test(D1): exchange, Jet and fine-surface checks run on all-blocked storage

- test_exchange_window and test_raw_exchange_gpu fine(): dense root becomes
  a (12, 8) blocked root (fine(None) stays only for the dense reference test
  that D2 retires);
- jets: G1/G3 compare blocked layouts against the FINE_BLOCKS[0] reference
  (blocked root); the dense run is a separate leg D2 removes; G2 legacy
  actuation vs Jet on the blocked reference; CPU fixtures use a blocked root;
- test_refinement_blocked_fine.raw reads every level through its layout;
- test_fine_surface_gpu: new test_fine_surface_blocked_matches_uniform_four_hops
  (fine ROI after each of four hops vs the uniform 48^2 solver with the same
  tolerances, blowing/suction mass sign from raw exchange, zero actuation ==
  passive bit for bit); the cross-interface donor rejection passes tiles.
  The dense hook/ledger/audit probes in test_fine_surface_four_hops and the
  dense SoA fault tests stay until D2 retires them by name.

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

tests/validation/multigrid/test_exchange_window.py

@@ -109,7 +109,7 @@ class ExchangeWindowGPU(unittest.TestCase):
         from test_raw_exchange_gpu import BOUNDS as RB, CFG, circle as rc, jet, state
         from CelerisLab.cuda.context import CudaContext
         with CudaContext(0) as ctx:
-            kw = dict(context=ctx, refinements=[RB], level_blocks=[None, (16, 20)],
+            kw = dict(context=ctx, refinements=[RB], level_blocks=[(12, 8), (16, 20)],
                       level_objects=[(), (rc(),)], level_jets=[(), (jet(),)],
                       observe_raw_exchange=True)
             with FixedPatchStepper(CFG, observe_window_exchange=True, **kw) as windowed, \

tests/validation/multigrid/test_fine_surface_gpu.py

@@ -279,6 +279,61 @@ def test_fine_surface_four_hops(speed):
             np.testing.assert_array_equal(got, wanted)
 
 
+def run_blocked_patch(speed):
+    """All-blocked run (whole-level tiles) of the same case; fine ROI after each fine hop."""
+    cfg = LBMConfig(nx=24, ny=24, viscosity=.05, velocity=0, rho=1)
+    source = object_spec(speed)
+    with CudaContext(0) as context:
+        with implementation.FixedPatchStepper(cfg, context=context, refinements=[BOUNDS],
+                                              level_blocks=[(24, 24), (40, 40)],
+                                              level_objects=[(), (source,)],
+                                              observe_raw_exchange=True) as backend:
+            c, f = backend._levels
+            flags = backend.fine_boundary.flags.reshape(40, 40)
+            active = np.zeros((40, 40), bool)
+            active[4:-4, 4:-4] = True
+            active &= ((flags & FLUID) != 0) & ((flags & SOLID) == 0)
+            yy, xx = np.mgrid[:40, :40]
+            roi = active & ((xx-20)**2+(yy-20)**2 <= 9**2)
+            backend.initialize(rest((24, 24)), rest((40, 40)))
+            traces, original = [], f.kernel
+
+            def fine(*args, **kwargs):
+                original(*args, **kwargs)
+                backend._stream.synchronize()
+                advanced = f.layout.to_dense(read(args[1], f.layout.blocked_shape))
+                assert np.isfinite(advanced[:, active]).all()
+                assert np.isnan(advanced[:, (flags & SOLID) != 0]).all()
+                traces.append(advanced[:, roi].copy())
+
+            with patch.object(f, "kernel", fine):
+                backend.step(2)
+            assert len(traces) == 4 and backend.step_count == 2
+            mass = backend.read_raw_exchange(0)["total"]["mass_into_fluid"]
+            return traces, roi, backend.download(), float(mass)
+
+
+@pytest.mark.parametrize("speed", [0., .002, -.002])
+def test_fine_surface_blocked_matches_uniform_four_hops(speed):
+    """Blocked storage: fine ROI after each of four hops against the uniform 48^2 solver
+    (same tolerances as the dense test); blowing/suction drives mass with its sign;
+    zero actuation equals the passive circle bit for bit."""
+    actual, roi, final, mass = run_blocked_patch(speed)
+    reference = uniform_reference(speed, roi)
+    for got, wanted in zip(actual, reference):
+        np.testing.assert_allclose(got, wanted, atol=3e-6, rtol=0)
+        np.testing.assert_allclose(MOMENTS@got.astype(float), MOMENTS@wanted.astype(float),
+                                   atol=2e-5, rtol=0)
+    if speed:
+        assert speed*mass > 0
+    else:
+        passive, passive_roi, passive_final, _ = run_blocked_patch(None)
+        np.testing.assert_array_equal(passive_roi, roi)
+        for got, wanted in zip(actual+list(final), passive+list(passive_final)):
+            np.testing.assert_array_equal(got, wanted)
+    print(f"blocked speed={speed} four hops vs uniform PASS mass_into_fluid={mass:+.9g}", flush=True)
+
+
 @pytest.mark.parametrize("missing", ["surface_velocity_gpu", "link_exchange_gpu"])
 def test_missing_surface_device_buffer_rejects_before_launch(missing):
     backend = implementation.FixedPatchStepper.__new__(implementation.FixedPatchStepper)
@@ -328,6 +383,7 @@ def test_cross_interface_donor_rejects_before_cuda():
     with patch.object(cuda, "mem_alloc") as allocate:
         with pytest.raises(ValueError, match="donor crosses"):
             implementation.FixedPatchStepper(cfg, context=context, refinements=[BOUNDS],
+                                             level_blocks=[(24, 24), (40, 40)],
                                              level_objects=[(), (source,)])
         context.activate.assert_not_called()
         allocate.assert_not_called()

tests/validation/multigrid/test_raw_exchange_gpu.py

@@ -92,8 +92,9 @@ def census(run, index, buffer):
 
 
 def fine(block):
-    """Two-level circle+jet construction on BOUNDS; block None keeps the dense fine level."""
-    return dict(refinements=[BOUNDS], level_blocks=[None, block],
+    """Two-level circle+jet construction on BOUNDS, all-blocked with a (12, 8) root tile;
+    block None is the dense route (dense reference test only, retired with the dense route)."""
+    return dict(refinements=[BOUNDS], level_blocks=[None if block is None else (12, 8), block],
                 level_objects=[(), (circle(),)], level_jets=[(), (jet(),)])
 
 class RawExchangeGPU(unittest.TestCase):

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
-test_refinement_mrt_gpu::test_dense_versus_blocked_mrt_two_level), shifted and
-depth >= 3 (no dense route to compare), several bodies (later slice).
+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
@@ -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)
-        p = FixedPatch2D.from_config(cfg, BOUNDS)
+        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:
+            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()

