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

Build refinement transfers from union and pin physical body support

Browse Source
This commit is contained in:
Codex
2026-10-16 20:51:29 -07:00
parent 1e8e2c9a69
commit fbe73779c1
13 changed files with 357 additions and 106 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/CelerisLab/lbm/refinement.py
+104 -58
View File
@@ -151,20 +151,17 @@ class FixedPatch2D:
raise ValueError("levels are not synchronized: fine_steps must equal 2*coarse_steps")
MIN_PATCH_GAP = 4 # coarse cells: the full 4-fine-cell supports stay disjoint
MIN_PATCH_GAP = 4 # Legacy exported constant; union transfers impose no gap policy.
@dataclass(frozen=True, eq=False)
class FixedPatchSet2D:
"""Several separated FixedPatch2D rectangles sharing one 2:1 fine level.
"""Rectangles describing one 2:1 fine union, with unique transfer ownership.
The fine level is the whole (2*ny, 2*nx) grid in level coordinates; only
the patch supports (active region + HALO ghost band) are stored, by a
sparse BlockLayout2D from fine_layout(). Tables have FixedPatch2D meaning
and row/child order, re-addressed to level cells: ghost_fine are level
fine cells, coalesce_fine are q*N+cell level slots; coarse ids unchanged.
Patches are kept in canonical (y0, x0) order, so input order is irrelevant.
Active regions must be >= MIN_PATCH_GAP coarse cells apart (Chebyshev).
Input patches remain canonical metadata. Numerical kinds and tables are
built once from their union; adjacent, overlapping and concave unions are
accepted when their actual transfer dependencies close. Sparse storage owns
every whole tile intersecting the support, independent of rectangle edges.
"""
coarse_shape: tuple[int, int]
patches: tuple
@@ -182,42 +179,43 @@ class FixedPatchSet2D:
if not rectangles:
raise ValueError("at least one patch rectangle is required")
patches = sorted((FixedPatch2D.from_config(cfg, r) for r in rectangles),
key=lambda p: (p.bounds[1], p.bounds[0]))
for i, a in enumerate(patches):
for b in patches[i + 1:]:
(ax0, ay0, ax1, ay1), (bx0, by0, bx1, by1) = a.bounds, b.bounds
if max(ay0 - by1, by0 - ay1, ax0 - bx1, bx0 - ax1) < MIN_PATCH_GAP:
raise ValueError(f"patches {a.bounds} and {b.bounds} are closer "
f"than {MIN_PATCH_GAP} coarse cells")
key=lambda p: (p.bounds[1], p.bounds[0], p.bounds[3], p.bounds[2]))
ny, nx = patches[0].coarse_shape
fy, fx = 2 * ny, 2 * nx
if fy * fx * 9 > np.iinfo(np.uint64).max // 2:
raise ValueError("fine level exceeds 64-bit indexing")
kind = np.full((fy, fx), OUTER_GHOST, np.uint8)
ghost_fine, coalesce_fine = [], []
for p in patches:
y0, x0 = p.fine_level_origin
ly, lx = p.fine_shape
kind[y0:y0 + ly, x0:x0 + lx] = p.fine_kind
gy, gx = np.divmod(p.ghost_fine.astype(np.int64), lx)
ghost_fine.append((gy + y0) * fx + gx + x0)
q, cell = np.divmod(p.coalesce_fine.astype(np.int64), ly * lx)
cy, cx = np.divmod(cell, lx)
coalesce_fine.append(q * fy * fx + (cy + y0) * fx + cx + x0)
ghost_coarse = np.concatenate([p.ghost_coarse for p in patches])
coalesce_coarse = np.concatenate([p.coalesce_coarse for p in patches])
covered = np.zeros((ny, nx), bool)
for p in patches:
covered |= p.covered
cells = ny * nx
if (covered.ravel()[ghost_coarse].any()
or covered.ravel()[coalesce_coarse % np.uint64(cells)].any()
or len(np.unique(coalesce_coarse)) != len(coalesce_coarse)):
raise ValueError("patch transfers overlap another patch")
return cls((ny, nx), tuple(patches), patches[0].tau_coarse, patches[0].tau_fine,
_readonly(kind, np.uint8), _readonly(np.concatenate(ghost_fine), np.uint64),
_readonly(ghost_coarse, np.uint32), _readonly(coalesce_coarse, np.uint64),
_readonly(np.concatenate(coalesce_fine), np.uint64).reshape(-1, 4))
fy, fx = 2 * ny, 2 * nx
if fy * fx * 9 > np.iinfo(np.uint64).max // 2:
raise ValueError("fine level exceeds 64-bit indexing")
active = np.kron(covered, np.ones((2, 2), bool))
inner = _dilate(active, 2)
support = _dilate(active, HALO)
kind = np.full((fy, fx), OUTER_GHOST, np.uint8)
kind[inner] = INNER_GHOST
kind[active] = ACTIVE
gy, gx = np.nonzero(support & ~active)
ghost_fine = (gy * fx + gx).astype(np.uint64)
ghost_coarse = ((gy >> 1) * nx + (gx >> 1)).astype(np.uint32)
targets, children = [], []
ry, rx = np.nonzero(_dilate(covered, 1) & ~covered)
for y, x in zip(ry.tolist(), rx.tolist()):
for q in range(1, 9):
if _closed_union_contains(covered, 2 * x + 1 - D2Q9_EX[q], 2 * y + 1 - D2Q9_EY[q]):
targets.append(q * ny * nx + y * nx + x)
children.append([q * fy * fx + (2 * y + dy) * fx + 2 * x + dx
for dy in (0, 1) for dx in (0, 1)])
targets = np.asarray(targets, np.uint64)
children = np.asarray(children, np.uint64).reshape(-1, 4)
cells = np.asarray(children % np.uint64(fy * fx), np.intp)
if not (kind.ravel()[cells] == INNER_GHOST).all():
raise ValueError("coalesce children must be INNER ghost cells")
pair = cls((ny, nx), tuple(patches), patches[0].tau_coarse, patches[0].tau_fine,
_readonly(kind, np.uint8), _readonly(ghost_fine, np.uint64),
_readonly(ghost_coarse, np.uint32), _readonly(targets, np.uint64),
_readonly(children, np.uint64))
_transfer_footprint(pair)
_check_parent_pull_closure(pair)
return pair
@property
def fine_shape(self) -> tuple[int, int]:
@@ -240,22 +238,61 @@ class FixedPatchSet2D:
return mask
def fine_layout(self, block):
"""Sparse non-periodic layout storing exactly the patch supports."""
"""Sparse non-periodic tiles covering union support; domain divisibility remains."""
from .blocks import BlockLayout2D
block = tuple(block)
if len(block) != 2 or any(isinstance(v, (bool, np.bool_)) or not isinstance(v, Integral)
or v <= 0 for v in block):
raise ValueError("fine block must be two positive integers (by, bx)")
by, bx = map(int, block)
origins = set()
for p in self.patches:
y0, x0 = p.fine_level_origin
ly, lx = p.fine_shape
if y0 % by or x0 % bx or ly % by or lx % bx:
raise ValueError(f"patch {p.bounds} support is not aligned to fine blocks {block}")
origins.update((y, x) for y in range(y0, y0 + ly, by) for x in range(x0, x0 + lx, bx))
support = _dilate(self.fine_kind == ACTIVE, HALO)
ty, tx = np.nonzero(support)
origins = sorted({(int(y) // by * by, int(x) // bx * bx) for y, x in zip(ty, tx)})
return BlockLayout2D.build(*self.fine_shape, by, bx, periodic=(False, False),
origins=sorted(origins))
origins=origins)
def _dilate(mask, steps):
"""Chebyshev dilation with a zero (non-periodic) exterior."""
out = np.asarray(mask, bool)
for _ in range(steps):
grown = out.copy()
grown[1:] |= out[:-1]
grown[:-1] |= out[1:]
wide = grown.copy()
wide[:, 1:] |= grown[:, :-1]
wide[:, :-1] |= grown[:, 1:]
out = wide
return out
def _closed_union_contains(covered, mx, my):
"""Does doubled point (mx, my) lie in a closed covered square [2X, 2X+2] x [2Y, 2Y+2]?"""
xs = (mx // 2 - 1, mx // 2) if mx % 2 == 0 else ((mx - 1) // 2,)
ys = (my // 2 - 1, my // 2) if my % 2 == 0 else ((my - 1) // 2,)
ny, nx = covered.shape
return any(0 <= y < ny and 0 <= x < nx and covered[y, x] for y in ys for x in xs)
def _check_parent_pull_closure(pair):
"""Every uncovered parent cell pulling q from a covered cell needs its (R, q) row; no row
may stage into an uncovered parent cell unless its link passes the closed union."""
ny, nx = pair.coarse_shape
rows = set(map(int, pair.coalesce_coarse))
if len(rows) != len(pair.coalesce_coarse):
raise ValueError("duplicate coalesce writer (receiver, q)")
covered = pair.covered
for y, x in zip(*np.nonzero(~covered)):
for q in range(1, 9):
sy, sx = y - D2Q9_EY[q], x - D2Q9_EX[q]
if 0 <= sy < ny and 0 <= sx < nx and covered[sy, sx]:
if q * ny * nx + y * nx + x not in rows:
raise ValueError("a surviving parent ACTIVE pull lacks its coalesce row")
for slot in rows:
q, cell = divmod(slot, ny * nx)
y, x = divmod(cell, nx)
if covered[y, x] or not _closed_union_contains(covered, 2*x + 1 - D2Q9_EX[q], 2*y + 1 - D2Q9_EY[q]):
raise ValueError("coalesce row outside the union boundary")
@dataclass(frozen=True, eq=False)
@@ -323,14 +360,14 @@ def nest_patch(cfg, parent, child_bounds):
def nest_patch_set(cfg, parent, child_rectangles_in_parent):
"""Separated child union in the parent's allocated L1 array coordinates.
"""Child rectangle union in the parent's allocated L1 array coordinates.
The root edge coordinate of L1 index (x, y) is
``parent.fine_origin + (x, y)/2``. For a set parent its origin is zero;
a scalar parent's window includes its ghost band and has a shifted origin.
The child set allocates a dense public extent twice that L1 array, while
its blocked storage omits holes. Both adjacent pairs retain gap >= 4 in
their respective parent-cell units. Validate support over the whole union.
its blocked storage omits holes. Validate the actual adjacent-pair support over the whole union;
no rectangle-gap policy is imposed.
"""
child = FixedPatchSet2D.from_config(
_middle_config(cfg, parent), child_rectangles_in_parent)
@@ -373,6 +410,17 @@ def _nested_kind(parent, child):
return _readonly(composed, np.uint8)
def _circle_in_coverage(covered, cx, cy, radius):
"""Exact disk versus uncovered parent squares; no containing-rectangle rule."""
ny, nx = covered.shape
if cx-radius < 0 or cy-radius < 0 or cx+radius > nx or cy+radius > ny:
return False
y, x = np.nonzero(~covered)
dx = np.maximum(np.maximum(x-cx, cx-(x+1)), 0.)
dy = np.maximum(np.maximum(y-cy, cy-(y+1)), 0.)
return not np.any(dx*dx + dy*dy < radius*radius)
def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=(),
parent_origin=(0., 0.), parent_spacing=1.):
"""Project fixed circles and reject boundary/transfer overlap before CUDA.
@@ -386,7 +434,7 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=(),
Original objects are not retained or mutated. No runtime geometry updates.
``patch`` may be a FixedPatchSet2D: objects are then projected once onto
the level-wide fine grid (fine_origin (0, 0)); each circle must lie inside
one patch, body ids are the positions in ``objects`` (never patch ids), and
the union, body ids are the positions in ``objects`` (never patch ids), and
jets and the (receiver, dir) writer checks run once over the whole level.
``jets`` use the same coarse cell-edge frame as the objects; their
selection is checked against true fine hits here. Jet amplitudes need no
@@ -416,7 +464,6 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=(),
flags = np.full(nx * ny, FLUID, np.uint16)
active = patch.fine_kind.ravel() == ACTIVE
local, links = [], []
rects = patch.bounds if isinstance(patch, FixedPatchSet2D) else (patch.bounds,)
covered = patch.covered
dx = dt = .5 * parent_spacing
fine_origin = origin + parent_spacing * np.asarray(patch.fine_origin)
@@ -430,9 +477,8 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=(),
px, py = (np.asarray((cx, cy)) - origin) / parent_spacing
pr = radius / parent_spacing
if (not np.isfinite([cx, cy, radius]).all() or radius <= 0
or not any(x0 <= px-pr and px+pr <= x1 and y0 <= py-pr
and py+pr <= y1 for x0, y0, x1, y1 in rects)):
raise ValueError("Entire fine circle must lie inside one physical patch")
or not _circle_in_coverage(covered, px, py, pr)):
raise ValueError("Entire fine circle must lie inside the physical patch union")
center = tuple((v-o)/dx-.5 for v, o in zip((cx, cy), fine_origin))
spec = obj.surface_actuation
if spec is not None:
src/CelerisLab/lbm/refinement_stepper.py
+7 -8
View File
@@ -170,13 +170,13 @@ class FixedPatchStepper:
re-uploads the per-link velocity column only.
Several patches: pass ``bounds`` as a sequence of (x0, y0, x1, y1)
rectangles (FixedPatchSet2D, >= 4 coarse cells apart). This requires
rectangles (FixedPatchSet2D, dependency-validated union). This requires
fine_block and base_block; fine_objects/jets are projected once onto the
level-wide fine grid (each circle inside one patch, body ids = object
level-wide fine grid (each circle inside the union, body ids = object
positions, jet and writer checks level-wide), and
uses one level-wide fine grid (2*ny, 2*nx): public fine arrays have that
shape, NaN outside the patches' ACTIVE cells. The step sequence is the
blocked one above, with every patch's transfers in one launch each.
blocked one above, with unique union transfers in one launch each.
Three-level separated unions: ``child_rectangles_in_parent=[...]`` is
mutually exclusive with scalar ``child`` and accepts either parent form.
@@ -507,8 +507,8 @@ class FixedPatchStepper:
Periodic three-level set children only; L0/L1 and block sizes stay fixed.
A body run requires its original root-edge-frame ``child_objects``:
stationary circles pinned in unchanged complete L2 rectangles.
Separated fluid-only rectangles may be added, moved or removed.
stationary circles with identical physical boundary and BFL support.
Other fluid coverage may change, even in the same connected union.
Current L2 Jets (including updated amplitudes)
are preserved. Local projected boundary objects are not valid inputs.
The lock and stream synchronization certify
@@ -547,10 +547,9 @@ class FixedPatchStepper:
raise ValueError("child_objects changed boundary ownership or surface velocity")
new = RemapTopology.from_config(cfg, pairs[0], rectangles,
child_objects=child_objects, child_jets=jets)
if (old.pinned_rectangles != new.pinned_rectangles
or not same_physical_boundary(old.boundary, new.boundary)
if (not same_physical_boundary(old.boundary, new.boundary)
or not np.array_equal(old.surface_velocity, new.surface_velocity)):
raise ValueError("remap requires the identical pinned body rectangles and boundary")
raise ValueError("remap requires the identical pinned physical boundary and BFL support")
base_block, fine_block, child_block = self._level_blocks
# Check sparse block alignment/launch size before any device read/allocation.
layout = new.child.fine_layout(child_block)
src/CelerisLab/lbm/synchronized_remap.py
+12 -10
View File
@@ -1,7 +1,7 @@
"""CPU-only whole-volume G transfer, not a runtime adaptive solver.
Periodic three-level grids with fixed L0/L1 topology are supported, including
stationary circles pinned inside unchanged L2 rectangles. The caller
stationary circles with unchanged physical boundary/BFL support. The caller
certifies that current owned leaves are quiescent,
unshifted D2Q9 postcollision G at the supplied completed base interval.
No device buffers, device boundary state, ghosts, collision or time advance
@@ -20,13 +20,14 @@ from ..body.jets import surface_link_velocity
@dataclass(frozen=True, eq=False)
class RemapTopology:
"""Validated fixed parent and separated L2 set; use ``from_config``.
"""Validated fixed parent and L2 union; use ``from_config``.
Child rectangles are in the parent's allocated L1 coordinates, including
its window offset/halo when the parent is a scalar patch. Dense masks are
rebuilt from geometry, never inferred from finite values in input arrays.
``pinned_rectangles`` canonically lists each body-containing rectangle once;
fluid-only layouts have no pins.
``pinned_rectangles`` is legacy metadata for single-rectangle containment.
Eligibility pins physical boundary/stencil data, not entire rectangles or
connected fluid regions; each topology independently validates transfer closure.
"""
parent: object
child: FixedPatchSet2D
@@ -59,9 +60,11 @@ class RemapTopology:
for obj in child_objects:
cx, cy = (np.asarray(obj.geometry.center) - parent.fine_origin) * 2
radius = obj.geometry._radius * 2
containing.add(next(r for r in child.bounds if
r[0] <= cx-radius and cx+radius <= r[2] and
r[1] <= cy-radius and cy+radius <= r[3]))
match = next((r for r in child.bounds if
r[0] <= cx-radius and cx+radius <= r[2] and
r[1] <= cy-radius and cy+radius <= r[3]), None)
if match is not None:
containing.add(match)
pinned = tuple(sorted(containing))
velocity = boundary_velocity(cfg, child, boundary, child_jets)
if velocity is not None:
@@ -167,10 +170,9 @@ def remap_synchronized(old, new, old_g, *, base_count, level_steps, phase,
or not _same_parent(old.parent, new.parent)
or old.child.tau_fine != new.child.tau_fine):
raise ValueError("remap requires unchanged L0/L1 geometry and relaxation")
if (old.pinned_rectangles != new.pinned_rectangles
or not same_physical_boundary(old.boundary, new.boundary)
if (not same_physical_boundary(old.boundary, new.boundary)
or not np.array_equal(old.surface_velocity, new.surface_velocity)):
raise ValueError("remap requires the identical pinned body rectangles and physical boundary")
raise ValueError("remap requires the identical pinned physical boundary and BFL support")
if not isinstance(level_steps, (tuple, list)):
raise ValueError("level_steps must be three completed integer counts")
counts = (base_count, *level_steps)