From 76555de4b894033d3e520ae2fd2d3e6068833442 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 16 Oct 2026 14:35:57 +0000 Subject: [PATCH] Three levels: nested body-free child patch on a single chain FixedPatchStepper(..., child=(x0,y0,x1,y1), child_block=(by,bx)) adds an L2 patch in L1 window cells (requires fine_block/base_block, one L1 patch, no bodies/jets). nest_patch() builds the L1->L2 FixedPatch2D (2*nu on L1) and the composed L1 kind (inherited INNER/OUTER kept, child-covered ACTIVE made OUTER), rejecting from the actual tables: coverage outside ACTIVE, donors or receivers not uncovered ACTIVE, surviving ACTIVE pulls of covered/OUTER cells without a child-supplied arrival, and overlap with L1's restriction footprint. Step order per L1 hop: explode L1(src)->L2, two L2 hops (halo, step), L1 halo, receiver-local L2->L1 staging into src, L1 step. The second hop thus feeds L2 from the new L1 G; L2 is never refilled from L0; L0 stays once per step (1:2:4). Existing kernels only; a small _Pair holds the child transfer bindings (explode from the child module, parent window stride). Two-level API unchanged. Co-Authored-By: Claude Opus 5.5 --- src/CelerisLab/lbm/refinement.py | 45 +++++++ src/CelerisLab/lbm/refinement_stepper.py | 124 +++++++++++++++--- tests/unit/test_refinement.py | 27 ++++ .../multigrid/test_refinement_three_level.py | 89 +++++++++++++ 4 files changed, 264 insertions(+), 21 deletions(-) create mode 100644 tests/validation/multigrid/test_refinement_three_level.py diff --git a/src/CelerisLab/lbm/refinement.py b/src/CelerisLab/lbm/refinement.py index 713960d..43ea6c8 100644 --- a/src/CelerisLab/lbm/refinement.py +++ b/src/CelerisLab/lbm/refinement.py @@ -299,6 +299,51 @@ def _transfer_footprint(patch): return fine, coarse, _readonly(sorted(reads), np.int64) +def nest_patch(cfg, parent, child_bounds): + """Child FixedPatch2D inside a fine patch window, plus the composed middle kind. + + ``cfg`` is the base (L0) configuration and ``parent`` the L0->L1 + FixedPatch2D; ``child_bounds`` (x0, y0, x1, y1) are in L1 window cells + (its dense fine grid, ghost band included). Returns (child, kind): the + L1->L2 FixedPatch2D (viscosity 2*nu on L1) and the L1 kind array with + inherited INNER/OUTER kept and only the child-covered ACTIVE cells made + OUTER (invalid). Rejects, from the actual tables before any allocation: + coverage outside inherited ACTIVE; explosion donors or coalescence + receivers that are not uncovered ACTIVE; any surviving ACTIVE pull that + would read a covered or OUTER cell without a child-supplied arrival; and + child coverage of L1's own two-hop restriction footprint toward L0. + """ + fy, fx = parent.fine_shape + middle = replace(cfg, nx=fx, ny=fy, viscosity=2 * cfg.viscosity, y_wall_bc="bounce_back", + top_wall_bc="velocity", wall_slot_width=0, wall_normal_velocity=0.0, + velocity=0.0) + child = FixedPatch2D.from_config(middle, child_bounds) + kind = parent.fine_kind + covered = child.covered + active = (kind == ACTIVE) & ~covered + if not (kind[covered] == ACTIVE).all(): + raise ValueError("child must cover only inherited ACTIVE cells") + cells = fx * fy + receivers = child.coalesce_coarse % np.uint64(cells) + if (not active.ravel()[child.ghost_coarse].all() + or not active.ravel()[receivers.astype(np.intp)].all()): + raise ValueError("child donors/receivers must be uncovered parent ACTIVE cells") + supplied = set(map(int, child.coalesce_coarse)) + invalid = covered | (kind == OUTER_GHOST) + for y, x in zip(*np.nonzero(active)): + for q in range(1, 9): + sy, sx = y - D2Q9_EY[q], x - D2Q9_EX[q] + inside = 0 <= sy < fy and 0 <= sx < fx + if (not inside or invalid[sy, sx]) and q * cells + y * fx + x not in supplied: + raise ValueError("a surviving parent ACTIVE pull would read an invalid cell") + footprint, _, _ = _transfer_footprint(parent) + if covered.ravel()[sorted(footprint)].any(): + raise ValueError("child covers the parent's restriction footprint") + composed = kind.copy() + composed[covered] = OUTER_GHOST + return child, _readonly(composed, np.uint8) + + def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=()): """Project fixed circles and reject boundary/transfer overlap before CUDA. diff --git a/src/CelerisLab/lbm/refinement_stepper.py b/src/CelerisLab/lbm/refinement_stepper.py index e245e50..bfa397c 100644 --- a/src/CelerisLab/lbm/refinement_stepper.py +++ b/src/CelerisLab/lbm/refinement_stepper.py @@ -19,7 +19,7 @@ from .curved_links import surface_link_velocity from ..body.jets import jet_link_shapes, jet_link_velocity from .block_stepper import MAX_GRID_DIM from .refinement import (ACTIVE, HALO, OUTER_GHOST, FixedPatch2D, FixedPatchSet2D, - prepare_fine_boundary) + nest_patch, prepare_fine_boundary) from .descriptors import SOLID, D2Q9_EX, D2Q9_EY @@ -31,15 +31,36 @@ def blocked_fine_tables(patch, layout, order): unused halo entries OUTER). Coarse donors/targets are unchanged. """ owner = layout.owner.ravel().astype(np.intp) - owned = layout.owned.ravel() # sparse multipatch levels leave holes unstored cells = patch.fine_shape[0] * patch.fine_shape[1] q, cell = np.divmod(patch.coalesce_fine[order], np.uint64(cells)) - kind = np.full(layout.nblocks * int(np.prod(layout.local_shape)), OUTER_GHOST, np.uint8) - kind[owner[owned]] = patch.fine_kind.ravel()[owned] + kind = blocked_kind(layout, patch.fine_kind) return (owner[patch.ghost_fine.astype(np.intp)].astype(np.uint64), layout.slot(owner[cell.astype(np.intp)], q), kind) +def blocked_kind(layout, dense_kind): + """Per-block-cell kind; unused halo entries and unstored holes stay OUTER.""" + owner = layout.owner.ravel().astype(np.intp) + owned = layout.owned.ravel() # sparse multipatch levels leave holes unstored + kind = np.full(layout.nblocks * int(np.prod(layout.local_shape)), OUTER_GHOST, np.uint8) + kind[owner[owned]] = np.asarray(dense_kind).ravel()[owned] + return kind + + +@dataclass +class _Pair: + """Device transfer tables of one nested parent->child relationship.""" + explode: object # BlockExplode from the CHILD module (child TOTAL_CELLS) + coalesce: object + ghosts: object + donors: object + targets: object # receiver-local staged parent slots + children: object + n_ghosts: int + n_rows: int + parent_cells: int # parent block window cells (explode donor stride) + + def blocked_base_tables(patch, layout, order, wall_aware=False): """Parent donor cells and receiver-local arrival slots, in transfer row order. @@ -119,7 +140,7 @@ class FixedPatchStepper: blocked one above, with every patch's transfers in one launch each. """ def __init__(self, cfg, bounds, *, context, coarse_boundary="periodic", fine_objects=(), - fine_block=None, base_block=None, jets=()): + fine_block=None, base_block=None, jets=(), child=None, child_block=None): if HALO != 4: raise RuntimeError("refinement.cu requires a four-cell fine halo") if coarse_boundary not in ("periodic", "slot_reservoir"): @@ -165,6 +186,22 @@ class FixedPatchStepper: if (self.fine_layout.nblocks > MAX_GRID_DIM or self.fine_layout.local_shape[0] > MAX_GRID_DIM): raise ValueError("fine blocks exceed launch limits") + self.child_patch = self.child_layout = None + self._middle_kind = None + if child is not None: + if multipatch or fine_block is None or base_block is None or child_block is None: + raise ValueError("a nested child needs one patch plus fine_block, base_block " + "and child_block") + if tuple(fine_objects) or tuple(jets): + raise ValueError("a nested child does not support fine_objects or jets yet") + if len(child_block) != 2: + raise ValueError("child_block must be (by, bx)") + self.child_patch, self._middle_kind = nest_patch(cfg, self.patch, child) + self.child_layout = BlockLayout2D.build(*self.child_patch.fine_shape, *child_block, + periodic=(False, False)) + if (self.child_layout.nblocks > MAX_GRID_DIM + or self.child_layout.local_shape[0] > MAX_GRID_DIM): + raise ValueError("child blocks exceed launch limits") self._jets = tuple(jets) self.fine_boundary = prepare_fine_boundary( cfg, self.patch, fine_objects, wall_aware=self._wall_aware, jets=self._jets) @@ -173,7 +210,11 @@ class FixedPatchStepper: 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 + self._fine_leaf = (self.patch.fine_kind if self._middle_kind is None + else self._middle_kind) == ACTIVE + self._child_leaf = (None if self.child_patch is None + else self.child_patch.fine_kind == ACTIVE) + self._child = None if self.fine_boundary is not None: self._fine_leaf &= (self.fine_boundary.flags.reshape(self.patch.fine_shape) & SOLID) == 0 self.ctx = context @@ -208,9 +249,12 @@ class FixedPatchStepper: with self.ctx.activate(): try: self._stream = cuda.Stream() - for index, (config, layout, tau, name) in enumerate(( - (replace(cfg), self.base_layout, self.patch.tau_coarse, "PatchCoarseStep"), - (fine_cfg, self.fine_layout, self.patch.tau_fine, "PatchFineStep"))): + specs = [(replace(cfg), self.base_layout, self.patch.tau_coarse, "PatchCoarseStep"), + (fine_cfg, self.fine_layout, self.patch.tau_fine, "PatchFineStep")] + if self.child_patch is not None: + specs.append((replace(fine_cfg, viscosity=4 * cfg.viscosity), self.child_layout, + self.child_patch.tau_fine, "PatchFineStep")) + for index, (config, layout, tau, name) in enumerate(specs): source = "refinement.cu" if layout is not None: config = replace(config, nx=layout.local_shape[1], ny=layout.local_shape[0], @@ -233,7 +277,8 @@ class FixedPatchStepper: if size != struct.calcsize("=ffffI"): raise RuntimeError("LBMParams host/device layout mismatch") cuda.memcpy_htod(ptr, struct.pack("=ffffI", config.omega, 0, 0, 0, n_objects)) - shape = (9, *(self.patch.coarse_shape if index == 0 else self.patch.fine_shape)) + shape = (9, *(self.patch.coarse_shape, self.patch.fine_shape, + self.child_patch and self.child_patch.fine_shape)[index]) size = int(np.prod(shape if layout is None else layout.blocked_shape)) * 4 grid = ((config.nx + self._block[0] - 1) // self._block[0], config.ny, 1 if layout is None else layout.nblocks) @@ -252,7 +297,7 @@ class FixedPatchStepper: cuda.memcpy_htod(ptr, struct.pack("=IIIf", cfg.nx, cfg.ny, cfg.wall_slot_width, cfg.wall_normal_velocity)) p = self.patch - c, f = self._levels + c, f = self._levels[:2] order = np.argsort(p.coalesce_coarse) fine_tables = (None if f.layout is None else blocked_fine_tables(p, f.layout, order)) @@ -268,10 +313,23 @@ class FixedPatchStepper: self._explode = f.module.get_function("PatchExplode" if fine_tables is None else "BlockExplode") if fine_tables is not None: - f.kind = self._upload(fine_tables[2]) + f.kind = self._upload(fine_tables[2] if self._middle_kind is None + else blocked_kind(f.layout, self._middle_kind)) if base_tables is not None: c.kind = self._upload(base_tables[2]) self._coalesce = c.module.get_function("BlockCoalesce") + if self.child_patch is not None: + g, cp = self._levels[2], self.child_patch + order = np.argsort(cp.coalesce_coarse) + ghosts, children, g_kind = blocked_fine_tables(cp, g.layout, order) + donors, targets, _ = blocked_base_tables(cp, f.layout, order) + g.kind = self._upload(g_kind) + g.halo_copy = g.module.get_function("BlockHaloCopy") + self._child = _Pair(g.module.get_function("BlockExplode"), + f.module.get_function("BlockCoalesce"), + self._upload(ghosts), self._upload(donors.astype(np.uint64)), + self._upload(targets), self._upload(children), + len(ghosts), len(targets), int(np.prod(f.layout.local_shape))) 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() @@ -314,24 +372,28 @@ class FixedPatchStepper: def step_count(self): return self._step_count - def _mask_nonleaves(self, coarse, fine): + def _mask_nonleaves(self, coarse, fine, child=None): coarse[:, self.patch.covered] = np.nan if self._wall_aware: coarse[:, (0, -1), :] = np.nan fine[:, ~self._fine_leaf] = np.nan + if child is not None: + child[:, ~self._child_leaf] = np.nan - def initialize(self, coarse_g, fine_g): - """Upload an explicit G_0; invalid non-leaves are discarded.""" + def initialize(self, coarse_g, fine_g, child_g=None): + """Upload an explicit G_0 (child_g only with a nested child); non-leaves are discarded.""" with self._lock, self.ctx.activate(): self._require("new", "ready") - arrays = [np.array(a, dtype=np.float32, order="C", copy=True) - for a in (coarse_g, fine_g)] + given = (coarse_g, fine_g) + (() if child_g is None else (child_g,)) + if len(given) != len(self._levels): + raise ValueError("Expected one Q-major array per level") + arrays = [np.array(a, dtype=np.float32, order="C", copy=True) for a in given] if any(a.shape != level.shape for a, level in zip(arrays, self._levels)): - raise ValueError("Expected Q-major (9,ny,nx) arrays for both levels") + raise ValueError("Expected Q-major (9,ny,nx) arrays for every level") leaf = ~self.patch.covered if self._wall_aware: leaf[[0, -1]] = False - for array, mask in zip(arrays, (leaf, self._fine_leaf)): + for array, mask in zip(arrays, (leaf, self._fine_leaf, self._child_leaf)): if not np.isfinite(array[:, mask]).all() or (array[:, mask].sum(axis=0) <= 0).any(): raise ValueError("Leaf populations must be finite with positive density") self._mask_nonleaves(*arrays) @@ -412,6 +474,16 @@ class FixedPatchStepper: 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 _advance_child(self, parent_buffer, **kw): + """Explode from the given L1 buffer, then two L2 hops (halo -> step).""" + pair, g = self._child, self._levels[2] + pair.explode(parent_buffer, g.current, pair.ghosts, pair.donors, np.uint32(pair.n_ghosts), + np.uint64(pair.parent_cells), + grid=((pair.n_ghosts + self._block[0] - 1) // self._block[0], 1), **kw) + for src, dst in ((g.current, g.spare), (g.spare, g.current)): + self._gather_halo(g, src, **kw) + g.kernel(src, dst, np.uint32(g.layout.halo), np.uintp(0), g.kind, grid=g.grid, **kw) + def _gather_halo(self, level, buffer, **kw): pairs = len(level.layout.halo_dst) if pairs: @@ -441,7 +513,17 @@ class FixedPatchStepper: # 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)): + if self._child is not None: + # L1 feeds L2 from its CURRENT G (fresh at the midpoint), + # L2 takes two hops, then is restricted into src. + self._advance_child(src, **kw) self._gather_halo(f, src, **kw) + if self._child is not None: + pair = self._child + pair.coalesce(self._levels[2].current, src, pair.targets, + pair.children, np.uint32(pair.n_rows), + grid=((pair.n_rows + self._block[0] - 1) + // self._block[0], 1), **kw) self._blocked_fine_boundary(src, **kw) f.kernel(src, dst, np.uint32(self.fine_layout.halo), flags, f.kind, grid=f.grid, **kw) @@ -477,7 +559,7 @@ class FixedPatchStepper: raise def download(self): - """Return two synchronized G arrays, NaN on all non-leaves; no host cache.""" + """Return synchronized G arrays (one per level), NaN on all non-leaves; no host cache.""" with self._lock, self.ctx.activate(): self._require("ready") arrays = [np.empty(level.shape if level.layout is None else level.layout.blocked_shape, @@ -508,7 +590,7 @@ class FixedPatchStepper: except Exception as exc: error = exc self._levels.clear() - self._explode = self._coalesce = None + self._explode = self._coalesce = self._child = None if getattr(self, "_fine_curved", None) is not None: try: self._fine_curved.free() diff --git a/tests/unit/test_refinement.py b/tests/unit/test_refinement.py index 4fcfe17..0f2df56 100644 --- a/tests/unit/test_refinement.py +++ b/tests/unit/test_refinement.py @@ -233,5 +233,32 @@ class FixedPatchSetTests(unittest.TestCase): self.build(self.FACE).fine_layout(block) +class NestPatchTests(unittest.TestCase): + """Three-level single chain: L2 bounds in L1 window cells (probe fixture).""" + def setUp(self): + self.cfg = config() + self.l1 = FixedPatch2D.from_config(self.cfg, (4, 4, 12, 12)) # 24x24 window + + def test_probe_fixture_composes_kinds(self): + child, kind = refinement.nest_patch(self.cfg, self.l1, (8, 8, 16, 16)) + self.assertEqual(child.fine_shape, (24, 24)) + self.assertAlmostEqual(child.tau_coarse, self.l1.tau_fine) + base = self.l1.fine_kind + for role in (INNER_GHOST, OUTER_GHOST): + np.testing.assert_array_equal(kind[base == role], role) + self.assertTrue((kind[child.covered] == OUTER_GHOST).all()) + self.assertEqual(int((kind == ACTIVE).sum()), 192) + self.assertEqual(int((kind == INNER_GHOST).sum()), 144) + self.assertEqual(int(child.covered.sum()), 64) + + def test_invalid_nesting_rejected(self): + for bounds in ((2, 8, 16, 16), # covers inherited INNER + (4, 8, 16, 16), # explosion donors in INNER + (8, 8, 20, 16), # donors reach OUTER/INNER at the far side + (5, 8, 16, 16)): # one-cell outward: restriction footprint/donors + with self.subTest(bounds=bounds), self.assertRaises(ValueError): + refinement.nest_patch(self.cfg, self.l1, bounds) + + if __name__ == "__main__": unittest.main(verbosity=2) diff --git a/tests/validation/multigrid/test_refinement_three_level.py b/tests/validation/multigrid/test_refinement_three_level.py new file mode 100644 index 0000000..8cb873d --- /dev/null +++ b/tests/validation/multigrid/test_refinement_three_level.py @@ -0,0 +1,89 @@ +"""Nested three-level FixedPatchStepper (child=..., single chain): CPU guards and GPU gates. + +Fixture = reviewed CPU probe: L0 16x16 periodic, L1 (4,4,12,12) -> 24x24 window, +L2 (8,8,16,16) in L1 window cells -> 24x24 window; 4x4 blocks on every level. +Run CPU and GPU classes in separate processes. GPU gates (frozen before any run): +- uniform moving equilibrium stays on all three levels' ACTIVE leaves, |dG| <= 1e-6; +- seeded non-equilibrium: every level's ACTIVE leaves finite, raw owner cores NaN on + OUTER/covered (INNER partly invalid by contract, not asserted); +- each L1 hop feeds L2 from the L1 buffer it is about to step (current, then the new + G written by the first hop), never from L0: 2 child advances per coarse step. +Numerical comparison against the CPU probe (G vs C(f)) is done by the acceptance run. +""" +import unittest +from unittest.mock import patch + +import numpy as np +import pycuda.driver as cuda + +from CelerisLab.config import LBMConfig +from CelerisLab.lbm.equilibrium import compute_feq_d2q9 +from CelerisLab.lbm.refinement import ACTIVE, OUTER_GHOST +from CelerisLab.lbm.refinement_stepper import FixedPatchStepper +from test_refinement_blocked_fine import state, circle +from test_refinement_blocked_base import raw + +CFG = LBMConfig(nx=16, ny=16, viscosity=.1) +L1, L2 = (4, 4, 12, 12), (8, 8, 16, 16) +KW = dict(fine_block=(4, 4), base_block=(4, 4), child=L2, child_block=(4, 4)) + + +class ThreeLevelCPU(unittest.TestCase): + def test_rejections_before_cuda(self): + context = object() # never activated: every case must fail first + bad = (dict(KW, child_block=None), dict(KW, fine_block=None), + dict(KW, fine_objects=[circle()]), dict(KW, child=(4, 8, 16, 16)), + dict(KW, child_block=(5, 4))) + for kw in bad: + with self.subTest(kw=kw), self.assertRaises(ValueError): + FixedPatchStepper(CFG, L1, context=context, **kw) + with self.assertRaises(ValueError): # several L1 patches + child: not yet + FixedPatchStepper(LBMConfig(nx=32, ny=16, viscosity=.1), ((4, 4, 12, 12), (20, 4, 28, 12)), + context=context, **KW) + + +class ThreeLevelGPU(unittest.TestCase): + def test_equilibrium_is_kept_on_every_level(self): + from CelerisLab.cuda.context import CudaContext + feq = np.asarray(compute_feq_d2q9(1.0, .03, -.02), np.float32)[:, None, None] + with CudaContext(0) as ctx, FixedPatchStepper(CFG, L1, context=ctx, **KW) as run: + shapes = ((16, 16), run.patch.fine_shape, run.child_patch.fine_shape) + run.initialize(*(np.broadcast_to(feq, (9, *s)).copy() for s in shapes)) + run.step(8) + self.assertEqual(run.step_count, 8) + for level, g in enumerate(run.download()): + leaf = np.isfinite(g[0]) + self.assertTrue(leaf.any(), f"level {level}") + np.testing.assert_allclose(g[:, leaf], np.broadcast_to(feq, g.shape)[:, leaf], + atol=1e-6, rtol=0, err_msg=f"level {level}") + + def test_nonequilibrium_masks_and_fresh_middle_feeds(self): + from CelerisLab.cuda.context import CudaContext + with CudaContext(0) as ctx, FixedPatchStepper(CFG, L1, context=ctx, **KW) as run: + p1, p2 = run.patch, run.child_patch + run.initialize(state(16, 16, 1), state(*p1.fine_shape, 2), state(*p2.fine_shape, 3)) + f = run._levels[1] + feeds, original = [], run._advance_child + def record(buffer, **kw): + feeds.append(int(buffer)) + original(buffer, **kw) + with patch.object(run, "_advance_child", side_effect=record): + cur, spare = int(f.current), int(f.spare) + run.step(1) + self.assertEqual(feeds, [cur, spare]) # current G, then the new G of hop 1 + run.step(7) + middle = run._middle_kind + for g, leaf in zip(run.download(), (~p1.covered, middle == ACTIVE, p2.fine_kind == ACTIVE)): + self.assertTrue(np.isfinite(g[:, leaf]).all()) + self.assertTrue(np.isnan(g[:, ~leaf]).all()) + rc, rf, rg = raw(run) + self.assertTrue(np.isnan(rc[:, p1.covered]).all()) + self.assertTrue(np.isfinite(rf[:, middle == ACTIVE]).all()) + self.assertTrue(np.isnan(rf[:, middle == OUTER_GHOST]).all()) # incl. child-covered + self.assertTrue(np.isnan(rf[:, p2.covered]).all()) + self.assertTrue(np.isfinite(rg[:, p2.fine_kind == ACTIVE]).all()) + self.assertTrue(np.isnan(rg[:, p2.fine_kind == OUTER_GHOST]).all()) + + +if __name__ == "__main__": + unittest.main()