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/lifecycle.md

@@ -6,8 +6,8 @@ it is; it adds no framework. Line references are to candidate 0f8fdbc3.
 | Class | Contents | How it changes | Owner rule |
 |---|---|---|---|
 | 1. Module and its globals | Compile-time (visible to the compiler, fixed per module): `config.h` macros incl. `NX`/`NY`, `COLLISION_MODEL`, `BLOCK_ROOT_BC`, `WALL_*`. Runtime per-module globals (not visible to the compiler; written by host code): `__constant__` `d_params`, `d_w`, `d_cx`, `d_cy`, `d_block_channel`, `d_block_wall`; plus bound kernel handles | Globals are written once by the owner right after `load_module`. Uniform grid: `Simulation._recompile` (construction, body-topology sync) loads a new module and re-uploads; `update_runtime_params` rewrites `d_params` (omega, force, n_objects) between runs. Steppers: never rewritten after construction, so they are immutable for the instance's lifetime although they are runtime values. The blocked slot/reservoir root reads its slot from the `d_block_wall` constant (the former dense route used `WALL_*` macros; both were bit-identical in the parity gates). | Every owner loads its own module (`compiler.load_module` = `module_from_file`, no cache). A loaded module is never shared between live owners: globals are per module, so sharing would share mutable `d_params`. |
-| 2. Geometry, layout, tables | flags, curved-link tables, block layouts and halo pairs, transfer (explode/coalesce) tables, block kinds | Immutable for the life of an instance. Uniform grid: rebuilt by body sync. `FixedPatchStepper`: a geometry change is a new instance (`remap_finest`), committed only after a successful transfer. | No in-place table edits. |
-| 3. Controls and dependent boundary columns | Jet amplitude -> per-link `surface_velocity`; body actions (uniform grid) | Between steps, without a table or module rebuild (`set_jet_amplitude`). | Input (index, level, amplitude) is validated before any device write; invalid input changes nothing and the instance stays usable. The host copy commits only after the upload returned. Any error after validation leaves the device column undefined and fails the instance: `FixedPatchStepper` status `failed`, `Simulation._failed`; every later device-touching public method (`Simulation`: `stream`, `set_jet_amplitude`, `initialize`, `sync_bodies`, `run`/`step`, `read_surface_exchange`, `read_body(s)`, `get_macroscopic`, `get_ddf`, `update_runtime_params`, `snapshot`, `restore`, save/load checkpoint) raises; cached host readers and host-side staging (e.g. `get_flags`, `add_body`/`remove_body` before a sync) do not touch the device and stay callable. `close()` still releases resources and stays retryable. |
+| 2. Geometry, layout, tables | flags, curved-link tables, block layouts and halo pairs, transfer (explode/coalesce) tables, block kinds | Immutable for the life of an instance, except translating circles on the uniform grid: their swept-band flags and cut links are rewritten on the device every advancing step (`MovingCircleCells`/`MovingCircleLinks`, docs/moving_bodies.md); the host flag mirror is then stale for those cells. Uniform grid otherwise: rebuilt by body sync. `FixedPatchStepper`: a geometry change is a new instance (`remap_finest`), committed only after a successful transfer. | No in-place table edits. |
+| 3. Controls and dependent boundary columns | Jet amplitude -> per-link `surface_velocity`; body actions (uniform grid) | Between steps, without a table or module rebuild (`set_jet_amplitude`). | Input (index, level, amplitude) is validated before any device write; invalid input changes nothing and the instance stays usable. The host copy commits only after the upload returned. Any error after validation leaves the device column undefined and fails the instance: `FixedPatchStepper` status `failed`, `Simulation._failed`; every later device-touching public method (`Simulation`: `stream`, `set_jet_amplitude`, `initialize`, `sync_bodies`, `run`/`step`, `read_surface_exchange`, `read_body(s)`, `get_macroscopic`, `get_ddf`, `update_runtime_params`, `snapshot`, `restore`, save/load checkpoint) raises; cached host readers and host-side staging (e.g. `get_flags` without translating circles, `add_body`/`remove_body` before a sync) do not touch the device and stay callable; with translating circles `get_flags` reads the device flags back and refuses after a failure. Translating-circle targets follow the same rule: `set_body` validates before mutation, the applied offset commits only after a successful synchronized step, and a failed moving step fails the instance (docs/moving_bodies.md). `close()` still releases resources and stays retryable. |
 | 4. Per-step state | population double buffers and their parity, raw-exchange ledgers, observations | Device-only during `step(n)` / `run(n)`; buffer swaps are host-handle swaps. | Valid on the host only after the single wait at the end of the call. |
 
 Storage routes are separate concerns:

