test_block_step_sync: locate block_stepper without importing it

The module-level import of CelerisLab.lbm.block_stepper loaded
pycuda.driver into the shared test process, unlike the other CPU host
tests (which load mocked copies via AST). find_spec gives the source path
without the import. Test content and assertions unchanged.

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

docs/compatibility.md

@@ -0,0 +1,35 @@
+# Compatibility and API changes
+
+Changes that can make a previously accepted setup or call behave differently.
+Each entry names the commit, the affected entry points and the observable result.
+
+## 87c3a31: cut-link ownership is validated for every body
+
+`validate_cut_link_ownership` (`body/preprocess/flag_overlay.py`) now runs for
+passive bodies as well as Jet/aperture bodies. Every cut link must have an
+active fluid receiver (FLUID set, SOLID clear), and each `(receiver, direction)`
+population may have only one writer.
+
+Rejected with `ValueError`, before flags, populations or device buffers change:
+
+1. Overlapping or touching passive solids whose cut links write the same
+   population ("Boundary requires unique cut-link writers").
+2. A body whose cut-link receiver is a non-fluid cell, e.g. a channel wall row
+   or the inlet/outlet column ("Cut links must have active fluid receivers").
+
+Both were previously accepted silently for passive bodies; the interpolated
+boundary then had undefined results. The rules apply to the cut links of all
+bodies together, so a passive body overlapping a Jet body is rejected as well.
+The check runs in `Simulation.initialize`, runtime body sync (`sync_bodies`),
+checkpoint save and checkpoint load (`common/checkpoint.py`). Load validates
+before restoring anything, so a checkpoint whose saved geometry violates either
+rule no longer restores; rebuild the geometry with separated bodies.
+
+## 87c3a31: `Simulation.read_surface_exchange` accepts Jet-driven bodies
+
+`read_surface_exchange(id)` previously required `SurfaceActuation` on the body.
+It now also accepts bodies that own Jet aperture links, including Jets with zero
+or negative amplitude (the aperture test is geometric). Bodies with neither
+`SurfaceActuation` nor aperture links still raise
+`ValueError("Body has no surface actuation")`. The value is unchanged: the latest
+step's wall/aperture exchange, not total nozzle force.

docs/validation/dfg_re20.md

@@ -0,0 +1,49 @@
+# DFG 2D-1 (Re 20) against the official reference intervals
+
+`tests/validation/run_dfg_re20.py` runs the steady DFG 2D-1 case on the refined
+channel root (D40 root, D80 wide fine patch) and on the uniform D80 grid. Geometry,
+inflow, initial state, layout and observation are those of the pinned Re100 runner
+(`run_refined_dfg.py`, `re100_helpers.py`, both unchanged); only the lattice
+viscosity is set for Re 20 (root tau 0.56, fine 0.62, uniform 0.62, lattice speed 0.01).
+
+Label: "compared against the official reference bands", not a full benchmark
+certification: the outlet is pressure-NEQ (`neq_extrap`), not do-nothing, and the
+inlet is local Zou-He.
+
+## Quantities
+
+- Cd = 2 Fx / (rho U^2 D), Cl = 2 Fy / (rho U^2 D), interval means over each sample.
+- p_diff = p(0.15, 0.2) - p(0.25, 0.2), linear extrapolation to the exact surface
+  points (`re100_helpers.pressure_probes`). The helper returns dp / (rho U^2) in
+  scaled units; p_diff = dp * 0.04 for U = 0.2 m/s, rho = 1. Both are reported.
+- Official intervals (Schaefer & Turek 1996): Cd [5.57, 5.59], Cl [0.0104, 0.0110],
+  p_diff [0.1172, 0.1176].
+
+## Bounds and steady rule (frozen before any run)
+
+- End time t* = 60 (30 s physical): 240,000 refined root steps, 480,000 uniform steps.
+- One sample every t* = 0.5 (0.25 s) at the end of a synchronized root step; Cd and
+  Cl are interval means. Trailing window t* 50..60 inclusive (21 samples).
+- Steady rule (our engineering choice, not official): in the window, for each of Cd,
+  Cl, p_diff, max - min <= 10 % and |least-squares slope x window| <= 2 % of its
+  official interval width (range 2e-3 / 6e-5 / 4e-5; slope 4e-4 / 1.2e-5 / 8e-6).
+  Absolute limits, so the small Cl is never a divisor. Never relaxed after a run.
+- The verdict uses the window mean; the distance to the nearest interval edge is
+  reported for each quantity (negative outside).
+- Force sign: fluid on cylinder, x along the inflow, y up from the bottom wall; Cd > 0
+  and the official Cl > 0 (`cd_positive` is reported).
+- Not steady at t* = 60: status `NOT_STEADY`, no interval verdict. Health limits
+  (Mach <= 0.2, rho in [0.5, 1.5]) stop the run with `UNHEALTHY`.
+- Steady: `WITHIN_INTERVALS` or `OUTSIDE_INTERVALS`, with the per-quantity verdict.
+
+## Run
+
+```
+python tests/validation/run_dfg_re20.py --output <new-dir> --plan
+CUDA_VISIBLE_DEVICES=<n> python tests/validation/run_dfg_re20.py --case both \
+    --collision SRT --output <new-dir> --max-wall-seconds <cap>
+```
+
+The uniform D80 control separates method error from refinement error. `--plan` writes `plan.json` only. A GPU run writes `plan.json` and `result.json`
+(per case: lattice parameters, every sample, window ranges, status, verdict).
+Host contract tests: `tests/unit/test_dfg_re20.py`.

