Document bounded multigrid and add public remap example

prepare_fine_boundary accepts a FixedPatchSet2D: objects are projected once
onto the level-wide fine grid (FixedPatchSet2D.fine_origin = (0, 0), same
coarse cell-edge -> fine cell transform), each circle must lie inside one
patch, coarse transfer cells are checked against the covered union, and body
ids stay object positions, so the existing jet one-body-per-jet, overlap and
(receiver, dir) writer checks run once over the whole level (caller obj_id
cannot bypass them). The stepper no longer rejects fine_objects/jets with
several patches; blocked link tables, BFL and the jet amplitude cache are
reused unchanged. No kernel change. README/docstring scope updated.

Tests: CPU translation/ids/rejections (outside one patch, one jet over two
bodies at zero amplitude); GPU one body via the new path == single-patch path
(both modes), two bodies (+jets, amplitude switch) == tiled single (bitwise).

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 patches, two or three levels with separated rectangles.
 - `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,87 @@ 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. Wall-touching patches are unsupported.
+  Three-level nesting uses the explicit child coordinates described below.
+
+  ```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. This scalar `child=` form uses one L1 patch.
+  For separated L1/L2 sets, use `child_rectangles_in_parent`; see
+  [coordinate frames and supported layouts](docs/fixed_three_level_multipatch.md).
+  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
+  ```
+- Manual synchronized L2 replacement: `replacement = current.remap_child(rectangles,
+  child_objects=original_objects)` preserves the absolute step count and live Jet
+  amplitudes. L0/L1 and whole body-containing L2 rectangles stay fixed; separated
+  fluid rectangles may change. The source remains independently owned. See
+  [ownership, pinning and transfer rules](docs/synchronized_leaf_remap.md).
+- Runnable three-level example: `PYTHONPATH=src python examples/multigrid_remap.py`.
+  It uses two separated L1 regions, an L2 stationary Circle/Jet and an explicit
+  fluid-rectangle move, retaining original root-frame body definitions. Add
+  `--check-geometry` for a CPU-only check or `--verify-jet` for an optional
+  next-step comparison with an explicitly driven control. No flowfiles are saved.
 - Not available: per-link surface exchange/force readout in blocked mode,
-  checkpoints, moving bodies, multiple patches or dynamic refinement, `Simulation`
-  integration. Validation so far is storage/operator equivalence and short runs,
-  not accuracy, convergence or long-time stability.
+  checkpoints, moving bodies, L1 bodies/Jets in a three-level hierarchy,
+  arbitrary depth, automatic AMR or `Simulation` integration.
+  Storage/operator equivalence, remap and short-run gates do not establish
+  convergence or long-time physical accuracy. The earlier shear-decay result is
+  a research accuracy limitation, not a demonstrated solver bug. Pinned-body
+  replay conditioned on actual L1 feeds is not an independent whole-body solver.
+  See [three-level validation scope](docs/fixed_three_level_multipatch.md) and
+  [remap validation scope](docs/synchronized_leaf_remap.md). The new usage example
+  has CPU geometry/syntax checks only until its GPU command is run on Horsea.
 
 ## Installation
 

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,185 @@
+# 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 periodic
+hierarchy, either body-free or with stationary Circles/Jets in unchanged pinned
+L2 rectangles (see below). 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 L0/L1 body collection; use `child_objects`/`child_jets` for pinned L2
+bodies. 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 undeclared bodies or 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
+periodic three-level runs whose existing L2 uses
+`child_rectangles_in_parent`, either body-free or under the pinned-body rules below.
+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. These are GPU execution gates,
+not CPU-mock results. The separate body-free 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 transfer/equivalence cases do not establish evolution accuracy. Historical
+pinned-body FP64 replay uses actual L1 feeds: it checks conditional L2 evolution,
+not an independent whole-body solver or whole-hierarchy accuracy. The body-free
+replay gate above does not establish that stronger claim either.
+
+## Small runnable usage example
+
+[`examples/multigrid_remap.py`](../examples/multigrid_remap.py) keeps the original
+root-frame Circle/Jet definitions, moves a fluid-only L2 rectangle, checks fluid
+leaves with geometry-derived masks and retires both old/new owners before their
+shared context. `--verify-jet` adds a temporary control solely to compare the
+next step against an explicitly set live amplitude; normal use needs no control.
+
+```sh
+# Host-only: no CUDA import or allocation.
+PYTHONPATH=src python examples/multigrid_remap.py --check-geometry
+# Horsea: select an approved idle physical GPU; it becomes logical device 0.
+PYTHONPATH="$PWD/src" CUDA_VISIBLE_DEVICES=0 python examples/multigrid_remap.py --verify-jet
+```
+
+The new example has passed host geometry and syntax/import checks only; its GPU
+run is pending Horsea. Expected scope: base counts 2 before/after remap, then 5;
+finite positive-density fluid leaves, NaN nonleaves, and exact next-step control
+agreement. This is a short public-API smoke run, not a physical accuracy test.
+No flowfiles are written.

examples/multigrid_remap.py

@@ -0,0 +1,126 @@
+"""Three levels, separated regions and an explicit fluid-region remap.
+
+Run on an approved idle GPU (logical device 0):
+    PYTHONPATH=src CUDA_VISIBLE_DEVICES=0 python examples/multigrid_remap.py --verify-jet
+Host-only geometry check, without CUDA:
+    PYTHONPATH=src python examples/multigrid_remap.py --check-geometry
+
+A short public-API smoke example, not an accuracy/stability benchmark. No files
+are written. Edit the rectangles below; this is manual remapping, not auto AMR.
+"""
+import argparse
+from contextlib import ExitStack
+
+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.config import LBMConfig
+from CelerisLab.lbm.equilibrium import compute_feq_d2q9
+from CelerisLab.lbm.refinement import FixedPatchSet2D
+from CelerisLab.lbm.synchronized_remap import RemapTopology
+
+CFG = LBMConfig(nx=48, ny=32, viscosity=.1)  # D2Q9/SRT/FP32, periodic
+PARENTS = ((4, 4, 20, 28), (24, 4, 44, 28))  # L0 cell-edge rectangles
+PIN = (16, 16, 32, 24)  # L1 allocated-grid indices; whole body region stays fixed
+BEFORE = (PIN, (56, 32, 72, 40))
+AFTER = (PIN, (60, 32, 76, 40))  # move only the separated fluid rectangle
+LIVE_AMPLITUDE = -.004
+
+
+def geometry():
+    # Retain these original L0 edge-frame definitions for every remap; never
+    # substitute the stepper's projected fine-local boundary objects.
+    objects = (SimObject(0, geometry=CircleGeometry(12., 10., .8),
+                         center=(12., 10.), radius=.8),)
+    jets = (Jet(Tube((10.5, 10.), 0., 3., .3), Parallel(0.), amplitude=.002),)
+    parent = FixedPatchSet2D.from_config(CFG, PARENTS)
+    meshes = tuple(RemapTopology.from_config(
+        CFG, parent, rectangles, child_objects=objects, child_jets=jets)
+        for rectangles in (BEFORE, AFTER))
+    if any(mesh.pinned_rectangles != (PIN,) for mesh in meshes):
+        raise RuntimeError("the body-containing rectangle must remain pinned")
+    return objects, jets, meshes
+
+
+def diagnostics(run, mesh):
+    """Masks come from geometry, never from whether downloaded values are finite."""
+    mass = 0.
+    for level, (g, leaf) in enumerate(zip(run.download(), mesh.leaf_masks)):
+        fluid = g[:, leaf].astype(np.float64)
+        rho = fluid.sum(axis=0)
+        if not np.isfinite(fluid).all() or np.any(rho <= 0):
+            raise RuntimeError(f"invalid L{level} fluid populations/density")
+        if not np.isnan(g[:, ~leaf]).all():
+            raise RuntimeError(f"unexpected L{level} nonleaf values")
+        mass += rho.sum() / 4**level  # area weights: 1, 1/4, 1/16
+    print(f"absolute base step {run.step_count}: finite fluid leaves, mass={mass:.9f}")
+    # Driven boundary injects/extracts mass; mass need not stay constant in time.
+
+
+def verify_live_jet(current, objects):
+    """Optional public-API check; normal continuation needs no control solver."""
+    # Identity remap copies the same state. Explicitly setting the known live
+    # amplitude on the control detects a stale amplitude on the continuation.
+    # No private fields or constructor-time Jet snapshots inspect the solver.
+    with current.remap_child(AFTER, child_objects=objects) as control:
+        control.set_jet_amplitude(0, LIVE_AMPLITUDE, level=2)
+        current.step()
+        control.step()
+        for actual, expected in zip(current.download(), control.download()):
+            if not np.array_equal(actual, expected, equal_nan=True):
+                raise RuntimeError("live Jet continuation differs from explicitly set control")
+    print("live Jet amplitude preserved: next step equals explicitly set control")
+
+
+def main():
+    parser = argparse.ArgumentParser(description=__doc__.splitlines()[0])
+    parser.add_argument("--check-geometry", action="store_true")
+    parser.add_argument("--verify-jet", action="store_true",
+                        help="compare one step with an explicitly driven control")
+    args = parser.parse_args()
+    objects, jets, (old_mesh, new_mesh) = geometry()
+    if args.check_geometry:
+        print("geometry OK: three levels, two L1 regions, fixed body pin, moved fluid region")
+        return
+
+    # Lazy imports keep --check-geometry useful on a CPU-only machine.
+    from CelerisLab.cuda.context import CudaContext
+    from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
+
+    # ExitStack owns every old/new stepper, including on failure. All are closed
+    # before the shared context. remap_child never consumes/closes its source.
+    with CudaContext(0) as context, ExitStack() as owned:
+        old = owned.enter_context(FixedPatchStepper(
+            CFG, PARENTS, context=context, base_block=(4, 4), fine_block=(4, 4),
+            child_block=(4, 4), child_rectangles_in_parent=BEFORE,
+            child_objects=objects, child_jets=jets))
+        rest = np.asarray(compute_feq_d2q9(1., 0., 0.), np.float32)[:, None, None]
+        old.initialize(*(np.broadcast_to(rest, (9, *mask.shape))
+                         for mask in old_mesh.leaf_masks))
+        old.step(2)  # one base step advances L1 twice and L2 four times
+        old.set_jet_amplitude(0, LIVE_AMPLITUDE, level=2)
+        diagnostics(old, old_mesh)
+        count = old.step_count
+        current = owned.enter_context(old.remap_child(AFTER, child_objects=objects))
+        if current is old or current.step_count != count:
+            raise RuntimeError("remap must return a distinct continuation at the same time")
+        diagnostics(current, new_mesh)
+
+        if args.verify_jet:
+            verify_live_jet(current, objects)  # optional control; advances one step
+        else:
+            current.step()
+        if old.step_count != count or current.step_count != count + 1:
+            raise RuntimeError("old/new absolute step ownership changed")
+        old.close()  # retire the old state only after successful replacement
+        current.step(2)
+        if current.step_count != count + 3:
+            raise RuntimeError("continuation lost its absolute step count")
+        diagnostics(current, new_mesh)
+        print("smoke example only: no accuracy or long-time stability claim")
+
+
+if __name__ == "__main__":
+    main()

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/lbm/refinement.py

@@ -114,6 +114,19 @@ class FixedPatch2D:
     def fine_origin(self) -> tuple[float, float]:
         return self.bounds[0] - HALO / 2, self.bounds[1] - HALO / 2
 
+    @property
+    def fine_level_origin(self) -> tuple[int, int]:
+        """(y, x) of this window's first cell on the level-wide 2:1 fine grid."""
+        return 2 * self.bounds[1] - HALO, 2 * self.bounds[0] - HALO
+
+    @property
+    def covered(self) -> np.ndarray:
+        """Dense coarse mask of cells replaced by the patch (not leaves)."""
+        x0, y0, x1, y1 = self.bounds
+        mask = np.zeros(self.coarse_shape, bool)
+        mask[y0:y1, x0:x1] = True
+        return mask
+
     @property
     def table_bytes(self) -> int:
         return sum(a.nbytes for a in (self.fine_kind, self.ghost_fine,
@@ -138,6 +151,113 @@ class FixedPatch2D:
             raise ValueError("levels are not synchronized: fine_steps must equal 2*coarse_steps")
 
 
+MIN_PATCH_GAP = 4  # coarse cells: the full 4-fine-cell supports stay disjoint
+
+
+@dataclass(frozen=True, eq=False)
+class FixedPatchSet2D:
+    """Several separated FixedPatch2D rectangles sharing one 2:1 fine level.
+
+    The fine level is the whole (2*ny, 2*nx) grid in level coordinates; only
+    the patch supports (active region + HALO ghost band) are stored, by a
+    sparse BlockLayout2D from fine_layout(). Tables have FixedPatch2D meaning
+    and row/child order, re-addressed to level cells: ghost_fine are level
+    fine cells, coalesce_fine are q*N+cell level slots; coarse ids unchanged.
+    Patches are kept in canonical (y0, x0) order, so input order is irrelevant.
+    Active regions must be >= MIN_PATCH_GAP coarse cells apart (Chebyshev).
+    """
+    coarse_shape: tuple[int, int]
+    patches: tuple
+    tau_coarse: float
+    tau_fine: float
+    fine_kind: np.ndarray
+    ghost_fine: np.ndarray
+    ghost_coarse: np.ndarray
+    coalesce_coarse: np.ndarray
+    coalesce_fine: np.ndarray
+
+    @classmethod
+    def from_config(cls, cfg: LBMConfig, rectangles) -> "FixedPatchSet2D":
+        rectangles = tuple(rectangles)
+        if not rectangles:
+            raise ValueError("at least one patch rectangle is required")
+        patches = sorted((FixedPatch2D.from_config(cfg, r) for r in rectangles),
+                         key=lambda p: (p.bounds[1], p.bounds[0]))
+        for i, a in enumerate(patches):
+            for b in patches[i + 1:]:
+                (ax0, ay0, ax1, ay1), (bx0, by0, bx1, by1) = a.bounds, b.bounds
+                if max(ay0 - by1, by0 - ay1, ax0 - bx1, bx0 - ax1) < MIN_PATCH_GAP:
+                    raise ValueError(f"patches {a.bounds} and {b.bounds} are closer "
+                                     f"than {MIN_PATCH_GAP} coarse cells")
+        ny, nx = patches[0].coarse_shape
+        fy, fx = 2 * ny, 2 * nx
+        if fy * fx * 9 > np.iinfo(np.uint64).max // 2:
+            raise ValueError("fine level exceeds 64-bit indexing")
+        kind = np.full((fy, fx), OUTER_GHOST, np.uint8)
+        ghost_fine, coalesce_fine = [], []
+        for p in patches:
+            y0, x0 = p.fine_level_origin
+            ly, lx = p.fine_shape
+            kind[y0:y0 + ly, x0:x0 + lx] = p.fine_kind
+            gy, gx = np.divmod(p.ghost_fine.astype(np.int64), lx)
+            ghost_fine.append((gy + y0) * fx + gx + x0)
+            q, cell = np.divmod(p.coalesce_fine.astype(np.int64), ly * lx)
+            cy, cx = np.divmod(cell, lx)
+            coalesce_fine.append(q * fy * fx + (cy + y0) * fx + cx + x0)
+        ghost_coarse = np.concatenate([p.ghost_coarse for p in patches])
+        coalesce_coarse = np.concatenate([p.coalesce_coarse for p in patches])
+        covered = np.zeros((ny, nx), bool)
+        for p in patches:
+            covered |= p.covered
+        cells = ny * nx
+        if (covered.ravel()[ghost_coarse].any()
+                or covered.ravel()[coalesce_coarse % np.uint64(cells)].any()
+                or len(np.unique(coalesce_coarse)) != len(coalesce_coarse)):
+            raise ValueError("patch transfers overlap another patch")
+        return cls((ny, nx), tuple(patches), patches[0].tau_coarse, patches[0].tau_fine,
+                   _readonly(kind, np.uint8), _readonly(np.concatenate(ghost_fine), np.uint64),
+                   _readonly(ghost_coarse, np.uint32), _readonly(coalesce_coarse, np.uint64),
+                   _readonly(np.concatenate(coalesce_fine), np.uint64).reshape(-1, 4))
+
+    @property
+    def fine_shape(self) -> tuple[int, int]:
+        return self.fine_kind.shape
+
+    @property
+    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)
+        for p in self.patches:
+            mask |= p.covered
+        return mask
+
+    def fine_layout(self, block):
+        """Sparse non-periodic layout storing exactly the patch supports."""
+        from .blocks import BlockLayout2D
+        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
+            ly, lx = p.fine_shape
+            if y0 % by or x0 % bx or ly % by or lx % bx:
+                raise ValueError(f"patch {p.bounds} support is not aligned to fine blocks {block}")
+            origins.update((y, x) for y in range(y0, y0 + ly, by) for x in range(x0, x0 + lx, bx))
+        return BlockLayout2D.build(*self.fine_shape, by, bx, periodic=(False, False),
+                                   origins=sorted(origins))
+
+
 @dataclass(frozen=True, eq=False)
 class FineBoundary2D:
     """Fixed local projection plus auditable, finite transfer dependency sets.
@@ -184,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")
@@ -215,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"))
@@ -224,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)
@@ -286,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.
@@ -294,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_stepper.py

@@ -3,7 +3,8 @@
 
 The caller owns the CudaContext and must close this backend before that context.
 Only synchronized postcollision_g leaf populations are public. No checkpoint,
-snapshot, forcing, runtime parameter mutation or dynamic refinement is offered.
+forcing or runtime parameter mutation is offered. Bounded manual periodic
+L2 set replacement is available at completed base intervals via remap_child.
 """
 from dataclasses import dataclass, replace
 from numbers import Integral
@@ -18,7 +19,8 @@ from .blocks import BlockLayout2D, blocked_link_tables
 from .curved_links import surface_link_velocity
 from ..body.jets import jet_link_shapes, jet_link_velocity
 from .block_stepper import MAX_GRID_DIM
-from .refinement import ACTIVE, HALO, OUTER_GHOST, FixedPatch2D, prepare_fine_boundary
+from .refinement import (ACTIVE, HALO, OUTER_GHOST, FixedPatch2D, FixedPatchSet2D,
+                         nest_patch, nest_patch_set, prepare_fine_boundary)
 from .descriptors import SOLID, D2Q9_EX, D2Q9_EY
 
 
@@ -32,31 +34,62 @@ def blocked_fine_tables(patch, layout, order):
     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 = np.full(layout.nblocks * int(np.prod(layout.local_shape)), OUTER_GHOST, np.uint8)
-    kind[owner] = patch.fine_kind.ravel()
+    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_base_tables(patch, layout, order, wall_aware=False):
+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
+
+
+@dataclass
+class _Pair:
+    """Device transfer tables of one nested parent->child relationship."""
+    explode: object      # BlockExplode from the CHILD module (child TOTAL_CELLS)
+    coalesce: object
+    ghosts: object
+    donors: object
+    targets: object      # receiver-local staged parent slots
+    children: object
+    n_ghosts: int
+    n_rows: int
+    parent_cells: int    # parent block window cells (explode donor stride)
+
+
+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[receiver.astype(np.intp)].astype(np.int64) - offset
-    kind = np.full(layout.nblocks * int(np.prod(layout.local_shape)), OUTER_GHOST, np.uint8)
-    leaves = np.ones(patch.coarse_shape, bool)
-    x0, y0, x1, y1 = patch.bounds
-    leaves[y0:y1, x0:x1] = False
-    if wall_aware:
-        leaves[[0, -1]] = False
-    kind[owner[leaves.ravel()].astype(np.intp)] = ACTIVE
-    return (owner[patch.ghost_coarse].astype(np.uint64), layout.slot(source, q), kind)
+    source = owner[receivers].astype(np.int64) - offset
+    return (owner[donors].astype(np.uint64), layout.slot(source, q),
+            blocked_kind(layout, dense_kind))
 
 
 @dataclass
@@ -72,6 +105,36 @@ class _Level:
     halo_copy: object = None
     halo_dst: object = None
     halo_src: object = None
+    boundary: object = None
+
+
+@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
+
+
+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)
+    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"], _ = 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)
 
 
 class FixedPatchStepper:
@@ -110,23 +173,65 @@ class FixedPatchStepper:
     ``jets`` (body.jets.Jet) use the same edge frame as ``fine_objects`` and
     drive fine cut links selected at their true hits; ``set_jet_amplitude``
     re-uploads the per-link velocity column only.
+
+    Several patches: pass ``bounds`` as a sequence of (x0, y0, x1, y1)
+    rectangles (FixedPatchSet2D, >= 4 coarse cells apart). This requires
+    fine_block and base_block; fine_objects/jets are projected once onto the
+    level-wide fine grid (each circle inside one patch, body ids = object
+    positions, jet and writer checks level-wide), and
+    uses one level-wide fine grid (2*ny, 2*nx): public fine arrays have that
+    shape, NaN outside the patches' ACTIVE cells. The step sequence is the
+    blocked one above, with every patch's transfers in one launch each.
+
+    Three-level separated unions: ``child_rectangles_in_parent=[...]`` is
+    mutually exclusive with scalar ``child`` and accepts either parent form.
+    Coordinates index the allocated L1 array, including ghost bands and holes:
+    root edge = patch.fine_origin + L1 index / 2. A parent set has origin zero;
+    a scalar parent retains its shifted window origin. The child set's public
+    extent is twice the allocated L1 shape, even though storage is sparse.
+    Each pair requires gap >= 4 in its own parent-cell units. All regions share
+    one buffer pair per level and the fixed 1:2:4 schedule; L2-only bodies/Jets
+    below apply unchanged, with object IDs assigned over the whole level.
+
+    Three-level single chain: ``child`` bounds use L1 window cells; all three
+    levels must be blocked. ``child_objects`` and ``child_jets`` use the L0
+    edge frame, with L2 dx=dt=.25 and unchanged lattice velocity. L1 bodies
+    and Jets remain unsupported in a chain. Every whole circle and its BFL
+    stencil must fit strictly away from the L2 transfer footprint. L2 runs
+    halo -> BFL -> step four times per base step. Runtime child Jet changes
+    require ``set_jet_amplitude(index, amplitude, level=2)``.
     """
     def __init__(self, cfg, bounds, *, context, coarse_boundary="periodic", fine_objects=(),
-                 fine_block=None, base_block=None, jets=()):
+                 fine_block=None, base_block=None, jets=(), child=None, child_block=None,
+                 child_objects=(), child_jets=(), child_rectangles_in_parent=None):
         if HALO != 4:
             raise RuntimeError("refinement.cu requires a four-cell fine halo")
         if coarse_boundary not in ("periodic", "slot_reservoir"):
             raise ValueError("coarse_boundary must be periodic or slot_reservoir")
-        self.patch = FixedPatch2D.from_config(cfg, bounds)
+        if child is not None and child_rectangles_in_parent is not None:
+            raise ValueError("child and child_rectangles_in_parent are mutually exclusive")
+        nested = child is not None or child_rectangles_in_parent is not None
+        bounds = tuple(bounds)
+        multipatch = len(bounds) > 0 and not isinstance(bounds[0], (Integral, np.integer))
+        if multipatch:
+            if fine_block is None or base_block is None:
+                raise ValueError("several patches require fine_block and base_block")
+            self.patch = FixedPatchSet2D.from_config(cfg, bounds)
+        else:
+            self.patch = FixedPatch2D.from_config(cfg, bounds)
         self._wall_aware = coarse_boundary == "slot_reservoir"
         if self._wall_aware:
             if (cfg.y_wall_bc != "on_site_velocity" or cfg.top_wall_bc != "reservoir"
                     or cfg.velocity != 0 or not 2 <= cfg.wall_slot_width < cfg.nx):
                 raise ValueError("slot_reservoir requires on_site_velocity, reservoir, "
                                  "velocity=0 and 2 <= wall_slot_width < nx")
-            if self.patch.bounds[1] < 4 or self.patch.bounds[3] > cfg.ny - 4:
+            rects = self.patch.bounds if multipatch else (self.patch.bounds,)
+            if any(r[1] < 4 or r[3] > cfg.ny - 4 for r in rects):
                 raise ValueError("slot_reservoir requires y0 >= 4 and y1 <= ny-4")
-        fine_objects = tuple(fine_objects)
+        fine_objects, jets = tuple(fine_objects), tuple(jets)
+        child_objects, child_jets = tuple(child_objects), tuple(child_jets)
+        if not nested and (child_objects or child_jets):
+            raise ValueError("child_objects/child_jets require a nested child")
         if base_block is not None and fine_block is None:
             raise ValueError("base_block requires fine_block")
         self.base_layout = None
@@ -143,21 +248,70 @@ class FixedPatchStepper:
             if len(fine_block) != 2:
                 raise ValueError("fine_block must be (by, bx)")
             fy, fx = self.patch.fine_shape
-            self.fine_layout = BlockLayout2D.build(fy, fx, *fine_block, periodic=(False, False))
+            self.fine_layout = (self.patch.fine_layout(fine_block) if multipatch else
+                                BlockLayout2D.build(fy, fx, *fine_block, periodic=(False, False)))
             if (self.fine_layout.nblocks > MAX_GRID_DIM
                     or self.fine_layout.local_shape[0] > MAX_GRID_DIM):
                 raise ValueError("fine blocks exceed launch limits")
+        self.child_patch = self.child_layout = None
+        self._middle_kind = None
+        if nested:
+            if child is not None and multipatch:
+                raise ValueError("scalar child requires one parent patch; use child_rectangles_in_parent")
+            if fine_block is None or base_block is None or child_block is None:
+                raise ValueError("a nested child needs fine_block, base_block and child_block")
+            if tuple(fine_objects) or tuple(jets):
+                raise ValueError("a nested child supports bodies/Jets only on L2; use child_objects/child_jets")
+            if len(child_block) != 2:
+                raise ValueError("child_block must be (by, bx)")
+            if child_rectangles_in_parent is not None:
+                self.child_patch, self._middle_kind = nest_patch_set(
+                    cfg, self.patch, child_rectangles_in_parent)
+                self.child_layout = self.child_patch.fine_layout(child_block)
+            else:
+                self.child_patch, self._middle_kind = nest_patch(cfg, self.patch, child)
+                self.child_layout = BlockLayout2D.build(*self.child_patch.fine_shape, *child_block,
+                                                       periodic=(False, False))
+            if (self.child_layout.nblocks > MAX_GRID_DIM
+                    or self.child_layout.local_shape[0] > MAX_GRID_DIM):
+                raise ValueError("child blocks exceed launch limits")
         self._jets = tuple(jets)
         self.fine_boundary = prepare_fine_boundary(
             cfg, self.patch, fine_objects, wall_aware=self._wall_aware, jets=self._jets)
-        self._fine_tables, self._actuated, self._jet_cache = None, False, None
-        if self.fine_layout is not None and self.fine_boundary is not None:
-            fluid_idx, dirs, q, _, _, _, _, fallback = self.fine_boundary.curved_host
-            self._fine_tables, _ = blocked_link_tables(
-                self.fine_layout, self.fine_boundary.flags, fluid_idx, dirs, q, fallback)
-        self._fine_leaf = self.patch.fine_kind == ACTIVE
+        fine_cfg = replace(cfg, nx=self.patch.fine_shape[1], ny=self.patch.fine_shape[0],
+                           viscosity=2 * cfg.viscosity)
+        if self._wall_aware:
+            fine_cfg = replace(fine_cfg, y_wall_bc="bounce_back", top_wall_bc="velocity",
+                               wall_slot_width=0, wall_normal_velocity=0.0, velocity=0.0)
+        self.child_boundary = (None if self.child_patch is None else prepare_fine_boundary(
+            fine_cfg, self.child_patch, child_objects, jets=child_jets,
+            parent_origin=self.patch.fine_origin, parent_spacing=.5))
+        child_cfg = (None if self.child_patch is None else replace(
+            fine_cfg, nx=self.child_patch.fine_shape[1], ny=self.child_patch.fine_shape[0],
+            viscosity=4 * cfg.viscosity))
+        self._boundaries = {
+            1: _prepare_blocked_boundary(self.fine_boundary, self.fine_layout, fine_cfg, self._jets),
+            2: _prepare_blocked_boundary(self.child_boundary, self.child_layout, child_cfg, child_jets),
+        }
+        # Compatibility aliases for the existing L1 inspection API.
+        b = self._boundaries[1]
+        self._fine_tables = None if b is None else b.tables
+        self._actuated = False if b is None else b.actuated
+        self._jet_cache = None if b is None else b.jet_cache
+        self._fine_leaf = (self.patch.fine_kind if self._middle_kind is None
+                           else self._middle_kind) == ACTIVE
+        self._child_leaf = (None if self.child_patch is None
+                            else self.child_patch.fine_kind == ACTIVE)
+        if self.child_boundary is not None:
+            self._child_leaf &= (self.child_boundary.flags.reshape(
+                self.child_patch.fine_shape) & SOLID) == 0
+        self._child = None
         if self.fine_boundary is not None:
             self._fine_leaf &= (self.fine_boundary.flags.reshape(self.patch.fine_shape) & SOLID) == 0
+        # Private constructor inputs for the bounded manual L2 continuation.
+        self._remap_cfg = replace(cfg)
+        self._remap_blocks = tuple(None if b is None else tuple(b)
+                                   for b in (base_block, fine_block, child_block))
         self.ctx = context
         self._lock = RLock()
         self._status = "new"
@@ -169,30 +323,15 @@ class FixedPatchStepper:
         self._coalesce = None
         self._link_tables = {}
         self._block = (cfg.threads_per_block, 1, 1)
-        fine_cfg = replace(cfg, nx=self.patch.fine_shape[1], ny=self.patch.fine_shape[0],
-                           viscosity=2 * cfg.viscosity)
-        if self._wall_aware:
-            fine_cfg = replace(fine_cfg, y_wall_bc="bounce_back", top_wall_bc="velocity",
-                               wall_slot_width=0, wall_normal_velocity=0.0, velocity=0.0)
-        if self._fine_tables is not None:
-            # Any actuation spec (even zero amplitude) requires the velocity buffer.
-            b = self.fine_boundary
-            self._actuated = bool(self._jets) or any(
-                obj.surface_actuation is not None for obj in b.objects)
-            if self._actuated:
-                fine_cfg.validate()  # derives nq etc.; the unique-writer key needs them
-                fluid_idx, dirs, _, body_id, rx, ry, _, _ = b.curved_host
-                self._fine_tables["surface_velocity"], _ = surface_link_velocity(
-                    b.objects, fine_cfg, b.flags, fluid_idx, dirs, body_id, rx, ry,
-                    jets=self._jets, hits=b.curved_hits)
-                if self._jets:
-                    self._jet_cache = jet_link_shapes(self._jets, b.curved_hits)
         with self.ctx.activate():
             try:
                 self._stream = cuda.Stream()
-                for index, (config, layout, tau, name) in enumerate((
-                        (replace(cfg), self.base_layout, self.patch.tau_coarse, "PatchCoarseStep"),
-                        (fine_cfg, self.fine_layout, self.patch.tau_fine, "PatchFineStep"))):
+                specs = [(replace(cfg), self.base_layout, self.patch.tau_coarse, "PatchCoarseStep"),
+                         (fine_cfg, self.fine_layout, self.patch.tau_fine, "PatchFineStep")]
+                if self.child_patch is not None:
+                    specs.append((child_cfg, self.child_layout,
+                                  self.child_patch.tau_fine, "PatchFineStep"))
+                for index, (config, layout, tau, name) in enumerate(specs):
                     source = "refinement.cu"
                     if layout is not None:
                         config = replace(config, nx=layout.local_shape[1], ny=layout.local_shape[0],
@@ -215,12 +354,14 @@ class FixedPatchStepper:
                     if size != struct.calcsize("=ffffI"):
                         raise RuntimeError("LBMParams host/device layout mismatch")
                     cuda.memcpy_htod(ptr, struct.pack("=ffffI", config.omega, 0, 0, 0, n_objects))
-                    shape = (9, *(self.patch.coarse_shape if index == 0 else self.patch.fine_shape))
+                    shape = (9, *(self.patch.coarse_shape, self.patch.fine_shape,
+                                  self.child_patch and self.child_patch.fine_shape)[index])
                     size = int(np.prod(shape if layout is None else layout.blocked_shape)) * 4
                     grid = ((config.nx + self._block[0] - 1) // self._block[0], config.ny,
                             1 if layout is None else layout.nblocks)
                     level = _Level(shape, module, self._allocate(size), self._allocate(size),
                                    module.get_function(name), grid, layout)
+                    level.boundary = self._boundaries.get(index)
                     self._levels.append(level)
                     if layout is not None:
                         level.halo_copy = module.get_function("BlockHaloCopy")
@@ -234,7 +375,7 @@ class FixedPatchStepper:
                             cuda.memcpy_htod(ptr, struct.pack("=IIIf", cfg.nx, cfg.ny,
                                                             cfg.wall_slot_width, cfg.wall_normal_velocity))
                 p = self.patch
-                c, f = self._levels
+                c, f = self._levels[:2]
                 order = np.argsort(p.coalesce_coarse)
                 fine_tables = (None if f.layout is None else
                                blocked_fine_tables(p, f.layout, order))
@@ -250,15 +391,33 @@ class FixedPatchStepper:
                 self._explode = f.module.get_function("PatchExplode" if fine_tables is None
                                                        else "BlockExplode")
                 if fine_tables is not None:
-                    f.kind = self._upload(fine_tables[2])
+                    f.kind = self._upload(fine_tables[2] if self._middle_kind is None
+                                          else blocked_kind(f.layout, self._middle_kind))
                 if base_tables is not None:
                     c.kind = self._upload(base_tables[2])
                     self._coalesce = c.module.get_function("BlockCoalesce")
-                if self._fine_tables is not None:
-                    self._fine_flags = self._upload(self._fine_tables["flags"])
-                    self._link_tables = {k: self._upload(v) for k, v in self._fine_tables.items()
-                                         if k != "flags"}
-                    self._curved = self._levels[1].module.get_function("BlockCurvedBoundary")
+                if self.child_patch is not None:
+                    g, cp = self._levels[2], self.child_patch
+                    order = np.argsort(cp.coalesce_coarse)
+                    ghosts, children, g_kind = blocked_fine_tables(cp, g.layout, order)
+                    donors, targets, _ = blocked_base_tables(
+                        cp, f.layout, order, composed_kind=self._middle_kind)
+                    g.kind = self._upload(g_kind)
+                    g.halo_copy = g.module.get_function("BlockHaloCopy")
+                    self._child = _Pair(g.module.get_function("BlockExplode"),
+                                        f.module.get_function("BlockCoalesce"),
+                                        self._upload(ghosts), self._upload(donors.astype(np.uint64)),
+                                        self._upload(targets), self._upload(children),
+                                        len(ghosts), len(targets), int(np.prod(f.layout.local_shape)))
+                for level in self._levels:
+                    b = level.boundary
+                    if b is not None:
+                        b.flags = self._upload(b.tables["flags"])
+                        b.device = {k: self._upload(v) for k, v in b.tables.items() if k != "flags"}
+                        b.kernel = level.module.get_function("BlockCurvedBoundary")
+                if self._boundaries[1] is not None:
+                    b = self._boundaries[1]
+                    self._fine_flags, self._link_tables, self._curved = b.flags, b.device, b.kernel
                 elif self.fine_boundary is not None:
                     from .curved_links import CurvedLinkSoA
                     boundary = self.fine_boundary
@@ -296,27 +455,28 @@ class FixedPatchStepper:
     def step_count(self):
         return self._step_count
 
-    def _mask_nonleaves(self, coarse, fine):
-        x0, y0, x1, y1 = self.patch.bounds
-        coarse[:, y0:y1, x0:x1] = np.nan
+    def _mask_nonleaves(self, coarse, fine, child=None):
+        coarse[:, self.patch.covered] = np.nan
         if self._wall_aware:
             coarse[:, (0, -1), :] = np.nan
         fine[:, ~self._fine_leaf] = np.nan
+        if child is not None:
+            child[:, ~self._child_leaf] = np.nan
 
-    def initialize(self, coarse_g, fine_g):
-        """Upload an explicit G_0; invalid non-leaves are discarded."""
+    def initialize(self, coarse_g, fine_g, child_g=None):
+        """Upload an explicit G_0 (child_g only with a nested child); non-leaves are discarded."""
         with self._lock, self.ctx.activate():
             self._require("new", "ready")
-            arrays = [np.array(a, dtype=np.float32, order="C", copy=True)
-                      for a in (coarse_g, fine_g)]
+            given = (coarse_g, fine_g) + (() if child_g is None else (child_g,))
+            if len(given) != len(self._levels):
+                raise ValueError("Expected one Q-major array per level")
+            arrays = [np.array(a, dtype=np.float32, order="C", copy=True) for a in given]
             if any(a.shape != level.shape for a, level in zip(arrays, self._levels)):
-                raise ValueError("Expected Q-major (9,ny,nx) arrays for both levels")
-            x0, y0, x1, y1 = self.patch.bounds
-            leaf = np.ones(self.patch.coarse_shape, dtype=bool)
-            leaf[y0:y1, x0:x1] = False
+                raise ValueError("Expected Q-major (9,ny,nx) arrays for every level")
+            leaf = ~self.patch.covered
             if self._wall_aware:
                 leaf[[0, -1]] = False
-            for array, mask in zip(arrays, (leaf, self._fine_leaf)):
+            for array, mask in zip(arrays, (leaf, self._fine_leaf, self._child_leaf)):
                 if not np.isfinite(array[:, mask]).all() or (array[:, mask].sum(axis=0) <= 0).any():
                     raise ValueError("Leaf populations must be finite with positive density")
             self._mask_nonleaves(*arrays)
@@ -334,6 +494,86 @@ class FixedPatchStepper:
                 self._status = "failed"
                 raise
 
+    def remap_child(self, rectangles, *, child_objects=()):
+        """Return a NEW ready stepper with an L2 set replaced at this base time.
+
+        Periodic three-level set children only; L0/L1 and block sizes stay fixed.
+        A body run requires its original root-edge-frame ``child_objects``:
+        stationary circles pinned in unchanged complete L2 rectangles.
+        Separated fluid-only rectangles may be added, moved or removed.
+        Current L2 Jets (including updated amplitudes)
+        are preserved. Local projected boundary objects are not valid inputs.
+        The lock and stream synchronization certify
+        current postcollision G and the completed 1:2:4 counts internally.
+        Both objects remain independently owned: the caller must close each,
+        and chooses when to stop using the unchanged old instance. Candidate
+        construction/upload errors never invalidate the old instance; failure
+        to synchronize/read its device state does, as with step/download.
+        """
+        from .synchronized_remap import (RemapTopology, remap_synchronized,
+                                         same_physical_boundary, boundary_velocity)
+
+        with self._lock:
+            self._require("ready")
+            if (self._wall_aware or len(self._levels) != 3
+                    or not isinstance(self.child_patch, FixedPatchSet2D)
+                    or self.fine_boundary is not None
+                    or any(level.boundary is not None for level in self._levels[:2])):
+                raise ValueError("remap_child requires periodic three-level L2 sets with no L0/L1 bodies")
+            cfg = self._remap_cfg
+            child_objects = tuple(child_objects)
+            actual = self._boundaries[2]
+            jets = () if actual is None else actual.jets
+            old = RemapTopology.from_config(cfg, self.patch, self.child_patch.bounds,
+                                           child_objects=child_objects, child_jets=jets)
+            if (not same_physical_boundary(old.boundary, self.child_boundary)
+                    or (self.child_boundary is not None and actual is None)):
+                raise ValueError("child_objects must reproduce the original root-edge-frame boundary")
+            if actual is not None:
+                velocity = actual.tables.get("surface_velocity")
+                if velocity is None:
+                    velocity = boundary_velocity(cfg, self.child_patch, actual.host, jets)
+                if (not same_physical_boundary(old.boundary, actual.host)
+                        or not np.array_equal(old.surface_velocity, velocity)):
+                    raise ValueError("child_objects changed boundary ownership or surface velocity")
+            new = RemapTopology.from_config(cfg, self.patch, rectangles,
+                                           child_objects=child_objects, child_jets=jets)
+            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 boundary")
+            base_block, fine_block, child_block = self._remap_blocks
+            # Check sparse block alignment/launch size before any device read/allocation.
+            layout = new.child.fine_layout(child_block)
+            if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM:
+                raise ValueError("child blocks exceed launch limits")
+            with self.ctx.activate():
+                try:
+                    self._stream.synchronize()
+                    source = self.download()
+                except BaseException:
+                    self._status = "failed"
+                    raise
+            count = self._step_count
+            transferred = remap_synchronized(
+                old, new, source, base_count=count, level_steps=(count, 2*count, 4*count),
+                phase="postcollision_g")
+            candidate = type(self)(
+                cfg, self.patch.bounds, context=self.ctx, base_block=base_block,
+                fine_block=fine_block, child_block=child_block,
+                child_rectangles_in_parent=new.child.bounds,
+                child_objects=child_objects, child_jets=jets)
+            try:
+                candidate.initialize(*transferred.g)
+                candidate._step_count = count
+            except BaseException:
+                try:
+                    candidate.close()
+                except Exception:
+                    pass  # Preserve the original candidate initialization failure.
+                raise
+            return candidate
+
     def _apply_fine_boundary(self, buffer, **kw):
         """Stage solid source slots; exchange buffer holds this fine substep only."""
         cl = self._fine_curved
@@ -354,48 +594,72 @@ class FixedPatchStepper:
                      cl.link_exchange_gpu if cl.link_exchange_gpu is not None else np.uintp(0),
                      grid=((cl.count + self._block[0] - 1) // self._block[0], 1), **kw)
 
-    def set_jet_amplitude(self, index, amplitude):
+    def set_jet_amplitude(self, index, amplitude, *, level=1):
         """Replace one jet's amplitude between steps; no geometry/table rebuild.
 
+        ``level=1`` preserves the L1 default; use ``level=2`` for child_jets.
         Invalid input raises before any change. A device/sync error marks the
         instance failed (as step/initialize do); host state commits only after
         a successful upload.
         """
         with self._lock, self.ctx.activate():
             self._require("new", "ready")
-            jets = list(self._jets)
+            if (isinstance(level, bool) or not isinstance(level, Integral)
+                    or level not in (1, 2) or level >= len(self._levels)):
+                raise ValueError("Jet level must be 1 or an existing child level 2")
+            boundary = self._boundaries[level]
+            jets = list(self._jets if level == 1 else (() if boundary is None else boundary.jets))
             jets[index] = jets[index].with_amplitude(amplitude)  # validates first
             velocity = None
-            if self._fine_curved is None and self._jet_cache is not None:
-                shape, jet_id = self._jet_cache
-                velocity = self._fine_tables["surface_velocity"].copy()
+            if boundary is not None and boundary.jet_cache is not None:
+                shape, jet_id = boundary.jet_cache
+                velocity = boundary.tables["surface_velocity"].copy()
                 driven = jet_id >= 0
                 velocity[driven] = jet_link_velocity(jets, shape, jet_id)[driven]
             try:
-                if self._fine_curved is not None:
+                if level == 1 and self._fine_curved is not None:
                     self._fine_curved.set_jet_amplitude(index, amplitude, self._stream)
                 elif velocity is not None:
                     self._stream.synchronize()  # no in-flight step reads the column
-                    cuda.memcpy_htod(self._link_tables["surface_velocity"], velocity)
+                    cuda.memcpy_htod(boundary.device["surface_velocity"], velocity)
             except BaseException:
                 self._status = "failed"
                 raise
             if velocity is not None:
-                self._fine_tables["surface_velocity"] = velocity
-            self._jets = tuple(jets)
+                boundary.tables["surface_velocity"] = velocity
+            if boundary is not None:
+                boundary.jets = tuple(jets)
+            if level == 1:
+                self._jets = tuple(jets)
 
     def _blocked_fine_boundary(self, buffer, **kw):
-        if self._fine_tables is None:
+        self._blocked_boundary(self._levels[1], buffer, **kw)
+
+    def _blocked_boundary(self, level, buffer, **kw):
+        b = level.boundary
+        if b is None:
             return
-        if self._curved is None or self._fine_flags is None:
-            raise RuntimeError("Fine boundary device buffers are missing")
-        t, count = self._link_tables, len(self._fine_tables["cell"])
+        if b.kernel is None or b.flags is None:
+            raise RuntimeError("Boundary device buffers are missing")
+        t, count = b.device, len(b.tables["cell"])
         velocity = t.get("surface_velocity")
-        if self._actuated and velocity is None:
+        if b.actuated and velocity is None:
             raise RuntimeError("Surface velocity device buffer is missing")
-        self._curved(buffer, t["block"], t["cell"], t["dir"], t["q"], t["fallback"],
-                     np.uint32(count), velocity if velocity is not None else np.uintp(0),
-                     grid=((count + self._block[0] - 1) // self._block[0], 1), **kw)
+        b.kernel(buffer, t["block"], t["cell"], t["dir"], t["q"], t["fallback"],
+                 np.uint32(count), velocity if velocity is not None else np.uintp(0),
+                 grid=((count + self._block[0] - 1) // self._block[0], 1), **kw)
+
+    def _advance_child(self, parent_buffer, **kw):
+        """Explode from the given L1 buffer, then two L2 hops (halo -> BFL -> step)."""
+        pair, g = self._child, self._levels[2]
+        pair.explode(parent_buffer, g.current, pair.ghosts, pair.donors, np.uint32(pair.n_ghosts),
+                     np.uint64(pair.parent_cells),
+                     grid=((pair.n_ghosts + self._block[0] - 1) // self._block[0], 1), **kw)
+        for src, dst in ((g.current, g.spare), (g.spare, g.current)):
+            self._gather_halo(g, src, **kw)
+            self._blocked_boundary(g, src, **kw)
+            flags = np.uintp(0) if g.boundary is None else g.boundary.flags
+            g.kernel(src, dst, np.uint32(g.layout.halo), flags, g.kind, grid=g.grid, **kw)
 
     def _gather_halo(self, level, buffer, **kw):
         pairs = len(level.layout.halo_dst)
@@ -409,7 +673,7 @@ class FixedPatchStepper:
             raise ValueError("n must be a nonnegative integer")
         with self._lock, self.ctx.activate():
             self._require("ready")
-            c, f = self._levels
+            c, f = self._levels[:2]  # a nested child (third level) runs inside the fine hops
             p = self.patch
             kw = dict(block=self._block, stream=self._stream)
             try:
@@ -426,7 +690,17 @@ class FixedPatchStepper:
                         # Same buffer per hop: halo -> receiver-local BFL -> step.
                         flags = self._fine_flags if self._fine_flags is not None else np.uintp(0)
                         for src, dst in ((f.current, f.spare), (f.spare, f.current)):
+                            if self._child is not None:
+                                # L1 feeds L2 from its CURRENT G (fresh at the midpoint),
+                                # L2 takes two hops, then is restricted into src.
+                                self._advance_child(src, **kw)
                             self._gather_halo(f, src, **kw)
+                            if self._child is not None:
+                                pair = self._child
+                                pair.coalesce(self._levels[2].current, src, pair.targets,
+                                              pair.children, np.uint32(pair.n_rows),
+                                              grid=((pair.n_rows + self._block[0] - 1)
+                                                    // self._block[0], 1), **kw)
                             self._blocked_fine_boundary(src, **kw)
                             f.kernel(src, dst, np.uint32(self.fine_layout.halo), flags,
                                      f.kind, grid=f.grid, **kw)
@@ -462,7 +736,7 @@ class FixedPatchStepper:
                 raise
 
     def download(self):
-        """Return two synchronized G arrays, NaN on all non-leaves; no host cache."""
+        """Return synchronized G arrays (one per level), NaN on all non-leaves; no host cache."""
         with self._lock, self.ctx.activate():
             self._require("ready")
             arrays = [np.empty(level.shape if level.layout is None else level.layout.blocked_shape,
@@ -492,8 +766,11 @@ class FixedPatchStepper:
                         self._stream.synchronize()
                 except Exception as exc:
                     error = exc
+                for b in getattr(self, "_boundaries", {}).values():
+                    if b is not None:
+                        b.flags = b.device = b.kernel = None
                 self._levels.clear()
-                self._explode = self._coalesce = None
+                self._explode = self._coalesce = self._child = None
                 if getattr(self, "_fine_curved", None) is not None:
                     try:
                         self._fine_curved.free()

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

@@ -163,5 +163,221 @@ class FixedPatchTests(unittest.TestCase):
                 p.check_synchronization(c, f, phase=phase)
 
 
+class FixedPatchSetTests(unittest.TestCase):
+    """Several separated patches: level-wide re-addressing of FixedPatch2D tables."""
+    FACE = ((4, 4, 8, 12), (12, 4, 16, 12))          # 16x24 coarse, face gap 4
+    S1 = ((4, 4, 12, 12), (20, 4, 28, 12))           # 16x32 coarse, translation copy
+
+    def build(self, rects, ny=16, nx=24):
+        return refinement.FixedPatchSet2D.from_config(LBMConfig(nx=nx, ny=ny, viscosity=.1), rects)
+
+    def test_tables_are_shifted_single_patch_tables(self):
+        s = self.build(self.FACE)
+        fy, fx = s.fine_shape
+        self.assertEqual((fy, fx), (32, 48))
+        rows = cells = 0
+        for p in s.patches:
+            y0, x0 = p.fine_level_origin
+            ly, lx = p.fine_shape
+            np.testing.assert_array_equal(s.fine_kind[y0:y0 + ly, x0:x0 + lx], p.fine_kind)
+            gy, gx = np.divmod(p.ghost_fine.astype(np.int64), lx)
+            np.testing.assert_array_equal(s.ghost_fine[cells:cells + len(gy)],
+                                          (gy + y0) * fx + gx + x0)
+            np.testing.assert_array_equal(s.ghost_coarse[cells:cells + len(gy)], p.ghost_coarse)
+            q, c = np.divmod(p.coalesce_fine.astype(np.int64), ly * lx)
+            cy, cx = np.divmod(c, lx)
+            np.testing.assert_array_equal(s.coalesce_fine[rows:rows + len(q)],
+                                          q * fy * fx + (cy + y0) * fx + cx + x0)
+            np.testing.assert_array_equal(s.coalesce_coarse[rows:rows + len(q)], p.coalesce_coarse)
+            rows, cells = rows + len(q), cells + len(gy)
+        self.assertEqual((rows, cells), (len(s.coalesce_coarse), len(s.ghost_fine)))
+        self.assertEqual(int((s.fine_kind == ACTIVE).sum()), 4 * int(s.covered.sum()))
+
+    def test_input_order_is_irrelevant(self):
+        a, b = self.build(self.FACE), self.build(self.FACE[::-1])
+        for name in ("fine_kind", "ghost_fine", "ghost_coarse", "coalesce_coarse", "coalesce_fine"):
+            np.testing.assert_array_equal(getattr(a, name), getattr(b, name))
+        self.assertEqual(a.bounds, b.bounds)
+
+    def test_one_patch_set_matches_fixed_patch(self):
+        p = FixedPatch2D.from_config(config(), (4, 4, 12, 12))
+        s = self.build(((4, 4, 12, 12),), nx=16)
+        np.testing.assert_array_equal(s.covered, p.covered)
+        self.assertEqual((s.tau_coarse, s.tau_fine), (p.tau_coarse, p.tau_fine))
+        np.testing.assert_array_equal(s.coalesce_coarse, p.coalesce_coarse)
+        np.testing.assert_array_equal(s.ghost_coarse, p.ghost_coarse)
+
+    def test_gap_and_overlap_rejected(self):
+        for gap in (0, 2, 3):
+            with self.subTest(gap=gap), self.assertRaisesRegex(ValueError, "closer"):
+                self.build(((4, 4, 8, 12), (8 + gap, 4, 12 + gap, 12)))
+        with self.assertRaises(ValueError):
+            self.build(())
+        with self.assertRaises(ValueError):  # FixedPatch2D margin still applies
+            self.build(((1, 4, 8, 12), (12, 4, 16, 12)))
+
+    def test_sparse_fine_layout_stores_exactly_the_supports(self):
+        for rects, nx in ((self.FACE, 24), (self.S1, 32)):
+            s = self.build(rects, nx=nx)
+            layout = s.fine_layout((4, 4))
+            support = np.zeros(s.fine_shape, bool)
+            for p in s.patches:
+                y0, x0 = p.fine_level_origin
+                support[y0:y0 + p.fine_shape[0], x0:x0 + p.fine_shape[1]] = True
+            np.testing.assert_array_equal(layout.owned, support)
+            self.assertTrue(layout.owned.ravel()[s.ghost_fine.astype(np.intp)].all())
+        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__":
     unittest.main(verbosity=2)

tests/unit/test_refinement_three_level_multipatch.py

@@ -0,0 +1,314 @@
+"""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._explode.side_effect=record('explode01')
+            run._child.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._child.coalesce.side_effect=lambda src,dst,*a,**kw: events.append(('restrict21',dst))
+            run._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_child
+            def checked_advance(buffer, **kw):
+                self.assertIs(buffer,(f.current,f.spare)[len(feeds)])
+                feeds.append(buffer)
+                advance(buffer,**kw)
+            with patch.object(run,'_advance_child',side_effect=lambda buffer,**kw:
+                              checked_advance(f.current,**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)
+
+
+if __name__=='__main__':
+    unittest.main()

tests/unit/test_runtime_child_remap.py

@@ -0,0 +1,119 @@
+"""CPU mocks execute production remap lifecycle; no GPU numerical claim."""
+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)
+        for attribute, value in (('_wall_aware', True), ('child_boundary', object()),
+                                 ('fine_boundary', object()), ('child_patch', 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_refinement_multipatch.py

@@ -0,0 +1,226 @@
+"""Opt-in multipatch FixedPatchStepper (FixedPatchSet2D): CPU guards plus GPU gates.
+
+Run CPU and GPU classes in separate processes. GPU gates (frozen before any run):
+- one rectangle through the multipatch interface is bit-identical to the old
+  single-patch blocked path (both boundary modes, raw owner masks included);
+- 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, 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
+import unittest
+
+import numpy as np
+import pycuda.driver as cuda
+
+from CelerisLab.config import LBMConfig
+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
+
+FACE = ((4, 4, 8, 12), (12, 4, 16, 12))
+S1 = ((4, 4, 12, 12), (20, 4, 28, 12))
+
+
+def level_fine(patch_set, windows):
+    """Level-wide fine array with each patch window filled from `windows`."""
+    out = np.full((9, *patch_set.fine_shape), np.nan, np.float32)
+    for p, w in zip(patch_set.patches, windows):
+        y0, x0 = p.fine_level_origin
+        out[:, y0:y0 + p.fine_shape[0], x0:x0 + p.fine_shape[1]] = w
+    return out
+
+
+def window(patch_set, index, fine):
+    p = patch_set.patches[index]
+    y0, x0 = p.fine_level_origin
+    return fine[:, y0:y0 + p.fine_shape[0], x0:x0 + p.fine_shape[1]]
+
+
+class MultipatchCPU(unittest.TestCase):
+    def test_rejections_before_cuda(self):
+        cfg = config("periodic")
+        context = object()  # never activated: every case must fail first
+        for kw in (dict(), dict(fine_block=(8, 4)), dict(base_block=(4, 4)),
+                   dict(fine_block=(8, 4), base_block=(4, 4), fine_objects=[circle()]),
+                   dict(fine_block=(5, 4), base_block=(4, 4))):
+            with self.subTest(kw=kw), self.assertRaises(ValueError):
+                FixedPatchStepper(cfg, ((4, 4, 8, 12), (14, 4, 18, 12)), context=context, **kw)
+        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):
+    STEPS = (0, 1, 2, 8)
+
+    def run_steps(self, context, cfg, bounds, coarse, fine, **kw):
+        out = []
+        with FixedPatchStepper(cfg, bounds, context=context, **kw) as backend:
+            backend.initialize(coarse, fine)
+            for step in self.STEPS:
+                backend.step(step - backend.step_count)
+                out.append(backend.download() + raw(backend))
+        return out
+
+    def test_one_rectangle_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))
+                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, b[2], a[2], label + " raw coarse")
+                    exact(self, window(s, 0, b[3]), a[3], label + " raw fine")
+
+    def test_s1_translation_copy_and_order(self):
+        from CelerisLab.cuda.context import CudaContext
+        cfg16 = LBMConfig(nx=16, ny=16, viscosity=.1)
+        cfg32 = replace(cfg16, nx=32)
+        one = FixedPatchSet2D.from_config(cfg16, (S1[0],))
+        two = FixedPatchSet2D.from_config(cfg32, S1)
+        coarse, fine = state(16, 16, 3), state(*one.patches[0].fine_shape, 4)
+        kw = dict(fine_block=(4, 4), base_block=(4, 4))
+        with CudaContext(0) as context:
+            ref = self.run_steps(context, cfg16, (S1[0],), coarse, level_fine(one, [fine]), **kw)
+            tiled = np.concatenate([coarse, coarse], axis=2)
+            for bounds in (S1, S1[::-1]):
+                got = self.run_steps(context, cfg32, bounds, tiled,
+                                     level_fine(two, [fine, fine]), **kw)
+                for step, a, b in zip(self.STEPS, ref, got):
+                    label = f"order {bounds[0]} 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}")
+
+    def test_face_gap4_finite_masks_and_order(self):
+        from CelerisLab.cuda.context import CudaContext
+        cfg = LBMConfig(nx=24, ny=16, viscosity=.1)
+        s = FixedPatchSet2D.from_config(cfg, FACE)
+        coarse = state(16, 24, 5)
+        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])]
+        for step, a, b in zip(self.STEPS, *runs):
+            for x, y in zip(a, b):
+                exact(self, y, x, f"reversed order step {step}")
+            c, f = a[0], a[1]
+            self.assertTrue(np.isfinite(c[:, ~s.covered]).all())
+            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()