test: shear-wave modal-extinction handling; TG mock verify_source flow arg

- shear_wave.fit_rate: finite late R=0 returns rate_fit_valid=false with
  null fit quantities (no clipping, no invented rate); own screens and
  budgets treat it as a failure (CONTROL_SCREEN_FAIL / refined FAIL),
  other finite metrics retained, JSON-safe.
- test_projected_pairs: dry-run/execution mocks expect the explicit
  'taylor_green' flow argument now passed to verify_source.
- run_case: SHEAR_WAVE_SHA256 repinned; extinction test added.

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

docs/validation/taylor_green.md

@@ -159,3 +159,139 @@ theoretical stress-consistency certification or production-kernel change.
 ## Repository integration
 
 The four Python harness files remain together under tests/validation/taylor_green, so runtime provenance hashes precisely this harness directory. Formatting changes produce a new harness identity; preparation and integrated receipts must not be mixed. All analytic definitions, numerical floors, assessment gates, sample times, and source-import assertions are unchanged. Production source is unchanged from the cd3 snapshot; these periodic tests do not waive the separate channel GPU gate.
+
+## Oblique shear wave (test-only diagnostic, not GPU-qualified)
+
+`shear_wave.py` reuses the frozen oracle's constants, cell centers and
+geometry. `oracle.py` is byte-identical, and `test_shear_wave.py` pins its
+eec4677 SHA-256. `run_case.py --flow shear_wave` selects the shear wave;
+the default flow and its receipt are the unchanged Taylor–Green path.
+
+### Convention
+
+Integer mode (2,1) on the L=N periodic root domain:
+- the wavevector `k = (2*pi/L)*(2,1)` points 26.565° from +x;
+- mode (1,2) would point 63.435°, and a wrong-mode test catches a swap.
+
+Polarization is `E = (-1,2)/sqrt(5)`, perpendicular to k:
+
+    u = U*E*sin(k.x + phi)*exp(-nu*|k|^2*t),  phi = 0,  p = const,  rho0 = 1
+
+`nu = U*N/100`, the Taylor–Green Re100 viscosity at the same N. This is an
+exact incompressible Navier–Stokes solution. The initial condition is FP64
+second-order equilibrium at the actual leaf centers, with no
+non-equilibrium part.
+
+### Physical caveats
+
+1. *Scope of the cubic-term argument.* The second-order D2Q9 equilibrium
+   lacks the cubic third moment `rho*u_a*u_b*u_g` of the Maxwellian (Dellar,
+   https://people.maths.ox.ac.uk/dellar/papers/Galilean.pdf, eq. 3). At
+   Navier–Stokes order this produces a stress error proportional to
+   `d_g(rho u_a u_b u_g)`.
+
+   For the ideal constant-density, zero-background transverse wave, k.E = 0
+   makes that contracted derivative vanish. That is the only claim made.
+   It does not remove other discrete or finite-Mach errors. It is not
+   guaranteed once AMR introduces density, longitudinal or higher-mode
+   contamination. The analytic NS field itself is exact.
+2. *Startup.* Equilibrium initialization is not the discrete shear
+   eigenmode (Mei et al., https://doi.org/10.1016/j.compfluid.2005.08.008).
+   - The decay rate is a free-intercept least-squares fit of ln R over the
+     samples with t* >= 0.12, where `R = hypot(a_hat, b_hat)`.
+   - The intercept cancels only a constant amplitude prefactor, and only
+     once a single exponential dominates.
+   - Per-interval rates and the maximum log-fit residual are kept, so a
+     lingering transient stays visible.
+   - Late sampling alone does not show that transients have gone.
+
+### Metrics
+
+Computed on ACTIVE leaf cells, each weighted by its area `4**-level` in
+root-cell units. The partition must cover the domain exactly once, and one
+global pressure gauge is used.
+
+- `a_hat`, `b_hat`: weighted least-squares fit of `u.E` onto sin/cos, using
+  the 2x2 Gram matrix of the actual composite quadrature. A wrong Fourier
+  mode need not project to exactly zero on this quadrature.
+- `amplitude_error = R/A_exact - 1`.
+- `phase_drift = atan2(b_hat, a_hat)`, in radians, reported separately.
+- Rates are positive (minus the fitted slope).
+  - Per root step they are compared with `nu*|k|^2`. Uniform 2N uses its
+    own `nu = 0.0768` at N=256 and doubled step counts.
+  - In t* they are compared with `5*(2*pi)^2/Re = 1.97392088`.
+  - `rate_error = |rate/exact - 1|`.
+- `velocity_l2`, `residual_l2` (velocity left after removing the fitted
+  mode) and `longitudinal_l2` (the `u.k_hat` part), all relative to the
+  exact-field norm. A longitudinal perturbation also raises the
+  full-vector residual and the velocity L2.
+- `spurious_pressure_rms`: the raw gauged rms of `p = (rho-1)/3`.
+- `spurious_pressure_scaled`: that rms divided by `A_exact(t)^2/4`, a
+  continuum reference scale. It is *not* the discrete, globally gauged
+  Taylor–Green norm on the composite grid. Exact pressure is constant, so
+  no relative-to-exact claim is made. Max and per-level values are kept.
+
+### Preregistered screens (`shear_wave.assess`)
+
+Run set: uniform N, uniform 2N, and refined depth-2 runs at root N.
+
+- `eps(N) = (|k| h_root)^2 = 5*(2*pi/N)^2`: dimensionless, coefficient
+  one, 0.01205 at N=128. This is a chosen screen, not a proven error bound.
+- Absolute screens:
+  - rate_error, max |phase_drift| (radians) and max residual_l2 must each
+    be <= eps;
+  - velocity_l2 <= 0.05, mass <= 1e-5, spurious_pressure_scaled <= 0.02.
+- Budgets keep the `assess.py` form, `refined <= max(2*control_N, floor)`:
+  - floors are 0.002 for rate and residual, 0.02 for pressure;
+  - the comparative pressure budget is always redundant:
+    `max(2*control, 0.02)` is never below the 0.02 absolute ceiling;
+  - raw ratios to both the root and finest controls, including raw
+    pressure rms, are always reported.
+- Every record is validated before any screen is applied:
+  - exact scheduled steps and t*, with no duplicates;
+  - every numeric field finite, including list fields and density;
+  - error quantities nonnegative and rho_min > 0;
+  - signed quantities (amplitudes, `amplitude_error`, `phase_drift`,
+    per-interval rate errors) keep their sign.
+- Outcomes are kept distinct:
+  - INVALID_CONTROL: a control or refined run is missing, a control is
+    incomplete or invalid, or the records come from mixed sources;
+  - INVALID_RUN: a refined record fails validation;
+  - CONTROL_SCREEN_FAIL: complete, valid controls fail their own screens.
+    Each refined run's `own_screens_pass` is still reported, but
+    `controlled_pass` is false, so it is never a controlled-comparison pass;
+  - PASS/FAIL: every refined run's `controlled_pass` (own screens and
+    budgets against valid, passing controls).
+  - Finite modal extinction (late R = 0) is a physical failure, not
+    malformed data. The run gets `rate_fit_valid=false` with null fit
+    quantities, and its other metrics are kept. Nothing is clipped and no
+    rate is invented. The result is CONTROL_SCREEN_FAIL for a control, or
+    FAIL for a refined run.
+  - Screens are never retuned.
+
+### Limitations
+
+One flow, one phase and one rectangle, SRT only. The rectangle's vertical
+faces span 0.625 tangential periods and its horizontal faces 1.5, so this
+is not placement robustness. Two resolutions show a trend, not an order.
+
+### Receipt assessment and runner integrity
+
+- Fixed four-case N128 contract:
+  `python assess_shear_wave.py U128.json U256.json R128_constant.json R128_projected.json --output report.json --expected-revision <sha>`.
+  It reuses the shared provenance, build, workload and sample-time checks
+  from `assess_projected_pairs.py`, and adds shear metadata checks (flow,
+  mode, polarization, angle, phase, eps, limits, disabled early trace, and
+  pinned `shear_wave.py` hash).
+  - Missing, extra, duplicate or mixed-source records give `INVALID_SET`,
+    as does any failed shared check.
+  - Otherwise the result comes from `shear_wave.assess`, which separates
+    `INVALID_CONTROL` from `CONTROL_SCREEN_FAIL`.
+- `run_case.verify_source(..., flow="shear_wave")` requires the imported
+  `shear_wave.py` to be the selected harness file with the pinned
+  `SHEAR_WAVE_SHA256`. The default TG flow does not import it.
+- Sample and READY finiteness use `finite_row`, which accepts scalar and
+  list diagnostics. For TG's all-scalar rows it behaves the same as before.
+- Merge note: native 92b8afe adds a `scaled_n256` receipt block. When the
+  two are merged, keep both blocks and make `--flow shear_wave` and
+  `scaled_n256` mutually exclusive. Not merged here.

tests/validation/taylor_green/assess_shear_wave.py

@@ -0,0 +1,81 @@
+"""Fixed N128 four-case shear-wave contract: U128, U256, constant R128, projected R128.
+
+Shared provenance, build, workload and sample-time checks are reused from
+assess_projected_pairs (same cases, source pins and geometry); only the
+Taylor-Green-specific checks are replaced by shear-wave metadata and physics.
+"""
+
+import argparse
+import json
+from pathlib import Path
+
+import assess_projected_pairs as tg
+import shear_wave as sw
+from run_case import SHEAR_WAVE_SHA256
+
+SHARED = ("complete_matrix", "no_extra_cases", "no_duplicate_cases", "finite_receipts",
+          "same_source_and_harness", "clean_a135252_numerical_source",
+          "same_recorded_runtime", "same_explicit_device", "same_early_sampling")
+PER_CASE = ("complete_times", "scales", "actual_builds", "reference", "runtime_source",
+            "patch", "workload")
+
+
+def metadata_ok(record, n):
+    return (record.get("flow") == "shear_wave" and record.get("mode") == list(sw.MODE)
+            and record.get("polarization") == sw.E.tolist()
+            and record.get("angle_deg") == sw.ANGLE_DEG and record.get("phase") == sw.PHI
+            and record.get("eps") == sw.eps(n)
+            and record.get("numerics") == tg.NUMERICS
+            and record.get("initialization") == sw.INITIALIZATION
+            and record.get("limits") == sw.LIMITS
+            and record.get("early_pressure") == dict(enabled=False, rows=[],
+                                                     sampled_max_pressure_l2=None)
+            and (record.get("source") or {}).get("shear_wave_sha256") == SHEAR_WAVE_SHA256)
+
+
+def assess(records, expected_revision=None):
+    if not isinstance(records, list) or not all(isinstance(r, dict) for r in records):
+        raise ValueError("records must be a list of receipt objects")
+    base_checks = tg.assess(records, expected_revision)["checks"]
+    checks = {name: base_checks.get(name, False) for name in SHARED}
+    indexed = {}
+    for r in records:
+        indexed.setdefault(tuple(r.get(k) for k in ("re_L", "n", "depth", "explosion")), r)
+    for key in tg.CASES:
+        _, n, depth, scheme = key
+        label = f"{n}/{depth}/{scheme}"
+        for name in PER_CASE:
+            checks[f"{label} {name}"] = base_checks.get(f"{label} {name}", False)
+        checks[f"{label} shear_metadata"] = key in indexed and metadata_ok(indexed[key], n)
+    contract_valid = all(checks.values())
+    result = dict(checks=checks, contract_valid=contract_valid,
+                  cases=[list(k) for k in tg.CASES], expected_revision=expected_revision,
+                  scope="2D periodic Re100 oblique (2,1) shear wave; one phase and one "
+                        "rectangle; no placement robustness or convergence order")
+    if not contract_valid:
+        result.update(status="INVALID_SET", passed=False)
+        return result
+    try:
+        physics = sw.assess(indexed[tg.CASES[0]], indexed[tg.CASES[1]],
+                            [indexed[tg.CASES[2]], indexed[tg.CASES[3]]])
+    except (KeyError, ValueError) as error:
+        result.update(status="INVALID_SET", passed=False, reason=str(error))
+        return result
+    result.update(physics=physics, status=physics["status"],
+                  passed=physics["status"] == "PASS")
+    return result
+
+
+if __name__ == "__main__":
+    p = argparse.ArgumentParser(description=__doc__)
+    p.add_argument("files", type=Path, nargs=4)
+    p.add_argument("--output", type=Path, required=True)
+    p.add_argument("--expected-revision", required=True)
+    args = p.parse_args()
+    if args.output.exists():
+        p.error("Output already exists; choose a new path to retain prior evidence")
+    result = assess([json.loads(path.read_text()) for path in args.files], args.expected_revision)
+    text = json.dumps(result, indent=2, allow_nan=False) + "\n"
+    args.output.write_text(text)
+    print(text, end="")
+    raise SystemExit(0 if result["passed"] else 1)

tests/validation/taylor_green/run_case.py

@@ -18,6 +18,30 @@ SOURCE_SHA256 = "a5dbfd804befce868067d79a23e5fd47dce7f505b3f9b1bb2145fb57cb5f217
 BUILD_CAPTURE_SHA256 = "3dbaa6547cb0a4cb1541054f3b64ad8b6e18cc70527ac6166e56b944e9bd464a"
 REFERENCE_HARNESS = "11901a7f77557a901e017a171eeae1758d584fa5"
 INITIALIZATION = "physical FP64 equilibrium G at actual cell centers; rho=1+3p"
+SHEAR_WAVE_SHA256 = "81d5c8ee60864372887cddf9239580c82fb6750d956038f1e3ee9de61aa4e877"
+
+
+def flow_functions(flow="taylor_green"):
+    """(initial, metrics, limits, initialization); the default is the frozen TG oracle."""
+    if flow == "taylor_green":
+        return initial, metrics, LIMITS, INITIALIZATION
+    if flow == "shear_wave":
+        import shear_wave
+        return shear_wave.initial, shear_wave.metrics, shear_wave.LIMITS, shear_wave.INITIALIZATION
+    raise ValueError(f"unknown flow {flow!r}")
+
+
+def finite_row(row):
+    """Every scalar or list-valued diagnostic is finite (TG rows are all scalars)."""
+    return all(np.isfinite(np.asarray(value, dtype=float)).all() for value in row.values())
+
+
+def sample(run, masks, step, n, viscosity, initial_mass, metrics, depth, label="sample"):
+    arrays = (run.download(),) if depth == 1 else run.download()
+    row = metrics(arrays, masks, step, n, viscosity, initial_mass)
+    if not finite_row(row):
+        raise FloatingPointError(f"nonfinite {label} diagnostics")
+    return row
 
 
 def sample_steps(n, early_pressure=False):
@@ -64,7 +88,7 @@ def provenance(repo):
     )
 
 
