research: one-shot Dütsch Re100/KC5 prescribed-oscillation driver (no production change)
protocol.py holds the pure arithmetic (phase/cycle bookkeeping from the right
turning point, 1-degree bins, interval-averaged Morison bases and fit with
added-mass c_i, solid-aware bilinear sampling, settling and common-window
rules). run.py drives the existing Simulation API (set_body + run(1) +
read_force + moving_diagnostics; macroscopic fields and flags only inside the
three phase bins) and writes raw data outside the checkout. analyze.py builds
the cycle ledger, K_case, the common window, Morison fits, Fy sanity and the
cross-resolution waveform difference. CPU tests cover the §6a arithmetic.
Carries the protocol doc from claude/duetsch-protocol 8c1f3c4d.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
ClaudeandClaude Opus 5.5 committed
1 parent
653fd997a6
commit
6f95c34c0b
5 files changed
+571
No files matched your search
@@ -0,0 +1,101 @@
|
||||
"""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):
|
||||
s = case["samples"]
|
||||
rows = s[s[:, 1] == k]
|
||||
return {(ph, xs, x2): (u, v) for _, _, ph, xs, x2, u, v in rows}
|
||||
|
||||
|
||||
def cycle_morison(case, ks):
|
||||
T = case["T"]
|
||||
idx = np.concatenate([np.arange(len(case["force"]))[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])
|
||||
|
||||
|
||||
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]))]
|
||||
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 ledger(case):
|
||||
ks = sorted(cycle_slices(case))
|
||||
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(paths):
|
||||
cases = [load(path) for path in paths]
|
||||
ledgers = [ledger(c) 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])
|
||||
report = dict(window=window, cases=[])
|
||||
for case, (pairs, k_case) in zip(cases, ledgers):
|
||||
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,
|
||||
no_donor=float(case["diag"][:, 13].sum()), eq_only=float(case["diag"][:, 12].sum()))
|
||||
if window:
|
||||
entry["morison"] = cycle_morison(case, window)
|
||||
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"])
|
||||
ref = np.concatenate([coarse_wave(fine, k) for k in window])
|
||||
for case, entry in zip(cases, report["cases"]):
|
||||
r = fine["T"] // case["T"]
|
||||
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)))
|
||||
print(json.dumps(report, indent=1, default=float))
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main(sys.argv[1:])
|
||||
@@ -0,0 +1,107 @@
|
||||
"""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
|
||||
SETTLE = dict(wave=.02, scalar=.01, profile=.03, morison=.02)
|
||||
|
||||
|
||||
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)."""
|
||||
return A * math.cos(2.0 * math.pi * s / T)
|
||||
|
||||
|
||||
def phase(s, T):
|
||||
return (360.0 * s / 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 = n + 0.5 * (x + 1.0)
|
||||
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 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).
|
||||
"""
|
||||
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, ()
|
||||
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
|
||||
|
||||
|
||||
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):
|
||||
"""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."""
|
||||
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)):
|
||||
return None
|
||||
return K, K + 1, K + 2
|
||||
@@ -0,0 +1,93 @@
|
||||
"""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
|
||||
|
||||
|
||||
def config(D, box, T):
|
||||
u = p.KC * D / T
|
||||
n = int(round(box * D))
|
||||
return {"grid": {"lattice_model": "D2Q9", "nx": n, "ny": n, "nz": 1},
|
||||
"physics": {"data_type": "FP32", "viscosity": u * D / p.RE, "velocity": 0.0, "rho": 1.0},
|
||||
"method": {"collision": "MRT", "streaming": "double_buffer", "store_precision": "FP32",
|
||||
"ddf_shifting": False, "y_wall_bc": "bounce_back"}}, 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)
|
||||
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()
|
||||
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"))
|
||||
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)
|
||||
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)))
|
||||
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))
|
||||
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)
|
||||
finally:
|
||||
sim.close()
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
main()
|
||||
Reference in new issue
Block a user