Document bounded multigrid and add public remap example

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,
+  checkpoints, moving bodies, L1 bodies/Jets in a three-level hierarchy,
-  several patches per level with nesting, dynamic refinement, `Simulation`
+  arbitrary depth, automatic AMR or `Simulation` integration.
-  integration. Validation so far is storage/operator equivalence and short runs,
+  Storage/operator equivalence, remap and short-run gates do not establish
-  not accuracy, convergence or long-time stability. The new L2-body GPU gates
+  convergence or long-time physical accuracy. The earlier shear-decay result is
-  (layouts, masks, amplitude reversal, legacy equivalence) and independent
+  a research accuracy limitation, not a demonstrated solver bug. Pinned-body
-  every-substep numerical replay remain required acceptance checks; CPU tests
+  replay conditioned on actual L1 feeds is not an independent whole-body solver.
-  alone do not establish GPU numerical correctness.
+  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
 

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,
+It changes only a separated L2 rectangle set in a three-level periodic
-periodic hierarchy. L0/L1 topology, relaxation parameters, direction order,
+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
+actual state has no undeclared bodies or forcing, that all device work and old
-quiescent, and that inputs are the **current owned-leaf postcollision G**.
+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
+periodic three-level runs whose existing L2 uses
-`child_rectangles_in_parent`. Scalar or union L1 parents are supported; scalar
+`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
+finite leaves, exact masks and absolute counts. These are GPU execution gates,
-host. The separate independent FP64 capture/replay gate is max population
+not CPU-mock results. The separate body-free independent FP64 capture/replay
-error **3e-6**, with exact masks; remap mass/momentum allow FP32 roundoff.
+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
+These transfer/equivalence cases do not establish evolution accuracy. Historical
-accuracy. A 3e-6 FP32-vs-FP64 evolution gate requires a separate independent
+pinned-body FP64 replay uses actual L1 feeds: it checks conditional L2 evolution,
-BFL reference; the body-free replay gate above does not establish that claim.
+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.

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()