research: VIV driver — bounded --smoke --steps N (<=64) --D plumbing mode replaces --max-steps

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

docs/validation/leontini_viv_predeclare.md

@@ -131,6 +131,14 @@ counted from the measured signal; T_n is used only for the seed check.
   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

research/viv_re150/run.py

@@ -40,23 +40,30 @@ def config(D=D, U=U):
 def main():
     ap = argparse.ArgumentParser()
     ap.add_argument("--out", required=True)
-    ap.add_argument("--max-steps", type=int, default=0,
+    ap.add_argument("--D", type=int, default=D)
-                    help="plumbing smoke only: stop after N steps; never scored")
+    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")
-    cfg, nx, ny = config()
+    Dn = args.D
-    s, cyc = viv.Structure(D, U), viv.Cycles(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.max_steps:
+    if args.smoke:
-        step_cap = min(step_cap, args.max_steps)
+        step_cap = args.steps
-    centre = (30.0 * D, 0.5 * (ny - 1))
+    centre = (30.0 * Dn, 0.5 * (ny - 1))
 
     from CelerisLab.body.motion import Translation
     from CelerisLab.simulation import Simulation
@@ -65,18 +72,19 @@ def main():
     sim = Simulation(fh.name)
     rows, status, seed_step, f_last = [], "running", None, 0.0
     try:
-        body = sim.add_body("circle", centre, 0.5 * D,
+        body = sim.add_body("circle", centre, 0.5 * Dn,
-                            translation=Translation((0.0, 0.0, -1.5 * D, 1.5 * D), 0.5))
+                            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 n == n_seed and max(abs(r[1]) for r in rows[n_seed - quiet - 1:n_seed - 1]) < viv.ELIGIBLE * D:
+            if not args.smoke and n == n_seed and \
-                s.v += s.seed_velocity(D)          # after the quiet check, before the update
+                    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 * D:
+            if not (math.isfinite(y1) and math.isfinite(v1)) or abs(y1) > viv.Y_LIMIT * Dn:
                 status = "fail_bound_or_nonfinite"
                 break
             try:
@@ -94,6 +102,8 @@ def main():
             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 args.smoke:
+                continue                           # no settle / score in the plumbing smoke
             if cyc.update(s.y) == "scored":
                 status = "scored"
                 break
@@ -101,17 +111,18 @@ def main():
                 status = "timeout_response_cycles"
                 break
         else:
-            status = "smoke_stopped" if args.max_steps else "timeout_native_steps"
+            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=bool(args.max_steps), steps=len(rows), t_U_over_D=len(rows) * U / D, T_n=T_n,
+        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.max_steps:
+        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)