tests/validation/multigrid/test_refinement_blocked_fine.py

@@ -88,15 +88,13 @@ def seam_cells(layout, mapped):
 
 
 def raw(backend):
-    """Unmasked device arrays; blocked fine storage unpacked from owner cores."""
-    coarse, fine = backend._levels
-    out = [np.empty(coarse.shape, np.float32)]
-    out.append(np.empty(fine.shape if backend.fine_layout is None
-                        else backend.fine_layout.blocked_shape, np.float32))
-    cuda.memcpy_dtoh(out[0], coarse.current)
-    cuda.memcpy_dtoh(out[1], fine.current)
-    if backend.fine_layout is not None:
-        out[1] = backend.fine_layout.to_dense(out[1])
+    """Unmasked device arrays of every level; blocked storage unpacked from owner cores."""
+    out = []
+    for level in backend._levels:
+        array = np.empty(level.shape if level.layout is None else level.layout.blocked_shape,
+                         np.float32)
+        cuda.memcpy_dtoh(array, level.current)
+        out.append(array if level.layout is None else level.layout.to_dense(array))
     return out
 
 

tests/validation/multigrid/test_refinement_followup.py

@@ -17,6 +17,7 @@ import pycuda.driver as cuda
 
 import CelerisLab
 from CelerisLab.lbm import refinement_stepper as implementation
+from CelerisLab.lbm.blocks import BlockLayout2D
 from test_refinement_gpu import (ROOT, POP_ATOL, MOMENT_ATOL, ref, postcollision,
                                  totals, LBMConfig, CudaContext, compiler,
                                  PeriodicPatchStepper)
@@ -83,6 +84,21 @@ def scalar_interval(coarse, fine, bounds, slot=None):
     return result, current
 
 
+def whole_tiles(cfg, bounds):
+    """One whole-level tile per level: (ny, nx) at the root, the fine window below."""
+    return [(cfg.ny, cfg.nx), implementation.FixedPatch2D.from_config(cfg, bounds).fine_shape]
+
+
+def owner_arrays(backend):
+    """Raw owner storage of every level as dense (9, ny, nx), before public masking."""
+    arrays = []
+    for level in backend._levels:
+        blocked = np.empty(level.layout.blocked_shape, np.float32)
+        cuda.memcpy_dtoh(blocked, level.current)
+        arrays.append(level.layout.to_dense(blocked))
+    return arrays
+
+
 class FaultGuardsCPU(unittest.TestCase):
     def test_boundary_mode_and_margin_contract(self):
         cfg = LBMConfig(nx=32, ny=24, viscosity=.1, velocity=0, wall_normal_velocity=.001,
@@ -95,7 +111,7 @@ class FaultGuardsCPU(unittest.TestCase):
             with self.subTest(mode=mode, bounds=bounds):
                 with self.assertRaises(ValueError):
                     implementation.FixedPatchStepper(cfg, context=context, coarse_boundary=mode,
-                                                     refinements=[bounds])
+                                                     refinements=[bounds], level_blocks=[(1, 1)] * 2)
                 context.activate.assert_not_called()
         # Both exact channel limits and the original periodic two-cell margin
         # must reach context activation, with no allocation or device required.
