diff --git a/docs/validation/leontini_viv_predeclare.md b/docs/validation/leontini_viv_predeclare.md index 4cbdfb2..aecb4c1 100644 --- a/docs/validation/leontini_viv_predeclare.md +++ b/docs/validation/leontini_viv_predeclare.md @@ -45,8 +45,8 @@ this page is jointly agreed. F_last → candidate (y_{n+1}, v_{n+1}) by the update below → `set_body` with offset y_{n+1} and wall velocity (y_{n+1} − y_n)/dt with lattice dt = 1 (interval average; a different quantity from the endpoint v_{n+1}) → `run(1)` → only on success: - commit y_n ← y_{n+1}, v_n ← v_{n+1}, F_last ← force read for that step. Any - failure leaves the Simulation failed (existing lifecycle). Whether the + commit y_n ← y_{n+1}, v_n ← v_{n+1}, F_last ← force read for that step (failure + semantics: see the window section). Whether the measured force represents the interval's start, midpoint or end is not asserted; the one-step lag is declared, not corrected. - Initial state: fluid at rest except the uniform inlet; y_0 = v_0 = 0; diff --git a/research/viv_re150/run.py b/research/viv_re150/run.py new file mode 100644 index 0000000..31457c6 --- /dev/null +++ b/research/viv_re150/run.py @@ -0,0 +1,110 @@ +"""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() diff --git a/research/viv_re150/viv.py b/research/viv_re150/viv.py new file mode 100644 index 0000000..a06ca9f --- /dev/null +++ b/research/viv_re150/viv.py @@ -0,0 +1,121 @@ +"""Pure logic of the predeclared Re 150 free transverse VIV case (stdlib only). + +Contract: docs/validation/leontini_viv_predeclare.md. Oscillator update with a +constant force over one lattice step, response-cycle tracking with fixed y = 0 +crossings, eligibility, 10-cycle blocks, settling/scoring with deterministic +restart, and the conditional seed. +""" +from __future__ import annotations + +import math + +RE, U_STAR, MASS_NUMBER = 150.0, 5.0, 2.0 # m/(rho D^2) = 2 -> m* = 8/pi +TARGET_A, BAND = 0.520, 0.10 # Qin §4.5 maximum amplitude, ±10 % +ELIGIBLE, BLOCK, LAST5 = 0.01, 10, 5 # in D; cycles; cycles +WINDOW_TOL, LAST5_TOL = 0.02, 0.01 +CYCLE_CAP, STEP_CAP_TN, SEED_TN, Y_LIMIT = 150, 200, 40, 1.0 + + +class Structure: + """Undamped spring-mass, exact harmonic update over dt = 1 with constant force.""" + + def __init__(self, D, u): + self.m = MASS_NUMBER * D * D # lattice rho = 1 + self.f_n = u / (U_STAR * D) + self.w = 2.0 * math.pi * self.f_n + self.k = self.m * self.w * self.w + self.y = self.v = 0.0 + + def candidate(self, force): + ye = force / self.k + c, s = math.cos(self.w), math.sin(self.w) + return (ye + (self.y - ye) * c + self.v / self.w * s, + -(self.y - ye) * self.w * s + self.v * c) + + def seed_velocity(self, D): + return self.w * ELIGIBLE * D # amplitude ~ .01 D + + +class Cycles: + """Upward y = 0 crossings, eligibility, blocks, settle/score with restart.""" + + def __init__(self, D): + self.D = D + self.reset_proof() + self.counted = 0 + self.log = [] # one dict per completed cycle + self._cur = None + self._prev = 0.0 + self._since = 0 + + def reset_proof(self): + self.amps = [] # eligible amplitudes since proof start + self.blocks = [] # block means + self.settled_at = None # block index where R held + self.window = None # completed scoring cycles (list of dicts) + + def discard_current(self): + self._cur, self._since = None, 0 + + def rule(self, k): + """R(k): block-mean change ≤ 2 % and last-five spread ≤ 1 %.""" + if k < 1: + return False + w0, w1 = self.blocks[k - 1], self.blocks[k] + last = self.amps[(k + 1) * BLOCK - LAST5:(k + 1) * BLOCK] + a5 = sum(last) / LAST5 + return abs(w1 - w0) / w0 <= WINDOW_TOL and max(abs(a - a5) for a in last) / a5 <= LAST5_TOL + + def update(self, y): + """Feed y after a committed step; returns 'scored' when a stable window completes.""" + out = None + crossing = self._prev < 0.0 <= y and self._since >= 2 + if crossing: + if self._cur is not None: + out = self._close() + self._cur, self._since = dict(ymax=y, ymin=y, ysum=0.0, n=0), 0 + if self._cur is not None: + c = self._cur + c["ymax"], c["ymin"] = max(c["ymax"], y), min(c["ymin"], y) + c["ysum"] += y + c["n"] += 1 + self._since += 1 + self._prev = y + return out + + def _close(self): + c = self._cur + a = 0.5 * (c["ymax"] - c["ymin"]) / self.D + cyc = dict(ymax=c["ymax"] / self.D, ymin=c["ymin"] / self.D, + offset=c["ysum"] / c["n"] / self.D, amp=a, steps=c["n"], eligible=a >= ELIGIBLE) + self.log.append(cyc) + if not cyc["eligible"]: + return None + self.counted += 1 + self.amps.append(a) + if self.settled_at is not None: + self.window.append(cyc) + if len(self.amps) % BLOCK: + return None + self.blocks.append(sum(self.amps[-BLOCK:]) / BLOCK) + k = len(self.blocks) - 1 + if self.settled_at is None: + if self.rule(k): + self.settled_at, self.window = k, [] + return None + if not self.rule(k): # scoring block failed: restart from it + amps, blocks = self.amps[-BLOCK:], self.blocks[-1:] + self.reset_proof() + self.amps, self.blocks = amps, blocks + return None + return "scored" if k == self.settled_at + 2 else None + + +def score(window): + """Target metric and reported quantities over the 20 scoring cycles.""" + a_max = max(max(abs(c["ymax"]), abs(c["ymin"])) for c in window) + return dict(A_max=a_max, mean_amp=sum(c["amp"] for c in window) / len(window), + ymax=max(c["ymax"] for c in window), ymin=min(c["ymin"] for c in window), + offset=sum(c["offset"] for c in window) / len(window), + T_resp=sum(c["steps"] for c in window) / len(window), + within_band=abs(a_max - TARGET_A) <= BAND * TARGET_A) diff --git a/tests/unit/test_viv_logic.py b/tests/unit/test_viv_logic.py new file mode 100644 index 0000000..c7513cd --- /dev/null +++ b/tests/unit/test_viv_logic.py @@ -0,0 +1,88 @@ +"""CPU tests of the predeclared VIV logic (research/viv_re150/viv.py; stdlib only).""" +import importlib.util +import math +from pathlib import Path + +_spec = importlib.util.spec_from_file_location( + "viv_logic", Path(__file__).resolve().parents[2] / "research/viv_re150/viv.py") +viv = importlib.util.module_from_spec(_spec) +_spec.loader.exec_module(viv) + +D, U = 40, 0.04 + + +def _feed(cyc, amps, period=100): + """Feed sine cycles of the given amplitudes (in D); return the first non-None result.""" + n = 0 + for a in amps: + for i in range(period): + out = cyc.update(a * D * math.sin(2 * math.pi * (i + 1) / period)) + n += 1 + if out: + return out, n + return None, n + + +def test_structure_exact_update_conserves_free_oscillation_and_equilibrium(): + s = viv.Structure(D, U) + assert math.isclose(1 / s.f_n, 5000.0) and math.isclose(s.m, 2 * D * D) + s.y, s.v = 0.3 * D, 0.0 + e0 = 0.5 * s.k * s.y ** 2 + for _ in range(12345): + s.y, s.v = s.candidate(0.0) + assert math.isclose(0.5 * s.k * s.y ** 2 + 0.5 * s.m * s.v ** 2, e0, rel_tol=1e-9) + s.y, s.v = 2.0 / s.k, 0.0 + assert all(math.isclose(a, b, abs_tol=1e-12) for a, b in zip(s.candidate(2.0), (2.0 / s.k, 0.0))) + assert math.isclose(s.seed_velocity(D) / s.w, viv.ELIGIBLE * D) + + +def test_constant_amplitude_settles_then_scores_twenty_cycles(): + cyc = viv.Cycles(D) + out, _ = _feed(cyc, [0.5] * 60) + assert out == "scored" and len(cyc.window) == 20 and cyc.settled_at == 1 + r = viv.score(cyc.window) + assert math.isclose(r["A_max"], 0.5, rel_tol=1e-3) and abs(r["offset"]) < 1e-3 + assert r["T_resp"] == 100 and r["within_band"] # |.500 - .520| <= .052 + + +def test_small_cycles_are_not_counted_and_y0_is_not_a_crossing(): + cyc = viv.Cycles(D) + assert cyc.update(0.0) is None and cyc.log == [] + out, _ = _feed(cyc, [0.005] * 50) + assert out is None and cyc.counted == 0 and cyc.log and not any(c["eligible"] for c in cyc.log) + + +def test_scoring_failure_restarts_from_the_failing_block(): + cyc = viv.Cycles(D) + amps = [0.5] * 20 + [0.55] * 10 + [0.55] * 40 # settle at block 1, block 2 jumps 10 % + out, _ = _feed(cyc, amps) + assert out == "scored" + assert all(math.isclose(c["amp"], 0.55, rel_tol=1e-3) for c in cyc.window) + + +def test_seed_reset_discards_incomplete_cycle_and_proof(): + cyc = viv.Cycles(D) + _feed(cyc, [0.5] * 15) + cyc.update(0.2 * D) # inside an incomplete cycle + counted = cyc.counted + cyc.discard_current() + cyc.reset_proof() + assert cyc.amps == [] and cyc.blocks == [] and cyc.counted == counted + + +def test_viv_config_passes_the_real_loader(tmp_path): + import json + import sys + sys.path.insert(0, str(Path(__file__).resolve().parents[2] / "research/viv_re150")) + try: + import run + finally: + sys.path.pop(0) + from CelerisLab.config import load_lbm_config + cfg_dict, nx, ny = run.config() + path = tmp_path / "config_lbm.json" + path.write_text(json.dumps(cfg_dict)) + cfg = load_lbm_config(str(path)) + assert (cfg.nx, cfg.ny) == (nx, ny) == (3000, 2400) + assert cfg.collision == "MRT" and cfg.y_wall_bc == "free_slip" and cfg.velocity == .04 + assert math.isclose(cfg.viscosity, .04 * 40 / 150)