refactor: remove redundant periodic patch stepper alias

Retired with the dense refinement route (results recorded at 1573942 and
2469af7 or named as accepted losses):
- test_refinement_blocked_base: BlockedBaseGPU exact parity (3 tests),
  DenseBlockedRouteParity (6 tests incl. MRT and two bodies),
  DenseBlockedRouteCPU (capability matrix, ledger comparator);
- test_refinement_blocked_fine: BlockedFineGPU dense-vs-blocked (3 tests);
- test_refinement_mrt_gpu: test_dense_versus_blocked_mrt_two_level;
- test_raw_exchange_gpu: test_dense_fine_storage_reference (private
  link_exchange cut-link oracle, accepted loss);
- test_refinement_gpu: test_bounded_equivalence_and_ownership (dense
  CELERIS_REFINEMENT_AUDIT ownership counters, accepted loss); its physics
  stays in test_bounded_equivalence_blocked; eight_steps loses the audit
  branch and the misleading 'ownership/counts PASS' label;
- test_fine_surface_gpu: test_fine_surface_four_hops (dense per-hop flux
  ledger, CurvedLinkSoA reflection oracle, audit counts: accepted losses;
  raw exchange covers the boundary part only, not full flux closure) and the
  dense SoA fault tests (missing surface buffer, close frees SoA); blocked
  equivalents remain in test_refinement_blocked_fine/_base;
- test_raw_exchange.test_dense_observation_rejected_before_device_access.
New TileContractCPU: level_blocks omitted/None/wrong length/None entry/
wrong arity/non-integer/bool/non-positive/non-tuple/not dividing all raise
before CUDA; whole-level tiles reach CUDA. Rejection tests and fixtures pass
tiles; launch-order and kernel-text tests follow the blocked root;
examples/fixed_periodic_patch uses whole-level tiles (+ host tile test).

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

docs/compatibility.md

@@ -3,6 +3,13 @@
 Changes that can make a previously accepted setup or call behave differently.
 Each entry names the commit, the affected entry points and the observable result.
 
+## Redundant periodic stepper alias removed
+
+`PeriodicPatchStepper` is no longer exported from `lbm.refinement_stepper`.
+Use `FixedPatchStepper` with the same arguments, including explicit
+`level_blocks`. The old name was an alias for that exact class; periodic
+boundaries remain the default. Existing imports of the old name must be updated.
+
 ## 87c3a31: cut-link ownership is validated for every body
 
 `validate_cut_link_ownership` (`body/preprocess/flag_overlay.py`) now runs for

docs/lifecycle.md

@@ -5,7 +5,7 @@ it is; it adds no framework. Line references are to candidate 0f8fdbc3.
 
 | Class | Contents | How it changes | Owner rule |
 |---|---|---|---|
-| 1. Module and its globals | Compile-time (visible to the compiler, fixed per module): `config.h` macros incl. `NX`/`NY`, `COLLISION_MODEL`, `BLOCK_ROOT_BC`, `WALL_*`. Runtime per-module globals (not visible to the compiler; written by host code): `__constant__` `d_params`, `d_w`, `d_cx`, `d_cy`, `d_block_channel`, `d_block_wall`; plus bound kernel handles | Globals are written once by the owner right after `load_module`. Uniform grid: `Simulation._recompile` (construction, body-topology sync) loads a new module and re-uploads; `update_runtime_params` rewrites `d_params` (omega, force, n_objects) between runs. Steppers: never rewritten after construction, so they are immutable for the instance's lifetime although they are runtime values. The same quantity can be a macro in one route and a global in another (dense slot root: `WALL_SLOT_WIDTH`/`WALL_NORMAL_VELOCITY` macros; blocked slot root: `d_block_wall`), which can change generated arithmetic. | Every owner loads its own module (`compiler.load_module` = `module_from_file`, no cache). A loaded module is never shared between live owners: globals are per module, so sharing would share mutable `d_params`. |
+| 1. Module and its globals | Compile-time (visible to the compiler, fixed per module): `config.h` macros incl. `NX`/`NY`, `COLLISION_MODEL`, `BLOCK_ROOT_BC`, `WALL_*`. Runtime per-module globals (not visible to the compiler; written by host code): `__constant__` `d_params`, `d_w`, `d_cx`, `d_cy`, `d_block_channel`, `d_block_wall`; plus bound kernel handles | Globals are written once by the owner right after `load_module`. Uniform grid: `Simulation._recompile` (construction, body-topology sync) loads a new module and re-uploads; `update_runtime_params` rewrites `d_params` (omega, force, n_objects) between runs. Steppers: never rewritten after construction, so they are immutable for the instance's lifetime although they are runtime values. The blocked slot/reservoir root reads its slot from the `d_block_wall` constant (the former dense route used `WALL_*` macros; both were bit-identical in the parity gates). | Every owner loads its own module (`compiler.load_module` = `module_from_file`, no cache). A loaded module is never shared between live owners: globals are per module, so sharing would share mutable `d_params`. |
 | 2. Geometry, layout, tables | flags, curved-link tables, block layouts and halo pairs, transfer (explode/coalesce) tables, block kinds | Immutable for the life of an instance. Uniform grid: rebuilt by body sync. `FixedPatchStepper`: a geometry change is a new instance (`remap_child`), committed only after a successful transfer. | No in-place table edits. |
 | 3. Controls and dependent boundary columns | Jet amplitude -> per-link `surface_velocity`; body actions (uniform grid) | Between steps, without a table or module rebuild (`set_jet_amplitude`). | Host copy commits only after the device upload succeeded; a device error fails the instance. |
 | 4. Per-step state | population double buffers and their parity, raw-exchange ledgers, observations | Device-only during `step(n)` / `run(n)`; buffer swaps are host-handle swaps. | Valid on the host only after the single wait at the end of the call. |
@@ -15,10 +15,11 @@ Storage routes are separate concerns:
 - Original uniform grid (`uniform.cu`, `Simulation`): whole-domain `NX`/`NY`.
 - Blocked levels (`blocked.cu`; every `FixedPatchStepper` level, including the
   single level of `refinements=[]`): `NX`/`NY` are one block window including halo.
-- Dense refinement levels (`refinement.cu`, `level_blocks[l] is None`).
+- Dense refinement levels (`refinement.cu`) were removed: `level_blocks` is
+  required and every refinement level is blocked.
 
 A single-rectangle refinement (`FixedPatch2D`) is stored in its own fine
-window frame in both dense and blocked form (`refinement_stepper.py:222-224`,
+window frame (`FixedPatchStepper.__init__`,
 `BlockLayout2D.build(*pair.fine_shape, ...)`). Only a rectangle set
 (`FixedPatchSet2D`) uses the level-wide frame `(2*ny, 2*nx)`, whose
 `fine_layout` requires that full grid to be divisible by the tile

docs/supported_settings.md

@@ -10,12 +10,13 @@ measured; they are not general bounds.
 
 | Setting | Supported | Guard |
 |---|---|---|
-| Dimension / lattice | 2D D2Q9 only (3D refinement is future work) | `lbm/refinement.py` `FixedPatch2D.from_config`; `#error` guards in `lbm/kernels/blocks.cu` and `lbm/kernels/refinement.cu` |
+| Dimension / lattice | 2D D2Q9 only (3D refinement is future work) | `lbm/refinement.py` `FixedPatch2D.from_config`; `#error` guards in `lbm/kernels/blocked.cu` |
 | Precision / streaming | FP32 storage, double buffer | same guards |
 | Collision on refinement levels | SRT or MRT; TRT and LES rejected | `lbm/refinement.py` `FixedPatch2D.from_config` (re-checked for fine boundaries in `prepare_fine_boundary`) |
-| Shifted (ddf_shifting) storage | only when every level is blocked; dense or mixed rejected | `lbm/refinement_stepper.py` `FixedPatchStepper.__init__`, `FixedPatchStepper._check_storage`; `#error` guard in `lbm/kernels/refinement.cu` (dense patch kernel unshifted only) |
+| Shifted (ddf_shifting) storage | supported (every level is blocked); one representation on all levels | `lbm/refinement_stepper.py` `FixedPatchStepper._check_storage` |
 | Fine relaxation | 1/tau_fine must lie inside the configured omega guards | `lbm/refinement.py` `FixedPatch2D.from_config` |
 | Single blocked level `FixedPatchStepper(refinements=[], level_blocks=[(by, bx)], level_objects=[circles])` (replaces `PeriodicBlockStepper`, removed in 8e1feff) | doubly periodic, SRT only, no LES/Jets/projection/exchange observation, stationary passive circles; `download()` returns a 1-tuple | `lbm/refinement_stepper.py` `_check_levels`, `FixedPatchStepper.__init__`; `lbm/blocks.py` `static_circle_tables` |
+| Refinement storage | blocked on every level: `level_blocks` is required, one tile `(by, bx)` per level dividing that level; a whole-level tile (`(ny, nx)` at the root, the fine window shape for one rectangle) always divides. Dense and mixed dense/blocked storage were removed (parity of D/D, D/B and B/B recorded at 1573942 and 2469af7) | `lbm/refinement_stepper.py` `_check_levels`, `FixedPatchStepper.__init__` |
 
 Tests: `tests/unit/test_refinement_mrt.py`, `tests/validation/multigrid/test_refinement_mrt_gpu.py`,
 `tests/unit/test_shifted_storage.py`, `tests/unit/test_shifted_block_operators.py`.

examples/fixed_periodic_patch.py

@@ -16,7 +16,7 @@ from CelerisLab.config import load_lbm_config
 from CelerisLab.cuda.context import CudaContext
 from CelerisLab.lbm.equilibrium import compute_feq_d2q9, compute_macro_from_ddf
 from CelerisLab.lbm.refinement import ACTIVE
-from CelerisLab.lbm.refinement_stepper import PeriodicPatchStepper
+from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
 
 
 def shear_equilibrium(shape, origin, dx, period):
@@ -48,7 +48,8 @@ def main():
                   nx=16, ny=16, viscosity=.1, velocity=1e-3, rho=1.0,
                   lattice_model="D2Q9", collision="SRT", streaming="double_buffer",
                   store_precision="FP32", ddf_shifting=False, les_enabled=False)
-    with CudaContext(0) as context, PeriodicPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)]) as run:
+    with CudaContext(0) as context, FixedPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)],
+                                                           level_blocks=[(16, 16), (24, 24)]) as run:
         p = run.patch
         coarse = shear_equilibrium(p.coarse_shape, (0, 0), 1.0, cfg.ny)
         fine = shear_equilibrium(p.fine_shape, p.fine_origin, .5, cfg.ny)

src/CelerisLab/cuda/compiler.py

@@ -7,8 +7,7 @@ build_module runs, in a private directory it removes afterwards:
   1. build_dir = generate_config(lbm_cfg, n_objects) — isolated source snapshot
   2. compile_kernel(arch, source=entry, build_dir=build_dir) — nvcc -ptx
   3. load_module(ptx_path)              — PyCUDA module from PTX