@@ -104,18 +120,19 @@ class FaultGuardsCPU(unittest.TestCase):
             with self.subTest(mode=mode, bounds=bounds):
                 with self.assertRaisesRegex(RuntimeError, "validated before CUDA"):
                     implementation.FixedPatchStepper(cfg, context=context, coarse_boundary=mode,
-                                                     refinements=[bounds])
+                                                     refinements=[bounds], level_blocks=[(1, 1)] * 2)
 
     def fake_backend(self):
         backend = PeriodicPatchStepper.__new__(PeriodicPatchStepper)
-        backend._layouts = [None, None]  # dense storage on both levels
+        layout = BlockLayout2D.build(1, 1, 1, 1, periodic=(False, False))
+        backend._layouts = [layout, layout]  # blocked storage on both levels
         backend._ddf_shifting = False  # unshifted on every level, as _check_storage expects
         backend._level_cfgs = [LBMConfig(nx=16, ny=16, viscosity=.1)] * 2  # ddf_shifting off
         backend.ctx = Mock(activate=nullcontext)
         backend._lock = RLock()
         backend._status = "ready"
         backend._stream = Mock()
-        backend._levels = [Mock(shape=(9, 1, 1), layout=None), Mock(shape=(9, 1, 1), layout=None)]
+        backend._levels = [Mock(shape=(9, 1, 1), layout=layout), Mock(shape=(9, 1, 1), layout=layout)]
         backend._allocations = [Mock(), Mock()]
         return backend
 
@@ -158,7 +175,8 @@ class FaultGuardsCPU(unittest.TestCase):
         context = Mock(activate=nullcontext, sm_arch="sm_70")
         with patch.object(implementation, "HALO", 2):
             with self.assertRaisesRegex(RuntimeError, "four-cell"):
-                PeriodicPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)])
+                PeriodicPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)],
+                                     level_blocks=[(16, 16), (24, 24)])
         p = implementation.FixedPatch2D.from_config(cfg, (4, 4, 12, 12))
         from CelerisLab.lbm import hierarchy
         bad = patch.object(implementation.FixedPatch2D, "from_config",
@@ -169,10 +187,12 @@ class FaultGuardsCPU(unittest.TestCase):
         with patch.object(cuda, "Stream"), patch.object(cuda, "mem_get_info", return_value=(1 << 40, 1 << 40)), bad, patch.object(cuda, "mem_alloc") as alloc, \
                 patch.object(cuda, "memcpy_htod") as upload:
             with self.assertRaisesRegex(RuntimeError, "config chain is inconsistent"):
-                PeriodicPatchStepper(cfg, context=context, refinements=[p.bounds])
+                PeriodicPatchStepper(cfg, context=context, refinements=[p.bounds],
+                                     level_blocks=[(16, 16), (24, 24)])
             with patch.object(hierarchy, "_check_chain"), \
                     self.assertRaisesRegex(RuntimeError, "omega/tau"):
-                PeriodicPatchStepper(cfg, context=context, refinements=[p.bounds])
+                PeriodicPatchStepper(cfg, context=context, refinements=[p.bounds],
+                                     level_blocks=[(16, 16), (24, 24)])
             alloc.assert_not_called()
             upload.assert_not_called()
         module = Mock()
@@ -182,7 +202,8 @@ class FaultGuardsCPU(unittest.TestCase):
                 patch.object(compiler, "load_module", return_value=module), \
                 patch.object(cuda, "memcpy_htod") as upload:
             with self.assertRaisesRegex(RuntimeError, "layout mismatch"):
-                PeriodicPatchStepper(cfg, context=context, refinements=[p.bounds])
+                PeriodicPatchStepper(cfg, context=context, refinements=[p.bounds],
+                                     level_blocks=[(16, 16), (24, 24)])
             upload.assert_not_called()
 
 
@@ -250,23 +271,27 @@ class NormalBuildGPU(unittest.TestCase):
                 parts = line.split()
                 if len(parts) == 3 and parts[0] == "#define":
                     definitions[parts[1]] = parts[2]
-            generated.append((metadata["config"], definitions))
+            generated.append((metadata["config"], definitions, metadata))
             return build
 
         max_error = max_bc_error = 0.0
         with CudaContext(0) as context:
             with patch.object(compiler, "generate_config", side_effect=record_config):
                 backend = implementation.FixedPatchStepper(cfg, context=context, refinements=[bounds],
-                                                            coarse_boundary="slot_reservoir")
+                                                            coarse_boundary="slot_reservoir",
+                                                            level_blocks=whole_tiles(cfg, bounds))
             with backend:
                 self.assertEqual(len(generated), 2)