src/CelerisLab/lbm/block_stepper.py

@@ -204,11 +204,21 @@ class PeriodicBlockStepper:
                 raise
 
     def step(self, n=1):
-        """Advance n steps: halo refresh, then core pull/collide, per step."""
+        """Advance n steps: halo refresh, then core pull/collide, per step.
+
+        All n steps are enqueued in order on the one stream and the host waits
+        once before returning (same launches and bits as n step(1) calls).
+        Asynchronous device errors surface at that single wait, so they are not
+        attributed to the step that caused them; call step(1) repeatedly to
+        localize one. On any failure the instance is failed and step_count keeps
+        its entry value.
+        """
         if isinstance(n, bool) or not isinstance(n, Integral) or n < 0:
             raise ValueError("n must be a nonnegative integer")
         with self._lock, self.ctx.activate():
             self._require("ready")
+            if n == 0:
+                return
             kw = dict(block=self._block, stream=self._stream)
             try:
                 self._status = "advancing"
@@ -224,9 +234,9 @@ class PeriodicBlockStepper:
                     self._step(self._current, self._spare, np.uint32(self.layout.halo),
                                self._tables.get("flags", np.uintp(0)), np.uintp(0),  # no kind
                                grid=self._grid, **kw)
-                    self._stream.synchronize()
-                    self._current, self._spare = self._spare, self._current
-                    self._step_count += 1
+                    self._current, self._spare = self._spare, self._current  # host handles only
+                self._stream.synchronize()
+                self._step_count += n
                 self._status = "ready"
             except BaseException:
                 self._status = "failed"

src/CelerisLab/lbm/refinement_stepper.py

