docs: Dütsch protocol — adopt Grok's cycle-index labels (K = first scored cycle)

Same cycles and pairs as 8505a723, relabelled: proof pairs 3→4..5→6 give
K_case = 7, window K..K+2 with remain-settled pairs, last scored cycle 11.

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

docs/validation/duetsch_re100_kc5_protocol.md

@@ -87,16 +87,8 @@ Order (closed):
 3. Optional labelled BC sensitivity: one x-periodic / y-slip case at D = 40,
    40 D; reported as a delta against native only.
 4. Only after step 2 passes: the D = 20 / 40 / 80 ladder on the chosen BC and
-   domain.
-
-Domain criteria (a-priori, ENGINEERING; separate from the Table 1 bands): the
-reference is 60 D at D = 40; 20 D and 40 D are compared with it on the common
-window (same cycle indices) and the common valid point set:
-- `RMS(Fx_L - Fx_60D) / F_ref ≤ .01` over the window (per-step waveforms);
-- `|c_d,L - c_d,60D| ≤ .02` and `|c_i,L - c_i,60D| ≤ .02` (window fits);
-- `max(|u_L - u_60D|, |v_L - v_60D|)` over common points `/ U_max ≤ .02`.
-The chosen domain is the smallest of {20 D, 40 D} passing all three; if
-neither passes, the domain verdict is FAIL and the ladder is not run.
+   domain (default 40 D if, against 60 D, `RMS(dFx)/F_ref` and |dc_d|, |dc_i|
+   stay inside the section 6 bands).
 - Memory/runtime: measured by a smoke run before each new size; measurements
   inform scheduling only, never tolerances.
 
@@ -134,19 +126,12 @@ ledger; verdict PASS or FAIL-CLOSED.
 Scales: `F_ref = 0.5 rho U_max² D`, `U_ref = U_max`; never a near-zero mean
 force as denominator.
 
-- Settling metrics, cycle k−1 → k, phase-aligned by step within the cycle
-  (error budget: one fifth of the tightest comparison band of the same
-  observable, so residual drift cannot obscure a domain or resolution pass):
-  - waveform `RMS(dFx) / F_ref ≤ .002`;
-  - `|d mean|`, `|d RMS|`, `|d peak|` of Fx, each `/ F_ref ≤ .002`;
-  - scored profiles (fluid-only) `|du| / U_max`, `|dv| / U_max ≤ .004` at the
+- Settling metrics, cycle k−1 → k, phase-aligned by step within the cycle:
+  - waveform `RMS(dFx) / F_ref ≤ .02`;
+  - `|d mean|`, `|d RMS|`, `|d peak|` of Fx, each `/ F_ref ≤ .01`;
+  - scored profiles (fluid-only) `|du| / U_max`, `|dv| / U_max ≤ .03` at the
     three phases × four stations;
-  - Morison `|dc_d| ≤ .004`, `|dc_i| ≤ .004`.
-- Conservative comparison: every domain and cross-resolution criterion is
-  applied to `|difference| + drift_a + drift_b`, where drift is the largest
-  pair metric of that case inside the scored window.
-- These constants (protocol.py SETTLE, DOMAIN, CROSS) are frozen before any
-  physical result; no post-result loosening.
+  - Morison `|dc_d| ≤ .02`, `|dc_i| ≤ .02`.
 - Rule: all metrics for N = 3 consecutive whole paper cycles → settled.
 
 Cycle indexing (worked example, binding). Paper cycle k spans
@@ -159,20 +144,21 @@ quarter `s/T < 0.25`. Pair p_k compares cycle k−1 with cycle k.
   first cycle AFTER the last of those pairs. Earliest possible: pairs 3→4,
   4→5, 5→6 pass ⇒ K_case = 7.
 - Common window: K = max over compared cases of K_case; every case scores the
-  same three complete cycles K, K+1, K+2 and must stay settled continuously:
-  every pair p_j for j = its own K_case .. K+2 must pass (including the wait
-  between its local settle and the common K), otherwise that case is
-  FAIL-CLOSED.
+  same three complete cycles K, K+1, K+2 and must stay settled throughout:
+  pairs (K−1→K), (K→K+1), (K+1→K+2) must pass on every case, otherwise that
+  case is FAIL-CLOSED (an earlier local settle does not waive this).
 - Cap: cycle k ends at s/T = k + 0.25, so the last admissible scored cycle is
   11 (ends 11.25 T) and K ≤ 9. Examples: K = 7 → cycles 7–9 (ends 9.25 T);
   K = 9 → cycles 9–11 (ends 11.25 T); any K_case > 9 → FAIL-CLOSED.
-- Run length: cases run sequentially; each runs to the end of cycle 11
-  (11.25 T: 28,350 / 56,700 / 113,400 steps for D = 20 / 40 / 80; the 12 T
-  cap is 30,240 / 60,480 / 120,960) and keeps its scalar and profile data, then
-  the common K is chosen. An unscorable partial cycle 12 is not run.
-- Profiles: an empty common valid point set is FAIL-CLOSED; coverage, missing
-  points, bin time weights (equal per step) and the original fluid corner
-  weight sums are reported.
+- Every case is run to the end of cycle K + 2 (proportionally more native
+  steps on finer grids).
+- Discard: the initial quarter cycle plus the first 2 full paper cycles
+  (s/T < 2.25). Settling is searched from cycle 3; cap 12 periods from s = 0;
+  unsettled → FAIL-CLOSED (no literature, cross-resolution or VIV statistic).
+- Common window: K = the largest first-settled cycle index over all compared
+  cases; every case is scored on the same 3 cycles K .. K + 2 and must reach
+  them within the cap (proportionally more native steps on finer grids), or
+  the comparison is FAIL-CLOSED.
 - The earlier G4 timing check (100 coarse vs 400 diffusive-fine steps) was a
   matched short transient, not periodic evidence.
 
@@ -204,25 +190,6 @@ Superseded (earlier drafts): the velocity-jump start, guessed stations and
 phases, "C_M" for c_i, and the uncertainty wording; cycles 4–6, the 30-cycle
 cap and the earlier tolerance numbers.
 