-Entries: uniform.cu (uniform grid), blocked.cu (blocked level), refinement.cu
-(dense refinement level).
+Entries: uniform.cu (uniform grid), blocked.cu (every blocked level, single or refined).
 """
 
 from __future__ import annotations

src/CelerisLab/lbm/kernels/refinement.cu

@@ -1,40 +0,0 @@
-// Fixed, unforced internal 2:1 patch. Separate compilation entry;
-// default OneStep/channel kernels are deliberately not included.
-#include <stdint.h>
-#include <math.h>
-#include "config.h"
-#include "core/precision.cuh"
-#include "core/descriptors.cuh"
-#include "core/layout.cuh"
-#include "core/params.cuh"
-#include "core/flags.cuh"
-#include "core/obs.cuh"
-#include "core/cell_kind.cuh"
-#include "operators/macro.cuh"
-#include "operators/collision_srt.cuh"
-#if COLLISION_MODEL == 2
-#include "operators/collision_mrt.cuh"
-#endif
-#include "streaming/pull_double_buffer.cuh"
-#include "boundary/planar_velocity.cuh"
-#include "boundary/curved_boundary.cuh"
-
-extern "C" {
-#include "step/curved_double.cu"
-}
-
-#if DIM != 2 || NQ != 9 || (COLLISION_MODEL != 0 && COLLISION_MODEL != 2) || STORE_PRECISION != 0 || USE_DDF_SHIFTING || USE_LES || STREAMING_MODEL != 0
-#error "Refinement requires unshifted FP32 D2Q9 SRT or MRT double buffers without LES"
-#endif
-
-// Enabled only in the validation build; no diagnostic arrays in normal builds.
-#ifdef CELERIS_REFINEMENT_AUDIT
-extern "C" {
-__device__ unsigned int collision_counts[TOTAL_CELLS];
-__device__ unsigned int explosion_counts[TOTAL_CELLS];
-__device__ unsigned int incoming_counts[NQ * TOTAL_CELLS];
-}
-#endif
-
-// Kernels: patch collision helper, PatchExplode, PatchFineStep, PatchCoarseStep.
-#include "step/patch_step.cu"

src/CelerisLab/lbm/kernels/step/block_transfer.cu

@@ -16,8 +16,8 @@ extern "C" __global__ void BlockHaloCopy(
 }
 
 // Fixed-patch explosion into blocked fine storage: target = owner block cell
-// 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.
+// id of a patch-ghost cell; donors are parent block cell ids and coarse_cells
+// is the parent's block window size.
 extern "C" __global__ void BlockExplode(
     const fpxx* coarse, fpxx* fine, const uint64_t* target, const uint64_t* donor,
     unsigned int count, unsigned long coarse_cells)
@@ -66,7 +66,7 @@ extern "C" __global__ void BlockExplodeProjected(
 }
 
 // Targets are receiver-local SOURCE slots for the immediately following pull.
-// Keep child summation order identical to PatchCoarseStep; no collision here.
+// Child summation order ((c0 + c1) + c2) + c3 is part of the numerics; no collision here.
 extern "C" __global__ void BlockCoalesce(
     const fpxx* fine, fpxx* coarse, const uint64_t* target, const uint64_t* children,
     unsigned int count)

src/CelerisLab/lbm/kernels/step/patch_step.cu

@@ -1,95 +0,0 @@
-// CelerisLab/lbm/kernels/step/patch_step.cu
-// Dense fixed-patch collision, explosion and the fine/coarse steps.
-// Part of the refinement.cu compilation entry; relies on its includes, guards
-// and the optional CELERIS_REFINEMENT_AUDIT counters declared there.
-
-__device__ __forceinline__ void patch_collide(float* f, unsigned long k) {
-#if COLLISION_MODEL == 2
-    collide_mrt_unforced(f);
-#else
-    float rho, ux, uy, feq[NQ];
-    compute_rho_u(f, rho, ux, uy);
-    compute_feq(rho, ux, uy, feq);
-    collide_srt_no_force(f, feq, d_params.omega);
-#endif
-#ifdef CELERIS_REFINEMENT_AUDIT
-    ++collision_counts[k];
-#endif
-}
-
-extern "C" __global__ void PatchExplode(
-    const float* coarse, float* fine, const uint32_t* targets,
-    const uint32_t* donors, unsigned int count, unsigned long coarse_cells)
-{
-    unsigned int row = blockIdx.x * blockDim.x + threadIdx.x;
-    if (row >= count) return;
-    for (unsigned int q = 0; q < NQ; ++q)
-        fine[index_f(targets[row], q)] = coarse[q * coarse_cells + donors[row]];
-#ifdef CELERIS_REFINEMENT_AUDIT
-    ++explosion_counts[targets[row]];
-#endif
-}
-
-extern "C" __global__ void PatchFineStep(const float* in, float* out, const uint16_t* flags) {
-    unsigned int x, y; unsigned long k;
-    index_from_thread(x, y, k);
-    if (x >= NX || y >= NY) return;
-    float f[NQ];
-    // Outer halo is invalid after EACH hop, including when reusing the buffer
-    // that held explosion data. Inner halo only streams; it never collides.
-    if (x < 2 || y < 2 || x >= NX - 2 || y >= NY - 2
-            || (flags != nullptr && is_solid(flags[k]))) {
-        fill_invalid(f);
-    } else {
-        unsigned long j[NQ];
-        compute_neighbors(x, y, j);  // These interior targets never wrap.
-        stream_pull_load(k, f, in, j);
-        if (x >= 4 && y >= 4 && x < NX - 4 && y < NY - 4)
-            patch_collide(f, k);
-    }
-    stream_pull_store(k, f, out);
-}
-
-extern "C" __global__ void PatchCoarseStep(
-    const float* in, float* out, const float* fine,
-    const uint64_t* targets, const uint64_t* children, unsigned int count,
-    unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1,
-    unsigned int wall_aware)
-{
-    unsigned int x, y; unsigned long k;
-    index_from_thread(x, y, k);
-    if (x >= NX || y >= NY) return;
-    // Both buffers were poisoned at initialization; covered cells stay invalid.
-    if (x >= x0 && x < x1 && y >= y0 && y < y1) return;
-    if (wall_aware && (y == 0 || y == NY - 1)) return;
-    unsigned long j[NQ];
-    compute_neighbors(x, y, j);
-    float f[NQ];
-    stream_pull_load(k, f, in, j);
-    if (x >= x0 - 1 && x <= x1 && y >= y0 - 1 && y <= y1) {
-        for (unsigned int q = 1; q < NQ; ++q) {
-            uint64_t slot = index_f(k, q);
-            unsigned int lo = 0, hi = count;
-            while (lo < hi) {
-                unsigned int mid = lo + (hi - lo) / 2;
-                if (targets[mid] < slot) lo = mid + 1;
-                else hi = mid;
-            }
-            if (lo < count && targets[lo] == slot) {
-                const uint64_t* s = children + 4ul * lo;
-                f[q] = 0.25f * (fine[s[0]] + fine[s[1]] + fine[s[2]] + fine[s[3]]);
-#ifdef CELERIS_REFINEMENT_AUDIT
-                ++incoming_counts[slot];
-#endif
-            }
-        }
-    }
-    // Includes the eight grazing-corner links with UNCOVERED coarse donors.
-    // All replacement arrivals must be in place before computing rho/u.
-    if (wall_aware) {
-        if (y == 1) reconstruct_bottom_velocity(f, 0.0f, bottom_wall_normal_velocity(x));
-        else if (y == NY - 2) reconstruct_top_reservoir(f, RHO);
-    }
-    patch_collide(f, k);
-    stream_pull_store(k, f, out);
-}

src/CelerisLab/lbm/refinement_stepper.py

@@ -19,7 +19,7 @@ from .blocks import (BlockLayout2D, launch_curved_boundary, launch_halo_copy, le
                      static_circle_tables)
 from .channel import channel_masks, validate_channel_config, validate_channel_transfers
 from ..body.jets import jet_link_velocity
-from .hierarchy import build_hierarchy, is_rectangle
+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)
                                 blocked_kind, projected_explosion_tables, _prepare_blocked_boundary,
@@ -75,21 +75,19 @@ def _check_levels(refinements, blocks, objects, jets):
         if len(values) != depth:
             raise ValueError(f"{name} needs one entry per level ({depth})")
     for level, block in enumerate(blocks):
-        if block is not None and len(block) != 2:
-            raise ValueError(f"level_blocks[{level}] must be (by, bx)")
+        if (not isinstance(block, tuple) or len(block) != 2
+                or any(isinstance(v, (bool, np.bool_)) or not isinstance(v, Integral) or v <= 0
+                       for v in block)):
+            raise ValueError(f"level_blocks[{level}] must be a tile (by, bx) of positive integers; "
+                             "a whole-level tile ((ny, nx) at the root, the fine window shape "
+                             "for one rectangle) always divides its level")
     if depth == 1:  # one blocked level: the root carries the static circles
-        if blocks[0] is None or jets[0]:
-            raise ValueError("a single level must be blocked and has no Jets")
+        if jets[0]:
+            raise ValueError("a single level has no Jets")
         return
     if objects[0] or jets[0]:
         raise ValueError("base-level bodies/Jets are not supported")
-    if any(blocks[level - 1] is not None and blocks[level] is None for level in range(1, depth)):
-        raise ValueError("a blocked level requires blocked refined levels below it")
-    if not is_rectangle(refinements[0]) and (blocks[0] is None or blocks[1] is None):
-        raise ValueError("several L1 patches require level_blocks[0] and level_blocks[1]")
     if depth > 2:
-        if any(block is None for block in blocks):
-            raise ValueError("three or more levels must all be blocked")
         if any(objects[level] or jets[level] for level in range(1, depth - 1)):
             raise ValueError("bodies/Jets are supported only on the finest level of a nested hierarchy")
 
@@ -101,15 +99,16 @@ class FixedPatchStepper:
     (x0, y0, x1, y1) (FixedPatch2D, a fine window with its own origin) or a
     sequence of rectangles (FixedPatchSet2D, a dependency-validated union on one
     level-wide grid, public fine arrays (2*ny, 2*nx), NaN outside ACTIVE cells);
-    see lbm.hierarchy. ``level_blocks``, ``level_objects`` and ``level_jets``
-    hold one entry per level, L0 included (L0 bodies/Jets are not supported).
+    see lbm.hierarchy. ``level_blocks`` (required), ``level_objects`` and
+    ``level_jets`` hold one entry per level, L0 included (L0 bodies/Jets are not
+    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).
 
     ``coarse_boundary`` is ``periodic`` (default) or ``slot_reservoir``:
     a bottom on-site slot and top fixed-density reservoir from ``cfg``.
-    ``channel`` requires all-blocked levels and matches uniform SRT Zou-He
+    ``channel`` matches uniform SRT Zou-He
     inlet, pressure-NEQ outlet and stationary bounce-back y walls. Boundary
     ghosts evolve internally but public downloads exclude them; initialization
     must supply their finite positive-density state. Transfers must avoid the
@@ -132,14 +131,13 @@ class FixedPatchStepper:
     stream. On a device/launch failure the instance becomes unusable: close and
     recreate it, rather than pretending a partially advanced buffer rolled back.
 
-    A ``level_blocks`` entry (by, bx) stores that level in blocked.cu storage:
-    the ENTIRE level shape, including the four-cell patch ghost band, is
-    partitioned by a non-periodic BlockLayout2D with one-cell same-level halos.
-    Patch kinds stay per cell (OUTER invalid, INNER stream only, ACTIVE
-    collide); halos are gathered before each hop. Public arrays keep the dense
-    shapes; only the allocation differs. A blocked level requires blocked
-    refined levels below it; several L1 rectangles and depth >= 3 require all
-    levels blocked. Blocked bodies are stationary circles, optionally
+    Every level is stored in blocked.cu storage with its ``level_blocks`` tile
+    (by, bx): the ENTIRE level shape, including the four-cell patch ghost band,
+    is partitioned by a BlockLayout2D with one-cell same-level halos; the tile
+    must divide the level (a whole-level tile always does). Patch kinds stay
+    per cell (OUTER invalid, INNER stream only, ACTIVE collide); halos are
+    gathered before each hop. Public arrays are dense (9, ny, nx); only the
+    allocation is blocked. Bodies are stationary circles, optionally
     surface-actuated: prepare_fine_boundary decides legality on the dense grid,
     then cut links are staged receiver-local per hop: halo -> BFL -> step, same
     buffer. A blocked parent gathers halos before receiver-local coalescence
@@ -148,8 +146,7 @@ class FixedPatchStepper:
     ``observe_raw_exchange=True`` records blocked physical cut-link exchange
     without changing populations. read_raw_exchange reports the last complete
     base interval or last substep in base-lattice units per unit depth. This
-    is boundary momentum reaction, not a validated nozzle/stress load. Dense
-    body levels are not supported by this diagnostic.
+    is boundary momentum reaction, not a validated nozzle/stress load.
     ``observe_window_exchange=True`` (requires observe_raw_exchange) adds a second
     per-link FP64 ledger accumulated by the same kernel between
     start_exchange_window() and read_exchange_window(); disabled by default.
@@ -180,30 +177,32 @@ class FixedPatchStepper:
                  level_blocks=None, level_objects=None, level_jets=None,
                  observe_raw_exchange=False, observe_window_exchange=False, explosion="constant"):
         if HALO != 4:
-            raise RuntimeError("refinement.cu requires a four-cell fine halo")
+            raise RuntimeError("the refinement transfers require a four-cell fine halo")
         if coarse_boundary not in ("periodic", "slot_reservoir", "channel"):
             raise ValueError("coarse_boundary must be periodic, slot_reservoir or channel")
         depth = len(refinements) + 1
-        level_blocks = [None] * depth if level_blocks is None else list(level_blocks)
+        if level_blocks is None:
+            raise ValueError("level_blocks is required: one tile (by, bx) per level; a whole-level "
+                             "tile ((ny, nx) at the root, the fine window shape for one rectangle) "
+                             "always divides its level")
+        level_blocks = list(level_blocks)
         level_objects = [()] * depth if level_objects is None else list(level_objects)
         level_jets = [()] * depth if level_jets is None else list(level_jets)
         level_objects = [tuple(o) for o in level_objects]
         level_jets = [tuple(j) for j in level_jets]
-        level_blocks = [None if b is None else tuple(b) for b in level_blocks]
+        level_blocks = [tuple(b) if isinstance(b, (list, tuple)) else b for b in level_blocks]
         _check_levels(refinements, level_blocks, level_objects, level_jets)
         depth = len(refinements) + 1
         if explosion not in ("constant", "projected_linear"):
             raise ValueError("explosion must be constant or projected_linear")
         if explosion == "projected_linear" and (
-                coarse_boundary != "periodic" or any(b is None for b in level_blocks)
+                coarse_boundary != "periodic"
                 or any(level_objects) or any(level_jets)
                 or cfg.collision != "SRT" or cfg.les_enabled):
-            raise ValueError("projected explosion requires periodic all-blocked SRT levels "
+            raise ValueError("projected explosion requires periodic SRT levels "
                              "without bodies/Jets or LES")
         self._explosion = explosion
         self._ddf_shifting = bool(cfg.ddf_shifting)
-        if self._ddf_shifting and any(b is None for b in level_blocks):
-            raise ValueError("shifted populations require all-blocked levels; dense/mixed is unsupported")
         if depth == 1 and (coarse_boundary != "periodic" or cfg.collision != "SRT" or cfg.les_enabled
                            or explosion != "constant"
                            or observe_raw_exchange or observe_window_exchange):
@@ -225,18 +224,16 @@ class FixedPatchStepper:
                 raise ValueError("slot_reservoir requires y0 >= 4 and y1 <= ny-4")
         self._layouts = []
         for level, block in enumerate(level_blocks):
-            layout = None
-            if block is not None:
-                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_name(level)} blocks exceed launch limits")
+            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_name(level)} blocks exceed launch limits")
             self._layouts.append(layout)
         if self._channel:
             base = blocked_base_tables(pairs[0], self._layouts[0], np.argsort(pairs[0].coalesce_coarse))
@@ -274,9 +271,6 @@ class FixedPatchStepper:
             spacings.append(spacing)
         if not isinstance(observe_raw_exchange, bool):
             raise ValueError("observe_raw_exchange must be bool")
-        if observe_raw_exchange and any(b is not None and layout is None
-                                      for b, layout in zip(self._boundary_hosts, self._layouts)):
-            raise ValueError("raw exchange currently requires blocked body levels")
         self._observe_raw_exchange = observe_raw_exchange
         if not isinstance(observe_window_exchange, bool):
             raise ValueError("observe_window_exchange must be bool")
@@ -319,8 +313,6 @@ class FixedPatchStepper:
         self._allocations = []
         self._levels, self._pairs = [], []
         self._stream = None
-        self._fine_flags = self._fine_curved = self._curved = None
-        self._link_tables = {}
         self._block = (cfg.threads_per_block, 1, 1)
         with self.ctx.activate():
             try:
@@ -348,23 +340,6 @@ class FixedPatchStepper:
                                 b.window_integral = self._allocate(count * 4 * 8)
                             b.accumulate = level.module.get_function("BlockAccumulateExchange")
                         b.kernel = level.module.get_function("BlockCurvedBoundary")
-                if depth == 1:
-                    pass
-                elif self._level_boundaries[1] is not None:
-                    b = self._level_boundaries[1]
-                    self._fine_flags, self._link_tables, self._curved = b.flags, b.device, b.kernel
-                elif self._boundary_hosts[1] is not None:
-                    from .curved_links import CurvedLinkSoA
-                    boundary = self._boundary_hosts[1]
-                    self._fine_flags = self._upload(boundary.flags)
-                    # SoA owns its allocations; never also register them in _allocations.
-                    self._fine_curved = CurvedLinkSoA()
-                    self._fine_curved.assign_host(*boundary.curved_host)
-                    self._fine_curved.assign_surface(boundary.objects, fine_cfg, boundary.flags,
-                                                     jets=self._level_jets[1],
-                                                     hits=boundary.curved_hits)
-                    self._fine_curved.upload()
-                    self._curved = self._levels[1].module.get_function("CurvedBoundaryKernel")
             except BaseException:
                 try:
                     self.close()
@@ -380,11 +355,8 @@ class FixedPatchStepper:
         """
         total = 0
         for index, layout in enumerate(self._layouts):
