# 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 _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), } ] } from CelerisLab import Simulation with tempfile.TemporaryDirectory(prefix="celeris_sah04_mrt_") as tmpd: 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) with Simulation(lbm_config_path=lbm_tmp, body_config_path=body_tmp) as sim, sim.ctx.activate(): sim.initialize() stream = sim.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 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) 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) failed = evaluation["cases_failed"] > 0 if not args.smoke: failed |= not evaluation["pass_gate_all_completed"] return int(failed) if __name__ == "__main__": raise SystemExit(main())