Use committed Taylor body IDs after topology sync

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

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

README.md

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

docs/continuous_viv_video.md

@@ -0,0 +1,86 @@
+# Continuous VIV videos
+
+Research-only, two separate settings, each one uninterrupted OFF -> ON -> OFF
+fluid/oscillator trajectory. No production changes, topology switching, moving
+checkpoint feature, retraining, force filtering, or parameter search.
+
+The user selected separate videos because runtime topology changes are unsupported
+with translating circles. This is deliberately different from the old Legacy
+multi-setting video, which rebuilt solvers and transferred full state.
+
+## Frozen physical and control contract
+
+- Source baseline `73d3f6214e0d28ced4744665aa3bb05522aa2fa8`; D40, U=.005,
+  nu=.004, MRT/FP32/unshifted double-buffer. Same parabolic regularized inlet,
+  NEQ outlet and bounce-back walls as the retained m*=4 pair.
+- Common transverse assembly: m*=4 (all three disks), fn D/U=.2, zeta=0.
+  Structure uses previous native-step sum CME_y+Wen_y. No policy-force averaging
+  is applied to this RHS. Same +/-D travel, speed, finite-state and no-donor stops.
+- Each phase: Karman400 complete SI of3200, steady320 SI of4000 native steps
+  =160 D/U. Steady's SI follows the actual old video1000-step interval under
+  diffusive D40 scaling, rather than borrowing PPO's800-step interval. Total480 D/U.
+  No extension, retuning or replacement if steady motion/control is weak.
+- Steady: canonical new-reproduce constant surface velocities
+  (0,-5.1,+5.1) U. This is NOT a PPO transfer. No model/normalizer is loaded.
+  Existing 2,048,000-step at-rest warmup precedes the free response.
+- Karman: actual pinball Re_D50, historically named `karman_re100`, frozen
+  model `148240c1dcb8b11d8e5a0a9e991f6874c4e799c9744707bf0fb32a8efeb0468d`
+  and norm `534fda818c030a56960e95d0f3ac1e350d887d03db7f80f1d6453195a502b7f0`.
+  Use the verified static release once, at video start (clock4576000).
+- Both switches occur on SI boundaries. The oscillator y/v, previous RHS,
+  both DDF buffers, moving geometry and solver clock are never reset. OFF still
+  gathers sensor/CME observations. First PPO action uses the immediately previous
+  complete OFF interval, unlike an independent arm's initial all-zero observation.
+  Deterministic feed-forward PPO has no recurrent state. EMA resets to zero each
+  SI, as in training/reproduction; original video's SplitAdvance also reset the
+  constant-steady actuation every interval. Physical OFF is exactly zero wall speed.
+  There is no extra transition blend. The first1 D/U after either switch is
+  explicitly labeled a switch transient. ON from an already moving state differs
+  from the previous independent pair; +/-D stops are retained without retuning.
+  A visual ON label means enabled control,
+  not a proven suppression/cloaking result.
+
+## Video provenance and presentation
+
+Read-only source: DynamisLab `src/drl_pinball/reproduce/continuous_video.py`, SHA256
+`3bd9c777171c9e8932c78298accbb6f61a8323b0eed5938e0f56e167b3df0743`.
+Approved original MP4 SHA256
+`578d2cf5113429006b6919bbb8ef8ff6ce65ebcd0c056ac478f8ce2515e4cbab`.
+The actual new-reproduce conversation confirms fixed-position plain-language
+labels, white bars, OFF/ON toggle, rotation sectors, and 1/2/5/10/5/2/1 gears.
+
+Only the display idiom/gear schedule are reused. Legacy population decoding,
+physics, topology handoff and policy runner are not imported. Modern macroscopic
+fields feed centered vorticity with all stencil nodes fluid. Fixed color limits
+[-3,+3], 21 RdBu_r bins; per-frame clipping fraction retained. Box downsampling
+is display-only. Both main field and enlarged inset have fixed lab coordinates;
+the inset shows the equilibrium line. Actual moving centers and applied surface
+speeds drive the body overlays. y/D trace has fixed +/-1 limits. No motion gain.
+
+1280x800,30fps; first80 D/U of each phase uses the original display gear ramps,
+last80 D/U uses normal cadence (800 native steps/frame =.1D/U). Every native
+step and every SI is still executed. Streaming ffmpeg avoids stored PNG sequences.
+Keep raw scalar histories, SI observations/actions, per-frame time/motion records,
+selected keyframes, source/input hashes and terminal status. These are new
+continuous runs, not an exact replay of the old independent paired histories.
+
+## Gates and limits
+
+CPU: switch-state/clock continuity, first ON observation, exactly zero OFF action,
+SI timing despite sampling, physical scaling, schedule monotonicity, vorticity
+sign/mask, fixed-coordinate mapping. Real smoke:32 warmup +3x32 response steps
+on D40 (smoke SI32), all switches/encoder; not a physical comparison. Grok reviews
+the exact research diff before full runs. Capture no-donor, finite state, native
+RHS and all stop outcomes. Report per-phase late-half motion statistics and
+raw-force limits; appearance and successful duration are not physics qualification.
+
+Example (outputs must be new, outside checkout; select a safe GPU):
+
+```sh
+CUDA_VISIBLE_DEVICES=0 PYTHONPATH=src python research/continuous_viv/run.py --scene steady --out /tmp/celeris-results-continuous-viv/steady
+CUDA_VISIBLE_DEVICES=0 PYTHONPATH=src python research/continuous_viv/run.py --scene karman --model /tmp/celeris-karman-viv-assets/model.zip --norm /tmp/celeris-karman-viv-assets/norm.json --release-dir /tmp/celeris-results-karman-viv-516975e/pair --out /tmp/celeris-results-continuous-viv/karman
+```
+
+Append `--smoke` and use separate output directories for the bounded gates.
+The Karman smoke does not read/install the physical release. Source and review
+SHA, interpreter/version and GPU resource readings belong in the execution receipt.

docs/fig2_collection.md

@@ -0,0 +1,109 @@
+# Fig. 2 material collection: modern moving plant
+
+Research-only child of `76873d911c4e2caba1bacc25aaa1495ea768453f`.
+No production/kernel/manuscript or DynamisLab edits. These are fresh common-release
+paired experiments, separate from the continuous OFF/ON/OFF demonstration videos.
+
+## Fixed plant and cases
+
+All cases use D40, Uref=.005, nu=.004, MRT/FP32/unshifted double buffering,
+the existing mean-parabolic/regularized inlet, NEQ outlet and bounce-back walls.
+One transverse assembly DOF: m*=4 based on all three disks, fnD/U=.2, zeta=0.
+No force filter or retuning. Existing +/-D, speed.025, nonfinite and no-donor stops.
+Structure RHS is previous native-step sum of pinball CME_y+Wen_y. Policy uses
+SI-mean packed CME and fixed lab probe signals; these are different force paths.
+SI=3200; EMA alpha=.9**.25 resets every SI as in the acquisition contract.
+Steady video used SI4000; its descriptive numbers are not the new collection.
+
+| case | prelude/release clock | controller | collection |
+|---|---|---|---|
+| steady | new static pinball, 2,048,000 steps | constant wall speeds (0,-5.1,5.1)U, no PPO | 400 SI,160 D/U |
+| Karman | verified saved static release at4,576,000 | historical `karman_re100` PPO; actual pinball Re_D50, scale8/bias(0,-4,4) | 400 SI,160 D/U |
+| Taylor y0 | clean512,000 then static pinball512,000; clock1,024,000 | `vortex_taylor` PPO, scale4/samebias | 150 SI,60 D/U |
+
+Zero/control share exact current/spare populations, flags, clock and zero initial
+structural state. Target is a separate no-pinball reference at the SAME absolute
+clocks, not a copy of a solid-containing release. Its different geometry is the
+definition of the reference, not identical microscopic initial conditions.
+Karman target must run to4,576,000; old target.npy is from earlier clocks and is
+not reused for same-clock field errors. Model/norm files are verified by SHA256
+in the driver and written with all source hashes into case.json.
+
+## Taylor injection
+
+Exactly one event at release in every role, inside the owning synchronized CUDA
+context: centre(600,510.5), R80, strength **.000075**. Reuse the existing initializer
+unchanged. It rewrites ALL interior cells to equilibrium and uploads BOTH buffers;
+this removes existing nonequilibrium globally. The impulse and first two SI of
+native force/displacement remain in the record. This is a declared limitation,
+not an independently qualified physical vortex generator.
+
+Strength has inverse-time units: u_theta=(s*pi/2)r exp(-r²/R²). D20 R40/U.01/s.0003
+maps to D40 R80/U.005/s.000075 under diffusive scaling. Using .03*U at D40 would
+double u/U and quadruple p/U²; actual-function CPU positive/negative controls lock it.
+
+Current Horsea `acquisition.py` SHA
+`ea674080501407d00ded6d23467423c516cd37e04d166d5adda16465100204d3`
+defines EIGHT offsets(-20,-15,-10,-5,0,5,10,15). `legacy_test/acquire.py` SHA
+`91cd29a368dfc53a0804ec37397d2e5f25f0f8721f01b560a0f4328b8fd1dd9c`
+selects the controlled centre-uy event: median first20, peak over38..87 inclusive.
+All roles use that common event index; own target/role peak delays are separately
+reported. The earlier five-offset remote summary is superseded by actual source.
+
+## Estimands and storage
+
+Primary: unaligned same-clock/event velocity-vector errors, Eq1
+`sqrt(mean_ROI((dux²+duy²)/Uref²))`. No target-energy denominator. Lab ROI is
+front-equilibrium-centre-relative x/D=[4,24], y/D=[-5,5], equivalent to sensor-plane
+x/D=[-6,14]; D40 inclusive columns1360..2160 and rows311..710. All ROI cells must
+be fluid and finite in every compared field; record support coverage.
+
+Motion uses native last80 D/U for steady/Karman, whole60 D/U for Taylor; mean,
+centered RMS and half-range remain distinct. Karman complete-cycle means and
+eight independently phased nearest snapshots reuse actual acquisition helpers,
+centre-uy column3 and minimum crossing gap. Missing crossings yield unavailable
+phase diagnostics, never fabricated phase slots. Steady retains its last snapshot
+and separately labeled fixed-window mean. Field collection for steady/Karman is
+the last160 SI boundaries (96.4..160 D/U); motion window remains80..160 D/U.
+Taylor retains its eight common event snapshots. Fixed-window mean error differs
+from mean instantaneous error; both are labeled. Native six-channel DTW is a
+supplementary training-space metric, never field recovery or physical delay;
+its actual legacy lag uses flat sensor column1, unlike the phase helper's column3.
+
+Use one FP32 full-grid RAM cache per role, plus compact ROI RAM histories for
+same-clock comparisons. Admit only with >=16GiB MemAvailable; no3.35GB disk field
+cache. Save only selected fields, means, explicit masks, scalars and source receipts.
+Nonfluid placeholders are zero-filled for helpers but always masked in plots and
+scores; nonfinite FLUID is a hard failure. Serialize compilation/GPU work. Verify
+compressed native histories before removing only redundant task scratch copies.
+
+Plot standards: vorticity curl(u)D/U, five-point fluid stencil, common RdBu_r21bins
+with fixed[-3,3] limits and nearest interpolation; lab-fixed coordinates and ROI.
+Mean plots exclude cells that were not fluid at every collected boundary. Selected
+snapshots use actual step geometry. Raw arrays and scalar errors are not clipped.
+
+## V5 boundary
+
+Standard V5 `kar_re100` seeds41..45 DO exist (actual pinball Re_D50), but use
+2000x600 uniform/free-slip and different absolute geometry/normalization/action
+contracts. These are not silently transferred onto the present Legacy-contract
+plant. No verified retained trained V5 steady/Taylor bundle was found; a vortex
+environment prototype is not a trained policy. Existing V5 within-case results
+may be indexed separately; cross-generation scores cannot imply an identical plant.
+
+## Execution
+
+First run focused CPU contracts, independent source review, then each scene with
+`--smoke` (three32-step SI per role; no physical scores). Only the reviewed exact
+SHA may run physical cases, one at a time, new outputs outside the checkout:
+
+```
+CUDA_VISIBLE_DEVICES=3 PYTHONPATH=src python research/fig2_collection/run.py --scene karman --out /tmp/celeris-results-fig2-collection/karman
+python research/fig2_collection/analyze.py /tmp/celeris-results-fig2-collection/karman
+python research/fig2_collection/plot.py /tmp/celeris-results-fig2-collection/karman
+```
+
+Use the same invocation for steady/taylor with new output directories. No automatic
+extension or parameter sweep. All results are developmental material, not asserted
+settled/physically qualified. Preserve partial failures and distinguish missing
+evidence from passing gates.

docs/karman_viv.md

@@ -0,0 +1,120 @@
+# Re_D = 50 Kármán inflow: paired pinball VIV / frozen PPO demo
+
+Development experiment, frozen before results, 2026-10-31. No new production
+method, structural damping, force filter, policy training, or main merge.
+[Force-spike research index](research/moving_boundary_force_spikes.md).
+
+## Plant and structure
+
+Base `e9213dd2cc9cfed0e8487a3cf4ebdd94d23dc64d`; fixed-policy contract from
+`580397eaf83c85b233495791ecc51662a374a247`. Reuse `steady_pinball` config,
+coordinate map and exact held-force oscillator. D40, U=.005, nu=.004, MRT,
+2559 x 1022, FP32 unshifted double buffers. Regularized parabolic inlet,
+NEQ extrapolation outlet, bounce-back y walls (.5,1020.5); no BC changes.
+D20 coordinates map x'=2x, y'=.5+2(y-.5).
+Upstream cylinder (400,510.5), radius40, is static. Pinballs radius20 at
+(1200,510.5), (1252,540.5), (1252,480.5) share only transverse translation.
+Lab probes radius10 at (1600,590.5), (1600,510.5), (1600,430.5) stay fixed.
+
+One shared `Translation`, travel ±D, max wall speed .025, existing clearance
+and per-step move guards unchanged. M2: m*=2 based on all three disks,
+M=7539.8223686, fn=.000025, K=.000186037660, C=zeta=0; release y=v=0.
+Each native step uses previous-step sum(CME_y+Wen_y), then exact oscillator
+update with held force. Wall speed is interval displacement y[n+1]-y[n].
+Cover/uncover impulses are logged, never added to the selected RHS. This
+estimator remains provisional; existing spike/convergence limitations remain.
+
+## Same release state, two clocks
+
+A static precursor retains the accepted Kármán sequence: upstream+probes,
+2,048,000-step warmup, 150 x 3,200-step no-pinball target samples, add static
+pinballs/sync, second 2,048,000-step warmup. Save both DDF buffers, flags,
+clock and last instantaneous force. Construct each moving arm fresh at the
+same zero-offset geometry, check flags, copy both buffers and clock, verify
+exact copy hashes. This one-time research initialization is not a moving
+checkpoint or FSI rollback. Initial policy observation and spin are zero.
+
+Frozen model SHA256 `148240c1dcb8b11d8e5a0a9e991f6874c4e799c9744707bf0fb32a8efeb0468d`;
+normalizer `534fda818c030a56960e95d0f3ac1e350d887d03db7f80f1d6453195a502b7f0`.
+Only these explicit supplied assets are accepted. Policy is CPU deterministic.
+Native policy signals: packed CME (same `read_force` path), lab area-mean
+velocity, each averaged over SI3200. Convert to original D20 training units:
+forces x2, area-mean velocity x2 x78, then frozen force-first normalization
+and clipping. Fine-grid probe pixel footprint is not exactly identical;
+using the original78 cells avoids a hidden area-count normalization change.
+The moving plant, motion-induced force, lab-frame wake displacement and D40
+resolution remain out of training distribution. No transfer guarantee.
+
+A=(8a+(0,-4,4))*.005; EMA alpha=.9^(1/4), reset0 each SI; omega=+A_eff/20.
+This preserves historical physical actuator timescale under dt/4.
+The zero arm substitutes a=(0,.5,-.5), so all spins are identically zero.
+Policy interval averaging does **not** average the native structural RHS.
+
+## Bound, output, and interpretation
+
+One serial pair, zero then controlled: 400SI=1,280,000 native steps=160D/U
+=32Tn per arm. No retuning/automatic extension. Stop on nonfinite state,
+no_donor, travel/speed/move rejection, or four-hour whole-pair wall cap.
+A failure stops the remaining run; retain valid rows and completed intervals.
+NaN trailing rows are unexecuted, not zero-valued observations.
+
+Report first/last80D/U windows, mean offset separately from centered RMS,
+half-range, peak, frequency and drift; native versus SI-averaged force,
+CME/Wen/node/crossing terms and action/observation clipping. Last80D/U is
+the principal descriptive window, not an asserted settled window. If zero
+arm half-range <.1D there, classify weak response and do not claim vibration
+suppression. Report continuous ratios regardless. Cloaking uses retained
+canonical last30 lagged sensor DTW against this precursor's no-pinball target,
+with fixed-body D20 results contextual only. These are development criteria,
+not paper uncertainty, grid convergence, or proof of validated FSI.
+
+Save native18-column rows (~184MB/arm), small SI signals, target, config,
+release two buffers (~188MiB), hashes and terminal reports outside checkout.
+No per-step field output. Expected new data <650MB, check actual resources.
+
+```
+CUDA_VISIBLE_DEVICES=0 PYTHONPATH=src python research/karman_viv/run.py \
+  --model /tmp/celeris-karman-viv-assets/model.zip \
+  --norm /tmp/celeris-karman-viv-assets/norm.json --out /tmp/NEW_OUTPUT
+```
+
+First `--smoke`: same D40 geometry, warmup32/target2x8/warmup32,
+4x8 response steps each arm; strictly plumbing only. CPU source checks and
+independent exact-tip review precede GPU. Preserve original DynamisLab files.
+
+## One predeclared m*=4 exploration (2026-10-31)
+
+The m*=2 attempt at `516975ef6482f983ace347cb10de7450d2fbd680` stopped
+at the declared travel guard in both arms: zero 10.469375 D/U, PPO 27.5525 D/U.
+Its principal 80–160 D/U paired window is unavailable. Early attenuation is
+not sustained suppression; control-only DTW is not a paired cloaking result.
+Those results remain unchanged and retained.
+
+The sole next candidate doubles M and K together: m*=4, M=15079.644737231,
+K=0.000372075320164, fn=0.000025, C=zeta=0. This increases inertia and stiffness,
+not damping. Half response under an identical prescribed external force is
+the selection rationale only; coupled flow feedback can invalidate it. This
+is an exploratory development case, not independent validation or a sweep.
+
+Reuse the exact saved static release, config, target, model and normalizer.
+`--release-dir` checks precursor metadata, config equality, both DDF/flags
+hashes and clock against the original transfer receipt. Each fresh arm still
+checks its upload against that seed. No warmup, moving restart or new release
+file is generated. All fluid/geometry/EMA/SI/RHS settings, 400SI per arm,
+zero damping, ±D guard and other stop conditions above remain unchanged.
+The original failure stops the remaining run rule also remains unchanged.
+
+```
+CUDA_VISIBLE_DEVICES=0 PYTHONPATH=src python research/karman_viv/run.py \
+  --model /tmp/celeris-karman-viv-assets/model.zip \
+  --norm /tmp/celeris-karman-viv-assets/norm.json --out /tmp/NEW_M4_OUTPUT \
+  --mass-ratio 4 --release-dir /tmp/celeris-results-karman-viv-516975e/pair
+```
+
+Report matched valid windows only, separately retaining mean offset,
+centered RMS, half-range, interval/native forces, crossing/diagnostic counts
+and paired DTW only where both arms have sufficient data. The larger
+controlled force excursions remain unresolved: retained m*=2 data strongly
+associate them with the post-action-reset ramp, not exclusively current-step
+node switching. No force smoothing, action-contract change, enlarged travel
+band or further parameter adjustment follows an unsuccessful outcome.

docs/lifecycle.md

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

docs/moving_bodies.md

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

docs/research/moving_boundary_force_spikes.md

