158 lines
7.2 KiB
Python
158 lines
7.2 KiB
Python
"""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()
|