-def verify_source(repo, revision):
+def verify_source(repo, revision, flow="taylor_green"):
     """Exact clean harness commit, fixed production tree, and reviewed capture helper."""
     import oracle
     from tests.validation.multigrid import check_projected_pairs as gate
@@ -83,6 +107,13 @@ def verify_source(repo, revision):
     if result["source_tree_sha256"] != SOURCE_SHA256:
         raise RuntimeError("Wrong numerical source fingerprint")
     result["build_capture_sha256"] = BUILD_CAPTURE_SHA256
+    if flow == "shear_wave":
+        import shear_wave
+        module = harness / "shear_wave.py"
+        if (Path(shear_wave.__file__).resolve() != module
+                or hashlib.sha256(module.read_bytes()).hexdigest() != SHEAR_WAVE_SHA256):
+            raise RuntimeError("Wrong shear_wave module checkout or content")
+        result["shear_wave_sha256"] = SHEAR_WAVE_SHA256
     return result
 
 
@@ -151,7 +182,11 @@ def main():
     p.add_argument("--explosion", choices=["constant", "projected_linear"], default="constant")
     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()
+    if args.flow == "shear_wave" and (args.re != 100 or args.depth == 3 or args.early_pressure):
+        p.error("shear_wave requires --re 100, depth 1 or 2 and no --early-pressure")
+    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:
@@ -161,7 +196,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 = verify_source(repo, args.expected_revision)
+    source = verify_source(repo, args.expected_revision, args.flow)
     runtime = runtime_source(repo)
     if args.dry_run:
         print(json.dumps(dict(source=source, runtime_source=runtime, dry_run=True)))
