Granular.jl

Julia package for granular dynamics simulation
git clone git://src.adamsgaard.dk/Granular.jl
Log | Files | Refs | README | LICENSE

commit cdd40dffc47470cc2d0513fb561b6275eed0f98b
parent 1d1b5457fd62a9751f2983af9432ade2cfe9b8f5
Author: Anders Damsgaard <andersd@riseup.net>
Date:   Thu, 27 Apr 2017 15:45:36 -0400

add functionality to create idealized oceans

Diffstat:
Msrc/ocean.jl | 59++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-
Atest/ocean.jl | 27+++++++++++++++++++++++++++
Mtest/runtests.jl | 1+
3 files changed, 86 insertions(+), 1 deletion(-)

diff --git a/src/ocean.jl b/src/ocean.jl @@ -1,3 +1,4 @@ +export createEmptyOcean "Returns empty ocean type for initialization purposes." function createEmptyOcean() return Ocean(false, @@ -18,6 +19,7 @@ function createEmptyOcean() zeros(1,1,1,1)) end +export readOceanNetCDF """ Read ocean NetCDF files generated by MOM6 from disk and return as `Ocean` data structure. @@ -59,6 +61,7 @@ function readOceanNetCDF(velocity_file::String, grid_file::String) return ocean end +export readOceanStateNetCDF """ Read NetCDF file with ocean state generated by MOM6 (e.g. `prog__####_###.nc` or `########.ocean_month.nc`) from disk and return time stamps, velocity fields, @@ -97,6 +100,7 @@ function readOceanStateNetCDF(filename::String) return time, u, v, h, e, zl, zi end +export readOceanGridNetCDF """ Read NetCDF file with ocean *supergrid* information generated by MOM6 (e.g. `ocean_hrid.nc`) from disk and return as `Ocean` data structure. This file is @@ -125,6 +129,7 @@ function readOceanGridNetCDF(filename::String) return xh, yh, xq, yq end +export interpolateOceanVelocitiesToCorners """ Convert gridded data from Arakawa-C type (decomposed velocities at faces) to Arakawa-B type (velocities at corners) through interpolation. @@ -158,6 +163,7 @@ function interpolateOceanVelocitiesToCorners(u_in::Array{float, 4}, return u, v end +export interpolateOceanState """ Ocean data is containted in `Ocean` type at discrete times (`Ocean.time`). This function performs linear interpolation between time steps to get the approximate @@ -194,8 +200,46 @@ function interpolateOceanState(ocean::Ocean, t::float) ocean.e[:,:,:,i]*(1. - rel_time) + ocean.e[:,:,:,i+1]*rel_time end +export createRegularOceanGrid """ -Add Stokes-type drag from velocity difference between ocean and ice floe. +Initialize and return a regular, Cartesian `Ocean` grid with `n[1]` by `n[2]` +cells in the horizontal dimension, and `n[3]` vertical cells. The cell corner +and center coordinates will be set according to the grid spatial dimensions +`L[1]`, `L[2]`, and `L[3]`. The grid `u`, `v`, `h`, and `e` fields will contain +one 4-th dimension matrix per `time` step. +""" +function createRegularOceanGrid(n::Array{Int, 1}, + L::Array{float, 1}; + origo::Array{float, 1} = zeros(3), + time::Array{float, 1} = zeros(1), + name::String = "unnamed") + + xq = repmat(linspace(origo[1], L[1], n[1] + 1), 1, n[2] + 1) + yq = repmat(linspace(origo[2], L[2], n[2] + 1)', n[1] + 1, 1) + + dx = L./n + xh = repmat(linspace(origo[1] + .5*dx[1], L[1] - .5*dx[1], n[1]), 1, n[2]) + yh = repmat(linspace(origo[2] + .5*dx[2], L[2] - .5*dx[2], n[2])', n[1], 1) + + zl = linspace(origo[3] + .5*dx[3], L[3] - .5*dx[3], n[3]) + zi = linspace(origo[3], L[3], n[3] + 1) + + u = zeros(n[1] + 1, n[2] + 1, n[3], length(time)) + v = zeros(n[1] + 1, n[2] + 1, n[3], length(time)) + h = zeros(n[1] + 1, n[2] + 1, n[3], length(time)) + e = zeros(n[1] + 1, n[2] + 1, n[3], length(time)) + + return Ocean("unnamed", + time, + xq, yq, + xh, yh, + zl, zi, + u, v, h, e) +end + +export addOceanDrag +""" +Add Stokes-type drag from velocity difference between ocean and all ice floes. """ function addOceanDrag!(simulation::Simulation) if !simulation.ocean.id @@ -204,4 +248,17 @@ function addOceanDrag!(simulation::Simulation) u, v, h, e = interpolateOceanState(simulation.ocean, simulation.time) + for ice_floe in simulation.ice_floes + applyOceanDragToIceFloe!(ice_floe, u, v) + end +end + +export applyOceanDragToIceFloe +""" +Add Stokes-type drag from velocity difference between ocean and a single ice +floe. +""" +function applyOceanDragToIceFloe!(ice_floe::IceFloeCylindrical, + u::float, v::float) + end diff --git a/test/ocean.jl b/test/ocean.jl @@ -0,0 +1,27 @@ +#!/usr/bin/env julia + +# Check if ocean-specific functions and grid operations are functioning +# correctly + +info("#### $(basename(@__FILE__)) ####") + +info("Testing regular grid generation") +ocean = SeaIce.createRegularOceanGrid([6, 6, 6], [1., 1., 1.]) +@test size(ocean.xq) == (7, 7) +@test size(ocean.yq) == (7, 7) +@test size(ocean.xh) == (6, 6) +@test size(ocean.yh) == (6, 6) +@test ocean.xq[1, :, 1] ≈ zeros(7) +@test ocean.xq[4, :, 1] ≈ .5*ones(7) +@test ocean.xq[end, :, 1] ≈ 1.*ones(7) +@test ocean.yq[:, 1, 1] ≈ zeros(7) +@test ocean.yq[:, 4, 1] ≈ .5*ones(7) +@test ocean.yq[:, end, 1] ≈ 1.*ones(7) +@test size(ocean.u) == (7, 7, 6, 1) +@test size(ocean.v) == (7, 7, 6, 1) +@test size(ocean.h) == (7, 7, 6, 1) +@test size(ocean.e) == (7, 7, 6, 1) +@test ocean.u ≈ zeros(7, 7, 6, 1) +@test ocean.v ≈ zeros(7, 7, 6, 1) +@test ocean.h ≈ zeros(7, 7, 6, 1) +@test ocean.e ≈ zeros(7, 7, 6, 1) diff --git a/test/runtests.jl b/test/runtests.jl @@ -6,3 +6,4 @@ include("collision-2floes-normal.jl") include("netcdf.jl") include("vtk.jl") include("grid.jl") +include("ocean.jl")