timplement tangential elasticity - 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 09f45f8ba9a400b62df5ac97affca06066b3707c
(DIR) parent 5e16f290d31e8d0d24032fc60ec40a5dfbcf52f6
(HTM) Author: Anders Damsgaard <andersd@riseup.net>
Date: Wed, 10 May 2017 13:26:59 -0400
implement tangential elasticity
Diffstat:
M examples/nares_strait.jl | 13 ++++++-------
M src/contact_search.jl | 33 ++++++++++++++++++++++---------
M src/datatypes.jl | 8 +++++---
M src/icefloe.jl | 26 +++++++++++++++++++++-----
M src/interaction.jl | 84 ++++++++++++++++---------------
M src/simulation.jl | 18 +++---------------
M test/collision-2floes-normal.jl | 92 +++++++++++++++++++++++++++++++
M test/collision-2floes-oblique.jl | 337 ++++++++++++++++++++++---------
M test/contact-search-and-geometry.jl | 134 ++++++++++++++++++++++++++-----
M test/vtk.jl | 2 +-
10 files changed, 549 insertions(+), 198 deletions(-)
---
(DIR) diff --git a/examples/nares_strait.jl b/examples/nares_strait.jl
t@@ -1,11 +1,11 @@
#!/usr/bin/env julia
import SeaIce
-sim = SeaIce.createSimulation(id="nares_strait")
-n = [25, 25, 2]
+#sim = SeaIce.createSimulation(id="nares_strait")
+#n = [25, 25, 2]
-#sim = SeaIce.createSimulation(id="nares_strait_coarse")
-#n = [6, 6, 2]
+sim = SeaIce.createSimulation(id="nares_strait_coarse")
+n = [6, 6, 2]
# Initialize ocean
Lx = 50.e3
t@@ -110,7 +110,7 @@ info("added $(n) ice floes")
SeaIce.removeSimulationFiles(sim)
k_n = 1e6 # N/m
-k_t = k_n
+k_t = 0.
gamma_t = 1e7 # N/(m/s)
mu_d = 0.7
rotating = true
t@@ -130,8 +130,7 @@ SeaIce.setTimeStep!(sim)
# Run simulation for 10 time steps, then add new icefloes from the top
while sim.time < sim.time_total
for it=1:10
- SeaIce.run!(sim, status_interval=1, single_step=true,
- contact_tangential_rheology="Linear Viscous Frictional")
+ SeaIce.run!(sim, status_interval=1, single_step=true)
end
for i=1:size(sim.ocean.xh, 1)
if sim.ocean.ice_floe_list[i, end] == []
(DIR) diff --git a/src/contact_search.jl b/src/contact_search.jl
t@@ -120,23 +120,38 @@ written to `simulation.contact_parallel_displacement`.
* `i::Int`: index of the first ice floe.
* `j::Int`: index of the second ice floe.
"""
-function checkAndAddContact!(simulation::Simulation, i::Int, j::Int)
+function checkAndAddContact!(sim::Simulation, i::Int, j::Int)
if i < j
- if (simulation.ice_floes[i].fixed && simulation.ice_floes[j].fixed) ||
- !simulation.ice_floes[i].enabled ||
- !simulation.ice_floes[j].enabled
+ if (sim.ice_floes[i].fixed && sim.ice_floes[j].fixed) ||
+ !sim.ice_floes[i].enabled ||
+ !sim.ice_floes[j].enabled
return
end
# Inter-grain position vector and grain overlap
- position_ij = interIceFloePositionVector(simulation, i, j)
- overlap_ij = findOverlap(simulation, i, j, position_ij)
+ position_ij = interIceFloePositionVector(sim, i, j)
+ overlap_ij = findOverlap(sim, i, j, position_ij)
# Check if grains overlap (overlap when negative)
- if overlap_ij < 0.0
- push!(simulation.contact_pairs, [i, j])
- push!(simulation.contact_parallel_displacement, zeros(2))
+ if overlap_ij < 0.
+ for ic=1:(Nc_max + 1)
+ if ic == (Nc_max + 1)
+ error("contact $i-$j exceeds max. number of contacts " *
+ "(Nc_max = $Nc_max) for ice floe $i")
+
+ elseif 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
+ sim.ice_floes[i].contact_parallel_displacement[ic] =
+ zeros(2)
+ break
+ end
+ end
end
end
end
(DIR) diff --git a/src/datatypes.jl b/src/datatypes.jl
t@@ -44,9 +44,12 @@ type IceFloeCylindrical
contact_static_friction::float
contact_dynamic_friction::float
+ # Interaction
pressure::float
-
+ n_contacts::Int
ocean_grid_pos::Array{Int, 1}
+ contacts::Array{Int, 1}
+ contact_parallel_displacement::Array{Array{Float64, 1}, 1}
end
# Type for gathering data from ice floe objects into single arrays
t@@ -92,6 +95,7 @@ type IceFloeArrays
contact_dynamic_friction
pressure
+ n_contacts
end
#=
t@@ -171,8 +175,6 @@ type Simulation
file_time_since_output_file::Float64
ice_floes::Array{IceFloeCylindrical, 1}
- contact_pairs::Array{Array{Int, 1}, 1}
- contact_parallel_displacement::Array{Array{Float64, 1}, 1}
ocean::Ocean
end
(DIR) diff --git a/src/icefloe.jl b/src/icefloe.jl
t@@ -1,5 +1,7 @@
## Manage icefloes in the model
+Nc_max = 32 # max. no. of contacts per ice floe
+
export addIceFloeCylindrical
"""
Adds a grain to the simulation. Example:
t@@ -20,7 +22,7 @@ function addIceFloeCylindrical(simulation::Simulation,
torque::float = 0.,
density::float = 934.,
contact_stiffness_normal::float = 1.e6,
- contact_stiffness_tangential::float = 1.e6,
+ contact_stiffness_tangential::float = 0.,
contact_viscosity_normal::float = 0.,
contact_viscosity_tangential::float = 0.,
contact_static_friction::float = 0.4,
t@@ -30,7 +32,13 @@ function addIceFloeCylindrical(simulation::Simulation,
rotating::Bool = true,
enabled::Bool = true,
verbose::Bool = true,
- ocean_grid_pos::Array{Int, 1} = [0, 0])
+ ocean_grid_pos::Array{Int, 1} = [0, 0],
+ n_contacts::Int = 0,
+ contacts::Array{Int, 1} = zeros(Int, Nc_max),
+ contact_parallel_displacement::
+ Array{Array{Float64, 1}, 1}
+ =
+ Array{Array{Float64, 1}, 1}(Nc_max))
# Check input values
if length(lin_pos) != 2
t@@ -65,6 +73,9 @@ function addIceFloeCylindrical(simulation::Simulation,
areal_radius = contact_radius
end
+ for i=1:Nc_max
+ contact_parallel_displacement[i] = zeros(2)
+ end
# Create icefloe object with placeholder values for surface area, volume,
# mass, and moment of inertia.
t@@ -101,8 +112,10 @@ function addIceFloeCylindrical(simulation::Simulation,
contact_dynamic_friction,
pressure,
-
- ocean_grid_pos
+ n_contacts,
+ ocean_grid_pos,
+ contacts,
+ contact_parallel_displacement
)
# Overwrite previous placeholder values
t@@ -194,7 +207,8 @@ function convertIceFloeDataToArrays(simulation::Simulation)
Array(Float64, length(simulation.ice_floes)),
Array(Float64, length(simulation.ice_floes)),
- Array(Float64, length(simulation.ice_floes))
+ Array(Float64, length(simulation.ice_floes)),
+ Array(Int, length(simulation.ice_floes))
)
# fill arrays
t@@ -240,6 +254,8 @@ function convertIceFloeDataToArrays(simulation::Simulation)
simulation.ice_floes[i].contact_dynamic_friction
ifarr.pressure[i] = simulation.ice_floes[i].pressure
+
+ ifarr.n_contacts[i] = simulation.ice_floes[i].n_contacts
end
return ifarr
(DIR) diff --git a/src/interaction.jl b/src/interaction.jl
t@@ -2,22 +2,28 @@
export interact!
"""
+ interact!(simulation::Simulation)
+
Resolve mechanical interaction between all particle pairs.
"""
-function interact!(simulation::Simulation;
- contact_normal_rheology::String = "Linear Elastic",
- contact_tangential_rheology::String = "None")
-
- # IceFloe to grain collisions
- while !isempty(simulation.contact_pairs)
- contact_pair = pop!(simulation.contact_pairs)
- contact_parallel_displacement =
- pop!(simulation.contact_parallel_displacement)
- interactIceFloes!(simulation, contact_pair[1], contact_pair[2],
- contact_parallel_displacement,
- contact_normal_rheology=contact_normal_rheology,
- contact_tangential_rheology=
- contact_tangential_rheology)
+function interact!(simulation::Simulation)
+ for i=1:Int(ceil(length(simulation.ice_floes)/2.)) # i <= Int(N/2)
+ for ic=1:simulation.ice_floes[i].n_contacts
+ #for ic=1:Nc_max
+
+ j = simulation.ice_floes[i].contacts[ic]
+
+ 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].contacts[ic] = 0 # remove contact
+ simulation.ice_floes[i].n_contacts -= 1
+ else
+ interactIceFloes!(simulation, i, j, ic)
+ end
+ end
end
end
t@@ -27,11 +33,7 @@ Resolve an grain-to-grain interaction using a prescibed contact law. This
function adds the compressive force of the interaction to the ice floe
`pressure` field of mean compressive stress on the ice floe sides.
"""
-function interactIceFloes!(simulation::Simulation,
- i::Int, j::Int,
- contact_parallel_displacement::vector;
- contact_normal_rheology::String = "Linear Elastic",
- contact_tangential_rheology::String = "None")
+function interactIceFloes!(simulation::Simulation, i::Int, j::Int, ic::Int)
force_n = 0. # Contact-normal force
force_t = 0. # Contact-parallel (tangential) force
t@@ -61,8 +63,8 @@ function interactIceFloes!(simulation::Simulation,
simulation.ice_floes[j].ang_vel)
# Correct old tangential displacement for contact rotation and add new
- delta_t = dot(t, contact_parallel_displacement -
- (n*dot(n, contact_parallel_displacement))) +
+ delta_t0 =simulation.ice_floes[i].contact_parallel_displacement[ic]
+ delta_t = dot(t, delta_t0 - (n*dot(n, delta_t0))) +
vel_t*simulation.time_step
# Effective radius
t@@ -91,35 +93,27 @@ function interactIceFloes!(simulation::Simulation,
simulation.ice_floes[j].contact_dynamic_friction)
# Determine contact forces
- if contact_normal_rheology == "None"
+ if k_n_harmonic_mean ≈ 0. && gamma_n_harmonic_mean ≈ 0.
force_n = 0.
- elseif contact_normal_rheology == "Linear Elastic"
+ elseif k_n_harmonic_mean > 0. && gamma_n_harmonic_mean ≈ 0.
force_n = -k_n_harmonic_mean*delta_n
- elseif contact_normal_rheology == "Linear Viscous Elastic"
+ elseif k_n_harmonic_mean > 0. && gamma_n_harmonic_mean > 0.
force_n = -k_n_harmonic_mean*delta_n - gamma_n_harmonic_mean*vel_n
if force_n < 0.
force_n = 0.
end
else
- error("unknown contact_normal_rheology '$contact_normal_rheology'")
+ error("unknown contact_normal_rheology (k_n = $k_n_harmonic_mean," *
+ " gamma_n = $gamma_n_harmonic_mean")
end
- if contact_tangential_rheology == "None"
+ if k_t_harmonic_mean ≈ 0. && gamma_t_harmonic_mean ≈ 0.
# do nothing
- elseif contact_tangential_rheology ==
- "Linear Elastic Viscous Frictional"
- force_t = -k_t_harmonic_mean*delta_t - gamma_t_harmonic_mean*vel_t
- if abs(force_t) > mu_d_minimum*abs(force_n)
- force_t = mu_d_minimum*abs(force_n)*force_t/abs(force_t)
- delta_t = (-force_t - gamma_t_harmonic_mean*vel_t)/
- k_t_harmonic_mean
- end
-
- elseif contact_tangential_rheology == "Linear Viscous Frictional"
+ elseif k_t_harmonic_mean ≈ 0. && gamma_t_harmonic_mean > 0.
force_t = abs(gamma_t_harmonic_mean * vel_t)
if force_t > mu_d_minimum*abs(force_n)
force_t = mu_d_minimum*abs(force_n)
t@@ -128,17 +122,27 @@ function interactIceFloes!(simulation::Simulation,
force_t = -force_t
end
+ elseif k_t_harmonic_mean > 0.
+
+ force_t = -k_t_harmonic_mean*delta_t -
+ gamma_t_harmonic_mean*vel_t
+
+ if abs(force_t) > mu_d_minimum*abs(force_n)
+ force_t = mu_d_minimum*abs(force_n)*force_t/abs(force_t)
+ delta_t = (-force_t - gamma_t_harmonic_mean*vel_t)/
+ k_t_harmonic_mean
+ end
+
else
- error("unknown contact_tangential_rheology ",
- "'$contact_tangential_rheology'")
+ error("unknown contact_tangential_rheology (k_t = " *
+ "$k_t_harmonic_mean, gamma_t = $gamma_t_harmonic_mean")
end
else
error("function called to process non-existent contact between ice " *
"floes $i and $j")
end
-
- contact_parallel_displacement = delta_t*t
+ simulation.ice_floes[i].contact_parallel_displacement[ic] = delta_t*t
simulation.ice_floes[i].force += force_n*n + force_t*t;
simulation.ice_floes[j].force -= force_n*n + force_t*t;
(DIR) diff --git a/src/simulation.jl b/src/simulation.jl
t@@ -9,10 +9,7 @@ export createSimulation
time_step::Float64=-1.,
file_time_step::Float64=-1.,
file_number::Int=0,
- ice_floes=Array{IceFloeCylindrical, 1}[],
- contact_pairs=Array{Int64, 1}[],
- contact_parallel_displacement=
- Array{Array{Float64, 1}, 1}[])
+ ice_floes=Array{IceFloeCylindrical, 1}[])
Create a simulation object containing all relevant variables such as temporal
parameters, and lists of ice floes and contacts.
t@@ -29,9 +26,6 @@ function createSimulation(;id::String="unnamed",
file_number::Int=0,
file_time_since_output_file::Float64=0.,
ice_floes=Array{IceFloeCylindrical, 1}[],
- contact_pairs=Array{Int64, 1}[],
- contact_parallel_displacement=
- Array{Array{Float64, 1}, 1}[],
ocean::Ocean=createEmptyOcean())
return Simulation(id,
t@@ -43,8 +37,6 @@ function createSimulation(;id::String="unnamed",
file_number,
file_time_since_output_file,
ice_floes,
- contact_pairs,
- contact_parallel_displacement,
ocean)
end
t@@ -85,9 +77,7 @@ function run!(simulation::Simulation;
status_interval::Int=100,
show_file_output::Bool=true,
single_step::Bool=false,
- temporal_integration_method::String="Three-term Taylor",
- contact_normal_rheology::String = "Linear Elastic",
- contact_tangential_rheology::String = "None")
+ temporal_integration_method::String="Three-term Taylor")
if single_step && simulation.time >= simulation.time_total
simulation.time_total += simulation.time_step
t@@ -119,9 +109,7 @@ function run!(simulation::Simulation;
else
findContacts!(simulation, method="all to all")
end
- interact!(simulation,
- contact_normal_rheology=contact_normal_rheology,
- contact_tangential_rheology=contact_tangential_rheology)
+ interact!(simulation)
if typeof(simulation.ocean.input_file) != Bool
addOceanDrag!(simulation)
end
(DIR) diff --git a/test/collision-2floes-normal.jl b/test/collision-2floes-normal.jl
t@@ -115,3 +115,95 @@ E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
@test_approx_eq_eps E_kin_lin_init E_kin_lin_final E_kin_lin_init*tol
@test_approx_eq E_kin_rot_init E_kin_rot_final
+
+
+info("# Adding contact-normal viscosity")
+info("# One ice floe fixed")
+sim = SeaIce.createSimulation(id="test")
+SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., verbose=verbose)
+SeaIce.addIceFloeCylindrical(sim, [20.05, 0.], 10., 1., verbose=verbose)
+sim.ice_floes[1].lin_vel[1] = 0.1
+sim.ice_floes[1].contact_viscosity_normal = 1e4
+sim.ice_floes[2].contact_viscosity_normal = 1e4
+sim.ice_floes[2].fixed = true
+
+E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+
+# With decreasing timestep (epsilon towards 0), the explicit integration scheme
+# should become more correct
+
+SeaIce.setTotalTime!(sim, 10.0)
+sim_init = deepcopy(sim)
+
+
+info("Testing kinetic energy conservation with Two-term Taylor scheme")
+sim = deepcopy(sim_init)
+SeaIce.setTimeStep!(sim, epsilon=0.007)
+tol = 0.02
+info("Relative tolerance: $(tol*100.)%")
+SeaIce.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbose)
+
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test E_kin_lin_init > E_kin_lin_final
+@test_approx_eq E_kin_rot_init E_kin_rot_final
+
+
+info("Testing kinetic energy conservation with Three-term Taylor scheme")
+sim = deepcopy(sim_init)
+SeaIce.setTimeStep!(sim, epsilon=0.07)
+tol = 0.01
+info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
+SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
+ verbose=verbose)
+
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test E_kin_lin_init > E_kin_lin_final
+@test_approx_eq E_kin_rot_init E_kin_rot_final
+
+
+info("# Ice floes free to move")
+
+sim = SeaIce.createSimulation(id="test")
+SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., verbose=verbose)
+SeaIce.addIceFloeCylindrical(sim, [20.05, 0.], 10., 1., verbose=verbose)
+sim.ice_floes[1].lin_vel[1] = 0.1
+sim.ice_floes[1].contact_viscosity_normal = 1e4
+sim.ice_floes[2].contact_viscosity_normal = 1e4
+
+E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+
+# With decreasing timestep (epsilon towards 0), the explicit integration scheme
+# should become more correct
+
+SeaIce.setTotalTime!(sim, 10.0)
+sim_init = deepcopy(sim)
+
+info("Testing kinetic energy conservation with Two-term Taylor scheme")
+sim = deepcopy(sim_init)
+SeaIce.setTimeStep!(sim, epsilon=0.007)
+tol = 0.02
+info("Relative tolerance: $(tol*100.)%")
+SeaIce.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbose)
+
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test E_kin_lin_init > E_kin_lin_final
+@test_approx_eq E_kin_rot_init E_kin_rot_final
+
+
+info("Testing kinetic energy conservation with Three-term Taylor scheme")
+sim = deepcopy(sim_init)
+SeaIce.setTimeStep!(sim, epsilon=0.07)
+tol = 0.01
+info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
+SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
+ verbose=verbose)
+
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test E_kin_lin_init > E_kin_lin_final
+@test_approx_eq E_kin_rot_init E_kin_rot_final
(DIR) diff --git a/test/collision-2floes-oblique.jl b/test/collision-2floes-oblique.jl
t@@ -29,8 +29,7 @@ info("Testing kinetic energy conservation with Two-term Taylor scheme")
SeaIce.setTimeStep!(sim, epsilon=0.07)
tol = 0.1
info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
-SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
- contact_tangential_rheology="None", verbose=verbose)
+SeaIce.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbose)
E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
t@@ -43,8 +42,7 @@ sim = deepcopy(sim_init)
SeaIce.setTimeStep!(sim, epsilon=0.007)
tol = 0.01
info("Relative tolerance: $(tol*100.)%")
-SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
- contact_tangential_rheology="None", verbose=verbose)
+SeaIce.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbose)
E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
t@@ -57,8 +55,7 @@ sim = deepcopy(sim_init)
SeaIce.setTimeStep!(sim, epsilon=0.07)
tol = 0.01
info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
-SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
- contact_tangential_rheology="None", verbose=verbose)
+SeaIce.run!(sim, temporal_integration_method="Three-term Taylor", verbose=verbose)
E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
t@@ -85,8 +82,7 @@ info("Testing kinetic energy conservation with Two-term Taylor scheme")
SeaIce.setTimeStep!(sim, epsilon=0.07)
tol = 0.1
info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
-SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
- contact_tangential_rheology="None", verbose=verbose)
+SeaIce.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbose)
E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
t@@ -99,8 +95,7 @@ sim = deepcopy(sim_init)
SeaIce.setTimeStep!(sim, epsilon=0.007)
tol = 0.01
info("Relative tolerance: $(tol*100.)%")
-SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
- contact_tangential_rheology="None", verbose=verbose)
+SeaIce.run!(sim, temporal_integration_method="Two-term Taylor", verbose=verbose)
E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
t@@ -114,7 +109,7 @@ SeaIce.setTimeStep!(sim, epsilon=0.07)
tol = 0.01
info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
- contact_tangential_rheology="None", verbose=verbose)
+ verbose=verbose)
E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
t@@ -123,85 +118,19 @@ E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
info("## Contact-normal elasticity and tangential viscosity and friction")
-info("Testing gamma_t = 0.")
-sim = SeaIce.createSimulation(id="test")
-SeaIce.addIceFloeCylindrical(sim, [0., 10.], 10., 1., verbose=verbose)
-SeaIce.addIceFloeCylindrical(sim, [19.0, 0.], 10., 1., verbose=verbose)
-sim.ice_floes[1].lin_vel[1] = 0.1
-
-E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
-E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
-
-# With decreasing timestep (epsilon towards 0), the explicit integration scheme
-# should become more correct
-
SeaIce.setTotalTime!(sim, 30.0)
-sim_init = deepcopy(sim)
-
-info("Testing kinetic energy conservation with Two-term Taylor scheme")
-SeaIce.setTimeStep!(sim, epsilon=0.07)
-tol = 0.1
-info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
-E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
-E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
-SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
- contact_tangential_rheology="Linear Viscous Frictional",
- verbose=verbose)
-
-@test sim.ice_floes[1].ang_pos ≈ 0.
-@test sim.ice_floes[1].ang_vel ≈ 0.
-@test sim.ice_floes[2].ang_pos ≈ 0.
-E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
-E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
-@test_approx_eq_eps E_kin_lin_init E_kin_lin_final E_kin_lin_init*tol
-@test_approx_eq E_kin_rot_init E_kin_rot_final
-
-info("Testing kinetic energy conservation with Three-term Taylor scheme")
+sim_init.ice_floes[1].contact_viscosity_tangential = 1e6
+sim_init.ice_floes[2].contact_viscosity_tangential = 1e6
+sim_init.ice_floes[1].contact_dynamic_friction = 1e2 # no Coulomb sliding
+sim_init.ice_floes[2].contact_dynamic_friction = 1e2 # no Coulomb sliding
+sim_init.ice_floes[2].fixed = true
sim = deepcopy(sim_init)
-SeaIce.setTimeStep!(sim, epsilon=0.07)
-tol = 0.01
-info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
-E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
-E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
-SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
- contact_tangential_rheology="Linear Viscous Frictional",
- verbose=verbose)
-
-@test sim.ice_floes[1].ang_pos ≈ 0.
-@test sim.ice_floes[1].ang_vel ≈ 0.
-@test sim.ice_floes[2].ang_pos ≈ 0.
-E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
-E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
-@test_approx_eq_eps E_kin_lin_init E_kin_lin_final E_kin_lin_init*tol
-@test_approx_eq E_kin_rot_init E_kin_rot_final
-
-
-info("## Contact-normal elasticity and tangential viscosity and friction")
-info("# One ice floe fixed")
-sim = SeaIce.createSimulation(id="test")
-SeaIce.addIceFloeCylindrical(sim, [0., 10.], 10., 1., verbose=verbose)
-SeaIce.addIceFloeCylindrical(sim, [19., 0.], 10., 1., verbose=verbose)
-sim.ice_floes[1].lin_vel[1] = 0.1
-sim.ice_floes[2].fixed = true
-sim.ice_floes[1].contact_viscosity_tangential = 1e4
-sim.ice_floes[2].contact_viscosity_tangential = 1e4
-
-E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
-E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
-
-# With decreasing timestep (epsilon towards 0), the explicit integration scheme
-# should become more correct
-
-SeaIce.setTotalTime!(sim, 30.0)
-#sim.file_time_step = 1.
-sim_init = deepcopy(sim)
info("Testing kinetic energy conservation with Two-term Taylor scheme")
SeaIce.setTimeStep!(sim, epsilon=0.07)
tol = 0.1
info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
- contact_tangential_rheology="Linear Viscous Frictional",
verbose=verbose)
@test sim.ice_floes[1].ang_pos < 0.
t@@ -221,7 +150,6 @@ info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
- contact_tangential_rheology="Linear Viscous Frictional",
verbose=verbose)
@test sim.ice_floes[1].ang_pos ≈ 0.
t@@ -235,10 +163,9 @@ E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
info("Testing kinetic energy conservation with Two-term Taylor scheme")
sim = deepcopy(sim_init)
SeaIce.setTimeStep!(sim, epsilon=0.007)
-tol = 0.02
+tol = 0.1
info("Relative tolerance: $(tol*100.)%")
SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
- contact_tangential_rheology="Linear Viscous Frictional",
verbose=verbose)
@test sim.ice_floes[1].ang_pos < 0.
t@@ -253,10 +180,9 @@ E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
info("Testing kinetic energy conservation with Three-term Taylor scheme")
sim = deepcopy(sim_init)
SeaIce.setTimeStep!(sim, epsilon=0.07)
-tol = 0.03
+tol = 0.09
info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
- contact_tangential_rheology="Linear Viscous Frictional",
verbose=verbose)
@test sim.ice_floes[1].ang_pos < 0.
t@@ -290,7 +216,6 @@ SeaIce.setTimeStep!(sim, epsilon=0.07)
tol = 0.1
info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
- contact_tangential_rheology="Linear Viscous Frictional",
verbose=verbose)
@test sim.ice_floes[1].ang_pos < 0.
t@@ -307,7 +232,6 @@ SeaIce.setTimeStep!(sim, epsilon=0.007)
tol = 0.02
info("Relative tolerance: $(tol*100.)%")
SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
- contact_tangential_rheology="Linear Viscous Frictional",
verbose=verbose)
E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
t@@ -321,7 +245,6 @@ SeaIce.setTimeStep!(sim, epsilon=0.07)
tol = 0.02
info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
- contact_tangential_rheology="Linear Viscous Frictional",
verbose=verbose)
@test sim.ice_floes[1].ang_pos < 0.
t@@ -356,7 +279,6 @@ SeaIce.setTimeStep!(sim, epsilon=0.07)
tol = 0.1
info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
- contact_tangential_rheology="Linear Viscous Frictional",
verbose=verbose)
@test sim.ice_floes[1].ang_pos > 0.
t@@ -373,7 +295,6 @@ SeaIce.setTimeStep!(sim, epsilon=0.007)
tol = 0.02
info("Relative tolerance: $(tol*100.)%")
SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
- contact_tangential_rheology="Linear Viscous Frictional",
verbose=verbose)
E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
t@@ -387,7 +308,6 @@ SeaIce.setTimeStep!(sim, epsilon=0.07)
tol = 0.02
info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
- contact_tangential_rheology="Linear Viscous Frictional",
verbose=verbose)
@test sim.ice_floes[1].ang_pos > 0.
t@@ -422,7 +342,6 @@ SeaIce.setTimeStep!(sim, epsilon=0.07)
tol = 0.1
info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
- contact_tangential_rheology="Linear Viscous Frictional",
verbose=verbose)
@test sim.ice_floes[1].ang_pos < 0.
t@@ -439,7 +358,6 @@ SeaIce.setTimeStep!(sim, epsilon=0.007)
tol = 0.02
info("Relative tolerance: $(tol*100.)%")
SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
- contact_tangential_rheology="Linear Viscous Frictional",
verbose=verbose)
E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
t@@ -453,7 +371,6 @@ SeaIce.setTimeStep!(sim, epsilon=0.07)
tol = 0.02
info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
- contact_tangential_rheology="Linear Viscous Frictional",
verbose=verbose)
@test sim.ice_floes[1].ang_pos < 0.
t@@ -463,3 +380,231 @@ SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
@test_approx_eq_eps E_kin_lin_init+E_kin_rot_init E_kin_lin_final+E_kin_rot_final E_kin_lin_init*tol
+
+
+info("# Tangential elasticity, no tangential viscosity, no Coulomb slip")
+
+sim = SeaIce.createSimulation(id="test")
+SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., verbose=verbose)
+SeaIce.addIceFloeCylindrical(sim, [19.0, -10.], 10., 1., verbose=verbose)
+sim.ice_floes[2].lin_vel[1] = -0.1
+sim.ice_floes[1].contact_dynamic_friction = 1e3 # disable Coulomb slip
+sim.ice_floes[2].contact_dynamic_friction = 1e3 # disable Coulomb slip
+sim.ice_floes[1].contact_viscosity_tangential = 0. # disable tan. viscosity
+sim.ice_floes[2].contact_viscosity_tangential = 0. # disable tan. viscosity
+sim.ice_floes[1].contact_stiffness_tangential =
+ sim.ice_floes[1].contact_stiffness_normal # enable tangential elasticity
+sim.ice_floes[2].contact_stiffness_tangential =
+ sim.ice_floes[2].contact_stiffness_normal # enable tangential elasticity
+
+E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+
+# With decreasing timestep (epsilon towards 0), the explicit integration scheme
+# should become more correct
+
+SeaIce.setTotalTime!(sim, 30.0)
+sim_init = deepcopy(sim)
+
+info("Testing kinetic energy conservation with Two-term Taylor scheme")
+SeaIce.setTimeStep!(sim, epsilon=0.07)
+tol = 0.1
+info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
+SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
+ verbose=verbose)
+
+@test sim.ice_floes[1].ang_pos < 0.
+@test sim.ice_floes[1].ang_vel < 0.
+@test sim.ice_floes[2].ang_pos < 0.
+@test sim.ice_floes[2].ang_vel < 0.
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test_approx_eq_eps E_kin_lin_init+E_kin_rot_init E_kin_lin_final+E_kin_rot_final E_kin_lin_init*tol
+
+info("Testing kinetic energy conservation with Two-term Taylor scheme")
+sim = deepcopy(sim_init)
+SeaIce.setTimeStep!(sim, epsilon=0.007)
+tol = 0.02
+info("Relative tolerance: $(tol*100.)%")
+SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
+ verbose=verbose)
+
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test_approx_eq_eps E_kin_lin_init+E_kin_rot_init E_kin_lin_final+E_kin_rot_final E_kin_lin_init*tol
+
+
+info("Testing kinetic energy conservation with Three-term Taylor scheme")
+sim = deepcopy(sim_init)
+SeaIce.setTimeStep!(sim, epsilon=0.07)
+tol = 0.03
+info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
+SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
+ verbose=verbose)
+
+@test sim.ice_floes[1].ang_pos < 0.
+@test sim.ice_floes[1].ang_vel < 0.
+@test sim.ice_floes[2].ang_pos < 0.
+@test sim.ice_floes[2].ang_vel < 0.
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test_approx_eq_eps E_kin_lin_init+E_kin_rot_init E_kin_lin_final+E_kin_rot_final E_kin_lin_init*tol
+
+
+info("# Tangential elasticity, no tangential viscosity, Coulomb slip")
+
+sim = SeaIce.createSimulation(id="test")
+SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., verbose=verbose)
+SeaIce.addIceFloeCylindrical(sim, [19.0, -10.], 10., 1., verbose=verbose)
+sim.ice_floes[2].lin_vel[1] = -0.1
+sim.ice_floes[1].contact_dynamic_friction = 0.1 # enable Coulomb slip
+sim.ice_floes[2].contact_dynamic_friction = 0.1 # enable Coulomb slip
+sim.ice_floes[1].contact_viscosity_tangential = 0. # disable tan. viscosity
+sim.ice_floes[2].contact_viscosity_tangential = 0. # disable tan. viscosity
+sim.ice_floes[1].contact_stiffness_tangential =
+ sim.ice_floes[1].contact_stiffness_normal # enable tangential elasticity
+sim.ice_floes[2].contact_stiffness_tangential =
+ sim.ice_floes[2].contact_stiffness_normal # enable tangential elasticity
+
+E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+
+# With decreasing timestep (epsilon towards 0), the explicit integration scheme
+# should become more correct
+
+SeaIce.setTotalTime!(sim, 30.0)
+sim_init = deepcopy(sim)
+
+info("Testing kinetic energy conservation with Two-term Taylor scheme")
+sim = deepcopy(sim_init)
+SeaIce.setTimeStep!(sim, epsilon=0.007)
+tol = 0.02
+info("Relative tolerance: $(tol*100.)%")
+SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
+ verbose=verbose)
+
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test E_kin_lin_init+E_kin_rot_init > E_kin_lin_final+E_kin_rot_final
+
+info("Testing kinetic energy conservation with Three-term Taylor scheme")
+sim = deepcopy(sim_init)
+SeaIce.setTimeStep!(sim, epsilon=0.07)
+tol = 0.03
+info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
+SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
+ verbose=verbose)
+
+@test sim.ice_floes[1].ang_pos < 0.
+@test sim.ice_floes[1].ang_vel < 0.
+@test sim.ice_floes[2].ang_pos < 0.
+@test sim.ice_floes[2].ang_vel < 0.
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test E_kin_lin_init+E_kin_rot_init > E_kin_lin_final+E_kin_rot_final
+
+
+info("# Tangential elasticity, tangential viscosity, no Coulomb slip")
+
+sim = SeaIce.createSimulation(id="test")
+SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., verbose=verbose)
+SeaIce.addIceFloeCylindrical(sim, [19.0, -10.], 10., 1., verbose=verbose)
+sim.ice_floes[2].lin_vel[1] = -0.1
+sim.ice_floes[1].contact_dynamic_friction = 1e3 # disable Coulomb slip
+sim.ice_floes[2].contact_dynamic_friction = 1e3 # disable Coulomb slip
+sim.ice_floes[1].contact_viscosity_tangential = 1e4 # enable tan. viscosity
+sim.ice_floes[2].contact_viscosity_tangential = 1e4 # enable tan. viscosity
+sim.ice_floes[1].contact_stiffness_tangential =
+ sim.ice_floes[1].contact_stiffness_normal # enable tangential elasticity
+sim.ice_floes[2].contact_stiffness_tangential =
+ sim.ice_floes[2].contact_stiffness_normal # enable tangential elasticity
+
+E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+
+# With decreasing timestep (epsilon towards 0), the explicit integration scheme
+# should become more correct
+
+SeaIce.setTotalTime!(sim, 30.0)
+sim_init = deepcopy(sim)
+
+info("Testing kinetic energy conservation with Two-term Taylor scheme")
+sim = deepcopy(sim_init)
+SeaIce.setTimeStep!(sim, epsilon=0.007)
+tol = 0.02
+info("Relative tolerance: $(tol*100.)%")
+SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
+ verbose=verbose)
+
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test E_kin_lin_init+E_kin_rot_init > E_kin_lin_final+E_kin_rot_final
+
+info("Testing kinetic energy conservation with Three-term Taylor scheme")
+sim = deepcopy(sim_init)
+SeaIce.setTimeStep!(sim, epsilon=0.07)
+tol = 0.03
+info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
+SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
+ verbose=verbose)
+
+@test sim.ice_floes[1].ang_pos < 0.
+@test sim.ice_floes[1].ang_vel < 0.
+@test sim.ice_floes[2].ang_pos < 0.
+@test sim.ice_floes[2].ang_vel < 0.
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test E_kin_lin_init+E_kin_rot_init > E_kin_lin_final+E_kin_rot_final
+
+
+info("# Tangential elasticity, tangential viscosity, Coulomb slip")
+
+sim = SeaIce.createSimulation(id="test")
+SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., verbose=verbose)
+SeaIce.addIceFloeCylindrical(sim, [19.0, -10.], 10., 1., verbose=verbose)
+sim.ice_floes[2].lin_vel[1] = -0.1
+sim.ice_floes[1].contact_dynamic_friction = 0.1 # enable Coulomb slip
+sim.ice_floes[2].contact_dynamic_friction = 0.1 # enable Coulomb slip
+sim.ice_floes[1].contact_viscosity_tangential = 1e4 # enable tan. viscosity
+sim.ice_floes[2].contact_viscosity_tangential = 1e4 # enable tan. viscosity
+sim.ice_floes[1].contact_stiffness_tangential =
+ sim.ice_floes[1].contact_stiffness_normal # enable tangential elasticity
+sim.ice_floes[2].contact_stiffness_tangential =
+ sim.ice_floes[2].contact_stiffness_normal # enable tangential elasticity
+
+E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+
+# With decreasing timestep (epsilon towards 0), the explicit integration scheme
+# should become more correct
+
+SeaIce.setTotalTime!(sim, 30.0)
+sim_init = deepcopy(sim)
+
+info("Testing kinetic energy conservation with Two-term Taylor scheme")
+sim = deepcopy(sim_init)
+SeaIce.setTimeStep!(sim, epsilon=0.007)
+tol = 0.02
+info("Relative tolerance: $(tol*100.)%")
+SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
+ verbose=verbose)
+
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test E_kin_lin_init+E_kin_rot_init > E_kin_lin_final+E_kin_rot_final
+
+info("Testing kinetic energy conservation with Three-term Taylor scheme")
+sim = deepcopy(sim_init)
+SeaIce.setTimeStep!(sim, epsilon=0.07)
+tol = 0.03
+info("Relative tolerance: $(tol*100.)% with time step: $(sim.time_step)")
+SeaIce.run!(sim, temporal_integration_method="Three-term Taylor",
+ verbose=verbose)
+
+@test sim.ice_floes[1].ang_pos < 0.
+@test sim.ice_floes[1].ang_vel < 0.
+@test sim.ice_floes[2].ang_pos < 0.
+@test sim.ice_floes[2].ang_vel < 0.
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+@test E_kin_lin_init+E_kin_rot_init > E_kin_lin_final+E_kin_rot_final
(DIR) diff --git a/test/contact-search-and-geometry.jl b/test/contact-search-and-geometry.jl
t@@ -20,24 +20,40 @@ info("Testing findContactsAllToAll(...)")
sim_copy = deepcopy(sim)
SeaIce.findContactsAllToAll!(sim)
-@test 1 == length(sim.contact_pairs)
-@test_approx_eq [1, 2] sim.contact_pairs[1]
-
info("Testing findContacts(...)")
sim = deepcopy(sim_copy)
SeaIce.findContacts!(sim)
sim.ice_floes[1].fixed = true
-@test 1 == length(sim.contact_pairs)
-@test_approx_eq [1, 2] sim.contact_pairs[1]
+# The contact should be registered in ice floe 1, but not ice floe 2
+@test 2 == sim.ice_floes[1].contacts[1]
+for ic=2:32
+ @test 0 == sim.ice_floes[1].contacts[ic]
+ @test [0., 0.] ≈ sim.ice_floes[1].contact_parallel_displacement[ic]
+end
+for ic=1:32
+ @test 0 == sim.ice_floes[2].contacts[ic]
+ @test [0., 0.] ≈ sim.ice_floes[2].contact_parallel_displacement[ic]
+end
+@test 1 == sim.ice_floes[1].n_contacts
+@test 1 == sim.ice_floes[2].n_contacts
info("Testing findContacts(...)")
sim = deepcopy(sim_copy)
SeaIce.findContacts!(sim)
-@test 1 == length(sim.contact_pairs)
-@test_approx_eq [1, 2] sim.contact_pairs[1]
+@test 2 == sim.ice_floes[1].contacts[1]
+for ic=2:32
+ @test 0 == sim.ice_floes[1].contacts[ic]
+ @test [0., 0.] ≈ sim.ice_floes[1].contact_parallel_displacement[ic]
+end
+for ic=1:32
+ @test 0 == sim.ice_floes[2].contacts[ic]
+ @test [0., 0.] ≈ sim.ice_floes[2].contact_parallel_displacement[ic]
+end
+@test 1 == sim.ice_floes[1].n_contacts
+@test 1 == sim.ice_floes[2].n_contacts
@test_throws ErrorException SeaIce.findContacts!(sim, method="")
t@@ -45,18 +61,47 @@ sim = deepcopy(sim_copy)
sim.ice_floes[1].fixed = true
sim.ice_floes[2].fixed = true
SeaIce.findContacts!(sim)
-@test 0 == length(sim.contact_pairs)
+for ic=1:32
+ @test 0 == sim.ice_floes[1].contacts[ic]
+ @test [0., 0.] ≈ sim.ice_floes[1].contact_parallel_displacement[ic]
+end
+for ic=1:32
+ @test 0 == sim.ice_floes[2].contacts[ic]
+ @test [0., 0.] ≈ sim.ice_floes[2].contact_parallel_displacement[ic]
+end
+@test 0 == sim.ice_floes[1].n_contacts
+@test 0 == sim.ice_floes[2].n_contacts
+
sim = deepcopy(sim_copy)
SeaIce.disableIceFloe!(sim, 1)
SeaIce.findContacts!(sim)
-@test 0 == length(sim.contact_pairs)
+for ic=1:32
+ @test 0 == sim.ice_floes[1].contacts[ic]
+ @test [0., 0.] ≈ sim.ice_floes[1].contact_parallel_displacement[ic]
+end
+for ic=1:32
+ @test 0 == sim.ice_floes[2].contacts[ic]
+ @test [0., 0.] ≈ sim.ice_floes[2].contact_parallel_displacement[ic]
+end
+@test 0 == sim.ice_floes[1].n_contacts
+@test 0 == sim.ice_floes[2].n_contacts
+
sim = deepcopy(sim_copy)
SeaIce.disableIceFloe!(sim, 1)
SeaIce.disableIceFloe!(sim, 2)
SeaIce.findContacts!(sim)
-@test 0 == length(sim.contact_pairs)
+for ic=1:32
+ @test 0 == sim.ice_floes[1].contacts[ic]
+ @test [0., 0.] ≈ sim.ice_floes[1].contact_parallel_displacement[ic]
+end
+for ic=1:32
+ @test 0 == sim.ice_floes[2].contacts[ic]
+ @test [0., 0.] ≈ sim.ice_floes[2].contact_parallel_displacement[ic]
+end
+@test 0 == sim.ice_floes[1].n_contacts
+@test 0 == sim.ice_floes[2].n_contacts
info("Testing if interact(...) removes contacts correctly")
sim = deepcopy(sim_copy)
t@@ -64,8 +109,18 @@ SeaIce.findContacts!(sim)
SeaIce.interact!(sim)
SeaIce.findContacts!(sim)
-@test 1 == length(sim.contact_pairs)
-@test_approx_eq [1, 2] sim.contact_pairs[1]
+@test 2 == sim.ice_floes[1].contacts[1]
+for ic=2:32
+ @test 0 == sim.ice_floes[1].contacts[ic]
+ @test [0., 0.] ≈ sim.ice_floes[1].contact_parallel_displacement[ic]
+end
+for ic=1:32
+ @test 0 == sim.ice_floes[2].contacts[ic]
+ @test [0., 0.] ≈ sim.ice_floes[2].contact_parallel_displacement[ic]
+end
+@test 1 == sim.ice_floes[1].n_contacts
+@test 1 == sim.ice_floes[2].n_contacts
+
info("Testing findContactsOceanGrid(...)")
sim = deepcopy(sim_copy)
t@@ -73,8 +128,18 @@ sim.ocean = SeaIce.createRegularOceanGrid([4, 4, 2], [80., 80., 2.])
SeaIce.sortIceFloesInOceanGrid!(sim)
SeaIce.findContactsOceanGrid!(sim)
-@test 1 == length(sim.contact_pairs)
-@test_approx_eq [1, 2] sim.contact_pairs[1]
+@test 2 == sim.ice_floes[1].contacts[1]
+for ic=2:32
+ @test 0 == sim.ice_floes[1].contacts[ic]
+ @test [0., 0.] ≈ sim.ice_floes[1].contact_parallel_displacement[ic]
+end
+for ic=1:32
+ @test 0 == sim.ice_floes[2].contacts[ic]
+ @test [0., 0.] ≈ sim.ice_floes[2].contact_parallel_displacement[ic]
+end
+@test 1 == sim.ice_floes[1].n_contacts
+@test 1 == sim.ice_floes[2].n_contacts
+
sim = deepcopy(sim_copy)
sim.ocean = SeaIce.createRegularOceanGrid([4, 4, 2], [80., 80., 2.])
t@@ -82,8 +147,18 @@ sim.ice_floes[1].fixed = true
SeaIce.sortIceFloesInOceanGrid!(sim)
SeaIce.findContactsOceanGrid!(sim)
-@test 1 == length(sim.contact_pairs)
-@test_approx_eq [1, 2] sim.contact_pairs[1]
+@test 2 == sim.ice_floes[1].contacts[1]
+for ic=2:32
+ @test 0 == sim.ice_floes[1].contacts[ic]
+ @test [0., 0.] ≈ sim.ice_floes[1].contact_parallel_displacement[ic]
+end
+for ic=1:32
+ @test 0 == sim.ice_floes[2].contacts[ic]
+ @test [0., 0.] ≈ sim.ice_floes[2].contact_parallel_displacement[ic]
+end
+@test 1 == sim.ice_floes[1].n_contacts
+@test 1 == sim.ice_floes[2].n_contacts
+
sim = deepcopy(sim_copy)
sim.ocean = SeaIce.createRegularOceanGrid([4, 4, 2], [80., 80., 2.])
t@@ -92,7 +167,16 @@ sim.ice_floes[2].fixed = true
SeaIce.sortIceFloesInOceanGrid!(sim)
SeaIce.findContactsOceanGrid!(sim)
-@test 0 == length(sim.contact_pairs)
+for ic=1:32
+ @test 0 == sim.ice_floes[1].contacts[ic]
+ @test [0., 0.] ≈ sim.ice_floes[1].contact_parallel_displacement[ic]
+end
+for ic=1:32
+ @test 0 == sim.ice_floes[2].contacts[ic]
+ @test [0., 0.] ≈ sim.ice_floes[2].contact_parallel_displacement[ic]
+end
+@test 0 == sim.ice_floes[1].n_contacts
+@test 0 == sim.ice_floes[2].n_contacts
info("Testing findContacts(...)")
sim = deepcopy(sim_copy)
t@@ -100,10 +184,16 @@ sim.ocean = SeaIce.createRegularOceanGrid([4, 4, 2], [80., 80., 2.])
SeaIce.sortIceFloesInOceanGrid!(sim)
SeaIce.findContacts!(sim)
-@test 1 == length(sim.contact_pairs)
-@test_approx_eq [1, 2] sim.contact_pairs[1]
-
-@test 1 == length(sim.contact_pairs)
-@test_approx_eq [1, 2] sim.contact_pairs[1]
+@test 2 == sim.ice_floes[1].contacts[1]
+for ic=2:32
+ @test 0 == sim.ice_floes[1].contacts[ic]
+ @test [0., 0.] ≈ sim.ice_floes[1].contact_parallel_displacement[ic]
+end
+for ic=1:32
+ @test 0 == sim.ice_floes[2].contacts[ic]
+ @test [0., 0.] ≈ sim.ice_floes[2].contact_parallel_displacement[ic]
+end
+@test 1 == sim.ice_floes[1].n_contacts
+@test 1 == sim.ice_floes[2].n_contacts
@test_throws ErrorException SeaIce.findContacts!(sim, method="")
(DIR) diff --git a/test/vtk.jl b/test/vtk.jl
t@@ -24,7 +24,7 @@ else
end
@test readstring(`$(cmd) test.icefloes.1.vtu$(cmd_post)`) ==
-"72f4e4b854d7e92afd8cde0b79a4af6a29e49714b751ffc30a4ff3867f44b505 test.icefloes.1.vtu\n"
+"a01d322026a56b1332c2174e4b513015c63ad44e2a28140bd2c2cccf7df38a13 test.icefloes.1.vtu\n"
@test readstring(`$(cmd) test.ocean.1.vts$(cmd_post)`) ==
"f0117e414c4e71a0c55980f63865eb03b6c597fa2546983258b8a57eb4ff2a25 test.ocean.1.vts\n"