docs: checkpoint save/load scope and read_surface_exchange wording

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

docs/compatibility.md

@@ -0,0 +1,35 @@
+# Compatibility and API changes
+
+Changes that can make a previously accepted setup or call behave differently.
+Each entry names the commit, the affected entry points and the observable result.
+
+## 87c3a31: cut-link ownership is validated for every body
+
+`validate_cut_link_ownership` (`body/preprocess/flag_overlay.py`) now runs for
+passive bodies as well as Jet/aperture bodies. Every cut link must have an
+active fluid receiver (FLUID set, SOLID clear), and each `(receiver, direction)`
+population may have only one writer.
+
+Rejected with `ValueError`, before flags, populations or device buffers change:
+
+1. Overlapping or touching passive solids whose cut links write the same
+   population ("Boundary requires unique cut-link writers").
+2. A body whose cut-link receiver is a non-fluid cell, e.g. a channel wall row
+   or the inlet/outlet column ("Cut links must have active fluid receivers").
+
+Both were previously accepted silently for passive bodies; the interpolated
+boundary then had undefined results. The rules apply to the cut links of all
+bodies together, so a passive body overlapping a Jet body is rejected as well.
+The check runs in `Simulation.initialize`, runtime body sync (`sync_bodies`),
+checkpoint save and checkpoint load (`common/checkpoint.py`). Load validates
+before restoring anything, so a checkpoint whose saved geometry violates either
+rule no longer restores; rebuild the geometry with separated bodies.
+
+## 87c3a31: `Simulation.read_surface_exchange` accepts Jet-driven bodies
+
+`read_surface_exchange(id)` previously required `SurfaceActuation` on the body.
+It now also accepts bodies that own Jet aperture links, including Jets with zero
+or negative amplitude (the aperture test is geometric). Bodies with neither
+`SurfaceActuation` nor aperture links still raise
+`ValueError("Body has no surface actuation")`. The value is unchanged: the latest
+step's wall/aperture exchange, not total nozzle force.

docs/validation/dfg_re20.md

@@ -1,49 +0,0 @@
-# DFG 2D-1 (Re 20) against the official reference intervals
-
-`tests/validation/run_dfg_re20.py` runs the steady DFG 2D-1 case on the refined
-channel root (D40 root, D80 wide fine patch) and on the uniform D80 grid. Geometry,
-inflow, initial state, layout and observation are those of the pinned Re100 runner
-(`run_refined_dfg.py`, `re100_helpers.py`, both unchanged); only the lattice
-viscosity is set for Re 20 (root tau 0.56, fine 0.62, uniform 0.62, lattice speed 0.01).
-
-Label: "compared against the official reference bands", not a full benchmark
-certification: the outlet is pressure-NEQ (`neq_extrap`), not do-nothing, and the
-inlet is local Zou-He.
-
-## Quantities
-
-- Cd = 2 Fx / (rho U^2 D), Cl = 2 Fy / (rho U^2 D), interval means over each sample.
-- p_diff = p(0.15, 0.2) - p(0.25, 0.2), linear extrapolation to the exact surface
-  points (`re100_helpers.pressure_probes`). The helper returns dp / (rho U^2) in
-  scaled units; p_diff = dp * 0.04 for U = 0.2 m/s, rho = 1. Both are reported.
-- Official intervals (Schaefer & Turek 1996): Cd [5.57, 5.59], Cl [0.0104, 0.0110],
-  p_diff [0.1172, 0.1176].
-
-## Bounds and steady rule (frozen before any run)
-
-- End time t* = 60 (30 s physical): 240,000 refined root steps, 480,000 uniform steps.
-- One sample every t* = 0.5 (0.25 s) at the end of a synchronized root step; Cd and
-  Cl are interval means. Trailing window t* 50..60 inclusive (21 samples).
-- Steady rule (our engineering choice, not official): in the window, for each of Cd,
-  Cl, p_diff, max - min <= 10 % and |least-squares slope x window| <= 2 % of its
-  official interval width (range 2e-3 / 6e-5 / 4e-5; slope 4e-4 / 1.2e-5 / 8e-6).
-  Absolute limits, so the small Cl is never a divisor. Never relaxed after a run.
-- The verdict uses the window mean; the distance to the nearest interval edge is
-  reported for each quantity (negative outside).
-- Force sign: fluid on cylinder, x along the inflow, y up from the bottom wall; Cd > 0
-  and the official Cl > 0 (`cd_positive` is reported).
-- Not steady at t* = 60: status `NOT_STEADY`, no interval verdict. Health limits
-  (Mach <= 0.2, rho in [0.5, 1.5]) stop the run with `UNHEALTHY`.
-- Steady: `WITHIN_INTERVALS` or `OUTSIDE_INTERVALS`, with the per-quantity verdict.
-
-## Run
-
-```
-python tests/validation/run_dfg_re20.py --output <new-dir> --plan
-CUDA_VISIBLE_DEVICES=<n> python tests/validation/run_dfg_re20.py --case both \
-    --collision SRT --output <new-dir> --max-wall-seconds <cap>
-```
-
-The uniform D80 control separates method error from refinement error. `--plan` writes `plan.json` only. A GPU run writes `plan.json` and `result.json`
-(per case: lattice parameters, every sample, window ranges, status, verdict).
-Host contract tests: `tests/unit/test_dfg_re20.py`.

