From 5304525e394897610e16664bbfc2cceec576c966 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 15 Oct 2026 05:06:19 +0000 Subject: [PATCH] Add fixed-patch circle example and README scope section examples/fixed_patch_circle.py: 32x24 coarse, patch (8,6,24,18), circle (16,12) R2.15 in the cell-edge frame with the tested 0.7/0.85 aperture; --mode periodic|slot_reservoir, --speed (+/-, 0 passive), optional --fine-block BY BX, default 128 coarse steps; rest-equilibrium postcollision initialization, caller-owned context, fluid-only leaf mass/density/speed printed every 32 steps; no output files. README: experimental backend usage and scope (limits, cell-edge frame, whole-allocation divisibility, fine_block as AMR infrastructure without speedup claim, unavailable features). No solver changes. Co-Authored-By: Claude Opus 5.5 --- README.md | 42 ++++++++++++++ examples/fixed_patch_circle.py | 102 +++++++++++++++++++++++++++++++++ 2 files changed, 144 insertions(+) create mode 100644 examples/fixed_patch_circle.py diff --git a/README.md b/README.md index 115172c..e25d9f2 100644 --- a/README.md +++ b/README.md @@ -91,6 +91,48 @@ sim.run(100, stream=stream, upload_act=False, sync_obs=False) stream.synchronize() obs = sim.read_bodies(stream=stream) # sync already done, just read pinned buffer +### Fixed 2:1 refinement patch (experimental backend) + +`FixedPatchStepper` advances one fixed 2:1 refinement patch inside a coarse +grid, optionally with stationary circles in the fine region and a prescribed +normal surface velocity (blowing/suction). It is a standalone backend, not part +of `Simulation`. Runnable demo: `examples/fixed_patch_circle.py`. + +```python +from CelerisLab.body.geometry.circle import CircleGeometry +from CelerisLab.body.objects import SimObject, SurfaceActuation +from CelerisLab.config import LBMConfig +from CelerisLab.cuda.context import CudaContext +from CelerisLab.lbm.refinement_stepper import FixedPatchStepper + +cfg = LBMConfig(nx=32, ny=24, viscosity=.1) +body = SimObject(0, geometry=CircleGeometry(16.0, 12.0, 2.15), center=(16.0, 12.0), radius=2.15, + surface_actuation=SurfaceActuation(0.7, 0.85, 0.002)) # angle, half-width, speed +with CudaContext(0) as ctx, FixedPatchStepper(cfg, (8, 6, 24, 18), context=ctx, + fine_objects=[body], fine_block=(16, 20)) as run: + run.initialize(coarse_g, fine_g) # (9, ny, nx) postcollision populations per level + run.step(128) # coarse steps; fine level takes two per step + coarse_g, fine_g = run.download() # NaN on covered, ghost and solid cells +``` + +Scope and limits: +- D2Q9, SRT, FP32, double buffer, unshifted, no LES or forcing; one fixed 2:1 patch. +- `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)`), + unlike `Simulation`'s integer-node frame. Surface speeds are in coarse lattice units. +- Bodies must be stationary circles whose whole boundary stencil lies in the fine + active region, clear of the coarse/fine transfer cells; violations are rejected + before any GPU allocation. +- `fine_block=(by, bx)` stores the fine level in fixed-size blocks; `by`/`bx` must + divide the WHOLE fine allocation, i.e. `2*(y1-y0)+8` by `2*(x1-x0)+8` including the + four-cell ghost band. It gives results identical to dense storage and is + infrastructure for future block AMR, not a speedup. +- Not available: per-link surface exchange/force readout in blocked mode, + checkpoints, moving bodies, multiple patches or dynamic refinement, `Simulation` + integration. Validation so far is storage/operator equivalence and short runs, + not accuracy, convergence or long-time stability. + ## Installation ### Prerequisites diff --git a/examples/fixed_patch_circle.py b/examples/fixed_patch_circle.py new file mode 100644 index 0000000..ee0f034 --- /dev/null +++ b/examples/fixed_patch_circle.py @@ -0,0 +1,102 @@ +"""Public API demo: fixed 2:1 patch with a stationary, optionally actuated circle. + +From this checkout, after checking the selected physical GPU is idle: + PYTHONPATH="$PWD/src" CUDA_VISIBLE_DEVICES=0 python examples/fixed_patch_circle.py + ... --speed -0.002 --fine-block 16 20 --mode slot_reservoir + +A short smoke run (default 128 coarse steps): it shows setup, stepping and +fluid-only diagnostics, not a developed flow, accuracy or long-time stability. +Normal blowing/suction injects/extracts mass, so leaf mass is expected to change. +No files, plots or checkpoints are written. +""" +import argparse +from pathlib import Path + +import numpy as np + +import CelerisLab +from CelerisLab.body.geometry.circle import CircleGeometry +from CelerisLab.body.objects import SimObject, SurfaceActuation +from CelerisLab.config import LBMConfig +from CelerisLab.cuda.context import CudaContext +from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY +from CelerisLab.lbm.equilibrium import compute_feq_d2q9 +from CelerisLab.lbm.refinement_stepper import FixedPatchStepper + +NX, NY, BOUNDS = 32, 24, (8, 6, 24, 18) # coarse cells; patch is half-open (x0, y0, x1, y1) +CENTER, RADIUS = (16.0, 12.0), 2.15 # coarse CELL-EDGE frame: cell I spans [I, I+1) + + +def config(mode): + if mode == "periodic": # x and y periodic coarse level + return LBMConfig(nx=NX, ny=NY, viscosity=.1) + # Bottom on-site velocity slot, top fixed-density reservoir, x periodic. + return LBMConfig(nx=NX, ny=NY, viscosity=.1, velocity=0, wall_normal_velocity=.001, + wall_slot_width=8, y_wall_bc="on_site_velocity", top_wall_bc="reservoir") + + +def leaf_diagnostics(coarse, fine): + """Fluid leaves only (downloads are NaN on covered/ghost/solid cells). + + Mass counts a fine cell as 1/4 of a coarse cell (2D, 2:1).""" + out = {} + for name, g, area in (("coarse", coarse, 1.0), ("fine", fine, .25)): + leaf = np.isfinite(g).all(axis=0) + f = g[:, leaf].astype(np.float64) + rho = f.sum(axis=0) + ux = (np.array(D2Q9_EX)[:, None] * f).sum(axis=0) / rho + uy = (np.array(D2Q9_EY)[:, None] * f).sum(axis=0) / rho + out[name] = (area * rho.sum(), rho.min(), rho.max(), np.hypot(ux, uy).max()) + return out + + +def main(): + parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + parser.add_argument("--mode", choices=("periodic", "slot_reservoir"), default="periodic") + parser.add_argument("--speed", type=float, default=.002, + help="normal surface velocity, coarse lattice units (+ blows, - sucks, 0 passive)") + parser.add_argument("--fine-block", type=int, nargs=2, metavar=("BY", "BX"), + help="blocked fine storage; must divide the whole fine allocation (32 x 40 here)") + parser.add_argument("--steps", type=int, default=128) + args = parser.parse_args() + root = Path(__file__).resolve().parents[1] + if Path(CelerisLab.__file__).resolve() != root / "src/CelerisLab/__init__.py": + raise RuntimeError(f"Set PYTHONPATH to {root / 'src'}") + + # Aperture centred at angle 0.7 rad, half-width 0.85 rad (the tested fixture). + spec = SurfaceActuation(0.7, 0.85, args.speed) if args.speed else None + body = SimObject(0, geometry=CircleGeometry(*CENTER, RADIUS), center=CENTER, radius=RADIUS, + surface_actuation=spec) + cfg = config(args.mode) + fine_block = tuple(args.fine_block) if args.fine_block else None + # The caller owns the CUDA context; the stepper must be closed first (the + # nested with-statement does that). + with CudaContext(0) as context, FixedPatchStepper( + cfg, BOUNDS, context=context, coarse_boundary=args.mode, + fine_objects=[body], fine_block=fine_block) as run: + p = run.patch + # Rest equilibrium is already postcollision G, so no initial collision is + # needed. Non-leaf entries (covered, ghost, solid) are discarded by initialize. + rest = np.asarray(compute_feq_d2q9(1.0, 0.0, 0.0), np.float32)[:, None, None] + run.initialize(np.broadcast_to(rest, (9, *p.coarse_shape)), + np.broadcast_to(rest, (9, *p.fine_shape))) + storage = f"blocked {fine_block}" if fine_block else "dense" + print(f"mode={args.mode} fine storage={storage} speed={args.speed:+g} " + f"(fine allocation {p.fine_shape[0]}x{p.fine_shape[1]}, fine origin {p.fine_origin})") + start = leaf_diagnostics(*run.download()) + for _ in range(0, args.steps, 32): + run.step(min(32, args.steps - run.step_count)) + coarse, fine = run.download() # synchronized postcollision G, NaN on non-leaves + d = leaf_diagnostics(coarse, fine) + mass = d["coarse"][0] + d["fine"][0] + rho_lo = min(d["coarse"][1], d["fine"][1]) + rho_hi = max(d["coarse"][2], d["fine"][2]) + speed = max(d["coarse"][3], d["fine"][3]) + print(f"coarse step {run.step_count:4d}: leaf mass {mass:.9f} " + f"(change {mass - start['coarse'][0] - start['fine'][0]:+.3e}) " + f"rho [{rho_lo:.6f}, {rho_hi:.6f}] max |u| {speed:.3e}") + print("smoke run only: storage/stepping demonstration, not accuracy or stability evidence") + + +if __name__ == "__main__": + main()