#!/usr/bin/env python3 """Run a real finite-volume GPU primitive over exported OpenFOAM arrays.""" from __future__ import annotations import argparse import os import shutil import sys import time from dataclasses import dataclass 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 verify_gpu_step_timing import json_ready, nvidia_device_identity, prepare_case, time_openfoam_step, write_report ROOT = Path(__file__).resolve().parents[1] DEFAULT_WORK = ROOT / "tmp/gpu_algorithm_check" DEFAULT_REPORT = DEFAULT_WORK / "report.json" ALGORITHM_NAME = "cell_flux_imbalance" @dataclass(frozen=True) class FluxInputs: owner: np.ndarray neighbour: np.ndarray phi: np.ndarray n_cells: int n_internal_faces: int source_summary: dict[str, Any] class VerificationFailure(Exception): """Verifier failure with JSON-reportable details.""" def __init__(self, message: str, *, details: Mapping[str, Any] | None = None): super().__init__(message) self.details = dict(details or {}) @qd.kernel def zero_cell_flux_imbalance(n_cells: int, out: qd.types.NDArray[qd.f64, 1]) -> None: for cell in range(n_cells): out[cell] = 0.0 @qd.kernel def cell_flux_imbalance( n_internal_faces: int, owner: qd.types.NDArray[qd.i32, 1], neighbour: qd.types.NDArray[qd.i32, 1], phi_internal: qd.types.NDArray[qd.f64, 1], out: qd.types.NDArray[qd.f64, 1], ) -> None: for face in range(n_internal_faces): flux = phi_internal[face] qd.atomic_add(out[owner[face]], flux) qd.atomic_add(out[neighbour[face]], -flux) def _require_mapping(value: Any, path: str) -> Mapping[str, Any]: if not isinstance(value, Mapping): raise ValueError(f"{path}: expected mapping, got {type(value).__name__}") return value def _contiguous_1d_array(value: Any, path: str) -> np.ndarray: array = np.asarray(value) if array.ndim != 1: raise ValueError(f"{path}: expected 1-D array, got shape {array.shape}") return np.ascontiguousarray(array) def load_openfoam_flux_inputs(case: Path) -> FluxInputs: import foam_stepper as foam state = foam.Case(case).make_stepper().export_state(required_fields=("phi",)) mesh = _require_mapping(state.get("mesh"), "mesh") sizes = _require_mapping(mesh.get("sizes"), "mesh.sizes") connectivity = _require_mapping(mesh.get("connectivity"), "mesh.connectivity") fields = _require_mapping(state.get("fields"), "fields") phi_field = _require_mapping(fields.get("phi"), "fields.phi") try: n_cells = int(sizes["n_cells"]) n_internal_faces = int(sizes["n_internal_faces"]) except KeyError as exc: raise ValueError(f"mesh.sizes.{exc.args[0]}: missing required size") from exc owner_raw = _contiguous_1d_array(connectivity.get("owner"), "mesh.connectivity.owner") neighbour_raw = _contiguous_1d_array(connectivity.get("neighbour"), "mesh.connectivity.neighbour") phi_raw = _contiguous_1d_array(phi_field.get("internal"), "fields.phi.internal") expected_shape = (n_internal_faces,) for path, array in ( ("mesh.connectivity.owner", owner_raw), ("mesh.connectivity.neighbour", neighbour_raw), ("fields.phi.internal", phi_raw), ): if array.shape != expected_shape: raise ValueError(f"{path}: expected shape {expected_shape}, got {array.shape}") if not np.issubdtype(owner_raw.dtype, np.integer): raise ValueError(f"mesh.connectivity.owner: expected integer dtype, got {owner_raw.dtype}") if not np.issubdtype(neighbour_raw.dtype, np.integer): raise ValueError(f"mesh.connectivity.neighbour: expected integer dtype, got {neighbour_raw.dtype}") if not np.issubdtype(phi_raw.dtype, np.floating): raise ValueError(f"fields.phi.internal: expected floating dtype, got {phi_raw.dtype}") owner_min = int(owner_raw.min(initial=0)) if owner_raw.size else 0 owner_max = int(owner_raw.max(initial=0)) if owner_raw.size else -1 neighbour_min = int(neighbour_raw.min(initial=0)) if neighbour_raw.size else 0 neighbour_max = int(neighbour_raw.max(initial=0)) if neighbour_raw.size else -1 if owner_min < 0 or neighbour_min < 0 or owner_max >= n_cells or neighbour_max >= n_cells: raise ValueError( "mesh.connectivity owner/neighbour indices out of cell range: " f"owner=[{owner_min}, {owner_max}], neighbour=[{neighbour_min}, {neighbour_max}], n_cells={n_cells}" ) if owner_max > np.iinfo(np.int32).max or neighbour_max > np.iinfo(np.int32).max: raise ValueError("mesh.connectivity owner/neighbour exceed int32 GPU index range") owner = np.ascontiguousarray(owner_raw.astype(np.int32, copy=False)) neighbour = np.ascontiguousarray(neighbour_raw.astype(np.int32, copy=False)) phi = np.ascontiguousarray(phi_raw.astype(np.float64, copy=False)) return FluxInputs( owner=owner, neighbour=neighbour, phi=phi, n_cells=n_cells, n_internal_faces=n_internal_faces, source_summary={ "case": str(case), "mesh": { "n_cells": n_cells, "n_internal_faces": n_internal_faces, "owner_shape": list(owner_raw.shape), "owner_dtype": str(owner_raw.dtype), "neighbour_shape": list(neighbour_raw.shape), "neighbour_dtype": str(neighbour_raw.dtype), }, "fields": { "phi": { "entity_kind": phi_field.get("entity_kind"), "entity_count": int(phi_field.get("entity_count", -1)), "internal_shape": list(phi_raw.shape), "internal_dtype": str(phi_raw.dtype), } }, "gpu_input_dtypes": { "owner": str(owner.dtype), "neighbour": str(neighbour.dtype), "phi_internal": str(phi.dtype), }, }, ) def run_gpu_flux_imbalance(inputs: FluxInputs, *, 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) owner_gpu = qd.ndarray(qd.i32, shape=inputs.owner.shape) neighbour_gpu = qd.ndarray(qd.i32, shape=inputs.neighbour.shape) phi_gpu = qd.ndarray(qd.f64, shape=inputs.phi.shape) out_gpu = qd.ndarray(qd.f64, shape=(inputs.n_cells,)) owner_gpu.from_numpy(inputs.owner) neighbour_gpu.from_numpy(inputs.neighbour) phi_gpu.from_numpy(inputs.phi) zero_cell_flux_imbalance(inputs.n_cells, out_gpu) cell_flux_imbalance(inputs.n_internal_faces, owner_gpu, neighbour_gpu, phi_gpu, out_gpu) qd.sync() qd.profiler.clear_kernel_profiler_info() start = time.perf_counter() for _ in range(repeats): zero_cell_flux_imbalance(inputs.n_cells, out_gpu) cell_flux_imbalance(inputs.n_internal_faces, owner_gpu, neighbour_gpu, phi_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(ALGORITHM_NAME in name for name in generated_kernel_names): raise AssertionError( f"Quadrants CUDA profiler recorded no {ALGORITHM_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 = np.asarray(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": [ALGORITHM_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 actual, evidence def compute_cpu_flux_imbalance_reference(inputs: FluxInputs) -> np.ndarray: reference = np.zeros(inputs.n_cells, dtype=np.float64) np.add.at(reference, inputs.owner, inputs.phi) np.add.at(reference, inputs.neighbour, -inputs.phi) return reference def compare_flux_imbalance(actual: np.ndarray, expected: np.ndarray, *, rtol: float, atol: float) -> dict[str, Any]: actual_array = np.asarray(actual) expected_array = np.asarray(expected) if actual_array.shape != expected_array.shape: report = { "allclose": False, "reason": "shape_mismatch", "actual_shape": list(actual_array.shape), "expected_shape": list(expected_array.shape), "actual_dtype": str(actual_array.dtype), "expected_dtype": str(expected_array.dtype), "rtol": rtol, "atol": atol, } raise VerificationFailure("cell_flux_imbalance shape mismatch", details={"comparison": report}) if not np.issubdtype(actual_array.dtype, np.floating): report = { "allclose": False, "reason": "actual_dtype_not_floating", "shape": list(actual_array.shape), "actual_dtype": str(actual_array.dtype), "expected_dtype": str(expected_array.dtype), "rtol": rtol, "atol": atol, } raise VerificationFailure("GPU output dtype is not floating", details={"comparison": report}) if not np.issubdtype(expected_array.dtype, np.floating): report = { "allclose": False, "reason": "expected_dtype_not_floating", "shape": list(actual_array.shape), "actual_dtype": str(actual_array.dtype), "expected_dtype": str(expected_array.dtype), "rtol": rtol, "atol": atol, } raise VerificationFailure("CPU reference dtype is not floating", details={"comparison": report}) actual64 = actual_array.astype(np.float64, copy=False) expected64 = expected_array.astype(np.float64, copy=False) abs_diff = np.abs(actual64 - expected64) max_abs = float(np.max(abs_diff)) if abs_diff.size else 0.0 denominator = np.maximum(np.abs(expected64), atol) rel_diff = np.divide(abs_diff, denominator, out=np.zeros_like(abs_diff), where=denominator > 0.0) max_rel = float(np.max(rel_diff)) if rel_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(actual64, expected64, rtol=rtol, atol=atol)) report = { "allclose": allclose, "shape": list(actual_array.shape), "dtype": str(actual_array.dtype), "actual_shape": list(actual_array.shape), "expected_shape": list(expected_array.shape), "actual_dtype": str(actual_array.dtype), "expected_dtype": str(expected_array.dtype), "rtol": rtol, "atol": atol, "max_abs_error": max_abs, "max_rel_error": max_rel, "largest_difference_index": largest_index, } if not allclose: raise VerificationFailure("cell_flux_imbalance numerical mismatch", details={"comparison": report}) return report 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" case = args.work / "gpu_algorithm_case" prepare_case(openfoam_case) prepare_case(case) openfoam_step = time_openfoam_step(openfoam_case) inputs = load_openfoam_flux_inputs(case) cpu_reference_start = time.perf_counter() expected = compute_cpu_flux_imbalance_reference(inputs) cpu_reference_wall_ms = (time.perf_counter() - cpu_reference_start) * 1000.0 actual, gpu_evidence = run_gpu_flux_imbalance(inputs, repeats=args.repeats) comparison = compare_flux_imbalance(actual, expected, rtol=args.rtol, atol=args.atol) conservation_sum = float(np.sum(actual, dtype=np.float64)) if actual.size else 0.0 gpu_wall_per_step = gpu_evidence["wall_ms_per_step"] cpu_reference_speedup = cpu_reference_wall_ms / gpu_wall_per_step if gpu_wall_per_step > 0.0 else float("inf") openfoam_step_speedup = openfoam_step["wall_ms"] / gpu_wall_per_step if gpu_wall_per_step > 0.0 else float("inf") return { "status": "passed", "backend_requested": "gpu", "backend_selected": "gpu", "device_kind": gpu_evidence["device_kind"], "device_name": gpu_evidence["device_name"], "device_uuid": gpu_evidence["device_uuid"], "used_cpu_fallback": False, "algorithm": { "name": ALGORITHM_NAME, "kind": "internal_face_finite_volume_primitive", "formula": "cell_flux_imbalance[cell] = sum(owner phi_internal) - sum(neighbour phi_internal)", "inputs": [ "mesh.connectivity.owner", "mesh.connectivity.neighbour", "mesh.sizes.n_internal_faces", "fields.phi.internal", ], "cpu_reference_complete": True, }, "openfoam_step": openfoam_step, "cpu_reference": { "operation": ALGORITHM_NAME, "implementation": "numpy.add.at owner(+phi) and neighbour(-phi)", "output": "cell_flux_imbalance", "output_shape": list(expected.shape), "output_dtype": str(expected.dtype), "wall_ms": cpu_reference_wall_ms, }, "openfoam_inputs": inputs.source_summary, "gpu_algorithm": { "operation": ALGORITHM_NAME, "output": "cell_flux_imbalance", "output_shape": list(actual.shape), "output_dtype": str(actual.dtype), "repeats": args.repeats, "evidence": gpu_evidence, }, "comparison": comparison, "timing": { "openfoam_step_wall_ms": openfoam_step["wall_ms"], "openfoam_operation": openfoam_step["operation"], "cpu_numpy_reference_wall_ms": cpu_reference_wall_ms, "gpu_algorithm_wall_ms_total": gpu_evidence["wall_ms_total"], "gpu_algorithm_wall_ms_per_step": gpu_evidence["wall_ms_per_step"], "gpu_algorithm_device_ms_total": gpu_evidence["device_time_ms_total"], "gpu_algorithm_device_ms_per_step": gpu_evidence["device_time_ms_per_step"], "gpu_repeats": args.repeats, "speedup_vs_openfoam_step_wall": openfoam_step_speedup, "speedup_vs_cpu_numpy_reference_wall": cpu_reference_speedup, }, "conservation_check": { "description": "owner and neighbour scatter signs should make the global internal-face imbalance sum cancel", "sum": conservation_sum, "abs_sum": abs(conservation_sum), }, "work": args.work, } 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=1e-10) parser.add_argument("--atol", type=float, default=1e-12) 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", "backend_requested": "gpu", "backend_selected": None, "used_cpu_fallback": False, "algorithm": {"name": ALGORITHM_NAME}, "work": args.work, "failure": {"type": type(exc).__name__, "message": str(exc)}, } details = getattr(exc, "details", None) if details: failure["failure"]["details"] = details write_report(failure, report_path) print(f"gpu algorithm verification failed: {exc}", file=sys.stderr) print(f"report={report_path}", file=sys.stderr) return 1 write_report(report, report_path) evidence = report["gpu_algorithm"]["evidence"] print("gpu algorithm verification passed") print(f"report={report_path}") print(f"algorithm={report['algorithm']['name']}") print(f"backend_selected={report['backend_selected']}") print(f"profile_record_count={evidence['profile_record_count']}") print(f"gpu_wall_ms_per_step={evidence['wall_ms_per_step']:.6f}") print(f"gpu_device_ms_per_step={evidence['device_time_ms_per_step']:.6f}") print(f"openfoam_step_wall_ms={report['timing']['openfoam_step_wall_ms']:.3f}") print(f"cpu_numpy_reference_wall_ms={report['timing']['cpu_numpy_reference_wall_ms']:.6f}") print(f"output_shape={report['gpu_algorithm']['output_shape']}") return 0 if __name__ == "__main__": raise SystemExit(main())