Blocked fine storage: velocity-only surface actuation on fixed circles

- curved_links.surface_link_velocity: pure host helper extracted from
  CurvedLinkSoA.assign_surface (which now calls it; dense path unchanged).
  Rows follow physical cut-link order.
- blocks.cu BlockCurvedBoundary: one nullable trailing surface_velocity
  argument (2 floats/link), same convention as CurvedBoundaryKernel; no
  exchange output. Both call sites pass it explicitly (PeriodicBlockStepper
  passes null; single-level behaviour unchanged).
- FixedPatchStepper blocked mode accepts actuated stationary circles,
  validated against the physical fine_cfg; uploads the velocity table only
  when any actuation spec exists and fails before BFL if that buffer is
  missing (even at zero amplitude). Blocked exchange diagnostics are
  unavailable; dense exchange behaviour unchanged. No force aggregation,
  no motion, no new numerical formula.
Tests: 1x1 seam (0,0), receivers in 4 blocks, all q<.5 links fallback 0;
dense SoA never built in blocked mode; launch order incl. null vs velocity
pointer; missing-buffer guard at 0 and 0.002; helper == assign_surface;
aperture crosses a seam; GPU: periodic and slot_reservoir, +/-0.002,
8 coarse steps, dense vs 1x1/2x2/4x2 exact values + NaN masks.

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

src/CelerisLab/lbm/block_stepper.py

@@ -210,7 +210,8 @@ class PeriodicBlockStepper:
                     if self._curved is not None:
                         t = self._tables
                         self._curved(self._current, t["block"], t["cell"], t["dir"], t["q"],
-                                     t["fallback"], np.uint32(self._links), grid=self._link_grid, **kw)
+                                     t["fallback"], np.uint32(self._links), np.uintp(0),  # no actuation
+                                     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
                                grid=self._grid, **kw)

src/CelerisLab/lbm/curved_links.py

@@ -42,6 +42,33 @@ def _upload_columns(columns, stream=None):
         stream.synchronize()
 
 
+def surface_link_velocity(objects, cfg, flags, fluid_idx, dir, body_id, rx, ry):
+    """Per-link prescribed surface velocity (n, 2) FP32 and aperture mask.
+
+    Host-only and storage-neutral: rows follow the given physical cut-link
+    order. Returns (None, None) when no object carries surface actuation.
+    """
+    enabled = [obj for obj in objects if obj.surface_actuation is not None]
+    if not enabled:
+        return None, None
+    from .descriptors import FLUID, SOLID
+    if np.any((flags[fluid_idx] & FLUID) == 0) or np.any(flags[fluid_idx] & SOLID):
+        raise ValueError("Surface cut links must have active fluid receivers")
+    if len(np.unique(fluid_idx.astype(np.uint64)*cfg.nq+dir)) != len(dir):
+        raise ValueError("Surface boundary requires unique cut-link writers")
+    surface_velocity = np.zeros((len(dir), 2), np.float32)
+    surface_aperture = np.zeros(len(dir), bool)
+    for obj in enabled:
+        obj.validate_surface_actuation(cfg)
+        selected = body_id == obj.obj_id
+        velocity, aperture = obj.surface_actuation.project(rx[selected], ry[selected])
+        if not aperture.any():
+            raise ValueError("Surface aperture hits no resolved cut links")
+        surface_velocity[selected] = velocity
+        surface_aperture[selected] = aperture
+    return surface_velocity, surface_aperture
+
+
 @dataclass
 class CurvedLinkSoA:
     """Cut-link columns on host with matching ``*_gpu`` device buffers."""
@@ -83,25 +110,9 @@ class CurvedLinkSoA:
 
     def assign_surface(self, objects, cfg, flags):
         """Keep the original eight-column geometry/checkpoint contract unchanged."""
-        enabled = [obj for obj in objects if obj.surface_actuation is not None]
         self.surface_velocity = self.aperture = None
