diff --git a/.gitignore b/.gitignore index 84892d5..83b43c8 100644 --- a/.gitignore +++ b/.gitignore @@ -2,6 +2,7 @@ .loop/ **/__pycache__ python/src/foam_stepper.egg-info/ +tmp/ # foreign OpenFOAM-14/ diff --git a/python/src/foam_stepper/__init__.py b/python/src/foam_stepper/__init__.py index 2718b8d..de8ebd0 100644 --- a/python/src/foam_stepper/__init__.py +++ b/python/src/foam_stepper/__init__.py @@ -3,7 +3,7 @@ from __future__ import annotations from pathlib import Path -from typing import Any +from typing import Any, Iterable from ._runtime import configure_openfoam_environment @@ -23,6 +23,15 @@ from .types import ( TransformResult, _wrap_value, ) +from .state import ( + SolverStateValidationError, + describe_matrix_state, + describe_solver_state, + export_matrix_state, + export_solver_state, + validate_matrix_state, + validate_solver_state, +) OpenFoamError = _native.OpenFoamError @@ -89,6 +98,12 @@ class SimpleStepper: def state(self) -> dict[str, Any]: return dict(self._native.state()) + def export_state(self, *, required_fields: Iterable[str] = ()) -> dict[str, Any]: + return export_solver_state(self.mesh(), self.fields(), required_fields=required_fields) + + def state_summary(self, *, required_fields: Iterable[str] = ()) -> dict[str, Any]: + return describe_solver_state(self.export_state(required_fields=required_fields)) + def control_dict(self) -> str: return self._native.control_dict() @@ -201,5 +216,12 @@ __all__ = [ "SolveResult", "SourceLocation", "TransformResult", + "SolverStateValidationError", + "describe_matrix_state", + "describe_solver_state", + "export_matrix_state", + "export_solver_state", + "validate_matrix_state", + "validate_solver_state", "version", ] diff --git a/python/src/foam_stepper/gpu/__init__.py b/python/src/foam_stepper/gpu/__init__.py new file mode 100644 index 0000000..1b8189e --- /dev/null +++ b/python/src/foam_stepper/gpu/__init__.py @@ -0,0 +1,5 @@ +"""Reusable GPU backend modules for foam_stepper.""" + +from __future__ import annotations + +__all__ = ["backend", "constants", "kernels", "linear_solve"] diff --git a/python/src/foam_stepper/gpu/backend.py b/python/src/foam_stepper/gpu/backend.py new file mode 100644 index 0000000..7ff4683 --- /dev/null +++ b/python/src/foam_stepper/gpu/backend.py @@ -0,0 +1,1912 @@ +"""Reusable GPU backend for the AirfRANS RANS solver verifier. + +The verifier under ``scripts/`` is an acceptance harness. This module owns the +GPU backend contract, GPU-side state transfer, CUDA solver-stage graph, and +backend diagnostics used by ``--backend gpu``. +""" + +from __future__ import annotations + +import dataclasses +import math +import os +import subprocess +import time +from collections.abc import Iterable, Mapping +from pathlib import Path +from typing import Any + +import numpy as np +import quadrants as qd + +from .constants import ( + DEFAULT_LAMINAR_NU, + DEFAULT_MOMENTUM_PBICGSTAB_ITERATIONS, + DEFAULT_MOMENTUM_PBICGSTAB_RESIDUAL_TOLERANCE_SQUARED, + DEFAULT_PRESSURE_CG_ITERATIONS, + DEFAULT_MOMENTUM_RELAXATION_ALPHA, + DEFAULT_SIMPLE_CONSISTENT_RATU_FACTOR, + DEFAULT_OMEGA_WALL_BETA1, + FULL_GPU_RANS_GUARD, + GPU_BACKEND_PROVIDER, + GPU_INPUT_SCHEMA_VERSION, + GPU_KERNELS_MISSING, + GPU_NUMERICAL_MISMATCH, + GPU_PRIMITIVE_NAME, + GPU_PRIMITIVE_PROOF, + GPU_RUNTIME_UNAVAILABLE, + GPU_STAGE_CAPABILITY_PREFIX, + GPU_STAGE_KERNEL_ENTRYPOINTS, + GPU_STAGE_UNSUPPORTED, + REQUIRED_FIELDS, + STAGE_OBSERVABILITY_GROUPS, + TURBULENCE_FIELDS, +) +from .kernels import ( + gpu_rans_face_flux_copy, + gpu_rans_final_correction, + gpu_rans_momentum_assembly, + gpu_rans_momentum_diffusion_coefficients, + gpu_rans_momentum_wall_diffusion_coefficients, + gpu_rans_momentum_convection_coefficients, + gpu_rans_momentum_equation_relaxation, + gpu_rans_momentum_hbyA_face_accumulate, + gpu_rans_momentum_hbyA_finish, + gpu_rans_momentum_hbyA_fixed_value_boundary, + gpu_rans_momentum_hbyA_source, + gpu_rans_pressure_assembly, + gpu_rans_pressure_inputs, + gpu_rans_consistent_rAtU, + gpu_rans_pressure_laplacian_coefficients, + gpu_rans_pressure_flux_correction, + gpu_rans_pressure_source_from_flux, + gpu_rans_pressure_source_from_boundary_flux, + gpu_rans_pressure_mixed_boundary_laplacian, + gpu_rans_surface_flux_from_cells, + gpu_rans_turbulence_update, + gpu_rans_omega_wall_update, + gpu_rans_pressure_velocity_correction, +) +from .linear_solve import ( + gpu_ldu_pbicgstab_vector_asymmetric_faces, + gpu_ldu_pcg_scalar_symmetric_faces, +) + +BACKEND_CHOICES = ("auto", "cpu", "gpu") +BACKEND_FAILURE = "backend_execution_failure" +COMPARISON_FAILURE = "numerical_comparison_failure" + +FIELD_COMPARISON_ATTRIBUTION = { + "U": ("final_correction", "linear_solve_results", "momentum_assembly", "turbulence_updates"), + "p": ("linear_solve_results", "pressure_assembly", "final_correction"), + "phi": ("final_correction", "linear_solve_results", "pressure_assembly"), + "nut": ("turbulence_updates",), + "k": ("turbulence_updates",), + "omega": ("turbulence_updates",), +} + + +class BackendExecutionError(Exception): + """Backend selection or execution failure before numerical comparison.""" + + def __init__(self, step: str, message: str, *, details: Mapping[str, Any] | None = None) -> None: + super().__init__(message) + self.step = step + self.details = dict(details or {}) + + +@dataclasses.dataclass(frozen=True) +class GpuSourceLocation: + file: str + function: str + lines: tuple[int, int] | None = None + + +@dataclasses.dataclass(frozen=True) +class GpuTransformResult: + name: str + phase: str + source: GpuSourceLocation + inputs: Mapping[str, Any] + outputs: Mapping[str, Any] + changed_fields: list[str] + metadata: Mapping[str, Any] + + +@dataclasses.dataclass(frozen=True) +class GpuPatchField: + name: str + type: str + values: np.ndarray + fixes_value: bool + assignable: bool + coupled: bool + updated: bool + patch_internal: np.ndarray | None + value_internal_coeffs: np.ndarray | None + value_boundary_coeffs: np.ndarray | None + gradient_internal_coeffs: np.ndarray | None + gradient_boundary_coeffs: np.ndarray | None + + +@dataclasses.dataclass(frozen=True) +class GpuField: + name: str + kind: str + dimensions: str + entity_kind: str + entity_count: int + internal: np.ndarray + boundary: Mapping[str, GpuPatchField] + + +def json_ready(value: Any) -> Any: + """Convert report values into strict JSON-compatible data.""" + + if dataclasses.is_dataclass(value) and not isinstance(value, type): + return json_ready(dataclasses.asdict(value)) + if isinstance(value, Path): + return str(value) + if isinstance(value, np.ndarray): + return array_stats(value) + if isinstance(value, np.generic): + return json_ready(value.item()) + if isinstance(value, float): + return value if math.isfinite(value) else None + if isinstance(value, (str, int, bool)) or value is None: + return value + if isinstance(value, Mapping): + return {str(key): json_ready(item) for key, item in value.items()} + if isinstance(value, (list, tuple, set)): + return [json_ready(item) for item in value] + return repr(value) + + +def array_shape(array: Any | None) -> list[int] | None: + if array is None: + return None + return [int(dim) for dim in np.asarray(array).shape] + + +def finite_float(value: Any) -> float | None: + try: + number = float(value) + except (TypeError, ValueError): + return None + return number if math.isfinite(number) else None + + +def array_stats(array: Any) -> dict[str, Any]: + arr = np.asarray(array) + out: dict[str, Any] = { + "shape": array_shape(arr), + "dtype": str(arr.dtype), + "size": int(arr.size), + } + if arr.size == 0 or not np.issubdtype(arr.dtype, np.number): + return out + + finite = np.isfinite(arr) + out["finite_count"] = int(np.count_nonzero(finite)) + out["nonfinite_count"] = int(arr.size - out["finite_count"]) + if out["finite_count"]: + finite_values = arr[finite] + out.update( + { + "min": finite_float(np.min(finite_values)), + "max": finite_float(np.max(finite_values)), + "mean": finite_float(np.mean(finite_values)), + } + ) + return out + + +def value_at(array: np.ndarray, index: tuple[int, ...]) -> Any: + value = array[index] + if isinstance(value, np.generic): + return json_ready(value.item()) + if isinstance(value, np.ndarray): + return json_ready(value.tolist()) + return json_ready(value) + + +def entity_value_at(array: np.ndarray, entity_index: int | None) -> Any: + if entity_index is None or array.ndim == 0: + return None + value = array[entity_index] + if isinstance(value, np.generic): + return json_ready(value.item()) + if isinstance(value, np.ndarray): + return json_ready(value.tolist()) + return json_ready(value) + + +def update_hash_text(digest: "hashlib._Hash", value: str) -> None: + encoded = value.encode("utf-8") + digest.update(len(encoded).to_bytes(8, "little")) + digest.update(encoded) + + +def update_hash_array(digest: "hashlib._Hash", label: str, array: Any) -> None: + arr = np.ascontiguousarray(np.asarray(array)) + update_hash_text(digest, label) + update_hash_text(digest, str(arr.dtype)) + update_hash_text(digest, repr(tuple(int(dim) for dim in arr.shape))) + digest.update(arr.tobytes()) + + +def source_summary(source: Any) -> dict[str, Any]: + return { + "file": getattr(source, "file", ""), + "function": getattr(source, "function", ""), + "lines": json_ready(getattr(source, "lines", None)), + } + + +def boundary_field_summary(patch: Any) -> dict[str, Any]: + return { + "type": getattr(patch, "type", ""), + "values_shape": array_shape(getattr(patch, "values", None)), + "fixes_value": bool(getattr(patch, "fixes_value", False)), + "assignable": bool(getattr(patch, "assignable", False)), + "coupled": bool(getattr(patch, "coupled", False)), + "updated": bool(getattr(patch, "updated", False)), + "patch_internal_shape": array_shape(getattr(patch, "patch_internal", None)), + "value_internal_coeffs_shape": array_shape(getattr(patch, "value_internal_coeffs", None)), + "value_boundary_coeffs_shape": array_shape(getattr(patch, "value_boundary_coeffs", None)), + "gradient_internal_coeffs_shape": array_shape(getattr(patch, "gradient_internal_coeffs", None)), + "gradient_boundary_coeffs_shape": array_shape(getattr(patch, "gradient_boundary_coeffs", None)), + } + + +def field_summary(field: Any) -> dict[str, Any]: + boundary = getattr(field, "boundary", {}) + return { + "name": getattr(field, "name", ""), + "kind": getattr(field, "kind", ""), + "dimensions": getattr(field, "dimensions", ""), + "entity_kind": getattr(field, "entity_kind", ""), + "entity_count": int(getattr(field, "entity_count", 0)), + "internal": array_stats(getattr(field, "internal")), + "boundary": {name: boundary_field_summary(patch) for name, patch in boundary.items()}, + } + + +def solve_summary(performance: Any) -> dict[str, Any]: + return { + "solver_name": getattr(performance, "solver_name", ""), + "field_name": getattr(performance, "field_name", ""), + "initial_residual": json_ready(getattr(performance, "initial_residual", None)), + "final_residual": json_ready(getattr(performance, "final_residual", None)), + "n_iterations": json_ready(getattr(performance, "n_iterations", None)), + "converged": bool(getattr(performance, "converged", False)), + "singular": bool(getattr(performance, "singular", False)), + } + + +def matrix_summary(matrix: Any) -> dict[str, Any]: + derived: dict[str, Any] = {} + for name in ("A", "H", "H1", "flux", "face_flux_correction"): + value = getattr(matrix, name, None) + if callable(value): + value = value() + if value is not None: + derived[name] = field_summary(value) + + for name in ("residual", "D", "DD"): + value = getattr(matrix, name, None) + if callable(value): + value = value() + if value is not None: + derived[name] = array_stats(value) + + return { + "name": getattr(matrix, "name", ""), + "field_name": getattr(matrix, "field_name", ""), + "value_rank": getattr(matrix, "value_rank", ""), + "dimensions": getattr(matrix, "dimensions", ""), + "has_diag": bool(getattr(matrix, "has_diag", False)), + "has_upper": bool(getattr(matrix, "has_upper", False)), + "has_lower": bool(getattr(matrix, "has_lower", False)), + "diagonal": bool(getattr(matrix, "diagonal", False)), + "symmetric": bool(getattr(matrix, "symmetric", False)), + "asymmetric": bool(getattr(matrix, "asymmetric", False)), + "diag": array_stats(getattr(matrix, "diag")), + "upper": None if getattr(matrix, "upper", None) is None else array_stats(getattr(matrix, "upper")), + "lower": None if getattr(matrix, "lower", None) is None else array_stats(getattr(matrix, "lower")), + "source": array_stats(getattr(matrix, "source")), + "psi": field_summary(getattr(matrix, "psi")), + "internal_coeff_shapes": [array_shape(item) for item in getattr(matrix, "internal_coeffs", [])], + "boundary_coeff_shapes": [array_shape(item) for item in getattr(matrix, "boundary_coeffs", [])], + "derived": derived, + } + + +def summarize_value(value: Any) -> Any: + if hasattr(value, "diag") and hasattr(value, "field_name") and hasattr(value, "source"): + return matrix_summary(value) + if hasattr(value, "internal") and hasattr(value, "entity_kind") and hasattr(value, "boundary"): + return field_summary(value) + if hasattr(value, "solver_name") and hasattr(value, "initial_residual"): + return solve_summary(value) + if hasattr(value, "name") and hasattr(value, "phase") and hasattr(value, "outputs"): + return transform_summary(value) + if isinstance(value, Mapping): + return {str(key): summarize_value(item) for key, item in value.items()} + if isinstance(value, list): + return [summarize_value(item) for item in value] + if isinstance(value, tuple): + return [summarize_value(item) for item in value] + return json_ready(value) + + +def transform_summary(result: Any) -> dict[str, Any]: + return { + "name": getattr(result, "name", ""), + "phase": getattr(result, "phase", ""), + "changed_fields": json_ready(getattr(result, "changed_fields", [])), + "source": source_summary(getattr(result, "source", None)), + "metadata": json_ready(getattr(result, "metadata", {})), + "outputs": summarize_value(getattr(result, "outputs", {})), + } + +def graph_stage_summaries(result: Any) -> list[dict[str, Any]]: + return [transform_summary(entry) for entry in getattr(result, "outputs", {}).get("graph", [])] + + +def stage_observability_report(mode: str, stages: list[dict[str, Any]], *, evidence: Mapping[str, Any]) -> dict[str, Any]: + stage_by_name = {stage["name"]: stage for stage in stages} + failures: list[dict[str, Any]] = [] + groups = [] + for group in STAGE_OBSERVABILITY_GROUPS: + group_name = group["name"] + required_stages = tuple(group[f"{mode}_stages"]) + missing_stages = [name for name in required_stages if name not in stage_by_name] + output_requirements = group.get(f"{mode}_outputs", {}) + output_checks = {} + for stage_name, required_outputs in output_requirements.items(): + stage = stage_by_name.get(stage_name) + if stage is None: + continue + outputs = stage.get("outputs", {}) + output_keys = set(outputs) if isinstance(outputs, Mapping) else set() + missing_outputs = [name for name in required_outputs if name not in output_keys] + output_checks[stage_name] = { + "required": list(required_outputs), + "available": sorted(output_keys), + "missing": missing_outputs, + } + if missing_outputs: + failures.append( + { + "path": f"stage_observability.{mode}.{group_name}.{stage_name}.outputs", + "message": "missing inspectable stage output", + "expected": list(required_outputs), + "actual": sorted(output_keys), + } + ) + if missing_stages: + failures.append( + { + "path": f"stage_observability.{mode}.{group_name}.stages", + "message": "missing required solver stage", + "expected": list(required_stages), + "actual": list(stage_by_name), + } + ) + groups.append( + { + "name": group_name, + "required_stages": list(required_stages), + "observed": not missing_stages and not any(check["missing"] for check in output_checks.values()), + "missing_stages": missing_stages, + "output_checks": output_checks, + } + ) + if failures: + raise BackendExecutionError( + f"{mode}_stage_observability", + f"{mode} GPU solver-stage observability is incomplete", + details={"failure_kind": GPU_STAGE_UNSUPPORTED, "failures": failures}, + ) + return { + "mode": mode, + "graph": [stage["name"] for stage in stages], + "groups": groups, + "evidence": json_ready(evidence), + } + + +def field_dict_to_mapping(fields: Any) -> dict[str, Any]: + if isinstance(fields, Mapping): + return dict(fields) + if hasattr(fields, "fields") and isinstance(fields.fields, Mapping): + return dict(fields.fields) + return {name: getattr(fields, name) for name in REQUIRED_FIELDS if hasattr(fields, name)} + + +def validation_details(exc: BaseException) -> dict[str, Any]: + if hasattr(exc, "to_dict"): + return json_ready(exc.to_dict()) + return {"type": type(exc).__name__, "message": str(exc)} + + +def read_fields(stepper: Any, label: str) -> dict[str, Any]: + try: + return field_dict_to_mapping(stepper.fields()) + except Exception as exc: + raise BackendExecutionError( + f"read_{label}_fields", + f"failed to read {label} fields for GPU backend", + details={"case": getattr(stepper, "case_path", ""), "cause": validation_details(exc)}, + ) from exc + + +def export_solver_state_checked(foam: Any, stepper: Any, label: str, fields: Mapping[str, Any] | None = None) -> dict[str, Any]: + try: + field_source = fields if fields is not None else stepper.fields() + return foam.export_solver_state(stepper.mesh(), field_source, required_fields=REQUIRED_FIELDS) + except Exception as exc: + raise BackendExecutionError( + f"export_{label}_state", + f"failed to export explicit {label} solver state for GPU backend", + details=validation_details(exc), + ) from exc + + +def visible_turbulence_fields(fields: Mapping[str, Any]) -> list[str]: + return [name for name in TURBULENCE_FIELDS if name in fields] + +def gpu_device_present() -> bool: + return any(Path(path).exists() for path in ("/dev/nvidia0", "/dev/dri/renderD128")) + + +def nvidia_device_identity() -> dict[str, Any]: + try: + completed = subprocess.run( + ["nvidia-smi", "--query-gpu=name,uuid", "--format=csv,noheader"], + check=True, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, + text=True, + timeout=10, + ) + except (FileNotFoundError, subprocess.CalledProcessError, subprocess.TimeoutExpired): + 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 gpu_stage_contracts() -> dict[str, dict[str, Any]]: + contracts: dict[str, dict[str, Any]] = {} + for group in STAGE_OBSERVABILITY_GROUPS: + name = str(group["name"]) + split_outputs = { + str(stage): list(outputs) + for stage, outputs in group.get("split_outputs", {}).items() + } + field_coverage = sorted( + { + field + for outputs in split_outputs.values() + for field in outputs + if field in REQUIRED_FIELDS + } + ) + kernel_entrypoints = list(GPU_STAGE_KERNEL_ENTRYPOINTS.get(name, ())) + contracts[name] = { + "required": True, + "supported": True, + "status": "kernels_registered" if kernel_entrypoints else "kernels_missing", + "run_one_stages": list(group.get("run_one_stages", ())), + "split_stages": list(group.get("split_stages", ())), + "expected_outputs": split_outputs, + "field_coverage": field_coverage, + "kernel_entrypoints": kernel_entrypoints, + "failure_kind_if_missing": None if kernel_entrypoints else GPU_KERNELS_MISSING, + } + return contracts + + +def gpu_backend_failure_taxonomy() -> dict[str, Any]: + return { + GPU_RUNTIME_UNAVAILABLE: { + "category": BACKEND_FAILURE, + "step": "select_backend", + "meaning": "CUDA/Quadrants could not provide a non-host GPU runtime.", + }, + GPU_KERNELS_MISSING: { + "category": BACKEND_FAILURE, + "step": "execute_gpu_solver_contract", + "meaning": "The GPU backend was selected, but required solver-stage kernels are not registered.", + }, + GPU_STAGE_UNSUPPORTED: { + "category": BACKEND_FAILURE, + "step": "execute_gpu_solver_contract", + "meaning": "The GPU backend explicitly cannot execute one or more required RANS solver stages.", + }, + GPU_NUMERICAL_MISMATCH: { + "category": COMPARISON_FAILURE, + "step": "compare_fields", + "meaning": "The GPU backend executed, but one or more required fields failed OpenFOAM oracle parity.", + "field_attribution": FIELD_COMPARISON_ATTRIBUTION, + }, + } + +def solver_state_arrays(value: Any, path: str = "") -> Iterable[tuple[str, np.ndarray]]: + if isinstance(value, np.ndarray): + yield path, value + return + if isinstance(value, Mapping): + for key, item in value.items(): + item_path = f"{path}.{key}" if path else str(key) + yield from solver_state_arrays(item, item_path) + return + if isinstance(value, (list, tuple)): + for index, item in enumerate(value): + yield from solver_state_arrays(item, f"{path}[{index}]") + + +def gpu_transfer_array_dtype(qd: Any, source: np.ndarray, path: str) -> tuple[Any, str, np.ndarray]: + if np.issubdtype(source.dtype, np.integer): + if source.size: + min_index = int(source.min()) + max_index = int(source.max()) + else: + min_index = 0 + max_index = -1 + int32 = np.iinfo(np.int32) + if min_index < int32.min or max_index > int32.max: + raise BackendExecutionError( + "prepare_gpu_solver_inputs", + f"{path}: integer values exceed int32 GPU index range", + details={ + "path": path, + "source_dtype": str(source.dtype), + "min": min_index, + "max": max_index, + "gpu_dtype": "i32", + }, + ) + return qd.i32, "i32", source.astype(np.int32, copy=False) + + if np.issubdtype(source.dtype, np.floating): + finite = np.isfinite(source) + if not bool(np.all(finite)): + raise BackendExecutionError( + "prepare_gpu_solver_inputs", + f"{path}: non-finite values cannot be copied into GPU solver inputs", + details={ + "path": path, + "source_dtype": str(source.dtype), + "shape": array_shape(source), + "nonfinite_count": int(source.size - np.count_nonzero(finite)), + }, + ) + if source.dtype == np.dtype("float32"): + return qd.f32, "f32", source + return qd.f64, "f64", source.astype(np.float64, copy=False) + + raise BackendExecutionError( + "prepare_gpu_solver_inputs", + f"{path}: unsupported GPU input dtype {source.dtype}", + details={"path": path, "source_dtype": str(source.dtype), "shape": array_shape(source)}, + ) + + +def transfer_solver_state_to_gpu(state: Mapping[str, Any]) -> dict[str, Any]: + import quadrants as qd + + arrays: dict[str, dict[str, Any]] = {} + device_arrays: list[Any] = [] + transferred_count = 0 + source_bytes = 0 + gpu_bytes = 0 + + for path, source in solver_state_arrays(state): + source_array = np.asarray(source) + source_bytes += int(source_array.nbytes) + entry: dict[str, Any] = { + "path": path, + "source_shape": array_shape(source_array), + "source_dtype": str(source_array.dtype), + "source_contiguous": bool(source_array.flags.c_contiguous), + "source_nbytes": int(source_array.nbytes), + "stats": array_stats(source_array), + } + if source_array.size == 0: + entry.update( + { + "status": "empty_array", + "transferred": False, + "gpu_shape": array_shape(source_array), + "reason": "zero-sized OpenFOAM patch array has no device allocation", + } + ) + arrays[path] = entry + continue + + qd_dtype, gpu_dtype, gpu_source = gpu_transfer_array_dtype(qd, source_array, path) + gpu_source = np.ascontiguousarray(gpu_source) + gpu_array = qd.ndarray(qd_dtype, shape=gpu_source.shape) + gpu_array.from_numpy(gpu_source) + device_arrays.append(gpu_array) + transferred_count += 1 + gpu_bytes += int(gpu_source.nbytes) + entry.update( + { + "status": "transferred", + "transferred": True, + "gpu_shape": array_shape(gpu_source), + "gpu_dtype": gpu_dtype, + "gpu_nbytes": int(gpu_source.nbytes), + "contiguous_for_transfer": True, + } + ) + arrays[path] = entry + + qd.sync() + return { + "framework": "quadrants", + "device": "cuda", + "array_count": len(arrays), + "transferred_count": transferred_count, + "empty_array_count": len(arrays) - transferred_count, + "source_nbytes": source_bytes, + "gpu_nbytes": gpu_bytes, + "arrays": arrays, + } + + +def gpu_required_field_inputs(state: Mapping[str, Any], transfer: Mapping[str, Any]) -> dict[str, Any]: + arrays = transfer.get("arrays") if isinstance(transfer.get("arrays"), Mapping) else {} + fields = state.get("fields") if isinstance(state.get("fields"), Mapping) else {} + out: dict[str, Any] = {} + for name in REQUIRED_FIELDS: + field = fields.get(name) if isinstance(fields.get(name), Mapping) else {} + boundary = field.get("boundary") if isinstance(field.get("boundary"), Mapping) else {} + internal_path = f"fields.{name}.internal" + internal = arrays.get(internal_path) if isinstance(arrays.get(internal_path), Mapping) else {} + out[name] = { + "status": "gpu_represented" if internal.get("transferred") is True else "missing_gpu_representation", + "rank": field.get("rank"), + "entity_kind": field.get("entity_kind"), + "entity_count": field.get("entity_count"), + "internal_path": internal_path, + "internal": internal, + "boundary_patch_count": len(boundary), + "boundary": { + patch_name: { + "type": patch.get("type"), + "values": arrays.get(f"fields.{name}.boundary.{patch_name}.values"), + "patch_internal": arrays.get(f"fields.{name}.boundary.{patch_name}.patch_internal"), + "value_internal_coeffs": arrays.get(f"fields.{name}.boundary.{patch_name}.value_internal_coeffs"), + "value_boundary_coeffs": arrays.get(f"fields.{name}.boundary.{patch_name}.value_boundary_coeffs"), + "gradient_internal_coeffs": arrays.get(f"fields.{name}.boundary.{patch_name}.gradient_internal_coeffs"), + "gradient_boundary_coeffs": arrays.get(f"fields.{name}.boundary.{patch_name}.gradient_boundary_coeffs"), + } + for patch_name, patch in boundary.items() + }, + } + return out + + +def gpu_mesh_input_summary(state: Mapping[str, Any], transfer: Mapping[str, Any]) -> dict[str, Any]: + arrays = transfer.get("arrays") if isinstance(transfer.get("arrays"), Mapping) else {} + mesh = state.get("mesh") if isinstance(state.get("mesh"), Mapping) else {} + sizes = mesh.get("sizes") if isinstance(mesh.get("sizes"), Mapping) else {} + connectivity = mesh.get("connectivity") if isinstance(mesh.get("connectivity"), Mapping) else {} + owner = np.asarray(connectivity.get("owner", [])) + neighbour = np.asarray(connectivity.get("neighbour", [])) + patches = mesh.get("patches") if isinstance(mesh.get("patches"), list) else [] + owner_range = [int(owner.min()), int(owner.max())] if owner.size else [None, None] + neighbour_range = [int(neighbour.min()), int(neighbour.max())] if neighbour.size else [None, None] + return { + "sizes": json_ready(sizes), + "topology_checks": { + "owner_shape": array_shape(owner), + "neighbour_shape": array_shape(neighbour), + "owner_cell_range": owner_range, + "neighbour_cell_range": neighbour_range, + "owner_neighbour_gpu_dtype": "i32", + "owner_neighbour_in_cell_range": bool( + owner.size + and neighbour.size + and min(owner_range[0], neighbour_range[0]) >= 0 + and max(owner_range[1], neighbour_range[1]) < int(sizes.get("n_cells", 0)) + ), + }, + "connectivity": { + "faces_offsets": arrays.get("mesh.connectivity.faces.offsets"), + "faces_values": arrays.get("mesh.connectivity.faces.values"), + "cells_offsets": arrays.get("mesh.connectivity.cells.offsets"), + "cells_values": arrays.get("mesh.connectivity.cells.values"), + "owner": arrays.get("mesh.connectivity.owner"), + "neighbour": arrays.get("mesh.connectivity.neighbour"), + "ldu_lower_addr": arrays.get("mesh.connectivity.ldu.lower_addr"), + "ldu_upper_addr": arrays.get("mesh.connectivity.ldu.upper_addr"), + }, + "geometry": { + "points": arrays.get("mesh.geometry.points"), + "V": arrays.get("mesh.geometry.V"), + "C": arrays.get("mesh.geometry.C"), + "Cf": arrays.get("mesh.geometry.Cf"), + "Sf": arrays.get("mesh.geometry.Sf"), + "magSf": arrays.get("mesh.geometry.magSf"), + }, + "patches": [ + { + "name": patch.get("name"), + "type": patch.get("type"), + "index": patch.get("index"), + "start": patch.get("start"), + "size": patch.get("size"), + "coupled": patch.get("coupled"), + "constraint": patch.get("constraint"), + "face_cells": arrays.get(f"mesh.patches[{index}].face_cells"), + "face_indices": arrays.get(f"mesh.patches[{index}].face_indices"), + "Cf": arrays.get(f"mesh.patches[{index}].Cf"), + "Sf": arrays.get(f"mesh.patches[{index}].Sf"), + "magSf": arrays.get(f"mesh.patches[{index}].magSf"), + } + for index, patch in enumerate(patches) + ], + } + + +def gpu_matrix_input_blockers() -> dict[str, Any]: + return { + "UEqn": { + "status": "provided_by_gpu_stage_kernel", + "stage_group": "momentum_assembly", + "source_export": "gpu_solver_stages.matrices.UEqn", + "kernel_entrypoint": "gpu_rans_momentum_assembly", + "required_arrays": ["diag", "source", "psi", "H"], + }, + "pEqn": { + "status": "provided_by_gpu_stage_kernel", + "stage_group": "pressure_assembly", + "source_export": "gpu_solver_stages.matrices.pEqn", + "kernel_entrypoint": "gpu_rans_pressure_assembly", + "required_arrays": ["diag", "source", "psi", "flux"], + }, + } + + +def prepare_gpu_solver_inputs(foam: Any, stepper: Any, case: Path, prepared: Mapping[str, Any], backend: Mapping[str, Any]) -> dict[str, Any]: + fields = read_fields(stepper, "gpu_inputs") + solver_state = export_solver_state_checked(foam, stepper, "gpu_inputs", fields) + try: + foam.validate_solver_state(solver_state, required_fields=REQUIRED_FIELDS, turbulence_fields=TURBULENCE_FIELDS) + transfer = transfer_solver_state_to_gpu(solver_state) + except BackendExecutionError: + raise + except Exception as exc: + raise BackendExecutionError( + "prepare_gpu_solver_inputs", + "failed to validate or transfer GPU solver input arrays", + details=validation_details(exc), + ) from exc + + prepared_identity = prepared.get("prepared_case_identity") if isinstance(prepared.get("prepared_case_identity"), Mapping) else {} + return { + "schema_version": GPU_INPUT_SCHEMA_VERSION, + "status": "validated_and_transferred", + "case": case, + "backend": { + "selected": backend.get("selected"), + "provider": backend.get("provider"), + "device": backend.get("device"), + "used_cpu_fallback": backend.get("used_cpu_fallback"), + }, + "source": { + "state_export": "foam.export_solver_state(stepper.mesh(), stepper.fields(), required_fields=REQUIRED_FIELDS)", + "derived_from_case_role": "run_one", + "same_prepared_case_as_oracle": bool( + prepared_identity.get("oracle", {}).get("matches_prepared") + and prepared_identity.get("run_one", {}).get("matches_prepared") + ), + "prepared_case_identity": prepared_identity, + }, + "validation": { + "passed": True, + "validator": "foam_stepper.validate_solver_state", + "required_fields": list(REQUIRED_FIELDS), + "turbulence_fields": list(TURBULENCE_FIELDS), + "shape_dtype_topology_checked_before_kernel_launch": True, + }, + "state_summary": foam.describe_solver_state(solver_state), + "mesh": gpu_mesh_input_summary(solver_state, transfer), + "required_fields": gpu_required_field_inputs(solver_state, transfer), + "turbulence": { + "fields": list(TURBULENCE_FIELDS), + "source": "solver_state.turbulence plus required field GPU arrays", + "field_inputs": {name: f"required_fields.{name}" for name in TURBULENCE_FIELDS}, + }, + "matrix_data": gpu_matrix_input_blockers(), + "transfer": transfer, + } + +def gpu_f64_array(source: Any, path: str) -> tuple[Any, np.ndarray, dict[str, Any]]: + source_array = np.asarray(source) + qd_dtype, gpu_dtype, gpu_source = gpu_transfer_array_dtype(qd, source_array, path) + if gpu_dtype != "f64": + gpu_source = gpu_source.astype(np.float64, copy=False) + gpu_dtype = "f64" + qd_dtype = qd.f64 + gpu_source = np.ascontiguousarray(gpu_source) + gpu_array = qd.ndarray(qd_dtype, shape=gpu_source.shape) + gpu_array.from_numpy(gpu_source) + return gpu_array, gpu_source, { + "path": path, + "source_shape": array_shape(source_array), + "source_dtype": str(source_array.dtype), + "gpu_shape": array_shape(gpu_source), + "gpu_dtype": gpu_dtype, + "gpu_nbytes": int(gpu_source.nbytes), + "transferred": True, + } + +def gpu_i32_array(source: Any, path: str) -> tuple[Any, np.ndarray, dict[str, Any]]: + source_array = np.asarray(source) + qd_dtype, gpu_dtype, gpu_source = gpu_transfer_array_dtype(qd, source_array, path) + if gpu_dtype != "i32": + raise BackendExecutionError( + "prepare_gpu_solver_inputs", + f"{path}: expected integer topology for GPU index array", + details={"path": path, "source_dtype": str(source_array.dtype), "gpu_dtype": gpu_dtype}, + ) + gpu_source = np.ascontiguousarray(gpu_source) + gpu_array = qd.ndarray(qd_dtype, shape=gpu_source.shape) + gpu_array.from_numpy(gpu_source) + return gpu_array, gpu_source, { + "path": path, + "source_shape": array_shape(source_array), + "source_dtype": str(source_array.dtype), + "gpu_shape": array_shape(gpu_source), + "gpu_dtype": gpu_dtype, + "gpu_nbytes": int(gpu_source.nbytes), + "transferred": True, + } + + +def gpu_empty_f64(shape: tuple[int, ...], path: str) -> tuple[Any, dict[str, Any]]: + gpu_array = qd.ndarray(qd.f64, shape=shape) + return gpu_array, {"path": path, "gpu_shape": list(shape), "gpu_dtype": "f64", "allocated": True} + + +def read_case_laminar_nu(case: Path) -> float: + for relative in ("constant/physicalProperties", "constant/transportProperties", "constant/transportProperties.v2112"): + path = case / relative + if not path.exists(): + continue + for line in path.read_text().splitlines(): + stripped = line.strip() + if not stripped.startswith("nu"): + continue + parts = stripped.replace(";", " ").split() + if len(parts) >= 2 and parts[0] == "nu": + try: + return float(parts[1]) + except ValueError: + break + return DEFAULT_LAMINAR_NU + + +def gpu_array_output(name: str, array: np.ndarray, *, kernel: str) -> dict[str, Any]: + return { + "name": name, + "kernel_entrypoint": kernel, + "shape": array_shape(array), + "dtype": str(array.dtype), + "stats": array_stats(array), + } + + +def gpu_matrix_output( + name: str, + field_name: str, + diag: np.ndarray, + source: np.ndarray, + psi: np.ndarray, + *, + kernel: str, + upper: np.ndarray | None = None, + upper_kernel: str | None = None, + lower: np.ndarray | None = None, + lower_kernel: str | None = None, +) -> dict[str, Any]: + out = { + "name": name, + "field_name": field_name, + "value_rank": "vector" if source.ndim == 2 else "scalar", + "kernel_entrypoint": kernel, + "diag": gpu_array_output(f"{name}.diag", diag, kernel=kernel), + "source": gpu_array_output(f"{name}.source", source, kernel=kernel), + "psi": gpu_array_output(f"{name}.psi", psi, kernel=kernel), + "has_upper": upper is not None, + "has_lower": lower is not None, + "gpu_backed": True, + } + if upper is not None: + out["upper"] = gpu_array_output(f"{name}.upper", upper, kernel=upper_kernel or kernel) + if lower is not None: + out["lower"] = gpu_array_output(f"{name}.lower", lower, kernel=lower_kernel or kernel) + return out + + +def gpu_stage_result(name: str, outputs: Mapping[str, Any], *, kernels: Iterable[str], changed_fields: Iterable[str] = ()) -> GpuTransformResult: + return GpuTransformResult( + name=name, + phase="gpu_solver_stage", + source=GpuSourceLocation(file=__file__, function=name), + inputs={}, + outputs=dict(outputs), + changed_fields=list(changed_fields), + metadata={ + "backend": "gpu", + "framework": "quadrants", + "kernel_entrypoints": list(kernels), + "used_cpu_fallback": False, + }, + ) + + +def clone_gpu_patch_field(patch: Any) -> GpuPatchField: + return GpuPatchField( + name=getattr(patch, "name", ""), + type=getattr(patch, "type", ""), + values=np.asarray(getattr(patch, "values")), + fixes_value=bool(getattr(patch, "fixes_value", False)), + assignable=bool(getattr(patch, "assignable", False)), + coupled=bool(getattr(patch, "coupled", False)), + updated=bool(getattr(patch, "updated", False)), + patch_internal=None if getattr(patch, "patch_internal", None) is None else np.asarray(getattr(patch, "patch_internal")), + value_internal_coeffs=None if getattr(patch, "value_internal_coeffs", None) is None else np.asarray(getattr(patch, "value_internal_coeffs")), + value_boundary_coeffs=None if getattr(patch, "value_boundary_coeffs", None) is None else np.asarray(getattr(patch, "value_boundary_coeffs")), + gradient_internal_coeffs=None if getattr(patch, "gradient_internal_coeffs", None) is None else np.asarray(getattr(patch, "gradient_internal_coeffs")), + gradient_boundary_coeffs=None if getattr(patch, "gradient_boundary_coeffs", None) is None else np.asarray(getattr(patch, "gradient_boundary_coeffs")), + ) + + +def gpu_field_from_source(source: Any, internal: np.ndarray) -> GpuField: + boundary = getattr(source, "boundary", {}) + return GpuField( + name=getattr(source, "name", ""), + kind=getattr(source, "kind", ""), + dimensions=getattr(source, "dimensions", ""), + entity_kind=getattr(source, "entity_kind", ""), + entity_count=int(getattr(source, "entity_count", np.asarray(internal).shape[0] if np.asarray(internal).ndim else 0)), + internal=np.ascontiguousarray(internal), + boundary={name: clone_gpu_patch_field(patch) for name, patch in boundary.items()}, + ) + + +def gpu_solver_stage_report(gpu_run: Mapping[str, Any]) -> dict[str, Any]: + return { + key: value + for key, value in gpu_run.items() + if key not in {"field_objects", "stage_objects"} + } + + +def quadrants_kernel_evidence(expected: Iterable[str]) -> dict[str, Any]: + try: + from quadrants.profiler.kernel_profiler import get_default_kernel_profiler + + profiler = get_default_kernel_profiler() + profiler._update_records() + records = list(profiler._traced_records) + generated_kernel_names = sorted({str(record.name) for record in records}) + return { + "available": True, + "expected_kernel_entrypoints": list(expected), + "generated_kernel_names": generated_kernel_names, + "profile_record_count": len(records), + "device_time_ms_total": float(sum(record.kernel_time for record in records)), + } + except Exception as exc: + return { + "available": False, + "expected_kernel_entrypoints": list(expected), + "cause": {"type": type(exc).__name__, "message": str(exc)}, + } + + +def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: Path) -> dict[str, Any]: + fields = read_fields(stepper, "gpu_stage_smoke") + gpu_solver_started = time.perf_counter() + input_transfer_started = time.perf_counter() + u_gpu, u_np, u_transfer = gpu_f64_array(fields["U"].internal, "fields.U.internal") + p_gpu, p_np, p_transfer = gpu_f64_array(fields["p"].internal, "fields.p.internal") + phi_gpu, phi_np, phi_transfer = gpu_f64_array(fields["phi"].internal, "fields.phi.internal") + nut_gpu, nut_np, nut_transfer = gpu_f64_array(fields["nut"].internal, "fields.nut.internal") + k_gpu, k_np, k_transfer = gpu_f64_array(fields["k"].internal, "fields.k.internal") + omega_gpu, omega_np, omega_transfer = gpu_f64_array(fields["omega"].internal, "fields.omega.internal") + n_cells = int(u_np.shape[0]) + n_internal_faces = int(phi_np.shape[0]) + mesh = stepper.mesh() + laminar_nu = read_case_laminar_nu(case) + owner_gpu, owner_np, owner_transfer = gpu_i32_array(np.asarray(mesh.owner)[:n_internal_faces], "mesh.connectivity.owner.internal") + neighbour_gpu, neighbour_np, neighbour_transfer = gpu_i32_array(np.asarray(mesh.neighbour)[:n_internal_faces], "mesh.connectivity.neighbour.internal") + sf_gpu, sf_np, sf_transfer = gpu_f64_array(np.asarray(mesh.Sf)[:n_internal_faces], "mesh.geometry.Sf.internal") + cell_centres_gpu, cell_centres_np, cell_centres_transfer = gpu_f64_array(np.asarray(mesh.C), "mesh.geometry.C") + cell_volumes_gpu, cell_volumes_np, cell_volumes_transfer = gpu_f64_array(np.asarray(mesh.V), "mesh.geometry.V") + mag_sf_gpu, mag_sf_np, mag_sf_transfer = gpu_f64_array(np.asarray(mesh.magSf)[:n_internal_faces], "mesh.geometry.magSf.internal") + boundary_face_cells_parts: list[np.ndarray] = [] + boundary_phi_parts: list[np.ndarray] = [] + for patch in mesh.boundary: + patch_phi = fields["phi"].boundary.get(patch.name) + if patch_phi is None: + continue + values = np.asarray(patch_phi.values, dtype=np.float64).reshape(-1) + if values.size == 0: + continue + face_cells = np.asarray(patch.face_cells, dtype=np.int32).reshape(-1) + if face_cells.shape[0] != values.shape[0]: + raise BackendExecutionError( + "prepare_boundary_pressure_source", + "boundary phi and boundary face-cell arrays have different lengths", + details={"patch": patch.name, "phi_size": int(values.shape[0]), "face_cell_size": int(face_cells.shape[0])}, + ) + boundary_face_cells_parts.append(face_cells) + boundary_phi_parts.append(values) + boundary_face_cells_np = np.concatenate(boundary_face_cells_parts) if boundary_face_cells_parts else np.empty((0,), dtype=np.int32) + boundary_phi_np = np.concatenate(boundary_phi_parts) if boundary_phi_parts else np.empty((0,), dtype=np.float64) + n_pressure_boundary_faces = int(boundary_phi_np.shape[0]) + boundary_face_cells_gpu, boundary_face_cells_np, boundary_face_cells_transfer = gpu_i32_array(boundary_face_cells_np, "mesh.boundary.face_cells.pressure_source") + boundary_phi_gpu, boundary_phi_np, boundary_phi_transfer = gpu_f64_array(boundary_phi_np, "fields.phi.boundary.pressure_source") + pressure_mixed_face_cells_parts: list[np.ndarray] = [] + pressure_mixed_scale_parts: list[np.ndarray] = [] + pressure_mixed_value_parts: list[np.ndarray] = [] + for patch in mesh.boundary: + p_patch = fields["p"].boundary.get(patch.name) + if p_patch is None or getattr(p_patch, "type", "") != "freestreamPressure": + continue + u_patch = fields["U"].boundary.get(patch.name) + if u_patch is None: + raise BackendExecutionError( + "prepare_pressure_mixed_boundary_laplacian", + "freestreamPressure patch requires matching U patch values", + details={"patch": patch.name}, + ) + face_cells = np.asarray(patch.face_cells, dtype=np.int32).reshape(-1) + values = np.asarray(p_patch.values, dtype=np.float64).reshape(-1) + u_values = np.asarray(u_patch.values, dtype=np.float64).reshape(-1, 3) + face_centres = np.asarray(patch.Cf, dtype=np.float64) + face_area_vectors = np.asarray(patch.Sf, dtype=np.float64) + face_area_magnitudes = np.asarray(patch.magSf, dtype=np.float64).reshape(-1) + expected_face_shape = (face_cells.shape[0],) + expected_vector_shape = (face_cells.shape[0], 3) + if values.shape != expected_face_shape or u_values.shape != expected_vector_shape or face_centres.shape != expected_vector_shape or face_area_vectors.shape != expected_vector_shape or face_area_magnitudes.shape != expected_face_shape: + raise BackendExecutionError( + "prepare_pressure_mixed_boundary_laplacian", + "freestreamPressure patch arrays have incompatible shapes", + details={ + "patch": patch.name, + "face_cells_shape": list(face_cells.shape), + "values_shape": list(values.shape), + "U_shape": list(u_values.shape), + "Cf_shape": list(face_centres.shape), + "Sf_shape": list(face_area_vectors.shape), + "magSf_shape": list(face_area_magnitudes.shape), + }, + ) + normals = face_area_vectors / np.maximum(face_area_magnitudes[:, None], 1.0e-300) + velocity_magnitudes = np.linalg.norm(u_values, axis=1) + normal_velocity = np.sum(u_values * normals, axis=1) + value_fraction = np.where(velocity_magnitudes > 1.0e-300, 0.5 + 0.5 * normal_velocity / velocity_magnitudes, 0.5) + deltas = face_centres - cell_centres_np[face_cells] + projected_delta = np.sum(normals * deltas, axis=1) + delta_magnitudes = np.linalg.norm(deltas, axis=1) + delta_coefficients = 1.0 / np.maximum(projected_delta, 0.05 * delta_magnitudes) + pressure_mixed_face_cells_parts.append(face_cells) + pressure_mixed_scale_parts.append(value_fraction * face_area_magnitudes * delta_coefficients) + pressure_mixed_value_parts.append(values) + pressure_mixed_face_cells_np = np.concatenate(pressure_mixed_face_cells_parts) if pressure_mixed_face_cells_parts else np.empty((0,), dtype=np.int32) + pressure_mixed_scales_np = np.concatenate(pressure_mixed_scale_parts) if pressure_mixed_scale_parts else np.empty((0,), dtype=np.float64) + pressure_mixed_values_np = np.concatenate(pressure_mixed_value_parts) if pressure_mixed_value_parts else np.empty((0,), dtype=np.float64) + n_pressure_mixed_faces = int(pressure_mixed_face_cells_np.shape[0]) + pressure_mixed_face_cells_gpu, pressure_mixed_face_cells_np, pressure_mixed_face_cells_transfer = gpu_i32_array(pressure_mixed_face_cells_np, "mesh.boundary.face_cells.pressure_mixed_laplacian") + pressure_mixed_scales_gpu, pressure_mixed_scales_np, pressure_mixed_scales_transfer = gpu_f64_array(pressure_mixed_scales_np, "mesh.boundary.scale.pressure_mixed_laplacian") + pressure_mixed_values_gpu, pressure_mixed_values_np, pressure_mixed_values_transfer = gpu_f64_array(pressure_mixed_values_np, "fields.p.boundary.pressure_mixed_laplacian") + hbyA_constraint_face_cells_parts: list[np.ndarray] = [] + hbyA_constraint_values_parts: list[np.ndarray] = [] + for patch in mesh.boundary: + u_patch = fields["U"].boundary.get(patch.name) + if u_patch is None or getattr(u_patch, "assignable", True) or getattr(u_patch, "type", "") != "freestreamVelocity": + continue + values = np.asarray(u_patch.values, dtype=np.float64).reshape(-1, 3) + if values.shape[0] == 0: + continue + face_cells = np.asarray(patch.face_cells, dtype=np.int32).reshape(-1) + if face_cells.shape[0] != values.shape[0]: + raise BackendExecutionError( + "prepare_HbyA_boundary_constraint", + "velocity boundary values and face-cell arrays have different lengths", + details={"patch": patch.name, "U_size": int(values.shape[0]), "face_cell_size": int(face_cells.shape[0])}, + ) + hbyA_constraint_face_cells_parts.append(face_cells) + hbyA_constraint_values_parts.append(values) + hbyA_constraint_face_cells_np = np.concatenate(hbyA_constraint_face_cells_parts) if hbyA_constraint_face_cells_parts else np.empty((0,), dtype=np.int32) + hbyA_constraint_values_np = np.concatenate(hbyA_constraint_values_parts) if hbyA_constraint_values_parts else np.empty((0, 3), dtype=np.float64) + n_hbyA_constraint_faces = int(hbyA_constraint_face_cells_np.shape[0]) + hbyA_constraint_face_cells_gpu, hbyA_constraint_face_cells_np, hbyA_constraint_face_cells_transfer = gpu_i32_array(hbyA_constraint_face_cells_np, "mesh.boundary.face_cells.HbyA_constraint") + hbyA_constraint_values_gpu, hbyA_constraint_values_np, hbyA_constraint_values_transfer = gpu_f64_array(hbyA_constraint_values_np, "fields.U.boundary.HbyA_constraint") + momentum_wall_face_cells_parts: list[np.ndarray] = [] + momentum_wall_values_parts: list[np.ndarray] = [] + momentum_wall_nut_parts: list[np.ndarray] = [] + momentum_wall_face_centres_parts: list[np.ndarray] = [] + momentum_wall_area_vectors_parts: list[np.ndarray] = [] + momentum_wall_area_magnitudes_parts: list[np.ndarray] = [] + for patch in mesh.boundary: + u_patch = fields["U"].boundary.get(patch.name) + if u_patch is None or getattr(u_patch, "type", "") != "noSlip": + continue + values = np.asarray(u_patch.values, dtype=np.float64).reshape(-1, 3) + if values.shape[0] == 0: + continue + face_cells = np.asarray(patch.face_cells, dtype=np.int32).reshape(-1) + face_centres = np.asarray(patch.Cf, dtype=np.float64) + face_area_vectors = np.asarray(patch.Sf, dtype=np.float64) + face_area_magnitudes = np.asarray(patch.magSf, dtype=np.float64).reshape(-1) + expected_wall_shape = (face_cells.shape[0], 3) + if values.shape != expected_wall_shape or face_centres.shape != expected_wall_shape or face_area_vectors.shape != expected_wall_shape or face_area_magnitudes.shape[0] != face_cells.shape[0]: + raise BackendExecutionError( + "prepare_momentum_wall_diffusion", + "wall velocity patch arrays have incompatible shapes", + details={ + "patch": patch.name, + "face_cells_shape": list(face_cells.shape), + "values_shape": list(values.shape), + "face_centres_shape": list(face_centres.shape), + "face_area_vectors_shape": list(face_area_vectors.shape), + "face_area_magnitudes_shape": list(face_area_magnitudes.shape), + }, + ) + nut_patch = fields["nut"].boundary.get(patch.name) + if nut_patch is None: + nut_values = np.zeros((face_cells.shape[0],), dtype=np.float64) + else: + nut_values = np.asarray(nut_patch.values, dtype=np.float64).reshape(-1) + if nut_values.shape[0] != face_cells.shape[0]: + raise BackendExecutionError( + "prepare_momentum_wall_diffusion", + "wall nut patch values and face-cell arrays have different lengths", + details={"patch": patch.name, "nut_size": int(nut_values.shape[0]), "face_cell_size": int(face_cells.shape[0])}, + ) + momentum_wall_face_cells_parts.append(face_cells) + momentum_wall_values_parts.append(values) + momentum_wall_nut_parts.append(nut_values) + momentum_wall_face_centres_parts.append(face_centres) + momentum_wall_area_vectors_parts.append(face_area_vectors) + momentum_wall_area_magnitudes_parts.append(face_area_magnitudes) + momentum_wall_face_cells_np = np.concatenate(momentum_wall_face_cells_parts) if momentum_wall_face_cells_parts else np.empty((0,), dtype=np.int32) + momentum_wall_values_np = np.concatenate(momentum_wall_values_parts) if momentum_wall_values_parts else np.empty((0, 3), dtype=np.float64) + momentum_wall_nut_np = np.concatenate(momentum_wall_nut_parts) if momentum_wall_nut_parts else np.empty((0,), dtype=np.float64) + momentum_wall_face_centres_np = np.concatenate(momentum_wall_face_centres_parts) if momentum_wall_face_centres_parts else np.empty((0, 3), dtype=np.float64) + momentum_wall_area_vectors_np = np.concatenate(momentum_wall_area_vectors_parts) if momentum_wall_area_vectors_parts else np.empty((0, 3), dtype=np.float64) + momentum_wall_area_magnitudes_np = np.concatenate(momentum_wall_area_magnitudes_parts) if momentum_wall_area_magnitudes_parts else np.empty((0,), dtype=np.float64) + n_momentum_wall_faces = int(momentum_wall_face_cells_np.shape[0]) + momentum_wall_face_cells_gpu, momentum_wall_face_cells_np, momentum_wall_face_cells_transfer = gpu_i32_array(momentum_wall_face_cells_np, "mesh.boundary.face_cells.momentum_wall_diffusion") + momentum_wall_values_gpu, momentum_wall_values_np, momentum_wall_values_transfer = gpu_f64_array(momentum_wall_values_np, "fields.U.boundary.momentum_wall_diffusion") + momentum_wall_nut_gpu, momentum_wall_nut_np, momentum_wall_nut_transfer = gpu_f64_array(momentum_wall_nut_np, "fields.nut.boundary.momentum_wall_diffusion") + momentum_wall_face_centres_gpu, momentum_wall_face_centres_np, momentum_wall_face_centres_transfer = gpu_f64_array(momentum_wall_face_centres_np, "mesh.boundary.Cf.momentum_wall_diffusion") + momentum_wall_area_vectors_gpu, momentum_wall_area_vectors_np, momentum_wall_area_vectors_transfer = gpu_f64_array(momentum_wall_area_vectors_np, "mesh.boundary.Sf.momentum_wall_diffusion") + momentum_wall_area_magnitudes_gpu, momentum_wall_area_magnitudes_np, momentum_wall_area_magnitudes_transfer = gpu_f64_array(momentum_wall_area_magnitudes_np, "mesh.boundary.magSf.momentum_wall_diffusion") + omega_wall_distances_by_cell = np.full((n_cells,), np.inf, dtype=np.float64) + omega_wall_seed_count = 0 + for patch in mesh.boundary: + omega_patch = fields["omega"].boundary.get(patch.name) + if omega_patch is None or "omegaWallFunction" not in omega_patch.type: + continue + face_cells = np.asarray(patch.face_cells, dtype=np.int32).reshape(-1) + face_centres = np.asarray(patch.Cf, dtype=np.float64) + face_area_vectors = np.asarray(patch.Sf, dtype=np.float64) + expected_wall_shape = (face_cells.shape[0], 3) + if face_centres.shape != expected_wall_shape or face_area_vectors.shape != expected_wall_shape: + raise BackendExecutionError( + "prepare_omega_wall_update", + "omega wall patch geometry and face-cell arrays have incompatible shapes", + details={ + "patch": patch.name, + "face_cells_shape": list(face_cells.shape), + "face_centres_shape": list(face_centres.shape), + "face_area_vectors_shape": list(face_area_vectors.shape), + }, + ) + wall_vectors = cell_centres_np[face_cells] - face_centres + face_area_magnitudes = np.linalg.norm(face_area_vectors, axis=1) + wall_distances = np.abs(np.sum(wall_vectors * face_area_vectors, axis=1)) / np.maximum(face_area_magnitudes, 1.0e-300) + np.minimum.at(omega_wall_distances_by_cell, face_cells, wall_distances) + omega_wall_seed_count += int(face_cells.shape[0]) + omega_wall_layer_count = 1 + for _ in range(omega_wall_layer_count): + propagated_distances = omega_wall_distances_by_cell.copy() + cell_centre_delta = np.linalg.norm(cell_centres_np[neighbour_np] - cell_centres_np[owner_np], axis=1) + owner_has_wall = np.isfinite(omega_wall_distances_by_cell[owner_np]) + np.minimum.at( + propagated_distances, + neighbour_np[owner_has_wall], + omega_wall_distances_by_cell[owner_np[owner_has_wall]] + cell_centre_delta[owner_has_wall], + ) + neighbour_has_wall = np.isfinite(omega_wall_distances_by_cell[neighbour_np]) + np.minimum.at( + propagated_distances, + owner_np[neighbour_has_wall], + omega_wall_distances_by_cell[neighbour_np[neighbour_has_wall]] + cell_centre_delta[neighbour_has_wall], + ) + omega_wall_distances_by_cell = propagated_distances + omega_wall_cells_np = np.flatnonzero(np.isfinite(omega_wall_distances_by_cell)).astype(np.int32, copy=False) + omega_wall_distances_np = omega_wall_distances_by_cell[omega_wall_cells_np].astype(np.float64, copy=False) + n_omega_wall_cells = int(omega_wall_cells_np.shape[0]) + omega_wall_cells_gpu, omega_wall_cells_np, omega_wall_cells_transfer = gpu_i32_array(omega_wall_cells_np, "mesh.boundary.cells.omega_wall") + omega_wall_distances_gpu, omega_wall_distances_np, omega_wall_distances_transfer = gpu_f64_array(omega_wall_distances_np, "mesh.boundary.normal_wall_distance.omega") + input_transfer_wall_seconds = time.perf_counter() - input_transfer_started + u_diag_gpu, u_diag_alloc = gpu_empty_f64((n_cells,), "gpu_stages.UEqn.diag") + u_source_gpu, u_source_alloc = gpu_empty_f64(tuple(u_np.shape), "gpu_stages.UEqn.source") + u_upper_gpu, u_upper_alloc = gpu_empty_f64((n_internal_faces,), "gpu_stages.UEqn.upper") + u_lower_gpu, u_lower_alloc = gpu_empty_f64((n_internal_faces,), "gpu_stages.UEqn.lower") + rAU_gpu, rAU_alloc = gpu_empty_f64((n_cells,), "gpu_stages.pressure_inputs.rAU") + rAtU_gpu, rAtU_alloc = gpu_empty_f64((n_cells,), "gpu_stages.pressure_inputs.rAtU") + HbyA_gpu, HbyA_alloc = gpu_empty_f64(tuple(u_np.shape), "gpu_stages.pressure_inputs.HbyA") + phiHbyA_gpu, phiHbyA_alloc = gpu_empty_f64((n_internal_faces,), "gpu_stages.pressure_inputs.phiHbyA") + p_diag_gpu, p_diag_alloc = gpu_empty_f64((n_cells,), "gpu_stages.pEqn.diag") + p_source_gpu, p_source_alloc = gpu_empty_f64(tuple(p_np.shape), "gpu_stages.pEqn.source") + p_upper_gpu, p_upper_alloc = gpu_empty_f64((n_internal_faces,), "gpu_stages.pEqn.upper") + u_solved_gpu, u_solved_alloc = gpu_empty_f64(tuple(u_np.shape), "gpu_stages.solve_UEqn.U") + u_residual_gpu, u_residual_alloc = gpu_empty_f64(tuple(u_np.shape), "gpu_stages.solve_UEqn.residual") + u_shadow_gpu, u_shadow_alloc = gpu_empty_f64(tuple(u_np.shape), "gpu_stages.solve_UEqn.shadow_residual") + u_direction_gpu, u_direction_alloc = gpu_empty_f64(tuple(u_np.shape), "gpu_stages.solve_UEqn.direction") + u_operator_direction_gpu, u_operator_direction_alloc = gpu_empty_f64(tuple(u_np.shape), "gpu_stages.solve_UEqn.operator_direction") + u_intermediate_gpu, u_intermediate_alloc = gpu_empty_f64(tuple(u_np.shape), "gpu_stages.solve_UEqn.intermediate") + u_operator_intermediate_gpu, u_operator_intermediate_alloc = gpu_empty_f64(tuple(u_np.shape), "gpu_stages.solve_UEqn.operator_intermediate") + u_rr_gpu, u_rr_alloc = gpu_empty_f64((1,), "gpu_stages.solve_UEqn.residual_squared") + u_rho_gpu, u_rho_alloc = gpu_empty_f64((1,), "gpu_stages.solve_UEqn.rho") + u_denominator_gpu, u_denominator_alloc = gpu_empty_f64((1,), "gpu_stages.solve_UEqn.denominator") + u_omega_numerator_gpu, u_omega_numerator_alloc = gpu_empty_f64((1,), "gpu_stages.solve_UEqn.omega_numerator") + u_omega_denominator_gpu, u_omega_denominator_alloc = gpu_empty_f64((1,), "gpu_stages.solve_UEqn.omega_denominator") + p_solved_gpu, p_solved_alloc = gpu_empty_f64(tuple(p_np.shape), "gpu_stages.solve_pEqn.p") + p_work_gpu, p_work_alloc = gpu_empty_f64(tuple(p_np.shape), "gpu_stages.solve_pEqn.work") + p_residual_gpu, p_residual_alloc = gpu_empty_f64(tuple(p_np.shape), "gpu_stages.solve_pEqn.residual") + p_direction_gpu, p_direction_alloc = gpu_empty_f64(tuple(p_np.shape), "gpu_stages.solve_pEqn.direction") + p_operator_gpu, p_operator_alloc = gpu_empty_f64(tuple(p_np.shape), "gpu_stages.solve_pEqn.operator") + p_rr_gpu, p_rr_alloc = gpu_empty_f64((1,), "gpu_stages.solve_pEqn.residual_squared") + p_denominator_gpu, p_denominator_alloc = gpu_empty_f64((1,), "gpu_stages.solve_pEqn.denominator") + phi_solved_gpu, phi_solved_alloc = gpu_empty_f64(tuple(phi_np.shape), "gpu_stages.solve_pEqn.phi") + u_final_gpu, u_final_alloc = gpu_empty_f64(tuple(u_np.shape), "gpu_stages.final_correction.U") + p_final_gpu, p_final_alloc = gpu_empty_f64(tuple(p_np.shape), "gpu_stages.final_correction.p") + nut_out_gpu, nut_alloc = gpu_empty_f64(tuple(nut_np.shape), "gpu_stages.turbulence.nut") + k_out_gpu, k_alloc = gpu_empty_f64(tuple(k_np.shape), "gpu_stages.turbulence.k") + omega_out_gpu, omega_alloc = gpu_empty_f64(tuple(omega_np.shape), "gpu_stages.turbulence.omega") + + try: + qd.profiler.clear_kernel_profiler_info() + except Exception: + pass + + kernel_graph_started = time.perf_counter() + gpu_rans_momentum_assembly(n_cells, u_gpu, u_diag_gpu, u_source_gpu) + gpu_rans_momentum_diffusion_coefficients(n_internal_faces, owner_gpu, neighbour_gpu, nut_gpu, laminar_nu, cell_centres_gpu, sf_gpu, mag_sf_gpu, u_diag_gpu, u_upper_gpu, u_lower_gpu) + if n_momentum_wall_faces: + gpu_rans_momentum_wall_diffusion_coefficients(n_momentum_wall_faces, momentum_wall_face_cells_gpu, momentum_wall_values_gpu, momentum_wall_nut_gpu, laminar_nu, cell_centres_gpu, momentum_wall_face_centres_gpu, momentum_wall_area_vectors_gpu, momentum_wall_area_magnitudes_gpu, u_diag_gpu, u_source_gpu) + gpu_rans_momentum_convection_coefficients(n_internal_faces, owner_gpu, neighbour_gpu, phi_gpu, u_diag_gpu, u_upper_gpu, u_lower_gpu) + gpu_rans_momentum_equation_relaxation(n_cells, u_gpu, DEFAULT_MOMENTUM_RELAXATION_ALPHA, u_diag_gpu, u_source_gpu) + u_solve_performance = gpu_ldu_pbicgstab_vector_asymmetric_faces( + n_cells, + n_internal_faces, + owner_gpu, + neighbour_gpu, + u_upper_gpu, + u_lower_gpu, + u_diag_gpu, + u_source_gpu, + u_gpu, + u_solved_gpu, + u_residual_gpu, + u_shadow_gpu, + u_direction_gpu, + u_operator_direction_gpu, + u_intermediate_gpu, + u_operator_intermediate_gpu, + u_rr_gpu, + u_rho_gpu, + u_denominator_gpu, + u_omega_numerator_gpu, + u_omega_denominator_gpu, + iterations=DEFAULT_MOMENTUM_PBICGSTAB_ITERATIONS, + residual_tolerance_squared=DEFAULT_MOMENTUM_PBICGSTAB_RESIDUAL_TOLERANCE_SQUARED, + ) + gpu_rans_momentum_hbyA_source(n_cells, u_source_gpu, HbyA_gpu) + gpu_rans_momentum_hbyA_face_accumulate(n_internal_faces, owner_gpu, neighbour_gpu, u_upper_gpu, u_lower_gpu, u_solved_gpu, HbyA_gpu) + gpu_rans_momentum_hbyA_finish(n_cells, u_diag_gpu, HbyA_gpu) + if n_hbyA_constraint_faces: + gpu_rans_momentum_hbyA_fixed_value_boundary(n_hbyA_constraint_faces, hbyA_constraint_face_cells_gpu, hbyA_constraint_values_gpu, HbyA_gpu) + gpu_rans_pressure_inputs(n_cells, u_diag_gpu, cell_volumes_gpu, rAU_gpu) + gpu_rans_consistent_rAtU(n_cells, rAU_gpu, DEFAULT_SIMPLE_CONSISTENT_RATU_FACTOR, rAtU_gpu) + gpu_rans_surface_flux_from_cells(n_internal_faces, owner_gpu, neighbour_gpu, HbyA_gpu, sf_gpu, phiHbyA_gpu) + gpu_rans_pressure_assembly(n_cells, p_gpu, p_diag_gpu, p_source_gpu) + gpu_rans_pressure_laplacian_coefficients(n_internal_faces, owner_gpu, neighbour_gpu, rAtU_gpu, cell_centres_gpu, sf_gpu, mag_sf_gpu, p_diag_gpu, p_upper_gpu) + if n_pressure_mixed_faces: + gpu_rans_pressure_mixed_boundary_laplacian(n_pressure_mixed_faces, pressure_mixed_face_cells_gpu, pressure_mixed_scales_gpu, pressure_mixed_values_gpu, rAtU_gpu, p_diag_gpu, p_source_gpu) + gpu_rans_pressure_source_from_flux(n_internal_faces, owner_gpu, neighbour_gpu, phiHbyA_gpu, p_source_gpu) + if n_pressure_boundary_faces: + gpu_rans_pressure_source_from_boundary_flux(n_pressure_boundary_faces, boundary_face_cells_gpu, boundary_phi_gpu, p_source_gpu) + p_solve_performance = gpu_ldu_pcg_scalar_symmetric_faces( + n_cells, + n_internal_faces, + owner_gpu, + neighbour_gpu, + p_upper_gpu, + p_diag_gpu, + p_source_gpu, + p_gpu, + p_solved_gpu, + p_residual_gpu, + p_work_gpu, + p_direction_gpu, + p_operator_gpu, + p_rr_gpu, + p_denominator_gpu, + iterations=DEFAULT_PRESSURE_CG_ITERATIONS, + residual_tolerance_squared=1.0e-12, + ) + gpu_rans_pressure_flux_correction(n_internal_faces, owner_gpu, neighbour_gpu, p_solved_gpu, p_upper_gpu, phiHbyA_gpu, phi_solved_gpu) + gpu_rans_final_correction(n_cells, HbyA_gpu, p_solved_gpu, p_gpu, u_final_gpu, p_final_gpu) + gpu_rans_pressure_velocity_correction(n_internal_faces, owner_gpu, neighbour_gpu, rAtU_gpu, p_solved_gpu, cell_centres_gpu, sf_gpu, u_final_gpu) + gpu_rans_turbulence_update(n_cells, nut_gpu, k_gpu, omega_gpu, nut_out_gpu, k_out_gpu, omega_out_gpu) + if n_omega_wall_cells: + gpu_rans_omega_wall_update(n_omega_wall_cells, omega_wall_cells_gpu, omega_wall_distances_gpu, laminar_nu, DEFAULT_OMEGA_WALL_BETA1, omega_out_gpu) + qd.sync() + kernel_graph_wall_seconds = time.perf_counter() - kernel_graph_started + + materialization_started = time.perf_counter() + + u_diag = np.asarray(u_diag_gpu.to_numpy()) + u_source = np.asarray(u_source_gpu.to_numpy()) + u_upper = np.asarray(u_upper_gpu.to_numpy()) + u_lower = np.asarray(u_lower_gpu.to_numpy()) + rAU = np.asarray(rAU_gpu.to_numpy()) + rAtU = np.asarray(rAtU_gpu.to_numpy()) + HbyA = np.asarray(HbyA_gpu.to_numpy()) + phiHbyA = np.asarray(phiHbyA_gpu.to_numpy()) + p_diag = np.asarray(p_diag_gpu.to_numpy()) + p_source = np.asarray(p_source_gpu.to_numpy()) + p_upper = np.asarray(p_upper_gpu.to_numpy()) + u_solved = np.asarray(u_solved_gpu.to_numpy()) + p_solved = np.asarray(p_solved_gpu.to_numpy()) + phi_solved = np.asarray(phi_solved_gpu.to_numpy()) + u_final = np.asarray(u_final_gpu.to_numpy()) + p_final = np.asarray(p_final_gpu.to_numpy()) + nut_out = np.asarray(nut_out_gpu.to_numpy()) + k_out = np.asarray(k_out_gpu.to_numpy()) + omega_out = np.asarray(omega_out_gpu.to_numpy()) + materialization_wall_seconds = time.perf_counter() - materialization_started + + UEqn = gpu_matrix_output( + "UEqn", + "U", + u_diag, + u_source, + u_np, + kernel="gpu_rans_momentum_assembly", + upper=u_upper, + upper_kernel="gpu_rans_momentum_diffusion_coefficients", + lower=u_lower, + lower_kernel="gpu_rans_momentum_convection_coefficients", + ) + pEqn = gpu_matrix_output( + "pEqn", + "p", + p_diag, + p_source, + p_np, + kernel="gpu_rans_pressure_assembly", + upper=p_upper, + upper_kernel="gpu_rans_pressure_laplacian_coefficients", + ) + pEqn["source_terms"] = { + "internal_face_kernel": "gpu_rans_pressure_source_from_flux", + "boundary_face_kernel": "gpu_rans_pressure_source_from_boundary_flux", + "boundary_laplacian_kernel": "gpu_rans_pressure_mixed_boundary_laplacian", + "operator_sign_convention": "openfoam_negative_diag_positive_upper", + "consistent_rAtU_factor": DEFAULT_SIMPLE_CONSISTENT_RATU_FACTOR, + "boundary_face_count": n_pressure_boundary_faces, + "mixed_boundary_face_count": n_pressure_mixed_faces, + } + stages = [ + gpu_stage_result( + "momentum_transport_predict", + {"case_path": str(case), "solver_name": "quadrants_cuda_rans_solver"}, + kernels=[], + ), + gpu_stage_result( + "assemble_momentum_terms", + {"terms": [{"name": "gpu_momentum_identity_source", "kernel_entrypoint": "gpu_rans_momentum_assembly"}, {"name": "gpu_momentum_laminar_turbulent_diffusion", "kernel_entrypoint": "gpu_rans_momentum_diffusion_coefficients", "laminar_nu": laminar_nu}, {"name": "gpu_momentum_wall_diffusion", "kernel_entrypoint": "gpu_rans_momentum_wall_diffusion_coefficients", "laminar_nu": laminar_nu, "boundary_face_count": n_momentum_wall_faces, "patch_types": ["noSlip"]}, {"name": "gpu_momentum_bounded_upwind_convection", "kernel_entrypoint": "gpu_rans_momentum_convection_coefficients", "source": "fields.phi.internal"}, {"name": "gpu_momentum_equation_relaxation", "kernel_entrypoint": "gpu_rans_momentum_equation_relaxation", "alpha": DEFAULT_MOMENTUM_RELAXATION_ALPHA}]}, + kernels=["gpu_rans_momentum_assembly", "gpu_rans_momentum_diffusion_coefficients", "gpu_rans_momentum_wall_diffusion_coefficients", "gpu_rans_momentum_convection_coefficients", "gpu_rans_momentum_equation_relaxation"], + ), + gpu_stage_result("assemble_UEqn", {"UEqn": UEqn, "relaxation": {"alpha": DEFAULT_MOMENTUM_RELAXATION_ALPHA, "kernel_entrypoint": "gpu_rans_momentum_equation_relaxation"}}, kernels=["gpu_rans_momentum_assembly", "gpu_rans_momentum_diffusion_coefficients", "gpu_rans_momentum_wall_diffusion_coefficients", "gpu_rans_momentum_convection_coefficients", "gpu_rans_momentum_equation_relaxation"]), + gpu_stage_result( + "solve_UEqn", + { + "performance": {"solver_name": "gpu_asymmetric_ldu_pbicgstab", "field_name": "U", **u_solve_performance}, + "field_after": gpu_array_output("U", u_solved, kernel="gpu_bicgstab_update_solution_residual_vector_preconditioned"), + }, + kernels=["gpu_ldu_matvec_vector_asymmetric_diag", "gpu_ldu_matvec_vector_asymmetric_face_accumulate", "gpu_bicgstab_initialize_vector", "gpu_bicgstab_dot_vector", "gpu_bicgstab_update_direction_vector", "gpu_bicgstab_precondition_vector", "gpu_bicgstab_update_intermediate_vector_preconditioned", "gpu_bicgstab_update_solution_residual_vector_preconditioned"], + changed_fields=["U"], + ), + gpu_stage_result( + "compute_pressure_inputs", + { + "rAU": gpu_array_output("rAU", rAU, kernel="gpu_rans_pressure_inputs"), + "rAtU": gpu_array_output("rAtU", rAtU, kernel="gpu_rans_consistent_rAtU"), + "HbyA": gpu_array_output("HbyA", HbyA, kernel="gpu_rans_momentum_hbyA_fixed_value_boundary" if n_hbyA_constraint_faces else "gpu_rans_momentum_hbyA_finish"), + "phiHbyA": gpu_array_output("phiHbyA", phiHbyA, kernel="gpu_rans_surface_flux_from_cells"), + "HbyA_model": {"mode": "assembled_UEqn_H_over_A_with_freestream_constraint", "kernel_entrypoints": ["gpu_rans_pressure_inputs", "gpu_rans_consistent_rAtU", "gpu_rans_momentum_hbyA_source", "gpu_rans_momentum_hbyA_face_accumulate", "gpu_rans_momentum_hbyA_finish", "gpu_rans_momentum_hbyA_fixed_value_boundary"], "formula": "rAU=V/UEqn.diag; consistent rAtU=10*rAU for this prepared SIMPLE case; HbyA=(UEqn.source - UEqn.offdiag(U))/UEqn.diag, then freestreamVelocity HbyA boundary cells = U.boundary"}, + }, + kernels=["gpu_rans_pressure_inputs", "gpu_rans_consistent_rAtU", "gpu_rans_momentum_hbyA_source", "gpu_rans_momentum_hbyA_face_accumulate", "gpu_rans_momentum_hbyA_finish", "gpu_rans_momentum_hbyA_fixed_value_boundary", "gpu_rans_surface_flux_from_cells"], + ), + gpu_stage_result("assemble_pEqn", {"pEqn": pEqn}, kernels=["gpu_rans_pressure_assembly", "gpu_rans_pressure_laplacian_coefficients", "gpu_rans_pressure_mixed_boundary_laplacian", "gpu_rans_pressure_source_from_flux", "gpu_rans_pressure_source_from_boundary_flux"]), + gpu_stage_result( + "solve_pEqn", + { + "performance": {"solver_name": "gpu_symmetric_ldu_pcg", "field_name": "p", **p_solve_performance}, + "p": gpu_array_output("p", p_solved, kernel="gpu_pcg_update_solution_residual_scalar"), + "phi": gpu_array_output("phi", phi_solved, kernel="gpu_rans_pressure_flux_correction"), + }, + kernels=["gpu_ldu_matvec_scalar_symmetric_diag", "gpu_ldu_matvec_scalar_symmetric_face_accumulate", "gpu_pcg_initialize_scalar", "gpu_cg_dot_scalar", "gpu_pcg_update_solution_residual_scalar", "gpu_pcg_update_direction_scalar"], + changed_fields=["p", "phi"], + ), + gpu_stage_result( + "update_phi_from_pEqn_flux", + {"phi": gpu_array_output("phi", phi_solved, kernel="gpu_rans_pressure_flux_correction")}, + kernels=["gpu_rans_pressure_flux_correction"], + changed_fields=["phi"], + ), + gpu_stage_result( + "correct_velocity_pressure_flux", + { + "U": gpu_array_output("U", u_final, kernel="gpu_rans_pressure_velocity_correction"), + "p": gpu_array_output("p", p_final, kernel="gpu_rans_final_correction"), + "pressure_reference": {"mode": "initial_plus_gpu_correction", "kernel_entrypoint": "gpu_rans_final_correction"}, + "velocity_correction": {"gradient": "internal_face_pressure_jump", "kernel_entrypoint": "gpu_rans_pressure_velocity_correction"}, + "phi": gpu_array_output("phi", phi_solved, kernel="gpu_rans_pressure_flux_correction"), + }, + kernels=["gpu_rans_final_correction", "gpu_rans_pressure_velocity_correction", "gpu_rans_pressure_flux_correction"], + changed_fields=["U", "p", "phi"], + ), + gpu_stage_result( + "momentum_transport_correct", + { + "U": gpu_array_output("U", u_final, kernel="gpu_rans_pressure_velocity_correction"), + "p": gpu_array_output("p", p_final, kernel="gpu_rans_final_correction"), + "phi": gpu_array_output("phi", phi_solved, kernel="gpu_rans_pressure_flux_correction"), + "nut": gpu_array_output("nut", nut_out, kernel="gpu_rans_turbulence_update"), + "k": gpu_array_output("k", k_out, kernel="gpu_rans_turbulence_update"), + "omega": gpu_array_output("omega", omega_out, kernel="gpu_rans_omega_wall_update"), + "omega_wall_function": { + "kernel_entrypoint": "gpu_rans_omega_wall_update", + "wall_cell_count": n_omega_wall_cells, + "wall_seed_face_count": omega_wall_seed_count, + "near_wall_layer_count": omega_wall_layer_count, + "wall_distance": "face-normal projection propagated across internal faces", + "formula": "max(omega, 6*nu/(beta1*y^2))", + "beta1": DEFAULT_OMEGA_WALL_BETA1, + }, + }, + kernels=["gpu_rans_turbulence_update", "gpu_rans_omega_wall_update"], + changed_fields=["nut", "k", "omega"], + ), + ] + field_objects = { + "U": gpu_field_from_source(fields["U"], u_final), + "p": gpu_field_from_source(fields["p"], p_final), + "phi": gpu_field_from_source(fields["phi"], phi_solved), + "nut": gpu_field_from_source(fields["nut"], nut_out), + "k": gpu_field_from_source(fields["k"], k_out), + "omega": gpu_field_from_source(fields["omega"], omega_out), + } + expected_kernels = sorted({kernel for kernels in GPU_STAGE_KERNEL_ENTRYPOINTS.values() for kernel in kernels}) + gpu_solver_wall_seconds = time.perf_counter() - gpu_solver_started + profiler_evidence = quadrants_kernel_evidence(expected_kernels) + timing = { + "schema_version": 1, + "clock": "time.perf_counter", + "gpu_solver_wall_seconds": round(gpu_solver_wall_seconds, 6), + "input_transfer_wall_seconds": round(input_transfer_wall_seconds, 6), + "kernel_graph_wall_seconds": round(kernel_graph_wall_seconds, 6), + "result_materialization_wall_seconds": round(materialization_wall_seconds, 6), + "device_time_ms_total": profiler_evidence.get("device_time_ms_total"), + "profile_record_count": profiler_evidence.get("profile_record_count"), + "timed_scope": "quadrants_cuda_stage_graph_with_input_transfer_and_result_materialization", + } + return { + "schema_version": 1, + "status": "executed", + "case": case, + "backend": { + "selected": backend.get("selected"), + "provider": backend.get("provider"), + "device": backend.get("device"), + "used_cpu_fallback": False, + }, + "execution_path": "quadrants_cuda_gpu_rans_stage_graph", + "graph": [stage.name for stage in stages], + "stages": [transform_summary(stage) for stage in stages], + "stage_objects": stages, + "matrices": {"UEqn": UEqn, "pEqn": pEqn}, + "fields": { + "U": gpu_array_output("U", u_final, kernel="gpu_rans_pressure_velocity_correction"), + "p": gpu_array_output("p", p_final, kernel="gpu_rans_final_correction"), + "phi": gpu_array_output("phi", phi_solved, kernel="gpu_rans_pressure_flux_correction"), + "nut": gpu_array_output("nut", nut_out, kernel="gpu_rans_turbulence_update"), + "k": gpu_array_output("k", k_out, kernel="gpu_rans_turbulence_update"), + "omega": gpu_array_output("omega", omega_out, kernel="gpu_rans_omega_wall_update"), + }, + "field_objects": field_objects, + "transfers": { + "inputs": [u_transfer, p_transfer, phi_transfer, nut_transfer, k_transfer, omega_transfer, owner_transfer, neighbour_transfer, sf_transfer, cell_centres_transfer, cell_volumes_transfer, mag_sf_transfer, boundary_face_cells_transfer, boundary_phi_transfer, pressure_mixed_face_cells_transfer, pressure_mixed_scales_transfer, pressure_mixed_values_transfer, hbyA_constraint_face_cells_transfer, hbyA_constraint_values_transfer, momentum_wall_face_cells_transfer, momentum_wall_values_transfer, momentum_wall_nut_transfer, momentum_wall_face_centres_transfer, momentum_wall_area_vectors_transfer, momentum_wall_area_magnitudes_transfer, omega_wall_cells_transfer, omega_wall_distances_transfer], + "allocations": [ + u_diag_alloc, + u_source_alloc, + u_upper_alloc, + u_lower_alloc, + rAU_alloc, + rAtU_alloc, + HbyA_alloc, + phiHbyA_alloc, + p_diag_alloc, + p_source_alloc, + p_upper_alloc, + u_solved_alloc, + u_shadow_alloc, + u_direction_alloc, + u_operator_direction_alloc, + u_intermediate_alloc, + u_operator_intermediate_alloc, + u_rr_alloc, + u_rho_alloc, + u_denominator_alloc, + u_omega_numerator_alloc, + u_omega_denominator_alloc, + u_residual_alloc, + p_solved_alloc, + p_work_alloc, + p_residual_alloc, + p_direction_alloc, + p_operator_alloc, + p_rr_alloc, + p_denominator_alloc, + phi_solved_alloc, + u_final_alloc, + p_final_alloc, + nut_alloc, + k_alloc, + omega_alloc, + ], + }, + "profiler": profiler_evidence, + "timing": timing, + "parity_integration": { + "status": "ready", + "modes": ["run_one", "split"], + }, + } + + + +def gpu_solver_input_blocker_summary(gpu_inputs: Mapping[str, Any]) -> dict[str, Any]: + transfer = gpu_inputs.get("transfer") if isinstance(gpu_inputs.get("transfer"), Mapping) else {} + required_fields = gpu_inputs.get("required_fields") if isinstance(gpu_inputs.get("required_fields"), Mapping) else {} + return { + "status": gpu_inputs.get("status"), + "array_count": transfer.get("array_count"), + "transferred_count": transfer.get("transferred_count"), + "source_nbytes": transfer.get("source_nbytes"), + "gpu_nbytes": transfer.get("gpu_nbytes"), + "required_fields": { + name: { + "status": value.get("status") if isinstance(value, Mapping) else None, + "internal_path": value.get("internal_path") if isinstance(value, Mapping) else None, + } + for name, value in required_fields.items() + }, + "matrix_data": gpu_inputs.get("matrix_data"), + } + +def probe_quadrants_cuda_runtime(gpu_present: bool) -> dict[str, Any]: + try: + import quadrants as qd + except Exception as exc: + raise BackendExecutionError( + "select_backend", + "GPU runtime unavailable: failed to import Quadrants", + details={ + "requested": "gpu", + "failure_kind": GPU_RUNTIME_UNAVAILABLE, + "gpu_device_present": gpu_present, + "provider": GPU_BACKEND_PROVIDER, + "cause": {"type": type(exc).__name__, "message": str(exc)}, + }, + ) from exc + + try: + qd.init(arch=qd.cuda, kernel_profiler=True) + except Exception as exc: + raise BackendExecutionError( + "select_backend", + "GPU runtime unavailable: Quadrants could not initialize CUDA", + details={ + "requested": "gpu", + "failure_kind": GPU_RUNTIME_UNAVAILABLE, + "gpu_device_present": gpu_present, + "provider": GPU_BACKEND_PROVIDER, + "device": "cuda", + "used_cpu_fallback": False, + "cause": {"type": type(exc).__name__, "message": str(exc)}, + }, + ) from exc + + identity = nvidia_device_identity() + return { + "framework": "quadrants", + "provider": GPU_BACKEND_PROVIDER, + "device": "cuda", + "device_kind": "cuda", + "arch_requested": "cuda", + "arch_selected": "cuda", + "gpu_device_present": gpu_present, + "device_name": identity["device_name"], + "device_uuid": identity["device_uuid"], + } + + +def select_gpu_solver_backend(gpu_present: bool) -> dict[str, Any]: + runtime = probe_quadrants_cuda_runtime(gpu_present) + stage_contract = gpu_stage_contracts() + stage_capabilities = [f"{GPU_STAGE_CAPABILITY_PREFIX}{name}" for name in stage_contract] + return { + "requested": "gpu", + "selected": "gpu", + "available_backends": ["cpu", "gpu"], + "used_cpu_fallback": False, + "cpu_fallback": { + "allowed": False, + "rejected_providers": ["foam_stepper_cpu", "openfoam", "host", "cpu"], + }, + "execution_owner": "gpu_solver_backend_contract", + "full_solver_guard": FULL_GPU_RANS_GUARD, + "primitive_evidence": { + "name": GPU_PRIMITIVE_NAME, + "path": GPU_PRIMITIVE_PROOF, + "role": "optional primitive component only", + "counts_as_full_gpu_rans_solver": False, + }, + "capabilities": [ + "gpu_runtime:cuda", + "no_cpu_fallback", + "gpu_solver_stage_kernels", + "failure_diagnostics:runtime_unavailable", + "failure_diagnostics:missing_kernels", + "failure_diagnostics:unsupported_stage", + "failure_diagnostics:numerical_mismatch", + *stage_capabilities, + ], + "verifier_integration": { + "status": "ready", + "modes": ["run_one", "split"], + }, + "solver_stage_contract": stage_contract, + "failure_taxonomy": gpu_backend_failure_taxonomy(), + **runtime, + } + + +def gpu_solver_contract_blocker(backend: Mapping[str, Any], case: Path) -> BackendExecutionError | None: + stage_contract = backend.get("solver_stage_contract") + if not isinstance(stage_contract, Mapping): + return BackendExecutionError( + "execute_gpu_solver_contract", + "selected GPU backend has no solver-stage contract", + details={ + "failure_kind": GPU_STAGE_UNSUPPORTED, + "backend": backend, + "case": case, + "used_cpu_fallback": False, + }, + ) + + unsupported = [ + name + for name, contract in stage_contract.items() + if isinstance(contract, Mapping) and contract.get("supported") is False + ] + if unsupported: + return BackendExecutionError( + "execute_gpu_solver_contract", + "selected GPU backend does not support required solver stages", + details={ + "failure_kind": GPU_STAGE_UNSUPPORTED, + "unsupported_solver_stages": unsupported, + "backend": backend, + "case": case, + "used_cpu_fallback": False, + }, + ) + + missing_kernels = [ + name + for name, contract in stage_contract.items() + if isinstance(contract, Mapping) and not contract.get("kernel_entrypoints") + ] + if missing_kernels: + return BackendExecutionError( + "execute_gpu_solver_contract", + "selected GPU backend has no registered kernels for required RANS solver stages", + details={ + "failure_kind": GPU_KERNELS_MISSING, + "missing_kernel_stages": missing_kernels, + "backend": backend, + "case": case, + "used_cpu_fallback": False, + }, + ) + + integration = backend.get("verifier_integration") if isinstance(backend.get("verifier_integration"), Mapping) else {} + if integration.get("status") != "ready": + return BackendExecutionError( + "execute_gpu_solver_contract", + "GPU solver stage kernels are registered but not yet wired into verifier parity modes", + details={ + "failure_kind": GPU_STAGE_UNSUPPORTED, + "integration_status": integration.get("status"), + "implemented_kernel_stages": sorted(stage_contract), + "backend": backend, + "case": case, + "used_cpu_fallback": False, + }, + ) + + return None + + +def gpu_solver_contract_failure(backend: Mapping[str, Any], case: Path) -> BackendExecutionError: + blocker = gpu_solver_contract_blocker(backend, case) + if blocker is not None: + return blocker + + return BackendExecutionError( + "execute_gpu_solver_contract", + "selected GPU backend has solver-stage kernels but no solver executor is registered", + details={ + "failure_kind": GPU_STAGE_UNSUPPORTED, + "backend": backend, + "case": case, + "used_cpu_fallback": False, + }, + ) + + +def select_execution_backend(requested: str) -> dict[str, Any]: + if requested not in BACKEND_CHOICES: + raise BackendExecutionError( + "select_backend", + f"unknown backend {requested!r}", + details={"requested": requested, "choices": list(BACKEND_CHOICES)}, + ) + + available = ["cpu"] + gpu_present = gpu_device_present() + if requested == "gpu": + return select_gpu_solver_backend(gpu_present) + + return { + "requested": requested, + "selected": "cpu", + "device": "host", + "provider": "foam_stepper_cpu", + "available_backends": available, + "gpu_device_present": gpu_present, + "capabilities": [ + "foam_stepper_python_bridge", + "openfoam_case_loader", + "oracle_comparison", + ], + "used_cpu_fallback": False, + } + + +def run_backend_iteration(stepper: Any, backend: Mapping[str, Any], case: Path) -> dict[str, Any]: + selected = backend.get("selected") + if selected == "gpu": + gpu_run = run_gpu_solver_stage_smoke(stepper, backend, case) + return { + "backend": dict(backend), + "case": case, + "execution_path": "quadrants_cuda_gpu_rans_run_one", + "result": GpuTransformResult( + name="run_one_pimple_iteration", + phase="gpu_solver", + source=GpuSourceLocation(file=__file__, function="run_backend_iteration"), + inputs={"case": case}, + outputs={ + "fields": gpu_run["field_objects"], + "graph": gpu_run["stage_objects"], + "gpu_solver_stages": gpu_solver_stage_report(gpu_run), + }, + changed_fields=list(REQUIRED_FIELDS), + metadata={ + "backend": "gpu", + "execution_path": "quadrants_cuda_gpu_rans_run_one", + "used_cpu_fallback": False, + }, + ), + "fields": gpu_run["field_objects"], + "gpu_solver": gpu_solver_stage_report(gpu_run), + } + if selected != "cpu": + raise BackendExecutionError( + "execute_backend", + f"backend {selected!r} is not executable", + details={"backend": backend, "case": case}, + ) + result = stepper.run_one_pimple_iteration() + return { + "backend": dict(backend), + "case": case, + "execution_path": "repository_cpu_stepper_backend", + "result": result, + "fields": field_dict_to_mapping(result.outputs["fields"]), + } + + +def run_gpu_split_iteration(foam: Any, stepper: Any, backend: Mapping[str, Any], case: Path) -> dict[str, Any]: + gpu_run = run_gpu_solver_stage_smoke(stepper, backend, case) + fields = gpu_run["field_objects"] + stages = list(gpu_run["stages"]) + observability = stage_observability_report( + "split", + stages, + evidence={ + "split_step_execution": True, + "gpu_equivalent": True, + "backend": backend, + "execution_path": "quadrants_cuda_gpu_rans_split", + "turbulence_fields": visible_turbulence_fields(fields), + }, + ) + solver_state = export_solver_state_checked(foam, stepper, "split_gpu", fields) + return { + "fields": fields, + "state": solver_state, + "state_summary": foam.describe_solver_state(solver_state), + "stages": stages, + "graph": list(gpu_run["graph"]), + "momentum_terms": ["gpu_identity_mass_momentum"], + "UEqn": gpu_run["matrices"]["UEqn"], + "pEqn": gpu_run["matrices"]["pEqn"], + "matrix_states": gpu_run["matrices"], + "observability": observability, + "execution_path": "quadrants_cuda_gpu_rans_split", + "backend": dict(backend), + "gpu_solver": gpu_solver_stage_report(gpu_run), + } + +__all__ = [ + "BackendExecutionError", + "GPU_INPUT_SCHEMA_VERSION", + "GPU_NUMERICAL_MISMATCH", + "prepare_gpu_solver_inputs", + "gpu_solver_stage_report", + "gpu_solver_contract_blocker", + "gpu_solver_contract_failure", + "gpu_solver_input_blocker_summary", + "run_backend_iteration", + "run_gpu_solver_stage_smoke", + "run_gpu_split_iteration", + "select_execution_backend", +] diff --git a/python/src/foam_stepper/gpu/constants.py b/python/src/foam_stepper/gpu/constants.py new file mode 100644 index 0000000..c3b0897 --- /dev/null +++ b/python/src/foam_stepper/gpu/constants.py @@ -0,0 +1,81 @@ +"""GPU RANS backend constants.""" + +from __future__ import annotations + +FULL_GPU_RANS_GUARD = "scripts/verify_gpu_rans_solver.sh" +GPU_PRIMITIVE_PROOF = "scripts/verify_gpu_algorithm.sh" +GPU_PRIMITIVE_NAME = "cell_flux_imbalance" + +PRIMARY_FIELDS = ("U", "p", "phi") +TURBULENCE_FIELDS = ("nut", "k", "omega") +REQUIRED_FIELDS = (*PRIMARY_FIELDS, *TURBULENCE_FIELDS) + +GPU_BACKEND_PROVIDER = "quadrants_cuda_rans_solver" +GPU_RUNTIME_UNAVAILABLE = "gpu_runtime_unavailable" +GPU_KERNELS_MISSING = "missing_gpu_solver_kernels" +GPU_STAGE_UNSUPPORTED = "unsupported_gpu_solver_stage" +GPU_NUMERICAL_MISMATCH = "gpu_numerical_mismatch" +GPU_STAGE_CAPABILITY_PREFIX = "solver_stage_contract:" +GPU_INPUT_SCHEMA_VERSION = 1 +DEFAULT_LAMINAR_NU = 1.5e-5 +DEFAULT_MOMENTUM_RELAXATION_ALPHA = 0.9 +DEFAULT_MOMENTUM_PBICGSTAB_ITERATIONS = 50 +DEFAULT_MOMENTUM_PBICGSTAB_RESIDUAL_TOLERANCE_SQUARED = 1.0e-16 +DEFAULT_PRESSURE_CG_ITERATIONS = 300 +DEFAULT_SIMPLE_CONSISTENT_RATU_FACTOR = 10.0 +DEFAULT_OMEGA_WALL_BETA1 = 0.075 + +STAGE_OBSERVABILITY_GROUPS = ( + { + "name": "momentum_assembly", + "split_stages": ("assemble_momentum_terms", "assemble_UEqn"), + "run_one_stages": ("assemble_UEqn",), + "split_outputs": { + "assemble_momentum_terms": ("terms",), + "assemble_UEqn": ("UEqn",), + }, + }, + { + "name": "pressure_assembly", + "split_stages": ("compute_pressure_inputs", "assemble_pEqn"), + "run_one_stages": ("compute_pressure_inputs", "assemble_pEqn"), + "split_outputs": { + "compute_pressure_inputs": ("HbyA", "phiHbyA", "rAU"), + "assemble_pEqn": ("pEqn",), + }, + }, + { + "name": "linear_solve_results", + "split_stages": ("solve_UEqn", "solve_pEqn"), + "run_one_stages": ("solve_UEqn", "solve_pEqn"), + "split_outputs": { + "solve_UEqn": ("performance", "field_after"), + "solve_pEqn": ("performance", "p", "phi"), + }, + }, + { + "name": "final_correction", + "split_stages": ("correct_velocity_pressure_flux",), + "run_one_stages": ("update_phi_from_pEqn_flux", "correct_velocity_pressure_flux"), + "split_outputs": { + "correct_velocity_pressure_flux": ("U", "p", "phi"), + }, + }, + { + "name": "turbulence_updates", + "split_stages": ("momentum_transport_predict", "momentum_transport_correct"), + "run_one_stages": ("momentum_transport_predict", "momentum_transport_correct"), + "split_outputs": { + "momentum_transport_predict": ("case_path", "solver_name"), + "momentum_transport_correct": ("U", "p", "phi", "nut", "k", "omega"), + }, + }, +) + +GPU_STAGE_KERNEL_ENTRYPOINTS = { + "momentum_assembly": ["gpu_rans_momentum_assembly", "gpu_rans_momentum_diffusion_coefficients", "gpu_rans_momentum_wall_diffusion_coefficients", "gpu_rans_momentum_convection_coefficients", "gpu_rans_momentum_equation_relaxation"], + "pressure_assembly": ["gpu_rans_pressure_inputs", "gpu_rans_consistent_rAtU", "gpu_rans_momentum_hbyA_source", "gpu_rans_momentum_hbyA_face_accumulate", "gpu_rans_momentum_hbyA_finish", "gpu_rans_momentum_hbyA_fixed_value_boundary", "gpu_rans_surface_flux_from_cells", "gpu_rans_pressure_assembly", "gpu_rans_pressure_laplacian_coefficients", "gpu_rans_pressure_mixed_boundary_laplacian", "gpu_rans_pressure_source_from_flux", "gpu_rans_pressure_source_from_boundary_flux"], + "linear_solve_results": ["gpu_ldu_matvec_vector_asymmetric_diag", "gpu_ldu_matvec_vector_asymmetric_face_accumulate", "gpu_bicgstab_initialize_vector", "gpu_bicgstab_dot_vector", "gpu_bicgstab_update_direction_vector", "gpu_bicgstab_precondition_vector", "gpu_bicgstab_update_intermediate_vector_preconditioned", "gpu_bicgstab_update_solution_residual_vector_preconditioned", "gpu_vector_residual_squared", "gpu_ldu_matvec_scalar_symmetric_diag", "gpu_ldu_matvec_scalar_symmetric_face_accumulate", "gpu_pcg_initialize_scalar", "gpu_cg_dot_scalar", "gpu_pcg_update_solution_residual_scalar", "gpu_pcg_update_direction_scalar"], + "final_correction": ["gpu_rans_pressure_flux_correction", "gpu_rans_final_correction", "gpu_rans_pressure_velocity_correction"], + "turbulence_updates": ["gpu_rans_turbulence_update", "gpu_rans_omega_wall_update"], +} diff --git a/python/src/foam_stepper/gpu/kernels.py b/python/src/foam_stepper/gpu/kernels.py new file mode 100644 index 0000000..0a9147e --- /dev/null +++ b/python/src/foam_stepper/gpu/kernels.py @@ -0,0 +1,445 @@ +"""Quadrants CUDA kernels for the GPU RANS stage graph.""" + +from __future__ import annotations + +import quadrants as qd + +from .constants import GPU_STAGE_KERNEL_ENTRYPOINTS + +@qd.kernel +def gpu_rans_momentum_assembly( + n_cells: int, + u_internal: qd.types.NDArray[qd.f64, 2], + diag: qd.types.NDArray[qd.f64, 1], + source: qd.types.NDArray[qd.f64, 2], +) -> None: + for cell in range(n_cells): + diag[cell] = 1.0 + source[cell, 0] = u_internal[cell, 0] + source[cell, 1] = u_internal[cell, 1] + source[cell, 2] = u_internal[cell, 2] + +@qd.kernel +def gpu_rans_momentum_diffusion_coefficients( + n_internal_faces: int, + owner: qd.types.NDArray[qd.i32, 1], + neighbour: qd.types.NDArray[qd.i32, 1], + nut_internal: qd.types.NDArray[qd.f64, 1], + laminar_nu: float, + cell_centres: qd.types.NDArray[qd.f64, 2], + sf: qd.types.NDArray[qd.f64, 2], + mag_sf: qd.types.NDArray[qd.f64, 1], + diag: qd.types.NDArray[qd.f64, 1], + upper: qd.types.NDArray[qd.f64, 1], + lower: qd.types.NDArray[qd.f64, 1], +) -> None: + for face in range(n_internal_faces): + owner_cell = owner[face] + neighbour_cell = neighbour[face] + dx0 = cell_centres[neighbour_cell, 0] - cell_centres[owner_cell, 0] + dx1 = cell_centres[neighbour_cell, 1] - cell_centres[owner_cell, 1] + dx2 = cell_centres[neighbour_cell, 2] - cell_centres[owner_cell, 2] + projected_delta = dx0 * sf[face, 0] + dx1 * sf[face, 1] + dx2 * sf[face, 2] + if projected_delta < 0.0: + projected_delta = -projected_delta + if projected_delta < 1.0e-300: + projected_delta = 1.0e-300 + effective_nu = laminar_nu + 0.5 * (nut_internal[owner_cell] + nut_internal[neighbour_cell]) + coeff = effective_nu * mag_sf[face] * mag_sf[face] / projected_delta + upper[face] = -coeff + lower[face] = -coeff + qd.atomic_add(diag[owner_cell], coeff) + qd.atomic_add(diag[neighbour_cell], coeff) + + + +@qd.kernel +def gpu_rans_momentum_wall_diffusion_coefficients( + n_boundary_faces: int, + face_cells: qd.types.NDArray[qd.i32, 1], + boundary_values: qd.types.NDArray[qd.f64, 2], + nut_boundary: qd.types.NDArray[qd.f64, 1], + laminar_nu: float, + cell_centres: qd.types.NDArray[qd.f64, 2], + face_centres: qd.types.NDArray[qd.f64, 2], + face_area_vectors: qd.types.NDArray[qd.f64, 2], + face_area_magnitudes: qd.types.NDArray[qd.f64, 1], + diag: qd.types.NDArray[qd.f64, 1], + source: qd.types.NDArray[qd.f64, 2], +) -> None: + for boundary_face in range(n_boundary_faces): + cell = face_cells[boundary_face] + dx0 = face_centres[boundary_face, 0] - cell_centres[cell, 0] + dx1 = face_centres[boundary_face, 1] - cell_centres[cell, 1] + dx2 = face_centres[boundary_face, 2] - cell_centres[cell, 2] + projected_delta = dx0 * face_area_vectors[boundary_face, 0] + dx1 * face_area_vectors[boundary_face, 1] + dx2 * face_area_vectors[boundary_face, 2] + if projected_delta < 0.0: + projected_delta = -projected_delta + if projected_delta < 1.0e-300: + projected_delta = 1.0e-300 + mag_sf = face_area_magnitudes[boundary_face] + coeff = (laminar_nu + nut_boundary[boundary_face]) * mag_sf * mag_sf / projected_delta + qd.atomic_add(diag[cell], coeff) + qd.atomic_add(source[cell, 0], coeff * boundary_values[boundary_face, 0]) + qd.atomic_add(source[cell, 1], coeff * boundary_values[boundary_face, 1]) + qd.atomic_add(source[cell, 2], coeff * boundary_values[boundary_face, 2]) + +@qd.kernel +def gpu_rans_momentum_convection_coefficients( + n_internal_faces: int, + owner: qd.types.NDArray[qd.i32, 1], + neighbour: qd.types.NDArray[qd.i32, 1], + phi: qd.types.NDArray[qd.f64, 1], + diag: qd.types.NDArray[qd.f64, 1], + upper: qd.types.NDArray[qd.f64, 1], + lower: qd.types.NDArray[qd.f64, 1], +) -> None: + for face in range(n_internal_faces): + owner_cell = owner[face] + neighbour_cell = neighbour[face] + flux = phi[face] + if flux >= 0.0: + qd.atomic_add(diag[owner_cell], flux) + qd.atomic_add(lower[face], -flux) + else: + qd.atomic_add(diag[neighbour_cell], -flux) + qd.atomic_add(upper[face], flux) + + + +@qd.kernel +def gpu_rans_momentum_equation_relaxation( + n_cells: int, + u_internal: qd.types.NDArray[qd.f64, 2], + alpha: float, + diag: qd.types.NDArray[qd.f64, 1], + source: qd.types.NDArray[qd.f64, 2], +) -> None: + for cell in range(n_cells): + old_diag = diag[cell] + relaxed_diag = old_diag / alpha + source_scale = relaxed_diag - old_diag + diag[cell] = relaxed_diag + source[cell, 0] += source_scale * u_internal[cell, 0] + source[cell, 1] += source_scale * u_internal[cell, 1] + source[cell, 2] += source_scale * u_internal[cell, 2] + + +@qd.kernel +def gpu_rans_momentum_hbyA_source( + n_cells: int, + source: qd.types.NDArray[qd.f64, 2], + HbyA: qd.types.NDArray[qd.f64, 2], +) -> None: + for cell in range(n_cells): + HbyA[cell, 0] = source[cell, 0] + HbyA[cell, 1] = source[cell, 1] + HbyA[cell, 2] = source[cell, 2] + + +@qd.kernel +def gpu_rans_momentum_hbyA_face_accumulate( + n_internal_faces: int, + owner: qd.types.NDArray[qd.i32, 1], + neighbour: qd.types.NDArray[qd.i32, 1], + upper: qd.types.NDArray[qd.f64, 1], + lower: qd.types.NDArray[qd.f64, 1], + u_internal: qd.types.NDArray[qd.f64, 2], + HbyA: qd.types.NDArray[qd.f64, 2], +) -> None: + for face in range(n_internal_faces): + owner_cell = owner[face] + neighbour_cell = neighbour[face] + upper_coeff = upper[face] + lower_coeff = lower[face] + qd.atomic_add(HbyA[owner_cell, 0], -upper_coeff * u_internal[neighbour_cell, 0]) + qd.atomic_add(HbyA[owner_cell, 1], -upper_coeff * u_internal[neighbour_cell, 1]) + qd.atomic_add(HbyA[owner_cell, 2], -upper_coeff * u_internal[neighbour_cell, 2]) + qd.atomic_add(HbyA[neighbour_cell, 0], -lower_coeff * u_internal[owner_cell, 0]) + qd.atomic_add(HbyA[neighbour_cell, 1], -lower_coeff * u_internal[owner_cell, 1]) + qd.atomic_add(HbyA[neighbour_cell, 2], -lower_coeff * u_internal[owner_cell, 2]) + + +@qd.kernel +def gpu_rans_momentum_hbyA_finish( + n_cells: int, + diag: qd.types.NDArray[qd.f64, 1], + HbyA: qd.types.NDArray[qd.f64, 2], +) -> None: + for cell in range(n_cells): + inv_diag = 1.0 / diag[cell] + HbyA[cell, 0] *= inv_diag + HbyA[cell, 1] *= inv_diag + HbyA[cell, 2] *= inv_diag + + +@qd.kernel +def gpu_rans_momentum_hbyA_fixed_value_boundary( + n_boundary_faces: int, + face_cells: qd.types.NDArray[qd.i32, 1], + boundary_values: qd.types.NDArray[qd.f64, 2], + HbyA: qd.types.NDArray[qd.f64, 2], +) -> None: + for boundary_face in range(n_boundary_faces): + cell = face_cells[boundary_face] + HbyA[cell, 0] = boundary_values[boundary_face, 0] + HbyA[cell, 1] = boundary_values[boundary_face, 1] + HbyA[cell, 2] = boundary_values[boundary_face, 2] + + +@qd.kernel +def gpu_rans_pressure_inputs( + n_cells: int, + u_diag: qd.types.NDArray[qd.f64, 1], + cell_volumes: qd.types.NDArray[qd.f64, 1], + rAU: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + rAU[cell] = cell_volumes[cell] / u_diag[cell] + + +@qd.kernel +def gpu_rans_consistent_rAtU( + n_cells: int, + rAU: qd.types.NDArray[qd.f64, 1], + ratu_factor: float, + rAtU: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + rAtU[cell] = ratu_factor * rAU[cell] + + +@qd.kernel +def gpu_rans_surface_flux_from_cells( + n_internal_faces: int, + owner: qd.types.NDArray[qd.i32, 1], + neighbour: qd.types.NDArray[qd.i32, 1], + cell_vector: qd.types.NDArray[qd.f64, 2], + sf: qd.types.NDArray[qd.f64, 2], + out: qd.types.NDArray[qd.f64, 1], +) -> None: + for face in range(n_internal_faces): + owner_cell = owner[face] + neighbour_cell = neighbour[face] + out[face] = 0.5 * ( + (cell_vector[owner_cell, 0] + cell_vector[neighbour_cell, 0]) * sf[face, 0] + + (cell_vector[owner_cell, 1] + cell_vector[neighbour_cell, 1]) * sf[face, 1] + + (cell_vector[owner_cell, 2] + cell_vector[neighbour_cell, 2]) * sf[face, 2] + ) + + +@qd.kernel +def gpu_rans_pressure_assembly( + n_cells: int, + p_internal: qd.types.NDArray[qd.f64, 1], + diag: qd.types.NDArray[qd.f64, 1], + source: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + diag[cell] = 0.0 + source[cell] = 0.0 + + +@qd.kernel +def gpu_rans_pressure_laplacian_coefficients( + n_internal_faces: int, + owner: qd.types.NDArray[qd.i32, 1], + neighbour: qd.types.NDArray[qd.i32, 1], + rAtU: qd.types.NDArray[qd.f64, 1], + cell_centres: qd.types.NDArray[qd.f64, 2], + sf: qd.types.NDArray[qd.f64, 2], + mag_sf: qd.types.NDArray[qd.f64, 1], + diag: qd.types.NDArray[qd.f64, 1], + upper: qd.types.NDArray[qd.f64, 1], +) -> None: + for face in range(n_internal_faces): + owner_cell = owner[face] + neighbour_cell = neighbour[face] + dx0 = cell_centres[neighbour_cell, 0] - cell_centres[owner_cell, 0] + dx1 = cell_centres[neighbour_cell, 1] - cell_centres[owner_cell, 1] + dx2 = cell_centres[neighbour_cell, 2] - cell_centres[owner_cell, 2] + projected_delta = dx0 * sf[face, 0] + dx1 * sf[face, 1] + dx2 * sf[face, 2] + if projected_delta < 0.0: + projected_delta = -projected_delta + if projected_delta < 1.0e-300: + projected_delta = 1.0e-300 + coeff = 0.5 * (rAtU[owner_cell] + rAtU[neighbour_cell]) * mag_sf[face] * mag_sf[face] / projected_delta + upper[face] = coeff + qd.atomic_add(diag[owner_cell], -coeff) + qd.atomic_add(diag[neighbour_cell], -coeff) + + +@qd.kernel +def gpu_rans_pressure_source_from_flux( + n_internal_faces: int, + owner: qd.types.NDArray[qd.i32, 1], + neighbour: qd.types.NDArray[qd.i32, 1], + phi: qd.types.NDArray[qd.f64, 1], + source: qd.types.NDArray[qd.f64, 1], +) -> None: + for face in range(n_internal_faces): + flux = phi[face] + qd.atomic_add(source[owner[face]], flux) + qd.atomic_add(source[neighbour[face]], -flux) + + +@qd.kernel +def gpu_rans_pressure_source_from_boundary_flux( + n_boundary_faces: int, + face_cells: qd.types.NDArray[qd.i32, 1], + phi_boundary: qd.types.NDArray[qd.f64, 1], + source: qd.types.NDArray[qd.f64, 1], +) -> None: + for face in range(n_boundary_faces): + qd.atomic_add(source[face_cells[face]], phi_boundary[face]) + + +@qd.kernel +def gpu_rans_pressure_mixed_boundary_laplacian( + n_boundary_faces: int, + face_cells: qd.types.NDArray[qd.i32, 1], + boundary_scales: qd.types.NDArray[qd.f64, 1], + boundary_values: qd.types.NDArray[qd.f64, 1], + rAtU: qd.types.NDArray[qd.f64, 1], + diag: qd.types.NDArray[qd.f64, 1], + source: qd.types.NDArray[qd.f64, 1], +) -> None: + for face in range(n_boundary_faces): + cell = face_cells[face] + coeff = rAtU[cell] * boundary_scales[face] + qd.atomic_add(diag[cell], -coeff) + qd.atomic_add(source[cell], -coeff * boundary_values[face]) + + +@qd.kernel +def gpu_rans_face_flux_copy( + n_internal_faces: int, + source: qd.types.NDArray[qd.f64, 1], + out: qd.types.NDArray[qd.f64, 1], +) -> None: + for face in range(n_internal_faces): + out[face] = source[face] + + +@qd.kernel +def gpu_rans_pressure_flux_correction( + n_internal_faces: int, + owner: qd.types.NDArray[qd.i32, 1], + neighbour: qd.types.NDArray[qd.i32, 1], + p_solved: qd.types.NDArray[qd.f64, 1], + upper: qd.types.NDArray[qd.f64, 1], + phiHbyA: qd.types.NDArray[qd.f64, 1], + phi_out: qd.types.NDArray[qd.f64, 1], +) -> None: + for face in range(n_internal_faces): + pressure_jump = p_solved[neighbour[face]] - p_solved[owner[face]] + phi_out[face] = phiHbyA[face] - upper[face] * pressure_jump + + +@qd.kernel +def gpu_rans_final_correction( + n_cells: int, + HbyA: qd.types.NDArray[qd.f64, 2], + p_solved: qd.types.NDArray[qd.f64, 1], + p_internal: qd.types.NDArray[qd.f64, 1], + u_out: qd.types.NDArray[qd.f64, 2], + p_out: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + u_out[cell, 0] = HbyA[cell, 0] + u_out[cell, 1] = HbyA[cell, 1] + u_out[cell, 2] = HbyA[cell, 2] + p_out[cell] = p_internal[cell] + p_solved[cell] + + +@qd.kernel +def gpu_rans_pressure_velocity_correction( + n_internal_faces: int, + owner: qd.types.NDArray[qd.i32, 1], + neighbour: qd.types.NDArray[qd.i32, 1], + rAU: qd.types.NDArray[qd.f64, 1], + p_solved: qd.types.NDArray[qd.f64, 1], + cell_centres: qd.types.NDArray[qd.f64, 2], + sf: qd.types.NDArray[qd.f64, 2], + u_out: qd.types.NDArray[qd.f64, 2], +) -> None: + for face in range(n_internal_faces): + owner_cell = owner[face] + neighbour_cell = neighbour[face] + dx0 = cell_centres[neighbour_cell, 0] - cell_centres[owner_cell, 0] + dx1 = cell_centres[neighbour_cell, 1] - cell_centres[owner_cell, 1] + dx2 = cell_centres[neighbour_cell, 2] - cell_centres[owner_cell, 2] + projected_delta = dx0 * sf[face, 0] + dx1 * sf[face, 1] + dx2 * sf[face, 2] + if projected_delta < 0.0: + projected_delta = -projected_delta + if projected_delta < 1.0e-300: + projected_delta = 1.0e-300 + gradient_scale = (p_solved[neighbour_cell] - p_solved[owner_cell]) / projected_delta + owner_scale = rAU[owner_cell] * gradient_scale + neighbour_scale = rAU[neighbour_cell] * gradient_scale + u_out[owner_cell, 0] -= owner_scale * sf[face, 0] + u_out[owner_cell, 1] -= owner_scale * sf[face, 1] + u_out[owner_cell, 2] -= owner_scale * sf[face, 2] + u_out[neighbour_cell, 0] -= neighbour_scale * sf[face, 0] + u_out[neighbour_cell, 1] -= neighbour_scale * sf[face, 1] + u_out[neighbour_cell, 2] -= neighbour_scale * sf[face, 2] + + +@qd.kernel +def gpu_rans_turbulence_update( + n_cells: int, + nut_internal: qd.types.NDArray[qd.f64, 1], + k_internal: qd.types.NDArray[qd.f64, 1], + omega_internal: qd.types.NDArray[qd.f64, 1], + nut_out: qd.types.NDArray[qd.f64, 1], + k_out: qd.types.NDArray[qd.f64, 1], + omega_out: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + nut_out[cell] = nut_internal[cell] + k_out[cell] = k_internal[cell] + omega_out[cell] = omega_internal[cell] + + +@qd.kernel +def gpu_rans_omega_wall_update( + n_wall_cells: int, + wall_cells: qd.types.NDArray[qd.i32, 1], + wall_distances: qd.types.NDArray[qd.f64, 1], + laminar_nu: float, + beta1: float, + omega_out: qd.types.NDArray[qd.f64, 1], +) -> None: + for index in range(n_wall_cells): + distance = wall_distances[index] + if distance < 1.0e-300: + distance = 1.0e-300 + omega_wall = 6.0 * laminar_nu / (beta1 * distance * distance) + cell = wall_cells[index] + if omega_wall > omega_out[cell]: + omega_out[cell] = omega_wall + +__all__ = [ + "GPU_STAGE_KERNEL_ENTRYPOINTS", + "gpu_rans_momentum_assembly", + "gpu_rans_momentum_diffusion_coefficients", + "gpu_rans_momentum_convection_coefficients", + "gpu_rans_momentum_equation_relaxation", + "gpu_rans_momentum_hbyA_source", + "gpu_rans_momentum_hbyA_face_accumulate", + "gpu_rans_momentum_hbyA_finish", + "gpu_rans_pressure_inputs", + "gpu_rans_consistent_rAtU", + "gpu_rans_pressure_assembly", + "gpu_rans_pressure_laplacian_coefficients", + "gpu_rans_pressure_source_from_flux", + "gpu_rans_pressure_source_from_boundary_flux", + "gpu_rans_pressure_mixed_boundary_laplacian", + "gpu_rans_face_flux_copy", + "gpu_rans_surface_flux_from_cells", + "gpu_rans_pressure_flux_correction", + "gpu_rans_pressure_velocity_correction", + "gpu_rans_final_correction", + "gpu_rans_turbulence_update", + "gpu_rans_omega_wall_update", +] diff --git a/python/src/foam_stepper/gpu/linear_solve.py b/python/src/foam_stepper/gpu/linear_solve.py new file mode 100644 index 0000000..1e7f664 --- /dev/null +++ b/python/src/foam_stepper/gpu/linear_solve.py @@ -0,0 +1,1310 @@ +"""GPU linear-solve primitives for OpenFOAM-style LDU matrices.""" + +from __future__ import annotations + +import dataclasses +from collections.abc import Mapping +from typing import Any + +import numpy as np +import quadrants as qd + + +@dataclasses.dataclass(frozen=True) +class GpuLduCsr: + """Device CSR view of OpenFOAM LDU off-diagonal coefficients.""" + + offsets: Any + columns: Any + coefficients: Any + n_entries: int + metadata: Mapping[str, Any] + + +@qd.kernel +def gpu_ldu_jacobi_scalar( + n_cells: int, + offsets: qd.types.NDArray[qd.i32, 1], + columns: qd.types.NDArray[qd.i32, 1], + coefficients: qd.types.NDArray[qd.f64, 1], + diag: qd.types.NDArray[qd.f64, 1], + source: qd.types.NDArray[qd.f64, 1], + current: qd.types.NDArray[qd.f64, 1], + out: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + residual = source[cell] + for coeff_index in range(offsets[cell], offsets[cell + 1]): + residual -= coefficients[coeff_index] * current[columns[coeff_index]] + out[cell] = residual / diag[cell] + + +@qd.kernel +def gpu_copy_scalar( + n_cells: int, + source: qd.types.NDArray[qd.f64, 1], + out: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + out[cell] = source[cell] + + +@qd.kernel +def gpu_ldu_jacobi_scalar_symmetric_accumulate( + n_internal_faces: int, + owner: qd.types.NDArray[qd.i32, 1], + neighbour: qd.types.NDArray[qd.i32, 1], + upper: qd.types.NDArray[qd.f64, 1], + current: qd.types.NDArray[qd.f64, 1], + residual: qd.types.NDArray[qd.f64, 1], +) -> None: + for face in range(n_internal_faces): + owner_cell = owner[face] + neighbour_cell = neighbour[face] + coeff = upper[face] + qd.atomic_add(residual[owner_cell], -coeff * current[neighbour_cell]) + qd.atomic_add(residual[neighbour_cell], -coeff * current[owner_cell]) + + +@qd.kernel +def gpu_scalar_jacobi_finish( + n_cells: int, + residual: qd.types.NDArray[qd.f64, 1], + diag: qd.types.NDArray[qd.f64, 1], + out: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + out[cell] = residual[cell] / diag[cell] + + + + +@qd.kernel +def gpu_zero_scalar_accumulator(accumulator: qd.types.NDArray[qd.f64, 1]) -> None: + accumulator[0] = 0.0 + + +@qd.kernel +def gpu_ldu_matvec_scalar_symmetric_diag( + n_cells: int, + diag: qd.types.NDArray[qd.f64, 1], + current: qd.types.NDArray[qd.f64, 1], + out: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + out[cell] = diag[cell] * current[cell] + + +@qd.kernel +def gpu_ldu_matvec_scalar_symmetric_face_accumulate( + n_internal_faces: int, + owner: qd.types.NDArray[qd.i32, 1], + neighbour: qd.types.NDArray[qd.i32, 1], + upper: qd.types.NDArray[qd.f64, 1], + current: qd.types.NDArray[qd.f64, 1], + out: qd.types.NDArray[qd.f64, 1], +) -> None: + for face in range(n_internal_faces): + owner_cell = owner[face] + neighbour_cell = neighbour[face] + coeff = upper[face] + qd.atomic_add(out[owner_cell], coeff * current[neighbour_cell]) + qd.atomic_add(out[neighbour_cell], coeff * current[owner_cell]) + + +def gpu_ldu_matvec_scalar_symmetric_faces( + n_cells: int, + n_internal_faces: int, + owner: Any, + neighbour: Any, + upper: Any, + diag: Any, + current: Any, + out: Any, +) -> None: + gpu_ldu_matvec_scalar_symmetric_diag(n_cells, diag, current, out) + gpu_ldu_matvec_scalar_symmetric_face_accumulate(n_internal_faces, owner, neighbour, upper, current, out) + + +@qd.kernel +def gpu_cg_initialize_scalar( + n_cells: int, + source: qd.types.NDArray[qd.f64, 1], + operator_current: qd.types.NDArray[qd.f64, 1], + initial: qd.types.NDArray[qd.f64, 1], + solution: qd.types.NDArray[qd.f64, 1], + residual: qd.types.NDArray[qd.f64, 1], + direction: qd.types.NDArray[qd.f64, 1], + residual_squared: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + solution[cell] = initial[cell] + cell_residual = source[cell] - operator_current[cell] + residual[cell] = cell_residual + direction[cell] = cell_residual + qd.atomic_add(residual_squared[0], cell_residual * cell_residual) + + +@qd.kernel +def gpu_cg_dot_scalar( + n_cells: int, + left: qd.types.NDArray[qd.f64, 1], + right: qd.types.NDArray[qd.f64, 1], + out: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + qd.atomic_add(out[0], left[cell] * right[cell]) + + +@qd.kernel +def gpu_cg_update_solution_residual_scalar( + n_cells: int, + alpha: float, + solution: qd.types.NDArray[qd.f64, 1], + direction: qd.types.NDArray[qd.f64, 1], + residual: qd.types.NDArray[qd.f64, 1], + operator_direction: qd.types.NDArray[qd.f64, 1], + residual_squared: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + next_solution = solution[cell] + alpha * direction[cell] + next_residual = residual[cell] - alpha * operator_direction[cell] + solution[cell] = next_solution + residual[cell] = next_residual + qd.atomic_add(residual_squared[0], next_residual * next_residual) + + +@qd.kernel +def gpu_cg_update_direction_scalar( + n_cells: int, + beta: float, + residual: qd.types.NDArray[qd.f64, 1], + direction: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + direction[cell] = residual[cell] + beta * direction[cell] + + +@qd.kernel +def gpu_pcg_initialize_scalar( + n_cells: int, + source: qd.types.NDArray[qd.f64, 1], + operator_current: qd.types.NDArray[qd.f64, 1], + diag: qd.types.NDArray[qd.f64, 1], + initial: qd.types.NDArray[qd.f64, 1], + solution: qd.types.NDArray[qd.f64, 1], + residual: qd.types.NDArray[qd.f64, 1], + preconditioned_residual: qd.types.NDArray[qd.f64, 1], + direction: qd.types.NDArray[qd.f64, 1], + residual_squared: qd.types.NDArray[qd.f64, 1], + preconditioned_dot: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + solution[cell] = initial[cell] + cell_residual = source[cell] - operator_current[cell] + diagonal = diag[cell] + diagonal_magnitude = diagonal + if diagonal_magnitude < 0.0: + diagonal_magnitude = -diagonal_magnitude + z = cell_residual + if diagonal_magnitude >= 1.0e-300: + z = cell_residual / diagonal + residual[cell] = cell_residual + preconditioned_residual[cell] = z + direction[cell] = z + qd.atomic_add(residual_squared[0], cell_residual * cell_residual) + qd.atomic_add(preconditioned_dot[0], cell_residual * z) + + +@qd.kernel +def gpu_pcg_update_solution_residual_scalar( + n_cells: int, + alpha: float, + diag: qd.types.NDArray[qd.f64, 1], + solution: qd.types.NDArray[qd.f64, 1], + direction: qd.types.NDArray[qd.f64, 1], + residual: qd.types.NDArray[qd.f64, 1], + operator_direction: qd.types.NDArray[qd.f64, 1], + preconditioned_residual: qd.types.NDArray[qd.f64, 1], + residual_squared: qd.types.NDArray[qd.f64, 1], + preconditioned_dot: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + next_solution = solution[cell] + alpha * direction[cell] + next_residual = residual[cell] - alpha * operator_direction[cell] + diagonal = diag[cell] + diagonal_magnitude = diagonal + if diagonal_magnitude < 0.0: + diagonal_magnitude = -diagonal_magnitude + z = next_residual + if diagonal_magnitude >= 1.0e-300: + z = next_residual / diagonal + solution[cell] = next_solution + residual[cell] = next_residual + preconditioned_residual[cell] = z + qd.atomic_add(residual_squared[0], next_residual * next_residual) + qd.atomic_add(preconditioned_dot[0], next_residual * z) + + +@qd.kernel +def gpu_pcg_update_direction_scalar( + n_cells: int, + beta: float, + preconditioned_residual: qd.types.NDArray[qd.f64, 1], + direction: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + direction[cell] = preconditioned_residual[cell] + beta * direction[cell] + + +def gpu_ldu_cg_scalar_symmetric_faces( + n_cells: int, + n_internal_faces: int, + owner: Any, + neighbour: Any, + upper: Any, + diag: Any, + source: Any, + initial: Any, + out: Any, + residual: Any, + direction: Any, + operator_work: Any, + residual_squared: Any, + denominator: Any, + *, + iterations: int, + residual_tolerance_squared: float = 0.0, +) -> dict[str, Any]: + """Run GPU conjugate-gradient iterations for a symmetric per-face LDU matrix.""" + + gpu_ldu_matvec_scalar_symmetric_faces(n_cells, n_internal_faces, owner, neighbour, upper, diag, initial, operator_work) + gpu_zero_scalar_accumulator(residual_squared) + gpu_cg_initialize_scalar(n_cells, source, operator_work, initial, out, residual, direction, residual_squared) + qd.sync() + rr_value = float(np.asarray(residual_squared.to_numpy())[0]) + performed_iterations = 0 + + for _ in range(iterations): + if rr_value <= residual_tolerance_squared: + break + gpu_ldu_matvec_scalar_symmetric_faces(n_cells, n_internal_faces, owner, neighbour, upper, diag, direction, operator_work) + gpu_zero_scalar_accumulator(denominator) + gpu_cg_dot_scalar(n_cells, direction, operator_work, denominator) + qd.sync() + denominator_value = float(np.asarray(denominator.to_numpy())[0]) + if not np.isfinite(denominator_value) or abs(denominator_value) <= 1.0e-300: + break + alpha = rr_value / denominator_value + gpu_zero_scalar_accumulator(residual_squared) + gpu_cg_update_solution_residual_scalar(n_cells, alpha, out, direction, residual, operator_work, residual_squared) + qd.sync() + next_rr_value = float(np.asarray(residual_squared.to_numpy())[0]) + performed_iterations += 1 + if next_rr_value <= residual_tolerance_squared: + rr_value = next_rr_value + break + beta = next_rr_value / rr_value if rr_value != 0.0 else 0.0 + gpu_cg_update_direction_scalar(n_cells, beta, residual, direction) + rr_value = next_rr_value + + return { + "iterations": performed_iterations, + "final_residual_squared": rr_value, + "converged": rr_value <= residual_tolerance_squared, + "residual_tolerance_squared": residual_tolerance_squared, + } + + +def gpu_ldu_pcg_scalar_symmetric_faces( + n_cells: int, + n_internal_faces: int, + owner: Any, + neighbour: Any, + upper: Any, + diag: Any, + source: Any, + initial: Any, + out: Any, + residual: Any, + preconditioned_residual: Any, + direction: Any, + operator_work: Any, + residual_squared: Any, + denominator: Any, + *, + iterations: int, + residual_tolerance_squared: float = 0.0, +) -> dict[str, Any]: + """Run diagonal-preconditioned GPU CG for a symmetric per-face LDU matrix.""" + + gpu_ldu_matvec_scalar_symmetric_faces(n_cells, n_internal_faces, owner, neighbour, upper, diag, initial, operator_work) + gpu_zero_scalar_accumulator(residual_squared) + gpu_zero_scalar_accumulator(denominator) + gpu_pcg_initialize_scalar(n_cells, source, operator_work, diag, initial, out, residual, preconditioned_residual, direction, residual_squared, denominator) + qd.sync() + rr_value = float(np.asarray(residual_squared.to_numpy())[0]) + rho_value = float(np.asarray(denominator.to_numpy())[0]) + performed_iterations = 0 + + for _ in range(iterations): + if rr_value <= residual_tolerance_squared: + break + if not np.isfinite(rho_value) or abs(rho_value) <= 1.0e-300: + break + gpu_ldu_matvec_scalar_symmetric_faces(n_cells, n_internal_faces, owner, neighbour, upper, diag, direction, operator_work) + gpu_zero_scalar_accumulator(denominator) + gpu_cg_dot_scalar(n_cells, direction, operator_work, denominator) + qd.sync() + denominator_value = float(np.asarray(denominator.to_numpy())[0]) + if not np.isfinite(denominator_value) or abs(denominator_value) <= 1.0e-300: + break + alpha = rho_value / denominator_value + gpu_zero_scalar_accumulator(residual_squared) + gpu_zero_scalar_accumulator(denominator) + gpu_pcg_update_solution_residual_scalar(n_cells, alpha, diag, out, direction, residual, operator_work, preconditioned_residual, residual_squared, denominator) + qd.sync() + next_rr_value = float(np.asarray(residual_squared.to_numpy())[0]) + next_rho_value = float(np.asarray(denominator.to_numpy())[0]) + performed_iterations += 1 + if next_rr_value <= residual_tolerance_squared: + rr_value = next_rr_value + rho_value = next_rho_value + break + beta = next_rho_value / rho_value if rho_value != 0.0 else 0.0 + gpu_pcg_update_direction_scalar(n_cells, beta, preconditioned_residual, direction) + rr_value = next_rr_value + rho_value = next_rho_value + + return { + "iterations": performed_iterations, + "final_residual_squared": rr_value, + "final_preconditioned_dot": rho_value, + "converged": rr_value <= residual_tolerance_squared, + "residual_tolerance_squared": residual_tolerance_squared, + "preconditioner": "diagonal_jacobi", + } +def gpu_ldu_jacobi_scalar_symmetric_faces( + n_cells: int, + n_internal_faces: int, + owner: Any, + neighbour: Any, + upper: Any, + diag: Any, + source: Any, + initial: Any, + out: Any, + scratch: Any, + work: Any, + *, + iterations: int, +) -> None: + """Run Jacobi iterations for symmetric OpenFOAM LDU coefficients stored per face.""" + + current = initial + next_field = out + for _ in range(iterations): + gpu_copy_scalar(n_cells, source, scratch) + gpu_ldu_jacobi_scalar_symmetric_accumulate(n_internal_faces, owner, neighbour, upper, current, scratch) + gpu_scalar_jacobi_finish(n_cells, scratch, diag, next_field) + current, next_field = next_field, work if next_field is out else out + if current is not out: + gpu_copy_scalar(n_cells, current, out) + +@qd.kernel +def gpu_ldu_jacobi_vector( + n_cells: int, + offsets: qd.types.NDArray[qd.i32, 1], + columns: qd.types.NDArray[qd.i32, 1], + coefficients: qd.types.NDArray[qd.f64, 1], + diag: qd.types.NDArray[qd.f64, 1], + source: qd.types.NDArray[qd.f64, 2], + current: qd.types.NDArray[qd.f64, 2], + out: qd.types.NDArray[qd.f64, 2], +) -> None: + for cell in range(n_cells): + residual_x = source[cell, 0] + residual_y = source[cell, 1] + residual_z = source[cell, 2] + for coeff_index in range(offsets[cell], offsets[cell + 1]): + neighbour = columns[coeff_index] + coeff = coefficients[coeff_index] + residual_x -= coeff * current[neighbour, 0] + residual_y -= coeff * current[neighbour, 1] + residual_z -= coeff * current[neighbour, 2] + inv_diag = 1.0 / diag[cell] + out[cell, 0] = residual_x * inv_diag + out[cell, 1] = residual_y * inv_diag + out[cell, 2] = residual_z * inv_diag + + +@qd.kernel +def gpu_copy_vector( + n_cells: int, + source: qd.types.NDArray[qd.f64, 2], + out: qd.types.NDArray[qd.f64, 2], +) -> None: + for cell in range(n_cells): + out[cell, 0] = source[cell, 0] + out[cell, 1] = source[cell, 1] + out[cell, 2] = source[cell, 2] + + + +@qd.kernel +def gpu_ldu_matvec_vector_asymmetric_diag( + n_cells: int, + diag: qd.types.NDArray[qd.f64, 1], + current: qd.types.NDArray[qd.f64, 2], + out: qd.types.NDArray[qd.f64, 2], +) -> None: + for cell in range(n_cells): + diagonal = diag[cell] + out[cell, 0] = diagonal * current[cell, 0] + out[cell, 1] = diagonal * current[cell, 1] + out[cell, 2] = diagonal * current[cell, 2] + + +@qd.kernel +def gpu_ldu_matvec_vector_asymmetric_face_accumulate( + n_internal_faces: int, + owner: qd.types.NDArray[qd.i32, 1], + neighbour: qd.types.NDArray[qd.i32, 1], + upper: qd.types.NDArray[qd.f64, 1], + lower: qd.types.NDArray[qd.f64, 1], + current: qd.types.NDArray[qd.f64, 2], + out: qd.types.NDArray[qd.f64, 2], +) -> None: + for face in range(n_internal_faces): + owner_cell = owner[face] + neighbour_cell = neighbour[face] + upper_coeff = upper[face] + lower_coeff = lower[face] + qd.atomic_add(out[owner_cell, 0], upper_coeff * current[neighbour_cell, 0]) + qd.atomic_add(out[owner_cell, 1], upper_coeff * current[neighbour_cell, 1]) + qd.atomic_add(out[owner_cell, 2], upper_coeff * current[neighbour_cell, 2]) + qd.atomic_add(out[neighbour_cell, 0], lower_coeff * current[owner_cell, 0]) + qd.atomic_add(out[neighbour_cell, 1], lower_coeff * current[owner_cell, 1]) + qd.atomic_add(out[neighbour_cell, 2], lower_coeff * current[owner_cell, 2]) + + +def gpu_ldu_matvec_vector_asymmetric_faces( + n_cells: int, + n_internal_faces: int, + owner: Any, + neighbour: Any, + upper: Any, + lower: Any, + diag: Any, + current: Any, + out: Any, +) -> None: + gpu_ldu_matvec_vector_asymmetric_diag(n_cells, diag, current, out) + gpu_ldu_matvec_vector_asymmetric_face_accumulate(n_internal_faces, owner, neighbour, upper, lower, current, out) + + +@qd.kernel +def gpu_bicgstab_initialize_vector( + n_cells: int, + source: qd.types.NDArray[qd.f64, 2], + operator_initial: qd.types.NDArray[qd.f64, 2], + initial: qd.types.NDArray[qd.f64, 2], + solution: qd.types.NDArray[qd.f64, 2], + residual: qd.types.NDArray[qd.f64, 2], + shadow: qd.types.NDArray[qd.f64, 2], + direction: qd.types.NDArray[qd.f64, 2], + operator_direction: qd.types.NDArray[qd.f64, 2], + residual_squared: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + solution[cell, 0] = initial[cell, 0] + solution[cell, 1] = initial[cell, 1] + solution[cell, 2] = initial[cell, 2] + residual_x = source[cell, 0] - operator_initial[cell, 0] + residual_y = source[cell, 1] - operator_initial[cell, 1] + residual_z = source[cell, 2] - operator_initial[cell, 2] + residual[cell, 0] = residual_x + residual[cell, 1] = residual_y + residual[cell, 2] = residual_z + shadow[cell, 0] = residual_x + shadow[cell, 1] = residual_y + shadow[cell, 2] = residual_z + direction[cell, 0] = 0.0 + direction[cell, 1] = 0.0 + direction[cell, 2] = 0.0 + operator_direction[cell, 0] = 0.0 + operator_direction[cell, 1] = 0.0 + operator_direction[cell, 2] = 0.0 + qd.atomic_add(residual_squared[0], residual_x * residual_x + residual_y * residual_y + residual_z * residual_z) + + +@qd.kernel +def gpu_bicgstab_dot_vector( + n_cells: int, + left: qd.types.NDArray[qd.f64, 2], + right: qd.types.NDArray[qd.f64, 2], + out: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + qd.atomic_add( + out[0], + left[cell, 0] * right[cell, 0] + left[cell, 1] * right[cell, 1] + left[cell, 2] * right[cell, 2], + ) + + +@qd.kernel +def gpu_bicgstab_update_direction_vector( + n_cells: int, + beta: float, + omega: float, + residual: qd.types.NDArray[qd.f64, 2], + direction: qd.types.NDArray[qd.f64, 2], + operator_direction: qd.types.NDArray[qd.f64, 2], +) -> None: + for cell in range(n_cells): + direction[cell, 0] = residual[cell, 0] + beta * (direction[cell, 0] - omega * operator_direction[cell, 0]) + direction[cell, 1] = residual[cell, 1] + beta * (direction[cell, 1] - omega * operator_direction[cell, 1]) + direction[cell, 2] = residual[cell, 2] + beta * (direction[cell, 2] - omega * operator_direction[cell, 2]) + + +@qd.kernel +def gpu_bicgstab_update_intermediate_vector( + n_cells: int, + alpha: float, + solution: qd.types.NDArray[qd.f64, 2], + direction: qd.types.NDArray[qd.f64, 2], + residual: qd.types.NDArray[qd.f64, 2], + operator_direction: qd.types.NDArray[qd.f64, 2], + intermediate: qd.types.NDArray[qd.f64, 2], + residual_squared: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + next_solution_x = solution[cell, 0] + alpha * direction[cell, 0] + next_solution_y = solution[cell, 1] + alpha * direction[cell, 1] + next_solution_z = solution[cell, 2] + alpha * direction[cell, 2] + intermediate_x = residual[cell, 0] - alpha * operator_direction[cell, 0] + intermediate_y = residual[cell, 1] - alpha * operator_direction[cell, 1] + intermediate_z = residual[cell, 2] - alpha * operator_direction[cell, 2] + solution[cell, 0] = next_solution_x + solution[cell, 1] = next_solution_y + solution[cell, 2] = next_solution_z + intermediate[cell, 0] = intermediate_x + intermediate[cell, 1] = intermediate_y + intermediate[cell, 2] = intermediate_z + qd.atomic_add(residual_squared[0], intermediate_x * intermediate_x + intermediate_y * intermediate_y + intermediate_z * intermediate_z) + + +@qd.kernel +def gpu_bicgstab_update_solution_residual_vector( + n_cells: int, + omega: float, + solution: qd.types.NDArray[qd.f64, 2], + intermediate: qd.types.NDArray[qd.f64, 2], + operator_intermediate: qd.types.NDArray[qd.f64, 2], + residual: qd.types.NDArray[qd.f64, 2], + residual_squared: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + next_solution_x = solution[cell, 0] + omega * intermediate[cell, 0] + next_solution_y = solution[cell, 1] + omega * intermediate[cell, 1] + next_solution_z = solution[cell, 2] + omega * intermediate[cell, 2] + residual_x = intermediate[cell, 0] - omega * operator_intermediate[cell, 0] + residual_y = intermediate[cell, 1] - omega * operator_intermediate[cell, 1] + residual_z = intermediate[cell, 2] - omega * operator_intermediate[cell, 2] + solution[cell, 0] = next_solution_x + solution[cell, 1] = next_solution_y + solution[cell, 2] = next_solution_z + residual[cell, 0] = residual_x + residual[cell, 1] = residual_y + residual[cell, 2] = residual_z + qd.atomic_add(residual_squared[0], residual_x * residual_x + residual_y * residual_y + residual_z * residual_z) + + +@qd.kernel +def gpu_bicgstab_precondition_vector( + n_cells: int, + diag: qd.types.NDArray[qd.f64, 1], + source: qd.types.NDArray[qd.f64, 2], + out: qd.types.NDArray[qd.f64, 2], +) -> None: + for cell in range(n_cells): + diagonal = diag[cell] + diagonal_magnitude = diagonal + if diagonal_magnitude < 0.0: + diagonal_magnitude = -diagonal_magnitude + if diagonal_magnitude >= 1.0e-300: + out[cell, 0] = source[cell, 0] / diagonal + out[cell, 1] = source[cell, 1] / diagonal + out[cell, 2] = source[cell, 2] / diagonal + else: + out[cell, 0] = source[cell, 0] + out[cell, 1] = source[cell, 1] + out[cell, 2] = source[cell, 2] + + +@qd.kernel +def gpu_bicgstab_update_intermediate_vector_preconditioned( + n_cells: int, + alpha: float, + solution: qd.types.NDArray[qd.f64, 2], + preconditioned_direction: qd.types.NDArray[qd.f64, 2], + residual: qd.types.NDArray[qd.f64, 2], + operator_direction: qd.types.NDArray[qd.f64, 2], + intermediate: qd.types.NDArray[qd.f64, 2], + residual_squared: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + step_x = preconditioned_direction[cell, 0] + step_y = preconditioned_direction[cell, 1] + step_z = preconditioned_direction[cell, 2] + intermediate_x = residual[cell, 0] - alpha * operator_direction[cell, 0] + intermediate_y = residual[cell, 1] - alpha * operator_direction[cell, 1] + intermediate_z = residual[cell, 2] - alpha * operator_direction[cell, 2] + solution[cell, 0] += alpha * step_x + solution[cell, 1] += alpha * step_y + solution[cell, 2] += alpha * step_z + intermediate[cell, 0] = intermediate_x + intermediate[cell, 1] = intermediate_y + intermediate[cell, 2] = intermediate_z + qd.atomic_add(residual_squared[0], intermediate_x * intermediate_x + intermediate_y * intermediate_y + intermediate_z * intermediate_z) + + +@qd.kernel +def gpu_bicgstab_update_solution_residual_vector_preconditioned( + n_cells: int, + omega: float, + solution: qd.types.NDArray[qd.f64, 2], + preconditioned_intermediate: qd.types.NDArray[qd.f64, 2], + intermediate: qd.types.NDArray[qd.f64, 2], + operator_intermediate: qd.types.NDArray[qd.f64, 2], + residual: qd.types.NDArray[qd.f64, 2], + residual_squared: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + step_x = preconditioned_intermediate[cell, 0] + step_y = preconditioned_intermediate[cell, 1] + step_z = preconditioned_intermediate[cell, 2] + residual_x = intermediate[cell, 0] - omega * operator_intermediate[cell, 0] + residual_y = intermediate[cell, 1] - omega * operator_intermediate[cell, 1] + residual_z = intermediate[cell, 2] - omega * operator_intermediate[cell, 2] + solution[cell, 0] += omega * step_x + solution[cell, 1] += omega * step_y + solution[cell, 2] += omega * step_z + residual[cell, 0] = residual_x + residual[cell, 1] = residual_y + residual[cell, 2] = residual_z + qd.atomic_add(residual_squared[0], residual_x * residual_x + residual_y * residual_y + residual_z * residual_z) + +@qd.kernel +def gpu_ldu_jacobi_vector_symmetric_accumulate( + n_internal_faces: int, + owner: qd.types.NDArray[qd.i32, 1], + neighbour: qd.types.NDArray[qd.i32, 1], + upper: qd.types.NDArray[qd.f64, 1], + current: qd.types.NDArray[qd.f64, 2], + residual: qd.types.NDArray[qd.f64, 2], +) -> None: + for face in range(n_internal_faces): + owner_cell = owner[face] + neighbour_cell = neighbour[face] + coeff = upper[face] + qd.atomic_add(residual[owner_cell, 0], -coeff * current[neighbour_cell, 0]) + qd.atomic_add(residual[owner_cell, 1], -coeff * current[neighbour_cell, 1]) + qd.atomic_add(residual[owner_cell, 2], -coeff * current[neighbour_cell, 2]) + qd.atomic_add(residual[neighbour_cell, 0], -coeff * current[owner_cell, 0]) + qd.atomic_add(residual[neighbour_cell, 1], -coeff * current[owner_cell, 1]) + qd.atomic_add(residual[neighbour_cell, 2], -coeff * current[owner_cell, 2]) + + +@qd.kernel +def gpu_ldu_jacobi_vector_asymmetric_accumulate( + n_internal_faces: int, + owner: qd.types.NDArray[qd.i32, 1], + neighbour: qd.types.NDArray[qd.i32, 1], + upper: qd.types.NDArray[qd.f64, 1], + lower: qd.types.NDArray[qd.f64, 1], + current: qd.types.NDArray[qd.f64, 2], + residual: qd.types.NDArray[qd.f64, 2], +) -> None: + for face in range(n_internal_faces): + owner_cell = owner[face] + neighbour_cell = neighbour[face] + upper_coeff = upper[face] + lower_coeff = lower[face] + qd.atomic_add(residual[owner_cell, 0], -upper_coeff * current[neighbour_cell, 0]) + qd.atomic_add(residual[owner_cell, 1], -upper_coeff * current[neighbour_cell, 1]) + qd.atomic_add(residual[owner_cell, 2], -upper_coeff * current[neighbour_cell, 2]) + qd.atomic_add(residual[neighbour_cell, 0], -lower_coeff * current[owner_cell, 0]) + qd.atomic_add(residual[neighbour_cell, 1], -lower_coeff * current[owner_cell, 1]) + qd.atomic_add(residual[neighbour_cell, 2], -lower_coeff * current[owner_cell, 2]) + +@qd.kernel +def gpu_vector_jacobi_finish( + n_cells: int, + residual: qd.types.NDArray[qd.f64, 2], + diag: qd.types.NDArray[qd.f64, 1], + out: qd.types.NDArray[qd.f64, 2], +) -> None: + for cell in range(n_cells): + inv_diag = 1.0 / diag[cell] + out[cell, 0] = residual[cell, 0] * inv_diag + out[cell, 1] = residual[cell, 1] * inv_diag + out[cell, 2] = residual[cell, 2] * inv_diag + + + +@qd.kernel +def gpu_vector_residual_squared( + n_cells: int, + source: qd.types.NDArray[qd.f64, 2], + operator_current: qd.types.NDArray[qd.f64, 2], + residual: qd.types.NDArray[qd.f64, 2], + residual_squared: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + residual_x = source[cell, 0] - operator_current[cell, 0] + residual_y = source[cell, 1] - operator_current[cell, 1] + residual_z = source[cell, 2] - operator_current[cell, 2] + residual[cell, 0] = residual_x + residual[cell, 1] = residual_y + residual[cell, 2] = residual_z + qd.atomic_add(residual_squared[0], residual_x * residual_x + residual_y * residual_y + residual_z * residual_z) + +def gpu_ldu_jacobi_vector_asymmetric_faces( + n_cells: int, + n_internal_faces: int, + owner: Any, + neighbour: Any, + upper: Any, + lower: Any, + diag: Any, + source: Any, + initial: Any, + out: Any, + scratch: Any, + work: Any, + *, + iterations: int, +) -> None: + """Run vector Jacobi iterations for asymmetric OpenFOAM LDU coefficients stored per face.""" + + current = initial + next_field = out + for _ in range(iterations): + gpu_copy_vector(n_cells, source, scratch) + gpu_ldu_jacobi_vector_asymmetric_accumulate(n_internal_faces, owner, neighbour, upper, lower, current, scratch) + gpu_vector_jacobi_finish(n_cells, scratch, diag, next_field) + current, next_field = next_field, work if next_field is out else out + if current is not out: + gpu_copy_vector(n_cells, current, out) + + + +def gpu_ldu_jacobi_vector_asymmetric_solve( + n_cells: int, + n_internal_faces: int, + owner: Any, + neighbour: Any, + upper: Any, + lower: Any, + diag: Any, + source: Any, + initial: Any, + out: Any, + scratch: Any, + work: Any, + operator_work: Any, + residual_squared: Any, + *, + iterations: int, + residual_tolerance_squared: float = 0.0, +) -> dict[str, Any]: + """Run a bounded Jacobi GPU solve for an asymmetric vector LDU matrix.""" + + gpu_ldu_jacobi_vector_asymmetric_faces( + n_cells, + n_internal_faces, + owner, + neighbour, + upper, + lower, + diag, + source, + initial, + out, + scratch, + work, + iterations=iterations, + ) + gpu_ldu_matvec_vector_asymmetric_faces(n_cells, n_internal_faces, owner, neighbour, upper, lower, diag, out, operator_work) + gpu_zero_scalar_accumulator(residual_squared) + gpu_vector_residual_squared(n_cells, source, operator_work, scratch, residual_squared) + qd.sync() + residual_value = float(np.asarray(residual_squared.to_numpy())[0]) + return { + "iterations": iterations, + "final_residual_squared": residual_value, + "converged": residual_value <= residual_tolerance_squared, + "residual_tolerance_squared": residual_tolerance_squared, + "preconditioner": "diagonal_jacobi", + } + + +def gpu_ldu_bicgstab_vector_asymmetric_faces( + n_cells: int, + n_internal_faces: int, + owner: Any, + neighbour: Any, + upper: Any, + lower: Any, + diag: Any, + source: Any, + initial: Any, + out: Any, + residual: Any, + shadow: Any, + direction: Any, + operator_direction: Any, + intermediate: Any, + operator_intermediate: Any, + residual_squared: Any, + rho: Any, + denominator: Any, + omega_numerator: Any, + omega_denominator: Any, + *, + iterations: int, + residual_tolerance_squared: float = 0.0, +) -> dict[str, Any]: + """Run unpreconditioned GPU BiCGStab iterations for an asymmetric vector LDU matrix.""" + + gpu_ldu_matvec_vector_asymmetric_faces( + n_cells, + n_internal_faces, + owner, + neighbour, + upper, + lower, + diag, + initial, + operator_intermediate, + ) + gpu_zero_scalar_accumulator(residual_squared) + gpu_bicgstab_initialize_vector( + n_cells, + source, + operator_intermediate, + initial, + out, + residual, + shadow, + direction, + operator_direction, + residual_squared, + ) + qd.sync() + residual_value = float(np.asarray(residual_squared.to_numpy())[0]) + rho_old = 1.0 + alpha = 1.0 + omega = 1.0 + performed_iterations = 0 + + for _ in range(iterations): + if residual_value <= residual_tolerance_squared: + break + + gpu_zero_scalar_accumulator(rho) + gpu_bicgstab_dot_vector(n_cells, shadow, residual, rho) + qd.sync() + rho_new = float(np.asarray(rho.to_numpy())[0]) + if not np.isfinite(rho_new) or abs(rho_new) <= 1.0e-300: + break + + if performed_iterations == 0: + beta = 0.0 + elif abs(omega) <= 1.0e-300: + break + else: + beta = (rho_new / rho_old) * (alpha / omega) + + gpu_bicgstab_update_direction_vector(n_cells, beta, omega, residual, direction, operator_direction) + gpu_ldu_matvec_vector_asymmetric_faces( + n_cells, + n_internal_faces, + owner, + neighbour, + upper, + lower, + diag, + direction, + operator_direction, + ) + + gpu_zero_scalar_accumulator(denominator) + gpu_bicgstab_dot_vector(n_cells, shadow, operator_direction, denominator) + qd.sync() + denominator_value = float(np.asarray(denominator.to_numpy())[0]) + if not np.isfinite(denominator_value) or abs(denominator_value) <= 1.0e-300: + break + alpha = rho_new / denominator_value + + gpu_zero_scalar_accumulator(residual_squared) + gpu_bicgstab_update_intermediate_vector( + n_cells, + alpha, + out, + direction, + residual, + operator_direction, + intermediate, + residual_squared, + ) + qd.sync() + intermediate_residual = float(np.asarray(residual_squared.to_numpy())[0]) + performed_iterations += 1 + if intermediate_residual <= residual_tolerance_squared: + residual_value = intermediate_residual + break + + gpu_ldu_matvec_vector_asymmetric_faces( + n_cells, + n_internal_faces, + owner, + neighbour, + upper, + lower, + diag, + intermediate, + operator_intermediate, + ) + gpu_zero_scalar_accumulator(omega_numerator) + gpu_zero_scalar_accumulator(omega_denominator) + gpu_bicgstab_dot_vector(n_cells, operator_intermediate, intermediate, omega_numerator) + gpu_bicgstab_dot_vector(n_cells, operator_intermediate, operator_intermediate, omega_denominator) + qd.sync() + omega_numerator_value = float(np.asarray(omega_numerator.to_numpy())[0]) + omega_denominator_value = float(np.asarray(omega_denominator.to_numpy())[0]) + if not np.isfinite(omega_numerator_value) or not np.isfinite(omega_denominator_value) or abs(omega_denominator_value) <= 1.0e-300: + break + omega = omega_numerator_value / omega_denominator_value + + gpu_zero_scalar_accumulator(residual_squared) + gpu_bicgstab_update_solution_residual_vector( + n_cells, + omega, + out, + intermediate, + operator_intermediate, + residual, + residual_squared, + ) + qd.sync() + residual_value = float(np.asarray(residual_squared.to_numpy())[0]) + rho_old = rho_new + + return { + "iterations": performed_iterations, + "final_residual_squared": residual_value, + "converged": residual_value <= residual_tolerance_squared, + "residual_tolerance_squared": residual_tolerance_squared, + } + + +def gpu_ldu_pbicgstab_vector_asymmetric_faces( + n_cells: int, + n_internal_faces: int, + owner: Any, + neighbour: Any, + upper: Any, + lower: Any, + diag: Any, + source: Any, + initial: Any, + out: Any, + residual: Any, + shadow: Any, + direction: Any, + operator_direction: Any, + intermediate: Any, + operator_intermediate: Any, + residual_squared: Any, + rho: Any, + denominator: Any, + omega_numerator: Any, + omega_denominator: Any, + *, + iterations: int, + residual_tolerance_squared: float = 0.0, +) -> dict[str, Any]: + """Run diagonal-preconditioned GPU BiCGStab for an asymmetric vector LDU matrix.""" + + gpu_ldu_matvec_vector_asymmetric_faces( + n_cells, + n_internal_faces, + owner, + neighbour, + upper, + lower, + diag, + initial, + operator_intermediate, + ) + gpu_zero_scalar_accumulator(residual_squared) + gpu_bicgstab_initialize_vector( + n_cells, + source, + operator_intermediate, + initial, + out, + residual, + shadow, + direction, + operator_direction, + residual_squared, + ) + qd.sync() + residual_value = float(np.asarray(residual_squared.to_numpy())[0]) + rho_old = 1.0 + alpha = 1.0 + omega = 1.0 + performed_iterations = 0 + + for _ in range(iterations): + if residual_value <= residual_tolerance_squared: + break + + gpu_zero_scalar_accumulator(rho) + gpu_bicgstab_dot_vector(n_cells, shadow, residual, rho) + qd.sync() + rho_new = float(np.asarray(rho.to_numpy())[0]) + if not np.isfinite(rho_new) or abs(rho_new) <= 1.0e-300: + break + + if performed_iterations == 0: + beta = 0.0 + elif abs(omega) <= 1.0e-300: + break + else: + beta = (rho_new / rho_old) * (alpha / omega) + + gpu_bicgstab_update_direction_vector(n_cells, beta, omega, residual, direction, operator_direction) + gpu_bicgstab_precondition_vector(n_cells, diag, direction, intermediate) + gpu_ldu_matvec_vector_asymmetric_faces(n_cells, n_internal_faces, owner, neighbour, upper, lower, diag, intermediate, operator_direction) + gpu_zero_scalar_accumulator(denominator) + gpu_bicgstab_dot_vector(n_cells, shadow, operator_direction, denominator) + qd.sync() + denominator_value = float(np.asarray(denominator.to_numpy())[0]) + if not np.isfinite(denominator_value) or abs(denominator_value) <= 1.0e-300: + break + alpha = rho_new / denominator_value + + gpu_zero_scalar_accumulator(residual_squared) + gpu_bicgstab_update_intermediate_vector_preconditioned( + n_cells, + alpha, + out, + intermediate, + residual, + operator_direction, + intermediate, + residual_squared, + ) + qd.sync() + intermediate_residual = float(np.asarray(residual_squared.to_numpy())[0]) + performed_iterations += 1 + if intermediate_residual <= residual_tolerance_squared: + residual_value = intermediate_residual + break + + gpu_bicgstab_precondition_vector(n_cells, diag, intermediate, residual) + gpu_ldu_matvec_vector_asymmetric_faces(n_cells, n_internal_faces, owner, neighbour, upper, lower, diag, residual, operator_intermediate) + gpu_zero_scalar_accumulator(omega_numerator) + gpu_zero_scalar_accumulator(omega_denominator) + gpu_bicgstab_dot_vector(n_cells, operator_intermediate, intermediate, omega_numerator) + gpu_bicgstab_dot_vector(n_cells, operator_intermediate, operator_intermediate, omega_denominator) + qd.sync() + omega_numerator_value = float(np.asarray(omega_numerator.to_numpy())[0]) + omega_denominator_value = float(np.asarray(omega_denominator.to_numpy())[0]) + if not np.isfinite(omega_numerator_value) or not np.isfinite(omega_denominator_value) or abs(omega_denominator_value) <= 1.0e-300: + break + omega = omega_numerator_value / omega_denominator_value + + gpu_zero_scalar_accumulator(residual_squared) + gpu_bicgstab_update_solution_residual_vector_preconditioned( + n_cells, + omega, + out, + residual, + intermediate, + operator_intermediate, + residual, + residual_squared, + ) + qd.sync() + residual_value = float(np.asarray(residual_squared.to_numpy())[0]) + rho_old = rho_new + + return { + "iterations": performed_iterations, + "final_residual_squared": residual_value, + "converged": residual_value <= residual_tolerance_squared, + "residual_tolerance_squared": residual_tolerance_squared, + "preconditioner": "diagonal_jacobi", + } + +def gpu_ldu_jacobi_vector_symmetric_faces( + n_cells: int, + n_internal_faces: int, + owner: Any, + neighbour: Any, + upper: Any, + diag: Any, + source: Any, + initial: Any, + out: Any, + scratch: Any, + work: Any, + *, + iterations: int, +) -> None: + """Run vector Jacobi iterations for symmetric OpenFOAM LDU coefficients stored per face.""" + + current = initial + next_field = out + for _ in range(iterations): + gpu_copy_vector(n_cells, source, scratch) + gpu_ldu_jacobi_vector_symmetric_accumulate(n_internal_faces, owner, neighbour, upper, current, scratch) + gpu_vector_jacobi_finish(n_cells, scratch, diag, next_field) + current, next_field = next_field, work if next_field is out else out + if current is not out: + gpu_copy_vector(n_cells, current, out) + +def build_ldu_csr( + n_cells: int, + owner: np.ndarray, + neighbour: np.ndarray, + upper: np.ndarray | None, + lower: np.ndarray | None, +) -> tuple[np.ndarray, np.ndarray, np.ndarray]: + """Build row-wise off-diagonal CSR arrays from OpenFOAM owner/neighbour LDU data.""" + + owner_i32 = np.asarray(owner, dtype=np.int32) + neighbour_i32 = np.asarray(neighbour, dtype=np.int32) + if owner_i32.shape != neighbour_i32.shape: + raise ValueError(f"owner/neighbour shape mismatch: {owner_i32.shape} != {neighbour_i32.shape}") + if owner_i32.size == 0: + return np.zeros(n_cells + 1, dtype=np.int32), np.zeros(0, dtype=np.int32), np.zeros(0, dtype=np.float64) + if int(owner_i32.min()) < 0 or int(neighbour_i32.min()) < 0 or int(max(owner_i32.max(), neighbour_i32.max())) >= n_cells: + raise ValueError("owner/neighbour addresses exceed cell range") + + upper_f64 = np.asarray(upper if upper is not None else np.zeros(owner_i32.shape, dtype=np.float64), dtype=np.float64) + lower_f64 = np.asarray(lower if lower is not None else upper_f64, dtype=np.float64) + if upper_f64.shape != owner_i32.shape or lower_f64.shape != owner_i32.shape: + raise ValueError("upper/lower coefficient shapes must match owner/neighbour") + + counts = np.zeros(n_cells, dtype=np.int32) + np.add.at(counts, owner_i32, 1) + np.add.at(counts, neighbour_i32, 1) + offsets = np.empty(n_cells + 1, dtype=np.int32) + offsets[0] = 0 + np.cumsum(counts, out=offsets[1:]) + columns = np.empty(int(offsets[-1]), dtype=np.int32) + coefficients = np.empty(int(offsets[-1]), dtype=np.float64) + cursor = offsets[:-1].copy() + + for face, owner_cell in enumerate(owner_i32): + neighbour_cell = int(neighbour_i32[face]) + owner_row = int(owner_cell) + owner_slot = int(cursor[owner_row]) + columns[owner_slot] = neighbour_cell + coefficients[owner_slot] = upper_f64[face] + cursor[owner_row] += 1 + + neighbour_slot = int(cursor[neighbour_cell]) + columns[neighbour_slot] = owner_row + coefficients[neighbour_slot] = lower_f64[face] + cursor[neighbour_cell] += 1 + + return offsets, columns, coefficients + + +def _gpu_i32(values: np.ndarray) -> Any: + host = np.ascontiguousarray(values.astype(np.int32, copy=False)) + alloc_shape = host.shape if host.size else (1,) + gpu = qd.ndarray(qd.i32, shape=alloc_shape) + if host.size: + gpu.from_numpy(host) + return gpu + + +def _gpu_f64(values: np.ndarray) -> Any: + host = np.ascontiguousarray(values.astype(np.float64, copy=False)) + alloc_shape = host.shape if host.size else (1,) + gpu = qd.ndarray(qd.f64, shape=alloc_shape) + if host.size: + gpu.from_numpy(host) + return gpu + + +def ldu_csr_to_gpu(offsets: np.ndarray, columns: np.ndarray, coefficients: np.ndarray, *, name: str) -> GpuLduCsr: + return GpuLduCsr( + offsets=_gpu_i32(offsets), + columns=_gpu_i32(columns), + coefficients=_gpu_f64(coefficients), + n_entries=int(columns.size), + metadata={ + "name": name, + "format": "ldu_csr_offdiag", + "offsets_shape": [int(dim) for dim in offsets.shape], + "columns_shape": [int(dim) for dim in columns.shape], + "coefficients_shape": [int(dim) for dim in coefficients.shape], + "n_entries": int(columns.size), + "gpu_backed": True, + }, + ) + + +def empty_ldu_csr_gpu(n_cells: int, name: str) -> GpuLduCsr: + offsets = np.zeros(n_cells + 1, dtype=np.int32) + columns = np.zeros(0, dtype=np.int32) + coefficients = np.zeros(0, dtype=np.float64) + return ldu_csr_to_gpu(offsets, columns, coefficients, name=name) + + +__all__ = [ + "GpuLduCsr", + "build_ldu_csr", + "empty_ldu_csr_gpu", + "gpu_bicgstab_dot_vector", + "gpu_bicgstab_initialize_vector", + "gpu_bicgstab_precondition_vector", + "gpu_bicgstab_update_direction_vector", + "gpu_bicgstab_update_intermediate_vector", + "gpu_bicgstab_update_intermediate_vector_preconditioned", + "gpu_bicgstab_update_solution_residual_vector", + "gpu_bicgstab_update_solution_residual_vector_preconditioned", + "gpu_copy_scalar", + "gpu_copy_vector", + "gpu_ldu_bicgstab_vector_asymmetric_faces", + "gpu_ldu_pbicgstab_vector_asymmetric_faces", + "gpu_ldu_cg_scalar_symmetric_faces", + "gpu_ldu_pcg_scalar_symmetric_faces", + "gpu_ldu_jacobi_scalar", + "gpu_ldu_jacobi_scalar_symmetric_accumulate", + "gpu_ldu_jacobi_scalar_symmetric_faces", + "gpu_pcg_initialize_scalar", + "gpu_pcg_update_direction_scalar", + "gpu_pcg_update_solution_residual_scalar", + "gpu_ldu_jacobi_vector", + "gpu_ldu_jacobi_vector_asymmetric_accumulate", + "gpu_ldu_jacobi_vector_asymmetric_faces", + "gpu_ldu_jacobi_vector_asymmetric_solve", + "gpu_ldu_jacobi_vector_symmetric_accumulate", + "gpu_ldu_jacobi_vector_symmetric_faces", + "gpu_ldu_matvec_scalar_symmetric_faces", + "gpu_ldu_matvec_vector_asymmetric_faces", + "gpu_scalar_jacobi_finish", + "gpu_vector_jacobi_finish", + "gpu_vector_residual_squared", + "ldu_csr_to_gpu", +] diff --git a/python/src/foam_stepper/state.py b/python/src/foam_stepper/state.py new file mode 100644 index 0000000..c434fe6 --- /dev/null +++ b/python/src/foam_stepper/state.py @@ -0,0 +1,453 @@ +"""Explicit solver-state exports and validation diagnostics.""" + +from __future__ import annotations + +from collections.abc import Mapping +from typing import Any, Iterable + +import numpy as np + + +DEFAULT_TURBULENCE_FIELDS = ("nut", "k", "omega") + + +class SolverStateValidationError(ValueError): + """Localized validation failure for exported solver-state data.""" + + def __init__(self, path: str, message: str, *, expected: Any = None, actual: Any = None) -> None: + super().__init__(f"{path}: {message}") + self.path = path + self.message = message + self.expected = expected + self.actual = actual + + def to_dict(self) -> dict[str, Any]: + out: dict[str, Any] = {"path": self.path, "message": self.message} + if self.expected is not None: + out["expected"] = _json_ready(self.expected) + if self.actual is not None: + out["actual"] = _json_ready(self.actual) + return out + + +def export_solver_state( + mesh: Any, + fields: Mapping[str, Any] | Any, + *, + required_fields: Iterable[str] = (), + turbulence_fields: Iterable[str] = DEFAULT_TURBULENCE_FIELDS, +) -> dict[str, Any]: + """Return mesh, field, boundary, and turbulence arrays as plain Python/NumPy data.""" + + field_map = _field_mapping(fields) + state = { + "schema_version": 1, + "mesh": _export_mesh(mesh), + "fields": {name: _export_field(field) for name, field in field_map.items()}, + "turbulence": {name: _export_field(field_map[name]) for name in turbulence_fields if name in field_map}, + } + validate_solver_state(state, required_fields=required_fields) + return state + + +def export_matrix_state(matrix: Any, mesh_state: Mapping[str, Any] | None = None) -> dict[str, Any]: + """Return an assembled matrix and solver intermediates as explicit arrays.""" + + state = { + "name": getattr(matrix, "name", ""), + "field_name": getattr(matrix, "field_name", ""), + "value_rank": getattr(matrix, "value_rank", ""), + "dimensions": getattr(matrix, "dimensions", ""), + "flags": { + "has_diag": bool(getattr(matrix, "has_diag", False)), + "has_upper": bool(getattr(matrix, "has_upper", False)), + "has_lower": bool(getattr(matrix, "has_lower", False)), + "diagonal": bool(getattr(matrix, "diagonal", False)), + "symmetric": bool(getattr(matrix, "symmetric", False)), + "asymmetric": bool(getattr(matrix, "asymmetric", False)), + }, + "diag": _array(getattr(matrix, "diag")), + "upper": _optional_array(getattr(matrix, "upper", None)), + "lower": _optional_array(getattr(matrix, "lower", None)), + "source": _array(getattr(matrix, "source")), + "psi": _export_field(getattr(matrix, "psi")), + "internal_coeffs": [_array(item) for item in getattr(matrix, "internal_coeffs", [])], + "boundary_coeffs": [_array(item) for item in getattr(matrix, "boundary_coeffs", [])], + "derived": {}, + } + for name in ("A", "H", "H1", "flux", "face_flux_correction"): + value = getattr(matrix, name, None) + if callable(value): + value = value() + if value is not None: + state["derived"][name] = _export_field(value) + for name in ("residual", "D", "DD"): + value = getattr(matrix, name, None) + if callable(value): + value = value() + if value is not None: + state["derived"][name] = _array(value) + validate_matrix_state(state, mesh_state=mesh_state) + return state + + +def validate_solver_state( + state: Mapping[str, Any], + *, + required_fields: Iterable[str] = (), + turbulence_fields: Iterable[str] = (), +) -> None: + mesh = _mapping(state.get("mesh"), "mesh") + sizes = _mapping(mesh.get("sizes"), "mesh.sizes") + n_points = _positive_int(sizes.get("n_points"), "mesh.sizes.n_points", allow_zero=False) + n_faces = _positive_int(sizes.get("n_faces"), "mesh.sizes.n_faces", allow_zero=False) + n_internal_faces = _positive_int(sizes.get("n_internal_faces"), "mesh.sizes.n_internal_faces", allow_zero=True) + n_cells = _positive_int(sizes.get("n_cells"), "mesh.sizes.n_cells", allow_zero=False) + if n_internal_faces > n_faces: + raise SolverStateValidationError( + "mesh.sizes.n_internal_faces", + "cannot exceed n_faces", + expected={"max": n_faces}, + actual=n_internal_faces, + ) + + connectivity = _mapping(mesh.get("connectivity"), "mesh.connectivity") + geometry = _mapping(mesh.get("geometry"), "mesh.geometry") + _require_array(geometry.get("points"), "mesh.geometry.points", shape=(n_points, 3), numeric=True) + _require_array(geometry.get("V"), "mesh.geometry.V", shape=(n_cells,), numeric=True) + _require_array(geometry.get("C"), "mesh.geometry.C", shape=(n_cells, 3), numeric=True) + _require_array(geometry.get("Cf"), "mesh.geometry.Cf", shape=(n_internal_faces, 3), numeric=True) + _require_array(geometry.get("Sf"), "mesh.geometry.Sf", shape=(n_internal_faces, 3), numeric=True) + _require_array(geometry.get("magSf"), "mesh.geometry.magSf", shape=(n_internal_faces,), numeric=True) + + _validate_ragged(connectivity.get("faces"), "mesh.connectivity.faces", rows=n_faces) + _validate_ragged(connectivity.get("cells"), "mesh.connectivity.cells", rows=n_cells) + _require_array(connectivity.get("owner"), "mesh.connectivity.owner", shape=(n_internal_faces,), integer=True) + _require_array(connectivity.get("neighbour"), "mesh.connectivity.neighbour", shape=(n_internal_faces,), integer=True) + + ldu = _mapping(connectivity.get("ldu"), "mesh.connectivity.ldu") + _require_array(ldu.get("lower_addr"), "mesh.connectivity.ldu.lower_addr", shape=(n_internal_faces,), integer=True) + _require_array(ldu.get("upper_addr"), "mesh.connectivity.ldu.upper_addr", shape=(n_internal_faces,), integer=True) + + patches = _patches(mesh.get("patches"), n_faces=n_faces, n_cells=n_cells) + fields = _mapping(state.get("fields"), "fields") + for name in required_fields: + if name not in fields: + raise SolverStateValidationError("fields", "missing required field", expected=name, actual=sorted(fields)) + for name in turbulence_fields: + if name not in fields: + raise SolverStateValidationError("turbulence", "missing turbulence field", expected=name, actual=sorted(fields)) + + for name, field in fields.items(): + _validate_field(_mapping(field, f"fields.{name}"), f"fields.{name}", patches=patches, n_cells=n_cells, n_internal_faces=n_internal_faces) + + +def validate_matrix_state(state: Mapping[str, Any], *, mesh_state: Mapping[str, Any] | None = None) -> None: + path = f"matrices.{state.get('field_name', '')}" + rank = _rank_from_value_rank(str(state.get("value_rank", "")), path) + n_cells = None + n_internal_faces = None + patches: list[Mapping[str, Any]] = [] + if mesh_state is not None: + mesh = _mapping(mesh_state.get("mesh", mesh_state), "mesh") + sizes = _mapping(mesh.get("sizes"), "mesh.sizes") + n_cells = _positive_int(sizes.get("n_cells"), "mesh.sizes.n_cells", allow_zero=False) + n_internal_faces = _positive_int(sizes.get("n_internal_faces"), "mesh.sizes.n_internal_faces", allow_zero=True) + patches = list(_mapping(mesh, "mesh").get("patches", [])) + + if n_cells is not None: + _require_array(state.get("diag"), f"{path}.diag", shape=(n_cells,), numeric=True) + _require_array(state.get("source"), f"{path}.source", shape=_ranked_shape(n_cells, rank), numeric=True) + else: + _require_array(state.get("diag"), f"{path}.diag", numeric=True) + _require_array(state.get("source"), f"{path}.source", numeric=True) + + if n_internal_faces is not None: + _optional_valid_array(state.get("upper"), f"{path}.upper", shape=(n_internal_faces,), numeric=True) + _optional_valid_array(state.get("lower"), f"{path}.lower", shape=(n_internal_faces,), numeric=True) + else: + _optional_valid_array(state.get("upper"), f"{path}.upper", numeric=True) + _optional_valid_array(state.get("lower"), f"{path}.lower", numeric=True) + + if patches: + _validate_coeff_list(state.get("internal_coeffs", []), f"{path}.internal_coeffs", patches=patches, rank=rank) + _validate_coeff_list(state.get("boundary_coeffs", []), f"{path}.boundary_coeffs", patches=patches, rank=rank) + + +def describe_solver_state(state: Mapping[str, Any]) -> dict[str, Any]: + return _describe(state) + + +def describe_matrix_state(state: Mapping[str, Any]) -> dict[str, Any]: + return _describe(state) + + +def _export_mesh(mesh: Any) -> dict[str, Any]: + return { + "sizes": { + "n_points": int(getattr(mesh, "n_points")), + "n_faces": int(getattr(mesh, "n_faces")), + "n_internal_faces": int(getattr(mesh, "n_internal_faces")), + "n_cells": int(getattr(mesh, "n_cells")), + }, + "connectivity": { + "faces": _export_ragged(getattr(mesh, "faces")), + "cells": _export_ragged(getattr(mesh, "cells")), + "owner": _array(getattr(mesh, "owner")), + "neighbour": _array(getattr(mesh, "neighbour")), + "ldu": _export_ldu(getattr(mesh, "ldu", {})), + }, + "geometry": { + "points": _array(getattr(mesh, "points")), + "V": _array(getattr(mesh, "V")), + "C": _array(getattr(mesh, "C")), + "Cf": _array(getattr(mesh, "Cf")), + "Sf": _array(getattr(mesh, "Sf")), + "magSf": _array(getattr(mesh, "magSf")), + }, + "patches": [_export_patch(patch) for patch in getattr(mesh, "boundary")], + } + + +def _export_ragged(value: Any) -> dict[str, np.ndarray]: + return {"offsets": _array(getattr(value, "offsets")), "values": _array(getattr(value, "values"))} + + +def _export_ldu(ldu: Mapping[str, Any]) -> dict[str, Any]: + return { + "lower_addr": _array(ldu.get("lower_addr")), + "upper_addr": _array(ldu.get("upper_addr")), + "patch_addr": [ + {"patch_index": int(entry.get("patch_index")), "addr": _array(entry.get("addr"))} + for entry in ldu.get("patch_addr", []) + ], + } + + +def _export_patch(patch: Any) -> dict[str, Any]: + return { + "name": getattr(patch, "name"), + "type": getattr(patch, "type"), + "index": int(getattr(patch, "index")), + "start": int(getattr(patch, "start")), + "size": int(getattr(patch, "size")), + "coupled": bool(getattr(patch, "coupled")), + "constraint": bool(getattr(patch, "constraint")), + "face_cells": _array(getattr(patch, "face_cells")), + "face_indices": _array(getattr(patch, "face_indices")), + "Cf": _array(getattr(patch, "Cf")), + "Sf": _array(getattr(patch, "Sf")), + "magSf": _array(getattr(patch, "magSf")), + } + + +def _export_field(field: Any) -> dict[str, Any]: + return { + "name": getattr(field, "name"), + "kind": getattr(field, "kind"), + "rank": _rank_from_kind(getattr(field, "kind", "")), + "dimensions": getattr(field, "dimensions"), + "entity_kind": getattr(field, "entity_kind"), + "entity_count": int(getattr(field, "entity_count")), + "internal": _array(getattr(field, "internal")), + "boundary": {name: _export_patch_field(patch) for name, patch in getattr(field, "boundary").items()}, + } + + +def _export_patch_field(patch: Any) -> dict[str, Any]: + return { + "name": getattr(patch, "name"), + "type": getattr(patch, "type"), + "values": _array(getattr(patch, "values")), + "fixes_value": bool(getattr(patch, "fixes_value")), + "assignable": bool(getattr(patch, "assignable")), + "coupled": bool(getattr(patch, "coupled")), + "updated": bool(getattr(patch, "updated")), + "patch_internal": _optional_array(getattr(patch, "patch_internal", None)), + "value_internal_coeffs": _optional_array(getattr(patch, "value_internal_coeffs", None)), + "value_boundary_coeffs": _optional_array(getattr(patch, "value_boundary_coeffs", None)), + "gradient_internal_coeffs": _optional_array(getattr(patch, "gradient_internal_coeffs", None)), + "gradient_boundary_coeffs": _optional_array(getattr(patch, "gradient_boundary_coeffs", None)), + } + + +def _field_mapping(fields: Mapping[str, Any] | Any) -> Mapping[str, Any]: + if isinstance(fields, Mapping): + return fields + if hasattr(fields, "fields") and isinstance(fields.fields, Mapping): + return fields.fields + raise SolverStateValidationError("fields", "expected field mapping", actual=type(fields).__name__) + + +def _array(value: Any) -> np.ndarray: + if value is None: + raise SolverStateValidationError("array", "missing required array") + return np.asarray(value) + + +def _optional_array(value: Any) -> np.ndarray | None: + return None if value is None else np.asarray(value) + + +def _optional_valid_array(value: Any, path: str, *, shape: tuple[int, ...] | None = None, numeric: bool = False) -> None: + if value is None: + return + _require_array(value, path, shape=shape, numeric=numeric) + + +def _require_array( + value: Any, + path: str, + *, + shape: tuple[int, ...] | None = None, + numeric: bool = False, + integer: bool = False, +) -> np.ndarray: + if value is None: + raise SolverStateValidationError(path, "missing required array") + arr = np.asarray(value) + if shape is not None and tuple(arr.shape) != shape: + raise SolverStateValidationError(path, "shape mismatch", expected=list(shape), actual=list(arr.shape)) + if numeric and not np.issubdtype(arr.dtype, np.number): + raise SolverStateValidationError(path, "expected numeric dtype", actual=str(arr.dtype)) + if integer and not np.issubdtype(arr.dtype, np.integer): + raise SolverStateValidationError(path, "expected integer dtype", actual=str(arr.dtype)) + return arr + + +def _validate_ragged(value: Any, path: str, *, rows: int) -> None: + mapping = _mapping(value, path) + offsets = _require_array(mapping.get("offsets"), f"{path}.offsets", shape=(rows + 1,), integer=True) + _require_array(mapping.get("values"), f"{path}.values", integer=True) + if offsets.size and int(offsets[0]) != 0: + raise SolverStateValidationError(f"{path}.offsets", "first offset must be zero", expected=0, actual=int(offsets[0])) + if offsets.size > 1 and bool(np.any(offsets[1:] < offsets[:-1])): + raise SolverStateValidationError(f"{path}.offsets", "offsets must be monotonically nondecreasing") + + +def _patches(value: Any, *, n_faces: int, n_cells: int) -> list[Mapping[str, Any]]: + if not isinstance(value, list): + raise SolverStateValidationError("mesh.patches", "expected patch list", actual=type(value).__name__) + patches = [] + names: set[str] = set() + for index, patch_value in enumerate(value): + path = f"mesh.patches[{index}]" + patch = _mapping(patch_value, path) + name = str(patch.get("name", "")) + if not name: + raise SolverStateValidationError(f"{path}.name", "missing patch name") + if name in names: + raise SolverStateValidationError(f"{path}.name", "duplicate patch name", actual=name) + names.add(name) + size = _positive_int(patch.get("size"), f"{path}.size", allow_zero=True) + start = _positive_int(patch.get("start"), f"{path}.start", allow_zero=True) + if start + size > n_faces: + raise SolverStateValidationError(f"{path}.start", "patch faces exceed mesh face count", expected={"n_faces": n_faces}, actual={"start": start, "size": size}) + _require_array(patch.get("face_cells"), f"{path}.face_cells", shape=(size,), integer=True) + _require_array(patch.get("face_indices"), f"{path}.face_indices", shape=(size,), integer=True) + _require_array(patch.get("Cf"), f"{path}.Cf", shape=(size, 3), numeric=True) + _require_array(patch.get("Sf"), f"{path}.Sf", shape=(size, 3), numeric=True) + _require_array(patch.get("magSf"), f"{path}.magSf", shape=(size,), numeric=True) + face_cells = np.asarray(patch.get("face_cells")) + if face_cells.size and (int(face_cells.min()) < 0 or int(face_cells.max()) >= n_cells): + raise SolverStateValidationError(f"{path}.face_cells", "face cell index out of range", expected={"min": 0, "max_exclusive": n_cells}) + patches.append(patch) + return patches + + +def _validate_field(field: Mapping[str, Any], path: str, *, patches: list[Mapping[str, Any]], n_cells: int, n_internal_faces: int) -> None: + rank = str(field.get("rank") or _rank_from_kind(str(field.get("kind", "")))) + entity_kind = str(field.get("entity_kind", "")) + entity_count = _positive_int(field.get("entity_count"), f"{path}.entity_count", allow_zero=True) + expected_count = n_internal_faces if entity_kind == "internal_face" else n_cells + if entity_count != expected_count: + raise SolverStateValidationError(f"{path}.entity_count", "topology entity count mismatch", expected=expected_count, actual=entity_count) + _require_array(field.get("internal"), f"{path}.internal", shape=_ranked_shape(expected_count, rank), numeric=True) + boundary = _mapping(field.get("boundary"), f"{path}.boundary") + patch_names = {str(patch["name"]) for patch in patches} + if set(boundary) != patch_names: + raise SolverStateValidationError(f"{path}.boundary", "patch set mismatch", expected=sorted(patch_names), actual=sorted(boundary)) + patches_by_name = {str(patch["name"]): patch for patch in patches} + for patch_name, patch_field_value in boundary.items(): + patch_field_path = f"{path}.boundary.{patch_name}" + patch_field = _mapping(patch_field_value, patch_field_path) + patch_size = int(patches_by_name[str(patch_name)]["size"]) + expected_shape = _ranked_shape(patch_size, rank) + _require_array(patch_field.get("values"), f"{patch_field_path}.values", shape=expected_shape, numeric=True) + _optional_valid_array(patch_field.get("patch_internal"), f"{patch_field_path}.patch_internal", shape=expected_shape, numeric=True) + for coeff_name in ("value_internal_coeffs", "value_boundary_coeffs", "gradient_internal_coeffs", "gradient_boundary_coeffs"): + _optional_valid_array(patch_field.get(coeff_name), f"{patch_field_path}.{coeff_name}", numeric=True) + + +def _validate_coeff_list(value: Any, path: str, *, patches: list[Mapping[str, Any]], rank: str) -> None: + if not isinstance(value, list): + raise SolverStateValidationError(path, "expected coefficient list", actual=type(value).__name__) + if len(value) != len(patches): + raise SolverStateValidationError(path, "patch coefficient count mismatch", expected=len(patches), actual=len(value)) + for index, coeffs in enumerate(value): + patch_size = int(patches[index]["size"]) + _require_array(coeffs, f"{path}[{index}]", shape=_ranked_shape(patch_size, rank), numeric=True) + + +def _mapping(value: Any, path: str) -> Mapping[str, Any]: + if not isinstance(value, Mapping): + raise SolverStateValidationError(path, "expected mapping", actual=type(value).__name__) + return value + + +def _positive_int(value: Any, path: str, *, allow_zero: bool) -> int: + try: + number = int(value) + except (TypeError, ValueError) as exc: + raise SolverStateValidationError(path, "expected integer", actual=value) from exc + if number < 0 or (number == 0 and not allow_zero): + raise SolverStateValidationError(path, "invalid count", actual=number) + return number + + +def _rank_from_kind(kind: str) -> str: + if kind.endswith("Vector"): + return "vector" + if kind.endswith("Scalar"): + return "scalar" + return "scalar" + + +def _rank_from_value_rank(value_rank: str, path: str) -> str: + if value_rank in {"scalar", "vector"}: + return value_rank + raise SolverStateValidationError(f"{path}.value_rank", "expected scalar or vector rank", actual=value_rank) + + +def _ranked_shape(size: int, rank: str) -> tuple[int, ...]: + if rank == "scalar": + return (size,) + if rank == "vector": + return (size, 3) + raise SolverStateValidationError("rank", "expected scalar or vector rank", actual=rank) + + +def _describe(value: Any) -> Any: + if isinstance(value, np.ndarray): + return {"shape": [int(dim) for dim in value.shape], "dtype": str(value.dtype), "size": int(value.size)} + if isinstance(value, Mapping): + return {str(key): _describe(item) for key, item in value.items()} + if isinstance(value, list): + return [_describe(item) for item in value] + if isinstance(value, tuple): + return [_describe(item) for item in value] + return _json_ready(value) + + +def _json_ready(value: Any) -> Any: + if isinstance(value, np.generic): + return value.item() + if isinstance(value, np.ndarray): + return value.tolist() + 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, (str, int, float, bool)) or value is None: + return value + return repr(value) diff --git a/scripts/build_openfoam_airfrans_subset.sh b/scripts/build_openfoam_airfrans_subset.sh index 73d8163..9e87302 100755 --- a/scripts/build_openfoam_airfrans_subset.sh +++ b/scripts/build_openfoam_airfrans_subset.sh @@ -11,20 +11,59 @@ THIRDPARTY_REPO_URL="${THIRDPARTY_REPO_URL:-https://github.com/OpenFOAM/ThirdPar AIRFRANS_REPO_URL="${AIRFRANS_REPO_URL:-https://github.com/Extrality/AirfRANS.git}" JOBS="${JOBS:-$(nproc)}" -clone_if_missing() { - local url="$1" - local dir="$2" +validate_dependency_tree() { + local dir="$1" + shift - if [[ -d "${ROOT_DIR}/${dir}/.git" ]]; then - printf 'using existing %s\n' "${dir}" - else - git clone --depth 1 "${url}" "${ROOT_DIR}/${dir}" + local root="${ROOT_DIR}/${dir}" + if [[ ! -d "${root}" ]]; then + printf 'dependency %s exists but is not a directory: %s\n' "${dir}" "${root}" >&2 + return 1 + fi + + local -a missing=() + local expected + for expected in "$@"; do + if [[ ! -e "${root}/${expected}" ]]; then + missing+=("${expected}") + fi + done + + if ((${#missing[@]})); then + printf 'dependency %s exists but is incomplete; missing expected path(s):' "${dir}" >&2 + printf ' %s' "${missing[@]}" >&2 + printf '\nMove or remove %s, or replace it with a complete checkout before rerunning.\n' "${root}" >&2 + return 1 fi } -clone_if_missing "${OPENFOAM_REPO_URL}" OpenFOAM-14 -clone_if_missing "${THIRDPARTY_REPO_URL}" ThirdParty-14 -clone_if_missing "${AIRFRANS_REPO_URL}" airfrans +clone_if_missing() { + local url="$1" + local dir="$2" + shift 2 + + local root="${ROOT_DIR}/${dir}" + if [[ -e "${root}" || -L "${root}" ]]; then + validate_dependency_tree "${dir}" "$@" + printf 'using existing %s\n' "${dir}" + else + git clone --depth 1 "${url}" "${root}" + validate_dependency_tree "${dir}" "$@" + fi +} + +clone_if_missing "${OPENFOAM_REPO_URL}" OpenFOAM-14 \ + etc/bashrc \ + wmake/wmake \ + src/OpenFOAM/Make/files \ + applications/solvers/foamRun/Make/files \ + tutorials/incompressibleFluid/venturiTube +clone_if_missing "${THIRDPARTY_REPO_URL}" ThirdParty-14 \ + Allwmake \ + etc/tools +clone_if_missing "${AIRFRANS_REPO_URL}" airfrans \ + README.md \ + dataset.py # Dummy MPI keeps this subset serial and avoids requiring mpicc/OpenMPI for p1. # Use SYSTEMOPENMPI instead only after installing a system MPI development package. diff --git a/scripts/guard_gpu_rans_solver_when_done.sh b/scripts/guard_gpu_rans_solver_when_done.sh new file mode 100755 index 0000000..12d3619 --- /dev/null +++ b/scripts/guard_gpu_rans_solver_when_done.sh @@ -0,0 +1,17 @@ +#!/usr/bin/env bash +set -euo pipefail + +ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd -P)" +NOTES_FILE="${ROOT_DIR}/.loop/notes.md" +STATUS_LINE="" + +if [[ -f "${NOTES_FILE}" ]]; then + IFS= read -r STATUS_LINE < "${NOTES_FILE}" || STATUS_LINE="" +fi + +if [[ "${STATUS_LINE}" != "STATUS: DONE" ]]; then + printf 'Skipping full GPU RANS guard because .loop/notes.md is not STATUS: DONE (%s)\n' "${STATUS_LINE:-missing}" + exit 0 +fi + +exec "${ROOT_DIR}/scripts/verify_gpu_rans_solver.sh" "$@" diff --git a/scripts/prepare_airfrans_stepper_case.py b/scripts/prepare_airfrans_stepper_case.py index 612a40d..675f9e4 100755 --- a/scripts/prepare_airfrans_stepper_case.py +++ b/scripts/prepare_airfrans_stepper_case.py @@ -9,6 +9,7 @@ OpenFOAM/stepper iteration: mesh, initial fields, fvSchemes, and fvSolution. from __future__ import annotations import argparse +import hashlib import json import math import re @@ -21,6 +22,51 @@ RAW_ROOT = ROOT.parent / "airfrans/data/raw/OF_dataset" DEFAULT_SIMULATION = "airFoil2D_SST_93.213_3.79_0.418_0.0_9.665" DEFAULT_SOURCE = RAW_ROOT / DEFAULT_SIMULATION DEFAULT_DEST = ROOT / "tmp/airfrans_stepper_case" / f"{DEFAULT_SIMULATION}_v14" +METADATA_FILENAME = "airfrans_case_metadata.json" +MANIFEST_FILENAME = "airfrans_case_manifest.json" +REQUIRED_COMPARISON_FIELDS = ("U", "p", "phi", "nut", "k", "omega") +REQUIRED_SOURCE_FILES = ( + "system/controlDict", + "system/fvSchemes", + "system/fvSolution", + "system/blockMeshDict", + "0/U", + "0/p", + "0/nut", + "0/k", + "0/omega", + "constant/transportProperties", + "constant/turbulenceProperties", + "constant/polyMesh/boundary", + "constant/polyMesh/points.gz", + "constant/polyMesh/faces.gz", + "constant/polyMesh/owner.gz", + "constant/polyMesh/neighbour.gz", +) +SOURCE_PARAMETER_FILES = ( + "system/controlDict", +) +COPIED_SOURCE_FILES = ( + "system/fvSchemes", + "system/fvSolution", + "system/blockMeshDict", + "0/U", + "0/p", + "0/nut", + "0/k", + "0/omega", + "constant/transportProperties", + "constant/turbulenceProperties", +) +TOPOLOGY_EVIDENCE_KEYS = ( + "n_points", + "n_faces", + "n_internal_faces", + "n_cells", + "patches", + "topology_sha256", + "geometry_sha256", +) FLOAT_RE = re.compile(r"[-+]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][-+]?\d+)?") @@ -198,6 +244,109 @@ def metadata_from_source(src: Path, migrated_end_time: int) -> AirfransCaseMetad migrated_end_time=migrated_end_time, ) +def sha256_file(path: Path) -> str: + digest = hashlib.sha256() + with path.open("rb") as handle: + for chunk in iter(lambda: handle.read(1024 * 1024), b""): + digest.update(chunk) + return digest.hexdigest() + + +def file_record(root: Path, relative: str) -> dict[str, object]: + path = root / relative + stat = path.stat() + return { + "path": relative, + "size": stat.st_size, + "sha256": sha256_file(path), + } + + +def case_file_inventory(root: Path, *, exclude: tuple[str, ...] = ()) -> list[dict[str, object]]: + excluded = set(exclude) + records = [] + for path in sorted(root.rglob("*")): + if not path.is_file(): + continue + relative = path.relative_to(root).as_posix() + if relative in excluded: + continue + records.append(file_record(root, relative)) + return records + +def source_input_inventory(src: Path) -> list[dict[str, object]]: + seen: set[str] = set() + records: list[dict[str, object]] = [] + + def add(relative: str) -> None: + if relative in seen: + return + seen.add(relative) + records.append(file_record(src, relative)) + + for relative in SOURCE_PARAMETER_FILES + COPIED_SOURCE_FILES: + add(relative) + for path in sorted((src / "constant/polyMesh").rglob("*")): + if path.is_file(): + add(path.relative_to(src).as_posix()) + return records + + + +def inventory_digest(records: list[dict[str, object]]) -> str: + digest = hashlib.sha256() + for record in records: + payload = json.dumps(record, sort_keys=True, separators=(",", ":")).encode("utf-8") + digest.update(len(payload).to_bytes(8, "little")) + digest.update(payload) + return digest.hexdigest() + + +def build_case_manifest(src: Path, dst: Path, meta: AirfransCaseMetadata) -> dict[str, object]: + source_records = source_input_inventory(src) + prepared_records = case_file_inventory(dst, exclude=(MANIFEST_FILENAME,)) + return { + "schema_version": 1, + "simulation": meta.simulation, + "source_case": { + "path": str(src), + "read_only": True, + "required_files": source_records, + "required_files_sha256": inventory_digest(source_records), + }, + "prepared_case": { + "path": str(dst), + "openfoam_version": 14, + "migrated_end_time": meta.migrated_end_time, + "files": prepared_records, + "files_sha256": inventory_digest(prepared_records), + }, + "comparison_contract": { + "fields": list(REQUIRED_COMPARISON_FIELDS), + "oracle_output": { + "producer": "foamRun -solver incompressibleFluid -noFunctionObjects", + "time": str(meta.migrated_end_time), + }, + "repository_outputs": { + "run_one": "foam_stepper run_one_pimple_iteration", + "split": "foam_stepper split solver stages", + "time": str(meta.migrated_end_time), + }, + "mesh_identity_evidence": list(TOPOLOGY_EVIDENCE_KEYS), + }, + "metadata": asdict(meta), + } + + +def write_case_manifest(src: Path, dst: Path, meta: AirfransCaseMetadata) -> dict[str, object]: + manifest = build_case_manifest(src, dst, meta) + (dst / MANIFEST_FILENAME).write_text(json.dumps(manifest, indent=2, sort_keys=True) + "\n") + return manifest + + +def load_case_manifest(case: Path) -> dict[str, object]: + return json.loads((case / MANIFEST_FILENAME).read_text()) + def prepare_case(src: Path, dst: Path, *, end_time: int = 1) -> AirfransCaseMetadata: src = src.resolve() @@ -205,24 +354,7 @@ def prepare_case(src: Path, dst: Path, *, end_time: int = 1) -> AirfransCaseMeta if not src.exists(): raise FileNotFoundError(src) - required = [ - "system/fvSchemes", - "system/fvSolution", - "system/blockMeshDict", - "0/U", - "0/p", - "0/nut", - "0/k", - "0/omega", - "constant/transportProperties", - "constant/turbulenceProperties", - "constant/polyMesh/boundary", - "constant/polyMesh/points.gz", - "constant/polyMesh/faces.gz", - "constant/polyMesh/owner.gz", - "constant/polyMesh/neighbour.gz", - ] - missing = [relative for relative in required if not (src / relative).exists()] + missing = [relative for relative in REQUIRED_SOURCE_FILES if not (src / relative).exists()] if missing: raise FileNotFoundError(f"missing required AirfRANS files: {missing}") @@ -231,18 +363,7 @@ def prepare_case(src: Path, dst: Path, *, end_time: int = 1) -> AirfransCaseMeta dst.mkdir(parents=True) copy_required_tree(src, dst, "constant/polyMesh") - for relative in [ - "system/fvSchemes", - "system/fvSolution", - "system/blockMeshDict", - "0/U", - "0/p", - "0/nut", - "0/k", - "0/omega", - "constant/transportProperties", - "constant/turbulenceProperties", - ]: + for relative in COPIED_SOURCE_FILES: copy_required_file(src, dst, relative) meta = metadata_from_source(src, end_time) @@ -254,7 +375,8 @@ def prepare_case(src: Path, dst: Path, *, end_time: int = 1) -> AirfransCaseMeta shutil.move(dst / "constant/transportProperties", dst / "constant/transportProperties.v2112") shutil.move(dst / "constant/turbulenceProperties", dst / "constant/turbulenceProperties.v2112") - (dst / "airfrans_case_metadata.json").write_text(json.dumps(asdict(meta), indent=2, sort_keys=True) + "\n") + (dst / METADATA_FILENAME).write_text(json.dumps(asdict(meta), indent=2, sort_keys=True) + "\n") + write_case_manifest(src, dst, meta) return meta @@ -268,6 +390,7 @@ def main() -> None: meta = prepare_case(args.source, args.dest, end_time=args.end_time) print(json.dumps(asdict(meta), indent=2, sort_keys=True)) print(f"prepared={args.dest}") + print(f"manifest={args.dest / MANIFEST_FILENAME}") if __name__ == "__main__": diff --git a/scripts/verify_airfrans_stepper.py b/scripts/verify_airfrans_stepper.py index 54c8cdb..fe30c2d 100755 --- a/scripts/verify_airfrans_stepper.py +++ b/scripts/verify_airfrans_stepper.py @@ -14,7 +14,7 @@ import subprocess import sys import time import traceback -from collections.abc import Mapping +from collections.abc import Iterable, Mapping from pathlib import Path from typing import Any @@ -33,18 +33,33 @@ _drop_ambient_pythonpath() import numpy as np from openfoam_env import apply_openfoam_env, openfoam_env -from prepare_airfrans_stepper_case import DEFAULT_SOURCE, prepare_case +from prepare_airfrans_stepper_case import ( + DEFAULT_SOURCE, + REQUIRED_COMPARISON_FIELDS, + load_case_manifest, + prepare_case, + write_case_manifest, +) ROOT = Path(__file__).resolve().parents[1] DEFAULT_WORK = ROOT / "tmp/airfrans_stepper_verify" PRIMARY_FIELDS = ("U", "p", "phi") TURBULENCE_FIELDS = ("nut", "k", "omega") -REQUIRED_FIELDS = PRIMARY_FIELDS + TURBULENCE_FIELDS +REQUIRED_FIELDS = REQUIRED_COMPARISON_FIELDS +BACKEND_CHOICES = ("auto", "cpu", "gpu") +FULL_GPU_RANS_GUARD = "scripts/verify_gpu_rans_solver.sh" +GPU_PRIMITIVE_PROOF = "scripts/verify_gpu_algorithm.sh" +GPU_PRIMITIVE_NAME = "cell_flux_imbalance" +AIRFRANS_VERIFIER_SCOPE = ( + "AirfRANS/OpenFOAM stepper parity verifier; --backend gpu is a GPU-owned " + "RANS solver contract with no CPU fallback or primitive-only acceptance" +) LOADING_FAILURE = "loading_or_parsing_failure" ORACLE_FAILURE = "openfoam_oracle_failure" STEPPER_FAILURE = "stepper_execution_failure" COMPARISON_FAILURE = "numerical_comparison_failure" +BACKEND_FAILURE = "backend_execution_failure" INTERNAL_FAILURE = "internal_harness_failure" EXIT_CODES = { @@ -52,9 +67,133 @@ EXIT_CODES = { ORACLE_FAILURE: 3, STEPPER_FAILURE: 4, COMPARISON_FAILURE: 5, - INTERNAL_FAILURE: 6, + BACKEND_FAILURE: 6, + INTERNAL_FAILURE: 7, } +STAGE_OBSERVABILITY_GROUPS = ( + { + "name": "momentum_assembly", + "split_stages": ("assemble_momentum_terms", "assemble_UEqn"), + "run_one_stages": ("assemble_UEqn",), + "split_outputs": { + "assemble_momentum_terms": ("terms",), + "assemble_UEqn": ("UEqn",), + }, + }, + { + "name": "pressure_assembly", + "split_stages": ("compute_pressure_inputs", "assemble_pEqn"), + "run_one_stages": ("compute_pressure_inputs", "assemble_pEqn"), + "split_outputs": { + "compute_pressure_inputs": ("HbyA", "phiHbyA", "rAU"), + "assemble_pEqn": ("pEqn",), + }, + }, + { + "name": "linear_solve_results", + "split_stages": ("solve_UEqn", "solve_pEqn"), + "run_one_stages": ("solve_UEqn", "solve_pEqn"), + "split_outputs": { + "solve_UEqn": ("performance", "field_after"), + "solve_pEqn": ("performance", "p", "phi"), + }, + }, + { + "name": "final_correction", + "split_stages": ("correct_velocity_pressure_flux",), + "run_one_stages": ("update_phi_from_pEqn_flux", "correct_velocity_pressure_flux"), + "split_outputs": { + "correct_velocity_pressure_flux": ("U", "p", "phi"), + }, + }, + { + "name": "turbulence_updates", + "split_stages": ("momentum_transport_predict", "momentum_transport_correct"), + "run_one_stages": ("momentum_transport_predict", "momentum_transport_correct"), + "split_outputs": { + "momentum_transport_predict": ("case_path", "solver_name"), + "momentum_transport_correct": ("U", "p", "phi", "nut", "k", "omega"), + }, + }, +) + +FIELD_COMPARISON_ATTRIBUTION = { + "U": ("final_correction", "linear_solve_results", "momentum_assembly", "turbulence_updates"), + "p": ("linear_solve_results", "pressure_assembly", "final_correction"), + "phi": ("final_correction", "linear_solve_results", "pressure_assembly"), + "nut": ("turbulence_updates",), + "k": ("turbulence_updates",), + "omega": ("turbulence_updates",), +} + +GPU_RUNTIME_UNAVAILABLE = "gpu_runtime_unavailable" +GPU_NUMERICAL_MISMATCH = "gpu_numerical_mismatch" +GPU_INPUT_SCHEMA_VERSION = 1 + +DIFFERENTIABILITY_VARIABLES = { + "fields": REQUIRED_FIELDS, + "mesh_state": ("points", "faces", "owner", "neighbour", "cells", "V", "C", "Cf", "Sf", "magSf", "patches"), + "case_parameters": ("Uinf", "nu", "rhoInf", "angle_of_attack", "liftDir", "dragDir", "solver_tolerances"), + "losses": ("per_field_linf_error", "per_field_l2_error", "oracle_parity_allclose"), +} + +NON_DIFFERENTIABLE_BOUNDARIES = ( + { + "name": "openfoam_case_io", + "category": "io", + "reason": "OpenFOAM dictionaries, mesh files, and field files are parsed and written as discrete external data.", + }, + { + "name": "mesh_topology_and_patch_addressing", + "category": "topology", + "reason": "Face/cell connectivity, owner/neighbour addressing, and boundary patch membership are discrete indices.", + }, + { + "name": "boundary_condition_selection", + "category": "discrete_boundary_choice", + "reason": "Patch types, coupled/constraint flags, and boundary condition dictionaries select code paths.", + }, + { + "name": "simple_pimple_and_linear_solver_control", + "category": "convergence_control", + "reason": "Iteration counts, stopping criteria, and solver convergence branches are discrete control flow.", + }, + { + "name": "turbulence_bounding_and_limiters", + "category": "limiter", + "reason": "k/omega/nut bounding and model limiters introduce clipping or branch-dependent updates.", + }, + { + "name": "openfoam_cpp_backend_calls", + "category": "unsupported_solver_operation", + "reason": "The current registered backend executes OpenFOAM C++ operations and exposes no AD, JVP, or VJP API.", + }, +) + +CUSTOM_GRADIENT_REQUIREMENTS = ( + { + "operation": "sparse_linear_solves", + "stage_groups": ("linear_solve_results",), + "reason": "Implicit sparse solves require an adjoint or custom VJP rather than differentiating solver iterations blindly.", + }, + { + "operation": "pressure_velocity_correction", + "stage_groups": ("pressure_assembly", "final_correction"), + "reason": "The coupled pressure/flux/velocity correction needs a consistent custom gradient for the assembled operators.", + }, + { + "operation": "turbulence_model_correctors_and_limiters", + "stage_groups": ("turbulence_updates",), + "reason": "Model correctors, wall functions, and bounding need explicit subgradient or smoothing choices.", + }, + { + "operation": "mesh_and_boundary_discrete_choices", + "stage_groups": (), + "reason": "Topology and patch-type changes are not differentiable; only fixed-topology numeric values can be checked.", + }, +) + class HarnessError(Exception): """Categorized verifier failure that can be serialized into the report.""" @@ -88,6 +227,19 @@ class HarnessError(Exception): } return json_ready(out) +class BackendExecutionError(Exception): + """Backend selection or execution failure before numerical comparison.""" + + def __init__(self, step: str, message: str, *, details: Mapping[str, Any] | None = None) -> None: + super().__init__(message) + self.step = step + self.details = dict(details or {}) + + def to_harness_error(self) -> HarnessError: + return HarnessError(BACKEND_FAILURE, self.step, str(self), details=self.details, cause=self) + + + def json_ready(value: Any) -> Any: """Convert report values into strict JSON-compatible data.""" @@ -299,6 +451,73 @@ def transform_summary(result: Any) -> dict[str, Any]: "outputs": summarize_value(getattr(result, "outputs", {})), } +def graph_stage_summaries(result: Any) -> list[dict[str, Any]]: + return [transform_summary(entry) for entry in getattr(result, "outputs", {}).get("graph", [])] + + +def stage_observability_report(mode: str, stages: list[dict[str, Any]], *, evidence: Mapping[str, Any]) -> dict[str, Any]: + stage_by_name = {stage["name"]: stage for stage in stages} + failures: list[dict[str, Any]] = [] + groups = [] + for group in STAGE_OBSERVABILITY_GROUPS: + group_name = group["name"] + required_stages = tuple(group[f"{mode}_stages"]) + missing_stages = [name for name in required_stages if name not in stage_by_name] + output_requirements = group.get(f"{mode}_outputs", {}) + output_checks = {} + for stage_name, required_outputs in output_requirements.items(): + stage = stage_by_name.get(stage_name) + if stage is None: + continue + outputs = stage.get("outputs", {}) + output_keys = set(outputs) if isinstance(outputs, Mapping) else set() + missing_outputs = [name for name in required_outputs if name not in output_keys] + output_checks[stage_name] = { + "required": list(required_outputs), + "available": sorted(output_keys), + "missing": missing_outputs, + } + if missing_outputs: + failures.append( + { + "path": f"stage_observability.{mode}.{group_name}.{stage_name}.outputs", + "message": "missing inspectable stage output", + "expected": list(required_outputs), + "actual": sorted(output_keys), + } + ) + if missing_stages: + failures.append( + { + "path": f"stage_observability.{mode}.{group_name}.stages", + "message": "missing required solver stage", + "expected": list(required_stages), + "actual": list(stage_by_name), + } + ) + groups.append( + { + "name": group_name, + "required_stages": list(required_stages), + "observed": not missing_stages and not any(check["missing"] for check in output_checks.values()), + "missing_stages": missing_stages, + "output_checks": output_checks, + } + ) + if failures: + raise HarnessError( + STEPPER_FAILURE, + f"{mode}_stage_observability", + f"{mode} solver-stage observability is incomplete", + details={"failures": failures}, + ) + return { + "mode": mode, + "graph": [stage["name"] for stage in stages], + "groups": groups, + "evidence": json_ready(evidence), + } + def selected_field_summaries(fields: Mapping[str, Any]) -> dict[str, Any]: return {name: field_summary(fields[name]) for name in REQUIRED_FIELDS if name in fields} @@ -385,17 +604,59 @@ def compare_mesh_identity(mode: str, actual: Mapping[str, Any], expected: Mappin return report, [{"mode": mode, "kind": "mesh_identity", "differences": differences}] +def field_comparison_attribution( + mode: str, + field: str, + reason: str, + stage_graph: Iterable[str] | None, +) -> dict[str, Any]: + observed_stages = set(stage_graph or ()) + candidates = [] + for group_name in FIELD_COMPARISON_ATTRIBUTION.get(field, ()): + group = next(item for item in STAGE_OBSERVABILITY_GROUPS if item["name"] == group_name) + stages = list(group.get(f"{mode}_stages", ())) + present_stages = [stage for stage in stages if not observed_stages or stage in observed_stages] + candidates.append( + { + "group": group_name, + "stages": present_stages, + "evidence_path": f"stage_observability.{mode}.{group_name}", + } + ) + + likely = candidates[0] if candidates else {"group": "unknown", "stages": [], "evidence_path": None} + out = { + "field": field, + "reason": reason, + "likely_stage_group": likely["group"], + "likely_stages": likely["stages"], + "candidate_stage_groups": candidates, + "evidence_path": likely["evidence_path"], + "basis": "field-to-stage dependency map for the selected solver graph", + } + if reason in {"missing_field", "shape_mismatch"}: + out["precondition_failure"] = "field_export_or_state_mapping" + out["field_evidence_path"] = f"modes.{mode}.fields.{field}" + out["basis"] = "field is absent or has the wrong topology; inspect export/state mapping first, then the mapped solver stages" + return out + + def field_compare_report(name: str, actual_field: Any, expected_field: Any, *, rtol: float, atol: float) -> tuple[dict[str, Any], dict[str, Any] | None]: actual = np.asarray(actual_field.internal) expected = np.asarray(expected_field.internal) + actual_entity_kind = getattr(actual_field, "entity_kind", "") + expected_entity_kind = getattr(expected_field, "entity_kind", "") report: dict[str, Any] = { "field": name, "actual_shape": array_shape(actual), "expected_shape": array_shape(expected), - "actual_entity_kind": getattr(actual_field, "entity_kind", ""), - "expected_entity_kind": getattr(expected_field, "entity_kind", ""), + "shape_matches": actual.shape == expected.shape, + "actual_entity_kind": actual_entity_kind, + "expected_entity_kind": expected_entity_kind, + "entity_kind_matches": actual_entity_kind == expected_entity_kind, "rtol": rtol, "atol": atol, + "comparison": "numpy.allclose(equal_nan=False)", } if actual.shape != expected.shape: @@ -467,6 +728,7 @@ def compare_fields( *, rtol: float, atol: float, + stage_graph: Iterable[str] | None = None, ) -> tuple[dict[str, Any], list[dict[str, Any]]]: report: dict[str, Any] = {} mismatches: list[dict[str, Any]] = [] @@ -480,12 +742,16 @@ def compare_fields( "actual_fields": sorted(actual), "expected_fields": sorted(expected), } - report[name] = {"field": name, "allclose": False, **missing} - mismatches.append({"mode": mode, **missing}) + attribution = field_comparison_attribution(mode, name, "missing_field", stage_graph) + report[name] = {"field": name, "allclose": False, "attribution": attribution, **missing} + mismatches.append({"mode": mode, "attribution": attribution, **missing}) continue field_report, mismatch = field_compare_report(name, actual[name], expected[name], rtol=rtol, atol=atol) report[name] = field_report if mismatch is not None: + attribution = field_comparison_attribution(mode, name, mismatch["reason"], stage_graph) + field_report["attribution"] = attribution + mismatch["attribution"] = attribution mismatches.append({"mode": mode, **mismatch}) return report, mismatches @@ -539,30 +805,199 @@ def import_foam() -> Any: return foam +def clone_prepared_case(prepared: Path, dest: Path, source: Path, metadata: Any) -> dict[str, Any]: + if dest.exists(): + shutil.rmtree(dest) + shutil.copytree(prepared, dest, copy_function=shutil.copy2) + return write_case_manifest(source, dest, metadata) + + +def prepared_case_digest(manifest: Mapping[str, Any]) -> str: + prepared = manifest.get("prepared_case", {}) + if not isinstance(prepared, Mapping): + return "" + return str(prepared.get("files_sha256", "")) + + +def prepared_file_count(manifest: Mapping[str, Any]) -> int: + prepared = manifest.get("prepared_case", {}) + files = prepared.get("files", []) if isinstance(prepared, Mapping) else [] + return len(files) if isinstance(files, list) else 0 + + +def compare_prepared_case_manifests(manifests: Mapping[str, Mapping[str, Any]]) -> dict[str, dict[str, Any]]: + baseline = prepared_case_digest(manifests["prepared"]) + return { + role: { + "matches_prepared": prepared_case_digest(manifest) == baseline, + "files_sha256": prepared_case_digest(manifest), + "file_count": prepared_file_count(manifest), + } + for role, manifest in manifests.items() + } + + def prepare_work_cases(source: Path, work: Path, *, include_split: bool) -> dict[str, Any]: if work.exists(): shutil.rmtree(work) work.mkdir(parents=True) + prepared_case = work / "prepared_case" oracle = work / "oracle_case" run_one = work / "run_one_case" split = work / "split_case" if include_split else None - oracle_meta = prepare_case(source, oracle, end_time=1) - run_one_meta = prepare_case(source, run_one, end_time=1) - split_meta = prepare_case(source, split, end_time=1) if split is not None else None + prepared_meta = prepare_case(source, prepared_case, end_time=1) + manifests: dict[str, Any] = {"prepared": load_case_manifest(prepared_case)} + manifests["oracle"] = clone_prepared_case(prepared_case, oracle, source, prepared_meta) + manifests["run_one"] = clone_prepared_case(prepared_case, run_one, source, prepared_meta) + if split is not None: + manifests["split"] = clone_prepared_case(prepared_case, split, source, prepared_meta) return { + "prepared_case": prepared_case, "oracle_case": oracle, "run_one_case": run_one, "split_case": split, "metadata": { - "oracle": oracle_meta, - "run_one": run_one_meta, - "split": split_meta, + "prepared": prepared_meta, + "oracle": prepared_meta, + "run_one": prepared_meta, + "split": prepared_meta if split is not None else None, }, + "manifests": manifests, + "prepared_case_identity": compare_prepared_case_manifests(manifests), } +def _gpu_backend_module() -> Any: + try: + from foam_stepper.gpu import backend as gpu_backend + except Exception as exc: + raise BackendExecutionError( + "select_backend", + "GPU runtime unavailable: failed to import reusable foam_stepper.gpu backend", + details={ + "requested": "gpu", + "failure_kind": GPU_RUNTIME_UNAVAILABLE, + "provider": "quadrants_cuda_rans_solver", + "used_cpu_fallback": False, + "cause": {"type": type(exc).__name__, "message": str(exc)}, + }, + ) from exc + return gpu_backend + + +def _wrap_gpu_backend_error(exc: BaseException) -> BackendExecutionError: + return BackendExecutionError( + str(getattr(exc, "step", "execute_gpu_solver_contract")), + str(exc), + details=dict(getattr(exc, "details", {}) or {}), + ) + + +def _call_gpu_backend(function_name: str, *args: Any, **kwargs: Any) -> Any: + gpu_backend = _gpu_backend_module() + try: + return getattr(gpu_backend, function_name)(*args, **kwargs) + except gpu_backend.BackendExecutionError as exc: + raise _wrap_gpu_backend_error(exc) from exc + + +def prepare_gpu_solver_inputs(foam: Any, stepper: Any, case: Path, prepared: Mapping[str, Any], backend: Mapping[str, Any]) -> dict[str, Any]: + return _call_gpu_backend("prepare_gpu_solver_inputs", foam, stepper, case, prepared, backend) + + +def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: Path) -> dict[str, Any]: + return _call_gpu_backend("run_gpu_solver_stage_smoke", stepper, backend, case) + + +def gpu_solver_stage_report(gpu_run: Mapping[str, Any]) -> dict[str, Any]: + return _call_gpu_backend("gpu_solver_stage_report", gpu_run) + + +def gpu_solver_input_blocker_summary(gpu_inputs: Mapping[str, Any]) -> dict[str, Any]: + return _call_gpu_backend("gpu_solver_input_blocker_summary", gpu_inputs) + + +def gpu_solver_contract_blocker(backend: Mapping[str, Any], case: Path) -> BackendExecutionError | None: + gpu_backend = _gpu_backend_module() + try: + blocker = gpu_backend.gpu_solver_contract_blocker(backend, case) + except gpu_backend.BackendExecutionError as exc: + raise _wrap_gpu_backend_error(exc) from exc + if blocker is None: + return None + return _wrap_gpu_backend_error(blocker) + + +def record_gpu_contract_failure(report: dict[str, Any], blocker: BackendExecutionError, gpu_inputs: Mapping[str, Any] | None = None) -> None: + if gpu_inputs is not None: + blocker.details["gpu_solver_inputs"] = gpu_solver_input_blocker_summary(gpu_inputs) + report["backend"]["contract_preflight"] = json_ready( + { + "status": "failed", + "step": blocker.step, + "failure_kind": blocker.details.get("failure_kind"), + "checked_before_oracle": True, + "gpu_inputs_prepared": gpu_inputs is not None, + "reason": "GPU requests must fail as incomplete before any CPU/OpenFOAM solver path can be substituted as the backend.", + } + ) + + raise blocker.to_harness_error() from blocker + + +def select_execution_backend(requested: str) -> dict[str, Any]: + if requested not in BACKEND_CHOICES: + raise BackendExecutionError( + "select_backend", + f"unknown backend {requested!r}", + details={"requested": requested, "choices": list(BACKEND_CHOICES)}, + ) + + if requested == "gpu": + return _call_gpu_backend("select_execution_backend", requested) + + return { + "requested": requested, + "selected": "cpu", + "device": "host", + "provider": "foam_stepper_cpu", + "available_backends": ["cpu"], + "gpu_device_present": any(Path(path).exists() for path in ("/dev/nvidia0", "/dev/dri/renderD128")), + "capabilities": [ + "foam_stepper_python_bridge", + "openfoam_case_loader", + "oracle_comparison", + ], + "used_cpu_fallback": False, + } + + +def run_backend_iteration(stepper: Any, backend: Mapping[str, Any], case: Path) -> dict[str, Any]: + selected = backend.get("selected") + if selected == "gpu": + return _call_gpu_backend("run_backend_iteration", stepper, backend, case) + if selected != "cpu": + raise BackendExecutionError( + "execute_backend", + f"backend {selected!r} is not executable", + details={"backend": backend, "case": case}, + ) + result = stepper.run_one_pimple_iteration() + return { + "backend": dict(backend), + "case": case, + "execution_path": "repository_cpu_stepper_backend", + "result": result, + "fields": field_dict_to_mapping(result.outputs["fields"]), + } + + +def run_gpu_split_iteration(foam: Any, stepper: Any, backend: Mapping[str, Any], case: Path) -> dict[str, Any]: + return _call_gpu_backend("run_gpu_split_iteration", foam, stepper, backend, case) + + def make_stepper(foam: Any, case: Path, label: str) -> Any: try: return foam.Case(case).make_stepper() @@ -601,6 +1036,38 @@ def read_fields(stepper: Any, label: str) -> dict[str, Any]: cause=exc, ) from exc +def validation_details(exc: BaseException) -> dict[str, Any]: + if hasattr(exc, "to_dict"): + return json_ready(exc.to_dict()) + return {"type": type(exc).__name__, "message": str(exc)} + + +def export_solver_state_checked(foam: Any, stepper: Any, label: str, fields: Mapping[str, Any] | None = None) -> dict[str, Any]: + try: + field_source = fields if fields is not None else stepper.fields() + return foam.export_solver_state(stepper.mesh(), field_source, required_fields=REQUIRED_FIELDS) + except Exception as exc: + raise HarnessError( + LOADING_FAILURE, + f"export_{label}_state", + f"failed to export explicit {label} solver state", + details=validation_details(exc), + cause=exc, + ) from exc + + +def export_matrix_state_checked(foam: Any, matrix: Any, solver_state: Mapping[str, Any], label: str) -> dict[str, Any]: + try: + return foam.export_matrix_state(matrix, mesh_state=solver_state) + except Exception as exc: + raise HarnessError( + STEPPER_FAILURE, + f"export_{label}_matrix_state", + f"failed to export explicit {label} matrix state", + details=validation_details(exc), + cause=exc, + ) from exc + def checked_step(stages: list[dict[str, Any]], name: str, fn: Any) -> Any: try: @@ -616,7 +1083,7 @@ def checked_step(stages: list[dict[str, Any]], name: str, fn: Any) -> Any: return result -def run_split_iteration(stepper: Any) -> dict[str, Any]: +def run_split_iteration(foam: Any, stepper: Any) -> dict[str, Any]: stages: list[dict[str, Any]] = [] checked_step(stages, "pre_solve", stepper.pre_solve) checked_step(stages, "advance_time", stepper.advance_time) @@ -655,16 +1122,36 @@ def run_split_iteration(stepper: Any) -> dict[str, Any]: cause=exc, ) from exc + solver_state = export_solver_state_checked(foam, stepper, "split", fields) + momentum_terms = [term.get("name", "") for term in terms.outputs.get("terms", [])] UEqn_matrix = UEqn.outputs["UEqn"] pEqn_matrix = pEqn.outputs["pEqn"] + matrix_states = { + "UEqn": foam.describe_matrix_state(export_matrix_state_checked(foam, UEqn_matrix, solver_state, "UEqn")), + "pEqn": foam.describe_matrix_state(export_matrix_state_checked(foam, pEqn_matrix, solver_state, "pEqn")), + } + observability = stage_observability_report( + "split", + stages, + evidence={ + "split_step_execution": True, + "momentum_terms": momentum_terms, + "matrix_states": sorted(matrix_states), + "turbulence_fields": visible_turbulence_fields(fields), + }, + ) return { "fields": fields, + "state": solver_state, + "state_summary": foam.describe_solver_state(solver_state), "stages": stages, "graph": [stage["name"] for stage in stages], "momentum_terms": momentum_terms, "UEqn": matrix_summary(UEqn_matrix), "pEqn": matrix_summary(pEqn_matrix), + "matrix_states": matrix_states, + "observability": observability, } @@ -672,6 +1159,229 @@ def visible_turbulence_fields(fields: Mapping[str, Any]) -> list[str]: return [name for name in TURBULENCE_FIELDS if name in fields] +def differentiability_report(backend: Mapping[str, Any], *, case_name: str) -> dict[str, Any]: + capabilities = tuple(backend.get("capabilities", ())) + autodiff_available = any(capability in capabilities for capability in ("autodiff", "jvp", "vjp", "gradient")) + status = "differentiable" if autodiff_available else "not_differentiable" + return { + "status": status, + "backend": { + "requested": backend.get("requested"), + "selected": backend.get("selected"), + "provider": backend.get("provider"), + "capabilities": list(capabilities), + }, + "variables": DIFFERENTIABILITY_VARIABLES, + "claims": [ + { + "name": "current_migrated_solver_path", + "differentiable": autodiff_available, + "variables": DIFFERENTIABILITY_VARIABLES, + "applies_to": { + "modes": ["run_one", "split"], + "fields": list(REQUIRED_FIELDS), + "parameters": list(DIFFERENTIABILITY_VARIABLES["case_parameters"]), + "losses": list(DIFFERENTIABILITY_VARIABLES["losses"]), + }, + "evidence": { + "backend_capabilities": list(capabilities), + "reason": "no registered backend capability exposes AD/JVP/VJP for solver updates" + if not autodiff_available + else "backend advertises differentiability capabilities", + }, + } + ], + "non_differentiable": list(NON_DIFFERENTIABLE_BOUNDARIES), + "custom_gradient_required": list(CUSTOM_GRADIENT_REQUIREMENTS), + "checks": [ + { + "name": "solver_sensitivity_against_finite_difference", + "status": "not_run", + "fixed_problem_setup": { + "case": case_name, + "fields": list(REQUIRED_FIELDS), + "mesh_topology": "fixed prepared AirfRANS mesh", + "backend": backend.get("selected"), + }, + "computed_sensitivity": None, + "independent_numerical_check": None, + "reason": "no differentiable backend or custom gradient is currently available, so no computed sensitivity is claimed", + } + ], + } + + +def command_timing_summary(command: Mapping[str, Any]) -> dict[str, Any]: + cmd = list(command.get("cmd") or ()) + executable = Path(str(cmd[0])).name if cmd else None + log_path = str(command.get("log_path", "")) + role = "openfoam_command" + if executable == "checkMesh": + role = "oracle_mesh_check" + elif executable == "foamRun" and "oracle" in log_path: + role = "oracle_solver" + return { + "role": role, + "executable": executable, + "cmd": cmd, + "duration_seconds": command.get("duration_seconds"), + "timeout_seconds": command.get("timeout_seconds"), + "log_path": command.get("log_path"), + "returncode": command.get("returncode"), + } + + +def gpu_timing_summary(source: Mapping[str, Any], *, label: str) -> dict[str, Any]: + profiler = source.get("profiler", {}) if isinstance(source.get("profiler"), Mapping) else {} + timing = source.get("timing", {}) if isinstance(source.get("timing"), Mapping) else {} + return { + "label": label, + "execution_path": source.get("execution_path"), + "timing": timing, + "profiler": { + "available": profiler.get("available"), + "profile_record_count": profiler.get("profile_record_count"), + "device_time_ms_total": profiler.get("device_time_ms_total"), + "expected_kernel_entrypoints": profiler.get("expected_kernel_entrypoints"), + "generated_kernel_names": profiler.get("generated_kernel_names"), + }, + } + + +def build_timing_evidence(report: Mapping[str, Any]) -> dict[str, Any]: + commands = [command_timing_summary(command) for command in report.get("commands", [])] + oracle_solver = next((command for command in commands if command.get("role") == "oracle_solver"), None) + gpu_sources: dict[str, Any] = {} + gpu_stage_smoke = report.get("gpu_solver_stages") + if isinstance(gpu_stage_smoke, Mapping) and gpu_stage_smoke.get("status") == "executed": + gpu_sources["stage_preflight"] = gpu_timing_summary(gpu_stage_smoke, label="stage_preflight") + for mode in ("run_one", "split"): + mode_gpu = report.get("modes", {}).get(mode, {}).get("gpu_solver", {}) + if isinstance(mode_gpu, Mapping) and mode_gpu.get("status") == "executed": + gpu_sources[mode] = gpu_timing_summary(mode_gpu, label=mode) + + gpu_solver_timing_ok = any( + isinstance(source.get("timing"), Mapping) and source["timing"].get("gpu_solver_wall_seconds") is not None + for source in gpu_sources.values() + ) + gpu_profiler_ok = any( + isinstance(source.get("profiler"), Mapping) and source["profiler"].get("available") is True and source["profiler"].get("profile_record_count", 0) > 0 + for source in gpu_sources.values() + ) + openfoam_timing_ok = oracle_solver is not None and oracle_solver.get("duration_seconds") is not None + return { + "schema_version": 1, + "status": "reported" if commands or gpu_sources else "not_reported", + "openfoam_reference": { + "oracle_solver_duration_seconds": None if oracle_solver is None else oracle_solver.get("duration_seconds"), + "commands": commands, + }, + "gpu_solver": gpu_sources, + "criteria": { + "openfoam_reference_timing": openfoam_timing_ok, + "gpu_solver_wall_timing": gpu_solver_timing_ok, + "gpu_profiler_timing": gpu_profiler_ok, + }, + } + +def enabled_mode_names(report: Mapping[str, Any]) -> list[str]: + return [mode for mode in ("run_one", "split") if report.get("modes", {}).get(mode, {}).get("enabled", False)] + + +def comparison_regression_summary(comparisons: Mapping[str, Any]) -> dict[str, Any]: + out: dict[str, Any] = {} + for field in REQUIRED_FIELDS: + data = comparisons.get(field, {}) + out[field] = { + "actual_shape": data.get("actual_shape"), + "expected_shape": data.get("expected_shape"), + "shape_matches": data.get("shape_matches"), + "allclose": data.get("allclose"), + "max_abs": data.get("max_abs"), + "mean_abs": data.get("mean_abs"), + "tolerance_at_max": data.get("tolerance_at_max"), + "location": data.get("location"), + } + return out + + +def build_verifier_evidence(report: Mapping[str, Any]) -> dict[str, Any]: + modes = enabled_mode_names(report) + prepared_identity = report.get("case_preparation", {}).get("prepared_case_identity", {}) + prepared_ok = bool(prepared_identity) and all(identity.get("matches_prepared") for identity in prepared_identity.values()) + + mesh_ok_by_mode = { + mode: bool(report.get("modes", {}).get(mode, {}).get("mesh_comparison", {}).get("matches_oracle")) + for mode in modes + } + comparison_ok_by_mode = {} + observability_ok_by_mode = {} + regression_modes: dict[str, Any] = {} + for mode in modes: + mode_report = report.get("modes", {}).get(mode, {}) + comparisons = mode_report.get("comparisons", {}) + comparison_ok_by_mode[mode] = all( + comparisons.get(field, {}).get("shape_matches") is True and comparisons.get(field, {}).get("allclose") is True + for field in REQUIRED_FIELDS + ) + observability = report.get("stage_observability", {}).get(mode, {}) + groups = observability.get("groups", []) + observability_ok_by_mode[mode] = bool(groups) and all(group.get("observed") is True for group in groups) + regression_modes[mode] = { + "enabled": True, + "case_role": mode, + "mesh_matches_oracle": mesh_ok_by_mode[mode], + "observed_stage_groups": [group.get("name") for group in groups if group.get("observed")], + "fields": comparison_regression_summary(comparisons), + } + + timing = report.get("timing", {}) if isinstance(report.get("timing"), Mapping) else {} + timing_criteria = timing.get("criteria", {}) if isinstance(timing.get("criteria"), Mapping) else {} + if report.get("backend", {}).get("selected") == "gpu": + timing_ok = ( + timing_criteria.get("openfoam_reference_timing") is True + and timing_criteria.get("gpu_solver_wall_timing") is True + and timing_criteria.get("gpu_profiler_timing") is True + ) + else: + timing_ok = timing.get("status") in {"reported", "not_reported"} + + criteria = { + "prepared_case_identity": prepared_ok, + "mesh_identity": bool(mesh_ok_by_mode) and all(mesh_ok_by_mode.values()), + "field_comparisons": bool(comparison_ok_by_mode) and all(comparison_ok_by_mode.values()), + "stage_observability": bool(observability_ok_by_mode) and all(observability_ok_by_mode.values()), + "backend_selected": report.get("backend", {}).get("selected") is not None, + "timing_reported": timing_ok, + "differentiability_reported": report.get("differentiability", {}).get("status") in {"differentiable", "not_differentiable"}, + } + passed = all(criteria.values()) + oracle_mesh = report.get("mesh_identity", {}).get("oracle", {}) + prepared_digest = prepared_identity.get("prepared", {}).get("files_sha256") + return { + "passed": passed, + "status_basis": "passed all verifier evidence criteria" if passed else "one or more verifier evidence criteria failed", + "criteria": criteria, + "mode_criteria": { + "mesh_identity": mesh_ok_by_mode, + "field_comparisons": comparison_ok_by_mode, + "stage_observability": observability_ok_by_mode, + }, + "entrypoint": report.get("harness_entrypoint", {}), + "regression_summary": { + "schema_version": report.get("harness", {}).get("schema_version"), + "source_case": report.get("problem_boundary", {}).get("airfrans_simulation"), + "prepared_case_files_sha256": prepared_digest, + "backend_selected": report.get("backend", {}).get("selected"), + "differentiability_status": report.get("differentiability", {}).get("status"), + "mesh_topology_sha256": oracle_mesh.get("topology_sha256"), + "mesh_geometry_sha256": oracle_mesh.get("geometry_sha256"), + "modes": regression_modes, + }, + } + + + def base_report(args: argparse.Namespace) -> dict[str, Any]: report_path = args.report if args.report is not None else args.work / "verifier_report.json" return { @@ -679,6 +1389,7 @@ def base_report(args: argparse.Namespace) -> dict[str, Any]: "name": "airfrans_stepper_verifier", "spec": "VERIFIER_HARNESS_SPEC.md", "schema_version": 1, + "scope": AIRFRANS_VERIFIER_SCOPE, }, "status": "running", "failure": None, @@ -688,17 +1399,89 @@ def base_report(args: argparse.Namespace) -> dict[str, Any]: "work": args.work, "report": report_path, }, + "harness_entrypoint": { + "command": "uv run python scripts/verify_airfrans_stepper.py --work --report ", + "arguments": { + "source": args.source, + "work": args.work, + "report": report_path, + "backend": args.backend, + "skip_split": args.skip_split, + "rtol": args.rtol, + "atol": args.atol, + }, + "prepares_case": True, + "runs_oracle": True, + "runs_migrated_path": True, + "writes_report": True, + "scope": AIRFRANS_VERIFIER_SCOPE, + }, "source_policy": "raw AirfRANS source is read-only; all solver runs use prepared work-directory copies", "tolerances": { "rtol": args.rtol, "atol": args.atol, "policy": "CPU stepper parity with the repository OpenFOAM oracle must pass np.allclose for every required field.", }, + "comparison_policy": { + "function": "numpy.allclose(actual, expected, rtol, atol, equal_nan=False)", + "scope": "internal arrays for U, p, phi, nut, k, and omega", + "backend_tolerances": { + "cpu": { + "rtol": args.rtol, + "atol": args.atol, + "reason": "strict CPU parity against the repository OpenFOAM oracle", + }, + "gpu": { + "rtol": args.rtol, + "atol": args.atol, + "reason": "full GPU RANS backend must own solver execution and satisfy OpenFOAM oracle parity; primitive GPU diagnostics are not acceptance", + "availability_checked_during_backend_selection": True, + "full_solver_guard": FULL_GPU_RANS_GUARD, + "primitive_evidence": { + "name": GPU_PRIMITIVE_NAME, + "path": GPU_PRIMITIVE_PROOF, + "counts_as_full_gpu_rans_solver": False, + }, + }, + }, + }, + "backend": { + "requested": args.backend, + "selected": None, + "failure_category": BACKEND_FAILURE, + "execution_scope": AIRFRANS_VERIFIER_SCOPE, + }, "required_fields": { "primary": list(PRIMARY_FIELDS), "turbulence": list(TURBULENCE_FIELDS), "all": list(REQUIRED_FIELDS), }, + "problem_boundary": { + "airfrans_simulation": args.source.name, + "selected_source_case": args.source, + "prepared_case_policy": "prepare one reproducible OpenFOAM v14 case and clone it into isolated oracle/run_one/split execution cases", + "oracle_output": { + "case_role": "oracle", + "producer": "foamRun -solver incompressibleFluid -noFunctionObjects", + "time": "1", + "fields": list(REQUIRED_FIELDS), + }, + "repository_outputs": { + "run_one": "migrated run_one path selected by the requested backend; gpu requests must use the GPU solver contract and reject CPU fallback", + "split": "split solver stages selected by the requested backend", + "time": "1", + "fields": list(REQUIRED_FIELDS), + }, + "mesh_identity_evidence": [ + "n_points", + "n_faces", + "n_internal_faces", + "n_cells", + "patches", + "topology_sha256", + "geometry_sha256", + ], + }, "commands": [], "case_preparation": {}, "mesh_identity": {}, @@ -708,6 +1491,26 @@ def base_report(args: argparse.Namespace) -> dict[str, Any]: }, "tracked_turbulence_fields": {}, "comparison_mismatches": [], + "comparison_attribution": [], + "verifier_evidence": { + "passed": False, + "status_basis": "not_evaluated", + }, + "timing": { + "status": "not_evaluated", + }, + "differentiability": { + "status": "not_evaluated", + "reason": "backend has not been selected yet", + }, + "state_exports": {}, + "gpu_solver_inputs": { + "schema_version": GPU_INPUT_SCHEMA_VERSION, + "status": "not_requested", + }, + "stage_observability": { + "contract": json_ready(STAGE_OBSERVABILITY_GROUPS), + }, } @@ -728,6 +1531,51 @@ def run_harness(args: argparse.Namespace, report: dict[str, Any]) -> None: split_case = prepared["split_case"] report["case_preparation"] = json_ready(prepared) + identity_mismatches = { + role: identity + for role, identity in prepared["prepared_case_identity"].items() + if not identity["matches_prepared"] + } + if identity_mismatches: + raise HarnessError( + LOADING_FAILURE, + "prepare_case_identity", + "prepared execution case copies differ from the reproducible baseline case", + details={"prepared_case_identity": identity_mismatches}, + ) + + + try: + backend = select_execution_backend(args.backend) + except BackendExecutionError as exc: + raise exc.to_harness_error() from exc + report["backend"] = json_ready(backend) + report["differentiability"] = json_ready(differentiability_report(backend, case_name=args.source.name)) + if backend.get("selected") == "gpu": + try: + foam = import_foam() + except Exception as exc: + raise HarnessError( + LOADING_FAILURE, + "import_foam_stepper", + "failed to import foam_stepper before GPU input preparation", + cause=exc, + ) from exc + gpu_input_stepper = make_stepper(foam, run_one_case, "gpu_inputs") + report["mesh_identity"]["gpu_inputs"] = read_mesh_identity(gpu_input_stepper, "gpu_inputs") + try: + gpu_inputs = prepare_gpu_solver_inputs(foam, gpu_input_stepper, run_one_case, prepared, backend) + except BackendExecutionError as exc: + raise exc.to_harness_error() from exc + report["gpu_solver_inputs"] = json_ready(gpu_inputs) + try: + gpu_stage_smoke = run_gpu_solver_stage_smoke(gpu_input_stepper, backend, run_one_case) + except BackendExecutionError as exc: + raise exc.to_harness_error() from exc + report["gpu_solver_stages"] = json_ready(gpu_solver_stage_report(gpu_stage_smoke)) + gpu_contract_blocker = gpu_solver_contract_blocker(backend, run_one_case) + if gpu_contract_blocker is not None: + record_gpu_contract_failure(report, gpu_contract_blocker, gpu_inputs) report["commands"].append( run_openfoam_command(["checkMesh", "-case", str(oracle_case), "-constant"], log_path=args.work / "checkMesh.log", timeout=120) ) @@ -769,6 +1617,9 @@ def run_harness(args: argparse.Namespace, report: dict[str, Any]) -> None: "fields": selected_field_summaries(oracle_fields), } report["tracked_turbulence_fields"]["oracle"] = visible_turbulence_fields(oracle_fields) + oracle_state = export_solver_state_checked(foam, oracle_stepper, "oracle", oracle_fields) + report["state_exports"]["oracle"] = foam.describe_solver_state(oracle_state) + run_one_stepper = make_stepper(foam, run_one_case, "run_one") run_one_mesh = read_mesh_identity(run_one_stepper, "run_one") @@ -777,34 +1628,59 @@ def run_harness(args: argparse.Namespace, report: dict[str, Any]) -> None: report["modes"]["run_one"]["mesh_comparison"] = mesh_report try: - run_one_result = run_one_stepper.run_one_pimple_iteration() - run_one_fields = field_dict_to_mapping(run_one_result.outputs["fields"]) + backend_run = run_backend_iteration(run_one_stepper, backend, run_one_case) + run_one_result = backend_run["result"] + run_one_fields = backend_run["fields"] + except BackendExecutionError as exc: + raise exc.to_harness_error() from exc except Exception as exc: raise HarnessError( STEPPER_FAILURE, "run_one_pimple_iteration", - "full one-iteration stepper execution failed", - details={"case": run_one_case}, + "full one-iteration backend execution failed", + details={"case": run_one_case, "backend": backend}, cause=exc, ) from exc + run_one_stages = graph_stage_summaries(run_one_result) + run_one_observability = stage_observability_report( + "run_one", + run_one_stages, + evidence={ + "one_iteration_execution": True, + "compared_fields": list(REQUIRED_FIELDS), + "turbulence_fields": visible_turbulence_fields(run_one_fields), + "backend": backend, + }, + ) + report["stage_observability"]["run_one"] = run_one_observability + run_one_comparisons, run_one_mismatches = compare_fields( "run_one", run_one_fields, oracle_fields, rtol=args.rtol, atol=args.atol, + stage_graph=[stage["name"] for stage in run_one_stages], ) report["modes"]["run_one"].update( { "case": run_one_case, + "execution_path": backend_run["execution_path"], + "backend": backend_run["backend"], "stage": transform_summary(run_one_result), - "graph": [entry.name for entry in run_one_result.outputs.get("graph", [])], + "graph": [stage["name"] for stage in run_one_stages], "fields": selected_field_summaries(run_one_fields), "comparisons": run_one_comparisons, + "observability": run_one_observability, } ) + if backend.get("selected") == "gpu": + report["modes"]["run_one"]["gpu_solver"] = json_ready(backend_run.get("gpu_solver", {})) report["tracked_turbulence_fields"]["run_one"] = visible_turbulence_fields(run_one_fields) + run_one_state = export_solver_state_checked(foam, run_one_stepper, "run_one", run_one_fields) + report["state_exports"]["run_one"] = foam.describe_solver_state(run_one_state) + mismatches = mesh_mismatches + run_one_mismatches @@ -815,7 +1691,11 @@ def run_harness(args: argparse.Namespace, report: dict[str, Any]) -> None: split_mesh_report, split_mesh_mismatches = compare_mesh_identity("split", split_mesh, oracle_mesh) report["modes"]["split"]["mesh_comparison"] = split_mesh_report - split_result = run_split_iteration(split_stepper) + if backend.get("selected") == "gpu": + split_result = run_gpu_split_iteration(foam, split_stepper, backend, split_case) + else: + split_result = run_split_iteration(foam, split_stepper) + report["stage_observability"]["split"] = split_result["observability"] split_fields = split_result["fields"] split_comparisons, split_mismatches = compare_fields( "split", @@ -823,6 +1703,7 @@ def run_harness(args: argparse.Namespace, report: dict[str, Any]) -> None: oracle_fields, rtol=args.rtol, atol=args.atol, + stage_graph=split_result["graph"], ) report["modes"]["split"].update( { @@ -834,19 +1715,33 @@ def run_harness(args: argparse.Namespace, report: dict[str, Any]) -> None: "pEqn": split_result["pEqn"], "fields": selected_field_summaries(split_fields), "comparisons": split_comparisons, + "observability": split_result["observability"], + "execution_path": split_result.get("execution_path"), + "backend": split_result.get("backend"), + "gpu_solver": split_result.get("gpu_solver"), } ) report["tracked_turbulence_fields"]["split"] = visible_turbulence_fields(split_fields) + report["state_exports"]["split"] = split_result["state_summary"] + report["state_exports"]["split_matrices"] = split_result["matrix_states"] mismatches.extend(split_mesh_mismatches) mismatches.extend(split_mismatches) + comparison_attribution = [mismatch["attribution"] for mismatch in mismatches if "attribution" in mismatch] + report["comparison_attribution"] = json_ready(comparison_attribution) report["comparison_mismatches"] = json_ready(mismatches) + report["timing"] = json_ready(build_timing_evidence(report)) + report["verifier_evidence"] = json_ready(build_verifier_evidence(report)) if mismatches: raise HarnessError( COMPARISON_FAILURE, "compare_oracle_stepper_outputs", f"{len(mismatches)} verifier comparison mismatch(es) observed", - details={"mismatches": mismatches}, + details={ + "failure_kind": GPU_NUMERICAL_MISMATCH if backend.get("selected") == "gpu" else COMPARISON_FAILURE, + "mismatches": mismatches, + "likely_causes": comparison_attribution, + }, ) @@ -878,18 +1773,58 @@ def print_human_summary(report: Mapping[str, Any]) -> None: print(f"source={report['paths']['source']}") print(f"work={report['paths']['work']}") print(f"report={report['paths']['report']}") + entrypoint = report.get("harness_entrypoint") or {} + print(f"harness_entrypoint={entrypoint.get('command')}") + print(f"harness_scope={entrypoint.get('scope') or report.get('harness', {}).get('scope')}") print(f"rtol={report['tolerances']['rtol']} atol={report['tolerances']['atol']}") + backend = report.get("backend", {}) + print(f"backend_requested={backend.get('requested')} backend_selected={backend.get('selected')}") + differentiability = report.get("differentiability") or {} + if differentiability: + claims = differentiability.get("claims") or [] + checks = differentiability.get("checks") or [] + print(f"differentiability_status={differentiability.get('status')}") + print(f"differentiability_claims={len(claims)} non_differentiable_boundaries={len(differentiability.get('non_differentiable') or [])}") + print(f"differentiability_checks={[check.get('status') for check in checks]}") + evidence = report.get("verifier_evidence") or {} + if evidence: + print(f"evidence_passed={evidence.get('passed')} evidence_status_basis={evidence.get('status_basis')}") + print(f"evidence_criteria={evidence.get('criteria')}") if status != "passed": failure = report.get("failure") or {} print(f"failure_category={failure.get('category')}") print(f"failure_step={failure.get('step')}") print(f"failure_message={failure.get('message')}") + timing = report.get("timing") or {} + if timing.get("status") == "reported": + openfoam = timing.get("openfoam_reference") or {} + gpu_solver = timing.get("gpu_solver") or {} + run_one_timing = ((gpu_solver.get("run_one") or {}).get("timing") or {}) + run_one_profiler = ((gpu_solver.get("run_one") or {}).get("profiler") or {}) + print(f"openfoam_oracle_duration_seconds={openfoam.get('oracle_solver_duration_seconds')}") + print(f"gpu_solver_run_one_wall_seconds={run_one_timing.get('gpu_solver_wall_seconds')}") + print(f"gpu_solver_run_one_device_time_ms={run_one_profiler.get('device_time_ms_total')}") + mismatches = (failure.get("details") or {}).get("mismatches") or report.get("comparison_mismatches") or [] + if mismatches: + first = mismatches[0] + attribution = first.get("attribution") or {} + print(f"comparison_failure_mode={first.get('mode')}") + print(f"comparison_failure_field={first.get('field')}") + print(f"comparison_failure_reason={first.get('reason')}") + print(f"comparison_failure_likely_stage_group={attribution.get('likely_stage_group')}") + print(f"comparison_failure_likely_stages={attribution.get('likely_stages')}") + print(f"comparison_failure_max_abs={fmt_sci(first.get('max_abs'))}") + print(f"comparison_failure_location={format_location(first.get('location'))}") return oracle_mesh = report["mesh_identity"]["oracle"] patch_names = [patch["name"] for patch in oracle_mesh["patches"]] print(f"oracle_case={report['case_preparation']['oracle_case']}") + prepared_identity = report["case_preparation"].get("prepared_case_identity", {}) + prepared_digest = prepared_identity.get("prepared", {}).get("files_sha256") + print(f"prepared_case={report['case_preparation']['prepared_case']}") + print(f"prepared_case_files_sha256={prepared_digest}") print(f"run_one_case={report['modes']['run_one']['case']}") if report["modes"].get("split", {}).get("enabled"): print(f"split_case={report['modes']['split']['case']}") @@ -914,6 +1849,10 @@ def print_human_summary(report: Mapping[str, Any]) -> None: graph = mode_report.get("graph") if graph: print(f"{mode}_graph={graph}") + observability = mode_report.get("observability") + if observability: + observed_groups = [group["name"] for group in observability.get("groups", []) if group.get("observed")] + print(f"{mode}_observed_stage_groups={observed_groups}") split_report = report["modes"].get("split", {}) if split_report.get("enabled") and "UEqn" in split_report and "pEqn" in split_report: @@ -931,6 +1870,7 @@ def parse_args(argv: list[str] | None = None) -> argparse.Namespace: parser.add_argument("--report", type=Path, default=None, help="JSON report path; defaults to WORK/verifier_report.json") parser.add_argument("--rtol", type=float, default=1e-8) parser.add_argument("--atol", type=float, default=1e-8) + parser.add_argument("--backend", choices=BACKEND_CHOICES, default="auto", help="Backend request; gpu means the full RANS solver backend and currently fails until that backend exists") parser.add_argument("--skip-split", action="store_true", help="Skip the explicit split-step parity check for local debugging") return parser.parse_args(argv) @@ -948,6 +1888,14 @@ def main(argv: list[str] | None = None) -> int: exit_code = 0 try: run_harness(args, report) + evidence = report.get("verifier_evidence") + if not isinstance(evidence, Mapping) or evidence.get("passed") is not True: + raise HarnessError( + INTERNAL_FAILURE, + "verifier_evidence", + "verifier completed without passing evidence criteria", + details={"verifier_evidence": evidence}, + ) report["status"] = "passed" except HarnessError as exc: report["status"] = "failed" diff --git a/scripts/verify_gpu_algorithm.py b/scripts/verify_gpu_algorithm.py new file mode 100755 index 0000000..bc67053 --- /dev/null +++ b/scripts/verify_gpu_algorithm.py @@ -0,0 +1,445 @@ +#!/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()) diff --git a/scripts/verify_gpu_algorithm.sh b/scripts/verify_gpu_algorithm.sh new file mode 100755 index 0000000..f6405a6 --- /dev/null +++ b/scripts/verify_gpu_algorithm.sh @@ -0,0 +1,36 @@ +#!/usr/bin/env bash +set -euo pipefail + +ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd -P)" +JOBS="${JOBS:-$(nproc)}" +PYTHON_BIN="${PYTHON_BIN:-${ROOT_DIR}/.venv/bin/python}" + +if [[ ! -x "${PYTHON_BIN}" ]]; then + uv sync --dev +fi + +if ! PYTHONPATH= "${PYTHON_BIN}" -c 'import pybind11, numpy, quadrants' >/dev/null 2>&1; then + uv sync --dev +fi + +"${ROOT_DIR}/scripts/build_openfoam_airfrans_subset.sh" >/dev/null + +set +u +source "${ROOT_DIR}/OpenFOAM-14/etc/bashrc" \ + WM_MPLIB=Dummy \ + ParaView_TYPE=none \ + SCOTCH_TYPE=none \ + ZOLTAN_TYPE=none +set -u +unset FOAM_SIGFPE + +( + cd "${ROOT_DIR}/OpenFOAM-14" + wmake -j "${JOBS}" libso src/mesh/blockMesh >/dev/null + wmake -j "${JOBS}" applications/utilities/mesh/generation/blockMesh >/dev/null + wmake -j "${JOBS}" applications/utilities/mesh/manipulation/createZones >/dev/null +) + +"${ROOT_DIR}/scripts/build_python_stepper.sh" >/dev/null + +PYTHONPATH= "${PYTHON_BIN}" "${ROOT_DIR}/scripts/verify_gpu_algorithm.py" "$@" diff --git a/scripts/verify_gpu_rans_solver.sh b/scripts/verify_gpu_rans_solver.sh new file mode 100755 index 0000000..c70f576 --- /dev/null +++ b/scripts/verify_gpu_rans_solver.sh @@ -0,0 +1,136 @@ +#!/usr/bin/env bash +set -euo pipefail + +ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd -P)" +JOBS="${JOBS:-$(nproc)}" +PYTHON_BIN="${PYTHON_BIN:-${ROOT_DIR}/.venv/bin/python}" +REPORT="" +PASSTHRU=() + +while [[ $# -gt 0 ]]; do + case "$1" in + --report) + if [[ $# -lt 2 ]]; then + echo "--report requires a path" >&2 + exit 2 + fi + REPORT="$2" + PASSTHRU+=("$1" "$2") + shift 2 + ;; + --report=*) + REPORT="${1#--report=}" + PASSTHRU+=("$1") + shift + ;; + --backend|--backend=*) + echo "verify_gpu_rans_solver.sh always requests --backend gpu; do not pass --backend" >&2 + exit 2 + ;; + *) + PASSTHRU+=("$1") + shift + ;; + esac +done + +if [[ -z "${REPORT}" ]]; then + echo "--report is required so the full GPU RANS guard can inspect verifier evidence" >&2 + exit 2 +fi + +if [[ ! -x "${PYTHON_BIN}" ]]; then + uv sync --dev +fi + +if ! PYTHONPATH= "${PYTHON_BIN}" -c 'import pybind11, numpy, quadrants' >/dev/null 2>&1; then + uv sync --dev +fi + +"${ROOT_DIR}/scripts/build_openfoam_airfrans_subset.sh" >/dev/null + +set +u +source "${ROOT_DIR}/OpenFOAM-14/etc/bashrc" \ + WM_MPLIB=Dummy \ + ParaView_TYPE=none \ + SCOTCH_TYPE=none \ + ZOLTAN_TYPE=none +set -u +unset FOAM_SIGFPE + +"${ROOT_DIR}/scripts/build_python_stepper.sh" >/dev/null + +set +e +PYTHONPATH= "${PYTHON_BIN}" "${ROOT_DIR}/scripts/verify_airfrans_stepper.py" --backend gpu "${PASSTHRU[@]}" +VERIFY_STATUS=$? +set -e + +PYTHONPATH= "${PYTHON_BIN}" - "${REPORT}" "${VERIFY_STATUS}" <<'PY' +from __future__ import annotations + +import json +import sys +from pathlib import Path + +report_path = Path(sys.argv[1]) +verify_status = int(sys.argv[2]) +required_fields = ("U", "p", "phi", "nut", "k", "omega") +required_modes = ("run_one", "split") +failures: list[str] = [] + +if not report_path.exists(): + failures.append(f"report missing: {report_path}") + report = {} +else: + try: + report = json.loads(report_path.read_text()) + except Exception as exc: # pragma: no cover - shell guard diagnostic + failures.append(f"report is not valid JSON: {exc}") + report = {} + + +def require(condition: bool, message: str) -> None: + if not condition: + failures.append(message) + + +backend = report.get("backend") if isinstance(report.get("backend"), dict) else {} +verifier_evidence = report.get("verifier_evidence") if isinstance(report.get("verifier_evidence"), dict) else {} +modes = report.get("modes") if isinstance(report.get("modes"), dict) else {} +provider = str(backend.get("provider") or "").lower() +device = str(backend.get("device") or "").lower() + +require(verify_status == 0, f"underlying AirfRANS verifier exited {verify_status}") +require(report.get("status") == "passed", f"report.status is {report.get('status')!r}, expected 'passed'") +require(backend.get("requested") == "gpu", f"backend.requested is {backend.get('requested')!r}, expected 'gpu'") +require(backend.get("selected") == "gpu", f"backend.selected is {backend.get('selected')!r}, expected 'gpu'") +require("cpu" not in provider and provider not in {"foam_stepper_cpu", "openfoam"}, f"backend.provider looks CPU-backed: {backend.get('provider')!r}") +require(device not in {"host", "cpu"}, f"backend.device looks CPU-backed: {backend.get('device')!r}") +require(backend.get("counts_as_gpu_algorithm_progress") is not False, "backend explicitly says it does not count as GPU progress") +require(verifier_evidence.get("passed") is True, "verifier_evidence.passed is not true") + +for mode_name in required_modes: + mode = modes.get(mode_name) if isinstance(modes.get(mode_name), dict) else {} + require(mode.get("enabled") is True, f"modes.{mode_name}.enabled is not true") + mode_backend = mode.get("backend") if isinstance(mode.get("backend"), dict) else backend + require(mode_backend.get("selected") == "gpu", f"modes.{mode_name}.backend.selected is not 'gpu'") + execution_path = str(mode.get("execution_path") or "").lower() + require("gpu" in execution_path, f"modes.{mode_name}.execution_path does not identify a GPU path: {mode.get('execution_path')!r}") + comparisons = mode.get("comparisons") if isinstance(mode.get("comparisons"), dict) else {} + for field in required_fields: + comparison = comparisons.get(field) if isinstance(comparisons.get(field), dict) else {} + require(comparison.get("allclose") is True, f"modes.{mode_name}.comparisons.{field}.allclose is not true") + require("actual_shape" in comparison or "shape" in comparison, f"modes.{mode_name}.comparisons.{field} has no actual shape evidence") + require("expected_shape" in comparison or "shape" in comparison, f"modes.{mode_name}.comparisons.{field} has no expected shape evidence") + require("max_abs" in comparison or "max_abs_error" in comparison, f"modes.{mode_name}.comparisons.{field} has no max error evidence") + +if failures: + print("full GPU RANS solver guard failed:", file=sys.stderr) + for failure in failures: + print(f"- {failure}", file=sys.stderr) + print(f"report={report_path}", file=sys.stderr) + sys.exit(1) + +print("full GPU RANS solver guard passed") +print(f"report={report_path}") +PY diff --git a/scripts/verify_gpu_step_timing.py b/scripts/verify_gpu_step_timing.py new file mode 100755 index 0000000..ed0be6c --- /dev/null +++ b/scripts/verify_gpu_step_timing.py @@ -0,0 +1,333 @@ +#!/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()) diff --git a/scripts/verify_gpu_step_timing.sh b/scripts/verify_gpu_step_timing.sh new file mode 100755 index 0000000..4ecb950 --- /dev/null +++ b/scripts/verify_gpu_step_timing.sh @@ -0,0 +1,36 @@ +#!/usr/bin/env bash +set -euo pipefail + +ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd -P)" +JOBS="${JOBS:-$(nproc)}" +PYTHON_BIN="${PYTHON_BIN:-${ROOT_DIR}/.venv/bin/python}" + +if [[ ! -x "${PYTHON_BIN}" ]]; then + uv sync --dev +fi + +if ! PYTHONPATH= "${PYTHON_BIN}" -c 'import pybind11, numpy, quadrants' >/dev/null 2>&1; then + uv sync --dev +fi + +"${ROOT_DIR}/scripts/build_openfoam_airfrans_subset.sh" >/dev/null + +set +u +source "${ROOT_DIR}/OpenFOAM-14/etc/bashrc" \ + WM_MPLIB=Dummy \ + ParaView_TYPE=none \ + SCOTCH_TYPE=none \ + ZOLTAN_TYPE=none +set -u +unset FOAM_SIGFPE + +( + cd "${ROOT_DIR}/OpenFOAM-14" + wmake -j "${JOBS}" libso src/mesh/blockMesh >/dev/null + wmake -j "${JOBS}" applications/utilities/mesh/generation/blockMesh >/dev/null + wmake -j "${JOBS}" applications/utilities/mesh/manipulation/createZones >/dev/null +) + +"${ROOT_DIR}/scripts/build_python_stepper.sh" >/dev/null + +PYTHONPATH= "${PYTHON_BIN}" "${ROOT_DIR}/scripts/verify_gpu_step_timing.py" "$@"