research: karman reproduce uses committed pinball ids after sync_bodies

add_body after initialize returns a pending negative id; sync_bodies commits
the new bodies in order, so the pinball ids are n0..n0+2 (f2f88c phase3 did
the same). Fixes the KeyError(-1) in the GPU smoke.

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

README.md

@@ -394,6 +394,12 @@ before `run()`, the force will be reset to zero. For the common usage pattern
 | `sim.load_checkpoint(path)` | Restore from HDF5 (config must match) |
 | `sim.snapshot()` / `sim.restore()` | In-memory field snapshot |
 
+In-memory snapshots save populations only. Restore keeps the current clock,
+controls, observations, geometry and flags; use it only while geometry is
+unchanged. This is not a restart or an FSI rollback. In particular, do not
+restore a snapshot taken before a translating circle changed position.
+This geometry restriction is a caller precondition, not a runtime check.
+
 #### Low-level access
 
 | Attribute | Description |

docs/moving_bodies.md

@@ -30,6 +30,10 @@ cells (M1 `MovingCircleCells`), moving-wall Bouzidi with wall velocity
   mass, eq_only, no_donor, crossings) accumulated since obs were zeroed;
   `obs`/`read_force` keep their contract (CME on the body).
 - `get_flags()` reads device flags back when moving circles exist.
+- `snapshot()` / `restore()` save and restore populations only, keeping the
+  current clock, controls, observations, geometry and flags. Use them only
+  with unchanged geometry; restoring across a translation is unsupported
+  and is not rejected at runtime. They do not provide FSI rollback.
 
 *Supported:* dense double buffer, D2Q9, FP32 unshifted, SRT/MRT, circles,
 spin. *Rejected before mutation:* EsoPull, ddf shifting, D3Q19, refinement /

docs/validation/viv_re150_predeclare.md

@@ -29,10 +29,17 @@ Claude, Grok and Astra at 1abfe954; implemented in research/viv_re150/.
 
 - D = 40, MRT, dense double buffer; U_lat = .04, nu = U D / 150 = .010667,
   tau = .532.
-- Box [0, 75 D] × [0, 60 D] (3000 × 2400) with the cylinder centre at
-  (30 D, 30 D): 30 D from the inlet, 45 D from the outlet, centred in y
-  (Kolahdouz §4.1.1). Qin's own box ([0, 40 D] × [0, 20 D], centre (10 D, 10 D))
-  is the target paper's, not ours.
+- Nominal 75 D × 60 D lattice allocation (3000 × 2400), based on the
+  Kolahdouz §4.1.1 layout. Coordinate clarification of the existing run:
+  the driver uses `(30 D, (ny - 1)/2)` = `(1200, 1199.5)` in lattice
+  coordinates, or `(30 D, 29.9875 D)`. The earlier `(30 D, 30 D)` wording
+  described the nominal box, not the exact executed centre.
+  Free-slip reflection planes are at `y = .5` and `y = ny - 1.5 = 2398.5`;
+  their separation is 2398 cells = 59.95 D, with the centre exactly midway,
+  29.975 D from each wall. Node counts are not physical wall-to-wall lengths.
+  This corrects the description, without changing the historical geometry
+  or results. Qin's own box ([0, 40 D] × [0, 20 D], centre (10 D, 10 D)) is
+  the target paper's, not ours.
 - Boundaries: uniform Zou–He inlet, NEQ outlet (rho = 1, backflow clamp),
   free-slip y walls. Kolahdouz uses a traction outlet: declared mismatch.
 - Qin's own VIV box (40 D × 20 D) is not used; a smaller box would be a

research/duetsch_re100_kc5/protocol.py

@@ -8,6 +8,7 @@ s/T in [0.25 + (k-1), 0.25 + k); cycle 0 is the startup quarter.
 from __future__ import annotations
 
 import math
+from functools import lru_cache
 
 import numpy as np
 
@@ -53,10 +54,16 @@ def bin_phase(s, T):
     return None
 
 
+@lru_cache(maxsize=8)
+def _quadrature_rule(points):
+    """Reuse the nodes and weights, which depend only on quadrature order."""
+    return np.polynomial.legendre.leggauss(points)
+
+
 def interval_bases(n, T, points=16):
     """Single-step interval averages over [n, n + 1] of cos phi |cos phi| and sin phi
     (phi in radians = 2 pi s / T - pi / 2), Gauss-Legendre quadrature."""
-    x, w = np.polynomial.legendre.leggauss(points)
+    x, w = _quadrature_rule(points)
     s = int(n) % T + 0.5 * (x + 1.0)  # periodic: identical bases in every cycle
     phi = 2.0 * np.pi * s / T - 0.5 * np.pi
     c = np.cos(phi)