-## 6b. Amendment: duration-only diagnostic (declared before the run)
-
-- Background: the D40 domain stage (20 D / 40 D / 60 D, 11.25 T each) failed
-  settling in every box (no window). That result stands.
-- One declared diagnostic: D = 40, 60 D box, from rest, through cycle 31
-  (157,500 steps = 31.25 T). Physics, settling tolerances, bands and
-  production unchanged; three proof pairs plus a continuously settled
-  K..K+2 with K ≤ 29.
-- Invocation: `run.py --D 40 --box 60 --cycles 31` then
-  `analyze.py --horizon 31`; the JSON records `horizon` and `horizon_mode`.
-  Without `--horizon` the analyzer is the original cap-11 scoring.
-- Failure at 31.25 T: report and stop; no automatic extension, matrix, D80
-  or VIV. A pass shows only that this box settles with time; it is not a
-  domain or literature pass.
-- The phase-locked high-frequency force component (residual of the mean
-  waveform after its Morison reconstruction, above 20/T) is reported per case
-  as unqualified; it needs separate qualification before any Morison or
-  literature statement.
-
 ## 6a. Driver unit tests (CPU, before any run)
 
 - Phase/cycle arithmetic: s ↔ (k, φ), cycle boundaries, bin membership,

research/duetsch_re100_kc5/analyze.py

