Reserve a pressure condition on Jet, velocity remains the only implemented one

Jet gains a last field condition='velocity'. 'pressure' raises
NotImplementedError in __post_init__ before profile/amplitude validation, so a
pressure jet never reaches surface_link_velocity, tables or kernels; unknown
values raise ValueError. Amplitude stays a lattice velocity (not a flow rate);
no pressure units, field reuse, tables or kernel changes are introduced.

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

README.md

@@ -123,7 +123,7 @@ with CudaContext(0) as ctx, FixedPatchStepper(cfg, (8, 6, 24, 18), context=ctx,
 ```
 
 Scope and limits:
-- D2Q9, SRT, FP32, double buffer, unshifted, no LES or forcing; one fixed 2:1 patch.
+- D2Q9, SRT, FP32, double buffer, unshifted, no LES or forcing; fixed 2:1 patch(es), two levels (or one nested third).
 - `coarse_boundary="periodic"` (default) or `"slot_reservoir"` (bottom slot, top
   reservoir; patch needs `y0 >= 4`, `y1 <= ny-4`). Periodic mode ignores the wall/slot config fields.
 - Object centers/radii use the coarse **cell-edge** frame (cell `I` spans `[I, I+1)`),
@@ -135,10 +135,71 @@ Scope and limits:
   divide the WHOLE fine allocation, i.e. `2*(y1-y0)+8` by `2*(x1-x0)+8` including the
   four-cell ghost band. It gives results identical to dense storage and is
   infrastructure for future block AMR, not a speedup.
+- Several separated patches: pass a sequence of `(x0, y0, x1, y1)`
+  rectangles as `bounds`, at least 4 coarse cells apart, with both `fine_block`
+  and `base_block`. 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
+  patch windows (active region plus 4-cell ghost band, origin `(2*y0-4, 2*x0-4)`)
+  are stored; each window must align with `fine_block`. Order of rectangles does
+  not matter. `fine_objects`/`jets` work as for one patch: each circle (same
+  coarse cell-edge frame) must lie wholly inside one patch's active region with
+  its boundary stencil clear of transfer cells; a jet may not select links of two
+  different bodies, even across patches. No wall-touching patches or nesting yet.
+
+  ```python
+  cfg = LBMConfig(nx=32, ny=16, viscosity=.1)
+  rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None]
+  with CudaContext(0) as ctx, FixedPatchStepper(
+          cfg, ((4, 4, 12, 12), (16, 4, 24, 12)), context=ctx,
+          fine_block=(4, 4), base_block=(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.step(32)
+      coarse_g, fine_g = run.download()  # coarse: NaN on every covered cell;
+                                         # fine: NaN on every non-ACTIVE cell
+  ```
+- Three levels (one nested chain): add `child=(x0, y0, x1, y1)` in
+  **L1 window cells** (the patch's fine grid including its 4-cell ghost band) and
+  `child_block`, with `fine_block` and `base_block`. L2 must sit inside L1's
+  ACTIVE cells, away from L1's ghost band and transfer footprint (checked before
+  any GPU allocation). `initialize` and `download` take/return three arrays; each
+  level is NaN on its non-ACTIVE (and child-covered) cells. One coarse step = 2 L1
+  steps = 4 L2 steps. One L1 patch only. Optional stationary circles and Jets
+  live on L2 via `child_objects` and `child_jets`; `fine_objects`/`jets` retain
+  their L1 meaning and must be empty in this mode. All object and Jet geometry
+  uses the **L0 edge frame**, not the child bounds' L1 window frame. L2 has
+  dx=dt=0.25, so prescribed lattice velocity is unchanged. Each whole body and
+  its BFL stencil must be inside L2 ACTIVE cells, outside the transfer footprint.
+  A Jet may select one physical body across multiple blocks, never two bodies.
+
+  ```python
+  from CelerisLab.body.geometry.circle import CircleGeometry
+  from CelerisLab.body.objects import SimObject
+  from CelerisLab.body.jets import Jet, Tube, Parallel
+
+  cfg = LBMConfig(nx=16, ny=16, viscosity=.1)
+  body = SimObject(0, geometry=CircleGeometry(8., 8., .8), center=(8., 8.), radius=.8)
+  jet = Jet(Tube((6.5, 8.), 0., 3., .3), Parallel(0.), amplitude=.002)
+  rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None]
+  with CudaContext(0) as ctx, FixedPatchStepper(
+          cfg, (4, 4, 12, 12), context=ctx, fine_block=(4, 4), base_block=(4, 4),
+          child=(8, 8, 16, 16), child_block=(4, 4),
+          child_objects=[body], child_jets=[jet]) as run:
+      shapes = (run.patch.coarse_shape, run.patch.fine_shape, run.child_patch.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.step(4)
+      coarse_g, l1_g, l2_g = run.download()  # solids/covered/nonleaves are NaN
+  ```
 - Not available: per-link surface exchange/force readout in blocked mode,
-  checkpoints, moving bodies, multiple patches or dynamic refinement, `Simulation`
+  checkpoints, moving bodies, L1 bodies/Jets in a three-level chain,
+  several patches per level with nesting, dynamic refinement, `Simulation`
   integration. Validation so far is storage/operator equivalence and short runs,
-  not accuracy, convergence or long-time stability.
+  not accuracy, convergence or long-time stability. The new L2-body GPU gates
+  (layouts, masks, amplitude reversal, legacy equivalence) and independent
+  every-substep numerical replay remain required acceptance checks; CPU tests
+  alone do not establish GPU numerical correctness.
 
 ## Installation
 

docs/development/raw_exchange.md

@@ -0,0 +1,73 @@
+# Raw multilevel boundary exchange (opt-in candidate)
+
+Use `FixedPatchStepper(..., observe_raw_exchange=True)` and, after a successful
+base step, `read_raw_exchange(body_id, level=l, window="last_base_step")`.
+The alternative window is `"last_substep"`. Body IDs are object positions within
+the named level, not the object's original user label. The torque reference is
+that circle's geometry.center (the center used to construct BFL arms), reported
+in the root edge frame. A differing SimObject.center does not change this reference. Current geometry
+restrictions are unchanged: fixed complete bodies on the finest level, isolated
+from the transfer footprint. New observation currently requires blocked body
+storage; dense-body configurations raise before device allocation. Existing
+dense diagnostics are unchanged.
+
+Each result contains `wall`, `aperture`, `total`, with `force`, `torque`,
+`impulse`, `angular_impulse`, `mass_into_fluid`, and `mass_rate`. These are raw
+boundary momentum reaction, not validated stress or whole-nozzle mechanical
+load. No GI/advective correction or guessed thrust is added. Apertures retain
+their geometric classification when their prescribed velocity is zero.
+
+All results are base-lattice units per unit out-of-plane depth. For local
+spacing and timestep h=dt, raw `(jx,jy,tz,dm)` converts to impulse
+`(h*dt*jx,h*dt*jy,h*h*dt*tz)` and mass `h*h*dm`.
+The accumulator intentionally accepts h and dt separately to make the spatial
+and temporal conversion explicit; the current acoustic hierarchy supplies dt=h.
+The stored rx/ry arms run from geometry.center to the actual wall hit, not to
+the fluid node. For this raw exchange j=c_i(a+b), moving the arm from the fluid
+node to the hit adds q*c_i, whose cross product with j is zero. Thus either arm
+gives the same raw torque (up to floating-point rounding); this equivalence
+must not be assumed for a future nonparallel stress/advective correction.
+
+Force/torque are impulse/angular impulse divided by the reported window length.
+`dm=b-a` is positive into the exterior fluid, while j is the opposite of the
+fluid's momentum change across that cut link. Interpolated BFL can introduce
+mass drift; dm is an observation, not an assertion of impermeability.
+
+`last_base_step` means [t-1,t] even after step(n); `last_substep` means [t-h,t]
+and is not a common instantaneous sample across levels. A Jet amplitude set
+before step(1) applies to all its substeps. Reads are repeatable; step(0) leaves
+them unchanged. Initialization and a remap candidate invalidate observation.
+Only successful stream synchronization publishes a window; failed instances
+reject reads rather than returning partially advanced measurements.
+
+Storage is one FP32 [N,4] raw array and one FP64 [N,4] integral per measured
+level, plus row metadata. No substep history is stored. Same-stream sequencing
+is halo/restriction -> BFL -> independent row accumulation -> step. There are
+no accumulation atomics or body obs offsets (blocked n_objects remains zero).
+Disabled observation allocates neither observation buffer and launches no
+accumulation kernel. The standalone periodic block caller passes explicit
+null optional arguments.
+
+## Gates
+
+CPU (standard library unittest, plus existing NumPy dependency):
+
+    PYTHONPATH=src python -m unittest discover -s tests/unit -p test_raw_exchange.py -v
+
+GPU, explicitly on a CUDA host with the existing environment:
+
+    PYTHONPATH=src python -m unittest discover -s tests/validation/multigrid -p test_raw_exchange_gpu.py -v
+
+The GPU suite includes an independent actual-population cut-link census and
+ACTIVE-fluid momentum/mass ledger with non-body boundary flux, observation
+on/off population equivalence, dense/blocked row comparison, three block
+partitions, zero/step/reversed/disabled Jet amplitudes, depth-four eight-sample
+accumulation, and a short-time whole-domain uniform-fine comparison with the
+same physical geometry, viscosity and Jet velocity. The uniform reference
+uses a test-only adapter to the standalone kernel's optional inputs.
+
+These are implementation checks. They do not establish physical Jet/stress
+accuracy, long-time AMR error, or convergence order. Resolution/long-time
+studies and independently validated continuum load benchmarks remain separate.
+GPU gates were prepared but not run in the authoring environment (no CUDA,
+PyCUDA or nvcc); run them before merging this candidate.

docs/fixed_three_level_multipatch.md

@@ -0,0 +1,76 @@
+# Fixed three-level separated patches
+
+`FixedPatchStepper` supports exactly three static, blocked levels with separated
+rectangle sets on either adjacent pair. This is not general AMR: there is no
+per-region timestep, moving geometry, dynamic refinement or L1 body support.
+All regions share one current/spare buffer pair per level. The existing schedule
+advances L0 once, L1 twice and L2 four times per base step, with aggregate transfers.
+No CUDA kernel changes are involved.
+
+## Explicit coordinate frames
+
+Keep the existing scalar `child=(x0,y0,x1,y1)` for a single scalar parent patch.
+Its coordinates retain their original meaning: half-open integer indices in the
+allocated L1 window, including its ghost band.
+
+The new keyword `child_rectangles_in_parent=[...]` is mutually exclusive with
+`child`. It accepts a scalar or a set parent and interprets every rectangle in
+the actual allocated L1 array, including ghost bands and holes. It never infers
+region-local coordinates. The physical root cell edge corresponding to an L1
+index `(x,y)` is `parent.fine_origin + (x,y)/2`:
+
+- A scalar L1 parent has a shifted window origin.
+- A set L1 parent has origin `(0,0)` and public shape `(2*ny,2*nx)`.
+- A child set has public shape twice the allocated L1 shape. Its root origin is
+  the L1 allocated origin, even when the parent is a scalar window.
+
+For a 16x16 root, scalar parent `(4,4,12,12)` and scalar child `(8,8,16,16)`
+are physically equivalent to parent set `[(4,4,12,12)]` with child set
+`[(12,12,20,20)]`. A child set on the *scalar* parent instead uses
+`[(8,8,16,16)]`. Reusing the same numbers after changing the parent frame can
+move the physical region.
+
+Both adjacent pairs require a Chebyshev gap of at least four in their own
+parent-cell units. The whole child union must preserve inherited ACTIVE
+coverage, uncovered ACTIVE donors/receivers, surviving pulls, and the upstream
+restriction footprint. Sparse holes remain unowned; inherited INNER remains
+stream-only and OUTER remains invalid. Block sizes must align with every
+support edge, including its four-cell ghost band.
+
+## Example
+
+```python
+run = FixedPatchStepper(
+    cfg, [(4, 4, 20, 28), (24, 4, 44, 28)], context=ctx,
+    base_block=(4, 4), fine_block=(4, 4), child_block=(4, 4),
+    child_rectangles_in_parent=[(16, 16, 32, 40), (56, 16, 80, 40)],
+)
+```
+
+This example uses a 48x32 root. `initialize` and `download` use three Q-major
+dense arrays. Public arrays and masks still scale with the full allocated
+extent; sparse device buffers do not remove that host-memory limit.
+
+Stationary whole-circle `child_objects` and `child_jets` use root edge
+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.
+
+## Verification
+
+CPU geometry, table, pre-CUDA rejection, schedule and cleanup tests:
+
+```sh
+PYTHONPATH=src python -m unittest discover -s tests/unit -p 'test_refinement*.py' -v
+```
+
+Explicit GPU gates are in
+`tests/validation/multigrid/test_refinement_three_level_multipatch_gpu.py`.
+They cover physical-frame equivalence, twins, asymmetric refinement, multiple
+children, reversed order, independently varied aligned layouts, public/raw
+masks and two L2 bodies with static/zero/driven/switched Jets. Snapshots are
+0/1/2/8. CPU mocks do not execute CUDA and do not establish numerical accuracy.
+GPU equivalence gates supplement the separately frozen independent dense
+three-level numerical replay; passing either alone is not a conservation,
+performance or physical-accuracy claim.

docs/synchronized_leaf_remap.md

@@ -0,0 +1,160 @@
+# Synchronized leaf G remap and manual continuation
+
+This is a CPU state-transfer component, **not a runtime adaptive solver**.
+It changes only a separated L2 rectangle set in a three-level, body-free,
+periodic hierarchy. L0/L1 topology, relaxation parameters, direction order,
+physical frame and completed time remain unchanged. Existing geometry checks
+validate child support/nesting and separation. The CPU helper itself adds no CUDA, scheduler, allocator, boundary handling,
+collision or ghost refill. A bounded runtime adapter is described below.
+
+```python
+from CelerisLab.lbm.refinement import FixedPatchSet2D
+from CelerisLab.lbm.synchronized_remap import RemapTopology, remap_synchronized
+
+# cfg: D2Q9/SRT/FP32/unshifted/double_buffer, nz=1, no LES
+parent = FixedPatchSet2D.from_config(cfg, [(4, 4, 20, 28), (24, 4, 44, 28)])
+old = RemapTopology.from_config(cfg, parent, [(16, 16, 32, 40)])
+new = RemapTopology.from_config(cfg, parent, [(20, 16, 36, 40)])
+candidate = remap_synchronized(
+    old, new, current_g, base_count=7, level_steps=(7, 14, 28),
+    phase="postcollision_g", origin=(0., 0.), spacing=1.,
+)
+```
+
+The example has a 48×32 root; both child supports allow 4×4 blocks.
+Child bounds use the parent's allocated L1 coordinates, not rectangle-local
+indices. A scalar FixedPatch2D parent is also supported: its shifted L1 window
+origin is retained. Physical lower edges are `origin + spacing*(x,y)` on L0,
+`origin + spacing*(parent.fine_origin + (x,y)/2)` on L1 and
+`origin + spacing*(parent.fine_origin + (x,y)/4)` on L2.
+
+Factory options `coarse_boundary` and `objects` must declare `"periodic"` and
+an empty body collection (all levels). 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
+actual state has no bodies/forcing, that all device work and old readers are
+quiescent, and that inputs are the **current owned-leaf postcollision G**.
+Count/phase stamps check consistency, not freshness or active-buffer identity.
+
+`current_g` is three NumPy Q-major arrays `(9,ny,nx)`. For blocked inputs,
+pass `layouts=(base_layout, fine_layout, child_layout)`; each optional `None`
+entry denotes dense input. Every leaf must have an owner. Only owner cores
+are read, never redundant block halos. Mixed precision is rejected. FP64 is
+accepted for CPU reference tests; this does not enable FP64 runtime storage.
+Outputs `candidate.g` are separate dense arrays regardless of input layout.
+Memory scales with public grid extents, not only sparse allocated blocks.
+
+Same-level leaf overlap is copied bitwise in physical coordinates. Newly
+refined L2 copies all nine values from its valid old L1 parent. Newly exposed
+L1 receives the complete four-child mean, including q0. Covered parents,
+ghosts and holes never serve as donors; NaN/Inf outside leaves is ignored,
+while nonfinite required leaves or nonpositive density reject the operation.
+Restriction uses FP64 quarter-weight accumulation and rounds to the input
+dtype. New nonleaves are NaN. Source arrays/metadata are not changed.
+
+`candidate` retains `base_count`, `level_steps`, `phase`, `origin`, `spacing`
+and new `topology.leaf_masks`. Existing `FixedPatchStepper.initialize` resets
+the counter: using it is **not** a continuation adapter. This component neither
+publishes a candidate nor retires old state. Runtime table rebuild/halo refill
+and correct first-step scheduling are handled only by the explicit adapter below.
+
+Tests accumulate leaf-area integrals in FP64 for each population, mass and
+momentum. The declared absolute tolerance is `2e-13 * sum(abs(old population
+integrals))` for FP64, `2e-7 * scale` for FP32, accounting for reduction and
+restriction rounding (not a solver accuracy bound). Coarse→fine→coarse copy
+round-trips preserve state; coarsening arbitrary fine structure loses it and
+is not invertible. No kinetic accuracy or stability claim follows.
+
+```sh
+PYTHONPATH=src python -m unittest discover -s tests/unit -p 'test_synchronized_remap.py' -v
+```
+
+## Explicit runtime continuation
+
+`replacement = current.remap_child(rectangles)` returns a new ready
+`FixedPatchStepper` at the **same absolute step_count**. It supports only
+body-free periodic three-level runs whose existing L2 uses
+`child_rectangles_in_parent`. Scalar or union L1 parents are supported; scalar
+`child=` is deliberately rejected. L0/L1, relaxation and block sizes stay fixed.
+No tagger, scheduler, in-place replacement, or general rollback is provided.
+
+The method holds the stepper lock, validates geometry, synchronizes its stream,
+and downloads the actual current owned leaves. It certifies postcollision G
+and completed 1:2:4 counts internally, applies the CPU helper, rebuilds all
+geometry/transfer/kind/halo buffers, and uploads both current and spare buffers.
+No old halo is reused. The regular step schedule fills transit ghosts/halos
+before use. `initialize` itself still resets the count to zero.
+
+The old stepper remains ready, unchanged and independently usable. The caller
+owns **both** objects and must close both before closing their shared CUDA
+context. Candidate construction or upload failure releases only candidate
+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)])
+# Explicitly choose the returned continuation only after success.
+current.close()
+current = replacement
+current.step()  # advances the preserved absolute count by one
+```
+
+CPU lifecycle mocks cover fresh current-buffer reads, exact transfer, separate
+ownership, time preservation, ordinary initialization reset, unsupported states
+and geometry, and candidate/source failure cleanup. GPU tests are predeclared
+in `tests/validation/multigrid/test_runtime_child_remap_gpu.py`: identity next
+step bitwise, shifted A→B→A with two intervals between changes, uniform moving
+and nonequilibrium states, exact CPU remap and overlap, both buffers reset,
+finite leaves, exact masks and absolute counts. They were not run on the CPU
+host. The separate independent FP64 capture/replay gate is max population
+error **3e-6**, with exact masks; remap mass/momentum allow FP32 roundoff.
+Per-population conservation is not an advancement target under collision.
+
+### Bounded pinned-body extension
+
+`replacement = current.remap_child(rectangles, child_objects=original_objects)`
+also accepts stationary Circles in separated L2 rectangles. Every body’s
+**complete containing rectangle must remain identical**. Several bodies may
+share one pinned rectangle; there need not be any fluid-only rectangles.
+Separated fluid-only rectangles may be added, removed or moved. Existing
+nesting, four-cell gap, full-circle containment, reverse cut-link 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.
+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
+IDs, q, lever arms, hits, BFL support,
+legacy SurfaceActuation specification and current per-link velocity against
+the existing boundary. Omission, changed geometry/actuation and changed body
+rectangle fail without device reads or allocations; reordered bodies fail
+too. Root object IDs are still assigned by list position, as in the constructor.
+
+Jets are taken from the live L2 boundary, including the most recent successful
+`set_jet_amplitude(..., level=2)`. There is no replacement-Jets argument or
+constructor-time Jet snapshot. New blocked addresses may differ; the physical
+boundary and velocities must match exactly.
+
+`RemapTopology.from_config(..., child_objects=..., child_jets=...)` prepares a
+typed `FineBoundary2D`; it derives L2 fluid leaves as ACTIVE minus SOLID. No
+caller mask is accepted. Solids, ghosts and holes are never donors and remain
+NaN, physical-fluid overlap is exact, and exposed L1 still restricts all four
+old fluid children. The helper checks the old/new body pins and boundary too.
+Its immutable `pinned_rectangles` is a canonical, deduplicated tuple, empty for body-free layouts; rectangle input
+order does not affect it. Scalar-parent window offsets remain part of the
+coordinate conversion.
+
+Focused CPU mocks: `tests/unit/test_pinned_body_remap.py`. GPU declarations:
+`tests/validation/multigrid/test_pinned_body_remap_gpu.py`, covering no-op
+next-step bitwise identity, moved-fluid A→B→A, current-amplitude preservation,
+physical boundary/velocity equality, source independence and absolute counts.
+Multiple-body CPU coverage: `tests/unit/test_multiple_pinned_remap.py`.
+`tests/validation/multigrid/test_multiple_pinned_remap_gpu.py` adds two pins,
+two independently updated Jets, variable fluid-region counts and add/move/remove
+immediate comparisons with CPU transfer. All nine solid-excluded area-weighted
+population integrals, exact overlap and nonleaf NaNs are checked.
+These new GPU cases have **not been run here**. They do not establish evolution
+accuracy. A 3e-6 FP32-vs-FP64 evolution gate requires a separate independent
+BFL reference; the body-free replay gate above does not establish that claim.

research/shear_refinement/README.md

@@ -0,0 +1,70 @@
+# One measured-strain refinement example
+
+Research only: one completed base-time observation, one possible growth-only
+fluid L2 enlargement, then two more base steps. L0/L1 and the Circle-containing
+L2 rectangle stay fixed. No scheduler, coarsening, packing, CUDA changes, or
+claim of accuracy, stable long-time adaptation, or optimal mesh selection.
+
+The periodic initial shear wave is evolved for two base steps. The downloaded
+postcollision G moments supply actual L1 velocity, at dx=.5 and the parent's
+physical origin. The tagger measures |S|=sqrt(2 S:S), not vorticity. We use
+s*=|S|/0.05, cutoff .001, no wall band/history/buffer. This deliberately small
+cutoff exercises the mechanism; it is not an error estimator or calibrated
+physical criterion. A high cutoff yields no change with the same observation.
+
+The eligible corridor is L1 indices [72,76) x [20,36), immediately right of the
+fluid child. Only L1 fluid leaves are valid gradient donors and owned voters;
+covered cells, holes and ghosts cannot supply gradients. L1 has no physical
+solid cells: the Circle is entirely in the fixed L2 patch. The cropped sparse
+arrays are never treated as a full periodic domain. One-sided derivatives may
+be used at missing-neighbor edges. Unknown gradients, unavailable corridor
+cells, and tags outside the corridor are reported separately. Ignored tags do
+not certify adequate resolution elsewhere. Old coverage is always retained.
+
+`grow_rectangle` knows only a rectangle, boolean requests, alignment and area
+cap. It returns an unchanged/enlarged tuple or raises ValueError before any
+solver mutation. The caller validates the entire proposed mesh with the
+existing RemapTopology and sparse-layout guards; remap_child repeats its full
+runtime checks. This example targets exactly the original mesh once, not a
+reusable adaptation loop. A rejection propagates explicitly; no partial
+proposal is installed. All solver instances are independently closed.
+
+## Run
+
+From repository root, with dependencies already available:
+
+    PYTHONPATH=src:. python -m unittest discover -s tests/unit -p test_shear_refinement_example.py -v
+
+Only on an approved idle CUDA device:
+
+    PYTHONPATH=src:. python -m research.shear_refinement.example
+    PYTHONPATH=src:.:tests/validation/multigrid python -m unittest discover -s tests/validation/multigrid -p test_shear_refinement_example_gpu.py -v
+
+The GPU gate compares automatic selection to an explicitly identical manual
+rectangle after short continuation, requires finite leaves, exact NaN masks,
+and an unchanged source. It does not replace independent numerical replay.
+GPU execution has not been performed in the implementation cloud environment.
+
+## Tagger provenance and license
+
+`boundary_shear_tagging.py` is the unchanged accepted Grok standalone v0.2
+preparation (2026-10-16), SHA256
+`12bc70ab989953501a752a2dc23e3712b6e69bfd47e67664c1771ff2dfa4da0b`.
+The delivered `README_boundary_shear_tagging_v0.2.md` identifies it as original
+code under MIT, the same license as CelerisLab (see repository LICENSE), with
+nothing copied from upstream examples. Its source ledger/review is in the
+preparation/grok_amr_review/v0.2 handoff. Vendoring keeps this example runnable
+without external preparation paths; it is not installed as production API.
+The separate reviewed level-request policy v0.2.1 is intentionally not copied:
+this example observes only L1 and can request only the existing L2 level.
+
+## Local implementation checks (2026-10-16)
+
+Seven focused example tests passed. The combined CPU suite below passed all
+60 tests, including existing topology, pinned-body and synchronized remap gates:
+
+    PYTHONPATH=src:. python -m unittest tests.unit.test_refinement tests.unit.test_refinement_three_level_multipatch tests.unit.test_pinned_body_remap tests.unit.test_runtime_child_remap tests.unit.test_synchronized_remap tests.unit.test_shear_refinement_example -q
+
+Broad unit discovery was also attempted; 14 legacy module imports failed
+because PyCUDA is absent. That broad run is not a pass. No installation or GPU
+execution was attempted. The example uses a passive Circle, with no Jet drive.

research/shear_refinement/boundary_shear_tagging.py

@@ -0,0 +1,285 @@
+"""Boundary-proximity + physically scaled shear tagging with hysteresis (CPU prototype), v0.2.
+
+Status: standalone preparation for CelerisLab AMR (baseline dev 93d1cdd). Original code, not
+copied from any project. It only *proposes* cells for refinement: it returns masks and
+diagnostics and never builds patches, regrids, applies a max_level, or touches solver state.
+Not a block/AMR integration: callers must still legalise/align proposals and run the existing
+BlockLayout2D / FixedPatch2D / prepare_fine_boundary guards.
+
+v0.2 changes (2026-10-16 review): explicit periodic axes for distance AND dilation (minimum
+image); three separate masks (physical solid / valid stencil samples / owned cells); core tags
+only on owned-eligible fluid; hold previous core where the strain is unevaluable.
+
+Masks (all bool, shape (ny, nx)):
+  * solid  - physical solid geometry. The ONLY source of walls. Never fill holes with it.
+  * valid  - cells holding a usable, finite macroscopic velocity sample (may include ghosts).
+             Derivative stencils read only valid cells. Default: ~solid. Must not overlap solid.
+             Cells that are neither solid nor valid are "absent" (sparse holes): never walls,
+             never derivative donors, never tagged; their velocity is never read.
+  * owned  - cells this caller owns and may propose (e.g. ACTIVE fine cells, owned blocks).
+             Default: all True. eligible = owned & ~solid, and must be a subset of valid.
+  Core tags exist only on eligible cells; the buffer is clipped to owned cells.
+
+Units (match CelerisLab refinement.py: coarse dx = dt = 1, level l has dx = dt = 2**-l):
+  * Lengths (origin, wall distance, L_ref, periods) are in COARSE cells.
+  * Velocities are lattice velocities; under acoustic scaling these equal coarse units.
+  * Strain rate |S| = sqrt(2 S_ij S_ij) in 1/coarse-time-step; u = (g*y, 0) gives |S| = |g|.
+  * Shear indicator s* = |S| * L_ref / U_ref: dimensionless; the physical definition is
+    level-independent (its finite-difference estimate is exactly level-invariant only for
+    affine fields; otherwise it changes by truncation error).
+  * cell_strain = |S| * dx_level / U_ref: an analogous dx-scaled STRAIN diagnostic, returned
+    but not used to tag. It is not the quantity waLBerla/AMReX threshold: those use
+    vorticity-related indicators (rigid rotation: zero strain, nonzero vorticity). We borrow
+    only their level-scaling pattern. A dx-scaled indicator is resolution-sensitive; it is not
+    an accuracy certificate and not a substitute for a hard max_level in the policy layer.
+  * Velocity must be MACROSCOPIC (moments of one consistent population phase).
+
+Tag rule per call, on eligible cells:
+  wall:  d <= band_on or (prev and d <= band_off)               band_off >= band_on (units L_ref)
+  shear: where evaluable: s* >= shear_on or (prev and s* >= shear_off)
+         where unevaluable (no valid stencil): HOLD prev (never switches a tag off)
+  core = wall | shear;  tagged = core dilated n_buffer cells (Chebyshev, periodic-aware) & owned.
+Feed ``core`` (never ``tagged``) back as ``prev_core``. prev_core outside eligible cells is
+rejected; buffered history on an all-fluid grid cannot be detected, so the contract stays
+core-only. History is only meaningful on unchanged topology: reset or remap after a regrid.
+An untagged unevaluable cell (reported in ``unevaluable``) is NOT evidence that it may coarsen.
+"""
+from dataclasses import dataclass
+from numbers import Integral, Real
+
+import numpy as np
+
+WALL, SHEAR, HELD = 1, 2, 4  # bits in TagResult.reason (HELD: shear state held, unevaluable)
+_CHUNK = 4096
+
+
+def _ro(a):
+    a = np.ascontiguousarray(a)
+    a.setflags(write=False)
+    return a
+
+
+def _num(name, v, *, positive=False, nonneg=False):
+    if isinstance(v, (bool, np.bool_)) or not isinstance(v, Real) or not np.isfinite(v):
+        raise ValueError(f"{name} must be a finite real number")
+    v = float(v)
+    if positive and v <= 0:
+        raise ValueError(f"{name} must be > 0")
+    if nonneg and v < 0:
+        raise ValueError(f"{name} must be >= 0")
+    return v
+
+
+def _periodic(periodic):
+    if (not isinstance(periodic, (tuple, list)) or len(periodic) != 2
+            or not all(isinstance(p, (bool, np.bool_)) for p in periodic)):
+        raise ValueError("periodic must be two bools (x, y)")
+    return bool(periodic[0]), bool(periodic[1])
+
+
+def _mask(name, m, shape, default):
+    if m is None:
+        return default
+    m = np.asarray(m)
+    if m.dtype != np.bool_ or m.shape != shape:
+        raise ValueError(f"{name} must be a bool array of the grid shape")
+    return m
+
+
+def _neighbour_or(m, periodic):
+    """OR over the 3x3 Chebyshev neighbourhood; wraps on periodic axes, False-pads others."""
+    px, py = periodic
+    p = np.pad(m, ((0 if py else 1,) * 2, (0 if px else 1,) * 2), constant_values=False)
+    out = np.zeros(m.shape, bool)
+    ny, nx = m.shape
+    for dy in (-1, 0, 1):
+        for dx in (-1, 0, 1):
+            if px or py:
+                s = np.roll(p, (-dy if py else 0, -dx if px else 0), axis=(0, 1))
+            else:
+                s = p
+            oy = 0 if py else 1 + dy
+            ox = 0 if px else 1 + dx
+            out |= s[oy:oy + ny, ox:ox + nx]
+    return out
+
+
+def dilate(mask, n, periodic=(False, False)):
+    """Chebyshev dilation by n cells; wraps across periodic axes only."""
+    periodic = _periodic(periodic)
+    m = np.asarray(mask, bool).copy()
+    for _ in range(int(n)):
+        m = _neighbour_or(m, periodic)
+    return m
+
+
+@dataclass(frozen=True, eq=False)
+class TagResult:
+    core: np.ndarray          # bool: hysteresis state on eligible cells; feed back as prev_core
+    tagged: np.ndarray        # bool: core + buffer, clipped to owned; proposal for a regridder
+    reason: np.ndarray        # uint8 bits WALL|SHEAR|HELD on core cells
+    strain_star: np.ndarray   # float64 s*; NaN where unevaluable or not valid
+    cell_strain: np.ndarray   # float64 |S| dx_level / U_ref (diagnostic only)
+    wall_distance: np.ndarray # float64 coarse cells from non-solid centres; inf if no wall / solid
+    unevaluable: np.ndarray   # bool: eligible cells with no usable gradient stencil
+    eligible: np.ndarray      # bool: owned & ~solid (the only cells that can hold core tags)
+    n_unevaluable_fluid: int  # count of `unevaluable`
+
+
+def strain_rate(ux, uy, fluid, dx_level, periodic=(False, False)):
+    """|S| = sqrt(2 S:S) in 1/coarse-time-step on `fluid` (= valid sample) cells; NaN elsewhere.
+
+    Each derivative uses a central difference when both neighbours are valid (inside the
+    domain or across a periodic seam), else a one-sided difference (first order) to one valid
+    neighbour, else NaN. Values in non-valid cells are never read.
+    """
+    fluid = np.asarray(fluid, bool)
+    px, py = _periodic(periodic)
+
+    def ddir(f, axis, per):
+        n = f.shape[axis]
+        idx = np.arange(n)
+        ip, im = idx + 1, idx - 1
+        okp, okm = ip < n, im >= 0
+        if per:
+            ip, im = ip % n, im % n
+            okp = okm = np.ones(n, bool)
+        ip, im = np.clip(ip, 0, n - 1), np.clip(im, 0, n - 1)
+        take = lambda a, i: np.take(a, i, axis=axis)
+        shape = [1, 1]; shape[axis] = n
+        fp = take(fluid, ip) & okp.reshape(shape)
+        fm = take(fluid, im) & okm.reshape(shape)
+        fz = np.where(fluid, f, 0.0)
+        vp, vm = take(fz, ip), take(fz, im)
+        out = np.full(f.shape, np.nan)
+        c = fluid & fp & fm
+        out[c] = 0.5 * (vp - vm)[c]
+        o = fluid & fp & ~fm
+        out[o] = (vp - fz)[o]
+        o = fluid & ~fp & fm
+        out[o] = (fz - vm)[o]
+        return out / dx_level
+
+    dudx, dudy = ddir(ux, 1, px), ddir(ux, 0, py)
+    dvdx, dvdy = ddir(uy, 1, px), ddir(uy, 0, py)
+    sxy = 0.5 * (dudy + dvdx)
+    return np.sqrt(2.0 * (dudx**2 + dvdy**2 + 2.0 * sxy**2))
+
+
+def wall_distance(fluid, dx_level, origin=(0.0, 0.0), distance_fn=None,
+                  periodic=(False, False), solid=None):
+    """Distance (coarse cells) from each non-solid cell centre to the wall; inf if no wall.
+
+    `solid` (default ~fluid) is the physical geometry; only solid cells with a non-solid
+    Chebyshev neighbour (seam-aware) are wall cells, so absent cells are never walls.
+    distance_fn(x, y) -> signed distance (coarse units, negative inside bodies) is preferred
+    when geometry is known; it is evaluated at cell centres as given, so for periodic axes it
+    must itself return the minimum-image distance. Otherwise a raster HEURISTIC is used:
+    distance to the nearest wall-cell centre (minimum image on periodic axes) minus dx/2.
+    It is exact only for a planar wall lying on cell faces; it has no general error bound
+    (a physical wall off the cell face, curved, sub-cell or off-patch bodies differ).
+    """
+    fluid = np.asarray(fluid, bool)
+    solid = ~fluid if solid is None else np.asarray(solid, bool)
+    px, py = _periodic(periodic)
+    ny, nx = fluid.shape
+    x = origin[0] + (np.arange(nx) + 0.5) * dx_level
+    y = origin[1] + (np.arange(ny) + 0.5) * dx_level
+    X, Y = np.meshgrid(x, y)
+    target = ~solid
+    d = np.full(fluid.shape, np.inf)
+    if distance_fn is not None:
+        dv = np.asarray(distance_fn(X, Y), dtype=np.float64)
+        if dv.shape != fluid.shape or not np.all(np.isfinite(dv[target])):
+            raise ValueError("distance_fn must return finite values of the grid shape")
+        d[target] = np.maximum(dv[target], 0.0)
+        return d
+    wall = solid & _neighbour_or(~solid, (px, py))
+    if not wall.any() or not target.any():
+        return d
+    Lx, Ly = nx * dx_level, ny * dx_level
+    wx, wy = X[wall], Y[wall]
+    fx, fy = X[target], Y[target]
+    best = np.empty(fx.size)
+    for s in range(0, fx.size, _CHUNK):
+        ex, ey = np.abs(fx[s:s + _CHUNK, None] - wx), np.abs(fy[s:s + _CHUNK, None] - wy)
+        if px:
+            ex = np.minimum(ex, Lx - ex)
+        if py:
+            ey = np.minimum(ey, Ly - ey)
+        best[s:s + _CHUNK] = np.sqrt((ex * ex + ey * ey).min(axis=1))
+    d[target] = np.maximum(best - 0.5 * dx_level, 0.0)
+    return d
+
+
+def tag_cells(ux, uy, solid, *, dx_level, U_ref, L_ref, shear_on, shear_off,
+              band_on, band_off, n_buffer=2, prev_core=None, origin=(0.0, 0.0),
+              periodic=(False, False), distance_fn=None, valid=None, owned=None):
+    """Propose refinement tags for one level. Pure function: inputs are not modified.
+
+    band_on/band_off are in units of L_ref. See module docstring for the three masks.
+    """
+    solid = np.asarray(solid)
+    if solid.dtype != np.bool_ or solid.ndim != 2 or min(solid.shape) < 3:
+        raise ValueError("solid must be a 2-D bool array with both sides >= 3")
+    shape = solid.shape
+    ux = np.asarray(ux, dtype=np.float64)
+    uy = np.asarray(uy, dtype=np.float64)
+    if ux.shape != shape or uy.shape != shape:
+        raise ValueError("ux, uy and solid must have the same (ny, nx) shape")
+    valid = _mask("valid", valid, shape, ~solid)
+    owned = _mask("owned", owned, shape, np.ones(shape, bool))
+    if np.any(valid & solid):
+        raise ValueError("valid must not overlap solid (solid cells are not velocity samples)")
+    eligible = owned & ~solid
+    if np.any(eligible & ~valid):
+        raise ValueError("owned non-solid cells must be valid (have a velocity sample)")
+    if not (np.all(np.isfinite(ux[valid])) and np.all(np.isfinite(uy[valid]))):
+        raise ValueError("velocity must be finite on every valid cell")
+    dx_level = _num("dx_level", dx_level, positive=True)
+    lvl = -np.log2(dx_level)
+    if abs(lvl - round(lvl)) > 1e-12 or round(lvl) < 0:
+        raise ValueError("dx_level must be 2**-level coarse cells (1, 0.5, 0.25, ...)")
+    U_ref = _num("U_ref", U_ref, positive=True)
+    L_ref = _num("L_ref", L_ref, positive=True)
+    shear_on = _num("shear_on", shear_on, positive=True)
+    shear_off = _num("shear_off", shear_off, nonneg=True)
+    if not shear_off < shear_on:
+        raise ValueError("hysteresis needs shear_off < shear_on")
+    band_on = _num("band_on", band_on, nonneg=True)
+    band_off = _num("band_off", band_off, nonneg=True)
+    if not band_off >= band_on:
+        raise ValueError("hysteresis needs band_off >= band_on")
+    if isinstance(n_buffer, (bool, np.bool_)) or not isinstance(n_buffer, Integral) or n_buffer < 0:
+        raise ValueError("n_buffer must be a non-negative integer")
+    if len(origin) != 2:
+        raise ValueError("origin must be (x0, y0) in coarse units")
+    origin = (_num("origin[0]", origin[0]), _num("origin[1]", origin[1]))
+    periodic = _periodic(periodic)
+    if prev_core is None:
+        prev = np.zeros(shape, bool)
+    else:
+        prev = np.asarray(prev_core)
+        if prev.dtype != np.bool_ or prev.shape != shape:
+            raise ValueError("prev_core must be a bool array of the grid shape")
+        if np.any(prev & ~eligible):
+            raise ValueError("prev_core marks solid/unowned cells: pass TagResult.core, not tagged")
+
+    S = strain_rate(ux, uy, valid, dx_level, periodic)
+    s_star = S * L_ref / U_ref
+    cell_strain = S * dx_level / U_ref
+    d = wall_distance(~solid, dx_level, origin, distance_fn, periodic, solid=solid)
+    dn = d / L_ref
+
+    ok = np.isfinite(s_star)
+    sv = np.where(ok, s_star, -np.inf)
+    unevaluable = eligible & ~ok
+    held = unevaluable & prev
+    shear = eligible & (((sv >= shear_on) | (prev & (sv >= shear_off))) | held)
+    wall = eligible & ((dn <= band_on) | (prev & (dn <= band_off)))
+    core = shear | wall
+    reason = (wall * WALL + shear * SHEAR + held * HELD).astype(np.uint8)
+    tagged = dilate(core, n_buffer, periodic) & owned
+    return TagResult(_ro(core), _ro(tagged), _ro(reason), _ro(s_star), _ro(cell_strain),
+                     _ro(d), _ro(unevaluable), _ro(eligible),
+                     int(np.count_nonzero(unevaluable)))

research/shear_refinement/example.py

@@ -0,0 +1,140 @@
+"""One measured L1 strain observation, one growth-only L2 remap, no scheduler.
+
+Run from the repository root: PYTHONPATH=src:. python -m research.shear_refinement.example
+Requires CUDA; CPU helpers and tests never import CUDA.
+"""
+from contextlib import ExitStack
+from numbers import Integral
+import numpy as np
+
+from CelerisLab.config import LBMConfig
+from CelerisLab.body.geometry.circle import CircleGeometry
+from CelerisLab.body.objects import SimObject
+from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
+from CelerisLab.lbm.refinement import FixedPatchSet2D
+from CelerisLab.lbm.synchronized_remap import RemapTopology
+from .boundary_shear_tagging import tag_cells
+
+CFG = LBMConfig(nx=48, ny=32, viscosity=.1)
+PARENTS = ((4, 4, 20, 28), (24, 4, 44, 28))
+CHILDREN = ((16, 16, 32, 40), (56, 16, 72, 40))
+CORRIDOR = (72, 20, 76, 36)  # L1 allocated indices; bounded to right of fluid L2
+
+
+def grow_rectangle(current, request, *, alignment=4, area_cap=480):
+    """Bounding-box growth only, in request-array cell-edge coordinates.
+
+    No body or solver knowledge. Reject invalid/budget-exceeding proposals with
+    ValueError; never clip a request, shrink old coverage or mutate an input.
+    Alignment is edge divisibility, not a substitute for full topology guards.
+    """
+    request = np.asarray(request)
+    if request.dtype != bool or request.ndim != 2:
+        raise ValueError('request must be a 2D boolean mask')
+    values = (*current, alignment, area_cap)
+    if len(current) != 4 or any(isinstance(v, (bool, np.bool_)) or not isinstance(v, Integral) for v in values):
+        raise ValueError('rectangle, alignment and area cap must be integers')
+    x0, y0, x1, y1 = current
+    ny, nx = request.shape
+    if alignment <= 0 or area_cap <= 0 or not (0 <= x0 < x1 <= nx and 0 <= y0 < y1 <= ny):
+        raise ValueError('invalid rectangle, alignment or area cap')
+    if any(v % alignment for v in current):
+        raise ValueError('current rectangle must be aligned')
+    y, x = np.nonzero(request)
+    if x.size:
+        x0, y0 = min(x0, int(x.min())), min(y0, int(y.min()))
+        x1, y1 = max(x1, int(x.max())+1), max(y1, int(y.max())+1)
+    a = alignment
+    result = (x0//a*a, y0//a*a, (x1+a-1)//a*a, (y1+a-1)//a*a)
+    if result[2] > nx or result[3] > ny:
+        raise ValueError('aligned growth exceeds observation frame')
+    if (result[2]-result[0])*(result[3]-result[1]) > area_cap:
+        raise ValueError('growth exceeds area cap')
+    return result
+
+
+def objects():
+    return (SimObject(71, geometry=CircleGeometry(12., 14., .8),
+                      center=(12., 14.), radius=.8),)
+
+
+def topology(rectangles=CHILDREN):
+    return RemapTopology.from_config(CFG, FixedPatchSet2D.from_config(CFG, PARENTS),
+                                     rectangles, child_objects=objects())
+
+
+def observe(g, mesh, *, shear_on=.001):
+    """Moments of completed-time postcollision G; only real L1 fluid leaves.
+
+    No reconstructed covered values, ghosts or holes are derivative donors.
+    solid is independent geometry (empty on L1: the body is pinned at L2).
+    The array is cropped/sparse, so neither derivative axis wraps.
+    """
+    valid = mesh.leaf_masks[1]
+    f = g[:, valid].astype(np.float64)
+    rho = f.sum(axis=0)
+    if not np.isfinite(f).all() or np.any(rho <= 0):
+        raise ValueError('invalid completed-time L1 fluid populations')
+    ux, uy = (np.full(valid.shape, np.nan) for _ in range(2))
+    ux[valid] = (np.asarray(D2Q9_EX)[:, None]*f).sum(axis=0)/rho
+    uy[valid] = (np.asarray(D2Q9_EY)[:, None]*f).sum(axis=0)/rho
+    tags = tag_cells(ux, uy, np.zeros(valid.shape, bool), valid=valid, owned=valid,
+                     dx_level=.5, origin=mesh.parent.fine_origin, periodic=(False, False),
+                     U_ref=.05, L_ref=1., shear_on=shear_on, shear_off=0.,
+                     band_on=0., band_off=0., n_buffer=0)
+    corridor = np.zeros(valid.shape, bool)
+    x0, y0, x1, y1 = CORRIDOR
+    corridor[y0:y1, x0:x1] = True
+    request = tags.tagged & corridor
+    proposed = (CHILDREN[0], grow_rectangle(CHILDREN[1], request))
+    # Host guards cover nesting, donor footprints, complete body and separation.
+    candidate = topology(proposed)
+    candidate.child.fine_layout((4, 4))
+    report = dict(requested=int(request.sum()), unknown=int(tags.unevaluable.sum()),
+                  corridor_unknown=int((tags.unevaluable & corridor).sum()),
+                  corridor_unavailable=int((corridor & ~valid).sum()),
+                  out_of_scope=int((tags.tagged & ~corridor).sum()), rectangles=proposed)
+    return proposed, report
+
+
+def initial_fields(mesh):
+    """Smooth periodic shear wave, sampled in each level's physical frame."""
+    from CelerisLab.lbm.equilibrium import compute_feq_d2q9
+    frames = ((np.zeros(2), 1.), (np.asarray(mesh.parent.fine_origin), .5),
+              (np.asarray(mesh.parent.fine_origin)+.5*np.asarray(mesh.child.fine_origin), .25))
+    arrays = []
+    for leaf, (origin, dx) in zip(mesh.leaf_masks, frames):
+        y = origin[1]+dx*(np.arange(leaf.shape[0])+.5)
+        ux = .01*np.sin(2*np.pi*y/CFG.ny)
+        rows = np.stack([compute_feq_d2q9(1., u, 0.) for u in ux], axis=1)
+        arrays.append(np.broadcast_to(rows[:, :, None], (9, *leaf.shape)).copy())
+    return tuple(arrays)
+
+
+def initialize(run):
+    run.initialize(*initial_fields(topology()))
+
+
+def make(context):
+    from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+    return FixedPatchStepper(CFG, PARENTS, context=context, base_block=(4, 4),
+                             fine_block=(4, 4), child_block=(4, 4),
+                             child_rectangles_in_parent=CHILDREN, child_objects=objects())
+
+
+def main():
+    from CelerisLab.cuda.context import CudaContext
+    with CudaContext(0) as ctx, ExitStack() as owned:
+        run = owned.enter_context(make(ctx))
+        initialize(run)
+        run.step(2)  # Returns only at completed synchronized base time.
+        proposed, report = observe(run.download()[1], topology())
+        print(report)
+        if proposed != CHILDREN:
+            run = owned.enter_context(run.remap_child(proposed, child_objects=objects()))
+        run.step(2)
+        print(f'completed base step {run.step_count}; body rectangle {run.child_patch.bounds[0]}')
+
+
+if __name__ == '__main__':
+    main()

src/CelerisLab/body/jets.py

@@ -9,6 +9,11 @@ whose links are evaluated; callers convert frames/units before evaluation.
 One physical link (receiver, direction) may be driven by at most one jet, and
 one jet may drive the links of only one body (checked at assignment, any amplitude).
 
+``condition`` names the prescribed boundary type. Only ``"velocity"`` (the
+default: amplitude is a signed lattice velocity, not a flow rate) is
+implemented; ``"pressure"`` is reserved and rejected at construction, before
+any amplitude validation, so it never reaches assignment or a kernel.
+
 Edge rule: the lateral bound |s| < 1 is exclusive (the profile vanishes
 there for parabolic/quartic); Sector r_min/r_max and Tube 0 <= along <= length
 are inclusive. Jets are immutable; an amplitude change replaces the Jet.
@@ -25,6 +30,9 @@ from dataclasses import dataclass, replace
 
 import numpy as np
 
+CONDITIONS = ("velocity",)
+RESERVED_CONDITIONS = ("pressure",)
+
 PROFILES = {
     "uniform": lambda s: np.ones_like(s),
     "parabolic": lambda s: 1 - s*s,
@@ -122,7 +130,8 @@ class Parallel:
 
 @dataclass(frozen=True)
 class Jet:
-    """Region (Sector|Tube), field (Radial|Parallel), profile name, signed amplitude.
+    """Region (Sector|Tube), field (Radial|Parallel), profile name, signed amplitude,
+    prescribed boundary condition (velocity only; pressure reserved).
 
     Frozen: change amplitude with ``with_amplitude`` (validated, new object);
     the cached per-link shapes stay valid because region/field/profile cannot change.
@@ -131,8 +140,14 @@ class Jet:
     field: object
     profile: str = "uniform"
     amplitude: float = 0.0
+    condition: str = "velocity"
 
     def __post_init__(self):
+        # The condition decides what amplitude and field mean: check it first.
+        if self.condition in RESERVED_CONDITIONS:
+            raise NotImplementedError(f"{self.condition} jets are reserved but not implemented")
+        if self.condition not in CONDITIONS:
+            raise ValueError(f"Jet condition must be one of {CONDITIONS + RESERVED_CONDITIONS}")
         if self.profile not in PROFILES:
             raise ValueError(f"Jet profile must be one of {sorted(PROFILES)}")
         _finite(self.amplitude)

src/CelerisLab/lbm/block_stepper.py

@@ -211,6 +211,7 @@ class PeriodicBlockStepper:
                         t = self._tables
                         self._curved(self._current, t["block"], t["cell"], t["dir"], t["q"],
                                      t["fallback"], np.uint32(self._links), np.uintp(0),  # no actuation
+                                     np.uintp(0), np.uintp(0), np.uintp(0),  # no observation
                                      grid=self._link_grid, **kw)
                     self._step(self._current, self._spare, np.uint32(self.layout.halo),
                                self._tables.get("flags", np.uintp(0)), np.uintp(0),  # no kind

src/CelerisLab/lbm/hierarchy.py

@@ -0,0 +1,100 @@
+# CelerisLab/lbm/hierarchy.py
+"""Host-only geometry of a fixed 2:1 refinement hierarchy of any depth (numpy, no CUDA).
+
+Loops the existing adjacent-pair logic of refinement.py; it adds no transfer rule.
+``refinements[l]`` describes level l+1 in level-l integer cells: a scalar
+(x0, y0, x1, y1) is one FixedPatch2D window (its level keeps the window shape), a
+sequence of rectangles is a FixedPatchSet2D on the level-wide grid (origin 0).
+A scalar child is only accepted under a scalar parent. Every pair after the first
+is built by nest_patch/nest_patch_set, whose _nested_kind validates coverage,
+donors/receivers, surviving pulls and the parent's upward restriction footprint.
+"""
+from dataclasses import dataclass, replace
+from numbers import Integral
+
+import numpy as np
+
+from .refinement import (ACTIVE, OUTER_GHOST, FixedPatch2D, FixedPatchSet2D,
+                         _middle_config, nest_patch, nest_patch_set)
+
+
+@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)."""
+    configs: tuple
+    pairs: tuple
+    kinds: tuple
+
+    @property
+    def depth(self):
+        return len(self.kinds)
+
+
+def is_rectangle(value):
+    value = tuple(value)
+    return len(value) == 4 and all(isinstance(v, (Integral, np.integer))
+                                   and not isinstance(v, (bool, np.bool_)) for v in value)
+
+
+def build_hierarchy(cfg, refinements, *, wall_aware=False):
+    """Geometry for len(refinements) + 1 levels; wall_aware marks base rows 0 and -1 OUTER."""
+    configs, pairs, kinds = [_validated(cfg)], [], []
+    for level, entry in enumerate(refinements):
+        scalar = is_rectangle(entry)
+        if not scalar and not tuple(entry):
+            raise ValueError("every refinement level needs at least one rectangle")
+        if level == 0:
+            pair = (FixedPatch2D if scalar else FixedPatchSet2D).from_config(configs[0], entry)
+            kind = np.where(pair.covered, OUTER_GHOST, ACTIVE).astype(np.uint8)
+            if wall_aware:
+                kind[[0, -1]] = OUTER_GHOST
+        elif scalar:
+            if not isinstance(pairs[-1], FixedPatch2D):
+                raise ValueError("a scalar child requires a scalar parent; use a rectangle list")
+            pair, kind = nest_patch(configs[level - 1], pairs[-1], entry)
+        else:
+            pair, kind = nest_patch_set(configs[level - 1], pairs[-1], entry)
+        kinds.append(np.asarray(kind))
+        pairs.append(pair)
+        configs.append(_validated(_middle_config(configs[level], pair)))
+        _check_chain(configs[level], configs[level + 1], pair)
+    if pairs:
+        kinds.append(pairs[-1].fine_kind)
+    else:
+        kind = np.full((cfg.ny, cfg.nx), ACTIVE, np.uint8)
+        if wall_aware:
+            kind[[0, -1]] = OUTER_GHOST
+        kinds.append(kind)
+    return Hierarchy(tuple(configs), tuple(pairs), tuple(kinds))
+
+
+def _validated(cfg):
+    """LBMConfig.omega is computed only by validate(); replace() copies it stale."""
+    cfg = replace(cfg)
+    cfg.validate()
+    return cfg
+
+
+def _check_chain(coarse, fine, pair):
+    if not (np.isclose(pair.tau_coarse, 1 / coarse.omega, rtol=1e-12, atol=0)
+            and np.isclose(pair.tau_fine, 1 / fine.omega, rtol=1e-12, atol=0)
+            and tuple(pair.fine_shape) == (fine.ny, fine.nx)):
+        raise RuntimeError("level config chain is inconsistent with the pair tables")
+
+
+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.
+    """
+    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))
+            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)
+    return base

src/CelerisLab/lbm/kernels/blocks.cu

@@ -76,11 +76,12 @@ extern "C" __global__ void BlockCoalesce(
 // cell, so the existing BFL writes the solid-source window image that this
 // block's next pull reads (no halo copy between this kernel and BlockStep).
 // surface_velocity (nullable, 2 floats per link): prescribed fluid velocity on
-// fixed geometry, same convention as CurvedBoundaryKernel. No exchange output.
+// fixed geometry, same convention as CurvedBoundaryKernel. Optional raw rows
+// (jx,jy,torque,dm) are unique physical links, never body-observation slots.
 extern "C" __global__ void BlockCurvedBoundary(
     fpxx* f, const uint32_t* block, const uint32_t* cell, const unsigned char* dir,
     const float* q, const unsigned char* fallback, unsigned int count,
-    const float* surface_velocity)
+    const float* surface_velocity, const float* rx, const float* ry, float* raw)
 {
     unsigned int row = blockIdx.x * blockDim.x + threadIdx.x;
     if (row >= count) return;
@@ -90,7 +91,24 @@ extern "C" __global__ void BlockCurvedBoundary(
         if (u != 0.0f || v != 0.0f) { Uw = u; Vw = v; }
     }
     apply_bouzidi_link(dir[row], q[row], cell[row], f + (unsigned long)block[row] * NQ * TOTAL_CELLS,
-                       Uw, Vw, 0.0f, 0.0f, fallback[row], nullptr, -1);
+                       Uw, Vw, rx ? rx[row] : 0.0f, ry ? ry[row] : 0.0f,
+                       fallback[row], nullptr, -1, raw ? raw + 4ul * row : nullptr,
+                       raw ? raw + 4ul * row + 3 : nullptr);
+}
+
+// Same-stream BFL -> this row-owned accumulation -> pull step. No atomics,
+// no cross-row reduction and no history allocation proportional to substeps.
+extern "C" __global__ void BlockAccumulateExchange(
+    const float* raw, double* integral, unsigned int count, double h, double dt)
+{
+    unsigned int row = blockIdx.x * blockDim.x + threadIdx.x;
+    if (row >= count) return;
+    const unsigned long k = 4ul * row;
+    const double impulse = h * dt;
+    integral[k] += impulse * (double)raw[k];
+    integral[k + 1] += impulse * (double)raw[k + 1];
+    integral[k + 2] += h * impulse * (double)raw[k + 2];
+    integral[k + 3] += h * h * (double)raw[k + 3];
 }
 
 // Pull/collide/store core cells only; halo cells of out are left untouched.

