docs: Dütsch protocol — adopt Grok's cycle-index labels (K = first scored cycle)

Same cycles and pairs as 8505a723, relabelled: proof pairs 3→4..5→6 give
K_case = 7, window K..K+2 with remain-settled pairs, last scored cycle 11.

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

@@ -40,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,116 @@
+# 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). All three received parent/candidate
+generated-config and host/GPU equality checks. Only the collision slice ran a new PTX
+matrix (a4fc46d, whose src is identical to c5a3cdb); 65339ab and 442979b reused prior
+PTX evidence after unchanged CUDA/compiler and generated-input checks. In the collision
+matrix, only the two channel-root SRT modules changed PTX (omega load/clamp moved);
+their parent/candidate outputs matched exactly in the exercised cases.
+
+## 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

@@ -6,8 +6,8 @@ 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 blocked slot/reservoir root reads its slot from the `d_block_wall` constant (the former dense route used `WALL_*` macros; both were bit-identical in the parity gates). | Every owner loads its own module (`compiler.load_module` = `module_from_file`, no cache). A loaded module is never shared between live owners: globals are per module, so sharing would share mutable `d_params`. |
-| 2. Geometry, layout, tables | flags, curved-link tables, block layouts and halo pairs, transfer (explode/coalesce) tables, block kinds | Immutable for the life of an instance. Uniform grid: rebuilt by body sync. `FixedPatchStepper`: a geometry change is a new instance (`remap_child`), committed only after a successful transfer. | No in-place table edits. |
-| 3. Controls and dependent boundary columns | Jet amplitude -> per-link `surface_velocity`; body actions (uniform grid) | Between steps, without a table or module rebuild (`set_jet_amplitude`). | Host copy commits only after the device upload succeeded; a device error fails the instance. |
+| 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:

docs/moving_bodies.md

@@ -0,0 +1,403 @@
+# Moving circles on the uniform dense grid — S0 DRAFT contract
+
+Status: **DRAFT, not frozen** for the force RHS, kernels and physical
+acceptance. Host-only S1 (state, validation, guards; no kernels) is on branch
+`claude/moving-bodies-s1`; positive-step runs with a moving body are rejected
+until S2. When S2 lifts that guard, `run(n > 1)` with a pending target must
+raise `ValueError` before mutation (4a), and `run(0)` must not upload a motion
+record or consume the pending target. Base: main
+`b2034f87e03656c5f67ab88562a70d6c338b9918`. Scope approved by Frank (via TRON,
+thread C0C78UFKGRM/1793335557.460509): uniform-grid prescribed translation and
+single-cylinder free transverse VIV. Out of scope: refinement/blocked path,
+pinball assembly, EsoPull, shifted storage, D3Q19, DynamisLab, manuscript.
+
+Items marked **OPEN** are undecided; items marked **CANDIDATE** are proposals.
+
+## 1. Working standards
+
+- Lean, readable kernels in the existing idiom (`kernels/step/*.cu`, shared
+  `apply_bouzidi_link`); Python manages state and validation. No framework,
+  no compatibility layer, no new dependency.
+- The existing stationary/rotating path (`CurvedBoundaryKernel`, host link
+  packing, `sync_bodies`) stays byte-identical. Moving circles are a separate,
+  explicitly enabled path with its own numerical contract (this file).
+- One authority per state (section 4). No per-step full-domain download,
+  rebuild or recompile; local swept-band device updates are in scope.
+- Small slices (S1 host spec, S2 kernels without motion, S3 motion + refill,
+  S4 oscillator); each reviewed SHA is frozen before targeted runs.
+- Tolerances, benchmark parameters and the force composition used in the
+  structural RHS are declared before the runs that judge them. Fresh vs carried
+  evidence and existing baseline failures are recorded; smoke, convergence and
+  physical validation are reported separately. No undeclared force filtering.
+- External code (waLBerla GPLv3, FluidX3D source-available, VIVSIM) is a
+  concept reference only; nothing is copied.
+
+## 2. Current update order (main b2034f87, dense double buffer)
+
+`Simulation.run(n)`: upload host action (if `upload_act`) → zero obs (if
+`zero_obs`) → `DoubleBufferStepper.step(n)` → download obs → synchronize.
+
+Per step (`lbm/stepper.py` `step`):
+
+1. `CurvedBoundaryKernel` over the packed cut links: reads post-collision
+   `G^{n-1}` at `k_f` (and `k_f - c` for q < 1/2), adds the moving-wall term
+   (`Uw = -omega*ry, Vw = omega*rx` + optional per-link `surface_velocity`),
+   writes the reflected value into slot `opp` of the solid node `k_s`, and
+   accumulates `c(f_toward + f_reflected)` and torque into `obs`.
+2. `OneStep` pull/collide `ddf → temp` (collides FLUID cells; obstacle
+   interiors flagged `SOLID|OBSTACLE|BC_CURVED`), then buffer swap.
+3. `ForceRegionKernel`, `SensorKernel`.
+
+Raw `obs` force is the force **on the body** (link derivation in section 6;
+DFG docs report `Cd > 0` from it). Translation `vx, vy` in the action slots is
+currently stored but not read by any curved kernel.
+
+## 3. Proposed update order (moving circles)
+
+Per step n, before the existing step 1:
+
+- **M1 `MovingCircleCells`** (one thread per swept-band cell): classify with the
+  analytic predicates `in_old = inside(x, c^{n-1})`, `in_new = inside(x, c^n)`;
+  write `flag = in_new ? SOLID|OBSTACLE|BC_CURVED : FLUID`; covered cells
+  (`!in_old && in_new`) accumulate `node_cover`; fresh cells
+  (`in_old && !in_new`) are refilled (section 5) and accumulate `node_uncover`.
+- **M2 `MovingCircleLinks`** (one thread per band cell × direction): compute
+  q and hit point at `c^n` (FP64, same formula as `find_circle_intersection`,
+  same clamp [0.01, 1] and fallback classes as `CircleGeometry.build_cut_links`),
+  skip non-cut slots, call `apply_bouzidi_link` with
+  `u_w = (vx^n, vy^n) + omega × (hit - c^n)`, accumulate the link terms
+  (section 6).
+
+Then the unchanged steps 1–3. M2 must follow M1 because a q < 1/2 link may
+read `k_f - c` at a fresh cell.
+
+Exclusivity: moving bodies are never packed into the `CurvedBoundaryKernel`
+link arrays (`n_curved` covers stationary/rotating bodies only), so M2 is the
+sole writer of their solid slots and the sole source of their link force; no
+double write and no double `obs`. The stationary path's packing, launch and
+arithmetic are unchanged.
+
+Completeness: after M1, every (fluid cell, direction) pair in the swept band
+whose neighbour is `in_new` is visited by M2. This holds because the band is
+the circle AABB at the declared extreme offsets plus 2 cells, and the host
+rejects per-step moves with `|c_new - c_old|_2 > 0.5` and bands touching the domain edge,
+other bodies, sensors or force regions. Each M2 write lands at
+`index_f(k_s, opp)`, the exact slot the adjacent fluid cell pulls, so every
+fluid pull from a solid slot of a moving body is BFL-written.
+
+## 4. State, masks and authority
+
+- Host `ObjectState` is the single authority for one displacement
+  `offset = (dx, dy)` relative to the immutable base centre, plus
+  `velocity = (vx, vy)` and `omega`. Centres are always derived
+  (`c = base + offset`), never stored as independent state. Velocity is
+  supplied, never differenced from offsets.
+- `c^{n-1}` is derived from the last *committed* offset, `c^n` from the
+  offset being applied. Per body a small device **motion record**
+  `{c^{n-1}, c^n, u^n}` is uploaded with the action buffer; action-slot meaning
+  is not overloaded. (CANDIDATE layout; OPEN.)
+- Fail-closed commit: an enqueued upload is not a completed one. The new
+  offset becomes the committed offset only after the stream synchronization
+  that follows the upload and the step group succeeds. Any upload or later
+  synchronization failure leaves host and device motion state unknown: the
+  simulation is marked unusable (existing `_require_usable` lifecycle), and
+  no API may report a host/device match afterwards.
+
+### 4a. Applied geometry vs pending target (CANDIDATE, narrow API)
+
+Two distinct facts per moving body, one authority each:
+`applied` = offset of the geometry last applied on the device, changed only
+by a successful commit; `pending` = an optional requested target
+`(offset, velocity)`, changed only by `set_body`. A setter returning proves
+nothing about the device. The uploaded motion record is derived each run:
+`{c_old = base + applied, c_new = base + (pending.offset if pending is not
+None else applied), u = pending.velocity if pending is not None else 0}`
+(explicit `is not None`, so a zero offset or zero velocity is a valid target). Translational wall velocity exists only on a step
+that advances to a pending target; a hold step (`c_old == c_new`) has
+translational `u_w = 0` (spin unchanged), so a non-advancing circle never
+acts as a sliding wall. One moving step consumes one target.
+
+| situation | contract |
+|---|---|
+| `set_body` target | validated against `applied` (box, `|Δc|_2 ≤ 0.5` cell, speed bound) before mutation; on success replaces `pending` |
+| target with `offset == applied` and nonzero translational velocity | `ValueError` before mutation (would be a hold step with sliding wall velocity) |
+| two setters before a run | last wins; each validated against `applied`, not against the earlier pending target |
+| `run(1)` with pending | upload record, step, sync; on success `applied := pending.offset` (applied stores only the offset), `pending := None` |
+| `run(1)` without pending | hold: `c_old == c_new`, no cover/uncover, links at `applied`, translational `u_w = 0` |
+| `run(n > 1)` with pending | `ValueError` before mutation (single-step-only moving API) |
+| `run(n > 1)`, no pending | allowed: hold for n steps, translational `u_w = 0`, spin only |
+| `run(0)` | no motion upload, no commit; `pending` kept |
+| validation failure | `ValueError`, nothing mutated, retry allowed |
+| upload/launch/sync failure after enqueue | simulation unusable; `applied` and `pending` are not reported as device truth; no re-entry |
+
+Displacement norm: Euclidean per-step `|c_new - c_old|_2 ≤ 0.5` lattice
+cells (declared bound; CANDIDATE value).
+- Membership is never read from a flag that M1 is overwriting. The host
+  verifies at `add_body` that every swept-band cell has base flag FLUID, so the
+  band's base state is known without reading `flag_gpu`.
+- `flag_gpu` is the device authority for band cells while a moving body
+  exists; the host flag mirror is stale for those cells and is not used.
+
+## 5. Refill of fresh cells (CANDIDATE)
+
+- Stored time level: everything in `ddf` at M1 is post-collision `G^{n-1}`; the
+  refill writes all 9 slots at that same level, so every slot of a fresh cell
+  that the following pull or a q < 1/2 link may read is defined. Donor values
+  are read at the same level `G^{n-1}` (old-fluid cells, untouched by M1).
+- `G = f_eq(rho_bar, u_w(x)) + G_neq(x_d)` where `rho_bar` = mean density of
+  valid neighbours and `x_d` = the valid neighbour best aligned with the outward
+  normal `(x - c^n)/|x - c^n|`; `G_neq(x_d) = G(x_d) - f_eq(rho(x_d), u(x_d))`.
+- Valid donor: inside the domain, base FLUID, `!in_old && !in_new`.
+- Sets: `V(x)` = the 8 lattice neighbours that are in the domain, base FLUID
+  and `!in_old && !in_new`. `rho_bar` = mean density over `V(x)`. NEQ donor =
+  the member of `V(x)` with the largest `c·n` among those with `c·n > 0`
+  (`n` = outward normal at `c^n`); none → NEQ part 0 (counted as `eq_only`).
+- `V(x)` empty → no defined density: the cell is refilled with
+  `f_eq(RHO, u_w(x))` so all 9 slots stay defined, the `no_donor` counter is
+  incremented, and the host treats `no_donor > 0` after the run as a hard
+  failure (simulation unusable). It is not a silent fallback.
+- Non-emptiness is NOT assumed from the displacement bound. Geometry
+  assumptions: circle radius `R ≥ R_min` (CANDIDATE 4 cells), `|Δc|_2 ≤ 0.5`.
+  S1 includes a host enumeration over sub-cell centre positions and
+  displacement directions on a fine grid within these bounds that checks every
+  fresh cell has nonempty `V(x)`; it is evidence for those bounds, not a
+  proof, and runtime keeps the counter.
+- Fresh cells only write, old-fluid cells only read: no race.
+
+## 6. Discrete momentum ledger (to review; not a selection)
+
+Notation per step n: `F^n` = fluid set after M1; links L = (k_f, dir) at
+`c^n`; `f_t = G_dir(k_f)`; `f_r` = value M2 writes into `(k_s, opp)`.
+Fluid momentum before streaming:
+
+    P_pre = P(G^{n-1} on F^{n-1}) - sum_cover sum_q c_q G_q + sum_fresh sum_q c_q G_q^refill
+
+During the pull, every population entering a fluid cell comes from a fluid
+cell except the BFL-written slots; `f_t` leaves into solid and is never pulled
+by a fluid cell. Collision conserves momentum for SRT/MRT without forcing.
+Hence, in exact arithmetic,
+
+    P_fluid^n - P_fluid^{n-1} = - sum_L c (f_t + f_r) - I_cover + I_uncover + B^n + S^n
+
+with `B^n` external-boundary momentum (inlet/outlet/walls) and `S^n` forcing.
+So **`F_CME + node_cover - node_uncover`** is, by construction, the body
+momentum that balances the fluid ledger of *this* storage/update order.
+Branch details: halfway fallback links (`f_r = f_t` plus the moving-wall
+addend) and q < 1/2 links enter the link accumulators exactly like Bouzidi
+links, with `c(f_t + f_r)`. In the q < 1/2 branch `G_dir(k_f - c)` is an
+interpolation input only, not a fluid-momentum sink. Pull-slot completeness
+is the section 3 contract. The scope excludes `ddf_shifting` (rejected), so
+the unshifted force form applies. The balance is a fluid-side identity; it
+does not close body+fluid globally in a domain with inlet/outlet/forcing
+(`B^n`, `S^n`), and the solid-slot writes are boundary plumbing, not physical
+streaming packets.
+
+The Wen 2014 relative-velocity link sum adds `-u_w (f_t - f_r) = u_w dm` per
+link (`dm = f_r - f_t`, the BFL mass defect). That term has no counterpart in
+the fluid ledger above. Peng et al. (PRE 2016 submitted, §II.C, Eqs. 20–21)
+show for a straight wall that the bounce-back error accumulated before a node
+switch cancels with the switching momentum. Whether `sum u_w dm` approximates
+`I_cover - I_uncover` here (implying double counting if both are used) is
+**NOT demonstrated** for this implementation: it needs a bounded experiment
+(prescribed slow translation; cumulative `sum u_w dm` vs cumulative node
+impulse, per crossing cycle). Momentum conservation is also distinct from
+Galilean invariance, which G4 tests.
+
+External reference note: waLBerla's covered/refill momentum switches
+(`conserveMomentumWhenMapping`, `conserveMomentumWhenRestoring`) default to
+false (Astra), so their existence does not establish which combination its
+validation cases use. A bookkeeping identity alone does not select the
+physical force either.
+
+Mass ledger (same proposition status). With `M^n` the active-fluid mass
+summed over `F^n` at the sample point,
+
+    M^n - M^{n-1} = sum_L dm - m_cover + m_uncover + B_m^n
+
+(`dm = f_r - f_t`, `m_cover`/`m_uncover` = summed `sum_q G` of covered and
+refilled cells, `B_m` external-boundary mass flux; collision conserves mass).
+
+Sample point (fixed): post-collision storage after the swap of step n, i.e.
+the same level `G` that M1 reads at step n+1; both ledgers are evaluated
+there. Ledger residual = left side minus every listed term. Zero mass flux
+does not imply zero momentum exchange: closed walls (bounce-back, free-slip)
+still exchange momentum, so "no inlet/outlet, no forcing" does NOT give
+`B = S = 0`. The rounding-level momentum check therefore uses either a
+momentum-isolated case (periodic in both directions, no forcing) or counts
+the wall impulse explicitly as part of `B`; the mass check may use a closed
+domain. Otherwise `B`, `S` are not separately measured, the residual is
+reported as their estimate and is not a pass criterion.
+Both ledger identities are **proposed, unreviewed**: they are not called
+proven until review and a host replay agree.
+
+Accumulators (separate, per body, per step window; all reported raw):
+`link_cme`, `link_wen_correction` (= `sum u_w dm`), `node_cover`,
+`node_uncover` (momentum), `link_dm`, `mass_cover`, `mass_uncover` (mass),
+plus counters of crossings, `eq_only` and `no_donor`, so force spikes can be
+correlated with crossing events.
+
+**OPEN — structural RHS candidates (not selected):**
+(a) `link_cme + node_cover - node_uncover`; (b) `link_cme + link_wen_correction`
+(no node impulse). Selection follows the ledger review and the bounded
+experiment above.
+
+## 7. Rejection semantics
+
+- Before any mutation (`ValueError`, no device write, host state unchanged):
+  moving body with EsoPull, `ddf_shifting`, refinement/blocked stepper, D3Q19,
+  Jets or surface actuation on that body, `sync_bodies`/checkpoint while a
+  moving body exists; swept band touching domain boundaries, other bodies,
+  sensors or force regions; `set_body` offset outside the declared box,
+  per-step move `|c_new - c_old|_2 > 0.5` cell (the single bound of 4a), speed
+above the declared bound.
+- Host motion state is committed only after the upload *and* the following
+  synchronization of the step group succeed (section 4); any failure in
+  between makes the simulation unusable and no host/device match is reported.
+- The only soft device diagnostic is `eq_only` (NEQ part omitted, counted).
+  `no_donor > 0` is a hard failure: it is checked after the step group and
+  before the step is committed or any results are returned; on failure the
+  simulation is unusable (section 5).
+
+## 8. Coupling and units (CANDIDATE)
+
+Lattice units, rho0 = 1, per unit span. Step n uses `c^n`, `u^n` and yields
+`F^n`. Weak staggered update
+`v^{n+1} = v^n + (F^n - c v^n - k y^n)/m`, `y^{n+1} = y^n + v^{n+1}`
+(Newmark is an alternative; OPEN). Mass, natural frequency, damping and `U_r`
+follow the chosen benchmark's own definitions. Leontini, Thompson & Hourigan,
+JFS 22:857–864 (2006), DOI 10.1016/j.jfluidstructs.2006.04.003, as read by
+Grok from the full text: Re = 200, `m* = m/m_f = 10` (m_f displaced fluid
+mass), `zeta = c/c_crit = 0.01`, `U* = U/(f_N D)` with `f_N` the natural
+frequency in fluid assuming an inviscid added-mass term; the paper also
+reports frequencies normalised by the in-vacuo value. Our runner derives `k`
+from that in-fluid `f_N` (added-mass coefficient as stated in the paper) and
+never from the vacuum formula; `c_crit` convention and the original
+comparison data are still OPEN. Stability/convergence is gated by fluid
+time/space refinement of the coupled run, not assumed from the mass ratio;
+structural sub-stepping does not substitute for fluid-coupling time
+refinement.
+
+## 9. Gates (map; tolerances OPEN, declared before runs)
+
+G0 zero-motion regression, split: G0a populations of the existing stationary
+and spin-only path bitwise unchanged vs main; G0b moving path at rest (with
+and without spin) reproduces the stationary link set and populations; force
+estimator outputs are compared separately from populations, since the moving
+path adds estimator terms · G1 raw force sign (`Cd > 0`) ·
+G2/G3 prescribed translation incl. ledger accumulators · G4 matched
+relative-flow (moving body in quiescent fluid vs fixed body in uniform
+stream), optional boost, co-moving sanity · G5 Dütsch et al. 1998 Re 100,
+KC 5 in-line (needs 2-component offset) · G6 Leontini free transverse VIV
+(m* = 10, ζ = 0.01, Re = 200, definitions pending) · G7 time-level contract
+frozen before the oscillator.
+
+Candidate pre-run engineering tolerances (Astra; for discussion, not
+literature limits or measured outcomes; norms, averaging window, convergence
+protocol and physical reference uncertainty must be fixed before adoption, and
+none may be weakened after seeing runs):
+normalised momentum-ledger residual FP32 ≤ 5e-6 (FP64 oracle ≤ 2e-13);
+10-cycle relative mass drift ≤ 1e-4; prescribed-motion grid study D = 40/80/160
+with change in the force-response amplitude (observed quantity and norm to be
+named before adoption) ≤ 2 % and phase ≤ 2°; enlarged domain
+amplitude change ≤ 1 %; matched-frame force RMS ≤ 2 % and improving with
+refinement; VIV refinement amplitude ≤ 2 %, frequency ≤ 1 %, five consecutive
+cycle amplitudes within 1 %.
+
+Conditions stated before any moving kernel (S2/S3), not only at the
+oscillator stage: (i) fluid refinement protocol for the prescribed-motion
+gates — diffusive scaling (dx halved, dt quartered, physical Re, KC and
+amplitude fixed), the observed quantities and norms named in advance;
+(ii) Galilean constant-boost G4-boost requires boundaries compatible with a
+uniform boost (periodic in the boost direction or boundaries that move with
+it); if the current dense configuration cannot represent that, G4-boost is
+recorded as not runnable rather than replaced by another test.
+
+## 10. Proposed minimal S1 scope (host only; held until cross-review)
+
+- `add_body` moving-circle declaration: base centre, declared offset box
+  (x/y), speed bound; swept-band computation and band-placement validation.
+- `set_body(offset=..., velocity=..., omega=...)` validation before mutation
+  (section 7), applied/pending state and the run-transition table of 4a,
+  and the fail-closed commit of section 4.
+- Host enumeration check of refill-donor non-emptiness (section 5) within
+  the declared `R_min` and displacement bound.
+- Motion-record host layout and upload path; exclusion of moving bodies from
+  the packed `CurvedBoundaryKernel` links (stationary packing unchanged).
+- Explicit rejections for unsupported modes.
+- S1 execution guard: until the S2 moving kernels exist, any `run(n)` with
+  `n > 0` (and `step(n)`) while a moving body is declared raises before any
+  device mutation (no action/motion upload, no launch), including the
+  nominal at-rest case, because the moving body is excluded from the static
+  packed links and nothing else would apply its boundary. `run(0)` keeps its
+  existing behaviour. Simulations without moving bodies are unaffected.
+  Preparation/commit helpers are exercised by host tests with mocks; no
+  operational unsupported stepping path is exposed.
+- Host tests (Astra runs them): rejection-before-mutation, band geometry,
+  derived centres, fail-closed behaviour under an injected upload/sync
+  failure, stationary packing byte-identical, execution guard (positive-step
+  run with a moving body raises before upload/launch; stationary unaffected). No kernels, no GPU runs, no
+  y-only VIV oscillator (S4), no force RHS.
+
+## 10a. S2a — fixed-centre moving links (implemented on `claude/moving-bodies-s2a`)
+
+- `MovingCircleLinks` (`kernels/step/moving_circle.cu`, dense double buffer)
+  runs before the static curved kernel when moving tables are attached: one
+  thread per (band cell, direction) slot, FP64 `{cx, cy, r}` at the applied
+  offset, q/hit/clamp/fallback in the order of `build_cut_links` with explicit
+  `_rn` operations (host/device bit parity is a test target, not a claim),
+  then the unchanged `apply_bouzidi_link` with spin wall velocity only.
+- `obs` keeps its contract (CME and torque from `apply_bouzidi_link`). A
+  separate per-body diagnostics buffer receives `cme`, the Wen term
+  `u_w * dm`, `link_dm` and the link count from the same link outputs; it is
+  zeroed with obs and only exists when moving circles do. Optional
+  `slot_debug` (q, fallback, rx, ry per slot) is null in production.
+- Displacement is disabled: `set_body(offset=..., vx=..., vy=...)` on a
+  translating circle raises before mutation, before and after initialize, so
+  no pending target can exist and `run`/`step`/direct `LBMStepper.step` all
+  hold at the applied geometry. `initialize()` never reads or consumes
+  `pending`. `stepper.blocked` is cleared only when the moving tables are
+  attached (initialize, stepper recreation).
+- S3 rule (not implemented): once targets are re-enabled, any pending target
+  must raise in `run`/`step` and in `LBMStepper.step` alike before any launch,
+  swap or step-count change; zero-step calls never consume it.
+- Not validation: a pass of the zero-motion gates says nothing about moving
+  geometry. GPU parity (G0a/G0b) is unverified until run on a free device.
+
+## 10b. S3a — translation with refill (implemented on `claude/moving-bodies-s3a`)
+
+- Targets are enabled again (4a). `set_body(offset, vx, vy)` validates against
+  `applied`; a pending target blocks direct `LBMStepper.step(n > 0)` before any
+  launch; `run(n > 1)` with a target raises before mutation; `run(0)` never
+  consumes it; `run(1)` advances it.
+- Per advancing step, on one stream: upload the immutable FP64 records
+  `{c_old, c_new, r, u}` (only when changed) → M1 `MovingCircleCells` → M2
+  `MovingCircleLinks` (`u_w = u + omega x r` at the hit point) → static curved →
+  `OneStep` → swap → force region → sensors → synchronize → read `no_donor` →
+  commit (`applied := pending.offset`) → upload the hold record
+  (`c_old = c_new = applied`, `u = 0`) → synchronize. Hold steps launch no M1
+  and are the S2a path. Commit point: the advancing step-group sync plus a
+  zero `no_donor` count; a failure before it commits nothing. The hold-record
+  upload stages the next step after the commit: if it (or its sync) fails,
+  `applied` already matches the moved device geometry, but the Simulation is
+  still marked failed and refuses reuse. Any failure from the first upload on
+  marks the circles and the Simulation failed.
+- M1 classifies with the old/new predicates only (no flag reads). Covered
+  cells add raw `sum c G` / `sum G` of `G^{n-1}` to `cover_*` and become
+  `SOLID|OBSTACLE|BC_CURVED`. Uncovered cells: donors are fluid in both
+  geometries; `rho_bar` over all of them; NEQ from the donor with the largest
+  unnormalized `c·n > 0` (first direction wins ties, `n` from `c_new`);
+  `u_w` (translation + spin) at the cell position; all nine slots written;
+  raw `sum c G`, `sum G` of the refill go to `uncover_*`; no donor →
+  `f_eq(RHO, u_w)` plus `no_donor` (hard failure before commit); no NEQ donor
+  → `eq_only`.
+- `get_flags()` reads the device flags back (synchronized) when moving tables
+  exist; otherwise it returns the host array as before.
+- Not selected: the structural force. CME+node vs CME+Wen stay OPEN until X1
+  and G4 have run.
+
+## 11. Remaining decisions
+
+Structural RHS (section 6); motion-record layout; refill details; coupling
+scheme; benchmark definitions and digitized reference data; tolerances;
+execution/resource feasibility (Astra, pending); whether a quiescent domain
+(zero inlet velocity) runs cleanly in the current config.

