ClaudeandClaude Opus 5.5
48ee0982c9
Tests load research drivers via importlib under unique module names instead of a bare 'import run' (shared-process collision between the Dütsch and VIV drivers). lifecycle.md qualifies immutable geometry and the host-only get_flags for translating circles; supported_settings.md points to the dense-only moving-circle feature; leontini_viv_predeclare.md is renamed viv_re150_predeclare.md (Re 150 Kolahdouz/Qin case) with references updated. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
127 lines
4.9 KiB
Python
127 lines
4.9 KiB
Python
"""Pure logic of the predeclared Re 150 free transverse VIV case (stdlib only).
|
|
|
|
Contract: docs/validation/viv_re150_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
|
|
|
|
|
|
def finite(values):
|
|
"""True iff every recorded value is finite (fail-closed output check)."""
|
|
return all(math.isfinite(float(v)) for v in values)
|
|
|
|
|
|
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)
|