ADR-0011 — Non-positive density is a hard error, not a value to floor¶

Status: Accepted

Context¶

Per-cell source terms recover primitives from the conserved state by dividing by ρ — u = ρu/ρ, T = p/(ρR). The codebase had four different answers for what to do when ρ is not strictly positive, chosen ad hoc per closure:

convention	sites	behaviour at ρ ≤ 0
floor ρ	`filter.substrate_wall_into`, `doc_sulfur`, `hnco_hydrolysis`	continue on a fabricated 1e-30
skip the cell	the `@njit` kernels in `filter.py` (`if rho <= 0.0: continue`)	source term silently dropped
gate, then fall back	`porous_source` (`rho > 0.0 and p > 0.0` → constant μ)	partially effective (below)
nothing	`heat_source`, `catalyst_apply_cells`	divide through

A 2026-07-31 sweep for invariants guarded in one sibling but not the others flagged the inconsistency, but deferred the fix: floors and fail-loud serve opposite goals, so the convention had to be decided before it could be applied. The deferral rested on one premise — that leaving a site unguarded is safe because "under njit the NaN is caught one step later by the existing non-finite-dt guard."

That premise is false, and measuring it is what settled this ADR.

pipe_min_dt takes a max of max_wave_speed over cells. At ρ ≤ 0 the sound speed sqrt(γp/ρ) is NaN — and NaN loses every > comparison, so the bad cell simply fails to update smax. On a 200-cell pipe with one cell at ρ = −1, the guard returned a dt bit-identical to the fully healthy pipe:

all healthy       s[7]=371.565   dt=1.345659e-06
ONE cell rho=-1   s[7]=nan       dt=1.345659e-06   <-- NOT CAUGHT

Only a wholly degenerate pipe was ever caught, and then only because max_wave_speed happens to raise ZeroDivisionError at exactly ρ = 0 — the guard's own docstring says "an all-NaN … pipe". A localised divergence, which is what actually happens when a solution starts to go bad, passed straight through. Downstream, each closure then applied its own convention, so the observable behaviour depended on which closures the case happened to enable: a DPF floored and marched on, a catalyst divided through, porous_source returned −inf, filter silently dropped the term.

Two further measurements shaped the decision:

friction_source was never in this class. It multiplies by ρ (s_mom = -(f/2D_h)·ρ·|u|·u) and is well-behaved at ρ = 0. The earlier note listing it alongside heat_source was wrong.
porous_source's rho > 0.0 gate does not guard what it appears to. It prevents the division by zero but not the garbage: at ρ = 1e−12 the recovered T is enormous, so Sutherland's law returns a μ six orders too large (s_mom = −8.6e9), and at ρ = 1e−300 it returns −inf.

Decision¶

ρ ≤ 0 in any cell is an unphysical state and terminates the run with a diagnostic that names the pipe, the cell, and the value. It is never floored, skipped, or divided through as a matter of policy.

The check lives in the three dt kernels — kernels.pipe_min_dt, packed_pipes._packed_min_dt_{par,ser}, hllc_eos_njit.pipe_min_dt_eos_jit — which already visit every cell each step, so it costs one comparison in a loop that already does a sqrt. A bad cell routes to the established 1e300 degenerate sentinel and pipe_min_dt_checked re-diagnoses it into nonpositive_density_error.

Three details are deliberate:

A plain rho <= 0.0 comparison, not an isnan test. These kernels compile under fastmath, whose nnan flag makes NaN-based tests unreliable. An ordinary comparison against a finite value needs no NaN semantics: verified to flag −1.0, 0.0 and −inf identically under fastmath=True and fastmath=False. This is also why the guard's Python half stays in Python — matching the existing note that "the guard lives at Python level on purpose: the @njit kernel itself stays bit-identical".
The check must be in the kernel, not in numpy around it. The obvious U[:, 0].min() wrapper was measured against the kernel it guards:

N=50    pipe_min_dt= 1.344us   rho.min()=3.438us  (256% of the kernel)
N=200   pipe_min_dt= 4.072us   rho.min()=3.718us  ( 91%)
N=1000  pipe_min_dt=19.002us   rho.min()=4.733us  ( 25%)

Pipes in real cases are small — the EPW diesel manifold runners are n_cells = 15 — so numpy call overhead dominates and the wrapper would have been a serious per-step regression. In the kernel it is free.

Tiny-but-positive ρ is NOT thresholded. ρ = 1e−300 still passes. Any cutoff is arbitrary and risks rejecting a legitimate near-vacuum state, so the line is drawn at the unambiguous unphysical set. The existing closure-local floors remain as defence-in-depth for that grey zone — they are now belt-and-braces, not policy, and are commented as such.

Consequences¶

A localised divergence now stops the run naming the cause (Non-positive density in pipe 'ex_run_3': 1 of 200 cells have rho <= 0 (worst rho = -3.1e-04, cell index 7)) instead of completing with a plausible-looking wrong answer. This is the same silent-wrong-answer class as the wrap-around valve cam and the junction port-area mismatch.
Healthy runs are bit-identical. Every cell has ρ > 0, so the branch never fires and smax is unchanged. Confirmed under AEOLUS_FASTMATH=0.
A second, pre-existing defect is fixed in passing. Three kernels documented that "routing it to the 1e300 sentinel makes the caller raise the non-finite-dt error" — but 1e300 is finite, so pipe_min_dt_checked's isfinite test passed it straight through. network_timeloop and timeloop each re-tested >= 1e299 themselves and were safe; the two implicit_network call sites did not, and both take a min() across pipes, where a sentinel from one dead pipe loses to the healthy ones and vanishes. The test now lives in pipe_min_dt_checked, covering every caller.
The four conventions are not mechanically unified — that would have meant either adding floors (hiding divergence) or removing them (perturbing 15 validated suites). Instead the reachable case is now impossible, which makes the remaining differences unobservable rather than merely tolerated.

Alternatives considered¶

Floor ρ everywhere. Rejected: it converts a diverged solution into a finite, plausible, wrong answer that a user cannot distinguish from a good one — precisely the failure mode this project keeps having to fix.
Skip the cell everywhere (the filter.py convention). Rejected for the same reason: a silently dropped source term is still a silent wrong answer.
Detect the NaN wave speed inside the kernel rather than testing ρ. This would also catch tiny-positive ρ, which produces a NaN speed via an overflowed p. Rejected because it depends on NaN semantics inside a fastmath kernel.

Amended 2026-09-29 — non-positive pressure joins the refusal. A cell with ρ > 0 but ρ·(ρE) − ½(ρu)² <= 0 (p ≤ 0, including the tiny-ρ states whose pressure overflows) used to get a NaN wave speed, which loses every > and dropped the cell from the dt minimum silently. pipe_min_dt, its var-γ and packed twins now return the same sentinel for it, and pipe_min_dt_checked / the network driver raise a named "non-positive pressure" error. The test is the fastmath-safe product above, not a NaN check, so this ADR's rule stands. The EOS kernel (pipe_min_dt_eos_jit) is unchanged: under a real EOS, e ≤ 0 does not imply p ≤ 0. * A numpy positivity pass in pipe_min_dt_checked. Rejected on the measurement in Decision §2.