"""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()