test: compare fixed centered reconstruction on audited square faces

docs/validation/taylor_green.md

@@ -459,3 +459,148 @@ What's kept and enforced:
   a 600 kB cap. If invalid it is dropped and an error is recorded, keeping
   the scalar receipt.
 - There are no new pass gates; this is a mechanism diagnostic only.
+
+### Single x-gap-width contrast (`wide_x_shear_n128`)
+
+This opt-in admits only Re100/N128/projected depth3 with `--probe-trace`.
+The outer boxes, three grid spacings, viscosities, initialization, five main
+sample times and trace windows are unchanged. L1→L2 parent-frame rectangles
+are `[(56,44,96,212),(160,60,200,228)]`: both x-normal middle gaps widen from
+6 to 12 root cells, while y gaps remain 6. Active leaves are
+`[7936,20352,53760]`, hops `[2048,4096,8192]`, work 540,016,640 active updates.
+This also moves the L2 x faces and reduces its volume; it is not a pure
+travel-distance experiment. No new qualification matrix, assessor, threshold
+or snapshot is added. Existing scalar screens and all failures stay visible.
+
+The same 24 targets and actual-center metadata are retained. The horizontal
+probes stop before the new inner face x = 28: they do not bracket that interface
+or sample its fine side. The first 8 horizontal L0/L1 probes keep identical
+centers, and the y-normal x = 38 transect is unchanged. Their temporal residuals
+can support a descriptive timing contrast, but cannot establish a crossing
+or reflection coefficient. Changed-level probes must not be treated as
+collocated samples. `rho-1` remains raw density deviation; do not substitute
+the 24-probe mean for the global pressure gauge.
+
+### POST-G initialization counterfactual (`postg_shear_n128`)
+
+Exactly three fresh Re100 cases are admitted: U128, U512, and projected R128
+depth 3 with the original nested rectangles. Each requires `--flow shear_wave
+--diagnostic postg_shear_n128 --probe-trace`; snapshots and early-pressure mode
+are excluded. Geometry, physical velocity/density, viscosities, collisions,
+five main times and the 24 signed probes/windows are unchanged. Total work is
+2,952,265,728 active updates, with 15 main rows, 774 trace rows and five builds.
+
+Only initialization changes. At level spacing h, k_level = h k and
+A_i = 3 w_i U cos(theta) (c_i dot e) (c_i dot k_level). Initialize physical
+FP64 **POST G = feq - (tau - 1) A**, then use the existing shifted FP32 codec.
+This is the first-gradient, low-Mach Chapman-Enskog approximation, not an
+exact discrete eigenstate. The PRE coefficient is -tau; the physical viscous
+stress coefficient tau - 1/2 must not replace the raw POST-G coefficient.
+The correction has zero mass/momentum before encoding; FP32 introduces normal
+roundoff. U512 uses its native k and tau, matching nested L2 physical spacing.
+
+Receipts pin the new initializer module and identify `postg_first_gradient_v1`,
+with `qualification_eligible=false`. Existing scalar screens remain visible
+as descriptive outcomes, including failures; no new pass gate or assessor is
+introduced. Historical assessors reject the new diagnostic and initialization
+identities. Compare raw rows and signed traces with the retained equilibrium
+receipts as a startup-state counterfactual. The host first-interval pressure
+RMS reduction is not a prediction of long-time improvement, and individual
+probe errors can increase. All equilibrium paths and historical verdicts stay
+unchanged.
+
+### Short square orientation diagnostic (`axis_square_shear_n128`)
+
+Four fresh cases only: U128, U256, constant R128 and projected R128 depth 2.
+Use `--flow shear_wave --diagnostic axis_square_shear_n128 --probe-trace`.
+The supported union is `[[(32,32,96,96)]]`; a full-height periodic fine strip
+is not supported by the unchanged production geometry. All four cases use
+POST G = feq - (tau - 1) A, mode (2,0), polarization (0,1), phase pi/4,
+U=.03 and Re100. Root/fine nu are .0384/.0768 and tau are .6152/.7304.
+This first-gradient low-Mach initialization is not an exact discrete mode.
+
+Stop at root step 8: two global rows at 0/8, nine local rows at 0..8;
+U256 advances two native steps per root step. Refined leaves are [12288,16384],
+hops [8,16], work 360,448 per refined case. Four-case work is 1,900,544 ACTIVE
+updates, six build sidecars and 720 native cell-time samples. Existing runtime,
+build/source guards and receipt/trace caps apply. No snapshots or late fit.
+`screening_pass=null`, `qualification_eligible=false`: exit zero means valid
+finite completion, with no physical pass threshold or legacy assessor.
+
+Eight root-unit-cell groups (i,j): inner L/R/B/T = (32,64)/(95,63)/(64,32)/(63,95),
+outer L/R/B/T = (31,64)/(96,63)/(64,31)/(63,96). Refined inner groups contain four
+fine cells, outer groups one coarse cell; U128 uses one and U256 four per group.
+Retain raw rho/ux/uy, native indices, physical centers/areas and parity. Group
+residuals are signed means of rho-1 and rho*u-u_exact, evaluated at each actual
+center; normalization uses that group's mean analytic d(uy)/dx. Only inner
+L/B and R/T pairs have exactly opposite equal-magnitude gradients. Controls
+have their own native quadrature. Probe density is ungauged; global pressure
+RMS uses one domain gauge.
+
+Actual-table local corridors match tangentially extended straight faces beyond
+22.5 root cells. The conservative paired-map radius is 2.5 root cells per root
+step, so these selected groups through step 8 are locally straight-face
+equivalent. This is not a claim about global corner errors. For projected
+explosion with actual missing-normal-gradient fallback, the stationary
+linearized affine reference after the first complete root interval predicts
+zero conserved-moment error on vertical faces. At the
+bottom/top inner fine quartets its mean delta jx / d(uy)/dx is
+-/+ 1/(576*tau_f), with opposite adjacent-coarse redistribution.
+Keep raw values, step-zero rounding and fresh controls: sinusoidal curvature
+and finite-Mach effects prevent treating this coefficient as an exact gate.
+
+### Test-only centered straight-face comparison
+
+`centered_face_square_n128` reuses the square axis-wave POST-G initializer,
+geometry, native probes and root steps 0–8. It has four fresh cases: U128,
+U256, current projected R128, and experimental centered R128. This is a
+qualification-ineligible diagnostic; there is no legacy assessor or new physics gate.
+
+Only the candidate's private generated fine-level `blocks.cu` changes. A
+2,112-byte constant mask selects 1,792 rows (448 complete donor quartets).
+The diagonal projected correction is doubled on straight faces with exactly
+one unsupported normal axis and both legal tangent neighbors. Donor centers
+within four root cells of a corner retain the current rule, as do all other
+patterns, axial populations and rest. The mask is symmetric under transpose
+and reflection. Its hash and the actual modified source/PTX hashes are retained;
+production files, restrictions, collision, scheduling and launch arguments stay fixed.
+A complete quartet has zero summed offset, so each population's quartet mean
+is preserved. The perpendicular gradient also retains equal characteristic
+packets; the existing two-hop conservation routing therefore remains applicable.
+Host tests additionally check the mixed mask through complete two-hop/collision
+cycles, including corner grazing, and execute the actual edited kernel via C++.
+The central probes are more than 20 root cells from mask transitions; the
+reviewed conservative dependency radius is 2.5 root cells per root interval.
+
+The separate comparison subtracts each group's step-zero residual and the
+fresh uniform control at matching quadrature (U256 inner, U128 outer). It retains
+signed density and both momentum components, actual-center gradient normalization,
+parities, and sign flips. Step 1 requires bottom-inner cross-momentum magnitude
+improvement greater than its deadband, and top-inner nonworsening within it.
+Each horizontal face's step-8 RMS combines its inner and outer group means;
+both must be nonworse. All four vertical groups must remain unchanged within
+the component deadbands at every sampled root step. The norm
+`sqrt(drho²/3 + djx² + djy²)` must not worsen beyond its corresponding deadband.
+Sampled raw-pressure envelope and step-8 global gauged pressure RMS must not
+worsen beyond the pressure deadband. Nonfinite/nonpositive sampled density,
+invalid traces, source/build mismatch, or exceeded observed population envelope
+reject the comparison. A small top signal or sign reversal remains explicit.
+
+The predeclared shifted population envelope is 1/32. From per-population oracle
+allowance `6e-8 + 2e-7*1/32`, pairwise moment deadbands are
+`[1.1925e-6, 7.95e-7, 7.95e-7]`; pressure uses the density deadband divided by 3.
+Area averaging cannot enlarge this local allowance. These are practical diagnostic
+deadbands, not certified accumulated root-8 rounding bounds. The envelope is
+checked on downloaded synchronized ACTIVE populations, not hidden substeps.
+Known quadratic fine-cell O(h) terms, corner transitions, repeated coupling and
+stability remain unresolved. Passing permits review of a candidate, not promotion.
+
+Run the existing `run_case.py` with `--diagnostic centered_face_square_n128
+--flow shear_wave --re 100 --probe-trace --test-reconstruction current` for
+U128/U256 (depth 1, constant) and current R128 (depth 2, projected_linear).
+For candidate R128 use `--test-reconstruction centered_face_v1`. All require
+an exact clean `--repo`, `--expected-revision`, fresh `--output`, and either
+`--dry-run` or explicit single-device `--run-gpu`. Compare ordered receipts with
+`shear_centered.py --repo REPO --expected-revision SHA --output comparison.json
+U128.json U256.json current.json candidate.json`. Total work is 1,900,544 ACTIVE
+updates, six build sidecars, and 720 cell-time samples; scalar cap stays 1 MiB/case.

tests/validation/taylor_green/run_case.py

@@ -1,6 +1,7 @@
 """One bounded GPU case, with explicit opt-in experimental explosion."""
 
 import argparse
+from contextlib import nullcontext
 import hashlib
 import json
 import os
@@ -25,7 +26,15 @@ NESTED_SHEAR = "nested_shear_n128"  # fixed three-level/two-region shear preset
 SHEAR_NESTED_SHA256 = "fa6cfe87fb4df8c6d732c637c885276833c81973b950c7724d734f8fc350bcb0"
 TWIN_SHEAR = "twin_shear_n128"  # depth-2 twin: nested L0->L1 boxes without L2
 SHEAR_TWIN_SHA256 = "1cbf37f18e89aea577581631c4565f51e68d22a3739b698816ef64acbf2b2bb4"
-SHEAR_PROBES_SHA256 = "53ba8af8ac9c8078c40f84488c7ade4c11ccc91dd1fc275ff806042e9af30166"  # signed interface-probe mechanism diagnostic
+WIDE_X_SHEAR = "wide_x_shear_n128"  # one projected x-gap-width contrast
+SHEAR_WIDE_X_SHA256 = "13d16ba355117e2d074edccd9dc5d694c2fc13c673a0abdfca6401f6ed917733"
+AXIS_SQUARE = "axis_square_shear_n128"  # short four-face orientation diagnostic
+SHEAR_AXIS_SQUARE_SHA256 = "0baffa6eb99c268ce998d90300c993d45bb3b7b27c524b636874d75808bd0a8c"
+CENTERED_SQUARE = "centered_face_square_n128"  # private-build reconstruction comparison
+SHEAR_CENTERED_SHA256 = "78bd7360e266e5e28ef4f03496413375551fcdc10ceba61d43bd88c33b4a8f5c"
+POSTG_SHEAR = "postg_shear_n128"  # initialization-only three-case counterfactual
+SHEAR_POSTG_SHA256 = "9c909fa013c092490ee21532e3125b2a78a4c5102c7a1bdc0e76986d3f0ef903"
+SHEAR_PROBES_SHA256 = "e2e6490d281cc275f448dc6d6ea7b9b4407cc74a51e9256b86c08c8626b89b76"  # signed interface-probe mechanism diagnostic
 
 
 def flow_functions(flow="taylor_green"):
@@ -121,6 +130,34 @@ def verify_source(repo, revision, flow="taylor_green", probes=False):
                 or hashlib.sha256(module.read_bytes()).hexdigest() != SHEAR_PROBES_SHA256):
             raise RuntimeError("Wrong shear_probes module checkout or content")
         result["shear_probes_sha256"] = SHEAR_PROBES_SHA256
+    if flow == "shear_centered":
+        import shear_centered
+        module = harness / "shear_centered.py"
+        if (Path(shear_centered.__file__).resolve() != module
+                or hashlib.sha256(module.read_bytes()).hexdigest() != SHEAR_CENTERED_SHA256):
+            raise RuntimeError("Wrong shear_centered module checkout or content")
+        result["shear_centered_sha256"] = SHEAR_CENTERED_SHA256
+    if flow in ("shear_axis_square", "shear_centered"):
+        import shear_axis_square
+        module = harness / "shear_axis_square.py"
+        if (Path(shear_axis_square.__file__).resolve() != module
+                or hashlib.sha256(module.read_bytes()).hexdigest() != SHEAR_AXIS_SQUARE_SHA256):
+            raise RuntimeError("Wrong shear_axis_square module checkout or content")
+        result["shear_axis_square_sha256"] = SHEAR_AXIS_SQUARE_SHA256
+    if flow == "shear_postg":
+        import shear_postg
+        module = harness / "shear_postg.py"
+        if (Path(shear_postg.__file__).resolve() != module
+                or hashlib.sha256(module.read_bytes()).hexdigest() != SHEAR_POSTG_SHA256):
+            raise RuntimeError("Wrong shear_postg module checkout or content")
+        result["shear_postg_sha256"] = SHEAR_POSTG_SHA256
+    if flow == "shear_wide_x":
+        import shear_wide_x
+        module = harness / "shear_wide_x.py"
+        if (Path(shear_wide_x.__file__).resolve() != module
+                or hashlib.sha256(module.read_bytes()).hexdigest() != SHEAR_WIDE_X_SHA256):
+            raise RuntimeError("Wrong shear_wide_x module checkout or content")
+        result["shear_wide_x_sha256"] = SHEAR_WIDE_X_SHA256
     if flow == "shear_twin":
         import shear_twin
         module = harness / "shear_twin.py"
