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>

tests/validation/multigrid/test_refinement_blocked_base.py

@@ -2,6 +2,37 @@
 
 Run CPU and GPU classes in separate processes. GPU callers must first check
 resources/device availability; this file never selects or reserves a GPU.
+
+Dense/blocked route parity (DenseBlockedRouteParity): the same configuration,
+initial state, body/Jet schedule and tilings on three storage routes, only
+level_blocks varying: D/D [None, None], D/B [None, FINE_BLOCK] (dense root,
+the mixed pattern of current examples) and B/B [(12, 8), FINE_BLOCK]. Gate:
+every public and owner array identical in FP32 bits and NaN masks at steps
+0, 1, 2, 8 and after re-initialize; D/B and B/B raw-exchange ledgers equal.
+No tolerance. A mismatch is reported with route pair/level/step and diagnosed
+(boundary semantics vs arithmetic ordering), never widened.
+
+Route capabilities (rejections checked before CUDA by test_route_capability_matrix):
+
+  combination                                   D/D       D/B       B/B
+  periodic / slot_reservoir root, fine circle   yes       yes       yes   (parity gate)
+  actuated circle, velocity Jet + amplitude     yes       yes       yes   (parity gate)
+  raw/window exchange on a body level           rejected  yes       yes
+  shifted storage                               rejected  rejected  yes
+  projected explosion                           rejected  rejected  yes
+  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. MRT body-free is also covered by
+test_refinement_mrt_gpu::test_dense_versus_blocked_mrt_two_level. Parked:
+shifted and depth >= 3 (no dense route to compare).
 """
 from contextlib import ExitStack, contextmanager, nullcontext
 from dataclasses import replace
@@ -20,6 +51,7 @@ from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
 from CelerisLab.lbm.refinement import ACTIVE, OUTER_GHOST, FixedPatch2D
 from CelerisLab.lbm.refinement_stepper import FixedPatchStepper, blocked_base_tables
 from test_refinement_blocked_fine import NX, NY, BOUNDS, config, state, circle, actuation
+from test_refinement_jets import SPEED, jet
 
 BASE_BLOCKS = ((24, 32), (12, 8), (3, 4), (1, 1))
 FINE_BLOCK = (16, 20)
@@ -324,5 +356,167 @@ class BlockedBaseGPU(unittest.TestCase):
                             print(f"PASS {label}: full finite bits/NaN masks at 0,1,2,8 + reset", flush=True)
 
 
+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):
+    """Exact equality of nested exchange summaries (dicts, sequences, arrays, numbers)."""
+    if isinstance(a, dict):
+        return a.keys() == b.keys() and all(same_summary(a[k], b[k]) for k in a)
+    if isinstance(a, (list, tuple)):
+        return len(a) == len(b) and all(same_summary(x, y) for x, y in zip(a, b))
+    a, b = np.asarray(a), np.asarray(b)
+    if a.dtype.kind in "biufc" and b.dtype.kind in "biufc":
+        return a.shape == b.shape and np.array_equal(a, b, equal_nan=True)
+    return a.dtype == b.dtype and np.array_equal(a, b)  # labels such as the window name
+
+
+class DenseBlockedRouteCPU(unittest.TestCase):
+    def test_same_summary_compares_labels_and_every_number(self):
+        ledger = {"window": "last_base_step", "time_start": 3, "total": {"force": [1.0, -2.0],
+                  "torque": 0.5}, "aperture": {"force": np.array([np.nan, 1.0])}}
+        self.assertTrue(same_summary(ledger, replace_value(ledger, ())))
+        for path, value in ((("window",), "last_substep"), (("time_start",), 4),
+                            (("total", "force"), [1.0, -2.000001]), (("total", "torque"), 0.25),
+                            (("aperture", "force"), np.array([0.0, 1.0]))):
+            with self.subTest(path=path):
+                self.assertFalse(same_summary(ledger, replace_value(ledger, path, value)))
+                self.assertFalse(same_summary(replace_value(ledger, path, value), ledger))
+
+    def test_route_capability_matrix(self):
+        context = Mock(activate=Mock(side_effect=AssertionError("CUDA reached")))
+        periodic = config("periodic")
+        inner = (BOUNDS[0] + 2, BOUNDS[1] + 2, BOUNDS[2] - 2, BOUNDS[3] - 2)
+        second = (BOUNDS[0], BOUNDS[1], BOUNDS[0] + 2, BOUNDS[1] + 2)
+        cases = (
+            ("raw exchange on dense body level", periodic, dict(
+                level_blocks=[None, None], level_objects=[(), (circle(),)], observe_raw_exchange=True)),
+            ("shifted D/D", replace(periodic, ddf_shifting=True), dict(level_blocks=[None, None])),
+            ("shifted D/B", replace(periodic, ddf_shifting=True), dict(level_blocks=[None, FINE_BLOCK])),
+            ("projected D/D", periodic, dict(level_blocks=[None, None], explosion="projected_linear")),
+            ("projected D/B", periodic, dict(level_blocks=[None, FINE_BLOCK], explosion="projected_linear")),
+            ("channel D/B", periodic, dict(level_blocks=[None, FINE_BLOCK], coarse_boundary="channel")),
+            ("several L1 rectangles", periodic, dict(refinements=[[BOUNDS, second]], level_blocks=[None, None])),
+            ("depth 3 dense", periodic, dict(refinements=[BOUNDS, inner], level_blocks=[None, None, None])),
+        )
+        for label, cfg, kw in cases:
+            kw.setdefault("refinements", [BOUNDS])
+            with self.subTest(label), self.assertRaises(ValueError):
+                FixedPatchStepper(cfg, context=context, **kw)
+        context.activate.assert_not_called()
+
+
+def replace_value(tree, path, value=None):
+    """Deep copy of a nested dict with the item at ``path`` replaced (no change for ())."""
+    if isinstance(tree, dict):
+        return {k: (value if path == (k,) else replace_value(v, path[1:], value) if path[:1] == (k,)
+                    else replace_value(v, (), None)) for k, v in tree.items()}
+    return np.array(tree, copy=True) if isinstance(tree, np.ndarray) else tree
+
+
+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, cfg=None, bounds=BOUNDS, routes=ROUTES):
+        from CelerisLab.cuda.context import CudaContext
+        cfg = cfg or config(mode)
+        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:
+                observe = bool(objects) and blocks[1] is not None
+                snapshots, ledgers = [], []
+                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:
+                    backend.initialize(*arrays)
+                    for step in (0, 1, 2, 8):
+                        if step == 8 and switch is not None:
+                            backend.set_jet_amplitude(0, switch)
+                        backend.step(step - backend.step_count)
+                        snapshots.append(backend.download() + raw(backend))
+                        if observe and step:
+                            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))
+                results[name] = snapshots, ledgers
+        return results
+
+    def check(self, label, results):
+        reference, _ = results["D/D"]
+        for name in ("D/B", "B/B"):
+            for step, actual, expected in zip((0, 1, 2, 8, "reset"), results[name][0], reference):
+                for index, (a, b) in enumerate(zip(actual, expected)):
+                    exact(self, a, b, f"{label} {name} vs D/D step={step} array={index}")
+        if results["D/B"][1]:
+            for step, a, b in zip((1, 2, 8), results["D/B"][1], results["B/B"][1]):
+                self.assertTrue(same_summary(a, b), f"{label} raw exchange D/B vs B/B step={step}")
+        for name, (snapshots, _) in results.items():
+            masses = [float(np.nansum(level.astype(np.float64))) for level in snapshots[3][:2]]
+            print(f"{label} {name}: step 8 mass by level {masses}", flush=True)
+        print(f"PASS {label}: D/D = D/B = B/B bits/NaN masks at 0,1,2,8 + reset", flush=True)
+
+    def test_body_free(self):
+        for mode in ("periodic", "slot_reservoir"):
+            self.check(f"{mode} body-free", self.run_routes(mode))
+
+    def test_passive_circle(self):
+        for mode in ("periodic", "slot_reservoir"):
+            self.check(f"{mode} circle", self.run_routes(mode, objects=[circle()]))
+
+    def test_actuated_circle(self):
+        for mode in ("periodic", "slot_reservoir"):
+            for speed in (SPEED, -SPEED):
+                self.check(f"{mode} actuated {speed}",
+                           self.run_routes(mode, objects=[circle(actuation(speed))]))
+
+    def test_velocity_jet_with_amplitude_switch(self):
+        for mode in ("periodic", "slot_reservoir"):
+            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()