Mid-run livelock: DoDynamicsThenSync wedges (GPU 100%, no progress) during dense settle (bare-metal Linux, CUDA 12.8)

[Ward-Vandepitte]

Summary

On a single-GPU run, DoDynamicsThenSync intermittently never returns. The process stays alive with
the GPU at 99–100% SM utilization, but the simulation makes no further progress — no error, no CUDA
Xid, no exception. It is a livelock, not a crash.

Onset is stochastic and concentrated in the contact-ramp regime — the moment a gravity-settled
granular pack lands on the floor and the contact count ramps sharply. In a controlled repeat of a
minimal scene (below) it wedged in 3 of 7 runs, always at chunk ~30–40 (the ramp), where the
confirmed-contact count is climbing through ~2k → ~12k; the other 4 runs completed cleanly. It also
occasionally wedges much later (I have seen it at ~chunk 1449 in a longer run).

I originally hit this with a custom JIT force model and suspected my own kernel, but it reproduces
identically with the stock UseFrictionalHertzianModel() and no custom code, so it appears to be
in the async engine itself rather than in user force code.

Correction (updated): an earlier version of this report listed the OS as WSL2. That was wrong —
this reproduces on bare-metal Linux with a native (non-virtualized) GPU. See the Environment
section. Re-confirmed on bare-metal: rolling the minimal repro below, one run wedged at chunk
k=19 (nc≈2.5k, mid contact-ramp) with the GPU at 100% SM, one thread R (33% CPU) and 23
in futex_wait_queue — the exact signature below; a separate run of the same scene completed
3000/3000 chunks cleanly, i.e. the stochasticity holds off-WSL. This matters for the diagnosis: it
is not a WSL GPU-scheduling artifact, and it is therefore likely a different root cause from
#70 (which is genuinely on WSL2), or a shared root cause that is platform-independent.

Environment

• DEME 2.3.3 (PyPI sdist, built from source, CUDAARCHS=89)
• GPU: NVIDIA RTX 2000 Ada (compute 8.9), 16 GB
• Driver 580.159.03, CUDA toolkit 12.8, GCC 13.3
• OS: bare-metal Ubuntu 24.04 (native GPU — not a VM, not WSL; nvidia-smi -q reports
  Virtualization Mode: None, GPU on PCI 01:00.0), Python 3.12, DEMSolver(nGPUs=1)

Signature (how to recognize it)

• DoDynamicsThenSync does not return; process alive.
• nvidia-smi: the process at 99–100% SM, memory-controller ~0%.
• Thread states (ps -L -o stat,wchan): 1–2 threads R (running), the rest (~20+) in
  futex_wait_queue — the kinematic thread appears to spin while the dynamics thread and workers wait.
• faulthandler.dump_traceback_later(...) shows the main Python thread blocked inside
  DoDynamicsThenSync (a bare C-call frame; no Python-level loop).
• Distinct from a GPU-idle hang: here the GPU is busy.

Minimal reproduction

Gravity-settle ~600 three-sphere clumps (sub-mm grains, the bundled clumps/3_clump.csv scaled so
r_max ≈ 0.35 mm) into a 16 mm walled box, stock frictional-Hertzian model, dt = 2.66e-6, stepping in
chunks of 200 substeps. Because onset is stochastic, run it a handful of times — it wedged in ~3/7
here, typically within the first few tens of chunks (the contact ramp).

import faulthandler, os, time
import numpy as np
import DEME

RHO, GRAIN_MAX_R, FOOT, H, C, MU = 2650.0, 0.35e-3, 0.016, 0.05, 200, 0.50
YOUNG, POISSON, COR = 1.0e7, 0.30, 0.40
N_TARGET, MAX_CHUNKS, DT = 600, 400, 2.662e-6

faulthandler.enable()
faulthandler.dump_traceback_later(120.0, repeat=True)  # dumps thread stacks if wedged

rows = []
for line in open(DEME.GetDEMEDataFile("clumps/3_clump.csv")):
    line = line.strip()
    if not line or line.startswith("#") or line.lower().startswith("x"):
        continue
    rows.append([float(v) for v in line.split(",")[:4]])
comp0 = np.asarray(rows)
scale = GRAIN_MAX_R / comp0[:, 3].max()
comp = comp0 * scale
mass = RHO * 5.5886717 * scale**3
moi = [RHO * scale**5 * I for I in (1.8327927, 2.1580013, 0.77010059)]
bound_r = float(np.max(np.linalg.norm(comp[:, :3], axis=1) + comp[:, 3]))
sep = 2.2 * bound_r

tmp = "/tmp/deme_repro_clump.csv"
with open(tmp, "w") as f:
    f.write("x,y,z,r\n")
    for row in comp:
        f.write("%f,%f,%f,%f\n" % tuple(row[:4]))

