From 20fb16ded5d7d5e0550eb187f62c7371b52ae35c Mon Sep 17 00:00:00 2001 From: TRON Date: Fri, 16 Oct 2026 08:49:08 -0700 Subject: [PATCH] Support stationary bodies and geometry Jets on nested L2 --- README.md | 32 +++- src/CelerisLab/lbm/refinement.py | 27 ++- src/CelerisLab/lbm/refinement_stepper.py | 166 +++++++++++------ tests/unit/test_refinement.py | 69 +++++++ .../multigrid/test_refinement_three_level.py | 173 +++++++++++++++++- 5 files changed, 400 insertions(+), 67 deletions(-) diff --git a/README.md b/README.md index 3e9caf8..ee0b029 100644 --- a/README.md +++ b/README.md @@ -158,30 +158,48 @@ Scope and limits: coarse_g, fine_g = run.download() # coarse: NaN on every covered cell; # fine: NaN on every non-ACTIVE cell ``` -- Three levels (one nested chain, body-free): add `child=(x0, y0, x1, y1)` in +- Three levels (one nested chain): add `child=(x0, y0, x1, y1)` in **L1 window cells** (the patch's fine grid including its 4-cell ghost band) and `child_block`, with `fine_block` and `base_block`. L2 must sit inside L1's ACTIVE cells, away from L1's ghost band and transfer footprint (checked before any GPU allocation). `initialize` and `download` take/return three arrays; each level is NaN on its non-ACTIVE (and child-covered) cells. One coarse step = 2 L1 - steps = 4 L2 steps. One L1 patch only, no bodies/Jets. + steps = 4 L2 steps. One L1 patch only. Optional stationary circles and Jets + live on L2 via `child_objects` and `child_jets`; `fine_objects`/`jets` retain + their L1 meaning and must be empty in this mode. All object and Jet geometry + uses the **L0 edge frame**, not the child bounds' L1 window frame. L2 has + dx=dt=0.25, so prescribed lattice velocity is unchanged. Each whole body and + its BFL stencil must be inside L2 ACTIVE cells, outside the transfer footprint. + A Jet may select one physical body across multiple blocks, never two bodies. ```python + from CelerisLab.body.geometry.circle import CircleGeometry + from CelerisLab.body.objects import SimObject + from CelerisLab.body.jets import Jet, Tube, Parallel + cfg = LBMConfig(nx=16, ny=16, viscosity=.1) + body = SimObject(0, geometry=CircleGeometry(8., 8., .8), center=(8., 8.), radius=.8) + jet = Jet(Tube((6.5, 8.), 0., 3., .3), Parallel(0.), amplitude=.002) rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None] with CudaContext(0) as ctx, FixedPatchStepper( cfg, (4, 4, 12, 12), context=ctx, fine_block=(4, 4), base_block=(4, 4), - child=(8, 8, 16, 16), child_block=(4, 4)) as run: + child=(8, 8, 16, 16), child_block=(4, 4), + child_objects=[body], child_jets=[jet]) as run: shapes = (run.patch.coarse_shape, run.patch.fine_shape, run.child_patch.fine_shape) run.initialize(*(np.broadcast_to(rest, (9, *s)) for s in shapes)) # 16², 24², 24² - run.step(32) - coarse_g, l1_g, l2_g = run.download() + run.step(4) + run.set_jet_amplitude(0, -.002, level=2) # default level=1 still addresses L1 + run.step(4) + coarse_g, l1_g, l2_g = run.download() # solids/covered/nonleaves are NaN ``` - Not available: per-link surface exchange/force readout in blocked mode, - checkpoints, moving bodies, bodies/Jets with three levels, + checkpoints, moving bodies, L1 bodies/Jets in a three-level chain, several patches per level with nesting, dynamic refinement, `Simulation` integration. Validation so far is storage/operator equivalence and short runs, - not accuracy, convergence or long-time stability. + not accuracy, convergence or long-time stability. The new L2-body GPU gates + (layouts, masks, amplitude reversal, legacy equivalence) and independent + every-substep numerical replay remain required acceptance checks; CPU tests + alone do not establish GPU numerical correctness. ## Installation diff --git a/src/CelerisLab/lbm/refinement.py b/src/CelerisLab/lbm/refinement.py index 5a769a0..288fffa 100644 --- a/src/CelerisLab/lbm/refinement.py +++ b/src/CelerisLab/lbm/refinement.py @@ -349,12 +349,16 @@ def nest_patch(cfg, parent, child_bounds): return child, _readonly(composed, np.uint8) -def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=()): +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. - Input centers/radii are in coarse physical coordinates (dx=dt=1). + Input centers/radii and Jets are in the root edge frame. ``parent_origin`` + and ``parent_spacing`` map this patch's parent grid into that frame; the + default preserves the L0->L1 convention. cfg must describe the parent grid. SurfaceActuation retains its lattice-velocity meaning: input coarse units - are converted into independent fine-local objects with dx=dt=.5. + are converted into independent fine-local objects with + dx=dt=parent_spacing/2. 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 @@ -364,6 +368,10 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=()): selection is checked against true fine hits here. Jet amplitudes need no rescaling (lattice velocity scales by dt/dx = 1 under acoustic 2:1). """ + origin = np.asarray(parent_origin, dtype=np.float64) + if (origin.shape != (2,) or not np.isfinite(origin).all() + or not np.isfinite(parent_spacing) or parent_spacing <= 0): + raise ValueError("Parent frame must have a finite origin and positive spacing") if not objects: if jets: raise ValueError("Fine jets need fine boundary objects") @@ -386,7 +394,8 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=()): local, links = [], [] rects = patch.bounds if isinstance(patch, FixedPatchSet2D) else (patch.bounds,) covered = patch.covered - dx = dt = .5 + dx = dt = .5 * parent_spacing + fine_origin = origin + parent_spacing * np.asarray(patch.fine_origin) for bid, obj in enumerate(objects): if (type(obj.geometry) is not CircleGeometry or obj.is_sensor or obj.is_force_region or any(getattr(obj.state, name) != 0 for name in ("vx", "vy", "omega")) @@ -394,11 +403,13 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=()): raise ValueError("Fine boundaries require stationary circles without other roles") cx, cy = obj.geometry.center radius = obj.geometry._radius + 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 <= cx-radius and cx+radius <= x1 and y0 <= cy-radius - and cy+radius <= y1 for x0, y0, x1, y1 in rects)): + 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") - center = tuple((v-o)/dx-.5 for v, o in zip((cx, cy), patch.fine_origin)) + center = tuple((v-o)/dx-.5 for v, o in zip((cx, cy), fine_origin)) spec = obj.surface_actuation if spec is not None: spec = replace(spec, normal_velocity=spec.normal_velocity * dt / dx) @@ -465,7 +476,7 @@ def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False, jets=()): raise ValueError("Transfer footprint touches the coarse wall reconstruction stencil") from ..body.jets import fine_hits_to_edge_frame, surface_link_velocity hits = fine_hits_to_edge_frame([(link.hit_x, link.hit_y) for link in links], - patch.fine_origin, dx) + fine_origin, dx) packed = pack_cut_links_to_soa(links) if jets: # Same assignment checks as the stepper: empty/overlapping jets and diff --git a/src/CelerisLab/lbm/refinement_stepper.py b/src/CelerisLab/lbm/refinement_stepper.py index 4385bac..0fd2fcd 100644 --- a/src/CelerisLab/lbm/refinement_stepper.py +++ b/src/CelerisLab/lbm/refinement_stepper.py @@ -93,6 +93,36 @@ class _Level: halo_copy: object = None halo_dst: object = None halo_src: object = None + boundary: object = None + + +@dataclass +class _BlockedBoundary: + """One level's stationary BFL data and mutable Jet velocity column.""" + host: object + tables: dict + jets: tuple + actuated: bool + jet_cache: object = None + flags: object = None + device: object = None + kernel: object = None + + +def _prepare_blocked_boundary(boundary, layout, cfg, jets): + if boundary is None or layout is None: + return None + fluid_idx, dirs, q, body_id, rx, ry, _, fallback = boundary.curved_host + tables, _ = blocked_link_tables(layout, boundary.flags, fluid_idx, dirs, q, fallback) + actuated = bool(jets) or any(o.surface_actuation is not None for o in boundary.objects) + if actuated: + cfg = replace(cfg) + cfg.validate() + tables["surface_velocity"], _ = surface_link_velocity( + boundary.objects, cfg, boundary.flags, fluid_idx, dirs, body_id, rx, ry, + jets=jets, hits=boundary.curved_hits) + return _BlockedBoundary(boundary, tables, jets, actuated, + jet_link_shapes(jets, boundary.curved_hits) if jets else None) class FixedPatchStepper: @@ -140,9 +170,18 @@ class FixedPatchStepper: 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. + + Three-level single chain: ``child`` bounds use L1 window cells; all three + levels must be blocked. ``child_objects`` and ``child_jets`` use the L0 + edge frame, with L2 dx=dt=.25 and unchanged lattice velocity. L1 bodies + and Jets remain unsupported in a chain. Every whole circle and its BFL + stencil must fit strictly away from the L2 transfer footprint. L2 runs + halo -> BFL -> step four times per base step. Runtime child Jet changes + require ``set_jet_amplitude(index, amplitude, level=2)``. """ def __init__(self, cfg, bounds, *, context, coarse_boundary="periodic", fine_objects=(), - fine_block=None, base_block=None, jets=(), child=None, child_block=None): + fine_block=None, base_block=None, jets=(), child=None, child_block=None, + child_objects=(), child_jets=()): if HALO != 4: raise RuntimeError("refinement.cu requires a four-cell fine halo") if coarse_boundary not in ("periodic", "slot_reservoir"): @@ -164,7 +203,10 @@ class FixedPatchStepper: rects = self.patch.bounds if multipatch else (self.patch.bounds,) if any(r[1] < 4 or r[3] > cfg.ny - 4 for r in rects): raise ValueError("slot_reservoir requires y0 >= 4 and y1 <= ny-4") - fine_objects = tuple(fine_objects) + fine_objects, jets = tuple(fine_objects), tuple(jets) + child_objects, child_jets = tuple(child_objects), tuple(child_jets) + if child is None and (child_objects or child_jets): + raise ValueError("child_objects/child_jets require a nested child") if base_block is not None and fine_block is None: raise ValueError("base_block requires fine_block") self.base_layout = None @@ -193,7 +235,7 @@ class FixedPatchStepper: 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") + raise ValueError("a nested child supports bodies/Jets only on L2; use child_objects/child_jets") 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) @@ -205,15 +247,33 @@ class FixedPatchStepper: self._jets = tuple(jets) self.fine_boundary = prepare_fine_boundary( cfg, self.patch, fine_objects, wall_aware=self._wall_aware, jets=self._jets) - self._fine_tables, self._actuated, self._jet_cache = None, False, None - if self.fine_layout is not None and self.fine_boundary is not None: - 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) + fine_cfg = replace(cfg, nx=self.patch.fine_shape[1], ny=self.patch.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.child_boundary = (None if self.child_patch is None else prepare_fine_boundary( + fine_cfg, self.child_patch, child_objects, jets=child_jets, + parent_origin=self.patch.fine_origin, parent_spacing=.5)) + child_cfg = (None if self.child_patch is None else replace( + fine_cfg, nx=self.child_patch.fine_shape[1], ny=self.child_patch.fine_shape[0], + viscosity=4 * cfg.viscosity)) + self._boundaries = { + 1: _prepare_blocked_boundary(self.fine_boundary, self.fine_layout, fine_cfg, self._jets), + 2: _prepare_blocked_boundary(self.child_boundary, self.child_layout, child_cfg, child_jets), + } + # Compatibility aliases for the existing L1 inspection API. + b = self._boundaries[1] + self._fine_tables = None if b is None else b.tables + self._actuated = False if b is None else b.actuated + self._jet_cache = None if b is None else b.jet_cache 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) + if self.child_boundary is not None: + self._child_leaf &= (self.child_boundary.flags.reshape( + self.child_patch.fine_shape) & SOLID) == 0 self._child = None if self.fine_boundary is not None: self._fine_leaf &= (self.fine_boundary.flags.reshape(self.patch.fine_shape) & SOLID) == 0 @@ -228,31 +288,13 @@ class FixedPatchStepper: self._coalesce = None self._link_tables = {} self._block = (cfg.threads_per_block, 1, 1) - fine_cfg = replace(cfg, nx=self.patch.fine_shape[1], ny=self.patch.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) - if self._fine_tables is not None: - # Any actuation spec (even zero amplitude) requires the velocity buffer. - b = self.fine_boundary - self._actuated = bool(self._jets) or any( - obj.surface_actuation is not None for obj in b.objects) - if self._actuated: - fine_cfg.validate() # derives nq etc.; the unique-writer key needs them - fluid_idx, dirs, _, body_id, rx, ry, _, _ = b.curved_host - self._fine_tables["surface_velocity"], _ = surface_link_velocity( - b.objects, fine_cfg, b.flags, fluid_idx, dirs, body_id, rx, ry, - jets=self._jets, hits=b.curved_hits) - if self._jets: - self._jet_cache = jet_link_shapes(self._jets, b.curved_hits) with self.ctx.activate(): try: self._stream = cuda.Stream() 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, + specs.append((child_cfg, self.child_layout, self.child_patch.tau_fine, "PatchFineStep")) for index, (config, layout, tau, name) in enumerate(specs): source = "refinement.cu" @@ -284,6 +326,7 @@ class FixedPatchStepper: 1 if layout is None else layout.nblocks) level = _Level(shape, module, self._allocate(size), self._allocate(size), module.get_function(name), grid, layout) + level.boundary = self._boundaries.get(index) self._levels.append(level) if layout is not None: level.halo_copy = module.get_function("BlockHaloCopy") @@ -330,11 +373,15 @@ class FixedPatchStepper: 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() - if k != "flags"} - self._curved = self._levels[1].module.get_function("BlockCurvedBoundary") + 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 != "flags"} + b.kernel = level.module.get_function("BlockCurvedBoundary") + if self._boundaries[1] is not None: + b = self._boundaries[1] + self._fine_flags, self._link_tables, self._curved = b.flags, b.device, b.kernel elif self.fine_boundary is not None: from .curved_links import CurvedLinkSoA boundary = self.fine_boundary @@ -431,58 +478,72 @@ class FixedPatchStepper: cl.link_exchange_gpu if cl.link_exchange_gpu is not None else np.uintp(0), grid=((cl.count + self._block[0] - 1) // self._block[0], 1), **kw) - def set_jet_amplitude(self, index, amplitude): + def set_jet_amplitude(self, index, amplitude, *, level=1): """Replace one jet's amplitude 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. """ with self._lock, self.ctx.activate(): self._require("new", "ready") - jets = list(self._jets) + if (isinstance(level, bool) or not isinstance(level, Integral) + or level not in (1, 2) or level >= len(self._levels)): + raise ValueError("Jet level must be 1 or an existing child level 2") + boundary = self._boundaries[level] + jets = list(self._jets if level == 1 else (() if boundary is None else boundary.jets)) jets[index] = jets[index].with_amplitude(amplitude) # validates first velocity = None - if self._fine_curved is None and self._jet_cache is not None: - shape, jet_id = self._jet_cache - velocity = self._fine_tables["surface_velocity"].copy() + if boundary is not None and boundary.jet_cache is not None: + shape, jet_id = boundary.jet_cache + velocity = boundary.tables["surface_velocity"].copy() driven = jet_id >= 0 velocity[driven] = jet_link_velocity(jets, shape, jet_id)[driven] try: - if self._fine_curved is not None: + if level == 1 and self._fine_curved is not None: self._fine_curved.set_jet_amplitude(index, amplitude, self._stream) elif velocity is not None: self._stream.synchronize() # no in-flight step reads the column - cuda.memcpy_htod(self._link_tables["surface_velocity"], velocity) + cuda.memcpy_htod(boundary.device["surface_velocity"], velocity) except BaseException: self._status = "failed" raise if velocity is not None: - self._fine_tables["surface_velocity"] = velocity - self._jets = tuple(jets) + boundary.tables["surface_velocity"] = velocity + if boundary is not None: + boundary.jets = tuple(jets) + if level == 1: + self._jets = tuple(jets) def _blocked_fine_boundary(self, buffer, **kw): - if self._fine_tables is None: + self._blocked_boundary(self._levels[1], buffer, **kw) + + def _blocked_boundary(self, level, buffer, **kw): + b = level.boundary + if b is None: return - if self._curved is None or self._fine_flags is None: - raise RuntimeError("Fine boundary device buffers are missing") - t, count = self._link_tables, len(self._fine_tables["cell"]) + if b.kernel is None or b.flags is None: + raise RuntimeError("Boundary device buffers are missing") + t, count = b.device, len(b.tables["cell"]) velocity = t.get("surface_velocity") - if self._actuated and velocity is None: + if b.actuated and velocity is None: raise RuntimeError("Surface velocity device buffer is missing") - self._curved(buffer, t["block"], t["cell"], t["dir"], t["q"], t["fallback"], - np.uint32(count), velocity if velocity is not None else np.uintp(0), - grid=((count + self._block[0] - 1) // self._block[0], 1), **kw) + b.kernel(buffer, t["block"], t["cell"], t["dir"], t["q"], t["fallback"], + 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).""" + """Explode from the given L1 buffer, then two L2 hops (halo -> BFL -> 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) + self._blocked_boundary(g, src, **kw) + flags = np.uintp(0) if g.boundary is None else g.boundary.flags + g.kernel(src, dst, np.uint32(g.layout.halo), flags, g.kind, grid=g.grid, **kw) def _gather_halo(self, level, buffer, **kw): pairs = len(level.layout.halo_dst) @@ -589,6 +650,9 @@ class FixedPatchStepper: self._stream.synchronize() except Exception as exc: error = exc + for b in getattr(self, "_boundaries", {}).values(): + if b is not None: + b.flags = b.device = b.kernel = None self._levels.clear() self._explode = self._coalesce = self._child = None if getattr(self, "_fine_curved", None) is not None: diff --git a/tests/unit/test_refinement.py b/tests/unit/test_refinement.py index a993e35..5dd80c2 100644 --- a/tests/unit/test_refinement.py +++ b/tests/unit/test_refinement.py @@ -310,5 +310,74 @@ class NestPatchTests(unittest.TestCase): refinement.nest_patch(self.cfg, self.l1, bounds) +class ChildBoundaryTests(unittest.TestCase): + @staticmethod + def body(center=(8., 8.), radius=.8, obj_id=99, spec=None): + from CelerisLab.body.geometry.circle import CircleGeometry + from CelerisLab.body.objects import SimObject + return SimObject(obj_id, geometry=CircleGeometry(*center, radius), center=center, + radius=radius, surface_actuation=spec) + + def prepare(self, objects, jets=(), large=False): + from dataclasses import replace + cfg = LBMConfig(nx=32 if large else 16, ny=32 if large else 16, viscosity=.1) + parent = FixedPatch2D.from_config(cfg, (4, 4, 28, 28) if large else (4, 4, 12, 12)) + child, _ = refinement.nest_patch(cfg, parent, (8, 8, 40, 40) if large else (8, 8, 16, 16)) + middle = replace(cfg, nx=parent.fine_shape[1], ny=parent.fine_shape[0], viscosity=.2) + return child, refinement.prepare_fine_boundary(middle, child, objects, jets=jets, + parent_origin=parent.fine_origin, parent_spacing=.5) + + def test_root_frame_matches_independent_parent_frame_projection(self): + from CelerisLab.body.jets import Jet, Parallel, Tube, jet_link_shapes + jet = Jet(Tube((6.5, 8.), 0., 3., .3), Parallel(0.), amplitude=.002) + original = self.body() + child, b = self.prepare([original], [jet]) + middle = LBMConfig(nx=24, ny=24, viscosity=.2) + transformed = self.body((12., 12.), 1.6) + local_jet = Jet(Tube((9., 12.), 0., 6., .6), Parallel(0.), amplitude=.002) + other = refinement.prepare_fine_boundary(middle, child, [transformed], jets=[local_jet]) + for a, c in zip(b.curved_host, other.curved_host): + np.testing.assert_array_equal(a, c) + np.testing.assert_array_equal(b.flags, other.flags) + np.testing.assert_allclose(b.curved_hits, 2 + .5 * other.curved_hits, atol=0, rtol=0) + for a, c in zip(jet_link_shapes([jet], b.curved_hits), + jet_link_shapes([local_jet], other.curved_hits)): + np.testing.assert_allclose(a, c, atol=1e-14, rtol=0) + self.assertEqual(len(b.curved_host[0]), 60) + self.assertAlmostEqual(child.tau_fine, .5 + 12 * .1) + self.assertEqual(original.geometry.center, (8., 8.)) + self.assertEqual(original.obj_id, 99) + + def test_invalid_body_and_jet_contracts(self): + from CelerisLab.body.jets import Jet, Parallel, Tube + from CelerisLab.body.objects import SurfaceActuation + jet = Jet(Tube((6.5, 8.), 0., 3., .3), Parallel(0.), amplitude=.002) + for body in (self.body((6., 8.)), self.body((6.9, 8.)), self.body(radius=2.1)): + with self.assertRaises(ValueError): + self.prepare([body]) + moving = self.body() + moving.state.vx = .01 + with self.assertRaisesRegex(ValueError, "stationary"): + self.prepare([moving]) + with self.assertRaises(ValueError): + self.prepare([self.body()], [jet, jet]) + empty = Jet(Tube((0., 0.), 0., 1., .1), Parallel(0.), amplitude=0.) + with self.assertRaises(ValueError): + self.prepare([self.body()], [empty]) + spec = SurfaceActuation(.7, .85, .002) + with self.assertRaisesRegex(ValueError, "SurfaceActuation"): + self.prepare([self.body(spec=spec)], [spec.as_jet((8., 8.), .8)]) + + def test_one_jet_two_real_bodies_same_caller_id_rejected(self): + from CelerisLab.body.jets import Jet, Parallel, Tube + bodies = [self.body((12., 14.), 1., 7), self.body((17., 14.), 1., 7)] + _, boundary = self.prepare(bodies, large=True) # establish legal geometry first + self.assertEqual(np.unique(boundary.curved_host[3]).tolist(), [0, 1]) + for amplitude in (0., .002): + jet = Jet(Tube((10., 14.), 0., 9., .4), Parallel(0.), amplitude=amplitude) + with self.assertRaisesRegex(ValueError, "one jet may drive only one body"): + self.prepare(bodies, [jet], large=True) + + 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 index 8cb873d..ecc91c2 100644 --- a/tests/validation/multigrid/test_refinement_three_level.py +++ b/tests/validation/multigrid/test_refinement_three_level.py @@ -10,8 +10,10 @@ Run CPU and GPU classes in separate processes. GPU gates (frozen before any run) 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. """ +from contextlib import nullcontext +from dataclasses import replace import unittest -from unittest.mock import patch +from unittest.mock import Mock, patch import numpy as np import pycuda.driver as cuda @@ -28,6 +30,32 @@ 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)) +def child_body(spec=None): + from CelerisLab.body.geometry.circle import CircleGeometry + from CelerisLab.body.objects import SimObject + return SimObject(99, geometry=CircleGeometry(8., 8., .8), center=(8., 8.), radius=.8, + surface_actuation=spec) + + +def child_jet(speed=.002): + from CelerisLab.body.jets import Jet, Tube, Parallel + return Jet(Tube((6.5, 8.), 0., 3., .3), Parallel(0.), amplitude=speed) + + +def mocked_child(): + from CelerisLab.cuda import compiler_v2 as compiler + context = Mock(activate=nullcontext, sm_arch="sm_70") + module = Mock(get_function=Mock(side_effect=lambda name: Mock(name=name))) + module.get_global.return_value = (Mock(), 20) + with patch.object(cuda, "Stream"), patch.object(compiler, "generate_config"), \ + patch.object(compiler, "compile_kernel"), patch.object(compiler, "cleanup_build_dir"), \ + patch.object(compiler, "load_module", return_value=module), \ + patch.object(cuda, "mem_alloc", side_effect=lambda size: Mock()), \ + patch.object(cuda, "memcpy_htod"): + return FixedPatchStepper(CFG, L1, context=context, child_objects=[child_body()], + child_jets=[child_jet()], **KW) + + class ThreeLevelCPU(unittest.TestCase): def test_rejections_before_cuda(self): context = object() # never activated: every case must fail first @@ -37,11 +65,92 @@ class ThreeLevelCPU(unittest.TestCase): for kw in bad: with self.subTest(kw=kw), self.assertRaises(ValueError): FixedPatchStepper(CFG, L1, context=context, **kw) + for kw in (dict(child_objects=[child_body()]), dict(child_jets=[child_jet()]), + dict(KW, child_jets=[child_jet()])): + with 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) + def test_four_halo_bfl_steps_and_fresh_middle_feeds(self): + run = mocked_child() + run._status = "ready" + f, g = run._levels[1:] + b = g.boundary + events, feeds = [], [] + advance, gather = run._advance_child, run._gather_halo + def record_feed(buffer, **kw): + self.assertIs(buffer, (f.current, f.spare)[len(feeds)]) + feeds.append(buffer) + advance(buffer, **kw) + def record_halo(level, buffer, **kw): + if level is g: + events.append(("halo", buffer)) + gather(level, buffer, **kw) + b.kernel.side_effect = lambda buffer, *a, **k: events.append(("bfl", buffer)) + def record_step(src, dst, halo, flags, kind, **kw): + self.assertIs(flags, b.flags) + self.assertIs(kind, g.kind) + events.append(("step", src)) + g.kernel.side_effect = record_step + with patch.object(run, "_advance_child", side_effect=record_feed), \ + patch.object(run, "_gather_halo", side_effect=record_halo): + run.step() + self.assertEqual(feeds, [f.current, f.spare]) + expected = [(kind, buffer) for buffer in (g.current, g.spare, g.current, g.spare) + for kind in ("halo", "bfl", "step")] + self.assertEqual(events, expected) + self.assertGreater(len(np.unique(b.tables["block"])), 1) + # Negative control: reusing L1's old source at the midpoint must fail. + feeds.clear() + with patch.object(run, "_advance_child", side_effect=lambda buffer, **kw: + record_feed(f.current, **kw)): + with self.assertRaises(AssertionError): + run.step() + run.close() + + def test_child_update_isolated_and_failure_atomic(self): + run = mocked_child() + b = run._levels[2].boundary + before = b.tables["surface_velocity"].copy() + with self.assertRaises(IndexError): + run.set_jet_amplitude(0, -.002) # default stays L1 + for level in (0, 3, True): + with self.assertRaises(ValueError): + run.set_jet_amplitude(0, -.002, level=level) + with patch.object(cuda, "memcpy_htod") as upload: + run.set_jet_amplitude(0, -.002, level=2) + self.assertEqual(upload.call_count, 1) + self.assertIs(upload.call_args.args[0], b.device["surface_velocity"]) + np.testing.assert_array_equal(b.tables["surface_velocity"], -before) + self.assertEqual(run._jets, ()) + with patch.object(cuda, "memcpy_htod", side_effect=RuntimeError("copy failed")): + with self.assertRaisesRegex(RuntimeError, "copy failed"): + run.set_jet_amplitude(0, .002, level=2) + self.assertEqual(run._status, "failed") + self.assertEqual(b.jets[0].amplitude, -.002) + np.testing.assert_array_equal(b.tables["surface_velocity"], -before) + run.close() + + def test_child_public_masks_include_solids(self): + from CelerisLab.lbm.descriptors import SOLID + run = mocked_child() + solid = run.child_boundary.flags.reshape(run.child_patch.fine_shape) & SOLID != 0 + self.assertTrue(solid.any()) + self.assertFalse(run._child_leaf[solid].any()) + arrays = [np.ones(level.shape, np.float32) for level in run._levels] + run._mask_nonleaves(*arrays) + self.assertTrue(np.isnan(arrays[2][:, solid]).all()) + with patch.object(cuda, "memcpy_htod") as upload: + run.initialize(*arrays) + sent = upload.call_args_list[-1].args[1] + dense = run.child_layout.to_dense(sent) + self.assertTrue(np.isnan(dense[:, ~run._child_leaf]).all()) + run.close() + + class ThreeLevelGPU(unittest.TestCase): def test_equilibrium_is_kept_on_every_level(self): from CelerisLab.cuda.context import CudaContext @@ -85,5 +194,67 @@ class ThreeLevelGPU(unittest.TestCase): self.assertTrue(np.isnan(rg[:, p2.fine_kind == OUTER_GHOST]).all()) +class ThreeLevelBodyGPU(unittest.TestCase): + """Frozen gates; CPU mocks do not establish these numerical results.""" + def run_case(self, mode, block, *, jets=(), spec=None, switch=False): + from CelerisLab.cuda.context import CudaContext + from CelerisLab.lbm.descriptors import SOLID + cfg = CFG if mode == "periodic" else replace( + CFG, velocity=0., wall_normal_velocity=.001, wall_slot_width=4, + y_wall_bc="on_site_velocity", top_wall_bc="reservoir") + with CudaContext(0) as ctx, FixedPatchStepper( + cfg, L1, context=ctx, coarse_boundary=mode, child_objects=[child_body(spec)], + child_jets=jets, **dict(KW, child_block=block)) as run: + run.initialize(*(state(*level.shape[1:], i + 1) for i, level in enumerate(run._levels))) + snapshots = [] + for stop in (0, 1, 2, 4, 8): + run.step(stop - run.step_count) + if stop == 4 and switch: + run.set_jet_amplitude(0, -.002, level=2) + public, stored = run.download(), raw(run) + leaves = [~run.patch.covered, run._fine_leaf, run._child_leaf] + if mode != "periodic": + leaves[0][[0, -1]] = False + for g, mask in zip(public, leaves): + self.assertTrue(np.isfinite(g[:, mask]).all()) + self.assertTrue(np.isnan(g[:, ~mask]).all()) + invalid = [run.patch.covered.copy(), run._middle_kind == OUTER_GHOST, + (run.child_patch.fine_kind == OUTER_GHOST) | + ((run.child_boundary.flags.reshape(run.child_patch.fine_shape) & SOLID) != 0)] + if mode != "periodic": + invalid[0][[0, -1]] = True + for g, mask in zip(stored, invalid): + self.assertTrue(np.isnan(g[:, mask]).all()) + snapshots.append(public + tuple(stored)) + return snapshots + + def test_layouts_static_driven_switched_and_zero(self): + for mode in ("periodic", "slot_reservoir"): + cases = [self.run_case(mode, (4, 4), jets=jets, switch=switch) + for jets, switch in (((), False), ([child_jet()], False), + ([child_jet()], True), ([child_jet(0.)], False))] + for a, b in zip(cases[0], cases[3]): + for x, y in zip(a, b): + np.testing.assert_array_equal(x, y) + self.assertFalse(np.array_equal(cases[1][-1][2], cases[2][-1][2], equal_nan=True)) + for block in ((8, 4), (4, 8)): + for expected, (jets, switch) in zip(cases, (((), False), ([child_jet()], False), + ([child_jet()], True), ([child_jet(0.)], False))): + actual = self.run_case(mode, block, jets=jets, switch=switch) + for a, b in zip(actual, expected): + for x, y in zip(a, b): + np.testing.assert_array_equal(x, y) + + def test_legacy_surface_equivalence(self): + from CelerisLab.body.objects import SurfaceActuation + spec = SurfaceActuation(.7, .85, .002) + for mode in ("periodic", "slot_reservoir"): + legacy = self.run_case(mode, (4, 4), spec=spec)[-1][2] + driven = self.run_case(mode, (4, 4), jets=[spec.as_jet((8., 8.), .8)])[-1][2] + mask = np.isfinite(legacy) + np.testing.assert_array_equal(mask, np.isfinite(driven)) + self.assertLessEqual(np.max(np.abs(legacy[mask] - driven[mask])), 1e-6) + + if __name__ == "__main__": unittest.main()