From f75d543586315b6272cfdf23046aa7d4550f5fbe Mon Sep 17 00:00:00 2001 From: TRON Date: Fri, 16 Oct 2026 11:54:43 -0700 Subject: [PATCH] Document bounded multigrid and add public remap example --- README.md | 36 ++++++--- docs/synchronized_leaf_remap.md | 51 +++++++++---- examples/multigrid_remap.py | 126 ++++++++++++++++++++++++++++++++ 3 files changed, 190 insertions(+), 23 deletions(-) create mode 100644 examples/multigrid_remap.py diff --git a/README.md b/README.md index ee0b029..c3c0b85 100644 --- a/README.md +++ b/README.md @@ -123,7 +123,7 @@ with CudaContext(0) as ctx, FixedPatchStepper(cfg, (8, 6, 24, 18), context=ctx, ``` Scope and limits: -- D2Q9, SRT, FP32, double buffer, unshifted, no LES or forcing; fixed 2:1 patch(es), two levels (or one nested third). +- D2Q9, SRT, FP32, double buffer, unshifted, no LES or forcing; fixed 2:1 patches, two or three levels with separated rectangles. - `coarse_boundary="periodic"` (default) or `"slot_reservoir"` (bottom slot, top reservoir; patch needs `y0 >= 4`, `y1 <= ny-4`). Periodic mode ignores the wall/slot config fields. - Object centers/radii use the coarse **cell-edge** frame (cell `I` spans `[I, I+1)`), @@ -144,7 +144,8 @@ Scope and limits: not matter. `fine_objects`/`jets` work as for one patch: each circle (same coarse cell-edge frame) must lie wholly inside one patch's active region with its boundary stencil clear of transfer cells; a jet may not select links of two - different bodies, even across patches. No wall-touching patches or nesting yet. + different bodies, even across patches. Wall-touching patches are unsupported. + Three-level nesting uses the explicit child coordinates described below. ```python cfg = LBMConfig(nx=32, ny=16, viscosity=.1) @@ -164,7 +165,10 @@ Scope and limits: 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. Optional stationary circles and Jets + steps = 4 L2 steps. This scalar `child=` form uses one L1 patch. + For separated L1/L2 sets, use `child_rectangles_in_parent`; see + [coordinate frames and supported layouts](docs/fixed_three_level_multipatch.md). + 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 @@ -192,14 +196,26 @@ Scope and limits: run.step(4) coarse_g, l1_g, l2_g = run.download() # solids/covered/nonleaves are NaN ``` +- Manual synchronized L2 replacement: `replacement = current.remap_child(rectangles, + child_objects=original_objects)` preserves the absolute step count and live Jet + amplitudes. L0/L1 and whole body-containing L2 rectangles stay fixed; separated + fluid rectangles may change. The source remains independently owned. See + [ownership, pinning and transfer rules](docs/synchronized_leaf_remap.md). +- Runnable three-level example: `PYTHONPATH=src python examples/multigrid_remap.py`. + It uses two separated L1 regions, an L2 stationary Circle/Jet and an explicit + fluid-rectangle move, retaining original root-frame body definitions. Add + `--check-geometry` for a CPU-only check or `--verify-jet` for an optional + next-step comparison with an explicitly driven control. No flowfiles are saved. - Not available: per-link surface exchange/force readout in blocked mode, - 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. 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. + checkpoints, moving bodies, L1 bodies/Jets in a three-level hierarchy, + arbitrary depth, automatic AMR or `Simulation` integration. + Storage/operator equivalence, remap and short-run gates do not establish + convergence or long-time physical accuracy. The earlier shear-decay result is + a research accuracy limitation, not a demonstrated solver bug. Pinned-body + replay conditioned on actual L1 feeds is not an independent whole-body solver. + See [three-level validation scope](docs/fixed_three_level_multipatch.md) and + [remap validation scope](docs/synchronized_leaf_remap.md). The new usage example + has CPU geometry/syntax checks only until its GPU command is run on Horsea. ## Installation diff --git a/docs/synchronized_leaf_remap.md b/docs/synchronized_leaf_remap.md index 639cf53..bf6763e 100644 --- a/docs/synchronized_leaf_remap.md +++ b/docs/synchronized_leaf_remap.md @@ -1,8 +1,9 @@ # Synchronized leaf G remap and manual continuation This is a CPU state-transfer component, **not a runtime adaptive solver**. -It changes only a separated L2 rectangle set in a three-level, body-free, -periodic hierarchy. L0/L1 topology, relaxation parameters, direction order, +It changes only a separated L2 rectangle set in a three-level periodic +hierarchy, either body-free or with stationary Circles/Jets in unchanged pinned +L2 rectangles (see below). L0/L1 topology, relaxation parameters, direction order, physical frame and completed time remain unchanged. Existing geometry checks validate child support/nesting and separation. The CPU helper itself adds no CUDA, scheduler, allocator, boundary handling, collision or ghost refill. A bounded runtime adapter is described below. @@ -29,11 +30,12 @@ origin is retained. Physical lower edges are `origin + spacing*(x,y)` on L0, `origin + spacing*(parent.fine_origin + (x,y)/4)` on L2. Factory options `coarse_boundary` and `objects` must declare `"periodic"` and -an empty body collection (all levels). LBMConfig's channel settings do not +an empty L0/L1 body collection; use `child_objects`/`child_jets` for pinned L2 +bodies. LBMConfig's channel settings do not select this operator's boundary mode. Geometry metadata is factory-created; manual mutation/forging of metadata is unsupported. Caller certifies the -actual state has no bodies/forcing, that all device work and old readers are -quiescent, and that inputs are the **current owned-leaf postcollision G**. +actual state has no undeclared bodies or forcing, that all device work and old +readers are quiescent, and that inputs are the **current owned-leaf postcollision G**. Count/phase stamps check consistency, not freshness or active-buffer identity. `current_g` is three NumPy Q-major arrays `(9,ny,nx)`. For blocked inputs, @@ -73,8 +75,9 @@ PYTHONPATH=src python -m unittest discover -s tests/unit -p 'test_synchronized_r `replacement = current.remap_child(rectangles)` returns a new ready `FixedPatchStepper` at the **same absolute step_count**. It supports only -body-free periodic three-level runs whose existing L2 uses -`child_rectangles_in_parent`. Scalar or union L1 parents are supported; scalar +periodic three-level runs whose existing L2 uses +`child_rectangles_in_parent`, either body-free or under the pinned-body rules below. +Scalar or union L1 parents are supported; scalar `child=` is deliberately rejected. L0/L1, relaxation and block sizes stay fixed. No tagger, scheduler, in-place replacement, or general rollback is provided. @@ -105,9 +108,9 @@ and geometry, and candidate/source failure cleanup. GPU tests are predeclared in `tests/validation/multigrid/test_runtime_child_remap_gpu.py`: identity next step bitwise, shifted A→B→A with two intervals between changes, uniform moving and nonequilibrium states, exact CPU remap and overlap, both buffers reset, -finite leaves, exact masks and absolute counts. They were not run on the CPU -host. The separate independent FP64 capture/replay gate is max population -error **3e-6**, with exact masks; remap mass/momentum allow FP32 roundoff. +finite leaves, exact masks and absolute counts. These are GPU execution gates, +not CPU-mock results. The separate body-free independent FP64 capture/replay +gate is max population error **3e-6**, with exact masks; remap mass/momentum allow FP32 roundoff. Per-population conservation is not an advancement target under collision. ### Bounded pinned-body extension @@ -155,6 +158,28 @@ Multiple-body CPU coverage: `tests/unit/test_multiple_pinned_remap.py`. two independently updated Jets, variable fluid-region counts and add/move/remove immediate comparisons with CPU transfer. All nine solid-excluded area-weighted population integrals, exact overlap and nonleaf NaNs are checked. -These new GPU cases have **not been run here**. They do not establish evolution -accuracy. A 3e-6 FP32-vs-FP64 evolution gate requires a separate independent -BFL reference; the body-free replay gate above does not establish that claim. +These transfer/equivalence cases do not establish evolution accuracy. Historical +pinned-body FP64 replay uses actual L1 feeds: it checks conditional L2 evolution, +not an independent whole-body solver or whole-hierarchy accuracy. The body-free +replay gate above does not establish that stronger claim either. + +## Small runnable usage example + +[`examples/multigrid_remap.py`](../examples/multigrid_remap.py) keeps the original +root-frame Circle/Jet definitions, moves a fluid-only L2 rectangle, checks fluid +leaves with geometry-derived masks and retires both old/new owners before their +shared context. `--verify-jet` adds a temporary control solely to compare the +next step against an explicitly set live amplitude; normal use needs no control. + +```sh +# Host-only: no CUDA import or allocation. +PYTHONPATH=src python examples/multigrid_remap.py --check-geometry +# Horsea: select an approved idle physical GPU; it becomes logical device 0. +PYTHONPATH="$PWD/src" CUDA_VISIBLE_DEVICES=0 python examples/multigrid_remap.py --verify-jet +``` + +The new example has passed host geometry and syntax/import checks only; its GPU +run is pending Horsea. Expected scope: base counts 2 before/after remap, then 5; +finite positive-density fluid leaves, NaN nonleaves, and exact next-step control +agreement. This is a short public-API smoke run, not a physical accuracy test. +No flowfiles are written. diff --git a/examples/multigrid_remap.py b/examples/multigrid_remap.py new file mode 100644 index 0000000..24e1c1d --- /dev/null +++ b/examples/multigrid_remap.py @@ -0,0 +1,126 @@ +"""Three levels, separated regions and an explicit fluid-region remap. + +Run on an approved idle GPU (logical device 0): + PYTHONPATH=src CUDA_VISIBLE_DEVICES=0 python examples/multigrid_remap.py --verify-jet +Host-only geometry check, without CUDA: + PYTHONPATH=src python examples/multigrid_remap.py --check-geometry + +A short public-API smoke example, not an accuracy/stability benchmark. No files +are written. Edit the rectangles below; this is manual remapping, not auto AMR. +""" +import argparse +from contextlib import ExitStack + +import numpy as np + +from CelerisLab.body.geometry.circle import CircleGeometry +from CelerisLab.body.jets import Jet, Parallel, Tube +from CelerisLab.body.objects import SimObject +from CelerisLab.config import LBMConfig +from CelerisLab.lbm.equilibrium import compute_feq_d2q9 +from CelerisLab.lbm.refinement import FixedPatchSet2D +from CelerisLab.lbm.synchronized_remap import RemapTopology + +CFG = LBMConfig(nx=48, ny=32, viscosity=.1) # D2Q9/SRT/FP32, periodic +PARENTS = ((4, 4, 20, 28), (24, 4, 44, 28)) # L0 cell-edge rectangles +PIN = (16, 16, 32, 24) # L1 allocated-grid indices; whole body region stays fixed +BEFORE = (PIN, (56, 32, 72, 40)) +AFTER = (PIN, (60, 32, 76, 40)) # move only the separated fluid rectangle +LIVE_AMPLITUDE = -.004 + + +def geometry(): + # Retain these original L0 edge-frame definitions for every remap; never + # substitute the stepper's projected fine-local boundary objects. + objects = (SimObject(0, geometry=CircleGeometry(12., 10., .8), + center=(12., 10.), radius=.8),) + jets = (Jet(Tube((10.5, 10.), 0., 3., .3), Parallel(0.), amplitude=.002),) + parent = FixedPatchSet2D.from_config(CFG, PARENTS) + meshes = tuple(RemapTopology.from_config( + CFG, parent, rectangles, child_objects=objects, child_jets=jets) + for rectangles in (BEFORE, AFTER)) + if any(mesh.pinned_rectangles != (PIN,) for mesh in meshes): + raise RuntimeError("the body-containing rectangle must remain pinned") + return objects, jets, meshes + + +def diagnostics(run, mesh): + """Masks come from geometry, never from whether downloaded values are finite.""" + mass = 0. + for level, (g, leaf) in enumerate(zip(run.download(), mesh.leaf_masks)): + fluid = g[:, leaf].astype(np.float64) + rho = fluid.sum(axis=0) + if not np.isfinite(fluid).all() or np.any(rho <= 0): + raise RuntimeError(f"invalid L{level} fluid populations/density") + if not np.isnan(g[:, ~leaf]).all(): + raise RuntimeError(f"unexpected L{level} nonleaf values") + mass += rho.sum() / 4**level # area weights: 1, 1/4, 1/16 + print(f"absolute base step {run.step_count}: finite fluid leaves, mass={mass:.9f}") + # Driven boundary injects/extracts mass; mass need not stay constant in time. + + +def verify_live_jet(current, objects): + """Optional public-API check; normal continuation needs no control solver.""" + # Identity remap copies the same state. Explicitly setting the known live + # amplitude on the control detects a stale amplitude on the continuation. + # No private fields or constructor-time Jet snapshots inspect the solver. + with current.remap_child(AFTER, child_objects=objects) as control: + control.set_jet_amplitude(0, LIVE_AMPLITUDE, level=2) + current.step() + control.step() + for actual, expected in zip(current.download(), control.download()): + if not np.array_equal(actual, expected, equal_nan=True): + raise RuntimeError("live Jet continuation differs from explicitly set control") + print("live Jet amplitude preserved: next step equals explicitly set control") + + +def main(): + parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + parser.add_argument("--check-geometry", action="store_true") + parser.add_argument("--verify-jet", action="store_true", + help="compare one step with an explicitly driven control") + args = parser.parse_args() + objects, jets, (old_mesh, new_mesh) = geometry() + if args.check_geometry: + print("geometry OK: three levels, two L1 regions, fixed body pin, moved fluid region") + return + + # Lazy imports keep --check-geometry useful on a CPU-only machine. + from CelerisLab.cuda.context import CudaContext + from CelerisLab.lbm.refinement_stepper import FixedPatchStepper + + # ExitStack owns every old/new stepper, including on failure. All are closed + # before the shared context. remap_child never consumes/closes its source. + with CudaContext(0) as context, ExitStack() as owned: + old = owned.enter_context(FixedPatchStepper( + CFG, PARENTS, context=context, base_block=(4, 4), fine_block=(4, 4), + child_block=(4, 4), child_rectangles_in_parent=BEFORE, + child_objects=objects, child_jets=jets)) + rest = np.asarray(compute_feq_d2q9(1., 0., 0.), np.float32)[:, None, None] + old.initialize(*(np.broadcast_to(rest, (9, *mask.shape)) + for mask in old_mesh.leaf_masks)) + old.step(2) # one base step advances L1 twice and L2 four times + old.set_jet_amplitude(0, LIVE_AMPLITUDE, level=2) + diagnostics(old, old_mesh) + count = old.step_count + current = owned.enter_context(old.remap_child(AFTER, child_objects=objects)) + if current is old or current.step_count != count: + raise RuntimeError("remap must return a distinct continuation at the same time") + diagnostics(current, new_mesh) + + if args.verify_jet: + verify_live_jet(current, objects) # optional control; advances one step + else: + current.step() + if old.step_count != count or current.step_count != count + 1: + raise RuntimeError("old/new absolute step ownership changed") + old.close() # retire the old state only after successful replacement + current.step(2) + if current.step_count != count + 3: + raise RuntimeError("continuation lost its absolute step count") + diagnostics(current, new_mesh) + print("smoke example only: no accuracy or long-time stability claim") + + +if __name__ == "__main__": + main()