diff --git a/src/CelerisLab/lbm/block_stepper.py b/src/CelerisLab/lbm/block_stepper.py index d687426..f9c4c39 100644 --- a/src/CelerisLab/lbm/block_stepper.py +++ b/src/CelerisLab/lbm/block_stepper.py @@ -18,61 +18,13 @@ import numpy as np import pycuda.driver as cuda from ..cuda import compiler -from .blocks import (NQ, BlockLayout2D, blocked_link_tables, launch_curved_boundary, - launch_halo_copy, level_config) -from .descriptors import D2Q9_EX, D2Q9_EY, FLUID, SOLID +from .blocks import (NQ, BlockLayout2D, launch_curved_boundary, launch_halo_copy, level_config, + static_circle_tables) from .population_storage import _StoredPopulations, _encode_physical, _decode_stored MAX_GRID_DIM = 65535 # CUDA gridDim.y/z limit (rows, blocks) -def _static_circles(layout, objects): - """Host tables for stationary circles; rejects before any CUDA work. - - Every periodic neighbor of every fluid core is inspected, including domain - edges that CircleGeometry does not enumerate, so clipped or incomplete - geometry cannot pass. Each physical link has one receiver-owning block. - """ - from ..body.geometry.circle import CircleGeometry - - ny, nx = layout.shape - flags = np.full(ny * nx, FLUID, np.uint16) - for obj in objects: - if (type(obj.geometry) is not CircleGeometry or obj.is_sensor or obj.is_force_region - or obj.surface_actuation is not None - or any(getattr(obj.state, name) != 0 for name in ("vx", "vy", "omega")) - or any(vars(obj.control).values())): - raise ValueError("Block bodies must be stationary passive circles") - cx, cy = obj.geometry.center - radius = obj.geometry._radius - if not np.isfinite([cx, cy, radius]).all() or radius <= 0: - raise ValueError("Block circle center/radius must be finite with radius > 0") - mask = obj.geometry.build_flag_mask(nx, ny) - solid = (mask & SOLID) != 0 - if not solid.any() or np.any(flags[solid] & SOLID): - raise ValueError("Block circles must be resolved and disjoint") - flags[solid] = mask[solid] - solid = ((flags & SOLID) != 0).reshape(ny, nx) - if solid[[0, -1]].any() or solid[:, [0, -1]].any(): - raise ValueError("Block circles must not touch the domain edge") - links = [link for obj in objects - for link in obj.geometry.build_cut_links(nx, ny, domain_flags=flags)] - def source(r, q): - y, x = divmod(r, nx) - return (y + D2Q9_EY[q]) % ny, (x + D2Q9_EX[q]) % nx - receivers = {(int(link.fluid_idx), int(link.dir)) for link in links} - required = {(int(r), q) for r in np.flatnonzero(~solid.ravel()) - 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") - 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: """G->G stepping of Q-major (9, ny, nx) postcollision populations. @@ -96,7 +48,7 @@ class PeriodicBlockStepper: if layout.nblocks > MAX_GRID_DIM or ly > MAX_GRID_DIM or not 0 < pairs < 2**31: raise ValueError("block count, block rows or halo pairs exceed launch limits") objects = tuple(objects) # an empty iterable takes the body-free path - self.boundary = _static_circles(layout, objects) if objects else None + self.boundary = static_circle_tables(layout, objects) if objects else None links = 0 if self.boundary is None else len(self.boundary["cell"]) if links >= 2**31: raise ValueError("cut links exceed launch limits") diff --git a/src/CelerisLab/lbm/blocks.py b/src/CelerisLab/lbm/blocks.py index 3626926..6a9c7bb 100644 --- a/src/CelerisLab/lbm/blocks.py +++ b/src/CelerisLab/lbm/blocks.py @@ -407,3 +407,51 @@ def launch_curved_boundary(kernel, buffer, tables, count, *, block, stream, np.uint32(count), none if velocity is None else velocity, none if rx is None else rx, none if ry is None else ry, none if raw is None else raw, grid=((count + block[0] - 1) // block[0], 1), block=block, stream=stream) + + +def static_circle_tables(layout, objects): + """Host tables for stationary circles; rejects before any CUDA work. + + Every periodic neighbor of every fluid core is inspected, including domain + edges that CircleGeometry does not enumerate, so clipped or incomplete + geometry cannot pass. Each physical link has one receiver-owning block. + """ + from ..body.geometry.circle import CircleGeometry + from .descriptors import FLUID, SOLID + + ny, nx = layout.shape + flags = np.full(ny * nx, FLUID, np.uint16) + for obj in objects: + if (type(obj.geometry) is not CircleGeometry or obj.is_sensor or obj.is_force_region + or obj.surface_actuation is not None + or any(getattr(obj.state, name) != 0 for name in ("vx", "vy", "omega")) + or any(vars(obj.control).values())): + raise ValueError("Block bodies must be stationary passive circles") + cx, cy = obj.geometry.center + radius = obj.geometry._radius + if not np.isfinite([cx, cy, radius]).all() or radius <= 0: + raise ValueError("Block circle center/radius must be finite with radius > 0") + mask = obj.geometry.build_flag_mask(nx, ny) + solid = (mask & SOLID) != 0 + if not solid.any() or np.any(flags[solid] & SOLID): + raise ValueError("Block circles must be resolved and disjoint") + flags[solid] = mask[solid] + solid = ((flags & SOLID) != 0).reshape(ny, nx) + if solid[[0, -1]].any() or solid[:, [0, -1]].any(): + raise ValueError("Block circles must not touch the domain edge") + links = [link for obj in objects + for link in obj.geometry.build_cut_links(nx, ny, domain_flags=flags)] + def source(r, q): + y, x = divmod(r, nx) + return (y + D2Q9_EY[q]) % ny, (x + D2Q9_EX[q]) % nx + receivers = {(int(link.fluid_idx), int(link.dir)) for link in links} + required = {(int(r), q) for r in np.flatnonzero(~solid.ravel()) + 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") + 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) diff --git a/src/CelerisLab/lbm/refinement_stepper.py b/src/CelerisLab/lbm/refinement_stepper.py index 3c23880..8c9c3c9 100644 --- a/src/CelerisLab/lbm/refinement_stepper.py +++ b/src/CelerisLab/lbm/refinement_stepper.py @@ -15,7 +15,8 @@ import numpy as np import pycuda.driver as cuda from ..cuda import compiler -from .blocks import BlockLayout2D, launch_curved_boundary, launch_halo_copy, level_config +from .blocks import (BlockLayout2D, launch_curved_boundary, launch_halo_copy, level_config, + static_circle_tables) from .channel import channel_masks, validate_channel_config, validate_channel_transfers from ..body.jets import jet_link_velocity from .block_stepper import MAX_GRID_DIM @@ -23,6 +24,7 @@ from .hierarchy import build_hierarchy, is_rectangle from .refinement import ACTIVE, HALO, FixedPatch2D, FixedPatchSet2D, prepare_fine_boundary from .refinement_tables import (blocked_base_tables, blocked_fine_tables, # noqa: F401 (re-export) blocked_kind, projected_explosion_tables, _prepare_blocked_boundary, + _BlockedBoundary, BLOCKED_BOUNDARY_DEVICE_COLUMNS) from .descriptors import SOLID from .raw_exchange import summarize_exchange @@ -68,14 +70,16 @@ def _level_name(level): def _check_levels(refinements, blocks, objects, jets): """Candidate limits of the level lists, checked before any device work.""" depth = len(refinements) + 1 - if depth < 2: - raise ValueError("at least one refinement level is required") for name, values in (("level_blocks", blocks), ("level_objects", objects), ("level_jets", jets)): if len(values) != depth: raise ValueError(f"{name} needs one entry per level ({depth})") for level, block in enumerate(blocks): if block is not None and len(block) != 2: raise ValueError(f"level_blocks[{level}] must be (by, bx)") + if depth == 1: # one blocked level: the root carries the static circles + if blocks[0] is None or jets[0]: + raise ValueError("a single level must be blocked and has no Jets") + return if objects[0] or jets[0]: raise ValueError("base-level bodies/Jets are not supported") if any(blocks[level - 1] is not None and blocks[level] is None for level in range(1, depth)): @@ -98,6 +102,9 @@ class FixedPatchStepper: level-wide grid, public fine arrays (2*ny, 2*nx), NaN outside ACTIVE cells); see lbm.hierarchy. ``level_blocks``, ``level_objects`` and ``level_jets`` hold one entry per level, L0 included (L0 bodies/Jets are not supported). + ``refinements=[]`` is one blocked level: periodic SRT root without LES, + Jets or exchange observation; its ``level_objects[0]`` are stationary + passive circles (static BFL, as PeriodicBlockStepper). ``coarse_boundary`` is ``periodic`` (default) or ``slot_reservoir``: a bottom on-site slot and top fixed-density reservoir from ``cfg``. @@ -196,6 +203,10 @@ class FixedPatchStepper: self._ddf_shifting = bool(cfg.ddf_shifting) if self._ddf_shifting and any(b is None for b in level_blocks): raise ValueError("shifted populations require all-blocked levels; dense/mixed is unsupported") + if depth == 1 and (coarse_boundary != "periodic" or cfg.collision != "SRT" or cfg.les_enabled + or observe_raw_exchange or observe_window_exchange): + raise ValueError("a single level requires a periodic SRT root without LES or exchange " + "observation") self._channel = coarse_boundary == "channel" if self._channel: validate_channel_config(cfg, level_blocks) @@ -229,12 +240,14 @@ class FixedPatchStepper: base = blocked_base_tables(pairs[0], self._layouts[0], np.argsort(pairs[0].coalesce_coarse)) validate_channel_transfers(pairs[0], self._layouts[0], base[1]) # Kernel configurations per level (geometry configs live in self.hierarchy). - fine_cfg = replace(cfg, nx=pairs[0].fine_shape[1], ny=pairs[0].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._level_cfgs = [replace(cfg), fine_cfg] + self._level_cfgs = [replace(cfg)] + if pairs: + fine_cfg = replace(cfg, nx=pairs[0].fine_shape[1], ny=pairs[0].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._level_cfgs.append(fine_cfg) for pair in pairs[1:]: parent = self._level_cfgs[-1] self._level_cfgs.append(replace(parent, nx=pair.fine_shape[1], ny=pair.fine_shape[0], @@ -275,12 +288,18 @@ class FixedPatchStepper: self._exchange_centers = tuple(tuple(tuple(o.geometry.center) for o in objects) for objects in level_objects) self._level_jets = level_jets - self._level_boundaries = [None] + [ + root = None + if depth == 1 and level_objects[0]: + tables = static_circle_tables(self._layouts[0], level_objects[0]) + root = _BlockedBoundary(None, tables, (), False) + self._level_boundaries = [root] + [ _prepare_blocked_boundary(self._boundary_hosts[level], self._layouts[level], self._level_cfgs[level], level_jets[level], _shifted_blocked=self._ddf_shifting) for level in range(1, depth)] self._leaf = [kind == ACTIVE for kind in self.hierarchy.kinds] + if root is not None: + self._leaf[0] &= ~root.tables["solid"] for level in range(1, depth): host = self._boundary_hosts[level] if host is not None: @@ -324,7 +343,9 @@ class FixedPatchStepper: b.window_integral = self._allocate(count * 4 * 8) b.accumulate = level.module.get_function("BlockAccumulateExchange") b.kernel = level.module.get_function("BlockCurvedBoundary") - if self._level_boundaries[1] is not None: + if depth == 1: + pass + elif self._level_boundaries[1] is not None: b = self._level_boundaries[1] self._fine_flags, self._link_tables, self._curved = b.flags, b.device, b.kernel elif self._boundary_hosts[1] is not None: @@ -350,7 +371,8 @@ class FixedPatchStepper: """Compile one level's step module and allocate its two G buffers.""" layout, config = self._layouts[index], self._level_cfgs[index] pairs = self.hierarchy.pairs - tau = pairs[0].tau_coarse if index == 0 else pairs[index - 1].tau_fine + tau = (pairs[0].tau_coarse if pairs else 3 * config.viscosity + .5) if index == 0 \ + else pairs[index - 1].tau_fine name, source = ("PatchCoarseStep" if index == 0 else "PatchFineStep"), "refinement.cu" if layout is not None: config = level_config(config, layout, tau) @@ -732,7 +754,8 @@ class FixedPatchStepper: np.float64(h), np.float64(h), grid=((count + self._block[0] - 1) // self._block[0], 1), **kw) flags = np.uintp(0) if level.boundary is None else level.boundary.flags - level.kernel(src, dst, np.uint32(level.layout.halo), flags, level.kind, grid=level.grid, **kw) + kind = np.uintp(0) if level.kind is None else level.kind # single level: all ACTIVE + level.kernel(src, dst, np.uint32(level.layout.halo), flags, kind, grid=level.grid, **kw) def _gather_halo(self, level, buffer, **kw): pairs = len(level.layout.halo_dst) diff --git a/tests/validation/multigrid/test_block_stepper.py b/tests/validation/multigrid/test_block_stepper.py index 12939d2..648e796 100644 --- a/tests/validation/multigrid/test_block_stepper.py +++ b/tests/validation/multigrid/test_block_stepper.py @@ -25,6 +25,7 @@ from CelerisLab.cuda import compiler from CelerisLab.lbm import block_stepper as implementation from CelerisLab.lbm.blocks import NO_DONOR, BlockLayout2D from CelerisLab.lbm.block_stepper import PeriodicBlockStepper +from CelerisLab.lbm.refinement_stepper import FixedPatchStepper from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY ROOT = Path(__file__).resolve().parents[3] @@ -118,6 +119,23 @@ class CpuContracts(unittest.TestCase): def forbidden_context(self): return Mock(activate=Mock(side_effect=AssertionError("CUDA reached")), sm_arch="sm_70") + def test_single_level_fixed_patch_contract_before_cuda(self): + def single(cfg=None, **kw): + options = dict(context=self.forbidden_context(), refinements=[], level_blocks=[(12, 16)]) + options.update(kw) + return FixedPatchStepper(cfg or config(), **options) + for cfg, kw in ((None, dict(level_blocks=[None])), + (None, dict(coarse_boundary="channel")), + (None, dict(coarse_boundary="slot_reservoir")), + (config(collision="MRT"), {}), + (None, dict(observe_raw_exchange=True)), + (None, dict(level_jets=[(object(),)])), + (None, dict(level_objects=[(circle((15.5, 4.0), 4.3),)]))): # touches edge + with self.subTest(kw=kw), self.assertRaises(ValueError): + single(cfg, **kw) + with self.assertRaisesRegex(AssertionError, "CUDA reached"): + single(level_objects=[(circle(),)]) + def test_unsupported_inputs_rejected_before_cuda(self): layout = BlockLayout2D.build(NY, NX, 12, 16) cases = ((config(), BlockLayout2D.build(NY, NX, 12, 16, periodic=(True, False))), @@ -139,7 +157,7 @@ class CpuContracts(unittest.TestCase): with self.subTest(center=center, radius=radius), self.assertRaises(ValueError): PeriodicBlockStepper(config(), layout, context=self.forbidden_context(), objects=[circle(center, radius)]) - with patch.object(implementation, "_static_circles") as prepare: + with patch.object(implementation, "static_circle_tables") as prepare: with self.assertRaisesRegex(AssertionError, "CUDA reached"): # body-free path PeriodicBlockStepper(config(), layout, context=self.forbidden_context(), objects=iter(())) prepare.assert_not_called() @@ -177,7 +195,7 @@ class CpuContracts(unittest.TestCase): def test_static_circle_tables(self): layout = BlockLayout2D.build(NY, NX, 12, 16) - tables = implementation._static_circles(layout, (circle(),)) + tables = implementation.static_circle_tables(layout, (circle(),)) links = CircleGeometry(*CENTER, RADIUS).build_cut_links(NX, NY, domain_flags=dense_flags()) self.assertEqual(len(tables["cell"]), len(links)) np.testing.assert_array_equal(tables["q"], np.array([l.q for l in links], np.float32)) @@ -331,6 +349,29 @@ class PeriodicBlockGPU(unittest.TestCase): *(a - b for a, b in zip(after, before))), flush=True) self.assertLessEqual(error, ORACLE_ATOL) + def test_single_level_fixed_patch_is_bit_identical(self): + """FixedPatchStepper(refinements=[]) reproduces PeriodicBlockStepper bit for bit.""" + from CelerisLab.cuda.context import CudaContext + g0 = initial_state() + with CudaContext(0) as context: + for (by, bx), objects in [(split, objects) for split in SPLITS + for objects in ((), (circle(),))]: + with self.subTest(split=(by, bx), bodies=len(objects)): + layout = BlockLayout2D.build(NY, NX, by, bx) + with PeriodicBlockStepper(config(), layout, context=context, + objects=objects) as backend: + backend.initialize(g0) + backend.step(STEPS) + expected = backend.download() + with FixedPatchStepper(config(), context=context, refinements=[], + level_blocks=[(by, bx)], + level_objects=[objects]) as backend: + backend.initialize(g0) + backend.step(STEPS) + self.assertEqual(backend.step_count, STEPS) + (actual,) = backend.download() + np.testing.assert_array_equal(actual, expected) + if __name__ == "__main__": unittest.main()