From 2bbf82121ef7070dfa2a5fb234932a0ef8f2ce71 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 15 Oct 2026 04:25:09 +0000 Subject: [PATCH] Blocked fine storage: passive stationary circles via receiver-local BFL FixedPatchStepper(fine_block=..., fine_objects=[circle]) now accepts passive stationary circles; any surface actuation is rejected before CUDA (dense fine path unchanged). prepare_fine_boundary still decides legality on the dense fine grid; its SoA is mapped by the new array-based storage adapter blocks.blocked_link_tables (map_cut_links: unique receiver owner and local writer, row order kept; blocked flags sized nblocks*L), which PeriodicBlockStepper._static_circles now also uses (single-level legality builder unchanged). Per fine hop on the same buffer: halo -> existing BlockCurvedBoundary -> BlockStep(flags, kind). No CUDA changes. Tests: fixture 76 links / 68 q<.5 and 8 seam writes + 8 seam donors per decomposed layout (both modes); adapter flags sizing; mocked launch-order contract; GPU: periodic and slot_reservoir, 8 coarse steps, dense fine body vs 1x1/2x2/4x2 fine blocks, exact values + NaN masks on masked downloads and raw final arrays, ACTIVE fluid finite, solid/OUTER NaN. Co-Authored-By: Claude Opus 5.5 --- src/CelerisLab/lbm/block_stepper.py | 17 ++- src/CelerisLab/lbm/blocks.py | 28 +++++ src/CelerisLab/lbm/refinement_stepper.py | 39 +++++- .../multigrid/test_refinement_blocked_fine.py | 118 +++++++++++++++++- 4 files changed, 181 insertions(+), 21 deletions(-) diff --git a/src/CelerisLab/lbm/block_stepper.py b/src/CelerisLab/lbm/block_stepper.py index 4d6a203..7baf73d 100644 --- a/src/CelerisLab/lbm/block_stepper.py +++ b/src/CelerisLab/lbm/block_stepper.py @@ -19,7 +19,7 @@ import numpy as np import pycuda.driver as cuda from ..cuda import compiler_v2 as compiler -from .blocks import NQ, BlockLayout2D +from .blocks import NQ, BlockLayout2D, blocked_link_tables from .descriptors import D2Q9_EX, D2Q9_EY, FLUID, SOLID MAX_GRID_DIM = 65535 # CUDA gridDim.y/z limit (rows, blocks) @@ -64,15 +64,12 @@ def _static_circles(layout, objects): for q in range(1, NQ) if solid[source(r, q)]} if len(receivers) != len(links) or receivers != required: raise ValueError("Cut links must cover every fluid pull from a solid exactly once") - mapped = layout.map_cut_links(links) # unique receiver owner and local writer - cells = int(np.prod(layout.local_shape)) - blocked_flags = np.zeros(layout.nblocks * cells, np.uint16) # halo entries unused - blocked_flags[layout.owner.ravel().astype(np.intp)] = flags - return dict(solid=solid, flags=blocked_flags, block=mapped.block, - cell=(mapped.receiver % cells).astype(np.uint32), - dir=np.array([link.dir for link in links], np.uint8), - q=np.array([link.q for link in links], np.float32), - fallback=np.array([link.fallback for link in links], np.uint8)) + tables, _ = blocked_link_tables( + layout, flags, np.array([link.fluid_idx for link in links], np.uint32), + np.array([link.dir for link in links], np.uint8), + np.array([link.q for link in links], np.float32), + np.array([link.fallback for link in links], np.uint8)) + return dict(solid=solid, **tables) class PeriodicBlockStepper: diff --git a/src/CelerisLab/lbm/blocks.py b/src/CelerisLab/lbm/blocks.py index 4fc9975..57d8269 100644 --- a/src/CelerisLab/lbm/blocks.py +++ b/src/CelerisLab/lbm/blocks.py @@ -27,6 +27,7 @@ No moving geometry is implied. """ from dataclasses import dataclass from numbers import Integral +from types import SimpleNamespace import numpy as np @@ -266,3 +267,30 @@ class BlockLayout2D: or not blocked.flags.c_contiguous: raise ValueError(f"expected C-contiguous blocked array {self.blocked_shape}") return blocked + + +def blocked_link_tables(layout, flags, fluid_idx, dir, q, fallback): + """Storage adapter: dense static cut-link SoA + dense flags -> blocked tables. + + Legality (geometry, coverage, transfer isolation) stays with the caller's + own builder; this only maps storage via map_cut_links (unique receiver + owner and local writer, input row order kept). Returns (tables, mapped): + tables holds device-ready arrays, with flags sized nblocks*L (unused halo + entries 0) and cell = the receiver's local core cell. + """ + flags = np.asarray(flags) + if flags.shape != (int(np.prod(layout.shape)),) or flags.dtype != np.uint16: + raise ValueError("flags must be dense uint16 with one entry per grid cell") + columns = [np.asarray(fluid_idx), np.asarray(dir), np.asarray(q), np.asarray(fallback)] + if len({len(c) for c in columns}) != 1: + raise ValueError("cut-link columns must have equal length") + mapped = layout.map_cut_links(SimpleNamespace(fluid_idx=r, dir=d, q=v, fallback=f) + for r, d, v, f in zip(*columns)) + cells = int(np.prod(layout.local_shape)) + blocked_flags = np.zeros(layout.nblocks * cells, np.uint16) + blocked_flags[layout.owner.ravel().astype(np.intp)] = flags + tables = dict(flags=blocked_flags, block=mapped.block, + cell=(mapped.receiver % cells).astype(np.uint32), + dir=columns[1].astype(np.uint8), q=columns[2].astype(np.float32), + fallback=columns[3].astype(np.uint8)) + return tables, mapped diff --git a/src/CelerisLab/lbm/refinement_stepper.py b/src/CelerisLab/lbm/refinement_stepper.py index f63eb87..57663ff 100644 --- a/src/CelerisLab/lbm/refinement_stepper.py +++ b/src/CelerisLab/lbm/refinement_stepper.py @@ -14,7 +14,7 @@ import numpy as np import pycuda.driver as cuda from ..cuda import compiler_v2 as compiler -from .blocks import BlockLayout2D +from .blocks import BlockLayout2D, blocked_link_tables from .block_stepper import MAX_GRID_DIM from .refinement import ACTIVE, HALO, OUTER_GHOST, FixedPatch2D, prepare_fine_boundary from .descriptors import SOLID @@ -65,12 +65,15 @@ class FixedPatchStepper: On a device/launch failure the instance becomes unusable: close and recreate it, rather than pretending the partially advanced fine buffer rolled back. - ``fine_block=(by, bx)`` (opt-in, body-free only) stores the fine level in + ``fine_block=(by, bx)`` (opt-in) stores the fine level in blocks.cu storage: the ENTIRE fine_shape, including the four-cell patch ghost band, is partitioned by a non-periodic BlockLayout2D with one-cell 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): + prepare_fine_boundary decides legality on the dense fine grid, then cut + links are staged receiver-local per hop: halo -> BFL -> step, same buffer. """ def __init__(self, cfg, bounds, *, context, coarse_boundary="periodic", fine_objects=(), fine_block=None): @@ -90,8 +93,8 @@ class FixedPatchStepper: fine_objects = tuple(fine_objects) self.fine_layout = None if fine_block is not None: - if fine_objects: - raise ValueError("fine_block storage is body-free; fine_objects need the dense fine level") + 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 @@ -101,6 +104,11 @@ 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 + 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) self._fine_leaf = self.patch.fine_kind == ACTIVE if self.fine_boundary is not None: self._fine_leaf &= (self.fine_boundary.flags.reshape(self.patch.fine_shape) & SOLID) == 0 @@ -113,6 +121,7 @@ class FixedPatchStepper: self._stream = None self._fine_flags = self._fine_curved = self._curved = None self._kind = self._halo_copy = self._halo_dst = self._halo_src = 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) @@ -173,7 +182,12 @@ class FixedPatchStepper: module = self._levels[1].module self._explode = module.get_function("BlockExplode") self._halo_copy = module.get_function("BlockHaloCopy") - if self.fine_boundary is not None: + 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") + elif self.fine_boundary is not None: from .curved_links import CurvedLinkSoA boundary = self.fine_boundary self._fine_flags = self._upload(boundary.flags) @@ -266,6 +280,15 @@ 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 _blocked_fine_boundary(self, buffer, **kw): + if self._fine_tables 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"]) + 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) + def _fine_halo(self, buffer, **kw): pairs = len(self.fine_layout.halo_dst) if pairs: @@ -290,9 +313,12 @@ class FixedPatchStepper: **kw) if self.fine_layout is not None: # Gather one-cell block halos before EACH fine hop. + # Same buffer per hop: halo -> receiver-local BFL -> step. + flags = self._fine_flags if self._fine_flags is not None else np.uintp(0) for src, dst in ((f.current, f.spare), (f.spare, f.current)): self._fine_halo(src, **kw) - f.kernel(src, dst, np.uint32(self.fine_layout.halo), np.uintp(0), + self._blocked_fine_boundary(src, **kw) + f.kernel(src, dst, np.uint32(self.fine_layout.halo), flags, self._kind, grid=f.grid, **kw) else: flags = self._fine_flags if self._fine_flags is not None else np.uintp(0) @@ -355,6 +381,7 @@ class FixedPatchStepper: error = error or exc self._fine_curved = None self._curved = self._fine_flags = None + self._link_tables = {} for pointer in self._allocations: try: pointer.free() diff --git a/tests/validation/multigrid/test_refinement_blocked_fine.py b/tests/validation/multigrid/test_refinement_blocked_fine.py index 6d14e29..c2c0d1a 100644 --- a/tests/validation/multigrid/test_refinement_blocked_fine.py +++ b/tests/validation/multigrid/test_refinement_blocked_fine.py @@ -13,12 +13,14 @@ import numpy as np import pycuda.driver as cuda from CelerisLab.body.geometry.circle import CircleGeometry -from CelerisLab.body.objects import SimObject +from CelerisLab.body.objects import SimObject, SurfaceActuation from CelerisLab.config import LBMConfig from CelerisLab.cuda import compiler_v2 as compiler -from CelerisLab.lbm.blocks import BlockLayout2D +from CelerisLab.lbm.blocks import NO_DONOR, BlockLayout2D, blocked_link_tables from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY -from CelerisLab.lbm.refinement import ACTIVE, INNER_GHOST, OUTER_GHOST, FixedPatch2D +from CelerisLab.lbm.descriptors import SOLID +from CelerisLab.lbm.refinement import (ACTIVE, INNER_GHOST, OUTER_GHOST, FixedPatch2D, + prepare_fine_boundary) from CelerisLab.lbm.refinement_stepper import FixedPatchStepper, blocked_fine_tables NX, NY, BOUNDS, STEPS = 32, 24, (8, 6, 24, 18), 8 @@ -45,6 +47,27 @@ def state(ny, nx, seed): return (feq + noise).astype(np.float32) +def circle(actuation=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) + + +def fine_circle(mode): + cfg = config(mode) + patch = FixedPatch2D.from_config(cfg, BOUNDS) + return patch, prepare_fine_boundary(cfg, patch, (circle(),), wall_aware=mode == "slot_reservoir") + + +def seam_cells(layout, mapped): + """Distinct non-core cells written / read as q<.5 donors (receiver-local halos).""" + L = int(np.prod(layout.local_shape)) + core = set(layout.owner.ravel().tolist()) + writes = {int(w) // (9 * L) * L + int(w) % L for w in mapped.write} + donors = {int(d) for d in mapped.donor if d != NO_DONOR} + return len(writes - core), len(donors - core) + + def raw(backend): """Unmasked device arrays; blocked fine storage unpacked from owner cores.""" coarse, fine = backend._levels @@ -63,9 +86,9 @@ class BlockedFineCPU(unittest.TestCase): return Mock(activate=Mock(side_effect=AssertionError("CUDA reached")), sm_arch="sm_70") def test_rejections_before_cuda(self): - circle = SimObject(0, geometry=CircleGeometry(16.0, 12.0, 2.0), center=(16.0, 12.0), radius=2.0) + 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": [circle]}): + {"fine_block": (16, 20), "fine_objects": [actuated]}): with self.subTest(kw=kw), self.assertRaises(ValueError): FixedPatchStepper(config("periodic"), BOUNDS, context=self.forbidden_context(), **kw) @@ -94,6 +117,64 @@ class BlockedFineCPU(unittest.TestCase): backend.close() self.assertTrue(sizes) + def test_circle_fixture_counts_and_seam_staging(self): + for mode in ("periodic", "slot_reservoir"): + patch, boundary = fine_circle(mode) + 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:]: + 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)) + + def test_storage_adapter_flags_sizing(self): + patch, boundary = fine_circle("periodic") + columns = [boundary.curved_host[i] for i in (0, 1, 2, 7)] + layout = BlockLayout2D.build(*patch.fine_shape, 16, 20, periodic=(False, False)) + tables, mapped = blocked_link_tables(layout, boundary.flags, *columns) + L = int(np.prod(layout.local_shape)) + self.assertEqual(tables["flags"].shape, (layout.nblocks * L,)) + np.testing.assert_array_equal(tables["flags"][layout.owner.ravel().astype(np.intp)], + boundary.flags) + np.testing.assert_array_equal(tables["cell"], mapped.receiver % np.uint64(L)) + for bad in (boundary.flags[:-1], boundary.flags.astype(np.uint32)): + with self.assertRaises(ValueError): + blocked_link_tables(layout, bad, *columns) + + def test_same_buffer_halo_bfl_step_launch_order(self): + context = Mock(activate=nullcontext, sm_arch="sm_70") + calls, functions = [], {} + + def get_function(name): + return functions.setdefault(name, Mock(side_effect=lambda *a, **k: calls.append((name, a)))) + module = Mock(get_function=Mock(side_effect=get_function)) + 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"): + 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) + c, f = backend._levels + coarse, cur, spare = c.current, f.current, f.spare + flags, kind = backend._fine_flags, backend._kind + 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) + def test_tables_map_back_to_dense_slots_in_order(self): patch = FixedPatch2D.from_config(config("periodic"), BOUNDS) fy, fx = patch.fine_shape @@ -151,6 +232,33 @@ class BlockedFineGPU(unittest.TestCase): print(f"{mode}: dense vs {len(blocked)} fine-block layouts bit-identical " f"after {STEPS} coarse steps", flush=True) + def test_passive_circle_dense_and_blocked_fine_storage_are_bit_identical(self): + from CelerisLab.cuda.context import CudaContext + for mode in ("periodic", "slot_reservoir"): + cfg = config(mode) + patch, boundary = fine_circle(mode) + 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()], 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(), mode) + self.assertTrue(np.isnan(fine[:, solid]).all(), mode) + self.assertTrue(np.isnan(fine[:, patch.fine_kind == OUTER_GHOST]).all(), mode) + for fine_block, (public, raw_arrays) in zip(FINE_BLOCKS, blocked): + msg = f"circle {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"circle {mode}: dense vs {len(blocked)} fine-block layouts bit-identical " + f"after {STEPS} coarse steps", flush=True) + if __name__ == "__main__": unittest.main()