From 6f95c34c0bf88408416655fc3baeee02d1455096 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 30 Oct 2026 10:54:59 +0000 Subject: [PATCH] =?UTF-8?q?research:=20one-shot=20D=C3=BCtsch=20Re100/KC5?= =?UTF-8?q?=20prescribed-oscillation=20driver=20(no=20production=20change)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 --- docs/validation/duetsch_re100_kc5_protocol.md | 209 ++++++++++++++++++ research/duetsch_re100_kc5/analyze.py | 101 +++++++++ research/duetsch_re100_kc5/protocol.py | 107 +++++++++ research/duetsch_re100_kc5/run.py | 93 ++++++++ tests/unit/test_duetsch_protocol.py | 61 +++++ 5 files changed, 571 insertions(+) create mode 100644 docs/validation/duetsch_re100_kc5_protocol.md create mode 100644 research/duetsch_re100_kc5/analyze.py create mode 100644 research/duetsch_re100_kc5/protocol.py create mode 100644 research/duetsch_re100_kc5/run.py create mode 100644 tests/unit/test_duetsch_protocol.py diff --git a/docs/validation/duetsch_re100_kc5_protocol.md b/docs/validation/duetsch_re100_kc5_protocol.md new file mode 100644 index 0000000..9e90c73 --- /dev/null +++ b/docs/validation/duetsch_re100_kc5_protocol.md @@ -0,0 +1,209 @@ +# Prescribed in-line oscillation: Dütsch Re 100 / KC 5 — PROTOCOL (review draft) + +Production under test: `claude/moving-bodies-s3a` at `653fd997` (unchanged). +Primary source: Dütsch, Durst, Becker & Lienhart, J. Fluid Mech. 360, 249–271 +(1998), doi:10.1017/S002211209800860X. Source facts below were verified from +the paper (Slack C0C78UFKGRM/1793335557.460509, TRON 1793356317–1793356724, +page numbers given). Everything marked ENGINEERING is our a-priori choice, +fixed before any benchmark result; it is not a paper uncertainty. + +## 0. Source facts (with provenance) + +- Motion (p.251): `x_c = -A sin(2 pi f t)`; paper phase phi = 0 is the mean + position moving left at `-U_max`. KC = 5 gives `A/D = 5/(2 pi) = .7957747`. +- Start (p.254): the computation starts a quarter cycle before phi = 0, at + `x = +A` with zero velocity (phi = -90°). +- Experiment (pp.252, 254–256): D = 10 mm, L = .65 m, tank 1.2 × .7 × .7 m, + rho = 998.2 kg/m³, mu = 1.00818 mPa·s, Re = 100, KC = 5; LDA with 1° phase + bins and 200 samples per node and degree; no quantified experimental error + recovered. +- Profiles (Fig. 8, p.260): LAB-FIXED stations x1 = -6, 0, 6, 12 mm, i.e. + x/D = -.6, 0, .6, 1.2, at phi = 180°, 210°, 330°; u1 and u2 versus x2. +- Original numerics (§2): cylinder-fixed accelerated frame (`x_i = x~_i + + x_i,s`, outer velocity `-u_c`, inertial source `-rho du_c/dt`), O-grid outer + diameter 120 D, illustrated 192 circumferential nodes, 720 steps/cycle; tank + influence neglected. +- Force fit (verified): force ON THE BODY per unit span + `F1 = -0.5 rho D c_d xdot|xdot| - (pi/4) rho D² c_i xddot` with `x = -A sin phi`; + normalized by `0.5 rho D U_max²`: `c_d cos phi |cos phi| - (pi² c_i / KC) sin phi`. + `c_i` is the ADDED-MASS coefficient; total inertia `c_m = c_i + 1` is a + different convention and is not used. +- Table 1 computed Morison coefficients (c_d, c_i) for 96 / 192 / 384 perimeter + nodes at dt/T = 1/720: (2.15, 1.42), (2.10, 1.45), (2.09, 1.45). These are + fitted computed coefficients, not a measured force and not an error bar. The + paper states the coarse direct force is closer to the finest direct force + than the Morison fit, but gives no waveform-error norm. + +## 1. Kinematics, time levels and phase bookkeeping + +- Our run is in the LAB frame: the circle moves, the fluid is initially at rest + (rho = 1, u = 0); no fictitious inertial source is added. Lab-fixed profile + stations are compared directly (no frame transform needed). +- Geometry reference: the circle's base centre is the INITIAL position, the + right turning point `x_mean + A`; the declared offset box is x in [-2A, 0], + y = 0. The body is never teleported: offset 0 at s = 0. +- Elapsed time s = n dt from the start. Analytic position relative to the mean + `x(s) = A cos(2 pi s / T)`, velocity `u(s) = -U_max sin(2 pi s / T)`; paper + phase `phi = (360 s / T - 90) mod 360`. +- Step n moves from s_n to s_{n+1}: target offset `x(s_{n+1}) - A`; wall + velocity = the single-step interval average `(x(s_{n+1}) - x(s_n)) / dt` + (analytic, supplied, not differenced from stored offsets). The per-step + force of `run(1)` / `read_force(normalize=True)` is attributed to that same + interval. +- Paper cycle k = 1, 2, ... spans `s in [T/4 + (k-1) T, T/4 + k T)`; the + initial quarter cycle (s < T/4) is cycle 0 and never scored. Every logged + sample carries (k, phi). Example: phi = 330° of cycle 1 is at s = 7T/6. +- T is chosen as a multiple of 360 native steps, so phi = 180°, 210°, 330° fall + exactly on step boundaries and 1° bins hold whole steps. + +## 2. Resolution ladder (uniform grid, fixed Mach), ENGINEERING + +| D (cells) | T (steps) | U_max = 5 D / T | nu_lat = U_max D / 100 | tau (3 nu + .5) | +|---|---|---|---|---| +| 20 | 2,520 | .039683 | .0079365 | .5238 | +| 40 | 5,040 | .039683 | .015873 | .5476 | +| 80 | 10,080 | .039683 | .031746 | .5952 | + +- MRT collision. Mach fixed (Ma ≈ .0687); a separate Mach check at D = 40 + with T = 10,080 (U_max halved, tau .5238) is labelled MACH and not part of the + ladder. +- Per-step displacement ≤ U_max ≈ .04 cells (bound .5). + +## 3. Domain and boundaries (finite-domain reproduction), ENGINEERING + +The paper's 120 D O-grid is not reproduced; our rectangular box is a +reproduction choice. The native dense channel (zero-velocity inlet, +pressure-NEQ outlet, bounce-back y walls) is a finite-domain approximation, +not infinite still fluid. The G4 x-periodic / y-slip override is a separately +labelled sensitivity family and never mixed into the primary ladder. The +runner owner confirms the zero-velocity inlet is supported as configured +before any run. + +Order (closed): +1. Smoke + ledger at D = 40, native BC, 40 D × 40 D (finite box). +2. Domain sensitivity before any literature statement: same BC, D = 40, + L ∈ {20 D, 40 D, 60 D}; record the change in Fx waveform, c_d, c_i and the + scored profiles. +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). +- Memory/runtime: measured by a smoke run before each new size; measurements + inform scheduling only, never tolerances. + +## 4. Observables (raw, per step, preserved) + +- `obs` Fx, Fy (CME) and all `MovingTables.DIAG` columns per step: CME, + CME + Wen and CME + cover − uncover reported separately; crossings, eq_only, + no_donor (no_donor must stay 0). +- Velocity u1, u2 on lab-fixed lines x/D = -.6, 0, .6, 1.2 (relative to the + mean cylinder position) versus x2, at the scored phases. Lab frame + throughout; if any body-frame intermediate is used, position and velocity + are both transformed. +- Sampling cost: raw force and diagnostics counters are read every step + (small). Flags and velocity fields are read only at the steps inside the + three prescribed 1° phase bins (φ ∈ [180, 181), [210, 211), [330, 331), + T/360 steps each) of the cycles needed for settling and scoring, plus + sparse declared audits; never every step. +- Point set: fixed x2/D sample points per station, identical on every grid + and declared before running; per point the interpolation weights are + recorded; the scored set is the common subset valid on all compared grids + at that phase, with its coverage recorded; point sets never change silently. +- Interpolation: bilinear on cell centres using only cells that are fluid at + that step (flags read back via `get_flags()`); weights renormalized over the + fluid corners; a sample with no fluid corner, or inside the circle at that + step, is MISSING and logged, never set to zero. One-degree bins keep valid + fluid weights. The numerical average need not mimic the 200 LDA samples. + +## 5. Periodic-regime gate (every case, before any statistic), ENGINEERING + +Ledger rows, same schema for every case: dx, dt (in D and T units); U_max, nu, +Re, KC (lattice and nondimensional); T (native steps, nondimensional); +duration (native steps, s/T, paper cycles); discard; accepted window; cycle +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: + - 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| ≤ .02`, `|dc_i| ≤ .02`. +- Rule: all metrics for N = 3 consecutive whole paper cycles → settled. + +Cycle indexing (worked example, binding). Paper cycle k spans +`s/T in [0.25 + (k-1), 0.25 + k)`, i.e. φ0 → φ0; cycle 0 is the startup +quarter `s/T < 0.25`. Pair p_k compares cycle k−1 with cycle k. + +- Discarded: cycles 0, 1, 2 (s/T < 2.25). Pairs are formed only between + cycles ≥ 3, so the first pair is p_4 (3 → 4). +- Settling proof for a case: three consecutive pairs all pass; K_case = the + 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). +- 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. +- The earlier G4 timing check (100 coarse vs 400 diffusive-fine steps) was a + matched short transient, not periodic evidence. + +## 6. Metrics and per-observable criteria (ENGINEERING, not paper uncertainty) + +Three criteria stay separate; the G4 same-grid 2% target is not reused. + +- Morison fit: least squares of the normalized per-step force over the + accepted window against the two basis functions, each averaged over its + single step interval, fixed paper phase, no free phase shift; c_d and + added-mass c_i per force candidate. +- Cross-resolution: every scored observable's error against its reference must + decrease monotonically D20 → D40 → D80 (necessary, not sufficient). +- Against Table 1 (computed fit, not measured force): at D ≥ 40, + `|c_d - 2.09| ≤ .15` and `|c_i - 1.45| ≤ .15`, and monotone. +- Force waveform against its own Morison reconstruction: residual reported + only, no pass/fail (the paper reports fit error above its coarse numerical + error and gives no percentage). +- Fig. 8 LDA profiles: QUALITATIVE ONLY. The original figure could not be + accessed or rendered; no digitized ordinates and no measurement uncertainty + exist. After a digitization with figure number, tool and uncertainty note: + `RMS(|du|, |dv|) / U_max ≤ .15` at D ≥ 40 and monotone; LDA error stays + UNQUANTIFIED and is stated with every comparison. +- Fy sanity (not literature): `|mean Fy| / F_ref ≤ .05` in the accepted window. +- No force candidate is selected for an RHS by these results; each is + reported. + +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. + +## 6a. Driver unit tests (CPU, before any run) + +- Phase/cycle arithmetic: s ↔ (k, φ), cycle boundaries, bin membership, + window/cap selection for K = 7, 8, 9, 10. +- Solid-aware interpolation: all-fluid, partial-solid (renormalized weights), + no-fluid (MISSING) and inside-circle cases. +- Interval-averaged Morison bases and the two-parameter fit on a synthetic + waveform with known c_d, c_i. + +## 7. Order and scope + +Protocol review → one-shot Dütsch driver (analytic `set_body` + `run(1)`, +cycle ledger, lab-fixed fluid-only station sampler, two-basis Morison fit; no +reusable runner until the first smoke shows shared helpers pay off; no +production change) → smoke (memory/runtime) → domain sensitivity → +ladder → physical-validation verdict → only then single-cylinder VIV. No AMR, +no pinball, no manuscript work. diff --git a/research/duetsch_re100_kc5/analyze.py b/research/duetsch_re100_kc5/analyze.py new file mode 100644 index 0000000..68bf284 --- /dev/null +++ b/research/duetsch_re100_kc5/analyze.py @@ -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:]) diff --git a/research/duetsch_re100_kc5/protocol.py b/research/duetsch_re100_kc5/protocol.py new file mode 100644 index 0000000..bb287b2 --- /dev/null +++ b/research/duetsch_re100_kc5/protocol.py @@ -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 diff --git a/research/duetsch_re100_kc5/run.py b/research/duetsch_re100_kc5/run.py new file mode 100644 index 0000000..a0f2747 --- /dev/null +++ b/research/duetsch_re100_kc5/run.py @@ -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() diff --git a/tests/unit/test_duetsch_protocol.py b/tests/unit/test_duetsch_protocol.py new file mode 100644 index 0000000..3b367d9 --- /dev/null +++ b/tests/unit/test_duetsch_protocol.py @@ -0,0 +1,61 @@ +"""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 = dict(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, dict(passes, **{8: False})]) is None # must stay settled + + +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)) + fluid[1, 1] = False # one solid corner + value, used = p.bilinear_fluid(field, fluid, 1.25, 1.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, ()) + + +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)