diff --git a/docs/validation/dfg_re20.md b/docs/validation/dfg_re20.md index 255475b..4662f26 100644 --- a/docs/validation/dfg_re20.md +++ b/docs/validation/dfg_re20.md @@ -6,7 +6,7 @@ inflow, initial state, layout and observation are those of the pinned Re100 runn (`run_refined_dfg.py`, `re100_helpers.py`, both unchanged); only the lattice viscosity is set for Re 20 (root tau 0.56, fine 0.62, uniform 0.62, lattice speed 0.01). -This is a comparison with the official intervals, not a full benchmark +Label: "compared against the official reference bands", not a full benchmark certification: the outlet is pressure-NEQ (`neq_extrap`), not do-nothing, and the inlet is local Zou-He. @@ -22,10 +22,13 @@ inlet is local Zou-He. ## Bounds and steady rule (frozen before any run) - End time t* = 60 (30 s physical): 240,000 refined root steps, 480,000 uniform steps. -- One sample every t* = 0.5; trailing window t* = 10 (5 s, 20 samples). -- Steady: in the window, max - min of each of Cd, Cl, p_diff is at most one tenth of - its official interval width (0.002, 0.00006, 0.00004). Absolute ranges, so the - small Cl is never used as a divisor. +- One sample every t* = 1 (0.5 s) at the end of a synchronized root step; trailing + window t* = 10 (last 5 s, 11 samples). +- Steady: in the window, both max - min and |least-squares slope x window| of each + of Cd, Cl, p_diff are at most 1 % of its official interval width (2e-4, 6e-6, + 4e-6). Absolute limits, so the small Cl is never a divisor. Never relaxed after a run. +- Force sign: fluid on cylinder, x along the inflow, y up from the bottom wall; Cd > 0 + and the official Cl > 0 (`cd_positive` is reported). - Not steady at t* = 60: status `NOT_STEADY`, no interval verdict. Health limits (Mach <= 0.2, rho in [0.5, 1.5]) stop the run with `UNHEALTHY`. - Steady: `WITHIN_INTERVALS` or `OUTSIDE_INTERVALS`, with the per-quantity verdict. @@ -34,10 +37,10 @@ inlet is local Zou-He. ``` python tests/validation/run_dfg_re20.py --output --plan -CUDA_VISIBLE_DEVICES= python tests/validation/run_dfg_re20.py --case both \ +CUDA_VISIBLE_DEVICES= python tests/validation/run_dfg_re20.py --case refined \ --collision SRT --output --max-wall-seconds ``` -`--plan` writes `plan.json` only. A GPU run writes `plan.json` and `result.json` +Run `--case uniform` only if the refined result is outside a band. `--plan` writes `plan.json` only. A GPU run writes `plan.json` and `result.json` (per case: lattice parameters, every sample, window ranges, status, verdict). Host contract tests: `tests/unit/test_dfg_re20.py`. diff --git a/tests/unit/test_dfg_re20.py b/tests/unit/test_dfg_re20.py index e90423b..6a3de8c 100644 --- a/tests/unit/test_dfg_re20.py +++ b/tests/unit/test_dfg_re20.py @@ -21,8 +21,8 @@ class Re20Contract(unittest.TestCase): root, uniform = v.lattice_parameters(40), v.lattice_parameters(80) self.assertAlmostEqual(root["tau"], 0.56, places=12) self.assertAlmostEqual(uniform["tau"], 0.62, places=12) - self.assertEqual((root["steps"], root["interval"]), (240000, 2000)) - self.assertEqual((uniform["steps"], uniform["interval"]), (480000, 4000)) + self.assertEqual((root["steps"], root["interval"]), (240000, 4000)) + self.assertEqual((uniform["steps"], uniform["interval"]), (480000, 8000)) self.assertAlmostEqual(v.plan()["refined_fine_tau"], 0.62, places=12) def test_config_only_changes_viscosity(self): @@ -47,17 +47,22 @@ class Re20Contract(unittest.TestCase): def test_steady_and_verdict(self): inside = dict(cd=5.58, cl=0.0106, p_diff=0.1174) out = v.summarize(rows(lambda i: inside)) - self.assertEqual((out["status"], out["steady"]), ("WITHIN_INTERVALS", True)) + self.assertEqual((out["status"], out["steady"], out["window_samples"]), + ("WITHIN_INTERVALS", True, 11)) + self.assertTrue(out["cd_positive"]) self.assertTrue(all(out["verdict"].values())) low = v.summarize(rows(lambda i: dict(inside, cl=0.0100))) self.assertEqual(low["status"], "OUTSIDE_INTERVALS") self.assertFalse(low["verdict"]["cl"]) def test_drift_gives_no_verdict(self): - # Cl drifts by more than a tenth of its interval width inside the window. - out = v.summarize(rows(lambda i: dict(cd=5.58, cl=0.0106 + 1e-5 * i, p_diff=0.1174))) + # Cl drifts by more than 1 % of its interval width (6e-6) inside the window. + out = v.summarize(rows(lambda i: dict(cd=5.58, cl=0.0106 + 1e-6 * i, p_diff=0.1174))) self.assertEqual((out["status"], out["verdict"]), ("NOT_STEADY", None)) self.assertGreater(out["ranges"]["cl"], out["steady_limits"]["cl"]) + self.assertGreater(out["drifts"]["cl"], out["steady_limits"]["cl"]) + small = v.summarize(rows(lambda i: dict(cd=5.58, cl=0.0106 + 2e-6 * (i % 2), p_diff=0.1174))) + self.assertTrue(small["steady"]) # 2e-6 alternation, no trend: within both limits def test_incomplete_run_gives_no_verdict(self): out = v.summarize(rows(lambda i: dict(cd=5.58, cl=0.0106, p_diff=0.1174), tstar_end=30)) diff --git a/tests/validation/run_dfg_re20.py b/tests/validation/run_dfg_re20.py index 8cd08a0..414a43b 100644 --- a/tests/validation/run_dfg_re20.py +++ b/tests/validation/run_dfg_re20.py @@ -11,16 +11,24 @@ Cd = 2 Fx / (rho U^2 D) and Cl = 2 Fy / (rho U^2 D) are scale-free, and ``pressure_probes`` returns dp / (rho U^2). The official p_diff = p(0.15, 0.2) - p(0.25, 0.2) uses U = 0.2 m/s, rho = 1, so p_diff = dp * 0.04. Both are reported; only p_diff is compared with the official interval. Forces use the existing sign -convention (positive Cd along the inflow, positive Cl towards +y). +convention: F is the force of the fluid on the cylinder, i.e. the refined raw +exchange impulse (window ledger, ``total.impulse``) or the uniform ``read_force``, +with x along the inflow and y up from the bottom wall. Cd must be positive, and the +official Cl is positive (centre 0.20 below the centreline 0.205), so a flipped sign +shows up as Cl outside its interval. -Sampling. One sample every SAMPLE_TSTAR (t* = U t / D): Cd and Cl are the mean -over that interval (window exchange ledger refined, per-batch forces uniform); -p_diff is the end-of-interval value at the exact surface points. +Sampling. One sample every SAMPLE_TSTAR (t* = U t / D; 0.5 s), taken at the end of a +completed root step with all levels synchronized; the step index is recorded. Cd +and Cl are the mean over the sample interval (smooths step-parity/acoustic +oscillation; equal to the instantaneous value at steady state); p_diff is the +end-of-interval value at the exact surface points. -Steady rule (frozen before any run). Over the trailing WINDOW_TSTAR, the range -(max - min) of each of Cd, Cl and p_diff must not exceed one tenth of the width of -its official interval. Absolute ranges avoid dividing by the small Cl. If the rule -fails at TSTAR_MAX the run reports "not steady" and gives no interval verdict. +Steady rule (frozen before any run). Over the trailing WINDOW_TSTAR (last 5 s, 11 +samples), both the range (max - min) and |least-squares slope * window| of each of +Cd, Cl and p_diff must not exceed 1 % of the width of its official interval. These +are absolute limits, so the small Cl is never a divisor, and the slope term catches +slow monotone drift. If the rule fails at TSTAR_MAX the run reports "not steady" and +gives no interval verdict; the limits are never relaxed after a run. """ import argparse @@ -41,11 +49,11 @@ import run_refined_dfg as rr # noqa: E402 RE = 20.0 U_PHYS = 0.2 # m/s mean inflow; rho = 1 TSTAR_MAX = 60.0 # 30 s physical -SAMPLE_TSTAR = 0.5 +SAMPLE_TSTAR = 1.0 # 0.5 s WINDOW_TSTAR = 10.0 # trailing 5 s # Schaefer & Turek (1996), 2D-1 lower/upper bounds. BANDS = dict(cd=(5.57, 5.59), cl=(0.0104, 0.0110), p_diff=(0.1172, 0.1176)) -STEADY_RANGE = {k: (hi - lo) / 10 for k, (lo, hi) in BANDS.items()} +STEADY_LIMIT = {k: (hi - lo) / 100 for k, (lo, hi) in BANDS.items()} DEVIATIONS = ("pressure-NEQ outlet (neq_extrap), not do-nothing", "local Zou-He inlet") _re100_config = rr.config @@ -78,11 +86,14 @@ def summarize(rows): if not rows or rows[-1]["tstar"] < TSTAR_MAX - 1e-9: return dict(status="INCOMPLETE", steady=False, verdict=None) window = [r for r in rows if r["tstar"] >= TSTAR_MAX - WINDOW_TSTAR - 1e-9] + t = np.array([r["tstar"] for r in window]) ranges = {k: max(r[k] for r in window) - min(r[k] for r in window) for k in BANDS} - steady = all(ranges[k] <= STEADY_RANGE[k] for k in BANDS) + drifts = {k: abs(np.polyfit(t, [r[k] for r in window], 1)[0]) * WINDOW_TSTAR for k in BANDS} + steady = all(ranges[k] <= STEADY_LIMIT[k] and drifts[k] <= STEADY_LIMIT[k] for k in BANDS) final = {k: rows[-1][k] for k in BANDS} - result = dict(window_samples=len(window), ranges=ranges, steady_limits=STEADY_RANGE, - steady=steady, final=final, bands=BANDS, deviations=DEVIATIONS) + result = dict(window_samples=len(window), ranges=ranges, drifts=drifts, + steady_limits=STEADY_LIMIT, steady=steady, final=final, bands=BANDS, + deviations=DEVIATIONS, cd_positive=final["cd"] > 0) if not steady: return dict(result, status="NOT_STEADY", verdict=None) inside = {k: BANDS[k][0] <= final[k] <= BANDS[k][1] for k in BANDS} @@ -171,16 +182,16 @@ def run_case(a, h, case): def plan(collision="SRT"): return dict(re=RE, collision=collision, tstar_max=TSTAR_MAX, sample_tstar=SAMPLE_TSTAR, - window_tstar=WINDOW_TSTAR, bands=BANDS, steady_limits=STEADY_RANGE, + window_tstar=WINDOW_TSTAR, bands=BANDS, steady_limits=STEADY_LIMIT, deviations=DEVIATIONS, pressure_scale=U_PHYS**2, refined_root=lattice_parameters(40), refined_fine_tau=0.5 + 3 * 2 * rr.SPEED * 40 / RE, uniform=lattice_parameters(80), - label="compared with official reference intervals; not a full benchmark certification") + label="compared against the official reference bands (non-canonical neq_extrap outlet, zou_he_local inlet); not a full benchmark certification") def main(): parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) - parser.add_argument("--case", choices=["both", "refined", "uniform"], default="both") + parser.add_argument("--case", choices=["refined", "uniform", "both"], default="refined") parser.add_argument("--collision", choices=["SRT", "MRT"], default="SRT") parser.add_argument("--source", type=Path, default=rr.DEFAULT_SOURCE) parser.add_argument("--output", type=Path, required=True)