Add bounded pinned-circle manual L2 remapping

src/CelerisLab/lbm/refinement_stepper.py

@@ -3,7 +3,7 @@
 
 The caller owns the CudaContext and must close this backend before that context.
 Only synchronized postcollision_g leaf populations are public. No checkpoint,
-forcing or runtime parameter mutation is offered. Manual body-free periodic
+forcing or runtime parameter mutation is offered. Bounded manual periodic
 L2 set replacement is available at completed base intervals via remap_child.
 """
 from dataclasses import dataclass, replace
@@ -494,29 +494,54 @@ class FixedPatchStepper:
                 self._status = "failed"
                 raise
 
-    def remap_child(self, rectangles):
+    def remap_child(self, rectangles, *, child_objects=()):
         """Return a NEW ready stepper with an L2 set replaced at this base time.
 
-        Only body-free periodic three-level set children are supported; L0/L1
-        and block sizes stay fixed. The lock and stream synchronization certify
+        Periodic three-level set children only; L0/L1 and block sizes stay fixed.
+        A body run requires its original root-edge-frame ``child_objects``:
+        one stationary circle pinned in one of exactly two L2 rectangles.
+        The complete body rectangle stays identical; only the other fluid
+        rectangle may move. Current L2 Jets (including updated amplitudes)
+        are preserved. Local projected boundary objects are not valid inputs.
+        The lock and stream synchronization certify
         current postcollision G and the completed 1:2:4 counts internally.
         Both objects remain independently owned: the caller must close each,
         and chooses when to stop using the unchanged old instance. Candidate
         construction/upload errors never invalidate the old instance; failure
         to synchronize/read its device state does, as with step/download.
         """
-        from .synchronized_remap import RemapTopology, remap_synchronized
+        from .synchronized_remap import (RemapTopology, remap_synchronized,
+                                         same_physical_boundary, boundary_velocity)
 
         with self._lock:
             self._require("ready")
             if (self._wall_aware or len(self._levels) != 3
                     or not isinstance(self.child_patch, FixedPatchSet2D)
-                    or self.fine_boundary is not None or self.child_boundary is not None
-                    or any(level.boundary is not None for level in self._levels)):
-                raise ValueError("remap_child requires body-free periodic three-level L2 sets")
+                    or self.fine_boundary is not None
+                    or any(level.boundary is not None for level in self._levels[:2])):
+                raise ValueError("remap_child requires periodic three-level L2 sets with no L0/L1 bodies")
             cfg = self._remap_cfg
-            old = RemapTopology.from_config(cfg, self.patch, self.child_patch.bounds)
-            new = RemapTopology.from_config(cfg, self.patch, rectangles)
+            child_objects = tuple(child_objects)
+            actual = self._boundaries[2]
+            jets = () if actual is None else actual.jets
+            old = RemapTopology.from_config(cfg, self.patch, self.child_patch.bounds,
+                                           child_objects=child_objects, child_jets=jets)
+            if (not same_physical_boundary(old.boundary, self.child_boundary)
+                    or (self.child_boundary is not None and actual is None)):
+                raise ValueError("child_objects must reproduce the original root-edge-frame boundary")
+            if actual is not None:
+                velocity = actual.tables.get("surface_velocity")
+                if velocity is None:
+                    velocity = boundary_velocity(cfg, self.child_patch, actual.host, jets)
+                if (not same_physical_boundary(old.boundary, actual.host)
+                        or not np.array_equal(old.surface_velocity, velocity)):
+                    raise ValueError("child_objects changed boundary ownership or surface velocity")
+            new = RemapTopology.from_config(cfg, self.patch, rectangles,
+                                           child_objects=child_objects, child_jets=jets)
+            if (old.pinned_rectangle != new.pinned_rectangle
+                    or not same_physical_boundary(old.boundary, new.boundary)
+                    or not np.array_equal(old.surface_velocity, new.surface_velocity)):
+                raise ValueError("remap requires the identical pinned body rectangle and boundary")
             base_block, fine_block, child_block = self._remap_blocks
             # Check sparse block alignment/launch size before any device read/allocation.
             layout = new.child.fine_layout(child_block)
@@ -536,7 +561,8 @@ class FixedPatchStepper:
             candidate = type(self)(
                 cfg, self.patch.bounds, context=self.ctx, base_block=base_block,
                 fine_block=fine_block, child_block=child_block,
-                child_rectangles_in_parent=new.child.bounds)
+                child_rectangles_in_parent=new.child.bounds,
+                child_objects=child_objects, child_jets=jets)
             try:
                 candidate.initialize(*transferred.g)
                 candidate._step_count = count

src/CelerisLab/lbm/synchronized_remap.py

@@ -1,16 +1,21 @@
 """CPU-only whole-volume G transfer, not a runtime adaptive solver.
 
-Only body-free periodic three-level grids with a fixed L0/L1 topology are
-supported. The caller certifies that current owned leaves are quiescent,
+Periodic three-level grids with fixed L0/L1 topology are supported, including
+one stationary circle pinned inside an unchanged L2 rectangle. The caller
+certifies that current owned leaves are quiescent,
 unshifted D2Q9 postcollision G at the supplied completed base interval.
