research: reuse quadrature and clarify snapshot and VIV coordinates

Index values are unchanged (1, 3, 7, 10, 13); the run recorded at fb507aa1
used the same columns. Research-only.

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
+- Nominal 75 D × 60 D lattice allocation (3000 × 2400), based on the
-  (30 D, 30 D): 30 D from the inlet, 45 D from the outlet, centred in y
+  Kolahdouz §4.1.1 layout. Coordinate clarification of the existing run:
-  (Kolahdouz §4.1.1). Qin's own box ([0, 40 D] × [0, 20 D], centre (10 D, 10 D))
+  the driver uses `(30 D, (ny - 1)/2)` = `(1200, 1199.5)` in lattice
-  is the target paper's, not ours.
+  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/viv_re150/README.md

@@ -0,0 +1,23 @@
+# Research VIV driver — coupling contract (not production)
+
+`run.py` couples a 1-DOF undamped spring–mass (`viv.py`) to the production
+moving-circle API. Everything specific to fluid–structure coupling lives here,
+not in `src/`.
+
+- Per step: known (y_n, v_n) and the last completed force F_last → exact
+  harmonic candidate (y_{n+1}, v_{n+1}) → `set_body(offset=(0, y_{n+1}),
+  vy=y_{n+1} − y_n)` (interval wall velocity, lattice dt = 1) → `run(1)` →
+  commit y, v and F_last only after the step succeeded. One step of explicit
+  force lag (weak coupling); no sub-iteration.
+- Structural RHS = CME + Wen (y), read from `moving_diagnostics()` by the
+  names in `MovingTables.DIAG`; raw CME, the `obs` force and node terms are
+  logged as diagnostics. This is a stage choice, not a validated force law.
+- Production guarantees used: fail-closed target commit (synchronized step,
+  `no_donor == 0`), input rejection before mutation, a failed device step
+  marks the Simulation failed.
+- Not provided, by design: device-resident structure state or batched
+  (multi-step) coupling, strong/sub-iterated coupling (would need state
+  rollback), more than one DOF. These are future work and would be a separate
+  execution path, not an abstraction over this driver.
+- Recorded result and limits: `docs/validation/moving_bodies_stage1.md`;
+  predeclaration: `docs/validation/viv_re150_predeclare.md`.

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
+Re 150, U* = 5, m/(rho D^2) = 2, zeta = 0, D = 40 on a nominal 75 D x 60 D
-channel (uniform inlet, NEQ outlet, free-slip y walls), cylinder at
+lattice (uniform inlet, NEQ outlet, free-slip y walls). The cylinder is at
-(30 D, 30 D). Structural RHS = CME + Wen (diagnostics: raw CME, node terms).
+(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
@@ -65,7 +66,8 @@ def main():
         step_cap = args.steps
     centre = (30.0 * Dn, 0.5 * (ny - 1))
 
-    from CelerisLab.body.motion import Translation
+    from CelerisLab.body.motion import MovingTables, Translation
+    col = {k: MovingTables.DIAG.index(k) for k in ("cme_y", "wen_y", "cover_jy", "uncover_jy", "no_donor")}
     from CelerisLab.simulation import Simulation
     with tempfile.NamedTemporaryFile("w", suffix=".json", delete=False) as fh:
         json.dump(cfg, fh)
@@ -95,9 +97,9 @@ def main():
             sim.run(1)                             # device failures raise and mark the Simulation failed
             d = sim.moving_diagnostics()[0]
             f_obs = sim.read_force(body)[1]
-            f_rhs = float(d[1] + d[3])             # CME + Wen (y)
+            f_rhs = float(d[col["cme_y"]] + d[col["wen_y"]])   # CME + Wen (y)
-            row = (n, y1, v1, f_rhs, float(d[1]), float(d[3]), float(f_obs),
+            row = (n, y1, v1, f_rhs, float(d[col["cme_y"]]), float(d[col["wen_y"]]), float(f_obs),
-                   float(d[7] - d[10]), float(d[13]))
+                   float(d[col["cover_jy"]] - d[col["uncover_jy"]]), float(d[col["no_donor"]]))
             if not viv.finite(row):                # every recorded output, not only the RHS
                 status = "fail_nonfinite_output"
                 break

src/CelerisLab/body/motion.py

@@ -1,10 +1,10 @@
 # CelerisLab/body/motion.py
-"""Host state of a translating circle (dense double buffer, D2Q9 only).
+"""Translation targets and device tables (dense double buffer, D2Q9 only).
 
 Contract: docs/moving_bodies.md, sections 4a, 5 and 7. The applied geometry
 (last successfully committed offset) and an optional pending target are two
-separate facts; centres are always ``base + offset``. Nothing here touches
+separate facts; centres are always ``base + offset``. MovingCircle owns host
-the device: the caller supplies upload/step/synchronize callables.
+targets, MovingTables owns device records, and Simulation commits completed steps.
 """
 from __future__ import annotations
 
@@ -18,10 +18,6 @@ R_MIN = 4.0        # smallest radius covered by the donor enumeration
 
 Vec = Tuple[float, float]
 
-# D2Q9 non-rest directions, paired order of lbm/descriptors.
-_CX = (1, -1, 0, 0, 1, -1, 1, -1)
-_CY = (0, 0, 1, -1, 1, -1, -1, 1)
-
 
 def _vec(value, name: str) -> Vec:
     x, y = (float(v) for v in value)
@@ -98,28 +94,6 @@ class MovingCircle:
             offset, velocity = self.applied, (0.0, 0.0)
         return self.centre(self.applied), self.centre(offset), velocity
 
-    def advance(self, upload: Callable, step: Callable, synchronize: Callable) -> None:
-        """Single-circle transaction helper (unit tests): upload, step, sync, commit.
-
-        Not the production path: ``Simulation._run_moving`` batches all circles,
-        adds the ``no_donor`` gate and uses ``commit``/``fail`` directly.
-
-        Upload the record, run one step and commit only after a successful sync.
-
-        Any exception after the upload is attempted leaves device geometry
-        unknown: the circle is marked failed, never reports applied state again,
-        and the owner's ``on_failure`` hook is called.
-        """
-        record = self.record()
-        try:
-            upload(record)
-            step()
-            synchronize()
-        except BaseException as exc:
-            self.fail(exc)
-            raise
-        self.commit()
-
     def commit(self) -> None:
         """Applied := pending offset (caller has synchronized a successful step)."""
         if self.pending is not None:
@@ -143,43 +117,6 @@ class MovingCircle:
             raise RuntimeError("moving circle state is unknown after a failed step")
 
 
-def _inside(x: int, y: int, c: Vec, r: float) -> bool:
-    return (x - c[0]) ** 2 + (y - c[1]) ** 2 < r * r
-
-
-def empty_donor_cells(c_old: Vec, c_new: Vec, radius: float):
-    """Fresh cells (inside c_old, outside c_new) whose donor set V(x) is empty."""
-    r = radius
-    xs = range(math.floor(min(c_old[0], c_new[0]) - r) - 1, math.ceil(max(c_old[0], c_new[0]) + r) + 2)
-    ys = range(math.floor(min(c_old[1], c_new[1]) - r) - 1, math.ceil(max(c_old[1], c_new[1]) + r) + 2)
-    empty = []
-    for x in xs:
-        for y in ys:
-            if not _inside(x, y, c_old, r) or _inside(x, y, c_new, r):
-                continue
-            if not any(not _inside(x + cx, y + cy, c_old, r) and not _inside(x + cx, y + cy, c_new, r)
-                       for cx, cy in zip(_CX, _CY)):
-                empty.append((x, y))
-    return empty
-
-
-def enumerate_donor_coverage(radius: float, positions: int = 8, directions: int = 16) -> int:
-    """Count empty-donor fresh cells over sub-cell centres and moves of MAX_STEP.
-
-    Evidence for the declared bounds (radius >= R_MIN, |dc|_2 <= MAX_STEP),
-    not a proof; the device keeps the ``no_donor`` counter.
-    """
-    count = 0
-    for i in range(positions):
-        for j in range(positions):
-            c_old = (i / positions, j / positions)
-            for k in range(directions):
-                a = 2.0 * math.pi * k / directions
-                c_new = (c_old[0] + MAX_STEP * math.cos(a), c_old[1] + MAX_STEP * math.sin(a))
-                count += len(empty_donor_cells(c_old, c_new, radius))
-    return count
-
-
 class MovingTables:
     """Device tables of the translating circles for ``MovingCircleCells``/``Links``.
 

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

tests/unit/test_moving_bodies.py

@@ -1,14 +1,13 @@
-"""Translating-circle host contract, S1 (docs/moving_bodies.md 4a, 5, 7).
+"""Translating-circle host contract (docs/moving_bodies.md 4a, 5, 7).
 
-Host doubles only: no kernels exist yet, so every positive-step run with a
+Real Simulation entrypoints with device doubles; no CUDA context or GPU work.
-moving body must be rejected before any device write.
 """
+import math
 from types import SimpleNamespace
 
 import pytest
 
-from CelerisLab.body.motion import (MAX_STEP, R_MIN, MovingCircle, Translation,
+from CelerisLab.body.motion import MAX_STEP, R_MIN, MovingCircle, Translation
-                                    enumerate_donor_coverage)
 
 BAD = (float("nan"), float("inf"))
 
@@ -46,46 +45,96 @@ def test_targets_validate_against_applied_and_last_wins():
     assert c.applied == (0., 0.)
 
 
-def test_record_commit_and_hold():
+def test_record_commit_and_hold(monkeypatch):
-    c = circle()
+    sim, writes, c = routed_advancing_simulation(monkeypatch)
-    assert c.record() == ((50., 40.), (50., 40.), (0., 0.))  # hold: no translation velocity
+    assert c.record() == ((50., 40.), (50., 40.), (0., 0.))
-    c.set_target((0., .25), (0., .05))
+    sim.set_body(0, offset=(0., .25), vy=.05)
     assert c.record() == ((50., 40.), (50., 40.25), (0., .05))
-    uploads = []
+    seen_at_sync = []
-    c.advance(uploads.append, lambda: None, lambda: None)
+    sim.run(1, stream=SimpleNamespace(synchronize=lambda: seen_at_sync.append(c.applied)))
-    assert uploads == [((50., 40.), (50., 40.25), (0., .05))]
+    assert seen_at_sync == [(0., 0.), (0., .25)]  # commit between step sync and hold sync
+    uploads = [w[1][0] for w in writes if isinstance(w, tuple) and w[0] == "upload"]
+    assert uploads == [((50., 40.), (50., 40.25), (0., .05)),
+                       ((50., 40.25), (50., 40.25), (0., 0.))]
     assert c.applied == (0., .25) and c.pending is None  # offset-only commit
-    c.advance(uploads.append, lambda: None, lambda: None)  # no new target: hold step
+    sim.run(1, stream=STREAM)  # no new target: hold step
+    uploads = [w[1][0] for w in writes if isinstance(w, tuple) and w[0] == "upload"]
     assert uploads[-1] == ((50., 40.25), (50., 40.25), (0., 0.)) and c.applied == (0., .25)
-    c.set_target((0., 0.), (0., -.05))  # return to zero offset is a real target
+    sim.set_body(0, offset=(0., 0.), vy=-.05)  # zero offset is a real target
     assert c.pending == ((0., 0.), (0., -.05))
-    c.advance(uploads.append, lambda: None, lambda: None)
+    sim.run(1, stream=STREAM)
     assert c.applied == (0., 0.) and c.pending is None
 
 
 @pytest.mark.parametrize("stage", ("upload", "step", "synchronize"))
-def test_failure_after_upload_attempt_is_fail_closed(stage):
+def test_failure_after_upload_attempt_is_fail_closed(monkeypatch, stage):
-    c = circle()
+    sim, writes, c = routed_advancing_simulation(monkeypatch, fail="upload" if stage == "upload" else None)
-    c.set_target((0., .25), (0., .05))
+    sim.set_body(0, offset=(0., .25), vy=.05)
 
-    def call(name):
+    def fail(*a, **k):
-        def run(*_):
+        raise RuntimeError(f"{stage} failed")
-            if name == stage:
+    if stage == "step":
-                raise RuntimeError(f"{name} failed")
+        sim.stepper._step = fail
-        return run
+    stream = SimpleNamespace(synchronize=fail) if stage == "synchronize" else STREAM
     with pytest.raises(RuntimeError, match=f"{stage} failed"):
-        c.advance(call("upload"), call("step"), call("synchronize"))
+        sim.run(1, stream=stream)
-    assert c.failed and c.applied == (0., 0.)
+    assert c.failed and c.applied == (0., 0.) and c.pending == ((0., .25), (0., .05))
-    for later in (c.record, lambda: c.set_target((0., .1), (0., 0.)),
+    assert sim._failed is not None and sim.stepper.blocked == sim._failed
-                  lambda: c.advance(*(lambda *_: None,) * 3)):
+    for later in (c.record, lambda: sim.set_body(0, offset=(0., .1), vy=0.)):
         with pytest.raises(RuntimeError, match="unknown"):
             later()
+    with pytest.raises(RuntimeError, match="Simulation failed"):
+        sim.run(1, stream=STREAM)
+    with pytest.raises(RuntimeError, match="moving circle step failed"):
+        sim.stepper.step(1, action_gpu=None, obs_gpu=None)
 
 
 def test_band_covers_offset_box_with_margin():
     assert circle().band() == (50 - 2 - 10 - 2, 50 + 2 + 10 + 2, 40 - 3 - 10 - 2, 40 + 3 + 10 + 2)
 
 
+# D2Q9 non-rest directions, paired order of lbm/descriptors.
+_CX = (1, -1, 0, 0, 1, -1, 1, -1)
+_CY = (0, 0, 1, -1, 1, -1, -1, 1)
+
+
+def _inside(x: int, y: int, c, r: float) -> bool:
+    return (x - c[0]) ** 2 + (y - c[1]) ** 2 < r * r
+
+
+def empty_donor_cells(c_old, c_new, radius: float):
+    """Fresh cells (inside c_old, outside c_new) whose donor set V(x) is empty."""
+    r = radius
+    xs = range(math.floor(min(c_old[0], c_new[0]) - r) - 1, math.ceil(max(c_old[0], c_new[0]) + r) + 2)
+    ys = range(math.floor(min(c_old[1], c_new[1]) - r) - 1, math.ceil(max(c_old[1], c_new[1]) + r) + 2)
+    empty = []
+    for x in xs:
+        for y in ys:
+            if not _inside(x, y, c_old, r) or _inside(x, y, c_new, r):
+                continue
+            if not any(not _inside(x + cx, y + cy, c_old, r) and not _inside(x + cx, y + cy, c_new, r)
+                       for cx, cy in zip(_CX, _CY)):
+                empty.append((x, y))
+    return empty
+
+
+def enumerate_donor_coverage(radius: float, positions: int = 8, directions: int = 16) -> int:
+    """Count empty-donor fresh cells over sub-cell centres and moves of MAX_STEP.
+
+    Evidence for the declared bounds (radius >= R_MIN, |dc|_2 <= MAX_STEP),
+    not a proof; the device keeps the ``no_donor`` counter.
+    """
+    count = 0
+    for i in range(positions):
+        for j in range(positions):
+            c_old = (i / positions, j / positions)
+            for k in range(directions):
+                a = 2.0 * math.pi * k / directions
+                c_new = (c_old[0] + MAX_STEP * math.cos(a), c_old[1] + MAX_STEP * math.sin(a))
+                count += len(empty_donor_cells(c_old, c_new, radius))
+    return count
+
+
 @pytest.mark.parametrize("radius", (R_MIN, 4.3, 7.5, 12.0))
 def test_refill_donor_sets_nonempty_within_bounds(radius):
     # Evidence for R >= R_MIN and |dc|_2 <= MAX_STEP on a fixed sample, not a proof.
@@ -278,6 +327,26 @@ def advancing_simulation(monkeypatch, **tables):
     return sim, writes, sim.bodies.objects[0].motion
 
 
+def routed_advancing_simulation(monkeypatch, **tables):
+    """Real body registration, target routing, run and owner hook; double device I/O."""
+    sim = routed_simulation(monkeypatch)
+    bid = sim.add_body("circle", (50., 40.), 10., translation=Translation((-2., 2., -3., 3.), .1))
+    sim._initialized = True
+    sim._assert_runtime_contracts = lambda: None
+    writes = []
+    monkeypatch.setattr(sim.bodies, "_upload_action_async", lambda s: writes.append("action"))
+    monkeypatch.setattr(sim.bodies, "zero_obs_async", lambda s: writes.append("obs"))
+    monkeypatch.setattr(sim.bodies, "download_obs_full_async", lambda s: writes.append("download"))
+    sim.stepper = SimpleNamespace(
+        moving=Tables(writes, **tables), blocked=None,
+        _step=lambda n, advance=False, **k: writes.append(("step", n, advance)))
+    from CelerisLab.lbm.stepper import LBMStepper
+    sim.stepper.step = LBMStepper.step.__get__(sim.stepper)
+    sim._moving_tables = sim.stepper.moving
+    sim._attach_moving()
+    return sim, writes, sim.bodies.get(bid).motion
+
+
 def test_run_one_advances_and_commits_after_sync(monkeypatch):
     sim, writes, motion = advancing_simulation(monkeypatch)
     motion.set_target((0., .25), (0., .05))
@@ -545,23 +614,21 @@ def test_moving_circle_excluded_from_static_links():
 
 
 def test_failed_moving_step_marks_simulation_failed(monkeypatch):
-    sim = routed_simulation(monkeypatch)
+    sim, writes, motion = routed_advancing_simulation(monkeypatch)
     sim.ctx = ClosableOwner()
-    bid = sim.add_body("circle", (50., 40.), 10., translation=Translation((0., 0., -1., 1.), .1))
+    sim.set_body(0, offset=(0., .2), vy=.05)
-    motion = sim.bodies.get(bid).motion
-    motion.set_target((0., .2), (0., .05))
 
     def fail(*_):
         raise RuntimeError("sync failed")
     with pytest.raises(RuntimeError, match="sync failed"):
-        motion.advance(lambda *_: None, lambda: None, fail)
+        sim.run(1, stream=SimpleNamespace(synchronize=fail))
     assert sim._failed is not None and motion.applied == (0., 0.)
     for call in (lambda: sim.run(0), lambda: sim.run(1), sim.step, sim.sync_bodies,
                  sim.save_checkpoint, lambda: sim.load_checkpoint("x.h5")):
         with pytest.raises(RuntimeError, match="Simulation failed"):
             call()
     with pytest.raises(RuntimeError, match="unknown"):  # no retry on the failed circle
-        motion.advance(lambda *_: None, lambda: None, lambda: None)
+        sim.set_body(0, offset=(0., .1), vy=.05)
     with sim.ctx.activate():  # a refused close keeps the failed state and resources
         with pytest.raises(RuntimeError, match="active CUDA"):
             sim.close()