-                fine_config, fine_defines = generated[1]
-                self.assertEqual((fine_config["nx"], fine_config["ny"]), (fx, fy))
+                fine_config, fine_defines, _ = generated[1]
+                # Blocked kernels see one tile window including the 1-cell halo.
+                self.assertEqual((fine_config["nx"], fine_config["ny"]), (fx+2, fy+2))
                 for name in ("Y_WALL_BC", "TOP_WALL_BC", "WALL_SLOT_WIDTH",
                              "WALL_NORMAL_VELOCITY", "TOP_WALL_NORMAL_VELOCITY"):
                     self.assertEqual(float(fine_defines[name].removesuffix("f")), 0, name)
-                self.assertEqual(int(generated[0][1]["Y_WALL_BC"]), 2)
-                self.assertEqual(int(generated[0][1]["TOP_WALL_BC"]), 1)
+                # The root's physical slot/reservoir walls come from BLOCK_ROOT_BC=1
+                # (module constant d_block_wall); the fine level has none.
+                root_meta, fine_meta = generated[0][2], generated[1][2]  # read before cleanup
+                self.assertEqual((root_meta["block_root_bc"], fine_meta["block_root_bc"]), (1, 0))
                 # Independently derive all target links, including the eight
                 # whose coarse sources are uncovered, from physical midpoints.
                 targets, corners = set(), set()
@@ -286,9 +311,7 @@ class NormalBuildGPU(unittest.TestCase):
                     expected = scalar_interval(*expected, bounds, slot=(amplitude, width))
                     backend.step()
                     self.assertEqual(backend.step_count, step)
-                    raw = [np.empty(a.shape, np.float32) for a in initial]
-                    for level, array in zip(backend._levels, raw):
-                        cuda.memcpy_dtoh(array, level.current)
+                    raw = owner_arrays(backend)
                     # Read raw storage before masking: active NaNs cannot be
                     # hidden by the public download's non-leaf cleanup.
                     self.assertTrue(np.isfinite(raw[0][:, coarse_leaf]).all())
@@ -327,7 +350,7 @@ class NormalBuildGPU(unittest.TestCase):
         max_pop = max_drift = max_moment = 0.0
         with CudaContext(0) as context, PeriodicPatchStepper(
                 LBMConfig(nx=16, ny=16, viscosity=.1), context=context,
-                refinements=[(4, 4, 12, 12)]) as backend:
+                refinements=[(4, 4, 12, 12)], level_blocks=[(16, 16), (24, 24)]) as backend:
             backend.initialize(*postcollision(state, spec, maps))
             before = totals(backend.download(), maps)
             for step in range(1, 9):
@@ -358,17 +381,17 @@ class NormalBuildGPU(unittest.TestCase):
         fine_mask[4:-4, 4:-4] = True
         initial[1][:, ~fine_mask] = np.nan
         expected = scalar_interval(*[a.astype(np.float64) for a in initial], bounds)
-        with CudaContext(0) as context, PeriodicPatchStepper(cfg, context=context, refinements=[bounds]) as backend:
+        with CudaContext(0) as context, PeriodicPatchStepper(cfg, context=context, refinements=[bounds],
+                                                             level_blocks=[(13, 17), (14, 16)]) as backend:
             self.assertEqual(backend.patch.fine_shape, (14, 16))
             p = backend.patch
-            allocated = 2 * 9 * (17*13+16*14) * 4 + sum(a.nbytes for a in
-                        (p.ghost_fine, p.ghost_coarse, p.coalesce_coarse, p.coalesce_fine))
+            allocated = 2 * 4 * sum(int(np.prod(level.layout.blocked_shape)) for level in backend._levels) \
+                + sum(a.nbytes for a in (p.ghost_fine, p.ghost_coarse, p.coalesce_coarse, p.coalesce_fine))
             self.assertLess(allocated, 128*1024)
             backend.initialize(*initial)
             backend.step()
             coarse, fine = backend.download()
-            raw_fine = np.empty((9, 14, 16), np.float32)
-            cuda.memcpy_dtoh(raw_fine, backend._levels[1].current)
+            raw_fine = owner_arrays(backend)[1]
             max_error = 0.0
             for got, want in zip((coarse, raw_fine), expected):
                 np.testing.assert_array_equal(np.isnan(got), np.isnan(want))

