Predeclare zero-damping pinball development cases

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

@@ -0,0 +1,104 @@
+# Uncontrolled steady pinball: development protocol
+
+Frozen before results: 1280 x 512 D2Q9/FP32/MRT, unshifted double buffer;
+rho=1, mean U=.01, nu=.004, D=20, Re_D=50. Centres are (600,255.5),
+(626,270.5), (626,240.5), radius10. No upstream cylinder. Passive probes stay
+at (800,295.5), (800,255.5), (800,215.5), radius5. Every omega is zero.
+Parabolic inlet, regularized NEQ damping1; y bounce-back planes .5/510.5;
+NEQ-extrapolated outlet with backflow clamp (the configured blend is inactive in
+this MRT route). Initial fluid is
+parabolic equilibrium. Fixed warmup is 4 NX/U=512000 steps.
+
+Geometry and fluid parameters follow DynamisLab 609116ec26f19766f5e9ce2b0df427d4e24faae8,
+`legacy_test/core/legacy_env_builder.py:build_steady_cloak` and the retained
+`config_lbm_pinball_legacy_re50d_u001.json`. This is a modern Celeris reconstruction,
+not a bitwise Legacy reproduction: Legacy uses shifted incompressible MRT, a
+different curved interpolation and prescribed parabolic outlet velocity. Current
+Celeris uses standard MRT, two-branch BFL and extrapolated outlet velocity.
+The fixed-hold control uses the SAME modern settings as the free case.
+
+## Shared translation and ownership
+
+`Translation(..., shared=True)` is explicitly transverse only. Reusing one
+envelope INSTANCE declares common motion; equal but separate envelopes are not
+linked. Before initialization and every positive Simulation.run, group members
+must have identical applied offsets, effective target offsets and interval wall
+velocities. Stage each member with set_body, then call run(1) once. An incomplete
+or inconsistent target group is rejected before device uploads and can be fixed.
+Existing pending-target direct-stepper restrictions remain.
+
+Only members of that group may have overlapping whole swept bands. Their fixed
+surface clearance must exceed4 lattice cells. This conservative bound exceeds
+2 sqrt(2)+MAX_STEP = 3.329: M1 changed cells lie in old/new disks and read one-link
+neighbors; M2 reads at most two links outside its new disk. Old/new relative
+displacement is at most .5. Thus reads cannot reach another member's old/new
+solid or M1 write region. M1 writes are disjoint. M2 writes opposite population
+slots inside its own new disk, also disjoint. Separate first_slot ranges still
+allocate each body's logical band/debug storage. Shared band cells alone do not
+write anything: M1 returns unless membership changes; M2 returns unless a cut
+link exists. Unrelated moving/static body and probe overlaps still reject.
+
+This is a sufficient restricted-group condition, not support for near-contact
+or independent motion. The r4, distance8.3, dy=.5 counterexample remains rejected.
+CPU tests enumerate candidate donors/write footprints over subcell phases;
+GPU tests compare fixed vs hold and member-order permutations through crossings
+at wall speed .01 over the full +/-D envelope (8000 steps). The initial .5-speed
+fixture failed with NaNs; a single-circle control also failed at step69. That
+Mach~.866 fixture was unsuitable as a low-Mach numerical gate. Its failure is
+retained; the geometric MAX_STEP proof is unchanged, not a stability guarantee.
+CUDA kernels and arithmetic are unchanged.
+
+## New structural parameters (not inherited from the paper)
+
+One transverse state y,v; total assembly M=10 rho (3 pi D^2/4)=9424.77796,
+fn=U/(5D)=.0001, Tn=10000, zeta=.01, K=M(2 pi fn)^2,
+C=2 zeta M(2 pi fn). Initial y=v=0, no seed or later kick. Force is the sum
+of the three instantaneous CME_y+Wen_y values from the previous fluid step;
+node impulses are logged separately and never added. The damped oscillator is
+integrated exactly for that piecewise-constant force over dt=1. Wall velocity
+is y[n+1]-y[n], not the endpoint oscillator velocity. No filtering or retuning.
+
+Observe600000 steps=60Tn=300 D/U after warmup, in each of free and fixed-hold
+cases. Final20Tn is a descriptive comparison window; ten-Tn blocks expose drift.
+No settled/periodic claim follows merely from reaching the cap. Stop on |y|>D,
+nonfinite values, no_donor, rejected displacement or speed>.05. Report any quiet
+response honestly. No automatic extension. --smoke uses32 warmup+32 steps only.
+
+Native-step rows retain y,v, each body's CME/Wen/node terms and refill counters;
+each600000x15 FP64 file is72MB. No full-field output or policy training. Compile
+and run on one device serially, keep original outputs and dirty trees untouched.
+This stage makes no claim of physical grid convergence or vibration suppression.
+
+Run from checkout with PYTHONPATH=src and the existing Python environment:
+`python research/steady_pinball/run.py --smoke --out <new-external-directory>`
+then separately `--out <new-free-directory>` and `--fixed --out <new-fixed-directory>`.
+
+## Zero-damping development follow-up (predeclared 2026-10-31)
+
+The original d0bcefa run and its zeta=.01 results remain historical evidence.
+The current runner exposes three fixed `--case` choices, all with C=0 as requested:
+
+| Case | M / (rho 3 pi D^2/4) | fn | M | K | Purpose |
+|---|---:|---:|---:|---:|---|
+| Z0 | 10 | .0001 | 9424.77796 | .00372075320 | Remove damping only |
+| M2 | 2 | .0001 | 1884.95559 | .000744150640 | Reduce mass and stiffness together at fixed fn |
+| R2 | 2 | .00006 | 1884.95559 | .000267894231 | Change stiffness only relative to M2 |
+
+R2 deliberately places fn at the retained fixed-body lift peak, f=.00006
+(St=.12 on mean U). This is an informed development choice, not independent
+physical validation or a promise of lock-in. All geometries, fluid settings,
+512000-step fixed warmup, y0=v0=0, force time level and unfiltered RHS are unchanged.
+No run result may select a replacement parameter or an automatic extension.
+
+Each case has a 600000-step cap, the same 300 D/U physical duration. This equals
+60 natural periods for Z0/M2 and 36 for R2. Compare the common final200000 native
+steps, and six100000-step blocks, retaining actual natural/response frequencies
+and ongoing drift. Do not call those blocks ten R2 periods. No settled-VIV claim
+follows from a larger displacement. Report all three outcomes, including quiet
+responses or stopped cases, alongside the existing fixed and damped-free records.
+
+Unchanged stops: |y|>D, invalid move/speed, nonfinite values or no_donor. A 1800s
+per-case wall cap and GPU headroom check bound resource use. Preserve partial rows;
+do not add damping, clip trajectories, filter forces, or relax any stop. Each
+case adds72MB of scalar history, no full fields. The separate Legacy/modern
+baseline diagnosis does not change this FSI protocol.

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

