tadded test of kinetic energy magnitude, corrected velocity - 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 624061bd7a2a2b9424853e0848aabc2942de873d
(DIR) parent 64e194a6c52ee909e414d69c73567254683d0544
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Fri, 20 Jun 2014 09:49:36 +0200
added test of kinetic energy magnitude, corrected velocity
Diffstat:
M src/contactmodels.cuh | 17 ++---------------
M tests/contactmodel_wall.py | 6 +++++-
2 files changed, 7 insertions(+), 16 deletions(-)
---
(DIR) diff --git a/src/contactmodels.cuh b/src/contactmodels.cuh
t@@ -32,7 +32,7 @@ __device__ Float contactLinear_wall(Float3* F, Float3* T, Float* es_dot,
// Contact velocity is the sum of the linear and
// rotational components
- //Float3 vel = linvel + radius_a * cross(n, angvel) + wvel;
+ //Float3 vel = vel_linear + radius_a*cross(n, angvel) + wvel;
Float3 vel = vel_linear + (radius_a + delta/2.0) * cross(n, angvel) + wvel;
// Normal component of the contact velocity
t@@ -42,18 +42,13 @@ __device__ Float contactLinear_wall(Float3* F, Float3* T, Float* es_dot,
// The tangential velocity is the contact velocity
// with the normal component subtracted
- //Float3 vel_t = vel - n * (dot(vel, n));
const Float3 vel_t = vel - n * (dot(n, vel));
const Float vel_t_length = length(vel_t);
- // Calculate elastic normal component
- //Float3 f_n = -devC_params.k_n * delta * n;
-
// Normal force component: Elastic - viscous damping
- //Float3 f_n = (-devC_params.k_n * delta - devC_params.gamma_wn * vel_n) * n;
- //Float3 f_n = (-devC_params.k_n * delta + devC_params.gamma_wn * vel_n) * n;
Float3 f_n = fmax(0.0, -devC_params.k_n*delta
- devC_params.gamma_wn*vel_n) * n;
+ const Float f_n_length = length(f_n); // Save length for later use
// Print data for contact model validation
/*printf("f_n_elast = %f\tgamma_wn = %f\tf_n_visc = %f\n",
t@@ -65,14 +60,6 @@ __device__ Float contactLinear_wall(Float3* F, Float3* T, Float* es_dot,
// contactLinear()
*ev_dot += devC_params.gamma_wn * vel_n * vel_n;
- // Make sure the viscous damping doesn't exceed the elastic component,
- // i.e. the damping factor doesn't exceed the critical damping:
- // 2*sqrt(m*k_n)
- if (dot(f_n, n) < 0.0)
- f_n = MAKE_FLOAT3(0.0, 0.0, 0.0);
-
- const Float f_n_length = length(f_n); // Save length for later use
-
// Initialize vectors
Float3 f_t = MAKE_FLOAT3(0.0, 0.0, 0.0);
(DIR) diff --git a/tests/contactmodel_wall.py b/tests/contactmodel_wall.py
t@@ -73,8 +73,12 @@ orig.run(verbose=False)
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:", tolerance=0.03)
+ "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')