src/CelerisLab/lbm/kernels/boundary/curved_boundary.cuh

@@ -82,7 +82,8 @@ __device__ inline void apply_bouzidi_link(
     float Uw, float Vw,
     float rx, float ry,
     unsigned char fallback_class,
-    float* __restrict__ obs, int body_id, float* link_exchange = nullptr)
+    float* __restrict__ obs, int body_id, float* link_exchange = nullptr,
+    float* link_mass = nullptr)
 {
     const unsigned int dir_opp = (unsigned int)opp_dir((int)dir);
 
@@ -120,6 +121,10 @@ __device__ inline void apply_bouzidi_link(
     // corrected incoming population into the adjacent fluid node.
     store_ddf(fi, index_f(k_s, dir_opp), f_reflected);
 
+    // Equal opposite weights cancel even with shifted DDFs.
+    if (link_mass != nullptr)
+        *link_mass = f_reflected - f_toward;
+
     if (obs != nullptr || link_exchange != nullptr) {
 #if USE_DDF_SHIFTING
         const float w_i = (float)d_w[dir];

src/CelerisLab/lbm/raw_exchange.py

@@ -0,0 +1,36 @@
+"""Host-only reduction of raw cut-link observations; no stress/GI correction.
+
+Rows are unique physical links. All quantities are 2D, per unit out-of-plane
+thickness; centers are in the root edge frame, torque is about the named body.
+"""
+import numpy as np
+
+
+def summarize_exchange(rows, body_ids, aperture, *, body_id, level, h,
+                       time_end, window, reference_center):
+    """Reduce FP64 interval integrals or convert FP32 last-substep raw rows.
+
+    dt_l=h_l under acoustic 2:1 refinement. Local momentum j -> base force
+    h*j -> impulse h*dt*j; torque adds one h; mass is h**2*dm.
+    """
+    if window not in ("last_base_step", "last_substep"):
+        raise ValueError("unknown exchange window")
+    values = np.array(rows, dtype=np.float64, copy=True)
+    duration = 1.0 if window == "last_base_step" else h
+    if window == "last_substep":
+        values *= np.array([h*h, h*h, h*h*h, h*h])
+    selected = np.asarray(body_ids) == body_id
+    aperture = np.asarray(aperture, bool)
+
+    def reduce(mask):
+        integral = values[mask].sum(axis=0, dtype=np.float64)
+        return dict(impulse=integral[:2].copy(), angular_impulse=float(integral[2]),
+                    mass_into_fluid=float(integral[3]), force=integral[:2] / duration,
+                    torque=float(integral[2] / duration), mass_rate=float(integral[3] / duration))
+
+    return dict(wall=reduce(selected & ~aperture), aperture=reduce(selected & aperture),
+                total=reduce(selected), body_id=int(body_id), level=int(level),
+                time_start=time_end-duration, time_end=time_end, dt=h,
+                sample_count=round(1/h) if window == "last_base_step" else 1,
+                window=window, units="base_lattice_per_unit_depth",
+                definition="raw_boundary_exchange", reference_center=tuple(reference_center), valid=True)

src/CelerisLab/lbm/refinement.py

@@ -227,6 +227,11 @@ class FixedPatchSet2D:
     def bounds(self) -> tuple:
         return tuple(p.bounds for p in self.patches)
 
+    @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 covered(self) -> np.ndarray:
         mask = np.zeros(self.coarse_shape, bool)
@@ -237,7 +242,11 @@ class FixedPatchSet2D:
     def fine_layout(self, block):
         """Sparse non-periodic layout storing exactly the patch supports."""
         from .blocks import BlockLayout2D
-        by, bx = block
+        block = tuple(block)
+        if len(block) != 2 or any(isinstance(v, (bool, np.bool_)) or not isinstance(v, Integral)
+                                  or v <= 0 for v in block):
+            raise ValueError("fine block must be two positive integers (by, bx)")
+        by, bx = map(int, block)
         origins = set()
         for p in self.patches:
             y0, x0 = p.fine_level_origin
@@ -295,17 +304,98 @@ def _transfer_footprint(patch):
     return fine, coarse, _readonly(sorted(reads), np.int64)
 
 
-def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=()):
+def nest_patch(cfg, parent, child_bounds):
+    """Child FixedPatch2D inside a fine patch window, plus the composed middle kind.
+
+    ``cfg`` is the base (L0) configuration and ``parent`` the L0->L1
+    FixedPatch2D; ``child_bounds`` (x0, y0, x1, y1) are in L1 window cells
+    (its dense fine grid, ghost band included). Returns (child, kind): the
+    L1->L2 FixedPatch2D (viscosity 2*nu on L1) and the L1 kind array with
+    inherited INNER/OUTER kept and only the child-covered ACTIVE cells made
+    OUTER (invalid). Rejects, from the actual tables before any allocation:
+    coverage outside inherited ACTIVE; explosion donors or coalescence
+    receivers that are not uncovered ACTIVE; any surviving ACTIVE pull that
+    would read a covered or OUTER cell without a child-supplied arrival; and
+    child coverage of L1's own two-hop restriction footprint toward L0.
+    """
+    child = FixedPatch2D.from_config(_middle_config(cfg, parent), child_bounds)
+    return child, _nested_kind(parent, child)
+
+
+def nest_patch_set(cfg, parent, child_rectangles_in_parent):
+    """Separated child union in the parent's allocated L1 array coordinates.
+
+    The root edge coordinate of L1 index (x, y) is
+    ``parent.fine_origin + (x, y)/2``. For a set parent its origin is zero;
+    a scalar parent's window includes its ghost band and has a shifted origin.
+    The child set allocates a dense public extent twice that L1 array, while
+    its blocked storage omits holes. Both adjacent pairs retain gap >= 4 in
+    their respective parent-cell units. Validate support over the whole union.
+    """
+    child = FixedPatchSet2D.from_config(
+        _middle_config(cfg, parent), child_rectangles_in_parent)
+    return child, _nested_kind(parent, child)
+
+
+def _middle_config(cfg, parent):
+    fy, fx = parent.fine_shape
+    return replace(cfg, nx=fx, ny=fy, viscosity=2 * cfg.viscosity, y_wall_bc="bounce_back",
+                   top_wall_bc="velocity", wall_slot_width=0, wall_normal_velocity=0.0,
+                   velocity=0.0)
+
+
+def _nested_kind(parent, child):
+    """Validate composed transfer support and preserve inherited ghost kinds."""
+    fy, fx = parent.fine_shape
+    kind = parent.fine_kind
+    covered = child.covered
+    active = (kind == ACTIVE) & ~covered
+    if not (kind[covered] == ACTIVE).all():
+        raise ValueError("child must cover only inherited ACTIVE cells")
+    cells = fx * fy
+    receivers = child.coalesce_coarse % np.uint64(cells)
+    if (not active.ravel()[child.ghost_coarse].all()
+            or not active.ravel()[receivers.astype(np.intp)].all()):
+        raise ValueError("child donors/receivers must be uncovered parent ACTIVE cells")
+    supplied = set(map(int, child.coalesce_coarse))
+    invalid = covered | (kind == OUTER_GHOST)
+    for y, x in zip(*np.nonzero(active)):
+        for q in range(1, 9):
+            sy, sx = y - D2Q9_EY[q], x - D2Q9_EX[q]
+            inside = 0 <= sy < fy and 0 <= sx < fx
+            if (not inside or invalid[sy, sx]) and q * cells + y * fx + x not in supplied:
+                raise ValueError("a surviving parent ACTIVE pull would read an invalid cell")
+    footprint, _, _ = _transfer_footprint(parent)
+    if covered.ravel()[sorted(footprint)].any():
+        raise ValueError("child covers the parent's restriction footprint")
+    composed = kind.copy()
+    composed[covered] = OUTER_GHOST
+    return _readonly(composed, np.uint8)
+
+
+def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=(),
+                          parent_origin=(0., 0.), parent_spacing=1.):
     """Project fixed circles and reject boundary/transfer overlap before CUDA.
 
-    Input centers/radii are in coarse physical coordinates (dx=dt=1).
+    Input centers/radii and Jets are in the root edge frame. ``parent_origin``
+    and ``parent_spacing`` map this patch's parent grid into that frame; the
+    default preserves the L0->L1 convention. cfg must describe the parent grid.
     SurfaceActuation retains its lattice-velocity meaning: input coarse units
-    are converted into independent fine-local objects with dx=dt=.5.
+    are converted into independent fine-local objects with
+    dx=dt=parent_spacing/2.
     Original objects are not retained or mutated. No runtime geometry updates.
+    ``patch`` may be a FixedPatchSet2D: objects are then projected once onto
+    the level-wide fine grid (fine_origin (0, 0)); each circle must lie inside
+    one patch, body ids are the positions in ``objects`` (never patch ids), and
+    jets and the (receiver, dir) writer checks run once over the whole level.
     ``jets`` use the same coarse cell-edge frame as the objects; their
     selection is checked against true fine hits here. Jet amplitudes need no
     rescaling (lattice velocity scales by dt/dx = 1 under acoustic 2:1).
     """
+    origin = np.asarray(parent_origin, dtype=np.float64)
+    if (origin.shape != (2,) or not np.isfinite(origin).all()
+            or not np.isfinite(parent_spacing) or parent_spacing <= 0):
+        raise ValueError("Parent frame must have a finite origin and positive spacing")
     if not objects:
         if jets:
             raise ValueError("Fine jets need fine boundary objects")
@@ -326,8 +416,10 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=()):
     flags = np.full(nx * ny, FLUID, np.uint16)
     active = patch.fine_kind.ravel() == ACTIVE
     local, links = [], []