docs/moving_bodies.md

@@ -1,17 +1,65 @@
-# Moving circles on the uniform dense grid — S0 DRAFT contract
+# Moving circles on the uniform dense grid
 
-Status: **DRAFT, not frozen** for the force RHS, kernels and physical
-acceptance. Host-only S1 (state, validation, guards; no kernels) is on branch
-`claude/moving-bodies-s1`; positive-step runs with a moving body are rejected
-until S2. When S2 lifts that guard, `run(n > 1)` with a pending target must
-raise `ValueError` before mutation (4a), and `run(0)` must not upload a motion
-record or consume the pending target. Base: main
-`b2034f87e03656c5f67ab88562a70d6c338b9918`. Scope approved by Frank (via TRON,
-thread C0C78UFKGRM/1793335557.460509): uniform-grid prescribed translation and
-single-cylinder free transverse VIV. Out of scope: refinement/blocked path,
-pinball assembly, EsoPull, shifted storage, D3Q19, DynamisLab, manuscript.
+## Status and user guide (current)
 
-Items marked **OPEN** are undecided; items marked **CANDIDATE** are proposals.
+Development-stage feature, integrated from the reviewed slices S1–S3a
+(production code identical to `653fd997`). Physical evidence and its limits:
+`docs/validation/moving_bodies_stage1.md`. The force used as a structural RHS
+is a declared stage choice, not an established result.
+
+*What it does.* Translating circles (with optional spin) on the dense
+double-buffer D2Q9 path: per-step cut links from an FP64 geometry record
+(M2 `MovingCircleLinks`), reclassification and refill of covered/uncovered
+cells (M1 `MovingCircleCells`), moving-wall Bouzidi with wall velocity
+`u + omega x r`, and separate per-body diagnostics.
+
+*API.*
+- `Simulation.add_body("circle", center=(x, y), radius=r,
+  translation=Translation(offset_box, max_speed))` before `initialize()`;
+  `offset_box = (x_min, x_max, y_min, y_max)` in lattice cells around the
+  base centre (must contain 0), `0 < max_speed ≤ .5`.
+- `set_body(id, offset=(dx, dy), vx=..., vy=..., omega=...)` sets a pending
+  target (validated against the applied offset: box, Euclidean step ≤ .5
+  cell, speed, no sliding hold); `vx`/`vy` are the interval wall velocity.
+- `run(1)` advances a pending target and commits it only after the step
+  synchronizes with `no_donor == 0`; `run(n > 1)` with a pending target
+  raises; without a target `run(n)` holds the geometry (spin allowed).
+  Direct `LBMStepper.step` refuses while a target is pending.
+- `moving_diagnostics()` returns per body the columns of
+  `MovingTables.DIAG` (CME, Wen term, link mass, cover/uncover momentum and
+  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 /
+blocked steppers, Jets or surface actuation on a moving circle,
+`sync_bodies` and checkpoints while a moving circle exists, swept bands that
+touch the domain edge, other bodies, sensors or force regions.
+
+*Failure lifecycle.* Device-side failures (upload, step, sync, no_donor)
+mark the Simulation failed; later run/step/sync/checkpoint calls and direct
+stepper calls refuse until `close()`.
+
+*Known limitations.* No mass correction at curved walls (static or moving);
+cover/uncover and Wen terms are diagnostics, not a selected force law; the
+host loop costs a fixed per-step synchronization (≈ 1 ms/step measured on a
+7.2 M-cell grid, run(1) class); outer boundaries are the single dense
+channel package (see `configs/CONFIG.md`).
+
+## Design record (staged contract, kept as written)
+
+The sections below are the contract the slices were reviewed against
+(S0 → S3a). They are kept unchanged for provenance; where they say "draft",
+"proposed" or "until S2/S3", the status above supersedes them.
+
+Original scope (Frank via TRON, thread C0C78UFKGRM/1793335557.460509):
+uniform-grid prescribed translation and single-cylinder free transverse VIV;
+base main `b2034f87e03656c5f67ab88562a70d6c338b9918`.
 
 ## 1. Working standards
 

docs/supported_settings.md

@@ -82,6 +82,13 @@ Tests: `tests/unit/test_jets.py`, `test_jets_review.py`, `test_shifted_jets.py`,
 - DFG 2D-1 Re 20 (V4, SRT, t* 60, pressure-NEQ outlet and local Zou-He inlet): refined and uniform
   runs were NOT_STEADY under the frozen steady rule, so no interval verdict was given.
 
+## Moving circles (dense `Simulation`, uniform grid only)
+
+Translating circles with optional spin are supported on the dense double
+buffer D2Q9 FP32 path only (not on `FixedPatchStepper` or refined levels).
+API, supported/rejected modes and limits: `docs/moving_bodies.md`; stage-1
+evidence: `docs/validation/moving_bodies_stage1.md`.
+
 ## Not supported in this release
 
 3D, dynamic/automatic AMR (tagging, time-varying hierarchies beyond `remap_finest`),

docs/validation/moving_bodies_stage1.md

@@ -0,0 +1,46 @@
+# Moving circles — stage-1 evidence (development level)
+
+Owner criterion (Yanqi via TRON): prescribed motion and one free VIV case both
+"reasonably close" to a predeclared reference = development-stage success.
+This is practical closeness, not a formal validation claim. Every number
+below is from Astra's receipts as posted in Slack C0C78UFKGRM thread
+1793335557.460509; raw outputs stay on Horsea under `/tmp/celeris-results-*`
+and the Mac evidence directory (not in git).
+
+## Code
+
+- Production: `src/` identical to `653fd997` (slices S1 `983020b6`,
+  S2a `b920a303`, S3a `653fd997`). CPU suites and GPU hold/crossing gates
+  passed on those slices (G0a/G0b population-bitwise; X1 crossing ledger).
+- Research drivers: `research/duetsch_re100_kc5/` (`dd37883d`, capture
+  `a0aea3b9`) and `research/viv_re150/` (`fb507aa1`), run as committed.
+
+## Prescribed in-line oscillation (Dütsch et al. 1998, Re 100, KC 5)
+
+- Protocol: `docs/validation/duetsch_re100_kc5_protocol.md`.
+- Frozen settling gate (.002 F_ref waveform): **FAILED** and kept on record:
+  D40 domain stage (20/40/60 D, 11.25 T) and the 60 D duration diagnostic
+  (31.25 T, 0/28 pairs; last waveform drift .00613 F_ref).
+- Unscored late-cycle CME fits c_d ≈ 2.02–2.05, c_i ≈ 1.46–1.48 vs the
+  paper's computed (2.09, 1.45); fields clean and symmetric (owner view).
+  The owner accepts ~0.61 % cycle drift at D40 for this stage; the gate is
+  not re-scored.
+- Fig. 8 LDA profiles: qualitative only (no digitized ordinates).
+
+## Free transverse VIV (Kolahdouz §4.1.1 setup, Qin §4.5 target)
+
+- Predeclaration: `docs/validation/viv_re150_predeclare.md` (agreed at
+  `1abfe954` under its former name `leontini_viv_predeclare.md`); run at `fb507aa1`, single case, D = 40, 75 D × 60 D.
+- Result: status `scored`, 336,384 steps (t U/D = 336.4); A_max = .5464 D vs
+  .520 D (+5.1 %, inside the declared ±10 % band); mean half peak-to-peak
+  .5448 D; f/f_n = .955; block drifts .78/.29/.18 %, last-five spreads
+  .15/.11/.10 %; finite, no_donor = 0, no seed; 337.7 s, ~1.0 ms/step.
+
+## Limitations (binding)
+
+- One resolution (D40) per case; no grid or time convergence.
+- Structural RHS = CME + Wen is a theory-motivated stage choice; CME-only
+  and node-impulse variants were not run as RHS; no force law is validated.
+- Boundary mismatches: NEQ outlet with backflow clamp vs traction outlet;
+  finite box; the strict Dütsch settling failure is unexplained.
+- No mass correction at curved walls; ~1 ms/step host-coupled overhead.

docs/validation/viv_re150_predeclare.md

@@ -1,4 +1,6 @@
-# Single-cylinder free transverse VIV — PREDECLARATION (agreed at 1abfe954)
+# Single-cylinder free transverse VIV, Re 150 (Kolahdouz/Qin) — PREDECLARATION (agreed at 1abfe954)
+
+(File renamed from `leontini_viv_predeclare.md`; the Leontini Re 200 case is only the queued second case below.)
 
 Stage goal (Yanqi via TRON): prescribed Dütsch case + one VIV case both
 "reasonably close" to a predeclared reference = development-stage success, then
@@ -27,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/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/leontini_viv_predeclare.md).
+"""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
@@ -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)
-            row = (n, y1, v1, f_rhs, float(d[1]), float(d[3]), float(f_obs),
-                   float(d[7] - d[10]), float(d[13]))
+            f_rhs = float(d[col["cme_y"]] + d[col["wen_y"]])   # CME + Wen (y)
+            row = (n, y1, v1, f_rhs, float(d[col["cme_y"]]), float(d[col["wen_y"]]), float(f_obs),
+                   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

research/viv_re150/viv.py

@@ -1,6 +1,6 @@
 """Pure logic of the predeclared Re 150 free transverse VIV case (stdlib only).
 