-            cells = int(np.prod(layout.blocked_shape if layout is not None
-                                else (9, *self.hierarchy.kinds[index].shape)))
-            total += 2 * 4 * cells
-            if layout is not None:
-                total += layout.halo_dst.nbytes + layout.halo_src.nbytes
+            total += 2 * 4 * int(np.prod(layout.blocked_shape))
+            total += layout.halo_dst.nbytes + layout.halo_src.nbytes
             b = self._level_boundaries[index]
             if b is not None:
                 total += sum(np.asarray(v).nbytes for k, v in b.tables.items()
@@ -397,83 +369,62 @@ class FixedPatchStepper:
         pairs = self.hierarchy.pairs
         tau = (pairs[0].tau_coarse if pairs else 3 * config.viscosity + .5) if index == 0 \
             else pairs[index - 1].tau_fine
-        name, source = ("PatchCoarseStep" if index == 0 else "PatchFineStep"), "refinement.cu"
-        if layout is not None:
-            config = level_config(config, layout, tau)
-            name, source = ("BlockChannelStep" if index == 0 and self._channel else "BlockStep"), "blocked.cu"
-        else:
-            config.validate()
-            if not np.isclose(config.omega, 1 / tau, rtol=1e-14, atol=0):
-                raise RuntimeError("level omega/tau mismatch")
-        host = self._boundary_hosts[index]
-        n_objects = len(host.objects) if name == "PatchFineStep" and host is not None else 0
-        block_root_bc = 0
-        if index == 0 and layout is not None:
-            block_root_bc = 2 if self._channel else 1 if self._wall_aware else 0
-        module = compiler.build_module(config, source, arch=self.ctx.sm_arch, n_objects=n_objects,
+        config = level_config(config, layout, tau)
+        name = "BlockChannelStep" if index == 0 and self._channel else "BlockStep"
+        block_root_bc = (2 if self._channel else 1 if self._wall_aware else 0) if index == 0 else 0
+        module = compiler.build_module(config, "blocked.cu", arch=self.ctx.sm_arch,
                                        block_root_bc=block_root_bc)
-        compiler.set_params(module, config.omega, n_objects=n_objects)
+        compiler.set_params(module, config.omega)
         shape = (9, *self.hierarchy.kinds[index].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)
+        size = int(np.prod(layout.blocked_shape)) * 4
+        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]
-        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._channel:
-                cfg = self._remap_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
-                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))
+        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
+            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
+            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))
         return level
 
     def _build_pair(self, index, pair):
         """Upload the transfer tables of levels (index, index+1) and the child's block kinds."""
         parent, child = self._levels[index], self._levels[index + 1]
         order = np.argsort(pair.coalesce_coarse)
-        fine = None if child.layout is None else blocked_fine_tables(pair, child.layout, order)
-        base = None
-        if parent.layout is not None:
-            if index == 0:
-                base = blocked_base_tables(pair, parent.layout, order, self._wall_aware)
-            else:
-                # Nested pairs keep inherited INNER/OUTER and sparse holes of the parent.
-                base = blocked_base_tables(pair, parent.layout, order,
-                                           composed_kind=self.hierarchy.kinds[index])
-        ghosts = pair.ghost_fine if fine is None else fine[0]
-        donors = pair.ghost_coarse if base is None else base[0]
-        if fine is not None:
-            donors = donors.astype(np.uint64)
-        targets = pair.coalesce_coarse[order] if base is None else base[1]
-        children = pair.coalesce_fine[order] if fine is None else fine[1]
-        result = _Pair(child.module.get_function("PatchExplode" if fine is None else "BlockExplode"),
-                       None if parent.layout is None else parent.module.get_function("BlockCoalesce"),
-                       self._upload(ghosts), self._upload(donors), self._upload(targets),
-                       self._upload(children), len(ghosts), len(targets),
-                       int(np.prod(pair.coarse_shape if parent.layout is None
-                                   else parent.layout.local_shape)))
+        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,
+                                                               self._wall_aware)
+        else:
+            # Nested pairs keep inherited INNER/OUTER and sparse holes of the parent.
+            donors, targets, parent_kind = blocked_base_tables(
+                pair, parent.layout, order, composed_kind=self.hierarchy.kinds[index])
+        result = _Pair(child.module.get_function("BlockExplode"),
+                       parent.module.get_function("BlockCoalesce"),
+                       self._upload(ghosts), self._upload(donors.astype(np.uint64)),
+                       self._upload(targets), self._upload(children), len(ghosts), len(targets),
+                       int(np.prod(parent.layout.local_shape)))
         if self._explosion == "projected_linear":
             neighbors, offsets = projected_explosion_tables(
                 pair, parent.layout, self.hierarchy.kinds[index])
             result.explode = child.module.get_function("BlockExplodeProjected")
             result.neighbors, result.offsets = self._upload(neighbors), self._upload(offsets)
-        if child.layout is not None:
-            child.kind = self._upload(blocked_kind(child.layout, self.hierarchy.kinds[index + 1]))
-        if index == 0 and base is not None:
-            parent.kind = self._upload(base[2])
+        child.kind = self._upload(blocked_kind(child.layout, self.hierarchy.kinds[index + 1]))
+        if index == 0:
+            parent.kind = self._upload(parent_kind)
         return result
 
     def _allocate(self, size):
@@ -506,8 +457,6 @@ class FixedPatchStepper:
     def _check_storage(self):
         if any(bool(c.ddf_shifting) != self._ddf_shifting for c in self._level_cfgs):
             raise ValueError("mixed population representations are unsupported")
-        if self._ddf_shifting and any(layout is None for layout in self._layouts):
-            raise ValueError("shifted populations require all-blocked levels")
 
     def initialize(self, *level_g):
         """Upload physical postcollision G, encoding before any FP32 rounding.
@@ -552,16 +501,13 @@ class FixedPatchStepper:
                         or (values.sum(axis=0, dtype=np.float64) + int(self._ddf_shifting) <= 0).any()):
                     raise ValueError("Leaf populations must be finite with positive density")
             self._mask_nonstate(*arrays)
-            arrays = [a if level.layout is None else level.layout.to_blocked(a)
-                      for a, level in zip(arrays, self._levels)]
+            arrays = [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)
                     # A single level starts its spare NaN so stale halos cannot hide.
                     cuda.memcpy_htod(level.spare, np.full_like(array, np.nan)
                                      if len(self._levels) == 1 else array)
-                if self._fine_curved is not None and self._fine_curved.link_exchange_gpu is not None:
-                    cuda.memset_d8(self._fine_curved.link_exchange_gpu, 0, self._fine_curved.count * 12)
                 self._step_count = 0
                 self._exchange_valid = False
                 self._window_valid = False
@@ -660,26 +606,6 @@ class FixedPatchStepper:
                 raise
             return candidate
 
-    def _apply_fine_boundary(self, buffer, **kw):
-        """Stage solid source slots; exchange buffer holds this fine substep only."""
-        cl = self._fine_curved
-        if cl is None:
-            if self._boundary_hosts[1] is not None:
-                raise RuntimeError("Fine boundary device buffers are missing")
-            return
-        if self._curved is None or cl.fluid_idx_gpu is None:
-            raise RuntimeError("Fine boundary device buffers are missing")
-        driven = bool(self._level_jets[1]) or any(
-            obj.surface_actuation is not None for obj in self._boundary_hosts[1].objects)
-        if driven and (cl.surface_velocity_gpu is None or cl.link_exchange_gpu is None):
-            raise RuntimeError("Surface boundary device buffers are missing")
-        self._curved(buffer, cl.fluid_idx_gpu, cl.dir_gpu, cl.q_gpu,
-                     cl.rx_gpu, cl.ry_gpu, cl.rz_gpu, cl.fallback_class_gpu, cl.body_id_gpu,
-                     np.uintp(0), np.uintp(0), np.uint32(cl.count),
-                     cl.surface_velocity_gpu if cl.surface_velocity_gpu is not None else np.uintp(0),
-                     cl.link_exchange_gpu if cl.link_exchange_gpu is not None else np.uintp(0),
-                     grid=((cl.count + self._block[0] - 1) // self._block[0], 1), **kw)
-
     def set_jet_amplitude(self, index, amplitude, *, level=1):
         """Replace one jet's amplitude between steps; no geometry/table rebuild.
 
@@ -703,9 +629,7 @@ class FixedPatchStepper:
                 driven = jet_id >= 0
                 velocity[driven] = jet_link_velocity(jets, shape, jet_id)[driven]
             try:
-                if level == 1 and self._fine_curved is not None:
-                    self._fine_curved.set_jet_amplitude(index, amplitude, self._stream)
-                elif velocity is not None:
+                if velocity is not None:
                     self._stream.synchronize()  # no in-flight step reads the column
                     cuda.memcpy_htod(boundary.device["surface_velocity"], velocity)
             except BaseException:
@@ -747,19 +671,6 @@ class FixedPatchStepper:
         nested = index + 1 < len(self._levels)
         if nested:
             self._advance(index + 1, src, **kw)
-        if level.layout is None:
-            if index == 0:  # dense root: PatchCoarseStep fuses restriction and collision
-                pair = self._pairs[0]
-                level.kernel(src, dst, self._levels[1].current, pair.targets, pair.children,
-                             np.uint32(pair.n_rows), *(np.uint32(v) for v in self.hierarchy.pairs[0].bounds),
-                             np.uint32(self._wall_aware), grid=level.grid, **kw)
-                return
-            if self._boundary_hosts[index] is not None and self._fine_flags is None:
-                raise RuntimeError("Fine boundary flags are missing")
-            self._apply_fine_boundary(src, **kw)
-            flags = self._fine_flags if self._fine_flags is not None else np.uintp(0)
-            level.kernel(src, dst, flags, grid=level.grid, **kw)
-            return
         # Child donors were read above; no halo gather may overwrite staged arrivals below.
         self._gather_halo(level, src, **kw)
         if nested:
@@ -924,7 +835,7 @@ class FixedPatchStepper:
         with self._lock, self.ctx.activate():
             self._require("ready")
             self._check_storage()
