test: signed interface-probe trace diagnostic (--probe-trace)

Test-only mechanism diagnostic, no pass gates. New pinned shear_probes.py:
24 fixed root-physical points on vertical/horizontal interface-normal
transects (>=8 root cells from corners), ACTIVE-leaf selection storing
level/i/j/actual centers, root-sync rho/ux/uy at root steps 0-128 and
1920-2048 (U512 every 4 lattice steps), analytic residuals at actual
centers, schedule/finite/600 kB validation. run_case: --probe-trace for
nested U128/U512/projected d3 and twin projected d2 (not with
--final-snapshot), shear_probes.py pin, probe-only steps skip metrics,
five qualification samples unchanged. shear_wave/shear_nested/shear_twin/
oracle byte-identical.

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

docs/validation/taylor_green.md

@@ -418,3 +418,39 @@ Labels:
 - Neither contrast is additive by construction.
 - The thin-band/acoustic explanation for the nested L1 pressure is a
   hypothesis, not established.
+
+### Signed interface-probe traces (opt-in `--probe-trace`, mechanism diagnostic)
+
+`shear_probes.py` (pinned) defines 24 fixed points on two interface-normal
+transects, all at least 8 root cells from any box corner:
+- V: y = 64.1, crossing the vertical faces x = 16 (L0/L1) and x = 22 (L1/L2);
+- H: x = 90.1, crossing the horizontal faces y = 24 and y = 30;
+- offsets from the outer face: −3.9, −1.9, −0.4, 0.1, 0.6, 2.1, 4.1, 5.6, 6.1,
+  6.6, 8.1, 10.1.
+
+Allowed cases:
+- nested U128, nested U512, nested projected depth 3;
+- twin projected depth 2.
+
+The option can't be combined with `--final-snapshot`.
+
+Sampling:
+- Each point is read from the ACTIVE leaf cell containing it. The receipt
+  stores the leaf level, `i`/`j`, and the actual cell center, all in
+  root-physical units (N=128).
+- Grids differ between cases, so centers are *not* collocated. Analytic
+  fields are evaluated at each actual center (`shear_probes.residuals`), and
+  comparisons across cases are between residuals at nearby but distinct
+  centers.
+- Times are root-physical steps 0–128 and 1920–2048. U512 is sampled every 4
+  lattice steps to match.
+- Each sample stores `rho`, `ux`, `uy` from the FP64 physical populations,
+  POST-collision G at root synchronization only. No substep content is
+  observed, and grid-scale alternation can alias.
+
+What's kept and enforced:
+- The five qualification samples stay unchanged.
+- No field dumps. The trace is validated for exact schedule, finiteness and
+  a 600 kB cap. If invalid it is dropped and an error is recorded, keeping
+  the scalar receipt.
+- There are no new pass gates; this is a mechanism diagnostic only.

tests/validation/taylor_green/run_case.py

@@ -25,6 +25,7 @@ NESTED_SHEAR = "nested_shear_n128"  # fixed three-level/two-region shear preset
 SHEAR_NESTED_SHA256 = "fa6cfe87fb4df8c6d732c637c885276833c81973b950c7724d734f8fc350bcb0"
 TWIN_SHEAR = "twin_shear_n128"  # depth-2 twin: nested L0->L1 boxes without L2
 SHEAR_TWIN_SHA256 = "1cbf37f18e89aea577581631c4565f51e68d22a3739b698816ef64acbf2b2bb4"
+SHEAR_PROBES_SHA256 = "98ea8090ace48f04b83dfb015db10e681e414e57eb748b34c0fbb90fb790d13b"  # signed interface-probe mechanism diagnostic
 
 
 def flow_functions(flow="taylor_green"):
@@ -94,7 +95,7 @@ def provenance(repo):
     )
 
 
-def verify_source(repo, revision, flow="taylor_green"):
+def verify_source(repo, revision, flow="taylor_green", probes=False):
     """Exact clean harness commit, fixed production tree, and reviewed capture helper."""
     import oracle
     from tests.validation.multigrid import check_projected_pairs as gate
@@ -113,6 +114,13 @@ def verify_source(repo, revision, flow="taylor_green"):
     if result["source_tree_sha256"] != SOURCE_SHA256:
         raise RuntimeError("Wrong numerical source fingerprint")
     result["build_capture_sha256"] = BUILD_CAPTURE_SHA256
