docs: sync src_coherence with the landed slices

Stage history and verification, historical label on the migration map,
host-preparation/device order, pair accessors, per-level state split into
fixed plan structure / mutable boundary data / runtime _Level, the two
block_collide families, and the deferred items with their evidence needs.
Doc-only.

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

README.md

@@ -117,8 +117,8 @@ with CudaContext(0) as ctx, FixedPatchStepper(cfg, context=ctx, refinements=[(8,
                                               level_objects=[(), (body,)]) as run:
     rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None]
     # (9, ny, nx) postcollision populations per level; non-leaf cells are discarded.
-    run.initialize(np.broadcast_to(rest, (9, *run.patch.coarse_shape)),
-                   np.broadcast_to(rest, (9, *run.patch.fine_shape)))
+    run.initialize(np.broadcast_to(rest, (9, *run.hierarchy.pairs[0].coarse_shape)),
+                   np.broadcast_to(rest, (9, *run.hierarchy.pairs[0].fine_shape)))
     run.step(128)                      # coarse steps; fine level takes two per step
     coarse_g, fine_g = run.download()  # NaN on covered, ghost and solid cells
 ```
@@ -143,7 +143,7 @@ Scope and limits:
   from their union; adjacent, overlapping and concave unions are accepted when
   their actual dependencies close. There is no fixed rectangle-gap rule.
   The fine level is then one level-wide grid `(2*ny, 2*nx)`:
-  `run.patch.fine_shape` gives it, `initialize`/`download` use it, and only the
+  `run.hierarchy.pairs[0].fine_shape` gives it, `initialize`/`download` use it, and only the
   grid-aligned tiles intersecting the union plus its four-cell ghost band are stored.
   Block sizes must divide the level-wide grid; rectangle edges need not align
   with tiles. Order of rectangles does not matter. L1 bodies/Jets (`level_objects[1]`/`level_jets[1]`) work
@@ -159,8 +159,8 @@ Scope and limits:
   with CudaContext(0) as ctx, FixedPatchStepper(
           cfg, context=ctx, refinements=[((4, 4, 12, 12), (16, 4, 24, 12))],
           level_blocks=[(4, 4), (4, 4)]) as run:
-      run.initialize(np.broadcast_to(rest, (9, *run.patch.coarse_shape)),
-                     np.broadcast_to(rest, (9, *run.patch.fine_shape)))  # (9, 32, 64)
+      run.initialize(np.broadcast_to(rest, (9, *run.hierarchy.pairs[0].coarse_shape)),
+                     np.broadcast_to(rest, (9, *run.hierarchy.pairs[0].fine_shape)))  # (9, 32, 64)
       run.step(32)
       coarse_g, fine_g = run.download()  # coarse: NaN on every covered cell;
                                          # fine: NaN on every non-ACTIVE cell
@@ -193,10 +193,10 @@ Scope and limits:
           cfg, context=ctx, refinements=[(4, 4, 12, 12), (8, 8, 16, 16)],
           level_blocks=[(4, 4)] * 3, level_objects=[(), (), (body,)],
           level_jets=[(), (), (jet,)]) as run:
-      shapes = (run.patch.coarse_shape, run.patch.fine_shape, run.child_patch.fine_shape)
+      shapes = (run.hierarchy.pairs[0].coarse_shape, run.hierarchy.pairs[0].fine_shape, run.hierarchy.pairs[1].fine_shape)
       run.initialize(*(np.broadcast_to(rest, (9, *s)) for s in shapes))  # 16², 24², 24²
       run.step(4)
-      run.set_jet_amplitude(0, -.002, level=2)  # default level=1 still addresses L1
+      run.set_jet_amplitude(0, -.002, level=2)  # explicit level selects the driven level
       run.step(4)
       coarse_g, l1_g, l2_g = run.download()  # solids/covered/nonleaves are NaN
   ```

docs/compatibility.md

@@ -3,6 +3,13 @@
 Changes that can make a previously accepted setup or call behave differently.
 Each entry names the commit, the affected entry points and the observable result.
 
+## Redundant periodic stepper alias removed
+
+`PeriodicPatchStepper` is no longer exported from `lbm.refinement_stepper`.
+Use `FixedPatchStepper` with the same arguments, including explicit
+`level_blocks`. The old name was an alias for that exact class; periodic
+boundaries remain the default. Existing imports of the old name must be updated.
+
 ## 87c3a31: cut-link ownership is validated for every body
 
 `validate_cut_link_ownership` (`body/preprocess/flag_overlay.py`) now runs for
@@ -33,3 +40,18 @@ or negative amplitude (the aperture test is geometric). Bodies with neither
 `SurfaceActuation` nor aperture links still raise
 `ValueError("Body has no surface actuation")`. The value is unchanged: the latest
 step's wall/aperture exchange, not total nozzle force.
+
+## 7afb5ad: `Simulation` fails closed after a Jet amplitude device-write failure
+
+`Simulation.set_jet_amplitude` previously only documented that a device error
+leaves the per-link velocity column undefined; the instance stayed usable.
+Now any error after input validation (internal stream creation or the upload)
+marks the Simulation failed. Every later device-touching public method,
+including the `stream` property (guard added in the follow-up repair commit on
+`claude/jet-update-consistency`), raises
+`RuntimeError("Simulation failed (...); close() and recreate it")`. The host
+amplitude keeps its previous value. `close()` still works and stays retryable;
+cached host readers and host staging calls are not blocked. Invalid input
+(`ValueError`/`IndexError`) is unchanged: it raises before any change and the
+instance stays usable. `FixedPatchStepper.set_jet_amplitude` already failed the
+instance (status `failed`) and is unchanged.

docs/development/projected_pair_explosion.md

@@ -4,7 +4,7 @@
 G-to-G explosion reconstruction. The omitted/default `"constant"` selector retains
 the same kernel, launch arguments and device allocations. Projected mode adds two
 compact tables per adjacent pair; recursion, collision, BFL, restriction, buffer
-ownership and synchronized publication are unchanged. `remap_child` preserves it.
+ownership and synchronized publication are unchanged. `remap_finest` preserves it.
 
 Runtime eligibility is periodic-root, fully blocked D2Q9/SRT/FP32, body/Jet-free,
 without LES; the fixed-patch backend has no forcing interface. Shifted and unshifted

docs/development/src_coherence.md

@@ -0,0 +1,113 @@
+# src coherence stage (from main 914906e)
+
+Owner direction 2026-10-28: unified style over historical compatibility, inside `src`.
+No numerical change is intended. Stages after 5f69eba: 65339ab (one boundary/Jet record,
+host preparation helpers, pair accessors), c5a3cdb (`block_collide` families), 442979b
+(exchange metadata on the boundary record). Each was checked by parent/candidate
+config.h/PTX comparison and host/GPU equality runs; c5a3cdb changed the PTX of the two
+channel-root SRT modules only (omega load/clamp moved, outputs exact).
+
+## Conventions (defined once, in `lbm.hierarchy.Hierarchy`)
+
+- Level `l = 0 .. depth-1` (0 = root). Pair `hierarchy.pairs[l-1]` builds level `l`.
+- Array shapes are `(ny, nx)`; rectangles, origins and centers are `(x, y)`.
+- Coarse/fine name the two members (parent/child) of ANY pair, not absolute L0/L1:
+  `fine_shape`, `fine_origin`, `prepare_fine_boundary`, `blocked_fine_tables`,
+  `blocked_coarse_tables`.
+- Refined levels (l >= 1) of FixedPatchStepper place bodies/Jets in the root EDGE frame
+  (root cell `i` = `[i, i+1)`). Uniform `Simulation` objects AND the static circles of a
+  single-level (depth 1) FixedPatchStepper use the integer-NODE frame. Nothing converts
+  between them. This concerns input placement only: within the EDGE frame, refined
+  levels still project root-edge coordinates to fine-local cell indices with the
+  half-cell shift (`prepare_fine_boundary`: `(v - o)/dx - .5`; hits back via `fine_hits_to_edge_frame`
+  `origin + dx*(hit + .5)`), using `edge_frame`/`fine_origin` below.
+- `Hierarchy.edge_frame(l)` -> `((ox, oy), h)`: cell `(i, j)` of level `l`'s allocated
+  array (ghost band and holes included) covers `[o + h*i, o + h*(i+1))`. A scalar
+  FixedPatch2D window carries the half-coarse-cell ghost shift `bounds - HALO/2` in
+  `fine_origin`; a set's level-wide grid has origin 0. Level 1 boundaries are built from
+  the root config (the wall-aware fields live there); deeper ones from the parent frame.
+- `FixedPatch2D.fine_level_offset` (was `fine_level_origin`) is an array INDEX offset
+  `(y, x)`, not an `(x, y)` origin.
+
+## FixedPatchStepper migration map
+
+Historical: left column = names at or after 914906e that no longer exist; right column =
+the current name.
+
+| removed / changed | replacement |
+|---|---|
+| `patch`, `child_patch` | `hierarchy.pairs[0]`, `hierarchy.pairs[1]` |
+| `base_layout`, `fine_layout`, `child_layout` (stepper properties) | `layouts[0]`, `layouts[1]`, `layouts[2]` (`FixedPatchSet2D.fine_layout(block)` is unchanged) |
+| `fine_boundary`, `child_boundary` | `boundaries[1]`, `boundaries[2]` (host geometry; `None` = no bodies/Jets) |
+| `_jets` | `jets[1]` |
+| `_boundaries[l]` (dict from L1, level keys), `_level_boundaries` (b4eacf7), `_device_boundaries` (65339ab) | `_boundaries[l]`: the one combined per-level `_BlockedBoundary` record (list from L0, `None` at L0 unless single-level static circles): host geometry `.host`, `.tables`, `.jets`, `.jet_cache`, exchange `.spacing`/`.centers`, device handles |
+| `_fine_tables`, `_actuated`, `_jet_cache` | `_boundaries[1].tables` / `.actuated` / `.jet_cache` |
+| `_level_jets` (private) | `_boundaries[l].jets` (public `jets[l]` is derived) |
+| `_exchange_spacing[l]`, `_exchange_centers[l]` (private) | `_boundaries[l].spacing` (h_l, root units) / `.centers` (root edge frame); only levels with a boundary report exchange |
+| `_level_blocks` (private) | `_layouts[l].block` |
+| `_remap_cfg` (private) | `_root_cfg` |
+| `refinement_tables.blocked_base_tables` (also via `refinement_stepper`) | `refinement_tables.blocked_coarse_tables` (import from `refinement_tables`; the stepper no longer re-exports) |
+| `_fine_leaf`, `_child_leaf` | `leaf[1]`, `leaf[2]` (copies) |
+| `_middle_kind` | `hierarchy.kinds[1]` |
+| `_boundary_hosts` (private) | public `boundaries[l]` = `_boundaries[l].host` (host geometry only; the private record is the combined one) |
+| `remap_child(rects, child_objects=...)` | `remap_finest(rects, objects=...)` |
+| `set_jet_amplitude(i, a)` (implicit level 1) | `set_jet_amplitude(i, a, level=l)` (required) |
+| `read_raw_exchange(b)` / `read_exchange_window(b)` (implicit level 1) | `read_raw_exchange(b, level=l)` / `read_exchange_window(b, level=l)` (required) |
+| `_level_name(l)` "base/fine/child" messages | `"level l"` |
+| `refinement.MIN_PATCH_GAP` | deleted (no gap policy exists) |
+| `params.cuh HaloPlan`, `collide_trt_no_force` | deleted (zero references) |
+
+Kept: `step_count`, constructor keywords (`refinements`, `level_blocks`, `level_objects`,
+`level_jets`, `coarse_boundary`, `explosion`), `Simulation.set_jet_amplitude(i, a)` (one level),
+`RemapTopology.from_config(..., child_objects=, child_jets=)` (internal remap helper; "child" =
+the replaced level relative to its parent).
+
+## Responsibilities
+
+- Construction: `_check_levels` (inputs); host preparation, numpy only, before any device
+  work: `_level_layouts`, `_kernel_configs` (per-level kernel configs), `_host_boundaries`,
+  `_prepare_blocked_boundary`, `_fluid_leaves`; `Hierarchy.edge_frame` (geometry frames);
+  then `_planned_bytes()` (memory preflight) and `_build_level` / `_build_pair` /
+  `_upload_boundary` (device). A construction failure closes the partial instance.
+  Observation reads share `_observed_boundary` + `_summarize`.
+- Pair accessors: `rectangles` (a scalar FixedPatch2D gives `(bounds,)`) and
+  `fine_level_offset` (`(0, 0)` for a set) on both pair types; the types and their
+  storage frames stay distinct.
+- Per-level state, by lifecycle:
+  - Fixed plan structure (built once in `__init__`, read by `_planned_bytes()` before any
+    allocation, readable after `close()` through the public views): `_layouts`,
+    `_level_cfgs`, `_leaf`, `_state`, and in each `_boundaries[l]` record the host geometry
+    `.host`, link tables, `.jet_cache`, `.spacing`, `.centers`.
+  - Mutable boundary data in the same record: `.jets` and `.tables["surface_velocity"]`,
+    changed only by `set_jet_amplitude` (host prepare -> device write -> commit; an uncertain
+    write fails the instance); device handles (`.flags`, `.device`, `.kernel`, ledgers) are
+    set by `_upload_boundary` and cleared by `close()`.
+  - Runtime: `_levels` (`_Level`) holds compiled modules, kernels, the two population
+    buffers (swapped per step) and transfer/kind handles. It is created only after the
+    memory preflight and cleared on `close()`; `.layout`/`.boundary` reference the plan for
+    `_hop`. Folding the plan into `_Level` would mix both lifecycles in one record.
+- Blocked kernels: `operators/block_collision.cuh block_collide` is the one collision extension
+  point of `BlockStep` and `BlockChannelStep`. It computes rho/u once and reads no runtime
+  force (`d_params.fx/fy`). Two compile-time families: uniform-equivalent (MRT on any level,
+  SRT on the channel root `BLOCK_ROOT_BC == 2`): clamped omega, `collide_mrt_no_force` or
+  `collide_srt` with its zero Guo term; and level SRT (`BLOCK_ROOT_BC` 0/1): unclamped omega,
+  `collide_srt_no_force`, with the shifted direct-momentum feq under `USE_DDF_SHIFTING`.
+  The uniform path keeps `collide_dispatch` (runtime force, LES, TRT).
+- `explosion="projected_linear"` stays an optional capability; `"constant"` stays default.
+
+## Not done / open
+
+- Two per-level config derivations stay distinct on purpose: `Hierarchy.configs`
+  (`_middle_config`, wall fields reset on every refined level, geometry validation) vs.
+  `FixedPatchStepper._kernel_configs` (reset only on level 1 of a slot_reservoir root;
+  generates each level's config.h). Folding them needs a config.h/PTX diff gate.
+- No single-`Simulation` facade over FixedPatchStepper: the two own different state models
+  (uniform node frame, body sync, forcing, checkpoint vs. fixed blocked hierarchy, edge
+  frame, remap-as-new-instance). A facade would only forward calls; no concrete need found.
+- Deferred, each needing its own evidence first: sharing the channel-root cell/donor selector
+  with OneStep (the closure functions are already shared; a shared selector changes uniform
+  signatures and needs a uniform PTX gate); merging FixedPatch2D into FixedPatchSet2D
+  (different origin, ghost band, layout size and addressing; needs a geometry/parity table and
+  would need tiles that are not grid-aligned); one `clamp_omega` helper for uniform and
+  blocked paths (it lives in the hash-pinned LES header `turbulence_smag.cuh`).
+- Tests, examples and validation callers are migrated separately (pinned runners stay frozen).

docs/fixed_three_level_multipatch.md

@@ -61,7 +61,7 @@ Stationary whole-circle L2 bodies and Jets (`level_objects[2]`, `level_jets[2]`)
 coordinates and a single L2 boundary object. Body IDs are physical object-list
 positions, not region IDs or caller `obj_id`. Complete boundary stencils must
 avoid the transfer footprint. Change a child Jet with
-`set_jet_amplitude(index, value, level=2)`; L1 bodies/Jets remain unsupported.
+`set_jet_amplitude(index, value, level=2)` (`level` is required); L1 bodies/Jets remain unsupported.
 
 ## Verification
 

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`. |
-| 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. |
+| 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_finest`), 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`). | 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 device-touching public method (`Simulation`: `stream`, `set_jet_amplitude`, `initialize`, `sync_bodies`, `run`/`step`, `read_surface_exchange`, `read_body(s)`, `get_macroscopic`, `get_ddf`, `update_runtime_params`, `snapshot`, `restore`, save/load checkpoint) raises; cached host readers and host-side staging (e.g. `get_flags`, `add_body`/`remove_body` before a sync) do not touch the device and stay callable. `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`.
@@ -34,7 +35,7 @@ Tests: `tests/unit/test_refinement_mrt.py`, `tests/validation/multigrid/test_ref
 | `slot_reservoir` | on-site velocity bottom slot, reservoir top, velocity 0, 2 <= slot width < nx, patch rows 4..ny-4 | `lbm/refinement_stepper.py` `FixedPatchStepper.__init__`; `lbm/refinement.py` `prepare_fine_boundary` |
 | `channel` | all levels blocked; Zou-He inlet, pressure-NEQ outlet, stationary bounce-back walls; transfers away from the closure | `lbm/channel.py` `validate_channel_config`, `validate_channel_transfers` |
 | Explosion | `constant` (default); `projected_linear` only periodic, all-blocked SRT, no bodies/Jets, no LES (experimental) | `lbm/refinement_stepper.py` `FixedPatchStepper.__init__` |
-| L2 runtime remap (`remap_child`) | periodic three-level L2 sets only, no L0/L1 bodies, identical pinned boundary | `lbm/refinement_stepper.py` `FixedPatchStepper.remap_child` |
+| L2 runtime remap (`remap_finest`) | periodic three-level L2 sets only, no L0/L1 bodies, identical pinned boundary | `lbm/refinement_stepper.py` `FixedPatchStepper.remap_finest` |
 
 Tests: `tests/unit/test_refinement_three_level_multipatch.py`, `tests/validation/multigrid/test_refinement_multipatch.py`,
 `test_refinement_three_level.py`, `test_refinement_three_level_multipatch_gpu.py`, `test_refinement_blocked_base.py`,
@@ -83,6 +84,6 @@ Tests: `tests/unit/test_jets.py`, `test_jets_review.py`, `test_shifted_jets.py`,
 
 ## Not supported in this release
 
-3D, dynamic/automatic AMR (tagging, time-varying hierarchies beyond `remap_child`),
+3D, dynamic/automatic AMR (tagging, time-varying hierarchies beyond `remap_finest`),
 moving bodies on refined levels, pressure Jets, Jet checkpoints, LES or TRT on refined
 levels, partial (non-dividing) tiles.

docs/synchronized_leaf_remap.md

@@ -33,7 +33,7 @@ origin is retained. Physical lower edges are `origin + spacing*(x,y)` on L0,
 `origin + spacing*(parent.fine_origin + (x,y)/4)` on L2.
 
 Factory options `coarse_boundary` and `objects` must declare `"periodic"` and
-an empty collection; L2 bodies use `child_objects` as described below.
+an empty collection; L2 bodies use `objects` as described below.
 LBMConfig's channel settings do not select this operator's boundary mode.
 Geometry metadata is factory-created;
 manual mutation/forging of metadata is unsupported. Caller certifies the
@@ -77,7 +77,7 @@ PYTHONPATH=src python -m unittest discover -s tests/unit -p 'test_synchronized_r
 
 ## Explicit runtime continuation
 
-`replacement = current.remap_child(rectangles)` returns a new ready
+`replacement = current.remap_finest(rectangles)` returns a new ready
 `FixedPatchStepper` at the **same absolute step_count**. It supports only
 periodic three-level runs with no L0/L1 bodies and an existing L2
 `FixedPatchSet2D`; stationary L2 bodies/Jets follow the support contract below.
@@ -99,7 +99,7 @@ resources; an actual old-device synchronization/read failure marks the old
 stepper failed, consistent with its existing failure semantics.
 
 ```python
-replacement = current.remap_child([(20, 16, 36, 40)])
+replacement = current.remap_finest([(20, 16, 36, 40)])
 # Explicitly choose the returned continuation only after success.
 current.close()
 current = replacement
@@ -119,7 +119,7 @@ Per-population conservation is not an advancement target under collision.
 
 ### Bounded pinned-body extension
 
-`replacement = current.remap_child(rectangles, child_objects=original_objects)`
+`replacement = current.remap_finest(rectangles, objects=original_objects)`
 also accepts stationary Circles inside the L2 union. Each body's **physical
 boundary, complete BFL support and current per-link velocity must remain
 identical**. Rectangles are not pinned: they may be split, merged or resized,
@@ -130,7 +130,7 @@ coverage and two-hop BFL transfer-footprint checks still apply. This is manual
 remapping, not adaptive body tracking or runtime body editing.
 
 Supply all original root edge-frame definitions in their construction order.
-Do not pass `current.child_boundary.objects`: those are projected fine-local copies.
+Do not pass `current.boundaries[2].objects`: those are projected fine-local copies.
 There is no retained mutable original-object registry. Before reading device
 state, remapping reprojects the supplied definition onto the old topology and
 compares exact ordered circle centers/radii, flags, physical cut links, body

examples/fixed_patch_circle.py

@@ -1,9 +1,9 @@
 """Public API demo: fixed 2:1 patch with a stationary, optionally actuated circle.
 
 From this checkout, after checking the selected physical GPU is idle:
-    PYTHONPATH="$PWD/src" CUDA_VISIBLE_DEVICES=0 python examples/fixed_patch_circle.py --fine-block 16 20
-    ... --speed -0.002 --fine-block 16 20 --mode slot_reservoir
-Omit --fine-block for the (equivalent) dense fine storage.
+    PYTHONPATH="$PWD/src" CUDA_VISIBLE_DEVICES=0 python examples/fixed_patch_circle.py
+    ... --speed -0.002 --mode slot_reservoir
+Both levels use blocked storage (--base-block, --fine-block tiles).
 
 A short smoke run (default 128 coarse steps): it shows setup, stepping and
 fluid-only diagnostics, not a developed flow, accuracy or long-time stability.
@@ -40,13 +40,13 @@ def config(mode):
 def leaf_masks(run, mode):
     """Known fluid leaves: uncovered coarse cells (minus wall rows in slot mode)
     and fine ACTIVE cells outside the body."""
-    p = run.patch
+    p = run.hierarchy.pairs[0]
     x0, y0, x1, y1 = p.bounds
     coarse = np.ones(p.coarse_shape, bool)
     coarse[y0:y1, x0:x1] = False
     if mode == "slot_reservoir":
         coarse[[0, -1]] = False
-    solid = (run.fine_boundary.flags.reshape(p.fine_shape) & SOLID) != 0
+    solid = (run.boundaries[1].flags.reshape(p.fine_shape) & SOLID) != 0
     return coarse, (p.fine_kind == ACTIVE) & ~solid
 
 
@@ -71,8 +71,10 @@ def main():
     parser.add_argument("--mode", choices=("periodic", "slot_reservoir"), default="periodic")
     parser.add_argument("--speed", type=float, default=.002,
                         help="normal surface velocity, coarse lattice units (+ blows, - sucks, 0 passive)")
-    parser.add_argument("--fine-block", type=int, nargs=2, metavar=("BY", "BX"),
-                        help="blocked fine storage; must divide the whole fine allocation (32 x 40 here)")
+    parser.add_argument("--base-block", type=int, nargs=2, metavar=("BY", "BX"), default=(12, 8),
+                        help="root tile; must divide the 24 x 32 root")
+    parser.add_argument("--fine-block", type=int, nargs=2, metavar=("BY", "BX"), default=(16, 20),
+                        help="fine tile; must divide the whole fine allocation (32 x 40 here)")
     parser.add_argument("--steps", type=int, default=128)
     args = parser.parse_args()
     if args.steps < 0:
@@ -86,20 +88,20 @@ def main():
     body = SimObject(0, geometry=CircleGeometry(*CENTER, RADIUS), center=CENTER, radius=RADIUS,
                      surface_actuation=spec)
     cfg = config(args.mode)
-    fine_block = tuple(args.fine_block) if args.fine_block else None
     # The caller owns the CUDA context; the stepper must be closed first (the
     # nested with-statement does that).
     with CudaContext(0) as context, FixedPatchStepper(
             cfg, context=context, coarse_boundary=args.mode, refinements=[BOUNDS],
-            level_blocks=[None, fine_block], level_objects=[(), (body,)]) as run:
-        p = run.patch
+            level_blocks=[tuple(args.base_block), tuple(args.fine_block)],
+            level_objects=[(), (body,)]) as run:
+        p = run.hierarchy.pairs[0]
         # Rest equilibrium is already postcollision G, so no initial collision is
         # needed. Non-leaf entries (covered, ghost, solid) are discarded by initialize.
         rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None]
         run.initialize(np.broadcast_to(rest, (9, *p.coarse_shape)),
                        np.broadcast_to(rest, (9, *p.fine_shape)))
-        storage = f"blocked {fine_block}" if fine_block else "dense"
-        print(f"mode={args.mode} fine storage={storage} speed={args.speed:+g} "
+        storage = f"blocked root {tuple(args.base_block)} fine {tuple(args.fine_block)}"
+        print(f"mode={args.mode} storage={storage} speed={args.speed:+g} "
               f"(fine allocation {p.fine_shape[0]}x{p.fine_shape[1]}, fine origin {p.fine_origin})")
         masks = leaf_masks(run, args.mode)
         start = leaf_diagnostics(*run.download(), masks)

examples/fixed_patch_jet.py

@@ -31,6 +31,7 @@ from CelerisLab.lbm.refinement import ACTIVE
 from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
 
 BOUNDS, AMPLITUDE, WIDTH = (16, 4, 32, 28), .02, 12
+LEVEL_BLOCKS = [(12, 8), (8, 8)]  # blocked storage: 48x48 root, 56x40 fine window
 TIMES, MAX_MACH, MAX_SECONDS = (0, 300, 1200), .1, 120
 
 
@@ -132,7 +133,8 @@ def main():
                   lattice_model='D2Q9', collision='SRT', streaming='double_buffer',
                   store_precision='FP32', ddf_shifting=False, les_enabled=False)
     report = dict(config=dict(nx=48, ny=48, bounds=BOUNDS, W=WIDTH, A=AMPLITUDE, nu=.1,
-                              rho=1, tau_coarse=.8, tau_fine=1.1, coarse_boundary='slot_reservoir'),
+                              rho=1, tau_coarse=.8, tau_fine=1.1, coarse_boundary='slot_reservoir',
+                              level_blocks=LEVEL_BLOCKS),
                   guards=dict(check_every=50, positive_density=True, finite=True,
                               max_Mach=MAX_MACH, max_seconds=MAX_SECONDS,
                               timing='cooperative checks, not an exact interruption deadline'),
@@ -149,8 +151,9 @@ def main():
     start, run = time.monotonic(), None
     try:
         with CudaContext(0) as context, FixedPatchStepper(cfg, context=context,
-                coarse_boundary='slot_reservoir', refinements=[BOUNDS]) as run:
-            p = run.patch
+                coarse_boundary='slot_reservoir', refinements=[BOUNDS],
+                level_blocks=LEVEL_BLOCKS) as run:
+            p = run.hierarchy.pairs[0]
             masks = [np.ones(p.coarse_shape, bool), p.fine_kind == ACTIVE]
             masks[0][[0, -1]] = False
             masks[0][4:28, 16:32] = False

examples/fixed_periodic_patch.py

@@ -16,7 +16,7 @@ from CelerisLab.config import load_lbm_config
 from CelerisLab.cuda.context import CudaContext
 from CelerisLab.lbm.equilibrium import compute_feq_d2q9, compute_macro_from_ddf
 from CelerisLab.lbm.refinement import ACTIVE
-from CelerisLab.lbm.refinement_stepper import PeriodicPatchStepper
+from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
 
 
 def shear_equilibrium(shape, origin, dx, period):
@@ -48,8 +48,9 @@ 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:
-        p = run.patch
+    with CudaContext(0) as context, FixedPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)],
+                                                           level_blocks=[(16, 16), (24, 24)]) as run:
+        p = run.hierarchy.pairs[0]
         coarse = shear_equilibrium(p.coarse_shape, (0, 0), 1.0, cfg.ny)
         fine = shear_equilibrium(p.fine_shape, p.fine_origin, .5, cfg.ny)
         coarse_leaf = np.ones(p.coarse_shape, dtype=bool)

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/blocks.py

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

src/CelerisLab/lbm/hierarchy.py

@@ -21,7 +21,24 @@ from .refinement import (ACTIVE, OUTER_GHOST, FixedPatch2D, FixedPatchSet2D,
 @dataclass(frozen=True, eq=False)
 class Hierarchy:
     """configs[l]: validated LBMConfig of level l; pairs[l]: tables of (l, l+1);
-    kinds[l]: composed per-cell kind (child coverage OUTER on every non-finest level)."""
+    kinds[l]: composed per-cell kind (child coverage OUTER on every non-finest level).
+
+    Index conventions (the single place they are defined):
+      * level l = 0 (root) .. depth-1; pair index p = l-1 builds level l from level p;
+      * array shapes are (ny, nx); rectangles, origins and centers are (x, y);
+      * a pair's two members are its COARSE (parent) and FINE (child) level, for every
+        pair, not absolute L0/L1 (fine_shape, fine_origin, blocked_coarse_tables, ...);
+      * refined levels (l >= 1) of FixedPatchStepper place bodies and Jets in the root
+        EDGE frame: root cell i covers [i, i+1), center i+.5. Uniform Simulation objects
+        and the static circles of a single-level (depth 1) FixedPatchStepper use the
+        integer-NODE frame (CircleGeometry.build_flag_mask at integer (x, y)). No
+        NODE<->EDGE input conversion is applied anywhere; inside the EDGE frame, refined
+        boundaries still project root-edge coordinates to fine-local cell indices with
+        the half-cell shift ((v - o)/h - .5, see prepare_fine_boundary);
+      * cell (i, j) of level l's allocated array (ghost band and holes included) covers
+        [o + h*i, o + h*(i+1)) with (o, h) = edge_frame(l); level 1 is built from the root
+        config directly, deeper levels from their parent frame (see edge_frame).
+    """
     configs: tuple
     pairs: tuple
     kinds: tuple
@@ -30,6 +47,18 @@ class Hierarchy:
     def depth(self):
         return len(self.kinds)
 
+    def edge_frame(self, level):
+        """(origin (x, y), spacing) of level ``level``'s array in the root edge frame.
+
+        Each pair's fine_origin is in parent cells: zero for a set (level-wide grid),
+        the shifted window origin for a scalar FixedPatch2D.
+        """
+        origin, spacing = (0.0, 0.0), 1.0
+        for pair in self.pairs[:level]:
+            origin = tuple(o + spacing * f for o, f in zip(origin, pair.fine_origin))
+            spacing /= 2
+        return origin, spacing
+
 
 def is_rectangle(value):
     value = tuple(value)
@@ -86,14 +115,14 @@ def _check_chain(coarse, fine, pair):
 def leaf_partition(h):
     """Base-cell fraction published as ACTIVE over all levels (exact: powers of two).
 
-    Window (scalar) levels are offset by their fine_level_origin first.
+    Window (scalar) levels are offset by their fine_level_offset first.
     """
     base = np.zeros(h.kinds[0].shape)
     oy = ox = 0
     for level, kind in enumerate(h.kinds):
         if level:
             pair = h.pairs[level - 1]
-            py, px = (pair.fine_level_origin if isinstance(pair, FixedPatch2D) else (0, 0))
+            py, px = pair.fine_level_offset
             oy, ox = 2 * oy + py, 2 * ox + px
         y, x = np.nonzero(kind == ACTIVE)
         np.add.at(base, ((y + oy) >> level, (x + ox) >> level), 0.25 ** level)

src/CelerisLab/lbm/kernels/blocked.cu

@@ -34,9 +34,7 @@
 #if BLOCK_ROOT_BC == 1
 #include "boundary/block_wall.cuh"
 #endif
-#if BLOCK_ROOT_BC != 2
 #include "operators/block_collision.cuh"
-#endif
 
 // Kernels, in definition order: transfers, boundary staging, then the one full
 // step selected for this module (BlockStep, or BlockChannelStep for a channel root).

src/CelerisLab/lbm/kernels/core/params.cuh

@@ -44,16 +44,4 @@ struct RigidBodyControl2D {
     int   control_mode;           // 0=force, 1=velocity, 2=displacement
 };
 
-// ============================================================================
-// Halo plan (multi-GPU / AMR 预留)
-// ============================================================================
-struct HaloPlan {
-    int       face_id;
-    int       peer_gpu;
-    int       level;
-    size_t    count;
-    int*      send_idx;
-    int*      recv_idx;
-};
-
 #endif // CELERIS_CORE_PARAMS_CUH

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

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

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

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

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

@@ -57,25 +57,4 @@ __device__ __forceinline__ void collide_trt(float* __restrict__ f,
     }
 }
 
-__device__ __forceinline__ void collide_trt_no_force(float* __restrict__ f,
-                                                     const float* __restrict__ feq,
-                                                     float omega)
-{
-    const float wp = omega;
-    const float wm = compute_omega_minus(wp);
-
-    f[0] = f[0] - wp * (f[0] - feq[0]);
-
-    #pragma unroll
-    for (int i = 1; i < NQ; i += 2) {
-        const int ib = i + 1;
-        const float fi = f[i], fb = f[ib];
-        const float feqi = feq[i], feqb = feq[ib];
-        const float delta_s = 0.5f * ((fi - feqi) + (fb - feqb));
-        const float delta_a = 0.5f * ((fi - feqi) - (fb - feqb));
-        f[i]  = fi  - wp * delta_s - wm * delta_a;
-        f[ib] = fb  - wp * delta_s + wm * delta_a;
-    }
-}
-
 #endif // CELERIS_OPERATORS_COLLISION_TRT_CUH

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_channel.cu

@@ -48,18 +48,7 @@ extern "C" __global__ void BlockChannelStep(
         } else {
             apply_wall_bb_d2q9_rows(gy, d_block_channel.ny - 2, f, in + base, k);
         }
-        if (!wall && !outlet) {
-#if COLLISION_MODEL == 2
-            collide_mrt_unforced(f);
-#else
-            // Match uniform SRT arithmetic, including its shifted equilibrium path.
-            float rho, ux, uy, feq[NQ], forcing[NQ] = {};
-            compute_rho_u(f, rho, ux, uy);
-            compute_feq(rho, ux, uy, feq);
-            const float omega = fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega));
-            collide_srt(f, feq, forcing, omega);
-#endif
-        }
+        if (!wall && !outlet) block_collide(f);
     }
     stream_pull_store(k, f, out + base);
 }

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

@@ -115,10 +115,16 @@ class FixedPatch2D:
         return self.bounds[0] - HALO / 2, self.bounds[1] - HALO / 2
 
     @property
-    def fine_level_origin(self) -> tuple[int, int]:
-        """(y, x) of this window's first cell on the level-wide 2:1 fine grid."""
+    def fine_level_offset(self) -> tuple[int, int]:
+        """Array index (y, x) of this window's first cell on the level-wide 2:1 fine grid
+        (an index offset, not an (x, y) origin)."""
         return 2 * self.bounds[1] - HALO, 2 * self.bounds[0] - HALO
 
+    @property
+    def rectangles(self) -> tuple:
+        """The one parent rectangle, as a sequence like FixedPatchSet2D.rectangles."""
+        return (self.bounds,)
+
     @property
     def covered(self) -> np.ndarray:
         """Dense coarse mask of cells replaced by the patch (not leaves)."""
@@ -151,7 +157,6 @@ class FixedPatch2D:
             raise ValueError("levels are not synchronized: fine_steps must equal 2*coarse_steps")
 
 
-MIN_PATCH_GAP = 4  # Legacy exported constant; union transfers impose no gap policy.
 
 
 @dataclass(frozen=True, eq=False)
@@ -225,11 +230,18 @@ class FixedPatchSet2D:
     def bounds(self) -> tuple:
         return tuple(p.bounds for p in self.patches)
 
+    rectangles = bounds
+
     @property
     def fine_origin(self) -> tuple[float, float]:
         """Lower edge of the level-wide fine grid in coarse physical units."""
         return 0.0, 0.0
 
+    @property
+    def fine_level_offset(self) -> tuple[int, int]:
+        """The level-wide grid itself: zero (y, x) index offset."""
+        return 0, 0
+
     @property
     def covered(self) -> np.ndarray:
         mask = np.zeros(self.coarse_shape, bool)

src/CelerisLab/lbm/refinement_stepper.py

@@ -1,10 +1,11 @@
 # CelerisLab/lbm/refinement_stepper.py
-"""Opt-in fixed patch backend with optional stationary fine-interior circles.
+"""Blocked fixed 2:1 refinement hierarchy (any depth) with optional stationary circles/Jets.
 
 The caller owns the CudaContext and must close this backend before that context.
 Only synchronized postcollision_g leaf populations are public. No checkpoint,
-forcing or runtime parameter mutation is offered. Bounded manual periodic
-L2 set replacement is available at completed base intervals via remap_child.
+forcing or runtime parameter mutation besides Jet amplitudes is offered. A bounded
+manual periodic level-2 set replacement is available at completed base intervals
+via remap_finest. Level/pair/shape/frame conventions: lbm.hierarchy.Hierarchy.
 """
 from dataclasses import dataclass, replace
 from numbers import Integral
@@ -19,9 +20,9 @@ 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)
+from .refinement_tables import (blocked_coarse_tables, blocked_fine_tables,
                                 blocked_kind, projected_explosion_tables, _prepare_blocked_boundary,
                                 _BlockedBoundary,
                                 BLOCKED_BOUNDARY_DEVICE_COLUMNS)
@@ -62,8 +63,26 @@ class _Level:
     boundary: object = None
 
 
-def _level_name(level):
-    return ("base", "fine", "child")[level] if level < 3 else f"level {level}"
+def _kernel_configs(cfg, pairs, wall_aware):
+    """Per-level kernel LBMConfig (window shape, viscosity doubling per level).
+
+    Deliberately distinct from Hierarchy.configs (refinement._middle_config), which
+    resets the wall fields on EVERY refined level for geometry validation. Here only
+    level 1 of a slot_reservoir root resets them (internal bounce-back/velocity ghost
+    setup); other levels inherit their parent's fields. These configs generate each
+    level's compiled config.h, so folding the two derivations is a numerical-parity
+    change (config.h/PTX diff gate), not a style one.
+    """
+    configs = [replace(cfg)]
+    for index, pair in enumerate(pairs):
+        parent = configs[-1]
+        level = replace(parent, nx=pair.fine_shape[1], ny=pair.fine_shape[0],
+                        viscosity=2 * parent.viscosity)
+        if index == 0 and wall_aware:
+            level = replace(level, y_wall_bc="bounce_back", top_wall_bc="velocity",
+                            wall_slot_width=0, wall_normal_velocity=0.0, velocity=0.0)
+        configs.append(level)
+    return configs
 
 
 MAX_GRID_DIM = 65535  # CUDA gridDim.y/z limit (rows, blocks)
@@ -75,25 +94,74 @@ 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")
 
 
+# Host preparation (numpy only, before any device work); each returns one entry per level.
+
+def _level_layouts(cfg, pairs, level_blocks, wall_aware):
+    """Blocked storage layout of every level, checked against the launch limits."""
+    layouts = []
+    for level, block in enumerate(level_blocks):
+        if level == 0:
+            layout = BlockLayout2D.build(cfg.ny, cfg.nx, *block, periodic=(True, not wall_aware))
+        elif isinstance(pairs[level - 1], FixedPatchSet2D):
+            layout = pairs[level - 1].fine_layout(block)
+        else:
+            layout = BlockLayout2D.build(*pairs[level - 1].fine_shape, *block,
+                                         periodic=(False, False))
+        if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM:
+            raise ValueError(f"level {level} blocks exceed launch limits")
+        layouts.append(layout)
+    return layouts
+
+
+def _host_boundaries(hierarchy, level_cfgs, level_objects, level_jets, wall_aware, shifted):
+    """Host boundary geometry of every refined level in the root edge frame (None at L0)."""
+    hosts = [None]
+    for level in range(1, hierarchy.depth):
+        pair = hierarchy.pairs[level - 1]
+        if level == 1:  # the root config carries the wall-aware boundary fields
+            host = prepare_fine_boundary(level_cfgs[0], pair, level_objects[1],
+                                         wall_aware=wall_aware, jets=level_jets[1],
+                                         _shifted_blocked=shifted)
+        else:
+            origin, spacing = hierarchy.edge_frame(level - 1)
+            host = prepare_fine_boundary(level_cfgs[level - 1], pair, level_objects[level],
+                                         jets=level_jets[level], parent_origin=origin,
+                                         parent_spacing=spacing, _shifted_blocked=shifted)
+        hosts.append(host)
+    return hosts
+
+
+def _fluid_leaves(hierarchy, boundaries):
+    """ACTIVE cells of every level that no static circle or body makes solid."""
+    leaves = [kind == ACTIVE for kind in hierarchy.kinds]
+    for level, b in enumerate(boundaries):
+        if b is None:
+            continue
+        if b.host is None:  # single level: static root circles
+            leaves[level] &= ~b.tables["solid"]
+        else:
+            fine_shape = hierarchy.pairs[level - 1].fine_shape
+            leaves[level] &= (b.host.flags.reshape(fine_shape) & SOLID) == 0
+    return leaves
+
+
 class FixedPatchStepper:
     """D2Q9 SRT/MRT FP32 fixed 2:1 refinement hierarchy, built from level lists.
 
@@ -101,15 +169,17 @@ 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).
+    passive circles (static BFL) in the integer-NODE frame, like Simulation
+    objects (see Hierarchy for frames).
 
     ``coarse_boundary`` is ``periodic`` (default) or ``slot_reservoir``:
     a bottom on-site slot and top fixed-density reservoir from ``cfg``.
-    ``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
@@ -117,7 +187,7 @@ class FixedPatchStepper:
     Periodic mode does not apply the wall/slot boundary fields in ``cfg``.
     Outer walls remain coarse; fine halo donors must be ordinary fluid.
 
-    Bodies are fixed circles in root (L0) physical coordinates; velocities are
+    Bodies of refined levels are fixed circles in root (L0) physical coordinates; velocities are
     in root lattice units (dx=dt=1). Independent local copies are projected at
     fine cell centers; their complete stencil must be active and disjoint from
     the actual two-hop transfer footprint. These use the edge/cell-center frame
@@ -132,14 +202,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 +217,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.
@@ -166,11 +234,18 @@ class FixedPatchStepper:
     lattice velocity); every whole circle and its BFL stencil must fit strictly
     away from that level's transfer footprint.
 
-    ``explosion="projected_linear"`` selects experimental transverse linear
-    reconstruction on each adjacent pair. Default "constant" is unchanged.
-    Projected mode requires periodic root, all-blocked SRT levels without
-    bodies/Jets or LES; this backend supplies no forcing. Geometry is generic,
-    but union/multilevel GPU and matched-physics qualification remain pending.
+    ``explosion`` selects the parent->child reconstruction of every pair:
+    ``"constant"`` (default) or ``"projected_linear"``, an optional transverse
+    linear reconstruction (Schornbaum-Ruede/Chen style) that is GPU-confirmed
+    more accurate than constant only on periodic two-level single-rectangle
+    Taylor-Green cases; union/multilevel GPU and matched-physics qualification
+    remain pending. It requires a periodic root and SRT levels without
+    bodies/Jets or LES; this backend supplies no forcing.
+
+    Per-level read-only state: ``hierarchy.pairs`` (pair l-1 builds level l),
+    ``layouts``, ``boundaries`` (host geometry), ``jets`` and ``leaf`` masks,
+    each indexed by level. ``set_jet_amplitude``, ``read_raw_exchange`` and
+    ``read_exchange_window`` take an explicit ``level=``.
 
     One recursive hop advances every level: a child explodes from its parent's
     current hop source, takes two hops, then the parent gathers halos, stages
@@ -180,30 +255,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):
@@ -220,63 +297,17 @@ class FixedPatchStepper:
         self.hierarchy = build_hierarchy(cfg, refinements, wall_aware=self._wall_aware)
         pairs = self.hierarchy.pairs
         if self._wall_aware:
-            rects = (pairs[0].bounds,) if isinstance(pairs[0], FixedPatch2D) else pairs[0].bounds
-            if any(r[1] < 4 or r[3] > cfg.ny - 4 for r in rects):
+            if any(r[1] < 4 or r[3] > cfg.ny - 4 for r in pairs[0].rectangles):
                 raise ValueError("slot_reservoir requires y0 >= 4 and y1 <= ny-4")
-        self._layouts = []
-        for level, block in enumerate(level_blocks):
-            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")
-            self._layouts.append(layout)
+        self._layouts = _level_layouts(cfg, pairs, level_blocks, self._wall_aware)
         if self._channel:
-            base = blocked_base_tables(pairs[0], self._layouts[0], np.argsort(pairs[0].coalesce_coarse))
+            base = blocked_coarse_tables(pairs[0], self._layouts[0], np.argsort(pairs[0].coalesce_coarse))
             validate_channel_transfers(pairs[0], self._layouts[0], base[1])
-        # Kernel configurations per level (geometry configs live in self.hierarchy).
-        self._level_cfgs = [replace(cfg)]
-        if pairs:
-            fine_cfg = replace(cfg, nx=pairs[0].fine_shape[1], ny=pairs[0].fine_shape[0],
-                               viscosity=2 * cfg.viscosity)
-            if self._wall_aware:
-                fine_cfg = replace(fine_cfg, y_wall_bc="bounce_back", top_wall_bc="velocity",
-                                   wall_slot_width=0, wall_normal_velocity=0.0, velocity=0.0)
-            self._level_cfgs.append(fine_cfg)
-        for pair in pairs[1:]:
-            parent = self._level_cfgs[-1]
-            self._level_cfgs.append(replace(parent, nx=pair.fine_shape[1], ny=pair.fine_shape[0],
-                                            viscosity=2 * parent.viscosity))
-        # Per-level boundary geometry in the root edge frame (level l spacing 2**-l).
-        self._boundary_hosts = [None]
-        origin, spacing = (0.0, 0.0), 1.0
-        spacings = [spacing]
-        for level in range(1, depth):
-            pair = pairs[level - 1]
-            if level == 1:
-                host = prepare_fine_boundary(cfg, pair, level_objects[1],
-                                             wall_aware=self._wall_aware, jets=level_jets[1],
-                                             _shifted_blocked=self._ddf_shifting)
-            else:
-                host = prepare_fine_boundary(self._level_cfgs[level - 1], pair, level_objects[level],
-                                             jets=level_jets[level], parent_origin=origin,
-                                             parent_spacing=spacing, _shifted_blocked=self._ddf_shifting)
-            self._boundary_hosts.append(host)
-            origin = tuple(o + spacing * f for o, f in zip(origin, pair.fine_origin))
-            spacing /= 2
-            spacings.append(spacing)
+        self._level_cfgs = _kernel_configs(cfg, pairs, self._wall_aware)
+        hosts = _host_boundaries(self.hierarchy, self._level_cfgs, level_objects, level_jets,
+                                 self._wall_aware, self._ddf_shifting)
         if not isinstance(observe_raw_exchange, bool):
             raise ValueError("observe_raw_exchange must be bool")
-        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")
@@ -286,32 +317,25 @@ class FixedPatchStepper:
         self._window_valid = False
         self._window_start = None
         self._exchange_valid = False
-        self._exchange_spacing = tuple(spacings)
-        self._exchange_centers = tuple(tuple(tuple(o.geometry.center) for o in objects)
-                                       for objects in level_objects)
-        self._level_jets = level_jets
         root = None
         if depth == 1 and level_objects[0]:
             tables = static_circle_tables(self._layouts[0], level_objects[0])
             root = _BlockedBoundary(None, tables, (), False)
-        self._level_boundaries = [root] + [
-            _prepare_blocked_boundary(self._boundary_hosts[level], self._layouts[level],
+        # The one boundary/Jet record per level; host geometry and Jets are views of it.
+        self._boundaries = [root] + [
+            _prepare_blocked_boundary(hosts[level], self._layouts[level],
                                       self._level_cfgs[level], level_jets[level],
                                       _shifted_blocked=self._ddf_shifting)
             for level in range(1, depth)]
-        self._leaf = [kind == ACTIVE for kind in self.hierarchy.kinds]
-        if root is not None:
-            self._leaf[0] &= ~root.tables["solid"]
-        for level in range(1, depth):
-            host = self._boundary_hosts[level]
-            if host is not None:
-                self._leaf[level] &= (host.flags.reshape(pairs[level - 1].fine_shape) & SOLID) == 0
+        for level, b in enumerate(self._boundaries):
+            if b is not None:  # exchange scaling (dt_l = h_l) and root-edge-frame moment centers
+                b.spacing = self.hierarchy.edge_frame(level)[1]
+                b.centers = tuple(tuple(o.geometry.center) for o in level_objects[level])
+        self._leaf = _fluid_leaves(self.hierarchy, self._boundaries)
         self._state = [mask.copy() for mask in self._leaf]
         if self._channel:
             self._leaf[0] &= channel_masks((cfg.ny, cfg.nx))[0]
-        # Private constructor inputs for the bounded manual L2 continuation.
-        self._remap_cfg = replace(cfg)
-        self._level_blocks = level_blocks
+        self._root_cfg = replace(cfg)  # private constructor input for remap_finest
         self.ctx = context
         self._lock = RLock()
         self._status = "new"
@@ -319,8 +343,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:
@@ -333,38 +355,8 @@ class FixedPatchStepper:
                 for index, pair in enumerate(pairs):
                     self._pairs.append(self._build_pair(index, pair))
                 for level in self._levels:
-                    b = level.boundary
-                    if b is not None:
-                        b.flags = self._upload(b.tables["flags"])
-                        b.device = {k: self._upload(v) for k, v in b.tables.items()
-                                    if k in BLOCKED_BOUNDARY_DEVICE_COLUMNS}
-                        if self._observe_raw_exchange:
-                            b.device.update(rx=self._upload(b.tables["rx"]),
-                                            ry=self._upload(b.tables["ry"]))
-                            count = len(b.tables["cell"])
-                            b.last_raw = self._allocate(count * 4 * 4)
-                            b.interval_integral = self._allocate(count * 4 * 8)
-                            if self._observe_window:
-                                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")
+                    if level.boundary is not None:
+                        self._upload_boundary(level)
             except BaseException:
                 try:
                     self.close()
@@ -372,6 +364,22 @@ class FixedPatchStepper:
                     pass  # Preserve the original construction failure.
                 raise
 
+    def _upload_boundary(self, level):
+        """Device columns, BFL kernel and (when observed) the exchange ledgers of one level."""
+        b = level.boundary
+        b.flags = self._upload(b.tables["flags"])
+        b.device = {k: self._upload(v) for k, v in b.tables.items()
+                    if k in BLOCKED_BOUNDARY_DEVICE_COLUMNS}
+        if self._observe_raw_exchange:
+            b.device.update(rx=self._upload(b.tables["rx"]), ry=self._upload(b.tables["ry"]))
+            count = len(b.tables["cell"])
+            b.last_raw = self._allocate(count * 4 * 4)
+            b.interval_integral = self._allocate(count * 4 * 8)
+            if self._observe_window:
+                b.window_integral = self._allocate(count * 4 * 8)
+            b.accumulate = level.module.get_function("BlockAccumulateExchange")
+        b.kernel = level.module.get_function("BlockCurvedBoundary")
+
     def _planned_bytes(self):
         """Both G buffers, halo pairs and boundary columns of every level (transfer tables excluded).
 
@@ -380,12 +388,9 @@ 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
-            b = self._level_boundaries[index]
+            total += 2 * 4 * int(np.prod(layout.blocked_shape))
+            total += layout.halo_dst.nbytes + layout.halo_src.nbytes
+            b = self._boundaries[index]
             if b is not None:
                 total += sum(np.asarray(v).nbytes for k, v in b.tables.items()
                              if k == "flags" or k in BLOCKED_BOUNDARY_DEVICE_COLUMNS)
@@ -397,83 +402,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.boundary = self._boundaries[index]
+        level.halo_copy = module.get_function("BlockHaloCopy")
+        if len(layout.halo_dst):
+            level.halo_dst = self._upload(layout.halo_dst)
+            level.halo_src = self._upload(layout.halo_src)
+        if index == 0 and self._channel:
+            cfg = self._root_cfg
+            ptr, size = module.get_global("d_block_channel")
+            if size != struct.calcsize("=IIf"):
+                raise RuntimeError("Block channel host/device layout mismatch")
+            cuda.memcpy_htod(ptr, struct.pack("=IIf", cfg.nx, cfg.ny, cfg.velocity))
+        if index == 0 and self._wall_aware:
+            cfg = self._root_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_coarse_tables(pair, parent.layout, order,
+                                                               self._wall_aware)
+        else:
+            # Nested pairs keep inherited INNER/OUTER and sparse holes of the parent.
+            donors, targets, parent_kind = blocked_coarse_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 +490,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 +534,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
@@ -570,14 +549,15 @@ class FixedPatchStepper:
                 self._status = "failed"
                 raise
 
-    def remap_child(self, rectangles, *, child_objects=()):
-        """Return a NEW ready stepper with an L2 set replaced at this base time.
+    def remap_finest(self, rectangles, *, objects=()):
+        """Return a NEW ready stepper with the level-2 set replaced at this base time.
 
-        Periodic three-level set children only; L0/L1 and block sizes stay fixed.
-        A body run requires its original root-edge-frame ``child_objects``:
+        Periodic three-level hierarchies whose level 2 is a set only; levels 0/1 and
+        block sizes stay fixed. A body run requires its original root-edge-frame
+        level-2 ``objects``:
         stationary circles with identical physical boundary and BFL support.
         Other fluid coverage may change, even in the same connected union.
-        Current L2 Jets (including updated amplitudes)
+        Current level-2 Jets (including updated amplitudes)
         are preserved. Local projected boundary objects are not valid inputs.
         The lock and stream synchronization certify
         current postcollision G and the completed 1:2:4 counts internally.
@@ -595,19 +575,18 @@ class FixedPatchStepper:
             pairs = self.hierarchy.pairs
             if (self._wall_aware or self._channel or self.hierarchy.depth != 3
                     or not isinstance(pairs[1], FixedPatchSet2D)
-                    or self._boundary_hosts[1] is not None
-                    or any(b is not None for b in self._level_boundaries[:2])):
-                raise ValueError("remap_child requires periodic three-level L2 sets with no L0/L1 bodies")
-            cfg = self._remap_cfg
-            child_objects = tuple(child_objects)
-            actual = self._level_boundaries[2]
+                    or any(b is not None for b in self._boundaries[:2])):
+                raise ValueError("remap_finest requires a periodic three-level hierarchy with a "
+                                 "level-2 set and no level-0/1 bodies")
+            cfg = self._root_cfg
+            objects = tuple(objects)
+            actual = self._boundaries[2]
             jets = () if actual is None else actual.jets
             old = RemapTopology.from_config(cfg, pairs[0], pairs[1].bounds,
-                                           child_objects=child_objects, child_jets=jets,
+                                           child_objects=objects, child_jets=jets,
                                            _shifted_blocked=self._ddf_shifting)
-            if (not same_physical_boundary(old.boundary, self._boundary_hosts[2])
-                    or (self._boundary_hosts[2] is not None and actual is None)):
-                raise ValueError("child_objects must reproduce the original root-edge-frame boundary")
+            if not same_physical_boundary(old.boundary, None if actual is None else actual.host):
+                raise ValueError("objects must reproduce the original root-edge-frame level-2 boundary")
             if actual is not None:
                 velocity = actual.tables.get("surface_velocity")
                 if velocity is None:
@@ -615,18 +594,17 @@ class FixedPatchStepper:
                                                  _shifted_blocked=self._ddf_shifting)
                 if (not same_physical_boundary(old.boundary, actual.host)
                         or not np.array_equal(old.surface_velocity, velocity)):
-                    raise ValueError("child_objects changed boundary ownership or surface velocity")
+                    raise ValueError("objects changed boundary ownership or surface velocity")
             new = RemapTopology.from_config(cfg, pairs[0], rectangles,
-                                           child_objects=child_objects, child_jets=jets,
+                                           child_objects=objects, child_jets=jets,
                                            _shifted_blocked=self._ddf_shifting)
             if (not same_physical_boundary(old.boundary, new.boundary)
                     or not np.array_equal(old.surface_velocity, new.surface_velocity)):
                 raise ValueError("remap requires the identical pinned physical boundary and BFL support")
-            base_block, fine_block, child_block = self._level_blocks
             # Check sparse block alignment/launch size before any device read/allocation.
-            layout = new.child.fine_layout(child_block)
+            layout = new.child.fine_layout(self._layouts[2].block)
             if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM:
-                raise ValueError("child blocks exceed launch limits")
+                raise ValueError("level 2 blocks exceed launch limits")
             with self.ctx.activate():
                 try:
                     self._stream.synchronize()
@@ -642,8 +620,8 @@ class FixedPatchStepper:
                 phase="postcollision_g", **options)
             candidate = type(self)(
                 cfg, context=self.ctx, refinements=[pairs[0].bounds, new.child.bounds],
-                level_blocks=[base_block, fine_block, child_block],
-                level_objects=[(), (), child_objects], level_jets=[(), (), jets],
+                level_blocks=[layout.block for layout in self._layouts],
+                level_objects=[(), (), objects], level_jets=[(), (), jets],
                 observe_raw_exchange=self._observe_raw_exchange,
                 observe_window_exchange=self._observe_window, explosion=self._explosion)
             try:
@@ -660,30 +638,9 @@ 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):
+        """Replace Jet ``index`` of level ``level`` (``jets[level]``) between steps; no table rebuild.
 
-    def set_jet_amplitude(self, index, amplitude, *, level=1):
-        """Replace one jet's amplitude between steps; no geometry/table rebuild.
-
-        ``level=1`` preserves the L1 default; use ``level=2`` for child_jets.
         Invalid input raises before any change. A device/sync error marks the
         instance failed (as step/initialize do); host state commits only after
         a successful upload.
@@ -692,30 +649,25 @@ class FixedPatchStepper:
             self._require("new", "ready")
             if (isinstance(level, bool) or not isinstance(level, Integral)
                     or not 1 <= level < self.hierarchy.depth):
-                raise ValueError("Jet level must name an existing refined level (1 = first)")
-            boundary = self._level_boundaries[level]
-            jets = list(self._level_jets[level])
+                raise ValueError("Jet level must name an existing refined level (1 .. depth-1)")
+            boundary = self._boundaries[level]
+            # Host preparation: validate and build the new column without touching state.
+            jets = list(() if boundary is None else boundary.jets)
             jets[index] = jets[index].with_amplitude(amplitude)  # validates first
-            velocity = None
-            if boundary is not None and boundary.jet_cache is not None:
-                shape, jet_id = boundary.jet_cache
-                velocity = boundary.tables["surface_velocity"].copy()
-                driven = jet_id >= 0
-                velocity[driven] = jet_link_velocity(jets, shape, jet_id)[driven]
+            shape, jet_id = boundary.jet_cache
+            velocity = boundary.tables["surface_velocity"].copy()
+            driven = jet_id >= 0
+            velocity[driven] = jet_link_velocity(jets, shape, jet_id)[driven]
+            # Device action: an uncertain write fails the instance closed.
             try:
-                if level == 1 and self._fine_curved is not None:
-                    self._fine_curved.set_jet_amplitude(index, amplitude, self._stream)
-                elif velocity is not None:
-                    self._stream.synchronize()  # no in-flight step reads the column
-                    cuda.memcpy_htod(boundary.device["surface_velocity"], velocity)
+                self._stream.synchronize()  # no in-flight step reads the column
+                cuda.memcpy_htod(boundary.device["surface_velocity"], velocity)
             except BaseException:
                 self._status = "failed"
                 raise
-            if velocity is not None:
-                boundary.tables["surface_velocity"] = velocity
-            if boundary is not None:
-                boundary.jets = tuple(jets)
-            self._level_jets[level] = tuple(jets)
+            # Commit: the one boundary record holds both the column and the Jets.
+            boundary.tables["surface_velocity"] = velocity
+            boundary.jets = tuple(jets)
 
     def _blocked_boundary(self, level, buffer, **kw):
         b = level.boundary
@@ -747,19 +699,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:
@@ -771,7 +710,7 @@ class FixedPatchStepper:
         b = level.boundary
         if b is not None and b.last_raw is not None:
             count = len(b.tables["cell"])
-            h = self._exchange_spacing[index]  # acoustic scaling: dt_l = h_l
+            h = b.spacing  # acoustic scaling: dt_l = h_l
             b.accumulate(b.last_raw, b.interval_integral, np.uint32(count),
                          np.float64(h), np.float64(h),
                          grid=((count + self._block[0] - 1) // self._block[0], 1), **kw)
@@ -815,7 +754,7 @@ class FixedPatchStepper:
                 self._exchange_valid = False
                 for _ in range(n):
                     if self._observe_raw_exchange:  # stream-ordered: holds the last interval only
-                        for b in self._level_boundaries:
+                        for b in self._boundaries:
                             if b is not None:
                                 cuda.memset_d8_async(b.interval_integral, 0,
                                                      len(b.tables["cell"]) * 32, self._stream)
@@ -838,7 +777,7 @@ class FixedPatchStepper:
                 raise RuntimeError("window exchange observation is disabled")
             self._window_valid = False
             try:
-                for b in self._level_boundaries:
+                for b in self._boundaries:
                     if b is not None and b.window_integral is not None:
                         cuda.memset_d8(b.window_integral, 0, len(b.tables["cell"]) * 32)
             except BaseException:
@@ -847,7 +786,23 @@ class FixedPatchStepper:
             self._window_start = self._step_count
             self._window_valid = True
 
-    def read_exchange_window(self, body_id, *, level=1):
+    def _observed_boundary(self, level, body_id):
+        if (isinstance(level, bool) or not isinstance(level, Integral)
+                or not 0 <= level < len(self._levels)):
+            raise ValueError("level must name an existing level")
+        b = self._boundaries[level]
+        if (b is None or isinstance(body_id, bool) or not isinstance(body_id, Integral)
+                or not 0 <= body_id < len(b.host.objects)):
+            raise ValueError("body_id must name a body on this level")
+        return b
+
+    def _summarize(self, b, rows, level, body_id, window, **kw):
+        return summarize_exchange(rows, b.tables["body_id"], b.tables["aperture"],
+                                  body_id=body_id, level=level, h=b.spacing,
+                                  time_end=self._step_count, window=window,
+                                  reference_center=b.centers[body_id], **kw)
+
+    def read_exchange_window(self, body_id, *, level):
         """Non-destructive read of [window start, current step count] in base-lattice units."""
         with self._lock, self.ctx.activate():
             self._require("ready")
@@ -858,26 +813,16 @@ class FixedPatchStepper:
             duration = self._step_count - self._window_start
             if duration < 1:
                 raise RuntimeError("exchange window requires at least one completed base interval")
-            if (isinstance(level, bool) or not isinstance(level, Integral)
-                    or not 0 <= level < len(self._levels)):
-                raise ValueError("level must name an existing level")
-            b = self._level_boundaries[level]
-            if (b is None or isinstance(body_id, bool) or not isinstance(body_id, Integral)
-                    or not 0 <= body_id < len(b.host.objects)):
-                raise ValueError("body_id must name a body on this level")
+            b = self._observed_boundary(level, body_id)
             rows = np.empty((len(b.tables["cell"]), 4), np.float64)
             try:
                 cuda.memcpy_dtoh(rows, b.window_integral)
             except BaseException:
                 self._status = "failed"
                 raise
-            return summarize_exchange(rows, b.tables["body_id"], b.tables["aperture"],
-                                      body_id=body_id, level=level, h=self._exchange_spacing[level],
-                                      time_end=self._step_count, window="window",
-                                      reference_center=self._exchange_centers[level][body_id],
-                                      duration=duration)
+            return self._summarize(b, rows, level, body_id, "window", duration=duration)
 
-    def read_raw_exchange(self, body_id, *, level=1, window="last_base_step"):
+    def read_raw_exchange(self, body_id, *, level, window="last_base_step"):
         """Read raw boundary reaction about its geometric center in the root edge frame.
 
         Repeated reads and step(0) preserve the window. last_substep samples
@@ -891,13 +836,7 @@ class FixedPatchStepper:
                 raise RuntimeError("raw exchange observation is disabled")
             if not self._exchange_valid:
                 raise RuntimeError("raw exchange requires a completed base interval")
-            if (isinstance(level, bool) or not isinstance(level, Integral)
-                    or not 0 <= level < len(self._levels)):
-                raise ValueError("level must name an existing level")
-            b = self._level_boundaries[level]
-            if (b is None or isinstance(body_id, bool) or not isinstance(body_id, Integral)
-                    or not 0 <= body_id < len(b.host.objects)):
-                raise ValueError("body_id must name a body on this level")
+            b = self._observed_boundary(level, body_id)
             if window not in ("last_base_step", "last_substep"):
                 raise ValueError("window must be last_base_step or last_substep")
             count = len(b.tables["cell"])
@@ -907,10 +846,7 @@ class FixedPatchStepper:
             except BaseException:
                 self._status = "failed"
                 raise
-            return summarize_exchange(rows, b.tables["body_id"], b.tables["aperture"],
-                                      body_id=body_id, level=level, h=self._exchange_spacing[level],
-                                      time_end=self._step_count, window=window,
-                                      reference_center=self._exchange_centers[level][body_id])
+            return self._summarize(b, rows, level, body_id, window)
 
     def download(self):
         """Return physical G (FP64 for shifted storage), NaN on non-leaves."""
@@ -924,7 +860,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 +868,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)
@@ -951,21 +887,13 @@ class FixedPatchStepper:
                         self._stream.synchronize()
                 except Exception as exc:
                     error = exc
-                for b in getattr(self, "_level_boundaries", ()):
+                for b in getattr(self, "_boundaries", ()):
                     if b is not None:
                         b.flags = b.device = b.kernel = None
                         b.last_raw = b.interval_integral = b.accumulate = None
                         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()
@@ -986,28 +914,21 @@ class FixedPatchStepper:
             if exc[0] is None:
                 raise
 
-    # Read-only compatibility views over the level lists, kept on purpose for the existing
-    # tests and docs. Internal code must not use them; it reads hierarchy/_layouts/_level_*.
-    patch = property(lambda self: self.hierarchy.pairs[0])
-    child_patch = property(lambda self: self.hierarchy.pairs[1] if self.hierarchy.depth > 2 else None)
-    base_layout = property(lambda self: self._layouts[0])
-    fine_layout = property(lambda self: self._layouts[1])
-    child_layout = property(lambda self: self._layouts[2] if len(self._layouts) > 2 else None)
-    fine_boundary = property(lambda self: self._boundary_hosts[1])
-    child_boundary = property(lambda self: self._boundary_hosts[2] if len(self._boundary_hosts) > 2
-                              else None)
-    _jets = property(lambda self: self._level_jets[1])
-    _boundaries = property(lambda self: {i: b for i, b in enumerate(self._level_boundaries) if i})
-    _fine_tables = property(lambda self: None if self._level_boundaries[1] is None
-                            else self._level_boundaries[1].tables)
-    _actuated = property(lambda self: self._level_boundaries[1] is not None
-                         and self._level_boundaries[1].actuated)
-    _jet_cache = property(lambda self: None if self._level_boundaries[1] is None
-                          else self._level_boundaries[1].jet_cache)
-    _fine_leaf = property(lambda self: self._leaf[1])
-    _child_leaf = property(lambda self: self._leaf[2] if len(self._leaf) > 2 else None)
-    _middle_kind = property(lambda self: self.hierarchy.kinds[1] if self.hierarchy.depth > 2 else None)
+    # Read-only per-level state, indexed by level (0 = root). Pairs: hierarchy.pairs[l-1].
+    @property
+    def layouts(self):
+        return tuple(self._layouts)
 
+    @property
+    def boundaries(self):
+        """Host boundary geometry per level (None: no bodies/Jets on a refined level)."""
+        return tuple(None if b is None else b.host for b in self._boundaries)
 
-# Preserve the original entry point and its periodic default.
-PeriodicPatchStepper = FixedPatchStepper
+    @property
+    def jets(self):
+        return tuple(() if b is None else b.jets for b in self._boundaries)
+
+    @property
+    def leaf(self):
+        """Per-level (ny, nx) masks of public (non-NaN) download cells."""
+        return tuple(m.copy() for m in self._leaf)

src/CelerisLab/lbm/refinement_tables.py

@@ -36,7 +36,7 @@ def blocked_kind(layout, dense_kind):
     kind[owner[owned]] = np.asarray(dense_kind).ravel()[owned]
     return kind
 
-def blocked_base_tables(patch, layout, order, wall_aware=False, *, composed_kind=None):
+def blocked_coarse_tables(patch, layout, order, wall_aware=False, *, composed_kind=None):
     """Parent donor cells and receiver-local arrival slots, in transfer row order.
 
     Stage arrivals only after child reads and parent halo gather, immediately
@@ -130,6 +130,8 @@ class _BlockedBoundary:
     interval_integral: object = None
     accumulate: object = None
     window_integral: object = None  # optional observation-window FP64 ledger
+    spacing: float = None  # level spacing h_l in root units (exchange scaling)
+    centers: tuple = ()    # body reference centers in the root edge frame
 
 
 def _prepare_blocked_boundary(boundary, layout, cfg, jets, *, _shifted_blocked=False):

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)
@@ -89,6 +90,7 @@ class Simulation:
     @context_bound
     def stream(self) -> cuda.Stream:
         """Internal CUDA stream for async DTOH coordination."""
+        self._require_usable()
         if self._stream is None:
             self._stream = cuda.Stream()
         return self._stream
@@ -172,15 +174,27 @@ 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
+        (stream creation or upload) leaves the host amplitude unchanged but the
+        device column undefined, so the Simulation is marked failed: every later
+        device-touching public method (incl. ``stream``) raises; close() still
+        works. 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 +352,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 +380,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 +400,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 +463,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 +496,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 +592,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 +657,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 +676,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 +706,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/cpp/channel_root.cpp

@@ -42,6 +42,7 @@ inline float __fdiv_rn(float a,float b){return a/b;}
 #include "step/one_step_double.cu"
 #else
 #include "core/cell_kind.cuh"
+#include "operators/block_collision.cuh"
 #include "step/block_channel.cu"
 #endif
 int main(){

tests/unit/test_channel_root.py

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

tests/unit/test_example_tiles.py

@@ -0,0 +1,70 @@
+"""The public examples run all-blocked: their level tiles must reach the device unchanged.
+
+Host only: a context whose activation raises proves construction (layouts, tables,
+boundary checks) succeeded up to the first CUDA call.
+"""
+import ast
+from dataclasses import replace
+from pathlib import Path
+import unittest
+from unittest.mock import Mock
+
+from CelerisLab.body.geometry.circle import CircleGeometry
+from CelerisLab.body.objects import SimObject
+from CelerisLab.config import LBMConfig
+from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+
+EXAMPLES = Path(__file__).resolve().parents[2] / "examples"
+
+
+def literal(name, constant):
+    """Module-level literal from an example, read without importing it (no pycuda)."""
+    tree = ast.parse((EXAMPLES / name).read_text())
+    for node in tree.body:
+        if isinstance(node, ast.Assign):
+            for target in node.targets:
+                names = target.elts if isinstance(target, ast.Tuple) else [target]
+                for i, t in enumerate(names):
+                    if isinstance(t, ast.Name) and t.id == constant:
+                        value = node.value.elts[i] if isinstance(target, ast.Tuple) else node.value
+                        return ast.literal_eval(value)
+    raise KeyError(constant)
+
+
+def forbidden():
+    return Mock(activate=Mock(side_effect=AssertionError("CUDA reached")), sm_arch="sm_70")
+
+
+class ExampleTiles(unittest.TestCase):
+    def test_fixed_patch_jet_tiles(self):
+        bounds, blocks = literal("fixed_patch_jet.py", "BOUNDS"), literal("fixed_patch_jet.py", "LEVEL_BLOCKS")
+        width, amplitude = literal("fixed_patch_jet.py", "WIDTH"), literal("fixed_patch_jet.py", "AMPLITUDE")
+        cfg = LBMConfig(nx=48, ny=48, viscosity=.1, velocity=0., wall_slot_width=width,
+                        wall_normal_velocity=amplitude, y_wall_bc="on_site_velocity",
+                        top_wall_bc="reservoir")
+        with self.assertRaisesRegex(AssertionError, "CUDA reached"):
+            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")
+        center, radius = literal("fixed_patch_circle.py", "CENTER"), literal("fixed_patch_circle.py", "RADIUS")
+        body = SimObject(0, geometry=CircleGeometry(*center, radius), center=center, radius=radius)
+        periodic = LBMConfig(nx=nx, ny=ny, viscosity=.1)
+        slot = replace(periodic, velocity=0, wall_normal_velocity=.001, wall_slot_width=8,
+                       y_wall_bc="on_site_velocity", top_wall_bc="reservoir")
+        for mode, cfg in (("periodic", periodic), ("slot_reservoir", slot)):
+            with self.subTest(mode=mode), self.assertRaisesRegex(AssertionError, "CUDA reached"):
+                FixedPatchStepper(cfg, context=forbidden(), coarse_boundary=mode,
+                                  refinements=[bounds], level_blocks=[(12, 8), (16, 20)],
+                                  level_objects=[(), (body,)])
+
+
+if __name__ == "__main__":
+    unittest.main()

tests/unit/test_jet_update_failure.py

@@ -0,0 +1,149 @@
+"""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 import simulation as module
+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_stream_creation_failure_fails_closed(monkeypatch):
+    sim, jet = simulation(monkeypatch)
+    created = []
+    def stream():
+        created.append(1)
+        raise RuntimeError("stream creation failed")
+    monkeypatch.setattr(module.cuda, "Stream", stream, raising=False)
+    sim._stream = None
+    with pytest.raises(RuntimeError, match="stream creation failed"):
+        sim.set_jet_amplitude(0, 2e-3)
+    assert sim.bodies.jets[0] is jet and sim._failed is not None
+    assert sim.field.curved.writes == [] and created == [1]
+    for call in (lambda: sim.stream, lambda: sim.set_jet_amplitude(0, 2e-3)):
+        with pytest.raises(RuntimeError, match="Simulation failed"):
+            call()
+    assert created == [1]  # the public stream property does not retry creation
+    sim.ctx.detach_failure = True
+    with pytest.raises(RuntimeError, match="detach failed"):
+        sim.close()
+    assert not sim._closed
+    sim.close()  # close stays retryable after the failure
+    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._boundaries[2]
+    return boundary.tables["surface_velocity"].copy(), boundary.jets, run.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._boundaries[2].tables["surface_velocity"], velocity)
+        assert (run._boundaries[2].jets, run.jets[2]) == (jets, level_jets)
+        assert run._status == "ready"
+        run.set_jet_amplitude(0, .017, level=2)
+        assert run.jets[2][0].amplitude == .017
+    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._boundaries[2].tables["surface_velocity"], velocity)
+        assert (run._boundaries[2].jets, run.jets[2]) == (jets, level_jets)
+        assert run._status == "failed"
+        module.cuda.memcpy_htod.side_effect = None
+        run._stream.synchronize.side_effect = None
+        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_multiple_pinned_remap.py

@@ -94,7 +94,7 @@ class MultiplePinnedRemap(unittest.TestCase):
         run.set_jet_amplitude(0, .017, level=2)
         run.set_jet_amplitude(1, -.011, level=2)
         before = run.download()
-        clone = run.remap_child((PINS+MOVED+ADDED)[::-1], child_objects=bodies)
+        clone = run.remap_finest((PINS+MOVED+ADDED)[::-1], objects=bodies)
         try:
             self.assertEqual(tuple(j.amplitude for j in clone._boundaries[2].jets), (.017, -.011))
             np.testing.assert_array_equal(clone._boundaries[2].tables['surface_velocity'],
@@ -116,7 +116,7 @@ class MultiplePinnedRemap(unittest.TestCase):
         start = len(allocations)
         module.cuda.memcpy_dtoh.reset_mock()
         with self.assertRaisesRegex(ValueError, 'boundary'):
-            run.remap_child(PINS+FLUID, child_objects=changed)
+            run.remap_finest(PINS+FLUID, objects=changed)
         self.assertEqual(len(allocations), start)
         module.cuda.memcpy_dtoh.assert_not_called()
         with self.assertRaisesRegex(ValueError, 'pinned'):
@@ -128,7 +128,7 @@ class MultiplePinnedRemap(unittest.TestCase):
         cases = [(PINS+FLUID, bodies[:1]), (PINS+FLUID, bodies[::-1]),
                  (PINS+FLUID, (bodies[0], body(center=(32.1, 10.)))),
                  (PINS+FLUID, (bodies[0], body(SurfaceActuation(0., .7, 0.), (32., 10.)))),
-                 (PINS+FLUID, run.child_boundary.objects)]
+                 (PINS+FLUID, run.boundaries[2].objects)]
         for index in range(2):
             changed = list(PINS+FLUID)
             x0, y0, x1, y1 = changed[index]
@@ -141,7 +141,7 @@ class MultiplePinnedRemap(unittest.TestCase):
         run._stream.synchronize.reset_mock()
         for rects, supplied in cases:
             with self.subTest(rects=rects), self.assertRaises(ValueError):
-                run.remap_child(rects, child_objects=supplied)
+                run.remap_finest(rects, objects=supplied)
         self.assertEqual(len(allocations), start)
         module.cuda.memcpy_dtoh.assert_not_called()
         run._stream.synchronize.assert_not_called()

tests/unit/test_pinned_body_remap.py

@@ -61,11 +61,11 @@ class PinnedBodyRemap(unittest.TestCase):
         run,module,allocations,storage,objects=ready()
         run.set_jet_amplitude(0,.017,level=2)
         before=run.download()
-        clone=run.remap_child(B,child_objects=objects)
+        clone=run.remap_finest(B,objects=objects)
         self.assertEqual(clone._boundaries[2].jets[0].amplitude,.017)
         np.testing.assert_array_equal(clone._boundaries[2].tables['surface_velocity'],
                                       run._boundaries[2].tables['surface_velocity'])
-        self.assertTrue(same_physical_boundary(run.child_boundary,clone.child_boundary))
+        self.assertTrue(same_physical_boundary(run.boundaries[2],clone.boundaries[2]))
         clone.step()
         self.assertEqual((run.step_count,clone.step_count),(7,8))
         for got,want in zip(run.download(),before): np.testing.assert_array_equal(got,want)
@@ -73,7 +73,7 @@ class PinnedBodyRemap(unittest.TestCase):
 
     def test_static_body_without_actuation(self):
         run,module,allocations,storage,objects=ready(jets=())
-        clone=run.remap_child(A,child_objects=objects)
+        clone=run.remap_finest(A,objects=objects)
         self.assertFalse(clone._boundaries[2].actuated)
         self.assertNotIn('surface_velocity',clone._boundaries[2].tables)
         for a,b in zip(run.download(),clone.download()): np.testing.assert_array_equal(a,b)
@@ -82,13 +82,13 @@ class PinnedBodyRemap(unittest.TestCase):
     def test_reject_before_device_reads_or_allocation(self):
         run,module,allocations,storage,objects=ready()
         start=len(allocations)
-        cases=[(B,()),(B,run.child_boundary.objects),(B,(body(center=(12.1,14.)),)),
+        cases=[(B,()),(B,run.boundaries[2].objects),(B,(body(center=(12.1,14.)),)),
                (((24,16,36,40),B[1]),objects),(((16,28,32,40),B[1]),objects),
                ((B[1],),objects)]
         module.cuda.memcpy_dtoh.reset_mock()
         for rects,objs in cases:
             with self.subTest(rects=rects),self.assertRaises(ValueError):
-                run.remap_child(rects,child_objects=objs)
+                run.remap_finest(rects,objects=objs)
         self.assertEqual(len(allocations),start)
         module.cuda.memcpy_dtoh.assert_not_called()
         self.assertEqual((run._status,run.step_count),('ready',7))
@@ -100,9 +100,9 @@ class PinnedBodyRemap(unittest.TestCase):
                 (body(SurfaceActuation(0.,.7,speed)),),())
             for spec in (SurfaceActuation(0.,.7,.003),SurfaceActuation(2.,.7,speed)):
                 module.cuda.memcpy_dtoh.reset_mock()
-                with self.assertRaises(ValueError): run.remap_child(B,child_objects=(body(spec),))
+                with self.assertRaises(ValueError): run.remap_finest(B,objects=(body(spec),))
                 module.cuda.memcpy_dtoh.assert_not_called()
-            clone=run.remap_child(B,child_objects=objects)
+            clone=run.remap_finest(B,objects=objects)
             clone.close();run.close()
 
     def test_allocation_failure_keeps_old_ready_and_count(self):
@@ -117,7 +117,7 @@ class PinnedBodyRemap(unittest.TestCase):
             return allocate(size)
         with patch.object(module.cuda,'mem_alloc',side_effect=fail):
             with self.assertRaisesRegex(RuntimeError,'allocation'):
-                run.remap_child(B,child_objects=objects)
+                run.remap_finest(B,objects=objects)
         self.assertEqual((run._status,run.step_count),('ready',7))
         for got,want in zip(run.download(),before): np.testing.assert_array_equal(got,want)
         for a in saved: a.free.assert_not_called()

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()
 
@@ -166,7 +164,7 @@ class ProjectedPairDispatch(unittest.TestCase):
             run, host, _, _, storage = backend(refinements=FOUR[:2], explosion=explosion)
             run._step_count = 7
             old = run.download()
-            clone = run.remap_child(FOUR[1])
+            clone = run.remap_finest(FOUR[1])
             self.assertEqual(clone._explosion, explosion)
             self.assertEqual(clone.step_count, 7)
             for pair in clone._pairs:
@@ -184,7 +182,7 @@ class ProjectedPairDispatch(unittest.TestCase):
                 refinements=FOUR[:2], explosion='projected_linear')
             run._step_count = 3
             before = run.download()
-            clone = run.remap_child(changed)
+            clone = run.remap_finest(changed)
             self.assertEqual(clone._explosion, 'projected_linear')
             self.assertEqual(clone.hierarchy.pairs[1].bounds, tuple(changed))
             self.assertNotEqual(clone._pairs[1].n_ghosts, run._pairs[1].n_ghosts)

tests/unit/test_raw_exchange.py

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

@@ -524,6 +524,11 @@ def test_public_jet_reversal_boundary_and_windows():
         assert r.jet_interval(handle, "reverse", switch + 80, 80, diameter) == (80, None)
         assert r.jet_interval(handle, "fixed", switch, 80, diameter) == (80, None)
         assert not calls
+    refined_calls = []
+    refined = SimpleNamespace(set_jet_amplitude=lambda i, a, *, level: refined_calls.append((i, a, level)))
+    n, event = r.jet_interval(refined, "reverse", switch, 80, diameter, level=1)
+    assert n == 80 and event["step"] == switch and event["tstar"] == 300
+    assert refined_calls == [(0, -.001, 1), (1, .001, 1)]
     rows = [dict(tstar=t, tstar_force=t - .5, cd=2.0, cl=-.1) for t in range(251, 306)]
     windows = r.jet_diagnostics(rows, "reverse")
     assert len(windows) == 6 and all(w["complete"] for w in windows)

tests/unit/test_refinement.py

@@ -177,7 +177,7 @@ class FixedPatchSetTests(unittest.TestCase):
         self.assertEqual((fy, fx), (32, 48))
         rows = cells = 0
         for p in s.patches:
-            y0, x0 = p.fine_level_origin
+            y0, x0 = p.fine_level_offset
             ly, lx = p.fine_shape
             np.testing.assert_array_equal(s.fine_kind[y0:y0 + ly, x0:x0 + lx], p.fine_kind)
             gy, gx = np.divmod(p.ghost_fine.astype(np.int64), lx)
@@ -224,7 +224,7 @@ class FixedPatchSetTests(unittest.TestCase):
             layout = s.fine_layout((4, 4))
             support = np.zeros(s.fine_shape, bool)
             for p in s.patches:
-                y0, x0 = p.fine_level_origin
+                y0, x0 = p.fine_level_offset
                 support[y0:y0 + p.fine_shape[0], x0:x0 + p.fine_shape[1]] = True
             np.testing.assert_array_equal(layout.owned, support)
             self.assertTrue(layout.owned.ravel()[s.ghost_fine.astype(np.intp)].all())

tests/unit/test_refinement_mrt.py

@@ -1,5 +1,6 @@
 """Host guards and source contract for MRT reuse in refined/blocked stepping (no CUDA)."""
 from pathlib import Path
+import re
 import unittest
 
 from CelerisLab.config import LBMConfig
@@ -27,27 +28,30 @@ 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):
+        self.assertIn("collide_srt(f, feq, forcing, omega);", block_collision)
+        self.assertIn("#if COLLISION_MODEL == 2 || BLOCK_ROOT_BC == 2", block_collision)
+        self.assertEqual(channel.count("block_collide(f);"), 1)
+        # Executable text only: the header comment explains the absent force read.
+        code = re.sub(r"//[^\n]*", "", channel + "\n" + block_collision)
+        self.assertIn("d_params.omega", code)
+        for text in ("collide_mrt_unforced", "collide_dispatch(", "d_params.fx", "d_params.fy"):
+            self.assertNotIn(text, code)
+        for source in (blocks,):
             self.assertIn("(COLLISION_MODEL != 0 && COLLISION_MODEL != 2)", source)
             self.assertIn("USE_LES", source)
             self.assertIn("STORE_PRECISION != 0", source)
             self.assertIn('#if COLLISION_MODEL == 2\n#include "operators/collision_mrt.cuh"\n#endif', source)
-        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):
-            self.assertEqual(source.count("collide_mrt_unforced(f);"), 1)
-            self.assertNotIn("collide_mrt_no_force(", source)
-            self.assertNotIn("compute_rho_u(f, rho, ux, uy);\n            collide_mrt", source)
-            self.assertNotRegex(source, r"compute_rho_u\(f, rho, ux, uy\);\s*collide_mrt")
+        # One unforced MRT call, in block_collide only; the steps never collide directly.
+        self.assertEqual(block_collision.count("collide_mrt_no_force(f, rho, ux, uy, omega);"), 1)
+        self.assertEqual(block_collision.count("compute_rho_u(f, rho, ux, uy);"), 1)
+        self.assertEqual(block_collision.count("fmaxf(OMEGA_COLLISION_MIN, d_params.omega)"), 1)
+        for source in (blocks, channel):
+            self.assertNotIn("collide_mrt", source)
+            self.assertNotIn("collide_srt", source)
 
 
 if __name__ == "__main__": unittest.main()

tests/unit/test_refinement_three_level_multipatch.py

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

tests/unit/test_runtime_child_remap.py

@@ -1,5 +1,6 @@
 """CPU mocks execute production remap lifecycle; no GPU numerical claim."""
 from dataclasses import replace
+from types import SimpleNamespace
 import unittest
 from unittest.mock import patch
 
@@ -29,7 +30,7 @@ class RuntimeChildRemap(unittest.TestCase):
         saved = {key: value.copy() for key, value in storage.items()}
         old_allocations = list(allocations)
         module.cuda.memcpy_dtoh.reset_mock()
-        clone = run.remap_child(A)
+        clone = run.remap_finest(A)
         self.assertEqual(module.cuda.memcpy_dtoh.call_count, 3)
         for a, b in zip(clone.download(), before):
             np.testing.assert_array_equal(a, b)
@@ -48,7 +49,7 @@ class RuntimeChildRemap(unittest.TestCase):
         run, module, allocations, storage = ready_backend()
         before = run.download()
         expected = remap(topology(A), topology(B), before).g
-        clone = run.remap_child(B)
+        clone = run.remap_finest(B)
         for level, want, got in zip(clone._levels, expected, clone.download()):
             np.testing.assert_array_equal(got, want)
             for ptr in (level.current, level.spare):
@@ -65,17 +66,18 @@ class RuntimeChildRemap(unittest.TestCase):
         module.cuda.memcpy_dtoh.reset_mock()
         for status in ('new', 'failed', 'closed', 'advancing'):
             run._status = status
-            with self.assertRaises(RuntimeError): run.remap_child(B)
+            with self.assertRaises(RuntimeError): run.remap_finest(B)
         run._status = 'ready'
         for bad in ([], [(40,16,48,40)], [(16,16,32,24),(16,27,32,55)]):
-            with self.assertRaises(ValueError): run.remap_child(bad)
-        hosts, pairs = run._boundary_hosts, run.hierarchy.pairs
+            with self.assertRaises(ValueError): run.remap_finest(bad)
+        records, pairs = run._boundaries, run.hierarchy.pairs
         for attribute, value in (('_wall_aware', True),
-                                 ('_boundary_hosts', [*hosts[:2], object()]),        # child boundary
-                                 ('_boundary_hosts', [hosts[0], object(), hosts[2]]),  # fine boundary
+                                 # level-2 record whose host geometry differs from objects=()
+                                 ('_boundaries', [*records[:2], SimpleNamespace(host=object(), jets=())]),
+                                 ('_boundaries', [records[0], object(), records[2]]),  # level-1 body
                                  ('hierarchy', replace(run.hierarchy, pairs=(pairs[0], None)))):
             with patch.object(run, attribute, value):
-                with self.assertRaises(ValueError): run.remap_child(B)
+                with self.assertRaises(ValueError): run.remap_finest(B)
         self.assertEqual(len(allocations), start)
         module.cuda.memcpy_dtoh.assert_not_called()
         run.close()
@@ -102,7 +104,7 @@ class RuntimeChildRemap(unittest.TestCase):
                         initialize(candidate, *arrays)
                 operation = patch.object(type(run), 'initialize', fail)
             with operation, self.assertRaisesRegex(RuntimeError, 'candidate'):
-                run.remap_child(B)
+                run.remap_finest(B)
             self.assertEqual((run._status, run.step_count), ('ready', 7))
             for a, b in zip(run.download(), before): np.testing.assert_array_equal(a, b)
             for allocation in old_allocations: allocation.free.assert_not_called()
@@ -116,7 +118,7 @@ class RuntimeChildRemap(unittest.TestCase):
             target = run._stream.synchronize if operation == 'sync' else module.cuda.memcpy_dtoh
             target.side_effect = RuntimeError('old device failure')
             with self.assertRaisesRegex(RuntimeError, 'old device'):
-                run.remap_child(B)
+                run.remap_finest(B)
             self.assertEqual(run._status, 'failed')
             self.assertEqual(len(allocations), start)
             target.side_effect = None

tests/unit/test_shifted_jets.py

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

tests/unit/test_shifted_storage.py

@@ -95,7 +95,7 @@ class ShiftedStorage(unittest.TestCase):
             before = run._download_raw()
             expected = remap(topology(A),topology(B),before.values)
             with patch.object(run,'download',side_effect=AssertionError('public decode used')):
-                clone = run.remap_child(B)
+                clone = run.remap_finest(B)
             self.assertEqual((run.step_count,clone.step_count),(7,7))
             self.assertFalse(clone._exchange_valid)
             for level, want, got in zip(clone._levels,expected.g,clone._download_raw().values):
@@ -114,7 +114,7 @@ class ShiftedStorage(unittest.TestCase):
             physical = tuple(np.broadcast_to(_WEIGHTS, level.shape).copy() for level in run._levels)
             for a in physical: a[1] += 2**-35
             run.initialize(*physical)
-            for a,mask in zip(run._download_raw().values,run._leaf):
+            for a,mask in zip(run._download_raw().values,run.leaf):
                 self.assertTrue((a[1,mask] == np.float32(2**-35)).all())
             self.assertEqual(run.step_count,0)
             module.cuda.memcpy_htod.reset_mock()
@@ -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_single_level.py

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

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

@@ -35,10 +35,10 @@ class UnionBodySupport(unittest.TestCase):
         run,module,allocations,storage,objects=ready()
         try:
             run.set_jet_amplitude(0,.017,level=2)
-            clone=run.remap_child(((16,16,24,40),(24,16,32,40),A[1]),child_objects=objects)
+            clone=run.remap_finest(((16,16,24,40),(24,16,32,40),A[1]),objects=objects)
             try:
                 self.assertEqual(clone._boundaries[2].jets[0].amplitude,.017)
-                self.assertTrue(same_physical_boundary(run.child_boundary,clone.child_boundary))
+                self.assertTrue(same_physical_boundary(run.boundaries[2],clone.boundaries[2]))
                 for a,b in zip(run.download(),clone.download()):np.testing.assert_array_equal(a,b)
             finally:clone.close()
         finally:run.close()

tests/unit/test_union_refinement.py

@@ -40,7 +40,7 @@ class Equivalence(unittest.TestCase):
                          table_sets(FixedPatchSet2D.from_config(CFG, [(8, 8, 20, 22)])))
 
     def test_separated_rectangles_reproduce_current_tables(self):
-        rects = [(4, 4, 10, 10), (16, 16, 26, 26)]  # gap 6 >= MIN_PATCH_GAP
+        rects = [(4, 4, 10, 10), (16, 16, 26, 26)]  # separated rectangles; union dependencies determine legality
         self.assertEqual(table_sets(union(*rects)),
                          table_sets(FixedPatchSet2D.from_config(CFG, rects)))
 

tests/validation/multigrid/test_block_stepper.py

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

tests/validation/multigrid/test_channel_root_gpu.py

@@ -113,7 +113,7 @@ class ChannelRootGPU(unittest.TestCase):
                             arrays = [state(*level.shape[1:]) for level in backend._levels]
                             backend.initialize(*arrays)
                             backend.step(4)
-                            if bodies: backend.set_jet_amplitude(0, -.001)
+                            if bodies: backend.set_jet_amplitude(0, -.001, level=1)
                             backend.step(4)
                             raw = backend._download_raw()
                             physical = backend.download()
@@ -128,7 +128,7 @@ class ChannelRootGPU(unittest.TestCase):
                             else:
                                 for a,b in zip(reference, raw.values): np.testing.assert_array_equal(a,b)
                             with self.assertRaisesRegex(ValueError, "periodic"):
-                                backend.remap_child([])
+                                backend.remap_finest([])
 
 
 if __name__ == "__main__": unittest.main()

tests/validation/multigrid/test_exchange_window.py

@@ -59,44 +59,44 @@ class ExchangeWindowCPU(unittest.TestCase):
             stack.enter_context(patch.object(cuda, "memcpy_dtoh",
                                              side_effect=lambda out, src: out.fill(0)))
             backend = build(m)
-            b = backend._level_boundaries[1]
+            b = backend._boundaries[1]
             self.assertIsNotNone(b.window_integral)
             with self.assertRaises(RuntimeError):
                 backend.start_exchange_window()  # not initialized
             arrays = [np.full(level.shape, 1 / 9, np.float32) for level in backend._levels]
             backend.initialize(*arrays)
             with self.assertRaisesRegex(RuntimeError, "no valid exchange window"):
-                backend.read_exchange_window(0)
+                backend.read_exchange_window(0, level=1)
             backend.step(1)
             self.assertEqual(len(accumulations(m)), 2)  # closed window: interval ledger only
             backend.start_exchange_window()
             with self.assertRaisesRegex(RuntimeError, "at least one completed"):
-                backend.read_exchange_window(0)
+                backend.read_exchange_window(0, level=1)
             m["calls"].clear()
             backend.step(2)
             calls = accumulations(m)
             self.assertEqual(len(calls), 8)  # two fine hops x two ledgers x two base steps
             self.assertEqual({id(c[1][1]) for c in calls},
                              {id(b.interval_integral), id(b.window_integral)})
-            out = backend.read_exchange_window(0)
+            out = backend.read_exchange_window(0, level=1)
             self.assertEqual((out["time_start"], out["time_end"], out["sample_count"]), (1., 3, 4))
-            self.assertEqual(backend.read_exchange_window(0)["time_end"], 3)  # non-destructive
+            self.assertEqual(backend.read_exchange_window(0, level=1)["time_end"], 3)  # non-destructive
             backend.initialize(*arrays)  # invalidates
             with self.assertRaisesRegex(RuntimeError, "no valid exchange window"):
-                backend.read_exchange_window(0)
+                backend.read_exchange_window(0, level=1)
             backend.start_exchange_window()
             with patch.object(backend, "_hop", side_effect=RuntimeError("injected")):
                 with self.assertRaises(RuntimeError):
                     backend.step(1)
             self.assertFalse(backend._window_valid)
             with self.assertRaises(RuntimeError):
-                backend.read_exchange_window(0)  # failed status
+                backend.read_exchange_window(0, level=1)  # failed status
 
     def test_disabled_default_has_no_window_cost(self):
         with mocked_backend() as m, ExitStack() as stack:
             stack.enter_context(patch.object(cuda, "memset_d8_async"))
             backend = build(m, window=False)
-            self.assertIsNone(backend._level_boundaries[1].window_integral)
+            self.assertIsNone(backend._boundaries[1].window_integral)
             backend.initialize(*[np.full(l.shape, 1 / 9, np.float32) for l in backend._levels])
             backend.step(1)
             self.assertEqual(len(accumulations(m)), 2)
@@ -109,7 +109,7 @@ class ExchangeWindowGPU(unittest.TestCase):
         from test_raw_exchange_gpu import BOUNDS as RB, CFG, circle as rc, jet, state
         from CelerisLab.cuda.context import CudaContext
         with CudaContext(0) as ctx:
-            kw = dict(context=ctx, refinements=[RB], level_blocks=[None, (16, 20)],
+            kw = dict(context=ctx, refinements=[RB], level_blocks=[(12, 8), (16, 20)],
                       level_objects=[(), (rc(),)], level_jets=[(), (jet(),)],
                       observe_raw_exchange=True)
             with FixedPatchStepper(CFG, observe_window_exchange=True, **kw) as windowed, \
@@ -120,7 +120,7 @@ class ExchangeWindowGPU(unittest.TestCase):
                 windowed.step(2)
                 plain.step(2)
                 windowed.start_exchange_window()
-                b = windowed._level_boundaries[1]
+                b = windowed._boundaries[1]
                 count = len(b.tables["cell"])
                 factors = np.array([.25, .25, .125, .25])  # h*dt, h*h*dt, h*h for level-1 h=.5
                 scale = np.zeros(4)
@@ -140,12 +140,12 @@ class ExchangeWindowGPU(unittest.TestCase):
                 summed = {part: [0., 0., 0., 0.] for part in ("total", "wall", "aperture")}
                 steps = 0
                 for speed in (0., .002, .002, -.001, 0., 0.):  # amplitude updates stay in-window
-                    windowed.set_jet_amplitude(0, speed)
-                    plain.set_jet_amplitude(0, speed)
+                    windowed.set_jet_amplitude(0, speed, level=1)
+                    plain.set_jet_amplitude(0, speed, level=1)
                     windowed.step(1)
                     plain.step(1)
                     steps += 1
-                    obs = windowed.read_raw_exchange(0)
+                    obs = windowed.read_raw_exchange(0, level=1)
                     for part in summed:
                         o = obs[part]
                         summed[part] = [summed[part][0] + o["impulse"][0], summed[part][1] + o["impulse"][1],
@@ -153,7 +153,7 @@ class ExchangeWindowGPU(unittest.TestCase):
                                         summed[part][3] + o["mass_into_fluid"]]
                     for a, e in zip(windowed._download_raw().values, plain._download_raw().values):
                         same_bits(self, a, e)
-                out = windowed.read_exchange_window(0)
+                out = windowed.read_exchange_window(0, level=1)
                 self.assertEqual((out["time_end"] - out["time_start"], out["sample_count"]),
                                  (steps, 2 * steps))
                 self.assertEqual(substeps, 2 * steps)
@@ -165,11 +165,11 @@ class ExchangeWindowGPU(unittest.TestCase):
                     diff = np.abs(got - np.array(summed[part]))
                     self.assertTrue((diff <= bound).all(), f"{part}: {diff} > {bound}")
                     np.testing.assert_array_equal(o["force"], o["impulse"] / steps)
-                again = windowed.read_exchange_window(0)
+                again = windowed.read_exchange_window(0, level=1)
                 np.testing.assert_array_equal(again["total"]["impulse"], out["total"]["impulse"])
                 windowed.initialize(*arrays)
                 with self.assertRaises(RuntimeError):
-                    windowed.read_exchange_window(0)
+                    windowed.read_exchange_window(0, level=1)
 
 
 def audit_substeps(stepper, b, count, record):
@@ -212,7 +212,7 @@ class ExchangeWindowShiftedMRTGPU(unittest.TestCase):
                     run.step(2)
                     if run is not plain:
                         run.start_exchange_window()
-                b = windowed._level_boundaries[1]
+                b = windowed._boundaries[1]
                 count = len(b.tables["cell"])
                 record = dict(scale=np.zeros(4), substeps=0)
                 audit_substeps(windowed, b, count, record)
@@ -220,7 +220,7 @@ class ExchangeWindowShiftedMRTGPU(unittest.TestCase):
                 for _ in range(steps):
                     windowed.step(1)
                     plain.step(1)
-                    summed += components(windowed.read_raw_exchange(0))
+                    summed += components(windowed.read_raw_exchange(0, level=1))
                     for a, e in zip(windowed._download_raw().values, plain._download_raw().values):
                         same_bits(self, a, e)
                 split.step(steps)  # one call: same ledger arithmetic as five step(1) calls
@@ -228,31 +228,31 @@ class ExchangeWindowShiftedMRTGPU(unittest.TestCase):
                     same_bits(self, a, e)
                 rows = [np.empty((count, 4), np.float64) for _ in range(2)]
                 cuda.memcpy_dtoh(rows[0], b.window_integral)
-                cuda.memcpy_dtoh(rows[1], split._level_boundaries[1].window_integral)
+                cuda.memcpy_dtoh(rows[1], split._boundaries[1].window_integral)
                 np.testing.assert_array_equal(rows[0].view(np.uint64), rows[1].view(np.uint64))
 
-                out = windowed.read_exchange_window(0)
+                out = windowed.read_exchange_window(0, level=1)
                 self.assertEqual(record["substeps"], 2 * steps)
                 self.assertEqual((out["time_end"] - out["time_start"], out["sample_count"]),
                                  (steps, 2 * steps))
                 bound = 4 * (count * record["substeps"]) * EPS * record["scale"] + 1e-300
                 diff = np.abs(components(out) - summed)
                 self.assertTrue((diff <= bound).all(), f"{diff} > {bound}")
-                np.testing.assert_array_equal(split.read_exchange_window(0)["total"]["impulse"],
+                np.testing.assert_array_equal(split.read_exchange_window(0, level=1)["total"]["impulse"],
                                               out["total"]["impulse"])
 
                 windowed.step(0)  # no interval: window and its duration unchanged
-                again = windowed.read_exchange_window(0)
+                again = windowed.read_exchange_window(0, level=1)
                 self.assertEqual(again["time_end"], out["time_end"])
                 np.testing.assert_array_equal(again["total"]["impulse"], out["total"]["impulse"])
 
                 windowed.start_exchange_window()  # reset after a window
                 with self.assertRaisesRegex(RuntimeError, "at least one completed"):
-                    windowed.read_exchange_window(0)
+                    windowed.read_exchange_window(0, level=1)
                 record.update(scale=np.zeros(4), substeps=0)
                 windowed.step(1)
-                single = components(windowed.read_raw_exchange(0))
-                reset = windowed.read_exchange_window(0)
+                single = components(windowed.read_raw_exchange(0, level=1))
+                reset = windowed.read_exchange_window(0, level=1)
                 self.assertEqual(reset["time_end"] - reset["time_start"], 1)
                 bound = 4 * (count * record["substeps"]) * EPS * record["scale"] + 1e-300
                 self.assertTrue((np.abs(components(reset) - single) <= bound).all())

tests/validation/multigrid/test_fine_surface_gpu.py

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

tests/validation/multigrid/test_multiple_pinned_remap_gpu.py

@@ -56,7 +56,7 @@ def make(context):
 
 def topology(run, rectangles, objects):
     return RemapTopology.from_config(
-        CFG, run.patch, rectangles, child_objects=objects,
+        CFG, run.hierarchy.pairs[0], rectangles, child_objects=objects,
         child_jets=run._boundaries[2].jets,
     )
 
@@ -75,7 +75,7 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
         # Each jet selects a different body. Updating one must preserve the
         # other body's current per-link velocities, including its driven links.
         boundary = run._boundaries[2]
-        body_ids = run.child_boundary.curved_host[3]
+        body_ids = run.boundaries[2].curved_host[3]
         for index, amplitude in enumerate(amplitudes):
             before = boundary.tables['surface_velocity'].copy()
             run.set_jet_amplitude(index, amplitude, level=2)
@@ -87,8 +87,8 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
 
     def check_boundary(self, old, new, mesh, amplitudes):
         self.assertEqual(mesh.pinned_rectangles, PINS)
-        self.assertTrue(same_physical_boundary(old.child_boundary, new.child_boundary))
-        self.assertTrue(same_physical_boundary(new.child_boundary, mesh.boundary))
+        self.assertTrue(same_physical_boundary(old.boundaries[2], new.boundaries[2]))
+        self.assertTrue(same_physical_boundary(new.boundaries[2], mesh.boundary))
         # Compare physical cut-link rows, never layout-dependent blocked addresses.
         for run in (old, new):
             boundary = run._boundaries[2]
@@ -107,7 +107,7 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
             before, count = old.download(), old.step_count
             mesh = topology(old, START, objects)
             # Input order is irrelevant even with several pins and fluid regions.
-            new = owned.enter_context(old.remap_child(START[::-1], child_objects=objects))
+            new = owned.enter_context(old.remap_finest(START[::-1], objects=objects))
             self.assertIsNot(old, new)
             self.assertEqual((old.step_count, new.step_count), (count, count))
             self.check_boundary(old, new, mesh, amplitudes)
@@ -133,7 +133,7 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
                     previous = current
                     self.change_amplitudes(previous, amplitudes)
                     before, count = previous.download(), previous.step_count
-                    old_mesh = topology(previous, previous.child_patch.bounds, objects)
+                    old_mesh = topology(previous, previous.hierarchy.pairs[1].bounds, objects)
                     new_mesh = topology(previous, rectangles, objects)
                     expected = remap_synchronized(
                         old_mesh, new_mesh, before, base_count=count,
@@ -141,7 +141,7 @@ class MultiplePinnedRemapGPU(unittest.TestCase):
                         phase='postcollision_g',
                     ).g
                     current = owned.enter_context(
-                        previous.remap_child(rectangles, child_objects=objects))
+                        previous.remap_finest(rectangles, objects=objects))
                     self.assertIsNot(previous, current)
                     self.assertEqual((previous.step_count, current.step_count), (count, count))
                     self.check_boundary(previous, current, new_mesh, amplitudes)

tests/validation/multigrid/test_pinned_body_remap_gpu.py

@@ -36,8 +36,8 @@ class PinnedBodyRemapGPU(unittest.TestCase):
 
     def check_boundary(self, before, after):
         for name in ('flags','boundary_cells','curved_hits'):
-            exact(self,getattr(after.child_boundary,name),getattr(before.child_boundary,name),name)
-        for a,b in zip(after.child_boundary.curved_host,before.child_boundary.curved_host):
+            exact(self,getattr(after.boundaries[2],name),getattr(before.boundaries[2],name),name)
+        for a,b in zip(after.boundaries[2].curved_host,before.boundaries[2].curved_host):
             exact(self,a,b,'physical cut links')
         exact(self,after._boundaries[2].tables['surface_velocity'],
               before._boundaries[2].tables['surface_velocity'],'current per-link velocities')
@@ -49,7 +49,7 @@ class PinnedBodyRemapGPU(unittest.TestCase):
             old.initialize(*initial(old,False));old.step(2)
             old.set_jet_amplitude(0,.017,level=2)
             before=old.download()
-            new=owned.enter_context(old.remap_child(A,child_objects=objects))
+            new=owned.enter_context(old.remap_finest(A,objects=objects))
             self.check_boundary(old,new);self.check_reset(new,before)
             old.step();new.step()
             for a,b in zip(new.download(),old.download()): exact(self,a,b,'no-op next step bitwise')
@@ -64,13 +64,13 @@ class PinnedBodyRemapGPU(unittest.TestCase):
             for rects in (B,A):
                 previous=current;before=previous.download();count=previous.step_count
                 jets=previous._boundaries[2].jets
-                old=RemapTopology.from_config(CFG,previous.patch,previous.child_patch.bounds,
+                old=RemapTopology.from_config(CFG,previous.hierarchy.pairs[0],previous.hierarchy.pairs[1].bounds,
                                              child_objects=objects,child_jets=jets)
-                new=RemapTopology.from_config(CFG,previous.patch,rects,
+                new=RemapTopology.from_config(CFG,previous.hierarchy.pairs[0],rects,
                                              child_objects=objects,child_jets=jets)
                 expected=remap_synchronized(old,new,before,base_count=count,
                             level_steps=(count,2*count,4*count),phase='postcollision_g').g
-                current=owned.enter_context(previous.remap_child(rects,child_objects=objects))
+                current=owned.enter_context(previous.remap_finest(rects,objects=objects))
                 self.check_boundary(previous,current);self.check_reset(current,expected)
                 for a,b,om,nm in zip(current.download(),before,old.leaf_masks,new.leaf_masks):
                     exact(self,a[:,om & nm],b[:,om & nm],'physical-fluid overlap')

tests/validation/multigrid/test_raw_exchange_gpu.py

@@ -61,7 +61,7 @@ def census(run, index, buffer):
     layout, b = level.layout, level.boundary
     g = download(buffer, layout.blocked_shape)
     dense = layout.to_dense(g)
-    leaf = run._leaf[index]
+    leaf = run.leaf[index]
     ny, nx = leaf.shape
     solid = (b.host.flags.reshape(ny, nx) & SOLID) != 0
     y, x = np.nonzero(leaf)
@@ -92,8 +92,8 @@ def census(run, index, buffer):
 
 
 def fine(block):
-    """Two-level circle+jet construction on BOUNDS; block None keeps the dense fine level."""
-    return dict(refinements=[BOUNDS], level_blocks=[None, block],
+    """Two-level circle+jet construction on BOUNDS, blocked with a (12, 8) root tile."""
+    return dict(refinements=[BOUNDS], level_blocks=[(12, 8), block],
                 level_objects=[(), (circle(),)], level_jets=[(), (jet(),)])
 
 class RawExchangeGPU(unittest.TestCase):
@@ -105,7 +105,7 @@ class RawExchangeGPU(unittest.TestCase):
                     arrays = [state(l.shape) for l in measured._levels]
                     measured.initialize(*arrays)
                     plain.initialize(*arrays)
-                    b = measured._level_boundaries[1]
+                    b = measured._boundaries[1]
                     count = len(b.tables['cell'])
                     kernel = measured._levels[1].kernel
                     samples = []
@@ -119,7 +119,7 @@ class RawExchangeGPU(unittest.TestCase):
                         measured._stream.synchronize()
                         after = measured._levels[1].layout.to_dense(
                             download(dst, measured._levels[1].layout.blocked_shape))
-                        change = totals(after, measured._leaf[1])-before
+                        change = totals(after, measured.leaf[1])-before
                         expected = external + np.array([-raw[:, 0].astype(float).sum(),
                                                         -raw[:, 1].astype(float).sum(),
                                                          raw[:, 3].astype(float).sum()])
@@ -127,8 +127,8 @@ class RawExchangeGPU(unittest.TestCase):
                         samples.append(raw.copy())
                     measured._levels[1].kernel = audited_step
                     for speed in (0., .002, -.001, 0.):
-                        measured.set_jet_amplitude(0, speed)
-                        plain.set_jet_amplitude(0, speed)
+                        measured.set_jet_amplitude(0, speed, level=1)
+                        plain.set_jet_amplitude(0, speed, level=1)
                         samples.clear()
                         measured.step()
                         plain.step()
@@ -139,33 +139,13 @@ class RawExchangeGPU(unittest.TestCase):
                         np.testing.assert_array_equal(integral, expected)
                         last = download(b.last_raw, (count, 4))
                         np.testing.assert_array_equal(last, samples[-1])
-                        obs = measured.read_raw_exchange(0)
+                        obs = measured.read_raw_exchange(0, level=1)
                         assert obs['sample_count'] == 2
                         assert b.tables['aperture'].any()
                         measured.step(0)
-                        np.testing.assert_array_equal(measured.read_raw_exchange(0)['total']['impulse'],
+                        np.testing.assert_array_equal(measured.read_raw_exchange(0, level=1)['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)
@@ -178,7 +158,7 @@ class RawExchangeGPU(unittest.TestCase):
                  FixedPatchStepper(cfg, context=ctx, **settings) as plain:
                 arrays = [state(l.shape, ux=0.) for l in run._levels]
                 run.initialize(*arrays); plain.initialize(*arrays)
-                b = run._level_boundaries[3]
+                b = run._boundaries[3]
                 count = len(b.tables['cell'])
                 kernel = run._levels[3].kernel
                 samples, independent_samples = [], []
@@ -191,7 +171,7 @@ class RawExchangeGPU(unittest.TestCase):
                     run._stream.synchronize()
                     after = run._levels[3].layout.to_dense(
                         download(dst, run._levels[3].layout.blocked_shape))
-                    change = totals(after, run._leaf[3])-before
+                    change = totals(after, run.leaf[3])-before
                     expected = external + np.array([-independent[:, 0].sum(),
                                                     -independent[:, 1].sum(), independent[:, 3].sum()])
                     np.testing.assert_allclose(change, expected, rtol=2e-4, atol=8e-5)
@@ -233,7 +213,7 @@ class RawExchangeGPU(unittest.TestCase):
                  FixedPatchStepper(CFG, context=ctx, **fine((16, 20)), observe_raw_exchange=True) as amr:
                 amr.initialize(*(state(l.shape) for l in amr._levels))
                 uniform.initialize(state((9, 48, 64)))
-                b = amr._level_boundaries[1]
+                b = amr._boundaries[1]
                 count = len(b.tables['cell'])
                 links = obj.geometry.build_cut_links(64, 48)
                 np.testing.assert_array_equal([v.dir for v in links], b.tables['dir'])
@@ -241,7 +221,7 @@ class RawExchangeGPU(unittest.TestCase):
                 # Same fine-grid force arms after translation to the same center: the
                 # single level's static BFL gets the AMR Jet column and force arms; the raw
                 # output pointer is substituted at launch (no stepper-side accumulation).
-                root = uniform._level_boundaries[0]
+                root = uniform._boundaries[0]
                 velocity = uniform._upload(b.tables['surface_velocity'])
                 root.device.update(surface_velocity=velocity, rx=uniform._upload(b.tables['rx']),
                                    ry=uniform._upload(b.tables['ry']))
@@ -254,7 +234,7 @@ class RawExchangeGPU(unittest.TestCase):
                     samples.append(download(raw, (count, 4)))
                 root.kernel = measured_boundary
                 for speed in (0., .002, -.001, 0.):
-                    amr.set_jet_amplitude(0, speed)
+                    amr.set_jet_amplitude(0, speed, level=1)
                     cuda.memcpy_htod(velocity, b.tables['surface_velocity'])
                     samples.clear()
                     amr.step(); uniform.step(2)

tests/validation/multigrid/test_refinement_blocked_base.py

@@ -2,6 +2,12 @@
 
 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 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
@@ -18,7 +24,8 @@ from CelerisLab.lbm import refinement_stepper
 from CelerisLab.lbm.blocks import BlockLayout2D
 from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
 from CelerisLab.lbm.refinement import ACTIVE, OUTER_GHOST, FixedPatch2D
-from CelerisLab.lbm.refinement_stepper import FixedPatchStepper, blocked_base_tables
+from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+from CelerisLab.lbm.refinement_tables import blocked_coarse_tables
 from test_refinement_blocked_fine import NX, NY, BOUNDS, config, state, circle, actuation
 
 BASE_BLOCKS = ((24, 32), (12, 8), (3, 4), (1, 1))
@@ -102,7 +109,7 @@ class BlockedBaseCPU(unittest.TestCase):
                 leaves[[0, -1]] = False
             for by, bx in BASE_BLOCKS:
                 layout = BlockLayout2D.build(NY, NX, by, bx, periodic=(True, mode == "periodic"))
-                donors, targets, kind = blocked_base_tables(p, layout, order, mode != "periodic")
+                donors, targets, kind = blocked_coarse_tables(p, layout, order, mode != "periodic")
                 ly, lx = layout.local_shape
                 L = ly * lx
                 expected, source_cells = [], []
@@ -132,7 +139,7 @@ class BlockedBaseCPU(unittest.TestCase):
         p = FixedPatch2D.from_config(config("periodic"), BOUNDS)
         order = np.argsort(p.coalesce_coarse)
         layout = BlockLayout2D.build(NY, NX, 3, 4)
-        donors, targets, _ = blocked_base_tables(p, layout, order)
+        donors, targets, _ = blocked_coarse_tables(p, layout, order)
         dense = np.arange(9 * NY * NX, dtype=np.float32).reshape(9, NY, NX)
         initial = layout.to_blocked(dense)
         arrivals = np.arange(len(targets), dtype=np.float32) + 10000
@@ -172,7 +179,7 @@ class BlockedBaseCPU(unittest.TestCase):
         # With one parent block, all eight grazing sources are owner cores.
         # Each is read for four fine ghost children: late explosion corrupts32.
         single = BlockLayout2D.build(NY, NX, NY, NX)
-        donors, targets, _ = blocked_base_tables(p, single, order)
+        donors, targets, _ = blocked_coarse_tables(p, single, order)
         initial = single.to_blocked(dense)
         staged = initial.copy()
         staged.ravel()[targets.astype(np.intp)] = arrivals
@@ -198,7 +205,7 @@ class BlockedBaseCPU(unittest.TestCase):
                                         level_blocks=[(3, 4), FINE_BLOCK],
                                         level_objects=[(), (circle(actuation(.002)),)])
             try:
-                backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2))
+                backend.initialize(state(NY, NX, 1), state(*backend.hierarchy.pairs[0].fine_shape, 2))
                 coarse, fine = backend._levels
                 current, spare, fine_current = coarse.current, coarse.spare, fine.current
                 backend.step(2)
@@ -249,7 +256,7 @@ class BlockedBaseCPU(unittest.TestCase):
         with mocked_backend() as m:
             backend = FixedPatchStepper(config("periodic"), context=m["context"], refinements=[BOUNDS],
                                         level_blocks=[(3, 4), FINE_BLOCK])
-            backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2))
+            backend.initialize(state(NY, NX, 1), state(*backend.hierarchy.pairs[0].fine_shape, 2))
             m["functions"]["BlockCoalesce"].side_effect = RuntimeError("injected parent stage failure")
             current = backend._levels[0].current
             with self.assertRaisesRegex(RuntimeError, "injected parent stage"):
@@ -266,62 +273,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)
+    def test_level_blocks_required_and_validated(self):
+        context = Mock(activate=Mock(side_effect=AssertionError("CUDA reached")))
+        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, refinements=[BOUNDS], **kw)
+        with self.assertRaisesRegex(ValueError, "whole-level tile"):
+            FixedPatchStepper(cfg, context=context, refinements=[BOUNDS])
+        context.activate.assert_not_called()
+        # 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
@@ -22,7 +21,8 @@ from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
 from CelerisLab.lbm.descriptors import SOLID
 from CelerisLab.lbm.refinement import (ACTIVE, INNER_GHOST, OUTER_GHOST, FixedPatch2D,
                                        prepare_fine_boundary)
-from CelerisLab.lbm.refinement_stepper import FixedPatchStepper, blocked_fine_tables
+from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+from CelerisLab.lbm.refinement_tables import blocked_fine_tables
 
 NX, NY, BOUNDS, STEPS = 32, 24, (8, 6, 24, 18), 8
 FINE_BLOCKS = ((32, 40), (16, 20), (16, 10))  # fine window 32x40: 1x1, 2x2, 4x2
@@ -88,15 +88,13 @@ def seam_cells(layout, mapped):
 
 
 def raw(backend):
-    """Unmasked device arrays; blocked fine storage unpacked from owner cores."""
-    coarse, fine = backend._levels
-    out = [np.empty(coarse.shape, np.float32)]
-    out.append(np.empty(fine.shape if backend.fine_layout is None
-                        else backend.fine_layout.blocked_shape, np.float32))
-    cuda.memcpy_dtoh(out[0], coarse.current)
-    cuda.memcpy_dtoh(out[1], fine.current)
-    if backend.fine_layout is not None:
-        out[1] = backend.fine_layout.to_dense(out[1])
+    """Unmasked device arrays of every level; blocked storage unpacked from owner cores."""
+    out = []
+    for level in backend._levels:
+        array = np.empty(level.shape if level.layout is None else level.layout.blocked_shape,
+                         np.float32)
+        cuda.memcpy_dtoh(array, level.current)
+        out.append(array if level.layout is None else level.layout.to_dense(array))
     return out
 
 
@@ -127,8 +125,8 @@ 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]])
-            self.assertEqual(len(backend.fine_layout.halo_dst), 0)
+                                        level_blocks=[(12, 8), FINE_BLOCKS[0]])
+            self.assertEqual(len(backend.layouts[1].halo_dst), 0)
             self.assertIsNone(backend._levels[1].halo_dst)
             backend._gather_halo(backend._levels[1], Mock())  # no launch without pairs
             module.get_function.return_value.assert_not_called()
@@ -179,9 +177,9 @@ 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))
+            backend.initialize(state(NY, NX, 1), state(*backend.hierarchy.pairs[0].fine_shape, 2))
         return backend, calls
 
     def test_same_buffer_halo_bfl_step_launch_order(self):
@@ -189,8 +187,8 @@ class BlockedFineCPU(unittest.TestCase):
             backend, calls = self.mocked_backend([circle(spec)])
             c, f = backend._levels
             coarse, cur, spare = c.current, f.current, f.spare
-            flags, kind = backend._fine_flags, backend._levels[1].kind
-            velocity = backend._link_tables.get("surface_velocity")
+            flags, kind = backend._boundaries[1].flags, backend._levels[1].kind
+            velocity = backend._boundaries[1].device.get("surface_velocity")
             self.assertEqual(velocity is None, spec is None)
             backend.step(1)
             backend.close()
@@ -198,8 +196,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
@@ -214,8 +213,8 @@ class BlockedFineCPU(unittest.TestCase):
     def test_actuation_without_velocity_buffer_fails_before_bfl(self):
         for speed in (0.0, .002):  # any spec, even zero amplitude
             backend, calls = self.mocked_backend([circle(actuation(speed))])
-            self.assertIsNotNone(backend._fine_tables["surface_velocity"])
-            del backend._link_tables["surface_velocity"]
+            self.assertIsNotNone(backend._boundaries[1].tables["surface_velocity"])
+            del backend._boundaries[1].device["surface_velocity"]
             with self.assertRaisesRegex(RuntimeError, "Surface velocity device buffer is missing"):
                 backend.step(1)
             names = [name for name, _ in calls]
@@ -281,93 +280,5 @@ class BlockedFineCPU(unittest.TestCase):
             self.assertEqual(set(np.unique(kind).tolist()) - {ACTIVE, INNER_GHOST, OUTER_GHOST}, set())
 
 
-class BlockedFineGPU(unittest.TestCase):
-    def test_dense_and_blocked_fine_storage_are_bit_identical(self):
-        from CelerisLab.cuda.context import CudaContext
-        for mode in ("periodic", "slot_reservoir"):
-            cfg = config(mode)
-            patch = FixedPatch2D.from_config(cfg, BOUNDS)
-            g_coarse, g_fine = state(NY, NX, 1), state(*patch.fine_shape, 2)
-            results = []
-            with CudaContext(0) as context:
-                for fine_block in (None, *FINE_BLOCKS):
-                    with FixedPatchStepper(cfg, context=context, coarse_boundary=mode,
-                                           refinements=[BOUNDS],
-                                           level_blocks=[None, fine_block]) as backend:
-                        backend.initialize(g_coarse, g_fine)
-                        backend.step(STEPS)
-                        self.assertEqual(backend.step_count, STEPS)
-                        results.append((backend.download(), raw(backend)))
-            (dense_public, dense_raw), blocked = results[0], results[1:]
-            kind = patch.fine_kind
-            self.assertTrue(np.isnan(dense_raw[1][:, kind == OUTER_GHOST]).all(), mode)
-            # INNER populations may be NaN after hop 2 (pulled from poisoned OUTER);
-            # their exact pattern is pinned by the dense-vs-blocked comparison below.
-            self.assertTrue(np.isfinite(dense_raw[1][:, kind == ACTIVE]).all(), mode)
-            for fine_block, (public, raw_arrays) in zip(FINE_BLOCKS, blocked):
-                msg = f"{mode} fine_block={fine_block}"
-                for a, b in zip(public + tuple(raw_arrays), dense_public + tuple(dense_raw)):
-                    np.testing.assert_array_equal(a, b, err_msg=msg)  # values and NaN masks
-            print(f"{mode}: dense vs {len(blocked)} fine-block layouts bit-identical "
-                  f"after {STEPS} coarse steps", flush=True)
-
-    def test_passive_circle_dense_and_blocked_fine_storage_are_bit_identical(self):
-        from CelerisLab.cuda.context import CudaContext
-        for mode in ("periodic", "slot_reservoir"):
-            cfg = config(mode)
-            patch, boundary = fine_circle(mode)
-            solid = (boundary.flags.reshape(patch.fine_shape) & SOLID) != 0
-            g_coarse, g_fine = state(NY, NX, 1), state(*patch.fine_shape, 2)
-            results = []
-            with CudaContext(0) as context:
-                for fine_block in (None, *FINE_BLOCKS):
-                    with FixedPatchStepper(cfg, context=context, coarse_boundary=mode,
-                                           refinements=[BOUNDS], level_blocks=[None, fine_block],
-                                           level_objects=[(), (circle(),)]) as backend:
-                        backend.initialize(g_coarse, g_fine)
-                        backend.step(STEPS)
-                        results.append((backend.download(), raw(backend)))
-            (dense_public, dense_raw), blocked = results[0], results[1:]
-            fine = dense_raw[1]
-            self.assertTrue(np.isfinite(fine[:, (patch.fine_kind == ACTIVE) & ~solid]).all(), mode)
-            self.assertTrue(np.isnan(fine[:, solid]).all(), mode)
-            self.assertTrue(np.isnan(fine[:, patch.fine_kind == OUTER_GHOST]).all(), mode)
-            for fine_block, (public, raw_arrays) in zip(FINE_BLOCKS, blocked):
-                msg = f"circle {mode} fine_block={fine_block}"
-                for a, b in zip(public + tuple(raw_arrays), dense_public + tuple(dense_raw)):
-                    np.testing.assert_array_equal(a, b, err_msg=msg)  # values and NaN masks
-            print(f"circle {mode}: dense vs {len(blocked)} fine-block layouts bit-identical "
-                  f"after {STEPS} coarse steps", flush=True)
-
-    def test_actuated_circle_dense_and_blocked_fine_storage_are_bit_identical(self):
-        from CelerisLab.cuda.context import CudaContext
-        for mode in ("periodic", "slot_reservoir"):
-            for speed in (.002, -.002):
-                cfg = config(mode)
-                patch, boundary = fine_circle(mode, actuation(speed))
-                solid = (boundary.flags.reshape(patch.fine_shape) & SOLID) != 0
-                g_coarse, g_fine = state(NY, NX, 1), state(*patch.fine_shape, 2)
-                results = []
-                with CudaContext(0) as context:
-                    for fine_block in (None, *FINE_BLOCKS):
-                        with FixedPatchStepper(cfg, context=context, coarse_boundary=mode,
-                                               refinements=[BOUNDS], level_blocks=[None, fine_block],
-                                               level_objects=[(), (circle(actuation(speed)),)]) as backend:
-                            backend.initialize(g_coarse, g_fine)
-                            backend.step(STEPS)
-                            results.append((backend.download(), raw(backend)))
-                (dense_public, dense_raw), blocked = results[0], results[1:]
-                fine = dense_raw[1]
-                self.assertTrue(np.isfinite(fine[:, (patch.fine_kind == ACTIVE) & ~solid]).all())
-                self.assertTrue(np.isnan(fine[:, solid]).all())
-                self.assertTrue(np.isnan(fine[:, patch.fine_kind == OUTER_GHOST]).all())
-                for fine_block, (public, raw_arrays) in zip(FINE_BLOCKS, blocked):
-                    msg = f"actuated {speed} {mode} fine_block={fine_block}"
-                    for a, b in zip(public + tuple(raw_arrays), dense_public + tuple(dense_raw)):
-                        np.testing.assert_array_equal(a, b, err_msg=msg)  # values and NaN masks
-                print(f"actuated {speed:+.3f} {mode}: dense vs {len(blocked)} fine-block layouts "
-                      f"bit-identical after {STEPS} coarse steps", flush=True)
-
-
 if __name__ == "__main__":
     unittest.main()

tests/validation/multigrid/test_refinement_followup.py

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

tests/validation/multigrid/test_refinement_gpu.py

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

tests/validation/multigrid/test_refinement_jets.py

@@ -3,11 +3,13 @@
 Jets use the same coarse cell-edge frame as fine_objects. The equivalent jet
 of Astra's actuated fixture is SurfaceActuation(0.7, 0.85, U).as_jet((16, 12),
 2.15). Frozen GPU gates (8 coarse steps, both coarse boundary modes):
-  G1 jet run: dense vs 1x1/2x2/4x2 fine blocks bit-identical (values + NaN masks);
-  G2 jet vs legacy actuation (dense): max |dG| over finite fine ACTIVE fluid
+  G1 jet run: 1x1/2x2/4x2 fine blocks (blocked (12, 8) root) bit-identical to the
+     FINE_BLOCKS[0] reference (values + NaN masks);
+  G2 jet vs legacy actuation (blocked reference): max |dG| over finite fine ACTIVE fluid
      populations <= 1e-6 (velocity rows differ only by FP32 rounding; reported);
-  G3 amplitude change after 4 steps (+U -> -U), 4 more steps: dense vs blocked
-     bit-identical and the run differs from the unchanged run.
+  G3 amplitude change after 4 steps (+U -> -U), 4 more steps: blocked layouts
+     bit-identical to the FINE_BLOCKS[0] reference and the run differs from the
+     unchanged run.
 """
 from contextlib import nullcontext
 import unittest
@@ -26,6 +28,11 @@ from test_refinement_blocked_fine import (BOUNDS, FINE_BLOCKS, NX, NY, STEPS, ac
 SPEED = .002
 
 
+def root_and(fine_block):
+    """Blocked levels with a (12, 8) root tile."""
+    return [(12, 8), fine_block]
+
+
 def jet(speed=SPEED):
     return actuation(speed).as_jet((16.0, 12.0), 2.15)
 
@@ -42,7 +49,7 @@ def mocked(objects, jets, fine_block):
             patch.object(cuda, "mem_alloc", side_effect=lambda size: Mock()), \
             patch.object(cuda, "memcpy_htod", htod):
         backend = FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS],
-                                    level_blocks=[None, fine_block], level_objects=[(), tuple(objects)],
+                                    level_blocks=root_and(fine_block), level_objects=[(), tuple(objects)],
                                     level_jets=[(), tuple(jets)])
     return backend, htod
 
@@ -55,8 +62,8 @@ class JetStepperCPU(unittest.TestCase):
         C = (16.0, 12.0)
         whole = Jet(Sector(C, .1, math.pi, 3.), Radial(C, inward=False), amplitude=0.)
         backend, _ = mocked([circle()], [whole], FINE_BLOCKS[1])
-        _, jet_id = backend._jet_cache
-        self.assertGreater(len(np.unique(backend._fine_tables["block"][jet_id == 0])), 1)
+        _, jet_id = backend._boundaries[1].jet_cache
+        self.assertGreater(len(np.unique(backend._boundaries[1].tables["block"][jet_id == 0])), 1)
         backend.close()
 
     def test_fine_ownership_uses_remapped_physical_owner(self):
@@ -69,10 +76,10 @@ class JetStepperCPU(unittest.TestCase):
         source = SimObject(99, geometry=CircleGeometry(*C, 2.15), center=C, radius=2.15)
         whole = Jet(Sector(C, .1, math.pi, 3.), Radial(C, inward=False), amplitude=0.)
         backend, _ = mocked([source], [whole], FINE_BLOCKS[1])
-        _, jet_id = backend._jet_cache
-        owners = backend.fine_boundary.curved_host[3][jet_id == 0]
+        _, jet_id = backend._boundaries[1].jet_cache
+        owners = backend.boundaries[1].curved_host[3][jet_id == 0]
         self.assertEqual(np.unique(owners).tolist(), [0])
-        self.assertGreater(len(np.unique(backend._fine_tables["block"][jet_id == 0])), 1)
+        self.assertGreater(len(np.unique(backend._boundaries[1].tables["block"][jet_id == 0])), 1)
         self.assertEqual(source.obj_id, 99)
         backend.close()
         cfg = config("periodic")
@@ -90,8 +97,8 @@ class JetStepperCPU(unittest.TestCase):
     def test_blocked_jet_velocity_matches_legacy_rows(self):
         legacy, _ = mocked([circle(actuation(SPEED))], (), FINE_BLOCKS[1])
         driven, _ = mocked([circle()], [jet()], FINE_BLOCKS[1])
-        a = legacy._fine_tables["surface_velocity"]
-        b = driven._fine_tables["surface_velocity"]
+        a = legacy._boundaries[1].tables["surface_velocity"]
+        b = driven._boundaries[1].tables["surface_velocity"]
         diff = np.abs(a - b).max()
         print(f"fine jet vs legacy velocity rows: max |dU| = {diff:.3e} (exact={diff == 0})")
         self.assertLessEqual(diff, 4 * np.finfo(np.float32).eps * SPEED)
@@ -99,58 +106,53 @@ 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"]
+        before = backend._boundaries[1].tables["surface_velocity"].copy()
+        device = backend._boundaries[1].device["surface_velocity"]
         with patch.object(cuda, "memcpy_htod", htod):
             htod.reset_mock()
-            backend.set_jet_amplitude(0, -SPEED)
+            backend.set_jet_amplitude(0, -SPEED, level=1)
             (ptr, sent), _ = htod.call_args
             self.assertEqual(htod.call_count, 1)
             self.assertIs(ptr, device)
             np.testing.assert_allclose(sent, -before, rtol=1e-6, atol=0)
-            np.testing.assert_array_equal(backend._fine_tables["surface_velocity"], sent)
+            np.testing.assert_array_equal(backend._boundaries[1].tables["surface_velocity"], sent)
             with self.assertRaises(ValueError):
-                backend.set_jet_amplitude(0, 1.5)
+                backend.set_jet_amplitude(0, 1.5, level=1)
             self.assertEqual(htod.call_count, 1)
-        self.assertEqual(backend._jets[0].amplitude, -SPEED)
+        self.assertEqual(backend.jets[1][0].amplitude, -SPEED)
         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)):
                 with self.assertRaises(error):
-                    backend.set_jet_amplitude(0, bad)
+                    backend.set_jet_amplitude(0, bad, level=1)
                 self.assertEqual(backend._status, status)  # invalid input stays usable
             with self.assertRaises(IndexError):
-                backend.set_jet_amplitude(1, -SPEED)
+                backend.set_jet_amplitude(1, -SPEED, level=1)
             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._boundaries[1].tables["surface_velocity"].copy()
             boom = Mock(side_effect=RuntimeError("copy failed"))
             with patch.object(cuda, "memcpy_htod", boom), patch.object(cuda, "memcpy_htod_async", boom), \
                     patch.object(cuda, "pagelocked_empty", side_effect=lambda shape, dtype: np.empty(shape, dtype)):
                 with self.assertRaisesRegex(RuntimeError, "copy failed"):
-                    backend.set_jet_amplitude(0, -SPEED)
+                    backend.set_jet_amplitude(0, -SPEED, level=1)
             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)
+            self.assertEqual(backend.jets[1][0].amplitude, SPEED)
+            np.testing.assert_array_equal(backend._boundaries[1].tables["surface_velocity"], before)
             with self.assertRaisesRegex(RuntimeError, "failed"):
                 backend.step(1)
             with self.assertRaisesRegex(RuntimeError, "failed"):
-                backend.set_jet_amplitude(0, -SPEED)
+                backend.set_jet_amplitude(0, -SPEED, level=1)
             backend.close()
 
     def test_jet_overlapping_legacy_aperture_rejected_before_cuda(self):
         with patch.object(cuda, "Stream", side_effect=AssertionError("CUDA touched")):
             with self.assertRaisesRegex(ValueError, "SurfaceActuation"):
                 FixedPatchStepper(config("periodic"), context=Mock(), refinements=[BOUNDS],
+                                  level_blocks=root_and(FINE_BLOCKS[1]),
                                   level_objects=[(), (circle(actuation(SPEED)),)],
                                   level_jets=[(), (jet(),)])
 
@@ -161,14 +163,14 @@ def run(mode, fine_block, *, objects, jets, switch=None):
     patch_ = FixedPatch2D.from_config(cfg, BOUNDS)
     with CudaContext(0) as context, FixedPatchStepper(
             cfg, context=context, coarse_boundary=mode, refinements=[BOUNDS],
-            level_blocks=[None, fine_block], level_objects=[(), tuple(objects)],
+            level_blocks=root_and(fine_block), level_objects=[(), tuple(objects)],
             level_jets=[(), tuple(jets)]) as backend:
         backend.initialize(state(NY, NX, 1), state(*patch_.fine_shape, 2))
         if switch is None:
             backend.step(STEPS)
         else:
             backend.step(STEPS // 2)
-            backend.set_jet_amplitude(0, switch)
+            backend.set_jet_amplitude(0, switch, level=1)
             backend.step(STEPS - STEPS // 2)
         return backend.download(), raw(backend)
 
@@ -179,29 +181,29 @@ class JetStepperGPU(unittest.TestCase):
             patch_, boundary = fine_circle(mode)
             fluid = (patch_.fine_kind == ACTIVE) & (
                 (boundary.flags.reshape(patch_.fine_shape) & SOLID) == 0)
-            dense = run(mode, None, objects=[circle()], jets=[jet()])
-            self.assertTrue(np.isfinite(dense[1][1][:, fluid]).all(), mode)
-            for fine_block in FINE_BLOCKS:
+            reference = run(mode, FINE_BLOCKS[0], objects=[circle()], jets=[jet()])
+            self.assertTrue(np.isfinite(reference[1][1][:, fluid]).all(), mode)
+            for fine_block in FINE_BLOCKS[1:]:
                 public, raw_arrays = run(mode, fine_block, objects=[circle()], jets=[jet()])
-                for a, b in zip(public + tuple(raw_arrays), dense[0] + tuple(dense[1])):
+                for a, b in zip(public + tuple(raw_arrays), reference[0] + tuple(reference[1])):
                     np.testing.assert_array_equal(a, b, err_msg=f"G1 {mode} {fine_block}")
-            legacy = run(mode, None, objects=[circle(actuation(SPEED))], jets=())
-            diff = np.abs(dense[1][1][:, fluid] - legacy[1][1][:, fluid]).max()
+            legacy = run(mode, FINE_BLOCKS[0], objects=[circle(actuation(SPEED))], jets=())
+            diff = np.abs(reference[1][1][:, fluid] - legacy[1][1][:, fluid]).max()
             print(f"G2 {mode}: jet vs legacy fine max |dG| = {diff:.3e}", flush=True)
             self.assertLessEqual(diff, 1e-6, mode)
 
     def test_g3_amplitude_switch_layouts_identical(self):
         for mode in ("periodic", "slot_reservoir"):
-            dense = run(mode, None, objects=[circle()], jets=[jet()], switch=-SPEED)
-            steady = run(mode, None, objects=[circle()], jets=[jet()])
-            self.assertFalse(np.array_equal(dense[1][1], steady[1][1], equal_nan=True), mode)
-            for fine_block in FINE_BLOCKS:
+            reference = run(mode, FINE_BLOCKS[0], objects=[circle()], jets=[jet()], switch=-SPEED)
+            steady = run(mode, FINE_BLOCKS[0], objects=[circle()], jets=[jet()])
+            self.assertFalse(np.array_equal(reference[1][1], steady[1][1], equal_nan=True), mode)
+            for fine_block in FINE_BLOCKS[1:]:
                 public, raw_arrays = run(mode, fine_block, objects=[circle()], jets=[jet()],
                                          switch=-SPEED)
-                for a, b in zip(public + tuple(raw_arrays), dense[0] + tuple(dense[1])):
+                for a, b in zip(public + tuple(raw_arrays), reference[0] + tuple(reference[1])):
                     np.testing.assert_array_equal(a, b, err_msg=f"G3 {mode} {fine_block}")
-            print(f"G3 {mode}: amplitude switch dense vs {len(FINE_BLOCKS)} layouts identical",
-                  flush=True)
+            print(f"G3 {mode}: amplitude switch, {len(FINE_BLOCKS) - 1} layouts "
+                  f"identical to {FINE_BLOCKS[0]}", flush=True)
 
 
 if __name__ == "__main__":

tests/validation/multigrid/test_refinement_mrt_gpu.py

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

tests/validation/multigrid/test_refinement_multipatch.py

@@ -36,14 +36,14 @@ def level_fine(patch_set, windows):
     """Level-wide fine array with each patch window filled from `windows`."""
     out = np.full((9, *patch_set.fine_shape), np.nan, np.float32)
     for p, w in zip(patch_set.patches, windows):
-        y0, x0 = p.fine_level_origin
+        y0, x0 = p.fine_level_offset
         out[:, y0:y0 + p.fine_shape[0], x0:x0 + p.fine_shape[1]] = w
     return out
 
 
 def window(patch_set, index, fine):
     p = patch_set.patches[index]
-    y0, x0 = p.fine_level_origin
+    y0, x0 = p.fine_level_offset
     return fine[:, y0:y0 + p.fine_shape[0], x0:x0 + p.fine_shape[1]]
 
 
@@ -175,7 +175,7 @@ class MultipatchBodiesGPU(unittest.TestCase):
             for step in self.STEPS:
                 if switch is not None and step == 8:
                     for index, amplitude in switch:
-                        backend.set_jet_amplitude(index, amplitude)
+                        backend.set_jet_amplitude(index, amplitude, level=1)
                 backend.step(step - backend.step_count)
                 out.append(backend.download() + raw(backend))
         return out

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(),)))):
@@ -128,8 +128,10 @@ class ThreeLevelCPU(unittest.TestCase):
         run = mocked_child()
         b = run._levels[2].boundary
         before = b.tables["surface_velocity"].copy()
+        with self.assertRaises(TypeError):
+            run.set_jet_amplitude(0, -.002)  # level is required
         with self.assertRaises(IndexError):
-            run.set_jet_amplitude(0, -.002)  # default stays L1
+            run.set_jet_amplitude(0, -.002, level=1)  # level 1 has no Jets
         for level in (0, 3, True):
             with self.assertRaises(ValueError):
                 run.set_jet_amplitude(0, -.002, level=level)
@@ -138,7 +140,7 @@ class ThreeLevelCPU(unittest.TestCase):
             self.assertEqual(upload.call_count, 1)
             self.assertIs(upload.call_args.args[0], b.device["surface_velocity"])
         np.testing.assert_array_equal(b.tables["surface_velocity"], -before)
-        self.assertEqual(run._jets, ())
+        self.assertEqual(run.jets[1], ())
         with patch.object(cuda, "memcpy_htod", side_effect=RuntimeError("copy failed")):
             with self.assertRaisesRegex(RuntimeError, "copy failed"):
                 run.set_jet_amplitude(0, .002, level=2)
@@ -150,17 +152,17 @@ class ThreeLevelCPU(unittest.TestCase):
     def test_child_public_masks_include_solids(self):
         from CelerisLab.lbm.descriptors import SOLID
         run = mocked_child()
-        solid = run.child_boundary.flags.reshape(run.child_patch.fine_shape) & SOLID != 0
+        solid = run.boundaries[2].flags.reshape(run.hierarchy.pairs[1].fine_shape) & SOLID != 0
         self.assertTrue(solid.any())
-        self.assertFalse(run._child_leaf[solid].any())
+        self.assertFalse(run.leaf[2][solid].any())
         arrays = [np.ones(level.shape, np.float32) for level in run._levels]
         run._mask_nonleaves(*arrays)
         self.assertTrue(np.isnan(arrays[2][:, solid]).all())
         with patch.object(cuda, "memcpy_htod") as upload:
             run.initialize(*arrays)
         sent = upload.call_args_list[-1].args[1]
-        dense = run.child_layout.to_dense(sent)
-        self.assertTrue(np.isnan(dense[:, ~run._child_leaf]).all())
+        dense = run.layouts[2].to_dense(sent)
+        self.assertTrue(np.isnan(dense[:, ~run.leaf[2]]).all())
         run.close()
 
 
@@ -169,7 +171,7 @@ class ThreeLevelGPU(unittest.TestCase):
         from CelerisLab.cuda.context import CudaContext
         feq = np.asarray(compute_feq_d2q9(1.0, .03, -.02), np.float32)[:, None, None]
         with CudaContext(0) as ctx, FixedPatchStepper(CFG, context=ctx, **KW) as run:
-            shapes = ((16, 16), run.patch.fine_shape, run.child_patch.fine_shape)
+            shapes = ((16, 16), run.hierarchy.pairs[0].fine_shape, run.hierarchy.pairs[1].fine_shape)
             run.initialize(*(np.broadcast_to(feq, (9, *s)).copy() for s in shapes))
             run.step(8)
             self.assertEqual(run.step_count, 8)
@@ -182,7 +184,7 @@ class ThreeLevelGPU(unittest.TestCase):
     def test_nonequilibrium_masks_and_fresh_middle_feeds(self):
         from CelerisLab.cuda.context import CudaContext
         with CudaContext(0) as ctx, FixedPatchStepper(CFG, context=ctx, **KW) as run:
-            p1, p2 = run.patch, run.child_patch
+            p1, p2 = run.hierarchy.pairs[0], run.hierarchy.pairs[1]
             run.initialize(state(16, 16, 1), state(*p1.fine_shape, 2), state(*p2.fine_shape, 3))
             f = run._levels[1]
             feeds, original = [], run._advance
@@ -195,7 +197,7 @@ class ThreeLevelGPU(unittest.TestCase):
                 run.step(1)
             self.assertEqual(feeds, [cur, spare])  # current G, then the new G of hop 1
             run.step(7)
-            middle = run._middle_kind
+            middle = run.hierarchy.kinds[1]
             for g, leaf in zip(run.download(), (~p1.covered, middle == ACTIVE, p2.fine_kind == ACTIVE)):
                 self.assertTrue(np.isfinite(g[:, leaf]).all())
                 self.assertTrue(np.isnan(g[:, ~leaf]).all())
@@ -227,15 +229,15 @@ class ThreeLevelBodyGPU(unittest.TestCase):
                 if stop == 4 and switch:
                     run.set_jet_amplitude(0, -.002, level=2)
                 public, stored = run.download(), raw(run)
-                leaves = [~run.patch.covered, run._fine_leaf, run._child_leaf]
+                leaves = [~run.hierarchy.pairs[0].covered, run.leaf[1], run.leaf[2]]
                 if mode != "periodic":
                     leaves[0][[0, -1]] = False
                 for g, mask in zip(public, leaves):
                     self.assertTrue(np.isfinite(g[:, mask]).all())
                     self.assertTrue(np.isnan(g[:, ~mask]).all())
-                invalid = [run.patch.covered.copy(), run._middle_kind == OUTER_GHOST,
-                           (run.child_patch.fine_kind == OUTER_GHOST) |
-                           ((run.child_boundary.flags.reshape(run.child_patch.fine_shape) & SOLID) != 0)]
+                invalid = [run.hierarchy.pairs[0].covered.copy(), run.hierarchy.kinds[1] == OUTER_GHOST,
+                           (run.hierarchy.pairs[1].fine_kind == OUTER_GHOST) |
+                           ((run.boundaries[2].flags.reshape(run.hierarchy.pairs[1].fine_shape) & SOLID) != 0)]
                 if mode != "periodic":
                     invalid[0][[0, -1]] = True
                 for g, mask in zip(stored, invalid):

tests/validation/multigrid/test_refinement_three_level_multipatch_gpu.py

@@ -54,9 +54,9 @@ def jets(speed):
 
 
 def frames(run):
-    middle_origin = np.asarray(run.patch.fine_origin)
+    middle_origin = np.asarray(run.hierarchy.pairs[0].fine_origin)
     return ((np.zeros(2), 1.), (middle_origin, .5),
-            (middle_origin + .5*np.asarray(run.child_patch.fine_origin), .25))
+            (middle_origin + .5*np.asarray(run.hierarchy.pairs[1].fine_origin), .25))
 
 
 def physical_state(shape, origin, spacing):
@@ -92,15 +92,15 @@ class ThreeLevelMultipatchGPU(unittest.TestCase):
                     for index in range(len(drive)):
                         run.set_jet_amplitude(index, -.002, level=2)
                 public, stored = run.download(), raw(run)
-                masks = [~run.patch.covered, run._fine_leaf, run._child_leaf]
+                masks = [~run.hierarchy.pairs[0].covered, run.leaf[1], run.leaf[2]]
                 if mode != 'periodic':
                     masks[0][[0, -1]] = False
-                invalid = [run.patch.covered.copy(), run._middle_kind == OUTER_GHOST,
-                           run.child_patch.fine_kind == OUTER_GHOST]
+                invalid = [run.hierarchy.pairs[0].covered.copy(), run.hierarchy.kinds[1] == OUTER_GHOST,
+                           run.hierarchy.pairs[1].fine_kind == OUTER_GHOST]
                 if mode != 'periodic':
                     invalid[0][[0, -1]] = True
-                if run.child_boundary is not None:
-                    invalid[2] |= (run.child_boundary.flags.reshape(run.child_patch.fine_shape) & SOLID) != 0
+                if run.boundaries[2] is not None:
+                    invalid[2] |= (run.boundaries[2].flags.reshape(run.hierarchy.pairs[1].fine_shape) & SOLID) != 0
                 for level, (g, mask, dense, bad) in enumerate(zip(public, masks, stored, invalid)):
                     label = f'{mode} stop={stop} level={level}'
                     self.assertTrue(np.isfinite(g[:, mask]).all(), label)

tests/validation/multigrid/test_runtime_child_remap_gpu.py

@@ -39,7 +39,7 @@ def initial(run, uniform):
 
 class RuntimeChildRemapGPU(unittest.TestCase):
     def check_public(self, run):
-        for array, leaf in zip(run.download(), (~run.patch.covered, run._fine_leaf, run._child_leaf)):
+        for array, leaf in zip(run.download(), (~run.hierarchy.pairs[0].covered, run.leaf[1], run.leaf[2])):
             self.assertTrue(np.isfinite(array[:, leaf]).all())
             self.assertTrue(np.isnan(array[:, ~leaf]).all())
 
@@ -57,7 +57,7 @@ class RuntimeChildRemapGPU(unittest.TestCase):
             old.initialize(*initial(old, False))
             old.step(2)
             before = old.download()
-            new = owned.enter_context(old.remap_child(A))
+            new = owned.enter_context(old.remap_finest(A))
             self.assertEqual((old.step_count, new.step_count), (2, 2))
             self.check_reset(new, before)
             for got, want in zip(old.download(), before): exact(self, got, want, 'old unchanged')
@@ -75,12 +75,12 @@ class RuntimeChildRemapGPU(unittest.TestCase):
                 current.step(2)
                 for rectangles in (B, A):
                     before, count = current.download(), current.step_count
-                    old_mesh = RemapTopology.from_config(CFG, current.patch, current.child_patch.bounds)
-                    new_mesh = RemapTopology.from_config(CFG, current.patch, rectangles)
+                    old_mesh = RemapTopology.from_config(CFG, current.hierarchy.pairs[0], current.hierarchy.pairs[1].bounds)
+                    new_mesh = RemapTopology.from_config(CFG, current.hierarchy.pairs[0], rectangles)
                     expected = remap_synchronized(old_mesh, new_mesh, before, base_count=count,
                         level_steps=(count, 2*count, 4*count), phase='postcollision_g').g
                     previous = current
-                    current = owned.enter_context(previous.remap_child(rectangles))
+                    current = owned.enter_context(previous.remap_finest(rectangles))
                     self.assertEqual(current.step_count, count)
                     self.check_reset(current, expected)
                     for got, want, old_g, old_leaf, new_leaf in zip(

tests/validation/multigrid/test_shear_refinement_example_gpu.py

@@ -8,14 +8,17 @@ from contextlib import ExitStack
 import unittest
 import numpy as np
 from CelerisLab.cuda.context import CudaContext
-from research.shear_refinement.example import (CHILDREN, make, initialize, observe,
-                                                topology, objects)
+from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+from research.shear_refinement.example import (CFG, PARENTS, CHILDREN, initialize,
+                                                observe, topology, objects)
 
 
 class ShearRefinementGPU(unittest.TestCase):
     def test_measured_growth_matches_explicit_manual_continuation(self):
         with CudaContext(0) as ctx, ExitStack() as owned:
-            old = owned.enter_context(make(ctx))
+            old = owned.enter_context(FixedPatchStepper(
+                CFG, context=ctx, refinements=[PARENTS, CHILDREN],
+                level_blocks=[(4, 4)] * 3, level_objects=[(), (), objects()]))
             initialize(old)
             old.step(2)
             before = old.download()
@@ -25,10 +28,10 @@ class ShearRefinementGPU(unittest.TestCase):
             self.assertGreater(report['requested'], 0)
             self.assertEqual(report['corridor_unknown'], 0)
             self.assertEqual(report['corridor_unavailable'], 0)
-            auto = owned.enter_context(old.remap_child(chosen, child_objects=objects()))
-            explicit = owned.enter_context(old.remap_child(manual, child_objects=objects()))
+            auto = owned.enter_context(old.remap_finest(chosen, objects=objects()))
+            explicit = owned.enter_context(old.remap_finest(manual, objects=objects()))
             for run in (auto, explicit):
-                self.assertEqual(run.child_patch.bounds[0], CHILDREN[0])
+                self.assertEqual(run.hierarchy.pairs[1].bounds[0], CHILDREN[0])
                 self.assertEqual(run.step_count, 2)
                 run.step(2)
             for a, b, leaf in zip(auto.download(), explicit.download(), topology(manual).leaf_masks):

tests/validation/multigrid/test_shifted_jets_gpu.py

@@ -58,7 +58,7 @@ def census(run,index,buffer):
     """Independent physical link/face census of actual staged receiver windows."""
     level=run._levels[index];layout=level.layout;b=level.boundary
     stored=download(buffer,layout.blocked_shape);dense=layout.to_dense(stored)
-    leaf=run._leaf[index];ny,nx=leaf.shape
+    leaf=run.leaf[index];ny,nx=leaf.shape
     solid=(b.host.flags.reshape(ny,nx)&SOLID)!=0
     y,x=np.nonzero(leaf);owner=layout.owner[y,x].astype(np.intp)
     cells=int(np.prod(layout.local_shape));block,cell=np.divmod(owner,cells)
@@ -93,7 +93,7 @@ class Ledger:
             self.max_row=max(self.max_row,float(abs(raw.astype(float)-independent).max()))
             self.original(src,dst,*args,**kwargs);run._stream.synchronize()
             after=level.layout.to_dense(download(dst,level.layout.blocked_shape))
-            change=raw_totals(after,run._leaf[index])-before
+            change=raw_totals(after,run.leaf[index])-before
             expected=external+np.array([-independent[:,0].sum(),-independent[:,1].sum(),independent[:,3].sum()])
             np.testing.assert_allclose(change,expected,rtol=2e-4,atol=8e-5)
             self.max_balance=max(self.max_balance,float(abs(change-expected).max()))
@@ -102,7 +102,7 @@ class Ledger:
 
     def step(self):
         run=self.run;self.samples.clear();run.step()
-        b=run._level_boundaries[self.index];h=run._exchange_spacing[self.index]
+        b=run._boundaries[self.index];h=b.spacing
         self.test.assertEqual(len(self.samples),2**self.index)
         expected=np.zeros_like(self.samples[0],dtype=np.float64)
         weights=np.array([h*h,h*h,h*h*h,h*h])
@@ -113,7 +113,7 @@ class Ledger:
         for name in ('impulse','angular_impulse','mass_into_fluid'):
             np.testing.assert_allclose(obs['total'][name],np.asarray(obs['wall'][name])+obs['aperture'][name],rtol=0,atol=2e-13)
         self.test.assertEqual(obs['sample_count'],2**self.index)
-        for a,leaf in zip(run._download_raw().values,run._leaf):
+        for a,leaf in zip(run._download_raw().values,run.leaf):
             self.test.assertTrue(np.isfinite(a[:,leaf]).all())
             self.test.assertTrue((1+a[:,leaf].sum(axis=0,dtype=np.float64)>0).all())
         return obs
@@ -140,10 +140,10 @@ class ShiftedJetsGPU(unittest.TestCase):
                         level_objects=[(),(obj,)],level_jets=[(),(whole,)],observe_raw_exchange=True) as run:
                     run._initialize_raw(initial(run,rho,u=.013,v=0.))
                     expected=run._download_raw();ledger=Ledger(self,run,1)
-                    self.assertTrue(run._level_boundaries[1].tables['aperture'].all())
+                    self.assertTrue(run._boundaries[1].tables['aperture'].all())
                     for _ in range(4):ledger.step()
                     error=max(float(abs(a[:,leaf]-b[:,leaf]).max()) for a,b,leaf in
-                              zip(run._download_raw().values,expected.values,run._leaf))
+                              zip(run._download_raw().values,expected.values,run.leaf))
                     self.assertLessEqual(error,1e-6)  # inherited uniform-state cap
                     print(f'matched whole-circle rho={rho}: max_abs_raw_state={error:.9g}',flush=True)
                     ledger.report(f'matched whole-circle rho={rho}')
@@ -161,7 +161,7 @@ class ShiftedJetsGPU(unittest.TestCase):
                 ledger=Ledger(self,measured,1)
                 for speed in SPEEDS:
                     before=measured._download_raw();count=measured.step_count
-                    measured.set_jet_amplitude(0,speed);plain.set_jet_amplitude(0,speed)
+                    measured.set_jet_amplitude(0,speed, level=1);plain.set_jet_amplitude(0,speed, level=1)
                     self.exact_raw(before,measured._download_raw());self.assertEqual(measured.step_count,count)
                     ledger.step();plain.step()
                     self.exact_raw(measured._download_raw(),plain._download_raw())
@@ -193,15 +193,15 @@ class ShiftedJetsGPU(unittest.TestCase):
             old._initialize_raw(initial(old,.99));ledger=Ledger(self,old,2)
             for speed in (.002,-.001):old.set_jet_amplitude(0,speed,level=2);ledger.step()
             before=old._download_raw();count=old.step_count
-            topology=lambda rects:RemapTopology.from_config(cfg,old.patch,rects,child_objects=(obj,),
-                child_jets=old._level_jets[2],_shifted_blocked=True)
+            topology=lambda rects:RemapTopology.from_config(cfg,old.hierarchy.pairs[0],rects,child_objects=(obj,),
+                child_jets=old.jets[2],_shifted_blocked=True)
             expected=_remap_stored(topology(a),topology(b),before,shifted=True,base_count=count,
                 level_steps=(count,2*count,4*count),phase='postcollision_g')
-            new=owned.enter_context(old.remap_child(b,child_objects=(obj,)))
+            new=owned.enter_context(old.remap_finest(b,objects=(obj,)))
             self.exact_raw(before,old._download_raw())
             self.exact_raw(_StoredPopulations(expected.g,True),new._download_raw())
-            self.assertEqual(new._level_jets[2][0].amplitude,-.001)
-            np.testing.assert_array_equal(new._level_boundaries[2].tables['surface_velocity'],old._level_boundaries[2].tables['surface_velocity'])
+            self.assertEqual(new.jets[2][0].amplitude,-.001)
+            np.testing.assert_array_equal(new._boundaries[2].tables['surface_velocity'],old._boundaries[2].tables['surface_velocity'])
             newledger=Ledger(self,new,2);ledger.step();newledger.step()
             self.assertEqual((old.step_count,new.step_count),(count+1,count+1))
             ledger.report('L2 source');newledger.report('L2 remapped')

tests/validation/multigrid/test_shifted_remap_gpu.py

@@ -61,7 +61,7 @@ class ShiftedRemapGPU(unittest.TestCase):
     def check_finite(self,run):
         raw=run._download_raw()
         self.assertTrue(raw.shifted)
-        for a,leaf in zip(raw.values,run._leaf):
+        for a,leaf in zip(raw.values,run.leaf):
             self.assertTrue(np.isfinite(a[:,leaf]).all())
             self.assertTrue(np.isnan(a[:,~leaf]).all())
             self.assertTrue((1+a[:,leaf].sum(axis=0,dtype=np.float64)>0).all())
@@ -71,7 +71,7 @@ class ShiftedRemapGPU(unittest.TestCase):
             old=owned.enter_context(make(ctx))
             old._initialize_raw(initial(old));old.step(2)
             before=old._download_raw()
-            new=owned.enter_context(old.remap_child(A))
+            new=owned.enter_context(old.remap_finest(A))
             self.check_reset(new,before.values)
             self.snapshot_equal(old._download_raw(),before,'old unchanged')
             self.snapshot_equal(new._download_raw(),before,'identity raw copy')
@@ -85,11 +85,11 @@ class ShiftedRemapGPU(unittest.TestCase):
             old=owned.enter_context(make(ctx))
             old._initialize_raw(initial(old));old.step(2)
             before=old._download_raw();count=old.step_count
-            old_mesh=RemapTopology.from_config(CFG,old.patch,A)
-            new_mesh=RemapTopology.from_config(CFG,old.patch,B)
+            old_mesh=RemapTopology.from_config(CFG,old.hierarchy.pairs[0],A)
+            new_mesh=RemapTopology.from_config(CFG,old.hierarchy.pairs[0],B)
             expected=_remap_stored(old_mesh,new_mesh,before,shifted=True,
                 base_count=count,level_steps=(count,2*count,4*count),phase='postcollision_g')
-            new=owned.enter_context(old.remap_child(B))
+            new=owned.enter_context(old.remap_finest(B))
             self.check_reset(new,expected.g)
             self.snapshot_equal(new._download_raw(),_StoredPopulations(expected.g,True),'host raw fill')
             for a,b,om,nm in zip(before.values,new._download_raw().values,old_mesh.leaf_masks,new_mesh.leaf_masks):
@@ -103,7 +103,7 @@ class ShiftedRemapGPU(unittest.TestCase):
         with CudaContext(0) as ctx,make(ctx) as run:
             run._initialize_raw(initial(run));run.step(2)
             before=run._download_raw()
-            with self.assertRaises(ValueError):run.remap_child([])
+            with self.assertRaises(ValueError):run.remap_finest([])
             self.snapshot_equal(run._download_raw(),before,'rejection preserves source')
             self.assertEqual(run.step_count,2)
             run.step();self.check_finite(run)

tests/validation/run_refined_dfg.py

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

tests/validation/time_window_dfg.py

@@ -73,7 +73,7 @@ def run_variant(backend, variant, n_intervals, deadline):
             else:
                 for _ in range(INTERVAL):
                     backend.step(1)
-            totals.append(components(backend.read_exchange_window(0)))
+            totals.append(components(backend.read_exchange_window(0, level=1)))
         t.check_deadline(deadline)
     return np.array(totals), (np.array(scales) if scales else None)