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

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

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

README.md

@@ -123,7 +123,7 @@ with CudaContext(0) as ctx, FixedPatchStepper(cfg, (8, 6, 24, 18), context=ctx,
 ```
 
 Scope and limits:
-- D2Q9, SRT, FP32, double buffer, unshifted, no LES or forcing; fixed 2:1 patches, two or three levels with separated rectangles.
+- D2Q9, SRT, FP32, double buffer, unshifted, no LES or forcing; fixed 2:1 patch(es), two levels (or one nested third).
 - `coarse_boundary="periodic"` (default) or `"slot_reservoir"` (bottom slot, top
   reservoir; patch needs `y0 >= 4`, `y1 <= ny-4`). Periodic mode ignores the wall/slot config fields.
 - Object centers/radii use the coarse **cell-edge** frame (cell `I` spans `[I, I+1)`),
@@ -144,8 +144,7 @@ Scope and limits:
   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.
+  different bodies, even across patches. No wall-touching patches or nesting yet.
 
   ```python
   cfg = LBMConfig(nx=32, ny=16, viscosity=.1)
@@ -165,10 +164,7 @@ Scope and limits:
   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
+  steps = 4 L2 steps. One L1 patch only. Optional stationary circles and Jets
   live on L2 via `child_objects` and `child_jets`; `fine_objects`/`jets` retain
   their L1 meaning and must be empty in this mode. All object and Jet geometry
   uses the **L0 edge frame**, not the child bounds' L1 window frame. L2 has
@@ -196,26 +192,14 @@ Scope and limits:
       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, 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.
+  checkpoints, moving bodies, L1 bodies/Jets in a three-level chain,
+  several patches per level with nesting, dynamic refinement, `Simulation`
+  integration. Validation so far is storage/operator equivalence and short runs,
+  not accuracy, convergence or long-time stability. The new L2-body GPU gates
+  (layouts, masks, amplitude reversal, legacy equivalence) and independent
+  every-substep numerical replay remain required acceptance checks; CPU tests
+  alone do not establish GPU numerical correctness.
 
 ## Installation
 

docs/development/raw_exchange.md

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

docs/synchronized_leaf_remap.md

@@ -1,9 +1,8 @@
 # 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,
+It changes only a separated L2 rectangle set in a three-level, body-free,
+periodic hierarchy. L0/L1 topology, relaxation parameters, direction order,
 physical frame and completed time remain unchanged. Existing geometry checks
 validate child support/nesting and separation. The CPU helper itself adds no CUDA, scheduler, allocator, boundary handling,
 collision or ghost refill. A bounded runtime adapter is described below.
@@ -30,12 +29,11 @@ origin is retained. Physical lower edges are `origin + spacing*(x,y)` on L0,
 `origin + spacing*(parent.fine_origin + (x,y)/4)` on L2.
 
 Factory options `coarse_boundary` and `objects` must declare `"periodic"` and
-an empty L0/L1 body collection; use `child_objects`/`child_jets` for pinned L2
-bodies. LBMConfig's channel settings do not
+an empty body collection (all levels). LBMConfig's channel settings do not
 select this operator's boundary mode. Geometry metadata is factory-created;
 manual mutation/forging of metadata is unsupported. Caller certifies the
-actual state has no undeclared bodies or forcing, that all device work and old
-readers are quiescent, and that inputs are the **current owned-leaf postcollision G**.
+actual state has no bodies/forcing, that all device work and old readers are
+quiescent, and that inputs are the **current owned-leaf postcollision G**.
 Count/phase stamps check consistency, not freshness or active-buffer identity.
 
 `current_g` is three NumPy Q-major arrays `(9,ny,nx)`. For blocked inputs,
@@ -75,9 +73,8 @@ PYTHONPATH=src python -m unittest discover -s tests/unit -p 'test_synchronized_r
 
 `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
+body-free periodic three-level runs whose existing L2 uses
+`child_rectangles_in_parent`. Scalar or union L1 parents are supported; scalar
 `child=` is deliberately rejected. L0/L1, relaxation and block sizes stay fixed.
 No tagger, scheduler, in-place replacement, or general rollback is provided.
 
@@ -108,9 +105,9 @@ 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.
+finite leaves, exact masks and absolute counts. They were not run on the CPU
+host. The separate independent FP64 capture/replay gate is max population
+error **3e-6**, with exact masks; remap mass/momentum allow FP32 roundoff.
 Per-population conservation is not an advancement target under collision.
 
 ### Bounded pinned-body extension
@@ -158,28 +155,6 @@ Multiple-body CPU coverage: `tests/unit/test_multiple_pinned_remap.py`.
 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.
+These new GPU cases have **not been run here**. They do not establish evolution
+accuracy. A 3e-6 FP32-vs-FP64 evolution gate requires a separate independent
+BFL reference; the body-free replay gate above does not establish that claim.

examples/multigrid_remap.py

@@ -1,126 +0,0 @@
-"""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()

src/CelerisLab/body/jets.py

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

src/CelerisLab/lbm/block_stepper.py

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

src/CelerisLab/lbm/hierarchy.py

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

src/CelerisLab/lbm/kernels/blocks.cu

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

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

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

src/CelerisLab/lbm/raw_exchange.py

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

src/CelerisLab/lbm/refinement_tables.py

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

tests/__init__.py

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

tests/unit/test_hierarchy.py

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

tests/unit/test_jet_condition.py

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

tests/unit/test_raw_exchange.py

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

tests/unit/test_refinement_three_level_multipatch.py

@@ -251,14 +251,14 @@ class ConstructorAndSchedule(unittest.TestCase):
             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')
+            run._pairs[0].explode.side_effect=record('explode01')
+            run._pairs[1].explode.side_effect=record('explode12')
             for i,level in enumerate(run._levels):
                 level.kernel.side_effect=record(f'step{i}')
             if bodies:
                 g.boundary.kernel.side_effect=record('bfl2')
-            run._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._pairs[1].coalesce.side_effect=lambda src,dst,*a,**kw: events.append(('restrict21',dst))
+            run._pairs[0].coalesce.side_effect=lambda src,dst,*a,**kw: events.append(('restrict10',dst))
             run._status='ready'
             expected=[('explode01',initial_c)]
             for middle in (f.current,f.spare):
@@ -278,13 +278,14 @@ class ConstructorAndSchedule(unittest.TestCase):
             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)):
+            advance=run._advance
+            def checked_advance(index, buffer, **kw):
+                if index == 2:  # L1 -> L2 feed: must be the L1 hop source
+                    self.assertIs(buffer,(f.current,f.spare)[len(feeds)])
+                    feeds.append(buffer)
+                advance(index,buffer,**kw)
+            with patch.object(run,'_advance',side_effect=lambda index,buffer,**kw:
+                              checked_advance(index,f.current if index == 2 else buffer,**kw)):
                 with self.assertRaises(AssertionError):
                     run.step()
             # A per-region shared-parent advance must fail the 1:2:4 guard.
@@ -310,5 +311,24 @@ class ConstructorAndSchedule(unittest.TestCase):
                 mock_backend(**options)
 
 
+
+
+class InitialLevels(unittest.TestCase):
+    """initialize() entry adapter: positional per-level arrays or the legacy names, never both."""
+
+    def test_forms(self):
+        a, b, c = object(), object(), object()
+        f = host._initial_levels
+        self.assertEqual(f((a, b), None, None, None), (a, b))
+        self.assertEqual(f((a, b, c), None, None, None), (a, b, c))
+        self.assertEqual(f((a, b), None, None, c), (a, b, c))      # old initialize(c, f, child_g=g)
+        self.assertEqual(f((), a, b, None), (a, b))
+        self.assertEqual(f((), a, b, c), (a, b, c))
+        for args in (((a,), None, None, c), ((a, b, c), None, None, c), ((a,), a, b, None),
+                     ((), a, None, None), ((), None, b, c)):
+            with self.subTest(args=args), self.assertRaises(TypeError):
+                f(*args)
+
+
 if __name__=='__main__':
     unittest.main()

tests/unit/test_runtime_child_remap.py

@@ -1,4 +1,5 @@
 """CPU mocks execute production remap lifecycle; no GPU numerical claim."""
+from dataclasses import replace
 import unittest
 from unittest.mock import patch
 
@@ -68,8 +69,11 @@ class RuntimeChildRemap(unittest.TestCase):
         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)):
+        hosts, pairs = run._boundary_hosts, run.hierarchy.pairs
+        for attribute, value in (('_wall_aware', True),
+                                 ('_boundary_hosts', [*hosts[:2], object()]),        # child boundary
+                                 ('_boundary_hosts', [hosts[0], object(), hosts[2]]),  # fine boundary
+                                 ('hierarchy', replace(run.hierarchy, pairs=(pairs[0], None)))):
             with patch.object(run, attribute, value):
                 with self.assertRaises(ValueError): run.remap_child(B)
         self.assertEqual(len(allocations), start)

tests/validation/multigrid/test_raw_exchange_gpu.py

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

tests/validation/multigrid/test_refinement_gpu.py

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