test_refinement_mrt: check force reads in executable text, not comments

The block_collision.cuh header comment names the absent d_params.fx/fy
read, which tripped the source-text assertion. Strip // comments first
and also require the real d_params.omega read. Test-only; no kernel source
change.

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

docs/development/src_coherence.md

@@ -8,8 +8,16 @@ pending a GPU run).
 
 - Level `l = 0 .. depth-1` (0 = root). Pair `hierarchy.pairs[l-1]` builds level `l`.
 - Array shapes are `(ny, nx)`; rectangles, origins and centers are `(x, y)`.
-- FixedPatchStepper bodies/Jets use the root EDGE frame (root cell `i` = `[i, i+1)`);
-  uniform `Simulation` objects use the integer-NODE frame. Nothing converts between them.
+- Coarse/fine name the two members (parent/child) of ANY pair, not absolute L0/L1:
+  `fine_shape`, `fine_origin`, `prepare_fine_boundary`, `blocked_fine_tables`,
+  `blocked_coarse_tables`.
+- Refined levels (l >= 1) of FixedPatchStepper place bodies/Jets in the root EDGE frame
+  (root cell `i` = `[i, i+1)`). Uniform `Simulation` objects AND the static circles of a
+  single-level (depth 1) FixedPatchStepper use the integer-NODE frame. Nothing converts
+  between them. This concerns input placement only: within the EDGE frame, refined
+  levels still project root-edge coordinates to fine-local cell indices with the
+  half-cell shift (`prepare_fine_boundary`: `(v - o)/dx - .5`; hits back via `fine_hits_to_edge_frame`
+  `origin + dx*(hit + .5)`), using `edge_frame`/`fine_origin` below.
 - `Hierarchy.edge_frame(l)` -> `((ox, oy), h)`: cell `(i, j)` of level `l`'s allocated
   array (ghost band and holes included) covers `[o + h*i, o + h*(i+1))`. A scalar
   FixedPatch2D window carries the half-coarse-cell ghost shift `bounds - HALO/2` in
@@ -26,8 +34,11 @@ pending a GPU run).
 | `base_layout`, `fine_layout`, `child_layout` (stepper properties) | `layouts[0]`, `layouts[1]`, `layouts[2]` (`FixedPatchSet2D.fine_layout(block)` is unchanged) |
 | `fine_boundary`, `child_boundary` | `boundaries[1]`, `boundaries[2]` (host geometry; `None` = no bodies/Jets) |
 | `_jets` | `jets[1]` |
-| `_boundaries[l]` (dict from L1, level keys) | `_level_boundaries[l]` (list from L0, `None` at L0 unless single-level static circles) |
-| `_fine_tables`, `_actuated`, `_jet_cache` | `_level_boundaries[1].tables` / `.actuated` / `.jet_cache` |
+| `_boundaries[l]` (dict from L1, level keys) | `_device_boundaries[l]` (list from L0, `None` at L0 unless single-level static circles) |
+| `_level_boundaries` (private, b4eacf7) | `_device_boundaries` |
+| `_fine_tables`, `_actuated`, `_jet_cache` | `_device_boundaries[1].tables` / `.actuated` / `.jet_cache` |
+| `_remap_cfg` (private) | `_root_cfg` |
+| `refinement_tables.blocked_base_tables` (also via `refinement_stepper`) | `refinement_tables.blocked_coarse_tables` (import from `refinement_tables`; the stepper no longer re-exports) |
 | `_fine_leaf`, `_child_leaf` | `leaf[1]`, `leaf[2]` (copies) |
 | `_middle_kind` | `hierarchy.kinds[1]` |
 | `_boundary_hosts` (private) | `_boundaries` (private list behind `boundaries`) |
@@ -56,6 +67,10 @@ the replaced level relative to its parent).
 
 ## Not done / open
 
+- Two per-level config derivations stay distinct on purpose: `Hierarchy.configs`
+  (`_middle_config`, wall fields reset on every refined level, geometry validation) vs.
+  `FixedPatchStepper._kernel_configs` (reset only on level 1 of a slot_reservoir root;
+  generates each level's config.h). Folding them needs a config.h/PTX diff gate.
 - No single-`Simulation` facade over FixedPatchStepper: the two own different state models
   (uniform node frame, body sync, forcing, checkpoint vs. immutable blocked hierarchy, edge
   frame, remap-as-new-instance). A facade would only forward calls; no concrete need found.

src/CelerisLab/lbm/blocks.py

@@ -18,7 +18,9 @@ all NQ populations from an owner core cell into another block's halo cell.
 Periodicity is explicit per axis. Across a non-periodic exterior, halo cells
 have no pair and stay invalid (NaN after packing): no wrap, no implied BC.
 
-Preparation only: no CUDA, no boundary conditions, no refinement levels.
+Host preparation plus thin launch wrappers: no CUDA import or allocation
+(launch_* only call kernels and device tables the caller owns), no refinement
+levels; level_config and static_circle_tables prepare one blocked level.
 The one-cell same-level halo here is unrelated to the four-cell fine halo of
 the fixed refinement patch (refinement.py), whose semantics are unchanged.
 

src/CelerisLab/lbm/hierarchy.py

@@ -26,8 +26,15 @@ class Hierarchy:
     Index conventions (the single place they are defined):
       * level l = 0 (root) .. depth-1; pair index p = l-1 builds level l from level p;
       * array shapes are (ny, nx); rectangles, origins and centers are (x, y);
-      * the root EDGE frame puts root cell i at [i, i+1), centers at i+.5 (bodies and
-        Jets of FixedPatchStepper); uniform Simulation objects use the integer-NODE frame;
+      * a pair's two members are its COARSE (parent) and FINE (child) level, for every
+        pair, not absolute L0/L1 (fine_shape, fine_origin, blocked_coarse_tables, ...);
+      * refined levels (l >= 1) of FixedPatchStepper place bodies and Jets in the root
+        EDGE frame: root cell i covers [i, i+1), center i+.5. Uniform Simulation objects
+        and the static circles of a single-level (depth 1) FixedPatchStepper use the
+        integer-NODE frame (CircleGeometry.build_flag_mask at integer (x, y)). No
+        NODE<->EDGE input conversion is applied anywhere; inside the EDGE frame, refined
+        boundaries still project root-edge coordinates to fine-local cell indices with
+        the half-cell shift ((v - o)/h - .5, see prepare_fine_boundary);
       * cell (i, j) of level l's allocated array (ghost band and holes included) covers
         [o + h*i, o + h*(i+1)) with (o, h) = edge_frame(l); level 1 is built from the root
         config directly, deeper levels from their parent frame (see edge_frame).
@@ -115,7 +122,7 @@ def leaf_partition(h):
     for level, kind in enumerate(h.kinds):
         if level:
             pair = h.pairs[level - 1]
-            py, px = (pair.fine_level_offset if isinstance(pair, FixedPatch2D) else (0, 0))
+            py, px = pair.fine_level_offset
             oy, ox = 2 * oy + py, 2 * ox + px
         y, x = np.nonzero(kind == ACTIVE)
         np.add.at(base, ((y + oy) >> level, (x + ox) >> level), 0.25 ** level)

src/CelerisLab/lbm/kernels/operators/block_collision.cuh

@@ -1,10 +1,12 @@
+// CelerisLab – operators/block_collision.cuh
 // The one collision extension point of the blocked steps (BlockStep, BlockChannelStep).
-// Compile-time selection: MRT shares one path; SRT keeps two explicit arithmetic
-// variants, unchanged from their original kernels:
-//   level  (BLOCK_ROOT_BC 0/1): unclamped omega, no forcing term; shifted storage
-//          builds feq from the retained drho and direct momentum (numerical variant);
-//   channel (BLOCK_ROOT_BC 2): matches uniform OneStep SRT arithmetic (omega clamp,
-//          zero Guo forcing term, compute_feq).
+// Unforced by construction: FixedPatchStepper supplies no forcing, so no runtime
+// d_params.fx/fy branch is read. Two compile-time arithmetic families:
+//   uniform-equivalent (MRT on every level; SRT on the channel root, BLOCK_ROOT_BC 2):
+//          collide_dispatch arithmetic with zero force and no LES: clamped omega;
+//          SRT keeps the zero Guo forcing term of collide_srt;
+//   level SRT (BLOCK_ROOT_BC 0/1): unclamped omega, no forcing term; shifted storage
+//          builds feq from the retained drho and direct momentum (numerical variant).
 #ifndef CELERIS_OPERATORS_BLOCK_COLLISION_CUH
 #define CELERIS_OPERATORS_BLOCK_COLLISION_CUH
 
@@ -16,17 +18,19 @@
 #endif
 
 __device__ __forceinline__ void block_collide(float* f) {
-#if COLLISION_MODEL == 2
-    collide_mrt_unforced(f);
-#elif BLOCK_ROOT_BC == 2
-    float rho, ux, uy, feq[NQ], forcing[NQ] = {};
+    float rho, ux, uy;
     compute_rho_u(f, rho, ux, uy);
-    compute_feq(rho, ux, uy, feq);
+#if COLLISION_MODEL == 2 || BLOCK_ROOT_BC == 2
     const float omega = fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega));