-Contract: docs/validation/leontini_viv_predeclare.md. Oscillator update with a
+Contract: docs/validation/viv_re150_predeclare.md. Oscillator update with a
 constant force over one lattice step, response-cycle tracking with fixed y = 0
 crossings, eligibility, 10-cycle blocks, settling/scoring with deterministic
 restart, and the conditional seed.

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
-the device: the caller supplies upload/step/synchronize callables.
+separate facts; centres are always ``base + offset``. MovingCircle owns host
+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,23 +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:
-        """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:
@@ -138,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``.
 

src/CelerisLab/configs/CONFIG.md

@@ -86,6 +86,20 @@ step/one_step_*.cu         → kernel 编排
 - 转速 `omega` 属于运行时状态，可在步进期间更新；该更新不改变 `N_OBJS` 或 `OBS_*` 内存布局，因此**不需要**因 `omega` 变化触发重新编译。
 - 力矩观测从 `obs` 的 torque segment 读取；force / torque / sensor 段统一由 `config_obs.h` 宏描述。
 
+### 外边界能力（Outer boundaries, dense `Simulation`）
+
+The dense path supports one channel package, not four independently
+selectable sides:
+- left x = 0: velocity inlet (`method.inlet.scheme`, `profile`, `physics.velocity`);
+- right x = nx−1: outlet (`method.outlet.mode`; `neq_extrap` holds rho = 1 with
+  a `max(u, 0)` backflow clamp when `backflow_clamp` is set);
+- bottom and top: one shared `y_wall_bc` (`bounce_back` halfway or
+  `free_slip`); corners are walls.
+`y_wall_bc = "on_site_velocity"` switches to x-periodic with on-site y planes
+and requires an object-free domain. There is no x/y-periodic or symmetry
+option for domains with bodies; adding outer BC types is a separate design
+question, not covered by these keys.
+
 ### 稳定性默认窗口（Stability defaults）
 
 - 默认 `method.omega_guard = {min: 0.01, max: 1.99}`，用于约束碰撞频率与 LES 有效粘度路径。

src/CelerisLab/lbm/stepper.py

@@ -21,7 +21,7 @@ class LBMStepper:
 
         self._esopull = (cfg.streaming == "esopull")
         # Set by the owner when a step would be wrong (translating circles
-        # excluded from the packed links before the moving kernels exist).
+        # without attached moving tables, a pending target, or a failed step).
         self.blocked: str | None = None
         self.moving = None  # MovingTables attached by the owner (translating circles)
 
@@ -39,8 +39,10 @@ class LBMStepper:
         self.apply_esopull_curved_fn = module.get_function("ApplyEsoPullCurvedKernel")
         self.sensor_fn = module.get_function("SensorKernel")
         self.force_region_fn = module.get_function("ForceRegionKernel")
-        self.moving_fn = None if self._esopull else module.get_function("MovingCircleLinks")
-        self.moving_cells_fn = None if self._esopull else module.get_function("MovingCircleCells")
+        # Moving-circle kernels exist only in D2Q9 double-buffer modules.
+        moving = not self._esopull and cfg.is_d2q9
+        self.moving_fn = module.get_function("MovingCircleLinks") if moving else None
+        self.moving_cells_fn = module.get_function("MovingCircleCells") if moving else None
 
         tpb = cfg.threads_per_block
         self.block = (tpb, 1, 1)

src/CelerisLab/simulation.py

@@ -28,7 +28,7 @@ from .cuda.context import CudaContext, context_bound
 from .cuda import compiler
 from .lbm.descriptors import FLUID
 
-# Translating circles step only through their moving tables (fixed centre in S2a).
+# Translating circles step only through their moving tables (attached at initialize).
 PENDING_BLOCKED = "a pending translation target is advanced only by Simulation.run(1)"
 MOVING_BLOCKED = ("translating circles need the moving kernels and their tables (attached "
                   "at initialize); "
@@ -122,8 +122,9 @@ class Simulation:
                 D2Q9/FP32/unshifted/DB only; set before initialization or stage
                 a new body. No runtime motion or amplitude setter.
             translation: ``body.motion.Translation`` declaring a translating
-                circle (before initialize only; docs/moving_bodies.md). Until
-                the moving kernels exist, positive-step runs are rejected.
+                circle (before initialize only; docs/moving_bodies.md).
+                Positive-step runs are rejected until ``initialize()`` has
+                attached its moving tables.
             **kwargs: passed to the geometry constructor.
 
         Returns:
@@ -241,7 +242,8 @@ class Simulation:
         self._require_usable()
         tables = self.stepper.moving
         if tables is None:
-            return np.zeros((0, 6), np.float32)
+            from .body.motion import MovingTables
+            return np.zeros((0, MovingTables.DIAG_FLOATS), np.float32)
         self.stream.synchronize()
         return tables.download_diagnostics()
 
@@ -360,8 +362,11 @@ class Simulation:
         Args:
             id: body_id from add_body().
             omega: rotation speed (lattice units).  ``None`` = no change.
-            vx: translation velocity x (future).  ``None`` = no change.
-            vy: translation velocity y (future).  ``None`` = no change.
+            vx: translation velocity x. For a translating circle: the
+                interval wall velocity of the target, only together with
+                ``offset``; otherwise stored in the body state (not applied
+                to stationary walls).  ``None`` = no change.
+            vy: translation velocity y; same rules as ``vx``.
             offset: target centre offset of a translating circle; with vx/vy
                 (default 0) it replaces the pending target after validation
                 against the applied offset (docs/moving_bodies.md, 4a). Only
@@ -781,6 +786,7 @@ class Simulation:
         (docs/moving_bodies.md 4a, 5). Any failure leaves the Simulation unusable."""
         tables = self.stepper.moving
         circles = [(obj.motion, obj.obj_id) for obj in self._moving_bodies()]
