FixedPatchStepper: level-indexed API, drop 15 fine/child compatibility views

- Hierarchy documents level/pair/shape/frame conventions; edge_frame(l)
  replaces the stepper's inline origin/spacing loop.
- Public per-level read-only state: layouts, boundaries, jets, leaf.
- set_jet_amplitude/read_raw_exchange/read_exchange_window: level= required.
- remap_child(child_objects=) -> remap_finest(objects=).
- _kernel_configs, _upload_boundary, _observed_boundary/_summarize split
  construction and observation duplication; behavior unchanged.
- fine_level_origin -> fine_level_offset (a (y, x) index, not an origin).
- docs: migration map in docs/development/src_coherence.md.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Claude
2026-10-28 01:39:48 +00:00
co-authored by Claude Opus 5.5
parent 0c64f821a4
commit b4eacf7526
9 changed files with 212 additions and 135 deletions
+25 -3
View File
@@ -21,7 +21,17 @@ from .refinement import (ACTIVE, OUTER_GHOST, FixedPatch2D, FixedPatchSet2D,
@dataclass(frozen=True, eq=False)
class Hierarchy:
"""configs[l]: validated LBMConfig of level l; pairs[l]: tables of (l, l+1);
kinds[l]: composed per-cell kind (child coverage OUTER on every non-finest level)."""
kinds[l]: composed per-cell kind (child coverage OUTER on every non-finest level).
Index conventions (the single place they are defined):
* level l = 0 (root) .. depth-1; pair index p = l-1 builds level l from level p;
* array shapes are (ny, nx); rectangles, origins and centers are (x, y);
* the root EDGE frame puts root cell i at [i, i+1), centers at i+.5 (bodies and
Jets of FixedPatchStepper); uniform Simulation objects use the integer-NODE frame;
* cell (i, j) of level l's allocated array (ghost band and holes included) covers
[o + h*i, o + h*(i+1)) with (o, h) = edge_frame(l); level 1 is built from the root
config directly, deeper levels from their parent frame (see edge_frame).
"""
configs: tuple
pairs: tuple
kinds: tuple
@@ -30,6 +40,18 @@ class Hierarchy:
def depth(self):
return len(self.kinds)
def edge_frame(self, level):
"""(origin (x, y), spacing) of level ``level``'s array in the root edge frame.
Each pair's fine_origin is in parent cells: zero for a set (level-wide grid),
the shifted window origin for a scalar FixedPatch2D.
"""
origin, spacing = (0.0, 0.0), 1.0
for pair in self.pairs[:level]:
origin = tuple(o + spacing * f for o, f in zip(origin, pair.fine_origin))
spacing /= 2
return origin, spacing
def is_rectangle(value):
value = tuple(value)
@@ -86,14 +108,14 @@ def _check_chain(coarse, fine, pair):
def leaf_partition(h):
"""Base-cell fraction published as ACTIVE over all levels (exact: powers of two).
Window (scalar) levels are offset by their fine_level_origin first.
Window (scalar) levels are offset by their fine_level_offset first.
"""
base = np.zeros(h.kinds[0].shape)
oy = ox = 0
for level, kind in enumerate(h.kinds):
if level:
pair = h.pairs[level - 1]
py, px = (pair.fine_level_origin if isinstance(pair, FixedPatch2D) else (0, 0))
py, px = (pair.fine_level_offset if isinstance(pair, FixedPatch2D) else (0, 0))
oy, ox = 2 * oy + py, 2 * ox + px
y, x = np.nonzero(kind == ACTIVE)
np.add.at(base, ((y + oy) >> level, (x + ox) >> level), 0.25 ** level)
+3 -2
View File
@@ -115,8 +115,9 @@ class FixedPatch2D:
return self.bounds[0] - HALO / 2, self.bounds[1] - HALO / 2
@property
def fine_level_origin(self) -> tuple[int, int]:
"""(y, x) of this window's first cell on the level-wide 2:1 fine grid."""
def fine_level_offset(self) -> tuple[int, int]:
"""Array index (y, x) of this window's first cell on the level-wide 2:1 fine grid
(an index offset, not an (x, y) origin)."""
return 2 * self.bounds[1] - HALO, 2 * self.bounds[0] - HALO
@property
+113 -121
View File
@@ -1,10 +1,11 @@
# CelerisLab/lbm/refinement_stepper.py
"""Opt-in fixed patch backend with optional stationary fine-interior circles.
"""Blocked fixed 2:1 refinement hierarchy (any depth) with optional stationary circles/Jets.
The caller owns the CudaContext and must close this backend before that context.
Only synchronized postcollision_g leaf populations are public. No checkpoint,
forcing or runtime parameter mutation is offered. Bounded manual periodic
L2 set replacement is available at completed base intervals via remap_child.
forcing or runtime parameter mutation besides Jet amplitudes is offered. A bounded
manual periodic level-2 set replacement is available at completed base intervals
via remap_finest. Level/pair/shape/frame conventions: lbm.hierarchy.Hierarchy.
"""
from dataclasses import dataclass, replace
from numbers import Integral
@@ -62,8 +63,22 @@ class _Level:
boundary: object = None
def _level_name(level):
return ("base", "fine", "child")[level] if level < 3 else f"level {level}"
def _kernel_configs(cfg, pairs, wall_aware):
"""Per-level kernel LBMConfig (window shape, viscosity doubling per level).
Geometry configs live in Hierarchy.configs; level 1 of a slot_reservoir root
replaces the root wall fields by the internal bounce-back/velocity ghost setup.
"""
configs = [replace(cfg)]
for index, pair in enumerate(pairs):
parent = configs[-1]
level = replace(parent, nx=pair.fine_shape[1], ny=pair.fine_shape[0],
viscosity=2 * parent.viscosity)
if index == 0 and wall_aware:
level = replace(level, y_wall_bc="bounce_back", top_wall_bc="velocity",
wall_slot_width=0, wall_normal_velocity=0.0, velocity=0.0)
configs.append(level)
return configs
MAX_GRID_DIM = 65535 # CUDA gridDim.y/z limit (rows, blocks)
@@ -163,11 +178,18 @@ class FixedPatchStepper:
lattice velocity); every whole circle and its BFL stencil must fit strictly
away from that level's transfer footprint.
``explosion="projected_linear"`` selects experimental transverse linear
reconstruction on each adjacent pair. Default "constant" is unchanged.
Projected mode requires periodic root, all-blocked SRT levels without
bodies/Jets or LES; this backend supplies no forcing. Geometry is generic,
but union/multilevel GPU and matched-physics qualification remain pending.
``explosion`` selects the parent->child reconstruction of every pair:
``"constant"`` (default) or ``"projected_linear"``, an optional transverse
linear reconstruction (Schornbaum-Ruede/Chen style) that is GPU-confirmed
more accurate than constant only on periodic two-level single-rectangle
Taylor-Green cases; union/multilevel GPU and matched-physics qualification
remain pending. It requires a periodic root and SRT levels without
bodies/Jets or LES; this backend supplies no forcing.
Per-level read-only state: ``hierarchy.pairs`` (pair l-1 builds level l),
``layouts``, ``boundaries`` (host geometry), ``jets`` and ``leaf`` masks,
each indexed by level. ``set_jet_amplitude``, ``read_raw_exchange`` and
``read_exchange_window`` take an explicit ``level=``.
One recursive hop advances every level: a child explodes from its parent's
current hop source, takes two hops, then the parent gathers halos, stages
@@ -233,42 +255,26 @@ class FixedPatchStepper:
layout = BlockLayout2D.build(*pairs[level - 1].fine_shape, *block,
periodic=(False, False))
if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM:
raise ValueError(f"{_level_name(level)} blocks exceed launch limits")
raise ValueError(f"level {level} blocks exceed launch limits")
self._layouts.append(layout)
if self._channel:
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).
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],
viscosity=2 * parent.viscosity))
# Per-level boundary geometry in the root edge frame (level l spacing 2**-l).
self._boundary_hosts = [None]
origin, spacing = (0.0, 0.0), 1.0
spacings = [spacing]
self._level_cfgs = _kernel_configs(cfg, pairs, self._wall_aware)
# Host boundary geometry per level in the root edge frame (Hierarchy.edge_frame).
self._boundaries = [None]
for level in range(1, depth):
pair = pairs[level - 1]
if level == 1:
if level == 1: # the root config carries the wall-aware boundary fields
host = prepare_fine_boundary(cfg, pair, level_objects[1],
wall_aware=self._wall_aware, jets=level_jets[1],
_shifted_blocked=self._ddf_shifting)
else:
origin, spacing = self.hierarchy.edge_frame(level - 1)
host = prepare_fine_boundary(self._level_cfgs[level - 1], pair, level_objects[level],
jets=level_jets[level], parent_origin=origin,
parent_spacing=spacing, _shifted_blocked=self._ddf_shifting)
self._boundary_hosts.append(host)
origin = tuple(o + spacing * f for o, f in zip(origin, pair.fine_origin))
spacing /= 2
spacings.append(spacing)
self._boundaries.append(host)
if not isinstance(observe_raw_exchange, bool):
raise ValueError("observe_raw_exchange must be bool")
self._observe_raw_exchange = observe_raw_exchange
@@ -280,7 +286,7 @@ class FixedPatchStepper:
self._window_valid = False
self._window_start = None
self._exchange_valid = False
self._exchange_spacing = tuple(spacings)
self._exchange_spacing = tuple(self.hierarchy.edge_frame(l)[1] for l in range(depth))
self._exchange_centers = tuple(tuple(tuple(o.geometry.center) for o in objects)
for objects in level_objects)
self._level_jets = level_jets
@@ -289,7 +295,7 @@ class FixedPatchStepper:
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],
_prepare_blocked_boundary(self._boundaries[level], self._layouts[level],
self._level_cfgs[level], level_jets[level],
_shifted_blocked=self._ddf_shifting)
for level in range(1, depth)]
@@ -297,13 +303,13 @@ class FixedPatchStepper:
if root is not None:
self._leaf[0] &= ~root.tables["solid"]
for level in range(1, depth):
host = self._boundary_hosts[level]
host = self._boundaries[level]
if host is not None:
self._leaf[level] &= (host.flags.reshape(pairs[level - 1].fine_shape) & SOLID) == 0
self._state = [mask.copy() for mask in self._leaf]
if self._channel:
self._leaf[0] &= channel_masks((cfg.ny, cfg.nx))[0]
# Private constructor inputs for the bounded manual L2 continuation.
# Private constructor inputs for remap_finest.
self._remap_cfg = replace(cfg)
self._level_blocks = level_blocks
self.ctx = context
@@ -325,21 +331,8 @@ class FixedPatchStepper:
for index, pair in enumerate(pairs):
self._pairs.append(self._build_pair(index, pair))
for level in self._levels:
b = level.boundary
if b is not None:
b.flags = self._upload(b.tables["flags"])
b.device = {k: self._upload(v) for k, v in b.tables.items()
if k in BLOCKED_BOUNDARY_DEVICE_COLUMNS}
if self._observe_raw_exchange:
b.device.update(rx=self._upload(b.tables["rx"]),
ry=self._upload(b.tables["ry"]))
count = len(b.tables["cell"])
b.last_raw = self._allocate(count * 4 * 4)
b.interval_integral = self._allocate(count * 4 * 8)
if self._observe_window:
b.window_integral = self._allocate(count * 4 * 8)
b.accumulate = level.module.get_function("BlockAccumulateExchange")
b.kernel = level.module.get_function("BlockCurvedBoundary")
if level.boundary is not None:
self._upload_boundary(level)
except BaseException:
try:
self.close()
@@ -347,6 +340,22 @@ class FixedPatchStepper:
pass # Preserve the original construction failure.
raise
def _upload_boundary(self, level):
"""Device columns, BFL kernel and (when observed) the exchange ledgers of one level."""
b = level.boundary
b.flags = self._upload(b.tables["flags"])
b.device = {k: self._upload(v) for k, v in b.tables.items()
if k in BLOCKED_BOUNDARY_DEVICE_COLUMNS}
if self._observe_raw_exchange:
b.device.update(rx=self._upload(b.tables["rx"]), ry=self._upload(b.tables["ry"]))
count = len(b.tables["cell"])
b.last_raw = self._allocate(count * 4 * 4)
b.interval_integral = self._allocate(count * 4 * 8)
if self._observe_window:
b.window_integral = self._allocate(count * 4 * 8)
b.accumulate = level.module.get_function("BlockAccumulateExchange")
b.kernel = level.module.get_function("BlockCurvedBoundary")
def _planned_bytes(self):
"""Both G buffers, halo pairs and boundary columns of every level (transfer tables excluded).
@@ -516,14 +525,15 @@ class FixedPatchStepper:
self._status = "failed"
raise
def remap_child(self, rectangles, *, child_objects=()):
"""Return a NEW ready stepper with an L2 set replaced at this base time.
def remap_finest(self, rectangles, *, objects=()):
"""Return a NEW ready stepper with the level-2 set replaced at this base time.
Periodic three-level set children only; L0/L1 and block sizes stay fixed.
A body run requires its original root-edge-frame ``child_objects``:
Periodic three-level hierarchies whose level 2 is a set only; levels 0/1 and
block sizes stay fixed. A body run requires its original root-edge-frame
level-2 ``objects``:
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)
Current level-2 Jets (including updated amplitudes)
are preserved. Local projected boundary objects are not valid inputs.
The lock and stream synchronization certify
current postcollision G and the completed 1:2:4 counts internally.
@@ -541,19 +551,20 @@ class FixedPatchStepper:
pairs = self.hierarchy.pairs
if (self._wall_aware or self._channel or self.hierarchy.depth != 3
or not isinstance(pairs[1], FixedPatchSet2D)
or self._boundary_hosts[1] is not None
or self._boundaries[1] is not None
or any(b is not None for b in self._level_boundaries[:2])):
raise ValueError("remap_child requires periodic three-level L2 sets with no L0/L1 bodies")
raise ValueError("remap_finest requires a periodic three-level hierarchy with a "
"level-2 set and no level-0/1 bodies")
cfg = self._remap_cfg
child_objects = tuple(child_objects)
objects = tuple(objects)
actual = self._level_boundaries[2]
jets = () if actual is None else actual.jets
old = RemapTopology.from_config(cfg, pairs[0], pairs[1].bounds,
child_objects=child_objects, child_jets=jets,
child_objects=objects, child_jets=jets,
_shifted_blocked=self._ddf_shifting)
if (not same_physical_boundary(old.boundary, self._boundary_hosts[2])
or (self._boundary_hosts[2] is not None and actual is None)):
raise ValueError("child_objects must reproduce the original root-edge-frame boundary")
if (not same_physical_boundary(old.boundary, self._boundaries[2])
or (self._boundaries[2] is not None and actual is None)):
raise ValueError("objects must reproduce the original root-edge-frame level-2 boundary")
if actual is not None:
velocity = actual.tables.get("surface_velocity")
if velocity is None:
@@ -561,18 +572,17 @@ class FixedPatchStepper:
_shifted_blocked=self._ddf_shifting)
if (not same_physical_boundary(old.boundary, actual.host)
or not np.array_equal(old.surface_velocity, velocity)):
raise ValueError("child_objects changed boundary ownership or surface velocity")
raise ValueError("objects changed boundary ownership or surface velocity")
new = RemapTopology.from_config(cfg, pairs[0], rectangles,
child_objects=child_objects, child_jets=jets,
child_objects=objects, child_jets=jets,
_shifted_blocked=self._ddf_shifting)
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 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)
layout = new.child.fine_layout(self._level_blocks[2])
if layout.nblocks > MAX_GRID_DIM or layout.local_shape[0] > MAX_GRID_DIM:
raise ValueError("child blocks exceed launch limits")
raise ValueError("level 2 blocks exceed launch limits")
with self.ctx.activate():
try:
self._stream.synchronize()
@@ -588,8 +598,8 @@ class FixedPatchStepper:
phase="postcollision_g", **options)
candidate = type(self)(
cfg, context=self.ctx, refinements=[pairs[0].bounds, new.child.bounds],
level_blocks=[base_block, fine_block, child_block],
level_objects=[(), (), child_objects], level_jets=[(), (), jets],
level_blocks=list(self._level_blocks),
level_objects=[(), (), objects], level_jets=[(), (), jets],
observe_raw_exchange=self._observe_raw_exchange,
observe_window_exchange=self._observe_window, explosion=self._explosion)
try:
@@ -606,10 +616,9 @@ class FixedPatchStepper:
raise
return candidate
def set_jet_amplitude(self, index, amplitude, *, level=1):
"""Replace one jet's amplitude between steps; no geometry/table rebuild.
def set_jet_amplitude(self, index, amplitude, *, level):
"""Replace Jet ``index`` of ``level_jets[level]`` between steps; no geometry/table rebuild.
``level=1`` preserves the L1 default; use ``level=2`` for child_jets.
Invalid input raises before any change. A device/sync error marks the
instance failed (as step/initialize do); host state commits only after
a successful upload.
@@ -618,7 +627,7 @@ class FixedPatchStepper:
self._require("new", "ready")
if (isinstance(level, bool) or not isinstance(level, Integral)
or not 1 <= level < self.hierarchy.depth):
raise ValueError("Jet level must name an existing refined level (1 = first)")
raise ValueError("Jet level must name an existing refined level (1 .. depth-1)")
boundary = self._level_boundaries[level]
jets = list(self._level_jets[level])
jets[index] = jets[index].with_amplitude(amplitude) # validates first
@@ -758,7 +767,23 @@ class FixedPatchStepper:
self._window_start = self._step_count
self._window_valid = True
def read_exchange_window(self, body_id, *, level=1):
def _observed_boundary(self, level, body_id):
if (isinstance(level, bool) or not isinstance(level, Integral)
or not 0 <= level < len(self._levels)):
raise ValueError("level must name an existing level")
b = self._level_boundaries[level]
if (b is None or isinstance(body_id, bool) or not isinstance(body_id, Integral)
or not 0 <= body_id < len(b.host.objects)):
raise ValueError("body_id must name a body on this level")
return b
def _summarize(self, b, rows, level, body_id, window, **kw):
return summarize_exchange(rows, b.tables["body_id"], b.tables["aperture"],
body_id=body_id, level=level, h=self._exchange_spacing[level],
time_end=self._step_count, window=window,
reference_center=self._exchange_centers[level][body_id], **kw)
def read_exchange_window(self, body_id, *, level):
"""Non-destructive read of [window start, current step count] in base-lattice units."""
with self._lock, self.ctx.activate():
self._require("ready")
@@ -769,26 +794,16 @@ class FixedPatchStepper:
duration = self._step_count - self._window_start
if duration < 1:
raise RuntimeError("exchange window requires at least one completed base interval")
if (isinstance(level, bool) or not isinstance(level, Integral)
or not 0 <= level < len(self._levels)):
raise ValueError("level must name an existing level")
b = self._level_boundaries[level]
if (b is None or isinstance(body_id, bool) or not isinstance(body_id, Integral)
or not 0 <= body_id < len(b.host.objects)):
raise ValueError("body_id must name a body on this level")
b = self._observed_boundary(level, body_id)
rows = np.empty((len(b.tables["cell"]), 4), np.float64)
try:
cuda.memcpy_dtoh(rows, b.window_integral)
except BaseException:
self._status = "failed"
raise
return summarize_exchange(rows, b.tables["body_id"], b.tables["aperture"],
body_id=body_id, level=level, h=self._exchange_spacing[level],
time_end=self._step_count, window="window",
reference_center=self._exchange_centers[level][body_id],
duration=duration)
return self._summarize(b, rows, level, body_id, "window", duration=duration)
def read_raw_exchange(self, body_id, *, level=1, window="last_base_step"):
def read_raw_exchange(self, body_id, *, level, window="last_base_step"):
"""Read raw boundary reaction about its geometric center in the root edge frame.
Repeated reads and step(0) preserve the window. last_substep samples
@@ -802,13 +817,7 @@ class FixedPatchStepper:
raise RuntimeError("raw exchange observation is disabled")
if not self._exchange_valid:
raise RuntimeError("raw exchange requires a completed base interval")
if (isinstance(level, bool) or not isinstance(level, Integral)
or not 0 <= level < len(self._levels)):
raise ValueError("level must name an existing level")
b = self._level_boundaries[level]
if (b is None or isinstance(body_id, bool) or not isinstance(body_id, Integral)
or not 0 <= body_id < len(b.host.objects)):
raise ValueError("body_id must name a body on this level")
b = self._observed_boundary(level, body_id)
if window not in ("last_base_step", "last_substep"):
raise ValueError("window must be last_base_step or last_substep")
count = len(b.tables["cell"])
@@ -818,10 +827,7 @@ class FixedPatchStepper:
except BaseException:
self._status = "failed"
raise
return summarize_exchange(rows, b.tables["body_id"], b.tables["aperture"],
body_id=body_id, level=level, h=self._exchange_spacing[level],
time_end=self._step_count, window=window,
reference_center=self._exchange_centers[level][body_id])
return self._summarize(b, rows, level, body_id, window)
def download(self):
"""Return physical G (FP64 for shifted storage), NaN on non-leaves."""
@@ -889,24 +895,10 @@ class FixedPatchStepper:
if exc[0] is None:
raise
# Read-only compatibility views over the level lists, kept on purpose for the existing
# tests and docs. Internal code must not use them; it reads hierarchy/_layouts/_level_*.
patch = property(lambda self: self.hierarchy.pairs[0])
child_patch = property(lambda self: self.hierarchy.pairs[1] if self.hierarchy.depth > 2 else None)
base_layout = property(lambda self: self._layouts[0])
fine_layout = property(lambda self: self._layouts[1])
child_layout = property(lambda self: self._layouts[2] if len(self._layouts) > 2 else None)
fine_boundary = property(lambda self: self._boundary_hosts[1])
child_boundary = property(lambda self: self._boundary_hosts[2] if len(self._boundary_hosts) > 2
else None)
_jets = property(lambda self: self._level_jets[1])
_boundaries = property(lambda self: {i: b for i, b in enumerate(self._level_boundaries) if i})
_fine_tables = property(lambda self: None if self._level_boundaries[1] is None
else self._level_boundaries[1].tables)
_actuated = property(lambda self: self._level_boundaries[1] is not None
and self._level_boundaries[1].actuated)
_jet_cache = property(lambda self: None if self._level_boundaries[1] is None
else self._level_boundaries[1].jet_cache)
_fine_leaf = property(lambda self: self._leaf[1])
_child_leaf = property(lambda self: self._leaf[2] if len(self._leaf) > 2 else None)
_middle_kind = property(lambda self: self.hierarchy.kinds[1] if self.hierarchy.depth > 2 else None)
# Read-only per-level state, indexed by level (0 = root). Pairs: hierarchy.pairs[l-1].
layouts = property(lambda self: tuple(self._layouts))
boundaries = property(lambda self: tuple(self._boundaries),
doc="Host boundary geometry per level (None: no bodies/Jets).")
jets = property(lambda self: tuple(self._level_jets))
leaf = property(lambda self: tuple(m.copy() for m in self._leaf),
doc="Per-level (ny, nx) masks of public (non-NaN) download cells.")