docs: S0 draft — wall impulse in momentum ledger, no_donor is a hard failure

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

docs/development/src_coherence.md

@@ -1,8 +1,14 @@
 # src coherence stage (from main 914906e)
 
 Owner direction 2026-10-28: unified style over historical compatibility, inside `src`.
-No numerical change is intended; kernel text was moved verbatim (PTX/bit parity is
-pending a GPU run).
+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`)
 
@@ -28,20 +34,25 @@ pending a GPU run).
 
 ## 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) | `_device_boundaries[l]` (list from L0, `None` at L0 unless single-level static circles) |
-| `_level_boundaries` (private, b4eacf7) | `_device_boundaries` |
-| `_fine_tables`, `_actuated`, `_jet_cache` | `_device_boundaries[1].tables` / `.actuated` / `.jet_cache` |
+| `_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) | `_boundaries` (private list behind `boundaries`) |
+| `_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) |
@@ -56,13 +67,35 @@ the replaced level relative to its parent).
 
 ## Responsibilities
 
-- Construction: `_check_levels` (inputs), `_kernel_configs` (per-level kernel configs),
-  `Hierarchy.edge_frame` (geometry frames), `_build_level` / `_build_pair` / `_upload_boundary`
-  (device). Observation reads share `_observed_boundary` + `_summarize`.
+- 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`. MRT shares one path; SRT keeps its two explicit
-  arithmetic variants (level: unclamped, no-force, shifted direct-momentum feq; channel: uniform
-  OneStep arithmetic with omega clamp), selected by `BLOCK_ROOT_BC`.
+  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
@@ -72,6 +105,12 @@ the replaced level relative to its parent).
   `FixedPatchStepper._kernel_configs` (reset only on level 1 of a slot_reservoir root;
   generates each level's config.h). Folding them needs a config.h/PTX diff gate.
 - No single-`Simulation` facade over FixedPatchStepper: the two own different state models
-  (uniform node frame, body sync, forcing, checkpoint vs. immutable blocked hierarchy, edge
+  (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/moving_bodies.md

@@ -0,0 +1,339 @@
+# Moving circles on the uniform dense grid — S0 DRAFT contract
+
+Status: **DRAFT for review, not frozen.** No implementation exists. 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 of 1 cell or more 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 or applied),
+u = pending.velocity or 0}`. 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` |
+| 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`, `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.
+
+## 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.

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

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

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

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

src/CelerisLab/lbm/refinement_stepper.py

@@ -317,26 +317,25 @@ class FixedPatchStepper:
         self._window_valid = False
         self._window_start = None
         self._exchange_valid = False
-        self._exchange_spacing = tuple(self.hierarchy.edge_frame(l)[1] for l in range(depth))
-        self._exchange_centers = tuple(tuple(tuple(o.geometry.center) for o in objects)
-                                       for objects in level_objects)
         root = None
         if depth == 1 and level_objects[0]:
             tables = static_circle_tables(self._layouts[0], level_objects[0])
             root = _BlockedBoundary(None, tables, (), False)
         # The one boundary/Jet record per level; host geometry and Jets are views of it.
-        self._device_boundaries = [root] + [
+        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 = _fluid_leaves(self.hierarchy, self._device_boundaries)
+        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 remap_finest.
-        self._root_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"
@@ -391,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._device_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)
@@ -414,7 +413,7 @@ 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._device_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)
@@ -576,12 +575,12 @@ 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 any(b is not None for b in self._device_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._device_boundaries[2]
+            actual = self._boundaries[2]
             jets = () if actual is None else actual.jets
             old = RemapTopology.from_config(cfg, pairs[0], pairs[1].bounds,
                                            child_objects=objects, child_jets=jets,
@@ -603,7 +602,7 @@ class FixedPatchStepper:
                     or not np.array_equal(old.surface_velocity, new.surface_velocity)):
                 raise ValueError("remap requires the identical pinned physical boundary and BFL support")
             # Check sparse block alignment/launch size before any device read/allocation.
-            layout = new.child.fine_layout(self._level_blocks[2])
+            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("level 2 blocks exceed launch limits")
             with self.ctx.activate():
@@ -621,7 +620,7 @@ class FixedPatchStepper:
                 phase="postcollision_g", **options)
             candidate = type(self)(
                 cfg, context=self.ctx, refinements=[pairs[0].bounds, new.child.bounds],
-                level_blocks=list(self._level_blocks),
+                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)
@@ -651,7 +650,7 @@ class FixedPatchStepper:
             if (isinstance(level, bool) or not isinstance(level, Integral)
                     or not 1 <= level < self.hierarchy.depth):
                 raise ValueError("Jet level must name an existing refined level (1 .. depth-1)")