research/karman_reproduce/run.py

@@ -0,0 +1,116 @@
+"""Replay the frozen Legacy karman_re100 PPO policy on current CelerisLab.
+
+Research only; nothing in src/ is changed. Plant: Legacy Karman cloak scene
+(1280x512, D=20, parabolic inlet, mean U0 .01, nu .004, Re_D 50; upstream
+cylinder r=20 at x=200; pinball front (600, 255.5), rear (626, 255.5+-15);
+probes r=5 at x=800, y=255.5+40/0/-40). Interface contract = DynamisLab
+legacy_test (the executor behind the Legacy rerun), not f2f88c phase3:
+
+- command per interval: A = (8 a + b) U0, b = (0, -4, 4), a = deterministic
+  policy output in [-1, 1]^3; wall speed command, omega = +A / R (both
+  solvers impose u_w = omega x r, CCW positive);
+- smoothing per LATTICE step, A_eff <- .9 A_eff + .1 A, with A_eff = 0 at the
+  start of every 800-step interval (legacy_test.runtime.run_historical_interval);
+- observation: interval means of the three pinball forces (front, top, bottom)
+  and the three probe (ux, uy) sums, normalized by the frozen norm.json with
+  legacy_test.runtime.policy_observation; the first observation is zero.
+
+Output (--out, outside the checkout): signals.npz with sensors (N,6), forces
+(N,6), actions (N,3) in the f2f88c phase3 layout, so its comparator applies.
+
+    PYTHONPATH=src python research/karman_reproduce/run.py --dynamis DYN --model .../d1a3o12_re100.zip \\
+        --norm .../karman_re100/norm.json --config .../config_lbm_pinball_legacy_compat.json \\
+        --out /tmp/new-dir [--device 1]
+"""
+from __future__ import annotations
+
+import argparse
+import json
+import os
+import subprocess
+import sys
+
+import numpy as np
+
+L0, U0, NX, NY = 20.0, 0.01, 1280, 512
+R, YC = L0 / 2, (NY - 1) / 2
+SI, INTERVALS, WARMUP, FIFO_LEN = 800, 200, int(4 * NX / U0), 150
+BIAS = np.array([0.0, -4.0, 4.0])
+
+
+def git_head(path):
+    return subprocess.run(["git", "-C", path, "rev-parse", "HEAD"], capture_output=True,
+                          text=True).stdout.strip()
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--dynamis", required=True, help="DynamisLab checkout (imported read-only)")
+    ap.add_argument("--model", required=True, help="frozen d1a3o12_re100.zip")
+    ap.add_argument("--norm", required=True, help="frozen karman_re100 norm.json")
+    ap.add_argument("--config", required=True)
+    ap.add_argument("--out", required=True)
+    ap.add_argument("--device", type=int, default=0)
+    ap.add_argument("--intervals", type=int, default=INTERVALS)
+    args = ap.parse_args()
+    here = os.path.realpath(os.path.join(os.path.dirname(__file__), "..", ".."))
+    if os.path.realpath(args.out).startswith(here):
+        raise SystemExit("--out must be outside the checkout")
+    os.makedirs(args.out)  # fresh directory only; never overwrite results
+
+    sys.dont_write_bytecode = True  # import only: write nothing into the DynamisLab tree
+    sys.path.insert(0, os.path.join(args.dynamis, "src"))
+    from stable_baselines3 import PPO
+    from drl_pinball.legacy_test.core.model_loader import _dummy_env
+    from drl_pinball.legacy_test.runtime import load_policy_norm, policy_observation
+    model = PPO.load(args.model, env=_dummy_env(12), device="cpu")
+    norm = load_policy_norm(args.norm)
+
+    from CelerisLab.simulation import Simulation
+    sim = Simulation(args.config, device_id=args.device)
+    try:
+        sim.add_body("circle", center=(10 * L0, YC, 0.0), radius=L0)
+        probes = [sim.add_body("sensor", center=(40 * L0, YC + dy, 0.0), radius=5.0)
+                  for dy in (40.0, 0.0, -40.0)]
+        sim.initialize()
+        # Legacy build_karman_cloak order: stabilize disturbance + probes, then the
+        # 150-interval zero-action target sampling (kept as flow time), then the pinball.
+        sim.run(WARMUP + FIFO_LEN * SI)
+        n0 = sim.bodies.count
+        for c in ((30 * L0, YC, 0.0), (626.0, YC + 15, 0.0), (626.0, YC - 15, 0.0)):
+            sim.add_body("circle", center=c, radius=R)  # pending (negative) id until sync
+        sim.sync_bodies()
+        pinball = list(range(n0, n0 + 3))  # committed ids, as in f2f88c phase3
+        sim.run(WARMUP)
+        cells = float(len(sim.bodies.get(probes[0]).get_sensor_list(NX, NY)[0]))  # probe sum = mean * cells
+
+        obs = np.zeros(12, np.float32)
+        sensors, forces, actions = (np.zeros((args.intervals, n), np.float32) for n in (6, 6, 3))
+        for n in range(args.intervals):
+            a = np.asarray(model.predict(obs, deterministic=True)[0], np.float32).reshape(3)
+            target, eff = (8.0 * a + BIAS) * U0, np.zeros(3)
+            for k in range(SI):
+                eff = 0.9 * eff + 0.1 * target
+                for pid, w in zip(pinball, eff / R):
+                    sim.set_body(pid, omega=float(w))
+                sim.run(1, zero_obs=(k == 0), sync_obs=(k == SI - 1))
+            s = [v * cells for pid in probes for v in sim.read_sensor(pid)[:2]]
+            f = [v for pid in pinball for v in sim.read_force(pid)[:2]]
+            raw = np.asarray(s + f, np.float32)
+            if not np.all(np.isfinite(raw)):
+                raise RuntimeError(f"non-finite observation at interval {n}")
+            sensors[n], forces[n], actions[n] = raw[:6], raw[6:], a
+            obs = policy_observation(raw, norm)
+        np.savez_compressed(os.path.join(args.out, "signals.npz"), sensors=sensors,
+                            forces=forces, actions=actions)
+        with open(os.path.join(args.out, "case.json"), "w") as fh:
+            json.dump(dict(celeris=git_head(here), dynamis=git_head(args.dynamis),
+                           config=os.path.realpath(args.config), model=os.path.realpath(args.model),
+                           norm=os.path.realpath(args.norm), warmup=WARMUP, target_sampling=FIFO_LEN * SI, interval=SI, intervals=args.intervals,
+                           probe_cells=cells), fh, indent=1)
+    finally:
+        sim.close()
+
+
+if __name__ == "__main__":
+    main()

