commit b189f7b4e3f1e021691cc3a766b66bf2f2557538
parent 87def518c59023ee58581b32919e0aa9cea9ed6c
Author: Anders Damsgaard <andersd@riseup.net>
Date: Fri, 8 Sep 2017 12:03:06 -0400
Merge branch 'master' of github.com:anders-dc/SeaIce.jl
Diffstat:
22 files changed, 1031 insertions(+), 413 deletions(-)
diff --git a/REQUIRE b/REQUIRE
@@ -1,3 +1,4 @@
julia 0.6
WriteVTK
NetCDF
+PyPlot
diff --git a/src/SeaIce.jl b/src/SeaIce.jl
@@ -19,5 +19,6 @@ include("temporal_integration.jl")
include("io.jl")
include("ocean.jl")
include("atmosphere.jl")
+include("util.jl")
end # module end
diff --git a/src/atmosphere.jl b/src/atmosphere.jl
@@ -15,7 +15,7 @@ function createEmptyAtmosphere()
zeros(1,1,1,1),
zeros(1,1,1,1),
- Array{Array{Int, 1}}(1, 1),
+ Array{Vector{Int}}(1, 1),
false)
end
@@ -25,16 +25,16 @@ export interpolateAtmosphereVelocitiesToCorners
Convert gridded data from Arakawa-C type (decomposed velocities at faces) to
Arakawa-B type (velocities at corners) through interpolation.
"""
-function interpolateAtmosphereVelocitiesToCorners(u_in::Array{float, 4},
- v_in::Array{float, 4})
+function interpolateAtmosphereVelocitiesToCorners(u_in::Array{Float64, 4},
+ v_in::Array{Float64, 4})
if size(u_in) != size(v_in)
error("size of u_in ($(size(u_in))) must match v_in ($(size(v_in)))")
end
nx, ny, nz, nt = size(u_in)
- #u = Array{float}(nx+1, ny+1, nz, nt)
- #v = Array{float}(nx+1, ny+1, nz, nt)
+ #u = Array{Float64}(nx+1, ny+1, nz, nt)
+ #v = Array{Float64}(nx+1, ny+1, nz, nt)
u = zeros(nx+1, ny+1, nz, nt)
v = zeros(nx+1, ny+1, nz, nt)
for i=1:nx
@@ -62,7 +62,7 @@ steps to get the approximate atmosphere state at any point in time. If the
`Atmosphere` data set only contains a single time step, values from that time
are returned.
"""
-function interpolateAtmosphereState(atmosphere::Atmosphere, t::float)
+function interpolateAtmosphereState(atmosphere::Atmosphere, t::Float64)
if length(atmosphere.time) == 1
return atmosphere.u, atmosphere.v
elseif t < atmosphere.time[1] || t > atmosphere.time[end]
@@ -102,10 +102,10 @@ one 4-th dimension matrix per `time` step. Sea surface will be at `z=0.` with
the atmosphere spanning `z<0.`. Vertical indexing starts with `k=0` at the sea
surface, and increases downwards.
"""
-function createRegularAtmosphereGrid(n::Array{Int, 1},
- L::Array{float, 1};
- origo::Array{float, 1} = zeros(2),
- time::Array{float, 1} = zeros(1),
+function createRegularAtmosphereGrid(n::Vector{Int},
+ L::Vector{Float64};
+ origo::Vector{Float64} = zeros(2),
+ time::Array{Float64, 1} = zeros(1),
name::String = "unnamed")
xq = repmat(linspace(origo[1], L[1], n[1] + 1), 1, n[2] + 1)
@@ -141,6 +141,11 @@ function addAtmosphereDrag!(simulation::Simulation)
end
u, v = interpolateAtmosphereState(simulation.atmosphere, simulation.time)
+ uv_interp = Vector{Float64}(2)
+ sw = Vector{Float64}(2)
+ se = Vector{Float64}(2)
+ ne = Vector{Float64}(2)
+ nw = Vector{Float64}(2)
for ice_floe in simulation.ice_floes
@@ -163,17 +168,13 @@ function addAtmosphereDrag!(simulation::Simulation)
""")
end
- applyAtmosphereDragToIceFloe!(ice_floe,
- bilinearInterpolation(u,
- x_tilde, y_tilde,
- i, j, k, 1),
- bilinearInterpolation(v,
- x_tilde, y_tilde,
- i, j, k, 1))
+ bilinearInterpolation!(uv_interp, u, v, x_tilde, y_tilde, i, j, k, 1)
+ applyAtmosphereDragToIceFloe!(ice_floe, uv_interp[1], uv_interp[2])
applyAtmosphereVorticityToIceFloe!(ice_floe,
- curl(simulation.atmosphere,
- x_tilde, y_tilde, i, j, k, 1))
+ curl(simulation.atmosphere, x_tilde, y_tilde,
+ i, j, k, 1, sw, se, ne, nw))
end
+ nothing
end
export applyAtmosphereDragToIceFloe!
@@ -182,18 +183,19 @@ Add Stokes-type drag from velocity difference between atmosphere and a single
ice floe.
"""
function applyAtmosphereDragToIceFloe!(ice_floe::IceFloeCylindrical,
- u::float, v::float)
- freeboard = .1*ice_floe.thickness # height above water
+ u::Float64, v::Float64)
rho_a = 1.2754 # atmosphere density
length = ice_floe.areal_radius*2.
width = ice_floe.areal_radius*2.
- drag_force = rho_a * (.5*ice_floe.ocean_drag_coeff_vert*width*freeboard +
- ice_floe.atmosphere_drag_coeff_horiz*length*width) *
+ drag_force = rho_a *
+ (.5*ice_floe.ocean_drag_coeff_vert*width*.1*ice_floe.thickness +
+ ice_floe.atmosphere_drag_coeff_horiz*length*width) *
([u, v] - ice_floe.lin_vel)*norm([u, v] - ice_floe.lin_vel)
ice_floe.force += drag_force
ice_floe.atmosphere_stress = drag_force/ice_floe.horizontal_surface_area
+ nothing
end
export applyAtmosphereVorticityToIceFloe!
@@ -203,16 +205,16 @@ single ice floe. See Eq. 9.28 in "Introduction to Fluid Mechanics" by Nakayama
and Boucher, 1999.
"""
function applyAtmosphereVorticityToIceFloe!(ice_floe::IceFloeCylindrical,
- atmosphere_curl::float)
- freeboard = .1*ice_floe.thickness # height above water
+ atmosphere_curl::Float64)
rho_a = 1.2754 # atmosphere density
ice_floe.torque +=
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)*
+ .1*ice_floe.thickness*ice_floe.atmosphere_drag_coeff_vert)*
abs(.5*atmosphere_curl - ice_floe.ang_vel)*
(.5*atmosphere_curl - ice_floe.ang_vel)
+ nothing
end
export compareAtmospheres
@@ -240,4 +242,5 @@ function compareAtmospheres(atmosphere1::Atmosphere, atmosphere2::Atmosphere)
if isassigned(atmosphere1.ice_floe_list, 1)
Base.Test.@test atmosphere1.ice_floe_list == atmosphere2.ice_floe_list
end
+ nothing
end
diff --git a/src/contact_search.jl b/src/contact_search.jl
@@ -32,6 +32,7 @@ function findContacts!(simulation::Simulation;
else
error("Unknown contact search method '$method'")
end
+ nothing
end
export interIceFloePositionVector
@@ -47,18 +48,19 @@ Returns a `vector` pointing from ice floe `i` to ice floe `j` in the
* `j::Int`: index of the second ice floe.
"""
function interIceFloePositionVector(simulation::Simulation,
- i::Integer, j::Integer)
- return simulation.ice_floes[i].lin_pos - simulation.ice_floes[j].lin_pos
+ i::Int, j::Int)
+ @inbounds return simulation.ice_floes[i].lin_pos -
+ simulation.ice_floes[j].lin_pos
end
"""
position_ij is the inter-grain position vector, and can be found with
interIceFloePositionVector().
"""
-function findOverlap(simulation::Simulation, i::Integer, j::Integer,
- position_ij::vector)
- return norm(position_ij) - (simulation.ice_floes[i].contact_radius +
- simulation.ice_floes[j].contact_radius)
+function findOverlap(simulation::Simulation, i::Int, j::Int,
+ position_ij::Vector{Float64})
+ @inbounds return norm(position_ij) - (simulation.ice_floes[i].contact_radius
+ + simulation.ice_floes[j].contact_radius)
end
export findContactsAllToAll!
@@ -70,13 +72,14 @@ Perform an O(n^2) all-to-all contact search between all ice floes in the
"""
function findContactsAllToAll!(simulation::Simulation)
- for i = 1:length(simulation.ice_floes)
+ @inbounds for i = 1:length(simulation.ice_floes)
# Check contacts with other grains
for j = 1:length(simulation.ice_floes)
checkAndAddContact!(simulation, i, j)
end
end
+ nothing
end
export findContactsInGrid!
@@ -107,12 +110,13 @@ function findContactsInGrid!(simulation::Simulation, grid::Any)
continue
end
- for idx_j in grid.ice_floe_list[i, j]
+ @inbounds for idx_j in grid.ice_floe_list[i, j]
checkAndAddContact!(simulation, idx_i, idx_j)
end
end
end
end
+ nothing
end
export checkAndAddContact!
@@ -134,7 +138,7 @@ written to `simulation.contact_parallel_displacement`.
function checkAndAddContact!(sim::Simulation, i::Int, j::Int)
if i < j
- if (sim.ice_floes[i].fixed && sim.ice_floes[j].fixed) ||
+ @inbounds if (sim.ice_floes[i].fixed && sim.ice_floes[j].fixed) ||
!sim.ice_floes[i].enabled || !sim.ice_floes[j].enabled
return
end
@@ -143,28 +147,49 @@ function checkAndAddContact!(sim::Simulation, i::Int, j::Int)
position_ij = interIceFloePositionVector(sim, i, j)
overlap_ij = findOverlap(sim, i, j, position_ij)
+ contact_found = false
+
# Check if grains overlap (overlap when negative)
if overlap_ij < 0.
- for ic=1:(sim.Nc_max + 1)
- if ic == (sim.Nc_max + 1)
- error("contact $i-$j exceeds max. number of contacts " *
- "(sim.Nc_max = $(sim.Nc_max)) for ice floe $i")
-
- else
- if sim.ice_floes[i].contacts[ic] == j
- break # contact already registered
-
- elseif sim.ice_floes[i].contacts[ic] == 0 # empty
- sim.ice_floes[i].n_contacts += 1 # register new contact
- sim.ice_floes[j].n_contacts += 1
- sim.ice_floes[i].contacts[ic] = j
- fill!(sim.ice_floes[i].
+
+ # Check if contact is already registered
+ for ic=1:sim.Nc_max
+ @inbounds if sim.ice_floes[i].contacts[ic] == j
+ contact_found = true
+ break # contact already registered
+ end
+ end
+
+ # Register as new contact in first empty position
+ if !contact_found
+
+ for ic=1:(sim.Nc_max + 1)
+
+ # Test if this contact exceeds the number of contacts
+ if ic == (sim.Nc_max + 1)
+ for ic=1:sim.Nc_max
+ warn("ice_floes[$i].contacts[$ic] = " *
+ "$(sim.ice_floes[i].contacts[ic])")
+ warn("ice_floes[$i].contact_age[$ic] = " *
+ "$(sim.ice_floes[i].contact_age[ic])")
+ end
+ error("contact $i-$j exceeds max. number of contacts " *
+ "(sim.Nc_max = $(sim.Nc_max)) for ice floe $i")
+ end
+
+ # Register as new contact
+ @inbounds if sim.ice_floes[i].contacts[ic] == 0 # empty
+ @inbounds sim.ice_floes[i].n_contacts += 1
+ @inbounds sim.ice_floes[j].n_contacts += 1
+ @inbounds sim.ice_floes[i].contacts[ic] = j
+ @inbounds fill!(sim.ice_floes[i].
contact_parallel_displacement[ic] , 0.)
- sim.ice_floes[i].contact_age[ic] = 0.
+ @inbounds sim.ice_floes[i].contact_age[ic] = 0.
break
end
end
end
end
end
+ nothing
end
diff --git a/src/datatypes.jl b/src/datatypes.jl
@@ -1,35 +1,31 @@
-# Define floating point data types
-const float = Float64
-const vector = Array{Float64, 1}
-
## Particle composite types
type IceFloeCylindrical
# Material properties
- density::float
+ density::Float64
# Geometrical parameters
- thickness::float
- contact_radius::float
- areal_radius::float
- circumreference::float
- horizontal_surface_area::float
- side_surface_area::float
- volume::float
- mass::float
- moment_of_inertia::float
+ thickness::Float64
+ contact_radius::Float64
+ areal_radius::Float64
+ circumreference::Float64
+ horizontal_surface_area::Float64
+ side_surface_area::Float64
+ volume::Float64
+ mass::Float64
+ moment_of_inertia::Float64
# Linear kinematic degrees of freedom along horizontal plane
- lin_pos::vector
- lin_vel::vector
- lin_acc::vector
- force::vector
+ lin_pos::Vector{Float64}
+ lin_vel::Vector{Float64}
+ lin_acc::Vector{Float64}
+ force::Vector{Float64}
# Angular kinematic degrees of freedom for vertical rotation around center
- ang_pos::float
- ang_vel::float
- ang_acc::float
- torque::float
+ ang_pos::Float64
+ ang_vel::Float64
+ ang_acc::Float64
+ torque::Float64
# Kinematic constraint flags
fixed::Bool
@@ -37,97 +33,98 @@ type IceFloeCylindrical
enabled::Bool
# Rheological parameters
- contact_stiffness_normal::float
- contact_stiffness_tangential::float
- contact_viscosity_normal::float
- contact_viscosity_tangential::float
- contact_static_friction::float
- contact_dynamic_friction::float
-
- youngs_modulus::float
- poissons_ratio::float
- tensile_strength::float
- tensile_heal_rate::float
- compressive_strength_prefactor::float
+ contact_stiffness_normal::Float64
+ contact_stiffness_tangential::Float64
+ contact_viscosity_normal::Float64
+ contact_viscosity_tangential::Float64
+ contact_static_friction::Float64
+ contact_dynamic_friction::Float64
+
+ youngs_modulus::Float64
+ poissons_ratio::Float64
+ tensile_strength::Float64
+ tensile_heal_rate::Float64
+ compressive_strength_prefactor::Float64
# Ocean/atmosphere interaction parameters
- ocean_drag_coeff_vert::float
- ocean_drag_coeff_horiz::float
- atmosphere_drag_coeff_vert::float
- atmosphere_drag_coeff_horiz::float
+ ocean_drag_coeff_vert::Float64
+ ocean_drag_coeff_horiz::Float64
+ atmosphere_drag_coeff_vert::Float64
+ atmosphere_drag_coeff_horiz::Float64
# Interaction
- pressure::float
+ pressure::Float64
n_contacts::Int
- ocean_grid_pos::Array{Int, 1}
- atmosphere_grid_pos::Array{Int, 1}
- contacts::Array{Int, 1}
- contact_parallel_displacement::Array{Array{Float64, 1}, 1}
- contact_age::Array{Float64, 1}
-
- granular_stress::vector
- ocean_stress::vector
- atmosphere_stress::vector
+ ocean_grid_pos::Vector{Int}
+ atmosphere_grid_pos::Vector{Int}
+ contacts::Vector{Int}
+ contact_parallel_displacement::Vector{Vector{Float64}}
+ contact_age::Vector{Float64}
+
+ granular_stress::Vector{Float64}
+ ocean_stress::Vector{Float64}
+ atmosphere_stress::Vector{Float64}
end
# Type for gathering data from ice floe objects into single arrays
type IceFloeArrays
# Material properties
- density
+ density::Vector{Float64}
# Geometrical parameters
- thickness
- contact_radius
- areal_radius
- circumreference
- horizontal_surface_area
- side_surface_area
- volume
- mass
- moment_of_inertia
+ thickness::Vector{Float64}
+ contact_radius::Vector{Float64}
+ areal_radius::Vector{Float64}
+ circumreference::Vector{Float64}
+ horizontal_surface_area::Vector{Float64}
+ side_surface_area::Vector{Float64}
+ volume::Vector{Float64}
+ mass::Vector{Float64}
+ moment_of_inertia::Vector{Float64}
# Linear kinematic degrees of freedom along horizontal plane
- lin_pos
- lin_vel
- lin_acc
- force
+ lin_pos::Array{Float64, 2}
+ lin_vel::Array{Float64, 2}
+ lin_acc::Array{Float64, 2}
+ force::Array{Float64, 2}
# Angular kinematic degrees of freedom for vertical rotation around center
- ang_pos
- ang_vel
- ang_acc
- torque
+ ang_pos::Array{Float64, 2}
+ ang_vel::Array{Float64, 2}
+ ang_acc::Array{Float64, 2}
+ torque::Array{Float64, 2}
# Kinematic constraint flags
- fixed
- rotating
- enabled
+ fixed::Vector{Int}
+ rotating::Vector{Int}
+ enabled::Vector{Int}
# Rheological parameters
- contact_stiffness_normal
- contact_stiffness_tangential
- contact_viscosity_normal
- contact_viscosity_tangential
- contact_static_friction
- contact_dynamic_friction
-
- youngs_modulus
- poissons_ratio
- tensile_strength
- compressive_strength_prefactor
-
- ocean_drag_coeff_vert
- ocean_drag_coeff_horiz
- atmosphere_drag_coeff_vert
- atmosphere_drag_coeff_horiz
-
- pressure
- n_contacts
-
- granular_stress
- ocean_stress
- atmosphere_stress
+ contact_stiffness_normal::Vector{Float64}
+ contact_stiffness_tangential::Vector{Float64}
+ contact_viscosity_normal::Vector{Float64}
+ contact_viscosity_tangential::Vector{Float64}
+ contact_static_friction::Vector{Float64}
+ contact_dynamic_friction::Vector{Float64}
+
+ youngs_modulus::Vector{Float64}
+ poissons_ratio::Vector{Float64}
+ tensile_strength::Vector{Float64}
+ #tensile_heal_rate::Vector{Float64}
+ compressive_strength_prefactor::Vector{Float64}
+
+ ocean_drag_coeff_vert::Vector{Float64}
+ ocean_drag_coeff_horiz::Vector{Float64}
+ atmosphere_drag_coeff_vert::Vector{Float64}
+ atmosphere_drag_coeff_horiz::Vector{Float64}
+
+ pressure::Vector{Float64}
+ n_contacts::Vector{Int}
+
+ granular_stress::Array{Float64, 2}
+ ocean_stress::Array{Float64, 2}
+ atmosphere_stress::Array{Float64, 2}
end
#=
@@ -171,7 +168,7 @@ h-points. During read, the velocities are interpolated to the cell corners
type Ocean
input_file::Any
- time::Array{Float64, 1}
+ time::Vector{Float64}
# q-point (cell corner) positions
xq::Array{Float64, 2}
@@ -182,8 +179,8 @@ type Ocean
yh::Array{Float64, 2}
# Vertical positions
- zl::Array{Float64, 1}
- zi::Array{Float64, 1}
+ zl::Vector{Float64}
+ zi::Vector{Float64}
# Field values
u::Array{Float64, 4}
@@ -191,7 +188,7 @@ type Ocean
h::Array{Float64, 4}
e::Array{Float64, 4}
- ice_floe_list::Array{Array{Int, 1}, 2}
+ ice_floe_list::Array{Vector{Int}, 2}
end
#=
@@ -213,12 +210,12 @@ cell corners (q-points).
# Fields
* `input_file::String`: path to input NetCDF file
-* `time::Array{Float64, 1}`: time in days
-* `xq::Array{Float64, 1}`: nominal longitude of q-points [degrees_E]
-* `yq::Array{Float64, 1}`: nominal latitude of q-points [degrees_N]
-* `xh::Array{Float64, 1}`: nominal longitude of h-points [degrees_E]
-* `yh::Array{Float64, 1}`: nominal latitude of h-points [degrees_N]
-* `zl::Array{Float64, 1}`: vertical position [m]
+* `time::Vector{Float64}`: time in days
+* `xq::Array{Float64, 2}`: nominal longitude of q-points [degrees_E]
+* `yq::Array{Float64, 2}`: nominal latitude of q-points [degrees_N]
+* `xh::Array{Float64, 2}`: nominal longitude of h-points [degrees_E]
+* `yh::Array{Float64, 2}`: nominal latitude of h-points [degrees_N]
+* `zl::Vector{Float64}`: vertical position [m]
* `u::Array{Float64, Int}`: zonal velocity (positive towards west) [m/s],
dimensions correspond to placement in `[xq, yq, zl, time]`.
* `v::Array{Float64, Int}`: meridional velocity (positive towards north) [m/s],
@@ -229,7 +226,7 @@ cell corners (q-points).
type Atmosphere
input_file::Any
- time::Array{Float64, 1}
+ time::Vector{Float64}
# q-point (cell corner) positions
xq::Array{Float64, 2}
@@ -240,13 +237,13 @@ type Atmosphere
yh::Array{Float64, 2}
# Vertical positions
- zl::Array{Float64, 1}
+ zl::Vector{Float64}
# Field values
u::Array{Float64, 4}
v::Array{Float64, 4}
- ice_floe_list::Array{Array{Int, 1}, 2}
+ ice_floe_list::Array{Vector{Int}, 2}
# If true the grid positions are identical to the ocean grid
collocated_with_ocean_grid::Bool
@@ -264,7 +261,7 @@ type Simulation
file_number::Int
file_time_since_output_file::Float64
- ice_floes::Array{IceFloeCylindrical, 1}
+ ice_floes::Vector{IceFloeCylindrical}
ocean::Ocean
atmosphere::Atmosphere
diff --git a/src/grid.jl b/src/grid.jl
@@ -7,41 +7,65 @@ south-west (-x, -y)-facing corner.
# Arguments
* `field::Array{Float64, 4}`: a scalar field to interpolate from
-* `x_tilde::float`: x point position [0;1]
-* `y_tilde::float`: y point position [0;1]
+* `x_tilde::Float64`: x point position [0;1]
+* `y_tilde::Float64`: y point position [0;1]
* `i::Int`: i-index of cell containing point
* `j::Int`: j-index of scalar field to interpolate from
* `it::Int`: time step from scalar field to interpolate from
"""
-function bilinearInterpolation(field::Array{Float64, 4},
- x_tilde::Float64,
- y_tilde::Float64,
- i::Int,
- j::Int,
- k::Int,
- it::Int)
+@inline function bilinearInterpolation!(interp_val::Vector{Float64},
+ field_x::Array{Float64, 2},
+ field_y::Array{Float64, 2},
+ x_tilde::Float64,
+ y_tilde::Float64,
+ i::Int,
+ j::Int)
#if x_tilde < 0. || x_tilde > 1. || y_tilde < 0. || y_tilde > 1.
#error("relative coordinates outside bounds ($(x_tilde), $(y_tilde))")
#end
- return (field[i+1, j+1, k, it]*x_tilde +
- field[i, j+1, k, it]*(1. - x_tilde))*y_tilde +
- (field[i+1, j, k, it]*x_tilde +
- field[i, j, k, it]*(1. - x_tilde))*(1. - y_tilde)
+ x_tilde_inv = 1. - x_tilde
+
+ @views interp_val[1] =
+ (field_x[i+1, j+1]*x_tilde + field_x[i, j+1]*x_tilde_inv)*y_tilde +
+ (field_x[i+1, j]*x_tilde + field_x[i, j]*x_tilde_inv)*(1. - y_tilde)
+
+ @views interp_val[2] =
+ (field_y[i+1, j+1]*x_tilde + field_y[i, j+1]*x_tilde_inv)*y_tilde +
+ (field_y[i+1, j]*x_tilde + field_y[i, j]*x_tilde_inv)*(1. - y_tilde)
+
+ nothing
end
-function bilinearInterpolation(field::Array{Float64, 2},
- x_tilde::Float64,
- y_tilde::Float64,
- i::Int,
- j::Int)
-
- if x_tilde < 0. || x_tilde > 1. || y_tilde < 0. || y_tilde > 1.
- error("relative coordinates outside bounds ($(x_tilde), $(y_tilde))")
- end
+@inline function bilinearInterpolation!(interp_val::Vector{Float64},
+ field_x::Array{Float64, 4},
+ field_y::Array{Float64, 4},
+ x_tilde::Float64,
+ y_tilde::Float64,
+ i::Int,
+ j::Int,
+ k::Int,
+ it::Int)
+
+ #if x_tilde < 0. || x_tilde > 1. || y_tilde < 0. || y_tilde > 1.
+ #error("relative coordinates outside bounds ($(x_tilde), $(y_tilde))")
+ #end
- return (field[i+1, j+1]*x_tilde + field[i, j+1]*(1. - x_tilde))*y_tilde +
- (field[i+1, j]*x_tilde + field[i, j]*(1. - x_tilde))*(1. - y_tilde)
+ x_tilde_inv = 1. - x_tilde
+
+ @views interp_val[1] =
+ (field_x[i+1, j+1, k, it]*x_tilde +
+ field_x[i, j+1, k, it]*x_tilde_inv)*y_tilde +
+ (field_x[i+1, j, k, it]*x_tilde +
+ field_x[i, j, k, it]*x_tilde_inv)*(1. - y_tilde)
+
+ @views interp_val[2] =
+ (field_y[i+1, j+1, k, it]*x_tilde +
+ field_y[i, j+1, k, it]*x_tilde_inv)*y_tilde +
+ (field_y[i+1, j, k, it]*x_tilde +
+ field_y[i, j, k, it]*x_tilde_inv)*(1. - y_tilde)
+
+ nothing
end
"""
@@ -53,8 +77,8 @@ grid to an arbitrary position in a cell. Assumes south-west convention, i.e.
# Arguments
* `grid::Any`: grid for which to determine curl
-* `x_tilde::float`: x point position [0;1]
-* `y_tilde::float`: y point position [0;1]
+* `x_tilde::Float64`: x point position [0;1]
+* `y_tilde::Float64`: y point position [0;1]
* `i::Int`: i-index of cell containing point
* `j::Int`: j-index of scalar field to interpolate from
* `it::Int`: time step from scalar field to interpolate from
@@ -65,15 +89,27 @@ function curl(grid::Any,
i::Int,
j::Int,
k::Int,
- it::Int)
-
- sw, se, ne, nw = getCellCornerCoordinates(grid.xq, grid.yq, i, j)
+ it::Int,
+ sw::Vector{Float64} = Vector{Float64}(2),
+ se::Vector{Float64} = Vector{Float64}(2),
+ ne::Vector{Float64} = Vector{Float64}(2),
+ nw::Vector{Float64} = Vector{Float64}(2))
+
+ #sw, se, ne, nw = getCellCornerCoordinates(grid.xq, grid.yq, i, j)
+ sw[1] = grid.xq[ i, j]
+ sw[2] = grid.yq[ i, j]
+ se[1] = grid.xq[i+1, j]
+ se[2] = grid.yq[i+1, j]
+ ne[1] = grid.xq[i+1, j+1]
+ ne[2] = grid.yq[i+1, j+1]
+ nw[1] = grid.xq[ i, j+1]
+ nw[2] = grid.yq[ i, j+1]
sw_se = norm(sw - se)
se_ne = norm(se - ne)
nw_ne = norm(nw - ne)
sw_nw = norm(sw - nw)
- return (
+ @views @inbounds return (
((grid.v[i+1, j , k,it] - grid.v[i , j , k,it])/sw_se*(1. - y_tilde) +
((grid.v[i+1, j+1, k,it] - grid.v[i , j+1, k,it])/nw_ne)*y_tilde) -
((grid.u[i , j+1, k,it] - grid.u[i , j , k,it])/sw_nw*(1. - x_tilde) +
@@ -92,36 +128,42 @@ function sortIceFloesInGrid!(simulation::Simulation, grid::Any; verbose=false)
for i=1:size(grid.xh, 1)
for j=1:size(grid.xh, 2)
- grid.ice_floe_list[i, j] = Int[]
+ @inbounds grid.ice_floe_list[i, j] = Int[]
end
end
else
for i=1:size(grid.xh, 1)
for j=1:size(grid.xh, 2)
- empty!(grid.ice_floe_list[i, j])
+ @inbounds empty!(grid.ice_floe_list[i, j])
end
end
end
- for idx in 1:length(simulation.ice_floes)
+ sw = Vector{Float64}(2)
+ se = Vector{Float64}(2)
+ ne = Vector{Float64}(2)
+ nw = Vector{Float64}(2)
+
+ for idx=1:length(simulation.ice_floes)
- if !simulation.ice_floes[idx].enabled
+ @inbounds if !simulation.ice_floes[idx].enabled
continue
end
# After first iteration, check if ice floe is in same cell before
# traversing entire grid
if typeof(grid) == Ocean
- i_old, j_old = simulation.ice_floes[idx].ocean_grid_pos
+ @inbounds i_old, j_old = simulation.ice_floes[idx].ocean_grid_pos
elseif typeof(grid) == Atmosphere
- i_old, j_old = simulation.ice_floes[idx].atmosphere_grid_pos
+ @inbounds i_old, j_old =
+ simulation.ice_floes[idx].atmosphere_grid_pos
else
error("grid type not understood.")
end
if simulation.time > 0. &&
i_old > 0 && j_old > 0 &&
isPointInCell(grid, i_old, j_old,
- simulation.ice_floes[idx].lin_pos)
+ simulation.ice_floes[idx].lin_pos, sw, se, ne, nw)
i = i_old
j = j_old
@@ -139,8 +181,9 @@ function sortIceFloesInGrid!(simulation::Simulation, grid::Any; verbose=false)
i_t = max(min(i_old + i_rel, nx), 1)
j_t = max(min(j_old + j_rel, ny), 1)
- if isPointInCell(grid, i_t, j_t,
- simulation.ice_floes[idx].lin_pos)
+ @inbounds if isPointInCell(grid, i_t, j_t,
+ simulation.ice_floes[idx].lin_pos,
+ sw, se, ne, nw)
i = i_t
j = j_t
found = true
@@ -166,17 +209,18 @@ function sortIceFloesInGrid!(simulation::Simulation, grid::Any; verbose=false)
# add cell to ice floe
if typeof(grid) == Ocean
- simulation.ice_floes[idx].ocean_grid_pos = [i, j]
+ @inbounds simulation.ice_floes[idx].ocean_grid_pos = [i, j]
elseif typeof(grid) == Atmosphere
- simulation.ice_floes[idx].atmosphere_grid_pos = [i, j]
+ @inbounds simulation.ice_floes[idx].atmosphere_grid_pos = [i, j]
else
error("grid type not understood.")
end
end
# add ice floe to cell
- push!(grid.ice_floe_list[i, j], idx)
+ @inbounds push!(grid.ice_floe_list[i, j], idx)
end
+ nothing
end
export findCellContainingPoint
@@ -191,16 +235,23 @@ found the function returns `(0,0)`.
# Arguments
* `grid::Any`: grid object containing ocean or atmosphere data.
-* `point::Array{float, 1}`: two-dimensional vector of point to check.
+* `point::Vector{Float64}`: 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(grid::Any, point::Array{float, 1};
+function findCellContainingPoint(grid::Any, point::Vector{Float64};
method::String="Conformal")
+ sw = Vector{Float64}(2)
+ se = Vector{Float64}(2)
+ ne = Vector{Float64}(2)
+ nw = Vector{Float64}(2)
+
for i=1:size(grid.xh, 1)
for j=1:size(grid.yh, 2)
- if isPointInCell(grid, i, j, point, method=method)
+ if isPointInCell(grid, i, j, point,
+ sw, se, ne, nw,
+ method=method)
return i, j
end
end
@@ -217,7 +268,7 @@ This function is a wrapper for `getCellCornerCoordinates()` and
`conformalQuadrilateralCoordinates()`.
"""
function getNonDimensionalCellCoordinates(grid::Any, i::Int, j::Int,
- point::Array{float, 1})
+ point::Vector{Float64})
sw, se, ne, nw = getCellCornerCoordinates(grid.xq, grid.yq, i, j)
return conformalQuadrilateralCoordinates(sw, se, ne, nw, point)
@@ -230,10 +281,22 @@ 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`.
"""
-function isPointInCell(grid::Any, i::Int, j::Int, point::Array{float, 1};
+function isPointInCell(grid::Any, i::Int, j::Int, point::Vector{Float64},
+ sw::Vector{Float64} = Vector{Float64}(2),
+ se::Vector{Float64} = Vector{Float64}(2),
+ ne::Vector{Float64} = Vector{Float64}(2),
+ nw::Vector{Float64} = Vector{Float64}(2);
method::String="Conformal")
- sw, se, ne, nw = getCellCornerCoordinates(grid.xq, grid.yq, i, j)
+ #sw, se, ne, nw = getCellCornerCoordinates(grid.xq, grid.yq, i, j)
+ @views sw[1] = grid.xq[ i, j]
+ @views sw[2] = grid.yq[ i, j]
+ @views se[1] = grid.xq[i+1, j]
+ @views se[2] = grid.yq[i+1, j]
+ @views ne[1] = grid.xq[i+1, j+1]
+ @views ne[2] = grid.yq[i+1, j+1]
+ @views nw[1] = grid.xq[ i, j+1]
+ @views nw[2] = grid.yq[ i, j+1]
if method == "Area"
if areaOfQuadrilateral(sw, se, ne, nw) ≈
@@ -275,7 +338,7 @@ south-east corner, north-east corner, north-west corner).
@inline function getCellCornerCoordinates(xq::Array{Float64, 2},
yq::Array{Float64, 2},
i::Int, j::Int)
- return Float64[xq[ i, j], yq[ i, j]],
+ @inbounds return Float64[xq[ i, j], yq[ i, j]],
Float64[xq[i+1, j], yq[i+1, j]],
Float64[xq[i+1, j+1], yq[i+1, j+1]],
Float64[xq[ i, j+1], yq[ i, j+1]]
@@ -300,9 +363,9 @@ end
export areaOfTriangle
"Returns the area of an triangle with corner coordinates `a`, `b`, and `c`."
-function areaOfTriangle(a::Array{float, 1},
- b::Array{float, 1},
- c::Array{float, 1})
+function areaOfTriangle(a::Vector{Float64},
+ b::Vector{Float64},
+ c::Vector{Float64})
return abs(
(a[1]*(b[2] - c[2]) +
b[1]*(c[2] - a[2]) +
@@ -317,10 +380,10 @@ Returns the area of a quadrilateral with corner coordinates `a`, `b`, `c`, and
true for `b` and `d`. This is true if the four corners are passed as arguments
in a "clockwise" or "counter-clockwise" manner.
"""
-function areaOfQuadrilateral(a::Array{float, 1},
- b::Array{float, 1},
- c::Array{float, 1},
- d::Array{float, 1})
+function areaOfQuadrilateral(a::Vector{Float64},
+ b::Vector{Float64},
+ c::Vector{Float64},
+ d::Vector{Float64})
return areaOfTriangle(a, b, c) + areaOfTriangle(c, d, a)
end
@@ -331,11 +394,11 @@ within a quadrilateral with corner coordinates `A`, `B`, `C`, and `D`.
Points must be ordered in counter-clockwise order, starting from south-west
corner.
"""
-function conformalQuadrilateralCoordinates(A::Array{float, 1},
- B::Array{float, 1},
- C::Array{float, 1},
- D::Array{float, 1},
- p::Array{float, 1})
+function conformalQuadrilateralCoordinates(A::Vector{Float64},
+ B::Vector{Float64},
+ C::Vector{Float64},
+ D::Vector{Float64},
+ p::Vector{Float64})
if !(A[1] < B[1] && B[2] < C[2] && C[1] > D[1])
error("corner coordinates are not passed in the correct order")
@@ -404,13 +467,13 @@ function findEmptyPositionInGridCell(simulation::Simulation,
grid::Any,
i::Int,
j::Int,
- r::float;
+ r::Float64;
n_iter::Int = 10,
seed::Int = 1,
verbose::Bool = false)
overlap_found = false
i_iter = 0
- pos = [NaN NaN]
+ pos = [NaN, NaN]
nx, ny = size(grid.xh)
@@ -421,8 +484,7 @@ function findEmptyPositionInGridCell(simulation::Simulation,
# generate random candidate position
x_tilde = Base.Random.rand()
y_tilde = Base.Random.rand()
- pos = [bilinearInterpolation(grid.xq, x_tilde, y_tilde, i, j)
- bilinearInterpolation(grid.yq, x_tilde, y_tilde, i, j)]
+ bilinearInterpolation!(pos, grid.xq, grid.yq, x_tilde, y_tilde, i, j)
if verbose
info("trying poisition $pos in cell $i,$j")
end
@@ -496,4 +558,5 @@ function copyGridSortingInfo!(ocean::Ocean, atmosphere::Atmosphere,
icefloe.atmosphere_grid_pos = deepcopy(icefloe.ocean_grid_pos)
end
atmosphere.ice_floe_list = deepcopy(ocean.ice_floe_list)
+ nothing
end
diff --git a/src/icefloe.jl b/src/icefloe.jl
@@ -1,4 +1,5 @@
## Manage icefloes in the model
+import PyPlot
export addIceFloeCylindrical!
"""
@@ -7,39 +8,39 @@ Adds a grain to the simulation. Example:
SeaIce.addIceFloeCylindrical([1.0, 2.0, 3.0], 1.0)
"""
function addIceFloeCylindrical!(simulation::Simulation,
- lin_pos::vector,
- contact_radius::float,
- thickness::float;
+ lin_pos::Vector{Float64},
+ contact_radius::Float64,
+ thickness::Float64;
areal_radius = false,
- lin_vel::vector = [0., 0.],
- lin_acc::vector = [0., 0.],
- force::vector = [0., 0.],
- ang_pos::float = 0.,
- ang_vel::float = 0.,
- ang_acc::float = 0.,
- torque::float = 0.,
- density::float = 934.,
- contact_stiffness_normal::float = 1e7,
- contact_stiffness_tangential::float = 0.,
- contact_viscosity_normal::float = 0.,
- contact_viscosity_tangential::float = 0.,
- contact_static_friction::float = 0.4,
- contact_dynamic_friction::float = 0.4,
- youngs_modulus::float = 2e7,
- #youngs_modulus::float = 2e9, # Hopkins 2004
- poissons_ratio::float = 0.185, # Hopkins 2004
- #tensile_strength::float = 500e3, # Hopkins2004
- tensile_strength::float = 0.,
- tensile_heal_rate::float = 0.,
- compressive_strength_prefactor::float = 1285e3,
+ lin_vel::Vector{Float64} = [0., 0.],
+ lin_acc::Vector{Float64} = [0., 0.],
+ force::Vector{Float64} = [0., 0.],
+ ang_pos::Float64 = 0.,
+ ang_vel::Float64 = 0.,
+ ang_acc::Float64 = 0.,
+ torque::Float64 = 0.,
+ density::Float64 = 934.,
+ contact_stiffness_normal::Float64 = 1e7,
+ contact_stiffness_tangential::Float64 = 0.,
+ contact_viscosity_normal::Float64 = 0.,
+ contact_viscosity_tangential::Float64 = 0.,
+ contact_static_friction::Float64 = 0.4,
+ contact_dynamic_friction::Float64 = 0.4,
+ youngs_modulus::Float64 = 2e7,
+ #youngs_modulus::Float64 = 2e9, # Hopkins 2004
+ poissons_ratio::Float64 = 0.185, # Hopkins 2004
+ #tensile_strength::Float64 = 500e3, # Hopkins2004
+ tensile_strength::Float64 = 0.,
+ tensile_heal_rate::Float64 = 0.,
+ compressive_strength_prefactor::Float64 = 1285e3,
# Hopkins 2004
- ocean_drag_coeff_vert::float = 0.85, # H&C 2011
- ocean_drag_coeff_horiz::float = 5e-4, # H&C 2011
- atmosphere_drag_coeff_vert::float = 0.4,
+ ocean_drag_coeff_vert::Float64 = 0.85, # H&C 2011
+ ocean_drag_coeff_horiz::Float64 = 5e-4, # H&C 2011
+ atmosphere_drag_coeff_vert::Float64 = 0.4,
# H&C 2011
- atmosphere_drag_coeff_horiz::float = 2.5e-4,
+ atmosphere_drag_coeff_horiz::Float64 = 2.5e-4,
# H&C2011
- pressure::float = 0.,
+ pressure::Float64 = 0.,
fixed::Bool = false,
rotating::Bool = true,
enabled::Bool = true,
@@ -47,9 +48,9 @@ function addIceFloeCylindrical!(simulation::Simulation,
ocean_grid_pos::Array{Int, 1} = [0, 0],
atmosphere_grid_pos::Array{Int, 1} = [0, 0],
n_contacts::Int = 0,
- granular_stress::vector = [0., 0.],
- ocean_stress::vector = [0., 0.],
- atmosphere_stress::vector = [0., 0.])
+ granular_stress::Vector{Float64} = [0., 0.],
+ ocean_stress::Vector{Float64} = [0., 0.],
+ atmosphere_stress::Vector{Float64} = [0., 0.])
# Check input values
if length(lin_pos) != 2
@@ -77,8 +78,8 @@ function addIceFloeCylindrical!(simulation::Simulation,
contacts::Array{Int, 1} = zeros(Int, simulation.Nc_max)
contact_parallel_displacement =
- Array{Array{Float64, 1}, 1}(simulation.Nc_max)
- contact_age::Array{Float64, 1} = zeros(Float64, simulation.Nc_max)
+ Vector{Vector{Float64}}(simulation.Nc_max)
+ contact_age::Vector{Float64} = zeros(Float64, simulation.Nc_max)
for i=1:simulation.Nc_max
contact_parallel_displacement[i] = zeros(2)
end
@@ -151,6 +152,7 @@ function addIceFloeCylindrical!(simulation::Simulation,
# Add to simulation object
addIceFloe!(simulation, icefloe, verbose)
+ nothing
end
export iceFloeCircumreference
@@ -318,56 +320,62 @@ function convertIceFloeDataToArrays(simulation::Simulation)
end
function deleteIceFloeArrays!(ifarr::IceFloeArrays)
- ifarr.density = 0
-
- ifarr.thickness = 0
- ifarr.contact_radius = 0
- ifarr.areal_radius = 0
- ifarr.circumreference = 0
- ifarr.horizontal_surface_area = 0
- ifarr.side_surface_area = 0
- ifarr.volume = 0
- ifarr.mass = 0
- ifarr.moment_of_inertia = 0
-
- ifarr.lin_pos = 0
- ifarr.lin_vel = 0
- ifarr.lin_acc = 0
- ifarr.force = 0
-
- ifarr.ang_pos = 0
- ifarr.ang_vel = 0
- ifarr.ang_acc = 0
- ifarr.torque = 0
-
- ifarr.fixed = 0
- ifarr.rotating = 0
- ifarr.enabled = 0
-
- ifarr.contact_stiffness_normal = 0
- ifarr.contact_stiffness_tangential = 0
- ifarr.contact_viscosity_normal = 0
- ifarr.contact_viscosity_tangential = 0
- ifarr.contact_static_friction = 0
- ifarr.contact_dynamic_friction = 0
-
- ifarr.youngs_modulus = 0
- ifarr.poissons_ratio = 0
- ifarr.tensile_strength = 0
- ifarr.compressive_strength_prefactor = 0
-
- ifarr.ocean_drag_coeff_vert = 0
- ifarr.ocean_drag_coeff_horiz = 0
- ifarr.atmosphere_drag_coeff_vert = 0
- ifarr.atmosphere_drag_coeff_horiz = 0
-
- ifarr.pressure = 0
- ifarr.n_contacts = 0
-
- ifarr.granular_stress = 0
- ifarr.ocean_stress = 0
- ifarr.atmosphere_stress = 0
+ f1 = zeros(1)
+ f2 = zeros(1,1)
+ i1 = zeros(Int, 1)
+
+ ifarr.density = f1
+
+ ifarr.thickness = f1
+ ifarr.contact_radius = f1
+ ifarr.areal_radius = f1
+ ifarr.circumreference = f1
+ ifarr.horizontal_surface_area = f1
+ ifarr.side_surface_area = f1
+ ifarr.volume = f1
+ ifarr.mass = f1
+ ifarr.moment_of_inertia = f1
+
+ ifarr.lin_pos = f2
+ ifarr.lin_vel = f2
+ ifarr.lin_acc = f2
+ ifarr.force = f2
+
+ ifarr.ang_pos = f2
+ ifarr.ang_vel = f2
+ ifarr.ang_acc = f2
+ ifarr.torque = f2
+
+ ifarr.fixed = i1
+ ifarr.rotating = i1
+ ifarr.enabled = i1
+
+ ifarr.contact_stiffness_normal = f1
+ ifarr.contact_stiffness_tangential = f1
+ ifarr.contact_viscosity_normal = f1
+ ifarr.contact_viscosity_tangential = f1
+ ifarr.contact_static_friction = f1
+ ifarr.contact_dynamic_friction = f1
+
+ ifarr.youngs_modulus = f1
+ ifarr.poissons_ratio = f1
+ ifarr.tensile_strength = f1
+ ifarr.compressive_strength_prefactor = f1
+
+ ifarr.ocean_drag_coeff_vert = f1
+ ifarr.ocean_drag_coeff_horiz = f1
+ ifarr.atmosphere_drag_coeff_vert = f1
+ ifarr.atmosphere_drag_coeff_horiz = f1
+
+ ifarr.pressure = f1
+ ifarr.n_contacts = i1
+
+ ifarr.granular_stress = f2
+ ifarr.ocean_stress = f2
+ ifarr.atmosphere_stress = f2
+
gc()
+ nothing
end
export printIceFloeInfo
@@ -427,6 +435,7 @@ function printIceFloeInfo(f::IceFloeCylindrical)
println(" granular_stress: $(f.granular_stress) Pa")
println(" ocean_stress: $(f.ocean_stress) Pa")
println(" atmosphere_stress: $(f.atmosphere_stress) Pa")
+ nothing
end
export iceFloeKineticTranslationalEnergy
@@ -535,4 +544,64 @@ function compareIceFloes(if1::IceFloeCylindrical, if2::IceFloeCylindrical)
Base.Test.@test if1.granular_stress ≈ if2.granular_stress
Base.Test.@test if1.ocean_stress ≈ if2.ocean_stress
Base.Test.@test if1.atmosphere_stress ≈ if2.atmosphere_stress
+ nothing
+end
+
+export plotIceFloeSizeDistribution
+"""
+ plotIceFloeSizeDistribution(simulation, [filename_postfix], [nbins],
+ [size_type], [figsize], [filetype])
+
+Plot the ice-floe size distribution as a histogram and save it to the disk. The
+plot is saved accoring to the simulation id, the optional `filename_postfix`
+string, and the `filetype`, and is written to the current folder.
+
+# Arguments
+* `simulation::Simulation`: the simulation object containing the ice floes.
+* `filename_postfix::String`: optional string for the output filename.
+* `nbins::Int`: number of bins in the histogram (default = 12).
+* `size_type::String`: specify whether to use the `contact` or `areal` radius
+ for the ice-floe size. The default is `contact`.
+* `figsize::Tuple`: the output figure size in inches (default = (6,4).
+* `filetype::String`: the output file type (default = "png").
+* `verbose::String`: show output file as info message in stdout (default =
+ true).
+* `skip_fixed::Bool`: ommit ice floes that are fixed in space from the size
+ distribution (default = true).
+* `logy::Bool`: plot y-axis in log scale.
+"""
+function plotIceFloeSizeDistribution(simulation::Simulation;
+ filename_postfix::String = "",
+ nbins::Int=12,
+ size_type::String = "contact",
+ figsize::Tuple = (6,4),
+ filetype::String = "png",
+ verbose::Bool = true,
+ skip_fixed::Bool = true,
+ log_y::Bool = true)
+
+ diameters = Float64[]
+ for i=1:length(simulation.ice_floes)
+ if simulation.ice_floes[i].fixed && skip_fixed
+ continue
+ end
+ if size_type == "contact"
+ push!(diameters, simulation.ice_floes[i].contact_radius*2.)
+ elseif size_type == "areal"
+ push!(diameters, simulation.ice_floes[i].areal_radius*2.)
+ else
+ error("size_type '$size_type' not understood")
+ end
+ end
+ PyPlot.pygui(false)
+ PyPlot.figure(figsize=figsize)
+ PyPlot.plt[:hist](diameters, nbins, log=log_y)
+ PyPlot.xlabel("Diameter [m]")
+ PyPlot.ylabel("Count [-]")
+ filename = string(simulation.id * filename_postfix *
+ "-ice-floe-size-distribution." * filetype)
+ PyPlot.savefig(filename)
+ if verbose
+ info(filename)
+ end
end
diff --git a/src/interaction.jl b/src/interaction.jl
@@ -16,26 +16,26 @@ function interact!(simulation::Simulation)
continue
end
- """
- if norm(simulation.ice_floes[i].lin_pos -
+ #=if norm(simulation.ice_floes[i].lin_pos -
simulation.ice_floes[j].lin_pos) -
- (simulation.ice_floes[i].contact_radius +
- simulation.ice_floes[j].contact_radius) > 0.
+ (simulation.ice_floes[i].contact_radius +
+ simulation.ice_floes[j].contact_radius) > 0.
simulation.ice_floes[i].contacts[ic] = 0 # remove contact
simulation.ice_floes[i].n_contacts -= 1
simulation.ice_floes[j].n_contacts -= 1
- else
- """
+ else=#
interactIceFloes!(simulation, i, j, ic)
#end
end
end
for i=1:length(simulation.ice_floes)
- simulation.ice_floes[i].granular_stress = simulation.ice_floes[i].force/
+ @inbounds simulation.ice_floes[i].granular_stress =
+ simulation.ice_floes[i].force/
simulation.ice_floes[i].horizontal_surface_area
end
+ nothing
end
export interactIceFloes!
@@ -197,12 +197,5 @@ function interactIceFloes!(simulation::Simulation, i::Int, j::Int, ic::Int)
force_n/simulation.ice_floes[i].side_surface_area;
simulation.ice_floes[j].pressure +=
force_n/simulation.ice_floes[j].side_surface_area;
-end
-
-function harmonicMean(a::Any, b::Any)
- if a ≈ 0. && b ≈ 0
- return 0.
- else
- return 2.*a*b/(a + b)
- end
+ nothing
end
diff --git a/src/io.jl b/src/io.jl
@@ -41,6 +41,7 @@ function writeSimulation(simulation::Simulation;
info("simulation written to $filename")
end
end
+ nothing
end
export readSimulation
@@ -56,6 +57,7 @@ function readSimulation(filename::String="";
warn("Package JLD not found. Simulation save/read not supported. " *
"Please install JLD and its " *
"requirements with `Pkg.add(\"JLD\")`.")
+ nothing
else
if verbose
info("reading simulation from $filename")
@@ -85,6 +87,7 @@ function writeSimulationStatus(simulation::Simulation;
if verbose
info("wrote status to $filename")
end
+ nothing
end
export readSimulationStatus
@@ -197,6 +200,7 @@ function status(folder::String=".";
repeat = false
end
end
+ nothing
end
export writeVTK
@@ -240,6 +244,7 @@ function writeVTK(simulation::Simulation;
simulation.file_number)
writeGridVTK(simulation.atmosphere, filename, verbose=verbose)
end
+ nothing
end
export writeIceFloeVTK
@@ -343,9 +348,8 @@ function writeIceFloeVTK(simulation::Simulation,
outfiles = WriteVTK.vtk_save(vtkfile)
if verbose
info("Output file: " * outfiles[1])
- else
- return nothing
end
+ nothing
end
export writeIceFloeInteractionVTK
@@ -363,14 +367,14 @@ function writeIceFloeInteractionVTK(simulation::Simulation,
i1 = Int64[]
i2 = Int64[]
- inter_particle_vector = Array{float, 1}[]
- force = float[]
- effective_radius = float[]
- contact_area = float[]
- contact_stiffness = float[]
- tensile_stress = float[]
- shear_displacement = Array{float, 1}[]
- contact_age = float[]
+ inter_particle_vector = Vector{Float64}[]
+ force = Float64[]
+ effective_radius = Float64[]
+ contact_area = Float64[]
+ contact_stiffness = Float64[]
+ tensile_stress = Float64[]
+ shear_displacement = Vector{Float64}[]
+ contact_age = Float64[]
for i=1:length(simulation.ice_floes)
for ic=1:simulation.Nc_max
if simulation.ice_floes[i].contacts[ic] > 0
@@ -457,8 +461,8 @@ function writeIceFloeInteractionVTK(simulation::Simulation,
"Name=\"Shear displacement [m]\" " *
"NumberOfComponents=\"3\" format=\"ascii\">\n")
for i=1:length(i1)
- write(f, "$(shear_displacement[i][1]) ")
- write(f, "$(shear_displacement[i][2]) ")
+ @inbounds write(f, "$(shear_displacement[i][1]) ")
+ @inbounds write(f, "$(shear_displacement[i][2]) ")
write(f, "0.0 ")
end
write(f, "\n")
@@ -467,7 +471,7 @@ function writeIceFloeInteractionVTK(simulation::Simulation,
write(f, " <DataArray type=\"Float32\" Name=\"Force [N]\" " *
"NumberOfComponents=\"1\" format=\"ascii\">\n")
for i=1:length(i1)
- write(f, "$(force[i]) ")
+ @inbounds write(f, "$(force[i]) ")
end
write(f, "\n")
write(f, " </DataArray>\n")
@@ -476,7 +480,7 @@ function writeIceFloeInteractionVTK(simulation::Simulation,
"Name=\"Effective radius [m]\" " *
"NumberOfComponents=\"1\" format=\"ascii\">\n")
for i=1:length(i1)
- write(f, "$(effective_radius[i]) ")
+ @inbounds write(f, "$(effective_radius[i]) ")
end
write(f, "\n")
write(f, " </DataArray>\n")
@@ -485,7 +489,7 @@ function writeIceFloeInteractionVTK(simulation::Simulation,
"Name=\"Contact area [m^2]\" " *
"NumberOfComponents=\"1\" format=\"ascii\">\n")
for i=1:length(i1)
- write(f, "$(contact_area[i]) ")
+ @inbounds write(f, "$(contact_area[i]) ")
end
write(f, "\n")
write(f, " </DataArray>\n")
@@ -494,7 +498,7 @@ function writeIceFloeInteractionVTK(simulation::Simulation,
"Name=\"Contact stiffness [N/m]\" " *
"NumberOfComponents=\"1\" format=\"ascii\">\n")
for i=1:length(i1)
- write(f, "$(contact_stiffness[i]) ")
+ @inbounds write(f, "$(contact_stiffness[i]) ")
end
write(f, "\n")
write(f, " </DataArray>\n")
@@ -503,7 +507,7 @@ function writeIceFloeInteractionVTK(simulation::Simulation,
"Name=\"Tensile stress [Pa]\" " *
"NumberOfComponents=\"1\" format=\"ascii\">\n")
for i=1:length(i1)
- write(f, "$(tensile_stress[i]) ")
+ @inbounds write(f, "$(tensile_stress[i]) ")
end
write(f, "\n")
write(f, " </DataArray>\n")
@@ -512,7 +516,7 @@ function writeIceFloeInteractionVTK(simulation::Simulation,
"Name=\"Contact age [s]\" NumberOfComponents=\"1\"
format=\"ascii\">\n")
for i=1:length(i1)
- write(f, "$(contact_age[i]) ")
+ @inbounds write(f, "$(contact_age[i]) ")
end
write(f, "\n")
write(f, " </DataArray>\n")
@@ -525,7 +529,7 @@ function writeIceFloeInteractionVTK(simulation::Simulation,
write(f, " <DataArray type=\"Float32\" Name=\"Points\" " *
"NumberOfComponents=\"3\" format=\"ascii\">\n")
for i in simulation.ice_floes
- write(f, "$(i.lin_pos[1]) $(i.lin_pos[2]) 0.0 ")
+ @inbounds write(f, "$(i.lin_pos[1]) $(i.lin_pos[2]) 0.0 ")
end
write(f, "\n")
write(f, " </DataArray>\n")
@@ -544,7 +548,7 @@ function writeIceFloeInteractionVTK(simulation::Simulation,
write(f, " <DataArray type=\"Int64\" Name=\"connectivity\" " *
"format=\"ascii\">\n")
for i=1:length(i1)
- write(f, "$(i1[i] - 1) $(i2[i] - 1) ")
+ @inbounds write(f, "$(i1[i] - 1) $(i2[i] - 1) ")
end
write(f, "\n")
write(f, " </DataArray>\n")
@@ -554,7 +558,7 @@ function writeIceFloeInteractionVTK(simulation::Simulation,
write(f, " <DataArray type=\"Int64\" Name=\"offsets\" " *
"format=\"ascii\">\n")
for i=1:length(i1)
- write(f, "$((i - 1)*2 + 2) ")
+ @inbounds write(f, "$((i - 1)*2 + 2) ")
end
write(f, "\n")
write(f, " </DataArray>\n")
@@ -580,6 +584,7 @@ function writeIceFloeInteractionVTK(simulation::Simulation,
write(f, " </PolyData>\n")
write(f, "</VTKFile>\n")
end
+ nothing
end
export writeOceanVTK
@@ -599,12 +604,12 @@ function writeGridVTK(grid::Any,
zq = similar(grid.u[:,:,:,1])
for iz=1:size(xq, 3)
- xq[:,:,iz] = grid.xq
- yq[:,:,iz] = grid.yq
+ @inbounds xq[:,:,iz] = grid.xq
+ @inbounds yq[:,:,iz] = grid.yq
end
for ix=1:size(xq, 1)
for iy=1:size(xq, 2)
- zq[ix,iy,:] = grid.zl
+ @inbounds zq[ix,iy,:] = grid.zl
end
end
@@ -620,8 +625,8 @@ function writeGridVTK(grid::Any,
for ix=1:size(xq, 1)
for iy=1:size(xq, 2)
for iz=1:size(xq, 3)
- vel[1, ix, iy, iz] = grid.u[ix, iy, iz, 1]
- vel[2, ix, iy, iz] = grid.v[ix, iy, iz, 1]
+ @inbounds vel[1, ix, iy, iz] = grid.u[ix, iy, iz, 1]
+ @inbounds vel[2, ix, iy, iz] = grid.v[ix, iy, iz, 1]
end
end
end
@@ -638,9 +643,8 @@ function writeGridVTK(grid::Any,
outfiles = WriteVTK.vtk_save(vtkfile)
if verbose
info("Output file: " * outfiles[1])
- else
- return nothing
end
+ nothing
end
export writeParaviewStateFile
@@ -18176,6 +18180,7 @@ function writeParaviewStateFile(simulation::Simulation;
if verbose
info("Output file: " * filename)
end
+ nothing
end
export removeSimulationFiles
@@ -18192,4 +18197,5 @@ function removeSimulationFiles(simulation::Simulation; folder::String=".")
run(`bash -c "rm -rf $(folder)/$(simulation.id).status.txt"`)
run(`bash -c "rm -rf $(folder)/$(simulation.id).*.jld"`)
run(`bash -c "rm -rf $(folder)"`)
+ nothing
end
diff --git a/src/ocean.jl b/src/ocean.jl
@@ -71,17 +71,17 @@ or `########.ocean_month.nc`) from disk and return time stamps, velocity fields,
layer thicknesses, interface heights, and vertical coordinates.
# Returns
-* `time::Array{Float, 1}`: Time [s]
-* `u::Array{Float, 2}`: Cell corner zonal velocity [m/s],
+* `time::Vector{Float64}`: Time [s]
+* `u::Array{Float64, 2}`: Cell corner zonal velocity [m/s],
dimensions correspond to placement in `[xq, yq, zl, time]`
-* `v::Array{Float, 2}`: Cell corner meridional velocity [m/s],
+* `v::Array{Float64, 2}`: Cell corner meridional velocity [m/s],
dimensions correspond to placement in `[xq, yq, zl, time]`
* `h::Array{Float64, 2}`: layer thickness [m], dimensions correspond to
placement in `[xh, yh, zl, time]`
* `e::Array{Float64, 2}`: interface height relative to mean sea level [m],
dimensions correspond to placement in `[xh, yh, zi, time]`
-* `zl::Array{Float64, 1}`: layer target potential density [kg m^-3]
-* `zi::Array{Float64, 1}`: interface target potential density [kg m^-3]
+* `zl::Vector{Float64}`: layer target potential density [kg m^-3]
+* `zi::Vector{Float64}`: interface target potential density [kg m^-3]
"""
function readOceanStateNetCDF(filename::String)
@@ -93,13 +93,13 @@ function readOceanStateNetCDF(filename::String)
v_staggered = convert(Array{Float64, 4}, NetCDF.ncread(filename, "v"))
u, v = interpolateOceanVelocitiesToCorners(u_staggered, v_staggered)
- time = convert(Array{Float64, 1},
+ time = convert(Vector{Float64},
NetCDF.ncread(filename, "time")*24.*60.*60.)
h = convert(Array{Float64, 4}, NetCDF.ncread(filename, "h"))
e = convert(Array{Float64, 4}, NetCDF.ncread(filename, "e"))
- zl = convert(Array{Float64, 1}, NetCDF.ncread(filename, "zl"))
- zi = convert(Array{Float64, 1}, NetCDF.ncread(filename, "zi"))
+ zl = convert(Vector{Float64}, NetCDF.ncread(filename, "zl"))
+ zi = convert(Vector{Float64}, NetCDF.ncread(filename, "zi"))
return time, u, v, h, e, zl, zi
end
@@ -172,7 +172,7 @@ function performs linear interpolation between time steps to get the approximate
ocean state at any point in time. If the `Ocean` data set only contains a
single time step, values from that time are returned.
"""
-function interpolateOceanState(ocean::Ocean, t::float)
+function interpolateOceanState(ocean::Ocean, t::Float64)
if length(ocean.time) == 1
return ocean.u, ocean.v, ocean.h, ocean.e
elseif t < ocean.time[1] || t > ocean.time[end]
@@ -213,9 +213,9 @@ the ocean spanning `z<0.`. Vertical indexing starts with `k=0` at the sea
surface, and increases downwards.
"""
function createRegularOceanGrid(n::Array{Int, 1},
- L::Array{float, 1};
- origo::Array{float, 1} = zeros(2),
- time::Array{float, 1} = zeros(1),
+ L::Vector{Float64};
+ origo::Vector{Float64} = zeros(2),
+ time::Vector{Float64} = zeros(1),
name::String = "unnamed")
xq = repmat(linspace(origo[1], L[1], n[1] + 1), 1, n[2] + 1)
@@ -253,6 +253,11 @@ function addOceanDrag!(simulation::Simulation)
end
u, v, h, e = interpolateOceanState(simulation.ocean, simulation.time)
+ uv_interp = Vector{Float64}(2)
+ sw = Vector{Float64}(2)
+ se = Vector{Float64}(2)
+ ne = Vector{Float64}(2)
+ nw = Vector{Float64}(2)
for ice_floe in simulation.ice_floes
@@ -274,15 +279,13 @@ function addOceanDrag!(simulation::Simulation)
""")
end
- applyOceanDragToIceFloe!(ice_floe,
- bilinearInterpolation(u, x_tilde, y_tilde,
- i, j, k, 1),
- bilinearInterpolation(v, x_tilde, y_tilde,
- i, j, k, 1))
+ bilinearInterpolation!(uv_interp, u, v, x_tilde, y_tilde, i, j, k, 1)
+ applyOceanDragToIceFloe!(ice_floe, uv_interp[1], uv_interp[2])
applyOceanVorticityToIceFloe!(ice_floe,
curl(simulation.ocean, x_tilde, y_tilde,
- i, j, k, 1))
+ i, j, k, 1, sw, se, ne, nw))
end
+ nothing
end
export applyOceanDragToIceFloe!
@@ -291,7 +294,7 @@ Add Stokes-type drag from velocity difference between ocean and a single ice
floe.
"""
function applyOceanDragToIceFloe!(ice_floe::IceFloeCylindrical,
- u::float, v::float)
+ u::Float64, v::Float64)
freeboard = .1*ice_floe.thickness # height above water
rho_o = 1000. # ocean density
draft = ice_floe.thickness - freeboard # height of submerged thickness
@@ -304,6 +307,7 @@ function applyOceanDragToIceFloe!(ice_floe::IceFloeCylindrical,
ice_floe.force += drag_force
ice_floe.ocean_stress = drag_force/ice_floe.horizontal_surface_area
+ nothing
end
export applyOceanVorticityToIceFloe!
@@ -313,7 +317,7 @@ single ice floe. See Eq. 9.28 in "Introduction to Fluid Mechanics" by Nakayama
and Boucher, 1999.
"""
function applyOceanVorticityToIceFloe!(ice_floe::IceFloeCylindrical,
- ocean_curl::float)
+ ocean_curl::Float64)
freeboard = .1*ice_floe.thickness # height above water
rho_o = 1000. # ocean density
draft = ice_floe.thickness - freeboard # height of submerged thickness
@@ -323,6 +327,7 @@ function applyOceanVorticityToIceFloe!(ice_floe::IceFloeCylindrical,
(ice_floe.areal_radius/5.*ice_floe.ocean_drag_coeff_horiz +
draft*ice_floe.ocean_drag_coeff_vert)*
abs(.5*ocean_curl - ice_floe.ang_vel)*(.5*ocean_curl - ice_floe.ang_vel)
+ nothing
end
export compareOceans
@@ -353,4 +358,5 @@ function compareOceans(ocean1::Ocean, ocean2::Ocean)
if isassigned(ocean1.ice_floe_list, 1)
Base.Test.@test ocean1.ice_floe_list == ocean2.ice_floe_list
end
+ nothing
end
diff --git a/src/simulation.jl b/src/simulation.jl
@@ -93,7 +93,7 @@ function run!(simulation::Simulation;
simulation.time_total += simulation.time_step
end
- checkTimeParameters(simulation)
+ checkTimeParameters(simulation, single_step=single_step)
# if both are enabled, check if the atmosphere grid spatial geometry is
# identical to the ocean grid
@@ -175,7 +175,7 @@ function run!(simulation::Simulation;
incrementCurrentTime!(simulation, simulation.time_step)
if single_step
- return
+ return nothing
end
end
@@ -191,6 +191,8 @@ function run!(simulation::Simulation;
reportSimulationTimeToStdout(simulation)
println()
end
+ gc()
+ nothing
end
export addIceFloe!
@@ -210,16 +212,17 @@ function addIceFloe!(simulation::Simulation,
if verbose
info("Added IceFloe $(length(simulation.ice_floes))")
end
+ nothing
end
export disableIceFloe!
"Disable ice floe with index `i` in the `simulation` object."
-function disableIceFloe!(simulation::Simulation, i::Integer)
+function disableIceFloe!(simulation::Simulation, i::Int)
if i < 1
error("Index must be greater than 0 (i = $i)")
end
-
simulation.ice_floes[i].enabled = false
+ nothing
end
export zeroForcesAndTorques!
@@ -230,6 +233,7 @@ function zeroForcesAndTorques!(simulation::Simulation)
icefloe.torque = 0.
icefloe.pressure = 0.
end
+ nothing
end
export reportSimulationTimeToStdout
@@ -237,6 +241,7 @@ export reportSimulationTimeToStdout
function reportSimulationTimeToStdout(simulation::Simulation)
print("\r t = ", simulation.time, '/', simulation.time_total,
" s ")
+ nothing
end
export compareSimulations
@@ -265,4 +270,46 @@ function compareSimulations(sim1::Simulation, sim2::Simulation)
compareAtmospheres(sim1.atmosphere, sim2.atmosphere)
Base.Test.@test sim1.Nc_max == sim2.Nc_max
+ nothing
+end
+
+export printMemoryUsage
+"""
+ printMemoryUsage(sim::Simulation)
+
+Shows the memory footprint of the simulation object.
+"""
+function printMemoryUsage(sim::Simulation)
+
+ reportMemory(sim, "sim")
+ println(" where:")
+
+ reportMemory(sim.ice_floes, " sim.ice_floes",
+ "(N=$(length(sim.ice_floes)))")
+
+ reportMemory(sim.ocean, " sim.ocean",
+ "($(size(sim.ocean.xh, 1))x" *
+ "$(size(sim.ocean.xh, 2))x" *
+ "$(size(sim.ocean.xh, 3)))")
+
+ reportMemory(sim.atmosphere, " sim.atmosphere",
+ "($(size(sim.atmosphere.xh, 1))x" *
+ "$(size(sim.atmosphere.xh, 2))x" *
+ "$(size(sim.atmosphere.xh, 3)))")
+ nothing
+end
+
+function reportMemory(variable, head::String, tail::String="")
+ bytes = Base.summarysize(variable)
+ if bytes < 10_000
+ size_str = @sprintf "%5d bytes" bytes
+ elseif bytes < 10_000 * 1024
+ size_str = @sprintf "%5d kb" bytes ÷ 1024
+ elseif bytes < 10_000 * 1024 * 1024
+ size_str = @sprintf "%5d Mb" bytes ÷ 1024 ÷ 1024
+ else
+ size_str = @sprintf "%5d Gb" bytes ÷ 1024 ÷ 1024 ÷ 1024
+ end
+ @printf("%-20s %s %s\n", head, size_str, tail)
+ nothing
end
diff --git a/src/temporal.jl b/src/temporal.jl
@@ -1,70 +1,77 @@
export setTotalTime!
"""
- setTotalTime!(simulation::Simulation, t::float)
+ setTotalTime!(simulation::Simulation, t::Float64)
Sets the total simulation time of the `simulation` object to `t`, with parameter
value checks.
"""
-function setTotalTime!(simulation::Simulation, t::float)
+function setTotalTime!(simulation::Simulation, t::Float64)
if t <= 0.0
error("Total time should be a positive value (t = $t s)")
end
simulation.time_total = t
+ nothing
end
export setCurrentTime!
"""
- setCurrentTime!(simulation::Simulation, t::float)
+ setCurrentTime!(simulation::Simulation, t::Float64)
Sets the current simulation time of the `simulation` object to `t`, with
parameter value checks.
"""
-function setCurrentTime!(simulation::Simulation, t::float)
+function setCurrentTime!(simulation::Simulation, t::Float64)
if t <= 0.0
error("Current time should be a positive value (t = $t s)")
end
simulation.time = t
+ nothing
end
export incrementCurrentTime!
"""
- incrementCurrentTime!(simulation::Simulation, t::float)
+ incrementCurrentTime!(simulation::Simulation, t::Float64)
Sets the current simulation time of the `simulation` object to `t`, with
parameter value checks.
"""
-function incrementCurrentTime!(simulation::Simulation, t::float)
+function incrementCurrentTime!(simulation::Simulation, t::Float64)
if t <= 0.0
error("Current time increment should be a positive value (t = $t s)")
end
simulation.time += t
simulation.file_time_since_output_file += t
+ nothing
end
export setOutputFileInterval!
"""
- setOutputFileInterval!(simulation::Simulation, t::float)
+ setOutputFileInterval!(simulation::Simulation, t::Float64)
Sets the simulation-time interval between output files are written to disk. If
this value is zero or negative, no output files will be written.
"""
-function setOutputFileInterval!(simulation::Simulation, t::float; verbose=true)
+function setOutputFileInterval!(simulation::Simulation, t::Float64;
+ verbose=true)
if t <= 0.0 && verbose
info("No output files will be written")
end
simulation.file_time_step = t
+ nothing
end
export disableOutputFiles!
"Disables the write of output files to disk during a simulation."
function disableOutputFiles!(simulation::Simulation)
simulation.file_time_step = 0.0
+ nothing
end
export checkTimeParameters
"Checks if simulation temporal parameters are of reasonable values."
-function checkTimeParameters(simulation::Simulation)
- if simulation.time_total <= 0.0 || simulation.time_total <= simulation.time
+function checkTimeParameters(simulation::Simulation; single_step::Bool=false)
+ if !single_step && (simulation.time_total <= 0.0 || simulation.time_total <=
+ simulation.time)
error("Total time should be positive and larger than current time.\n",
" t_total = ", simulation.time_total, " s, t = ", simulation.time,
" s")
@@ -72,6 +79,7 @@ function checkTimeParameters(simulation::Simulation)
if simulation.time_step <= 0.0
error("Time step should be positive (t = ", simulation.time_step, ")")
end
+ nothing
end
export findSmallestIceFloeMass
@@ -80,8 +88,8 @@ the optimal time step length."
function findSmallestIceFloeMass(simulation::Simulation)
m_min = 1e20
i_min = -1
- for i in length(simulation.ice_floes)
- icefloe = simulation.ice_floes[i]
+ for i=1:length(simulation.ice_floes)
+ @inbounds icefloe = simulation.ice_floes[i]
if icefloe.mass < m_min
m_min = icefloe.mass
i_min = i
@@ -98,9 +106,9 @@ function findLargestIceFloeStiffness(simulation::Simulation)
k_t_max = 0.
i_n_max = -1
i_t_max = -1
- for i in length(simulation.ice_floes)
+ for i=1:length(simulation.ice_floes)
- icefloe = simulation.ice_floes[i]
+ @inbounds icefloe = simulation.ice_floes[i]
if icefloe.youngs_modulus > 0.
k_n = icefloe.youngs_modulus*icefloe.thickness # A = h*r
@@ -129,7 +137,7 @@ Find the computational time step length suitable given the grain radii, contact
stiffnesses, and grain density. Uses the scheme by Radjaii et al. 2011.
"""
function setTimeStep!(simulation::Simulation;
- epsilon::float=0.07, verbose::Bool=true)
+ epsilon::Float64=0.07, verbose::Bool=true)
if length(simulation.ice_floes) < 1
error("At least 1 grain is needed to find the computational time step.")
@@ -147,4 +155,5 @@ function setTimeStep!(simulation::Simulation;
if verbose
info("Time step length t=", simulation.time_step, " s")
end
+ nothing
end
diff --git a/src/temporal_integration.jl b/src/temporal_integration.jl
@@ -36,6 +36,7 @@ function updateIceFloeKinematics!(simulation::Simulation;
else
error("Unknown integration method '$method'")
end
+ nothing
end
export updateIceFloeKinematicsTwoTermTaylor!
@@ -64,6 +65,7 @@ function updateIceFloeKinematicsTwoTermTaylor!(icefloe::IceFloeCylindrical,
icefloe.lin_vel += icefloe.lin_acc * simulation.time_step
icefloe.ang_vel += icefloe.ang_acc * simulation.time_step
+ nothing
end
export updateIceFloeKinematicsThreeTermTaylor!
@@ -109,4 +111,5 @@ function updateIceFloeKinematicsThreeTermTaylor!(icefloe::IceFloeCylindrical,
0.5 * d_lin_acc_dt * simulation.time_step^2.
icefloe.ang_vel += icefloe.ang_acc * simulation.time_step +
0.5 * d_ang_acc_dt * simulation.time_step^2.
+ nothing
end
diff --git a/src/util.jl b/src/util.jl
@@ -0,0 +1,43 @@
+#!/usr/bin/env julia
+
+export randpower
+"""
+ randpower([nvals], [distribution_power], [min_val], [max_val])
+
+Returns one or more random numbers from a power-law probability distribution.
+
+# Arguments
+* `dims::Any`: the dimensions of random values (default = 1)
+* `distribution_power::Number`: the distribution power (default = 1.)
+* `min_val::Number`: the lower bound of the distribution range (default = 0.)
+* `max_val::Number`: the upper bound of the distribution range (default = 1.)
+"""
+@inline function randpower(dims::Any = 1,
+ distribution_power::Number = 1.,
+ min_val::Number = 0.,
+ max_val::Number = 1.)
+
+ val = ((max_val^(distribution_power + 1.) -
+ min_val^(distribution_power + 1.))*Base.Random.rand(dims) +
+ min_val^(distribution_power + 1.)).^(1./(distribution_power + 1.))
+
+ if dims == 1
+ return val[1]
+ else
+ return val
+ end
+end
+
+export harmonicMean
+"""
+ harmonicMean(a, b)
+
+Returns the harmonic mean of two numbers `a::Number` and `b::Number`.
+"""
+function harmonicMean(a::Number, b::Number)::Number
+ if a ≈ 0. && b ≈ 0
+ return 0.
+ else
+ return 2.*a*b/(a + b)
+ end
+end
diff --git a/test/contact-search-and-geometry.jl b/test/contact-search-and-geometry.jl
@@ -197,3 +197,47 @@ end
@test 1 == sim.ice_floes[2].n_contacts
@test_throws ErrorException SeaIce.findContacts!(sim, method="")
+
+info("Testing contact registration with multiple contacts")
+sim = SeaIce.createSimulation(id="test")
+SeaIce.addIceFloeCylindrical!(sim, [2., 2.], 1.01, 1., verbose=false)
+SeaIce.addIceFloeCylindrical!(sim, [4., 2.], 1.01, 1., verbose=false)
+SeaIce.addIceFloeCylindrical!(sim, [6., 2.], 1.01, 1., verbose=false)
+SeaIce.addIceFloeCylindrical!(sim, [2., 4.], 1.01, 1., verbose=false)
+SeaIce.addIceFloeCylindrical!(sim, [4., 4.], 1.01, 1., verbose=false)
+SeaIce.addIceFloeCylindrical!(sim, [6., 4.], 1.01, 1., verbose=false)
+SeaIce.addIceFloeCylindrical!(sim, [2., 6.], 1.01, 1., verbose=false)
+SeaIce.addIceFloeCylindrical!(sim, [4., 6.], 1.01, 1., verbose=false)
+SeaIce.addIceFloeCylindrical!(sim, [6., 6.], 1.01, 1., verbose=false)
+sim.ocean = SeaIce.createRegularOceanGrid([4, 4, 2], [8., 8., 2.])
+SeaIce.sortIceFloesInGrid!(sim, sim.ocean)
+SeaIce.findContacts!(sim)
+@test 2 == sim.ice_floes[1].n_contacts
+@test 3 == sim.ice_floes[2].n_contacts
+@test 2 == sim.ice_floes[3].n_contacts
+@test 3 == sim.ice_floes[4].n_contacts
+@test 4 == sim.ice_floes[5].n_contacts
+@test 3 == sim.ice_floes[6].n_contacts
+@test 2 == sim.ice_floes[7].n_contacts
+@test 3 == sim.ice_floes[8].n_contacts
+@test 2 == sim.ice_floes[9].n_contacts
+SeaIce.interact!(sim)
+SeaIce.interact!(sim)
+SeaIce.interact!(sim)
+SeaIce.interact!(sim)
+@test 2 == sim.ice_floes[1].n_contacts
+@test 3 == sim.ice_floes[2].n_contacts
+@test 2 == sim.ice_floes[3].n_contacts
+@test 3 == sim.ice_floes[4].n_contacts
+@test 4 == sim.ice_floes[5].n_contacts
+@test 3 == sim.ice_floes[6].n_contacts
+@test 2 == sim.ice_floes[7].n_contacts
+@test 3 == sim.ice_floes[8].n_contacts
+@test 2 == sim.ice_floes[9].n_contacts
+for i=1:9
+ sim.ice_floes[i].contact_radius = 0.99
+end
+SeaIce.interact!(sim)
+for i=1:9
+ @test sim.ice_floes[i].n_contacts == 0
+end
diff --git a/test/grid.jl b/test/grid.jl
@@ -22,22 +22,23 @@ info("Testing area-determination methods")
@test SeaIce.areaOfQuadrilateral([1., 0.], [0., 1.], [0., 0.], [1., 1.]) ≈ 1.
info("Testing area-based cell content determination")
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.5, 53.5], sw, se, ne, nw) == true
@test SeaIce.isPointInCell(ocean, 1, 1, [6.5, 53.5]) == true
@test SeaIce.getNonDimensionalCellCoordinates(ocean, 1, 1, [6.5, 53.5]) ≈
[.5, .5]
-@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.5]) == true
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.5], sw, se, ne, nw) == true
@test SeaIce.getNonDimensionalCellCoordinates(ocean, 1, 1, [6.1, 53.5]) ≈
[.1, .5]
-@test SeaIce.isPointInCell(ocean, 1, 1, [6.0, 53.5]) == true
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.0, 53.5], sw, se, ne, nw) == true
@test SeaIce.getNonDimensionalCellCoordinates(ocean, 1, 1, [6.0, 53.5]) ≈
[.0, .5]
-@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.7]) == true
-@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.9]) == true
-@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.99999]) == true
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.7], sw, se, ne, nw) == true
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.9], sw, se, ne, nw) == true
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.99999], sw, se, ne, nw) == true
@test SeaIce.getNonDimensionalCellCoordinates(ocean, 1, 1, [6.1, 53.99999]) ≈
[.1, .99999]
-@test SeaIce.isPointInCell(ocean, 1, 1, [7.5, 53.5]) == false
-@test SeaIce.isPointInCell(ocean, 1, 1, [0.0, 53.5]) == false
+@test SeaIce.isPointInCell(ocean, 1, 1, [7.5, 53.5], sw, se, ne, nw) == false
+@test SeaIce.isPointInCell(ocean, 1, 1, [0.0, 53.5], sw, se, ne, nw) == false
x_tilde, _ = SeaIce.getNonDimensionalCellCoordinates(ocean, 1, 1, [0., 53.5])
@test x_tilde < 0.
@@ -64,16 +65,21 @@ info("Testing conformal mapping methods")
[7.5,-1.5])
info("Checking cell content using conformal mapping methods")
-@test SeaIce.isPointInCell(ocean, 1, 1, [6.4, 53.4], method="Conformal") == true
-@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.5], method="Conformal") == true
-@test SeaIce.isPointInCell(ocean, 1, 1, [6.0, 53.5], method="Conformal") == true
-@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.7], method="Conformal") == true
-@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.9], method="Conformal") == true
-@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.99999],
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.4, 53.4], sw, se, ne, nw,
method="Conformal") == true
-@test SeaIce.isPointInCell(ocean, 1, 1, [7.5, 53.5],
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.5], sw, se, ne, nw,
+ method="Conformal") == true
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.0, 53.5], sw, se, ne, nw,
+ method="Conformal") == true
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.7], sw, se, ne, nw,
+ method="Conformal") == true
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.9], sw, se, ne, nw,
+ method="Conformal") == true
+@test SeaIce.isPointInCell(ocean, 1, 1, [6.1, 53.99999], sw, se, ne, nw,
+ method="Conformal") == true
+@test SeaIce.isPointInCell(ocean, 1, 1, [7.5, 53.5], sw, se, ne, nw,
method="Conformal") == false
-@test SeaIce.isPointInCell(ocean, 1, 1, [0.0, 53.5],
+@test SeaIce.isPointInCell(ocean, 1, 1, [0.0, 53.5], sw, se, ne, nw,
method="Conformal") == false
info("Testing bilinear interpolation scheme on conformal mapping")
@@ -81,11 +87,29 @@ ocean.u[1, 1, 1, 1] = 1.0
ocean.u[2, 1, 1, 1] = 1.0
ocean.u[2, 2, 1, 1] = 0.0
ocean.u[1, 2, 1, 1] = 0.0
-@test SeaIce.bilinearInterpolation(ocean.u, .5, .5, 1, 1, 1, 1) ≈ .5
-@test SeaIce.bilinearInterpolation(ocean.u, 1., 1., 1, 1, 1, 1) ≈ .0
-@test SeaIce.bilinearInterpolation(ocean.u, 0., 0., 1, 1, 1, 1) ≈ 1.
-@test SeaIce.bilinearInterpolation(ocean.u, .25, .25, 1, 1, 1, 1) ≈ .75
-@test SeaIce.bilinearInterpolation(ocean.u, .75, .75, 1, 1, 1, 1) ≈ .25
+val = [NaN, NaN]
+SeaIce.bilinearInterpolation!(val, ocean.u[:,:,1,1], ocean.u[:,:,1,1],
+ .5, .5, 1, 1)
+@time SeaIce.bilinearInterpolation!(val, ocean.u[:,:,1,1], ocean.u[:,:,1,1], .5,
+ .5, 1, 1)
+@test val[1] ≈ .5
+@test val[2] ≈ .5
+SeaIce.bilinearInterpolation!(val, ocean.u[:,:,1,1], ocean.u[:,:,1,1], 1., 1.,
+1, 1)
+@test val[1] ≈ .0
+@test val[2] ≈ .0
+SeaIce.bilinearInterpolation!(val, ocean.u[:,:,1,1], ocean.u[:,:,1,1], 0., 0.,
+1, 1)
+@test val[1] ≈ 1.
+@test val[2] ≈ 1.
+SeaIce.bilinearInterpolation!(val, ocean.u[:,:,1,1], ocean.u[:,:,1,1], .25, .25,
+1, 1)
+@test val[1] ≈ .75
+@test val[2] ≈ .75
+SeaIce.bilinearInterpolation!(val, ocean.u[:,:,1,1], ocean.u[:,:,1,1], .75, .75,
+1, 1)
+@test val[1] ≈ .25
+@test val[2] ≈ .25
info("Testing cell binning - Area-based approach")
@test SeaIce.findCellContainingPoint(ocean, [6.2,53.4], method="Area") == (1, 1)
@@ -144,14 +168,18 @@ ocean.u[2, 1, 1, 1] = 1.0
ocean.u[2, 2, 1, 1] = 0.0
ocean.u[1, 2, 1, 1] = 0.0
ocean.v[:, :, 1, 1] = 0.0
-@test SeaIce.curl(ocean, .5, .5, 1, 1, 1, 1) > 0.
+sw = Vector{Float64}(2)
+se = Vector{Float64}(2)
+ne = Vector{Float64}(2)
+nw = Vector{Float64}(2)
+@test SeaIce.curl(ocean, .5, .5, 1, 1, 1, 1, sw, se, ne, nw) > 0.
ocean.u[1, 1, 1, 1] = 0.0
ocean.u[2, 1, 1, 1] = 0.0
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.
+@test SeaIce.curl(ocean, .5, .5, 1, 1, 1, 1, sw, se, ne, nw) < 0.
info("Testing atmosphere drag")
sim = SeaIce.createSimulation()
@@ -186,6 +214,7 @@ 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.
+@test SeaIce.curl(atmosphere, .5, .5, 1, 1, 1, 1, sw, se, ne, nw) > 0.
atmosphere.u[1, 1, 1, 1] = 0.0
atmosphere.u[2, 1, 1, 1] = 0.0
@@ -193,6 +222,7 @@ 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.
+@test SeaIce.curl(atmosphere, .5, .5, 1, 1, 1, 1, sw, se, ne, nw) < 0.
info("Testing findEmptyPositionInGridCell")
diff --git a/test/icefloe.jl b/test/icefloe.jl
@@ -28,3 +28,13 @@ sim = SeaIce.createSimulation(id="test")
SeaIce.addIceFloeCylindrical!(sim, [ 0., 0.], 10., 1., verbose=false)
SeaIce.addIceFloeCylindrical!(sim, [ 0., 0.], 10., 1., verbose=false)
SeaIce.compareIceFloes(sim.ice_floes[1], sim.ice_floes[2])
+
+SeaIce.plotIceFloeSizeDistribution(sim)
+rm("test-ice-floe-size-distribution.png")
+SeaIce.plotIceFloeSizeDistribution(sim, skip_fixed=false)
+rm("test-ice-floe-size-distribution.png")
+SeaIce.plotIceFloeSizeDistribution(sim, log_y=false)
+rm("test-ice-floe-size-distribution.png")
+SeaIce.plotIceFloeSizeDistribution(sim, size_type="areal")
+rm("test-ice-floe-size-distribution.png")
+@test_throws ErrorException SeaIce.plotIceFloeSizeDistribution(sim, size_type="asdf")
diff --git a/test/memory-management.jl b/test/memory-management.jl
@@ -0,0 +1,175 @@
+#!/usr/bin/env julia
+import SeaIce
+using Base.Test
+
+info("#### $(basename(@__FILE__)) ####")
+
+info("Testing memory footprint of SeaIce types")
+
+sim = SeaIce.createSimulation()
+empty_sim_size = 96
+empty_sim_size_recursive = 472
+
+@test sizeof(sim) == empty_sim_size
+@test Base.summarysize(sim) == empty_sim_size_recursive
+
+size_per_icefloe = 352
+size_per_icefloe_recursive = 1136
+
+info("Testing memory usage when adding ice floes")
+for i=1:100
+ SeaIce.addIceFloeCylindrical!(sim, [1., 1.], 1., 1., verbose=false)
+
+ @test sizeof(sim) == empty_sim_size
+
+ @test sizeof(sim.ice_floes) == sizeof(Int)*i
+ @test sizeof(sim.ice_floes[:]) == sizeof(Int)*i
+ @test Base.summarysize(sim.ice_floes) == size_per_icefloe_recursive*i +
+ sizeof(Int)*i
+
+ @test Base.summarysize(sim) == empty_sim_size_recursive + sizeof(Int)*i +
+ size_per_icefloe_recursive*i
+
+ @test Base.summarysize(sim.ice_floes[i]) == size_per_icefloe_recursive
+
+ for j=1:i
+ @test sizeof(sim.ice_floes[j]) == size_per_icefloe
+ @test Base.summarysize(sim.ice_floes[j]) == size_per_icefloe_recursive
+ end
+
+end
+
+info("Checking memory footprint when overwriting simulation object")
+sim = SeaIce.createSimulation()
+@test sizeof(sim) == empty_sim_size
+@test Base.summarysize(sim) == empty_sim_size_recursive
+
+info("Check memory usage when stepping time for empty simulation object")
+sim = SeaIce.createSimulation()
+sim.time_step = 1.0
+for i=1:10
+ SeaIce.run!(sim, single_step=true, verbose=false)
+ @test sizeof(sim) == empty_sim_size
+ @test Base.summarysize(sim) == empty_sim_size_recursive
+end
+
+info("Check memory when stepping time with single ice floe")
+sim = SeaIce.createSimulation()
+SeaIce.addIceFloeCylindrical!(sim, [1., 1.], 1., 1., verbose=false)
+sim.time_step = 1.0
+for i=1:10
+ SeaIce.run!(sim, single_step=true, verbose=false)
+ @test sizeof(sim) == empty_sim_size
+ @test Base.summarysize(sim) == empty_sim_size_recursive +
+ sizeof(Int)*length(sim.ice_floes) +
+ size_per_icefloe_recursive*length(sim.ice_floes)
+end
+
+info("Check memory when stepping time with two separate ice floes")
+sim = SeaIce.createSimulation()
+SeaIce.addIceFloeCylindrical!(sim, [1., 1.], 1., 1., verbose=false)
+SeaIce.addIceFloeCylindrical!(sim, [1., 1.], 3., 1., verbose=false)
+sim.time_step = 1.0
+for i=1:10
+ SeaIce.run!(sim, single_step=true, verbose=false)
+ @test sizeof(sim) == empty_sim_size
+ @test Base.summarysize(sim) == empty_sim_size_recursive +
+ sizeof(Int)*length(sim.ice_floes) +
+ size_per_icefloe_recursive*length(sim.ice_floes)
+end
+
+info("Check memory when stepping time with two interacting ice floes (all to all)")
+sim = SeaIce.createSimulation()
+SeaIce.addIceFloeCylindrical!(sim, [1., 1.], 1., 1., verbose=false)
+SeaIce.addIceFloeCylindrical!(sim, [1., 1.], 1.9, 1., verbose=false)
+sim.time_step = 1.0
+for i=1:10
+ SeaIce.run!(sim, single_step=true, verbose=false)
+ @test sizeof(sim) == empty_sim_size
+ @test Base.summarysize(sim) == empty_sim_size_recursive +
+ sizeof(Int)*length(sim.ice_floes) +
+ size_per_icefloe_recursive*length(sim.ice_floes)
+end
+
+info("Check memory when stepping time with two interacting ice floes (cell sorting)")
+sim = SeaIce.createSimulation()
+SeaIce.addIceFloeCylindrical!(sim, [1., 1.], 1., 1., verbose=false)
+SeaIce.addIceFloeCylindrical!(sim, [1., 1.], 1.9, 1., verbose=false)
+nx, ny, nz = 5, 5, 1
+sim.ocean = SeaIce.createRegularOceanGrid([nx, ny, nz], [10., 10., 10.])
+sim.time_step = 1e-6
+SeaIce.run!(sim, single_step=true, verbose=false)
+original_size_recursive = Base.summarysize(sim)
+for i=1:10
+ SeaIce.run!(sim, single_step=true, verbose=false)
+ @test Base.summarysize(sim) == original_size_recursive
+end
+
+info("Checking if memory is freed after ended collision (all to all)")
+sim = SeaIce.createSimulation(id="test")
+SeaIce.addIceFloeCylindrical!(sim, [0., 0.], 10., 1., verbose=false)
+SeaIce.addIceFloeCylindrical!(sim, [20.05, 0.], 10., 1., verbose=false)
+sim.ice_floes[1].lin_vel[1] = 0.1
+SeaIce.setTotalTime!(sim, 10.0)
+SeaIce.setTimeStep!(sim, epsilon=0.07, verbose=false)
+original_size_recursive = Base.summarysize(sim)
+SeaIce.run!(sim, verbose=false)
+@test Base.summarysize(sim) == original_size_recursive
+
+info("Checking if memory is freed after ended collision (cell sorting)")
+sim = SeaIce.createSimulation(id="test")
+SeaIce.addIceFloeCylindrical!(sim, [0., 0.], 10., 1., verbose=false)
+SeaIce.addIceFloeCylindrical!(sim, [20.05, 0.], 10., 1., verbose=false)
+sim.ocean = SeaIce.createRegularOceanGrid([2, 2, 2], [40., 40., 10.])
+sim.ice_floes[1].lin_vel[1] = 0.1
+SeaIce.setTotalTime!(sim, 10.0)
+SeaIce.setTimeStep!(sim, epsilon=0.07, verbose=false)
+SeaIce.run!(sim, single_step=true, verbose=false)
+original_sim_size_recursive = Base.summarysize(sim)
+original_icefloes_size_recursive = Base.summarysize(sim.ice_floes)
+original_ocean_size_recursive = Base.summarysize(sim.ocean)
+original_atmosphere_size_recursive = Base.summarysize(sim.atmosphere)
+SeaIce.run!(sim, verbose=false)
+@test Base.summarysize(sim.ice_floes) == original_icefloes_size_recursive
+@test Base.summarysize(sim.ocean) == original_ocean_size_recursive
+@test Base.summarysize(sim.atmosphere) == original_atmosphere_size_recursive
+@test Base.summarysize(sim) == original_sim_size_recursive
+
+sim = SeaIce.createSimulation(id="test")
+SeaIce.addIceFloeCylindrical!(sim, [0., 0.], 10., 1., verbose=false)
+SeaIce.addIceFloeCylindrical!(sim, [20.05, 0.], 10., 1., verbose=false)
+sim.atmosphere = SeaIce.createRegularAtmosphereGrid([2, 2, 2], [40., 40., 10.])
+sim.ice_floes[1].lin_vel[1] = 0.1
+SeaIce.setTotalTime!(sim, 10.0)
+SeaIce.setTimeStep!(sim, epsilon=0.07, verbose=false)
+SeaIce.run!(sim, single_step=true, verbose=false)
+original_sim_size_recursive = Base.summarysize(sim)
+original_icefloes_size_recursive = Base.summarysize(sim.ice_floes)
+original_ocean_size_recursive = Base.summarysize(sim.ocean)
+original_atmosphere_size_recursive = Base.summarysize(sim.atmosphere)
+SeaIce.run!(sim, verbose=false)
+@test Base.summarysize(sim.ice_floes) == original_icefloes_size_recursive
+@test Base.summarysize(sim.ocean) == original_ocean_size_recursive
+@test Base.summarysize(sim.atmosphere) == original_atmosphere_size_recursive
+@test Base.summarysize(sim) == original_sim_size_recursive
+
+sim = SeaIce.createSimulation(id="test")
+SeaIce.addIceFloeCylindrical!(sim, [0., 0.], 10., 1., verbose=false)
+SeaIce.addIceFloeCylindrical!(sim, [20.05, 0.], 10., 1., verbose=false)
+sim.atmosphere = SeaIce.createRegularAtmosphereGrid([2, 2, 2], [40., 40., 10.])
+sim.ocean = SeaIce.createRegularOceanGrid([2, 2, 2], [40., 40., 10.])
+sim.ice_floes[1].lin_vel[1] = 0.1
+SeaIce.setTotalTime!(sim, 10.0)
+SeaIce.setTimeStep!(sim, epsilon=0.07, verbose=false)
+SeaIce.run!(sim, single_step=true, verbose=false)
+original_sim_size_recursive = Base.summarysize(sim)
+original_icefloes_size_recursive = Base.summarysize(sim.ice_floes)
+original_ocean_size_recursive = Base.summarysize(sim.ocean)
+original_atmosphere_size_recursive = Base.summarysize(sim.atmosphere)
+SeaIce.run!(sim, verbose=false)
+@test Base.summarysize(sim.ice_floes) == original_icefloes_size_recursive
+@test Base.summarysize(sim.ocean) == original_ocean_size_recursive
+@test Base.summarysize(sim.atmosphere) == original_atmosphere_size_recursive
+@test Base.summarysize(sim) == original_sim_size_recursive
+
+SeaIce.printMemoryUsage(sim)
diff --git a/test/profiling.jl b/test/profiling.jl
@@ -85,7 +85,7 @@ function timeSingleStepInDenseSimulation(nx::Int; verbose::Bool=true,
@test sim.ice_floes[nx].n_contacts == 0
@test sim.ice_floes[nx + 1].n_contacts == 1
@test sim.ice_floes[nx + 2].n_contacts == 4
- return t_elapsed
+ return t_elapsed, Base.summarysize(sim)
end
#nx = Int[4 8 16 32 64 96]
@@ -95,13 +95,16 @@ t_elapsed = zeros(length(nx))
t_elapsed_all_to_all = zeros(length(nx))
t_elapsed_cell_sorting = zeros(length(nx))
t_elapsed_cell_sorting2 = zeros(length(nx))
+memory_usage_all_to_all = zeros(length(nx))
+memory_usage_cell_sorting = zeros(length(nx))
+memory_usage_cell_sorting2 = zeros(length(nx))
for i=1:length(nx)
info("nx = $(nx[i])")
- t_elapsed_all_to_all[i] =
+ t_elapsed_all_to_all[i], memory_usage_all_to_all[i] =
timeSingleStepInDenseSimulation(Int(nx[i]), grid_sorting=false)
- t_elapsed_cell_sorting[i] =
+ t_elapsed_cell_sorting[i], memory_usage_cell_sorting[i] =
timeSingleStepInDenseSimulation(Int(nx[i]), grid_sorting=true)
- t_elapsed_cell_sorting2[i] =
+ t_elapsed_cell_sorting2[i], memory_usage_cell_sorting2[i] =
timeSingleStepInDenseSimulation(Int(nx[i]), grid_sorting=true,
include_atmosphere=true)
elements[i] = nx[i]*nx[i]
@@ -148,4 +151,46 @@ Plots.plot!(elements, fit_cell_sorting2(elements),
Plots.title!("Dense granular system " * "(host: $(gethostname()))")
Plots.xaxis!("Number of ice floes")
Plots.yaxis!("Wall time per time step [s]")
-Plots.savefig("profiling.pdf")
+Plots.savefig("profiling-cpu.pdf")
+
+Plots.scatter(elements, memory_usage_all_to_all .÷ 1024,
+ xscale=:log10,
+ yscale=:log10,
+ label="All to all")
+fit_all_to_all = CurveFit.curve_fit(CurveFit.PowerFit,
+ elements, memory_usage_all_to_all .÷ 1024)
+label_all_to_all = @sprintf "%1.3g n^%3.2f" fit_all_to_all.coefs[1] fit_all_to_all.coefs[2]
+Plots.plot!(elements, fit_all_to_all(elements),
+ xscale=:log10,
+ yscale=:log10,
+ label=label_all_to_all)
+
+Plots.scatter!(elements, memory_usage_cell_sorting .÷ 1024,
+ xscale=:log10,
+ yscale=:log10,
+ label="Cell-based spatial decomposition (ocean only)")
+fit_cell_sorting = CurveFit.curve_fit(CurveFit.PowerFit,
+ elements, memory_usage_cell_sorting .÷ 1024)
+label_cell_sorting = @sprintf "%1.3g n^%3.2f" fit_cell_sorting.coefs[1] fit_cell_sorting.coefs[2]
+Plots.plot!(elements, fit_cell_sorting(elements),
+ xscale=:log10,
+ yscale=:log10,
+ label=label_cell_sorting)
+
+Plots.scatter!(elements, memory_usage_cell_sorting2 .÷ 1024,
+ xscale=:log10,
+ yscale=:log10,
+ label="Cell-based spatial decomposition (ocean + atmosphere)")
+fit_cell_sorting2 = CurveFit.curve_fit(CurveFit.PowerFit,
+ elements,
+ memory_usage_cell_sorting2 .÷ 1024)
+label_cell_sorting2 = @sprintf "%1.3g n^%3.2f" fit_cell_sorting2.coefs[1] fit_cell_sorting2.coefs[2]
+Plots.plot!(elements, fit_cell_sorting2(elements),
+ xscale=:log10,
+ yscale=:log10,
+ label=label_cell_sorting2)
+
+Plots.title!("Dense granular system " * "(host: $(gethostname()))")
+Plots.xaxis!("Number of ice floes")
+Plots.yaxis!("Memory usage [kb]")
+Plots.savefig("profiling-memory-usage.pdf")
diff --git a/test/profiling.sh b/test/profiling.sh
@@ -0,0 +1,18 @@
+#!/bin/bash
+
+# Optionally use this script to launch profiling.jl with different Julia
+# optimization levels. Defaults are --optimize=2, --math-mode=ieee.
+
+declare -a arr=(\
+ "--procs 1 --optimize=1 --math-mode=ieee" \
+ "--procs 1 --optimize=2 --math-mode=ieee" \
+ "--procs 1 --optimize=3 --math-mode=ieee" \
+ "--procs 1 --optimize=1 --math-mode=fast" \
+ "--procs 1 --optimize=2 --math-mode=fast" \
+ "--procs 1 --optimize=3 --math-mode=fast")
+
+for flags in "${arr[@]}"; do
+ julia --color=yes $flags profiling.jl
+ mv profiling-cpu.pdf "profiling-cpu.$flags.pdf"
+ mv profiling-memory-usage.pdf "profiling-memory-usage.$flags.pdf"
+done
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -2,6 +2,7 @@ import SeaIce
using Base.Test
include("icefloe.jl")
+include("util.jl")
include("temporal.jl")
include("contact-search-and-geometry.jl")
include("collision-2floes-normal.jl")
@@ -14,3 +15,4 @@ include("jld.jl")
include("grid.jl")
include("ocean.jl")
include("atmosphere.jl")
+include("memory-management.jl")
diff --git a/test/util.jl b/test/util.jl
@@ -0,0 +1,28 @@
+#!/usr/bin/env julia
+import SeaIce
+using Base.Test
+
+info("#### $(basename(@__FILE__)) ####")
+
+info("Testing power-law RNG")
+
+@test 1 == length(SeaIce.randpower())
+@test () == size(SeaIce.randpower())
+@test 1 == length(SeaIce.randpower(1))
+@test () == size(SeaIce.randpower(1))
+@test 4 == length(SeaIce.randpower((2,2)))
+@test (2,2) == size(SeaIce.randpower((2,2)))
+@test 5 == length(SeaIce.randpower(5))
+@test (5,) == size(SeaIce.randpower(5))
+
+Base.Random.srand(1)
+for i=1:10^5
+ @test 0. <= SeaIce.randpower() <= 1.
+ @test 0. <= SeaIce.randpower(1, 1., 0., 1.) <= 1.
+ @test 0. <= SeaIce.randpower(1, 1., 0., .1) <= .1
+ @test 5. <= SeaIce.randpower(1, 1., 5., 6.) <= 6.
+ @test 0. <= minimum(SeaIce.randpower((2,2), 1., 0., 1.))
+ @test 1. >= maximum(SeaIce.randpower((2,2), 1., 0., 1.))
+ @test 0. <= minimum(SeaIce.randpower(5, 1., 0., 1.))
+ @test 1. >= minimum(SeaIce.randpower(5, 1., 0., 1.))
+end
