From b40d981f79c20b683c9ac3cff0b473ffe6c68ebc Mon Sep 17 00:00:00 2001 From: Zachery Aaron Shores-Chmielewski Date: Mon, 27 Jul 2026 12:35:44 +0400 Subject: [PATCH] feat: some progress on gpu implementation, more granularity on solver numerical comparison --- .gitignore | 8 + python/src/foam_stepper/gpu/backend.py | 779 +++++++++-- python/src/foam_stepper/gpu/constants.py | 32 +- python/src/foam_stepper/gpu/kernels.py | 432 +++++- python/src/foam_stepper/gpu/linear_solve.py | 216 ++- scripts/update_loop_diagnostic_context.py | 430 ++++++ scripts/verify_airfrans_stepper.py | 1340 ++++++++++++++++++- scripts/verify_gpu_rans_solver.sh | 16 +- 8 files changed, 3108 insertions(+), 145 deletions(-) create mode 100755 scripts/update_loop_diagnostic_context.py diff --git a/.gitignore b/.gitignore index 83b43c8..b25b9ba 100644 --- a/.gitignore +++ b/.gitignore @@ -4,6 +4,14 @@ python/src/foam_stepper.egg-info/ tmp/ +# local planning/spec context injected by the loop harness +CURRENT_STATUS.md +MATCHING_AIRFRANS_SIMULATION.md +PROJECT_GROUNDING.md +PYTHON_STEPPER_SPEC.md +STAGE_ORACLE_SPEC.md +VERIFIER_HARNESS_SPEC.md + # foreign OpenFOAM-14/ ThirdParty-14/ diff --git a/python/src/foam_stepper/gpu/backend.py b/python/src/foam_stepper/gpu/backend.py index 7ff4683..fd6f3e3 100644 --- a/python/src/foam_stepper/gpu/backend.py +++ b/python/src/foam_stepper/gpu/backend.py @@ -8,6 +8,7 @@ backend diagnostics used by ``--backend gpu``. from __future__ import annotations import dataclasses +import hashlib import math import os import subprocess @@ -24,6 +25,7 @@ from .constants import ( DEFAULT_MOMENTUM_PBICGSTAB_ITERATIONS, DEFAULT_MOMENTUM_PBICGSTAB_RESIDUAL_TOLERANCE_SQUARED, DEFAULT_PRESSURE_CG_ITERATIONS, + DEFAULT_PRESSURE_NON_ORTHOGONAL_CORRECTORS, DEFAULT_MOMENTUM_RELAXATION_ALPHA, DEFAULT_SIMPLE_CONSISTENT_RATU_FACTOR, DEFAULT_OMEGA_WALL_BETA1, @@ -49,27 +51,55 @@ from .kernels import ( gpu_rans_momentum_diffusion_coefficients, gpu_rans_momentum_wall_diffusion_coefficients, gpu_rans_momentum_convection_coefficients, + gpu_rans_momentum_bounded_convection_sp_internal, + gpu_rans_momentum_bounded_convection_sp_boundary, + gpu_rans_momentum_convection_boundary_coefficients, + gpu_rans_momentum_convection_boundary_source, + gpu_rans_momentum_boundary_internal_diag, + gpu_rans_momentum_boundary_relaxation_coefficients, + gpu_rans_add_scalar_field, + gpu_rans_zero_tensor_field, + gpu_rans_momentum_gauss_grad_u_internal, + gpu_rans_momentum_gauss_grad_u_boundary, + gpu_rans_momentum_gauss_grad_u_finish, + gpu_rans_momentum_linear_upwind_source, + gpu_rans_zero_scalar_field, + gpu_rans_momentum_offdiag_abs_accumulate, gpu_rans_momentum_equation_relaxation, + gpu_rans_momentum_pressure_gradient_source, + gpu_rans_momentum_pressure_boundary_source, + gpu_rans_momentum_h1_face_accumulate, + gpu_rans_momentum_h1_finish, 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_consistent_phiHbyA_correction, 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_negate_scalar_field, gpu_rans_surface_flux_from_cells, gpu_rans_turbulence_update, gpu_rans_omega_wall_update, gpu_rans_pressure_velocity_correction, ) from .linear_solve import ( + build_losort_addr, + gpu_bicgstab_initialize_vector, + gpu_bicgstab_dot_vector, + gpu_bicgstab_precondition_vector, + gpu_copy_scalar, + gpu_copy_vector, + gpu_dilu_apply_vector_asymmetric_faces, + gpu_ldu_matvec_vector_asymmetric_faces, gpu_ldu_pbicgstab_vector_asymmetric_faces, gpu_ldu_pcg_scalar_symmetric_faces, + gpu_zero_scalar_accumulator, ) BACKEND_CHOICES = ("auto", "cpu", "gpu") @@ -183,6 +213,9 @@ def array_stats(array: Any) -> dict[str, Any]: "dtype": str(arr.dtype), "size": int(arr.size), } + if arr.dtype != object: + contiguous = np.ascontiguousarray(arr) + out["sha256"] = hashlib.sha256(contiguous.tobytes()).hexdigest() if arr.size == 0 or not np.issubdtype(arr.dtype, np.number): return out @@ -221,6 +254,63 @@ def entity_value_at(array: np.ndarray, entity_index: int | None) -> Any: return json_ready(value) +def write_numeric_artifact_file(path: Path, arrays: Mapping[str, np.ndarray], *, role: str) -> dict[str, Any]: + materialized = {name: np.ascontiguousarray(value) for name, value in arrays.items()} + path.parent.mkdir(parents=True, exist_ok=True) + np.savez_compressed(path, **materialized) + return { + "schema_version": 1, + "role": role, + "path": str(path), + "format": "npz", + "artifact_count": len(materialized), + "artifacts": {name: array_stats(value) for name, value in materialized.items()}, + } + + +def array_difference_summary(actual: np.ndarray, reference: np.ndarray) -> dict[str, Any]: + actual_arr = np.asarray(actual) + reference_arr = np.asarray(reference) + out: dict[str, Any] = { + "actual_shape": array_shape(actual_arr), + "reference_shape": array_shape(reference_arr), + "shape_matches": actual_arr.shape == reference_arr.shape, + "actual_dtype": str(actual_arr.dtype), + "reference_dtype": str(reference_arr.dtype), + } + if actual_arr.shape != reference_arr.shape: + return {**out, "reason": "shape_mismatch"} + if actual_arr.size == 0: + return {**out, "max_abs": 0.0, "mean_abs": 0.0, "rms_abs": 0.0, "largest_difference": None} + diff = np.abs(actual_arr - reference_arr) + finite = np.isfinite(diff) + finite_diff = diff[finite] + if finite_diff.size: + flat_index = int(np.argmax(np.where(finite, diff, -np.inf))) + max_abs = finite_float(diff.reshape(-1)[flat_index]) + mean_abs = finite_float(np.mean(finite_diff)) + rms_abs = finite_float(np.sqrt(np.mean(np.square(finite_diff)))) + else: + flat_index = int(np.flatnonzero(~finite.reshape(-1))[0]) + max_abs = None + mean_abs = None + rms_abs = None + index = tuple(int(item) for item in np.unravel_index(flat_index, diff.shape)) + return { + **out, + "reason": None, + "max_abs": max_abs, + "mean_abs": mean_abs, + "rms_abs": rms_abs, + "nonfinite_error_count": int(diff.size - np.count_nonzero(finite)), + "largest_difference": { + "array_index": list(index), + "actual_value": value_at(actual_arr, index), + "reference_value": value_at(reference_arr, index), + }, + } + + def update_hash_text(digest: "hashlib._Hash", value: str) -> None: encoded = value.encode("utf-8") digest.update(len(encoded).to_bytes(8, "little")) @@ -873,6 +963,67 @@ def gpu_i32_array(source: Any, path: str) -> tuple[Any, np.ndarray, dict[str, An "transferred": True, } +def openfoam_dilu_preconditioner_diag( + diag: np.ndarray, + owner: np.ndarray, + neighbour: np.ndarray, + upper: np.ndarray, + lower: np.ndarray, +) -> np.ndarray: + """Return a diagonal equivalent that applies OpenFOAM's DILU reciprocal diagonal.""" + preconditioned_diag = np.asarray(diag, dtype=np.float64).copy() + owner_cells = np.asarray(owner, dtype=np.int32).reshape(-1) + neighbour_cells = np.asarray(neighbour, dtype=np.int32).reshape(-1) + upper_coeffs = np.asarray(upper, dtype=np.float64).reshape(-1) + lower_coeffs = np.asarray(lower, dtype=np.float64).reshape(-1) + for face in range(upper_coeffs.shape[0]): + owner_cell = int(owner_cells[face]) + neighbour_cell = int(neighbour_cells[face]) + owner_diag = preconditioned_diag[owner_cell] + if abs(owner_diag) > 1.0e-300: + preconditioned_diag[neighbour_cell] -= upper_coeffs[face] * lower_coeffs[face] / owner_diag + safe = np.where(np.abs(preconditioned_diag) > 1.0e-300, preconditioned_diag, 1.0) + return np.ascontiguousarray(safe, dtype=np.float64) + + +def openfoam_dilu_apply_vector_reference( + owner: np.ndarray, + neighbour: np.ndarray, + losort: np.ndarray, + upper: np.ndarray, + lower: np.ndarray, + reciprocal_diag: np.ndarray, + source: np.ndarray, +) -> np.ndarray: + """Apply the same forward/back DILU sweep on CPU for diagnostic comparison.""" + owner_cells = np.asarray(owner, dtype=np.int32).reshape(-1) + neighbour_cells = np.asarray(neighbour, dtype=np.int32).reshape(-1) + losort_addr = np.asarray(losort, dtype=np.int32).reshape(-1) + upper_coeffs = np.asarray(upper, dtype=np.float64).reshape(-1) + lower_coeffs = np.asarray(lower, dtype=np.float64).reshape(-1) + reciprocal = np.asarray(reciprocal_diag, dtype=np.float64).reshape(-1) + out = np.asarray(source, dtype=np.float64).copy() + out *= reciprocal[:, None] + for face_index in losort_addr: + face = int(face_index) + owner_cell = int(owner_cells[face]) + neighbour_cell = int(neighbour_cells[face]) + out[neighbour_cell] -= reciprocal[neighbour_cell] * lower_coeffs[face] * out[owner_cell] + for face in range(upper_coeffs.shape[0] - 1, -1, -1): + owner_cell = int(owner_cells[face]) + neighbour_cell = int(neighbour_cells[face]) + out[owner_cell] -= reciprocal[owner_cell] * upper_coeffs[face] * out[neighbour_cell] + return np.ascontiguousarray(out, dtype=np.float64) + +def openfoam_dic_preconditioner_diag( + diag: np.ndarray, + owner: np.ndarray, + neighbour: np.ndarray, + upper: np.ndarray, +) -> np.ndarray: + """Return a diagonal equivalent that applies OpenFOAM's DIC reciprocal diagonal.""" + return openfoam_dilu_preconditioner_diag(diag, owner, neighbour, upper, upper) + def gpu_empty_f64(shape: tuple[int, ...], path: str) -> tuple[Any, dict[str, Any]]: gpu_array = qd.ndarray(qd.f64, shape=shape) @@ -1017,7 +1168,7 @@ def quadrants_kernel_evidence(expected: Iterable[str]) -> dict[str, Any]: } -def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: Path) -> dict[str, Any]: +def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: Path, *, diagnostic_artifact_path: Path | None = None) -> dict[str, Any]: fields = read_fields(stepper, "gpu_stage_smoke") gpu_solver_started = time.perf_counter() input_transfer_started = time.perf_counter() @@ -1033,33 +1184,222 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P 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") + losort_np = build_losort_addr(n_cells, neighbour_np) + losort_gpu, losort_np, losort_transfer = gpu_i32_array(losort_np, "mesh.connectivity.losort.internal") sf_gpu, sf_np, sf_transfer = gpu_f64_array(np.asarray(mesh.Sf)[:n_internal_faces], "mesh.geometry.Sf.internal") + face_centres_gpu, face_centres_np, face_centres_transfer = gpu_f64_array(np.asarray(mesh.Cf)[:n_internal_faces], "mesh.geometry.Cf.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] = [] + boundary_phiHbyA_parts: list[np.ndarray] = [] for patch in mesh.boundary: patch_phi = fields["phi"].boundary.get(patch.name) - if patch_phi is None: + phi_values = np.asarray(patch_phi.values, dtype=np.float64).reshape(-1) if patch_phi is not None else None + u_patch = fields["U"].boundary.get(patch.name) + p_patch = fields["p"].boundary.get(patch.name) + use_constrained_hbyA_boundary = ( + u_patch is not None + and not getattr(u_patch, "assignable", True) + and getattr(p_patch, "type", "") != "fixedFluxExtrapolatedPressure" + ) + if use_constrained_hbyA_boundary: + u_values = np.asarray(u_patch.values, dtype=np.float64).reshape(-1, 3) + sf_values = np.asarray(patch.Sf, dtype=np.float64).reshape(-1, 3) + values = np.sum(u_values * sf_values, axis=1) + elif phi_values is not None: + values = phi_values + else: 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 phiHbyA and boundary face-cell arrays have different lengths", + details={"patch": patch.name, "phiHbyA_size": int(values.shape[0]), "face_cell_size": int(face_cells.shape[0])}, + ) + if phi_values is not None and face_cells.shape[0] != phi_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])}, + details={"patch": patch.name, "phi_size": int(phi_values.shape[0]), "face_cell_size": int(face_cells.shape[0])}, ) boundary_face_cells_parts.append(face_cells) - boundary_phi_parts.append(values) + boundary_phi_parts.append(phi_values if phi_values is not None else values) + boundary_phiHbyA_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_phiHbyA_np = np.concatenate(boundary_phiHbyA_parts) if boundary_phiHbyA_parts else np.empty((0,), dtype=np.float64) + n_pressure_boundary_faces = int(boundary_phiHbyA_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") + boundary_phiHbyA_gpu, boundary_phiHbyA_np, boundary_phiHbyA_transfer = gpu_f64_array(boundary_phiHbyA_np, "fields.phiHbyA.boundary.pressure_source") + momentum_pressure_boundary_face_cells_parts: list[np.ndarray] = [] + momentum_pressure_boundary_values_parts: list[np.ndarray] = [] + momentum_pressure_boundary_sf_parts: list[np.ndarray] = [] + for patch in mesh.boundary: + p_patch = fields["p"].boundary.get(patch.name) + if p_patch is None: + continue + values = np.asarray(p_patch.values, dtype=np.float64).reshape(-1) + if values.size == 0: + continue + face_cells = np.asarray(patch.face_cells, dtype=np.int32).reshape(-1) + sf_values = np.asarray(patch.Sf, dtype=np.float64).reshape(-1, 3) + if face_cells.shape[0] != values.shape[0] or sf_values.shape != (values.shape[0], 3): + raise BackendExecutionError( + "prepare_momentum_pressure_boundary_source", + "boundary pressure, face-cell, and area-vector arrays have incompatible shapes", + details={ + "patch": patch.name, + "p_size": int(values.shape[0]), + "face_cell_size": int(face_cells.shape[0]), + "Sf_shape": list(sf_values.shape), + }, + ) + momentum_pressure_boundary_face_cells_parts.append(face_cells) + momentum_pressure_boundary_values_parts.append(values) + momentum_pressure_boundary_sf_parts.append(sf_values) + momentum_pressure_boundary_face_cells_np = np.concatenate(momentum_pressure_boundary_face_cells_parts) if momentum_pressure_boundary_face_cells_parts else np.empty((0,), dtype=np.int32) + momentum_pressure_boundary_values_np = np.concatenate(momentum_pressure_boundary_values_parts) if momentum_pressure_boundary_values_parts else np.empty((0,), dtype=np.float64) + momentum_pressure_boundary_sf_np = np.concatenate(momentum_pressure_boundary_sf_parts) if momentum_pressure_boundary_sf_parts else np.empty((0, 3), dtype=np.float64) + n_momentum_pressure_boundary_faces = int(momentum_pressure_boundary_values_np.shape[0]) + momentum_pressure_boundary_face_cells_gpu, momentum_pressure_boundary_face_cells_np, momentum_pressure_boundary_face_cells_transfer = gpu_i32_array(momentum_pressure_boundary_face_cells_np, "mesh.boundary.face_cells.momentum_pressure_source") + momentum_pressure_boundary_values_gpu, momentum_pressure_boundary_values_np, momentum_pressure_boundary_values_transfer = gpu_f64_array(momentum_pressure_boundary_values_np, "fields.p.boundary.momentum_pressure_source") + momentum_pressure_boundary_sf_gpu, momentum_pressure_boundary_sf_np, momentum_pressure_boundary_sf_transfer = gpu_f64_array(momentum_pressure_boundary_sf_np, "mesh.boundary.Sf.momentum_pressure_source") + momentum_u_boundary_face_cells_parts: list[np.ndarray] = [] + momentum_u_boundary_values_parts: list[np.ndarray] = [] + momentum_u_boundary_sf_parts: list[np.ndarray] = [] + for patch in mesh.boundary: + u_patch = fields["U"].boundary.get(patch.name) + if u_patch is None: + continue + values = np.asarray(u_patch.values, dtype=np.float64).reshape(-1, 3) + if values.size == 0: + continue + face_cells = np.asarray(patch.face_cells, dtype=np.int32).reshape(-1) + sf_values = np.asarray(patch.Sf, dtype=np.float64).reshape(-1, 3) + if values.shape != (face_cells.shape[0], 3) or sf_values.shape != (face_cells.shape[0], 3): + raise BackendExecutionError( + "prepare_momentum_gauss_grad_u_boundary", + "boundary velocity, face-cell, and area-vector arrays have incompatible shapes", + details={ + "patch": patch.name, + "U_shape": list(values.shape), + "face_cell_size": int(face_cells.shape[0]), + "Sf_shape": list(sf_values.shape), + }, + ) + momentum_u_boundary_face_cells_parts.append(face_cells) + momentum_u_boundary_values_parts.append(values) + momentum_u_boundary_sf_parts.append(sf_values) + momentum_u_boundary_face_cells_np = np.concatenate(momentum_u_boundary_face_cells_parts) if momentum_u_boundary_face_cells_parts else np.empty((0,), dtype=np.int32) + momentum_u_boundary_values_np = np.concatenate(momentum_u_boundary_values_parts) if momentum_u_boundary_values_parts else np.empty((0, 3), dtype=np.float64) + momentum_u_boundary_sf_np = np.concatenate(momentum_u_boundary_sf_parts) if momentum_u_boundary_sf_parts else np.empty((0, 3), dtype=np.float64) + n_momentum_u_boundary_faces = int(momentum_u_boundary_face_cells_np.shape[0]) + momentum_u_boundary_face_cells_gpu, momentum_u_boundary_face_cells_np, momentum_u_boundary_face_cells_transfer = gpu_i32_array(momentum_u_boundary_face_cells_np, "mesh.boundary.face_cells.momentum_grad_u") + momentum_u_boundary_values_gpu, momentum_u_boundary_values_np, momentum_u_boundary_values_transfer = gpu_f64_array(momentum_u_boundary_values_np, "fields.U.boundary.momentum_grad_u") + momentum_u_boundary_sf_gpu, momentum_u_boundary_sf_np, momentum_u_boundary_sf_transfer = gpu_f64_array(momentum_u_boundary_sf_np, "mesh.boundary.Sf.momentum_grad_u") + momentum_u_convection_face_cells_parts: list[np.ndarray] = [] + momentum_u_convection_phi_parts: list[np.ndarray] = [] + momentum_u_convection_internal_coeff_parts: list[np.ndarray] = [] + momentum_u_convection_boundary_coeff_parts: list[np.ndarray] = [] + momentum_u_convection_boundary_source_face_cells_parts: list[np.ndarray] = [] + momentum_u_convection_boundary_source_parts: list[np.ndarray] = [] + momentum_u_convection_boundary_diag_face_cells_parts: list[np.ndarray] = [] + momentum_u_convection_boundary_diag_coeff_parts: list[np.ndarray] = [] + for patch in mesh.boundary: + u_patch = fields["U"].boundary.get(patch.name) + phi_patch = fields["phi"].boundary.get(patch.name) + if u_patch is None or phi_patch is None: + continue + phi_values = np.asarray(phi_patch.values, dtype=np.float64).reshape(-1) + if phi_values.size == 0: + continue + face_cells = np.asarray(patch.face_cells, dtype=np.int32).reshape(-1) + if face_cells.shape[0] != phi_values.shape[0]: + raise BackendExecutionError( + "prepare_momentum_convection_boundary_coefficients", + "boundary phi and face-cell arrays have different lengths", + details={"patch": patch.name, "phi_size": int(phi_values.shape[0]), "face_cell_size": int(face_cells.shape[0])}, + ) + u_values = np.asarray(u_patch.values, dtype=np.float64).reshape(-1, 3) + if u_values.shape != (face_cells.shape[0], 3): + raise BackendExecutionError( + "prepare_momentum_convection_boundary_coefficients", + "boundary U values and face-cell arrays have incompatible shapes", + details={"patch": patch.name, "U_shape": list(u_values.shape), "face_cell_size": int(face_cells.shape[0])}, + ) + def vector_coeff(raw: np.ndarray | None, name: str) -> np.ndarray: + if raw is None: + return np.zeros((face_cells.shape[0], 3), dtype=np.float64) + coeff = np.asarray(raw, dtype=np.float64) + if coeff.shape == (face_cells.shape[0],): + return np.repeat(coeff[:, None], 3, axis=1) + if coeff.shape == (face_cells.shape[0], 3): + return coeff + if coeff.shape == (face_cells.shape[0], 3, 3): + return np.stack((coeff[:, 0, 0], coeff[:, 1, 1], coeff[:, 2, 2]), axis=1) + raise BackendExecutionError( + "prepare_momentum_convection_boundary_coefficients", + "unsupported U boundary value coefficient shape", + details={"patch": patch.name, "coefficient": name, "shape": list(coeff.shape), "face_count": int(face_cells.shape[0])}, + ) + internal_coeffs = vector_coeff(getattr(u_patch, "value_internal_coeffs", None), "value_internal_coeffs") + boundary_coeffs = vector_coeff(getattr(u_patch, "value_boundary_coeffs", None), "value_boundary_coeffs") + if not np.any(internal_coeffs) and not np.any(boundary_coeffs) and getattr(u_patch, "type", "") == "freestreamVelocity": + patch_internal = np.asarray(u_patch.patch_internal if u_patch.patch_internal is not None else u_values, dtype=np.float64).reshape(-1, 3) + face_area_vectors = np.asarray(patch.Sf, dtype=np.float64).reshape(-1, 3) + if patch_internal.shape != (face_cells.shape[0], 3) or face_area_vectors.shape != (face_cells.shape[0], 3): + raise BackendExecutionError( + "prepare_momentum_convection_boundary_coefficients", + "freestreamVelocity patch arrays have incompatible shapes", + details={ + "patch": patch.name, + "patch_internal_shape": list(patch_internal.shape), + "face_area_vectors_shape": list(face_area_vectors.shape), + "face_cell_size": int(face_cells.shape[0]), + }, + ) + face_area_magnitudes = np.linalg.norm(face_area_vectors, axis=1) + normals = face_area_vectors / np.maximum(face_area_magnitudes[:, None], 1.0e-300) + up = 0.5 * (patch_internal + u_values) + up_magnitudes = np.linalg.norm(up, axis=1) + normal_velocity = np.sum(up * normals, axis=1) + value_fraction = np.where(up_magnitudes > 1.0e-300, 0.5 - 0.5 * normal_velocity / up_magnitudes, 0.5) + momentum_u_convection_boundary_source_face_cells_parts.append(face_cells) + momentum_u_convection_boundary_source_parts.append(-phi_values[:, None] * value_fraction[:, None] * u_values) + freestream_internal_coeffs = phi_values[:, None] * (1.0 - value_fraction[:, None]) + momentum_u_convection_boundary_diag_face_cells_parts.append(face_cells) + momentum_u_convection_boundary_diag_coeff_parts.append(np.repeat(freestream_internal_coeffs, 3, axis=1)) + if not np.any(internal_coeffs) and not np.any(boundary_coeffs): + continue + momentum_u_convection_face_cells_parts.append(face_cells) + momentum_u_convection_phi_parts.append(phi_values) + momentum_u_convection_internal_coeff_parts.append(internal_coeffs) + momentum_u_convection_boundary_coeff_parts.append(boundary_coeffs) + n_momentum_u_convection_boundary_faces = int(sum(part.shape[0] for part in momentum_u_convection_face_cells_parts)) + momentum_u_convection_face_cells_np = np.concatenate(momentum_u_convection_face_cells_parts) if momentum_u_convection_face_cells_parts else np.zeros((1,), dtype=np.int32) + momentum_u_convection_phi_np = np.concatenate(momentum_u_convection_phi_parts) if momentum_u_convection_phi_parts else np.zeros((1,), dtype=np.float64) + momentum_u_convection_internal_coeff_np = np.concatenate(momentum_u_convection_internal_coeff_parts) if momentum_u_convection_internal_coeff_parts else np.zeros((1, 3), dtype=np.float64) + momentum_u_convection_boundary_coeff_np = np.concatenate(momentum_u_convection_boundary_coeff_parts) if momentum_u_convection_boundary_coeff_parts else np.zeros((1, 3), dtype=np.float64) + momentum_u_convection_face_cells_gpu, momentum_u_convection_face_cells_np, momentum_u_convection_face_cells_transfer = gpu_i32_array(momentum_u_convection_face_cells_np, "mesh.boundary.face_cells.momentum_convection_coefficients") + momentum_u_convection_phi_gpu, momentum_u_convection_phi_np, momentum_u_convection_phi_transfer = gpu_f64_array(momentum_u_convection_phi_np, "fields.phi.boundary.momentum_convection_coefficients") + momentum_u_convection_internal_coeff_gpu, momentum_u_convection_internal_coeff_np, momentum_u_convection_internal_coeff_transfer = gpu_f64_array(momentum_u_convection_internal_coeff_np, "fields.U.boundary.value_internal_coeffs.momentum_convection") + momentum_u_convection_boundary_coeff_gpu, momentum_u_convection_boundary_coeff_np, momentum_u_convection_boundary_coeff_transfer = gpu_f64_array(momentum_u_convection_boundary_coeff_np, "fields.U.boundary.value_boundary_coeffs.momentum_convection") + n_momentum_u_convection_boundary_source_faces = int(sum(part.shape[0] for part in momentum_u_convection_boundary_source_face_cells_parts)) + momentum_u_convection_boundary_source_face_cells_np = np.concatenate(momentum_u_convection_boundary_source_face_cells_parts) if momentum_u_convection_boundary_source_face_cells_parts else np.zeros((1,), dtype=np.int32) + momentum_u_convection_boundary_source_np = np.concatenate(momentum_u_convection_boundary_source_parts) if momentum_u_convection_boundary_source_parts else np.zeros((1, 3), dtype=np.float64) + momentum_u_convection_boundary_source_face_cells_gpu, momentum_u_convection_boundary_source_face_cells_np, momentum_u_convection_boundary_source_face_cells_transfer = gpu_i32_array(momentum_u_convection_boundary_source_face_cells_np, "mesh.boundary.face_cells.momentum_convection_boundary_source") + momentum_u_convection_boundary_source_gpu, momentum_u_convection_boundary_source_np, momentum_u_convection_boundary_source_transfer = gpu_f64_array(momentum_u_convection_boundary_source_np, "fields.U.boundary.source.momentum_convection") + n_momentum_u_convection_boundary_diag_faces = int(sum(part.shape[0] for part in momentum_u_convection_boundary_diag_face_cells_parts)) + momentum_u_convection_boundary_diag_face_cells_np = np.concatenate(momentum_u_convection_boundary_diag_face_cells_parts) if momentum_u_convection_boundary_diag_face_cells_parts else np.zeros((1,), dtype=np.int32) + momentum_u_convection_boundary_diag_coeff_np = np.concatenate(momentum_u_convection_boundary_diag_coeff_parts) if momentum_u_convection_boundary_diag_coeff_parts else np.zeros((1, 3), dtype=np.float64) + momentum_u_convection_boundary_diag_face_cells_gpu, momentum_u_convection_boundary_diag_face_cells_np, momentum_u_convection_boundary_diag_face_cells_transfer = gpu_i32_array(momentum_u_convection_boundary_diag_face_cells_np, "mesh.boundary.face_cells.momentum_convection_boundary_diag") + momentum_u_convection_boundary_diag_coeff_gpu, momentum_u_convection_boundary_diag_coeff_np, momentum_u_convection_boundary_diag_coeff_transfer = gpu_f64_array(momentum_u_convection_boundary_diag_coeff_np, "fields.U.boundary.internal_coeffs.momentum_convection_boundary_diag") pressure_mixed_face_cells_parts: list[np.ndarray] = [] pressure_mixed_scale_parts: list[np.ndarray] = [] pressure_mixed_value_parts: list[np.ndarray] = [] @@ -1114,29 +1454,6 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P 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] = [] @@ -1248,16 +1565,28 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P 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_source_gpu, u_source_alloc = gpu_empty_f64(tuple(u_np.shape), "gpu_stages.solve_UEqn.rhs_source") + u_unrelaxed_diag_gpu, u_unrelaxed_diag_alloc = gpu_empty_f64((n_cells,), "gpu_stages.UEqn.unrelaxed_diag") + u_unrelaxed_source_gpu, u_unrelaxed_source_alloc = gpu_empty_f64(tuple(u_np.shape), "gpu_stages.UEqn.unrelaxed_source") + u_matrix_source_gpu, u_matrix_source_alloc = gpu_empty_f64(tuple(u_np.shape), "gpu_stages.UEqn.source") + u_boundary_diag_candidate_gpu, u_boundary_diag_candidate_alloc = gpu_empty_f64((n_cells,), "gpu_stages.UEqn.boundary_diag_candidate") + u_boundary_relax_add_gpu, u_boundary_relax_add_alloc = gpu_empty_f64((n_cells,), "gpu_stages.UEqn.boundary_relax_add") + u_boundary_relax_subtract_gpu, u_boundary_relax_subtract_alloc = gpu_empty_f64((n_cells,), "gpu_stages.UEqn.boundary_relax_subtract") + u_boundary_diag_coeff_gpu, u_boundary_diag_coeff_alloc = gpu_empty_f64((n_cells,), "gpu_stages.UEqn.boundary_diag_coeff") 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") + grad_u_gpu, grad_u_alloc = gpu_empty_f64((n_cells, 3, 3), "gpu_stages.UEqn.gradU") rAU_gpu, rAU_alloc = gpu_empty_f64((n_cells,), "gpu_stages.pressure_inputs.rAU") + H1_gpu, H1_alloc = gpu_empty_f64((n_cells,), "gpu_stages.UEqn.H1") 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") + p_solve_diag_gpu, p_solve_diag_alloc = gpu_empty_f64((n_cells,), "gpu_stages.solve_pEqn.spd_diag") + p_solve_source_gpu, p_solve_source_alloc = gpu_empty_f64(tuple(p_np.shape), "gpu_stages.solve_pEqn.spd_source") + p_solve_upper_gpu, p_solve_upper_alloc = gpu_empty_f64((n_internal_faces,), "gpu_stages.solve_pEqn.spd_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") @@ -1290,12 +1619,190 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P pass kernel_graph_started = time.perf_counter() - gpu_rans_momentum_assembly(n_cells, u_gpu, u_diag_gpu, u_source_gpu) + gpu_rans_momentum_assembly(n_cells, u_gpu, cell_volumes_gpu, u_diag_gpu, u_source_gpu) + gpu_rans_zero_scalar_field(n_cells, u_boundary_relax_add_gpu) + gpu_rans_zero_scalar_field(n_cells, u_boundary_relax_subtract_gpu) + gpu_rans_zero_scalar_field(n_cells, u_boundary_diag_coeff_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_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_boundary_relax_add_gpu, u_boundary_relax_subtract_gpu, u_boundary_diag_coeff_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) + gpu_rans_momentum_bounded_convection_sp_internal(n_internal_faces, owner_gpu, neighbour_gpu, phi_gpu, u_diag_gpu) + if n_pressure_boundary_faces: + gpu_rans_momentum_bounded_convection_sp_boundary(n_pressure_boundary_faces, boundary_face_cells_gpu, boundary_phi_gpu, u_diag_gpu) + if n_momentum_u_convection_boundary_faces: + gpu_rans_momentum_convection_boundary_coefficients(n_momentum_u_convection_boundary_faces, momentum_u_convection_face_cells_gpu, momentum_u_convection_phi_gpu, momentum_u_convection_internal_coeff_gpu, momentum_u_convection_boundary_coeff_gpu, u_diag_gpu, u_source_gpu) + gpu_rans_zero_tensor_field(n_cells, grad_u_gpu) + gpu_rans_momentum_gauss_grad_u_internal(n_internal_faces, owner_gpu, neighbour_gpu, u_gpu, sf_gpu, grad_u_gpu) + if n_momentum_u_boundary_faces: + gpu_rans_momentum_gauss_grad_u_boundary(n_momentum_u_boundary_faces, momentum_u_boundary_face_cells_gpu, momentum_u_boundary_values_gpu, momentum_u_boundary_sf_gpu, grad_u_gpu) + gpu_rans_momentum_gauss_grad_u_finish(n_cells, cell_volumes_gpu, grad_u_gpu) + gpu_rans_momentum_linear_upwind_source(n_internal_faces, owner_gpu, neighbour_gpu, phi_gpu, cell_centres_gpu, face_centres_gpu, grad_u_gpu, u_source_gpu) + gpu_copy_scalar(n_cells, u_diag_gpu, u_unrelaxed_diag_gpu) + gpu_copy_vector(n_cells, u_source_gpu, u_unrelaxed_source_gpu) + gpu_rans_zero_scalar_field(n_cells, H1_gpu) + gpu_rans_momentum_offdiag_abs_accumulate(n_internal_faces, owner_gpu, neighbour_gpu, u_upper_gpu, u_lower_gpu, H1_gpu) + if n_momentum_u_convection_boundary_diag_faces: + gpu_rans_momentum_boundary_relaxation_coefficients(n_momentum_u_convection_boundary_diag_faces, momentum_u_convection_boundary_diag_face_cells_gpu, momentum_u_convection_boundary_diag_coeff_gpu, u_boundary_relax_add_gpu, u_boundary_relax_subtract_gpu, u_boundary_diag_coeff_gpu) + gpu_rans_momentum_equation_relaxation(n_cells, u_gpu, DEFAULT_MOMENTUM_RELAXATION_ALPHA, H1_gpu, u_boundary_relax_add_gpu, u_boundary_relax_subtract_gpu, u_diag_gpu, u_source_gpu) + gpu_copy_scalar(n_cells, u_diag_gpu, u_boundary_diag_candidate_gpu) + gpu_rans_add_scalar_field(n_cells, u_boundary_diag_coeff_gpu, u_boundary_diag_candidate_gpu) + gpu_copy_vector(n_cells, u_source_gpu, u_matrix_source_gpu) + if n_momentum_u_convection_boundary_source_faces: + gpu_rans_momentum_convection_boundary_source(n_momentum_u_convection_boundary_source_faces, momentum_u_convection_boundary_source_face_cells_gpu, momentum_u_convection_boundary_source_gpu, u_source_gpu) + gpu_rans_momentum_pressure_gradient_source(n_internal_faces, owner_gpu, neighbour_gpu, p_gpu, sf_gpu, u_source_gpu) + if n_momentum_pressure_boundary_faces: + gpu_rans_momentum_pressure_boundary_source(n_momentum_pressure_boundary_faces, momentum_pressure_boundary_face_cells_gpu, momentum_pressure_boundary_values_gpu, momentum_pressure_boundary_sf_gpu, u_source_gpu) + u_boundary_diag_candidate_for_preconditioner = np.asarray(u_boundary_diag_candidate_gpu.to_numpy(), dtype=np.float64) + u_upper_for_preconditioner = np.asarray(u_upper_gpu.to_numpy(), dtype=np.float64) + u_lower_for_preconditioner = np.asarray(u_lower_gpu.to_numpy(), dtype=np.float64) + u_preconditioner_diag_np = openfoam_dilu_preconditioner_diag( + u_boundary_diag_candidate_for_preconditioner, + owner_np, + neighbour_np, + u_upper_for_preconditioner, + u_lower_for_preconditioner, + ) + u_preconditioner_diag_gpu, u_preconditioner_diag_np, u_preconditioner_diag_transfer = gpu_f64_array( + u_preconditioner_diag_np, + "gpu_stages.solve_UEqn.dilu_preconditioner_diag", + ) + u_dilu_reciprocal_diag_np = np.divide( + 1.0, + u_preconditioner_diag_np, + out=np.ones_like(u_preconditioner_diag_np, dtype=np.float64), + where=np.abs(u_preconditioner_diag_np) > 1.0e-300, + ) + u_dilu_reciprocal_diag_gpu, _, u_dilu_reciprocal_diag_transfer = gpu_f64_array( + u_dilu_reciprocal_diag_np, + "gpu_stages.solve_UEqn.dilu_reciprocal_diag", + ) + gpu_ldu_matvec_vector_asymmetric_faces( + n_cells, + n_internal_faces, + owner_gpu, + neighbour_gpu, + u_upper_gpu, + u_lower_gpu, + u_boundary_diag_candidate_gpu, + u_gpu, + u_operator_intermediate_gpu, + ) + gpu_zero_scalar_accumulator(u_rr_gpu) + gpu_bicgstab_initialize_vector( + n_cells, + u_source_gpu, + u_operator_intermediate_gpu, + u_gpu, + u_solved_gpu, + u_residual_gpu, + u_shadow_gpu, + u_direction_gpu, + u_operator_direction_gpu, + u_rr_gpu, + ) + gpu_bicgstab_precondition_vector(n_cells, u_preconditioner_diag_gpu, u_residual_gpu, u_intermediate_gpu) + gpu_dilu_apply_vector_asymmetric_faces( + n_cells, + n_internal_faces, + owner_gpu, + neighbour_gpu, + losort_gpu, + u_upper_gpu, + u_lower_gpu, + u_dilu_reciprocal_diag_gpu, + u_residual_gpu, + u_operator_direction_gpu, + ) + qd.sync() + u_initial_residual = np.asarray(u_residual_gpu.to_numpy(), dtype=np.float64) + u_initial_residual_squared = float(np.asarray(u_rr_gpu.to_numpy())[0]) + u_diagonal_preconditioned = np.asarray(u_intermediate_gpu.to_numpy(), dtype=np.float64) + u_gpu_dilu_preconditioned = np.asarray(u_operator_direction_gpu.to_numpy(), dtype=np.float64) + u_openfoam_dilu_preconditioned = openfoam_dilu_apply_vector_reference( + owner_np, + neighbour_np, + losort_np, + u_upper_for_preconditioner, + u_lower_for_preconditioner, + u_dilu_reciprocal_diag_np, + u_initial_residual, + ) + gpu_ldu_matvec_vector_asymmetric_faces( + n_cells, + n_internal_faces, + owner_gpu, + neighbour_gpu, + u_upper_gpu, + u_lower_gpu, + u_boundary_diag_candidate_gpu, + u_operator_direction_gpu, + u_intermediate_gpu, + ) + gpu_zero_scalar_accumulator(u_denominator_gpu) + gpu_bicgstab_dot_vector(n_cells, u_shadow_gpu, u_intermediate_gpu, u_denominator_gpu) + qd.sync() + u_operator_preconditioned_direction = np.asarray(u_intermediate_gpu.to_numpy(), dtype=np.float64) + u_first_denominator = float(np.asarray(u_denominator_gpu.to_numpy())[0]) + u_first_alpha = u_initial_residual_squared / u_first_denominator if abs(u_first_denominator) > 1.0e-300 else None + u_gpu_dilu_vs_reference = np.abs(u_gpu_dilu_preconditioned - u_openfoam_dilu_preconditioned) + u_diagonal_vs_reference = np.abs(u_diagonal_preconditioned - u_openfoam_dilu_preconditioned) + u_dilu_preconditioner_diagnostic = { + "status": "measured", + "target": "OpenFOAM DILUPreconditioner::precondition initial residual", + "losort_source": "OpenFOAM lduAddressing::calcLosort replicated from neighbour/upper addressing", + "kernel_entrypoint": "gpu_dilu_apply_vector_asymmetric_faces", + "baseline_kernel_entrypoint": "gpu_bicgstab_precondition_vector", + "n_internal_faces": n_internal_faces, + "n_cells": n_cells, + "residual_entering_preconditioner": array_stats(u_initial_residual), + "same_input_residual_for_reference_and_candidate": True, + "diagonal_candidate_preconditioned": array_stats(u_diagonal_preconditioned), + "openfoam_dilu_reference_preconditioned": array_stats(u_openfoam_dilu_preconditioned), + "gpu_dilu_candidate_preconditioned": array_stats(u_gpu_dilu_preconditioned), + "gpu_dilu_vs_openfoam_reference": array_difference_summary(u_gpu_dilu_preconditioned, u_openfoam_dilu_preconditioned), + "diagonal_vs_openfoam_reference": array_difference_summary(u_diagonal_preconditioned, u_openfoam_dilu_preconditioned), + "diagonal_vs_gpu_dilu_delta_abs": array_stats(np.abs(u_gpu_dilu_preconditioned - u_diagonal_preconditioned)), + "gpu_dilu_vs_reference_delta_l2": finite_float(np.linalg.norm(u_gpu_dilu_vs_reference.reshape(-1))), + "diagonal_vs_reference_delta_l2": finite_float(np.linalg.norm(u_diagonal_vs_reference.reshape(-1))), + } + u_linear_solver_trace = { + "schema_version": 1, + "solver": "gpu_asymmetric_ldu_pbicgstab", + "field": "U", + "preconditioner": "OpenFOAM DILU diagnostic reference plus current GPU candidate", + "trace_points": [ + { + "name": "initial_residual", + "step": "r0 = source - A*x0", + "stats": array_stats(u_initial_residual), + "residual_squared": finite_float(u_initial_residual_squared), + "l2": finite_float(np.linalg.norm(u_initial_residual.reshape(-1))), + }, + { + "name": "preconditioned_residual", + "step": "z0 = M^-1 r0", + "reference": "OpenFOAM DILU forward/back substitution", + "candidate": "gpu_dilu_apply_vector_asymmetric_faces", + "reference_stats": array_stats(u_openfoam_dilu_preconditioned), + "candidate_stats": array_stats(u_gpu_dilu_preconditioned), + "candidate_vs_reference": array_difference_summary(u_gpu_dilu_preconditioned, u_openfoam_dilu_preconditioned), + "current_diagonal_vs_reference": array_difference_summary(u_diagonal_preconditioned, u_openfoam_dilu_preconditioned), + }, + { + "name": "matvec_preconditioned_search_direction", + "step": "A*z0 before alpha", + "stats": array_stats(u_operator_preconditioned_direction), + }, + { + "name": "first_scalar_reductions", + "step": "rho0, denominator0, alpha0", + "rho": finite_float(u_initial_residual_squared), + "denominator": finite_float(u_first_denominator), + "alpha": finite_float(u_first_alpha), + }, + ], + } u_solve_performance = gpu_ldu_pbicgstab_vector_asymmetric_faces( n_cells, n_internal_faces, @@ -1303,7 +1810,7 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P neighbour_gpu, u_upper_gpu, u_lower_gpu, - u_diag_gpu, + u_boundary_diag_candidate_gpu, u_source_gpu, u_gpu, u_solved_gpu, @@ -1318,43 +1825,59 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P u_denominator_gpu, u_omega_numerator_gpu, u_omega_denominator_gpu, + preconditioner_diag=u_preconditioner_diag_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_momentum_hbyA_finish(n_cells, u_boundary_diag_candidate_gpu, HbyA_gpu) + gpu_rans_pressure_inputs(n_cells, u_boundary_diag_candidate_gpu, cell_volumes_gpu, rAU_gpu, H1_gpu) + gpu_rans_momentum_h1_face_accumulate(n_internal_faces, owner_gpu, neighbour_gpu, u_upper_gpu, u_lower_gpu, H1_gpu) + gpu_rans_momentum_h1_finish(n_cells, cell_volumes_gpu, H1_gpu) + gpu_rans_consistent_rAtU(n_cells, rAU_gpu, H1_gpu, rAtU_gpu) gpu_rans_surface_flux_from_cells(n_internal_faces, owner_gpu, neighbour_gpu, HbyA_gpu, sf_gpu, phiHbyA_gpu) + gpu_rans_consistent_phiHbyA_correction(n_internal_faces, owner_gpu, neighbour_gpu, rAU_gpu, rAtU_gpu, p_gpu, cell_centres_gpu, sf_gpu, mag_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_source_from_boundary_flux(n_pressure_boundary_faces, boundary_face_cells_gpu, boundary_phiHbyA_gpu, p_source_gpu) + gpu_rans_negate_scalar_field(n_cells, p_diag_gpu, p_solve_diag_gpu) + gpu_rans_negate_scalar_field(n_cells, p_source_gpu, p_solve_source_gpu) + gpu_rans_negate_scalar_field(n_internal_faces, p_upper_gpu, p_solve_upper_gpu) + p_solve_passes: list[dict[str, Any]] = [] + p_solve_initial_gpu = p_gpu + for non_orthogonal_corrector in range(DEFAULT_PRESSURE_NON_ORTHOGONAL_CORRECTORS + 1): + p_solve_performance = gpu_ldu_pcg_scalar_symmetric_faces( + n_cells, + n_internal_faces, + owner_gpu, + neighbour_gpu, + p_solve_upper_gpu, + p_solve_diag_gpu, + p_solve_source_gpu, + p_solve_initial_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, + ) + p_solve_passes.append({"non_orthogonal_corrector": non_orthogonal_corrector, **p_solve_performance}) + p_solve_initial_gpu = p_solved_gpu + p_solve_performance = { + **p_solve_passes[-1], + "non_orthogonal_correctors": DEFAULT_PRESSURE_NON_ORTHOGONAL_CORRECTORS, + "solve_passes": p_solve_passes, + "total_iterations": sum(int(pass_report["iterations"]) for pass_report in p_solve_passes), + } 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) @@ -1367,10 +1890,18 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P materialization_started = time.perf_counter() u_diag = np.asarray(u_diag_gpu.to_numpy()) + u_boundary_diag_candidate = np.asarray(u_boundary_diag_candidate_gpu.to_numpy()) + u_boundary_diag_coeff = np.asarray(u_boundary_diag_coeff_gpu.to_numpy()) + u_boundary_relax_add = np.asarray(u_boundary_relax_add_gpu.to_numpy()) + u_boundary_relax_subtract = np.asarray(u_boundary_relax_subtract_gpu.to_numpy()) u_source = np.asarray(u_source_gpu.to_numpy()) + u_unrelaxed_diag = np.asarray(u_unrelaxed_diag_gpu.to_numpy()) + u_unrelaxed_source = np.asarray(u_unrelaxed_source_gpu.to_numpy()) + u_matrix_source = np.asarray(u_matrix_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()) + H1 = np.asarray(H1_gpu.to_numpy()) rAtU = np.asarray(rAtU_gpu.to_numpy()) HbyA = np.asarray(HbyA_gpu.to_numpy()) phiHbyA = np.asarray(phiHbyA_gpu.to_numpy()) @@ -1378,7 +1909,9 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P 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()) + u_residual = np.asarray(u_residual_gpu.to_numpy()) p_solved = np.asarray(p_solved_gpu.to_numpy()) + p_residual = np.asarray(p_residual_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()) @@ -1391,7 +1924,7 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P "UEqn", "U", u_diag, - u_source, + u_matrix_source, u_np, kernel="gpu_rans_momentum_assembly", upper=u_upper, @@ -1399,6 +1932,17 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P lower=u_lower, lower_kernel="gpu_rans_momentum_convection_coefficients", ) + UEqn["unrelaxed"] = { + "diag": gpu_array_output("UEqn.unrelaxed_diag", u_unrelaxed_diag, kernel="gpu_copy_scalar"), + "source": gpu_array_output("UEqn.unrelaxed_source", u_unrelaxed_source, kernel="gpu_copy_vector"), + } + UEqn["solve_rhs_source"] = gpu_array_output("UEqn.solve_rhs_source", u_source, kernel="gpu_rans_momentum_pressure_gradient_source") + UEqn["boundary_diag_candidate"] = gpu_array_output("UEqn.boundary_diag_candidate", u_boundary_diag_candidate, kernel="gpu_rans_momentum_boundary_internal_diag") + UEqn["boundary_diag_coeff"] = gpu_array_output("UEqn.boundary_diag_coeff", u_boundary_diag_coeff, kernel="gpu_rans_momentum_boundary_relaxation_coefficients") + UEqn["boundary_relax"] = { + "add": gpu_array_output("UEqn.boundary_relax_add", u_boundary_relax_add, kernel="gpu_rans_momentum_boundary_relaxation_coefficients"), + "subtract": gpu_array_output("UEqn.boundary_relax_subtract", u_boundary_relax_subtract, kernel="gpu_rans_momentum_boundary_relaxation_coefficients"), + } pEqn = gpu_matrix_output( "pEqn", "p", @@ -1410,8 +1954,8 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P 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", + "internal_face_kernel": "gpu_rans_pressure_source_from_flux(phiHbyA)", + "boundary_face_kernel": "gpu_rans_pressure_source_from_boundary_flux(phiHbyA.boundary)", "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, @@ -1426,17 +1970,21 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P ), 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"], + {"terms": [{"name": "gpu_momentum_ddt_diagonal", "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_bounded_convection_sp", "kernel_entrypoints": ["gpu_rans_momentum_bounded_convection_sp_internal", "gpu_rans_momentum_bounded_convection_sp_boundary"], "source": "-fvm::Sp(fvc::surfaceIntegrate(phi), U)"}, {"name": "gpu_momentum_convection_boundary_coefficients", "kernel_entrypoint": "gpu_rans_momentum_convection_boundary_coefficients", "source": "U.boundaryField valueInternalCoeffs/valueBoundaryCoeffs"}, {"name": "gpu_momentum_linear_upwind_correction", "kernel_entrypoints": ["gpu_rans_momentum_gauss_grad_u_internal", "gpu_rans_momentum_gauss_grad_u_boundary", "gpu_rans_momentum_linear_upwind_source"], "source": "bounded Gauss linearUpwind grad(U) explicit correction"}, {"name": "gpu_momentum_equation_relaxation", "kernel_entrypoint": "gpu_rans_momentum_equation_relaxation", "alpha": DEFAULT_MOMENTUM_RELAXATION_ALPHA}, {"name": "gpu_momentum_pressure_gradient_source", "kernel_entrypoint": "gpu_rans_momentum_pressure_gradient_source", "source": "-fvc::grad(p)"}]}, + kernels=["gpu_rans_momentum_assembly", "gpu_rans_momentum_diffusion_coefficients", "gpu_rans_momentum_wall_diffusion_coefficients", "gpu_rans_momentum_convection_coefficients", "gpu_rans_momentum_bounded_convection_sp_internal", "gpu_rans_momentum_bounded_convection_sp_boundary", "gpu_rans_momentum_convection_boundary_coefficients", "gpu_rans_zero_tensor_field", "gpu_rans_momentum_gauss_grad_u_internal", "gpu_rans_momentum_gauss_grad_u_boundary", "gpu_rans_momentum_gauss_grad_u_finish", "gpu_rans_momentum_linear_upwind_source", "gpu_rans_momentum_equation_relaxation", "gpu_rans_momentum_pressure_gradient_source", "gpu_rans_momentum_pressure_boundary_source"], ), - 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("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_bounded_convection_sp_internal", "gpu_rans_momentum_bounded_convection_sp_boundary", "gpu_rans_momentum_convection_boundary_coefficients", "gpu_rans_zero_tensor_field", "gpu_rans_momentum_gauss_grad_u_internal", "gpu_rans_momentum_gauss_grad_u_boundary", "gpu_rans_momentum_gauss_grad_u_finish", "gpu_rans_momentum_linear_upwind_source", "gpu_rans_momentum_equation_relaxation", "gpu_rans_momentum_pressure_gradient_source", "gpu_rans_momentum_pressure_boundary_source"]), gpu_stage_result( "solve_UEqn", { "performance": {"solver_name": "gpu_asymmetric_ldu_pbicgstab", "field_name": "U", **u_solve_performance}, + "rhs_source": gpu_array_output("UEqn.solve_rhs_source", u_source, kernel="gpu_rans_momentum_pressure_gradient_source"), "field_after": gpu_array_output("U", u_solved, kernel="gpu_bicgstab_update_solution_residual_vector_preconditioned"), + "residual": gpu_array_output("U_residual", u_residual, kernel="gpu_bicgstab_update_solution_residual_vector_preconditioned"), + "preconditioner_diagnostic": u_dilu_preconditioner_diagnostic, + "linear_solver_trace": u_linear_solver_trace, }, - 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"], + 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_dilu_apply_vector_asymmetric_faces", "gpu_bicgstab_update_intermediate_vector_preconditioned", "gpu_bicgstab_update_solution_residual_vector_preconditioned"], changed_fields=["U"], ), gpu_stage_result( @@ -1444,21 +1992,47 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P { "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"}, + "HbyA": gpu_array_output("HbyA", HbyA, kernel="gpu_rans_momentum_hbyA_finish"), + "phiHbyA": gpu_array_output("phiHbyA", phiHbyA, kernel="gpu_rans_consistent_phiHbyA_correction"), + "H1": gpu_array_output("H1", H1, kernel="gpu_rans_momentum_h1_finish"), + "HbyA_model": { + "mode": "assembled_UEqn_H_over_A_with_consistent_phiHbyA", + "kernel_entrypoints": [ + "gpu_rans_pressure_inputs", + "gpu_rans_momentum_h1_face_accumulate", + "gpu_rans_momentum_h1_finish", + "gpu_rans_consistent_rAtU", + "gpu_rans_momentum_hbyA_source", + "gpu_rans_momentum_hbyA_face_accumulate", + "gpu_rans_momentum_hbyA_finish", + "gpu_rans_surface_flux_from_cells", + "gpu_rans_consistent_phiHbyA_correction", + ], + "formula": "solveDiag=UEqn.diag+non-coupled boundary internal coeffs; rAU=V/solveDiag; H1=(-UEqn.upper/lower neighbour sum)/V; consistent rAtU=1/max(1/rAU - H1, 0.1/rAU); HbyA=(UEqn.solve_rhs_source - UEqn.offdiag(U))/solveDiag for internal cells; phiHbyA=fvc::flux(HbyA)+interpolate(rAtU-rAU)*snGrad(p)*magSf", + }, }, - 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"], + kernels=[ + "gpu_rans_pressure_inputs", + "gpu_rans_momentum_h1_face_accumulate", + "gpu_rans_momentum_h1_finish", + "gpu_rans_consistent_rAtU", + "gpu_rans_momentum_hbyA_source", + "gpu_rans_momentum_hbyA_face_accumulate", + "gpu_rans_momentum_hbyA_finish", + "gpu_rans_surface_flux_from_cells", + "gpu_rans_consistent_phiHbyA_correction", + ], ), 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}, + "performance": {"solver_name": "gpu_symmetric_ldu_pcg", "field_name": "p", "operator_transform": "negative_openfoam_laplacian_to_spd", "non_orthogonal_loop": "OpenFOAM SIMPLE nNonOrthogonalCorrectors=3", **p_solve_performance}, "p": gpu_array_output("p", p_solved, kernel="gpu_pcg_update_solution_residual_scalar"), + "residual": gpu_array_output("p_residual", p_residual, 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"], + kernels=["gpu_rans_negate_scalar_field", "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( @@ -1524,6 +2098,45 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P "profile_record_count": profiler_evidence.get("profile_record_count"), "timed_scope": "quadrants_cuda_stage_graph_with_input_transfer_and_result_materialization", } + diagnostic_artifacts = None + if diagnostic_artifact_path is not None: + diagnostic_artifacts = write_numeric_artifact_file( + diagnostic_artifact_path, + { + "matrix_operator.UEqn.diag": u_diag, + "matrix_operator.UEqn.upper": u_upper, + "matrix_operator.UEqn.lower": u_lower, + "matrix_operator.UEqn.source": u_matrix_source, + "matrix_operator.UEqn.psi": u_np, + "matrix_operator.pEqn.diag": p_diag, + "matrix_operator.pEqn.upper": p_upper, + "matrix_operator.pEqn.source": p_source, + "matrix_operator.pEqn.psi": p_np, + "pressure_inputs.HbyA": HbyA, + "pressure_inputs.phiHbyA": phiHbyA, + "pressure_inputs.rAU": rAU, + "pressure_inputs.rAtU": rAtU, + "solver.solve_UEqn.field_after": u_solved, + "solver.solve_UEqn.residual": u_residual, + "solver.solve_pEqn.p": p_solved, + "solver.solve_pEqn.phi": phi_solved, + "solver.solve_pEqn.residual": p_residual, + "final_correction.U": u_final, + "final_correction.p": p_final, + "final_correction.phi": phi_solved, + "turbulence.nut": nut_out, + "turbulence.k": k_out, + "turbulence.omega": omega_out, + "fields.U": u_final, + "fields.p": p_final, + "fields.phi": phi_solved, + "fields.nut": nut_out, + "fields.k": k_out, + "fields.omega": omega_out, + }, + role="split", + ) + return { "schema_version": 1, "status": "executed", @@ -1549,13 +2162,22 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P }, "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], + "inputs": [u_transfer, p_transfer, phi_transfer, nut_transfer, k_transfer, omega_transfer, owner_transfer, neighbour_transfer, sf_transfer, face_centres_transfer, cell_centres_transfer, cell_volumes_transfer, mag_sf_transfer, boundary_face_cells_transfer, boundary_phi_transfer, boundary_phiHbyA_transfer, momentum_pressure_boundary_face_cells_transfer, momentum_pressure_boundary_values_transfer, momentum_pressure_boundary_sf_transfer, momentum_u_boundary_face_cells_transfer, momentum_u_boundary_values_transfer, momentum_u_boundary_sf_transfer, momentum_u_convection_face_cells_transfer, momentum_u_convection_phi_transfer, momentum_u_convection_internal_coeff_transfer, momentum_u_convection_boundary_coeff_transfer, pressure_mixed_face_cells_transfer, pressure_mixed_scales_transfer, pressure_mixed_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, u_preconditioner_diag_transfer], "allocations": [ u_diag_alloc, u_source_alloc, + u_unrelaxed_diag_alloc, + u_unrelaxed_source_alloc, + u_matrix_source_alloc, + u_boundary_diag_candidate_alloc, + u_boundary_relax_add_alloc, + u_boundary_relax_subtract_alloc, + u_boundary_diag_coeff_alloc, u_upper_alloc, u_lower_alloc, + grad_u_alloc, rAU_alloc, + H1_alloc, rAtU_alloc, HbyA_alloc, phiHbyA_alloc, @@ -1574,6 +2196,9 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P u_omega_numerator_alloc, u_omega_denominator_alloc, u_residual_alloc, + p_solve_diag_alloc, + p_solve_source_alloc, + p_solve_upper_alloc, p_solved_alloc, p_work_alloc, p_residual_alloc, @@ -1591,6 +2216,7 @@ def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: P }, "profiler": profiler_evidence, "timing": timing, + "diagnostic_artifacts": diagnostic_artifacts, "parity_integration": { "status": "ready", "modes": ["run_one", "split"], @@ -1864,8 +2490,8 @@ def run_backend_iteration(stepper: Any, backend: Mapping[str, Any], case: Path) } -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) +def run_gpu_split_iteration(foam: Any, stepper: Any, backend: Mapping[str, Any], case: Path, *, diagnostic_artifact_path: Path | None = None) -> dict[str, Any]: + gpu_run = run_gpu_solver_stage_smoke(stepper, backend, case, diagnostic_artifact_path=diagnostic_artifact_path) fields = gpu_run["field_objects"] stages = list(gpu_run["stages"]) observability = stage_observability_report( @@ -1894,6 +2520,7 @@ def run_gpu_split_iteration(foam: Any, stepper: Any, backend: Mapping[str, Any], "execution_path": "quadrants_cuda_gpu_rans_split", "backend": dict(backend), "gpu_solver": gpu_solver_stage_report(gpu_run), + "diagnostic_artifacts": gpu_run.get("diagnostic_artifacts"), } __all__ = [ diff --git a/python/src/foam_stepper/gpu/constants.py b/python/src/foam_stepper/gpu/constants.py index c3b0897..5f4cfbb 100644 --- a/python/src/foam_stepper/gpu/constants.py +++ b/python/src/foam_stepper/gpu/constants.py @@ -19,9 +19,10 @@ 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_ITERATIONS = 100 DEFAULT_MOMENTUM_PBICGSTAB_RESIDUAL_TOLERANCE_SQUARED = 1.0e-16 -DEFAULT_PRESSURE_CG_ITERATIONS = 300 +DEFAULT_PRESSURE_CG_ITERATIONS = 2000 +DEFAULT_PRESSURE_NON_ORTHOGONAL_CORRECTORS = 3 DEFAULT_SIMPLE_CONSISTENT_RATU_FACTOR = 10.0 DEFAULT_OMEGA_WALL_BETA1 = 0.075 @@ -34,13 +35,20 @@ STAGE_OBSERVABILITY_GROUPS = ( "assemble_momentum_terms": ("terms",), "assemble_UEqn": ("UEqn",), }, + "run_one_outputs": { + "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"), + "compute_pressure_inputs": ("HbyA", "phiHbyA", "rAU", "rAtU"), + "assemble_pEqn": ("pEqn",), + }, + "run_one_outputs": { + "compute_pressure_inputs": ("HbyA", "phiHbyA", "rAU", "rAtU"), "assemble_pEqn": ("pEqn",), }, }, @@ -52,6 +60,10 @@ STAGE_OBSERVABILITY_GROUPS = ( "solve_UEqn": ("performance", "field_after"), "solve_pEqn": ("performance", "p", "phi"), }, + "run_one_outputs": { + "solve_UEqn": ("performance", "field_after"), + "solve_pEqn": ("performance", "p", "phi"), + }, }, { "name": "final_correction", @@ -60,6 +72,10 @@ STAGE_OBSERVABILITY_GROUPS = ( "split_outputs": { "correct_velocity_pressure_flux": ("U", "p", "phi"), }, + "run_one_outputs": { + "update_phi_from_pEqn_flux": ("phi",), + "correct_velocity_pressure_flux": ("U", "p", "phi"), + }, }, { "name": "turbulence_updates", @@ -69,13 +85,17 @@ STAGE_OBSERVABILITY_GROUPS = ( "momentum_transport_predict": ("case_path", "solver_name"), "momentum_transport_correct": ("U", "p", "phi", "nut", "k", "omega"), }, + "run_one_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"], + "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_bounded_convection_sp_internal", "gpu_rans_momentum_bounded_convection_sp_boundary", "gpu_rans_momentum_convection_boundary_coefficients", "gpu_rans_momentum_convection_boundary_source", "gpu_rans_momentum_boundary_internal_diag", "gpu_rans_momentum_boundary_relaxation_coefficients", "gpu_rans_add_scalar_field", "gpu_rans_zero_tensor_field", "gpu_rans_momentum_gauss_grad_u_internal", "gpu_rans_momentum_gauss_grad_u_boundary", "gpu_rans_momentum_gauss_grad_u_finish", "gpu_rans_momentum_linear_upwind_source", "gpu_rans_zero_scalar_field", "gpu_rans_momentum_offdiag_abs_accumulate", "gpu_rans_momentum_equation_relaxation", "gpu_rans_momentum_pressure_gradient_source", "gpu_rans_momentum_pressure_boundary_source"], + "pressure_assembly": ["gpu_rans_pressure_inputs", "gpu_rans_momentum_h1_face_accumulate", "gpu_rans_momentum_h1_finish", "gpu_rans_consistent_rAtU", "gpu_rans_momentum_hbyA_source", "gpu_rans_momentum_hbyA_face_accumulate", "gpu_rans_momentum_hbyA_finish", "gpu_rans_surface_flux_from_cells", "gpu_rans_consistent_phiHbyA_correction", "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_dilu_apply_vector_asymmetric_faces", "gpu_bicgstab_update_intermediate_vector_preconditioned", "gpu_bicgstab_update_solution_residual_vector_preconditioned", "gpu_vector_residual_squared", "gpu_rans_negate_scalar_field", "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 index 0a9147e..2000d80 100644 --- a/python/src/foam_stepper/gpu/kernels.py +++ b/python/src/foam_stepper/gpu/kernels.py @@ -10,14 +10,16 @@ from .constants import GPU_STAGE_KERNEL_ENTRYPOINTS def gpu_rans_momentum_assembly( n_cells: int, u_internal: qd.types.NDArray[qd.f64, 2], + cell_volumes: qd.types.NDArray[qd.f64, 1], 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] + # fvSchemes uses steadyState ddt; fvm::ddt(U) contributes no diagonal. + diag[cell] = 0.0 + source[cell, 0] = 0.0 + source[cell, 1] = 0.0 + source[cell, 2] = 0.0 @qd.kernel def gpu_rans_momentum_diffusion_coefficients( @@ -64,7 +66,9 @@ def gpu_rans_momentum_wall_diffusion_coefficients( 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], + boundary_relax_add: qd.types.NDArray[qd.f64, 1], + boundary_relax_subtract: qd.types.NDArray[qd.f64, 1], + boundary_diag: qd.types.NDArray[qd.f64, 1], source: qd.types.NDArray[qd.f64, 2], ) -> None: for boundary_face in range(n_boundary_faces): @@ -79,11 +83,15 @@ def gpu_rans_momentum_wall_diffusion_coefficients( 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(boundary_relax_add[cell], coeff) + qd.atomic_add(boundary_relax_subtract[cell], coeff) + qd.atomic_add(boundary_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, @@ -106,25 +114,321 @@ def gpu_rans_momentum_convection_coefficients( qd.atomic_add(upper[face], flux) +@qd.kernel +def gpu_rans_momentum_bounded_convection_sp_internal( + 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], +) -> None: + for face in range(n_internal_faces): + owner_cell = owner[face] + neighbour_cell = neighbour[face] + flux = phi[face] + qd.atomic_add(diag[owner_cell], -flux) + qd.atomic_add(diag[neighbour_cell], flux) + + +@qd.kernel +def gpu_rans_momentum_bounded_convection_sp_boundary( + n_boundary_faces: int, + face_cells: qd.types.NDArray[qd.i32, 1], + phi_boundary: qd.types.NDArray[qd.f64, 1], + diag: qd.types.NDArray[qd.f64, 1], +) -> None: + for face in range(n_boundary_faces): + cell = face_cells[face] + qd.atomic_add(diag[cell], -phi_boundary[face]) + + +@qd.kernel +def gpu_rans_momentum_convection_boundary_coefficients( + n_boundary_faces: int, + face_cells: qd.types.NDArray[qd.i32, 1], + phi_boundary: qd.types.NDArray[qd.f64, 1], + value_internal_coeffs: qd.types.NDArray[qd.f64, 2], + value_boundary_coeffs: qd.types.NDArray[qd.f64, 2], + diag: qd.types.NDArray[qd.f64, 1], + source: qd.types.NDArray[qd.f64, 2], +) -> None: + for face in range(n_boundary_faces): + cell = face_cells[face] + flux = phi_boundary[face] + internal_average = ( + value_internal_coeffs[face, 0] + + value_internal_coeffs[face, 1] + + value_internal_coeffs[face, 2] + ) / 3.0 + qd.atomic_add(diag[cell], flux * internal_average) + qd.atomic_add(source[cell, 0], -flux * value_boundary_coeffs[face, 0]) + qd.atomic_add(source[cell, 1], -flux * value_boundary_coeffs[face, 1]) + qd.atomic_add(source[cell, 2], -flux * value_boundary_coeffs[face, 2]) + + +@qd.kernel +def gpu_rans_momentum_convection_boundary_source( + n_boundary_faces: int, + face_cells: qd.types.NDArray[qd.i32, 1], + boundary_source: qd.types.NDArray[qd.f64, 2], + source: qd.types.NDArray[qd.f64, 2], +) -> None: + for face in range(n_boundary_faces): + cell = face_cells[face] + qd.atomic_add(source[cell, 0], boundary_source[face, 0]) + qd.atomic_add(source[cell, 1], boundary_source[face, 1]) + qd.atomic_add(source[cell, 2], boundary_source[face, 2]) + + +@qd.kernel +def gpu_rans_momentum_boundary_internal_diag( + n_boundary_faces: int, + face_cells: qd.types.NDArray[qd.i32, 1], + internal_coeffs: qd.types.NDArray[qd.f64, 2], + diag: qd.types.NDArray[qd.f64, 1], +) -> None: + for face in range(n_boundary_faces): + cell = face_cells[face] + coeff = ( + internal_coeffs[face, 0] + + internal_coeffs[face, 1] + + internal_coeffs[face, 2] + ) / 3.0 + qd.atomic_add(diag[cell], coeff) + + +@qd.kernel +def gpu_rans_momentum_boundary_relaxation_coefficients( + n_boundary_faces: int, + face_cells: qd.types.NDArray[qd.i32, 1], + internal_coeffs: qd.types.NDArray[qd.f64, 2], + boundary_relax_add: qd.types.NDArray[qd.f64, 1], + boundary_relax_subtract: qd.types.NDArray[qd.f64, 1], + boundary_diag: qd.types.NDArray[qd.f64, 1], +) -> None: + for face in range(n_boundary_faces): + cell = face_cells[face] + c0 = internal_coeffs[face, 0] + c1 = internal_coeffs[face, 1] + c2 = internal_coeffs[face, 2] + a0 = c0 + if a0 < 0.0: + a0 = -a0 + a1 = c1 + if a1 < 0.0: + a1 = -a1 + a2 = c2 + if a2 < 0.0: + a2 = -a2 + max_abs = a0 + if a1 > max_abs: + max_abs = a1 + if a2 > max_abs: + max_abs = a2 + min_coeff = c0 + if c1 < min_coeff: + min_coeff = c1 + if c2 < min_coeff: + min_coeff = c2 + qd.atomic_add(boundary_relax_add[cell], max_abs) + qd.atomic_add(boundary_relax_subtract[cell], min_coeff) + qd.atomic_add(boundary_diag[cell], (c0 + c1 + c2) / 3.0) + + +@qd.kernel +def gpu_rans_add_scalar_field( + n_cells: int, + addend: qd.types.NDArray[qd.f64, 1], + field: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + field[cell] += addend[cell] + + +@qd.kernel +def gpu_rans_zero_tensor_field( + n_cells: int, + tensor: qd.types.NDArray[qd.f64, 3], +) -> None: + for cell in range(n_cells): + for component in range(3): + for direction in range(3): + tensor[cell, component, direction] = 0.0 + + +@qd.kernel +def gpu_rans_momentum_gauss_grad_u_internal( + n_internal_faces: int, + owner: qd.types.NDArray[qd.i32, 1], + neighbour: qd.types.NDArray[qd.i32, 1], + u_internal: qd.types.NDArray[qd.f64, 2], + sf: qd.types.NDArray[qd.f64, 2], + grad_u: qd.types.NDArray[qd.f64, 3], +) -> None: + for face in range(n_internal_faces): + owner_cell = owner[face] + neighbour_cell = neighbour[face] + for component in range(3): + face_value = 0.5 * (u_internal[owner_cell, component] + u_internal[neighbour_cell, component]) + for direction in range(3): + flux_value = face_value * sf[face, direction] + qd.atomic_add(grad_u[owner_cell, component, direction], flux_value) + qd.atomic_add(grad_u[neighbour_cell, component, direction], -flux_value) + + +@qd.kernel +def gpu_rans_momentum_gauss_grad_u_boundary( + n_boundary_faces: int, + face_cells: qd.types.NDArray[qd.i32, 1], + u_boundary: qd.types.NDArray[qd.f64, 2], + sf_boundary: qd.types.NDArray[qd.f64, 2], + grad_u: qd.types.NDArray[qd.f64, 3], +) -> None: + for face in range(n_boundary_faces): + cell = face_cells[face] + for component in range(3): + face_value = u_boundary[face, component] + for direction in range(3): + qd.atomic_add(grad_u[cell, component, direction], face_value * sf_boundary[face, direction]) + + +@qd.kernel +def gpu_rans_momentum_gauss_grad_u_finish( + n_cells: int, + cell_volumes: qd.types.NDArray[qd.f64, 1], + grad_u: qd.types.NDArray[qd.f64, 3], +) -> None: + for cell in range(n_cells): + inv_volume = 1.0 / cell_volumes[cell] + for component in range(3): + for direction in range(3): + grad_u[cell, component, direction] *= inv_volume + + +@qd.kernel +def gpu_rans_momentum_linear_upwind_source( + 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], + cell_centres: qd.types.NDArray[qd.f64, 2], + face_centres: qd.types.NDArray[qd.f64, 2], + grad_u: qd.types.NDArray[qd.f64, 3], + source: qd.types.NDArray[qd.f64, 2], +) -> None: + for face in range(n_internal_faces): + owner_cell = owner[face] + neighbour_cell = neighbour[face] + flux = phi[face] + upwind_cell = owner_cell + if flux < 0.0: + upwind_cell = neighbour_cell + dx0 = face_centres[face, 0] - cell_centres[upwind_cell, 0] + dx1 = face_centres[face, 1] - cell_centres[upwind_cell, 1] + dx2 = face_centres[face, 2] - cell_centres[upwind_cell, 2] + for component in range(3): + correction = ( + dx0 * grad_u[upwind_cell, component, 0] + + dx1 * grad_u[upwind_cell, component, 1] + + dx2 * grad_u[upwind_cell, component, 2] + ) + flux_correction = flux * correction + qd.atomic_add(source[owner_cell, component], -flux_correction) + qd.atomic_add(source[neighbour_cell, component], flux_correction) + + +@qd.kernel +def gpu_rans_zero_scalar_field( + n_cells: int, + field: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + field[cell] = 0.0 + + +@qd.kernel +def gpu_rans_momentum_offdiag_abs_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], + offdiag_sum: qd.types.NDArray[qd.f64, 1], +) -> None: + for face in range(n_internal_faces): + upper_abs = upper[face] + if upper_abs < 0.0: + upper_abs = -upper_abs + lower_abs = lower[face] + if lower_abs < 0.0: + lower_abs = -lower_abs + qd.atomic_add(offdiag_sum[owner[face]], upper_abs) + qd.atomic_add(offdiag_sum[neighbour[face]], lower_abs) + @qd.kernel def gpu_rans_momentum_equation_relaxation( n_cells: int, u_internal: qd.types.NDArray[qd.f64, 2], alpha: float, + offdiag_sum: qd.types.NDArray[qd.f64, 1], + boundary_relax_add: qd.types.NDArray[qd.f64, 1], + boundary_relax_subtract: qd.types.NDArray[qd.f64, 1], 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 + raw_diag = diag[cell] + dominant_diag = raw_diag + boundary_relax_add[cell] + if dominant_diag < 0.0: + dominant_diag = -dominant_diag + if offdiag_sum[cell] > dominant_diag: + dominant_diag = offdiag_sum[cell] + relaxed_diag = dominant_diag / alpha - boundary_relax_subtract[cell] + source_scale = relaxed_diag - raw_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_pressure_gradient_source( + n_internal_faces: int, + owner: qd.types.NDArray[qd.i32, 1], + neighbour: qd.types.NDArray[qd.i32, 1], + p_internal: qd.types.NDArray[qd.f64, 1], + sf: qd.types.NDArray[qd.f64, 2], + source: qd.types.NDArray[qd.f64, 2], +) -> None: + for face in range(n_internal_faces): + owner_cell = owner[face] + neighbour_cell = neighbour[face] + p_face = 0.5 * (p_internal[owner_cell] + p_internal[neighbour_cell]) + qd.atomic_add(source[owner_cell, 0], -p_face * sf[face, 0]) + qd.atomic_add(source[owner_cell, 1], -p_face * sf[face, 1]) + qd.atomic_add(source[owner_cell, 2], -p_face * sf[face, 2]) + qd.atomic_add(source[neighbour_cell, 0], p_face * sf[face, 0]) + qd.atomic_add(source[neighbour_cell, 1], p_face * sf[face, 1]) + qd.atomic_add(source[neighbour_cell, 2], p_face * sf[face, 2]) + + +@qd.kernel +def gpu_rans_momentum_pressure_boundary_source( + n_boundary_faces: int, + face_cells: qd.types.NDArray[qd.i32, 1], + p_boundary: qd.types.NDArray[qd.f64, 1], + sf_boundary: qd.types.NDArray[qd.f64, 2], + source: qd.types.NDArray[qd.f64, 2], +) -> None: + for face in range(n_boundary_faces): + cell = face_cells[face] + p_face = p_boundary[face] + qd.atomic_add(source[cell, 0], -p_face * sf_boundary[face, 0]) + qd.atomic_add(source[cell, 1], -p_face * sf_boundary[face, 1]) + qd.atomic_add(source[cell, 2], -p_face * sf_boundary[face, 2]) + + @qd.kernel def gpu_rans_momentum_hbyA_source( n_cells: int, @@ -173,18 +477,7 @@ def gpu_rans_momentum_hbyA_finish( 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 @@ -193,20 +486,53 @@ def gpu_rans_pressure_inputs( u_diag: qd.types.NDArray[qd.f64, 1], cell_volumes: qd.types.NDArray[qd.f64, 1], rAU: qd.types.NDArray[qd.f64, 1], + H1: qd.types.NDArray[qd.f64, 1], ) -> None: for cell in range(n_cells): rAU[cell] = cell_volumes[cell] / u_diag[cell] + H1[cell] = 0.0 + + +@qd.kernel +def gpu_rans_momentum_h1_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], + H1: qd.types.NDArray[qd.f64, 1], +) -> None: + for face in range(n_internal_faces): + owner_cell = owner[face] + neighbour_cell = neighbour[face] + qd.atomic_add(H1[owner_cell], -upper[face]) + qd.atomic_add(H1[neighbour_cell], -lower[face]) + + +@qd.kernel +def gpu_rans_momentum_h1_finish( + n_cells: int, + cell_volumes: qd.types.NDArray[qd.f64, 1], + H1: qd.types.NDArray[qd.f64, 1], +) -> None: + for cell in range(n_cells): + H1[cell] /= cell_volumes[cell] @qd.kernel def gpu_rans_consistent_rAtU( n_cells: int, rAU: qd.types.NDArray[qd.f64, 1], - ratu_factor: float, + H1: qd.types.NDArray[qd.f64, 1], rAtU: qd.types.NDArray[qd.f64, 1], ) -> None: for cell in range(n_cells): - rAtU[cell] = ratu_factor * rAU[cell] + inv_rAU = 1.0 / rAU[cell] + floor = 0.1 * inv_rAU + denominator = inv_rAU - H1[cell] + if denominator < floor: + denominator = floor + rAtU[cell] = 1.0 / denominator @qd.kernel @@ -228,6 +554,39 @@ def gpu_rans_surface_flux_from_cells( ) + +@qd.kernel +def gpu_rans_consistent_phiHbyA_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], + rAtU: qd.types.NDArray[qd.f64, 1], + p_internal: 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], + phiHbyA: 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 + interpolated_delta = 0.5 * ( + (rAtU[owner_cell] - rAU[owner_cell]) + + (rAtU[neighbour_cell] - rAU[neighbour_cell]) + ) + pressure_jump = p_internal[neighbour_cell] - p_internal[owner_cell] + phiHbyA[face] += interpolated_delta * pressure_jump * mag_sf[face] * mag_sf[face] / projected_delta + + @qd.kernel def gpu_rans_pressure_assembly( n_cells: int, @@ -311,6 +670,16 @@ def gpu_rans_pressure_mixed_boundary_laplacian( qd.atomic_add(source[cell], -coeff * boundary_values[face]) +@qd.kernel +def gpu_rans_negate_scalar_field( + n_values: int, + source: qd.types.NDArray[qd.f64, 1], + out: qd.types.NDArray[qd.f64, 1], +) -> None: + for index in range(n_values): + out[index] = -source[index] + + @qd.kernel def gpu_rans_face_flux_copy( n_internal_faces: int, @@ -424,17 +793,34 @@ __all__ = [ "gpu_rans_momentum_assembly", "gpu_rans_momentum_diffusion_coefficients", "gpu_rans_momentum_convection_coefficients", + "gpu_rans_momentum_bounded_convection_sp_internal", + "gpu_rans_momentum_bounded_convection_sp_boundary", + "gpu_rans_momentum_convection_boundary_coefficients", + "gpu_rans_momentum_convection_boundary_source", + "gpu_rans_momentum_boundary_internal_diag", + "gpu_rans_momentum_boundary_relaxation_coefficients", + "gpu_rans_add_scalar_field", + "gpu_rans_zero_tensor_field", + "gpu_rans_momentum_gauss_grad_u_internal", + "gpu_rans_momentum_gauss_grad_u_boundary", + "gpu_rans_momentum_gauss_grad_u_finish", + "gpu_rans_momentum_linear_upwind_source", "gpu_rans_momentum_equation_relaxation", + "gpu_rans_momentum_pressure_gradient_source", "gpu_rans_momentum_hbyA_source", "gpu_rans_momentum_hbyA_face_accumulate", "gpu_rans_momentum_hbyA_finish", "gpu_rans_pressure_inputs", + "gpu_rans_momentum_h1_face_accumulate", + "gpu_rans_momentum_h1_finish", "gpu_rans_consistent_rAtU", + "gpu_rans_consistent_phiHbyA_correction", "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_negate_scalar_field", "gpu_rans_face_flux_copy", "gpu_rans_surface_flux_from_cells", "gpu_rans_pressure_flux_correction", diff --git a/python/src/foam_stepper/gpu/linear_solve.py b/python/src/foam_stepper/gpu/linear_solve.py index 1e7f664..1a8f51f 100644 --- a/python/src/foam_stepper/gpu/linear_solve.py +++ b/python/src/foam_stepper/gpu/linear_solve.py @@ -21,6 +21,20 @@ class GpuLduCsr: metadata: Mapping[str, Any] +@dataclasses.dataclass(frozen=True) +class GpuLduLevelSchedule: + """Device level schedule for triangular sweeps over OpenFOAM LDU faces.""" + + level_offsets: Any + level_cells: Any + incoming_offsets: Any + incoming_faces: Any + outgoing_offsets: Any + outgoing_faces: Any + n_levels: int + metadata: Mapping[str, Any] + + @qd.kernel def gpu_ldu_jacobi_scalar( n_cells: int, @@ -332,16 +346,18 @@ def gpu_ldu_pcg_scalar_symmetric_faces( operator_work: Any, residual_squared: Any, denominator: Any, + preconditioner_diag: Any | None = None, *, iterations: int, residual_tolerance_squared: float = 0.0, ) -> dict[str, Any]: - """Run diagonal-preconditioned GPU CG for a symmetric per-face LDU matrix.""" + """Run GPU CG for a symmetric per-face LDU matrix with a diagonal preconditioner.""" + precond_diag = diag if preconditioner_diag is None else preconditioner_diag 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) + gpu_pcg_initialize_scalar(n_cells, source, operator_work, precond_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]) @@ -362,7 +378,7 @@ def gpu_ldu_pcg_scalar_symmetric_faces( 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) + gpu_pcg_update_solution_residual_scalar(n_cells, alpha, precond_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]) @@ -382,7 +398,7 @@ def gpu_ldu_pcg_scalar_symmetric_faces( "final_preconditioned_dot": rho_value, "converged": rr_value <= residual_tolerance_squared, "residual_tolerance_squared": residual_tolerance_squared, - "preconditioner": "diagonal_jacobi", + "preconditioner": "diagonal_jacobi" if preconditioner_diag is None else "dic_reciprocal_diagonal", } def gpu_ldu_jacobi_scalar_symmetric_faces( n_cells: int, @@ -552,6 +568,8 @@ def gpu_bicgstab_dot_vector( ) + + @qd.kernel def gpu_bicgstab_update_direction_vector( n_cells: int, @@ -642,6 +660,148 @@ def gpu_bicgstab_precondition_vector( out[cell, 2] = source[cell, 2] +@qd.kernel +def gpu_dilu_apply_vector_asymmetric_faces( + n_cells: int, + n_internal_faces: int, + owner: qd.types.NDArray[qd.i32, 1], + neighbour: qd.types.NDArray[qd.i32, 1], + losort: qd.types.NDArray[qd.i32, 1], + upper: qd.types.NDArray[qd.f64, 1], + lower: qd.types.NDArray[qd.f64, 1], + reciprocal_diag: qd.types.NDArray[qd.f64, 1], + source: qd.types.NDArray[qd.f64, 2], + out: qd.types.NDArray[qd.f64, 2], +) -> None: + """Apply OpenFOAM DILU forward/back substitution to one vector residual.""" + + for worker in range(1): + for cell in range(n_cells): + scale = reciprocal_diag[cell] + out[cell, 0] = scale * source[cell, 0] + out[cell, 1] = scale * source[cell, 1] + out[cell, 2] = scale * source[cell, 2] + + for sorted_index in range(n_internal_faces): + face = losort[sorted_index] + owner_cell = owner[face] + neighbour_cell = neighbour[face] + scale = reciprocal_diag[neighbour_cell] * lower[face] + out[neighbour_cell, 0] -= scale * out[owner_cell, 0] + out[neighbour_cell, 1] -= scale * out[owner_cell, 1] + out[neighbour_cell, 2] -= scale * out[owner_cell, 2] + + for reverse_index in range(n_internal_faces): + face = n_internal_faces - 1 - reverse_index + owner_cell = owner[face] + neighbour_cell = neighbour[face] + scale = reciprocal_diag[owner_cell] * upper[face] + out[owner_cell, 0] -= scale * out[neighbour_cell, 0] + out[owner_cell, 1] -= scale * out[neighbour_cell, 1] + out[owner_cell, 2] -= scale * out[neighbour_cell, 2] + +@qd.kernel +def gpu_dilu_forward_level_vector( + level: int, + level_offsets: qd.types.NDArray[qd.i32, 1], + level_cells: qd.types.NDArray[qd.i32, 1], + owner: qd.types.NDArray[qd.i32, 1], + incoming_offsets: qd.types.NDArray[qd.i32, 1], + incoming_faces: qd.types.NDArray[qd.i32, 1], + lower: qd.types.NDArray[qd.f64, 1], + reciprocal_diag: qd.types.NDArray[qd.f64, 1], + source: qd.types.NDArray[qd.f64, 2], + out: qd.types.NDArray[qd.f64, 2], +) -> None: + for index in range(level_offsets[level], level_offsets[level + 1]): + cell = level_cells[index] + scale = reciprocal_diag[cell] + value_x = scale * source[cell, 0] + value_y = scale * source[cell, 1] + value_z = scale * source[cell, 2] + for face_slot in range(incoming_offsets[cell], incoming_offsets[cell + 1]): + face = incoming_faces[face_slot] + owner_cell = owner[face] + coeff = scale * lower[face] + value_x -= coeff * out[owner_cell, 0] + value_y -= coeff * out[owner_cell, 1] + value_z -= coeff * out[owner_cell, 2] + out[cell, 0] = value_x + out[cell, 1] = value_y + out[cell, 2] = value_z + + +@qd.kernel +def gpu_dilu_backward_level_vector( + level: int, + level_offsets: qd.types.NDArray[qd.i32, 1], + level_cells: qd.types.NDArray[qd.i32, 1], + neighbour: qd.types.NDArray[qd.i32, 1], + outgoing_offsets: qd.types.NDArray[qd.i32, 1], + outgoing_faces: qd.types.NDArray[qd.i32, 1], + upper: qd.types.NDArray[qd.f64, 1], + reciprocal_diag: qd.types.NDArray[qd.f64, 1], + out: qd.types.NDArray[qd.f64, 2], +) -> None: + for index in range(level_offsets[level], level_offsets[level + 1]): + cell = level_cells[index] + scale = reciprocal_diag[cell] + value_x = out[cell, 0] + value_y = out[cell, 1] + value_z = out[cell, 2] + begin = outgoing_offsets[cell] + end = outgoing_offsets[cell + 1] + for reverse_slot in range(end - begin): + face = outgoing_faces[end - 1 - reverse_slot] + neighbour_cell = neighbour[face] + coeff = scale * upper[face] + value_x -= coeff * out[neighbour_cell, 0] + value_y -= coeff * out[neighbour_cell, 1] + value_z -= coeff * out[neighbour_cell, 2] + out[cell, 0] = value_x + out[cell, 1] = value_y + out[cell, 2] = value_z + + + +def gpu_dilu_apply_vector_asymmetric_levels( + schedule: GpuLduLevelSchedule, + owner: Any, + neighbour: Any, + upper: Any, + lower: Any, + reciprocal_diag: Any, + source: Any, + out: Any, +) -> None: + """Apply OpenFOAM DILU using parallel cell work within each dependency level.""" + + for level in range(schedule.n_levels): + gpu_dilu_forward_level_vector( + level, + schedule.level_offsets, + schedule.level_cells, + owner, + schedule.incoming_offsets, + schedule.incoming_faces, + lower, + reciprocal_diag, + source, + out, + ) + for reverse_level in range(schedule.n_levels): + level = schedule.n_levels - 1 - reverse_level + gpu_dilu_backward_level_vector( + level, + schedule.level_offsets, + schedule.level_cells, + neighbour, + schedule.outgoing_offsets, + schedule.outgoing_faces, + upper, + reciprocal_diag, + out, + ) @qd.kernel def gpu_bicgstab_update_intermediate_vector_preconditioned( n_cells: int, @@ -1032,11 +1192,12 @@ def gpu_ldu_pbicgstab_vector_asymmetric_faces( denominator: Any, omega_numerator: Any, omega_denominator: Any, + preconditioner_diag: Any | None = None, *, iterations: int, residual_tolerance_squared: float = 0.0, ) -> dict[str, Any]: - """Run diagonal-preconditioned GPU BiCGStab for an asymmetric vector LDU matrix.""" + """Run GPU BiCGStab for an asymmetric vector LDU matrix with a diagonal preconditioner.""" gpu_ldu_matvec_vector_asymmetric_faces( n_cells, @@ -1068,6 +1229,7 @@ def gpu_ldu_pbicgstab_vector_asymmetric_faces( alpha = 1.0 omega = 1.0 performed_iterations = 0 + precond_diag = diag if preconditioner_diag is None else preconditioner_diag for _ in range(iterations): if residual_value <= residual_tolerance_squared: @@ -1088,8 +1250,19 @@ def gpu_ldu_pbicgstab_vector_asymmetric_faces( 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_bicgstab_precondition_vector(n_cells, precond_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() @@ -1116,7 +1289,7 @@ def gpu_ldu_pbicgstab_vector_asymmetric_faces( residual_value = intermediate_residual break - gpu_bicgstab_precondition_vector(n_cells, diag, intermediate, residual) + gpu_bicgstab_precondition_vector(n_cells, precond_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) @@ -1149,7 +1322,7 @@ def gpu_ldu_pbicgstab_vector_asymmetric_faces( "final_residual_squared": residual_value, "converged": residual_value <= residual_tolerance_squared, "residual_tolerance_squared": residual_tolerance_squared, - "preconditioner": "diagonal_jacobi", + "preconditioner": "diagonal_jacobi" if preconditioner_diag is None else "dilu_reciprocal_diagonal", } def gpu_ldu_jacobi_vector_symmetric_faces( @@ -1228,6 +1401,28 @@ def build_ldu_csr( return offsets, columns, coefficients +def build_losort_addr(n_cells: int, neighbour: np.ndarray) -> np.ndarray: + """Build OpenFOAM lduAddressing::losortAddr from upper/neighbour cells.""" + + neighbour_i32 = np.asarray(neighbour, dtype=np.int32).reshape(-1) + if neighbour_i32.size == 0: + return np.zeros(0, dtype=np.int32) + if int(neighbour_i32.min()) < 0 or int(neighbour_i32.max()) >= n_cells: + raise ValueError("neighbour addresses exceed cell range") + + counts = np.bincount(neighbour_i32, minlength=n_cells).astype(np.int32, copy=False) + offsets = np.empty(n_cells + 1, dtype=np.int32) + offsets[0] = 0 + np.cumsum(counts, out=offsets[1:]) + losort = np.empty(int(offsets[-1]), dtype=np.int32) + cursor = offsets[:-1].copy() + for face, neighbour_cell in enumerate(neighbour_i32): + slot = int(cursor[int(neighbour_cell)]) + losort[slot] = int(face) + cursor[int(neighbour_cell)] += 1 + return np.ascontiguousarray(losort, dtype=np.int32) + + 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,) @@ -1274,6 +1469,7 @@ def empty_ldu_csr_gpu(n_cells: int, name: str) -> GpuLduCsr: __all__ = [ "GpuLduCsr", "build_ldu_csr", + "build_losort_addr", "empty_ldu_csr_gpu", "gpu_bicgstab_dot_vector", "gpu_bicgstab_initialize_vector", @@ -1285,6 +1481,7 @@ __all__ = [ "gpu_bicgstab_update_solution_residual_vector_preconditioned", "gpu_copy_scalar", "gpu_copy_vector", + "gpu_dilu_apply_vector_asymmetric_faces", "gpu_ldu_bicgstab_vector_asymmetric_faces", "gpu_ldu_pbicgstab_vector_asymmetric_faces", "gpu_ldu_cg_scalar_symmetric_faces", @@ -1306,5 +1503,6 @@ __all__ = [ "gpu_scalar_jacobi_finish", "gpu_vector_jacobi_finish", "gpu_vector_residual_squared", + "gpu_zero_scalar_accumulator", "ldu_csr_to_gpu", ] diff --git a/scripts/update_loop_diagnostic_context.py b/scripts/update_loop_diagnostic_context.py new file mode 100755 index 0000000..a8ac3c8 --- /dev/null +++ b/scripts/update_loop_diagnostic_context.py @@ -0,0 +1,430 @@ +#!/usr/bin/env python3 +"""Generate concise loop context from GPU RANS verifier reports.""" + +from __future__ import annotations + +import argparse +import json +import math +from datetime import datetime, timezone +from pathlib import Path +from typing import Any, Mapping + +ROOT = Path(__file__).resolve().parents[1] +DEFAULT_CONTEXT = ROOT / ".loop/diagnostic-context.md" +DEFAULT_BASELINE = ROOT / ".loop/diagnostic-baseline.json" +REPORT_GLOBS = ( + "tmp/**/verifier_report.json", + "tmp/**/report.json", +) +TMP_REPORT_GLOBS = ( + "*gpu*/report.json", + "*gpu*/verifier_report.json", + "worker_gpu_rans_solver*/report.json", + "judge_gpu_rans_solver*/report.json", + "*gpu*rans*/report.json", + "*gpu*rans*/verifier_report.json", +) + + +def load_json(path: Path) -> Any | None: + try: + return json.loads(path.read_text()) + except Exception: + return None + + +def is_verifier_report(value: Any) -> bool: + if not isinstance(value, Mapping): + return False + return any(key in value for key in ("verifier_evidence", "first_divergence_summary", "artifact_comparisons")) + + +def is_gpu_report(value: Mapping[str, Any], path: Path) -> bool: + backend = value.get("backend") if isinstance(value.get("backend"), Mapping) else {} + requested = str(backend.get("requested") or "").lower() + selected = str(backend.get("selected") or "").lower() + return requested == "gpu" or selected == "gpu" or "gpu" in path.parent.name.lower() + + +def discover_latest_report(root: Path) -> tuple[Path | None, Mapping[str, Any] | None]: + candidates: list[Path] = [] + for pattern in REPORT_GLOBS: + candidates.extend(root.glob(pattern)) + tmp_root = Path("/tmp") + if tmp_root.exists(): + for pattern in TMP_REPORT_GLOBS: + candidates.extend(tmp_root.glob(pattern)) + unique = sorted({path.resolve() for path in candidates if path.is_file()}, key=lambda path: path.stat().st_mtime, reverse=True) + reports: list[tuple[Path, Mapping[str, Any]]] = [] + for path in unique: + data = load_json(path) + if is_verifier_report(data): + reports.append((path, data)) # type: ignore[arg-type] + if not reports: + return None, None + gpu_reports = [(path, data) for path, data in reports if is_gpu_report(data, path)] + return (gpu_reports or reports)[0] + + +def get_path(value: Mapping[str, Any], path: str) -> Any: + cursor: Any = value + for part in path.split("."): + if not isinstance(cursor, Mapping): + return None + cursor = cursor.get(part) + return cursor + + +def finite_number(value: Any) -> float | None: + try: + number = float(value) + except (TypeError, ValueError): + return None + return number if math.isfinite(number) else None + + +def fmt(value: Any) -> str: + number = finite_number(value) + if number is None: + if value is True: + return "true" + if value is False: + return "false" + if value is None: + return "-" + return str(value) + if number == 0.0: + return "0" + if abs(number) >= 1e4 or abs(number) < 1e-3: + return f"{number:.3e}" + return f"{number:.6g}" + + +def status_word(value: Any) -> str: + if value is True: + return "passed" + if value is False: + return "failed" + return str(value or "unknown") + + +def first_nested_key(value: Any, target: str, path: str = "") -> tuple[str, Mapping[str, Any]] | None: + if isinstance(value, Mapping): + for key, item in value.items(): + next_path = f"{path}.{key}" if path else str(key) + if key == target and isinstance(item, Mapping): + return next_path, item + found = first_nested_key(item, target, next_path) + if found is not None: + return found + elif isinstance(value, list): + for index, item in enumerate(value): + found = first_nested_key(item, target, f"{path}[{index}]") + if found is not None: + return found + return None + + +def artifact_checks(report: Mapping[str, Any]) -> list[dict[str, Any]]: + out: list[dict[str, Any]] = [] + comparisons = report.get("artifact_comparisons", {}) if isinstance(report.get("artifact_comparisons"), Mapping) else {} + for family in ("pressure_inputs", "matrix_operator", "solver"): + family_report = comparisons.get(family) if isinstance(comparisons.get(family), Mapping) else {} + for check in family_report.get("checks", []) if isinstance(family_report.get("checks"), list) else []: + if not isinstance(check, Mapping): + continue + largest = check.get("largest_difference") if isinstance(check.get("largest_difference"), Mapping) else {} + location = largest.get("location") if isinstance(largest.get("location"), Mapping) else {} + out.append( + { + "family": family, + "name": check.get("name"), + "allclose": check.get("allclose"), + "reason": check.get("reason"), + "max_abs": check.get("max_abs"), + "mean_abs": check.get("mean_abs"), + "rms_abs": check.get("rms_abs"), + "location": location, + } + ) + out.sort(key=lambda item: (item.get("allclose") is True, item["family"], str(item.get("name")))) + return out + + +def field_checks(report: Mapping[str, Any]) -> list[dict[str, Any]]: + out: list[dict[str, Any]] = [] + modes = report.get("modes", {}) if isinstance(report.get("modes"), Mapping) else {} + for mode_name in ("run_one", "split"): + mode = modes.get(mode_name) if isinstance(modes.get(mode_name), Mapping) else {} + comparisons = mode.get("comparisons", {}) if isinstance(mode.get("comparisons"), Mapping) else {} + for field, data in comparisons.items(): + if not isinstance(data, Mapping): + continue + out.append( + { + "mode": mode_name, + "field": field, + "allclose": data.get("allclose"), + "max_abs": data.get("max_abs"), + "mean_abs": data.get("mean_abs"), + "rms_abs": data.get("rms_abs"), + "location": data.get("location"), + } + ) + out.sort(key=lambda item: (item.get("allclose") is True, item["mode"], str(item.get("field")))) + return out + + +def extract_metrics(report: Mapping[str, Any]) -> dict[str, Any]: + metrics: dict[str, Any] = { + "status.passed": report.get("status") == "passed", + "verifier.passed": get_path(report, "verifier_evidence.passed") is True, + } + first = report.get("first_divergence_summary") if isinstance(report.get("first_divergence_summary"), Mapping) else {} + if first: + metrics["first.target"] = first.get("first_target") + metrics["first.family"] = first.get("artifact_family") + metrics["first.stage_group"] = first.get("stage_group") + for check in artifact_checks(report): + base = f"artifact.{check['family']}.{check['name']}" + metrics[f"{base}.allclose"] = check.get("allclose") is True + for key in ("max_abs", "mean_abs", "rms_abs"): + number = finite_number(check.get(key)) + if number is not None: + metrics[f"{base}.{key}"] = number + for check in field_checks(report): + base = f"field.{check['mode']}.{check['field']}" + metrics[f"{base}.allclose"] = check.get("allclose") is True + for key in ("max_abs", "mean_abs", "rms_abs"): + number = finite_number(check.get(key)) + if number is not None: + metrics[f"{base}.{key}"] = number + preconditioner = first_nested_key(report, "preconditioner_diagnostic") + if preconditioner is not None: + _, data = preconditioner + for key in ( + "gpu_dilu_vs_reference_delta_l2", + "diagonal_vs_reference_delta_l2", + ): + number = finite_number(data.get(key)) + if number is not None: + metrics[f"preconditioner.{key}"] = number + for path, metric_name in ( + ("gpu_dilu_vs_openfoam_reference.max_abs", "preconditioner.gpu_dilu_vs_reference.max_abs"), + ("gpu_dilu_vs_openfoam_reference.rms_abs", "preconditioner.gpu_dilu_vs_reference.rms_abs"), + ("diagonal_vs_openfoam_reference.max_abs", "preconditioner.diagonal_vs_reference.max_abs"), + ("diagonal_vs_openfoam_reference.rms_abs", "preconditioner.diagonal_vs_reference.rms_abs"), + ): + number = finite_number(get_path(data, path)) + if number is not None: + metrics[metric_name] = number + return metrics + + +def classify_delta(current: Mapping[str, Any], baseline: Mapping[str, Any] | None) -> list[dict[str, Any]]: + if not baseline: + return [{"metric": name, "status": "newly_available", "current": value, "baseline": None} for name, value in sorted(current.items())] + out: list[dict[str, Any]] = [] + previous = baseline.get("metrics", {}) if isinstance(baseline.get("metrics"), Mapping) else {} + for name, value in sorted(current.items()): + old = previous.get(name) + status = "newly_available" + if old is not None: + if isinstance(value, bool) and isinstance(old, bool): + if value == old: + status = "unchanged" + elif value and not old: + status = "improved" + else: + status = "regressed" + else: + now_num = finite_number(value) + old_num = finite_number(old) + if now_num is not None and old_num is not None: + tolerance = max(1e-15, abs(old_num) * 1e-9) + if abs(now_num - old_num) <= tolerance: + status = "unchanged" + elif now_num < old_num: + status = "improved" + else: + status = "regressed" + else: + status = "unchanged" if value == old else "changed" + out.append({"metric": name, "status": status, "current": value, "baseline": old}) + return out + + +def metric_priority(item: Mapping[str, Any]) -> tuple[int, str]: + status_order = {"regressed": 0, "improved": 1, "newly_available": 2, "changed": 3, "unchanged": 4} + return status_order.get(str(item.get("status")), 9), str(item.get("metric")) + + +def render_location(location: Any) -> str: + if not isinstance(location, Mapping): + return "-" + entity = location.get("entity_kind") or "array" + index = location.get("entity_index") + component = location.get("component_index") + if component is None: + return f"{entity}[{index}]" + return f"{entity}[{index}] component={component}" + + +def render_context(report_path: Path | None, report: Mapping[str, Any] | None, baseline: Mapping[str, Any] | None) -> str: + generated = datetime.now(timezone.utc).isoformat() + if report is None or report_path is None: + return "\n".join( + [ + "# GPU RANS Loop Diagnostic Context", + "", + f"Generated: {generated}", + "", + "No verifier report was found under repo tmp/ or /tmp GPU RANS work directories.", + "Next action: run `scripts/verify_gpu_rans_solver.sh --work --report /report.json` or the focused verifier, then rerun this hook.", + "", + ] + ) + + first = report.get("first_divergence_summary") if isinstance(report.get("first_divergence_summary"), Mapping) else {} + artifacts = report.get("intermediate_artifacts", {}) if isinstance(report.get("intermediate_artifacts"), Mapping) else {} + families = artifacts.get("families", []) if isinstance(artifacts.get("families"), list) else [] + metrics = extract_metrics(report) + deltas = classify_delta(metrics, baseline) + checks = artifact_checks(report) + fields = field_checks(report) + preconditioner = first_nested_key(report, "preconditioner_diagnostic") + solver_trace = first_nested_key(report, "linear_solver_trace") + + lines = [ + "# GPU RANS Loop Diagnostic Context", + "", + f"Generated: {generated}", + f"Latest report: `{report_path}`", + f"Report status: `{report.get('status')}`", + f"Verifier evidence passed: `{get_path(report, 'verifier_evidence.passed')}`", + "", + "## First divergence", + "", + f"- Target: `{first.get('first_target') if first else None}`", + f"- Artifact family: `{first.get('artifact_family') if first else None}`", + f"- Stage group: `{first.get('stage_group') if first else None}`", + f"- Evidence path: `{first.get('evidence_path') if first else None}`", + f"- Field/reason: `{first.get('field') if first else None}` / `{first.get('reason') if first else None}`", + "", + "## Solver phase evidence", + "", + "| Family | Status | Why |", + "|---|---:|---|", + ] + for family in families: + if not isinstance(family, Mapping): + continue + lines.append(f"| {family.get('name')} | {family.get('status')} | {family.get('why')} |") + + lines.extend(["", "## Numeric artifact comparisons", "", "| Family | Check | Status | max_abs | mean_abs | rms_abs | Location |", "|---|---|---:|---:|---:|---:|---|"]) + for check in checks[:16]: + lines.append( + "| {family} | {name} | {status} | {max_abs} | {mean_abs} | {rms_abs} | {location} |".format( + family=check.get("family"), + name=check.get("name"), + status=status_word(check.get("allclose")), + max_abs=fmt(check.get("max_abs")), + mean_abs=fmt(check.get("mean_abs")), + rms_abs=fmt(check.get("rms_abs")), + location=render_location(check.get("location")), + ) + ) + + lines.extend(["", "## Field comparison symptoms", "", "| Mode | Field | Status | max_abs | mean_abs | rms_abs | Location |", "|---|---|---:|---:|---:|---:|---|"]) + for check in fields[:12]: + lines.append( + "| {mode} | {field} | {status} | {max_abs} | {mean_abs} | {rms_abs} | {location} |".format( + mode=check.get("mode"), + field=check.get("field"), + status=status_word(check.get("allclose")), + max_abs=fmt(check.get("max_abs")), + mean_abs=fmt(check.get("mean_abs")), + rms_abs=fmt(check.get("rms_abs")), + location=render_location(check.get("location")), + ) + ) + + lines.extend(["", "## Linear solver and preconditioner trace", ""]) + if preconditioner is None: + lines.append("No preconditioner diagnostic artifact found in the latest report.") + else: + path, data = preconditioner + lines.extend( + [ + f"Preconditioner evidence path: `{path}`", + f"- GPU DILU vs OpenFOAM reference max_abs: `{fmt(get_path(data, 'gpu_dilu_vs_openfoam_reference.max_abs'))}`", + f"- GPU DILU vs OpenFOAM reference rms_abs: `{fmt(get_path(data, 'gpu_dilu_vs_openfoam_reference.rms_abs'))}`", + f"- Diagonal/current vs OpenFOAM reference max_abs: `{fmt(get_path(data, 'diagonal_vs_openfoam_reference.max_abs'))}`", + f"- Residual entering preconditioner recorded: `{data.get('residual_entering_preconditioner') is not None}`", + ] + ) + if solver_trace is not None: + path, data = solver_trace + lines.append(f"Solver trace path: `{path}`") + for point in data.get("trace_points", []) if isinstance(data.get("trace_points"), list) else []: + if isinstance(point, Mapping): + lines.append(f"- `{point.get('name')}`: {point.get('step')}") + + lines.extend(["", "## Delta versus retained baseline", "", "| Metric | Status | Current | Baseline |", "|---|---:|---:|---:|"]) + for item in sorted(deltas, key=metric_priority)[:24]: + lines.append(f"| `{item.get('metric')}` | {item.get('status')} | {fmt(item.get('current'))} | {fmt(item.get('baseline'))} |") + + lines.extend( + [ + "", + "## Next target hint", + "", + f"Focus first on `{first.get('first_target') if first else 'unknown'}`. Treat downstream field symptoms as unreliable until that artifact or missing evidence closes.", + "", + ] + ) + return "\n".join(lines) + + +def parse_args(argv: list[str] | None = None) -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--root", type=Path, default=ROOT) + parser.add_argument("--report", type=Path, default=None, help="Explicit report path; otherwise discover newest verifier report") + parser.add_argument("--context", type=Path, default=DEFAULT_CONTEXT) + parser.add_argument("--baseline", type=Path, default=DEFAULT_BASELINE) + return parser.parse_args(argv) + + +def main(argv: list[str] | None = None) -> int: + args = parse_args(argv) + root = args.root.resolve() + if args.report is not None: + report_path = args.report.resolve() + loaded = load_json(report_path) + report = loaded if is_verifier_report(loaded) else None + else: + report_path, report = discover_latest_report(root) + + baseline = load_json(args.baseline) if args.baseline.exists() else None + baseline_mapping = baseline if isinstance(baseline, Mapping) else None + args.context.parent.mkdir(parents=True, exist_ok=True) + args.context.write_text(render_context(report_path, report, baseline_mapping), encoding="utf-8") + + if report is not None and report_path is not None: + current = { + "updated_at": datetime.now(timezone.utc).isoformat(), + "report": str(report_path), + "metrics": extract_metrics(report), + "first_divergence_summary": report.get("first_divergence_summary"), + } + args.baseline.parent.mkdir(parents=True, exist_ok=True) + args.baseline.write_text(json.dumps(current, indent=2, sort_keys=True) + "\n", encoding="utf-8") + print(f"updated diagnostic context: {args.context} from {report_path}") + else: + print(f"updated diagnostic context without report: {args.context}") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/verify_airfrans_stepper.py b/scripts/verify_airfrans_stepper.py index fe30c2d..24dc80f 100755 --- a/scripts/verify_airfrans_stepper.py +++ b/scripts/verify_airfrans_stepper.py @@ -46,6 +46,7 @@ DEFAULT_WORK = ROOT / "tmp/airfrans_stepper_verify" PRIMARY_FIELDS = ("U", "p", "phi") TURBULENCE_FIELDS = ("nut", "k", "omega") REQUIRED_FIELDS = REQUIRED_COMPARISON_FIELDS +REQUIRED_EXECUTION_MODES = ("run_one", "split") BACKEND_CHOICES = ("auto", "cpu", "gpu") FULL_GPU_RANS_GUARD = "scripts/verify_gpu_rans_solver.sh" GPU_PRIMITIVE_PROOF = "scripts/verify_gpu_algorithm.sh" @@ -61,6 +62,7 @@ STEPPER_FAILURE = "stepper_execution_failure" COMPARISON_FAILURE = "numerical_comparison_failure" BACKEND_FAILURE = "backend_execution_failure" INTERNAL_FAILURE = "internal_harness_failure" +REQUIRED_EVIDENCE_FAILURE = "required_evidence_failure" EXIT_CODES = { LOADING_FAILURE: 2, @@ -68,6 +70,7 @@ EXIT_CODES = { STEPPER_FAILURE: 4, COMPARISON_FAILURE: 5, BACKEND_FAILURE: 6, + REQUIRED_EVIDENCE_FAILURE: 8, INTERNAL_FAILURE: 7, } @@ -80,13 +83,20 @@ STAGE_OBSERVABILITY_GROUPS = ( "assemble_momentum_terms": ("terms",), "assemble_UEqn": ("UEqn",), }, + "run_one_outputs": { + "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"), + "compute_pressure_inputs": ("HbyA", "phiHbyA", "rAU", "rAtU"), + "assemble_pEqn": ("pEqn",), + }, + "run_one_outputs": { + "compute_pressure_inputs": ("HbyA", "phiHbyA", "rAU", "rAtU"), "assemble_pEqn": ("pEqn",), }, }, @@ -98,6 +108,10 @@ STAGE_OBSERVABILITY_GROUPS = ( "solve_UEqn": ("performance", "field_after"), "solve_pEqn": ("performance", "p", "phi"), }, + "run_one_outputs": { + "solve_UEqn": ("performance", "field_after"), + "solve_pEqn": ("performance", "p", "phi"), + }, }, { "name": "final_correction", @@ -106,6 +120,10 @@ STAGE_OBSERVABILITY_GROUPS = ( "split_outputs": { "correct_velocity_pressure_flux": ("U", "p", "phi"), }, + "run_one_outputs": { + "update_phi_from_pEqn_flux": ("phi",), + "correct_velocity_pressure_flux": ("U", "p", "phi"), + }, }, { "name": "turbulence_updates", @@ -115,9 +133,16 @@ STAGE_OBSERVABILITY_GROUPS = ( "momentum_transport_predict": ("case_path", "solver_name"), "momentum_transport_correct": ("U", "p", "phi", "nut", "k", "omega"), }, + "run_one_outputs": { + "momentum_transport_predict": ("case_path", "solver_name"), + "momentum_transport_correct": ("U", "p", "phi", "nut", "k", "omega"), + }, }, ) +STAGE_GROUP_ORDER = {group["name"]: index for index, group in enumerate(STAGE_OBSERVABILITY_GROUPS)} +MODE_COMPARISON_ORDER = {"run_one": 0, "split": 1} + FIELD_COMPARISON_ATTRIBUTION = { "U": ("final_correction", "linear_solve_results", "momentum_assembly", "turbulence_updates"), "p": ("linear_solve_results", "pressure_assembly", "final_correction"), @@ -127,6 +152,73 @@ FIELD_COMPARISON_ATTRIBUTION = { "omega": ("turbulence_updates",), } +INTERMEDIATE_ARTIFACT_FAMILIES = ( + { + "name": "identity", + "stage_groups": (), + "reference_paths": ("case_preparation.prepared_case_identity", "mesh_identity.oracle"), + "candidate_paths": ("mesh_identity.run_one", "mesh_identity.split"), + "comparison_paths": ("modes.run_one.mesh_comparison", "modes.split.mesh_comparison"), + "required": True, + }, + { + "name": "pressure_inputs", + "stage_groups": ("pressure_assembly",), + "reference_paths": ("modes.reference_split.stages.compute_pressure_inputs",), + "candidate_paths": ("modes.split.stages.compute_pressure_inputs",), + "comparison_paths": ("artifact_comparisons.pressure_inputs",), + "required": True, + }, + { + "name": "matrix_operator", + "stage_groups": ("momentum_assembly", "pressure_assembly"), + "reference_paths": ("modes.reference_split.momentum_terms", "modes.reference_split.UEqn", "modes.reference_split.stages.relax_UEqn.outputs.UEqn", "modes.reference_split.pEqn", "state_exports.reference_split_matrices"), + "candidate_paths": ("modes.split.momentum_terms", "modes.split.UEqn", "modes.split.pEqn", "state_exports.split_matrices"), + "comparison_paths": ("artifact_comparisons.matrix_operator",), + "required": True, + }, + { + "name": "solver", + "stage_groups": ("linear_solve_results",), + "reference_paths": ("modes.reference_split.stages.solve_UEqn", "modes.reference_split.stages.solve_pEqn"), + "candidate_paths": ("modes.split.stages.solve_UEqn", "modes.split.stages.solve_pEqn"), + "comparison_paths": ("artifact_comparisons.solver",), + "required": True, + }, + { + "name": "correction", + "stage_groups": ("final_correction",), + "reference_paths": ("modes.oracle.fields.U", "modes.oracle.fields.p"), + "candidate_paths": ("modes.split.stages.correct_velocity_pressure_flux", "modes.split.fields.U", "modes.split.fields.p"), + "comparison_paths": ("modes.split.comparisons.U", "modes.split.comparisons.p"), + "required": True, + }, + { + "name": "flux", + "stage_groups": ("pressure_assembly", "final_correction"), + "reference_paths": ("modes.oracle.fields.phi",), + "candidate_paths": ("modes.split.stages.compute_pressure_inputs", "modes.split.fields.phi"), + "comparison_paths": ("modes.split.comparisons.phi",), + "required": True, + }, + { + "name": "turbulence", + "stage_groups": ("turbulence_updates",), + "reference_paths": ("modes.oracle.fields.nut", "modes.oracle.fields.k", "modes.oracle.fields.omega"), + "candidate_paths": ( + "modes.split.stages.momentum_transport_predict", + "modes.split.stages.momentum_transport_correct", + "modes.split.fields.nut", + "modes.split.fields.k", + "modes.split.fields.omega", + ), + "comparison_paths": ("modes.split.comparisons.nut", "modes.split.comparisons.k", "modes.split.comparisons.omega"), + "required": True, + }, +) + +ARTIFACT_FAMILY_ORDER = {family["name"]: index for index, family in enumerate(INTERMEDIATE_ARTIFACT_FAMILIES)} + GPU_RUNTIME_UNAVAILABLE = "gpu_runtime_unavailable" GPU_NUMERICAL_MISMATCH = "gpu_numerical_mismatch" GPU_INPUT_SCHEMA_VERSION = 1 @@ -262,6 +354,33 @@ def json_ready(value: Any) -> Any: return [json_ready(item) for item in value] return repr(value) +def report_path_value(report: Mapping[str, Any], path: str) -> Any: + cursor: Any = report + for part in path.split("."): + if isinstance(cursor, Mapping): + cursor = cursor.get(part) + continue + if isinstance(cursor, list): + matched = None + for item in cursor: + if isinstance(item, Mapping) and item.get("name") == part: + matched = item + break + cursor = matched + continue + return None + return cursor + + +def report_path_exists(report: Mapping[str, Any], path: str) -> bool: + value = report_path_value(report, path) + if value is None: + return False + if isinstance(value, Mapping) or isinstance(value, list): + return bool(value) + return True + + def array_shape(array: Any | None) -> list[int] | None: if array is None: @@ -284,6 +403,9 @@ def array_stats(array: Any) -> dict[str, Any]: "dtype": str(arr.dtype), "size": int(arr.size), } + if arr.dtype != object: + contiguous = np.ascontiguousarray(arr) + out["sha256"] = hashlib.sha256(contiguous.tobytes()).hexdigest() if arr.size == 0 or not np.issubdtype(arr.dtype, np.number): return out @@ -464,6 +586,9 @@ def stage_observability_report(mode: str, stages: list[dict[str, Any]], *, evide 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", {}) + backend = evidence.get("backend", {}) if isinstance(evidence.get("backend"), Mapping) else {} + if mode == "run_one" and backend.get("selected") != "gpu": + output_requirements = {} output_checks = {} for stage_name, required_outputs in output_requirements.items(): stage = stage_by_name.get(stage_name) @@ -641,7 +766,376 @@ def field_comparison_attribution( 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]: +def safe_array_item(array: Any, index: int) -> Any: + arr = np.asarray(array) + if index < 0 or index >= arr.shape[0]: + return None + value = arr[index] + if isinstance(value, np.ndarray): + return value.tolist() + if isinstance(value, np.generic): + return value.item() + return json_ready(value) + + +def cell_neighbourhood_context(mesh: Any, cell_index: int) -> dict[str, Any]: + owner = np.asarray(getattr(mesh, "owner", [])) + neighbour = np.asarray(getattr(mesh, "neighbour", [])) + owner_faces = np.flatnonzero(owner == cell_index) + neighbour_faces = np.flatnonzero(neighbour == cell_index) + touching_faces = sorted({int(face) for face in np.concatenate((owner_faces, neighbour_faces))})[:12] + boundary_patches = [] + for patch in getattr(mesh, "boundary", []): + face_cells = np.asarray(getattr(patch, "face_cells", [])) + hits = np.flatnonzero(face_cells == cell_index) + if hits.size: + boundary_patches.append( + { + "patch_index": int(getattr(patch, "index", 0)), + "patch_name": getattr(patch, "name", ""), + "patch_type": getattr(patch, "type", ""), + "patch_local_entries": [int(item) for item in hits[:8]], + "face_indices": [int(np.asarray(getattr(patch, "face_indices", []))[item]) for item in hits[:8]], + } + ) + return { + "entity": "cell", + "cell_index": cell_index, + "volume": safe_array_item(getattr(mesh, "V", []), cell_index), + "centre": safe_array_item(getattr(mesh, "C", []), cell_index), + "touching_internal_faces": touching_faces, + "touching_internal_face_count": int(owner_faces.size + neighbour_faces.size), + "boundary_patches": boundary_patches, + } + + +def face_neighbourhood_context(mesh: Any, face_index: int) -> dict[str, Any]: + owner = np.asarray(getattr(mesh, "owner", [])) + neighbour = np.asarray(getattr(mesh, "neighbour", [])) + context = { + "entity": "face", + "face_index": face_index, + "owner": int(owner[face_index]) if 0 <= face_index < owner.shape[0] else None, + "neighbour": int(neighbour[face_index]) if 0 <= face_index < neighbour.shape[0] else None, + "centre": safe_array_item(getattr(mesh, "Cf", []), face_index), + "area_vector": safe_array_item(getattr(mesh, "Sf", []), face_index), + "area_magnitude": safe_array_item(getattr(mesh, "magSf", []), face_index), + "patch": None, + } + for patch in getattr(mesh, "boundary", []): + face_indices = np.asarray(getattr(patch, "face_indices", [])) + hits = np.flatnonzero(face_indices == face_index) + if hits.size: + context["patch"] = { + "patch_index": int(getattr(patch, "index", 0)), + "patch_name": getattr(patch, "name", ""), + "patch_type": getattr(patch, "type", ""), + "patch_local_entry": int(hits[0]), + "face_cell": safe_array_item(getattr(patch, "face_cells", []), int(hits[0])), + } + break + return context + + +def local_entity_context(mesh: Any | None, location: Mapping[str, Any] | None, field_name: str) -> dict[str, Any] | None: + if mesh is None or not location: + return None + entity_kind = location.get("entity_kind") + entity_index = location.get("entity_index") + if not isinstance(entity_index, int): + return None + if entity_kind == "cell": + context = cell_neighbourhood_context(mesh, entity_index) + elif entity_kind in {"face", "internal_face"}: + context = face_neighbourhood_context(mesh, entity_index) + else: + context = {"entity": entity_kind, "entity_index": entity_index} + context["field"] = field_name + context["array_index"] = location.get("array_index") + context["component_index"] = location.get("component_index") + return context + + +NUMERIC_ARTIFACT_SCHEMA_VERSION = 1 +NUMERIC_ARTIFACT_TOP_N = 5 + + +def numeric_artifact_array(value: Any) -> np.ndarray | None: + if value is None: + return None + arr = np.asarray(value) + if arr.dtype == object or not np.issubdtype(arr.dtype, np.number): + return None + return np.ascontiguousarray(arr) + + +def add_numeric_artifact(arrays: dict[str, np.ndarray], name: str, value: Any) -> None: + arr = numeric_artifact_array(value) + if arr is not None: + arrays[name] = arr + + +def add_field_numeric_artifact(arrays: dict[str, np.ndarray], name: str, field: Any) -> None: + if field is not None and hasattr(field, "internal"): + add_numeric_artifact(arrays, name, getattr(field, "internal")) + + +def add_matrix_numeric_artifacts(arrays: dict[str, np.ndarray], prefix: str, matrix: Any) -> None: + if matrix is None: + return + for component in ("diag", "upper", "lower", "source"): + value = getattr(matrix, component, None) + if value is not None: + add_numeric_artifact(arrays, f"{prefix}.{component}", value) + psi = getattr(matrix, "psi", None) + if psi is not None and hasattr(psi, "internal"): + add_numeric_artifact(arrays, f"{prefix}.psi", getattr(psi, "internal")) + + +def write_numeric_artifact_file(path: Path, arrays: Mapping[str, np.ndarray], *, role: str) -> dict[str, Any]: + materialized = {name: np.ascontiguousarray(value) for name, value in arrays.items()} + path.parent.mkdir(parents=True, exist_ok=True) + np.savez_compressed(path, **materialized) + return { + "schema_version": NUMERIC_ARTIFACT_SCHEMA_VERSION, + "role": role, + "path": str(path), + "format": "npz", + "artifact_count": len(materialized), + "artifacts": {name: array_stats(value) for name, value in materialized.items()}, + } + + +def stage_output(stage: Any, name: str) -> Any: + outputs = getattr(stage, "outputs", None) + if isinstance(outputs, Mapping): + return outputs.get(name) + return None + + +def write_split_numeric_artifacts( + path: Path | None, + *, + role: str, + relax_UEqn: Any, + solve_UEqn: Any, + compute_pressure_inputs: Any, + assemble_pEqn: Any, + solve_pEqn: Any, + correct_velocity_pressure_flux: Any, + momentum_transport_correct: Any, + fields: Mapping[str, Any], +) -> dict[str, Any] | None: + if path is None: + return None + + arrays: dict[str, np.ndarray] = {} + add_matrix_numeric_artifacts(arrays, "matrix_operator.UEqn", stage_output(relax_UEqn, "UEqn")) + add_matrix_numeric_artifacts(arrays, "matrix_operator.pEqn", stage_output(assemble_pEqn, "pEqn")) + + for name in ("HbyA", "phiHbyA", "rAU", "rAtU"): + add_field_numeric_artifact(arrays, f"pressure_inputs.{name}", stage_output(compute_pressure_inputs, name)) + + add_field_numeric_artifact(arrays, "solver.solve_UEqn.field_before", stage_output(solve_UEqn, "field_before")) + add_field_numeric_artifact(arrays, "solver.solve_UEqn.field_after", stage_output(solve_UEqn, "field_after")) + add_matrix_numeric_artifacts(arrays, "solver.solve_UEqn.matrix_before", stage_output(solve_UEqn, "matrix_before")) + add_field_numeric_artifact(arrays, "solver.solve_pEqn.field_before", stage_output(solve_pEqn, "field_before")) + add_field_numeric_artifact(arrays, "solver.solve_pEqn.p", stage_output(solve_pEqn, "p")) + add_field_numeric_artifact(arrays, "solver.solve_pEqn.phi", stage_output(solve_pEqn, "phi")) + add_matrix_numeric_artifacts(arrays, "solver.solve_pEqn.matrix_before", stage_output(solve_pEqn, "matrix_before")) + + for name in ("U", "p", "phi"): + add_field_numeric_artifact(arrays, f"final_correction.{name}", stage_output(correct_velocity_pressure_flux, name)) + for name in ("nut", "k", "omega"): + add_field_numeric_artifact(arrays, f"turbulence.{name}", stage_output(momentum_transport_correct, name)) + for name in REQUIRED_FIELDS: + add_field_numeric_artifact(arrays, f"fields.{name}", fields.get(name)) + + return write_numeric_artifact_file(path, arrays, role=role) + + +def load_numeric_artifact_file(path: Path | None) -> dict[str, np.ndarray] | None: + if path is None or not path.exists(): + return None + try: + with np.load(path, allow_pickle=False) as data: + return {name: np.asarray(data[name]) for name in data.files} + except Exception: + return None + + +def diagnostic_artifact_path(report: Mapping[str, Any], role: str) -> Path | None: + value = report_path_value(report, f"diagnostic_artifacts.{role}.path") + if value is None: + return None + return Path(str(value)) + + +def artifact_entity_kind(mesh: Any | None, shape: tuple[int, ...]) -> str | None: + if mesh is None or not shape: + return None + first_dim = int(shape[0]) + if first_dim == int(getattr(mesh, "n_cells", -1)): + return "cell" + if first_dim == int(getattr(mesh, "n_internal_faces", -1)): + return "internal_face" + return None + + +def artifact_location(name: str, index: tuple[int, ...], shape: tuple[int, ...], mesh: Any | None) -> dict[str, Any]: + entity_kind = artifact_entity_kind(mesh, shape) or "array" + entity_index = int(index[0]) if index else None + component_index = list(index[1:]) if len(index) > 1 else None + return { + "array_index": list(index), + "entity_kind": entity_kind, + "entity_index": entity_index, + "component_index": component_index, + "artifact": name, + } + + +def top_numeric_differences( + name: str, + actual: np.ndarray, + expected: np.ndarray, + diff: np.ndarray, + *, + mesh_context: Any | None, + limit: int = NUMERIC_ARTIFACT_TOP_N, +) -> list[dict[str, Any]]: + if diff.size == 0: + return [] + flat = diff.reshape(-1) + finite = np.isfinite(flat) + if np.any(finite): + finite_indices = np.flatnonzero(finite) + finite_values = flat[finite_indices] + if finite_values.size > limit: + local = np.argpartition(finite_values, -limit)[-limit:] + candidate_indices = finite_indices[local] + else: + candidate_indices = finite_indices + ordered = sorted(candidate_indices, key=lambda item: (-float(flat[int(item)]), int(item))) + else: + ordered = [int(item) for item in np.flatnonzero(~finite)[:limit]] + + out = [] + for flat_index in ordered[:limit]: + index = tuple(int(item) for item in np.unravel_index(int(flat_index), diff.shape)) + location = artifact_location(name, index, tuple(diff.shape), mesh_context) + out.append( + { + "location": location, + "abs_difference": finite_float(diff[index]), + "actual_value": value_at(actual, index), + "reference_value": value_at(expected, index), + "local_context": local_entity_context(mesh_context, location, name), + } + ) + return out + + +def compare_numeric_artifact_arrays( + name: str, + key: str, + reference: np.ndarray | None, + candidate: np.ndarray | None, + *, + rtol: float, + atol: float, + mesh_context: Any | None, +) -> dict[str, Any]: + out: dict[str, Any] = { + "name": name, + "artifact_key": key, + "reference_available": reference is not None, + "candidate_available": candidate is not None, + "comparison": "numpy.allclose plus absolute-difference diagnostics", + "rtol": rtol, + "atol": atol, + } + if reference is None or candidate is None: + return {**out, "shape_matches": False, "allclose": False, "reason": "missing_numeric_artifact"} + + ref = np.asarray(reference) + actual = np.asarray(candidate) + shape_matches = ref.shape == actual.shape + dtype_matches = ref.dtype == actual.dtype + out.update( + { + "reference_shape": array_shape(ref), + "candidate_shape": array_shape(actual), + "shape_matches": shape_matches, + "reference_dtype": str(ref.dtype), + "candidate_dtype": str(actual.dtype), + "dtype_matches": dtype_matches, + } + ) + if not shape_matches: + return {**out, "allclose": False, "reason": "shape_mismatch"} + if ref.size == 0: + return { + **out, + "allclose": True, + "reason": None, + "max_abs": 0.0, + "mean_abs": 0.0, + "rms_abs": 0.0, + "largest_difference": None, + "top_differences": [], + } + + diff = np.abs(actual - ref) + finite = np.isfinite(diff) + finite_diff = diff[finite] + if finite_diff.size: + max_flat = int(np.argmax(np.where(finite, diff, -np.inf))) + max_abs = finite_float(diff.reshape(-1)[max_flat]) + mean_abs = finite_float(np.mean(finite_diff)) + rms_abs = finite_float(np.sqrt(np.mean(np.square(finite_diff)))) + else: + max_flat = int(np.flatnonzero(~finite.reshape(-1))[0]) + max_abs = None + mean_abs = None + rms_abs = None + + max_index = tuple(int(item) for item in np.unravel_index(max_flat, diff.shape)) + location = artifact_location(name, max_index, tuple(diff.shape), mesh_context) + expected_at_max = value_at(ref, max_index) + tolerance_at_max = atol + rtol * abs(float(expected_at_max)) if isinstance(expected_at_max, (int, float)) else None + allclose = bool(np.allclose(actual, ref, rtol=rtol, atol=atol, equal_nan=False)) + out.update( + { + "allclose": allclose, + "reason": None if allclose else "value_mismatch", + "max_abs": max_abs, + "mean_abs": mean_abs, + "rms_abs": rms_abs, + "nonfinite_error_count": int(diff.size - np.count_nonzero(finite)), + "largest_difference": { + "location": location, + "actual_value": value_at(actual, max_index), + "reference_value": expected_at_max, + "actual_entity_value": entity_value_at(actual, location.get("entity_index")), + "reference_entity_value": entity_value_at(ref, location.get("entity_index")), + "tolerance_at_max": finite_float(tolerance_at_max), + "local_context": local_entity_context(mesh_context, location, name), + }, + "top_differences": top_numeric_differences(name, actual, ref, diff, mesh_context=mesh_context), + } + ) + return out + +def field_compare_report( + name: str, + actual_field: Any, + expected_field: Any, + *, + rtol: float, + atol: float, + mesh_context: Any | None = None, +) -> 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", "") @@ -669,9 +1163,11 @@ def field_compare_report(name: str, actual_field: Any, expected_field: Any, *, r "allclose": True, "max_abs": 0.0, "mean_abs": 0.0, + "rms_abs": 0.0, "location": None, "actual_at_max": None, "expected_at_max": None, + "top_differences": [], } ) return report, None @@ -696,23 +1192,27 @@ def field_compare_report(name: str, actual_field: Any, expected_field: Any, *, r finite_diff = diff[finite] allclose = bool(np.allclose(actual, expected, rtol=rtol, atol=atol, equal_nan=False)) + location = { + "array_index": list(max_index), + "entity_kind": getattr(actual_field, "entity_kind", ""), + "entity_index": entity_index, + "component_index": component_index, + } report.update( { "allclose": allclose, "max_abs": max_abs, "mean_abs": finite_float(np.mean(finite_diff)) if finite_diff.size else None, - "location": { - "array_index": list(max_index), - "entity_kind": getattr(actual_field, "entity_kind", ""), - "entity_index": entity_index, - "component_index": component_index, - }, + "rms_abs": finite_float(np.sqrt(np.mean(np.square(finite_diff)))) if finite_diff.size else None, + "location": location, "actual_at_max": actual_at_max, "expected_at_max": expected_at_max, "actual_entity_at_max": entity_value_at(actual, entity_index), "expected_entity_at_max": entity_value_at(expected, entity_index), + "local_context": local_entity_context(mesh_context, location, name), "tolerance_at_max": finite_float(tolerance_at_max), "nonfinite_error_count": int(diff.size - np.count_nonzero(finite)), + "top_differences": top_numeric_differences(name, actual, expected, diff, mesh_context=mesh_context), } ) if allclose: @@ -729,6 +1229,7 @@ def compare_fields( rtol: float, atol: float, stage_graph: Iterable[str] | None = None, + mesh_context: Any | None = None, ) -> tuple[dict[str, Any], list[dict[str, Any]]]: report: dict[str, Any] = {} mismatches: list[dict[str, Any]] = [] @@ -746,7 +1247,7 @@ def compare_fields( 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) + field_report, mismatch = field_compare_report(name, actual[name], expected[name], rtol=rtol, atol=atol, mesh_context=mesh_context) report[name] = field_report if mismatch is not None: attribution = field_comparison_attribution(mode, name, mismatch["reason"], stage_graph) @@ -756,6 +1257,209 @@ def compare_fields( return report, mismatches +def mismatch_sort_key(mismatch: Mapping[str, Any]) -> tuple[int, int, int]: + attribution = mismatch.get("attribution") if isinstance(mismatch.get("attribution"), Mapping) else {} + stage_group = attribution.get("likely_stage_group") + field = str(mismatch.get("field", "")) + try: + field_index = REQUIRED_FIELDS.index(field) + except ValueError: + field_index = len(REQUIRED_FIELDS) + return ( + MODE_COMPARISON_ORDER.get(str(mismatch.get("mode", "")), len(MODE_COMPARISON_ORDER)), + STAGE_GROUP_ORDER.get(str(stage_group), len(STAGE_GROUP_ORDER)), + field_index, + ) + + +def field_failure_context(mismatch: Mapping[str, Any]) -> dict[str, Any]: + return { + "field": mismatch.get("field"), + "reason": mismatch.get("reason"), + "shape": { + "actual": mismatch.get("actual_shape"), + "expected": mismatch.get("expected_shape"), + "matches": mismatch.get("shape_matches"), + }, + "entity_kind": { + "actual": mismatch.get("actual_entity_kind"), + "expected": mismatch.get("expected_entity_kind"), + "matches": mismatch.get("entity_kind_matches"), + }, + "tolerance": { + "rtol": mismatch.get("rtol"), + "atol": mismatch.get("atol"), + "at_largest_difference": mismatch.get("tolerance_at_max"), + }, + "error": { + "max_abs": mismatch.get("max_abs"), + "mean_abs": mismatch.get("mean_abs"), + "rms_abs": mismatch.get("rms_abs"), + "nonfinite_error_count": mismatch.get("nonfinite_error_count"), + }, + "largest_difference": { + "location": mismatch.get("location"), + "actual_value": mismatch.get("actual_at_max"), + "expected_value": mismatch.get("expected_at_max"), + "actual_entity_value": mismatch.get("actual_entity_at_max"), + "expected_entity_value": mismatch.get("expected_entity_at_max"), + }, + "local_context": mismatch.get("local_context"), + "top_differences": mismatch.get("top_differences"), + } + +def stage_group_observability(report: Mapping[str, Any], mode: Any, group_name: Any) -> dict[str, Any] | None: + observability = report.get("stage_observability", {}).get(mode, {}) + groups = observability.get("groups", []) if isinstance(observability, Mapping) else [] + for group in groups: + if isinstance(group, Mapping) and group.get("name") == group_name: + return dict(group) + return None + + +def finest_supported_target(mismatch: Mapping[str, Any], attribution: Mapping[str, Any], report: Mapping[str, Any]) -> dict[str, Any]: + mode = mismatch.get("mode") + stage_group = attribution.get("likely_stage_group", "unknown") + stages = list(attribution.get("likely_stages") or []) + observability = stage_group_observability(report, mode, stage_group) + output_targets: list[dict[str, Any]] = [] + if mismatch.get("kind") == "mesh_identity": + return { + "kind": "identity_artifact_mismatch", + "single_next_target": f"{mode}.mesh_identity", + "finest_granularity": "mesh identity digest/patch table", + "artifact_family": "mesh_identity", + "evidence_path": f"modes.{mode}.mesh_comparison", + "differences": mismatch.get("differences"), + "why": "mesh identity failed before field or solver-stage comparisons, so all downstream numerical evidence is untrustworthy", + } + if observability is not None: + output_checks = observability.get("output_checks", {}) + if isinstance(output_checks, Mapping): + for stage_name in stages: + check = output_checks.get(stage_name) + if not isinstance(check, Mapping): + continue + required = list(check.get("required") or []) + output_targets.append( + { + "stage": stage_name, + "artifact_outputs": required, + "available_outputs": list(check.get("available") or []), + "missing_outputs": list(check.get("missing") or []), + "evidence_path": f"stage_observability.{mode}.{stage_group}.{stage_name}.outputs", + } + ) + + if mismatch.get("reason") in {"missing_field", "shape_mismatch"}: + return { + "kind": "field_export_precondition", + "single_next_target": f"{mode}.{mismatch.get('field')} field export/state mapping", + "finest_granularity": "field availability/shape", + "evidence_path": attribution.get("field_evidence_path"), + "why": "the required comparison field is absent or has incompatible topology, so solver-stage attribution is not yet trustworthy", + } + + artifact_family = "solver" if stage_group == "linear_solve_results" else "matrix_operator" if stage_group in {"momentum_assembly", "pressure_assembly"} else None + artifact_comparison = report.get("artifact_comparisons", {}).get(artifact_family) if artifact_family else None + if isinstance(artifact_comparison, Mapping) and artifact_comparison.get("allclose") is False: + failed_checks = artifact_comparison.get("failed_checks", []) + first_failed = failed_checks[0] if failed_checks and isinstance(failed_checks[0], Mapping) else {} + return { + "kind": "failed_artifact_comparison", + "single_next_target": f"artifact_comparisons.{artifact_family}.{first_failed.get('name')}", + "finest_granularity": "reference-vs-candidate numeric artifact comparison", + "stage": first_failed.get("name"), + "artifact_family": artifact_family, + "evidence_path": f"artifact_comparisons.{artifact_family}", + "failed_check": first_failed, + "why": "a direct reference-vs-candidate artifact comparison now exists and failed; fix this artifact before downstream field mismatches", + } + + if output_targets: + first_stage = output_targets[0]["stage"] + return { + "kind": "missing_artifact_comparison", + "single_next_target": f"{mode}.{stage_group}.{first_stage} artifact comparison evidence", + "finest_granularity": "stage output family", + "stage": first_stage, + "artifact_family": stage_group, + "artifact_outputs": output_targets[0]["artifact_outputs"], + "evidence_path": output_targets[0]["evidence_path"], + "available_stage_outputs": output_targets, + "why": "final field mismatch is the first numerical failure; no finer reference-vs-candidate artifact comparison is recorded yet, so the missing artifact comparison is the blocking evidence gap", + } + + return { + "kind": "missing_stage_granularity", + "single_next_target": f"{mode}.{stage_group} observability/evidence", + "finest_granularity": "stage group", + "stage_names": stages, + "evidence_path": attribution.get("evidence_path"), + "why": "the harness cannot name a finer artifact from the current evidence; add or inspect finer stage artifact evidence before chasing downstream fields", + } + + + +def build_failure_diagnostics(mismatches: list[dict[str, Any]], report: Mapping[str, Any]) -> dict[str, Any]: + if not mismatches: + return { + "status": "passed", + "first_blocking_point": None, + "downstream_unreliable": [], + "basis": "all required field and mesh comparisons passed", + } + + ordered = sorted(mismatches, key=mismatch_sort_key) + first = ordered[0] + attribution = first.get("attribution") if isinstance(first.get("attribution"), Mapping) else {} + first_stage_group = attribution.get("likely_stage_group", "unknown") + first_mode = first.get("mode") + first_stage_index = STAGE_GROUP_ORDER.get(str(first_stage_group), len(STAGE_GROUP_ORDER)) + downstream = [] + for mismatch in ordered[1:]: + mismatch_attr = mismatch.get("attribution") if isinstance(mismatch.get("attribution"), Mapping) else {} + mismatch_stage = mismatch_attr.get("likely_stage_group", "unknown") + if mismatch.get("mode") != first_mode or STAGE_GROUP_ORDER.get(str(mismatch_stage), len(STAGE_GROUP_ORDER)) >= first_stage_index: + downstream.append( + { + "mode": mismatch.get("mode"), + "field": mismatch.get("field"), + "reason": mismatch.get("reason"), + "stage_group": mismatch_stage, + "why_unreliable": "blocked by the first failing comparison; fix earlier evidence before interpreting this mismatch", + } + ) + + observability_path = attribution.get("evidence_path") + observed_group = None + if observability_path: + cursor: Any = report + for part in str(observability_path).split("."): + cursor = cursor.get(part) if isinstance(cursor, Mapping) else None + if isinstance(cursor, Mapping): + observed_group = cursor + + first_target = finest_supported_target(first, attribution, report) + return { + "status": "failed", + "basis": "ordered by execution mode and solver-stage dependency before field declaration order; the first target is the finest granularity supported by recorded evidence", + "first_blocking_point": { + "mode": first_mode, + "category": COMPARISON_FAILURE, + "affected_artifact_family": "mesh_identity" if first.get("kind") == "mesh_identity" else first_target.get("artifact_family", "field_comparison"), + "stage_group": first_stage_group, + "stage_names": attribution.get("likely_stages", []), + "evidence_path": first_target.get("evidence_path") or observability_path, + "first_target": first_target, + "field_context": field_failure_context(first), + "stage_observability": observed_group, + }, + "downstream_unreliable": downstream, + "all_mismatch_count": len(mismatches), + } + + def run_openfoam_command(cmd: list[str], *, log_path: Path, timeout: int) -> dict[str, Any]: started = time.monotonic() env = openfoam_env() @@ -787,6 +1491,54 @@ def run_openfoam_command(cmd: list[str], *, log_path: Path, timeout: int) -> dic "timeout_seconds": timeout, } +def diagnostic_artifact_root(args: argparse.Namespace) -> Path: + return args.diagnostic_artifacts if args.diagnostic_artifacts is not None else args.work / "diagnostic_artifacts" + + +def run_reference_split_diagnostic(args: argparse.Namespace) -> tuple[dict[str, Any], dict[str, Any] | None]: + """Run the CPU split reference in a separate verifier process for GPU artifact diagnostics.""" + work = args.work / "reference_split_diagnostic" + report_path = work / "report.json" + reference_artifact_root = diagnostic_artifact_root(args) / "reference_split" + reference_artifact_path = reference_artifact_root / "split.npz" + cmd = [ + sys.executable, + str(Path(__file__).resolve()), + "--backend", + "cpu", + "--source", + str(args.source), + "--work", + str(work), + "--report", + str(report_path), + "--rtol", + str(args.rtol), + "--atol", + str(args.atol), + "--diagnostic-artifacts", + str(reference_artifact_root), + ] + started = time.monotonic() + env = dict(os.environ) + env.pop("PYTHONPATH", None) + completed = subprocess.run(cmd, cwd=ROOT, env=env, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True, timeout=900) + summary = { + "cmd": cmd, + "exit_code": completed.returncode, + "duration_seconds": round(time.monotonic() - started, 6), + "report": report_path, + "diagnostic_artifacts": reference_artifact_path, + "stdout_tail": completed.stdout[-4000:], + } + if not report_path.exists(): + return summary, None + try: + return summary, json.loads(report_path.read_text()) + except Exception as exc: + summary["report_parse_error"] = {"type": type(exc).__name__, "message": str(exc)} + return summary, None + def patch_start_from_latest(case: Path) -> None: path = case / "system/controlDict" @@ -845,6 +1597,7 @@ def prepare_work_cases(source: Path, work: Path, *, include_split: bool) -> dict oracle = work / "oracle_case" run_one = work / "run_one_case" split = work / "split_case" if include_split else None + reference_split = work / "reference_split_case" if include_split else None prepared_meta = prepare_case(source, prepared_case, end_time=1) manifests: dict[str, Any] = {"prepared": load_case_manifest(prepared_case)} @@ -852,17 +1605,20 @@ def prepare_work_cases(source: Path, work: Path, *, include_split: bool) -> dict 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) + manifests["reference_split"] = clone_prepared_case(prepared_case, reference_split, source, prepared_meta) return { "prepared_case": prepared_case, "oracle_case": oracle, "run_one_case": run_one, "split_case": split, + "reference_split_case": reference_split, "metadata": { "prepared": prepared_meta, "oracle": prepared_meta, "run_one": prepared_meta, "split": prepared_meta if split is not None else None, + "reference_split": prepared_meta if reference_split is not None else None, }, "manifests": manifests, "prepared_case_identity": compare_prepared_case_manifests(manifests), @@ -907,8 +1663,8 @@ def prepare_gpu_solver_inputs(foam: Any, stepper: Any, case: Path, prepared: Map 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 run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: Path, *, diagnostic_artifact_path: Path | None = None) -> dict[str, Any]: + return _call_gpu_backend("run_gpu_solver_stage_smoke", stepper, backend, case, diagnostic_artifact_path=diagnostic_artifact_path) def gpu_solver_stage_report(gpu_run: Mapping[str, Any]) -> dict[str, Any]: @@ -994,8 +1750,8 @@ def run_backend_iteration(stepper: Any, backend: Mapping[str, Any], case: Path) } -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 run_gpu_split_iteration(foam: Any, stepper: Any, backend: Mapping[str, Any], case: Path, *, diagnostic_artifact_path: Path | None = None) -> dict[str, Any]: + return _call_gpu_backend("run_gpu_split_iteration", foam, stepper, backend, case, diagnostic_artifact_path=diagnostic_artifact_path) def make_stepper(foam: Any, case: Path, label: str) -> Any: @@ -1083,7 +1839,7 @@ def checked_step(stages: list[dict[str, Any]], name: str, fn: Any) -> Any: return result -def run_split_iteration(foam: Any, stepper: Any) -> dict[str, Any]: +def run_split_iteration(foam: Any, stepper: Any, *, diagnostic_artifact_path: Path | None = None) -> dict[str, Any]: stages: list[dict[str, Any]] = [] checked_step(stages, "pre_solve", stepper.pre_solve) checked_step(stages, "advance_time", stepper.advance_time) @@ -1101,14 +1857,14 @@ def run_split_iteration(foam: Any, stepper: Any) -> dict[str, Any]: checked_step(stages, "momentum_transport_predict", stepper.momentum_transport_predictor) terms = checked_step(stages, "assemble_momentum_terms", stepper.assemble_momentum_terms) UEqn = checked_step(stages, "assemble_UEqn", stepper.assemble_momentum_matrix) - checked_step(stages, "relax_UEqn", stepper.relax_matrix) + relax_UEqn = checked_step(stages, "relax_UEqn", stepper.relax_matrix) checked_step(stages, "constrain_UEqn", stepper.constrain_matrix) - checked_step(stages, "solve_UEqn", stepper.solve_momentum) - checked_step(stages, "compute_pressure_inputs", stepper.compute_pressure_inputs) + solve_UEqn = checked_step(stages, "solve_UEqn", stepper.solve_momentum) + compute_pressure_inputs = checked_step(stages, "compute_pressure_inputs", stepper.compute_pressure_inputs) pEqn = checked_step(stages, "assemble_pEqn", stepper.assemble_pressure_matrix) - checked_step(stages, "solve_pEqn", stepper.solve_pressure) - checked_step(stages, "correct_velocity_pressure_flux", stepper.correct_velocity_pressure_flux) - checked_step(stages, "momentum_transport_correct", stepper.momentum_transport_corrector) + solve_pEqn = checked_step(stages, "solve_pEqn", stepper.solve_pressure) + correct_velocity_pressure_flux = checked_step(stages, "correct_velocity_pressure_flux", stepper.correct_velocity_pressure_flux) + momentum_transport_correct = checked_step(stages, "momentum_transport_correct", stepper.momentum_transport_corrector) checked_step(stages, "end_pimple_iteration", stepper.end_pimple_iteration) checked_step(stages, "post_solve", lambda: stepper.post_solve(write=False)) @@ -1131,6 +1887,18 @@ def run_split_iteration(foam: Any, stepper: Any) -> dict[str, Any]: "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")), } + diagnostic_artifacts = write_split_numeric_artifacts( + diagnostic_artifact_path, + role="split", + relax_UEqn=relax_UEqn, + solve_UEqn=solve_UEqn, + compute_pressure_inputs=compute_pressure_inputs, + assemble_pEqn=pEqn, + solve_pEqn=solve_pEqn, + correct_velocity_pressure_flux=correct_velocity_pressure_flux, + momentum_transport_correct=momentum_transport_correct, + fields=fields, + ) observability = stage_observability_report( "split", stages, @@ -1152,6 +1920,7 @@ def run_split_iteration(foam: Any, stepper: Any) -> dict[str, Any]: "pEqn": matrix_summary(pEqn_matrix), "matrix_states": matrix_states, "observability": observability, + "diagnostic_artifacts": diagnostic_artifacts, } @@ -1284,6 +2053,343 @@ def build_timing_evidence(report: Mapping[str, Any]) -> dict[str, Any]: }, } +def path_presence(report: Mapping[str, Any], paths: Iterable[str]) -> tuple[list[str], list[str]]: + available = [] + missing = [] + for path in paths: + (available if report_path_exists(report, path) else missing).append(path) + return available, missing + +def artifact_stats(value: Any) -> Mapping[str, Any] | None: + if not isinstance(value, Mapping): + return None + if "sha256" in value and "shape" in value: + return value + stats = value.get("stats") + if isinstance(stats, Mapping) and "sha256" in stats and "shape" in stats: + return stats + internal = value.get("internal") + if isinstance(internal, Mapping) and "sha256" in internal and "shape" in internal: + return internal + return None + + +def compare_artifact_stats(report: Mapping[str, Any], *, name: str, reference_path: str, candidate_path: str) -> dict[str, Any]: + reference = artifact_stats(report_path_value(report, reference_path)) + candidate = artifact_stats(report_path_value(report, candidate_path)) + out: dict[str, Any] = { + "name": name, + "reference_path": reference_path, + "candidate_path": candidate_path, + "reference_available": reference is not None, + "candidate_available": candidate is not None, + } + if reference is None or candidate is None: + return {**out, "shape_matches": False, "allclose": False, "reason": "missing_artifact_stats"} + shape_matches = reference.get("shape") == candidate.get("shape") + dtype_matches = reference.get("dtype") == candidate.get("dtype") + sha_matches = reference.get("sha256") == candidate.get("sha256") + return { + **out, + "shape_matches": shape_matches, + "dtype_matches": dtype_matches, + "allclose": bool(shape_matches and dtype_matches and sha_matches), + "reason": None if shape_matches and dtype_matches and sha_matches else "artifact_hash_mismatch", + "reference": {key: reference.get(key) for key in ("shape", "dtype", "size", "sha256", "min", "max", "mean")}, + "candidate": {key: candidate.get(key) for key in ("shape", "dtype", "size", "sha256", "min", "max", "mean")}, + } + + +def summarize_artifact_comparison(name: str, checks: list[dict[str, Any]]) -> dict[str, Any]: + failures = [check for check in checks if check.get("allclose") is not True] + return { + "name": name, + "shape_matches": all(check.get("shape_matches") is True for check in checks), + "allclose": not failures, + "reason": None if not failures else failures[0].get("reason"), + "checks": checks, + "failed_checks": failures, + } + + +def build_artifact_comparisons(report: Mapping[str, Any], *, mesh_context: Any | None = None) -> dict[str, Any]: + pressure_input_pairs = [ + ("HbyA", "pressure_inputs.HbyA", "modes.reference_split.stages.compute_pressure_inputs.outputs.HbyA.internal", "modes.split.stages.compute_pressure_inputs.outputs.HbyA"), + ("phiHbyA", "pressure_inputs.phiHbyA", "modes.reference_split.stages.compute_pressure_inputs.outputs.phiHbyA.internal", "modes.split.stages.compute_pressure_inputs.outputs.phiHbyA"), + ("rAU", "pressure_inputs.rAU", "modes.reference_split.stages.compute_pressure_inputs.outputs.rAU.internal", "modes.split.stages.compute_pressure_inputs.outputs.rAU"), + ("rAtU", "pressure_inputs.rAtU", "modes.reference_split.stages.compute_pressure_inputs.outputs.rAtU.internal", "modes.split.stages.compute_pressure_inputs.outputs.rAtU"), + ] + matrix_pairs = [ + ("UEqn.diag", "matrix_operator.UEqn.diag", "modes.reference_split.stages.relax_UEqn.outputs.UEqn.diag", "modes.split.UEqn.diag"), + ("UEqn.upper", "matrix_operator.UEqn.upper", "modes.reference_split.stages.relax_UEqn.outputs.UEqn.upper", "modes.split.UEqn.upper"), + ("UEqn.lower", "matrix_operator.UEqn.lower", "modes.reference_split.stages.relax_UEqn.outputs.UEqn.lower", "modes.split.UEqn.lower"), + ("UEqn.source", "matrix_operator.UEqn.source", "modes.reference_split.stages.relax_UEqn.outputs.UEqn.source", "modes.split.UEqn.source"), + ("UEqn.psi", "matrix_operator.UEqn.psi", "modes.reference_split.stages.relax_UEqn.outputs.UEqn.psi", "modes.split.UEqn.psi"), + ("pEqn.diag", "matrix_operator.pEqn.diag", "modes.reference_split.pEqn.diag", "modes.split.pEqn.diag"), + ("pEqn.upper", "matrix_operator.pEqn.upper", "modes.reference_split.pEqn.upper", "modes.split.pEqn.upper"), + ("pEqn.source", "matrix_operator.pEqn.source", "modes.reference_split.pEqn.source", "modes.split.pEqn.source"), + ("pEqn.psi", "matrix_operator.pEqn.psi", "modes.reference_split.pEqn.psi", "modes.split.pEqn.psi"), + ] + solver_pairs = [ + ("solve_UEqn.field_after", "solver.solve_UEqn.field_after", "modes.reference_split.stages.solve_UEqn.outputs.field_after", "modes.split.stages.solve_UEqn.outputs.field_after"), + ("solve_pEqn.p", "solver.solve_pEqn.p", "modes.reference_split.stages.solve_pEqn.outputs.p", "modes.split.stages.solve_pEqn.outputs.p"), + ("solve_pEqn.phi", "solver.solve_pEqn.phi", "modes.reference_split.stages.solve_pEqn.outputs.phi", "modes.split.stages.solve_pEqn.outputs.phi"), + ] + + reference_artifact_path = diagnostic_artifact_path(report, "reference_split") + candidate_artifact_path = diagnostic_artifact_path(report, "split") + reference_artifacts = load_numeric_artifact_file(reference_artifact_path) + candidate_artifacts = load_numeric_artifact_file(candidate_artifact_path) + tolerances = report.get("tolerances", {}) if isinstance(report.get("tolerances"), Mapping) else {} + rtol = float(tolerances.get("rtol", 0.0) or 0.0) + atol = float(tolerances.get("atol", 0.0) or 0.0) + backend_selected = report.get("backend", {}).get("selected") if isinstance(report.get("backend"), Mapping) else None + + + def compare_pair(name: str, key: str, reference_path: str, candidate_path: str) -> dict[str, Any]: + if backend_selected == "cpu" and reference_artifacts is None: + if candidate_artifacts is not None: + return compare_numeric_artifact_arrays( + name, + key, + candidate_artifacts.get(key), + candidate_artifacts.get(key), + rtol=rtol, + atol=atol, + mesh_context=mesh_context, + ) + return compare_artifact_stats(report, name=name, reference_path=candidate_path, candidate_path=candidate_path) + if reference_artifacts is not None and candidate_artifacts is not None: + return compare_numeric_artifact_arrays( + name, + key, + reference_artifacts.get(key), + candidate_artifacts.get(key), + rtol=rtol, + atol=atol, + mesh_context=mesh_context, + ) + return compare_artifact_stats(report, name=name, reference_path=reference_path, candidate_path=candidate_path) + + pressure_input_checks = [compare_pair(name, key, reference, candidate) for name, key, reference, candidate in pressure_input_pairs] + matrix_checks = [compare_pair(name, key, reference, candidate) for name, key, reference, candidate in matrix_pairs] + solver_checks = [compare_pair(name, key, reference, candidate) for name, key, reference, candidate in solver_pairs] + return { + "schema_version": 2, + "numeric_artifacts": { + "reference_path": str(reference_artifact_path) if reference_artifact_path is not None else None, + "candidate_path": str(candidate_artifact_path) if candidate_artifact_path is not None else None, + "reference_loaded": reference_artifacts is not None, + "candidate_loaded": candidate_artifacts is not None, + "comparison_basis": "npz_numeric_artifacts" if reference_artifacts is not None and candidate_artifacts is not None else "candidate_numeric_self_check" if backend_selected == "cpu" and candidate_artifacts is not None else "summary_hash_stats", + }, + "pressure_inputs": summarize_artifact_comparison("pressure_inputs", pressure_input_checks), + "matrix_operator": summarize_artifact_comparison("matrix_operator", matrix_checks), + "solver": summarize_artifact_comparison("solver", solver_checks), + } + + +def comparison_paths_failed(report: Mapping[str, Any], paths: Iterable[str]) -> list[dict[str, Any]]: + failed = [] + for path in paths: + value = report_path_value(report, path) + if isinstance(value, Mapping): + if value.get("matches_oracle") is False or value.get("allclose") is False or value.get("shape_matches") is False: + failed_checks = value.get("failed_checks", []) + first_failed = failed_checks[0] if failed_checks and isinstance(failed_checks[0], Mapping) else value + failed.append( + { + "path": path, + "status": "failed", + "reason": first_failed.get("reason") or value.get("reason"), + "first_failed_check": first_failed.get("name"), + "max_abs": first_failed.get("max_abs"), + "mean_abs": first_failed.get("mean_abs"), + "rms_abs": first_failed.get("rms_abs"), + } + ) + return failed + + +def family_downstream_of_first(family: Mapping[str, Any], first_blocker: Mapping[str, Any]) -> bool: + first_group = first_blocker.get("stage_group") + if not first_group or first_group == "unknown": + return False + first_index = STAGE_GROUP_ORDER.get(str(first_group)) + family_groups = family.get("stage_groups", ()) + family_indexes = [STAGE_GROUP_ORDER[group] for group in family_groups if group in STAGE_GROUP_ORDER] + return first_index is not None and bool(family_indexes) and min(family_indexes) > first_index + + +def intermediate_artifact_family_report(report: Mapping[str, Any]) -> dict[str, Any]: + diagnostics = report.get("failure_diagnostics", {}) if isinstance(report.get("failure_diagnostics"), Mapping) else {} + backend_selected = report.get("backend", {}).get("selected") + first_blocker = diagnostics.get("first_blocking_point", {}) if isinstance(diagnostics.get("first_blocking_point"), Mapping) else {} + families = [] + for family in INTERMEDIATE_ARTIFACT_FAMILIES: + reference_available, reference_missing = path_presence(report, family["reference_paths"]) + candidate_available, candidate_missing = path_presence(report, family["candidate_paths"]) + comparison_available, comparison_missing = path_presence(report, family["comparison_paths"]) + comparison_failures = comparison_paths_failed(report, comparison_available) + contract_gaps = [] + if backend_selected != "cpu": + if not family["reference_paths"]: + contract_gaps.append("reference_artifact_evidence") + if not family["comparison_paths"]: + contract_gaps.append("direct_reference_candidate_comparison") + downstream = family_downstream_of_first(family, first_blocker) + if downstream: + status = "downstream_unreliable" + elif comparison_failures: + status = "failed" + elif backend_selected == "cpu" and not candidate_missing and candidate_available and not comparison_failures: + status = "passed" + elif reference_missing or candidate_missing or comparison_missing or contract_gaps: + status = "missing" + elif comparison_available: + status = "passed" + else: + status = "available" + families.append( + { + "name": family["name"], + "required": family["required"], + "status": status, + "stage_groups": list(family["stage_groups"]), + "reference_evidence": {"available": reference_available, "missing": reference_missing}, + "candidate_evidence": {"available": candidate_available, "missing": candidate_missing}, + "comparison_evidence": {"available": comparison_available, "missing": comparison_missing, "failures": comparison_failures}, + "missing_core_evidence": contract_gaps, + "comparison_basis": "cpu_openfoam_artifact_self_check" if backend_selected == "cpu" and status == "passed" and not comparison_available else "direct_reference_candidate_comparison", + "blocked_by": first_blocker.get("first_target") if downstream else None, + "why": ( + "later artifact family is downstream of the first blocker" + if downstream + else "direct reference-vs-candidate comparison failed" + if comparison_failures + else "required reference, candidate, or direct comparison evidence is missing" + if status == "missing" + else "CPU/OpenFOAM artifact evidence is present on the trusted execution path" + if backend_selected == "cpu" and status == "passed" and not comparison_available + else "required direct comparison evidence is present and passed" + if status == "passed" + else "artifact evidence is present but has no direct comparison contract" + ), + } + ) + first_actionable = next((family for family in families if family["status"] in {"failed", "missing"}), None) + status = "complete" if all(family["status"] == "passed" for family in families if family["required"]) else "incomplete" + return { + "schema_version": 1, + "status": status, + "families": families, + "first_actionable_family": first_actionable, + "status_legend": { + "passed": "required reference, candidate, and direct comparison evidence exists and passed", + "failed": "direct comparison evidence exists and failed", + "missing": "required reference, candidate, or comparison evidence is absent", + "downstream_unreliable": "earlier blocker makes this family unsuitable for diagnosis", + "available": "artifacts are present but not directly compared", + }, + } + +def first_divergence_summary(report: Mapping[str, Any]) -> dict[str, Any]: + diagnostics = report.get("failure_diagnostics", {}) if isinstance(report.get("failure_diagnostics"), Mapping) else {} + first_blocker = diagnostics.get("first_blocking_point") if isinstance(diagnostics.get("first_blocking_point"), Mapping) else None + artifacts = report.get("intermediate_artifacts", {}) if isinstance(report.get("intermediate_artifacts"), Mapping) else {} + families = artifacts.get("families", []) if isinstance(artifacts.get("families"), list) else [] + family_statuses = { + family.get("name"): family.get("status") + for family in families + if isinstance(family, Mapping) + } + downstream_ignore = [ + { + "mode": item.get("mode"), + "field": item.get("field"), + "artifact_family": item.get("stage_group"), + "reason": item.get("why_unreliable"), + } + for item in diagnostics.get("downstream_unreliable", []) + if isinstance(item, Mapping) + ] + if first_blocker is not None: + field_context = first_blocker.get("field_context", {}) if isinstance(first_blocker.get("field_context"), Mapping) else {} + first_target = first_blocker.get("first_target", {}) if isinstance(first_blocker.get("first_target"), Mapping) else {} + return { + "schema_version": 1, + "status": "failed", + "category": first_blocker.get("category"), + "first_target": first_target.get("single_next_target"), + "target_kind": first_target.get("kind"), + "evidence_path": first_target.get("evidence_path") or first_blocker.get("evidence_path"), + "artifact_family": first_blocker.get("affected_artifact_family") or first_target.get("artifact_family"), + "mode": first_blocker.get("mode"), + "stage_group": first_blocker.get("stage_group"), + "stage_names": first_blocker.get("stage_names"), + "field": field_context.get("field"), + "reason": field_context.get("reason") or first_target.get("why"), + "local_context": field_context.get("local_context"), + "missing_evidence": first_target.get("missing_core_evidence") or first_target.get("artifact_outputs"), + "intermediate_artifact_statuses": family_statuses, + "downstream_symptoms_to_ignore": downstream_ignore, + } + + first_family = artifacts.get("first_actionable_family") if isinstance(artifacts.get("first_actionable_family"), Mapping) else None + if first_family is not None: + missing_evidence = { + "reference": first_family.get("reference_evidence", {}).get("missing"), + "candidate": first_family.get("candidate_evidence", {}).get("missing"), + "comparison": first_family.get("comparison_evidence", {}).get("missing"), + "core": first_family.get("missing_core_evidence"), + } + return { + "schema_version": 1, + "status": first_family.get("status") or "incomplete", + "category": REQUIRED_EVIDENCE_FAILURE, + "first_target": f"{first_family.get('name')} artifact evidence", + "target_kind": "intermediate_artifact_family", + "evidence_path": "intermediate_artifacts.first_actionable_family", + "artifact_family": first_family.get("name"), + "mode": None, + "stage_group": (first_family.get("stage_groups") or [None])[0], + "stage_names": first_family.get("stage_groups"), + "field": None, + "reason": first_family.get("why"), + "local_context": None, + "missing_evidence": missing_evidence, + "intermediate_artifact_statuses": family_statuses, + "downstream_symptoms_to_ignore": [ + { + "artifact_family": family.get("name"), + "reason": "ignore until the first missing or failed artifact family is closed", + } + for family in families + if isinstance(family, Mapping) and family.get("status") == "downstream_unreliable" + ], + } + + return { + "schema_version": 1, + "status": "passed", + "category": None, + "first_target": None, + "target_kind": None, + "evidence_path": None, + "artifact_family": None, + "mode": None, + "stage_group": None, + "stage_names": [], + "field": None, + "reason": "all required comparisons and artifact families passed", + "local_context": None, + "missing_evidence": None, + "intermediate_artifact_statuses": family_statuses, + "downstream_symptoms_to_ignore": [], + } + + + 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)] @@ -1305,11 +2411,77 @@ def comparison_regression_summary(comparisons: Mapping[str, Any]) -> dict[str, A return out +def backend_trust_evidence(report: Mapping[str, Any]) -> dict[str, Any]: + backend = report.get("backend", {}) if isinstance(report.get("backend"), Mapping) else {} + selected = backend.get("selected") + requested = backend.get("requested") + modes = { + name: mode + for name, mode in report.get("modes", {}).items() + if name in {"run_one", "split"} and isinstance(mode, Mapping) and mode.get("enabled") is True + } + primitive = backend.get("primitive_evidence") if isinstance(backend.get("primitive_evidence"), Mapping) else {} + capabilities = set(backend.get("capabilities", ())) if isinstance(backend.get("capabilities", ()), Iterable) else set() + + mode_gpu_ownership = {} + for name, mode in modes.items(): + mode_backend = mode.get("backend") if isinstance(mode.get("backend"), Mapping) else {} + gpu_solver = mode.get("gpu_solver") if isinstance(mode.get("gpu_solver"), Mapping) else {} + gpu_solver_backend = gpu_solver.get("backend") if isinstance(gpu_solver.get("backend"), Mapping) else {} + execution_path = str(mode.get("execution_path") or "").lower() + mode_gpu_ownership[name] = { + "backend_selected_gpu": mode_backend.get("selected") == "gpu", + "execution_path_gpu_owned": "gpu" in execution_path and "cpu" not in execution_path, + "gpu_solver_executed": gpu_solver.get("status") == "executed", + "gpu_solver_backend_selected_gpu": gpu_solver_backend.get("selected") == "gpu", + "used_cpu_fallback": mode_backend.get("used_cpu_fallback") is True or gpu_solver_backend.get("used_cpu_fallback") is True, + } + + gpu_requested = requested == "gpu" or selected == "gpu" + gpu_claim_checks = { + "selected_gpu_backend": selected == "gpu", + "provider_not_cpu": "cpu" not in str(backend.get("provider") or "").lower() + and str(backend.get("provider") or "").lower() not in {"openfoam", "foam_stepper_cpu"}, + "device_not_host": str(backend.get("device") or "").lower() not in {"", "host", "cpu"}, + "no_cpu_fallback_capability": "no_cpu_fallback" in capabilities, + "used_cpu_fallback_false": backend.get("used_cpu_fallback") is False, + "cpu_fallback_disallowed": (backend.get("cpu_fallback") if isinstance(backend.get("cpu_fallback"), Mapping) else {}).get("allowed") is False, + "primitive_not_acceptance": primitive.get("counts_as_full_gpu_rans_solver") is False, + "full_solver_guard_named": bool(backend.get("full_solver_guard")), + "modes_gpu_owned": bool(mode_gpu_ownership) and all( + item["backend_selected_gpu"] + and item["execution_path_gpu_owned"] + and item["gpu_solver_executed"] + and item["gpu_solver_backend_selected_gpu"] + and not item["used_cpu_fallback"] + for item in mode_gpu_ownership.values() + ), + } + cpu_checks = { + "selected_cpu_backend": selected == "cpu", + "used_cpu_fallback_false": backend.get("used_cpu_fallback") is False, + } + checks = gpu_claim_checks if gpu_requested else cpu_checks + return { + "status": "trusted" if all(checks.values()) else "untrusted", + "requested": requested, + "selected": selected, + "gpu_claimed": gpu_requested, + "checks": checks, + "mode_gpu_ownership": mode_gpu_ownership, + "policy": "GPU-backed success requires GPU-owned run_one and split executions, no CPU fallback, and primitive evidence marked non-acceptance.", + } + + def build_verifier_evidence(report: Mapping[str, Any]) -> dict[str, Any]: - modes = enabled_mode_names(report) + modes = REQUIRED_EXECUTION_MODES 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()) + mode_enabled = { + mode: report.get("modes", {}).get(mode, {}).get("enabled") is True + for mode in REQUIRED_EXECUTION_MODES + } mesh_ok_by_mode = { mode: bool(report.get("modes", {}).get(mode, {}).get("mesh_comparison", {}).get("matches_oracle")) for mode in modes @@ -1328,7 +2500,7 @@ def build_verifier_evidence(report: Mapping[str, Any]) -> dict[str, Any]: 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, + "enabled": mode_enabled[mode], "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")], @@ -1345,15 +2517,21 @@ def build_verifier_evidence(report: Mapping[str, Any]) -> dict[str, Any]: ) else: timing_ok = timing.get("status") in {"reported", "not_reported"} + backend_trust = backend_trust_evidence(report) + intermediate_artifacts = report.get("intermediate_artifacts", {}) if isinstance(report.get("intermediate_artifacts"), Mapping) else {} + + criteria = { + "required_execution_modes": all(mode_enabled.values()), "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, + "mesh_identity": all(mesh_ok_by_mode.values()), + "field_comparisons": all(comparison_ok_by_mode.values()), + "stage_observability": all(observability_ok_by_mode.values()), + "backend_selected": report.get("backend", {}).get("selected") is not None and backend_trust.get("status") == "trusted", "timing_reported": timing_ok, - "differentiability_reported": report.get("differentiability", {}).get("status") in {"differentiable", "not_differentiable"}, + "intermediate_artifacts_complete": intermediate_artifacts.get("status") == "complete", + "intermediate_artifacts_reported": bool(intermediate_artifacts.get("families")), } passed = all(criteria.values()) oracle_mesh = report.get("mesh_identity", {}).get("oracle", {}) @@ -1366,6 +2544,20 @@ def build_verifier_evidence(report: Mapping[str, Any]) -> dict[str, Any]: "mesh_identity": mesh_ok_by_mode, "field_comparisons": comparison_ok_by_mode, "stage_observability": observability_ok_by_mode, + "backend_trust": backend_trust, + "intermediate_artifacts": { + "status": intermediate_artifacts.get("status"), + "families": { + family.get("name"): family.get("status") + for family in intermediate_artifacts.get("families", []) + if isinstance(family, Mapping) + }, + "first_actionable_family": ( + intermediate_artifacts.get("first_actionable_family", {}).get("name") + if isinstance(intermediate_artifacts.get("first_actionable_family"), Mapping) + else None + ), + }, }, "entrypoint": report.get("harness_entrypoint", {}), "regression_summary": { @@ -1373,8 +2565,10 @@ def build_verifier_evidence(report: Mapping[str, Any]) -> dict[str, Any]: "source_case": report.get("problem_boundary", {}).get("airfrans_simulation"), "prepared_case_files_sha256": prepared_digest, "backend_selected": report.get("backend", {}).get("selected"), + "backend_trust_status": backend_trust.get("status"), "differentiability_status": report.get("differentiability", {}).get("status"), "mesh_topology_sha256": oracle_mesh.get("topology_sha256"), + "first_divergence_summary": report.get("first_divergence_summary"), "mesh_geometry_sha256": oracle_mesh.get("geometry_sha256"), "modes": regression_modes, }, @@ -1492,6 +2686,11 @@ def base_report(args: argparse.Namespace) -> dict[str, Any]: "tracked_turbulence_fields": {}, "comparison_mismatches": [], "comparison_attribution": [], + "failure_diagnostics": { + "status": "not_evaluated", + "first_blocking_point": None, + "downstream_unreliable": [], + }, "verifier_evidence": { "passed": False, "status_basis": "not_evaluated", @@ -1504,6 +2703,11 @@ def base_report(args: argparse.Namespace) -> dict[str, Any]: "reason": "backend has not been selected yet", }, "state_exports": {}, + "diagnostic_artifacts": { + "schema_version": NUMERIC_ARTIFACT_SCHEMA_VERSION, + "status": "not_evaluated", + "root": diagnostic_artifact_root(args), + }, "gpu_solver_inputs": { "schema_version": GPU_INPUT_SCHEMA_VERSION, "status": "not_requested", @@ -1551,6 +2755,24 @@ def run_harness(args: argparse.Namespace, report: dict[str, Any]) -> None: 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" and split_case is not None: + reference_summary, reference_report = run_reference_split_diagnostic(args) + report["reference_split_diagnostic"] = json_ready(reference_summary) + if reference_report is not None: + reference_split = reference_report.get("modes", {}).get("split", {}) + report["modes"]["reference_split"] = { + **reference_split, + "role": "separate_process_openfoam_split_reference_for_gpu_artifact_comparison", + "source_report": reference_summary["report"], + } + reference_state_exports = reference_report.get("state_exports", {}) + if "split" in reference_state_exports: + report["state_exports"]["reference_split"] = reference_state_exports["split"] + if "split_matrices" in reference_state_exports: + report["state_exports"]["reference_split_matrices"] = reference_state_exports["split_matrices"] + reference_artifacts = reference_report.get("diagnostic_artifacts", {}) + if isinstance(reference_artifacts, Mapping) and isinstance(reference_artifacts.get("split"), Mapping): + report["diagnostic_artifacts"]["reference_split"] = reference_artifacts["split"] if backend.get("selected") == "gpu": try: foam = import_foam() @@ -1621,6 +2843,7 @@ def run_harness(args: argparse.Namespace, report: dict[str, Any]) -> None: 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") report["mesh_identity"]["run_one"] = run_one_mesh @@ -1662,6 +2885,7 @@ def run_harness(args: argparse.Namespace, report: dict[str, Any]) -> None: rtol=args.rtol, atol=args.atol, stage_graph=[stage["name"] for stage in run_one_stages], + mesh_context=run_one_stepper.mesh(), ) report["modes"]["run_one"].update( { @@ -1684,6 +2908,8 @@ def run_harness(args: argparse.Namespace, report: dict[str, Any]) -> None: mismatches = mesh_mismatches + run_one_mismatches + split_artifact_mesh = None + if split_case is not None: split_stepper = make_stepper(foam, split_case, "split") split_mesh = read_mesh_identity(split_stepper, "split") @@ -1691,10 +2917,12 @@ 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_artifact_path = diagnostic_artifact_root(args) / "split.npz" + split_artifact_mesh = split_stepper.mesh() if backend.get("selected") == "gpu": - split_result = run_gpu_split_iteration(foam, split_stepper, backend, split_case) + split_result = run_gpu_split_iteration(foam, split_stepper, backend, split_case, diagnostic_artifact_path=split_artifact_path) else: - split_result = run_split_iteration(foam, split_stepper) + split_result = run_split_iteration(foam, split_stepper, diagnostic_artifact_path=split_artifact_path) report["stage_observability"]["split"] = split_result["observability"] split_fields = split_result["fields"] split_comparisons, split_mismatches = compare_fields( @@ -1704,6 +2932,7 @@ def run_harness(args: argparse.Namespace, report: dict[str, Any]) -> None: rtol=args.rtol, atol=args.atol, stage_graph=split_result["graph"], + mesh_context=split_stepper.mesh(), ) report["modes"]["split"].update( { @@ -1724,13 +2953,20 @@ def run_harness(args: argparse.Namespace, report: dict[str, Any]) -> None: 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"] + if isinstance(split_result.get("diagnostic_artifacts"), Mapping): + report["diagnostic_artifacts"]["split"] = split_result["diagnostic_artifacts"] + report["diagnostic_artifacts"]["status"] = "available" 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["artifact_comparisons"] = json_ready(build_artifact_comparisons(report, mesh_context=split_artifact_mesh)) + report["failure_diagnostics"] = json_ready(build_failure_diagnostics(mismatches, report)) report["timing"] = json_ready(build_timing_evidence(report)) + report["intermediate_artifacts"] = json_ready(intermediate_artifact_family_report(report)) + report["first_divergence_summary"] = json_ready(first_divergence_summary(report)) report["verifier_evidence"] = json_ready(build_verifier_evidence(report)) if mismatches: raise HarnessError( @@ -1741,9 +2977,19 @@ def run_harness(args: argparse.Namespace, report: dict[str, Any]) -> None: "failure_kind": GPU_NUMERICAL_MISMATCH if backend.get("selected") == "gpu" else COMPARISON_FAILURE, "mismatches": mismatches, "likely_causes": comparison_attribution, + "first_divergence_summary": report["first_divergence_summary"], + "diagnostics": report["failure_diagnostics"], }, ) + if report["verifier_evidence"].get("passed") is not True: + raise HarnessError( + REQUIRED_EVIDENCE_FAILURE, + "verifier_evidence", + "required verifier evidence is incomplete or failed", + details={"verifier_evidence": report["verifier_evidence"], "first_divergence_summary": report["first_divergence_summary"]}, + ) + def write_report(report: Mapping[str, Any], path: Path) -> None: path.parent.mkdir(parents=True, exist_ok=True) @@ -1790,6 +3036,20 @@ def print_human_summary(report: Mapping[str, Any]) -> None: if evidence: print(f"evidence_passed={evidence.get('passed')} evidence_status_basis={evidence.get('status_basis')}") print(f"evidence_criteria={evidence.get('criteria')}") + artifacts = report.get("intermediate_artifacts") or {} + if artifacts: + families = artifacts.get("families") or [] + print(f"intermediate_artifacts_status={artifacts.get('status')}") + print(f"intermediate_artifact_families={[(family.get('name'), family.get('status')) for family in families]}") + first_family = artifacts.get("first_actionable_family") or {} + print(f"intermediate_first_actionable_family={first_family.get('name')}") + compact = report.get("first_divergence_summary") or {} + if compact: + print(f"first_divergence_summary={compact}") + print(f"first_divergence_status={compact.get('status')}") + print(f"first_divergence_target={compact.get('first_target')}") + print(f"first_divergence_evidence_path={compact.get('evidence_path')}") + print(f"first_divergence_artifact_family={compact.get('artifact_family')}") if status != "passed": failure = report.get("failure") or {} @@ -1805,6 +3065,27 @@ def print_human_summary(report: Mapping[str, Any]) -> None: 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')}") + diagnostics = ((failure.get("details") or {}).get("diagnostics") or report.get("failure_diagnostics") or {}) + first_blocker = diagnostics.get("first_blocking_point") or {} + if first_blocker: + field_context = first_blocker.get("field_context") or {} + error_context = field_context.get("error") or {} + largest = field_context.get("largest_difference") or {} + print(f"first_blocking_mode={first_blocker.get('mode')}") + print(f"first_blocking_stage_group={first_blocker.get('stage_group')}") + print(f"first_blocking_stages={first_blocker.get('stage_names')}") + first_target = first_blocker.get("first_target") or {} + print(f"first_debug_target={first_target.get('single_next_target')}") + print(f"first_debug_target_kind={first_target.get('kind')}") + print(f"first_debug_target_evidence_path={first_target.get('evidence_path')}") + print(f"first_debug_target_why={first_target.get('why')}") + print(f"first_blocking_field={field_context.get('field')}") + print(f"first_blocking_reason={field_context.get('reason')}") + print(f"first_blocking_max_abs={fmt_sci(error_context.get('max_abs'))}") + print(f"first_blocking_location={format_location(largest.get('location'))}") + local_context = field_context.get("local_context") or {} + print(f"first_blocking_local_context={local_context}") + print(f"downstream_unreliable_count={len(diagnostics.get('downstream_unreliable') or [])}") mismatches = (failure.get("details") or {}).get("mismatches") or report.get("comparison_mismatches") or [] if mismatches: first = mismatches[0] @@ -1868,6 +3149,7 @@ def parse_args(argv: list[str] | None = None) -> argparse.Namespace: parser.add_argument("--source", type=Path, default=DEFAULT_SOURCE) parser.add_argument("--work", type=Path, default=DEFAULT_WORK) parser.add_argument("--report", type=Path, default=None, help="JSON report path; defaults to WORK/verifier_report.json") + parser.add_argument("--diagnostic-artifacts", type=Path, default=None, help="Directory for focused numeric split-stage NPZ artifacts; defaults to WORK/diagnostic_artifacts") 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") diff --git a/scripts/verify_gpu_rans_solver.sh b/scripts/verify_gpu_rans_solver.sh index c70f576..0c12a4c 100755 --- a/scripts/verify_gpu_rans_solver.sh +++ b/scripts/verify_gpu_rans_solver.sh @@ -96,17 +96,23 @@ def require(condition: bool, message: str) -> None: 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 {} +backend_trust = ((verifier_evidence.get("mode_criteria") or {}).get("backend_trust") or {}) if isinstance(verifier_evidence.get("mode_criteria"), 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() - +capabilities = set(backend.get("capabilities") or []) +primitive = backend.get("primitive_evidence") if isinstance(backend.get("primitive_evidence"), dict) else {} 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(backend.get("used_cpu_fallback") is False, f"backend.used_cpu_fallback is {backend.get('used_cpu_fallback')!r}, expected false") +require((backend.get("cpu_fallback") or {}).get("allowed") is False, "backend.cpu_fallback.allowed is not false") +require("no_cpu_fallback" in capabilities, "backend capabilities do not include no_cpu_fallback") +require(primitive.get("counts_as_full_gpu_rans_solver") is False, "primitive GPU evidence is not explicitly marked non-acceptance") +require(backend_trust.get("status") == "trusted", f"backend trust status is {backend_trust.get('status')!r}, expected 'trusted'") require(verifier_evidence.get("passed") is True, "verifier_evidence.passed is not true") for mode_name in required_modes: @@ -116,6 +122,12 @@ for mode_name in required_modes: 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}") + require(mode_backend.get("used_cpu_fallback") is False, f"modes.{mode_name}.backend.used_cpu_fallback is not false") + gpu_solver = mode.get("gpu_solver") if isinstance(mode.get("gpu_solver"), dict) else {} + gpu_solver_backend = gpu_solver.get("backend") if isinstance(gpu_solver.get("backend"), dict) else {} + require(gpu_solver.get("status") == "executed", f"modes.{mode_name}.gpu_solver.status is not executed") + require(gpu_solver_backend.get("selected") == "gpu", f"modes.{mode_name}.gpu_solver.backend.selected is not 'gpu'") + require(gpu_solver_backend.get("used_cpu_fallback") is False, f"modes.{mode_name}.gpu_solver.backend.used_cpu_fallback is not false") 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 {}