+    if probes:
+        import shear_probes
+        module = harness / "shear_probes.py"
+        if (Path(shear_probes.__file__).resolve() != module
+                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_twin":
         import shear_twin
         module = harness / "shear_twin.py"
@@ -206,6 +214,9 @@ def main():
                         "(scaled_n256: Taylor-Green; scaled_shear_n256: shear wave) or the "
                         "fixed three-level/two-region shear preset (nested_shear_n128) "
                         "or its depth-2 twin (twin_shear_n128)")
+    p.add_argument("--probe-trace", action="store_true",
+                   help="signed interface-probe trace (nested U128/U512/projected d3, "
+                        "twin projected d2); mechanism diagnostic, no gates")
     p.add_argument("--final-snapshot", action="store_true",
                    help="nested_shear_n128 projected depth-3 only: final macro snapshot")
     p.add_argument("--early-pressure", action="store_true",
@@ -216,6 +227,13 @@ def main():
     if args.flow == "shear_wave" and (args.re != 100 or (args.depth == 3 and not nested)
                                       or args.early_pressure):
         p.error("shear_wave requires --re 100, depth 1 or 2 and no --early-pressure")
+    probe_cases = {(NESTED_SHEAR, 128, 1, "constant"), (NESTED_SHEAR, 512, 1, "constant"),
+                   (NESTED_SHEAR, 128, 3, "projected_linear"),
+                   (TWIN_SHEAR, 128, 2, "projected_linear")}
+    if args.probe_trace and (args.final_snapshot or (
+            args.diagnostic, args.n, args.depth, args.explosion) not in probe_cases):
+        p.error("--probe-trace requires nested U128/U512/projected depth 3 or twin projected "
+                "depth 2, without --final-snapshot")
     if args.final_snapshot and not (nested and args.depth == 3
                                     and args.explosion == "projected_linear"):
         p.error("--final-snapshot requires nested_shear_n128 projected depth 3")
@@ -253,7 +271,8 @@ def main():
     repo = args.repo.resolve()
     sys.path[:0] = [str(repo / "src"), str(repo)]
     source = verify_source(repo, args.expected_revision,
-                           "shear_nested" if nested else "shear_twin" if twin else args.flow)
+                           "shear_nested" if nested else "shear_twin" if twin else args.flow,
+                           probes=args.probe_trace)
     runtime = runtime_source(repo)
     if args.dry_run:
         print(json.dumps(dict(source=source, runtime_source=runtime, dry_run=True,
@@ -351,8 +370,23 @@ def main():
                 ]
             )
             initial_mass = None
-            for step in sorted(set(main_steps + early_steps)):
+            probe_steps = set()
+            if args.probe_trace:
+                import shear_probes
+                probe_root = dict(zip(shear_probes.case_steps(args.n), shear_probes.STEPS))
+                probe_steps = set(probe_root)
+                probes = shear_probes.select(masks, args.n)
+                report["probe_trace"] = shear_probes.new_trace(
+                    probes, cfg.viscosity, [args.diagnostic, args.n, args.depth, args.explosion])
+            for step in sorted(set(main_steps + early_steps) | probe_steps):
                 run.step(step - run.step_count)
+                if step in probe_steps:
+                    trace = report["probe_trace"]
+                    trace["steps"].append(probe_root[step])  # root-physical step
+                    trace["values"].append(shear_probes.sample(
+                        (run.download(),) if args.depth == 1 else run.download(), probes))
+                    if step not in main_steps and step not in early_steps:
+                        continue  # probe-only step: no qualification metrics
                 row = sample(run, masks, step, args.n, cfg.viscosity, initial_mass,
                              metrics, args.depth)
                 row["step"] = step
@@ -405,6 +439,11 @@ def main():
                         ),
                         flush=True,
                     )
