tadded contact model tests - 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 a1bd00c040694d453ddfdc6742b1050d14ac2215
(DIR) parent f0459aec942c8b6de50bd358ce0dc1e5b588e4c4
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Fri, 20 Jun 2014 11:40:19 +0200
added contact model tests
Diffstat:
M python/sphere.py | 8 ++++----
M tests/CMakeLists.txt | 3 +++
A tests/contactmodel.py | 72 +++++++++++++++++++++++++++++++
M tests/contactmodel_wall.py | 8 ++++----
4 files changed, 83 insertions(+), 8 deletions(-)
---
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -3294,18 +3294,18 @@ class sim:
:param idx: The particle index
:type idx: int
:returns: The particle momentum [N*s]
+ :return type: numpy.array
'''
- return self.rho*V_sphere(self.radius[idx])\
- *numpy.linalg.norm(self.vel[idx,:])
+ return self.rho*V_sphere(self.radius[idx])*self.vel[idx,:]
def totalMomentum(self):
'''
Returns the sum of particle momentums.
:returns: The sum of particle momentums (m*v) [N*s]
- :return type: float
+ :return type: numpy.array
'''
- m_sum = 0.0
+ m_sum = numpy.zeros(3)
for i in range(self.np):
m_sum += self.momentum(i)
return m_sum
(DIR) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
t@@ -6,6 +6,9 @@ add_test(io_tests ${PYTHON_EXECUTABLE}
add_test(wall_contact_model_tests ${PYTHON_EXECUTABLE}
${CMAKE_CURRENT_BINARY_DIR}/contactmodel_wall.py)
+add_test(contact_model_tests ${PYTHON_EXECUTABLE}
+ ${CMAKE_CURRENT_BINARY_DIR}/contactmodel.py)
+
add_test(io_tests_fluid ${PYTHON_EXECUTABLE}
${CMAKE_CURRENT_BINARY_DIR}/io_tests_fluid.py)
(DIR) diff --git a/tests/contactmodel.py b/tests/contactmodel.py
t@@ -0,0 +1,72 @@
+#!/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]
+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(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.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.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):")
+
+
+
+
+#orig.cleanup()
(DIR) diff --git a/tests/contactmodel_wall.py b/tests/contactmodel_wall.py
t@@ -23,13 +23,13 @@ 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.totalMomentum()
+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.totalMomentum()
+moment_after = orig.totalKineticEnergy()
#print(moment_before)
#print(moment_after)
#print("time step: " + str(orig.time_dt[0]))
t@@ -47,11 +47,11 @@ 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)
-moment_before = orig.totalMomentum()
+moment_before = orig.totalKineticEnergy()
orig.run(verbose=False)
#orig.writeVTKall()
orig.readlast(verbose=False)
-moment_after = orig.totalMomentum()
+moment_after = orig.totalKineticEnergy()
pytestutils.compareFloats(moment_before, moment_after,\
" 45 deg. wall collision:\t")