tests/validation/multigrid/test_refinement_gpu.py

@@ -52,6 +52,17 @@ def totals(arrays, maps):
     return np.array([populations.sum(), populations @ ref.CX, populations @ ref.CY])
 
 
+def level_storage(level):
+    """Raw device storage of one level as dense (9, ny, nx), dense or blocked."""
+    if level.layout is None:
+        raw = np.empty(level.shape, np.float32)
+        cuda.memcpy_dtoh(raw, level.current)
+        return raw
+    blocked = np.empty(level.layout.blocked_shape, np.float32)
+    cuda.memcpy_dtoh(blocked, level.current)
+    return level.layout.to_dense(blocked)
+
+
 def global_array(module, name, shape):
     out = np.empty(shape, dtype=np.uint32)
     ptr, size = module.get_global(name)
@@ -82,7 +93,7 @@ class CoupledGPU(unittest.TestCase):
         expected = postcollision(next_state, spec, maps)
         arrivals = (next_state.coarse_f, next_state.fine_f)
 
-        cases = (("dense_unshifted", False, {}),
+        cases = (("blocked_unshifted", False, dict(level_blocks=[(16, 16), (24, 24)])),
                  ("blocked_shifted", True,
                   dict(level_blocks=[(4, 4), (6, 6)])))
         with CudaContext(0) as context:
@@ -148,7 +159,8 @@ class CoupledGPU(unittest.TestCase):
         maps = ref.build_maps(spec)
         with CudaContext(0) as context:
             with PeriodicPatchStepper(LBMConfig(nx=16, ny=16, viscosity=.1),
-                                      context=context, refinements=[(4, 4, 12, 12)]) as backend:
+                                      context=context, refinements=[(4, 4, 12, 12)],
+                                      level_blocks=[(16, 16), (24, 24)]) as backend:
                 backend.initialize(*postcollision(ref.make_state(spec), spec, maps))
                 coarse, fine = backend._levels
 
@@ -164,12 +176,94 @@ class CoupledGPU(unittest.TestCase):
                 for operation in (backend.download, backend.step):
                     with self.assertRaisesRegex(RuntimeError, "failed"):
                         operation()
-                raw = np.empty(fine.shape, np.float32)
-                cuda.memcpy_dtoh(raw, fine.current)
+                raw = level_storage(fine)
                 self.assertTrue(np.isfinite(raw.ravel()[backend.patch.coalesce_fine]).all())
             self.assertFalse(context.is_closed)
         print("partial failure: two fine hops completed, no coarse step published; cleanup PASS", flush=True)
 
