test: compare fixed centered reconstruction on audited square faces

docs/validation/taylor_green.md

@@ -548,3 +548,59 @@ 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
@@ -29,6 +30,8 @@ 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
@@ -127,7 +130,14 @@ 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_axis_square":
+    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
@@ -162,7 +172,7 @@ def verify_source(repo, revision, flow="taylor_green", probes=False):
                 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", "shear_wide_x", "shear_postg", "shear_axis_square"):
+    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
@@ -237,7 +247,7 @@ def main():
     p.add_argument("--explosion", choices=["constant", "projected_linear"], default="constant")
     p.add_argument("--diagnostic",
                    choices=[SCALED_DIAGNOSTIC, SCALED_SHEAR, NESTED_SHEAR, TWIN_SHEAR,
-                            WIDE_X_SHEAR, POSTG_SHEAR, AXIS_SQUARE],
+                            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) "
@@ -245,6 +255,8 @@ def main():
                         "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, wide-x projected d3, POST-G three-case set, "
@@ -256,7 +268,16 @@ def main():
                    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()
-    axis_square = args.diagnostic == AXIS_SQUARE
+    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
     wide_x = args.diagnostic == WIDE_X_SHEAR
     postg = args.diagnostic == POSTG_SHEAR
@@ -273,7 +294,7 @@ def main():
                    (POSTG_SHEAR, 128, 3, "projected_linear")}
     if axis_square:
         import shear_axis_square
-        probe_cases.update((AXIS_SQUARE, n, depth, scheme)
+        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):
@@ -325,7 +346,7 @@ 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_flow = ("shear_axis_square" if axis_square else "shear_postg" if postg else "shear_wide_x" if wide_x else "shear_nested" if nested
+    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)
@@ -401,6 +422,10 @@ def main():
         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")
@@ -415,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(
                 *[
@@ -452,6 +478,8 @@ def main():
                     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:

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/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)