-            arrays = [np.empty(level.shape if level.layout is None else level.layout.blocked_shape,
+            arrays = [np.empty(level.layout.blocked_shape,
                                np.float32) for level in self._levels]
             try:
                 for level, array in zip(self._levels, arrays):
@@ -932,7 +843,7 @@ class FixedPatchStepper:
             except BaseException:
                 self._status = "failed"
                 raise
-            arrays = [a if level.layout is None else level.layout.to_dense(a)
+            arrays = [level.layout.to_dense(a)
                       for a, level in zip(arrays, self._levels)]
             self._mask_nonstate(*arrays)
             return _StoredPopulations(tuple(arrays), self._ddf_shifting)
@@ -958,14 +869,6 @@ class FixedPatchStepper:
                         b.window_integral = None
                 self._levels.clear()
                 self._pairs = []
-                if getattr(self, "_fine_curved", None) is not None:
-                    try:
-                        self._fine_curved.free()
-                    except Exception as exc:
-                        error = error or exc
-                    self._fine_curved = None
-                self._curved = self._fine_flags = None
-                self._link_tables = {}
                 for pointer in self._allocations:
                     try:
                         pointer.free()
@@ -1007,7 +910,3 @@ class FixedPatchStepper:
     _fine_leaf = property(lambda self: self._leaf[1])
     _child_leaf = property(lambda self: self._leaf[2] if len(self._leaf) > 2 else None)
     _middle_kind = property(lambda self: self.hierarchy.kinds[1] if self.hierarchy.depth > 2 else None)
-
-
-# Preserve the original entry point and its periodic default.
-PeriodicPatchStepper = FixedPatchStepper

tests/unit/test_example_tiles.py

@@ -46,6 +46,11 @@ class ExampleTiles(unittest.TestCase):
             FixedPatchStepper(cfg, context=forbidden(), coarse_boundary="slot_reservoir",
                               refinements=[bounds], level_blocks=[tuple(b) for b in blocks])
 
+    def test_fixed_periodic_patch_tiles(self):
+        with self.assertRaisesRegex(AssertionError, "CUDA reached"):
+            FixedPatchStepper(LBMConfig(nx=16, ny=16, viscosity=.1, velocity=1e-3), context=forbidden(),
+                              refinements=[(4, 4, 12, 12)], level_blocks=[(16, 16), (24, 24)])
+
     def test_fixed_patch_circle_default_tiles(self):
         bounds = literal("fixed_patch_circle.py", "BOUNDS")
         nx, ny = literal("fixed_patch_circle.py", "NX"), literal("fixed_patch_circle.py", "NY")

tests/unit/test_projected_pair_stepper.py

@@ -146,7 +146,6 @@ class ProjectedPairDispatch(unittest.TestCase):
         for cfg, changes in [(CFG, {'explosion':'unknown'}),
                              (CFG, {'coarse_boundary':'slot_reservoir'}),
                              (CFG, {'coarse_boundary':'channel'}),
-                             (CFG, {'level_blocks':[None,None], 'refinements':[(6,6,26,26)]}),
                              (CFG, {'level_objects':[(),(object(),)]}),
                              (CFG, {'level_jets':[(),(object(),)]}),
                              (replace(CFG,collision='MRT'), {}),
@@ -155,9 +154,8 @@ class ProjectedPairDispatch(unittest.TestCase):
                 host.FixedPatchStepper(cfg, **dict(options, **changes))
         context.activate.assert_not_called()
         host.cuda.mem_alloc.assert_not_called()
-        # The legacy constant path still allows MRT and dense storage.
-        for options in ({'cfg':replace(CFG,collision='MRT')},
-                        {'level_blocks':[None,None], 'refinements':[(6,6,26,26)]}):
+        # The constant path still allows MRT.
+        for options in ({'cfg':replace(CFG,collision='MRT')},):
             run, _, _, _, _ = backend(**options)
             run.close()
 

tests/unit/test_raw_exchange.py

@@ -153,15 +153,6 @@ class RawExchange(unittest.TestCase):
         run._stream.synchronize.side_effect = None
         run.close()
 
-    def test_dense_observation_rejected_before_device_access(self):
-        host = host_module()
-        ctx = Mock()
-        with self.assertRaisesRegex(ValueError, 'blocked body'):
-            host.FixedPatchStepper(LBMConfig(nx=32, ny=32, viscosity=.1), context=ctx,
-                                   refinements=[(4, 4, 28, 28)],
-                                   level_objects=[(), (body(radius=2.15),)], observe_raw_exchange=True)
-        ctx.activate.assert_not_called()
-
     def test_all_depth_weights(self):
         for level in (1, 2, 3, 7):
             self._check_weights_split_and_window(level)

tests/unit/test_refinement_mrt.py

@@ -27,23 +27,20 @@ class RefinementMRTHost(unittest.TestCase):
         block_collision = (KERNELS / "operators/block_collision.cuh").read_text()
         self.assertIn('#include "operators/block_collision.cuh"', blocks)
         self.assertEqual(blocks.count("block_collide(f);"), 1)
-        refinement = (KERNELS / "refinement.cu").read_text() + (KERNELS / "step/patch_step.cu").read_text()
         for text in ("compute_feq_from_momentum(drho, jx, jy, ux, uy, feq);",
                      "compute_feq(rho, ux, uy, feq);",
                      "collide_srt_no_force(f, feq, d_params.omega);"):
             self.assertIn(text, block_collision)
         self.assertIn("collide_srt(f, feq, forcing, omega);", channel)
-        self.assertIn("collide_srt_no_force(f, feq, d_params.omega);", refinement)
-        for source in (blocks, refinement):
+        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)
-        self.assertIn("USE_DDF_SHIFTING ||", refinement)  # dense/mixed shifted stays rejected
         # 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, refinement):
+        for source in (blocks + "\n" + block_collision, channel):
             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)

tests/unit/test_shifted_storage.py

@@ -125,11 +125,11 @@ class ShiftedStorage(unittest.TestCase):
             with self.assertRaisesRegex(ValueError,'mixed'): run._download_raw()
         finally: run.close()
 
-    def test_dense_or_mixed_shifted_rejected_before_compilation(self):
+    def test_missing_tiles_rejected_before_compilation(self):
         module=backend.host_module()
         cfg=replace(backend.CFG,ddf_shifting=True)
         for blocks in ([None,None],[(4,4),None],[None,(4,4)]):
-            with self.assertRaisesRegex(ValueError,'all-blocked|blocked refined levels below'):
+            with self.assertRaisesRegex(ValueError,'must be a tile'):
                 module.FixedPatchStepper(cfg,context=None,refinements=[(4,4,20,28)],
                                          level_blocks=blocks)
         module.compiler.build_module.assert_not_called()

tests/validation/multigrid/test_exchange_window.py

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

tests/validation/multigrid/test_fine_surface_gpu.py

@@ -1,13 +1,17 @@
-"""Two coarse/four fine steps: isolated fine circle, no flux-accuracy claim.
+"""Two coarse/four fine steps: isolated fine circle on blocked storage, no flux-accuracy claim.
 
-Run only through the bounded GPU runner after a fresh resource check. The test
-hooks read actual kernel inputs/outputs; they never change the stepping order.
+Run only through the bounded GPU runner after a fresh resource check. The fine
+step hook reads the actual output buffer; it never changes the stepping order.
+
+Retired with the dense refinement route (accepted losses, see the D2 commit):
+the dense-kernel per-hop flux ledger (fine/coarse transfer + boundary budget),
+the CurvedLinkSoA reflection oracle and the CELERIS_REFINEMENT_AUDIT collision
+counts. The blocked test keeps the four-hop comparison with the uniform solver,
+the raw-exchange mass sign and passive == zero actuation; raw exchange covers
+the boundary part only, not full population flux closure.
 """
 import importlib.util
-from contextlib import nullcontext
 from pathlib import Path
-from threading import RLock
-from types import SimpleNamespace
 from unittest.mock import Mock, patch
 
 import numpy as np
@@ -18,7 +22,6 @@ import CelerisLab
 from CelerisLab.body.geometry.circle import CircleGeometry
 from CelerisLab.body.objects import SimObject
 from CelerisLab.config import LBMConfig
-from CelerisLab.cuda import compiler
 from CelerisLab.cuda.context import CudaContext
 from CelerisLab.lbm import refinement_stepper as implementation
 from CelerisLab.lbm.descriptors import FLUID, SOLID
@@ -27,10 +30,8 @@ from CelerisLab.lbm.descriptors import FLUID, SOLID
 ROOT = Path(__file__).resolve().parents[3]
 CX = np.array([0, 1, -1, 0, 0, 1, -1, 1, -1])
 CY = np.array([0, 0, 0, 1, -1, 1, -1, -1, 1])
-OPP = np.array([0, 2, 1, 4, 3, 6, 5, 8, 7])
 W = np.array([4/9]+[1/9]*4+[1/36]*4)
 MOMENTS = np.stack((np.ones(9), CX, CY))
-GEOMETRY = ("fluid_idx", "dir", "q", "body_id", "rx", "ry", "rz", "fallback_class")
 BOUNDS = (4, 4, 20, 20)
 
 
@@ -53,197 +54,6 @@ def object_spec(speed):
                      radius=2.15, surface_actuation=surface(speed))
 
 
-def audit_compiler(original):
-    def compile_private(*args, **kwargs):
-        entry = Path(kwargs["build_dir"])/"fine_surface_audit.cu"
-        entry.write_text('#define CELERIS_REFINEMENT_AUDIT\n#include "refinement.cu"\n')
-        kwargs["source"] = entry.name
-        return original(*args, **kwargs)
-    return compile_private
-
-
-def reflection_oracle(before, cl, speed, *, verify_projection=True):
-    """FP64 BFL from frozen spatial populations, independent arc projection."""
-    flat = before.astype(float).reshape(9, -1)
-    radius = np.hypot(cl.rx.astype(float), cl.ry.astype(float))
-    delta = (np.arctan2(cl.ry.astype(float), cl.rx.astype(float))-.7+np.pi) % (2*np.pi)-np.pi
-    amplitude = 0 if speed is None else speed
-    strength = amplitude*np.maximum(0, 1-(delta/.85)**2)**2
-    velocity = strength[:, None]*np.column_stack((cl.rx, cl.ry))/radius[:, None]
-    if verify_projection and cl.surface_velocity is not None:
-        np.testing.assert_allclose(cl.surface_velocity, velocity, atol=3e-9, rtol=0)
-    incoming = []
-    for i, (k, d, q, fallback) in enumerate(zip(cl.fluid_idx, cl.dir, cl.q, cl.fallback_class)):
-        k, d, q = int(k), int(d), float(q)
-        rho = flat[:, k].sum()
-        dot = CX[d]*velocity[i, 0]+CY[d]*velocity[i, 1]
-        if fallback:
-            value = flat[d, k]-6*W[d]*rho*dot
-        elif q < .5:
-            value = 2*q*flat[d, k]+(1-2*q)*flat[d, k-CX[d]-40*CY[d]]-6*W[d]*rho*dot
-        else:
-            value = flat[d, k]/(2*q)+(1-1/(2*q))*flat[OPP[d], k]-3*W[d]*rho*dot/q
-        incoming.append(value)
-    return np.array(incoming)
-
-
-def fine_flux(before, staged, active, solid):
-    """Population ledger across the fine leaf boundary: solid versus halo.
-
-    Internal fluid-to-fluid streaming cancels pairwise; solid and halo packets
-    are counted in their actual direction, before any collision.
-    """
-    boundary, transfer = np.zeros(9), np.zeros(9)
-    magnitude = 0.0
-    for d, (cx, cy) in enumerate(zip(CX, CY)):
-        incoming = active & ~np.roll(active, (cy, cx), (0, 1))
-        outgoing = active & ~np.roll(active, (-cy, -cx), (0, 1))
-        from_solid = np.roll(solid, (cy, cx), (0, 1))
-        into_solid = np.roll(solid, (-cy, -cx), (0, 1))
-        received = np.roll(staged[d], (cy, cx), (0, 1)).astype(float)
-        sent = before[d].astype(float)
-        for vector, inc, out in ((boundary, incoming & from_solid, outgoing & into_solid),
-                                 (transfer, incoming & ~from_solid, outgoing & ~into_solid)):
-            assert np.isfinite(received[inc]).all() and np.isfinite(sent[out]).all()
-            vector[d] = received[inc].sum()-sent[out].sum()
-            magnitude += np.abs(received[inc]).sum()+np.abs(sent[out]).sum()
-    return boundary, transfer, magnitude
-
-
-def coarse_flux(before, fine, active):
-    """Derive every replacement from physical midpoints, not production tables."""
-    x0, y0, x1, y1 = BOUNDS
-    incoming = np.stack([np.roll(before[d], (CY[d], CX[d]), (0, 1)) for d in range(9)]).astype(float)
-    magnitude = 0.0
-    for y, x in np.argwhere(active):
-        for d in range(1, 9):
-            if x0 <= x+.5-CX[d]/2 <= x1 and y0 <= y+.5-CY[d]/2 <= y1:
-                a, b = 4+2*(x-x0), 4+2*(y-y0)
-                children = fine[d, b:b+2, a:a+2].astype(float)
-                assert children.shape == (2, 2) and np.isfinite(children).all()
-                incoming[d, y, x] = children.sum()/4
-                magnitude += np.abs(children).sum()/4
-    assert np.isfinite(incoming[:, active]).all()
-    # This retains the eight grazing-corner replacements of uncovered donors.
-    flux = (incoming[:, active]-before[:, active].astype(float)).sum(axis=1)
-    magnitude += np.abs(incoming[:, active]).sum()+np.abs(before[:, active]).sum(dtype=float)
-    return flux, magnitude
-
-
-def run_patch(speed, *, audit=True):
-    cfg = LBMConfig(nx=24, ny=24, viscosity=.05, velocity=0, rho=1)
-    coarse_active = np.ones((24, 24), bool)
-    coarse_active[4:20, 4:20] = False
-    source = object_spec(speed)
-    state = vars(source.state).copy()
-    with CudaContext(0) as context:
-        build = (patch.object(compiler, "compile_kernel", audit_compiler(compiler.compile_kernel))
-                 if audit else nullcontext())
-        with build:
-            backend = implementation.FixedPatchStepper(cfg, context=context, refinements=[BOUNDS],
-                                             level_objects=[(), (source,)])
-        with backend, context.activate():
-            c, f = backend._levels
-            cl = backend._fine_curved
-            geometry = {name: getattr(cl, name).copy() for name in GEOMETRY}
-            flags = backend.fine_boundary.flags.reshape(40, 40)
-            solid = (flags & SOLID) != 0
-            active = np.zeros((40, 40), bool)
-            active[4:-4, 4:-4] = True
-            active &= (flags & FLUID) != 0
-            active &= ~solid
-            local = backend.fine_boundary.objects[0]
-            assert local.center == (20., 20.) and local.radius == 4.3
-            assert source.center == (12.25, 12.25) and source.radius == 2.15
-            yy, xx = np.mgrid[:40, :40]
-            roi = active & ((xx-20)**2+(yy-20)**2 <= 9**2)
-            backend.initialize(rest((24, 24)), rest((40, 40)))
-            if audit:
-                for level in backend._levels:
-                    for name in ("collision_counts", "explosion_counts", "incoming_counts"):
-                        pointer, size = level.module.get_global(name)
-                        cuda.memset_d8(pointer, 0, size)
-            initial = backend.download()
-            initial_pop = initial[0][:, coarse_active].sum(axis=1, dtype=float)
-            initial_pop += .25*initial[1][:, active].sum(axis=1, dtype=float)
-            boundary, fine_transfer, coarse_transfer = np.zeros(9), np.zeros(9), np.zeros(9)
-            budget_sum, reflection_error, drive_mass = 0., 0., 0.
-            traces = []
-            original_curved, original_fine, original_coarse = backend._curved, f.kernel, c.kernel
-
-            def curved(*args, **kwargs):
-                nonlocal boundary, fine_transfer, budget_sum, reflection_error, drive_mass
-                backend._stream.synchronize()
-                before = read(args[0], f.shape)
-                expected = reflection_oracle(before, cl, speed)
-                original_curved(*args, **kwargs)
-                backend._stream.synchronize()
-                staged = read(args[0], f.shape)
-                d, k = cl.dir.astype(int), cl.fluid_idx.astype(int)
-                values = staged.reshape(9, -1)[OPP[d], k+CX[d]+40*CY[d]]
-                passive = reflection_oracle(before, cl, None, verify_projection=False)
-                drive_mass += .25*float(np.sum(values-passive))
-                error = float(np.max(np.abs(values-expected)))
-                reflection_error = max(reflection_error, error)
-                assert error <= 1e-7
-                np.testing.assert_array_equal(staged[:, active], before[:, active])
-                b, transfer, magnitude = fine_flux(before, staged, active, solid)
-                boundary += .25*b
-                fine_transfer += .25*transfer
-                budget_sum += .25*magnitude
-
-            def fine(*args, **kwargs):
-                nonlocal budget_sum
-                backend._stream.synchronize()
-                staged = read(args[0], f.shape)
-                original_fine(*args, **kwargs)
-                backend._stream.synchronize()
-                advanced = read(args[1], f.shape)
-                assert np.isfinite(advanced[:, active]).all()
-                assert np.isnan(advanced[:, solid]).all()
-                # Include actual collision operands in the FP32 rounding budget.
-                budget_sum += .25*np.abs(staged[:, active]).sum(dtype=float)
-                traces.append(advanced[:, roi].copy())
-
-            def coarse(*args, **kwargs):
-                nonlocal coarse_transfer, budget_sum
-                backend._stream.synchronize()
-                before, fine_values = read(args[0], c.shape), read(args[2], f.shape)
-                transfer, magnitude = coarse_flux(before, fine_values, coarse_active)
-                coarse_transfer += transfer
-                budget_sum += magnitude
-                return original_coarse(*args, **kwargs)
-
-            with patch.object(backend, "_curved", curved), patch.object(f, "kernel", fine), \
-                    patch.object(c, "kernel", coarse):
-                backend.step(2)
-            actual = backend.download()
-            population = actual[0][:, coarse_active].sum(axis=1, dtype=float)
-            population += .25*actual[1][:, active].sum(axis=1, dtype=float)
-            ledger = boundary+fine_transfer+coarse_transfer
-            error = MOMENTS@(population-initial_pop-ledger)
-            budget = 8*np.finfo(np.float32).eps*budget_sum
-            assert np.isfinite(budget_sum) and budget_sum > 0
-            assert np.max(np.abs(error)) <= budget
-            if speed:
-                assert speed*drive_mass > 0
-            assert len(traces) == 4 and backend.step_count == 2
-            if audit:
-                for level, mask, count in ((c, coarse_active, 2), (f, active, 4)):
-                    pointer, size = level.module.get_global("collision_counts")
-                    counts = read(pointer, mask.shape, np.uint32)
-                    assert size == counts.nbytes
-                    np.testing.assert_array_equal(counts, count*mask)
-            for name, values in geometry.items():
-                np.testing.assert_array_equal(getattr(cl, name), values)
-            assert vars(source.state) == state and source.geometry.center == source.center
-            print(f"speed={speed} audit={audit} BFL={reflection_error:.9g} ledger={error.tolist()} "
-                  f"budget={budget:.9g} boundary={list(MOMENTS@boundary)} "
-                  f"transfer={list(MOMENTS@(fine_transfer+coarse_transfer))} "
-                  f"drive_mass={drive_mass:+.9g}; counts={'PASS' if audit else 'normal build'}", flush=True)
-            return traces, roi, actual
-
-
 def uniform_reference(speed, roi):
     # Reuse the existing explicit config helper without importing unrelated tests.
     spec = importlib.util.spec_from_file_location("fine_surface_cases", ROOT/"tests/integration/_cases.py")
@@ -261,62 +71,59 @@ def uniform_reference(speed, roi):
         return traces
 
 
+def run_blocked_patch(speed):
+    """All-blocked run (whole-level tiles) of the same case; fine ROI after each fine hop."""
+    cfg = LBMConfig(nx=24, ny=24, viscosity=.05, velocity=0, rho=1)
+    source = object_spec(speed)
+    with CudaContext(0) as context:
+        with implementation.FixedPatchStepper(cfg, context=context, refinements=[BOUNDS],
+                                              level_blocks=[(24, 24), (40, 40)],
+                                              level_objects=[(), (source,)],
+                                              observe_raw_exchange=True) as backend:
+            c, f = backend._levels
+            flags = backend.fine_boundary.flags.reshape(40, 40)
+            active = np.zeros((40, 40), bool)
+            active[4:-4, 4:-4] = True
+            active &= ((flags & FLUID) != 0) & ((flags & SOLID) == 0)
+            yy, xx = np.mgrid[:40, :40]
+            roi = active & ((xx-20)**2+(yy-20)**2 <= 9**2)
+            backend.initialize(rest((24, 24)), rest((40, 40)))
+            traces, original = [], f.kernel
+
+            def fine(*args, **kwargs):
+                original(*args, **kwargs)
+                backend._stream.synchronize()
+                advanced = f.layout.to_dense(read(args[1], f.layout.blocked_shape))
+                assert np.isfinite(advanced[:, active]).all()
+                assert np.isnan(advanced[:, (flags & SOLID) != 0]).all()
+                traces.append(advanced[:, roi].copy())
+
+            with patch.object(f, "kernel", fine):
+                backend.step(2)
+            assert len(traces) == 4 and backend.step_count == 2
+            mass = backend.read_raw_exchange(0)["total"]["mass_into_fluid"]
+            return traces, roi, backend.download(), float(mass)
+
+
 @pytest.mark.parametrize("speed", [0., .002, -.002])
-def test_fine_surface_four_hops(speed):
-    assert Path(CelerisLab.__file__).resolve() == ROOT/"src/CelerisLab/__init__.py"
-    assert Path(compiler._KERNEL_DIR).resolve() == ROOT/"src/CelerisLab/lbm/kernels"
-    assert Path(implementation.__file__).resolve() == ROOT/"src/CelerisLab/lbm/refinement_stepper.py"
-    actual, roi, final = run_patch(speed)
+def test_fine_surface_blocked_matches_uniform_four_hops(speed):
+    """Blocked storage: fine ROI after each of four hops against the uniform 48^2 solver
+    (same tolerances as the dense test); blowing/suction drives mass with its sign;
+    zero actuation equals the passive circle bit for bit."""
+    actual, roi, final, mass = run_blocked_patch(speed)
     reference = uniform_reference(speed, roi)
     for got, wanted in zip(actual, reference):
         np.testing.assert_allclose(got, wanted, atol=3e-6, rtol=0)
         np.testing.assert_allclose(MOMENTS@got.astype(float), MOMENTS@wanted.astype(float),
                                    atol=2e-5, rtol=0)
-    if speed == 0:
-        passive, passive_roi, passive_final = run_patch(None, audit=False)
+    if speed:
+        assert speed*mass > 0
+    else:
+        passive, passive_roi, passive_final, _ = run_blocked_patch(None)
         np.testing.assert_array_equal(passive_roi, roi)
         for got, wanted in zip(actual+list(final), passive+list(passive_final)):
             np.testing.assert_array_equal(got, wanted)
-
-
-@pytest.mark.parametrize("missing", ["surface_velocity_gpu", "link_exchange_gpu"])
-def test_missing_surface_device_buffer_rejects_before_launch(missing):
-    backend = implementation.FixedPatchStepper.__new__(implementation.FixedPatchStepper)
-    backend._boundary_hosts = [None, SimpleNamespace(objects=[object_spec(.002)])]
-    backend._level_jets = [(), ()]
-    backend._fine_curved = SimpleNamespace(fluid_idx_gpu=object(),
-                                          surface_velocity_gpu=object(), link_exchange_gpu=object())
-    setattr(backend._fine_curved, missing, None)
-    backend._curved = Mock()
-    with pytest.raises(RuntimeError, match="Surface boundary device buffers are missing"):
-        backend._apply_fine_boundary(object())
-    backend._curved.assert_not_called()
-
-
-def test_close_frees_fine_boundary_after_stream_failure_once():
-    backend = implementation.FixedPatchStepper.__new__(implementation.FixedPatchStepper)
-    backend.ctx = SimpleNamespace(activate=nullcontext)
-    backend._lock = RLock()
-    backend._status = "ready"
-    stream, soa, first, second = Mock(), Mock(), Mock(), Mock()
-    stream.synchronize.side_effect = RuntimeError("stream synchronization failed")
-    backend._stream = stream
-    backend._fine_curved = soa
-    backend._fine_flags = first
-    backend._curved = Mock()
-    backend._levels = [object(), object()]
-    backend._allocations = [first, second]
-    with pytest.raises(RuntimeError, match="stream synchronization failed"):
-        backend.close()
-    assert backend._status == "closed"
-    assert backend._levels == [] and backend._allocations == []
-    assert backend._fine_curved is None and backend._fine_flags is None
-    assert backend._curved is None and backend._stream is None
-    backend.close()
-    stream.synchronize.assert_called_once_with()
-    soa.free.assert_called_once_with()
-    first.free.assert_called_once_with()
-    second.free.assert_called_once_with()
+    print(f"blocked speed={speed} four hops vs uniform PASS mass_into_fluid={mass:+.9g}", flush=True)
 
 
 def test_cross_interface_donor_rejects_before_cuda():
@@ -328,6 +135,7 @@ def test_cross_interface_donor_rejects_before_cuda():
     with patch.object(cuda, "mem_alloc") as allocate:
         with pytest.raises(ValueError, match="donor crosses"):
             implementation.FixedPatchStepper(cfg, context=context, refinements=[BOUNDS],
+                                             level_blocks=[(24, 24), (40, 40)],
                                              level_objects=[(), (source,)])
         context.activate.assert_not_called()
         allocate.assert_not_called()

tests/validation/multigrid/test_raw_exchange_gpu.py

@@ -92,8 +92,8 @@ def census(run, index, buffer):
 
 
 def fine(block):
-    """Two-level circle+jet construction on BOUNDS; block None keeps the dense fine level."""
-    return dict(refinements=[BOUNDS], level_blocks=[None, block],
+    """Two-level circle+jet construction on BOUNDS, blocked with a (12, 8) root tile."""
+    return dict(refinements=[BOUNDS], level_blocks=[(12, 8), block],
                 level_objects=[(), (circle(),)], level_jets=[(), (jet(),)])
 
 class RawExchangeGPU(unittest.TestCase):
@@ -146,26 +146,6 @@ class RawExchangeGPU(unittest.TestCase):
                         np.testing.assert_array_equal(measured.read_raw_exchange(0)['total']['impulse'],
                                                       obs['total']['impulse'])
 
-    def test_dense_fine_storage_reference(self):
-        """Same fine resolution/initial data with dense storage: a numerical
-        reference for cut-link rows, not a uniform-whole-domain convergence study."""
-        with CudaContext(0) as ctx:
-            with FixedPatchStepper(CFG, context=ctx, **fine(None)) as dense, \
-                 FixedPatchStepper(CFG, context=ctx, **fine((16, 20)), observe_raw_exchange=True) as blocked:
-                arrays = [state(l.shape) for l in blocked._levels]
-                dense.initialize(*arrays)
-                blocked.initialize(*arrays)
-                for speed in (0., .002, 0.):
-                    dense.set_jet_amplitude(0, speed)
-                    blocked.set_jet_amplitude(0, speed)
-                    dense.step(); blocked.step()
-                    for a, e in zip(blocked.download(), dense.download()):
-                        np.testing.assert_array_equal(a, e)
-                    count = dense._fine_curved.count
-                    expected = download(dense._fine_curved.link_exchange_gpu, (count, 3))
-                    actual = download(blocked._level_boundaries[1].last_raw, (count, 4))
-                    np.testing.assert_array_equal(actual[:, :3], expected)
-
     def test_depth_four_interval_storage_and_sampling(self):
         cfg = LBMConfig(nx=48, ny=48, viscosity=.02)
         obj = SimObject(0, geometry=CircleGeometry(16., 16., .5), center=(16., 16.), radius=.5)

tests/validation/multigrid/test_refinement_blocked_base.py

@@ -3,33 +3,11 @@
 Run CPU and GPU classes in separate processes. GPU callers must first check
 resources/device availability; this file never selects or reserves a GPU.
 
-Dense/blocked route parity (DenseBlockedRouteParity): the same configuration,
-initial state, body/Jet schedule and tilings on three storage routes, only
-level_blocks varying: D/D [None, None], D/B [None, FINE_BLOCK] (dense root,
-the mixed pattern of current examples) and B/B [(12, 8), FINE_BLOCK]. Gate:
-every public and owner array identical in FP32 bits and NaN masks at steps
-0, 1, 2, 8 and after re-initialize; D/B and B/B raw-exchange ledgers equal.
-No tolerance. A mismatch is reported with route pair/level/step and diagnosed
-(boundary semantics vs arithmetic ordering), never widened.
-
-Route capabilities (rejections checked before CUDA by test_route_capability_matrix):
-
-  combination                                   D/D       D/B       B/B
-  periodic / slot_reservoir root, fine circle   yes       yes       yes   (parity gate)
-  actuated circle, velocity Jet + amplitude     yes       yes       yes   (parity gate)
-  raw/window exchange on a body level           rejected  yes       yes
-  shifted storage                               rejected  rejected  yes
-  projected explosion                           rejected  rejected  yes
-  channel root                                  rejected  rejected  yes
-  several L1 rectangles                         rejected  rejected  yes
-  depth >= 3 with a dense level                 rejected  rejected  yes
-Slot/reservoir walls are compile-time WALL_* macros on a dense root and the
-d_block_wall constant on a blocked root: BC-encoding parity, bit-identical by
-this gate. Dense fine bodies use CurvedLinkSoA + CurvedBoundaryKernel
-(n_objects > 0), blocked ones BlockCurvedBoundary; same links and physics,
-compared through the population rows. Parked here: MRT (covered by
-test_refinement_mrt_gpu::test_dense_versus_blocked_mrt_two_level), shifted and
-depth >= 3 (no dense route to compare), several bodies (later slice).
+Dense refinement storage was removed: every level is blocked and level_blocks
+is required (TileContractCPU). The retired dense/blocked route parity gates
+(BlockedBaseGPU exact parity, DenseBlockedRouteParity incl. MRT and two
+bodies, the route capability matrix) passed and are recorded at 1573942 and
+2469af7.
 """
 from contextlib import ExitStack, contextmanager, nullcontext
 from dataclasses import replace
@@ -48,7 +26,6 @@ from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
 from CelerisLab.lbm.refinement import ACTIVE, OUTER_GHOST, FixedPatch2D
 from CelerisLab.lbm.refinement_stepper import FixedPatchStepper, blocked_base_tables
 from test_refinement_blocked_fine import NX, NY, BOUNDS, config, state, circle, actuation
-from test_refinement_jets import SPEED, jet
 
 BASE_BLOCKS = ((24, 32), (12, 8), (3, 4), (1, 1))
 FINE_BLOCK = (16, 20)
@@ -295,185 +272,33 @@ class BlockedBaseCPU(unittest.TestCase):
                 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)])
