test: dense/blocked parity rows for MRT bodies and two fine bodies

Extends DenseBlockedRouteParity (test-only, on 1573942) with the residual
dense-supported combinations named before any dense-route deletion:
- test_mrt_bodies: MRT on a periodic root with a passive circle, an
  actuated circle (+/-) and a velocity Jet with an amplitude switch
  (MRT slot/reservoir is rejected on every route: on_site_velocity is
  SRT-only);
- test_two_bodies: two circles on a wider patch (4, 4, 28, 20), periodic
  and slot/reservoir, Jet on one + passive other, and two actuated (+/-);
  routes D/D, D/B [None, (8, 8)], B/B [(12, 8), (8, 8)].
Same gate: exact FP32 bits and NaN masks at 0, 1, 2, 8 and after
re-initialize; D/B vs B/B raw-exchange ledgers equal for every body.
Docstring table updated; no tolerance, no src change.

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

examples/fixed_patch_circle.py

@@ -1,9 +1,9 @@
 """Public API demo: fixed 2:1 patch with a stationary, optionally actuated circle.
 
 From this checkout, after checking the selected physical GPU is idle:
-    PYTHONPATH="$PWD/src" CUDA_VISIBLE_DEVICES=0 python examples/fixed_patch_circle.py
+    PYTHONPATH="$PWD/src" CUDA_VISIBLE_DEVICES=0 python examples/fixed_patch_circle.py --fine-block 16 20
-    ... --speed -0.002 --mode slot_reservoir
+    ... --speed -0.002 --fine-block 16 20 --mode slot_reservoir
-Both levels use blocked storage (--base-block, --fine-block tiles).
+Omit --fine-block for the (equivalent) dense fine storage.
 
 A short smoke run (default 128 coarse steps): it shows setup, stepping and
 fluid-only diagnostics, not a developed flow, accuracy or long-time stability.
@@ -71,10 +71,8 @@ def main():
     parser.add_argument("--mode", choices=("periodic", "slot_reservoir"), default="periodic")
     parser.add_argument("--speed", type=float, default=.002,
                         help="normal surface velocity, coarse lattice units (+ blows, - sucks, 0 passive)")
-    parser.add_argument("--base-block", type=int, nargs=2, metavar=("BY", "BX"), default=(12, 8),
+    parser.add_argument("--fine-block", type=int, nargs=2, metavar=("BY", "BX"),
-                        help="root tile; must divide the 24 x 32 root")
+                        help="blocked fine storage; must divide the whole fine allocation (32 x 40 here)")
-    parser.add_argument("--fine-block", type=int, nargs=2, metavar=("BY", "BX"), default=(16, 20),
-                        help="fine tile; must divide the whole fine allocation (32 x 40 here)")
     parser.add_argument("--steps", type=int, default=128)
     args = parser.parse_args()
     if args.steps < 0:
@@ -88,20 +86,20 @@ def main():
     body = SimObject(0, geometry=CircleGeometry(*CENTER, RADIUS), center=CENTER, radius=RADIUS,
                      surface_actuation=spec)
     cfg = config(args.mode)
+    fine_block = tuple(args.fine_block) if args.fine_block else None
     # The caller owns the CUDA context; the stepper must be closed first (the
     # nested with-statement does that).
     with CudaContext(0) as context, FixedPatchStepper(
             cfg, context=context, coarse_boundary=args.mode, refinements=[BOUNDS],
-            level_blocks=[tuple(args.base_block), tuple(args.fine_block)],
+            level_blocks=[None, fine_block], level_objects=[(), (body,)]) as run:
-            level_objects=[(), (body,)]) as run:
         p = run.patch
         # Rest equilibrium is already postcollision G, so no initial collision is
         # needed. Non-leaf entries (covered, ghost, solid) are discarded by initialize.
         rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None]
         run.initialize(np.broadcast_to(rest, (9, *p.coarse_shape)),
                        np.broadcast_to(rest, (9, *p.fine_shape)))
-        storage = f"blocked root {tuple(args.base_block)} fine {tuple(args.fine_block)}"
+        storage = f"blocked {fine_block}" if fine_block else "dense"
-        print(f"mode={args.mode} storage={storage} speed={args.speed:+g} "
+        print(f"mode={args.mode} fine storage={storage} speed={args.speed:+g} "
               f"(fine allocation {p.fine_shape[0]}x{p.fine_shape[1]}, fine origin {p.fine_origin})")
         masks = leaf_masks(run, args.mode)
         start = leaf_diagnostics(*run.download(), masks)

examples/fixed_patch_jet.py

@@ -31,7 +31,6 @@ from CelerisLab.lbm.refinement import ACTIVE
 from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
 
 BOUNDS, AMPLITUDE, WIDTH = (16, 4, 32, 28), .02, 12
-LEVEL_BLOCKS = [(12, 8), (8, 8)]  # blocked storage: 48x48 root, 56x40 fine window
 TIMES, MAX_MACH, MAX_SECONDS = (0, 300, 1200), .1, 120
 
 
@@ -133,8 +132,7 @@ def main():
                   lattice_model='D2Q9', collision='SRT', streaming='double_buffer',
                   store_precision='FP32', ddf_shifting=False, les_enabled=False)
     report = dict(config=dict(nx=48, ny=48, bounds=BOUNDS, W=WIDTH, A=AMPLITUDE, nu=.1,
-                              rho=1, tau_coarse=.8, tau_fine=1.1, coarse_boundary='slot_reservoir',
+                              rho=1, tau_coarse=.8, tau_fine=1.1, coarse_boundary='slot_reservoir'),
-                              level_blocks=LEVEL_BLOCKS),
                   guards=dict(check_every=50, positive_density=True, finite=True,
                               max_Mach=MAX_MACH, max_seconds=MAX_SECONDS,
                               timing='cooperative checks, not an exact interruption deadline'),
@@ -151,8 +149,7 @@ def main():
     start, run = time.monotonic(), None
     try:
         with CudaContext(0) as context, FixedPatchStepper(cfg, context=context,
-                coarse_boundary='slot_reservoir', refinements=[BOUNDS],
+                coarse_boundary='slot_reservoir', refinements=[BOUNDS]) as run:
-                level_blocks=LEVEL_BLOCKS) as run:
             p = run.patch
             masks = [np.ones(p.coarse_shape, bool), p.fine_kind == ACTIVE]
             masks[0][[0, -1]] = False

tests/unit/test_example_tiles.py

@@ -1,65 +0,0 @@
-"""The public examples run all-blocked: their level tiles must reach the device unchanged.
-
-Host only: a context whose activation raises proves construction (layouts, tables,
-boundary checks) succeeded up to the first CUDA call.
-"""
-import ast
-from dataclasses import replace
-from pathlib import Path
-import unittest
-from unittest.mock import Mock
-
-from CelerisLab.body.geometry.circle import CircleGeometry
-from CelerisLab.body.objects import SimObject
-from CelerisLab.config import LBMConfig
-from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
-
-EXAMPLES = Path(__file__).resolve().parents[2] / "examples"
-
-
-def literal(name, constant):
-    """Module-level literal from an example, read without importing it (no pycuda)."""
-    tree = ast.parse((EXAMPLES / name).read_text())
-    for node in tree.body:
-        if isinstance(node, ast.Assign):
-            for target in node.targets:
-                names = target.elts if isinstance(target, ast.Tuple) else [target]
-                for i, t in enumerate(names):
-                    if isinstance(t, ast.Name) and t.id == constant:
-                        value = node.value.elts[i] if isinstance(target, ast.Tuple) else node.value
-                        return ast.literal_eval(value)
-    raise KeyError(constant)
-
-
-def forbidden():
-    return Mock(activate=Mock(side_effect=AssertionError("CUDA reached")), sm_arch="sm_70")
-
-
-class ExampleTiles(unittest.TestCase):
-    def test_fixed_patch_jet_tiles(self):
-        bounds, blocks = literal("fixed_patch_jet.py", "BOUNDS"), literal("fixed_patch_jet.py", "LEVEL_BLOCKS")
-        width, amplitude = literal("fixed_patch_jet.py", "WIDTH"), literal("fixed_patch_jet.py", "AMPLITUDE")
-        cfg = LBMConfig(nx=48, ny=48, viscosity=.1, velocity=0., wall_slot_width=width,
-                        wall_normal_velocity=amplitude, y_wall_bc="on_site_velocity",
-                        top_wall_bc="reservoir")
-        with self.assertRaisesRegex(AssertionError, "CUDA reached"):
-            FixedPatchStepper(cfg, context=forbidden(), coarse_boundary="slot_reservoir",
-                              refinements=[bounds], level_blocks=[tuple(b) for b in blocks])
-
-    def test_fixed_patch_circle_default_tiles(self):
-        bounds = literal("fixed_patch_circle.py", "BOUNDS")
-        nx, ny = literal("fixed_patch_circle.py", "NX"), literal("fixed_patch_circle.py", "NY")
-        center, radius = literal("fixed_patch_circle.py", "CENTER"), literal("fixed_patch_circle.py", "RADIUS")
-        body = SimObject(0, geometry=CircleGeometry(*center, radius), center=center, radius=radius)
-        periodic = LBMConfig(nx=nx, ny=ny, viscosity=.1)
-        slot = replace(periodic, velocity=0, wall_normal_velocity=.001, wall_slot_width=8,
-                       y_wall_bc="on_site_velocity", top_wall_bc="reservoir")
-        for mode, cfg in (("periodic", periodic), ("slot_reservoir", slot)):
-            with self.subTest(mode=mode), self.assertRaisesRegex(AssertionError, "CUDA reached"):
-                FixedPatchStepper(cfg, context=forbidden(), coarse_boundary=mode,
-                                  refinements=[bounds], level_blocks=[(12, 8), (16, 20)],
-                                  level_objects=[(), (body,)])
-
-
-if __name__ == "__main__":
-    unittest.main()

tests/validation/multigrid/test_refinement_blocked_base.py

@@ -23,13 +23,16 @@ Route capabilities (rejections checked before CUDA by test_route_capability_matr
   channel root                                  rejected  rejected  yes
   several L1 rectangles                         rejected  rejected  yes
   depth >= 3 with a dense level                 rejected  rejected  yes
+  MRT, periodic root, circle/actuation/Jet      yes       yes       yes   (parity gate)
+  two fine bodies (Jet + passive, two actuated) yes       yes       yes   (parity gate)
+  MRT with slot/reservoir root                  rejected on every route (on_site_velocity is SRT-only)
 Slot/reservoir walls are compile-time WALL_* macros on a dense root and the
 d_block_wall constant on a blocked root: BC-encoding parity, bit-identical by
 this gate. Dense fine bodies use CurvedLinkSoA + CurvedBoundaryKernel
 (n_objects > 0), blocked ones BlockCurvedBoundary; same links and physics,
-compared through the population rows. Parked here: MRT (covered by
+compared through the population rows. MRT body-free is also covered by
-test_refinement_mrt_gpu::test_dense_versus_blocked_mrt_two_level), shifted and
+test_refinement_mrt_gpu::test_dense_versus_blocked_mrt_two_level. Parked:
-depth >= 3 (no dense route to compare), several bodies (later slice).
+shifted and depth >= 3 (no dense route to compare).
 """
 from contextlib import ExitStack, contextmanager, nullcontext
 from dataclasses import replace
