test: add opt-in shifted FP32 Jet CV control

docs/projected-explosion-experiment.md

@@ -1,62 +0,0 @@
-# Projected explosion: isolated spatial-interpolation experiment
-
-Historical baseline: `9805383f9c50d867e8555471baf5dbda8333a652`. This candidate does not replace or relabel its results. The TG runner records the actual candidate revision, dirty state, source hash and harness hashes, plus the historical baseline revision and selected explosion scheme.
-
-Default remains `constant`. Explicit `explosion="projected_linear"` supports only a periodic, body-free, two-level, fully blocked single rectangle (window or singleton union representation). Dense, multilevel, multi-rectangle and body/Jet cases are rejected. Constructor guards have been statically reviewed; GPU construction is not yet exercised.
-
-For each population, form a centered gradient only where both axis neighbors are original explosion donors, owned and uncovered. Otherwise set that axis derivative to zero. Project the gradient perpendicular to the population lattice direction (rest population retains its gradient), then evaluate at child offsets ±1/4 coarse cell. Read owner cores, never stale block halos. Repeating the donor for a missing axis avoids invalid reads. All neighbor reads remain within the original coarse transfer footprint.
-
-Timing, two ghost-only streaming hops, four-child coalescence, relaxation, shifted-storage convention, bulk collision and legacy constant kernel are unchanged. There is no NEQ rescaling, limiter, time interpolation or new framework. Projection commutes with subtracting the constant lattice weights.
-
-Method attribution: Schornbaum and Rüde, [arXiv:1508.07982](https://arxiv.org/pdf/1508.07982), Algorithm 2; Chen et al., Physica A 362 (2006) 158–167, [DOI 10.1016/j.physa.2005.09.036](https://doi.org/10.1016/j.physa.2005.09.036). Independently implemented mathematical stencil, with no GPL implementation copied. The Chen full text was not directly audited in this experiment.
-
-## CPU checks
-
-Run `PYTHONPATH=src python -m unittest discover -s tests/unit -p test_projected_explosion.py -v`.
-
-- Constant populations, affine available-axis gradients, characteristic constancy, quartet means and shift invariance.
-- Original donor footprint only; poison covered cells and check finite reconstruction.
-- Both window and singleton-union layouts, centered/off-axis/minimal near-edge rectangles.
-- Eight complete intervals, identity and SRT collision: actual fine pulls and coarse restriction tables, all eight grazing corner paths included; composite mass and both momentum components conserved within 2e-12 FP64 absolute tolerance.
-- Actual CUDA kernel body compiled with a CPU C++ shim: blocked address mapping, FP32 shifted/unshifted inputs, poisoned halo/covered storage, independent NumPy parity. This checks arithmetic and indexing, not CUDA compilation or execution.
-- Unprojected reconstruction is distinguished by failure of characteristic constancy. It need not fail global mass conservation in this geometry.
-
-## Next diagnostic, only on authorized Horsea GPU
-
-Use the same candidate checkout for paired Re100 N64 depth2 runs, differing only in `--explosion constant` versus `--explosion projected_linear`:
-
-`python tests/validation/taylor_green/run_case.py --repo "$PWD" --re 100 --n 64 --depth 2 --explosion projected_linear --output <new-result-path>`
-
-Select the physical CUDA device explicitly. Never overwrite historical reports. Compare frozen TG error gates and full trajectories, including mass, energy, pressure and velocity. Amplitude/nonmodal shape decomposition can inform mechanism but cannot substitute for the gates. First establish CUDA compilation and CPU/device parity, then paired end-to-end budgets. The CPU diagnostic below establishes improvement only in its FP64 reference trajectories; it does not establish production suitability or a candidate GPU pass.
-
-## Local verification limits
-
-Focused candidate/refinement/block-layout/build-contract checks pass (62 unittest tests). Full unittest discovery is not green: the attempted 200-test run reports 18 missing-dependency/import errors and three historical source-pin refusals (jet-circle and LES guards), plus one skipped test. Those pin refusals correctly reject the changed `blocks.cu`; their immutable fixtures were not weakened or rewritten. Python syntax checks and `git diff --check` pass. No candidate CUDA build, device execution, GPU performance or GPU TG improvement has been established.
-
-
-## Donor-support audit
-
-The current support is exactly the original two-coarse-cell-wide explosion ring. Centered normal derivatives therefore vanish along straight faces: the inner donor lacks an uncovered inward neighbor, and the outer donor lacks a third-ring outward neighbor. This is an explicit limitation, not evidence of a missing donor layer. Independent inspection of the authors' [LinearExplosion implementation](https://github.com/lssfau/walberla/blob/master/src/lbm/refinement/LinearExplosion.h) found the same two-cell directional slabs and both-neighbors-valid condition; their union equals all 368 distinct donors in the N64 TG geometry. No upstream code was copied. This does not claim identical block-local tangential fallback at every fine/fine block seam.
-
-A separate read-only dependency trace compared the ring to all uncovered periodic parent neighbors for four rectangles. Wider support changed some unused outer ghost populations but none consumed by final fine ACTIVE outputs or coarse restriction through the two hops, including full collision dependencies. Widening the support is not part of this candidate.
-
-## Reproducible CPU Taylor–Green A/B
-
-Run from this checkout, with NumPy and Python installed:
-
-`python tests/validation/taylor_green/run_cpu_explosion.py --repo "$PWD" --n 32 --output-dir /tmp/explosion-ab`
-
-Repeat with `--n 64` for the second retained case. These are the only supported grids. The script reuses unchanged CPUREF `advance_coarse_once`, existing TG initialization/metrics, actual singleton-rectangle transfer tables, and the tested projected reconstruction. The rectangular adapter verifies its coalescence children against the production table and includes all eight grazing-corner links. Only the explosion operation differs within each pair.
-
-The reference stores physical FP64 arrival populations; production stores shifted FP32 postcollision populations. Collision preserves the macroscopic moments compared here, and initialization is equilibrium. This supports a macro-trajectory comparison, not population-by-population/device equivalence. Outputs record actual checkout/harness provenance and hashes of both imported diagnostic helpers. The compact retained `tests/validation/taylor_green/cpu_explosion_summary.json` records original N32/N64 results and historical GPU report hashes, without importing those reports into this branch.
-
-At Re100, depth 2, t*=0.48:
-
-| Root N | Velocity L2, constant → projected | Pressure L2 | Orthogonal velocity residual |
-|---|---|---|---|
-| 32 | 3.24942% → 1.52264% | 8.98696% → 4.85601% | 2.28814% → 0.98615% |
-| 64 | 1.13113% → 0.49304% | 3.55050% → 1.65801% | 0.93699% → 0.40289% |
-
-Across sampled times, constant CPU versus historical constant GPU differs by at most 6.93e-8 in velocity L2 and 3.58e-6 in pressure L2 for N32; 1.74e-8 and 8.64e-7 for N64. Relative mass drift remains below 1.42e-13 and absolute momentum change below 2.11e-13 in these CPU pairs. Independent N32 rerunning reproduced trajectory/error metrics exactly; recomputed quadrature-factor metadata agrees to roundoff.
-
-Amplitude is the volume-weighted projection onto the exact velocity field; the orthogonal residual isolates the remaining shape error. The energy decomposition includes the discrete composite-grid exact-energy quadrature factor recorded per row, not simply amplitude squared plus residual squared. These are bounded diagnostic improvements, not proof of complete method fidelity, uniform-grid accuracy equivalence or a validated CUDA implementation. Frozen production validation gates are unchanged.

docs/validation/jet_circle_cv.md

@@ -113,3 +113,47 @@ These criteria are fixed before the continuation runs. There is no new physical
 Actual earlier scale-1 pilot cost was 52.3884s launch-to-exit, 46.1414s after READY for800 steps with full-field diagnostics every base step. A naive same-diagnostic3200-step scale-1 extrapolation is about210s; the new876 field samples reduce transfer count but add ray diagnostics and compression. Use conservative planning envelopes210–500s for scale1 and900–3200s for scale2, not timing guarantees. The scale2 conservative envelope accounts for4x cells and4x root steps; sparse sampling mitigates but does not eliminate that cost. Post-READY safety caps are1200s and5400s respectively, replacing the unsuitable fixed1800s assumption. Root owns compilation/startup timeout and serial launch decisions.
 
 Device-envelope estimates remain128/512MiB for s1/s2, at least1GiB free (or twice estimate). Output estimate is under64MiB per run including five primitive snapshots and876 sample records; no all-time populations are saved. Plan/metadata and sampled files are rewritten at each100-reference-time window to retain bounded partial evidence if interrupted. READY still follows compiled warmup and state/time reset, before measured stepping begins.
+
+## Dense-root temporal diagnostic: preparation only
+
+The audited continuation did not resolve time quadrature: all48 audited stride2 RMS changes exceeded their stride1 CV residuals, by factors1.55–11.46. That historical failure and its gates remain unchanged. The following is a measurement-only experiment to test this explanation, not a solver fix, a physics pass, or permission to extend GPU work. The separate periodic global stock-minus-raw mass residual also persists; denser contour quadrature cannot by itself resolve it.
+
+Future launch, only after the computer is online and root approves a serialized bounded run:
+
+```sh
+CUDA_VISIBLE_DEVICES="$VERIFIED_IDLE_GPU" python tests/validation/run_jet_circle_cv.py --repo "$PWD" --output "$NEW_DENSE_S1" --dense-root --steps 800 --scale 1
+CUDA_VISIBLE_DEVICES="$VERIFIED_IDLE_GPU" python tests/validation/run_jet_circle_cv.py --repo "$PWD" --output "$NEW_DENSE_S2" --dense-root --steps 800 --scale 2
+python tests/validation/analyze_jet_cv_continuation.py "$NEW_DENSE_S1" "$NEW_DENSE_S2" --output "$DENSE_COMPARISON_JSON"
+```
+
+This mode fixes Re100,beta.1 and reference t=0..800 (t*≈.9965), with800/3200 root intervals and801/3201 endpoint records. It is mutually exclusive with --continuation and rejects longer runs or other scales/Re. Physical geometry, source manifest, solver, velocity commands, kernel and raw observation path are unchanged. There are no new finest-hop hooks. Existing GPU population downloads occur every root interval for CPU primitive diagnostics; those distributions are discarded rather than retained.
+
+Dense scalar records contain CV endpoint inventories, each flux and split momentum terms, density/speed, physical composite mass, exact cumulative raw partitions and cumulative absolute aperture net mass. Rays remain on the exact old reference cadence:0..100 at dt1, then104..800 at dt4,276 samples for each scale. Schema2 maps their smaller arrays through ray_sample_indices into the dense records. Full primitive snapshots occur only at reference0,100,800. Saved data is capped at32MiB per case; dense arrays are scalar diagnostics, not all-time full fields. Schema/unit metadata, file hashes, complete-schedule checks and matched root time intervals remain mandatory.
+
+The analyzer derives root strides1,2,4 and the exact old cadence from the SAME dense primitive records, always keeping window endpoints and the SAME raw ledger. At scale2 these root strides correspond to reference dt.25,.5,1; the old cadence corresponds to root gaps4/16. Its fourth quadrature result and change from dense are separately identified. Existing residual normalizations and historical files are preserved. No automatic refinement-progress flags are issued for this temporal-only mode. Root-sampled agreement cannot certify finer-than-root time convergence or exclude aliasing; an unresolved outcome stays unresolved, with no automatic longer run.
+
+Three distinct throughput quantities must remain separate:
+
+- Prescribed unsigned nominal mass at rho_ref=1:2*Q_command*window_duration.
+- Observed raw absolute NET aperture mass: sum over root intervals of abs(dm_blowing)+abs(dm_suction), after spatial reduction within each aperture. This retains reversals between root intervals that abs(window-net) loses. Within-aperture or within-root/substep reversals remain unresolved by this ledger. The legacy raw_gross_mass_over_nominal field retains its old abs(window-net)/ONE-aperture-nominal definition; it is not this new quantity.
+- Sparse reconstructed absolute Q/M: true-arc quadrature of abs(u.n) and abs(rho*u.n), using each existing linear/quadratic ray estimate, followed by old-cadence time integration. This is estimator-dependent measured-field evidence, not exact wall truth, not dense-root quadrature, and not a prescribed nominal flow.
+
+Global stock-minus-raw mass balance is reported per root interval and summed per window, in reference mass units. Raw total-minus-wall/aperture partition discrepancies are also separate. Neither is bundled into the CV temporal error, and neither is called a proven raw-extractor defect. Source/geometry/state validity, both signed and gross budgets, old failures, and physics_validation_pass=false remain intact.
+
+Planning envelope from the earlier52-second800-step pilot: s1 60–240s, s2 600–2400s, with conservative post-READY caps1200/5400s. Scale2 has4x cells and4x root observations over the same reference horizon; sparse saving reduces output, not all CPU/PCIe observation cost. These are preparation estimates, not measured timings. Root owns startup/compilation timeout and GPU serialization. No GPU/remote action was taken to prepare this option.
+
+Budget enforcement preflights the projected retained file sizes before replacing each snapshot/window report or the sample+metadata pair. If32MiB would be exceeded, staged replacements are removed, the previous verified sample/hash-metadata pair remains intact, and execution stops. Temporary compression files need additional short-lived disk space and are distinct from the32MiB retained-evidence cap. This protects budget failures; it is not a filesystem crash-recovery transaction. Both raw partition identities, total−wall−aperture and aperture−blowing−suction, are saved separately from stock/raw and CV residuals.
+
+## Optional shifted-FP32 representation control
+
+The dense-root campaign accepts `--shifted-fp32` as an explicit numerical-representation control. Without it, the historical unshifted FP32 default remains unchanged. Use a fresh output directory and append the flag to the reviewed dense command, for example:
+
+```bash
+CUDA_VISIBLE_DEVICES="$VERIFIED_IDLE_GPU" python tests/validation/run_jet_circle_cv.py --repo "$PWD" --output "$NEW_SHIFTED_DENSE_S1" --dense-root --steps 800 --scale 1 --shifted-fp32
+```
+
+Root owns the future GPU launch and scale selection. This option requires the bounded dense campaign; it does not authorize another run. Geometry, Jet profile/amplitude, initial physical seed, viscosity, time horizon, root sampling, sparse rays/fields, normalization, thresholds, budgets and prior failed gates are unchanged. The `case.ddf_shifting` flag and `case.storage_variant` receipt appear in the plan, READY, window report and sample metadata. Preserve the unshifted evidence in its original directory. Compare like-scale records at identical root endpoints, including CV mass/momentum, global stock-minus-raw balance, raw partitions and measured aperture flux. A reduced mass residual would support representation sensitivity; it would not establish physical convergence or prove that every residual is roundoff.
+
+Both production blocked steppers already support this representation: initialization takes physical populations and encodes `g=f-w` to FP32; public download returns physical `f` decoded in FP64. The runner deliberately retains its historical FP32-rounded physical equilibrium seed for BOTH variants rather than also changing initial precision. Encoding introduces a small additional rounding difference: maximum per-population change 1.035e-10 and density change 3.201e-10 at both scales for this seed. Thus the decoded starting states are not bit-identical. The residual initial density offset is therefore not removed. All existing measurement routines receive physical populations and apply no extra weight shift. Snapshots keep the existing primitive-field precision.
+
+The existing `BlockCurvedBoundary` optional-input adapter invokes the same shifted-aware BFL helper as the AMR boundary: wall density restores the unit background, raw mass uses reflected-minus-outgoing (opposite weights cancel), and reaction restores twice the direction weight. No host correction is applied to raw rows. The shifted branch selects `compute_feq_from_momentum` rather than unshifted `compute_feq`, so this control changes existing numerical algebra as well as storage representation. It cannot isolate storage precision alone or identify collision as the sole cause. Block collision preserves the density perturbation in its shifted branch; blocked halo copies and AMR transfers already support shifted storage. This option changes no solver/kernel/API or source manifest. Existing shifted-storage/Jet host tests and the new runner codec tests are CPU evidence only; numerical GPU behavior remains to be checked by the reviewed control.

src/CelerisLab/lbm/kernels/blocks.cu

@@ -60,39 +60,6 @@ extern "C" __global__ void BlockExplode(
         to[index_f(t % TOTAL_CELLS, q)] = from[q * coarse_cells + d % coarse_cells];
 }
 
-// Experimental Schornbaum--Ruede (arXiv:1508.07982, Algorithm 2)
-// and Chen et al. (doi:10.1016/j.physa.2005.09.036) population interpolation.
-// Independent implementation of the mathematical stencil; no upstream source copied.
-// Stencil cells are owner cores; a disabled axis repeats its center twice.
-extern "C" __global__ void BlockExplodeProjected(
-    const fpxx* coarse, fpxx* fine, const uint64_t* target, const uint64_t* donor,
-    unsigned int count, unsigned long coarse_cells,
-    const uint64_t* neighbors, const float* offsets)
-{
-    unsigned int row = blockIdx.x * blockDim.x + threadIdx.x;
-    if (row >= count) return;
-    const uint64_t t = target[row], d = donor[row];
-    fpxx* to = fine + (t / TOTAL_CELLS) * NQ * TOTAL_CELLS;
-    for (unsigned int q = 0; q < NQ; ++q) {
-        float v[4];
-        for (unsigned int a = 0; a < 4; ++a) {
-            const uint64_t n = neighbors[4ul * row + a];
-            v[a] = coarse[(n / coarse_cells) * NQ * coarse_cells
-                          + q * coarse_cells + n % coarse_cells];
-        }
-        float gx = .5f * (v[1] - v[0]), gy = .5f * (v[3] - v[2]);
-        const float ex = d_cx[q], ey = d_cy[q];
-        if (q != 0) {
-            const float along = (ex * gx + ey * gy) / (ex * ex + ey * ey);
-            gx -= ex * along;
-            gy -= ey * along;
-        }
-        to[index_f(t % TOTAL_CELLS, q)] =
-            coarse[(d / coarse_cells) * NQ * coarse_cells + q * coarse_cells + d % coarse_cells]
-            + offsets[2ul * row] * gx + offsets[2ul * row + 1] * gy;
-    }
-}
-
 // Targets are receiver-local SOURCE slots for the immediately following pull.
 // Keep child summation order identical to PatchCoarseStep; no collision here.
 extern "C" __global__ void BlockCoalesce(

src/CelerisLab/lbm/refinement_stepper.py

@@ -22,7 +22,7 @@ from .block_stepper import MAX_GRID_DIM
 from .hierarchy import build_hierarchy, is_rectangle
 from .refinement import ACTIVE, HALO, FixedPatch2D, FixedPatchSet2D, prepare_fine_boundary
 from .refinement_tables import (blocked_base_tables, blocked_fine_tables,  # noqa: F401 (re-export)
-                                blocked_kind, projected_explosion_tables, _prepare_blocked_boundary,
+                                blocked_kind, _prepare_blocked_boundary,
                                 BLOCKED_BOUNDARY_DEVICE_COLUMNS)
 from .descriptors import SOLID
 from .raw_exchange import summarize_exchange
@@ -41,7 +41,6 @@ class _Pair:
     n_ghosts: int
     n_rows: int
     parent_cells: int    # parent block window cells (explode donor stride)
-    interpolation: tuple = ()  # optional neighbor/offset device tables
 
 
 @dataclass
@@ -217,7 +216,7 @@ class FixedPatchStepper:
                  fine_block=None, base_block=None, jets=(), child=None, child_block=None,
                  child_objects=(), child_jets=(), child_rectangles_in_parent=None,
                  refinements=None, level_blocks=None, level_objects=None, level_jets=None,
-                 observe_raw_exchange=False, explosion="constant"):
+                 observe_raw_exchange=False):
         if HALO != 4:
             raise RuntimeError("refinement.cu requires a four-cell fine halo")
         if coarse_boundary not in ("periodic", "slot_reservoir", "channel"):
