Granular.jl

commit bb73c093ae2450a3c1b80889a5edb356157fe705
parent 5951464686ab69926cfb9d4905f42ea898d38789
Author: Anders Damsgaard <andersd@riseup.net>
Date:   Wed, 15 Nov 2017 16:25:50 -0500

finish implementing dynamic walls

Diffstat:
M
src/interaction.jl
 | 
9
+++++++--
M
src/temporal_integration.jl
 | 
3
+--
M
test/runtests.jl
 | 
2
+-
M
test/wall.jl
 | 
53
++++++++++++++++++++++++++---------------------------

4 files changed, 35 insertions(+), 32 deletions(-)
diff --git a/src/interaction.jl b/src/interaction.jl
@@ -46,12 +46,16 @@ interaction is frictionless in the tangential direction.
 """
 function interactWalls!(sim::Simulation)
 
+    orientation::Float64 = 0.0
     δ_n::Float64 = 0.0
     k_n::Float64 = 0.0
 
     for iw=1:length(sim.walls)
         for i=1:length(sim.grains)
 
+            orientation = sign(dot(sim.walls[iw].normal,
+                                   sim.grains[i].lin_pos) - sim.walls[iw].pos)
+
             # get overlap distance by projecting grain position onto wall-normal
             # vector
             δ_n = abs(dot(sim.walls[iw].normal, sim.grains[i].lin_pos) -
@@ -64,8 +68,9 @@ function interactWalls!(sim::Simulation)
                     k_n = sim.grains[i].contact_stiffness_normal
                 end
 
-                sim.grains[i].force += k_n * abs(δ_n) * sim.walls[iw].normal
-                sim.walls[iw].force += k_n * δ_n
+                sim.grains[i].force += -k_n * δ_n .* sim.walls[iw].normal .*
+                    orientation
+                sim.walls[iw].force += k_n * δ_n * orientation
             end
         end
     end
diff --git a/src/temporal_integration.jl b/src/temporal_integration.jl
@@ -220,8 +220,7 @@ function updateWallKinematicsThreeTermTaylor!(wall::WallLinearFrictionless,
     if wall.bc == "velocity"
         wall.acc = 0.0
     else
-        f_n::Float64 = -wall.normal_stress*wall.surface_area
-        wall.acc = (wall.force + f_n)/wall.mass
+        wall.acc = (wall.force + wall.normal_stress*wall.surface_area)/wall.mass
     end
 
     # Temporal gradient in acceleration using backwards differences
diff --git a/test/runtests.jl b/test/runtests.jl
@@ -1,8 +1,8 @@
 using Compat.Test
 import Granular
 
-include("grain.jl")
 include("wall.jl")
+include("grain.jl")
 include("packing.jl")
 include("util.jl")
 include("temporal.jl")
diff --git a/test/wall.jl b/test/wall.jl
@@ -57,7 +57,7 @@ sim = Granular.createSimulation()
 sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
 Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
 Granular.addWallLinearFrictionless!(sim, [1., 0.], -1.01, verbose=false)
-Granular.setTimeStep!(sim)
+Granular.setTimeStep!(sim, verbose=false)
 Granular.interactWalls!(sim)
 @test sim.walls[1].force ≈ 0.
 @test sim.grains[1].force[1] ≈ 0.
@@ -68,7 +68,7 @@ sim = Granular.createSimulation()
 sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
 Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
 Granular.addWallLinearFrictionless!(sim, [1., 0.], +2.01, verbose=false)
-Granular.setTimeStep!(sim)
+Granular.setTimeStep!(sim, verbose=false)
 Granular.interactWalls!(sim)
 @test sim.walls[1].force ≈ 0.
 @test sim.grains[1].force[1] ≈ 0.
@@ -79,10 +79,10 @@ sim = Granular.createSimulation()
 sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
 Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
 Granular.addWallLinearFrictionless!(sim, [1., 0.], -1. + .01, verbose=false)
-Granular.setTimeStep!(sim)
+Granular.setTimeStep!(sim, verbose=false)
 Granular.interactWalls!(sim)
-@test sim.walls[1].force > 0.
-@test sim.grains[1].force[1] < 0.
+@test sim.walls[1].force < 0.
+@test sim.grains[1].force[1] > 0.
 @test sim.grains[1].force[2] ≈ 0.
 
 info("Wall at +x")
@@ -90,7 +90,7 @@ sim = Granular.createSimulation()
 sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
 Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
 Granular.addWallLinearFrictionless!(sim, [1., 0.], +1. - .01, verbose=false)
-Granular.setTimeStep!(sim)
+Granular.setTimeStep!(sim, verbose=false)
 Granular.interactWalls!(sim)
 @test sim.walls[1].force > 0.
 @test sim.grains[1].force[1] < 0.
@@ -101,7 +101,7 @@ sim = Granular.createSimulation()
 sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
 Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
 Granular.addWallLinearFrictionless!(sim, [0., 1.], -1. + .01, verbose=false)
-Granular.setTimeStep!(sim)
+Granular.setTimeStep!(sim, verbose=false)
 Granular.interactWalls!(sim)
 @test sim.walls[1].force < 0.
 @test sim.grains[1].force[1] ≈ 0.
@@ -112,7 +112,7 @@ sim = Granular.createSimulation()
 sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
 Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
 Granular.addWallLinearFrictionless!(sim, [0., 1.], +1. - .01, verbose=false)
-Granular.setTimeStep!(sim)
+Granular.setTimeStep!(sim, verbose=false)
 Granular.interactWalls!(sim)
 @test sim.walls[1].force > 0.
 @test sim.grains[1].force[1] ≈ 0.
@@ -126,7 +126,7 @@ sim = Granular.createSimulation()
 sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
 Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
 Granular.addWallLinearFrictionless!(sim, [1., 0.], -1.01, verbose=false)
-Granular.setTimeStep!(sim)
+Granular.setTimeStep!(sim, verbose=false)
 Granular.interactWalls!(sim)
 Granular.updateWallKinematics!(sim)
 @test sim.walls[1].force ≈ 0.
@@ -141,7 +141,7 @@ sim = Granular.createSimulation()
 sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
 Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
 Granular.addWallLinearFrictionless!(sim, [1., 0.], -1.01, verbose=false)
-Granular.setTimeStep!(sim)
+Granular.setTimeStep!(sim, verbose=false)
 Granular.interactWalls!(sim)
 Granular.updateWallKinematics!(sim, method="Two-term Taylor")
 @test sim.walls[1].force ≈ 0.
@@ -156,7 +156,7 @@ sim = Granular.createSimulation()
 sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
 Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
 Granular.addWallLinearFrictionless!(sim, [1., 0.], -1.01, verbose=false)
-Granular.setTimeStep!(sim)
+Granular.setTimeStep!(sim, verbose=false)
 Granular.interactWalls!(sim)
 @test_throws ErrorException Granular.updateWallKinematics!(sim, method="asdf")
 Granular.updateWallKinematics!(sim, method="Three-term Taylor")
@@ -172,7 +172,7 @@ sim = Granular.createSimulation()
 sim.ocean = Granular.createRegularOceanGrid([1, 1, 1], [10., 20., 1.0])
 Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
 Granular.addWallLinearFrictionless!(sim, [1., 0.], -1.01, verbose=false)
-Granular.setTimeStep!(sim)
+Granular.setTimeStep!(sim, verbose=false)
 Granular.interactWalls!(sim)
 Granular.updateWallKinematics!(sim)
 @test sim.walls[1].acc ≈ 0.
@@ -186,7 +186,7 @@ Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
 Granular.addWallLinearFrictionless!(sim, [1., 0.], -1.01,
                                     bc="velocity", vel=1.0,
                                     verbose=false)
-Granular.setTimeStep!(sim)
+Granular.setTimeStep!(sim, verbose=false)
 Granular.interactWalls!(sim)
 Granular.updateWallKinematics!(sim)
 @test sim.walls[1].force ≈ 0.
@@ -203,7 +203,7 @@ Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
 Granular.addWallLinearFrictionless!(sim, [1., 0.], -1.01,
                                     bc="velocity", vel=1.0,
                                     verbose=false)
-Granular.setTimeStep!(sim)
+Granular.setTimeStep!(sim, verbose=false)
 Granular.interactWalls!(sim)
 Granular.updateWallKinematics!(sim, method="Two-term Taylor")
 @test sim.walls[1].force ≈ 0.
@@ -220,7 +220,7 @@ Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
 Granular.addWallLinearFrictionless!(sim, [1., 0.], -1.01,
                                     bc="velocity", vel=1.0,
                                     verbose=false)
-Granular.setTimeStep!(sim)
+Granular.setTimeStep!(sim, verbose=false)
 Granular.interactWalls!(sim)
 @test_throws ErrorException Granular.updateWallKinematics!(sim, method="asdf")
 Granular.updateWallKinematics!(sim, method="Three-term Taylor")
@@ -238,7 +238,7 @@ Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
 Granular.addWallLinearFrictionless!(sim, [1., 0.], -0.9,
                                     bc="velocity", vel=1.0,
                                     verbose=false)
-Granular.setTimeStep!(sim)
+Granular.setTimeStep!(sim, verbose=false)
 Granular.interactWalls!(sim)
 Granular.updateWallKinematics!(sim, method="Two-term Taylor")
 @test sim.walls[1].bc == "velocity"
@@ -253,7 +253,7 @@ Granular.addGrainCylindrical!(sim, [ 0., 0.], 1., 2., verbose=false)
 Granular.addWallLinearFrictionless!(sim, [1., 0.], -0.9,
                                     bc="velocity", vel=1.0,
                                     verbose=false)
-Granular.setTimeStep!(sim)
+Granular.setTimeStep!(sim, verbose=false)
 Granular.interactWalls!(sim)
 Granular.updateWallKinematics!(sim, method="Two-term Taylor")
 @test sim.walls[1].acc ≈ 0.
@@ -268,7 +268,7 @@ Granular.addWallLinearFrictionless!(sim, [1., 0.], 2.,
                                     bc="normal stress",
                                     normal_stress=0.,
                                     verbose=false)
-Granular.setTimeStep!(sim)
+Granular.setTimeStep!(sim, verbose=false)
 Granular.interactWalls!(sim)
 Granular.updateWallKinematics!(sim)
 @test sim.walls[1].force ≈ 0.
@@ -284,9 +284,9 @@ Granular.addGrainCylindrical!(sim, [ 1., 1.], 1., 1., verbose=false)
 sim.grains[1].fixed = true
 Granular.addWallLinearFrictionless!(sim, [1., 0.], 2.,
                                     bc="normal stress",
-                                    normal_stress=1e4,
+                                    normal_stress=-1e4,
                                     verbose=false)
-Granular.setTimeStep!(sim)
+Granular.setTimeStep!(sim, verbose=false)
 Granular.interactWalls!(sim)
 Granular.updateWallKinematics!(sim)
 @test sim.walls[1].force ≈ 0.
@@ -298,11 +298,10 @@ for i=1:5
     Granular.interactWalls!(sim)
     Granular.updateWallKinematics!(sim)
     println(sim.walls[1].pos)
+    @test sim.walls[1].force > 0.
+    @test sim.walls[1].acc < 0.
+    @test sim.walls[1].vel < 0.
+    @test sim.walls[1].pos < 2.
+    @test sim.grains[1].force[1] < 0.
+    @test sim.grains[1].force[2] ≈ 0.
 end
-@test sim.walls[1].force < 0.
-@test sim.walls[1].acc < 0.
-@test sim.walls[1].vel < 0.
-@test sim.walls[1].pos < 1.
-@test sim.grains[1].force[1] > 0.
-@test sim.grains[1].force[2] ≈ 0.
-