@@ -128,14 +165,14 @@ def verify_source(repo, revision, flow="taylor_green", probes=False):
                 or hashlib.sha256(module.read_bytes()).hexdigest() != SHEAR_TWIN_SHA256):
             raise RuntimeError("Wrong shear_twin module checkout or content")
         result["shear_twin_sha256"] = SHEAR_TWIN_SHA256
-    if flow in ("shear_nested", "shear_twin"):  # the twin reuses the nested geometry
+    if flow in ("shear_nested", "shear_twin", "shear_wide_x", "shear_postg"):
         import shear_nested
         module = harness / "shear_nested.py"
         if (Path(shear_nested.__file__).resolve() != module
                 or hashlib.sha256(module.read_bytes()).hexdigest() != SHEAR_NESTED_SHA256):
             raise RuntimeError("Wrong shear_nested module checkout or content")
         result["shear_nested_sha256"] = SHEAR_NESTED_SHA256
-    if flow in ("shear_wave", "shear_nested", "shear_twin"):
+    if flow in ("shear_wave", "shear_nested", "shear_twin", "shear_wide_x", "shear_postg", "shear_axis_square", "shear_centered"):
         import shear_wave
         module = harness / "shear_wave.py"
         if (Path(shear_wave.__file__).resolve() != module
@@ -209,48 +246,87 @@ def main():
     p.add_argument("--output", type=Path, required=True)
     p.add_argument("--explosion", choices=["constant", "projected_linear"], default="constant")
     p.add_argument("--diagnostic",
-                   choices=[SCALED_DIAGNOSTIC, SCALED_SHEAR, NESTED_SHEAR, TWIN_SHEAR],
+                   choices=[SCALED_DIAGNOSTIC, SCALED_SHEAR, NESTED_SHEAR, TWIN_SHEAR,
+                            WIDE_X_SHEAR, POSTG_SHEAR, AXIS_SQUARE, CENTERED_SQUARE],
                    help="explicit four-case N256/fine512 same-placement diagnostic "
                         "(scaled_n256: Taylor-Green; scaled_shear_n256: shear wave) or the "
                         "fixed three-level/two-region shear preset (nested_shear_n128) "
-                        "or its depth-2 twin (twin_shear_n128)")
+                        "or its depth-2 twin (twin_shear_n128), or the single projected "
+                        "x-gap-width contrast (wide_x_shear_n128), or three-case POST-G "
+                        "initialization counterfactual (postg_shear_n128), or short square "
+                        "orientation diagnostic (axis_square_shear_n128)")
+    p.add_argument("--test-reconstruction", choices=["current", "centered_face_v1"],
+                   help="required only for centered_face_square_n128; private test build")
     p.add_argument("--probe-trace", action="store_true",
                    help="signed interface-probe trace (nested U128/U512/projected d3, "
-                        "twin projected d2); mechanism diagnostic, no gates")
+                        "twin projected d2, wide-x projected d3, POST-G three-case set, "
+                        "axis-square four-case set); "
+                        "mechanism diagnostic, no gates")
     p.add_argument("--final-snapshot", action="store_true",
                    help="nested_shear_n128 projected depth-3 only: final macro snapshot")
     p.add_argument("--early-pressure", action="store_true",
                    help="also save 65 scalar pressure samples at t*=k*.001875, k=0..64")
     p.add_argument("--flow", choices=["taylor_green", "shear_wave"], default="taylor_green")
     args = p.parse_args()
+    centered = args.diagnostic == CENTERED_SQUARE
+    if centered != (args.test_reconstruction is not None):
+        p.error("--test-reconstruction is required only for centered_face_square_n128")
+    if centered:
+        import shear_centered
+        if (args.n, args.depth, args.test_reconstruction) not in shear_centered.CASES:
+            p.error("centered comparison requires one of its four fixed cases")
+        if args.explosion != ("constant" if args.depth == 1 else "projected_linear"):
+            p.error("centered comparison uses constant uniforms and projected refined cases")
+    axis_square = args.diagnostic in (AXIS_SQUARE, CENTERED_SQUARE)
     nested = args.diagnostic == NESTED_SHEAR
-    if args.flow == "shear_wave" and (args.re != 100 or (args.depth == 3 and not nested)
+    wide_x = args.diagnostic == WIDE_X_SHEAR
+    postg = args.diagnostic == POSTG_SHEAR
+    if (wide_x or postg or axis_square) and not args.probe_trace:
+        p.error(f"{args.diagnostic} requires --probe-trace")
+    if args.flow == "shear_wave" and (args.re != 100 or (args.depth == 3 and not (nested or wide_x or postg))
                                       or args.early_pressure):
         p.error("shear_wave requires --re 100, depth 1 or 2 and no --early-pressure")
     probe_cases = {(NESTED_SHEAR, 128, 1, "constant"), (NESTED_SHEAR, 512, 1, "constant"),
                    (NESTED_SHEAR, 128, 3, "projected_linear"),
-                   (TWIN_SHEAR, 128, 2, "projected_linear")}
+                   (TWIN_SHEAR, 128, 2, "projected_linear"),
+                   (WIDE_X_SHEAR, 128, 3, "projected_linear"),
+                   (POSTG_SHEAR, 128, 1, "constant"), (POSTG_SHEAR, 512, 1, "constant"),
+                   (POSTG_SHEAR, 128, 3, "projected_linear")}
+    if axis_square:
+        import shear_axis_square
+        probe_cases.update((args.diagnostic, n, depth, scheme)
+                           for _, n, depth, scheme in shear_axis_square.CASES)
     if args.probe_trace and (args.final_snapshot or (
             args.diagnostic, args.n, args.depth, args.explosion) not in probe_cases):
-        p.error("--probe-trace requires nested U128/U512/projected depth 3 or twin projected "
-                "depth 2, without --final-snapshot")
+        p.error("--probe-trace requires a fixed nested, twin, wide-x, POST-G or axis-square case, "
+                "without --final-snapshot")
     if args.final_snapshot and not (nested and args.depth == 3
                                     and args.explosion == "projected_linear"):
         p.error("--final-snapshot requires nested_shear_n128 projected depth 3")
     initial, metrics, limits, initialization = flow_functions(args.flow)
     if args.output.exists():
         p.error("Output already exists; choose a new path to retain prior evidence")
-    if args.explosion != "constant" and args.depth != 2 and not nested:
+    if args.explosion != "constant" and args.depth != 2 and not (nested or wide_x or postg):
         p.error("projected_linear requires depth 2")
     twin = args.diagnostic == TWIN_SHEAR
-    if nested or twin:
+    if axis_square:
+        preset = shear_axis_square
+        if args.flow != "shear_wave" or (args.re, args.n, args.depth, args.explosion) not in preset.CASES:
+            p.error("axis_square_shear_n128 requires one of its four Re100 shear cases")
+        initial, metrics, limits, initialization = preset.initial, preset.metrics, {}, preset.INITIALIZATION
+    elif nested or twin or wide_x or postg:
         import shear_nested
         import shear_twin
-        preset = shear_nested if nested else shear_twin
+        import shear_wide_x
+        preset = shear_wide_x if wide_x else shear_nested if nested else shear_twin
+        if postg:
+            import shear_postg
+            preset = shear_postg
+            initial, initialization = preset.initial, preset.INITIALIZATION
         if args.flow != "shear_wave":
             p.error(f"--diagnostic {args.diagnostic} requires --flow shear_wave")
         if (args.re, args.n, args.depth, args.explosion) not in preset.CASES:
-            p.error(f"{args.diagnostic} requires one of its four Re100 cases")
+            p.error(f"{args.diagnostic} requires one of its fixed Re100 cases")
     elif args.diagnostic:
         cases = {(256, 1, "constant"), (512, 1, "constant"),
                  (256, 2, "constant"), (256, 2, "projected_linear")}
@@ -270,9 +346,9 @@ def main():
         p.error("--run-gpu and exactly one explicit CUDA_VISIBLE_DEVICES are required")
     repo = args.repo.resolve()
     sys.path[:0] = [str(repo / "src"), str(repo)]