+    if args.probe_trace:
+        try:
+            report["probe_trace"]["bytes"] = shear_probes.validate(report["probe_trace"])
+        except ValueError as error:  # keep the scalar receipt; drop an invalid trace
+            report["probe_trace"] = dict(error=str(error))
     report.update(
         completed=True,
         seconds=time.monotonic() - started,

tests/validation/taylor_green/shear_probes.py

@@ -0,0 +1,113 @@
+"""Signed interface-probe traces for the nested/twin shear presets (mechanism diagnostic).
+
+No pass gates. 24 fixed physical points on two interface-normal transects, away
+from corners, crossing L0->L1 and (nested only) L1->L2:
+  V: y = 64.1, across the vertical faces x = 16 (L0/L1) and x = 22 (L1/L2) of box 1;
+  H: x = 90.1, across the horizontal faces y = 24 (L0/L1) and y = 30 (L1/L2) of box 2.
+Each point is sampled at the ACTIVE leaf cell containing it, so cases with
+different grids (U128, U512, nested, twin) give different cell centers: the
+actual center and level are stored and analytic fields are evaluated there.
+Samples are POST-collision G at root synchronization only (root-step units);
+no substep content is observed, and grid-scale alternation can alias.
+"""
+
+import math
+
+import numpy as np
+
+import shear_wave as sw
+from oracle import CX, CY, U
+
+ROOT_N = 128
+OFFSETS = (-3.9, -1.9, -0.4, 0.1, 0.6, 2.1, 4.1, 5.6, 6.1, 6.6, 8.1, 10.1)  # from outer face
+TRANSECTS = {"V": dict(axis=0, face=16.0, other=64.1),   # x varies, y fixed
+             "H": dict(axis=1, face=24.0, other=90.1)}   # y varies, x fixed
+STEPS = tuple(range(0, 129)) + tuple(range(1920, 2049))  # root steps, both inclusive
+TRACE_MAX_BYTES = 600_000
+
+
+def points():
+    out = []
+    for name, t in TRANSECTS.items():
+        for offset in OFFSETS:
+            along = t["face"] + offset
+            xy = (along, t["other"]) if t["axis"] == 0 else (t["other"], along)
+            out.append((name, offset, xy))
+    return out
+
+
+def case_steps(n):
+    """Case lattice steps at the fixed root-physical sample times (U512: 4 per root step)."""
+    return [s * (n // ROOT_N) for s in STEPS]
+
+
+def select(masks, n=ROOT_N):
+    """Leaf cell containing each point: finest level whose ACTIVE mask covers it.
+
+    Points and stored centers are in root-physical units (root cell = 1); a case
+    with root N = n has cells of size ROOT_N / (n * 2**level) there.
+    """
+    if n % ROOT_N:
+        raise ValueError("case root N must be a multiple of 128")
+    probes = []
+    for name, offset, (x, y) in points():
+        for level in reversed(range(len(masks))):
+            scale = (n // ROOT_N) * 2 ** level
+            i, j = int(math.floor(x * scale)), int(math.floor(y * scale))
+            if masks[level][j, i]:
+                h = 1.0 / scale
+                probes.append(dict(transect=name, offset=offset, point=[x, y], level=level,
+                                   i=i, j=j, center=[h * (i + 0.5), h * (j + 0.5)]))
+                break
+        else:
+            raise ValueError(f"point {(x, y)} is not covered by an ACTIVE leaf")
+    return probes
+
+
+def sample(arrays, probes):
+    """Flat [rho, ux, uy] per probe from FP64 physical populations."""
+    row = []
+    for p in probes:
+        f = np.asarray(arrays[p["level"]][:, p["j"], p["i"]], dtype=float)
+        rho = float(f.sum())
+        row += [rho, float(CX @ f) / rho, float(CY @ f) / rho]
+    return row
+
+
+def new_trace(probes, case_viscosity, case):
+    """Times/coordinates are root-physical (N=128); nu_root = U*128/Re for all cases."""
+    return dict(schema=1, kind="mechanism diagnostic, no pass gates", case=case,
+                time_unit="root step (N=128)", sample_phase="postcollision G at root synchronization",
+                root_n=ROOT_N, viscosity_root=U * ROOT_N / sw.RE, case_viscosity=case_viscosity,
+                mode=list(sw.MODE), phase=sw.PHI,
+                amplitude=U, fields=["rho", "ux", "uy"], probes=probes, steps=[], values=[])
+
+
+def residuals(trace):
+    """Signed (drho, dux, duy) per step and probe, analytic field at each actual center."""
+    out = []
+    for step, row in zip(trace["steps"], trace["values"]):
+        res = []
+        for n, p in enumerate(trace["probes"]):
+            x, y = p["center"]
+            ux, uy, _ = sw.exact(x, y, step, trace["root_n"], trace["viscosity_root"])
+            rho, vx, vy = row[3 * n: 3 * n + 3]
+            res.append((rho - 1.0, vx - float(ux), vy - float(uy)))
+        out.append(res)
+    return out
+
+
+def validate(trace):
+    """Exact schedule, finite values, fixed probe count and byte cap."""
+    import json
+    if trace.get("steps") != list(STEPS) or len(trace["values"]) != len(STEPS):
+        raise ValueError("trace steps differ from the fixed schedule")
+    if len(trace["probes"]) != len(OFFSETS) * len(TRANSECTS):
+        raise ValueError("unexpected probe count")
+    width = 3 * len(trace["probes"])
+    if any(len(r) != width or not all(math.isfinite(v) for v in r) for r in trace["values"]):
+        raise ValueError("malformed or nonfinite trace row")
+    size = len(json.dumps(trace, allow_nan=False).encode())
+    if size > TRACE_MAX_BYTES:
+        raise ValueError(f"trace exceeds {TRACE_MAX_BYTES} bytes ({size})")
+    return size

tests/validation/taylor_green/test_shear_wave.py

@@ -930,3 +930,189 @@ class TwinPresetTests(unittest.TestCase):
             self.assertTrue(all(len(r["spurious_pressure_scaled_by_level"]) == 2 for r in rec["rows"]))
             self.assertEqual(selectors, [("projected_linear", [list(self.twin.REFINEMENTS[0])])])
             self.assertEqual(sorted(p.name for p in Path(tmp).iterdir()), ["projected_twin.json"])
+
+
+class ProbeTraceTests(unittest.TestCase):
+    """Signed interface-probe mechanism diagnostic (--probe-trace); no pass gates."""
+
+    def setUp(self):
+        import shear_nested
+        import shear_probes
+        import shear_twin
+        self.nested, self.probes, self.twin = shear_nested, shear_probes, shear_twin
+
+    def exact_arrays(self, masks, n, t):
+        """Equilibrium of the exact wave at root-physical time t on a case grid of root N=n."""
+        out = []
+        for level, m in enumerate(masks):
+            x, y = centers(m.shape, level)
+            x, y = x * 128 / n, y * 128 / n
+            ux, uy, _ = sw.exact(x, y, t, 128, U * 128 / sw.RE)
+            out.append(sw.equilibrium(np.ones_like(ux), ux, uy))
+        return out
+
+    def test_frozen_history_identities(self):
+        for name, digest in (("shear_wave.py", "81d5c8ee60864372887cddf9239580c82fb6750d956038f1e3ee9de61aa4e877"),
+                             ("shear_nested.py", "fa6cfe87fb4df8c6d732c637c885276833c81973b950c7724d734f8fc350bcb0"),
+                             ("shear_twin.py", "1cbf37f18e89aea577581631c4565f51e68d22a3739b698816ef64acbf2b2bb4"),
+                             ("oracle.py", ORACLE_SHA256_EEC4677)):
+            self.assertEqual(hashlib.sha256((HERE / name).read_bytes()).hexdigest(), digest, name)
+        self.assertEqual(run_case.SHEAR_PROBES_SHA256,
+                         hashlib.sha256((HERE / "shear_probes.py").read_bytes()).hexdigest())
+
+    def test_probe_selection_per_case(self):
+        sp = self.probes
+        self.assertEqual(len(sp.points()), 24)
+        corners = [(x, y) for level in (self.nested.REFINEMENTS[0],
+                                        [tuple(v / 2 for v in r) for r in self.nested.REFINEMENTS[1]])
+                   for x0, y0, x1, y1 in level for x in (x0, x1) for y in (y0, y1)]
+        for _, _, (x, y) in sp.points():
+            self.assertGreaterEqual(min(math.hypot(x - cx, y - cy) for cx, cy in corners), 8.0)
+        pattern = [0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2]
+        nested = sp.select(self.nested.leaf_masks())
+        self.assertEqual([p["level"] for p in nested], pattern * 2)
+        twin = sp.select(self.twin.leaf_masks())
+        self.assertEqual([p["level"] for p in twin], [min(v, 1) for v in pattern] * 2)
+        u128 = sp.select([np.ones((128, 128), bool)], 128)
+        u512 = sp.select([np.ones((512, 512), bool)], 512)
+        self.assertTrue(all(p["level"] == 0 for p in u128 + u512))
+        for case in (nested, twin, u128, u512):
+            for p in case:
+                # The stored center is the actual cell center containing the point (root units).
+                h = 1.0 / (2 ** p["level"]) if case is not u512 else 0.25
+                self.assertLessEqual(abs(p["center"][0] - p["point"][0]), h / 2)
+                self.assertLessEqual(abs(p["center"][1] - p["point"][1]), h / 2)
+        # Different grids give different centers: no pretended collocation.
+        self.assertNotEqual(u128[3]["center"], u512[3]["center"])
+        self.assertEqual(sp.case_steps(512)[:3], [0, 4, 8])
+        self.assertEqual(sp.case_steps(128), list(sp.STEPS))
+
+    def test_residuals_vanish_for_the_exact_field(self):
+        sp = self.probes
+        for masks, n in ((self.nested.leaf_masks(), 128), ([np.ones((512, 512), bool)], 512)):
+            probes = sp.select(masks, n)
+            trace = sp.new_trace(probes, U * n / sw.RE, ["x", n, len(masks), "c"])
+            for t in (0, 100, 2048):
+                trace["steps"].append(t)
+                trace["values"].append(sp.sample(self.exact_arrays(masks, n, t), probes))
+            for res in sp.residuals(trace):
+                self.assertLess(max(abs(v) for r in res for v in r), 1e-14)
+
+    def test_validate_schedule_values_and_cap(self):
+        sp = self.probes
+        probes = sp.select(self.nested.leaf_masks())
+        trace = sp.new_trace(probes, 0.0384, ["n", 128, 3, "projected_linear"])
+        trace["steps"] = list(sp.STEPS)
+        trace["values"] = [[1.0, 0.01, -0.02] * 24 for _ in sp.STEPS]
+        size = sp.validate(trace)
+        self.assertLess(size, sp.TRACE_MAX_BYTES)
+        realistic = dict(trace, values=[[1.0000012345678901, 0.012345678901234567,
+                                         -0.023456789012345678] * 24 for _ in sp.STEPS])
+        self.assertLess(sp.validate(realistic), sp.TRACE_MAX_BYTES)
+        for bad in (dict(trace, steps=trace["steps"][:-1]),
+                    dict(trace, values=trace["values"][:-1] + [[float("nan")] * 72]),
+                    dict(trace, probes=trace["probes"][:-1])):
+            with self.assertRaises(ValueError):
+                sp.validate(bad)
+        with patch.object(sp, "TRACE_MAX_BYTES", 1000):
+            with self.assertRaises(ValueError):
+                sp.validate(trace)
+
+    def test_probe_runner_options(self):
+        repo = Path(__file__).resolve().parents[3]
+        with tempfile.TemporaryDirectory() as tmp:
+            base = ["run_case.py", "--repo", str(repo), "--expected-revision", "f" * 40,
+                    "--output", str(Path(tmp) / "case.json"), "--dry-run", "--re", "100",
+                    "--flow", "shear_wave", "--probe-trace"]
+            nest = ["--diagnostic", run_case.NESTED_SHEAR]
+            tw = ["--diagnostic", run_case.TWIN_SHEAR]
+            proj = ["--explosion", "projected_linear"]
+            for extra, valid in (
+                    (["--n", "128", "--depth", "1"] + nest, True),
+                    (["--n", "512", "--depth", "1"] + nest, True),
+                    (["--n", "128", "--depth", "3"] + proj + nest, True),
+                    (["--n", "128", "--depth", "2"] + proj + tw, True),
+                    (["--n", "128", "--depth", "3"] + nest, False),          # constant d3
+                    (["--n", "128", "--depth", "1"] + tw, False),
+                    (["--n", "128", "--depth", "3"] + proj + nest + ["--final-snapshot"], False),
+                    (["--n", "128", "--depth", "2"], False)):
+                argv = base + extra
+                with self.subTest(argv=argv), patch.object(sys, "argv", argv), \
+                     patch.object(run_case, "verify_source", return_value={}) as verify, \
+                     patch.object(run_case, "runtime_source", return_value={}), \
+                     redirect_stdout(io.StringIO()), patch("sys.stderr", io.StringIO()):
+                    if valid:
+                        self.assertEqual(run_case.main(), 0)
+                        self.assertIs(verify.call_args.kwargs["probes"], True)
+                    else:
+                        with self.assertRaises(SystemExit):
+                            run_case.main()
+
+    def run_mocked(self, n, depth, diagnostic, scheme, masks):
+        from types import SimpleNamespace
+        from contextlib import nullcontext
+        arrays = self.exact_arrays(masks, n, 0.0)
+        counts = []
+
+        class FakeRun:
+            def __init__(self, *args, **kwargs):
+                self.step_count = 0
+
+            def __enter__(self):
+                return self
+
+            def __exit__(self, *args):
+                pass
+
+            def initialize(self, *given):
+                pass
+
+            def step(self, count):
+                self.step_count += count
+                counts.append(self.step_count)
+
+            def download(self):
+                return arrays[0] if depth == 1 else arrays
+
+        modules = {"CelerisLab.cuda.context": SimpleNamespace(CudaContext=FakeRun),
+                   "CelerisLab.lbm.block_stepper": SimpleNamespace(PeriodicBlockStepper=FakeRun),
+                   "CelerisLab.lbm.refinement_stepper": SimpleNamespace(FixedPatchStepper=FakeRun)}
+        repo = Path(__file__).resolve().parents[3]
+        with tempfile.TemporaryDirectory() as tmp:
+            output = Path(tmp) / "probe.json"
+            argv = ["run_case.py", "--repo", str(repo), "--re", "100", "--n", str(n),
+                    "--depth", str(depth), "--explosion", scheme, "--flow", "shear_wave",
+                    "--diagnostic", diagnostic, "--probe-trace", "--output", str(output),
+                    "--expected-revision", "f" * 40, "--run-gpu"]
+            with patch.dict(sys.modules, modules), patch.object(sys, "argv", argv), \
+                 patch.dict("os.environ", CUDA_VISIBLE_DEVICES="mock"), \
+                 patch.object(run_case, "verify_source", return_value=dict(revision="f" * 40)), \
+                 patch.object(run_case, "runtime_source", return_value={}), \
+                 patch.object(run_case, "runtime_identity", return_value={}), \
+                 patch.object(run_case, "build_capture", return_value=nullcontext()), \
+                 redirect_stdout(io.StringIO()):
+                run_case.main()
+            return json.loads(output.read_text()), counts, sorted(p.name for p in Path(tmp).iterdir())
+
+    def test_mocked_probe_traces(self):
+        sp = self.probes
+        rec, counts, files = self.run_mocked(128, 3, run_case.NESTED_SHEAR, "projected_linear",
+                                             self.nested.leaf_masks())
+        trace = rec["probe_trace"]
+        self.assertEqual(trace["steps"], list(sp.STEPS))
+        self.assertEqual(len(trace["values"]), len(sp.STEPS))
+        self.assertLessEqual(trace["bytes"], sp.TRACE_MAX_BYTES)
+        self.assertEqual(trace["case"], [run_case.NESTED_SHEAR, 128, 3, "projected_linear"])
+        self.assertEqual(len(rec["rows"]), 5)  # qualification samples preserved
+        self.assertEqual(files, ["probe.json"])  # no field dumps
+        self.assertEqual(sp.validate(trace) <= sp.TRACE_MAX_BYTES, True)
+        rec512, counts512, _ = self.run_mocked(512, 1, run_case.NESTED_SHEAR, "constant",
+                                               [np.ones((512, 512), bool)])
+        t512 = rec512["probe_trace"]
+        self.assertEqual(t512["steps"], list(sp.STEPS))  # root-physical times
+        self.assertEqual(t512["viscosity_root"], U * 128 / sw.RE)
+        self.assertEqual(t512["case_viscosity"], U * 512 / sw.RE)
+        self.assertIn(4 * 1920, counts512)  # U512 advanced to 4x lattice steps
+        self.assertTrue(all(p["level"] == 0 for p in t512["probes"]))
+        res = sp.residuals(t512)
+        self.assertLess(abs(res[0][0][1]), 1e-14)  # exact field at t=0 at actual centers