-        if not enabled:
-            return
-        from .descriptors import FLUID, SOLID
-        if np.any((flags[self.fluid_idx] & FLUID) == 0) or np.any(flags[self.fluid_idx] & SOLID):
-            raise ValueError("Surface cut links must have active fluid receivers")
-        if len(np.unique(self.fluid_idx.astype(np.uint64)*cfg.nq+self.dir)) != len(self.dir):
-            raise ValueError("Surface boundary requires unique cut-link writers")
-        self.surface_velocity = np.zeros((len(self.dir), 2), np.float32)
-        self.aperture = np.zeros(len(self.dir), bool)
-        for obj in enabled:
-            obj.validate_surface_actuation(cfg)
-            selected = self.body_id == obj.obj_id
-            velocity, aperture = obj.surface_actuation.project(self.rx[selected], self.ry[selected])
-            if not aperture.any():
-                raise ValueError("Surface aperture hits no resolved cut links")
-            self.surface_velocity[selected] = velocity
-            self.aperture[selected] = aperture
+        self.surface_velocity, self.aperture = surface_link_velocity(
+            objects, cfg, flags, self.fluid_idx, self.dir, self.body_id, self.rx, self.ry)
 
     def surface_exchange(self, body_id):
         """Latest boundary step's (Fx,Fy,Tz), split by aperture; not nozzle thrust.

src/CelerisLab/lbm/kernels/blocks.cu

@@ -54,14 +54,22 @@ extern "C" __global__ void BlockExplode(
 // Static cut links, receiver-local staging: cell is the receiver's local core
 // 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.
 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* q, const unsigned char* fallback, unsigned int count,
+    const float* surface_velocity)
 {
     unsigned int row = blockIdx.x * blockDim.x + threadIdx.x;
     if (row >= count) return;
+    float Uw = 0.0f, Vw = 0.0f;
+    if (surface_velocity != nullptr) {
+        const float u = surface_velocity[2 * row], v = surface_velocity[2 * row + 1];
+        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,
-                       0.0f, 0.0f, 0.0f, 0.0f, fallback[row], nullptr, -1);
+                       Uw, Vw, 0.0f, 0.0f, fallback[row], nullptr, -1);
 }
 
 // Pull/collide/store core cells only; halo cells of out are left untouched.

src/CelerisLab/lbm/refinement_stepper.py

@@ -15,6 +15,7 @@ import pycuda.driver as cuda
 
 from ..cuda import compiler_v2 as compiler
 from .blocks import BlockLayout2D, blocked_link_tables
+from .curved_links import surface_link_velocity
 from .block_stepper import MAX_GRID_DIM
 from .refinement import ACTIVE, HALO, OUTER_GHOST, FixedPatch2D, prepare_fine_boundary
 from .descriptors import SOLID
@@ -71,9 +72,11 @@ class FixedPatchStepper:
     same-level halos. Patch kinds stay per cell (OUTER invalid, INNER stream
     only, ACTIVE collide); halos are gathered before each fine hop. Public
     arrays keep the dense (9, ny, nx) fine shape; only the allocation differs.
-    Blocked fine_objects must be passive stationary circles (no actuation):
+    Blocked fine_objects are stationary circles, optionally surface-actuated:
     prepare_fine_boundary decides legality on the dense fine grid, then cut
     links are staged receiver-local per hop: halo -> BFL -> step, same buffer.
+    Blocked mode applies prescribed surface velocity only; its per-link
+    exchange diagnostic is unavailable (rows keep physical cut-link order).
     """
     def __init__(self, cfg, bounds, *, context, coarse_boundary="periodic", fine_objects=(),
                  fine_block=None):
@@ -93,8 +96,6 @@ class FixedPatchStepper:
         fine_objects = tuple(fine_objects)
         self.fine_layout = None
         if fine_block is not None:
-            if any(getattr(obj, "surface_actuation", None) is not None for obj in fine_objects):
-                raise ValueError("fine_block supports passive circles only; actuation needs the dense fine level")
             if len(fine_block) != 2:
                 raise ValueError("fine_block must be (by, bx)")
             fy, fx = self.patch.fine_shape
@@ -104,7 +105,7 @@ class FixedPatchStepper:
                 raise ValueError("fine blocks exceed launch limits")
         self.fine_boundary = prepare_fine_boundary(
             cfg, self.patch, fine_objects, wall_aware=self._wall_aware)
-        self._fine_tables = None
+        self._fine_tables, self._actuated = None, False
         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(
@@ -128,6 +129,14 @@ class FixedPatchStepper:
         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 = any(obj.surface_actuation is not None for obj in b.objects)
+            if self._actuated:
+                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)
         with self.ctx.activate():
             try:
                 self._stream = cuda.Stream()
