"""Uncontrolled paper-geometry steady Re_D=50, one transverse assembly DOF. One fixed protocol; --smoke performs 32 warmup + 32 observation steps only. 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 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()) 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 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): windows = [] for start in range(0, len(rows), 10*PERIOD): block = rows[start:start + 10*PERIOD] 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), 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') args = ap.parse_args() 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() (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) 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, 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, 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] sim.initialize() # At release obtain the last instantaneous force, not a warmup average. remaining = warmup - 1 while remaining: 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): y1, v1 = (0., 0.) if args.fixed else s.candidate(f_last) if not all(map(math.isfinite, (y1, v1))) or abs(y1) > 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='period', step=completed, t_Tn=completed/PERIOD, y_D=s.y/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]), 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()