@@ -0,0 +1,242 @@
+# Prescribed in-line oscillation: Dütsch Re 100 / KC 5 — PROTOCOL (review draft)
+
+Production under test: `claude/moving-bodies-s3a` at `653fd997` (unchanged).
+Primary source: Dütsch, Durst, Becker & Lienhart, J. Fluid Mech. 360, 249–271
+(1998), doi:10.1017/S002211209800860X. Source facts below were verified from
+the paper (Slack C0C78UFKGRM/1793335557.460509, TRON 1793356317–1793356724,
+page numbers given). Everything marked ENGINEERING is our a-priori choice,
+fixed before any benchmark result; it is not a paper uncertainty.
+
+## 0. Source facts (with provenance)
+
+- Motion (p.251): `x_c = -A sin(2 pi f t)`; paper phase phi = 0 is the mean
+  position moving left at `-U_max`. KC = 5 gives `A/D = 5/(2 pi) = .7957747`.
+- Start (p.254): the computation starts a quarter cycle before phi = 0, at
+  `x = +A` with zero velocity (phi = -90°).
+- Experiment (pp.252, 254–256): D = 10 mm, L = .65 m, tank 1.2 × .7 × .7 m,
+  rho = 998.2 kg/m³, mu = 1.00818 mPa·s, Re = 100, KC = 5; LDA with 1° phase
+  bins and 200 samples per node and degree; no quantified experimental error
+  recovered.
+- Profiles (Fig. 8, p.260): LAB-FIXED stations x1 = -6, 0, 6, 12 mm, i.e.
+  x/D = -.6, 0, .6, 1.2, at phi = 180°, 210°, 330°; u1 and u2 versus x2.
+- Original numerics (§2): cylinder-fixed accelerated frame (`x_i = x~_i +
+  x_i,s`, outer velocity `-u_c`, inertial source `-rho du_c/dt`), O-grid outer
+  diameter 120 D, illustrated 192 circumferential nodes, 720 steps/cycle; tank
+  influence neglected.
+- Force fit (verified): force ON THE BODY per unit span
+  `F1 = -0.5 rho D c_d xdot|xdot| - (pi/4) rho D² c_i xddot` with `x = -A sin phi`;
+  normalized by `0.5 rho D U_max²`: `c_d cos phi |cos phi| - (pi² c_i / KC) sin phi`.
+  `c_i` is the ADDED-MASS coefficient; total inertia `c_m = c_i + 1` is a
+  different convention and is not used.
+- Table 1 computed Morison coefficients (c_d, c_i) for 96 / 192 / 384 perimeter
+  nodes at dt/T = 1/720: (2.15, 1.42), (2.10, 1.45), (2.09, 1.45). These are
+  fitted computed coefficients, not a measured force and not an error bar. The
+  paper states the coarse direct force is closer to the finest direct force
+  than the Morison fit, but gives no waveform-error norm.
+
+## 1. Kinematics, time levels and phase bookkeeping
+
+- Our run is in the LAB frame: the circle moves, the fluid is initially at rest
+  (rho = 1, u = 0); no fictitious inertial source is added. Lab-fixed profile
+  stations are compared directly (no frame transform needed).
+- Geometry reference: the circle's base centre is the INITIAL position, the
+  right turning point `x_mean + A`; the declared offset box is x in [-2A, 0],
+  y = 0. The body is never teleported: offset 0 at s = 0.
+- Elapsed time s = n dt from the start. Analytic position relative to the mean
+  `x(s) = A cos(2 pi s / T)`, velocity `u(s) = -U_max sin(2 pi s / T)`; paper
+  phase `phi = (360 s / T - 90) mod 360`.
+- Step n moves from s_n to s_{n+1}: target offset `x(s_{n+1}) - A`; wall
+  velocity = the single-step interval average `(x(s_{n+1}) - x(s_n)) / dt`
+  (analytic, supplied, not differenced from stored offsets). The per-step
+  force of `run(1)` / `read_force(normalize=True)` is attributed to that same
+  interval.
+- Paper cycle k = 1, 2, ... spans `s in [T/4 + (k-1) T, T/4 + k T)`; the
+  initial quarter cycle (s < T/4) is cycle 0 and never scored. Every logged
+  sample carries (k, phi). Example: phi = 330° of cycle 1 is at s = 7T/6.
+- T is chosen as a multiple of 360 native steps, so phi = 180°, 210°, 330° fall
+  exactly on step boundaries and 1° bins hold whole steps.
+
+## 2. Resolution ladder (uniform grid, fixed Mach), ENGINEERING
+
+| D (cells) | T (steps) | U_max = 5 D / T | nu_lat = U_max D / 100 | tau (3 nu + .5) |
+|---|---|---|---|---|
+| 20 | 2,520 | .039683 | .0079365 | .5238 |
+| 40 | 5,040 | .039683 | .015873 | .5476 |
+| 80 | 10,080 | .039683 | .031746 | .5952 |
+
+- MRT collision. Mach fixed (Ma ≈ .0687); a separate Mach check at D = 40
+  with T = 10,080 (U_max halved, tau .5238) is labelled MACH and not part of the
+  ladder.
+- Per-step displacement ≤ U_max ≈ .04 cells (bound .5).
+
+## 3. Domain and boundaries (finite-domain reproduction), ENGINEERING
+
+The paper's 120 D O-grid is not reproduced; our rectangular box is a
+reproduction choice. The native dense channel (zero-velocity inlet,
+pressure-NEQ outlet, bounce-back y walls) is a finite-domain approximation,
+not infinite still fluid. The G4 x-periodic / y-slip override is a separately
+labelled sensitivity family and never mixed into the primary ladder. The
+runner owner confirms the zero-velocity inlet is supported as configured
+before any run.
+
+Order (closed):
+1. Smoke + ledger at D = 40, native BC, 40 D × 40 D (finite box).
+2. Domain sensitivity before any literature statement: same BC, D = 40,
+   L ∈ {20 D, 40 D, 60 D}; record the change in Fx waveform, c_d, c_i and the
+   scored profiles.
+3. Optional labelled BC sensitivity: one x-periodic / y-slip case at D = 40,
+   40 D; reported as a delta against native only.
+4. Only after step 2 passes: the D = 20 / 40 / 80 ladder on the chosen BC and
+   domain.
+
+Domain criteria (a-priori, ENGINEERING; separate from the Table 1 bands): the
+reference is 60 D at D = 40; 20 D and 40 D are compared with it on the common
+window (same cycle indices) and the common valid point set:
+- `RMS(Fx_L - Fx_60D) / F_ref ≤ .01` over the window (per-step waveforms);
+- `|c_d,L - c_d,60D| ≤ .02` and `|c_i,L - c_i,60D| ≤ .02` (window fits);
+- `max(|u_L - u_60D|, |v_L - v_60D|)` over common points `/ U_max ≤ .02`.
+The chosen domain is the smallest of {20 D, 40 D} passing all three; if
+neither passes, the domain verdict is FAIL and the ladder is not run.
+- Memory/runtime: measured by a smoke run before each new size; measurements
+  inform scheduling only, never tolerances.
+
+## 4. Observables (raw, per step, preserved)
+
+- `obs` Fx, Fy (CME) and all `MovingTables.DIAG` columns per step: CME,
+  CME + Wen and CME + cover − uncover reported separately; crossings, eq_only,
+  no_donor (no_donor must stay 0).
+- Velocity u1, u2 on lab-fixed lines x/D = -.6, 0, .6, 1.2 (relative to the
+  mean cylinder position) versus x2, at the scored phases. Lab frame
+  throughout; if any body-frame intermediate is used, position and velocity
+  are both transformed.
+- Sampling cost: raw force and diagnostics counters are read every step
+  (small). Flags and velocity fields are read only at the steps inside the
+  three prescribed 1° phase bins (φ ∈ [180, 181), [210, 211), [330, 331),
+  T/360 steps each) of the cycles needed for settling and scoring, plus
+  sparse declared audits; never every step.
+- Point set: fixed x2/D sample points per station, identical on every grid
+  and declared before running; per point the interpolation weights are
+  recorded; the scored set is the common subset valid on all compared grids
+  at that phase, with its coverage recorded; point sets never change silently.
+- Interpolation: bilinear on cell centres using only cells that are fluid at
+  that step (flags read back via `get_flags()`); weights renormalized over the
+  fluid corners; a sample with no fluid corner, or inside the circle at that
+  step, is MISSING and logged, never set to zero. One-degree bins keep valid
+  fluid weights. The numerical average need not mimic the 200 LDA samples.
+
+## 5. Periodic-regime gate (every case, before any statistic), ENGINEERING
+
+Ledger rows, same schema for every case: dx, dt (in D and T units); U_max, nu,
+Re, KC (lattice and nondimensional); T (native steps, nondimensional);
+duration (native steps, s/T, paper cycles); discard; accepted window; cycle
+ledger; verdict PASS or FAIL-CLOSED.
+
+Scales: `F_ref = 0.5 rho U_max² D`, `U_ref = U_max`; never a near-zero mean
+force as denominator.
+
+- Settling metrics, cycle k−1 → k, phase-aligned by step within the cycle
+  (error budget: one fifth of the tightest comparison band of the same
+  observable, so residual drift cannot obscure a domain or resolution pass):
+  - waveform `RMS(dFx) / F_ref ≤ .002`;
+  - `|d mean|`, `|d RMS|`, `|d peak|` of Fx, each `/ F_ref ≤ .002`;
+  - scored profiles (fluid-only) `|du| / U_max`, `|dv| / U_max ≤ .004` at the
+    three phases × four stations;
+  - Morison `|dc_d| ≤ .004`, `|dc_i| ≤ .004`.
+- Conservative comparison: every domain and cross-resolution criterion is
+  applied to `|difference| + drift_a + drift_b`, where drift is the largest
+  pair metric of that case inside the scored window.
+- These constants (protocol.py SETTLE, DOMAIN, CROSS) are frozen before any
+  physical result; no post-result loosening.
+- Rule: all metrics for N = 3 consecutive whole paper cycles → settled.
+
+Cycle indexing (worked example, binding). Paper cycle k spans
+`s/T in [0.25 + (k-1), 0.25 + k)`, i.e. φ0 → φ0; cycle 0 is the startup
+quarter `s/T < 0.25`. Pair p_k compares cycle k−1 with cycle k.
+
+- Discarded: cycles 0, 1, 2 (s/T < 2.25). Pairs are formed only between
+  cycles ≥ 3, so the first pair is p_4 (3 → 4).
+- Settling proof for a case: three consecutive pairs all pass; K_case = the
+  first cycle AFTER the last of those pairs. Earliest possible: pairs 3→4,
+  4→5, 5→6 pass ⇒ K_case = 7.
+- Common window: K = max over compared cases of K_case; every case scores the
+  same three complete cycles K, K+1, K+2 and must stay settled continuously:
+  every pair p_j for j = its own K_case .. K+2 must pass (including the wait
+  between its local settle and the common K), otherwise that case is
+  FAIL-CLOSED.
+- Cap: cycle k ends at s/T = k + 0.25, so the last admissible scored cycle is
+  11 (ends 11.25 T) and K ≤ 9. Examples: K = 7 → cycles 7–9 (ends 9.25 T);
+  K = 9 → cycles 9–11 (ends 11.25 T); any K_case > 9 → FAIL-CLOSED.
+- Run length: cases run sequentially; each runs to the end of cycle 11
+  (11.25 T: 28,350 / 56,700 / 113,400 steps for D = 20 / 40 / 80; the 12 T
+  cap is 30,240 / 60,480 / 120,960) and keeps its scalar and profile data, then
+  the common K is chosen. An unscorable partial cycle 12 is not run.
+- Profiles: an empty common valid point set is FAIL-CLOSED; coverage, missing
+  points, bin time weights (equal per step) and the original fluid corner
+  weight sums are reported.
+- The earlier G4 timing check (100 coarse vs 400 diffusive-fine steps) was a
+  matched short transient, not periodic evidence.
+
+## 6. Metrics and per-observable criteria (ENGINEERING, not paper uncertainty)
+
+Three criteria stay separate; the G4 same-grid 2% target is not reused.
+
+- Morison fit: least squares of the normalized per-step force over the
+  accepted window against the two basis functions, each averaged over its
+  single step interval, fixed paper phase, no free phase shift; c_d and
+  added-mass c_i per force candidate.
+- Cross-resolution: every scored observable's error against its reference must
+  decrease monotonically D20 → D40 → D80 (necessary, not sufficient).
+- Against Table 1 (computed fit, not measured force): at D ≥ 40,
+  `|c_d - 2.09| ≤ .15` and `|c_i - 1.45| ≤ .15`, and monotone.
+- Force waveform against its own Morison reconstruction: residual reported
+  only, no pass/fail (the paper reports fit error above its coarse numerical
+  error and gives no percentage).
+- Fig. 8 LDA profiles: QUALITATIVE ONLY. The original figure could not be
+  accessed or rendered; no digitized ordinates and no measurement uncertainty
+  exist. After a digitization with figure number, tool and uncertainty note:
+  `RMS(|du|, |dv|) / U_max ≤ .15` at D ≥ 40 and monotone; LDA error stays
+  UNQUANTIFIED and is stated with every comparison.
+- Fy sanity (not literature): `|mean Fy| / F_ref ≤ .05` in the accepted window.
+- No force candidate is selected for an RHS by these results; each is
+  reported.
+
+Superseded (earlier drafts): the velocity-jump start, guessed stations and
+phases, "C_M" for c_i, and the uncertainty wording; cycles 4–6, the 30-cycle
+cap and the earlier tolerance numbers.
+
+## 6b. Amendment: duration-only diagnostic (declared before the run)
+
+- Background: the D40 domain stage (20 D / 40 D / 60 D, 11.25 T each) failed
+  settling in every box (no window). That result stands.
+- One declared diagnostic: D = 40, 60 D box, from rest, through cycle 31
+  (157,500 steps = 31.25 T). Physics, settling tolerances, bands and
+  production unchanged; three proof pairs plus a continuously settled
+  K..K+2 with K ≤ 29.
+- Invocation: `run.py --D 40 --box 60 --cycles 31` then
+  `analyze.py --horizon 31`; the JSON records `horizon` and `horizon_mode`.
+  Without `--horizon` the analyzer is the original cap-11 scoring.
+- Failure at 31.25 T: report and stop; no automatic extension, matrix, D80
+  or VIV. A pass shows only that this box settles with time; it is not a
+  domain or literature pass.
+- The phase-locked high-frequency force component (residual of the mean
+  waveform after its Morison reconstruction, above 20/T) is reported per case
+  as unqualified; it needs separate qualification before any Morison or
+  literature statement.
+
+## 6a. Driver unit tests (CPU, before any run)
+
+- Phase/cycle arithmetic: s ↔ (k, φ), cycle boundaries, bin membership,
+  window/cap selection for K = 7, 8, 9, 10.
+- Solid-aware interpolation: all-fluid, partial-solid (renormalized weights),
+  no-fluid (MISSING) and inside-circle cases.
+- Interval-averaged Morison bases and the two-parameter fit on a synthetic
+  waveform with known c_d, c_i.
+
+## 7. Order and scope
+
+Protocol review → one-shot Dütsch driver (analytic `set_body` + `run(1)`,
+cycle ledger, lab-fixed fluid-only station sampler, two-basis Morison fit; no
+reusable runner until the first smoke shows shared helpers pay off; no
+production change) → smoke (memory/runtime) → domain sensitivity →
+ladder → physical-validation verdict → only then single-cylinder VIV. No AMR,
+no pinball, no manuscript work.

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