+    def eight_steps(self, backend, spec, maps, audit):
+        """Rest and seeded nonequilibrium, 8 steps against the scalar reference.
+
+        audit=True also checks the dense CELERIS_REFINEMENT_AUDIT ownership counters.
+        """
+        p = backend.patch
+        for case in ("rest", "nonequilibrium_seed1"):
+            with self.subTest(case=case):
+                state = ref.make_state(spec, velocity=(0, 0) if case == "rest" else (.03, -.02))
+                if case != "rest":
+                    rng = np.random.default_rng(1)
+                    state = ref.State(*(a * (1 + 1e-3 * rng.standard_normal(a.shape))
+                                        for a in (state.coarse_f, state.fine_f)))
+                initial = postcollision(state, spec, maps)
+                backend.initialize(*initial)
+                if audit:  # dense audit build only
+                    for level in backend._levels:
+                        for name in ("collision_counts", "explosion_counts", "incoming_counts"):
+                            ptr, size = level.module.get_global(name)
+                            cuda.memset_d8(ptr, 0, size)
+                before = totals(backend.download(), maps)
+                max_error = max_drift = max_moment_error = max_corner = 0.0
+                for step in range(1, 9):
+                    state, diagnostic = ref.advance_coarse_once(spec, maps, state)
+                    backend.step()
+                    actual = backend.download()
+                    expected = postcollision(state, spec, maps)
+                    self.assertEqual(backend.step_count, step)
+                    p.check_synchronization(step, state.tick, phase="postcollision_g")
+                    for got, want in zip(actual, expected):
+                        np.testing.assert_array_equal(np.isnan(got), np.isnan(want))
+                        error = float(np.nanmax(np.abs(got - want)))
+                        self.assertLessEqual(error, POP_ATOL)
+                        max_error = max(max_error, error)
+                    drift = float(np.max(np.abs(totals(actual, maps) - before)))
+                    moment_error = float(np.max(np.abs(totals(actual, maps) - totals(expected, maps))))
+                    self.assertLessEqual(drift, MOMENT_ATOL)
+                    self.assertLessEqual(moment_error, MOMENT_ATOL)
+                    max_drift, max_moment_error = max(max_drift, drift), max(max_moment_error, moment_error)
+                    c, f = backend._levels
+                    if audit:
+                        np.testing.assert_array_equal(global_array(c.module, "collision_counts", (16, 16)),
+                                                      step * maps.coarse_active)
+                        np.testing.assert_array_equal(global_array(f.module, "collision_counts", (24, 24)),
+                                                      2 * step * maps.fine_active)
+                        np.testing.assert_array_equal(global_array(f.module, "explosion_counts", (24, 24)),
+                                                      step * ~maps.fine_active)
+                        owners = np.zeros((9, 16, 16), np.uint32)
+                        owners.ravel()[p.coalesce_coarse] = step
+                        np.testing.assert_array_equal(global_array(c.module, "incoming_counts", owners.shape), owners)
+                    raw = level_storage(f)
+                    np.testing.assert_array_equal(np.isfinite(raw), diagnostic["hop2_valid"])
+                    self.assertTrue(np.isfinite(raw.ravel()[p.coalesce_fine]).all())
+                    self.assertEqual(len(diagnostic["corner_packets"]), 8)
+                    for packet in diagnostic["corner_packets"].values():
+                        sx, sy = packet["source_coarse_xy"]
+                        self.assertTrue(maps.coarse_active[sy, sx])
+                        values = [raw[packet["direction"], y, x]
+                                  for y, x in packet["children_fine_yx"]]
+                        error = float(np.max(np.abs(np.array(values) - packet["packets"])))
+                        self.assertLessEqual(error, POP_ATOL)
+                        max_corner = max(max_corner, error)
+                print(f"{case}: steps=8 max_population={max_error:.9g} "
+                      f"max_leaf_drift={max_drift:.9g} max_leaf_error={max_moment_error:.9g} "
+                      f"max_corner_packet={max_corner:.9g}; ownership/counts PASS", flush=True)
+
+    def test_bounded_equivalence_blocked(self):
+        """The same eight-step reference and launch-failure checks on whole-level tiles."""
+        spec = ref.make_spec(collision="srt")
+        maps = ref.build_maps(spec)
+        cfg = LBMConfig(nx=16, ny=16, viscosity=.1)
+        with CudaContext(0) as context:
+            with PeriodicPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)],
+                                      level_blocks=[(16, 16), (24, 24)]) as backend:
+                self.eight_steps(backend, spec, maps, audit=False)
+                with patch.object(backend._pairs[0], "explode", side_effect=RuntimeError("injected launch failure")):
+                    with self.assertRaisesRegex(RuntimeError, "injected"):
+                        backend.step()
+                for operation in (backend.download, backend.step):
+                    with self.assertRaisesRegex(RuntimeError, "failed"):
+                        operation()
+            self.assertFalse(context.is_closed)
+
     def test_bounded_equivalence_and_ownership(self):
         self.assertEqual(Path(CelerisLab.__file__).resolve(), ROOT / "src/CelerisLab/__init__.py")
         self.assertEqual(Path(compiler._KERNEL_DIR).resolve(), ROOT / "src/CelerisLab/lbm/kernels")
@@ -198,64 +292,7 @@ class CoupledGPU(unittest.TestCase):
                                   for level in backend._levels)
                 self.assertLess(device_bytes + audit_bytes, 128 * 1024)
                 print(f"device_arrays={device_bytes} audit_arrays={audit_bytes} bytes", flush=True)
