research: VIV driver — fail closed on any non-finite recorded output

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

docs/validation/leontini_viv_predeclare.md

@@ -0,0 +1,146 @@
+# Single-cylinder free transverse VIV — PREDECLARATION (agreed at 1abfe954)
+
+Stage goal (Yanqi via TRON): prescribed Dütsch case + one VIV case both
+"reasonably close" to a predeclared reference = development-stage success, then
+cleanup. The frozen Dütsch settling failure (0/28 pairs at 31.25 T) stays on
+record. Practical closeness is not a formal validation claim. Agreed by
+Claude, Grok and Astra at 1abfe954; implemented in research/viv_re150/.
+
+## Primary case (one case)
+
+- Reference: Qin et al., arXiv:2003.12186, §4.5 (maximum amplitudes at
+  U_red = 3…8: .075, .551, .520, .459, .370, .084 D in the text, as quoted by
+  Grok/Astra) for the target; setup from Kolahdouz et al., arXiv:2003.12046,
+  §4.1.1 (read 2026-10-30 via ar5iv): U* sweep 3–8 with amplitudes only in
+  Fig. 3(b) (no text values), domain [−30, 45] D × [−30, 30] D with the
+  cylinder "initially at rest and centered at the origin", uniform inflow on
+  the left, zero normal traction and tangential velocity at the outlet, zero
+  normal velocity and tangential traction on y = ±30 D; no initial
+  perturbation and no total simulated time stated.
+- Re = 150, U* = U/(f_n D) = 5, m/(rho D²) = 2 (m* = rho_s/rho_f = 8/pi),
+  zeta = 0, f_n = sqrt(K/M)/(2 pi) in vacuum. Transverse DOF only.
+- Target: A_max = .520 D (peak amplitude, Qin §4.5). Multi-code reference
+  (Kolahdouz ILE-NS, Qin II-LBM). Reference uncertainty: not stated; the
+  spread between codes is not quantified here.
+
+## Our reproduction (declared differences)
+
+- D = 40, MRT, dense double buffer; U_lat = .04, nu = U D / 150 = .010667,
+  tau = .532.
+- Box [0, 75 D] × [0, 60 D] (3000 × 2400) with the cylinder centre at
+  (30 D, 30 D): 30 D from the inlet, 45 D from the outlet, centred in y
+  (Kolahdouz §4.1.1). Qin's own box ([0, 40 D] × [0, 20 D], centre (10 D, 10 D))
+  is the target paper's, not ours.
+- Boundaries: uniform Zou–He inlet, NEQ outlet (rho = 1, backflow clamp),
+  free-slip y walls. Kolahdouz uses a traction outlet: declared mismatch.
+- Qin's own VIV box (40 D × 20 D) is not used; a smaller box would be a
+  declared cost cut with more blockage, not "the same case".
+- Allowed y offset ±1.5 D; peak transverse speed ≈ .026 cells/step (bound .5).
+
+## Structure, coupling and time levels
+
+- M = 2 rho D² (lattice rho = 1), f_n = U/(5 D) = 2e-4 per step, K = M (2 pi f_n)²,
+  zeta = 0. Equation of motion uses the structural mass only.
+- Per-step sequence (binding): known (y_n, v_n) and the last completed force
+  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 (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 initialized by `Simulation.initialize()` to the uniform
+  inlet velocity U everywhere (the code's actual initial condition; the
+  source gives no fluid initial condition), body at rest at the declared
+  centre, y_0 = v_0 = 0, F_last = 0.
+- Structural update over one step with constant force F = F_last (exact
+  harmonic solution, zeta = 0): y_e = F/K, w = 2 pi f_n,
+  y^{n+1} = y_e + (y^n − y_e) cos w + (v^n / w) sin w,
+  v^{n+1} = −(y^n − y_e) w sin w + v^n cos w.
+- Primary RHS (single coupled run): CME + Wen (theory-motivated stage choice:
+  Galilean-invariant link correction; G4 could not separate it from CME).
+  Raw CME and the node terms are logged as diagnostics of that trajectory;
+  the node impulse is not added. A run with CME-only as RHS is a separate,
+  labelled sensitivity case, decided separately. Neither RHS is claimed proven.
+
+## Window, failure semantics, acceptance (engineering, development stage)
+
+Clocks: lattice dt = 1; T_n = 1/f_n = 5,000 steps. Response cycles are
+counted from the measured signal; T_n is used only for the seed check.
+
+- Start: body at rest at the declared centre (the only initial condition the
+  source states), y_0 = v_0 = 0, F_last = 0; fluid uniform at U (as
+  initialized by the code; our choice, not the source's).
+- Seed (conditional, one time, declared): at step n_seed = 200,000
+  (t = 40 T_n), if max |y_t| over states t in [175,000, 200,000) is below
+  .01 D, then
+  in that step, after the check and before the candidate update, set
+  v_n += 2 pi f_n (.01 D) (F_last unchanged) and proceed normally; the
+  geometry moves continuously through the update. `seed_applied`, step and
+  size are logged; never repeated or changed.
+- Cycle boundary (discrete, fixed threshold): step n closes a cycle iff
+  y[n−1] < 0 and y[n] ≥ 0 (upward crossing of the equilibrium y = 0);
+  y_0 = 0 is not a crossing; no boundary fewer than 2 samples after the
+  previous one.
+- Per cycle: positive peak max(y), negative peak min(y), mean offset
+  (time-mean of y), amplitude a = (max − min)/2 (the only amplitude used for
+  settling).
+- Eligibility: a completed cycle is counted (toward the settling proof and
+  the 150 budget) only if a ≥ .01 D; ineligible cycles are logged, never
+  counted, so near-zero motion cannot settle or pass.
+- Seed reset: if the seed is applied, the incomplete cycle at n_seed is
+  discarded and the settling proof starts empty; the 150 budget and the
+  native cap keep running from t = 0.
+- Windows: adjacent non-overlapping blocks of 10 counted cycles;
+  W_k = mean amplitude of block k (W_k ≥ .01 D by eligibility, so no zero
+  denominator).
+- Block rule R(k), k ≥ 2: |W_k − W_{k−1}| / W_{k−1} ≤ .02 and, for the last 5
+  cycles of block k, max_i |a_i − ā_5| / ā_5 ≤ .01 with ā_5 their mean.
+- Settled at the end of block k when R(k) holds. Scoring = the next two
+  blocks k+1, k+2 (20 counted cycles); both must satisfy R. If one fails,
+  that block becomes the new reference: the settling proof restarts from it
+  and continues within the caps (deterministic; no other retry). A pass
+  requires one complete stable 20-cycle scoring window before a cap; success
+  stops after its 20th cycle.
+- Caps (both always armed; the first hit wins and is reported): 150 counted
+  response cycles including the 20 scoring cycles (so settling must complete
+  by counted cycle 130) → `timeout_response_cycles`; 1,000,000 native steps
+  (200 T_n) → `timeout_native_steps`.
+- Durations reported: total steps, t U/D, counted and uncounted cycles,
+  measured T_resp (mean boundary interval); N_resp / f_n only as a labelled
+  secondary clock.
+- Growth is not a failure: initial VIV growth is physical. Checked every
+  step (FAIL-CLOSED): non-finite y, v or force; |y| > 1.0 D (the tighter
+  gate; about twice the largest reference peak, .551 D) or the ±1.5 D offset
+  box reached; no_donor > 0; a failed upload/step/sync.
+- On any failure the driver stops and reports; it never retries,
+  sub-iterates, retunes the RHS or parameters, or starts a second case.
+  Device-side failures (upload/step/sync, no_donor) mark the Simulation
+  failed through the existing lifecycle; input rejection by `set_body`
+  (offset box, speed, step bound) and timeouts end the driver run without
+  that mark.
+- Target metric over the 20 scoring cycles: A_max/D = max over those cycles
+  of max(|positive peak|, |negative peak|) about y = 0 (Qin reports "maximum
+  amplitude"); also reported: mean amplitude, positive and negative peaks,
+  mean offset, f/f_n (no target).
+- Acceptance: |A_max − .520 D| ≤ .052 D (±10 %), stated together with the
+  boundary mismatch and the absence of a quoted reference uncertainty. A pass
+  is "reasonably close" for the stage: not grid convergence, not force-RHS
+  validation, and it does not touch the Dütsch settling record.
+- Driver reads force and the small diagnostics per step only; field reads
+  only at a few predeclared snapshots. No runtime claim: there is no measured
+  timing on this grid yet (the D40/60D receipt is a different, sampled case).
+
+## Plumbing smoke (not the benchmark)
+
+`run.py --smoke --D 20 --steps 32` runs exactly 32 steps of the same BC
+package and coupling loop (no seed, settle or score; never scored);
+status `smoke_ok` iff every step is finite, no_donor stays 0 and the body
+stays inside its box. The physical case is only the default invocation
+(`--D 40`, no `--smoke`).
+
+## Queued, not in this stage
+
+Leontini et al. JFS 22 (2006) Re = 200, m* = 10, zeta = 0.01, U* = 5.0 (needs
+the owners' c_crit convention and a Fig. 1 digitization). D80, literature
+claims, pinball, AMR, BC redesign, threshold changes.

research/viv_re150/run.py

@@ -0,0 +1,134 @@
+"""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("--D", type=int, default=D)
+    ap.add_argument("--smoke", action="store_true",
+                    help="plumbing smoke only: exactly --steps steps (<= 64), never scored")
+    ap.add_argument("--steps", type=int, default=0)
+    args = ap.parse_args()
+    if args.smoke != bool(args.steps) or args.steps > 64 or args.steps < 0:
+        raise SystemExit("--smoke requires --steps N with 1 <= N <= 64 (and only with --smoke)")
+    if args.D != D and not args.smoke:
+        raise SystemExit("the physical case is predeclared at D = 40")
+    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")
+    Dn = args.D
+    cfg, nx, ny = config(Dn)
+    s, cyc = viv.Structure(Dn, U), viv.Cycles(Dn)
+    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.smoke:
+        step_cap = args.steps
+    centre = (30.0 * Dn, 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 * Dn,
+                            translation=Translation((0.0, 0.0, -1.5 * Dn, 1.5 * Dn), 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 not args.smoke and n == n_seed and \
+                    max(abs(r[1]) for r in rows[n_seed - quiet - 1:n_seed - 1]) < viv.ELIGIBLE * Dn:
+                s.v += s.seed_velocity(Dn)         # 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 * Dn:
+                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)
+            row = (n, y1, v1, f_rhs, float(d[1]), float(d[3]), float(f_obs),
+                   float(d[7] - d[10]), float(d[13]))
+            if not viv.finite(row):                # every recorded output, not only the RHS
+                status = "fail_nonfinite_output"
+                break
+            s.y, s.v, f_last = y1, v1, f_rhs       # commit only after a successful step
+            rows.append(row)
+            if args.smoke:
+                continue                           # no settle / score in the plumbing smoke
+            if cyc.update(s.y) == "scored":
+                status = "scored"
+                break
+            if cyc.counted >= viv.CYCLE_CAP:
+                status = "timeout_response_cycles"
+                break
+        else:
+            status = ("smoke_ok" if all(r[8] == 0.0 for r in rows) else "smoke_no_donor") \
+                if args.smoke 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=args.smoke, D=Dn, steps=len(rows), t_U_over_D=len(rows) * U / Dn, 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.smoke:
+            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()

research/viv_re150/viv.py

@@ -0,0 +1,126 @@
+"""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
+
+
+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)

tests/unit/test_viv_logic.py

@@ -0,0 +1,96 @@
+"""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)
+
+
+def test_finite_rejects_any_nonfinite_recorded_output():
+    row = (3, .1, .2, .3, .4, .5, .6, .7, 0.0)
+    assert viv.finite(row)
+    for i in range(len(row)):
+        for bad in (float("nan"), float("inf"), -float("inf")):
+            assert not viv.finite(row[:i] + (bad,) + row[i + 1:])