-    x0, y0, x1, y1 = patch.bounds
-    dx = dt = .5
+    rects = patch.bounds if isinstance(patch, FixedPatchSet2D) else (patch.bounds,)
+    covered = patch.covered
+    dx = dt = .5 * parent_spacing
+    fine_origin = origin + parent_spacing * np.asarray(patch.fine_origin)
     for bid, obj in enumerate(objects):
         if (type(obj.geometry) is not CircleGeometry or obj.is_sensor or obj.is_force_region
                 or any(getattr(obj.state, name) != 0 for name in ("vx", "vy", "omega"))
@@ -335,10 +427,13 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=()):
             raise ValueError("Fine boundaries require stationary circles without other roles")
         cx, cy = obj.geometry.center
         radius = obj.geometry._radius
+        px, py = (np.asarray((cx, cy)) - origin) / parent_spacing
+        pr = radius / parent_spacing
         if (not np.isfinite([cx, cy, radius]).all() or radius <= 0
-                or cx-radius < x0 or cx+radius > x1 or cy-radius < y0 or cy+radius > y1):
-            raise ValueError("Entire fine circle must lie inside the physical patch")
-        center = tuple((v-o)/dx-.5 for v, o in zip((cx, cy), patch.fine_origin))
+                or not any(x0 <= px-pr and px+pr <= x1 and y0 <= py-pr
+                           and py+pr <= y1 for x0, y0, x1, y1 in rects)):
+            raise ValueError("Entire fine circle must lie inside one physical patch")
+        center = tuple((v-o)/dx-.5 for v, o in zip((cx, cy), fine_origin))
         spec = obj.surface_actuation
         if spec is not None:
             spec = replace(spec, normal_velocity=spec.normal_velocity * dt / dx)
@@ -397,7 +492,7 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=()):
         raise ValueError("Boundary stencil overlaps the two-hop transfer footprint")
     for k in coarse:
         y, x = divmod(k, cfg.nx)
-        if x0 <= x < x1 and y0 <= y < y1:
+        if covered[y, x]:
             raise ValueError("Transfer coarse donor/target is not an uncovered bulk cell")
         # Rows 0/ny-1 are never computed and rows 1/ny-2 are reconstructed after
         # coalescence; rows 2/ny-3 are ordinary bulk cells only read as donors.
@@ -405,7 +500,7 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=()):
             raise ValueError("Transfer footprint touches the coarse wall reconstruction stencil")
     from ..body.jets import fine_hits_to_edge_frame, surface_link_velocity
     hits = fine_hits_to_edge_frame([(link.hit_x, link.hit_y) for link in links],
-                                   patch.fine_origin, dx)
+                                   fine_origin, dx)
     packed = pack_cut_links_to_soa(links)
     if jets:
         # Same assignment checks as the stepper: empty/overlapping jets and

src/CelerisLab/lbm/refinement_tables.py

@@ -0,0 +1,107 @@
+# CelerisLab/lbm/refinement_tables.py
+"""Host-only (numpy) mapping of fixed-patch transfer and boundary tables to block storage.
+
+No CUDA import: CPU tests import these directly. refinement_stepper re-exports the names.
+"""
+from dataclasses import dataclass, replace
+
+import numpy as np
+
+from .blocks import blocked_link_tables
+from ..body.jets import jet_link_shapes, surface_link_velocity
+from .descriptors import D2Q9_EX, D2Q9_EY
+from .refinement import ACTIVE, OUTER_GHOST
+
+
+def blocked_fine_tables(patch, layout, order):
+    """Remap fixed-patch fine tables to owner-core blocked storage (host only).
+
+    Returns (explosion target block cells, coalescence child slots in the same
+    rows/child order as ``patch.coalesce_fine[order]``, per-block-cell kind with
+    unused halo entries OUTER). Coarse donors/targets are unchanged.
+    """
+    owner = layout.owner.ravel().astype(np.intp)
+    cells = patch.fine_shape[0] * patch.fine_shape[1]
+    q, cell = np.divmod(patch.coalesce_fine[order], np.uint64(cells))
+    kind = blocked_kind(layout, patch.fine_kind)
+    return (owner[patch.ghost_fine.astype(np.intp)].astype(np.uint64),
+            layout.slot(owner[cell.astype(np.intp)], q), kind)
+
+
+def blocked_kind(layout, dense_kind):
+    """Per-block-cell kind; unused halo entries and unstored holes stay OUTER."""
+    owner = layout.owner.ravel().astype(np.intp)
+    owned = layout.owned.ravel()  # sparse multipatch levels leave holes unstored
+    kind = np.full(layout.nblocks * int(np.prod(layout.local_shape)), OUTER_GHOST, np.uint8)
+    kind[owner[owned]] = np.asarray(dense_kind).ravel()[owned]
+    return kind
+
+def blocked_base_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
+    before parent pull. Source *owner* slots are wrong at block seams.
+    Nested pairs must pass composed_kind to retain inherited INNER/OUTER and
+    sparse holes; donors and receivers must be owned, uncovered ACTIVE cells.
+    """
+    owner = layout.owner.ravel()
+    q, receiver = np.divmod(patch.coalesce_coarse[order], np.uint64(owner.size))
+    if composed_kind is None:
+        dense_kind = np.where(patch.covered, OUTER_GHOST, ACTIVE).astype(np.uint8)
+        if wall_aware:
+            dense_kind[[0, -1]] = OUTER_GHOST
+    else:
+        dense_kind = np.asarray(composed_kind)
+        if dense_kind.shape != patch.coarse_shape:
+            raise ValueError("composed parent kind has the wrong shape")
+    donors = patch.ghost_coarse.astype(np.intp)
+    receivers = receiver.astype(np.intp)
+    legal = layout.owned.ravel() & (dense_kind.ravel() == ACTIVE) & ~patch.covered.ravel()
+    if not legal[donors].all() or not legal[receivers].all():
+        raise ValueError("child donors/receivers need owned uncovered parent ACTIVE cells")
+    offset = (np.asarray(D2Q9_EY)[q.astype(np.intp)] * layout.local_shape[1]
+              + np.asarray(D2Q9_EX)[q.astype(np.intp)])
+    source = owner[receivers].astype(np.int64) - offset
+    return (owner[donors].astype(np.uint64), layout.slot(source, q),
+            blocked_kind(layout, dense_kind))
+
+
+# Required device columns plus optional actuation; observation arms are uploaded
+# separately only when enabled. body_id/aperture remain host reduction metadata.
+BLOCKED_BOUNDARY_DEVICE_COLUMNS = ("block", "cell", "dir", "q", "fallback", "surface_velocity")
+
+
+@dataclass
+class _BlockedBoundary:
+    """One level's stationary BFL data and mutable Jet velocity column."""
+    host: object
+    tables: dict
+    jets: tuple
+    actuated: bool
+    jet_cache: object = None
+    flags: object = None
+    device: object = None
+    kernel: object = None
+    last_raw: object = None
+    interval_integral: object = None
+    accumulate: object = None
+
+
+def _prepare_blocked_boundary(boundary, layout, cfg, jets):
+    if boundary is None or layout is None:
+        return None
+    fluid_idx, dirs, q, body_id, rx, ry, _, fallback = boundary.curved_host
+    tables, _ = blocked_link_tables(layout, boundary.flags, fluid_idx, dirs, q, fallback)
+    # blocked_link_tables/map_cut_links already rejects repeated (fluid_idx, dir),
+    # independently of body_id, and repeated local writers.
+    tables.update(body_id=np.asarray(body_id, np.int32), rx=np.asarray(rx, np.float32),
+                  ry=np.asarray(ry, np.float32), aperture=np.zeros(len(dirs), bool))
+    actuated = bool(jets) or any(o.surface_actuation is not None for o in boundary.objects)
+    if actuated:
+        cfg = replace(cfg)
+        cfg.validate()
+        tables["surface_velocity"], tables["aperture"] = surface_link_velocity(
+            boundary.objects, cfg, boundary.flags, fluid_idx, dirs, body_id, rx, ry,
+            jets=jets, hits=boundary.curved_hits)
+    return _BlockedBoundary(boundary, tables, jets, actuated,
+                            jet_link_shapes(jets, boundary.curved_hits) if jets else None)

src/CelerisLab/lbm/synchronized_remap.py