@@ -286,8 +295,12 @@ class FixedPatchStepper:
         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"])
+        velocity = t.get("surface_velocity")
+        if self._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), grid=((count + self._block[0] - 1) // self._block[0], 1), **kw)
+                     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 _fine_halo(self, buffer, **kw):
         pairs = len(self.fine_layout.halo_dst)

tests/validation/multigrid/test_refinement_blocked_fine.py

@@ -16,6 +16,7 @@ from CelerisLab.body.geometry.circle import CircleGeometry
 from CelerisLab.body.objects import SimObject, SurfaceActuation
 from CelerisLab.config import LBMConfig
 from CelerisLab.cuda import compiler_v2 as compiler
+from CelerisLab.lbm import curved_links
 from CelerisLab.lbm.blocks import NO_DONOR, BlockLayout2D, blocked_link_tables
 from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
 from CelerisLab.lbm.descriptors import SOLID
@@ -47,16 +48,31 @@ def state(ny, nx, seed):
     return (feq + noise).astype(np.float32)
 
 
-def circle(actuation=None):
+def circle(spec=None):
     """Astra fixture: coarse cell-edge C=(16,12), R=2.15 -> fine C=(19.5,15.5), R=4.3."""
     return SimObject(0, geometry=CircleGeometry(16.0, 12.0, 2.15), center=(16.0, 12.0),
-                     radius=2.15, surface_actuation=actuation)
+                     radius=2.15, surface_actuation=spec)
 
 
-def fine_circle(mode):
+def actuation(speed):
+    """Aperture centred on +x: hits links straddling the horizontal y 15|16 seam."""
+    return SurfaceActuation(0.0, 0.6, speed)
+
+
+def fine_circle(mode, spec=None):
     cfg = config(mode)
     patch = FixedPatch2D.from_config(cfg, BOUNDS)
-    return patch, prepare_fine_boundary(cfg, patch, (circle(),), wall_aware=mode == "slot_reservoir")
+    return patch, prepare_fine_boundary(cfg, patch, (circle(spec),),
+                                        wall_aware=mode == "slot_reservoir")
+
+
+def fine_cfg(mode, patch):
+    from dataclasses import replace
+    cfg = replace(config(mode), nx=patch.fine_shape[1], ny=patch.fine_shape[0], viscosity=.2)
+    if mode == "slot_reservoir":
+        cfg = replace(cfg, y_wall_bc="bounce_back", top_wall_bc="velocity",
+                      wall_slot_width=0, wall_normal_velocity=0.0, velocity=0.0)
+    return cfg
 
 
 def seam_cells(layout, mapped):
@@ -86,9 +102,7 @@ class BlockedFineCPU(unittest.TestCase):
         return Mock(activate=Mock(side_effect=AssertionError("CUDA reached")), sm_arch="sm_70")
 
     def test_rejections_before_cuda(self):
-        actuated = circle(SurfaceActuation(0.0, 0.5, 0.002))
-        for kw in ({"fine_block": (7, 7)}, {"fine_block": (16,)},
-                   {"fine_block": (16, 20), "fine_objects": [actuated]}):
+        for kw in ({"fine_block": (7, 7)}, {"fine_block": (16,)}):
             with self.subTest(kw=kw), self.assertRaises(ValueError):
                 FixedPatchStepper(config("periodic"), BOUNDS, context=self.forbidden_context(), **kw)
 
@@ -123,10 +137,13 @@ class BlockedFineCPU(unittest.TestCase):
             fluid_idx, dirs, q, _, _, _, _, fallback = boundary.curved_host
             self.assertEqual(len(fluid_idx), 76, mode)
             self.assertEqual(int(np.count_nonzero(q < .5)), 68, mode)
-            for by, bx in FINE_BLOCKS[1:]:
+            self.assertTrue((fallback[q < .5] == 0).all(), mode)  # every q<.5 link uses its donor
+            for by, bx in FINE_BLOCKS:
                 layout = BlockLayout2D.build(*patch.fine_shape, by, bx, periodic=(False, False))
                 _, mapped = blocked_link_tables(layout, boundary.flags, fluid_idx, dirs, q, fallback)
