Add matched diffusive resolution cases for pinball response

README.md

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

docs/moving_bodies.md

@@ -30,6 +30,10 @@ cells (M1 `MovingCircleCells`), moving-wall Bouzidi with wall velocity
   mass, eq_only, no_donor, crossings) accumulated since obs were zeroed;
   `obs`/`read_force` keep their contract (CME on the body).
 - `get_flags()` reads device flags back when moving circles exist.
+- `snapshot()` / `restore()` save and restore populations only, keeping the
+  current clock, controls, observations, geometry and flags. Use them only
+  with unchanged geometry; restoring across a translation is unsupported
+  and is not rejected at runtime. They do not provide FSI rollback.
 
 *Supported:* dense double buffer, D2Q9, FP32 unshifted, SRT/MRT, circles,
 spin. *Rejected before mutation:* EsoPull, ddf shifting, D3Q19, refinement /

docs/steady_pinball_fiv.md

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

docs/validation/viv_re150_predeclare.md

@@ -29,10 +29,17 @@ Claude, Grok and Astra at 1abfe954; implemented in research/viv_re150/.
 
 - D = 40, MRT, dense double buffer; U_lat = .04, nu = U D / 150 = .010667,
   tau = .532.
-- Box [0, 75 D] × [0, 60 D] (3000 × 2400) with the cylinder centre at
-  (30 D, 30 D): 30 D from the inlet, 45 D from the outlet, centred in y
-  (Kolahdouz §4.1.1). Qin's own box ([0, 40 D] × [0, 20 D], centre (10 D, 10 D))
-  is the target paper's, not ours.
+- Nominal 75 D × 60 D lattice allocation (3000 × 2400), based on the
+  Kolahdouz §4.1.1 layout. Coordinate clarification of the existing run:
+  the driver uses `(30 D, (ny - 1)/2)` = `(1200, 1199.5)` in lattice
+  coordinates, or `(30 D, 29.9875 D)`. The earlier `(30 D, 30 D)` wording
+  described the nominal box, not the exact executed centre.
+  Free-slip reflection planes are at `y = .5` and `y = ny - 1.5 = 2398.5`;
+  their separation is 2398 cells = 59.95 D, with the centre exactly midway,
+  29.975 D from each wall. Node counts are not physical wall-to-wall lengths.
+  This corrects the description, without changing the historical geometry
+  or results. Qin's own box ([0, 40 D] × [0, 20 D], centre (10 D, 10 D)) is
+  the target paper's, not ours.
 - Boundaries: uniform Zou–He inlet, NEQ outlet (rho = 1, backflow clamp),
   free-slip y walls. Kolahdouz uses a traction outlet: declared mismatch.
 - Qin's own VIV box (40 D × 20 D) is not used; a smaller box would be a

research/duetsch_re100_kc5/protocol.py

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

research/steady_pinball/config.json

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

research/steady_pinball/run.py

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

research/viv_re150/run.py

@@ -1,8 +1,9 @@
 """One-shot predeclared VIV run (docs/validation/viv_re150_predeclare.md).
 
-Re 150, U* = 5, m/(rho D^2) = 2, zeta = 0, D = 40 on a 75 D x 60 D native
-channel (uniform inlet, NEQ outlet, free-slip y walls), cylinder at
-(30 D, 30 D). Structural RHS = CME + Wen (diagnostics: raw CME, node terms).
+Re 150, U* = 5, m/(rho D^2) = 2, zeta = 0, D = 40 on a nominal 75 D x 60 D
+lattice (uniform inlet, NEQ outlet, free-slip y walls). The cylinder is at
+(30 D, (ny - 1)/2); y-wall planes are .5 and ny - 1.5 in lattice coordinates.
+Structural RHS = CME + Wen (diagnostics: raw CME, node terms).
 Research driver only: Simulation API, no production change.
 
     PYTHONPATH=src python research/viv_re150/run.py --out /path/outside/checkout

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

@@ -82,6 +82,35 @@ def test_interval_bases_and_morison_fit_recover_known_coefficients():
     assert p.morison_fit(f, b[:, 0], b[:, 1]) == pytest.approx((2.09, 1.45), abs=1e-9)
 
 
+def test_quadrature_rule_reused_across_steps_and_periods(monkeypatch):
+    original = np.polynomial.legendre.leggauss
+    orders = []
+
+    def counted(points):
+        orders.append(points)
+        return original(points)
+
+    p._quadrature_rule.cache_clear()
+    monkeypatch.setattr(np.polynomial.legendre, "leggauss", counted)
+    try:
+        for points in (16, 8, 16):
+            for period in (T, 2 * T):
+                for n in (0, 123, period + 123):
+                    p.interval_bases(n, period, points)
+        assert orders == [16, 8]
+    finally:
+        p._quadrature_rule.cache_clear()
+
+
+def test_interval_sine_basis_matches_exact_integral():
+    for period in (T, 2 * T):
+        half_width = math.pi / period
+        for n in (0, period // 4, period // 2, period - 1, 7 * period + 123):
+            midpoint = 2 * math.pi * ((n % period) + .5) / period - math.pi / 2
+            expected = math.sin(midpoint) * math.sin(half_width) / half_width
+            assert p.interval_bases(n, period)[1] == pytest.approx(expected, rel=0, abs=2e-15)
+
+
 def test_geometry_repeats_exactly_across_cycles():
     A = p.amplitude(20)
     for m in (0, 1, T // 4, T // 2 - 1, 7 * T // 12):

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