tadd tests for contact model based on Young's modulus - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
(HTM) git clone git://src.adamsgaard.dk/sphere
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
(DIR) commit f67134dcadd498588f5689329190fa3940efff8d
(DIR) parent 29ce3482c183769fdafcfbadf4a957831bee4e04
(HTM) Author: Anders Damsgaard Christensen <adc@geo.au.dk>
Date: Tue, 16 Aug 2016 12:12:44 -0700
add tests for contact model based on Young's modulus
Diffstat:
M tests/CMakeLists.txt | 6 ++++++
A tests/contactmodel_wall_young.py | 115 +++++++++++++++++++++++++++++++
A tests/contactmodel_young.py | 232 ++++++++++++++++++++++++++++++
3 files changed, 353 insertions(+), 0 deletions(-)
---
(DIR) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
t@@ -9,6 +9,12 @@ add_test(wall_contact_model_tests ${PYTHON_EXECUTABLE}
add_test(contact_model_tests ${PYTHON_EXECUTABLE}
${CMAKE_CURRENT_BINARY_DIR}/contactmodel.py)
+add_test(wall_contact_model_young_tests ${PYTHON_EXECUTABLE}
+ ${CMAKE_CURRENT_BINARY_DIR}/contactmodel_wall_young.py)
+
+add_test(contact_model_young_tests ${PYTHON_EXECUTABLE}
+ ${CMAKE_CURRENT_BINARY_DIR}/contactmodel_young.py)
+
add_test(io_tests_fluid ${PYTHON_EXECUTABLE}
${CMAKE_CURRENT_BINARY_DIR}/io_tests_fluid.py)
(DIR) diff --git a/tests/contactmodel_wall_young.py b/tests/contactmodel_wall_young.py
t@@ -0,0 +1,115 @@
+#!/usr/bin/env python
+'''
+Validate the implemented contact models by observing the behavior of one or two
+particles.
+'''
+
+import sphere
+import numpy
+import pytestutils
+
+### Wall-particle interaction ##################################################
+
+## Linear elastic collisions
+
+# Normal impact: Check for conservation of momentum (sum(v_i*m_i))
+orig = sphere.sim(np=1, nw=0, sid='contactmodeltest')
+sphere.cleanup(orig)
+orig.radius[:] = 1.0
+orig.x[0,:] = [5.0, 5.0, 1.05]
+orig.setYoungsModulus(70.0e9)
+orig.vel[0,2] = -0.1
+orig.defineWorldBoundaries(L=[10,10,10])
+orig.gamma_wn[0] = 0.0 # Disable wall viscosity
+orig.gamma_wt[0] = 0.0 # Disable wall viscosity
+orig.initTemporal(total = 1.0, file_dt = 0.01)
+#orig.time_dt = orig.time_dt*0.1
+moment_before = orig.totalKineticEnergy()
+orig.run(verbose=False)
+#orig.writeVTKall()
+orig.readlast(verbose=False)
+pytestutils.compareFloats(orig.vel[0,2], 0.1,\
+ "Elastic normal wall collision (1/2):")
+moment_after = orig.totalKineticEnergy()
+#print(moment_before)
+#print(moment_after)
+#print("time step: " + str(orig.time_dt[0]))
+#print(str((moment_after[0]-moment_before[0])/moment_before[0]*100.0) + " %")
+pytestutils.compareFloats(moment_before, moment_after,\
+ "Elastic normal wall collision (2/2):")
+
+# Oblique impact: Check for conservation of momentum (sum(v_i*m_i))
+orig = sphere.sim(np=1, sid='contactmodeltest')
+orig.radius[:] = 1.0
+orig.x[0,:] = [5.0, 5.0, 1.05]
+orig.setYoungsModulus(70.0e9)
+orig.vel[0,2] = -0.1
+orig.vel[0,0] = 0.1
+orig.defineWorldBoundaries(L=[10,10,10])
+orig.gamma_wn[0] = 0.0 # Disable wall viscosity
+orig.gamma_wt[0] = 0.0 # Disable wall viscosity
+orig.initTemporal(total = 0.3, file_dt = 0.01)
+moment_before = orig.totalKineticEnergy()
+orig.run(verbose=False)
+#orig.writeVTKall()
+orig.readlast(verbose=False)
+moment_after = orig.totalKineticEnergy()
+pytestutils.compareFloats(moment_before, moment_after,\
+ " 45 deg. wall collision:\t")
+
+## Visco-elastic collisions
+
+# Normal impact with normal viscous damping. Test that the lost kinetic energy
+# is saved as dissipated viscous energy
+orig = sphere.sim(np=1, sid='contactmodeltest')
+orig.radius[:] = 1.0
+orig.x[0,:] = [5.0, 5.0, 1.05]
+orig.setYoungsModulus(70.0e9)
+orig.vel[0,2] = -0.1
+orig.defineWorldBoundaries(L=[10,10,10])
+orig.gamma_wn[0] = 1.0e6
+orig.gamma_wt[0] = 0.0
+orig.initTemporal(total = 1.0, file_dt = 0.01)
+Ekin_before = orig.energy('kin')
+orig.run(verbose=False)
+#orig.writeVTKall()
+orig.readlast(verbose=False)
+Ekin_after = orig.energy('kin')
+Ev_after = orig.energy('visc_n')
+#print("Ekin_before = " + str(Ekin_before) + " J")
+#print("Ekin_after = " + str(Ekin_after) + " J")
+pytestutils.test(Ekin_before > Ekin_after,
+ "Viscoelastic normal wall collision (1/2):")
+pytestutils.compareFloats(Ekin_before, Ekin_after+Ev_after,\
+ "Viscoelastic normal wall collision (2/2):", tolerance=0.05)
+
+# Oblique impact: Check for conservation of momentum (sum(v_i*m_i))
+orig = sphere.sim(np=1, sid='contactmodeltest')
+orig.radius[:] = 1.0
+orig.x[0,:] = [5.0, 5.0, 1.05]
+orig.setYoungsModulus(70.0e9)
+orig.vel[0,2] = -0.1
+orig.vel[0,0] = 0.1
+orig.defineWorldBoundaries(L=[10,10,10])
+orig.gamma_wn[0] = 1.0e6
+orig.gamma_wt[0] = 1.0e6
+orig.initTemporal(total = 1.0, file_dt = 0.01)
+E_kin_before = orig.energy('kin')
+orig.run(verbose=False)
+#orig.writeVTKall()
+orig.readlast(verbose=False)
+#Ekin_after = orig.energy('kin')
+#Erot_after = orig.energy('rot')
+#Es_after = orig.energy('shear')
+#pytestutils.compareFloats(Ekin_before,\
+ #Ekin_after+Erot_after+Es_after,\
+ #" 45 deg. wall collision:", tolerance=0.03)
+pytestutils.test(Ekin_before > Ekin_after,
+ " 45 deg. wall collision (1/2):")
+pytestutils.test((orig.angvel[0,0] == 0.0 and orig.angvel[0,1] > 0.0 \
+ and orig.angvel[0,2] == 0.0),
+ " 45 deg. wall collision (2/2):")
+
+
+
+sphere.cleanup(orig)
(DIR) diff --git a/tests/contactmodel_young.py b/tests/contactmodel_young.py
t@@ -0,0 +1,232 @@
+#!/usr/bin/env python
+'''
+Validate the implemented contact models by observing the behavior of two
+particles.
+'''
+
+import sphere
+import numpy
+import pytestutils
+
+### Particle-particle interaction ##############################################
+
+## Linear elastic collisions
+
+# Normal impact: Check for conservation of momentum (sum(v_i*m_i))
+orig = sphere.sim(np=2, sid='contactmodeltest')
+after = sphere.sim(np=2, sid='contactmodeltest')
+sphere.cleanup(orig)
+#orig.radius[:] = [1.0, 2.0]
+orig.radius[:] = [1.0, 1.0]
+orig.x[0,:] = [5.0, 5.0, 2.0]
+orig.x[1,:] = [5.0, 5.0, 4.05]
+orig.setYoungsModulus(70.0e9)
+v_orig = 1
+orig.vel[0,2] = v_orig
+orig.defineWorldBoundaries(L=[10,10,10])
+orig.initTemporal(total = 0.1, file_dt = 0.01)
+
+orig.run(dry=True)
+orig.run(verbose=False)
+after.readlast(verbose=False)
+pytestutils.compareFloats(orig.vel[0,2], after.vel[1,2],\
+ "Elastic normal collision (1/4):")
+#print(orig.totalKineticEnergy())
+#print(after.totalKineticEnergy())
+pytestutils.compareFloats(orig.totalKineticEnergy(), after.totalKineticEnergy(),\
+ "Elastic normal collision (2/4):")
+
+# Normal impact with different sizes: Check for conservation of momentum
+orig = sphere.sim(np=2, sid='contactmodeltest')
+after = sphere.sim(np=2, sid='contactmodeltest')
+sphere.cleanup(orig)
+orig.radius[:] = [2.0, 1.0]
+orig.x[0,:] = [5.0, 5.0, 2.0]
+orig.x[1,:] = [5.0, 5.0, 5.05]
+orig.setYoungsModulus(70.0e9)
+orig.vel[0,2] = 1.0
+orig.defineWorldBoundaries(L=[10,10,10])
+orig.initTemporal(total = 0.1, file_dt = 0.01)
+
+orig.run(verbose=False)
+after.readlast(verbose=False)
+pytestutils.compareFloats(orig.totalKineticEnergy(), after.totalKineticEnergy(),\
+ "Elastic normal collision (3/4):")
+
+# Normal impact with different sizes: Check for conservation of momentum
+orig = sphere.sim(np=2, sid='contactmodeltest')
+after = sphere.sim(np=2, sid='contactmodeltest')
+sphere.cleanup(orig)
+orig.radius[:] = [1.0, 2.0]
+orig.x[0,:] = [5.0, 5.0, 2.0]
+orig.x[1,:] = [5.0, 5.0, 5.05]
+orig.setYoungsModulus(70.0e9)
+orig.vel[0,2] = 1.0
+orig.defineWorldBoundaries(L=[10,10,10])
+orig.initTemporal(total = 0.1, file_dt = 0.01)
+
+orig.run(verbose=False)
+after.readlast(verbose=False)
+pytestutils.compareFloats(orig.totalKineticEnergy(), after.totalKineticEnergy(),\
+ "Elastic normal collision (4/4):")
+
+
+## Linear viscous-elastic collisions
+
+# Normal impact: Check for conservation of momentum (sum(v_i*m_i))
+orig = sphere.sim(np=2, sid='contactmodeltest')
+after = sphere.sim(np=2, sid='contactmodeltest')
+sphere.cleanup(orig)
+orig.radius[:] = [1.0, 1.0]
+orig.x[0,:] = [5.0, 5.0, 2.0]
+orig.x[1,:] = [5.0, 5.0, 4.05]
+orig.setYoungsModulus(70.0e9)
+v_orig = 1
+orig.vel[0,2] = v_orig
+orig.defineWorldBoundaries(L=[10,10,10])
+orig.initTemporal(total = 0.1, file_dt = 0.01)
+orig.gamma_n[0] = 1.0e6
+
+orig.run(verbose=False)
+after.readlast(verbose=False)
+#print(orig.totalKineticEnergy())
+#print(after.totalKineticEnergy())
+#print(after.totalViscousEnergy())
+pytestutils.test(orig.vel[0,2] > after.vel[1,2],\
+ "Viscoelastic normal collision (1/4):")
+pytestutils.compareFloats(orig.totalKineticEnergy(),
+ after.totalKineticEnergy()
+ + after.totalViscousEnergy(),
+ "Viscoelastic normal collision (2/4):", tolerance=0.05)
+
+# Normal impact with different sizes: Check for conservation of momentum
+orig = sphere.sim(np=2, sid='contactmodeltest')
+after = sphere.sim(np=2, sid='contactmodeltest')
+sphere.cleanup(orig)
+orig.radius[:] = [2.0, 1.0]
+orig.x[0,:] = [5.0, 5.0, 2.0]
+orig.x[1,:] = [5.0, 5.0, 5.05]
+orig.setYoungsModulus(70.0e9)
+orig.vel[0,2] = 1.0
+orig.defineWorldBoundaries(L=[10,10,10])
+orig.initTemporal(total = 0.1, file_dt = 0.01)
+orig.gamma_n[0] = 1.0e6
+
+orig.run(verbose=False)
+after.readlast(verbose=False)
+pytestutils.compareFloats(orig.totalKineticEnergy(),
+ after.totalKineticEnergy()
+ + after.totalViscousEnergy(),
+ "Viscoelastic normal collision (3/4):", tolerance=0.05)
+
+# Normal impact with different sizes: Check for conservation of momentum
+orig = sphere.sim(np=2, sid='contactmodeltest')
+after = sphere.sim(np=2, sid='contactmodeltest')
+sphere.cleanup(orig)
+orig.radius[:] = [1.0, 2.0]
+orig.x[0,:] = [5.0, 5.0, 2.0]
+orig.x[1,:] = [5.0, 5.0, 5.05]
+orig.setYoungsModulus(70.0e9)
+orig.vel[0,2] = 1.0
+orig.defineWorldBoundaries(L=[10,10,10])
+orig.initTemporal(total = 0.1, file_dt = 0.01)
+orig.gamma_n[0] = 1.0e6
+
+orig.run(verbose=False)
+after.readlast(verbose=False)
+pytestutils.compareFloats(orig.totalKineticEnergy(),
+ after.totalKineticEnergy()
+ + after.totalViscousEnergy(),
+ "Viscoelastic normal collision (4/4):", tolerance=0.05)
+
+
+
+## Oblique elastic collisions
+
+# Normal impact, low angle, no slip
+orig = sphere.sim(np=2, sid='contactmodeltest')
+after = sphere.sim(np=2, sid='contactmodeltest')
+sphere.cleanup(orig)
+#orig.radius[:] = [1.0, 2.0]
+orig.radius[:] = [1.0, 1.0]
+orig.x[0,:] = [5.0, 5.0, 2.0]
+orig.x[1,:] = [5.0, 5.0, 4.05]
+orig.setYoungsModulus(70.0e9)
+orig.vel[0,2] = 1
+orig.vel[0,0] = 1
+orig.mu_s[0] = 1e9 # no slip
+orig.mu_d[0] = 1e9 # no slip
+orig.defineWorldBoundaries(L=[10,10,10])
+orig.initTemporal(total = 0.1, file_dt = 0.01)
+
+orig.run(verbose=False)
+after.readlast(verbose=False)
+pytestutils.test((after.angvel[:,1] < 0.0).all(),
+ "Oblique normal collision (1/8):")
+pytestutils.compareFloats(orig.totalKineticEnergy(),
+ after.totalKineticEnergy()
+ + after.totalRotationalEnergy(),
+ "Oblique normal collision (2/8):", tolerance=0.05)
+
+# Normal impact, low angle, slip
+orig = sphere.sim(np=2, sid='contactmodeltest')
+after = sphere.sim(np=2, sid='contactmodeltest')
+sphere.cleanup(orig)
+#orig.radius[:] = [1.0, 2.0]
+orig.radius[:] = [1.0, 1.0]
+orig.x[0,:] = [5.0, 5.0, 2.0]
+orig.x[1,:] = [5.0, 5.0, 4.05]
+orig.setYoungsModulus(70.0e9)
+orig.vel[0,2] = 1
+orig.vel[0,0] = 1
+orig.mu_s[0] = 0.3
+orig.mu_d[0] = 0.3
+orig.defineWorldBoundaries(L=[10,10,10])
+orig.initTemporal(total = 0.1, file_dt = 0.01)
+
+orig.run(verbose=False)
+after.readlast(verbose=False)
+pytestutils.compareFloats(orig.totalKineticEnergy(),
+ after.totalKineticEnergy()
+ + after.totalRotationalEnergy()
+ + after.totalFrictionalEnergy(),
+ "Oblique normal collision (3/8):", tolerance=0.05)
+pytestutils.test((after.angvel[:,1] < 0.0).all(),
+ "Oblique normal collision (4/8):")
+pytestutils.test(after.totalFrictionalEnergy() > 0.0,
+ "Oblique normal collision (5/8):")
+
+# Normal impact, low angle, slip, viscous damping tangentially
+orig = sphere.sim(np=2, sid='contactmodeltest')
+after = sphere.sim(np=2, sid='contactmodeltest')
+sphere.cleanup(orig)
+#orig.radius[:] = [1.0, 2.0]
+orig.radius[:] = [1.0, 1.0]
+orig.x[0,:] = [5.0, 5.0, 2.0]
+orig.x[1,:] = [5.0, 5.0, 4.05]
+orig.setYoungsModulus(70.0e9)
+orig.vel[0,2] = 1
+orig.vel[0,0] = 1
+orig.mu_s[0] = 0.3
+orig.mu_d[0] = 0.3
+orig.gamma_t[0] = 1.0e3
+orig.defineWorldBoundaries(L=[10,10,10])
+orig.initTemporal(total = 0.1, file_dt = 0.01)
+
+orig.run(verbose=False)
+after.readlast(verbose=False)
+print(after.totalViscousEnergy())
+pytestutils.compareFloats(orig.totalKineticEnergy(),
+ after.totalKineticEnergy()
+ + after.totalRotationalEnergy()
+ + after.totalFrictionalEnergy()
+ + after.totalViscousEnergy(),
+ "Oblique normal collision (6/8):", tolerance=0.05)
+pytestutils.test((after.angvel[:,1] < 0.0).all(),
+ "Oblique normal collision (7/8):")
+pytestutils.test(after.totalFrictionalEnergy() > 0.0,
+ "Oblique normal collision (8/8):")
+pytestutils.test(after.totalFrictionalEnergy() > 0.0,
+ "Oblique normal collision (8/8):")
+
+orig.cleanup()