test: verify Simulation Jet amplitude uploads on GPU

Simulation.set_jet_amplitude validated before upload and committed host
state after it, but after an uncertain device write nothing stopped reuse.
Record the failure in Simulation._failed and refuse every public device
call except close() (FixedPatchStepper already uses status 'failed').
Invalid input stays side-effect free and retryable. Adds host fault tests
for both owners and states the rule in docs/lifecycle.md.

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

docs/lifecycle.md

@@ -5,9 +5,9 @@ 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. |
+| 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`). | Input (index, level, amplitude) is validated before any device write; invalid input changes nothing and the instance stays usable. The host copy commits only after the upload returned. Any error after validation leaves the device column undefined and fails the instance: `FixedPatchStepper` status `failed`, `Simulation._failed`; every later call except `close()` raises, `close()` still releases resources and stays retryable. |
 | 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. |
 
 Storage routes are separate concerns:
@@ -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

@@ -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, PeriodicPatchStepper(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()

src/CelerisLab/simulation.py

@@ -47,6 +47,7 @@ class Simulation:
         # Internal stream for async DTOH coordination
         self._stream: Optional[cuda.Stream] = None
         self._closed = False
+        self._failed = None  # reason after an uncertain device write; only close() remains
 
         # Load configs
         self.lbm_cfg = load_lbm_config(lbm_config_path)
@@ -172,15 +173,26 @@ class Simulation:
     def set_jet_amplitude(self, index: int, amplitude: float) -> None:
         """Change one jet's signed amplitude; re-uploads per-link velocity only.
 
-        Invalid input changes nothing. A device error leaves the host amplitude
-        unchanged but the device column undefined: treat it like any other CUDA
-        exception here and close/recreate the Simulation before stepping again.
+        Invalid input raises before any change. Any error after validation
+        leaves the host amplitude unchanged but the device column undefined,
+        so the Simulation is marked failed: every later call except close()
+        raises; close and recreate it.
         """
+        self._require_usable()
         jet = self.bodies.jets[index].with_amplitude(amplitude)
         if self._initialized:
-            self.field.curved.set_jet_amplitude(index, amplitude, self.stream)
+            try:
+                self.field.curved.set_jet_amplitude(index, amplitude, self.stream)
+            except BaseException as exc:
+                self._failed = f"jet amplitude upload failed: {exc!r}"
+                raise
         self.bodies.jets[index] = jet  # later geometry syncs keep the new amplitude
 
+    def _require_usable(self):
+        failed = getattr(self, "_failed", None)
+        if failed is not None:
+            raise RuntimeError(f"Simulation failed ({failed}); close() and recreate it")
+
     def remove_body(self, id: int) -> None:
         """Remove a body from the simulation.
 
@@ -338,6 +350,7 @@ class Simulation:
         Unlike read_body this is not accumulated/time-averaged. It resets after
         geometry sync or checkpoint load; advance before interpreting it.
         """
+        self._require_usable()
         obj = self.bodies.get(id)
         curved = self.field.curved
         has_aperture = (curved.aperture is not None
@@ -365,6 +378,7 @@ class Simulation:
             BodyTelemetry with fields ``force``, ``torque``, ``sensor``.
             ``sensor`` is a dim-length array of zeros for non-sensor bodies.
         """
+        self._require_usable()
         if self.bodies.obs_pinned is None:
             raise RuntimeError("No obs buffer. Call run() first.")
         if stream is None:
@@ -384,6 +398,7 @@ class Simulation:
         This is the most efficient readback for DRL control loops
         (single contiguous D2H copy, no Python per-body overhead).
         """
+        self._require_usable()
         if self.bodies.obs_pinned is None:
             raise RuntimeError("No obs buffer. Call run() first.")
         if stream is None:
@@ -446,6 +461,7 @@ class Simulation:
     @context_bound
     def initialize(self):
         """Initialize flow field and sync objects to GPU."""
+        self._require_usable()
         self._assert_allocation_contract()
         # Validate the complete geometry before flags, populations or modules
         # change. Compact-list construction publishes no outputs on failure.
@@ -478,6 +494,7 @@ class Simulation:
         In EsoPull mode, the backing layout is decoded to physical DDF on
         GPU for the patched cells only, patched on host, then encoded back.
         """
+        self._require_usable()
         if not self._initialized:
             raise RuntimeError("Call initialize() before sync_bodies()")
 
@@ -573,6 +590,7 @@ class Simulation:
                 sensor) on GPU and reset the step accumulator before stepping.
             checkpoint_interval: Nonnegative integer; if >0, save every N steps.
         """
+        self._require_usable()
         if not self._initialized:
             raise RuntimeError("Call initialize() first")
         steps = operator.index(steps)
@@ -637,10 +655,12 @@ class Simulation:
     @context_bound
     def get_macroscopic(self) -> Dict[str, np.ndarray]:
         """Download DDF and return rho, ux, uy [, uz]."""
+        self._require_usable()
         return self.field.get_macroscopic(step_id=self.stepper.step_count)
 
     @context_bound
     def get_ddf(self) -> np.ndarray:
+        self._require_usable()
         self.field.download_ddf(step_id=self.stepper.step_count)
         return self.field.ddf.copy()
 
@@ -654,23 +674,27 @@ class Simulation:
 
         Accepted: omega, fx, fy, fz, n_objects.
         """
+        self._require_usable()
         self.field.update_params(**kwargs)
 
     # -- Snapshots -----------------------------------------------------------
     @context_bound
     def snapshot(self):
         """Snapshot the flow field; restore keeps the current clock and controls."""
+        self._require_usable()
         self.field.snapshot(step_id=self.stepper.step_count)
 
     @context_bound
     def restore(self):
         """Restore the saved flow field at the current streaming parity."""
+        self._require_usable()
         self.field.restore(step_id=self.stepper.step_count)
 
     # -- Checkpoint / Restore (HDF5) ----------------------------------------
     @context_bound
     def save_checkpoint(self, path: str = None):
         """Save full simulation state to an HDF5 file."""
+        self._require_usable()
         from .common.checkpoint import save_checkpoint
         if self.bodies.jets:
             raise NotImplementedError("Jets are not checkpointed yet")
@@ -680,6 +704,7 @@ class Simulation:
     @context_bound
     def load_checkpoint(self, path: str):
         """Restore simulation state from an HDF5 checkpoint."""
+        self._require_usable()
         from .common.checkpoint import load_checkpoint
         if not self._initialized:
             self.initialize()

tests/integration/test_jet_update_gpu.py

@@ -0,0 +1,53 @@
+"""Direct Simulation Jet control upload; boundary-law accuracy is tested elsewhere."""
+import numpy as np
+import pycuda.driver as cuda
+
+from CelerisLab.body.jets import Jet, Radial, Sector
+from CelerisLab.lbm.descriptors import FLUID
+from ._cases import make_sim
+
+
+def test_simulation_jet_amplitude_upload_and_continued_steps():
+    # Same grid/circle/sector as the existing uniform surface-actuation fixture.
+    center = (16., 12.)
+    with make_sim(nx=32, ny=24, velocity=0) as sim:
+        sim.add_body("circle", center=center, radius=4.3)
+        index = sim.add_jet(Jet(Sector(center, .7, .85, 4.4),
+                                Radial(center, inward=False), "quartic", .002))
+        sim.initialize()
+        curved = sim.field.curved
+        shape, jet_id = curved.jet_shape.copy(), curved.jet_id.copy()
+        flags = sim.get_flags()
+        driven = jet_id == index
+        assert driven.any() and (~driven).any()
+        velocity_pointer = int(curved.surface_velocity_gpu)
+        previous = None
+        for interval, amplitude in enumerate((.002, -.001)):
+            if interval:
+                count = sim.stepper.step_count
+                sim.set_jet_amplitude(index, amplitude)
+                assert sim.stepper.step_count == count  # update does not step
+            expected = np.zeros_like(curved.surface_velocity)
+            expected[driven] = (amplitude * shape[driven]).astype(np.float32)
+            actual = np.empty_like(expected)
+            with sim.ctx.activate():
+                cuda.memcpy_dtoh(actual, curved.surface_velocity_gpu)
+            np.testing.assert_array_equal(actual, expected)  # actual device column
+            np.testing.assert_array_equal(curved.surface_velocity, expected)
+            assert np.any(actual[driven] != 0)
+            if previous is not None:
+                assert not np.array_equal(actual[driven], previous[driven])
+            previous = actual.copy()
+            assert sim.bodies.jets[index].amplitude == amplitude
+            assert curved.jets[index].amplitude == amplitude
+            assert int(curved.surface_velocity_gpu) == velocity_pointer
+            np.testing.assert_array_equal(curved.jet_shape, shape)
+            np.testing.assert_array_equal(curved.jet_id, jet_id)
+            np.testing.assert_array_equal(sim.get_flags(), flags)
+            sim.run(2)
+            assert sim.stepper.step_count == 2 * (interval + 1)
+            assert sim._failed is None
+            macro = sim.get_macroscopic()
+            fluid = (flags.reshape(24, 32) & FLUID) != 0
+            assert all(np.isfinite(value[fluid]).all() for value in macro.values())
+            assert np.all(macro["rho"][fluid] > 0)

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_jet_update_failure.py

@@ -0,0 +1,124 @@
+"""Jet amplitude update failure contract (host doubles only, no CUDA numerics).
+
+Both owners validate before any device write, commit host state only after the
+upload returns, and refuse every later call except close() after an uncertain
+device write.
+"""
+from types import SimpleNamespace
+
+import numpy as np
+import pytest
+
+pytest.importorskip("pycuda.driver")
+from CelerisLab.body.jets import Jet, Radial, Sector
+from tests.unit.test_pinned_body_remap import ready
+from tests.unit.test_simulation_lifecycle import make_simulation
+
+BAD = (float("nan"), float("inf"), 1.)
+
+
+class Curved:
+    """Records uploads; ``fail`` raises once the write was attempted."""
+    def __init__(self):
+        self.writes, self.fail = [], False
+
+    def set_jet_amplitude(self, index, amplitude, stream):
+        self.writes.append((index, amplitude))
+        if self.fail:
+            raise RuntimeError("upload failed")
+
+
+def simulation(monkeypatch):
+    sim = make_simulation(monkeypatch)
+    jet = Jet(Sector((6., 5.), 0., .5, 4.), Radial((6., 5.)), "uniform", 1e-3)
+    sim._failed, sim._initialized = None, True
+    sim.bodies, sim.field = SimpleNamespace(jets=[jet]), SimpleNamespace(curved=Curved())
+    return sim, jet
+
+
+def test_simulation_invalid_input_writes_nothing_and_stays_usable(monkeypatch):
+    sim, jet = simulation(monkeypatch)
+    for amplitude in BAD:
+        with pytest.raises(ValueError):
+            sim.set_jet_amplitude(0, amplitude)
+    with pytest.raises(IndexError):
+        sim.set_jet_amplitude(1, 2e-3)
+    assert sim.field.curved.writes == [] and sim.bodies.jets[0] is jet and sim._failed is None
+    sim.set_jet_amplitude(0, 2e-3)  # retry after a rejected input
+    assert sim.field.curved.writes == [(0, 2e-3)] and sim.bodies.jets[0].amplitude == 2e-3
+    sim.close()
+
+
+def test_simulation_upload_failure_keeps_host_and_blocks_reuse(monkeypatch):
+    sim, jet = simulation(monkeypatch)
+    sim.field.curved.fail = True
+    with pytest.raises(RuntimeError, match="upload failed"):
+        sim.set_jet_amplitude(0, 2e-3)
+    assert sim.bodies.jets[0] is jet and sim._failed is not None
+    sim.field.curved.fail = False
+    for call in (lambda: sim.set_jet_amplitude(0, 2e-3), lambda: sim.run(1), sim.step,
+                 sim.sync_bodies, sim.initialize, sim.get_ddf, sim.get_macroscopic,
+                 lambda: sim.read_surface_exchange(0), lambda: sim.read_body(0),
+                 sim.read_bodies, sim.snapshot, sim.restore, sim.save_checkpoint,
+                 lambda: sim.load_checkpoint("x.h5"), lambda: sim.update_runtime_params(omega=1.)):
+        with pytest.raises(RuntimeError, match="Simulation failed"):
+            call()
+    assert sim.field.curved.writes == [(0, 2e-3)]
+    sim.close()  # cleanup still works
+    assert sim._closed and sim.ctx.is_closed
+
+
+def test_simulation_before_initialize_commits_host_only(monkeypatch):
+    sim, _ = simulation(monkeypatch)
+    sim._initialized = False
+    sim.set_jet_amplitude(0, -2e-3)
+    assert sim.field.curved.writes == [] and sim.bodies.jets[0].amplitude == -2e-3
+    sim.close()
+
+
+def stepper_state(run):
+    boundary = run._level_boundaries[2]
+    return boundary.tables["surface_velocity"].copy(), boundary.jets, run._level_jets[2]
+
+
+def test_stepper_invalid_input_writes_nothing_and_stays_ready():
+    run, module, _, _, _ = ready()
+    try:
+        velocity, jets, level_jets = stepper_state(run)
+        module.cuda.memcpy_htod.reset_mock()
+        for args, kw, error in [((0, a), {"level": 2}, ValueError) for a in BAD] + [
+                ((1, .01), {"level": 2}, IndexError), ((0, .01), {"level": 3}, ValueError)]:
+            with pytest.raises(error):
+                run.set_jet_amplitude(*args, **kw)
+        module.cuda.memcpy_htod.assert_not_called()
+        np.testing.assert_array_equal(run._level_boundaries[2].tables["surface_velocity"], velocity)
+        assert (run._level_boundaries[2].jets, run._level_jets[2]) == (jets, level_jets)
+        assert run._status == "ready"
+        run.set_jet_amplitude(0, .017, level=2)
+        assert run._level_jets[2][0].amplitude == .017
+    finally:
+        run.close()
+
+
+@pytest.mark.parametrize("where", ["upload", "synchronize"])
+def test_stepper_device_failure_keeps_host_and_blocks_reuse(where):
+    run, module, _, _, _ = ready()
+    try:
+        velocity, jets, level_jets = stepper_state(run)
+        if where == "upload":
+            module.cuda.memcpy_htod.side_effect = RuntimeError("upload failed")
+        else:
+            run._stream.synchronize.side_effect = RuntimeError("upload failed")
+        with pytest.raises(RuntimeError, match="upload failed"):
+            run.set_jet_amplitude(0, .017, level=2)
+        np.testing.assert_array_equal(run._level_boundaries[2].tables["surface_velocity"], velocity)
+        assert (run._level_boundaries[2].jets, run._level_jets[2]) == (jets, level_jets)
+        assert run._status == "failed"
+        module.cuda.memcpy_htod.side_effect = None
+        run._stream.synchronize.side_effect = None
+        for call in (lambda: run.set_jet_amplitude(0, .017, level=2), run.step, run.download):
+            with pytest.raises(RuntimeError, match="failed"):
+                call()
+    finally:
+        run.close()
+    assert run._status == "closed"

tests/unit/test_jets.py

@@ -194,6 +194,7 @@ def test_public_surface_exchange_for_jet_preserves_aperture_at_zero_and_reverse(
     field, bodies = jet_case(device, [jet])
     sim = SimpleNamespace(bodies=bodies, field=field, _initialized=True,
                           stream=SimpleNamespace(synchronize=device.synchronize))
+    sim._require_usable = Simulation._require_usable.__get__(sim)
     cl = field.curved
     assert bodies.get(0).surface_actuation is None
     aperture, velocity = cl.aperture, cl.surface_velocity.copy()
@@ -218,6 +219,7 @@ def test_public_surface_exchange_keeps_legacy_and_passive_guards(device):
     bodies.sync_to_gpu(field)
     sim = SimpleNamespace(bodies=bodies, field=field,
                           stream=SimpleNamespace(synchronize=device.synchronize))
+    sim._require_usable = Simulation._require_usable.__get__(sim)
     for value in Simulation.read_surface_exchange(sim, 0).values():
         np.testing.assert_array_equal(value, 0.)
     field.curved.free_gpu_columns()

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/unit/test_state_regressions.py

@@ -295,6 +295,7 @@ def test_passive_overlap_failure_keeps_live_state(device, entry, second_x):
                           _initialized=entry == "runtime_edit",
                           _assert_allocation_contract=lambda: None,
                           stream=SimpleNamespace(synchronize=device.synchronize))
+    sim._require_usable = Simulation._require_usable.__get__(sim)
     with pytest.raises(ValueError, match="active fluid receivers|unique cut-link writers"):
         if entry == "sync":
             bodies.sync_to_gpu(field)

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,24 +71,6 @@ def uniform_reference(speed, roi):
         return traces
 
 
-@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)
-    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)
-        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)
-
-
 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)
@@ -334,46 +126,6 @@ def test_fine_surface_blocked_matches_uniform_four_hops(speed):
     print(f"blocked speed={speed} four hops vs uniform PASS mass_into_fluid={mass:+.9g}", flush=True)
 
 
-@pytest.mark.parametrize("missing", ["surface_velocity_gpu", "link_exchange_gpu"])
-def test_missing_surface_device_buffer_rejects_before_launch(missing):
-    backend = implementation.FixedPatchStepper.__new__(implementation.FixedPatchStepper)
-    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()
-
-
 def test_cross_interface_donor_rejects_before_cuda():
     cfg = LBMConfig(nx=24, ny=24, viscosity=.05)
     source = SimObject(0, CircleGeometry(5.05, 12.25, .7),

tests/validation/multigrid/test_raw_exchange_gpu.py

@@ -92,9 +92,8 @@ def census(run, index, buffer):
 
 
 def fine(block):
-    """Two-level circle+jet construction on BOUNDS, all-blocked with a (12, 8) root tile;
-    block None is the dense route (dense reference test only, retired with the dense route)."""
-    return dict(refinements=[BOUNDS], level_blocks=[None if block is None else (12, 8), block],
+    """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):
@@ -147,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,36 +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
-  MRT, periodic root, circle/actuation/Jet      yes       yes       yes   (parity gate)
-  two fine bodies (Jet + passive, two actuated) yes       yes       yes   (parity gate)
-  MRT with slot/reservoir root                  rejected on every route (on_site_velocity is SRT-only)
-Slot/reservoir walls are compile-time WALL_* macros on a dense root and the
-d_block_wall constant on a blocked root: BC-encoding parity, bit-identical by
-this gate. Dense fine bodies use CurvedLinkSoA + CurvedBoundaryKernel
-(n_objects > 0), blocked ones BlockCurvedBoundary; same links and physics,
-compared through the population rows. MRT body-free is also covered by
-test_refinement_mrt_gpu::test_dense_versus_blocked_mrt_two_level. Parked:
-shifted and depth >= 3 (no dense route to compare).
+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
@@ -51,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)
@@ -298,224 +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]))
-# Two-body rows: patch (4, 4, 28, 20) gives a 40 x 56 fine window, tiled (8, 8).
-WIDE_BOUNDS = (4, 4, 28, 20)
-WIDE_ROUTES = (("D/D", [None, None]), ("D/B", [None, (8, 8)]), ("B/B", [(12, 8), (8, 8)]))
-LEFT, RIGHT, RADIUS = (12.0, 12.0), (20.0, 12.0), 2.15
-
-
-def body(ident, center, spec=None):
-    from CelerisLab.body.geometry.circle import CircleGeometry
-    from CelerisLab.body.objects import SimObject
-    return SimObject(ident, geometry=CircleGeometry(*center, RADIUS), center=center, radius=RADIUS,
-                     surface_actuation=spec)
-
-
-TWO_BODY_CASES = (
-    ("jet + passive", lambda: [body(0, LEFT), body(1, RIGHT)],
-     lambda: [actuation(SPEED).as_jet(LEFT, RADIUS)], -SPEED),
-    ("actuated +/-", lambda: [body(0, LEFT, actuation(SPEED)), body(1, RIGHT, actuation(-SPEED))],
-     lambda: [], None),
-)
-
-
-def same_summary(a, b):
-    """Exact equality of nested exchange summaries (dicts, sequences, arrays, numbers)."""
-    if isinstance(a, dict):
-        return a.keys() == b.keys() and all(same_summary(a[k], b[k]) for k in a)
-    if isinstance(a, (list, tuple)):
-        return len(a) == len(b) and all(same_summary(x, y) for x, y in zip(a, b))
-    a, b = np.asarray(a), np.asarray(b)
-    if a.dtype.kind in "biufc" and b.dtype.kind in "biufc":
-        return a.shape == b.shape and np.array_equal(a, b, equal_nan=True)
-    return a.dtype == b.dtype and np.array_equal(a, b)  # labels such as the window name
-
-
-class DenseBlockedRouteCPU(unittest.TestCase):
-    def test_same_summary_compares_labels_and_every_number(self):
-        ledger = {"window": "last_base_step", "time_start": 3, "total": {"force": [1.0, -2.0],
-                  "torque": 0.5}, "aperture": {"force": np.array([np.nan, 1.0])}}
-        self.assertTrue(same_summary(ledger, replace_value(ledger, ())))
-        for path, value in ((("window",), "last_substep"), (("time_start",), 4),
-                            (("total", "force"), [1.0, -2.000001]), (("total", "torque"), 0.25),
-                            (("aperture", "force"), np.array([0.0, 1.0]))):
-            with self.subTest(path=path):
-                self.assertFalse(same_summary(ledger, replace_value(ledger, path, value)))
-                self.assertFalse(same_summary(replace_value(ledger, path, value), ledger))
-
-    def test_route_capability_matrix(self):
+    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, cfg=None, bounds=BOUNDS, routes=ROUTES):
-        from CelerisLab.cuda.context import CudaContext
-        cfg = cfg or config(mode)
-        p = FixedPatch2D.from_config(cfg, bounds)
-        arrays = state(NY, NX, 1), state(*p.fine_shape, 2)
-        bodies = sorted({obj.obj_id for obj in objects})
-        results = {}
-        with CudaContext(0) as context:
-            for name, blocks in routes:
-                observe = bool(objects) and blocks[1] is not None
-                snapshots, ledgers = [], []
-                with FixedPatchStepper(cfg, context=context, coarse_boundary=mode, refinements=[bounds],
-                                       level_blocks=blocks, level_objects=[(), tuple(objects)],
-                                       level_jets=[(), tuple(jets)],
-                                       observe_raw_exchange=observe) as backend:
-                    backend.initialize(*arrays)
-                    for step in (0, 1, 2, 8):
-                        if step == 8 and switch is not None:
-                            backend.set_jet_amplitude(0, switch)
-                        backend.step(step - backend.step_count)
-                        snapshots.append(backend.download() + raw(backend))
-                        if observe and step:
-                            ledgers.append(tuple(backend.read_raw_exchange(body, window=window)
-                                                 for body in bodies
-                                                 for window in ("last_base_step", "last_substep")))
-                    backend.initialize(*arrays)
-                    snapshots.append(backend.download() + raw(backend))
-                results[name] = snapshots, ledgers
-        return results
-
-    def check(self, label, results):
-        reference, _ = results["D/D"]
-        for name in ("D/B", "B/B"):
-            for step, actual, expected in zip((0, 1, 2, 8, "reset"), results[name][0], reference):
-                for index, (a, b) in enumerate(zip(actual, expected)):
-                    exact(self, a, b, f"{label} {name} vs D/D step={step} array={index}")
-        if results["D/B"][1]:
-            for step, a, b in zip((1, 2, 8), results["D/B"][1], results["B/B"][1]):
-                self.assertTrue(same_summary(a, b), f"{label} raw exchange D/B vs B/B step={step}")
-        for name, (snapshots, _) in results.items():
-            masses = [float(np.nansum(level.astype(np.float64))) for level in snapshots[3][:2]]
-            print(f"{label} {name}: step 8 mass by level {masses}", flush=True)
-        print(f"PASS {label}: D/D = D/B = B/B bits/NaN masks at 0,1,2,8 + reset", flush=True)
-
-    def test_body_free(self):
-        for mode in ("periodic", "slot_reservoir"):
-            self.check(f"{mode} body-free", self.run_routes(mode))
-
-    def test_passive_circle(self):
-        for mode in ("periodic", "slot_reservoir"):
-            self.check(f"{mode} circle", self.run_routes(mode, objects=[circle()]))
-
-    def test_actuated_circle(self):
-        for mode in ("periodic", "slot_reservoir"):
-            for speed in (SPEED, -SPEED):
-                self.check(f"{mode} actuated {speed}",
-                           self.run_routes(mode, objects=[circle(actuation(speed))]))
-
-    def test_velocity_jet_with_amplitude_switch(self):
-        for mode in ("periodic", "slot_reservoir"):
-            self.check(f"{mode} jet switch",
-                       self.run_routes(mode, objects=[circle()], jets=[jet()], switch=-SPEED))
-
-    def test_mrt_bodies(self):
-        """MRT, periodic root (on_site_velocity slot/reservoir is SRT-only on every route)."""
-        cfg = replace(config("periodic"), collision="MRT")
-        for label, objects, jets, switch in (("circle", [circle()], (), None),
-                                             ("actuated +", [circle(actuation(SPEED))], (), None),
-                                             ("actuated -", [circle(actuation(-SPEED))], (), None),
-                                             ("jet switch", [circle()], [jet()], -SPEED)):
-            self.check(f"MRT periodic {label}",
-                       self.run_routes("periodic", objects, jets, switch, cfg=cfg))
-
-    def test_two_bodies(self):
-        """Two circles on a wider patch (same tiles on every route): Jet + passive, two actuated."""
-        for mode in ("periodic", "slot_reservoir"):
-            for label, objects, jets, switch in TWO_BODY_CASES:
-                self.check(f"{mode} two bodies {label}",
-                           self.run_routes(mode, objects(), jets(), switch, bounds=WIDE_BOUNDS,
-                                           routes=WIDE_ROUTES))
+        # 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
@@ -125,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
@@ -177,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
@@ -187,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()
@@ -196,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
@@ -213,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]
@@ -279,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_gpu.py

@@ -53,24 +53,12 @@ def totals(arrays, maps):
 
 
 def level_storage(level):
-    """Raw device storage of one level as dense (9, ny, nx), dense or blocked."""
-    if level.layout is None:
-        raw = np.empty(level.shape, np.float32)
-        cuda.memcpy_dtoh(raw, level.current)
-        return raw
+    """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)
 
 
-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
-
-
 class CoupledGPU(unittest.TestCase):
     def test_actual_tg_rates_one_interval(self):
         """Algebraic rate/phase check, not a new physical Taylor--Green benchmark.
@@ -181,11 +169,8 @@ class CoupledGPU(unittest.TestCase):
             self.assertFalse(context.is_closed)
         print("partial failure: two fine hops completed, no coarse step published; cleanup PASS", flush=True)
 
-    def eight_steps(self, backend, spec, maps, audit):
-        """Rest and seeded nonequilibrium, 8 steps against the scalar reference.
-
-        audit=True also checks the dense CELERIS_REFINEMENT_AUDIT ownership counters.
-        """
+    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):
@@ -196,11 +181,6 @@ class CoupledGPU(unittest.TestCase):
                                         for a in (state.coarse_f, state.fine_f)))
                 initial = postcollision(state, spec, maps)
                 backend.initialize(*initial)
-                if audit:  # dense audit build only
-                    for level in backend._levels:
-                        for name in ("collision_counts", "explosion_counts", "incoming_counts"):
-                            ptr, size = level.module.get_global(name)
-                            cuda.memset_d8(ptr, 0, size)
                 before = totals(backend.download(), maps)
                 max_error = max_drift = max_moment_error = max_corner = 0.0
                 for step in range(1, 9):
@@ -220,17 +200,7 @@ class CoupledGPU(unittest.TestCase):
                     self.assertLessEqual(drift, MOMENT_ATOL)
                     self.assertLessEqual(moment_error, MOMENT_ATOL)
                     max_drift, max_moment_error = max(max_drift, drift), max(max_moment_error, moment_error)
-                    c, f = backend._levels
-                    if audit:
-                        np.testing.assert_array_equal(global_array(c.module, "collision_counts", (16, 16)),
-                                                      step * maps.coarse_active)
-                        np.testing.assert_array_equal(global_array(f.module, "collision_counts", (24, 24)),
-                                                      2 * step * maps.fine_active)
-                        np.testing.assert_array_equal(global_array(f.module, "explosion_counts", (24, 24)),
-                                                      step * ~maps.fine_active)
-                        owners = np.zeros((9, 16, 16), np.uint32)
-                        owners.ravel()[p.coalesce_coarse] = step
-                        np.testing.assert_array_equal(global_array(c.module, "incoming_counts", owners.shape), owners)
+                    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())
@@ -245,7 +215,7 @@ class CoupledGPU(unittest.TestCase):
                         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)
+                      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."""
@@ -255,7 +225,7 @@ class CoupledGPU(unittest.TestCase):
         with CudaContext(0) as context:
             with PeriodicPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)],
                                       level_blocks=[(16, 16), (24, 24)]) as backend:
-                self.eight_steps(backend, spec, maps, audit=False)
+                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()
@@ -264,59 +234,5 @@ class CoupledGPU(unittest.TestCase):
                         operation()
             self.assertFalse(context.is_closed)
 
-    def test_bounded_equivalence_and_ownership(self):
-        self.assertEqual(Path(CelerisLab.__file__).resolve(), ROOT / "src/CelerisLab/__init__.py")
-        self.assertEqual(Path(compiler._KERNEL_DIR).resolve(), ROOT / "src/CelerisLab/lbm/kernels")
-        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)
-                self.eight_steps(backend, spec, maps, audit=True)
-
-                # Fail before any additional time step; no stale synchronized export.
-                with patch.object(backend._pairs[0], "explode", side_effect=RuntimeError("injected launch failure")):
-                    with self.assertRaisesRegex(RuntimeError, "injected"):
-                        backend.step()
-                with self.assertRaisesRegex(RuntimeError, "failed"):
-                    backend.download()
-                with self.assertRaisesRegex(RuntimeError, "failed"):
-                    backend.step()
-            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

@@ -4,12 +4,12 @@ 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: 1x1/2x2/4x2 fine blocks (blocked (12, 8) root) bit-identical to the
-     FINE_BLOCKS[0] reference (values + NaN masks); dense leg kept until dense retirement;
+     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: blocked layouts
-     bit-identical to the FINE_BLOCKS[0] reference (dense leg until retirement) and the
-     run differs from the unchanged run.
+     bit-identical to the FINE_BLOCKS[0] reference and the run differs from the
+     unchanged run.
 """
 from contextlib import nullcontext
 import unittest
@@ -29,8 +29,8 @@ SPEED = .002
 
 
 def root_and(fine_block):
-    """All-blocked levels with a (12, 8) root tile; fine_block None is the dense route."""
-    return [None if fine_block is None else (12, 8), fine_block]
+    """Blocked levels with a (12, 8) root tile."""
+    return [(12, 8), fine_block]
 
 
 def jet(speed=SPEED):
@@ -107,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)
@@ -123,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)):
@@ -133,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)):
@@ -143,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"):
@@ -158,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(),)])
 
@@ -192,9 +187,6 @@ class JetStepperGPU(unittest.TestCase):
                 public, raw_arrays = run(mode, fine_block, objects=[circle()], jets=[jet()])
                 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}")
-            dense = run(mode, None, objects=[circle()], jets=[jet()])  # dense leg (D2 retires)
-            for a, b in zip(dense[0] + tuple(dense[1]), reference[0] + tuple(reference[1])):
-                np.testing.assert_array_equal(a, b, err_msg=f"G1 {mode} dense")
             legacy = run(mode, FINE_BLOCKS[0], objects=[circle(actuation(SPEED))], jets=())
             diff = np.abs(reference[1][1][:, fluid] - legacy[1][1][:, fluid]).max()
             print(f"G2 {mode}: jet vs legacy fine max |dG| = {diff:.3e}", flush=True)
@@ -205,12 +197,12 @@ class JetStepperGPU(unittest.TestCase):
             reference = run(mode, FINE_BLOCKS[0], objects=[circle()], jets=[jet()], switch=-SPEED)
             steady = run(mode, FINE_BLOCKS[0], objects=[circle()], jets=[jet()])
             self.assertFalse(np.array_equal(reference[1][1], steady[1][1], equal_nan=True), mode)
-            for fine_block in (*FINE_BLOCKS[1:], None):  # None: dense leg (D2 retires)
+            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), 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, {len(FINE_BLOCKS) - 1} layouts + dense "
+            print(f"G3 {mode}: amplitude switch, {len(FINE_BLOCKS) - 1} layouts "
                   f"identical to {FINE_BLOCKS[0]}", flush=True)
 
 

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(),)))):