From dac9dce752d0f011452ddd4b336d39b6729f4d7e Mon Sep 17 00:00:00 2001 From: DanNegrut Date: Fri, 4 Sep 2026 12:03:57 -0500 Subject: [PATCH] dT: check the integrated state is finite, covering angular velocity as well kT has two NaN tests already, and neither can see angular velocity. That is not a placement mistake: angular velocity never crosses to kT. kT is handed the translational speed array and the orientation quaternions, so no check placed in kT can cover this at any strength. Hence a check on dT. The two existing tests, stated precisely so this change is not mistaken for a duplicate of either: 1. unpackMyBuffer() tests isfinite() on a max-reduced velocity. This one is weak even on the translational side, because cub::DeviceReduce::Max compares with >, every comparison against NaN is false, so NaN operands are skipped and an all-NaN array reduces to the initial value -FLT_MAX, which is finite and passes. Measured: a run whose whole state was NaN reported "kinetic energy is nan while max speed is -3.40282e+38". 2. computeMarginFromAbsv() does a per-element isfinite() with DEME_ABORT_KERNEL. This one is strong, and it is the real guard on the translational side. But it reads marginSize, i.e. translational speed, and it is skipped entirely when isExpandFactorFixed, since the margin is then not derived from speed at all. Angular velocity is where the corruption arrives first, which is what makes the gap matter rather than merely exist. The defect fixed by #72 and #74 produced NaN through cross(cntPnt, F) with F exactly zero, and that lands in angular acceleration via atomicAdd on alphaX (DEMCollectForceKernels_Compact.cu:51 and :97), which integrates into angular velocity, then through the orientation quaternion into positions, contact geometry, forces, and only last into translational velocity. Measured on a pre-fix tree with this check in place, 40 repetitions, one Blackwell GPU, CUDA 12.9, one process per repetition, machine otherwise idle. The reproducer sets no user expand factor, so the automatic margin route was active and test 2 above was on duty throughout: 36 ran clean 3 this check reported angular velocity (x) as first non-finite 0 this check reported translational velocity as first non-finite 1 kT wedged in the contact-detection livelock before dT got another turn Test 2 did not abort in any of those 40 runs, nor in 20 further sanitized runs on current main. It was awake and the failures went past it, because all of them were angular. In each of the three reports the last completed step was healthy, at 0.570, 0.746 and 0.589 m/s with finite kinetic energy, so the abort lands on the first step where the state goes bad rather than once the corruption has surfaced. All four affected runs failed inside the contact-array growth window, at 4366 to 9043 contacts. checkStateIsFinite() probes the translational speed array and the three angular components with sum reductions, which propagate NaN where a max does not, and names which array failed. cubSumReduce is already instantiated, so no new template is introduced. Cost is at most four device reductions per contact-detection cycle rather than per time step, short-circuiting on the first bad array. Verified not to fire on healthy runs: 40 repetitions on current main, zero firings and zero failures. That property matters because the check aborts. An injected NaN in angular velocity is reported by name, where without this change the same injection dies with an unspecified launch failure out of DEMCubContactDetection. An injected NaN in translational velocity is not improved by this change and does not need to be, since test 2 catches it first. Known limitation: if kT wedges in the contact-detection livelock first, dT waits on it and this check is starved rather than blind. That is the one repetition above. The abort on non-finite bin coordinates added at DEMBinSphereKernels.cu:63 covers that path from the kT side; the two are complementary. Not done here, deliberately: giving cubDEMMax a NaN-propagating functor would also fix the clump_max_absv inspector, which still returns a finite number for an all-NaN system, but that changes public inspector behaviour. Co-Authored-By: Claude Opus 5 (1M context) --- src/DEM/dT.cpp | 50 ++++++++++++++++++++++++++++++++++++++++++++++++++ src/DEM/dT.h | 6 ++++++ 2 files changed, 56 insertions(+) diff --git a/src/DEM/dT.cpp b/src/DEM/dT.cpp index 86e272a3..1fe9cad9 100644 --- a/src/DEM/dT.cpp +++ b/src/DEM/dT.cpp @@ -7,6 +7,7 @@ #include #include #include +#include #ifdef DEME_USE_CHPF #include @@ -2289,10 +2290,59 @@ inline void DEMDynamicThread::ifProduceFreshThenUseIt() { } } + +inline void DEMDynamicThread::checkStateIsFinite() { + // WHY THIS EXISTS, AND WHY IT IS HERE RATHER THAN IN kT. + // kT has two NaN tests already, and neither one can see angular velocity, because angular + // velocity never crosses to kT: it lives here on dT. kT is handed the translational speed + // array and the orientation quaternions, so no check placed in kT can cover this at any + // strength. That is the whole reason this function is on dT. + // The two kT tests, for the record, so nobody re-adds a third one there: + // 1. unpackMyBuffer() tests isfinite() on a max-reduced velocity. Weak even on the + // translational side: cub::DeviceReduce::Max compares with >, every comparison against + // NaN is false, so NaN operands are skipped and an all-NaN array reduces to the initial + // value -FLT_MAX, which is finite and passes. + // 2. computeMarginFromAbsv() does a per-element isfinite() with DEME_ABORT_KERNEL, which + // is strong, but it reads marginSize, i.e. translational speed. It is also skipped + // entirely when isExpandFactorFixed, since then the margin is not derived from speed. + // Measured on a settling bed with the pre-fix tree, 40 runs, automatic margin route active + // so test 2 was on duty throughout: three real corruptions, all three naming angular + // velocity (x) here, with translational speed still a healthy 0.570, 0.746 and 0.589 m/s and + // kinetic energy finite. Test 2 did not fire in any of those 40 runs, nor in 20 more. + // A SUM propagates NaN where a MAX does not, so probing each array with one is enough. The + // arrays are checked separately so the message can say which part of the state went bad. + const size_t n = simParams->nOwnerBodies; + const char* which = nullptr; + float* const arrays[4] = {pCycleVel, granData->omgBarX, granData->omgBarY, granData->omgBarZ}; + const char* const names[4] = {"velocity", "angular velocity (x)", "angular velocity (y)", + "angular velocity (z)"}; + for (int i = 0; i < 4; i++) { + if (!arrays[i]) + continue; + cubSumReduce(arrays[i], &(stateFiniteProbe), n, streamInfo.stream, solverScratchSpace); + stateFiniteProbe.toHost(); + if (!std::isfinite(*(stateFiniteProbe))) { + which = names[i]; + break; + } + } + if (which) { + DEME_ERROR( + "The system state is not finite: the %s array contains NaN or Inf.\nThis usually means the " + "simulation diverged. Decreasing the step size, softening material properties, or checking for " + "elements initialized inside walls may help.\nNote that a reported max velocity can look " + "entirely reasonable while this fires, because a max reduction cannot propagate NaN, and because " + "angular velocity is not part of that reduction at all.", + which); + } +} inline void DEMDynamicThread::calibrateParams() { // Unpacking is done; now we can use temp arrays again to derive max velocity and send to kT pCycleVel = determineSysVel(); + // Refuse to hand kT a corrupted state, and say so, rather than continuing silently. + checkStateIsFinite(); + if (solverFlags.autoUpdateFreq) { unsigned int comfortable_drift; if (accumStepUpdater.Query(comfortable_drift)) { diff --git a/src/DEM/dT.h b/src/DEM/dT.h index c5fa9851..04e58243 100644 --- a/src/DEM/dT.h +++ b/src/DEM/dT.h @@ -81,6 +81,10 @@ class DEMDynamicThread { // dT believes this amount of future drift is ideal DualStruct perhapsIdealFutureDrift = DualStruct(0); + // Finiteness probes for the system state, one per array checked. SUM reductions are used + // rather than MAX because a max reduction cannot propagate NaN: IEEE comparisons against + // NaN are false, so cub::DeviceReduce::Max silently skips NaN operands. + DualStruct stateFiniteProbe = DualStruct(0.f); // Buffer arrays for storing info from the dT side. // kT modifies these arrays; dT uses them only. @@ -693,6 +697,8 @@ class DEMDynamicThread { // Change sim params based on dT's experience, if needed inline void calibrateParams(); + // Verifies the integrated state is finite, and reports which component is not. + inline void checkStateIsFinite(); // Determine the max vel for this cycle, kT needs it inline float* determineSysVel();