@@ -0,0 +1,213 @@
+"""CPU-only whole-volume G transfer, not a runtime adaptive solver.
+
+Periodic three-level grids with fixed L0/L1 topology are supported, including
+stationary circles pinned inside unchanged L2 rectangles. The caller
+certifies that current owned leaves are quiescent,
+unshifted D2Q9 postcollision G at the supplied completed base interval.
+No device buffers, device boundary state, ghosts, collision or time advance
+occur; boundary preparation is host-only.
+"""
+from dataclasses import dataclass, replace
+from numbers import Integral
+
+import numpy as np
+
+from .refinement import (ACTIVE, FixedPatch2D, FixedPatchSet2D, FineBoundary2D,
+                         nest_patch_set, prepare_fine_boundary)
+from .descriptors import SOLID
+from ..body.jets import surface_link_velocity
+
+
+@dataclass(frozen=True, eq=False)
+class RemapTopology:
+    """Validated fixed parent and separated L2 set; use ``from_config``.
+
+    Child rectangles are in the parent's allocated L1 coordinates, including
+    its window offset/halo when the parent is a scalar patch. Dense masks are
+    rebuilt from geometry, never inferred from finite values in input arrays.
+    ``pinned_rectangles`` canonically lists each body-containing rectangle once;
+    fluid-only layouts have no pins.
+    """
+    parent: object
+    child: FixedPatchSet2D
+    leaf_masks: tuple
+    boundary: FineBoundary2D | None = None
+    pinned_rectangles: tuple = ()
+    surface_velocity: np.ndarray | None = None
+
+    @classmethod
+    def from_config(cls, cfg, parent, child_rectangles, *, coarse_boundary="periodic", objects=(),
+                    child_objects=(), child_jets=()):
+        if coarse_boundary != "periodic" or tuple(objects) or cfg.nz != 1:
+            raise ValueError("remap requires periodic 2D state; use child_objects for a pinned L2 body")
+        if type(parent) not in (FixedPatch2D, FixedPatchSet2D):
+            raise ValueError("parent must be a fixed patch or fixed patch set")
+        parent = type(parent).from_config(cfg, parent.bounds)
+        child, middle_kind = nest_patch_set(cfg, parent, child_rectangles)
+        child_objects, child_jets = tuple(child_objects), tuple(child_jets)
+        fine_cfg = replace(cfg, nx=parent.fine_shape[1], ny=parent.fine_shape[0],
+                           viscosity=2 * cfg.viscosity)
+        boundary = prepare_fine_boundary(
+            fine_cfg, child, child_objects, jets=child_jets,
+            parent_origin=parent.fine_origin, parent_spacing=.5)
+        pinned, velocity = (), None
+        child_leaf = child.fine_kind == ACTIVE
+        if boundary is not None:
+            child_leaf &= (boundary.flags.reshape(child.fine_shape) & SOLID) == 0
+            # Whole circle containment, not merely the resolved solid cells.
+            containing = set()
+            for obj in child_objects:
+                cx, cy = (np.asarray(obj.geometry.center) - parent.fine_origin) * 2
+                radius = obj.geometry._radius * 2
+                containing.add(next(r for r in child.bounds if
+                                    r[0] <= cx-radius and cx+radius <= r[2] and
+                                    r[1] <= cy-radius and cy+radius <= r[3]))
+            pinned = tuple(sorted(containing))
+            velocity = boundary_velocity(cfg, child, boundary, child_jets)
+            if velocity is not None:
+                velocity.setflags(write=False)
+        masks = (~parent.covered, middle_kind == ACTIVE, child_leaf)
+        for mask in masks:
+            mask.setflags(write=False)
+        return cls(parent, child, masks, boundary, pinned, velocity)
+
+
+def boundary_velocity(cfg, child, boundary, jets=()):
+    """Physical cut-link order; never blocked addresses or stale Jet snapshots."""
+    local_cfg = replace(cfg, nx=child.fine_shape[1], ny=child.fine_shape[0],
+                        viscosity=4 * cfg.viscosity)
+    local_cfg.validate()
+    fluid, dirs, _, ids, rx, ry, _, _ = boundary.curved_host
+    velocity, _ = surface_link_velocity(
+        boundary.objects, local_cfg, boundary.flags, fluid, dirs, ids, rx, ry,
+        jets=jets, hits=boundary.curved_hits)
+    return velocity
+
+
+def same_physical_boundary(a, b):
+    """Compare geometry, ownership, hits and the entire BFL support stencil."""
+    if a is None or b is None:
+        return a is b
+    if not isinstance(a, FineBoundary2D) or not isinstance(b, FineBoundary2D):
+        return False
+    # Exact ordered geometry also catches changes hidden by cut-table rounding.
+    def signatures(boundary):
+        return tuple((o.geometry.center, o.geometry._radius, o.surface_actuation)
+                     for o in boundary.objects)
+    if signatures(a) != signatures(b):
+        return False
+    return all(np.array_equal(x, y) for x, y in zip(
+        (a.flags, a.boundary_cells, a.curved_hits, *a.curved_host),
+        (b.flags, b.boundary_cells, b.curved_hits, *b.curved_host)))
+
+
+@dataclass(frozen=True, eq=False)
+class RemapCandidate:
+    """Separate dense arrays plus explicit synchronization/frame metadata.
+
+    Passing these arrays to the existing stepper.initialize resets its count;
+    this candidate is deliberately NOT an installed runtime continuation.
+    """
+    topology: RemapTopology
+    g: tuple
+    base_count: int
+    level_steps: tuple
+    phase: str
+    origin: tuple
+    spacing: float
+
+
+def _same_parent(a, b):
+    return (type(a) is type(b) and a.bounds == b.bounds
+            and a.coarse_shape == b.coarse_shape
+            and a.tau_coarse == b.tau_coarse and a.tau_fine == b.tau_fine)
+
+
+def _read_leaves(topology, arrays, layouts):
+    """Unpack only owned cores; nonleaf values are not evidence of validity."""
+    if len(arrays) != 3 or len(layouts) != 3:
+        raise ValueError("expected three G arrays and three layouts")
+    result = []
+    for a, layout, leaf in zip(arrays, layouts, topology.leaf_masks):
+        if not isinstance(a, np.ndarray) or a.dtype not in (np.dtype('float32'), np.dtype('float64')):
+            raise ValueError("G arrays must be FP32 or FP64 NumPy arrays")
+        if layout is not None:
+            if layout.shape != leaf.shape or not layout.owned[leaf].all():
+                raise ValueError("layout must own every required leaf donor")
+            a = layout.to_dense(a)
+        if a.shape != (9, *leaf.shape):
+            raise ValueError("G must have Q-major (9, ny, nx) shape")
+        values = a[:, leaf]
+        if not np.isfinite(values).all() or (values.sum(axis=0, dtype=np.float64) <= 0).any():
+            raise ValueError("required leaf donors must be finite with positive density")
+        result.append(a)
+    if len({a.dtype for a in result}) != 1:
+        raise ValueError("all levels must use the same population dtype")
+    return result
+
+
+def remap_synchronized(old, new, old_g, *, base_count, level_steps, phase,
+                       origin=(0., 0.), spacing=1., layouts=None):
+    """Return a non-mutating candidate at the SAME physical time and frame.
+
+    ``old``/``new`` are RemapTopology.from_config results. L0/L1 geometry and
+    relaxation parameters must match; only the L2 rectangle set may change.
+    ``old_g`` contains current owned-leaf G, dense or (with ``layouts``) blocked.
+    ``level_steps`` must equal (base_count, 2*base_count, 4*base_count), and
+    ``phase`` must be 'postcollision_g'. These stamps are caller assertions,
+    not proof of device quiescence or active-buffer freshness.
+
+    Preserve same-level leaf overlap bitwise. Copy all nine parent G values
+    to newly refined children; average ALL four old children (including q0)
+    for newly exposed L1. Restriction accumulates in FP64, then rounds once to
+    the source dtype. No equilibrium reconstruction or nonequilibrium scaling.
+    All destination nonleaves are NaN. Inputs and their write flags are intact.
+    """
+    if (not isinstance(old, RemapTopology) or not isinstance(new, RemapTopology)
+            or not _same_parent(old.parent, new.parent)
+            or old.child.tau_fine != new.child.tau_fine):
+        raise ValueError("remap requires unchanged L0/L1 geometry and relaxation")
+    if (old.pinned_rectangles != new.pinned_rectangles
+            or 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 body rectangles and physical boundary")
+    if not isinstance(level_steps, (tuple, list)):
+        raise ValueError("level_steps must be three completed integer counts")
+    counts = (base_count, *level_steps)
+    if (len(counts) != 4 or any(isinstance(n, (bool, np.bool_))
+                              or not isinstance(n, Integral) or n < 0 for n in counts)
+            or tuple(level_steps) != (base_count, 2*base_count, 4*base_count)
+            or phase != 'postcollision_g'):
+        raise ValueError("requires synchronized completed counts and postcollision_g")
+    try:
+        frame = np.asarray(origin, dtype=np.float64)
+        spacing = float(spacing)
+    except (TypeError, ValueError):
+        raise ValueError("frame needs a finite (x,y) origin and positive spacing") from None
+    if (frame.shape != (2,) or not np.isfinite(frame).all()
+            or not np.isscalar(spacing) or not np.isfinite(spacing) or spacing <= 0):
+        raise ValueError("frame needs a finite (x,y) origin and positive spacing")
+    source = _read_leaves(old, old_g, (None,)*3 if layouts is None else layouts)
+    candidate = [np.full((9, *mask.shape), np.nan, dtype=source[0].dtype)
+                 for mask in new.leaf_masks]
+    for level, (before, after) in enumerate(zip(old.leaf_masks, new.leaf_masks)):
+        overlap = before & after
+        candidate[level][:, overlap] = source[level][:, overlap]
+
+    # A child set uses a level-wide grid in its parent's allocated coordinates.
+    # Thus child (y,x) belongs to parent (y//2,x//2), even for a shifted L1 window.
+    y, x = np.nonzero(new.leaf_masks[2] & ~old.leaf_masks[2])
+    if not old.leaf_masks[1][y//2, x//2].all():
+        raise ValueError("new L2 requires an old owned L1 leaf")
+    candidate[2][:, y, x] = source[1][:, y//2, x//2]
+
+    y, x = np.nonzero(new.leaf_masks[1] & ~old.leaf_masks[1])
+    children = [(2*y+dy, 2*x+dx) for dy in (0, 1) for dx in (0, 1)]
+    if any(not old.leaf_masks[2][cy, cx].all() for cy, cx in children):
+        raise ValueError("reactivated L1 requires all four old owned L2 leaves")
+    # Scale first to avoid an overflowing intermediate sum for finite FP64 G.
+    mean = sum(source[2][:, cy, cx].astype(np.float64) * .25 for cy, cx in children)
+    candidate[1][:, y, x] = mean
+    _read_leaves(new, candidate, (None,)*3)
+    return RemapCandidate(new, tuple(candidate), int(base_count), tuple(map(int, level_steps)),
+                          phase, tuple(map(float, frame)), float(spacing))

tests/__init__.py

@@ -0,0 +1 @@
+"""CelerisLab test package; keeps shared test helpers independent of discovery paths."""

tests/unit/test_hierarchy.py

@@ -0,0 +1,97 @@
+"""CPU geometry tests for lbm.hierarchy (no CUDA). Geometry construction only, not LBM numerics."""
+import unittest
+
+import numpy as np
+from CelerisLab.config import LBMConfig
+from CelerisLab.lbm import refinement as r
+from CelerisLab.lbm.hierarchy import build_hierarchy, leaf_partition
+
+CFG = LBMConfig(nx=48, ny=48, viscosity=.02)
+L1 = [(6, 6, 26, 26), (32, 30, 42, 42)]   # L0 cells; L1 ACTIVE x/y 12..52 and x 64..84, y 60..84
+L2 = [(20, 20, 44, 44), (72, 68, 76, 76)]  # L1 cells; L2 ACTIVE 40..88 and x 144..152, y 136..152
+L3 = [(56, 56, 72, 72)]                    # L2 cells
+FOUR = [L1, L2, L3]
+SMALL = LBMConfig(nx=16, ny=16, viscosity=.1)
+
+
+class FourLevels(unittest.TestCase):
+    @classmethod
+    def setUpClass(cls):
+        cls.h = build_hierarchy(CFG, FOUR)
+
+    def test_config_chain(self):
+        self.assertEqual(self.h.depth, 4)
+        for level, c in enumerate(self.h.configs):
+            self.assertEqual((c.ny, c.nx), (48 << level, 48 << level))
+            self.assertAlmostEqual(c.viscosity, .02 * 2 ** level, places=12)
+            self.assertEqual(self.h.kinds[level].shape, (c.ny, c.nx))
+        for level, p in enumerate(self.h.pairs):
+            self.assertAlmostEqual(p.tau_coarse, 1 / self.h.configs[level].omega, places=12)
+            self.assertAlmostEqual(p.tau_fine, 1 / self.h.configs[level + 1].omega, places=12)
+
+    def test_partition_and_composition(self):
+        h = self.h
+        np.testing.assert_array_equal(leaf_partition(h), np.ones((48, 48)))
+        for level in (1, 2):
+            covered, raw = h.pairs[level].covered, h.pairs[level - 1].fine_kind
+            self.assertTrue((h.kinds[level][covered] == r.OUTER_GHOST).all())
+            np.testing.assert_array_equal(h.kinds[level][~covered], raw[~covered])
+        np.testing.assert_array_equal(h.kinds[3], h.pairs[2].fine_kind)
+
+    def test_depth_is_list_length(self):
+        for depth in range(1, 5):
+            h = build_hierarchy(CFG, FOUR[:depth - 1])
+            self.assertEqual(h.depth, depth)
+            np.testing.assert_array_equal(leaf_partition(h), np.ones((48, 48)))
+            for level in range(depth - 1):
+                np.testing.assert_array_equal(h.kinds[level], self.h.kinds[level])
+
+    def test_illegal_coverage_reports_first_guard(self):
+        cases = [([L1, [(10, 20, 30, 44)]], "inherited ACTIVE"),
+                 ([L1, L2[:1] + [(54, 54, 58, 58)]], "inherited ACTIVE"),
+                 ([L1, L2, [(40, 40, 60, 60)]], "donors/receivers"),
+                 ([L1, L2, [(41, 41, 60, 60)]], "donors/receivers"),
+                 ([L1, L2, [(56, 56, 64, 64), (66, 56, 72, 64)]], "closer than")]
+        for levels, reason in cases:
+            with self.subTest(levels=levels), self.assertRaisesRegex(ValueError, reason):
+                build_hierarchy(CFG, levels)
+        build_hierarchy(CFG, [L1, L2, [(42, 42, 60, 60)]])  # minimal legal margin
+
+
+class LegacyForms(unittest.TestCase):
+    def test_scalar_parent_and_child_match_existing_calls(self):
+        local = r.FixedPatch2D.from_config(SMALL, (4, 4, 12, 12))
+        child, kind = r.nest_patch(SMALL, local, (8, 8, 16, 16))
+        h = build_hierarchy(SMALL, [(4, 4, 12, 12), (8, 8, 16, 16)])
+        np.testing.assert_array_equal(h.kinds[0], np.where(local.covered, r.OUTER_GHOST, r.ACTIVE))
+        np.testing.assert_array_equal(h.kinds[1], kind)
+        np.testing.assert_array_equal(h.kinds[2], child.fine_kind)
+        for name in ("ghost_fine", "ghost_coarse", "coalesce_coarse", "coalesce_fine"):
+            np.testing.assert_array_equal(getattr(h.pairs[1], name), getattr(child, name))
+        np.testing.assert_array_equal(leaf_partition(h), np.ones((16, 16)))
+
+    def test_set_child_under_scalar_parent(self):
+        local = r.FixedPatch2D.from_config(SMALL, (4, 4, 12, 12))
+        child, kind = r.nest_patch_set(SMALL, local, [(8, 8, 16, 16)])
+        h = build_hierarchy(SMALL, [(4, 4, 12, 12), [(8, 8, 16, 16)]])
+        np.testing.assert_array_equal(h.kinds[1], kind)
+        np.testing.assert_array_equal(h.pairs[1].coalesce_fine, child.coalesce_fine)
+        np.testing.assert_array_equal(leaf_partition(h), np.ones((16, 16)))
+
+    def test_scalar_child_under_set_parent_is_rejected(self):
+        with self.assertRaisesRegex(ValueError, "scalar parent"):
+            build_hierarchy(SMALL, [[(4, 4, 12, 12)], (12, 12, 20, 20)])
+
+    def test_wall_rows_and_single_level(self):
+        h = build_hierarchy(SMALL, [(4, 4, 12, 12)], wall_aware=True)
+        self.assertTrue((h.kinds[0][[0, -1]] == r.OUTER_GHOST).all())
+        one = build_hierarchy(SMALL, [])
+        self.assertEqual(one.depth, 1)
+        self.assertTrue((one.kinds[0] == r.ACTIVE).all())
+        walled = build_hierarchy(SMALL, [], wall_aware=True)  # regression: no refinement
+        self.assertTrue((walled.kinds[0][[0, -1]] == r.OUTER_GHOST).all())
+        self.assertTrue((walled.kinds[0][1:-1] == r.ACTIVE).all())
+
+
+if __name__ == "__main__":
+    unittest.main()

tests/unit/test_jet_condition.py

@@ -0,0 +1,32 @@
+"""Jet boundary-condition hook: velocity implemented, pressure reserved (host only, no CUDA)."""
+import pytest
+
+from CelerisLab.body.jets import Jet, Parallel, Tube
+
+
+def region():
+    return Tube((16., 12.), 0., 3., .7)
+
+
+def test_default_is_velocity_and_survives_amplitude_change():
+    jet = Jet(region(), Parallel(0.), amplitude=.01)
+    assert jet.condition == "velocity"
+    assert jet.with_amplitude(-.02).condition == "velocity"
+    # positional calls keep their meaning; condition is the new last field
+    assert (Jet(region(), Parallel(0.), "parabolic", .01, "velocity")
+            == Jet(region(), Parallel(0.), "parabolic", .01))
+
+
+def test_pressure_is_rejected_before_amplitude_validation():
+    # An amplitude that the velocity check would reject must not decide the error.
+    for amplitude in (0., 5., float("nan")):
+        with pytest.raises(NotImplementedError, match="pressure jets are reserved"):
+            Jet(region(), Parallel(0.), amplitude=amplitude, condition="pressure")
+    with pytest.raises(NotImplementedError):
+        Jet(region(), None, profile="no-such-profile", condition="pressure")
+
+
+def test_unknown_condition_is_a_value_error():
+    for condition in ("Velocity", "flow_rate", "", None):
+        with pytest.raises(ValueError, match="condition"):
+            Jet(region(), Parallel(0.), amplitude=.01, condition=condition)

tests/unit/test_multiple_pinned_remap.py

@@ -0,0 +1,153 @@
+"""Multiple body pins and variable fluid sets; mocked runtime is host-only."""
+import unittest
+from unittest.mock import patch
+
+import numpy as np
+
+from CelerisLab.config import LBMConfig
+from CelerisLab.body.jets import Jet, Tube, Parallel
+from CelerisLab.body.objects import SurfaceActuation
+from CelerisLab.lbm.refinement import FixedPatch2D
+from CelerisLab.lbm.synchronized_remap import RemapTopology
+from tests.unit.test_pinned_body_remap import body, mesh
+from tests.unit.test_refinement_three_level_multipatch import mock_backend
+from tests.unit.test_synchronized_remap import fields, integrals, remap
+
+PINS = ((16, 16, 32, 24), (56, 16, 72, 24))
+FLUID = ((16, 32, 32, 40), (56, 32, 72, 40))
+ADDED = ((16, 44, 32, 48),)
+MOVED = ((20, 32, 36, 40), (60, 32, 76, 40))
+
+
+def objects():
+    return (body(center=(12., 10.)), body(center=(32., 10.)))
+
+
+def ready():
+    bodies = objects()
+    jets = tuple(Jet(Tube((x-1.5, 10.), 0., 3., .3), Parallel(0.),
+                     'parabolic', .003) for x in (12., 32.))
+    with patch('tests.unit.test_refinement_three_level_multipatch.bodies_and_jets',
+               return_value=(bodies, jets)):
+        run, module, allocations = mock_backend(children=PINS+FLUID, bodies=True)
+    storage = {}
+    module.cuda.memcpy_htod.side_effect = lambda p, a: storage.__setitem__(id(p), np.array(a, copy=True))
+    module.cuda.memcpy_dtoh.side_effect = lambda a, p: np.copyto(a, storage[id(p)])
+    run.initialize(*fields(mesh(PINS+FLUID, bodies, jets), np.float32))
+    run._step_count = 7
+    return run, module, allocations, bodies
+
+
+class MultiplePinnedRemap(unittest.TestCase):
+    def check_transfer(self, old, new, arrays, tolerance):
+        out = remap(old, new, arrays)
+        before = integrals(old, arrays)
+        np.testing.assert_allclose(integrals(new, out.g), before, rtol=0,
+                                   atol=tolerance*np.abs(before[:9]).sum())
+        for a, b, om, nm in zip(arrays, out.g, old.leaf_masks, new.leaf_masks):
+            self.assertEqual(a[:, om & nm].tobytes(), b[:, om & nm].tobytes())
+            self.assertTrue(np.isnan(b[:, ~nm]).all())
+        return out.g
+
+    def test_add_move_remove_and_body_only_all_populations(self):
+        for dtype, tolerance in ((np.float64, 2e-13), (np.float32, 2e-7)):
+            old = mesh(PINS+FLUID, objects())
+            arrays = fields(old, dtype)
+            for rects in (PINS+FLUID+ADDED, PINS+MOVED+ADDED, PINS, PINS+FLUID):
+                new = mesh(rects, objects())
+                self.assertEqual(new.pinned_rectangles, PINS)
+                arrays = self.check_transfer(old, new, arrays, tolerance)
+                old = new
+
+    def test_shared_pin_order_invariance_and_body_free(self):
+        bodies = (body(center=(11., 14.)), body(center=(14., 14.)))
+        rects = ((16, 16, 32, 40), (56, 16, 72, 40))
+        old, new = mesh(rects, bodies), mesh(rects[::-1], bodies)
+        self.assertEqual(old.pinned_rectangles, (rects[0],))
+        self.assertEqual(new.pinned_rectangles, old.pinned_rectangles)
+        self.check_transfer(old, new, fields(old), 2e-13)
+        only = mesh((rects[0],), bodies)
+        self.check_transfer(old, only, fields(old), 2e-13)
+        self.assertEqual(mesh(rects, ()).pinned_rectangles, ())
+
+    def test_shifted_scalar_parent_coordinates(self):
+        cfg = LBMConfig(nx=48, ny=32, viscosity=.1)
+        parent = FixedPatch2D.from_config(cfg, (4, 4, 44, 28))
+        self.assertEqual(parent.fine_origin, (2., 2.))
+        bodies = (body(center=(14., 12.)), body(center=(34., 12.)))
+        def topology(rects):
+            return RemapTopology.from_config(cfg, parent, rects, child_objects=bodies)
+        old, new = topology(PINS+FLUID), topology(PINS+MOVED)
+        self.assertEqual(old.pinned_rectangles, PINS)
+        self.check_transfer(old, new, fields(old), 2e-13)
+
+    def test_direct_helper_rejects_each_changed_pin_and_body_order(self):
+        old = mesh(PINS+FLUID, objects())
+        candidates = [mesh(PINS+FLUID, objects()[::-1])]
+        for index in range(2):
+            changed = list(PINS+FLUID)
+            x0, y0, x1, y1 = changed[index]
+            changed[index] = (x0-4, y0, x1, y1)
+            candidates.append(mesh(changed, objects()))
+        for new in candidates:
+            # Invalid donors ensure boundary validation precedes array reads.
+            with self.assertRaisesRegex(ValueError, 'pinned'):
+                remap(old, new, None)
+
+    def test_live_independent_jets_and_source_ownership(self):
+        run, module, allocations, bodies = ready()
+        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)
+        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'],
+                                          run._boundaries[2].tables['surface_velocity'])
+            clone.set_jet_amplitude(1, .005, level=2)
+            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):
+                np.testing.assert_array_equal(got, want)
+        finally:
+            clone.close()
+            run.close()
+
+    def test_sub_ulp_cut_table_geometry_change_rejected(self):
+        run, module, allocations, bodies = ready()
+        changed = objects()
+        changed[1].geometry._radius = np.nextafter(.8, np.inf)
+        start = len(allocations)
+        module.cuda.memcpy_dtoh.reset_mock()
+        with self.assertRaisesRegex(ValueError, 'boundary'):
+            run.remap_child(PINS+FLUID, child_objects=changed)
+        self.assertEqual(len(allocations), start)
+        module.cuda.memcpy_dtoh.assert_not_called()
+        with self.assertRaisesRegex(ValueError, 'pinned'):
+            remap(mesh(PINS+FLUID, bodies), mesh(PINS+FLUID, changed), None)
+        run.close()
+
+    def test_each_pin_and_body_identity_rejected_before_device_access(self):
+        run, module, allocations, bodies = ready()
+        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)]
+        for index in range(2):
+            changed = list(PINS+FLUID)
+            x0, y0, x1, y1 = changed[index]
+            changed[index] = (x0-4, y0, x1, y1)
+            cases.append((changed, bodies))
+            cases.append((tuple(r for i, r in enumerate(PINS+FLUID) if i != index), bodies))
+        start = len(allocations)
+        module.cuda.memcpy_dtoh.reset_mock()
+        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)
+        self.assertEqual(len(allocations), start)
+        module.cuda.memcpy_dtoh.assert_not_called()
+        run._stream.synchronize.assert_not_called()
+        self.assertEqual((run._status, run.step_count), ('ready', 7))
+        run.close()

tests/unit/test_pinned_body_remap.py

@@ -0,0 +1,130 @@
+"""Bounded pinned-body transfer and host lifecycle; no CUDA numerical claim."""
+import unittest
+from unittest.mock import patch
+import numpy as np
+
+from CelerisLab.body.geometry.circle import CircleGeometry
+from CelerisLab.body.objects import SimObject, SurfaceActuation
+from CelerisLab.body.jets import Jet, Tube, Parallel
+from CelerisLab.lbm.synchronized_remap import RemapTopology, same_physical_boundary
+from tests.unit.test_synchronized_remap import CFG, PARENTS, fields, remap
+from tests.unit.test_refinement_three_level_multipatch import mock_backend
+from CelerisLab.lbm.refinement import FixedPatchSet2D
+
+A = ((16,16,32,40), (56,16,72,40))
+B = ((16,16,32,40), (60,16,76,40))
+
+
+def body(spec=None, center=(12.,14.)):
+    return SimObject(71, geometry=CircleGeometry(*center,.8), center=center,
+                     radius=.8, surface_actuation=spec)
+
+
+def jet():
+    return Jet(Tube((10.5,14.),0.,3.,.3), Parallel(0.), 'parabolic', .003)
+
+
+def mesh(rects, objects, jets=()):
+    return RemapTopology.from_config(CFG, FixedPatchSet2D.from_config(CFG,PARENTS),
+                                    rects, child_objects=objects, child_jets=jets)
+
+
+def ready(objects=None, jets=None):
+    objects = (body(),) if objects is None else objects
+    jets = (jet(),) if jets is None else jets
+    with patch('tests.unit.test_refinement_three_level_multipatch.bodies_and_jets',
+               return_value=(objects,jets)):
+        run, module, allocations = mock_backend(children=A,bodies=True)
+    storage = {}
+    module.cuda.memcpy_htod.side_effect = lambda p,a: storage.__setitem__(id(p),np.array(a,copy=True))
+    module.cuda.memcpy_dtoh.side_effect = lambda a,p: np.copyto(a,storage[id(p)])
+    run.initialize(*fields(mesh(A,objects,jets),np.float32))
+    run._step_count = 7
+    return run,module,allocations,storage,objects
+
+
+class PinnedBodyRemap(unittest.TestCase):
+    def test_fluid_only_transfer_poison_and_round_trip(self):
+        old,new = mesh(A,(body(),)),mesh(B,(body(),))
+        arrays=fields(old,np.float32)
+        for a,leaf in zip(arrays,old.leaf_masks): a[:,~leaf]=12345
+        out=remap(old,new,arrays)
+        for a,b,before,after in zip(arrays,out.g,old.leaf_masks,new.leaf_masks):
+            np.testing.assert_array_equal(a[:,before & after],b[:,before & after])
+            self.assertTrue(np.isnan(b[:,~after]).all())
+        self.assertTrue(same_physical_boundary(old.boundary,new.boundary))
+        back=remap(new,old,out.g)
+        self.assertTrue(np.isnan(back.g[2][:,~old.leaf_masks[2]]).all())
+        self.assertEqual(old.pinned_rectangles,(A[0],))
+
+    def test_current_jets_and_source_independence(self):
+        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)
+        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))
+        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)
+        clone.close();run.close()
+
+    def test_static_body_without_actuation(self):
+        run,module,allocations,storage,objects=ready(jets=())
+        clone=run.remap_child(A,child_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)
+        clone.close();run.close()
+
+    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.)),)),
+               (((20,16,36,40),B[1]),objects),(((16,20,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)
+        self.assertEqual(len(allocations),start)
+        module.cuda.memcpy_dtoh.assert_not_called()
+        self.assertEqual((run._status,run.step_count),('ready',7))
+        run.close()
+
+    def test_legacy_velocity_and_zero_velocity_ownership_rejected(self):
+        for speed in (.002,0.):
+            run,module,allocations,storage,objects=ready(
+                (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),))
+                module.cuda.memcpy_dtoh.assert_not_called()
+            clone=run.remap_child(B,child_objects=objects)
+            clone.close();run.close()
+
+    def test_allocation_failure_keeps_old_ready_and_count(self):
+        run,module,allocations,storage,objects=ready()
+        before=run.download(); saved=list(allocations)
+        allocate=module.cuda.mem_alloc.side_effect
+        calls=0
+        def fail(size):
+            nonlocal calls
+            calls+=1
+            if calls == 3: raise RuntimeError('candidate allocation')
+            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)
+        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()
+        for a in allocations[len(saved):]: a.free.assert_called_once()
+        run.close()
+
+    def test_direct_helper_rejects_changed_pinned_rectangle(self):
+        old=mesh(A,(body(),))
+        new=mesh(((12,16,32,40),A[1]),(body(),))
+        with self.assertRaisesRegex(ValueError,'pinned'): remap(old,new,fields(old))

tests/unit/test_raw_exchange.py

@@ -0,0 +1,170 @@
+"""CPU reduction and real host scheduling with CUDA doubles, not GPU physics."""
+from contextlib import nullcontext
+from dataclasses import replace
+from unittest.mock import Mock
+import numpy as np
+import unittest
+from CelerisLab.body.geometry.circle import CircleGeometry
+from CelerisLab.body.objects import SimObject, SurfaceActuation
+from CelerisLab.config import LBMConfig
+from CelerisLab.lbm.blocks import BlockLayout2D
+from CelerisLab.lbm.refinement import FixedPatch2D, prepare_fine_boundary
+from CelerisLab.lbm.refinement_tables import _prepare_blocked_boundary
+from CelerisLab.lbm.raw_exchange import summarize_exchange
+from tests.unit.test_refinement_three_level_multipatch import host_module
+
+
+def body(center=(16., 16.), radius=.5, speed=0.):
+    return SimObject(91, geometry=CircleGeometry(*center, radius), center=center, radius=radius,
+                     surface_actuation=SurfaceActuation(.7, .85, speed))
+
+
+def backend(observe=True, obj=None):
+    host = host_module()
+    calls = []
+    module = Mock()
+    module.get_global.return_value = (Mock(), 20)
+    module.get_function.side_effect = lambda name: Mock(side_effect=lambda *a, **k: calls.append((name, a, k)))
+    host.compiler.load_module.return_value = module
+    host.cuda.mem_alloc.side_effect = lambda n: Mock(size=n)
+    host.cuda.memcpy_dtoh.side_effect = lambda array, pointer: array.fill(0)
+    ctx = Mock(activate=nullcontext, sm_arch='sm_70')
+    run = host.FixedPatchStepper(LBMConfig(nx=48, ny=48, viscosity=.02), context=ctx,
+        refinements=[[(6, 6, 26, 26)], [(20, 20, 44, 44)], [(56, 56, 72, 72)]],
+        level_blocks=[(4, 4)]*4, level_objects=[(), (), (), (body() if obj is None else obj,)],
+        observe_raw_exchange=observe)
+    arrays = [np.full(l.shape, 1/9, np.float32) for l in run._levels]
+    run.initialize(*arrays)
+    return run, host, calls, arrays
+
+
+class RawExchange(unittest.TestCase):
+    def _check_weights_split_and_window(self, level):
+        h = 2. ** -level
+        rows = np.array([[2, -3, 5, 7], [-1, 4, -2, -6], [99, 99, 99, 99]], np.float32)
+        kw = dict(body_id=0, level=level, h=h, time_end=3, reference_center=(7, 8))
+        last = summarize_exchange(rows, [0, 0, 1], [False, True, False], window='last_substep', **kw)
+        integral = rows.astype(float)*[h*h, h*h, h**3, h*h]/h
+        full = summarize_exchange(integral, [0, 0, 1], [False, True, False], window='last_base_step', **kw)
+        np.testing.assert_array_equal(last['total']['force'], np.array([1, 1])*h)
+        self.assertEqual(last['total']['torque'], 3 * h * h)
+        self.assertEqual(last['total']['mass_into_fluid'], h * h)
+        for name in ('force', 'torque', 'mass_rate'):
+            np.testing.assert_array_equal(last['total'][name], full['total'][name])
+        for name in ('force', 'impulse', 'torque', 'angular_impulse', 'mass_into_fluid'):
+            np.testing.assert_array_equal(last['total'][name],
+                                         np.asarray(last['wall'][name])+last['aperture'][name])
+        self.assertEqual(last['time_start'], 3 - h)
+        self.assertEqual(full['time_start'], 2)
+        self.assertEqual(last['sample_count'], 1)
+        self.assertEqual(full['sample_count'], 2 ** level)
+        self.assertEqual(full['reference_center'], (7, 8))
+
+    def test_final_substep_is_not_interval_mean(self):
+        kw = dict(body_id=0, level=1, h=.5, time_end=1, reference_center=(0, 0))
+        last = summarize_exchange([[4, 0, 0, 0]], [0], [False], window='last_substep', **kw)
+        full = summarize_exchange([[1.5, 0, 0, 0]], [0], [False], window='last_base_step', **kw)
+        self.assertEqual(last['total']['force'][0], 2)
+        self.assertEqual(full['total']['force'][0], 1.5)
+
+    def test_host_metadata_zero_jet_and_rows(self):
+        cfg = LBMConfig(nx=32, ny=32, viscosity=.1)
+        pair = FixedPatch2D.from_config(cfg, (4, 4, 28, 28))
+        boundary = prepare_fine_boundary(cfg, pair, [body(radius=2.15)])
+        local = replace(cfg, nx=pair.fine_shape[1], ny=pair.fine_shape[0], viscosity=.2)
+        local.validate()
+        layout = BlockLayout2D.build(*pair.fine_shape, 8, 8, periodic=(False, False))
+        b = _prepare_blocked_boundary(boundary, layout, local, ())
+        for key, idx in [('body_id', 3), ('rx', 4), ('ry', 5)]:
+            np.testing.assert_array_equal(b.tables[key], boundary.curved_host[idx])
+        self.assertTrue(b.tables['aperture'].any())
+        self.assertTrue((~b.tables['aperture']).any())
+        self.assertFalse(b.tables['surface_velocity'].any())
+        # Classification does not disappear when prescribed velocity is zero.
+        self.assertEqual(len(b.tables['cell']), len(boundary.curved_host[0]))
+
+    def test_four_level_schedule_publication_and_fixed_storage(self):
+        run, host, calls, arrays = backend()
+        b = run._level_boundaries[3]
+        rows = len(b.tables['cell'])
+        self.assertEqual(b.last_raw.size, rows * 16)
+        self.assertEqual(b.interval_integral.size, rows * 32)
+        with self.assertRaisesRegex(RuntimeError, 'completed'):
+            run.read_raw_exchange(0, level=3)
+        run.step()
+        names = [name for name, _, _ in calls]
+        self.assertEqual(names.count('BlockAccumulateExchange'), 8)
+        for i, (name, args, kw) in enumerate(calls):
+            if name == 'BlockAccumulateExchange':
+                self.assertEqual(names[i - 1], 'BlockCurvedBoundary')
+                self.assertEqual(names[i + 1], 'BlockStep')
+                self.assertIs(args[0], b.last_raw)
+                self.assertIs(args[1], b.interval_integral)
+                self.assertEqual(args[3:], (0.125, 0.125))
+                self.assertIs(kw['stream'], run._stream)
+        self.assertEqual(host.cuda.memset_d8_async.call_count, 1)
+        value = run.read_raw_exchange(0, level=3)
+        self.assertEqual(value['sample_count'], 8)
+        self.assertEqual(value['time_start'], 0)
+        run.step(0)
+        self.assertEqual(run.read_raw_exchange(0, level=3)['time_end'], 1)
+        self.assertEqual(run.read_raw_exchange(0, level=3, window='last_substep')['time_start'], 0.875)
+        run.step(2)
+        self.assertEqual(run.read_raw_exchange(0, level=3)['time_start'], 2)
+        run.initialize(*arrays)
+        with self.assertRaisesRegex(RuntimeError, 'completed'):
+            run.read_raw_exchange(0, level=3)
+        run.close()
+
+    def test_torque_reference_uses_geometry_center_not_optional_object_center(self):
+        for center in (None, (3., 4.)):
+            obj = SimObject(91, geometry=CircleGeometry(16., 16., .5), center=center, radius=.5)
+            run, _, _, _ = backend(obj=obj)
+            run.step()
+            result = run.read_raw_exchange(0, level=3)
+            self.assertEqual(result['reference_center'], (16.0, 16.0))
+            self.assertNotEqual(result['reference_center'], tuple(obj.center))
+            run.close()
+
+    def test_disabled_has_nulls_and_no_observation_allocation_or_launch(self):
+        run, host, calls, _ = backend(False)
+        b = run._level_boundaries[3]
+        self.assertIs(b.last_raw, b.interval_integral)
+        self.assertIs(b.interval_integral, None)
+        run.step()
+        for name, args, _ in calls:
+            if name == 'BlockCurvedBoundary':
+                self.assertEqual(tuple(args[8:]), (0, 0, 0))
+        self.assertTrue(all((n != 'BlockAccumulateExchange' for n, _, _ in calls)))
+        host.cuda.memset_d8_async.assert_not_called()
+        with self.assertRaisesRegex(RuntimeError, 'disabled'):
+            run.read_raw_exchange(0, level=3)
+        run.close()
+
+    def test_failed_sync_never_publishes_partial_interval(self):
+        run, host, calls, _ = backend()
+        run.step()
+        run._stream.synchronize.side_effect = RuntimeError('device failure')
+        with self.assertRaisesRegex(RuntimeError, 'device failure'):
+            run.step()
+        self.assertFalse(run._exchange_valid)
+        with self.assertRaisesRegex(RuntimeError, 'failed'):
+            run.read_raw_exchange(0, level=3)
+        run._stream.synchronize.side_effect = None
+        run.close()
+
+    def test_dense_observation_rejected_before_device_access(self):
+        host = host_module()
+        ctx = Mock()
+        with self.assertRaisesRegex(ValueError, 'blocked body'):
+            host.FixedPatchStepper(LBMConfig(nx=32, ny=32, viscosity=.1), (4, 4, 28, 28),
+                                   context=ctx, fine_objects=[body(radius=2.15)], observe_raw_exchange=True)
+        ctx.activate.assert_not_called()
+
+    def test_all_depth_weights(self):
+        for level in (1, 2, 3, 7):
+            self._check_weights_split_and_window(level)
+
+
+if __name__ == "__main__":
+    unittest.main()

