docs: drop 'rigorous' wording for the window tolerance (no threshold change)

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

docs/development/les_normalization_candidate.md

@@ -1,70 +0,0 @@
-# Minimal Smagorinsky normalization integration candidate
-
-This local candidate changes only the two closed-form D2Q9/D3Q19 coefficients
-from 18 to 18 sqrt(2), plus their convention comments. The intended convention
-is |S| = sqrt(2 S:S), with self-consistent effective tau. No Cs default,
-iterative formula/cap/iteration count, collision dispatch, AMR guard, inlet,
-channel-parity rule, or unrelated source file changes.
-
-The resulting turbulence_smag.cuh SHA256 is
-`a3db46c2e7d27530b2d1436e65cdbed0b72cc7e047bfbeef2a640877271b51aa`,
-exactly the independently reviewed temporary correction used in the paired
-scalar and end-to-end flow experiments. Its uncorrected baseline remains at
-validation checkpoint `9805383f9c50d867e8555471baf5dbda8333a652`.
-
-## Evidence and limits
-
-The separately audited scalar-helper GPU diagnostic on the cd3 source campaign
-covered 48 cases per variant. All 32 paired iterative results were bit-identical;
-maximum patched closed-form error was approximately 1.19e-7, versus 0.06454 for
-the baseline against the same standard convention. This is scalar evidence,
-not full-collision or turbulence qualification.
-
-The independently reviewed wall-bounded, x-periodic SRT/shifted-FP32/double-buffer
-shear experiment distinguished baseline from corrected source under identical
-predeclared criteria. Baseline Cs=1 SGS decrement errors at the last two
-checkpoints were approximately 28.50% and 28.05%; corrected errors were 0.0628%
-and 0.00612%. Corrected maximum Cs=1 profile error was approximately 0.11344%,
-and the Cs=0 observations were identical between variants. Cs=1 is an amplified
-normalization diagnostic, not a recommended physical coefficient. These data
-qualify this bounded normalization comparison, not general turbulence accuracy.
-
-The current direct-header host regression uses the existing independent
-Hermite tensor construction and quadratic bisection oracle. It covers 120
-parameter subcases across D2/D3, Cs=0/.16, rho=.7/1/1.3, tau0=.56/.8,
-equilibrium, both signs of shear and trace-free extension, including a 3D axis.
-It also checks converged iterative agreement, Cs=0 omega guard bounds, and
-explicitly proves that removing sqrt(2) in temporary headers fails the oracle.
-No stored candidate header or full solver fixture is added.
-
-MRT full-trace/bulk coupling, forcing/Guo correction, other streaming/storage
-paths, three-dimensional turbulent flows, and AMR LES remain uncertified.
-Full-stress closed form and deviatoric/capped iteration intentionally need not
-agree on isotropic or cap-active stresses. Existing guards remain unchanged.
-
-## Run and retain provenance
-
-From this candidate checkout:
-
-    PYTHONPATH=src python tests/unit/test_smagorinsky_normalization.py -v
-    PYTHONPATH=src python tests/unit/test_les_flow_runner.py -v
-
-The flow runner already supports `--variant corrected` against the exact
-reviewed helper hash. Its source-pin tests now exercise that current variant;
-they still verify that the archived baseline variant rejects the corrected
-source and that unrelated source tampering is rejected.
-
-Historical Re100, Jet, matched-DFG and paired-scalar manifests/runners are not
-repinned or weakened. Their baseline-source guards are expected to reject this
-changed production tree. Run those archived source-specific checks on the clean
-9805383 checkpoint; their old GPU evidence does not become evidence for this
-candidate. The new direct-current-header test is the regression for corrected
-production. Do not describe the historical full suite as passing on this tree.
-
-Before publication or main integration, obtain final independent diff review
-and rerun the applicable scalar/flow GPU commands on the exact final candidate
-commit after the executor reconnects. Retain new checkout/import/source-hash
-receipts and compare against the frozen baseline using separate processes.
-A new commit receipt is required even though its corrected header matches the
-previously tested temporary correction. No GPU rerun, publication, or main
-merge is claimed by this candidate preparation.

docs/development/refinement_mrt.md

@@ -0,0 +1,34 @@
+# MRT reuse in blocked/refined stepping (candidate)
+
+Base `9805383`. Opt-in by `LBMConfig(collision="MRT")`; SRT remains the default and its
+kernel statements, including the shifted retained-delta/direct-momentum path, are unchanged
+under the alternate compile branch.
+
+## What changes
+- `blocks.cu` (`BlockStep`), `step/block_channel.cu` (`BlockChannelStep`) and `refinement.cu`
+  (`patch_collide`) gain a compile-time `COLLISION_MODEL == 2` branch that repeats the uniform
+  `collide_dispatch` sequence with zero forcing: `compute_rho_u`, omega clamp to
+  `[OMEGA_COLLISION_MIN, OMEGA_COLLISION_MAX]`, `collide_mrt_no_force`.
+- Kernel and host guards accept SRT or MRT only. TRT, LES, non-FP32 storage, non-double-buffer
+  streaming, and shifted dense/mixed refinement remain rejected.
+- Transfers, subcycling, BFL, raw exchange ledgers and channel closures are untouched. The
+  existing pressure-NEQ outlet selects its MRT branch (unknowns-only NEQ extrapolation) from
+  `COLLISION_MODEL`, so MRT channel parity must use uniform MRT as reference.
+
+## Contracts (existing operator, not redefined here)
+- Shifted storage uses the existing MRT representation (equilibrium moments from `rho - 1`),
+  not the SRT retained-delta/direct-momentum precision. Mass and momentum are conserved
+  moments (relaxation rate 0).
+- Kinetic rates are the operator's fixed `s_e = s_eps = s_q = 1.2` on every level; the shear
+  rate is the existing per-level omega. Under acoustic 2:1 scaling the physical bulk viscosity
+  therefore halves on each finer level. No rate rescaling is applied.
+- No LES, no external forcing (`d_params` force fields stay zero).
+- Shifted Jets stay SRT-only: `body/jets.py` keeps its existing guard. This extension does not
+  provide shifted MRT Jets.
+- `PeriodicBlockStepper` (uniform blocked stepper) remains SRT-only.
+
+## Evidence status
+CPU: `tests/unit/test_refinement_mrt.py` (host guards, source contract). GPU:
+`tests/validation/multigrid/test_refinement_mrt_gpu.py` (uniform-MRT channel parity, AMR MRT
+smoke) and unchanged SRT GPU tests. Strict parity may be subject to cross-kernel FMA
+contraction differences already documented for SRT. Not a DFG/benchmark validation.

docs/validation/refined_dfg_mrt.md

@@ -0,0 +1,71 @@
+# Matched MRT DFG pilot preparation
+
+Preparation only: this test-only extension reuses the existing runner and existing MRT
+operator. GPU admission awaits reviewed MRT capability results; no pilot/full execution
+is implied by CPU preflight. No production, historical SRT manifest, reference data,
+periodic gate, geometry, initialization or observation algorithm is changed here.
+
+## Selection and provenance
+
+Use `--collision MRT` for **both** uniform D80 and refined root D40/fine D80. There are
+no independently selectable collision models per case. Default remains SRT.
+
+MRT verifies `data/refined_dfg/source_manifest_mrt.json`, pinned to reviewed numerical
+source snapshot `3ea131e92c057843034f11a1c574301a2430dfce`, including the shared geometry,
+pressure/cycle helpers and official reference. Test-only commits may change HEAD without
+changing this numerical source pin. Plan, READY, case summaries and comparison separately
+record `checkout_revision` and `checkout_dirty` using existing source-provenance helper;
+the runner hash remains separate. Collision is recorded in every receipt, and full-mode
+validation rejects missing or different collision labels in plan, either summary or either
+READY receipt, in addition to existing source/runner/health/time checks.
+
+The historical SRT manifest remains byte-for-byte unchanged. SRT on the MRT candidate
+checkout deliberately fails that source pin; use a matching historical checkout via
+`--source` for SRT. No fallback or automatic repinning is permitted. Historical SRT
+failure/diagnostic evidence remains attached to its original source and runner hashes;
+it cannot authorize an MRT full run. The existing failed-observation step-reporting
+undercount is not changed in this narrowly scoped selection patch.
+
+## Frozen physical comparison
+
+All settings remain those documented in `refined_dfg.md`: U_lattice=.01, physical mean
+U=1 m/s, D=.1 m, shear ν=.001 m²/s, fixed physical domain/body/wall planes, patch
+(40,24,320,140), shifted FP32 double buffers, no LES/forcing/Jets, finest dx=.00125 m,
+finest dt=.0000125 s, unchanged half-cell frame conversion and axis-specific phase.
+Force interval accumulation, physical leaf mass, pressure interpolation, flux sampling,
+safety thresholds, t*=5 pilot and t*=300 cap are unchanged. Existing accuracy/periodicity
+criteria are imported unchanged. Both MRT cases are fresh runs; old MRT or SRT results
+are not substituted as controls.
+
+Both cases use the **existing MRT pressure-NEQ branch**: only unknown directions q2/6/8
+are extrapolated. This differs from SRT's all-nine .5-damped branch; the retained
+`srt_neq_damp=.5` configuration value is inactive under MRT. Backflow clamp remains off.
+The outlet is still noncanonical relative to exact DFG do-nothing conditions.
+
+The reused MRT operator keeps fixed kinetic rates `s_e=s_eps=s_q=1.2` on all levels;
+only the shear rate follows existing per-level viscosity. Physical bulk relaxation
+therefore differs between coarse and fine levels. Transfer operators and subcycling
+are unchanged. This is the candidate's declared behavior, not a new fully mode-scaled
+MRT transfer scheme. Do not tune rates, the patch or gates after seeing outcomes.
+Capability smoke/parity gates do not establish DFG accuracy or speedup.
+
+## CPU and admission commands
+
+From this checkout:
+
+```sh
+python tests/unit/test_refined_dfg.py
+python tests/validation/run_refined_dfg.py --collision MRT --mode plan --output /tmp/new-mrt-dfg-plan
+```
+
+After root verifies the candidate MRT root/parity/refinement capability results and admits
+an idle device, the paired pilot may use the existing bounds and external watchdog:
+
+```sh
+CUDA_VISIBLE_DEVICES=<admitted-device> timeout --signal=TERM --kill-after=30s <hard-budget>s python tests/validation/run_refined_dfg.py --collision MRT --mode pilot --gpu-gates-approved --max-wall-seconds <soft-budget> --output <new-mrt-pilot-directory>
+```
+
+The flag is caller attestation, not evidence. An unhealthy first case prevents the second.
+Existing strict pilot receipt validation and calibrated wall-time controls remain; full
+mode is not being admitted as part of this preparation. Root must review both healthy
+pilot outcomes and separately decide any full-run admission.

docs/validation/refined_dfg_timing.md

