Support blocked parent storage in fixed patch stepping

src/CelerisLab/lbm/kernels/blocks.cu

@@ -16,6 +16,7 @@
 #include "operators/collision_srt.cuh"
 #include "streaming/pull_double_buffer.cuh"
 #include "boundary/curved_boundary.cuh"
+#include "boundary/planar_velocity.cuh"
 
 #if DIM != 2 || NQ != 9 || COLLISION_MODEL != 0 || STORE_PRECISION != 0 || USE_DDF_SHIFTING || USE_LES || STREAMING_MODEL != 0
 #error "Blocks require unshifted FP32 D2Q9 SRT double buffers without LES"
@@ -24,6 +25,11 @@
 // Fixed-patch cell kinds, matching refinement.py ACTIVE/INNER_GHOST/OUTER_GHOST.
 constexpr unsigned char KIND_ACTIVE = 0u, KIND_INNER = 1u, KIND_OUTER = 2u;
 
+// Zero-initialized for ordinary blocks/fine levels. Set only for parent walls;
+// NX/NY remain local storage strides, never physical wall coordinates.
+struct BlockWall { unsigned int nx, ny, slot_width; float amplitude; };
+extern "C" { __constant__ BlockWall d_block_wall; }
+
 // Halo refresh: block cell ids b*TOTAL_CELLS+c, sources are always core cells.
 extern "C" __global__ void BlockHaloCopy(
     fpxx* f, const uint64_t* dst, const uint64_t* src, unsigned int count)
@@ -38,17 +44,32 @@ extern "C" __global__ void BlockHaloCopy(
 }
 
 // Fixed-patch explosion into blocked fine storage: target = owner block cell
-// id of a patch-ghost cell, donor = dense coarse cell (coarse is Q-major dense).
+// id of a patch-ghost cell; donors are parent block cell ids. Dense parent
+// storage is the one-window case with coarse_cells = full domain size.
 extern "C" __global__ void BlockExplode(
-    const fpxx* coarse, fpxx* fine, const uint64_t* target, const uint32_t* donor,
+    const fpxx* coarse, fpxx* fine, const uint64_t* target, const uint64_t* donor,
     unsigned int count, unsigned long coarse_cells)
 {
     unsigned int row = blockIdx.x * blockDim.x + threadIdx.x;
     if (row >= count) return;
     const uint64_t t = target[row];
     fpxx* to = fine + (t / TOTAL_CELLS) * NQ * TOTAL_CELLS;
+    const uint64_t d = donor[row];
+    const fpxx* from = coarse + (d / coarse_cells) * NQ * coarse_cells;
     for (unsigned int q = 0; q < NQ; ++q)
-        to[index_f(t % TOTAL_CELLS, q)] = coarse[q * coarse_cells + donor[row]];
+        to[index_f(t % TOTAL_CELLS, q)] = from[q * coarse_cells + d % coarse_cells];
+}
+
+// Targets are receiver-local SOURCE slots for the immediately following pull.
+// Keep child summation order identical to PatchCoarseStep; no collision here.
+extern "C" __global__ void BlockCoalesce(
+    const fpxx* fine, fpxx* coarse, const uint64_t* target, const uint64_t* children,
+    unsigned int count)
+{
+    unsigned int row = blockIdx.x * blockDim.x + threadIdx.x;
+    if (row >= count) return;
+    const uint64_t* s = children + 4ul * row;
+    coarse[target[row]] = 0.25f * (fine[s[0]] + fine[s[1]] + fine[s[2]] + fine[s[3]]);
 }
 
 // Static cut links, receiver-local staging: cell is the receiver's local core