tests/unit/test_refinement.py

@@ -228,6 +228,155 @@ class FixedPatchSetTests(unittest.TestCase):
             self.assertTrue(layout.owned.ravel()[s.ghost_fine.astype(np.intp)].all())
         with self.assertRaisesRegex(ValueError, "aligned"):
             self.build(self.FACE).fine_layout((5, 4))
+        for block in ((0, 4), (4, 0), (-4, 4), (4.0, 4), (True, 4), (4,), (4, 4, 4), (np.bool_(1), 4)):
+            with self.subTest(block=block), self.assertRaisesRegex(ValueError, "positive integers"):
+                self.build(self.FACE).fine_layout(block)
+
+
+class PatchSetBoundaryTests(unittest.TestCase):
+    """fine_objects/jets on a level-wide patch set (16x32, patches 4 cells apart)."""
+    PATCHES = ((2, 2, 14, 14), (18, 2, 30, 14))
+
+    @staticmethod
+    def body(x, y, obj_id=0):
+        from CelerisLab.body.geometry.circle import CircleGeometry
+        from CelerisLab.body.objects import SimObject
+        return SimObject(obj_id, geometry=CircleGeometry(x, y, 2.15), center=(x, y), radius=2.15)
+
+    def prepare(self, nx, rects, objects, jets=()):
+        cfg = LBMConfig(nx=nx, ny=16, viscosity=.1)
+        s = refinement.FixedPatchSet2D.from_config(cfg, rects)
+        return s, refinement.prepare_fine_boundary(cfg, s, objects, jets=jets)
+
+    def test_two_bodies_are_the_one_body_case_translated(self):
+        _, one = self.prepare(16, self.PATCHES[:1], [self.body(8., 8.)])
+        _, two = self.prepare(32, self.PATCHES, [self.body(8., 8., 7), self.body(24., 8., 7)])
+        f1, f2 = one.flags.reshape(32, 32), two.flags.reshape(32, 64)
+        np.testing.assert_array_equal(f2[:, :32], f1)
+        np.testing.assert_array_equal(f2[:, 32:], f1)
+        n = len(one.curved_host[0])
+        self.assertEqual(len(two.curved_host[0]), 2 * n)
+        body_id = two.curved_host[3]
+        self.assertEqual(sorted(set(body_id.tolist())), [0, 1])  # positions, not obj_id 7
+        # Translation must reproduce identical link geometry (prerequisite of the
+        # bitwise GPU tiling gate): dirs, q and fallback equal, receivers shifted.
+        fidx, dirs, q, fallback = (two.curved_host[i] for i in (0, 1, 2, 7))
+        for column in (dirs, q, fallback):
+            np.testing.assert_array_equal(column[:n], column[n:])
+        y, x = np.divmod(fidx.astype(np.int64), 64)
+        np.testing.assert_array_equal(y[:n], y[n:])
+        np.testing.assert_array_equal(x[:n] + 32, x[n:])
+
+    def test_body_outside_a_single_patch_rejected(self):
+        for x in (16., 13.):  # in the gap / crossing the patch edge
+            with self.subTest(x=x), self.assertRaises(ValueError):
+                self.prepare(32, self.PATCHES, [self.body(x, 8.)])
+
+    def test_one_jet_over_two_bodies_rejected_even_at_zero_amplitude(self):
+        from CelerisLab.body.jets import Jet, Parallel, Tube
+        bodies = [self.body(8., 8.), self.body(24., 8.)]
+        span = Jet(Tube((5., 8.), 0., 22., .6), Parallel(0.), amplitude=0.)
+        with self.assertRaisesRegex(ValueError, "one jet"):
+            self.prepare(32, self.PATCHES, bodies, jets=[span])
+        own = [Jet(Tube((5., 8.), 0., 6., .6), Parallel(0.), amplitude=.002),
+               Jet(Tube((21., 8.), 0., 6., .6), Parallel(0.), amplitude=.002)]
+        self.prepare(32, self.PATCHES, bodies, jets=own)  # one jet per body: accepted
+
+
+class NestPatchTests(unittest.TestCase):
+    """Three-level single chain: L2 bounds in L1 window cells (probe fixture)."""
+    def setUp(self):
+        self.cfg = config()
+        self.l1 = FixedPatch2D.from_config(self.cfg, (4, 4, 12, 12))  # 24x24 window
+
+    def test_probe_fixture_composes_kinds(self):
+        child, kind = refinement.nest_patch(self.cfg, self.l1, (8, 8, 16, 16))
+        self.assertEqual(child.fine_shape, (24, 24))
+        self.assertAlmostEqual(child.tau_coarse, self.l1.tau_fine)
+        base = self.l1.fine_kind
+        for role in (INNER_GHOST, OUTER_GHOST):
+            np.testing.assert_array_equal(kind[base == role], role)
+        self.assertTrue((kind[child.covered] == OUTER_GHOST).all())
+        self.assertEqual(int((kind == ACTIVE).sum()), 192)
+        self.assertEqual(int((kind == INNER_GHOST).sum()), 144)
+        self.assertEqual(int(child.covered.sum()), 64)
+
+    def test_invalid_nesting_rejected(self):
+        for bounds in ((2, 8, 16, 16),   # covers inherited INNER
+                       (4, 8, 16, 16),   # explosion donors in INNER
+                       (8, 8, 20, 16),   # donors reach OUTER/INNER at the far side
+                       (5, 8, 16, 16)):  # one-cell outward: restriction footprint/donors
+            with self.subTest(bounds=bounds), self.assertRaises(ValueError):
+                refinement.nest_patch(self.cfg, self.l1, bounds)
+
+
+class ChildBoundaryTests(unittest.TestCase):
+    @staticmethod
+    def body(center=(8., 8.), radius=.8, obj_id=99, spec=None):
+        from CelerisLab.body.geometry.circle import CircleGeometry
+        from CelerisLab.body.objects import SimObject
+        return SimObject(obj_id, geometry=CircleGeometry(*center, radius), center=center,
+                         radius=radius, surface_actuation=spec)
+
+    def prepare(self, objects, jets=(), large=False):
+        from dataclasses import replace
+        cfg = LBMConfig(nx=32 if large else 16, ny=32 if large else 16, viscosity=.1)
+        parent = FixedPatch2D.from_config(cfg, (4, 4, 28, 28) if large else (4, 4, 12, 12))
+        child, _ = refinement.nest_patch(cfg, parent, (8, 8, 40, 40) if large else (8, 8, 16, 16))
+        middle = replace(cfg, nx=parent.fine_shape[1], ny=parent.fine_shape[0], viscosity=.2)
+        return child, refinement.prepare_fine_boundary(middle, child, objects, jets=jets,
+                         parent_origin=parent.fine_origin, parent_spacing=.5)
+
+    def test_root_frame_matches_independent_parent_frame_projection(self):
+        from CelerisLab.body.jets import Jet, Parallel, Tube, jet_link_shapes
+        jet = Jet(Tube((6.5, 8.), 0., 3., .3), Parallel(0.), amplitude=.002)
+        original = self.body()
+        child, b = self.prepare([original], [jet])
+        middle = LBMConfig(nx=24, ny=24, viscosity=.2)
+        transformed = self.body((12., 12.), 1.6)
+        local_jet = Jet(Tube((9., 12.), 0., 6., .6), Parallel(0.), amplitude=.002)
+        other = refinement.prepare_fine_boundary(middle, child, [transformed], jets=[local_jet])
+        for a, c in zip(b.curved_host, other.curved_host):
+            np.testing.assert_array_equal(a, c)
+        np.testing.assert_array_equal(b.flags, other.flags)
+        np.testing.assert_allclose(b.curved_hits, 2 + .5 * other.curved_hits, atol=0, rtol=0)
+        for a, c in zip(jet_link_shapes([jet], b.curved_hits),
+                        jet_link_shapes([local_jet], other.curved_hits)):
+            np.testing.assert_allclose(a, c, atol=1e-14, rtol=0)
+        self.assertEqual(len(b.curved_host[0]), 60)
+        self.assertAlmostEqual(child.tau_fine, .5 + 12 * .1)
+        self.assertEqual(original.geometry.center, (8., 8.))
+        self.assertEqual(original.obj_id, 99)
+
+    def test_invalid_body_and_jet_contracts(self):
+        from CelerisLab.body.jets import Jet, Parallel, Tube
+        from CelerisLab.body.objects import SurfaceActuation
+        jet = Jet(Tube((6.5, 8.), 0., 3., .3), Parallel(0.), amplitude=.002)
+        for body in (self.body((6., 8.)), self.body((6.9, 8.)), self.body(radius=2.1)):
+            with self.assertRaises(ValueError):
+                self.prepare([body])
+        moving = self.body()
+        moving.state.vx = .01
+        with self.assertRaisesRegex(ValueError, "stationary"):
+            self.prepare([moving])
+        with self.assertRaises(ValueError):
+            self.prepare([self.body()], [jet, jet])
+        empty = Jet(Tube((0., 0.), 0., 1., .1), Parallel(0.), amplitude=0.)
+        with self.assertRaises(ValueError):
+            self.prepare([self.body()], [empty])
+        spec = SurfaceActuation(.7, .85, .002)
+        with self.assertRaisesRegex(ValueError, "SurfaceActuation"):
+            self.prepare([self.body(spec=spec)], [spec.as_jet((8., 8.), .8)])
+
+    def test_one_jet_two_real_bodies_same_caller_id_rejected(self):
+        from CelerisLab.body.jets import Jet, Parallel, Tube
+        bodies = [self.body((12., 14.), 1., 7), self.body((17., 14.), 1., 7)]
+        _, boundary = self.prepare(bodies, large=True)  # establish legal geometry first
+        self.assertEqual(np.unique(boundary.curved_host[3]).tolist(), [0, 1])
+        for amplitude in (0., .002):
+            jet = Jet(Tube((10., 14.), 0., 9., .4), Parallel(0.), amplitude=amplitude)
+            with self.assertRaisesRegex(ValueError, "one jet may drive only one body"):
+                self.prepare(bodies, [jet], large=True)
 
 
 if __name__ == "__main__":

tests/unit/test_refinement_three_level_multipatch.py

@@ -0,0 +1,334 @@
+"""CPU-only union geometry and production host-control-flow tests.
+
+The AST loader omits CUDA imports and supplies explicit mocks; it executes
+unchanged production definitions, never CUDA kernels or a numerical emulator.
+"""
+import ast
+from contextlib import nullcontext
+from dataclasses import replace
+from pathlib import Path
+import types
+import unittest
+from unittest.mock import Mock, patch
+
+import numpy as np
+from CelerisLab.config import LBMConfig
+from CelerisLab.lbm import refinement as r
+
+
+def host_module():
+    path=Path(r.__file__).with_name('refinement_stepper.py')
+    tree=ast.parse(path.read_text())
+    tree.body=[node for node in tree.body if not (
+        isinstance(node,ast.Import) and any(a.name=='pycuda.driver' for a in node.names)
+        or isinstance(node,ast.ImportFrom) and node.module in (
+            'cuda','block_stepper','curved_links'))]
+    module=types.ModuleType('CelerisLab.lbm._test_host_refinement')
+    module.__package__='CelerisLab.lbm'
+    from CelerisLab.body.jets import surface_link_velocity
+    module.__dict__.update(cuda=Mock(),compiler=Mock(),MAX_GRID_DIM=65535,
+                           surface_link_velocity=surface_link_velocity)
+    exec(compile(tree,str(path),'exec'),module.__dict__)
+    return module
+
+
+host=host_module()
+CFG = LBMConfig(nx=48, ny=32, viscosity=.1)
+PARENTS = ((4,4,20,28), (24,4,44,28))  # exact gap=4 coarse cells
+CHILDREN = ((16,16,32,40), (56,16,80,40))  # parent-allocated L1 cells
+
+
+class Gate(unittest.TestCase):
+    def setUp(self):
+        self.parent = r.FixedPatchSet2D.from_config(CFG, PARENTS)
+
+    def test_one_element_equivalence(self):
+        cfg = LBMConfig(nx=16,ny=16,viscosity=.1)
+        local = r.FixedPatch2D.from_config(cfg,(4,4,12,12))
+        global_parent = r.FixedPatchSet2D.from_config(cfg,[(4,4,12,12)])
+        a, ak = r.nest_patch(cfg,local,(8,8,16,16))
+        b, bk = r.nest_patch_set(cfg,global_parent,[(12,12,20,20)])
+        local_set, local_kind = r.nest_patch_set(cfg,local,[(8,8,16,16)])
+        np.testing.assert_array_equal(ak,local_kind)
+        np.testing.assert_array_equal(a.fine_kind,local_set.patches[0].fine_kind)
+        np.testing.assert_array_equal(a.ghost_coarse,local_set.ghost_coarse)
+        np.testing.assert_array_equal(a.coalesce_coarse,local_set.coalesce_coarse)
+        np.testing.assert_array_equal(ak, bk[4:28,4:28])
+        # L2 node-index offsets in global L2 grid: local parent origin=2 L0.
+        ay,ax=np.nonzero(a.fine_kind == r.ACTIVE)
+        by,bx=np.nonzero(b.fine_kind == r.ACTIVE)
+        np.testing.assert_array_equal(np.stack([ay+20,ax+20],-1),np.stack([by,bx],-1))
+        np.testing.assert_array_equal(a.fine_kind,b.fine_kind[20:44,20:44])
+        # Compare complete per-region maps after undoing set addressing.
+        bp=b.patches[0]
+        np.testing.assert_array_equal(a.coalesce_fine,bp.coalesce_fine)
+        aq,ac=np.divmod(a.coalesce_coarse,24*24)
+        bq,bc=np.divmod(bp.coalesce_coarse,32*32)
+        np.testing.assert_array_equal(aq,bq)
+        np.testing.assert_array_equal(np.stack(np.divmod(ac,24),-1)+4,np.stack(np.divmod(bc,32),-1))
+        ay,ax=np.divmod(a.ghost_coarse,24)
+        by,bx=np.divmod(bp.ghost_coarse,32)
+        np.testing.assert_array_equal(np.stack([ay+4,ax+4],-1),np.stack([by,bx],-1))
+
+    def test_twin_tile_chains_and_order(self):
+        children=((16,16,32,40),(56,16,72,40))
+        child,kind=r.nest_patch_set(CFG,self.parent,children)
+        other,ok=r.nest_patch_set(CFG,r.FixedPatchSet2D.from_config(CFG,PARENTS[::-1]),children[::-1])
+        for name in ('fine_kind','ghost_fine','ghost_coarse','coalesce_coarse','coalesce_fine'):
+            np.testing.assert_array_equal(getattr(child,name),getattr(other,name))
+        np.testing.assert_array_equal(kind,ok)
+        for p in child.patches:
+            self.assertEqual(p.fine_shape,child.patches[0].fine_shape)
+        np.testing.assert_array_equal(child.patches[0].fine_kind,child.patches[1].fine_kind)
+        self.assertEqual(child.tau_coarse,self.parent.tau_fine)
+        self.assertAlmostEqual(child.tau_fine,1.7)
+
+    def test_asymmetric_no_child_region(self):
+        child,kind=r.nest_patch_set(CFG,self.parent,[CHILDREN[0]])
+        np.testing.assert_array_equal(kind[:,48:],self.parent.fine_kind[:,48:])
+        self.assertEqual(int(child.covered.sum()),16*24)
+
+    def test_two_children_one_parent(self):
+        child,kind=r.nest_patch_set(CFG,self.parent,[(16,16,32,24),(16,28,32,40)])
+        self.assertEqual(len(child.patches),2)
+        self.assertEqual(int(child.covered.sum()),16*20)
+        self.assertEqual(int((kind==r.ACTIVE).sum()),int((self.parent.fine_kind==r.ACTIVE).sum())-16*20)
+
+    def test_holes_donors_receivers_footprint_and_invalid_pull(self):
+        child,_=r.nest_patch_set(CFG,self.parent,CHILDREN)
+        cases=[('inherited ACTIVE', lambda: r.nest_patch_set(CFG,self.parent,[(40,16,48,40)])),
+               ('donors/receivers', lambda: r.nest_patch_set(CFG,self.parent,[(8,16,32,40)]))]
+        for text,call in cases:
+            with self.subTest(text=text),self.assertRaisesRegex(ValueError,text):call()
+        # Synthetic upstream dependency inside child coverage isolates the final guard.
+        upstream=self.parent.coalesce_fine.copy()
+        upstream[0,:]=20*self.parent.fine_shape[1]+20  # q=0, deep inherited ACTIVE
+        with self.assertRaisesRegex(ValueError,'restriction footprint'):
+            r._nested_kind(replace(self.parent,coalesce_fine=upstream),child)
+        hole=int(np.flatnonzero(~self.parent.fine_layout((4,4)).owned.ravel())[0])
+        donors=child.ghost_coarse.copy();donors[0]=hole
+        with self.assertRaisesRegex(ValueError,'donors/receivers'):
+            r._nested_kind(self.parent,replace(child,ghost_coarse=donors))
+        receivers=child.coalesce_coarse.copy();receivers[0]=hole
+        with self.assertRaisesRegex(ValueError,'donors/receivers'):
+            r._nested_kind(self.parent,replace(child,coalesce_coarse=receivers))
+        # Remove all supplying rows; covered-to-uncovered pulls must be rejected.
+        with self.assertRaisesRegex(ValueError,'surviving parent ACTIVE pull'):
+            r._nested_kind(self.parent,replace(child,coalesce_coarse=child.coalesce_coarse[:0],coalesce_fine=child.coalesce_fine[:0]))
+
+    def test_existing_l2_body_jet_preparation_on_union(self):
+        from CelerisLab.body.geometry.circle import CircleGeometry
+        from CelerisLab.body.objects import SimObject
+        from CelerisLab.body.jets import Jet, Tube, Parallel, jet_link_shapes
+        child,_=r.nest_patch_set(CFG,self.parent,CHILDREN)
+        centers=((12.,14.),(34.,14.))
+        objects=[SimObject(99,geometry=CircleGeometry(*c,.8),center=c,radius=.8) for c in centers]
+        jets=[Jet(Tube((c[0]-1.5,c[1]),0.,3.,.3),Parallel(0.),amplitude=.002) for c in centers]
+        middle=replace(CFG,nx=self.parent.fine_shape[1],ny=self.parent.fine_shape[0],viscosity=.2)
+        boundary=r.prepare_fine_boundary(middle,child,objects,jets=jets,
+                                         parent_origin=self.parent.fine_origin,parent_spacing=.5)
+        ids=boundary.curved_host[3]
+        np.testing.assert_array_equal(np.unique(ids),[0,1])
+        self.assertEqual(int((ids==0).sum()),int((ids==1).sum()))
+        np.testing.assert_allclose(boundary.curved_hits[ids==1]-boundary.curved_hits[ids==0],
+                                   np.broadcast_to([22.,0.],boundary.curved_hits[ids==0].shape),atol=1e-14,rtol=0)
+        _,jid=jet_link_shapes(jets,boundary.curved_hits)
+        for i in (0,1):
+            self.assertTrue((ids[jid==i]==i).all())
+            self.assertTrue((jid==i).any())
+
+    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)
+        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])
+        inherited=self.parent.fine_kind != r.ACTIVE
+        np.testing.assert_array_equal(kind[inherited],self.parent.fine_kind[inherited])
+        self.assertTrue((kind[child.covered] == r.OUTER_GHOST).all())
+        l2=child.fine_layout((4,4))
+        ghosts,children,ckind=host.blocked_fine_tables(child,l2,order)
+        self.assertTrue((ghosts<l2.hole).all())
+        self.assertTrue((children<np.prod(l2.blocked_shape)).all())
+        for p,lay in ((self.parent,layout),(child,l2)):
+            self.assertTrue(lay.owned.ravel()[p.ghost_fine].all())
+            self.assertTrue(lay.owned.ravel()[p.coalesce_fine % np.uint64(np.prod(p.fine_shape))].all())
+        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)
+        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)
+
+
+def bodies_and_jets(speed=.002):
+    from CelerisLab.body.geometry.circle import CircleGeometry
+    from CelerisLab.body.objects import SimObject
+    from CelerisLab.body.jets import Jet, Tube, Parallel
+    centers=((12.,14.),(34.,14.))
+    bodies=[SimObject(99,geometry=CircleGeometry(*c,.8),center=c,radius=.8) for c in centers]
+    jets=[Jet(Tube((c[0]-1.5,c[1]),0.,3.,.3),Parallel(0.),amplitude=speed) for c in centers]
+    return bodies,jets
+
+
+KW=dict(base_block=(4,4),fine_block=(4,4),child_block=(4,4),
+        child_rectangles_in_parent=CHILDREN)
+
+
+def mock_backend(*, children=CHILDREN, bodies=False, fail_compile=False, fail_upload=False):
+    module=host_module()
+    context=Mock(activate=nullcontext,sm_arch='sm_70')
+    allocations=[]
+    def allocate(size):
+        pointer=Mock(name=f'allocation_{len(allocations)}')
+        allocations.append(pointer)
+        return pointer
+    module.cuda.mem_alloc.side_effect=allocate
+    def load(ptx):
+        gpu=Mock()
+        gpu.get_function.side_effect=lambda name: Mock(name=name)
+        gpu.get_global.side_effect=lambda name: (Mock(),20 if name=='d_params' else 16)
+        return gpu
+    module.compiler.load_module.side_effect=load
+    if fail_compile:
+        module.compiler.compile_kernel.side_effect=['first','second',RuntimeError('L2 compile failed')]
+    if fail_upload:
+        def upload(*args):
+            if len(allocations)>20:
+                raise RuntimeError('table upload failed')
+        module.cuda.memcpy_htod.side_effect=upload
+    objects,jets=bodies_and_jets() if bodies else ((),())
+    try:
+        run=module.FixedPatchStepper(CFG,PARENTS,context=context,
+            child_objects=objects,child_jets=jets,
+            **dict(KW,child_rectangles_in_parent=children))
+    except RuntimeError:
+        for allocation in allocations:
+            allocation.free.assert_called_once_with()
+        raise
+    return run,module,allocations
+
+
+class ConstructorAndSchedule(unittest.TestCase):
+    def test_invalid_inputs_fail_before_context_or_allocations(self):
+        context=Mock()
+        cases=[dict(child=(8,8,16,16)),dict(child_rectangles_in_parent=[]),
+               dict(child_rectangles_in_parent=[(16,16,32)]),
+               dict(child_rectangles_in_parent=[(16.5,16,32,40)]),
+               dict(child_rectangles_in_parent=[(16,16,32,24),(16,27,32,40)]),
+               dict(child_rectangles_in_parent=[(40,16,48,40)]),
+               dict(child_rectangles_in_parent=[(8,16,32,40)]),
+               dict(child_block=None),dict(fine_block=None),dict(base_block=None),
+               dict(child_block=(5,4)),dict(fine_block=(5,4)),
+               dict(child_block=(4,)),dict(child_rectangles_in_parent=(16,16,32,40))]
+        for changes in cases:
+            with self.subTest(changes=changes),self.assertRaises((ValueError,TypeError)):
+                host.FixedPatchStepper(CFG,PARENTS,context=context,**dict(KW,**changes))
+        with self.assertRaisesRegex(ValueError,'closer'):
+            host.FixedPatchStepper(CFG,[(4,4,20,28),(23,4,44,28)],context=context,**KW)
+        objects,jets=bodies_and_jets()
+        with self.assertRaisesRegex(ValueError,'only on L2'):
+            host.FixedPatchStepper(CFG,PARENTS,context=context,fine_objects=objects,**KW)
+        from CelerisLab.body.jets import Jet,Tube,Parallel
+        spanning=Jet(Tube((10.5,14.),0.,25.,.3),Parallel(0.),amplitude=0.)
+        with self.assertRaisesRegex(ValueError,'one jet'):
+            host.FixedPatchStepper(CFG,PARENTS,context=context,child_objects=objects,
+                                  child_jets=[spanning],**KW)
+        context.activate.assert_not_called()
+        host.cuda.mem_alloc.assert_not_called()
+
+    def test_aggregate_schedule_independent_of_region_count(self):
+        for children,bodies in ((CHILDREN,True),(CHILDREN[:1],False),
+                                (((16,16,32,24),(16,28,32,40)),False)):
+            run,_,allocations=mock_backend(children=children,bodies=bodies)
+            self.assertEqual(len(run._levels),3)
+            c,f,g=run._levels
+            self.assertEqual(len({id(level.current) for level in run._levels} |
+                                 {id(level.spare) for level in run._levels}),6)
+            events=[]
+            initial_c=c.current
+            def record(name):
+                return lambda source,*a,**kw: events.append((name,source))
+            run._pairs[0].explode.side_effect=record('explode01')
+            run._pairs[1].explode.side_effect=record('explode12')
+            for i,level in enumerate(run._levels):
+                level.kernel.side_effect=record(f'step{i}')
+            if bodies:
+                g.boundary.kernel.side_effect=record('bfl2')
+            run._pairs[1].coalesce.side_effect=lambda src,dst,*a,**kw: events.append(('restrict21',dst))
+            run._pairs[0].coalesce.side_effect=lambda src,dst,*a,**kw: events.append(('restrict10',dst))
+            run._status='ready'
+            expected=[('explode01',initial_c)]
+            for middle in (f.current,f.spare):
+                expected.append(('explode12',middle))
+                for leaf in (g.current,g.spare):
+                    expected.append(('halo2',leaf))
+                    if bodies:
+                        expected.append(('bfl2',leaf))
+                    expected.append(('step2',leaf))
+                expected.extend([('halo1',middle),('restrict21',middle),('step1',middle)])
+            expected.extend([('halo0',initial_c),('restrict10',initial_c),('step0',initial_c)])
+            def gather(level,buffer,**kw):
+                events.append((f'halo{next(i for i,l in enumerate(run._levels) if l is level)}',buffer))
+            with patch.object(run,'_gather_halo',side_effect=gather):
+                run.step()
+            self.assertEqual(events,expected)
+            self.assertEqual([l.kernel.call_count for l in run._levels],[1,2,4])
+            # A stale-feed implementation must fail the source-identity guard.
+            feeds=[]
+            advance=run._advance
+            def checked_advance(index, buffer, **kw):
+                if index == 2:  # L1 -> L2 feed: must be the L1 hop source
+                    self.assertIs(buffer,(f.current,f.spare)[len(feeds)])
+                    feeds.append(buffer)
+                advance(index,buffer,**kw)
+            with patch.object(run,'_advance',side_effect=lambda index,buffer,**kw:
+                              checked_advance(index,f.current if index == 2 else buffer,**kw)):
+                with self.assertRaises(AssertionError):
+                    run.step()
+            # A per-region shared-parent advance must fail the 1:2:4 guard.
+            run._status='ready'
+            counts=[0,0,0]
+            def count(index, multiple=1):
+                def step(*a,**kw):
+                    counts[index]+=multiple
+                return step
+            c.kernel.side_effect=count(0)
+            f.kernel.side_effect=count(1,2)  # mutant: twice per shared L1 hop
+            g.kernel.side_effect=count(2)
+            run.step()
+            with self.assertRaises(AssertionError):
+                self.assertEqual(counts,[1,2,4])
+            run.close();run.close()
+            for allocation in allocations:
+                allocation.free.assert_called_once_with()
+
+    def test_partial_construction_frees_every_allocation(self):
+        for options in (dict(fail_compile=True),dict(fail_upload=True)):
+            with self.subTest(options=options),self.assertRaisesRegex(RuntimeError,'failed'):
+                mock_backend(**options)
+
+
+
+
+class InitialLevels(unittest.TestCase):
+    """initialize() entry adapter: positional per-level arrays or the legacy names, never both."""
+
+    def test_forms(self):
+        a, b, c = object(), object(), object()
+        f = host._initial_levels
+        self.assertEqual(f((a, b), None, None, None), (a, b))
+        self.assertEqual(f((a, b, c), None, None, None), (a, b, c))
+        self.assertEqual(f((a, b), None, None, c), (a, b, c))      # old initialize(c, f, child_g=g)
+        self.assertEqual(f((), a, b, None), (a, b))
+        self.assertEqual(f((), a, b, c), (a, b, c))
+        for args in (((a,), None, None, c), ((a, b, c), None, None, c), ((a,), a, b, None),
+                     ((), a, None, None), ((), None, b, c)):
+            with self.subTest(args=args), self.assertRaises(TypeError):
+                f(*args)
+
+
+if __name__=='__main__':
+    unittest.main()