+  #if COLLISION_MODEL == 2
+    collide_mrt_no_force(f, rho, ux, uy, omega);
+  #else
+    float feq[NQ], forcing[NQ] = {};
+    compute_feq(rho, ux, uy, feq);
     collide_srt(f, feq, forcing, omega);
+  #endif
 #else
-    float rho, ux, uy, feq[NQ];
-    compute_rho_u(f, rho, ux, uy);
+    float feq[NQ];
 #if USE_DDF_SHIFTING
     // Preserve the density perturbation before adding the unit reference.
     const float drho = f[0]+f[1]+f[2]+f[3]+f[4]+f[5]+f[6]+f[7]+f[8];

src/CelerisLab/lbm/kernels/operators/collision_mrt.cuh

@@ -114,16 +114,6 @@ __device__ __forceinline__ void collide_mrt_no_force(float* __restrict__ g,
     collide_mrt(g, rho, ux, uy, Fin, omega);
 }
 
-// Unforced MRT with the clamped device omega, as in collide_dispatch with zero
-// forcing. Shared by the blocked, channel-root and dense-patch collisions.
-__device__ __forceinline__ void collide_mrt_unforced(float* __restrict__ g)
-{
-    float rho, ux, uy;
-    compute_rho_u(g, rho, ux, uy);
-    collide_mrt_no_force(g, rho, ux, uy,
-        fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega)));
-}
-
 #elif LATTICE_MODEL == LATTICE_D3Q19
 
 __device__ __forceinline__ void collide_mrt(float* __restrict__ g,

src/CelerisLab/lbm/refinement.py

@@ -120,6 +120,11 @@ class FixedPatch2D:
         (an index offset, not an (x, y) origin)."""
         return 2 * self.bounds[1] - HALO, 2 * self.bounds[0] - HALO
 
+    @property
+    def rectangles(self) -> tuple:
+        """The one parent rectangle, as a sequence like FixedPatchSet2D.rectangles."""
+        return (self.bounds,)
+
     @property
     def covered(self) -> np.ndarray:
         """Dense coarse mask of cells replaced by the patch (not leaves)."""
@@ -225,11 +230,18 @@ class FixedPatchSet2D:
     def bounds(self) -> tuple:
         return tuple(p.bounds for p in self.patches)
 
+    rectangles = bounds
+
     @property
     def fine_origin(self) -> tuple[float, float]:
         """Lower edge of the level-wide fine grid in coarse physical units."""
         return 0.0, 0.0
 
+    @property
+    def fine_level_offset(self) -> tuple[int, int]:
+        """The level-wide grid itself: zero (y, x) index offset."""
+        return 0, 0
+
     @property
     def covered(self) -> np.ndarray:
         mask = np.zeros(self.coarse_shape, bool)

src/CelerisLab/lbm/refinement_stepper.py

@@ -22,7 +22,7 @@ from .channel import channel_masks, validate_channel_config, validate_channel_tr
 from ..body.jets import jet_link_velocity
 from .hierarchy import build_hierarchy
 from .refinement import ACTIVE, HALO, FixedPatch2D, FixedPatchSet2D, prepare_fine_boundary
-from .refinement_tables import (blocked_base_tables, blocked_fine_tables,  # noqa: F401 (re-export)
+from .refinement_tables import (blocked_coarse_tables, blocked_fine_tables,
                                 blocked_kind, projected_explosion_tables, _prepare_blocked_boundary,
                                 _BlockedBoundary,
                                 BLOCKED_BOUNDARY_DEVICE_COLUMNS)
@@ -66,8 +66,12 @@ class _Level:
 def _kernel_configs(cfg, pairs, wall_aware):
     """Per-level kernel LBMConfig (window shape, viscosity doubling per level).
 
-    Geometry configs live in Hierarchy.configs; level 1 of a slot_reservoir root
-    replaces the root wall fields by the internal bounce-back/velocity ghost setup.
+    Deliberately distinct from Hierarchy.configs (refinement._middle_config), which
+    resets the wall fields on EVERY refined level for geometry validation. Here only
+    level 1 of a slot_reservoir root resets them (internal bounce-back/velocity ghost
+    setup); other levels inherit their parent's fields. These configs generate each
+    level's compiled config.h, so folding the two derivations is a numerical-parity
+    change (config.h/PTX diff gate), not a style one.
     """
     configs = [replace(cfg)]
     for index, pair in enumerate(pairs):
@@ -107,6 +111,57 @@ def _check_levels(refinements, blocks, objects, jets):
             raise ValueError("bodies/Jets are supported only on the finest level of a nested hierarchy")
 
 
+# Host preparation (numpy only, before any device work); each returns one entry per level.
+
+def _level_layouts(cfg, pairs, level_blocks, wall_aware):
+    """Blocked storage layout of every level, checked against the launch limits."""
+    layouts = []
+    for level, block in enumerate(level_blocks):
+        if level == 0:
+            layout = BlockLayout2D.build(cfg.ny, cfg.nx, *block, periodic=(True, not wall_aware))
+        elif isinstance(pairs[level - 1], FixedPatchSet2D):
+            layout = pairs[level - 1].fine_layout(block)
+        else:
+            layout = BlockLayout2D.build(*pairs[level - 1].fine_shape, *block,
+                                         periodic=(False, False))
+        if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM:
+            raise ValueError(f"level {level} blocks exceed launch limits")
+        layouts.append(layout)
+    return layouts
+
+
+def _host_boundaries(hierarchy, level_cfgs, level_objects, level_jets, wall_aware, shifted):
+    """Host boundary geometry of every refined level in the root edge frame (None at L0)."""
+    hosts = [None]
+    for level in range(1, hierarchy.depth):
+        pair = hierarchy.pairs[level - 1]
+        if level == 1:  # the root config carries the wall-aware boundary fields
+            host = prepare_fine_boundary(level_cfgs[0], pair, level_objects[1],
+                                         wall_aware=wall_aware, jets=level_jets[1],
+                                         _shifted_blocked=shifted)
+        else:
+            origin, spacing = hierarchy.edge_frame(level - 1)
+            host = prepare_fine_boundary(level_cfgs[level - 1], pair, level_objects[level],
+                                         jets=level_jets[level], parent_origin=origin,
+                                         parent_spacing=spacing, _shifted_blocked=shifted)
+        hosts.append(host)
+    return hosts
+
+
+def _fluid_leaves(hierarchy, boundaries):
+    """ACTIVE cells of every level that no static circle or body makes solid."""
+    leaves = [kind == ACTIVE for kind in hierarchy.kinds]
+    for level, b in enumerate(boundaries):
+        if b is None:
+            continue
+        if b.host is None:  # single level: static root circles
+            leaves[level] &= ~b.tables["solid"]
+        else:
+            fine_shape = hierarchy.pairs[level - 1].fine_shape
+            leaves[level] &= (b.host.flags.reshape(fine_shape) & SOLID) == 0
+    return leaves
+
+
 class FixedPatchStepper:
     """D2Q9 SRT/MRT FP32 fixed 2:1 refinement hierarchy, built from level lists.
 
@@ -119,7 +174,8 @@ class FixedPatchStepper:
     supported).
     ``refinements=[]`` is one blocked level: periodic SRT root without LES,
     Jets or exchange observation; its ``level_objects[0]`` are stationary
-    passive circles (static BFL).
+    passive circles (static BFL) in the integer-NODE frame, like Simulation
+    objects (see Hierarchy for frames).
 
     ``coarse_boundary`` is ``periodic`` (default) or ``slot_reservoir``:
     a bottom on-site slot and top fixed-density reservoir from ``cfg``.
@@ -131,7 +187,7 @@ class FixedPatchStepper:
     Periodic mode does not apply the wall/slot boundary fields in ``cfg``.
     Outer walls remain coarse; fine halo donors must be ordinary fluid.
 
-    Bodies are fixed circles in root (L0) physical coordinates; velocities are
+    Bodies of refined levels are fixed circles in root (L0) physical coordinates; velocities are
     in root lattice units (dx=dt=1). Independent local copies are projected at
     fine cell centers; their complete stencil must be active and disjoint from
     the actual two-hop transfer footprint. These use the edge/cell-center frame
@@ -241,40 +297,15 @@ class FixedPatchStepper:
         self.hierarchy = build_hierarchy(cfg, refinements, wall_aware=self._wall_aware)
         pairs = self.hierarchy.pairs
         if self._wall_aware:
