tadd debugging info to stdout - 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 44093ce85fa78cecda72209b859fb2a3f3f610ed
(DIR) parent feec5954462d3a34b0a2aa9e142373e9a85e0f9f
(HTM) Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Fri, 10 Feb 2023 14:33:01 +0100
add debugging info to stdout
Diffstat:
M src/contactsearch.cuh | 85 ++++++++++++++++---------------
1 file changed, 44 insertions(+), 41 deletions(-)
---
(DIR) diff --git a/src/contactsearch.cuh b/src/contactsearch.cuh
t@@ -89,21 +89,21 @@ __device__ int findDistMod(int3* targetCell, Float3* distmod)
// Kernel executed on device, and callable from device only.
// Function is called from interact().
__device__ void findAndProcessContactsInCell(
- int3 targetCell,
- const unsigned int idx_a,
+ int3 targetCell,
+ const unsigned int idx_a,
const Float4 x_a,
const Float radius_a,
Float3* F,
- Float3* T,
+ Float3* T,
Float* es_dot,
Float* ev_dot,
Float* p,
- const Float4* __restrict__ dev_x_sorted,
- const Float4* __restrict__ dev_vel_sorted,
+ const Float4* __restrict__ dev_x_sorted,
+ const Float4* __restrict__ dev_vel_sorted,
const Float4* __restrict__ dev_angvel_sorted,
- const unsigned int* __restrict__ dev_cellStart,
+ const unsigned int* __restrict__ dev_cellStart,
const unsigned int* __restrict__ dev_cellEnd,
- const Float4* __restrict__ dev_walls_nx,
+ const Float4* __restrict__ dev_walls_nx,
Float4* __restrict__ dev_walls_mvfd)
//uint4 bonds)
{
t@@ -118,7 +118,7 @@ __device__ void findAndProcessContactsInCell(
// Calculate linear cell ID
unsigned int cellID = targetCell.x + targetCell.y * devC_grid.num[0]
- + (devC_grid.num[0] * devC_grid.num[1]) * targetCell.z;
+ + (devC_grid.num[0] * devC_grid.num[1]) * targetCell.z;
// Lowest particle index in cell
unsigned int startIdx = dev_cellStart[cellID];
t@@ -139,8 +139,8 @@ __device__ void findAndProcessContactsInCell(
Float kappa = devC_params.kappa;
// Distance between particle centers (Float4 -> Float3)
- Float3 x_ab = MAKE_FLOAT3(x_a.x - x_b.x,
- x_a.y - x_b.y,
+ Float3 x_ab = MAKE_FLOAT3(x_a.x - x_b.x,
+ x_a.y - x_b.y,
x_a.z - x_b.z);
// Adjust interparticle vector if periodic boundary/boundaries
t@@ -149,21 +149,21 @@ __device__ void findAndProcessContactsInCell(
Float x_ab_length = length(x_ab);
// Distance between particle perimeters
- Float delta_ab = x_ab_length - (radius_a + radius_b);
+ Float delta_ab = x_ab_length - (radius_a + radius_b);
// Check for particle overlap
if (delta_ab < 0.0f) {
- contactLinearViscous(F, T, es_dot, ev_dot, p,
+ contactLinearViscous(F, T, es_dot, ev_dot, p,
idx_a, idx_b,
- dev_vel_sorted,
+ dev_vel_sorted,
dev_angvel_sorted,
- radius_a, radius_b,
+ radius_a, radius_b,
x_ab, x_ab_length,
delta_ab, kappa);
- } else if (delta_ab < devC_params.db) {
+ } else if (delta_ab < devC_params.db) {
// Check wether particle distance satisfies the
// capillary bond distance
- capillaryCohesion_exp(F, radius_a, radius_b, delta_ab,
+ capillaryCohesion_exp(F, radius_a, radius_b, delta_ab,
x_ab, x_ab_length, kappa);
}
t@@ -192,12 +192,12 @@ __device__ void findAndProcessContactsInCell(
// Kernel executed on device, and callable from device only.
// Function is called from topology().
__device__ void findContactsInCell(
- int3 targetCell,
- const unsigned int idx_a,
+ int3 targetCell,
+ const unsigned int idx_a,
const Float4 x_a,
const Float radius_a,
- const Float4* __restrict__ dev_x_sorted,
- const unsigned int* __restrict__ dev_cellStart,
+ const Float4* __restrict__ dev_x_sorted,
+ const unsigned int* __restrict__ dev_cellStart,
const unsigned int* __restrict__ dev_cellEnd,
const unsigned int* __restrict__ dev_gridParticleIndex,
int* nc,
t@@ -215,7 +215,7 @@ __device__ void findContactsInCell(
// Calculate linear cell ID
unsigned int cellID = targetCell.x + targetCell.y * devC_grid.num[0]
- + (devC_grid.num[0] * devC_grid.num[1]) * targetCell.z;
+ + (devC_grid.num[0] * devC_grid.num[1]) * targetCell.z;
// Lowest particle index in cell
unsigned int startIdx = dev_cellStart[cellID];
t@@ -246,8 +246,8 @@ __device__ void findContactsInCell(
//__syncthreads();
// Distance between particle centers (Float4 -> Float3)
- Float3 x_ab = MAKE_FLOAT3(x_a.x - x_b.x,
- x_a.y - x_b.y,
+ Float3 x_ab = MAKE_FLOAT3(x_a.x - x_b.x,
+ x_a.y - x_b.y,
x_a.z - x_b.z);
// Adjust interparticle vector if periodic boundary/boundaries
t@@ -257,7 +257,7 @@ __device__ void findContactsInCell(
Float x_ab_length = length(x_ab);
// Distance between particle perimeters
- Float delta_ab = x_ab_length - (radius_a + radius_b);
+ Float delta_ab = x_ab_length - (radius_a + radius_b);
// Check for particle overlap
if (delta_ab < 0.0f) {
t@@ -323,11 +323,11 @@ __device__ void findContactsInCell(
// For a single particle:
-// Search for neighbors to particle 'idx' inside the 27 closest cells,
+// Search for neighbors to particle 'idx' inside the 27 closest cells,
// and save the contact pairs in global memory.
// Function is called from mainGPU loop.
__global__ void topology(
- const unsigned int* __restrict__ dev_cellStart,
+ const unsigned int* __restrict__ dev_cellStart,
const unsigned int* __restrict__ dev_cellEnd,
const unsigned int* __restrict__ dev_gridParticleIndex,
const Float4* __restrict__ dev_x_sorted,
t@@ -353,14 +353,17 @@ __global__ void topology(
gridPos.y = floor((x_a.y - devC_grid.origo[1]) / (devC_grid.L[1]/devC_grid.num[1]));
gridPos.z = floor((x_a.z - devC_grid.origo[2]) / (devC_grid.L[2]/devC_grid.num[2]));
+ printf("\ntopology: idx_a = %d, gridpos = [%d, %d, %d]\n", idx_a, gridPos.x, gridPos.y, gridPos.z);
+
// Find overlaps between particle no. idx and all particles
// from its own cell + 26 neighbor cells
for (int z_dim=-1; z_dim<2; ++z_dim) { // z-axis
for (int y_dim=-1; y_dim<2; ++y_dim) { // y-axis
for (int x_dim=-1; x_dim<2; ++x_dim) { // x-axis
targetPos = gridPos + make_int3(x_dim, y_dim, z_dim);
+ printf("\n targetPos = [%d, %d, %d]\n", targetPos.x, targetPos.y, targetPos.z);
findContactsInCell(targetPos, idx_a, x_a, radius_a,
- dev_x_sorted,
+ dev_x_sorted,
dev_cellStart, dev_cellEnd,
dev_gridParticleIndex,
&nc, dev_contacts, dev_distmod);
t@@ -373,7 +376,7 @@ __global__ void topology(
// For a single particle:
// If contactmodel=1:
-// Search for neighbors to particle 'idx' inside the 27 closest cells,
+// Search for neighbors to particle 'idx' inside the 27 closest cells,
// and compute the resulting force and torque on it.
// If contactmodel=2:
// Process contacts saved in dev_contacts by topology(), and compute
t@@ -417,11 +420,11 @@ __global__ void interact(
Float radius_a = x_a.w;
// Fetch world dimensions in constant memory read
- Float3 origo = MAKE_FLOAT3(devC_grid.origo[0],
- devC_grid.origo[1],
- devC_grid.origo[2]);
- Float3 L = MAKE_FLOAT3(devC_grid.L[0],
- devC_grid.L[1],
+ Float3 origo = MAKE_FLOAT3(devC_grid.origo[0],
+ devC_grid.origo[1],
+ devC_grid.origo[2]);
+ Float3 L = MAKE_FLOAT3(devC_grid.L[0],
+ devC_grid.L[1],
devC_grid.L[2]);
// Fetch wall data in global read
t@@ -470,7 +473,7 @@ __global__ void interact(
//uint4 bonds = dev_bonds_sorted[idx_a];
// Initiate shear friction loss rate at 0.0
- Float es_dot = 0.0;
+ Float es_dot = 0.0;
Float ev_dot = 0.0;
// Initiate pressure on particle at 0.0
t@@ -518,28 +521,28 @@ __global__ void interact(
// Process collision if the particles are overlapping
if (delta_n < 0.0) {
if (devC_params.contactmodel == 2) {
- contactLinear(&F, &T, &es_dot, &ev_dot, &p,
+ contactLinear(&F, &T, &es_dot, &ev_dot, &p,
idx_a_orig,
idx_b_orig,
vel_a,
dev_vel,
angvel_a,
dev_angvel,
- radius_a, radius_b,
+ radius_a, radius_b,
x_ab, x_ab_length,
- delta_n, dev_delta_t,
+ delta_n, dev_delta_t,
mempos);
} else if (devC_params.contactmodel == 3) {
- contactHertz(&F, &T, &es_dot, &ev_dot, &p,
+ contactHertz(&F, &T, &es_dot, &ev_dot, &p,
idx_a_orig,
idx_b_orig,
vel_a,
dev_vel,
angvel_a,
dev_angvel,
- radius_a, radius_b,
+ radius_a, radius_b,
x_ab, x_ab_length,
- delta_n, dev_delta_t,
+ delta_n, dev_delta_t,
mempos);
}
} else {
t@@ -555,7 +558,7 @@ __global__ void interact(
__syncthreads();
// Zero sum of shear displacement in this position
dev_delta_t[mempos] = MAKE_FLOAT4(0.0, 0.0, 0.0, 0.0);
- }
+ }
} // Contact loop end