research: predeclared Re150 VIV driver and logic (research-only)
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>
This commit is contained in:
Claude
co-authored by
Claude Opus 5.5
Claude Opus 5.5
parent
1abfe954cd
commit
707fa8b5fc
@@ -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;
|
||||
|
||||
@@ -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()
|
||||
@@ -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)
|
||||
@@ -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)
|
||||
Reference in New Issue
Block a user