ClaudeandClaude Opus 5.5
707fa8b5fc
viv.py: exact undamped harmonic update under the last completed force, y = 0 response cycles with eligibility, 10-cycle blocks, settle/score with deterministic restart, conditional velocity seed, scoring metrics. run.py: Simulation-API driver (set_body + run(1)), RHS = CME + Wen from the moving diagnostics, raw CME/obs/node terms logged, fail-closed bounds, caps and input rejection, outputs outside the checkout. CPU tests for the logic and the real config loader. Predeclare: stale failure sentence removed. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
111 lines
4.6 KiB
Python
111 lines
4.6 KiB
Python
"""One-shot predeclared VIV run (docs/validation/leontini_viv_predeclare.md).
|
|
|
|
Re 150, U* = 5, m/(rho D^2) = 2, zeta = 0, D = 40 on a 75 D x 60 D native
|
|
channel (uniform inlet, NEQ outlet, free-slip y walls), cylinder at
|
|
(30 D, 30 D). Structural RHS = CME + Wen (diagnostics: raw CME, node terms).
|
|
Research driver only: Simulation API, no production change.
|
|
|
|
PYTHONPATH=src python research/viv_re150/run.py --out /path/outside/checkout
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
import argparse
|
|
import json
|
|
import math
|
|
import os
|
|
import tempfile
|
|
|
|
import numpy as np
|
|
|
|
import viv
|
|
|
|
D, U = 40, 0.04
|
|
TEMPLATE = os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "..", "src",
|
|
"CelerisLab", "configs", "config_lbm.json")
|
|
|
|
|
|
def config(D=D, U=U):
|
|
nx, ny = 75 * D, 60 * D
|
|
with open(TEMPLATE) as fh:
|
|
cfg = json.load(fh)
|
|
cfg["grid"].update(lattice_model="D2Q9", nx=nx, ny=ny, nz=1)
|
|
cfg["physics"].update(data_type="FP32", viscosity=U * D / viv.RE, velocity=U, rho=1.0)
|
|
cfg["method"].update(collision="MRT", streaming="double_buffer", store_precision="FP32",
|
|
ddf_shifting=False, y_wall_bc="free_slip")
|
|
cfg["method"]["les"]["enabled"] = False
|
|
cfg["method"]["inlet"]["profile"] = "uniform"
|
|
return cfg, nx, ny
|
|
|
|
|
|
def main():
|
|
ap = argparse.ArgumentParser()
|
|
ap.add_argument("--out", required=True)
|
|
args = ap.parse_args()
|
|
if os.path.realpath(args.out).startswith(os.path.realpath(os.path.join(os.path.dirname(__file__), "..", ".."))):
|
|
raise SystemExit("--out must be outside the checkout")
|
|
os.makedirs(args.out, exist_ok=True)
|
|
cfg, nx, ny = config()
|
|
s, cyc = viv.Structure(D, U), viv.Cycles(D)
|
|
T_n = 1.0 / s.f_n
|
|
step_cap, n_seed = int(round(viv.STEP_CAP_TN * T_n)), int(round(viv.SEED_TN * T_n))
|
|
centre = (30.0 * D, 0.5 * (ny - 1))
|
|
|
|
from CelerisLab.body.motion import Translation
|
|
from CelerisLab.simulation import Simulation
|
|
with tempfile.NamedTemporaryFile("w", suffix=".json", delete=False) as fh:
|
|
json.dump(cfg, fh)
|
|
sim = Simulation(fh.name)
|
|
rows, status, seed_step, f_last = [], "running", None, 0.0
|
|
try:
|
|
body = sim.add_body("circle", centre, 0.5 * D,
|
|
translation=Translation((0.0, 0.0, -1.5 * D, 1.5 * D), 0.5))
|
|
sim.initialize()
|
|
for n in range(step_cap):
|
|
if n == n_seed and max((abs(r[1]) for r in rows[-int(5 * T_n):]), default=0.0) < viv.ELIGIBLE * D:
|
|
s.v += s.seed_velocity(D) # after the quiet check, before the update
|
|
seed_step = n
|
|
cyc.discard_current()
|
|
cyc.reset_proof()
|
|
y1, v1 = s.candidate(f_last)
|
|
if not (math.isfinite(y1) and math.isfinite(v1)) or abs(y1) > viv.Y_LIMIT * D:
|
|
status = "fail_bound_or_nonfinite"
|
|
break
|
|
try:
|
|
sim.set_body(body, offset=(0.0, y1), vy=y1 - s.y)
|
|
except ValueError as exc: # input rejection: stop, no retry
|
|
status = f"rejected: {exc}"
|
|
break
|
|
sim.run(1) # device failures raise and mark the Simulation failed
|
|
d = sim.moving_diagnostics()[0]
|
|
f_obs = sim.read_force(body)[1]
|
|
f_rhs = float(d[1] + d[3]) # CME + Wen (y)
|
|
if not math.isfinite(f_rhs):
|
|
status = "fail_nonfinite_force"
|
|
break
|
|
s.y, s.v, f_last = y1, v1, f_rhs # commit only after a successful step
|
|
rows.append((n, s.y, s.v, f_rhs, float(d[1]), float(d[3]), float(f_obs),
|
|
float(d[7] - d[10]), float(d[13])))
|
|
if cyc.update(s.y) == "scored":
|
|
status = "scored"
|
|
break
|
|
if cyc.counted >= viv.CYCLE_CAP:
|
|
status = "timeout_response_cycles"
|
|
break
|
|
else:
|
|
status = "timeout_native_steps"
|
|
finally:
|
|
np.save(os.path.join(args.out, "steps.npy"), np.asarray(rows, float))
|
|
report = dict(status=status, steps=len(rows), t_U_over_D=len(rows) * U / D, T_n=T_n,
|
|
seed_step=seed_step, counted_cycles=cyc.counted, cycles=cyc.log,
|
|
columns=["n", "y", "v", "F_rhs_cme_wen", "F_cme", "F_wen", "F_obs",
|
|
"node_cover_minus_uncover_y", "no_donor"])
|
|
if status == "scored":
|
|
report["score"] = viv.score(cyc.window)
|
|
with open(os.path.join(args.out, "report.json"), "w") as fh:
|
|
json.dump(report, fh, indent=1)
|
|
sim.close()
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|