@@ -246,15 +245,6 @@ class FixedPatchStepper:
         level_blocks = [None if b is None else tuple(b) for b in level_blocks]
         _check_levels(refinements, level_blocks, level_objects, level_jets)
         depth = len(refinements) + 1
-        if explosion not in ("constant", "projected_linear"):
-            raise ValueError("explosion must be constant or projected_linear")
-        if explosion == "projected_linear" and (
-                depth != 2 or coarse_boundary != "periodic"
-                or any(b is None for b in level_blocks)
-                or any(level_objects) or any(level_jets)
-                or not (is_rectangle(refinements[0]) or len(refinements[0]) == 1)):
-            raise ValueError("projected explosion requires periodic two-level blocked single rectangle without bodies/Jets")
-        self._explosion = explosion
         self._ddf_shifting = bool(cfg.ddf_shifting)
         if self._ddf_shifting and any(b is None for b in level_blocks):
             raise ValueError("shifted populations require all-blocked levels; dense/mixed is unsupported")
@@ -477,10 +467,6 @@ class FixedPatchStepper:
                        self._upload(children), len(ghosts), len(targets),
                        int(np.prod(pair.coarse_shape if parent.layout is None
                                    else parent.layout.local_shape)))
-        if self._explosion == "projected_linear":
-            neighbors, offsets = projected_explosion_tables(pair, parent.layout)
-            result.explode = child.module.get_function("BlockExplodeProjected")
-            result.interpolation = (self._upload(neighbors), self._upload(offsets))
         if child.layout is not None:
             child.kind = self._upload(blocked_kind(child.layout, self.hierarchy.kinds[index + 1]))
         if index == 0 and base is not None:
