Frank14f/CelerisLab
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Three levels: nested body-free child patch on a single chain

Browse Source 
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 <noreply@anthropic.com>
This commit is contained in:
Claude
2026-10-16 14:35:57 +00:00
co-authored by Claude Opus 5.5
parent 6301887093
commit 76555de4b8
4 changed files with 264 additions and 21 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/CelerisLab/lbm/refinement.py
+45
View File
@@ -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.
src/CelerisLab/lbm/refinement_stepper.py
+103 -21
View File
@@ -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()