@@ -1,224 +0,0 @@
-"""Ledger, settling, common window and Morison fits for Dütsch cases (numpy only).
-
-    python analyze.py CASE_DIR [CASE_DIR ...]   # cases of one comparison
-"""
-from __future__ import annotations
-
-import json
-import math
-import os
-import sys
-
-import numpy as np
-
-import protocol as p
-
-
-def load(path):
-    with open(os.path.join(path, "case.json")) as fh:
-        case = json.load(fh)
-    raw = np.load(os.path.join(path, "raw.npz"))
-    case.update(force=raw["force"], diag=raw["diag"], samples=raw["samples"], path=path)
-    return case
-
-
-def cycle_slices(case):
-    """{k: slice of step indices m whose interval [m, m+1] lies in paper cycle k}."""
-    T, q = case["T"], case["T"] // 4
-    return {k: slice(q + (k - 1) * T, q + k * T)
-            for k in range(1, case["cycles"] + 1) if q + k * T <= len(case["force"])}
-
-
-def coarse_wave(case, k, r=1):
-    """Normalized per-step Fx of cycle k, averaged over r native steps per sample."""
-    f = case["force"][cycle_slices(case)[k], 0] / case["F_ref"]
-    return f.reshape(-1, r).mean(axis=1)
-
-
-def profiles(case, k):
-    """{(phase, station, x2): (u, v)} of cycle k, time-averaged over the bin steps
-    (equal weights per step) from valid samples only."""
-    s = case["samples"]
-    acc = {}
-    for _, _, ph, xs, x2, u, v, _ in s[s[:, 1] == k]:
-        if not (math.isnan(u) or math.isnan(v)):
-            acc.setdefault((ph, xs, x2), []).append((u, v))
-    return {key: tuple(np.mean(vals, axis=0)) for key, vals in acc.items()}
-
-
-def coverage(case, k):
-    s = case["samples"]
-    rows = s[s[:, 1] == k]
-    missing = int(np.isnan(rows[:, 5]).sum())
-    return dict(samples=len(rows), missing=missing, bin_steps=len(set(rows[:, 0])),
-                min_fluid_weight=float(np.nanmin(np.where(rows[:, 7] > 0, rows[:, 7], np.nan)))
-                if len(rows) else math.nan)
-
-
-def common_points(dicts):
-    keys = set(dicts[0])
-    for d in dicts[1:]:
-        keys &= set(d)
-    return sorted(keys)
-
-
-def candidates(case):
-    """Per-step Fx of each reported force candidate (none is selected as an RHS)."""
-    d = case["diag"]
-    return {"cme": d[:, 0], "cme_wen": d[:, 0] + d[:, 2], "cme_node": d[:, 0] + d[:, 6] - d[:, 9]}
-
-
-def cycle_morison(case, ks, fx=None):
-    T = case["T"]
-    fx = case["force"][:, 0] if fx is None else fx
-    idx = np.concatenate([np.arange(len(fx))[cycle_slices(case)[k]] for k in ks])
-    b = np.array([p.interval_bases(int(m), T) for m in idx])
-    return p.morison_fit(fx[idx] / case["F_ref"], b[:, 0], b[:, 1])
-
-
-def pair_metrics(case, k):
-    a, b = coarse_wave(case, k - 1), coarse_wave(case, k)
-    pa, pb = profiles(case, k - 1), profiles(case, k)
-    common = common_points([pa, pb])  # empty set -> inf -> pair fails (FAIL-CLOSED)
-    prof = max((max(abs(pa[c][0] - pb[c][0]), abs(pa[c][1] - pb[c][1])) for c in common),
-               default=math.inf) / case["u_max"]
-    (cd0, ci0), (cd1, ci1) = cycle_morison(case, [k - 1]), cycle_morison(case, [k])
-    m = dict(wave=float(np.sqrt(np.mean((b - a) ** 2))),
-             scalar=float(max(abs(b.mean() - a.mean()), abs(np.sqrt(np.mean(b * b)) - np.sqrt(np.mean(a * a))),
-                              abs(np.abs(b).max() - np.abs(a).max()))),
-             profile=float(prof), morison=float(max(abs(cd1 - cd0), abs(ci1 - ci0))),
-             profile_points=len(common))
-    m["pass"] = all(m[key] <= tol for key, tol in p.SETTLE.items())
-    return m
-
-
-def hf_component(case, ks):
-    """Unqualified phase-locked high-frequency force component (reported, not checked):
-    mean CME waveform over cycles ks minus its two-basis Morison reconstruction,
-    split at 20/T; cycle differencing does not remove it."""
-    T = case["T"]
-    if not ks:
-        return dict(label="unqualified phase-locked HF force component, not checked", cycles=[])
-    mean = np.mean([coarse_wave(case, k) for k in ks], axis=0)
-    c_d, c_i = cycle_morison(case, ks)
-    q = T // 4
-    b = np.array([p.interval_bases(q + m, T) for m in range(T)])
-    resid = mean - (c_d * b[:, 0] - (np.pi ** 2 * c_i / p.KC) * b[:, 1])
-    spec = np.fft.rfft(resid)
-    high = spec.copy()
-    high[:21] = 0.0
-    energy = one_sided_energy(spec, T)
-    return dict(label="unqualified phase-locked HF force component, not checked",
-                cycles=list(ks), hf_rms_over_F_ref=float(np.sqrt(np.mean(np.fft.irfft(high, T) ** 2))),
-                hf_energy_fraction=float(energy[21:].sum() / max(energy.sum(), 1e-300)),
-                top_harmonics=[int(j) for j in np.argsort(energy[21:])[::-1][:5] + 21])
-
-
-def one_sided_energy(spec, n):
-    """Parseval-consistent energy per rfft bin: interior bins count twice, DC and
-    (for even n) Nyquist once."""
-    w = np.full(len(spec), 2.0)
-    w[0] = 1.0
-    if n % 2 == 0:
-        w[-1] = 1.0
-    return w * np.abs(spec) ** 2
-
-
-def ledger(case, horizon=p.LAST_SCORED):
-    ks = [k for k in sorted(cycle_slices(case)) if k <= horizon]
-    pairs = {k: pair_metrics(case, k) for k in ks if k - 1 > p.DISCARD}
-    return pairs, p.settled_start({k: v["pass"] for k, v in pairs.items()})
-
-
-def main(argv):
-    """argv: CASE_DIR ... [--reference CASE_DIR]; the reference defaults to the
-    largest D (resolution ladder); domain runs name the 60 D case explicitly."""
-    ref_path = None
-    if "--reference" in argv:
-        i = argv.index("--reference")
-        ref_path, argv = argv[i + 1], argv[:i] + argv[i + 2:]
-    horizon, explicit = p.LAST_SCORED, "--horizon" in argv  # default: original scoring
-    if explicit:
-        i = argv.index("--horizon")
-        horizon, argv = int(argv[i + 1]), argv[:i] + argv[i + 2:]
-    paths = argv if ref_path is None or ref_path in argv else argv + [ref_path]
-    cases = [load(path) for path in paths]
-    for c in cases if explicit else ():  # complete cycles 1..N; never pad or shift
-        if c["cycles"] < horizon or len(c["force"]) < round((0.25 + horizon) * c["T"]):
-            raise SystemExit(f"{c['path']}: run does not contain complete cycles 1..{horizon}")
-    ledgers = [ledger(c, horizon) for c in cases]
-    window = p.common_window([k for _, k in ledgers],
-                             [{k: v["pass"] for k, v in pairs.items()} for pairs, _ in ledgers],
-                             last_scored=horizon)
-    report = dict(window=window, horizon=horizon, cases=[],
-                  horizon_mode="default-cap" if horizon == p.LAST_SCORED else
-                  "extended duration diagnostic (declared amendment, D40/60D only)")
-    for case, (pairs, k_case) in zip(cases, ledgers):
-        limit = round((0.25 + horizon) * case["T"])  # events inside the horizon only
-        entry = {key: case[key] for key in ("path", "D", "T", "box", "n", "u_max", "nu", "tau", "A",
-                                            "steps", "cycles", "F_ref")}
-        entry.update(duration_T=case["steps"] / case["T"], K_case=k_case, pairs=pairs,
-                     coverage={k: coverage(case, k) for k in sorted(cycle_slices(case))
-                               if p.DISCARD < k <= horizon},
-                     no_donor=float(case["diag"][:limit, 13].sum()),
-                     eq_only=float(case["diag"][:limit, 12].sum()),
-                     hf_component=hf_component(case, window or tuple(
-                         [k for k in sorted(cycle_slices(case)) if p.DISCARD < k <= horizon][-3:])))
-        if window:  # residual drift inside the scored window (pairs K..K+2)
-            entry["window_drift"] = {key: max(pairs[k][key] for k in window)
-                                     for key in ("wave", "profile", "morison")}
-            entry["morison_obs"] = cycle_morison(case, window)
-            entry["morison"] = {name: cycle_morison(case, window, fx)
-                                for name, fx in candidates(case).items()}
-            entry["fy_mean"] = float(np.mean(np.concatenate(
-                [case["force"][cycle_slices(case)[k], 1] for k in window])) / case["F_ref"])
-        report["cases"].append(entry)
-    if window and len(cases) > 1:  # cross-resolution waveform against the finest case
-        fine = (next(c for c in cases if c["path"] == ref_path) if ref_path
-                else max(cases, key=lambda c: c["D"]))
-        report["reference"] = fine["path"]
-        ref = np.concatenate([coarse_wave(fine, k) for k in window])
-        for case, entry in zip(cases, report["cases"]):
-            r = fine["T"] // case["T"]  # 1 for equal-D domain runs
-            ref_c = ref.reshape(-1, r).mean(axis=1)
-            own = np.concatenate([coarse_wave(case, k) for k in window])
-            entry["wave_rms_vs_finest"] = float(np.sqrt(np.mean((own - ref_c) ** 2)))
-        for k in window:  # profiles on the common valid point set of all compared cases
-            per_case = [profiles(c, k) for c in cases]
-            keys = common_points(per_case)
-            report.setdefault("common_profile_points", {})[k] = len(keys)
-            if not keys:
-                report["verdict"] = "FAIL-CLOSED: empty common profile set"
-                continue
-            fine_p = per_case[cases.index(fine)]
-            for case, entry, prof in zip(cases, report["cases"], per_case):
-                d = max(max(abs(prof[key][0] - fine_p[key][0]), abs(prof[key][1] - fine_p[key][1]))
-                        for key in keys) / case["u_max"]
-                entry.setdefault("profile_max_vs_finest", {})[k] = float(d)
-        ref_entry = report["cases"][cases.index(fine)]
-        for entry in report["cases"]:  # conservative: difference + both drifts
-            if entry is ref_entry:
-                continue
-            drift = {key: entry["window_drift"][key] + ref_entry["window_drift"][key]
-                     for key in ("wave", "profile", "morison")}
-            dm = max(abs(entry["morison"]["cme"][0] - ref_entry["morison"]["cme"][0]),
-                     abs(entry["morison"]["cme"][1] - ref_entry["morison"]["cme"][1]))
-            entry["conservative"] = dict(
-                wave=entry["wave_rms_vs_finest"] + drift["wave"],
-                profile=max(entry.get("profile_max_vs_finest", {}).values(), default=math.inf)
-                + drift["profile"],
-                morison_cme=dm + drift["morison"])
-            bands = p.DOMAIN if ref_path else p.CROSS  # domain runs name the 60 D reference
-            c = entry["conservative"]
-            entry["pass_scope"] = (
-                ("wave+profile+morison_cme" if "morison" in bands else "wave+profile")
-                + ": conservative (|diff| + both window drifts) on the obs/CME Fx waveform and "
-                "simulation-to-simulation profiles; CME+Wen and CME+node reported, not checked; "
-                "no Fig. 8 / literature verdict")
-            entry["pass"] = (c["wave"] <= bands["wave"] and c["profile"] <= bands["profile"]
-                             and ("morison" not in bands or c["morison_cme"] <= bands["morison"]))
-    print(json.dumps(report, indent=1, default=float))
-
-
-if __name__ == "__main__":
-    main(sys.argv[1:])