tests/unit/test_runtime_child_remap.py

@@ -0,0 +1,123 @@
+"""CPU mocks execute production remap lifecycle; no GPU numerical claim."""
+from dataclasses import replace
+import unittest
+from unittest.mock import patch
+
+import numpy as np
+
+from tests.unit.test_refinement_three_level_multipatch import mock_backend
+from tests.unit.test_synchronized_remap import A, B, fields, topology, remap
+
+
+def ready_backend():
+    run, module, allocations = mock_backend(children=A)
+    storage = {}
+    module.cuda.memcpy_htod.side_effect = lambda ptr, a: storage.__setitem__(id(ptr), np.array(a, copy=True))
+    module.cuda.memcpy_dtoh.side_effect = lambda a, ptr: np.copyto(a, storage[id(ptr)])
+    run.initialize(*fields(topology(A), np.float32))
+    run._step_count = 7
+    return run, module, allocations, storage
+
+
+class RuntimeChildRemap(unittest.TestCase):
+    def test_current_snapshot_identity_and_count_continuation(self):
+        run, module, allocations, storage = ready_backend()
+        # Make current state distinct from initialization and spare, as after a step.
+        for level in run._levels:
+            storage[id(level.current)] *= np.float32(1.01)
+        before = run.download()
+        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)
+        self.assertEqual(module.cuda.memcpy_dtoh.call_count, 3)
+        for a, b in zip(clone.download(), before):
+            np.testing.assert_array_equal(a, b)
+        self.assertEqual((run.step_count, clone.step_count), (7, 7))
+        for key, value in saved.items():
+            np.testing.assert_array_equal(storage[key], value)
+        for allocation in old_allocations:
+            allocation.free.assert_not_called()
+        clone.step()
+        self.assertEqual((run.step_count, clone.step_count), (7, 8))
+        run.step()
+        self.assertEqual(run.step_count, 8)
+        clone.close(); run.close()
+
+    def test_shift_uses_cpu_remap_and_resets_both_buffers(self):
+        run, module, allocations, storage = ready_backend()
+        before = run.download()
+        expected = remap(topology(A), topology(B), before).g
+        clone = run.remap_child(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):
+                np.testing.assert_array_equal(storage[id(ptr)], level.layout.to_blocked(want))
+        for a, b in zip(run.download(), before):
+            np.testing.assert_array_equal(a, b)
+        clone.initialize(*expected)
+        self.assertEqual(clone.step_count, 0)  # Existing public initialization semantics.
+        clone.close(); run.close()
+
+    def test_invalid_state_geometry_and_body_reject_before_reads_allocations(self):
+        run, module, allocations, storage = ready_backend()
+        start = len(allocations)
+        module.cuda.memcpy_dtoh.reset_mock()
+        for status in ('new', 'failed', 'closed', 'advancing'):
+            run._status = status
+            with self.assertRaises(RuntimeError): run.remap_child(B)
+        run._status = 'ready'
+        for bad in ([], [(40,16,48,40)], [(16,16,32,24),(16,27,32,40)]):
+            with self.assertRaises(ValueError): run.remap_child(bad)
+        hosts, pairs = run._boundary_hosts, 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
+                                 ('hierarchy', replace(run.hierarchy, pairs=(pairs[0], None)))):
+            with patch.object(run, attribute, value):
+                with self.assertRaises(ValueError): run.remap_child(B)
+        self.assertEqual(len(allocations), start)
+        module.cuda.memcpy_dtoh.assert_not_called()
+        run.close()
+
+    def test_candidate_construction_or_upload_failure_preserves_source(self):
+        for failure in ('construction', 'initialize'):
+            run, module, allocations, storage = ready_backend()
+            before = run.download()
+            old_allocations = list(allocations)
+            if failure == 'construction':
+                operation = patch.object(module.compiler, 'compile_kernel', side_effect=RuntimeError('candidate compile'))
+            else:
+                initialize = type(run).initialize
+                def fail(candidate, *arrays):
+                    upload = module.cuda.memcpy_htod.side_effect
+                    calls = 0
+                    def partial_upload(pointer, array):
+                        nonlocal calls
+                        calls += 1
+                        if calls == 2:
+                            raise RuntimeError('candidate upload')
+                        upload(pointer, array)
+                    with patch.object(module.cuda, 'memcpy_htod', side_effect=partial_upload):
+                        initialize(candidate, *arrays)
+                operation = patch.object(type(run), 'initialize', fail)
+            with operation, self.assertRaisesRegex(RuntimeError, 'candidate'):
+                run.remap_child(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()
+            for allocation in allocations[len(old_allocations):]: allocation.free.assert_called_once()
+            run.close()
+
+    def test_old_sync_or_read_failure_marks_only_old_failed(self):
+        for operation in ('sync', 'read'):
+            run, module, allocations, storage = ready_backend()
+            start = len(allocations)
+            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)
+            self.assertEqual(run._status, 'failed')
+            self.assertEqual(len(allocations), start)
+            target.side_effect = None
+            run.close()

tests/unit/test_shear_refinement_example.py

@@ -0,0 +1,84 @@
+"""CPU-only research-example contracts; no solver/GPU import required."""
+import unittest
+import numpy as np
+from research.shear_refinement.boundary_shear_tagging import tag_cells, strain_rate
+from research.shear_refinement.example import (CHILDREN, grow_rectangle, topology,
+                                                observe, initial_fields)
+
+
+class ShearRefinementExample(unittest.TestCase):
+    def test_affine_shear_physical_spacing_and_rotation(self):
+        y, x = np.indices((12, 20), dtype=float)
+        x, y = 7.+.5*(x+.5), 3.+.5*(y+.5)
+        valid = np.ones(x.shape, bool)
+        np.testing.assert_allclose(strain_rate(.003*y, 0*x, valid, .5), .003, atol=1e-17)
+        np.testing.assert_allclose(strain_rate(-.003*y, .003*x, valid, .5), 0., atol=1e-17)
+
+    def test_holes_are_unknown_not_walls_or_zero(self):
+        valid = np.zeros((7, 7), bool); valid[3, 3] = True
+        u = np.full(valid.shape, np.nan); u[3, 3] = .1
+        result = tag_cells(u, u, np.zeros_like(valid), valid=valid, owned=valid,
+                           dx_level=.5, origin=(7., 3.), U_ref=.05, L_ref=1.,
+                           shear_on=.01, shear_off=0., band_on=0., band_off=0., n_buffer=0)
+        self.assertTrue(result.unevaluable[3, 3])
+        self.assertTrue(np.isnan(result.strain_star[3, 3]))
+        self.assertFalse(result.tagged.any())
+        self.assertTrue(np.isinf(result.wall_distance).all())
+
+    def test_masked_hole_never_supplies_a_gradient(self):
+        y, x = np.indices((9, 9), dtype=float)
+        valid = np.ones(x.shape, bool); valid[:, 4] = False
+        ux = .003*(5.+.5*y); uy = np.zeros_like(ux)
+        ux[~valid] = np.inf; uy[~valid] = np.nan
+        strain = strain_rate(ux, uy, valid, .5)
+        np.testing.assert_allclose(strain[valid], .003, atol=1e-17)
+        self.assertTrue(np.isnan(strain[~valid]).all())
+
+    def test_predictable_growth_and_input_unchanged(self):
+        mask = np.zeros((64, 96), bool); mask[22, 74] = True
+        before = mask.copy()
+        self.assertEqual(grow_rectangle(CHILDREN[1], mask), (56, 16, 76, 40))
+        np.testing.assert_array_equal(mask, before)
+        mask[:] = False
+        self.assertEqual(grow_rectangle(CHILDREN[1], mask), CHILDREN[1])
+
+    def test_budget_and_geometry_rejection(self):
+        mask = np.zeros((64, 96), bool); mask[22, 74] = True
+        with self.assertRaisesRegex(ValueError, 'area cap'):
+            grow_rectangle(CHILDREN[1], mask, area_cap=479)
+        with self.assertRaises(ValueError):
+            grow_rectangle((57, 16, 72, 40), mask)
+        mask = np.zeros((9, 9), bool); mask[8, 8] = True
+        with self.assertRaisesRegex(ValueError, 'frame'):
+            grow_rectangle((0, 0, 4, 4), mask)
+        with self.assertRaises(ValueError):
+            topology((CHILDREN[0], (56, 16, 88, 40)))
+
+    def test_uniform_velocity_observation_does_not_grow(self):
+        from CelerisLab.lbm.equilibrium import compute_feq_d2q9
+        mesh = topology()
+        g = np.broadcast_to(compute_feq_d2q9(1., .01, 0.)[:, None, None],
+                            (9, *mesh.leaf_masks[1].shape)).copy()
+        g[:, ~mesh.leaf_masks[1]] = np.nan
+        rectangles, report = observe(g, mesh)
+        self.assertEqual(rectangles, CHILDREN)
+        self.assertEqual(report['requested'], 0)
+
+    def test_observation_growth_and_low_strain_no_change(self):
+        mesh = topology()
+        g = initial_fields(mesh)[1]
+        g[:, ~mesh.leaf_masks[1]] = np.nan
+        before = g.copy()
+        rectangles, report = observe(g, mesh)
+        self.assertEqual(rectangles, (CHILDREN[0], (56, 16, 76, 40)))
+        self.assertEqual(report['requested'], 64)
+        self.assertEqual(report['corridor_unknown'], 0)
+        self.assertEqual(report['corridor_unavailable'], 0)
+        self.assertGreater(report['out_of_scope'], 0)
+        np.testing.assert_array_equal(g, before)
+        rectangles, report = observe(g, mesh, shear_on=1.)
+        self.assertEqual(rectangles, CHILDREN)
+        self.assertEqual(report['requested'], 0)
+
+
+if __name__ == '__main__': unittest.main()

tests/unit/test_synchronized_remap.py

@@ -0,0 +1,205 @@
+"""CPU remap tests: full-population conservation, never a runtime AMR claim."""
+from dataclasses import replace
+import unittest
+
+import numpy as np
+
+from CelerisLab.config import LBMConfig
+from CelerisLab.lbm.blocks import BlockLayout2D
+from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
+from CelerisLab.lbm.refinement import FixedPatch2D, FixedPatchSet2D
+from CelerisLab.lbm.synchronized_remap import RemapTopology, remap_synchronized
+
+
+CFG = LBMConfig(nx=48, ny=32, viscosity=.1)
+PARENTS = [(4, 4, 20, 28), (24, 4, 44, 28)]
+A, B = [(16, 16, 32, 40)], [(20, 16, 36, 40)]
+
+
+def topology(rectangles, parent=None, cfg=CFG):
+    if parent is None:
+        parent = FixedPatchSet2D.from_config(cfg, PARENTS)
+    return RemapTopology.from_config(cfg, parent, rectangles)
+
+
+def fields(mesh, dtype=np.float64):
+    rng = np.random.default_rng(8473)
+    result = []
+    for leaf in mesh.leaf_masks:
+        # Direction-dependent and non-equilibrium, with independently varying
+        # children: q0 and momentum cannot be recovered from equilibrium rho/u.
+        g = rng.uniform(.03, .15, (9, *leaf.shape)).astype(dtype)
+        g += np.arange(9, dtype=dtype)[:, None, None] * .01
+        g[:, ~leaf] = np.nan
+        result.append(g)
+    return result
+
+
+def remap(old, new, g, **kwargs):
+    options = dict(base_count=7, level_steps=(7, 14, 28), phase='postcollision_g',
+                   origin=(-3., 2.), spacing=1.25)
+    options.update(kwargs)
+    return remap_synchronized(old, new, g, **options)
+
+
+def integrals(mesh, g, spacing=1.25):
+    populations = sum(a[:, mask].sum(axis=1, dtype=np.float64) * spacing**2 / 4**level
+                      for level, (a, mask) in enumerate(zip(g, mesh.leaf_masks)))
+    return np.r_[populations, populations.sum(), populations @ D2Q9_EX,
+                 populations @ D2Q9_EY]
+
+
+class SynchronizedRemapTests(unittest.TestCase):
+    def setUp(self):
+        self.old, self.new = topology(A), topology(B)
+        self.g = fields(self.old)
+
+    def test_identity_bitwise_and_no_mutation_or_count_advance(self):
+        copies = [a.copy() for a in self.g]
+        flags = [a.flags.writeable for a in self.g]
+        out = remap(self.old, self.old, self.g)
+        self.assertEqual((out.base_count, out.level_steps, out.phase, out.origin, out.spacing),
+                         (7, (7, 14, 28), 'postcollision_g', (-3., 2.), 1.25))
+        for a, b, saved, leaf, flag in zip(self.g, out.g, copies, self.old.leaf_masks, flags):
+            self.assertEqual(a[:, leaf].tobytes(), b[:, leaf].tobytes())
+            np.testing.assert_array_equal(a, saved)
+            self.assertEqual(a.flags.writeable, flag)
+            self.assertFalse(np.shares_memory(a, b))
+            self.assertTrue(np.isnan(b[:, ~leaf]).all())
+
+    def test_shift_preserves_physical_overlap_not_local_patch_index(self):
+        out = remap(self.old, self.new, self.g)
+        for level in range(3):
+            overlap = self.old.leaf_masks[level] & self.new.leaf_masks[level]
+            self.assertEqual(self.g[level][:, overlap].tobytes(), out.g[level][:, overlap].tobytes())
+        y, x = np.nonzero(self.new.leaf_masks[2] & ~self.old.leaf_masks[2])
+        np.testing.assert_array_equal(out.g[2][:, y, x], self.g[1][:, y//2, x//2])
+        y, x = np.nonzero(self.new.leaf_masks[1] & ~self.old.leaf_masks[1])
+        expected = sum(self.g[2][:, 2*y+dy, 2*x+dx] * .25
+                       for dy in (0, 1) for dx in (0, 1))
+        np.testing.assert_array_equal(out.g[1][:, y, x], expected)
+        self.assertGreater(np.max(np.abs(expected[0] - self.g[2][0, 2*y, 2*x])), .001)
+
+    def test_fp64_and_fp32_area_population_mass_momentum_conservation(self):
+        for dtype, tolerance in ((np.float64, 2e-13), (np.float32, 2e-7)):
+            g = fields(self.old, dtype)
+            out = remap(self.old, self.new, g)
+            before, after = integrals(self.old, g), integrals(self.new, out.g)
+            # FP64 accumulation; absolute scale handles near-cancelled momentum.
+            scale = np.sum(np.abs(before[:9]))
+            np.testing.assert_allclose(after, before, rtol=0, atol=tolerance * scale)
+            for array, leaf in zip(out.g, self.new.leaf_masks):
+                self.assertEqual(array.dtype, dtype)
+                self.assertTrue(np.isfinite(array[:, leaf]).all())
+                self.assertTrue(np.isnan(array[:, ~leaf]).all())
+
+    def test_layout_and_rectangle_permutation_independence(self):
+        old = topology(A + [(56, 16, 80, 40)])
+        new = topology(B + [(56, 16, 80, 40)])
+        reordered = topology(list(reversed(B + [(56, 16, 80, 40)])),
+                             FixedPatchSet2D.from_config(CFG, list(reversed(PARENTS))))
+        g = fields(old)
+        expected = remap(old, new, g)
+        for block in ((4, 4), (2, 2)):
+            layouts = (BlockLayout2D.build(CFG.ny, CFG.nx, *block),
+                       old.parent.fine_layout(block), old.child.fine_layout(block))
+            blocked = [layout.to_blocked(a) for layout, a in zip(layouts, g)]
+            # Corrupt all block halos; only owner-core populations are admissible.
+            for layout, a in zip(layouts, blocked):
+                h = layout.halo
+                a[:, :, :h, :] = np.inf
+                a[:, :, -h:, :] = np.inf
+                a[:, :, :, :h] = np.inf
+                a[:, :, :, -h:] = np.inf
+            saved = [a.copy() for a in blocked]
+            actual = remap(old, reordered, blocked, layouts=layouts)
+            for a, b in zip(expected.g, actual.g):
+                np.testing.assert_array_equal(a, b)
+            for a, b in zip(saved, blocked):
+                np.testing.assert_array_equal(a, b)
+
+    def test_shifted_parent_physical_frame(self):
+        # L1 window origin=(2,2) in root units; child coordinates remain parent-local.
+        cfg = LBMConfig(nx=32, ny=32, viscosity=.1)
+        parent = FixedPatch2D.from_config(cfg, (4, 4, 28, 28))
+        old, new = topology(A, parent, cfg), topology(B, parent, cfg)
+        g = fields(old)
+        out = remap(old, new, g, origin=(10., -7.))
+        self.assertEqual(out.topology.parent.fine_origin, (2., 2.))
+        overlap = old.leaf_masks[2] & new.leaf_masks[2]
+        np.testing.assert_array_equal(out.g[2][:, overlap], g[2][:, overlap])
+        np.testing.assert_allclose(integrals(old, g), integrals(new, out.g), rtol=2e-13, atol=2e-12)
+
+    def test_outside_leaf_nan_inf_ignored_required_leaf_rejected(self):
+        for level in range(3):
+            dirty = [a.copy() for a in self.g]
+            dirty[level][:, ~self.old.leaf_masks[level]] = np.inf
+            actual = remap(self.old, self.new, dirty)
+            expected = remap(self.old, self.new, self.g)
+            for a, b in zip(expected.g, actual.g):
+                np.testing.assert_array_equal(a, b)
+        # Specifically test restriction children, prolongation donors, and overlap.
+        required = [(0, self.old.leaf_masks[0]),
+                    (1, self.old.leaf_masks[1] & ~self.new.leaf_masks[1]),
+                    (2, self.old.leaf_masks[2] & ~self.new.leaf_masks[2])]
+        for level, mask in required:
+            y, x = np.argwhere(mask)[0]
+            for bad in (np.nan, np.inf, -np.inf):
+                dirty = [a.copy() for a in self.g]
+                dirty[level][0, y, x] = bad
+                with self.assertRaisesRegex(ValueError, 'required leaf donors'):
+                    remap(self.old, self.new, dirty)
+        dirty = [a.copy() for a in self.g]
+        dirty[0][:, 0, 0] = -1
+        with self.assertRaisesRegex(ValueError, 'positive density'):
+            remap(self.old, self.new, dirty)
+
+    def test_absent_required_block_rejected(self):
+        layouts = (None, None, BlockLayout2D.build(*self.old.leaf_masks[2].shape, 4, 4,
+                                                  periodic=(False, False), origins=[(0, 0)]))
+        arrays = [*self.g[:2], layouts[2].to_blocked(self.g[2])]
+        with self.assertRaisesRegex(ValueError, 'own every required leaf'):
+            remap(self.old, self.new, arrays, layouts=layouts)
+
+    def test_coarse_fine_coarse_identity_not_general_inverse(self):
+        # Grow then shrink: newly refined cells copy one parent value, so means
+        # recover it exactly. Shrink then grow loses arbitrary child variation.
+        wide = topology([(16, 16, 36, 40)])
+        refined = remap(self.old, wide, self.g)
+        restored = remap(wide, self.old, refined.g)
+        for a, b in zip(self.g, restored.g):
+            np.testing.assert_array_equal(a, b)
+        shifted = remap(self.old, self.new, self.g)
+        returned = remap(self.new, self.old, shifted.g)
+        lost = self.old.leaf_masks[2] & ~self.new.leaf_masks[2]
+        self.assertFalse(np.array_equal(self.g[2][:, lost], returned.g[2][:, lost]))
+
+    def test_invalid_geometry_counts_representation_frame(self):
+        for rectangles in ([], [(2, 2, 10, 10)], [(16, 16, 32, 40)] * 2,
+                           [(16, 16, 33.5, 40)], [(16, 16, 32, 40), (33, 16, 36, 40)]):
+            with self.assertRaises(ValueError):
+                topology(rectangles)
+        parent = self.old.parent
+        for changes in (dict(lattice_model='D3Q19'), dict(ddf_shifting=True),
+                        dict(collision='TRT'), dict(les_enabled=True), dict(nz=2)):
+            with self.assertRaises(ValueError):
+                RemapTopology.from_config(replace(CFG, **changes), parent, A)
+        for kwargs in (dict(coarse_boundary='slot_reservoir'), dict(objects=(object(),))):
+            with self.assertRaisesRegex(ValueError, 'periodic 2D'):
+                RemapTopology.from_config(CFG, parent, A, **kwargs)
+        different = topology(A, cfg=replace(CFG, viscosity=.12))
+        with self.assertRaisesRegex(ValueError, 'unchanged'):
+            remap(self.old, different, self.g)
+        for kwargs in (dict(base_count=True), dict(base_count=-1), dict(level_steps=(7, 14, 27)),
+                       dict(level_steps=(7, 14)), dict(level_steps=None), dict(phase='arrival_f'),
+                       dict(origin=(np.nan, 0)), dict(spacing=0), dict(spacing='bad')):
+            with self.assertRaises(ValueError):
+                remap(self.old, self.new, self.g, **kwargs)
+        for g in (self.g[:2], [self.g[0].astype(np.float32), *self.g[1:]],
+                  [self.g[0][0], *self.g[1:]]):
+            with self.assertRaises(ValueError):
+                remap(self.old, self.new, g)
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/validation/multigrid/test_multiple_pinned_remap_gpu.py

@@ -0,0 +1,175 @@
+"""Manual GPU gates for two pinned bodies and several fluid-only L2 regions.
+
+Run explicitly on an approved idle GPU. The CPU helper checks immediate
+whole-population transfer at the same physical time, not independent evolution
+accuracy. Continuation checks cover finite leaves and no-op bitwise identity;
+they do not establish an independent BFL replay or a numerical error bound.
+"""
+from contextlib import ExitStack
+import unittest
+
+import numpy as np
+
+from CelerisLab.body.geometry.circle import CircleGeometry
+from CelerisLab.body.jets import Jet, Parallel, Tube
+from CelerisLab.body.objects import SimObject
+from CelerisLab.cuda.context import CudaContext
+from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+from CelerisLab.lbm.synchronized_remap import (
+    RemapTopology, remap_synchronized, same_physical_boundary,
+)
+from test_refinement_blocked_base import exact
+from test_refinement_three_level_multipatch_gpu import CFG, PARENTS
+from test_runtime_child_remap_gpu import RuntimeChildRemapGPU, initial
+
+
+# All bounds use the allocated L1 frame, with support aligned to 4x4 blocks.
+PINS = ((16, 16, 32, 24), (56, 16, 72, 24))
+FLUID = ((16, 32, 32, 40), (56, 32, 72, 40))
+EXTRA = ((16, 44, 32, 48),)
+MOVED = ((20, 32, 36, 40), (60, 32, 76, 40))
+START = PINS + FLUID
+CHANGES = (
+    ('add', START + EXTRA, (.017, -.011)),
+    ('move', PINS + MOVED + EXTRA, (-.009, .005)),
+    ('remove', PINS + MOVED, (.013, -.007)),
+)
+
+
+def make(context):
+    centers = ((12., 10.), (32., 10.))
+    objects = tuple(
+        SimObject(71 + index, geometry=CircleGeometry(*center, .8),
+                  center=center, radius=.8)
+        for index, center in enumerate(centers)
+    )
+    jets = tuple(
+        Jet(Tube((x - 1.5, y), 0., 3., .3), Parallel(0.), 'parabolic', amplitude)
+        for (x, y), amplitude in zip(centers, (.003, -.002))
+    )
+    run = FixedPatchStepper(
+        CFG, PARENTS, context=context, base_block=(4, 4), fine_block=(4, 4),
+        child_block=(4, 4), child_rectangles_in_parent=START,
+        child_objects=objects, child_jets=jets,
+    )
+    return run, objects
+
+
+def topology(run, rectangles, objects):
+    return RemapTopology.from_config(
+        CFG, run.patch, rectangles, child_objects=objects,
+        child_jets=run._boundaries[2].jets,
+    )
+
+
+def population_integrals(mesh, arrays):
+    """Nine separate area integrals in the root frame, accumulated in FP64."""
+    return sum(array[:, leaf].sum(axis=1, dtype=np.float64) / 4**level
+               for level, (array, leaf) in enumerate(zip(arrays, mesh.leaf_masks)))
+
+
+class MultiplePinnedRemapGPU(unittest.TestCase):
+    check_public = RuntimeChildRemapGPU.check_public
+    check_reset = RuntimeChildRemapGPU.check_reset
+
+    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._boundaries[2]
+        body_ids = run.child_boundary.curved_host[3]
+        for index, amplitude in enumerate(amplitudes):
+            before = boundary.tables['surface_velocity'].copy()
+            run.set_jet_amplitude(index, amplitude, level=2)
+            after = boundary.tables['surface_velocity']
+            exact(self, after[body_ids != index], before[body_ids != index],
+                  'other body velocity unchanged')
+            self.assertTrue(np.any(after[body_ids == index] != before[body_ids == index]))
+        self.assertEqual(tuple(jet.amplitude for jet in boundary.jets), amplitudes)
+
+    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))
+        # Compare physical cut-link rows, never layout-dependent blocked addresses.
+        for run in (old, new):
+            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')
+
+    def test_noop_next_step_bitwise(self):
+        with CudaContext(0) as context, ExitStack() as owned:
+            old, objects = make(context)
+            owned.enter_context(old)
+            old.initialize(*initial(old, False))
+            old.step(2)
+            amplitudes = (.017, -.011)
+            self.change_amplitudes(old, amplitudes)
+            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))
+            self.assertIsNot(old, new)
+            self.assertEqual((old.step_count, new.step_count), (count, count))
+            self.check_boundary(old, new, mesh, amplitudes)
+            self.check_reset(new, before)
+            for got, want in zip(old.download(), before):
+                exact(self, got, want, 'no-op source unchanged')
+            old.step()
+            new.step()
+            self.assertEqual((old.step_count, new.step_count), (count + 1, count + 1))
+            for got, want in zip(new.download(), old.download()):
+                exact(self, got, want, 'no-op next step bitwise')
+            self.check_public(old)
+            self.check_public(new)
+
+    def test_add_move_remove_immediate_transfer_and_source_independence(self):
+        with CudaContext(0) as context, ExitStack() as owned:
+            current, objects = make(context)
+            owned.enter_context(current)
+            current.initialize(*initial(current, False))
+            current.step(2)
+            for action, rectangles, amplitudes in CHANGES:
+                with self.subTest(action=action):
+                    previous = current
+                    self.change_amplitudes(previous, amplitudes)
+                    before, count = previous.download(), previous.step_count
+                    old_mesh = topology(previous, previous.child_patch.bounds, objects)
+                    new_mesh = topology(previous, rectangles, objects)
+                    expected = remap_synchronized(
+                        old_mesh, new_mesh, before, base_count=count,
+                        level_steps=(count, 2 * count, 4 * count),
+                        phase='postcollision_g',
+                    ).g
+                    current = owned.enter_context(
+                        previous.remap_child(rectangles, child_objects=objects))
+                    self.assertIsNot(previous, current)
+                    self.assertEqual((previous.step_count, current.step_count), (count, count))
+                    self.check_boundary(previous, current, new_mesh, amplitudes)
+                    self.check_reset(current, expected)
+                    actual = current.download()
+                    for got, want, source, old_leaf, new_leaf in zip(
+                            actual, expected, before, old_mesh.leaf_masks, new_mesh.leaf_masks):
+                        exact(self, got, want, 'immediate CPU transfer parity')
+                        exact(self, got[:, old_leaf & new_leaf], source[:, old_leaf & new_leaf],
+                              'same-level physical fluid overlap')
+                        self.assertTrue(np.isfinite(got[:, new_leaf]).all())
+                        self.assertTrue(np.isnan(got[:, ~new_leaf]).all())
+                    # Conservation is asserted only across remap, not collision.
+                    np.testing.assert_allclose(
+                        population_integrals(new_mesh, actual),
+                        population_integrals(old_mesh, before),
+                        rtol=2 * np.finfo(np.float32).eps, atol=0,
+                        err_msg='all nine area-weighted population integrals',
+                    )
+                    self.check_public(current)
+                    current.step()
+                    self.assertEqual((previous.step_count, current.step_count), (count, count + 1))
+                    for got, want in zip(previous.download(), before):
+                        exact(self, got, want, 'source unchanged after candidate advances')
+                    self.check_boundary(previous, current, new_mesh, amplitudes)
+                    self.check_public(current)
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/validation/multigrid/test_pinned_body_remap_gpu.py