@@ -822,27 +822,41 @@ class FixedPatchStepper:
                             grid=((pairs + self._block[0] - 1) // self._block[0], 1), **kw)
 
     def step(self, n=1):
-        """Advance G→G by n base intervals: level l takes 2**l hops per interval."""
+        """Advance G→G by n base intervals: level l takes 2**l hops per interval.
+
+        All n intervals are enqueued in order on the one stream, and the host waits
+        once, before returning; step(0) launches and waits for nothing. On return
+        the device work is complete, step_count and an open window cover all n
+        intervals, and the raw-exchange ledger holds the final interval only.
+
+        Asynchronous device errors surface at that single wait, so they are not
+        attributed to the interval that caused them; call step(1) repeatedly to
+        localize one (same launches, same bits). On any failure the instance is
+        failed, the exchange ledger and window are invalid, and step_count keeps
+        its value from entry to the failing call.
+        """
         if isinstance(n, bool) or not isinstance(n, Integral) or n < 0:
             raise ValueError("n must be a nonnegative integer")
         with self._lock, self.ctx.activate():
             self._require("ready")
+            if n == 0:
+                return
             root = self._levels[0]
             kw = dict(block=self._block, stream=self._stream)
             try:
                 self._status = "advancing"
+                self._exchange_valid = False
                 for _ in range(n):
-                    if self._observe_raw_exchange:
-                        self._exchange_valid = False
+                    if self._observe_raw_exchange:  # stream-ordered: holds the last interval only
                         for b in self._level_boundaries:
                             if b is not None:
                                 cuda.memset_d8_async(b.interval_integral, 0,
                                                      len(b.tables["cell"]) * 32, self._stream)
                     self._hop(0, root.current, root.spare, **kw)
-                    self._stream.synchronize()
-                    root.current, root.spare = root.spare, root.current
-                    self._step_count += 1
-                    self._exchange_valid = self._observe_raw_exchange
+                    root.current, root.spare = root.spare, root.current  # host handles only
+                self._stream.synchronize()
+                self._step_count += n
+                self._exchange_valid = self._observe_raw_exchange
                 self._status = "ready"
             except BaseException:
                 self._window_valid = False

tests/unit/test_block_step_sync.py

@@ -0,0 +1,104 @@
+"""Host contract of PeriodicBlockStepper.step(n): one stream wait per call (mocked driver)."""
+import ast
+from contextlib import nullcontext
+import importlib.util
+from pathlib import Path
+import types
+from unittest.mock import Mock
+import unittest
+
+import numpy as np
+
+from CelerisLab.config import LBMConfig
+from CelerisLab.lbm.blocks import BlockLayout2D
+
+
+def host_module():
+    """block_stepper with pycuda/compiler replaced by mocks; no GPU context."""
+    # Locate the source without importing it: importing block_stepper loads pycuda.driver,
+    # which CPU-only import-boundary tests (test_build_contracts) must never see.
+    path = Path(importlib.util.find_spec("CelerisLab.lbm.block_stepper").origin)
+    tree = ast.parse(path.read_text())
+    tree.body = [node for node in tree.body if not (
+        isinstance(node, ast.Import) and any(a.name == "pycuda.driver" for a in node.names)
+        or isinstance(node, ast.ImportFrom) and node.module == "cuda")]
+    module = types.ModuleType("CelerisLab.lbm._test_host_block_stepper")
+    module.__package__ = "CelerisLab.lbm"
+    module.__dict__.update(cuda=Mock(), compiler=Mock())
+    exec(compile(tree, str(path), "exec"), module.__dict__)
+    return module
+
+
+def backend():
+    host, log = host_module(), []
+    host.cuda.mem_get_info.return_value = (1 << 40, 1 << 40)
+    host.cuda.mem_alloc.side_effect = lambda size: Mock(size=size)  # distinct buffers
+    module = Mock()
+    module.get_global.return_value = (Mock(), 20)
+    module.get_function.side_effect = lambda name: Mock(
+        side_effect=lambda *a, **k: log.append((name, a[:2])))
+    host.compiler.load_module.return_value = module
+    cfg = LBMConfig(nx=16, ny=16, viscosity=.1)
+    run = host.PeriodicBlockStepper(cfg, BlockLayout2D.build(16, 16, 8, 8),
+                                    context=Mock(activate=nullcontext, sm_arch="sm_70"))
+    run._stream.synchronize.side_effect = lambda: log.append(("sync", ()))
+    run.initialize(np.full((9, 16, 16), 1 / 9))
+    log.clear()
+    return run, log
+
+
+class BlockStepSync(unittest.TestCase):
+    def test_step_n_waits_once_and_matches_n_single_steps(self):
+        many, log_many = backend()
+        many.step(3)
+        single, log_single = backend()
+        for _ in range(3):
+            single.step(1)
+        names = [name for name, _ in log_many]
+        self.assertEqual(names.count("sync"), 1)
+        self.assertEqual(names[-1], "sync")
+        self.assertEqual(names.count("BlockStep"), 3)
+        # Same launch order and buffer roles as three step(1) calls.
+        strip = lambda log: [(n, tuple(id(x) for x in a)) for n, a in log if n != "sync"]
+        roles = lambda run, log: [(n, tuple({id(run._current): "c", id(run._spare): "s"}.get(x, "-")
+                                            for x in ids)) for n, ids in strip(log)]
+        self.assertEqual([n for n, _ in strip(log_many)], [n for n, _ in strip(log_single)])
+        self.assertEqual(roles(many, log_many), roles(single, log_single))
+        self.assertIsNot(many._current, many._spare)
+        self.assertEqual((many.step_count, single.step_count), (3, 3))
+
+    def test_step_zero_launches_and_waits_for_nothing(self):
+        run, log = backend()
+        run.step(0)
+        self.assertEqual((log, run.step_count), ([], 0))
+
+    def test_failure_keeps_entry_step_count(self):
+        run, log = backend()
+        run.step(1)
+        run._stream.synchronize.side_effect = RuntimeError("device fault")
+        with self.assertRaises(RuntimeError):
+            run.step(2)
+        self.assertEqual((run.step_count, run._status), (1, "failed"))
+
+    def test_launch_failure_mid_batch_fails_closed_without_waiting(self):
+        run, log = backend()
+        calls = []
+        step = run._step
+
+        def failing(*args, **kwargs):
+            calls.append(1)
+            if len(calls) == 2:
+                raise RuntimeError("launch failure")
+            return step(*args, **kwargs)
+
+        run._step = failing
+        with self.assertRaises(RuntimeError):
+            run.step(3)
+        self.assertNotIn("sync", [name for name, _ in log])
+        self.assertEqual((run.step_count, run._status), (0, "failed"))
+        with self.assertRaises(RuntimeError):
+            run.download()
+
+
+if __name__ == "__main__":
+    unittest.main()

tests/unit/test_dfg_re20.py

@@ -0,0 +1,79 @@
+"""CPU-only contract tests for run_dfg_re20. Not CUDA or DFG acceptance evidence."""
+
+import sys
+from pathlib import Path
+import unittest
+
+sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "validation"))
+import run_dfg_re20 as v  # noqa: E402
+import run_refined_dfg as rr  # noqa: E402
+
+h, _ = rr.helpers(rr.DEFAULT_SOURCE)
+
+
+def rows(values, tstar_end=v.TSTAR_MAX):
+    n = int(round(tstar_end / v.SAMPLE_TSTAR))
+    return [dict(tstar=(i + 1) * v.SAMPLE_TSTAR, **values(i)) for i in range(n)]
+
+
+class Re20Contract(unittest.TestCase):
+    def test_lattice_parameters_and_bounds(self):
+        root, uniform = v.lattice_parameters(40), v.lattice_parameters(80)
+        self.assertAlmostEqual(root["tau"], 0.56, places=12)
+        self.assertAlmostEqual(uniform["tau"], 0.62, places=12)
+        self.assertEqual((root["steps"], root["interval"]), (240000, 2000))
+        self.assertEqual((uniform["steps"], uniform["interval"]), (480000, 4000))
+        self.assertAlmostEqual(v.plan()["refined_fine_tau"], 0.62, places=12)
+
+    def test_config_only_changes_viscosity(self):
+        base = v._re100_config(h, 40, "SRT", "none")
+        cfg, body, p = v.config(h, 40, "SRT")
+        self.assertAlmostEqual(cfg["physics"]["viscosity"], 0.01 * 40 / 20, places=15)
+        self.assertEqual((p["re"], p["viscosity"]), (20.0, cfg["physics"]["viscosity"]))
+        base[0]["physics"]["viscosity"] = cfg["physics"]["viscosity"]
+        self.assertEqual(base[0], cfg)
+        self.assertEqual(base[1], body)
+        with self.assertRaisesRegex(ValueError, "no Jets"):
+            v.config(h, 40, "MRT", "zero")
+
+    def test_units(self):
+        s = v.sample([2.0, -0.5], 4, 40, 3.0)
+        den = rr.SPEED**2 * 40 / 2
+        self.assertAlmostEqual(s["cd"], 2.0 / 4 / den)
+        self.assertAlmostEqual(s["cl"], -0.5 / 4 / den)
+        self.assertAlmostEqual(s["p_diff"], 3.0 * 0.04)
+        self.assertEqual(s["dp_scaled"], 3.0)
+
+    def test_steady_and_verdict(self):
+        inside = dict(cd=5.58, cl=0.0106, p_diff=0.1174)
+        out = v.summarize(rows(lambda i: inside))
+        self.assertEqual((out["status"], out["steady"], out["window_samples"]),
+                         ("WITHIN_INTERVALS", True, 21))
+        self.assertTrue(out["cd_positive"])
+        self.assertAlmostEqual(out["edge_distance"]["cl"], 0.0002)
+        self.assertTrue(all(out["verdict"].values()))
+        low = v.summarize(rows(lambda i: dict(inside, cl=0.0100)))
+        self.assertEqual(low["status"], "OUTSIDE_INTERVALS")
+        self.assertFalse(low["verdict"]["cl"])
+
+    def test_slow_drift_gives_no_verdict(self):
+        # Range 2e-5 passes (limit 6e-5); slope * window 2e-5 exceeds 1.2e-5.
+        out = v.summarize(rows(lambda i: dict(cd=5.58, cl=0.0106 + 1e-6 * i, p_diff=0.1174)))
+        self.assertEqual((out["status"], out["verdict"]), ("NOT_STEADY", None))
+        self.assertLessEqual(out["ranges"]["cl"], out["range_limits"]["cl"])
+        self.assertGreater(out["drifts"]["cl"], out["slope_limits"]["cl"])
+
+    def test_noise_range(self):
+        quiet = v.summarize(rows(lambda i: dict(cd=5.58, cl=0.0106 + 2e-5 * (i % 2), p_diff=0.1174)))
+        self.assertTrue(quiet["steady"])  # within the range limit, no trend
+        noisy = v.summarize(rows(lambda i: dict(cd=5.58, cl=0.0106 + 1e-4 * (i % 2), p_diff=0.1174)))
+        self.assertEqual(noisy["status"], "NOT_STEADY")
+        self.assertGreater(noisy["ranges"]["cl"], noisy["range_limits"]["cl"])
+
+    def test_incomplete_run_gives_no_verdict(self):
+        out = v.summarize(rows(lambda i: dict(cd=5.58, cl=0.0106, p_diff=0.1174), tstar_end=30))
+        self.assertEqual((out["status"], out["verdict"]), ("INCOMPLETE", None))
+
+
+if __name__ == "__main__":
+    unittest.main()

