CelerisLab/research/viv_re150/run.py

"""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)
    ap.add_argument("--max-steps", type=int, default=0,
                    help="plumbing smoke only: stop after N steps; never scored")
    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)
    for name in ("report.json", "steps.npy"):  # never overwrite earlier evidence
        if os.path.exists(os.path.join(args.out, name)):
            raise SystemExit(f"{args.out} already holds {name}; use a fresh output directory")
    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))
    quiet = int(round(5 * T_n))
    if args.max_steps:
        step_cap = min(step_cap, args.max_steps)
    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):
            # rows[k] holds the state y_{k+1}; the declared window is y_t, t in [n_seed - 5 T_n, n_seed)
            if n == n_seed and max(abs(r[1]) for r in rows[n_seed - quiet - 1:n_seed - 1]) < 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 = "smoke_stopped" if args.max_steps else "timeout_native_steps"
    except BaseException as exc:
        status = f"error: {exc!r}"
        raise
    finally:
        np.save(os.path.join(args.out, "steps.npy"), np.asarray(rows, float))
        report = dict(status=status, smoke=bool(args.max_steps), 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" and not args.max_steps:
            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()