-                for case in ("rest", "nonequilibrium_seed1"):
-                    with self.subTest(case=case):
-                        state = ref.make_state(spec, velocity=(0, 0) if case == "rest" else (.03, -.02))
-                        if case != "rest":
-                            rng = np.random.default_rng(1)
-                            state = ref.State(*(a * (1 + 1e-3 * rng.standard_normal(a.shape))
-                                                for a in (state.coarse_f, state.fine_f)))
-                        initial = postcollision(state, spec, maps)
-                        backend.initialize(*initial)
-                        for level in backend._levels:
-                            for name in ("collision_counts", "explosion_counts", "incoming_counts"):
-                                ptr, size = level.module.get_global(name)
-                                cuda.memset_d8(ptr, 0, size)
-                        before = totals(backend.download(), maps)
-                        max_error = max_drift = max_moment_error = max_corner = 0.0
-                        for step in range(1, 9):
-                            state, diagnostic = ref.advance_coarse_once(spec, maps, state)
-                            backend.step()
-                            actual = backend.download()
-                            expected = postcollision(state, spec, maps)
-                            self.assertEqual(backend.step_count, step)
-                            p.check_synchronization(step, state.tick, phase="postcollision_g")
-                            for got, want in zip(actual, expected):
-                                np.testing.assert_array_equal(np.isnan(got), np.isnan(want))
-                                error = float(np.nanmax(np.abs(got - want)))
-                                self.assertLessEqual(error, POP_ATOL)
-                                max_error = max(max_error, error)
-                            drift = float(np.max(np.abs(totals(actual, maps) - before)))
-                            moment_error = float(np.max(np.abs(totals(actual, maps) - totals(expected, maps))))
-                            self.assertLessEqual(drift, MOMENT_ATOL)
-                            self.assertLessEqual(moment_error, MOMENT_ATOL)
-                            max_drift, max_moment_error = max(max_drift, drift), max(max_moment_error, moment_error)
-                            c, f = backend._levels
-                            np.testing.assert_array_equal(global_array(c.module, "collision_counts", (16, 16)),
-                                                          step * maps.coarse_active)
-                            np.testing.assert_array_equal(global_array(f.module, "collision_counts", (24, 24)),
-                                                          2 * step * maps.fine_active)
-                            np.testing.assert_array_equal(global_array(f.module, "explosion_counts", (24, 24)),
-                                                          step * ~maps.fine_active)
-                            owners = np.zeros((9, 16, 16), np.uint32)
-                            owners.ravel()[p.coalesce_coarse] = step
-                            np.testing.assert_array_equal(global_array(c.module, "incoming_counts", owners.shape), owners)
-                            raw = np.empty(f.shape, np.float32)
-                            cuda.memcpy_dtoh(raw, f.current)
-                            np.testing.assert_array_equal(np.isfinite(raw), diagnostic["hop2_valid"])
-                            self.assertTrue(np.isfinite(raw.ravel()[p.coalesce_fine]).all())
-                            self.assertEqual(len(diagnostic["corner_packets"]), 8)
-                            for packet in diagnostic["corner_packets"].values():
-                                sx, sy = packet["source_coarse_xy"]
-                                self.assertTrue(maps.coarse_active[sy, sx])
-                                values = [raw[packet["direction"], y, x]
-                                          for y, x in packet["children_fine_yx"]]
-                                error = float(np.max(np.abs(np.array(values) - packet["packets"])))
-                                self.assertLessEqual(error, POP_ATOL)
-                                max_corner = max(max_corner, error)
-                        print(f"{case}: steps=8 max_population={max_error:.9g} "
-                              f"max_leaf_drift={max_drift:.9g} max_leaf_error={max_moment_error:.9g} "
-                              f"max_corner_packet={max_corner:.9g}; ownership/counts PASS", flush=True)
+                self.eight_steps(backend, spec, maps, audit=True)
 
                 # Fail before any additional time step; no stale synchronized export.
                 with patch.object(backend._pairs[0], "explode", side_effect=RuntimeError("injected launch failure")):

tests/validation/multigrid/test_refinement_jets.py

@@ -3,11 +3,13 @@
 Jets use the same coarse cell-edge frame as fine_objects. The equivalent jet
 of Astra's actuated fixture is SurfaceActuation(0.7, 0.85, U).as_jet((16, 12),
 2.15). Frozen GPU gates (8 coarse steps, both coarse boundary modes):
-  G1 jet run: dense vs 1x1/2x2/4x2 fine blocks bit-identical (values + NaN masks);
-  G2 jet vs legacy actuation (dense): max |dG| over finite fine ACTIVE fluid
+  G1 jet run: 1x1/2x2/4x2 fine blocks (blocked (12, 8) root) bit-identical to the
+     FINE_BLOCKS[0] reference (values + NaN masks); dense leg kept until dense retirement;
+  G2 jet vs legacy actuation (blocked reference): max |dG| over finite fine ACTIVE fluid
      populations <= 1e-6 (velocity rows differ only by FP32 rounding; reported);