tests/unit/test_dfg_re20.py

@@ -1,79 +0,0 @@
-"""CPU-only contract tests for run_dfg_re20. Not CUDA or DFG acceptance evidence."""
-
-import sys
-from pathlib import Path
-import unittest
-
-sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "validation"))
-import run_dfg_re20 as v  # noqa: E402
-import run_refined_dfg as rr  # noqa: E402
-
-h, _ = rr.helpers(rr.DEFAULT_SOURCE)
-
-
-def rows(values, tstar_end=v.TSTAR_MAX):
-    n = int(round(tstar_end / v.SAMPLE_TSTAR))
-    return [dict(tstar=(i + 1) * v.SAMPLE_TSTAR, **values(i)) for i in range(n)]
-
-
-class Re20Contract(unittest.TestCase):
-    def test_lattice_parameters_and_bounds(self):
-        root, uniform = v.lattice_parameters(40), v.lattice_parameters(80)
-        self.assertAlmostEqual(root["tau"], 0.56, places=12)
-        self.assertAlmostEqual(uniform["tau"], 0.62, places=12)
-        self.assertEqual((root["steps"], root["interval"]), (240000, 2000))
-        self.assertEqual((uniform["steps"], uniform["interval"]), (480000, 4000))
-        self.assertAlmostEqual(v.plan()["refined_fine_tau"], 0.62, places=12)
-
-    def test_config_only_changes_viscosity(self):
-        base = v._re100_config(h, 40, "SRT", "none")
-        cfg, body, p = v.config(h, 40, "SRT")
-        self.assertAlmostEqual(cfg["physics"]["viscosity"], 0.01 * 40 / 20, places=15)
-        self.assertEqual((p["re"], p["viscosity"]), (20.0, cfg["physics"]["viscosity"]))
-        base[0]["physics"]["viscosity"] = cfg["physics"]["viscosity"]
-        self.assertEqual(base[0], cfg)
-        self.assertEqual(base[1], body)
-        with self.assertRaisesRegex(ValueError, "no Jets"):
-            v.config(h, 40, "MRT", "zero")
-
-    def test_units(self):
-        s = v.sample([2.0, -0.5], 4, 40, 3.0)
-        den = rr.SPEED**2 * 40 / 2
-        self.assertAlmostEqual(s["cd"], 2.0 / 4 / den)
-        self.assertAlmostEqual(s["cl"], -0.5 / 4 / den)
-        self.assertAlmostEqual(s["p_diff"], 3.0 * 0.04)
-        self.assertEqual(s["dp_scaled"], 3.0)
-
-    def test_steady_and_verdict(self):
-        inside = dict(cd=5.58, cl=0.0106, p_diff=0.1174)
-        out = v.summarize(rows(lambda i: inside))
-        self.assertEqual((out["status"], out["steady"], out["window_samples"]),
-                         ("WITHIN_INTERVALS", True, 21))
-        self.assertTrue(out["cd_positive"])
-        self.assertAlmostEqual(out["edge_distance"]["cl"], 0.0002)
-        self.assertTrue(all(out["verdict"].values()))
-        low = v.summarize(rows(lambda i: dict(inside, cl=0.0100)))
-        self.assertEqual(low["status"], "OUTSIDE_INTERVALS")
-        self.assertFalse(low["verdict"]["cl"])
-
-    def test_slow_drift_gives_no_verdict(self):
-        # Range 2e-5 passes (limit 6e-5); slope * window 2e-5 exceeds 1.2e-5.
-        out = v.summarize(rows(lambda i: dict(cd=5.58, cl=0.0106 + 1e-6 * i, p_diff=0.1174)))
-        self.assertEqual((out["status"], out["verdict"]), ("NOT_STEADY", None))
-        self.assertLessEqual(out["ranges"]["cl"], out["range_limits"]["cl"])
-        self.assertGreater(out["drifts"]["cl"], out["slope_limits"]["cl"])
-
-    def test_noise_range(self):
-        quiet = v.summarize(rows(lambda i: dict(cd=5.58, cl=0.0106 + 2e-5 * (i % 2), p_diff=0.1174)))
-        self.assertTrue(quiet["steady"])  # within the range limit, no trend
-        noisy = v.summarize(rows(lambda i: dict(cd=5.58, cl=0.0106 + 1e-4 * (i % 2), p_diff=0.1174)))
-        self.assertEqual(noisy["status"], "NOT_STEADY")
-        self.assertGreater(noisy["ranges"]["cl"], noisy["range_limits"]["cl"])
-
-    def test_incomplete_run_gives_no_verdict(self):
-        out = v.summarize(rows(lambda i: dict(cd=5.58, cl=0.0106, p_diff=0.1174), tstar_end=30))
-        self.assertEqual((out["status"], out["verdict"]), ("INCOMPLETE", None))
-
-
-if __name__ == "__main__":
-    unittest.main()

