tzero forces and torques at the beginning of each time step and expand integration tests - 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 1c2e3c2a2abc0125d658dfb267281e664bf32641
(DIR) parent fd905feff759e602d9239e0f0b64bb0a05206678
(HTM) Author: Anders Damsgaard <andersd@riseup.net>
Date: Thu, 20 Apr 2017 11:15:42 -0400
zero forces and torques at the beginning of each time step and expand integration tests
Diffstat:
M src/interaction.jl | 2 --
M src/simulation.jl | 13 +++++++++++--
M src/temporal_integration.jl | 14 ++++++++++++++
M test/collision-2floes-normal.jl | 49 ++++++++++++++++++++++++++-----
M test/contact-search-and-geometry.jl | 12 ++++++++++++
5 files changed, 79 insertions(+), 11 deletions(-)
---
(DIR) diff --git a/src/interaction.jl b/src/interaction.jl
t@@ -4,8 +4,6 @@
Resolve mechanical interaction between all grain pairs and walls.
"""
function interact!(simulation::Simulation)
- contact_pair = Array{Integer, 1}
- #overlap_ij = float
# IceFloe to grain collisions
while !isempty(simulation.contact_pairs)
(DIR) diff --git a/src/simulation.jl b/src/simulation.jl
t@@ -29,7 +29,8 @@ function run!(simulation::Simulation;
verbose::Bool = true,
status_interval = 100.,
show_file_output = true,
- single_step=false)
+ single_step=false,
+ temporal_integration_method="Three-term Taylor")
if single_step && simulation.time >= simulation.time_total
simulation.time_total += simulation.time_step
t@@ -56,9 +57,10 @@ function run!(simulation::Simulation;
" s ")
end
+ zeroForcesAndTorques!(simulation)
findContacts!(simulation)
interact!(simulation)
- updateIceFloeKinematics!(simulation)
+ updateIceFloeKinematics!(simulation, method=temporal_integration_method)
# Update time variables
simulation.time_iteration += 1
t@@ -95,3 +97,10 @@ function removeIceFloe!(simulation::Simulation, i::Integer)
delete!(simulation.ice_floes, i)
end
+
+function zeroForcesAndTorques!(simulation::Simulation)
+ for icefloe in simulation.ice_floes
+ icefloe.force = zeros(2)
+ icefloe.torque = 0.
+ end
+end
(DIR) diff --git a/src/temporal_integration.jl b/src/temporal_integration.jl
t@@ -23,6 +23,13 @@ function updateIceFloeKinematicsTwoTermTaylor(icefloe::IceFloeCylindrical,
icefloe.lin_acc = icefloe.force/icefloe.mass
icefloe.ang_acc = icefloe.torque/icefloe.moment_of_inertia
+ if icefloe.fixed
+ icefloe.lin_acc = zeros(2)
+ icefloe.ang_acc = 0.
+ elseif !icefloe.rotating
+ icefloe.ang_acc = 0.
+ end
+
icefloe.lin_pos +=
icefloe.lin_vel * simulation.time_step +
0.5*icefloe.lin_acc * simulation.time_step^2.0
t@@ -48,6 +55,13 @@ function updateIceFloeKinematicsThreeTermTaylor(icefloe::IceFloeCylindrical,
icefloe.lin_acc = icefloe.force/icefloe.mass
icefloe.ang_acc = icefloe.torque/icefloe.moment_of_inertia
+ if icefloe.fixed
+ icefloe.lin_acc = zeros(2)
+ icefloe.ang_acc = 0.
+ elseif !icefloe.rotating
+ icefloe.ang_acc = 0.
+ end
+
# Temporal gradient in acceleration using backwards differences
d_lin_acc_dt = (icefloe.lin_acc - lin_acc_0)/simulation.time_step
d_ang_acc_dt = (icefloe.ang_acc - ang_acc_0)/simulation.time_step
(DIR) diff --git a/test/collision-2floes-normal.jl b/test/collision-2floes-normal.jl
t@@ -8,20 +8,55 @@ import SeaIce
info("#### $(basename(@__FILE__)) ####")
+
sim = SeaIce.createSimulation(id="test")
-SeaIce.addIceFloeCylindrical(sim, [ 0., 0.], 10., 1., verbose=false)
-SeaIce.addIceFloeCylindrical(sim, [20., 0.], 10., 1., verbose=false)
-sim.ice_floes[1].lin_vel[1] = 1.
+SeaIce.addIceFloeCylindrical(sim, [0., 0.], 10., 1., verbose=false)
+SeaIce.addIceFloeCylindrical(sim, [20.05, 0.], 10., 1., verbose=false)
+sim.ice_floes[1].lin_vel[1] = 0.1
+sim.ice_floes[2].fixed = true
E_kin_lin_init = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
E_kin_rot_init = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
-SeaIce.setTimeStep!(sim)
-SeaIce.setTotalTime!(sim, 1.0)
-SeaIce.run!(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")
+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")
+
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+Base.Test.@test_approx_eq_eps E_kin_lin_init E_kin_lin_final E_kin_lin_init*tol
+Base.Test.@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.02
+info("Relative tolerance: $(tol*100.)%")
+SeaIce.run!(sim, temporal_integration_method="Two-term Taylor")
E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+Base.Test.@test_approx_eq_eps E_kin_lin_init E_kin_lin_final E_kin_lin_init*tol
+Base.Test.@test_approx_eq E_kin_rot_init E_kin_rot_final
+
-Base.Test.@test_approx_eq E_kin_lin_init E_kin_lin_final
+info("Testing kinetic energy conservation with Two-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")
+
+E_kin_lin_final = SeaIce.totalIceFloeKineticTranslationalEnergy(sim)
+E_kin_rot_final = SeaIce.totalIceFloeKineticRotationalEnergy(sim)
+Base.Test.@test_approx_eq_eps E_kin_lin_init E_kin_lin_final E_kin_lin_init*tol
Base.Test.@test_approx_eq E_kin_rot_init E_kin_rot_final
(DIR) diff --git a/test/contact-search-and-geometry.jl b/test/contact-search-and-geometry.jl
t@@ -36,3 +36,15 @@ Base.Test.@test 1 == length(sim.overlaps)
Base.Test.@test 1 == length(sim.contact_pairs)
Base.Test.@test_approx_eq [1, 2] sim.contact_pairs[1]
Base.Test.@test_approx_eq [-2., 0.] sim.overlaps[1]
+
+
+info("Testing if interact(...) removes contacts correctly")
+sim = deepcopy(sim_copy)
+SeaIce.findContacts!(sim)
+SeaIce.interact!(sim)
+SeaIce.findContacts!(sim)
+
+Base.Test.@test 1 == length(sim.overlaps)
+Base.Test.@test 1 == length(sim.contact_pairs)
+Base.Test.@test_approx_eq [1, 2] sim.contact_pairs[1]
+Base.Test.@test_approx_eq [-2., 0.] sim.overlaps[1]