tcleanup contactmodels.cuh, label in/out arrays, verbose runs - 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 f924918c5ac9717d5ebc6cb8143e854c7c445e95
(DIR) parent cebda4c5601b577b6535bd0f11bdccfd606a76d2
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Mon, 30 Jun 2014 10:13:07 +0200
cleanup contactmodels.cuh, label in/out arrays, verbose runs
Diffstat:
M src/contactmodels.cuh | 145 +------------------------------
M src/contactsearch.cuh | 45 ++++++++++++++++---------------
M tests/io_tests.py | 5 +++--
3 files changed, 29 insertions(+), 166 deletions(-)
---
(DIR) diff --git a/src/contactmodels.cuh b/src/contactmodels.cuh
t@@ -50,12 +50,6 @@ __device__ Float contactLinear_wall(Float3* F, Float3* T, Float* es_dot,
- 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",
- devC_params.k_n*delta,
- devC_params.gamma_wn,
- devC_params.gamma_wn*vel_n);*/
-
// Store the energy lost by viscous damping. See derivation in
// contactLinear()
*ev_dot += devC_params.gamma_wn * vel_n * vel_n;
t@@ -93,16 +87,6 @@ __device__ Float contactLinear_wall(Float3* F, Float3* T, Float* es_dot,
}
}
- /* if (angvel_length > 0.f) {
- // Apply rolling resistance (Zhou et al. 1999)
- //T_res = -angvel_a/angvel_length * devC_params.mu_r * radius_a * f_n_length;
-
- // New rolling resistance model
- T_res = -1.0f * fmin(devC_params.gamma_r * radius_a * angvel_length,
- devC_params.mu_r * radius_a * f_n_length)
- * angvel_a/angvel_length;
- }*/
-
// Total force from wall
*F += f_n + f_t;
t@@ -113,7 +97,6 @@ __device__ Float contactLinear_wall(Float3* F, Float3* T, Float* es_dot,
*p += f_n_length / (4.0f * PI * radius_a*radius_a);
// Return force excerted onto the wall
- //return -dot(*F, n);
return dot(f_n, n);
}
t@@ -142,7 +125,6 @@ __device__ void contactLinearViscous(Float3* F, Float3* T,
// Force between grain pair decomposed into normal- and tangential part
Float3 f_n, f_t, f_c;
- //Float3 T_res;
// Normal vector of contact
Float3 n_ab = x_ab/x_ab_length;
t@@ -170,31 +152,6 @@ __device__ void contactLinearViscous(Float3* F, Float3* T,
Float3 vel_t_ab = vel_ab - (n_ab * dot(vel_ab, n_ab));
Float vel_t_ab_length = length(vel_t_ab);
- // Compute the normal stiffness of the contact
- //Float k_n_ab = k_n_a * k_n_b / (k_n_a + k_n_b);
-
- // Calculate rolling radius
- //Float R_bar = (radius_a + radius_b) / 2.0f;
-
- // Normal force component: linear-elastic approximation (Augier 2009, eq. 3)
- // with velocity dependant damping
- // Damping coefficient: alpha = 0.8
- //f_n = (-k_n_ab * delta_ab + 2.0f * 0.8f * sqrtf(m_eff*k_n_ab) * vel_ab) * n_ab;
-
- // Linear spring for normal component (Renzo 2004, eq. 35)
- // Dissipation due to plastic deformation is modelled by using a different
- // unloading spring constant (Walton and Braun 1986)
- // Here the factor in the second term determines the relative strength of the
- // unloading spring relative to the loading spring.
- /* if (vel_n_ab > 0.0f) { // Loading
- f_n = (-k_n_ab * delta_ab) * n_ab;
- } else { // Unloading
- f_n = (-k_n_ab * 0.90f * delta_ab) * n_ab;
- } // f_n is OK! */
-
- // Normal force component: Elastic
- //f_n = -devC_params.k_n * delta_ab * n_ab;
-
// Normal force component: Elastic - viscous damping
f_n = fmax(0.0, -devC_params.k_n * delta_ab
- devC_params.gamma_n * vel_n_ab) * n_ab;
t@@ -210,8 +167,7 @@ __device__ void contactLinearViscous(Float3* F, Float3* T,
f_c = -kappa * sqrtf(radius_a * radius_b) * n_ab;
// Initialize force vectors to zero
- f_t = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
- //T_res = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
+ f_t = MAKE_FLOAT3(0.0, 0.0, 0.0);
// Shear force component: Nonlinear relation
// Coulomb's law of friction limits the tangential force to less or equal
t@@ -242,20 +198,8 @@ __device__ void contactLinearViscous(Float3* F, Float3* T,
}
}
- /* if (angvel_ab_length > 0.f) {
- // Apply rolling resistance (Zhou et al. 1999)
- //T_res = -angvel_ab/angvel_ab_length * devC_params.mu_r * R_bar * length(f_n);
-
- // New rolling resistance model
- T_res = -1.0f * fmin(devC_params.gamma_r * R_bar * angvel_ab_length,
- devC_params.mu_r * R_bar * f_n_length)
- * angvel_ab/angvel_ab_length;
- }
- */
-
// Add force components from this collision to total force for particle
*F += f_n + f_t + f_c;
- //*T += -(radius_a + delta_ab/2.0f) * cross(n_ab, f_t) + T_res;
*T += -(radius_a + delta_ab/2.0f) * cross(n_ab, f_t);
// Pressure excerted onto the particle from this contact
t@@ -265,32 +209,6 @@ __device__ void contactLinearViscous(Float3* F, Float3* T,
// Linear elastic contact model for particle-particle interactions
-/*__device__ void contactLinear_bck(Float3* F, Float3* T,
- Float* es_dot, Float* ev_dot, Float* p,
- unsigned int idx_a_orig,
- unsigned int idx_b_orig,
- Float4 vel_a,
- Float4* dev_vel,
- Float3 angvel_a,
- Float4* dev_angvel,
- Float radius_a, Float radius_b,
- Float3 x_ab, Float x_ab_length,
- Float delta_ab, Float4* dev_delta_t,
- unsigned int mempos)
- {
- Float4 vel_b = dev_vel[idx_b_orig];
- Float4 angvel4_b = dev_vel[idx_b_orig];
-
-// Fe
-
-
-
-
-}*/
-
-
-
-// Linear elastic contact model for particle-particle interactions
__device__ void contactLinear(Float3* F, Float3* T,
Float* es_dot, Float* ev_dot, Float* p,
unsigned int idx_a_orig,
t@@ -325,7 +243,6 @@ __device__ void contactLinear(Float3* F, Float3* T,
// Force between grain pair decomposed into normal- and tangential part
Float3 f_n, f_t, f_c;
- //Float3 T_res;
// Normal vector of contact
Float3 n = x / x_length;
t@@ -348,33 +265,22 @@ __device__ void contactLinear(Float3* F, Float3* T,
Float angvel_length = length(angvel);
// Normal component of the relative contact interface velocity
- //Float vel_n = dot(vel_linear, n);
Float vel_n = -dot(vel_linear, n);
// Tangential component of the relative contact interface velocity
// Hinrichsen and Wolf 2004, eq. 13.9
- //Float3 vel_t = vel - vel_n * n;
Float3 vel_t = vel - n * dot(n, vel);
Float vel_t_length = length(vel_t);
// Correct tangential displacement vector, which is
// necessary if the tangential plane rotated
- //Float3 delta_t0 = delta_t0_uncor - (n * dot(delta_t0_uncor, n));
Float3 delta_t0 = delta_t0_uncor - (n * dot(n, delta_t0_uncor));
- //cuPrintf("delta_t0: %f\t%f\t%f\n", delta_t0.x, delta_t0.y, delta_t0.z);
Float delta_t0_length = length(delta_t0);
// New tangential displacement vector
Float3 delta_t;
- // Compute the normal stiffness of the contact
- //Float k_n_ab = k_n_a * k_n_b / (k_n_a + k_n_b);
-
- // Normal force component: Elastic
- //f_n = -devC_params.k_n * delta * n_ab;
-
// Normal force component: Elastic - viscous damping
- //f_n = (-devC_params.k_n * delta - devC_params.gamma_n * vel_n) * n;
f_n = fmax(0.0, -devC_params.k_n*delta + devC_params.gamma_n * vel_n) * n;
Float f_n_length = length(f_n);
t@@ -390,9 +296,6 @@ __device__ void contactLinear(Float3* F, Float3* T,
// Initialize force vectors to zero
f_t = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
- //T_res = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
-
- //cuPrintf("mu_s = %f\n", devC_params.mu_s);
// Apply a tangential force if the previous tangential displacement
// is non-zero, or the current sliding velocity is non-zero.
t@@ -402,7 +305,6 @@ __device__ void contactLinear(Float3* F, Float3* T,
delta_t = delta_t0 + vel_t * devC_dt;
// Tangential force: Visco-Elastic, before limitation criterion
- //Float3 f_t_elast = -devC_params.k_t * delta_t0;
Float3 f_t_elast = -devC_params.k_t * delta_t;
Float3 f_t_visc = -devC_params.gamma_t * vel_t;
f_t = f_t_elast + f_t_visc;
t@@ -422,27 +324,11 @@ __device__ void contactLinear(Float3* F, Float3* T,
// resulting in a slip and energy dissipation
if (f_t_length > f_t_limit) { // Static friciton exceeded: Dynamic case
- //cuPrintf("slip! %f > %f\n", f_t_length, f_t_limit);
-
// tangential vector
Float3 t = f_t/length(f_t);
// Frictional force is reduced to equal the dynamic limit
- //f_t *= (devC_params.mu_d * length(f_n-f_c))/f_t_length;
f_t = f_t_limit * t;
- //f_t = f_t * (devC_params.mu_d * f_n_length)/f_t;
-
- // A slip event zeros the displacement vector
- //delta_t = MAKE_FLOAT3(0.0f, 0.0f, 0.0f);
-
- // In a slip event, the tangential spring is adjusted to a
- // length which is consistent with Coulomb's equation
- // (Hinrichsen and Wolf, 2004)
- //delta_t = -1.0f/devC_params.k_t * (f_t + devC_params.gamma_t * vel_t_ab);
- //delta_t = -1.0f/devC_params.k_t * f_t;
- //delta_t = -1.0/devC_params.k_t * f_t + devC_params.gamma_t * vel_t_ab;
- //delta_t = -1.0/devC_params.k_t * devC_params.mu_d * t +
- //+ devC_params.gamma_t * vel_t;
// In the sliding friction case, the tangential spring is adjusted
// to a length consistent with Coulombs (dynamic) condition (Luding
t@@ -453,41 +339,17 @@ __device__ void contactLinear(Float3* F, Float3* T,
// Shear friction heat production rate:
// The energy lost from the tangential spring is dissipated as heat
- //*es_dot += -dot(vel_t_ab, f_t);
- //*es_dot += length(delta_t0 - delta_t) * devC_params.k_t / devC_dt;
*es_dot += 0.5*length(length(f_t) * vel_t * devC_dt) / devC_dt; // Seen in ESyS-Particle
- //*es_dot += fabs(dot(delta_t0 - delta_t, f_t)) / devC_dt;
-
- } //else { // Static case
- //cuPrintf("no slip: %f < %f\n", f_t_length, f_t_limit);
- // No correction of f_t is required
-
- // Add tangential displacement to total tangential displacement
- //delta_t = delta_t0 + vel_t_ab * devC_dt;
- //}
+ }
}
- //if (angvel_ab_length > 0.f) {
- // Apply rolling resistance (Zhou et al. 1999)
- //T_res = -angvel_ab/angvel_ab_length * devC_params.mu_r * R_bar * length(f_n);
-
- // New rolling resistance model
- /*T_res = -1.0f * fmin(devC_params.gamma_r * R_bar * angvel_ab_length,
- devC_params.mu_r * R_bar * f_n_length)
- * angvel_ab/angvel_ab_length;*/
- //T_res = -1.0f * fmin(devC_params.gamma_r * radius_a * angvel_ab_length,
- // devC_params.mu_r * radius_a * f_n_length)
- // * angvel_ab/angvel_ab_length;
- //}
-
// Add force components from this collision to total force for particle
*F += f_n + f_t + f_c;
+
// Add torque components from this collision to total torque for particle
// Comment out the line below to disable rotation
- //*T += -(radius_a + delta_ab/2.0f) * cross(n_ab, f_t) + T_res;
- //*T += cross(-(radius_a + delta*0.5) * n_ab, f_t) + T_res;
*T += cross(-(radius_a + delta*0.5) * n, f_t);
// Pressure excerted onto the particle from this contact
t@@ -583,7 +445,6 @@ __device__ void contactHertz(Float3* F, Float3* T,
// watt = N*m/s = N*s/m * m/s * m/s * s / s
*ev_dot += devC_params.gamma_n * vel_n_ab * vel_n_ab;
-
// 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_ab) < 0.0f)
(DIR) diff --git a/src/contactsearch.cuh b/src/contactsearch.cuh
t@@ -371,28 +371,29 @@ __global__ void topology(unsigned int* dev_cellStart,
// Collide with top- and bottom walls, save resulting force on upper wall.
// Kernel is executed on device, and is callable from host only.
// Function is called from mainGPU loop.
-__global__ void interact(unsigned int* dev_gridParticleIndex, // Input: Unsorted-sorted key
- unsigned int* dev_cellStart,
- unsigned int* dev_cellEnd,
- Float4* dev_x,
- Float4* dev_x_sorted,
- Float4* dev_vel_sorted,
- Float4* dev_angvel_sorted,
- Float4* dev_vel,
- Float4* dev_angvel,
- Float4* dev_force,
- Float4* dev_torque,
- Float* dev_es_dot,
- Float* dev_ev_dot,
- Float* dev_es,
- Float* dev_ev,
- Float* dev_p,
- Float4* dev_walls_nx,
- Float4* dev_walls_mvfd,
- Float* dev_walls_force_pp, //uint4* dev_bonds_sorted,
- unsigned int* dev_contacts,
- Float4* dev_distmod,
- Float4* dev_delta_t)
+__global__ void interact(
+ unsigned int* dev_gridParticleIndex, // in
+ unsigned int* dev_cellStart, // in
+ unsigned int* dev_cellEnd, // in
+ Float4* dev_x, // in
+ Float4* dev_x_sorted, // in
+ Float4* dev_vel_sorted, // in
+ Float4* dev_angvel_sorted, // in
+ Float4* dev_vel, // in
+ Float4* dev_angvel, // in
+ Float4* dev_force, // out
+ Float4* dev_torque, // out
+ Float* dev_es_dot, // out
+ Float* dev_ev_dot, // out
+ Float* dev_es, // out
+ Float* dev_ev, // out
+ Float* dev_p, // out
+ Float4* dev_walls_nx, // in
+ Float4* dev_walls_mvfd, // in
+ Float* dev_walls_force_pp, // out
+ unsigned int* dev_contacts, // in
+ Float4* dev_distmod, // in
+ Float4* dev_delta_t) // in
{
// Thread index equals index of particle A
unsigned int idx_a = blockIdx.x * blockDim.x + threadIdx.x;
(DIR) diff --git a/tests/io_tests.py b/tests/io_tests.py
t@@ -25,8 +25,9 @@ py.readbin("../input/" + orig.sid + ".bin", verbose=False)
compare(orig, py, "Python IO:")
# Test C++ IO routines
-#orig.run(verbose=False, hideinputfile=True)
-orig.run(verbose=True, hideinputfile=True)
+#orig.run(verbose=True, hideinputfile=True)
+orig.run(dry=True)
+orig.run()
cpp = sphere.sim()
cpp.readbin("../output/" + orig.sid + ".output00000.bin", verbose=False)
compare(orig, cpp, "C++ IO: ")