tests/unit/test_step_sync.py

@@ -0,0 +1,170 @@
+"""Host contract of FixedPatchStepper.step(n): one stream wait per call (mocked driver).
+
+All n base intervals are enqueued in order on the backend stream and the host
+waits once before step returns. Launch order and buffer arguments must equal n
+step(1) calls; completion metadata, raw-exchange and window ledgers, and failure
+invalidation are published only after that wait.
+"""
+from contextlib import nullcontext
+from unittest.mock import Mock
+import unittest
+
+import numpy as np
+
+from CelerisLab.config import LBMConfig
+from tests.unit.test_raw_exchange import body
+from tests.unit.test_refinement_three_level_multipatch import host_module
+
+BODY_LEVEL, HOPS = 3, 8  # body on L3: 2**3 hops per base interval
+
+
+def backend(window=False):
+    """Four-level host backend; every launch, memset and stream wait goes to one log."""
+    host = host_module()
+    log, allocations = [], []
+
+    def allocate(size):
+        allocations.append(Mock(size=size))
+        return allocations[-1]
+
+    module = Mock()
+    module.get_global.return_value = (Mock(), 20)
+    module.get_function.side_effect = lambda name: Mock(
+        side_effect=lambda *a, **k: log.append((name, a, k)))
+    host.compiler.load_module.return_value = module
+    host.cuda.mem_alloc.side_effect = allocate
+    host.cuda.memcpy_dtoh.side_effect = lambda array, pointer: array.fill(0)
+    host.cuda.memset_d8_async.side_effect = lambda *a: log.append(("memset", a, {}))
+    run = host.FixedPatchStepper(
+        LBMConfig(nx=48, ny=48, viscosity=.02), context=Mock(activate=nullcontext, sm_arch="sm_70"),
+        refinements=[[(6, 6, 26, 26)], [(20, 20, 44, 44)], [(56, 56, 72, 72)]],
+        level_blocks=[(4, 4)] * 4, level_objects=[(), (), (), (body(),)],
+        observe_raw_exchange=True, observe_window_exchange=window)
+    run._stream.synchronize.side_effect = lambda: log.append(("sync", (), {}))
+    run.initialize(*[np.full(level.shape, 1 / 9, np.float32) for level in run._levels])
+    log.clear()
+    return run, log, allocations
+
+
+def signature(log, allocations):
+    """Event names with device-pointer arguments replaced by allocation index."""
+    index = {id(a): i for i, a in enumerate(allocations)}
+    return [(name, tuple(index.get(id(a), "-") for a in args)) for name, args, _ in log]
+
+
+def names(log):
+    return [name for name, _, _ in log]
+
+
+class StepSync(unittest.TestCase):
+    def test_step1_waits_once_after_all_launches(self):
+        run, log, _ = backend()
+        run.step(1)
+        self.assertEqual(names(log).count("sync"), 1)
+        self.assertEqual(names(log)[-1], "sync")
+        self.assertEqual(names(log)[0], "memset")
+        self.assertEqual(names(log).count("BlockStep"), 1 + 2 + 4 + 8)
+        self.assertEqual(run.step_count, 1)
+        run.close()
+
+    def test_step_n_equals_n_step1_launches_with_one_wait(self):
+        single, single_log, single_alloc = backend()
+        for _ in range(3):
+            single.step(1)
+        batched, batched_log, batched_alloc = backend()
+        batched.step(3)
+        expected = [e for e in signature(single_log, single_alloc) if e[0] != "sync"]
+        got = signature(batched_log, batched_alloc)
+        self.assertEqual(got[-1][0], "sync")
+        self.assertEqual(names(batched_log).count("sync"), 1)
+        self.assertEqual(got[:-1], expected)  # same order, same ping-pong buffers
+        self.assertEqual(names(batched_log).count("memset"), 3)
+        self.assertEqual(batched.step_count, single.step_count)
+        index = [{id(a): i for i, a in enumerate(alloc)} for alloc in (single_alloc, batched_alloc)]
+        self.assertEqual(index[1][id(batched._levels[0].current)],
+                         index[0][id(single._levels[0].current)])
+        single.close()
+        batched.close()
+
+    def test_step0_launches_and_waits_for_nothing(self):
+        run, log, _ = backend()
+        run.step(2)
+        log.clear()
+        run.step(0)
+        self.assertEqual(log, [])
+        self.assertEqual(run.step_count, 2)
+        value = run.read_raw_exchange(0, level=BODY_LEVEL)  # step(0) keeps the ledger valid
+        self.assertEqual((value["time_start"], value["time_end"]), (1, 2))
+        run.close()
+
+    def test_completion_metadata_covers_all_intervals(self):
+        run, _, _ = backend()
+        with self.assertRaisesRegex(RuntimeError, "completed"):
+            run.read_raw_exchange(0, level=BODY_LEVEL)
+        run.step(3)
+        value = run.read_raw_exchange(0, level=BODY_LEVEL)
+        self.assertEqual((value["time_start"], value["time_end"], value["sample_count"]),
+                         (2, 3, HOPS))  # step(n) retains only the final base interval
+        run.step(2)
+        self.assertEqual(run.read_raw_exchange(0, level=BODY_LEVEL)["time_end"], 5)
+        run.close()
+
+    def test_window_accounting_across_batched_calls(self):
+        windows = []
+        for calls in ((1, 3), (4,), (1, 1, 1, 1)):
+            run, log, _ = backend(window=True)
+            run.step(1)
+            run.start_exchange_window()
+            log.clear()
+            for n in calls:
+                run.step(n)
+            self.assertEqual(names(log).count("BlockAccumulateExchange"), 2 * 4 * HOPS)
+            self.assertEqual(names(log).count("sync"), len(calls))
+            out = run.read_exchange_window(0, level=BODY_LEVEL)
+            windows.append((out["time_start"], out["time_end"], out["sample_count"]))
+            run.close()
+        self.assertEqual(windows, [(1., 5, 4 * HOPS)] * 3)
+
+    def test_failed_wait_publishes_nothing(self):
+        run, _, _ = backend(window=True)
+        run.step(2)
+        run.start_exchange_window()
+        run._stream.synchronize.side_effect = RuntimeError("device failure")
+        with self.assertRaisesRegex(RuntimeError, "device failure"):
+            run.step(3)
+        self.assertEqual(run.step_count, 2)  # value at entry to the failing call
+        self.assertFalse(run._exchange_valid)
+        self.assertFalse(run._window_valid)
+        for read in (lambda: run.read_raw_exchange(0, level=BODY_LEVEL),
+                     lambda: run.read_exchange_window(0, level=BODY_LEVEL), run.download,
+                     lambda: run.step(1)):
+            with self.assertRaisesRegex(RuntimeError, "failed"):
+                read()
+        run._stream.synchronize.side_effect = None
+        run.close()
+
+    def test_launch_error_mid_batch_skips_wait_and_fails(self):
+        run, log, _ = backend()
+        run.step(1)
+        log.clear()
+        hop, calls = run._hop, []
+
+        def failing(index, *args, **kw):
+            if index == 0:
+                calls.append(1)
+                if len(calls) == 2:
+                    raise RuntimeError("launch failure")
+            return hop(index, *args, **kw)
+        run._hop = failing
+        with self.assertRaisesRegex(RuntimeError, "launch failure"):
+            run.step(3)
+        self.assertNotIn("sync", names(log))
+        self.assertEqual(run.step_count, 1)
+        self.assertFalse(run._exchange_valid)
+        with self.assertRaisesRegex(RuntimeError, "failed"):
+            run.read_raw_exchange(0, level=BODY_LEVEL)
+        run.close()
+
+
+if __name__ == "__main__":
+    unittest.main()