@@ -0,0 +1,84 @@
+"""Pinned-body GPU lifecycle gates. NOT executed by CPU test collection.
+
+No evolution accuracy tolerance is asserted here: an independent BFL replay
+reference is required before a 3e-6 FP32-vs-FP64 gate can be claimed. These
+checks require exact physical boundary/velocity, remap overlap and no-op next
+step identity; they never compare blocked boundary addresses across layouts.
+"""
+from contextlib import ExitStack
+import unittest
+import numpy as np
+from CelerisLab.cuda.context import CudaContext
+from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+from CelerisLab.lbm.synchronized_remap import RemapTopology, remap_synchronized
+from CelerisLab.body.geometry.circle import CircleGeometry
+from CelerisLab.body.objects import SimObject
+from CelerisLab.body.jets import Jet, Tube, Parallel
+from test_runtime_child_remap_gpu import RuntimeChildRemapGPU, initial
+from test_refinement_three_level_multipatch_gpu import CFG, PARENTS
+from test_refinement_blocked_base import exact
+
+A = ((16,16,32,40), (56,16,72,40))
+B = ((16,16,32,40), (60,16,76,40))
+
+
+def make(ctx):
+    body = SimObject(71, geometry=CircleGeometry(12.,14.,.8),center=(12.,14.),radius=.8)
+    jet = Jet(Tube((10.5,14.),0.,3.,.3),Parallel(0.),'parabolic',.003)
+    run = FixedPatchStepper(CFG,PARENTS,context=ctx,base_block=(4,4),fine_block=(4,4),
+                           child_block=(4,4),child_rectangles_in_parent=A,
+                           child_objects=(body,),child_jets=(jet,))
+    return run,(body,)
+
+
+class PinnedBodyRemapGPU(unittest.TestCase):
+    check_public = RuntimeChildRemapGPU.check_public
+    check_reset = RuntimeChildRemapGPU.check_reset
+
+    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,a,b,'physical cut links')
+        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:
+            old,objects=make(ctx);owned.enter_context(old)
+            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))
+            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')
+            self.assertEqual((old.step_count,new.step_count),(3,3))
+            self.check_public(new)
+
+    def test_moved_fluid_return_and_source_independence(self):
+        with CudaContext(0) as ctx, ExitStack() as owned:
+            current,objects=make(ctx);owned.enter_context(current)
+            current.initialize(*initial(current,False));current.step(2)
+            current.set_jet_amplitude(0,.017,level=2)
+            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,
+                                             child_objects=objects,child_jets=jets)
+                new=RemapTopology.from_config(CFG,previous.patch,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))
+                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')
+                current.step(2)
+                self.assertEqual((previous.step_count,current.step_count),(count,count+2))
+                for a,b in zip(previous.download(),before): exact(self,a,b,'source independent')
+                self.check_public(current)
+
+
+if __name__ == '__main__': unittest.main()

tests/validation/multigrid/test_raw_exchange_gpu.py

@@ -0,0 +1,264 @@
+"""GPU gates for raw exchange. Run explicitly; requires PyCUDA and CUDA.
+
+Independent census reads actual staged incoming/outgoing populations, including
+non-body ACTIVE-boundary flux, then checks postcollision fluid momentum/mass.
+This is a discrete implementation gate, NOT a physical nozzle-load validation.
+"""
+from dataclasses import replace
+import unittest
+import numpy as np
+import pycuda.driver as cuda
+
+from CelerisLab.body.geometry.circle import CircleGeometry
+from CelerisLab.body.objects import SimObject
+from CelerisLab.body.jets import Jet, Tube, Parallel
+from CelerisLab.config import LBMConfig
+from CelerisLab.cuda.context import CudaContext
+from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY, SOLID
+
+CX, CY = np.asarray(D2Q9_EX), np.asarray(D2Q9_EY)
+OPP = np.array([np.flatnonzero((CX == -x) & (CY == -y))[0] for x, y in zip(CX, CY)])
+W = np.array([4/9] + [1/9]*4 + [1/36]*4, np.float32)
+CFG = LBMConfig(nx=32, ny=24, viscosity=.1)
+BOUNDS = (8, 6, 24, 18)
+
+
+def circle():
+    return SimObject(17, geometry=CircleGeometry(16., 12., 2.15), center=(16., 12.), radius=2.15)
+
+
+def jet():
+    # Asymmetric aperture avoids hiding exchange errors by paired-jet cancellation.
+    return Jet(Tube((16., 12.), 0., 3., .7), Parallel(0.), amplitude=0.)
+
+
+def state(shape, ux=.013):
+    cu = CX[:, None, None]*ux
+    values = W[:, None, None]*(1+3*cu+4.5*cu*cu-1.5*ux*ux)
+    return np.broadcast_to(values, shape).astype(np.float32).copy()
+
+
+def download(pointer, shape, dtype=np.float32):
+    result = np.empty(shape, dtype)
+    cuda.memcpy_dtoh(result, pointer)
+    return result
+
+
+def totals(g, leaf):
+    values = g[:, leaf].astype(np.float64)
+    return np.array([(values*CX[:, None]).sum(), (values*CY[:, None]).sum(), values.sum()])
+
+
+def census(run, index, buffer):
+    """Independent incoming/outgoing face census on ACTIVE fluid cells.
+
+    Reads the receiver's actual block window, not a freshly gathered owner image
+    (that would erase BFL/restriction staging). Classifies external flux solely
+    from geometry, without consulting raw observation arrays.
+    """
+    level = run._levels[index]
+    layout, b = level.layout, level.boundary
+    g = download(buffer, layout.blocked_shape)
+    dense = layout.to_dense(g)
+    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)
+    flat = g.reshape(layout.nblocks, 9, cells)
+    external = np.zeros(3, np.float64)
+    for q, (cx, cy) in enumerate(zip(CX, CY)):
+        sy, sx = y-cy, x-cx
+        incoming = flat[block, q, cell-cy*layout.local_shape[1]-cx].astype(np.float64)
+        entering = ~leaf[sy, sx] & ~solid[sy, sx]
+        external += np.array([cx, cy, 1])*incoming[entering].sum()
+        dy, dx = y+cy, x+cx
+        leaving = ~leaf[dy, dx] & ~solid[dy, dx]
+        outgoing = flat[block, q, cell].astype(np.float64)
+        external -= np.array([cx, cy, 1])*outgoing[leaving].sum()
+    receiver, direction, _, _, rx, ry, _, _ = b.host.curved_host
+    owner = layout.owner.ravel()[receiver].astype(np.intp)
+    block, cell = np.divmod(owner, cells)
+    a = flat[block, direction, cell].astype(np.float64)
+    source = cell+CY[direction]*layout.local_shape[1]+CX[direction]
+    incoming = flat[block, OPP[direction], source].astype(np.float64)
+    jx, jy = CX[direction]*(a+incoming), CY[direction]*(a+incoming)
+    independent = np.column_stack((jx, jy, rx*jy-ry*jx, incoming-a))
+    return totals(dense, leaf), external, independent
+
+
+class RawExchangeGPU(unittest.TestCase):
+    def test_row_ledger_interval_jet_switch_and_no_population_effect(self):
+        with CudaContext(0) as ctx:
+            for block in ((32, 40), (16, 20), (8, 10)):
+                with FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_block=block,
+                        fine_objects=[circle()], jets=[jet()], observe_raw_exchange=True) as measured, \
+                     FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_block=block,
+                        fine_objects=[circle()], jets=[jet()]) as plain:
+                    arrays = [state(l.shape) for l in measured._levels]
+                    measured.initialize(*arrays)
+                    plain.initialize(*arrays)
+                    b = measured._level_boundaries[1]
+                    count = len(b.tables['cell'])
+                    kernel = measured._levels[1].kernel
+                    samples = []
+
+                    def audited_step(src, dst, *args, **kw):
+                        measured._stream.synchronize()
+                        before, external, independent = census(measured, 1, src)
+                        raw = download(b.last_raw, (count, 4))
+                        np.testing.assert_allclose(raw, independent, rtol=2e-6, atol=2e-7)
+                        kernel(src, dst, *args, **kw)
+                        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
+                        expected = external + np.array([-raw[:, 0].astype(float).sum(),
+                                                        -raw[:, 1].astype(float).sum(),
+                                                         raw[:, 3].astype(float).sum()])
+                        np.testing.assert_allclose(change, expected, rtol=2e-4, atol=8e-5)
+                        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)
+                        samples.clear()
+                        measured.step()
+                        plain.step()
+                        for a, e in zip(measured.download(), plain.download()):
+                            np.testing.assert_array_equal(a, e)
+                        integral = download(b.interval_integral, (count, 4), np.float64)
+                        expected = sum(s.astype(float)*[.25, .25, .125, .25] for s in samples)
+                        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)
+                        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'],
+                                                      obs['total']['impulse'])
+
+    def test_dense_fine_storage_reference(self):
+        """Same fine resolution/initial data with dense storage: a numerical
+        reference for cut-link rows, not a uniform-whole-domain convergence study."""
+        with CudaContext(0) as ctx:
+            with FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_objects=[circle()], jets=[jet()]) as dense, \
+                 FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_block=(16, 20),
+                    fine_objects=[circle()], jets=[jet()], observe_raw_exchange=True) as blocked:
+                arrays = [state(l.shape) for l in blocked._levels]
+                dense.initialize(*arrays)
+                blocked.initialize(*arrays)
+                for speed in (0., .002, 0.):
+                    dense.set_jet_amplitude(0, speed)
+                    blocked.set_jet_amplitude(0, speed)
+                    dense.step(); blocked.step()
+                    for a, e in zip(blocked.download(), dense.download()):
+                        np.testing.assert_array_equal(a, e)
+                    count = dense._fine_curved.count
+                    expected = download(dense._fine_curved.link_exchange_gpu, (count, 3))
+                    actual = download(blocked._level_boundaries[1].last_raw, (count, 4))
+                    np.testing.assert_array_equal(actual[:, :3], expected)
+
+    def test_depth_four_interval_storage_and_sampling(self):
+        cfg = LBMConfig(nx=48, ny=48, viscosity=.02)
+        obj = SimObject(0, geometry=CircleGeometry(16., 16., .5), center=(16., 16.), radius=.5)
+        drive = Jet(Tube((16., 16.), 0., .75, .2), Parallel(0.), amplitude=0.)
+        settings = dict(refinements=[[(6, 6, 26, 26)], [(20, 20, 44, 44)], [(56, 56, 72, 72)]],
+                        level_blocks=[(4, 4)]*4, level_objects=[(), (), (), (obj,)],
+                        level_jets=[(), (), (), (drive,)])
+        with CudaContext(0) as ctx:
+            with FixedPatchStepper(cfg, context=ctx, observe_raw_exchange=True, **settings) as run, \
+                 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]
+                count = len(b.tables['cell'])
+                kernel = run._levels[3].kernel
+                samples, independent_samples = [], []
+                def audited_step(src, dst, *args, **kw):
+                    run._stream.synchronize()
+                    before, external, independent = census(run, 3, src)
+                    raw = download(b.last_raw, (count, 4))
+                    np.testing.assert_allclose(raw, independent, rtol=2e-6, atol=2e-7)
+                    kernel(src, dst, *args, **kw)
+                    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
+                    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)
+                    samples.append(raw)
+                    independent_samples.append(independent)
+                run._levels[3].kernel = audited_step
+                for speed in (0., .002, 0.):
+                    run.set_jet_amplitude(0, speed, level=3)
+                    plain.set_jet_amplitude(0, speed, level=3)
+                    samples.clear(); independent_samples.clear()
+                    run.step(); plain.step()
+                    self.assertEqual(len(samples), 8)
+                    weights = np.array([1/64, 1/64, 1/512, 1/64])
+                    integral = download(b.interval_integral, (count, 4), np.float64)
+                    np.testing.assert_array_equal(integral, sum(v.astype(float)*weights for v in samples))
+                    np.testing.assert_allclose(integral, sum(v*weights for v in independent_samples),
+                                               rtol=2e-6, atol=2e-7)
+                    np.testing.assert_array_equal(download(b.last_raw, (count, 4)), samples[-1])
+                    for actual, expected in zip(run.download(), plain.download()):
+                        np.testing.assert_array_equal(actual, expected)
+                    result = run.read_raw_exchange(0, level=3)
+                    self.assertEqual(result['sample_count'], 8)
+                    self.assertEqual(result['time_start'], run.step_count-1)
+                    if run.step_count == 1:
+                        np.testing.assert_allclose(result['total']['force'], 0, atol=2e-7)
+                        self.assertAlmostEqual(result['total']['torque'], 0, delta=2e-7)
+
+    def test_uniform_whole_domain_fine_short_time_reference(self):
+        """Same root geometry, nu and Jet speed on a whole-domain fine mesh.
+        Eight local hops are a short-time equivalence check, not convergence.
+        The test adapter only enables the standalone kernel's optional inputs.
+        """
+        from CelerisLab.lbm.blocks import BlockLayout2D
+        from CelerisLab.lbm.block_stepper import PeriodicBlockStepper
+        fine_cfg = replace(CFG, nx=64, ny=48, viscosity=.2)
+        obj = SimObject(0, geometry=CircleGeometry(31.5, 23.5, 4.3),
+                        center=(31.5, 23.5), radius=4.3)
+        layout = BlockLayout2D.build(48, 64, 8, 8)
+        with CudaContext(0) as ctx:
+            with PeriodicBlockStepper(fine_cfg, layout, context=ctx, objects=[obj]) as uniform, \
+                 FixedPatchStepper(CFG, BOUNDS, context=ctx, fine_block=(16, 20),
+                    fine_objects=[circle()], jets=[jet()], 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]
+                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'])
+                np.testing.assert_allclose([v.q for v in links], b.tables['q'], atol=2e-7, rtol=0)
+                # Same fine-grid force arms after translation to the same center.
+                rx = uniform._upload(b.tables['rx'])
+                ry = uniform._upload(b.tables['ry'])
+                velocity = uniform._upload(b.tables['surface_velocity'])
+                raw = uniform._allocate(count*16)
+                original = uniform._curved
+                samples = []
+                def measured_boundary(*args, **kw):
+                    original(*args[:7], velocity, rx, ry, raw, **kw)
+                    uniform._stream.synchronize()
+                    samples.append(download(raw, (count, 4)))
+                uniform._curved = measured_boundary
+                for speed in (0., .002, -.001, 0.):
+                    amr.set_jet_amplitude(0, speed)
+                    cuda.memcpy_htod(velocity, b.tables['surface_velocity'])
+                    samples.clear()
+                    amr.step(); uniform.step(2)
+                    expected = sum(v.astype(float)*[.25, .25, .125, .25] for v in samples)
+                    actual = download(b.interval_integral, (count, 4), np.float64)
+                    np.testing.assert_allclose(actual, expected, rtol=1e-5, atol=2e-6)
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/validation/multigrid/test_refinement_gpu.py

@@ -194,7 +194,7 @@ class CoupledGPU(unittest.TestCase):
                               f"max_corner_packet={max_corner:.9g}; ownership/counts PASS", flush=True)
 
                 # Fail before any additional time step; no stale synchronized export.
-                with patch.object(backend, "_explode", side_effect=RuntimeError("injected launch failure")):
+                with patch.object(backend._pairs[0], "explode", side_effect=RuntimeError("injected launch failure")):
                     with self.assertRaisesRegex(RuntimeError, "injected"):
                         backend.step()
                 with self.assertRaisesRegex(RuntimeError, "failed"):

tests/validation/multigrid/test_refinement_multipatch.py

@@ -6,7 +6,8 @@ Run CPU and GPU classes in separate processes. GPU gates (frozen before any run)
 - S1 translation copy: two patches on 16x32 equal the 16x16 one-patch run tiled
   twice, bit for bit, every step;
 - reversed rectangle order gives bit-identical results;
