From 73c2f507f41201e42059821764d66e5732ede473 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 30 Oct 2026 10:58:48 +0000 Subject: [PATCH] =?UTF-8?q?research:=20D=C3=BCtsch=20analysis=20=E2=80=94?= =?UTF-8?q?=20Morison=20fit=20per=20force=20candidate,=20explicit=20--refe?= =?UTF-8?q?rence?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Opus 5.5 --- research/duetsch_re100_kc5/analyze.py | 32 +++++++++++++++++++++------ 1 file changed, 25 insertions(+), 7 deletions(-) diff --git a/research/duetsch_re100_kc5/analyze.py b/research/duetsch_re100_kc5/analyze.py index 4a5075a..40f55a3 100644 --- a/research/duetsch_re100_kc5/analyze.py +++ b/research/duetsch_re100_kc5/analyze.py @@ -62,11 +62,18 @@ def common_points(dicts): return sorted(keys) -def cycle_morison(case, ks): +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"] - idx = np.concatenate([np.arange(len(case["force"]))[cycle_slices(case)[k]] for k in ks]) + 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(case["force"][idx, 0] / case["F_ref"], b[:, 0], b[:, 1]) + return p.morison_fit(fx[idx] / case["F_ref"], b[:, 0], b[:, 1]) def pair_metrics(case, k): @@ -91,7 +98,14 @@ def ledger(case): return pairs, p.settled_start({k: v["pass"] for k, v in pairs.items()}) -def main(paths): +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:] + paths = argv if ref_path is None or ref_path in argv else argv + [ref_path] cases = [load(path) for path in paths] ledgers = [ledger(c) for c in cases] window = p.common_window([k for _, k in ledgers], @@ -104,15 +118,19 @@ def main(paths): 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) + 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 = max(cases, key=lambda c: c["D"]) + 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"] + 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)))