tests/validation/multigrid/test_refinement_three_level.py

@@ -80,11 +80,12 @@ class ThreeLevelCPU(unittest.TestCase):
         f, g = run._levels[1:]
         b = g.boundary
         events, feeds = [], []
-        advance, gather = run._advance_child, run._gather_halo
-        def record_feed(buffer, **kw):
-            self.assertIs(buffer, (f.current, f.spare)[len(feeds)])
-            feeds.append(buffer)
-            advance(buffer, **kw)
+        advance, gather = run._advance, run._gather_halo
+        def record_feed(index, buffer, **kw):  # index 2: L1 hop feeds L2
+            if index == 2:
+                self.assertIs(buffer, (f.current, f.spare)[len(feeds)])
+                feeds.append(buffer)
+            advance(index, buffer, **kw)
         def record_halo(level, buffer, **kw):
             if level is g:
                 events.append(("halo", buffer))
@@ -95,7 +96,7 @@ class ThreeLevelCPU(unittest.TestCase):
             self.assertIs(kind, g.kind)
             events.append(("step", src))
         g.kernel.side_effect = record_step
-        with patch.object(run, "_advance_child", side_effect=record_feed), \
+        with patch.object(run, "_advance", side_effect=record_feed), \
                 patch.object(run, "_gather_halo", side_effect=record_halo):
             run.step()
         self.assertEqual(feeds, [f.current, f.spare])
