Use committed Taylor body IDs after topology sync

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/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>`.

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/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()

src/CelerisLab/body/motion.py

@@ -14,6 +14,7 @@ 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]
@@ -31,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)
@@ -41,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."""
 
@@ -72,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}")

src/CelerisLab/simulation.py

@@ -259,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"):
@@ -277,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:
@@ -719,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:

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_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_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) == []