Add continuous OFF ON OFF VIV video research driver

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

@@ -102,3 +102,49 @@ per-case wall cap and GPU headroom check bound resource use. Preserve partial ro
 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/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/run.py

@@ -24,19 +24,31 @@ COLUMNS = ('step', 'y', 'v', 'cme_F', 'cme_U', 'cme_L', 'wen_F', 'wen_U', 'wen_L
 ROOT = Path(__file__).resolve().parents[2]
 
 
-def config():
-    return json.loads(Path(__file__).with_name('config.json').read_text())
+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.):
+    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 * D * D / 4.
+        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
@@ -55,15 +67,15 @@ class Structure:
                 -self.omega**2*self.es*z + (self.ec - self.a*self.es)*self.v)
 
 
-def summarize(rows):
+def summarize(rows, diameter=D, block_steps=10*PERIOD):
     windows = []
-    for start in range(0, len(rows), 10*PERIOD):
-        block = rows[start:start + 10*PERIOD]
+    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()/D),
-                            y_rms_D=float(np.sqrt(np.mean(block[:, 1]**2))/D),
-                            y_peak_D=float(np.abs(block[:, 1]).max()/D),
+                            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
@@ -76,16 +88,26 @@ def main():
     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()
+    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()
-    s = Structure(*STRUCTURES[args.case], zeta=0.)
-    warmup, count = (32, 32) if args.smoke else (WARMUP, STEPS)
+    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
@@ -99,21 +121,26 @@ def main():
                   structure=dict(mass=s.mass, stiffness=s.stiffness, damping=s.damping,
                                  mass_ratio_three_disks=s.mass_ratio,
                                  zeta=s.zeta, fn=s.frequency, Tn=1./s.frequency, y0=0., v0=0.),
-                  warmup_steps=warmup, requested_steps=count, dt=1., D=D, U=U, nu=NU,
-                  duration_Tn=count*s.frequency, duration_tU_D=count*U/D,
-                  summary_block_steps=10*PERIOD,
+                  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.)
-            envelope = Translation((0., 0., -D, D), .05, shared=True)
-            bodies = [sim.add_body('circle', c, D/2, translation=envelope) for c in CENTRES]
+            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
@@ -123,8 +150,10 @@ def main():
             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) > D:
+                if not all(map(math.isfinite, (y1, v1))) or abs(y1) > args.D:
                     status = 'stopped_nonfinite_or_travel_bound'
                     break
                 if not args.fixed:
@@ -143,10 +172,10 @@ def main():
                 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:
+                if completed % period == 0:
                     rows.flush()
                     print(json.dumps(dict(event='progress', step=completed, t_Tn=completed*s.frequency,
-                                          y_D=s.y/D, elapsed_s=time.monotonic()-started)), flush=True)
+                                          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:
@@ -155,7 +184,7 @@ def main():
     finally:
         rows.flush()
         report.update(status=status, completed_steps=completed, elapsed_s=time.monotonic()-started,
-                      windows=summarize(rows[:completed]),
+                      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.')

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_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_steady_pinball.py

@@ -82,3 +82,47 @@ def test_predeclared_comparisons_change_structure_only():
     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) == []