@@ -354,6 +357,25 @@ class BlockedBaseGPU(unittest.TestCase):
 
 
 ROUTES = (("D/D", [None, None]), ("D/B", [None, FINE_BLOCK]), ("B/B", [(12, 8), FINE_BLOCK]))
+# Two-body rows: patch (4, 4, 28, 20) gives a 40 x 56 fine window, tiled (8, 8).
+WIDE_BOUNDS = (4, 4, 28, 20)
+WIDE_ROUTES = (("D/D", [None, None]), ("D/B", [None, (8, 8)]), ("B/B", [(12, 8), (8, 8)]))
+LEFT, RIGHT, RADIUS = (12.0, 12.0), (20.0, 12.0), 2.15
+
+
+def body(ident, center, spec=None):
+    from CelerisLab.body.geometry.circle import CircleGeometry
+    from CelerisLab.body.objects import SimObject
+    return SimObject(ident, geometry=CircleGeometry(*center, RADIUS), center=center, radius=RADIUS,
+                     surface_actuation=spec)
+
+
+TWO_BODY_CASES = (
+    ("jet + passive", lambda: [body(0, LEFT), body(1, RIGHT)],
+     lambda: [actuation(SPEED).as_jet(LEFT, RADIUS)], -SPEED),
+    ("actuated +/-", lambda: [body(0, LEFT, actuation(SPEED)), body(1, RIGHT, actuation(-SPEED))],
+     lambda: [], None),
+)
 
 
 def same_summary(a, b):
