# CelerisLab/tests/validation/run_sah04_st_matrix.py """Sah04 MRT-only Strouhal validation on S1-S4 anchors. This runner implements the current validation contract in ``tests/Sah04_validation.md``: - Cases: S1-S4 only (hard periodic anchors). - Collision: MRT only. - Inlet: parabolic + channel_stabilized. - Walls: no-slip channel (from base config + confined geometry). - Grid policy: S3/S4 can use configurable refined diameter for diagnostics. Usage:: conda run -n pycuda_3_10 python tests/run_sah04_st_matrix.py conda run -n pycuda_3_10 python tests/run_sah04_st_matrix.py --case S3 --smoke conda run -n pycuda_3_10 python tests/run_sah04_st_matrix.py --gate-pct 10 --json-out tests/output/sah04_mrt/summary.json """ from __future__ import annotations import argparse 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 _PKG_ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..")) _DEFAULT_LBM = os.path.join(_PKG_ROOT, "src", "CelerisLab", "configs", "config_lbm.json") _BASE_D = 30.0 @dataclass(frozen=True) class Sah04Case: """One hard benchmark case from Sah04_validation.md.""" case_id: str beta_nominal: float re_nominal: float target_st: float h_fluid: int steps: int burn: int @dataclass(frozen=True) class CaseGeometry: """Resolved lattice geometry for one case.""" diameter: float h_fluid: int nx: int ny: int center_x: float center_y: float radius: float beta_real: float wall_gap_cells: float CASES: Tuple[Sah04Case, ...] = ( Sah04Case("S1", 0.3, 100.0, 0.2115, 100, 120_000, 45_000), Sah04Case("S2", 0.5, 200.0, 0.3513, 60, 120_000, 45_000), Sah04Case("S3", 0.8, 160.0, 0.5537, 38, 220_000, 99_000), Sah04Case("S4", 0.9, 200.0, 0.5314, 33, 220_000, 99_000), ) 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 rfft_power_spectrum(samples: np.ndarray, *, sample_dt: float) -> Tuple[np.ndarray, np.ndarray]: """Mean-subtracted Hanning-windowed signal to rFFT power spectrum.""" x = np.asarray(samples, dtype=np.float64) x = x - np.mean(x) n = x.size if n < 64: return np.zeros(0, dtype=np.float64), np.zeros(0, dtype=np.float64) win = np.hanning(n) spec = np.abs(np.fft.rfft(x * win)) ** 2 freqs = np.fft.rfftfreq(n, d=float(sample_dt)) return freqs.astype(np.float64), spec.astype(np.float64) def _parabolic_peak_freq(freqs: np.ndarray, spec: np.ndarray, idx: int) -> float: """Sub-bin frequency estimate with local log-parabolic interpolation.""" 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 _strouhal_from_lift( lift: np.ndarray, *, diameter: float, u_max: float, sample_dt: float, f_hz_min: float, f_hz_max: float, ) -> Tuple[float, float]: """Return guided Strouhal and guided dominant frequency.""" freqs, spec = rfft_power_spectrum(lift, sample_dt=sample_dt) if freqs.size == 0: return float("nan"), float("nan") band = (freqs >= float(f_hz_min)) & (freqs <= float(f_hz_max)) if not np.any(band): return float("nan"), float("nan") f0 = 0.5 * (float(f_hz_min) + float(f_hz_max)) sigma = max(1e-12, 0.18 * f0) weight = np.exp(-((freqs - f0) / sigma) ** 2) idx = int(np.argmax(spec * band.astype(np.float64) * weight)) f_peak = _parabolic_peak_freq(freqs, spec, idx) return float(f_peak * diameter / u_max), float(f_peak) def shedding_freq_band_hz( target_st: float, u_max: float, diameter: float, *, half_width: float = 0.42, ) -> Tuple[float, float]: """Frequency band around target shedding frequency for robust FFT pick.""" f0 = float(target_st) * float(u_max) / float(diameter) return max(1e-8, f0 * (1.0 - half_width)), f0 * (1.0 + half_width) def vorticity_z_from_velocity(ux: np.ndarray, uy: np.ndarray) -> np.ndarray: """Return z-vorticity for 2D velocity fields.""" 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_final_vorticity_png(path: str, ux: np.ndarray, uy: np.ndarray, *, title: str) -> None: """Save final-step vorticity image; requires matplotlib.""" try: import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt except ImportError as exc: raise RuntimeError("save_final_vorticity_png requires matplotlib.") from exc omega = vorticity_z_from_velocity(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 _build_case_geometry( case: Sah04Case, *, refine_high_beta: bool, high_beta_diameter: float, diameter_override: Optional[float], h_fluid_override: Optional[int], ) -> CaseGeometry: """Map case to geometry, with optional high-beta refinement.""" diameter = _BASE_D h_fluid = int(case.h_fluid) if diameter_override is not None or h_fluid_override is not None: if diameter_override is None or h_fluid_override is None: raise ValueError("Both --diameter-override and --h-fluid-override must be set together.") diameter = float(diameter_override) h_fluid = int(h_fluid_override) elif refine_high_beta and case.case_id in ("S3", "S4"): # Diagnostic default keeps high-blockage runs around D~80 for faster sweeps. diameter = float(high_beta_diameter) h_fluid = int(round(float(diameter) / float(case.beta_nominal))) nx = int(80.0 * diameter + 2.0) ny = int(h_fluid + 2) center_x = 40.0 * diameter + 0.5 center_y = 0.5 * float(h_fluid) + 0.5 radius = 0.5 * diameter beta_real = float(diameter / float(h_fluid)) wall_gap_cells = 0.5 * float(h_fluid - diameter) return CaseGeometry( diameter=diameter, h_fluid=h_fluid, nx=nx, ny=ny, center_x=center_x, center_y=center_y, radius=radius, beta_real=beta_real, wall_gap_cells=wall_gap_cells, ) def _relative_error(measured: float, target: float) -> Optional[float]: if not np.isfinite(measured) or target <= 0.0: return None return abs(float(measured) - float(target)) / float(target) def _realized_umax(ux: np.ndarray, *, probe_x: int) -> float: """Estimate developed centerline maximum from one downstream vertical profile.""" # Skip top/bottom walls and sample at a downstream x station. prof = np.asarray(ux[1:-1, int(probe_x)], dtype=np.float64) if prof.size == 0: return float("nan") return float(np.max(prof)) def run_one_simulation( case: Sah04Case, geometry: CaseGeometry, *, collision: str, outlet_mode: str, inlet_profile: str, inlet_scheme: str, u_max_nominal: float, steps: int, burn: int, record_every: int, f_hz_min: float, f_hz_max: float, dump_npz_path: Optional[str] = None, final_vorticity_png_path: Optional[str] = None, crash_dump_dir: Optional[str] = None, ) -> Dict[str, Any]: """Build config, run simulation, and return measured metrics.""" nu = float(u_max_nominal * geometry.diameter / case.re_nominal) u0_mean = float(u_max_nominal / 1.5) cfg = _load_json(_DEFAULT_LBM) cfg["grid"]["nx"] = int(geometry.nx) cfg["grid"]["ny"] = int(geometry.ny) cfg["grid"]["nz"] = 1 cfg["physics"]["viscosity"] = float(nu) cfg["physics"]["velocity"] = float(u0_mean) cfg["physics"]["rho"] = 1.0 cfg["method"]["collision"] = str(collision).upper() cfg["method"]["streaming"] = "double_buffer" cfg["method"]["les"]["enabled"] = False cfg["method"]["inlet"]["profile"] = str(inlet_profile) cfg["method"]["inlet"]["scheme"] = str(inlet_scheme) cfg["method"]["outlet"]["mode"] = str(outlet_mode) body_doc = { "objects": [ { "type": "cylinder", "center": [float(geometry.center_x), float(geometry.center_y)], "radius": float(geometry.radius), } ] } tmpd = tempfile.mkdtemp(prefix="celeris_sah04_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_doc) from CelerisLab import Simulation # noqa: WPS433 sim = Simulation(lbm_config_path=lbm_tmp, body_config_path=body_tmp) sim.initialize() stream = cuda.Stream() rec_every = max(1, int(record_every)) lift_hist: List[float] = [] fx_hist: List[float] = [] step_hist: List[int] = [] n_curved = int(sim.field.n_curved) fallback_links = int(sim.bodies.fallback_link_count()) low_q_links = int(sim.bodies.low_q_link_count()) for step in range(1, int(steps) + 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_every == 0 or step == int(steps): stream.synchronize() sim.bodies.download_obs_full_async(stream) stream.synchronize() fvec = sim.bodies.read_force(0, normalize=False) lift = float(fvec[1]) drag = float(fvec[0]) if not np.isfinite(lift) or not np.isfinite(drag): crash_npz_path: Optional[str] = None crash_png_path: Optional[str] = None if crash_dump_dir: os.makedirs(crash_dump_dir, exist_ok=True) stream.synchronize() macro_bad = sim.get_macroscopic() ux_bad = np.asarray(macro_bad["ux"], dtype=np.float64).reshape(geometry.ny, geometry.nx) uy_bad = np.asarray(macro_bad["uy"], dtype=np.float64).reshape(geometry.ny, geometry.nx) rho_bad = np.asarray(macro_bad["rho"], dtype=np.float64).reshape(geometry.ny, geometry.nx) prefix = f"{case.case_id.lower()}_{str(collision).lower()}_crash_step{step}" crash_npz_path = os.path.join(crash_dump_dir, f"{prefix}.npz") np.savez_compressed( crash_npz_path, rho=rho_bad.astype(np.float32), ux=ux_bad.astype(np.float32), uy=uy_bad.astype(np.float32), step=np.array([int(step)], dtype=np.int64), case_id=np.array([case.case_id]), collision=np.array([str(collision)]), inlet_scheme=np.array([str(inlet_scheme)]), outlet_mode=np.array([str(outlet_mode)]), ) crash_png_path = os.path.join(crash_dump_dir, f"{prefix}.png") try: save_final_vorticity_png( crash_png_path, ux_bad, uy_bad, title=( f"Sah04 {case.case_id} {collision} crash@{step} " f"D={int(geometry.diameter)} H={geometry.h_fluid}" ), ) except Exception: # noqa: BLE001 crash_png_path = None sim.close() msg = f"NaN/Inf force at step {step}" if crash_npz_path: msg += f"; crash_npz={crash_npz_path}" if crash_png_path: msg += f"; crash_png={crash_png_path}" raise RuntimeError(msg) lift_hist.append(lift) fx_hist.append(drag) step_hist.append(step) stream.synchronize() macro_last = sim.get_macroscopic() ux_last = np.asarray(macro_last["ux"], dtype=np.float64).reshape(geometry.ny, geometry.nx) uy_last = np.asarray(macro_last["uy"], dtype=np.float64).reshape(geometry.ny, geometry.nx) rho_last = np.asarray(macro_last["rho"], dtype=np.float64).reshape(geometry.ny, geometry.nx) sim.close() if final_vorticity_png_path: out_dir = os.path.dirname(os.path.abspath(final_vorticity_png_path)) if out_dir: os.makedirs(out_dir, exist_ok=True) save_final_vorticity_png( final_vorticity_png_path, ux_last, uy_last, title=f"Sah04 {case.case_id} MRT Re={case.re_nominal:.0f} beta_nom={case.beta_nominal:.1f}", ) lift_arr = np.asarray(lift_hist, dtype=np.float64) fx_arr = np.asarray(fx_hist, dtype=np.float64) step_arr = np.asarray(step_hist, dtype=np.int64) burn_idx = min(int(burn) // rec_every, max(0, lift_arr.size - 16)) lift_tail = lift_arr[burn_idx:] st, f_peak = _strouhal_from_lift( lift_tail, diameter=float(geometry.diameter), u_max=float(u_max_nominal), sample_dt=float(rec_every), f_hz_min=float(f_hz_min), f_hz_max=float(f_hz_max), ) mean_cd = ( float(np.mean(fx_arr[burn_idx:]) * 2.0 / (u_max_nominal**2 * geometry.diameter)) if fx_arr.size else float("nan") ) beta_real = float(geometry.beta_real) probe_x = min(geometry.nx - 2, max(2, geometry.nx - 10)) u_max_real = _realized_umax(ux_last, probe_x=probe_x) re_real = float(u_max_real * geometry.diameter / nu) if np.isfinite(u_max_real) and nu > 0.0 else float("nan") if dump_npz_path: freqs, power = rfft_power_spectrum(lift_tail, sample_dt=float(rec_every)) out_dir = os.path.dirname(os.path.abspath(dump_npz_path)) if out_dir: os.makedirs(out_dir, exist_ok=True) np.savez_compressed( dump_npz_path, lift_samples=lift_arr, drag_samples=fx_arr, sample_lbm_step=step_arr, burn_index_samples=int(burn_idx), record_every_lbm_steps=int(rec_every), freqs_hz_post_burn=freqs, power_post_burn=power, rho_final=rho_last.astype(np.float32), ux_final=ux_last.astype(np.float32), uy_final=uy_last.astype(np.float32), st=np.array([st], dtype=np.float64), f_peak=np.array([f_peak], dtype=np.float64), re_real=np.array([re_real], dtype=np.float64), beta_real=np.array([beta_real], dtype=np.float64), ) return { "collision": str(collision).upper(), "inlet_profile": inlet_profile, "inlet_scheme": inlet_scheme, "St": float(st), "f_peak_per_step": float(f_peak), "mean_Cd": float(mean_cd), "Re_real": float(re_real), "U_max_real": float(u_max_real), "beta_real": float(beta_real), "n_curved": n_curved, "fallback_links": fallback_links, "low_q_links": low_q_links, "rho_min_final": float(np.min(rho_last)), "rho_max_final": float(np.max(rho_last)), "n_lift_samples": int(lift_arr.size), } def evaluate_rows(rows: Sequence[Dict[str, Any]], *, gate_pct: float) -> Dict[str, Any]: """Aggregate pass/fail summary for S1-S4 hard anchors.""" valid_rows = [r for r in rows if "error" not in r] st_errs = [r.get("St_error_pct") for r in valid_rows if r.get("St_error_pct") is not None] pass_count = sum(1 for v in st_errs if float(v) <= float(gate_pct)) return { "gate_pct": float(gate_pct), "cases_total": len(rows), "cases_completed": len(valid_rows), "cases_failed": sum(1 for r in rows if "error" in r), "cases_within_gate": int(pass_count), "pass_gate_all_completed": bool(len(valid_rows) > 0 and pass_count == len(valid_rows)), } def main() -> int: ap = argparse.ArgumentParser(description="Sah04 MRT-only S1-S4 validation runner") ap.add_argument("--case", default="all", help='S1-S4 or "all"') ap.add_argument("--collision", default="MRT", choices=("SRT", "TRT", "MRT")) ap.add_argument("--outlet", default="neq_extrap", choices=("neq_extrap", "zero_gradient", "blended")) ap.add_argument("--inlet-profile", default="parabolic", choices=("parabolic",)) ap.add_argument( "--inlet-scheme", default="channel_stabilized", choices=("channel_stabilized", "regularized", "zou_he_local", "equilibrium"), ) ap.add_argument("--record-every", type=int, default=5) ap.add_argument("--smoke", action="store_true", help="Short run for wiring checks.") ap.add_argument("--steps", type=int, default=None, help="Override case steps (ignored with --smoke).") ap.add_argument("--burn", type=int, default=None, help="Override case burn (ignored with --smoke).") ap.add_argument("--gate-pct", type=float, default=5.0, help="Pass gate for St relative error percent.") ap.add_argument("--json-out", type=str, default=None, help="Write summary JSON.") ap.add_argument("--dump-npz-dir", type=str, default=None, help="Optional directory for case NPZ dumps.") ap.add_argument("--final-vorticity-dir", type=str, default=None, help="Optional directory for final vorticity PNG.") ap.add_argument( "--crash-dump-dir", type=str, default=None, help="Optional directory to dump full flowfield NPZ/PNG immediately before crash.", ) ap.add_argument( "--no-refine-high-beta", action="store_true", help="Disable default refined geometry for S3/S4 (debug only).", ) ap.add_argument( "--high-beta-diameter", type=float, default=80.0, help="Refined diameter for S3/S4 when high-beta refinement is enabled.", ) ap.add_argument("--diameter-override", type=float, default=None, help="Override cylinder diameter for selected case.") ap.add_argument("--h-fluid-override", type=int, default=None, help="Override fluid height H for selected case.") args = ap.parse_args() if not os.path.isfile(_DEFAULT_LBM): print(f"Missing base config: {_DEFAULT_LBM}", file=sys.stderr) return 2 selected_case = str(args.case).upper() if selected_case != "ALL" and selected_case not in {c.case_id for c in CASES}: print("--case must be one of S1,S2,S3,S4,all", file=sys.stderr) return 2 cases_to_run = [c for c in CASES if selected_case == "ALL" or c.case_id == selected_case] if args.dump_npz_dir: os.makedirs(args.dump_npz_dir, exist_ok=True) if args.final_vorticity_dir: os.makedirs(args.final_vorticity_dir, exist_ok=True) rows: List[Dict[str, Any]] = [] for case in cases_to_run: geometry = _build_case_geometry( case, refine_high_beta=not bool(args.no_refine_high_beta), high_beta_diameter=float(args.high_beta_diameter), diameter_override=args.diameter_override, h_fluid_override=args.h_fluid_override, ) steps = 5000 if args.smoke else (int(args.steps) if args.steps is not None else case.steps) burn = 1500 if args.smoke else (int(args.burn) if args.burn is not None else case.burn) f_lo, f_hi = shedding_freq_band_hz(case.target_st, 0.1, geometry.diameter) npz_path = os.path.join(args.dump_npz_dir, f"{case.case_id.lower()}_mrt.npz") if args.dump_npz_dir else None vort_path = ( os.path.join(args.final_vorticity_dir, f"{case.case_id.lower()}_mrt_laststep.png") if args.final_vorticity_dir else None ) print( f"--- {case.case_id} {args.collision} beta_nom={case.beta_nominal:.1f} Re_nom={case.re_nominal:.0f} " f"D={int(geometry.diameter)} H={geometry.h_fluid} gap~{geometry.wall_gap_cells:.2f} " f"steps={steps} burn={burn} inlet={args.inlet_scheme}/{args.inlet_profile} ---", flush=True, ) try: out = run_one_simulation( case, geometry, collision=args.collision, outlet_mode=args.outlet, inlet_profile=args.inlet_profile, inlet_scheme=args.inlet_scheme, u_max_nominal=0.1, steps=steps, burn=burn, record_every=int(args.record_every), f_hz_min=f_lo, f_hz_max=f_hi, dump_npz_path=npz_path, final_vorticity_png_path=vort_path, crash_dump_dir=args.crash_dump_dir, ) except Exception as exc: # noqa: BLE001 rows.append( { "case_id": case.case_id, "collision": str(args.collision).upper(), "inlet_scheme": args.inlet_scheme, "inlet_profile": args.inlet_profile, "outlet": args.outlet, "grid": {"nx": geometry.nx, "ny": geometry.ny, "diameter": int(geometry.diameter), "h_fluid": geometry.h_fluid}, "error": str(exc), } ) print(f"FAILED: {exc}", flush=True) continue rel_err = _relative_error(out["St"], case.target_st) st_err_pct = (100.0 * rel_err) if rel_err is not None else None row = { "case_id": case.case_id, "collision": out["collision"], "inlet_scheme": out["inlet_scheme"], "inlet_profile": out["inlet_profile"], "grid": {"nx": geometry.nx, "ny": geometry.ny, "diameter": int(geometry.diameter), "h_fluid": geometry.h_fluid}, "steps": int(steps), "burn_in": int(burn), "Re_nominal": float(case.re_nominal), "Re_real": out["Re_real"], "beta_nominal": float(case.beta_nominal), "beta_real": out["beta_real"], "wall_gap_cells": geometry.wall_gap_cells, "target_St": float(case.target_st), "St": out["St"], "St_error_pct": float(st_err_pct) if st_err_pct is not None and np.isfinite(st_err_pct) else None, "gate_pct": float(args.gate_pct), "gate_pass": bool(st_err_pct is not None and st_err_pct <= float(args.gate_pct)), "mean_Cd": out["mean_Cd"], "U_max_real": out["U_max_real"], "rho_min_final": out["rho_min_final"], "rho_max_final": out["rho_max_final"], "n_curved": out["n_curved"], "fallback_links": out["fallback_links"], "low_q_links": out["low_q_links"], "n_lift_samples": out["n_lift_samples"], } rows.append(row) st_err_txt = f"{row['St_error_pct']:.2f}%" if row["St_error_pct"] is not None else "n/a" re_real_txt = f"{row['Re_real']:.2f}" if np.isfinite(row["Re_real"]) else "nan" print( f" St={row['St']:.5f} target={row['target_St']:.5f} err={st_err_txt} " f"Re_real={re_real_txt} beta_real={row['beta_real']:.4f} " f"[{'PASS' if row['gate_pass'] else 'CHECK'}]", flush=True, ) evaluation = evaluate_rows(rows, gate_pct=float(args.gate_pct)) print("\n=== Sah04 MRT S1-S4 summary ===", flush=True) print(json.dumps(evaluation, indent=2), flush=True) if args.json_out: json_out_path = os.path.abspath(args.json_out) json_out_dir = os.path.dirname(json_out_path) if json_out_dir: os.makedirs(json_out_dir, exist_ok=True) _write_json( json_out_path, { "requested": { "case": args.case, "outlet": args.outlet, "inlet_profile": args.inlet_profile, "inlet_scheme": args.inlet_scheme, "record_every": int(args.record_every), "smoke": bool(args.smoke), "steps_override": args.steps, "burn_override": args.burn, "gate_pct": float(args.gate_pct), }, "rows": rows, "evaluation": evaluation, }, ) print(f"Wrote: {json_out_path}", flush=True) return 0 if __name__ == "__main__": raise SystemExit(main())