Performance Update Mesh - Mesh

- kT much faster via Fp32 (fp64 fallback)
- dT early memory improvments

Author: fl0rianr

src/DEM/dT.cpp

@@ -514,7 +514,7 @@ void DEMDynamicThread::changeOwnerSizes(const std::vector<bodyID_t>& IDs, const
     size_t ownerFactorSize = (size_t)simParams->nOwnerBodies * sizeof(float);
     // Bool table for whether this owner should change
     notStupidBool_t* idBool = (notStupidBool_t*)solverScratchSpace.allocateTempVector("idBool", idBoolSize);
-    DEME_GPU_CALL(cudaMemset(idBool, 0, idBoolSize));
+    DEME_GPU_CALL(cudaMemsetAsync(idBool, 0, idBoolSize, streamInfo.stream));
     float* ownerFactors = (float*)solverScratchSpace.allocateTempVector("ownerFactors", ownerFactorSize);
 
     // Mark on the bool array those owners that need a change
@@ -632,7 +632,7 @@ void DEMDynamicThread::allocateGPUArrays(size_t nOwnerBodies,
     // maxTriTriPenetration usually keeps the max tri--tri penetration during the on-going simulation. But after
     // initialization, when it stores no meaningful values, dT will send a work order to kT, so maxTriTriPenetration's
     // value has to be initialized.
-    DEME_GPU_CALL(cudaMemset(maxTriTriPenetration.getDevicePointer(), 0, sizeof(double)));
+    DEME_GPU_CALL(cudaMemsetAsync(maxTriTriPenetration.getDevicePointer(), 0, sizeof(double), streamInfo.stream));
 
     // Resize to the number of analytical geometries
     DEME_DUAL_ARRAY_RESIZE(ownerAnalBody, nAnalGM, 0);
@@ -2495,8 +2495,26 @@ inline void DEMDynamicThread::dispatchPatchBasedForceCorrections(
     const ContactTypeMap<std::pair<contactPairs_t, contactPairs_t>>& typeStartCountPrimitiveMap,
     const ContactTypeMap<std::pair<contactPairs_t, contactPairs_t>>& typeStartCountPatchMap,
     const ContactTypeMap<std::vector<std::pair<std::shared_ptr<JitHelper::CachedProgram>, std::string>>>& typeKernelMap) {
-    // Reset max tri-tri penetration for this timestep on device (kT may need this info)
-    DEME_GPU_CALL(cudaMemset(maxTriTriPenetration.getDevicePointer(), 0, sizeof(double)));
+    // Reset max tri-tri penetration only when mesh contacts are present (avoid needless clears), async to overlap
+    bool needTriPenReset = false;
+    auto markIfTri = [&](contact_t c) {
+        auto itP = typeStartCountPrimitiveMap.find(c);
+        if (itP != typeStartCountPrimitiveMap.end() && itP->second.second > 0) {
+            needTriPenReset = true;
+            return;
+        }
+        auto itPatch = typeStartCountPatchMap.find(c);
+        if (itPatch != typeStartCountPatchMap.end() && itPatch->second.second > 0) {
+            needTriPenReset = true;
+        }
+    };
+    markIfTri(SPHERE_TRIANGLE_CONTACT);
+    markIfTri(TRIANGLE_TRIANGLE_CONTACT);
+    markIfTri(TRIANGLE_ANALYTICAL_CONTACT);
+    if (needTriPenReset) {
+        DEME_GPU_CALL(
+            cudaMemsetAsync(maxTriTriPenetration.getDevicePointer(), 0, sizeof(double), streamInfo.stream));
+    }
 
     // For each contact type that exists, check if it is patch(mesh)-related type...
     for (size_t i = 0; i < m_numExistingTypes; i++) {
@@ -3375,7 +3393,7 @@ float* DEMDynamicThread::inspectCall(const std::shared_ptr<JitHelper::CachedProg
     // If this boolArrExclude is 1 at an element, that means this element is exluded in the reduction
     notStupidBool_t* boolArrExclude =
         (notStupidBool_t*)solverScratchSpace.allocateTempVector("boolArrExclude", regionTempSize);
-    DEME_GPU_CALL(cudaMemset(boolArrExclude, 0, regionTempSize));
+    DEME_GPU_CALL(cudaMemsetAsync(boolArrExclude, 0, regionTempSize, streamInfo.stream));
 
     // We may actually have 2 reduced returns: in regional reduction, key 0 and 1 give one return each.
     size_t returnSize = sizeof(float) * 2;

src/DEM/kT.cpp

@@ -608,7 +608,7 @@ void DEMKinematicThread::changeOwnerSizes(const std::vector<bodyID_t>& IDs, cons
     size_t ownerFactorSize = (size_t)simParams->nOwnerBodies * sizeof(float);
     // Bool table for whether this owner should change
     notStupidBool_t* idBool = (notStupidBool_t*)solverScratchSpace.allocateTempVector("idBool", idBoolSize);
-    DEME_GPU_CALL(cudaMemset(idBool, 0, idBoolSize));
+    DEME_GPU_CALL(cudaMemsetAsync(idBool, 0, idBoolSize, streamInfo.stream));
     float* ownerFactors = (float*)solverScratchSpace.allocateTempVector("ownerFactors", ownerFactorSize);
 
     // Mark on the bool array those owners that need a change

src/kernel/DEMBinTriangleKernels.cu

@@ -56,7 +56,62 @@ __global__ void makeTriangleSandwich(deme::DEMSimParams* simParams,
     }
 }
 
-inline __device__ void figureOutNodeAndBoundingBox(deme::DEMSimParams* simParams,
+
+// Compute triangle AABB -> bin index bounds using mixed precision (FP32 fast path + FP64 fallback).
+// This mirrors the sphere binning approach (axis_bounds) to avoid precision issues when a bound lies
+// close to a bin boundary.
+inline __device__ bool boundingBoxIntersectBinAxisBounds(deme::binID_t* L,
+                                                         deme::binID_t* U,
+                                                         const float3& vA,
+                                                         const float3& vB,
+                                                         const float3& vC,
+                                                         deme::DEMSimParams* simParams) {
+    float3 min_pt;
+    min_pt.x = DEME_MIN(vA.x, DEME_MIN(vB.x, vC.x));
+    min_pt.y = DEME_MIN(vA.y, DEME_MIN(vB.y, vC.y));
+    min_pt.z = DEME_MIN(vA.z, DEME_MIN(vB.z, vC.z));
+
+    float3 max_pt;
+    max_pt.x = DEME_MAX(vA.x, DEME_MAX(vB.x, vC.x));
+    max_pt.y = DEME_MAX(vA.y, DEME_MAX(vB.y, vC.y));
+    max_pt.z = DEME_MAX(vA.z, DEME_MAX(vB.z, vC.z));
+
+    // Enlarge bounding box, so that no triangle lies right between 2 layers of bins
+    const float enlarge = (float)DEME_BIN_ENLARGE_RATIO_FOR_FACETS * (float)simParams->dyn.binSize;
+    min_pt -= enlarge;
+    max_pt += enlarge;
+
+    const double invBinSize = simParams->dyn.inv_binSize;
+    const int nbX = (int)simParams->nbX;
+    const int nbY = (int)simParams->nbY;
+    const int nbZ = (int)simParams->nbZ;
+
+    // Convert [min,max] to (center, half-range) and use axis_bounds (FP32 fast path with FP64 fallback).
+    const double cx = 0.5 * ((double)min_pt.x + (double)max_pt.x);
+    const double rx = 0.5 * ((double)max_pt.x - (double)min_pt.x);
+    const deme::AxisBounds bx = axis_bounds(cx, rx, nbX, invBinSize);
+    if (bx.imax < bx.imin) return false;
+
+    const double cy = 0.5 * ((double)min_pt.y + (double)max_pt.y);
+    const double ry = 0.5 * ((double)max_pt.y - (double)min_pt.y);
+    const deme::AxisBounds by = axis_bounds(cy, ry, nbY, invBinSize);
+    if (by.imax < by.imin) return false;
+
+    const double cz = 0.5 * ((double)min_pt.z + (double)max_pt.z);
+    const double rz = 0.5 * ((double)max_pt.z - (double)min_pt.z);
+    const deme::AxisBounds bz = axis_bounds(cz, rz, nbZ, invBinSize);
+    if (bz.imax < bz.imin) return false;
+
+    L[0] = (deme::binID_t)bx.imin;
+    U[0] = (deme::binID_t)bx.imax;
+    L[1] = (deme::binID_t)by.imin;
+    U[1] = (deme::binID_t)by.imax;
+    L[2] = (deme::binID_t)bz.imin;
+    U[2] = (deme::binID_t)bz.imax;
+    return true;
+}
+
+inline __device__ bool figureOutNodeAndBoundingBox(deme::DEMSimParams* simParams,
                                                    deme::DEMDataKT* granData,
                                                    const deme::bodyID_t& triID,
                                                    float3& vA,
@@ -85,7 +140,7 @@ inline __device__ void figureOutNodeAndBoundingBox(deme::DEMSimParams* simParams
     vB = ownerXYZ + loc_vB;
     vC = ownerXYZ + loc_vC;
 
-    boundingBoxIntersectBin(L, U, vA, vB, vC, simParams);
+    return boundingBoxIntersectBinAxisBounds(L, U, vA, vB, vC, simParams);
 }
 
 __global__ void getNumberOfBinsEachTriangleTouches(deme::DEMSimParams* simParams,
@@ -106,38 +161,64 @@ __global__ void getNumberOfBinsEachTriangleTouches(deme::DEMSimParams* simParams
         // bin-based locations don't need that)
         float3 vA1, vB1, vC1, vA2, vB2, vC2;
         deme::binID_t L1[3], L2[3], U1[3], U2[3];
-        figureOutNodeAndBoundingBox(simParams, granData, triID, vA1, vB1, vC1, L1, U1, nodeA1[triID], nodeB1[triID],
-                                    nodeC1[triID]);
-        figureOutNodeAndBoundingBox(simParams, granData, triID, vA2, vB2, vC2, L2, U2, nodeA2[triID], nodeB2[triID],
-                                    nodeC2[triID]);
-        L1[0] = DEME_MIN(L1[0], L2[0]);
-        L1[1] = DEME_MIN(L1[1], L2[1]);
-        L1[2] = DEME_MIN(L1[2], L2[2]);
-        U1[0] = DEME_MAX(U1[0], U2[0]);
-        U1[1] = DEME_MAX(U1[1], U2[1]);
-        U1[2] = DEME_MAX(U1[2], U2[2]);
+        const bool ok1 = figureOutNodeAndBoundingBox(simParams, granData, triID, vA1, vB1, vC1, L1, U1, nodeA1[triID],
+                                                     nodeB1[triID], nodeC1[triID]);
+        const bool ok2 = figureOutNodeAndBoundingBox(simParams, granData, triID, vA2, vB2, vC2, L2, U2, nodeA2[triID],
+                                                     nodeB2[triID], nodeC2[triID]);
+
+        // If neither triangle sandwich intersects the bin grid, it cannot touch any bin.
+        if (!ok1 && !ok2) {
+            numBinsTriTouches[triID] = 0;
+            if (meshUniversalContact) {
+                numAnalGeoTriTouches[triID] = 0;
+            }
+            return;
+        }
+
+        // Merge bounds (or take the valid one, if only one is valid).
+        if (ok1 && ok2) {
+            L1[0] = DEME_MIN(L1[0], L2[0]);
+            L1[1] = DEME_MIN(L1[1], L2[1]);
+            L1[2] = DEME_MIN(L1[2], L2[2]);
+            U1[0] = DEME_MAX(U1[0], U2[0]);
+            U1[1] = DEME_MAX(U1[1], U2[1]);
+            U1[2] = DEME_MAX(U1[2], U2[2]);
+        } else if (!ok1) {
+            L1[0] = L2[0];
+            L1[1] = L2[1];
+            L1[2] = L2[2];
+            U1[0] = U2[0];
+            U1[1] = U2[1];
+            U1[2] = U2[2];
+        }
 
         unsigned int numSDsTouched = 0;
         // Triangle may span a collection of bins...
         // BTW, I don't know why Chrono::GPU had to check the so-called 3 cases, and create thread divergence like that.
         // Just sweep through all potential bins and you are fine.
         float BinCenter[3];
-        float BinHalfSizes[3];
-        BinHalfSizes[0] = simParams->dyn.binSize / 2. + DEME_BIN_ENLARGE_RATIO_FOR_FACETS * simParams->dyn.binSize;
-        BinHalfSizes[1] = simParams->dyn.binSize / 2. + DEME_BIN_ENLARGE_RATIO_FOR_FACETS * simParams->dyn.binSize;
-        BinHalfSizes[2] = simParams->dyn.binSize / 2. + DEME_BIN_ENLARGE_RATIO_FOR_FACETS * simParams->dyn.binSize;
-        for (deme::binID_t i = L1[0]; i <= U1[0]; i++) {
+        const float binSizeF = (float)simParams->dyn.binSize;
+        const float binHalfSpan = binSizeF * (0.5f + (float)DEME_BIN_ENLARGE_RATIO_FOR_FACETS);
+        float BinHalfSizes[3] = {binHalfSpan, binHalfSpan, binHalfSpan};
+        const float startX = binSizeF * (float)L1[0] + 0.5f * binSizeF;
+        const float startY = binSizeF * (float)L1[1] + 0.5f * binSizeF;
+        const float startZ = binSizeF * (float)L1[2] + 0.5f * binSizeF;
+        for (deme::binID_t i = L1[0], ix = 0; i <= U1[0]; i++, ix++) {
+            float cy0 = startY;
+            BinCenter[0] = startX + ix * binSizeF;
             for (deme::binID_t j = L1[1]; j <= U1[1]; j++) {
+                float cz = startZ;
+                BinCenter[1] = cy0;
                 for (deme::binID_t k = L1[2]; k <= U1[2]; k++) {
-                    BinCenter[0] = simParams->dyn.binSize * i + simParams->dyn.binSize / 2.;
-                    BinCenter[1] = simParams->dyn.binSize * j + simParams->dyn.binSize / 2.;
-                    BinCenter[2] = simParams->dyn.binSize * k + simParams->dyn.binSize / 2.;
+                    BinCenter[2] = cz;
 
                     if (check_TriangleBoxOverlap(BinCenter, BinHalfSizes, vA1, vB1, vC1) ||
                         check_TriangleBoxOverlap(BinCenter, BinHalfSizes, vA2, vB2, vC2)) {
                         numSDsTouched++;
                     }
+                    cz += binSizeF;
                 }
+                cy0 += binSizeF;
             }
         }
         numBinsTriTouches[triID] = numSDsTouched;
@@ -221,36 +302,56 @@ __global__ void populateBinTriangleTouchingPairs(deme::DEMSimParams* simParams,
         // 3 vertices of the triangle
         float3 vA1, vB1, vC1, vA2, vB2, vC2;
         deme::binID_t L1[3], L2[3], U1[3], U2[3];
-        figureOutNodeAndBoundingBox(simParams, granData, triID, vA1, vB1, vC1, L1, U1, nodeA1[triID], nodeB1[triID],
-                                    nodeC1[triID]);
-        figureOutNodeAndBoundingBox(simParams, granData, triID, vA2, vB2, vC2, L2, U2, nodeA2[triID], nodeB2[triID],
-                                    nodeC2[triID]);
-        L1[0] = DEME_MIN(L1[0], L2[0]);
-        L1[1] = DEME_MIN(L1[1], L2[1]);
-        L1[2] = DEME_MIN(L1[2], L2[2]);
-        U1[0] = DEME_MAX(U1[0], U2[0]);
-        U1[1] = DEME_MAX(U1[1], U2[1]);
-        U1[2] = DEME_MAX(U1[2], U2[2]);
+        const bool ok1 = figureOutNodeAndBoundingBox(simParams, granData, triID, vA1, vB1, vC1, L1, U1, nodeA1[triID],
+                                                     nodeB1[triID], nodeC1[triID]);
+        const bool ok2 = figureOutNodeAndBoundingBox(simParams, granData, triID, vA2, vB2, vC2, L2, U2, nodeA2[triID],
+                                                     nodeB2[triID], nodeC2[triID]);
+
+        // If neither triangle sandwich intersects the bin grid, it cannot touch any bin.
+        if (!ok1 && !ok2) {
+            return;
+        }
+
+        // Merge bounds (or take the valid one, if only one is valid).
+        if (ok1 && ok2) {
+            L1[0] = DEME_MIN(L1[0], L2[0]);
+            L1[1] = DEME_MIN(L1[1], L2[1]);
+            L1[2] = DEME_MIN(L1[2], L2[2]);
+            U1[0] = DEME_MAX(U1[0], U2[0]);
+            U1[1] = DEME_MAX(U1[1], U2[1]);
+            U1[2] = DEME_MAX(U1[2], U2[2]);
+        } else if (!ok1) {
+            L1[0] = L2[0];
+            L1[1] = L2[1];
+            L1[2] = L2[2];
+            U1[0] = U2[0];
+            U1[1] = U2[1];
+            U1[2] = U2[2];
+        }
 
         deme::binsTriangleTouchPairs_t myReportOffset = numBinsTriTouchesScan[triID];
         // In case this sweep does not agree with the previous one, we need to intercept such potential segfaults
         const deme::binsTriangleTouchPairs_t myReportOffset_end = numBinsTriTouchesScan[triID + 1];
 
         // Triangle may span a collection of bins...
         float BinCenter[3];
-        float BinHalfSizes[3];
-        BinHalfSizes[0] = simParams->dyn.binSize / 2. + DEME_BIN_ENLARGE_RATIO_FOR_FACETS * simParams->dyn.binSize;
-        BinHalfSizes[1] = simParams->dyn.binSize / 2. + DEME_BIN_ENLARGE_RATIO_FOR_FACETS * simParams->dyn.binSize;
-        BinHalfSizes[2] = simParams->dyn.binSize / 2. + DEME_BIN_ENLARGE_RATIO_FOR_FACETS * simParams->dyn.binSize;
-        for (deme::binID_t i = L1[0]; i <= U1[0]; i++) {
+        const float binSizeF = (float)simParams->dyn.binSize;
+        const float binHalfSpan = binSizeF * (0.5f + (float)DEME_BIN_ENLARGE_RATIO_FOR_FACETS);
+        float BinHalfSizes[3] = {binHalfSpan, binHalfSpan, binHalfSpan};
+        const float startX = binSizeF * (float)L1[0] + 0.5f * binSizeF;
+        const float startY = binSizeF * (float)L1[1] + 0.5f * binSizeF;
+        const float startZ = binSizeF * (float)L1[2] + 0.5f * binSizeF;
+        for (deme::binID_t i = L1[0], ix = 0; i <= U1[0]; i++, ix++) {
+            BinCenter[0] = startX + ix * binSizeF;
+            float cy0 = startY;
             for (deme::binID_t j = L1[1]; j <= U1[1]; j++) {
+                BinCenter[1] = cy0;
+                float cz = startZ;
                 for (deme::binID_t k = L1[2]; k <= U1[2]; k++) {
                     if (myReportOffset >= myReportOffset_end) {
                         continue;  // Don't step on the next triangle's domain
                     }
-                    BinCenter[0] = simParams->dyn.binSize * i + simParams->dyn.binSize / 2.;
-                    BinCenter[1] = simParams->dyn.binSize * j + simParams->dyn.binSize / 2.;
-                    BinCenter[2] = simParams->dyn.binSize * k + simParams->dyn.binSize / 2.;
+                    BinCenter[2] = cz;
 
                     if (check_TriangleBoxOverlap(BinCenter, BinHalfSizes, vA1, vB1, vC1) ||
                         check_TriangleBoxOverlap(BinCenter, BinHalfSizes, vA2, vB2, vC2)) {
@@ -259,7 +360,9 @@ __global__ void populateBinTriangleTouchingPairs(deme::DEMSimParams* simParams,
                         triIDsEachBinTouches[myReportOffset] = triID;
                         myReportOffset++;
                     }
+                    cz += binSizeF;
                 }
+                cy0 += binSizeF;
             }
         }
         // This can happen for like 1 in 10^9 chance, for the tri--bin contact algorithm has stochasticity on GPU