Support stationary bodies and geometry Jets on nested L2

README.md

@@ -158,30 +158,48 @@ Scope and limits:
       coarse_g, fine_g = run.download()  # coarse: NaN on every covered cell;
                                          # fine: NaN on every non-ACTIVE cell
   ```
-- Three levels (one nested chain, body-free): add `child=(x0, y0, x1, y1)` in
+- Three levels (one nested chain): add `child=(x0, y0, x1, y1)` in
   **L1 window cells** (the patch's fine grid including its 4-cell ghost band) and
   `child_block`, with `fine_block` and `base_block`. L2 must sit inside L1's
   ACTIVE cells, away from L1's ghost band and transfer footprint (checked before
   any GPU allocation). `initialize` and `download` take/return three arrays; each
   level is NaN on its non-ACTIVE (and child-covered) cells. One coarse step = 2 L1
-  steps = 4 L2 steps. One L1 patch only, no bodies/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
+  dx=dt=0.25, so prescribed lattice velocity is unchanged. Each whole body and
+  its BFL stencil must be inside L2 ACTIVE cells, outside the transfer footprint.
+  A Jet may select one physical body across multiple blocks, never two bodies.
 
   ```python
+  from CelerisLab.body.geometry.circle import CircleGeometry
+  from CelerisLab.body.objects import SimObject
+  from CelerisLab.body.jets import Jet, Tube, Parallel
+
   cfg = LBMConfig(nx=16, ny=16, viscosity=.1)
+  body = SimObject(0, geometry=CircleGeometry(8., 8., .8), center=(8., 8.), radius=.8)
+  jet = Jet(Tube((6.5, 8.), 0., 3., .3), Parallel(0.), amplitude=.002)
   rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None]
   with CudaContext(0) as ctx, FixedPatchStepper(
           cfg, (4, 4, 12, 12), context=ctx, fine_block=(4, 4), base_block=(4, 4),
-          child=(8, 8, 16, 16), child_block=(4, 4)) as run:
+          child=(8, 8, 16, 16), child_block=(4, 4),
+          child_objects=[body], child_jets=[jet]) as run:
       shapes = (run.patch.coarse_shape, run.patch.fine_shape, run.child_patch.fine_shape)
       run.initialize(*(np.broadcast_to(rest, (9, *s)) for s in shapes))  # 16², 24², 24²