tests/validation/run_dfg_re20.py

@@ -1,222 +0,0 @@
-"""Steady DFG 2D-1 (Re 20): refined and uniform runs against the official reference intervals.
-
-Not a full benchmark certification. Declared deviations from the official setup:
-pressure-NEQ outlet (``neq_extrap``) instead of do-nothing, local Zou-He inlet.
-Geometry, inflow profile, initial state, refined layout and observation reuse the
-pinned Re100 runner (``run_refined_dfg``) and ``re100_helpers`` unchanged; only the
-lattice viscosity is set for Re 20 (nu = U_lattice * D / 20 at the root level).
-
-Units. The helpers use scaled DFG units with mean inflow speed 1, so
-Cd = 2 Fx / (rho U^2 D) and Cl = 2 Fy / (rho U^2 D) are scale-free, and
-``pressure_probes`` returns dp / (rho U^2). The official p_diff = p(0.15, 0.2) -
-p(0.25, 0.2) uses U = 0.2 m/s, rho = 1, so p_diff = dp * 0.04. Both are reported;
-only p_diff is compared with the official interval. Forces use the existing sign
-convention: F is the force of the fluid on the cylinder, i.e. the refined raw
-exchange impulse (window ledger, ``total.impulse``) or the uniform ``read_force``,
-with x along the inflow and y up from the bottom wall. Cd must be positive, and the
-official Cl is positive (centre 0.20 below the centreline 0.205), so a flipped sign
-shows up as Cl outside its interval.
-
-Sampling. One sample every SAMPLE_TSTAR (t* = U t / D; 0.25 s), taken at the end of
-a completed root step with all levels synchronized; the step index is recorded. Cd
-and Cl are means over the sample interval (window exchange ledger refined, per-batch
-forces uniform), which filters acoustic and step-parity oscillation; p_diff is the
-end-of-interval value at the exact surface points.
-
-Steady rule (our engineering choice, frozen before any run; not official). Over the
-trailing WINDOW_TSTAR (t* 50..60 inclusive, 21 samples), for each of Cd, Cl and
-p_diff: range (max - min) <= 10 % and |least-squares slope * window| <= 2 % of the
-width of its official interval. Absolute limits, so the small Cl is never a divisor;
-the slope term bounds slow monotone drift that the range alone would pass. If the
-rule fails at TSTAR_MAX the run reports "not steady" and gives no interval verdict;
-limits are never relaxed after a run. The verdict uses the window mean, and the
-distance of each mean to its nearest interval edge is reported (negative outside).
-"""
-
-import argparse
-import json
-import math
-import os
-from pathlib import Path
-import sys
-import time
-from tempfile import TemporaryDirectory
-
-import numpy as np
-
-HERE = Path(__file__).resolve().parent
-sys.path.insert(0, str(HERE))
-import run_refined_dfg as rr  # noqa: E402
-
-RE = 20.0
-U_PHYS = 0.2                   # m/s mean inflow; rho = 1
-TSTAR_MAX = 60.0               # 30 s physical
-SAMPLE_TSTAR = 0.5             # 0.25 s
-WINDOW_TSTAR = 10.0            # trailing 5 s, 21 samples
-# Schaefer & Turek (1996), 2D-1 lower/upper bounds.
-BANDS = dict(cd=(5.57, 5.59), cl=(0.0104, 0.0110), p_diff=(0.1172, 0.1176))
-STEADY_RANGE = {k: (hi - lo) / 10 for k, (lo, hi) in BANDS.items()}
-STEADY_SLOPE = {k: (hi - lo) / 50 for k, (lo, hi) in BANDS.items()}
-DEVIATIONS = ("pressure-NEQ outlet (neq_extrap), not do-nothing", "local Zou-He inlet")
-
-_re100_config = rr.config
-
-
-def config(h, diameter, collision="SRT", jet_scenario="none"):
-    """The pinned Re100 configuration with the lattice viscosity set for Re 20."""
-    if jet_scenario != "none":
-        raise ValueError("DFG 2D-1 runs have no Jets")
-    cfg, body, p = _re100_config(h, diameter, collision, jet_scenario)
-    nu = p["mean_speed"] * p["diameter"] / RE
-    cfg["physics"]["viscosity"] = nu
-    return cfg, body, dict(p, viscosity=nu, re=RE)
-
-
-def lattice_parameters(diameter, speed=rr.SPEED):
-    nu = speed * diameter / RE
-    return dict(viscosity=nu, tau=0.5 + 3 * nu, steps=math.ceil(TSTAR_MAX * diameter / speed),
-                interval=round(SAMPLE_TSTAR * diameter / speed))
-
-
-def sample(impulse, n, diameter, dp):
-    den = rr.SPEED**2 * diameter / 2
-    return dict(cd=float(impulse[0] / n / den), cl=float(impulse[1] / n / den),
-                dp_scaled=float(dp), p_diff=float(dp) * U_PHYS**2)
-
-
-def summarize(rows):
-    """Steady rule over the trailing window, then the interval verdict (or none)."""
-    if not rows or rows[-1]["tstar"] < TSTAR_MAX - 1e-9:
-        return dict(status="INCOMPLETE", steady=False, verdict=None)
-    window = [r for r in rows if r["tstar"] >= TSTAR_MAX - WINDOW_TSTAR - 1e-9]
-    t = np.array([r["tstar"] for r in window])
-    ranges = {k: max(r[k] for r in window) - min(r[k] for r in window) for k in BANDS}
-    drifts = {k: abs(np.polyfit(t, [r[k] for r in window], 1)[0]) * WINDOW_TSTAR for k in BANDS}
-    steady = all(ranges[k] <= STEADY_RANGE[k] and drifts[k] <= STEADY_SLOPE[k] for k in BANDS)
-    final = {k: float(np.mean([r[k] for r in window])) for k in BANDS}
-    edge = {k: min(final[k] - BANDS[k][0], BANDS[k][1] - final[k]) for k in BANDS}
-    result = dict(window_samples=len(window), ranges=ranges, drifts=drifts,
-                  range_limits=STEADY_RANGE, slope_limits=STEADY_SLOPE, steady=steady,
-                  window_mean=final, edge_distance=edge, bands=BANDS,
-                  deviations=DEVIATIONS, cd_positive=final["cd"] > 0)
-    if not steady:
-        return dict(result, status="NOT_STEADY", verdict=None)
-    inside = {k: BANDS[k][0] <= final[k] <= BANDS[k][1] for k in BANDS}
-    return dict(result, status="WITHIN_INTERVALS" if all(inside.values()) else "OUTSIDE_INTERVALS",
-                verdict=inside)
-
-
-def healthy(row):
-    return row["max_mach"] <= 0.2 and row["rho_min"] >= 0.5 and row["rho_max"] <= 1.5
-
-
-def run_case(a, h, case):
-    from contextlib import ExitStack
-    from CelerisLab.cuda.context import CudaContext
-    from CelerisLab.lbm.descriptors import FLUID
-    from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
-    from CelerisLab.simulation import Simulation
-
-    refined = case == "refined"
-    cfgdict, bodydict, p = config(h, 40 if refined else 80, a.collision)
-    lp = lattice_parameters(p["diameter"])
-    rows = []
-    with ExitStack() as stack:
-        if refined:
-            cfg, pair, body, host = rr.host_geometry(a.source, h, a.collision)
-            context = stack.enter_context(CudaContext(0))
-            backend = stack.enter_context(FixedPatchStepper(
-                cfg, rr.PLACEMENTS["wide"], context=context, coarse_boundary="channel",
-                base_block=rr.BASE_BLOCK, fine_block=rr.FINE_BLOCK, fine_objects=(body,),
-                observe_raw_exchange=True, observe_window_exchange=True))
-            from CelerisLab.lbm.channel import channel_masks
-            from CelerisLab.lbm.refinement import ACTIVE
-            from CelerisLab.lbm.descriptors import SOLID
-            masks = [channel_masks(pair.coarse_shape)[0] & ~pair.covered,
-                     (pair.fine_kind == ACTIVE)
-                     & ((host.flags.reshape(pair.fine_shape) & SOLID) == 0)]
-            backend.initialize(rr.initial_state(pair.coarse_shape, (0.0, 0.0), 1.0, root=True),
-                               rr.initial_state(pair.fine_shape, pair.fine_origin, 0.5))
-            probe_p = dict(p, center=[(80.5 - pair.fine_origin[0]) / 0.5 - 0.5,
-                                      (81.0 - pair.fine_origin[1]) / 0.5 - 0.5], radius=40.0)
-        else:
-            tmp = Path(stack.enter_context(TemporaryDirectory()))
-            (tmp / "lbm.json").write_text(json.dumps(cfgdict))
-            (tmp / "body.json").write_text(json.dumps(bodydict))
-            sim = stack.enter_context(Simulation(str(tmp / "lbm.json"), str(tmp / "body.json"),
-                                                 device_id=0))
-            sim.initialize()
-            probe_p = p
-            fluid = (sim.field.flag.reshape(p["ny"], p["nx"]) & FLUID) != 0
-        begin = 0
-        while begin < lp["steps"]:
-            rr.check_deadline(a.deadline)
-            n = min(lp["interval"], lp["steps"] - begin)
-            if refined:
-                impulse = rr.advance_refined_interval(backend, n, a.deadline, "window")
-                fields = [rr.macros(g) for g in backend.download()]
-                rho, ux, uy = fields[1]
-                fluid = masks[1]
-                vals = [tuple(v[m] for v in f) for m, f in zip(masks, fields)]
-            else:
-                impulse, done = np.zeros(2), 0
-                while done < n:
-                    batch = min(100, n - done)
-                    sim.run(batch)
-                    impulse += sim.read_force(0).astype(np.float64) * batch
-                    done += batch
-                m = sim.get_macroscopic()
-                rho, ux, uy = [np.asarray(m[k], np.float64).reshape(p["ny"], p["nx"])
-                               for k in ("rho", "ux", "uy")]
-                vals = [(rho[fluid], ux[fluid], uy[fluid])]
-            begin += n
-            row = dict(step=begin, tstar=begin * rr.SPEED / p["diameter"],
-                       **sample(impulse, n, p["diameter"],
-                                h.pressure_probes(rho, fluid, probe_p)["dp"]),
-                       rho_min=float(min(v[0].min() for v in vals)),
-                       rho_max=float(max(v[0].max() for v in vals)),
-                       max_mach=float(max(np.sqrt(3 * (v[1]**2 + v[2]**2)).max() for v in vals)))
-            rr.ensure_finite(row)
-            rows.append(row)
-            if not healthy(row):
-                break
-    summary = summarize(rows) if healthy(rows[-1]) else dict(status="UNHEALTHY", steady=False,
-                                                             verdict=None)
-    return dict(case=case, collision=a.collision, lattice=lp, rows=rows, **summary)
-
-
-def plan(collision="SRT"):
-    return dict(re=RE, collision=collision, tstar_max=TSTAR_MAX, sample_tstar=SAMPLE_TSTAR,
-                window_tstar=WINDOW_TSTAR, bands=BANDS, range_limits=STEADY_RANGE, slope_limits=STEADY_SLOPE,
-                deviations=DEVIATIONS, pressure_scale=U_PHYS**2,
-                refined_root=lattice_parameters(40), refined_fine_tau=0.5 + 3 * 2 * rr.SPEED * 40 / RE,
-                uniform=lattice_parameters(80),
-                label="compared against the official reference bands (non-canonical neq_extrap outlet, zou_he_local inlet); not a full benchmark certification")
-
-
-def main():
-    parser = argparse.ArgumentParser(description=__doc__.splitlines()[0])
-    parser.add_argument("--case", choices=["refined", "uniform", "both"], default="refined")
-    parser.add_argument("--collision", choices=["SRT", "MRT"], default="SRT")
-    parser.add_argument("--source", type=Path, default=rr.DEFAULT_SOURCE)
-    parser.add_argument("--output", type=Path, required=True)
-    parser.add_argument("--plan", action="store_true", help="Write the frozen plan; no CUDA")
-    parser.add_argument("--max-wall-seconds", type=float)
-    a = parser.parse_args()
-    a.source = a.source.resolve()
-    a.output.mkdir(parents=True, exist_ok=False)
-    (a.output / "plan.json").write_text(json.dumps(plan(a.collision), indent=2, allow_nan=False))
-    if a.plan:
-        return
-    if not a.max_wall_seconds or a.max_wall_seconds <= 0 or not os.environ.get("CUDA_VISIBLE_DEVICES"):
-        raise SystemExit("GPU run needs --max-wall-seconds > 0 and explicit CUDA_VISIBLE_DEVICES")
-    h, _ = rr.helpers(a.source)
-    rr.config = config  # host_geometry builds its root configuration through this name
-    a.deadline = time.monotonic() + a.max_wall_seconds
-    cases = ["refined", "uniform"] if a.case == "both" else [a.case]
-    results = [run_case(a, h, case) for case in cases]
-    (a.output / "result.json").write_text(json.dumps(results, indent=2, allow_nan=False))
-
-
-if __name__ == "__main__":
-    main()