+class TileContractCPU(unittest.TestCase):
+    """level_blocks is required: one positive (by, bx) tile per level, checked before CUDA."""
 
-    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, context=context, coarse_boundary=mode,
-                                               refinements=[BOUNDS], level_blocks=[base_block, FINE_BLOCK],
-                                               level_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)
-
-
-ROUTES = (("D/D", [None, None]), ("D/B", [None, FINE_BLOCK]), ("B/B", [(12, 8), FINE_BLOCK]))
-
-
-def same_summary(a, b):
-    """Exact equality of nested exchange summaries (dicts, sequences, arrays, numbers)."""
-    if isinstance(a, dict):
-        return a.keys() == b.keys() and all(same_summary(a[k], b[k]) for k in a)
-    if isinstance(a, (list, tuple)):
-        return len(a) == len(b) and all(same_summary(x, y) for x, y in zip(a, b))
-    a, b = np.asarray(a), np.asarray(b)
-    if a.dtype.kind in "biufc" and b.dtype.kind in "biufc":
-        return a.shape == b.shape and np.array_equal(a, b, equal_nan=True)
-    return a.dtype == b.dtype and np.array_equal(a, b)  # labels such as the window name
-
-
-class DenseBlockedRouteCPU(unittest.TestCase):
-    def test_same_summary_compares_labels_and_every_number(self):
-        ledger = {"window": "last_base_step", "time_start": 3, "total": {"force": [1.0, -2.0],
-                  "torque": 0.5}, "aperture": {"force": np.array([np.nan, 1.0])}}
-        self.assertTrue(same_summary(ledger, replace_value(ledger, ())))
-        for path, value in ((("window",), "last_substep"), (("time_start",), 4),
-                            (("total", "force"), [1.0, -2.000001]), (("total", "torque"), 0.25),
-                            (("aperture", "force"), np.array([0.0, 1.0]))):
-            with self.subTest(path=path):
-                self.assertFalse(same_summary(ledger, replace_value(ledger, path, value)))
-                self.assertFalse(same_summary(replace_value(ledger, path, value), ledger))
-
-    def test_route_capability_matrix(self):
+    def test_level_blocks_required_and_validated(self):
         context = Mock(activate=Mock(side_effect=AssertionError("CUDA reached")))
