From b31a6141b3230501e16e95572d30b0bba5c4b229 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 30 Oct 2026 10:58:05 +0000 Subject: [PATCH] =?UTF-8?q?research:=20D=C3=BCtsch=20driver=20=E2=80=94=20?= =?UTF-8?q?continuous=20settling,=20domain=20criteria,=20profile=20coverag?= =?UTF-8?q?e?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit common_window requires every pair from each case's own K_case to K+2; a-priori domain criteria against the 60 D reference (separate from Table 1); bilinear sampling returns the original fluid corner weight sum; analysis time-averages bin samples, reports coverage/missing points and fails closed on an empty common profile set; cases run to the end of cycle 11. Co-Authored-By: Claude Opus 5.5 --- docs/validation/duetsch_re100_kc5_protocol.md | 35 ++++++++++------- research/duetsch_re100_kc5/analyze.py | 38 ++++++++++++++++++- research/duetsch_re100_kc5/protocol.py | 17 +++++---- research/duetsch_re100_kc5/run.py | 10 ++--- tests/unit/test_duetsch_protocol.py | 11 ++++-- 5 files changed, 79 insertions(+), 32 deletions(-) diff --git a/docs/validation/duetsch_re100_kc5_protocol.md b/docs/validation/duetsch_re100_kc5_protocol.md index 9e90c73..85f1851 100644 --- a/docs/validation/duetsch_re100_kc5_protocol.md +++ b/docs/validation/duetsch_re100_kc5_protocol.md @@ -87,8 +87,16 @@ 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 (default 40 D if, against 60 D, `RMS(dFx)/F_ref` and |dc_d|, |dc_i| - stay inside the section 6 bands). + 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. - Memory/runtime: measured by a smoke run before each new size; measurements inform scheduling only, never tolerances. @@ -144,21 +152,20 @@ 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 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). + 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. - 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. -- 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. +- 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. - The earlier G4 timing check (100 coarse vs 400 diffusive-fine steps) was a matched short transient, not periodic evidence. diff --git a/research/duetsch_re100_kc5/analyze.py b/research/duetsch_re100_kc5/analyze.py index 68bf284..4a5075a 100644 --- a/research/duetsch_re100_kc5/analyze.py +++ b/research/duetsch_re100_kc5/analyze.py @@ -36,9 +36,30 @@ def coarse_wave(case, k, r=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] - return {(ph, xs, x2): (u, v) for _, _, ph, xs, x2, u, v in rows} + 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 cycle_morison(case, ks): @@ -51,7 +72,7 @@ def cycle_morison(case, ks): 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 = [key for key in pa if key in pb and not any(map(math.isnan, pa[key] + pb[key]))] + 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]) @@ -80,6 +101,7 @@ def main(paths): 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 k > p.DISCARD}, no_donor=float(case["diag"][:, 13].sum()), eq_only=float(case["diag"][:, 12].sum())) if window: entry["morison"] = cycle_morison(case, window) @@ -94,6 +116,18 @@ def main(paths): 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) print(json.dumps(report, indent=1, default=float)) diff --git a/research/duetsch_re100_kc5/protocol.py b/research/duetsch_re100_kc5/protocol.py index bb287b2..df65e35 100644 --- a/research/duetsch_re100_kc5/protocol.py +++ b/research/duetsch_re100_kc5/protocol.py @@ -17,6 +17,7 @@ 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 SETTLE = dict(wave=.02, scalar=.01, profile=.03, morison=.02) +DOMAIN = dict(wave=.01, morison=.02, profile=.02) # a-priori, vs the 60 D reference def amplitude(D): @@ -65,8 +66,9 @@ def morison_fit(f_norm, b1, b2): def bilinear_fluid(field, fluid, x, y): """Bilinear value at (x, y) in cell coordinates from fluid cells only. - Returns (value, ((ix, iy, weight), ...)) with weights renormalized over the - fluid corners, or (None, ()) when no corner is fluid (MISSING, never zero). + 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 @@ -75,9 +77,9 @@ def bilinear_fluid(field, fluid, x, y): 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, () + 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 + return float(sum(w * field[j, i] for i, j, w in used)), used, float(total) def settled_start(pairs_pass): @@ -95,13 +97,14 @@ def settled_start(pairs_pass): def common_window(k_cases, pairs_pass_by_case): """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 a window pair.""" + 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 pairs in pairs_pass_by_case: - if not all(pairs.get(k, False) for k in (K, K + 1, K + 2)): + 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 diff --git a/research/duetsch_re100_kc5/run.py b/research/duetsch_re100_kc5/run.py index a0f2747..f50cce9 100644 --- a/research/duetsch_re100_kc5/run.py +++ b/research/duetsch_re100_kc5/run.py @@ -73,18 +73,18 @@ def main(): for xs in p.STATIONS: for x2 in p.X2: x, y = xm + xs * D, ym + x2 * D - u, wu = p.bilinear_fluid(ux, fluid, x, y) - v, _ = p.bilinear_fluid(uy, fluid, x, y) + u, wu, wsum = p.bilinear_fluid(ux, fluid, x, y) + v, _, _ = p.bilinear_fluid(uy, fluid, x, y) 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, - json.dumps(wu))) + wsum, json.dumps(wu))) np.savez(os.path.join(args.out, "raw.npz"), force=force, diag=np.asarray(diag), - samples=np.asarray([r[:7] for r in samples], float)) + 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, - weights=[r[7] for r in samples]), fh) + weights=[r[8] for r in samples]), fh) finally: sim.close() diff --git a/tests/unit/test_duetsch_protocol.py b/tests/unit/test_duetsch_protocol.py index 3b367d9..2de5e34 100644 --- a/tests/unit/test_duetsch_protocol.py +++ b/tests/unit/test_duetsch_protocol.py @@ -38,19 +38,22 @@ def test_settling_and_common_window(): 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, dict(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], [dict(passes, **{8: False}), passes]) is None def test_solid_aware_bilinear(): field = np.arange(16.0).reshape(4, 4) fluid = np.ones((4, 4), bool) - value, used = p.bilinear_fluid(field, fluid, 1.25, 1.5) - assert value == pytest.approx(0.75 * 0.5 * (5 + 9) + 0.25 * 0.5 * (6 + 10)) + 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 = p.bilinear_fluid(field, fluid, 1.25, 1.5) + 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, ()) + assert p.bilinear_fluid(field, fluid, 1.25, 1.5) == (None, (), 0.0) def test_interval_bases_and_morison_fit_recover_known_coefficients():