-                self.assertEqual(seam_cells(layout, mapped), (8, 8), (mode, by, bx))
+                single = layout.nblocks == 1
+                self.assertEqual(seam_cells(layout, mapped), (0, 0) if single else (8, 8), (mode, by, bx))
+                self.assertEqual(len(set(mapped.block.tolist())), 1 if single else 4, (mode, by, bx))
 
     def test_storage_adapter_flags_sizing(self):
         patch, boundary = fine_circle("periodic")
@@ -142,7 +159,8 @@ class BlockedFineCPU(unittest.TestCase):
             with self.assertRaises(ValueError):
                 blocked_link_tables(layout, bad, *columns)
 
-    def test_same_buffer_halo_bfl_step_launch_order(self):
+    def mocked_backend(self, objects):
+        """Real constructor/initialize with driver doubles; kernels record (name, args)."""
         context = Mock(activate=nullcontext, sm_arch="sm_70")
         calls, functions = [], {}
 
@@ -154,26 +172,77 @@ class BlockedFineCPU(unittest.TestCase):
                 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"):
+                patch.object(cuda, "memcpy_htod"), \
+                patch.object(curved_links, "CurvedLinkSoA", side_effect=AssertionError("dense SoA built")):
             backend = FixedPatchStepper(config("periodic"), BOUNDS, context=context,
-                                        fine_objects=[circle()], fine_block=FINE_BLOCKS[1])
-            g_fine = state(*backend.patch.fine_shape, 2)
-            backend.initialize(state(NY, NX, 1), g_fine)
+                                        fine_objects=objects, fine_block=FINE_BLOCKS[1])
+            backend.initialize(state(NY, NX, 1), state(*backend.patch.fine_shape, 2))
+        return backend, calls
+
+    def test_same_buffer_halo_bfl_step_launch_order(self):
+        for spec in (None, actuation(.002)):
+            backend, calls = self.mocked_backend([circle(spec)])
             c, f = backend._levels
             coarse, cur, spare = c.current, f.current, f.spare
             flags, kind = backend._fine_flags, backend._kind
+            velocity = backend._link_tables.get("surface_velocity")
+            self.assertEqual(velocity is None, spec is None)
             backend.step(1)
             backend.close()
