Add matched diffusive resolution cases for pinball response

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