research/viv_re150/run.py

@@ -1,8 +1,9 @@
 """One-shot predeclared VIV run (docs/validation/viv_re150_predeclare.md).
 
-Re 150, U* = 5, m/(rho D^2) = 2, zeta = 0, D = 40 on a 75 D x 60 D native
-channel (uniform inlet, NEQ outlet, free-slip y walls), cylinder at
-(30 D, 30 D). Structural RHS = CME + Wen (diagnostics: raw CME, node terms).
+Re 150, U* = 5, m/(rho D^2) = 2, zeta = 0, D = 40 on a nominal 75 D x 60 D
+lattice (uniform inlet, NEQ outlet, free-slip y walls). The cylinder is at
+(30 D, (ny - 1)/2); y-wall planes are .5 and ny - 1.5 in lattice coordinates.
+Structural RHS = CME + Wen (diagnostics: raw CME, node terms).
 Research driver only: Simulation API, no production change.
 
     PYTHONPATH=src python research/viv_re150/run.py --out /path/outside/checkout

tests/unit/test_duetsch_protocol.py

@@ -82,6 +82,35 @@ def test_interval_bases_and_morison_fit_recover_known_coefficients():
     assert p.morison_fit(f, b[:, 0], b[:, 1]) == pytest.approx((2.09, 1.45), abs=1e-9)
 
 
+def test_quadrature_rule_reused_across_steps_and_periods(monkeypatch):
+    original = np.polynomial.legendre.leggauss
+    orders = []
+
+    def counted(points):
+        orders.append(points)
+        return original(points)
+
+    p._quadrature_rule.cache_clear()
+    monkeypatch.setattr(np.polynomial.legendre, "leggauss", counted)
+    try:
+        for points in (16, 8, 16):
+            for period in (T, 2 * T):
+                for n in (0, 123, period + 123):
+                    p.interval_bases(n, period, points)
+        assert orders == [16, 8]
+    finally:
+        p._quadrature_rule.cache_clear()
+
+
+def test_interval_sine_basis_matches_exact_integral():
+    for period in (T, 2 * T):
+        half_width = math.pi / period
+        for n in (0, period // 4, period // 2, period - 1, 7 * period + 123):
+            midpoint = 2 * math.pi * ((n % period) + .5) / period - math.pi / 2
+            expected = math.sin(midpoint) * math.sin(half_width) / half_width
+            assert p.interval_bases(n, period)[1] == pytest.approx(expected, rel=0, abs=2e-15)
+
+
 def test_geometry_repeats_exactly_across_cycles():
     A = p.amplitude(20)
     for m in (0, 1, T // 4, T // 2 - 1, 7 * T // 12):