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>

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.

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