@@ -197,7 +232,7 @@ def main():
         mach=U * np.sqrt(3),
         viscosity=cfg.viscosity,
         tau_by_level=[0.5 + 3 * c.viscosity for c in h.configs],
-        initialization=INITIALIZATION,
+        initialization=initialization,
         numerics=dict(lattice_model=cfg.lattice_model, collision=cfg.collision,
                       store_precision=cfg.store_precision, ddf_shifting=cfg.ddf_shifting,
                       streaming=cfg.streaming, les_enabled=cfg.les_enabled,
@@ -205,13 +240,18 @@ def main():
                       sample_phase="postcollision G at root synchronization"),
         cuda_visible_devices=selected,
         device_id=0,
-        limits=LIMITS,
+        limits=limits,
         workload=workload(masks, round(TIMES[-1] * args.n / U)),
         rows=[],
         early_pressure=dict(enabled=args.early_pressure, rows=[],
                             sampled_max_pressure_l2=None),
         completed=False,
     )
+    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))
     args.output.parent.mkdir(parents=True, exist_ok=True)
     def save():
         temporary = args.output.with_suffix(args.output.suffix + ".tmp")
@@ -248,10 +288,8 @@ def main():
             initial_mass = None
             for step in sorted(set(main_steps + early_steps)):
                 run.step(step - run.step_count)
-                arrays = (run.download(),) if args.depth == 1 else run.download()
-                row = metrics(arrays, masks, step, args.n, cfg.viscosity, initial_mass)
-                if not all(np.isfinite(value) for value in row.values()):
-                    raise FloatingPointError("nonfinite sample diagnostics")
+                row = sample(run, masks, step, args.n, cfg.viscosity, initial_mass,
+                             metrics, args.depth)
                 row["step"] = step
                 if initial_mass is None:
                     initial_mass = row["mass"]
@@ -278,10 +316,8 @@ def main():
                     # Exercise real kernels/halos/transfers on every level before READY.
                     # This is the first measured trajectory step, not discarded warmup.
                     run.step(1)
-                    probe = (run.download(),) if args.depth == 1 else run.download()
-                    probe_metrics = metrics(probe, masks, 1, args.n, cfg.viscosity, initial_mass)
-                    if not all(np.isfinite(value) for value in probe_metrics.values()):
-                        raise FloatingPointError("nonfinite READY diagnostics")
+                    sample(run, masks, 1, args.n, cfg.viscosity, initial_mass,
+                           metrics, args.depth, "READY")
                     report["ready_after_step"] = 1
                     save()
                     print(
@@ -301,7 +337,7 @@ def main():
         completed=True,
         seconds=time.monotonic() - started,
         screening_pass=all(
-            row[key] <= limit for row in report["rows"] for key, limit in LIMITS.items()
+            row[key] <= limit for row in report["rows"] for key, limit in limits.items()
         ),
     )
     save()

tests/validation/taylor_green/shear_wave.py

