tfix interaction scheme - 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 7515db2c7960dea61383549159b3b94308a9fc75
(DIR) parent dadd1cadc476c6f083ecf13f38110e453b465257
(HTM) Author: Anders Damsgaard <andersd@riseup.net>
Date: Wed, 10 May 2017 15:33:57 -0400
fix interaction scheme
Diffstat:
M src/contact_search.jl | 25 +++++++++++++------------
M src/icefloe.jl | 45 +++++++++++++++++++++++++++++++
M src/interaction.jl | 18 ++++++++----------
M test/collision-5floes-normal.jl | 47 ++++++++++++++++++++++++++++++-
4 files changed, 112 insertions(+), 23 deletions(-)
---
(DIR) diff --git a/src/contact_search.jl b/src/contact_search.jl
t@@ -124,8 +124,7 @@ function checkAndAddContact!(sim::Simulation, i::Int, j::Int)
if i < j
if (sim.ice_floes[i].fixed && sim.ice_floes[j].fixed) ||
- !sim.ice_floes[i].enabled ||
- !sim.ice_floes[j].enabled
+ !sim.ice_floes[i].enabled || !sim.ice_floes[j].enabled
return
end
t@@ -140,16 +139,18 @@ function checkAndAddContact!(sim::Simulation, i::Int, j::Int)
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
+ 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
+ sim.ice_floes[i].contact_parallel_displacement[ic] =
+ zeros(2)
+ break
+ end
end
end
end
(DIR) diff --git a/src/icefloe.jl b/src/icefloe.jl
t@@ -261,6 +261,51 @@ function convertIceFloeDataToArrays(simulation::Simulation)
return ifarr
end
+export printIceFloeInfo
+"""
+ printIceFloeInfo(icefloe::IceFloeCylindrical)
+
+Prints the contents of an ice floe to stdout in a formatted style.
+"""
+function printIceFloeInfo(f::IceFloeCylindrical)
+ println(" density: $(f.density) kg/m^3")
+ println(" thickness: $(f.thickness) m")
+ println(" contact_radius: $(f.contact_radius) m")
+ println(" areal_radius: $(f.areal_radius) m")
+ println(" circumreference: $(f.circumreference) m")
+ println(" horizontal_surface_area: $(f.horizontal_surface_area) m^2")
+ println(" side_surface_area: $(f.side_surface_area) m^2")
+ println(" volume: $(f.volume) m^3")
+ println(" mass: $(f.mass) kg")
+ println(" moment_of_inertia: $(f.moment_of_inertia) kg*m^2\n")
+
+ println(" lin_pos: $(f.lin_pos) m")
+ println(" lin_vel: $(f.lin_vel) m/s")
+ println(" lin_acc: $(f.lin_acc) m/s^2")
+ println(" force: $(f.force) N\n")
+
+ println(" ang_pos: $(f.ang_pos) rad")
+ println(" ang_vel: $(f.ang_vel) rad/s")
+ println(" ang_acc: $(f.ang_acc) rad/s^2")
+ println(" torque: $(f.torque) N*m\n")
+
+ println(" fixed: $(f.fixed)")
+ println(" rotating: $(f.rotating)")
+ println(" enabled: $(f.enabled)\n")
+
+ println(" k_n: $(f.contact_stiffness_normal) N/m")
+ println(" k_t: $(f.contact_stiffness_tangential) N/m")
+ println(" gamma_n: $(f.contact_viscosity_normal) N/(m/s)")
+ println(" gamma_t: $(f.contact_viscosity_tangential) N/(m/s)")
+ println(" mu_s: $(f.contact_static_friction)")
+ println(" mu_d: $(f.contact_dynamic_friction)\n")
+
+ println(" pressure: $(f.pressure) Pa")
+ println(" n_contacts: $(f.n_contacts)")
+ println(" contacts: $(f.contacts)")
+ println(" delta_t: $(f.contact_parallel_displacement)")
+end
+
export iceFloeKineticTranslationalEnergy
"Returns the translational kinetic energy of the ice floe"
function iceFloeKineticTranslationalEnergy(icefloe::IceFloeCylindrical)
(DIR) diff --git a/src/interaction.jl b/src/interaction.jl
t@@ -7,13 +7,13 @@ export interact!
Resolve mechanical interaction between all particle pairs.
"""
function interact!(simulation::Simulation)
- for i=1:Int(ceil(length(simulation.ice_floes)/2.)) # i <= Int(N/2)
+ for i=1:length(simulation.ice_floes)
for ic=1:Nc_max
j = simulation.ice_floes[i].contacts[ic]
- if j == 0
- break # end of contact list reached
+ if i > j # skip i > j and j == 0
+ continue
end
if norm(simulation.ice_floes[i].lin_pos -
t@@ -52,7 +52,10 @@ function interactIceFloes!(simulation::Simulation, i::Int, j::Int, ic::Int)
# Floe distance
delta_n = dist - (r_i + r_j)
- if delta_n < 0. # Contact (this should always occur)
+ if delta_n > 0. # Double-check contact
+ error("function called to process non-existent contact between ice " *
+ "floes $i and $j")
+ else
# Local axes
n = p/dist
t@@ -128,8 +131,7 @@ function interactIceFloes!(simulation::Simulation, i::Int, j::Int, ic::Int)
elseif k_t_harmonic_mean > 0.
- force_t = -k_t_harmonic_mean*delta_t -
- gamma_t_harmonic_mean*vel_t
+ 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)
t@@ -141,10 +143,6 @@ function interactIceFloes!(simulation::Simulation, i::Int, j::Int, ic::Int)
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
simulation.ice_floes[i].contact_parallel_displacement[ic] = delta_t*t
(DIR) diff --git a/test/collision-5floes-normal.jl b/test/collision-5floes-normal.jl
t@@ -39,6 +39,11 @@ 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
+@test 0. < norm(sim.ice_floes[1].lin_vel)
+for i=2:5
+ info("testing ice floe $i")
+ @test 0. ≈ norm(sim.ice_floes[i].lin_vel)
+end
info("Testing kinetic energy conservation with Two-term Taylor scheme")
t@@ -52,6 +57,11 @@ 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
+@test 0. < norm(sim.ice_floes[1].lin_vel)
+for i=2:5
+ info("testing ice floe $i")
+ @test 0. ≈ norm(sim.ice_floes[i].lin_vel)
+end
info("Testing kinetic energy conservation with Three-term Taylor scheme")
t@@ -66,6 +76,11 @@ 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
+@test 0. < norm(sim.ice_floes[1].lin_vel)
+for i=2:5
+ info("testing ice floe $i")
+ @test 0. ≈ norm(sim.ice_floes[i].lin_vel)
+end
info("# Ice floes free to move")
t@@ -84,7 +99,7 @@ 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)
+SeaIce.setTotalTime!(sim, 40.0)
sim_init = deepcopy(sim)
info("Testing kinetic energy conservation with Two-term Taylor scheme")
t@@ -97,6 +112,10 @@ 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
+for i=1:5
+ info("testing ice floe $i")
+ @test 0. < norm(sim.ice_floes[i].lin_vel)
+end
info("Testing kinetic energy conservation with Two-term Taylor scheme")
t@@ -110,6 +129,10 @@ 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
+for i=1:5
+ info("testing ice floe $i")
+ @test 0. < norm(sim.ice_floes[i].lin_vel)
+end
info("Testing kinetic energy conservation with Three-term Taylor scheme")
t@@ -124,6 +147,10 @@ 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
+for i=1:5
+ info("testing ice floe $i")
+ @test 0. < norm(sim.ice_floes[i].lin_vel)
+end
info("# Adding contact-normal viscosity")
t@@ -166,6 +193,11 @@ 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
+@test 0. < norm(sim.ice_floes[1].lin_vel)
+for i=2:5
+ info("testing ice floe $i")
+ @test 0. ≈ norm(sim.ice_floes[i].lin_vel)
+end
info("Testing kinetic energy conservation with Three-term Taylor scheme")
t@@ -180,6 +212,11 @@ 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
+@test 0. < norm(sim.ice_floes[1].lin_vel)
+for i=2:5
+ info("testing ice floe $i")
+ @test 0. ≈ norm(sim.ice_floes[i].lin_vel)
+end
info("# Ice floes free to move")
t@@ -217,6 +254,10 @@ 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
+for i=1:5
+ info("testing ice floe $i")
+ @test 0. < norm(sim.ice_floes[i].lin_vel)
+end
info("Testing kinetic energy conservation with Three-term Taylor scheme")
t@@ -231,3 +272,7 @@ 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
+for i=1:5
+ info("testing ice floe $i")
+ @test 0. < norm(sim.ice_floes[i].lin_vel)
+end