@@ -745,7 +731,7 @@ class FixedPatchStepper:
         """Explode from the given parent buffer, then two hops of level ``index``."""
         pair, level = self._pairs[index - 1], self._levels[index]
         pair.explode(parent_buffer, level.current, pair.ghosts, pair.donors, np.uint32(pair.n_ghosts),
-                     np.uint64(pair.parent_cells), *pair.interpolation,
+                     np.uint64(pair.parent_cells),
                      grid=((pair.n_ghosts + self._block[0] - 1) // self._block[0], 1), **kw)
         self._hop(index, level.current, level.spare, **kw)
         self._hop(index, level.spare, level.current, **kw)

src/CelerisLab/lbm/refinement_tables.py

@@ -66,40 +66,6 @@ def blocked_base_tables(patch, layout, order, wall_aware=False, *, composed_kind
             blocked_kind(layout, dense_kind))
 
 
-def projected_explosion_tables(patch, layout):
-    """Experimental periodic single-rectangle stencil, owner-core cell IDs.
-
-    Schornbaum & Rüde (2016), arXiv:1508.07982, Algorithm 2 explosion
-    interpolation, following Chen et al., Physica A 362 (2006), 158-167,
-    doi:10.1016/j.physa.2005.09.036: centered gradients, perpendicular projection,
-    and child offsets +/-1/4. Independently implemented from the mathematics.
-    Only original explosion donors may supply gradients: this preserves the
-    existing coarse transfer footprint. Missing axis pairs use center twice.
-    Caller must restrict use to two periodic, body-free blocked levels.
-    """
-    if hasattr(patch, "patches") and len(patch.patches) != 1:
-        raise ValueError("projected explosion requires one rectangle")
-    ny, nx = patch.coarse_shape
-    donor = patch.ghost_coarse.astype(np.intp)
-    legal = np.zeros(ny * nx, bool)
-    legal[donor] = True
-    legal &= layout.owned.ravel() & ~patch.covered.ravel()
-    if not legal[donor].all():
-        raise ValueError("projected explosion donors must be owned uncovered cells")
-    y, x = np.divmod(donor, nx)
-    neighbors = np.repeat(donor[:, None], 4, axis=1)
-    for axis, (dx, dy) in enumerate(((1, 0), (0, 1))):
-        minus = ((y - dy) % ny) * nx + (x - dx) % nx
-        plus = ((y + dy) % ny) * nx + (x + dx) % nx
-        valid = legal[minus] & legal[plus]
-        neighbors[valid, 2 * axis] = minus[valid]
-        neighbors[valid, 2 * axis + 1] = plus[valid]
-    fy, fx = np.divmod(patch.ghost_fine, patch.fine_shape[1])
-    offsets = np.column_stack((fx % 2, fy % 2)).astype(np.float32) * .5 - .25
-    return (np.ascontiguousarray(layout.owner.ravel()[neighbors], dtype=np.uint64),
-            np.ascontiguousarray(offsets))
-
-
 # Required device columns plus optional actuation; observation arms are uploaded
 # separately only when enabled. body_id/aperture remain host reduction metadata.
 BLOCKED_BOUNDARY_DEVICE_COLUMNS = ("block", "cell", "dir", "q", "fallback", "surface_velocity")

tests/unit/test_jet_circle_cv.py

@@ -28,6 +28,41 @@ def eq(rho, u, v):
 
 
 class Host(unittest.TestCase):
+    def test_shifted_storage_receipt_preserves_physical_case(self):
+        for scale in (1, 2):
+            base = case(scale)
+            shifted = case(scale, ddf_shifting=True)
+            self.assertFalse(base.pop("ddf_shifting"))
+            self.assertTrue(shifted.pop("ddf_shifting"))
+            self.assertEqual(base.pop("storage_variant")["name"], "unshifted_fp32")
+            self.assertEqual(shifted.pop("storage_variant")["name"], "shifted_fp32")
+            self.assertEqual(base, shifted)
+
+    def test_shifted_physical_download_roundtrip(self):
+        from CelerisLab.lbm.population_storage import _encode_physical, _decode_stored
+        from run_jet_circle_cv import moments
+        from jet_cv_measurements import primitives
+
+        for speed in (0.03, 0.015):
+            seed = state((9, 48, 64), speed)
+            self.assertEqual(seed.dtype, np.float32)
+            plain = _decode_stored(_encode_physical(seed, False), False)
+            decoded = _decode_stored(_encode_physical(seed, True), True)
+            self.assertEqual(decoded.dtype, np.float64)
+            np.testing.assert_array_equal(plain, seed)
+            # Same FP32 physical seed; only encoding quantization may differ.
+            np.testing.assert_allclose(decoded, seed, rtol=0, atol=3e-10)
+            rho, momentum, velocity = moments(decoded)
+            rho2, velocity2 = primitives(decoded, CX, CY)
+            np.testing.assert_array_equal(rho, rho2)
+            np.testing.assert_array_equal(velocity, velocity2)
+            np.testing.assert_allclose(momentum, rho * velocity, atol=1e-17)
+            fluid = np.ones(seed.shape[1:], bool)
+            inv, flux = cv(decoded, fluid, (0, 0), (5, 4, 25, 20), 0.0072)
+            expected, _ = cv(seed, fluid, (0, 0), (5, 4, 25, 20), 0.0072)
+            np.testing.assert_allclose(inv, expected, rtol=0, atol=1e-6)
+            np.testing.assert_allclose(flux, 0, atol=1e-14)
+
     def test_descriptor_and_uniform_closed_flux(self):
         from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
 

tests/unit/test_jet_cv_measurements.py

@@ -31,6 +31,111 @@ from analyze_jet_cv_continuation import integrate, analyze
 
 
 class Measurements(unittest.TestCase):
+    def test_budget_failure_preserves_verified_pair(self):
+        from jet_cv_measurements import install_evidence
+
+        with tempfile.TemporaryDirectory() as tmp:
+            path = Path(tmp)
+            data = path / "samples.npz"
+            meta = path / "samples.json"
+            data.write_bytes(b"old-data")
+            meta.write_text("old-meta")
+            staged = path / ".samples.pending.npz"
+            staged_meta = path / ".meta.pending"
+            staged.write_bytes(b"x" * 100)
+            staged_meta.write_text("new-meta")
+            with self.assertRaisesRegex(RuntimeError, "last verified files preserved"):
+                install_evidence({data: staged, meta: staged_meta}, limit=20)
+            self.assertEqual(data.read_bytes(), b"old-data")
+            self.assertEqual(meta.read_text(), "old-meta")
+            self.assertEqual(
+                sorted(p.name for p in path.iterdir()), ["samples.json", "samples.npz"]
+            )
+
+    def test_dense_root_and_exact_old_cadence(self):
+        from jet_cv_measurements import old_sample_steps
+
+        for scale, count in ((1, 801), (2, 3201)):
+            c = case(scale, steps=800)
+            c["dense_root"] = True
+            dense = sample_steps(c)
+            old = old_sample_steps(c)
+            self.assertEqual(len(dense), count)
+            self.assertEqual(len(old), 276)
+            np.testing.assert_array_equal(dense, np.arange(800 * scale**2 + 1))
+            self.assertTrue(np.isin(old, dense).all())
+            np.testing.assert_array_equal(old[:101] / scale**2, np.arange(101))
+            np.testing.assert_array_equal(old[101:] / scale**2, np.arange(104, 801, 4))
+
+    def test_root_interval_absolute_net_mass_reversals(self):
+        from jet_cv_measurements import absolute_aperture_interval_mass
+
+        masks = (np.array([True, False]), np.array([False, True]))
+        rows = np.zeros((2, 4))
+        rows[0, 3] = 1
+        first = absolute_aperture_interval_mass(rows, *masks)
+        rows[0, 3] = -1
+        self.assertEqual(first + absolute_aperture_interval_mass(rows, *masks), 2)
+        self.assertEqual(abs(1 - 1), 0)
+
+    def test_dense_offline_oscillation_preserves_ledger(self):
+        import json
+        from jet_cv_measurements import old_sample_steps
+
+        for s in (1, 2):
+            c = case(s, steps=800)
+            c["dense_root"] = True
+            steps = sample_steps(c)
+            t = steps / s**2
+            n = len(t)
+            ridx = old_sample_steps(c)
+            # Nonzero steady reference fluxes cancel exact stock changes.
+            inv = np.zeros((n, 2, 2, 3))
+            inv[..., 0] = (100 - 0.001 * t[:, None, None]) * s**2
+            inv[..., 1] = -t[:, None, None] * 0.002 * s
+            flux = np.zeros_like(inv)
+            flux[..., 0] = 0.001 + 1e-4 * np.cos(np.pi / 2 * t[:, None, None])
+            flux[..., 1] = 0.002 / s
+            raw = np.zeros((n, 2, 5, 4))
+            with tempfile.TemporaryDirectory() as tmp:
+                path = Path(tmp)
+                (path / f"samples_s{s}.json").write_text(
+                    json.dumps(dict(case=c, sample_file="samples.npz"))
+                )
+                np.savez(
+                    path / "samples.npz",
+                    t_reference=t,
+                    cv_inventory=inv,
+                    cv_flux=flux,
+                    raw_cumulative=raw,
+                    ray_sample_indices=ridx,
+                    arc_estimates=np.zeros((len(ridx), 2, 2, 2, 2)),
+                    arc_absolute_estimates=np.zeros((len(ridx), 2, 2, 2, 2)),
+                    raw_absolute_aperture_mass_cumulative=np.broadcast_to(
+                        2 * t[:, None] * s**2, (n, 2)
+                    ),
+                    global_mass_reference=np.ones((n, 2)) * 100,
+                    level_stats=np.broadcast_to([1, 1, 0, 0.03, 100], (n, 3, 5)),
+                    velocity_error=np.zeros(n),
+                )
+                result = analyze(path)
+                w = result["windows"][1]
+                self.assertTrue(result["complete_schedule"])
+                np.testing.assert_allclose(
+                    np.array(w["normalized_cv_residual"])[0], 0, atol=1e-11
+                )
+                self.assertGreater(
+                    abs(np.array(w["old_cadence_change_from_dense"])[..., 0]).max(),
+                    0.005,
+                )
+                np.testing.assert_allclose(
+                    w["global_stock_raw_balance_sum_reference"], 0
+                )
+                np.testing.assert_allclose(
+                    w["raw_absolute_aperture_net_interval_mass_reference"], 200
+                )
+                np.testing.assert_allclose(w["raw_reaction_coefficient"], 0)
+
     def test_common_physical_schedule(self):
         a = sample_steps(case(1, steps=3200))
         b = sample_steps(case(2, steps=3200))
@@ -134,6 +239,57 @@ class Measurements(unittest.TestCase):
                 snap["uniform_velocity_reference"][0][fluid], 0.01, rtol=1e-7
             )
 
+    def test_dense_recorder_retains_sparse_rays(self):
+        h = 0.5
+        y, x = np.mgrid[:64, :64]
+        x = (x + 0.5) * h
+        y = (y + 0.5) * h
+        fluid = (x - 16) ** 2 + (y - 16) ** 2 > 25
+        fine_leaf = fluid.copy()
+        fine_leaf[:, :4] = False
+        root_leaf = np.zeros((32, 32), bool)
+        root_leaf[:, :2] = True
+
+        def eq(ny, nx):
+            cu = CX[:, None, None] * 0.01
+            return np.broadcast_to(
+                W[:, None, None] * (1 + 3 * cu + 4.5 * cu * cu - 1.5 * 0.01**2),
+                (9, ny, nx),
+            ).copy()
+
+        ug = eq(64, 64)
+        levels = [eq(32, 32), ug.copy()]
+        c = dict(
+            scale=1,
+            dense_root=True,
+            steps=800,
+            center=[16, 16],
+            diameter=10,
+            u=0.01,
+            nu=0.0072,
+            contours=[[9, 9, 23, 23], [7, 7, 25, 25]],
+        )
+        data = {
+            name: [cv(ug, mask, (0, 0), r, c["nu"]) for r in c["contours"]]
+            for name, mask in [("uniform_adapter", fluid), ("amr_api", fine_leaf)]
+        }
+        with tempfile.TemporaryDirectory() as tmp:
+            rec = MeasurementRecorder(
+                c, tmp, fluid, [root_leaf, fine_leaf], (0, 0), CX, CY
+            )
+            for step in (0, 101):
+                rec.capture(step, ug, levels, data, {k: np.zeros(2) for k in data}, 0.0)
+                rec.records[-1]["raw_cumulative"] = np.zeros((2, 5, 4))
+                rec.records[-1]["raw_absolute_aperture_mass_cumulative"] = np.zeros(2)
+            rec.flush()
+            a = np.load(Path(tmp) / "samples_s1.npz")
+            np.testing.assert_array_equal(a["step"], [0, 101])
+            np.testing.assert_array_equal(a["ray_sample_indices"], [0])
+            self.assertEqual(a["arc_estimates"].shape, (1, 2, 2, 2, 2))
+            self.assertTrue(
+                np.all(a["arc_absolute_estimates"] >= abs(a["arc_estimates"]) - 1e-14)
+            )
+
     def test_offline_mass_force_normalization(self):
         import json
 

tests/unit/test_projected_explosion.py

@@ -1,204 +0,0 @@
-"""CPU stencil and actual closed-interface packet budgets; no CUDA claims."""
-import numpy as np
-import unittest
-from pathlib import Path
-import shutil
-import subprocess
-import tempfile
-
-from CelerisLab.config import LBMConfig
-from CelerisLab.lbm.blocks import BlockLayout2D
-from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
-from CelerisLab.lbm.refinement import FixedPatch2D, FixedPatchSet2D, ACTIVE, OUTER_GHOST
-from CelerisLab.lbm.refinement_tables import projected_explosion_tables
-
-E = np.array([D2Q9_EX, D2Q9_EY]).T
-W = np.array([4/9] + [1/9]*4 + [1/36]*4)
-
-
-def setup(union=False, bounds=(4, 4, 12, 12)):
-    cfg = LBMConfig(nx=16, ny=16, viscosity=.1)
-    p = (FixedPatchSet2D.from_config(cfg, [bounds]) if union else
-         FixedPatch2D.from_config(cfg, bounds))
-    layout = BlockLayout2D.build(16, 16, 4, 4, periodic=(True, True))
-    neighbors, offsets = projected_explosion_tables(p, layout)
-    inverse = np.full(layout.nblocks * np.prod(layout.local_shape), -1, dtype=int)
-    inverse[layout.owner.ravel()] = np.arange(256)
-    return p, inverse[neighbors], offsets
-
-
-def reconstruct(c, p, neighbors, offsets, project=True):
-    # Independent NumPy evaluation of the documented stencil; neighbors have
-    # been decoded from actual blocked owner IDs, not recomputed geometry.
-    values = c.reshape(9, -1)[:, neighbors]
-    g = .5 * (values[:, :, 1::2] - values[:, :, ::2])
-    if project:
-        g[1:] -= (np.einsum('qra,qa->qr', g[1:], E[1:]) /
-                  np.sum(E[1:]**2, axis=1)[:, None])[:, :, None] * E[1:, None, :]
-    out = c.reshape(9, -1)[:, p.ghost_coarse] + np.einsum('qra,ra->qr', g, offsets)
-    return out, g
-
-
-def check_legal_footprint_mean_and_characteristic(union, bounds):
-    p, neighbors, offsets = setup(union, bounds)
-    assert set(neighbors.ravel()) <= set(p.ghost_coarse)
-    assert not p.covered.ravel()[neighbors].any()
-    assert (neighbors >= 0).all()
-    c = np.random.default_rng(72).normal(size=(9,16,16))
-    c[:, p.covered] = np.nan  # poison every forbidden donor
-    out, g = reconstruct(c, p, neighbors, offsets)
-    assert np.isfinite(out).all()
-    np.testing.assert_allclose(np.einsum('qra,qa->qr', g, E), 0, atol=5e-16)
-    for d in np.unique(p.ghost_coarse):
-        rows = p.ghost_coarse == d
-        assert rows.sum() == 4
-        np.testing.assert_allclose(out[:, rows].mean(axis=1), c.reshape(9,-1)[:,d], atol=5e-16)
-        for q in range(1,9):
-            delta = offsets[rows, None, :] - offsets[None, rows, :]
-            # Parallel displacements have zero 2D cross product.
-            parallel = np.isclose(delta[...,0]*E[q,1]-delta[...,1]*E[q,0], 0)
-            diff = out[q, rows][:,None]-out[q, rows][None,:]
-            np.testing.assert_allclose(diff[parallel], 0, atol=5e-16)
-    const = np.broadcast_to(W[:,None,None], (9,16,16)).copy()
-    constant, _ = reconstruct(const,p,neighbors,offsets)
-    np.testing.assert_array_equal(constant, np.broadcast_to(W[:,None], constant.shape))
-
-
-def collide(f, tau):
-    rho = f.sum(axis=0)
-    u = np.einsum('qa,q...->a...', E, f) / rho
-    cu = np.einsum('qa,a...->q...', E, u)
-    eq = W[:,None] * rho * (1 + 3*cu + 4.5*cu**2 - 1.5*np.sum(u*u,axis=0))
-    return f + (eq-f)/tau
-
-
-def totals(c, f, p):
-    q = c[:,~p.covered].sum(axis=1)+.25*f[:,p.fine_kind==ACTIVE].sum(axis=1)
-    return np.r_[q.sum(), E.T@q]
-
-
-def check_full_two_hop_closed_corner_budget(union, collision, bounds):
-    p, neighbors, offsets = setup(union,bounds)
-    rng=np.random.default_rng(17)
-    c=W[:,None,None]*(1+rng.normal(0,.002,(9,16,16)))
-    f=W[:,None,None]*(1+rng.normal(0,.002,(9,*p.fine_shape)))
-    start=totals(c,f,p)
-    # Count exactly the eight coarse diagonal links whose source is also
-    # uncovered, but whose midpoint lies on the closed rectangular interface.
-    q, k=np.divmod(p.coalesce_coarse.astype(np.int64),256)
-    y,x=np.divmod(k,16)
-    grazing=~p.covered[y-E[q,1],x-E[q,0]]
-    assert grazing.sum()==8
-    c[:,p.covered]=np.nan
-    for _ in range(8):
-        f.reshape(9,-1)[:,p.ghost_fine]=reconstruct(c,p,neighbors,offsets)[0]
-        for _ in range(2):
-            f=np.array([np.roll(f[q],(ey,ex),(0,1)) for q,(ex,ey) in enumerate(E)])
-            f[:,p.fine_kind==OUTER_GHOST]=np.nan
-            if collision:
-                f[:,p.fine_kind==ACTIVE]=collide(f[:,p.fine_kind==ACTIVE],p.tau_fine)
-        out=np.array([np.roll(c[q],(ey,ex),(0,1)) for q,(ex,ey) in enumerate(E)])
-        packets=f.reshape(-1)[p.coalesce_fine].mean(axis=1)
-        assert np.isfinite(packets[grazing]).all()
-        out.reshape(-1)[p.coalesce_coarse]=packets
-        if collision:out[:,~p.covered]=collide(out[:,~p.covered],p.tau_coarse)
-        out[:,p.covered]=np.nan
-        c=out
-        np.testing.assert_allclose(totals(c,f,p),start,rtol=0,atol=2e-12)
-
-
-def check_projection_is_distinguished_by_characteristic_not_global_mass():
-    p,n,o=setup()
-    c=np.random.default_rng(8).normal(size=(9,16,16))
-    _,g=reconstruct(c,p,n,o,project=False)
-    assert np.max(abs(np.einsum('qra,qa->qr',g,E))) > .1
-
-
-def check_multiple_rectangles_rejected():
-    cfg=LBMConfig(nx=24,ny=16,viscosity=.1)
-    p=FixedPatchSet2D.from_config(cfg, [(4,4,8,12),(12,4,16,12)])
-    layout=BlockLayout2D.build(16,24,4,4,periodic=(True,True))
-    with unittest.TestCase().assertRaisesRegex(ValueError, 'one rectangle'):
-        projected_explosion_tables(p,layout)
-
-
-class ProjectedExplosionTests(unittest.TestCase):
-    def test_stencil(self):
-        for union in (False, True):
-            for bounds in ((4,4,12,12),(2,3,9,10),(2,2,3,3)):
-                with self.subTest(union=union,bounds=bounds):
-                    check_legal_footprint_mean_and_characteristic(union,bounds)
-
-    def test_closed_interface(self):
-        for union in (False,True):
-            for collision in (False,True):
-                for bounds in ((4,4,12,12),(2,3,14,12),(2,2,3,3)):
-                    with self.subTest(union=union,collision=collision,bounds=bounds):
-                        check_full_two_hop_closed_corner_budget(union,collision,bounds)
-
-    def test_projection(self):
-        check_projection_is_distinguished_by_characteristic_not_global_mass()
-
-    def test_affine_and_shift_invariance(self):
-        p,n,o=setup()
-        y,x=np.indices((16,16))
-        c=np.array([.7+q*.01+x*.03+y*.05 for q in range(9)])
-        out,g=reconstruct(c,p,n,o)
-        available=n[:,::2]!=n[:,1::2]
-        gradient=np.array([.03,.05])*available
-        expected=np.broadcast_to(gradient,(9,*gradient.shape)).copy()
-        expected[1:]-=(np.einsum('qra,qa->qr',expected[1:],E[1:])/
-                       np.sum(E[1:]**2,axis=1)[:,None])[:,:,None]*E[1:,None,:]
-        np.testing.assert_allclose(g,expected,atol=2e-16)
-        shifted,_=reconstruct(c-W[:,None,None],p,n,o)
-        np.testing.assert_allclose(out-W[:,None],shifted,atol=5e-16)
-
-    def test_reject(self):
-        check_multiple_rectangles_rejected()
-
-
-class ProjectedKernelTests(unittest.TestCase):
-    @unittest.skipUnless(shutil.which("g++"), "g++ is required for actual-kernel CPU arithmetic")
-    def test_actual_kernel_fp32_shifted_and_unshifted(self):
-        p, neighbors, offsets = setup(True)
-        parent = BlockLayout2D.build(16,16,4,4,periodic=(True,True))
-        child = p.fine_layout((8,4))  # deliberately unequal parent/child strides
-        ids, offsets = projected_explosion_tables(p,parent)
-        cells = int(np.prod(parent.local_shape))
-        targets = child.owner.ravel()[p.ghost_fine]
-        donors = parent.owner.ravel()[p.ghost_coarse]
-        source = (Path(__file__).resolve().parents[2] /
-                  "src/CelerisLab/lbm/kernels/blocks.cu").read_text()
-        kernel = source.split('extern "C" __global__ void BlockExplodeProjected(',1)[1]
-        kernel = 'void BlockExplodeProjected(' + kernel.split('// Targets are receiver-local SOURCE slots',1)[0]
-        def array(name, values, kind):
-            v = np.asarray(values).ravel()
-            text = ','.join(str(int(x)) if kind=='uint64_t' else float(x).hex()+'f' for x in v)
-            return f"{kind} {name}[] = {{{text}}};\n"
-        for shifted in (False,True):
-            c=np.random.default_rng(11).normal(0,.001,(9,16,16))
-            if not shifted:c+=W[:,None,None]
-            c=c.astype(np.float32)
-            packed=np.full(parent.nblocks*9*cells,np.nan,np.float32)
-            for q in range(9):packed[parent.slot(parent.owner.ravel(),q)]=c[q].ravel()
-            # All halo/covered inputs are poisoned; kernel must use legal cores.
-            for q in range(9):packed[parent.slot(parent.owner.ravel()[p.covered.ravel()],q)]=np.nan
-            # Hex initializer needs a portable expression for poisoned NaNs.
-            data=array('coarse',packed,'float').replace('nanf','NAN')
-            code = '#include <cstdint>\n#include <cmath>\n#include <cstdio>\nusing fpxx=float;\n'
-            code += f'constexpr unsigned NQ=9, TOTAL_CELLS={int(np.prod(child.local_shape))};\n'
-            code += 'struct Id {unsigned x;}; Id blockIdx{0},blockDim{1},threadIdx{0};\n'
-            code += 'unsigned index_f(unsigned k,unsigned q){return q*TOTAL_CELLS+k;}\n'
-            code += 'int d_cx[]={0,1,-1,0,0,1,-1,1,-1},d_cy[]={0,0,0,1,-1,1,-1,-1,1};\n'+kernel
-            code += data+array('target',targets,'uint64_t')+array('donor',donors,'uint64_t')
-            code += array('neighbor',ids,'uint64_t')+array('offset',offsets,'float')
-            code += f'float fine[{child.nblocks*9*int(np.prod(child.local_shape))}]={{}};\n'
-            code += f'int main(){{for(unsigned r=0;r<{len(targets)};++r){{threadIdx.x=r;BlockExplodeProjected(coarse,fine,target,donor,{len(targets)},{cells},neighbor,offset);}}'
-            code += f'for(unsigned q=0;q<9;++q)for(unsigned r=0;r<{len(targets)};++r){{auto t=target[r];printf("%.9g\\n",fine[(t/TOTAL_CELLS)*9*TOTAL_CELLS+q*TOTAL_CELLS+t%TOTAL_CELLS]);}}}}'
-            with tempfile.TemporaryDirectory() as folder:
-                src=Path(folder)/'kernel.cpp'; binary=Path(folder)/'kernel'
-                src.write_text(code)
-                subprocess.run(['g++','-std=c++17','-O2',str(src),'-o',str(binary)],check=True,capture_output=True)
-                result=np.fromstring(subprocess.check_output([str(binary)],text=True),sep=' ').reshape(9,-1)
-            expected,_=reconstruct(c.astype(float),p,neighbors,offsets)
-            np.testing.assert_allclose(result,expected,atol=6e-8,rtol=2e-7)

tests/validation/analyze_jet_cv_continuation.py

@@ -32,7 +32,7 @@ def analyze(directory):
     t = a["t_reference"]
     if not len(t) or t[0] != 0 or np.any(np.diff(t) <= 0):
         raise ValueError("Sample reference times must increase strictly from zero")
-    from jet_cv_measurements import sample_steps
+    from jet_cv_measurements import sample_steps, old_sample_steps
 
     complete_schedule = np.array_equal(t, sample_steps(c) / s**2)
     records = []
@@ -45,7 +45,15 @@ def analyze(directory):
     )  # angular impulse reference = root/s^2 (not used in force tests)
     q0 = 2 * np.pi / 12 * (c["diameter"] / 2) * c["jet_speed"]
     f0 = 0.5 * (c["u"] * s) ** 2 * (c["diameter"] / s)