research/duetsch_re100_kc5/capture.py

@@ -1,52 +0,0 @@
-"""Research-only field capture helpers for the Dütsch driver (numpy only).
-
-Snapshots are taken at the first step of the scored phase bins in declared
-cycles: a near-body crop at full resolution and a fluid-weighted block-mean
-downsample of the whole domain; plus scalar trends at every bin read.
-"""
-from __future__ import annotations
-
-import numpy as np
-
-CROP_X, CROP_Y, DOWN = 8.0, 6.0, 8      # crop half-widths in D; downsample factor
-STRIP = 2.0                             # inlet/outlet strip width in D
-
-
-def crop_box(xm, ym, D, n, body_extent):
-    """Inclusive integer crop (x0, x1, y0, y1); fail closed if it leaves the domain
-    interior or does not contain the body's swept extent (x_lo, x_hi, y_lo, y_hi)."""
-    x0, x1 = int(round(xm - CROP_X * D)), int(round(xm + CROP_X * D))
-    y0, y1 = int(round(ym - CROP_Y * D)), int(round(ym + CROP_Y * D))
-    if x0 < 1 or y0 < 1 or x1 > n - 2 or y1 > n - 2:
-        raise ValueError("capture crop leaves the domain interior")
-    bx0, bx1, by0, by1 = body_extent
-    if not (x0 <= bx0 and bx1 <= x1 and y0 <= by0 and by1 <= y1):
-        raise ValueError("capture crop clips the body's swept extent")
-    return x0, x1, y0, y1
-
-
-def downsample(field, fluid, f=DOWN):
-    """Fluid-weighted f x f block means; NaN where a block has no fluid cell."""
-    ny, nx = field.shape
-    h, w = ny // f, nx // f
-    v = np.where(fluid, field, 0.0)[:h * f, :w * f].reshape(h, f, w, f).sum(axis=(1, 3))
-    c = fluid[:h * f, :w * f].reshape(h, f, w, f).sum(axis=(1, 3))
-    with np.errstate(invalid="ignore", divide="ignore"):
-        return np.where(c > 0, v / np.maximum(c, 1), np.nan).astype(np.float32)
-
-
-def fluid_mean(field, fluid, xs=slice(None), ys=slice(None)):
-    m = fluid[ys, xs]
-    return float(field[ys, xs][m].mean()) if m.any() else float("nan")
-
-
-def trend_row(s, k, ph, rho, ux, uy, fluid, D, n, crop):
-    """(s, cycle, phase, <rho>, <ux>, <uy>, <rho>_crop, inlet <rho>, <ux>, outlet <rho>, <ux>)."""
-    w = max(1, int(round(STRIP * D)))
-    x0, x1, y0, y1 = crop
-    cx, cy = slice(x0, x1 + 1), slice(y0, y1 + 1)
-    inlet, outlet, inner = slice(1, 1 + w), slice(n - 1 - w, n - 1), slice(1, n - 1)
-    return (s, k, ph, fluid_mean(rho, fluid), fluid_mean(ux, fluid), fluid_mean(uy, fluid),
-            fluid_mean(rho, fluid, cx, cy),
-            fluid_mean(rho, fluid, inlet, inner), fluid_mean(ux, fluid, inlet, inner),
-            fluid_mean(rho, fluid, outlet, inner), fluid_mean(ux, fluid, outlet, inner))

research/duetsch_re100_kc5/protocol.py

