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

[Ward-Vandepitte]

Body

Summary

The stochastic "kT wedged at 100% GPU" livelock reported in #71 is not a
synchronization bug. It is a finite work explosion in contact detection, triggered by
NaN owner quaternions, which in turn originate from uninitialized device memory in
the contact recording arrays. This PR fixes the root cause and adds two defensive
layers.

Causal chain (each link verified with env-gated host-side probes on captured

wedges; ~90 instrumented repro runs)

1. During contact-count growth, contactEventArraysResize resizes
   contactForces / contactTorque_convToForce / contactPointGeometryA/B with
   DEME_DUAL_ARRAY_RESIZE(..., make_float3(0)). DualArray::resize(n, val) fills
   host values only (its own comment says so); resizeDevice raw-allocates and
   copies the old prefix — the grown device tail is recycled, uninitialized memory.
2. prepareForceArrays clears the two force arrays every step, but
   contactPointGeometryA/B are never cleared, and calculateContactForces writes
   them only for ContactType != NOT_A_CONTACT. Margin-only pairs (the majority of
   the pair list, by design of the CD margin) therefore keep garbage contact points
   for their entire lifetime.
3. forceToAngAcc computes cross(cntPnt, (F + torque_inForceForm) * mod) / MOI
   over ALL contact entries. For a margin-only entry the force is zero, and
   cross(finite_garbage, 0) = 0 — silently fine. But cross(NaN/inf, 0) = NaN,
   so one NaN-holding persistent margin pair injects NaN into its owner's angular
   acceleration every step. The linear collect (acc = F/mass) never reads
   cntPnt, which is why only the angular channel dies (observed: velocities stay
   finite through the whole event, so SetErrorOutVelocity never fires).
4. NaN alpha -> omgBar -> quaternion (the integrator multiplies by
   (0.5*h*omgBar, 1) and normalizes) -> NaN sphere positions. The clamp in
   getNumberOfBinsEachSphereTouches is NaN-blind (every ternary comparison is
   false -> lo=0, hi=nb-1 in all three dims), so each NaN sphere claims the ENTIRE
   bin grid. Measured at the wedge: 73,728,000 sphere-bin pairs = 1800 spheres x
   40960 bins exactly; maxGeoInBin = all spheres.
5. The per-bin n(n-1)/2 sweep (getNumberOfSphereContactsEachBin) then faces
   ~6.6e10 pair tests in one launch. One captured instance completed in 56.66 s and
   the sibling populate kernel hung next; ordinarily the launch simply never returns
   within any reasonable watchdog window. GPU pegged at 100%, kT inside one kernel:
   the reported livelock signature.

Notes that may help maintainers reconcile this with prior reports: the failure is
concentrated at the contact ramp because that is when the arrays grow; it is
stochastic because it depends on the recycled memory content (usually finite garbage
= a silent zero-effect; occasionally NaN/inf); and it is config-independent — a
pre-registered campaign falsified tuner on/off, drift clamps, fixed expand factor,
CD frequency (only CDUpdateFreq=1 suppressed the rate, by shrinking the margin pair
population), and a CUDA 12.6/12.8 toolkit swap. When the garbage lands in the force
tail instead, the run dies loudly at SetErrorOutVelocity — likely some of the
"diverged for no reason" reports are the same defect's other face.

Fix (three layers)

• (a) Root: contactEventArraysResize AND allocateGPUArrays now cudaMemset the
  device ranges of contactPointGeometryA/B after resizing (the second site is
  reachable from every Initialize() and mid-run UpdateClumps(), and its garbage
  is also user-visible through GetContactDetailedInfo() default output). The two
  force arrays need no such treatment: prepareForceArrays clears them every step.
  The deeper hazard is that DualArray::resize(n, val) applies val to host
  memory only -- happy to follow up with a patch there instead/as well if you
  prefer that semantics change.
• (b) Defense: forceToAngAcc returns zero early when F + torque_inForceForm is
  exactly zero, so never-written cntPnt slots are unconditionally harmless. These
  entries are the majority of the list, so this also skips their MOI and quaternion
  loads (and removes a silent finite-garbage torque perturbation that nobody would
  otherwise notice).
• (c) Belt: getNumberOfBinsEachSphereTouches aborts loudly (same style as the
  errOutBinSphNum check) on a non-finite position or radius span instead of
  silently claiming the whole grid. This is an INTENTIONAL new fail-fast: a
  non-finite sphere position is never a legitimate simulation state, so no opt-out
  is provided -- any future NaN state becomes a clean error rather than a phantom
  hang.

Validation

On the #71 reproduction scene (600 clumps / 1800 spheres, frictional Hertzian,
stock settings; stochastic livelock rate 27%+ per truncated run, 95% CI LB 0.269
over the campaign):

• Fixed build: 17/17 truncated runs clean (p ~ 0.005 against the stock rate) and
  4/4 full-length runs clean.
• Physics pins unchanged: 5/5 (bridge stiffness, normal stiffness, Coulomb
  friction, restitution, wall friction) pass on the fixed build.
• 24/24 production-shaped settling runs clean under a no-retry watchdog harness.
• Post-review final shape (this exact PR content, rebuilt): physics pins 5/5,
  17/17 truncated runs clean (arm prshape_final, 2026-07-22).