research: VIV driver — fail closed on any non-finite recorded output

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

docs/moving_bodies.md

@@ -338,6 +338,63 @@ recorded as not runnable rather than replaced by another test.
   run with a moving body raises before upload/launch; stationary unaffected). No kernels, no GPU runs, no
   y-only VIV oscillator (S4), no force RHS.
 
+## 10a. S2a — fixed-centre moving links (implemented on `claude/moving-bodies-s2a`)
+
+- `MovingCircleLinks` (`kernels/step/moving_circle.cu`, dense double buffer)
+  runs before the static curved kernel when moving tables are attached: one
+  thread per (band cell, direction) slot, FP64 `{cx, cy, r}` at the applied
+  offset, q/hit/clamp/fallback in the order of `build_cut_links` with explicit
+  `_rn` operations (host/device bit parity is a test target, not a claim),
+  then the unchanged `apply_bouzidi_link` with spin wall velocity only.
+- `obs` keeps its contract (CME and torque from `apply_bouzidi_link`). A
+  separate per-body diagnostics buffer receives `cme`, the Wen term
+  `u_w * dm`, `link_dm` and the link count from the same link outputs; it is
+  zeroed with obs and only exists when moving circles do. Optional
+  `slot_debug` (q, fallback, rx, ry per slot) is null in production.
+- Displacement is disabled: `set_body(offset=..., vx=..., vy=...)` on a
+  translating circle raises before mutation, before and after initialize, so
+  no pending target can exist and `run`/`step`/direct `LBMStepper.step` all
+  hold at the applied geometry. `initialize()` never reads or consumes
+  `pending`. `stepper.blocked` is cleared only when the moving tables are
+  attached (initialize, stepper recreation).
+- S3 rule (not implemented): once targets are re-enabled, any pending target
+  must raise in `run`/`step` and in `LBMStepper.step` alike before any launch,
+  swap or step-count change; zero-step calls never consume it.
+- Not validation: a pass of the zero-motion gates says nothing about moving
+  geometry. GPU parity (G0a/G0b) is unverified until run on a free device.
+
+## 10b. S3a — translation with refill (implemented on `claude/moving-bodies-s3a`)
+
+- Targets are enabled again (4a). `set_body(offset, vx, vy)` validates against
+  `applied`; a pending target blocks direct `LBMStepper.step(n > 0)` before any
+  launch; `run(n > 1)` with a target raises before mutation; `run(0)` never
+  consumes it; `run(1)` advances it.
+- Per advancing step, on one stream: upload the immutable FP64 records
+  `{c_old, c_new, r, u}` (only when changed) → M1 `MovingCircleCells` → M2
+  `MovingCircleLinks` (`u_w = u + omega x r` at the hit point) → static curved →
+  `OneStep` → swap → force region → sensors → synchronize → read `no_donor` →
+  commit (`applied := pending.offset`) → upload the hold record
+  (`c_old = c_new = applied`, `u = 0`) → synchronize. Hold steps launch no M1
+  and are the S2a path. Commit point: the advancing step-group sync plus a
+  zero `no_donor` count; a failure before it commits nothing. The hold-record
+  upload stages the next step after the commit: if it (or its sync) fails,
+  `applied` already matches the moved device geometry, but the Simulation is
+  still marked failed and refuses reuse. Any failure from the first upload on
+  marks the circles and the Simulation failed.
+- M1 classifies with the old/new predicates only (no flag reads). Covered
+  cells add raw `sum c G` / `sum G` of `G^{n-1}` to `cover_*` and become
+  `SOLID|OBSTACLE|BC_CURVED`. Uncovered cells: donors are fluid in both
+  geometries; `rho_bar` over all of them; NEQ from the donor with the largest
+  unnormalized `c·n > 0` (first direction wins ties, `n` from `c_new`);
+  `u_w` (translation + spin) at the cell position; all nine slots written;
+  raw `sum c G`, `sum G` of the refill go to `uncover_*`; no donor →
+  `f_eq(RHO, u_w)` plus `no_donor` (hard failure before commit); no NEQ donor
+  → `eq_only`.
+- `get_flags()` reads the device flags back (synchronized) when moving tables
+  exist; otherwise it returns the host array as before.
+- Not selected: the structural force. CME+node vs CME+Wen stay OPEN until X1
+  and G4 have run.
+
 ## 11. Remaining decisions
 
 Structural RHS (section 6); motion-record layout; refill details; coupling

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

@@ -0,0 +1,146 @@
+# Single-cylinder free transverse VIV — PREDECLARATION (agreed at 1abfe954)
+
+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.
+- Box [0, 75 D] × [0, 60 D] (3000 × 2400) with the cylinder centre at
+  (30 D, 30 D): 30 D from the inlet, 45 D from the outlet, centred in y
+  (Kolahdouz §4.1.1). Qin's own box ([0, 40 D] × [0, 20 D], centre (10 D, 10 D))
+  is the target paper's, not ours.
+- 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,129 @@
+"""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
+
+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
+
+
+def interval_bases(n, T, points=16):
+    """Single-step interval averages over [n, n + 1] of cos phi |cos phi| and sin phi
+    (phi in radians = 2 pi s / T - pi / 2), Gauss-Legendre quadrature."""
+    x, w = np.polynomial.legendre.leggauss(points)
+    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/viv_re150/run.py

@@ -0,0 +1,134 @@
+"""One-shot predeclared VIV run (docs/validation/leontini_viv_predeclare.md).
+
+Re 150, U* = 5, m/(rho D^2) = 2, zeta = 0, D = 40 on a 75 D x 60 D native
+channel (uniform inlet, NEQ outlet, free-slip y walls), cylinder at
+(30 D, 30 D). Structural RHS = CME + Wen (diagnostics: raw CME, node terms).
+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 Translation
+    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[1] + d[3])             # CME + Wen (y)
+            row = (n, y1, v1, f_rhs, float(d[1]), float(d[3]), float(f_obs),
+                   float(d[7] - d[10]), float(d[13]))
+            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/leontini_viv_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