-- face gap-4 pair (16x24): ACTIVE leaves finite, covered/hole masks exact.
+- face gap-4 pair (16x24): ACTIVE leaves finite, covered/hole masks exact, also on
+  raw owner cores per patch (INNER partly invalid by contract, not asserted).
 Weighted leaf mass is printed as a diagnostic only.
 """
 from dataclasses import replace
@@ -16,7 +17,7 @@ import numpy as np
 import pycuda.driver as cuda
 
 from CelerisLab.config import LBMConfig
-from CelerisLab.lbm.refinement import ACTIVE, FixedPatch2D, FixedPatchSet2D
+from CelerisLab.lbm.refinement import ACTIVE, OUTER_GHOST, FixedPatch2D, FixedPatchSet2D
 from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
 from test_refinement_blocked_fine import BOUNDS, NX, NY, config, state, circle
 from test_refinement_blocked_base import exact, raw
@@ -52,6 +53,13 @@ class MultipatchCPU(unittest.TestCase):
         with self.assertRaisesRegex(ValueError, "closer"):
             FixedPatchStepper(cfg, ((4, 4, 8, 12), (10, 4, 14, 12)), context=context,
                               fine_block=(8, 4), base_block=(4, 4))
+        from CelerisLab.body.jets import Jet, Parallel, Tube
+        cfg32 = LBMConfig(nx=32, ny=16, viscosity=.1)
+        span = Jet(Tube((5., 8.), 0., 22., .6), Parallel(0.), amplitude=0.)
+        with self.assertRaisesRegex(ValueError, "one jet"):  # same obj_id cannot bypass it
+            FixedPatchStepper(cfg32, ((2, 2, 14, 14), (18, 2, 30, 14)), context=context,
+                              fine_block=(4, 4), base_block=(4, 4),
+                              fine_objects=[body(8., 8.), body(24., 8.)], jets=[span])
 
 
 class MultipatchGPU(unittest.TestCase):
@@ -112,6 +120,9 @@ class MultipatchGPU(unittest.TestCase):
         fine = level_fine(s, [state(*p.fine_shape, 6 + i) for i, p in enumerate(s.patches)])
         kw = dict(fine_block=(4, 4), base_block=(4, 4))
         leaf = s.fine_kind == ACTIVE
+        support = s.fine_layout((4, 4)).owned
+        outer = support & (s.fine_kind == OUTER_GHOST)
+        holes = ~support
         with CudaContext(0) as context:
             runs = [self.run_steps(context, cfg, bounds, coarse, fine, **kw)
                     for bounds in (FACE, FACE[::-1])]
@@ -123,9 +134,93 @@ class MultipatchGPU(unittest.TestCase):
             self.assertTrue(np.isnan(c[:, s.covered]).all())
             self.assertTrue(np.isfinite(f[:, leaf]).all())
             self.assertTrue(np.isnan(f[:, ~leaf]).all())
+            # Raw owner cores, independent of the download mask (every patch).
+            # INNER may be partly NaN after hop 2 by contract: not asserted.
+            rc, rf = a[2], a[3]
+            label = f"raw step {step}"
+            self.assertTrue(np.isfinite(rc[:, ~s.covered]).all(), label + " coarse leaves")
+            if step:  # covered parent cores are poisoned by the parent step
+                self.assertTrue(np.isnan(rc[:, s.covered]).all(), label + " covered")
+            self.assertTrue(np.isfinite(rf[:, leaf]).all(), label + " fine ACTIVE")
+            self.assertTrue(np.isnan(rf[:, outer]).all(), label + " fine OUTER")
+            self.assertTrue(np.isnan(rf[:, holes]).all(), label + " fine holes")
+            for i, p in enumerate(s.patches):
+                active = window(s, i, rf)[:, p.fine_kind == ACTIVE]
+                self.assertTrue(np.isfinite(active).all(), f"{label} patch {i} ACTIVE")
             mass = float(np.nansum(c, dtype=np.float64) + .25 * np.nansum(f, dtype=np.float64))
             print(f"face gap4 step {step}: weighted leaf mass {mass:.9g} (diagnostic)", flush=True)
 
 
+def body(x, y, obj_id=0):
+    from CelerisLab.body.geometry.circle import CircleGeometry
+    from CelerisLab.body.objects import SimObject
+    return SimObject(obj_id, geometry=CircleGeometry(x, y, 2.15), center=(x, y), radius=2.15)
+
+
+class MultipatchBodiesGPU(unittest.TestCase):
+    """fine_objects/jets on several patches (frozen gates, bitwise)."""
+    PATCHES = ((2, 2, 14, 14), (18, 2, 30, 14))
+    STEPS = (0, 1, 2, 8)
+
+    def run_steps(self, context, cfg, bounds, coarse, fine, switch=None, **kw):
+        out = []
+        with FixedPatchStepper(cfg, bounds, context=context, **kw) as backend:
+            backend.initialize(coarse, fine)
+            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.step(step - backend.step_count)
+                out.append(backend.download() + raw(backend))
+        return out
+
+    def test_one_body_matches_single_patch_path(self):
+        from CelerisLab.cuda.context import CudaContext
+        with CudaContext(0) as context:
+            for mode in ("periodic", "slot_reservoir"):
+                cfg = config(mode)
+                p = FixedPatch2D.from_config(cfg, BOUNDS)
+                s = FixedPatchSet2D.from_config(cfg, (BOUNDS,))
+                coarse, fine = state(NY, NX, 1), state(*p.fine_shape, 2)
+                kw = dict(coarse_boundary=mode, fine_block=(8, 4), base_block=(3, 4),
+                          fine_objects=[circle()])
+                old = self.run_steps(context, cfg, BOUNDS, coarse, fine, **kw)
+                new = self.run_steps(context, cfg, (BOUNDS,), coarse, level_fine(s, [fine]), **kw)
+                for step, a, b in zip(self.STEPS, old, new):
+                    label = f"{mode} step {step}"
+                    exact(self, b[0], a[0], label + " coarse")
+                    exact(self, window(s, 0, b[1]), a[1], label + " fine")
+                    exact(self, window(s, 0, b[3]), a[3], label + " raw fine")
+
+    def test_two_bodies_and_jets_equal_tiled_single(self):
+        """Translation copy across a block seam, with jets and an amplitude switch."""
+        from CelerisLab.body.jets import Jet, Parallel, Tube
+        from CelerisLab.cuda.context import CudaContext
+        cfg16 = LBMConfig(nx=16, ny=16, viscosity=.1)
+        cfg32 = replace(cfg16, nx=32)
+        one = FixedPatchSet2D.from_config(cfg16, self.PATCHES[:1])
+        two = FixedPatchSet2D.from_config(cfg32, self.PATCHES)
+        coarse, fine = state(16, 16, 3), state(*one.patches[0].fine_shape, 4)
+        jet = lambda x: Jet(Tube((x, 8.), 0., 6., .6), Parallel(0.), amplitude=.002)
+        kw = dict(fine_block=(4, 4), base_block=(4, 4))
+        with CudaContext(0) as context:
+            for jets in ((), (jet(5.),)):
+                ref = self.run_steps(context, cfg16, self.PATCHES[:1], coarse, level_fine(one, [fine]),
+                                     switch=[(0, -.001)] if jets else None,
+                                     fine_objects=[body(8., 8.)], jets=jets, **kw)
+                both = jets and (jet(5.), jet(21.))
+                got = self.run_steps(context, cfg32, self.PATCHES, np.concatenate([coarse, coarse], 2),
+                                     level_fine(two, [fine, fine]),
+                                     switch=[(0, -.001), (1, -.001)] if jets else None,
+                                     fine_objects=[body(8., 8., 5), body(24., 8., 5)],
+                                     jets=both or (), **kw)
+                for step, a, b in zip(self.STEPS, ref, got):
+                    label = f"jets={len(jets)} step {step}"
+                    exact(self, b[0], np.concatenate([a[0], a[0]], axis=2), label + " coarse")
+                    for i in (0, 1):
+                        exact(self, window(two, i, b[1]), window(one, 0, a[1]), f"{label} fine {i}")
+                        exact(self, window(two, i, b[3]), window(one, 0, a[3]), f"{label} raw fine {i}")
+
+
 if __name__ == "__main__":
     unittest.main()

tests/validation/multigrid/test_refinement_three_level.py

@@ -0,0 +1,260 @@
+"""Nested three-level FixedPatchStepper (child=..., single chain): CPU guards and GPU gates.
+
+Fixture = reviewed CPU probe: L0 16x16 periodic, L1 (4,4,12,12) -> 24x24 window,
+L2 (8,8,16,16) in L1 window cells -> 24x24 window; 4x4 blocks on every level.
+Run CPU and GPU classes in separate processes. GPU gates (frozen before any run):
+- uniform moving equilibrium stays on all three levels' ACTIVE leaves, |dG| <= 1e-6;
+- seeded non-equilibrium: every level's ACTIVE leaves finite, raw owner cores NaN on
+  OUTER/covered (INNER partly invalid by contract, not asserted);
+- each L1 hop feeds L2 from the L1 buffer it is about to step (current, then the new
+  G written by the first hop), never from L0: 2 child advances per coarse step.
+Numerical comparison against the CPU probe (G vs C(f)) is done by the acceptance run.
+"""
+from contextlib import nullcontext
+from dataclasses import replace
+import unittest
+from unittest.mock import Mock, patch
+
+import numpy as np
+import pycuda.driver as cuda
+
+from CelerisLab.config import LBMConfig
+from CelerisLab.lbm.equilibrium import compute_feq_d2q9
+from CelerisLab.lbm.refinement import ACTIVE, OUTER_GHOST
+from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+from test_refinement_blocked_fine import state, circle
+from test_refinement_blocked_base import raw
+
+CFG = LBMConfig(nx=16, ny=16, viscosity=.1)
+L1, L2 = (4, 4, 12, 12), (8, 8, 16, 16)
+KW = dict(fine_block=(4, 4), base_block=(4, 4), child=L2, child_block=(4, 4))
+
+
+def child_body(spec=None):
+    from CelerisLab.body.geometry.circle import CircleGeometry
+    from CelerisLab.body.objects import SimObject
+    return SimObject(99, geometry=CircleGeometry(8., 8., .8), center=(8., 8.), radius=.8,
+                     surface_actuation=spec)
+
+
+def child_jet(speed=.002):
+    from CelerisLab.body.jets import Jet, Tube, Parallel
+    return Jet(Tube((6.5, 8.), 0., 3., .3), Parallel(0.), amplitude=speed)
+
+
+def mocked_child():
+    from CelerisLab.cuda import compiler_v2 as compiler
+    context = Mock(activate=nullcontext, sm_arch="sm_70")
+    module = Mock(get_function=Mock(side_effect=lambda name: Mock(name=name)))
+    module.get_global.return_value = (Mock(), 20)
+    with patch.object(cuda, "Stream"), patch.object(compiler, "generate_config"), \
+            patch.object(compiler, "compile_kernel"), patch.object(compiler, "cleanup_build_dir"), \
+            patch.object(compiler, "load_module", return_value=module), \
+            patch.object(cuda, "mem_alloc", side_effect=lambda size: Mock()), \
+            patch.object(cuda, "memcpy_htod"):
+        return FixedPatchStepper(CFG, L1, context=context, child_objects=[child_body()],
+                                 child_jets=[child_jet()], **KW)
+
+
+class ThreeLevelCPU(unittest.TestCase):
+    def test_rejections_before_cuda(self):
+        context = object()  # never activated: every case must fail first
+        bad = (dict(KW, child_block=None), dict(KW, fine_block=None),
+               dict(KW, fine_objects=[circle()]), dict(KW, child=(4, 8, 16, 16)),
+               dict(KW, child_block=(5, 4)))
+        for kw in bad:
+            with self.subTest(kw=kw), self.assertRaises(ValueError):
+                FixedPatchStepper(CFG, L1, context=context, **kw)
+        for kw in (dict(child_objects=[child_body()]), dict(child_jets=[child_jet()]),
+                   dict(KW, child_jets=[child_jet()])):
+            with self.assertRaises(ValueError):
+                FixedPatchStepper(CFG, L1, context=context, **kw)
+        with self.assertRaises(ValueError):  # several L1 patches + child: not yet
+            FixedPatchStepper(LBMConfig(nx=32, ny=16, viscosity=.1), ((4, 4, 12, 12), (20, 4, 28, 12)),
+                              context=context, **KW)
+
+
+    def test_four_halo_bfl_steps_and_fresh_middle_feeds(self):
+        run = mocked_child()
+        run._status = "ready"
+        f, g = run._levels[1:]
+        b = g.boundary
+        events, feeds = [], []
+        advance, gather = run._advance_child, run._gather_halo
+        def record_feed(buffer, **kw):
+            self.assertIs(buffer, (f.current, f.spare)[len(feeds)])
+            feeds.append(buffer)
+            advance(buffer, **kw)
+        def record_halo(level, buffer, **kw):
+            if level is g:
+                events.append(("halo", buffer))
+            gather(level, buffer, **kw)
+        b.kernel.side_effect = lambda buffer, *a, **k: events.append(("bfl", buffer))
+        def record_step(src, dst, halo, flags, kind, **kw):
+            self.assertIs(flags, b.flags)
+            self.assertIs(kind, g.kind)
+            events.append(("step", src))
+        g.kernel.side_effect = record_step
+        with patch.object(run, "_advance_child", side_effect=record_feed), \
+                patch.object(run, "_gather_halo", side_effect=record_halo):
+            run.step()
+        self.assertEqual(feeds, [f.current, f.spare])
+        expected = [(kind, buffer) for buffer in (g.current, g.spare, g.current, g.spare)
+                    for kind in ("halo", "bfl", "step")]
+        self.assertEqual(events, expected)
+        self.assertGreater(len(np.unique(b.tables["block"])), 1)
+        # Negative control: reusing L1's old source at the midpoint must fail.
+        feeds.clear()
+        with patch.object(run, "_advance_child", side_effect=lambda buffer, **kw:
+                          record_feed(f.current, **kw)):
+            with self.assertRaises(AssertionError):
+                run.step()
+        run.close()
+
+    def test_child_update_isolated_and_failure_atomic(self):
+        run = mocked_child()
+        b = run._levels[2].boundary
+        before = b.tables["surface_velocity"].copy()
+        with self.assertRaises(IndexError):
+            run.set_jet_amplitude(0, -.002)  # default stays L1
+        for level in (0, 3, True):
+            with self.assertRaises(ValueError):
+                run.set_jet_amplitude(0, -.002, level=level)
+        with patch.object(cuda, "memcpy_htod") as upload:
+            run.set_jet_amplitude(0, -.002, level=2)
+            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, ())
+        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)
+        self.assertEqual(run._status, "failed")
+        self.assertEqual(b.jets[0].amplitude, -.002)
+        np.testing.assert_array_equal(b.tables["surface_velocity"], -before)
+        run.close()
+
+    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
+        self.assertTrue(solid.any())
+        self.assertFalse(run._child_leaf[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())
+        run.close()
+
+
+class ThreeLevelGPU(unittest.TestCase):
+    def test_equilibrium_is_kept_on_every_level(self):
+        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, L1, context=ctx, **KW) as run:
+            shapes = ((16, 16), run.patch.fine_shape, run.child_patch.fine_shape)
+            run.initialize(*(np.broadcast_to(feq, (9, *s)).copy() for s in shapes))
+            run.step(8)
+            self.assertEqual(run.step_count, 8)
+            for level, g in enumerate(run.download()):
+                leaf = np.isfinite(g[0])
+                self.assertTrue(leaf.any(), f"level {level}")
+                np.testing.assert_allclose(g[:, leaf], np.broadcast_to(feq, g.shape)[:, leaf],
+                                           atol=1e-6, rtol=0, err_msg=f"level {level}")
+
+    def test_nonequilibrium_masks_and_fresh_middle_feeds(self):
+        from CelerisLab.cuda.context import CudaContext
+        with CudaContext(0) as ctx, FixedPatchStepper(CFG, L1, context=ctx, **KW) as run:
+            p1, p2 = run.patch, run.child_patch
+            run.initialize(state(16, 16, 1), state(*p1.fine_shape, 2), state(*p2.fine_shape, 3))
+            f = run._levels[1]
+            feeds, original = [], run._advance_child
+            def record(buffer, **kw):
+                feeds.append(int(buffer))
+                original(buffer, **kw)
+            with patch.object(run, "_advance_child", side_effect=record):
+                cur, spare = int(f.current), int(f.spare)
+                run.step(1)
+            self.assertEqual(feeds, [cur, spare])  # current G, then the new G of hop 1
+            run.step(7)
+            middle = run._middle_kind
+            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())
+            rc, rf, rg = raw(run)
+            self.assertTrue(np.isnan(rc[:, p1.covered]).all())
+            self.assertTrue(np.isfinite(rf[:, middle == ACTIVE]).all())
+            self.assertTrue(np.isnan(rf[:, middle == OUTER_GHOST]).all())  # incl. child-covered
+            self.assertTrue(np.isnan(rf[:, p2.covered]).all())
+            self.assertTrue(np.isfinite(rg[:, p2.fine_kind == ACTIVE]).all())
+            self.assertTrue(np.isnan(rg[:, p2.fine_kind == OUTER_GHOST]).all())
+
+
+class ThreeLevelBodyGPU(unittest.TestCase):
+    """Frozen gates; CPU mocks do not establish these numerical results."""
+    def run_case(self, mode, block, *, jets=(), spec=None, switch=False):
+        from CelerisLab.cuda.context import CudaContext
+        from CelerisLab.lbm.descriptors import SOLID
+        cfg = CFG if mode == "periodic" else replace(
+            CFG, velocity=0., wall_normal_velocity=.001, wall_slot_width=4,
+            y_wall_bc="on_site_velocity", top_wall_bc="reservoir")
+        with CudaContext(0) as ctx, FixedPatchStepper(
+                cfg, L1, context=ctx, coarse_boundary=mode, child_objects=[child_body(spec)],
+                child_jets=jets, **dict(KW, child_block=block)) as run:
+            run.initialize(*(state(*level.shape[1:], i + 1) for i, level in enumerate(run._levels)))
+            snapshots = []
+            for stop in (0, 1, 2, 4, 8):
+                run.step(stop - run.step_count)
+                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]
+                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)]
+                if mode != "periodic":
+                    invalid[0][[0, -1]] = True
+                for g, mask in zip(stored, invalid):
+                    self.assertTrue(np.isnan(g[:, mask]).all())
+                snapshots.append(public + tuple(stored))
+            return snapshots
+
+    def test_layouts_static_driven_switched_and_zero(self):
+        for mode in ("periodic", "slot_reservoir"):
+            cases = [self.run_case(mode, (4, 4), jets=jets, switch=switch)
+                     for jets, switch in (((), False), ([child_jet()], False),
+                                          ([child_jet()], True), ([child_jet(0.)], False))]
+            for a, b in zip(cases[0], cases[3]):
+                for x, y in zip(a, b):
+                    np.testing.assert_array_equal(x, y)
+            self.assertFalse(np.array_equal(cases[1][-1][2], cases[2][-1][2], equal_nan=True))
+            for block in ((8, 4), (4, 8)):
+                for expected, (jets, switch) in zip(cases, (((), False), ([child_jet()], False),
+                                                          ([child_jet()], True), ([child_jet(0.)], False))):
+                    actual = self.run_case(mode, block, jets=jets, switch=switch)
+                    for a, b in zip(actual, expected):
+                        for x, y in zip(a, b):
+                            np.testing.assert_array_equal(x, y)
+
+    def test_legacy_surface_equivalence(self):
+        from CelerisLab.body.objects import SurfaceActuation
+        spec = SurfaceActuation(.7, .85, .002)
+        for mode in ("periodic", "slot_reservoir"):
+            legacy = self.run_case(mode, (4, 4), spec=spec)[-1][2]
+            driven = self.run_case(mode, (4, 4), jets=[spec.as_jet((8., 8.), .8)])[-1][2]
+            mask = np.isfinite(legacy)
+            np.testing.assert_array_equal(mask, np.isfinite(driven))
+            self.assertLessEqual(np.max(np.abs(legacy[mask] - driven[mask])), 1e-6)
+
+
+if __name__ == "__main__":
+    unittest.main()

tests/validation/multigrid/test_refinement_three_level_multipatch_gpu.py

@@ -0,0 +1,178 @@
+"""Frozen GPU gates for the bounded three-level separated-union backend.
+
+Run explicitly on an available GPU; CPU geometry/mocks live in tests/unit.
+Snapshots 0/1/2/8: finite public leaves, exact NaN masks, raw invalid cores;
+reversed order and aligned layout changes must be bit-identical. The scalar
+window comparison uses root physical coordinates, never equal array indices.
+These equivalence checks supplement, not replace, independent dense CPU replay.
+"""
+from dataclasses import replace
+import unittest
+
+import numpy as np
+
+from CelerisLab.config import LBMConfig
+from CelerisLab.lbm.descriptors import SOLID
+from CelerisLab.lbm.equilibrium import compute_feq_d2q9
+from CelerisLab.lbm.refinement import OUTER_GHOST
+from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+from test_refinement_blocked_base import raw, exact
+
+CFG = LBMConfig(nx=48, ny=32, viscosity=.1)
+PARENTS = ((4, 4, 20, 28), (24, 4, 44, 28))
+CHILDREN = ((16, 16, 32, 40), (56, 16, 80, 40))
+TOPOLOGIES = {
+    'twins': ((16, 16, 32, 40), (56, 16, 72, 40)),
+    'asymmetric': (CHILDREN[0],),
+    'multiple_children': ((16, 16, 32, 24), (16, 28, 32, 40)),
+}
+# 8-wide blocks require support edges divisible by 8, i.e. rectangle edges
+# equal 2 modulo 4 in each pair's parent frame. The default fixture is 4-only.
+ALIGNED_CFG = LBMConfig(nx=56, ny=40, viscosity=.1)
+ALIGNED_PARENTS = ((6, 6, 22, 34), (26, 6, 50, 34))
+ALIGNED_CHILDREN = ((18, 18, 34, 50), (58, 18, 90, 50))
+STOPS = (0, 1, 2, 8)
+
+
+def config(cfg, mode):
+    return cfg if mode == 'periodic' else replace(
+        cfg, velocity=0., wall_normal_velocity=.001, wall_slot_width=4,
+        y_wall_bc='on_site_velocity', top_wall_bc='reservoir')
+
+
+def bodies():
+    from CelerisLab.body.geometry.circle import CircleGeometry
+    from CelerisLab.body.objects import SimObject
+    return [SimObject(99, geometry=CircleGeometry(x, 14., .8), center=(x, 14.), radius=.8)
+            for x in (12., 34.)]
+
+
+def jets(speed):
+    from CelerisLab.body.jets import Jet, Tube, Parallel
+    return [Jet(Tube((x-1.5, 14.), 0., 3., .3), Parallel(0.), amplitude=speed)
+            for x in (12., 34.)]
+
+
+def frames(run):
+    middle_origin = np.asarray(run.patch.fine_origin)
+    return ((np.zeros(2), 1.), (middle_origin, .5),
+            (middle_origin + .5*np.asarray(run.child_patch.fine_origin), .25))
+
+
+def physical_state(shape, origin, spacing):
+    """Same smooth non-equilibrium G at the same physical cell centers."""
+    y, x = np.indices(shape, dtype=np.float64)
+    x, y = origin[0] + spacing*(x+.5), origin[1] + spacing*(y+.5)
+    equilibrium = np.asarray(compute_feq_d2q9(1., .013, -.007), np.float64)[:, None, None]
+    q = np.arange(9)[:, None, None]
+    return (equilibrium + 1e-4*np.sin(.17*x + .23*y + q)).astype(np.float32)
+
+
+def physical_leaves(array, origin, spacing):
+    y, x = np.nonzero(np.isfinite(array[0]))
+    coordinates = np.stack((origin[1]+spacing*(y+.5), origin[0]+spacing*(x+.5)), axis=-1)
+    return coordinates, array[:, y, x]
+
+
+class ThreeLevelMultipatchGPU(unittest.TestCase):
+    def run_case(self, ctx, cfg=CFG, parents=PARENTS, children=CHILDREN,
+                 mode='periodic', blocks=((4, 4), (4, 4)), objects=(), drive=(), switch=False,
+                 scalar_child=None):
+        kwargs = dict(child_rectangles_in_parent=children) if scalar_child is None else dict(child=scalar_child)
+        with FixedPatchStepper(config(cfg, mode), parents, context=ctx, coarse_boundary=mode,
+                base_block=(4, 4), fine_block=blocks[0], child_block=blocks[1],
+                child_objects=objects, child_jets=drive, **kwargs) as run:
+            frame = frames(run)
+            run.initialize(*(physical_state(level.shape[1:], *f)
+                             for level, f in zip(run._levels, frame)))
+            snapshots = []
+            for stop in STOPS:
+                run.step(stop-run.step_count)
+                if switch and stop == 2:
+                    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]
+                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]
+                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
+                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)
+                    self.assertTrue(np.isnan(g[:, ~mask]).all(), label)
+                    self.assertTrue(np.isnan(dense[:, bad]).all(), label)
+                snapshots.append((public, stored, frame))
+            return snapshots
+
+    def compare(self, actual, expected):
+        for stop, a, b in zip(STOPS, actual, expected):
+            for kind in (0, 1):
+                for level, (x, y) in enumerate(zip(a[kind], b[kind])):
+                    exact(self, x, y, f'stop={stop} kind={kind} level={level}')
+
+    def test_topologies_reversed_order_and_masks(self):
+        from CelerisLab.cuda.context import CudaContext
+        with CudaContext(0) as ctx:
+            for mode in ('periodic', 'slot_reservoir'):
+                for name, children in TOPOLOGIES.items():
+                    with self.subTest(mode=mode, topology=name):
+                        expected = self.run_case(ctx, children=children, mode=mode)
+                        actual = self.run_case(ctx, parents=PARENTS[::-1], children=children[::-1], mode=mode)
+                        self.compare(actual, expected)
+
+    def test_independently_varied_aligned_layouts(self):
+        from CelerisLab.cuda.context import CudaContext
+        with CudaContext(0) as ctx:
+            for mode in ('periodic', 'slot_reservoir'):
+                kwargs = dict(cfg=ALIGNED_CFG, parents=ALIGNED_PARENTS,
+                              children=ALIGNED_CHILDREN, mode=mode)
+                expected = self.run_case(ctx, **kwargs)
+                for middle in ((4, 4), (8, 4), (4, 8)):
+                    for leaf in ((4, 4), (8, 4), (4, 8)):
+                        self.compare(self.run_case(ctx, blocks=(middle, leaf), **kwargs), expected)
+
+    def test_one_element_equivalent_physical_frames(self):
+        from CelerisLab.cuda.context import CudaContext
+        cfg = LBMConfig(nx=16, ny=16, viscosity=.1)
+        parent = (4, 4, 12, 12)
+        with CudaContext(0) as ctx:
+            for mode in ('periodic', 'slot_reservoir'):
+                old = self.run_case(ctx, cfg=cfg, parents=parent, mode=mode, scalar_child=(8, 8, 16, 16))
+                # Test both set-on-scalar and set-on-set allocated frames.
+                for parents, child in ((parent, (8, 8, 16, 16)), ((parent,), (12, 12, 20, 20))):
+                    new = self.run_case(ctx, cfg=cfg, parents=parents, children=[child], mode=mode)
+                    for a, b in zip(new, old):
+                        for level in range(3):
+                            ax, av = physical_leaves(a[0][level], *a[2][level])
+                            bx, bv = physical_leaves(b[0][level], *b[2][level])
+                            np.testing.assert_array_equal(ax, bx)
+                            exact(self, av, bv, f'physical level={level}')
+
+    def test_two_l2_bodies_zero_driven_switched_order_and_layout(self):
+        from CelerisLab.cuda.context import CudaContext
+        with CudaContext(0) as ctx:
+            for mode in ('periodic', 'slot_reservoir'):
+                kwargs = dict(cfg=ALIGNED_CFG, parents=ALIGNED_PARENTS,
+                              children=ALIGNED_CHILDREN, mode=mode, objects=bodies())
+                static = self.run_case(ctx, **kwargs)
+                zero = self.run_case(ctx, drive=jets(0.), **kwargs)
+                self.compare(zero, static)
+                driven = self.run_case(ctx, drive=jets(.002), **kwargs)
+                switched = self.run_case(ctx, drive=jets(.002), switch=True, **kwargs)
+                self.assertFalse(np.array_equal(driven[-1][0][2], static[-1][0][2], equal_nan=True))
+                self.assertFalse(np.array_equal(driven[-1][0][2], switched[-1][0][2], equal_nan=True))
+                for drive, switch, expected in (((), False, static), (jets(0.), False, zero),
+                                                (jets(.002), False, driven), (jets(.002), True, switched)):
+                    for middle, leaf in (((8, 4), (4, 8)), ((4, 8), (8, 4))):
+                        actual = self.run_case(ctx, blocks=(middle, leaf), drive=drive, switch=switch,
+                            **dict(kwargs, parents=ALIGNED_PARENTS[::-1], children=ALIGNED_CHILDREN[::-1]))
+                        self.compare(actual, expected)
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/validation/multigrid/test_runtime_child_remap_gpu.py

@@ -0,0 +1,106 @@
+"""Manual L2 continuation gates; run explicitly on an approved idle GPU.
+
+Frozen numerical gate: max |G - independent FP64 replay| <= 3e-6 on all
+ACTIVE leaves, with exact masks and exact same-level overlap at remap. The
+separate independent capture/replay acceptance harness supplies that numerical
+comparison; these tests cover lifecycle, exact transfer and GPU continuation.
+Remap mass/momentum must agree to FP32 roundoff, not per-q conservation during
+advancement (collision changes individual populations). No GPU run is implied.
+"""
+from contextlib import ExitStack
+import unittest
+
+import numpy as np
+import pycuda.driver as cuda
+
+from CelerisLab.cuda.context import CudaContext
+from CelerisLab.lbm.equilibrium import compute_feq_d2q9
+from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+from CelerisLab.lbm.synchronized_remap import RemapTopology, remap_synchronized
+from test_refinement_blocked_base import exact
+from test_refinement_three_level_multipatch_gpu import CFG, PARENTS, frames, physical_state
+
+A = ((16, 16, 32, 40),)
+B = ((20, 16, 36, 40),)
+FP64_ATOL = 3e-6
+
+
+def make(ctx):
+    return FixedPatchStepper(CFG, PARENTS, context=ctx, base_block=(4, 4),
+                             fine_block=(4, 4), child_block=(4, 4),
+                             child_rectangles_in_parent=A)
+
+
+def initial(run, uniform):
+    equilibrium = np.asarray(compute_feq_d2q9(1., .013, -.007), np.float32)[:, None, None]
+    return tuple(np.broadcast_to(equilibrium, level.shape).copy() if uniform else
+                 physical_state(level.shape[1:], *frame)
+                 for level, frame in zip(run._levels, frames(run)))
+
+
+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)):
+            self.assertTrue(np.isfinite(array[:, leaf]).all())
+            self.assertTrue(np.isnan(array[:, ~leaf]).all())
+
+    def check_reset(self, run, expected):
+        for level, want in zip(run._levels, expected):
+            blocked = level.layout.to_blocked(want)
+            for pointer in (level.current, level.spare):
+                actual = np.empty(level.layout.blocked_shape, np.float32)
+                cuda.memcpy_dtoh(actual, pointer)
+                exact(self, actual, blocked, 'fresh current and spare, including halos')
+
+    def test_identity_next_step_bitwise_and_old_still_usable(self):
+        with CudaContext(0) as ctx, ExitStack() as owned:
+            old = owned.enter_context(make(ctx))
+            old.initialize(*initial(old, False))
+            old.step(2)
+            before = old.download()
+            new = owned.enter_context(old.remap_child(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')
+            old.step(); new.step()
+            self.assertEqual((old.step_count, new.step_count), (3, 3))
+            for got, want in zip(new.download(), old.download()):
+                exact(self, got, want, 'identity next step')
+            self.check_public(new)
+
+    def test_shift_then_return_uniform_and_nonequilibrium(self):
+        for uniform in (True, False):
+            with self.subTest(uniform=uniform), CudaContext(0) as ctx, ExitStack() as owned:
+                current = owned.enter_context(make(ctx))
+                current.initialize(*initial(current, uniform))
+                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)
+                    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))
+                    self.assertEqual(current.step_count, count)
+                    self.check_reset(current, expected)
+                    for got, want, old_g, old_leaf, new_leaf in zip(
+                            current.download(), expected, before, old_mesh.leaf_masks, new_mesh.leaf_masks):
+                        exact(self, got, want, 'CPU remap exact')
+                        overlap = old_leaf & new_leaf
+                        exact(self, got[:, overlap], old_g[:, overlap], 'same-level overlap')
+                    for got, want in zip(previous.download(), before):
+                        exact(self, got, want, 'source unchanged')
+                    current.step(2)
+                    self.assertEqual(current.step_count, count+2)
+                    self.assertEqual(previous.step_count, count)
+                    self.check_public(current)
+                    if uniform:
+                        equilibrium = np.asarray(compute_feq_d2q9(1., .013, -.007), np.float64)
+                        for array in current.download():
+                            leaf = np.isfinite(array[0])
+                            self.assertLessEqual(np.max(np.abs(array[:, leaf]-equilibrium[:, None])), FP64_ATOL)
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/validation/multigrid/test_shear_refinement_example_gpu.py

@@ -0,0 +1,43 @@
+"""Run explicitly on an idle GPU; automatic/manual identical-layout equivalence.
+
+This is not independent flow-accuracy validation. CPU collection never runs it.
+PYTHONPATH=src:.:tests/validation/multigrid python -m unittest discover \
+    -s tests/validation/multigrid -p test_shear_refinement_example_gpu.py -v
+"""
+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)
+
+
+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))
+            initialize(old)
+            old.step(2)
+            before = old.download()
+            chosen, report = observe(before[1], topology())
+            manual = (CHILDREN[0], (56, 16, 76, 40))
+            self.assertEqual(chosen, manual)
+            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()))
+            for run in (auto, explicit):
+                self.assertEqual(run.child_patch.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):
+                np.testing.assert_array_equal(a, b)
+                self.assertTrue(np.isfinite(a[:, leaf]).all())
+                self.assertTrue(np.isnan(a[:, ~leaf]).all())
+            for a, b in zip(old.download(), before):
+                np.testing.assert_array_equal(a, b)
+            self.assertEqual((old.step_count, auto.step_count, explicit.step_count), (2, 4, 4))
+
+
+if __name__ == '__main__': unittest.main()