+        pending = any(motion.pending is not None for motion, _ in circles)
         try:
             advancing = tables.upload_records(circles, stream)
             self.stepper._step(steps, action_gpu=self.bodies.action_gpu,
@@ -793,10 +799,12 @@ class Simulation:
                 no_donor = tables.download_diagnostics()[:, tables.DIAG.index("no_donor")]
                 if np.any(no_donor > 0):
                     raise RuntimeError("translating circle refill found no donor cell")
+            if pending:  # a target equal to the applied offset is consumed as a hold step
                 for motion, _ in circles:
                     motion.commit()
-                tables.upload_records(circles, stream)  # first hold: old = new = applied, u = 0
-                stream.synchronize()
+                if advancing:
+                    tables.upload_records(circles, stream)  # first hold: old = new = applied, u = 0
+                    stream.synchronize()
                 self._attach_moving()
         except BaseException as exc:
             for motion, _ in circles:  # device geometry/state unknown after any failure

tests/unit/test_duetsch_protocol.py

@@ -11,6 +11,20 @@ _spec = importlib.util.spec_from_file_location(
 p = importlib.util.module_from_spec(_spec)
 _spec.loader.exec_module(p)
 
+def _load(rel, name):
+    """Load a research driver module under a unique name (no bare ``import run``)."""
+    import sys
+    path = Path(__file__).resolve().parents[2] / rel
+    sys.path.insert(0, str(path.parent))  # for the driver's own sibling imports
+    try:
+        spec = importlib.util.spec_from_file_location(name, path)
+        mod = importlib.util.module_from_spec(spec)
+        spec.loader.exec_module(mod)
+    finally:
+        sys.path.remove(str(path.parent))
+    return mod
+
+
 T = 2520
 
 
@@ -68,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):
@@ -93,12 +136,7 @@ def test_exact_point_inside_body_is_missing_even_with_fluid_corners():
 
 def test_generated_config_passes_the_real_loader(tmp_path):
     import json
-    import sys
-    sys.path.insert(0, str(Path(__file__).resolve().parents[2] / "research/duetsch_re100_kc5"))
-    try:
-        import run
-    finally:
-        sys.path.pop(0)
+    run = _load("research/duetsch_re100_kc5/run.py", "duetsch_run")
     from CelerisLab.config import load_lbm_config
     cfg_dict, u_max, n = run.config(20, 10.0, 2520)
     path = tmp_path / "config_lbm.json"
@@ -111,13 +149,7 @@ def test_generated_config_passes_the_real_loader(tmp_path):
 
 
 def _analyze():
-    import sys
-    sys.path.insert(0, str(Path(__file__).resolve().parents[2] / "research/duetsch_re100_kc5"))
-    try:
-        import analyze
-    finally:
-        sys.path.pop(0)
-    return analyze
+    return _load("research/duetsch_re100_kc5/analyze.py", "duetsch_analyze")
 
 
 def test_hf_energy_is_parseval_weighted():

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
-moving body must be rejected before any device write.
+Real Simulation entrypoints with device doubles; no CUDA context or GPU work.
 """
+import math
 from types import SimpleNamespace
 
 import pytest
 
-from CelerisLab.body.motion import (MAX_STEP, R_MIN, MovingCircle, Translation,
-                                    enumerate_donor_coverage)
+from CelerisLab.body.motion import MAX_STEP, R_MIN, MovingCircle, Translation
 
 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():
-    c = circle()
-    assert c.record() == ((50., 40.), (50., 40.), (0., 0.))  # hold: no translation velocity
-    c.set_target((0., .25), (0., .05))
+def test_record_commit_and_hold(monkeypatch):
+    sim, writes, c = routed_advancing_simulation(monkeypatch)
+    assert c.record() == ((50., 40.), (50., 40.), (0., 0.))
+    sim.set_body(0, offset=(0., .25), vy=.05)
     assert c.record() == ((50., 40.), (50., 40.25), (0., .05))
-    uploads = []
-    c.advance(uploads.append, lambda: None, lambda: None)
-    assert uploads == [((50., 40.), (50., 40.25), (0., .05))]
+    seen_at_sync = []
+    sim.run(1, stream=SimpleNamespace(synchronize=lambda: seen_at_sync.append(c.applied)))
+    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):
-    c = circle()
-    c.set_target((0., .25), (0., .05))
+def test_failure_after_upload_attempt_is_fail_closed(monkeypatch, stage):
+    sim, writes, c = routed_advancing_simulation(monkeypatch, fail="upload" if stage == "upload" else None)
+    sim.set_body(0, offset=(0., .25), vy=.05)
 
-    def call(name):
-        def run(*_):
-            if name == stage:
-                raise RuntimeError(f"{name} failed")
-        return run
+    def fail(*a, **k):
+        raise RuntimeError(f"{stage} failed")
+    if stage == "step":
+        sim.stepper._step = fail
+    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"))
-    assert c.failed and c.applied == (0., 0.)
-    for later in (c.record, lambda: c.set_target((0., .1), (0., 0.)),
-                  lambda: c.advance(*(lambda *_: None,) * 3)):
+        sim.run(1, stream=stream)
+    assert c.failed and c.applied == (0., 0.) and c.pending == ((0., .25), (0., .05))
+    assert sim._failed is not None and sim.stepper.blocked == sim._failed
+    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))
@@ -297,6 +366,44 @@ def test_run_one_advances_and_commits_after_sync(monkeypatch):
     assert ("step", 3, False) in writes and motion.applied == (0., .25)
 
 
+def test_zero_motion_target_is_consumed_by_run_one(monkeypatch):
+    sim, writes, motion = advancing_simulation(monkeypatch)
+    motion.set_target((0., 0.), (0., 0.))          # equal to applied: a hold target
+    sim._attach_moving()
+    assert sim.stepper.blocked == module.PENDING_BLOCKED
+    sim.run(1, stream=STREAM)
+    assert motion.pending is None and motion.applied == (0., 0.)
+    assert sim.stepper.blocked is None and ("step", 1, False) in writes
+    writes.clear()
+    sim.run(3, stream=STREAM)                      # no stale latch
+    assert ("step", 3, False) in writes
+
+
+@pytest.mark.parametrize("streaming,d2q9,bound", [("double_buffer", True, True),
+                                                  ("double_buffer", False, False),
+                                                  ("esopull", True, False)])
+def test_moving_kernels_bound_only_for_d2q9_double_buffer(monkeypatch, streaming, d2q9, bound):
+    from CelerisLab.lbm import stepper as stepper_module
+    monkeypatch.setattr(stepper_module.cuda, "mem_alloc", lambda n: object())
+    moving = {"MovingCircleLinks", "MovingCircleCells"}
+
+    def get_function(name):
+        if name in moving and not (d2q9 and streaming == "double_buffer"):
+            raise KeyError(name)                 # absent from D3Q19 / EsoPull modules
+        return lambda *a, **k: None
+    cfg = SimpleNamespace(streaming=streaming, threads_per_block=128, is_d2q9=d2q9,
+                          nx=16, ny=8, nz=4, nq=9 if d2q9 else 19)
+    stepper = stepper_module.LBMStepper(SimpleNamespace(), SimpleNamespace(get_function=get_function), cfg)
+    assert (stepper.moving_fn is not None) == bound and (stepper.moving_cells_fn is not None) == bound
+
+
+def test_moving_diagnostics_empty_shape_matches_layout(monkeypatch):
+    from CelerisLab.body.motion import MovingTables
+    sim = make_simulation(monkeypatch)
+    sim._failed, sim.stepper = None, SimpleNamespace(moving=None)
+    assert sim.moving_diagnostics().shape == (0, MovingTables.DIAG_FLOATS)
+
+
 @pytest.mark.parametrize("failure", ("upload", "first_sync", "no_donor", "hold_upload",
                                      "hold_sync"))
 def test_failed_advance_fails_simulation(monkeypatch, failure):
@@ -507,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))
-    motion = sim.bodies.get(bid).motion
-    motion.set_target((0., .2), (0., .05))
+    sim.set_body(0, offset=(0., .2), vy=.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()

tests/unit/test_viv_logic.py

@@ -8,6 +8,20 @@ _spec = importlib.util.spec_from_file_location(
 viv = importlib.util.module_from_spec(_spec)
 _spec.loader.exec_module(viv)
 
+def _load(rel, name):
+    """Load a research driver module under a unique name (no bare ``import run``)."""
+    import sys
+    path = Path(__file__).resolve().parents[2] / rel
+    sys.path.insert(0, str(path.parent))  # for the driver's own sibling imports
+    try:
+        spec = importlib.util.spec_from_file_location(name, path)
+        mod = importlib.util.module_from_spec(spec)
+        spec.loader.exec_module(mod)
+    finally:
+        sys.path.remove(str(path.parent))
+    return mod
+
+
 D, U = 40, 0.04
 
 
@@ -72,12 +86,7 @@ def test_seed_reset_discards_incomplete_cycle_and_proof():
 
 def test_viv_config_passes_the_real_loader(tmp_path):
     import json
-    import sys
-    sys.path.insert(0, str(Path(__file__).resolve().parents[2] / "research/viv_re150"))
-    try:
-        import run
-    finally:
-        sys.path.pop(0)
+    run = _load("research/viv_re150/run.py", "viv_run")
     from CelerisLab.config import load_lbm_config
     cfg_dict, nx, ny = run.config()
     path = tmp_path / "config_lbm.json"