-No device buffers, boundary state, ghosts, collision or time advance occur.
+No device buffers, device boundary state, ghosts, collision or time advance
+occur; boundary preparation is host-only.
 """
-from dataclasses import dataclass
+from dataclasses import dataclass, replace
 from numbers import Integral
 
 import numpy as np
 
-from .refinement import ACTIVE, FixedPatch2D, FixedPatchSet2D, nest_patch_set
+from .refinement import (ACTIVE, FixedPatch2D, FixedPatchSet2D, FineBoundary2D,
+                         nest_patch_set, prepare_fine_boundary)
+from .descriptors import SOLID
+from ..body.jets import surface_link_velocity
 
 
 @dataclass(frozen=True, eq=False)
@@ -24,19 +29,71 @@ class RemapTopology:
     parent: object
     child: FixedPatchSet2D
     leaf_masks: tuple
+    boundary: FineBoundary2D | None = None
+    pinned_rectangle: tuple | None = None
+    surface_velocity: np.ndarray | None = None
 
     @classmethod
-    def from_config(cls, cfg, parent, child_rectangles, *, coarse_boundary="periodic", objects=()):
+    def from_config(cls, cfg, parent, child_rectangles, *, coarse_boundary="periodic", objects=(),
+                    child_objects=(), child_jets=()):
         if coarse_boundary != "periodic" or tuple(objects) or cfg.nz != 1:
-            raise ValueError("remap supports only body-free periodic 2D state")
+            raise ValueError("remap requires periodic 2D state; use child_objects for a pinned L2 body")
         if type(parent) not in (FixedPatch2D, FixedPatchSet2D):
             raise ValueError("parent must be a fixed patch or fixed patch set")
         parent = type(parent).from_config(cfg, parent.bounds)
         child, middle_kind = nest_patch_set(cfg, parent, child_rectangles)
-        masks = (~parent.covered, middle_kind == ACTIVE, child.fine_kind == ACTIVE)
+        child_objects, child_jets = tuple(child_objects), tuple(child_jets)
+        if child_objects and (len(child_objects) != 1 or len(child.bounds) != 2):
+            raise ValueError("body remap requires one stationary circle and two L2 rectangles")
+        fine_cfg = replace(cfg, nx=parent.fine_shape[1], ny=parent.fine_shape[0],
+                           viscosity=2 * cfg.viscosity)
+        boundary = prepare_fine_boundary(
+            fine_cfg, child, child_objects, jets=child_jets,
+            parent_origin=parent.fine_origin, parent_spacing=.5)
+        pinned, velocity = None, None
+        child_leaf = child.fine_kind == ACTIVE
+        if boundary is not None:
+            child_leaf &= (boundary.flags.reshape(child.fine_shape) & SOLID) == 0
+            # Whole circle containment, not merely the resolved solid cells.
+            obj = child_objects[0]
+            cx, cy = (np.asarray(obj.geometry.center) - parent.fine_origin) * 2
+            radius = obj.geometry._radius * 2
+            pinned = next(r for r in child.bounds if
+                          r[0] <= cx-radius and cx+radius <= r[2] and
+                          r[1] <= cy-radius and cy+radius <= r[3])
+            velocity = boundary_velocity(cfg, child, boundary, child_jets)
+            if velocity is not None:
+                velocity.setflags(write=False)
+        masks = (~parent.covered, middle_kind == ACTIVE, child_leaf)
         for mask in masks:
             mask.setflags(write=False)
-        return cls(parent, child, masks)
+        return cls(parent, child, masks, boundary, pinned, velocity)
+
+
+def boundary_velocity(cfg, child, boundary, jets=()):
+    """Physical cut-link order; never blocked addresses or stale Jet snapshots."""
+    local_cfg = replace(cfg, nx=child.fine_shape[1], ny=child.fine_shape[0],
+                        viscosity=4 * cfg.viscosity)
+    local_cfg.validate()
+    fluid, dirs, _, ids, rx, ry, _, _ = boundary.curved_host
+    velocity, _ = surface_link_velocity(
+        boundary.objects, local_cfg, boundary.flags, fluid, dirs, ids, rx, ry,
+        jets=jets, hits=boundary.curved_hits)
+    return velocity
+
+
+def same_physical_boundary(a, b):
+    """Compare geometry, ownership, hits and the entire BFL support stencil."""
+    if a is None or b is None:
+        return a is b
+    if not isinstance(a, FineBoundary2D) or not isinstance(b, FineBoundary2D):
+        return False
+    if (tuple(o.surface_actuation for o in a.objects)
+            != tuple(o.surface_actuation for o in b.objects)):
+        return False
+    return all(np.array_equal(x, y) for x, y in zip(
+        (a.flags, a.boundary_cells, a.curved_hits, *a.curved_host),
+        (b.flags, b.boundary_cells, b.curved_hits, *b.curved_host)))
 
 
 @dataclass(frozen=True, eq=False)
@@ -105,6 +162,10 @@ def remap_synchronized(old, new, old_g, *, base_count, level_steps, phase,
             or not _same_parent(old.parent, new.parent)
             or old.child.tau_fine != new.child.tau_fine):
         raise ValueError("remap requires unchanged L0/L1 geometry and relaxation")
+    if (old.pinned_rectangle != new.pinned_rectangle
+            or not same_physical_boundary(old.boundary, new.boundary)
+            or not np.array_equal(old.surface_velocity, new.surface_velocity)):
+        raise ValueError("remap requires the identical pinned body rectangle and physical boundary")
     if not isinstance(level_steps, (tuple, list)):
         raise ValueError("level_steps must be three completed integer counts")
     counts = (base_count, *level_steps)