@@ -111,14 +111,21 @@ class MovingCircle:
             step()
             synchronize()
         except BaseException as exc:
-            self.failed = True
-            if self._on_failure is not None:
-                self._on_failure(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:
             self.applied = self.pending[0]
             self.pending = None
 
+    def fail(self, exc: BaseException) -> None:
+        self.failed = True
+        if self._on_failure is not None:
+            self._on_failure(exc)
+
     def band(self) -> Tuple[int, int, int, int]:
         """Inclusive cell AABB (x0, x1, y0, y1) swept over the whole offset box."""
         x0, x1, y0, y1 = self.translation.offset_box
@@ -166,3 +173,87 @@ def enumerate_donor_coverage(radius: float, positions: int = 8, directions: int
                 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``.
+
+    Band table (x0, y0, width, height, first_slot, body_id), immutable per run;
+    FP64 motion records (cx_old, cy_old, cx_new, cy_new, r, ux, uy), immutable
+    during a step; separate per-body diagnostics (see DIAG), never part of obs.
+    """
+    BAND_INTS, GEOM_DOUBLES = 6, 7
+    DIAG = ("cme_x", "cme_y", "wen_x", "wen_y", "link_dm", "n_links",
+            "cover_jx", "cover_jy", "cover_mass", "uncover_jx", "uncover_jy", "uncover_mass",
+            "eq_only", "no_donor", "crossings")
+    DIAG_FLOATS = len(DIAG)
+
+    def __init__(self, band, geom):
+        import numpy as np
+        import pycuda.driver as cuda
+        self.n_bodies = len(geom) // self.GEOM_DOUBLES
+        self.n_slots = int(band[-2] + band[-4] * band[-3] * 8)
+        self.n_cells = self.n_slots // 8
+        self.geom = list(geom)  # host copy of the uploaded records
+        self.band_gpu = cuda.to_device(np.asarray(band, np.int32))
+        self.geom_gpu = cuda.to_device(np.asarray(geom, np.float64))
+        # Owned page-locked staging for async record uploads (allocated once).
+        self._geom_pinned = cuda.pagelocked_empty(len(geom), np.float64)
+        self.diag_gpu = cuda.to_device(np.zeros(self.n_bodies * self.DIAG_FLOATS, np.float32))
+        self.slot_debug_gpu = None  # optional (q, fallback, rx, ry) per slot for parity tests
+
+    @classmethod
+    def upload(cls, circles):
+        """``circles``: [(MovingCircle, body_id)] in body order."""
+        return cls(*cls.layout(circles))
+
+    @staticmethod
+    def layout(circles):
+        """Host band and geometry tables (no device work)."""
+        band, geom, first = [], [], 0
+        for motion, body_id in circles:
+            x0, x1, y0, y1 = motion.band()
+            width, height = x1 - x0 + 1, y1 - y0 + 1
+            band += [x0, y0, width, height, first, int(body_id)]
+            first += width * height * 8
+            geom += MovingTables.record(motion)
+        return band, geom
+
+    @staticmethod
+    def record(motion):
+        """FP64 motion record of one circle for its next step (hold if nothing pending)."""
+        c_old, c_new, u = motion.record()
+        return [*c_old, *c_new, motion.radius, *u]
+
+    def upload_records(self, circles, stream) -> bool:
+        """Enqueue the records when they differ from the device copy; True if advancing.
+
+        The pinned staging buffer is reused: callers synchronize ``stream`` after every
+        upload before the next one (``Simulation._run_moving``), so no copy is in flight
+        when it is rewritten. ``geom`` is the enqueued state; after a failed copy the
+        owner is failed and the cache is never trusted again.
+        """
+        import pycuda.driver as cuda
+        geom = [v for motion, _ in circles for v in self.record(motion)]
+        if geom != self.geom:
+            self._geom_pinned[:] = geom
+            cuda.memcpy_htod_async(self.geom_gpu, self._geom_pinned, stream)
+            self.geom = geom
+        g = self.GEOM_DOUBLES
+        return any(geom[i:i + 2] != geom[i + 2:i + 4] for i in range(0, len(geom), g))
+
+    def enable_slot_debug(self):
+        import numpy as np
+        import pycuda.driver as cuda
+        self.slot_debug_gpu = cuda.to_device(np.zeros(self.n_slots * 4, np.float32))
+
+    def zero_diagnostics(self, stream):
+        import pycuda.driver as cuda
+        cuda.memset_d32_async(self.diag_gpu, 0, self.n_bodies * self.DIAG_FLOATS, stream)
+
+    def download_diagnostics(self):
+        import numpy as np
+        import pycuda.driver as cuda
+        out = np.empty((self.n_bodies, self.DIAG_FLOATS), np.float32)
+        cuda.memcpy_dtoh(out, self.diag_gpu)
+        return out

src/CelerisLab/lbm/kernels/step/moving_circle.cu

@@ -0,0 +1,183 @@
+// Translating circles (docs/moving_bodies.md): M1 MovingCircleCells reclassifies
+// band cells and refills uncovered ones; M2 MovingCircleLinks applies the cut links.
+// Geometry comes from an immutable FP64 record per body
+// {cx_old, cy_old, cx_new, cy_new, r, ux, uy}; the M2 operation order mirrors
+// CircleGeometry.build_cut_links / find_circle_intersection with explicit _rn
+// operations (bit parity with the host builder is a test target, not a claim).
+// Moving circles are never in the packed static links, so this kernel is the
+// only writer of their solid slots.
+#ifndef CELERIS_MOVING_CIRCLE_CU
+#define CELERIS_MOVING_CIRCLE_CU
+
+#if DIM == 2
+constexpr int MOVING_BAND_INTS = 6;   // x0, y0, width, height, first_slot, body_id
+constexpr int MOVING_GEOM_DOUBLES = 7;
+// cme_x, cme_y, wen_x, wen_y, link_dm, n_links, cover_jx, cover_jy, cover_mass,
+// uncover_jx, uncover_jy, uncover_mass, eq_only, no_donor, crossings
+constexpr int MOVING_DIAG_FLOATS = 15;
+
+__device__ __forceinline__ double mc_sq(double v) { return __dmul_rn(v, v); }
+
+__device__ __forceinline__ bool mc_inside(int x, int y, double cx, double cy, double r)
+{
+    return __dadd_rn(mc_sq(__dsub_rn((double)x, cx)), mc_sq(__dsub_rn((double)y, cy))) < mc_sq(r);
+}
+
+// Smallest root t in [0, 1] of |p + t c - centre| = r (find_circle_intersection order).
+__device__ __forceinline__ bool mc_hit(int x, int y, int ex, int ey,
+                                       double cx, double cy, double r, double& t)
+{
+    const double dx = (double)ex, dy = (double)ey;
+    const double px = __dsub_rn((double)x, cx), py = __dsub_rn((double)y, cy);
+    const double a = __dadd_rn(mc_sq(dx), mc_sq(dy));
+    const double b = __dmul_rn(2.0, __dadd_rn(__dmul_rn(dx, px), __dmul_rn(dy, py)));
+    const double c = __dsub_rn(__dadd_rn(mc_sq(px), mc_sq(py)), mc_sq(r));
+    const double det = __dsub_rn(mc_sq(b), __dmul_rn(__dmul_rn(4.0, a), c));
+    if (det < 0.0) return false;
+    const double s = __dsqrt_rn(det), a2 = __dmul_rn(2.0, a);
+    const double t1 = __ddiv_rn(__dsub_rn(-b, s), a2);
+    const double t2 = __ddiv_rn(__dadd_rn(-b, s), a2);
+    if (0.0 <= t1 && t1 <= 1.0) { t = t1; return true; }
+    if (0.0 <= t2 && t2 <= 1.0) { t = t2; return true; }
+    return false;
+}
+
+__global__ void MovingCircleLinks(
+    fpxx* fi, const uint16_t* flag,
+    const int* band, const double* geom,
+    unsigned int n_bodies, unsigned int n_slots,
+    const float* action, float* obs, float* diag, float* slot_debug)
+{
+    const unsigned int tid = threadIdx.x + blockIdx.x * blockDim.x;
+    if (tid >= n_slots) return;
+    unsigned int m = 0;
+    while (m + 1 < n_bodies && tid >= (unsigned int)band[(m + 1) * MOVING_BAND_INTS + 4]) ++m;
+    const int* b = band + m * MOVING_BAND_INTS;
+    const unsigned int local = tid - (unsigned int)b[4];
+    const unsigned int dir = local % 8u + 1u;
+    const int x = b[0] + (int)((local / 8u) % (unsigned int)b[2]);
+    const int y = b[1] + (int)((local / 8u) / (unsigned int)b[2]);
+    const double* g = geom + MOVING_GEOM_DOUBLES * m;
+    const double cx = g[2], cy = g[3], r = g[4];
+    const int ex = d_cx[dir], ey = d_cy[dir];
+    float* debug = slot_debug == nullptr ? nullptr : slot_debug + 4 * tid;
+    if (debug != nullptr) debug[0] = -1.0f;
+
+    double t;
+    if (mc_inside(x, y, cx, cy, r) || !mc_inside(x + ex, y + ey, cx, cy, r)
+            || !mc_hit(x, y, ex, ey, cx, cy, r, t))
+        return;
+    const double hx = __dadd_rn((double)x, __dmul_rn(t, (double)ex));
+    const double hy = __dadd_rn((double)y, __dmul_rn(t, (double)ey));
+    const double hdx = __dsub_rn(hx, (double)x), hdy = __dsub_rn(hy, (double)y);
+    const double norm = __dsqrt_rn((double)(ex * ex + ey * ey));
+    double qd = __ddiv_rn(__dsqrt_rn(__dadd_rn(mc_sq(hdx), mc_sq(hdy))), norm);
+    qd = fmax(0.01, fmin(qd, 1.0));
+    const float q = (float)qd;
+    const float rx = (float)__dsub_rn(hx, cx), ry = (float)__dsub_rn(hy, cy);
+
+    // Same two-stage order as build_cut_links + pack_cut_links_to_soa: the fallback
+    // class is chosen from the clamped FP64 q (q < 1/2 needs a fluid donor at x - c);
+    // only then are q, rx, ry rounded to FP32 for the BFL branch and spin velocity.
+    unsigned char fallback = CURVED_FALLBACK_BOUZIDI;
+    if (qd < 0.5) {
+        const int xd = x - ex, yd = y - ey;
+        if (xd < 0 || xd >= (int)NX || yd < 0 || yd >= (int)NY || mc_inside(xd, yd, cx, cy, r)
+                || (flag[linear_index((unsigned int)xd, (unsigned int)yd)] & FLAG_FLUID) == 0)
+            fallback = CURVED_FALLBACK_HALFWAY;
+    }
+    if (debug != nullptr) { debug[0] = q; debug[1] = (float)fallback; debug[2] = rx; debug[3] = ry; }
+
+    const int bid = b[5];
+    const float omega = action_omega(action, bid);
+    const float Uw = (float)g[5] - omega * ry, Vw = (float)g[6] + omega * rx;  // u + omega x r
+    float exchange[3], dm;
+    apply_bouzidi_link(dir, q, linear_index((unsigned int)x, (unsigned int)y), fi,
+                       Uw, Vw, rx, ry, fallback, obs, bid, exchange, &dm);
+    if (diag != nullptr) {
+        float* d = diag + m * MOVING_DIAG_FLOATS;
+        atomicAdd(&d[0], exchange[0]);
+        atomicAdd(&d[1], exchange[1]);
+        atomicAdd(&d[2], Uw * dm);   // Wen 2014 correction -u_w (f_t - f_r)
+        atomicAdd(&d[3], Vw * dm);
+        atomicAdd(&d[4], dm);
+        atomicAdd(&d[5], 1.0f);
+    }
+}
+// M1: one thread per band cell, launched only on steps with c_old != c_new.
+// Membership comes from the immutable old/new predicates, never from flags, so
+// classification, refill and the flag write of one cell cannot race with another.
+// ddf holds post-collision G^{n-1}; refilled cells get all nine slots at that level.
+__global__ void MovingCircleCells(
+    fpxx* fi, uint16_t* flag, const int* band, const double* geom,
+    unsigned int n_bodies, unsigned int n_cells, const float* action, float* diag)
+{
+    const unsigned int tid = threadIdx.x + blockIdx.x * blockDim.x;
+    if (tid >= n_cells) return;
+    unsigned int m = 0;
+    while (m + 1 < n_bodies && tid >= (unsigned int)band[(m + 1) * MOVING_BAND_INTS + 4] / 8u) ++m;
+    const int* b = band + m * MOVING_BAND_INTS;
+    const unsigned int local = tid - (unsigned int)b[4] / 8u;
+    const int x = b[0] + (int)(local % (unsigned int)b[2]);
+    const int y = b[1] + (int)(local / (unsigned int)b[2]);
+    const double* g = geom + MOVING_GEOM_DOUBLES * m;
+    const bool in_old = mc_inside(x, y, g[0], g[1], g[4]);
+    const bool in_new = mc_inside(x, y, g[2], g[3], g[4]);
+    if (in_old == in_new) return;
+    const unsigned long k = linear_index((unsigned int)x, (unsigned int)y);
+    float* d = diag + m * MOVING_DIAG_FLOATS;
+    float rho = 0.0f, jx = 0.0f, jy = 0.0f;
+    if (in_new) {  // covered: its fluid content leaves the active set
+        for (unsigned int i = 0; i < NQ; ++i) {
+            const float f = load_ddf(fi, index_f(k, i));
+            rho += f; jx += (float)d_cx[i] * f; jy += (float)d_cy[i] * f;
+        }
+        flag[k] = FLAG_SOLID | FLAG_OBSTACLE | FLAG_BC_CURVED;
+        atomicAdd(&d[6], jx); atomicAdd(&d[7], jy); atomicAdd(&d[8], rho);
+        atomicAdd(&d[14], 1.0f);
+        return;
+    }
+    // Uncovered: donors are fluid in BOTH geometries; NEQ from the neighbour with the
+    // largest unnormalized c.n > 0 (first direction wins ties), n from c_new.
+    const float nx = (float)__dsub_rn((double)x, g[2]), ny = (float)__dsub_rn((double)y, g[3]);
+    float rho_sum = 0.0f, best_dot = 0.0f;
+    int n_valid = 0, best = -1;
+    for (unsigned int i = 1; i < NQ; ++i) {
+        const int xn = x + d_cx[i], yn = y + d_cy[i];
+        if (mc_inside(xn, yn, g[0], g[1], g[4]) || mc_inside(xn, yn, g[2], g[3], g[4])) continue;
+        const unsigned long kn = linear_index((unsigned int)xn, (unsigned int)yn);
+        for (unsigned int q = 0; q < NQ; ++q) rho_sum += load_ddf(fi, index_f(kn, q));
+        ++n_valid;
+        const float dot = (float)d_cx[i] * nx + (float)d_cy[i] * ny;
+        if (dot > best_dot) { best_dot = dot; best = (int)i; }
+    }
+    const float omega = action_omega(action, b[5]);
+    const float uwx = (float)g[5] - omega * ny, uwy = (float)g[6] + omega * nx;  // at the cell
+    float feq[NQ], neq[NQ];
+    compute_feq(n_valid > 0 ? rho_sum / (float)n_valid : (float)RHO, uwx, uwy, feq);
+    for (unsigned int i = 0; i < NQ; ++i) neq[i] = 0.0f;
+    if (best >= 0) {
+        const unsigned long kd = linear_index((unsigned int)(x + d_cx[best]), (unsigned int)(y + d_cy[best]));
+        float gd[NQ], rd = 0.0f, mx = 0.0f, my = 0.0f, fd[NQ];
+        for (unsigned int i = 0; i < NQ; ++i) {
+            gd[i] = load_ddf(fi, index_f(kd, i));
+            rd += gd[i]; mx += (float)d_cx[i] * gd[i]; my += (float)d_cy[i] * gd[i];
+        }
+        compute_feq(rd, mx / rd, my / rd, fd);
+        for (unsigned int i = 0; i < NQ; ++i) neq[i] = gd[i] - fd[i];
+    } else {
+        atomicAdd(&d[12], 1.0f);
+    }
+    if (n_valid == 0) atomicAdd(&d[13], 1.0f);  // hard failure, checked before commit
+    for (unsigned int i = 0; i < NQ; ++i) {
+        const float f = feq[i] + neq[i];
+        store_ddf(fi, index_f(k, i), f);
+        rho += f; jx += (float)d_cx[i] * f; jy += (float)d_cy[i] * f;
+    }
+    flag[k] = FLAG_FLUID;
+    atomicAdd(&d[9], jx); atomicAdd(&d[10], jy); atomicAdd(&d[11], rho);
+    atomicAdd(&d[14], 1.0f);
+}
+#endif  // DIM == 2
+
+#endif  // CELERIS_MOVING_CIRCLE_CU

src/CelerisLab/lbm/kernels/uniform.cu

@@ -71,6 +71,7 @@ extern "C"
 #include "step/one_step_double.cu"
 #include "step/one_step_esopull.cu"
 #include "step/aux_kernels.cu"
+#include "step/moving_circle.cu"
 #include "step/esopull_macro.cu"
 #include "step/esopull_sync.cu"
 

src/CelerisLab/lbm/stepper.py

@@ -23,6 +23,7 @@ class LBMStepper:
         # Set by the owner when a step would be wrong (translating circles
         # excluded from the packed links before the moving kernels exist).
         self.blocked: str | None = None
+        self.moving = None  # MovingTables attached by the owner (translating circles)
 
         if self._esopull:
             self.step_fn = module.get_function("EsoPullStep")
@@ -38,6 +39,8 @@ 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")
 
         tpb = cfg.threads_per_block
         self.block = (tpb, 1, 1)
@@ -75,6 +78,11 @@ class LBMStepper:
         """Advance *n* time steps."""
         if n > 0 and self.blocked is not None:
             raise RuntimeError(self.blocked)
+        self._step(n, action_gpu=action_gpu, obs_gpu=obs_gpu, stream=stream)
+
+    def _step(self, n, *, action_gpu, obs_gpu, stream=None, advance=False):
+        """Unguarded steps; ``advance`` runs one step that moves translating circles."""
+        assert not advance or n == 1
         f = self.field
         launch_kw = {}
         if stream is not None:
@@ -86,6 +94,8 @@ class LBMStepper:
                     f.ddf_gpu, self.boundary_values_gpu, np.uint64(self._step_count),
                     block=self.block, grid=self.boundary_grid, **launch_kw,
                 )
+            if self.moving is not None:
+                self._launch_moving(action_gpu, obs_gpu, advance, **launch_kw)
             self._launch_curved(action_gpu, obs_gpu, **launch_kw)
 
             if self._esopull:
@@ -111,6 +121,21 @@ class LBMStepper:
             # Any completed GPU step invalidates host DDF mirror.
             f.invalidate_host_ddf_cache()
 
+    def _launch_moving(self, action_gpu, obs_gpu, advance, **launch_kw):
+        m, tpb = self.moving, self.cfg.threads_per_block
+        if advance:  # M1 completes before M2 on the same stream
+            self.moving_cells_fn(
+                self.field.ddf_gpu, self.field.flag_gpu, m.band_gpu, m.geom_gpu,
+                np.uint32(m.n_bodies), np.uint32(m.n_cells), action_gpu, m.diag_gpu,
+                block=(tpb, 1, 1), grid=((m.n_cells + tpb - 1) // tpb, 1, 1), **launch_kw,
+            )
+        self.moving_fn(
+            self.field.ddf_gpu, self.field.flag_gpu, m.band_gpu, m.geom_gpu,
+            np.uint32(m.n_bodies), np.uint32(m.n_slots), action_gpu, obs_gpu,
+            m.diag_gpu, m.slot_debug_gpu if m.slot_debug_gpu is not None else np.uintp(0),
+            block=(tpb, 1, 1), grid=((m.n_slots + tpb - 1) // tpb, 1, 1), **launch_kw,
+        )
+
     def _launch_curved(self, action_gpu, obs_gpu, **launch_kw):
         f = self.field
         if f.n_curved == 0:

src/CelerisLab/simulation.py

@@ -28,8 +28,10 @@ from .cuda.context import CudaContext, context_bound
 from .cuda import compiler
 from .lbm.descriptors import FLUID
 
-# S1: translating circles are declared and validated, but no moving kernel exists yet.
-MOVING_BLOCKED = ("translating circles need the moving kernels (not implemented); "
+# Translating circles step only through their moving tables (fixed centre in S2a).
+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); "
                   "positive-step runs are rejected")
 from .lbm.field import LBMField
 from .lbm.stepper import LBMStepper
@@ -211,6 +213,37 @@ class Simulation:
     def _fail_moving(self, exc: BaseException) -> None:
         # Device geometry is unknown after a failed moving step: same lifecycle as Jets.
         self._failed = f"moving circle step failed: {exc!r}"
+        self.stepper.blocked = self._failed  # direct stepper calls are refused as well
+
+    def _attach_moving(self) -> None:
+        # Direct stepper calls stay blocked until the moving tables exist.
+        tables = getattr(self, "_moving_tables", None)
+        self.stepper.moving = tables
+        moving = self._moving_bodies()
+        if moving and tables is None:
+            self.stepper.blocked = MOVING_BLOCKED
+        elif any(obj.motion.pending is not None for obj in moving):
+            self.stepper.blocked = PENDING_BLOCKED  # direct steps never consume a target
+        else:
+            self.stepper.blocked = None
+
+    def _build_moving_tables(self):
+        moving = self._moving_bodies()
+        if not moving:
+            return None
+        from .body.motion import MovingTables
+        return MovingTables.upload([(obj.motion, obj.obj_id) for obj in moving])
+
+    @context_bound
+    def moving_diagnostics(self) -> np.ndarray:
+        """Per translating circle (in body order), columns ``MovingTables.DIAG``,
+        accumulated since obs were last zeroed (diagnostics only, not a force RHS)."""
+        self._require_usable()
+        tables = self.stepper.moving
+        if tables is None:
+            return np.zeros((0, 6), np.float32)
+        self.stream.synchronize()
+        return tables.download_diagnostics()
 
     def _moving_bodies(self):
         return [obj for obj in self.bodies.objects if obj.motion is not None]
@@ -331,7 +364,8 @@ class Simulation:
             vy: translation velocity y (future).  ``None`` = no change.
             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).
+                against the applied offset (docs/moving_bodies.md, 4a). Only
+                ``run(1)`` advances it; direct stepper calls are blocked meanwhile.
         """
         obj = self.bodies.get(id)
         motion = obj.motion
@@ -360,6 +394,7 @@ class Simulation:
             changed = True
         if motion is not None and target is not None:
             motion.pending = target
+            self._attach_moving()
         if changed:
             self.bodies._refresh_action_from_objects()
 
@@ -516,7 +551,7 @@ class Simulation:
             self.field, self._module, self.lbm_cfg,
         )
         self.stepper._step_count = _prev_step_count
-        self.stepper.blocked = MOVING_BLOCKED if self._moving_bodies() else None
+        self._attach_moving()
 
         self._assert_object_count_contract(expected_count=expected_n_objects)
 
@@ -543,7 +578,7 @@ class Simulation:
         # change. Compact-list construction publishes no outputs on failure.
         next_flags = self.bodies.build_flags(self.field.build_channel_flags())
         self._validate_moving_bodies()
-        self.stepper.blocked = MOVING_BLOCKED if self._moving_bodies() else None
+        self._attach_moving()
         if self.bodies.jets and self.bodies.count == 0:
             raise ValueError("Jets need bodies with cut links")
         self.bodies.build_compact_lists_for(self.bodies.objects, next_flags)
@@ -560,6 +595,8 @@ class Simulation:
             self.bodies.sync_to_gpu(self.field, rebuild_flags=False)
         else:
             self.bodies.ensure_packed_buffers()
+        self._moving_tables = self._build_moving_tables()
+        self._attach_moving()
         self._assert_runtime_contracts()
         self._initialized = True
 
@@ -676,8 +713,13 @@ class Simulation:
         checkpoint_interval = operator.index(checkpoint_interval)
         if steps < 0 or checkpoint_interval < 0:
             raise ValueError("steps and checkpoint_interval must be nonnegative")
-        if steps > 0 and self._moving_bodies():
+        moving = self._moving_bodies()
+        if steps > 0 and moving and self.stepper.moving is None:
             raise RuntimeError(MOVING_BLOCKED)
+        if moving and checkpoint_interval > 0:
+            raise ValueError("checkpoints are not supported with translating circles")
+        if steps > 1 and any(obj.motion.pending is not None for obj in moving):
+            raise ValueError("a pending translation target needs run(1)")
         self._assert_runtime_contracts()
         # Discard any uncommitted body edits before stepping.
         if self.bodies.has_pending_edit():
@@ -694,6 +736,8 @@ class Simulation:
         # Zero obs segments before step group
         if zero_obs:
             self.bodies.zero_obs_async(stream)
+            if self.stepper.moving is not None:
+                self.stepper.moving.zero_diagnostics(stream)
 
         if checkpoint_interval > 0:
             done = 0
@@ -710,6 +754,9 @@ class Simulation:
                 if done < steps or done % checkpoint_interval == 0:
                     stream.synchronize()
                     self.save_checkpoint()
+        elif self.stepper.moving is not None and steps > 0:
+            self._run_moving(steps, stream, sync_obs)
+            return
         else:
             self.stepper.step(
                 steps,
@@ -728,6 +775,34 @@ class Simulation:
 
         stream.synchronize()
 
+    def _run_moving(self, steps, stream, sync_obs):
+        """Steps with translating circles: hold, or one advancing step that commits
+        each pending target only after synchronization and a zero ``no_donor`` count
+        (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()]
+        try:
+            advancing = tables.upload_records(circles, stream)
+            self.stepper._step(steps, action_gpu=self.bodies.action_gpu,
+                               obs_gpu=self.bodies.obs_gpu, stream=stream, advance=advancing)
+            self.bodies.increment_obs_steps(steps)
+            if sync_obs:
+                self.bodies.download_obs_full_async(stream)
+            stream.synchronize()
+            if advancing:
+                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")
+                for motion, _ in circles:
+                    motion.commit()
+                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
+                motion.fail(exc)
+            raise
+
     def step(self, n: int = 1):
         """Advance *n* steps (convenience for interactive use)."""
         self.run(n)
@@ -745,8 +820,16 @@ class Simulation:
         self.field.download_ddf(step_id=self.stepper.step_count)
         return self.field.ddf.copy()
 
+    @context_bound
     def get_flags(self) -> np.ndarray:
-        return self.field.flag.copy()
+        """Flag array; read back from the device when translating circles may have moved."""
+        if getattr(self.stepper, "moving", None) is None:
+            return self.field.flag.copy()
+        self._require_usable()
+        self.stream.synchronize()
+        flags = np.empty_like(self.field.flag)
+        cuda.memcpy_dtoh(flags, self.field.flag_gpu)
+        return flags
 
     # -- Runtime parameter updates -------------------------------------------
     @context_bound
@@ -807,6 +890,7 @@ class Simulation:
                 cuda.Context.synchronize()
                 # Release owned GPU resources while their context is current.
                 self.stepper = None
+                self._moving_tables = None
                 self.bodies = None
                 self.field = None
                 self._module = None

tests/unit/test_duetsch_protocol.py

@@ -0,0 +1,191 @@
+"""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)
+
+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_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
+    import sys
+    sys.path.insert(0, str(Path(__file__).resolve().parents[2] / "research/duetsch_re100_kc5"))
+    try:
+        import run
+    finally:
+        sys.path.pop(0)
+    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():
+    import sys
+    sys.path.insert(0, str(Path(__file__).resolve().parents[2] / "research/duetsch_re100_kc5"))
+    try:
+        import analyze
+    finally:
+        sys.path.pop(0)
+    return analyze
+
+
+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

@@ -117,7 +117,8 @@ def moving_simulation(monkeypatch, *, moving=True):
                                  zero_obs_async=lambda s: writes.append("obs"),
                                  increment_obs_steps=lambda n: writes.append("count"),
                                  download_obs_full_async=lambda s: writes.append("download"))