-    source = verify_source(repo, args.expected_revision,
-                           "shear_nested" if nested else "shear_twin" if twin else args.flow,
-                           probes=args.probe_trace)
+    source_flow = ("shear_centered" if centered else "shear_axis_square" if axis_square else "shear_postg" if postg else "shear_wide_x" if wide_x else "shear_nested" if nested
+                   else "shear_twin" if twin else args.flow)
+    source = verify_source(repo, args.expected_revision, source_flow, probes=args.probe_trace)
     runtime = runtime_source(repo)
     if args.dry_run:
         print(json.dumps(dict(source=source, runtime_source=runtime, dry_run=True,
@@ -288,7 +364,7 @@ def main():
     cfg = config(args.n, args.re)
     if args.depth == 1:
         refinements = []
-    elif nested or twin:
+    elif nested or twin or wide_x or postg or axis_square:
         refinements = [list(level) for level in preset.REFINEMENTS]
     else:
         refinements = rectangles(args.n, args.depth)
@@ -296,6 +372,8 @@ def main():
     np.testing.assert_array_equal(leaf_partition(h), np.ones((args.n, args.n)))
     masks = [k == ACTIVE for k in h.kinds]
     main_steps, early_steps = sample_steps(args.n, args.early_pressure)
+    if axis_square:
+        main_steps = [0, preset.case_steps(args.n)[-1]]
     report = dict(
         schema=1,
         source=source,
@@ -323,7 +401,7 @@ def main():
         cuda_visible_devices=selected,
         device_id=0,
         limits=limits,
-        workload=workload(masks, round(TIMES[-1] * args.n / U)),
+        workload=workload(masks, main_steps[-1]),
         rows=[],
         early_pressure=dict(enabled=args.early_pressure, rows=[],
                             sampled_max_pressure_l2=None),
@@ -331,11 +409,23 @@ def main():
     )
     if args.flow != "taylor_green":
         import shear_wave
-        report.update(flow=args.flow, mode=list(shear_wave.MODE),
-                      polarization=shear_wave.E.tolist(), angle_deg=shear_wave.ANGLE_DEG,
-                      phase=shear_wave.PHI, eps=shear_wave.eps(args.n))
+        wave = preset if axis_square else shear_wave
+        report.update(flow=args.flow, mode=list(wave.MODE),
+                      polarization=wave.E.tolist(), angle_deg=wave.ANGLE_DEG,
+                      phase=wave.PHI, eps=wave.eps(args.n))
     if args.diagnostic:
         report["diagnostic"] = args.diagnostic
+    if postg:
+        report.update(initialization_id=preset.INITIALIZATION_ID, qualification_eligible=False,
+                      screening_scope="descriptive scalar screens only; no qualification assessor")
+    if axis_square:
+        report.update(initialization_id=preset.INITIALIZATION_ID, qualification_eligible=False,
+                      reference_harness=None, historical_baseline_revision=None,
+                      screening_scope="no numerical gates; short orientation diagnostic")
+    if centered:
+        report.update(test_reconstruction=args.test_reconstruction, observed_shifted_max=0.,
+                      reconstruction=dict(variant="current"),
+                      screening_scope="predeclared centered-witness diagnostic; separate four-case comparison")
     args.output.parent.mkdir(parents=True, exist_ok=True)
     def save():
         temporary = args.output.with_suffix(args.output.suffix + ".tmp")
@@ -350,18 +440,19 @@ def main():
     with CudaContext(0) as ctx, build_capture(report, save):
         report["runtime_identity"] = runtime_identity(ctx)
         save()
-        if args.depth == 1:
-            layout = BlockLayout2D.build(args.n, args.n, 8, 8, periodic=(True, True))
-            run = PeriodicBlockStepper(cfg, layout, context=ctx)
-        else:
-            run = FixedPatchStepper(
-                cfg,
-                context=ctx,
-                coarse_boundary="periodic",
-                refinements=refinements,
-                level_blocks=[(8, 8)] * args.depth,
-                explosion=args.explosion,
-            )
+        with (shear_centered.private_build(h, report, args.test_reconstruction) if centered else nullcontext()):
+            if args.depth == 1:
+                layout = BlockLayout2D.build(args.n, args.n, 8, 8, periodic=(True, True))
+                run = PeriodicBlockStepper(cfg, layout, context=ctx)
+            else:
+                run = FixedPatchStepper(
+                    cfg,
+                    context=ctx,
+                    coarse_boundary="periodic",
+                    refinements=refinements,
+                    level_blocks=[(8, 8)] * args.depth,
+                    explosion=args.explosion,
+                )
         with run:
             run.initialize(
                 *[
@@ -373,6 +464,8 @@ def main():
             probe_steps = set()
             if args.probe_trace:
                 import shear_probes
+                if axis_square:
+                    shear_probes = preset
                 probe_root = dict(zip(shear_probes.case_steps(args.n), shear_probes.STEPS))
                 probe_steps = set(probe_root)
                 probes = shear_probes.select(masks, args.n)
@@ -382,9 +475,18 @@ def main():
                 run.step(step - run.step_count)
                 if step in probe_steps:
                     trace = report["probe_trace"]
-                    trace["steps"].append(probe_root[step])  # root-physical step
-                    trace["values"].append(shear_probes.sample(
-                        (run.download(),) if args.depth == 1 else run.download(), probes))
+                    if "error" not in trace:
+                        arrays = (run.download(),) if args.depth == 1 else run.download()
+                        try:
+                            if centered:
+                                report["observed_shifted_max"] = max(report["observed_shifted_max"], shear_centered.observe(arrays))
+                            values = shear_probes.sample(arrays, probes)
+                            shear_probes.validate_row(values, len(probes))
+                        except (ValueError, ZeroDivisionError) as error:
+                            report["probe_trace"] = dict(error=str(error))
+                        else:
+                            trace["steps"].append(probe_root[step])  # root-physical step
+                            trace["values"].append(values)
                     if step not in main_steps and step not in early_steps:
                         continue  # probe-only step: no qualification metrics
                 row = sample(run, masks, step, args.n, cfg.viscosity, initial_mass,
@@ -439,7 +541,7 @@ def main():
                         ),
                         flush=True,
                     )
-    if args.probe_trace:
+    if args.probe_trace and "error" not in report["probe_trace"]:
         try:
             report["probe_trace"]["bytes"] = shear_probes.validate(report["probe_trace"])
         except ValueError as error:  # keep the scalar receipt; drop an invalid trace
@@ -447,7 +549,7 @@ def main():
     report.update(
         completed=True,
         seconds=time.monotonic() - started,
-        screening_pass=all(
+        screening_pass=None if axis_square else all(
             row[key] <= limit for row in report["rows"] for key, limit in limits.items()
         ),
     )
@@ -458,6 +560,8 @@ def main():
         ),
         flush=True,
     )
+    if axis_square:
+        return int("error" in report["probe_trace"])
     return 0 if report["screening_pass"] else 1
 
 

tests/validation/taylor_green/shear_axis_square.py

@@ -0,0 +1,136 @@
+"""Short axis-shear orientation diagnostic; no physical acceptance gates."""
+
+import json
+import numpy as np
+
+import shear_wave as sw
+from shear_probes import sample, validate_row, TRACE_MAX_BYTES
+from oracle import CX, CY, W, U, centers
+
+ROOT_N = 128
+MODE, E, PHI, ANGLE_DEG = (2, 0), np.array([0., 1.]), np.pi / 4, 0.
+REFINEMENTS = [[(32, 32, 96, 96)]]
+CASES = ((100, 128, 1, "constant"), (100, 256, 1, "constant"),
+         (100, 128, 2, "constant"), (100, 128, 2, "projected_linear"))
+STEPS = tuple(range(9))
+INITIALIZATION_ID = "postg_first_gradient_v1"
+INITIALIZATION = "physical FP64 axis-shear POST G=feq-(tau-1)A; first-gradient low-Mach CE"
+GROUPS = dict(inner_left=(32, 64), inner_right=(95, 63),
+              inner_bottom=(64, 32), inner_top=(63, 95),
+              outer_left=(31, 64), outer_right=(96, 63),
+              outer_bottom=(64, 31), outer_top=(63, 96))
+
+
+def eps(n):
+    return (4 * np.pi / n)**2
+
+
+def exact(x, t, length, viscosity):
+    """Transverse velocity and physical x derivative in the supplied length units."""
+    k = 4 * np.pi / length
+    amplitude = U * np.exp(-viscosity * k * k * t)
+    return amplitude * np.sin(k * x + PHI), amplitude * k * np.cos(k * x + PHI)
+
+
+def initial(shape, level, length, viscosity):
+    x, _ = centers(shape, level)
+    uy, gradient = exact(x, 0, length, viscosity)
+    tau = .5 + 3 * viscosity * 2**level
+    a = (3 * W * CX * CY)[:, None, None] * gradient * 2.**-level
+    return sw.equilibrium(np.ones_like(uy), np.zeros_like(uy), uy) - (tau - 1) * a
+
+
+def metrics(arrays, masks, t, length, viscosity, initial_mass=None):
+    rho, ux, uy, x, _, w, _ = sw._leaves(arrays, masks)
+    if not np.isclose(w.sum(), length**2, rtol=0, atol=1e-10):
+        raise ValueError("leaf partition does not cover domain once")
+    ey, _ = exact(x, t, length, viscosity)
+    p = (rho - 1) / 3
+    p -= w @ p / w.sum()  # one domain gauge, never the probe mean
+    mass = float(w @ rho)
+    return dict(tstar=t * U / length, mass=mass,
+                mass_relative=0. if initial_mass is None else abs(mass-initial_mass)/initial_mass,
+                velocity_l2=float(np.sqrt(w @ (ux**2 + (uy-ey)**2) / (w @ ey**2))),
+                spurious_pressure_rms=float(np.sqrt(w @ p**2 / w.sum())),
+                rho_min=float(rho.min()), rho_max=float(rho.max()),
+                max_speed=float(np.hypot(ux, uy).max()), leaf_cells=len(rho))
+
+
+def case_steps(n):
+    return [r * (n // ROOT_N) for r in STEPS]
+
+
+def select(masks, n=ROOT_N):
+    """Exact native cells filling eight unit-area groups; no interpolated samples."""
+    if (n, len(masks)) not in ((128, 1), (256, 1), (128, 2)):
+        raise ValueError("unsupported square probe layout")
+    probes = []
+    for name, (x, y) in GROUPS.items():
+        level = int(len(masks) == 2 and name.startswith("inner"))
+        factor = (n // ROOT_N) * 2**level
+        for dy in range(factor):
+            for dx in range(factor):
+                i, j = factor*x + dx, factor*y + dy
+                if not masks[level][j, i]:
+                    raise ValueError("square probe is not ACTIVE")
+                probes.append(dict(group=name, level=level, i=i, j=j,
+                                   center=[(i+.5)/factor, (j+.5)/factor],
+                                   width=1/factor, area=1/factor**2, parity=[i % 2, j % 2]))
+    if len({(p["level"], p["i"], p["j"]) for p in probes}) != len(probes):
+        raise ValueError("duplicate square probes")
+    return probes
+
+
+def new_trace(probes, case_viscosity, case):
+    return dict(schema=1, kind="axis-square orientation diagnostic; no pass gates", case=case,
+                time_unit="root step (N=128)", sample_phase="postcollision G at root synchronization",
+                root_n=ROOT_N, viscosity_root=U*ROOT_N/100, case_viscosity=case_viscosity,
+                mode=list(MODE), phase=PHI, amplitude=U, fields=["rho", "ux", "uy"],
+                probes=probes, steps=[], values=[])
+
+
+def group_moments(trace):
+    """Signed unit-area means; normalize by each group's mean actual-center derivative.
+
+    Momentum is rho*u. Raw density is ungauged. Per-cell parity stays in the trace.
+    U128 samples one center, whereas U256/refined-inner use four different centers.
+    """
+    out = []
+    for step, row in zip(trace["steps"], trace["values"]):
+        values, groups = np.asarray(row).reshape(-1, 3), {}
+        for name in GROUPS:
+            indices = [i for i, p in enumerate(trace["probes"]) if p["group"] == name]
+            probes = [trace["probes"][i] for i in indices]
+            weights = np.array([p["area"] for p in probes])
+            if weights.sum() != 1:
+                raise ValueError("probe group does not have unit area")
+            rho, ux, uy = values[indices].T
+            ey, gradient = exact(np.array([p["center"][0] for p in probes]), step,
+                                 ROOT_N, trace["viscosity_root"])
+            drho, djx, djy = (float(weights @ v) for v in (rho-1, rho*ux, rho*uy-ey))
+            derivative = float(weights @ gradient)
+            groups[name] = dict(drho=drho, djx=djx, djy=djy, mean_duy_dx=derivative,
+                                djx_over_duy_dx=djx/derivative, djy_over_duy_dx=djy/derivative)
+        out.append(dict(root_step=step, groups=groups))
+    return out
+
+
+def validate(trace):
+    if trace["steps"] != list(STEPS) or len(trace["values"]) != len(STEPS):
+        raise ValueError("square trace schedule mismatch")
+    _, n, depth, scheme = trace["case"]
+    if (100, n, depth, scheme) not in CASES or len(trace["probes"]) != (
+            20 if depth == 2 else 8*(n//ROOT_N)**2):
+        raise ValueError("square trace case/count mismatch")
+    for row in trace["values"]:
+        validate_row(row, len(trace["probes"]))
+    trace["group_moments"] = group_moments(trace)
+    size = 0
+    while True:
+        final = len(json.dumps(dict(trace, bytes=size), allow_nan=False).encode())
+        if final == size:
+            break
+        size = final
+    if size > TRACE_MAX_BYTES:
+        raise ValueError("square trace exceeds byte cap")
+    return size

tests/validation/taylor_green/shear_centered.py

@@ -0,0 +1,229 @@
+"""Test-only straight-face centered witness; private CUDA snapshot, frozen source intact."""
+from contextlib import contextmanager
+import hashlib
+import json
+from pathlib import Path
+import numpy as np
+
+import shear_axis_square as axis
+from oracle import CX, CY, W
+
+DIAGNOSTIC = "centered_face_square_n128"
+VARIANT = "centered_face_v1"
+BUILD_SHA256 = "81e880017495190209bd867d4b523a2c243be496ed4681892e4cd789b7ba1f21"
+CASES = ((128, 1, "current"), (256, 1, "current"),
+         (128, 2, "current"), (128, 2, VARIANT))
+# Local arithmetic-derived diagnostic deadbands, NOT accumulated trajectory bounds.
+POP_ENVELOPE, OP_ATOL, OP_RTOL = 1/32, 6e-8, 2e-7
+MARGIN = 2 * np.array([9, 6, 6]) * (OP_ATOL + OP_RTOL * POP_ENVELOPE)
+
+
+def routing(hierarchy):
+    """Square-only mask: full donor quartets, one missing normal axis, legal tangent."""
+    from CelerisLab.lbm.blocks import BlockLayout2D
+    from CelerisLab.lbm.refinement_tables import projected_explosion_tables
+    if hierarchy.depth != 2 or hierarchy.kinds[0].shape != (128, 128):
+        raise ValueError("centered witness requires the fixed two-level square")
+    pair = hierarchy.pairs[0]
+    covered = np.zeros((128, 128), bool); covered[32:96, 32:96] = True
+    if pair.fine_origin != (0., 0.) or not np.array_equal(pair.covered, covered):
+        raise ValueError("wrong centered-witness geometry")
+    layout = BlockLayout2D.build(128, 128, 8, 8, periodic=(True, True))
+    neighbors, offsets = projected_explosion_tables(pair, layout, hierarchy.kinds[0])
+    donors = layout.owner.ravel()[pair.ghost_coarse]
+    off = neighbors.reshape(-1, 2, 2) == donors[:, None, None]
+    disabled = off.all(axis=2)
+    y, x = np.divmod(pair.ghost_coarse, 128)
+    vertical = np.isin(x, [30, 31, 96, 97]) & (y >= 36) & (y < 92)
+    horizontal = np.isin(y, [30, 31, 96, 97]) & (x >= 36) & (x < 92)
+    mask = ((vertical & disabled[:, 0] & ~off[:, 1].any(axis=1)) |
+            (horizontal & disabled[:, 1] & ~off[:, 0].any(axis=1))).astype(np.uint8)
+    for donor in np.unique(pair.ghost_coarse):
+        rows = pair.ghost_coarse == donor
+        if rows.sum() != 4 or len(set(mask[rows])) != 1 or sorted(map(tuple, offsets[rows])) != [(-.25,-.25), (-.25,.25), (.25,-.25), (.25,.25)]:
+            raise ValueError("routing must preserve complete +/- quarter donor quartets")
+    return mask
+
+
+def metadata(mask):
+    return dict(variant=VARIANT, mask_sha256=hashlib.sha256(mask.tobytes()).hexdigest(),
+                ghost_rows=len(mask), selected_rows=int(mask.sum()),
+                rule="diagonal projected correction x2; audited straight faces only")
+
+
+def patch_source(source, mask):
+    """Change only the projected function in the private generated build directory."""
+    anchor = 'extern "C" __global__ void BlockExplodeProjected('
+    position = "        to[index_f(t % TOTAL_CELLS, q)] ="
+    before, kernel = source.split(anchor)
+    body, after = kernel.split("// Targets are receiver-local SOURCE slots", 1)
+    if body.count(position) != 1:
+        raise ValueError("unexpected projected kernel source")
+    flag = "__device__ __constant__ unsigned char centered_face_mask[] = {" + ",".join(map(str, mask)) + "};\n"
+    edit = "        if (centered_face_mask[row] && ex != 0.f && ey != 0.f) { gx *= 2.f; gy *= 2.f; }\n"
+    return before + flag + anchor + body.replace(position, edit + position) + "// Targets are receiver-local SOURCE slots" + after
+
+
+def source_digest(build):
+    digest = hashlib.sha256()
+    for path in sorted(Path(build).rglob("*")):
+        if path.is_file() and path.name != ".celeris-build":
+            digest.update(str(path.relative_to(build)).encode()); digest.update(path.read_bytes())
+    return digest.hexdigest()
+
+
+@contextmanager
+def private_build(hierarchy, report, variant):
+    """No kernel upload/launch monkey-patch; same production function and arguments."""
+    if variant == "current":
+        yield
+        return
+    from CelerisLab.cuda import compiler_v2 as compiler
+    mask = routing(hierarchy)
+    report["reconstruction"] = metadata(mask)
+    original, calls = compiler.generate_config, []
+    def generate(cfg, *args, **kwargs):
+        build = Path(original(cfg, *args, **kwargs))
+        calls.append(cfg.viscosity)
+        if len(calls) == 2:
+            if calls != [.0384, .0768]:
+                raise ValueError("unexpected candidate level configuration")
+            path = build / "blocks.cu"
+            path.write_text(patch_source(path.read_text(), mask))
+            sidecar = build / ".celeris-build"
+            data = json.loads(sidecar.read_text())
+            data.update(source_sha256=source_digest(build), experimental_reconstruction=metadata(mask))
+            sidecar.write_text(json.dumps(data, sort_keys=True))
+        elif len(calls) > 2:
+            raise ValueError("candidate is one adjacent pair only")
+        return str(build)
+    compiler.generate_config = generate
+    try:
+        yield
+        if calls != [.0384, .0768]:
+            raise ValueError("candidate did not build its exact two levels")
+    finally:
+        compiler.generate_config = original
+
+
+def observe(arrays):
+    """Bound observed synchronized ACTIVE values only, not intermediate GPU states."""
+    maximum = max(float(np.nanmax(abs(a - W[:, None, None]))) for a in arrays)
+    if not np.isfinite(maximum) or maximum > POP_ENVELOPE:
+        raise ValueError("observed shifted-population envelope exceeded")
+    return maximum
+
+
+def compare(records, expected_source):
+    """Four fresh receipts; diagnostic screens only, never a physical qualification."""
+    import run_case as rc
+    from assess_projected_pairs import BUILD_SOURCES, NUMERICS, valid_builds, finite_tree
+    from CelerisLab.lbm.hierarchy import build_hierarchy
+    from CelerisLab.lbm.refinement import ACTIVE
+    if not finite_tree(records) or not expected_source.get("revision"):
+        raise ValueError("nonfinite receipt or missing exact source identity")
+    keys = [(r['n'], r['depth'], r.get('test_reconstruction')) for r in records]
+    if keys != list(CASES):
+        raise ValueError("expected ordered fresh U128/U256/current/candidate cases")
+    for r in records:
+        for key in ('source', 'runtime_identity', 'runtime_source', 'cuda_visible_devices'):
+            if not r.get(key) or r[key] != records[0][key]:
+                raise ValueError(f"mixed or missing {key}")
+        if (r['source'] != expected_source or r.get('numerical_revision') != rc.NUMERICAL_REVISION
+                or r.get('numerics') != NUMERICS or r.get('re_L') != 100 or r.get('U') != .03
+                or r.get('polarization') != [0., 1.] or r.get('initialization_id') != axis.INITIALIZATION_ID
+                or r['source'].get('dirty') is not False or r['source'].get('source_tree_sha256') != rc.SOURCE_SHA256
+                or r.get('diagnostic') != DIAGNOSTIC or not r.get('completed')
+                or r.get('qualification_eligible') is not False or r.get('initialization') != axis.INITIALIZATION
+                or r.get('mode') != [2, 0] or r.get('phase') != axis.PHI
+                or r.get('explosion') != ('constant' if r['depth'] == 1 else 'projected_linear')
+                or not 0 <= r.get('observed_shifted_max', np.inf) <= POP_ENVELOPE):
+            raise ValueError("wrong candidate receipt identity or observed envelope")
+        h = build_hierarchy(rc.config(r['n'], 100), [] if r['depth'] == 1 else axis.REFINEMENTS)
+        if r.get('refinements') != ([] if r['depth'] == 1 else [[list(axis.REFINEMENTS[0][0])]]):
+            raise ValueError("wrong geometry")
+        if (r.get('viscosity') != .03*r['n']/100
+                or r.get('tau_by_level') != [.5+3*c.viscosity for c in h.configs]
+                or r.get('workload') != rc.workload([k == ACTIVE for k in h.kinds], 8*(r['n']//128))):
+            raise ValueError("wrong physical scales or workload")
+        builds = r.get('actual_build_sidecars', [])
+        canonical = []
+        for level, build in enumerate(builds):
+            modified = r['test_reconstruction'] == VARIANT and level == 1
+            expected = BUILD_SHA256 if modified else BUILD_SOURCES[r['n']*2**level]
+            if (build.get('source_sha256') != expected
+                    or build.get('experimental_reconstruction') != (metadata(routing(h)) if modified else None)):
+                raise ValueError("wrong actual build or reconstruction identity")
+            # Reuse the frozen build-metadata checker after checking the private-source hash above.
+            canonical.append(dict(build, source_sha256=BUILD_SOURCES[r['n']*2**level]))
+        if not valid_builds(dict(actual_build_sidecars=canonical), r['n'], r['depth']):
+            raise ValueError("invalid architecture, flags, config or PTX digest")
+        if r.get('reconstruction') != (metadata(routing(h)) if r['test_reconstruction'] == VARIANT else dict(variant='current')):
+            raise ValueError("wrong receipt mask identity")
+        trace = r['probe_trace']
+        declared = axis.new_trace(trace['probes'], r['viscosity'], [DIAGNOSTIC, r['n'], r['depth'], r['explosion']])
+        if any(trace.get(k) != v for k,v in declared.items() if k not in ('steps','values')):
+            raise ValueError("wrong trace phase, clock or physical identity")
+        if trace['probes'] != axis.select([k == ACTIVE for k in h.kinds], r['n']):
+            raise ValueError("wrong actual probe cells")
+        if trace.get("bytes") != axis.validate(trace):
+            raise ValueError("trace serialized byte count mismatch")
+        if np.min(np.asarray(trace["values"])[:,0::3]) <= 0:
+            raise ValueError("nonpositive sampled density")
+        if [row['step'] for row in r['rows']] != [0, 8*(r['n']//128)] or not all(rc.finite_row(v) for v in r['rows']):
+            raise ValueError("invalid scalar rows")
+        if min(row['rho_min'] for row in r['rows']) <= 0:
+            raise ValueError("nonpositive density")
+    current, candidate = records[2:]
+    if current['probe_trace']['values'][0] != candidate['probe_trace']['values'][0]:
+        raise ValueError("candidate/current initial samples differ; cancellation invalid")
+    vectors = []
+    for r in records:
+        groups = axis.group_moments(r['probe_trace'])
+        vectors.append(np.array([[[row['groups'][name][key] for key in ('drho','djx','djy')]
+                                  for name in axis.GROUPS] for row in groups]))
+    # Inner quadrature matches U256; outer quadrature matches U128.
+    uniform = vectors[1].copy(); uniform[:, 4:] = vectors[0][:, 4:]
+    residual = [v-v[0]-uniform+uniform[0] for v in vectors[2:]]
+    old, new = residual
+    checks = dict(bottom_step1_clear=abs(old[1,2,1])-abs(new[1,2,1]) > MARGIN[1],
+                  top_step1_nonworse=abs(new[1,3,1]) <= abs(old[1,3,1])+MARGIN[1])
+    for face, indices in [('bottom',[2,6]), ('top',[3,7])]:
+        rms = lambda a: float(np.sqrt(np.mean(a[8,indices,1]**2)))
+        checks[face+'_step8_nonworse'] = rms(new) <= rms(old)+MARGIN[1]
+    checks['vertical_null'] = bool(np.all(abs(new[:,[0,1,4,5]]-old[:,[0,1,4,5]]) <= MARGIN))
+    acoustic = lambda a: np.sqrt(a[...,0]**2/3 + a[...,1]**2 + a[...,2]**2)
+    acoustic_margin = float(np.sqrt(MARGIN[0]**2/3 + MARGIN[1]**2 + MARGIN[2]**2))
+    checks['acoustic_nonworse'] = bool(np.all(acoustic(new) <= acoustic(old)+acoustic_margin))
+    raw_pressure = lambda r: float(np.max(abs(np.asarray(r['probe_trace']['values'])[:,0::3]-1))/3)
+    checks['sampled_pressure_nonworse'] = raw_pressure(candidate) <= raw_pressure(current)+MARGIN[0]/3
+    checks['step8_pressure_nonworse'] = candidate['rows'][-1]['spurious_pressure_rms'] <= current['rows'][-1]['spurious_pressure_rms']+MARGIN[0]/3
+    return dict(qualified=False, advance_to_review=all(checks.values()), checks=checks,
+                moment_deadband=MARGIN.tolist(), deadband_scope="local arithmetic-derived diagnostic; no accumulated error guarantee",
+                groups=list(axis.GROUPS), components=['drho','djx','djy'],
+                current_signed_residual=old.tolist(), candidate_signed_residual=new.tolist(),
+                raw_current=vectors[2].tolist(), raw_candidate=vectors[3].tolist(),
+                group_moments=[axis.group_moments(r['probe_trace']) for r in records],
+                significant_sign_flips=((old*new < 0) & (abs(old)>MARGIN) & (abs(new)>MARGIN)).tolist(),
+                top_step1_below_deadband=bool(abs(old[1,3,1]) <= MARGIN[1]),
+                note="Absolute reduction with sign retained; top below deadband may be inconclusive.")
+
+
+def main():
+    import argparse
+    p = argparse.ArgumentParser(description=__doc__)
+    p.add_argument("--repo", required=True, type=Path)
+    p.add_argument("--expected-revision", required=True)
+    p.add_argument("receipts", nargs=4, type=Path, help="U128 U256 current candidate, in order")
+    p.add_argument("--output", required=True, type=Path)
+    args = p.parse_args()
+    import run_case as rc
+    source = rc.verify_source(args.repo.resolve(), args.expected_revision, "shear_centered", probes=True)
+    result = compare([json.loads(path.read_text()) for path in args.receipts], source)
+    with args.output.open("x") as stream:
+        json.dump(result, stream, indent=2, allow_nan=False)
+    return 0 if result["advance_to_review"] else 1
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())

tests/validation/taylor_green/shear_postg.py

@@ -0,0 +1,34 @@
+"""First-gradient POST-G initialization counterfactual, not a qualification preset.
+
+PRE F = feq - tau*A; POST G = feq - (tau - 1)*A. This low-Mach,
+first-order Chapman-Enskog approximation is not an exact discrete eigenstate.
+The physical viscous-stress factor tau - 1/2 is not the raw POST-G factor.
+"""
+
+import numpy as np
+
+import shear_nested as nested
+import shear_wave as sw
+from oracle import CX, CY, W, U, centers
+
+REFINEMENTS = nested.REFINEMENTS
+CASES = ((100, 128, 1, "constant"), (100, 512, 1, "constant"),
+         (100, 128, 3, "projected_linear"))
+INITIALIZATION_ID = "postg_first_gradient_v1"
+INITIALIZATION = (
+    "physical FP64 POST G=feq-(tau-1)A at actual cell centers; "
+    "A=3*w*U*cos(theta)*(c.e)*(c.k_level); rho=1; first-gradient low-Mach CE"
+)
+
+
+def initial(shape, level, length, viscosity):
+    """Return physical FP64 G; production initialize owns shifted FP32 encoding."""
+    x, y = centers(shape, level)
+    k = sw.wavevector(length)
+    theta = k[0] * x + k[1] * y + sw.PHI
+    k_level = k * 2.0**-level
+    tau = 0.5 + 3 * viscosity * 2**level
+    ce = CX * sw.E[0] + CY * sw.E[1]
+    ck = CX * k_level[0] + CY * k_level[1]
+    correction = (3 * W * ce * ck)[:, None, None] * U * np.cos(theta)
+    return sw.initial(shape, level, length, viscosity) - (tau - 1) * correction

tests/validation/taylor_green/shear_probes.py

@@ -107,6 +107,12 @@ def residuals(trace):
     return out
 
 
+def validate_row(row, probe_count):
+    """Reject an invalid optional sample before it can enter a scalar receipt."""
+    if len(row) != 3 * probe_count or not all(math.isfinite(v) for v in row):
+        raise ValueError("malformed or nonfinite trace row")
+
+
 def validate(trace):
     """Exact schedule, finite values, fixed probe count and byte cap."""
     import json
@@ -114,10 +120,15 @@ def validate(trace):
         raise ValueError("trace steps differ from the fixed schedule")
     if len(trace["probes"]) != len(NORMAL_COORDS) * len(TRANSECTS):
         raise ValueError("unexpected probe count")
-    width = 3 * len(trace["probes"])
-    if any(len(r) != width or not all(math.isfinite(v) for v in r) for r in trace["values"]):
-        raise ValueError("malformed or nonfinite trace row")
-    size = len(json.dumps(trace, allow_nan=False).encode())
+    for row in trace["values"]:
+        validate_row(row, len(trace["probes"]))
+    # Include the final byte-count field itself, including its number of digits.
+    size = 0
+    while True:
+        final_size = len(json.dumps(dict(trace, bytes=size), allow_nan=False).encode())
+        if final_size == size:
+            break
+        size = final_size
     if size > TRACE_MAX_BYTES:
         raise ValueError(f"trace exceeds {TRACE_MAX_BYTES} bytes ({size})")
     return size

tests/validation/taylor_green/shear_wide_x.py

@@ -0,0 +1,20 @@
+"""One projected x-gap-width diagnostic (opt-in wide_x_shear_n128).
+
+Keep the nested outer boxes and all three spacings. Shrink only L2 in x,
+so each x-normal middle gap grows from 6 to 12 root cells; y gaps stay 6.
+This also reduces L2 volume and moves its x faces. It is a descriptive
+geometry contrast with no new qualification matrix, assessor or pass gate.
+"""
+
+import shear_nested as nested
+
+ROOT_N = nested.ROOT_N
+DEPTH = nested.DEPTH
+REFINEMENTS = [nested.REFINEMENTS[0],
+               [(56, 44, 96, 212), (160, 60, 200, 228)]]  # parent L1 coordinates
+CASES = ((100, ROOT_N, DEPTH, "projected_linear"),)
+LEAVES = [7936, 20352, 53760]
+HOPS = nested.HOPS
+WORK = sum(c * h for c, h in zip(LEAVES, HOPS))  # 540016640 active updates
+NU = nested.NU
+TAU = nested.TAU

tests/validation/taylor_green/test_axis_square.py

@@ -0,0 +1,172 @@
+"""Fixed square geometry, native sampling, POST phase and runner admission; CPU only."""
+import hashlib
+import io
+import json
+import sys
+import tempfile
+import unittest
+from copy import deepcopy
+from pathlib import Path
+from unittest.mock import patch
+from contextlib import redirect_stdout
+
+import numpy as np
+import run_case as rc
+import shear_axis_square as axis
+import shear_wave as sw
+import test_shear_wave as legacy
+from oracle import CX, CY, W, U, centers
+from CelerisLab.lbm.hierarchy import build_hierarchy, leaf_partition
+from CelerisLab.lbm.refinement import ACTIVE
+
+REPO = Path(__file__).resolve().parents[3]
+
+
+def masks(n, depth):
+    h = build_hierarchy(rc.config(n, 100), [] if depth == 1 else axis.REFINEMENTS)
+    np.testing.assert_array_equal(leaf_partition(h), np.ones((n, n)))
+    return [k == ACTIVE for k in h.kinds]
+
+
+class AxisSquareTests(unittest.TestCase):
+    def test_geometry_native_groups_and_work(self):
+        work = 0
+        for _, n, depth, _ in axis.CASES:
+            m = masks(n, depth)
+            p = axis.select(m, n)
+            self.assertEqual(len(p), 20 if depth == 2 else 8*(n//128)**2)
+            self.assertEqual(len({(v['level'], v['i'], v['j']) for v in p}), len(p))
+            roots = [(32, 64), (95, 63), (64, 32), (63, 95),
+                     (31, 64), (96, 63), (64, 31), (63, 96)]
+            expected = []
+            for group, (x, y) in enumerate(roots):
+                level = 1 if depth == 2 and group < 4 else 0
+                factor = (n//128) * 2**level
+                expected += [(level, factor*x+dx, factor*y+dy,
+                              [x+(dx+.5)/factor, y+(dy+.5)/factor])
+                             for dy in range(factor) for dx in range(factor)]
+            self.assertEqual([(v['level'], v['i'], v['j'], v['center']) for v in p], expected)
+            for name in axis.GROUPS:
+                group = [v for v in p if v['group'] == name]
+                self.assertEqual(sum(v['area'] for v in group), 1)
+                np.testing.assert_array_equal(np.mean([v['center'] for v in group], axis=0),
+                                              np.array(axis.GROUPS[name])+.5)
+                for v in group:
+                    self.assertTrue(m[v['level']][v['j'], v['i']])
+            w = rc.workload(m, axis.case_steps(n)[-1])
+            work += w['active_collision_updates']
+            if depth == 2:
+                self.assertEqual(w['active_cells_by_level'], [12288, 16384])
+                self.assertEqual(w['level_hops'], [8, 16])
+            self.assertEqual(axis.case_steps(n), [r*(n//128) for r in range(9)])
+        self.assertEqual(work, 1900544)
+
+    def test_postg_moments_stress_and_matching_fine_control(self):
+        for level in (0, 1):
+            x, _ = centers((33, 41), level)
+            k, h, tau = 4*np.pi/128, 2.**-level, .5+3*.0384*2**level
+            uy = U*np.sin(k*x+np.pi/4)
+            eq = sw.equilibrium(np.ones_like(x), np.zeros_like(x), uy)
+            g = axis.initial(x.shape, level, 128, .0384)
+            a = 3*W[:, None, None]*CX[:, None, None]*CY[:, None, None]*h*k*U*np.cos(k*x+np.pi/4)
+            np.testing.assert_allclose(g, eq-(tau-1)*a, rtol=0, atol=6e-17)
+            for v in (np.ones(9), CX, CY):
+                np.testing.assert_allclose(np.einsum('q,qyx->yx', v, g-eq), 0, rtol=0, atol=2e-16)
+            np.testing.assert_allclose(np.einsum('q,qyx->yx', CX*CY, g-eq),
+                                       (1-tau)*h*k*U*np.cos(k*x+np.pi/4)/3, atol=6e-17)
+            self.assertGreater(g.min(), 0)
+        np.testing.assert_array_equal(axis.initial((256, 256), 0, 256, .0768),
+                                      axis.initial((256, 256), 1, 128, .0384))
+
+    def test_signed_momentum_and_actual_center_time(self):
+        m = masks(128, 2)
+        probes = axis.select(m)
+        trace = axis.new_trace(probes, .0384, [rc.AXIS_SQUARE, 128, 2, 'projected_linear'])
+        for step in axis.STEPS:
+            arrays = []
+            for level, mask in enumerate(m):
+                x, _ = centers(mask.shape, level)
+                uy, _ = axis.exact(x, step, 128, .0384)
+                arrays.append(sw.equilibrium(np.full_like(x, 1.001), np.zeros_like(x), uy))
+            trace['steps'].append(step)
+            trace['values'].append(axis.sample(arrays, probes))
+        trace['bytes'] = axis.validate(trace)
+        self.assertEqual(trace['bytes'], len(json.dumps(trace, allow_nan=False).encode()))
+        for row in trace['group_moments']:
+            for name, v in row['groups'].items():
+                self.assertAlmostEqual(v['drho'], .001, places=14)
+                x = np.array([p['center'][0] for p in probes if p['group'] == name])
+                uy, gradient = axis.exact(x, row['root_step'], 128, .0384)
+                self.assertAlmostEqual(v['djy'], .001*uy.mean(), places=14)
+                self.assertAlmostEqual(v['mean_duy_dx'], gradient.mean(), places=14)
+            self.assertAlmostEqual(row['groups']['inner_left']['mean_duy_dx'],
+                                   -row['groups']['inner_bottom']['mean_duy_dx'], places=14)
+        with patch.object(axis, 'TRACE_MAX_BYTES', trace['bytes']-1), self.assertRaises(ValueError):
+            axis.validate(trace)
+        for mutate in (lambda t: t['steps'].pop(),
+                       lambda t: t['steps'].__setitem__(1, 0),
+                       lambda t: t['values'][1].__setitem__(0, float('nan')),
+                       lambda t: t['probes'].pop()):
+            bad = deepcopy(trace); mutate(bad)
+            with self.assertRaises(ValueError): axis.validate(bad)
+
+    def test_four_full_mock_schedules_and_old_rejection(self):
+        helper = legacy.ProbeTraceTests()
+        helper.setUp()
+        def arrays(m, n, t):
+            return [axis.initial(v.shape, level, n, U*n/100) for level, v in enumerate(m)]
+        for _, n, depth, scheme in axis.CASES:
+            with patch.object(helper, 'exact_arrays', side_effect=arrays):
+                r, counts, files = helper.run_mocked(n, depth, rc.AXIS_SQUARE, scheme,
+                                                     masks(n, depth), initializer=axis.initial)
+            self.assertEqual(counts[-1], 8*(n//128))
+            self.assertEqual([v['step'] for v in r['rows']], [0, 8*(n//128)])
+            self.assertEqual(r['probe_trace']['steps'], list(range(9)))
+            self.assertEqual(r['probe_trace']['probes'], axis.select(masks(n, depth), n))
+            self.assertEqual(r['mode'], [2, 0])
+            self.assertEqual(r['phase'], np.pi/4)
+            self.assertEqual(r['initialization'], axis.INITIALIZATION)
+            self.assertEqual(r['limits'], {})
+            self.assertIsNone(r['screening_pass'])
+            self.assertFalse(r['qualification_eligible'])
+            self.assertEqual(files, ['probe.json'])
+        old = legacy.NestedPresetTests(); old.setUp()
+        records = old.nested_set(); records[-1]['diagnostic'] = rc.AXIS_SQUARE
+        self.assertEqual(old.mod.assess_nested(records, 'r')['status'], 'INVALID_SET')
+        records = [legacy.receipt(128, 1), legacy.receipt(256, 1), legacy.receipt(128, 2),
+                   legacy.receipt(128, 2, 'projected_linear')]
+        records[-1].update(diagnostic=rc.AXIS_SQUARE, initialization=axis.INITIALIZATION)
+        self.assertEqual(old.mod.assess(records, 'r')['status'], 'INVALID_SET')
+        self.assertNotIn(rc.AXIS_SQUARE, old.mod.PRESETS)
+        with patch.object(axis, 'sample', return_value=[float('nan')]*24):
+            r, _, _ = helper.run_mocked(128, 1, rc.AXIS_SQUARE, 'constant', masks(128, 1),
+                                        initializer=axis.initial)
+        self.assertIn('error', r['probe_trace'])
+        self.assertTrue(r['completed'])
+        self.assertEqual(len(r['rows']), 2)
+
+    def test_admission_and_module_pin(self):
+        from tests.validation.multigrid import check_projected_pairs as gate
+        self.assertEqual(hashlib.sha256(Path(axis.__file__).read_bytes()).hexdigest(), rc.SHEAR_AXIS_SQUARE_SHA256)
+        with tempfile.TemporaryDirectory() as tmp:
+            base = ['run_case.py', '--repo', str(REPO), '--expected-revision', 'f'*40,
+                    '--output', str(Path(tmp)/'case.json'), '--dry-run', '--re', '100',
+                    '--flow', 'shear_wave', '--diagnostic', rc.AXIS_SQUARE, '--probe-trace']
+            for n, depth, scheme, valid in [(n, d, s, True) for _, n, d, s in axis.CASES] + [
+                    (512, 1, 'constant', False), (256, 2, 'projected_linear', False), (128, 3, 'constant', False)]:
+                for extra in ([], ['--final-snapshot'], ['--early-pressure'], ['--flow', 'taylor_green'], ['--re', '1000']):
+                    argv = base + ['--n', str(n), '--depth', str(depth), '--explosion', scheme] + extra
+                    with patch.object(sys, 'argv', argv), redirect_stdout(io.StringIO()), patch('sys.stderr', io.StringIO()), \
+                         patch.object(rc, 'verify_source', return_value={}) as verify, patch.object(rc, 'runtime_source', return_value={}):
+                        if valid and not extra:
+                            self.assertEqual(rc.main(), 0)
+                            verify.assert_called_once_with(REPO, 'f'*40, 'shear_axis_square', probes=True)
+                        else:
+                            with self.assertRaises(SystemExit): rc.main()
+                            verify.assert_not_called()
+        with patch.object(gate, 'preflight', return_value=rc.NUMERICAL_REVISION), \
+             patch.object(rc, 'provenance', return_value=dict(source_tree_sha256=rc.SOURCE_SHA256)):
+            rc.verify_source(REPO, 'r', 'shear_axis_square', probes=True)
+            for module, key, value in ((axis, '__file__', '/wrong/helper.py'), (rc, 'SHEAR_AXIS_SQUARE_SHA256', 'wrong')):
+                with patch.object(module, key, value), self.assertRaisesRegex(RuntimeError, 'shear_axis_square'):
+                    rc.verify_source(REPO, 'r', 'shear_axis_square', probes=True)

tests/validation/taylor_green/test_centered.py

@@ -0,0 +1,246 @@
+"""Independent coordinate oracle, exact kernel host shim and closed pair budgets."""
+import ctypes
+from contextlib import nullcontext
+from copy import deepcopy
+from dataclasses import replace
+import hashlib
+import json
+from pathlib import Path
+import subprocess
+import tempfile
+import unittest
+from unittest.mock import patch
+
+import numpy as np
+import run_case as rc
+import shear_axis_square as axis
+import shear_centered as candidate
+import test_shear_wave as legacy
+from test_axis_square import masks
+from CelerisLab.lbm.hierarchy import build_hierarchy
+from CelerisLab.lbm.blocks import BlockLayout2D
+from CelerisLab.lbm.refinement_tables import projected_explosion_tables
+from CelerisLab.lbm.refinement import ACTIVE, OUTER_GHOST
+from tests.unit import test_projected_explosion as old
+
+REPO = Path(__file__).resolve().parents[3]
+SOURCE = REPO/'src/CelerisLab/lbm/kernels/blocks.cu'
+
+
+def oracle(c, pair, centered=True):
+    """Uses only dense coordinates/legal donors, never production neighbor/mask arrays."""
+    out = np.empty((9, len(pair.ghost_fine)))
+    legal = set(map(int, pair.ghost_coarse))
+    for r, (donor, target) in enumerate(zip(pair.ghost_coarse, pair.ghost_fine)):
+        y, x = divmod(int(donor), 128); fy, fx = divmod(int(target), 256)
+        delta = np.array([(fx+.5)/2-x-.5, (fy+.5)/2-y-.5])
+        usable = [(128*(y-dy)+x-dx in legal and 128*(y+dy)+x+dx in legal)
+                  for dx,dy in [(1,0),(0,1)]]
+        straight = ((x in (30,31,96,97) and 36 <= y < 92 and usable == [False,True])
+                    or (y in (30,31,96,97) and 36 <= x < 92 and usable == [True,False]))
+        for q, e in enumerate(old.E):
+            grad = np.array([(c[q,y+dy,x+dx]-c[q,y-dy,x-dx])/2 if ok else 0.
+                             for (dx,dy),ok in zip([(1,0),(0,1)],usable)])
+            if q: grad -= e*(grad@e)/(e@e)
+            out[q,r] = c[q,y,x] + (2 if centered and straight and np.all(e) else 1)*(grad@delta)
+    return out
+
+
+class CenteredTests(unittest.TestCase):
+    @classmethod
+    def setUpClass(cls):
+        cls.h = build_hierarchy(rc.config(128,100), axis.REFINEMENTS)
+        cls.p = cls.h.pairs[0]
+        cls.mask = candidate.routing(cls.h)
+        cls.layout = BlockLayout2D.build(128,128,8,8,periodic=(True,True))
+        cls.ids, cls.offsets = projected_explosion_tables(cls.p,cls.layout,cls.h.kinds[0])
+        cls.neighbors = old.decode(cls.layout,cls.ids)
+
+    def test_mask_constants_quartets_and_transpose(self):
+        self.assertEqual(candidate.metadata(self.mask), dict(variant=candidate.VARIANT,
+            mask_sha256='112de63ce9bafd2e205a9e7029ef6ce5683fe14407b6ac04444579c84960b7ea',
+            ghost_rows=2112, selected_rows=1792,
+            rule='diagonal projected correction x2; audited straight faces only'))
+        c=np.random.default_rng(47).normal(0,.001,(9,128,128))
+        actual=oracle(c,self.p); baseline=oracle(c,self.p,False)
+        base,g=old.reconstruct(c,self.p,self.neighbors,self.offsets)
+        g[np.ix_(np.flatnonzero(np.all(old.E,axis=1)),np.flatnonzero(self.mask),[0,1])]*=2
+        expected=c.reshape(9,-1)[:,self.p.ghost_coarse]+np.einsum('qra,ra->qr',g,self.offsets)
+        np.testing.assert_allclose(actual,expected,atol=1e-18,rtol=0)
+        np.testing.assert_array_equal(actual[:,~self.mask.astype(bool)],baseline[:,~self.mask.astype(bool)])
+        for d in np.unique(self.p.ghost_coarse):
+            rows=self.p.ghost_coarse==d
+            np.testing.assert_allclose(actual[:,rows].mean(axis=1),c.reshape(9,-1)[:,d],atol=1e-18)
+        const=np.broadcast_to(old.W[:,None,None],c.shape)
+        np.testing.assert_array_equal(oracle(const,self.p),np.broadcast_to(old.W[:,None],actual.shape))
+        qswap=[list(map(tuple,old.E)).index(tuple(e[::-1])) for e in old.E]
+        index={int(t):r for r,t in enumerate(self.p.ghost_fine)}
+        swap=[index[int(t%256)*256+int(t//256)] for t in self.p.ghost_fine]
+        np.testing.assert_array_equal(self.mask,self.mask[swap])
+        np.testing.assert_allclose(oracle(c[qswap].transpose(0,2,1),self.p),actual[qswap][:,swap],atol=1e-18)
+        qreflect=[list(map(tuple,old.E)).index(tuple((-e[0],e[1]))) for e in old.E]
+        reflect=[index[int(t//256)*256+255-int(t%256)] for t in self.p.ghost_fine]
+        np.testing.assert_array_equal(self.mask,self.mask[reflect])
+        np.testing.assert_allclose(oracle(c[qreflect,:,::-1],self.p),actual[qreflect][:,reflect],atol=1e-18)
+        # Central dependency corridors remain in a constant mask regime beyond radius8*2.5.
+        for name,(x,y) in axis.GROUPS.items():
+            tangent = y if name.endswith(('left','right')) else x
+            self.assertGreater(min(tangent+.25-36,92-(tangent+.75)),20)
+
+    def test_actual_generated_kernel_shifted_unshifted(self):
+        child=self.p.fine_layout((8,4))  # unequal owner strides, unused slots poisoned
+        patched=candidate.patch_source(SOURCE.read_text(),self.mask)
+        kernel=patched.split('__device__ __constant__ unsigned char centered_face_mask',1)[1]
+        kernel='unsigned char centered_face_mask'+kernel.split('// Targets are receiver-local SOURCE slots',1)[0]
+        kernel=kernel.replace('extern "C" __global__','extern "C"')
+        code='#include <cstdint>\nusing fpxx=float;\n'
+        code+=f'constexpr unsigned NQ=9,TOTAL_CELLS={np.prod(child.local_shape)};\n'
+        code+='struct Id{unsigned x;}; Id blockIdx{0},blockDim{1},threadIdx{0};\n'
+        code+='unsigned index_f(unsigned k,unsigned q){return q*TOTAL_CELLS+k;}\n'
+        code+='int d_cx[]={0,1,-1,0,0,1,-1,1,-1},d_cy[]={0,0,0,1,-1,1,-1,-1,1};\n'+kernel
+        original=SOURCE.read_text().split('extern "C" __global__ void BlockExplodeProjected(',1)[1].split('// Targets are receiver-local SOURCE slots',1)[0]
+        code+='void Original('+original
+        code+='extern "C" void invoke(float*c,float*f,uint64_t*t,uint64_t*d,unsigned n,unsigned long l,uint64_t*b,float*o){for(unsigned r=0;r<n;r++){threadIdx.x=r;BlockExplodeProjected(c,f,t,d,n,l,b,o);}}'
+        code+=code[code.index('extern "C" void invoke'):].replace('void invoke(', 'void invoke_old(').replace('BlockExplodeProjected(c,f,','Original(c,f,')
+        with tempfile.TemporaryDirectory() as tmp:
+            src=Path(tmp)/'kernel.cpp'; lib=Path(tmp)/'kernel.so'; src.write_text(code)
+            subprocess.run(['g++','-std=c++17','-O2','-shared','-fPIC',str(src),'-o',str(lib)],check=True,capture_output=True)
+            library=ctypes.CDLL(str(lib));fn=library.invoke;fn_old=library.invoke_old
+            ptr=lambda a:a.ctypes.data_as(ctypes.c_void_p)
+            targets=np.ascontiguousarray(child.owner.ravel()[self.p.ghost_fine],np.uint64)
+            donors=np.ascontiguousarray(self.layout.owner.ravel()[self.p.ghost_coarse],np.uint64)
+            for shifted in (True,False):
+                c=np.random.default_rng(8).normal(0,.001,(9,128,128))
+                if not shifted:c+=old.W[:,None,None]
+                c=c.astype(np.float32); c[:,self.p.covered]=np.nan
+                packed=np.full(self.layout.nblocks*9*np.prod(self.layout.local_shape),np.nan,np.float32)
+                for q in range(9):packed[self.layout.slot(self.layout.owner.ravel(),q)]=c[q].ravel()
+                fine=np.full(child.nblocks*9*np.prod(child.local_shape),np.nan,np.float32)
+                fn(ptr(packed),ptr(fine),ptr(targets),ptr(donors),ctypes.c_uint(len(targets)),ctypes.c_ulong(np.prod(self.layout.local_shape)),ptr(self.ids),ptr(self.offsets))
+                original=np.full_like(fine,np.nan)
+                fn_old(ptr(packed),ptr(original),ptr(targets),ptr(donors),ctypes.c_uint(len(targets)),ctypes.c_ulong(np.prod(self.layout.local_shape)),ptr(self.ids),ptr(self.offsets))
+                output=np.array([fine[child.slot(targets,q)] for q in range(9)])
+                unchanged=np.array([original[child.slot(targets,q)] for q in range(9)])
+                np.testing.assert_array_equal(output[:,~self.mask.astype(bool)],unchanged[:,~self.mask.astype(bool)])
+                np.testing.assert_array_equal(output[:5],unchanged[:5])
+                np.testing.assert_allclose(output,oracle(c.astype(float),self.p),atol=6e-8,rtol=2e-7)
+
+    def test_mixed_mask_full_pair_mass_momentum(self):
+        rng=np.random.default_rng(12); p=self.p
+        c=old.W[:,None,None]*(1+rng.normal(0,.002,(9,128,128)))
+        f=old.W[:,None,None]*(1+rng.normal(0,.002,(9,256,256)))
+        start=old.totals(c,f,p); c[:,p.covered]=np.nan
+        for _ in range(2):
+            f.reshape(9,-1)[:,p.ghost_fine]=oracle(c,p)
+            for _ in range(2):
+                f=np.array([np.roll(f[q],(ey,ex),(0,1)) for q,(ex,ey) in enumerate(old.E)])
+                f[:,p.fine_kind==OUTER_GHOST]=np.nan
+                f[:,p.fine_kind==ACTIVE]=old.collide(f[:,p.fine_kind==ACTIVE],p.tau_fine)
+            out=np.array([np.roll(c[q],(ey,ex),(0,1)) for q,(ex,ey) in enumerate(old.E)])
+            out.reshape(-1)[p.coalesce_coarse]=f.reshape(-1)[p.coalesce_fine].mean(axis=1)
+            out[:,~p.covered]=old.collide(out[:,~p.covered],p.tau_coarse)
+            out[:,p.covered]=np.nan; c=out
+            np.testing.assert_allclose(old.totals(c,f,p),start,rtol=0,atol=4e-10)
+
+    def test_private_build_identity_cleanup_and_old_path(self):
+        from CelerisLab.cuda import compiler_v2 as compiler
+        from assess_projected_pairs import BUILD_SOURCES
+        original=compiler.generate_config; report={}
+        with candidate.private_build(self.h,report,candidate.VARIANT):
+            for level,cfg in enumerate(self.h.configs):
+                cfg=replace(cfg,nx=10,ny=10);cfg.validate(); build=compiler.generate_config(cfg)
+                try:
+                    meta=json.loads((Path(build)/'.celeris-build').read_text())
+                    self.assertEqual(meta['source_sha256'],candidate.source_digest(build))
+                    self.assertEqual(meta['source_sha256'],candidate.BUILD_SHA256 if level else BUILD_SOURCES[128])
+                    if level:self.assertEqual(meta['experimental_reconstruction'],candidate.metadata(self.mask))
+                finally:compiler.cleanup_build_dir(build)
+        self.assertIs(compiler.generate_config,original)
+        with candidate.private_build(None,{},'current'):
+            self.assertIs(compiler.generate_config,original)
+        self.assertEqual(hashlib.sha256(SOURCE.read_bytes()).hexdigest(),'668e5cf14e88bdcfa281c2fa6c75c0f8f8cdbbfb8ba2ca03064ae02ef7ad0057')
+
+    def mock_records(self):
+        import argparse
+        helper=legacy.ProbeTraceTests();helper.setUp()
+        parse=argparse.ArgumentParser.parse_args
+        records=[]
+        for n,depth,variant in candidate.CASES:
+            scheme='constant' if depth==1 else 'projected_linear'
+            def parse_extra(parser,*a,**k):
+                result=parse(parser,*a,**k);result.test_reconstruction=variant;return result
+            def arrays(m,n,t):return [axis.initial(v.shape,l,n,.03*n/100) for l,v in enumerate(m)]
+            with patch.object(argparse.ArgumentParser,'parse_args',parse_extra), \
+                 patch.object(helper,'exact_arrays',side_effect=arrays), \
+                 patch.object(candidate,'private_build',return_value=nullcontext()) as private:
+                r,counts,_=helper.run_mocked(n,depth,candidate.DIAGNOSTIC,scheme,masks(n,depth),initializer=axis.initial)
+            self.assertEqual(r['test_reconstruction'],variant)
+            self.assertEqual(r['probe_trace']['steps'],list(range(9)))
+            self.assertEqual([v['step'] for v in r['rows']],[0,8*(n//128)])
+            self.assertLess(r['observed_shifted_max'],candidate.POP_ENVELOPE)
+            self.assertEqual(counts[-1],8*(n//128));self.assertEqual(private.call_args.args[2],variant)
+            records.append(r)
+        return records
+
+    def test_four_mock_runs_admission_and_old_diagnostic_unchanged(self):
+        self.assertEqual(len(self.mock_records()),4)
+        self.assertEqual(hashlib.sha256(Path(candidate.__file__).read_bytes()).hexdigest(),rc.SHEAR_CENTERED_SHA256)
+        from tests.validation.multigrid import check_projected_pairs as gate
+        with patch.object(gate,'preflight',return_value=rc.NUMERICAL_REVISION), \
+             patch.object(rc,'provenance',return_value=dict(source_tree_sha256=rc.SOURCE_SHA256)):
+            rc.verify_source(REPO,'r','shear_centered',probes=True)
+            with patch.object(candidate,'__file__','/stale/shear_centered.py'), self.assertRaises(RuntimeError):
+                rc.verify_source(REPO,'r','shear_centered',probes=True)
+        import argparse, io, sys
+        from contextlib import redirect_stdout
+        with tempfile.TemporaryDirectory() as tmp:
+            base=['run_case.py','--repo',str(REPO),'--expected-revision','f'*40,
+                  '--output',str(Path(tmp)/'r.json'),'--re','100','--n','128','--depth','2',
+                  '--flow','shear_wave','--probe-trace','--explosion','projected_linear','--dry-run']
+            cases=[(['--diagnostic',candidate.DIAGNOSTIC],False),
+                   (['--diagnostic',rc.AXIS_SQUARE,'--test-reconstruction',candidate.VARIANT],False),
+                   (['--diagnostic',candidate.DIAGNOSTIC,'--test-reconstruction',candidate.VARIANT,'--depth','1'],False),
+                   (['--diagnostic',candidate.DIAGNOSTIC,'--test-reconstruction',candidate.VARIANT],True)]
+            for extra,valid in cases:
+                with patch.object(sys,'argv',base+extra), patch.object(rc,'verify_source',return_value={}) as verify, \
+                     patch.object(rc,'runtime_source',return_value={}), redirect_stdout(io.StringIO()), patch('sys.stderr',io.StringIO()):
+                    if valid:self.assertEqual(rc.main(),0)
+                    else:
+                        with self.assertRaises(SystemExit):rc.main()
+                        verify.assert_not_called()
+
+    def test_frozen_comparison_gates_and_corruptions(self):
+        import test_projected_pairs as fixture
+        records=self.mock_records(); source=fixture.source()
+        for r in records:
+            r.update(source=deepcopy(source),runtime_identity=deepcopy(fixture.RUNTIME),
+                     runtime_source=rc.runtime_source(REPO),cuda_visible_devices='0',
+                     actual_build_sidecars=[fixture.build_metadata(r['n']*2**l) for l in range(r['depth'])])
+        records[3]['reconstruction']=candidate.metadata(self.mask)
+        records[3]['actual_build_sidecars'][1].update(source_sha256=candidate.BUILD_SHA256,
+                                                    experimental_reconstruction=candidate.metadata(self.mask))
+        def change(r,step,group,component,value):
+            for i,p in enumerate(r['probe_trace']['probes']):
+                if p['group']==group:r['probe_trace']['values'][step][3*i+component]=value
+        for r,bottom,top in [(records[2],-1e-5,1e-6),(records[3],-1e-6,8e-7)]:
+            change(r,1,'inner_bottom',1,bottom);change(r,1,'inner_top',1,top)
+        for r in records:r["probe_trace"]["bytes"]=axis.validate(r["probe_trace"])
+        result=candidate.compare(records,source)
+        self.assertTrue(result['advance_to_review']);self.assertFalse(result['qualified'])
+        np.testing.assert_allclose(result['moment_deadband'],[1.1925e-6,7.95e-7,7.95e-7],rtol=2e-16,atol=0)
+        for group,step,value,check in [('inner_bottom',1,-9.8e-6,'bottom_step1_clear'),
+                ('inner_top',1,4e-6,'top_step1_nonworse'),('outer_top',8,4e-6,'top_step8_nonworse'),
+                ('inner_left',1,4e-6,'vertical_null')]:
+            bad=deepcopy(records);change(bad[3],step,group,1,value)
+            bad[3]['probe_trace']['bytes']=axis.validate(bad[3]['probe_trace'])
+            self.assertFalse(candidate.compare(bad,source)['checks'][check])
+        for mutation in [lambda r:r['source'].update(revision='wrong'),
+                         lambda r:r['runtime_identity'].update(device='wrong'),
+                         lambda r:r.update(viscosity=.1),lambda r:r.update(refinements=[]),
+                         lambda r:r.update(observed_shifted_max=.04),
+                         lambda r:r['actual_build_sidecars'][1].update(source_sha256='wrong'),
+                         lambda r:r['probe_trace']['values'][0].__setitem__(0,1.01),
+                         lambda r:r['probe_trace']['steps'].pop(),
+                         lambda r:r['probe_trace']['probes'][0].update(i=1)]:
+            bad=deepcopy(records);mutation(bad[3])
+            with self.assertRaises((ValueError,KeyError)):candidate.compare(bad,source)

tests/validation/taylor_green/test_projected_pairs.py

@@ -429,7 +429,7 @@ class ProjectedPairsTests(unittest.TestCase):
                  patch.object(run_case, "runtime_identity", side_effect=AssertionError("CUDA reached")), \
                  redirect_stdout(stdout):
                 self.assertEqual(run_case.main(), 0)
-            verify.assert_called_once_with(REPO, MOCK_REVISION, "taylor_green")
+            verify.assert_called_once_with(REPO, MOCK_REVISION, "taylor_green", probes=False)
             self.assertTrue(json.loads(stdout.getvalue())["dry_run"])
             self.assertFalse(output.exists())
 
@@ -531,7 +531,7 @@ class ProjectedPairsTests(unittest.TestCase):
                     else:
                         self.assertEqual(run_case.main(), 0)
                 receipt = json.loads(output.read_text())
-                verify.assert_called_once_with(REPO, MOCK_REVISION, "taylor_green")
+                verify.assert_called_once_with(REPO, MOCK_REVISION, "taylor_green", probes=False)
                 self.assertEqual(selectors, [scheme] if depth == 2 else [])
                 self.assertEqual(receipt["source"], source())
                 self.assertEqual(receipt["runtime_identity"], RUNTIME)

tests/validation/taylor_green/test_shear_wave.py

@@ -1012,6 +1012,14 @@ class ProbeTraceTests(unittest.TestCase):
         trace["values"] = [[1.0, 0.01, -0.02] * 24 for _ in sp.STEPS]
         size = sp.validate(trace)
         self.assertLess(size, sp.TRACE_MAX_BYTES)
+        final_trace = dict(trace, bytes=size)
+        self.assertEqual(size, len(json.dumps(final_trace, allow_nan=False).encode()))
+        self.assertEqual(sp.validate(final_trace), size)
+        with patch.object(sp, "TRACE_MAX_BYTES", size):
+            self.assertEqual(sp.validate(trace), size)
+        with patch.object(sp, "TRACE_MAX_BYTES", size - 1):
+            with self.assertRaises(ValueError):
+                sp.validate(trace)
         realistic = dict(trace, values=[[1.0000012345678901, 0.012345678901234567,
                                          -0.023456789012345678] * 24 for _ in sp.STEPS])
         self.assertLess(sp.validate(realistic), sp.TRACE_MAX_BYTES)
@@ -1054,15 +1062,17 @@ class ProbeTraceTests(unittest.TestCase):
                         with self.assertRaises(SystemExit):
                             run_case.main()
 
-    def run_mocked(self, n, depth, diagnostic, scheme, masks):
+    def run_mocked(self, n, depth, diagnostic, scheme, masks, initializer=sw.initial):
         from types import SimpleNamespace
         from contextlib import nullcontext
         arrays = self.exact_arrays(masks, n, 0.0)
-        counts = []
+        counts, selectors = [], []
 
         class FakeRun:
             def __init__(self, *args, **kwargs):
                 self.step_count = 0
+                if "refinements" in kwargs:
+                    selectors.append((kwargs["explosion"], kwargs["refinements"]))
 
             def __enter__(self):
                 return self
@@ -1071,7 +1081,9 @@ class ProbeTraceTests(unittest.TestCase):
                 pass
 
             def initialize(self, *given):
-                pass
+                assert len(given) == depth
+                for level, (values, mask) in enumerate(zip(given, masks)):
+                    np.testing.assert_array_equal(values, initializer(mask.shape, level, n, U * n / sw.RE))
 
             def step(self, count):
                 self.step_count += count
@@ -1098,7 +1110,11 @@ class ProbeTraceTests(unittest.TestCase):
                  patch.object(run_case, "build_capture", return_value=nullcontext()), \
                  redirect_stdout(io.StringIO()):
                 run_case.main()
-            return json.loads(output.read_text()), counts, sorted(p.name for p in Path(tmp).iterdir())
+            rec = json.loads(output.read_text())
+            expected = ([(scheme, [[tuple(r) for r in level] for level in rec["refinements"]])]
+                        if depth > 1 else [])
+            self.assertEqual(selectors, expected)
+            return rec, counts, sorted(p.name for p in Path(tmp).iterdir())
 
     def test_mocked_probe_traces(self):
         sp = self.probes
@@ -1108,6 +1124,7 @@ class ProbeTraceTests(unittest.TestCase):
         self.assertEqual(trace["steps"], list(sp.STEPS))
         self.assertEqual(len(trace["values"]), len(sp.STEPS))
         self.assertLessEqual(trace["bytes"], sp.TRACE_MAX_BYTES)
+        self.assertEqual(trace["bytes"], len(json.dumps(trace, allow_nan=False).encode()))
         self.assertEqual(trace["case"], [run_case.NESTED_SHEAR, 128, 3, "projected_linear"])
         self.assertEqual(len(rec["rows"]), 5)  # qualification samples preserved
         self.assertEqual(files, ["probe.json"])  # no field dumps
@@ -1122,3 +1139,278 @@ class ProbeTraceTests(unittest.TestCase):
         self.assertTrue(all(p["level"] == 0 for p in t512["probes"]))
         res = sp.residuals(t512)
         self.assertLess(abs(res[0][0][1]), 1e-14)  # exact field at t=0 at actual centers
+
+    def test_invalid_optional_trace_preserves_scalar_rows(self):
+        sp = self.probes
+        case = (128, 1, run_case.NESTED_SHEAR, "constant", [np.ones((128, 128), bool)])
+        baseline, _, _ = self.run_mocked(*case)
+        original = sp.sample
+        for bad in ([float("nan")] * 72, [float("inf")] * 72, [1.0],
+                    ZeroDivisionError("zero probe density")):
+            calls = []
+
+            def sample(arrays, probes):
+                calls.append(True)
+                if len(calls) == 2:  # probe-only step 1, before the next scalar save
+                    if isinstance(bad, Exception):
+                        raise bad
+                    return bad
+                return original(arrays, probes)
+
+            with self.subTest(bad=repr(bad)), patch.object(sp, "sample", side_effect=sample):
+                rec, counts, files = self.run_mocked(*case)
+            self.assertTrue(rec["completed"])
+            self.assertEqual(rec["rows"], baseline["rows"])
+            self.assertEqual(set(rec["probe_trace"]), {"error"})
+            self.assertEqual(len(calls), 2)  # no further probe sampling after rejection
+            self.assertEqual(counts[-1], 2048)
+            self.assertEqual(files, ["probe.json"])
+        with patch.object(sp, "TRACE_MAX_BYTES", 100):
+            rec, _, _ = self.run_mocked(*case)
+        self.assertTrue(rec["completed"])
+        self.assertEqual(rec["rows"], baseline["rows"])
+        self.assertIn("exceeds", rec["probe_trace"]["error"])
+
+
+class WideXDiagnosticTests(unittest.TestCase):
+    """One width contrast; historical qualification matrices remain unchanged."""
+
+    def setUp(self):
+        import shear_nested
+        import shear_probes
+        import shear_wide_x
+        from CelerisLab.lbm.hierarchy import build_hierarchy
+        self.wide, self.nested, self.probes = shear_wide_x, shear_nested, shear_probes
+        self.h = build_hierarchy(run_case.config(128, 100), self.wide.REFINEMENTS)
+
+    def test_geometry_counts_and_legal_pair_tables(self):
+        from CelerisLab.lbm.blocks import BlockLayout2D
+        from CelerisLab.lbm.hierarchy import build_hierarchy, leaf_partition
+        from CelerisLab.lbm.refinement import ACTIVE
+        from CelerisLab.lbm.refinement_tables import projected_explosion_tables
+        wide, h = self.wide, self.h
+        self.assertEqual(wide.REFINEMENTS[0], self.nested.REFINEMENTS[0])
+        self.assertEqual(wide.REFINEMENTS[1], [(56, 44, 96, 212), (160, 60, 200, 228)])
+        for outer, inner in zip(wide.REFINEMENTS[0], wide.REFINEMENTS[1]):
+            x0, y0, x1, y1 = (v / 2 for v in inner)
+            self.assertEqual((x0 - outer[0], outer[2] - x1), (12, 12))
+            self.assertEqual((y0 - outer[1], outer[3] - y1), (6, 6))
+        # Independent physical rectangle masks at each native spacing.
+        masks = []
+        for level, kind in enumerate(h.kinds):
+            x, y = centers(kind.shape, level)
+            outer = ((x >= 16) & (x < 60) & (y >= 16) & (y < 112)
+                     | (x >= 68) & (x < 112) & (y >= 24) & (y < 120))
+            inner = ((x >= 28) & (x < 48) & (y >= 22) & (y < 106)
+                     | (x >= 80) & (x < 100) & (y >= 30) & (y < 114))
+            masks.append(~outer if level == 0 else outer & ~inner if level == 1 else inner)
+            np.testing.assert_array_equal(kind == ACTIVE, masks[-1])
+        self.assertEqual([int(m.sum()) for m in masks], [7936, 20352, 53760])
+        np.testing.assert_array_equal(leaf_partition(h), np.ones((128, 128)))
+        self.assertEqual(run_case.workload(masks, 2048)["active_collision_updates"], 540016640)
+        self.assertEqual(wide.LEAVES, [7936, 20352, 53760])
+        self.assertEqual(wide.HOPS, [2048, 4096, 8192])
+        self.assertEqual(wide.WORK, 540016640)
+        np.testing.assert_allclose([c.viscosity for c in h.configs], wide.NU, rtol=0, atol=0)
+        np.testing.assert_allclose([1 / c.omega for c in h.configs], wide.TAU, rtol=1e-15)
+        old = build_hierarchy(run_case.config(128, 100), self.nested.REFINEMENTS)
+        for name in ("ghost_coarse", "ghost_fine", "coalesce_coarse", "coalesce_fine"):
+            np.testing.assert_array_equal(getattr(h.pairs[0], name), getattr(old.pairs[0], name))
+        root = BlockLayout2D.build(128, 128, 8, 8, periodic=(True, True))
+        for level, pair in enumerate(h.pairs):
+            layout = root if level == 0 else h.pairs[level - 1].fine_layout((8, 8))
+            neighbors, offsets = projected_explosion_tables(pair, layout, h.kinds[level])
+            legal = np.unique(pair.ghost_coarse)
+            legal = legal[(h.kinds[level].ravel()[legal] == ACTIVE) & ~pair.covered.ravel()[legal]]
+            self.assertTrue(np.isin(neighbors, layout.owner.ravel()[legal]).all())
+            self.assertEqual(neighbors.shape, (len(pair.ghost_fine), 4))
+            self.assertEqual(len(np.unique(pair.ghost_fine)), len(pair.ghost_fine))
+            self.assertTrue((np.abs(offsets) == .25).all())
+            if level == 0:
+                expected = projected_explosion_tables(old.pairs[0], root, old.kinds[0])
+                for a, b in zip((neighbors, offsets), expected):
+                    np.testing.assert_array_equal(a, b)
+
+    def test_fixed_probes_and_complete_runner_schedule(self):
+        from CelerisLab.lbm.refinement import ACTIVE
+        masks = [k == ACTIVE for k in self.h.kinds]
+        probes = self.probes.select(masks)
+        old = self.probes.select(self.nested.leaf_masks())
+        self.assertEqual([p["level"] for p in probes[:12]], [0] * 3 + [1] * 9)
+        self.assertLess(max(p["center"][0] for p in probes[:12]), 28)
+        self.assertEqual(probes[:8], old[:8])  # identical-center old L0/L1 probes
+        self.assertEqual(probes[12:], old[12:])  # unchanged y-normal transect
+        helper = ProbeTraceTests()
+        helper.setUp()
+        rec, counts, files = helper.run_mocked(128, 3, run_case.WIDE_X_SHEAR,
+                                               "projected_linear", masks)
+        self.assertTrue(rec["completed"])
+        self.assertEqual(rec["refinements"], [[list(r) for r in level] for level in self.wide.REFINEMENTS])
+        self.assertEqual(rec["workload"]["active_cells_by_level"], self.wide.LEAVES)
+        self.assertEqual(rec["workload"]["active_collision_updates"], self.wide.WORK)
+        self.assertEqual([r["step"] for r in rec["rows"]], [0, 512, 1024, 1536, 2048])
+        self.assertEqual(rec["probe_trace"]["steps"], list(self.probes.STEPS))
+        self.assertEqual(rec["probe_trace"]["probes"], probes)
+        self.assertEqual(counts[-1], 2048)
+        self.assertEqual(files, ["probe.json"])
+        self.assertNotIn("final_snapshot", rec)
+
+    def test_one_case_admission_and_pin(self):
+        self.assertEqual(run_case.SHEAR_WIDE_X_SHA256,
+                         hashlib.sha256((HERE / "shear_wide_x.py").read_bytes()).hexdigest())
+        repo = HERE.parents[2]
+        with tempfile.TemporaryDirectory() as tmp:
+            base = ["run_case.py", "--repo", str(repo), "--expected-revision", "f" * 40,
+                    "--output", str(Path(tmp) / "case.json"), "--dry-run", "--re", "100",
+                    "--flow", "shear_wave", "--diagnostic", run_case.WIDE_X_SHEAR,
+                    "--n", "128", "--depth", "3", "--explosion", "projected_linear", "--probe-trace"]
+            for extra in ([], ["--n", "512"], ["--depth", "2"], ["--explosion", "constant"],
+                          ["--re", "1000"], ["--flow", "taylor_green"], ["--final-snapshot"],
+                          ["--early-pressure"]):
+                with self.subTest(extra=extra), patch.object(sys, "argv", base + extra), \
+                     patch.object(run_case, "verify_source", return_value={}) as verify, \
+                     patch.object(run_case, "runtime_source", return_value={}), \
+                     redirect_stdout(io.StringIO()), patch("sys.stderr", io.StringIO()):
+                    if not extra:
+                        self.assertEqual(run_case.main(), 0)
+                        verify.assert_called_once_with(repo, "f" * 40, "shear_wide_x", probes=True)
+                    else:
+                        with self.assertRaises(SystemExit):
+                            run_case.main()
+                        verify.assert_not_called()
+            with patch.object(sys, "argv", base[:-1]), patch("sys.stderr", io.StringIO()):
+                with self.assertRaises(SystemExit):
+                    run_case.main()
+
+    def test_historical_matrix_rejects_width_case(self):
+        helper = NestedPresetTests()
+        helper.setUp()
+        records = helper.nested_set()
+        self.assertTrue(helper.mod.assess_nested(records, "r")["passed"])
+        records[-1].update(diagnostic=run_case.WIDE_X_SHEAR,
+                           refinements=self.wide.REFINEMENTS)
+        self.assertEqual(helper.mod.assess_nested(records, "r")["status"], "INVALID_SET")
+        self.assertNotIn(run_case.WIDE_X_SHEAR, helper.mod.PRESETS)
+
+
+class PostGDiagnosticTests(unittest.TestCase):
+    """Three initialization-only cases, excluded from historical qualification."""
+
+    def setUp(self):
+        import shear_postg
+        self.postg = shear_postg
+
+    def test_moments_phase_and_physical_scaling(self):
+        from CelerisLab.lbm.population_storage import _encode_physical, _decode_stored
+        velocities = np.array([CX, CY], dtype=float).T
+        for level in range(3):
+            shape, nu = (31, 47), .0384
+            h, tau = 2.0**-level, .5 + 3 * nu * 2**level
+            eq = sw.initial(shape, level, 128, nu)
+            g = self.postg.initial(shape, level, 128, nu)
+            x, y = centers(shape, level)
+            k = sw.wavevector(128)
+            strain = (.5 * U * np.cos(k[0] * x + k[1] * y + sw.PHI)[None, None]
+                      * (np.outer(sw.E, h * k) + np.outer(h * k, sw.E))[:, :, None, None])
+            tensor = np.einsum("qi,qj->qij", velocities, velocities) - np.eye(2)[None] / 3
+            a = 3 * oracle.W[:, None, None] * np.einsum("qij,ijyx->qyx", tensor, strain)
+            np.testing.assert_allclose(g, eq - (tau - 1) * a, rtol=0, atol=6e-17)
+            for weights in (np.ones(9), CX, CY):
+                np.testing.assert_allclose(np.einsum("q,qyx->yx", weights, g - eq),
+                                           0, rtol=0, atol=2e-16)
+            stress = np.einsum("qi,qj,qyx->ijyx", velocities, velocities, g - eq)
+            np.testing.assert_allclose(stress, (1 - tau) * 2 / 3 * strain, rtol=0, atol=6e-17)
+            pre = eq - tau * a
+            rho = pre.sum(axis=0)
+            pre_eq = sw.equilibrium(rho, np.einsum("q,qyx->yx", CX, pre) / rho,
+                                   np.einsum("q,qyx->yx", CY, pre) / rho)
+            np.testing.assert_allclose((1 - 1 / tau) * pre + pre_eq / tau, g, rtol=0, atol=4e-16)
+            self.assertAlmostEqual((tau - .5) * h, .1152, places=16)
+            self.assertGreater(1 - tau, 0)
+            encoded = _decode_stored(_encode_physical(g, True), True)
+            self.assertGreater(encoded.min(), 0)
+            for weights in (np.ones(9), CX, CY):
+                np.testing.assert_allclose(np.einsum("q,qyx->yx", weights, encoded - eq),
+                                           0, rtol=0, atol=3e-9)
+        # The U512 control is the same physical spacing/tau/phase as nested L2.
+        np.testing.assert_array_equal(self.postg.initial((512, 512), 0, 512, .1536),
+                                      self.postg.initial((512, 512), 2, 128, .0384))
+
+    def test_actual_initializers_and_three_complete_schedules(self):
+        helper = ProbeTraceTests()
+        helper.setUp()
+        work = 0
+        for _, n, depth, scheme in self.postg.CASES:
+            masks = helper.nested.leaf_masks() if depth == 3 else [np.ones((n, n), bool)]
+            rec, counts, files = helper.run_mocked(n, depth, run_case.POSTG_SHEAR, scheme, masks,
+                                                   initializer=self.postg.initial)
+            self.assertEqual(rec["initialization"], self.postg.INITIALIZATION)
+            self.assertEqual(rec["initialization_id"], "postg_first_gradient_v1")
+            self.assertFalse(rec["qualification_eligible"])
+            self.assertTrue(rec["completed"])
+            self.assertEqual(rec["limits"], sw.LIMITS)
+            self.assertEqual([r["step"] for r in rec["rows"]], run_case.sample_steps(n)[0])
+            self.assertEqual(rec["probe_trace"]["steps"], list(helper.probes.STEPS))
+            self.assertEqual(rec["probe_trace"]["probes"], helper.probes.select(masks, n))
+            self.assertEqual(counts[-1], 16 * n)
+            self.assertIn(1920 * (n // 128), counts)
+            self.assertEqual(files, ["probe.json"])
+            self.assertNotIn("final_snapshot", rec)
+            work += rec["workload"]["active_collision_updates"]
+        self.assertEqual(work, 2952265728)
+
+    def test_fixed_admission_and_source_pin(self):
+        from tests.validation.multigrid import check_projected_pairs as gate
+        repo = HERE.parents[2]
+        self.assertEqual(run_case.SHEAR_POSTG_SHA256,
+                         hashlib.sha256((HERE / "shear_postg.py").read_bytes()).hexdigest())
+        with tempfile.TemporaryDirectory() as tmp:
+            base = ["run_case.py", "--repo", str(repo), "--expected-revision", "f" * 40,
+                    "--output", str(Path(tmp) / "case.json"), "--dry-run", "--re", "100",
+                    "--flow", "shear_wave", "--diagnostic", run_case.POSTG_SHEAR, "--probe-trace"]
+            for n, depth, scheme, valid in ((128, 1, "constant", True), (512, 1, "constant", True),
+                                            (128, 3, "projected_linear", True),
+                                            (128, 3, "constant", False), (256, 1, "constant", False),
+                                            (128, 2, "projected_linear", False)):
+                argv = base + ["--n", str(n), "--depth", str(depth), "--explosion", scheme]
+                extras = ([], ["--final-snapshot"], ["--early-pressure"], ["--re", "1000"],
+                          ["--flow", "taylor_green"]) if valid else ([],)
+                for extra in extras:
+                    with self.subTest(argv=argv + extra), patch.object(sys, "argv", argv + extra), \
+                         patch.object(run_case, "verify_source", return_value={}) as verify, \
+                         patch.object(run_case, "runtime_source", return_value={}), \
+                         redirect_stdout(io.StringIO()), patch("sys.stderr", io.StringIO()):
+                        if valid and not extra:
+                            self.assertEqual(run_case.main(), 0)
+                            verify.assert_called_once_with(repo, "f" * 40, "shear_postg", probes=True)
+                        else:
+                            with self.assertRaises(SystemExit):
+                                run_case.main()
+                            verify.assert_not_called()
+            with patch.object(sys, "argv", base[:-1] + ["--n", "128", "--depth", "1"]), \
+                 patch("sys.stderr", io.StringIO()):
+                with self.assertRaises(SystemExit):
+                    run_case.main()
+        with patch.object(gate, "preflight", return_value=run_case.NUMERICAL_REVISION), \
+             patch.object(run_case, "provenance", return_value=dict(source_tree_sha256=run_case.SOURCE_SHA256)):
+            self.assertEqual(run_case.verify_source(repo, "r", "shear_postg")["shear_postg_sha256"],
+                             run_case.SHEAR_POSTG_SHA256)
+            for module, attr, value in ((self.postg, "__file__", "/wrong/shear_postg.py"),
+                                         (run_case, "SHEAR_POSTG_SHA256", "wrong")):
+                with patch.object(module, attr, value), self.assertRaisesRegex(RuntimeError, "shear_postg"):
+                    run_case.verify_source(repo, "r", "shear_postg")
+
+    def test_old_assessors_reject_identity_or_initializer(self):
+        from copy import deepcopy
+        helper = NestedPresetTests()
+        helper.setUp()
+        for base, assessor in ((helper.nested_set(), helper.mod.assess_nested),
+                               ([receipt(128, 1), receipt(256, 1), receipt(128, 2),
+                                 receipt(128, 2, "projected_linear")], helper.mod.assess)):
+            self.assertTrue(assessor(base, "r")["passed"])
+            for change in (dict(diagnostic=run_case.POSTG_SHEAR),
+                           dict(initialization=self.postg.INITIALIZATION)):
+                records = deepcopy(base)
+                records[-1].update(change)
+                self.assertEqual(assessor(records, "r")["status"], "INVALID_SET")
+        self.assertNotIn(run_case.POSTG_SHEAR, helper.mod.PRESETS)