@@ -92,6 +113,16 @@ extern "C" __global__ void BlockStep(const fpxx* in, fpxx* out, unsigned int hal
         compute_neighbors(x, y, j);  // Core neighbors stay inside the window: no wrap.
         stream_pull_load(k, f, in + base, j);
         if (cell_kind == KIND_ACTIVE) {
+            if (d_block_wall.nx != 0) {
+                const unsigned int bx = NX - 2 * halo, by = NY - 2 * halo;
+                const unsigned int nbx = d_block_wall.nx / bx;
+                const unsigned int gx = (blockIdx.z % nbx) * bx + x - halo;
+                const unsigned int gy = (blockIdx.z / nbx) * by + y - halo;
+                if (gy == 1)
+                    reconstruct_bottom_velocity(f, 0.0f, bottom_wall_normal_velocity(
+                        gx, d_block_wall.nx, d_block_wall.slot_width, d_block_wall.amplitude));
+                else if (gy == d_block_wall.ny - 2) reconstruct_top_reservoir(f, RHO);
+            }
             float rho, ux, uy, feq[NQ];
             compute_rho_u(f, rho, ux, uy);
             compute_feq(rho, ux, uy, feq);

src/CelerisLab/lbm/kernels/boundary/planar_velocity.cuh

@@ -3,16 +3,26 @@
 // the caller owns collision/store. No force, shifting or off-lattice walls.
 #ifndef CELERIS_BOUNDARY_PLANAR_VELOCITY_CUH
 #define CELERIS_BOUNDARY_PLANAR_VELOCITY_CUH
+#include <math.h>
+
+__device__ __forceinline__ float bottom_wall_normal_velocity(
+    unsigned int x, unsigned int width, unsigned int slot_width, float amplitude)
+{
+    // Width is the physical domain width, even when NX is a block window.
+    if (slot_width > 0 && slot_width < width) {
+        const float s = (2.0f*(float)x-(float)(width-1))/(float)slot_width;
+        const float w = fmaxf(0.0f, 1.0f-s*s);
+        return amplitude*w*w;
+    }
+    return amplitude; // Disabled or explicitly uniform full width.
+}
 
 __device__ __forceinline__ float bottom_wall_normal_velocity(unsigned int x)
 {
 #if WALL_SLOT_WIDTH > 0 && WALL_SLOT_WIDTH < NX
-    // C1 compact profile centered at (NX-1)/2; amplitude is its continuous peak.
-    const float s = (2.0f*(float)x-(float)(NX-1))/(float)WALL_SLOT_WIDTH;
-    const float w = fmaxf(0.0f, 1.0f-s*s);
-    return WALL_NORMAL_VELOCITY*w*w;
+    return bottom_wall_normal_velocity(x, NX, WALL_SLOT_WIDTH, WALL_NORMAL_VELOCITY);
 #else
-    return WALL_NORMAL_VELOCITY; // Disabled or explicitly uniform full width.
+    return WALL_NORMAL_VELOCITY;
 #endif
 }
 

src/CelerisLab/lbm/refinement_stepper.py

@@ -18,7 +18,7 @@ from .blocks import BlockLayout2D, blocked_link_tables
 from .curved_links import surface_link_velocity
 from .block_stepper import MAX_GRID_DIM
 from .refinement import ACTIVE, HALO, OUTER_GHOST, FixedPatch2D, prepare_fine_boundary
-from .descriptors import SOLID
+from .descriptors import SOLID, D2Q9_EX, D2Q9_EY
 
 
 def blocked_fine_tables(patch, layout, order):
@@ -37,6 +37,27 @@ def blocked_fine_tables(patch, layout, order):
             layout.slot(owner[cell.astype(np.intp)], q), kind)
 
 
+def blocked_base_tables(patch, layout, order, wall_aware=False):
+    """Parent donor cells and receiver-local arrival slots, in transfer row order.
+
+    Stage arrivals only after child reads and parent halo gather, immediately
+    before parent pull. Source *owner* slots are wrong at block seams.
+    """
+    owner = layout.owner.ravel()
+    q, receiver = np.divmod(patch.coalesce_coarse[order], np.uint64(owner.size))
+    offset = (np.asarray(D2Q9_EY)[q.astype(np.intp)] * layout.local_shape[1]
+              + np.asarray(D2Q9_EX)[q.astype(np.intp)])
+    source = owner[receiver.astype(np.intp)].astype(np.int64) - offset
+    kind = np.full(layout.nblocks * int(np.prod(layout.local_shape)), OUTER_GHOST, np.uint8)
+    leaves = np.ones(patch.coarse_shape, bool)
+    x0, y0, x1, y1 = patch.bounds
+    leaves[y0:y1, x0:x1] = False
+    if wall_aware:
+        leaves[[0, -1]] = False
+    kind[owner[leaves.ravel()].astype(np.intp)] = ACTIVE
+    return (owner[patch.ghost_coarse].astype(np.uint64), layout.slot(source, q), kind)
+
+
 @dataclass
 class _Level:
     shape: tuple
@@ -45,6 +66,11 @@ class _Level:
     spare: object
     kernel: object
     grid: tuple
+    layout: object = None
+    kind: object = None
+    halo_copy: object = None
+    halo_dst: object = None
+    halo_src: object = None
 
 
 class FixedPatchStepper:
@@ -75,11 +101,14 @@ class FixedPatchStepper:
     Blocked fine_objects are stationary circles, optionally surface-actuated:
     prepare_fine_boundary decides legality on the dense fine grid, then cut
     links are staged receiver-local per hop: halo -> BFL -> step, same buffer.
+    ``base_block=(by, bx)`` also blocks the parent; it requires ``fine_block``.
+    Parent halos precede receiver-local coalescence staging and pull/collision;
+    covered cells are explicitly poisoned on every step. Public shapes stay dense.
     Blocked mode applies prescribed surface velocity only; its per-link
     exchange diagnostic is unavailable (rows keep physical cut-link order).
     """
     def __init__(self, cfg, bounds, *, context, coarse_boundary="periodic", fine_objects=(),
-                 fine_block=None):
+                 fine_block=None, base_block=None):
         if HALO != 4:
             raise RuntimeError("refinement.cu requires a four-cell fine halo")
         if coarse_boundary not in ("periodic", "slot_reservoir"):
@@ -94,6 +123,17 @@ class FixedPatchStepper:
             if self.patch.bounds[1] < 4 or self.patch.bounds[3] > cfg.ny - 4:
                 raise ValueError("slot_reservoir requires y0 >= 4 and y1 <= ny-4")
         fine_objects = tuple(fine_objects)
+        if base_block is not None and fine_block is None:
+            raise ValueError("base_block requires fine_block")
+        self.base_layout = None
+        if base_block is not None:
+            if len(base_block) != 2:
+                raise ValueError("base_block must be (by, bx)")
+            self.base_layout = BlockLayout2D.build(
+                cfg.ny, cfg.nx, *base_block, periodic=(True, not self._wall_aware))
+            if (self.base_layout.nblocks > MAX_GRID_DIM
+                    or self.base_layout.local_shape[0] > MAX_GRID_DIM):
+                raise ValueError("base blocks exceed launch limits")
         self.fine_layout = None
         if fine_block is not None:
             if len(fine_block) != 2:
@@ -121,7 +161,7 @@ class FixedPatchStepper:
         self._levels = []
         self._stream = None
         self._fine_flags = self._fine_curved = self._curved = None
-        self._kind = self._halo_copy = self._halo_dst = self._halo_src = None
+        self._coalesce = None
         self._link_tables = {}
         self._block = (cfg.threads_per_block, 1, 1)
         fine_cfg = replace(cfg, nx=self.patch.fine_shape[1], ny=self.patch.fine_shape[0],
@@ -141,14 +181,15 @@ class FixedPatchStepper:
         with self.ctx.activate():
             try:
                 self._stream = cuda.Stream()
-                layout = self.fine_layout
-                fine = ((fine_cfg, "PatchFineStep", "refinement.cu") if layout is None else
-                        (replace(fine_cfg, nx=layout.local_shape[1], ny=layout.local_shape[0],
-                                 y_wall_bc="bounce_back", top_wall_bc="velocity", wall_slot_width=0,
-                                 wall_normal_velocity=0.0, velocity=0.0), "BlockStep", "blocks.cu"))
-                for config, name, tau, source in (
-                        (replace(cfg), "PatchCoarseStep", self.patch.tau_coarse, "refinement.cu"),
-                        (*fine[:2], self.patch.tau_fine, fine[2])):
+                for index, (config, layout, tau, name) in enumerate((
+                        (replace(cfg), self.base_layout, self.patch.tau_coarse, "PatchCoarseStep"),
+                        (fine_cfg, self.fine_layout, self.patch.tau_fine, "PatchFineStep"))):
+                    source = "refinement.cu"
+                    if layout is not None:
+                        config = replace(config, nx=layout.local_shape[1], ny=layout.local_shape[0],
+                                         y_wall_bc="bounce_back", top_wall_bc="velocity",
+                                         wall_slot_width=0, wall_normal_velocity=0.0, velocity=0.0)
+                        name, source = "BlockStep", "blocks.cu"
                     config.validate()
                     if not np.isclose(config.omega, 1 / tau, rtol=1e-14, atol=0):
                         raise RuntimeError("level omega/tau mismatch")
@@ -165,33 +206,45 @@ class FixedPatchStepper:
                     if size != struct.calcsize("=ffffI"):
                         raise RuntimeError("LBMParams host/device layout mismatch")
                     cuda.memcpy_htod(ptr, struct.pack("=ffffI", config.omega, 0, 0, 0, n_objects))
-                    shape = (9, config.ny, config.nx)
-                    size = int(np.prod(shape)) * 4
-                    grid = ((config.nx + self._block[0] - 1) // self._block[0], config.ny, 1)
-                    if name == "BlockStep":  # public dense shape, private blocked allocation
-                        shape = (9, *self.patch.fine_shape)
-                        size = int(np.prod(layout.blocked_shape)) * 4
-                        grid = grid[:2] + (layout.nblocks,)
-                    self._levels.append(_Level(shape, module, self._allocate(size),
-                                              self._allocate(size), module.get_function(name), grid))
+                    shape = (9, *(self.patch.coarse_shape if index == 0 else self.patch.fine_shape))
+                    size = int(np.prod(shape if layout is None else layout.blocked_shape)) * 4
+                    grid = ((config.nx + self._block[0] - 1) // self._block[0], config.ny,
+                            1 if layout is None else layout.nblocks)
+                    level = _Level(shape, module, self._allocate(size), self._allocate(size),
+                                   module.get_function(name), grid, layout)
+                    self._levels.append(level)
+                    if layout is not None:
+                        level.halo_copy = module.get_function("BlockHaloCopy")
+                        if len(layout.halo_dst):
+                            level.halo_dst = self._upload(layout.halo_dst)
+                            level.halo_src = self._upload(layout.halo_src)
+                        if index == 0 and self._wall_aware:
+                            ptr, size = module.get_global("d_block_wall")
+                            if size != struct.calcsize("=IIIf"):
+                                raise RuntimeError("Block wall host/device layout mismatch")
+                            cuda.memcpy_htod(ptr, struct.pack("=IIIf", cfg.nx, cfg.ny,
+                                                            cfg.wall_slot_width, cfg.wall_normal_velocity))
                 p = self.patch
+                c, f = self._levels
                 order = np.argsort(p.coalesce_coarse)
-                tables = None if layout is None else blocked_fine_tables(p, layout, order)
-                self._ghosts = self._upload(p.ghost_fine if layout is None else tables[0])
-                self._donors = self._upload(p.ghost_coarse)
-                self._targets = self._upload(p.coalesce_coarse[order])
-                if layout is None:
-                    self._children = self._upload(p.coalesce_fine[order])
-                    self._explode = self._levels[1].module.get_function("PatchExplode")
-                else:
-                    self._children = self._upload(tables[1])
-                    self._kind = self._upload(tables[2])
-                    if len(layout.halo_dst):  # a 1x1 non-periodic layout has no pairs
-                        self._halo_dst = self._upload(layout.halo_dst)
-                        self._halo_src = self._upload(layout.halo_src)
-                    module = self._levels[1].module
-                    self._explode = module.get_function("BlockExplode")
-                    self._halo_copy = module.get_function("BlockHaloCopy")
+                fine_tables = (None if f.layout is None else
+                               blocked_fine_tables(p, f.layout, order))
+                base_tables = (None if c.layout is None else
+                               blocked_base_tables(p, c.layout, order, self._wall_aware))
+                self._ghosts = self._upload(p.ghost_fine if fine_tables is None else fine_tables[0])
+                donors = p.ghost_coarse if base_tables is None else base_tables[0]
+                self._donors = self._upload(donors if fine_tables is None else donors.astype(np.uint64))
+                self._targets = self._upload(p.coalesce_coarse[order] if base_tables is None
+                                             else base_tables[1])
+                self._children = self._upload(p.coalesce_fine[order] if fine_tables is None
+                                              else fine_tables[1])
+                self._explode = f.module.get_function("PatchExplode" if fine_tables is None
+                                                       else "BlockExplode")
+                if fine_tables is not None:
+                    f.kind = self._upload(fine_tables[2])
+                if base_tables is not None:
+                    c.kind = self._upload(base_tables[2])
+                    self._coalesce = c.module.get_function("BlockCoalesce")
                 if self._fine_tables is not None:
                     self._fine_flags = self._upload(self._fine_tables["flags"])
                     self._link_tables = {k: self._upload(v) for k, v in self._fine_tables.items()
@@ -257,8 +310,8 @@ class FixedPatchStepper:
                 if not np.isfinite(array[:, mask]).all() or (array[:, mask].sum(axis=0) <= 0).any():
                     raise ValueError("Leaf populations must be finite with positive density")
             self._mask_nonleaves(*arrays)
-            if self.fine_layout is not None:
-                arrays[1] = self.fine_layout.to_blocked(arrays[1])
+            arrays = [a if level.layout is None else level.layout.to_blocked(a)
+                      for a, level in zip(arrays, self._levels)]
             try:
                 for level, array in zip(self._levels, arrays):
                     cuda.memcpy_htod(level.current, array)
@@ -303,10 +356,10 @@ class FixedPatchStepper:
                      np.uint32(count), velocity if velocity is not None else np.uintp(0),
                      grid=((count + self._block[0] - 1) // self._block[0], 1), **kw)
 
-    def _fine_halo(self, buffer, **kw):
-        pairs = len(self.fine_layout.halo_dst)
+    def _gather_halo(self, level, buffer, **kw):
+        pairs = len(level.layout.halo_dst)
         if pairs:
-            self._halo_copy(buffer, self._halo_dst, self._halo_src, np.uint32(pairs),
+            level.halo_copy(buffer, level.halo_dst, level.halo_src, np.uint32(pairs),
                             grid=((pairs + self._block[0] - 1) // self._block[0], 1), **kw)
 
     def step(self, n=1):
@@ -322,7 +375,9 @@ class FixedPatchStepper:
                 self._status = "advancing"
                 for _ in range(n):
                     self._explode(c.current, f.current, self._ghosts, self._donors,
-                                  np.uint32(len(p.ghost_fine)), np.uint64(np.prod(p.coarse_shape)),
+                                  np.uint32(len(p.ghost_fine)),
+                                  np.uint64(np.prod(p.coarse_shape if c.layout is None
+                                                    else c.layout.local_shape)),
                                   grid=((len(p.ghost_fine) + self._block[0] - 1) // self._block[0], 1),
                                   **kw)
                     if self.fine_layout is not None:
@@ -330,10 +385,10 @@ class FixedPatchStepper:
                         # Same buffer per hop: halo -> receiver-local BFL -> step.
                         flags = self._fine_flags if self._fine_flags is not None else np.uintp(0)
                         for src, dst in ((f.current, f.spare), (f.spare, f.current)):
-                            self._fine_halo(src, **kw)
+                            self._gather_halo(f, src, **kw)
                             self._blocked_fine_boundary(src, **kw)
                             f.kernel(src, dst, np.uint32(self.fine_layout.halo), flags,
-                                     self._kind, grid=f.grid, **kw)
+                                     f.kind, grid=f.grid, **kw)
                     else:
                         flags = self._fine_flags if self._fine_flags is not None else np.uintp(0)
                         if self.fine_boundary is not None and self._fine_flags is None:
@@ -342,10 +397,21 @@ class FixedPatchStepper:
                         f.kernel(f.current, f.spare, flags, grid=f.grid, **kw)
                         self._apply_fine_boundary(f.spare, **kw)
                         f.kernel(f.spare, f.current, flags, grid=f.grid, **kw)
-                    c.kernel(c.current, c.spare, f.current, self._targets, self._children,
-                             np.uint32(len(p.coalesce_coarse)),
-                             *(np.uint32(v) for v in p.bounds), np.uint32(self._wall_aware),
-                             grid=c.grid, **kw)
+                    if c.layout is None:
+                        c.kernel(c.current, c.spare, f.current, self._targets, self._children,
+                                 np.uint32(len(p.coalesce_coarse)),
+                                 *(np.uint32(v) for v in p.bounds), np.uint32(self._wall_aware),
+                                 grid=c.grid, **kw)
+                    else:
+                        # Child donors must be read before staging mutates parent G.
+                        # No halo gather may overwrite these receiver-local arrivals.
+                        self._gather_halo(c, c.current, **kw)
+                        self._coalesce(f.current, c.current, self._targets, self._children,
+                                       np.uint32(len(p.coalesce_coarse)),
+                                       grid=((len(p.coalesce_coarse) + self._block[0] - 1)
+                                             // self._block[0], 1), **kw)
+                        c.kernel(c.current, c.spare, np.uint32(c.layout.halo), np.uintp(0),
+                                 c.kind, grid=c.grid, **kw)
                     self._stream.synchronize()
                     c.current, c.spare = c.spare, c.current
                     self._step_count += 1
@@ -358,17 +424,16 @@ class FixedPatchStepper:
         """Return two synchronized G arrays, NaN on all non-leaves; no host cache."""
         with self._lock, self.ctx.activate():
             self._require("ready")
-            arrays = [np.empty(level.shape, np.float32) for level in self._levels]
-            if self.fine_layout is not None:
-                arrays[1] = np.empty(self.fine_layout.blocked_shape, np.float32)
+            arrays = [np.empty(level.shape if level.layout is None else level.layout.blocked_shape,
+                               np.float32) for level in self._levels]
             try:
                 for level, array in zip(self._levels, arrays):
                     cuda.memcpy_dtoh(array, level.current)
             except BaseException:
                 self._status = "failed"
                 raise
-            if self.fine_layout is not None:
-                arrays[1] = self.fine_layout.to_dense(arrays[1])
+            arrays = [a if level.layout is None else level.layout.to_dense(a)
+                      for a, level in zip(arrays, self._levels)]
             self._mask_nonleaves(*arrays)
             return tuple(arrays)
 
@@ -387,7 +452,7 @@ class FixedPatchStepper:
                 except Exception as exc:
                     error = exc
                 self._levels.clear()
-                self._explode = self._halo_copy = None
+                self._explode = self._coalesce = None
                 if getattr(self, "_fine_curved", None) is not None:
                     try:
                         self._fine_curved.free()

tests/validation/multigrid/test_refinement_blocked_base.py

@@ -0,0 +1,324 @@
+"""All-blocked fixed patch: transfer ordering, failure ownership and CUDA parity.
+
+Run CPU and GPU classes in separate processes. GPU callers must first check
+resources/device availability; this file never selects or reserves a GPU.
+"""
+from contextlib import ExitStack, contextmanager, nullcontext
+from dataclasses import replace
+from pathlib import Path
+import struct
+import unittest
+from unittest.mock import Mock, patch
+
+import numpy as np
+import pycuda.driver as cuda
+
+from CelerisLab.cuda import compiler_v2 as compiler
+from CelerisLab.lbm import refinement_stepper
+from CelerisLab.lbm.blocks import BlockLayout2D
+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
+
+BASE_BLOCKS = ((24, 32), (12, 8), (3, 4), (1, 1))
+FINE_BLOCK = (16, 20)
+
+
+def raw(backend):
+    arrays = []
+    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)
+        arrays.append(array if level.layout is None else level.layout.to_dense(array))
+    return tuple(arrays)
+
+
+def exact(test, actual, expected, label):
+    """Compare the complete invalid mask and every finite FP32 bit (incl. -0)."""
+    np.testing.assert_array_equal(np.isnan(actual), np.isnan(expected), err_msg=label)
+    test.assertFalse(np.isinf(actual).any(), label)
+    test.assertFalse(np.isinf(expected).any(), label)
+    finite = np.isfinite(expected)
+    np.testing.assert_array_equal(actual[finite].view(np.uint32),
+                                  expected[finite].view(np.uint32), err_msg=label)
+
+
+@contextmanager
+def mocked_backend(mode="periodic", fail_alloc=None):
+    """Use the real constructor, uploads and launch logic, recording driver calls."""
+    calls, allocations, uploads, functions = [], [], [], {}
+    stream = Mock()
+    context = Mock(activate=nullcontext, sm_arch="sm_70")
+
+    def function(name):
+        if name not in functions:
+            functions[name] = Mock(side_effect=lambda *a, **kw: calls.append((name, a, kw)))
+        return functions[name]
+
+    def allocate(size):
+        if len(allocations) == fail_alloc:
+            raise MemoryError("injected allocation failure")
+        if size <= 0:
+            raise AssertionError("empty allocation")
+        pointer = Mock()
+        allocations.append(pointer)
+        return pointer
+
+    module = Mock(get_function=Mock(side_effect=function))
+    module.get_global.side_effect = lambda name: (name, 16 if name == "d_block_wall" else 20)
+    with ExitStack() as stack:
+        for name in ("generate_config", "compile_kernel", "cleanup_build_dir"):
+            stack.enter_context(patch.object(compiler, name))
+        stack.enter_context(patch.object(compiler, "load_module", return_value=module))
+        stack.enter_context(patch.object(cuda, "Stream", return_value=stream))
+        stack.enter_context(patch.object(cuda, "mem_alloc", side_effect=allocate))
+        stack.enter_context(patch.object(cuda, "memcpy_htod",
+                                        side_effect=lambda p, a: uploads.append((p, a))))
+        yield dict(context=context, calls=calls, allocations=allocations, uploads=uploads,
+                   functions=functions, stream=stream)
+
+
+class BlockedBaseCPU(unittest.TestCase):
+    def test_reject_before_cuda(self):
+        context = Mock(activate=Mock(side_effect=AssertionError("CUDA reached")))
+        for kwargs in ({"base_block": (12, 8)},
+                       {"base_block": (7, 8), "fine_block": FINE_BLOCK},
+                       {"base_block": (12,), "fine_block": FINE_BLOCK}):
+            with self.subTest(kwargs=kwargs), self.assertRaises(ValueError):
+                FixedPatchStepper(config("periodic"), BOUNDS, context=context, **kwargs)
+        context.activate.assert_not_called()
+
+    def test_all_arrival_rows_donors_and_kind_masks(self):
+        p = FixedPatch2D.from_config(config("periodic"), BOUNDS)
+        order = np.argsort(p.coalesce_coarse)
+        self.assertEqual(len(order), 172)
+        for mode in ("periodic", "slot_reservoir"):
+            leaves = np.ones((NY, NX), bool)
+            x0, y0, x1, y1 = BOUNDS
+            leaves[y0:y1, x0:x1] = False
+            if mode != "periodic":
+                leaves[[0, -1]] = False
+            for by, bx in BASE_BLOCKS:
+                layout = BlockLayout2D.build(NY, NX, by, bx, periodic=(True, mode == "periodic"))
+                donors, targets, kind = blocked_base_tables(p, layout, order, mode != "periodic")
+                ly, lx = layout.local_shape
+                L = ly * lx
+                expected, source_cells = [], []
+                for slot in p.coalesce_coarse[order]:
+                    q, receiver = divmod(int(slot), NX * NY)
+                    y, x = divmod(receiver, NX)
+                    block = (y // by) * (NX // bx) + x // bx
+                    # Scalar coordinate construction independent of layout.owner / slot.
+                    sy, sx = y % by + 1 - D2Q9_EY[q], x % bx + 1 - D2Q9_EX[q]
+                    expected.append((block * 9 + q) * L + sy * lx + sx)
+                    source_cells.append((y - D2Q9_EY[q], x - D2Q9_EX[q]))
+                np.testing.assert_array_equal(targets, expected)
+                self.assertEqual(len(set(targets.tolist())), 172)
+                self.assertEqual(sum(not (y0 <= y < y1 and x0 <= x < x1)
+                                     for y, x in source_cells), 8)
+                expected_donors = []
+                for cell in p.ghost_coarse:
+                    y, x = divmod(int(cell), NX)
+                    block = (y // by) * (NX // bx) + x // bx
+                    expected_donors.append(block * L + (y % by + 1) * lx + x % bx + 1)
+                np.testing.assert_array_equal(donors, expected_donors)
+                expected_kind = np.full(layout.nblocks * L, OUTER_GHOST, np.uint8)
+                expected_kind[layout.owner[leaves].astype(np.intp)] = ACTIVE
+                np.testing.assert_array_equal(kind, expected_kind)
+
+    def test_staging_negative_controls_detect_seam_and_donor_corruption(self):
+        p = FixedPatch2D.from_config(config("periodic"), BOUNDS)
+        order = np.argsort(p.coalesce_coarse)
+        layout = BlockLayout2D.build(NY, NX, 3, 4)
+        donors, targets, _ = blocked_base_tables(p, layout, order)
+        dense = np.arange(9 * NY * NX, dtype=np.float32).reshape(9, NY, NX)
+        initial = layout.to_blocked(dense)
+        arrivals = np.arange(len(targets), dtype=np.float32) + 10000
+        staged = initial.copy()
+        staged.ravel()[targets.astype(np.intp)] = arrivals
+        def pull(storage):
+            out = np.empty_like(dense)
+            by, bx = layout.block
+            ly, lx = layout.local_shape
+            for y in range(NY):
+                for x in range(NX):
+                    b = (y // by) * (NX // bx) + x // bx
+                    for q in range(9):
+                        out[q, y, x] = storage[b, q, y % by + 1 - D2Q9_EY[q],
+                                               x % bx + 1 - D2Q9_EX[q]]
+            return out
+
+        expected = np.array([np.roll(dense[q], (D2Q9_EY[q], D2Q9_EX[q]), axis=(0, 1))
+                             for q in range(9)])
+        expected.ravel()[p.coalesce_coarse[order].astype(np.intp)] = arrivals
+        leaves = np.ones((NY, NX), bool)
+        leaves[6:18, 8:24] = False
+        # Compare ALL parent arrivals, not only the slots deliberately written.
+        np.testing.assert_array_equal(pull(staged)[:, leaves], expected[:, leaves])
+        # A late halo gather erases receiver-local arrivals at seams.
+        late_gather = staged.copy()
+        layout.exchange(late_gather)
+        self.assertGreater(np.count_nonzero(pull(late_gather)[:, leaves] != expected[:, leaves]), 0)
+        # Canonical source owners are not the slots read by each receiver's pull.
+        q, receiver = np.divmod(p.coalesce_coarse[order], NX * NY)
+        y, x = np.divmod(receiver.astype(np.intp), NX)
+        owners = layout.owner[y - np.asarray(D2Q9_EY)[q.astype(np.intp)],
+                              x - np.asarray(D2Q9_EX)[q.astype(np.intp)]]
+        wrong = initial.copy()
+        wrong.ravel()[layout.slot(owners, q).astype(np.intp)] = arrivals
+        self.assertGreater(np.count_nonzero(pull(wrong)[:, leaves] != expected[:, leaves]), 0)
+        # With one parent block, all eight grazing sources are owner cores.
+        # Each is read for four fine ghost children: late explosion corrupts32.
+        single = BlockLayout2D.build(NY, NX, NY, NX)
+        donors, targets, _ = blocked_base_tables(p, single, order)
+        initial = single.to_blocked(dense)
+        staged = initial.copy()
+        staged.ravel()[targets.astype(np.intp)] = arrivals
+        donor_slots = single.slot(donors[None, :], np.arange(9)[:, None]).astype(np.intp)
+        self.assertEqual(np.count_nonzero(staged.ravel()[donor_slots]
+                                          != initial.ravel()[donor_slots]), 32)
+
+    def test_raw_mask_gate_detects_leak_hidden_by_download_mask(self):
+        valid = state(NY, NX, 1)
+        valid[:, 6:18, 8:24] = np.nan
+        leaked = valid.copy()
+        leaked[5, 6, 8] = 42  # skipped invalid-cell store leaves stale staging data
+        public = leaked.copy()
+        public[:, 6:18, 8:24] = np.nan
+        exact(self, public, valid, "masked download cannot detect stale staging")
+        with self.assertRaises(AssertionError):
+            exact(self, leaked, valid, "raw owner-core must detect stale staging")
+
+    def test_launch_order_stream_wall_constants_and_buffer_swaps(self):
+        with mocked_backend("slot_reservoir") as m:
+            backend = FixedPatchStepper(config("slot_reservoir"), BOUNDS, context=m["context"],
+                                        coarse_boundary="slot_reservoir", fine_block=FINE_BLOCK,
+                                        base_block=(3, 4), fine_objects=[circle(actuation(.002))])
+            try:
+                backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2))
+                coarse, fine = backend._levels
+                current, spare, fine_current = coarse.current, coarse.spare, fine.current
+                backend.step(2)
+                names = [name for name, _, _ in m["calls"]]
+                interval = ["BlockExplode", "BlockHaloCopy", "BlockCurvedBoundary", "BlockStep",
+                            "BlockHaloCopy", "BlockCurvedBoundary", "BlockStep",
+                            "BlockHaloCopy", "BlockCoalesce", "BlockStep"]
+                self.assertEqual(names, interval * 2)
+                self.assertEqual(backend.step_count, 2)
+                self.assertIs(coarse.current, current)
+                self.assertIs(coarse.spare, spare)
+                self.assertIs(fine.current, fine_current)
+                for _, _, kw in m["calls"]:
+                    self.assertIs(kw["stream"], m["stream"])
+                for offset, src, dst in ((0, current, spare), (10, spare, current)):
+                    explode, gather, stage, collide = [m["calls"][offset + i][1] for i in (0, 7, 8, 9)]
+                    self.assertIs(explode[0], src)
+                    self.assertIs(gather[0], src)
+                    self.assertIs(stage[0], fine_current)
+                    self.assertIs(stage[1], src)
+                    self.assertIs(collide[0], src)
+                    self.assertIs(collide[1], dst)
+                    self.assertIs(collide[4], coarse.kind)
+                walls = [a for ptr, a in m["uploads"] if isinstance(ptr, str) and ptr == "d_block_wall"]
+                self.assertEqual(len(walls), 1)
+                width, height, slot, speed = struct.unpack("=IIIf", walls[0])
+                self.assertEqual((width, height, slot), (NX, NY, 8))
+                self.assertAlmostEqual(speed, .001)
+            finally:
+                backend.close()
+            for allocation in m["allocations"]:
+                allocation.free.assert_called_once_with()
+
+    def test_allocation_failure_releases_every_owned_buffer(self):
+        # Fail during the first level, second level and transfer-table uploads.
+        for fail_at in (1, 6, 10):
+            with self.subTest(fail_at=fail_at), mocked_backend(fail_alloc=fail_at) as m:
+                with self.assertRaisesRegex(MemoryError, "injected allocation"):
+                    FixedPatchStepper(config("periodic"), BOUNDS, context=m["context"],
+                                      fine_block=FINE_BLOCK, base_block=(3, 4))
+                self.assertEqual(len(m["allocations"]), fail_at)
+                for allocation in m["allocations"]:
+                    allocation.free.assert_called_once_with()
+
+    def test_failed_parent_stage_cannot_publish_or_resume(self):
+        with mocked_backend() as m:
+            backend = FixedPatchStepper(config("periodic"), BOUNDS, context=m["context"],
+                                        fine_block=FINE_BLOCK, base_block=(3, 4))
+            backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2))
+            m["functions"]["BlockCoalesce"].side_effect = RuntimeError("injected parent stage failure")
+            current = backend._levels[0].current
+            with self.assertRaisesRegex(RuntimeError, "injected parent stage"):
+                backend.step()
+            self.assertEqual(backend.step_count, 0)
+            self.assertEqual(backend._status, "failed")
+            self.assertIs(backend._levels[0].current, current)
+            with self.assertRaisesRegex(RuntimeError, "failed"):
+                backend.download()
+            with self.assertRaisesRegex(RuntimeError, "failed"):
+                backend.step()
+            backend.close()
+            for allocation in m["allocations"]:
+                allocation.free.assert_called_once_with()
+
+
+class BlockedBaseGPU(unittest.TestCase):
+    def test_bodyfree_exact_parity_and_raw_invalid_slots(self):
+        self.compare([(None, ((24, 32), (3, 4), (1, 1)), 8)])
+
+    def test_actuated_circle_exact_parity_and_raw_invalid_slots(self):
+        self.compare([(.002, ((12, 8),), 8), (-.002, ((12, 8),), 8)])
+
+    def test_non_power_two_slot_width_exact_parity(self):
+        self.compare([(None, ((3, 4),), 7), (None, ((3, 4),), 10)], ("slot_reservoir",))
+
+    def compare(self, cases, modes=("periodic", "slot_reservoir")):
+        from CelerisLab.cuda.context import CudaContext
+        expected_source = Path(__file__).resolve().parents[3] / "src/CelerisLab/lbm/refinement_stepper.py"
+        self.assertEqual(Path(refinement_stepper.__file__).resolve(), expected_source)
+        with CudaContext(0) as context:
+            for mode in modes:
+                for speed, layouts, width in cases:
+                    cfg = config(mode)
+                    if mode == "slot_reservoir":
+                        cfg = replace(cfg, wall_slot_width=width)
+                    p = FixedPatch2D.from_config(cfg, BOUNDS)
+                    arrays = state(NY, NX, 1), state(*p.fine_shape, 2)
+                    reference = None
+                    for base_block in (None, *layouts):
+                        label = f"{mode} width={width} speed={speed} base={base_block}"
+                        snapshots = []
+                        with FixedPatchStepper(cfg, BOUNDS, context=context, coarse_boundary=mode,
+                                               fine_block=FINE_BLOCK, base_block=base_block,
+                                               fine_objects=() if speed is None else [circle(actuation(speed))]) as backend:
+                            backend.initialize(*arrays)
+                            for step in (0, 1, 2, 8):
+                                if step == 1 and base_block is not None:
+                                    # Blocked staging dirties invalid input slots; its output
+                                    # must poison them. Dense parent does not stage there.
+                                    coarse = backend._levels[0]
+                                    spare = coarse.layout.to_blocked(np.full(coarse.shape, 42, np.float32))
+                                    cuda.memcpy_htod(coarse.spare, spare)
+                                backend.step(step - backend.step_count)
+                                public, owner = backend.download(), raw(backend)
+                                self.assertEqual(backend.step_count, step)
+                                self.assertTrue(np.isnan(owner[0][:, 6:18, 8:24]).all(), label)
+                                if mode != "periodic":
+                                    self.assertTrue(np.isnan(owner[0][:, (0, -1), :]).all(), label)
+                                snapshots.append(public + owner)
+                            backend.initialize(*arrays)
+                            self.assertEqual(backend.step_count, 0)
+                            for actual, expected in zip(backend.download() + raw(backend), snapshots[0]):
+                                exact(self, actual, expected, label + " reset")
+                        if reference is None:
+                            reference = snapshots
+                        else:
+                            for step, actual, expected in zip((0, 1, 2, 8), snapshots, reference):
+                                for index, (a, b) in enumerate(zip(actual, expected)):
+                                    exact(self, a, b, f"{label} step={step} array={index}")
+                            print(f"PASS {label}: full finite bits/NaN masks at 0,1,2,8 + reset", flush=True)
+
+
+if __name__ == "__main__":
+    unittest.main()

tests/validation/multigrid/test_refinement_blocked_fine.py

@@ -128,8 +128,8 @@ class BlockedFineCPU(unittest.TestCase):
             backend = FixedPatchStepper(config("periodic"), BOUNDS, context=context,
                                         fine_block=FINE_BLOCKS[0])
             self.assertEqual(len(backend.fine_layout.halo_dst), 0)
-            self.assertIsNone(backend._halo_dst)
-            backend._fine_halo(Mock())  # no launch without pairs
+            self.assertIsNone(backend._levels[1].halo_dst)
+            backend._gather_halo(backend._levels[1], Mock())  # no launch without pairs
             module.get_function.return_value.assert_not_called()
             backend.close()
         self.assertTrue(sizes)
@@ -187,7 +187,7 @@ class BlockedFineCPU(unittest.TestCase):
             backend, calls = self.mocked_backend([circle(spec)])
             c, f = backend._levels
             coarse, cur, spare = c.current, f.current, f.spare
-            flags, kind = backend._fine_flags, backend._kind
+            flags, kind = backend._fine_flags, backend._levels[1].kind
             velocity = backend._link_tables.get("surface_velocity")
             self.assertEqual(velocity is None, spec is None)
             backend.step(1)

tests/validation/multigrid/test_refinement_followup.py

@@ -111,7 +111,7 @@ class FaultGuardsCPU(unittest.TestCase):
         backend._lock = RLock()
         backend._status = "ready"
         backend._stream = Mock()
-        backend._levels = [Mock(shape=(9, 1, 1)), Mock(shape=(9, 1, 1))]
+        backend._levels = [Mock(shape=(9, 1, 1), layout=None), Mock(shape=(9, 1, 1), layout=None)]
         backend._allocations = [Mock(), Mock()]
         return backend