-            rects = (pairs[0].bounds,) if isinstance(pairs[0], FixedPatch2D) else pairs[0].bounds
-            if any(r[1] < 4 or r[3] > cfg.ny - 4 for r in rects):
+            if any(r[1] < 4 or r[3] > cfg.ny - 4 for r in pairs[0].rectangles):
                 raise ValueError("slot_reservoir requires y0 >= 4 and y1 <= ny-4")
-        self._layouts = []
-        for level, block in enumerate(level_blocks):
-            if level == 0:
-                layout = BlockLayout2D.build(cfg.ny, cfg.nx, *block,
-                                             periodic=(True, not self._wall_aware))
-            elif isinstance(pairs[level - 1], FixedPatchSet2D):
-                layout = pairs[level - 1].fine_layout(block)
-            else:
-                layout = BlockLayout2D.build(*pairs[level - 1].fine_shape, *block,
-                                             periodic=(False, False))
-            if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM:
-                raise ValueError(f"level {level} blocks exceed launch limits")
-            self._layouts.append(layout)
+        self._layouts = _level_layouts(cfg, pairs, level_blocks, self._wall_aware)
         if self._channel:
-            base = blocked_base_tables(pairs[0], self._layouts[0], np.argsort(pairs[0].coalesce_coarse))
+            base = blocked_coarse_tables(pairs[0], self._layouts[0], np.argsort(pairs[0].coalesce_coarse))
             validate_channel_transfers(pairs[0], self._layouts[0], base[1])
         self._level_cfgs = _kernel_configs(cfg, pairs, self._wall_aware)
-        # Host boundary geometry per level in the root edge frame (Hierarchy.edge_frame).
-        self._boundaries = [None]
-        for level in range(1, depth):
-            pair = pairs[level - 1]
-            if level == 1:  # the root config carries the wall-aware boundary fields
-                host = prepare_fine_boundary(cfg, pair, level_objects[1],
-                                             wall_aware=self._wall_aware, jets=level_jets[1],
-                                             _shifted_blocked=self._ddf_shifting)
-            else:
-                origin, spacing = self.hierarchy.edge_frame(level - 1)
-                host = prepare_fine_boundary(self._level_cfgs[level - 1], pair, level_objects[level],
-                                             jets=level_jets[level], parent_origin=origin,
-                                             parent_spacing=spacing, _shifted_blocked=self._ddf_shifting)
-            self._boundaries.append(host)
+        hosts = _host_boundaries(self.hierarchy, self._level_cfgs, level_objects, level_jets,
+                                 self._wall_aware, self._ddf_shifting)
         if not isinstance(observe_raw_exchange, bool):
             raise ValueError("observe_raw_exchange must be bool")
         self._observe_raw_exchange = observe_raw_exchange
@@ -289,28 +320,22 @@ class FixedPatchStepper:
         self._exchange_spacing = tuple(self.hierarchy.edge_frame(l)[1] for l in range(depth))
         self._exchange_centers = tuple(tuple(tuple(o.geometry.center) for o in objects)
                                        for objects in level_objects)
-        self._level_jets = level_jets
         root = None
         if depth == 1 and level_objects[0]:
             tables = static_circle_tables(self._layouts[0], level_objects[0])
             root = _BlockedBoundary(None, tables, (), False)
