From 7bf6101460921f78fef1c411a7a1deb3d67f59a8 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 15 Oct 2026 04:45:17 +0000 Subject: [PATCH] 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 --- src/CelerisLab/lbm/block_stepper.py | 3 +- src/CelerisLab/lbm/curved_links.py | 47 +++--- src/CelerisLab/lbm/kernels/blocks.cu | 12 +- src/CelerisLab/lbm/refinement_stepper.py | 23 ++- .../multigrid/test_refinement_blocked_fine.py | 146 +++++++++++++++--- 5 files changed, 181 insertions(+), 50 deletions(-) diff --git a/src/CelerisLab/lbm/block_stepper.py b/src/CelerisLab/lbm/block_stepper.py index 7baf73d..384d396 100644 --- a/src/CelerisLab/lbm/block_stepper.py +++ b/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) diff --git a/src/CelerisLab/lbm/curved_links.py b/src/CelerisLab/lbm/curved_links.py index 139de98..7f15e89 100644 --- a/src/CelerisLab/lbm/curved_links.py +++ b/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. diff --git a/src/CelerisLab/lbm/kernels/blocks.cu b/src/CelerisLab/lbm/kernels/blocks.cu index 543f07d..c504a34 100644 --- a/src/CelerisLab/lbm/kernels/blocks.cu +++ b/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. diff --git a/src/CelerisLab/lbm/refinement_stepper.py b/src/CelerisLab/lbm/refinement_stepper.py index 57663ff..58ec5d3 100644 --- a/src/CelerisLab/lbm/refinement_stepper.py +++ b/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) diff --git a/tests/validation/multigrid/test_refinement_blocked_fine.py b/tests/validation/multigrid/test_refinement_blocked_fine.py index c2c0d1a..95876a8 100644 --- a/tests/validation/multigrid/test_refinement_blocked_fine.py +++ b/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()