docs/supported_settings.md

@@ -35,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`,
@@ -84,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

docs/validation/duetsch_re100_kc5_protocol.md

@@ -0,0 +1,209 @@
+# Prescribed in-line oscillation: Dütsch Re 100 / KC 5 — PROTOCOL (review draft)
+
+Production under test: `claude/moving-bodies-s3a` at `653fd997` (unchanged).
+Primary source: Dütsch, Durst, Becker & Lienhart, J. Fluid Mech. 360, 249–271
+(1998), doi:10.1017/S002211209800860X. Source facts below were verified from
+the paper (Slack C0C78UFKGRM/1793335557.460509, TRON 1793356317–1793356724,
+page numbers given). Everything marked ENGINEERING is our a-priori choice,
+fixed before any benchmark result; it is not a paper uncertainty.
+
+## 0. Source facts (with provenance)
+
+- Motion (p.251): `x_c = -A sin(2 pi f t)`; paper phase phi = 0 is the mean
+  position moving left at `-U_max`. KC = 5 gives `A/D = 5/(2 pi) = .7957747`.
+- Start (p.254): the computation starts a quarter cycle before phi = 0, at
+  `x = +A` with zero velocity (phi = -90°).
+- Experiment (pp.252, 254–256): D = 10 mm, L = .65 m, tank 1.2 × .7 × .7 m,
+  rho = 998.2 kg/m³, mu = 1.00818 mPa·s, Re = 100, KC = 5; LDA with 1° phase
+  bins and 200 samples per node and degree; no quantified experimental error
+  recovered.
+- Profiles (Fig. 8, p.260): LAB-FIXED stations x1 = -6, 0, 6, 12 mm, i.e.
+  x/D = -.6, 0, .6, 1.2, at phi = 180°, 210°, 330°; u1 and u2 versus x2.
+- Original numerics (§2): cylinder-fixed accelerated frame (`x_i = x~_i +
+  x_i,s`, outer velocity `-u_c`, inertial source `-rho du_c/dt`), O-grid outer
+  diameter 120 D, illustrated 192 circumferential nodes, 720 steps/cycle; tank
+  influence neglected.
+- Force fit (verified): force ON THE BODY per unit span
+  `F1 = -0.5 rho D c_d xdot|xdot| - (pi/4) rho D² c_i xddot` with `x = -A sin phi`;
+  normalized by `0.5 rho D U_max²`: `c_d cos phi |cos phi| - (pi² c_i / KC) sin phi`.
+  `c_i` is the ADDED-MASS coefficient; total inertia `c_m = c_i + 1` is a
+  different convention and is not used.
+- Table 1 computed Morison coefficients (c_d, c_i) for 96 / 192 / 384 perimeter
+  nodes at dt/T = 1/720: (2.15, 1.42), (2.10, 1.45), (2.09, 1.45). These are
+  fitted computed coefficients, not a measured force and not an error bar. The
+  paper states the coarse direct force is closer to the finest direct force
+  than the Morison fit, but gives no waveform-error norm.
+
+## 1. Kinematics, time levels and phase bookkeeping
+
+- Our run is in the LAB frame: the circle moves, the fluid is initially at rest
+  (rho = 1, u = 0); no fictitious inertial source is added. Lab-fixed profile
+  stations are compared directly (no frame transform needed).
+- Geometry reference: the circle's base centre is the INITIAL position, the
+  right turning point `x_mean + A`; the declared offset box is x in [-2A, 0],
+  y = 0. The body is never teleported: offset 0 at s = 0.
+- Elapsed time s = n dt from the start. Analytic position relative to the mean
+  `x(s) = A cos(2 pi s / T)`, velocity `u(s) = -U_max sin(2 pi s / T)`; paper
+  phase `phi = (360 s / T - 90) mod 360`.
+- Step n moves from s_n to s_{n+1}: target offset `x(s_{n+1}) - A`; wall
+  velocity = the single-step interval average `(x(s_{n+1}) - x(s_n)) / dt`
+  (analytic, supplied, not differenced from stored offsets). The per-step
+  force of `run(1)` / `read_force(normalize=True)` is attributed to that same
+  interval.
+- Paper cycle k = 1, 2, ... spans `s in [T/4 + (k-1) T, T/4 + k T)`; the
+  initial quarter cycle (s < T/4) is cycle 0 and never scored. Every logged
+  sample carries (k, phi). Example: phi = 330° of cycle 1 is at s = 7T/6.
+- T is chosen as a multiple of 360 native steps, so phi = 180°, 210°, 330° fall
+  exactly on step boundaries and 1° bins hold whole steps.
+
+## 2. Resolution ladder (uniform grid, fixed Mach), ENGINEERING
+
+| D (cells) | T (steps) | U_max = 5 D / T | nu_lat = U_max D / 100 | tau (3 nu + .5) |
+|---|---|---|---|---|
+| 20 | 2,520 | .039683 | .0079365 | .5238 |
+| 40 | 5,040 | .039683 | .015873 | .5476 |
+| 80 | 10,080 | .039683 | .031746 | .5952 |
+
+- MRT collision. Mach fixed (Ma ≈ .0687); a separate Mach check at D = 40
+  with T = 10,080 (U_max halved, tau .5238) is labelled MACH and not part of the
+  ladder.
+- Per-step displacement ≤ U_max ≈ .04 cells (bound .5).
+
+## 3. Domain and boundaries (finite-domain reproduction), ENGINEERING
+
+The paper's 120 D O-grid is not reproduced; our rectangular box is a
+reproduction choice. The native dense channel (zero-velocity inlet,
+pressure-NEQ outlet, bounce-back y walls) is a finite-domain approximation,
+not infinite still fluid. The G4 x-periodic / y-slip override is a separately
+labelled sensitivity family and never mixed into the primary ladder. The
+runner owner confirms the zero-velocity inlet is supported as configured
+before any run.
+
+Order (closed):
+1. Smoke + ledger at D = 40, native BC, 40 D × 40 D (finite box).
+2. Domain sensitivity before any literature statement: same BC, D = 40,
+   L ∈ {20 D, 40 D, 60 D}; record the change in Fx waveform, c_d, c_i and the
+   scored profiles.
+3. Optional labelled BC sensitivity: one x-periodic / y-slip case at D = 40,
+   40 D; reported as a delta against native only.
+4. Only after step 2 passes: the D = 20 / 40 / 80 ladder on the chosen BC and
+   domain (default 40 D if, against 60 D, `RMS(dFx)/F_ref` and |dc_d|, |dc_i|
+   stay inside the section 6 bands).
+- Memory/runtime: measured by a smoke run before each new size; measurements
+  inform scheduling only, never tolerances.
+
+## 4. Observables (raw, per step, preserved)
+
+- `obs` Fx, Fy (CME) and all `MovingTables.DIAG` columns per step: CME,
+  CME + Wen and CME + cover − uncover reported separately; crossings, eq_only,
+  no_donor (no_donor must stay 0).
+- Velocity u1, u2 on lab-fixed lines x/D = -.6, 0, .6, 1.2 (relative to the
+  mean cylinder position) versus x2, at the scored phases. Lab frame
+  throughout; if any body-frame intermediate is used, position and velocity
+  are both transformed.
+- Sampling cost: raw force and diagnostics counters are read every step
+  (small). Flags and velocity fields are read only at the steps inside the
+  three prescribed 1° phase bins (φ ∈ [180, 181), [210, 211), [330, 331),
+  T/360 steps each) of the cycles needed for settling and scoring, plus
+  sparse declared audits; never every step.
+- Point set: fixed x2/D sample points per station, identical on every grid
+  and declared before running; per point the interpolation weights are
+  recorded; the scored set is the common subset valid on all compared grids
+  at that phase, with its coverage recorded; point sets never change silently.
+- Interpolation: bilinear on cell centres using only cells that are fluid at
+  that step (flags read back via `get_flags()`); weights renormalized over the
+  fluid corners; a sample with no fluid corner, or inside the circle at that
+  step, is MISSING and logged, never set to zero. One-degree bins keep valid
+  fluid weights. The numerical average need not mimic the 200 LDA samples.
+
+## 5. Periodic-regime gate (every case, before any statistic), ENGINEERING
+
+Ledger rows, same schema for every case: dx, dt (in D and T units); U_max, nu,
+Re, KC (lattice and nondimensional); T (native steps, nondimensional);
+duration (native steps, s/T, paper cycles); discard; accepted window; cycle
+ledger; verdict PASS or FAIL-CLOSED.
+
+Scales: `F_ref = 0.5 rho U_max² D`, `U_ref = U_max`; never a near-zero mean
+force as denominator.
+
+- Settling metrics, cycle k−1 → k, phase-aligned by step within the cycle:
+  - waveform `RMS(dFx) / F_ref ≤ .02`;
+  - `|d mean|`, `|d RMS|`, `|d peak|` of Fx, each `/ F_ref ≤ .01`;
+  - scored profiles (fluid-only) `|du| / U_max`, `|dv| / U_max ≤ .03` at the
+    three phases × four stations;
+  - Morison `|dc_d| ≤ .02`, `|dc_i| ≤ .02`.
+- Rule: all metrics for N = 3 consecutive whole paper cycles → settled.
+
+Cycle indexing (worked example, binding). Paper cycle k spans
+`s/T in [0.25 + (k-1), 0.25 + k)`, i.e. φ0 → φ0; cycle 0 is the startup
+quarter `s/T < 0.25`. Pair p_k compares cycle k−1 with cycle k.
+
+- Discarded: cycles 0, 1, 2 (s/T < 2.25). Pairs are formed only between
+  cycles ≥ 3, so the first pair is p_4 (3 → 4).
+- Settling proof for a case: three consecutive pairs all pass; K_case = the
+  first cycle AFTER the last of those pairs. Earliest possible: pairs 3→4,
+  4→5, 5→6 pass ⇒ K_case = 7.
+- Common window: K = max over compared cases of K_case; every case scores the
+  same three complete cycles K, K+1, K+2 and must stay settled throughout:
+  pairs (K−1→K), (K→K+1), (K+1→K+2) must pass on every case, otherwise that
+  case is FAIL-CLOSED (an earlier local settle does not waive this).
+- Cap: cycle k ends at s/T = k + 0.25, so the last admissible scored cycle is
+  11 (ends 11.25 T) and K ≤ 9. Examples: K = 7 → cycles 7–9 (ends 9.25 T);
+  K = 9 → cycles 9–11 (ends 11.25 T); any K_case > 9 → FAIL-CLOSED.
+- Every case is run to the end of cycle K + 2 (proportionally more native
+  steps on finer grids).
+- Discard: the initial quarter cycle plus the first 2 full paper cycles
+  (s/T < 2.25). Settling is searched from cycle 3; cap 12 periods from s = 0;
+  unsettled → FAIL-CLOSED (no literature, cross-resolution or VIV statistic).
+- Common window: K = the largest first-settled cycle index over all compared
+  cases; every case is scored on the same 3 cycles K .. K + 2 and must reach
+  them within the cap (proportionally more native steps on finer grids), or
+  the comparison is FAIL-CLOSED.
+- The earlier G4 timing check (100 coarse vs 400 diffusive-fine steps) was a
+  matched short transient, not periodic evidence.
+
+## 6. Metrics and per-observable criteria (ENGINEERING, not paper uncertainty)
+
+Three criteria stay separate; the G4 same-grid 2% target is not reused.
+
+- Morison fit: least squares of the normalized per-step force over the
+  accepted window against the two basis functions, each averaged over its
+  single step interval, fixed paper phase, no free phase shift; c_d and
+  added-mass c_i per force candidate.
+- Cross-resolution: every scored observable's error against its reference must
+  decrease monotonically D20 → D40 → D80 (necessary, not sufficient).
+- Against Table 1 (computed fit, not measured force): at D ≥ 40,
+  `|c_d - 2.09| ≤ .15` and `|c_i - 1.45| ≤ .15`, and monotone.
+- Force waveform against its own Morison reconstruction: residual reported
+  only, no pass/fail (the paper reports fit error above its coarse numerical
+  error and gives no percentage).
+- Fig. 8 LDA profiles: QUALITATIVE ONLY. The original figure could not be
+  accessed or rendered; no digitized ordinates and no measurement uncertainty
+  exist. After a digitization with figure number, tool and uncertainty note:
+  `RMS(|du|, |dv|) / U_max ≤ .15` at D ≥ 40 and monotone; LDA error stays
+  UNQUANTIFIED and is stated with every comparison.
+- Fy sanity (not literature): `|mean Fy| / F_ref ≤ .05` in the accepted window.
+- No force candidate is selected for an RHS by these results; each is
+  reported.
+
+Superseded (earlier drafts): the velocity-jump start, guessed stations and
+phases, "C_M" for c_i, and the uncertainty wording; cycles 4–6, the 30-cycle
+cap and the earlier tolerance numbers.
+
+## 6a. Driver unit tests (CPU, before any run)
+
+- Phase/cycle arithmetic: s ↔ (k, φ), cycle boundaries, bin membership,
+  window/cap selection for K = 7, 8, 9, 10.
+- Solid-aware interpolation: all-fluid, partial-solid (renormalized weights),
+  no-fluid (MISSING) and inside-circle cases.
+- Interval-averaged Morison bases and the two-parameter fit on a synthetic
+  waveform with known c_d, c_i.
+
+## 7. Order and scope
+
+Protocol review → one-shot Dütsch driver (analytic `set_body` + `run(1)`,
+cycle ledger, lab-fixed fluid-only station sampler, two-basis Morison fit; no
+reusable runner until the first smoke shows shared helpers pay off; no
+production change) → smoke (memory/runtime) → domain sensitivity →
+ladder → physical-validation verdict → only then single-cylinder VIV. No AMR,
+no pinball, no manuscript work.

examples/fixed_patch_circle.py

@@ -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
 
 
@@ -94,7 +94,7 @@ def main():
             cfg, context=context, coarse_boundary=args.mode, refinements=[BOUNDS],
             level_blocks=[tuple(args.base_block), tuple(args.fine_block)],
             level_objects=[(), (body,)]) as run:
-        p = run.patch
+        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]

examples/fixed_patch_jet.py

@@ -153,7 +153,7 @@ def main():
         with CudaContext(0) as context, FixedPatchStepper(cfg, context=context,
                 coarse_boundary='slot_reservoir', refinements=[BOUNDS],
                 level_blocks=LEVEL_BLOCKS) as run:
-            p = run.patch
+            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

@@ -50,7 +50,7 @@ def main():
                   store_precision="FP32", ddf_shifting=False, les_enabled=False)
     with CudaContext(0) as context, FixedPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)],
                                                            level_blocks=[(16, 16), (24, 24)]) as run:
-        p = run.patch
+        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/body/manager.py

@@ -405,7 +405,9 @@ class ObjectManager:
 
         for obj in objects:
             obj.validate_surface_actuation(self.cfg)
-            if (not obj.is_sensor) and (not obj.is_force_region) and hasattr(obj, 'get_curved_list'):
+            # Moving circles never enter the packed static links (docs/moving_bodies.md, section 3).
+            if (not obj.is_sensor) and (not obj.is_force_region) and obj.motion is None \
+                    and hasattr(obj, 'get_curved_list'):
                 (
                     fluid_idx, dirs, q_vals, body_ids,
                     rx_vals, ry_vals, rz_vals, fallback_vals

src/CelerisLab/body/motion.py

@@ -0,0 +1,259 @@
+# CelerisLab/body/motion.py
+"""Host state of a translating circle (dense double buffer, D2Q9 only).
+
+Contract: docs/moving_bodies.md, sections 4a, 5 and 7. The applied geometry
+(last successfully committed offset) and an optional pending target are two
+separate facts; centres are always ``base + offset``. Nothing here touches
+the device: the caller supplies upload/step/synchronize callables.
+"""
+from __future__ import annotations
+
+import math
+from dataclasses import dataclass
+from typing import Callable, Optional, Tuple
+
+MAX_STEP = 0.5     # Euclidean centre displacement per step (lattice cells)
+BAND_MARGIN = 2    # cells added around the swept circle AABB
+R_MIN = 4.0        # smallest radius covered by the donor enumeration
+
+Vec = Tuple[float, float]
+
+# D2Q9 non-rest directions, paired order of lbm/descriptors.
+_CX = (1, -1, 0, 0, 1, -1, 1, -1)
+_CY = (0, 0, 1, -1, 1, -1, -1, 1)
+
+
+def _vec(value, name: str) -> Vec:
+    x, y = (float(v) for v in value)
+    if not (math.isfinite(x) and math.isfinite(y)):
+        raise ValueError(f"{name} must be finite, got {value!r}")
+    return x, y
+
+
+@dataclass(frozen=True)
+class Translation:
+    """Declared motion envelope: offset box (x_min, x_max, y_min, y_max) and speed bound."""
+    offset_box: Tuple[float, float, float, float]
+    max_speed: float
+
+    def __post_init__(self):
+        box = tuple(float(v) for v in self.offset_box)
+        speed = float(self.max_speed)
+        if len(box) != 4 or not all(math.isfinite(v) for v in box + (speed,)):
+            raise ValueError("offset_box needs 4 finite values and max_speed must be finite")
+        if not (box[0] <= 0.0 <= box[1] and box[2] <= 0.0 <= box[3]):
+            raise ValueError("offset_box must contain the zero offset")
+        if not 0.0 < speed <= MAX_STEP:
+            raise ValueError(f"max_speed must satisfy 0 < max_speed <= {MAX_STEP}")
+        object.__setattr__(self, "offset_box", box)
+        object.__setattr__(self, "max_speed", speed)
+
+
+class MovingCircle:
+    """Applied offset, pending target and the fail-closed commit of one circle."""
+
+    def __init__(self, base, radius: float, translation: Translation,
+                 on_failure: Optional[Callable[[BaseException], None]] = None):
+        self.base = _vec(base, "base")
+        self.radius = float(radius)
+        if not (math.isfinite(self.radius) and self.radius >= R_MIN):
+            raise ValueError(f"moving circle radius must be >= {R_MIN}")
+        if not isinstance(translation, Translation):
+            raise TypeError("translation must be a Translation")
+        self.translation = translation
+        self.applied: Vec = (0.0, 0.0)
+        self.pending: Optional[Tuple[Vec, Vec]] = None
+        self.failed = False
+        self._on_failure = on_failure  # owner lifecycle hook (Simulation marks itself failed)
+
+    def centre(self, offset: Vec) -> Vec:
+        return self.base[0] + offset[0], self.base[1] + offset[1]
+
+    def check_target(self, offset, velocity) -> Tuple[Vec, Vec]:
+        """Validate a target against ``applied``; returns it, mutates nothing."""
+        self._require_usable()
+        offset, velocity = _vec(offset, "offset"), _vec(velocity, "velocity")
+        x0, x1, y0, y1 = self.translation.offset_box
+        if not (x0 <= offset[0] <= x1 and y0 <= offset[1] <= y1):
+            raise ValueError(f"offset {offset} outside declared box {self.translation.offset_box}")
+        step = math.hypot(offset[0] - self.applied[0], offset[1] - self.applied[1])
+        if step > MAX_STEP:
+            raise ValueError(f"per-step displacement {step} exceeds {MAX_STEP}")
+        if math.hypot(*velocity) > self.translation.max_speed:
+            raise ValueError(f"velocity {velocity} exceeds max_speed {self.translation.max_speed}")
+        if offset == self.applied and velocity != (0.0, 0.0):
+            raise ValueError("a hold target (offset == applied) needs zero translational velocity")
+        return offset, velocity
+
+    def set_target(self, offset, velocity) -> None:
+        """Replace the pending target (last wins); invalid input keeps the old one."""
+        self.pending = self.check_target(offset, velocity)
+
+    def record(self) -> Tuple[Vec, Vec, Vec]:
+        """Motion record (c_old, c_new, u) for the next single moving step."""
+        self._require_usable()
+        if self.pending is not None:
+            offset, velocity = self.pending
+        else:
+            offset, velocity = self.applied, (0.0, 0.0)
+        return self.centre(self.applied), self.centre(offset), velocity
+
+    def advance(self, upload: Callable, step: Callable, synchronize: Callable) -> None:
+        """Upload the record, run one step and commit only after a successful sync.
+
+        Any exception after the upload is attempted leaves device geometry
+        unknown: the circle is marked failed, never reports applied state again,
+        and the owner's ``on_failure`` hook is called.
+        """
+        record = self.record()
+        try:
+            upload(record)
+            step()
+            synchronize()
+        except BaseException as exc:
+            self.fail(exc)
+            raise
+        self.commit()
+
+    def commit(self) -> None:
+        """Applied := pending offset (caller has synchronized a successful step)."""
+        if self.pending is not None:
+            self.applied = self.pending[0]
+            self.pending = None
+
+    def fail(self, exc: BaseException) -> None:
+        self.failed = True
+        if self._on_failure is not None:
+            self._on_failure(exc)
+
+    def band(self) -> Tuple[int, int, int, int]:
+        """Inclusive cell AABB (x0, x1, y0, y1) swept over the whole offset box."""
+        x0, x1, y0, y1 = self.translation.offset_box
+        r, m = self.radius, BAND_MARGIN
+        return (math.floor(self.base[0] + x0 - r) - m, math.ceil(self.base[0] + x1 + r) + m,
+                math.floor(self.base[1] + y0 - r) - m, math.ceil(self.base[1] + y1 + r) + m)
+
+    def _require_usable(self) -> None:
+        if self.failed:
+            raise RuntimeError("moving circle state is unknown after a failed step")
+
+
+def _inside(x: int, y: int, c: Vec, r: float) -> bool:
+    return (x - c[0]) ** 2 + (y - c[1]) ** 2 < r * r
+
+
+def empty_donor_cells(c_old: Vec, c_new: Vec, radius: float):
+    """Fresh cells (inside c_old, outside c_new) whose donor set V(x) is empty."""
+    r = radius
+    xs = range(math.floor(min(c_old[0], c_new[0]) - r) - 1, math.ceil(max(c_old[0], c_new[0]) + r) + 2)
+    ys = range(math.floor(min(c_old[1], c_new[1]) - r) - 1, math.ceil(max(c_old[1], c_new[1]) + r) + 2)
+    empty = []
+    for x in xs:
+        for y in ys:
+            if not _inside(x, y, c_old, r) or _inside(x, y, c_new, r):
+                continue
+            if not any(not _inside(x + cx, y + cy, c_old, r) and not _inside(x + cx, y + cy, c_new, r)
+                       for cx, cy in zip(_CX, _CY)):
+                empty.append((x, y))
+    return empty
+
+
+def enumerate_donor_coverage(radius: float, positions: int = 8, directions: int = 16) -> int:
+    """Count empty-donor fresh cells over sub-cell centres and moves of MAX_STEP.
+
+    Evidence for the declared bounds (radius >= R_MIN, |dc|_2 <= MAX_STEP),
+    not a proof; the device keeps the ``no_donor`` counter.
+    """
+    count = 0
+    for i in range(positions):
+        for j in range(positions):
+            c_old = (i / positions, j / positions)
+            for k in range(directions):
+                a = 2.0 * math.pi * k / directions
+                c_new = (c_old[0] + MAX_STEP * math.cos(a), c_old[1] + MAX_STEP * math.sin(a))
+                count += len(empty_donor_cells(c_old, c_new, radius))
+    return count
+
+
+class MovingTables:
+    """Device tables of the translating circles for ``MovingCircleCells``/``Links``.
+
+    Band table (x0, y0, width, height, first_slot, body_id), immutable per run;
+    FP64 motion records (cx_old, cy_old, cx_new, cy_new, r, ux, uy), immutable
+    during a step; separate per-body diagnostics (see DIAG), never part of obs.
+    """
+    BAND_INTS, GEOM_DOUBLES = 6, 7
+    DIAG = ("cme_x", "cme_y", "wen_x", "wen_y", "link_dm", "n_links",
+            "cover_jx", "cover_jy", "cover_mass", "uncover_jx", "uncover_jy", "uncover_mass",
+            "eq_only", "no_donor", "crossings")
+    DIAG_FLOATS = len(DIAG)
+
+    def __init__(self, band, geom):
+        import numpy as np
+        import pycuda.driver as cuda
+        self.n_bodies = len(geom) // self.GEOM_DOUBLES
+        self.n_slots = int(band[-2] + band[-4] * band[-3] * 8)
+        self.n_cells = self.n_slots // 8
+        self.geom = list(geom)  # host copy of the uploaded records
+        self.band_gpu = cuda.to_device(np.asarray(band, np.int32))
+        self.geom_gpu = cuda.to_device(np.asarray(geom, np.float64))
+        # Owned page-locked staging for async record uploads (allocated once).
+        self._geom_pinned = cuda.pagelocked_empty(len(geom), np.float64)
+        self.diag_gpu = cuda.to_device(np.zeros(self.n_bodies * self.DIAG_FLOATS, np.float32))
+        self.slot_debug_gpu = None  # optional (q, fallback, rx, ry) per slot for parity tests
+
+    @classmethod
+    def upload(cls, circles):
+        """``circles``: [(MovingCircle, body_id)] in body order."""
+        return cls(*cls.layout(circles))
+
+    @staticmethod
+    def layout(circles):
+        """Host band and geometry tables (no device work)."""
+        band, geom, first = [], [], 0
+        for motion, body_id in circles:
+            x0, x1, y0, y1 = motion.band()
+            width, height = x1 - x0 + 1, y1 - y0 + 1
+            band += [x0, y0, width, height, first, int(body_id)]
+            first += width * height * 8
+            geom += MovingTables.record(motion)
+        return band, geom
+
+    @staticmethod
+    def record(motion):
+        """FP64 motion record of one circle for its next step (hold if nothing pending)."""
+        c_old, c_new, u = motion.record()
+        return [*c_old, *c_new, motion.radius, *u]
+
+    def upload_records(self, circles, stream) -> bool:
+        """Enqueue the records when they differ from the device copy; True if advancing.
+
+        The pinned staging buffer is reused: callers synchronize ``stream`` after every
+        upload before the next one (``Simulation._run_moving``), so no copy is in flight
+        when it is rewritten. ``geom`` is the enqueued state; after a failed copy the
+        owner is failed and the cache is never trusted again.
+        """
+        import pycuda.driver as cuda
+        geom = [v for motion, _ in circles for v in self.record(motion)]
+        if geom != self.geom:
+            self._geom_pinned[:] = geom
+            cuda.memcpy_htod_async(self.geom_gpu, self._geom_pinned, stream)
+            self.geom = geom
+        g = self.GEOM_DOUBLES
+        return any(geom[i:i + 2] != geom[i + 2:i + 4] for i in range(0, len(geom), g))
+
+    def enable_slot_debug(self):
+        import numpy as np
+        import pycuda.driver as cuda
+        self.slot_debug_gpu = cuda.to_device(np.zeros(self.n_slots * 4, np.float32))
+
+    def zero_diagnostics(self, stream):
+        import pycuda.driver as cuda
+        cuda.memset_d32_async(self.diag_gpu, 0, self.n_bodies * self.DIAG_FLOATS, stream)
+
+    def download_diagnostics(self):
+        import numpy as np
+        import pycuda.driver as cuda
+        out = np.empty((self.n_bodies, self.DIAG_FLOATS), np.float32)
+        cuda.memcpy_dtoh(out, self.diag_gpu)
+        return out

src/CelerisLab/body/objects.py

@@ -140,6 +140,7 @@ class SimObject:
             x=float(c[0]), y=float(c[1]), z=z0,
         )
         self.control = ObjectControl()
+        self.motion = None  # MovingCircle for a translating circle (body/motion.py)
 
     @property
     def surface_actuation(self):

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/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/moving_circle.cu

@@ -0,0 +1,183 @@
+// Translating circles (docs/moving_bodies.md): M1 MovingCircleCells reclassifies
+// band cells and refills uncovered ones; M2 MovingCircleLinks applies the cut links.
+// Geometry comes from an immutable FP64 record per body
+// {cx_old, cy_old, cx_new, cy_new, r, ux, uy}; the M2 operation order mirrors
+// CircleGeometry.build_cut_links / find_circle_intersection with explicit _rn
+// operations (bit parity with the host builder is a test target, not a claim).
+// Moving circles are never in the packed static links, so this kernel is the
+// only writer of their solid slots.
+#ifndef CELERIS_MOVING_CIRCLE_CU
+#define CELERIS_MOVING_CIRCLE_CU
+
+#if DIM == 2
+constexpr int MOVING_BAND_INTS = 6;   // x0, y0, width, height, first_slot, body_id
+constexpr int MOVING_GEOM_DOUBLES = 7;
+// cme_x, cme_y, wen_x, wen_y, link_dm, n_links, cover_jx, cover_jy, cover_mass,
+// uncover_jx, uncover_jy, uncover_mass, eq_only, no_donor, crossings
+constexpr int MOVING_DIAG_FLOATS = 15;
+
+__device__ __forceinline__ double mc_sq(double v) { return __dmul_rn(v, v); }
+
+__device__ __forceinline__ bool mc_inside(int x, int y, double cx, double cy, double r)
+{
+    return __dadd_rn(mc_sq(__dsub_rn((double)x, cx)), mc_sq(__dsub_rn((double)y, cy))) < mc_sq(r);
+}
+
+// Smallest root t in [0, 1] of |p + t c - centre| = r (find_circle_intersection order).
+__device__ __forceinline__ bool mc_hit(int x, int y, int ex, int ey,
+                                       double cx, double cy, double r, double& t)
+{
+    const double dx = (double)ex, dy = (double)ey;
+    const double px = __dsub_rn((double)x, cx), py = __dsub_rn((double)y, cy);
+    const double a = __dadd_rn(mc_sq(dx), mc_sq(dy));
+    const double b = __dmul_rn(2.0, __dadd_rn(__dmul_rn(dx, px), __dmul_rn(dy, py)));
+    const double c = __dsub_rn(__dadd_rn(mc_sq(px), mc_sq(py)), mc_sq(r));
+    const double det = __dsub_rn(mc_sq(b), __dmul_rn(__dmul_rn(4.0, a), c));
+    if (det < 0.0) return false;
+    const double s = __dsqrt_rn(det), a2 = __dmul_rn(2.0, a);
+    const double t1 = __ddiv_rn(__dsub_rn(-b, s), a2);
+    const double t2 = __ddiv_rn(__dadd_rn(-b, s), a2);
+    if (0.0 <= t1 && t1 <= 1.0) { t = t1; return true; }
+    if (0.0 <= t2 && t2 <= 1.0) { t = t2; return true; }
+    return false;
+}
+
+__global__ void MovingCircleLinks(
+    fpxx* fi, const uint16_t* flag,
+    const int* band, const double* geom,
+    unsigned int n_bodies, unsigned int n_slots,
+    const float* action, float* obs, float* diag, float* slot_debug)
+{
+    const unsigned int tid = threadIdx.x + blockIdx.x * blockDim.x;
+    if (tid >= n_slots) return;
+    unsigned int m = 0;
+    while (m + 1 < n_bodies && tid >= (unsigned int)band[(m + 1) * MOVING_BAND_INTS + 4]) ++m;
+    const int* b = band + m * MOVING_BAND_INTS;
+    const unsigned int local = tid - (unsigned int)b[4];
+    const unsigned int dir = local % 8u + 1u;
+    const int x = b[0] + (int)((local / 8u) % (unsigned int)b[2]);
+    const int y = b[1] + (int)((local / 8u) / (unsigned int)b[2]);
+    const double* g = geom + MOVING_GEOM_DOUBLES * m;
+    const double cx = g[2], cy = g[3], r = g[4];
+    const int ex = d_cx[dir], ey = d_cy[dir];
+    float* debug = slot_debug == nullptr ? nullptr : slot_debug + 4 * tid;
+    if (debug != nullptr) debug[0] = -1.0f;
+
+    double t;
+    if (mc_inside(x, y, cx, cy, r) || !mc_inside(x + ex, y + ey, cx, cy, r)
+            || !mc_hit(x, y, ex, ey, cx, cy, r, t))
+        return;
+    const double hx = __dadd_rn((double)x, __dmul_rn(t, (double)ex));
+    const double hy = __dadd_rn((double)y, __dmul_rn(t, (double)ey));
+    const double hdx = __dsub_rn(hx, (double)x), hdy = __dsub_rn(hy, (double)y);
+    const double norm = __dsqrt_rn((double)(ex * ex + ey * ey));
+    double qd = __ddiv_rn(__dsqrt_rn(__dadd_rn(mc_sq(hdx), mc_sq(hdy))), norm);
+    qd = fmax(0.01, fmin(qd, 1.0));
+    const float q = (float)qd;
+    const float rx = (float)__dsub_rn(hx, cx), ry = (float)__dsub_rn(hy, cy);
+
+    // Same two-stage order as build_cut_links + pack_cut_links_to_soa: the fallback
+    // class is chosen from the clamped FP64 q (q < 1/2 needs a fluid donor at x - c);
+    // only then are q, rx, ry rounded to FP32 for the BFL branch and spin velocity.
+    unsigned char fallback = CURVED_FALLBACK_BOUZIDI;
+    if (qd < 0.5) {
+        const int xd = x - ex, yd = y - ey;
+        if (xd < 0 || xd >= (int)NX || yd < 0 || yd >= (int)NY || mc_inside(xd, yd, cx, cy, r)
+                || (flag[linear_index((unsigned int)xd, (unsigned int)yd)] & FLAG_FLUID) == 0)
+            fallback = CURVED_FALLBACK_HALFWAY;
+    }
+    if (debug != nullptr) { debug[0] = q; debug[1] = (float)fallback; debug[2] = rx; debug[3] = ry; }
+
+    const int bid = b[5];
+    const float omega = action_omega(action, bid);
+    const float Uw = (float)g[5] - omega * ry, Vw = (float)g[6] + omega * rx;  // u + omega x r
+    float exchange[3], dm;
+    apply_bouzidi_link(dir, q, linear_index((unsigned int)x, (unsigned int)y), fi,
+                       Uw, Vw, rx, ry, fallback, obs, bid, exchange, &dm);
+    if (diag != nullptr) {
+        float* d = diag + m * MOVING_DIAG_FLOATS;
+        atomicAdd(&d[0], exchange[0]);
+        atomicAdd(&d[1], exchange[1]);
+        atomicAdd(&d[2], Uw * dm);   // Wen 2014 correction -u_w (f_t - f_r)
+        atomicAdd(&d[3], Vw * dm);
+        atomicAdd(&d[4], dm);
+        atomicAdd(&d[5], 1.0f);
+    }
+}
+// M1: one thread per band cell, launched only on steps with c_old != c_new.
+// Membership comes from the immutable old/new predicates, never from flags, so
+// classification, refill and the flag write of one cell cannot race with another.
+// ddf holds post-collision G^{n-1}; refilled cells get all nine slots at that level.
+__global__ void MovingCircleCells(
+    fpxx* fi, uint16_t* flag, const int* band, const double* geom,
+    unsigned int n_bodies, unsigned int n_cells, const float* action, float* diag)
+{
+    const unsigned int tid = threadIdx.x + blockIdx.x * blockDim.x;
+    if (tid >= n_cells) return;
+    unsigned int m = 0;
+    while (m + 1 < n_bodies && tid >= (unsigned int)band[(m + 1) * MOVING_BAND_INTS + 4] / 8u) ++m;
+    const int* b = band + m * MOVING_BAND_INTS;
+    const unsigned int local = tid - (unsigned int)b[4] / 8u;
+    const int x = b[0] + (int)(local % (unsigned int)b[2]);
+    const int y = b[1] + (int)(local / (unsigned int)b[2]);
+    const double* g = geom + MOVING_GEOM_DOUBLES * m;
+    const bool in_old = mc_inside(x, y, g[0], g[1], g[4]);
+    const bool in_new = mc_inside(x, y, g[2], g[3], g[4]);
+    if (in_old == in_new) return;
+    const unsigned long k = linear_index((unsigned int)x, (unsigned int)y);
+    float* d = diag + m * MOVING_DIAG_FLOATS;
+    float rho = 0.0f, jx = 0.0f, jy = 0.0f;
+    if (in_new) {  // covered: its fluid content leaves the active set
+        for (unsigned int i = 0; i < NQ; ++i) {
+            const float f = load_ddf(fi, index_f(k, i));
+            rho += f; jx += (float)d_cx[i] * f; jy += (float)d_cy[i] * f;
+        }
+        flag[k] = FLAG_SOLID | FLAG_OBSTACLE | FLAG_BC_CURVED;
+        atomicAdd(&d[6], jx); atomicAdd(&d[7], jy); atomicAdd(&d[8], rho);
+        atomicAdd(&d[14], 1.0f);
+        return;
+    }
+    // Uncovered: donors are fluid in BOTH geometries; NEQ from the neighbour with the
+    // largest unnormalized c.n > 0 (first direction wins ties), n from c_new.
+    const float nx = (float)__dsub_rn((double)x, g[2]), ny = (float)__dsub_rn((double)y, g[3]);
+    float rho_sum = 0.0f, best_dot = 0.0f;
+    int n_valid = 0, best = -1;
+    for (unsigned int i = 1; i < NQ; ++i) {
+        const int xn = x + d_cx[i], yn = y + d_cy[i];
+        if (mc_inside(xn, yn, g[0], g[1], g[4]) || mc_inside(xn, yn, g[2], g[3], g[4])) continue;
+        const unsigned long kn = linear_index((unsigned int)xn, (unsigned int)yn);
+        for (unsigned int q = 0; q < NQ; ++q) rho_sum += load_ddf(fi, index_f(kn, q));
+        ++n_valid;
+        const float dot = (float)d_cx[i] * nx + (float)d_cy[i] * ny;
+        if (dot > best_dot) { best_dot = dot; best = (int)i; }
+    }
+    const float omega = action_omega(action, b[5]);
+    const float uwx = (float)g[5] - omega * ny, uwy = (float)g[6] + omega * nx;  // at the cell
+    float feq[NQ], neq[NQ];
+    compute_feq(n_valid > 0 ? rho_sum / (float)n_valid : (float)RHO, uwx, uwy, feq);
+    for (unsigned int i = 0; i < NQ; ++i) neq[i] = 0.0f;
+    if (best >= 0) {
+        const unsigned long kd = linear_index((unsigned int)(x + d_cx[best]), (unsigned int)(y + d_cy[best]));
+        float gd[NQ], rd = 0.0f, mx = 0.0f, my = 0.0f, fd[NQ];
+        for (unsigned int i = 0; i < NQ; ++i) {
+            gd[i] = load_ddf(fi, index_f(kd, i));
+            rd += gd[i]; mx += (float)d_cx[i] * gd[i]; my += (float)d_cy[i] * gd[i];
+        }
+        compute_feq(rd, mx / rd, my / rd, fd);
+        for (unsigned int i = 0; i < NQ; ++i) neq[i] = gd[i] - fd[i];
+    } else {
+        atomicAdd(&d[12], 1.0f);
+    }
+    if (n_valid == 0) atomicAdd(&d[13], 1.0f);  // hard failure, checked before commit
+    for (unsigned int i = 0; i < NQ; ++i) {
+        const float f = feq[i] + neq[i];
+        store_ddf(fi, index_f(k, i), f);
+        rho += f; jx += (float)d_cx[i] * f; jy += (float)d_cy[i] * f;
+    }
+    flag[k] = FLAG_FLUID;
+    atomicAdd(&d[9], jx); atomicAdd(&d[10], jy); atomicAdd(&d[11], rho);
+    atomicAdd(&d[14], 1.0f);
+}
+#endif  // DIM == 2
+
+#endif  // CELERIS_MOVING_CIRCLE_CU

src/CelerisLab/lbm/kernels/uniform.cu

@@ -71,6 +71,7 @@ extern "C"
 #include "step/one_step_double.cu"
 #include "step/one_step_esopull.cu"
 #include "step/aux_kernels.cu"
+#include "step/moving_circle.cu"
 #include "step/esopull_macro.cu"
 #include "step/esopull_sync.cu"
 

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
@@ -21,7 +22,7 @@ from .channel import channel_masks, validate_channel_config, validate_channel_tr
 from ..body.jets import jet_link_velocity
 from .hierarchy import build_hierarchy
 from .refinement import ACTIVE, HALO, FixedPatch2D, FixedPatchSet2D, prepare_fine_boundary
-from .refinement_tables import (blocked_base_tables, blocked_fine_tables,  # noqa: F401 (re-export)
+from .refinement_tables import (blocked_coarse_tables, blocked_fine_tables,
                                 blocked_kind, projected_explosion_tables, _prepare_blocked_boundary,
                                 _BlockedBoundary,
                                 BLOCKED_BOUNDARY_DEVICE_COLUMNS)
@@ -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)
@@ -92,6 +111,57 @@ def _check_levels(refinements, blocks, objects, jets):
             raise ValueError("bodies/Jets are supported only on the finest level of a nested hierarchy")
 
 
+# Host preparation (numpy only, before any device work); each returns one entry per level.
+
+def _level_layouts(cfg, pairs, level_blocks, wall_aware):
+    """Blocked storage layout of every level, checked against the launch limits."""
+    layouts = []
+    for level, block in enumerate(level_blocks):
+        if level == 0:
+            layout = BlockLayout2D.build(cfg.ny, cfg.nx, *block, periodic=(True, not wall_aware))
+        elif isinstance(pairs[level - 1], FixedPatchSet2D):
+            layout = pairs[level - 1].fine_layout(block)
+        else:
+            layout = BlockLayout2D.build(*pairs[level - 1].fine_shape, *block,
+                                         periodic=(False, False))
+        if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM:
+            raise ValueError(f"level {level} blocks exceed launch limits")
+        layouts.append(layout)
+    return layouts
+
+
+def _host_boundaries(hierarchy, level_cfgs, level_objects, level_jets, wall_aware, shifted):
+    """Host boundary geometry of every refined level in the root edge frame (None at L0)."""
+    hosts = [None]
+    for level in range(1, hierarchy.depth):
+        pair = hierarchy.pairs[level - 1]
+        if level == 1:  # the root config carries the wall-aware boundary fields
+            host = prepare_fine_boundary(level_cfgs[0], pair, level_objects[1],
+                                         wall_aware=wall_aware, jets=level_jets[1],
+                                         _shifted_blocked=shifted)
+        else:
+            origin, spacing = hierarchy.edge_frame(level - 1)
+            host = prepare_fine_boundary(level_cfgs[level - 1], pair, level_objects[level],
+                                         jets=level_jets[level], parent_origin=origin,
+                                         parent_spacing=spacing, _shifted_blocked=shifted)
+        hosts.append(host)
+    return hosts
+
+
+def _fluid_leaves(hierarchy, boundaries):
+    """ACTIVE cells of every level that no static circle or body makes solid."""
+    leaves = [kind == ACTIVE for kind in hierarchy.kinds]
+    for level, b in enumerate(boundaries):
+        if b is None:
+            continue
+        if b.host is None:  # single level: static root circles
+            leaves[level] &= ~b.tables["solid"]
+        else:
+            fine_shape = hierarchy.pairs[level - 1].fine_shape
+            leaves[level] &= (b.host.flags.reshape(fine_shape) & SOLID) == 0
+    return leaves
+
+
 class FixedPatchStepper:
     """D2Q9 SRT/MRT FP32 fixed 2:1 refinement hierarchy, built from level lists.
 
@@ -104,7 +174,8 @@ class FixedPatchStepper:
     supported).
     ``refinements=[]`` is one blocked level: periodic SRT root without LES,
     Jets or exchange observation; its ``level_objects[0]`` are stationary
-    passive circles (static BFL).
+    passive circles (static BFL) in the integer-NODE frame, like Simulation
+    objects (see Hierarchy for frames).
 
     ``coarse_boundary`` is ``periodic`` (default) or ``slot_reservoir``:
     a bottom on-site slot and top fixed-density reservoir from ``cfg``.
@@ -116,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
@@ -163,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
@@ -219,56 +297,15 @@ class FixedPatchStepper:
         self.hierarchy = build_hierarchy(cfg, refinements, wall_aware=self._wall_aware)
         pairs = self.hierarchy.pairs
         if self._wall_aware:
-            rects = (pairs[0].bounds,) if isinstance(pairs[0], FixedPatch2D) else pairs[0].bounds
-            if any(r[1] < 4 or r[3] > cfg.ny - 4 for r in rects):
+            if any(r[1] < 4 or r[3] > cfg.ny - 4 for r in pairs[0].rectangles):
                 raise ValueError("slot_reservoir requires y0 >= 4 and y1 <= ny-4")
-        self._layouts = []
-        for level, block in enumerate(level_blocks):
-            if level == 0:
-                layout = BlockLayout2D.build(cfg.ny, cfg.nx, *block,
-                                             periodic=(True, not self._wall_aware))
-            elif isinstance(pairs[level - 1], FixedPatchSet2D):
-                layout = pairs[level - 1].fine_layout(block)
-            else:
-                layout = BlockLayout2D.build(*pairs[level - 1].fine_shape, *block,
-                                             periodic=(False, False))
-            if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM:
-                raise ValueError(f"{_level_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")
         self._observe_raw_exchange = observe_raw_exchange
@@ -280,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"
@@ -325,21 +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 level.boundary is not None:
+                        self._upload_boundary(level)
             except BaseException:
                 try:
                     self.close()
@@ -347,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).
 
@@ -357,7 +390,7 @@ class FixedPatchStepper:
         for index, layout in enumerate(self._layouts):
             total += 2 * 4 * int(np.prod(layout.blocked_shape))
             total += layout.halo_dst.nbytes + layout.halo_src.nbytes
-            b = self._level_boundaries[index]
+            b = self._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)
@@ -380,19 +413,19 @@ class FixedPatchStepper:
         grid = ((config.nx + self._block[0] - 1) // self._block[0], config.ny, layout.nblocks)
         level = _Level(shape, module, self._allocate(size), self._allocate(size),
                        module.get_function(name), grid, layout)
-        level.boundary = self._level_boundaries[index]
+        level.boundary = self._boundaries[index]
         level.halo_copy = module.get_function("BlockHaloCopy")
         if len(layout.halo_dst):
             level.halo_dst = self._upload(layout.halo_dst)
             level.halo_src = self._upload(layout.halo_src)
         if index == 0 and self._channel:
-            cfg = self._remap_cfg
+            cfg = self._root_cfg
             ptr, size = module.get_global("d_block_channel")
             if size != struct.calcsize("=IIf"):
                 raise RuntimeError("Block channel host/device layout mismatch")
             cuda.memcpy_htod(ptr, struct.pack("=IIf", cfg.nx, cfg.ny, cfg.velocity))
         if index == 0 and self._wall_aware:
-            cfg = self._remap_cfg
+            cfg = self._root_cfg
             ptr, size = module.get_global("d_block_wall")
             if size != struct.calcsize("=IIIf"):
                 raise RuntimeError("Block wall host/device layout mismatch")
@@ -406,11 +439,11 @@ class FixedPatchStepper:
         order = np.argsort(pair.coalesce_coarse)
         ghosts, children, child_kind = blocked_fine_tables(pair, child.layout, order)
         if index == 0:
-            donors, targets, parent_kind = blocked_base_tables(pair, parent.layout, order,
+            donors, targets, parent_kind = blocked_coarse_tables(pair, parent.layout, order,
                                                                self._wall_aware)
         else:
             # Nested pairs keep inherited INNER/OUTER and sparse holes of the parent.
-            donors, targets, parent_kind = blocked_base_tables(
+            donors, targets, parent_kind = blocked_coarse_tables(
                 pair, parent.layout, order, composed_kind=self.hierarchy.kinds[index])
         result = _Pair(child.module.get_function("BlockExplode"),
                        parent.module.get_function("BlockCoalesce"),
@@ -516,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.
@@ -541,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:
@@ -561,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()
@@ -588,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:
@@ -606,10 +638,9 @@ class FixedPatchStepper:
                 raise
             return candidate
 
-    def set_jet_amplitude(self, index, amplitude, *, level=1):
-        """Replace one jet's amplitude between steps; no geometry/table rebuild.
+    def set_jet_amplitude(self, index, amplitude, *, level):
+        """Replace Jet ``index`` of level ``level`` (``jets[level]``) between steps; no 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.
@@ -618,28 +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 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
@@ -682,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)
@@ -726,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)
@@ -749,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:
@@ -758,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")
@@ -769,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
@@ -802,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"])
@@ -818,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."""
@@ -862,7 +887,7 @@ 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
@@ -889,24 +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)
+
+    @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/lbm/stepper.py

@@ -20,6 +20,10 @@ class LBMStepper:
         self.cfg = cfg
 
         self._esopull = (cfg.streaming == "esopull")
+        # Set by the owner when a step would be wrong (translating circles
+        # excluded from the packed links before the moving kernels exist).
+        self.blocked: str | None = None
+        self.moving = None  # MovingTables attached by the owner (translating circles)
 
         if self._esopull:
             self.step_fn = module.get_function("EsoPullStep")
@@ -35,6 +39,8 @@ class LBMStepper:
         self.apply_esopull_curved_fn = module.get_function("ApplyEsoPullCurvedKernel")
         self.sensor_fn = module.get_function("SensorKernel")
         self.force_region_fn = module.get_function("ForceRegionKernel")
+        self.moving_fn = None if self._esopull else module.get_function("MovingCircleLinks")
+        self.moving_cells_fn = None if self._esopull else module.get_function("MovingCircleCells")
 
         tpb = cfg.threads_per_block
         self.block = (tpb, 1, 1)
@@ -70,6 +76,13 @@ class LBMStepper:
         stream: cuda.Stream | None = None,
     ):
         """Advance *n* time steps."""
+        if n > 0 and self.blocked is not None:
+            raise RuntimeError(self.blocked)
+        self._step(n, action_gpu=action_gpu, obs_gpu=obs_gpu, stream=stream)
+
+    def _step(self, n, *, action_gpu, obs_gpu, stream=None, advance=False):
+        """Unguarded steps; ``advance`` runs one step that moves translating circles."""
+        assert not advance or n == 1
         f = self.field
         launch_kw = {}
         if stream is not None:
@@ -81,6 +94,8 @@ class LBMStepper:
                     f.ddf_gpu, self.boundary_values_gpu, np.uint64(self._step_count),
                     block=self.block, grid=self.boundary_grid, **launch_kw,
                 )
+            if self.moving is not None:
+                self._launch_moving(action_gpu, obs_gpu, advance, **launch_kw)
             self._launch_curved(action_gpu, obs_gpu, **launch_kw)
 
             if self._esopull:
@@ -106,6 +121,21 @@ class LBMStepper:
             # Any completed GPU step invalidates host DDF mirror.
             f.invalidate_host_ddf_cache()
 
+    def _launch_moving(self, action_gpu, obs_gpu, advance, **launch_kw):
+        m, tpb = self.moving, self.cfg.threads_per_block
+        if advance:  # M1 completes before M2 on the same stream
+            self.moving_cells_fn(
+                self.field.ddf_gpu, self.field.flag_gpu, m.band_gpu, m.geom_gpu,
+                np.uint32(m.n_bodies), np.uint32(m.n_cells), action_gpu, m.diag_gpu,
+                block=(tpb, 1, 1), grid=((m.n_cells + tpb - 1) // tpb, 1, 1), **launch_kw,
+            )
+        self.moving_fn(
+            self.field.ddf_gpu, self.field.flag_gpu, m.band_gpu, m.geom_gpu,
+            np.uint32(m.n_bodies), np.uint32(m.n_slots), action_gpu, obs_gpu,
+            m.diag_gpu, m.slot_debug_gpu if m.slot_debug_gpu is not None else np.uintp(0),
+            block=(tpb, 1, 1), grid=((m.n_slots + tpb - 1) // tpb, 1, 1), **launch_kw,
+        )
+
     def _launch_curved(self, action_gpu, obs_gpu, **launch_kw):
         f = self.field
         if f.n_curved == 0:

src/CelerisLab/simulation.py

@@ -26,6 +26,13 @@ import pycuda.driver as cuda
 from .config import LBMConfig, BodyConfig, load_lbm_config, load_body_config
 from .cuda.context import CudaContext, context_bound
 from .cuda import compiler
+from .lbm.descriptors import FLUID
+
+# Translating circles step only through their moving tables (fixed centre in S2a).
+PENDING_BLOCKED = "a pending translation target is advanced only by Simulation.run(1)"
+MOVING_BLOCKED = ("translating circles need the moving kernels and their tables (attached "
+                  "at initialize); "
+                  "positive-step runs are rejected")
 from .lbm.field import LBMField
 from .lbm.stepper import LBMStepper
 from .body.objects import SimObject
@@ -47,6 +54,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,13 +97,15 @@ 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
 
     # -- Object management ---------------------------------------------------
     def add_body(self, type: str = "circle", center=None,
-                 radius: float = 0.0, *, surface_actuation=None, **kwargs) -> int:
+                 radius: float = 0.0, *, surface_actuation=None, translation=None,
+                 **kwargs) -> int:
         """Add a simulation body. Returns body_id.
 
         Before ``initialize()`` the object is added to the formal registry
@@ -111,6 +121,9 @@ class Simulation:
                 angle, half_width (radians), normal_velocity (positive outward).
                 D2Q9/FP32/unshifted/DB only; set before initialization or stage
                 a new body. No runtime motion or amplitude setter.
+            translation: ``body.motion.Translation`` declaring a translating
+                circle (before initialize only; docs/moving_bodies.md). Until
+                the moving kernels exist, positive-step runs are rejected.
             **kwargs: passed to the geometry constructor.
 
         Returns:
@@ -121,6 +134,10 @@ class Simulation:
         type_lower = str(type).strip().lower()
         if surface_actuation is not None and type_lower != "circle":
             raise ValueError("Surface actuation is only supported on circles")
+        if translation is not None and (type_lower != "circle" or surface_actuation is not None
+                                        or self._initialized):
+            raise ValueError("translation needs a circle without surface actuation, "
+                             "added before initialize()")
 
         if type_lower == "circle":
             from .body.geometry.circle import CircleGeometry
@@ -129,6 +146,11 @@ class Simulation:
             obj = SimObject(obj_id=-1, geometry=geom,
                             center=center, radius=radius,
                             is_sensor=False, surface_actuation=surface_actuation)
+            if translation is not None:
+                from .body.motion import MovingCircle
+                obj.motion = MovingCircle((cx, cy), radius, translation,
+                                          on_failure=self._fail_moving)
+                self.stepper.blocked = MOVING_BLOCKED
         elif type_lower == "sensor":
             from .body.geometry.circle import CircleGeometry
             cx, cy = float(center[0]), float(center[1])
@@ -172,15 +194,99 @@ 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 _fail_moving(self, exc: BaseException) -> None:
+        # Device geometry is unknown after a failed moving step: same lifecycle as Jets.
+        self._failed = f"moving circle step failed: {exc!r}"
+        self.stepper.blocked = self._failed  # direct stepper calls are refused as well
+
+    def _attach_moving(self) -> None:
+        # Direct stepper calls stay blocked until the moving tables exist.
+        tables = getattr(self, "_moving_tables", None)
+        self.stepper.moving = tables
+        moving = self._moving_bodies()
+        if moving and tables is None:
+            self.stepper.blocked = MOVING_BLOCKED
+        elif any(obj.motion.pending is not None for obj in moving):
+            self.stepper.blocked = PENDING_BLOCKED  # direct steps never consume a target
+        else:
+            self.stepper.blocked = None
+
+    def _build_moving_tables(self):
+        moving = self._moving_bodies()
+        if not moving:
+            return None
+        from .body.motion import MovingTables
+        return MovingTables.upload([(obj.motion, obj.obj_id) for obj in moving])
+
+    @context_bound
+    def moving_diagnostics(self) -> np.ndarray:
+        """Per translating circle (in body order), columns ``MovingTables.DIAG``,
+        accumulated since obs were last zeroed (diagnostics only, not a force RHS)."""
+        self._require_usable()
+        tables = self.stepper.moving
+        if tables is None:
+            return np.zeros((0, 6), np.float32)
+        self.stream.synchronize()
+        return tables.download_diagnostics()
+
+    def _moving_bodies(self):
+        return [obj for obj in self.bodies.objects if obj.motion is not None]
+
+    def _reject_moving_bodies(self, operation: str) -> None:
+        if self._moving_bodies():
+            raise ValueError(f"{operation} is not supported with translating circles")
+
+    def _validate_moving_bodies(self) -> None:
+        """Reject unsupported modes and swept-band placement before any mutation."""
+        moving = self._moving_bodies()
+        if not moving:
+            return
+        cfg = self.lbm_cfg
+        if (not cfg.is_d2q9 or cfg.streaming != "double_buffer" or cfg.ddf_shifting
+                or cfg.store_precision != "FP32"):
+            raise ValueError("translating circles need D2Q9/FP32/unshifted/double_buffer")
+        if self.bodies.jets:
+            raise ValueError("Jets are not supported together with translating circles")
+        nx, ny = self.field.nx, self.field.ny
+        base = self.field.build_channel_flags().reshape(ny, nx)
+        for obj in moving:
+            x0, x1, y0, y1 = obj.motion.band()
+            if x0 < 1 or y0 < 1 or x1 > nx - 2 or y1 > ny - 2:
+                raise ValueError(f"swept band {(x0, x1, y0, y1)} touches the domain edge")
+            if not np.all((base[y0:y1 + 1, x0:x1 + 1] & FLUID) != 0):
+                raise ValueError(f"swept band {(x0, x1, y0, y1)} contains non-fluid base cells")
+            for other in self.bodies.objects:
+                if other is obj:
+                    continue
+                (cx, cy), r = other.center[:2], float(other.radius)
+                if other.motion is not None:
+                    ox0, ox1, oy0, oy1 = other.motion.band()
+                else:
+                    ox0, ox1, oy0, oy1 = cx - r - 2, cx + r + 2, cy - r - 2, cy + r + 2
+                if ox0 <= x1 and x0 <= ox1 and oy0 <= y1 and y0 <= oy1:
+                    raise ValueError(f"swept band {(x0, x1, y0, y1)} overlaps body {other.obj_id}")
+
+    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.
 
@@ -246,7 +352,7 @@ class Simulation:
 
     # -- Runtime control -----------------------------------------------------
     def set_body(self, id: int, omega: float = None,
-                 vx: float = None, vy: float = None) -> None:
+                 vx: float = None, vy: float = None, offset=None) -> None:
         """Set runtime body state (host action array only, no H2D).
 
         The GPU action buffer is updated at the next ``run()`` call.
@@ -256,8 +362,21 @@ class Simulation:
             omega: rotation speed (lattice units).  ``None`` = no change.
             vx: translation velocity x (future).  ``None`` = no change.
             vy: translation velocity y (future).  ``None`` = no change.
+            offset: target centre offset of a translating circle; with vx/vy
+                (default 0) it replaces the pending target after validation
+                against the applied offset (docs/moving_bodies.md, 4a). Only
+                ``run(1)`` advances it; direct stepper calls are blocked meanwhile.
         """
         obj = self.bodies.get(id)
+        motion = obj.motion
+        if motion is not None:
+            if offset is None and (vx is not None or vy is not None):
+                raise ValueError("a translating circle takes vx/vy only with an offset")
+            target = None if offset is None else motion.check_target(
+                offset, (0.0 if vx is None else vx, 0.0 if vy is None else vy))
+            vx = vy = None  # the motion target, not ObjectState, owns translation
+        elif offset is not None:
+            raise ValueError(f"body {id} was not declared with a translation")
         proposed = dict(vars(obj.state))
         proposed.update({k: float(v) for k, v in (("omega", omega), ("vx", vx), ("vy", vy))
                          if v is not None})
@@ -273,6 +392,9 @@ class Simulation:
         if vy is not None:
             obj.state.vy = float(vy)
             changed = True
+        if motion is not None and target is not None:
+            motion.pending = target
+            self._attach_moving()
         if changed:
             self.bodies._refresh_action_from_objects()
 
@@ -338,6 +460,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 +488,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 +508,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:
@@ -426,6 +551,7 @@ class Simulation:
             self.field, self._module, self.lbm_cfg,
         )
         self.stepper._step_count = _prev_step_count
+        self._attach_moving()
 
         self._assert_object_count_contract(expected_count=expected_n_objects)
 
@@ -446,10 +572,13 @@ 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.
         next_flags = self.bodies.build_flags(self.field.build_channel_flags())
+        self._validate_moving_bodies()
+        self._attach_moving()
         if self.bodies.jets and self.bodies.count == 0:
             raise ValueError("Jets need bodies with cut links")
         self.bodies.build_compact_lists_for(self.bodies.objects, next_flags)
@@ -466,6 +595,8 @@ class Simulation:
             self.bodies.sync_to_gpu(self.field, rebuild_flags=False)
         else:
             self.bodies.ensure_packed_buffers()
+        self._moving_tables = self._build_moving_tables()
+        self._attach_moving()
         self._assert_runtime_contracts()
         self._initialized = True
 
@@ -478,8 +609,10 @@ 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()")
+        self._reject_moving_bodies("sync_bodies()")
 
         self._assert_allocation_contract()
         # 1. Synchronize GPU
@@ -573,12 +706,20 @@ 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)
         checkpoint_interval = operator.index(checkpoint_interval)
         if steps < 0 or checkpoint_interval < 0:
             raise ValueError("steps and checkpoint_interval must be nonnegative")
+        moving = self._moving_bodies()
+        if steps > 0 and moving and self.stepper.moving is None:
+            raise RuntimeError(MOVING_BLOCKED)
+        if moving and checkpoint_interval > 0:
+            raise ValueError("checkpoints are not supported with translating circles")
+        if steps > 1 and any(obj.motion.pending is not None for obj in moving):
+            raise ValueError("a pending translation target needs run(1)")
         self._assert_runtime_contracts()
         # Discard any uncommitted body edits before stepping.
         if self.bodies.has_pending_edit():
@@ -595,6 +736,8 @@ class Simulation:
         # Zero obs segments before step group
         if zero_obs:
             self.bodies.zero_obs_async(stream)
+            if self.stepper.moving is not None:
+                self.stepper.moving.zero_diagnostics(stream)
 
         if checkpoint_interval > 0:
             done = 0
@@ -611,6 +754,9 @@ class Simulation:
                 if done < steps or done % checkpoint_interval == 0:
                     stream.synchronize()
                     self.save_checkpoint()
+        elif self.stepper.moving is not None and steps > 0:
+            self._run_moving(steps, stream, sync_obs)
+            return
         else:
             self.stepper.step(
                 steps,
@@ -629,6 +775,34 @@ class Simulation:
 
         stream.synchronize()
 
+    def _run_moving(self, steps, stream, sync_obs):
+        """Steps with translating circles: hold, or one advancing step that commits
+        each pending target only after synchronization and a zero ``no_donor`` count
+        (docs/moving_bodies.md 4a, 5). Any failure leaves the Simulation unusable."""
+        tables = self.stepper.moving
+        circles = [(obj.motion, obj.obj_id) for obj in self._moving_bodies()]
+        try:
+            advancing = tables.upload_records(circles, stream)
+            self.stepper._step(steps, action_gpu=self.bodies.action_gpu,
+                               obs_gpu=self.bodies.obs_gpu, stream=stream, advance=advancing)
+            self.bodies.increment_obs_steps(steps)
+            if sync_obs:
+                self.bodies.download_obs_full_async(stream)
+            stream.synchronize()
+            if advancing:
+                no_donor = tables.download_diagnostics()[:, tables.DIAG.index("no_donor")]
+                if np.any(no_donor > 0):
+                    raise RuntimeError("translating circle refill found no donor cell")
+                for motion, _ in circles:
+                    motion.commit()
+                tables.upload_records(circles, stream)  # first hold: old = new = applied, u = 0
+                stream.synchronize()
+                self._attach_moving()
+        except BaseException as exc:
+            for motion, _ in circles:  # device geometry/state unknown after any failure
+                motion.fail(exc)
+            raise
+
     def step(self, n: int = 1):
         """Advance *n* steps (convenience for interactive use)."""
         self.run(n)
@@ -637,15 +811,25 @@ 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()
 
+    @context_bound
     def get_flags(self) -> np.ndarray:
-        return self.field.flag.copy()
+        """Flag array; read back from the device when translating circles may have moved."""
+        if getattr(self.stepper, "moving", None) is None:
+            return self.field.flag.copy()
+        self._require_usable()
+        self.stream.synchronize()
+        flags = np.empty_like(self.field.flag)
+        cuda.memcpy_dtoh(flags, self.field.flag_gpu)
+        return flags
 
     # -- Runtime parameter updates -------------------------------------------
     @context_bound
@@ -654,24 +838,29 @@ 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
+        self._reject_moving_bodies("save_checkpoint()")
         if self.bodies.jets:
             raise NotImplementedError("Jets are not checkpointed yet")
         return save_checkpoint(self.field, self.stepper, self.lbm_cfg,
@@ -680,7 +869,9 @@ 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
+        self._reject_moving_bodies("load_checkpoint()")
         if not self._initialized:
             self.initialize()
         load_checkpoint(path, self.field, self.stepper, self.lbm_cfg,
@@ -699,6 +890,7 @@ class Simulation:
                 cuda.Context.synchronize()
                 # Release owned GPU resources while their context is current.
                 self.stepper = None
+                self._moving_tables = None
                 self.bodies = None
                 self.field = None
                 self._module = None

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

@@ -0,0 +1,536 @@
+"""Translating-circle host contract, S1 (docs/moving_bodies.md 4a, 5, 7).
+
+Host doubles only: no kernels exist yet, so every positive-step run with a
+moving body must be rejected before any device write.
+"""
+from types import SimpleNamespace
+
+import pytest
+
+from CelerisLab.body.motion import (MAX_STEP, R_MIN, MovingCircle, Translation,
+                                    enumerate_donor_coverage)
+
+BAD = (float("nan"), float("inf"))
+
+
+def circle(box=(-2., 2., -3., 3.), speed=.1):
+    return MovingCircle((50., 40.), 10., Translation(box, speed))
+
+
+def test_translation_and_radius_validation():
+    for box in ((0., 1., -1., 1.)[:3], (.5, 1., -1., 1.), (-1., 1., BAD[0], 1.)):
+        with pytest.raises(ValueError):
+            Translation(box, .1)
+    for speed in (0., MAX_STEP + 1e-9) + BAD:
+        with pytest.raises(ValueError):
+            Translation((0., 0., -1., 1.), speed)
+    with pytest.raises(ValueError):
+        MovingCircle((0., 0.), R_MIN - .1, Translation((0., 0., -1., 1.), .1))
+
+
+def test_targets_validate_against_applied_and_last_wins():
+    c = circle()
+    c.set_target((0., .3), (0., .05))
+    c.set_target((0., .5), (0., .05))  # second setter replaces, checked vs applied (0, 0)
+    assert c.pending == ((0., .5), (0., .05))
+    for offset, velocity in (((0., .51), (0., 0.)),           # Euclidean bound
+                             ((.36, .36), (0., 0.)),          # |dc|_2 = .509
+                             ((0., 3.1), (0., 0.)),           # outside box
+                             ((0., .1), (0., .11)),           # speed bound
+                             ((0., 0.), (0., .01)),           # hold with sliding velocity
+                             ((BAD[0], 0.), (0., 0.)), ((0., 0.), (0., BAD[1]))):
+        with pytest.raises(ValueError):
+            c.set_target(offset, velocity)
+        assert c.pending == ((0., .5), (0., .05))  # invalid replacement keeps the old target
+    c.set_target((.35, .35), (0., 0.))  # |dc|_2 = .495 inside the bound
+    assert c.applied == (0., 0.)
+
+
+def test_record_commit_and_hold():
+    c = circle()
+    assert c.record() == ((50., 40.), (50., 40.), (0., 0.))  # hold: no translation velocity
+    c.set_target((0., .25), (0., .05))
+    assert c.record() == ((50., 40.), (50., 40.25), (0., .05))
+    uploads = []
+    c.advance(uploads.append, lambda: None, lambda: None)
+    assert uploads == [((50., 40.), (50., 40.25), (0., .05))]
+    assert c.applied == (0., .25) and c.pending is None  # offset-only commit
+    c.advance(uploads.append, lambda: None, lambda: None)  # no new target: hold step
+    assert uploads[-1] == ((50., 40.25), (50., 40.25), (0., 0.)) and c.applied == (0., .25)
+    c.set_target((0., 0.), (0., -.05))  # return to zero offset is a real target
+    assert c.pending == ((0., 0.), (0., -.05))
+    c.advance(uploads.append, lambda: None, lambda: None)
+    assert c.applied == (0., 0.) and c.pending is None
+
+
+@pytest.mark.parametrize("stage", ("upload", "step", "synchronize"))
+def test_failure_after_upload_attempt_is_fail_closed(stage):
+    c = circle()
+    c.set_target((0., .25), (0., .05))
+
+    def call(name):
+        def run(*_):
+            if name == stage:
+                raise RuntimeError(f"{name} failed")
+        return run
+    with pytest.raises(RuntimeError, match=f"{stage} failed"):
+        c.advance(call("upload"), call("step"), call("synchronize"))
+    assert c.failed and c.applied == (0., 0.)
+    for later in (c.record, lambda: c.set_target((0., .1), (0., 0.)),
+                  lambda: c.advance(*(lambda *_: None,) * 3)):
+        with pytest.raises(RuntimeError, match="unknown"):
+            later()
+
+
+def test_band_covers_offset_box_with_margin():
+    assert circle().band() == (50 - 2 - 10 - 2, 50 + 2 + 10 + 2, 40 - 3 - 10 - 2, 40 + 3 + 10 + 2)
+
+
+@pytest.mark.parametrize("radius", (R_MIN, 4.3, 7.5, 12.0))
+def test_refill_donor_sets_nonempty_within_bounds(radius):
+    # Evidence for R >= R_MIN and |dc|_2 <= MAX_STEP on a fixed sample, not a proof.
+    assert enumerate_donor_coverage(radius) == 0
+
+
+pytest.importorskip("pycuda.driver")
+from CelerisLab import simulation as module  # noqa: E402
+from tests.unit.test_simulation_lifecycle import Owner, make_simulation  # noqa: E402
+
+
+class ClosableOwner(Owner):
+    """Owner double that refuses activation after close, like CudaContext.activate."""
+    def activate(self):
+        if self.is_closed:
+            raise RuntimeError("CUDA context is closed")
+        return super().activate()
+
+
+def moving_simulation(monkeypatch, *, moving=True):
+    sim = make_simulation(monkeypatch)
+    obj = SimpleNamespace(motion=circle() if moving else None, obj_id=0,
+                          center=(50., 40.), radius=10.)
+    writes = []
+    sim._failed, sim._initialized = None, True
+    sim.bodies = SimpleNamespace(objects=[obj], count=1, jets=[], action_gpu=None, obs_gpu=None,
+                                 has_pending_edit=lambda: False,
+                                 _upload_action_async=lambda s: writes.append("action"),
+                                 zero_obs_async=lambda s: writes.append("obs"),
+                                 increment_obs_steps=lambda n: writes.append("count"),
+                                 download_obs_full_async=lambda s: writes.append("download"))
+    sim.stepper = SimpleNamespace(step=lambda n, **k: writes.append(("step", n)), blocked=None,
+                                  moving=None)
+    sim._assert_runtime_contracts = lambda: None
+    return sim, writes
+
+
+STREAM = SimpleNamespace(synchronize=lambda: None)
+
+
+def test_positive_steps_rejected_before_device_writes(monkeypatch):
+    sim, writes = moving_simulation(monkeypatch)
+    sim.bodies.objects[0].motion.set_target((0., .1), (0., .01))
+    sim.bodies.objects[0].motion.applied = (0., 0.)
+    for call in (lambda: sim.run(1, stream=STREAM), lambda: sim.run(5, stream=STREAM),
+                 lambda: sim.step(), lambda: sim.run(3, stream=STREAM, checkpoint_interval=1)):
+        with pytest.raises(RuntimeError, match="moving kernels"):
+            call()
+    with pytest.raises(RuntimeError, match="moving kernels"):  # nominal at rest as well
+        sim.bodies.objects[0].motion.pending = None
+        sim.run(1)
+    sim.bodies.objects[0].state = SimpleNamespace(omega=.02)  # spin does not lift the guard
+    with pytest.raises(RuntimeError, match="moving kernels"):
+        sim.run(2, stream=STREAM)
+    assert writes == [] and sim.bodies.objects[0].motion.applied == (0., 0.)
+
+
+def test_run_zero_keeps_pending_and_static_runs_proceed(monkeypatch):
+    sim, writes = moving_simulation(monkeypatch)
+    motion = sim.bodies.objects[0].motion
+    motion.set_target((0., .1), (0., .01))
+    sim.run(0, stream=STREAM)  # existing action/obs work only; no motion upload or commit
+    assert motion.pending == ((0., .1), (0., .01)) and motion.applied == (0., 0.)
+    assert ("step", 0) in writes
+    static, writes = moving_simulation(monkeypatch, moving=False)
+    static.save_checkpoint = lambda *a: writes.append("checkpoint")
+    static.run(3, stream=STREAM)
+    static.run(4, stream=STREAM, checkpoint_interval=2)
+    assert writes.count(("step", 3)) == 1 and writes.count(("step", 2)) == 2
+
+
+def test_direct_stepper_entry_is_guarded():
+    from CelerisLab.lbm.stepper import LBMStepper
+    launches = []
+    module = SimpleNamespace(get_function=lambda name: lambda *a, **k: launches.append(name))
+    cfg = SimpleNamespace(streaming="double_buffer", threads_per_block=128, is_d2q9=True,
+                          nx=64, ny=32, nz=1)
+    field = SimpleNamespace(ddf_gpu="ddf", temp_gpu="temp", n_curved=0)
+    stepper = LBMStepper(field, module, cfg)
+    assert stepper.blocked is None
+    stepper.blocked = "translating circles need the moving kernels"
+    for n in (1, 3):
+        with pytest.raises(RuntimeError, match="moving kernels"):
+            stepper.step(n, action_gpu=None, obs_gpu=None)
+    assert launches == [] and (field.ddf_gpu, field.temp_gpu) == ("ddf", "temp")
+    assert stepper.step_count == 0
+    stepper.step(0, action_gpu=None, obs_gpu=None)  # zero steps stay a no-op
+    assert launches == [] and stepper.step_count == 0
+
+
+def test_moving_body_rejects_topology_and_checkpoint(monkeypatch):
+    sim, _ = moving_simulation(monkeypatch)
+    for call in (sim.sync_bodies, sim.save_checkpoint, lambda: sim.load_checkpoint("x.h5")):
+        with pytest.raises(ValueError, match="translating circles"):
+            call()
+
+
+# -- Real-object routing, initialize rejections, packing, lifecycle -------------
+import numpy as np  # noqa: E402
+
+from CelerisLab.body.geometry.circle import CircleGeometry  # noqa: E402
+from CelerisLab.body.manager import ObjectManager  # noqa: E402
+from CelerisLab.body.objects import SimObject  # noqa: E402
+from CelerisLab.lbm.descriptors import FLUID, SOLID  # noqa: E402
+
+
+class Cfg:
+    dim, is_d3q19, is_d2q9 = 2, False, True
+    streaming, ddf_shifting, store_precision, y_wall_bc = "double_buffer", False, "FP32", "bounce_back"
+
+
+def body(x, y, r, translation=None, **kw):
+    obj = SimObject(obj_id=-1, geometry=CircleGeometry(x, y, r), center=(x, y), radius=r, **kw)
+    if translation is not None:
+        obj.motion = MovingCircle((x, y), r, translation)
+    return obj
+
+
+def routed_simulation(monkeypatch):
+    sim = make_simulation(monkeypatch)
+    sim._failed, sim._initialized, sim.lbm_cfg = None, False, Cfg()
+    sim.stepper = SimpleNamespace(blocked=None)
+    sim.bodies = ObjectManager(nx=120, ny=80, nz=1, nq=9, cfg=Cfg())
+    return sim
+
+
+def test_add_body_and_set_body_route_targets(monkeypatch):
+    sim = routed_simulation(monkeypatch)
+    bid = sim.add_body("circle", (50., 40.), 10., translation=Translation((-2., 2., -3., 3.), .1))
+    still = sim.add_body("circle", (100., 40.), 5.)
+    motion = sim.bodies.get(bid).motion
+    sim._moving_tables = "tables"
+    sim._attach_moving()
+    assert sim.stepper.blocked is None
+    sim.set_body(bid, offset=(0., .3), vy=.05)
+    assert sim.stepper.blocked == module.PENDING_BLOCKED  # direct steps never consume a target
+    sim.set_body(bid, offset=(0., .4), vy=.05)  # last wins
+    for kwargs in (dict(offset=(0., .6)), dict(vy=.05), dict(vx=.05),
+                   dict(offset=(0., .1), vy=.2), dict(offset=(0., 0.), vy=.01)):
+        with pytest.raises(ValueError):
+            sim.set_body(bid, **kwargs)
+        assert motion.pending == ((0., .4), (0., .05))  # invalid input keeps the target
+    sim.set_body(bid, omega=.01)  # omega-only keeps the pending target
+    obj = sim.bodies.get(bid)
+    assert obj.state.omega == .01 and (obj.state.vx, obj.state.vy) == (0., 0.)
+    assert motion.pending == ((0., .4), (0., .05)) and motion.applied == (0., 0.)
+    with pytest.raises(ValueError, match="not declared"):
+        sim.set_body(still, offset=(0., .1))
+    sim.set_body(still, vx=.01)  # stationary API unchanged
+    assert sim.bodies.get(still).state.vx == .01 and sim.bodies.get(still).motion is None
+    for kwargs in (dict(type="sensor"), dict(surface_actuation=dict(angle=0., half_width=.2,
+                                                                    normal_velocity=.01))):
+        kind = kwargs.pop("type", "circle")
+        with pytest.raises(ValueError, match="translation"):
+            sim.add_body(kind, (20., 20.), 5., translation=Translation((0., 0., 0., 0.), .1),
+                         **kwargs)
+
+
+class Tables:
+    """MovingTables double: records uploads, steps and diagnostics reads."""
+    DIAG = ("no_donor",)
+
+    def __init__(self, log, no_donor=0., fail=None):
+        self.log, self.no_donor, self.fail = log, no_donor, fail
+
+    def upload_records(self, circles, stream):
+        records = [m.record() for m, _ in circles]
+        self.log.append(("upload", records))
+        uploads = sum(1 for w in self.log if isinstance(w, tuple) and w[0] == "upload")
+        if self.fail == "upload" or (self.fail == "hold_upload" and uploads == 2):
+            raise RuntimeError("upload failed")
+        return any(r[0] != r[1] for r in records)
+
+    def zero_diagnostics(self, stream):
+        self.log.append("zero_diag")
+
+    def download_diagnostics(self):
+        return np.array([[self.no_donor]], np.float32)
+
+
+def advancing_simulation(monkeypatch, **tables):
+    sim, writes = moving_simulation(monkeypatch)
+    sim.bodies.objects[0].motion = circle()
+    sim.stepper = SimpleNamespace(
+        moving=Tables(writes, **tables), blocked=None,
+        _step=lambda n, advance=False, **k: writes.append(("step", n, advance)))
+    from CelerisLab.lbm.stepper import LBMStepper
+    sim.stepper.step = LBMStepper.step.__get__(sim.stepper)  # real public entry
+    sim._moving_tables = sim.stepper.moving
+    return sim, writes, sim.bodies.objects[0].motion
+
+
+def test_run_one_advances_and_commits_after_sync(monkeypatch):
+    sim, writes, motion = advancing_simulation(monkeypatch)
+    motion.set_target((0., .25), (0., .05))
+    with pytest.raises(ValueError, match="run\\(1\\)"):  # multi-step with a target
+        sim.run(2, stream=STREAM)
+    assert writes == [] and motion.pending is not None
+    sim.run(0, stream=STREAM)  # zero steps never consume the target
+    assert motion.pending is not None
+    sim.run(1, stream=STREAM)
+    uploads = [w for w in writes if isinstance(w, tuple) and w[0] == "upload"]
+    assert uploads[0][1] == [((50., 40.), (50., 40.25), (0., .05))]
+    assert uploads[1][1] == [((50., 40.25), (50., 40.25), (0., 0.))]  # first hold after move
+    assert ("step", 1, True) in writes
+    assert motion.applied == (0., .25) and motion.pending is None
+    writes.clear()
+    sim.run(3, stream=STREAM)  # hold steps: no advance
+    assert ("step", 3, False) in writes and motion.applied == (0., .25)
+
+
+@pytest.mark.parametrize("failure", ("upload", "first_sync", "no_donor", "hold_upload",
+                                     "hold_sync"))
+def test_failed_advance_fails_simulation(monkeypatch, failure):
+    kwargs = dict(no_donor=1.) if failure == "no_donor" else dict(fail=failure)
+    sim, writes, motion = advancing_simulation(monkeypatch, **kwargs)
+    motion._on_failure = sim._fail_moving
+    motion.set_target((0., .25), (0., .05))
+    syncs = []
+
+    def synchronize():  # stands in for an asynchronous failure surfacing at sync
+        syncs.append(1)
+        if (failure, len(syncs)) in (("first_sync", 1), ("hold_sync", 2)):
+            raise RuntimeError("sync failed")
+    with pytest.raises(RuntimeError):
+        sim.run(1, stream=SimpleNamespace(synchronize=synchronize))
+    # Before the step-group sync and no_donor check nothing is committed; the hold
+    # record is staged after the commit (applied already matches the moved device).
+    committed = failure in ("hold_upload", "hold_sync")
+    assert motion.applied == ((0., .25) if committed else (0., 0.))
+    assert (motion.pending is None) == committed
+    assert motion.failed and sim._failed is not None
+    assert sim.stepper.blocked == sim._failed  # direct stepper reuse refused too
+    with pytest.raises(RuntimeError, match="Simulation failed"):
+        sim.run(1, stream=STREAM)
+
+
+def test_failed_hold_step_blocks_the_direct_stepper(monkeypatch):
+    # No target is pending, so only the failure latch can refuse direct steps.
+    from CelerisLab.lbm.stepper import LBMStepper
+    sim, writes, motion = advancing_simulation(monkeypatch, fail="upload")
+    motion._on_failure = sim._fail_moving
+    launches = []
+    kernels = SimpleNamespace(get_function=lambda name: lambda *a, **k: launches.append(name))
+    cfg = SimpleNamespace(streaming="double_buffer", threads_per_block=128, is_d2q9=True,
+                          nx=64, ny=32, nz=1)
+    field = SimpleNamespace(ddf_gpu="ddf", temp_gpu="temp", flag_gpu="flag", n_curved=0)
+    stepper = LBMStepper(field, kernels, cfg)
+    stepper.moving = sim.stepper.moving
+    sim.stepper = stepper
+    assert motion.pending is None and stepper.blocked is None
+    with pytest.raises(RuntimeError, match="upload failed"):
+        sim.run(2, stream=STREAM)  # hold step
+    assert sim._failed is not None and stepper.blocked == sim._failed
+    with pytest.raises(RuntimeError, match="moving circle step failed"):
+        stepper.step(1, action_gpu=None, obs_gpu=None)
+    assert launches == [] and stepper.step_count == 0
+    assert (field.ddf_gpu, field.temp_gpu) == ("ddf", "temp")
+
+
+def test_moving_tables_layout_and_attachment(monkeypatch):
+    from CelerisLab.body.motion import MovingTables
+    a = MovingCircle((50., 40.), 10., Translation((-1., 1., -3., 3.), .1))
+    b = MovingCircle((90.25, 30.5), 5., Translation((0., 0., 0., 0.), .1))
+    band, geom = MovingTables.layout([(a, 2), (b, 5)])
+    wa, ha = 63 - 37 + 1, 55 - 25 + 1
+    assert band == [37, 25, wa, ha, 0, 2, 83, 23, 16, 16, wa * ha * 8, 5]
+    assert geom == [50., 40., 50., 40., 10., 0., 0., 90.25, 30.5, 90.25, 30.5, 5., 0., 0.]  # hold
+    sim = routed_simulation(monkeypatch)
+    sim.add_body("circle", (50., 40.), 10., translation=Translation((0., 0., 0., 0.), .1))
+    sim._moving_tables = "tables"
+    sim._attach_moving()
+    assert sim.stepper.moving == "tables" and sim.stepper.blocked is None
+    sim._moving_tables = None
+    sim._attach_moving()
+    assert sim.stepper.moving is None and sim.stepper.blocked is not None
+
+
+class CountingOwner(Owner):
+    """Owner double that counts nested activations like CudaContext."""
+    def __init__(self):
+        super().__init__()
+        self.active = 0
+
+    def activate(self):
+        from contextlib import contextmanager
+
+        @contextmanager
+        def scope():
+            self.active += 1
+            try:
+                yield self
+            finally:
+                self.active -= 1
+        return scope()
+
+
+def test_public_diagnostics_keep_owner_active_through_sync_and_download(monkeypatch):
+    sim = make_simulation(monkeypatch)
+    sim.ctx, sim._failed = CountingOwner(), None
+    seen = []
+    sim._stream = SimpleNamespace(synchronize=lambda: seen.append(("sync", sim.ctx.active)))
+    expected = np.arange(6, dtype=np.float32).reshape(1, 6)
+
+    def download():
+        seen.append(("download", sim.ctx.active))
+        return expected
+    sim.stepper = SimpleNamespace(moving=SimpleNamespace(download_diagnostics=download))
+    assert sim.moving_diagnostics() is expected
+    assert [name for name, _ in seen] == ["sync", "download"]
+    assert all(active > 0 for _, active in seen) and sim.ctx.active == 0
+
+
+def test_stepper_launches_moving_links_before_static_curved():
+    from CelerisLab.lbm.stepper import LBMStepper
+    launches = []
+    module = SimpleNamespace(get_function=lambda name: lambda *a, **k: launches.append((name, a, k)))
+    cfg = SimpleNamespace(streaming="double_buffer", threads_per_block=128, is_d2q9=True,
+                          nx=64, ny=32, nz=1)
+    curved = SimpleNamespace(**{f"{name}_gpu": name for name in (
+        "fluid_idx", "dir", "q", "rx", "ry", "rz", "fallback_class", "body_id")},
+        surface_velocity=None, surface_velocity_gpu=None, link_exchange_gpu=None)
+    field = SimpleNamespace(ddf_gpu="ddf", temp_gpu="temp", flag_gpu="flag", n_curved=5,
+                            curved=curved, invalidate_host_ddf_cache=lambda: None)
+    stepper = LBMStepper(field, module, cfg)
+    stepper._launch_force_region = stepper._launch_sensor = lambda *a, **k: None
+    stepper.moving = SimpleNamespace(band_gpu="band", geom_gpu="geom", diag_gpu="diag",
+                                     slot_debug_gpu=None, n_bodies=1, n_slots=300)
+    stepper.step(2, action_gpu="act", obs_gpu="obs")
+    names = [name for name, _, _ in launches]
+    assert names == ["MovingCircleLinks", "CurvedBoundaryKernel", "OneStep"] * 2
+    assert stepper.step_count == 2
+    assert launches[1][1][0] == "ddf" and launches[3][1][0] == "temp"  # swap between steps
+    args, kw = launches[0][1], launches[0][2]
+    assert args[:4] == ("ddf", "flag", "band", "geom") and int(args[5]) == 300
+    assert kw["grid"] == (3, 1, 1) and kw["block"] == (128, 1, 1)
+
+
+def test_recompile_installs_a_guarded_stepper(monkeypatch):
+    from CelerisLab.lbm.stepper import LBMStepper
+    sim = routed_simulation(monkeypatch)
+    sim.add_body("circle", (50., 40.), 10., translation=Translation((0., 0., -1., 1.), .1))
+    launches = []
+    kernel_module = SimpleNamespace(get_function=lambda name: lambda *a, **k: launches.append(name))
+    monkeypatch.setattr(module.compiler, "build_module", lambda *a, **k: kernel_module)
+    sim.lbm_cfg = SimpleNamespace(validate_runtime=lambda: None, viscosity=.01, omega=1.,
+                                  streaming="double_buffer", threads_per_block=128,
+                                  is_d2q9=True, nx=120, ny=80, nz=1)
+    sim._compiled_viscosity = .01
+    sim._assert_allocation_contract = lambda: None
+    sim._resolve_compile_arch = lambda: "sm_70"
+    sim._assert_object_count_contract = lambda **k: None
+    sim.field = SimpleNamespace(module=None, invalidate_module_constant_cache=lambda: None,
+                                _upload_params=lambda: None, _runtime_params={},
+                                ddf_gpu="ddf", temp_gpu="temp", n_curved=0)
+    sim.stepper = SimpleNamespace(_step_count=7, blocked=None)
+    sim._recompile()
+    assert isinstance(sim.stepper, LBMStepper) and sim.stepper.step_count == 7
+    with pytest.raises(RuntimeError, match="moving kernels"):
+        sim.stepper.step(1, action_gpu=None, obs_gpu=None)
+    assert launches == [] and (sim.field.ddf_gpu, sim.field.temp_gpu) == ("ddf", "temp")
+    assert sim.stepper.step_count == 7
+
+
+def validate(monkeypatch, objects, flags=None, **cfg):
+    sim = routed_simulation(monkeypatch)
+    for obj in objects:
+        sim.bodies.add(obj)
+    for name, value in cfg.items():
+        setattr(sim.lbm_cfg, name, value)
+    base = np.full(120 * 80, FLUID, np.uint16) if flags is None else flags
+    sim.field = SimpleNamespace(nx=120, ny=80, build_channel_flags=lambda: base.copy())
+    sim._validate_moving_bodies()
+
+
+def test_initialize_rejects_band_and_modes(monkeypatch):
+    t = Translation((-1., 1., -3., 3.), .1)
+    validate(monkeypatch, [body(50., 40., 10., t), body(100., 40., 5.)])  # accepted
+    walls = np.full(120 * 80, FLUID, np.uint16)
+    walls[60 + 30 * 120] = SOLID
+    cases = [([body(10., 40., 10., t)], None, {}, "domain edge"),
+             ([body(50., 40., 10., t)], walls, {}, "non-fluid"),
+             ([body(50., 40., 10., t), body(68., 40., 4.)], None, {}, "overlaps"),
+             ([body(50., 40., 10., t), body(50., 60., 3., is_sensor=True)], None, {}, "overlaps"),
+             ([body(50., 40., 10., t), body(64., 40., 3., is_force_region=True)], None, {}, "overlaps"),
+             ([body(50., 40., 10., t)], None, {"streaming": "esopull"}, "double_buffer"),
+             ([body(50., 40., 10., t)], None, {"ddf_shifting": True}, "unshifted"),
+             ([body(50., 40., 10., t)], None, {"is_d2q9": False}, "D2Q9")]
+    for objects, flags, cfg, match in cases:
+        with pytest.raises(ValueError, match=match):
+            validate(monkeypatch, objects, flags, **cfg)
+
+
+def test_initialize_validates_before_mutation(monkeypatch):
+    sim = routed_simulation(monkeypatch)
+    sim.bodies.add(body(10., 40., 10., Translation((0., 0., 0., 0.), .1)))
+    uploads = []
+    sim.field = SimpleNamespace(nx=120, ny=80, build_channel_flags=lambda: np.full(9600, FLUID, np.uint16),
+                                upload_flags=lambda: uploads.append("flags"))
+    sim._assert_allocation_contract = lambda: None
+    with pytest.raises(ValueError, match="domain edge"):
+        sim.initialize()
+    assert uploads == [] and not sim._initialized
+
+
+def test_moving_circle_excluded_from_static_links():
+    def manager(*objects):
+        mgr = ObjectManager(nx=120, ny=80, nz=1, nq=9, cfg=Cfg())
+        ids = [mgr.add(obj) for obj in objects]
+        return mgr, ids
+    reference_mgr, _ = manager(body(100., 40., 5.), body(20., 40., 6.))
+    mixed_mgr, ids = manager(body(100., 40., 5.), body(20., 40., 6.),
+                             body(55., 40., 10., Translation((0., 0., -1., 1.), .1)))
+    reference = reference_mgr.build_compact_lists_for(reference_mgr.objects)
+    mixed = mixed_mgr.build_compact_lists_for(mixed_mgr.objects)
+    for a, b in zip(reference[:8], mixed[:8]):  # static link columns byte-identical
+        assert np.asarray(a).tobytes() == np.asarray(b).tobytes()
+    assert len(reference[0]) > 0 and ids[2] not in set(np.asarray(mixed[3]).tolist())
+
+
+def test_failed_moving_step_marks_simulation_failed(monkeypatch):
+    sim = routed_simulation(monkeypatch)
+    sim.ctx = ClosableOwner()
+    bid = sim.add_body("circle", (50., 40.), 10., translation=Translation((0., 0., -1., 1.), .1))
+    motion = sim.bodies.get(bid).motion
+    motion.set_target((0., .2), (0., .05))
+
+    def fail(*_):
+        raise RuntimeError("sync failed")
+    with pytest.raises(RuntimeError, match="sync failed"):
+        motion.advance(lambda *_: None, lambda: None, fail)
+    assert sim._failed is not None and motion.applied == (0., 0.)
+    for call in (lambda: sim.run(0), lambda: sim.run(1), sim.step, sim.sync_bodies,
+                 sim.save_checkpoint, lambda: sim.load_checkpoint("x.h5")):
+        with pytest.raises(RuntimeError, match="Simulation failed"):
+            call()
+    with pytest.raises(RuntimeError, match="unknown"):  # no retry on the failed circle
+        motion.advance(lambda *_: None, lambda: None, lambda: None)
+    with sim.ctx.activate():  # a refused close keeps the failed state and resources
+        with pytest.raises(RuntimeError, match="active CUDA"):
+            sim.close()
+    assert sim._failed is not None and sim.bodies is not None and not sim._closed
+    sim.close()  # retry outside the activation succeeds
+    assert sim._closed and sim._failed is not None
+    with pytest.raises(RuntimeError, match="CUDA context is closed"):  # no re-entry after close
+        sim.run(1)
+    with pytest.raises(RuntimeError, match="unknown"):  # the failed circle stays refused
+        motion.record()

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

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

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
@@ -31,20 +32,26 @@ class RefinementMRTHost(unittest.TestCase):
                      "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(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)
-        # 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):
-            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()

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

@@ -294,7 +294,11 @@ def test_passive_overlap_failure_keeps_live_state(device, entry, second_x):
     sim = SimpleNamespace(bodies=bodies, field=field, lbm_cfg=field.cfg,
                           _initialized=entry == "runtime_edit",
                           _assert_allocation_contract=lambda: None,
-                          stream=SimpleNamespace(synchronize=device.synchronize))
+                          stream=SimpleNamespace(synchronize=device.synchronize),
+                          stepper=SimpleNamespace(blocked=None))
+    for name in ("_require_usable", "_moving_bodies", "_reject_moving_bodies",
+                 "_validate_moving_bodies", "_attach_moving"):
+        setattr(sim, name, getattr(Simulation, name).__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)
@@ -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

@@ -81,7 +81,7 @@ def run_blocked_patch(speed):
                                               level_objects=[(), (source,)],
                                               observe_raw_exchange=True) as backend:
             c, f = backend._levels
-            flags = backend.fine_boundary.flags.reshape(40, 40)
+            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)
@@ -101,7 +101,7 @@ def run_blocked_patch(speed):
             with patch.object(f, "kernel", fine):
                 backend.step(2)
             assert len(traces) == 4 and backend.step_count == 2
-            mass = backend.read_raw_exchange(0)["total"]["mass_into_fluid"]
+            mass = backend.read_raw_exchange(0, level=1)["total"]["mass_into_fluid"]
             return traces, roi, backend.download(), float(mass)
 
 

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

@@ -24,7 +24,8 @@ from CelerisLab.lbm import refinement_stepper
 from CelerisLab.lbm.blocks import BlockLayout2D
 from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
 from CelerisLab.lbm.refinement import ACTIVE, OUTER_GHOST, FixedPatch2D
-from CelerisLab.lbm.refinement_stepper import FixedPatchStepper, blocked_base_tables
+from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+from CelerisLab.lbm.refinement_tables import blocked_coarse_tables
 from test_refinement_blocked_fine import NX, NY, BOUNDS, config, state, circle, actuation
 
 BASE_BLOCKS = ((24, 32), (12, 8), (3, 4), (1, 1))
@@ -108,7 +109,7 @@ class BlockedBaseCPU(unittest.TestCase):
                 leaves[[0, -1]] = False
             for by, bx in BASE_BLOCKS:
                 layout = BlockLayout2D.build(NY, NX, by, bx, periodic=(True, mode == "periodic"))
-                donors, targets, kind = blocked_base_tables(p, layout, order, mode != "periodic")
+                donors, targets, kind = blocked_coarse_tables(p, layout, order, mode != "periodic")
                 ly, lx = layout.local_shape
                 L = ly * lx
                 expected, source_cells = [], []
@@ -138,7 +139,7 @@ class BlockedBaseCPU(unittest.TestCase):
         p = FixedPatch2D.from_config(config("periodic"), BOUNDS)
         order = np.argsort(p.coalesce_coarse)
         layout = BlockLayout2D.build(NY, NX, 3, 4)
-        donors, targets, _ = blocked_base_tables(p, layout, order)
+        donors, targets, _ = blocked_coarse_tables(p, layout, order)
         dense = np.arange(9 * NY * NX, dtype=np.float32).reshape(9, NY, NX)
         initial = layout.to_blocked(dense)
         arrivals = np.arange(len(targets), dtype=np.float32) + 10000
@@ -178,7 +179,7 @@ class BlockedBaseCPU(unittest.TestCase):
         # With one parent block, all eight grazing sources are owner cores.
         # Each is read for four fine ghost children: late explosion corrupts32.
         single = BlockLayout2D.build(NY, NX, NY, NX)
-        donors, targets, _ = blocked_base_tables(p, single, order)
+        donors, targets, _ = blocked_coarse_tables(p, single, order)
         initial = single.to_blocked(dense)
         staged = initial.copy()
         staged.ravel()[targets.astype(np.intp)] = arrivals
@@ -204,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)
@@ -255,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"):

tests/validation/multigrid/test_refinement_blocked_fine.py

@@ -21,7 +21,8 @@ from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
 from CelerisLab.lbm.descriptors import SOLID
 from CelerisLab.lbm.refinement import (ACTIVE, INNER_GHOST, OUTER_GHOST, FixedPatch2D,
                                        prepare_fine_boundary)
-from CelerisLab.lbm.refinement_stepper import FixedPatchStepper, blocked_fine_tables
+from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+from CelerisLab.lbm.refinement_tables import blocked_fine_tables
 
 NX, NY, BOUNDS, STEPS = 32, 24, (8, 6, 24, 18), 8
 FINE_BLOCKS = ((32, 40), (16, 20), (16, 10))  # fine window 32x40: 1x1, 2x2, 4x2
@@ -125,7 +126,7 @@ class BlockedFineCPU(unittest.TestCase):
                 patch.object(cuda, "mem_alloc", side_effect=alloc), patch.object(cuda, "memcpy_htod"):
             backend = FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS],
                                         level_blocks=[(12, 8), FINE_BLOCKS[0]])
-            self.assertEqual(len(backend.fine_layout.halo_dst), 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()
@@ -178,7 +179,7 @@ class BlockedFineCPU(unittest.TestCase):
             backend = FixedPatchStepper(config("periodic"), context=context, refinements=[BOUNDS],
                                         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):
@@ -186,8 +187,8 @@ class BlockedFineCPU(unittest.TestCase):
             backend, calls = self.mocked_backend([circle(spec)])
             c, f = backend._levels
             coarse, cur, spare = c.current, f.current, f.spare
-            flags, kind = backend._level_boundaries[1].flags, backend._levels[1].kind
-            velocity = backend._level_boundaries[1].device.get("surface_velocity")
+            flags, kind = backend._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()
@@ -212,8 +213,8 @@ class BlockedFineCPU(unittest.TestCase):
     def test_actuation_without_velocity_buffer_fails_before_bfl(self):
         for speed in (0.0, .002):  # any spec, even zero amplitude
             backend, calls = self.mocked_backend([circle(actuation(speed))])
-            self.assertIsNotNone(backend._fine_tables["surface_velocity"])
-            del backend._level_boundaries[1].device["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]

tests/validation/multigrid/test_refinement_followup.py

@@ -304,7 +304,7 @@ 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):
@@ -383,8 +383,8 @@ class NormalBuildGPU(unittest.TestCase):
         expected = scalar_interval(*[a.astype(np.float64) for a in initial], bounds)
         with CudaContext(0) as context, FixedPatchStepper(cfg, context=context, refinements=[bounds],
                                                              level_blocks=[(13, 17), (14, 16)]) as backend:
-            self.assertEqual(backend.patch.fine_shape, (14, 16))
-            p = backend.patch
+            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)

tests/validation/multigrid/test_refinement_gpu.py

@@ -91,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.
@@ -101,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))
@@ -165,13 +165,13 @@ class CoupledGPU(unittest.TestCase):
                     with self.assertRaisesRegex(RuntimeError, "failed"):
                         operation()
                 raw = level_storage(fine)
-                self.assertTrue(np.isfinite(raw.ravel()[backend.patch.coalesce_fine]).all())
+                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 eight_steps(self, backend, spec, maps):
         """Rest and seeded nonequilibrium, 8 steps against the scalar reference."""
-        p = backend.patch
+        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))

tests/validation/multigrid/test_refinement_jets.py

@@ -62,8 +62,8 @@ class JetStepperCPU(unittest.TestCase):
         C = (16.0, 12.0)
         whole = Jet(Sector(C, .1, math.pi, 3.), Radial(C, inward=False), amplitude=0.)
         backend, _ = mocked([circle()], [whole], FINE_BLOCKS[1])
-        _, jet_id = backend._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):
@@ -76,10 +76,10 @@ class JetStepperCPU(unittest.TestCase):
         source = SimObject(99, geometry=CircleGeometry(*C, 2.15), center=C, radius=2.15)
         whole = Jet(Sector(C, .1, math.pi, 3.), Radial(C, inward=False), amplitude=0.)
         backend, _ = mocked([source], [whole], FINE_BLOCKS[1])
-        _, jet_id = backend._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")
@@ -97,8 +97,8 @@ class JetStepperCPU(unittest.TestCase):
     def test_blocked_jet_velocity_matches_legacy_rows(self):
         legacy, _ = mocked([circle(actuation(SPEED))], (), FINE_BLOCKS[1])
         driven, _ = mocked([circle()], [jet()], FINE_BLOCKS[1])
-        a = legacy._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)
@@ -106,20 +106,20 @@ 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._level_boundaries[1].device["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):
@@ -128,24 +128,24 @@ class JetStepperCPU(unittest.TestCase):
             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)
-            before = 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)
-            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):
@@ -170,7 +170,7 @@ def run(mode, fine_block, *, objects, jets, switch=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)
 

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

@@ -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)