-  G3 amplitude change after 4 steps (+U -> -U), 4 more steps: dense vs blocked
-     bit-identical and the run differs from the unchanged run.
+  G3 amplitude change after 4 steps (+U -> -U), 4 more steps: blocked layouts
+     bit-identical to the FINE_BLOCKS[0] reference (dense leg until retirement) and the
+     run differs from the unchanged run.
 """
 from contextlib import nullcontext
 import unittest
@@ -26,6 +28,11 @@ from test_refinement_blocked_fine import (BOUNDS, FINE_BLOCKS, NX, NY, STEPS, ac
 SPEED = .002
 
 
+def root_and(fine_block):
+    """All-blocked levels with a (12, 8) root tile; fine_block None is the dense route."""
+    return [None if fine_block is None else (12, 8), fine_block]
+
+
 def jet(speed=SPEED):
     return actuation(speed).as_jet((16.0, 12.0), 2.15)
 
@@ -42,7 +49,7 @@ def mocked(objects, jets, fine_block):
             patch.object(cuda, "mem_alloc", side_effect=lambda size: Mock()), \
             patch.object(cuda, "memcpy_htod", htod):
         backend = FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS],
-                                    level_blocks=[None, fine_block], level_objects=[(), tuple(objects)],
+                                    level_blocks=root_and(fine_block), level_objects=[(), tuple(objects)],
                                     level_jets=[(), tuple(jets)])
     return backend, htod
 
@@ -161,7 +168,7 @@ def run(mode, fine_block, *, objects, jets, switch=None):
     patch_ = FixedPatch2D.from_config(cfg, BOUNDS)
     with CudaContext(0) as context, FixedPatchStepper(
             cfg, context=context, coarse_boundary=mode, refinements=[BOUNDS],
-            level_blocks=[None, fine_block], level_objects=[(), tuple(objects)],
+            level_blocks=root_and(fine_block), level_objects=[(), tuple(objects)],
             level_jets=[(), tuple(jets)]) as backend:
         backend.initialize(state(NY, NX, 1), state(*patch_.fine_shape, 2))
         if switch is None:
@@ -179,29 +186,32 @@ class JetStepperGPU(unittest.TestCase):
             patch_, boundary = fine_circle(mode)
             fluid = (patch_.fine_kind == ACTIVE) & (
                 (boundary.flags.reshape(patch_.fine_shape) & SOLID) == 0)
-            dense = run(mode, None, objects=[circle()], jets=[jet()])
-            self.assertTrue(np.isfinite(dense[1][1][:, fluid]).all(), mode)
-            for fine_block in FINE_BLOCKS:
+            reference = run(mode, FINE_BLOCKS[0], objects=[circle()], jets=[jet()])
+            self.assertTrue(np.isfinite(reference[1][1][:, fluid]).all(), mode)
+            for fine_block in FINE_BLOCKS[1:]:
                 public, raw_arrays = run(mode, fine_block, objects=[circle()], jets=[jet()])
-                for a, b in zip(public + tuple(raw_arrays), dense[0] + tuple(dense[1])):
+                for a, b in zip(public + tuple(raw_arrays), reference[0] + tuple(reference[1])):
                     np.testing.assert_array_equal(a, b, err_msg=f"G1 {mode} {fine_block}")
-            legacy = run(mode, None, objects=[circle(actuation(SPEED))], jets=())
-            diff = np.abs(dense[1][1][:, fluid] - legacy[1][1][:, fluid]).max()
+            dense = run(mode, None, objects=[circle()], jets=[jet()])  # dense leg (D2 retires)
+            for a, b in zip(dense[0] + tuple(dense[1]), reference[0] + tuple(reference[1])):
+                np.testing.assert_array_equal(a, b, err_msg=f"G1 {mode} dense")
+            legacy = run(mode, FINE_BLOCKS[0], objects=[circle(actuation(SPEED))], jets=())
+            diff = np.abs(reference[1][1][:, fluid] - legacy[1][1][:, fluid]).max()
             print(f"G2 {mode}: jet vs legacy fine max |dG| = {diff:.3e}", flush=True)
             self.assertLessEqual(diff, 1e-6, mode)
 
     def test_g3_amplitude_switch_layouts_identical(self):
         for mode in ("periodic", "slot_reservoir"):
-            dense = run(mode, None, objects=[circle()], jets=[jet()], switch=-SPEED)
-            steady = run(mode, None, objects=[circle()], jets=[jet()])
-            self.assertFalse(np.array_equal(dense[1][1], steady[1][1], equal_nan=True), mode)
-            for fine_block in FINE_BLOCKS:
+            reference = run(mode, FINE_BLOCKS[0], objects=[circle()], jets=[jet()], switch=-SPEED)
+            steady = run(mode, FINE_BLOCKS[0], objects=[circle()], jets=[jet()])
+            self.assertFalse(np.array_equal(reference[1][1], steady[1][1], equal_nan=True), mode)
+            for fine_block in (*FINE_BLOCKS[1:], None):  # None: dense leg (D2 retires)
                 public, raw_arrays = run(mode, fine_block, objects=[circle()], jets=[jet()],
                                          switch=-SPEED)
-                for a, b in zip(public + tuple(raw_arrays), dense[0] + tuple(dense[1])):
+                for a, b in zip(public + tuple(raw_arrays), reference[0] + tuple(reference[1])):
                     np.testing.assert_array_equal(a, b, err_msg=f"G3 {mode} {fine_block}")
-            print(f"G3 {mode}: amplitude switch dense vs {len(FINE_BLOCKS)} layouts identical",
-                  flush=True)
+            print(f"G3 {mode}: amplitude switch, {len(FINE_BLOCKS) - 1} layouts + dense "
+                  f"identical to {FINE_BLOCKS[0]}", flush=True)
 
 
 if __name__ == "__main__":