#!/usr/bin/env python3 """Compare OpenFOAM step timing with a proven GPU kernel over OpenFOAM-derived data. This is deliberately not a RANS-solver speedup claim. The GPU step here is the smallest honest executable GPU operation we can verify today: a Quadrants CUDA kernel consuming the OpenFOAM velocity field and producing per-cell squared speed. The report records OpenFOAM one-iteration timing next to the GPU kernel timing and fails unless the GPU kernel actually ran on the CUDA backend. """ from __future__ import annotations import argparse import json import os import shutil import subprocess import sys import time from pathlib import Path from typing import Any, Mapping def _drop_ambient_pythonpath() -> None: pythonpath = os.environ.pop("PYTHONPATH", "") if not pythonpath: return for entry in pythonpath.split(os.pathsep): if entry and entry in sys.path: sys.path.remove(entry) _drop_ambient_pythonpath() import numpy as np import quadrants as qd from quadrants.profiler.kernel_profiler import get_default_kernel_profiler from openfoam_env import apply_openfoam_env, openfoam_env ROOT = Path(__file__).resolve().parents[1] TUTORIAL = ROOT / "OpenFOAM-14/tutorials/incompressibleFluid/venturiTube" DEFAULT_WORK = ROOT / "tmp/gpu_step_timing_check" DEFAULT_REPORT = DEFAULT_WORK / "report.json" KERNEL_NAME = "squared_speed" HARNESS_SCOPE = ( "squared_speed is only a CUDA proof/timing harness over OpenFOAM-derived " "velocity data; it is not the target finite-volume or RANS GPU algorithm" ) @qd.kernel def squared_speed(n: int, u: qd.types.NDArray[qd.f32, 2], out: qd.types.NDArray[qd.f32, 1]) -> None: for i in range(n): ux = u[i, 0] uy = u[i, 1] uz = u[i, 2] out[i] = ux * ux + uy * uy + uz * uz def json_ready(value: Any) -> Any: if isinstance(value, Mapping): return {str(key): json_ready(item) for key, item in value.items()} if isinstance(value, (list, tuple)): return [json_ready(item) for item in value] if isinstance(value, np.generic): return value.item() if isinstance(value, np.ndarray): return value.tolist() if isinstance(value, Path): return str(value) return value def patch_control_dict(case: Path) -> None: path = case / "system/controlDict" text = path.read_text() replacements = { "startFrom startTime;": "startFrom startTime;", "endTime 1000;": "endTime 1;", "writeInterval 50;": "writeInterval 1;", } for old, new in replacements.items(): text = text.replace(old, new) path.write_text(text) def run_openfoam_command(cmd: list[str]) -> None: subprocess.run(cmd, cwd=ROOT, check=True, stdout=subprocess.DEVNULL, stderr=subprocess.STDOUT, env=openfoam_env()) def prepare_case(dst: Path) -> None: if dst.exists(): shutil.rmtree(dst) shutil.copytree(TUTORIAL, dst, ignore=shutil.ignore_patterns("processor*", "postProcessing", "*.log")) for orig in (dst / "0").glob("*.orig"): shutil.copyfile(orig, orig.with_suffix("")) patch_control_dict(dst) run_openfoam_command(["blockMesh", "-case", str(dst)]) run_openfoam_command(["createZones", "-case", str(dst)]) def import_foam() -> Any: apply_openfoam_env() import foam_stepper as foam return foam def time_openfoam_step(case: Path) -> dict[str, Any]: foam = import_foam() stepper = foam.Case(case).make_stepper() start = time.perf_counter() result = stepper.run_one_pimple_iteration() wall_ms = (time.perf_counter() - start) * 1000.0 fields = result.outputs["fields"] return { "operation": "foam_stepper.run_one_pimple_iteration", "backend": "OpenFOAM C++ CPU stepper", "wall_ms": wall_ms, "output_shapes": { "U": list(np.asarray(fields["U"].internal).shape), "p": list(np.asarray(fields["p"].internal).shape), "phi": list(np.asarray(fields["phi"].internal).shape), }, } def load_openfoam_velocity(case: Path) -> np.ndarray: foam = import_foam() fields = foam.Case(case).make_stepper().fields() velocity = np.asarray(fields.U.internal, dtype=np.float32) if velocity.ndim != 2 or velocity.shape[1] != 3: raise AssertionError(f"expected vector U field with shape (cells, 3), got {velocity.shape}") return np.ascontiguousarray(velocity) def nvidia_device_identity() -> dict[str, str | None]: try: completed = subprocess.run( ["nvidia-smi", "--query-gpu=name,uuid", "--format=csv,noheader,nounits"], check=True, stdout=subprocess.PIPE, stderr=subprocess.DEVNULL, text=True, ) except (FileNotFoundError, subprocess.CalledProcessError): return {"device_name": None, "device_uuid": None} first = completed.stdout.strip().splitlines()[0] if completed.stdout.strip() else "" if not first: return {"device_name": None, "device_uuid": None} parts = [part.strip() for part in first.split(",", 1)] return {"device_name": parts[0], "device_uuid": parts[1] if len(parts) > 1 else None} def compare_arrays(actual: np.ndarray, expected: np.ndarray, *, rtol: float, atol: float) -> dict[str, Any]: if actual.shape != expected.shape: raise AssertionError(f"GPU output shape mismatch: {actual.shape} != {expected.shape}") abs_diff = np.abs(actual - expected) max_abs = float(np.max(abs_diff)) if abs_diff.size else 0.0 max_rel = float(np.max(abs_diff / np.maximum(np.abs(expected), atol))) if abs_diff.size else 0.0 largest_index = None if abs_diff.size: largest_index = [int(index) for index in np.unravel_index(np.argmax(abs_diff), abs_diff.shape)] allclose = bool(np.allclose(actual, expected, rtol=rtol, atol=atol)) report = { "allclose": allclose, "shape": list(actual.shape), "dtype": str(actual.dtype), "rtol": rtol, "atol": atol, "max_abs_error": max_abs, "max_rel_error": max_rel, "largest_difference_index": largest_index, } if not allclose: raise AssertionError(f"GPU output mismatch: {json.dumps(report, sort_keys=True)}") return report def run_gpu_velocity_step(velocity: np.ndarray, *, repeats: int) -> tuple[np.ndarray, dict[str, Any]]: if repeats < 1: raise ValueError("repeats must be >= 1") qd.init(arch=qd.cuda, kernel_profiler=True) cells = int(velocity.shape[0]) u_gpu = qd.ndarray(qd.f32, shape=velocity.shape) out_gpu = qd.ndarray(qd.f32, shape=(cells,)) u_gpu.from_numpy(velocity) squared_speed(cells, u_gpu, out_gpu) qd.sync() qd.profiler.clear_kernel_profiler_info() start = time.perf_counter() for _ in range(repeats): squared_speed(cells, u_gpu, out_gpu) qd.sync() wall_ms = (time.perf_counter() - start) * 1000.0 profiler = get_default_kernel_profiler() profiler._update_records() records = list(profiler._traced_records) generated_kernel_names = sorted({str(record.name) for record in records}) if not any(KERNEL_NAME in name for name in generated_kernel_names): raise AssertionError(f"Quadrants CUDA profiler recorded no {KERNEL_NAME!r} kernel; recorded {generated_kernel_names}") device_time_ms_total = float(sum(record.kernel_time for record in records)) if device_time_ms_total <= 0.0: raise AssertionError(f"Quadrants CUDA profiler recorded nonpositive device time: {device_time_ms_total}") actual = out_gpu.to_numpy() identity = nvidia_device_identity() evidence = { "backend_requested": "gpu", "backend_selected": "gpu", "framework": "quadrants", "arch_requested": "cuda", "arch_selected": "cuda", "device_kind": "cuda", "device_name": identity["device_name"], "device_uuid": identity["device_uuid"], "kernel_names": [KERNEL_NAME], "generated_kernel_names": generated_kernel_names, "requested_kernel_calls": repeats, "profile_record_count": len(records), "used_cpu_fallback": False, "synchronized_before_timing": True, "synchronized_after_timing": True, "wall_ms_total": wall_ms, "wall_ms_per_step": wall_ms / repeats, "device_time_ms_total": device_time_ms_total, "device_time_ms_per_step": device_time_ms_total / repeats, "device_time_ms_min_record": float(min(record.kernel_time for record in records)), "device_time_ms_max_record": float(max(record.kernel_time for record in records)), } return np.asarray(actual), evidence def write_report(report: Mapping[str, Any], path: Path) -> None: path.parent.mkdir(parents=True, exist_ok=True) path.write_text(json.dumps(json_ready(report), indent=2, sort_keys=True, allow_nan=False) + "\n") def run_check(args: argparse.Namespace) -> dict[str, Any]: if args.work.exists(): shutil.rmtree(args.work) args.work.mkdir(parents=True) openfoam_case = args.work / "openfoam_step_case" gpu_case = args.work / "gpu_input_case" prepare_case(openfoam_case) prepare_case(gpu_case) openfoam_step = time_openfoam_step(openfoam_case) velocity = load_openfoam_velocity(gpu_case) expected = np.einsum("ij,ij->i", velocity, velocity).astype(np.float32, copy=False) actual, gpu_evidence = run_gpu_velocity_step(velocity, repeats=args.repeats) comparison = compare_arrays(actual, expected, rtol=args.rtol, atol=args.atol) speedup_vs_openfoam_wall = openfoam_step["wall_ms"] / gpu_evidence["wall_ms_per_step"] if gpu_evidence["wall_ms_per_step"] > 0 else float("inf") report = { "status": "passed", "scope": HARNESS_SCOPE, "algorithm": { "name": KERNEL_NAME, "role": "proof_timing_harness", "target_algorithm_complete": False, }, "work": args.work, "openfoam_step": openfoam_step, "gpu_step": { "operation": KERNEL_NAME, "input_field": "U", "input_shape": list(velocity.shape), "output": "squared_speed_per_cell", "output_shape": list(actual.shape), "output_dtype": str(actual.dtype), "repeats": args.repeats, "evidence": gpu_evidence, }, "comparison": comparison, "timing": { "openfoam_wall_ms": openfoam_step["wall_ms"], "gpu_wall_ms_total": gpu_evidence["wall_ms_total"], "gpu_wall_ms_per_step": gpu_evidence["wall_ms_per_step"], "gpu_device_ms_per_step": gpu_evidence["device_time_ms_per_step"], "speedup_vs_openfoam_wall": speedup_vs_openfoam_wall, }, } return report def parse_args(argv: list[str] | None = None) -> argparse.Namespace: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--work", type=Path, default=DEFAULT_WORK) parser.add_argument("--report", type=Path, default=None) parser.add_argument("--repeats", type=int, default=20) parser.add_argument("--rtol", type=float, default=5e-6) parser.add_argument("--atol", type=float, default=1e-6) return parser.parse_args(argv) def main(argv: list[str] | None = None) -> int: args = parse_args(argv) report_path = args.report if args.report is not None else args.work / "report.json" try: report = run_check(args) except Exception as exc: failure = { "status": "failed", "scope": HARNESS_SCOPE, "work": args.work, "failure": {"type": type(exc).__name__, "message": str(exc)}, } write_report(failure, report_path) print(f"gpu step timing verification failed: {exc}", file=sys.stderr) print(f"report={report_path}", file=sys.stderr) return 1 write_report(report, report_path) print("gpu step timing verification passed") print(f"report={report_path}") print(f"openfoam_wall_ms={report['timing']['openfoam_wall_ms']:.3f}") print(f"gpu_wall_ms_per_step={report['timing']['gpu_wall_ms_per_step']:.6f}") print(f"gpu_device_ms_per_step={report['timing']['gpu_device_ms_per_step']:.6f}") print(f"speedup_vs_openfoam_wall={report['timing']['speedup_vs_openfoam_wall']:.3f}") print(f"backend_selected={report['gpu_step']['evidence']['backend_selected']}") print(f"profile_record_count={report['gpu_step']['evidence']['profile_record_count']}") return 0 if __name__ == "__main__": raise SystemExit(main())