-    for lo in range(0, 3200, 100):
+    dense_root = c.get("dense_root", False)
+    old_times = old_sample_steps(c) / s**2
+    ray_indices = (
+        a["ray_sample_indices"] if "ray_sample_indices" in a else np.arange(len(t))
+    )
+    ray_times = t[ray_indices]
+    if dense_root and complete_schedule and not np.array_equal(ray_times, old_times):
+        raise ValueError("Sparse ray records do not match the declared old cadence")
+    for lo in range(0, c.get("steps", 3200 * s * s) // (s * s), 100):
         hi = lo + 100
         ids = np.flatnonzero((t >= lo) & (t <= hi))
         if len(ids) < 2 or t[ids[0]] != lo or t[ids[-1]] != hi:
@@ -54,8 +62,14 @@ def analyze(directory):
         delta = inv[last] - inv[first]
         ledger = raw[last] - raw[first]
         cv_variants = []
-        for stride in (1, 2, 4):
-            selected = np.unique(np.r_[ids[::stride], ids[-1]])
+        selections = [np.unique(np.r_[ids[::stride], ids[-1]]) for stride in (1, 2, 4)]
+        if dense_root:
+            selections.append(ids[np.isin(t[ids], old_times)])
+        for selected in selections:
+            if t[selected[0]] != lo or t[selected[-1]] != hi:
+                raise ValueError(
+                    "Every quadrature selection must retain both window endpoints"
+                )
             total = integrate(t[selected], flux[selected])
             m = delta[:, :, 0] + total[:, :, 0]
             j = -delta[:, :, 1:] - total[:, :, 1:]
@@ -68,16 +82,61 @@ def analyze(directory):
                     axis=-1,
                 )
             )
