tBond model passes stability, normal, shear, and bend tests, but not twist - 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 c2a5b8e5819c09eea3dba1b97fcade073e984bd9
(DIR) parent 11c8ccc15c6f107407771020cf27c7bd34233df4
(HTM) Author: Anders Damsgaard <adc@geo.au.dk>
Date: Thu, 7 Mar 2013 18:36:25 +0100
Bond model passes stability, normal, shear, and bend tests, but not twist
Diffstat:
M src/cohesion.cuh | 44 ++++++++++++++++++++-----------
M tests/CMakeLists.txt | 3 +++
A tests/bond_tests.py | 187 +++++++++++++++++++++++++++++++
3 files changed, 219 insertions(+), 15 deletions(-)
---
(DIR) diff --git a/src/cohesion.cuh b/src/cohesion.cuh
t@@ -53,15 +53,18 @@ __global__ void bondsLinear(
const Float4 x_j = dev_x[bond.y];
const Float4 vel_i = dev_vel[bond.x];
const Float4 vel_j = dev_vel[bond.y];
- const Float4 angvel_i = dev_angvel[bond.x];
- const Float4 angvel_j = dev_angvel[bond.y];
+ const Float4 angvel4_i = dev_angvel[bond.x];
+ const Float4 angvel4_j = dev_angvel[bond.y];
+
+ const Float3 angvel_i = MAKE_FLOAT3(angvel4_i.x, angvel4_i.y, angvel4_i.z);
+ const Float3 angvel_j = MAKE_FLOAT3(angvel4_j.x, angvel4_j.y, angvel4_j.z);
// Initialize force- and torque vectors
Float3 f, t, f_n, f_t, t_n, t_t;
//// Bond geometry and inertia
-
+
// Parallel-bond radius (Potyondy and Cundall 2004, eq. 12)
const Float R_bar = devC_params.lambda_bar * fmin(x_i.w, x_j.w);
t@@ -72,7 +75,8 @@ __global__ void bondsLinear(
const Float I = 0.25 * PI * R_bar*R_bar*R_bar*R_bar;
// Bond polar moment of inertia (Potyondy and Cundall 2004, eq. 15)
- const Float J = 0.50 * PI * R_bar*R_bar*R_bar*R_bar;
+ //const Float J = 0.50 * PI * R_bar*R_bar*R_bar*R_bar;
+ const Float J = I*2.0;
// Inter-particle vector
const Float3 x = MAKE_FLOAT3(
t@@ -80,6 +84,9 @@ __global__ void bondsLinear(
x_i.y - x_j.y,
x_i.z - x_j.z);
const Float x_length = length(x);
+
+ // Find overlap (negative value if overlapping)
+ const Float overlap = fmin(0.0, x_length - (x_i.w + x_j.w));
// Normal vector of contact (points from i to j)
const Float3 n = x/x_length;
t@@ -91,11 +98,15 @@ __global__ void bondsLinear(
//const Float3 delta_t0 = delta_t0_uncor - dot(delta_t0_uncor, n);
const Float3 delta_t0 = delta_t0_uncor - (n * dot(n, delta_t0_uncor));
- // Contact displacement (should this include rolling?)
- const Float3 ddelta = MAKE_FLOAT3(
- vel_i.x - vel_j.x,
- vel_i.y - vel_j.y,
- vel_i.z - vel_j.z) * devC_dt;
+ // Contact displacement
+ const Float3 ddelta = (
+ MAKE_FLOAT3(
+ vel_i.x - vel_j.x,
+ vel_i.y - vel_j.y,
+ vel_i.z - vel_j.z)
+ + (x_i.w + overlap/2.0) * cross(n, angvel_i)
+ + (x_j.w + overlap/2.0) * cross(n, angvel_j)
+ ) * devC_dt;
// Normal component of the displacement increment
//const Float ddelta_n = dot(ddelta, n);
t@@ -111,11 +122,14 @@ __global__ void bondsLinear(
// Tangential component of the total displacement
const Float3 delta_t = delta_t0 + ddelta_t;
- // Normal force: Elastic contact model
+ // Normal force: Visco-elastic contact model
+ // The elastic component caused by the overlap is subtracted.
//f_n = devC_params.k_n * A * delta_n * n;
f_n = (devC_params.k_n * A * delta_n + devC_params.gamma_n * ddelta_n/devC_dt) * n;
+ //f_n += devC_params.k_n * overlap * n;
+
- // Tangential force: Elastic contact model
+ // Tangential force: Visco-elastic contact model
//f_t = -devC_params.k_t * A * delta_t;
f_t = -devC_params.k_t * A * delta_t - devC_params.gamma_t * ddelta_t/devC_dt;
t@@ -149,11 +163,11 @@ __global__ void bondsLinear(
// Tangential component of the total displacement
const Float3 omega_t = omega_t0 + domega_t;
- // Twisting torque: Elastic contact model
- //t_n = -devC_params.k_t * J * omega_n * n;
- t_n = (devC_params.k_t * J * omega_n + devC_params.gamma_t * domega_n/devC_dt) * n;
+ // Twisting torque: Visco-elastic contact model
+ t_n = -devC_params.k_t * J * omega_n * n;
+ //t_n = (devC_params.k_t * J * omega_n + devC_params.gamma_t * domega_n/devC_dt) * n;
- // Bending torque: Elastic contact model
+ // Bending torque: Visco-elastic contact model
//t_t = -devC_params.k_n * I * omega_t;
//t_t = -devC_params.k_n * I * omega_t - devC_params.gamma_n * domega_t/devC_dt;
t_t = devC_params.k_n * I * omega_t - devC_params.gamma_n * domega_t/devC_dt;
(DIR) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
t@@ -9,3 +9,6 @@ add_test(porosity_tests ${PYTHON_EXECUTABLE}
add_test(memory_tests ${PYTHON_EXECUTABLE}
${CMAKE_CURRENT_BINARY_DIR}/memcheck_tests.py)
+
+add_test(bond_tests ${PYTHON_EXECUTABLE}
+ ${CMAKE_CURRENT_BINARY_DIR}/bond_tests.py)
(DIR) diff --git a/tests/bond_tests.py b/tests/bond_tests.py
t@@ -0,0 +1,187 @@
+#!/usr/bin/env python
+from pytestutils import *
+from sphere import *
+
+def printKinematics(sb):
+ print('bonds_delta_n'); print(sb.bonds_delta_n)
+ print('bonds_delta_t'); print(sb.bonds_delta_t)
+ print('bonds_omega_n'); print(sb.bonds_omega_n)
+ print('bonds_omega_t'); print(sb.bonds_omega_t)
+ print('force'); print(sb.force)
+ print('torque'); print(sb.torque)
+ print('vel'); print(sb.vel)
+ print('angvel'); print(sb.angvel)
+
+#### Bond tests ####
+print("### Bond tests ###")
+
+# Zero arrays
+z2_1 = numpy.zeros((2,1))
+z2_2 = numpy.zeros((2,2))
+z1_3 = numpy.zeros((1,3))
+z2_3 = numpy.zeros((2,3))
+
+# Small value arrays
+smallval = 1e-8
+s2_1 = numpy.ones((2,1))*smallval
+
+# Inter-particle distances to try (neg. for overlap)
+distances = [0.2, 0.0, -0.2]
+#distances = [-0.2]
+
+for d in distances:
+
+ radii = 0.5
+ print("## Inter-particle distance: " + str(d/radii) + " radii")
+
+ sb = Spherebin(np=2, sid='bondtest')
+ cleanup(sb)
+
+ # setup particles, bond, and simulation
+ sb.x[0,:] = numpy.array((2.0, 2.0, 2.0))
+ #sb.x[1,:] = numpy.array((2.2, 2.0, 2.0))
+ sb.x[1,:] = numpy.array((2.0+2.0*radii+d, 2.0, 2.0))
+ sb.radius = numpy.ones(sb.np)*radii
+ sb.initGridAndWorldsize(margin = 10, periodic = 1, contactmodel = 2, g = numpy.array([0.0, 0.0, 0.0]))
+ sb.bond(0, 1)
+ sb.defaultParams(gamma_n = 0.0, gamma_t = 0.0)
+ #sb.initTemporal(total=0.5, file_dt=0.01)
+
+ #sb.render(verbose=False)
+ #visualize(sb.sid, "energy")
+
+ print("# Stability test")
+ sb.zeroKinematics()
+ sb.initTemporal(total=0.2, file_dt=0.01)
+ #sb.initTemporal(total=0.01, file_dt=0.0001)
+ sb.writebin(verbose=False)
+ sb.run(verbose=False)
+ #sb.run()
+ sb.readlast(verbose=False)
+ compareFloats(0.0, sb.energy("kin") + sb.energy("rot") + sb.energy("bondpot"), "Energy cons.")
+ compareNumpyArrays(sb.vel, z2_3, "vel\t")
+ compareNumpyArrays(sb.angvel, z2_3, "angvel\t")
+ #printKinematics(sb)
+ #visualize(sb.sid, "energy")
+ #sb.readbin('../output/' + sb.sid + '.output00001.bin')
+ #printKinematics(sb)
+
+ print("# Normal expansion")
+ sb.zeroKinematics()
+ #sb.initTemporal(total=0.2, file_dt=0.01)
+ sb.initTemporal(total=0.2, file_dt=0.001)
+ sb.vel[1,0] = 1e-4
+ Ekinrot0 = sb.energy("kin") + sb.energy("rot")
+ sb.writebin(verbose=False)
+ sb.run(verbose=False)
+ sb.readlast(verbose=False)
+ compareFloats(Ekinrot0, sb.energy("kin") + sb.energy("rot") + sb.energy("bondpot"), "Energy cons.")
+ print("vel[:,0]"),
+ if ((sb.vel[0,0] <= 0.0) or (sb.vel[1,0] <= 0.0)):
+ print(failed())
+ else :
+ print(passed())
+ compareNumpyArrays(sb.vel[:,1:2], z2_2, "vel[:,1:2]")
+ compareNumpyArrays(sb.angvel, z2_3, "angvel\t")
+ #printKinematics(sb)
+ #visualize(sb.sid, "energy")
+
+ print("# Normal contraction")
+ sb.zeroKinematics()
+ sb.initTemporal(total=0.2, file_dt=0.01)
+ sb.vel[1,0] = -1e-4
+ Ekinrot0 = sb.energy("kin") + sb.energy("rot")
+ sb.writebin(verbose=False)
+ sb.run(verbose=False)
+ sb.readlast(verbose=False)
+ compareFloats(Ekinrot0, sb.energy("kin") + sb.energy("rot") + sb.energy("bondpot"), "Energy cons.")
+ print("vel[:,0]"),
+ if ((sb.vel[0,0] >= 0.0) or (sb.vel[1,0] >= 0.0)):
+ print(failed())
+ else :
+ print(passed())
+ compareNumpyArrays(sb.vel[:,1:2], z2_2, "vel[:,1:2]")
+ compareNumpyArrays(sb.angvel, z2_3, "angvel\t")
+ #printKinematics(sb)
+
+ print("# Shear")
+ sb.zeroKinematics()
+ sb.initTemporal(total=0.2, file_dt=0.01)
+ sb.vel[1,2] = 1e-4
+ Ekinrot0 = sb.energy("kin") + sb.energy("rot")
+ sb.writebin(verbose=False)
+ sb.run(verbose=False)
+ sb.readlast(verbose=False)
+ compareFloats(Ekinrot0, sb.energy("kin") + sb.energy("rot") + sb.energy("bondpot"), "Energy cons.")
+ #print("vel[:,0]"),
+ #if ((sb.vel[0,0] <= 0.0) or (sb.vel[1,0] >= 0.0)):
+ # print(failed())
+ #else :
+ # print(passed())
+ compareNumpyArrays(sb.vel[:,1], z2_1, "vel[:,1]")
+ print("vel[:,2]"),
+ if ((sb.vel[0,2] <= 0.0) or (sb.vel[1,2] <= 0.0)):
+ print(failed())
+ else :
+ print(passed())
+ #compareNumpyArrays(sb.angvel[:,0:2:2], z2_2, "angvel[:,0:2:2]")
+ #print("angvel[:,1]"),
+ #if ((sb.angvel[0,1] >= 0.0) or (sb.angvel[1,1] >= 0.0)):
+ # print(failed())
+ #else :
+ # print(passed())
+ compareNumpyArrays(sb.angvel, z2_3, "angvel\t")
+ #printKinematics(sb)
+ #visualize(sb.sid, "energy")
+
+
+ print("# Bend")
+ sb.zeroKinematics()
+ sb.initTemporal(total=0.2, file_dt=0.01)
+ sb.angvel[1,1] = 1e-4
+ Ekinrot0 = sb.energy("kin") + sb.energy("rot")
+ sb.writebin(verbose=False)
+ sb.run(verbose=False)
+ sb.readlast(verbose=False)
+ compareFloats(Ekinrot0, sb.energy("kin") + sb.energy("rot") + sb.energy("bondpot"), "Energy cons.")
+ print("vel[:,0]"),
+ if (sb.vel[0,0] <= 0.0) or (sb.vel[1,0] >= 0.0):
+ print(failed())
+ else :
+ print(passed())
+ compareNumpyArrays(sb.vel[:,1], z2_1, "vel[:,1]")
+ print("vel[:,2]"),
+ if ((sb.vel[0,2] <= 0.0) or (sb.vel[1,2] >= 0.0)):
+ print(failed())
+ else :
+ print(passed())
+ compareNumpyArrays(sb.angvel[:,0:2:2], z2_2, "angvel[:,0:2:2]")
+ print("angvel[:,1]"),
+ if ((sb.angvel[0,1] <= 0.0) or (sb.angvel[1,1] <= 0.0)):
+ print(failed())
+ else :
+ print(passed())
+ #printKinematics(sb)
+ #visualize(sb.sid, "energy")
+
+
+ print("# Twist")
+ sb.zeroKinematics()
+ sb.initTemporal(total=0.2, file_dt=0.01)
+ sb.angvel[1,0] = 1e-4
+ Ekinrot0 = sb.energy("kin") + sb.energy("rot")
+ sb.writebin(verbose=False)
+ sb.run(verbose=False)
+ sb.readlast(verbose=False)
+ compareFloats(Ekinrot0, sb.energy("kin") + sb.energy("rot") + sb.energy("bondpot"), "Energy cons.")
+ compareNumpyArrays(sb.vel, z2_3, "vel\t")
+ print("angvel[:,1]"),
+ if ((sb.angvel[0,1] <= 0.0) or (sb.angvel[0,1] <= 0.0)):
+ print(failed())
+ else :
+ print(passed())
+ compareNumpyArrays(sb.angvel[:,0:2:2], z2_2, "angvel[:,0:2:2]")
+ printKinematics(sb)
+ visualize(sb.sid, "energy")
+
+#cleanup(sb)