-            boundary = self._device_boundaries[level]
+            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
@@ -711,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)
@@ -755,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._device_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)
@@ -778,7 +777,7 @@ class FixedPatchStepper:
                 raise RuntimeError("window exchange observation is disabled")
             self._window_valid = False
             try:
-                for b in self._device_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:
@@ -791,7 +790,7 @@ class FixedPatchStepper:
         if (isinstance(level, bool) or not isinstance(level, Integral)
                 or not 0 <= level < len(self._levels)):
             raise ValueError("level must name an existing level")
-        b = self._device_boundaries[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")
@@ -799,9 +798,9 @@ class FixedPatchStepper:
 
     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=self._exchange_spacing[level],
+                                  body_id=body_id, level=level, h=b.spacing,
                                   time_end=self._step_count, window=window,
-                                  reference_center=self._exchange_centers[level][body_id], **kw)
+                                  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."""
@@ -888,7 +887,7 @@ class FixedPatchStepper:
                         self._stream.synchronize()
                 except Exception as exc:
                     error = exc
-                for b in getattr(self, "_device_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
@@ -923,11 +922,11 @@ class FixedPatchStepper:
     @property
     def boundaries(self):
         """Host boundary geometry per level (None: no bodies/Jets on a refined level)."""
-        return tuple(None if b is None else b.host for b in self._device_boundaries)
+        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._device_boundaries)
+        return tuple(() if b is None else b.jets for b in self._boundaries)
 
     @property
     def leaf(self):

src/CelerisLab/lbm/refinement_tables.py

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

tests/unit/test_jet_update_failure.py

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

tests/unit/test_multiple_pinned_remap.py

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

tests/unit/test_pinned_body_remap.py

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

tests/unit/test_raw_exchange.py

@@ -85,7 +85,7 @@ class RawExchange(unittest.TestCase):
 
     def test_four_level_schedule_publication_and_fixed_storage(self):
         run, host, calls, arrays = backend()
-        b = run._device_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._device_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_refinement_mrt.py

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

tests/unit/test_runtime_child_remap.py

@@ -70,11 +70,11 @@ class RuntimeChildRemap(unittest.TestCase):
         run._status = 'ready'
         for bad in ([], [(40,16,48,40)], [(16,16,32,24),(16,27,32,55)]):
             with self.assertRaises(ValueError): run.remap_finest(bad)
-        records, pairs = run._device_boundaries, run.hierarchy.pairs
+        records, pairs = run._boundaries, run.hierarchy.pairs
         for attribute, value in (('_wall_aware', True),
                                  # level-2 record whose host geometry differs from objects=()
-                                 ('_device_boundaries', [*records[:2], SimpleNamespace(host=object(), jets=())]),
-                                 ('_device_boundaries', [records[0], object(), records[2]]),  # level-1 body
+                                 ('_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_finest(B)

tests/unit/test_shifted_jets.py

@@ -67,7 +67,7 @@ class ShiftedJets(unittest.TestCase):
     def test_switch_upload_only_then_raw_remap_keeps_live_velocity(self):
         run,module,allocations,storage,objects=ready();clone=None
         try:
-            boundary=run._device_boundaries[2]
+            boundary=run._boundaries[2]
             before=run._download_raw()
             original=boundary.tables['surface_velocity'].copy()
             module.cuda.memcpy_htod.reset_mock()
@@ -79,7 +79,7 @@ class ShiftedJets(unittest.TestCase):
             for a,b in zip(before.values,run._download_raw().values):self.assertEqual(a.tobytes(),b.tobytes())
             clone=run.remap_finest(B,objects=objects)
             self.assertEqual(clone.jets[2][0].amplitude,-jet().amplitude)
-            np.testing.assert_array_equal(clone._device_boundaries[2].tables['surface_velocity'],
+            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()
@@ -92,12 +92,12 @@ class ShiftedJets(unittest.TestCase):
     def test_invalid_amplitude_does_not_upload_or_commit(self):
         run,module,_,_,_=ready()
         try:
-            before=run._device_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._device_boundaries[2].tables['surface_velocity'])
+            np.testing.assert_array_equal(before,run._boundaries[2].tables['surface_velocity'])
         finally:run.close()
 
 

tests/unit/test_single_level.py

@@ -73,7 +73,7 @@ class SingleLevel(unittest.TestCase):
 
     def test_circle_solids_are_not_state(self):
         run, _, _ = backend((circle(),))
-        solid = run._device_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_union_body_support.py

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

tests/validation/multigrid/test_block_stepper.py

@@ -240,7 +240,7 @@ class SingleLevelGPU(unittest.TestCase):
         with CudaContext(0) as context:
             for split in SPLITS:
                 with single(context=context, block=split, objects=[circle()]) as backend:
-                    tables = backend._device_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_exchange_window.py

@@ -59,7 +59,7 @@ class ExchangeWindowCPU(unittest.TestCase):
             stack.enter_context(patch.object(cuda, "memcpy_dtoh",
                                              side_effect=lambda out, src: out.fill(0)))
             backend = build(m)
-            b = backend._device_boundaries[1]
+            b = backend._boundaries[1]
             self.assertIsNotNone(b.window_integral)
             with self.assertRaises(RuntimeError):
                 backend.start_exchange_window()  # not initialized
@@ -96,7 +96,7 @@ class ExchangeWindowCPU(unittest.TestCase):
         with mocked_backend() as m, ExitStack() as stack:
             stack.enter_context(patch.object(cuda, "memset_d8_async"))
             backend = build(m, window=False)
-            self.assertIsNone(backend._device_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._device_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)
@@ -212,7 +212,7 @@ class ExchangeWindowShiftedMRTGPU(unittest.TestCase):
                     run.step(2)
                     if run is not plain:
                         run.start_exchange_window()
-                b = windowed._device_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)
@@ -228,7 +228,7 @@ class ExchangeWindowShiftedMRTGPU(unittest.TestCase):
                     same_bits(self, a, e)
                 rows = [np.empty((count, 4), np.float64) for _ in range(2)]
                 cuda.memcpy_dtoh(rows[0], b.window_integral)
-                cuda.memcpy_dtoh(rows[1], split._device_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, level=1)

tests/validation/multigrid/test_multiple_pinned_remap_gpu.py

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

tests/validation/multigrid/test_pinned_body_remap_gpu.py

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

tests/validation/multigrid/test_raw_exchange_gpu.py

@@ -105,7 +105,7 @@ class RawExchangeGPU(unittest.TestCase):
                     arrays = [state(l.shape) for l in measured._levels]
                     measured.initialize(*arrays)
                     plain.initialize(*arrays)
-                    b = measured._device_boundaries[1]
+                    b = measured._boundaries[1]
                     count = len(b.tables['cell'])
                     kernel = measured._levels[1].kernel
                     samples = []
@@ -158,7 +158,7 @@ class RawExchangeGPU(unittest.TestCase):
                  FixedPatchStepper(cfg, context=ctx, **settings) as plain:
                 arrays = [state(l.shape, ux=0.) for l in run._levels]
                 run.initialize(*arrays); plain.initialize(*arrays)
-                b = run._device_boundaries[3]
+                b = run._boundaries[3]
                 count = len(b.tables['cell'])
                 kernel = run._levels[3].kernel
                 samples, independent_samples = [], []
@@ -213,7 +213,7 @@ class RawExchangeGPU(unittest.TestCase):
                  FixedPatchStepper(CFG, context=ctx, **fine((16, 20)), observe_raw_exchange=True) as amr:
                 amr.initialize(*(state(l.shape) for l in amr._levels))
                 uniform.initialize(state((9, 48, 64)))
-                b = amr._device_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._device_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']))

tests/validation/multigrid/test_refinement_blocked_fine.py

@@ -187,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._device_boundaries[1].flags, backend._levels[1].kind
-            velocity = backend._device_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()
@@ -213,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._device_boundaries[1].tables["surface_velocity"])
-            del backend._device_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_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._device_boundaries[1].jet_cache
-        self.assertGreater(len(np.unique(backend._device_boundaries[1].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._device_boundaries[1].jet_cache
+        _, 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._device_boundaries[1].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._device_boundaries[1].tables["surface_velocity"]
-        b = driven._device_boundaries[1].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,8 +106,8 @@ class JetStepperCPU(unittest.TestCase):
 
     def test_blocked_amplitude_update_rewrites_only_velocity_column(self):
         backend, htod = mocked([circle()], [jet()], FINE_BLOCKS[1])
-        before = backend._device_boundaries[1].tables["surface_velocity"].copy()
-        device = backend._device_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, level=1)
@@ -115,7 +115,7 @@ class JetStepperCPU(unittest.TestCase):
             self.assertEqual(htod.call_count, 1)
             self.assertIs(ptr, device)
             np.testing.assert_allclose(sent, -before, rtol=1e-6, atol=0)
-            np.testing.assert_array_equal(backend._device_boundaries[1].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, level=1)
             self.assertEqual(htod.call_count, 1)
@@ -133,7 +133,7 @@ class JetStepperCPU(unittest.TestCase):
             with self.assertRaises(IndexError):
                 backend.set_jet_amplitude(1, -SPEED, level=1)
             self.assertEqual(backend._status, status)
-            before = backend._device_boundaries[1].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)):
@@ -141,7 +141,7 @@ class JetStepperCPU(unittest.TestCase):
                     backend.set_jet_amplitude(0, -SPEED, level=1)
             self.assertEqual(backend._status, "failed", fine_block)
             self.assertEqual(backend.jets[1][0].amplitude, SPEED)
-            np.testing.assert_array_equal(backend._device_boundaries[1].tables["surface_velocity"], before)
+            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"):

tests/validation/multigrid/test_shifted_jets_gpu.py

@@ -102,7 +102,7 @@ class Ledger:
 
     def step(self):
         run=self.run;self.samples.clear();run.step()
-        b=run._device_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])
@@ -140,7 +140,7 @@ class ShiftedJetsGPU(unittest.TestCase):
                         level_objects=[(),(obj,)],level_jets=[(),(whole,)],observe_raw_exchange=True) as run:
                     run._initialize_raw(initial(run,rho,u=.013,v=0.))
                     expected=run._download_raw();ledger=Ledger(self,run,1)
-                    self.assertTrue(run._device_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))
@@ -201,7 +201,7 @@ class ShiftedJetsGPU(unittest.TestCase):
             self.exact_raw(before,old._download_raw())
             self.exact_raw(_StoredPopulations(expected.g,True),new._download_raw())
             self.assertEqual(new.jets[2][0].amplitude,-.001)
-            np.testing.assert_array_equal(new._device_boundaries[2].tables['surface_velocity'],old._device_boundaries[2].tables['surface_velocity'])
+            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')