-        periodic = config("periodic")
-        inner = (BOUNDS[0] + 2, BOUNDS[1] + 2, BOUNDS[2] - 2, BOUNDS[3] - 2)
-        second = (BOUNDS[0], BOUNDS[1], BOUNDS[0] + 2, BOUNDS[1] + 2)
-        cases = (
-            ("raw exchange on dense body level", periodic, dict(
-                level_blocks=[None, None], level_objects=[(), (circle(),)], observe_raw_exchange=True)),
-            ("shifted D/D", replace(periodic, ddf_shifting=True), dict(level_blocks=[None, None])),
-            ("shifted D/B", replace(periodic, ddf_shifting=True), dict(level_blocks=[None, FINE_BLOCK])),
-            ("projected D/D", periodic, dict(level_blocks=[None, None], explosion="projected_linear")),
-            ("projected D/B", periodic, dict(level_blocks=[None, FINE_BLOCK], explosion="projected_linear")),
-            ("channel D/B", periodic, dict(level_blocks=[None, FINE_BLOCK], coarse_boundary="channel")),
-            ("several L1 rectangles", periodic, dict(refinements=[[BOUNDS, second]], level_blocks=[None, None])),
-            ("depth 3 dense", periodic, dict(refinements=[BOUNDS, inner], level_blocks=[None, None, None])),
-        )
-        for label, cfg, kw in cases:
-            kw.setdefault("refinements", [BOUNDS])
+        cfg = config("periodic")
+        for label, kw in (("omitted", {}),
+                          ("None list", dict(level_blocks=None)),
+                          ("too few levels", dict(level_blocks=[(12, 8)])),
+                          ("too many levels", dict(level_blocks=[(12, 8), FINE_BLOCK, FINE_BLOCK])),
+                          ("None root", dict(level_blocks=[None, FINE_BLOCK])),
+                          ("None fine", dict(level_blocks=[(12, 8), None])),
+                          ("wrong arity", dict(level_blocks=[(12,), FINE_BLOCK])),
+                          ("non-integer", dict(level_blocks=[(12, 8.0), FINE_BLOCK])),
+                          ("bool", dict(level_blocks=[(True, 8), FINE_BLOCK])),
+                          ("non-positive", dict(level_blocks=[(0, 8), FINE_BLOCK])),
+                          ("non-tuple", dict(level_blocks=[12, FINE_BLOCK])),
+                          ("not dividing", dict(level_blocks=[(7, 8), FINE_BLOCK]))):
             with self.subTest(label), self.assertRaises(ValueError):
-                FixedPatchStepper(cfg, context=context, **kw)
+                FixedPatchStepper(cfg, context=context, refinements=[BOUNDS], **kw)
+        with self.assertRaisesRegex(ValueError, "whole-level tile"):
+            FixedPatchStepper(cfg, context=context, refinements=[BOUNDS])
         context.activate.assert_not_called()
-
-
-def replace_value(tree, path, value=None):
-    """Deep copy of a nested dict with the item at ``path`` replaced (no change for ())."""
-    if isinstance(tree, dict):
-        return {k: (value if path == (k,) else replace_value(v, path[1:], value) if path[:1] == (k,)
-                    else replace_value(v, (), None)) for k, v in tree.items()}
-    return np.array(tree, copy=True) if isinstance(tree, np.ndarray) else tree
-
-
-class DenseBlockedRouteParity(unittest.TestCase):
-    """D/D vs D/B vs B/B, unshifted SRT, one circle or one velocity Jet; see module docstring."""
-
-    def run_routes(self, mode, objects=(), jets=(), switch=None):
-        from CelerisLab.cuda.context import CudaContext
-        cfg = config(mode)
-        p = FixedPatch2D.from_config(cfg, BOUNDS)
-        arrays = state(NY, NX, 1), state(*p.fine_shape, 2)
-        results = {}
-        with CudaContext(0) as context:
-            for name, blocks in ROUTES:
-                observe = bool(objects) and blocks[1] is not None
-                snapshots, ledgers = [], []
-                with FixedPatchStepper(cfg, context=context, coarse_boundary=mode, refinements=[BOUNDS],
-                                       level_blocks=blocks, level_objects=[(), tuple(objects)],
-                                       level_jets=[(), tuple(jets)],
-                                       observe_raw_exchange=observe) as backend:
-                    backend.initialize(*arrays)
-                    for step in (0, 1, 2, 8):
-                        if step == 8 and switch is not None:
-                            backend.set_jet_amplitude(0, switch)
-                        backend.step(step - backend.step_count)
-                        snapshots.append(backend.download() + raw(backend))
-                        if observe and step:
-                            ledgers.append(tuple(backend.read_raw_exchange(0, window=window)
-                                                 for window in ("last_base_step", "last_substep")))
-                    backend.initialize(*arrays)
-                    snapshots.append(backend.download() + raw(backend))
-                results[name] = snapshots, ledgers
-        return results
-
-    def check(self, label, results):
-        reference, _ = results["D/D"]
-        for name in ("D/B", "B/B"):
-            for step, actual, expected in zip((0, 1, 2, 8, "reset"), results[name][0], reference):
-                for index, (a, b) in enumerate(zip(actual, expected)):
-                    exact(self, a, b, f"{label} {name} vs D/D step={step} array={index}")
-        if results["D/B"][1]:
-            for step, a, b in zip((1, 2, 8), results["D/B"][1], results["B/B"][1]):
-                self.assertTrue(same_summary(a, b), f"{label} raw exchange D/B vs B/B step={step}")
-        for name, (snapshots, _) in results.items():
-            masses = [float(np.nansum(level.astype(np.float64))) for level in snapshots[3][:2]]
-            print(f"{label} {name}: step 8 mass by level {masses}", flush=True)
-        print(f"PASS {label}: D/D = D/B = B/B bits/NaN masks at 0,1,2,8 + reset", flush=True)
-
-    def test_body_free(self):
-        for mode in ("periodic", "slot_reservoir"):
-            self.check(f"{mode} body-free", self.run_routes(mode))
-
-    def test_passive_circle(self):
-        for mode in ("periodic", "slot_reservoir"):
-            self.check(f"{mode} circle", self.run_routes(mode, objects=[circle()]))
-
-    def test_actuated_circle(self):
-        for mode in ("periodic", "slot_reservoir"):
-            for speed in (SPEED, -SPEED):
-                self.check(f"{mode} actuated {speed}",
-                           self.run_routes(mode, objects=[circle(actuation(speed))]))
-
-    def test_velocity_jet_with_amplitude_switch(self):
-        for mode in ("periodic", "slot_reservoir"):
-            self.check(f"{mode} jet switch",
-                       self.run_routes(mode, objects=[circle()], jets=[jet()], switch=-SPEED))
+        # The documented whole-level tiles are always valid (construction reaches CUDA).
+        whole = [(NY, NX), FixedPatch2D.from_config(cfg, BOUNDS).fine_shape]
+        with self.assertRaisesRegex(AssertionError, "CUDA reached"):
+            FixedPatchStepper(cfg, context=context, refinements=[BOUNDS], level_blocks=whole)
 
 
 if __name__ == "__main__":

tests/validation/multigrid/test_refinement_blocked_fine.py

@@ -1,9 +1,8 @@
-"""Opt-in blocked fine storage for FixedPatchStepper: CPU tables plus GPU gate.
+"""Blocked fine storage for FixedPatchStepper: CPU tables, staging and launch contracts.
 
