commit e667d41a9b5f93780dc7d493620590183acdec10
parent 8c6f1f493f21bb6512ce390951af2aa1e3d33bfd
Author: Anders Damsgaard <andersd@riseup.net>
Date: Tue, 2 May 2017 15:48:20 -0400
check for nan in harmonic mean operation, return 0 if true. Add more tests
Diffstat:
4 files changed, 90 insertions(+), 12 deletions(-)
diff --git a/examples/nares_strait.jl b/examples/nares_strait.jl
@@ -116,7 +116,8 @@ 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)
+ SeaIce.run!(sim, status_interval=1, single_step=true,
+ contact_tangential_rheology="Linear Viscous Frictional")
end
for i=1:size(sim.ocean.xh, 1)
if sim.ocean.ice_floe_list[i, end] == []
@@ -137,7 +138,3 @@ while sim.time < sim.time_total
end
end
end
-
-
-
-
diff --git a/src/interaction.jl b/src/interaction.jl
@@ -122,7 +122,11 @@ function interactTangentialLinearViscousFrictional(simulation::Simulation,
end
function harmonicMean(a::Any, b::Any)
- return 2.*a*b/(a + b)
+ hm = 2.*a*b/(a + b)
+ if isnan(hm)
+ return 0.
+ end
+ return hm
end
function findTorque(simulation::Simulation, overlap_vector::vector,
diff --git a/src/temporal_integration.jl b/src/temporal_integration.jl
@@ -73,7 +73,7 @@ freedom for an `icefloe`.
"""
function updateIceFloeKinematicsThreeTermTaylor!(icefloe::IceFloeCylindrical,
simulation::Simulation)
-
+
if simulation.time_iteration == 0
lin_acc_0 = zeros(2)
ang_acc_0 = 0.
diff --git a/test/collision-2floes-oblique.jl b/test/collision-2floes-oblique.jl
@@ -123,12 +123,68 @@ 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.2
+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 = deepcopy(sim_init)
+SeaIce.setTimeStep!(sim, epsilon=0.07)
+tol = 0.02
+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)
@@ -142,7 +198,7 @@ sim_init = deepcopy(sim)
info("Testing kinetic energy conservation with Two-term Taylor scheme")
SeaIce.setTimeStep!(sim, epsilon=0.07)
-tol = 0.5
+tol = 0.2
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",
@@ -156,11 +212,30 @@ 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("mu_d = 0.")
+sim = deepcopy(sim_init)
+sim.ice_floes[1].contact_dynamic_friction = 0.
+SeaIce.setTimeStep!(sim, epsilon=0.07)
+tol = 0.02
+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("Testing kinetic energy conservation with Two-term Taylor scheme")
sim = deepcopy(sim_init)
SeaIce.setTimeStep!(sim, epsilon=0.007)
-tol = 0.2
+tol = 0.06
info("Relative tolerance: $(tol*100.)%")
SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
contact_tangential_rheology="Linear Viscous Frictional",
@@ -178,7 +253,7 @@ 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.2
+tol = 0.07
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",
@@ -198,6 +273,8 @@ 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
+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)
@@ -227,7 +304,7 @@ 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.1
+tol = 0.05
info("Relative tolerance: $(tol*100.)%")
SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
contact_tangential_rheology="Linear Viscous Frictional",
@@ -241,7 +318,7 @@ 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.1
+tol = 0.05
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",
