tadd atmosphere tests, fix many grid-related issues - Granular.jl - Julia package for granular dynamics simulation
(HTM) git clone git://src.adamsgaard.dk/Granular.jl
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
(DIR) commit a86d2ec05f76ac24aabe853af5b7da9296365c23
(DIR) parent f8aa00babd9ab63c2f2454e05f18043c6f62a1d8
(HTM) Author: Anders Damsgaard <andersd@riseup.net>
Date: Sat, 3 Jun 2017 10:08:38 -0400
add atmosphere tests, fix many grid-related issues
Diffstat:
M src/atmosphere.jl | 12 ++++--------
M src/contact_search.jl | 26 ++++++++++++++++----------
M src/grid.jl | 22 +++++++++++-----------
M src/simulation.jl | 2 +-
A test/atmosphere.jl | 127 +++++++++++++++++++++++++++++++
M test/contact-search-and-geometry.jl | 6 +++---
M test/grid.jl | 41 +++++++++++++++++++++++++++++++
M test/runtests.jl | 1 +
8 files changed, 204 insertions(+), 33 deletions(-)
---
(DIR) diff --git a/src/atmosphere.jl b/src/atmosphere.jl
t@@ -114,19 +114,16 @@ function createRegularAtmosphereGrid(n::Array{Int, 1},
yh = repmat(linspace(origo[2] + .5*dx[2], L[2] - .5*dx[2], n[2])', n[1], 1)
zl = -linspace(.5*dx[3], L[3] - .5*dx[3], n[3])
- zi = -linspace(0., 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 Atmosphere(name,
time,
xq, yq,
xh, yh,
- zl, zi,
- u, v, h, e,
+ zl,
+ u, v,
Array{Array{Int, 1}}(size(xh, 1), size(xh, 2)))
end
t@@ -140,8 +137,7 @@ function addAtmosphereDrag!(simulation::Simulation)
error("no atmosphere data read")
end
- u, v, h, e = interpolateAtmosphereState(simulation.atmosphere,
- simulation.time)
+ u, v = interpolateAtmosphereState(simulation.atmosphere, simulation.time)
for ice_floe in simulation.ice_floes
t@@ -203,7 +199,7 @@ function applyAtmosphereVorticityToIceFloe!(ice_floe::IceFloeCylindrical,
rho_a = 1.2754 # atmosphere density
ice_floe.torque +=
- pi*ice_floe.areal_radius^4.*rho_o*
+ pi*ice_floe.areal_radius^4.*rho_a*
(ice_floe.areal_radius/5.*ice_floe.atmosphere_drag_coeff_horiz +
freeboard*ice_floe.atmosphere_drag_coeff_vert)*
abs(.5*atmosphere_curl - ice_floe.ang_vel)*
(DIR) diff --git a/src/contact_search.jl b/src/contact_search.jl
t@@ -24,10 +24,10 @@ function findContacts!(simulation::Simulation;
findContactsAllToAll!(simulation)
elseif method == "ocean grid"
- findContactsOceanGrid!(simulation)
+ findContactsInGrid!(simulation, simulation.ocean)
elseif method == "atmosphere grid"
- error("not yet implemented")
+ findContactsInGrid!(simulation, simulation.atmosphere)
else
error("Unknown contact search method '$method'")
t@@ -79,19 +79,25 @@ function findContactsAllToAll!(simulation::Simulation)
end
end
-export findContactsOceanGrid!
+export findContactsInGrid!
"""
- findContactsOceanGrid!(simulation)
+ findContactsInGrid!(simulation)
Perform an O(n*log(n)) cell-based contact search between all ice floes in the
-`simulation` object. Contacts between fixed ice floes are ignored.
+`simulation` object. Contacts between fixed or disabled ice floes are ignored.
"""
-function findContactsOceanGrid!(simulation::Simulation)
+function findContactsInGrid!(simulation::Simulation, grid::Any)
for idx_i = 1:length(simulation.ice_floes)
- grid_pos = simulation.ice_floes[idx_i].ocean_grid_pos
- nx, ny = size(simulation.ocean.xh)
+ if typeof(grid) == Ocean
+ grid_pos = simulation.ice_floes[idx_i].ocean_grid_pos
+ elseif typeof(grid) == Atmosphere
+ grid_pos = simulation.ice_floes[idx_i].atmosphere_grid_pos
+ else
+ error("grid type not understood")
+ end
+ nx, ny = size(grid.xh)
for i=(grid_pos[1] - 1):(grid_pos[1] + 1)
for j=(grid_pos[2] - 1):(grid_pos[2] + 1)
t@@ -101,7 +107,7 @@ function findContactsOceanGrid!(simulation::Simulation)
continue
end
- for idx_j in simulation.ocean.ice_floe_list[i, j]
+ for idx_j in grid.ice_floe_list[i, j]
checkAndAddContact!(simulation, idx_i, idx_j)
end
end
t@@ -109,7 +115,7 @@ function findContactsOceanGrid!(simulation::Simulation)
end
end
-export addContact!
+export checkAndAddContact!
"""
checkAndAddContact!(simulation, i, j)
(DIR) diff --git a/src/grid.jl b/src/grid.jl
t@@ -158,26 +158,26 @@ end
export findCellContainingPoint
"""
- findCellContainingPoint(ocean, point[, method])
+ findCellContainingPoint(grid, point[, method])
-Returns the `i`, `j` index of the ocean grid cell containing the `point`.
+Returns the `i`, `j` index of the grid cell containing the `point`.
The function uses either an area-based approach (`method = "Area"`), or a
conformal mapping approach (`method = "Conformal"`). The area-based approach is
more robust. This function returns the coordinates of the cell. If no match is
found the function returns `(0,0)`.
# Arguments
-* `ocean::Ocean`: object containing ocean data.
+* `grid::Any`: grid object containing ocean or atmosphere data.
* `point::Array{float, 1}`: two-dimensional vector of point to check.
* `method::String`: approach to use for determining if point is inside cell or
not, can be "Conformal" (default) or "Areal".
"""
-function findCellContainingPoint(ocean::Ocean, point::Array{float, 1};
+function findCellContainingPoint(grid::Any, point::Array{float, 1};
method::String="Conformal")
- for i=1:size(ocean.xh, 1)
- for j=1:size(ocean.yh, 2)
- if isPointInCell(ocean, i, j, point, method=method)
+ for i=1:size(grid.xh, 1)
+ for j=1:size(grid.yh, 2)
+ if isPointInCell(grid, i, j, point, method=method)
return i, j
end
end
t@@ -188,22 +188,22 @@ end
export getNonDimensionalCellCoordinates
"""
Returns the non-dimensional conformal mapped coordinates for point `point` in
-cell `i,j`, based off the coordinates in the `ocean` grid.
+cell `i,j`, based off the coordinates in the grid.
This function is a wrapper for `getCellCornerCoordinates()` and
`conformalQuadrilateralCoordinates()`.
"""
-function getNonDimensionalCellCoordinates(ocean::Ocean, i::Int, j::Int,
+function getNonDimensionalCellCoordinates(grid::Any, i::Int, j::Int,
point::Array{float, 1})
- sw, se, ne, nw = getCellCornerCoordinates(ocean, i, j)
+ sw, se, ne, nw = getCellCornerCoordinates(grid, i, j)
x_tilde, y_tilde = conformalQuadrilateralCoordinates(sw, se, ne, nw, point)
return [x_tilde, y_tilde]
end
export isPointInCell
"""
-Check if a 2d point is contained inside a cell from the ocean grid.
+Check if a 2d point is contained inside a cell from the supplied grid.
The function uses either an area-based approach (`method = "Area"`), or a
conformal mapping approach (`method = "Conformal"`). The area-based approach is
more robust. This function returns `true` or `false`.
(DIR) diff --git a/src/simulation.jl b/src/simulation.jl
t@@ -111,7 +111,7 @@ function run!(simulation::Simulation;
zeroForcesAndTorques!(simulation)
if typeof(simulation.atmosphere.input_file) != Bool
- sortIceFloesInGrid!(sim, sim.atmosphere)
+ sortIceFloesInGrid!(simulation, simulation.atmosphere)
end
if typeof(simulation.ocean.input_file) != Bool
(DIR) diff --git a/test/atmosphere.jl b/test/atmosphere.jl
t@@ -0,0 +1,127 @@
+#!/usr/bin/env julia
+
+# Check if atmosphere-specific functions and grid operations are functioning
+# correctly
+
+info("#### $(basename(@__FILE__)) ####")
+
+info("Testing regular grid generation")
+sim = SeaIce.createSimulation()
+sim.atmosphere = SeaIce.createRegularAtmosphereGrid([6, 6, 6], [1., 1., 1.])
+@test size(sim.atmosphere.xq) == (7, 7)
+@test size(sim.atmosphere.yq) == (7, 7)
+@test size(sim.atmosphere.xh) == (6, 6)
+@test size(sim.atmosphere.yh) == (6, 6)
+@test sim.atmosphere.xq[1, :, 1] ≈ zeros(7)
+@test sim.atmosphere.xq[4, :, 1] ≈ .5*ones(7)
+@test sim.atmosphere.xq[end, :, 1] ≈ 1.*ones(7)
+@test sim.atmosphere.yq[:, 1, 1] ≈ zeros(7)
+@test sim.atmosphere.yq[:, 4, 1] ≈ .5*ones(7)
+@test sim.atmosphere.yq[:, end, 1] ≈ 1.*ones(7)
+@test size(sim.atmosphere.u) == (7, 7, 6, 1)
+@test size(sim.atmosphere.v) == (7, 7, 6, 1)
+@test sim.atmosphere.u ≈ zeros(7, 7, 6, 1)
+@test sim.atmosphere.v ≈ zeros(7, 7, 6, 1)
+
+info("Testing velocity drag interaction (static atmosphere)")
+SeaIce.addIceFloeCylindrical(sim, [.5, .5], .25, .1)
+SeaIce.setTotalTime!(sim, 5.)
+SeaIce.setTimeStep!(sim)
+sim_init = deepcopy(sim)
+sim.ice_floes[1].lin_vel[1] = 0.1
+E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+SeaIce.run!(sim, verbose=false)
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test E_kin_rot_init ≈ E_kin_rot_final # no rotation before or after
+@test E_kin_lin_init > E_kin_lin_final # linear velocity lost due to atmos drag
+
+info("Testing velocity drag interaction (static ice floe)")
+sim = deepcopy(sim_init)
+sim.atmosphere.v[:, :, 1, 1] = 0.1
+E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+SeaIce.run!(sim, verbose=false)
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test E_kin_rot_init ≈ E_kin_rot_final # no rotation before or after
+@test E_kin_lin_init < E_kin_lin_final # linear vel. gained due to atmos drag
+
+info("Testing vortex interaction (static atmosphere)")
+sim = deepcopy(sim_init)
+sim.ice_floes[1].ang_vel = 0.1
+E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+SeaIce.run!(sim, verbose=false)
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test E_kin_rot_init > E_kin_rot_final # energy lost to atmosphere
+@test sim.ice_floes[1].ang_vel > 0. # check angular velocity orientation
+@test sim.ice_floes[1].ang_pos > 0. # check angular position orientation
+@test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained
+
+info("Testing vortex interaction (static ice floe)")
+sim = deepcopy(sim_init)
+sim.atmosphere = SeaIce.createRegularAtmosphereGrid([1, 1, 1], [1., 1., 1.])
+sim.ice_floes[1].lin_pos[1] = 0.5
+sim.ice_floes[1].lin_pos[2] = 0.5
+sim.atmosphere.v[1, :, 1, 1] = -0.1
+sim.atmosphere.v[2, :, 1, 1] = 0.1
+E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+SeaIce.run!(sim, verbose=false)
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test sim.ice_floes[1].ang_vel > 0. # check angular velocity orientation
+@test sim.ice_floes[1].ang_pos > 0. # check angular position orientation
+@test E_kin_rot_init < E_kin_rot_final # rotation after due to atm vortex
+@test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained
+
+sim = deepcopy(sim_init)
+sim.atmosphere = SeaIce.createRegularAtmosphereGrid([1, 1, 1], [1., 1., 1.])
+sim.ice_floes[1].lin_pos[1] = 0.5
+sim.ice_floes[1].lin_pos[2] = 0.5
+sim.atmosphere.v[1, :, 1, 1] = 0.1
+sim.atmosphere.v[2, :, 1, 1] = -0.1
+E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+SeaIce.run!(sim, verbose=false)
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test sim.ice_floes[1].ang_vel < 0. # check angular velocity orientation
+@test sim.ice_floes[1].ang_pos < 0. # check angular position orientation
+@test E_kin_rot_init < E_kin_rot_final # rotation after due to atm vortex
+@test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained
+
+sim = deepcopy(sim_init)
+sim.atmosphere = SeaIce.createRegularAtmosphereGrid([1, 1, 1], [1., 1., 1.])
+sim.ice_floes[1].lin_pos[1] = 0.5
+sim.ice_floes[1].lin_pos[2] = 0.5
+sim.atmosphere.u[:, 1, 1, 1] = -0.1
+sim.atmosphere.u[:, 2, 1, 1] = 0.1
+E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+SeaIce.run!(sim, verbose=false)
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test sim.ice_floes[1].ang_vel < 0. # check angular velocity orientation
+@test sim.ice_floes[1].ang_pos < 0. # check angular position orientation
+@test E_kin_rot_init < E_kin_rot_final # rotation after due to atm vortex
+@test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained
+
+sim = deepcopy(sim_init)
+sim.atmosphere = SeaIce.createRegularAtmosphereGrid([1, 1, 1], [1., 1., 1.])
+sim.ice_floes[1].lin_pos[1] = 0.5
+sim.ice_floes[1].lin_pos[2] = 0.5
+sim.atmosphere.u[:, 1, 1, 1] = 0.1
+sim.atmosphere.u[:, 2, 1, 1] = -0.1
+E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+SeaIce.run!(sim, verbose=false)
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test sim.ice_floes[1].ang_vel > 0. # check angular velocity orientation
+@test sim.ice_floes[1].ang_pos > 0. # check angular position orientation
+@test E_kin_rot_init < E_kin_rot_final # rotation after due to atm vortex
+@test E_kin_lin_init ≈ E_kin_lin_final # no linear velocity gained
(DIR) diff --git a/test/contact-search-and-geometry.jl b/test/contact-search-and-geometry.jl
t@@ -126,7 +126,7 @@ info("Testing findContactsGrid(...)")
sim = deepcopy(sim_copy)
sim.ocean = SeaIce.createRegularOceanGrid([4, 4, 2], [80., 80., 2.])
SeaIce.sortIceFloesInGrid!(sim, sim.ocean)
-SeaIce.findContactsOceanGrid!(sim)
+SeaIce.findContactsInGrid!(sim, sim.ocean)
@test 2 == sim.ice_floes[1].contacts[1]
for ic=2:32
t@@ -145,7 +145,7 @@ sim = deepcopy(sim_copy)
sim.ocean = SeaIce.createRegularOceanGrid([4, 4, 2], [80., 80., 2.])
sim.ice_floes[1].fixed = true
SeaIce.sortIceFloesInGrid!(sim, sim.ocean)
-SeaIce.findContactsOceanGrid!(sim)
+SeaIce.findContactsInGrid!(sim, sim.ocean)
@test 2 == sim.ice_floes[1].contacts[1]
for ic=2:32
t@@ -165,7 +165,7 @@ sim.ocean = SeaIce.createRegularOceanGrid([4, 4, 2], [80., 80., 2.])
sim.ice_floes[1].fixed = true
sim.ice_floes[2].fixed = true
SeaIce.sortIceFloesInGrid!(sim, sim.ocean)
-SeaIce.findContactsOceanGrid!(sim)
+SeaIce.findContactsInGrid!(sim, sim.ocean)
for ic=1:32
@test 0 == sim.ice_floes[1].contacts[ic]
(DIR) diff --git a/test/grid.jl b/test/grid.jl
t@@ -152,3 +152,44 @@ ocean.u[2, 2, 1, 1] = 1.0
ocean.u[1, 2, 1, 1] = 1.0
ocean.v[:, :, 1, 1] = 0.0
@test SeaIce.curl(ocean, .5, .5, 1, 1, 1, 1) < 0.
+
+info("Testing atmosphere drag")
+sim = SeaIce.createSimulation()
+sim.atmosphere = SeaIce.createRegularAtmosphereGrid([4, 4, 2], [4., 4., 2.])
+atmosphere = SeaIce.createRegularAtmosphereGrid([4, 4, 2], [4., 4., 2.])
+sim.atmosphere.u[:,:,1,1] = 5.
+SeaIce.addIceFloeCylindrical(sim, [2.5, 3.5], 1., 1., verbose=verbose)
+SeaIce.addIceFloeCylindrical(sim, [2.6, 2.5], 1., 1., verbose=verbose)
+SeaIce.sortIceFloesInGrid!(sim, sim.atmosphere, verbose=verbose)
+sim.time = ocean.time[1]
+SeaIce.addAtmosphereDrag!(sim)
+@test sim.ice_floes[1].force[1] > 0.
+@test sim.ice_floes[1].force[2] ≈ 0.
+@test sim.ice_floes[2].force[1] > 0.
+@test sim.ice_floes[2].force[2] ≈ 0.
+sim.atmosphere.u[:,:,1,1] = -5.
+sim.atmosphere.v[:,:,1,1] = 5.
+SeaIce.addIceFloeCylindrical(sim, [2.5, 3.5], 1., 1., verbose=verbose)
+SeaIce.addIceFloeCylindrical(sim, [2.6, 2.5], 1., 1., verbose=verbose)
+SeaIce.sortIceFloesInGrid!(sim, sim.atmosphere, verbose=verbose)
+sim.time = ocean.time[1]
+SeaIce.addAtmosphereDrag!(sim)
+@test sim.ice_floes[1].force[1] < 0.
+@test sim.ice_floes[1].force[2] > 0.
+@test sim.ice_floes[2].force[1] < 0.
+@test sim.ice_floes[2].force[2] > 0.
+
+info("Testing curl function")
+atmosphere.u[1, 1, 1, 1] = 1.0
+atmosphere.u[2, 1, 1, 1] = 1.0
+atmosphere.u[2, 2, 1, 1] = 0.0
+atmosphere.u[1, 2, 1, 1] = 0.0
+atmosphere.v[:, :, 1, 1] = 0.0
+@test SeaIce.curl(atmosphere, .5, .5, 1, 1, 1, 1) > 0.
+
+atmosphere.u[1, 1, 1, 1] = 0.0
+atmosphere.u[2, 1, 1, 1] = 0.0
+atmosphere.u[2, 2, 1, 1] = 1.0
+atmosphere.u[1, 2, 1, 1] = 1.0
+atmosphere.v[:, :, 1, 1] = 0.0
+@test SeaIce.curl(atmosphere, .5, .5, 1, 1, 1, 1) < 0.
(DIR) diff --git a/test/runtests.jl b/test/runtests.jl
t@@ -12,3 +12,4 @@ include("netcdf.jl")
include("vtk.jl")
include("grid.jl")
include("ocean.jl")
+include("atmosphere.jl")