+        ray_ids = np.flatnonzero((ray_times >= lo) & (ray_times <= hi))
         arc = (
-            integrate(t[ids], a["arc_estimates"][ids]) / 100
+            integrate(ray_times[ray_ids], a["arc_estimates"][ray_ids]) / 100
         )  # Q/M rate conversion is unity
         # Each estimator, aperture, Q/M; compare against signed prescribed Q (rho_ref=1 for M only).
         arc_ratio = arc / np.array([q0, -q0])[None, None, :, None]
         global_delta = (
             a["global_mass_reference"][last] - a["global_mass_reference"][first]
         )
+        dense_details = {}
+        if dense_root:
+            absolute = (
+                a["raw_absolute_aperture_mass_cumulative"][last]
+                - a["raw_absolute_aperture_mass_cumulative"][first]
+            ) / s**2
+            stock_raw = np.diff(a["global_mass_reference"][ids], axis=0) - np.diff(
+                raw[ids, :, 2, 3], axis=0
+            )
+            partition = ledger[:, 2] - ledger[:, 0] - ledger[:, 1]
+            dense_details = dict(
+                quadrature_names=[
+                    "root_stride1",
+                    "root_stride2",
+                    "root_stride4",
+                    "exact_old_cadence",
+                ],
+                root_dt_reference=1 / s**2,
+                old_cadence_change_from_dense=(
+                    cv_variants[3] - cv_variants[0]
+                ).tolist(),
+                root_temporal_sensitivity_max=float(
+                    np.max(abs(np.array(cv_variants)[1:3] - cv_variants[0]))
+                ),
+                raw_absolute_aperture_net_interval_mass_reference=absolute.tolist(),
+                prescribed_unsigned_nominal_mass_reference=float(2 * q0 * 100),
+                legacy_raw_gross_field_definition="raw_gross_mass_over_nominal retains abs(window net blowing)+abs(window net suction), divided by ONE nominal aperture throughput; it is not root-interval absolute throughput.",
+                raw_absolute_to_prescribed_unsigned_nominal=(
+                    absolute / (2 * q0 * 100)
+                ).tolist(),
+                global_stock_raw_balance_max_abs_root_interval_reference=np.max(
+                    abs(stock_raw), axis=0
+                ).tolist(),
+                global_stock_raw_balance_sum_reference=stock_raw.sum(axis=0).tolist(),
+                raw_total_minus_wall_aperture_reference=partition.tolist(),
+                raw_aperture_minus_blowing_suction_reference=(
+                    ledger[:, 1] - ledger[:, 3] - ledger[:, 4]
+                ).tolist(),
+                measured_ray_absolute_volume_mass_reference=integrate(
+                    ray_times[ray_ids], a["arc_absolute_estimates"][ray_ids]
+                ).tolist(),
+                measured_ray_cadence="Sparse old cadence; estimator-dependent, not dense-root or exact wall throughput",
+            )
         records.append(
             dict(
+                **dense_details,
                 t_reference=[lo, hi],
                 startup=hi <= 800,
                 normalized_cv_residual=np.array(cv_variants).tolist(),
@@ -97,6 +156,12 @@ def analyze(directory):
         )
     return dict(
         scale=s,
+        dense_root=dense_root,
+        previous_temporal_failure_retained=(
+            "Prior continuation stride2 changes exceeded baseline residual in all48 audited comparisons; this diagnostic does not retroactively change those gates."
+            if dense_root
+            else None
+        ),
         complete_schedule=bool(complete_schedule),
         sample_checksum_verified=checksum_verified,
         windows=records,
@@ -122,8 +187,10 @@ def main():
         results=results,
         physics_validation_pass=False,
     )
-    if sorted(x["scale"] for x in results) == [1, 2] and all(
-        x["complete_schedule"] for x in results
+    if (
+        not any(x.get("dense_root", False) for x in results)
+        and sorted(x["scale"] for x in results) == [1, 2]
+        and all(x["complete_schedule"] for x in results)
     ):
         ordered = sorted(results, key=lambda x: x["scale"])
         summary = []

tests/validation/jet_cv_measurements.py

@@ -10,13 +10,63 @@ LINEAR = np.array([2.5, -1.5, 0.0])
 QUADRATIC = np.array([4.375, -5.25, 1.875])
 
 
+EVIDENCE_LIMIT = 32 * 1024**2
+
+
+def install_evidence(replacements, limit=None):
+    """Preflight retained bytes before replacing files; cap failure preserves old pair.
+
+    Staged compression bytes are temporary disk overhead, not retained evidence.
+    This guards budget failures, not a transactional crash-recovery filesystem.
+    """
+    temporary = set(replacements.values())
+    try:
+        directory = next(iter(replacements)).parent
+        retained = sum(
+            p.stat().st_size
+            for p in directory.iterdir()
+            if p.is_file() and p not in temporary
+        )
+        projected = retained - sum(p.stat().st_size for p in replacements if p.exists())
+        projected += sum(p.stat().st_size for p in temporary)
+        if limit is not None and projected > limit:
+            raise RuntimeError(
+                "Dense-root 32 MiB retained-evidence budget would be exceeded; last verified files preserved"
+            )
+        for target, staged in replacements.items():
+            staged.replace(target)
+    finally:
+        for p in temporary:
+            if p.exists():
+                p.unlink()
+
+
+def write_evidence_text(path, text, *, dense_root=False):
+    staged = path.with_name("." + path.name + ".pending")
+    staged.write_text(text)
+    install_evidence({path: staged}, EVIDENCE_LIMIT if dense_root else None)
+
+
 def sample_steps(c):
     """Same reference times at both scales; dense startup retained."""
+    if c.get("dense_root", False):
+        return np.arange(c["steps"] + 1, dtype=int)
+    return old_sample_steps(c)
+
+
+def old_sample_steps(c):
+    """Exact previous physical cadence, clipped to this campaign's endpoint."""
     s2 = c["scale"] ** 2
-    reference = sorted(set(range(0, 101)) | set(range(100, 3201, 4)))
+    end = c.get("steps", 3200 * s2) // s2
+    reference = sorted(set(range(0, min(100, end) + 1)) | set(range(100, end + 1, 4)))
     return np.array(reference, int) * s2
 
 
+def absolute_aperture_interval_mass(rows, blowing, suction):
+    """Sum absolute aperture NET masses for one root interval, not per-link flux."""
+    return float(abs(rows[blowing, 3].sum()) + abs(rows[suction, 3].sum()))
+
+
 def factors(scale):
     """Lref=one s=1 root length, Tref=one s=1 root time, rho_ref=1."""
     return dict(
@@ -105,7 +155,7 @@ def ray_geometry(fluid, origin, center, radius, h=0.5):
     )
 
 
-def reconstruct_aperture(rho, vel, geo):
+def reconstruct_aperture(rho, vel, geo, *, return_absolute=False):
     field = np.concatenate((rho[None], vel), axis=0)
     x, y, a, b = geo["ix"], geo["iy"], geo["a"], geo["b"]
     ray = (
@@ -115,12 +165,24 @@ def reconstruct_aperture(rho, vel, geo):
         + field[:, y + 1, x + 1] * a * b
     )
     estimates = []
+    absolute = []
     for coeff in (LINEAR, QUADRATIC):
         wall = np.einsum("cajk,k->caj", ray, coeff)
         un = np.einsum("caj,ajc->aj", wall[1:], geo["normal"])
         q = np.sum(un * geo["ds"], axis=1)
         mass = np.sum(wall[0] * un * geo["ds"], axis=1)
         estimates.append(np.stack((q, mass), axis=-1))
+        absolute.append(
+            np.stack(
+                (
+                    np.sum(abs(un) * geo["ds"], axis=1),
+                    np.sum(abs(wall[0] * un) * geo["ds"], axis=1),
+                ),
+                axis=-1,
+            )
+        )
+    if return_absolute:
+        return ray, np.array(estimates), np.array(absolute)
     return ray, np.array(estimates)
 
 
@@ -136,6 +198,8 @@ class MeasurementRecorder:
         self.cx = cx
         self.cy = cy
         self.records = []
+        self.sparse_rays = []
+        self.ray_steps = set(old_sample_steps(c).tolist())
         self.rays = [
             ray_geometry(fluid_u, (0.0, 0.0), np.array(c["center"]), c["diameter"] / 2),
             ray_geometry(leaves[1], origin_a, np.array(c["center"]), c["diameter"] / 2),
@@ -179,8 +243,16 @@ class MeasurementRecorder:
         record["cv_inventory"] = []
         record["cv_flux"] = []
         record["momentum_flux_terms"] = []
-        record["arc_ray_values"] = []
-        record["arc_estimates"] = []
+        save_ray = not c.get("dense_root", False) or step in self.ray_steps
+        ray_record = dict(
+            sample_index=len(self.records),
+            arc_ray_values=[],
+            arc_estimates=[],
+            arc_absolute_estimates=[],
+        )
+        if not c.get("dense_root", False):
+            record["arc_ray_values"] = []
+            record["arc_estimates"] = []
         record["nearest_mass_flux"] = [
             nearest_flux[k] for k in ("uniform_adapter", "amr_api")
         ]
@@ -203,9 +275,18 @@ class MeasurementRecorder:
                 rtol=2e-10,
             )
             record["momentum_flux_terms"].append(terms)
-            ray, estimates = reconstruct_aperture(rho, vel, self.rays[idx])
-            record["arc_ray_values"].append(ray)
-            record["arc_estimates"].append(estimates)
+            if save_ray:
+                ray, estimates, absolute = reconstruct_aperture(
+                    rho, vel, self.rays[idx], return_absolute=True
+                )
+                ray_record["arc_ray_values"].append(ray)
+                ray_record["arc_estimates"].append(estimates)
+                ray_record["arc_absolute_estimates"].append(absolute)
+                if not c.get("dense_root", False):
+                    record["arc_ray_values"].append(ray)
+                    record["arc_estimates"].append(estimates)
+        if c.get("dense_root", False) and save_ray:
+            self.sparse_rays.append(ray_record)
         self.records.append(record)
         if step // s**2 in (0, 100, 800, 1600, 3200) and step % s**2 == 0:
             fields = {}
@@ -223,7 +304,11 @@ class MeasurementRecorder:
                 fields[name + "_origin_reference"] = np.asarray(origin) / s
                 fields[name + "_h_reference"] = np.array(h / s)
             file = self.output / f"fields_s{s}_tref{step//s**2}.npz"
-            np.savez_compressed(file, **fields)
+            staged = file.with_name("." + file.name + ".pending.npz")
+            np.savez_compressed(staged, **fields)
+            install_evidence(
+                {file: staged}, EVIDENCE_LIMIT if c.get("dense_root", False) else None
+            )
             self.snapshots.append(
                 dict(
                     file=file.name,
@@ -236,19 +321,32 @@ class MeasurementRecorder:
     def flush(self):
         path = self.output / f"samples_s{self.c['scale']}.npz"
         arrays = {k: np.array([r[k] for r in self.records]) for k in self.records[0]}
+        if self.c.get("dense_root", False):
+            arrays["ray_sample_indices"] = np.array(
+                [r["sample_index"] for r in self.sparse_rays], dtype=int
+            )
+            for key in ("arc_ray_values", "arc_estimates", "arc_absolute_estimates"):
+                arrays[key] = np.array([r[key] for r in self.sparse_rays])
         for i, g in enumerate(self.rays):
             for k in ("normal", "ds", "positions"):
                 arrays[f"ray_{i}_{k}"] = g[k]
         arrays["ray_offsets_fine_cells"] = OFFSETS
         arrays["linear_coefficients"] = LINEAR
         arrays["quadratic_coefficients"] = QUADRATIC
-        np.savez_compressed(path, **arrays)
+        staged = path.with_name("." + path.name + ".pending.npz")
+        np.savez_compressed(staged, **arrays)
         metadata = dict(
-            schema=1,
+            schema=2 if self.c.get("dense_root", False) else 1,
             case=self.c,
+            dense_root=bool(self.c.get("dense_root", False)),
+            ray_sample_index_key=(
+                "ray_sample_indices" if self.c.get("dense_root", False) else None
+            ),
+            raw_absolute_throughput_definition="Cumulative sum of abs(blowing NET dm)+abs(suction NET dm) per root interval; unresolved within-aperture/substep reversals; not abs(window net).",
+            absolute_throughput_definition="Quadrature integral of abs(reconstructed u.n), abs(reconstructed rho*u.n); not abs(signed integral), not prescribed nominal, not exact wall truth.",
             reference_conversion=factors(self.c["scale"]),
             sample_file=path.name,
-            sample_sha256=hashlib.sha256(path.read_bytes()).hexdigest(),
+            sample_sha256=hashlib.sha256(staged.read_bytes()).hexdigest(),
             sample_count=len(self.records),
             snapshots=self.snapshots,
             level_stats_columns=[
@@ -280,6 +378,10 @@ class MeasurementRecorder:
             ],
             warning="Extrapolation and nearest-node values are estimator-dependent, not exact wall values. Fields/ray arrays use root-lattice units except explicitly named reference arrays.",
         )
-        (self.output / f"samples_s{self.c['scale']}.json").write_text(
-            json.dumps(metadata, indent=2) + "\n"
+        meta_path = self.output / f"samples_s{self.c['scale']}.json"
+        meta_staged = meta_path.with_name("." + meta_path.name + ".pending")
+        meta_staged.write_text(json.dumps(metadata, indent=2) + "\n")
+        install_evidence(
+            {path: staged, meta_path: meta_staged},
+            EVIDENCE_LIMIT if self.c.get("dense_root", False) else None,
         )

tests/validation/run_jet_circle_cv.py

@@ -18,11 +18,24 @@ CY = np.array([0, 0, 0, 1, -1, 1, -1, -1, 1])
 W = np.array([4 / 9] + [1 / 9] * 4 + [1 / 36] * 4)
 
 
-def case(scale=1, re=100, steps=800):
+def case(scale=1, re=100, steps=800, ddf_shifting=False):
     u = 0.03 / scale
     d = 24.084 * scale
     return dict(
         scale=scale,
+        ddf_shifting=bool(ddf_shifting),
+        storage_variant=dict(
+            name="shifted_fp32" if ddf_shifting else "unshifted_fp32",
+            device_populations="g=f-w" if ddf_shifting else "f",
+            device_dtype="float32",
+            equilibrium_path=(
+                "compute_feq_from_momentum" if ddf_shifting else "compute_feq"
+            ),
+            initial_physical_seed="historical FP32-rounded equilibrium, identical in both variants",
+            downloaded_populations="physical f, decoded by public stepper.download()",
+            downloaded_dtype="float64" if ddf_shifting else "float32",
+            raw_exchange="physical reaction and mass; no host weight correction",
+        ),
         re_d=re,
         beta=0.1,
         nx=192 * scale,
@@ -226,7 +239,7 @@ def run_case(c, ctx, output, provenance, continuation=False):
         collision="SRT",
         streaming="double_buffer",
         store_precision="FP32",
-        ddf_shifting=False,
+        ddf_shifting=c.get("ddf_shifting", False),
         les_enabled=False,
     )
     fine_cfg = LBMConfig(
@@ -237,7 +250,7 @@ def run_case(c, ctx, output, provenance, continuation=False):
         collision="SRT",
         streaming="double_buffer",
         store_precision="FP32",
-        ddf_shifting=False,
+        ddf_shifting=c.get("ddf_shifting", False),
         les_enabled=False,
     )
     obj = SimObject(
@@ -350,6 +363,10 @@ def run_case(c, ctx, output, provenance, continuation=False):
             for name in ("uniform_adapter", "amr_api")
         }
 
+        from jet_cv_measurements import absolute_aperture_interval_mass
+
+        cumulative_absolute_aperture_mass = np.zeros(2)
+
         def capture(step):
             ug = uniform.download()
             levels = amr.download()
@@ -381,6 +398,10 @@ def run_case(c, ctx, output, provenance, continuation=False):
                         for name in ("uniform_adapter", "amr_api")
                     ]
                 )
+                if c.get("dense_root", False):
+                    recorder.records[-1][
+                        "raw_absolute_aperture_mass_cumulative"
+                    ] = cumulative_absolute_aperture_mass.copy()
             return data, error, aperture
 
         previous, _, previous_aperture = capture(0)
@@ -427,7 +448,17 @@ def run_case(c, ctx, output, provenance, continuation=False):
             np.testing.assert_allclose(
                 ar.sum(axis=0)[:2], obs["total"]["impulse"], atol=1e-13, rtol=1e-13
             )
-            for name, rr in [("uniform_adapter", ur), ("amr_api", ar)]:
+            np.testing.assert_allclose(
+                ar[:, 3].sum(), obs["total"]["mass_into_fluid"], atol=1e-13, rtol=1e-13
+            )
+            for backend_index, (name, rr) in enumerate(
+                [("uniform_adapter", ur), ("amr_api", ar)]
+            ):
+                cumulative_absolute_aperture_mass[
+                    backend_index
+                ] += absolute_aperture_interval_mass(
+                    rr, masks["blowing"], masks["suction"]
+                )
                 for k, mask in masks.items():
                     subtotal = rr[mask].sum(axis=0)
                     accumulated[name][k] += subtotal
@@ -519,8 +550,12 @@ def run_case(c, ctx, output, provenance, continuation=False):
                     physics_validation_pass=False,
                     precision="FP32 evolution, FP64 host diagnostics",
                 )
-                (output / f"re{c['re_d']}_s{c['scale']}.json").write_text(
-                    json.dumps(payload, indent=2) + "\n"
+                from jet_cv_measurements import write_evidence_text
+
+                write_evidence_text(
+                    output / f"re{c['re_d']}_s{c['scale']}.json",
+                    json.dumps(payload, indent=2) + "\n",
+                    dense_root=c.get("dense_root", False),
                 )
                 print(
                     json.dumps(
@@ -552,16 +587,33 @@ def main():
         action="store_true",
         help="Fixed measurement campaign: Re100, scale1/2, 3200 reference steps",
     )
+    p.add_argument(
+        "--dense-root",
+        action="store_true",
+        help="Bounded ref0..800 scalar diagnostics at every root interval; sparse rays/fields",
+    )
+    p.add_argument(
+        "--shifted-fp32",
+        action="store_true",
+        help="Opt-in shifted numerical-representation control for --dense-root; historical unshifted default unchanged",
+    )
     p.add_argument("--plan", action="store_true")
     a = p.parse_args()
+    if a.shifted_fp32 and not a.dense_root:
+        p.error("--shifted-fp32 requires the bounded --dense-root campaign")
     if a.steps <= 0 or a.steps > 3200:
         p.error("steps must be 1..3200")
     if a.re == 400 and not a.allow_re400:
         p.error("Re400 requires explicit conditional-campaign opt-in")
+    if a.dense_root:
+        if a.continuation or a.scale not in (1, 2) or a.re != 100 or a.steps != 800:
+            p.error(
+                "--dense-root requires --re 100 --scale 1/2 --steps 800, without --continuation"
+            )
     if a.continuation:
         if a.scale not in (1, 2) or a.re != 100 or a.steps != 3200:
             p.error("--continuation requires --re 100 --scale 1/2 --steps 3200")
-    c = case(a.scale, a.re, a.steps)
+    c = case(a.scale, a.re, a.steps, ddf_shifting=a.shifted_fp32)
     if a.continuation:
         from jet_cv_measurements import sample_steps, factors
 
@@ -575,12 +627,28 @@ def main():
             post_ready_budget_seconds=1200 if a.scale == 1 else 5400,
             runtime_estimate_seconds=[210, 500] if a.scale == 1 else [900, 3200],
         )
+    if a.dense_root:
+        from jet_cv_measurements import sample_steps, old_sample_steps, factors
+
+        c.update(
+            dense_root=True,
+            measurement_schema=2,
+            reference_conversion=factors(a.scale),
+            post_ready_budget_seconds=1200 if a.scale == 1 else 5400,
+            runtime_estimate_seconds=[60, 240] if a.scale == 1 else [600, 2400],
+        )
+        c.update(
+            sample_count=len(sample_steps(c)),
+            sparse_ray_sample_count=len(old_sample_steps(c)),
+            scalar_sample_dt_reference=1 / a.scale**2,
+            raw_interval_dt_reference=1 / a.scale**2,
+        )
     plan = dict(
         case=c,
         required_revision=PIN,
         required_free_device_mib=1024,
         estimated_peak_device_mib=128 * a.scale**2,
-        estimated_output_mib=64 if a.continuation else 2,
+        estimated_output_mib=32 if a.dense_root else (64 if a.continuation else 2),
         max_steps=c["steps"],
         gpu_launch_owner="root; this runner does not select or reserve a GPU",
     )
@@ -612,7 +680,9 @@ def main():
             free, total = cuda.mem_get_info()
         if free < max(1024, 2 * plan["estimated_peak_device_mib"]) * 1024**2:
             raise RuntimeError("Insufficient free device memory")
-        run_case(c, ctx, a.output, provenance, continuation=a.continuation)
+        run_case(
+            c, ctx, a.output, provenance, continuation=a.continuation or a.dense_root
+        )
     print(
         json.dumps(
             dict(status="COMPLETE_DIAGNOSTIC_ONLY", physics_validation_pass=False)

tests/validation/taylor_green/check_gpu_explosion.py

@@ -1,145 +0,0 @@
-"""Explicit tiny CUDA explosion parity and constructor/one-interval smoke check.
-
-No installs or device selection defaults. Run only on the authorized GPU host.
-This is operator/runtime evidence, not Taylor-Green accuracy qualification.
-"""
-import argparse
-import hashlib
-import json
-import os
-from pathlib import Path
-import sys
-
-import numpy as np
-
-
-def check_case(ctx, shifted, projected):
-    import pycuda.driver as cuda
-    from CelerisLab.config import LBMConfig
-    from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
-    from CelerisLab.lbm.refinement_tables import projected_explosion_tables
-    from CelerisLab.lbm.refinement import ACTIVE
-    from test_projected_explosion import reconstruct, W
-
-    cfg = LBMConfig(nx=16, ny=16, viscosity=.1, velocity=0,
-                    ddf_shifting=shifted)
-    cfg.validate()
-    options = {"explosion": "projected_linear"} if projected else {}
-    # Omitting explosion in the constant case tests the production default.
-    with FixedPatchStepper(cfg, context=ctx, refinements=[[(4, 4, 12, 12)]],
-                           level_blocks=[(4, 4), (8, 4)], **options) as run:
-        parent, child = run._levels
-        pair = run._pairs[0]
-        patch = run.hierarchy.pairs[0]
-        coarse_cells = int(np.prod(parent.layout.local_shape))
-        fine_cells = int(np.prod(child.layout.local_shape))
-        assert coarse_cells != fine_cells
-        assert pair.parent_cells == coarse_cells
-        assert len(pair.interpolation) == (2 if projected else 0)
-        donor = patch.ghost_coarse.astype(np.intp)
-        # Random variation must exercise the gradient, not just a constant field.
-        c = np.random.default_rng(11).normal(0, .001, (9, 16, 16))
-        if not shifted:
-            c += W[:, None, None]
-        c = c.astype(np.float32)
-        legal = np.zeros(256, bool)
-        legal[donor] = True
-        c.reshape(9, -1)[:, ~legal] = np.nan
-        coarse = np.full(np.prod(parent.layout.blocked_shape), np.nan, np.float32)
-        for q in range(9):
-            coarse[parent.layout.slot(parent.layout.owner.ravel()[legal], q)] = c.reshape(9, -1)[q, legal]
-        fine = np.full(np.prod(child.layout.blocked_shape), np.nan, np.float32)
-        cuda.memcpy_htod(parent.current, coarse)
-        cuda.memcpy_htod(child.current, fine)
-        calls = []
-        real_hop = run._hop
-        try:
-            # Execute the production explosion launch exactly; isolate it from
-            # subsequent streaming/collision without replacing the kernel/tables.
-            run._hop = lambda *args, **kwargs: calls.append(args)
-            run._advance(1, parent.current, block=run._block, stream=run._stream)
-            run._stream.synchronize()
-        finally:
-            run._hop = real_hop
-        assert calls == [(1, child.current, child.spare), (1, child.spare, child.current)]
-        cuda.memcpy_dtoh(fine, child.current)
-        targets = child.layout.owner.ravel()[patch.ghost_fine]
-        slots = np.array([child.layout.slot(targets, q) for q in range(9)])
-        actual = fine[slots]
-        neighbors, offsets = projected_explosion_tables(patch, parent.layout)
-        inverse = np.full(parent.layout.nblocks * coarse_cells, -1, dtype=int)
-        inverse[parent.layout.owner.ravel()] = np.arange(256)
-        constant = c.reshape(9, -1)[:, donor].astype(float)
-        expected = reconstruct(c.astype(float), patch, inverse[neighbors], offsets)[0] if projected else constant
-        assert np.isfinite(actual).all()
-        if projected:
-            assert np.max(np.abs(expected - constant)) > 1e-5
-        np.testing.assert_allclose(actual, expected, atol=6e-8, rtol=2e-7)
-        untouched = np.ones(fine.size, bool)
-        untouched[slots.ravel()] = False
-        assert np.isnan(fine[untouched]).all(), "explosion wrote outside target slots"
-        error = float(np.max(np.abs(actual - expected)))
-
-        # Restore ordinary execution and exercise real constructor wiring,
-        # stream-only ghost hops, restriction, halo copies and collision once.
-        initial = [np.broadcast_to(W[:, None, None], (9, *k.shape)).copy()
-                   for k in run.hierarchy.kinds]
-        run.initialize(*initial)
-        run.step(1)
-        for values, kind in zip(run.download(), run.hierarchy.kinds):
-            active = values[:, kind == ACTIVE]
-            assert np.isfinite(active).all()
-            np.testing.assert_allclose(active, np.broadcast_to(W[:, None], active.shape),
-                                       atol=3e-7, rtol=1e-6)
-        assert run.step_count == 1
-        return dict(shifted=shifted, explosion="projected_linear" if projected else "constant_default",
-                    parent_stride=coarse_cells, child_stride=fine_cells,
-                    compared_slots=int(actual.size), max_abs_error=error,
-                    poisoned_nondonors=True, untouched_slots_preserved=True, equilibrium_step=True)
-
-
-def main():
-    parser = argparse.ArgumentParser(description=__doc__)
-    parser.add_argument("--repo", type=Path, required=True)
-    parser.add_argument("--output", type=Path, required=True)
-    parser.add_argument("--run-gpu", action="store_true", help="Explicitly opt into GPU execution")
-    args = parser.parse_args()
-    selected = os.environ.get("CUDA_VISIBLE_DEVICES", "")
-    if not args.run_gpu or not selected or "," in selected or selected == "-1":
-        parser.error("--run-gpu and exactly one explicit physical CUDA_VISIBLE_DEVICES are required")
-    if args.output.exists():
-        parser.error("output already exists; choose a new receipt path")
-    root = args.repo.resolve()
-    sys.path[:0] = [str(root / "src"), str(root / "tests/unit"),
-                   str(root / "tests/validation/taylor_green")]
-    import CelerisLab
-    if Path(CelerisLab.__file__).resolve() != root / "src/CelerisLab/__init__.py":
-        raise RuntimeError("Wrong imported source checkout")
-    from run_case import provenance
-    report = dict(source=provenance(root),
-                  numpy_oracle_sha256=hashlib.sha256((root / "tests/unit/test_projected_explosion.py").read_bytes()).hexdigest(),
-                  cuda_visible_devices=selected,
-                  logical_device=0, completed=False, cases=[])
-    args.output.parent.mkdir(parents=True, exist_ok=True)
-    def save():
-        args.output.write_text(json.dumps(report, indent=2, allow_nan=False) + "\n")
-    save()
-    try:
-        from CelerisLab.cuda.context import CudaContext
-        with CudaContext(0) as ctx:
-            report["device"] = ctx.device.name()
-            for shifted in (False, True):
-                for projected in (False, True):
-                    report["cases"].append(check_case(ctx, shifted, projected))
-                    save()
-        report["completed"] = True
-        save()
-    except Exception as exc:
-        report["error"] = f"{type(exc).__name__}: {exc}"
-        save()
-        raise
-    print(json.dumps(report, indent=2))
-
-
-if __name__ == "__main__":
-    main()

tests/validation/taylor_green/cpu_explosion_summary.json

@@ -1,101 +0,0 @@
-{
-  "schema": 1,
-  "implementation_revision": "43c527b",
-  "historical_baseline_revision": "9805383f9c50d867e8555471baf5dbda8333a652",
-  "scope": "Re100 depth2; FP64 CPU arrival_f macro diagnostics, not candidate CUDA execution or performance",
-  "tstar": 0.48,
-  "cases": [
-    {
-      "n": 32,
-      "historical_gpu_report": "CelerisLab_9805383_tg_matched_re100_n32_d2.json",
-      "historical_gpu_sha256": "6e84837781f2c4eb001721d6f828ea0d927b255fd4c46efdef443a9ae59af730",
-      "constant_max_abs_cpu_gpu_metric_difference": {
-        "velocity_l2": 6.928351422436085e-08,
-        "pressure_l2": 3.5717157327030957e-06,
-        "energy_relative": 2.1743219684455806e-07
-      },
-      "runs": [
-        {
-          "explosion": "constant",
-          "final": {
-            "velocity_l2": 0.03249421441793122,
-            "pressure_l2": 0.08986961517973119,
-            "energy_relative": 0.044970292392050686,
-            "amplitude_bias": -0.02307193176461997,
-            "orthogonal_residual_l2": 0.022881432107435998,
-            "mass_signed_relative": -7.616129948928574e-14,
-            "momentum_change": [
-              -3.26405569239796e-14,
-              -2.7000103541841014e-14
-            ]
-          },
-          "max_abs_mass_relative": 7.616129948928574e-14,
-          "max_abs_momentum_change": 3.608224830031759e-14
-        },
-        {
-          "explosion": "projected_linear",
-          "final": {
-            "velocity_l2": 0.015226371368103809,
-            "pressure_l2": 0.048560123195715434,
-            "energy_relative": 0.02285059695090361,
-            "amplitude_bias": -0.011601431258836059,
-            "orthogonal_residual_l2": 0.009861499773660909,
-            "mass_signed_relative": -7.616129948928574e-14,
-            "momentum_change": [
-              -4.030109579389318e-14,
-              -2.3444787777826548e-15
-            ]
-          },
-          "max_abs_mass_relative": 7.616129948928574e-14,
-          "max_abs_momentum_change": 4.030109579389318e-14
-        }
-      ]
-    },
-    {
-      "n": 64,
-      "historical_gpu_report": "CelerisLab_9805383_tg_re100_n64_d2.json",
-      "historical_gpu_sha256": "251f743d6adc33dff1d829db91fcc5f982c9ac299c177eb2c94e313249b25b9f",
-      "constant_max_abs_cpu_gpu_metric_difference": {
-        "velocity_l2": 1.7330029265794078e-08,
-        "pressure_l2": 8.634212284933107e-07,
-        "energy_relative": 8.74280938097105e-08
-      },
-      "runs": [
-        {
-          "explosion": "constant",
-          "final": {
-            "velocity_l2": 0.01131134574353476,
-            "pressure_l2": 0.03550495463085949,
-            "energy_relative": 0.012514764591008386,
-            "amplitude_bias": -0.006336459299048802,
-            "orthogonal_residual_l2": 0.009369942693596348,
-            "mass_signed_relative": -1.412203687323199e-13,
-            "momentum_change": [
-              -2.1049828546892968e-13,
-              -1.632548263241773e-14
-            ]
-          },
-          "max_abs_mass_relative": 1.412203687323199e-13,
-          "max_abs_momentum_change": 2.1049828546892968e-13
-        },
-        {
-          "explosion": "projected_linear",
-          "final": {
-            "velocity_l2": 0.004930396085095945,
-            "pressure_l2": 0.016580125701009935,
-            "energy_relative": 0.0056291365148441195,
-            "amplitude_bias": -0.0028419317091766816,
-            "orthogonal_residual_l2": 0.004028924138812507,
-            "mass_signed_relative": -1.4099832412739488e-13,
-            "momentum_change": [
-              -1.709743457922741e-13,
-              7.568945470382005e-14
-            ]
-          },
-          "max_abs_mass_relative": 1.4099832412739488e-13,
-          "max_abs_momentum_change": 1.709743457922741e-13
-        }
-      ]
-    }
-  ]
-}

tests/validation/taylor_green/run_case.py

@@ -1,4 +1,4 @@
-"""One bounded GPU case, with explicit opt-in experimental explosion."""
+"""One bounded GPU case. No production edits, plots, checkpoints or field dumps."""
 
 import argparse
 import hashlib
@@ -75,10 +75,7 @@ def main():
     p.add_argument("--n", type=int, choices=[32, 64, 128, 256], required=True)
     p.add_argument("--depth", type=int, choices=[1, 2, 3], required=True)
     p.add_argument("--output", type=Path, required=True)
-    p.add_argument("--explosion", choices=["constant", "projected_linear"], default="constant")
     args = p.parse_args()
-    if args.explosion != "constant" and args.depth != 2:
-        p.error("projected_linear requires depth 2")
     if not os.environ.get("CUDA_VISIBLE_DEVICES"):
         p.error("Explicit physical CUDA_VISIBLE_DEVICES selection required")
     repo = args.repo.resolve()
@@ -102,8 +99,6 @@ def main():
     report = dict(
         schema=1,
         source=provenance(repo),
-        historical_baseline_revision="9805383f9c50d867e8555471baf5dbda8333a652",
-        explosion=args.explosion,
         n=args.n,
         re_L=args.re,
         depth=args.depth,
@@ -137,7 +132,6 @@ def main():
                 coarse_boundary="periodic",
                 refinements=refinements,
                 level_blocks=[(8, 8)] * args.depth,
-                explosion=args.explosion,
             )
         with run:
             run.initialize(

tests/validation/taylor_green/run_cpu_explosion.py

@@ -1,85 +0,0 @@
-"""Diagnostic-only rectangular adapter to CPUREF-0.3. No production changes."""
-import sys,json,time,argparse,hashlib
-from pathlib import Path
-import numpy as np
-parser=argparse.ArgumentParser(description=__doc__)
-parser.add_argument('--repo',type=Path,required=True,help='Candidate checkout at 43c527b')
-parser.add_argument('--n',type=int,choices=(32,64),default=32)
-parser.add_argument('--output-dir',type=Path,required=True)
-args=parser.parse_args()
-ROOT=args.repo.resolve()
-sys.path[:0]=[str(ROOT/'src'),str(ROOT/'tests/validation/multigrid'),str(ROOT/'tests/validation/taylor_green'),str(ROOT/'tests/unit')]
-import cpu_reference as cpu
-import oracle
-from run_case import provenance
-from test_projected_explosion import reconstruct
-from CelerisLab.config import LBMConfig
-from CelerisLab.lbm.refinement import FixedPatchSet2D
-from CelerisLab.lbm.refinement_tables import projected_explosion_tables
-from CelerisLab.lbm.blocks import BlockLayout2D
-class RectSpec(cpu.Spec):
-    @property
-    def n_fine(self): return 2*self.n_coarse
-
-def setup(n):
-    nu=oracle.U*n/100
-    cfg=LBMConfig(nx=n,ny=n,viscosity=nu)
-    bounds=oracle.rectangles(n,2)[0]
-    p=FixedPatchSet2D.from_config(cfg,bounds)
-    spec=RectSpec(n_coarse=n,patch_lo=n//8,patch_hi=7*n//8,halo=n//4,tau_coarse=p.tau_coarse,tau_fine=p.tau_fine,collision='srt')
-    q,k=np.divmod(p.coalesce_coarse.astype(int),n*n); y,x=np.divmod(k,n)
-    rows=np.column_stack((q,y,x)); links=[]
-    for qi,yi,xi in rows:
-        sx,sy=xi-cpu.CX[qi],yi-cpu.CY[qi]
-        if not p.covered[sy,sx]: links.append(cpu.CornerLink(f'{xi}_{yi}_q{qi}',int(qi),(int(sx),int(sy)),(int(xi),int(yi))))
-    assert len(links)==8
-    maps=cpu.Maps(~p.covered,p.covered,p.fine_kind==0,p.fine_kind==1,p.fine_kind==2,p.fine_kind,np.arange(2*n)//2,(np.arange(2*n)+.5)/2,rows,tuple(links))
-    # Verify the CPU coalescence child locations equal actual production tables.
-    children=np.array([[qi*(2*n)**2+(2*yi+dy)*2*n+2*xi+dx for dy in (0,1) for dx in (0,1)] for qi,yi,xi in rows])
-    np.testing.assert_array_equal(children,p.coalesce_fine)
-    layout=BlockLayout2D.build(n,n,8,8,periodic=(True,True))
-    neighbors,offsets=projected_explosion_tables(p,layout)
-    inverse=np.full(layout.nblocks*np.prod(layout.local_shape),-1,dtype=int)
-    inverse[layout.owner.ravel()]=np.arange(n*n)
-    neighbors=inverse[neighbors]
-    def explode(gc,cv,gf,fv,m,writes):
-        # Only actual transfer footprint, including inner and outer ghost slabs.
-        gy,gx=np.divmod(p.ghost_fine.astype(int),2*n)
-        if mode[0]=='constant': values=gc.reshape(9,-1)[:,p.ghost_coarse]
-        else: values=reconstruct(gc,p,neighbors,offsets)[0]
-        assert np.isfinite(values).all()
-        gf[:,gy,gx]=values; fv[:,gy,gx]=True; writes[gy,gx]+=1
-    mode=['constant']; cpu._explode=explode
-    return spec,maps,p,mode,nu
-
-def run(n,which):
-    spec,m,p,mode,nu=setup(n); mode[0]=which
-    masks=[m.coarse_active,m.fine_active]
-    arrays=[oracle.initial(mask.shape,l,n,nu) for l,mask in enumerate(masks)]
-    for a,mask in zip(arrays,masks):a[:,~mask]=np.nan
-    s=cpu.State(*arrays); rows=[]; initmass=None; momentum0=None; start=time.monotonic()
-    for tstar in oracle.TIMES:
-        step=round(tstar*n/oracle.U)
-        while s.tick<2*step: s,_=cpu.advance_coarse_once(spec,m,s)
-        arrays=[s.coarse_f,s.fine_f]
-        row=oracle.metrics(arrays,masks,step,n,nu,initmass)
-        if initmass is None:initmass=row['mass']
-        dot=den=res=0.; parts=[]; momentum=np.zeros(2)
-        for l,(a,mask) in enumerate(zip(arrays,masks)):
-            v=a[:,mask];rho=v.sum(0);u=cpu.CX@v/rho;w=cpu.CY@v/rho
-            x,y=oracle.centers(mask.shape,l);ex,ey,_=oracle.exact(x[mask],y[mask],step,n,nu);vol=.25**l
-            dot+=vol*np.sum(u*ex+w*ey);den+=vol*np.sum(ex*ex+ey*ey)
-            parts.append((u,w,ex,ey,vol));momentum+=vol*np.array([np.sum(cpu.CX@v),np.sum(cpu.CY@v)])
-        amp=dot/den
-        res=sum(vol*np.sum((u-amp*ex)**2+(w-amp*ey)**2) for u,w,ex,ey,vol in parts)
-        if momentum0 is None:momentum0=momentum.copy()
-        row.update(step=step,exact_energy_quadrature_factor=den/(2*n*n*row['exact_energy']),amplitude_ratio=amp,amplitude_bias=amp-1,orthogonal_residual_l2=float(np.sqrt(res/den)),energy_signed_relative=row['energy']/row['exact_energy']-1,mass_signed_relative=row['mass']/initmass-1,momentum=momentum.tolist(),momentum_change=(momentum-momentum0).tolist())
-        rows.append(row)
-    return dict(n=n,explosion=which,seconds=time.monotonic()-start,rows=rows)
-if __name__=='__main__':
-    a=args
-    a.output_dir.mkdir(parents=True,exist_ok=True)
-    result=dict(implementation_revision='43c527b',source=provenance(ROOT),diagnostic_input_sha256={name:hashlib.sha256((ROOT/name).read_bytes()).hexdigest() for name in ('tests/validation/multigrid/cpu_reference.py','tests/unit/test_projected_explosion.py')},precision='FP64 physical populations',phase='arrival_f; macro moments invariant under collision, initialization equilibrium',runs=[])
-    out=a.output_dir/f'n{a.n}_ab.json'
-    for mode in ('constant','projected_linear'):
-        r=run(a.n,mode);result['runs'].append(r);out.write_text(json.dumps(result,indent=2)+'\n');print(json.dumps(r),flush=True)