@@ -1,129 +0,0 @@
-"""Pure arithmetic of the Dütsch Re 100 / KC 5 protocol (numpy only, no CUDA).
-
-docs/validation/duetsch_re100_kc5_protocol.md is the contract: elapsed time s
-in native steps from the right turning point, paper phase
-phi = (360 s / T - 90) mod 360, paper cycle k >= 1 spans
-s/T in [0.25 + (k-1), 0.25 + k); cycle 0 is the startup quarter.
-"""
-from __future__ import annotations
-
-import math
-
-import numpy as np
-
-KC, RE = 5.0, 100.0
-STATIONS = (-0.6, 0.0, 0.6, 1.2)          # lab-fixed x/D from the mean position
-PHASES = (180.0, 210.0, 330.0)            # degrees; bins [p, p + 1)
-X2 = tuple(np.round(np.arange(-1.5, 1.5001, 0.05), 10))  # declared x2/D points
-DISCARD, LAST_SCORED, N_PAIRS = 2, 11, 3   # cycles 0..2 discarded, cap 12 T
-# A-priori bands. Settling is one fifth of the tightest comparison band of the
-# same observable, and comparisons add both cases' residual window drift
-# (conservative: |difference| + drift_a + drift_b <= band).
-DOMAIN = dict(wave=.01, morison=.02, profile=.02)  # vs the 60 D reference at D = 40
-CROSS = dict(wave=.02, profile=.02)                # D40 vs D80
-SETTLE = dict(wave=.002, scalar=.002, profile=.004, morison=.004)
-
-
-def amplitude(D):
-    return KC * D / (2.0 * math.pi)
-
-
-def position(s, T, A):
-    """Centre position relative to the mean: x = A cos(2 pi s / T), evaluated from the
-    native step modulo T so every cycle repeats the same geometry bit for bit."""
-    return A * math.cos(2.0 * math.pi * (int(s) % T) / T)
-
-
-def phase(s, T):
-    """Paper phase in degrees from the native step modulo T (the unwrapped cycle
-    index is kept separately by ``cycle``)."""
-    return (360.0 * (int(s) % T) / T - 90.0) % 360.0
-
-
-def cycle(s, T):
-    return 0 if s < 0.25 * T else int(math.floor(s / T - 0.25)) + 1
-
-
-def bin_phase(s, T):
-    """Scored phase whose 1-degree bin [p, p + 1) contains time s, else None."""
-    phi = phase(s, T)
-    for p in PHASES:
-        if p <= phi < p + 1.0:
-            return p
-    return None
-
-
-def interval_bases(n, T, points=16):
-    """Single-step interval averages over [n, n + 1] of cos phi |cos phi| and sin phi
-    (phi in radians = 2 pi s / T - pi / 2), Gauss-Legendre quadrature."""
-    x, w = np.polynomial.legendre.leggauss(points)
-    s = int(n) % T + 0.5 * (x + 1.0)  # periodic: identical bases in every cycle
-    phi = 2.0 * np.pi * s / T - 0.5 * np.pi
-    c = np.cos(phi)
-    return 0.5 * float(np.sum(w * c * np.abs(c))), 0.5 * float(np.sum(w * np.sin(phi)))
-
-
-def morison_fit(f_norm, b1, b2):
-    """Least squares of f_norm = c_d b1 - (pi^2 c_i / KC) b2 (paper phase, no shift)."""
-    a = np.column_stack([b1, -(math.pi ** 2 / KC) * np.asarray(b2)])
-    (c_d, c_i), *_ = np.linalg.lstsq(a, np.asarray(f_norm), rcond=None)
-    return float(c_d), float(c_i)
-
-
-def inside_circle(x, y, centre, radius):
-    return (x - centre[0]) ** 2 + (y - centre[1]) ** 2 < radius * radius
-
-
-def sample(field, fluid, x, y, centre, radius):
-    """Profile sample: MISSING if the exact point is inside the body (checked before
-    any corner interpolation), else ``bilinear_fluid``."""
-    if inside_circle(x, y, centre, radius):
-        return None, (), 0.0
-    return bilinear_fluid(field, fluid, x, y)
-
-
-def bilinear_fluid(field, fluid, x, y):
-    """Bilinear value at (x, y) in cell coordinates from fluid cells only.
-
-    Returns (value, ((ix, iy, weight), ...), fluid_weight_sum) with weights
-    renormalized over the fluid corners and the original sum of their bilinear
-    weights, or (None, (), 0.0) when no corner is fluid (MISSING, never zero).
-    """
-    x0, y0 = int(math.floor(x)), int(math.floor(y))
-    fx, fy = x - x0, y - y0
-    corners = ((x0, y0, (1 - fx) * (1 - fy)), (x0 + 1, y0, fx * (1 - fy)),
-               (x0, y0 + 1, (1 - fx) * fy), (x0 + 1, y0 + 1, fx * fy))
-    used = [(i, j, w) for i, j, w in corners if w > 0.0 and fluid[j, i]]
-    total = sum(w for _, _, w in used)
-    if total == 0.0:
-        return None, (), 0.0
-    used = tuple((i, j, w / total) for i, j, w in used)
-    return float(sum(w * field[j, i] for i, j, w in used)), used, float(total)
-
-
-def settled_start(pairs_pass):
-    """K_case from {k: pass of pair (k-1 -> k)}: first cycle after three consecutive
-    passing pairs among cycles > DISCARD, or None."""
-    run = 0
-    for k in sorted(pairs_pass):
-        if k - 1 <= DISCARD:
-            continue
-        run = run + 1 if pairs_pass[k] else 0
-        if run == N_PAIRS:
-            return k + 1
-    return None
-
-
-def common_window(k_cases, pairs_pass_by_case, last_scored=LAST_SCORED):
-    """Common scored cycles K..K+2 with K = max K_case, or None (FAIL-CLOSED) when a
-    case never settled, K + 2 exceeds the cap, or a case fails any pair p_j for
-    j = its own K_case .. K + 2 (continuous settling, including the wait for K)."""
-    if any(k is None for k in k_cases):
-        return None
-    K = max(k_cases)
-    if K + 2 > last_scored:
-        return None
-    for k_case, pairs in zip(k_cases, pairs_pass_by_case):
-        if not all(pairs.get(j, False) for j in range(k_case, K + 3)):
-            return None
-    return K, K + 1, K + 2

research/duetsch_re100_kc5/run.py

