# CelerisLab/tests/validation/run_kan99b_rotating_cylinder.py """Kan99b MRT-only rotating-cylinder validation runner. This script follows ``docs/validation_specs/Kan99b_validation.md`` for the current round: - Primary matrix: K1-K5 with collision fixed to MRT. - Primary inlet: regularized (uniform profile). - Extra control: K2 with ``zou_he_local`` inlet for sensitivity only. """ from __future__ import annotations import argparse import csv import json import os import sys import tempfile from dataclasses import dataclass from typing import Any, Dict, List, Optional, Sequence, Tuple import numpy as np import pycuda.driver as cuda _REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..")) _DEFAULT_LBM = os.path.join(_REPO, "src", "CelerisLab", "configs", "config_lbm.json") U_INF = 0.03 D_LATTICE = 30.0 R_LATTICE = 15.0 _STORE_PRECISION = "FP32" _DDF_SHIFTING = False KAN99B_ANCHOR = { "St": 0.1655, "mean_cl": -2.4881, "mean_cd": 1.1040, "amp_cl": 0.3631, "amp_cd": 0.0993, } ANCHOR_BANDS = { "St": 0.03, "mean_cl": 0.04, "mean_cd": 0.05, "amp_cl": 0.08, "amp_cd": 0.10, } @dataclass(frozen=True) class DomainSpec: key: str nx: int ny: int center: Tuple[float, float] @dataclass(frozen=True) class KanCase: case_id: str re: float alpha: float steps: int burn: int @dataclass(frozen=True) class RunSpec: case_id: str variant: str domain: str re: float alpha: float inlet_scheme: str steps: int burn: int CASES: Tuple[KanCase, ...] = ( KanCase("K1", 100.0, 0.5, 200_000, 80_000), KanCase("K2", 100.0, 1.0, 200_000, 80_000), KanCase("K3", 60.0, 1.6, 240_000, 120_000), KanCase("K4", 100.0, 2.0, 240_000, 120_000), KanCase("K5", 160.0, 2.0, 240_000, 120_000), ) CASE_K0 = KanCase("K0", 100.0, 0.0, 180_000, 72_000) def _domain_specs() -> Dict[str, DomainSpec]: return { "S": DomainSpec("S", 1081, 481, (360.0, 240.0)), "M": DomainSpec("M", 1351, 601, (450.0, 300.0)), "L": DomainSpec("L", 1801, 721, (600.0, 360.0)), } def _load_json(path: str) -> dict: with open(path, "r", encoding="utf-8") as f: return json.load(f) def _write_json(path: str, payload: dict) -> None: with open(path, "w", encoding="utf-8") as f: json.dump(payload, f, indent=2) def _nu_from_re(re: float) -> float: return U_INF * D_LATTICE / float(re) def _omega_body(alpha: float) -> float: return 2.0 * float(alpha) * U_INF / D_LATTICE def _run_id(spec: RunSpec) -> str: a = f"{spec.alpha:.3f}".replace(".", "p") return ( f"{spec.case_id.lower()}_{spec.variant}_dom{spec.domain}_re{int(spec.re)}_a{a}_" f"{spec.inlet_scheme.lower()}_mrt" ) def _build_cfg( base_cfg: dict, *, nx: int, ny: int, re: float, inlet_scheme: str, ) -> dict: cfg = json.loads(json.dumps(base_cfg)) cfg["grid"]["nx"] = int(nx) cfg["grid"]["ny"] = int(ny) cfg["grid"]["nz"] = 1 cfg["physics"]["velocity"] = float(U_INF) cfg["physics"]["viscosity"] = float(_nu_from_re(re)) cfg["physics"]["rho"] = 1.0 cfg["method"]["collision"] = "MRT" cfg["method"]["streaming"] = "double_buffer" cfg["method"]["store_precision"] = _STORE_PRECISION cfg["method"]["ddf_shifting"] = _DDF_SHIFTING cfg["method"]["les"]["enabled"] = False cfg["method"]["inlet"]["profile"] = "uniform" cfg["method"]["inlet"]["scheme"] = str(inlet_scheme) cfg["method"]["outlet"]["mode"] = "neq_extrap" cfg["method"]["y_wall_bc"] = "free_slip" return cfg def _body_doc(center: Tuple[float, float], alpha: float) -> dict: return { "objects": [ { "type": "cylinder", "center": [float(center[0]), float(center[1])], "radius": float(R_LATTICE), "omega": float(_omega_body(alpha)), } ] } def _rfft_spectrum(x: np.ndarray, sample_dt: float) -> Tuple[np.ndarray, np.ndarray]: arr = np.asarray(x, dtype=np.float64) if arr.size < 64: return np.zeros(0, dtype=np.float64), np.zeros(0, dtype=np.float64) arr = arr - np.mean(arr) spec = np.abs(np.fft.rfft(arr * np.hanning(arr.size))) ** 2 freqs = np.fft.rfftfreq(arr.size, d=float(sample_dt)) return freqs.astype(np.float64), spec.astype(np.float64) def _peak_freq_parabolic(freqs: np.ndarray, spec: np.ndarray, idx: int) -> float: i = int(np.clip(idx, 0, spec.size - 1)) if i <= 0 or i + 1 >= spec.size: return float(freqs[i]) y0 = np.log(spec[i - 1] + 1e-30) y1 = np.log(spec[i] + 1e-30) y2 = np.log(spec[i + 1] + 1e-30) den = y0 - 2.0 * y1 + y2 if abs(den) < 1e-20: return float(freqs[i]) delta = float(np.clip(0.5 * (y0 - y2) / den, -1.0, 1.0)) return float(freqs[i]) + delta * float(freqs[i + 1] - freqs[i]) def _st_from_lift(lift: np.ndarray, sample_dt: float) -> float: freqs, spec = _rfft_spectrum(lift, sample_dt=sample_dt) if freqs.size <= 1: return float("nan") idx = int(np.argmax(spec[1:])) + 1 f_peak = _peak_freq_parabolic(freqs, spec, idx) return float(f_peak * D_LATTICE / U_INF) def _cycle_half_p2p(y: np.ndarray) -> float: arr = np.asarray(y, dtype=np.float64) if arr.size < 8: return float("nan") centered = arr - np.mean(arr) crossing = np.where((centered[:-1] <= 0.0) & (centered[1:] > 0.0))[0] if crossing.size >= 2: amps: List[float] = [] for i in range(crossing.size - 1): seg = arr[crossing[i] + 1 : crossing[i + 1] + 1] if seg.size >= 3: amps.append(0.5 * (float(np.max(seg)) - float(np.min(seg)))) if amps: return float(np.mean(amps)) return 0.5 * (float(np.max(arr)) - float(np.min(arr))) def _vorticity_z(ux: np.ndarray, uy: np.ndarray) -> np.ndarray: ux = np.asarray(ux, dtype=np.float64) uy = np.asarray(uy, dtype=np.float64) return np.gradient(uy, axis=1) - np.gradient(ux, axis=0) def _save_vorticity_png(path: str, ux: np.ndarray, uy: np.ndarray, title: str) -> None: try: import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt except ImportError: return omega = _vorticity_z(ux, uy) abs_o = np.abs(omega[np.isfinite(omega)]) vmax = float(np.percentile(abs_o, 99.5)) if abs_o.size else 1.0 if vmax <= 0.0: vmax = 1.0 ny, nx = omega.shape fig, ax = plt.subplots(figsize=(min(18.0, max(8.0, nx / 100.0)), min(12.0, max(3.0, ny / 40.0)))) im = ax.imshow( omega, origin="lower", aspect="equal", cmap="RdBu_r", vmin=-vmax, vmax=vmax, extent=(0, nx - 1, 0, ny - 1), ) ax.set_xlabel("x (lattice)") ax.set_ylabel("y (lattice)") ax.set_title(title) fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label="omega_z") fig.tight_layout() fig.savefig(path, dpi=150, bbox_inches="tight") plt.close(fig) def _run_one( spec: RunSpec, *, domain: DomainSpec, base_cfg: dict, out_dir: str, record_every: int, field_every: int, save_vorticity: bool, ) -> Dict[str, Any]: cfg = _build_cfg( base_cfg, nx=domain.nx, ny=domain.ny, re=spec.re, inlet_scheme=spec.inlet_scheme, ) body = _body_doc(domain.center, alpha=spec.alpha) tmpd = tempfile.mkdtemp(prefix="celeris_kan99b_mrt_") lbm_tmp = os.path.join(tmpd, "config_lbm.json") body_tmp = os.path.join(tmpd, "config_body.json") _write_json(lbm_tmp, cfg) _write_json(body_tmp, body) from CelerisLab import Simulation # noqa: WPS433 sim = Simulation(lbm_config_path=lbm_tmp, body_config_path=body_tmp) if sim.bodies.count < 1: sim.close() raise RuntimeError("Expected one cylinder in body config.") sim.bodies.get(0).state.omega = np.float32(_omega_body(spec.alpha)) sim.initialize() stream = cuda.Stream() rec = max(1, int(record_every)) total = int(spec.burn) + int(spec.steps) if total < 1: sim.close() raise ValueError("burn + steps must be >= 1") step_hist: List[int] = [] fx_hist: List[float] = [] fy_hist: List[float] = [] field_snapshots: List[str] = [] run_id = _run_id(spec) snap_dir = os.path.join(out_dir, "fields", run_id) if field_every > 0: os.makedirs(snap_dir, exist_ok=True) for step in range(1, total + 1): sim.bodies.zero_force_segment_async(stream) sim.stepper.step( 1, action_gpu=sim.bodies.action_gpu, obs_gpu=sim.bodies.obs_gpu, stream=stream, ) if step % rec == 0 or step == total: stream.synchronize() sim.bodies.download_obs_full_async(stream) stream.synchronize() force = sim.bodies.read_force(0) fx = float(force[0]) fy = float(force[1]) if not np.isfinite(fx) or not np.isfinite(fy): sim.close() raise RuntimeError(f"NaN/Inf force at step {step}") step_hist.append(step) fx_hist.append(fx) fy_hist.append(fy) if field_every > 0 and (step % int(field_every) == 0 or step == total): stream.synchronize() macro = sim.get_macroscopic() snap_path = os.path.join(snap_dir, f"macro_step{step:08d}.npz") np.savez_compressed( snap_path, rho=np.asarray(macro["rho"], dtype=np.float32), ux=np.asarray(macro["ux"], dtype=np.float32), uy=np.asarray(macro["uy"], dtype=np.float32), ) field_snapshots.append(snap_path) stream.synchronize() macro_last = sim.get_macroscopic() ux_last = np.asarray(macro_last["ux"], dtype=np.float64).reshape(domain.ny, domain.nx) uy_last = np.asarray(macro_last["uy"], dtype=np.float64).reshape(domain.ny, domain.nx) rho_last = np.asarray(macro_last["rho"], dtype=np.float64).reshape(domain.ny, domain.nx) sim.close() step_arr = np.asarray(step_hist, dtype=np.int64) fx_arr = np.asarray(fx_hist, dtype=np.float64) fy_arr = np.asarray(fy_hist, dtype=np.float64) burn_mask = step_arr >= int(spec.burn) if not np.any(burn_mask): burn_mask = np.ones_like(step_arr, dtype=bool) cl = 2.0 * fy_arr / (U_INF**2 * D_LATTICE) cd = 2.0 * fx_arr / (U_INF**2 * D_LATTICE) cl_tail = cl[burn_mask] cd_tail = cd[burn_mask] st = _st_from_lift(cl_tail, sample_dt=float(rec)) amp_cl = _cycle_half_p2p(cl_tail) amp_cd = _cycle_half_p2p(cd_tail) csv_dir = os.path.join(out_dir, "force_csv") os.makedirs(csv_dir, exist_ok=True) csv_path = os.path.join(csv_dir, f"{run_id}.csv") with open(csv_path, "w", newline="", encoding="utf-8") as f: w = csv.writer(f) w.writerow(["step", "fx", "fy", "cd", "cl"]) for i, s in enumerate(step_arr.tolist()): w.writerow([s, fx_arr[i], fy_arr[i], cd[i], cl[i]]) if save_vorticity: vdir = os.path.join(out_dir, "vorticity") os.makedirs(vdir, exist_ok=True) _save_vorticity_png( os.path.join(vdir, f"{run_id}.png"), ux_last, uy_last, title=( f"Kan99b {spec.case_id} {spec.variant} MRT dom={spec.domain} " f"Re={spec.re:.0f} alpha={spec.alpha:.3f} inlet={spec.inlet_scheme}" ), ) return { "run_id": run_id, "case_id": spec.case_id, "variant": spec.variant, "collision": "MRT", "inlet_scheme": spec.inlet_scheme, "inlet_profile": "uniform", "domain": spec.domain, "re": float(spec.re), "alpha": float(spec.alpha), "omega_body": float(_omega_body(spec.alpha)), "nu": float(_nu_from_re(spec.re)), "steps": int(spec.steps), "burn_in": int(spec.burn), "total_steps": int(total), "record_every": int(rec), "n_samples": int(step_arr.size), "St": float(st), "st": float(st), "mean_cl": float(np.mean(cl_tail)), "mean_cd": float(np.mean(cd_tail)), "amp_cl": float(amp_cl), "amp_cd": float(amp_cd), "rho_min_final": float(np.min(rho_last)), "rho_max_final": float(np.max(rho_last)), "grid": {"nx": int(domain.nx), "ny": int(domain.ny), "diameter": int(D_LATTICE)}, "beta_real": None, "Re_real": None, "re_real": None, "force_csv": csv_path, "field_snapshots": field_snapshots, } def _rel_err(measured: float, ref: float) -> Optional[float]: if not np.isfinite(measured) or ref == 0.0: return None return abs(float(measured) - float(ref)) / abs(float(ref)) def _k2_anchor_gate(rows: Sequence[Dict[str, Any]]) -> List[Dict[str, Any]]: """Evaluate K2 rows against Kan99b anchor tolerances.""" out: List[Dict[str, Any]] = [] for row in rows: if row.get("case_id") != "K2" or "error" in row: continue rel = { "St": _rel_err(row["St"], KAN99B_ANCHOR["St"]), "mean_cl": _rel_err(row["mean_cl"], KAN99B_ANCHOR["mean_cl"]), "mean_cd": _rel_err(row["mean_cd"], KAN99B_ANCHOR["mean_cd"]), "amp_cl": _rel_err(row["amp_cl"], KAN99B_ANCHOR["amp_cl"]), "amp_cd": _rel_err(row["amp_cd"], KAN99B_ANCHOR["amp_cd"]), } pass_bands = { key: (rel[key] is not None and rel[key] <= ANCHOR_BANDS[key]) for key in ANCHOR_BANDS } out.append( { "run_id": row["run_id"], "variant": row["variant"], "inlet_scheme": row["inlet_scheme"], "rel_err": rel, "pass_bands": pass_bands, "pass_all": bool(all(pass_bands.values())), } ) return out def _build_runs( cases: Sequence[KanCase], *, domain: str, include_k2_control: bool, steps_override: int, burn_override: int, ) -> List[RunSpec]: runs: List[RunSpec] = [] for case in cases: steps = int(steps_override) if steps_override > 0 else int(case.steps) burn = int(burn_override) if burn_override > 0 else int(case.burn) runs.append( RunSpec( case_id=case.case_id, variant="baseline", domain=domain, re=case.re, alpha=case.alpha, inlet_scheme="regularized", steps=steps, burn=burn, ) ) if include_k2_control and case.case_id == "K2": runs.append( RunSpec( case_id=case.case_id, variant="k2_inlet_control", domain=domain, re=case.re, alpha=case.alpha, inlet_scheme="zou_he_local", steps=steps, burn=burn, ) ) return runs def main() -> int: ap = argparse.ArgumentParser(description="Kan99b MRT-only primary matrix runner") ap.add_argument("--case", default="all", help='Case id K1-K5/K0 or "all"') ap.add_argument("--include-k0", action="store_true", help="Include optional K0 baseline.") ap.add_argument("--no-k2-control", action="store_true", help="Disable K2 zou_he_local control run.") ap.add_argument("--domain", default="M", choices=("S", "M", "L")) ap.add_argument("--steps", type=int, default=0, help="Override run steps for all selected runs.") ap.add_argument("--burn", type=int, default=0, help="Override burn steps for all selected runs.") ap.add_argument("--record-every", type=int, default=100) ap.add_argument("--field-every", type=int, default=0, help="Dump macro field .npz every N steps (0 disables).") ap.add_argument("--out-dir", type=str, default=os.path.join(_REPO, "tests", "output", "kan99b_validation")) ap.add_argument("--smoke", action="store_true", help="Very short run for wiring checks.") ap.add_argument("--save-vorticity", action="store_true", help="Save final vorticity PNG per run.") ap.add_argument("--store-precision", type=str, default="FP32", choices=("FP32", "FP16S"), help="DDF store precision (FP32, FP16S).") ap.add_argument("--ddf-shifting", action="store_true", help="Enable DDF shifting mode (f - w storage).") ap.add_argument("--json-out", type=str, default="", help="Optional explicit summary JSON output path.") args = ap.parse_args() # Global for _build_cfg() access global _STORE_PRECISION, _DDF_SHIFTING _STORE_PRECISION = str(args.store_precision).upper() _DDF_SHIFTING = bool(args.ddf_shifting) if not os.path.isfile(_DEFAULT_LBM): print(f"Missing base config: {_DEFAULT_LBM}", file=sys.stderr) return 2 base_cfg = _load_json(_DEFAULT_LBM) domains = _domain_specs() sel = str(args.case).upper() allowed = {case.case_id for case in CASES} | {"K0", "ALL"} if sel not in allowed: print("--case must be K0,K1,K2,K3,K4,K5,all", file=sys.stderr) return 2 selected: List[KanCase] = [] if sel == "ALL": selected.extend(CASES) if args.include_k0: selected.insert(0, CASE_K0) elif sel == "K0": selected.append(CASE_K0) else: selected.extend(case for case in CASES if case.case_id == sel) if not selected: print("No runs selected.", file=sys.stderr) return 2 steps_override = 2000 if args.smoke else max(0, int(args.steps)) burn_override = 800 if args.smoke else max(0, int(args.burn)) runs = _build_runs( selected, domain=args.domain, include_k2_control=not bool(args.no_k2_control), steps_override=steps_override, burn_override=burn_override, ) out_dir = os.path.abspath(args.out_dir) os.makedirs(out_dir, exist_ok=True) rows: List[Dict[str, Any]] = [] for spec in runs: print( f"--- {spec.case_id} {spec.variant} MRT dom={spec.domain} Re={spec.re:.0f} " f"alpha={spec.alpha:.3f} inlet={spec.inlet_scheme} burn={spec.burn} steps={spec.steps} ---", flush=True, ) try: row = _run_one( spec, domain=domains[spec.domain], base_cfg=base_cfg, out_dir=out_dir, record_every=max(1, int(args.record_every)), field_every=max(0, int(args.field_every)), save_vorticity=bool(args.save_vorticity), ) if spec.case_id == "K2": rel = _rel_err(row["St"], KAN99B_ANCHOR["St"]) row["St_error_pct"] = 100.0 * rel if rel is not None else None else: row["St_error_pct"] = None rows.append(row) print( f" St={row['St']:.5f} mean_CL={row['mean_cl']:.4f} mean_CD={row['mean_cd']:.4f} " f"C'L={row['amp_cl']:.4f} C'D={row['amp_cd']:.4f}", flush=True, ) except Exception as exc: # noqa: BLE001 rows.append( { "run_id": _run_id(spec), "case_id": spec.case_id, "variant": spec.variant, "collision": "MRT", "inlet_scheme": spec.inlet_scheme, "inlet_profile": "uniform", "domain": spec.domain, "re": float(spec.re), "alpha": float(spec.alpha), "steps": int(spec.steps), "burn_in": int(spec.burn), "error": str(exc), } ) print(f"FAILED: {exc}", flush=True) k2_gate = _k2_anchor_gate(rows) print("\n=== Kan99b K2 gate summary ===", flush=True) print(json.dumps({"k2_runs": k2_gate}, indent=2), flush=True) summary_csv = os.path.join(out_dir, "summary_runs.csv") csv_keys = [ "run_id", "case_id", "variant", "collision", "inlet_scheme", "inlet_profile", "domain", "re", "alpha", "omega_body", "nu", "burn_in", "steps", "total_steps", "record_every", "n_samples", "St", "st", "St_error_pct", "mean_cl", "mean_cd", "amp_cl", "amp_cd", "rho_min_final", "rho_max_final", "force_csv", "error", ] with open(summary_csv, "w", newline="", encoding="utf-8") as f: writer = csv.DictWriter(f, fieldnames=csv_keys) writer.writeheader() for row in rows: writer.writerow({k: row.get(k, "") for k in csv_keys}) summary = { "contract": { "collision": "MRT", "primary_inlet_scheme": "regularized", "k2_control_inlet_scheme": "zou_he_local", "inlet_profile": "uniform", "y_wall_bc": "free_slip", "outlet_mode": "neq_extrap", "streaming": "double_buffer", "store_precision": "FP32", "les_enabled": False, }, "requested": { "case": args.case, "include_k0": bool(args.include_k0), "include_k2_control": not bool(args.no_k2_control), "domain": args.domain, "smoke": bool(args.smoke), "steps_override": int(steps_override), "burn_override": int(burn_override), "record_every": int(args.record_every), "field_every": int(args.field_every), "save_vorticity": bool(args.save_vorticity), }, "counts": { "requested_runs": len(runs), "completed_runs": sum(1 for r in rows if "error" not in r), "failed_runs": sum(1 for r in rows if "error" in r), }, "k2_gate": k2_gate, "rows": rows, } json_out = ( os.path.abspath(args.json_out) if args.json_out.strip() else os.path.join(out_dir, "summary_runs.json") ) json_out_dir = os.path.dirname(json_out) if json_out_dir: os.makedirs(json_out_dir, exist_ok=True) _write_json(json_out, summary) print(f"Wrote: {summary_csv}", flush=True) print(f"Wrote: {json_out}", flush=True) return 0 if __name__ == "__main__": raise SystemExit(main())