@@ -0,0 +1,99 @@
+# Moving-boundary force spikes: evidence and source index
+
+Updated 2026-10-31.
+
+This note indexes supplied audited observations and primary-source extractions.
+It does not report a newly measured cure or implement an alternative method.
+Paper details below were supplied by the parent investigation; no new source
+retrieval or independent paper verification was performed for this note.
+
+## Current implementation and interpretation
+
+- Celeris M1 reconstructs a fresh node as `feq(rho_bar, u_wall) + donor NEQ`.
+- In exact arithmetic, the donor non-equilibrium component has zero mass and
+  momentum when formed against the donor's own equilibrium moments.
+- Re-projecting that same reconstruction onto the same wall velocity therefore
+  does not, by itself, change its mass or momentum or supply a cure for spikes.
+- Floating-point residuals should be measured separately from this algebraic fact.
+- Peng's alternative differs upstream: extrapolate PDFs from three donors, then
+  replace selected moments. It is not simply another projection of current M1.
+- PPO interval-force observations and native per-step FSI right-hand-side forces
+  are distinct signals; comparisons must identify which signal is being used.
+
+## Audited run observations supplied to this note
+
+- D40 revision: `e9213dd2cc9cfed0e8487a3cf4ebdd94d23dc64d`.
+- Window: first `97.303 D/U` of physical time.
+- Maximum displacement: `max |y|/D = 0.304893`.
+- Maximum native force-coefficient increment: `max |ΔCy| = 0.66688`.
+- D20 over the same physical-time window: `max |ΔCy| = 1.80735`.
+- The largest D40 event involved 30 nodes.
+- These are supplied audited observations, not measurements made for this note.
+- Time steps and trajectories differ between runs. This comparison establishes
+  neither spatial convergence nor a unique cause for the force spikes.
+
+Evidence locations supplied by the investigation:
+
+- Horsea: `/tmp/celeris-results-pinball-resolution`.
+- Mac: `validation-evidence/pinball-force-diagnosis`.
+
+These locations index run evidence, not local copies of the cited papers.
+
+## Tao, Hu, and Guo (2016): time-step decoupling
+
+Computers & Fluids **133**, 1–14. DOI:
+[10.1016/j.compfluid.2016.04.009](https://doi.org/10.1016/j.compfluid.2016.04.009).
+[Publisher abstract](https://www.sciencedirect.com/science/article/pii/S0045793016301116).
+
+- Abstract-level support only: non-equilibrium refill reduces oscillations.
+- With local iteration, the compared refill methods give similar oscillation magnitudes.
+- Time-step decoupling reduces oscillations further in the reported sedimentation,
+  Couette, and Poiseuille examples.
+- The exact decoupling schedule and time-step ratios have not been verified.
+- This abstract does not establish equivalence to an arbitrary force filter.
+- Do not implement a claimed reproduction until the actual methods are available.
+
+## Marson et al. (2021): local iterations
+
+Physical Review E **103**, 053308. DOI:
+[10.1103/PhysRevE.103.053308](https://doi.org/10.1103/PhysRevE.103.053308).
+[Open repository record](https://archive-ouverte.unige.ch/unige:151868).
+
+- §VI, p. 18 discusses pressure waves associated with fresh nodes.
+- Appendix E.3, p. 25 and Fig. 18, p. 24 describe the local procedure.
+- Copy a local sublattice, perform local collision and streaming, write the fresh
+  nodes, and repeat before the global collision and streaming operation.
+- One to three iterations performed best in the reported tests; this is not a
+  universal iteration-count recommendation.
+- The reported smoothed torque-noise assessment is qualitative here.
+- Its thin-shell setting differs from the current moving-body configuration.
+
+## Peng (2016): extrapolation and moment replacement
+
+Computers & Mathematics with Applications **72**, 349–374. DOI:
+[10.1016/j.camwa.2015.08.027](https://doi.org/10.1016/j.camwa.2015.08.027).
+
+- §5.1.2, Eq. (50), p. 368: extrapolate PDFs along an outward normal direction:
+  `f_new(x) = 3 f(x + e) - 3 f(x + 2e) + f(x + 3e)`.
+- §5.1.3, Eq. (51): wall equilibrium plus donor non-equilibrium is the family
+  containing the current M1 reconstruction.
+- §5.1.4, Eqs. (52)–(54): in MRT space, replace momentum with `j = rho0 u_wall`,
+  preserve the other moments, then apply the inverse transform.
+- A compressible Celeris adaptation must choose the appropriate density;
+  copying the paper's `rho0` blindly is not justified.
+- Three valid donors require new eligibility checks and fallback behavior beyond
+  the current one-link shared-body clearance rule.
+- §§5.2–5.3 retain nonzero fluctuations in specific configurations; these results
+  do not promise elimination of force spikes in the present problem.
+
+## Bounded next investigation
+
+1. Retain current M1 as the baseline and compare matched prescribed trajectories
+   over matched physical-time windows.
+2. Record raw force RHS, node crossings, donor moments, mass and momentum changes,
+   impulse, spectra, and trajectory response together.
+3. Evaluate local iterations, Peng-style reconstruction, and force filtering as
+   separate interventions, with explicit definitions and diagnostics.
+4. Obtain the actual time-step-decoupling methods before implementing that method.
+
+The current task indexes evidence only; it makes no alternative-method change.

docs/steady_pinball_fiv.md

@@ -0,0 +1,150 @@
+# Uncontrolled steady pinball: development protocol
+
+Frozen before results: 1280 x 512 D2Q9/FP32/MRT, unshifted double buffer;
+rho=1, mean U=.01, nu=.004, D=20, Re_D=50. Centres are (600,255.5),
+(626,270.5), (626,240.5), radius10. No upstream cylinder. Passive probes stay
+at (800,295.5), (800,255.5), (800,215.5), radius5. Every omega is zero.
+Parabolic inlet, regularized NEQ damping1; y bounce-back planes .5/510.5;
+NEQ-extrapolated outlet with backflow clamp (the configured blend is inactive in
+this MRT route). Initial fluid is
+parabolic equilibrium. Fixed warmup is 4 NX/U=512000 steps.
+
+Geometry and fluid parameters follow DynamisLab 609116ec26f19766f5e9ce2b0df427d4e24faae8,
+`legacy_test/core/legacy_env_builder.py:build_steady_cloak` and the retained
+`config_lbm_pinball_legacy_re50d_u001.json`. This is a modern Celeris reconstruction,
+not a bitwise Legacy reproduction: Legacy uses shifted incompressible MRT, a
+different curved interpolation and prescribed parabolic outlet velocity. Current
+Celeris uses standard MRT, two-branch BFL and extrapolated outlet velocity.
+The fixed-hold control uses the SAME modern settings as the free case.
+
+## Shared translation and ownership
+
+`Translation(..., shared=True)` is explicitly transverse only. Reusing one
+envelope INSTANCE declares common motion; equal but separate envelopes are not
+linked. Before initialization and every positive Simulation.run, group members
+must have identical applied offsets, effective target offsets and interval wall
+velocities. Stage each member with set_body, then call run(1) once. An incomplete
+or inconsistent target group is rejected before device uploads and can be fixed.
+Existing pending-target direct-stepper restrictions remain.
+
+Only members of that group may have overlapping whole swept bands. Their fixed
+surface clearance must exceed4 lattice cells. This conservative bound exceeds
+2 sqrt(2)+MAX_STEP = 3.329: M1 changed cells lie in old/new disks and read one-link
+neighbors; M2 reads at most two links outside its new disk. Old/new relative
+displacement is at most .5. Thus reads cannot reach another member's old/new
+solid or M1 write region. M1 writes are disjoint. M2 writes opposite population
+slots inside its own new disk, also disjoint. Separate first_slot ranges still
+allocate each body's logical band/debug storage. Shared band cells alone do not
+write anything: M1 returns unless membership changes; M2 returns unless a cut
+link exists. Unrelated moving/static body and probe overlaps still reject.
+
+This is a sufficient restricted-group condition, not support for near-contact
+or independent motion. The r4, distance8.3, dy=.5 counterexample remains rejected.
+CPU tests enumerate candidate donors/write footprints over subcell phases;
+GPU tests compare fixed vs hold and member-order permutations through crossings
+at wall speed .01 over the full +/-D envelope (8000 steps). The initial .5-speed
+fixture failed with NaNs; a single-circle control also failed at step69. That
+Mach~.866 fixture was unsuitable as a low-Mach numerical gate. Its failure is
+retained; the geometric MAX_STEP proof is unchanged, not a stability guarantee.
+CUDA kernels and arithmetic are unchanged.
+
+## New structural parameters (not inherited from the paper)
+
+One transverse state y,v; total assembly M=10 rho (3 pi D^2/4)=9424.77796,
+fn=U/(5D)=.0001, Tn=10000, zeta=.01, K=M(2 pi fn)^2,
+C=2 zeta M(2 pi fn). Initial y=v=0, no seed or later kick. Force is the sum
+of the three instantaneous CME_y+Wen_y values from the previous fluid step;
+node impulses are logged separately and never added. The damped oscillator is
+integrated exactly for that piecewise-constant force over dt=1. Wall velocity
+is y[n+1]-y[n], not the endpoint oscillator velocity. No filtering or retuning.
+
+Observe600000 steps=60Tn=300 D/U after warmup, in each of free and fixed-hold
+cases. Final20Tn is a descriptive comparison window; ten-Tn blocks expose drift.
+No settled/periodic claim follows merely from reaching the cap. Stop on |y|>D,
+nonfinite values, no_donor, rejected displacement or speed>.05. Report any quiet
+response honestly. No automatic extension. --smoke uses32 warmup+32 steps only.
+
+Native-step rows retain y,v, each body's CME/Wen/node terms and refill counters;
+each600000x15 FP64 file is72MB. No full-field output or policy training. Compile
+and run on one device serially, keep original outputs and dirty trees untouched.
+This stage makes no claim of physical grid convergence or vibration suppression.
+
+Run from checkout with PYTHONPATH=src and the existing Python environment:
+`python research/steady_pinball/run.py --smoke --out <new-external-directory>`
+then separately `--out <new-free-directory>` and `--fixed --out <new-fixed-directory>`.
+
+## Zero-damping development follow-up (predeclared 2026-10-31)
+
+The original d0bcefa run and its zeta=.01 results remain historical evidence.
+The current runner exposes three fixed `--case` choices, all with C=0 as requested:
+
+| Case | M / (rho 3 pi D^2/4) | fn | M | K | Purpose |
+|---|---:|---:|---:|---:|---|
+| Z0 | 10 | .0001 | 9424.77796 | .00372075320 | Remove damping only |
+| M2 | 2 | .0001 | 1884.95559 | .000744150640 | Reduce mass and stiffness together at fixed fn |
+| R2 | 2 | .00006 | 1884.95559 | .000267894231 | Change stiffness only relative to M2 |
+
+R2 deliberately places fn at the retained fixed-body lift peak, f=.00006
+(St=.12 on mean U). This is an informed development choice, not independent
+physical validation or a promise of lock-in. All geometries, fluid settings,
+512000-step fixed warmup, y0=v0=0, force time level and unfiltered RHS are unchanged.
+No run result may select a replacement parameter or an automatic extension.
+
+Each case has a 600000-step cap, the same 300 D/U physical duration. This equals
+60 natural periods for Z0/M2 and 36 for R2. Compare the common final200000 native
+steps, and six100000-step blocks, retaining actual natural/response frequencies
+and ongoing drift. Do not call those blocks ten R2 periods. No settled-VIV claim
+follows from a larger displacement. Report all three outcomes, including quiet
+responses or stopped cases, alongside the existing fixed and damped-free records.
+
+Unchanged stops: |y|>D, invalid move/speed, nonfinite values or no_donor. A 1800s
+per-case wall cap and GPU headroom check bound resource use. Preserve partial rows;
+do not add damping, clip trajectories, filter forces, or relax any stop. Each
+case adds72MB of scalar history, no full fields. The separate Legacy/modern
+baseline diagnosis does not change this FSI protocol.
+
+## Matched D40 resolution comparison (predeclared 2026-10-31)
+
+The first finer case is M2, with exactly the same mass ratio2, reduced velocity5,
+zero structural damping, initial rest state and physical geometry as the retained
+D20 M2. Targeting a development response of0.5--1D does not imply this case will
+reach it. Do not change the structure within this comparison or discard a smaller
+response. D60 is a supported scaling choice, not an automatically scheduled run.
+
+For s=D/20 use dx'=dx/s, dt'=dt/s^2, U'=.01/s, nu'=.004, rho'=1,
+fn'=fn/s^2, M'=s^2 M and K'=K/s^2. The same physical force per unit span
+has F'_lat=F_lat/s. For D40 this gives U=.005, fn=.000025,
+M=7539.82237, K=.000186037660, Re_D=50, and Tn=40000 steps.
+The y walls are at .5 and ny-1.5, while x inlet/outlet source nodes are at0
+and nx-1. Therefore map x'=s*x and y'=.5+s*(y-.5), with
+nx'=s*(1280-1)+1, ny'=s*(512-2)+2. D40 uses2559x1022 and centres
+(1200,510.5),(1252,540.5),(1252,480.5); every radius, probe and envelope scales
+with D. The parabolic profile retains the same normalized wall coordinate.
+This preserves the existing ghost-source boundary convention, without asserting
+that it is an exact continuum open-boundary location.
+
+Warmup=2048000 and response=2400000 native steps at D40: the SAME256 and300
+D/U durations as D20, including60 natural periods after release. Compare the
+six50-D/U blocks and the common final100-D/U window. Report max absolute y/D,
+half-range, RMS, frequency and drift; peak-to-peak is not amplitude. No settling
+claim follows from completion. Keep raw individual CME/Wen/node terms and all
+crossing/fallback counters. Normalize net force by .5*rho*U^2*D in each case.
+Report native crossing/non-crossing force increments separately, plus D40
+four-step block means aligned with D20 intervals (diagnostic only, never replace
+the raw spikes). Report final-window high-frequency variance above f*D/U=1,
+using demeaned FFT Parseval weights, with the same physical cutoff for both cases.
+Different motion histories confound spike-size comparisons; this is simultaneous
+space/time/Mach refinement, not proof of a unique coarse-grid cause.
+
+Only one full D40 M2 is initially scheduled, after CPU checks and a32+32-step
+smoke. Stop on the unchanged geometry/speed/nonfinite/no_donor checks or a7200s
+whole-case wall cap; save completed rows and status. The velocity cap scales to
+.025 so its physical bound is unchanged. No force filtering, trajectory clipping,
+automatic parameter search or extension. Scalar history is288MB; no full fields
+or DDF snapshots are written. Two FP32 population buffers require about188MiB,
+before flags, moving tables and context overhead; measure actual peak usage in
+the smoke and full run. Keep all D20 outcomes. No production source changes.
+
+Commands (one GPU and serial compilation, existing environment):
+`python research/steady_pinball/run.py --case M2 --D 40 --smoke --out <new-smoke-dir>`
+then `python research/steady_pinball/run.py --case M2 --D 40 --wall-seconds 7200 --out <new-full-dir>`.

docs/supported_settings.md

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

docs/validation/moving_bodies_stage1.md

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

docs/validation/viv_re150_predeclare.md

@@ -1,4 +1,6 @@
-# Single-cylinder free transverse VIV — PREDECLARATION (agreed at 1abfe954)
+# Single-cylinder free transverse VIV, Re 150 (Kolahdouz/Qin) — PREDECLARATION (agreed at 1abfe954)
+
+(File renamed from `leontini_viv_predeclare.md`; the Leontini Re 200 case is only the queued second case below.)
 
 Stage goal (Yanqi via TRON): prescribed Dütsch case + one VIV case both
 "reasonably close" to a predeclared reference = development-stage success, then
@@ -27,10 +29,17 @@ Claude, Grok and Astra at 1abfe954; implemented in research/viv_re150/.
 
 - D = 40, MRT, dense double buffer; U_lat = .04, nu = U D / 150 = .010667,
   tau = .532.
-- Box [0, 75 D] × [0, 60 D] (3000 × 2400) with the cylinder centre at
-  (30 D, 30 D): 30 D from the inlet, 45 D from the outlet, centred in y
-  (Kolahdouz §4.1.1). Qin's own box ([0, 40 D] × [0, 20 D], centre (10 D, 10 D))
-  is the target paper's, not ours.
+- Nominal 75 D × 60 D lattice allocation (3000 × 2400), based on the
+  Kolahdouz §4.1.1 layout. Coordinate clarification of the existing run:
+  the driver uses `(30 D, (ny - 1)/2)` = `(1200, 1199.5)` in lattice
+  coordinates, or `(30 D, 29.9875 D)`. The earlier `(30 D, 30 D)` wording
+  described the nominal box, not the exact executed centre.
+  Free-slip reflection planes are at `y = .5` and `y = ny - 1.5 = 2398.5`;
+  their separation is 2398 cells = 59.95 D, with the centre exactly midway,
+  29.975 D from each wall. Node counts are not physical wall-to-wall lengths.
+  This corrects the description, without changing the historical geometry
+  or results. Qin's own box ([0, 40 D] × [0, 20 D], centre (10 D, 10 D)) is
+  the target paper's, not ours.
 - Boundaries: uniform Zou–He inlet, NEQ outlet (rho = 1, backflow clamp),
   free-slip y walls. Kolahdouz uses a traction outlet: declared mismatch.
 - Qin's own VIV box (40 D × 20 D) is not used; a smaller box would be a

research/continuous_viv/render.py

@@ -0,0 +1,220 @@
+"""Sparse scientific frames and the approved new-reproduce presentation idiom.
+
+Frame gears/UI colors/rotation sectors derive from DynamisLab continuous_video.py
+(read-only SHA in docs). Modern velocity decoding and moving coordinates are used.
+All camera limits, vorticity limits, and motion scale are fixed, never following y.
+"""
+from functools import lru_cache
+import json
+from pathlib import Path
+import subprocess
+import time
+
+import numpy as np
+
+D,U,CY=40.,.005,510.5
+SIZE=(1280,800)
+CW,CCW='#8854B8','#138F78'
+MAIN=(0,120,1280,440)
+INSET=(960,492,1240,772)
+MAIN_LIMITS=(6.4,57.6,-6.4,6.4)
+INSET_LIMITS=(28.75,33.75,-2.5,2.5)
+
+
+def frame_schedule(steps,base_interval,accelerated=False):
+    """Original 15-frame gear holds; only readback cadence changes."""
+    if accelerated:
+        ramp=15*(1+2+5)*base_interval
+        if steps<=2*ramp:
+            raise ValueError('too short for both display ramps')
+        blocks=[(s,15*s*base_interval) for s in (1,2,5)]
+        blocks += [(10,steps-2*ramp)]
+        blocks += [(s,15*s*base_interval) for s in (5,2,1)]
+    else:
+        blocks=[(1,steps)]
+    result,done=[],0
+    for speed,duration in blocks:
+        end=done+duration
+        while done<end:
+            done=min(done+speed*base_interval,end)
+            result.append((done,speed))
+    return result
+
+
+def video_schedule(phase_steps,smoke=False):
+    if smoke:
+        phase=frame_schedule(phase_steps,8)
+    else:
+        # First 80 D/U of each phase is fast display; last 80 is normal.
+        half=phase_steps//2
+        phase=frame_schedule(half,800,True)
+        phase += [(half+n,s) for n,s in frame_schedule(phase_steps-half,800)]
+    return {p*phase_steps+n:s for p in range(3) for n,s in phase}
+
+
+def vorticity(ux,uy,flags):
+    """Centered dx=1 curl; mask every stencil touching nonfluid, no filling."""
+    ux,uy=np.asarray(ux,dtype=float),np.asarray(uy,dtype=float)
+    fluid=(flags&1)!=0
+    valid=(fluid[1:-1,1:-1]&fluid[1:-1,:-2]&fluid[1:-1,2:]&
+           fluid[:-2,1:-1]&fluid[2:,1:-1])
+    curl=((uy[1:-1,2:]-uy[1:-1,:-2])-(ux[2:,1:-1]-ux[:-2,1:-1]))/2
+    result=np.full(ux.shape,np.nan)
+    result[1:-1,1:-1]=np.where(valid,curl*D/U,np.nan)
+    return result
+
+
+def pixel(x,y,box,limits):
+    left,top,right,bottom=box
+    x0,x1,y0,y1=limits
+    return left+(x-x0)/(x1-x0)*(right-left),bottom-(y-y0)/(y1-y0)*(bottom-top)
+
+
+@lru_cache(maxsize=12)
+def font(size,bold=False):
+    from PIL import ImageFont
+    name='DejaVuSans-Bold.ttf' if bold else 'DejaVuSans.ttf'
+    return ImageFont.truetype('/usr/share/fonts/truetype/dejavu/'+name,size)
+
+
+def draw_circle(draw,x,y,r,wall_speed,box,limits):
+    cx,cy=pixel(x,y,box,limits)
+    radius=r*(box[2]-box[0])/(limits[1]-limits[0])
+    rect=(cx-radius,cy-radius,cx+radius,cy+radius)
+    draw.ellipse(rect,fill='white',outline='#263D4B',width=2)
+    if wall_speed is None:
+        return
+    fraction=abs(wall_speed)/(12*U)
+    if not np.isfinite(fraction) or fraction>1.00001:
+        raise RuntimeError('applied rotation exceeds fixed display reference')
+    r2=radius*.72
+    inner=(cx-r2,cy-r2,cx+r2,cy+r2)
+    draw.ellipse(inner,fill='#E4E9ED')
+    angle=min(fraction,1)*360
+    start,end=(-90,-90+angle) if wall_speed<0 else (-90-angle,-90)
+    if angle>0:
+        draw.pieslice(inner,start,end,fill=CW if wall_speed<0 else CCW)
+
+
+class Video:
+    def __init__(self,out,scene,bodies,phase_steps,smoke=False,source='unknown'):
+        from PIL import Image
+        from matplotlib import colormaps
+        self.out,self.scene,self.bodies=out,scene,bodies
+        self.schedule=video_schedule(phase_steps,smoke)
+        self.phase_steps=phase_steps
+        self.source=source
+        self.records=[]
+        self.palette=np.rint(colormaps['RdBu_r'](np.linspace(0,1,21))[:,:3]*255).astype('uint8')
+        self.log=open(out/'ffmpeg.log','wb')
+        self.proc=subprocess.Popen(['ffmpeg','-hide_banner','-loglevel','warning','-n',
+            '-f','rawvideo','-pix_fmt','rgb24','-s','1280x800','-r','30','-i','pipe:0',
+            '-an','-c:v','libx264','-threads','4','-preset','fast','-crf','18',
+            '-pix_fmt','yuv420p','-movflags','+faststart',str(out/'continuous_viv.mp4')],
+            stdin=subprocess.PIPE,stderr=self.log)
+        self.closed=False
+        self.seconds=0.
+
+    def capture(self,sim,step,phase,eff,rows):
+        speed=self.schedule.get(step)
+        if speed is None:
+            return
+        from PIL import Image,ImageDraw
+        start=time.monotonic()
+        macro=sim.get_macroscopic()
+        flags=sim.get_flags().reshape(sim.field.ny,sim.field.nx)
+        fluid=(flags&1)!=0
+        if any(not np.isfinite(macro[k][fluid]).all() for k in ('rho','ux','uy')):
+            raise RuntimeError('nonfinite fluid field at video sample')
+        omega=vorticity(macro['ux'],macro['uy'],flags)
+        finite=np.isfinite(omega)
+        index=np.floor((np.clip(np.nan_to_num(omega),-3,3)+3)/6*21).astype(int).clip(0,20)
+        rgb=self.palette[index]
+        rgb[~finite]=[232,235,237]
+        field=Image.fromarray(rgb[::-1])
+        canvas=Image.new('RGB',SIZE,'white')
+        draw=ImageDraw.Draw(canvas)
+        for box,limits in ((MAIN,MAIN_LIMITS),(INSET,INSET_LIMITS)):
+            x0,x1,y0,y1=limits
+            # Lattice cell-centre coordinates -> image pixel-edge coordinates.
+            crop=(x0*D+.5,sim.field.ny-(CY+y1*D+.5),
+                  x1*D+.5,sim.field.ny-(CY+y0*D+.5))
+            picture=field.crop(crop).resize((box[2]-box[0],box[3]-box[1]),Image.Resampling.BOX)
+            canvas.paste(picture,box[:2])
+            if box==INSET:
+                _,eq=pixel(30,0,box,limits)
+                for x in range(box[0],box[2],10):
+                    draw.line((x,eq,min(x+5,box[2]),eq),fill='#657887',width=1)
+            if self.scene=='karman':
+                draw_circle(draw,10,0,1,None,box,limits) if box==MAIN else None
+            for bid,w in zip(self.bodies,eff):
+                obj=sim.bodies.objects[bid]
+                x,y=obj.motion.centre(obj.motion.applied)
+                draw_circle(draw,x/D,(y-CY)/D,obj.radius/D,float(w),box,limits)
+        x0,y0=pixel(INSET_LIMITS[0],INSET_LIMITS[2],MAIN,MAIN_LIMITS)
+        x1,y1=pixel(INSET_LIMITS[1],INSET_LIMITS[3],MAIN,MAIN_LIMITS)
+        draw.rectangle((x0,y1,x1,y0),outline='#647C8C',width=1)
+        draw.rectangle(INSET,outline='#647C8C',width=1)
+        def text(x,y,value,size=16,color='#263D4B',bold=False):
+            draw.text((x,y),value,font=font(size,bold),fill=color)
+        on=phase=='on'
+        text(28,15,'INCOMING FLOW',12,'#667887',True)
+        text(28,39,'Steady inflow' if self.scene=='steady' else 'Periodic cylinder wake',29,bold=True)
+        text(28,83,'One continuous flow and oscillator trajectory | Re_D = 50 | zero structural damping',15)
+        text(850,18,'CYLINDER ROTATION',12,bold=True)
+        text(850,43,'Clockwise',14,CW);text(1020,43,'Counterclockwise',14,CCW)
+        text(850,70,'Sector angle = actual surface speed',13,'#667887')
+        text(850,90,'Full circle = 12 x reference inflow speed',13,'#667887')
+        draw.line((24,453,1256,453),fill='#E5EBEF')
+        text(28,465,'CLOAKING CONTROL',12,'#667887',True)
+        draw.rounded_rectangle((28,490,92,524),radius=17,fill='#137F73' if on else '#A8B4BC')
+        cx=75 if on else 45
+        draw.ellipse((cx-12,495,cx+12,519),fill='white')
+        text(108,491,'ON' if on else 'OFF',22,'#137F73' if on else '#667887',True)
+        text(195,493,('Fixed rotations' if self.scene=='steady' else 'Frozen PPO') if on else 'Zero rotation',18)
+        text(570,467,f't U/D = {step*U/D:.2f}',18,bold=True)
+        local=(step-1)%self.phase_steps+1
+        transition=(phase!='off_before' and local*U/D<=1.)
+        text(570,495,f'Playback {speed}x | '+('switch transient' if transition else phase.replace('_',' ')),16)
+        text(960,467,f'Fixed lab-frame inset | y/D={rows[-1,1]/D:+.3f}',13,bold=True)
+        # Motion trace: fixed +/-D limits, a fixed 60 D/U trailing time span.
+        box=(55,560,895,732)
+        draw.rectangle(box,outline='#CCD6DD')
+        middle=(box[1]+box[3])/2
+        draw.line((box[0],middle,box[2],middle),fill='#A8B4BC')
+        t=step*U/D
+        left=max(0.,t-60)
+        first=max(0,int(left*D/U))
+        sample=rows[first::max(1,(len(rows)-first)//1200)]
+        xy=[(box[0]+(r[0]*U/D-left)/60*(box[2]-box[0]),middle-r[1]/D*(box[3]-box[1])/2) for r in sample]
+        if len(xy)>1:draw.line(xy,fill='#263D4B',width=2)
+        text(28,535,'Transverse displacement y/D (fixed scale; dashed inset line = equilibrium)',14)
+        text(24,554,'+1',12);text(25,642,'0',12);text(24,722,'-1',12)
+        text(55,739,f'{left:.1f}',12);text(835,739,f'{left+60:.1f}',12)
+        text(350,739,'response time t U/D',12)
+        text(28,774,f'Vorticity D/U: fixed [-3, +3] | exploratory, not convergence proof | source {self.source[:8]}',12)
+        self.proc.stdin.write(canvas.tobytes())
+        local=step%self.phase_steps
+        if not self.records or local in (0,self.phase_steps//2) or self.records[-1]['phase']!=phase:
+            canvas.save(self.out/f'keyframe-{step:09d}.png')
+        elapsed=time.monotonic()-start
+        self.seconds+=elapsed
+        self.records.append(dict(frame=len(self.records),step=step,absolute_clock=sim.stepper.step_count,
+            tU_D=t,phase=phase,playback_speed=speed,y_D=float(rows[-1,1]/D),
+            surface_speeds=eff.tolist(),render_s=elapsed,
+            fraction_finite_outside_color_limits=float(np.mean(np.abs(omega[finite])>3))))
+
+    def close(self):
+        if self.closed:return
+        self.closed=True
+        self.proc.stdin.close()
+        code=self.proc.wait(timeout=60)
+        self.log.close()
+        receipt=dict(frames=self.records,frame_count=len(self.records),fps=30,size=SIZE,
+            duration_s=len(self.records)/30,render_total_s=self.seconds,ffmpeg_exit=code,
+            normalization='centered dx=1 vorticity D/U, all 5 stencil cells fluid',
+            color_limits=[-3,3],color_bins=21,main_limits_D=MAIN_LIMITS,inset_limits_D=INSET_LIMITS,
+            y_origin_lattice=CY,trace_limits_D=[-1,1],inset_camera='fixed lab frame',
+            acceleration='frame cadence only; native solver and SI unchanged')
+        (self.out/'video.json').write_text(json.dumps(receipt,indent=2)+'\n')
+        if code:raise RuntimeError(f'ffmpeg exit {code}')

research/continuous_viv/run.py

@@ -0,0 +1,222 @@
+"""Continuous OFF -> ON -> OFF VIV illustration; docs/continuous_viv_video.md.
+
+Each invocation owns ONE Simulation and ONE oscillator. Rendering never advances
+the solver or changes the policy clock. No topology changes or moving restart.
+"""
+from __future__ import annotations
+
+import argparse
+import json
+from pathlib import Path
+import subprocess
+import sys
+import time
+
+import numpy as np
+
+ROOT = Path(__file__).resolve().parents[2]
+sys.path.insert(0, str(ROOT))
+from research.karman_viv import run as base
+from research.continuous_viv.render import Video, frame_schedule
+
+D, U, SCALE, SI = 40., .005, 2., 3200
+PHASES = ('off_before', 'on', 'off_after')
+STEADY_WALL = np.array([0., -5.1, 5.1]) * U
+
+
+def timing(scene, smoke):
+    if smoke:
+        return 1,32,32
+    # Prior video's steady interval was 1000 at D20; diffusive D40 uses 4000.
+    return (320,4000,2048000) if scene=='steady' else (400,SI,2048000)
+
+
+def command(scene, phase, obs, model):
+    """Physical surface speeds; OFF is not a zero normalized PPO action."""
+    if phase != 'on':
+        return np.zeros(3), np.array([0., .5, -.5])
+    if scene == 'steady':
+        return STEADY_WALL.copy(), np.full(3, np.nan)
+    action = np.asarray(model.predict(obs, deterministic=True)[0], np.float32)
+    if action.shape != (3,) or not np.isfinite(action).all() or (np.abs(action)>1).any():
+        raise RuntimeError('invalid PPO action')
+    return base.wall_target(action, SCALE), action
+
+
+def add_scene(sim, scene):
+    if scene == 'karman':
+        return base.add_scene(sim, SCALE, moving=True)
+    from CelerisLab.body.motion import Translation
+    probes = [sim.add_body('sensor', c, 10.) for c in base.scaled_points(base.PROBES, SCALE)]
+    envelope = Translation((0., 0., -D, D), .025, shared=True)
+    bodies = [sim.add_body('circle', c, D/2, translation=envelope)
+              for c in base.scaled_points(base.CENTRES, SCALE)]
+    return probes, bodies
+
+
+def evolve(sim, probes, bodies, structure, force, *, scene, model, norm,
+           intervals, si, out, render, start, wall_cap):
+    """One native-step loop across all phases; completed SI feeds the next SI."""
+    from CelerisLab.body.motion import MovingTables
+    cols = {k:i for i,k in enumerate(MovingTables.DIAG)}
+    count = len(PHASES)*intervals*si
+    rows = np.lib.format.open_memmap(out/'steps.npy', mode='w+', dtype=np.float64,
+                                   shape=(count, len(base.STEP_COLUMNS)))
+    # The completed count is authoritative after any interrupted run.
+    completed = 0
+    obs = np.zeros(12, np.float32)
+    signals, switches = [], []
+    origin_clock = sim.stepper.step_count
+    status = 'running'
+    try:
+        for phase in PHASES:
+            switches.append(dict(phase=phase, step=completed, absolute_clock=sim.stepper.step_count,
+                                 y=structure.y, v=structure.v, previous_rhs=force,
+                                 input_observation=obs.tolist()))
+            for i in range(intervals):
+                target, action = command(scene, phase, obs, model)
+                # Match the frozen inference contract: reset the EMA every SI.
+                eff = np.zeros(3)
+                sensor_sum, force_sum = np.zeros(6), np.zeros(6)
+                for k in range(si):
+                    if k % 1000 == 0:
+                        base.check_time(start, wall_cap)
+                    y1, v1 = structure.candidate(force)
+                    if not np.isfinite([y1,v1]).all() or abs(y1)>D:
+                        raise RuntimeError('nonfinite structure or |y| > D; no retuning')
+                    eff = base.ema(eff, target, SCALE)
+                    omega = eff/(D/2)
+                    for bid,w in zip(bodies,omega):
+                        sim.set_body(bid, omega=float(w), offset=(0.,y1), vy=y1-structure.y)
+                    sim.run(1, zero_obs=True, sync_obs=True)
+                    d = sim.moving_diagnostics()
+                    cme,wen = d[:,cols['cme_y']], d[:,cols['wen_y']]
+                    node = d[:,cols['cover_jy']]-d[:,cols['uncover_jy']]
+                    row = (completed+1,y1,v1,*cme,*wen,*node,
+                           d[:,cols['eq_only']].sum(),d[:,cols['no_donor']].sum(),
+                           d[:,cols['crossings']].sum(),*omega)
+                    if not np.isfinite(row).all() or row[13] != 0:
+                        raise RuntimeError('nonfinite moving diagnostic or no_donor')
+                    f=np.array([v for p in bodies for v in sim.read_force(p)[:2]])
+                    s=np.array([v for p in probes for v in sim.read_sensor(p)[:2]])
+                    if not np.isfinite([f,s]).all():
+                        raise RuntimeError('nonfinite force/sensor observation')
+                    sensor_sum += s
+                    force_sum += f
+                    force=float(cme.sum(dtype=np.float64)+wen.sum(dtype=np.float64))
+                    rows[completed]=row
+                    completed += 1
+                    structure.y,structure.v=y1,v1
+                    if sim.stepper.step_count != origin_clock+completed:
+                        raise RuntimeError('noncontinuous solver clock')
+                    render(sim, completed, phase, eff, rows[:completed])
+                if scene == 'karman':
+                    raw, obs, clipped = base.observation(sensor_sum/si,force_sum/si,SCALE,norm)
+                else:
+                    # Steady uses prescribed rotations, not a model or borrowed norm.
+                    raw=np.r_[sensor_sum/si,force_sum/si]
+                    clipped=0
+                signals.append(dict(phase=phase, step=completed, raw=raw.tolist(),
+                                    observation=obs.tolist(), clipped=clipped,
+                                    action=action.tolist() if np.isfinite(action).all() else None,
+                                    target_surface_speed=target.tolist()))
+                if (i+1)%10==0 or i+1==intervals:
+                    rows.flush()
+                    progress=dict(phase=phase, interval=i+1, completed_steps=completed,
+                                  tU_D=completed*U/D, y_D=structure.y/D,
+                                  elapsed_s=time.monotonic()-start)
+                    (out/'progress.json').write_text(json.dumps(progress)+'\n')
+                    print(json.dumps(progress),flush=True)
+        status='duration_completed_not_physics_qualified'
+    except BaseException as exc:
+        status=f'error: {type(exc).__name__}: {exc}'
+        raise
+    finally:
+        rows.flush()
+        (out/'signals.json').write_text(json.dumps(signals,allow_nan=False)+'\n')
+        report=dict(status=status,completed_steps=completed,origin_clock=origin_clock,
+                    final_clock=sim.stepper.step_count,switches=switches,
+                    elapsed_s=time.monotonic()-start,columns=base.STEP_COLUMNS,
+                    M=structure.mass,K=structure.stiffness,fn=structure.frequency,damping=0.,
+                    max_abs_y_D=float(np.abs(rows[:completed,1]).max()/D) if completed else None)
+        (out/'report.json').write_text(json.dumps(report,indent=2,allow_nan=False)+'\n')
+
+
+def main():
+    ap=argparse.ArgumentParser(description=__doc__)
+    ap.add_argument('--scene',choices=('steady','karman'),required=True)
+    ap.add_argument('--out',type=Path,required=True)
+    ap.add_argument('--model',type=Path)
+    ap.add_argument('--norm',type=Path)
+    ap.add_argument('--release-dir',type=Path)
+    ap.add_argument('--smoke',action='store_true',help='3 x 32 steps, all switches, no physical result')
+    ap.add_argument('--wall-seconds',type=float,default=10800.)
+    args=ap.parse_args()
+    if not np.isfinite(args.wall_seconds) or args.wall_seconds<=0:
+        ap.error('finite positive wall cap required')
+    if args.scene=='karman':
+        if args.model is None or args.norm is None:
+            ap.error('Karman requires the frozen model and norm')
+        if base.sha(args.model)!=base.MODEL_SHA or base.sha(args.norm)!=base.NORM_SHA:
+            ap.error('frozen Karman model/normalizer hash mismatch')
+        if not args.smoke and args.release_dir is None:
+            ap.error('physical Karman requires its verified static release')
+    elif any(p is not None for p in (args.model,args.norm,args.release_dir)):
+        ap.error('steady uses constant control and its own warmup, no borrowed policy or release')
+    out=args.out.resolve()
+    if out==ROOT or ROOT in out.parents:
+        ap.error('output must be outside checkout')
+    out.mkdir(parents=True,exist_ok=False)
+    intervals,si,warm=timing(args.scene,args.smoke)
+    cfgpath=out/'config.json'
+    cfgpath.write_text(json.dumps(base.config(40),indent=2)+'\n')
+    source=subprocess.check_output(['git','-C',str(ROOT),'rev-parse','HEAD'],text=True).strip()
+    case=dict(source=source,scene=args.scene,smoke=args.smoke,D=D,U=U,nu=.004,
+              phases=PHASES,SI=si,intervals_per_phase=intervals,
+              phase_duration_tU_D=intervals*si*U/D, mass_ratio=4.,fn_D_U=.2,zeta=0.,
+              first_ppo_observation='previous complete OFF interval; no reset at switches',
+              rhs='previous native step CME_y+Wen_y; no filtering; policy uses SI-mean CME',
+              control='fixed (0,-5.1,+5.1)U' if args.scene=='steady' else 'frozen karman_re100 PPO; actual Re_D50',
+              model_sha256=base.sha(args.model) if args.model else None,
+              norm_sha256=base.sha(args.norm) if args.norm else None,
+              config_sha256=base.sha(cfgpath),render='sparse post-step snapshots only; no step skipping')
+    model=base.load_policy(args.model) if args.model else None
+    norm=json.loads(args.norm.read_text()) if args.norm else None
+    seed=None
+    if args.scene=='karman' and not args.smoke:
+        seed,force=base.saved_release(args.release_dir)
+        case['release_sha256']=base.sha(args.release_dir/'release.npz')
+    (out/'case.json').write_text(json.dumps(case,indent=2)+'\n')
+    from CelerisLab import Simulation
+    from CelerisLab.body.motion import MovingTables
+    started=time.monotonic()
+    structure=base.Structure(4.,2.5e-5,0.,D)
+    video=None
+    try:
+        with Simulation(str(cfgpath)) as sim:
+            probes,bodies=add_scene(sim,args.scene)
+            sim.initialize()
+            if seed is not None:
+                case['release_transfer']=base.install_seed(sim,seed)
+                del seed
+            else:
+                base.warmup(sim,warm-1,started,args.wall_seconds)
+                sim.run(1)
+                d=sim.moving_diagnostics()
+                force=float(d[:,MovingTables.DIAG.index('cme_y')].sum(dtype=np.float64)
+                            +d[:,MovingTables.DIAG.index('wen_y')].sum(dtype=np.float64))
+            case.update(warmup_executed=0 if 'release_transfer' in case else warm,
+                        origin_clock=sim.stepper.step_count)
+            (out/'case.json').write_text(json.dumps(case,indent=2)+'\n')
+            phase_steps=intervals*si
+            video=Video(out,args.scene,bodies,phase_steps,smoke=args.smoke,source=source)
+            evolve(sim,probes,bodies,structure,force,scene=args.scene,model=model,norm=norm,
+                   intervals=intervals,si=si,out=out,render=video.capture,
+                   start=started,wall_cap=args.wall_seconds)
+    finally:
+        if video is not None:
+            video.close()
+
+
+if __name__=='__main__':
+    main()

research/duetsch_re100_kc5/protocol.py

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

research/fig2_collection/analyze.py

@@ -0,0 +1,130 @@
+"""Offline Fig2 estimands; preserve raw same-clock errors beside registration.
+
+No thresholds are fitted here. Probe DTW is the actual Legacy helper, not a
+field metric; independent wake phase is a diagnostic, not a time correction.
+"""
+import argparse
+import gc
+import importlib.util
+import json
+from pathlib import Path
+import sys
+import types
+
+import numpy as np
+
+ROOT=Path(__file__).resolve().parents[2]
+sys.path.insert(0,str(ROOT))
+from research.fig2_collection import run
+from research.fig2_materials.metrics import field_error, motion_summary
+
+
+def frequency(signal,dt):
+    """Largest non-DC rectangular-window FFT bin, descriptive resolution only."""
+    a=np.asarray(signal,float);a=a-a.mean()
+    if len(a)<3 or np.max(np.abs(a))<1e-12:return None
+    power=np.abs(np.fft.rfft(a))**2
+    i=1+int(np.argmax(power[1:]))
+    return dict(f_D_U=float(np.fft.rfftfreq(len(a),dt)[i]),bin_width_D_U=1/(len(a)*dt),
+                method='demeaned rectangular window; dominant non-DC bin; not a confidence interval')
+
+
+def native_summary(rows,scene):
+    start=max(0,len(rows)-640000) if scene!='taylor' else 0
+    a=rows[start:]
+    y=a[:,1]/40.;cy=a[:,3:9].sum(axis=1)/(.5*.005**2*40)
+    return dict(window_D_U=[start*.005/40,len(rows)*.005/40],
+                motion=motion_summary(y),motion_frequency=frequency(y,.005/40),
+                structural_RHS_Cy_RMS=float(np.sqrt(np.mean(cy**2))),
+                structural_RHS_Cy_mean=float(cy.mean()),
+                max_native_delta_Cy=float(np.abs(np.diff(cy)).max()) if len(cy)>1 else None,
+                force_frequency=frequency(cy,.005/40),
+                no_donor=float(rows[:,13].sum()),eq_only=float(rows[:,12].sum()),
+                crossings=float(rows[:,14].sum()),
+                force_convention='body force, sum pinball CME_y+Wen_y; native previous-step RHS, no filtering')
+
+
+def compare_selected(target,current,u0=.005):
+    """Fixed lab ROI; all actual compared roles must be fluid and finite here."""
+    ys,xs=slice(311,711),slice(1360,2161)
+    t,c=target['ux'][:,ys,xs],current['ux'][:,ys,xs]
+    tv,cv=target['uy'][:,ys,xs],current['uy'][:,ys,xs]
+    if t.shape!=c.shape or tv.shape!=cv.shape:raise ValueError('selected shape mismatch')
+    for fields in (target,current):
+        if not np.all((fields['flags'][:,ys,xs]&1)!=0):raise ValueError('nonfluid selected ROI')
+    if not all(np.isfinite(a).all() for a in (t,c,tv,cv)):raise ValueError('nonfinite selected ROI')
+    e2=((c.astype(float)-t)/u0)**2+((cv.astype(float)-tv)/u0)**2
+    return dict(per_slot_E=np.sqrt(e2.mean(axis=(1,2))),aggregate_E=float(np.sqrt(e2.mean())))
+
+
+def main():
+    ap=argparse.ArgumentParser(description=__doc__);ap.add_argument('directory',type=Path)
+    args=ap.parse_args();p=args.directory
+    case=json.loads((p/'case.json').read_text());term=json.loads((p/'terminal.json').read_text())
+    if case['plan']['smoke'] or term['status']!='three_roles_completed_not_qualified':
+        raise ValueError('physical analysis requires all three complete roles; smoke/partial is unscored')
+    helpers=run.load_helpers()
+    # Reuse the canonical pure module under the same isolated read-only package.
+    native=run.load_file('_fig2_readonly_drl.legacy_test.metrics',
+                         run.ACQUISITION.parent/'legacy_test/metrics.py')
+    reports={k:json.loads((p/k/'report.json').read_text()) for k in ('target','zero','controlled')}
+    sig={}
+    result=dict(scene=case['scene'],source=case['source'],Re_D=50.,roles={},
+                primary='unaligned same-clock/event U-normalized vector errors; common lab ROI',
+                periodic_phase='independently registered centre-uy phase diagnostic, not a shared-clock comparison',
+                notes=['Modern moving-plant data; no Legacy numerical equivalence or physical qualification',
+                       'Motion, native structural force, policy SI means and probe DTW are different estimands'])
+    for role in ('target','zero','controlled'):
+        with np.load(p/role/'signals.npz') as z:sig[role]={k:z[k] for k in z.files}
+        if role!='target':
+            with np.load(p/role/'native.npz') as z:
+                result['roles'][role]=native_summary(z['rows'],case['scene'])
+    t=sig['target']['sensors']
+    for role in ('zero','controlled'):
+        s=sig[role]['sensors']
+        if not np.array_equal(sig[role]['absolute_steps'],sig['target']['absolute_steps']):
+            raise ValueError('probe absolute clocks mismatch')
+        result['roles'][role]['unaligned_probe_area_velocity_RMS_over_U']=float(np.sqrt(np.mean(((s/156-t/156)/.005)**2)))
+        result['roles'][role]['probe_channel_correlations']=native.channel_correlations(t,s)
+        if case['scene']=='taylor':
+            spec=types.SimpleNamespace(scene='vortex',conv_len=30)
+            dtw=native.native_similarity(spec,t,s,len(s)-1)
+            result['roles'][role]['native_probe_DTW']=dict(similarity=dtw,definition='source vortex training-window roll; supplementary, not event field recovery',conv_len=30)
+            own=helpers.select_vortex_event(s,'vortex_taylor')
+            target_event=helpers.select_vortex_event(t,'vortex_taylor')
+            result['roles'][role]['own_event_peak']=int(own['peak_index'])
+            result['roles'][role]['event_peak_delay_D_U']=(int(own['peak_index'])-int(target_event['peak_index']))*.4
+        else:
+            result['roles'][role]['native_probe_DTW']=native.cycle_similarity(t,s,30)
+            result['roles'][role]['native_probe_DTW'].update(definition='canonical six-channel legacy cycle DTW, training D20 units',
+                lag_column=1,target_template_indices=[30,60],state_window_indices=[len(s)-30,len(s)],
+                registration_only=True,lag_identifiable=bool(case['scene']=='karman' and np.std(t[:,1])>1e-10))
+    data={}
+    for role in ('target','zero','controlled'):
+        with np.load(p/role/'fields.npz') as z:data[role]={k:z[k] for k in z.files}
+    for role in ('zero','controlled'):
+        name={'karman':'independent_phase8_diagnostic','steady':'late_snapshot','taylor':'common_event8'}[case['scene']]
+        if case['scene']!='karman' or all(len(data[k]['ux'])==8 for k in (role,'target')):
+            result['roles'][role][name]=compare_selected(data['target'],data[role])
+        else:
+            result['roles'][role][name]=dict(status='unavailable_insufficient_periodic_crossings')
+        result['roles'][role]['selection_metadata']=reports[role]['selection']
+        result['roles'][role]['field_window_tU_D']=reports[role]['field_window_tU_D']
+        if case['scene']=='karman' and all('cycle_mean_ux' in data[k] for k in (role,'target')):
+            support=np.zeros(data[role]['cycle_mean_ux'].shape,bool);support[311:711,1360:2161]=True
+            result['roles'][role]['E_complete_cycle_mean']=field_error(
+                data[role]['cycle_mean_ux'],data[role]['cycle_mean_uy'],
+                data['target']['cycle_mean_ux'],data['target']['cycle_mean_uy'],.005,support)
+    if case['scene']=='taylor':
+        clocks=[data[k]['absolute_steps'] for k in data]
+        if not all(np.array_equal(clocks[0],v) for v in clocks[1:]):raise ValueError('event snapshot clocks differ')
+        result['Taylor_injection']=dict(operator='all interior equilibrium rewrite; both buffers, identical event clock',
+            strength=.000075,center=[600.,510.5],radius=80.,impulse_window='first two SI; retained native history',
+            source_offsets=list(helpers.VORTEX_EVENT_OFFSETS))
+    result['same_clock_field_errors']=json.loads((p/'same_clock_errors.json').read_text())
+    run.write_json(p/'statistics.json',result)
+    del data;gc.collect()
+    print(json.dumps(dict(scene=case['scene'],statistics=str(p/'statistics.json'))))
+
+
+if __name__=='__main__':main()

research/fig2_collection/plot.py

@@ -0,0 +1,93 @@
+"""Shared scientific plots of retained fields; no solver or synthetic fields."""
+import argparse
+import json
+from pathlib import Path
+import sys
+
+import numpy as np
+import matplotlib
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+from matplotlib.colors import BoundaryNorm
+from matplotlib.patches import Circle,Rectangle
+
+ROOT=Path(__file__).resolve().parents[2]
+sys.path.insert(0,str(ROOT))
+from research.continuous_viv.render import vorticity
+
+
+def decorate(ax,scene,role,y,*,mean=False):
+    if not mean:
+        if scene=='karman':ax.add_patch(Circle((-20,0),1,facecolor='white',edgecolor='#253544',lw=.7))
+        if role!='target':
+            for x,dy in ((0,0),(1.3,.75),(1.3,-.75)):
+                ax.add_patch(Circle((x,dy+y),.5,facecolor='white',edgecolor='#253544',lw=.7))
+    ax.add_patch(Rectangle((4,-5),20,10,fill=False,ls='--',lw=.6,edgecolor='#555'))
+    ax.set_xlim(-23.6,27.6);ax.set_ylim(-6.4,6.4)
+    ax.set_xlabel('(x - pinball-front equilibrium centre) / D',fontsize=8)
+    ax.set_ylabel('(y - channel centre) / D',fontsize=8)
+    ax.tick_params(labelsize=8)
+
+
+def main():
+    ap=argparse.ArgumentParser(description=__doc__);ap.add_argument('directory',type=Path)
+    args=ap.parse_args();p=args.directory
+    case=json.loads((p/'case.json').read_text());scene=case['scene']
+    cmap=plt.get_cmap('RdBu_r',21).copy();cmap.set_bad('#DDE2E5')
+    norm=BoundaryNorm(np.linspace(-3,3,22),21,clip=True)
+    prefix='SMOKE — plumbing only | ' if case['plan']['smoke'] else ''
+    roles=('target','zero','controlled')
+    data={};motion={}
+    for role in roles:
+        with np.load(p/role/'fields.npz') as z:data[role]={k:z[k] for k in z.files}
+        with np.load(p/role/'native.npz') as z:
+            motion[role]=z['rows'][:,1].copy()/40.
+    shape=data['target']['ux'].shape[-2:]
+    extent=((-0.5-1200)/40,(shape[1]-.5-1200)/40,
+            (-.5-510.5)/40,(shape[0]-.5-510.5)/40)
+    fig,axs=plt.subplots(3,1,figsize=(13,10),constrained_layout=True)
+    for ax,role in zip(axs,roles):
+        d=data[role];flags=d['common_fluid_all_boundaries'].astype(np.uint16)
+        w=vorticity(d['fixed_window_mean_ux'],d['fixed_window_mean_uy'],flags)
+        image=ax.imshow(w,origin='lower',extent=extent,cmap=cmap,norm=norm,interpolation='nearest')
+        decorate(ax,scene,role,0,mean=True)
+        label='no pinball reference' if role=='target' else ('zero rotation' if role=='zero' else ('constant rotations' if scene=='steady' else 'frozen PPO'))
+        ax.set_title(f'{role}: {label} | arithmetic fixed-window mean',fontsize=10)
+    fig.suptitle(prefix+f'{scene}: modern Re_D=50, m*=4, zero structural damping\n'
+                 'Gray: cells outside the common fluid support; mean fields are not instantaneous geometry',fontsize=12)
+    fig.colorbar(image,ax=list(axs),shrink=.65,label='vorticity D/U, fixed limits [-3,3]')
+    fig.savefig(p/'mean-flow.png',dpi=160);plt.close(fig)
+    slots=max(len(data[r]['ux']) for r in roles)
+    fig,axs=plt.subplots(3,slots,figsize=(max(8,3*slots),7),squeeze=False,constrained_layout=True)
+    for row,role in enumerate(roles):
+        d=data[role]
+        for col,ax in enumerate(axs[row]):
+            if col>=len(d['ux']):ax.axis('off');continue
+            w=vorticity(d['ux'][col],d['uy'][col],d['flags'][col])
+            image=ax.imshow(w,origin='lower',extent=extent,cmap=cmap,norm=norm,interpolation='nearest')
+            step=int(d['selected_response_steps'][col]);y=float(motion[role][step-1])
+            decorate(ax,scene,role,y)
+            ax.set_xlim(-3,14);ax.set_ylim(-3,3)
+            ax.set_xlabel('x/D relative to equilibrium front',fontsize=7)
+            ax.set_title(f'{role} | slot {col+1}\ntU/D={step*.005/40:.2f}',fontsize=8)
+    family={'karman':'independently phased diagnostic snapshots (different real times)',
+            'steady':'late snapshots at matching real times',
+            'taylor':'common event-relative snapshots at matching real times'}[scene]
+    fig.suptitle(prefix+f'{scene}: {family}\nFixed lab coordinates; same vorticity normalization/limits; source '+case['source'][:8],fontsize=12)
+    fig.colorbar(image,ax=list(axs.ravel()),shrink=.7,label='vorticity D/U [-3,3]')
+    fig.savefig(p/'selected-flow.png',dpi=160);plt.close(fig)
+    fig,ax=plt.subplots(figsize=(11,3.5),constrained_layout=True)
+    for role,color in (('zero','#5A6570'),('controlled','#147D73')):
+        y=motion[role];stride=max(1,len(y)//12000)
+        time=(np.arange(len(y))+1)*.005/40
+        ax.plot(time[::stride],y[::stride],lw=.8,label=role,color=color)
+    ax.set(xlabel='response time tU/D',ylabel='transverse y/D',ylim=(-1,1))
+    ax.legend();ax.grid(alpha=.2)
+    ax.set_title(prefix+scene+': independent paired releases; full retained duration')
+    fig.savefig(p/'motion.png',dpi=160);plt.close(fig)
+    print(json.dumps(dict(scene=scene,plots=['mean-flow.png','selected-flow.png','motion.png'],
+                          vorticity_limits=[-3,3],normalization='curl(u) D/U; five-node fluid stencil',
+                          mean_mask='intersection over collected boundary fields',display_interpolation='nearest')))
+
+
+if __name__=='__main__':main()

research/fig2_collection/run.py

@@ -0,0 +1,557 @@
+"""Research-only Fig. 2 acquisition on the modern moving plant.
+
+These paired releases are new experiments, not Legacy numerical reproductions.
+Nothing here qualifies a steady state or automatically extends/retunes a run.
+GPU dependencies are lazy so the collection contracts can be checked on CPU.
+"""
+from __future__ import annotations
+
+import argparse
+import ast
+import gc
+import hashlib
+import importlib.util
+import json
+import os
+from pathlib import Path
+import platform
+import subprocess
+import sys
+import time
+import types
+
+import numpy as np
+
+sys.dont_write_bytecode = True
+ROOT = Path(__file__).resolve().parents[2]
+D, U, NU, SI = 40., .005, .004, 3200
+BASELINE = '76873d911c4e2caba1bacc25aaa1495ea768453f'
+CENTRES = ((1200., 510.5), (1252., 540.5), (1252., 480.5))
+PROBES = ((1600., 590.5), (1600., 510.5), (1600., 430.5))
+BIAS = np.array([0., -4., 4.])
+STEADY_SURFACE = np.array([0., -5.1, 5.1])*U
+MIN_AVAILABLE = 16*1024**3
+ACQUISITION = Path('/home/frank14f/DynamisLab/src/drl_pinball/acquisition.py')
+K_RELEASE = Path('/tmp/celeris-results-karman-viv-516975e/pair')
+ASSETS = {
+    'karman': (Path('/tmp/celeris-karman-viv-assets/model.zip'),
+               '148240c1dcb8b11d8e5a0a9e991f6874c4e799c9744707bf0fb32a8efeb0468d',
+               Path('/tmp/celeris-karman-viv-assets/norm.json'),
+               '534fda818c030a56960e95d0f3ac1e350d887d03db7f80f1d6453195a502b7f0'),
+    'taylor': (Path('/home/frank14f/DynamisLab/models/old/vortex_taylor.zip'),
+               '58882721900f1eb0779f64aea0a2e233d307a94b640b9aaf33efa036db3c79e3',
+               Path('/home/frank14f/DynamisLab/src/SR_analysis/data/vortex/vortex_taylor/norm.json'),
+               'bd38b4cf2b6da56b5cf378e646d673b97a606531c68163460313b330c9957b86'),
+}
+NATIVE_COLUMNS = ('response_step', 'y', 'v', 'cme_F', 'cme_U', 'cme_L',
+                  'wen_F', 'wen_U', 'wen_L', 'node_F', 'node_U', 'node_L',
+                  'eq_only', 'no_donor', 'crossings', 'omega_F', 'omega_U', 'omega_L')
+
+
+def sha(path):
+    digest = hashlib.sha256()
+    with Path(path).open('rb') as stream:
+        for block in iter(lambda: stream.read(1024*1024), b''):
+            digest.update(block)
+    return digest.hexdigest()
+
+
+def json_value(value):
+    if isinstance(value, dict):
+        return {str(k): json_value(v) for k, v in value.items()}
+    if isinstance(value, (tuple, list)):
+        return [json_value(v) for v in value]
+    if isinstance(value, np.ndarray):
+        return json_value(value.tolist())
+    if isinstance(value, np.generic):
+        return value.item()
+    if isinstance(value, Path):
+        return str(value.resolve())
+    return value
+
+
+def write_json(path, value):
+    Path(path).write_text(json.dumps(json_value(value), indent=2, allow_nan=False)+'\n')
+
+
+def load_file(name, path):
+    spec = importlib.util.spec_from_file_location(name, path)
+    module = importlib.util.module_from_spec(spec)
+    sys.modules[name] = module
+    spec.loader.exec_module(module)
+    return module
+
+
+def load_helpers(path=ACQUISITION):
+    """Load read-only CPU helpers under an isolated package; no original imports.
+
+    Empty namespace packages avoid executing application __init__ modules. The
+    acquisition module and its pure DTW dependency execute directly from source.
+    Event functions are AST-selected verbatim to avoid acquire.py's heavy imports;
+    offsets come from that acquisition module, never a copied/fabricated selector.
+    """
+    path = Path(path)
+    package = '_fig2_readonly_drl'
+    for suffix, directory in (('', path.parent), ('.legacy_test', path.parent/'legacy_test'),
+                              ('.legacy_test.core', path.parent/'legacy_test/core')):
+        module = types.ModuleType(package+suffix)
+        module.__path__ = [str(directory)]
+        sys.modules[package+suffix] = module
+    acquisition = load_file(package+'.acquisition', path)
+    event_path = path.parent/'legacy_test/acquire.py'
+    source = ast.parse(event_path.read_text(), filename=str(event_path))
+    names = {'_vortex_event_channel', '_select_vortex_event'}
+    nodes = [node for node in source.body if isinstance(node, ast.FunctionDef) and node.name in names]
+    if {node.name for node in nodes} != names:
+        raise ValueError('authoritative event selector functions missing')
+    namespace = dict(np=np, VORTEX_EVENT_OFFSETS=acquisition.VORTEX_EVENT_OFFSETS)
+    exec(compile(ast.Module(body=nodes, type_ignores=[]), str(event_path), 'exec'), namespace)
+    acquisition.select_vortex_event = namespace['_select_vortex_event']
+    return acquisition
+
+
+def schedule(scene, smoke=False):
+    if scene not in ('karman', 'steady', 'taylor'):
+        raise ValueError('unknown scene')
+    intervals = 150 if scene == 'taylor' else 400
+    release = {'karman': 4576000, 'steady': 2048000, 'taylor': 1024000}[scene]
+    return dict(si=32 if smoke else SI, intervals=3 if smoke else intervals,
+                field_start=0 if smoke or scene == 'taylor' else 240,
+                release_clock=(64 if scene == 'taylor' else 32) if smoke else release,
+                smoke=bool(smoke))
+
+
+def available_memory(path=Path('/proc/meminfo')):
+    for line in Path(path).read_text().splitlines():
+        if line.startswith('MemAvailable:'):
+            return int(line.split()[1])*1024
+    raise RuntimeError('MemAvailable unavailable; RAM admission cannot be verified')
+
+
+def admit_memory(smoke=False):
+    if smoke:
+        return None
+    available = available_memory()
+    if available < MIN_AVAILABLE:
+        raise MemoryError(f'MemAvailable {available} < required {MIN_AVAILABLE}')
+    return available
+
+
+def wall_surface(scene, role, action=None):
+    if role == 'target' or role == 'zero':
+        return np.zeros(3, np.float64)  # exactly zero physical rotation
+    if scene == 'steady':
+        return STEADY_SURFACE.copy()
+    action = np.asarray(action)
+    if action.shape != (3,) or not np.isfinite(action).all() or (np.abs(action) > 1).any():
+        raise ValueError('invalid deterministic policy action')
+    return ((8. if scene == 'karman' else 4.)*action+BIAS)*U
+
+
+def check_structure(y, v, previous_y):
+    if not np.isfinite([y, v]).all() or abs(y) > D or abs(v) > .025 or abs(y-previous_y) > .025:
+        raise RuntimeError('nonfinite structure, travel > D or speed > .025; no retuning')
+
+
+def check_clock(stepper, release_clock, completed):
+    if stepper.step_count != release_clock+completed:
+        raise RuntimeError('observed solver clock differs from recorded response clock')
+
+
+def clean_field(fields, flags):
+    """Canonical modern bit-0 fluid mask; hard-fail before nonfluid zero fill."""
+    ux, uy = np.asarray(fields['ux']), np.asarray(fields['uy'])
+    rho = np.asarray(fields['rho'])
+    if ux.ndim != 2 or uy.shape != ux.shape or rho.shape != ux.shape:
+        raise ValueError('canonical shared (ny,nx) macroscopic arrays required')
+    flags = np.asarray(flags)
+    if flags.size != ux.size or flags.dtype != np.uint16:
+        raise ValueError('modern uint16 flags matching field required')
+    flags = flags.reshape(ux.shape)
+    fluid = (flags & 1) != 0
+    if any(not np.isfinite(value[fluid]).all() for value in (rho, ux, uy)):
+        raise RuntimeError('nonfinite FLUID macroscopic field')
+    result = np.zeros((2, *ux.shape), np.float32)
+    result[0, fluid], result[1, fluid] = ux[fluid], uy[fluid]
+    return result, flags.copy(), fluid
+
+
+def select_fields(scene, plan, helpers, times, sensors, cache, event=None):
+    """Fixed-window means remain separate from independently phased diagnostics."""
+    mean = helpers.accumulate_mean_fields(cache[:, 0], cache[:, 1], np.arange(len(cache)))
+    result = dict(fixed_window_mean_ux=mean['mean_ux'], fixed_window_mean_uy=mean['mean_uy'])
+    metadata = dict(fixed_window_count=len(cache), periodic_phase_status='not_requested')
+    indices = np.array([len(cache)-1], np.int64)
+    if plan['smoke']:
+        metadata['selection'] = 'smoke_last_snapshot_not_physical'
+    elif scene == 'taylor':
+        if event is None:
+            event = helpers.select_vortex_event(sensors, 'vortex_taylor')
+        indices = np.asarray(event['field_indices'], np.int64)
+        if (indices < 0).any() or (indices >= len(cache)).any():
+            raise ValueError('authoritative event offsets exceed collected fields')
+        metadata.update(selection='source_vortex_event', event=event)
+    elif scene == 'karman':
+        try:
+            phase = helpers.assign_periodic_phase(times, sensors, minimum_crossing_interval=.5*30*SI)
+            complete = helpers.complete_cycle_field_indices(times, phase['crossing_times'])
+            chosen = helpers.select_phase_fields(times, phase['crossing_times'])
+            cycle_mean = helpers.accumulate_mean_fields(cache[:, 0], cache[:, 1], complete['field_indices'])
+            result.update(cycle_mean_ux=cycle_mean['mean_ux'], cycle_mean_uy=cycle_mean['mean_uy'])
+            indices = chosen['field_indices']
+            metadata.update(periodic_phase_status='available_diagnostic_only', phase_selection=chosen,
+                            complete_cycles=complete, crossing_times=phase['crossing_times'],
+                            rejected_crossing_times=phase['rejected_crossing_times'])
+        except ValueError as exc:
+            metadata.update(periodic_phase_status='insufficient_crossings', phase_reason=str(exc))
+        metadata['selection'] = 'independent_center_probe_uy_diagnostic'
+    else:
+        metadata.update(selection='last_snapshot', periodic_phase_status='not_assumed_steady')
+    result.update(ux=cache[indices, 0], uy=cache[indices, 1], selected_index=indices,
+                  selected_response_steps=times[indices])
+    return result, metadata, event
+
+
+def add_scene(sim, scene, *, moving=False, pinball=True):
+    from CelerisLab.body.motion import Translation
+    if scene == 'karman':
+        sim.add_body('circle', (400., 510.5), 40.)
+    probes = [sim.add_body('sensor', point, 10.) for point in PROBES]
+    envelope = Translation((0., 0., -D, D), .025, shared=True) if moving else None
+    bodies = [sim.add_body('circle', point, 20., translation=envelope) for point in CENTRES] if pinball else []
+    return probes, bodies
+
+
+def inject_taylor(sim):
+    from CelerisLab.lbm.initializers import add_vortex
+    clock = int(sim.stepper.step_count)
+    with sim.ctx.activate():
+        sim.stream.synchronize()
+        add_vortex(sim.field, (600., 510.5), 80., .000075, 'taylor')
+        sim.stream.synchronize()
+    if sim.stepper.step_count != clock:
+        raise RuntimeError('vortex initializer advanced response clock')
+    return dict(absolute_clock=clock, response_clock=0, center=[600., 510.5], radius=80.,
+                strength=.000075, type='taylor', count=1,
+                operator='unmodified add_vortex: all interior reset to equilibrium; both buffers uploaded',
+                impulse='full equilibrium reset discards pre-existing non-equilibrium on all interior cells')
+
+
+def make_release(scene, plan, cfgpath, base, out, started, cap, release_dir):
+    from CelerisLab import Simulation
+    if scene == 'karman' and not plan['smoke']:
+        seed, force = base.saved_release(release_dir)
+        return seed, force, dict(kind='verified_saved_static_release', path=release_dir,
+                                 release_npz_sha256=sha(release_dir/'release.npz'))
+    with Simulation(str(cfgpath)) as sim:
+        _, bodies = add_scene(sim, scene, pinball=scene != 'taylor')
+        sim.initialize()
+        warm = 32 if plan['smoke'] else (512000 if scene == 'taylor' else 2048000)
+        if scene == 'taylor':
+            base.warmup(sim, warm, started, cap)
+            first = sim.bodies.count
+            for point in CENTRES:
+                sim.add_body('circle', point, 20.)  # runtime additions return staged ID -1
+            sim.sync_bodies()
+            bodies = list(range(first, first+len(CENTRES)))
+        base.warmup(sim, warm-1, started, cap)
+        sim.run(1, zero_obs=True, sync_obs=True)
+        force = sum(float(sim.read_force(body)[1]) for body in bodies)
+        seed = base.capture_seed(sim)
+        if seed[3] != plan['release_clock'] or not np.isfinite(force):
+            raise RuntimeError('unexpected static release clock/force')
+        np.savez_compressed(out/'release.npz', current=seed[0], spare=seed[1], flags=seed[2],
+                            clock=seed[3], force=force)
+        return seed, force, dict(kind='fresh_static_full_two_buffer_release',
+                                 release_npz_sha256=sha(out/'release.npz'))
+
+
+def compress_native(path, rows, completed):
+    """Bounded disk memmap; preserve raw file until parent elects removal."""
+    rows.flush()
+    compressed = path.with_suffix('.npz')
+    np.savez_compressed(compressed, rows=rows[:completed])
+    with np.load(compressed) as data:
+        retained = data['rows']
+        if retained.shape != (completed, len(NATIVE_COLUMNS)):
+            raise RuntimeError('compressed native shape mismatch')
+        for start in range(0, completed, 10000):
+            end = min(start+10000, completed)
+            if not np.array_equal(retained[start:end], rows[start:end]):
+                raise RuntimeError('compressed native verification failed')
+    return dict(raw_path=path.name, compressed_path=compressed.name, verified=True,
+                compressed_sha256=sha(compressed), raw_preserved_for_parent_cleanup=True)
+
+
+def run_role(scene, role, plan, cfgpath, base, helpers, metrics, seed, initial_force,
+             model, norm, out, started, cap, event):
+    from CelerisLab import Simulation
+    from CelerisLab.body.motion import MovingTables
+    directory = out/role
+    directory.mkdir()
+    steps = plan['si']*plan['intervals']
+    rows = np.lib.format.open_memmap(directory/'native.npy', mode='w+', dtype=np.float64,
+                                    shape=(steps, len(NATIVE_COLUMNS)))
+    signals = {key: np.zeros((plan['intervals'], width), np.float32) for key, width in
+               [('sensor_area_mean', 6), ('sensors', 6), ('forces', 6), ('actions', 3),
+                ('observations', 12), ('policy_input', 12), ('raw_observation', 12)]}
+    structure = base.Structure(4., 2.5e-5, 0., D)
+    columns = {name: i for i, name in enumerate(MovingTables.DIAG)}
+    completed = completed_si = 0
+    duration = steps*U/D
+    report = dict(status='running', role=role, target_is_free_body=False if role == 'target' else None,
+                  columns=NATIVE_COLUMNS,
+                  motion_window=dict(start_tU_D=max(0., duration-80.), end_tU_D=duration,
+                                     interpretation='smoke_short_not_physical' if plan['smoke'] else
+                                     ('entire_60DU_event_not_settled' if scene == 'taylor' else 'declared_last_80DU')),
+                  structure=None if role == 'target' else dict(mass=structure.mass, fn=structure.frequency,
+                                                             mass_ratio=4., zeta=0., y0=0., v0=0.))
+    cache, common_fluid, selected_flags, roi_cache = None, None, {}, None
+    obs, force = np.zeros(12, np.float32), initial_force
+    try:
+        report['MemAvailable_admission_bytes'] = admit_memory(plan['smoke'])
+        with Simulation(str(cfgpath)) as sim:
+            probes, bodies = add_scene(sim, scene, moving=role != 'target', pinball=role != 'target')
+            sim.initialize()
+            if role == 'target':
+                base.warmup(sim, plan['release_clock'], started, cap)
+                report['initialization'] = 'own no-pinball initialization at matching absolute clock'
+            else:
+                report['release_transfer'] = base.install_seed(sim, seed)
+            if sim.stepper.step_count != plan['release_clock']:
+                raise RuntimeError('role release clocks differ')
+            if scene == 'taylor':
+                report['injection'] = inject_taylor(sim)
+            count = plan['intervals']-plan['field_start']
+            shape = (sim.field.ny, sim.field.nx)
+            cache = np.empty((count, 2, *shape), np.float32)
+            common_fluid = np.ones(shape, bool)
+            roi = metrics.roi_slices(shape, D=D, xs=1600., y0=510.5,
+                                     downstream_solid_extent=1272.)
+            # Only compact ROI histories survive a role; full-grid histories never reach disk.
+            roi_shape = cache[(slice(None), slice(None), *roi)].shape
+            roi_cache = np.empty(roi_shape, np.float32)
+            retain = set(event['field_indices'].tolist()) if event is not None else set()
+            if scene == 'taylor' and event is None and not plan['smoke']:
+                # Authoritative selector searches 38..87; any imported offset can be selected.
+                retain = {i+int(offset) for i in range(38, 88) for offset in helpers.VORTEX_EVENT_OFFSETS}
+            # K phase unknown until the final signal: retain flags for every boundary in RAM.
+            flag_cache = np.empty((count, *shape), np.uint16) if scene == 'karman' else None
+            for interval in range(plan['intervals']):
+                base.check_time(started, cap)
+                action = np.zeros(3, np.float32)
+                signals['policy_input'][interval] = obs
+                if role == 'controlled' and scene != 'steady':
+                    action = np.asarray(model.predict(obs, deterministic=True)[0], np.float32)
+                target, effective = wall_surface(scene, role, action), np.zeros(3)
+                sensor_sum, force_sum = np.zeros(6), np.zeros(6)
+                for native in range(plan['si']):
+                    if native % 1000 == 0:
+                        base.check_time(started, cap)
+                    y, v = (0., 0.) if role == 'target' else structure.candidate(force)
+                    check_structure(y, v, structure.y)
+                    effective = .9**.25*effective + (1.-.9**.25)*target
+                    omega = effective/20.
+                    for body, w in zip(bodies, omega):
+                        sim.set_body(body, omega=float(w), offset=(0., y), vy=y-structure.y)
+                    sim.run(1, zero_obs=True, sync_obs=True)
+                    if bodies:
+                        diag = sim.moving_diagnostics()
+                        cme, wen = diag[:, columns['cme_y']], diag[:, columns['wen_y']]
+                        node = diag[:, columns['cover_jy']]-diag[:, columns['uncover_jy']]
+                        counters = [diag[:, columns[key]].sum() for key in ('eq_only', 'no_donor', 'crossings')]
+                        packed = np.array([value for body in bodies for value in sim.read_force(body)[:2]])
+                    else:
+                        cme, wen, node, counters, packed = np.zeros(3), np.zeros(3), np.zeros(3), [0., 0., 0.], np.zeros(6)
+                    sensors = np.array([value for probe in probes for value in sim.read_sensor(probe)[:2]])
+                    row = np.array([completed+1, y, v, *cme, *wen, *node, *counters, *omega])
+                    if not np.isfinite(row).all() or counters[1] != 0 or not np.isfinite([packed, sensors]).all():
+                        raise RuntimeError('nonfinite native diagnostics/observations or no_donor')
+                    rows[completed] = row
+                    completed += 1
+                    structure.y, structure.v = y, v
+                    force = float(cme.sum(dtype=np.float64)+wen.sum(dtype=np.float64))
+                    sensor_sum += sensors
+                    force_sum += packed
+                check_clock(sim.stepper, plan['release_clock'], completed)
+                sensors, packed = sensor_sum/plan['si'], force_sum/plan['si']
+                raw = np.r_[sensors*156., packed*2.].astype(np.float32)
+                if norm is not None:
+                    raw, obs, clipped = base.observation(sensors, packed, 2., norm)
+                    report['clipped_channels'] = report.get('clipped_channels', 0)+clipped
+                for key, value in [('sensor_area_mean', sensors), ('sensors', sensors*156.),
+                                   ('forces', packed*2.), ('actions', action), ('observations', obs),
+                                   ('raw_observation', raw)]:
+                    signals[key][interval] = value
+                if interval >= plan['field_start']:
+                    index = interval-plan['field_start']
+                    clean, flags, fluid = clean_field(sim.get_macroscopic(), sim.get_flags())
+                    if not fluid[roi].all():
+                        raise RuntimeError('prescribed downstream ROI contains nonfluid cells')
+                    cache[index] = clean
+                    common_fluid &= fluid
+                    roi_cache[index] = clean[(slice(None), *roi)]
+                    if flag_cache is not None:
+                        flag_cache[index] = flags
+                    elif index in retain or index == count-1:
+                        selected_flags[index] = flags
+                completed_si += 1
+                if completed_si % 10 == 0 or completed_si == plan['intervals']:
+                    rows.flush()
+                    print(json.dumps(dict(event='progress', scene=scene, role=role, interval=completed_si,
+                                          absolute_clock=sim.stepper.step_count, y_D=y/D)), flush=True)
+            times = np.arange(plan['field_start']+1, plan['intervals']+1, dtype=np.int64)*plan['si']
+            selected, selection, event = select_fields(scene, plan, helpers, times,
+                                                       signals['sensors'][plan['field_start']:], cache, event)
+            indices = selected['selected_index']
+            selected['flags'] = flag_cache[indices] if flag_cache is not None else np.stack([selected_flags[int(i)] for i in indices])
+            selected['common_fluid_all_boundaries'] = common_fluid
+            selected['absolute_steps'] = selected['selected_response_steps']+plan['release_clock']
+            np.savez_compressed(directory/'fields.npz', **selected)
+            report.update(selection=selection, field_response_steps=times,
+                          field_absolute_steps=times+plan['release_clock'],
+                          field_window_tU_D=[int(times[0])*U/D, int(times[-1])*U/D],
+                          status='smoke_completed_not_physical' if plan['smoke'] else 'duration_completed_not_qualified')
+            roi_mask = common_fluid[roi].copy()
+        report['motion_summary'] = metrics.motion_summary(rows[max(0, completed-640000):completed, 1]/D)
+        return dict(roi=roi_cache, mask=roi_mask, times=times, event=event)
+    except BaseException as exc:
+        report['status'] = f'error: {type(exc).__name__}: {exc}'
+        raise
+    finally:
+        report.update(completed_native_steps=completed, completed_intervals=completed_si)
+        np.savez_compressed(directory/'signals.npz', response_steps=np.arange(1, completed_si+1)*plan['si'],
+                            absolute_steps=plan['release_clock']+np.arange(1, completed_si+1)*plan['si'],
+                            **{key: value[:completed_si] for key, value in signals.items()})
+        if completed:
+            report['native_archive'] = compress_native(directory/'native.npy', rows, completed)
+        write_json(directory/'report.json', report)
+        del cache, rows
+        gc.collect()
+
+
+def score_roles(roles, metrics):
+    target = roles['target']
+    common = np.logical_and.reduce([role['mask'] for role in roles.values()])
+    if not common.all():
+        raise ValueError('prescribed downstream ROI must be entirely common fluid')
+    result = dict(definition='sqrt(mean common ROI (dux^2+duy^2)/U^2)',
+                  alignment='same absolute native clock; no lag/phase alignment',
+                  roi_xy_inclusive=[1360, 2160, 311, 710], support_cells=int(common.sum()),
+                  support_coverage_fraction=1.)
+    for name in ('zero', 'controlled'):
+        role = roles[name]
+        if not np.array_equal(role['times'], target['times']):
+            raise ValueError('unmatched response clocks')
+        errors = np.array([metrics.field_error(*a, *b, U, common) for a, b in zip(role['roi'], target['roi'])])
+        # FP64 reduction without duplicating entire history.
+        mean = [a.mean(axis=0, dtype=np.float64) for a in (role['roi'], target['roi'])]
+        result[name] = dict(same_clock_instantaneous_E=errors,
+                            mean_instantaneous_E=float(errors.mean()),
+                            E_of_fixed_window_mean=metrics.field_error(*mean[0], *mean[1], U, common))
+    return result, common
+
+
+def source_manifest(base, helpers):
+    files = [Path(__file__), ROOT/'research/fig2_materials/metrics.py', Path(base.__file__),
+             ROOT/'research/steady_pinball/run.py', ROOT/'research/steady_pinball/config.json',
+             Path(helpers.__file__), Path(helpers.__file__).parent/'legacy_test/acquire.py',
+             Path(helpers.__file__).parent/'legacy_test/core/dtw_metrics.py']
+    files += sorted((ROOT/'src/CelerisLab').rglob('*.py'))
+    for suffix in ('*.cu', '*.cuh', '*.h'):
+        files += sorted((ROOT/'src/CelerisLab').rglob(suffix))
+    return {str(path.resolve()): sha(path) for path in files}
+
+
+def main(argv=None):
+    ap = argparse.ArgumentParser(description=__doc__)
+    ap.add_argument('--scene', choices=('karman', 'steady', 'taylor'), required=True)
+    ap.add_argument('--out', type=Path, required=True)
+    ap.add_argument('--smoke', action='store_true', help='32-step warmups; 3 SI of 32; no physical score')
+    ap.add_argument('--D', type=int, choices=(40,), default=40)
+    ap.add_argument('--model', type=Path)
+    ap.add_argument('--norm', type=Path)
+    ap.add_argument('--release-dir', type=Path, default=K_RELEASE)
+    ap.add_argument('--acquisition', type=Path, default=ACQUISITION)
+    ap.add_argument('--wall-seconds', type=float, default=43200.)
+    args = ap.parse_args(argv)
+    if not np.isfinite(args.wall_seconds) or args.wall_seconds <= 0:
+        ap.error('positive finite wall cap required')
+    if args.scene == 'steady' and (args.model or args.norm):
+        ap.error('steady constant control has no PPO or normalizer')
+    out = args.out.resolve()
+    if out == ROOT or ROOT in out.parents:
+        ap.error('fresh output must be outside checkout')
+    sys.path.insert(0, str(ROOT))
+    from research.karman_viv import run as base
+    metrics = load_file('_fig2_metrics', ROOT/'research/fig2_materials/metrics.py')
+    helpers = load_helpers(args.acquisition)
+    model = norm = None
+    inputs = {}
+    if args.scene != 'steady':
+        model_path, model_sha, norm_path, norm_sha = ASSETS[args.scene]
+        model_path, norm_path = args.model or model_path, args.norm or norm_path
+        if sha(model_path) != model_sha or sha(norm_path) != norm_sha:
+            ap.error('model/normalizer hashes differ from frozen inputs')
+        norm = json.loads(norm_path.read_text())
+        inputs = dict(model_path=model_path, model_sha256=model_sha,
+                      norm_path=norm_path, norm_sha256=norm_sha)
+        model = base.load_policy(model_path)
+    sources = source_manifest(base, helpers)
+    source = subprocess.check_output(['git', '-C', str(ROOT), 'rev-parse', 'HEAD'], text=True).strip()
+    dirty = subprocess.check_output(['git', '-C', str(ROOT), 'status', '--porcelain'], text=True)
+    plan = schedule(args.scene, args.smoke)
+    admit_memory(args.smoke)
+    out.mkdir(parents=True, exist_ok=False)
+    cfgpath = out/'config.json'
+    write_json(cfgpath, base.config(40))
+    receipt = dict(scene=args.scene, plan=plan, D=D, U=U, nu=NU, actual_Re_D=50.,
+                   historical_case='karman_re100' if args.scene == 'karman' else args.scene,
+                   expected_baseline=BASELINE, source=source, git_dirty=bool(dirty), git_status=dirty,
+                   source_sha256=sources, inputs=inputs, centres=CENTRES, probes=PROBES,
+                   environment=dict(python=sys.version, platform=platform.platform(), numpy=np.__version__,
+                                    conda_env=os.environ.get('CONDA_DEFAULT_ENV'),
+                                    cuda_visible_devices=os.environ.get('CUDA_VISIBLE_DEVICES')),
+                   acquisition_event_offsets=list(helpers.VORTEX_EVENT_OFFSETS),
+                   event_selector='source _select_vortex_event AST verbatim; actual acquisition offsets',
+                   structure=dict(mass_ratio=4., fn_native=2.5e-5, fn_D_U=.2, zeta=0., travel_D=1., speed=.025),
+                   rhs='previous native-step sum CME_y+Wen_y; initial force last static native step',
+                   policy='SI mean packed CME *2; area-mean sensors *156; frozen norm; first input zeros',
+                   control='steady constant (0,-5.1,5.1)*U; K scale8 or Taylor scale4, bias(0,-4,4)*U',
+                   ema=dict(alpha=.9**.25, reset='every SI'),
+                   storage='one full-grid float32 RAM cache at a time; compact ROI histories in RAM; no full history field disk',
+                   limitation='Modern moving-plant experiment; no Legacy numerical equivalence or steady/convergence claim')
+    write_json(out/'case.json', receipt)
+    started, status = time.monotonic(), 'running'
+    try:
+        seed, force, release = make_release(args.scene, plan, cfgpath, base, out, started,
+                                            args.wall_seconds, args.release_dir)
+        receipt['release'] = release
+        receipt['release'].update(clock=int(seed[3]), initial_force=force,
+                                  current_sha256=base.array_sha(seed[0]),
+                                  spare_sha256=base.array_sha(seed[1]), flags_sha256=base.array_sha(seed[2]))
+        write_json(out/'case.json', receipt)
+        roles, event = {}, None
+        order = ('controlled', 'zero', 'target') if args.scene == 'taylor' else ('target', 'zero', 'controlled')
+        for role in order:
+            roles[role] = run_role(args.scene, role, plan, cfgpath, base, helpers, metrics,
+                                    seed, force, model, norm, out, started, args.wall_seconds, event)
+            if args.scene == 'taylor':
+                event = roles[role]['event']
+        if not args.smoke:
+            scores, support = score_roles(roles, metrics)
+            write_json(out/'same_clock_errors.json', scores)
+            np.savez_compressed(out/'common_roi_support.npz', support=support)
+        if source_manifest(base, helpers) != sources:
+            raise RuntimeError('source files changed during collection')
+        status = 'smoke_completed_not_physical' if args.smoke else 'three_roles_completed_not_qualified'
+    except BaseException as exc:
+        status = f'error: {type(exc).__name__}: {exc}'
+        raise
+    finally:
+        write_json(out/'terminal.json', dict(status=status, elapsed_s=time.monotonic()-started,
+                                            physical_score_claim=False if args.smoke else None))
+
+
+if __name__ == '__main__':
+    main()

research/fig2_materials/metrics.py

@@ -0,0 +1,94 @@
+"""Pure NumPy Fig2 primitives; no solver, phase selection, or probe DTW.
+
+Adapted formula and inclusive ROI bounds from DynamisLab
+src/drl_pinball/eval/wake_l2.py::mean_field_metric, roi_slices, inspected
+2026-11-02 at HEAD 609116ec26f19766f5e9ce2b0df427d4e24faae8
+(working-tree source, not asserted clean). The explicit common-fluid/finite
+mask and motion summaries are additions. legacy_test/metrics.py::cycle_similarity
+was inspected only: its six probe channels are not a field recovery metric.
+"""
+import math
+import numpy as np
+
+
+def roi_slices(shape, *, D, xs, y0=None, downstream_solid_extent):
+    """Lab-fixed ROI: sensor-relative x/D=-6..14, y/D=-5..5, inclusive.
+
+    Equivalent to pinball-front-relative x/D=4..24 when xs=front+10D.
+    Supply actual geometry, e.g. D=40, xs=1600, y0=510.5.
+    """
+    ny, nx = shape
+    y0 = (ny - 1) / 2 if y0 is None else float(y0)
+    if D <= 0 or not all(np.isfinite(v) for v in (D, xs, y0, downstream_solid_extent)):
+        raise ValueError("invalid geometry")
+    x0, x1 = math.ceil(xs - 6 * D), math.floor(xs + 14 * D)
+    lo, hi = math.ceil(y0 - 5 * D), math.floor(y0 + 5 * D)
+    if x0 < 0 or lo < 0 or x1 >= nx or hi >= ny:
+        raise ValueError("ROI outside grid")
+    if x0 <= downstream_solid_extent:
+        raise ValueError("ROI not downstream of solid extent")
+    return slice(lo, hi + 1), slice(x0, x1 + 1)
+
+
+def common_fluid_roi_mask(fluid_masks, fields, roi):
+    """Intersect explicit boolean fluid masks, ROI, and finite 2-D fields.
+
+    Pass every compared role/component/sample to freeze a common support over
+    time. Zero velocity never implies solid. Masks must share the lab grid.
+    """
+    fields = [np.asarray(a, dtype=float) for a in fields]
+    masks = [np.asarray(a) for a in fluid_masks]
+    if not fields or not masks:
+        raise ValueError("fields and explicit fluid masks required")
+    shape = fields[0].shape
+    if len(shape) != 2 or any(a.shape != shape for a in fields + masks):
+        raise ValueError("equal-shape 2-D arrays required")
+    if any(a.dtype != np.dtype(bool) for a in masks):
+        raise ValueError("fluid masks must be boolean")
+    support = np.zeros(shape, dtype=bool)
+    ys, xs = roi
+    for s, n in zip((ys, xs), shape):
+        if not isinstance(s, slice) or s.step not in (None, 1) or not 0 <= (s.start or 0) < (n if s.stop is None else s.stop) <= n:
+            raise ValueError("invalid ROI slices")
+    support[ys, xs] = True
+    for a in masks:
+        support &= a
+    for a in fields:
+        support &= np.isfinite(a)
+    if not support.any():
+        raise ValueError("empty common finite fluid ROI")
+    return support
+
+
+def field_error(rux, ruy, tux, tuy, u0, support):
+    """Eq1: sqrt(mean_support((dux**2 + duy**2)/U0**2)).
+
+    For E(mean field), pass arithmetic mean fields from matched windows.
+    Mean instantaneous E requires separate calls per snapshot and then a mean;
+    these operations are not interchangeable. No target-energy denominator.
+    """
+    a = [np.asarray(v, dtype=float) for v in (rux, ruy, tux, tuy)]
+    support = np.asarray(support)
+    if a[0].ndim != 2 or any(v.shape != a[0].shape for v in a) or support.shape != a[0].shape or support.dtype != np.dtype(bool):
+        raise ValueError("equal-shape 2-D fields and boolean support required")
+    if not np.isfinite(u0) or u0 <= 0 or not support.any():
+        raise ValueError("positive finite U0 and nonempty support required")
+    if any(not np.isfinite(v[support]).all() for v in a):
+        raise ValueError("nonfinite field on fixed support")
+    du = (a[0][support] - a[2][support]) / u0
+    dv = (a[1][support] - a[3][support]) / u0
+    return float(np.sqrt(np.mean(du * du + dv * dv)))
+
+
+def motion_summary(displacement):
+    """Arithmetic mean, centered RMS, half peak-to-peak along time axis 0.
+
+    Input is displacement in caller-declared units; no implicit D scaling.
+    Use the same retained window as field evaluation. No sinusoid assumption.
+    """
+    a = np.asarray(displacement, dtype=float)
+    if a.ndim == 0 or a.shape[0] == 0 or not np.isfinite(a).all():
+        raise ValueError("nonempty finite time samples required")
+    mean = a.mean(axis=0)
+    return {"mean": mean, "centered_rms": np.sqrt(np.mean((a - mean) ** 2, axis=0)),
+            "half_range": (a.max(axis=0) - a.min(axis=0)) / 2}

research/karman_viv/run.py

@@ -0,0 +1,332 @@
+"""One paired zero/PPO Karman pinball VIV experiment; see docs/karman_viv.md.
+
+Research-only adapter: static precursor -> two fresh, identical moving plants.
+No moving checkpoint/rollback, policy retraining, or force filtering.
+"""
+from __future__ import annotations
+
+import argparse
+import hashlib
+import json
+from pathlib import Path
+import subprocess
+import sys
+import time
+
+import numpy as np
+
+ROOT = Path(__file__).resolve().parents[2]
+sys.path.insert(0, str(ROOT))
+from research.steady_pinball.run import Structure, config, scaled_points, CENTRES, PROBES, COLUMNS
+
+BIAS = np.array([0., -4., 4.])
+MODEL_SHA = '148240c1dcb8b11d8e5a0a9e991f6874c4e799c9744707bf0fb32a8efeb0468d'
+NORM_SHA = '534fda818c030a56960e95d0f3ac1e350d887d03db7f80f1d6453195a502b7f0'
+RELEASE_SOURCE = '516975ef6482f983ace347cb10de7450d2fbd680'
+STEP_COLUMNS = (*COLUMNS, 'omega_F', 'omega_U', 'omega_L')
+
+
+def sha(path):
+    h = hashlib.sha256()
+    with open(path, 'rb') as f:
+        for block in iter(lambda: f.read(1024*1024), b''):
+            h.update(block)
+    return h.hexdigest()
+
+
+def array_sha(a):
+    return hashlib.sha256(memoryview(np.ascontiguousarray(a)).cast('B')).hexdigest()
+
+
+def observation(sensor_mean, force_mean, scale, norm):
+    """Fine-grid area means/CME -> trained D20 sensor sums/forces -> frozen norm."""
+    sensors = np.asarray(sensor_mean, np.float64)*scale*78.
+    forces = np.asarray(force_mean, np.float64)*scale
+    raw = np.asarray([*sensors, *forces], np.float32)
+    preclip = np.hstack((raw[6:]/norm['force_norm_fact'],
+                        (raw[:6]-np.asarray(norm['sens_deviation'], np.float32)) /
+                        np.asarray(norm['sens_norm_fact'], np.float32)))
+    if not np.isfinite(preclip).all():
+        raise RuntimeError('nonfinite policy observation')
+    return raw, np.clip(preclip, -1, 1).astype(np.float32), int((np.abs(preclip)>1).sum())
+
+
+def wall_target(action, scale):
+    return (8.*np.asarray(action)+BIAS)*(.01/scale)
+
+
+def ema(eff, target, scale):
+    alpha = .9**(1./scale**2)
+    return alpha*eff + (1.-alpha)*target
+
+
+def load_policy(model):
+    # Same inference spaces as legacy_test.core.model_loader, without importing it.
+    import gymnasium as gym
+    from stable_baselines3 import PPO
+    class InferenceEnv(gym.Env):
+        observation_space = gym.spaces.Box(-1., 1., (12,), np.float32)
+        action_space = gym.spaces.Box(-1., 1., (3,), np.float32)
+        def reset(self, *, seed=None, options=None):
+            super().reset(seed=seed)
+            return np.zeros(12, np.float32), {}
+        def step(self, action):
+            return np.zeros(12, np.float32), 0., False, False, {}
+    return PPO.load(str(model), env=InferenceEnv(), device='cpu')
+
+
+def add_scene(sim, scale, moving=False, pinball=True):
+    from CelerisLab.body.motion import Translation
+    sim.add_body('circle', scaled_points(((200., 255.5),), scale)[0], 20.*scale)
+    probes = [sim.add_body('sensor', c, 5.*scale) for c in scaled_points(PROBES, scale)]
+    bodies = []
+    if pinball:
+        envelope = Translation((0., 0., -20.*scale, 20.*scale), .05/scale, shared=True) if moving else None
+        bodies = [sim.add_body('circle', c, 10.*scale, translation=envelope)
+                  for c in scaled_points(CENTRES, scale)]
+    return probes, bodies
+
+
+def capture_seed(sim):
+    """Capture both FP32 buffers at fixed geometry; never a moving-state checkpoint."""
+    import pycuda.driver as cuda
+    with sim.ctx.activate():
+        sim.stream.synchronize()
+        buffers = [np.empty_like(sim.field.ddf) for _ in range(2)]
+        for dst, src in zip(buffers, (sim.field.ddf_gpu, sim.field.temp_gpu)):
+            cuda.memcpy_dtoh(dst, src)
+    return (*buffers, sim.get_flags(), sim.stepper.step_count)
+
+
+def install_seed(sim, seed):
+    """Only a new, zero-offset instance may receive the static release state."""
+    import pycuda.driver as cuda
+    current, spare, flags, clock = seed
+    if sim.stepper.step_count != 0 or any(a.shape != sim.field.ddf.shape or a.dtype != sim.field.ddf.dtype
+                                         for a in (current, spare)):
+        raise ValueError('release transfer requires a fresh matching field')
+    if not np.array_equal(sim.get_flags(), flags):
+        raise ValueError('static/moving release flags differ')
+    for obj in sim.bodies.objects:
+        if obj.motion is not None and (obj.motion.applied != (0., 0.) or obj.motion.pending is not None):
+            raise ValueError('release transfer requires zero geometry and no target')
+    with sim.ctx.activate():
+        sim.stream.synchronize()
+        cuda.memcpy_htod(sim.field.ddf_gpu, current)
+        cuda.memcpy_htod(sim.field.temp_gpu, spare)
+    sim.stepper._step_count = clock  # one-time research initialization, not restore()
+    sim.field.invalidate_host_ddf_cache()
+    copied = capture_seed(sim)
+    if copied[3] != clock or any(not np.array_equal(a, b, equal_nan=True) for a, b in zip(copied[:3], seed[:3])):
+        raise RuntimeError('release transfer did not preserve both buffers/flags')
+    return dict(current=array_sha(copied[0]), spare=array_sha(copied[1]), flags=array_sha(copied[2]), clock=copied[3])
+
+
+def check_time(start, cap):
+    if time.monotonic()-start > cap:
+        raise TimeoutError('paired experiment wall cap')
+
+
+def warmup(sim, count, start, cap):
+    while count:
+        check_time(start, cap)
+        n = min(count, 10000)
+        sim.run(n)
+        count -= n
+
+
+def precursor(cfgpath, scale, warm, si, targets, out, start, cap):
+    from CelerisLab import Simulation
+    with Simulation(str(cfgpath)) as sim:
+        probes, _ = add_scene(sim, scale, pinball=False)
+        sim.initialize()
+        warmup(sim, warm, start, cap)
+        target = []
+        for _ in range(targets):
+            check_time(start, cap)
+            sim.run(si)
+            target.append([v*scale*78. for p in probes for v in sim.read_sensor(p)[:2]])
+        np.save(out/'target.npy', np.asarray(target, np.float32))
+        first = sim.bodies.count
+        for c in scaled_points(CENTRES, scale):
+            sim.add_body('circle', c, 10.*scale)
+        sim.sync_bodies()
+        bodies = list(range(first, first+3))
+        warmup(sim, warm-1, start, cap)
+        sim.run(1)
+        force = sum(float(sim.read_force(p)[1]) for p in bodies)
+        seed = capture_seed(sim)
+        np.savez(out/'release.npz', current=seed[0], spare=seed[1], flags=seed[2], clock=seed[3], force=force)
+        print(json.dumps(dict(event='precursor_complete', clock=seed[3], force=force,
+                              elapsed_s=time.monotonic()-start)), flush=True)
+        return seed, force
+
+
+def run_arm(name, cfgpath, scale, si, intervals, seed, initial_force, model, norm, out, start, cap, *, mass_ratio=2.):
+    from CelerisLab import Simulation
+    from CelerisLab.body.motion import MovingTables
+    cols = {k:i for i, k in enumerate(MovingTables.DIAG)}
+    diameter = 20.*scale
+    structure = Structure(mass_ratio, 1.e-4/scale**2, 0., diameter)
+    arm = out/name
+    arm.mkdir()
+    rows = np.lib.format.open_memmap(arm/'steps.npy', mode='w+', dtype=np.float64,
+                                    shape=(si*intervals, len(STEP_COLUMNS)))
+    rows[:] = np.nan
+    signals = {k:np.full((intervals, width), np.nan, np.float32) for k, width in
+               [('sensors',6), ('forces',6), ('actions',3), ('observations',12)]}
+    completed = completed_intervals = clipped = 0
+    status, transfer = 'running', None
+    force = initial_force
+    obs = np.zeros(12, np.float32)
+    arm_start = time.monotonic()
+    try:
+        with Simulation(str(cfgpath)) as sim:
+            probes, bodies = add_scene(sim, scale, moving=True)
+            sim.initialize()
+            transfer = install_seed(sim, seed)
+            for i in range(intervals):
+                action = -BIAS/8 if name == 'zero' else np.asarray(model.predict(obs, deterministic=True)[0], np.float32)
+                if action.shape != (3,) or not np.isfinite(action).all() or (np.abs(action)>1).any():
+                    raise RuntimeError('invalid policy action')
+                target, eff = wall_target(action, scale), np.zeros(3)
+                sensor_sum, force_sum = np.zeros(6), np.zeros(6)
+                for k in range(si):
+                    if k % 1000 == 0:
+                        check_time(start, cap)
+                    y1, v1 = structure.candidate(force)
+                    if not np.isfinite([y1,v1]).all() or abs(y1)>diameter:
+                        raise RuntimeError('nonfinite structure or |y| > D; no retuning')
+                    eff = ema(eff, target, scale)
+                    omega = eff/(diameter/2)
+                    for p, w in zip(bodies, omega):
+                        sim.set_body(p, omega=float(w), offset=(0., y1), vy=y1-structure.y)
+                    sim.run(1, zero_obs=True, sync_obs=True)
+                    d = sim.moving_diagnostics()
+                    cme, wen = d[:,cols['cme_y']], d[:,cols['wen_y']]
+                    node = d[:,cols['cover_jy']]-d[:,cols['uncover_jy']]
+                    row = (completed+1,y1,v1,*cme,*wen,*node,
+                           d[:,cols['eq_only']].sum(),d[:,cols['no_donor']].sum(),d[:,cols['crossings']].sum(),*omega)
+                    if not np.isfinite(row).all() or row[13] != 0:
+                        raise RuntimeError('nonfinite moving diagnostic or no_donor')
+                    f = np.array([v for p in bodies for v in sim.read_force(p)[:2]])
+                    s = np.array([v for p in probes for v in sim.read_sensor(p)[:2]])
+                    if not np.isfinite([f,s]).all():
+                        raise RuntimeError('nonfinite force/sensor observation')
+                    # Policy sees packed CME as in fixed-body training; structure sees CME+Wen.
+                    force_sum += f
+                    sensor_sum += s
+                    force = float(cme.sum(dtype=np.float64)+wen.sum(dtype=np.float64))
+                    rows[completed] = row
+                    completed += 1
+                    structure.y, structure.v = y1,v1
+                raw, obs, nclip = observation(sensor_sum/si, force_sum/si, scale, norm)
+                clipped += nclip
+                for key, val in [('sensors',raw[:6]),('forces',raw[6:]),('actions',action),('observations',obs)]:
+                    signals[key][i] = val
+                completed_intervals = i+1
+                if (i+1)%10 == 0 or i+1 == intervals:
+                    rows.flush()
+                    np.savez(arm/'signals.npz', **signals)
+                    print(json.dumps(dict(event='progress', arm=name, interval=i+1, step=completed,
+                                          tU_D=completed*.01/scale/diameter, y_D=y1/diameter,
+                                          elapsed_s=time.monotonic()-arm_start)), flush=True)
+            status = 'duration_completed_not_qualified'
+    except BaseException as exc:
+        status = f'error: {type(exc).__name__}: {exc}'
+        raise
+    finally:
+        rows.flush()
+        np.savez(arm/'signals.npz', **signals)
+        report = dict(status=status, completed_steps=completed, completed_intervals=completed_intervals,
+                      elapsed_s=time.monotonic()-arm_start, release_transfer=transfer, clipped_observation_channels=clipped,
+                      columns=STEP_COLUMNS, M=structure.mass,K=structure.stiffness,fn=structure.frequency,damping=0.,
+                      max_abs_y_D=float(np.max(np.abs(rows[:completed,1]))/diameter) if completed else None)
+        (arm/'report.json').write_text(json.dumps(report,indent=2)+'\n')
+
+
+def saved_release(directory):
+    """Read the retained m*=2 static release; not a moving-state restart."""
+    case = json.loads((directory/'case.json').read_text())
+    expected = dict(source=RELEASE_SOURCE, model_sha256=MODEL_SHA, norm_sha256=NORM_SHA,
+                    smoke=False, D=40, U=.005, nu=.004, warmup_each=2048000,
+                    target_intervals=150, SI=3200)
+    if any(case.get(k) != v for k, v in expected.items()):
+        raise ValueError('saved release case does not match the frozen precursor')
+    if json.loads((directory/'config.json').read_text()) != config(40):
+        raise ValueError('saved release config differs')
+    with np.load(directory/'release.npz') as data:
+        seed = (data['current'], data['spare'], data['flags'], int(data['clock']))
+        force = float(data['force'])
+    transfer = json.loads((directory/'zero/report.json').read_text())['release_transfer']
+    if seed[3] != transfer['clock'] or seed[3] != 4576000 or not np.isfinite(force):
+        raise ValueError('saved release clock/force differs')
+    if [array_sha(a) for a in seed[:3]] != [transfer[k] for k in ('current','spare','flags')]:
+        raise ValueError('saved release buffer/flags hash differs')
+    target = np.load(directory/'target.npy')
+    if target.shape != (150,6) or not np.isfinite(target).all():
+        raise ValueError('saved target is incomplete or nonfinite')
+    return seed, force
+
+
+def main():
+    ap = argparse.ArgumentParser(description=__doc__)
+    ap.add_argument('--model', type=Path, required=True)
+    ap.add_argument('--norm', type=Path, required=True)
+    ap.add_argument('--out', type=Path, required=True)
+    ap.add_argument('--smoke', action='store_true', help='short plumbing only; no physical conclusion')
+    ap.add_argument('--mass-ratio', type=float, choices=(2.,4.), default=2.,
+                    help='4 is the single declared exploratory candidate; fn and damping unchanged')
+    ap.add_argument('--release-dir', type=Path, help='retained m*=2 precursor; skip warmup')
+    ap.add_argument('--wall-seconds', type=float, default=14400.)
+    args = ap.parse_args()
+    if not np.isfinite(args.wall_seconds) or args.wall_seconds<=0:
+        ap.error('positive finite wall cap required')
+    if args.release_dir is not None and args.smoke:
+        ap.error('saved physical release cannot be used for the plumbing smoke')
+    if sha(args.model)!=MODEL_SHA or sha(args.norm)!=NORM_SHA:
+        ap.error('model/normalizer must match frozen provenance')
+    out = args.out.resolve()
+    if out == ROOT or ROOT in out.parents:
+        ap.error('output must be outside checkout')
+    out.mkdir(parents=True, exist_ok=False)
+    norm = json.loads(args.norm.read_text())
+    model = load_policy(args.model)
+    scale = 2.
+    warm, si, targets, intervals = (32,8,2,4) if args.smoke else (2048000,3200,150,400)
+    cfgpath = out/'config.json'
+    cfgpath.write_text(json.dumps(config(40),indent=2)+'\n')
+    receipt = dict(source=subprocess.check_output(['git','-C',str(ROOT),'rev-parse','HEAD'],text=True).strip(),
+                   model_sha256=sha(args.model),norm_sha256=sha(args.norm),smoke=args.smoke,D=40,U=.005,nu=.004,
+                   warmup_each=warm,target_intervals=targets,SI=si,response_intervals=intervals,
+                   response_steps=si*intervals,response_tU_D=si*intervals*.005/40,
+                   rhs='previous native step sum CME_y+Wen_y; policy uses interval packed CME',
+                   probe_frame='lab-fixed; fine area mean * 2 * trained 78 cells',
+                   ema_native_alpha=.9**.25, structural_damping=0., mass_ratio=args.mass_ratio,
+                   exploratory=args.mass_ratio==4., columns=STEP_COLUMNS)
+    loaded = None
+    if args.release_dir is not None:
+        loaded = saved_release(args.release_dir)
+        cfgpath.write_bytes((args.release_dir/'config.json').read_bytes())
+        (out/'target.npy').write_bytes((args.release_dir/'target.npy').read_bytes())
+        receipt.update(release_dir=str(args.release_dir.resolve()), release_source=RELEASE_SOURCE,
+                       release_sha256=sha(args.release_dir/'release.npz'),
+                       warmup_executed=0, target_sha256=sha(out/'target.npy'))
+    (out/'case.json').write_text(json.dumps(receipt,indent=2)+'\n')
+    start = time.monotonic()
+    status = 'running'
+    try:
+        seed, force = loaded if loaded is not None else precursor(cfgpath,scale,warm,si,targets,out,start,args.wall_seconds)
+        for name in ('zero','controlled'):
+            run_arm(name,cfgpath,scale,si,intervals,seed,force,model,norm,out,start,args.wall_seconds,
+                    mass_ratio=args.mass_ratio)
+        status = 'smoke_completed' if args.smoke else 'pair_completed_not_qualified'
+    except BaseException as exc:
+        status = f'error: {type(exc).__name__}: {exc}'
+        raise
+    finally:
+        (out/'terminal.json').write_text(json.dumps(dict(status=status,elapsed_s=time.monotonic()-start),indent=2)+'\n')
+
+
+if __name__ == '__main__':
+    main()

research/steady_pinball/config.json

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

research/steady_pinball/run.py

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

research/viv_re150/README.md

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

research/viv_re150/run.py

@@ -1,8 +1,9 @@
-"""One-shot predeclared VIV run (docs/validation/leontini_viv_predeclare.md).
+"""One-shot predeclared VIV run (docs/validation/viv_re150_predeclare.md).
 
-Re 150, U* = 5, m/(rho D^2) = 2, zeta = 0, D = 40 on a 75 D x 60 D native
-channel (uniform inlet, NEQ outlet, free-slip y walls), cylinder at
-(30 D, 30 D). Structural RHS = CME + Wen (diagnostics: raw CME, node terms).
+Re 150, U* = 5, m/(rho D^2) = 2, zeta = 0, D = 40 on a nominal 75 D x 60 D
+lattice (uniform inlet, NEQ outlet, free-slip y walls). The cylinder is at
+(30 D, (ny - 1)/2); y-wall planes are .5 and ny - 1.5 in lattice coordinates.
+Structural RHS = CME + Wen (diagnostics: raw CME, node terms).
 Research driver only: Simulation API, no production change.
 
     PYTHONPATH=src python research/viv_re150/run.py --out /path/outside/checkout
@@ -65,7 +66,8 @@ def main():
         step_cap = args.steps
     centre = (30.0 * Dn, 0.5 * (ny - 1))
 
-    from CelerisLab.body.motion import Translation
+    from CelerisLab.body.motion import MovingTables, Translation
+    col = {k: MovingTables.DIAG.index(k) for k in ("cme_y", "wen_y", "cover_jy", "uncover_jy", "no_donor")}
     from CelerisLab.simulation import Simulation
     with tempfile.NamedTemporaryFile("w", suffix=".json", delete=False) as fh:
         json.dump(cfg, fh)
@@ -95,9 +97,9 @@ def main():
             sim.run(1)                             # device failures raise and mark the Simulation failed
             d = sim.moving_diagnostics()[0]
             f_obs = sim.read_force(body)[1]
-            f_rhs = float(d[1] + d[3])             # CME + Wen (y)
-            row = (n, y1, v1, f_rhs, float(d[1]), float(d[3]), float(f_obs),
-                   float(d[7] - d[10]), float(d[13]))
+            f_rhs = float(d[col["cme_y"]] + d[col["wen_y"]])   # CME + Wen (y)
+            row = (n, y1, v1, f_rhs, float(d[col["cme_y"]]), float(d[col["wen_y"]]), float(f_obs),
+                   float(d[col["cover_jy"]] - d[col["uncover_jy"]]), float(d[col["no_donor"]]))
             if not viv.finite(row):                # every recorded output, not only the RHS
                 status = "fail_nonfinite_output"
                 break

research/viv_re150/viv.py

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

src/CelerisLab/body/motion.py

@@ -1,10 +1,10 @@
 # CelerisLab/body/motion.py
-"""Host state of a translating circle (dense double buffer, D2Q9 only).
+"""Translation targets and device tables (dense double buffer, D2Q9 only).
 
 Contract: docs/moving_bodies.md, sections 4a, 5 and 7. The applied geometry
 (last successfully committed offset) and an optional pending target are two
-separate facts; centres are always ``base + offset``. Nothing here touches
-the device: the caller supplies upload/step/synchronize callables.
+separate facts; centres are always ``base + offset``. MovingCircle owns host
+targets, MovingTables owns device records, and Simulation commits completed steps.
 """
 from __future__ import annotations
 
@@ -14,14 +14,11 @@ from typing import Callable, Optional, Tuple
 
 MAX_STEP = 0.5     # Euclidean centre displacement per step (lattice cells)
 BAND_MARGIN = 2    # cells added around the swept circle AABB
+SHARED_CLEARANCE = 4.0  # > two diagonal links + MAX_STEP; docs/steady_pinball_fiv.md
 R_MIN = 4.0        # smallest radius covered by the donor enumeration
 
 Vec = Tuple[float, float]
 
-# D2Q9 non-rest directions, paired order of lbm/descriptors.
-_CX = (1, -1, 0, 0, 1, -1, 1, -1)
-_CY = (0, 0, 1, -1, 1, -1, -1, 1)
-
 
 def _vec(value, name: str) -> Vec:
     x, y = (float(v) for v in value)
@@ -35,6 +32,7 @@ class Translation:
     """Declared motion envelope: offset box (x_min, x_max, y_min, y_max) and speed bound."""
     offset_box: Tuple[float, float, float, float]
     max_speed: float
+    shared: bool = False  # same envelope instance: synchronous transverse motion
 
     def __post_init__(self):
         box = tuple(float(v) for v in self.offset_box)
@@ -45,10 +43,29 @@ class Translation:
             raise ValueError("offset_box must contain the zero offset")
         if not 0.0 < speed <= MAX_STEP:
             raise ValueError(f"max_speed must satisfy 0 < max_speed <= {MAX_STEP}")
+        if self.shared and box[:2] != (0.0, 0.0):
+            raise ValueError("shared translation is transverse only (x offsets must be zero)")
         object.__setattr__(self, "offset_box", box)
         object.__setattr__(self, "max_speed", speed)
 
 
+def validate_shared_translations(circles):
+    """Check common old/target offsets and wall velocities before device work.
+
+    Per-body setters may stage members in any order. Envelope identity (not
+    value equality) defines a group; a missing target means hold at applied.
+    """
+    groups = {}
+    for motion in circles:
+        if not motion.translation.shared:
+            continue
+        target = motion.pending if motion.pending is not None else (motion.applied, (0.0, 0.0))
+        state = (motion.applied, target)
+        key = id(motion.translation)
+        if groups.setdefault(key, state) != state:
+            raise ValueError("shared translation needs identical applied/target offsets and velocities")
+
+
 class MovingCircle:
     """Applied offset, pending target and the fail-closed commit of one circle."""
 
@@ -76,6 +93,8 @@ class MovingCircle:
         x0, x1, y0, y1 = self.translation.offset_box
         if not (x0 <= offset[0] <= x1 and y0 <= offset[1] <= y1):
             raise ValueError(f"offset {offset} outside declared box {self.translation.offset_box}")
+        if self.translation.shared and velocity[0] != 0.0:
+            raise ValueError("shared translation needs zero x velocity")
         step = math.hypot(offset[0] - self.applied[0], offset[1] - self.applied[1])
         if step > MAX_STEP:
             raise ValueError(f"per-step displacement {step} exceeds {MAX_STEP}")
@@ -98,23 +117,6 @@ class MovingCircle:
             offset, velocity = self.applied, (0.0, 0.0)
         return self.centre(self.applied), self.centre(offset), velocity
 
-    def advance(self, upload: Callable, step: Callable, synchronize: Callable) -> None:
-        """Upload the record, run one step and commit only after a successful sync.
-
-        Any exception after the upload is attempted leaves device geometry
-        unknown: the circle is marked failed, never reports applied state again,
-        and the owner's ``on_failure`` hook is called.
-        """
-        record = self.record()
-        try:
-            upload(record)
-            step()
-            synchronize()
-        except BaseException as exc:
-            self.fail(exc)
-            raise
-        self.commit()
-
     def commit(self) -> None:
         """Applied := pending offset (caller has synchronized a successful step)."""
         if self.pending is not None:
@@ -138,43 +140,6 @@ class MovingCircle:
             raise RuntimeError("moving circle state is unknown after a failed step")
 
 
-def _inside(x: int, y: int, c: Vec, r: float) -> bool:
-    return (x - c[0]) ** 2 + (y - c[1]) ** 2 < r * r
-
-
-def empty_donor_cells(c_old: Vec, c_new: Vec, radius: float):
-    """Fresh cells (inside c_old, outside c_new) whose donor set V(x) is empty."""
-    r = radius
-    xs = range(math.floor(min(c_old[0], c_new[0]) - r) - 1, math.ceil(max(c_old[0], c_new[0]) + r) + 2)
-    ys = range(math.floor(min(c_old[1], c_new[1]) - r) - 1, math.ceil(max(c_old[1], c_new[1]) + r) + 2)
-    empty = []
-    for x in xs:
-        for y in ys:
-            if not _inside(x, y, c_old, r) or _inside(x, y, c_new, r):
-                continue
-            if not any(not _inside(x + cx, y + cy, c_old, r) and not _inside(x + cx, y + cy, c_new, r)
-                       for cx, cy in zip(_CX, _CY)):
-                empty.append((x, y))
-    return empty
-
-
-def enumerate_donor_coverage(radius: float, positions: int = 8, directions: int = 16) -> int:
-    """Count empty-donor fresh cells over sub-cell centres and moves of MAX_STEP.
-
-    Evidence for the declared bounds (radius >= R_MIN, |dc|_2 <= MAX_STEP),
-    not a proof; the device keeps the ``no_donor`` counter.
-    """
-    count = 0
-    for i in range(positions):
-        for j in range(positions):
-            c_old = (i / positions, j / positions)
-            for k in range(directions):
-                a = 2.0 * math.pi * k / directions
-                c_new = (c_old[0] + MAX_STEP * math.cos(a), c_old[1] + MAX_STEP * math.sin(a))
-                count += len(empty_donor_cells(c_old, c_new, radius))
-    return count
-
-
 class MovingTables:
     """Device tables of the translating circles for ``MovingCircleCells``/``Links``.
 

src/CelerisLab/configs/CONFIG.md

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

src/CelerisLab/lbm/stepper.py

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

src/CelerisLab/simulation.py

@@ -28,7 +28,7 @@ from .cuda.context import CudaContext, context_bound
 from .cuda import compiler
 from .lbm.descriptors import FLUID
 
-# Translating circles step only through their moving tables (fixed centre in S2a).
+# Translating circles step only through their moving tables (attached at initialize).
 PENDING_BLOCKED = "a pending translation target is advanced only by Simulation.run(1)"
 MOVING_BLOCKED = ("translating circles need the moving kernels and their tables (attached "
                   "at initialize); "
@@ -122,8 +122,9 @@ class Simulation:
                 D2Q9/FP32/unshifted/DB only; set before initialization or stage
                 a new body. No runtime motion or amplitude setter.
             translation: ``body.motion.Translation`` declaring a translating
-                circle (before initialize only; docs/moving_bodies.md). Until
-                the moving kernels exist, positive-step runs are rejected.
+                circle (before initialize only; docs/moving_bodies.md).
+                Positive-step runs are rejected until ``initialize()`` has
+                attached its moving tables.
             **kwargs: passed to the geometry constructor.
 
         Returns:
@@ -241,7 +242,8 @@ class Simulation:
         self._require_usable()
         tables = self.stepper.moving
         if tables is None:
-            return np.zeros((0, 6), np.float32)
+            from .body.motion import MovingTables
+            return np.zeros((0, MovingTables.DIAG_FLOATS), np.float32)
         self.stream.synchronize()
         return tables.download_diagnostics()
 
@@ -257,6 +259,8 @@ class Simulation:
         moving = self._moving_bodies()
         if not moving:
             return
+        from .body.motion import SHARED_CLEARANCE, validate_shared_translations
+        validate_shared_translations([obj.motion for obj in moving])
         cfg = self.lbm_cfg
         if (not cfg.is_d2q9 or cfg.streaming != "double_buffer" or cfg.ddf_shifting
                 or cfg.store_precision != "FP32"):
@@ -275,6 +279,13 @@ class Simulation:
                 if other is obj:
                     continue
                 (cx, cy), r = other.center[:2], float(other.radius)
+                if (other.motion is not None and obj.motion.translation.shared
+                        and obj.motion.translation is other.motion.translation):
+                    # Whole envelopes may overlap; simultaneous local stencils may not.
+                    gap = np.hypot(obj.center[0] - cx, obj.center[1] - cy) - obj.radius - r
+                    if gap <= SHARED_CLEARANCE:
+                        raise ValueError("shared translation needs surface clearance > 4 cells")
+                    continue
                 if other.motion is not None:
                     ox0, ox1, oy0, oy1 = other.motion.band()
                 else:
@@ -360,8 +371,11 @@ class Simulation:
         Args:
             id: body_id from add_body().
             omega: rotation speed (lattice units).  ``None`` = no change.
-            vx: translation velocity x (future).  ``None`` = no change.
-            vy: translation velocity y (future).  ``None`` = no change.
+            vx: translation velocity x. For a translating circle: the
+                interval wall velocity of the target, only together with
+                ``offset``; otherwise stored in the body state (not applied
+                to stationary walls).  ``None`` = no change.
+            vy: translation velocity y; same rules as ``vx``.
             offset: target centre offset of a translating circle; with vx/vy
                 (default 0) it replaces the pending target after validation
                 against the applied offset (docs/moving_bodies.md, 4a). Only
@@ -714,6 +728,9 @@ class Simulation:
         if steps < 0 or checkpoint_interval < 0:
             raise ValueError("steps and checkpoint_interval must be nonnegative")
         moving = self._moving_bodies()
+        if steps > 0 and moving:
+            from .body.motion import validate_shared_translations
+            validate_shared_translations([obj.motion for obj in moving])
         if steps > 0 and moving and self.stepper.moving is None:
             raise RuntimeError(MOVING_BLOCKED)
         if moving and checkpoint_interval > 0:
@@ -781,6 +798,7 @@ class Simulation:
         (docs/moving_bodies.md 4a, 5). Any failure leaves the Simulation unusable."""
         tables = self.stepper.moving
         circles = [(obj.motion, obj.obj_id) for obj in self._moving_bodies()]
+        pending = any(motion.pending is not None for motion, _ in circles)
         try:
             advancing = tables.upload_records(circles, stream)
             self.stepper._step(steps, action_gpu=self.bodies.action_gpu,
@@ -793,10 +811,12 @@ class Simulation:
                 no_donor = tables.download_diagnostics()[:, tables.DIAG.index("no_donor")]
                 if np.any(no_donor > 0):
                     raise RuntimeError("translating circle refill found no donor cell")
+            if pending:  # a target equal to the applied offset is consumed as a hold step
                 for motion, _ in circles:
                     motion.commit()
-                tables.upload_records(circles, stream)  # first hold: old = new = applied, u = 0
-                stream.synchronize()
+                if advancing:
+                    tables.upload_records(circles, stream)  # first hold: old = new = applied, u = 0
+                    stream.synchronize()
                 self._attach_moving()
         except BaseException as exc:
             for motion, _ in circles:  # device geometry/state unknown after any failure

tests/integration/test_shared_translation_gpu.py

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

tests/unit/test_continuous_viv_video.py

@@ -0,0 +1,148 @@
+"""Research video contracts: continuous state and display-only sampling."""
+import importlib.util
+import json
+from pathlib import Path
+from types import SimpleNamespace
+
+import numpy as np
+import pytest
+
+ROOT=Path(__file__).resolve().parents[2]
+spec=importlib.util.spec_from_file_location('continuous_viv_video_driver',ROOT/'research/continuous_viv/run.py')
+m=importlib.util.module_from_spec(spec)
+spec.loader.exec_module(m)
+from research.continuous_viv import render as r
+
+
+def test_physical_zero_and_independent_steady_control():
+    model=SimpleNamespace(predict=lambda *a,**k:pytest.fail('steady/OFF must not invoke PPO'))
+    for scene in ('steady','karman'):
+        for phase in ('off_before','off_after'):
+            target,action=m.command(scene,phase,None,model)
+            np.testing.assert_array_equal(target,[0,0,0])
+            np.testing.assert_array_equal(m.base.wall_target(action,2),target)
+    target,action=m.command('steady','on',None,model)
+    np.testing.assert_array_equal(target,np.array([0,-5.1,5.1])*.005)
+    assert np.isnan(action).all()
+
+
+@pytest.mark.parametrize('bad',[[2,0,0],[np.nan,0,0],[0,0]])
+def test_invalid_policy_output(bad):
+    model=SimpleNamespace(predict=lambda *a,**k:(bad,None))
+    with pytest.raises(RuntimeError,match='invalid PPO'):
+        m.command('karman','on',np.zeros(12),model)
+
+
+def test_continuous_switches_previous_observation_and_native_rhs(tmp_path):
+    from CelerisLab.body.motion import MovingTables
+    calls,forces,seen,frames=[],[],[],[]
+    class Structure:
+        y=v=0.;mass=2.;stiffness=3.;frequency=.1
+        def candidate(self,f):
+            forces.append(f)
+            return self.y+.001,self.v+.002
+    class Sim:
+        stepper=SimpleNamespace(step_count=100)
+        n=0
+        def set_body(self,p,**kw):calls.append((p,kw))
+        def run(self,n,**kw):
+            assert n==1 and kw==dict(zero_obs=True,sync_obs=True)
+            self.n+=1;self.stepper.step_count+=1
+        def moving_diagnostics(self):
+            d=np.zeros((3,len(MovingTables.DIAG)),np.float32)
+            d[:,MovingTables.DIAG.index('cme_y')]=self.n
+            d[:,MovingTables.DIAG.index('wen_y')]=10*self.n
+            d[:,MovingTables.DIAG.index('cover_jy')]=1000
+            return d
+        def read_force(self,p):return [0.,float(self.n),0.]
+        def read_sensor(self,p):return [.001,0.,0.]
+    model=SimpleNamespace(predict=lambda obs,**kw:(seen.append(obs.copy()) or np.ones(3)*.1,None))
+    norm=dict(force_norm_fact=100.,sens_deviation=np.zeros(6),sens_norm_fact=np.ones(6))
+    sim,structure=Sim(),Structure()
+    def render(sim,step,phase,eff,rows):
+        if step%3==0:frames.append((step,sim.stepper.step_count,eff.copy()))
+    m.evolve(sim,[1,2,3],[4,5,6],structure,7.,scene='karman',model=model,norm=norm,
+             intervals=2,si=2,out=tmp_path,render=render,start=m.time.monotonic(),wall_cap=60)
+    assert sim.n==12 and sim.stepper.step_count==112
+    assert forces==[7.,*list(33.*np.arange(1,12))]  # no reset, SI averaging, or node impulse in RHS
+    rows=np.load(tmp_path/'steps.npy')
+    np.testing.assert_allclose(rows[:,1],.001*np.arange(1,13))
+    np.testing.assert_allclose(rows[:,2],.002*np.arange(1,13))
+    np.testing.assert_array_equal(rows[:4,15:],0)
+    np.testing.assert_array_equal(rows[8:,15:],0)
+    np.testing.assert_array_equal(rows[4,15:],rows[6,15:])  # per-SI EMA reset, not per-frame
+    signals=json.loads((tmp_path/'signals.json').read_text())
+    np.testing.assert_array_equal(seen[0],np.asarray(signals[1]['observation'],np.float32))
+    np.testing.assert_array_equal(seen[1],np.asarray(signals[2]['observation'],np.float32))
+    assert len(seen)==2 and seen[0].any()
+    report=json.loads((tmp_path/'report.json').read_text())
+    assert [x['absolute_clock'] for x in report['switches']]==[100,104,108]
+    assert [x['previous_rhs'] for x in report['switches']]==[7.,132.,264.]
+    np.testing.assert_allclose([x['y'] for x in report['switches']],[0,.004,.008])
+    assert [x[0] for x in frames]==[3,6,9,12]
+    assert all(x[1]==100+x[0] for x in frames)
+    assert len(calls)==36 and all(c[1]['vy']==pytest.approx(.001) for c in calls)
+
+
+def test_gear_schedule_and_complete_phase_endpoints():
+    n=1280000
+    schedule=r.video_schedule(n)
+    times=list(schedule)
+    assert times==sorted(set(times)) and times[-1]==3*n
+    assert all(p*n in schedule for p in (1,2,3))
+    assert set(schedule.values())=={1,2,5,10}
+    assert all(speed==1 for step,speed in schedule.items() if step%n>n//2)
+    part=r.frame_schedule(n//2,800,True)
+    assert [s for _,s in part[:15]]==[1]*15
+    assert [s for _,s in part[15:30]]==[2]*15
+    assert [s for _,s in part[-15:]]==[1]*15
+    assert part[-1][0]==n//2
+
+
+def test_smoke_schedule_covers_all_switches():
+    schedule=r.video_schedule(32,True)
+    assert list(schedule)==list(range(8,97,8))
+    assert set(schedule.values())=={1}
+
+
+def test_vorticity_sign_scale_and_solid_stencil_mask():
+    y,x=np.mgrid[:9,:11]
+    ux=-.0001*y;uy=.0001*x
+    flags=np.ones(x.shape,np.uint16)
+    omega=r.vorticity(ux,uy,flags)
+    np.testing.assert_allclose(omega[1:-1,1:-1],.0002*40/.005)
+    assert np.isnan(omega[0]).all()
+    flags[4,5]=2
+    omega=r.vorticity(ux,uy,flags)
+    for yy,xx in ((4,5),(3,5),(5,5),(4,4),(4,6)):
+        assert np.isnan(omega[yy,xx])
+    assert np.isfinite(omega[3,4])
+
+
+def test_fixed_inset_coordinate_and_no_motion_gain():
+    a=np.array(r.pixel(30,0,r.INSET,r.INSET_LIMITS))
+    b=np.array(r.pixel(30,1,r.INSET,r.INSET_LIMITS))
+    c=np.array(r.pixel(31,0,r.INSET,r.INSET_LIMITS))
+    assert b[0]==a[0] and a[1]-b[1]==56
+    assert c[0]-a[0]==56
+    assert r.pixel(r.INSET_LIMITS[0],r.INSET_LIMITS[3],r.INSET,r.INSET_LIMITS)==r.INSET[:2]
+
+
+def test_config_real_loader_and_structure_unchanged(tmp_path):
+    from CelerisLab.config import load_lbm_config
+    p=tmp_path/'config.json';p.write_text(json.dumps(m.base.config(40)))
+    cfg=load_lbm_config(str(p));cfg.validate()
+    assert (cfg.nx,cfg.ny,cfg.velocity,cfg.viscosity)==(2559,1022,.005,.004)
+    s=m.base.Structure(4.,2.5e-5,0.,40.)
+    assert s.mass==pytest.approx(15079.644737231007)
+    assert s.stiffness==pytest.approx(.00037207532016359776)
+    assert s.damping==0 and s.frequency*40/.005==.2
+
+
+def test_scene_timing_keeps_actual_legacy_intervals_and_same_physical_duration():
+    assert m.timing('steady',False)==(320,4000,2048000)
+    assert m.timing('karman',False)==(400,3200,2048000)
+    for scene in ('steady','karman'):
+        n,si,warm=m.timing(scene,False)
+        assert n*si*.005/40==160
+        assert m.timing(scene,True)==(1,32,32)

tests/unit/test_duetsch_protocol.py

@@ -11,6 +11,20 @@ _spec = importlib.util.spec_from_file_location(
 p = importlib.util.module_from_spec(_spec)
 _spec.loader.exec_module(p)
 
+def _load(rel, name):
+    """Load a research driver module under a unique name (no bare ``import run``)."""
+    import sys
+    path = Path(__file__).resolve().parents[2] / rel
+    sys.path.insert(0, str(path.parent))  # for the driver's own sibling imports
+    try:
+        spec = importlib.util.spec_from_file_location(name, path)
+        mod = importlib.util.module_from_spec(spec)
+        spec.loader.exec_module(mod)
+    finally:
+        sys.path.remove(str(path.parent))
+    return mod
+
+
 T = 2520
 
 
@@ -68,6 +82,35 @@ def test_interval_bases_and_morison_fit_recover_known_coefficients():
     assert p.morison_fit(f, b[:, 0], b[:, 1]) == pytest.approx((2.09, 1.45), abs=1e-9)
 
 
+def test_quadrature_rule_reused_across_steps_and_periods(monkeypatch):
+    original = np.polynomial.legendre.leggauss
+    orders = []
+
+    def counted(points):
+        orders.append(points)
+        return original(points)
+
+    p._quadrature_rule.cache_clear()
+    monkeypatch.setattr(np.polynomial.legendre, "leggauss", counted)
+    try:
+        for points in (16, 8, 16):
+            for period in (T, 2 * T):
+                for n in (0, 123, period + 123):
+                    p.interval_bases(n, period, points)
+        assert orders == [16, 8]
+    finally:
+        p._quadrature_rule.cache_clear()
+
+
+def test_interval_sine_basis_matches_exact_integral():
+    for period in (T, 2 * T):
+        half_width = math.pi / period
+        for n in (0, period // 4, period // 2, period - 1, 7 * period + 123):
+            midpoint = 2 * math.pi * ((n % period) + .5) / period - math.pi / 2
+            expected = math.sin(midpoint) * math.sin(half_width) / half_width
+            assert p.interval_bases(n, period)[1] == pytest.approx(expected, rel=0, abs=2e-15)
+
+
 def test_geometry_repeats_exactly_across_cycles():
     A = p.amplitude(20)
     for m in (0, 1, T // 4, T // 2 - 1, 7 * T // 12):
@@ -93,12 +136,7 @@ def test_exact_point_inside_body_is_missing_even_with_fluid_corners():
 
 def test_generated_config_passes_the_real_loader(tmp_path):
     import json
-    import sys
-    sys.path.insert(0, str(Path(__file__).resolve().parents[2] / "research/duetsch_re100_kc5"))
-    try:
-        import run
-    finally:
-        sys.path.pop(0)
+    run = _load("research/duetsch_re100_kc5/run.py", "duetsch_run")
     from CelerisLab.config import load_lbm_config
     cfg_dict, u_max, n = run.config(20, 10.0, 2520)
     path = tmp_path / "config_lbm.json"
@@ -111,13 +149,7 @@ def test_generated_config_passes_the_real_loader(tmp_path):
 
 
 def _analyze():
-    import sys
-    sys.path.insert(0, str(Path(__file__).resolve().parents[2] / "research/duetsch_re100_kc5"))
-    try:
-        import analyze
-    finally:
-        sys.path.pop(0)
-    return analyze
+    return _load("research/duetsch_re100_kc5/analyze.py", "duetsch_analyze")
 
 
 def test_hf_energy_is_parseval_weighted():

tests/unit/test_fig2_analysis.py

@@ -0,0 +1,36 @@
+"""Independent analytic data checks for offline estimands, no CFD context."""
+import unittest
+import numpy as np
+from research.fig2_collection.analyze import frequency, native_summary, compare_selected
+
+
+class AnalysisTests(unittest.TestCase):
+    def test_frequency_units_and_no_dc(self):
+        dt=.005/40;n=640000
+        t=np.arange(n)*dt
+        result=frequency(4+np.sin(2*np.pi*.2*t),dt)
+        self.assertAlmostEqual(result['f_D_U'],.2)
+        self.assertAlmostEqual(result['bin_width_D_U'],1/80)
+        self.assertIsNone(frequency(np.ones(20),dt))
+
+    def test_motion_and_native_rhs_are_distinct_from_node_impulse(self):
+        a=np.zeros((10,18));a[:,0]=np.arange(1,11);a[:,1]=40*np.tile([.5,-.5],5)
+        a[:,3]=.001;a[:,6]=.002;a[:,9]=999
+        result=native_summary(a,'taylor')
+        self.assertEqual(result['motion']['centered_rms'],.5)
+        self.assertEqual(result['motion']['half_range'],.5)
+        self.assertAlmostEqual(result['structural_RHS_Cy_RMS'],6.)
+        self.assertEqual(result['window_D_U'],[0.,10*.005/40])
+
+    def test_selected_vector_error_uses_common_roi_and_actual_flags(self):
+        shape=(1,712,2162)
+        target=dict(ux=np.zeros(shape,np.float32),uy=np.zeros(shape,np.float32),flags=np.ones(shape,np.uint16))
+        current={k:v.copy() for k,v in target.items()}
+        current['ux'][:,311:711,1360:2161]=.003
+        current['uy'][:,311:711,1360:2161]=.004
+        self.assertAlmostEqual(compare_selected(target,current)['aggregate_E'],1.,places=6)
+        current['flags'][0,400,1500]=2
+        with self.assertRaises(ValueError):compare_selected(target,current)
+
+
+if __name__=='__main__':unittest.main()

tests/unit/test_fig2_collection.py

@@ -0,0 +1,260 @@
+"""CPU contracts for Fig. 2 acquisition; no CUDA or simulation import required."""
+import importlib.util
+import json
+from pathlib import Path
+import sys
+import tempfile
+import types
+import unittest
+from unittest.mock import patch
+
+import numpy as np
+
+ROOT = Path(__file__).resolve().parents[2]
+SPEC = importlib.util.spec_from_file_location('fig2_collection_tested', ROOT/'research/fig2_collection/run.py')
+driver = importlib.util.module_from_spec(SPEC)
+SPEC.loader.exec_module(driver)
+METRIC_SPEC = importlib.util.spec_from_file_location('fig2_collection_metrics', ROOT/'research/fig2_materials/metrics.py')
+metrics = importlib.util.module_from_spec(METRIC_SPEC)
+METRIC_SPEC.loader.exec_module(metrics)
+
+
+class CollectionContracts(unittest.TestCase):
+    def test_frozen_clocks_and_windows(self):
+        for scene, release, count, first in [('karman', 4576000, 400, 240),
+                                             ('steady', 2048000, 400, 240),
+                                             ('taylor', 1024000, 150, 0)]:
+            plan = driver.schedule(scene)
+            self.assertEqual((plan['release_clock'], plan['intervals'], plan['field_start']),
+                             (release, count, first))
+            self.assertEqual(plan['si'], 3200)
+        self.assertEqual((400-240)*3200*driver.U/driver.D, 64.)
+        self.assertEqual(400*3200*driver.U/driver.D, 160.)
+        self.assertEqual(150*3200*driver.U/driver.D, 60.)
+
+    def test_smoke_has_no_physical_timing(self):
+        for scene in ('karman', 'steady', 'taylor'):
+            plan = driver.schedule(scene, True)
+            self.assertEqual((plan['si'], plan['intervals'], plan['field_start']), (32, 3, 0))
+            self.assertEqual(plan['release_clock'], 64 if scene == 'taylor' else 32)
+
+    def test_zero_is_exact_physical_zero_despite_bias(self):
+        for scene in ('karman', 'steady', 'taylor'):
+            np.testing.assert_array_equal(driver.wall_surface(scene, 'zero'), np.zeros(3))
+            np.testing.assert_array_equal(driver.wall_surface(scene, 'target'), np.zeros(3))
+        np.testing.assert_allclose(driver.wall_surface('karman', 'controlled', [.5, -.5, 1]),
+                                   [.02, -.04, .06])
+        np.testing.assert_allclose(driver.wall_surface('taylor', 'controlled', [.5, -.5, 1]),
+                                   [.01, -.03, .04])
+        np.testing.assert_allclose(driver.wall_surface('steady', 'controlled'), [0, -.0255, .0255])
+        with self.assertRaises(ValueError):
+            driver.wall_surface('karman', 'controlled', [0, 0, np.nan])
+
+    def test_structure_hard_stops(self):
+        driver.check_structure(40., 0., 40.)
+        for y, v, old in [(40.01, 0, 40), (0, .026, 0), (.026, 0, 0), (np.nan, 0, 0)]:
+            with self.subTest(y=y, v=v), self.assertRaises(RuntimeError):
+                driver.check_structure(y, v, old)
+
+    def test_fluid_nonfinite_fails_nonfluid_only_zero_filled(self):
+        flags = np.array([[1, 257, 2]], np.uint16)
+        fields = dict(rho=np.array([[1., 1., np.nan]]), ux=np.array([[0., .1, np.nan]]),
+                      uy=np.array([[0., -.1, np.inf]]))
+        clean, retained, fluid = driver.clean_field(fields, flags)
+        np.testing.assert_array_equal(fluid, [[True, True, False]])
+        np.testing.assert_array_equal(retained, flags)
+        self.assertTrue(np.isfinite(clean).all())
+        self.assertEqual(clean[0, 0, 0], 0.)  # zero flow remains fluid
+        fields['ux'][0, 1] = np.nan
+        with self.assertRaisesRegex(RuntimeError, 'FLUID'):
+            driver.clean_field(fields, flags)
+
+    def test_memory_admission_uses_available_not_free(self):
+        with tempfile.TemporaryDirectory() as directory:
+            path = Path(directory)/'meminfo'
+            path.write_text('MemFree: 1 kB\nMemAvailable: 16777216 kB\n')
+            self.assertEqual(driver.available_memory(path), driver.MIN_AVAILABLE)
+        with patch.object(driver, 'available_memory', return_value=driver.MIN_AVAILABLE-1):
+            with self.assertRaises(MemoryError):
+                driver.admit_memory()
+        self.assertIsNone(driver.admit_memory(True))
+
+    def test_roi_bounds_match_full_downstream_support(self):
+        roi = metrics.roi_slices((1022, 2559), D=40, xs=1600, y0=510.5,
+                                downstream_solid_extent=1272.)
+        self.assertEqual(roi, (slice(311, 711), slice(1360, 2161)))
+
+    def test_same_clock_error_and_mean_are_distinct(self):
+        shape = (2, 2, 2, 3)
+        target = np.zeros(shape, np.float32)
+        changed = target.copy()
+        changed[0, 0], changed[1, 0] = driver.U, -driver.U
+        role = lambda fields: dict(roi=fields, mask=np.ones((2, 3), bool), times=np.array([3200, 6400]))
+        roles = dict(target=role(target), zero=role(changed), controlled=role(target))
+        result, mask = driver.score_roles(roles, metrics)
+        self.assertAlmostEqual(result['zero']['mean_instantaneous_E'], 1., places=7)
+        self.assertEqual(result['zero']['E_of_fixed_window_mean'], 0.)
+        roles['zero']['times'] = np.array([6400, 9600])
+        with self.assertRaisesRegex(ValueError, 'clocks'):
+            driver.score_roles(roles, metrics)
+        roles['zero']['times'] = roles['target']['times']
+        roles['zero']['mask'][0, 0] = False
+        with self.assertRaisesRegex(ValueError, 'entirely common fluid'):
+            driver.score_roles(roles, metrics)
+
+    def test_native_archive_verified_raw_retained(self):
+        with tempfile.TemporaryDirectory() as directory:
+            path = Path(directory)/'native.npy'
+            rows = np.lib.format.open_memmap(path, mode='w+', dtype=np.float64,
+                                            shape=(9, len(driver.NATIVE_COLUMNS)))
+            rows[:] = np.arange(rows.size).reshape(rows.shape)
+            report = driver.compress_native(path, rows, 7)
+            self.assertTrue(report['verified'])
+            self.assertTrue(path.exists())
+            with np.load(path.with_suffix('.npz')) as data:
+                np.testing.assert_array_equal(data['rows'], rows[:7])
+
+    def test_readonly_import_reuses_actual_offsets_and_functions(self):
+        with tempfile.TemporaryDirectory() as directory:
+            root = Path(directory)
+            (root/'legacy_test/core').mkdir(parents=True)
+            # Package __init__ must never be executed by the isolated source loader.
+            for folder in (root, root/'legacy_test', root/'legacy_test/core'):
+                (folder/'__init__.py').write_text("raise RuntimeError('application import forbidden')\n")
+            (root/'legacy_test/core/dtw_metrics.py').write_text('calc_dtw_sim = None\ncalc_lag = None\n')
+            (root/'acquisition.py').write_text(
+                'from .legacy_test.core.dtw_metrics import calc_dtw_sim, calc_lag\n'
+                'VORTEX_EVENT_OFFSETS = (-20,-15,-10,-5,0,5,10,15)\n')
+            (root/'legacy_test/acquire.py').write_text(
+                "raise RuntimeError('heavy imports forbidden')\n"
+                'def _vortex_event_channel(case_name):\n    return 3, "center_sensor_uy"\n'
+                'def _select_vortex_event(sensors, case_name):\n'
+                '    return dict(offsets=VORTEX_EVENT_OFFSETS, channel=_vortex_event_channel(case_name))\n')
+            helper = driver.load_helpers(root/'acquisition.py')
+            result = helper.select_vortex_event(None, 'vortex_taylor')
+            self.assertEqual(result['offsets'], (-20,-15,-10,-5,0,5,10,15))
+            self.assertEqual(result['channel'], (3, 'center_sensor_uy'))
+            self.assertFalse(list(root.rglob('*.pyc')))
+
+    def test_geometry_and_shared_translation_envelope(self):
+        calls, translations = [], []
+        def translation(*args, **kwargs):
+            item = types.SimpleNamespace(args=args, kwargs=kwargs)
+            translations.append(item)
+            return item
+        def add_body(*args, **kwargs):
+            calls.append((args, kwargs))
+            return len(calls)-1
+        motion = types.ModuleType('CelerisLab.body.motion')
+        motion.Translation = translation
+        with patch.dict(sys.modules, {'CelerisLab.body.motion': motion}):
+            probes, bodies = driver.add_scene(types.SimpleNamespace(add_body=add_body), 'karman', moving=True)
+        self.assertEqual(calls[0][0], ('circle', (400.,510.5),40.))
+        self.assertEqual(probes, [1,2,3])
+        self.assertEqual(bodies, [4,5,6])
+        self.assertEqual(translations[0].args, ((0.,0.,-40.,40.), .025))
+        self.assertTrue(translations[0].kwargs['shared'])
+        self.assertTrue(all(calls[i][1]['translation'] is translations[0] for i in bodies))
+        calls.clear()
+        with patch.dict(sys.modules, {'CelerisLab.body.motion': motion}):
+            driver.add_scene(types.SimpleNamespace(add_body=add_body), 'steady', pinball=False)
+        self.assertEqual(len(calls), 3)
+        self.assertTrue(all(call[0][0] == 'sensor' for call in calls))
+
+    def test_taylor_operator_called_once_inside_context_at_response_zero(self):
+        events = []
+        class Context:
+            def __enter__(self):
+                events.append('enter')
+            def __exit__(self, *args):
+                events.append('exit')
+        initializer = types.ModuleType('CelerisLab.lbm.initializers')
+        initializer.add_vortex = lambda *args: events.append(args)
+        field = object()
+        sim = types.SimpleNamespace(field=field, stepper=types.SimpleNamespace(step_count=1024000),
+                                    ctx=types.SimpleNamespace(activate=Context),
+                                    stream=types.SimpleNamespace(synchronize=lambda: events.append('sync')))
+        with patch.dict(sys.modules, {'CelerisLab.lbm.initializers': initializer}):
+            receipt = driver.inject_taylor(sim)
+        self.assertEqual(events, ['enter', 'sync', (field, (600.,510.5),80.,.000075,'taylor'), 'sync', 'exit'])
+        self.assertEqual((receipt['response_clock'], receipt['absolute_clock'], receipt['count']),
+                         (0,1024000,1))
+
+
+class SelectionContracts(unittest.TestCase):
+    def test_taylor_release_uses_committed_ids_not_staged_minus_one(self):
+        reads=[]
+        class Sim:
+            def __init__(self,*args):
+                self.bodies=types.SimpleNamespace(count=0)
+                self.stepper=types.SimpleNamespace(step_count=0)
+                self.initialized=False;self.pending=0
+            def __enter__(self):return self
+            def __exit__(self,*args):pass
+            def add_body(self,*args,**kwargs):
+                if self.initialized:self.pending+=1;return -1
+                index=self.bodies.count;self.bodies.count+=1;return index
+            def initialize(self):self.initialized=True
+            def sync_bodies(self):self.bodies.count+=self.pending;self.pending=0
+            def run(self,n,**kwargs):self.stepper.step_count+=n
+            def read_force(self,body):
+                if body<0:raise IndexError('staged ID')
+                reads.append(body);return [0.,float(body),0.]
+        pkg=types.ModuleType('CelerisLab');pkg.Simulation=Sim
+        motion=types.ModuleType('CelerisLab.body.motion');motion.Translation=None
+        base=types.SimpleNamespace(warmup=lambda sim,n,*args:sim.run(n),
+            capture_seed=lambda sim:(np.zeros(18,np.float32),np.ones(18,np.float32),
+                                      np.ones(2,np.uint16),sim.stepper.step_count))
+        with tempfile.TemporaryDirectory() as directory,patch.dict(sys.modules,{'CelerisLab':pkg,'CelerisLab.body.motion':motion}):
+            seed,force,_=driver.make_release('taylor',driver.schedule('taylor',True),'unused',base,
+                                             Path(directory),0,60,None)
+        self.assertEqual(reads,[3,4,5])
+        self.assertEqual(force,12.)
+        self.assertEqual(seed[3],64)
+
+    def test_recorded_clock_must_equal_observed_solver(self):
+        driver.check_clock(types.SimpleNamespace(step_count=4579200),4576000,3200)
+        with self.assertRaisesRegex(RuntimeError,'observed solver clock'):
+            driver.check_clock(types.SimpleNamespace(step_count=4579201),4576000,3200)
+
+    @staticmethod
+    def mean(ux, uy, indices):
+        return dict(mean_ux=ux[indices].mean(axis=0), mean_uy=uy[indices].mean(axis=0))
+
+    def test_karman_missing_cycles_preserves_fixed_mean_and_last(self):
+        helper = types.SimpleNamespace(accumulate_mean_fields=self.mean)
+        def no_crossing(*args, **kwargs):
+            raise ValueError('at least two accepted rising crossings are required')
+        helper.assign_periodic_phase = no_crossing
+        cache = np.arange(24, dtype=np.float32).reshape(3, 2, 2, 2)
+        arrays, metadata, event = driver.select_fields('karman', dict(smoke=False), helper,
+                                                      np.array([3200, 6400, 9600]), np.zeros((3, 6)), cache)
+        self.assertEqual(metadata['periodic_phase_status'], 'insufficient_crossings')
+        self.assertNotIn('cycle_mean_ux', arrays)
+        np.testing.assert_array_equal(arrays['ux'][0], cache[-1, 0])
+        np.testing.assert_allclose(arrays['fixed_window_mean_ux'], cache[:, 0].mean(axis=0))
+
+    def test_taylor_common_controlled_event_not_reselected(self):
+        def forbidden(*args):
+            raise AssertionError('other roles must not select their own event')
+        helper = types.SimpleNamespace(accumulate_mean_fields=self.mean, select_vortex_event=forbidden)
+        cache = np.zeros((150, 2, 2, 2), np.float32)
+        event = dict(field_indices=np.array([20,25,30,35,40,45,50,55]))
+        arrays, meta, returned = driver.select_fields('taylor', dict(smoke=False), helper,
+                                                     np.arange(1,151)*3200, np.zeros((150,6)), cache, event)
+        self.assertIs(returned, event)
+        np.testing.assert_array_equal(arrays['selected_index'], event['field_indices'])
+        self.assertEqual(len(arrays['ux']), 8)
+
+    def test_smoke_selection_has_no_event_or_phase_claim(self):
+        helper = types.SimpleNamespace(accumulate_mean_fields=self.mean)
+        for scene in ('karman', 'steady', 'taylor'):
+            arrays, metadata, event = driver.select_fields(scene, dict(smoke=True), helper,
+                                                          np.array([32,64,96]), np.zeros((3,6)),
+                                                          np.zeros((3,2,2,2), np.float32))
+            self.assertEqual(metadata['selection'], 'smoke_last_snapshot_not_physical')
+            self.assertIsNone(event)
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/unit/test_fig2_metrics.py

@@ -0,0 +1,49 @@
+"""CPU-only analytic checks, runnable with stdlib unittest (or pytest)."""
+import unittest
+import numpy as np
+from research.fig2_materials.metrics import (
+    common_fluid_roi_mask, field_error, motion_summary, roi_slices,
+)
+
+
+class MetricsTests(unittest.TestCase):
+    def test_u0_normalization_not_target_energy(self):
+        zero = np.zeros((2, 2)); support = np.ones((2, 2), bool)
+        self.assertEqual(field_error(zero + 3, zero + 4, zero, zero, 2, support), 2.5)
+        self.assertEqual(field_error(zero + 13, zero + 4, zero + 10, zero, 2, support), 2.5)
+        for u0 in (0, -1, np.nan):
+            with self.assertRaises(ValueError):
+                field_error(zero, zero, zero, zero, u0, support)
+
+    def test_common_mask_roi_finite_and_fluid(self):
+        a = np.zeros((3, 4)); b = a.copy(); b[1, 2] = np.nan
+        m = np.ones(a.shape, bool); n = m.copy(); n[1, 1] = False
+        support = common_fluid_roi_mask([m, n], [a, b], (slice(1, 3), slice(1, 4)))
+        self.assertEqual(int(support.sum()), 4)
+        self.assertFalse(support[1, 1]); self.assertFalse(support[1, 2])
+        self.assertTrue(support[2, 1])  # zero velocity is fluid
+        self.assertEqual(field_error(a, a, b, a, 1, support), 0)
+        with self.assertRaises(ValueError):
+            common_fluid_roi_mask([m & False], [a], (slice(None), slice(None)))
+        with self.assertRaises(ValueError):
+            field_error(a, a, b, a, 1, m)  # fixed support must be finite
+
+    def test_mean_field_differs_from_mean_instantaneous_error(self):
+        snapshots = np.array([[[2.]], [[-2.]]]); z = np.zeros((1, 1)); mask = np.ones((1, 1), bool)
+        self.assertEqual(field_error(snapshots.mean(0), z, z, z, 1, mask), 0)
+        self.assertEqual(np.mean([field_error(s, z, z, z, 1, mask) for s in snapshots]), 2)
+
+    def test_actual_geometry_and_motion_bias(self):
+        roi = roi_slices((1022, 2400), D=40, xs=1600, y0=510.5, downstream_solid_extent=1350)
+        self.assertEqual(roi, (slice(311, 711), slice(1360, 2161)))
+        m = motion_summary([8, 10, 12])
+        self.assertEqual(m['mean'], 10)
+        self.assertAlmostEqual(m['centered_rms'], np.sqrt(8 / 3))
+        self.assertEqual(m['half_range'], 2)
+        for a in ([], [np.nan]):
+            with self.assertRaises(ValueError):
+                motion_summary(a)
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/unit/test_karman_viv.py

@@ -0,0 +1,183 @@
+"""Research interface contracts; no CUDA context or physical validation."""
+import importlib.util
+from pathlib import Path
+from types import SimpleNamespace
+from unittest.mock import patch
+
+import numpy as np
+import pytest
+
+ROOT = Path(__file__).resolve().parents[2]
+spec = importlib.util.spec_from_file_location('karman_viv_research_driver', ROOT/'research/karman_viv/run.py')
+m = importlib.util.module_from_spec(spec)
+spec.loader.exec_module(m)
+
+
+def test_observation_unit_mapping_and_order():
+    norm = dict(force_norm_fact=2.,sens_deviation=np.arange(6)/10.,sens_norm_fact=np.ones(6)*3)
+    s=np.arange(6)*.001; f=np.arange(6)*.01
+    raw,obs,nclip=m.observation(s,f,2,norm)
+    np.testing.assert_array_equal(raw,np.asarray([*(s*156),*(f*2)],np.float32))
+    expected=np.hstack((raw[6:]/2,(raw[:6]-np.asarray(norm['sens_deviation'],np.float32))/3))
+    np.testing.assert_array_equal(obs,expected.astype(np.float32))
+    assert nclip==0
+
+
+def test_same_physical_signal_d20_d40():
+    norm=dict(force_norm_fact=1.,sens_deviation=np.zeros(6),sens_norm_fact=np.ones(6))
+    a=m.observation(np.ones(6)*.001,np.arange(6)*.01,1,norm)
+    b=m.observation(np.ones(6)*.0005,np.arange(6)*.005,2,norm)
+    np.testing.assert_array_equal(a[0],b[0]); np.testing.assert_array_equal(a[1],b[1])
+
+
+def test_zero_action_and_ccw_scale():
+    np.testing.assert_array_equal(m.wall_target(-m.BIAS/8,2),np.zeros(3))
+    np.testing.assert_allclose(m.wall_target(np.ones(3),2),[.04,.02,.06])
+    assert m.wall_target(np.ones(3),2)[0]/20 == .002
+
+
+def test_ema_equal_physical_time_and_interval_reset():
+    target=np.array([.01,-.02,.03]); coarse=np.zeros(3); fine=np.zeros(3)
+    for _ in range(13): coarse=m.ema(coarse,target,1)
+    for _ in range(52): fine=m.ema(fine,target,2)
+    np.testing.assert_allclose(fine,coarse,rtol=2e-15,atol=1e-17)
+    np.testing.assert_allclose(m.ema(np.zeros(3),target,2),(1-.9**.25)*target)
+
+
+def test_config_real_parser(tmp_path):
+    import json
+    from CelerisLab.config import load_lbm_config
+    path=tmp_path/'config.json';path.write_text(json.dumps(m.config(40)))
+    cfg=load_lbm_config(str(path));cfg.validate()
+    assert (cfg.nx,cfg.ny,cfg.velocity,cfg.viscosity)==(2559,1022,.005,.004)
+
+
+def test_structure_physical_scaling_and_no_damping():
+    coarse=m.Structure(2,1e-4,0,20);fine=m.Structure(2,2.5e-5,0,40)
+    assert fine.mass==4*coarse.mass
+    assert fine.stiffness==coarse.stiffness/4
+    assert fine.damping==0
+    assert fine.candidate(0)==(0.,0.)
+
+
+def test_observation_clipping_and_nonfinite():
+    norm=dict(force_norm_fact=1.,sens_deviation=np.zeros(6),sens_norm_fact=np.ones(6))
+    _,obs,n=m.observation(np.ones(6),np.ones(6),2,norm)
+    assert n==12 and np.array_equal(obs,np.ones(12,np.float32))
+    with pytest.raises(RuntimeError,match='nonfinite'):
+        m.observation(np.full(6,np.nan),np.zeros(6),2,norm)
+
+
+@pytest.mark.parametrize('mass_ratio', [2.,4.])
+def test_arm_step_order_and_native_rhs(tmp_path, mass_ratio):
+    import CelerisLab
+    from CelerisLab.body.motion import MovingTables
+    from contextlib import nullcontext
+    log=[];candidates=[];structure_args=[]
+    class Structure:
+        y=v=0.;mass=2.;stiffness=3.;frequency=.1
+        def __init__(self,*args): structure_args.append(args)
+        def candidate(self,f):
+            candidates.append(f)
+            return self.y+.001,self.v+.002
+    class Sim:
+        def __init__(self,*args): self.n=0
+        def __enter__(self): return self
+        def __exit__(self,*args): pass
+        def initialize(self): pass
+        def set_body(self,p,**kw): log.append(('set',p,kw))
+        def run(self,n,**kw): self.n+=1;log.append(('run',n,kw))
+        def moving_diagnostics(self):
+            a=np.zeros((3,len(MovingTables.DIAG)),np.float32)
+            a[:,MovingTables.DIAG.index('cme_y')]=self.n
+            a[:,MovingTables.DIAG.index('wen_y')]=10*self.n
+            a[:,MovingTables.DIAG.index('cover_jy')]=1000
+            return a
+        def read_force(self,p): return [0.,float(self.n),0.]
+        def read_sensor(self,p): return [.001,0.,0.]
+    seen=[]
+    model=SimpleNamespace(predict=lambda obs,**kw:(seen.append(obs.copy()) or np.ones(3)*.1,None))
+    norm=dict(force_norm_fact=1.,sens_deviation=np.zeros(6),sens_norm_fact=np.ones(6))
+    with patch.object(CelerisLab,'Simulation',Sim),patch.object(m,'Structure',Structure), \
+         patch.object(m,'add_scene',return_value=([1,2,3],[4,5,6])),patch.object(m,'install_seed',return_value={}), \
+         patch.object(m,'check_time'):
+        m.run_arm('controlled','cfg',2,2,2,None,7.,model,norm,tmp_path,0,1,mass_ratio=mass_ratio)
+    assert structure_args == [(mass_ratio,2.5e-5,0.,40.)]
+    assert candidates==[7.,33.,66.,99.]  # previous native force, NOT SI mean or node impulse
+    assert not seen[0].any() and len(seen)==2
+    assert [e[0] for e in log]==['set','set','set','run']*4
+    for e in log:
+        if e[0]=='set': assert e[2]['vy']==pytest.approx(.001) and e[2]['omega']!=0
+        else: assert e[2]==dict(zero_obs=True,sync_obs=True)
+    rows=np.load(tmp_path/'controlled/steps.npy')
+    np.testing.assert_allclose(rows[:,1],[.001,.002,.003,.004])
+    assert rows.shape==(4,18)
+    np.testing.assert_array_equal(rows[0,15:],rows[2,15:])  # EMA resets at SI boundary
+    signals=np.load(tmp_path/'controlled/signals.npz')
+    np.testing.assert_array_equal(signals['forces'][:,1],[3.,7.])  # interval CME *2, not CME+Wen
+
+
+def test_mass4_changes_inertia_and_stiffness_not_frequency_or_damping():
+    a=m.Structure(2,2.5e-5,0,40);b=m.Structure(4,2.5e-5,0,40)
+    assert b.mass==2*a.mass and b.stiffness==2*a.stiffness
+    assert b.frequency==a.frequency and b.damping==a.damping==0
+    np.testing.assert_allclose(b.candidate(.001),np.asarray(a.candidate(.001))/2,rtol=0,atol=0)
+
+
+def make_release(path):
+    import json
+    (path/'zero').mkdir()
+    case=dict(source=m.RELEASE_SOURCE,model_sha256=m.MODEL_SHA,norm_sha256=m.NORM_SHA,
+              smoke=False,D=40,U=.005,nu=.004,warmup_each=2048000,target_intervals=150,SI=3200)
+    (path/'case.json').write_text(json.dumps(case))
+    (path/'config.json').write_text(json.dumps(m.config(40)))
+    seed=(np.arange(18,dtype=np.float32),np.arange(18,dtype=np.float32)+1,np.zeros(2,np.int32))
+    np.savez(path/'release.npz',current=seed[0],spare=seed[1],flags=seed[2],clock=4576000,force=.003)
+    np.save(path/'target.npy',np.zeros((150,6),np.float32))
+    transfer=dict(zip(('current','spare','flags'),map(m.array_sha,seed)),clock=4576000)
+    (path/'zero/report.json').write_text(json.dumps(dict(release_transfer=transfer)))
+    return seed
+
+
+@pytest.mark.parametrize('damage', [None,'source','config','buffer','target','clock'])
+def test_saved_release_checks_provenance_and_exact_inputs(tmp_path, damage):
+    import json
+    seed=make_release(tmp_path)
+    if damage in ('source','config'):
+        file=tmp_path/('case.json' if damage=='source' else 'config.json')
+        data=json.loads(file.read_text());data['source' if damage=='source' else 'extra']='changed'
+        file.write_text(json.dumps(data))
+    elif damage in ('buffer','clock'):
+        with np.load(tmp_path/'release.npz') as z: data={k:z[k] for k in z.files}
+        if damage=='buffer': data['current'][0]+=1
+        else: data['clock']=4576001
+        np.savez(tmp_path/'release.npz',**data)
+    elif damage=='target': np.save(tmp_path/'target.npy',np.zeros((149,6),np.float32))
+    if damage:
+        with pytest.raises(ValueError): m.saved_release(tmp_path)
+    else:
+        loaded,force=m.saved_release(tmp_path)
+        assert loaded[3]==4576000 and force==.003
+        for actual,expected in zip(loaded[:3],seed): np.testing.assert_array_equal(actual,expected)
+
+
+def test_mass4_cli_uses_saved_release_without_warmup(tmp_path, monkeypatch):
+    import json
+    original=tmp_path/'original';original.mkdir();make_release(original)
+    norm=tmp_path/'norm.json';norm.write_text('{}');model=tmp_path/'model.zip'
+    out=tmp_path/'output';calls=[]
+    monkeypatch.setattr(m.sys,'argv',['run.py','--model',str(model),'--norm',str(norm),
+                                    '--out',str(out),'--mass-ratio','4','--release-dir',str(original)])
+    monkeypatch.setattr(m,'sha',lambda path:m.MODEL_SHA if path==model else m.NORM_SHA if path==norm else 'fixture-sha')
+    monkeypatch.setattr(m,'load_policy',lambda path:object())
+    monkeypatch.setattr(m.subprocess,'check_output',lambda *a,**k:'candidate\n')
+    monkeypatch.setattr(m,'precursor',lambda *a:pytest.fail('saved release must skip precursor'))
+    monkeypatch.setattr(m,'run_arm',lambda *a,**k:calls.append((a,k)))
+    m.main()
+    assert [a[0] for a,k in calls]==['zero','controlled']
+    assert all(a[2:5]==(2.,3200,400) and k=={'mass_ratio':4.} for a,k in calls)
+    receipt=json.loads((out/'case.json').read_text())
+    assert receipt['mass_ratio']==4 and receipt['structural_damping']==0 and receipt['warmup_executed']==0
+    assert receipt['exploratory'] and receipt['response_tU_D']==160
+    assert (out/'target.npy').read_bytes()==(original/'target.npy').read_bytes()
+    assert (out/'config.json').read_bytes()==(original/'config.json').read_bytes()

tests/unit/test_moving_bodies.py

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

tests/unit/test_shared_translation.py

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

tests/unit/test_steady_pinball.py

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

tests/unit/test_viv_logic.py

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