s = DEME.DEMSolver(nGPUs=1)
s.SetVerbosity("WARNING")
mat = s.LoadMaterial({"E": YOUNG, "nu": POISSON, "CoR": COR, "mu": MU, "Crr": 0.0})
s.UseFrictionalHertzianModel()
clump = s.LoadClumpType(mass, moi, tmp, mat)

z_floor = -H / 2
half_xy = FOOT / 2 - sep
per_layer = max(1, int((2 * half_xy / sep + 1) ** 2))
n_layers = int(np.ceil(N_TARGET / per_layer)) + 2
col_h = n_layers * sep * np.sqrt(2.0 / 3.0)
z_bot = z_floor + 1.02 * bound_r
pts = DEME.HCPSampler(sep).SampleBox([0., 0., z_bot + col_h / 2.0],
                                     [half_xy, half_xy, col_h / 2.0])[:N_TARGET]
s.AddClumps(clump, pts)
s.InstructBoxDomainDimension(FOOT, FOOT, H)
s.InstructBoxDomainBoundingBC("top_open", mat)
s.SetGravitationalAcceleration([0., 0., -9.81])
s.SetInitTimeStep(DT)
s.SetMaxVelocity(0.5)
s.SetCDUpdateFreq(20)
s.Initialize()

print("SETUP done, n=%d" % len(pts), flush=True)
t0 = time.time()
for k in range(MAX_CHUNKS):
    s.DoDynamicsThenSync(C * DT)     # <- wedges here, typically k~30-40
    if k % 10 == 0:
        print("k=%d t=%.1f nc=%d" % (k, time.time() - t0, s.GetNumContacts()), flush=True)
print("DONE", flush=True)

What I ruled out (each falsified by a subsequent run that behaved oppositely)

• Adaptive kT update-frequency tuner — SetCDUpdateFreq(20) is initial-only; the tuner climbs it
  (20 → ~37 over ~30 chunks) and the CD margin / potential-pair count grow with it (~3× at the same
  stage vs a fixed frequency). DisableAdaptiveUpdateFreq() let one run complete 3000/3000 — but a
  later run stalled with the tuner disabled, so it is at most an aggravator.
• Adaptive bin sizing — GetBinSize()/GetBinNum() were stable up to the wedge; disabling it
  changed nothing.
• Host-side bulk getters racing the async protocol — a run calling GetOwnerPosition/
  GetOwnerVelocity every N chunks wedged at such a call, but a fully blind run (only
  DoDynamicsThenSync + GetNumContacts) also wedged. And the repro above wedges with neither.
• Custom force model — reproduces with the stock frictional-Hertzian model, as above.
• WSL2 / GPU virtualization — reproduces on bare-metal Linux with a native GPU (see Environment),
  so it is not a WSL GPU-scheduling artifact.

The through-line: no host-side call pattern reliably prevents it, which points at a race in the async
kT/dT engine, hazard concentrated at the contact ramp.

Workaround (in case it helps others landing here)

A supervisor that watchdogs stdout and, on >60 s of silence after Initialize(), kills and relaunches
the run — runs that clear the contact ramp almost always complete. Not a fix.

Questions

1. Is there a known race in the kT/dT handshake that a sharp contact-count increase can trigger?
2. Is there a supported way to run the solver fully synchronously (no async kinematic thread) as a
   diagnostic, to confirm the async engine is the locus?
3. Given this is bare-metal (not WSL2), is this the same root cause as Initialize hangs in pyDEME after setup on WSL Ubuntu 22.04 (CUDA 12.8, Blackwell) #70 or a separate issue?

Happy to run additional instrumentation (compute-sanitizer --tool racecheck, per-chunk timing, thread
dumps) if useful.

[Ward-Vandepitte]

Correction to the original report: I mislabeled the environment as WSL2 — this is actually bare-metal Ubuntu 24.04 with a native (non-virtualized) RTX 2000 Ada (nvidia-smi -q -> Virtualization Mode: None). The explorer.exe I saw on the GPU that led me to assume WSL was a Wine desktop stub from unrelated software on the same box, not a Windows host bleeding through a vGPU. I've updated the title and body accordingly.

Why it matters: the livelock is not a WSL GPU-scheduling artifact. Re-confirmed today on bare-metal — one roll of the minimal repro wedged at chunk k=19 (mid contact-ramp, nc~2.5k) with the GPU at 100% SM, one thread running (~33% CPU) and 23 in futex_wait_queue; a separate roll of the same scene completed 3000/3000 cleanly. So if this shares a root cause with #70 (WSL2), that cause is platform-independent; otherwise they are separate issues.