-Frozen GPU gate (same stream, 8 coarse steps, nonequilibrium start), in both
-periodic and slot_reservoir modes: dense fine storage vs 1x1, 2x2 and 4x2
-fine blocks give bit-identical masked downloads AND bit-identical raw coarse
-and unpacked owner-core fine arrays, including patch-ghost values/NaN masks.
+The former GPU gate (dense fine storage vs 1x1, 2x2 and 4x2 fine blocks,
+bit-identical) was retired with the dense route; its result and the
+dense/blocked route parity are recorded at 1573942 and 2469af7.
 """
 from contextlib import nullcontext
 import unittest
@@ -88,15 +87,13 @@ def seam_cells(layout, mapped):
 
 
 def raw(backend):
-    """Unmasked device arrays; blocked fine storage unpacked from owner cores."""
-    coarse, fine = backend._levels
-    out = [np.empty(coarse.shape, np.float32)]
-    out.append(np.empty(fine.shape if backend.fine_layout is None
-                        else backend.fine_layout.blocked_shape, np.float32))
-    cuda.memcpy_dtoh(out[0], coarse.current)
-    cuda.memcpy_dtoh(out[1], fine.current)
-    if backend.fine_layout is not None:
-        out[1] = backend.fine_layout.to_dense(out[1])
+    """Unmasked device arrays of every level; blocked storage unpacked from owner cores."""
+    out = []
+    for level in backend._levels:
+        array = np.empty(level.shape if level.layout is None else level.layout.blocked_shape,
+                         np.float32)
+        cuda.memcpy_dtoh(array, level.current)
+        out.append(array if level.layout is None else level.layout.to_dense(array))
     return out
 
 
@@ -127,7 +124,7 @@ class BlockedFineCPU(unittest.TestCase):
                 patch.object(compiler, "load_module", return_value=module), \
                 patch.object(cuda, "mem_alloc", side_effect=alloc), patch.object(cuda, "memcpy_htod"):
             backend = FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS],
-                                        level_blocks=[None, FINE_BLOCKS[0]])
+                                        level_blocks=[(12, 8), FINE_BLOCKS[0]])
             self.assertEqual(len(backend.fine_layout.halo_dst), 0)
             self.assertIsNone(backend._levels[1].halo_dst)
             backend._gather_halo(backend._levels[1], Mock())  # no launch without pairs
@@ -179,7 +176,7 @@ class BlockedFineCPU(unittest.TestCase):
                 patch.object(cuda, "memcpy_htod"), \
                 patch.object(curved_links, "CurvedLinkSoA", side_effect=AssertionError("dense SoA built")):
             backend = FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS],
-                                        level_blocks=[None, FINE_BLOCKS[1]],
+                                        level_blocks=[(12, 8), FINE_BLOCKS[1]],
                                         level_objects=[(), tuple(objects)])
             backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2))
         return backend, calls
@@ -189,8 +186,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._fine_flags, backend._levels[1].kind
-            velocity = backend._link_tables.get("surface_velocity")
+            flags, kind = backend._level_boundaries[1].flags, backend._levels[1].kind
+            velocity = backend._level_boundaries[1].device.get("surface_velocity")
             self.assertEqual(velocity is None, spec is None)
             backend.step(1)
             backend.close()
@@ -198,8 +195,9 @@ class BlockedFineCPU(unittest.TestCase):
             self.assertEqual(order, [("BlockExplode", coarse),
                                      ("BlockHaloCopy", cur), ("BlockCurvedBoundary", cur), ("BlockStep", cur),
                                      ("BlockHaloCopy", spare), ("BlockCurvedBoundary", spare),
-                                     ("BlockStep", spare), ("PatchCoarseStep", coarse)])
-            steps = [a for name, a in calls if name == "BlockStep"]
+                                     ("BlockStep", spare), ("BlockHaloCopy", coarse),
+                                     ("BlockCoalesce", cur), ("BlockStep", coarse)])
+            steps = [a for name, a in calls if name == "BlockStep"][:2]  # the two fine hops
             self.assertEqual([a[1] for a in steps], [spare, cur])
             for a in steps:
                 self.assertIs(a[3], flags)  # blocked flags, not null
@@ -215,7 +213,7 @@ class BlockedFineCPU(unittest.TestCase):
         for speed in (0.0, .002):  # any spec, even zero amplitude
             backend, calls = self.mocked_backend([circle(actuation(speed))])
             self.assertIsNotNone(backend._fine_tables["surface_velocity"])
-            del backend._link_tables["surface_velocity"]
+            del backend._level_boundaries[1].device["surface_velocity"]
             with self.assertRaisesRegex(RuntimeError, "Surface velocity device buffer is missing"):
                 backend.step(1)
             names = [name for name, _ in calls]
@@ -281,93 +279,5 @@ class BlockedFineCPU(unittest.TestCase):
             self.assertEqual(set(np.unique(kind).tolist()) - {ACTIVE, INNER_GHOST, OUTER_GHOST}, set())
 
 
-class BlockedFineGPU(unittest.TestCase):
-    def test_dense_and_blocked_fine_storage_are_bit_identical(self):
-        from CelerisLab.cuda.context import CudaContext
-        for mode in ("periodic", "slot_reservoir"):
-            cfg = config(mode)
-            patch = FixedPatch2D.from_config(cfg, BOUNDS)
-            g_coarse, g_fine = state(NY, NX, 1), state(*patch.fine_shape, 2)
-            results = []
-            with CudaContext(0) as context:
-                for fine_block in (None, *FINE_BLOCKS):
-                    with FixedPatchStepper(cfg, context=context, coarse_boundary=mode,
-                                           refinements=[BOUNDS],
-                                           level_blocks=[None, fine_block]) as backend:
-                        backend.initialize(g_coarse, g_fine)
-                        backend.step(STEPS)
-                        self.assertEqual(backend.step_count, STEPS)
-                        results.append((backend.download(), raw(backend)))
-            (dense_public, dense_raw), blocked = results[0], results[1:]
-            kind = patch.fine_kind
-            self.assertTrue(np.isnan(dense_raw[1][:, kind == OUTER_GHOST]).all(), mode)
-            # INNER populations may be NaN after hop 2 (pulled from poisoned OUTER);
-            # their exact pattern is pinned by the dense-vs-blocked comparison below.
-            self.assertTrue(np.isfinite(dense_raw[1][:, kind == ACTIVE]).all(), mode)
-            for fine_block, (public, raw_arrays) in zip(FINE_BLOCKS, blocked):
-                msg = f"{mode} fine_block={fine_block}"
-                for a, b in zip(public + tuple(raw_arrays), dense_public + tuple(dense_raw)):
-                    np.testing.assert_array_equal(a, b, err_msg=msg)  # values and NaN masks
-            print(f"{mode}: dense vs {len(blocked)} fine-block layouts bit-identical "
-                  f"after {STEPS} coarse steps", flush=True)
-
-    def test_passive_circle_dense_and_blocked_fine_storage_are_bit_identical(self):
-        from CelerisLab.cuda.context import CudaContext
-        for mode in ("periodic", "slot_reservoir"):
-            cfg = config(mode)
-            patch, boundary = fine_circle(mode)
-            solid = (boundary.flags.reshape(patch.fine_shape) & SOLID) != 0
-            g_coarse, g_fine = state(NY, NX, 1), state(*patch.fine_shape, 2)
-            results = []
-            with CudaContext(0) as context:
-                for fine_block in (None, *FINE_BLOCKS):
-                    with FixedPatchStepper(cfg, context=context, coarse_boundary=mode,
-                                           refinements=[BOUNDS], level_blocks=[None, fine_block],
-                                           level_objects=[(), (circle(),)]) as backend:
-                        backend.initialize(g_coarse, g_fine)
-                        backend.step(STEPS)
-                        results.append((backend.download(), raw(backend)))
-            (dense_public, dense_raw), blocked = results[0], results[1:]
-            fine = dense_raw[1]
-            self.assertTrue(np.isfinite(fine[:, (patch.fine_kind == ACTIVE) & ~solid]).all(), mode)
-            self.assertTrue(np.isnan(fine[:, solid]).all(), mode)
-            self.assertTrue(np.isnan(fine[:, patch.fine_kind == OUTER_GHOST]).all(), mode)
-            for fine_block, (public, raw_arrays) in zip(FINE_BLOCKS, blocked):
-                msg = f"circle {mode} fine_block={fine_block}"
-                for a, b in zip(public + tuple(raw_arrays), dense_public + tuple(dense_raw)):
-                    np.testing.assert_array_equal(a, b, err_msg=msg)  # values and NaN masks
-            print(f"circle {mode}: dense vs {len(blocked)} fine-block layouts bit-identical "
-                  f"after {STEPS} coarse steps", flush=True)
-
-    def test_actuated_circle_dense_and_blocked_fine_storage_are_bit_identical(self):
-        from CelerisLab.cuda.context import CudaContext
-        for mode in ("periodic", "slot_reservoir"):
-            for speed in (.002, -.002):
-                cfg = config(mode)
-                patch, boundary = fine_circle(mode, actuation(speed))
-                solid = (boundary.flags.reshape(patch.fine_shape) & SOLID) != 0
-                g_coarse, g_fine = state(NY, NX, 1), state(*patch.fine_shape, 2)
-                results = []
-                with CudaContext(0) as context:
-                    for fine_block in (None, *FINE_BLOCKS):
-                        with FixedPatchStepper(cfg, context=context, coarse_boundary=mode,
-                                               refinements=[BOUNDS], level_blocks=[None, fine_block],
-                                               level_objects=[(), (circle(actuation(speed)),)]) as backend:
-                            backend.initialize(g_coarse, g_fine)
-                            backend.step(STEPS)
-                            results.append((backend.download(), raw(backend)))
-                (dense_public, dense_raw), blocked = results[0], results[1:]
-                fine = dense_raw[1]
-                self.assertTrue(np.isfinite(fine[:, (patch.fine_kind == ACTIVE) & ~solid]).all())
-                self.assertTrue(np.isnan(fine[:, solid]).all())
-                self.assertTrue(np.isnan(fine[:, patch.fine_kind == OUTER_GHOST]).all())
-                for fine_block, (public, raw_arrays) in zip(FINE_BLOCKS, blocked):
-                    msg = f"actuated {speed} {mode} fine_block={fine_block}"
-                    for a, b in zip(public + tuple(raw_arrays), dense_public + tuple(dense_raw)):
-                        np.testing.assert_array_equal(a, b, err_msg=msg)  # values and NaN masks
-                print(f"actuated {speed:+.3f} {mode}: dense vs {len(blocked)} fine-block layouts "
-                      f"bit-identical after {STEPS} coarse steps", flush=True)
-
-
 if __name__ == "__main__":
     unittest.main()

tests/validation/multigrid/test_refinement_followup.py

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

tests/validation/multigrid/test_refinement_gpu.py

@@ -22,7 +22,7 @@ import CelerisLab
 from CelerisLab.config import LBMConfig
 from CelerisLab.cuda import compiler
 from CelerisLab.cuda.context import CudaContext
-from CelerisLab.lbm.refinement_stepper import FixedPatchStepper, PeriodicPatchStepper
+from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
 
 ROOT = Path(__file__).resolve().parents[3]
 POP_ATOL = 3e-6
@@ -52,12 +52,11 @@ def totals(arrays, maps):
     return np.array([populations.sum(), populations @ ref.CX, populations @ ref.CY])
 
 
-def global_array(module, name, shape):
-    out = np.empty(shape, dtype=np.uint32)
-    ptr, size = module.get_global(name)
-    assert size == out.nbytes
-    cuda.memcpy_dtoh(out, ptr)
-    return out
+def level_storage(level):
+    """Raw device storage of one blocked level as dense (9, ny, nx)."""
+    blocked = np.empty(level.layout.blocked_shape, np.float32)
+    cuda.memcpy_dtoh(blocked, level.current)
+    return level.layout.to_dense(blocked)
 
 
 class CoupledGPU(unittest.TestCase):
@@ -82,7 +81,7 @@ class CoupledGPU(unittest.TestCase):
         expected = postcollision(next_state, spec, maps)
         arrivals = (next_state.coarse_f, next_state.fine_f)
 
-        cases = (("dense_unshifted", False, {}),
+        cases = (("blocked_unshifted", False, dict(level_blocks=[(16, 16), (24, 24)])),
                  ("blocked_shifted", True,
                   dict(level_blocks=[(4, 4), (6, 6)])))
         with CudaContext(0) as context:
@@ -125,7 +124,7 @@ class CoupledGPU(unittest.TestCase):
                               "wrong-phase controls PASS", flush=True)
 
     def test_cleanup_attempts_every_allocation(self):
-        backend = PeriodicPatchStepper.__new__(PeriodicPatchStepper)
+        backend = FixedPatchStepper.__new__(FixedPatchStepper)
         backend.ctx = Mock(activate=nullcontext)
         backend._lock = RLock()
         backend._status = "failed"
@@ -147,8 +146,9 @@ class CoupledGPU(unittest.TestCase):
         spec = ref.make_spec(collision="srt")
         maps = ref.build_maps(spec)
         with CudaContext(0) as context:
-            with PeriodicPatchStepper(LBMConfig(nx=16, ny=16, viscosity=.1),
-                                      context=context, refinements=[(4, 4, 12, 12)]) as backend:
+            with FixedPatchStepper(LBMConfig(nx=16, ny=16, viscosity=.1),
+                                      context=context, refinements=[(4, 4, 12, 12)],
+                                      level_blocks=[(16, 16), (24, 24)]) as backend:
                 backend.initialize(*postcollision(ref.make_state(spec), spec, maps))
                 coarse, fine = backend._levels
 
@@ -164,122 +164,75 @@ class CoupledGPU(unittest.TestCase):
                 for operation in (backend.download, backend.step):
                     with self.assertRaisesRegex(RuntimeError, "failed"):
                         operation()
-                raw = np.empty(fine.shape, np.float32)
-                cuda.memcpy_dtoh(raw, fine.current)
+                raw = level_storage(fine)
                 self.assertTrue(np.isfinite(raw.ravel()[backend.patch.coalesce_fine]).all())
             self.assertFalse(context.is_closed)
         print("partial failure: two fine hops completed, no coarse step published; cleanup PASS", flush=True)
 
-    def test_bounded_equivalence_and_ownership(self):
-        self.assertEqual(Path(CelerisLab.__file__).resolve(), ROOT / "src/CelerisLab/__init__.py")
-        self.assertEqual(Path(compiler._KERNEL_DIR).resolve(), ROOT / "src/CelerisLab/lbm/kernels")
+    def eight_steps(self, backend, spec, maps):
+        """Rest and seeded nonequilibrium, 8 steps against the scalar reference."""
+        p = backend.patch
+        for case in ("rest", "nonequilibrium_seed1"):
+            with self.subTest(case=case):
+                state = ref.make_state(spec, velocity=(0, 0) if case == "rest" else (.03, -.02))
+                if case != "rest":
+                    rng = np.random.default_rng(1)
+                    state = ref.State(*(a * (1 + 1e-3 * rng.standard_normal(a.shape))
+                                        for a in (state.coarse_f, state.fine_f)))
+                initial = postcollision(state, spec, maps)
+                backend.initialize(*initial)
+                before = totals(backend.download(), maps)
+                max_error = max_drift = max_moment_error = max_corner = 0.0
+                for step in range(1, 9):
+                    state, diagnostic = ref.advance_coarse_once(spec, maps, state)
+                    backend.step()
+                    actual = backend.download()
+                    expected = postcollision(state, spec, maps)
+                    self.assertEqual(backend.step_count, step)
+                    p.check_synchronization(step, state.tick, phase="postcollision_g")
+                    for got, want in zip(actual, expected):
+                        np.testing.assert_array_equal(np.isnan(got), np.isnan(want))
+                        error = float(np.nanmax(np.abs(got - want)))
+                        self.assertLessEqual(error, POP_ATOL)
+                        max_error = max(max_error, error)
+                    drift = float(np.max(np.abs(totals(actual, maps) - before)))
+                    moment_error = float(np.max(np.abs(totals(actual, maps) - totals(expected, maps))))
+                    self.assertLessEqual(drift, MOMENT_ATOL)
+                    self.assertLessEqual(moment_error, MOMENT_ATOL)
+                    max_drift, max_moment_error = max(max_drift, drift), max(max_moment_error, moment_error)
+                    f = backend._levels[1]
+                    raw = level_storage(f)
+                    np.testing.assert_array_equal(np.isfinite(raw), diagnostic["hop2_valid"])
+                    self.assertTrue(np.isfinite(raw.ravel()[p.coalesce_fine]).all())
+                    self.assertEqual(len(diagnostic["corner_packets"]), 8)
+                    for packet in diagnostic["corner_packets"].values():
+                        sx, sy = packet["source_coarse_xy"]
+                        self.assertTrue(maps.coarse_active[sy, sx])
+                        values = [raw[packet["direction"], y, x]
+                                  for y, x in packet["children_fine_yx"]]
+                        error = float(np.max(np.abs(np.array(values) - packet["packets"])))
+                        self.assertLessEqual(error, POP_ATOL)
+                        max_corner = max(max_corner, error)
+                print(f"{case}: steps=8 max_population={max_error:.9g} "
+                      f"max_leaf_drift={max_drift:.9g} max_leaf_error={max_moment_error:.9g} "
+                      f"max_corner_packet={max_corner:.9g}; PASS", flush=True)
+
+    def test_bounded_equivalence_blocked(self):
+        """The same eight-step reference and launch-failure checks on whole-level tiles."""
         spec = ref.make_spec(collision="srt")
         maps = ref.build_maps(spec)
         cfg = LBMConfig(nx=16, ny=16, viscosity=.1)
-        compile_original = compiler.compile_kernel
-
-        def audit_compile(*args, **kwargs):
-            # Only this private build's entry differs. Production source unchanged.
-            entry = Path(kwargs["build_dir"]) / "refinement_audit.cu"
-            entry.write_text('#define CELERIS_REFINEMENT_AUDIT\n#include "refinement.cu"\n')
-            kwargs["source"] = entry.name
-            return compile_original(*args, **kwargs)
-
         with CudaContext(0) as context:
-            with patch.object(compiler, "compile_kernel", audit_compile):
-                backend = PeriodicPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)])
-            with backend:
-                # Production allocations: four DDF arrays and four compact tables.
-                device_bytes = 2 * 9 * (16**2 + 24**2) * 4
-                p = backend.patch
-                device_bytes += sum(a.nbytes for a in (p.ghost_fine, p.ghost_coarse,
-                                                       p.coalesce_coarse, p.coalesce_fine))
-                audit_bytes = sum((2 + 9) * np.prod(level.shape[1:]) * 4
-                                  for level in backend._levels)
-                self.assertLess(device_bytes + audit_bytes, 128 * 1024)
-                print(f"device_arrays={device_bytes} audit_arrays={audit_bytes} bytes", flush=True)
-                for case in ("rest", "nonequilibrium_seed1"):
-                    with self.subTest(case=case):
-                        state = ref.make_state(spec, velocity=(0, 0) if case == "rest" else (.03, -.02))
-                        if case != "rest":
-                            rng = np.random.default_rng(1)
-                            state = ref.State(*(a * (1 + 1e-3 * rng.standard_normal(a.shape))
-                                                for a in (state.coarse_f, state.fine_f)))
-                        initial = postcollision(state, spec, maps)
-                        backend.initialize(*initial)
-                        for level in backend._levels:
-                            for name in ("collision_counts", "explosion_counts", "incoming_counts"):
-                                ptr, size = level.module.get_global(name)
-                                cuda.memset_d8(ptr, 0, size)
-                        before = totals(backend.download(), maps)
-                        max_error = max_drift = max_moment_error = max_corner = 0.0
-                        for step in range(1, 9):
-                            state, diagnostic = ref.advance_coarse_once(spec, maps, state)
-                            backend.step()
-                            actual = backend.download()
-                            expected = postcollision(state, spec, maps)
-                            self.assertEqual(backend.step_count, step)
-                            p.check_synchronization(step, state.tick, phase="postcollision_g")
-                            for got, want in zip(actual, expected):
-                                np.testing.assert_array_equal(np.isnan(got), np.isnan(want))
-                                error = float(np.nanmax(np.abs(got - want)))
-                                self.assertLessEqual(error, POP_ATOL)
-                                max_error = max(max_error, error)
-                            drift = float(np.max(np.abs(totals(actual, maps) - before)))
-                            moment_error = float(np.max(np.abs(totals(actual, maps) - totals(expected, maps))))
-                            self.assertLessEqual(drift, MOMENT_ATOL)
-                            self.assertLessEqual(moment_error, MOMENT_ATOL)
-                            max_drift, max_moment_error = max(max_drift, drift), max(max_moment_error, moment_error)
-                            c, f = backend._levels
-                            np.testing.assert_array_equal(global_array(c.module, "collision_counts", (16, 16)),
-                                                          step * maps.coarse_active)
-                            np.testing.assert_array_equal(global_array(f.module, "collision_counts", (24, 24)),
-                                                          2 * step * maps.fine_active)
-                            np.testing.assert_array_equal(global_array(f.module, "explosion_counts", (24, 24)),
-                                                          step * ~maps.fine_active)
-                            owners = np.zeros((9, 16, 16), np.uint32)
-                            owners.ravel()[p.coalesce_coarse] = step
-                            np.testing.assert_array_equal(global_array(c.module, "incoming_counts", owners.shape), owners)
-                            raw = np.empty(f.shape, np.float32)
-                            cuda.memcpy_dtoh(raw, f.current)
-                            np.testing.assert_array_equal(np.isfinite(raw), diagnostic["hop2_valid"])
-                            self.assertTrue(np.isfinite(raw.ravel()[p.coalesce_fine]).all())
-                            self.assertEqual(len(diagnostic["corner_packets"]), 8)
-                            for packet in diagnostic["corner_packets"].values():
-                                sx, sy = packet["source_coarse_xy"]
-                                self.assertTrue(maps.coarse_active[sy, sx])
-                                values = [raw[packet["direction"], y, x]
-                                          for y, x in packet["children_fine_yx"]]
-                                error = float(np.max(np.abs(np.array(values) - packet["packets"])))
-                                self.assertLessEqual(error, POP_ATOL)
-                                max_corner = max(max_corner, error)
-                        print(f"{case}: steps=8 max_population={max_error:.9g} "
-                              f"max_leaf_drift={max_drift:.9g} max_leaf_error={max_moment_error:.9g} "
-                              f"max_corner_packet={max_corner:.9g}; ownership/counts PASS", flush=True)
-
-                # Fail before any additional time step; no stale synchronized export.
+            with FixedPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)],
+                                      level_blocks=[(16, 16), (24, 24)]) as backend:
+                self.eight_steps(backend, spec, maps)
                 with patch.object(backend._pairs[0], "explode", side_effect=RuntimeError("injected launch failure")):
                     with self.assertRaisesRegex(RuntimeError, "injected"):
                         backend.step()
-                with self.assertRaisesRegex(RuntimeError, "failed"):
-                    backend.download()
-                with self.assertRaisesRegex(RuntimeError, "failed"):
-                    backend.step()
+                for operation in (backend.download, backend.step):
+                    with self.assertRaisesRegex(RuntimeError, "failed"):
+                        operation()
             self.assertFalse(context.is_closed)
-            with self.assertRaisesRegex(RuntimeError, "closed"):
-                backend.download()
-
-            # Compile and execute the actual non-audited entry as well: one step.
-            with PeriodicPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)]) as plain:
-                rest = postcollision(ref.make_state(spec), spec, maps)
-                plain.initialize(*rest)
-                plain.step()
-                for got, expected in zip(plain.download(), rest):
-                    np.testing.assert_allclose(got, expected, rtol=0, atol=POP_ATOL, equal_nan=True)
-                print("normal-build rest: steps=1 PASS; source/kernel paths pinned", flush=True)
-        rss = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024
-        print(f"host_peak_rss={rss:.3f} MiB", flush=True)
-
 
 if __name__ == "__main__":
     unittest.main(verbosity=2)

tests/validation/multigrid/test_refinement_jets.py

@@ -3,11 +3,13 @@
 Jets use the same coarse cell-edge frame as fine_objects. The equivalent jet
 of Astra's actuated fixture is SurfaceActuation(0.7, 0.85, U).as_jet((16, 12),
 2.15). Frozen GPU gates (8 coarse steps, both coarse boundary modes):
-  G1 jet run: dense vs 1x1/2x2/4x2 fine blocks bit-identical (values + NaN masks);
-  G2 jet vs legacy actuation (dense): max |dG| over finite fine ACTIVE fluid
+  G1 jet run: 1x1/2x2/4x2 fine blocks (blocked (12, 8) root) bit-identical to the
+     FINE_BLOCKS[0] reference (values + NaN masks);
+  G2 jet vs legacy actuation (blocked reference): max |dG| over finite fine ACTIVE fluid
      populations <= 1e-6 (velocity rows differ only by FP32 rounding; reported);
-  G3 amplitude change after 4 steps (+U -> -U), 4 more steps: dense vs blocked
-     bit-identical and the run differs from the unchanged run.
+  G3 amplitude change after 4 steps (+U -> -U), 4 more steps: blocked layouts
+     bit-identical to the FINE_BLOCKS[0] reference and the run differs from the
+     unchanged run.
 """
 from contextlib import nullcontext
 import unittest
