Add matched diffusive resolution cases for pinball response

docs/steady_pinball_fiv.md

@@ -5,7 +5,8 @@ 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 existing backflow clamp/blending. Initial fluid is
+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,
@@ -71,3 +72,79 @@ This stage makes no claim of physical grid convergence or vibration suppression.
 Run from checkout with PYTHONPATH=src and the existing Python environment:
 `python research/steady_pinball/run.py --smoke --out <new-external-directory>`
 then separately `--out <new-free-directory>` and `--fixed --out <new-fixed-directory>`.
+
+## Zero-damping development follow-up (predeclared 2026-10-31)
+
+The original d0bcefa run and its zeta=.01 results remain historical evidence.
+The current runner exposes three fixed `--case` choices, all with C=0 as requested:
+
+| Case | M / (rho 3 pi D^2/4) | fn | M | K | Purpose |
+|---|---:|---:|---:|---:|---|
+| Z0 | 10 | .0001 | 9424.77796 | .00372075320 | Remove damping only |
+| M2 | 2 | .0001 | 1884.95559 | .000744150640 | Reduce mass and stiffness together at fixed fn |
+| R2 | 2 | .00006 | 1884.95559 | .000267894231 | Change stiffness only relative to M2 |
+
+R2 deliberately places fn at the retained fixed-body lift peak, f=.00006
+(St=.12 on mean U). This is an informed development choice, not independent
+physical validation or a promise of lock-in. All geometries, fluid settings,
+512000-step fixed warmup, y0=v0=0, force time level and unfiltered RHS are unchanged.
+No run result may select a replacement parameter or an automatic extension.
+
+Each case has a 600000-step cap, the same 300 D/U physical duration. This equals
+60 natural periods for Z0/M2 and 36 for R2. Compare the common final200000 native
+steps, and six100000-step blocks, retaining actual natural/response frequencies
+and ongoing drift. Do not call those blocks ten R2 periods. No settled-VIV claim
+follows from a larger displacement. Report all three outcomes, including quiet
+responses or stopped cases, alongside the existing fixed and damped-free records.
+
+Unchanged stops: |y|>D, invalid move/speed, nonfinite values or no_donor. A 1800s
+per-case wall cap and GPU headroom check bound resource use. Preserve partial rows;
+do not add damping, clip trajectories, filter forces, or relax any stop. Each
+case adds72MB of scalar history, no full fields. The separate Legacy/modern
+baseline diagnosis does not change this FSI protocol.
+
+## Matched D40 resolution comparison (predeclared 2026-10-31)
+
+The first finer case is M2, with exactly the same mass ratio2, reduced velocity5,
+zero structural damping, initial rest state and physical geometry as the retained
+D20 M2. Targeting a development response of0.5--1D does not imply this case will
+reach it. Do not change the structure within this comparison or discard a smaller
+response. D60 is a supported scaling choice, not an automatically scheduled run.
+
+For s=D/20 use dx'=dx/s, dt'=dt/s^2, U'=.01/s, nu'=.004, rho'=1,
+fn'=fn/s^2, M'=s^2 M and K'=K/s^2. The same physical force per unit span
+has F'_lat=F_lat/s. For D40 this gives U=.005, fn=.000025,
+M=7539.82237, K=.000186037660, Re_D=50, and Tn=40000 steps.
+The y walls are at .5 and ny-1.5, while x inlet/outlet source nodes are at0
+and nx-1. Therefore map x'=s*x and y'=.5+s*(y-.5), with
+nx'=s*(1280-1)+1, ny'=s*(512-2)+2. D40 uses2559x1022 and centres
+(1200,510.5),(1252,540.5),(1252,480.5); every radius, probe and envelope scales
+with D. The parabolic profile retains the same normalized wall coordinate.
+This preserves the existing ghost-source boundary convention, without asserting
+that it is an exact continuum open-boundary location.
+
+Warmup=2048000 and response=2400000 native steps at D40: the SAME256 and300
+D/U durations as D20, including60 natural periods after release. Compare the
+six50-D/U blocks and the common final100-D/U window. Report max absolute y/D,
+half-range, RMS, frequency and drift; peak-to-peak is not amplitude. No settling
+claim follows from completion. Keep raw individual CME/Wen/node terms and all
+crossing/fallback counters. Normalize net force by .5*rho*U^2*D in each case.
+Report native crossing/non-crossing force increments separately, plus D40
+four-step block means aligned with D20 intervals (diagnostic only, never replace
+the raw spikes). Report final-window high-frequency variance above f*D/U=1,
+using demeaned FFT Parseval weights, with the same physical cutoff for both cases.
+Different motion histories confound spike-size comparisons; this is simultaneous
+space/time/Mach refinement, not proof of a unique coarse-grid cause.
+
+Only one full D40 M2 is initially scheduled, after CPU checks and a32+32-step
+smoke. Stop on the unchanged geometry/speed/nonfinite/no_donor checks or a7200s
+whole-case wall cap; save completed rows and status. The velocity cap scales to
+.025 so its physical bound is unchanged. No force filtering, trajectory clipping,
+automatic parameter search or extension. Scalar history is288MB; no full fields
+or DDF snapshots are written. Two FP32 population buffers require about188MiB,
+before flags, moving tables and context overhead; measure actual peak usage in
+the smoke and full run. Keep all D20 outcomes. No production source changes.
+
+Commands (one GPU and serial compilation, existing environment):
+`python research/steady_pinball/run.py --case M2 --D 40 --smoke --out <new-smoke-dir>`
+then `python research/steady_pinball/run.py --case M2 --D 40 --wall-seconds 7200 --out <new-full-dir>`.

research/steady_pinball/run.py

@@ -1,6 +1,6 @@
 """Uncontrolled paper-geometry steady Re_D=50, one transverse assembly DOF.
 
-One fixed protocol; --smoke performs 32 warmup + 32 observation steps only.
+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
@@ -18,21 +18,39 @@ 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():
-    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):
-        self.mass = 10. * 3. * math.pi * D * D / 4.
-        self.omega = 2. * math.pi / PERIOD
-        self.zeta = .01
+    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.
@@ -49,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
@@ -68,16 +86,28 @@ def main():
     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()
+    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()
-    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
@@ -87,23 +117,30 @@ def main():
     col = {k: names.index(k) for k in names}
     status, completed, f_last = 'running', 0, 0.
     started = time.monotonic()
-    report = dict(source=source, fixed=args.fixed, smoke=args.smoke, columns=COLUMNS,
+    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,
-                                 zeta=s.zeta, fn=1./PERIOD, y0=0., v0=0.),
-                  warmup_steps=warmup, requested_steps=count, dt=1., D=D, U=U, nu=NU,
-                  duration_Tn=count/PERIOD, duration_tU_D=count*U/D,
+                                 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.)
-            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
@@ -113,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:
@@ -133,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='period', step=completed, t_Tn=completed/PERIOD,
-                                          y_D=s.y/D, elapsed_s=time.monotonic()-started)), flush=True)
+                    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:
@@ -145,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_steady_pinball.py

@@ -26,7 +26,7 @@ def test_actual_loader_and_paper_geometry(tmp_path):
 
 
 def test_structure_rest_equilibrium_and_energy_decay():
-    s = runner.Structure()
+    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.
@@ -52,3 +52,77 @@ def test_report_uses_completed_rows_and_net_force():
     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) == []