@@ -1,133 +0,0 @@
-"""One-shot Dütsch Re 100 / KC 5 prescribed in-line oscillation driver.
-
-Protocol: docs/validation/duetsch_re100_kc5_protocol.md. Native dense channel
-(zero-velocity inlet, pressure-NEQ outlet, bounce-back y walls), MRT, lab frame.
-Writes raw per-step force/diagnostics and phase-bin profile samples to --out
-(outside the checkout); analysis is in analyze.py. No production change.
-
-    python run.py --D 40 --box 40 --cycles 11 --out /path/outside/checkout
-"""
-from __future__ import annotations
-
-import argparse
-import json
-import math
-import os
-import tempfile
-
-import numpy as np
-
-import protocol as p
-
-
-TEMPLATE = os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "..", "src",
-                        "CelerisLab", "configs", "config_lbm.json")
-
-
-def config(D, box, T):
-    """Shipped config_lbm.json (full parser schema) with only the case fields
-    overridden: grid, still fluid (velocity 0), nu from Re, MRT, uniform inlet
-    profile, bounce-back y walls. Inlet scheme and outlet mode stay the
-    template's (zou_he_local, neq_extrap) and are recorded in case.json."""
-    u = p.KC * D / T
-    n = int(round(box * D))
-    with open(TEMPLATE) as fh:
-        cfg = json.load(fh)
-    cfg["grid"].update(lattice_model="D2Q9", nx=n, ny=n, nz=1)
-    cfg["physics"].update(data_type="FP32", viscosity=u * D / p.RE, velocity=0.0, rho=1.0)
-    cfg["method"].update(collision="MRT", streaming="double_buffer", store_precision="FP32",
-                         ddf_shifting=False, y_wall_bc="bounce_back")
-    cfg["method"]["les"]["enabled"] = False
-    cfg["method"]["inlet"]["profile"] = "uniform"
-    return cfg, u, n
-
-
-def main():
-    ap = argparse.ArgumentParser()
-    ap.add_argument("--D", type=int, required=True)
-    ap.add_argument("--box", type=float, default=40.0, help="box side in D")
-    ap.add_argument("--T", type=int, default=None, help="period in steps (multiple of 360)")
-    ap.add_argument("--cycles", type=int, default=p.LAST_SCORED, help="run to the end of this cycle")
-    ap.add_argument("--out", required=True)
-    ap.add_argument("--dump-fields", default="",
-                    help="research capture: comma-separated paper cycles, e.g. 11,21,30,31")
-    args = ap.parse_args()
-    D, T = args.D, args.T or 126 * args.D
-    if T % 360:
-        raise SystemExit("T must be a multiple of 360 steps")
-    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")
-    cfg, u_max, n = config(D, args.box, T)
-    A, R = p.amplitude(D), 0.5 * D
-    xm = ym = 0.5 * (n - 1)
-    os.makedirs(args.out, exist_ok=True)
-
-    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)
-    try:
-        body = sim.add_body("circle", (xm + A, ym), R,
-                            translation=Translation((-2.0 * A, 0.0, 0.0, 0.0), 1.01 * u_max))
-        sim.initialize()
-        dump = {int(c) for c in args.dump_fields.split(",") if c}
-        if dump:
-            import capture
-            if max(dump) > args.cycles or min(dump) <= p.DISCARD:
-                raise SystemExit("--dump-fields cycles must lie in the scored range of this run")
-            crop = capture.crop_box(xm, ym, D, n, (xm - A - R, xm + A + R, ym - R, ym + R))
-            trend = []
-        steps = int(round((0.25 + args.cycles) * T))
-        force = np.zeros((steps, 2)); diag = []; samples = []
-        for m in range(steps):
-            x0, x1 = p.position(m, T, A), p.position(m + 1, T, A)
-            sim.set_body(body, offset=(max(x1 - A, -2.0 * A), 0.0), vx=x1 - x0)
-            sim.run(1)
-            force[m] = sim.read_force(body)[:2]
-            diag.append(sim.moving_diagnostics()[0])
-            s = m + 1  # state after the step
-            ph = p.bin_phase(s, T)
-            if ph is not None and p.cycle(s, T) > p.DISCARD:
-                fields = sim.get_macroscopic()
-                fluid = (sim.get_flags().reshape(n, n) & 1) != 0
-                ux, uy = (np.asarray(fields[k]).reshape(n, n) for k in ("ux", "uy"))
-                if dump:  # research capture only; the default path is unchanged
-                    rho = np.asarray(fields["rho"]).reshape(n, n)
-                    trend.append(capture.trend_row(s, p.cycle(s, T), ph, rho, ux, uy, fluid, D, n, crop))
-                    if p.cycle(s, T) in dump and p.phase(s, T) == ph:  # first step of the bin
-                        x0c, x1c, y0c, y1c = crop
-                        cy_, cx_ = slice(y0c, y1c + 1), slice(x0c, x1c + 1)
-                        np.savez_compressed(
-                            os.path.join(args.out, f"fields_c{p.cycle(s, T)}_p{int(ph)}.npz"),
-                            s=s, cycle=p.cycle(s, T), phase=ph, crop=np.asarray(crop),
-                            force=force[m], centre=np.asarray((xm + p.position(s, T, A), ym)),
-                            ux=ux[cy_, cx_].astype(np.float32), uy=uy[cy_, cx_].astype(np.float32),
-                            rho=rho[cy_, cx_].astype(np.float32), fluid=fluid[cy_, cx_].astype(np.uint8),
-                            down_ux=capture.downsample(ux, fluid), down_uy=capture.downsample(uy, fluid),
-                            down_rho=capture.downsample(rho, fluid))
-                centre = (xm + p.position(s, T, A), ym)  # analytic centre after the step
-                for xs in p.STATIONS:
-                    for x2 in p.X2:
-                        x, y = xm + xs * D, ym + x2 * D
-                        u, wu, wsum = p.sample(ux, fluid, x, y, centre, R)
-                        v, _, _ = p.sample(uy, fluid, x, y, centre, R)
-                        samples.append((s, p.cycle(s, T), ph, xs, x2,
-                                        math.nan if u is None else u, math.nan if v is None else v,
-                                        wsum, json.dumps(wu)))
-        if dump:
-            np.save(os.path.join(args.out, "trend.npy"), np.asarray(trend, float))
-        np.savez(os.path.join(args.out, "raw.npz"), force=force, diag=np.asarray(diag),
-                 samples=np.asarray([r[:8] for r in samples], float))
-        with open(os.path.join(args.out, "case.json"), "w") as fh:
-            json.dump(dict(D=D, T=T, box=args.box, n=n, u_max=u_max, nu=cfg["physics"]["viscosity"],
-                           tau=3 * cfg["physics"]["viscosity"] + 0.5, A=A, steps=steps,
-                           cycles=args.cycles, F_ref=0.5 * u_max ** 2 * D,
-                           inlet=cfg["method"]["inlet"], outlet=cfg["method"]["outlet"],
-                           weights=[r[8] for r in samples]), fh)
-    finally:
-        sim.close()
-
-
-if __name__ == "__main__":
-    main()