-      run.step(32)
-      coarse_g, l1_g, l2_g = run.download()
+      run.step(4)
+      run.set_jet_amplitude(0, -.002, level=2)  # default level=1 still addresses L1
+      run.step(4)
+      coarse_g, l1_g, l2_g = run.download()  # solids/covered/nonleaves are NaN
   ```
 - Not available: per-link surface exchange/force readout in blocked mode,
-  checkpoints, moving bodies, bodies/Jets with three levels,
+  checkpoints, moving bodies, L1 bodies/Jets in a three-level chain,
   several patches per level with nesting, dynamic refinement, `Simulation`
   integration. Validation so far is storage/operator equivalence and short runs,
-  not accuracy, convergence or long-time stability.
+  not accuracy, convergence or long-time stability. The new L2-body GPU gates
+  (layouts, masks, amplitude reversal, legacy equivalence) and independent
+  every-substep numerical replay remain required acceptance checks; CPU tests
+  alone do not establish GPU numerical correctness.
 
 ## Installation
 

src/CelerisLab/lbm/refinement.py

@@ -349,12 +349,16 @@ def nest_patch(cfg, parent, child_bounds):
     return child, _readonly(composed, np.uint8)
 
 
-def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=()):
+def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=(),
+                          parent_origin=(0., 0.), parent_spacing=1.):
     """Project fixed circles and reject boundary/transfer overlap before CUDA.
 
-    Input centers/radii are in coarse physical coordinates (dx=dt=1).
+    Input centers/radii and Jets are in the root edge frame. ``parent_origin``
+    and ``parent_spacing`` map this patch's parent grid into that frame; the
+    default preserves the L0->L1 convention. cfg must describe the parent grid.
     SurfaceActuation retains its lattice-velocity meaning: input coarse units
-    are converted into independent fine-local objects with dx=dt=.5.
+    are converted into independent fine-local objects with
+    dx=dt=parent_spacing/2.
     Original objects are not retained or mutated. No runtime geometry updates.
     ``patch`` may be a FixedPatchSet2D: objects are then projected once onto
     the level-wide fine grid (fine_origin (0, 0)); each circle must lie inside
@@ -364,6 +368,10 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=()):
     selection is checked against true fine hits here. Jet amplitudes need no
     rescaling (lattice velocity scales by dt/dx = 1 under acoustic 2:1).
     """
+    origin = np.asarray(parent_origin, dtype=np.float64)
+    if (origin.shape != (2,) or not np.isfinite(origin).all()
+            or not np.isfinite(parent_spacing) or parent_spacing <= 0):
+        raise ValueError("Parent frame must have a finite origin and positive spacing")
     if not objects:
         if jets:
             raise ValueError("Fine jets need fine boundary objects")
@@ -386,7 +394,8 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=()):
     local, links = [], []
     rects = patch.bounds if isinstance(patch, FixedPatchSet2D) else (patch.bounds,)
     covered = patch.covered
-    dx = dt = .5
+    dx = dt = .5 * parent_spacing
+    fine_origin = origin + parent_spacing * np.asarray(patch.fine_origin)
     for bid, obj in enumerate(objects):
         if (type(obj.geometry) is not CircleGeometry or obj.is_sensor or obj.is_force_region
                 or any(getattr(obj.state, name) != 0 for name in ("vx", "vy", "omega"))
@@ -394,11 +403,13 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=()):
             raise ValueError("Fine boundaries require stationary circles without other roles")
         cx, cy = obj.geometry.center
         radius = obj.geometry._radius
+        px, py = (np.asarray((cx, cy)) - origin) / parent_spacing
+        pr = radius / parent_spacing
         if (not np.isfinite([cx, cy, radius]).all() or radius <= 0
-                or not any(x0 <= cx-radius and cx+radius <= x1 and y0 <= cy-radius
-                           and cy+radius <= y1 for x0, y0, x1, y1 in rects)):
+                or not any(x0 <= px-pr and px+pr <= x1 and y0 <= py-pr
+                           and py+pr <= y1 for x0, y0, x1, y1 in rects)):
             raise ValueError("Entire fine circle must lie inside one physical patch")
-        center = tuple((v-o)/dx-.5 for v, o in zip((cx, cy), patch.fine_origin))
+        center = tuple((v-o)/dx-.5 for v, o in zip((cx, cy), fine_origin))
         spec = obj.surface_actuation
         if spec is not None:
             spec = replace(spec, normal_velocity=spec.normal_velocity * dt / dx)
@@ -465,7 +476,7 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=()):
             raise ValueError("Transfer footprint touches the coarse wall reconstruction stencil")
     from ..body.jets import fine_hits_to_edge_frame, surface_link_velocity
     hits = fine_hits_to_edge_frame([(link.hit_x, link.hit_y) for link in links],
-                                   patch.fine_origin, dx)
+                                   fine_origin, dx)
     packed = pack_cut_links_to_soa(links)
     if jets:
         # Same assignment checks as the stepper: empty/overlapping jets and

src/CelerisLab/lbm/refinement_stepper.py

@@ -93,6 +93,36 @@ class _Level:
     halo_copy: object = None
     halo_dst: object = None
     halo_src: object = None
+    boundary: object = None
+
+
+@dataclass
+class _BlockedBoundary:
+    """One level's stationary BFL data and mutable Jet velocity column."""
+    host: object
+    tables: dict
+    jets: tuple
+    actuated: bool
+    jet_cache: object = None
+    flags: object = None
+    device: object = None
+    kernel: object = None
+
+
+def _prepare_blocked_boundary(boundary, layout, cfg, jets):
+    if boundary is None or layout is None:
+        return None
+    fluid_idx, dirs, q, body_id, rx, ry, _, fallback = boundary.curved_host
+    tables, _ = blocked_link_tables(layout, boundary.flags, fluid_idx, dirs, q, fallback)
+    actuated = bool(jets) or any(o.surface_actuation is not None for o in boundary.objects)
+    if actuated:
+        cfg = replace(cfg)
+        cfg.validate()
+        tables["surface_velocity"], _ = surface_link_velocity(
+            boundary.objects, cfg, boundary.flags, fluid_idx, dirs, body_id, rx, ry,
+            jets=jets, hits=boundary.curved_hits)
+    return _BlockedBoundary(boundary, tables, jets, actuated,
+                            jet_link_shapes(jets, boundary.curved_hits) if jets else None)
 
 
 class FixedPatchStepper:
@@ -140,9 +170,18 @@ class FixedPatchStepper:
     uses one level-wide fine grid (2*ny, 2*nx): public fine arrays have that
     shape, NaN outside the patches' ACTIVE cells. The step sequence is the
     blocked one above, with every patch's transfers in one launch each.
+
+    Three-level single chain: ``child`` bounds use L1 window cells; all three
+    levels must be blocked. ``child_objects`` and ``child_jets`` use the L0
+    edge frame, with L2 dx=dt=.25 and unchanged lattice velocity. L1 bodies
+    and Jets remain unsupported in a chain. Every whole circle and its BFL
+    stencil must fit strictly away from the L2 transfer footprint. L2 runs
+    halo -> BFL -> step four times per base step. Runtime child Jet changes
+    require ``set_jet_amplitude(index, amplitude, level=2)``.
     """
     def __init__(self, cfg, bounds, *, context, coarse_boundary="periodic", fine_objects=(),
-                 fine_block=None, base_block=None, jets=(), child=None, child_block=None):
+                 fine_block=None, base_block=None, jets=(), child=None, child_block=None,
+                 child_objects=(), child_jets=()):
         if HALO != 4:
             raise RuntimeError("refinement.cu requires a four-cell fine halo")
         if coarse_boundary not in ("periodic", "slot_reservoir"):
@@ -164,7 +203,10 @@ class FixedPatchStepper:
             rects = self.patch.bounds if multipatch else (self.patch.bounds,)
             if any(r[1] < 4 or r[3] > cfg.ny - 4 for r in rects):
                 raise ValueError("slot_reservoir requires y0 >= 4 and y1 <= ny-4")
-        fine_objects = tuple(fine_objects)
+        fine_objects, jets = tuple(fine_objects), tuple(jets)
+        child_objects, child_jets = tuple(child_objects), tuple(child_jets)
+        if child is None and (child_objects or child_jets):
+            raise ValueError("child_objects/child_jets require a nested child")
         if base_block is not None and fine_block is None:
             raise ValueError("base_block requires fine_block")
         self.base_layout = None
@@ -193,7 +235,7 @@ class FixedPatchStepper:
                 raise ValueError("a nested child needs one patch plus fine_block, base_block "
                                  "and child_block")
             if tuple(fine_objects) or tuple(jets):
-                raise ValueError("a nested child does not support fine_objects or jets yet")
+                raise ValueError("a nested child supports bodies/Jets only on L2; use child_objects/child_jets")
             if len(child_block) != 2:
                 raise ValueError("child_block must be (by, bx)")
             self.child_patch, self._middle_kind = nest_patch(cfg, self.patch, child)
@@ -205,15 +247,33 @@ class FixedPatchStepper:
         self._jets = tuple(jets)
         self.fine_boundary = prepare_fine_boundary(
             cfg, self.patch, fine_objects, wall_aware=self._wall_aware, jets=self._jets)
-        self._fine_tables, self._actuated, self._jet_cache = None, False, None
-        if self.fine_layout is not None and self.fine_boundary is not None:
-            fluid_idx, dirs, q, _, _, _, _, fallback = self.fine_boundary.curved_host
-            self._fine_tables, _ = blocked_link_tables(
-                self.fine_layout, self.fine_boundary.flags, fluid_idx, dirs, q, fallback)
+        fine_cfg = replace(cfg, nx=self.patch.fine_shape[1], ny=self.patch.fine_shape[0],
+                           viscosity=2 * cfg.viscosity)
+        if self._wall_aware:
+            fine_cfg = replace(fine_cfg, y_wall_bc="bounce_back", top_wall_bc="velocity",
+                               wall_slot_width=0, wall_normal_velocity=0.0, velocity=0.0)
+        self.child_boundary = (None if self.child_patch is None else prepare_fine_boundary(
+            fine_cfg, self.child_patch, child_objects, jets=child_jets,
+            parent_origin=self.patch.fine_origin, parent_spacing=.5))
+        child_cfg = (None if self.child_patch is None else replace(
+            fine_cfg, nx=self.child_patch.fine_shape[1], ny=self.child_patch.fine_shape[0],
+            viscosity=4 * cfg.viscosity))
+        self._boundaries = {
+            1: _prepare_blocked_boundary(self.fine_boundary, self.fine_layout, fine_cfg, self._jets),
+            2: _prepare_blocked_boundary(self.child_boundary, self.child_layout, child_cfg, child_jets),
+        }
+        # Compatibility aliases for the existing L1 inspection API.
+        b = self._boundaries[1]
+        self._fine_tables = None if b is None else b.tables
+        self._actuated = False if b is None else b.actuated
+        self._jet_cache = None if b is None else b.jet_cache
         self._fine_leaf = (self.patch.fine_kind if self._middle_kind is None
                            else self._middle_kind) == ACTIVE
         self._child_leaf = (None if self.child_patch is None
                             else self.child_patch.fine_kind == ACTIVE)
+        if self.child_boundary is not None:
+            self._child_leaf &= (self.child_boundary.flags.reshape(
+                self.child_patch.fine_shape) & SOLID) == 0
         self._child = None
         if self.fine_boundary is not None:
             self._fine_leaf &= (self.fine_boundary.flags.reshape(self.patch.fine_shape) & SOLID) == 0
@@ -228,31 +288,13 @@ class FixedPatchStepper:
         self._coalesce = None
         self._link_tables = {}
         self._block = (cfg.threads_per_block, 1, 1)
-        fine_cfg = replace(cfg, nx=self.patch.fine_shape[1], ny=self.patch.fine_shape[0],
-                           viscosity=2 * cfg.viscosity)
-        if self._wall_aware:
-            fine_cfg = replace(fine_cfg, y_wall_bc="bounce_back", top_wall_bc="velocity",
-                               wall_slot_width=0, wall_normal_velocity=0.0, velocity=0.0)
-        if self._fine_tables is not None:
-            # Any actuation spec (even zero amplitude) requires the velocity buffer.
-            b = self.fine_boundary
-            self._actuated = bool(self._jets) or any(
-                obj.surface_actuation is not None for obj in b.objects)
-            if self._actuated:
-                fine_cfg.validate()  # derives nq etc.; the unique-writer key needs them
-                fluid_idx, dirs, _, body_id, rx, ry, _, _ = b.curved_host
-                self._fine_tables["surface_velocity"], _ = surface_link_velocity(
-                    b.objects, fine_cfg, b.flags, fluid_idx, dirs, body_id, rx, ry,
-                    jets=self._jets, hits=b.curved_hits)
-                if self._jets:
-                    self._jet_cache = jet_link_shapes(self._jets, b.curved_hits)
         with self.ctx.activate():
             try:
                 self._stream = cuda.Stream()
                 specs = [(replace(cfg), self.base_layout, self.patch.tau_coarse, "PatchCoarseStep"),
                          (fine_cfg, self.fine_layout, self.patch.tau_fine, "PatchFineStep")]
                 if self.child_patch is not None:
-                    specs.append((replace(fine_cfg, viscosity=4 * cfg.viscosity), self.child_layout,
+                    specs.append((child_cfg, self.child_layout,
                                   self.child_patch.tau_fine, "PatchFineStep"))
                 for index, (config, layout, tau, name) in enumerate(specs):
                     source = "refinement.cu"
@@ -284,6 +326,7 @@ class FixedPatchStepper:
                             1 if layout is None else layout.nblocks)
                     level = _Level(shape, module, self._allocate(size), self._allocate(size),
                                    module.get_function(name), grid, layout)
+                    level.boundary = self._boundaries.get(index)
                     self._levels.append(level)
                     if layout is not None:
                         level.halo_copy = module.get_function("BlockHaloCopy")
@@ -330,11 +373,15 @@ class FixedPatchStepper:
                                         self._upload(ghosts), self._upload(donors.astype(np.uint64)),
                                         self._upload(targets), self._upload(children),
                                         len(ghosts), len(targets), int(np.prod(f.layout.local_shape)))
-                if self._fine_tables is not None:
-                    self._fine_flags = self._upload(self._fine_tables["flags"])
-                    self._link_tables = {k: self._upload(v) for k, v in self._fine_tables.items()
-                                         if k != "flags"}
-                    self._curved = self._levels[1].module.get_function("BlockCurvedBoundary")
+                for level in self._levels:
+                    b = level.boundary
+                    if b is not None:
+                        b.flags = self._upload(b.tables["flags"])
+                        b.device = {k: self._upload(v) for k, v in b.tables.items() if k != "flags"}
+                        b.kernel = level.module.get_function("BlockCurvedBoundary")
+                if self._boundaries[1] is not None:
+                    b = self._boundaries[1]
+                    self._fine_flags, self._link_tables, self._curved = b.flags, b.device, b.kernel
                 elif self.fine_boundary is not None:
                     from .curved_links import CurvedLinkSoA
                     boundary = self.fine_boundary
@@ -431,58 +478,72 @@ class FixedPatchStepper:
                      cl.link_exchange_gpu if cl.link_exchange_gpu is not None else np.uintp(0),
                      grid=((cl.count + self._block[0] - 1) // self._block[0], 1), **kw)
 
-    def set_jet_amplitude(self, index, amplitude):
+    def set_jet_amplitude(self, index, amplitude, *, level=1):
         """Replace one jet's amplitude between steps; no geometry/table rebuild.
 
+        ``level=1`` preserves the L1 default; use ``level=2`` for child_jets.
         Invalid input raises before any change. A device/sync error marks the
         instance failed (as step/initialize do); host state commits only after
         a successful upload.
         """
         with self._lock, self.ctx.activate():
             self._require("new", "ready")
-            jets = list(self._jets)
+            if (isinstance(level, bool) or not isinstance(level, Integral)
+                    or level not in (1, 2) or level >= len(self._levels)):
+                raise ValueError("Jet level must be 1 or an existing child level 2")
+            boundary = self._boundaries[level]
+            jets = list(self._jets if level == 1 else (() if boundary is None else boundary.jets))
             jets[index] = jets[index].with_amplitude(amplitude)  # validates first
             velocity = None
-            if self._fine_curved is None and self._jet_cache is not None:
-                shape, jet_id = self._jet_cache
-                velocity = self._fine_tables["surface_velocity"].copy()
+            if boundary is not None and boundary.jet_cache is not None:
+                shape, jet_id = boundary.jet_cache
+                velocity = boundary.tables["surface_velocity"].copy()
                 driven = jet_id >= 0
                 velocity[driven] = jet_link_velocity(jets, shape, jet_id)[driven]
             try:
-                if self._fine_curved is not None:
+                if level == 1 and self._fine_curved is not None:
                     self._fine_curved.set_jet_amplitude(index, amplitude, self._stream)
                 elif velocity is not None:
                     self._stream.synchronize()  # no in-flight step reads the column
-                    cuda.memcpy_htod(self._link_tables["surface_velocity"], velocity)
+                    cuda.memcpy_htod(boundary.device["surface_velocity"], velocity)
             except BaseException:
                 self._status = "failed"
                 raise
             if velocity is not None:
-                self._fine_tables["surface_velocity"] = velocity
-            self._jets = tuple(jets)
+                boundary.tables["surface_velocity"] = velocity
+            if boundary is not None:
+                boundary.jets = tuple(jets)
+            if level == 1:
+                self._jets = tuple(jets)
 
     def _blocked_fine_boundary(self, buffer, **kw):
-        if self._fine_tables is None:
+        self._blocked_boundary(self._levels[1], buffer, **kw)
+
+    def _blocked_boundary(self, level, buffer, **kw):
+        b = level.boundary
+        if b is None:
             return
-        if self._curved is None or self._fine_flags is None:
-            raise RuntimeError("Fine boundary device buffers are missing")
-        t, count = self._link_tables, len(self._fine_tables["cell"])
+        if b.kernel is None or b.flags is None:
+            raise RuntimeError("Boundary device buffers are missing")
+        t, count = b.device, len(b.tables["cell"])
         velocity = t.get("surface_velocity")
-        if self._actuated and velocity is None:
+        if b.actuated and velocity is None:
             raise RuntimeError("Surface velocity device buffer is missing")
-        self._curved(buffer, t["block"], t["cell"], t["dir"], t["q"], t["fallback"],
-                     np.uint32(count), velocity if velocity is not None else np.uintp(0),
-                     grid=((count + self._block[0] - 1) // self._block[0], 1), **kw)
+        b.kernel(buffer, t["block"], t["cell"], t["dir"], t["q"], t["fallback"],
+                 np.uint32(count), velocity if velocity is not None else np.uintp(0),
+                 grid=((count + self._block[0] - 1) // self._block[0], 1), **kw)
 
     def _advance_child(self, parent_buffer, **kw):
-        """Explode from the given L1 buffer, then two L2 hops (halo -> step)."""
+        """Explode from the given L1 buffer, then two L2 hops (halo -> BFL -> step)."""
         pair, g = self._child, self._levels[2]
         pair.explode(parent_buffer, g.current, pair.ghosts, pair.donors, np.uint32(pair.n_ghosts),
                      np.uint64(pair.parent_cells),
                      grid=((pair.n_ghosts + self._block[0] - 1) // self._block[0], 1), **kw)
         for src, dst in ((g.current, g.spare), (g.spare, g.current)):
             self._gather_halo(g, src, **kw)
-            g.kernel(src, dst, np.uint32(g.layout.halo), np.uintp(0), g.kind, grid=g.grid, **kw)
+            self._blocked_boundary(g, src, **kw)
+            flags = np.uintp(0) if g.boundary is None else g.boundary.flags
+            g.kernel(src, dst, np.uint32(g.layout.halo), flags, g.kind, grid=g.grid, **kw)
 
     def _gather_halo(self, level, buffer, **kw):
         pairs = len(level.layout.halo_dst)
@@ -589,6 +650,9 @@ class FixedPatchStepper:
                         self._stream.synchronize()
                 except Exception as exc:
                     error = exc
+                for b in getattr(self, "_boundaries", {}).values():
+                    if b is not None:
+                        b.flags = b.device = b.kernel = None
                 self._levels.clear()
                 self._explode = self._coalesce = self._child = None
                 if getattr(self, "_fine_curved", None) is not None:

tests/unit/test_refinement.py

@@ -310,5 +310,74 @@ class NestPatchTests(unittest.TestCase):
                 refinement.nest_patch(self.cfg, self.l1, bounds)
 
 
+class ChildBoundaryTests(unittest.TestCase):
+    @staticmethod
+    def body(center=(8., 8.), radius=.8, obj_id=99, spec=None):
+        from CelerisLab.body.geometry.circle import CircleGeometry
+        from CelerisLab.body.objects import SimObject
+        return SimObject(obj_id, geometry=CircleGeometry(*center, radius), center=center,
+                         radius=radius, surface_actuation=spec)
+
+    def prepare(self, objects, jets=(), large=False):
+        from dataclasses import replace
+        cfg = LBMConfig(nx=32 if large else 16, ny=32 if large else 16, viscosity=.1)
+        parent = FixedPatch2D.from_config(cfg, (4, 4, 28, 28) if large else (4, 4, 12, 12))
+        child, _ = refinement.nest_patch(cfg, parent, (8, 8, 40, 40) if large else (8, 8, 16, 16))
+        middle = replace(cfg, nx=parent.fine_shape[1], ny=parent.fine_shape[0], viscosity=.2)
+        return child, refinement.prepare_fine_boundary(middle, child, objects, jets=jets,
+                         parent_origin=parent.fine_origin, parent_spacing=.5)
+
+    def test_root_frame_matches_independent_parent_frame_projection(self):
+        from CelerisLab.body.jets import Jet, Parallel, Tube, jet_link_shapes
+        jet = Jet(Tube((6.5, 8.), 0., 3., .3), Parallel(0.), amplitude=.002)
+        original = self.body()
+        child, b = self.prepare([original], [jet])
+        middle = LBMConfig(nx=24, ny=24, viscosity=.2)
+        transformed = self.body((12., 12.), 1.6)
+        local_jet = Jet(Tube((9., 12.), 0., 6., .6), Parallel(0.), amplitude=.002)
+        other = refinement.prepare_fine_boundary(middle, child, [transformed], jets=[local_jet])
+        for a, c in zip(b.curved_host, other.curved_host):
+            np.testing.assert_array_equal(a, c)
+        np.testing.assert_array_equal(b.flags, other.flags)
+        np.testing.assert_allclose(b.curved_hits, 2 + .5 * other.curved_hits, atol=0, rtol=0)
+        for a, c in zip(jet_link_shapes([jet], b.curved_hits),
+                        jet_link_shapes([local_jet], other.curved_hits)):
+            np.testing.assert_allclose(a, c, atol=1e-14, rtol=0)
+        self.assertEqual(len(b.curved_host[0]), 60)
+        self.assertAlmostEqual(child.tau_fine, .5 + 12 * .1)
+        self.assertEqual(original.geometry.center, (8., 8.))
+        self.assertEqual(original.obj_id, 99)
+
+    def test_invalid_body_and_jet_contracts(self):
+        from CelerisLab.body.jets import Jet, Parallel, Tube
+        from CelerisLab.body.objects import SurfaceActuation
+        jet = Jet(Tube((6.5, 8.), 0., 3., .3), Parallel(0.), amplitude=.002)
+        for body in (self.body((6., 8.)), self.body((6.9, 8.)), self.body(radius=2.1)):
+            with self.assertRaises(ValueError):
+                self.prepare([body])
+        moving = self.body()
+        moving.state.vx = .01
+        with self.assertRaisesRegex(ValueError, "stationary"):
+            self.prepare([moving])
+        with self.assertRaises(ValueError):
+            self.prepare([self.body()], [jet, jet])
+        empty = Jet(Tube((0., 0.), 0., 1., .1), Parallel(0.), amplitude=0.)
+        with self.assertRaises(ValueError):
+            self.prepare([self.body()], [empty])
+        spec = SurfaceActuation(.7, .85, .002)
+        with self.assertRaisesRegex(ValueError, "SurfaceActuation"):
+            self.prepare([self.body(spec=spec)], [spec.as_jet((8., 8.), .8)])
+
+    def test_one_jet_two_real_bodies_same_caller_id_rejected(self):
+        from CelerisLab.body.jets import Jet, Parallel, Tube
+        bodies = [self.body((12., 14.), 1., 7), self.body((17., 14.), 1., 7)]
+        _, boundary = self.prepare(bodies, large=True)  # establish legal geometry first
+        self.assertEqual(np.unique(boundary.curved_host[3]).tolist(), [0, 1])
+        for amplitude in (0., .002):
+            jet = Jet(Tube((10., 14.), 0., 9., .4), Parallel(0.), amplitude=amplitude)
+            with self.assertRaisesRegex(ValueError, "one jet may drive only one body"):
+                self.prepare(bodies, [jet], large=True)
+
+
 if __name__ == "__main__":
     unittest.main(verbosity=2)

tests/validation/multigrid/test_refinement_three_level.py

@@ -10,8 +10,10 @@ Run CPU and GPU classes in separate processes. GPU gates (frozen before any run)
   G written by the first hop), never from L0: 2 child advances per coarse step.
 Numerical comparison against the CPU probe (G vs C(f)) is done by the acceptance run.
 """
+from contextlib import nullcontext
+from dataclasses import replace
 import unittest
-from unittest.mock import patch
+from unittest.mock import Mock, patch
 
 import numpy as np
 import pycuda.driver as cuda
@@ -28,6 +30,32 @@ L1, L2 = (4, 4, 12, 12), (8, 8, 16, 16)
 KW = dict(fine_block=(4, 4), base_block=(4, 4), child=L2, child_block=(4, 4))
 
 
+def child_body(spec=None):
+    from CelerisLab.body.geometry.circle import CircleGeometry
+    from CelerisLab.body.objects import SimObject
+    return SimObject(99, geometry=CircleGeometry(8., 8., .8), center=(8., 8.), radius=.8,
+                     surface_actuation=spec)
+
+
+def child_jet(speed=.002):
+    from CelerisLab.body.jets import Jet, Tube, Parallel
+    return Jet(Tube((6.5, 8.), 0., 3., .3), Parallel(0.), amplitude=speed)
+
+
+def mocked_child():
+    from CelerisLab.cuda import compiler_v2 as compiler
+    context = Mock(activate=nullcontext, sm_arch="sm_70")
+    module = Mock(get_function=Mock(side_effect=lambda name: Mock(name=name)))
+    module.get_global.return_value = (Mock(), 20)
+    with patch.object(cuda, "Stream"), patch.object(compiler, "generate_config"), \
+            patch.object(compiler, "compile_kernel"), patch.object(compiler, "cleanup_build_dir"), \
+            patch.object(compiler, "load_module", return_value=module), \
+            patch.object(cuda, "mem_alloc", side_effect=lambda size: Mock()), \
+            patch.object(cuda, "memcpy_htod"):
+        return FixedPatchStepper(CFG, L1, context=context, child_objects=[child_body()],
+                                 child_jets=[child_jet()], **KW)
+
+
 class ThreeLevelCPU(unittest.TestCase):
     def test_rejections_before_cuda(self):
         context = object()  # never activated: every case must fail first
@@ -37,11 +65,92 @@ class ThreeLevelCPU(unittest.TestCase):
         for kw in bad:
             with self.subTest(kw=kw), self.assertRaises(ValueError):
                 FixedPatchStepper(CFG, L1, context=context, **kw)
+        for kw in (dict(child_objects=[child_body()]), dict(child_jets=[child_jet()]),
+                   dict(KW, child_jets=[child_jet()])):
+            with self.assertRaises(ValueError):
+                FixedPatchStepper(CFG, L1, context=context, **kw)
         with self.assertRaises(ValueError):  # several L1 patches + child: not yet
             FixedPatchStepper(LBMConfig(nx=32, ny=16, viscosity=.1), ((4, 4, 12, 12), (20, 4, 28, 12)),
                               context=context, **KW)
 
 
+    def test_four_halo_bfl_steps_and_fresh_middle_feeds(self):
+        run = mocked_child()
+        run._status = "ready"
+        f, g = run._levels[1:]
+        b = g.boundary
+        events, feeds = [], []
+        advance, gather = run._advance_child, run._gather_halo
+        def record_feed(buffer, **kw):
+            self.assertIs(buffer, (f.current, f.spare)[len(feeds)])
+            feeds.append(buffer)
+            advance(buffer, **kw)
+        def record_halo(level, buffer, **kw):
+            if level is g:
+                events.append(("halo", buffer))
+            gather(level, buffer, **kw)
+        b.kernel.side_effect = lambda buffer, *a, **k: events.append(("bfl", buffer))
+        def record_step(src, dst, halo, flags, kind, **kw):
+            self.assertIs(flags, b.flags)
+            self.assertIs(kind, g.kind)
+            events.append(("step", src))
+        g.kernel.side_effect = record_step
+        with patch.object(run, "_advance_child", side_effect=record_feed), \
+                patch.object(run, "_gather_halo", side_effect=record_halo):
+            run.step()
+        self.assertEqual(feeds, [f.current, f.spare])
+        expected = [(kind, buffer) for buffer in (g.current, g.spare, g.current, g.spare)
+                    for kind in ("halo", "bfl", "step")]
+        self.assertEqual(events, expected)
+        self.assertGreater(len(np.unique(b.tables["block"])), 1)
+        # Negative control: reusing L1's old source at the midpoint must fail.
+        feeds.clear()
+        with patch.object(run, "_advance_child", side_effect=lambda buffer, **kw:
+                          record_feed(f.current, **kw)):
+            with self.assertRaises(AssertionError):
+                run.step()
+        run.close()
+
+    def test_child_update_isolated_and_failure_atomic(self):
+        run = mocked_child()
+        b = run._levels[2].boundary
+        before = b.tables["surface_velocity"].copy()
+        with self.assertRaises(IndexError):
+            run.set_jet_amplitude(0, -.002)  # default stays L1
+        for level in (0, 3, True):
+            with self.assertRaises(ValueError):
+                run.set_jet_amplitude(0, -.002, level=level)
+        with patch.object(cuda, "memcpy_htod") as upload:
+            run.set_jet_amplitude(0, -.002, level=2)
+            self.assertEqual(upload.call_count, 1)
+            self.assertIs(upload.call_args.args[0], b.device["surface_velocity"])
+        np.testing.assert_array_equal(b.tables["surface_velocity"], -before)
+        self.assertEqual(run._jets, ())
+        with patch.object(cuda, "memcpy_htod", side_effect=RuntimeError("copy failed")):
+            with self.assertRaisesRegex(RuntimeError, "copy failed"):
+                run.set_jet_amplitude(0, .002, level=2)
+        self.assertEqual(run._status, "failed")
+        self.assertEqual(b.jets[0].amplitude, -.002)
+        np.testing.assert_array_equal(b.tables["surface_velocity"], -before)
+        run.close()
+
+    def test_child_public_masks_include_solids(self):
+        from CelerisLab.lbm.descriptors import SOLID
+        run = mocked_child()
+        solid = run.child_boundary.flags.reshape(run.child_patch.fine_shape) & SOLID != 0
+        self.assertTrue(solid.any())
+        self.assertFalse(run._child_leaf[solid].any())
+        arrays = [np.ones(level.shape, np.float32) for level in run._levels]
+        run._mask_nonleaves(*arrays)
+        self.assertTrue(np.isnan(arrays[2][:, solid]).all())
+        with patch.object(cuda, "memcpy_htod") as upload:
+            run.initialize(*arrays)
+        sent = upload.call_args_list[-1].args[1]
+        dense = run.child_layout.to_dense(sent)
+        self.assertTrue(np.isnan(dense[:, ~run._child_leaf]).all())
+        run.close()
+
+
 class ThreeLevelGPU(unittest.TestCase):
     def test_equilibrium_is_kept_on_every_level(self):
         from CelerisLab.cuda.context import CudaContext
@@ -85,5 +194,67 @@ class ThreeLevelGPU(unittest.TestCase):
             self.assertTrue(np.isnan(rg[:, p2.fine_kind == OUTER_GHOST]).all())
 
 
+class ThreeLevelBodyGPU(unittest.TestCase):
+    """Frozen gates; CPU mocks do not establish these numerical results."""
+    def run_case(self, mode, block, *, jets=(), spec=None, switch=False):
+        from CelerisLab.cuda.context import CudaContext
+        from CelerisLab.lbm.descriptors import SOLID
+        cfg = CFG if mode == "periodic" else replace(
+            CFG, velocity=0., wall_normal_velocity=.001, wall_slot_width=4,
+            y_wall_bc="on_site_velocity", top_wall_bc="reservoir")
+        with CudaContext(0) as ctx, FixedPatchStepper(
+                cfg, L1, context=ctx, coarse_boundary=mode, child_objects=[child_body(spec)],
+                child_jets=jets, **dict(KW, child_block=block)) as run:
+            run.initialize(*(state(*level.shape[1:], i + 1) for i, level in enumerate(run._levels)))
+            snapshots = []
+            for stop in (0, 1, 2, 4, 8):
+                run.step(stop - run.step_count)
+                if stop == 4 and switch:
+                    run.set_jet_amplitude(0, -.002, level=2)
+                public, stored = run.download(), raw(run)
+                leaves = [~run.patch.covered, run._fine_leaf, run._child_leaf]
+                if mode != "periodic":
+                    leaves[0][[0, -1]] = False
+                for g, mask in zip(public, leaves):
+                    self.assertTrue(np.isfinite(g[:, mask]).all())
+                    self.assertTrue(np.isnan(g[:, ~mask]).all())
+                invalid = [run.patch.covered.copy(), run._middle_kind == OUTER_GHOST,
+                           (run.child_patch.fine_kind == OUTER_GHOST) |
+                           ((run.child_boundary.flags.reshape(run.child_patch.fine_shape) & SOLID) != 0)]
+                if mode != "periodic":
+                    invalid[0][[0, -1]] = True
+                for g, mask in zip(stored, invalid):
+                    self.assertTrue(np.isnan(g[:, mask]).all())
+                snapshots.append(public + tuple(stored))
+            return snapshots
+
+    def test_layouts_static_driven_switched_and_zero(self):
+        for mode in ("periodic", "slot_reservoir"):
+            cases = [self.run_case(mode, (4, 4), jets=jets, switch=switch)
+                     for jets, switch in (((), False), ([child_jet()], False),
+                                          ([child_jet()], True), ([child_jet(0.)], False))]
+            for a, b in zip(cases[0], cases[3]):
+                for x, y in zip(a, b):
+                    np.testing.assert_array_equal(x, y)
+            self.assertFalse(np.array_equal(cases[1][-1][2], cases[2][-1][2], equal_nan=True))
+            for block in ((8, 4), (4, 8)):
+                for expected, (jets, switch) in zip(cases, (((), False), ([child_jet()], False),
+                                                          ([child_jet()], True), ([child_jet(0.)], False))):
+                    actual = self.run_case(mode, block, jets=jets, switch=switch)
+                    for a, b in zip(actual, expected):
+                        for x, y in zip(a, b):
+                            np.testing.assert_array_equal(x, y)
+
+    def test_legacy_surface_equivalence(self):
+        from CelerisLab.body.objects import SurfaceActuation
+        spec = SurfaceActuation(.7, .85, .002)
+        for mode in ("periodic", "slot_reservoir"):
+            legacy = self.run_case(mode, (4, 4), spec=spec)[-1][2]
+            driven = self.run_case(mode, (4, 4), jets=[spec.as_jet((8., 8.), .8)])[-1][2]
+            mask = np.isfinite(legacy)
+            np.testing.assert_array_equal(mask, np.isfinite(driven))
+            self.assertLessEqual(np.max(np.abs(legacy[mask] - driven[mask])), 1e-6)
+
+
 if __name__ == "__main__":
     unittest.main()