@@ -414,17 +436,18 @@ def replace_value(tree, path, value=None):
 class DenseBlockedRouteParity(unittest.TestCase):
     """D/D vs D/B vs B/B, unshifted SRT, one circle or one velocity Jet; see module docstring."""
 
-    def run_routes(self, mode, objects=(), jets=(), switch=None):
+    def run_routes(self, mode, objects=(), jets=(), switch=None, cfg=None, bounds=BOUNDS, routes=ROUTES):
         from CelerisLab.cuda.context import CudaContext
-        cfg = config(mode)
+        cfg = cfg or config(mode)
-        p = FixedPatch2D.from_config(cfg, BOUNDS)
+        p = FixedPatch2D.from_config(cfg, bounds)
         arrays = state(NY, NX, 1), state(*p.fine_shape, 2)
+        bodies = sorted({obj.obj_id for obj in objects})
         results = {}
         with CudaContext(0) as context:
-            for name, blocks in ROUTES:
+            for name, blocks in routes:
                 observe = bool(objects) and blocks[1] is not None
                 snapshots, ledgers = [], []
-                with FixedPatchStepper(cfg, context=context, coarse_boundary=mode, refinements=[BOUNDS],
+                with FixedPatchStepper(cfg, context=context, coarse_boundary=mode, refinements=[bounds],
                                        level_blocks=blocks, level_objects=[(), tuple(objects)],
                                        level_jets=[(), tuple(jets)],
                                        observe_raw_exchange=observe) as backend:
@@ -435,7 +458,8 @@ class DenseBlockedRouteParity(unittest.TestCase):
                         backend.step(step - backend.step_count)
                         snapshots.append(backend.download() + raw(backend))
                         if observe and step:
-                            ledgers.append(tuple(backend.read_raw_exchange(0, window=window)
+                            ledgers.append(tuple(backend.read_raw_exchange(body, window=window)
+                                                 for body in bodies
                                                  for window in ("last_base_step", "last_substep")))
                     backend.initialize(*arrays)
                     snapshots.append(backend.download() + raw(backend))
@@ -475,6 +499,24 @@ class DenseBlockedRouteParity(unittest.TestCase):
             self.check(f"{mode} jet switch",
                        self.run_routes(mode, objects=[circle()], jets=[jet()], switch=-SPEED))
 
+    def test_mrt_bodies(self):
+        """MRT, periodic root (on_site_velocity slot/reservoir is SRT-only on every route)."""
+        cfg = replace(config("periodic"), collision="MRT")
+        for label, objects, jets, switch in (("circle", [circle()], (), None),
+                                             ("actuated +", [circle(actuation(SPEED))], (), None),
+                                             ("actuated -", [circle(actuation(-SPEED))], (), None),
+                                             ("jet switch", [circle()], [jet()], -SPEED)):
+            self.check(f"MRT periodic {label}",
+                       self.run_routes("periodic", objects, jets, switch, cfg=cfg))
+
+    def test_two_bodies(self):
+        """Two circles on a wider patch (same tiles on every route): Jet + passive, two actuated."""
+        for mode in ("periodic", "slot_reservoir"):
+            for label, objects, jets, switch in TWO_BODY_CASES:
+                self.check(f"{mode} two bodies {label}",
+                           self.run_routes(mode, objects(), jets(), switch, bounds=WIDE_BOUNDS,
+                                           routes=WIDE_ROUTES))
+
 
 if __name__ == "__main__":
     unittest.main()