research: Dütsch driver — continuous settling, domain criteria, profile coverage

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 <noreply@anthropic.com>

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))
 
 

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

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()