@@ -0,0 +1,155 @@
+# 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
+cleanup. The frozen Dütsch settling failure (0/28 pairs at 31.25 T) stays on
+record. Practical closeness is not a formal validation claim. Agreed by
+Claude, Grok and Astra at 1abfe954; implemented in research/viv_re150/.
+
+## Primary case (one case)
+
+- Reference: Qin et al., arXiv:2003.12186, §4.5 (maximum amplitudes at
+  U_red = 3…8: .075, .551, .520, .459, .370, .084 D in the text, as quoted by
+  Grok/Astra) for the target; setup from Kolahdouz et al., arXiv:2003.12046,
+  §4.1.1 (read 2026-10-30 via ar5iv): U* sweep 3–8 with amplitudes only in
+  Fig. 3(b) (no text values), domain [−30, 45] D × [−30, 30] D with the
+  cylinder "initially at rest and centered at the origin", uniform inflow on
+  the left, zero normal traction and tangential velocity at the outlet, zero
+  normal velocity and tangential traction on y = ±30 D; no initial
+  perturbation and no total simulated time stated.
+- Re = 150, U* = U/(f_n D) = 5, m/(rho D²) = 2 (m* = rho_s/rho_f = 8/pi),
+  zeta = 0, f_n = sqrt(K/M)/(2 pi) in vacuum. Transverse DOF only.
+- Target: A_max = .520 D (peak amplitude, Qin §4.5). Multi-code reference
+  (Kolahdouz ILE-NS, Qin II-LBM). Reference uncertainty: not stated; the
+  spread between codes is not quantified here.
+
+## Our reproduction (declared differences)
+
+- D = 40, MRT, dense double buffer; U_lat = .04, nu = U D / 150 = .010667,
+  tau = .532.
+- 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
+  declared cost cut with more blockage, not "the same case".
+- Allowed y offset ±1.5 D; peak transverse speed ≈ .026 cells/step (bound .5).
+
+## Structure, coupling and time levels
+
+- M = 2 rho D² (lattice rho = 1), f_n = U/(5 D) = 2e-4 per step, K = M (2 pi f_n)²,
+  zeta = 0. Equation of motion uses the structural mass only.
+- Per-step sequence (binding): known (y_n, v_n) and the last completed force
+  F_last → candidate (y_{n+1}, v_{n+1}) by the update below → `set_body` with
+  offset y_{n+1} and wall velocity (y_{n+1} − y_n)/dt with lattice dt = 1 (interval average; a
+  different quantity from the endpoint v_{n+1}) → `run(1)` → only on success:
+  commit y_n ← y_{n+1}, v_n ← v_{n+1}, F_last ← force read for that step (failure
+  semantics: see the window section). Whether the
+  measured force represents the interval's start, midpoint or end is not
+  asserted; the one-step lag is declared, not corrected.
+- Initial state: fluid initialized by `Simulation.initialize()` to the uniform
+  inlet velocity U everywhere (the code's actual initial condition; the
+  source gives no fluid initial condition), body at rest at the declared
+  centre, y_0 = v_0 = 0, F_last = 0.
+- Structural update over one step with constant force F = F_last (exact
+  harmonic solution, zeta = 0): y_e = F/K, w = 2 pi f_n,
+  y^{n+1} = y_e + (y^n − y_e) cos w + (v^n / w) sin w,
+  v^{n+1} = −(y^n − y_e) w sin w + v^n cos w.
+- Primary RHS (single coupled run): CME + Wen (theory-motivated stage choice:
+  Galilean-invariant link correction; G4 could not separate it from CME).
+  Raw CME and the node terms are logged as diagnostics of that trajectory;
+  the node impulse is not added. A run with CME-only as RHS is a separate,
+  labelled sensitivity case, decided separately. Neither RHS is claimed proven.
+
+## Window, failure semantics, acceptance (engineering, development stage)
+
+Clocks: lattice dt = 1; T_n = 1/f_n = 5,000 steps. Response cycles are
+counted from the measured signal; T_n is used only for the seed check.
+
+- Start: body at rest at the declared centre (the only initial condition the
+  source states), y_0 = v_0 = 0, F_last = 0; fluid uniform at U (as
+  initialized by the code; our choice, not the source's).
+- Seed (conditional, one time, declared): at step n_seed = 200,000
+  (t = 40 T_n), if max |y_t| over states t in [175,000, 200,000) is below
+  .01 D, then
+  in that step, after the check and before the candidate update, set
+  v_n += 2 pi f_n (.01 D) (F_last unchanged) and proceed normally; the
+  geometry moves continuously through the update. `seed_applied`, step and
+  size are logged; never repeated or changed.
+- Cycle boundary (discrete, fixed threshold): step n closes a cycle iff
+  y[n−1] < 0 and y[n] ≥ 0 (upward crossing of the equilibrium y = 0);
+  y_0 = 0 is not a crossing; no boundary fewer than 2 samples after the
+  previous one.
+- Per cycle: positive peak max(y), negative peak min(y), mean offset
+  (time-mean of y), amplitude a = (max − min)/2 (the only amplitude used for
+  settling).
+- Eligibility: a completed cycle is counted (toward the settling proof and
+  the 150 budget) only if a ≥ .01 D; ineligible cycles are logged, never
+  counted, so near-zero motion cannot settle or pass.
+- Seed reset: if the seed is applied, the incomplete cycle at n_seed is
+  discarded and the settling proof starts empty; the 150 budget and the
+  native cap keep running from t = 0.
+- Windows: adjacent non-overlapping blocks of 10 counted cycles;
+  W_k = mean amplitude of block k (W_k ≥ .01 D by eligibility, so no zero
+  denominator).
+- Block rule R(k), k ≥ 2: |W_k − W_{k−1}| / W_{k−1} ≤ .02 and, for the last 5
+  cycles of block k, max_i |a_i − ā_5| / ā_5 ≤ .01 with ā_5 their mean.
+- Settled at the end of block k when R(k) holds. Scoring = the next two
+  blocks k+1, k+2 (20 counted cycles); both must satisfy R. If one fails,
+  that block becomes the new reference: the settling proof restarts from it
+  and continues within the caps (deterministic; no other retry). A pass
+  requires one complete stable 20-cycle scoring window before a cap; success
+  stops after its 20th cycle.
+- Caps (both always armed; the first hit wins and is reported): 150 counted
+  response cycles including the 20 scoring cycles (so settling must complete
+  by counted cycle 130) → `timeout_response_cycles`; 1,000,000 native steps
+  (200 T_n) → `timeout_native_steps`.
+- Durations reported: total steps, t U/D, counted and uncounted cycles,
+  measured T_resp (mean boundary interval); N_resp / f_n only as a labelled
+  secondary clock.
+- Growth is not a failure: initial VIV growth is physical. Checked every
+  step (FAIL-CLOSED): non-finite y, v or force; |y| > 1.0 D (the tighter
+  gate; about twice the largest reference peak, .551 D) or the ±1.5 D offset
+  box reached; no_donor > 0; a failed upload/step/sync.
+- On any failure the driver stops and reports; it never retries,
+  sub-iterates, retunes the RHS or parameters, or starts a second case.
+  Device-side failures (upload/step/sync, no_donor) mark the Simulation
+  failed through the existing lifecycle; input rejection by `set_body`
+  (offset box, speed, step bound) and timeouts end the driver run without
+  that mark.
+- Target metric over the 20 scoring cycles: A_max/D = max over those cycles
+  of max(|positive peak|, |negative peak|) about y = 0 (Qin reports "maximum
+  amplitude"); also reported: mean amplitude, positive and negative peaks,
+  mean offset, f/f_n (no target).
+- Acceptance: |A_max − .520 D| ≤ .052 D (±10 %), stated together with the
+  boundary mismatch and the absence of a quoted reference uncertainty. A pass
+  is "reasonably close" for the stage: not grid convergence, not force-RHS
+  validation, and it does not touch the Dütsch settling record.
+- Driver reads force and the small diagnostics per step only; field reads
+  only at a few predeclared snapshots. No runtime claim: there is no measured
+  timing on this grid yet (the D40/60D receipt is a different, sampled case).
+
+## Plumbing smoke (not the benchmark)
+
+`run.py --smoke --D 20 --steps 32` runs exactly 32 steps of the same BC
+package and coupling loop (no seed, settle or score; never scored);
+status `smoke_ok` iff every step is finite, no_donor stays 0 and the body
+stays inside its box. The physical case is only the default invocation
+(`--D 40`, no `--smoke`).
+
+## Queued, not in this stage
+
+Leontini et al. JFS 22 (2006) Re = 200, m* = 10, zeta = 0.01, U* = 5.0 (needs
+the owners' c_crit convention and a Fig. 1 digitization). D80, literature
+claims, pinball, AMR, BC redesign, threshold changes.

research/duetsch_re100_kc5/analyze.py

@@ -0,0 +1,224 @@
+"""Ledger, settling, common window and Morison fits for Dütsch cases (numpy only).
+
+    python analyze.py CASE_DIR [CASE_DIR ...]   # cases of one comparison
+"""
+from __future__ import annotations
+
+import json
+import math
+import os
+import sys
+
+import numpy as np
+
+import protocol as p
+
+
+def load(path):
+    with open(os.path.join(path, "case.json")) as fh:
+        case = json.load(fh)
+    raw = np.load(os.path.join(path, "raw.npz"))
+    case.update(force=raw["force"], diag=raw["diag"], samples=raw["samples"], path=path)
+    return case
+
+
+def cycle_slices(case):
+    """{k: slice of step indices m whose interval [m, m+1] lies in paper cycle k}."""
+    T, q = case["T"], case["T"] // 4
+    return {k: slice(q + (k - 1) * T, q + k * T)
+            for k in range(1, case["cycles"] + 1) if q + k * T <= len(case["force"])}
+
+
+def coarse_wave(case, k, r=1):
+    """Normalized per-step Fx of cycle k, averaged over r native steps per sample."""
+    f = case["force"][cycle_slices(case)[k], 0] / case["F_ref"]
+    return f.reshape(-1, r).mean(axis=1)
+
+
+def profiles(case, k):
+    """{(phase, station, x2): (u, v)} of cycle k, time-averaged over the bin steps
+    (equal weights per step) from valid samples only."""
+    s = case["samples"]
+    acc = {}
+    for _, _, ph, xs, x2, u, v, _ in s[s[:, 1] == k]:
+        if not (math.isnan(u) or math.isnan(v)):
+            acc.setdefault((ph, xs, x2), []).append((u, v))
+    return {key: tuple(np.mean(vals, axis=0)) for key, vals in acc.items()}
+
+
+def coverage(case, k):
+    s = case["samples"]
+    rows = s[s[:, 1] == k]
+    missing = int(np.isnan(rows[:, 5]).sum())
+    return dict(samples=len(rows), missing=missing, bin_steps=len(set(rows[:, 0])),
+                min_fluid_weight=float(np.nanmin(np.where(rows[:, 7] > 0, rows[:, 7], np.nan)))
+                if len(rows) else math.nan)
+
+
+def common_points(dicts):
+    keys = set(dicts[0])
+    for d in dicts[1:]:
+        keys &= set(d)
+    return sorted(keys)
+
+
+def candidates(case):
+    """Per-step Fx of each reported force candidate (none is selected as an RHS)."""
+    d = case["diag"]
+    return {"cme": d[:, 0], "cme_wen": d[:, 0] + d[:, 2], "cme_node": d[:, 0] + d[:, 6] - d[:, 9]}
+
+
+def cycle_morison(case, ks, fx=None):
+    T = case["T"]
+    fx = case["force"][:, 0] if fx is None else fx
+    idx = np.concatenate([np.arange(len(fx))[cycle_slices(case)[k]] for k in ks])
+    b = np.array([p.interval_bases(int(m), T) for m in idx])
+    return p.morison_fit(fx[idx] / case["F_ref"], b[:, 0], b[:, 1])
+
+
+def pair_metrics(case, k):
+    a, b = coarse_wave(case, k - 1), coarse_wave(case, k)
+    pa, pb = profiles(case, k - 1), profiles(case, k)
+    common = common_points([pa, pb])  # empty set -> inf -> pair fails (FAIL-CLOSED)
+    prof = max((max(abs(pa[c][0] - pb[c][0]), abs(pa[c][1] - pb[c][1])) for c in common),
+               default=math.inf) / case["u_max"]
+    (cd0, ci0), (cd1, ci1) = cycle_morison(case, [k - 1]), cycle_morison(case, [k])
+    m = dict(wave=float(np.sqrt(np.mean((b - a) ** 2))),
+             scalar=float(max(abs(b.mean() - a.mean()), abs(np.sqrt(np.mean(b * b)) - np.sqrt(np.mean(a * a))),
+                              abs(np.abs(b).max() - np.abs(a).max()))),
+             profile=float(prof), morison=float(max(abs(cd1 - cd0), abs(ci1 - ci0))),
+             profile_points=len(common))
+    m["pass"] = all(m[key] <= tol for key, tol in p.SETTLE.items())
+    return m
+
+
+def hf_component(case, ks):
+    """Unqualified phase-locked high-frequency force component (reported, not checked):
+    mean CME waveform over cycles ks minus its two-basis Morison reconstruction,
+    split at 20/T; cycle differencing does not remove it."""
+    T = case["T"]
+    if not ks:
+        return dict(label="unqualified phase-locked HF force component, not checked", cycles=[])
+    mean = np.mean([coarse_wave(case, k) for k in ks], axis=0)
+    c_d, c_i = cycle_morison(case, ks)
+    q = T // 4
+    b = np.array([p.interval_bases(q + m, T) for m in range(T)])
+    resid = mean - (c_d * b[:, 0] - (np.pi ** 2 * c_i / p.KC) * b[:, 1])
+    spec = np.fft.rfft(resid)
+    high = spec.copy()
+    high[:21] = 0.0
+    energy = one_sided_energy(spec, T)
+    return dict(label="unqualified phase-locked HF force component, not checked",
+                cycles=list(ks), hf_rms_over_F_ref=float(np.sqrt(np.mean(np.fft.irfft(high, T) ** 2))),
+                hf_energy_fraction=float(energy[21:].sum() / max(energy.sum(), 1e-300)),
+                top_harmonics=[int(j) for j in np.argsort(energy[21:])[::-1][:5] + 21])
+
+
+def one_sided_energy(spec, n):
+    """Parseval-consistent energy per rfft bin: interior bins count twice, DC and
+    (for even n) Nyquist once."""
+    w = np.full(len(spec), 2.0)
+    w[0] = 1.0
+    if n % 2 == 0:
+        w[-1] = 1.0
+    return w * np.abs(spec) ** 2
+
+
+def ledger(case, horizon=p.LAST_SCORED):
+    ks = [k for k in sorted(cycle_slices(case)) if k <= horizon]
+    pairs = {k: pair_metrics(case, k) for k in ks if k - 1 > p.DISCARD}
+    return pairs, p.settled_start({k: v["pass"] for k, v in pairs.items()})
+
+
+def main(argv):
+    """argv: CASE_DIR ... [--reference CASE_DIR]; the reference defaults to the
+    largest D (resolution ladder); domain runs name the 60 D case explicitly."""
+    ref_path = None
+    if "--reference" in argv:
+        i = argv.index("--reference")
+        ref_path, argv = argv[i + 1], argv[:i] + argv[i + 2:]
+    horizon, explicit = p.LAST_SCORED, "--horizon" in argv  # default: original scoring
+    if explicit:
+        i = argv.index("--horizon")
+        horizon, argv = int(argv[i + 1]), argv[:i] + argv[i + 2:]
+    paths = argv if ref_path is None or ref_path in argv else argv + [ref_path]
+    cases = [load(path) for path in paths]
+    for c in cases if explicit else ():  # complete cycles 1..N; never pad or shift
+        if c["cycles"] < horizon or len(c["force"]) < round((0.25 + horizon) * c["T"]):
+            raise SystemExit(f"{c['path']}: run does not contain complete cycles 1..{horizon}")
+    ledgers = [ledger(c, horizon) for c in cases]
+    window = p.common_window([k for _, k in ledgers],
+                             [{k: v["pass"] for k, v in pairs.items()} for pairs, _ in ledgers],
+                             last_scored=horizon)
+    report = dict(window=window, horizon=horizon, cases=[],
+                  horizon_mode="default-cap" if horizon == p.LAST_SCORED else
+                  "extended duration diagnostic (declared amendment, D40/60D only)")
+    for case, (pairs, k_case) in zip(cases, ledgers):
+        limit = round((0.25 + horizon) * case["T"])  # events inside the horizon only
+        entry = {key: case[key] for key in ("path", "D", "T", "box", "n", "u_max", "nu", "tau", "A",
+                                            "steps", "cycles", "F_ref")}
+        entry.update(duration_T=case["steps"] / case["T"], K_case=k_case, pairs=pairs,
+                     coverage={k: coverage(case, k) for k in sorted(cycle_slices(case))
+                               if p.DISCARD < k <= horizon},
+                     no_donor=float(case["diag"][:limit, 13].sum()),
+                     eq_only=float(case["diag"][:limit, 12].sum()),
+                     hf_component=hf_component(case, window or tuple(
+                         [k for k in sorted(cycle_slices(case)) if p.DISCARD < k <= horizon][-3:])))
+        if window:  # residual drift inside the scored window (pairs K..K+2)
+            entry["window_drift"] = {key: max(pairs[k][key] for k in window)
+                                     for key in ("wave", "profile", "morison")}
+            entry["morison_obs"] = cycle_morison(case, window)
+            entry["morison"] = {name: cycle_morison(case, window, fx)
+                                for name, fx in candidates(case).items()}
+            entry["fy_mean"] = float(np.mean(np.concatenate(
+                [case["force"][cycle_slices(case)[k], 1] for k in window])) / case["F_ref"])
+        report["cases"].append(entry)
+    if window and len(cases) > 1:  # cross-resolution waveform against the finest case
+        fine = (next(c for c in cases if c["path"] == ref_path) if ref_path
+                else max(cases, key=lambda c: c["D"]))
+        report["reference"] = fine["path"]
+        ref = np.concatenate([coarse_wave(fine, k) for k in window])
+        for case, entry in zip(cases, report["cases"]):
+            r = fine["T"] // case["T"]  # 1 for equal-D domain runs
+            ref_c = ref.reshape(-1, r).mean(axis=1)
+            own = np.concatenate([coarse_wave(case, k) for k in window])
+            entry["wave_rms_vs_finest"] = float(np.sqrt(np.mean((own - ref_c) ** 2)))
+        for k in window:  # profiles on the common valid point set of all compared cases
+            per_case = [profiles(c, k) for c in cases]
+            keys = common_points(per_case)
+            report.setdefault("common_profile_points", {})[k] = len(keys)
+            if not keys:
+                report["verdict"] = "FAIL-CLOSED: empty common profile set"
+                continue
+            fine_p = per_case[cases.index(fine)]
+            for case, entry, prof in zip(cases, report["cases"], per_case):
+                d = max(max(abs(prof[key][0] - fine_p[key][0]), abs(prof[key][1] - fine_p[key][1]))
+                        for key in keys) / case["u_max"]
+                entry.setdefault("profile_max_vs_finest", {})[k] = float(d)
+        ref_entry = report["cases"][cases.index(fine)]
+        for entry in report["cases"]:  # conservative: difference + both drifts
+            if entry is ref_entry:
+                continue
+            drift = {key: entry["window_drift"][key] + ref_entry["window_drift"][key]
+                     for key in ("wave", "profile", "morison")}
+            dm = max(abs(entry["morison"]["cme"][0] - ref_entry["morison"]["cme"][0]),
+                     abs(entry["morison"]["cme"][1] - ref_entry["morison"]["cme"][1]))
+            entry["conservative"] = dict(
+                wave=entry["wave_rms_vs_finest"] + drift["wave"],
+                profile=max(entry.get("profile_max_vs_finest", {}).values(), default=math.inf)
+                + drift["profile"],
+                morison_cme=dm + drift["morison"])
+            bands = p.DOMAIN if ref_path else p.CROSS  # domain runs name the 60 D reference
+            c = entry["conservative"]
+            entry["pass_scope"] = (
+                ("wave+profile+morison_cme" if "morison" in bands else "wave+profile")
+                + ": conservative (|diff| + both window drifts) on the obs/CME Fx waveform and "
+                "simulation-to-simulation profiles; CME+Wen and CME+node reported, not checked; "
+                "no Fig. 8 / literature verdict")
+            entry["pass"] = (c["wave"] <= bands["wave"] and c["profile"] <= bands["profile"]
+                             and ("morison" not in bands or c["morison_cme"] <= bands["morison"]))
+    print(json.dumps(report, indent=1, default=float))
+
+
+if __name__ == "__main__":
+    main(sys.argv[1:])

research/duetsch_re100_kc5/capture.py

@@ -0,0 +1,52 @@
+"""Research-only field capture helpers for the Dütsch driver (numpy only).
+
+Snapshots are taken at the first step of the scored phase bins in declared
+cycles: a near-body crop at full resolution and a fluid-weighted block-mean
+downsample of the whole domain; plus scalar trends at every bin read.
+"""
+from __future__ import annotations
+
+import numpy as np
+
+CROP_X, CROP_Y, DOWN = 8.0, 6.0, 8      # crop half-widths in D; downsample factor
+STRIP = 2.0                             # inlet/outlet strip width in D
+
+
+def crop_box(xm, ym, D, n, body_extent):
+    """Inclusive integer crop (x0, x1, y0, y1); fail closed if it leaves the domain
+    interior or does not contain the body's swept extent (x_lo, x_hi, y_lo, y_hi)."""
+    x0, x1 = int(round(xm - CROP_X * D)), int(round(xm + CROP_X * D))
+    y0, y1 = int(round(ym - CROP_Y * D)), int(round(ym + CROP_Y * D))
+    if x0 < 1 or y0 < 1 or x1 > n - 2 or y1 > n - 2:
+        raise ValueError("capture crop leaves the domain interior")
+    bx0, bx1, by0, by1 = body_extent
+    if not (x0 <= bx0 and bx1 <= x1 and y0 <= by0 and by1 <= y1):
+        raise ValueError("capture crop clips the body's swept extent")
+    return x0, x1, y0, y1
+
+
+def downsample(field, fluid, f=DOWN):
+    """Fluid-weighted f x f block means; NaN where a block has no fluid cell."""
+    ny, nx = field.shape
+    h, w = ny // f, nx // f
+    v = np.where(fluid, field, 0.0)[:h * f, :w * f].reshape(h, f, w, f).sum(axis=(1, 3))
+    c = fluid[:h * f, :w * f].reshape(h, f, w, f).sum(axis=(1, 3))
+    with np.errstate(invalid="ignore", divide="ignore"):
+        return np.where(c > 0, v / np.maximum(c, 1), np.nan).astype(np.float32)
+
+
+def fluid_mean(field, fluid, xs=slice(None), ys=slice(None)):
+    m = fluid[ys, xs]
+    return float(field[ys, xs][m].mean()) if m.any() else float("nan")
+
+
+def trend_row(s, k, ph, rho, ux, uy, fluid, D, n, crop):
+    """(s, cycle, phase, <rho>, <ux>, <uy>, <rho>_crop, inlet <rho>, <ux>, outlet <rho>, <ux>)."""
+    w = max(1, int(round(STRIP * D)))
+    x0, x1, y0, y1 = crop
+    cx, cy = slice(x0, x1 + 1), slice(y0, y1 + 1)
+    inlet, outlet, inner = slice(1, 1 + w), slice(n - 1 - w, n - 1), slice(1, n - 1)
+    return (s, k, ph, fluid_mean(rho, fluid), fluid_mean(ux, fluid), fluid_mean(uy, fluid),
+            fluid_mean(rho, fluid, cx, cy),
+            fluid_mean(rho, fluid, inlet, inner), fluid_mean(ux, fluid, inlet, inner),
+            fluid_mean(rho, fluid, outlet, inner), fluid_mean(ux, fluid, outlet, inner))

research/duetsch_re100_kc5/protocol.py

@@ -0,0 +1,136 @@
+"""Pure arithmetic of the Dütsch Re 100 / KC 5 protocol (numpy only, no CUDA).
+
+docs/validation/duetsch_re100_kc5_protocol.md is the contract: elapsed time s
+in native steps from the right turning point, paper phase
+phi = (360 s / T - 90) mod 360, paper cycle k >= 1 spans
+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
+
+KC, RE = 5.0, 100.0
+STATIONS = (-0.6, 0.0, 0.6, 1.2)          # lab-fixed x/D from the mean position
+PHASES = (180.0, 210.0, 330.0)            # degrees; bins [p, p + 1)
+X2 = tuple(np.round(np.arange(-1.5, 1.5001, 0.05), 10))  # declared x2/D points
+DISCARD, LAST_SCORED, N_PAIRS = 2, 11, 3   # cycles 0..2 discarded, cap 12 T
+# A-priori bands. Settling is one fifth of the tightest comparison band of the
+# same observable, and comparisons add both cases' residual window drift
+# (conservative: |difference| + drift_a + drift_b <= band).
+DOMAIN = dict(wave=.01, morison=.02, profile=.02)  # vs the 60 D reference at D = 40
+CROSS = dict(wave=.02, profile=.02)                # D40 vs D80
+SETTLE = dict(wave=.002, scalar=.002, profile=.004, morison=.004)
+
+
+def amplitude(D):
+    return KC * D / (2.0 * math.pi)
+
+
+def position(s, T, A):
+    """Centre position relative to the mean: x = A cos(2 pi s / T), evaluated from the
+    native step modulo T so every cycle repeats the same geometry bit for bit."""
+    return A * math.cos(2.0 * math.pi * (int(s) % T) / T)
+
+
+def phase(s, T):
+    """Paper phase in degrees from the native step modulo T (the unwrapped cycle
+    index is kept separately by ``cycle``)."""
+    return (360.0 * (int(s) % T) / T - 90.0) % 360.0
+
+
+def cycle(s, T):
+    return 0 if s < 0.25 * T else int(math.floor(s / T - 0.25)) + 1
+
+
+def bin_phase(s, T):
+    """Scored phase whose 1-degree bin [p, p + 1) contains time s, else None."""
+    phi = phase(s, T)
+    for p in PHASES:
+        if p <= phi < p + 1.0:
+            return p
+    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 = _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)
+    return 0.5 * float(np.sum(w * c * np.abs(c))), 0.5 * float(np.sum(w * np.sin(phi)))
+
+
+def morison_fit(f_norm, b1, b2):
+    """Least squares of f_norm = c_d b1 - (pi^2 c_i / KC) b2 (paper phase, no shift)."""
+    a = np.column_stack([b1, -(math.pi ** 2 / KC) * np.asarray(b2)])
+    (c_d, c_i), *_ = np.linalg.lstsq(a, np.asarray(f_norm), rcond=None)
+    return float(c_d), float(c_i)
+
+
+def inside_circle(x, y, centre, radius):
+    return (x - centre[0]) ** 2 + (y - centre[1]) ** 2 < radius * radius
+
+
+def sample(field, fluid, x, y, centre, radius):
+    """Profile sample: MISSING if the exact point is inside the body (checked before
+    any corner interpolation), else ``bilinear_fluid``."""
+    if inside_circle(x, y, centre, radius):
+        return None, (), 0.0
+    return bilinear_fluid(field, fluid, x, y)
+
+
+def bilinear_fluid(field, fluid, x, y):
+    """Bilinear value at (x, y) in cell coordinates from fluid cells only.
+
+    Returns (value, ((ix, iy, weight), ...), fluid_weight_sum) with weights
+    renormalized over the fluid corners and the original sum of their bilinear
+    weights, or (None, (), 0.0) when no corner is fluid (MISSING, never zero).
+    """
+    x0, y0 = int(math.floor(x)), int(math.floor(y))
+    fx, fy = x - x0, y - y0
+    corners = ((x0, y0, (1 - fx) * (1 - fy)), (x0 + 1, y0, fx * (1 - fy)),
+               (x0, y0 + 1, (1 - fx) * fy), (x0 + 1, y0 + 1, fx * fy))
+    used = [(i, j, w) for i, j, w in corners if w > 0.0 and fluid[j, i]]
+    total = sum(w for _, _, w in used)
+    if total == 0.0:
+        return None, (), 0.0
+    used = tuple((i, j, w / total) for i, j, w in used)
+    return float(sum(w * field[j, i] for i, j, w in used)), used, float(total)
+
+
+def settled_start(pairs_pass):
+    """K_case from {k: pass of pair (k-1 -> k)}: first cycle after three consecutive
+    passing pairs among cycles > DISCARD, or None."""
+    run = 0
+    for k in sorted(pairs_pass):
+        if k - 1 <= DISCARD:
+            continue
+        run = run + 1 if pairs_pass[k] else 0
+        if run == N_PAIRS:
+            return k + 1
+    return None
+
+
+def common_window(k_cases, pairs_pass_by_case, last_scored=LAST_SCORED):
+    """Common scored cycles K..K+2 with K = max K_case, or None (FAIL-CLOSED) when a
+    case never settled, K + 2 exceeds the cap, or a case fails any pair p_j for
+    j = its own K_case .. K + 2 (continuous settling, including the wait for K)."""
+    if any(k is None for k in k_cases):
+        return None
+    K = max(k_cases)
+    if K + 2 > last_scored:
+        return None
+    for k_case, pairs in zip(k_cases, pairs_pass_by_case):
+        if not all(pairs.get(j, False) for j in range(k_case, K + 3)):
+            return None
+    return K, K + 1, K + 2

research/duetsch_re100_kc5/run.py

@@ -0,0 +1,133 @@
+"""One-shot Dütsch Re 100 / KC 5 prescribed in-line oscillation driver.
+
+Protocol: docs/validation/duetsch_re100_kc5_protocol.md. Native dense channel
+(zero-velocity inlet, pressure-NEQ outlet, bounce-back y walls), MRT, lab frame.
+Writes raw per-step force/diagnostics and phase-bin profile samples to --out
+(outside the checkout); analysis is in analyze.py. No production change.
+
+    python run.py --D 40 --box 40 --cycles 11 --out /path/outside/checkout
+"""
+from __future__ import annotations
+
+import argparse
+import json
+import math
+import os
+import tempfile
+
+import numpy as np
+
+import protocol as p
+
+
+TEMPLATE = os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "..", "src",
+                        "CelerisLab", "configs", "config_lbm.json")
+
+
+def config(D, box, T):
+    """Shipped config_lbm.json (full parser schema) with only the case fields
+    overridden: grid, still fluid (velocity 0), nu from Re, MRT, uniform inlet
+    profile, bounce-back y walls. Inlet scheme and outlet mode stay the
+    template's (zou_he_local, neq_extrap) and are recorded in case.json."""
+    u = p.KC * D / T
+    n = int(round(box * D))
+    with open(TEMPLATE) as fh:
+        cfg = json.load(fh)
+    cfg["grid"].update(lattice_model="D2Q9", nx=n, ny=n, nz=1)
+    cfg["physics"].update(data_type="FP32", viscosity=u * D / p.RE, velocity=0.0, rho=1.0)
+    cfg["method"].update(collision="MRT", streaming="double_buffer", store_precision="FP32",
+                         ddf_shifting=False, y_wall_bc="bounce_back")
+    cfg["method"]["les"]["enabled"] = False
+    cfg["method"]["inlet"]["profile"] = "uniform"
+    return cfg, u, n
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--D", type=int, required=True)
+    ap.add_argument("--box", type=float, default=40.0, help="box side in D")
+    ap.add_argument("--T", type=int, default=None, help="period in steps (multiple of 360)")
+    ap.add_argument("--cycles", type=int, default=p.LAST_SCORED, help="run to the end of this cycle")
+    ap.add_argument("--out", required=True)
+    ap.add_argument("--dump-fields", default="",
+                    help="research capture: comma-separated paper cycles, e.g. 11,21,30,31")
+    args = ap.parse_args()
+    D, T = args.D, args.T or 126 * args.D
+    if T % 360:
+        raise SystemExit("T must be a multiple of 360 steps")
+    if os.path.realpath(args.out).startswith(os.path.realpath(os.path.join(os.path.dirname(__file__), "..", ".."))):
+        raise SystemExit("--out must be outside the checkout")
+    cfg, u_max, n = config(D, args.box, T)
+    A, R = p.amplitude(D), 0.5 * D
+    xm = ym = 0.5 * (n - 1)
+    os.makedirs(args.out, exist_ok=True)
+
+    from CelerisLab.body.motion import Translation
+    from CelerisLab.simulation import Simulation
+    with tempfile.NamedTemporaryFile("w", suffix=".json", delete=False) as fh:
+        json.dump(cfg, fh)
+    sim = Simulation(fh.name)
+    try:
+        body = sim.add_body("circle", (xm + A, ym), R,
+                            translation=Translation((-2.0 * A, 0.0, 0.0, 0.0), 1.01 * u_max))
+        sim.initialize()
+        dump = {int(c) for c in args.dump_fields.split(",") if c}
+        if dump:
+            import capture
+            if max(dump) > args.cycles or min(dump) <= p.DISCARD:
+                raise SystemExit("--dump-fields cycles must lie in the scored range of this run")
+            crop = capture.crop_box(xm, ym, D, n, (xm - A - R, xm + A + R, ym - R, ym + R))
+            trend = []
+        steps = int(round((0.25 + args.cycles) * T))
+        force = np.zeros((steps, 2)); diag = []; samples = []
+        for m in range(steps):
+            x0, x1 = p.position(m, T, A), p.position(m + 1, T, A)
+            sim.set_body(body, offset=(max(x1 - A, -2.0 * A), 0.0), vx=x1 - x0)
+            sim.run(1)
+            force[m] = sim.read_force(body)[:2]
+            diag.append(sim.moving_diagnostics()[0])
+            s = m + 1  # state after the step
+            ph = p.bin_phase(s, T)
+            if ph is not None and p.cycle(s, T) > p.DISCARD:
+                fields = sim.get_macroscopic()
+                fluid = (sim.get_flags().reshape(n, n) & 1) != 0
+                ux, uy = (np.asarray(fields[k]).reshape(n, n) for k in ("ux", "uy"))
+                if dump:  # research capture only; the default path is unchanged
+                    rho = np.asarray(fields["rho"]).reshape(n, n)
+                    trend.append(capture.trend_row(s, p.cycle(s, T), ph, rho, ux, uy, fluid, D, n, crop))
+                    if p.cycle(s, T) in dump and p.phase(s, T) == ph:  # first step of the bin
+                        x0c, x1c, y0c, y1c = crop
+                        cy_, cx_ = slice(y0c, y1c + 1), slice(x0c, x1c + 1)
+                        np.savez_compressed(
+                            os.path.join(args.out, f"fields_c{p.cycle(s, T)}_p{int(ph)}.npz"),
+                            s=s, cycle=p.cycle(s, T), phase=ph, crop=np.asarray(crop),
+                            force=force[m], centre=np.asarray((xm + p.position(s, T, A), ym)),
+                            ux=ux[cy_, cx_].astype(np.float32), uy=uy[cy_, cx_].astype(np.float32),
+                            rho=rho[cy_, cx_].astype(np.float32), fluid=fluid[cy_, cx_].astype(np.uint8),
+                            down_ux=capture.downsample(ux, fluid), down_uy=capture.downsample(uy, fluid),
+                            down_rho=capture.downsample(rho, fluid))
+                centre = (xm + p.position(s, T, A), ym)  # analytic centre after the step
+                for xs in p.STATIONS:
+                    for x2 in p.X2:
+                        x, y = xm + xs * D, ym + x2 * D
+                        u, wu, wsum = p.sample(ux, fluid, x, y, centre, R)
+                        v, _, _ = p.sample(uy, fluid, x, y, centre, R)
+                        samples.append((s, p.cycle(s, T), ph, xs, x2,
+                                        math.nan if u is None else u, math.nan if v is None else v,
+                                        wsum, json.dumps(wu)))
+        if dump:
+            np.save(os.path.join(args.out, "trend.npy"), np.asarray(trend, float))
+        np.savez(os.path.join(args.out, "raw.npz"), force=force, diag=np.asarray(diag),
+                 samples=np.asarray([r[:8] for r in samples], float))
+        with open(os.path.join(args.out, "case.json"), "w") as fh:
+            json.dump(dict(D=D, T=T, box=args.box, n=n, u_max=u_max, nu=cfg["physics"]["viscosity"],
+                           tau=3 * cfg["physics"]["viscosity"] + 0.5, A=A, steps=steps,
+                           cycles=args.cycles, F_ref=0.5 * u_max ** 2 * D,
+                           inlet=cfg["method"]["inlet"], outlet=cfg["method"]["outlet"],
+                           weights=[r[8] for r in samples]), fh)
+    finally:
+        sim.close()
+
+
+if __name__ == "__main__":
+    main()

research/steady_pinball/config.json

@@ -0,0 +1,49 @@
+{
+    "_doc": "Legacy geometry Re_D=50, U0=0.01, SI=800 baseline (kar_legacy_re50d_u001)",
+    "grid": {
+        "lattice_model": "D2Q9",
+        "nx": 1280,
+        "ny": 512,
+        "nz": 1
+    },
+    "physics": {
+        "data_type": "FP32",
+        "viscosity": 0.004,
+        "velocity": 0.01,
+        "rho": 1.0
+    },
+    "method": {
+        "collision": "MRT",
+        "streaming": "double_buffer",
+        "store_precision": "FP32",
+        "ddf_shifting": false,
+        "les": {
+            "enabled": false,
+            "cs": 0.16,
+            "closed_form": true
+        },
+        "trt": {
+            "magic_param": 0.1875
+        },
+        "inlet": {
+            "profile": "parabolic",
+            "scheme": "regularized",
+            "regularized_neq_damp": 1.0
+        },
+        "outlet": {
+            "mode": "neq_extrap",
+            "backflow_clamp": true,
+            "blend_alpha": 0.7,
+            "srt_neq_damp": 0.5
+        },
+        "y_wall_bc": "bounce_back",
+        "omega_guard": {
+            "min": 0.01,
+            "max": 1.99
+        }
+    },
+    "cuda": {
+        "threads_per_block": 256,
+        "compute_capability": "auto"
+    }
+}

research/steady_pinball/run.py

@@ -0,0 +1,167 @@
+"""Uncontrolled paper-geometry steady Re_D=50, one transverse assembly DOF.
+
+Three predeclared zero-damping structures; --smoke uses only 32 + 32 steps.
+See docs/steady_pinball_fiv.md for historical solver differences and new structure.
+"""
+from __future__ import annotations
+
+import argparse
+import json
+import math
+from pathlib import Path
+import subprocess
+import time
+
+import numpy as np
+
+D, U, NU = 20., .01, .004
+CENTRES = ((600., 255.5), (626., 270.5), (626., 240.5))
+PROBES = ((800., 295.5), (800., 255.5), (800., 215.5))
+WARMUP, STEPS, PERIOD = 512000, 600000, 10000
+STRUCTURES = {'Z0': (10., 1.e-4), 'M2': (2., 1.e-4), 'R2': (2., 6.e-5)}
+COLUMNS = ('step', 'y', 'v', 'cme_F', 'cme_U', 'cme_L', 'wen_F', 'wen_U', 'wen_L',
+           'node_F', 'node_U', 'node_L', 'eq_only', 'no_donor', 'crossings')
+ROOT = Path(__file__).resolve().parents[2]
+
+
+def config():
+    return json.loads(Path(__file__).with_name('config.json').read_text())
+
+
+class Structure:
+    """Exact damped oscillator update for a force held over dt=1; no force filter."""
+    def __init__(self, mass_ratio=10., frequency=1.e-4, zeta=0.):
+        if not (math.isfinite(mass_ratio) and mass_ratio > 0. and
+                math.isfinite(frequency) and frequency > 0. and
+                math.isfinite(zeta) and 0. <= zeta < 1.):
+            raise ValueError('positive finite mass ratio/frequency and 0 <= zeta < 1 required')
+        self.mass_ratio, self.frequency = mass_ratio, frequency
+        self.mass = mass_ratio * 3. * math.pi * D * D / 4.
+        self.omega = 2. * math.pi * frequency
+        self.zeta = zeta
+        self.stiffness = self.mass * self.omega**2
+        self.damping = 2. * self.zeta * self.mass * self.omega
+        self.y = self.v = 0.
+        a = self.zeta * self.omega
+        wd = self.omega * math.sqrt(1. - self.zeta**2)
+        self.ec = math.exp(-a) * math.cos(wd)
+        self.es = math.exp(-a) * math.sin(wd) / wd
+        self.a = a
+
+    def candidate(self, force):
+        equilibrium = force / self.stiffness
+        z = self.y - equilibrium
+        return (equilibrium + (self.ec + self.a*self.es)*z + self.es*self.v,
+                -self.omega**2*self.es*z + (self.ec - self.a*self.es)*self.v)
+
+
+def summarize(rows):
+    windows = []
+    for start in range(0, len(rows), 10*PERIOD):
+        block = rows[start:start + 10*PERIOD]
+        fy = block[:, 3:9].sum(axis=1)
+        windows.append(dict(start_step=int(block[0, 0]), end_step=int(block[-1, 0]),
+                            y_mean_D=float(block[:, 1].mean()/D),
+                            y_rms_D=float(np.sqrt(np.mean(block[:, 1]**2))/D),
+                            y_peak_D=float(np.abs(block[:, 1]).max()/D),
+                            fy_mean=float(fy.mean()), fy_rms=float(np.sqrt(np.mean(fy**2))),
+                            mean_power=float(np.mean(fy*block[:, 2]))))
+    return windows
+
+
+def main():
+    ap = argparse.ArgumentParser(description=__doc__)
+    ap.add_argument('--out', required=True, type=Path)
+    ap.add_argument('--fixed', action='store_true', help='same modern moving-path geometry held fixed')
+    ap.add_argument('--smoke', action='store_true', help='32+32 steps, no physical result')
+    ap.add_argument('--case', choices=STRUCTURES, default='Z0',
+                    help='predeclared structure; all cases have zero structural damping')
+    args = ap.parse_args()
+    out = args.out.resolve()
+    if out == ROOT or ROOT in out.parents:
+        ap.error('output must be outside the checkout')
+    out.mkdir(parents=True, exist_ok=False)
+    cfg = config()
+    (out/'config.json').write_text(json.dumps(cfg, indent=2)+'\n')
+    source = subprocess.check_output(['git', '-C', str(ROOT), 'rev-parse', 'HEAD'], text=True).strip()
+    s = Structure(*STRUCTURES[args.case], zeta=0.)
+    warmup, count = (32, 32) if args.smoke else (WARMUP, STEPS)
+    rows = np.lib.format.open_memmap(out/'steps.npy', mode='w+', dtype=np.float64,
+                                     shape=(count, len(COLUMNS)))
+    rows[:] = np.nan
+    from CelerisLab import Simulation
+    from CelerisLab.body.motion import MovingTables, Translation
+    names = MovingTables.DIAG
+    col = {k: names.index(k) for k in names}
+    status, completed, f_last = 'running', 0, 0.
+    started = time.monotonic()
+    report = dict(source=source, case=args.case, fixed=args.fixed, smoke=args.smoke, columns=COLUMNS,
+                  structure=dict(mass=s.mass, stiffness=s.stiffness, damping=s.damping,
+                                 mass_ratio_three_disks=s.mass_ratio,
+                                 zeta=s.zeta, fn=s.frequency, Tn=1./s.frequency, y0=0., v0=0.),
+                  warmup_steps=warmup, requested_steps=count, dt=1., D=D, U=U, nu=NU,
+                  duration_Tn=count*s.frequency, duration_tU_D=count*U/D,
+                  summary_block_steps=10*PERIOD,
+                  rhs='sum of three instantaneous CME_y + Wen_y, one-step lag, no filtering')
+    try:
+        with Simulation(str(out/'config.json')) as sim:
+            # Match Legacy object order; probes remain passive and fixed in the lab.
+            for c in PROBES:
+                sim.add_body('sensor', c, 5.)
+            envelope = Translation((0., 0., -D, D), .05, shared=True)
+            bodies = [sim.add_body('circle', c, D/2, translation=envelope) for c in CENTRES]
+            sim.initialize()
+            # At release obtain the last instantaneous force, not a warmup average.
+            remaining = warmup - 1
+            while remaining:
+                chunk = min(10000, remaining)
+                sim.run(chunk)
+                remaining -= chunk
+            sim.run(1)
+            d = sim.moving_diagnostics()
+            f_last = float(d[:, col['cme_y']].sum(dtype=np.float64) + d[:, col['wen_y']].sum(dtype=np.float64))
+            report['warmup_elapsed_s'] = time.monotonic()-started
+            print(json.dumps(dict(event='release', **report)), flush=True)
+            for n in range(count):
+                y1, v1 = (0., 0.) if args.fixed else s.candidate(f_last)
+                if not all(map(math.isfinite, (y1, v1))) or abs(y1) > D:
+                    status = 'stopped_nonfinite_or_travel_bound'
+                    break
+                if not args.fixed:
+                    for bid in bodies:
+                        sim.set_body(bid, offset=(0., y1), vy=y1-s.y)
+                sim.run(1)
+                d = sim.moving_diagnostics()
+                cme, wen = d[:, col['cme_y']], d[:, col['wen_y']]
+                node = d[:, col['cover_jy']] - d[:, col['uncover_jy']]
+                row = (n+1, y1, v1, *cme, *wen, *node,
+                       d[:, col['eq_only']].sum(), d[:, col['no_donor']].sum(), d[:, col['crossings']].sum())
+                if not np.isfinite(row).all() or row[-2] != 0:
+                    status = 'stopped_nonfinite_or_no_donor'
+                    break
+                rows[n] = row
+                completed = n+1
+                s.y, s.v = y1, v1
+                f_last = float(cme.sum(dtype=np.float64)+wen.sum(dtype=np.float64))
+                if completed % PERIOD == 0:
+                    rows.flush()
+                    print(json.dumps(dict(event='progress', step=completed, t_Tn=completed*s.frequency,
+                                          y_D=s.y/D, elapsed_s=time.monotonic()-started)), flush=True)
+            else:
+                status = 'smoke_completed' if args.smoke else 'duration_completed_not_qualified'
+    except BaseException as exc:
+        status = f'error: {type(exc).__name__}: {exc}'
+        raise
+    finally:
+        rows.flush()
+        report.update(status=status, completed_steps=completed, elapsed_s=time.monotonic()-started,
+                      windows=summarize(rows[:completed]),
+                      totals={key: float(rows[:completed, i].sum()) for key, i in
+                              [('eq_only', 12), ('no_donor', 13), ('crossings', 14)]},
+                      limitation='Uncontrolled development response only. No automatic settling, grid convergence, control or historical-solver equivalence claim.')
+        (out/'report.json').write_text(json.dumps(report, indent=2)+'\n')
+        print(json.dumps(report), flush=True)
+
+
+if __name__ == '__main__':
+    main()

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

@@ -0,0 +1,136 @@
+"""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 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
+"""
+from __future__ import annotations
+
+import argparse
+import json
+import math
+import os
+import tempfile
+
+import numpy as np
+
+import viv
+
+D, U = 40, 0.04
+TEMPLATE = os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "..", "src",
+                        "CelerisLab", "configs", "config_lbm.json")
+
+
+def config(D=D, U=U):
+    nx, ny = 75 * D, 60 * D
+    with open(TEMPLATE) as fh:
+        cfg = json.load(fh)
+    cfg["grid"].update(lattice_model="D2Q9", nx=nx, ny=ny, nz=1)
+    cfg["physics"].update(data_type="FP32", viscosity=U * D / viv.RE, velocity=U, rho=1.0)
+    cfg["method"].update(collision="MRT", streaming="double_buffer", store_precision="FP32",
+                         ddf_shifting=False, y_wall_bc="free_slip")
+    cfg["method"]["les"]["enabled"] = False
+    cfg["method"]["inlet"]["profile"] = "uniform"
+    return cfg, nx, ny
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--out", required=True)
+    ap.add_argument("--D", type=int, default=D)
+    ap.add_argument("--smoke", action="store_true",
+                    help="plumbing smoke only: exactly --steps steps (<= 64), never scored")
+    ap.add_argument("--steps", type=int, default=0)
+    args = ap.parse_args()
+    if args.smoke != bool(args.steps) or args.steps > 64 or args.steps < 0:
+        raise SystemExit("--smoke requires --steps N with 1 <= N <= 64 (and only with --smoke)")
+    if args.D != D and not args.smoke:
+        raise SystemExit("the physical case is predeclared at D = 40")
+    if os.path.realpath(args.out).startswith(os.path.realpath(os.path.join(os.path.dirname(__file__), "..", ".."))):
+        raise SystemExit("--out must be outside the checkout")
+    os.makedirs(args.out, exist_ok=True)
+    for name in ("report.json", "steps.npy"):  # never overwrite earlier evidence
+        if os.path.exists(os.path.join(args.out, name)):
+            raise SystemExit(f"{args.out} already holds {name}; use a fresh output directory")
+    Dn = args.D
+    cfg, nx, ny = config(Dn)
+    s, cyc = viv.Structure(Dn, U), viv.Cycles(Dn)
+    T_n = 1.0 / s.f_n
+    step_cap, n_seed = int(round(viv.STEP_CAP_TN * T_n)), int(round(viv.SEED_TN * T_n))
+    quiet = int(round(5 * T_n))
+    if args.smoke:
+        step_cap = args.steps
+    centre = (30.0 * Dn, 0.5 * (ny - 1))
+
+    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)
+    sim = Simulation(fh.name)
+    rows, status, seed_step, f_last = [], "running", None, 0.0
+    try:
+        body = sim.add_body("circle", centre, 0.5 * Dn,
+                            translation=Translation((0.0, 0.0, -1.5 * Dn, 1.5 * Dn), 0.5))
+        sim.initialize()
+        for n in range(step_cap):
+            # rows[k] holds the state y_{k+1}; the declared window is y_t, t in [n_seed - 5 T_n, n_seed)
+            if not args.smoke and n == n_seed and \
+                    max(abs(r[1]) for r in rows[n_seed - quiet - 1:n_seed - 1]) < viv.ELIGIBLE * Dn:
+                s.v += s.seed_velocity(Dn)         # after the quiet check, before the update
+                seed_step = n
+                cyc.discard_current()
+                cyc.reset_proof()
+            y1, v1 = s.candidate(f_last)
+            if not (math.isfinite(y1) and math.isfinite(v1)) or abs(y1) > viv.Y_LIMIT * Dn:
+                status = "fail_bound_or_nonfinite"
+                break
+            try:
+                sim.set_body(body, offset=(0.0, y1), vy=y1 - s.y)
+            except ValueError as exc:              # input rejection: stop, no retry
+                status = f"rejected: {exc}"
+                break
+            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[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
+            s.y, s.v, f_last = y1, v1, f_rhs       # commit only after a successful step
+            rows.append(row)
+            if args.smoke:
+                continue                           # no settle / score in the plumbing smoke
+            if cyc.update(s.y) == "scored":
+                status = "scored"
+                break
+            if cyc.counted >= viv.CYCLE_CAP:
+                status = "timeout_response_cycles"
+                break
+        else:
+            status = ("smoke_ok" if all(r[8] == 0.0 for r in rows) else "smoke_no_donor") \
+                if args.smoke else "timeout_native_steps"
+    except BaseException as exc:
+        status = f"error: {exc!r}"
+        raise
+    finally:
+        np.save(os.path.join(args.out, "steps.npy"), np.asarray(rows, float))
+        report = dict(status=status, smoke=args.smoke, D=Dn, steps=len(rows), t_U_over_D=len(rows) * U / Dn, T_n=T_n,
+                      seed_step=seed_step, counted_cycles=cyc.counted, cycles=cyc.log,
+                      columns=["n", "y", "v", "F_rhs_cme_wen", "F_cme", "F_wen", "F_obs",
+                               "node_cover_minus_uncover_y", "no_donor"])
+        if status == "scored" and not args.smoke:
+            report["score"] = viv.score(cyc.window)
+        with open(os.path.join(args.out, "report.json"), "w") as fh:
+            json.dump(report, fh, indent=1)
+        sim.close()
+
+
+if __name__ == "__main__":
+    main()

research/viv_re150/viv.py

@@ -0,0 +1,126 @@
+"""Pure logic of the predeclared Re 150 free transverse VIV case (stdlib only).
+
+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.
+"""
+from __future__ import annotations
+
+import math
+
+RE, U_STAR, MASS_NUMBER = 150.0, 5.0, 2.0      # m/(rho D^2) = 2 -> m* = 8/pi
+TARGET_A, BAND = 0.520, 0.10                   # Qin §4.5 maximum amplitude, ±10 %
+ELIGIBLE, BLOCK, LAST5 = 0.01, 10, 5           # in D; cycles; cycles
+WINDOW_TOL, LAST5_TOL = 0.02, 0.01
+CYCLE_CAP, STEP_CAP_TN, SEED_TN, Y_LIMIT = 150, 200, 40, 1.0
+
+
+def finite(values):
+    """True iff every recorded value is finite (fail-closed output check)."""
+    return all(math.isfinite(float(v)) for v in values)
+
+
+class Structure:
+    """Undamped spring-mass, exact harmonic update over dt = 1 with constant force."""
+
+    def __init__(self, D, u):
+        self.m = MASS_NUMBER * D * D           # lattice rho = 1
+        self.f_n = u / (U_STAR * D)
+        self.w = 2.0 * math.pi * self.f_n
+        self.k = self.m * self.w * self.w
+        self.y = self.v = 0.0
+
+    def candidate(self, force):
+        ye = force / self.k
+        c, s = math.cos(self.w), math.sin(self.w)
+        return (ye + (self.y - ye) * c + self.v / self.w * s,
+                -(self.y - ye) * self.w * s + self.v * c)
+
+    def seed_velocity(self, D):
+        return self.w * ELIGIBLE * D           # amplitude ~ .01 D
+
+
+class Cycles:
+    """Upward y = 0 crossings, eligibility, blocks, settle/score with restart."""
+
+    def __init__(self, D):
+        self.D = D
+        self.reset_proof()
+        self.counted = 0
+        self.log = []                          # one dict per completed cycle
+        self._cur = None
+        self._prev = 0.0
+        self._since = 0
+
+    def reset_proof(self):
+        self.amps = []                         # eligible amplitudes since proof start
+        self.blocks = []                       # block means
+        self.settled_at = None                 # block index where R held
+        self.window = None                     # completed scoring cycles (list of dicts)
+
+    def discard_current(self):
+        self._cur, self._since = None, 0
+
+    def rule(self, k):
+        """R(k): block-mean change ≤ 2 % and last-five spread ≤ 1 %."""
+        if k < 1:
+            return False
+        w0, w1 = self.blocks[k - 1], self.blocks[k]
+        last = self.amps[(k + 1) * BLOCK - LAST5:(k + 1) * BLOCK]
+        a5 = sum(last) / LAST5
+        return abs(w1 - w0) / w0 <= WINDOW_TOL and max(abs(a - a5) for a in last) / a5 <= LAST5_TOL
+
+    def update(self, y):
+        """Feed y after a committed step; returns 'scored' when a stable window completes."""
+        out = None
+        crossing = self._prev < 0.0 <= y and self._since >= 2
+        if crossing:
+            if self._cur is not None:
+                out = self._close()
+            self._cur, self._since = dict(ymax=y, ymin=y, ysum=0.0, n=0), 0
+        if self._cur is not None:
+            c = self._cur
+            c["ymax"], c["ymin"] = max(c["ymax"], y), min(c["ymin"], y)
+            c["ysum"] += y
+            c["n"] += 1
+        self._since += 1
+        self._prev = y
+        return out
+
+    def _close(self):
+        c = self._cur
+        a = 0.5 * (c["ymax"] - c["ymin"]) / self.D
+        cyc = dict(ymax=c["ymax"] / self.D, ymin=c["ymin"] / self.D,
+                   offset=c["ysum"] / c["n"] / self.D, amp=a, steps=c["n"], eligible=a >= ELIGIBLE)
+        self.log.append(cyc)
+        if not cyc["eligible"]:
+            return None
+        self.counted += 1
+        self.amps.append(a)
+        if self.settled_at is not None:
+            self.window.append(cyc)
+        if len(self.amps) % BLOCK:
+            return None
+        self.blocks.append(sum(self.amps[-BLOCK:]) / BLOCK)
+        k = len(self.blocks) - 1
+        if self.settled_at is None:
+            if self.rule(k):
+                self.settled_at, self.window = k, []
+            return None
+        if not self.rule(k):                   # scoring block failed: restart from it
+            amps, blocks = self.amps[-BLOCK:], self.blocks[-1:]
+            self.reset_proof()
+            self.amps, self.blocks = amps, blocks
+            return None
+        return "scored" if k == self.settled_at + 2 else None
+
+
+def score(window):
+    """Target metric and reported quantities over the 20 scoring cycles."""
+    a_max = max(max(abs(c["ymax"]), abs(c["ymin"])) for c in window)
+    return dict(A_max=a_max, mean_amp=sum(c["amp"] for c in window) / len(window),
+                ymax=max(c["ymax"] for c in window), ymin=min(c["ymin"] for c in window),
+                offset=sum(c["offset"] for c in window) / len(window),
+                T_resp=sum(c["steps"] for c in window) / len(window),
+                within_band=abs(a_max - TARGET_A) <= BAND * TARGET_A)

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
 
@@ -14,14 +14,11 @@ from typing import Callable, Optional, Tuple
 
 MAX_STEP = 0.5     # Euclidean centre displacement per step (lattice cells)
 BAND_MARGIN = 2    # cells added around the swept circle AABB
+SHARED_CLEARANCE = 4.0  # > two diagonal links + MAX_STEP; docs/steady_pinball_fiv.md
 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)
@@ -35,6 +32,7 @@ class Translation:
     """Declared motion envelope: offset box (x_min, x_max, y_min, y_max) and speed bound."""
     offset_box: Tuple[float, float, float, float]
     max_speed: float
+    shared: bool = False  # same envelope instance: synchronous transverse motion
 
     def __post_init__(self):
         box = tuple(float(v) for v in self.offset_box)
@@ -45,10 +43,29 @@ class Translation:
             raise ValueError("offset_box must contain the zero offset")
         if not 0.0 < speed <= MAX_STEP:
             raise ValueError(f"max_speed must satisfy 0 < max_speed <= {MAX_STEP}")
+        if self.shared and box[:2] != (0.0, 0.0):
+            raise ValueError("shared translation is transverse only (x offsets must be zero)")
         object.__setattr__(self, "offset_box", box)
         object.__setattr__(self, "max_speed", speed)
 
 
+def validate_shared_translations(circles):
+    """Check common old/target offsets and wall velocities before device work.
+
+    Per-body setters may stage members in any order. Envelope identity (not
+    value equality) defines a group; a missing target means hold at applied.
+    """
+    groups = {}
+    for motion in circles:
+        if not motion.translation.shared:
+            continue
+        target = motion.pending if motion.pending is not None else (motion.applied, (0.0, 0.0))
+        state = (motion.applied, target)
+        key = id(motion.translation)
+        if groups.setdefault(key, state) != state:
+            raise ValueError("shared translation needs identical applied/target offsets and velocities")
+
+
 class MovingCircle:
     """Applied offset, pending target and the fail-closed commit of one circle."""
 
@@ -76,6 +93,8 @@ class MovingCircle:
         x0, x1, y0, y1 = self.translation.offset_box
         if not (x0 <= offset[0] <= x1 and y0 <= offset[1] <= y1):
             raise ValueError(f"offset {offset} outside declared box {self.translation.offset_box}")
+        if self.translation.shared and velocity[0] != 0.0:
+            raise ValueError("shared translation needs zero x velocity")
         step = math.hypot(offset[0] - self.applied[0], offset[1] - self.applied[1])
         if step > MAX_STEP:
             raise ValueError(f"per-step displacement {step} exceeds {MAX_STEP}")
@@ -98,23 +117,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 +140,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()
 
@@ -257,6 +259,8 @@ class Simulation:
         moving = self._moving_bodies()
         if not moving:
             return
+        from .body.motion import SHARED_CLEARANCE, validate_shared_translations
+        validate_shared_translations([obj.motion for obj in moving])
         cfg = self.lbm_cfg
         if (not cfg.is_d2q9 or cfg.streaming != "double_buffer" or cfg.ddf_shifting
                 or cfg.store_precision != "FP32"):
@@ -275,6 +279,13 @@ class Simulation:
                 if other is obj:
                     continue
                 (cx, cy), r = other.center[:2], float(other.radius)
+                if (other.motion is not None and obj.motion.translation.shared
+                        and obj.motion.translation is other.motion.translation):
+                    # Whole envelopes may overlap; simultaneous local stencils may not.
+                    gap = np.hypot(obj.center[0] - cx, obj.center[1] - cy) - obj.radius - r
+                    if gap <= SHARED_CLEARANCE:
+                        raise ValueError("shared translation needs surface clearance > 4 cells")
+                    continue
                 if other.motion is not None:
                     ox0, ox1, oy0, oy1 = other.motion.band()
                 else:
@@ -360,8 +371,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
@@ -714,6 +728,9 @@ class Simulation:
         if steps < 0 or checkpoint_interval < 0:
             raise ValueError("steps and checkpoint_interval must be nonnegative")
         moving = self._moving_bodies()
+        if steps > 0 and moving:
+            from .body.motion import validate_shared_translations
+            validate_shared_translations([obj.motion for obj in moving])
         if steps > 0 and moving and self.stepper.moving is None:
             raise RuntimeError(MOVING_BLOCKED)
         if moving and checkpoint_interval > 0:
@@ -781,6 +798,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 +811,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/integration/test_shared_translation_gpu.py

@@ -0,0 +1,62 @@
+"""Small serial GPU gates for shared transverse circles; no physical benchmark."""
+import json
+from pathlib import Path
+
+import numpy as np
+
+from CelerisLab import Simulation
+from CelerisLab.body.motion import MovingTables, Translation
+
+CENTRES = ((70., 79.5), (96., 94.5), (96., 64.5))
+
+
+def simulate(tmp_path, *, moving, offsets, reverse=False):
+    config = json.loads((Path(__file__).resolve().parents[2]/'research/steady_pinball/config.json').read_text())
+    config['grid'].update(nx=192, ny=160)
+    p = tmp_path/'config.json'
+    p.write_text(json.dumps(config))
+    with Simulation(str(p)) as sim:
+        envelope = Translation((0., 0., -20., 20.), .5, shared=True)
+        centres = list(reversed(CENTRES)) if reverse else list(CENTRES)
+        ids = [sim.add_body('circle', c, 10., translation=envelope if moving else None) for c in centres]
+        sim.initialize()
+        previous = 0.
+        crossings = 0.
+        for y in offsets:
+            if moving:
+                for bid in ids:
+                    sim.set_body(bid, offset=(0., float(y)), vy=float(y-previous))
+            sim.run(1)
+            if moving:
+                diag = sim.moving_diagnostics()
+                assert np.isfinite(diag).all()
+                assert np.all(diag[:, MovingTables.DIAG.index('no_donor')] == 0)
+                crossings += diag[:, MovingTables.DIAG.index('crossings')].sum()
+            previous = y
+        forces = np.array([sim.read_force(bid) for bid in ids])
+        if reverse:
+            forces = forces[::-1]
+        return sim.get_ddf(), sim.get_flags(), forces, crossings
+
+
+def test_fixed_vs_shared_hold(tmp_path):
+    fixed = simulate(tmp_path, moving=False, offsets=[0.]*32)
+    held = simulate(tmp_path, moving=True, offsets=[0.]*32)
+    assert np.array_equal(fixed[0], held[0], equal_nan=True)
+    assert np.array_equal(fixed[1], held[1])
+    # FP32 atomic force reduction order only; populations remain exactly equal.
+    np.testing.assert_allclose(fixed[2], held[2], rtol=0., atol=1e-4)
+    assert held[3] == 0
+
+
+def test_common_motion_member_order_has_no_population_effect(tmp_path):
+    # MAX_STEP=.5 is a geometric bound, not a low-Mach stability guarantee.
+    # Keep the complete +/-D excursion at wall speed .01, equal to mean inflow.
+    offsets = np.concatenate((np.arange(1, 2001), np.arange(1999, -2001, -1),
+                              np.arange(-1999, 1))) * .01
+    forward = simulate(tmp_path, moving=True, offsets=offsets)
+    reverse = simulate(tmp_path, moving=True, offsets=offsets, reverse=True)
+    assert forward[3] > 0 and forward[3] == reverse[3]
+    assert np.array_equal(forward[0], reverse[0], equal_nan=True)
+    assert np.array_equal(forward[1], reverse[1])
+    np.testing.assert_allclose(forward[2], reverse[2], rtol=0., atol=1e-4)

tests/unit/test_duetsch_protocol.py

@@ -0,0 +1,223 @@
+"""CPU tests of the Dütsch protocol arithmetic (docs/validation/duetsch_re100_kc5_protocol.md §6a)."""
+import importlib.util
+import math
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+_spec = importlib.util.spec_from_file_location(
+    "duetsch_protocol", Path(__file__).resolve().parents[2] / "research/duetsch_re100_kc5/protocol.py")
+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
+
+
+def test_phase_and_cycle_bookkeeping():
+    assert p.phase(0, T) == 270.0 and p.cycle(0, T) == 0          # +A, at rest, phi = -90
+    assert p.phase(T // 4, T) == 0.0 and p.cycle(T // 4, T) == 1  # mean position, u = -U_max
+    assert p.position(0, T, 1.0) == 1.0
+    assert p.cycle(T // 4 - 1, T) == 0 and p.cycle(T // 4 + T - 1, T) == 1
+    assert p.cycle(T // 4 + T, T) == 2
+    assert p.phase(T // 6, T) == pytest.approx(330.0) and p.cycle(T // 6, T) == 0  # startup quarter
+    assert p.phase(7 * T // 6, T) == pytest.approx(330.0) and p.cycle(7 * T // 6, T) == 1
+    for ph, frac in ((180.0, .75), (210.0, 5 / 6)):
+        s = round(frac * T) + T
+        assert p.phase(s, T) == pytest.approx(ph) and p.bin_phase(s, T) == ph
+    assert p.bin_phase(round(.75 * T) + T - 1, T) is None
+    assert sum(p.bin_phase(s, T) == 180.0 for s in range(T, 2 * T)) == T // 360
+
+
+def test_settling_and_common_window():
+    passes = {k: True for k in range(3, 12)}
+    assert p.settled_start(passes) == 7                          # pairs 3->4, 4->5, 5->6
+    late = {**passes, 5: False}
+    assert p.settled_start(late) == 9                            # pairs 6->7, 7->8, 8->9
+    assert p.common_window([7, 9], [passes, passes]) == (9, 10, 11)
+    assert p.common_window([7, 10], [passes, passes]) is None    # K + 2 = 12 beyond the cap
+    assert p.common_window([7, None], [passes, passes]) is None
+    assert p.common_window([7, 7], [passes, {**passes, 8: False}]) is None  # must stay settled
+    # a case settled early must keep passing while waiting for the common K
+    assert p.common_window([7, 9], [{**passes, 8: False}, passes]) is None
+    long = {k: True for k in range(3, 32)}
+    assert p.common_window([10], [long]) is None                 # default cap 11
+    assert p.common_window([29], [long], last_scored=31) == (29, 30, 31)
+    assert p.common_window([30], [long], last_scored=31) is None
+
+
+def test_solid_aware_bilinear():
+    field = np.arange(16.0).reshape(4, 4)
+    fluid = np.ones((4, 4), bool)
+    value, used, total = p.bilinear_fluid(field, fluid, 1.25, 1.5)
+    assert total == pytest.approx(1.0) and value == pytest.approx(0.75 * 0.5 * (5 + 9) + 0.25 * 0.5 * (6 + 10))
+    fluid[1, 1] = False                                           # one solid corner
+    value, used, total = p.bilinear_fluid(field, fluid, 1.25, 1.5)
+    assert total == pytest.approx(1.0 - 0.75 * 0.5)
+    assert all(not (i == 1 and j == 1) for i, j, _ in used)
+    assert sum(w for *_, w in used) == pytest.approx(1.0)
+    fluid[:] = False
+    assert p.bilinear_fluid(field, fluid, 1.25, 1.5) == (None, (), 0.0)
+
+
+def test_interval_bases_and_morison_fit_recover_known_coefficients():
+    b = np.array([p.interval_bases(n, T) for n in range(T // 4, T // 4 + 3 * T)])
+    phi = 2 * np.pi * (np.arange(T // 4, T // 4 + 3 * T) + .5) / T - np.pi / 2
+    assert np.allclose(b[:, 1], np.sin(phi), atol=1e-5)          # midpoint ~ interval mean
+    f = 2.09 * b[:, 0] - (math.pi ** 2 * 1.45 / p.KC) * b[:, 1]
+    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):
+        for c in (1, 5, 11):
+            assert p.position(m + c * T, T, A) == p.position(m, T, A)
+            assert p.phase(m + c * T, T) == p.phase(m, T)
+            v0 = p.position(m + 1, T, A) - p.position(m, T, A)
+            assert p.position(m + 1 + c * T, T, A) - p.position(m + c * T, T, A) == v0
+    assert p.cycle(T // 4 + 5 * T, T) == 6  # unwrapped index stays separate
+    assert p.interval_bases(123 + 7 * T, T) == p.interval_bases(123, T)
+
+
+def test_exact_point_inside_body_is_missing_even_with_fluid_corners():
+    field = np.zeros((6, 6))
+    fluid = np.ones((6, 6), bool)
+    fluid[3, 3] = False                         # only one corner solid
+    centre, radius = (3.0, 3.0), 0.3
+    value, used, total = p.sample(field, fluid, 3.2, 3.1, centre, radius)   # inside the circle
+    assert value is None and used == () and total == 0.0
+    value, used, total = p.sample(field, fluid, 3.35, 3.0, centre, radius)  # just outside
+    assert value is not None and total > 0.0
+
+
+def test_generated_config_passes_the_real_loader(tmp_path):
+    import json
+    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"
+    path.write_text(json.dumps(cfg_dict))
+    cfg = load_lbm_config(str(path))
+    assert (cfg.nx, cfg.ny) == (n, n) == (200, 200)
+    assert cfg.collision == "MRT" and cfg.streaming == "double_buffer" and not cfg.ddf_shifting
+    assert cfg.velocity == 0.0 and cfg.viscosity == pytest.approx(u_max * 20 / 100)
+    assert cfg.y_wall_bc == "bounce_back"
+
+
+def _analyze():
+    return _load("research/duetsch_re100_kc5/analyze.py", "duetsch_analyze")
+
+
+def test_hf_energy_is_parseval_weighted():
+    a = _analyze()
+    n = 360
+    t = np.arange(n)
+    x = 1.0 + np.cos(2 * np.pi * 5 * t / n) + 2.0 * np.cos(2 * np.pi * 30 * t / n) + 0.5 * np.cos(np.pi * t)
+    e = a.one_sided_energy(np.fft.rfft(x), n) / n
+    assert e.sum() == pytest.approx(np.sum(x ** 2))               # Parseval
+    expected = (2.0 ** 2 / 2 + 0.5 ** 2) / np.mean(x ** 2)       # bins >= 21 (30 and Nyquist)
+    assert e[21:].sum() / e.sum() == pytest.approx(expected)
+
+
+def test_horizon_excludes_future_coverage_and_events(monkeypatch):
+    a = _analyze()
+    T, horizon = 360, 4
+    steps = round((0.25 + 6) * T)                                 # cycles 1..6 recorded
+    diag = np.zeros((steps, 15))
+    diag[round((0.25 + horizon) * T):, 13] = 1.0                  # synthetic future no_donor
+    rows = [(s, p.cycle(s, T), 180.0, 0.0, 0.0, 0.1, 0.0, 1.0) for s in range(steps)
+            if p.bin_phase(s, T) is not None and p.cycle(s, T) > p.DISCARD]
+    case = dict(path="c", D=20, T=T, box=10, n=200, u_max=.04, nu=.008, tau=.524, A=1.0,
+                steps=steps, cycles=6, F_ref=1.0, force=np.zeros((steps, 2)), diag=diag,
+                samples=np.asarray(rows, float))
+    monkeypatch.setattr(a, "load", lambda path: case)
+    out = []
+    monkeypatch.setattr("builtins.print", lambda text: out.append(text))
+    a.main(["c", "--horizon", str(horizon)])
+    import json
+    entry = json.loads(out[0])["cases"][0]
+    assert sorted(int(k) for k in entry["coverage"]) == [3, 4]   # nothing from cycles 5, 6
+    assert entry["no_donor"] == 0.0
+
+
+
+def _capture():
+    spec = importlib.util.spec_from_file_location(
+        "duetsch_capture", Path(__file__).resolve().parents[2] / "research/duetsch_re100_kc5/capture.py")
+    mod = importlib.util.module_from_spec(spec)
+    spec.loader.exec_module(mod)
+    return mod
+
+
+def test_capture_crop_downsample_and_trends():
+    c = _capture()
+    n, D = 2400, 40
+    xm = ym = 0.5 * (n - 1)
+    A, R = p.amplitude(D), 20.0
+    box = c.crop_box(xm, ym, D, n, (xm - A - R, xm + A + R, ym - R, ym + R))
+    assert box == (880, 1520, 960, 1440) and box[1] - box[0] + 1 == 641 and box[3] - box[2] + 1 == 481
+    with pytest.raises(ValueError):
+        c.crop_box(xm, ym, D, 600, (0, 1, 0, 1))                    # leaves the domain
+    with pytest.raises(ValueError):
+        c.crop_box(xm, ym, D, n, (xm - 9 * D, xm, ym, ym))           # clips the swept body
+    field = np.arange(64.0).reshape(8, 8)
+    fluid = np.ones((8, 8), bool)
+    fluid[:4, :4] = False
+    down = c.downsample(field, fluid, 4)
+    assert np.isnan(down[0, 0]) and down[1, 1] == field[4:, 4:].mean()
+    rho = np.ones((n, n)); rho[:, 1:81] = 1.01
+    row = c.trend_row(10, 3, 180.0, rho, np.zeros((n, n)), np.zeros((n, n)),
+                      np.ones((n, n), bool), D, n, box)
+    assert row[7] == pytest.approx(1.01) and row[9] == pytest.approx(1.0)
+
+
+def test_dump_steps_are_the_first_step_of_each_bin():
+    T = 5040
+    for k in (11, 21, 30, 31):
+        starts = [s for s in range(T // 4 + (k - 1) * T, T // 4 + k * T)
+                  if p.bin_phase(s, T) is not None and p.phase(s, T) == p.bin_phase(s, T)]
+        assert len(starts) == 3 and all(p.cycle(s, T) == k for s in starts)

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_shared_translation.py

@@ -0,0 +1,122 @@
+"""Host admission and independent lattice-stencil checks; no CUDA context."""
+from types import SimpleNamespace
+
+import numpy as np
+import pytest
+
+from CelerisLab.body.motion import MovingCircle, MovingTables, Translation, validate_shared_translations
+from CelerisLab.lbm.descriptors import FLUID
+from tests.unit.test_moving_bodies import STREAM, body, moving_simulation, routed_simulation
+
+
+def scene(monkeypatch, envelope=None, gap=10.):
+    sim = routed_simulation(monkeypatch)
+    envelope = envelope or Translation((0., 0., -6., 6.), .5, shared=True)
+    for y in (40. - (20. + gap) / 2., 40. + (20. + gap) / 2.):
+        sim.bodies.add(body(60., y, 10., envelope))
+    sim.field = SimpleNamespace(nx=120, ny=80, build_channel_flags=lambda: np.full((80, 120), FLUID))
+    return sim
+
+
+def test_shared_bands_are_admitted_but_other_groups_are_not(monkeypatch):
+    sim = scene(monkeypatch)
+    sim._validate_moving_bodies()
+    a, b = [o.motion for o in sim.bodies.objects]
+    assert a.band()[3] >= b.band()[2]  # whole envelopes genuinely overlap
+    band, geom = MovingTables.layout([(a, 0), (b, 1)])
+    assert band[10] == band[2] * band[3] * 8
+    assert len(geom) == 14
+    b.translation = Translation((0., 0., -6., 6.), .5, shared=True)
+    with pytest.raises(ValueError, match='overlaps body'):
+        sim._validate_moving_bodies()  # equal values do not link independent motion
+    b.translation = a.translation
+    sim.bodies.add(body(60., 40., 2.))
+    with pytest.raises(ValueError, match='overlaps body'):
+        sim._validate_moving_bodies()
+
+
+@pytest.mark.parametrize('gap', (.3, 4.))
+def test_small_surface_gaps_rejected(monkeypatch, gap):
+    with pytest.raises(ValueError, match='clearance'):
+        scene(monkeypatch, gap=gap)._validate_moving_bodies()
+
+
+def test_shared_envelope_is_transverse_only():
+    with pytest.raises(ValueError, match='transverse'):
+        Translation((-.1, .1, -1., 1.), .5, shared=True)
+
+    envelope = Translation((0., 0., -1., 1.), .5, shared=True)
+    with pytest.raises(ValueError, match="zero x velocity"):
+        MovingCircle((0., 0.), 10., envelope).check_target((0., .1), (.1, .1))
+
+def test_group_rejection_is_before_writes_and_can_be_corrected(monkeypatch):
+    sim, writes = moving_simulation(monkeypatch)
+    env = Translation((0., 0., -2., 2.), .5, shared=True)
+    sim.bodies.objects = [body(40., y, 10., env) for y in (20., 50.)]
+    a, b = [o.motion for o in sim.bodies.objects]
+    sim.stepper.moving = object()
+    a.set_target((0., .2), (0., .2))
+    with pytest.raises(ValueError, match='identical'):
+        sim.run(1, stream=STREAM)
+    assert writes == [] and a.applied == b.applied == (0., 0.)
+    b.set_target((0., .2), (0., .1))
+    with pytest.raises(ValueError, match='identical'):
+        sim.run(1, stream=STREAM)
+    assert writes == []
+    b.set_target((0., .2), (0., .2))
+    validate_shared_translations([a, b])
+    b.applied = (0., -.1)
+    with pytest.raises(ValueError, match='identical'):
+        validate_shared_translations([a, b])
+
+
+C = np.array([(x, y) for x in (-1, 0, 1) for y in (-1, 0, 1) if (x, y) != (0, 0)])
+
+
+def footprint_conflicts(centres, radius, dy):
+    """Enumerate geometric write/read supports, not refill/BFL formulas.
+
+    Tests every candidate donor, conservatively including M2's second link.
+    Body-indexed slot numbers are separate from physical (cell,direction) writes.
+    """
+    centres = np.asarray(centres, float)
+    lo = np.floor(centres.min(0) - radius - 4).astype(int)
+    hi = np.ceil(centres.max(0) + radius + 4).astype(int)
+    yy, xx = np.mgrid[lo[1]:hi[1]+1, lo[0]:hi[0]+1]
+    points = np.column_stack((xx.ravel(), yy.ravel()))
+    def inside(p, c):
+        return ((p - c)**2).sum(axis=-1) < radius**2
+    old = np.array([inside(points, c) for c in centres])
+    new_centres = centres + (0., dy)
+    new = np.array([inside(points, c) for c in new_centres])
+    changed = old ^ new
+    conflicts = int(np.count_nonzero(changed.sum(0) > 1))
+    writes = set()
+    for i, (oc, nc) in enumerate(zip(centres, new_centres)):
+        others = [j for j in range(len(centres)) if j != i]
+        fresh = points[old[i] & ~new[i]]
+        for e in C:
+            donor = fresh + e
+            valid = ~inside(donor, oc) & ~inside(donor, nc)
+            for j in others:
+                conflicts += int(np.count_nonzero(valid & (inside(donor, centres[j]) | inside(donor, new_centres[j]))))
+            link = ~new[i] & inside(points + e, nc)
+            fluid = points[link]
+            solid = fluid + e
+            for p in solid:
+                slot = (int(p[0]), int(p[1]), int(e[0]), int(e[1]))
+                conflicts += int(slot in writes)
+                writes.add(slot)
+            for read in (fluid, fluid - e):
+                for j in others:
+                    conflicts += int(np.count_nonzero(inside(read, centres[j]) | inside(read, new_centres[j])))
+    return conflicts
+
+
+def test_counterexample_and_pinball_local_ownership():
+    # Close but disjoint circles are insufficient: parent/Grok independent example.
+    assert footprint_conflicts(((.25, .25), (.25, 8.55)), 4., .5) > 0
+    for phase in (0., .25, .5, .75):
+        for dy in (-.5, -.1, .1, .5):
+            centres = np.array(((0., 0.), (26., 15.), (26., -15.))) + phase
+            assert footprint_conflicts(centres, 10., dy) == 0

tests/unit/test_steady_pinball.py

@@ -0,0 +1,84 @@
+"""Real config loading and pure structural update for the frozen development case."""
+import importlib.util
+import math
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+path = Path(__file__).resolve().parents[2] / 'research/steady_pinball/run.py'
+spec = importlib.util.spec_from_file_location('steady_pinball_run', path)
+runner = importlib.util.module_from_spec(spec)
+spec.loader.exec_module(runner)
+
+
+def test_actual_loader_and_paper_geometry(tmp_path):
+    import json
+    from CelerisLab.config import load_lbm_config
+    p = tmp_path/'config.json'
+    p.write_text(json.dumps(runner.config()))
+    cfg = load_lbm_config(str(p))
+    cfg.validate()
+    assert (cfg.nx, cfg.ny, cfg.viscosity, cfg.velocity) == (1280, 512, .004, .01)
+    assert cfg.collision == 'MRT' and cfg.y_wall_bc == 'bounce_back'
+    assert runner.U * runner.D / runner.NU == 50.
+    assert runner.CENTRES == ((600., 255.5), (626., 270.5), (626., 240.5))
+
+
+def test_structure_rest_equilibrium_and_energy_decay():
+    s = runner.Structure(zeta=.01)
+    assert s.candidate(0.) == (0., 0.)
+    assert s.mass == pytest.approx(10*3*math.pi*100)
+    s.y, s.v = .2, 0.
+    # One damped period: independent continuous solution at total elapsed time.
+    for _ in range(runner.PERIOD):
+        s.y, s.v = s.candidate(0.)
+    t = runner.PERIOD
+    a = s.zeta*s.omega
+    w = s.omega*math.sqrt(1-s.zeta**2)
+    expected = .2*math.exp(-a*t)*(math.cos(w*t)+a/w*math.sin(w*t))
+    assert s.y == pytest.approx(expected, abs=2e-12)
+    assert .5*s.mass*s.v*s.v+.5*s.stiffness*s.y*s.y < .5*s.stiffness*.2**2
+    s.y, s.v = .3, 0.
+    assert s.candidate(s.stiffness*.3) == pytest.approx((.3, 0.), abs=1e-15)
+
+
+def test_report_uses_completed_rows_and_net_force():
+    rows = np.zeros((4, 15))
+    rows[:, 0] = np.arange(1, 5)
+    rows[:, 1] = 2.
+    rows[:, 3:9] = 1.
+    got = runner.summarize(rows)
+    assert got[0]['end_step'] == 4
+    assert got[0]['y_rms_D'] == .1
+    assert got[0]['fy_rms'] == 6.
+
+
+@pytest.mark.parametrize('case', ['Z0', 'M2', 'R2'])
+def test_zero_damping_cases_conserve_unforced_energy_and_match_forced_solution(case):
+    s = runner.Structure(*runner.STRUCTURES[case])
+    assert s.zeta == s.damping == 0.
+    s.y, s.v = .2, 1.e-5
+    y0, v0 = s.y, s.v
+    energy = .5*s.mass*v0**2 + .5*s.stiffness*y0**2
+    for _ in range(12345):
+        s.y, s.v = s.candidate(0.)
+    assert .5*s.mass*s.v**2 + .5*s.stiffness*s.y**2 == pytest.approx(energy, rel=2e-11)
+    force, elapsed = 1.e-4, 3210
+    s.y, s.v = y0, v0
+    for _ in range(elapsed):
+        s.y, s.v = s.candidate(force)
+    phase = s.omega*elapsed
+    equilibrium = force/s.stiffness
+    expected_y = equilibrium+(y0-equilibrium)*math.cos(phase)+v0*math.sin(phase)/s.omega
+    expected_v = -(y0-equilibrium)*s.omega*math.sin(phase)+v0*math.cos(phase)
+    assert (s.y, s.v) == pytest.approx((expected_y, expected_v), abs=2e-11)
+
+
+def test_predeclared_comparisons_change_structure_only():
+    z, m, r = [runner.Structure(*runner.STRUCTURES[k]) for k in ['Z0', 'M2', 'R2']]
+    assert m.mass/z.mass == pytest.approx(.2)
+    assert m.stiffness/z.stiffness == pytest.approx(.2)
+    assert r.mass == m.mass
+    assert r.stiffness/m.stiffness == pytest.approx(.36)
+    assert [runner.STEPS*s.frequency for s in [z, m, r]] == pytest.approx([60., 60., 36.])

tests/unit/test_viv_logic.py

@@ -0,0 +1,105 @@
+"""CPU tests of the predeclared VIV logic (research/viv_re150/viv.py; stdlib only)."""
+import importlib.util
+import math
+from pathlib import Path
+
+_spec = importlib.util.spec_from_file_location(
+    "viv_logic", Path(__file__).resolve().parents[2] / "research/viv_re150/viv.py")
+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
+
+
+def _feed(cyc, amps, period=100):
+    """Feed sine cycles of the given amplitudes (in D); return the first non-None result."""
+    n = 0
+    for a in amps:
+        for i in range(period):
+            out = cyc.update(a * D * math.sin(2 * math.pi * (i + 1) / period))
+            n += 1
+            if out:
+                return out, n
+    return None, n
+
+
+def test_structure_exact_update_conserves_free_oscillation_and_equilibrium():
+    s = viv.Structure(D, U)
+    assert math.isclose(1 / s.f_n, 5000.0) and math.isclose(s.m, 2 * D * D)
+    s.y, s.v = 0.3 * D, 0.0
+    e0 = 0.5 * s.k * s.y ** 2
+    for _ in range(12345):
+        s.y, s.v = s.candidate(0.0)
+    assert math.isclose(0.5 * s.k * s.y ** 2 + 0.5 * s.m * s.v ** 2, e0, rel_tol=1e-9)
+    s.y, s.v = 2.0 / s.k, 0.0
+    assert all(math.isclose(a, b, abs_tol=1e-12) for a, b in zip(s.candidate(2.0), (2.0 / s.k, 0.0)))
+    assert math.isclose(s.seed_velocity(D) / s.w, viv.ELIGIBLE * D)
+
+
+def test_constant_amplitude_settles_then_scores_twenty_cycles():
+    cyc = viv.Cycles(D)
+    out, _ = _feed(cyc, [0.5] * 60)
+    assert out == "scored" and len(cyc.window) == 20 and cyc.settled_at == 1
+    r = viv.score(cyc.window)
+    assert math.isclose(r["A_max"], 0.5, rel_tol=1e-3) and abs(r["offset"]) < 1e-3
+    assert r["T_resp"] == 100 and r["within_band"]          # |.500 - .520| <= .052
+
+
+def test_small_cycles_are_not_counted_and_y0_is_not_a_crossing():
+    cyc = viv.Cycles(D)
+    assert cyc.update(0.0) is None and cyc.log == []
+    out, _ = _feed(cyc, [0.005] * 50)
+    assert out is None and cyc.counted == 0 and cyc.log and not any(c["eligible"] for c in cyc.log)
+
+
+def test_scoring_failure_restarts_from_the_failing_block():
+    cyc = viv.Cycles(D)
+    amps = [0.5] * 20 + [0.55] * 10 + [0.55] * 40      # settle at block 1, block 2 jumps 10 %
+    out, _ = _feed(cyc, amps)
+    assert out == "scored"
+    assert all(math.isclose(c["amp"], 0.55, rel_tol=1e-3) for c in cyc.window)
+
+
+def test_seed_reset_discards_incomplete_cycle_and_proof():
+    cyc = viv.Cycles(D)
+    _feed(cyc, [0.5] * 15)
+    cyc.update(0.2 * D)                                  # inside an incomplete cycle
+    counted = cyc.counted
+    cyc.discard_current()
+    cyc.reset_proof()
+    assert cyc.amps == [] and cyc.blocks == [] and cyc.counted == counted
+
+
+def test_viv_config_passes_the_real_loader(tmp_path):
+    import json
+    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"
+    path.write_text(json.dumps(cfg_dict))
+    cfg = load_lbm_config(str(path))
+    assert (cfg.nx, cfg.ny) == (nx, ny) == (3000, 2400)
+    assert cfg.collision == "MRT" and cfg.y_wall_bc == "free_slip" and cfg.velocity == .04
+    assert math.isclose(cfg.viscosity, .04 * 40 / 150)
+
+
+def test_finite_rejects_any_nonfinite_recorded_output():
+    row = (3, .1, .2, .3, .4, .5, .6, .7, 0.0)
+    assert viv.finite(row)
+    for i in range(len(row)):
+        for bad in (float("nan"), float("inf"), -float("inf")):
+            assert not viv.finite(row[:i] + (bad,) + row[i + 1:])