@@ -0,0 +1,312 @@
+"""CPU-only oblique shear-wave reference and composite-leaf diagnostics.
+
+Convention: integer mode MODE=(2,1) on the L x L periodic root domain, so the
+wavevector k = (2*pi/L)*(2,1) points 26.565 deg from +x (mode (1,2) would be
+63.435 deg). Velocity is polarized along E = (-1,2)/sqrt(5), perpendicular to k.
+The frozen Taylor-Green oracle is reused for constants and geometry only.
+"""
+
+import math
+
+import numpy as np
+from oracle import CX, CY, W, U, TIMES, centers
+
+MODE = (2, 1)
+PHI = 0.0
+E = np.array([-MODE[1], MODE[0]], dtype=float) / math.hypot(*MODE)
+ANGLE_DEG = math.degrees(math.atan2(MODE[1], MODE[0]))
+RE = 100  # Re = U L / nu, the same nu as Taylor-Green Re100 at equal N
+FIT_FROM_TSTAR = 0.12  # late-time samples only; reduces startup sensitivity
+INITIALIZATION = (
+    "physical FP64 equilibrium G at actual cell centers; rho=1 (exact p const)"
+)
+# Per-row screens. phase_drift_abs, residual_l2 and the run-level rate screen
+# use eps(n) and are preregistered diagnostic choices, not error bounds.
+LIMITS = dict(velocity_l2=0.05, spurious_pressure_scaled=0.02, mass_relative=1e-5)
+# The comparative pressure budget max(2*control, 0.02) is never below the 0.02
+# absolute ceiling, so it is always redundant; raw control ratios are kept.
+FLOORS = dict(rate_error=0.002, residual_l2=0.002, spurious_pressure_scaled=0.02)
+BASELINE_FACTOR = 2.0
+
+
+def wavevector(length):
+    return 2 * np.pi / length * np.array(MODE, dtype=float)
+
+
+def decay_rate(length, viscosity):
+    """Exact amplitude decay rate per root step: nu |k|^2."""
+    k = wavevector(length)
+    return viscosity * float(k @ k)
+
+
+def eps(n):
+    """(|k| h_root)^2 with h_root = 1: dimensionless, coefficient one."""
+    k = wavevector(n)
+    return float(k @ k)
+
+
+def exact(x, y, t, length, viscosity, amplitude=U):
+    k = wavevector(length)
+    a = amplitude * np.exp(-decay_rate(length, viscosity) * t)
+    s = np.sin(k[0] * x + k[1] * y + PHI)
+    p = np.zeros_like(s)
+    return a * E[0] * s, a * E[1] * s, p  # kinematic pressure, constant
+
+
+def equilibrium(rho, ux, uy):
+    # Same second-order D2Q9 equilibrium as oracle.initial.
+    cu = CX[:, None, None] * ux + CY[:, None, None] * uy
+    return (
+        W[:, None, None]
+        * rho
+        * (1 + 3 * cu + 4.5 * cu * cu - 1.5 * (ux * ux + uy * uy))
+    )
+
+
+def initial(shape, level, length, viscosity):
+    x, y = centers(shape, level)
+    ux, uy, p = exact(x, y, 0.0, length, viscosity)
+    return equilibrium(1 + 3 * p, ux, uy)
+
+
+def _leaves(arrays, masks):
+    accum = []
+    for level, (g, mask) in enumerate(zip(arrays, masks)):
+        if g.shape != (9, *mask.shape):
+            raise ValueError("population/mask shape mismatch")
+        values = g[:, mask].astype(float)
+        if not np.isfinite(values).all():
+            raise FloatingPointError("nonfinite leaf populations")
+        rho = values.sum(axis=0)
+        if (rho <= 0).any():
+            raise FloatingPointError("nonpositive leaf density")
+        x, y = centers(mask.shape, level)
+        area = np.full(rho.size, 0.25**level)
+        lev = np.full(rho.size, level)
+        accum.append((rho, CX @ values / rho, CY @ values / rho, x[mask], y[mask], area, lev))
+    return [np.concatenate(v) for v in zip(*accum)]
+
+
+def metrics(arrays, masks, t, length, viscosity, initial_mass=None):
+    rho, ux, uy, x, y, w, lev = _leaves(arrays, masks)
+    volume = w.sum()
+    if not np.isclose(volume, length**2, rtol=0, atol=1e-10):
+        raise ValueError("leaf partition does not cover domain once")
+    k = wavevector(length)
+    theta = k[0] * x + k[1] * y + PHI
+    s, c = np.sin(theta), np.cos(theta)
+    v = ux * E[0] + uy * E[1]
+    # Weighted least squares on the actual composite quadrature (no
+    # continuous-orthogonality assumption).
+    gram = np.array([[w @ (s * s), w @ (s * c)], [w @ (s * c), w @ (c * c)]])
+    a_hat, b_hat = np.linalg.solve(gram, [w @ (v * s), w @ (v * c)])
+    ex, ey, _ = exact(x, y, t, length, viscosity)
+    exact_norm = w @ (ex * ex + ey * ey)
+    fit = a_hat * s + b_hat * c
+    longitudinal = (ux * k[0] + uy * k[1]) / math.sqrt(k @ k)
+    amplitude = U * math.exp(-decay_rate(length, viscosity) * t)
+    scale = amplitude * amplitude / 4  # continuum reference scale, not the TG norm
+    modal = math.hypot(a_hat, b_hat)
+    p = (rho - 1) / 3
+    p = p - w @ p / volume  # one global gauge
+    mass = float(w @ rho)
+    return dict(
+        tstar=t * U / length,
+        mass=mass,
+        mass_relative=(
+            0.0 if initial_mass is None else abs(mass - initial_mass) / initial_mass
+        ),
+        a_hat=float(a_hat),
+        b_hat=float(b_hat),
+        amplitude_exact=amplitude,
+        amplitude_ratio=float(a_hat / amplitude),
+        modal_amplitude=modal,
+        amplitude_error=modal / amplitude - 1,
+        phase_drift=float(math.atan2(b_hat, a_hat)),
+        phase_drift_abs=abs(float(math.atan2(b_hat, a_hat))),
+        velocity_l2=float(
+            np.sqrt(w @ ((ux - ex) ** 2 + (uy - ey) ** 2) / exact_norm)
+        ),
+        residual_l2=float(
+            np.sqrt(w @ ((ux - fit * E[0]) ** 2 + (uy - fit * E[1]) ** 2) / exact_norm)
+        ),
+        longitudinal_l2=float(np.sqrt(w @ (longitudinal**2) / exact_norm)),
+        spurious_pressure_rms=float(np.sqrt(w @ (p * p) / volume)),
+        spurious_pressure_scaled=float(np.sqrt(w @ (p * p) / volume) / scale),
+        spurious_pressure_max_scaled=float(np.max(np.abs(p)) / scale),
+        spurious_pressure_scaled_by_level=[
+            float(np.sqrt(w[lev == i] @ (p[lev == i] ** 2) / w[lev == i].sum()) / scale)
+            for i in range(len(masks))
+        ],
+        rho_min=float(rho.min()),
+        rho_max=float(rho.max()),
+        max_speed=float(np.max(np.hypot(ux, uy))),
+        leaf_cells=int(len(rho)),
+    )
+
+
+def fit_rate(rows, length, viscosity):
+    """Free-intercept least squares of ln R, R = hypot(a_hat, b_hat), late samples.
+
+    The intercept absorbs only a constant amplitude prefactor once a single
+    exponential dominates; per-interval rates and the fit residual are kept so
+    a lingering transient stays visible. Rates are per root step of the record
+    (uniform 2N uses its own nu and doubled step counts).
+    """
+    late = [r for r in rows if r["tstar"] >= FIT_FROM_TSTAR - 1e-12]
+    amp = np.array([math.hypot(r["a_hat"], r["b_hat"]) for r in late])
+    if len(late) < 3 or not np.isfinite(amp).all():
+        raise ValueError("need >= 3 late samples with finite R")
+    exact_rate = decay_rate(length, viscosity)
+    to_tstar = length / U  # root steps per unit t*
+    if not (amp > 0).all():
+        # Finite modal extinction: a physical failure, not malformed data. No
+        # clipping and no invented rate; fit quantities are unavailable.
+        return dict(rate_fit_valid=False, rate=None, exact_rate=exact_rate,
+                    rate_tstar=None, exact_rate_tstar=exact_rate * to_tstar,
+                    rate_error=None, interval_rate_errors=None, log_fit_residual_max=None)
+    t = np.array([r["step"] for r in late], dtype=float)
+    (slope, intercept) = np.polyfit(t, np.log(amp), 1)
+    residual = np.log(amp) - (slope * t + intercept)
+    return dict(
+        rate_fit_valid=True, rate=float(-slope), exact_rate=exact_rate,
+        rate_tstar=float(-slope * to_tstar),
+        exact_rate_tstar=exact_rate * to_tstar,  # = 5 (2 pi)^2 / Re
+        rate_error=abs(-slope / exact_rate - 1),
+        interval_rate_errors=[
+            float(-np.log(amp[i + 1] / amp[i]) / (t[i + 1] - t[i]) / exact_rate - 1)
+            for i in range(len(t) - 1)
+        ],
+        log_fit_residual_max=float(np.max(np.abs(residual))),
+    )
+
+
+NONNEGATIVE = ("tstar", "mass", "mass_relative", "amplitude_exact", "modal_amplitude",
+               "phase_drift_abs", "velocity_l2", "residual_l2", "longitudinal_l2",
+               "spurious_pressure_rms", "spurious_pressure_scaled",
+               "spurious_pressure_max_scaled", "rho_min", "rho_max", "max_speed")
+SIGNED = ("a_hat", "b_hat", "amplitude_ratio", "amplitude_error", "phase_drift")
+
+
+def _number(value):
+    return (isinstance(value, (int, float)) and not isinstance(value, bool)
+            and math.isfinite(value))
+
+
+def validate_rows(record):
+    """Exact schedule, finite numbers (incl. lists) and nonnegative error fields.
+
+    Signed quantities (amplitudes, amplitude/phase errors) keep their sign.
+    Raises ValueError so callers can report an invalid run before any screen.
+    """
+    n, rows = record["n"], record["rows"]
+    steps = [round(tstar * n / U) for tstar in TIMES]
+    if not isinstance(rows, list) or len(rows) != len(steps):
+        raise ValueError("rows do not match the scheduled samples")
+    for row, step in zip(rows, steps):
+        if type(row.get("step")) is not int or row["step"] != step:
+            raise ValueError(f"unscheduled or duplicate step {row.get('step')!r}")
+        if not _number(row.get("tstar")) or abs(row["tstar"] - step * U / n) > 3e-17:
+            raise ValueError(f"wrong sample time at step {step}")
+        if not all(_number(row.get(key)) for key in (*NONNEGATIVE, *SIGNED)):
+            raise ValueError(f"missing or nonfinite diagnostic at step {step}")
+        if any(row[key] < 0 for key in NONNEGATIVE) or row["rho_min"] <= 0:
+            raise ValueError(f"negative error quantity or density at step {step}")
+        levels = row.get("spurious_pressure_scaled_by_level")
+        if (not isinstance(levels, list) or len(levels) != record["depth"]
+                or not all(_number(v) and v >= 0 for v in levels)):
+            raise ValueError(f"invalid per-level pressure at step {step}")
+
+
+def summarize(record):
+    validate_rows(record)
+    rows = record["rows"]
+    n = record["n"]
+    out = fit_rate(rows, n, record["viscosity"])
+    out.update(
+        eps=eps(n),
+        complete=record["completed"] is True,
+        residual_l2=max(r["residual_l2"] for r in rows),
+        phase_drift_abs=max(r["phase_drift_abs"] for r in rows),
+        spurious_pressure_scaled=max(r["spurious_pressure_scaled"] for r in rows),
+        spurious_pressure_rms=max(r["spurious_pressure_rms"] for r in rows),
+        amplitude_error_abs=max(abs(r["amplitude_error"]) for r in rows),
+        velocity_l2=max(r["velocity_l2"] for r in rows),
+        mass_relative=max(r["mass_relative"] for r in rows),
+    )
+    if out["rate_fit_valid"] and not all(
+            _number(v) for v in (out["rate"], out["rate_error"], *out["interval_rate_errors"],
+                                 out["log_fit_residual_max"])):
+        raise ValueError("nonfinite rate fit")
+    out["own_screens_pass"] = out["complete"] and out["rate_fit_valid"] and all(
+        out[key] <= out["eps"] for key in ("rate_error", "phase_drift_abs", "residual_l2")
+    ) and all(out[key] <= limit for key, limit in LIMITS.items())
+    return out
+
+
+def _usable(record):
+    try:
+        out = summarize(record)
+    except (KeyError, TypeError, ValueError):
+        return None
+    return out if out["complete"] else None
+
+
+def assess(base_control, fine_control, refined):
+    """Uniform root-spacing and finest-spacing controls versus refined runs.
+
+    INVALID_CONTROL: missing, incomplete, invalid-row or mixed-source controls.
+    INVALID_RUN: a refined record has invalid rows (schedule, NaN, sign).
+    CONTROL_SCREEN_FAIL: complete valid controls that fail their own screens;
+    refined own-screen results are kept, but no controlled comparison passes.
+    PASS/FAIL: every refined run's controlled_pass (own screens and budgets).
+    """
+    if base_control is None or fine_control is None or not refined:
+        return dict(status="INVALID_CONTROL", reason="missing control or refined run")
+    n = base_control["n"]
+    if (base_control["depth"], fine_control["depth"], fine_control["n"]) != (1, 1, 2 * n):
+        raise ValueError("controls must be uniform N and uniform 2N")
+    if any(r["depth"] != 2 or r["n"] != n for r in refined):
+        raise ValueError("refined runs must be depth 2 at root N")
+    if len({r.get("explosion") for r in refined}) != len(refined):
+        raise ValueError("duplicate refined explosion scheme")
+    sources = {repr(r.get("source")) for r in [base_control, fine_control, *refined]}
+    base, fine = _usable(base_control), _usable(fine_control)
+    if base is None or fine is None or len(sources) != 1:
+        return dict(status="INVALID_CONTROL",
+                    reason="incomplete/invalid control or mixed source")
+    controls_pass = base["own_screens_pass"] and fine["own_screens_pass"]
+    report = dict(base_control=base, fine_control=fine, refined={},
+                  controls_pass=controls_pass)
+    statuses = []
+    for r in refined:
+        try:
+            s = summarize(r)
+        except (KeyError, TypeError, ValueError) as error:
+            report["refined"][r["explosion"]] = dict(valid=False, reason=str(error))
+            statuses.append("INVALID_RUN")
+            continue
+        s["valid"] = True
+        # rate_error is None when a modal fit is unavailable; never a pass.
+        s["budget"] = {
+            key: (None if base[key] is None else max(BASELINE_FACTOR * base[key], floor))
+            for key, floor in FLOORS.items()
+        }
+        s["budget_pass"] = all(
+            s[key] is not None and s["budget"][key] is not None and s[key] <= s["budget"][key]
+            for key in FLOORS)
+        s["control_ratios"] = {
+            name: {key: (s[key] / c[key] if s[key] is not None and c[key] else None)
+                   for key in (*FLOORS, "spurious_pressure_rms")}
+            for name, c in (("root", base), ("finest", fine))
+        }
+        s["controlled_pass"] = controls_pass and s["own_screens_pass"] and s["budget_pass"]
+        report["refined"][r["explosion"]] = s
+        statuses.append("PASS" if s["controlled_pass"] else "FAIL")
+    if "INVALID_RUN" in statuses:
+        report["status"] = "INVALID_RUN"
+    elif not controls_pass:
+        report["status"] = "CONTROL_SCREEN_FAIL"
+    else:
+        report["status"] = "PASS" if all(v == "PASS" for v in statuses) else "FAIL"
+    return report

