Fix stochastic contact-ramp livelock: uninitialized contact-point device memory NaNs the angular channel and degenerates CD into an all-pairs sweep (fixes #71)

src/DEM/dT.cpp

@@ -532,6 +532,12 @@ void DEMDynamicThread::allocateGPUArrays(size_t nOwnerBodies,
             DEME_DUAL_ARRAY_RESIZE(contactTorque_convToForce, cnt_arr_size, make_float3(0));
             DEME_DUAL_ARRAY_RESIZE(contactPointGeometryA, cnt_arr_size, make_float3(0));
             DEME_DUAL_ARRAY_RESIZE(contactPointGeometryB, cnt_arr_size, make_float3(0));
+            // The contact point device ranges must be zeroed for the same reason as in
+            // contactEventArraysResize: resize fills host values only, and these two arrays
+            // are never per-step cleared, so entries that no kernel writes (margin-only
+            // pairs) would otherwise expose recycled device memory.
+            DEME_GPU_CALL(cudaMemset(contactPointGeometryA.data(), 0, cnt_arr_size * sizeof(float3)));
+            DEME_GPU_CALL(cudaMemset(contactPointGeometryB.data(), 0, cnt_arr_size * sizeof(float3)));
         }
         // Allocate memory for each wildcard array
         contactWildcards.resize(simParams->nContactWildcards);
@@ -1945,6 +1951,16 @@ inline void DEMDynamicThread::contactEventArraysResize(size_t nContactPairs) {
         DEME_DUAL_ARRAY_RESIZE(contactTorque_convToForce, nContactPairs, make_float3(0));
         DEME_DUAL_ARRAY_RESIZE(contactPointGeometryA, nContactPairs, make_float3(0));
         DEME_DUAL_ARRAY_RESIZE(contactPointGeometryB, nContactPairs, make_float3(0));
+        // DualArray::resize(n, val) fills host values only, so the grown device range holds
+        // whatever the allocator recycled. contactForces and contactTorque_convToForce are
+        // cleared every time step by prepareForceArrays before use, but the contact point
+        // arrays are only ever written for entries whose ContactType is a real contact, and
+        // collectContactForcesAccStyle reads them for ALL entries: cross(garbage, 0) is 0 for
+        // finite garbage but NaN for NaN/inf garbage, which poisons that owner's angular
+        // acceleration and cascades (quaternion -> position -> full-grid bin ranges in kT's
+        // contact detection, presenting as a GPU-pegged stall). Zero them on device here.
+        DEME_GPU_CALL(cudaMemset(contactPointGeometryA.data(), 0, nContactPairs * sizeof(float3)));
+        DEME_GPU_CALL(cudaMemset(contactPointGeometryB.data(), 0, nContactPairs * sizeof(float3)));
     }
 
     // Re-packing pointers now is automatic

src/kernel/DEMBinSphereKernels.cu

@@ -56,6 +56,18 @@ __global__ void getNumberOfBinsEachSphereTouches(deme::DEMSimParams* simParams,
                 double myBinZ = myPosXYZ.z / simParams->binSize;
                 // How many bins my radius spans (with fractions)?
                 double myRadiusSpan = myRadius / simParams->binSize;
+                // A non-finite position or span makes every comparison in the range clamping
+                // below evaluate false, which silently assigns this sphere the ENTIRE bin grid
+                // and turns the per-bin contact sweep into an all-pairs workload that presents
+                // as a GPU-pegged stall. Erroring out loudly is the only sane option here.
+                if (!isfinite(myBinX) || !isfinite(myBinY) || !isfinite(myBinZ) || !isfinite(myRadiusSpan)) {
+                    DEME_ABORT_KERNEL(
+                        "Sphere %u has a non-finite position or contact margin (bin coordinates %f, %f, %f, radius "
+                        "span %f).\nThis usually means the simulation has diverged (a NaN orientation quaternion or "
+                        "angular velocity is a common cause), and relaxing the physics may help, such as decreasing "
+                        "the step size and modifying material properties.\n",
+                        sphereID, myBinX, myBinY, myBinZ, myRadiusSpan);
+                }
                 // printf("myRadius: %f\n", myRadiusSpan);
                 // Now, figure out how many bins I touch in each direction
                 numX = ((myBinX + myRadiusSpan < (double)simParams->nbX) ? (unsigned int)(myBinX + myRadiusSpan)

src/kernel/DEMCollectForceKernels.cu

@@ -102,6 +102,19 @@ __global__ void forceToAngAcc(float3* angAcc,
                               deme::DEMDataDT* granData) {
     deme::contactPairs_t myID = blockIdx.x * blockDim.x + threadIdx.x;
     if (myID < n) {
+        // torque_inForceForm is usually the contribution of rolling resistance and it contributes to torque only, not
+        // linear velocity
+        float3 myF = F[myID] + torque_inForceForm[myID];
+        // Zero-force entries must contribute exactly zero angular acceleration without touching
+        // cntPnt: entries the force kernel never writes (pairs registered by kT's contact margin
+        // that have no physical overlap) can hold uninitialized device memory in their cntPnt
+        // slot, and cross(NaN, 0) is NaN, which would poison this owner's angular acceleration.
+        // These entries are the majority of the list, so the early-out also skips their MOI and
+        // quaternion loads.
+        if (myF.x == 0.0f && myF.y == 0.0f && myF.z == 0.0f) {
+            angAcc[myID] = make_float3(0, 0, 0);
+            return;
+        }
         const deme::bodyID_t myOwner = owner[myID];
         float3 myMOI;
         // Get my mass info from either jitified arrays or global memory
@@ -114,9 +127,7 @@ __global__ void forceToAngAcc(float3* angAcc,
         const deme::oriQ_t myOriQz = oriQz[myOwner];
 
         float3 myCntPnt = cntPnt[myID];
-        // torque_inForceForm is usually the contribution of rolling resistance and it contributes to torque only, not
-        // linear velocity
-        float3 myF = (F[myID] + torque_inForceForm[myID]) * modifier;
+        myF = myF * modifier;
         // F is in global frame, but it needs to be in local to coordinate with moi and cntPnt
         applyOriQToVector3<float, deme::oriQ_t>(myF.x, myF.y, myF.z, myOriQw, -myOriQx, -myOriQy, -myOriQz);
         angAcc[myID] = cross(myCntPnt, myF) / myMOI;