@@ -105,8 +106,8 @@ class ThreeLevelCPU(unittest.TestCase):
         self.assertGreater(len(np.unique(b.tables["block"])), 1)
         # Negative control: reusing L1's old source at the midpoint must fail.
         feeds.clear()
-        with patch.object(run, "_advance_child", side_effect=lambda buffer, **kw:
-                          record_feed(f.current, **kw)):
+        with patch.object(run, "_advance", side_effect=lambda index, buffer, **kw:
+                          record_feed(index, f.current if index == 2 else buffer, **kw)):
             with self.assertRaises(AssertionError):
                 run.step()
         run.close()
@@ -172,11 +173,12 @@ class ThreeLevelGPU(unittest.TestCase):
             p1, p2 = run.patch, run.child_patch
             run.initialize(state(16, 16, 1), state(*p1.fine_shape, 2), state(*p2.fine_shape, 3))
             f = run._levels[1]
-            feeds, original = [], run._advance_child
-            def record(buffer, **kw):
-                feeds.append(int(buffer))
-                original(buffer, **kw)
-            with patch.object(run, "_advance_child", side_effect=record):
+            feeds, original = [], run._advance
+            def record(index, buffer, **kw):  # index 2: L1 hop feeds L2
+                if index == 2:
+                    feeds.append(int(buffer))
+                original(index, buffer, **kw)
+            with patch.object(run, "_advance", side_effect=record):
                 cur, spare = int(f.current), int(f.spare)
                 run.step(1)
             self.assertEqual(feeds, [cur, spare])  # current G, then the new G of hop 1

tests/validation/run_dfg_re20.py

