diff --git a/docs/validation/dfg_re20.md b/docs/validation/dfg_re20.md new file mode 100644 index 0000000..255475b --- /dev/null +++ b/docs/validation/dfg_re20.md @@ -0,0 +1,43 @@ +# DFG 2D-1 (Re 20) against the official reference intervals + +`tests/validation/run_dfg_re20.py` runs the steady DFG 2D-1 case on the refined +channel root (D40 root, D80 wide fine patch) and on the uniform D80 grid. Geometry, +inflow, initial state, layout and observation are those of the pinned Re100 runner +(`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 +certification: the outlet is pressure-NEQ (`neq_extrap`), not do-nothing, and the +inlet is local Zou-He. + +## Quantities + +- Cd = 2 Fx / (rho U^2 D), Cl = 2 Fy / (rho U^2 D), interval means over each sample. +- p_diff = p(0.15, 0.2) - p(0.25, 0.2), linear extrapolation to the exact surface + points (`re100_helpers.pressure_probes`). The helper returns dp / (rho U^2) in + scaled units; p_diff = dp * 0.04 for U = 0.2 m/s, rho = 1. Both are reported. +- Official intervals (Schaefer & Turek 1996): Cd [5.57, 5.59], Cl [0.0104, 0.0110], + p_diff [0.1172, 0.1176]. + +## 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. +- 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. + +## Run + +``` +python tests/validation/run_dfg_re20.py --output --plan +CUDA_VISIBLE_DEVICES= python tests/validation/run_dfg_re20.py --case both \ + --collision SRT --output --max-wall-seconds +``` + +`--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 new file mode 100644 index 0000000..e90423b --- /dev/null +++ b/tests/unit/test_dfg_re20.py @@ -0,0 +1,68 @@ +"""CPU-only contract tests for run_dfg_re20. Not CUDA or DFG acceptance evidence.""" + +import sys +from pathlib import Path +import unittest + +sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "validation")) +import run_dfg_re20 as v # noqa: E402 +import run_refined_dfg as rr # noqa: E402 + +h, _ = rr.helpers(rr.DEFAULT_SOURCE) + + +def rows(values, tstar_end=v.TSTAR_MAX): + n = int(round(tstar_end / v.SAMPLE_TSTAR)) + return [dict(tstar=(i + 1) * v.SAMPLE_TSTAR, **values(i)) for i in range(n)] + + +class Re20Contract(unittest.TestCase): + def test_lattice_parameters_and_bounds(self): + 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.assertAlmostEqual(v.plan()["refined_fine_tau"], 0.62, places=12) + + def test_config_only_changes_viscosity(self): + base = v._re100_config(h, 40, "SRT", "none") + cfg, body, p = v.config(h, 40, "SRT") + self.assertAlmostEqual(cfg["physics"]["viscosity"], 0.01 * 40 / 20, places=15) + self.assertEqual((p["re"], p["viscosity"]), (20.0, cfg["physics"]["viscosity"])) + base[0]["physics"]["viscosity"] = cfg["physics"]["viscosity"] + self.assertEqual(base[0], cfg) + self.assertEqual(base[1], body) + with self.assertRaisesRegex(ValueError, "no Jets"): + v.config(h, 40, "MRT", "zero") + + def test_units(self): + s = v.sample([2.0, -0.5], 4, 40, 3.0) + den = rr.SPEED**2 * 40 / 2 + self.assertAlmostEqual(s["cd"], 2.0 / 4 / den) + self.assertAlmostEqual(s["cl"], -0.5 / 4 / den) + self.assertAlmostEqual(s["p_diff"], 3.0 * 0.04) + self.assertEqual(s["dp_scaled"], 3.0) + + 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.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))) + self.assertEqual((out["status"], out["verdict"]), ("NOT_STEADY", None)) + self.assertGreater(out["ranges"]["cl"], out["steady_limits"]["cl"]) + + 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)) + self.assertEqual((out["status"], out["verdict"]), ("INCOMPLETE", None)) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/validation/run_dfg_re20.py b/tests/validation/run_dfg_re20.py new file mode 100644 index 0000000..8cd08a0 --- /dev/null +++ b/tests/validation/run_dfg_re20.py @@ -0,0 +1,206 @@ +"""Steady DFG 2D-1 (Re 20): refined and uniform runs against the official reference intervals. + +Not a full benchmark certification. Declared deviations from the official setup: +pressure-NEQ outlet (``neq_extrap``) instead of do-nothing, local Zou-He inlet. +Geometry, inflow profile, initial state, refined layout and observation reuse the +pinned Re100 runner (``run_refined_dfg``) and ``re100_helpers`` unchanged; only the +lattice viscosity is set for Re 20 (nu = U_lattice * D / 20 at the root level). + +Units. The helpers use scaled DFG units with mean inflow speed 1, so +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). + +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. + +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. +""" + +import argparse +import json +import math +import os +from pathlib import Path +import sys +import time +from tempfile import TemporaryDirectory + +import numpy as np + +HERE = Path(__file__).resolve().parent +sys.path.insert(0, str(HERE)) +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 +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()} +DEVIATIONS = ("pressure-NEQ outlet (neq_extrap), not do-nothing", "local Zou-He inlet") + +_re100_config = rr.config + + +def config(h, diameter, collision="SRT", jet_scenario="none"): + """The pinned Re100 configuration with the lattice viscosity set for Re 20.""" + if jet_scenario != "none": + raise ValueError("DFG 2D-1 runs have no Jets") + cfg, body, p = _re100_config(h, diameter, collision, jet_scenario) + nu = p["mean_speed"] * p["diameter"] / RE + cfg["physics"]["viscosity"] = nu + return cfg, body, dict(p, viscosity=nu, re=RE) + + +def lattice_parameters(diameter, speed=rr.SPEED): + nu = speed * diameter / RE + return dict(viscosity=nu, tau=0.5 + 3 * nu, steps=math.ceil(TSTAR_MAX * diameter / speed), + interval=round(SAMPLE_TSTAR * diameter / speed)) + + +def sample(impulse, n, diameter, dp): + den = rr.SPEED**2 * diameter / 2 + return dict(cd=float(impulse[0] / n / den), cl=float(impulse[1] / n / den), + dp_scaled=float(dp), p_diff=float(dp) * U_PHYS**2) + + +def summarize(rows): + """Steady rule over the trailing window, then the interval verdict (or none).""" + 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] + 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) + 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) + 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} + return dict(result, status="WITHIN_INTERVALS" if all(inside.values()) else "OUTSIDE_INTERVALS", + verdict=inside) + + +def healthy(row): + return row["max_mach"] <= 0.2 and row["rho_min"] >= 0.5 and row["rho_max"] <= 1.5 + + +def run_case(a, h, case): + from contextlib import ExitStack + from CelerisLab.cuda.context import CudaContext + from CelerisLab.lbm.descriptors import FLUID + from CelerisLab.lbm.refinement_stepper import FixedPatchStepper + from CelerisLab.simulation import Simulation + + refined = case == "refined" + cfgdict, bodydict, p = config(h, 40 if refined else 80, a.collision) + lp = lattice_parameters(p["diameter"]) + rows = [] + with ExitStack() as stack: + if refined: + cfg, pair, body, host = rr.host_geometry(a.source, h, a.collision) + context = stack.enter_context(CudaContext(0)) + backend = stack.enter_context(FixedPatchStepper( + cfg, rr.PLACEMENTS["wide"], context=context, coarse_boundary="channel", + base_block=rr.BASE_BLOCK, fine_block=rr.FINE_BLOCK, fine_objects=(body,), + observe_raw_exchange=True, observe_window_exchange=True)) + from CelerisLab.lbm.channel import channel_masks + from CelerisLab.lbm.refinement import ACTIVE + from CelerisLab.lbm.descriptors import SOLID + masks = [channel_masks(pair.coarse_shape)[0] & ~pair.covered, + (pair.fine_kind == ACTIVE) + & ((host.flags.reshape(pair.fine_shape) & SOLID) == 0)] + backend.initialize(rr.initial_state(pair.coarse_shape, (0.0, 0.0), 1.0, root=True), + rr.initial_state(pair.fine_shape, pair.fine_origin, 0.5)) + probe_p = dict(p, center=[(80.5 - pair.fine_origin[0]) / 0.5 - 0.5, + (81.0 - pair.fine_origin[1]) / 0.5 - 0.5], radius=40.0) + else: + tmp = Path(stack.enter_context(TemporaryDirectory())) + (tmp / "lbm.json").write_text(json.dumps(cfgdict)) + (tmp / "body.json").write_text(json.dumps(bodydict)) + sim = stack.enter_context(Simulation(str(tmp / "lbm.json"), str(tmp / "body.json"), + device_id=0)) + sim.initialize() + probe_p = p + fluid = (sim.field.flag.reshape(p["ny"], p["nx"]) & FLUID) != 0 + begin = 0 + while begin < lp["steps"]: + rr.check_deadline(a.deadline) + n = min(lp["interval"], lp["steps"] - begin) + if refined: + impulse = rr.advance_refined_interval(backend, n, a.deadline, "window") + fields = [rr.macros(g) for g in backend.download()] + rho, ux, uy = fields[1] + fluid = masks[1] + vals = [tuple(v[m] for v in f) for m, f in zip(masks, fields)] + else: + impulse, done = np.zeros(2), 0 + while done < n: + batch = min(100, n - done) + sim.run(batch) + impulse += sim.read_force(0).astype(np.float64) * batch + done += batch + m = sim.get_macroscopic() + rho, ux, uy = [np.asarray(m[k], np.float64).reshape(p["ny"], p["nx"]) + for k in ("rho", "ux", "uy")] + vals = [(rho[fluid], ux[fluid], uy[fluid])] + begin += n + row = dict(step=begin, tstar=begin * rr.SPEED / p["diameter"], + **sample(impulse, n, p["diameter"], + h.pressure_probes(rho, fluid, probe_p)["dp"]), + rho_min=float(min(v[0].min() for v in vals)), + rho_max=float(max(v[0].max() for v in vals)), + max_mach=float(max(np.sqrt(3 * (v[1]**2 + v[2]**2)).max() for v in vals))) + rr.ensure_finite(row) + rows.append(row) + if not healthy(row): + break + summary = summarize(rows) if healthy(rows[-1]) else dict(status="UNHEALTHY", steady=False, + verdict=None) + return dict(case=case, collision=a.collision, lattice=lp, rows=rows, **summary) + + +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, + 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") + + +def main(): + parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + parser.add_argument("--case", choices=["both", "refined", "uniform"], default="both") + 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) + parser.add_argument("--plan", action="store_true", help="Write the frozen plan; no CUDA") + parser.add_argument("--max-wall-seconds", type=float) + a = parser.parse_args() + a.source = a.source.resolve() + a.output.mkdir(parents=True, exist_ok=False) + (a.output / "plan.json").write_text(json.dumps(plan(a.collision), indent=2, allow_nan=False)) + if a.plan: + return + if not a.max_wall_seconds or a.max_wall_seconds <= 0 or not os.environ.get("CUDA_VISIBLE_DEVICES"): + raise SystemExit("GPU run needs --max-wall-seconds > 0 and explicit CUDA_VISIBLE_DEVICES") + h, _ = rr.helpers(a.source) + rr.config = config # host_geometry builds its root configuration through this name + a.deadline = time.monotonic() + a.max_wall_seconds + cases = ["refined", "uniform"] if a.case == "both" else [a.case] + results = [run_case(a, h, case) for case in cases] + (a.output / "result.json").write_text(json.dumps(results, indent=2, allow_nan=False)) + + +if __name__ == "__main__": + main()