-    sim.stepper = SimpleNamespace(step=lambda n, **k: writes.append(("step", n)), blocked=None)
+    sim.stepper = SimpleNamespace(step=lambda n, **k: writes.append(("step", n)), blocked=None,
+                                  moving=None)
     sim._assert_runtime_contracts = lambda: None
     return sim, writes
 
@@ -216,8 +217,11 @@ def test_add_body_and_set_body_route_targets(monkeypatch):
     bid = sim.add_body("circle", (50., 40.), 10., translation=Translation((-2., 2., -3., 3.), .1))
     still = sim.add_body("circle", (100., 40.), 5.)
     motion = sim.bodies.get(bid).motion
-    assert sim.stepper.blocked is not None  # declaring a moving circle blocks direct steps
+    sim._moving_tables = "tables"
+    sim._attach_moving()
+    assert sim.stepper.blocked is None
     sim.set_body(bid, offset=(0., .3), vy=.05)
+    assert sim.stepper.blocked == module.PENDING_BLOCKED  # direct steps never consume a target
     sim.set_body(bid, offset=(0., .4), vy=.05)  # last wins
     for kwargs in (dict(offset=(0., .6)), dict(vy=.05), dict(vx=.05),
                    dict(offset=(0., .1), vy=.2), dict(offset=(0., 0.), vy=.01)):