-        order = [(name, a[0]) for name, a in calls]
-        self.assertEqual(order, [("BlockExplode", coarse),
-                                 ("BlockHaloCopy", cur), ("BlockCurvedBoundary", cur), ("BlockStep", cur),
-                                 ("BlockHaloCopy", spare), ("BlockCurvedBoundary", spare),
-                                 ("BlockStep", spare), ("PatchCoarseStep", coarse)])
-        steps = [a for name, a in calls if name == "BlockStep"]
-        self.assertEqual([a[1] for a in steps], [spare, cur])
-        for a in steps:
-            self.assertIs(a[3], flags)  # blocked flags, not null
-            self.assertIs(a[4], kind)
+            order = [(name, a[0]) for name, a in calls]
+            self.assertEqual(order, [("BlockExplode", coarse),
+                                     ("BlockHaloCopy", cur), ("BlockCurvedBoundary", cur), ("BlockStep", cur),
+                                     ("BlockHaloCopy", spare), ("BlockCurvedBoundary", spare),
+                                     ("BlockStep", spare), ("PatchCoarseStep", coarse)])
+            steps = [a for name, a in calls if name == "BlockStep"]
+            self.assertEqual([a[1] for a in steps], [spare, cur])
+            for a in steps:
+                self.assertIs(a[3], flags)  # blocked flags, not null
+                self.assertIs(a[4], kind)
+            for a in (a for name, a in calls if name == "BlockCurvedBoundary"):
+                if spec is None:  # explicit null at the raw launch ABI
+                    self.assertIsInstance(a[7], np.uintp)
+                    self.assertEqual(int(a[7]), 0)
+                else:
+                    self.assertIs(a[7], velocity)
+
+    def test_actuation_without_velocity_buffer_fails_before_bfl(self):
+        for speed in (0.0, .002):  # any spec, even zero amplitude
+            backend, calls = self.mocked_backend([circle(actuation(speed))])
+            self.assertIsNotNone(backend._fine_tables["surface_velocity"])
+            del backend._link_tables["surface_velocity"]
+            with self.assertRaisesRegex(RuntimeError, "Surface velocity device buffer is missing"):
+                backend.step(1)
+            names = [name for name, _ in calls]
+            self.assertNotIn("BlockCurvedBoundary", names)
+            self.assertNotIn("BlockStep", names)
+            self.assertEqual(backend._status, "failed")
+            backend.close()
+
+    def test_surface_velocity_helper_matches_dense_soa(self):
+        for mode in ("periodic", "slot_reservoir"):
+            patch, boundary = fine_circle(mode, actuation(-.002))
+            cfg = fine_cfg(mode, patch)
+            soa = curved_links.CurvedLinkSoA()
+            soa.assign_host(*boundary.curved_host)
+            soa.assign_surface(boundary.objects, cfg, boundary.flags)
+            fluid_idx, dirs, _, body_id, rx, ry, _, _ = boundary.curved_host
+            velocity, aperture = curved_links.surface_link_velocity(
+                boundary.objects, cfg, boundary.flags, fluid_idx, dirs, body_id, rx, ry)
+            np.testing.assert_array_equal(velocity, soa.surface_velocity)
+            np.testing.assert_array_equal(aperture, soa.aperture)
+            self.assertEqual(curved_links.surface_link_velocity(
+                (circle(),), cfg, boundary.flags, fluid_idx, dirs, body_id, rx, ry), (None, None))
+            # The aperture must actually cross a block seam in both decomposed layouts.
+            for by, bx in FINE_BLOCKS[1:]:
+                layout = BlockLayout2D.build(*patch.fine_shape, by, bx, periodic=(False, False))
+                _, mapped = blocked_link_tables(layout, boundary.flags, fluid_idx, dirs,
+                                                boundary.curved_host[2], boundary.curved_host[7])
+                L = int(np.prod(layout.local_shape))
+                core = set(layout.owner.ravel().tolist())
+                seam = [int(w) // (9 * L) * L + int(w) % L not in core for w in mapped.write]
+                self.assertTrue(np.any(np.array(seam) & aperture & np.any(velocity != 0, axis=1)),
+                                (mode, by, bx))
 
     def test_tables_map_back_to_dense_slots_in_order(self):
         patch = FixedPatch2D.from_config(config("periodic"), BOUNDS)
@@ -259,6 +328,35 @@ class BlockedFineGPU(unittest.TestCase):
             print(f"circle {mode}: dense vs {len(blocked)} fine-block layouts bit-identical "
                   f"after {STEPS} coarse steps", flush=True)
 
+    def test_actuated_circle_dense_and_blocked_fine_storage_are_bit_identical(self):
+        from CelerisLab.cuda.context import CudaContext
+        for mode in ("periodic", "slot_reservoir"):
+            for speed in (.002, -.002):
+                cfg = config(mode)
+                patch, boundary = fine_circle(mode, actuation(speed))
+                solid = (boundary.flags.reshape(patch.fine_shape) & SOLID) != 0
+                g_coarse, g_fine = state(NY, NX, 1), state(*patch.fine_shape, 2)
+                results = []
+                with CudaContext(0) as context:
+                    for fine_block in (None, *FINE_BLOCKS):
+                        with FixedPatchStepper(cfg, BOUNDS, context=context, coarse_boundary=mode,
+                                               fine_objects=[circle(actuation(speed))],
+                                               fine_block=fine_block) as backend:
+                            backend.initialize(g_coarse, g_fine)
+                            backend.step(STEPS)
+                            results.append((backend.download(), raw(backend)))
+                (dense_public, dense_raw), blocked = results[0], results[1:]
+                fine = dense_raw[1]
+                self.assertTrue(np.isfinite(fine[:, (patch.fine_kind == ACTIVE) & ~solid]).all())
+                self.assertTrue(np.isnan(fine[:, solid]).all())
+                self.assertTrue(np.isnan(fine[:, patch.fine_kind == OUTER_GHOST]).all())
+                for fine_block, (public, raw_arrays) in zip(FINE_BLOCKS, blocked):
+                    msg = f"actuated {speed} {mode} fine_block={fine_block}"
+                    for a, b in zip(public + tuple(raw_arrays), dense_public + tuple(dense_raw)):
+                        np.testing.assert_array_equal(a, b, err_msg=msg)  # values and NaN masks
+                print(f"actuated {speed:+.3f} {mode}: dense vs {len(blocked)} fine-block layouts "
+                      f"bit-identical after {STEPS} coarse steps", flush=True)
+
 
 if __name__ == "__main__":
     unittest.main()