commit 23fc4a029c9f5c4a8b4781eb1404f4049fda3add
parent 0440db5e1befd66a7f4c9fb98061d8f4604d56f1
Author: Anders Damsgaard <andersd@riseup.net>
Date: Mon, 8 May 2017 21:18:02 -0400
concatenate interaction function, fix rotation issues
Diffstat:
2 files changed, 91 insertions(+), 179 deletions(-)
diff --git a/src/interaction.jl b/src/interaction.jl
@@ -35,152 +35,101 @@ function interactIceFloes!(simulation::Simulation,
contact_normal_rheology::String = "Linear Elastic",
contact_tangential_rheology::String = "None")
- force_n = zeros(2)
- force_t = zeros(2)
-
- if contact_normal_rheology == "None"
- # do nothing
- elseif contact_normal_rheology == "Linear Elastic"
- force_n = interactNormalLinearElastic(simulation, i, j, overlap)
- else
- error("Unknown contact_normal_rheology '$contact_normal_rheology'")
- end
-
- if contact_tangential_rheology == "None"
- # do nothing
- elseif contact_tangential_rheology == "Linear Viscous Frictional"
- force_t = interactTangentialLinearViscousFrictional(simulation, i, j,
- overlap,
- force_n)
- else
- error("Unknown contact_tangential_rheology ",
- "'$contact_tangential_rheology'")
- end
-
- simulation.ice_floes[i].force += force_n + force_t;
- simulation.ice_floes[j].force -= force_n + force_t;
-
- if norm(force_t) > 0.
- torque = -findTorque(simulation, overlap, force_t, i, j)
- simulation.ice_floes[i].torque += torque
- simulation.ice_floes[j].torque += torque
- end
-
- simulation.ice_floes[i].pressure +=
- norm(force_n)/simulation.ice_floes[i].side_surface_area;
- simulation.ice_floes[j].pressure +=
- norm(force_n)/simulation.ice_floes[j].side_surface_area;
-end
-
-export interactNormalLinearElastic
-"""
-Resolves linear-elastic interaction between two ice floes in the contact-normal
-direction.
-"""
-function interactNormalLinearElastic(simulation::Simulation,
- i::Int, j::Int,
- 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
-end
-
-export interactTangentialLinearElastic
-"""
- 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,
- force_n::vector)
-
- """
- contact_parallel_velocity = findContactParallelVelocity(simulation, i, j,
- overlap)
-
-
-
- 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 gamma_t_harmonic_mean ≈ 0. || mu_d_minimum ≈ 0.
- return [0., 0.]
- end
-
- 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))
- """
+ force_n = 0.
+ force_t = 0.
p = simulation.ice_floes[i].lin_pos - simulation.ice_floes[j].lin_pos
+ dist = norm(p)
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.
+ if dn < 0. # Contact (this should always occur)
+
+ # Local axes
n = p/dist
t = [-n[2], n[1]]
+ # Contact kinematics
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;
+ R_ij = harmonicMean(simulation.ice_floes[i].contact_radius,
+ simulation.ice_floes[j].contact_radius) -
+ norm(overlap)/2.
+
+ # Contact mechanical parameters
+ k_n_harmonic_mean =
+ harmonicMean(simulation.ice_floes[i].contact_stiffness_normal,
+ simulation.ice_floes[j].contact_stiffness_normal)
+
+ k_t_harmonic_mean =
+ harmonicMean(simulation.ice_floes[i].contact_stiffness_tangential,
+ simulation.ice_floes[j].contact_stiffness_tangential)
+
+ gamma_n_harmonic_mean = harmonicMean(
+ simulation.ice_floes[i].contact_viscosity_normal,
+ simulation.ice_floes[j].contact_viscosity_normal)
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
+ # Determine contact forces
+ if contact_normal_rheology == "None"
+ force_n = 0.
+ elseif contact_normal_rheology == "Linear Elastic"
+ force_n = -k_n_harmonic_mean*dn
+ elseif contact_normal_rheology == "Linear Viscous Elastic"
+ force_n = -k_n_harmonic_mean*dn - gamma_n_harmonic_mean*vel_n
+ if force_n < 0.
+ force_n = 0.
+ end
+ else
+ error("unknown contact_normal_rheology '$contact_normal_rheology'")
end
- return force_t*t
+ if contact_tangential_rheology == "None"
+ # do nothing
+ elseif contact_tangential_rheology == "Linear Elastic Frictional"
+ error("not yet implemented")
+ elseif contact_tangential_rheology == "Linear Viscous Frictional"
+ force_t = abs(gamma_t_harmonic_mean * vel_t)
+ 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
+ else
+ error("unknown contact_tangential_rheology ",
+ "'$contact_tangential_rheology'")
+ end
+ else
+ error("function called to process non-existent contact between ice " *
+ "floes $i and $j")
end
+ simulation.ice_floes[i].force += force_n*n + force_t*t;
+ simulation.ice_floes[j].force -= force_n*n + force_t*t;
+
+ simulation.ice_floes[i].torque += -force_t*R_ij
+ simulation.ice_floes[j].torque += -force_t*R_ij
+
+ simulation.ice_floes[i].pressure +=
+ 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)
@@ -190,40 +139,3 @@ function harmonicMean(a::Any, b::Any)
end
return hm
end
-
-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)
- return harmonicMean(simulation.ice_floes[i].contact_radius,
- 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)
- 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
@@ -27,7 +27,7 @@ sim_init = deepcopy(sim)
info("Testing kinetic energy conservation with Two-term Taylor scheme")
SeaIce.setTimeStep!(sim, epsilon=0.07)
-tol = 0.2
+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)
@@ -41,7 +41,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.02
+tol = 0.01
info("Relative tolerance: $(tol*100.)%")
SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
contact_tangential_rheology="None", verbose=verbose)
@@ -83,7 +83,7 @@ sim_init = deepcopy(sim)
info("Testing kinetic energy conservation with Two-term Taylor scheme")
SeaIce.setTimeStep!(sim, epsilon=0.07)
-tol = 0.2
+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)
@@ -97,7 +97,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.02
+tol = 0.01
info("Relative tolerance: $(tol*100.)%")
SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
contact_tangential_rheology="None", verbose=verbose)
@@ -140,7 +140,7 @@ sim_init = deepcopy(sim)
info("Testing kinetic energy conservation with Two-term Taylor scheme")
SeaIce.setTimeStep!(sim, epsilon=0.07)
-tol = 0.2
+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)
@@ -159,7 +159,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.02
+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)
@@ -198,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.2
+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",
@@ -216,7 +216,7 @@ 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
+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)
@@ -235,7 +235,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.06
+tol = 0.02
info("Relative tolerance: $(tol*100.)%")
SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
contact_tangential_rheology="Linear Viscous Frictional",
@@ -253,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.07
+tol = 0.03
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",
@@ -287,7 +287,7 @@ sim_init = deepcopy(sim)
info("Testing kinetic energy conservation with Two-term Taylor scheme")
SeaIce.setTimeStep!(sim, epsilon=0.07)
-tol = 0.2
+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",
@@ -304,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.05
+tol = 0.02
info("Relative tolerance: $(tol*100.)%")
SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
contact_tangential_rheology="Linear Viscous Frictional",
@@ -318,16 +318,16 @@ 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
+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.
-@test sim.ice_floes[1].ang_vel > 0.
-@test sim.ice_floes[2].ang_pos > 0.
-@test sim.ice_floes[2].ang_vel > 0.
+@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
@@ -353,7 +353,7 @@ sim_init = deepcopy(sim)
info("Testing kinetic energy conservation with Two-term Taylor scheme")
SeaIce.setTimeStep!(sim, epsilon=0.07)
-tol = 0.2
+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",
@@ -370,7 +370,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.05
+tol = 0.02
info("Relative tolerance: $(tol*100.)%")
SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
contact_tangential_rheology="Linear Viscous Frictional",
@@ -384,16 +384,16 @@ 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
+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.
-@test sim.ice_floes[1].ang_vel < 0.
-@test sim.ice_floes[2].ang_pos < 0.
-@test sim.ice_floes[2].ang_vel < 0.
+@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
@@ -419,7 +419,7 @@ sim_init = deepcopy(sim)
info("Testing kinetic energy conservation with Two-term Taylor scheme")
SeaIce.setTimeStep!(sim, epsilon=0.07)
-tol = 0.2
+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",
@@ -436,7 +436,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.05
+tol = 0.02
info("Relative tolerance: $(tol*100.)%")
SeaIce.run!(sim, temporal_integration_method="Two-term Taylor",
contact_tangential_rheology="Linear Viscous Frictional",
@@ -450,7 +450,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.04
+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",
