commit 784fc626abb1dde66a6aed241677ad8d5d8b53ac
parent bb21ec79edd4b21786d862d659daca6b2cbacc34
Author: Anders Damsgaard <andersd@riseup.net>
Date: Mon, 8 May 2017 16:57:16 -0400
test additional collisions, rotational tests will fail
Diffstat:
4 files changed, 275 insertions(+), 51 deletions(-)
diff --git a/src/contact_search.jl b/src/contact_search.jl
@@ -30,9 +30,20 @@ function findContacts!(simulation::Simulation;
end
export interIceFloePositionVector
+"""
+ interIceFloePositionVector(simulation, i, j)
+
+Returns a `vector` pointing from ice floe `i` to ice floe `j` in the
+`simulation`.
+
+# Arguments
+* `simulation::Simulation`: the simulation object containing the ice floes.
+* `i::Int`: index of the first ice floe.
+* `j::Int`: index of the second ice floe.
+"""
function interIceFloePositionVector(simulation::Simulation,
i::Integer, j::Integer)
- return simulation.ice_floes[j].lin_pos - simulation.ice_floes[i].lin_pos
+ return simulation.ice_floes[i].lin_pos - simulation.ice_floes[j].lin_pos
end
"""
@@ -95,20 +106,24 @@ end
export addContact!
"""
-checkAndAddContact!(simulation, i, j)
+ checkAndAddContact!(simulation, i, j)
-Check for contact between two ice floes and register the interaction.
+Check for contact between two ice floes and register the interaction in the
+`simulation` object. The indexes of the two ice floes is stored in
+`simulation.contact_pairs` as `[i, j]`. The overlap vector is parallel to a
+straight line connecting the ice floe centers, points away from ice floe `i` and
+towards `j`, and is stored in `simulation.overlaps`. A zero-length vector is
+written to `simulation.contact_parallel_displacement`.
# Arguments
-* `simulation::Simulation`: the simulation object containing the ice floes
-* `i::Int`: index of the first ice floe
-* `j::Int`: index of the second ice floe
+* `simulation::Simulation`: the simulation object containing the ice floes.
+* `i::Int`: index of the first ice floe.
+* `j::Int`: index of the second ice floe.
"""
function checkAndAddContact!(simulation::Simulation, i::Int, j::Int)
if i < j
- if (simulation.ice_floes[i].fixed &&
- simulation.ice_floes[j].fixed) ||
+ if (simulation.ice_floes[i].fixed && simulation.ice_floes[j].fixed) ||
!simulation.ice_floes[i].enabled ||
!simulation.ice_floes[j].enabled
return
diff --git a/src/interaction.jl b/src/interaction.jl
@@ -11,11 +11,11 @@ function interact!(simulation::Simulation;
# IceFloe to grain collisions
while !isempty(simulation.contact_pairs)
contact_pair = pop!(simulation.contact_pairs)
- overlap_vector = pop!(simulation.overlaps)
+ overlap = pop!(simulation.overlaps)
contact_parallel_displacement =
pop!(simulation.contact_parallel_displacement)
interactIceFloes!(simulation, contact_pair[1], contact_pair[2],
- overlap_vector, contact_parallel_displacement,
+ overlap, contact_parallel_displacement,
contact_normal_rheology=contact_normal_rheology,
contact_tangential_rheology=
contact_tangential_rheology)
@@ -30,8 +30,8 @@ function adds the compressive force of the interaction to the ice floe
"""
function interactIceFloes!(simulation::Simulation,
i::Int, j::Int,
- overlap_vector::Array{Float64, 1},
- contact_parallel_displacement::Array{Float64, 1};
+ overlap::vector,
+ contact_parallel_displacement::vector;
contact_normal_rheology::String = "Linear Elastic",
contact_tangential_rheology::String = "None")
@@ -41,7 +41,7 @@ function interactIceFloes!(simulation::Simulation,
if contact_normal_rheology == "None"
# do nothing
elseif contact_normal_rheology == "Linear Elastic"
- force_n = interactNormalLinearElastic(simulation, i, j, overlap_vector)
+ force_n = interactNormalLinearElastic(simulation, i, j, overlap)
else
error("Unknown contact_normal_rheology '$contact_normal_rheology'")
end
@@ -50,7 +50,7 @@ function interactIceFloes!(simulation::Simulation,
# do nothing
elseif contact_tangential_rheology == "Linear Viscous Frictional"
force_t = interactTangentialLinearViscousFrictional(simulation, i, j,
- overlap_vector,
+ overlap,
force_n)
else
error("Unknown contact_tangential_rheology ",
@@ -61,7 +61,7 @@ function interactIceFloes!(simulation::Simulation,
simulation.ice_floes[j].force -= force_n + force_t;
if norm(force_t) > 0.
- torque = -findTorque(simulation, overlap_vector, force_t, i, j)
+ torque = -findTorque(simulation, overlap, force_t, i, j)
simulation.ice_floes[i].torque += torque
simulation.ice_floes[j].torque += torque
end
@@ -79,46 +79,108 @@ direction.
"""
function interactNormalLinearElastic(simulation::Simulation,
i::Int, j::Int,
- overlap_vector::vector)
+ overlap::vector)
k_n_harmonic_mean =
harmonicMean(simulation.ice_floes[i].contact_stiffness_normal,
simulation.ice_floes[j].contact_stiffness_normal)
- return k_n_harmonic_mean * overlap_vector
+ return -k_n_harmonic_mean*overlap
end
export interactTangentialLinearElastic
"""
-Resolves linear-elastic interaction between two ice floes in the
-contact-parallel (tangential) direction.
+ interactTangentialLinearViscousFrictional(simulation, i, j,
+ overlap, force_n)
+
+Resolves linear-viscous interaction between two ice floes in the contact-
+parallel (tangential) direction, with a frictional limit dependent on the
+Coulomb criterion.
+
+# Arguments
+* `simulation::Simulation`: the simulation object containing the ice floes.
+* `i::Int`: index of the first ice floe.
+* `j::Int`: index of the second ice floe.
+* `overlap::vector`: two-dimensional vector pointing from i to j.
+* `force_n::vector`: normal force from the interaction.
"""
function interactTangentialLinearViscousFrictional(simulation::Simulation,
i::Int, j::Integer,
- overlap_vector::vector,
+ overlap::vector,
force_n::vector)
+ """
contact_parallel_velocity = findContactParallelVelocity(simulation, i, j,
- overlap_vector)
+ overlap)
- if norm(contact_parallel_velocity) ≈ 0.
+
+
+ if contact_parallel_velocity ≈ 0.
return [0., 0.]
end
-
gamma_t_harmonic_mean = harmonicMean(
simulation.ice_floes[i].contact_viscosity_tangential,
simulation.ice_floes[j].contact_viscosity_tangential)
+ mu_d_minimum = min(simulation.ice_floes[i].contact_dynamic_friction,
+ simulation.ice_floes[j].contact_dynamic_friction)
- if norm(gamma_t_harmonic_mean) ≈ 0.
+ if gamma_t_harmonic_mean ≈ 0. || mu_d_minimum ≈ 0.
return [0., 0.]
end
- mu_d_minimum = min(simulation.ice_floes[i].contact_dynamic_friction,
- simulation.ice_floes[j].contact_dynamic_friction)
+ force_t = abs(gamma_t_harmonic_mean*contact_parallel_velocity)
+ if force_t > mu_d_minimum*norm(force_n)
+ force_t = mu_d_minimum*norm(force_n)
+ end
+ if contact_parallel_velocity > 0.
+ force_t = -force_t
+ end
+
+ return force_t*contactParallelVector(contactNormalVector(overlap))
+ """
+
+ p = simulation.ice_floes[i].lin_pos - simulation.ice_floes[j].lin_pos
+
+ r_i = simulation.ice_floes[i].contact_radius
+ r_j = simulation.ice_floes[j].contact_radius
+
+ dist = norm(p)
+ dn = dist - (r_i + r_j)
+
+ if dn < 0.
+ n = p/dist
+ t = [-n[2], n[1]]
+
+ vel_lin = simulation.ice_floes[i].lin_vel -
+ simulation.ice_floes[j].lin_vel
+
+ vel_n = dot(vel_lin, n)
+ vel_t = dot(vel_lin, t) -
+ harmonicMean(r_i, r_j)*(simulation.ice_floes[i].ang_vel +
+ simulation.ice_floes[j].ang_vel)
+
+ #force_n = -kn * dn - nu * vn;
+
+ gamma_t_harmonic_mean = harmonicMean(
+ simulation.ice_floes[i].contact_viscosity_tangential,
+ simulation.ice_floes[j].contact_viscosity_tangential)
+
+ force_t = abs(gamma_t_harmonic_mean * vel_t)
+
+ mu_d_minimum = min(simulation.ice_floes[i].contact_dynamic_friction,
+ simulation.ice_floes[j].contact_dynamic_friction)
+
+ if force_t > mu_d_minimum*norm(force_n)
+ force_t = mu_d_minimum*norm(force_n)
+ end
+ if vel_t > 0.
+ force_t = -force_t
+ end
+
+ return force_t*t
+
+ end
- return -min(norm(gamma_t_harmonic_mean*contact_parallel_velocity),
- mu_d_minimum*norm(force_n))*
- contact_parallel_velocity/norm(contact_parallel_velocity)
end
function harmonicMean(a::Any, b::Any)
@@ -129,24 +191,39 @@ function harmonicMean(a::Any, b::Any)
return hm
end
-function findTorque(simulation::Simulation, overlap_vector::vector,
- force_t::vector, i::Int, j::Int)
- n = overlap_vector/norm(overlap_vector)
- return -findEffectiveRadius(simulation, i, j, overlap_vector)*
- (n[1]*force_t[2] - n[2]*force_t[1])
+function findTorque(simulation::Simulation, overlap::vector, force_t::vector,
+ i::Int, j::Int)
+ return -findEffectiveRadius(simulation, i, j, overlap)*norm(force_t)
end
function findEffectiveRadius(simulation::Simulation, i::Int, j::Int,
- overlap_vector::vector)
+ overlap::vector)
return harmonicMean(simulation.ice_floes[i].contact_radius,
- simulation.ice_floes[j].contact_radius)
- - norm(overlap_vector)/2.
+ simulation.ice_floes[j].contact_radius) -
+ norm(overlap)/2.
+end
+
+function findContactNormalVelocity(simulation::Simulation, i::Int, j::Int,
+ overlap::vector)
+ n = contactNormalVector(overlap)
+ v_ij = simulation.ice_floes[i].lin_vel - simulation.ice_floes[j].lin_vel
+ return v_ij[1]*n[1] + v_ij[2]*n[2]
end
function findContactParallelVelocity(simulation::Simulation, i::Int, j::Int,
- overlap_vector::vector)
- return simulation.ice_floes[i].lin_vel -
- simulation.ice_floes[j].lin_vel +
- findEffectiveRadius(simulation, i, j, overlap_vector)*
- (simulation.ice_floes[i].ang_vel + simulation.ice_floes[j].ang_vel)
+ overlap::vector)
+ v_ij = simulation.ice_floes[i].lin_vel - simulation.ice_floes[j].lin_vel
+ n = contactNormalVector(overlap)
+ t = contactParallelVector(n)
+ return (v_ij[1]*t[1] + v_ij[2]*t[2] -
+ findEffectiveRadius(simulation, i, j, overlap)*
+ (simulation.ice_floes[i].ang_vel + simulation.ice_floes[j].ang_vel))
+end
+
+function contactNormalVector(overlap::vector)
+ return overlap/norm(overlap)
+end
+
+function contactParallelVector(n::vector)
+ return [-n[2], n[1]]
end
diff --git a/test/collision-2floes-oblique.jl b/test/collision-2floes-oblique.jl
@@ -318,7 +318,73 @@ 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.04
+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.
+@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("# Ice floes free to move, mirrored")
+
+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_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_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)")
+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.
+@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.05
+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)
+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.04
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",
@@ -331,3 +397,69 @@ 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("# Ice floes free to move, mirrored #2")
+
+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_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_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)")
+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.
+@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.05
+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)
+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.04
+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.
+@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
diff --git a/test/contact-search-and-geometry.jl b/test/contact-search-and-geometry.jl
@@ -12,7 +12,7 @@ SeaIce.addIceFloeCylindrical(sim, [18., 0.], 10., 1., verbose=false)
position_ij = SeaIce.interIceFloePositionVector(sim, 1, 2)
overlap_ij = SeaIce.findOverlap(sim, 1, 2, position_ij)
-@test_approx_eq [18., 0.] position_ij
+@test_approx_eq [-18., 0.] position_ij
@test_approx_eq -2. overlap_ij
@@ -23,7 +23,7 @@ SeaIce.findContactsAllToAll!(sim)
@test 1 == length(sim.overlaps)
@test 1 == length(sim.contact_pairs)
@test_approx_eq [1, 2] sim.contact_pairs[1]
-@test_approx_eq [-2., 0.] sim.overlaps[1]
+@test_approx_eq [2., 0.] sim.overlaps[1]
info("Testing findContacts(...)")
@@ -34,7 +34,7 @@ sim.ice_floes[1].fixed = true
@test 1 == length(sim.overlaps)
@test 1 == length(sim.contact_pairs)
@test_approx_eq [1, 2] sim.contact_pairs[1]
-@test_approx_eq [-2., 0.] sim.overlaps[1]
+@test_approx_eq [2., 0.] sim.overlaps[1]
info("Testing findContacts(...)")
sim = deepcopy(sim_copy)
@@ -43,7 +43,7 @@ SeaIce.findContacts!(sim)
@test 1 == length(sim.overlaps)
@test 1 == length(sim.contact_pairs)
@test_approx_eq [1, 2] sim.contact_pairs[1]
-@test_approx_eq [-2., 0.] sim.overlaps[1]
+@test_approx_eq [2., 0.] sim.overlaps[1]
@test_throws ErrorException SeaIce.findContacts!(sim, method="")
@@ -76,7 +76,7 @@ SeaIce.findContacts!(sim)
@test 1 == length(sim.overlaps)
@test 1 == length(sim.contact_pairs)
@test_approx_eq [1, 2] sim.contact_pairs[1]
-@test_approx_eq [-2., 0.] sim.overlaps[1]
+@test_approx_eq [2., 0.] sim.overlaps[1]
info("Testing findContactsOceanGrid(...)")
sim = deepcopy(sim_copy)
@@ -87,7 +87,7 @@ SeaIce.findContactsOceanGrid!(sim)
@test 1 == length(sim.overlaps)
@test 1 == length(sim.contact_pairs)
@test_approx_eq [1, 2] sim.contact_pairs[1]
-@test_approx_eq [-2., 0.] sim.overlaps[1]
+@test_approx_eq [2., 0.] sim.overlaps[1]
sim = deepcopy(sim_copy)
sim.ocean = SeaIce.createRegularOceanGrid([4, 4, 2], [80., 80., 2.])
@@ -98,7 +98,7 @@ SeaIce.findContactsOceanGrid!(sim)
@test 1 == length(sim.overlaps)
@test 1 == length(sim.contact_pairs)
@test_approx_eq [1, 2] sim.contact_pairs[1]
-@test_approx_eq [-2., 0.] sim.overlaps[1]
+@test_approx_eq [2., 0.] sim.overlaps[1]
sim = deepcopy(sim_copy)
sim.ocean = SeaIce.createRegularOceanGrid([4, 4, 2], [80., 80., 2.])
@@ -119,11 +119,11 @@ SeaIce.findContacts!(sim)
@test 1 == length(sim.overlaps)
@test 1 == length(sim.contact_pairs)
@test_approx_eq [1, 2] sim.contact_pairs[1]
-@test_approx_eq [-2., 0.] sim.overlaps[1]
+@test_approx_eq [2., 0.] sim.overlaps[1]
@test 1 == length(sim.overlaps)
@test 1 == length(sim.contact_pairs)
@test_approx_eq [1, 2] sim.contact_pairs[1]
-@test_approx_eq [-2., 0.] sim.overlaps[1]
+@test_approx_eq [2., 0.] sim.overlaps[1]
@test_throws ErrorException SeaIce.findContacts!(sim, method="")