tests/unit/test_duetsch_protocol.py

@@ -1,191 +0,0 @@
-"""CPU tests of the Dütsch protocol arithmetic (docs/validation/duetsch_re100_kc5_protocol.md §6a)."""
-import importlib.util
-import math
-from pathlib import Path
-
-import numpy as np
-import pytest
-
-_spec = importlib.util.spec_from_file_location(
-    "duetsch_protocol", Path(__file__).resolve().parents[2] / "research/duetsch_re100_kc5/protocol.py")
-p = importlib.util.module_from_spec(_spec)
-_spec.loader.exec_module(p)
-
-T = 2520
-
-
-def test_phase_and_cycle_bookkeeping():
-    assert p.phase(0, T) == 270.0 and p.cycle(0, T) == 0          # +A, at rest, phi = -90
-    assert p.phase(T // 4, T) == 0.0 and p.cycle(T // 4, T) == 1  # mean position, u = -U_max
-    assert p.position(0, T, 1.0) == 1.0
-    assert p.cycle(T // 4 - 1, T) == 0 and p.cycle(T // 4 + T - 1, T) == 1
-    assert p.cycle(T // 4 + T, T) == 2
-    assert p.phase(T // 6, T) == pytest.approx(330.0) and p.cycle(T // 6, T) == 0  # startup quarter
-    assert p.phase(7 * T // 6, T) == pytest.approx(330.0) and p.cycle(7 * T // 6, T) == 1
-    for ph, frac in ((180.0, .75), (210.0, 5 / 6)):
-        s = round(frac * T) + T
-        assert p.phase(s, T) == pytest.approx(ph) and p.bin_phase(s, T) == ph
-    assert p.bin_phase(round(.75 * T) + T - 1, T) is None
-    assert sum(p.bin_phase(s, T) == 180.0 for s in range(T, 2 * T)) == T // 360
-
-
-def test_settling_and_common_window():
-    passes = {k: True for k in range(3, 12)}
-    assert p.settled_start(passes) == 7                          # pairs 3->4, 4->5, 5->6
-    late = {**passes, 5: False}
-    assert p.settled_start(late) == 9                            # pairs 6->7, 7->8, 8->9
-    assert p.common_window([7, 9], [passes, passes]) == (9, 10, 11)
-    assert p.common_window([7, 10], [passes, passes]) is None    # K + 2 = 12 beyond the cap
-    assert p.common_window([7, None], [passes, passes]) is None
-    assert p.common_window([7, 7], [passes, {**passes, 8: False}]) is None  # must stay settled
-    # a case settled early must keep passing while waiting for the common K
-    assert p.common_window([7, 9], [{**passes, 8: False}, passes]) is None
-    long = {k: True for k in range(3, 32)}
-    assert p.common_window([10], [long]) is None                 # default cap 11
-    assert p.common_window([29], [long], last_scored=31) == (29, 30, 31)
-    assert p.common_window([30], [long], last_scored=31) is None
-
-
-def test_solid_aware_bilinear():
-    field = np.arange(16.0).reshape(4, 4)
-    fluid = np.ones((4, 4), bool)
-    value, used, total = p.bilinear_fluid(field, fluid, 1.25, 1.5)
-    assert total == pytest.approx(1.0) and value == pytest.approx(0.75 * 0.5 * (5 + 9) + 0.25 * 0.5 * (6 + 10))
-    fluid[1, 1] = False                                           # one solid corner
-    value, used, total = p.bilinear_fluid(field, fluid, 1.25, 1.5)
-    assert total == pytest.approx(1.0 - 0.75 * 0.5)
-    assert all(not (i == 1 and j == 1) for i, j, _ in used)
-    assert sum(w for *_, w in used) == pytest.approx(1.0)
-    fluid[:] = False
-    assert p.bilinear_fluid(field, fluid, 1.25, 1.5) == (None, (), 0.0)
-
-
-def test_interval_bases_and_morison_fit_recover_known_coefficients():
-    b = np.array([p.interval_bases(n, T) for n in range(T // 4, T // 4 + 3 * T)])
-    phi = 2 * np.pi * (np.arange(T // 4, T // 4 + 3 * T) + .5) / T - np.pi / 2
-    assert np.allclose(b[:, 1], np.sin(phi), atol=1e-5)          # midpoint ~ interval mean
-    f = 2.09 * b[:, 0] - (math.pi ** 2 * 1.45 / p.KC) * b[:, 1]
-    assert p.morison_fit(f, b[:, 0], b[:, 1]) == pytest.approx((2.09, 1.45), abs=1e-9)
-
-
-def test_geometry_repeats_exactly_across_cycles():
-    A = p.amplitude(20)
-    for m in (0, 1, T // 4, T // 2 - 1, 7 * T // 12):
-        for c in (1, 5, 11):
-            assert p.position(m + c * T, T, A) == p.position(m, T, A)
-            assert p.phase(m + c * T, T) == p.phase(m, T)
-            v0 = p.position(m + 1, T, A) - p.position(m, T, A)
-            assert p.position(m + 1 + c * T, T, A) - p.position(m + c * T, T, A) == v0
-    assert p.cycle(T // 4 + 5 * T, T) == 6  # unwrapped index stays separate
-    assert p.interval_bases(123 + 7 * T, T) == p.interval_bases(123, T)
-
-
-def test_exact_point_inside_body_is_missing_even_with_fluid_corners():
-    field = np.zeros((6, 6))
-    fluid = np.ones((6, 6), bool)
-    fluid[3, 3] = False                         # only one corner solid
-    centre, radius = (3.0, 3.0), 0.3
-    value, used, total = p.sample(field, fluid, 3.2, 3.1, centre, radius)   # inside the circle
-    assert value is None and used == () and total == 0.0
-    value, used, total = p.sample(field, fluid, 3.35, 3.0, centre, radius)  # just outside
-    assert value is not None and total > 0.0
-
-
-def test_generated_config_passes_the_real_loader(tmp_path):
-    import json
-    import sys
-    sys.path.insert(0, str(Path(__file__).resolve().parents[2] / "research/duetsch_re100_kc5"))
-    try:
-        import run
-    finally:
-        sys.path.pop(0)
-    from CelerisLab.config import load_lbm_config
-    cfg_dict, u_max, n = run.config(20, 10.0, 2520)
-    path = tmp_path / "config_lbm.json"
-    path.write_text(json.dumps(cfg_dict))
-    cfg = load_lbm_config(str(path))
-    assert (cfg.nx, cfg.ny) == (n, n) == (200, 200)
-    assert cfg.collision == "MRT" and cfg.streaming == "double_buffer" and not cfg.ddf_shifting
-    assert cfg.velocity == 0.0 and cfg.viscosity == pytest.approx(u_max * 20 / 100)
-    assert cfg.y_wall_bc == "bounce_back"
-
-
-def _analyze():
-    import sys
-    sys.path.insert(0, str(Path(__file__).resolve().parents[2] / "research/duetsch_re100_kc5"))
-    try:
-        import analyze
-    finally:
-        sys.path.pop(0)
-    return analyze
-
-
-def test_hf_energy_is_parseval_weighted():
-    a = _analyze()
-    n = 360
-    t = np.arange(n)
-    x = 1.0 + np.cos(2 * np.pi * 5 * t / n) + 2.0 * np.cos(2 * np.pi * 30 * t / n) + 0.5 * np.cos(np.pi * t)
-    e = a.one_sided_energy(np.fft.rfft(x), n) / n
-    assert e.sum() == pytest.approx(np.sum(x ** 2))               # Parseval
-    expected = (2.0 ** 2 / 2 + 0.5 ** 2) / np.mean(x ** 2)       # bins >= 21 (30 and Nyquist)
-    assert e[21:].sum() / e.sum() == pytest.approx(expected)
-
-
-def test_horizon_excludes_future_coverage_and_events(monkeypatch):
-    a = _analyze()
-    T, horizon = 360, 4
-    steps = round((0.25 + 6) * T)                                 # cycles 1..6 recorded
-    diag = np.zeros((steps, 15))
-    diag[round((0.25 + horizon) * T):, 13] = 1.0                  # synthetic future no_donor
-    rows = [(s, p.cycle(s, T), 180.0, 0.0, 0.0, 0.1, 0.0, 1.0) for s in range(steps)
-            if p.bin_phase(s, T) is not None and p.cycle(s, T) > p.DISCARD]
-    case = dict(path="c", D=20, T=T, box=10, n=200, u_max=.04, nu=.008, tau=.524, A=1.0,
-                steps=steps, cycles=6, F_ref=1.0, force=np.zeros((steps, 2)), diag=diag,
-                samples=np.asarray(rows, float))
-    monkeypatch.setattr(a, "load", lambda path: case)
-    out = []
-    monkeypatch.setattr("builtins.print", lambda text: out.append(text))
-    a.main(["c", "--horizon", str(horizon)])
-    import json
-    entry = json.loads(out[0])["cases"][0]
-    assert sorted(int(k) for k in entry["coverage"]) == [3, 4]   # nothing from cycles 5, 6
-    assert entry["no_donor"] == 0.0
-
-
-
-def _capture():
-    spec = importlib.util.spec_from_file_location(
-        "duetsch_capture", Path(__file__).resolve().parents[2] / "research/duetsch_re100_kc5/capture.py")
-    mod = importlib.util.module_from_spec(spec)
-    spec.loader.exec_module(mod)
-    return mod
-
-
-def test_capture_crop_downsample_and_trends():
-    c = _capture()
-    n, D = 2400, 40
-    xm = ym = 0.5 * (n - 1)
-    A, R = p.amplitude(D), 20.0
-    box = c.crop_box(xm, ym, D, n, (xm - A - R, xm + A + R, ym - R, ym + R))
-    assert box == (880, 1520, 960, 1440) and box[1] - box[0] + 1 == 641 and box[3] - box[2] + 1 == 481
-    with pytest.raises(ValueError):
-        c.crop_box(xm, ym, D, 600, (0, 1, 0, 1))                    # leaves the domain
-    with pytest.raises(ValueError):
-        c.crop_box(xm, ym, D, n, (xm - 9 * D, xm, ym, ym))           # clips the swept body
-    field = np.arange(64.0).reshape(8, 8)
-    fluid = np.ones((8, 8), bool)
-    fluid[:4, :4] = False
-    down = c.downsample(field, fluid, 4)
-    assert np.isnan(down[0, 0]) and down[1, 1] == field[4:, 4:].mean()
-    rho = np.ones((n, n)); rho[:, 1:81] = 1.01
-    row = c.trend_row(10, 3, 180.0, rho, np.zeros((n, n)), np.zeros((n, n)),
-                      np.ones((n, n), bool), D, n, box)
-    assert row[7] == pytest.approx(1.01) and row[9] == pytest.approx(1.0)
-
-
-def test_dump_steps_are_the_first_step_of_each_bin():
-    T = 5040
-    for k in (11, 21, 30, 31):
-        starts = [s for s in range(T // 4 + (k - 1) * T, T // 4 + k * T)
-                  if p.bin_phase(s, T) is not None and p.phase(s, T) == p.bin_phase(s, T)]
-        assert len(starts) == 3 and all(p.cycle(s, T) == k for s in starts)