@@ -0,0 +1,222 @@
+"""Steady DFG 2D-1 (Re 20): refined and uniform runs against the official reference intervals.
+
+Not a full benchmark certification. Declared deviations from the official setup:
+pressure-NEQ outlet (``neq_extrap``) instead of do-nothing, local Zou-He inlet.
+Geometry, inflow profile, initial state, refined layout and observation reuse the
+pinned Re100 runner (``run_refined_dfg``) and ``re100_helpers`` unchanged; only the
+lattice viscosity is set for Re 20 (nu = U_lattice * D / 20 at the root level).
+
+Units. The helpers use scaled DFG units with mean inflow speed 1, so
+Cd = 2 Fx / (rho U^2 D) and Cl = 2 Fy / (rho U^2 D) are scale-free, and
+``pressure_probes`` returns dp / (rho U^2). The official p_diff = p(0.15, 0.2) -
+p(0.25, 0.2) uses U = 0.2 m/s, rho = 1, so p_diff = dp * 0.04. Both are reported;
+only p_diff is compared with the official interval. Forces use the existing sign
+convention: F is the force of the fluid on the cylinder, i.e. the refined raw
+exchange impulse (window ledger, ``total.impulse``) or the uniform ``read_force``,
+with x along the inflow and y up from the bottom wall. Cd must be positive, and the
+official Cl is positive (centre 0.20 below the centreline 0.205), so a flipped sign
+shows up as Cl outside its interval.
+
+Sampling. One sample every SAMPLE_TSTAR (t* = U t / D; 0.25 s), taken at the end of
+a completed root step with all levels synchronized; the step index is recorded. Cd
+and Cl are means over the sample interval (window exchange ledger refined, per-batch
+forces uniform), which filters acoustic and step-parity oscillation; p_diff is the
+end-of-interval value at the exact surface points.
+
+Steady rule (our engineering choice, frozen before any run; not official). Over the
+trailing WINDOW_TSTAR (t* 50..60 inclusive, 21 samples), for each of Cd, Cl and
+p_diff: range (max - min) <= 10 % and |least-squares slope * window| <= 2 % of the
+width of its official interval. Absolute limits, so the small Cl is never a divisor;
+the slope term bounds slow monotone drift that the range alone would pass. If the
+rule fails at TSTAR_MAX the run reports "not steady" and gives no interval verdict;
+limits are never relaxed after a run. The verdict uses the window mean, and the
+distance of each mean to its nearest interval edge is reported (negative outside).
+"""
+
+import argparse
+import json
+import math
+import os
+from pathlib import Path
+import sys
+import time
+from tempfile import TemporaryDirectory
+
+import numpy as np
+
+HERE = Path(__file__).resolve().parent
+sys.path.insert(0, str(HERE))
+import run_refined_dfg as rr  # noqa: E402
+
+RE = 20.0
+U_PHYS = 0.2                   # m/s mean inflow; rho = 1
+TSTAR_MAX = 60.0               # 30 s physical
+SAMPLE_TSTAR = 0.5             # 0.25 s
+WINDOW_TSTAR = 10.0            # trailing 5 s, 21 samples
+# Schaefer & Turek (1996), 2D-1 lower/upper bounds.
+BANDS = dict(cd=(5.57, 5.59), cl=(0.0104, 0.0110), p_diff=(0.1172, 0.1176))
+STEADY_RANGE = {k: (hi - lo) / 10 for k, (lo, hi) in BANDS.items()}
+STEADY_SLOPE = {k: (hi - lo) / 50 for k, (lo, hi) in BANDS.items()}
+DEVIATIONS = ("pressure-NEQ outlet (neq_extrap), not do-nothing", "local Zou-He inlet")
+
+_re100_config = rr.config
+
+
+def config(h, diameter, collision="SRT", jet_scenario="none"):
+    """The pinned Re100 configuration with the lattice viscosity set for Re 20."""
+    if jet_scenario != "none":
+        raise ValueError("DFG 2D-1 runs have no Jets")
+    cfg, body, p = _re100_config(h, diameter, collision, jet_scenario)
+    nu = p["mean_speed"] * p["diameter"] / RE
+    cfg["physics"]["viscosity"] = nu
+    return cfg, body, dict(p, viscosity=nu, re=RE)
+
+
+def lattice_parameters(diameter, speed=rr.SPEED):
+    nu = speed * diameter / RE
+    return dict(viscosity=nu, tau=0.5 + 3 * nu, steps=math.ceil(TSTAR_MAX * diameter / speed),
+                interval=round(SAMPLE_TSTAR * diameter / speed))
+
+
+def sample(impulse, n, diameter, dp):
+    den = rr.SPEED**2 * diameter / 2
+    return dict(cd=float(impulse[0] / n / den), cl=float(impulse[1] / n / den),
+                dp_scaled=float(dp), p_diff=float(dp) * U_PHYS**2)
+
+
+def summarize(rows):
+    """Steady rule over the trailing window, then the interval verdict (or none)."""
+    if not rows or rows[-1]["tstar"] < TSTAR_MAX - 1e-9:
+        return dict(status="INCOMPLETE", steady=False, verdict=None)
+    window = [r for r in rows if r["tstar"] >= TSTAR_MAX - WINDOW_TSTAR - 1e-9]
+    t = np.array([r["tstar"] for r in window])
+    ranges = {k: max(r[k] for r in window) - min(r[k] for r in window) for k in BANDS}
+    drifts = {k: abs(np.polyfit(t, [r[k] for r in window], 1)[0]) * WINDOW_TSTAR for k in BANDS}
+    steady = all(ranges[k] <= STEADY_RANGE[k] and drifts[k] <= STEADY_SLOPE[k] for k in BANDS)
+    final = {k: float(np.mean([r[k] for r in window])) for k in BANDS}
+    edge = {k: min(final[k] - BANDS[k][0], BANDS[k][1] - final[k]) for k in BANDS}
+    result = dict(window_samples=len(window), ranges=ranges, drifts=drifts,
+                  range_limits=STEADY_RANGE, slope_limits=STEADY_SLOPE, steady=steady,
+                  window_mean=final, edge_distance=edge, bands=BANDS,
+                  deviations=DEVIATIONS, cd_positive=final["cd"] > 0)
+    if not steady:
+        return dict(result, status="NOT_STEADY", verdict=None)
+    inside = {k: BANDS[k][0] <= final[k] <= BANDS[k][1] for k in BANDS}
+    return dict(result, status="WITHIN_INTERVALS" if all(inside.values()) else "OUTSIDE_INTERVALS",
+                verdict=inside)
+
+
+def healthy(row):
+    return row["max_mach"] <= 0.2 and row["rho_min"] >= 0.5 and row["rho_max"] <= 1.5
+
+
+def run_case(a, h, case):
+    from contextlib import ExitStack
+    from CelerisLab.cuda.context import CudaContext
+    from CelerisLab.lbm.descriptors import FLUID
+    from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+    from CelerisLab.simulation import Simulation
+
+    refined = case == "refined"
+    cfgdict, bodydict, p = config(h, 40 if refined else 80, a.collision)
+    lp = lattice_parameters(p["diameter"])
+    rows = []
+    with ExitStack() as stack:
+        if refined:
+            cfg, pair, body, host = rr.host_geometry(a.source, h, a.collision)
+            context = stack.enter_context(CudaContext(0))
+            backend = stack.enter_context(FixedPatchStepper(
+                cfg, rr.PLACEMENTS["wide"], context=context, coarse_boundary="channel",
+                base_block=rr.BASE_BLOCK, fine_block=rr.FINE_BLOCK, fine_objects=(body,),
+                observe_raw_exchange=True, observe_window_exchange=True))
+            from CelerisLab.lbm.channel import channel_masks
+            from CelerisLab.lbm.refinement import ACTIVE
+            from CelerisLab.lbm.descriptors import SOLID
+            masks = [channel_masks(pair.coarse_shape)[0] & ~pair.covered,
+                     (pair.fine_kind == ACTIVE)
+                     & ((host.flags.reshape(pair.fine_shape) & SOLID) == 0)]
+            backend.initialize(rr.initial_state(pair.coarse_shape, (0.0, 0.0), 1.0, root=True),
+                               rr.initial_state(pair.fine_shape, pair.fine_origin, 0.5))
+            probe_p = dict(p, center=[(80.5 - pair.fine_origin[0]) / 0.5 - 0.5,
+                                      (81.0 - pair.fine_origin[1]) / 0.5 - 0.5], radius=40.0)
+        else:
+            tmp = Path(stack.enter_context(TemporaryDirectory()))
+            (tmp / "lbm.json").write_text(json.dumps(cfgdict))
+            (tmp / "body.json").write_text(json.dumps(bodydict))
+            sim = stack.enter_context(Simulation(str(tmp / "lbm.json"), str(tmp / "body.json"),
+                                                 device_id=0))
+            sim.initialize()
+            probe_p = p
+            fluid = (sim.field.flag.reshape(p["ny"], p["nx"]) & FLUID) != 0
+        begin = 0
+        while begin < lp["steps"]:
+            rr.check_deadline(a.deadline)
+            n = min(lp["interval"], lp["steps"] - begin)
+            if refined:
+                impulse = rr.advance_refined_interval(backend, n, a.deadline, "window")
+                fields = [rr.macros(g) for g in backend.download()]
+                rho, ux, uy = fields[1]
+                fluid = masks[1]
+                vals = [tuple(v[m] for v in f) for m, f in zip(masks, fields)]
+            else:
+                impulse, done = np.zeros(2), 0
+                while done < n:
+                    batch = min(100, n - done)
+                    sim.run(batch)
+                    impulse += sim.read_force(0).astype(np.float64) * batch
+                    done += batch
+                m = sim.get_macroscopic()
+                rho, ux, uy = [np.asarray(m[k], np.float64).reshape(p["ny"], p["nx"])
+                               for k in ("rho", "ux", "uy")]
+                vals = [(rho[fluid], ux[fluid], uy[fluid])]
+            begin += n
+            row = dict(step=begin, tstar=begin * rr.SPEED / p["diameter"],
+                       **sample(impulse, n, p["diameter"],
+                                h.pressure_probes(rho, fluid, probe_p)["dp"]),
+                       rho_min=float(min(v[0].min() for v in vals)),
+                       rho_max=float(max(v[0].max() for v in vals)),
+                       max_mach=float(max(np.sqrt(3 * (v[1]**2 + v[2]**2)).max() for v in vals)))
+            rr.ensure_finite(row)
+            rows.append(row)
+            if not healthy(row):
+                break
+    summary = summarize(rows) if healthy(rows[-1]) else dict(status="UNHEALTHY", steady=False,
+                                                             verdict=None)
+    return dict(case=case, collision=a.collision, lattice=lp, rows=rows, **summary)
+
+
+def plan(collision="SRT"):
+    return dict(re=RE, collision=collision, tstar_max=TSTAR_MAX, sample_tstar=SAMPLE_TSTAR,
+                window_tstar=WINDOW_TSTAR, bands=BANDS, range_limits=STEADY_RANGE, slope_limits=STEADY_SLOPE,
+                deviations=DEVIATIONS, pressure_scale=U_PHYS**2,
+                refined_root=lattice_parameters(40), refined_fine_tau=0.5 + 3 * 2 * rr.SPEED * 40 / RE,
+                uniform=lattice_parameters(80),
+                label="compared against the official reference bands (non-canonical neq_extrap outlet, zou_he_local inlet); not a full benchmark certification")
+
+
+def main():
+    parser = argparse.ArgumentParser(description=__doc__.splitlines()[0])
+    parser.add_argument("--case", choices=["refined", "uniform", "both"], default="refined")
+    parser.add_argument("--collision", choices=["SRT", "MRT"], default="SRT")
+    parser.add_argument("--source", type=Path, default=rr.DEFAULT_SOURCE)
+    parser.add_argument("--output", type=Path, required=True)
+    parser.add_argument("--plan", action="store_true", help="Write the frozen plan; no CUDA")
+    parser.add_argument("--max-wall-seconds", type=float)
+    a = parser.parse_args()
+    a.source = a.source.resolve()
+    a.output.mkdir(parents=True, exist_ok=False)
+    (a.output / "plan.json").write_text(json.dumps(plan(a.collision), indent=2, allow_nan=False))
+    if a.plan:
+        return
+    if not a.max_wall_seconds or a.max_wall_seconds <= 0 or not os.environ.get("CUDA_VISIBLE_DEVICES"):
+        raise SystemExit("GPU run needs --max-wall-seconds > 0 and explicit CUDA_VISIBLE_DEVICES")
+    h, _ = rr.helpers(a.source)
+    rr.config = config  # host_geometry builds its root configuration through this name
+    a.deadline = time.monotonic() + a.max_wall_seconds
+    cases = ["refined", "uniform"] if a.case == "both" else [a.case]
+    results = [run_case(a, h, case) for case in cases]
+    (a.output / "result.json").write_text(json.dumps(results, indent=2, allow_nan=False))
+
+
+if __name__ == "__main__":
+    main()