@@ -26,6 +28,11 @@ from test_refinement_blocked_fine import (BOUNDS, FINE_BLOCKS, NX, NY, STEPS, ac
 SPEED = .002
 
 
+def root_and(fine_block):
+    """Blocked levels with a (12, 8) root tile."""
+    return [(12, 8), fine_block]
+
+
 def jet(speed=SPEED):
     return actuation(speed).as_jet((16.0, 12.0), 2.15)
 
@@ -42,7 +49,7 @@ def mocked(objects, jets, fine_block):
             patch.object(cuda, "mem_alloc", side_effect=lambda size: Mock()), \
             patch.object(cuda, "memcpy_htod", htod):
         backend = FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS],
-                                    level_blocks=[None, fine_block], level_objects=[(), tuple(objects)],
+                                    level_blocks=root_and(fine_block), level_objects=[(), tuple(objects)],
                                     level_jets=[(), tuple(jets)])
     return backend, htod
 
@@ -100,7 +107,7 @@ class JetStepperCPU(unittest.TestCase):
     def test_blocked_amplitude_update_rewrites_only_velocity_column(self):
         backend, htod = mocked([circle()], [jet()], FINE_BLOCKS[1])
         before = backend._fine_tables["surface_velocity"].copy()
-        device = backend._link_tables["surface_velocity"]
+        device = backend._level_boundaries[1].device["surface_velocity"]
         with patch.object(cuda, "memcpy_htod", htod):
             htod.reset_mock()
             backend.set_jet_amplitude(0, -SPEED)
@@ -116,7 +123,7 @@ class JetStepperCPU(unittest.TestCase):
         backend.close()
 
     def test_copy_failure_marks_failed_and_keeps_host_amplitude(self):
-        for fine_block in (None, FINE_BLOCKS[1]):  # dense SoA and blocked tables
+        for fine_block in (FINE_BLOCKS[1],):  # blocked tables
             backend, _ = mocked([circle()], [jet()], fine_block)
             status = backend._status
             for bad, error in ((1.5, ValueError), (float("nan"), ValueError)):
@@ -126,9 +133,7 @@ class JetStepperCPU(unittest.TestCase):
             with self.assertRaises(IndexError):
                 backend.set_jet_amplitude(1, -SPEED)
             self.assertEqual(backend._status, status)
-            dense = backend._fine_curved
-            before = (dense.surface_velocity if dense is not None
-                      else backend._fine_tables["surface_velocity"]).copy()
+            before = backend._fine_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)):
@@ -136,11 +141,7 @@ class JetStepperCPU(unittest.TestCase):
                     backend.set_jet_amplitude(0, -SPEED)
             self.assertEqual(backend._status, "failed", fine_block)
             self.assertEqual(backend._jets[0].amplitude, SPEED)
-            if dense is not None:
-                self.assertEqual(dense.jets[0].amplitude, SPEED)
-                np.testing.assert_array_equal(dense.surface_velocity, before)
-            else:
-                np.testing.assert_array_equal(backend._fine_tables["surface_velocity"], before)
+            np.testing.assert_array_equal(backend._fine_tables["surface_velocity"], before)
             with self.assertRaisesRegex(RuntimeError, "failed"):
                 backend.step(1)
             with self.assertRaisesRegex(RuntimeError, "failed"):
@@ -151,6 +152,7 @@ class JetStepperCPU(unittest.TestCase):
         with patch.object(cuda, "Stream", side_effect=AssertionError("CUDA touched")):
             with self.assertRaisesRegex(ValueError, "SurfaceActuation"):
                 FixedPatchStepper(config("periodic"), context=Mock(), refinements=[BOUNDS],
+                                  level_blocks=root_and(FINE_BLOCKS[1]),
                                   level_objects=[(), (circle(actuation(SPEED)),)],
                                   level_jets=[(), (jet(),)])
 
@@ -161,7 +163,7 @@ def run(mode, fine_block, *, objects, jets, switch=None):
     patch_ = FixedPatch2D.from_config(cfg, BOUNDS)
     with CudaContext(0) as context, FixedPatchStepper(
             cfg, context=context, coarse_boundary=mode, refinements=[BOUNDS],
-            level_blocks=[None, fine_block], level_objects=[(), tuple(objects)],
+            level_blocks=root_and(fine_block), level_objects=[(), tuple(objects)],
             level_jets=[(), tuple(jets)]) as backend:
         backend.initialize(state(NY, NX, 1), state(*patch_.fine_shape, 2))
         if switch is None:
@@ -179,29 +181,29 @@ class JetStepperGPU(unittest.TestCase):
             patch_, boundary = fine_circle(mode)
             fluid = (patch_.fine_kind == ACTIVE) & (
                 (boundary.flags.reshape(patch_.fine_shape) & SOLID) == 0)
-            dense = run(mode, None, objects=[circle()], jets=[jet()])
-            self.assertTrue(np.isfinite(dense[1][1][:, fluid]).all(), mode)
-            for fine_block in FINE_BLOCKS:
+            reference = run(mode, FINE_BLOCKS[0], objects=[circle()], jets=[jet()])
+            self.assertTrue(np.isfinite(reference[1][1][:, fluid]).all(), mode)
+            for fine_block in FINE_BLOCKS[1:]:
                 public, raw_arrays = run(mode, fine_block, objects=[circle()], jets=[jet()])
-                for a, b in zip(public + tuple(raw_arrays), dense[0] + tuple(dense[1])):
+                for a, b in zip(public + tuple(raw_arrays), reference[0] + tuple(reference[1])):
                     np.testing.assert_array_equal(a, b, err_msg=f"G1 {mode} {fine_block}")
-            legacy = run(mode, None, objects=[circle(actuation(SPEED))], jets=())
-            diff = np.abs(dense[1][1][:, fluid] - legacy[1][1][:, fluid]).max()
+            legacy = run(mode, FINE_BLOCKS[0], objects=[circle(actuation(SPEED))], jets=())
+            diff = np.abs(reference[1][1][:, fluid] - legacy[1][1][:, fluid]).max()
             print(f"G2 {mode}: jet vs legacy fine max |dG| = {diff:.3e}", flush=True)
             self.assertLessEqual(diff, 1e-6, mode)
 
     def test_g3_amplitude_switch_layouts_identical(self):
         for mode in ("periodic", "slot_reservoir"):
-            dense = run(mode, None, objects=[circle()], jets=[jet()], switch=-SPEED)
-            steady = run(mode, None, objects=[circle()], jets=[jet()])
-            self.assertFalse(np.array_equal(dense[1][1], steady[1][1], equal_nan=True), mode)
-            for fine_block in FINE_BLOCKS:
+            reference = run(mode, FINE_BLOCKS[0], objects=[circle()], jets=[jet()], switch=-SPEED)
+            steady = run(mode, FINE_BLOCKS[0], objects=[circle()], jets=[jet()])
+            self.assertFalse(np.array_equal(reference[1][1], steady[1][1], equal_nan=True), mode)
+            for fine_block in FINE_BLOCKS[1:]:
                 public, raw_arrays = run(mode, fine_block, objects=[circle()], jets=[jet()],
                                          switch=-SPEED)
-                for a, b in zip(public + tuple(raw_arrays), dense[0] + tuple(dense[1])):
+                for a, b in zip(public + tuple(raw_arrays), reference[0] + tuple(reference[1])):
                     np.testing.assert_array_equal(a, b, err_msg=f"G3 {mode} {fine_block}")
-            print(f"G3 {mode}: amplitude switch dense vs {len(FINE_BLOCKS)} layouts identical",
-                  flush=True)
+            print(f"G3 {mode}: amplitude switch, {len(FINE_BLOCKS) - 1} layouts "
+                  f"identical to {FINE_BLOCKS[0]}", flush=True)
 
 
 if __name__ == "__main__":

tests/validation/multigrid/test_refinement_mrt_gpu.py

@@ -54,25 +54,6 @@ class RefinementMRTGPU(unittest.TestCase):
                             else:
                                 for a, b in zip(reference, raw.values): np.testing.assert_array_equal(a, b)
 
-    def test_dense_versus_blocked_mrt_two_level(self):
-        """Unshifted MRT: all-dense PatchCoarseStep/PatchFineStep vs mixed/all-blocked storage."""
-        from test_refinement_blocked_base import BASE_BLOCKS, FINE_BLOCK
-        from test_refinement_blocked_fine import NX, NY, BOUNDS
-        import test_refinement_blocked_fine as fixture
-        cfg = LBMConfig(nx=NX, ny=NY, viscosity=.1, collision="MRT")
-        with CudaContext(0) as context:
-            results = {}
-            layouts = ((None, None), *((base, FINE_BLOCK) for base in (None, *BASE_BLOCKS)))
-            for base_block, fine_block in layouts:
-                with FixedPatchStepper(cfg, context=context, refinements=[BOUNDS],
-                                       level_blocks=[base_block, fine_block]) as backend:
-                    backend.initialize(*[fixture.state(*level.shape[1:], seed) for seed, level in enumerate(backend._levels)])
-                    backend.step(4)
-                    results[base_block, fine_block] = raw(backend)
-            for layout in layouts[1:]:
-                for level, (a, b) in enumerate(zip(results[layout], results[None, None])):
-                    exact(self, a, b, f"layout={layout} level={level}")
-
     def test_shifted_all_blocked_three_level_mrt_smoke(self):
         """Per-level shear rate omega_l = 1/(3*2**l*nu + 1/2); finite ACTIVE leaves."""
         cfg = LBMConfig(nx=16, ny=16, viscosity=.1, collision="MRT", ddf_shifting=True)

tests/validation/multigrid/test_refinement_three_level.py

@@ -75,7 +75,7 @@ class ThreeLevelCPU(unittest.TestCase):
             with self.subTest(kw=kw), self.assertRaises(ValueError):
                 FixedPatchStepper(CFG, context=context, **kw)
         # L2 bodies/Jets without an L2 level (one entry too many), and an L2 Jet without a body.
-        two_level = dict(refinements=[L1], level_blocks=[None, None])
+        two_level = dict(refinements=[L1], level_blocks=[B, B])
         for kw in (dict(two_level, level_objects=[(), (), (child_body(),)]),
                    dict(two_level, level_jets=[(), (), (child_jet(),)]),
                    levels(jets=((), (), (child_jet(),)))):