tests/validation/taylor_green/test_projected_pairs.py

@@ -319,7 +319,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)
+            verify.assert_called_once_with(REPO, MOCK_REVISION, "taylor_green")
             self.assertTrue(json.loads(stdout.getvalue())["dry_run"])
             self.assertFalse(output.exists())
 
@@ -412,7 +412,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)
+                verify.assert_called_once_with(REPO, MOCK_REVISION, "taylor_green")
                 self.assertEqual(selectors, ["projected_linear" if early else "constant"] if depth == 2 else [])
                 self.assertEqual(receipt["source"], source())
                 self.assertEqual(receipt["runtime_identity"], RUNTIME)

tests/validation/taylor_green/test_shear_wave.py

@@ -0,0 +1,471 @@
+"""Standalone CPU manufactured checks for the oblique shear wave.
+
+Run from this directory with PYTHONPATH=<checkout>/src (same as test_oracle.py).
+"""
+
+import hashlib
+import json
+import math
+import unittest
+from pathlib import Path
+
+import numpy as np
+import oracle
+import run_case
+import shear_wave as sw
+from oracle import CX, CY, U, TIMES, centers, rectangles
+
+N = 32
+NU = U * N / sw.RE
+HERE = Path(__file__).resolve().parent
+ORACLE_SHA256_EEC4677 = "83e4d7425d45538278e9fb181af450a36358bd9a5ed6be57ff611dd9ae5eb35b"
+
+
+def composite_masks(n):
+    """Level-wide union masks for rectangles(n, 2): root complement + fine box."""
+    (x0, y0, x1, y1), = rectangles(n, 2)[0]
+    root = np.ones((n, n), dtype=bool)
+    root[y0:y1, x0:x1] = False
+    fine = np.zeros((2 * n, 2 * n), dtype=bool)
+    fine[2 * y0:2 * y1, 2 * x0:2 * x1] = True
+    return [root, fine]
+
+
+def populations(masks, field, t=0.0):
+    """Equilibrium G of an arbitrary (ux, uy, rho) field on every level."""
+    out = []
+    for level, m in enumerate(masks):
+        x, y = centers(m.shape, level)
+        ux, uy, rho = field(x, y, t)
+        out.append(sw.equilibrium(rho, ux, uy))
+    return out
+
+
+def exact_field(phase=0.0, rate_factor=1.0, mode=sw.MODE, longitudinal=0.0, ripple=0.0,
+                offset=0.0):
+    k = 2 * np.pi / N * np.array(mode, dtype=float)
+    e = np.array([-mode[1], mode[0]], dtype=float) / math.hypot(*mode)
+    rate = rate_factor * sw.decay_rate(N, NU)
+    kh = k / np.linalg.norm(k)
+
+    def field(x, y, t):
+        theta = k[0] * x + k[1] * y + phase
+        a = U * np.exp(-rate * t)
+        ux = a * e[0] * np.sin(theta) + longitudinal * kh[0] * np.sin(theta)
+        uy = a * e[1] * np.sin(theta) + longitudinal * kh[1] * np.sin(theta)
+        rho = 1 + offset + ripple * np.cos(theta)
+        return ux, uy, rho
+
+    return field
+
+
+class ConventionTests(unittest.TestCase):
+    def test_mode_angle_and_polarization(self):
+        k = sw.wavevector(N)
+        self.assertAlmostEqual(sw.ANGLE_DEG, 26.565051177077994)
+        self.assertAlmostEqual(math.degrees(math.atan2(k[1], k[0])), sw.ANGLE_DEG)
+        self.assertAlmostEqual(float(sw.E @ k), 0.0)
+        self.assertAlmostEqual(float(sw.E @ sw.E), 1.0)
+        self.assertAlmostEqual(sw.eps(128), 5 * (2 * math.pi / 128) ** 2)
+
+    def test_navier_stokes_and_vanishing_advection(self):
+        # Independent central differences: du/dt = nu lap u, (u.grad)u = 0, div u = 0.
+        x = np.array([2.7, 8.1, 13.3])
+        y = np.array([3.2, 7.4, 21.6])
+        t, h = 10.0, 1e-3
+        f = lambda x, y, t: np.array(sw.exact(x, y, t, N, NU)[:2])
+        z = f(x, y, t)
+        dt = (f(x, y, t + h) - f(x, y, t - h)) / (2 * h)
+        dx = (f(x + h, y, t) - f(x - h, y, t)) / (2 * h)
+        dy = (f(x, y + h, t) - f(x, y - h, t)) / (2 * h)
+        lap = (f(x + h, y, t) + f(x - h, y, t) + f(x, y + h, t) + f(x, y - h, t) - 4 * z) / h**2
+        np.testing.assert_allclose(dt, NU * lap, atol=1e-9)
+        # Analytic identities: div u = A cos(theta) (E.k) and (u.grad)u = A sin(theta) (E.k) u/A.
+        k = sw.wavevector(N)
+        self.assertLess(abs(float(sw.E @ k)), 1e-17)
+        # Central differences do not cancel exactly: each derivative carries
+        # -(k_j h)^2/6 relative truncation, so div_h = -A cos(theta) h^2/6 sum_j E_j k_j^3
+        # + O(h^4) (nonzero for an oblique k). Bound it with that leading term plus roundoff.
+        a = U * math.exp(-sw.decay_rate(N, NU) * t)
+        trunc = a * h * h / 6 * float(np.abs(sw.E * k**3).sum()) * 1.01
+        # Argument rounding: eps*max|x|/h*U*|k| ~ 2.2e-16*21.6/1e-3*0.03*0.44 ~ 6e-14.
+        roundoff = 1e-13
+        self.assertLess(np.max(np.abs(dx[0] + dy[1])), trunc + roundoff)
+        # Discrete advection: same truncation scaled by |u| <= a.
+        self.assertLess(np.max(np.abs(z[0] * dx + z[1] * dy)), a * trunc + roundoff)
+        # Halving h reduces the divergence truncation about fourfold.
+        hh = h / 2
+        div2 = ((f(x + hh, y, t) - f(x - hh, y, t))[0] + (f(x, y + hh, t) - f(x, y - hh, t))[1]) / (2 * hh)
+        np.testing.assert_allclose(div2 * 4, dx[0] + dy[1], rtol=1e-2, atol=6 * roundoff)
+        self.assertTrue(np.all(sw.exact(x, y, t, N, NU)[2] == 0))
+
+
+class InitialTests(unittest.TestCase):
+    def test_moments_on_each_level(self):
+        for level, m in enumerate(composite_masks(N)):
+            g = sw.initial(m.shape, level, N, NU)
+            x, y = centers(m.shape, level)
+            ux, uy, _ = sw.exact(x, y, 0.0, N, NU)
+            rho = g.sum(axis=0)
+            np.testing.assert_allclose(rho, 1.0, atol=1e-14)
+            np.testing.assert_allclose(np.tensordot(CX, g, 1) / rho, ux, atol=1e-15)
+            np.testing.assert_allclose(np.tensordot(CY, g, 1) / rho, uy, atol=1e-15)
+
+    def test_equilibrium_matches_frozen_oracle_formula(self):
+        # Same equilibrium as oracle.initial when fed the TG field.
+        shape = (8, 8)
+        x, y = centers(shape, 1)
+        ux, uy, p = oracle.exact(x, y, 0.0, 4.0, 0.01)
+        np.testing.assert_array_equal(
+            sw.equilibrium(1 + 3 * p, ux, uy), oracle.initial(shape, 1, 4.0, 0.01)
+        )
+
+
+class MetricTests(unittest.TestCase):
+    def setUp(self):
+        self.masks = composite_masks(N)
+
+    def sample(self, field, t):
+        return sw.metrics(populations(self.masks, field, t), self.masks, t, N, NU, None)
+
+    def test_exact_field_on_composite_partition(self):
+        for t in (0.0, 100.0, 512.0):
+            r = self.sample(exact_field(), t)
+            self.assertLess(r["velocity_l2"], 1e-12)
+            self.assertLess(r["residual_l2"], 1e-12)
+            self.assertLess(r["longitudinal_l2"], 1e-12)
+            self.assertLess(r["phase_drift_abs"], 1e-12)
+            self.assertAlmostEqual(r["amplitude_ratio"], 1.0, places=12)
+            self.assertLess(abs(r["amplitude_error"]), 1e-12)
+            self.assertLess(r["spurious_pressure_scaled"], 1e-10)
+            self.assertEqual(len(r["spurious_pressure_scaled_by_level"]), 2)
+        self.assertEqual(r["leaf_cells"], int(sum(m.sum() for m in self.masks)))
+
+    def test_phase_control(self):
+        r = self.sample(exact_field(phase=0.01), 0.0)
+        self.assertAlmostEqual(r["phase_drift"], 0.01, delta=1e-9)
+        self.assertLess(abs(r["amplitude_error"]), 1e-12)  # R = hypot(a, b) is phase-free
+        self.assertLess(r["residual_l2"], 1e-12)  # a pure phase shift stays modal
+
+    def test_rate_control_and_fit(self):
+        for factor, expected in ((1.0, 0.0), (1.01, 0.01)):
+            field = exact_field(rate_factor=factor)
+            rows = []
+            for tstar in TIMES:
+                step = round(tstar * N / U)
+                rows.append(dict(self.sample(field, step), step=step))
+            fit = sw.fit_rate(rows, N, NU)
+            self.assertAlmostEqual(fit["rate_error"], expected, delta=1e-9)
+            np.testing.assert_allclose(fit["interval_rate_errors"], expected, atol=1e-9)
+            self.assertLess(fit["log_fit_residual_max"], 1e-9)
+            self.assertAlmostEqual(fit["exact_rate_tstar"], 5 * (2 * math.pi) ** 2 / sw.RE)
+            self.assertAlmostEqual(fit["exact_rate_tstar"], 1.97392088, places=8)
+
+    def test_fit_ignores_startup_sample(self):
+        rows = [dict(tstar=0.0, step=0, a_hat=2 * U, b_hat=0.0)]
+        rows += [dict(tstar=s, step=round(s * N / U), b_hat=0.0,
+                      a_hat=U * math.exp(-sw.decay_rate(N, NU) * round(s * N / U)))
+                 for s in TIMES[1:]]
+        self.assertLess(sw.fit_rate(rows, N, NU)["rate_error"], 1e-12)
+        rows[-1]["a_hat"] = float("nan")
+        with self.assertRaises(ValueError):
+            sw.fit_rate(rows, N, NU)
+
+    def test_wrong_mode_control(self):
+        # Swapping (2,1) for (1,2), i.e. 26.6 vs 63.4 deg, is not silently modal.
+        # Composite quadrature need not give an exactly zero projection.
+        r = self.sample(exact_field(mode=(1, 2)), 0.0)
+        self.assertLess(abs(r["amplitude_ratio"]), 0.05)
+        self.assertGreater(r["residual_l2"], 0.95)
+
+    def test_longitudinal_control(self):
+        r = self.sample(exact_field(longitudinal=1e-3), 0.0)
+        self.assertAlmostEqual(r["longitudinal_l2"], 1e-3 / U, delta=1e-9)
+        # The full-vector residual and velocity error see it too.
+        self.assertAlmostEqual(r["residual_l2"], 1e-3 / U, delta=1e-9)
+        self.assertAlmostEqual(r["velocity_l2"], 1e-3 / U, delta=1e-9)
+        self.assertLess(r["phase_drift_abs"], 1e-12)
+        self.assertAlmostEqual(r["amplitude_ratio"], 1.0, places=12)
+
+    def test_pressure_gauge_and_ripple(self):
+        r = self.sample(exact_field(offset=1e-4), 0.0)
+        self.assertLess(r["spurious_pressure_scaled"], 1e-10)  # one global gauge
+        ripple = 3e-5
+        r = self.sample(exact_field(ripple=ripple), 0.0)
+        # p = ripple/3 cos(theta); rms over a whole number of periods = ripple/(3 sqrt 2).
+        # Composite midpoint quadrature is not exactly periodic: 3% tolerance.
+        expected = ripple / (3 * math.sqrt(2)) / (U * U / 4)
+        self.assertAlmostEqual(r["spurious_pressure_scaled"], expected, delta=0.03 * expected)
+        self.assertAlmostEqual(r["spurious_pressure_rms"], expected * U * U / 4,
+                               delta=0.03 * expected * U * U / 4)
+
+    def test_volume_weighting_and_double_count(self):
+        bad = [np.ones_like(m) for m in self.masks]
+        with self.assertRaises(ValueError):
+            sw.metrics(populations(bad, exact_field()), bad, 0.0, N, NU)
+        # Unweighted (count-based) averaging would change the fit; weights matter.
+        r = self.sample(exact_field(), 0.0)
+        self.assertEqual(r["leaf_cells"], N * N - self.masks[1].sum() // 4 + self.masks[1].sum())
+
+
+def record(n, depth, rows, explosion="constant", completed=True, source="s"):
+    for row in rows:
+        row["spurious_pressure_scaled_by_level"] = [row["spurious_pressure_scaled"]] * depth
+    return dict(n=n, depth=depth, viscosity=U * n / sw.RE, rows=rows, completed=completed,
+                explosion=explosion, source=source)
+
+
+def rows_for(n, rate_error=0.0, residual=1e-4, pressure=1e-3, phase=0.0):
+    rate = (1 + rate_error) * sw.decay_rate(n, U * n / sw.RE)
+    out = []
+    for tstar in TIMES:
+        step = round(tstar * n / U)
+        amp = U * math.exp(-rate * step)
+        out.append(dict(tstar=step * U / n, step=step, a_hat=amp, b_hat=0.0,
+                        amplitude_ratio=1.0, modal_amplitude=amp, amplitude_exact=amp,
+                        phase_drift=phase, mass=float(n * n), rho_min=0.999, rho_max=1.001,
+                        max_speed=amp, longitudinal_l2=0.0, spurious_pressure_max_scaled=pressure,
+                        amplitude_error=0.0, residual_l2=residual, phase_drift_abs=phase,
+                        spurious_pressure_scaled=pressure,
+                        spurious_pressure_rms=pressure * U * U / 4,
+                        velocity_l2=residual, mass_relative=0.0))
+    return out
+
+
+class AssessTests(unittest.TestCase):
+    def test_budget_form_is_max_not_sum(self):
+        base = record(128, 1, rows_for(128, rate_error=1e-3, pressure=0.015))
+        fine = record(256, 1, rows_for(256, rate_error=2e-4))
+        # rate budget = max(2*1e-3, 0.002) = 0.002; 0.0025 would pass a sum form.
+        refined = record(128, 2, rows_for(128, rate_error=0.0025))
+        report = sw.assess(base, fine, [refined])
+        self.assertEqual(report["status"], "FAIL")
+        self.assertAlmostEqual(report["refined"]["constant"]["budget"]["rate_error"], 0.002)
+        # pressure budget = max(2*0.015, 0.02) = 0.03 (control-dominated).
+        self.assertAlmostEqual(
+            report["refined"]["constant"]["budget"]["spurious_pressure_scaled"], 0.03)
+        self.assertIn("finest", report["refined"]["constant"]["control_ratios"])
+        ok = record(128, 2, rows_for(128, rate_error=0.0015), explosion="projected_linear")
+        self.assertEqual(sw.assess(base, fine, [ok])["status"], "PASS")
+
+    def test_control_outcomes_are_distinct(self):
+        fine = record(256, 1, rows_for(256))
+        refined = record(128, 2, rows_for(128))
+        failing = record(128, 1, rows_for(128, residual=0.5))  # complete, > eps(128)
+        report = sw.assess(failing, fine, [refined])
+        self.assertEqual(report["status"], "CONTROL_SCREEN_FAIL")
+        own = report["refined"]["constant"]
+        self.assertTrue(own["own_screens_pass"])  # refined own screens kept
+        self.assertFalse(own["controlled_pass"])  # but no controlled comparison passes
+        good = record(128, 1, rows_for(128))
+        self.assertEqual(sw.assess(None, fine, [refined])["status"], "INVALID_CONTROL")
+        incomplete = record(128, 1, rows_for(128), completed=False)
+        self.assertEqual(sw.assess(incomplete, fine, [refined])["status"], "INVALID_CONTROL")
+        mixed = record(256, 1, rows_for(256), source="other")
+        self.assertEqual(sw.assess(good, mixed, [refined])["status"], "INVALID_CONTROL")
+
+    def test_invalid_rows_are_rejected_before_screens(self):
+        base = record(128, 1, rows_for(128))
+        fine = record(256, 1, rows_for(256))
+        for mutate in (
+            lambda rows: rows[2].update(rho_min=float("nan")),
+            lambda rows: rows[2].update(spurious_pressure_scaled_by_level=[0.0, float("nan")]),
+            lambda rows: rows[2].update(step=rows[1]["step"], tstar=rows[1]["tstar"]),
+            lambda rows: rows[2].update(residual_l2=-1e-6),
+            lambda rows: rows[2].pop("longitudinal_l2"),
+        ):
+            refined = record(128, 2, rows_for(128))
+            mutate(refined["rows"])
+            self.assertEqual(sw.assess(base, fine, [refined])["status"], "INVALID_RUN")
+        # Legitimate signs stay allowed: negative amplitude and phase errors.
+        refined = record(128, 2, rows_for(128))
+        for row in refined["rows"]:
+            row.update(amplitude_error=-1e-3, phase_drift=-1e-4, phase_drift_abs=1e-4)
+        report = sw.assess(base, fine, [refined])
+        self.assertEqual(report["status"], "PASS")
+        self.assertEqual(sw.assess(base, fine, [])["status"], "INVALID_CONTROL")
+        with self.assertRaises(ValueError):
+            sw.assess(base, fine, [record(128, 2, rows_for(128)), record(128, 2, rows_for(128))])
+
+    def test_modal_extinction_is_a_failure_not_malformed(self):
+        def extinct(rows):
+            for row in rows[1:]:
+                row.update(a_hat=0.0, b_hat=0.0, modal_amplitude=0.0, amplitude_ratio=0.0,
+                           amplitude_error=-1.0)
+            return rows
+        base = record(128, 1, rows_for(128))
+        fine = record(256, 1, rows_for(256))
+        refined = record(128, 2, extinct(rows_for(128)))
+        report = sw.assess(base, fine, [refined])
+        self.assertEqual(report["status"], "FAIL")
+        own = report["refined"]["constant"]
+        self.assertTrue(own["valid"])
+        self.assertFalse(own["rate_fit_valid"])
+        self.assertIsNone(own["rate_error"])
+        self.assertIsNone(own["control_ratios"]["root"]["rate_error"])
+        self.assertFalse(own["own_screens_pass"] or own["controlled_pass"] or own["budget_pass"])
+        self.assertEqual(own["residual_l2"], 1e-4)  # other finite metrics retained
+        json.loads(json.dumps(report, allow_nan=False))
+        control = sw.assess(record(128, 1, extinct(rows_for(128))), fine,
+                            [record(128, 2, rows_for(128))])
+        self.assertEqual(control["status"], "CONTROL_SCREEN_FAIL")
+        self.assertIsNone(control["refined"]["constant"]["budget"]["rate_error"])
+        self.assertFalse(control["refined"]["constant"]["controlled_pass"])
+        json.loads(json.dumps(control, allow_nan=False))
+
+    def test_wrong_shapes(self):
+        base = record(128, 1, rows_for(128))
+        fine = record(256, 1, rows_for(256))
+        refined = record(128, 2, rows_for(128))
+        with self.assertRaises(ValueError):
+            sw.assess(base, record(128, 1, rows_for(128)), [refined])
+        with self.assertRaises(ValueError):
+            sw.assess(record(128, 1, rows_for(128)), fine, [record(256, 2, rows_for(256))])
+
+    def test_incomplete_record_fails(self):
+        base = record(128, 1, rows_for(128))
+        fine = record(256, 1, rows_for(256))
+        refined = record(128, 2, rows_for(128), completed=False)
+        self.assertEqual(sw.assess(base, fine, [refined])["status"], "FAIL")
+
+
+class FakeRun:
+    def __init__(self, arrays):
+        self.arrays = arrays
+
+    def download(self):
+        return self.arrays if len(self.arrays) > 1 else self.arrays[0]
+
+
+class RunnerSampleTests(unittest.TestCase):
+    def test_refined_shear_sample_with_list_diagnostics(self):
+        masks = composite_masks(N)
+        arrays = [sw.initial(m.shape, level, N, NU) for level, m in enumerate(masks)]
+        row = run_case.sample(FakeRun(arrays), masks, 0, N, NU, None, sw.metrics, 2)
+        self.assertIsInstance(row["spurious_pressure_scaled_by_level"], list)
+        self.assertTrue(run_case.finite_row(row))
+        bad = [a.copy() for a in arrays]
+        bad[1][0, masks[1]] = np.nan
+        with self.assertRaises(FloatingPointError):
+            run_case.sample(FakeRun(bad), masks, 0, N, NU, None, sw.metrics, 2, "READY")
+        self.assertFalse(run_case.finite_row(dict(a=1.0, b=[0.0, float("inf")])))
+
+    def test_taylor_green_sample_unchanged(self):
+        masks = composite_masks(N)
+        arrays = [oracle.initial(m.shape, level, N, NU) for level, m in enumerate(masks)]
+        row = run_case.sample(FakeRun(arrays), masks, 0, N, NU, None, oracle.metrics, 2)
+        self.assertEqual(row, oracle.metrics(arrays, masks, 0, N, NU, None))
+        uniform = [oracle.initial((N, N), 0, N, NU)]
+        row = run_case.sample(FakeRun(uniform), [np.ones((N, N), bool)], 0, N, NU, None,
+                              oracle.metrics, 1)
+        self.assertEqual(row["leaf_cells"], N * N)
+
+
+def receipt(n, depth, explosion="constant"):
+    """Synthetic receipt satisfying every contract field of the fixed four-case set."""
+    import assess_projected_pairs as tg
+    steps = run_case.sample_steps(n)[0]
+    rows = rows_for(n)
+    for row, step in zip(rows, steps):
+        row.update(step=step, tstar=step * U / n,
+                   spurious_pressure_scaled_by_level=[row["spurious_pressure_scaled"]] * depth)
+    cells = [n * n] if depth == 1 else [8704, 30720]
+    hops = [16 * n * 2**level for level in range(depth)]
+    builds = []
+    for level in range(depth):
+        size = n * 2**level
+        builds.append(dict(source_sha256=tg.BUILD_SOURCES[size], arch="sm_70",
+                           flags=["--use_fast_math"], n_objects=0, ptx_sha256="a" * 64,
+                           config=dict(nx=10, ny=10, nz=1, lattice_model="D2Q9", collision="SRT",
+                                       store_precision="FP32", streaming="double_buffer",
+                                       ddf_shifting=True, les_enabled=False, velocity=0.0,
+                                       viscosity=U * size / 100, rho=1.0,
+                                       omega=1 / (.5 + 3 * U * size / 100),
+                                       threads_per_block=256)))
+    checkout = "/x/CelerisLab"
+    return dict(
+        n=n, depth=depth, re_L=100, explosion=explosion, U=U, mach=U * math.sqrt(3),
+        viscosity=U * n / 100, tau_by_level=[0.5 + 3 * U * n / 100 * 2**l for l in range(depth)],
+        source=dict(revision="r", dirty=False, source_tree_sha256=run_case.SOURCE_SHA256,
+                    harness_sha256={"x.py": "h"}, build_capture_sha256=run_case.BUILD_CAPTURE_SHA256,
+                    shear_wave_sha256=run_case.SHEAR_WAVE_SHA256),
+        runtime_identity=dict(device=dict(name="V100", pci_bus_id="0", capability=[7, 0],
+                                          driver=1, runtime=[12, 0]),
+                              nvcc_path="/nvcc", nvcc_version="v", python="3", numpy="1",
+                              pycuda="2"),
+        runtime_source=dict(checkout=checkout, module=checkout + "/src/CelerisLab/__init__.py",
+                            compiler=checkout + "/src/CelerisLab/cuda/compiler_v2.py",
+                            kernel_root=checkout + "/src/CelerisLab/lbm/kernels"),
+        cuda_visible_devices="0", numerical_revision=run_case.NUMERICAL_REVISION,
+        reference_harness=run_case.REFERENCE_HARNESS, numerics=dict(tg.NUMERICS),
+        initialization=sw.INITIALIZATION, limits=dict(sw.LIMITS),
+        early_pressure=dict(enabled=False, rows=[], sampled_max_pressure_l2=None),
+        refinements=[] if depth == 1 else rectangles(n, depth),
+        workload=dict(active_cells_by_level=cells, level_hops=hops, root_intervals=16 * n,
+                      active_collision_updates_by_level=[c * h for c, h in zip(cells, hops)],
+                      active_collision_updates=sum(c * h for c, h in zip(cells, hops))),
+        actual_build_sidecars=builds, rows=rows, completed=True,
+        flow="shear_wave", mode=list(sw.MODE), polarization=sw.E.tolist(),
+        angle_deg=sw.ANGLE_DEG, phase=sw.PHI, eps=sw.eps(n))
+
+
+class ContractTests(unittest.TestCase):
+    def setUp(self):
+        import assess_shear_wave
+        self.assess = assess_shear_wave.assess
+        self.records = [receipt(128, 1), receipt(256, 1), receipt(128, 2),
+                        receipt(128, 2, "projected_linear")]
+
+    def test_complete_contract_passes(self):
+        result = self.assess(self.records, "r")
+        self.assertTrue(result["contract_valid"], [k for k, v in result["checks"].items() if not v])
+        self.assertEqual(result["status"], "PASS")
+
+    def test_reproduced_false_passes_are_rejected(self):
+        a = self.assess
+        self.assertEqual(a(self.records[:2], "r")["status"], "INVALID_SET")  # empty refined
+        dup = [*self.records[:3], receipt(128, 2)]
+        self.assertEqual(a(dup, "r")["status"], "INVALID_SET")  # duplicate constant-only
+        for mutate in (
+            lambda r: r.pop("source"),
+            lambda r: r.update(flow="taylor_green"),
+            lambda r: r.update(mode=[1, 2]),
+            lambda r: r.update(re_L=1000),
+            lambda r: r.update(refinements=[[(0, 0, 8, 8)]]),
+            lambda r: r["actual_build_sidecars"][0].update(arch="sm_80"),
+            lambda r: r["rows"][0].update(step=1, tstar=U / 128),  # wrong initial time
+            lambda r: r["source"].update(shear_wave_sha256="0" * 64),
+        ):
+            records = [dict(r, rows=[dict(x) for x in r["rows"]],
+                            source=dict(r["source"]),
+                            actual_build_sidecars=[dict(b) for b in r["actual_build_sidecars"]])
+                       for r in self.records]
+            mutate(records[2])
+            self.assertEqual(a(records, "r")["status"], "INVALID_SET")
+        mixed = [dict(r) for r in self.records]
+        mixed[1] = dict(mixed[1], source=dict(mixed[1]["source"], revision="other"))
+        self.assertEqual(a(mixed, "r")["status"], "INVALID_SET")
+        self.assertEqual(a(self.records, "wrong-revision")["status"], "INVALID_SET")
+
+
+class DefaultPathTests(unittest.TestCase):
+    def test_frozen_oracle_is_byte_identical(self):
+        digest = hashlib.sha256((HERE / "oracle.py").read_bytes()).hexdigest()
+        self.assertEqual(digest, ORACLE_SHA256_EEC4677)
+
+    def test_default_flow_is_frozen_taylor_green(self):
+        init, met, limits, text = run_case.flow_functions()
+        self.assertIs(init, oracle.initial)
+        self.assertIs(met, oracle.metrics)
+        self.assertIs(limits, oracle.LIMITS)
+        self.assertEqual(text, run_case.INITIALIZATION)
+        self.assertEqual(run_case.flow_functions("taylor_green"), (init, met, limits, text))
+        shear = run_case.flow_functions("shear_wave")
+        self.assertIs(shear[0], sw.initial)
+        self.assertIs(shear[1], sw.metrics)
+        with self.assertRaises(ValueError):
+            run_case.flow_functions("tg")
+
+
+if __name__ == "__main__":
+    unittest.main()