-        self._level_boundaries = [root] + [
-            _prepare_blocked_boundary(self._boundaries[level], self._layouts[level],
+        # The one boundary/Jet record per level; host geometry and Jets are views of it.
+        self._device_boundaries = [root] + [
+            _prepare_blocked_boundary(hosts[level], self._layouts[level],
                                       self._level_cfgs[level], level_jets[level],
                                       _shifted_blocked=self._ddf_shifting)
             for level in range(1, depth)]
-        self._leaf = [kind == ACTIVE for kind in self.hierarchy.kinds]
-        if root is not None:
-            self._leaf[0] &= ~root.tables["solid"]
-        for level in range(1, depth):
-            host = self._boundaries[level]
-            if host is not None:
-                self._leaf[level] &= (host.flags.reshape(pairs[level - 1].fine_shape) & SOLID) == 0
+        self._leaf = _fluid_leaves(self.hierarchy, self._device_boundaries)
         self._state = [mask.copy() for mask in self._leaf]
         if self._channel:
             self._leaf[0] &= channel_masks((cfg.ny, cfg.nx))[0]
         # Private constructor inputs for remap_finest.
-        self._remap_cfg = replace(cfg)
+        self._root_cfg = replace(cfg)
         self._level_blocks = level_blocks
         self.ctx = context
         self._lock = RLock()
@@ -366,7 +391,7 @@ class FixedPatchStepper:
         for index, layout in enumerate(self._layouts):
             total += 2 * 4 * int(np.prod(layout.blocked_shape))
             total += layout.halo_dst.nbytes + layout.halo_src.nbytes
-            b = self._level_boundaries[index]
+            b = self._device_boundaries[index]
             if b is not None:
                 total += sum(np.asarray(v).nbytes for k, v in b.tables.items()
                              if k == "flags" or k in BLOCKED_BOUNDARY_DEVICE_COLUMNS)
@@ -389,19 +414,19 @@ class FixedPatchStepper:
         grid = ((config.nx + self._block[0] - 1) // self._block[0], config.ny, layout.nblocks)
         level = _Level(shape, module, self._allocate(size), self._allocate(size),
                        module.get_function(name), grid, layout)
-        level.boundary = self._level_boundaries[index]
+        level.boundary = self._device_boundaries[index]
         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._channel:
-            cfg = self._remap_cfg
+            cfg = self._root_cfg
             ptr, size = module.get_global("d_block_channel")
             if size != struct.calcsize("=IIf"):
                 raise RuntimeError("Block channel host/device layout mismatch")
             cuda.memcpy_htod(ptr, struct.pack("=IIf", cfg.nx, cfg.ny, cfg.velocity))
         if index == 0 and self._wall_aware:
-            cfg = self._remap_cfg
+            cfg = self._root_cfg
             ptr, size = module.get_global("d_block_wall")
             if size != struct.calcsize("=IIIf"):
                 raise RuntimeError("Block wall host/device layout mismatch")
@@ -415,11 +440,11 @@ class FixedPatchStepper:
         order = np.argsort(pair.coalesce_coarse)
         ghosts, children, child_kind = blocked_fine_tables(pair, child.layout, order)
         if index == 0:
-            donors, targets, parent_kind = blocked_base_tables(pair, parent.layout, order,
+            donors, targets, parent_kind = blocked_coarse_tables(pair, parent.layout, order,
                                                                self._wall_aware)
         else:
             # Nested pairs keep inherited INNER/OUTER and sparse holes of the parent.
-            donors, targets, parent_kind = blocked_base_tables(
+            donors, targets, parent_kind = blocked_coarse_tables(
                 pair, parent.layout, order, composed_kind=self.hierarchy.kinds[index])
         result = _Pair(child.module.get_function("BlockExplode"),
                        parent.module.get_function("BlockCoalesce"),
@@ -551,19 +576,17 @@ class FixedPatchStepper:
             pairs = self.hierarchy.pairs
             if (self._wall_aware or self._channel or self.hierarchy.depth != 3
                     or not isinstance(pairs[1], FixedPatchSet2D)
-                    or self._boundaries[1] is not None
-                    or any(b is not None for b in self._level_boundaries[:2])):
+                    or any(b is not None for b in self._device_boundaries[:2])):
                 raise ValueError("remap_finest requires a periodic three-level hierarchy with a "
                                  "level-2 set and no level-0/1 bodies")
-            cfg = self._remap_cfg
+            cfg = self._root_cfg
             objects = tuple(objects)
-            actual = self._level_boundaries[2]
+            actual = self._device_boundaries[2]
             jets = () if actual is None else actual.jets
             old = RemapTopology.from_config(cfg, pairs[0], pairs[1].bounds,
                                            child_objects=objects, child_jets=jets,
                                            _shifted_blocked=self._ddf_shifting)
-            if (not same_physical_boundary(old.boundary, self._boundaries[2])
-                    or (self._boundaries[2] is not None and actual is None)):
+            if not same_physical_boundary(old.boundary, None if actual is None else actual.host):
                 raise ValueError("objects must reproduce the original root-edge-frame level-2 boundary")
             if actual is not None:
                 velocity = actual.tables.get("surface_velocity")
@@ -617,7 +640,7 @@ class FixedPatchStepper:
             return candidate
 
     def set_jet_amplitude(self, index, amplitude, *, level):
-        """Replace Jet ``index`` of ``level_jets[level]`` between steps; no geometry/table rebuild.
+        """Replace Jet ``index`` of level ``level`` (``jets[level]``) between steps; no table rebuild.
 
         Invalid input raises before any change. A device/sync error marks the
         instance failed (as step/initialize do); host state commits only after
@@ -628,27 +651,24 @@ class FixedPatchStepper:
             if (isinstance(level, bool) or not isinstance(level, Integral)
                     or not 1 <= level < self.hierarchy.depth):
                 raise ValueError("Jet level must name an existing refined level (1 .. depth-1)")
-            boundary = self._level_boundaries[level]
-            jets = list(self._level_jets[level])
+            boundary = self._device_boundaries[level]
+            # Host preparation: validate and build the new column without touching state.
+            jets = list(() if boundary is None else boundary.jets)
             jets[index] = jets[index].with_amplitude(amplitude)  # validates first
-            velocity = None
-            if boundary is not None and boundary.jet_cache is not None:
-                shape, jet_id = boundary.jet_cache
-                velocity = boundary.tables["surface_velocity"].copy()
-                driven = jet_id >= 0
-                velocity[driven] = jet_link_velocity(jets, shape, jet_id)[driven]
+            shape, jet_id = boundary.jet_cache
+            velocity = boundary.tables["surface_velocity"].copy()
+            driven = jet_id >= 0
+            velocity[driven] = jet_link_velocity(jets, shape, jet_id)[driven]
+            # Device action: an uncertain write fails the instance closed.
             try:
-                if velocity is not None:
-                    self._stream.synchronize()  # no in-flight step reads the column
-                    cuda.memcpy_htod(boundary.device["surface_velocity"], velocity)
+                self._stream.synchronize()  # no in-flight step reads the column
+                cuda.memcpy_htod(boundary.device["surface_velocity"], velocity)
             except BaseException:
                 self._status = "failed"
                 raise
-            if velocity is not None:
-                boundary.tables["surface_velocity"] = velocity
-            if boundary is not None:
-                boundary.jets = tuple(jets)
-            self._level_jets[level] = tuple(jets)
+            # Commit: the one boundary record holds both the column and the Jets.
+            boundary.tables["surface_velocity"] = velocity
+            boundary.jets = tuple(jets)
 
     def _blocked_boundary(self, level, buffer, **kw):
         b = level.boundary
@@ -735,7 +755,7 @@ class FixedPatchStepper:
                 self._exchange_valid = False
                 for _ in range(n):
                     if self._observe_raw_exchange:  # stream-ordered: holds the last interval only
-                        for b in self._level_boundaries:
+                        for b in self._device_boundaries:
                             if b is not None:
                                 cuda.memset_d8_async(b.interval_integral, 0,
                                                      len(b.tables["cell"]) * 32, self._stream)
@@ -758,7 +778,7 @@ class FixedPatchStepper:
                 raise RuntimeError("window exchange observation is disabled")
             self._window_valid = False
             try:
-                for b in self._level_boundaries:
+                for b in self._device_boundaries:
                     if b is not None and b.window_integral is not None:
                         cuda.memset_d8(b.window_integral, 0, len(b.tables["cell"]) * 32)
             except BaseException:
@@ -771,7 +791,7 @@ class FixedPatchStepper:
         if (isinstance(level, bool) or not isinstance(level, Integral)
                 or not 0 <= level < len(self._levels)):
             raise ValueError("level must name an existing level")
-        b = self._level_boundaries[level]
+        b = self._device_boundaries[level]
         if (b is None or isinstance(body_id, bool) or not isinstance(body_id, Integral)
                 or not 0 <= body_id < len(b.host.objects)):
             raise ValueError("body_id must name a body on this level")
@@ -868,7 +888,7 @@ class FixedPatchStepper:
                         self._stream.synchronize()
                 except Exception as exc:
                     error = exc
-                for b in getattr(self, "_level_boundaries", ()):
+                for b in getattr(self, "_device_boundaries", ()):
                     if b is not None:
                         b.flags = b.device = b.kernel = None
                         b.last_raw = b.interval_integral = b.accumulate = None
@@ -896,9 +916,20 @@ class FixedPatchStepper:
                 raise
 
     # Read-only per-level state, indexed by level (0 = root). Pairs: hierarchy.pairs[l-1].
-    layouts = property(lambda self: tuple(self._layouts))
-    boundaries = property(lambda self: tuple(self._boundaries),
-                          doc="Host boundary geometry per level (None: no bodies/Jets).")
-    jets = property(lambda self: tuple(self._level_jets))
-    leaf = property(lambda self: tuple(m.copy() for m in self._leaf),
-                    doc="Per-level (ny, nx) masks of public (non-NaN) download cells.")
+    @property
+    def layouts(self):
+        return tuple(self._layouts)
+
+    @property
+    def boundaries(self):
+        """Host boundary geometry per level (None: no bodies/Jets on a refined level)."""
+        return tuple(None if b is None else b.host for b in self._device_boundaries)
+
+    @property
+    def jets(self):
+        return tuple(() if b is None else b.jets for b in self._device_boundaries)
+
+    @property
+    def leaf(self):
+        """Per-level (ny, nx) masks of public (non-NaN) download cells."""
+        return tuple(m.copy() for m in self._leaf)

src/CelerisLab/lbm/refinement_tables.py

@@ -36,7 +36,7 @@ def blocked_kind(layout, dense_kind):
     kind[owner[owned]] = np.asarray(dense_kind).ravel()[owned]
     return kind
 
-def blocked_base_tables(patch, layout, order, wall_aware=False, *, composed_kind=None):
+def blocked_coarse_tables(patch, layout, order, wall_aware=False, *, composed_kind=None):
     """Parent donor cells and receiver-local arrival slots, in transfer row order.
 
     Stage arrivals only after child reads and parent halo gather, immediately

tests/unit/test_channel_root.py

@@ -7,7 +7,7 @@ from CelerisLab.config import LBMConfig
 from CelerisLab.lbm.blocks import BlockLayout2D
 from CelerisLab.lbm.channel import channel_masks, validate_channel_config, validate_channel_transfers
 from CelerisLab.lbm.refinement import FixedPatch2D
-from CelerisLab.lbm.refinement_tables import blocked_base_tables
+from CelerisLab.lbm.refinement_tables import blocked_coarse_tables
 from .test_numerical_operators import compile_run
 
 
@@ -31,7 +31,7 @@ def transfer(bounds, block=(4, 6)):
     cfg = LBMConfig(nx=24, ny=20, viscosity=.1)
     pair = FixedPatch2D.from_config(cfg, bounds)
     layout = BlockLayout2D.build(20, 24, *block)
-    targets = blocked_base_tables(pair, layout, np.argsort(pair.coalesce_coarse))[1]
+    targets = blocked_coarse_tables(pair, layout, np.argsort(pair.coalesce_coarse))[1]
     return pair, layout, targets
 
 

tests/unit/test_jet_update_failure.py

@@ -102,7 +102,7 @@ def test_simulation_before_initialize_commits_host_only(monkeypatch):
 
 
 def stepper_state(run):
-    boundary = run._level_boundaries[2]
+    boundary = run._device_boundaries[2]
     return boundary.tables["surface_velocity"].copy(), boundary.jets, run.jets[2]
 
 
@@ -116,8 +116,8 @@ def test_stepper_invalid_input_writes_nothing_and_stays_ready():
             with pytest.raises(error):
                 run.set_jet_amplitude(*args, **kw)
         module.cuda.memcpy_htod.assert_not_called()
-        np.testing.assert_array_equal(run._level_boundaries[2].tables["surface_velocity"], velocity)
-        assert (run._level_boundaries[2].jets, run.jets[2]) == (jets, level_jets)
+        np.testing.assert_array_equal(run._device_boundaries[2].tables["surface_velocity"], velocity)
+        assert (run._device_boundaries[2].jets, run.jets[2]) == (jets, level_jets)
         assert run._status == "ready"
         run.set_jet_amplitude(0, .017, level=2)
         assert run.jets[2][0].amplitude == .017
@@ -136,8 +136,8 @@ def test_stepper_device_failure_keeps_host_and_blocks_reuse(where):
             run._stream.synchronize.side_effect = RuntimeError("upload failed")
         with pytest.raises(RuntimeError, match="upload failed"):
             run.set_jet_amplitude(0, .017, level=2)
-        np.testing.assert_array_equal(run._level_boundaries[2].tables["surface_velocity"], velocity)
-        assert (run._level_boundaries[2].jets, run.jets[2]) == (jets, level_jets)
+        np.testing.assert_array_equal(run._device_boundaries[2].tables["surface_velocity"], velocity)
+        assert (run._device_boundaries[2].jets, run.jets[2]) == (jets, level_jets)
         assert run._status == "failed"
         module.cuda.memcpy_htod.side_effect = None
         run._stream.synchronize.side_effect = None

tests/unit/test_multiple_pinned_remap.py

@@ -96,11 +96,11 @@ class MultiplePinnedRemap(unittest.TestCase):
         before = run.download()
         clone = run.remap_finest((PINS+MOVED+ADDED)[::-1], objects=bodies)
         try:
-            self.assertEqual(tuple(j.amplitude for j in clone._level_boundaries[2].jets), (.017, -.011))
-            np.testing.assert_array_equal(clone._level_boundaries[2].tables['surface_velocity'],
-                                          run._level_boundaries[2].tables['surface_velocity'])
+            self.assertEqual(tuple(j.amplitude for j in clone._device_boundaries[2].jets), (.017, -.011))
+            np.testing.assert_array_equal(clone._device_boundaries[2].tables['surface_velocity'],
+                                          run._device_boundaries[2].tables['surface_velocity'])
             clone.set_jet_amplitude(1, .005, level=2)
-            self.assertEqual(run._level_boundaries[2].jets[1].amplitude, -.011)
+            self.assertEqual(run._device_boundaries[2].jets[1].amplitude, -.011)
             clone.step()
             self.assertEqual((run.step_count, clone.step_count), (7, 8))
             for got, want in zip(run.download(), before):

tests/unit/test_pinned_body_remap.py

@@ -62,9 +62,9 @@ class PinnedBodyRemap(unittest.TestCase):
         run.set_jet_amplitude(0,.017,level=2)
         before=run.download()
         clone=run.remap_finest(B,objects=objects)
-        self.assertEqual(clone._level_boundaries[2].jets[0].amplitude,.017)
-        np.testing.assert_array_equal(clone._level_boundaries[2].tables['surface_velocity'],
-                                      run._level_boundaries[2].tables['surface_velocity'])
+        self.assertEqual(clone._device_boundaries[2].jets[0].amplitude,.017)
+        np.testing.assert_array_equal(clone._device_boundaries[2].tables['surface_velocity'],
+                                      run._device_boundaries[2].tables['surface_velocity'])
         self.assertTrue(same_physical_boundary(run.boundaries[2],clone.boundaries[2]))
         clone.step()
         self.assertEqual((run.step_count,clone.step_count),(7,8))
@@ -74,8 +74,8 @@ class PinnedBodyRemap(unittest.TestCase):
     def test_static_body_without_actuation(self):
         run,module,allocations,storage,objects=ready(jets=())
         clone=run.remap_finest(A,objects=objects)
-        self.assertFalse(clone._level_boundaries[2].actuated)
-        self.assertNotIn('surface_velocity',clone._level_boundaries[2].tables)
+        self.assertFalse(clone._device_boundaries[2].actuated)
+        self.assertNotIn('surface_velocity',clone._device_boundaries[2].tables)
         for a,b in zip(run.download(),clone.download()): np.testing.assert_array_equal(a,b)
         clone.close();run.close()
 

tests/unit/test_raw_exchange.py

@@ -85,7 +85,7 @@ class RawExchange(unittest.TestCase):
 
     def test_four_level_schedule_publication_and_fixed_storage(self):
         run, host, calls, arrays = backend()
-        b = run._level_boundaries[3]
+        b = run._device_boundaries[3]
         rows = len(b.tables['cell'])
         self.assertEqual(b.last_raw.size, rows * 16)
         self.assertEqual(b.interval_integral.size, rows * 32)
@@ -128,7 +128,7 @@ class RawExchange(unittest.TestCase):
 
     def test_disabled_has_nulls_and_no_observation_allocation_or_launch(self):
         run, host, calls, _ = backend(False)
-        b = run._level_boundaries[3]
+        b = run._device_boundaries[3]
         self.assertIs(b.last_raw, b.interval_integral)
         self.assertIs(b.interval_integral, None)
         run.step()

tests/unit/test_refinement_mrt.py

@@ -1,5 +1,6 @@
 """Host guards and source contract for MRT reuse in refined/blocked stepping (no CUDA)."""
 from pathlib import Path
+import re
 import unittest
 
 from CelerisLab.config import LBMConfig
@@ -32,22 +33,25 @@ class RefinementMRTHost(unittest.TestCase):
                      "collide_srt_no_force(f, feq, d_params.omega);"):
             self.assertIn(text, block_collision)
         self.assertIn("collide_srt(f, feq, forcing, omega);", block_collision)
-        self.assertIn("#elif BLOCK_ROOT_BC == 2", block_collision)
+        self.assertIn("#if COLLISION_MODEL == 2 || BLOCK_ROOT_BC == 2", block_collision)
         self.assertEqual(channel.count("block_collide(f);"), 1)
-        self.assertNotIn("collide_mrt_unforced(f);", channel)
+        # Executable text only: the header comment explains the absent force read.
+        code = re.sub(r"//[^\n]*", "", channel + "\n" + block_collision)
+        self.assertIn("d_params.omega", code)
+        for text in ("collide_mrt_unforced", "collide_dispatch(", "d_params.fx", "d_params.fy"):
+            self.assertNotIn(text, code)
         for source in (blocks,):
             self.assertIn("(COLLISION_MODEL != 0 && COLLISION_MODEL != 2)", source)
             self.assertIn("USE_LES", source)
             self.assertIn("STORE_PRECISION != 0", source)
             self.assertIn('#if COLLISION_MODEL == 2\n#include "operators/collision_mrt.cuh"\n#endif', source)
-        # One shared unforced MRT body; each blocked/patch path calls it exactly once.
-        mrt = (KERNELS / "operators/collision_mrt.cuh").read_text()
-        self.assertEqual(mrt.count("void collide_mrt_unforced("), 1)
-        for source in (blocks + "\n" + block_collision, channel + "\n" + block_collision):
-            self.assertEqual(source.count("collide_mrt_unforced(f);"), 1)
-            self.assertNotIn("collide_mrt_no_force(", source)
-            self.assertNotIn("compute_rho_u(f, rho, ux, uy);\n            collide_mrt", source)
-            self.assertNotRegex(source, r"compute_rho_u\(f, rho, ux, uy\);\s*collide_mrt")
+        # One unforced MRT call, in block_collide only; the steps never collide directly.
+        self.assertEqual(block_collision.count("collide_mrt_no_force(f, rho, ux, uy, omega);"), 1)
+        self.assertEqual(block_collision.count("compute_rho_u(f, rho, ux, uy);"), 1)
+        self.assertEqual(block_collision.count("fmaxf(OMEGA_COLLISION_MIN, d_params.omega)"), 1)
+        for source in (blocks, channel):
+            self.assertNotIn("collide_mrt", source)
+            self.assertNotIn("collide_srt", source)
 
 
 if __name__ == "__main__": unittest.main()

tests/unit/test_refinement_three_level_multipatch.py

@@ -141,7 +141,7 @@ class Gate(unittest.TestCase):
     def test_sparse_parent_tables_preserve_inherited_kinds(self):
         child,kind=r.nest_patch_set(CFG,self.parent,CHILDREN)
         layout=self.parent.fine_layout((4,4));order=np.argsort(child.coalesce_coarse)
-        donors,targets,bkind=host.blocked_base_tables(child,layout,order,composed_kind=kind)
+        donors,targets,bkind=host.blocked_coarse_tables(child,layout,order,composed_kind=kind)
         self.assertTrue((donors < layout.hole).all())
         self.assertTrue((targets < np.prod(layout.blocked_shape)).all())
         np.testing.assert_array_equal(bkind[layout.owner[layout.owned]],kind[layout.owned])
@@ -158,10 +158,10 @@ class Gate(unittest.TestCase):
         hole=int(np.flatnonzero(~layout.owned.ravel())[0])
         bad=child.ghost_coarse.copy();bad[0]=hole
         with self.assertRaisesRegex(ValueError,'owned uncovered'):
-            host.blocked_base_tables(replace(child,ghost_coarse=bad),layout,order,composed_kind=kind)
+            host.blocked_coarse_tables(replace(child,ghost_coarse=bad),layout,order,composed_kind=kind)
         bad=child.coalesce_coarse.copy();bad[0]=hole
         with self.assertRaisesRegex(ValueError,'owned uncovered'):
-            host.blocked_base_tables(replace(child,coalesce_coarse=bad),layout,order,composed_kind=kind)
+            host.blocked_coarse_tables(replace(child,coalesce_coarse=bad),layout,order,composed_kind=kind)
 
 
 def bodies_and_jets(speed=.002):

tests/unit/test_runtime_child_remap.py

@@ -1,5 +1,6 @@
 """CPU mocks execute production remap lifecycle; no GPU numerical claim."""
 from dataclasses import replace
+from types import SimpleNamespace
 import unittest
 from unittest.mock import patch
 
@@ -69,10 +70,11 @@ class RuntimeChildRemap(unittest.TestCase):
         run._status = 'ready'
         for bad in ([], [(40,16,48,40)], [(16,16,32,24),(16,27,32,55)]):
             with self.assertRaises(ValueError): run.remap_finest(bad)
-        hosts, pairs = run._boundaries, run.hierarchy.pairs
+        records, pairs = run._device_boundaries, run.hierarchy.pairs
         for attribute, value in (('_wall_aware', True),
-                                 ('_boundaries', [*hosts[:2], object()]),        # child boundary
-                                 ('_boundaries', [hosts[0], object(), hosts[2]]),  # fine boundary
+                                 # level-2 record whose host geometry differs from objects=()
+                                 ('_device_boundaries', [*records[:2], SimpleNamespace(host=object(), jets=())]),
+                                 ('_device_boundaries', [records[0], object(), records[2]]),  # level-1 body
                                  ('hierarchy', replace(run.hierarchy, pairs=(pairs[0], None)))):
             with patch.object(run, attribute, value):
                 with self.assertRaises(ValueError): run.remap_finest(B)

tests/unit/test_shifted_jets.py

@@ -67,7 +67,7 @@ class ShiftedJets(unittest.TestCase):
     def test_switch_upload_only_then_raw_remap_keeps_live_velocity(self):
         run,module,allocations,storage,objects=ready();clone=None
         try:
-            boundary=run._level_boundaries[2]
+            boundary=run._device_boundaries[2]
             before=run._download_raw()
             original=boundary.tables['surface_velocity'].copy()
             module.cuda.memcpy_htod.reset_mock()
@@ -79,7 +79,7 @@ class ShiftedJets(unittest.TestCase):
             for a,b in zip(before.values,run._download_raw().values):self.assertEqual(a.tobytes(),b.tobytes())
             clone=run.remap_finest(B,objects=objects)
             self.assertEqual(clone.jets[2][0].amplitude,-jet().amplitude)
-            np.testing.assert_array_equal(clone._level_boundaries[2].tables['surface_velocity'],
+            np.testing.assert_array_equal(clone._device_boundaries[2].tables['surface_velocity'],
                                           boundary.tables['surface_velocity'])
             self.assertEqual(clone.step_count,7)
             saved=run._download_raw()
@@ -92,12 +92,12 @@ class ShiftedJets(unittest.TestCase):
     def test_invalid_amplitude_does_not_upload_or_commit(self):
         run,module,_,_,_=ready()
         try:
-            before=run._level_boundaries[2].tables['surface_velocity'].copy()
+            before=run._device_boundaries[2].tables['surface_velocity'].copy()
             module.cuda.memcpy_htod.reset_mock()
             for amplitude in (float('nan'),float('inf'),1.):
                 with self.assertRaises(ValueError):run.set_jet_amplitude(0,amplitude,level=2)
             module.cuda.memcpy_htod.assert_not_called()
-            np.testing.assert_array_equal(before,run._level_boundaries[2].tables['surface_velocity'])
+            np.testing.assert_array_equal(before,run._device_boundaries[2].tables['surface_velocity'])
         finally:run.close()
 
 

tests/unit/test_single_level.py

@@ -73,7 +73,7 @@ class SingleLevel(unittest.TestCase):
 
     def test_circle_solids_are_not_state(self):
         run, _, _ = backend((circle(),))
-        solid = run._level_boundaries[0].tables["solid"]
+        solid = run._device_boundaries[0].tables["solid"]
         self.assertTrue(solid.any())
         np.testing.assert_array_equal(run._state[0], ~solid)
 

tests/unit/test_union_body_support.py

@@ -37,7 +37,7 @@ class UnionBodySupport(unittest.TestCase):
             run.set_jet_amplitude(0,.017,level=2)
             clone=run.remap_finest(((16,16,24,40),(24,16,32,40),A[1]),objects=objects)
             try:
-                self.assertEqual(clone._level_boundaries[2].jets[0].amplitude,.017)
+                self.assertEqual(clone._device_boundaries[2].jets[0].amplitude,.017)
                 self.assertTrue(same_physical_boundary(run.boundaries[2],clone.boundaries[2]))
                 for a,b in zip(run.download(),clone.download()):np.testing.assert_array_equal(a,b)
             finally:clone.close()

tests/validation/multigrid/test_block_stepper.py

@@ -240,7 +240,7 @@ class SingleLevelGPU(unittest.TestCase):
         with CudaContext(0) as context:
             for split in SPLITS:
                 with single(context=context, block=split, objects=[circle()]) as backend:
-                    tables = backend._level_boundaries[0].tables
+                    tables = backend._device_boundaries[0].tables
                     for key, values, dtype in (("q", [l.q for l in links], np.float32),
                                                ("dir", [l.dir for l in links], np.uint8),
                                                ("fallback", [l.fallback for l in links], np.uint8)):

tests/validation/multigrid/test_exchange_window.py

@@ -59,7 +59,7 @@ class ExchangeWindowCPU(unittest.TestCase):
             stack.enter_context(patch.object(cuda, "memcpy_dtoh",
                                              side_effect=lambda out, src: out.fill(0)))
             backend = build(m)
-            b = backend._level_boundaries[1]
+            b = backend._device_boundaries[1]
             self.assertIsNotNone(b.window_integral)
             with self.assertRaises(RuntimeError):
                 backend.start_exchange_window()  # not initialized
@@ -96,7 +96,7 @@ class ExchangeWindowCPU(unittest.TestCase):
         with mocked_backend() as m, ExitStack() as stack:
             stack.enter_context(patch.object(cuda, "memset_d8_async"))
             backend = build(m, window=False)
-            self.assertIsNone(backend._level_boundaries[1].window_integral)
+            self.assertIsNone(backend._device_boundaries[1].window_integral)
             backend.initialize(*[np.full(l.shape, 1 / 9, np.float32) for l in backend._levels])
             backend.step(1)
             self.assertEqual(len(accumulations(m)), 2)
@@ -120,7 +120,7 @@ class ExchangeWindowGPU(unittest.TestCase):
                 windowed.step(2)
                 plain.step(2)
                 windowed.start_exchange_window()
-                b = windowed._level_boundaries[1]
+                b = windowed._device_boundaries[1]
                 count = len(b.tables["cell"])
                 factors = np.array([.25, .25, .125, .25])  # h*dt, h*h*dt, h*h for level-1 h=.5
                 scale = np.zeros(4)
@@ -212,7 +212,7 @@ class ExchangeWindowShiftedMRTGPU(unittest.TestCase):
                     run.step(2)
                     if run is not plain:
                         run.start_exchange_window()
-                b = windowed._level_boundaries[1]
+                b = windowed._device_boundaries[1]
                 count = len(b.tables["cell"])
                 record = dict(scale=np.zeros(4), substeps=0)
                 audit_substeps(windowed, b, count, record)
@@ -228,7 +228,7 @@ class ExchangeWindowShiftedMRTGPU(unittest.TestCase):
                     same_bits(self, a, e)
                 rows = [np.empty((count, 4), np.float64) for _ in range(2)]
                 cuda.memcpy_dtoh(rows[0], b.window_integral)
-                cuda.memcpy_dtoh(rows[1], split._level_boundaries[1].window_integral)
+                cuda.memcpy_dtoh(rows[1], split._device_boundaries[1].window_integral)
                 np.testing.assert_array_equal(rows[0].view(np.uint64), rows[1].view(np.uint64))
 
                 out = windowed.read_exchange_window(0, level=1)

tests/validation/multigrid/test_multiple_pinned_remap_gpu.py

@@ -57,7 +57,7 @@ def make(context):
 def topology(run, rectangles, objects):
     return RemapTopology.from_config(
         CFG, run.hierarchy.pairs[0], rectangles, child_objects=objects,
-        child_jets=run._level_boundaries[2].jets,
+        child_jets=run._device_boundaries[2].jets,
     )
 
 
@@ -74,7 +74,7 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
     def change_amplitudes(self, run, amplitudes):
         # Each jet selects a different body. Updating one must preserve the
         # other body's current per-link velocities, including its driven links.
-        boundary = run._level_boundaries[2]
+        boundary = run._device_boundaries[2]
         body_ids = run.boundaries[2].curved_host[3]
         for index, amplitude in enumerate(amplitudes):
             before = boundary.tables['surface_velocity'].copy()
@@ -91,7 +91,7 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
         self.assertTrue(same_physical_boundary(new.boundaries[2], mesh.boundary))
         # Compare physical cut-link rows, never layout-dependent blocked addresses.
         for run in (old, new):
-            boundary = run._level_boundaries[2]
+            boundary = run._device_boundaries[2]
             self.assertEqual(tuple(jet.amplitude for jet in boundary.jets), amplitudes)
             exact(self, boundary.tables['surface_velocity'], mesh.surface_velocity,
                   'latest physical per-link velocities')

tests/validation/multigrid/test_pinned_body_remap_gpu.py

@@ -39,9 +39,9 @@ class PinnedBodyRemapGPU(unittest.TestCase):
             exact(self,getattr(after.boundaries[2],name),getattr(before.boundaries[2],name),name)
         for a,b in zip(after.boundaries[2].curved_host,before.boundaries[2].curved_host):
             exact(self,a,b,'physical cut links')
-        exact(self,after._level_boundaries[2].tables['surface_velocity'],
-              before._level_boundaries[2].tables['surface_velocity'],'current per-link velocities')
-        self.assertEqual(after._level_boundaries[2].jets[0].amplitude,.017)
+        exact(self,after._device_boundaries[2].tables['surface_velocity'],
+              before._device_boundaries[2].tables['surface_velocity'],'current per-link velocities')
+        self.assertEqual(after._device_boundaries[2].jets[0].amplitude,.017)
 
     def test_noop_next_step_bitwise(self):
         with CudaContext(0) as ctx, ExitStack() as owned:
@@ -63,7 +63,7 @@ class PinnedBodyRemapGPU(unittest.TestCase):
             current.set_jet_amplitude(0,.017,level=2)
             for rects in (B,A):
                 previous=current;before=previous.download();count=previous.step_count
-                jets=previous._level_boundaries[2].jets
+                jets=previous._device_boundaries[2].jets
                 old=RemapTopology.from_config(CFG,previous.hierarchy.pairs[0],previous.hierarchy.pairs[1].bounds,
                                              child_objects=objects,child_jets=jets)
                 new=RemapTopology.from_config(CFG,previous.hierarchy.pairs[0],rects,

tests/validation/multigrid/test_raw_exchange_gpu.py

@@ -105,7 +105,7 @@ class RawExchangeGPU(unittest.TestCase):
                     arrays = [state(l.shape) for l in measured._levels]
                     measured.initialize(*arrays)
                     plain.initialize(*arrays)
-                    b = measured._level_boundaries[1]
+                    b = measured._device_boundaries[1]
                     count = len(b.tables['cell'])
                     kernel = measured._levels[1].kernel
                     samples = []
@@ -158,7 +158,7 @@ class RawExchangeGPU(unittest.TestCase):
                  FixedPatchStepper(cfg, context=ctx, **settings) as plain:
                 arrays = [state(l.shape, ux=0.) for l in run._levels]
                 run.initialize(*arrays); plain.initialize(*arrays)
-                b = run._level_boundaries[3]
+                b = run._device_boundaries[3]
                 count = len(b.tables['cell'])
                 kernel = run._levels[3].kernel
                 samples, independent_samples = [], []
@@ -213,7 +213,7 @@ class RawExchangeGPU(unittest.TestCase):
                  FixedPatchStepper(CFG, context=ctx, **fine((16, 20)), observe_raw_exchange=True) as amr:
                 amr.initialize(*(state(l.shape) for l in amr._levels))
                 uniform.initialize(state((9, 48, 64)))
-                b = amr._level_boundaries[1]
+                b = amr._device_boundaries[1]
                 count = len(b.tables['cell'])
                 links = obj.geometry.build_cut_links(64, 48)
                 np.testing.assert_array_equal([v.dir for v in links], b.tables['dir'])
@@ -221,7 +221,7 @@ class RawExchangeGPU(unittest.TestCase):
                 # Same fine-grid force arms after translation to the same center: the
                 # single level's static BFL gets the AMR Jet column and force arms; the raw
                 # output pointer is substituted at launch (no stepper-side accumulation).
-                root = uniform._level_boundaries[0]
+                root = uniform._device_boundaries[0]
                 velocity = uniform._upload(b.tables['surface_velocity'])
                 root.device.update(surface_velocity=velocity, rx=uniform._upload(b.tables['rx']),
                                    ry=uniform._upload(b.tables['ry']))

tests/validation/multigrid/test_refinement_blocked_base.py

@@ -24,7 +24,8 @@ 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 CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+from CelerisLab.lbm.refinement_tables import blocked_coarse_tables
 from test_refinement_blocked_fine import NX, NY, BOUNDS, config, state, circle, actuation
 
 BASE_BLOCKS = ((24, 32), (12, 8), (3, 4), (1, 1))
@@ -108,7 +109,7 @@ class BlockedBaseCPU(unittest.TestCase):
                 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")
+                donors, targets, kind = blocked_coarse_tables(p, layout, order, mode != "periodic")
                 ly, lx = layout.local_shape
                 L = ly * lx
                 expected, source_cells = [], []
@@ -138,7 +139,7 @@ class BlockedBaseCPU(unittest.TestCase):
         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)
+        donors, targets, _ = blocked_coarse_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
@@ -178,7 +179,7 @@ class BlockedBaseCPU(unittest.TestCase):
         # 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)
+        donors, targets, _ = blocked_coarse_tables(p, single, order)
         initial = single.to_blocked(dense)
         staged = initial.copy()
         staged.ravel()[targets.astype(np.intp)] = arrivals

tests/validation/multigrid/test_refinement_blocked_fine.py

@@ -21,7 +21,8 @@ from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
 from CelerisLab.lbm.descriptors import SOLID
 from CelerisLab.lbm.refinement import (ACTIVE, INNER_GHOST, OUTER_GHOST, FixedPatch2D,
                                        prepare_fine_boundary)
-from CelerisLab.lbm.refinement_stepper import FixedPatchStepper, blocked_fine_tables
+from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+from CelerisLab.lbm.refinement_tables import blocked_fine_tables
 
 NX, NY, BOUNDS, STEPS = 32, 24, (8, 6, 24, 18), 8
 FINE_BLOCKS = ((32, 40), (16, 20), (16, 10))  # fine window 32x40: 1x1, 2x2, 4x2
@@ -186,8 +187,8 @@ 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._level_boundaries[1].flags, backend._levels[1].kind
-            velocity = backend._level_boundaries[1].device.get("surface_velocity")
+            flags, kind = backend._device_boundaries[1].flags, backend._levels[1].kind
+            velocity = backend._device_boundaries[1].device.get("surface_velocity")
             self.assertEqual(velocity is None, spec is None)
             backend.step(1)
             backend.close()
@@ -212,8 +213,8 @@ class BlockedFineCPU(unittest.TestCase):
     def test_actuation_without_velocity_buffer_fails_before_bfl(self):
         for speed in (0.0, .002):  # any spec, even zero amplitude
             backend, calls = self.mocked_backend([circle(actuation(speed))])
-            self.assertIsNotNone(backend._level_boundaries[1].tables["surface_velocity"])
-            del backend._level_boundaries[1].device["surface_velocity"]
+            self.assertIsNotNone(backend._device_boundaries[1].tables["surface_velocity"])
+            del backend._device_boundaries[1].device["surface_velocity"]
             with self.assertRaisesRegex(RuntimeError, "Surface velocity device buffer is missing"):
                 backend.step(1)
             names = [name for name, _ in calls]

tests/validation/multigrid/test_refinement_jets.py

@@ -62,8 +62,8 @@ class JetStepperCPU(unittest.TestCase):
         C = (16.0, 12.0)
         whole = Jet(Sector(C, .1, math.pi, 3.), Radial(C, inward=False), amplitude=0.)
         backend, _ = mocked([circle()], [whole], FINE_BLOCKS[1])
-        _, jet_id = backend._level_boundaries[1].jet_cache
-        self.assertGreater(len(np.unique(backend._level_boundaries[1].tables["block"][jet_id == 0])), 1)
+        _, jet_id = backend._device_boundaries[1].jet_cache
+        self.assertGreater(len(np.unique(backend._device_boundaries[1].tables["block"][jet_id == 0])), 1)
         backend.close()
 
     def test_fine_ownership_uses_remapped_physical_owner(self):
@@ -76,10 +76,10 @@ class JetStepperCPU(unittest.TestCase):
         source = SimObject(99, geometry=CircleGeometry(*C, 2.15), center=C, radius=2.15)
         whole = Jet(Sector(C, .1, math.pi, 3.), Radial(C, inward=False), amplitude=0.)
         backend, _ = mocked([source], [whole], FINE_BLOCKS[1])
-        _, jet_id = backend._level_boundaries[1].jet_cache
+        _, jet_id = backend._device_boundaries[1].jet_cache
         owners = backend.boundaries[1].curved_host[3][jet_id == 0]
         self.assertEqual(np.unique(owners).tolist(), [0])
-        self.assertGreater(len(np.unique(backend._level_boundaries[1].tables["block"][jet_id == 0])), 1)
+        self.assertGreater(len(np.unique(backend._device_boundaries[1].tables["block"][jet_id == 0])), 1)
         self.assertEqual(source.obj_id, 99)
         backend.close()
         cfg = config("periodic")
@@ -97,8 +97,8 @@ class JetStepperCPU(unittest.TestCase):
     def test_blocked_jet_velocity_matches_legacy_rows(self):
         legacy, _ = mocked([circle(actuation(SPEED))], (), FINE_BLOCKS[1])
         driven, _ = mocked([circle()], [jet()], FINE_BLOCKS[1])
-        a = legacy._level_boundaries[1].tables["surface_velocity"]
-        b = driven._level_boundaries[1].tables["surface_velocity"]
+        a = legacy._device_boundaries[1].tables["surface_velocity"]
+        b = driven._device_boundaries[1].tables["surface_velocity"]
         diff = np.abs(a - b).max()
         print(f"fine jet vs legacy velocity rows: max |dU| = {diff:.3e} (exact={diff == 0})")
         self.assertLessEqual(diff, 4 * np.finfo(np.float32).eps * SPEED)
@@ -106,8 +106,8 @@ class JetStepperCPU(unittest.TestCase):
 
     def test_blocked_amplitude_update_rewrites_only_velocity_column(self):
         backend, htod = mocked([circle()], [jet()], FINE_BLOCKS[1])
-        before = backend._level_boundaries[1].tables["surface_velocity"].copy()
-        device = backend._level_boundaries[1].device["surface_velocity"]
+        before = backend._device_boundaries[1].tables["surface_velocity"].copy()
+        device = backend._device_boundaries[1].device["surface_velocity"]
         with patch.object(cuda, "memcpy_htod", htod):
             htod.reset_mock()
             backend.set_jet_amplitude(0, -SPEED, level=1)
@@ -115,7 +115,7 @@ class JetStepperCPU(unittest.TestCase):
             self.assertEqual(htod.call_count, 1)
             self.assertIs(ptr, device)
             np.testing.assert_allclose(sent, -before, rtol=1e-6, atol=0)
-            np.testing.assert_array_equal(backend._level_boundaries[1].tables["surface_velocity"], sent)
+            np.testing.assert_array_equal(backend._device_boundaries[1].tables["surface_velocity"], sent)
             with self.assertRaises(ValueError):
                 backend.set_jet_amplitude(0, 1.5, level=1)
             self.assertEqual(htod.call_count, 1)
@@ -133,7 +133,7 @@ class JetStepperCPU(unittest.TestCase):
             with self.assertRaises(IndexError):
                 backend.set_jet_amplitude(1, -SPEED, level=1)
             self.assertEqual(backend._status, status)
-            before = backend._level_boundaries[1].tables["surface_velocity"].copy()
+            before = backend._device_boundaries[1].tables["surface_velocity"].copy()
             boom = Mock(side_effect=RuntimeError("copy failed"))
             with patch.object(cuda, "memcpy_htod", boom), patch.object(cuda, "memcpy_htod_async", boom), \
                     patch.object(cuda, "pagelocked_empty", side_effect=lambda shape, dtype: np.empty(shape, dtype)):
@@ -141,7 +141,7 @@ class JetStepperCPU(unittest.TestCase):
                     backend.set_jet_amplitude(0, -SPEED, level=1)
             self.assertEqual(backend._status, "failed", fine_block)
             self.assertEqual(backend.jets[1][0].amplitude, SPEED)
-            np.testing.assert_array_equal(backend._level_boundaries[1].tables["surface_velocity"], before)
+            np.testing.assert_array_equal(backend._device_boundaries[1].tables["surface_velocity"], before)
             with self.assertRaisesRegex(RuntimeError, "failed"):
                 backend.step(1)
             with self.assertRaisesRegex(RuntimeError, "failed"):

tests/validation/multigrid/test_shifted_jets_gpu.py

@@ -102,7 +102,7 @@ class Ledger:
 
     def step(self):
         run=self.run;self.samples.clear();run.step()
-        b=run._level_boundaries[self.index];h=run._exchange_spacing[self.index]
+        b=run._device_boundaries[self.index];h=run._exchange_spacing[self.index]
         self.test.assertEqual(len(self.samples),2**self.index)
         expected=np.zeros_like(self.samples[0],dtype=np.float64)
         weights=np.array([h*h,h*h,h*h*h,h*h])
@@ -140,7 +140,7 @@ class ShiftedJetsGPU(unittest.TestCase):
                         level_objects=[(),(obj,)],level_jets=[(),(whole,)],observe_raw_exchange=True) as run:
                     run._initialize_raw(initial(run,rho,u=.013,v=0.))
                     expected=run._download_raw();ledger=Ledger(self,run,1)
-                    self.assertTrue(run._level_boundaries[1].tables['aperture'].all())
+                    self.assertTrue(run._device_boundaries[1].tables['aperture'].all())
                     for _ in range(4):ledger.step()
                     error=max(float(abs(a[:,leaf]-b[:,leaf]).max()) for a,b,leaf in
                               zip(run._download_raw().values,expected.values,run.leaf))
@@ -201,7 +201,7 @@ class ShiftedJetsGPU(unittest.TestCase):
             self.exact_raw(before,old._download_raw())
             self.exact_raw(_StoredPopulations(expected.g,True),new._download_raw())
             self.assertEqual(new.jets[2][0].amplitude,-.001)
-            np.testing.assert_array_equal(new._level_boundaries[2].tables['surface_velocity'],old._level_boundaries[2].tables['surface_velocity'])
+            np.testing.assert_array_equal(new._device_boundaries[2].tables['surface_velocity'],old._device_boundaries[2].tables['surface_velocity'])
             newledger=Ledger(self,new,2);ledger.step();newledger.step()
             self.assertEqual((old.step_count,new.step_count),(count+1,count+1))
             ledger.report('L2 source');newledger.report('L2 remapped')

tests/validation/run_refined_dfg.py

@@ -81,7 +81,7 @@ def public_jets(center, radius, scenario):
     )
 
 
-def jet_interval(handle, scenario, begin, n, diameter, *, level=None):
+def jet_interval(handle, scenario, begin, n, diameter, level=None):
     """Apply a reversal only at a completed interval boundary, before advancing."""
     if scenario != "none":
         boundaries = [round(t * diameter / SPEED) for t in (250, 260, 270, 280, 290, 300, 305)]
@@ -92,9 +92,10 @@ def jet_interval(handle, scenario, begin, n, diameter, *, level=None):
     if begin < switch:
         return min(n, switch - begin), None
     if begin == switch:
-        target = {} if level is None else dict(level=level)
-        handle.set_jet_amplitude(0, -JET_PEAK, **target)
-        handle.set_jet_amplitude(1, JET_PEAK, **target)
+        # FixedPatchStepper needs an explicit Jet level; Simulation has one level.
+        kw = {} if level is None else {"level": level}
+        handle.set_jet_amplitude(0, -JET_PEAK, **kw)
+        handle.set_jet_amplitude(1, JET_PEAK, **kw)
         return n, dict(step=begin, tstar=300.0, amplitudes=[-JET_PEAK, JET_PEAK])
     return n, None
 
@@ -257,7 +258,7 @@ def host_geometry(source, h, collision="SRT", jet_scenario="none", placement="wi
         validate_channel_config,
         validate_channel_transfers,
     )
-    from CelerisLab.lbm.refinement_tables import blocked_base_tables
+    from CelerisLab.lbm import refinement_tables
     from CelerisLab.body.geometry.circle import CircleGeometry
     from CelerisLab.body.objects import SimObject
 
@@ -269,7 +270,15 @@ def host_geometry(source, h, collision="SRT", jet_scenario="none", placement="wi
     pair = FixedPatch2D.from_config(cfg, PLACEMENTS[placement])
     layout = BlockLayout2D.build(cfg.ny, cfg.nx, *BASE_BLOCK)
     validate_channel_config(cfg, [BASE_BLOCK, FINE_BLOCK])
-    targets = blocked_base_tables(pair, layout, np.argsort(pair.coalesce_coarse))[1]
+    # Select by the imported source's own API, not by scenario: the pinned historical
+    # source (verify_source) defines blocked_base_tables, the current checkout (host
+    # tests, public Jet scenarios) blocked_coarse_tables. Exactly one must exist.
+    names = [n for n in ("blocked_base_tables", "blocked_coarse_tables")
+             if hasattr(refinement_tables, n)]
+    if len(names) != 1:
+        raise RuntimeError(f"Unexpected coarse transfer table API: {names}")
+    coarse_tables = getattr(refinement_tables, names[0])
+    targets = coarse_tables(pair, layout, np.argsort(pair.coalesce_coarse))[1]
     validate_channel_transfers(pair, layout, targets)
     BlockLayout2D.build(*pair.fine_shape, *FINE_BLOCK, periodic=(False, False))
     body = SimObject(