@@ -240,6 +244,186 @@ def test_add_body_and_set_body_route_targets(monkeypatch):
                          **kwargs)
 
 
+class Tables:
+    """MovingTables double: records uploads, steps and diagnostics reads."""
+    DIAG = ("no_donor",)
+
+    def __init__(self, log, no_donor=0., fail=None):
+        self.log, self.no_donor, self.fail = log, no_donor, fail
+
+    def upload_records(self, circles, stream):
+        records = [m.record() for m, _ in circles]
+        self.log.append(("upload", records))
+        uploads = sum(1 for w in self.log if isinstance(w, tuple) and w[0] == "upload")
+        if self.fail == "upload" or (self.fail == "hold_upload" and uploads == 2):
+            raise RuntimeError("upload failed")
+        return any(r[0] != r[1] for r in records)
+
+    def zero_diagnostics(self, stream):
+        self.log.append("zero_diag")
+
+    def download_diagnostics(self):
+        return np.array([[self.no_donor]], np.float32)
+
+
+def advancing_simulation(monkeypatch, **tables):
+    sim, writes = moving_simulation(monkeypatch)
+    sim.bodies.objects[0].motion = circle()
+    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)  # real public entry
+    sim._moving_tables = sim.stepper.moving
+    return sim, writes, sim.bodies.objects[0].motion
+
+
+def test_run_one_advances_and_commits_after_sync(monkeypatch):
+    sim, writes, motion = advancing_simulation(monkeypatch)
+    motion.set_target((0., .25), (0., .05))
+    with pytest.raises(ValueError, match="run\\(1\\)"):  # multi-step with a target
+        sim.run(2, stream=STREAM)
+    assert writes == [] and motion.pending is not None
+    sim.run(0, stream=STREAM)  # zero steps never consume the target
+    assert motion.pending is not None
+    sim.run(1, stream=STREAM)
+    uploads = [w for w in writes if isinstance(w, tuple) and w[0] == "upload"]
+    assert uploads[0][1] == [((50., 40.), (50., 40.25), (0., .05))]
+    assert uploads[1][1] == [((50., 40.25), (50., 40.25), (0., 0.))]  # first hold after move
+    assert ("step", 1, True) in writes
+    assert motion.applied == (0., .25) and motion.pending is None
+    writes.clear()
+    sim.run(3, stream=STREAM)  # hold steps: no advance
+    assert ("step", 3, False) in writes and motion.applied == (0., .25)
+
+
+@pytest.mark.parametrize("failure", ("upload", "first_sync", "no_donor", "hold_upload",
+                                     "hold_sync"))
+def test_failed_advance_fails_simulation(monkeypatch, failure):
+    kwargs = dict(no_donor=1.) if failure == "no_donor" else dict(fail=failure)
+    sim, writes, motion = advancing_simulation(monkeypatch, **kwargs)
+    motion._on_failure = sim._fail_moving
+    motion.set_target((0., .25), (0., .05))
+    syncs = []
+
+    def synchronize():  # stands in for an asynchronous failure surfacing at sync
+        syncs.append(1)
+        if (failure, len(syncs)) in (("first_sync", 1), ("hold_sync", 2)):
+            raise RuntimeError("sync failed")
+    with pytest.raises(RuntimeError):
+        sim.run(1, stream=SimpleNamespace(synchronize=synchronize))
+    # Before the step-group sync and no_donor check nothing is committed; the hold
+    # record is staged after the commit (applied already matches the moved device).
+    committed = failure in ("hold_upload", "hold_sync")
+    assert motion.applied == ((0., .25) if committed else (0., 0.))
+    assert (motion.pending is None) == committed
+    assert motion.failed and sim._failed is not None
+    assert sim.stepper.blocked == sim._failed  # direct stepper reuse refused too
+    with pytest.raises(RuntimeError, match="Simulation failed"):
+        sim.run(1, stream=STREAM)
+
+
+def test_failed_hold_step_blocks_the_direct_stepper(monkeypatch):
+    # No target is pending, so only the failure latch can refuse direct steps.
+    from CelerisLab.lbm.stepper import LBMStepper
+    sim, writes, motion = advancing_simulation(monkeypatch, fail="upload")
+    motion._on_failure = sim._fail_moving
+    launches = []
+    kernels = SimpleNamespace(get_function=lambda name: lambda *a, **k: launches.append(name))
+    cfg = SimpleNamespace(streaming="double_buffer", threads_per_block=128, is_d2q9=True,
+                          nx=64, ny=32, nz=1)
+    field = SimpleNamespace(ddf_gpu="ddf", temp_gpu="temp", flag_gpu="flag", n_curved=0)
+    stepper = LBMStepper(field, kernels, cfg)
+    stepper.moving = sim.stepper.moving
+    sim.stepper = stepper
+    assert motion.pending is None and stepper.blocked is None
+    with pytest.raises(RuntimeError, match="upload failed"):
+        sim.run(2, stream=STREAM)  # hold step
+    assert sim._failed is not None and stepper.blocked == sim._failed
+    with pytest.raises(RuntimeError, match="moving circle step failed"):
+        stepper.step(1, action_gpu=None, obs_gpu=None)
+    assert launches == [] and stepper.step_count == 0
+    assert (field.ddf_gpu, field.temp_gpu) == ("ddf", "temp")
+
+
+def test_moving_tables_layout_and_attachment(monkeypatch):
+    from CelerisLab.body.motion import MovingTables
+    a = MovingCircle((50., 40.), 10., Translation((-1., 1., -3., 3.), .1))
+    b = MovingCircle((90.25, 30.5), 5., Translation((0., 0., 0., 0.), .1))
+    band, geom = MovingTables.layout([(a, 2), (b, 5)])
+    wa, ha = 63 - 37 + 1, 55 - 25 + 1
+    assert band == [37, 25, wa, ha, 0, 2, 83, 23, 16, 16, wa * ha * 8, 5]
+    assert geom == [50., 40., 50., 40., 10., 0., 0., 90.25, 30.5, 90.25, 30.5, 5., 0., 0.]  # hold
+    sim = routed_simulation(monkeypatch)
+    sim.add_body("circle", (50., 40.), 10., translation=Translation((0., 0., 0., 0.), .1))
+    sim._moving_tables = "tables"
+    sim._attach_moving()
+    assert sim.stepper.moving == "tables" and sim.stepper.blocked is None
+    sim._moving_tables = None
+    sim._attach_moving()
+    assert sim.stepper.moving is None and sim.stepper.blocked is not None
+
+
+class CountingOwner(Owner):
+    """Owner double that counts nested activations like CudaContext."""
+    def __init__(self):
+        super().__init__()
+        self.active = 0
+
+    def activate(self):
+        from contextlib import contextmanager
+
+        @contextmanager
+        def scope():
+            self.active += 1
+            try:
+                yield self
+            finally:
+                self.active -= 1
+        return scope()
+
+
+def test_public_diagnostics_keep_owner_active_through_sync_and_download(monkeypatch):
+    sim = make_simulation(monkeypatch)
+    sim.ctx, sim._failed = CountingOwner(), None
+    seen = []
+    sim._stream = SimpleNamespace(synchronize=lambda: seen.append(("sync", sim.ctx.active)))
+    expected = np.arange(6, dtype=np.float32).reshape(1, 6)
+
+    def download():
+        seen.append(("download", sim.ctx.active))
+        return expected
+    sim.stepper = SimpleNamespace(moving=SimpleNamespace(download_diagnostics=download))
+    assert sim.moving_diagnostics() is expected
+    assert [name for name, _ in seen] == ["sync", "download"]
+    assert all(active > 0 for _, active in seen) and sim.ctx.active == 0
+
+
+def test_stepper_launches_moving_links_before_static_curved():
+    from CelerisLab.lbm.stepper import LBMStepper
+    launches = []
+    module = SimpleNamespace(get_function=lambda name: lambda *a, **k: launches.append((name, a, k)))
+    cfg = SimpleNamespace(streaming="double_buffer", threads_per_block=128, is_d2q9=True,
+                          nx=64, ny=32, nz=1)
+    curved = SimpleNamespace(**{f"{name}_gpu": name for name in (
+        "fluid_idx", "dir", "q", "rx", "ry", "rz", "fallback_class", "body_id")},
+        surface_velocity=None, surface_velocity_gpu=None, link_exchange_gpu=None)
+    field = SimpleNamespace(ddf_gpu="ddf", temp_gpu="temp", flag_gpu="flag", n_curved=5,
+                            curved=curved, invalidate_host_ddf_cache=lambda: None)
+    stepper = LBMStepper(field, module, cfg)
+    stepper._launch_force_region = stepper._launch_sensor = lambda *a, **k: None
+    stepper.moving = SimpleNamespace(band_gpu="band", geom_gpu="geom", diag_gpu="diag",
+                                     slot_debug_gpu=None, n_bodies=1, n_slots=300)
+    stepper.step(2, action_gpu="act", obs_gpu="obs")
+    names = [name for name, _, _ in launches]
+    assert names == ["MovingCircleLinks", "CurvedBoundaryKernel", "OneStep"] * 2
+    assert stepper.step_count == 2
+    assert launches[1][1][0] == "ddf" and launches[3][1][0] == "temp"  # swap between steps
+    args, kw = launches[0][1], launches[0][2]
+    assert args[:4] == ("ddf", "flag", "band", "geom") and int(args[5]) == 300
+    assert kw["grid"] == (3, 1, 1) and kw["block"] == (128, 1, 1)
+
+
 def test_recompile_installs_a_guarded_stepper(monkeypatch):
     from CelerisLab.lbm.stepper import LBMStepper
     sim = routed_simulation(monkeypatch)

tests/unit/test_state_regressions.py

@@ -297,7 +297,7 @@ def test_passive_overlap_failure_keeps_live_state(device, entry, second_x):
                           stream=SimpleNamespace(synchronize=device.synchronize),
                           stepper=SimpleNamespace(blocked=None))
     for name in ("_require_usable", "_moving_bodies", "_reject_moving_bodies",
-                 "_validate_moving_bodies"):
+                 "_validate_moving_bodies", "_attach_moving"):
         setattr(sim, name, getattr(Simulation, name).__get__(sim))
     with pytest.raises(ValueError, match="active fluid receivers|unique cut-link writers"):
         if entry == "sync":

tests/unit/test_viv_logic.py

@@ -0,0 +1,96 @@
+"""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)
+
+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
+    import sys
+    sys.path.insert(0, str(Path(__file__).resolve().parents[2] / "research/viv_re150"))
+    try:
+        import run
+    finally:
+        sys.path.pop(0)
+    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:])