@@ -0,0 +1,58 @@
+# Refined DFG MRT timing harness (measurement only)
+
+`tests/validation/time_refined_dfg.py` times the matched MRT refined/uniform DFG Re100 setup of
+`run_refined_dfg.py` (same configuration, geometry and initial state functions). No numerical,
+physics-gate or production change; no force or accuracy result is produced.
+
+| Variant | Loop | Force reads |
+|---|---|---|
+| A | refined `step(1)` + `read_raw_exchange` per base step | runner-like force-read workload (values discarded) |
+| B | refined `step(1)` only, no host reads (`observe_raw_exchange` stays on) | timing-only |
+| C | refined single `step(N)` | timing-only; ledger keeps the final base interval only |
+| E | uniform `sim.run(<=100)` + `read_force` | runner-like force-read workload (values discarded) |
+| E' | uniform single `sim.run(2N)` | timing-only; at most a whole-window average, different summation order |
+
+- No variant retains a valid force history; A/E omit the runner's interval aggregation and
+  validation, and no variant includes macro/probe/CSV costs. Results are not full benchmark
+  throughput.
+- Status lines: `CONSTRUCTED` (construction/JIT complete) and `WARM_READY` (first warm-up and
+  health check passed) per case.
+- Equal simulated duration: N refined base steps vs 2N uniform D80 steps.
+- Every repeat starts from the same fresh state (refined `initialize`, uniform `snapshot/restore`)
+  and the same warm-up; construction/JIT, warm-up, downloads and parity checks are untimed.
+- The actual stepping stream is synchronized explicitly before and after each timed loop; the
+  wall timer (`time.perf_counter`) is the only source. No CUDA-event span is called kernel time.
+- `--max-wall-seconds` is required for GPU runs (outer timeout, checked between variants and
+  every 100 steps in per-step loops).
+- Final populations and NaN masks must be bitwise identical to the first A repeat (A/B/C) or the
+  first E repeat (E/E'); refined and uniform grids are never compared. Leaf/fluid populations
+  must be finite. Physical health (finite leaf/fluid macros, rho in [0.5, 1.5], max Mach <= 0.2,
+  the existing runner limits) is checked after warm-up and after each timed window. Matching NaN
+  is accepted only as a storage mask; any infinity or level-count mismatch fails parity. All
+  parity/health checks are untimed.
+- Cases run sequentially (uniform, then refined) in one process with startup serialized;
+  variants are interleaved within each repeat. A-B and B-C are measured comparative
+  differences under GPU variability, not exact causal decompositions.
+
+Usage (root selects one verified idle GPU):
+
+    PYTHONPATH=src python -m unittest tests/unit/test_time_refined_dfg.py
+    python tests/validation/time_refined_dfg.py --output <new> --plan
+    CUDA_VISIBLE_DEVICES=<gpu> python tests/validation/time_refined_dfg.py --output <new> --max-wall-seconds <s>
+
+Defaults: 2000 base steps, 200 warm-up base steps, 3 repeats. Profiler capture and the
+`observe_raw_exchange=False` variant are deferred until these measurements exist.
+
+## A versus observation-window ledger (`tests/validation/time_window_dfg.py`)
+
+Measurement-only comparison on the same shifted-MRT refined setup: A (per-base `step(1)` +
+`read_raw_exchange`, summed per 80-base-step interval), W (`start_exchange_window`, `step(80)`,
+one window read) and W1 (window with 80 x `step(1)`, isolating the read pattern from call
+batching). Same fresh state, warm-up and simulated duration; populations must be bitwise equal
+to the first A repeat and per-interval total impulses within the documented engineering
+tolerance, all outside timing. The candidate source is recorded (checkout revision, dirty flag,
+file hashes); frozen refined-DFG manifests/results are not used or changed.
+
+    PYTHONPATH=src python -m unittest tests/unit/test_time_window_dfg.py
+    python tests/validation/time_window_dfg.py --output <new> --plan
+    CUDA_VISIBLE_DEVICES=<gpu> python tests/validation/time_window_dfg.py --output <new> --max-wall-seconds <s>

src/CelerisLab/lbm/channel.py

@@ -1,4 +1,4 @@
-"""Host contract for the blocked root's uniform-SRT channel boundary.
+"""Host contract for the blocked root's uniform SRT/MRT channel boundary.
 
 Boundary ghosts are evolved state, but are not physical fluid leaves. Halo
 ownership remains periodic to reproduce the uniform kernel's wrapped pulls.
@@ -29,7 +29,7 @@ def validate_channel_config(cfg, blocks):
 def validate_channel_transfers(pair, layout, targets):
     """Reject actual exchange/closure overlap, not a padded geometric envelope.
 
-    Explosion donors and restriction receivers must use ordinary bulk SRT
+    Explosion donors and restriction receivers must use ordinary bulk collision
     semantics. Restriction writes receiver-local source slots; these may not
     overwrite the channel's extra immutable-input reads (pressure donor and
     same-node halfway-wall corrections). The slot check deliberately respects

src/CelerisLab/lbm/kernels/blocks.cu

@@ -14,12 +14,15 @@
 #include "core/obs.cuh"
 #include "operators/macro.cuh"
 #include "operators/collision_srt.cuh"
+#if COLLISION_MODEL == 2
+#include "operators/collision_mrt.cuh"
+#endif
 #include "streaming/pull_double_buffer.cuh"
 #include "boundary/curved_boundary.cuh"
 #include "boundary/planar_velocity.cuh"
 
-#if DIM != 2 || NQ != 9 || COLLISION_MODEL != 0 || STORE_PRECISION != 0 || USE_LES || STREAMING_MODEL != 0
-#error "Blocks require FP32 D2Q9 SRT double buffers without LES"
+#if DIM != 2 || NQ != 9 || (COLLISION_MODEL != 0 && COLLISION_MODEL != 2) || STORE_PRECISION != 0 || USE_LES || STREAMING_MODEL != 0
+#error "Blocks require FP32 D2Q9 SRT or MRT double buffers without LES"
 #endif
 
 // Fixed-patch cell kinds, matching refinement.py ACTIVE/INNER_GHOST/OUTER_GHOST.
@@ -141,6 +144,13 @@ extern "C" __global__ void BlockStep(const fpxx* in, fpxx* out, unsigned int hal
                         gx, d_block_wall.nx, d_block_wall.slot_width, d_block_wall.amplitude));
                 else if (gy == d_block_wall.ny - 2) reconstruct_top_reservoir(f, RHO);
             }
+#if COLLISION_MODEL == 2
+            // Uniform MRT sequence (collide_dispatch, zero forcing); SRT path below unchanged.
+            float rho, ux, uy;
+            compute_rho_u(f, rho, ux, uy);
+            collide_mrt_no_force(f, rho, ux, uy,
+                fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega)));
+#else
             float rho, ux, uy, feq[NQ];
             compute_rho_u(f, rho, ux, uy);
 #if USE_DDF_SHIFTING
@@ -153,6 +163,7 @@ extern "C" __global__ void BlockStep(const fpxx* in, fpxx* out, unsigned int hal
             compute_feq(rho, ux, uy, feq);
 #endif
             collide_srt_no_force(f, feq, d_params.omega);
+#endif
         }
     }
     stream_pull_store(k, f, out + base);

src/CelerisLab/lbm/kernels/operators/turbulence_smag.cuh

@@ -5,9 +5,7 @@
 //   nu_t = (C_s * Delta)^2 * |S|,  Delta = 1 (lattice unit)
 //   omega_eff = 1 / (3*(nu + nu_t) + 0.5)
 //
-// |S| = sqrt(2 S_ij S_ij), with Delta = 1 in lattice units.
-// The stress estimate uses self-consistent tau_eff; tau0 is molecular.
-// Closed form uses full stress. Iteration uses deviatoric stress and a cap.
+// |S| is estimated from non-equilibrium second-order moments using tau0.
 // ============================================================================
 
 #ifndef CELERIS_OPERATORS_TURBULENCE_SMAG_CUH
@@ -67,12 +65,12 @@ __device__ __forceinline__ float compute_omega_smag(const float* __restrict__ f,
 
 #if LES_CLOSED_FORM
     // Closed-form Smagorinsky (no iteration):
-    //   ω = 2 / (τ₀ + √(τ₀² + 18·√2·Cs²·√Q / ρ))
+    //   ω = 2 / (τ₀ + √(τ₀² + 18·Cs²·√Q / ρ))
     // where Q = Πneq_ij · Πneq_ij (Frobenius norm squared)
     const float Q = pixx*pixx + piyy*piyy + 2.0f*pixy*pixy;
     const float sqrtQ = sqrtf(fmaxf(Q, 0.0f));
     const float smag_const_sq = LES_CS * LES_CS;
-    const float discriminant = tau0*tau0 + 18.0f * sqrtf(2.0f) * smag_const_sq * sqrtQ / rho_safe;
+    const float discriminant = tau0*tau0 + 18.0f * smag_const_sq * sqrtQ / rho_safe;
     const float tau_eff = 0.5f * (tau0 + sqrtf(fmaxf(discriminant, tau0*tau0)));
     return clamp_omega(1.0f / fmaxf(tau_eff, 0.500001f));
 #else
@@ -142,7 +140,7 @@ __device__ __forceinline__ float compute_omega_smag(const float* __restrict__ f,
                   + 2.0f*(pixy*pixy + pixz*pixz + piyz*piyz);
     const float sqrtQ = sqrtf(fmaxf(Q, 0.0f));
     const float smag_const_sq = LES_CS * LES_CS;
-    const float discriminant = tau0*tau0 + 18.0f * sqrtf(2.0f) * smag_const_sq * sqrtQ / rho_safe;
+    const float discriminant = tau0*tau0 + 18.0f * smag_const_sq * sqrtQ / rho_safe;
     const float tau_eff = 0.5f * (tau0 + sqrtf(fmaxf(discriminant, tau0*tau0)));
     return clamp_omega(1.0f / fmaxf(tau_eff, 0.500001f));
 #else

src/CelerisLab/lbm/kernels/refinement.cu

@@ -11,6 +11,9 @@
 #include "core/obs.cuh"
 #include "operators/macro.cuh"
 #include "operators/collision_srt.cuh"
+#if COLLISION_MODEL == 2
+#include "operators/collision_mrt.cuh"
+#endif
 #include "streaming/pull_double_buffer.cuh"
 #include "boundary/planar_velocity.cuh"
 #include "boundary/curved_boundary.cuh"
@@ -19,8 +22,8 @@ extern "C" {
 #include "step/curved_double.cu"
 }
 
-#if DIM != 2 || NQ != 9 || COLLISION_MODEL != 0 || STORE_PRECISION != 0 || USE_DDF_SHIFTING || USE_LES || STREAMING_MODEL != 0
-#error "Refinement requires unshifted FP32 D2Q9 SRT double buffers without LES"
+#if DIM != 2 || NQ != 9 || (COLLISION_MODEL != 0 && COLLISION_MODEL != 2) || STORE_PRECISION != 0 || USE_DDF_SHIFTING || USE_LES || STREAMING_MODEL != 0
+#error "Refinement requires unshifted FP32 D2Q9 SRT or MRT double buffers without LES"
 #endif
 
 // Enabled only in the validation build; no diagnostic arrays in normal builds.
@@ -33,10 +36,17 @@ __device__ unsigned int incoming_counts[NQ * TOTAL_CELLS];
 #endif
 
 __device__ __forceinline__ void patch_collide(float* f, unsigned long k) {
+#if COLLISION_MODEL == 2
+    float rho, ux, uy;
+    compute_rho_u(f, rho, ux, uy);
+    collide_mrt_no_force(f, rho, ux, uy,
+        fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega)));
+#else
     float rho, ux, uy, feq[NQ];
     compute_rho_u(f, rho, ux, uy);
     compute_feq(rho, ux, uy, feq);
     collide_srt_no_force(f, feq, d_params.omega);
+#endif
 #ifdef CELERIS_REFINEMENT_AUDIT
     ++collision_counts[k];
 #endif

src/CelerisLab/lbm/kernels/step/block_channel.cu

@@ -54,12 +54,20 @@ extern "C" __global__ void BlockChannelStep(
             f[7] = load_ddf(in + base, index_f(k, 8));
         }
         if (!wall && !outlet) {
+#if COLLISION_MODEL == 2
+            // Match uniform MRT arithmetic (collide_dispatch, zero forcing).
+            float rho, ux, uy;
+            compute_rho_u(f, rho, ux, uy);
+            collide_mrt_no_force(f, rho, ux, uy,
+                fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega)));
+#else
             // Match uniform SRT arithmetic, including its shifted equilibrium path.
             float rho, ux, uy, feq[NQ], forcing[NQ] = {};
             compute_rho_u(f, rho, ux, uy);
             compute_feq(rho, ux, uy, feq);
             const float omega = fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega));
             collide_srt(f, feq, forcing, omega);
+#endif
         }
     }
     stream_pull_store(k, f, out + base);

src/CelerisLab/lbm/raw_exchange.py

@@ -7,16 +7,26 @@ import numpy as np
 
 
 def summarize_exchange(rows, body_ids, aperture, *, body_id, level, h,
-                       time_end, window, reference_center):
+                       time_end, window, reference_center, duration=None):
     """Reduce FP64 interval integrals or convert FP32 last-substep raw rows.
 
     dt_l=h_l under acoustic 2:1 refinement. Local momentum j -> base force
     h*j -> impulse h*dt*j; torque adds one h; mass is h**2*dm.
+    window="window" reduces an FP64 observation-window ledger spanning
+    ``duration`` (integer >= 1) base intervals; rates divide by that duration.
     """
-    if window not in ("last_base_step", "last_substep"):
+    if window not in ("last_base_step", "last_substep", "window"):
         raise ValueError("unknown exchange window")
+    if window == "window":
+        if isinstance(duration, bool) or not isinstance(duration, (int, np.integer)) or duration < 1:
+            raise ValueError("window duration must be an integer number of base intervals >= 1")
+    elif duration is not None:
+        raise ValueError("duration applies only to window='window'")
     values = np.array(rows, dtype=np.float64, copy=True)
-    duration = 1.0 if window == "last_base_step" else h
+    if window == "window":
+        duration = float(duration)
+    else:
+        duration = 1.0 if window == "last_base_step" else h
     if window == "last_substep":
         values *= np.array([h*h, h*h, h*h*h, h*h])
     selected = np.asarray(body_ids) == body_id
@@ -31,6 +41,7 @@ def summarize_exchange(rows, body_ids, aperture, *, body_id, level, h,
     return dict(wall=reduce(selected & ~aperture), aperture=reduce(selected & aperture),
                 total=reduce(selected), body_id=int(body_id), level=int(level),
                 time_start=time_end-duration, time_end=time_end, dt=h,
-                sample_count=round(1/h) if window == "last_base_step" else 1,
+                sample_count=(round(duration/h) if window == "window"
+                              else round(1/h) if window == "last_base_step" else 1),
                 window=window, units="base_lattice_per_unit_depth",
                 definition="raw_boundary_exchange", reference_center=tuple(reference_center), valid=True)

src/CelerisLab/lbm/refinement.py

@@ -52,17 +52,17 @@ class FixedPatch2D:
     def from_config(cls, cfg: LBMConfig, bounds) -> "FixedPatch2D":
         """Prepare maps without changing cfg or allocating coarse-sized tables.
 
-        The initial backend contract is D2Q9/SRT/FP32/double_buffer,
+        The backend contract is D2Q9/SRT or MRT/FP32/double_buffer,
         without LES. Two coarse cells of margin keep all halo donors in bounds;
         this does NOT certify that those donors are fluid or periodic.
         """
         checked = replace(cfg)
         checked.validate()
-        if (not checked.is_d2q9 or checked.collision != "SRT"
+        if (not checked.is_d2q9 or checked.collision not in ("SRT", "MRT")
                 or checked.streaming != "double_buffer"
                 or checked.store_precision != "FP32"
                 or checked.les_enabled):
-            raise ValueError("FixedPatch2D requires D2Q9/SRT/FP32/double_buffer without LES")
+            raise ValueError("FixedPatch2D requires D2Q9/SRT or MRT/FP32/double_buffer without LES")
         bounds = tuple(bounds)
         if len(bounds) != 4 or any(isinstance(v, (bool, np.bool_))
                                    or not isinstance(v, Integral) for v in bounds):
@@ -451,7 +451,7 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=(),
     checked = replace(cfg)
     checked.validate()
     if (patch.coarse_shape != (cfg.ny, cfg.nx) or not checked.is_d2q9
-            or checked.collision != "SRT" or checked.streaming != "double_buffer"
+            or checked.collision not in ("SRT", "MRT") or checked.streaming != "double_buffer"
             or checked.store_precision != "FP32" or checked.les_enabled
             or patch.tau_coarse != 3*checked.viscosity+.5):
         raise ValueError("Fine boundary requires the matching fixed-patch configuration")

src/CelerisLab/lbm/refinement_stepper.py

@@ -172,6 +172,9 @@ class FixedPatchStepper:
     base interval or last substep in base-lattice units per unit depth. This
     is boundary momentum reaction, not a validated nozzle/stress load. Dense
     body layouts are not supported by this new diagnostic in this candidate.
+    ``observe_window_exchange=True`` (requires observe_raw_exchange) adds a second
+    per-link FP64 ledger accumulated by the same kernel between
+    start_exchange_window() and read_exchange_window(); disabled by default.
     ``jets`` (body.jets.Jet) use the same edge frame as ``fine_objects`` and
     drive fine cut links selected at their true hits; ``set_jet_amplitude``
     re-uploads the per-link velocity column only.
@@ -216,7 +219,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):
+                 observe_raw_exchange=False, observe_window_exchange=False):
         if HALO != 4:
             raise RuntimeError("refinement.cu requires a four-cell fine halo")
         if coarse_boundary not in ("periodic", "slot_reservoir", "channel"):
@@ -315,6 +318,13 @@ class FixedPatchStepper:
                                       for b, layout in zip(self._boundary_hosts, self._layouts)):
             raise ValueError("raw exchange currently requires blocked body levels")
         self._observe_raw_exchange = observe_raw_exchange
+        if not isinstance(observe_window_exchange, bool):
+            raise ValueError("observe_window_exchange must be bool")
+        if observe_window_exchange and not observe_raw_exchange:
+            raise ValueError("observe_window_exchange requires observe_raw_exchange")
+        self._observe_window = observe_window_exchange
+        self._window_valid = False
+        self._window_start = None
         self._exchange_valid = False
         self._exchange_spacing = tuple(spacings)
         self._exchange_centers = tuple(tuple(tuple(o.geometry.center) for o in objects)
@@ -365,6 +375,8 @@ class FixedPatchStepper:
                             count = len(b.tables["cell"])
                             b.last_raw = self._allocate(count * 4 * 4)
                             b.interval_integral = self._allocate(count * 4 * 8)
+                            if self._observe_window:
+                                b.window_integral = self._allocate(count * 4 * 8)
                             b.accumulate = level.module.get_function("BlockAccumulateExchange")
                         b.kernel = level.module.get_function("BlockCurvedBoundary")
                 if self._level_boundaries[1] is not None:
@@ -559,6 +571,7 @@ class FixedPatchStepper:
                     cuda.memset_d8(self._fine_curved.link_exchange_gpu, 0, self._fine_curved.count * 12)
                 self._step_count = 0
                 self._exchange_valid = False
+                self._window_valid = False
                 self._status = "ready"
             except BaseException:
                 self._status = "failed"
@@ -639,7 +652,8 @@ class FixedPatchStepper:
                 fine_block=fine_block, child_block=child_block,
                 child_rectangles_in_parent=new.child.bounds,
                 child_objects=child_objects, child_jets=jets,
-                observe_raw_exchange=self._observe_raw_exchange)
+                observe_raw_exchange=self._observe_raw_exchange,
+                observe_window_exchange=self._observe_window)
             try:
                 if self._ddf_shifting:
                     candidate._initialize_raw(_StoredPopulations(transferred.g, transferred.shifted))
@@ -771,6 +785,10 @@ class FixedPatchStepper:
             b.accumulate(b.last_raw, b.interval_integral, np.uint32(count),
                          np.float64(h), np.float64(h),
                          grid=((count + self._block[0] - 1) // self._block[0], 1), **kw)
+            if self._window_valid and b.window_integral is not None:
+                b.accumulate(b.last_raw, b.window_integral, np.uint32(count),
+                             np.float64(h), np.float64(h),
+                             grid=((count + self._block[0] - 1) // self._block[0], 1), **kw)
         flags = np.uintp(0) if level.boundary is None else level.boundary.flags
         level.kernel(src, dst, np.uint32(level.layout.halo), flags, level.kind, grid=level.grid, **kw)
 
@@ -804,9 +822,57 @@ class FixedPatchStepper:
                     self._exchange_valid = self._observe_raw_exchange
                 self._status = "ready"
             except BaseException:
+                self._window_valid = False
                 self._status = "failed"
                 raise
 
+    def start_exchange_window(self):
+        """Zero the window ledgers and open a window at the current completed step count."""
+        with self._lock, self.ctx.activate():
+            self._require("ready")
+            if not self._observe_window:
+                raise RuntimeError("window exchange observation is disabled")
+            self._window_valid = False
+            try:
+                for b in self._level_boundaries:
+                    if b is not None and b.window_integral is not None:
+                        cuda.memset_d8(b.window_integral, 0, len(b.tables["cell"]) * 32)
+            except BaseException:
+                self._status = "failed"
+                raise
+            self._window_start = self._step_count
+            self._window_valid = True
+
+    def read_exchange_window(self, body_id, *, level=1):
+        """Non-destructive read of [window start, current step count] in base-lattice units."""
+        with self._lock, self.ctx.activate():
+            self._require("ready")
+            if not self._observe_window:
+                raise RuntimeError("window exchange observation is disabled")
+            if not self._window_valid:
+                raise RuntimeError("no valid exchange window; call start_exchange_window()")
+            duration = self._step_count - self._window_start
+            if duration < 1:
+                raise RuntimeError("exchange window requires at least one completed base interval")
+            if (isinstance(level, bool) or not isinstance(level, Integral)
+                    or not 0 <= level < len(self._levels)):
+                raise ValueError("level must name an existing level")
+            b = self._level_boundaries[level]
+            if (b is None or isinstance(body_id, bool) or not isinstance(body_id, Integral)
+                    or not 0 <= body_id < len(b.host.objects)):
+                raise ValueError("body_id must name a body on this level")
+            rows = np.empty((len(b.tables["cell"]), 4), np.float64)
+            try:
+                cuda.memcpy_dtoh(rows, b.window_integral)
+            except BaseException:
+                self._status = "failed"
+                raise
+            return summarize_exchange(rows, b.tables["body_id"], b.tables["aperture"],
+                                      body_id=body_id, level=level, h=self._exchange_spacing[level],
+                                      time_end=self._step_count, window="window",
+                                      reference_center=self._exchange_centers[level][body_id],
+                                      duration=duration)
+
     def read_raw_exchange(self, body_id, *, level=1, window="last_base_step"):
         """Read raw boundary reaction about its geometric center in the root edge frame.
 
@@ -885,6 +951,7 @@ class FixedPatchStepper:
                     if b is not None:
                         b.flags = b.device = b.kernel = None
                         b.last_raw = b.interval_integral = b.accumulate = None
+                        b.window_integral = None
                 self._levels.clear()
                 self._pairs = []
                 if getattr(self, "_fine_curved", None) is not None:

src/CelerisLab/lbm/refinement_tables.py

@@ -85,6 +85,7 @@ class _BlockedBoundary:
     last_raw: object = None
     interval_integral: object = None
     accumulate: object = None
+    window_integral: object = None  # optional observation-window FP64 ledger
 
 
 def _prepare_blocked_boundary(boundary, layout, cfg, jets, *, _shifted_blocked=False):

tests/unit/test_les_flow_runner.py

@@ -41,9 +41,9 @@ class GuardTests(unittest.TestCase):
         with self.assertRaisesRegex(RuntimeError, "velocity/profile"):
             decode_macros(invalid, w, cx, cy)
 
-    def test_current_corrected_pins(self):
+    def test_baseline_pins(self):
         self.assertEqual(
-            source_guard(ROOT, "corrected")["base_commit"],
+            source_guard(ROOT, "baseline")["base_commit"],
             "cd3a657f9d3833e812c3cca173525f165b48710a",
         )
 
@@ -55,7 +55,20 @@ class GuardTests(unittest.TestCase):
                 dest = root / name
                 dest.parent.mkdir(parents=True, exist_ok=True)
                 shutil.copyfile(ROOT / name, dest)
-            # Current production is the corrected variant; archived baseline pins stay frozen.
+            # Rebuild the reviewed candidate from actual pinned headers, not fixtures.
+            candidate = root / "header_candidate"
+            subprocess.run(
+                [
+                    sys.executable,
+                    str(HERE / "les/prepare.py"),
+                    "--repo",
+                    str(ROOT),
+                    "--output-dir",
+                    str(candidate),
+                ],
+                check=True,
+            )
+            shutil.copyfile(candidate / "patched" / HEADER, root / HEADER)
             with patch(
                 "run_uniform_les_flow.subprocess.check_output",
                 return_value="isolated-fixture\n",

tests/unit/test_raw_exchange_window.py

@@ -0,0 +1,45 @@
+"""Host reduction contract for the optional observation-window exchange ledger."""
+import unittest
+
+import numpy as np
+
+from CelerisLab.lbm.raw_exchange import summarize_exchange
+
+
+def reduce(rows, window, duration=None, h=.5, time_end=10):
+    return summarize_exchange(rows, np.array([0, 0, 1]), np.array([False, True, False]),
+                              body_id=0, level=1, h=h, time_end=time_end, window=window,
+                              reference_center=(0., 0.), duration=duration)
+
+
+class ExchangeWindowReduction(unittest.TestCase):
+    def test_window_duration_scaling(self):
+        rows = np.array([[1., 2., 3., 4.], [.5, -1., 0., 1.], [9., 9., 9., 9.]])
+        out = reduce(rows, "window", duration=4)
+        np.testing.assert_array_equal(out["total"]["impulse"], [1.5, 1.])
+        np.testing.assert_array_equal(out["total"]["force"], [1.5 / 4, 1. / 4])
+        self.assertEqual(out["total"]["torque"], 3. / 4)
+        self.assertEqual(out["total"]["mass_rate"], 5. / 4)
+        np.testing.assert_array_equal(out["aperture"]["impulse"], [.5, -1.])
+        self.assertEqual((out["time_start"], out["time_end"]), (6., 10))
+        self.assertEqual(out["sample_count"], 8)  # duration / h fine substeps
+        self.assertEqual(out["window"], "window")
+
+    def test_existing_windows_unchanged(self):
+        rows = np.array([[1., 2., 3., 4.], [0., 0., 0., 0.], [0., 0., 0., 0.]])
+        base = reduce(rows, "last_base_step")
+        self.assertEqual((base["time_start"], base["sample_count"]), (9., 2))
+        np.testing.assert_array_equal(base["total"]["force"], [1., 2.])
+        sub = reduce(rows, "last_substep")
+        self.assertEqual(sub["sample_count"], 1)
+
+    def test_duration_validation(self):
+        rows = np.zeros((3, 4))
+        for bad in (None, 0, -1, 1.5, True):
+            with self.subTest(duration=bad), self.assertRaises(ValueError):
+                reduce(rows, "window", duration=bad)
+        with self.assertRaises(ValueError):
+            reduce(rows, "last_base_step", duration=2)
+
+
+if __name__ == "__main__": unittest.main()

tests/unit/test_refined_dfg.py

@@ -13,7 +13,11 @@ h, periodic = r.helpers(r.DEFAULT_SOURCE)
 
 
 def test_source_pin():
-    assert r.verify_source(r.DEFAULT_SOURCE)["commit"].startswith("cd3a657")
+    historical = json.loads((r.HERE / "data/refined_dfg/source_manifest.json").read_text())
+    assert historical["commit"].startswith("cd3a657")
+    assert r.verify_source(r.DEFAULT_SOURCE, "MRT")["commit"].startswith("3ea131e")
+    with unittest.TestCase().assertRaisesRegex(RuntimeError, "Source mismatch"):
+        r.verify_source(r.DEFAULT_SOURCE, "SRT")
 
 
 def test_matched_units_and_planes():
@@ -135,16 +139,17 @@ def test_axis_specific_sampling_phase():
 
 def test_full_requires_completed_pilot():
     with unittest.TestCase().assertRaisesRegex(ValueError, "pilot-results"):
-        r.verify_pilot(None, r.verify_source(r.DEFAULT_SOURCE))
+        r.verify_pilot(None, r.verify_source(r.DEFAULT_SOURCE, "MRT"))
 
 
 def pilot_fixture(path):
     """Synthetic contract fixture only; never a physical result or run receipt."""
-    manifest = r.verify_source(r.DEFAULT_SOURCE)
+    manifest = r.verify_source(r.DEFAULT_SOURCE, "MRT")
     (path / "plan.json").write_text(
         json.dumps(
             dict(
                 source_manifest=manifest,
+                collision="MRT",
                 bounds_root_edge=list(r.BOUNDS),
                 runner_sha256=r.runner_hash(),
             )
@@ -177,6 +182,7 @@ def pilot_fixture(path):
             json.dumps(
                 dict(
                     case=case,
+                    collision="MRT",
                     mode="pilot",
                     status="pilot_complete",
                     healthy=True,
@@ -192,6 +198,7 @@ def pilot_fixture(path):
             json.dumps(
                 dict(
                     case=case,
+                    collision="MRT",
                     mode="pilot",
                     step=interval,
                     runner_sha256=r.runner_hash(),
@@ -208,7 +215,7 @@ def test_pilot_validation_and_calibration():
     with TemporaryDirectory() as tmp:
         path = Path(tmp)
         manifest = pilot_fixture(path)
-        assert r.verify_pilot(path, manifest) == 3606.0
+        assert r.verify_pilot(path, manifest, "MRT") == 3606.0
 
 
 def test_reject_malformed_pilot_receipts():
@@ -238,7 +245,7 @@ def test_reject_malformed_pilot_receipts():
             mutate(bad)
             receipt.write_text(json.dumps(bad))
             with unittest.TestCase().assertRaises(ValueError):
-                r.verify_pilot(path, manifest)
+                r.verify_pilot(path, manifest, "MRT")
         # Exact reviewer regression: previously accepted NaN endpoint.
         receipt.write_text(
             json.dumps(
@@ -250,11 +257,11 @@ def test_reject_malformed_pilot_receipts():
             )
         )
         with unittest.TestCase().assertRaises(ValueError):
-            r.verify_pilot(path, manifest)
+            r.verify_pilot(path, manifest, "MRT")
         receipt.write_text(json.dumps(valid))
         (path / "uniform/ready.json").write_text("{}")
         with unittest.TestCase().assertRaises(ValueError):
-            r.verify_pilot(path, manifest)
+            r.verify_pilot(path, manifest, "MRT")
 
 
 def test_strict_finite_json_and_deadline():
@@ -283,6 +290,51 @@ def test_pressure_probes_actual_body_mask():
     assert h.pressure_probes(np.ones(pair.fine_shape), fluid, p)["dp"] == 0.0
 
 
+def test_matched_mrt_selection():
+    for diameter in (40, 80):
+        srt, sbody, sp = r.config(h, diameter)
+        mrt, mbody, mp = r.config(h, diameter, "MRT")
+        assert mrt["method"]["collision"] == "MRT"
+        mrt["method"]["collision"] = "SRT"
+        assert mrt == srt and mbody == sbody and mp == sp
+    srt_plan, mrt_plan = r.plan(h), r.plan(h, "MRT")
+    for key in ("root", "uniform", "bounds_root_edge", "finest_dx_m", "finest_dt_s"):
+        assert srt_plan[key] == mrt_plan[key]
+    assert "q2/6/8" in mrt_plan["outlet"] and "all nine" in srt_plan["outlet"]
+    cfg, pair, _, host = r.host_geometry(r.DEFAULT_SOURCE, h, "MRT")
+    assert cfg.collision == "MRT" and pair.fine_shape == (240, 568)
+    assert len(host.boundary_cells) > 0
+
+
+def test_reject_collision_receipts():
+    from tempfile import TemporaryDirectory
+
+    for rel in ("plan.json", "uniform/summary.json", "refined/summary.json",
+                "uniform/ready.json", "refined/ready.json"):
+        for invalid in ("SRT", None):
+            with TemporaryDirectory() as tmp:
+                path = Path(tmp)
+                manifest = pilot_fixture(path)
+                receipt = path / rel
+                record = json.loads(receipt.read_text())
+                if invalid is None:
+                    record.pop("collision")
+                else:
+                    record["collision"] = invalid
+                receipt.write_text(json.dumps(record))
+                with unittest.TestCase().assertRaises(ValueError):
+                    r.verify_pilot(path, manifest, "MRT")
+
+
+def test_record_actual_checkout_separately():
+    from source_provenance import checkout_provenance
+
+    actual = checkout_provenance(r.DEFAULT_SOURCE)
+    assert len(actual["checkout_revision"]) == 40
+    assert isinstance(actual["checkout_dirty"], bool)
+    assert r.verify_source(r.DEFAULT_SOURCE, "MRT")["production_origin_commit"].startswith("3ea131e")
+
+
 if __name__ == "__main__":
     suite = unittest.TestSuite(
         unittest.FunctionTestCase(fn)

tests/unit/test_refinement.py

@@ -142,7 +142,7 @@ class FixedPatchTests(unittest.TestCase):
                 FixedPatch2D.from_config(LBMConfig(nx=nx, ny=ny, viscosity=.1), bounds)
 
     def test_unsupported_backend_and_fine_relaxation_rejected(self):
-        for kwargs in ({"collision": "TRT"}, {"collision": "MRT"}, {"streaming": "esopull"},
+        for kwargs in ({"collision": "TRT"}, {"streaming": "esopull"},
                        {"les_enabled": True},
                        {"store_precision": "FP16S"}, {"lattice_model": "D3Q19"},
                        {"data_type": "FP64"}, {"omega_min": 1.0}):

tests/unit/test_refinement_mrt.py

@@ -0,0 +1,44 @@
+"""Host guards and source contract for MRT reuse in refined/blocked stepping (no CUDA)."""
+from pathlib import Path
+import unittest
+
+from CelerisLab.config import LBMConfig
+from CelerisLab.lbm.refinement import FixedPatch2D
+
+KERNELS = Path(__file__).resolve().parents[2] / "src/CelerisLab/lbm/kernels"
+BOUNDS = (8, 6, 16, 14)
+
+
+class RefinementMRTHost(unittest.TestCase):
+    def test_patch_accepts_implemented_models(self):
+        for collision in ("SRT", "MRT"):
+            with self.subTest(collision=collision):
+                FixedPatch2D.from_config(LBMConfig(nx=24, ny=20, viscosity=.1, collision=collision), BOUNDS)
+
+    def test_patch_rejects_unimplemented(self):
+        for options in (dict(collision="TRT"), dict(collision="MRT", les_enabled=True)):
+            with self.subTest(options=options), self.assertRaises(ValueError):
+                FixedPatch2D.from_config(LBMConfig(nx=24, ny=20, viscosity=.1, **options), BOUNDS)
+
+    def test_srt_statements_and_guards_preserved(self):
+        """MRT is an additive compile-time branch; SRT statements and other guards remain."""
+        blocks = (KERNELS / "blocks.cu").read_text()
+        channel = (KERNELS / "step/block_channel.cu").read_text()
+        refinement = (KERNELS / "refinement.cu").read_text()
+        for text in ("compute_feq_from_momentum(drho, jx, jy, ux, uy, feq);",
+                     "compute_feq(rho, ux, uy, feq);",
+                     "collide_srt_no_force(f, feq, d_params.omega);"):
+            self.assertIn(text, blocks)
+        self.assertIn("collide_srt(f, feq, forcing, omega);", channel)
+        self.assertIn("collide_srt_no_force(f, feq, d_params.omega);", refinement)
+        for source in (blocks, refinement):
+            self.assertIn("(COLLISION_MODEL != 0 && COLLISION_MODEL != 2)", source)
+            self.assertIn("USE_LES", source)
+            self.assertIn("STORE_PRECISION != 0", source)
+            self.assertIn('#if COLLISION_MODEL == 2\n#include "operators/collision_mrt.cuh"\n#endif', source)
+        self.assertIn("USE_DDF_SHIFTING ||", refinement)  # dense/mixed shifted stays rejected
+        for source in (blocks, channel, refinement):
+            self.assertEqual(source.count("collide_mrt_no_force("), 1)
+
+
+if __name__ == "__main__": unittest.main()

tests/unit/test_smagorinsky_normalization.py

@@ -1,99 +0,0 @@
-"""Direct-current-header regression; historical candidate manifests stay frozen."""
-
-import importlib.util
-from pathlib import Path
-import shutil
-import sys
-from tempfile import TemporaryDirectory
-import unittest
-
-ROOT = Path(__file__).resolve().parents[2]
-HELPERS = ROOT / "tests/validation/les"
-sys.path.insert(0, str(HELPERS))
-spec = importlib.util.spec_from_file_location(
-    "smagorinsky_oracle", HELPERS / "test_candidate.py"
-)
-oracle = importlib.util.module_from_spec(spec)
-spec.loader.exec_module(oracle)
-HEADER = Path("src/CelerisLab/lbm/kernels/operators/turbulence_smag.cuh")
-DESCRIPTORS = Path("src/CelerisLab/lbm/kernels/core/descriptors.cuh")
-
-
-class SmagorinskyNormalization(unittest.TestCase):
-    def test_current_header_matches_independent_tensor_oracle(self):
-        with TemporaryDirectory() as directory:
-            build_dir = Path(directory)
-            for dim in (2, 3):
-                for cs in (0.0, 0.16):
-                    # Compile directly against this checkout, not saved candidate headers.
-                    closed = oracle.build(
-                        ROOT.parent, ROOT.name, dim, cs, 1, 3, build_dir
-                    )
-                    converged = oracle.build(
-                        ROOT.parent, ROOT.name, dim, cs, 0, 200, build_dir
-                    )
-                    for rho in (0.7, 1.0, 1.3):
-                        for tau0 in (0.56, 0.8):
-                            for kind, magnitude in (
-                                ("equilibrium", 0.0),
-                                ("shear", -0.02),
-                                ("shear", 0.02),
-                                ("extension", -0.02),
-                                ("extension", 0.02),
-                            ):
-                                stress = oracle.tensor(
-                                    dim, kind, magnitude * rho, 0, dim - 1
-                                )
-                                f, eq = oracle.state(dim, rho, stress)
-                                expected = 1 / oracle.standard_tau(
-                                    tau0, cs, oracle.norm(stress), rho
-                                )
-                                with self.subTest(
-                                    dim=dim,
-                                    cs=cs,
-                                    rho=rho,
-                                    tau=tau0,
-                                    kind=kind,
-                                    sign=magnitude,
-                                ):
-                                    self.assertAlmostEqual(
-                                        closed(f, eq, rho, tau0), expected, delta=3e-6
-                                    )
-                                    self.assertAlmostEqual(
-                                        converged(f, eq, rho, tau0),
-                                        expected,
-                                        delta=3e-6,
-                                    )
-                    if cs == 0:
-                        f, eq = oracle.state(
-                            dim, 1.0, oracle.tensor(dim, "equilibrium", 0.0)
-                        )
-                        self.assertAlmostEqual(
-                            closed(f, eq, 1.0, 0.50001), 1.99, delta=3e-6
-                        )
-                        self.assertAlmostEqual(
-                            closed(f, eq, 1.0, 200.0), 0.01, delta=3e-6
-                        )
-
-    def test_removing_sqrt_two_is_detected(self):
-        with TemporaryDirectory() as directory:
-            root = Path(directory)
-            for relative in (HEADER, DESCRIPTORS):
-                target = root / "old_convention" / relative
-                target.parent.mkdir(parents=True, exist_ok=True)
-                shutil.copyfile(ROOT / relative, target)
-            target = root / "old_convention" / HEADER
-            text = target.read_text()
-            factor = "18.0f * sqrtf(2.0f) * smag_const_sq"
-            self.assertEqual(text.count(factor), 2)
-            target.write_text(text.replace(factor, "18.0f * smag_const_sq"))
-            for dim in (2, 3):
-                old = oracle.build(root, "old_convention", dim, 0.16, 1, 3, root)
-                stress = oracle.tensor(dim, "shear", 0.02)
-                f, eq = oracle.state(dim, 1.0, stress)
-                expected = 1 / oracle.standard_tau(0.8, 0.16, oracle.norm(stress), 1.0)
-                self.assertGreater(abs(old(f, eq, 1.0, 0.8) - expected), 1e-3)
-
-
-if __name__ == "__main__":
-    unittest.main()

tests/unit/test_time_refined_dfg.py

@@ -0,0 +1,86 @@
+"""CPU contract of the measurement-only refined DFG timing harness (no CUDA)."""
+from pathlib import Path
+import sys
+import unittest
+
+import numpy as np
+
+sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "validation"))
+import time_refined_dfg as t  # noqa: E402
+
+
+class TimingHarnessContract(unittest.TestCase):
+    def test_equal_simulated_duration(self):
+        s = t.schedule(2000, 200)
+        self.assertEqual(s["refined"], (200, 2000))
+        self.assertEqual(s["uniform"], (400, 4000))  # one refined base step = two D80 steps
+        self.assertAlmostEqual(s["tstar_window"], 2000 * .01 / 40)
+        for bad in ((0, 1), (1, -1)):
+            with self.assertRaises(ValueError):
+                t.schedule(*bad)
+
+    def test_timing_only_labels(self):
+        for name in ("B", "C", "E'"):
+            self.assertFalse(t.VARIANTS[name]["force_history_valid"])
+            self.assertTrue(t.VARIANTS[name]["timing_only"])
+        for name in ("A", "E"):
+            self.assertFalse(t.VARIANTS[name]["force_history_valid"])  # reads discarded
+            self.assertTrue(t.VARIANTS[name]["force_reads"])
+            self.assertIn("runner-like", t.VARIANTS[name]["force_note"])
+        self.assertEqual({v["case"] for v in t.VARIANTS.values()}, {"refined", "uniform"})
+        self.assertIn("final base interval", t.VARIANTS["C"]["force_note"])
+        self.assertIn("whole-window", t.VARIANTS["E'"]["force_note"])
+
+    def test_deadline(self):
+        import time
+        with self.assertRaises(TimeoutError):
+            t.check_deadline(time.monotonic() - 1)
+        t.check_deadline(time.monotonic() + 60)
+
+    def test_parity_detects_bits_and_nan_masks(self):
+        a = np.array([[1., np.nan], [2., 3.]], np.float32)
+        self.assertIsNone(t.same_state([a], [a.copy()]))
+        b = a.copy(); b[1, 1] = np.nextafter(np.float32(3), np.float32(4))
+        self.assertIn("bits", t.same_state([a], [b]))
+        c = a.copy(); c[0, 0] = np.nan
+        self.assertIn("NaN", t.same_state([a], [c]))
+        self.assertIn("level count", t.same_state([a], [a, a]))
+        inf = a.copy(); inf[1, 0] = np.inf
+        self.assertIn("infinite", t.same_state([inf], [inf.copy()]))  # matching inf still fails
+        self.assertIn("infinite", t.same_state([a], [inf]))
+
+    def test_flat_uniform_ddf(self):
+        fluid = np.array([[True, False, True], [True, True, False]])
+        flat = np.ones(9 * fluid.size, np.float32)  # Simulation.get_ddf() layout
+        self.assertTrue(t.uniform_populations_finite(flat, fluid))
+        flat[1 * fluid.size + 1] = np.nan  # q=1, solid cell 1: ignored
+        self.assertTrue(t.uniform_populations_finite(flat, fluid))
+        flat[4 * fluid.size + 0] = np.nan  # q=4, fluid cell 0: detected
+        self.assertFalse(t.uniform_populations_finite(flat, fluid))
+        self.assertTrue(t.uniform_populations_finite(np.ones((9,) + fluid.shape), fluid))
+        with self.assertRaises(RuntimeError):
+            t.uniform_populations_finite(np.ones(9 * fluid.size + 1), fluid)
+
+    def test_physical_health_limits(self):
+        ok = (np.array([1., 1.01]), np.array([.01, .02]), np.array([0., .01]))
+        health = t.physical_health([ok])
+        self.assertAlmostEqual(health["rho_max"], 1.01)
+        for bad in ((np.array([1., np.nan]),) + ok[1:],
+                    (np.array([.4, 1.]),) + ok[1:],
+                    (ok[0], np.array([.2, 0.]), ok[2]),
+                    (np.array([]), np.array([]), np.array([]))):
+            with self.assertRaises(RuntimeError):
+                t.physical_health([bad])
+
+    def test_summary_reports_differences_without_claims(self):
+        rows = [dict(variant=v, seconds=s) for v, s in
+                (("A", 3.), ("A", 5.), ("B", 2.), ("C", 1.), ("E", 2.), ("E'", 1.5))]
+        s = t.summarize(rows, .5)
+        self.assertEqual(s["A"]["median_seconds"], 4.)
+        self.assertEqual(s["A_minus_B_seconds"], 2.)
+        self.assertEqual(s["B_minus_C_seconds"], 1.)
+        self.assertEqual(s["A_over_E"], 2.)
+        self.assertIn("not exact causal", s["interpretation"])
+
+
+if __name__ == "__main__": unittest.main()

tests/unit/test_time_window_dfg.py

@@ -0,0 +1,38 @@
+"""CPU contract of the A-versus-window timing harness (no CUDA)."""
+from pathlib import Path
+import sys
+import unittest
+
+import numpy as np
+
+sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "validation"))
+import time_window_dfg as w  # noqa: E402
+
+
+class WindowTimingContract(unittest.TestCase):
+    def test_interval_and_variants(self):
+        self.assertEqual(w.INTERVAL, 80)
+        self.assertEqual(w.VARIANTS, ("A", "W", "W1"))
+
+    def test_components(self):
+        obs = dict(total=dict(impulse=np.array([1., 2.]), angular_impulse=3., mass_into_fluid=4.))
+        np.testing.assert_array_equal(w.components(obs), [1., 2., 3., 4.])
+
+    def test_interval_comparison(self):
+        ref = np.array([[1., -1., 0., 0.]])
+        scale = np.array([[10., 10., 1., 1.]])
+        self.assertEqual(w.compare_intervals(ref, scale, ref.copy()), 0.)
+        near = ref + np.array([[5e-9, 0., 0., 0.]])
+        self.assertLessEqual(w.compare_intervals(ref, scale, near), 1.)
+        far = ref + np.array([[1e-7, 0., 0., 0.]])
+        self.assertGreater(w.compare_intervals(ref, scale, far), 1.)
+        zero = np.zeros((1, 4))
+        self.assertGreater(w.compare_intervals(zero, zero, zero + 1e-12), 1.)  # floor, no div-by-0
+
+    def test_candidate_identity_records_files(self):
+        identity = w.candidate_identity(Path(__file__).resolve().parents[2])
+        self.assertIn("checkout_revision", identity)
+        self.assertEqual(set(identity["sha256"]), set(w.CANDIDATE_FILES))
+
+
+if __name__ == "__main__": unittest.main()

tests/validation/data/refined_dfg/source_manifest_mrt.json

@@ -0,0 +1,112 @@
+{
+  "commit": "3ea131e92c057843034f11a1c574301a2430dfce",
+  "production_origin_commit": "3ea131e92c057843034f11a1c574301a2430dfce",
+  "collision": "MRT",
+  "sha256": {
+    "src/CelerisLab/__init__.py": "21cc469bd6e970a1f9ff5970c8eea4602065e41fc881854759ee565ee55c10c6",
+    "src/CelerisLab/body/__init__.py": "e4e0e8f0769a94428d5b02fac790443997556cfc5515a5f1ed1595d2c694d69d",
+    "src/CelerisLab/body/action_smoother.py": "b1a23ee6bbd4850b172cb8738c5c3546ee465a54189e6af9d8044f023d4687dc",
+    "src/CelerisLab/body/coupling/__init__.py": "c931f7eb318d48408c3e2febcc6917c88c8dcc1697cbf4329b1f340d94c35763",
+    "src/CelerisLab/body/coupling/soa_packer.py": "8ec26fd55ef74801bbd5ab6bd3d0a7dabc750cbd498aaac2ca4c61a0108ad3c9",
+    "src/CelerisLab/body/ddf_patch.py": "60fcfb7441e73bf290bcff41ca3a83da3e7ddad765d7247464e8f02e3043b4e2",
+    "src/CelerisLab/body/geometry/__init__.py": "ada88980ab672959aca9aa3a4c4aade3219524b2ab68110cd2a2f41e3e0ca90a",
+    "src/CelerisLab/body/geometry/base.py": "cceca14361ca2f34d34a3e43730e1a79b8ca81be00f0dd62e884db32b479419f",
+    "src/CelerisLab/body/geometry/circle.py": "7605ca7c21103f1d801ac73ab8dd18dffbdac425b2dee7b8858b283628a9e436",
+    "src/CelerisLab/body/jets.py": "4104c24af4651f4bb4f8d8fd36fb6d88f41c6e4469822f36e22e74c9ede31dac",
+    "src/CelerisLab/body/manager.py": "01b653d25bd9b4333379c8bfd12383596b686d2029f2e91b1c91b618540ec6b3",
+    "src/CelerisLab/body/objects.py": "8e63cee8faa9cc6243bcff428be628dfad73a02a9a3b8ece7d78c0a63406a08d",
+    "src/CelerisLab/body/preprocess/__init__.py": "dfd9d58bc07f626ab23a70053e31b88a6d9165d2e4c76a4e9f48abe674d3d2ec",
+    "src/CelerisLab/body/preprocess/flag_overlay.py": "831fc9e9736b758847759a085ff133eb49ec91318ace4fda76f85960138ee300",
+    "src/CelerisLab/body/registry.py": "fbb7327780adbdc6da15ebdb799980daf32df22ee974af0893d7409f177bbb68",
+    "src/CelerisLab/body/sync_plan.py": "3a15bfe9ebeb91482dbfc68850761ef39c872cc12d2efc795babd0c01e392323",
+    "src/CelerisLab/common/__init__.py": "d4dbde78bc97fccc9eabac9abae8f0f15c112c10ff8e9b45bd9401d4723df1ef",
+    "src/CelerisLab/common/_types.py": "aad2ff1a9033d4fd68783ad7123afb47201ddc58785cb6dc0246947037072a5f",
+    "src/CelerisLab/common/checkpoint.py": "b2f960aecd5ac27a3a7130b809ffad4709c44d522edef0359ab952654625410b",
+    "src/CelerisLab/common/pathline.py": "a1b6efcc0a076d0f1dc8d49a595779426d7d5626b00072a59bac710ff02ad046",
+    "src/CelerisLab/common/preprocess.py": "7fc78a08857cf8e98519500b094a1d896af189b6b15658ff9b388ed51bed4515",
+    "src/CelerisLab/common/render.py": "0e055476257eb2bcd11e733766e52195de02f8fad228f570a3511a7c831bdc29",
+    "src/CelerisLab/common/streakline/__init__.py": "13ff1cfa416fdd721a08a659e76b5a5a99b8d3201add10a83b7b0bc654bc1eec",
+    "src/CelerisLab/common/streakline/_config.py": "a9185ef68e9c7b5d5a1dc0b04a28341e2ad91d13c145a7c106636a878cf0a605",
+    "src/CelerisLab/common/streakline/_integrate.py": "3fbd4368b2c8d0e7fda40dc65aac5df02b085008b4c68acd83cb7ae45d33a13b",
+    "src/CelerisLab/common/streakline/_render.py": "236cd637414ee318c4d66f9aed32110d8dc42a100d65fde06172f7cd0ba62c65",
+    "src/CelerisLab/common/streakline/_streakline.py": "1a3e08daa706b5f83cdeb337e9a945a36ddf7583c40af1d0bb9d825ac954df37",
+    "src/CelerisLab/config.py": "2f1bba7515f5bc08cb3a8d140bcafe7fa0689e01b3b06494697bb2956a71530e",
+    "src/CelerisLab/configs/CONFIG.md": "41cafdbe525366c08aa9d0416ea026ccf0026ba7123b8a9b65f968a80d928e3a",
+    "src/CelerisLab/configs/config_body.json": "055f261f0828ca1d4f2bb694bd15b72bbcabe31d552326208f28e37f761adab1",
+    "src/CelerisLab/configs/config_lbm.json": "17e4311f1f250c6ece59177a8b6ff856b69a016c7c10267e774b030a0a390865",
+    "src/CelerisLab/configs/config_lbm_flume_mrt_zouhe.json": "8f593321db05a563a42931603a228d61cbcba06dfbe113262d5cd0ad93c0fe79",
+    "src/CelerisLab/configs/config_lbm_flume_trt_regularized.json": "f4bba8e6a4339fc3c4eb14f6e2f57c4b22993c73703c773d853edfbe34736e56",
+    "src/CelerisLab/configs/config_lbm_three_cylinder_triangle.json": "eeb8e7149dbe4aa6b3a9562802735aa950750dfbe2d972a0727e788573afa817",
+    "src/CelerisLab/cuda/__init__.py": "0d01dff9f0399c8baba822aab5ca0122a135f011008ed8a499ebd9f7e5b06aea",
+    "src/CelerisLab/cuda/compiler_v2.py": "c4cc3ec4445c8d75e6a253dd87963877117a722bc73b9d4811806fc6045d08c1",
+    "src/CelerisLab/cuda/context.py": "013dd4ec6d01d2fbd11afdfb103eac6a273bd9f798d663d761b851eeb27885e9",
+    "src/CelerisLab/cuda/obs_layout.py": "5ffb2c3aeeaaaaeb64f2cd8febfbcf4ad9919e5414020ad753cfade0683b21c2",
+    "src/CelerisLab/lbm/__init__.py": "7f60fb3aaae81501b510d08c2fedff44fcb05657260eb10b6f5cc91e7204306b",
+    "src/CelerisLab/lbm/block_stepper.py": "00b313daddc6e63d8ab8423c9190ca7e92ad291de20a2f6b80eb53e3f2d85ec7",
+    "src/CelerisLab/lbm/blocks.py": "1e644c5719993f418185ad4a10b6950a30636fc912befd9eab41ff48b4cf85cd",
+    "src/CelerisLab/lbm/channel.py": "40d0354e5dc97bcdb7e7b5c1e7f9303a109b7134d8a3af88b8ff4cd75ab5a095",
+    "src/CelerisLab/lbm/curved_links.py": "b6957fa13d13f864dba4f7364e93a39905d062dbb3ee9b3f5191add81fff0732",
+    "src/CelerisLab/lbm/descriptors.py": "c9f025dfc14da43b989047b84cc22bd05d5cc92c9ec695c1da68fdfb788e3f2b",
+    "src/CelerisLab/lbm/equilibrium.py": "2fb2f2bb7ff326d23abce80fc2a4c0907b17740ea3e29fe6ae079c5044783826",
+    "src/CelerisLab/lbm/field.py": "0e727eacc3706add4d8d15c418e9b631050871f831698d3ac689dd5f2c701b9b",
+    "src/CelerisLab/lbm/hierarchy.py": "f234975074f6b6d289d573ff6ea6832feae0e27cb04f80b37fe81792b13a12e7",
+    "src/CelerisLab/lbm/initializers.py": "dacf9586963e862d6d0a98ef1a066036be6f91e9d423720661d8ef2c1fb2df6a",
+    "src/CelerisLab/lbm/kernels/blocks.cu": "2fa894b20a29f006911ac3e3aad0793cc4aaf9ccd7a67d01c6320a0836b1d8ef",
+    "src/CelerisLab/lbm/kernels/boundary/bounce_back.cuh": "aba3748c564ed91b161d268ef8b737b7c7be96a823063d201f3da684dd845233",
+    "src/CelerisLab/lbm/kernels/boundary/curved_boundary.cuh": "aeb16bfd483514fb03e15a2fed437d673c3b8b2fe68ab83edf4244ef593878ea",
+    "src/CelerisLab/lbm/kernels/boundary/esopull_curved_helpers.cuh": "b9b474bc4edaff05681c5443235472d04072d0d775fc0eb9c5967faef52b0258",
+    "src/CelerisLab/lbm/kernels/boundary/ibm_kernels.cuh": "e6909b547dc8d64cdd8a76a90660243a8f5fee535caa161b52b4dd8b7c29934b",
+    "src/CelerisLab/lbm/kernels/boundary/inlet/channel_stabilized.cuh": "efe077afcbb6b731256a987c2fd79a1dbc470c902f87070d60421ea6cb13d4e2",
+    "src/CelerisLab/lbm/kernels/boundary/inlet/common.cuh": "59b38a89439878e00de7de16b9c7b03feff443fcc5a421593f7d2536dcbb7ab8",
+    "src/CelerisLab/lbm/kernels/boundary/inlet/equilibrium.cuh": "f0964c811b2c41a5aaa1c6d99c490b06e850c634eb8749f1f05ef61f93c2de1d",
+    "src/CelerisLab/lbm/kernels/boundary/inlet/regularized.cuh": "943d4a571818bacb89c70bccfaac202691121e3e385803ce6309123a0a03ddaf",
+    "src/CelerisLab/lbm/kernels/boundary/inlet/zou_he_local.cuh": "e73e66b853b1e4e0e30d8f4185d832aa414ce670a1813ce2b5459f2887f71f89",
+    "src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh": "de8b213a6da42c2e5e9f8966199fc66ea4541b85144a90f3f5dbad3c6ab63b97",
+    "src/CelerisLab/lbm/kernels/boundary/outlet/pressure_neq.cuh": "c8b2adc66d12fe602076ddc9693c5c208711fcaddb58aa9a31401a269f51c4c4",
+    "src/CelerisLab/lbm/kernels/boundary/planar_velocity.cuh": "501b2759d0e2e4280cb6fdd127d66a50fa89f6de48fc2cae9a7eafb11ae90bed",
+    "src/CelerisLab/lbm/kernels/config/config_grid.h": "1db3e030ebdc116cd7f02362d1fd077daac55e8f0b4329f6285948c224d3ca5d",
+    "src/CelerisLab/lbm/kernels/config/config_method.h": "4374e04a130f90f97d17cc20c40d2f9ea0ed0dd8603817d07f215d02266e247c",
+    "src/CelerisLab/lbm/kernels/config/config_objects.h": "db793af052c756561a190a3ff2cb93ee5fe28462ee4f97c6af80882b4dc8b072",
+    "src/CelerisLab/lbm/kernels/config/config_obs.h": "c11e2249359145b5c4cbaaf0487506fa77c1f4b82c7c169ddc574aa0babcd727",
+    "src/CelerisLab/lbm/kernels/config/config_physics.h": "8ac0990b168e1058588676cc05c80df8634cfd671dbc92a89246c3c372d83a4e",
+    "src/CelerisLab/lbm/kernels/config.h": "7d1124125e0b354558297934460e64c760fb68cdd7bf6d01a1238bb2eaab8c37",
+    "src/CelerisLab/lbm/kernels/core/descriptors.cuh": "8440bc14857da9007982650f828a23608eec9f7cc4bf772e8c8dc1239cb45356",
+    "src/CelerisLab/lbm/kernels/core/flags.cuh": "0025705f24710330d9d9195f713879bfb650a85a3449baef0e196b831af104ab",
+    "src/CelerisLab/lbm/kernels/core/layout.cuh": "1fbd820021c89ac3efcde24145996cb0147769e607e498e61f308ed3e0b13706",
+    "src/CelerisLab/lbm/kernels/core/obs.cuh": "b742b77c9f61e1f020fe9922c3838ae9d855e17b381154958b58f10a22938922",
+    "src/CelerisLab/lbm/kernels/core/params.cuh": "6e18bf730d936d940e3d85452f0714a2eaae31bf41002f9fd520ec5b372878dd",
+    "src/CelerisLab/lbm/kernels/core/precision.cuh": "849df14b4cce50af52399e4b8f3e997196f3ae5f1163f68e28df3e6b2507a8c4",
+    "src/CelerisLab/lbm/kernels/kernel_v2.cu": "c7bd3fb705b77261ba2441c52a24a7e6a7865c165cba088a1b75201c4d29c408",
+    "src/CelerisLab/lbm/kernels/operators/collision_mrt.cuh": "aee444eb60c47a3fd518dca975d1af05425a0818a0ec6c3fb3f6150f60c2cb44",
+    "src/CelerisLab/lbm/kernels/operators/collision_srt.cuh": "6b7ff842f392c8f853d2d1a97cc69bfe50f21283f73e659404da07cf6c95f5e1",
+    "src/CelerisLab/lbm/kernels/operators/collision_trt.cuh": "9725b9840b9ad8b4ec7be61ff5d06d3151e86d1a87585c13f2fbf934d0d2c5f6",
+    "src/CelerisLab/lbm/kernels/operators/forcing_guo.cuh": "969a430c1af03f32e3bd0df95d805e64f56334717c44dd4ba208ba503062a24b",
+    "src/CelerisLab/lbm/kernels/operators/helpers.cuh": "0b3a4d95d7d19f5fec09e5423b8f4cc99f8a226b541967425967330da5a051ec",
+    "src/CelerisLab/lbm/kernels/operators/macro.cuh": "2be20eb2b7d7f5522daaad7f4b57b1ff8cdf883efe4a67846a96783b2fe846a3",
+    "src/CelerisLab/lbm/kernels/operators/turbulence_smag.cuh": "ff7f000297ee0817e255ad9fcb6a13b647e3ddd5b206518d35f38d2895a9f501",
+    "src/CelerisLab/lbm/kernels/refinement.cu": "b45aad9c022cecb64f7e968ec3cae2d8385a3cc430ee409700c404a3f4dd8cf8",
+    "src/CelerisLab/lbm/kernels/step/aux_kernels.cu": "578a5b20a157241ee678547059abb1372dd5a589ad06bb2c15b987eb01ca1013",
+    "src/CelerisLab/lbm/kernels/step/block_channel.cu": "ee8a31398b8e5233d977c816e6b8e060c1eb2ec71320e594226dbca55593c0d0",
+    "src/CelerisLab/lbm/kernels/step/curved_double.cu": "bfa7a1df1e939020314d49755045f7a85fe4a19cc1c68dbe62653c66c5e7c2b3",
+    "src/CelerisLab/lbm/kernels/step/esopull_macro.cu": "c73116020ca217ede1e4de27c8269c31f9d91e33433ac8252a21e9b14187c50f",
+    "src/CelerisLab/lbm/kernels/step/esopull_sync.cu": "c0ff3e72900c3ab89a05dd03684e96d3b740cf1672d565727fc2af165933d087",
+    "src/CelerisLab/lbm/kernels/step/init_flow.cu": "305be3bb12fc0d039acec1701baabcf0e8a35e683c49c34f2af8ca9867d2b29c",
+    "src/CelerisLab/lbm/kernels/step/one_step_double.cu": "f9f8242b38e8e0cf8425b75e2b0215bc15363c62eaa0fb43c761b8e0931bfdc9",
+    "src/CelerisLab/lbm/kernels/step/one_step_esopull.cu": "aa223b2744ae6c9125673df311a96ce588032c1ba58bebb11fd33ecb990f8ed4",
+    "src/CelerisLab/lbm/kernels/streaming/esopull_semantic_helpers.cuh": "dd5a9091628b079efd97e53fd49e52d09dcee88529235db77a5648e02937af6e",
+    "src/CelerisLab/lbm/kernels/streaming/esopull_single_buffer.cuh": "77ed4ab9a0ebc5ff62a91a7c39066ef395391ae36e9040a69fb2bd8b1626dc79",
+    "src/CelerisLab/lbm/kernels/streaming/pull_double_buffer.cuh": "b65d741ed78b7b62f99531ba0c0ce1936663d8dc9640e7d7d16c707f1baf6549",
+    "src/CelerisLab/lbm/population_storage.py": "a6a86c29c52c892edfe4802fb773c387212e3f4cc86708d80e154986ddf88044",
+    "src/CelerisLab/lbm/raw_exchange.py": "8694202466eba848563e298656353b0f5db94e21a01fa05f51aad6dc2b9fd311",
+    "src/CelerisLab/lbm/refinement.py": "100d5e46ac894f6008b61364cdad5d37c2ca7ed949d8d71405572d6366526417",
+    "src/CelerisLab/lbm/refinement_stepper.py": "f2bb34fa90805ad21c3e7cf9f9a67d88d3ae3e249b4517d3898567b4bf9e364d",
+    "src/CelerisLab/lbm/refinement_tables.py": "fdd344024b89df3a83174dc810ca5e32db5e19eb8b1609f8f73b0821f56546b6",
+    "src/CelerisLab/lbm/stepper.py": "0dcb7cd7ae87e33e3dc67a38b6aa43042fbe673e37d9ad165342ee1d695933ce",
+    "src/CelerisLab/lbm/synchronized_remap.py": "7400fc4beaffd4649d6a9c524d38ce5fd6fa7eeb32eb5332e4030a4b6648f197",
+    "src/CelerisLab/simulation.py": "0a03904dc6de88667d7af699005a289adf8b09f0976cb0d099b041f4830a8f4c",
+    "tests/validation/re100_helpers.py": "289c6e01075ff361573f862b90cfd61e5980c9a8efcbc78aef69979b4004fd60",
+    "tests/validation/periodic.py": "a2212c0e2705126327a5019e9cafe64eeca398da1b75061ff23630c9843697a6",
+    "tests/validation/data/dfg_re100/reference_cycle.json": "073d1bda1629697efebad95667bfe80598e217d708dbffd7d5f6c6999379c55f",
+    "tests/validation/source_provenance.py": "6dfe8f2d819dba0e240ba761d31162a5033685d492d09db75eb2b8aaab3a5a4c"
+  }
+}

tests/validation/multigrid/test_exchange_window.py

@@ -0,0 +1,260 @@
+"""Optional observation-window exchange ledger: host contract (mocked driver) and GPU gate.
+
+Run CPU and GPU classes in separate processes; GPU callers select the device.
+ExchangeWindowGPU uses the SRT/unshifted Jet fixture; ExchangeWindowShiftedMRTGPU covers a
+shifted MRT, all-blocked, stationary-body case without Jets (no shifted-MRT Jet claim).
+GPU tolerance (conservative engineering tolerance, not a universal proof): the window ledger
+sums h*dt*raw per link over fine substeps, while the reference sums per-base interval ledgers
+over links and then over base steps. Both are FP64 reorderings of the same FP32-exact substep
+contributions, so the scale is the sum of the ACTUAL absolute substep contributions (read from
+last_raw after every fine hop), not of per-base totals that could hide substep cancellation:
+|window - reference| <= 4 * m * eps * sum|substep terms| + 1e-300, m = links * substeps.
+"""
+from contextlib import ExitStack
+import unittest
+from unittest.mock import patch
+
+import numpy as np
+import pycuda.driver as cuda
+
+from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+from test_refinement_blocked_base import mocked_backend
+from test_refinement_blocked_fine import BOUNDS, circle, config
+
+EPS = np.finfo(np.float64).eps
+
+
+def build(m, window=True, raw=True):
+    return FixedPatchStepper(config("periodic"), BOUNDS, context=m["context"], base_block=(3, 4),
+                             fine_block=(16, 20), fine_objects=[circle()],
+                             observe_raw_exchange=raw, observe_window_exchange=window)
+
+
+def same_bits(test, a, b):
+    """Identical NaN masks, no infinities, identical finite FP32 bits (incl. signed zero)."""
+    a, b = np.asarray(a), np.asarray(b)
+    test.assertEqual(a.dtype, b.dtype)
+    test.assertFalse(np.isinf(a).any() or np.isinf(b).any())
+    np.testing.assert_array_equal(np.isnan(a), np.isnan(b))
+    finite = np.isfinite(a)
+    np.testing.assert_array_equal(a[finite].view(np.uint32), b[finite].view(np.uint32))
+
+
+def accumulations(m):
+    return [c for c in m["calls"] if c[0] == "BlockAccumulateExchange"]
+
+
+class ExchangeWindowCPU(unittest.TestCase):
+    def test_constructor_validation(self):
+        with mocked_backend() as m:
+            with self.assertRaisesRegex(ValueError, "requires observe_raw_exchange"):
+                build(m, window=True, raw=False)
+            with self.assertRaisesRegex(ValueError, "must be bool"):
+                build(m, window=1)
+
+    def test_state_machine_and_launches(self):
+        with mocked_backend() as m, ExitStack() as stack:
+            stack.enter_context(patch.object(cuda, "memset_d8"))
+            stack.enter_context(patch.object(cuda, "memset_d8_async"))
+            stack.enter_context(patch.object(cuda, "memcpy_dtoh",
+                                             side_effect=lambda out, src: out.fill(0)))
+            backend = build(m)
+            b = backend._level_boundaries[1]
+            self.assertIsNotNone(b.window_integral)
+            with self.assertRaises(RuntimeError):
+                backend.start_exchange_window()  # not initialized
+            arrays = [np.full(level.shape, 1 / 9, np.float32) for level in backend._levels]
+            backend.initialize(*arrays)
+            with self.assertRaisesRegex(RuntimeError, "no valid exchange window"):
+                backend.read_exchange_window(0)
+            backend.step(1)
+            self.assertEqual(len(accumulations(m)), 2)  # closed window: interval ledger only
+            backend.start_exchange_window()
+            with self.assertRaisesRegex(RuntimeError, "at least one completed"):
+                backend.read_exchange_window(0)
+            m["calls"].clear()
+            backend.step(2)
+            calls = accumulations(m)
+            self.assertEqual(len(calls), 8)  # two fine hops x two ledgers x two base steps
+            self.assertEqual({id(c[1][1]) for c in calls},
+                             {id(b.interval_integral), id(b.window_integral)})
+            out = backend.read_exchange_window(0)
+            self.assertEqual((out["time_start"], out["time_end"], out["sample_count"]), (1., 3, 4))
+            self.assertEqual(backend.read_exchange_window(0)["time_end"], 3)  # non-destructive
+            backend.initialize(*arrays)  # invalidates
+            with self.assertRaisesRegex(RuntimeError, "no valid exchange window"):
+                backend.read_exchange_window(0)
+            backend.start_exchange_window()
+            with patch.object(backend, "_hop", side_effect=RuntimeError("injected")):
+                with self.assertRaises(RuntimeError):
+                    backend.step(1)
+            self.assertFalse(backend._window_valid)
+            with self.assertRaises(RuntimeError):
+                backend.read_exchange_window(0)  # failed status
+
+    def test_disabled_default_has_no_window_cost(self):
+        with mocked_backend() as m, ExitStack() as stack:
+            stack.enter_context(patch.object(cuda, "memset_d8_async"))
+            backend = build(m, window=False)
+            self.assertIsNone(backend._level_boundaries[1].window_integral)
+            backend.initialize(*[np.full(l.shape, 1 / 9, np.float32) for l in backend._levels])
+            backend.step(1)
+            self.assertEqual(len(accumulations(m)), 2)
+            with self.assertRaisesRegex(RuntimeError, "disabled"):
+                backend.start_exchange_window()
+
+
+class ExchangeWindowGPU(unittest.TestCase):
+    def test_window_matches_summed_base_reads_and_populations_unchanged(self):
+        from test_raw_exchange_gpu import BOUNDS as RB, CFG, circle as rc, jet, state
+        from CelerisLab.cuda.context import CudaContext
+        with CudaContext(0) as ctx:
+            kw = dict(context=ctx, fine_block=(16, 20), fine_objects=[rc()], jets=[jet()],
+                      observe_raw_exchange=True)
+            with FixedPatchStepper(CFG, RB, observe_window_exchange=True, **kw) as windowed, \
+                 FixedPatchStepper(CFG, RB, **kw) as plain:
+                arrays = [state(level.shape) for level in windowed._levels]
+                windowed.initialize(*arrays)
+                plain.initialize(*arrays)
+                windowed.step(2)
+                plain.step(2)
+                windowed.start_exchange_window()
+                b = windowed._level_boundaries[1]
+                count = len(b.tables["cell"])
+                factors = np.array([.25, .25, .125, .25])  # h*dt, h*h*dt, h*h for level-1 h=.5
+                scale = np.zeros(4)
+                substeps = 0
+                kernel = windowed._levels[1].kernel
+
+                def audited(src, dst, *args, **kw):
+                    nonlocal scale, substeps
+                    windowed._stream.synchronize()  # BFL and both accumulations are complete
+                    raw = np.empty((count, 4), np.float32)
+                    cuda.memcpy_dtoh(raw, b.last_raw)
+                    scale = scale + (np.abs(raw.astype(np.float64)) * factors).sum(axis=0)
+                    substeps += 1
+                    kernel(src, dst, *args, **kw)
+                windowed._levels[1].kernel = audited
+                keys = ("impulse", "angular_impulse", "mass_into_fluid")
+                summed = {part: [0., 0., 0., 0.] for part in ("total", "wall", "aperture")}
+                steps = 0
+                for speed in (0., .002, .002, -.001, 0., 0.):  # amplitude updates stay in-window
+                    windowed.set_jet_amplitude(0, speed)
+                    plain.set_jet_amplitude(0, speed)
+                    windowed.step(1)
+                    plain.step(1)
+                    steps += 1
+                    obs = windowed.read_raw_exchange(0)
+                    for part in summed:
+                        o = obs[part]
+                        summed[part] = [summed[part][0] + o["impulse"][0], summed[part][1] + o["impulse"][1],
+                                        summed[part][2] + o["angular_impulse"],
+                                        summed[part][3] + o["mass_into_fluid"]]
+                    for a, e in zip(windowed._download_raw().values, plain._download_raw().values):
+                        same_bits(self, a, e)
+                out = windowed.read_exchange_window(0)
+                self.assertEqual((out["time_end"] - out["time_start"], out["sample_count"]),
+                                 (steps, 2 * steps))
+                self.assertEqual(substeps, 2 * steps)
+                bound = 4 * (count * substeps) * EPS * scale + 1e-300
+                for part in summed:
+                    o = out[part]
+                    got = np.array([o["impulse"][0], o["impulse"][1], o["angular_impulse"],
+                                    o["mass_into_fluid"]])
+                    diff = np.abs(got - np.array(summed[part]))
+                    self.assertTrue((diff <= bound).all(), f"{part}: {diff} > {bound}")
+                    np.testing.assert_array_equal(o["force"], o["impulse"] / steps)
+                again = windowed.read_exchange_window(0)
+                np.testing.assert_array_equal(again["total"]["impulse"], out["total"]["impulse"])
+                windowed.initialize(*arrays)
+                with self.assertRaises(RuntimeError):
+                    windowed.read_exchange_window(0)
+
+
+def audit_substeps(stepper, b, count, record):
+    """Wrap the level-1 step kernel: after BFL and accumulation, add |substep terms| to record."""
+    factors = np.array([.25, .25, .125, .25])
+    kernel = stepper._levels[1].kernel
+
+    def audited(src, dst, *args, **kw):
+        stepper._stream.synchronize()
+        raw = np.empty((count, 4), np.float32)
+        cuda.memcpy_dtoh(raw, b.last_raw)
+        record["scale"] = record["scale"] + (np.abs(raw.astype(np.float64)) * factors).sum(axis=0)
+        record["substeps"] += 1
+        kernel(src, dst, *args, **kw)
+    stepper._levels[1].kernel = audited
+
+
+def components(obs):
+    o = obs["total"]
+    return np.array([o["impulse"][0], o["impulse"][1], o["angular_impulse"], o["mass_into_fluid"]])
+
+
+class ExchangeWindowShiftedMRTGPU(unittest.TestCase):
+    """Shifted MRT, all-blocked, stationary body, NO Jet (scoped separately from SRT/Jet)."""
+
+    def test_window_step0_reset_and_split_step(self):
+        from test_raw_exchange_gpu import BOUNDS as RB, circle as rc, state
+        from CelerisLab.config import LBMConfig
+        from CelerisLab.cuda.context import CudaContext
+        cfg = LBMConfig(nx=32, ny=24, viscosity=.1, collision="MRT", ddf_shifting=True)
+        with CudaContext(0) as ctx:
+            kw = dict(context=ctx, base_block=(12, 8), fine_block=(16, 20), fine_objects=[rc()],
+                      observe_raw_exchange=True)
+            with FixedPatchStepper(cfg, RB, observe_window_exchange=True, **kw) as windowed, \
+                 FixedPatchStepper(cfg, RB, observe_window_exchange=True, **kw) as split, \
+                 FixedPatchStepper(cfg, RB, **kw) as plain:
+                arrays = [state(level.shape) for level in windowed._levels]
+                for run in (windowed, split, plain):
+                    run.initialize(*arrays)
+                    run.step(2)
+                    if run is not plain:
+                        run.start_exchange_window()
+                b = windowed._level_boundaries[1]
+                count = len(b.tables["cell"])
+                record = dict(scale=np.zeros(4), substeps=0)
+                audit_substeps(windowed, b, count, record)
+                steps, summed = 5, np.zeros(4)
+                for _ in range(steps):
+                    windowed.step(1)
+                    plain.step(1)
+                    summed += components(windowed.read_raw_exchange(0))
+                    for a, e in zip(windowed._download_raw().values, plain._download_raw().values):
+                        same_bits(self, a, e)
+                split.step(steps)  # one call: same ledger arithmetic as five step(1) calls
+                for a, e in zip(windowed._download_raw().values, split._download_raw().values):
+                    same_bits(self, a, e)
+                rows = [np.empty((count, 4), np.float64) for _ in range(2)]
+                cuda.memcpy_dtoh(rows[0], b.window_integral)
+                cuda.memcpy_dtoh(rows[1], split._level_boundaries[1].window_integral)
+                np.testing.assert_array_equal(rows[0].view(np.uint64), rows[1].view(np.uint64))
+
+                out = windowed.read_exchange_window(0)
+                self.assertEqual(record["substeps"], 2 * steps)
+                self.assertEqual((out["time_end"] - out["time_start"], out["sample_count"]),
+                                 (steps, 2 * steps))
+                bound = 4 * (count * record["substeps"]) * EPS * record["scale"] + 1e-300
+                diff = np.abs(components(out) - summed)
+                self.assertTrue((diff <= bound).all(), f"{diff} > {bound}")
+                np.testing.assert_array_equal(split.read_exchange_window(0)["total"]["impulse"],
+                                              out["total"]["impulse"])
+
+                windowed.step(0)  # no interval: window and its duration unchanged
+                again = windowed.read_exchange_window(0)
+                self.assertEqual(again["time_end"], out["time_end"])
+                np.testing.assert_array_equal(again["total"]["impulse"], out["total"]["impulse"])
+
+                windowed.start_exchange_window()  # reset after a window
+                with self.assertRaisesRegex(RuntimeError, "at least one completed"):
+                    windowed.read_exchange_window(0)
+                record.update(scale=np.zeros(4), substeps=0)
+                windowed.step(1)
+                single = components(windowed.read_raw_exchange(0))
+                reset = windowed.read_exchange_window(0)
+                self.assertEqual(reset["time_end"] - reset["time_start"], 1)
+                bound = 4 * (count * record["substeps"]) * EPS * record["scale"] + 1e-300
+                self.assertTrue((np.abs(components(reset) - single) <= bound).all())
+
+
+if __name__ == "__main__": unittest.main()

tests/validation/multigrid/test_refinement_mrt_gpu.py

@@ -0,0 +1,95 @@
+"""MRT reuse GPU gates. Caller owns GPU admission (CUDA_VISIBLE_DEVICES).
+
+Run: python -m unittest discover -s tests/validation/multigrid -p test_refinement_mrt_gpu.py
+Parity compares the blocked MRT channel root with uniform MRT OneStep: the existing
+pressure-NEQ outlet branch differs by collision model, so SRT is never the reference.
+AMR checks are finite-evolution/split-invariance smoke tests, not benchmark validation.
+"""
+from dataclasses import replace
+import unittest
+import numpy as np
+
+from CelerisLab.config import LBMConfig
+from CelerisLab.cuda.context import CudaContext
+from CelerisLab.lbm.population_storage import _encode_physical
+from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+from test_channel_root_gpu import run_root, state
+from test_refinement_blocked_base import raw, exact
+
+
+class RefinementMRTGPU(unittest.TestCase):
+    def test_uniform_matched_mrt_root(self):
+        with CudaContext(0) as context:
+            for shifted in (False, True):
+                for profile in ("uniform", "parabolic"):
+                    cfg = LBMConfig(nx=24, ny=20, viscosity=.1, velocity=.025, collision="MRT",
+                                    ddf_shifting=shifted, inlet_profile=profile)
+                    initial = _encode_physical(state(cfg.ny, cfg.nx), shifted)
+                    expected = run_root(cfg, context, initial)
+                    self.assertTrue(np.isfinite(expected).all())
+                    for block in ((20, 24), (5, 6), (4, 1)):
+                        np.testing.assert_array_equal(run_root(cfg, context, initial, block), expected)
+
+    def test_amr_mrt_smoke(self):
+        from CelerisLab.body.geometry.circle import CircleGeometry
+        from CelerisLab.body.objects import SimObject
+        center = (24., 16.)
+        body = SimObject(0, geometry=CircleGeometry(*center, 2.15), center=center, radius=2.15)
+        with CudaContext(0) as context:
+            for boundary in ("periodic", "channel"):
+                for shifted in (False, True):
+                    cfg = LBMConfig(nx=48, ny=32, viscosity=.1, velocity=.015 if boundary == "channel" else 0.,
+                                    collision="MRT", ddf_shifting=shifted)
+                    reference = None
+                    for block in ((32, 48), (8, 8), (4, 1)):
+                        with FixedPatchStepper(cfg, (8, 6, 40, 26), context=context, coarse_boundary=boundary,
+                                               base_block=block, fine_block=(8, 8), fine_objects=(body,)) as backend:
+                            backend.initialize(*[state(*level.shape[1:]) for level in backend._levels])
+                            backend.step(8)
+                            raw = backend._download_raw()
+                            for a, mask in zip(raw.values, backend._state):
+                                self.assertTrue(np.isfinite(a[:, mask]).all())
+                            if reference is None: reference = raw.values
+                            else:
+                                for a, b in zip(reference, raw.values): np.testing.assert_array_equal(a, b)
+
+    def test_dense_versus_blocked_mrt_two_level(self):
+        """Unshifted MRT: all-dense PatchCoarseStep/PatchFineStep vs mixed/all-blocked storage."""
+        from test_refinement_blocked_base import BASE_BLOCKS, FINE_BLOCK
+        from test_refinement_blocked_fine import NX, NY, BOUNDS
+        import test_refinement_blocked_fine as fixture
+        cfg = LBMConfig(nx=NX, ny=NY, viscosity=.1, collision="MRT")
+        with CudaContext(0) as context:
+            results = {}
+            layouts = ((None, None), *((base, FINE_BLOCK) for base in (None, *BASE_BLOCKS)))
+            for base_block, fine_block in layouts:
+                with FixedPatchStepper(cfg, BOUNDS, context=context, fine_block=fine_block,
+                                       base_block=base_block) as backend:
+                    backend.initialize(*[fixture.state(*level.shape[1:], seed) for seed, level in enumerate(backend._levels)])
+                    backend.step(4)
+                    results[base_block, fine_block] = raw(backend)
+            for layout in layouts[1:]:
+                for level, (a, b) in enumerate(zip(results[layout], results[None, None])):
+                    exact(self, a, b, f"layout={layout} level={level}")
+
+    def test_shifted_all_blocked_three_level_mrt_smoke(self):
+        """Per-level shear rate omega_l = 1/(3*2**l*nu + 1/2); finite ACTIVE leaves."""
+        cfg = LBMConfig(nx=16, ny=16, viscosity=.1, collision="MRT", ddf_shifting=True)
+        kw = dict(fine_block=(4, 4), base_block=(4, 4), child=(8, 8, 16, 16), child_block=(4, 4))
+        with CudaContext(0) as context:
+            with FixedPatchStepper(cfg, (4, 4, 12, 12), context=context, **kw) as backend:
+                self.assertEqual(len(backend._level_cfgs), 3)
+                for level, level_cfg in enumerate(backend._level_cfgs):
+                    checked = replace(level_cfg)
+                    checked.validate()  # omega is derived only by validate()
+                    self.assertEqual(checked.collision, "MRT")
+                    self.assertTrue(checked.ddf_shifting)
+                    self.assertTrue(np.isclose(checked.omega, 1 / (3 * 2**level * .1 + .5), rtol=1e-14))
+                backend.initialize(*[state(*level.shape[1:]) for level in backend._levels])
+                backend.step(4)
+                values = backend._download_raw().values
+                for a, mask in zip(values, backend._state):
+                    self.assertTrue(np.isfinite(a[:, mask]).all())
+
+
+if __name__ == "__main__": unittest.main()

tests/validation/run_refined_dfg.py

@@ -1,4 +1,4 @@
-"""Frozen matched SRT DFG comparison. GPU execution requires explicit root admission."""
+"""Frozen matched SRT/MRT DFG comparison. GPU execution requires explicit root admission."""
 
 import argparse
 import csv
@@ -39,20 +39,21 @@ def helpers(source):
     return h, periodic
 
 
-def config(h, diameter):
-    cfg, body, p = h.make_config(diameter, SPEED, "SRT", "neq_extrap", "zou_he_local")
+def config(h, diameter, collision="SRT"):
+    cfg, body, p = h.make_config(diameter, SPEED, collision, "neq_extrap", "zou_he_local")
     cfg["method"]["outlet"]["srt_neq_damp"] = 0.5
     return cfg, body, p
 
 
-def plan(h):
-    _, _, c = config(h, 40)
-    _, _, f = config(h, 80)
+def plan(h, collision="SRT"):
+    _, _, c = config(h, 40, collision)
+    _, _, f = config(h, 80, collision)
     x0, y0, x1, y1 = BOUNDS
     fine_active = 4 * (x1 - x0) * (y1 - y0)
     coarse_uncovered = c["nx"] * c["ny"] - fine_active // 4
     return dict(
         status="prepared_not_run",
+        collision=collision,
         root_gpu_gates_required=True,
         root=c,
         uniform=f,
@@ -82,7 +83,11 @@ def plan(h):
         / (2 * f["nx"] * f["ny"]),
         cost_caveat="Geometry-only estimates include ghosts/solids; exclude block halos, transfers, observation and launch overhead. Not measured speedup.",
         sampling_offset="For origin (38,22), fine physical x=(j+75.5)*dx80 versus uniform i*dx80: half-fine-cell x offset. Fine physical y=(k+42.5)*dx80 equals uniform (l-.5)*dx80 at l=k+43: y samples align. Body, wall and source planes coincide.",
-        outlet="SRT pressure NEQ, all nine fneq damped by .5, no backflow clamp; noncanonical DFG outlet",
+        outlet=(
+            "SRT pressure NEQ, all nine fneq damped by .5, no backflow clamp; noncanonical DFG outlet"
+            if collision == "SRT" else
+            "MRT pressure NEQ, unknowns q2/6/8 extrapolated; no backflow clamp; noncanonical DFG outlet"
+        ),
         modes={"pilot_tstar": 5.0, "full_tstar": 300.0},
         observation_dense_tstar=0.02,
     )
@@ -120,7 +125,7 @@ def initial_state(shape, origin, spacing, height=164.0, speed=SPEED, root=False)
     return w * (1 + 3 * cu + 4.5 * cu * cu - 1.5 * velocity * velocity)
 
 
-def host_geometry(source, h):
+def host_geometry(source, h, collision="SRT"):
     """Production donor/transfer and complete cut-link checks, without CUDA."""
     from CelerisLab.config import load_lbm_config
     from CelerisLab.lbm.refinement import FixedPatch2D, prepare_fine_boundary
@@ -133,7 +138,7 @@ def host_geometry(source, h):
     from CelerisLab.body.geometry.circle import CircleGeometry
     from CelerisLab.body.objects import SimObject
 
-    cfgdict, _, _ = config(h, 40)
+    cfgdict, _, _ = config(h, 40, collision)
     with TemporaryDirectory() as tmp:
         p = Path(tmp) / "cfg.json"
         p.write_text(json.dumps(cfgdict))
@@ -151,8 +156,11 @@ def host_geometry(source, h):
     return cfg, pair, body, host
 
 
-def verify_source(source):
-    expected = json.loads((HERE / "data/refined_dfg/source_manifest.json").read_text())
+def verify_source(source, collision="SRT"):
+    if collision not in ("SRT", "MRT"):
+        raise ValueError("Matched collision must be SRT or MRT")
+    name = "source_manifest.json" if collision == "SRT" else "source_manifest_mrt.json"
+    expected = json.loads((HERE / "data/refined_dfg" / name).read_text())
     for name, digest in expected["sha256"].items():
         if hashlib.sha256((source / name).read_bytes()).hexdigest() != digest:
             raise RuntimeError("Source mismatch: " + name)
@@ -167,7 +175,7 @@ def run_case(a, h, periodic, case):
     from contextlib import ExitStack
 
     refined = case == "refined"
-    cfgdict, bodydict, p = config(h, 40 if refined else 80)
+    cfgdict, bodydict, p = config(h, 40 if refined else 80, a.collision)
     cap = 5.0 if a.mode == "pilot" else 300.0
     max_steps = math.ceil(cap * p["diameter"] / SPEED)
     rows = []
@@ -179,7 +187,7 @@ def run_case(a, h, periodic, case):
     dest.mkdir()
     with ExitStack() as stack:
         if refined:
-            cfg, pair, body, host = host_geometry(a.source, h)
+            cfg, pair, body, host = host_geometry(a.source, h, a.collision)
             context = stack.enter_context(CudaContext(0))
             backend = stack.enter_context(
                 FixedPatchStepper(
@@ -336,6 +344,8 @@ def run_case(a, h, periodic, case):
                 if len(rows) == 1:
                     ready = dict(
                         case=case,
+                        collision=a.collision,
+                        **a.provenance,
                         mode=a.mode,
                         step=begin,
                         tstar=row["tstar"],
@@ -348,6 +358,8 @@ def run_case(a, h, periodic, case):
                     print("READY " + json.dumps(ready, allow_nan=False), flush=True)
     result = dict(
         case=case,
+        collision=a.collision,
+        **a.provenance,
         mode=a.mode,
         history=rows,
         healthy=bool(healthy),
@@ -416,12 +428,14 @@ def check_deadline(deadline):
         )
 
 
-def verify_pilot(path, manifest):
+def verify_pilot(path, manifest, collision="SRT"):
     if path is None:
         raise ValueError(
             "Full comparison requires --pilot-results from a completed matched pilot"
         )
     previous = read_finite_json(path / "plan.json")
+    if previous.get("collision") != collision:
+        raise ValueError("Pilot collision differs")
     if previous["source_manifest"] != manifest or previous["bounds_root_edge"] != list(
         BOUNDS
     ):
@@ -432,7 +446,8 @@ def verify_pilot(path, manifest):
     for case, diameter in (("uniform", 80), ("refined", 40)):
         summary = read_finite_json(path / case / "summary.json")
         if (
-            summary.get("case") != case
+            summary.get("collision") != collision
+            or summary.get("case") != case
             or summary.get("mode") != "pilot"
             or summary.get("status") != "pilot_complete"
             or summary.get("healthy") is not True
@@ -483,7 +498,8 @@ def verify_pilot(path, manifest):
                 raise ValueError("Unhealthy pilot observations")
         ready = read_finite_json(path / case / "ready.json")
         if (
-            ready.get("case") != case
+            ready.get("collision") != collision
+            or ready.get("case") != case
             or ready.get("mode") != "pilot"
             or ready.get("step") != interval
             or ready.get("runner_sha256") != runner_hash()
@@ -507,6 +523,7 @@ def verify_pilot(path, manifest):
 
 def main():
     parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--collision", choices=["SRT", "MRT"], default="SRT")
     parser.add_argument("--source", type=Path, default=DEFAULT_SOURCE)
     parser.add_argument("--output", type=Path, required=True)
     parser.add_argument("--mode", choices=["plan", "pilot", "full"], default="plan")
@@ -519,9 +536,15 @@ def main():
     )
     a = parser.parse_args()
     a.source = a.source.resolve()
-    manifest = verify_source(a.source)
+    manifest = verify_source(a.source, a.collision)
     h, periodic = helpers(a.source)
-    calibration = verify_pilot(a.pilot_results, manifest) if a.mode == "full" else None
+    from source_provenance import checkout_provenance
+
+    a.provenance = checkout_provenance(a.source)
+    calibration = (
+        verify_pilot(a.pilot_results, manifest, a.collision)
+        if a.mode == "full" else None
+    )
     if a.mode != "plan" and (
         a.max_wall_seconds is None
         or not math.isfinite(a.max_wall_seconds)
@@ -536,16 +559,17 @@ def main():
     if a.output.exists() and any(a.output.iterdir()):
         raise ValueError("Output must be empty/new")
     a.output.mkdir(parents=True, exist_ok=True)
-    specification = plan(h)
+    specification = plan(h, a.collision)
     specification.update(
         source_manifest=manifest,
+        **a.provenance,
         gates=periodic.GATES,
         reference_gates=periodic.REFERENCE_GATES,
         runner_sha256=runner_hash(),
         wall_time_budget_seconds=budget,
         full_pilot_calibration_seconds=calibration,
     )
-    host_geometry(a.source, h)
+    host_geometry(a.source, h, a.collision)
     (a.output / "plan.json").write_text(json.dumps(specification, indent=2))
     if a.mode == "plan":
         return 0
@@ -561,6 +585,8 @@ def main():
         except TimeoutError as exc:
             failure = dict(
                 status="wall_time_limit",
+                collision=a.collision,
+                **a.provenance,
                 case=case,
                 matched_case_pass=False,
                 error=str(exc),
@@ -574,6 +600,8 @@ def main():
         if not result["healthy"]:
             break
     comparison = dict(
+        collision=a.collision,
+        **a.provenance,
         status="pilot_only" if a.mode == "pilot" else "incomplete",
         matched_case_pass=False,
     )

tests/validation/time_refined_dfg.py

@@ -0,0 +1,326 @@
+"""Measurement-only wall timing for the matched MRT refined DFG case. No numerical edits.
+
+Variants over EQUAL simulated duration (one refined base step = two uniform D80 steps):
+  A  refined: step(1) + read_raw_exchange per base step  -- runner-like force-read workload
+  B  refined: step(1) only, no host reads        -- timing-only, no force history
+  C  refined: one step(N) call                   -- timing-only; ledger keeps final base interval only
+  E  uniform: sim.run(<=100) + read_force                -- runner-like force-read workload
+  E' uniform: one sim.run(2N) call                -- timing-only; at most a whole-window force
+                                                    average (different summation order), no history
+A/E read forces but discard them: no retained force history, no interval aggregation/validation,
+and no macro/probe/CSV cost, so they are not full benchmark throughput.
+observe_raw_exchange stays True for A/B/C: B/C omit host reads only.
+Status lines: CONSTRUCTED = backend/simulation built (JIT done); WARM_READY = first warm-up and
+health check passed for that case.
+Parity compares A/B/C with each other and E/E' with each other, never uniform with refined.
+A-B and B-C are measured comparative differences, not exact causal decompositions.
+Each repeat starts from the same fresh initial state with the same warm-up; construction,
+compilation, parity and health checks are outside timed regions. The actual stepping stream is
+synchronized explicitly immediately before the timer starts and after the timed loop. Wall time
+(perf_counter) is the only timing source; no CUDA-event span is reported as kernel time. An outer
+--max-wall-seconds deadline is checked between variants and inside per-step loops.
+"""
+
+import argparse
+import json
+import math
+import os
+from pathlib import Path
+import statistics
+import sys
+import time
+from tempfile import TemporaryDirectory
+
+import numpy as np
+
+HERE = Path(__file__).resolve().parent
+sys.path.insert(0, str(HERE))
+import run_refined_dfg as r  # noqa: E402  (published runner: same setup functions)
+
+VARIANTS = {
+    "A": dict(case="refined", force_history_valid=False, force_reads=True,
+              force_note="runner-like force-read workload; values discarded, no aggregation"),
+    "B": dict(case="refined", force_history_valid=False, timing_only=True,
+              force_note="no host reads"),
+    "C": dict(case="refined", force_history_valid=False, timing_only=True,
+              force_note="raw ledger retains only the final base interval"),
+    "E": dict(case="uniform", force_history_valid=False, force_reads=True,
+              force_note="runner-like force-read workload; values discarded, no aggregation"),
+    "E'": dict(case="uniform", force_history_valid=False, timing_only=True,
+               force_note="at most a whole-window accumulated average; different summation order"),
+}
+
+
+def schedule(n_base, warmup_base):
+    """Steps per case for equal simulated duration (refined dt_root = 2 dt_uniform)."""
+    if n_base <= 0 or warmup_base < 0:
+        raise ValueError("n_base must be positive and warmup_base nonnegative")
+    return dict(refined=(warmup_base, n_base), uniform=(2 * warmup_base, 2 * n_base),
+                tstar_window=n_base * r.SPEED / 40.0)
+
+
+# Existing run_refined_dfg safety limits (rho in [.5, 1.5], sampled max Mach <= .2).
+RHO_LIMITS, MAX_MACH = (0.5, 1.5), 0.2
+
+
+def same_state(reference, values):
+    """Bitwise populations and identical NaN masks; returns a short description or None.
+
+    Matching NaN is accepted only as a storage mask (inactive/ghost/solid slots); physical
+    leaf health is checked separately. Any infinity in either state is a mismatch.
+    """
+    reference, values = list(reference), list(values)
+    if len(reference) != len(values):
+        return f"level count {len(reference)} != {len(values)}"
+    for level, (a, b) in enumerate(zip(reference, values)):
+        a, b = np.asarray(a), np.asarray(b)
+        if a.shape != b.shape:
+            return f"level {level}: shape {a.shape} != {b.shape}"
+        if np.isinf(a).any() or np.isinf(b).any():
+            return f"level {level}: infinite population"
+        if not np.array_equal(np.isnan(a), np.isnan(b)):
+            return f"level {level}: NaN mask differs"
+        finite = np.isfinite(a)
+        if a.dtype != b.dtype or not np.array_equal(a[finite].view(np.uint8), b[finite].view(np.uint8)):
+            return f"level {level}: finite bits differ"
+    return None
+
+
+def uniform_populations_finite(ddf, fluid):
+    """Simulation.get_ddf() may be flat (9*ny*nx); view it as (Q, cells) before masking."""
+    ddf, fluid = np.asarray(ddf), np.asarray(fluid, bool).ravel()
+    if fluid.size == 0 or ddf.size % fluid.size:
+        raise RuntimeError(f"DDF size {ddf.size} incompatible with {fluid.size} cells")
+    return bool(np.isfinite(ddf.reshape(ddf.size // fluid.size, fluid.size)[:, fluid]).all())
+
+
+def physical_health(fields):
+    """fields: (rho, ux, uy) arrays already restricted to physical leaves/fluid."""
+    if not fields or any(not np.isfinite(v).all() or v.size == 0 for f in fields for v in f):
+        raise RuntimeError("nonfinite or empty physical macroscopic field")
+    rho_min = float(min(f[0].min() for f in fields))
+    rho_max = float(max(f[0].max() for f in fields))
+    max_mach = float(max(np.sqrt(3 * (f[1] ** 2 + f[2] ** 2)).max() for f in fields))
+    if not (RHO_LIMITS[0] <= rho_min <= rho_max <= RHO_LIMITS[1] and max_mach <= MAX_MACH):
+        raise RuntimeError(f"health limits violated: rho [{rho_min}, {rho_max}], Mach {max_mach}")
+    return dict(rho_min=rho_min, rho_max=rho_max, max_mach=max_mach)
+
+
+def refined_health(backend, masks):
+    arrays = backend.download()  # physical FP64 populations; non-leaves NaN by contract
+    if any(not np.isfinite(g[:, m]).all() for g, m in zip(arrays, masks)):
+        raise RuntimeError("nonfinite physical leaf population")
+    return physical_health([tuple(v[m] for v in r.macros(g)) for g, m in zip(arrays, masks)])
+
+
+def uniform_health(sim, fluid):
+    m = sim.get_macroscopic()
+    return physical_health([tuple(np.asarray(m[k], np.float64).reshape(fluid.shape)[fluid]
+                                  for k in ("rho", "ux", "uy"))])
+
+
+def check_deadline(deadline):
+    if time.monotonic() >= deadline:
+        raise TimeoutError("timing harness exceeded --max-wall-seconds")
+
+
+def sync_refined(backend):
+    with backend.ctx.activate():
+        backend._stream.synchronize()
+
+
+def sync_uniform(sim):
+    with sim.ctx.activate():
+        sim.stream.synchronize()
+
+
+def refined_loop(backend, variant, n, deadline):
+    if variant == "C":
+        backend.step(n)
+        return
+    for i in range(n):
+        backend.step(1)
+        if variant == "A":
+            backend.read_raw_exchange(0, level=1, window="last_base_step")
+        if i % 100 == 99:
+            check_deadline(deadline)
+
+
+def uniform_loop(sim, variant, n, deadline):
+    if variant == "E'":
+        sim.run(n)
+        return
+    done = 0
+    while done < n:
+        batch = min(100, n - done)
+        sim.run(batch)
+        sim.read_force(0)
+        done += batch
+        check_deadline(deadline)
+
+
+def time_refined(a, h, sched):
+    from contextlib import ExitStack
+    from CelerisLab.cuda.context import CudaContext
+    from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+    from CelerisLab.lbm.channel import channel_masks
+    from CelerisLab.lbm.refinement import ACTIVE
+    from CelerisLab.lbm.descriptors import SOLID
+
+    warm, n = sched["refined"]
+    rows, reference = [], {}
+    with ExitStack() as stack:
+        cfg, pair, body, host = r.host_geometry(a.source, h, "MRT")
+        context = stack.enter_context(CudaContext(0))
+        backend = stack.enter_context(FixedPatchStepper(
+            cfg, r.BOUNDS, context=context, coarse_boundary="channel", base_block=r.BASE_BLOCK,
+            fine_block=r.FINE_BLOCK, fine_objects=(body,), observe_raw_exchange=True))
+        masks = [channel_masks(pair.coarse_shape)[0] & ~pair.covered,  # as run_refined_dfg
+                 (pair.fine_kind == ACTIVE) & ((host.flags.reshape(pair.fine_shape) & SOLID) == 0)]
+        fresh = (r.initial_state(pair.coarse_shape, (0.0, 0.0), 1.0, root=True),
+                 r.initial_state(pair.fine_shape, pair.fine_origin, 0.5))
+        print(json.dumps(dict(status="CONSTRUCTED", case="refined")), flush=True)
+        for repeat in range(a.repeats):
+            for variant in ("A", "B", "C"):
+                check_deadline(a.deadline)
+                backend.initialize(*fresh)  # resets step count; same fresh state every time
+                refined_loop(backend, variant, warm, a.deadline)
+                warm_health = refined_health(backend, masks)  # untimed
+                if repeat == 0 and variant == "A":
+                    print(json.dumps(dict(status="WARM_READY", case="refined")), flush=True)
+                sync_refined(backend)
+                start = time.perf_counter()
+                refined_loop(backend, variant, n, a.deadline)
+                sync_refined(backend)
+                seconds = time.perf_counter() - start
+                end_health = refined_health(backend, masks)  # untimed
+                values = [np.array(v, copy=True) for v in backend._download_raw().values]
+                if not all(np.isfinite(v[:, m]).all() for v, m in zip(values, backend._state)):
+                    raise RuntimeError(f"refined {variant} repeat {repeat}: nonfinite leaf populations")
+                mismatch = None
+                if repeat == 0 and variant == "A":
+                    reference["A"] = values
+                else:
+                    mismatch = same_state(reference["A"], values)
+                rows.append(dict(variant=variant, repeat=repeat, seconds=seconds,
+                                 base_steps=n, step_count=backend.step_count,
+                                 health_after_warmup=warm_health, health_after_window=end_health,
+                                 parity_vs_A=mismatch or "identical", **VARIANTS[variant]))
+                if mismatch:
+                    raise RuntimeError(f"refined {variant} repeat {repeat}: {mismatch}")
+    return rows
+
+
+def time_uniform(a, h, sched):
+    from CelerisLab.simulation import Simulation
+    from CelerisLab.lbm.descriptors import FLUID
+
+    warm, n = sched["uniform"]
+    rows, reference = [], {}
+    cfgdict, bodydict, params = r.config(h, 80, "MRT")
+    with TemporaryDirectory() as tmp:
+        tmp = Path(tmp)
+        (tmp / "lbm.json").write_text(json.dumps(cfgdict))
+        (tmp / "body.json").write_text(json.dumps(bodydict))
+        with Simulation(str(tmp / "lbm.json"), str(tmp / "body.json"), device_id=0) as sim:
+            sim.initialize()
+            sim.snapshot()  # fresh state; restore keeps clock and streaming parity
+            fluid = (sim.get_flags() & FLUID) != 0
+            ny, nx = params["ny"], params["nx"]
+            fluid = fluid.reshape(ny, nx)
+            print(json.dumps(dict(status="CONSTRUCTED", case="uniform")), flush=True)
+            for repeat in range(a.repeats):
+                for variant in ("E", "E'"):
+                    check_deadline(a.deadline)
+                    sim.restore()
+                    if sim.stepper.step_count % 2:
+                        raise RuntimeError("uniform restore parity changed between variants")
+                    uniform_loop(sim, variant, warm, a.deadline)
+                    warm_health = uniform_health(sim, fluid)  # untimed
+                    if repeat == 0 and variant == "E":
+                        print(json.dumps(dict(status="WARM_READY", case="uniform")), flush=True)
+                    sync_uniform(sim)
+                    start = time.perf_counter()
+                    uniform_loop(sim, variant, n, a.deadline)
+                    sync_uniform(sim)
+                    seconds = time.perf_counter() - start
+                    end_health = uniform_health(sim, fluid)  # untimed
+                    values = [sim.get_ddf()]
+                    if not uniform_populations_finite(values[0], fluid):
+                        raise RuntimeError(f"uniform {variant} repeat {repeat}: nonfinite fluid populations")
+                    mismatch = None
+                    if repeat == 0 and variant == "E":
+                        reference["E"] = values
+                    else:
+                        mismatch = same_state(reference["E"], values)
+                    rows.append(dict(variant=variant, repeat=repeat, seconds=seconds,
+                                     uniform_steps=n, parity_vs_E=mismatch or "identical",
+                                     health_after_warmup=warm_health, health_after_window=end_health,
+                                     **VARIANTS[variant]))
+                    if mismatch:
+                        raise RuntimeError(f"uniform {variant} repeat {repeat}: {mismatch}")
+    return rows
+
+
+def summarize(rows, tstar_window):
+    summary = {}
+    for variant in VARIANTS:
+        seconds = [row["seconds"] for row in rows if row["variant"] == variant]
+        if seconds:
+            median = statistics.median(seconds)
+            summary[variant] = dict(median_seconds=median, seconds=seconds,
+                                    seconds_per_tstar=median / tstar_window,
+                                    **VARIANTS[variant])
+    if {"A", "B", "C"} <= summary.keys():
+        summary["A_minus_B_seconds"] = summary["A"]["median_seconds"] - summary["B"]["median_seconds"]
+        summary["B_minus_C_seconds"] = summary["B"]["median_seconds"] - summary["C"]["median_seconds"]
+    if {"A", "E"} <= summary.keys():
+        summary["A_over_E"] = summary["A"]["median_seconds"] / summary["E"]["median_seconds"]
+    if {"C", "E'"} <= summary.keys():
+        summary["C_over_Eprime"] = summary["C"]["median_seconds"] / summary["E'"]["median_seconds"]
+    summary["interpretation"] = ("Measured comparative differences under GPU variability; "
+                                 "not exact causal decompositions or performance claims.")
+    return summary
+
+
+def main():
+    p = argparse.ArgumentParser(description=__doc__)
+    p.add_argument("--source", type=Path, default=r.DEFAULT_SOURCE)
+    p.add_argument("--output", type=Path, required=True)
+    p.add_argument("--base-steps", type=int, default=2000)
+    p.add_argument("--warmup-base-steps", type=int, default=200)
+    p.add_argument("--repeats", type=int, default=3)
+    p.add_argument("--max-wall-seconds", type=float, help="Required outer timeout for GPU runs")
+    p.add_argument("--plan", action="store_true")
+    a = p.parse_args()
+    if a.repeats < 1 or a.base_steps > 20000:
+        p.error("repeats must be >= 1 and base-steps <= 20000")
+    a.source = a.source.resolve()
+    manifest = r.verify_source(a.source, "MRT")
+    sched = schedule(a.base_steps, a.warmup_base_steps)
+    plan = dict(collision="MRT", manifest_commit=manifest["commit"], schedule=sched,
+                repeats=a.repeats, variants=VARIANTS, timer="time.perf_counter wall",
+                excluded="construction/JIT, warm-up, state download and parity checks")
+    if a.plan:
+        print(json.dumps(plan, indent=2))
+        return 0
+    if a.max_wall_seconds is None or not math.isfinite(a.max_wall_seconds) or a.max_wall_seconds <= 0:
+        p.error("GPU timing requires a finite positive --max-wall-seconds")
+    a.deadline = time.monotonic() + a.max_wall_seconds
+    if not os.environ.get("CUDA_VISIBLE_DEVICES") or "," in os.environ["CUDA_VISIBLE_DEVICES"]:
+        p.error("Root must select exactly one verified idle GPU via CUDA_VISIBLE_DEVICES")
+    if a.output.exists() and any(a.output.iterdir()):
+        p.error("Output must be empty/new")
+    a.output.mkdir(parents=True, exist_ok=True)
+    h, _ = r.helpers(a.source)
+    rows = time_uniform(a, h, sched) + time_refined(a, h, sched)
+    result = dict(plan=plan, rows=rows, summary=summarize(rows, sched["tstar_window"]))
+    if not all(math.isfinite(row["seconds"]) for row in rows):
+        raise RuntimeError("Nonfinite timing")
+    (a.output / "timing.json").write_text(json.dumps(result, indent=2))
+    print(json.dumps(result["summary"], indent=2), flush=True)
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())

tests/validation/time_window_dfg.py

@@ -0,0 +1,183 @@
+"""Measurement-only A-versus-window timing for the matched shifted-MRT refined DFG setup.
+
+Candidate source (window ledger) is recorded by checkout revision, dirty flag and file hashes;
+the frozen refined-DFG manifests/results are NOT used or modified (they pin pre-window source).
+Equal simulated duration, same fresh state and warm-up per variant, one observation interval =
+80 root base steps (the runner's dense-phase interval):
+  A   per base step: step(1) + read_raw_exchange, summed per interval (current runner workload)
+  W   per interval: start_exchange_window(); step(80); read_exchange_window()
+  W1  per interval: start_exchange_window(); 80 x step(1); read_exchange_window()
+All variants keep observe_raw_exchange=True; the backend is built once with
+observe_window_exchange=True (a closed window adds no launch). Populations (bitwise) and
+per-interval total impulses are compared OUTSIDE timing; the impulse tolerance is a
+conservative engineering tolerance (see INTERVAL_RTOL); the tighter substep-scaled engineering
+tolerance is exercised by the capability test tests/validation/multigrid/test_exchange_window.py.
+Wall time (perf_counter) with explicit stream synchronization is the only timing source.
+"""
+
+import argparse
+import hashlib
+import json
+import math
+import os
+from pathlib import Path
+import statistics
+import sys
+import time
+
+import numpy as np
+
+HERE = Path(__file__).resolve().parent
+sys.path.insert(0, str(HERE))
+import time_refined_dfg as t  # noqa: E402  (reuses run_refined_dfg setup via t.r)
+
+r = t.r
+INTERVAL = 80
+VARIANTS = ("A", "W", "W1")
+# Interval impulses differ only by FP64 summation order (per link over time vs per step over
+# links). Engineering tolerance: 1e-9 of the interval's summed absolute per-step totals.
+INTERVAL_RTOL = 1e-9
+CANDIDATE_FILES = ("src/CelerisLab/lbm/refinement_stepper.py", "src/CelerisLab/lbm/raw_exchange.py",
+                   "src/CelerisLab/lbm/refinement_tables.py", "src/CelerisLab/lbm/kernels/blocks.cu")
+
+
+def components(obs):
+    o = obs["total"]
+    return np.array([o["impulse"][0], o["impulse"][1], o["angular_impulse"], o["mass_into_fluid"]])
+
+
+def candidate_identity(source):
+    from source_provenance import checkout_provenance
+    return dict(checkout_provenance(source),
+                sha256={f: hashlib.sha256((source / f).read_bytes()).hexdigest()
+                        for f in CANDIDATE_FILES})
+
+
+def run_variant(backend, variant, n_intervals, deadline):
+    """Returns per-interval totals and, for A, per-interval summed |per-step totals|."""
+    totals, scales = [], []
+    for _ in range(n_intervals):
+        if variant == "A":
+            impulse, scale = np.zeros(4), np.zeros(4)
+            for _ in range(INTERVAL):
+                backend.step(1)
+                c = components(backend.read_raw_exchange(0, level=1, window="last_base_step"))
+                impulse += c
+                scale += np.abs(c)
+            totals.append(impulse)
+            scales.append(scale)
+        else:
+            backend.start_exchange_window()
+            if variant == "W":
+                backend.step(INTERVAL)
+            else:
+                for _ in range(INTERVAL):
+                    backend.step(1)
+            totals.append(components(backend.read_exchange_window(0)))
+        t.check_deadline(deadline)
+    return np.array(totals), (np.array(scales) if scales else None)
+
+
+def compare_intervals(reference, scales, totals):
+    """Max violation ratio; <= 1 passes. Scale includes a floor so zero cancels cannot divide."""
+    bound = INTERVAL_RTOL * scales + 1e-300
+    return float(np.max(np.abs(totals - reference) / bound))
+
+
+def main():
+    p = argparse.ArgumentParser(description=__doc__)
+    p.add_argument("--source", type=Path, default=r.DEFAULT_SOURCE)
+    p.add_argument("--output", type=Path, required=True)
+    p.add_argument("--base-steps", type=int, default=2000)
+    p.add_argument("--warmup-base-steps", type=int, default=200)
+    p.add_argument("--repeats", type=int, default=3)
+    p.add_argument("--max-wall-seconds", type=float)
+    p.add_argument("--plan", action="store_true")
+    a = p.parse_args()
+    if a.base_steps <= 0 or a.base_steps % INTERVAL or a.base_steps > 20000 or a.repeats < 1:
+        p.error(f"base-steps must be a positive multiple of {INTERVAL} (<= 20000); repeats >= 1")
+    a.source = a.source.resolve()
+    identity = candidate_identity(a.source)
+    plan = dict(collision="MRT", shifted_storage="as run_refined_dfg MRT config",
+                candidate_source=identity, frozen_manifests_used=False, interval_base_steps=INTERVAL,
+                base_steps=a.base_steps, warmup_base_steps=a.warmup_base_steps, repeats=a.repeats,
+                variants=VARIANTS, tstar_window=a.base_steps * r.SPEED / 40.0,
+                timer="time.perf_counter wall with explicit stream sync",
+                excluded="construction/JIT, warm-up, health, downloads, parity/impulse comparison")
+    if a.plan:
+        print(json.dumps(plan, indent=2))
+        return 0
+    if a.max_wall_seconds is None or not math.isfinite(a.max_wall_seconds) or a.max_wall_seconds <= 0:
+        p.error("GPU timing requires a finite positive --max-wall-seconds")
+    if not os.environ.get("CUDA_VISIBLE_DEVICES") or "," in os.environ["CUDA_VISIBLE_DEVICES"]:
+        p.error("Root must select exactly one verified idle GPU via CUDA_VISIBLE_DEVICES")
+    if a.output.exists() and any(a.output.iterdir()):
+        p.error("Output must be empty/new")
+    a.output.mkdir(parents=True, exist_ok=True)
+    deadline = time.monotonic() + a.max_wall_seconds
+    h, _ = r.helpers(a.source)
+
+    from contextlib import ExitStack
+    from CelerisLab.cuda.context import CudaContext
+    from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+    from CelerisLab.lbm.channel import channel_masks
+    from CelerisLab.lbm.refinement import ACTIVE
+    from CelerisLab.lbm.descriptors import SOLID
+
+    rows, ref = [], {}
+    n_intervals = a.base_steps // INTERVAL
+    with ExitStack() as stack:
+        cfg, pair, body, host = r.host_geometry(a.source, h, "MRT")
+        context = stack.enter_context(CudaContext(0))
+        backend = stack.enter_context(FixedPatchStepper(
+            cfg, r.BOUNDS, context=context, coarse_boundary="channel", base_block=r.BASE_BLOCK,
+            fine_block=r.FINE_BLOCK, fine_objects=(body,), observe_raw_exchange=True,
+            observe_window_exchange=True))
+        masks = [channel_masks(pair.coarse_shape)[0] & ~pair.covered,
+                 (pair.fine_kind == ACTIVE) & ((host.flags.reshape(pair.fine_shape) & SOLID) == 0)]
+        fresh = (r.initial_state(pair.coarse_shape, (0.0, 0.0), 1.0, root=True),
+                 r.initial_state(pair.fine_shape, pair.fine_origin, 0.5))
+        print(json.dumps(dict(status="CONSTRUCTED", case="refined")), flush=True)
+        for repeat in range(a.repeats):
+            for variant in VARIANTS:
+                t.check_deadline(deadline)
+                backend.initialize(*fresh)
+                for _ in range(a.warmup_base_steps):
+                    backend.step(1)
+                warm = t.refined_health(backend, masks)
+                if repeat == 0 and variant == "A":
+                    print(json.dumps(dict(status="WARM_READY", case="refined")), flush=True)
+                t.sync_refined(backend)
+                start = time.perf_counter()
+                totals, scales = run_variant(backend, variant, n_intervals, deadline)
+                t.sync_refined(backend)
+                seconds = time.perf_counter() - start
+                end = t.refined_health(backend, masks)
+                values = [np.array(v, copy=True) for v in backend._download_raw().values]
+                if variant == "A" and repeat == 0:
+                    ref.update(values=values, totals=totals, scales=scales)
+                    parity, ratio = "identical", 0.0
+                else:
+                    parity = t.same_state(ref["values"], values) or "identical"
+                    ratio = compare_intervals(ref["totals"], ref["scales"], totals)
+                rows.append(dict(variant=variant, repeat=repeat, seconds=seconds,
+                                 intervals=n_intervals, population_parity_vs_A=parity,
+                                 interval_impulse_violation_ratio=ratio,
+                                 health_after_warmup=warm, health_after_window=end))
+                if parity != "identical" or not ratio <= 1:
+                    raise RuntimeError(f"{variant} repeat {repeat}: parity {parity}, impulse ratio {ratio}")
+    summary = {v: dict(median_seconds=statistics.median(x["seconds"] for x in rows if x["variant"] == v),
+                       seconds=[x["seconds"] for x in rows if x["variant"] == v]) for v in VARIANTS}
+    summary.update(W_over_A=summary["W"]["median_seconds"] / summary["A"]["median_seconds"],
+                   W1_over_A=summary["W1"]["median_seconds"] / summary["A"]["median_seconds"],
+                   interpretation="Measured comparative timings under GPU variability; W also "
+                                  "batches step calls, W1 isolates the read pattern. Not a "
+                                  "benchmark-throughput or causal claim.")
+    (a.output / "window_timing.json").write_text(json.dumps(dict(plan=plan, rows=rows, summary=summary),
+                                                           indent=2))
+    print(json.dumps(summary, indent=2), flush=True)
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())