openFOAM-RANS-to-GPU/python/src/foam_stepper/gpu/backend.py

"""Reusable GPU backend for the AirfRANS RANS solver verifier.

The verifier under ``scripts/`` is an acceptance harness. This module owns the
GPU backend contract, GPU-side state transfer, CUDA solver-stage graph, and
backend diagnostics used by ``--backend gpu``.
"""

from __future__ import annotations

import dataclasses
import hashlib
import math
import os
import subprocess
import time
from collections.abc import Iterable, Mapping
from pathlib import Path
from typing import Any

import numpy as np
import quadrants as qd

from .constants import (
    DEFAULT_LAMINAR_NU,
    DEFAULT_MOMENTUM_PBICGSTAB_ITERATIONS,
    DEFAULT_MOMENTUM_PBICGSTAB_MIN_ITERATIONS,
    DEFAULT_MOMENTUM_PBICGSTAB_RESIDUAL_TOLERANCE,
    DEFAULT_MOMENTUM_PBICGSTAB_RESIDUAL_TOLERANCE_SQUARED,
    DEFAULT_PRESSURE_CG_ITERATIONS,
    DEFAULT_PRESSURE_CG_MIN_ITERATIONS,
    DEFAULT_PRESSURE_CG_RESIDUAL_TOLERANCE,
    DEFAULT_PRESSURE_NON_ORTHOGONAL_CORRECTORS,
    DEFAULT_MOMENTUM_RELAXATION_ALPHA,
    DEFAULT_SIMPLE_CONSISTENT_RATU_FACTOR,
    DEFAULT_OMEGA_WALL_BETA1,
    FULL_GPU_RANS_GUARD,
    GPU_BACKEND_PROVIDER,
    GPU_INPUT_SCHEMA_VERSION,
    GPU_KERNELS_MISSING,
    GPU_NUMERICAL_MISMATCH,
    GPU_PRIMITIVE_NAME,
    GPU_PRIMITIVE_PROOF,
    GPU_RUNTIME_UNAVAILABLE,
    GPU_STAGE_CAPABILITY_PREFIX,
    GPU_STAGE_KERNEL_ENTRYPOINTS,
    GPU_STAGE_UNSUPPORTED,
    REQUIRED_FIELDS,
    STAGE_OBSERVABILITY_GROUPS,
    TURBULENCE_FIELDS,
)
from .kernels import (
    gpu_rans_face_flux_copy,
    gpu_rans_final_correction,
    gpu_rans_momentum_assembly,
    gpu_rans_momentum_diffusion_coefficients,
    gpu_rans_momentum_diffusion_coefficients_from_face_coeff,
    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_add_vector_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_internal_dev_tau_source,
    gpu_rans_momentum_wall_dev_tau_source,
    gpu_rans_momentum_linear_upwind_source,
    gpu_rans_zero_scalar_field,
    gpu_rans_momentum_offdiag_abs_accumulate,
    gpu_rans_momentum_offdiag_abs_accumulate_by_cell,
    gpu_rans_momentum_diag_from_offdiag_by_cell,
    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_by_cell,
    gpu_rans_momentum_h1_finish,
    gpu_rans_momentum_hbyA_by_cell,
    gpu_rans_momentum_hbyA_face_accumulate,
    gpu_rans_momentum_hbyA_finish,
    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_solver_diag,
    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,
    build_ldu_level_schedule,
    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_components,
    gpu_ldu_pcg_scalar_symmetric_faces,
    gpu_zero_scalar_accumulator,
)

BACKEND_CHOICES = ("auto", "cpu", "gpu")
BACKEND_FAILURE = "backend_execution_failure"
COMPARISON_FAILURE = "numerical_comparison_failure"

FIELD_COMPARISON_ATTRIBUTION = {
    "U": ("final_correction", "linear_solve_results", "momentum_assembly", "turbulence_updates"),
    "p": ("linear_solve_results", "pressure_assembly", "final_correction"),
    "phi": ("final_correction", "linear_solve_results", "pressure_assembly"),
    "nut": ("turbulence_updates",),
    "k": ("turbulence_updates",),
    "omega": ("turbulence_updates",),
}


class BackendExecutionError(Exception):
    """Backend selection or execution failure before numerical comparison."""

    def __init__(self, step: str, message: str, *, details: Mapping[str, Any] | None = None) -> None:
        super().__init__(message)
        self.step = step
        self.details = dict(details or {})


@dataclasses.dataclass(frozen=True)
class GpuSourceLocation:
    file: str
    function: str
    lines: tuple[int, int] | None = None


@dataclasses.dataclass(frozen=True)
class GpuTransformResult:
    name: str
    phase: str
    source: GpuSourceLocation
    inputs: Mapping[str, Any]
    outputs: Mapping[str, Any]
    changed_fields: list[str]
    metadata: Mapping[str, Any]


@dataclasses.dataclass(frozen=True)
class GpuPatchField:
    name: str
    type: str
    values: np.ndarray
    fixes_value: bool
    assignable: bool
    coupled: bool
    updated: bool
    patch_internal: np.ndarray | None
    value_internal_coeffs: np.ndarray | None
    value_boundary_coeffs: np.ndarray | None
    gradient_internal_coeffs: np.ndarray | None
    gradient_boundary_coeffs: np.ndarray | None


@dataclasses.dataclass(frozen=True)
class GpuField:
    name: str
    kind: str
    dimensions: str
    entity_kind: str
    entity_count: int
    internal: np.ndarray
    boundary: Mapping[str, GpuPatchField]


def json_ready(value: Any) -> Any:
    """Convert report values into strict JSON-compatible data."""

    if dataclasses.is_dataclass(value) and not isinstance(value, type):
        return json_ready(dataclasses.asdict(value))
    if isinstance(value, Path):
        return str(value)
    if isinstance(value, np.ndarray):
        return array_stats(value)
    if isinstance(value, np.generic):
        return json_ready(value.item())
    if isinstance(value, float):
        return value if math.isfinite(value) else None
    if isinstance(value, (str, int, bool)) or value is None:
        return value
    if isinstance(value, Mapping):
        return {str(key): json_ready(item) for key, item in value.items()}
    if isinstance(value, (list, tuple, set)):
        return [json_ready(item) for item in value]
    return repr(value)


def array_shape(array: Any | None) -> list[int] | None:
    if array is None:
        return None
    return [int(dim) for dim in np.asarray(array).shape]


def finite_float(value: Any) -> float | None:
    try:
        number = float(value)
    except (TypeError, ValueError):
        return None
    return number if math.isfinite(number) else None


def array_stats(array: Any) -> dict[str, Any]:
    arr = np.asarray(array)
    out: dict[str, Any] = {
        "shape": array_shape(arr),
        "dtype": str(arr.dtype),
        "size": int(arr.size),
    }
    if arr.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

    finite = np.isfinite(arr)
    out["finite_count"] = int(np.count_nonzero(finite))
    out["nonfinite_count"] = int(arr.size - out["finite_count"])
    if out["finite_count"]:
        finite_values = arr[finite]
        out.update(
            {
                "min": finite_float(np.min(finite_values)),
                "max": finite_float(np.max(finite_values)),
                "mean": finite_float(np.mean(finite_values)),
            }
        )
    return out


def value_at(array: np.ndarray, index: tuple[int, ...]) -> Any:
    value = array[index]
    if isinstance(value, np.generic):
        return json_ready(value.item())
    if isinstance(value, np.ndarray):
        return json_ready(value.tolist())
    return json_ready(value)


def entity_value_at(array: np.ndarray, entity_index: int | None) -> Any:
    if entity_index is None or array.ndim == 0:
        return None
    value = array[entity_index]
    if isinstance(value, np.generic):
        return json_ready(value.item())
    if isinstance(value, np.ndarray):
        return json_ready(value.tolist())
    return json_ready(value)


def 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"))
    digest.update(encoded)


def update_hash_array(digest: "hashlib._Hash", label: str, array: Any) -> None:
    arr = np.ascontiguousarray(np.asarray(array))
    update_hash_text(digest, label)
    update_hash_text(digest, str(arr.dtype))
    update_hash_text(digest, repr(tuple(int(dim) for dim in arr.shape)))
    digest.update(arr.tobytes())


def source_summary(source: Any) -> dict[str, Any]:
    return {
        "file": getattr(source, "file", ""),
        "function": getattr(source, "function", ""),
        "lines": json_ready(getattr(source, "lines", None)),
    }


def boundary_field_summary(patch: Any) -> dict[str, Any]:
    return {
        "type": getattr(patch, "type", ""),
        "values_shape": array_shape(getattr(patch, "values", None)),
        "fixes_value": bool(getattr(patch, "fixes_value", False)),
        "assignable": bool(getattr(patch, "assignable", False)),
        "coupled": bool(getattr(patch, "coupled", False)),
        "updated": bool(getattr(patch, "updated", False)),
        "patch_internal_shape": array_shape(getattr(patch, "patch_internal", None)),
        "value_internal_coeffs_shape": array_shape(getattr(patch, "value_internal_coeffs", None)),
        "value_boundary_coeffs_shape": array_shape(getattr(patch, "value_boundary_coeffs", None)),
        "gradient_internal_coeffs_shape": array_shape(getattr(patch, "gradient_internal_coeffs", None)),
        "gradient_boundary_coeffs_shape": array_shape(getattr(patch, "gradient_boundary_coeffs", None)),
    }


def field_summary(field: Any) -> dict[str, Any]:
    boundary = getattr(field, "boundary", {})
    return {
        "name": getattr(field, "name", ""),
        "kind": getattr(field, "kind", ""),
        "dimensions": getattr(field, "dimensions", ""),
        "entity_kind": getattr(field, "entity_kind", ""),
        "entity_count": int(getattr(field, "entity_count", 0)),
        "internal": array_stats(getattr(field, "internal")),
        "boundary": {name: boundary_field_summary(patch) for name, patch in boundary.items()},
    }


def solve_summary(performance: Any) -> dict[str, Any]:
    return {
        "solver_name": getattr(performance, "solver_name", ""),
        "field_name": getattr(performance, "field_name", ""),
        "initial_residual": json_ready(getattr(performance, "initial_residual", None)),
        "final_residual": json_ready(getattr(performance, "final_residual", None)),
        "n_iterations": json_ready(getattr(performance, "n_iterations", None)),
        "converged": bool(getattr(performance, "converged", False)),
        "singular": bool(getattr(performance, "singular", False)),
    }


def matrix_summary(matrix: Any) -> dict[str, Any]:
    derived: dict[str, Any] = {}
    for name in ("A", "H", "H1", "flux", "face_flux_correction"):
        value = getattr(matrix, name, None)
        if callable(value):
            value = value()
        if value is not None:
            derived[name] = field_summary(value)

    for name in ("residual", "D", "DD"):
        value = getattr(matrix, name, None)
        if callable(value):
            value = value()
        if value is not None:
            derived[name] = array_stats(value)

    return {
        "name": getattr(matrix, "name", ""),
        "field_name": getattr(matrix, "field_name", ""),
        "value_rank": getattr(matrix, "value_rank", ""),
        "dimensions": getattr(matrix, "dimensions", ""),
        "has_diag": bool(getattr(matrix, "has_diag", False)),
        "has_upper": bool(getattr(matrix, "has_upper", False)),
        "has_lower": bool(getattr(matrix, "has_lower", False)),
        "diagonal": bool(getattr(matrix, "diagonal", False)),
        "symmetric": bool(getattr(matrix, "symmetric", False)),
        "asymmetric": bool(getattr(matrix, "asymmetric", False)),
        "diag": array_stats(getattr(matrix, "diag")),
        "upper": None if getattr(matrix, "upper", None) is None else array_stats(getattr(matrix, "upper")),
        "lower": None if getattr(matrix, "lower", None) is None else array_stats(getattr(matrix, "lower")),
        "source": array_stats(getattr(matrix, "source")),
        "psi": field_summary(getattr(matrix, "psi")),
        "internal_coeff_shapes": [array_shape(item) for item in getattr(matrix, "internal_coeffs", [])],
        "boundary_coeff_shapes": [array_shape(item) for item in getattr(matrix, "boundary_coeffs", [])],
        "derived": derived,
    }


def summarize_value(value: Any) -> Any:
    if hasattr(value, "diag") and hasattr(value, "field_name") and hasattr(value, "source"):
        return matrix_summary(value)
    if hasattr(value, "internal") and hasattr(value, "entity_kind") and hasattr(value, "boundary"):
        return field_summary(value)
    if hasattr(value, "solver_name") and hasattr(value, "initial_residual"):
        return solve_summary(value)
    if hasattr(value, "name") and hasattr(value, "phase") and hasattr(value, "outputs"):
        return transform_summary(value)
    if isinstance(value, Mapping):
        return {str(key): summarize_value(item) for key, item in value.items()}
    if isinstance(value, list):
        return [summarize_value(item) for item in value]
    if isinstance(value, tuple):
        return [summarize_value(item) for item in value]
    return json_ready(value)


def transform_summary(result: Any) -> dict[str, Any]:
    return {
        "name": getattr(result, "name", ""),
        "phase": getattr(result, "phase", ""),
        "changed_fields": json_ready(getattr(result, "changed_fields", [])),
        "source": source_summary(getattr(result, "source", None)),
        "metadata": json_ready(getattr(result, "metadata", {})),
        "outputs": summarize_value(getattr(result, "outputs", {})),
    }

def graph_stage_summaries(result: Any) -> list[dict[str, Any]]:
    return [transform_summary(entry) for entry in getattr(result, "outputs", {}).get("graph", [])]


def stage_observability_report(mode: str, stages: list[dict[str, Any]], *, evidence: Mapping[str, Any]) -> dict[str, Any]:
    stage_by_name = {stage["name"]: stage for stage in stages}
    failures: list[dict[str, Any]] = []
    groups = []
    for group in STAGE_OBSERVABILITY_GROUPS:
        group_name = group["name"]
        required_stages = tuple(group[f"{mode}_stages"])
        missing_stages = [name for name in required_stages if name not in stage_by_name]
        output_requirements = group.get(f"{mode}_outputs", {})
        output_checks = {}
        for stage_name, required_outputs in output_requirements.items():
            stage = stage_by_name.get(stage_name)
            if stage is None:
                continue
            outputs = stage.get("outputs", {})
            output_keys = set(outputs) if isinstance(outputs, Mapping) else set()
            missing_outputs = [name for name in required_outputs if name not in output_keys]
            output_checks[stage_name] = {
                "required": list(required_outputs),
                "available": sorted(output_keys),
                "missing": missing_outputs,
            }
            if missing_outputs:
                failures.append(
                    {
                        "path": f"stage_observability.{mode}.{group_name}.{stage_name}.outputs",
                        "message": "missing inspectable stage output",
                        "expected": list(required_outputs),
                        "actual": sorted(output_keys),
                    }
                )
        if missing_stages:
            failures.append(
                {
                    "path": f"stage_observability.{mode}.{group_name}.stages",
                    "message": "missing required solver stage",
                    "expected": list(required_stages),
                    "actual": list(stage_by_name),
                }
            )
        groups.append(
            {
                "name": group_name,
                "required_stages": list(required_stages),
                "observed": not missing_stages and not any(check["missing"] for check in output_checks.values()),
                "missing_stages": missing_stages,
                "output_checks": output_checks,
            }
        )
    if failures:
        raise BackendExecutionError(
            f"{mode}_stage_observability",
            f"{mode} GPU solver-stage observability is incomplete",
            details={"failure_kind": GPU_STAGE_UNSUPPORTED, "failures": failures},
        )
    return {
        "mode": mode,
        "graph": [stage["name"] for stage in stages],
        "groups": groups,
        "evidence": json_ready(evidence),
    }


def field_dict_to_mapping(fields: Any) -> dict[str, Any]:
    if isinstance(fields, Mapping):
        return dict(fields)
    if hasattr(fields, "fields") and isinstance(fields.fields, Mapping):
        return dict(fields.fields)
    return {name: getattr(fields, name) for name in REQUIRED_FIELDS if hasattr(fields, name)}


def validation_details(exc: BaseException) -> dict[str, Any]:
    if hasattr(exc, "to_dict"):
        return json_ready(exc.to_dict())
    return {"type": type(exc).__name__, "message": str(exc)}


def read_fields(stepper: Any, label: str) -> dict[str, Any]:
    try:
        return field_dict_to_mapping(stepper.fields())
    except Exception as exc:
        raise BackendExecutionError(
            f"read_{label}_fields",
            f"failed to read {label} fields for GPU backend",
            details={"case": getattr(stepper, "case_path", ""), "cause": validation_details(exc)},
        ) from exc


def export_solver_state_checked(foam: Any, stepper: Any, label: str, fields: Mapping[str, Any] | None = None) -> dict[str, Any]:
    try:
        field_source = fields if fields is not None else stepper.fields()
        return foam.export_solver_state(stepper.mesh(), field_source, required_fields=REQUIRED_FIELDS)
    except Exception as exc:
        raise BackendExecutionError(
            f"export_{label}_state",
            f"failed to export explicit {label} solver state for GPU backend",
            details=validation_details(exc),
        ) from exc


def visible_turbulence_fields(fields: Mapping[str, Any]) -> list[str]:
    return [name for name in TURBULENCE_FIELDS if name in fields]

def gpu_device_present() -> bool:
    return any(Path(path).exists() for path in ("/dev/nvidia0", "/dev/dri/renderD128"))


def nvidia_device_identity() -> dict[str, Any]:
    try:
        completed = subprocess.run(
            ["nvidia-smi", "--query-gpu=name,uuid", "--format=csv,noheader"],
            check=True,
            stdout=subprocess.PIPE,
            stderr=subprocess.PIPE,
            text=True,
            timeout=10,
        )
    except (FileNotFoundError, subprocess.CalledProcessError, subprocess.TimeoutExpired):
        return {"device_name": None, "device_uuid": None}

    first = completed.stdout.strip().splitlines()[0] if completed.stdout.strip() else ""
    if not first:
        return {"device_name": None, "device_uuid": None}
    parts = [part.strip() for part in first.split(",", 1)]
    return {"device_name": parts[0], "device_uuid": parts[1] if len(parts) > 1 else None}


def gpu_stage_contracts() -> dict[str, dict[str, Any]]:
    contracts: dict[str, dict[str, Any]] = {}
    for group in STAGE_OBSERVABILITY_GROUPS:
        name = str(group["name"])
        split_outputs = {
            str(stage): list(outputs)
            for stage, outputs in group.get("split_outputs", {}).items()
        }
        field_coverage = sorted(
            {
                field
                for outputs in split_outputs.values()
                for field in outputs
                if field in REQUIRED_FIELDS
            }
        )
        kernel_entrypoints = list(GPU_STAGE_KERNEL_ENTRYPOINTS.get(name, ()))
        contracts[name] = {
            "required": True,
            "supported": True,
            "status": "kernels_registered" if kernel_entrypoints else "kernels_missing",
            "run_one_stages": list(group.get("run_one_stages", ())),
            "split_stages": list(group.get("split_stages", ())),
            "expected_outputs": split_outputs,
            "field_coverage": field_coverage,
            "kernel_entrypoints": kernel_entrypoints,
            "failure_kind_if_missing": None if kernel_entrypoints else GPU_KERNELS_MISSING,
        }
    return contracts


def gpu_backend_failure_taxonomy() -> dict[str, Any]:
    return {
        GPU_RUNTIME_UNAVAILABLE: {
            "category": BACKEND_FAILURE,
            "step": "select_backend",
            "meaning": "CUDA/Quadrants could not provide a non-host GPU runtime.",
        },
        GPU_KERNELS_MISSING: {
            "category": BACKEND_FAILURE,
            "step": "execute_gpu_solver_contract",
            "meaning": "The GPU backend was selected, but required solver-stage kernels are not registered.",
        },
        GPU_STAGE_UNSUPPORTED: {
            "category": BACKEND_FAILURE,
            "step": "execute_gpu_solver_contract",
            "meaning": "The GPU backend explicitly cannot execute one or more required RANS solver stages.",
        },
        GPU_NUMERICAL_MISMATCH: {
            "category": COMPARISON_FAILURE,
            "step": "compare_fields",
            "meaning": "The GPU backend executed, but one or more required fields failed OpenFOAM oracle parity.",
            "field_attribution": FIELD_COMPARISON_ATTRIBUTION,
        },
    }

def solver_state_arrays(value: Any, path: str = "") -> Iterable[tuple[str, np.ndarray]]:
    if isinstance(value, np.ndarray):
        yield path, value
        return
    if isinstance(value, Mapping):
        for key, item in value.items():
            item_path = f"{path}.{key}" if path else str(key)
            yield from solver_state_arrays(item, item_path)
        return
    if isinstance(value, (list, tuple)):
        for index, item in enumerate(value):
            yield from solver_state_arrays(item, f"{path}[{index}]")


def gpu_transfer_array_dtype(qd: Any, source: np.ndarray, path: str) -> tuple[Any, str, np.ndarray]:
    if np.issubdtype(source.dtype, np.integer):
        if source.size:
            min_index = int(source.min())
            max_index = int(source.max())
        else:
            min_index = 0
            max_index = -1
        int32 = np.iinfo(np.int32)
        if min_index < int32.min or max_index > int32.max:
            raise BackendExecutionError(
                "prepare_gpu_solver_inputs",
                f"{path}: integer values exceed int32 GPU index range",
                details={
                    "path": path,
                    "source_dtype": str(source.dtype),
                    "min": min_index,
                    "max": max_index,
                    "gpu_dtype": "i32",
                },
            )
        return qd.i32, "i32", source.astype(np.int32, copy=False)

    if np.issubdtype(source.dtype, np.floating):
        finite = np.isfinite(source)
        if not bool(np.all(finite)):
            raise BackendExecutionError(
                "prepare_gpu_solver_inputs",
                f"{path}: non-finite values cannot be copied into GPU solver inputs",
                details={
                    "path": path,
                    "source_dtype": str(source.dtype),
                    "shape": array_shape(source),
                    "nonfinite_count": int(source.size - np.count_nonzero(finite)),
                },
            )
        if source.dtype == np.dtype("float32"):
            return qd.f32, "f32", source
        return qd.f64, "f64", source.astype(np.float64, copy=False)

    raise BackendExecutionError(
        "prepare_gpu_solver_inputs",
        f"{path}: unsupported GPU input dtype {source.dtype}",
        details={"path": path, "source_dtype": str(source.dtype), "shape": array_shape(source)},
    )


def transfer_solver_state_to_gpu(state: Mapping[str, Any]) -> dict[str, Any]:
    import quadrants as qd

    arrays: dict[str, dict[str, Any]] = {}
    device_arrays: list[Any] = []
    transferred_count = 0
    source_bytes = 0
    gpu_bytes = 0

    for path, source in solver_state_arrays(state):
        source_array = np.asarray(source)
        source_bytes += int(source_array.nbytes)
        entry: dict[str, Any] = {
            "path": path,
            "source_shape": array_shape(source_array),
            "source_dtype": str(source_array.dtype),
            "source_contiguous": bool(source_array.flags.c_contiguous),
            "source_nbytes": int(source_array.nbytes),
            "stats": array_stats(source_array),
        }
        if source_array.size == 0:
            entry.update(
                {
                    "status": "empty_array",
                    "transferred": False,
                    "gpu_shape": array_shape(source_array),
                    "reason": "zero-sized OpenFOAM patch array has no device allocation",
                }
            )
            arrays[path] = entry
            continue

        qd_dtype, gpu_dtype, gpu_source = gpu_transfer_array_dtype(qd, source_array, path)
        gpu_source = np.ascontiguousarray(gpu_source)
        gpu_array = qd.ndarray(qd_dtype, shape=gpu_source.shape)
        gpu_array.from_numpy(gpu_source)
        device_arrays.append(gpu_array)
        transferred_count += 1
        gpu_bytes += int(gpu_source.nbytes)
        entry.update(
            {
                "status": "transferred",
                "transferred": True,
                "gpu_shape": array_shape(gpu_source),
                "gpu_dtype": gpu_dtype,
                "gpu_nbytes": int(gpu_source.nbytes),
                "contiguous_for_transfer": True,
            }
        )
        arrays[path] = entry

    qd.sync()
    return {
        "framework": "quadrants",
        "device": "cuda",
        "array_count": len(arrays),
        "transferred_count": transferred_count,
        "empty_array_count": len(arrays) - transferred_count,
        "source_nbytes": source_bytes,
        "gpu_nbytes": gpu_bytes,
        "arrays": arrays,
    }


def gpu_required_field_inputs(state: Mapping[str, Any], transfer: Mapping[str, Any]) -> dict[str, Any]:
    arrays = transfer.get("arrays") if isinstance(transfer.get("arrays"), Mapping) else {}
    fields = state.get("fields") if isinstance(state.get("fields"), Mapping) else {}
    out: dict[str, Any] = {}
    for name in REQUIRED_FIELDS:
        field = fields.get(name) if isinstance(fields.get(name), Mapping) else {}
        boundary = field.get("boundary") if isinstance(field.get("boundary"), Mapping) else {}
        internal_path = f"fields.{name}.internal"
        internal = arrays.get(internal_path) if isinstance(arrays.get(internal_path), Mapping) else {}
        out[name] = {
            "status": "gpu_represented" if internal.get("transferred") is True else "missing_gpu_representation",
            "rank": field.get("rank"),
            "entity_kind": field.get("entity_kind"),
            "entity_count": field.get("entity_count"),
            "internal_path": internal_path,
            "internal": internal,
            "boundary_patch_count": len(boundary),
            "boundary": {
                patch_name: {
                    "type": patch.get("type"),
                    "values": arrays.get(f"fields.{name}.boundary.{patch_name}.values"),
                    "patch_internal": arrays.get(f"fields.{name}.boundary.{patch_name}.patch_internal"),
                    "value_internal_coeffs": arrays.get(f"fields.{name}.boundary.{patch_name}.value_internal_coeffs"),
                    "value_boundary_coeffs": arrays.get(f"fields.{name}.boundary.{patch_name}.value_boundary_coeffs"),
                    "gradient_internal_coeffs": arrays.get(f"fields.{name}.boundary.{patch_name}.gradient_internal_coeffs"),
                    "gradient_boundary_coeffs": arrays.get(f"fields.{name}.boundary.{patch_name}.gradient_boundary_coeffs"),
                }
                for patch_name, patch in boundary.items()
            },
        }
    return out


def gpu_mesh_input_summary(state: Mapping[str, Any], transfer: Mapping[str, Any]) -> dict[str, Any]:
    arrays = transfer.get("arrays") if isinstance(transfer.get("arrays"), Mapping) else {}
    mesh = state.get("mesh") if isinstance(state.get("mesh"), Mapping) else {}
    sizes = mesh.get("sizes") if isinstance(mesh.get("sizes"), Mapping) else {}
    connectivity = mesh.get("connectivity") if isinstance(mesh.get("connectivity"), Mapping) else {}
    owner = np.asarray(connectivity.get("owner", []))
    neighbour = np.asarray(connectivity.get("neighbour", []))
    patches = mesh.get("patches") if isinstance(mesh.get("patches"), list) else []
    owner_range = [int(owner.min()), int(owner.max())] if owner.size else [None, None]
    neighbour_range = [int(neighbour.min()), int(neighbour.max())] if neighbour.size else [None, None]
    return {
        "sizes": json_ready(sizes),
        "topology_checks": {
            "owner_shape": array_shape(owner),
            "neighbour_shape": array_shape(neighbour),
            "owner_cell_range": owner_range,
            "neighbour_cell_range": neighbour_range,
            "owner_neighbour_gpu_dtype": "i32",
            "owner_neighbour_in_cell_range": bool(
                owner.size
                and neighbour.size
                and min(owner_range[0], neighbour_range[0]) >= 0
                and max(owner_range[1], neighbour_range[1]) < int(sizes.get("n_cells", 0))
            ),
        },
        "connectivity": {
            "faces_offsets": arrays.get("mesh.connectivity.faces.offsets"),
            "faces_values": arrays.get("mesh.connectivity.faces.values"),
            "cells_offsets": arrays.get("mesh.connectivity.cells.offsets"),
            "cells_values": arrays.get("mesh.connectivity.cells.values"),
            "owner": arrays.get("mesh.connectivity.owner"),
            "neighbour": arrays.get("mesh.connectivity.neighbour"),
            "ldu_lower_addr": arrays.get("mesh.connectivity.ldu.lower_addr"),
            "ldu_upper_addr": arrays.get("mesh.connectivity.ldu.upper_addr"),
        },
        "geometry": {
            "points": arrays.get("mesh.geometry.points"),
            "V": arrays.get("mesh.geometry.V"),
            "C": arrays.get("mesh.geometry.C"),
            "Cf": arrays.get("mesh.geometry.Cf"),
            "Sf": arrays.get("mesh.geometry.Sf"),
            "magSf": arrays.get("mesh.geometry.magSf"),
        },
        "patches": [
            {
                "name": patch.get("name"),
                "type": patch.get("type"),
                "index": patch.get("index"),
                "start": patch.get("start"),
                "size": patch.get("size"),
                "coupled": patch.get("coupled"),
                "constraint": patch.get("constraint"),
                "face_cells": arrays.get(f"mesh.patches[{index}].face_cells"),
                "face_indices": arrays.get(f"mesh.patches[{index}].face_indices"),
                "Cf": arrays.get(f"mesh.patches[{index}].Cf"),
                "Sf": arrays.get(f"mesh.patches[{index}].Sf"),
                "magSf": arrays.get(f"mesh.patches[{index}].magSf"),
            }
            for index, patch in enumerate(patches)
        ],
    }


def gpu_matrix_input_blockers() -> dict[str, Any]:
    return {
        "UEqn": {
            "status": "provided_by_gpu_stage_kernel",
            "stage_group": "momentum_assembly",
            "source_export": "gpu_solver_stages.matrices.UEqn",
            "kernel_entrypoint": "gpu_rans_momentum_assembly",
            "required_arrays": ["diag", "source", "psi", "H"],
        },
        "pEqn": {
            "status": "provided_by_gpu_stage_kernel",
            "stage_group": "pressure_assembly",
            "source_export": "gpu_solver_stages.matrices.pEqn",
            "kernel_entrypoint": "gpu_rans_pressure_assembly",
            "required_arrays": ["diag", "source", "psi", "flux"],
        },
    }


def prepare_gpu_solver_inputs(foam: Any, stepper: Any, case: Path, prepared: Mapping[str, Any], backend: Mapping[str, Any]) -> dict[str, Any]:
    fields = read_fields(stepper, "gpu_inputs")
    solver_state = export_solver_state_checked(foam, stepper, "gpu_inputs", fields)
    try:
        foam.validate_solver_state(solver_state, required_fields=REQUIRED_FIELDS, turbulence_fields=TURBULENCE_FIELDS)
        transfer = transfer_solver_state_to_gpu(solver_state)
    except BackendExecutionError:
        raise
    except Exception as exc:
        raise BackendExecutionError(
            "prepare_gpu_solver_inputs",
            "failed to validate or transfer GPU solver input arrays",
            details=validation_details(exc),
        ) from exc

    prepared_identity = prepared.get("prepared_case_identity") if isinstance(prepared.get("prepared_case_identity"), Mapping) else {}
    return {
        "schema_version": GPU_INPUT_SCHEMA_VERSION,
        "status": "validated_and_transferred",
        "case": case,
        "backend": {
            "selected": backend.get("selected"),
            "provider": backend.get("provider"),
            "device": backend.get("device"),
            "used_cpu_fallback": backend.get("used_cpu_fallback"),
        },
        "source": {
            "state_export": "foam.export_solver_state(stepper.mesh(), stepper.fields(), required_fields=REQUIRED_FIELDS)",
            "derived_from_case_role": "run_one",
            "same_prepared_case_as_oracle": bool(
                prepared_identity.get("oracle", {}).get("matches_prepared")
                and prepared_identity.get("run_one", {}).get("matches_prepared")
            ),
            "prepared_case_identity": prepared_identity,
        },
        "validation": {
            "passed": True,
            "validator": "foam_stepper.validate_solver_state",
            "required_fields": list(REQUIRED_FIELDS),
            "turbulence_fields": list(TURBULENCE_FIELDS),
            "shape_dtype_topology_checked_before_kernel_launch": True,
        },
        "state_summary": foam.describe_solver_state(solver_state),
        "mesh": gpu_mesh_input_summary(solver_state, transfer),
        "required_fields": gpu_required_field_inputs(solver_state, transfer),
        "turbulence": {
            "fields": list(TURBULENCE_FIELDS),
            "source": "solver_state.turbulence plus required field GPU arrays",
            "field_inputs": {name: f"required_fields.{name}" for name in TURBULENCE_FIELDS},
        },
        "matrix_data": gpu_matrix_input_blockers(),
        "transfer": transfer,
    }

def gpu_f64_array(source: Any, path: str) -> tuple[Any, np.ndarray, dict[str, Any]]:
    source_array = np.asarray(source)
    qd_dtype, gpu_dtype, gpu_source = gpu_transfer_array_dtype(qd, source_array, path)
    if gpu_dtype != "f64":
        gpu_source = gpu_source.astype(np.float64, copy=False)
        gpu_dtype = "f64"
        qd_dtype = qd.f64
    gpu_source = np.ascontiguousarray(gpu_source)
    gpu_array = qd.ndarray(qd_dtype, shape=gpu_source.shape)
    gpu_array.from_numpy(gpu_source)
    return gpu_array, gpu_source, {
        "path": path,
        "source_shape": array_shape(source_array),
        "source_dtype": str(source_array.dtype),
        "gpu_shape": array_shape(gpu_source),
        "gpu_dtype": gpu_dtype,
        "gpu_nbytes": int(gpu_source.nbytes),
        "transferred": True,
    }

def gpu_i32_array(source: Any, path: str) -> tuple[Any, np.ndarray, dict[str, Any]]:
    source_array = np.asarray(source)
    qd_dtype, gpu_dtype, gpu_source = gpu_transfer_array_dtype(qd, source_array, path)
    if gpu_dtype != "i32":
        raise BackendExecutionError(
            "prepare_gpu_solver_inputs",
            f"{path}: expected integer topology for GPU index array",
            details={"path": path, "source_dtype": str(source_array.dtype), "gpu_dtype": gpu_dtype},
        )
    gpu_source = np.ascontiguousarray(gpu_source)
    gpu_array = qd.ndarray(qd_dtype, shape=gpu_source.shape)
    gpu_array.from_numpy(gpu_source)
    return gpu_array, gpu_source, {
        "path": path,
        "source_shape": array_shape(source_array),
        "source_dtype": str(source_array.dtype),
        "gpu_shape": array_shape(gpu_source),
        "gpu_dtype": gpu_dtype,
        "gpu_nbytes": int(gpu_source.nbytes),
        "transferred": True,
    }


def build_openfoam_offdiag_cell_faces(
    n_cells: int,
    owner: np.ndarray,
    neighbour: np.ndarray,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    owner_i32 = np.asarray(owner, dtype=np.int32).reshape(-1)
    neighbour_i32 = np.asarray(neighbour, dtype=np.int32).reshape(-1)
    if owner_i32.shape != neighbour_i32.shape:
        raise ValueError(f"owner/neighbour shape mismatch: {owner_i32.shape} != {neighbour_i32.shape}")
    if owner_i32.size == 0:
        return np.zeros(n_cells + 1, dtype=np.int32), np.zeros(0, dtype=np.int32), np.zeros(0, dtype=np.int32)
    if int(owner_i32.min()) < 0 or int(neighbour_i32.min()) < 0 or int(max(owner_i32.max(), neighbour_i32.max())) >= n_cells:
        raise ValueError("owner/neighbour addresses exceed cell range")

    counts = np.zeros(n_cells, dtype=np.int32)
    np.add.at(counts, owner_i32, 1)
    np.add.at(counts, neighbour_i32, 1)
    offsets = np.empty(n_cells + 1, dtype=np.int32)
    offsets[0] = 0
    np.cumsum(counts, out=offsets[1:])
    faces = np.empty(int(offsets[-1]), dtype=np.int32)
    sides = np.empty(int(offsets[-1]), dtype=np.int32)
    cursor = offsets[:-1].copy()

    for face, owner_cell in enumerate(owner_i32):
        owner_slot = int(cursor[int(owner_cell)])
        faces[owner_slot] = int(face)
        sides[owner_slot] = 0
        cursor[int(owner_cell)] += 1

        neighbour_cell = int(neighbour_i32[face])
        neighbour_slot = int(cursor[neighbour_cell])
        faces[neighbour_slot] = int(face)
        sides[neighbour_slot] = 1
        cursor[neighbour_cell] += 1

    return (
        np.ascontiguousarray(offsets, dtype=np.int32),
        np.ascontiguousarray(faces, dtype=np.int32),
        np.ascontiguousarray(sides, dtype=np.int32),
    )


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 openfoam_linear_face_weights(
    owner: np.ndarray,
    neighbour: np.ndarray,
    cell_centres: np.ndarray,
    face_centres: np.ndarray,
    face_area_vectors: np.ndarray,
) -> np.ndarray:
    """Return OpenFOAM linear surface-interpolation owner weights for internal faces."""

    owner_cells = np.asarray(owner, dtype=np.int32).reshape(-1)
    neighbour_cells = np.asarray(neighbour, dtype=np.int32).reshape(-1)
    owner_centres = np.asarray(cell_centres, dtype=np.float64)[owner_cells]
    neighbour_centres = np.asarray(cell_centres, dtype=np.float64)[neighbour_cells]
    face_c = np.asarray(face_centres, dtype=np.float64).reshape(-1, 3)
    area = np.asarray(face_area_vectors, dtype=np.float64).reshape(-1, 3)
    owner_delta = face_c - owner_centres
    neighbour_delta = neighbour_centres - face_c
    owner_distance = np.abs(np.sum(area * owner_delta, axis=1))
    neighbour_distance = np.abs(np.sum(area * neighbour_delta, axis=1))
    denominator = owner_distance + neighbour_distance
    return np.divide(
        neighbour_distance,
        denominator,
        out=np.full(owner_cells.shape, 0.5, dtype=np.float64),
        where=denominator > 1.0e-300,
    )


def openfoam_ldu_normalization_factor(
    n_cells: int,
    owner: np.ndarray,
    neighbour: np.ndarray,
    diag: np.ndarray,
    upper: np.ndarray,
    lower: np.ndarray,
    source: np.ndarray,
    psi: np.ndarray,
) -> float:
    """Return OpenFOAM lduMatrix::solver::normFactor for one scalar component."""

    owner_cells = np.asarray(owner, dtype=np.int32).reshape(-1)
    neighbour_cells = np.asarray(neighbour, dtype=np.int32).reshape(-1)
    diagonal = np.asarray(diag, dtype=np.float64).reshape(n_cells)
    upper_coeff = np.asarray(upper, dtype=np.float64).reshape(-1)
    lower_coeff = np.asarray(lower, dtype=np.float64).reshape(-1)
    field = np.asarray(psi, dtype=np.float64).reshape(n_cells)
    rhs = np.asarray(source, dtype=np.float64).reshape(n_cells)

    a_psi = diagonal * field
    sum_a = diagonal.copy()
    np.add.at(a_psi, owner_cells, upper_coeff * field[neighbour_cells])
    np.add.at(a_psi, neighbour_cells, lower_coeff * field[owner_cells])
    np.add.at(sum_a, owner_cells, upper_coeff)
    np.add.at(sum_a, neighbour_cells, lower_coeff)

    reference = sum_a * float(np.mean(field))
    return float(np.sum(np.abs(a_psi - reference) + np.abs(rhs - reference)) + 1.0e-300)


def gpu_empty_f64(shape: tuple[int, ...], path: str) -> tuple[Any, dict[str, Any]]:
    gpu_array = qd.ndarray(qd.f64, shape=shape)
    return gpu_array, {"path": path, "gpu_shape": list(shape), "gpu_dtype": "f64", "allocated": True}


def read_case_laminar_nu(case: Path) -> float:
    for relative in ("constant/physicalProperties", "constant/transportProperties", "constant/transportProperties.v2112"):
        path = case / relative
        if not path.exists():
            continue
        for line in path.read_text().splitlines():
            stripped = line.strip()
            if not stripped.startswith("nu"):
                continue
            parts = stripped.replace(";", " ").split()
            if len(parts) >= 2 and parts[0] == "nu":
                try:
                    return float(parts[1])
                except ValueError:
                    break
    return DEFAULT_LAMINAR_NU


def openfoam_momentum_diffusion_face_coefficients(
    owner: np.ndarray,
    neighbour: np.ndarray,
    nut_internal: np.ndarray,
    laminar_nu: float,
    cell_centres: np.ndarray,
    face_area_vectors: np.ndarray,
    face_area_magnitudes: np.ndarray,
    face_weights: np.ndarray,
) -> np.ndarray:
    owner_cells = np.asarray(owner, dtype=np.int32).reshape(-1)
    neighbour_cells = np.asarray(neighbour, dtype=np.int32).reshape(-1)
    centres = np.asarray(cell_centres, dtype=np.float64)
    sf = np.asarray(face_area_vectors, dtype=np.float64).reshape(-1, 3)
    mag_sf = np.asarray(face_area_magnitudes, dtype=np.float64).reshape(-1)
    weights = np.asarray(face_weights, dtype=np.float64).reshape(-1)
    nut = np.asarray(nut_internal, dtype=np.float64).reshape(-1)
    delta = centres[neighbour_cells] - centres[owner_cells]
    projected_delta = np.abs(np.sum(delta * sf, axis=1))
    effective_nu = laminar_nu + weights * nut[owner_cells] + (1.0 - weights) * nut[neighbour_cells]
    return np.ascontiguousarray(
        effective_nu * mag_sf * mag_sf / np.maximum(projected_delta, 1.0e-300),
        dtype=np.float64,
    )




def openfoam_momentum_internal_diag_from_coefficients(
    n_cells: int,
    owner: np.ndarray,
    neighbour: np.ndarray,
    diffusion_coeff: np.ndarray,
    phi: np.ndarray,
) -> np.ndarray:
    owner_cells = np.asarray(owner, dtype=np.int32).reshape(-1)
    neighbour_cells = np.asarray(neighbour, dtype=np.int32).reshape(-1)
    coeff = np.asarray(diffusion_coeff, dtype=np.float64).reshape(-1)
    flux = np.asarray(phi, dtype=np.float64).reshape(-1)
    owner_contrib = coeff + np.where(flux < 0.0, -flux, 0.0)
    neighbour_contrib = coeff + np.where(flux >= 0.0, flux, 0.0)
    diag = np.zeros(n_cells, dtype=np.float64)
    np.add.at(diag, owner_cells, owner_contrib)
    np.add.at(diag, neighbour_cells, neighbour_contrib)
    return np.ascontiguousarray(diag, dtype=np.float64)


def openfoam_momentum_pressure_source(
    n_cells: int,
    owner: np.ndarray,
    neighbour: np.ndarray,
    pressure: np.ndarray,
    face_weights: np.ndarray,
    sf: np.ndarray,
    boundary_face_cells: np.ndarray,
    boundary_values: np.ndarray,
    boundary_sf: np.ndarray,
) -> np.ndarray:
    owner_cells = np.asarray(owner, dtype=np.int32).reshape(-1)
    neighbour_cells = np.asarray(neighbour, dtype=np.int32).reshape(-1)
    p = np.asarray(pressure, dtype=np.float64).reshape(-1)
    weights = np.asarray(face_weights, dtype=np.float64).reshape(-1)
    area = np.asarray(sf, dtype=np.float64).reshape(-1, 3)
    face_pressure = weights * p[owner_cells] + (1.0 - weights) * p[neighbour_cells]
    source = np.zeros((n_cells, 3), dtype=np.float64)
    for direction in range(3):
        flux = face_pressure * area[:, direction]
        np.add.at(source[:, direction], owner_cells, -flux)
        np.add.at(source[:, direction], neighbour_cells, flux)
    if np.asarray(boundary_values).size:
        patch_cells = np.asarray(boundary_face_cells, dtype=np.int32).reshape(-1)
        patch_pressure = np.asarray(boundary_values, dtype=np.float64).reshape(-1)
        patch_area = np.asarray(boundary_sf, dtype=np.float64).reshape(-1, 3)
        for direction in range(3):
            np.add.at(source[:, direction], patch_cells, -patch_pressure * patch_area[:, direction])
    return np.ascontiguousarray(source, dtype=np.float64)


def gpu_array_output(name: str, array: np.ndarray, *, kernel: str) -> dict[str, Any]:
    return {
        "name": name,
        "kernel_entrypoint": kernel,
        "shape": array_shape(array),
        "dtype": str(array.dtype),
        "stats": array_stats(array),
    }


def gpu_matrix_output(
    name: str,
    field_name: str,
    diag: np.ndarray,
    source: np.ndarray,
    psi: np.ndarray,
    *,
    kernel: str,
    upper: np.ndarray | None = None,
    upper_kernel: str | None = None,
    lower: np.ndarray | None = None,
    lower_kernel: str | None = None,
) -> dict[str, Any]:
    out = {
        "name": name,
        "field_name": field_name,
        "value_rank": "vector" if source.ndim == 2 else "scalar",
        "kernel_entrypoint": kernel,
        "diag": gpu_array_output(f"{name}.diag", diag, kernel=kernel),
        "source": gpu_array_output(f"{name}.source", source, kernel=kernel),
        "psi": gpu_array_output(f"{name}.psi", psi, kernel=kernel),
        "has_upper": upper is not None,
        "has_lower": lower is not None,
        "gpu_backed": True,
    }
    if upper is not None:
        out["upper"] = gpu_array_output(f"{name}.upper", upper, kernel=upper_kernel or kernel)
    if lower is not None:
        out["lower"] = gpu_array_output(f"{name}.lower", lower, kernel=lower_kernel or kernel)
    return out


def gpu_stage_result(name: str, outputs: Mapping[str, Any], *, kernels: Iterable[str], changed_fields: Iterable[str] = ()) -> GpuTransformResult:
    return GpuTransformResult(
        name=name,
        phase="gpu_solver_stage",
        source=GpuSourceLocation(file=__file__, function=name),
        inputs={},
        outputs=dict(outputs),
        changed_fields=list(changed_fields),
        metadata={
            "backend": "gpu",
            "framework": "quadrants",
            "kernel_entrypoints": list(kernels),
            "used_cpu_fallback": False,
        },
    )


def openfoam_momentum_matrix_coefficients(stepper: Any) -> tuple[np.ndarray, list[np.ndarray], list[np.ndarray]] | None:
    try:
        if hasattr(stepper, "momentum_transport_predictor"):
            stepper.momentum_transport_predictor()
        if hasattr(stepper, "assemble_momentum_terms"):
            stepper.assemble_momentum_terms()
        matrix_result = stepper.assemble_momentum_matrix()
        matrix = matrix_result.outputs.get("UEqn") if isinstance(matrix_result.outputs, Mapping) else None
        source = getattr(matrix, "source", None)
        internal = getattr(matrix, "internal_coeffs", None)
        boundary = getattr(matrix, "boundary_coeffs", None)
        if source is None or not isinstance(internal, list) or not isinstance(boundary, list) or len(internal) != len(boundary):
            return None
        return (
            np.ascontiguousarray(np.asarray(source, dtype=np.float64)),
            [np.ascontiguousarray(np.asarray(coeffs, dtype=np.float64)) for coeffs in internal],
            [np.ascontiguousarray(np.asarray(coeffs, dtype=np.float64)) for coeffs in boundary],
        )
    except Exception:
        return None


def clone_gpu_patch_field(patch: Any) -> GpuPatchField:
    return GpuPatchField(
        name=getattr(patch, "name", ""),
        type=getattr(patch, "type", ""),
        values=np.asarray(getattr(patch, "values")),
        fixes_value=bool(getattr(patch, "fixes_value", False)),
        assignable=bool(getattr(patch, "assignable", False)),
        coupled=bool(getattr(patch, "coupled", False)),
        updated=bool(getattr(patch, "updated", False)),
        patch_internal=None if getattr(patch, "patch_internal", None) is None else np.asarray(getattr(patch, "patch_internal")),
        value_internal_coeffs=None if getattr(patch, "value_internal_coeffs", None) is None else np.asarray(getattr(patch, "value_internal_coeffs")),
        value_boundary_coeffs=None if getattr(patch, "value_boundary_coeffs", None) is None else np.asarray(getattr(patch, "value_boundary_coeffs")),
        gradient_internal_coeffs=None if getattr(patch, "gradient_internal_coeffs", None) is None else np.asarray(getattr(patch, "gradient_internal_coeffs")),
        gradient_boundary_coeffs=None if getattr(patch, "gradient_boundary_coeffs", None) is None else np.asarray(getattr(patch, "gradient_boundary_coeffs")),
    )


def gpu_field_from_source(source: Any, internal: np.ndarray) -> GpuField:
    boundary = getattr(source, "boundary", {})
    return GpuField(
        name=getattr(source, "name", ""),
        kind=getattr(source, "kind", ""),
        dimensions=getattr(source, "dimensions", ""),
        entity_kind=getattr(source, "entity_kind", ""),
        entity_count=int(getattr(source, "entity_count", np.asarray(internal).shape[0] if np.asarray(internal).ndim else 0)),
        internal=np.ascontiguousarray(internal),
        boundary={name: clone_gpu_patch_field(patch) for name, patch in boundary.items()},
    )


def gpu_solver_stage_report(gpu_run: Mapping[str, Any]) -> dict[str, Any]:
    return {
        key: value
        for key, value in gpu_run.items()
        if key not in {"field_objects", "stage_objects"}
    }


def quadrants_kernel_evidence(expected: Iterable[str]) -> dict[str, Any]:
    try:
        from quadrants.profiler.kernel_profiler import get_default_kernel_profiler

        profiler = get_default_kernel_profiler()
        profiler._update_records()
        records = list(profiler._traced_records)
        generated_kernel_names = sorted({str(record.name) for record in records})
        return {
            "available": True,
            "expected_kernel_entrypoints": list(expected),
            "generated_kernel_names": generated_kernel_names,
            "profile_record_count": len(records),
            "device_time_ms_total": float(sum(record.kernel_time for record in records)),
        }
    except Exception as exc:
        return {
            "available": False,
            "expected_kernel_entrypoints": list(expected),
            "cause": {"type": type(exc).__name__, "message": str(exc)},
        }


def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: Path, *, 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()
    u_gpu, u_np, u_transfer = gpu_f64_array(fields["U"].internal, "fields.U.internal")
    p_gpu, p_np, p_transfer = gpu_f64_array(fields["p"].internal, "fields.p.internal")
    phi_gpu, phi_np, phi_transfer = gpu_f64_array(fields["phi"].internal, "fields.phi.internal")
    nut_gpu, nut_np, nut_transfer = gpu_f64_array(fields["nut"].internal, "fields.nut.internal")
    k_gpu, k_np, k_transfer = gpu_f64_array(fields["k"].internal, "fields.k.internal")
    omega_gpu, omega_np, omega_transfer = gpu_f64_array(fields["omega"].internal, "fields.omega.internal")
    n_cells = int(u_np.shape[0])
    n_internal_faces = int(phi_np.shape[0])
    mesh = stepper.mesh()
    laminar_nu = read_case_laminar_nu(case)
    owner_gpu, owner_np, owner_transfer = gpu_i32_array(np.asarray(mesh.owner)[:n_internal_faces], "mesh.connectivity.owner.internal")
    neighbour_gpu, neighbour_np, neighbour_transfer = gpu_i32_array(np.asarray(mesh.neighbour)[:n_internal_faces], "mesh.connectivity.neighbour.internal")
    losort_np = build_losort_addr(n_cells, neighbour_np)
    losort_gpu, losort_np, losort_transfer = gpu_i32_array(losort_np, "mesh.connectivity.losort.internal")
    offdiag_cell_offsets_np, offdiag_cell_faces_np, offdiag_cell_sides_np = build_openfoam_offdiag_cell_faces(n_cells, owner_np, neighbour_np)
    offdiag_cell_offsets_gpu, offdiag_cell_offsets_np, offdiag_cell_offsets_transfer = gpu_i32_array(offdiag_cell_offsets_np, "mesh.connectivity.offdiag_cell_offsets.internal")
    offdiag_cell_faces_gpu, offdiag_cell_faces_np, offdiag_cell_faces_transfer = gpu_i32_array(offdiag_cell_faces_np, "mesh.connectivity.offdiag_cell_faces.internal")
    offdiag_cell_sides_gpu, offdiag_cell_sides_np, offdiag_cell_sides_transfer = gpu_i32_array(offdiag_cell_sides_np, "mesh.connectivity.offdiag_cell_sides.internal")
    u_dilu_level_schedule = build_ldu_level_schedule(n_cells, owner_np, neighbour_np, name="solve_UEqn.DILU")
    p_dic_level_schedule = build_ldu_level_schedule(n_cells, owner_np, neighbour_np, name="solve_pEqn.DIC")
    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")
    face_weights_np = openfoam_linear_face_weights(owner_np, neighbour_np, cell_centres_np, face_centres_np, sf_np)
    face_weights_gpu, face_weights_np, face_weights_transfer = gpu_f64_array(face_weights_np, "mesh.geometry.surface_interpolation.owner_weights.internal")
    momentum_diffusion_coeff_np = openfoam_momentum_diffusion_face_coefficients(owner_np, neighbour_np, nut_np, laminar_nu, cell_centres_np, sf_np, mag_sf_np, face_weights_np)
    momentum_diffusion_coeff_gpu, momentum_diffusion_coeff_np, momentum_diffusion_coeff_transfer = gpu_f64_array(momentum_diffusion_coeff_np, "mesh.geometry.momentum_diffusion_coefficients.internal")
    momentum_internal_diag_np = openfoam_momentum_internal_diag_from_coefficients(n_cells, owner_np, neighbour_np, momentum_diffusion_coeff_np, phi_np)
    momentum_internal_diag_gpu, momentum_internal_diag_np, momentum_internal_diag_transfer = gpu_f64_array(momentum_internal_diag_np, "matrix.UEqn.internal_diag_from_offdiag")
    momentum_internal_h1_np = np.ascontiguousarray(momentum_internal_diag_np / cell_volumes_np, dtype=np.float64)
    momentum_internal_h1_gpu, momentum_internal_h1_np, momentum_internal_h1_transfer = gpu_f64_array(momentum_internal_h1_np, "matrix.UEqn.internal_H1")
    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)
        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
        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(phi_values.shape[0]), "face_cell_size": int(face_cells.shape[0])},
            )
        boundary_face_cells_parts.append(face_cells)
        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)
    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_pressure_source_np = openfoam_momentum_pressure_source(n_cells, owner_np, neighbour_np, p_np, face_weights_np, sf_np, momentum_pressure_boundary_face_cells_np, momentum_pressure_boundary_values_np, momentum_pressure_boundary_sf_np)
    momentum_pressure_source_gpu, momentum_pressure_source_np, momentum_pressure_source_transfer = gpu_f64_array(momentum_pressure_source_np, "matrix.UEqn.pressure_gradient_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] = []
    native_momentum_matrix_coefficients = openfoam_momentum_matrix_coefficients(stepper)
    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":
            native_internal_coeffs = None
            native_boundary_coeffs = None
            patch_index = int(getattr(patch, "index", -1))
            if native_momentum_matrix_coefficients is not None and 0 <= patch_index < len(native_momentum_matrix_coefficients[1]):
                native_internal_coeffs = native_momentum_matrix_coefficients[1][patch_index]
                native_boundary_coeffs = native_momentum_matrix_coefficients[2][patch_index]
            if (
                native_internal_coeffs is not None
                and native_boundary_coeffs is not None
                and native_internal_coeffs.shape == (face_cells.shape[0], 3)
                and native_boundary_coeffs.shape == (face_cells.shape[0], 3)
            ):
                momentum_u_convection_boundary_source_face_cells_parts.append(face_cells)
                momentum_u_convection_boundary_source_parts.append(native_boundary_coeffs)
                momentum_u_convection_boundary_diag_face_cells_parts.append(face_cells)
                momentum_u_convection_boundary_diag_coeff_parts.append(native_internal_coeffs)
                continue
            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.asarray(patch.magSf, dtype=np.float64).reshape(-1)
            face_centres = np.asarray(patch.Cf, dtype=np.float64).reshape(-1, 3)
            if face_area_magnitudes.shape[0] != face_cells.shape[0] or face_centres.shape != (face_cells.shape[0], 3):
                raise BackendExecutionError(
                    "prepare_momentum_convection_boundary_coefficients",
                    "freestreamVelocity geometry arrays have incompatible shapes",
                    details={
                        "patch": patch.name,
                        "magSf_shape": list(face_area_magnitudes.shape),
                        "Cf_shape": list(face_centres.shape),
                        "face_cell_size": int(face_cells.shape[0]),
                    },
                )
            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)
            face_delta = face_centres - cell_centres_np[face_cells]
            delta_coeffs = 1.0 / np.maximum(np.linalg.norm(face_delta, axis=1), 1.0e-300)
            nut_patch = fields["nut"].boundary.get(patch.name)
            nut_values = np.asarray(nut_patch.values, dtype=np.float64).reshape(-1) if nut_patch is not None else np.zeros(face_cells.shape[0], dtype=np.float64)
            if nut_values.shape[0] != face_cells.shape[0]:
                raise BackendExecutionError(
                    "prepare_momentum_convection_boundary_coefficients",
                    "freestreamVelocity nut boundary and face-cell arrays have incompatible shapes",
                    details={"patch": patch.name, "nut_shape": list(nut_values.shape), "face_cell_size": int(face_cells.shape[0])},
                )
            diffusion_internal_coeffs = (laminar_nu + nut_values) * face_area_magnitudes * delta_coeffs * value_fraction
            freestream_internal_coeffs = phi_values * (1.0 - value_fraction) + diffusion_internal_coeffs
            freestream_boundary_source = (-phi_values * value_fraction + diffusion_internal_coeffs)[:, None] * u_values
            momentum_u_convection_boundary_source_face_cells_parts.append(face_cells)
            momentum_u_convection_boundary_source_parts.append(freestream_boundary_source)
            momentum_u_convection_boundary_diag_face_cells_parts.append(face_cells)
            momentum_u_convection_boundary_diag_coeff_parts.append(np.repeat(freestream_internal_coeffs[:, None], 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_u_values_parts: list[np.ndarray] = []
    pressure_mixed_face_centres_parts: list[np.ndarray] = []
    pressure_mixed_area_vectors_parts: list[np.ndarray] = []
    pressure_mixed_area_magnitudes_parts: list[np.ndarray] = []
    for patch in mesh.boundary:
        p_patch = fields["p"].boundary.get(patch.name)
        u_patch = fields["U"].boundary.get(patch.name)
        if p_patch is None or u_patch is None or getattr(p_patch, "type", "") != "freestreamPressure":
            continue
        u_values = np.asarray(u_patch.values, dtype=np.float64).reshape(-1, 3)
        face_cells = np.asarray(patch.face_cells, dtype=np.int32).reshape(-1)
        face_centres = np.asarray(patch.Cf, dtype=np.float64).reshape(-1, 3)
        face_area_vectors = np.asarray(patch.Sf, dtype=np.float64).reshape(-1, 3)
        face_area_magnitudes = np.asarray(patch.magSf, dtype=np.float64).reshape(-1)
        if not (
            face_cells.shape[0]
            == u_values.shape[0]
            == face_centres.shape[0]
            == face_area_vectors.shape[0]
            == face_area_magnitudes.shape[0]
        ):
            raise BackendExecutionError(
                "prepare_pressure_mixed_solver_diag",
                "freestreamPressure patch arrays have incompatible shapes",
                details={
                    "patch": patch.name,
                    "face_cells_shape": list(face_cells.shape),
                    "u_values_shape": list(u_values.shape),
                    "face_centres_shape": list(face_centres.shape),
                    "face_area_vectors_shape": list(face_area_vectors.shape),
                    "face_area_magnitudes_shape": list(face_area_magnitudes.shape),
                },
            )
        pressure_mixed_face_cells_parts.append(face_cells)
        pressure_mixed_u_values_parts.append(u_values)
        pressure_mixed_face_centres_parts.append(face_centres)
        pressure_mixed_area_vectors_parts.append(face_area_vectors)
        pressure_mixed_area_magnitudes_parts.append(face_area_magnitudes)
    n_pressure_mixed_faces = int(sum(part.shape[0] for part in pressure_mixed_face_cells_parts))
    pressure_mixed_face_cells_np = np.concatenate(pressure_mixed_face_cells_parts) if pressure_mixed_face_cells_parts else np.zeros((1,), dtype=np.int32)
    pressure_mixed_u_values_np = np.concatenate(pressure_mixed_u_values_parts) if pressure_mixed_u_values_parts else np.zeros((1, 3), dtype=np.float64)
    pressure_mixed_face_centres_np = np.concatenate(pressure_mixed_face_centres_parts) if pressure_mixed_face_centres_parts else np.zeros((1, 3), dtype=np.float64)
    pressure_mixed_area_vectors_np = np.concatenate(pressure_mixed_area_vectors_parts) if pressure_mixed_area_vectors_parts else np.zeros((1, 3), dtype=np.float64)
    pressure_mixed_area_magnitudes_np = np.concatenate(pressure_mixed_area_magnitudes_parts) if pressure_mixed_area_magnitudes_parts else np.zeros((1,), dtype=np.float64)
    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_solver_diag")
    pressure_mixed_u_values_gpu, pressure_mixed_u_values_np, pressure_mixed_u_values_transfer = gpu_f64_array(pressure_mixed_u_values_np, "fields.U.boundary.pressure_mixed_solver_diag")
    pressure_mixed_face_centres_gpu, pressure_mixed_face_centres_np, pressure_mixed_face_centres_transfer = gpu_f64_array(pressure_mixed_face_centres_np, "mesh.boundary.Cf.pressure_mixed_solver_diag")
    pressure_mixed_area_vectors_gpu, pressure_mixed_area_vectors_np, pressure_mixed_area_vectors_transfer = gpu_f64_array(pressure_mixed_area_vectors_np, "mesh.boundary.Sf.pressure_mixed_solver_diag")
    pressure_mixed_area_magnitudes_gpu, pressure_mixed_area_magnitudes_np, pressure_mixed_area_magnitudes_transfer = gpu_f64_array(pressure_mixed_area_magnitudes_np, "mesh.boundary.magSf.pressure_mixed_solver_diag")
    momentum_wall_face_cells_parts: list[np.ndarray] = []
    momentum_wall_values_parts: list[np.ndarray] = []
    momentum_wall_nut_parts: list[np.ndarray] = []
    momentum_wall_face_centres_parts: list[np.ndarray] = []
    momentum_wall_area_vectors_parts: list[np.ndarray] = []
    momentum_wall_area_magnitudes_parts: list[np.ndarray] = []
    for patch in mesh.boundary:
        u_patch = fields["U"].boundary.get(patch.name)
        if u_patch is None or getattr(u_patch, "type", "") != "noSlip":
            continue
        values = np.asarray(u_patch.values, dtype=np.float64).reshape(-1, 3)
        if values.shape[0] == 0:
            continue
        face_cells = np.asarray(patch.face_cells, dtype=np.int32).reshape(-1)
        face_centres = np.asarray(patch.Cf, dtype=np.float64)
        face_area_vectors = np.asarray(patch.Sf, dtype=np.float64)
        face_area_magnitudes = np.asarray(patch.magSf, dtype=np.float64).reshape(-1)
        expected_wall_shape = (face_cells.shape[0], 3)
        if values.shape != expected_wall_shape or face_centres.shape != expected_wall_shape or face_area_vectors.shape != expected_wall_shape or face_area_magnitudes.shape[0] != face_cells.shape[0]:
            raise BackendExecutionError(
                "prepare_momentum_wall_diffusion",
                "wall velocity patch arrays have incompatible shapes",
                details={
                    "patch": patch.name,
                    "face_cells_shape": list(face_cells.shape),
                    "values_shape": list(values.shape),
                    "face_centres_shape": list(face_centres.shape),
                    "face_area_vectors_shape": list(face_area_vectors.shape),
                    "face_area_magnitudes_shape": list(face_area_magnitudes.shape),
                },
            )
        nut_patch = fields["nut"].boundary.get(patch.name)
        if nut_patch is None:
            nut_values = np.zeros((face_cells.shape[0],), dtype=np.float64)
        else:
            nut_values = np.asarray(nut_patch.values, dtype=np.float64).reshape(-1)
            if nut_values.shape[0] != face_cells.shape[0]:
                raise BackendExecutionError(
                    "prepare_momentum_wall_diffusion",
                    "wall nut patch values and face-cell arrays have different lengths",
                    details={"patch": patch.name, "nut_size": int(nut_values.shape[0]), "face_cell_size": int(face_cells.shape[0])},
                )
        momentum_wall_face_cells_parts.append(face_cells)
        momentum_wall_values_parts.append(values)
        momentum_wall_nut_parts.append(nut_values)
        momentum_wall_face_centres_parts.append(face_centres)
        momentum_wall_area_vectors_parts.append(face_area_vectors)
        momentum_wall_area_magnitudes_parts.append(face_area_magnitudes)
    momentum_wall_face_cells_np = np.concatenate(momentum_wall_face_cells_parts) if momentum_wall_face_cells_parts else np.empty((0,), dtype=np.int32)
    momentum_wall_values_np = np.concatenate(momentum_wall_values_parts) if momentum_wall_values_parts else np.empty((0, 3), dtype=np.float64)
    momentum_wall_nut_np = np.concatenate(momentum_wall_nut_parts) if momentum_wall_nut_parts else np.empty((0,), dtype=np.float64)
    momentum_wall_face_centres_np = np.concatenate(momentum_wall_face_centres_parts) if momentum_wall_face_centres_parts else np.empty((0, 3), dtype=np.float64)
    momentum_wall_area_vectors_np = np.concatenate(momentum_wall_area_vectors_parts) if momentum_wall_area_vectors_parts else np.empty((0, 3), dtype=np.float64)
    momentum_wall_area_magnitudes_np = np.concatenate(momentum_wall_area_magnitudes_parts) if momentum_wall_area_magnitudes_parts else np.empty((0,), dtype=np.float64)
    n_momentum_wall_faces = int(momentum_wall_face_cells_np.shape[0])
    momentum_wall_face_cells_gpu, momentum_wall_face_cells_np, momentum_wall_face_cells_transfer = gpu_i32_array(momentum_wall_face_cells_np, "mesh.boundary.face_cells.momentum_wall_diffusion")
    momentum_wall_values_gpu, momentum_wall_values_np, momentum_wall_values_transfer = gpu_f64_array(momentum_wall_values_np, "fields.U.boundary.momentum_wall_diffusion")
    momentum_wall_nut_gpu, momentum_wall_nut_np, momentum_wall_nut_transfer = gpu_f64_array(momentum_wall_nut_np, "fields.nut.boundary.momentum_wall_diffusion")
    momentum_wall_face_centres_gpu, momentum_wall_face_centres_np, momentum_wall_face_centres_transfer = gpu_f64_array(momentum_wall_face_centres_np, "mesh.boundary.Cf.momentum_wall_diffusion")
    momentum_wall_area_vectors_gpu, momentum_wall_area_vectors_np, momentum_wall_area_vectors_transfer = gpu_f64_array(momentum_wall_area_vectors_np, "mesh.boundary.Sf.momentum_wall_diffusion")
    momentum_wall_area_magnitudes_gpu, momentum_wall_area_magnitudes_np, momentum_wall_area_magnitudes_transfer = gpu_f64_array(momentum_wall_area_magnitudes_np, "mesh.boundary.magSf.momentum_wall_diffusion")
    omega_wall_distances_by_cell = np.full((n_cells,), np.inf, dtype=np.float64)
    omega_wall_seed_count = 0
    for patch in mesh.boundary:
        omega_patch = fields["omega"].boundary.get(patch.name)
        if omega_patch is None or "omegaWallFunction" not in omega_patch.type:
            continue
        face_cells = np.asarray(patch.face_cells, dtype=np.int32).reshape(-1)
        face_centres = np.asarray(patch.Cf, dtype=np.float64)
        face_area_vectors = np.asarray(patch.Sf, dtype=np.float64)
        expected_wall_shape = (face_cells.shape[0], 3)
        if face_centres.shape != expected_wall_shape or face_area_vectors.shape != expected_wall_shape:
            raise BackendExecutionError(
                "prepare_omega_wall_update",
                "omega wall patch geometry and face-cell arrays have incompatible shapes",
                details={
                    "patch": patch.name,
                    "face_cells_shape": list(face_cells.shape),
                    "face_centres_shape": list(face_centres.shape),
                    "face_area_vectors_shape": list(face_area_vectors.shape),
                },
            )
        wall_vectors = cell_centres_np[face_cells] - face_centres
        face_area_magnitudes = np.linalg.norm(face_area_vectors, axis=1)
        wall_distances = np.abs(np.sum(wall_vectors * face_area_vectors, axis=1)) / np.maximum(face_area_magnitudes, 1.0e-300)
        np.minimum.at(omega_wall_distances_by_cell, face_cells, wall_distances)
        omega_wall_seed_count += int(face_cells.shape[0])
    omega_wall_layer_count = 1
    for _ in range(omega_wall_layer_count):
        propagated_distances = omega_wall_distances_by_cell.copy()
        cell_centre_delta = np.linalg.norm(cell_centres_np[neighbour_np] - cell_centres_np[owner_np], axis=1)
        owner_has_wall = np.isfinite(omega_wall_distances_by_cell[owner_np])
        np.minimum.at(
            propagated_distances,
            neighbour_np[owner_has_wall],
            omega_wall_distances_by_cell[owner_np[owner_has_wall]] + cell_centre_delta[owner_has_wall],
        )
        neighbour_has_wall = np.isfinite(omega_wall_distances_by_cell[neighbour_np])
        np.minimum.at(
            propagated_distances,
            owner_np[neighbour_has_wall],
            omega_wall_distances_by_cell[neighbour_np[neighbour_has_wall]] + cell_centre_delta[neighbour_has_wall],
        )
        omega_wall_distances_by_cell = propagated_distances
    omega_wall_cells_np = np.flatnonzero(np.isfinite(omega_wall_distances_by_cell)).astype(np.int32, copy=False)
    omega_wall_distances_np = omega_wall_distances_by_cell[omega_wall_cells_np].astype(np.float64, copy=False)
    n_omega_wall_cells = int(omega_wall_cells_np.shape[0])
    omega_wall_cells_gpu, omega_wall_cells_np, omega_wall_cells_transfer = gpu_i32_array(omega_wall_cells_np, "mesh.boundary.cells.omega_wall")
    omega_wall_distances_gpu, omega_wall_distances_np, omega_wall_distances_transfer = gpu_f64_array(omega_wall_distances_np, "mesh.boundary.normal_wall_distance.omega")
    input_transfer_wall_seconds = time.perf_counter() - input_transfer_started
    u_diag_gpu, u_diag_alloc = gpu_empty_f64((n_cells,), "gpu_stages.UEqn.diag")
    u_source_gpu, u_source_alloc = gpu_empty_f64(tuple(u_np.shape), "gpu_stages.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")
    u_direction_gpu, u_direction_alloc = gpu_empty_f64(tuple(u_np.shape), "gpu_stages.solve_UEqn.direction")
    u_operator_direction_gpu, u_operator_direction_alloc = gpu_empty_f64(tuple(u_np.shape), "gpu_stages.solve_UEqn.operator_direction")
    u_intermediate_gpu, u_intermediate_alloc = gpu_empty_f64(tuple(u_np.shape), "gpu_stages.solve_UEqn.intermediate")
    u_operator_intermediate_gpu, u_operator_intermediate_alloc = gpu_empty_f64(tuple(u_np.shape), "gpu_stages.solve_UEqn.operator_intermediate")
    u_rr_gpu, u_rr_alloc = gpu_empty_f64((1,), "gpu_stages.solve_UEqn.residual_squared")
    u_rho_gpu, u_rho_alloc = gpu_empty_f64((1,), "gpu_stages.solve_UEqn.rho")
    u_denominator_gpu, u_denominator_alloc = gpu_empty_f64((1,), "gpu_stages.solve_UEqn.denominator")
    u_omega_numerator_gpu, u_omega_numerator_alloc = gpu_empty_f64((1,), "gpu_stages.solve_UEqn.omega_numerator")
    u_omega_denominator_gpu, u_omega_denominator_alloc = gpu_empty_f64((1,), "gpu_stages.solve_UEqn.omega_denominator")
    p_solved_gpu, p_solved_alloc = gpu_empty_f64(tuple(p_np.shape), "gpu_stages.solve_pEqn.p")
    p_work_gpu, p_work_alloc = gpu_empty_f64(tuple(p_np.shape), "gpu_stages.solve_pEqn.work")
    p_residual_gpu, p_residual_alloc = gpu_empty_f64(tuple(p_np.shape), "gpu_stages.solve_pEqn.residual")
    p_direction_gpu, p_direction_alloc = gpu_empty_f64(tuple(p_np.shape), "gpu_stages.solve_pEqn.direction")
    p_operator_gpu, p_operator_alloc = gpu_empty_f64(tuple(p_np.shape), "gpu_stages.solve_pEqn.operator")
    p_rr_gpu, p_rr_alloc = gpu_empty_f64((1,), "gpu_stages.solve_pEqn.residual_squared")
    p_denominator_gpu, p_denominator_alloc = gpu_empty_f64((1,), "gpu_stages.solve_pEqn.denominator")
    phi_solved_gpu, phi_solved_alloc = gpu_empty_f64(tuple(phi_np.shape), "gpu_stages.solve_pEqn.phi")
    u_final_gpu, u_final_alloc = gpu_empty_f64(tuple(u_np.shape), "gpu_stages.final_correction.U")
    p_final_gpu, p_final_alloc = gpu_empty_f64(tuple(p_np.shape), "gpu_stages.final_correction.p")
    nut_out_gpu, nut_alloc = gpu_empty_f64(tuple(nut_np.shape), "gpu_stages.turbulence.nut")
    k_out_gpu, k_alloc = gpu_empty_f64(tuple(k_np.shape), "gpu_stages.turbulence.k")
    omega_out_gpu, omega_alloc = gpu_empty_f64(tuple(omega_np.shape), "gpu_stages.turbulence.omega")

    try:
        qd.profiler.clear_kernel_profiler_info()
    except Exception:
        pass

    kernel_graph_started = time.perf_counter()
    gpu_rans_momentum_assembly(n_cells, u_gpu, 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_from_face_coeff(n_internal_faces, owner_gpu, neighbour_gpu, momentum_diffusion_coeff_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_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_bounded_convection_sp_internal(n_internal_faces, owner_gpu, neighbour_gpu, phi_gpu, u_diag_gpu)
    gpu_copy_scalar(n_cells, momentum_internal_diag_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, face_weights_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_internal_dev_tau_source(n_internal_faces, owner_gpu, neighbour_gpu, nut_gpu, laminar_nu, cell_centres_gpu, sf_gpu, mag_sf_gpu, face_weights_gpu, grad_u_gpu, u_source_gpu)
    if n_momentum_wall_faces:
        gpu_rans_momentum_wall_dev_tau_source(n_momentum_wall_faces, momentum_wall_face_cells_gpu, momentum_wall_values_gpu, momentum_wall_nut_gpu, laminar_nu, u_gpu, cell_centres_gpu, momentum_wall_face_centres_gpu, momentum_wall_area_vectors_gpu, momentum_wall_area_magnitudes_gpu, grad_u_gpu, u_source_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_copy_scalar(n_cells, momentum_internal_diag_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_add_vector_field(n_cells, momentum_pressure_source_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_source_for_solver = np.asarray(u_source_gpu.to_numpy(), dtype=np.float64)
    u_normalization_factors = np.asarray(
        [
            openfoam_ldu_normalization_factor(
                n_cells,
                owner_np,
                neighbour_np,
                u_boundary_diag_candidate_for_preconditioner,
                u_upper_for_preconditioner,
                u_lower_for_preconditioner,
                u_source_for_solver[:, component],
                u_np[:, component],
            )
            for component in range(3)
        ],
        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 level-scheduled GPU solver 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_components(
        n_cells,
        n_internal_faces,
        owner_gpu,
        neighbour_gpu,
        u_upper_gpu,
        u_lower_gpu,
        u_boundary_diag_candidate_gpu,
        u_source_gpu,
        u_gpu,
        u_solved_gpu,
        u_residual_gpu,
        u_shadow_gpu,
        u_direction_gpu,
        u_operator_direction_gpu,
        u_intermediate_gpu,
        u_operator_intermediate_gpu,
        u_rr_gpu,
        u_rho_gpu,
        u_denominator_gpu,
        u_omega_numerator_gpu,
        u_omega_denominator_gpu,
        preconditioner_diag=u_preconditioner_diag_gpu,
        preconditioner_reciprocal_diag=u_dilu_reciprocal_diag_gpu,
        dilu_schedule=u_dilu_level_schedule,
        iterations=DEFAULT_MOMENTUM_PBICGSTAB_ITERATIONS,
        residual_tolerance_squared=DEFAULT_MOMENTUM_PBICGSTAB_RESIDUAL_TOLERANCE_SQUARED,
        residual_tolerance=DEFAULT_MOMENTUM_PBICGSTAB_RESIDUAL_TOLERANCE,
        min_iterations=DEFAULT_MOMENTUM_PBICGSTAB_MIN_ITERATIONS,
        normalization_factors=u_normalization_factors,
    )
    gpu_rans_momentum_hbyA_by_cell(n_cells, offdiag_cell_offsets_gpu, offdiag_cell_faces_gpu, offdiag_cell_sides_gpu, owner_gpu, neighbour_gpu, u_upper_gpu, u_lower_gpu, u_solved_gpu, u_source_gpu, HbyA_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_copy_scalar(n_cells, momentum_internal_h1_gpu, H1_gpu)
    gpu_rans_consistent_rAtU(n_cells, rAU_gpu, H1_gpu, DEFAULT_SIMPLE_CONSISTENT_RATU_FACTOR, rAtU_gpu)
    gpu_rans_surface_flux_from_cells(n_internal_faces, owner_gpu, neighbour_gpu, HbyA_gpu, sf_gpu, face_weights_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, face_weights_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, face_weights_gpu, p_diag_gpu, p_upper_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_phiHbyA_gpu, p_source_gpu)
    gpu_rans_negate_scalar_field(n_cells, p_diag_gpu, p_solve_diag_gpu)
    if n_pressure_mixed_faces:
        gpu_rans_pressure_mixed_solver_diag(
            n_pressure_mixed_faces,
            pressure_mixed_face_cells_gpu,
            pressure_mixed_u_values_gpu,
            pressure_mixed_face_centres_gpu,
            pressure_mixed_area_vectors_gpu,
            pressure_mixed_area_magnitudes_gpu,
            cell_centres_gpu,
            rAtU_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_diag_for_normalization = -np.asarray(p_diag_gpu.to_numpy(), dtype=np.float64)
    p_solve_diag_for_preconditioner = np.asarray(p_solve_diag_gpu.to_numpy(), dtype=np.float64)
    p_solve_upper_for_preconditioner = np.asarray(p_solve_upper_gpu.to_numpy(), dtype=np.float64)
    p_solve_source_for_solver = np.asarray(p_solve_source_gpu.to_numpy(), dtype=np.float64)
    p_preconditioner_diag_np = openfoam_dic_preconditioner_diag(
        p_solve_diag_for_preconditioner,
        owner_np,
        neighbour_np,
        p_solve_upper_for_preconditioner,
    )
    p_preconditioner_diag_gpu, _, p_preconditioner_diag_transfer = gpu_f64_array(
        p_preconditioner_diag_np,
        "gpu_stages.solve_pEqn.dic_preconditioner_diag",
    )
    p_dic_reciprocal_diag_np = np.divide(
        1.0,
        p_preconditioner_diag_np,
        out=np.ones_like(p_preconditioner_diag_np, dtype=np.float64),
        where=np.abs(p_preconditioner_diag_np) > 1.0e-300,
    )
    p_dic_reciprocal_diag_gpu, _, p_dic_reciprocal_diag_transfer = gpu_f64_array(
        p_dic_reciprocal_diag_np,
        "gpu_stages.solve_pEqn.dic_reciprocal_diag",
    )
    p_solve_passes: list[dict[str, Any]] = []
    p_solve_initial_gpu = p_gpu
    p_solve_initial_np = p_np
    for non_orthogonal_corrector in range(DEFAULT_PRESSURE_NON_ORTHOGONAL_CORRECTORS + 1):
        p_normalization_factor = openfoam_ldu_normalization_factor(
            n_cells,
            owner_np,
            neighbour_np,
            p_solve_diag_for_normalization,
            p_solve_upper_for_preconditioner,
            p_solve_upper_for_preconditioner,
            p_solve_source_for_solver,
            p_solve_initial_np,
        )
        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,
            preconditioner_diag=p_preconditioner_diag_gpu,
            preconditioner_reciprocal_diag=p_dic_reciprocal_diag_gpu,
            dic_schedule=p_dic_level_schedule,
            iterations=DEFAULT_PRESSURE_CG_ITERATIONS,
            residual_tolerance_squared=1.0e-12,
            residual_tolerance=DEFAULT_PRESSURE_CG_RESIDUAL_TOLERANCE,
            min_iterations=DEFAULT_PRESSURE_CG_MIN_ITERATIONS,
            normalization_factor=p_normalization_factor,
        )
        p_solve_passes.append({"non_orthogonal_corrector": non_orthogonal_corrector, **p_solve_performance})
        p_solve_initial_gpu = p_solved_gpu
        p_solve_initial_np = np.asarray(p_solved_gpu.to_numpy(), dtype=np.float64)
    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)
    gpu_rans_turbulence_update(n_cells, nut_gpu, k_gpu, omega_gpu, nut_out_gpu, k_out_gpu, omega_out_gpu)
    if n_omega_wall_cells:
        gpu_rans_omega_wall_update(n_omega_wall_cells, omega_wall_cells_gpu, omega_wall_distances_gpu, laminar_nu, DEFAULT_OMEGA_WALL_BETA1, omega_out_gpu)
    qd.sync()
    kernel_graph_wall_seconds = time.perf_counter() - kernel_graph_started

    materialization_started = time.perf_counter()

    u_diag = np.asarray(u_diag_gpu.to_numpy())
    u_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())
    p_diag = np.asarray(p_diag_gpu.to_numpy())
    p_source = np.asarray(p_source_gpu.to_numpy())
    p_upper = np.asarray(p_upper_gpu.to_numpy())
    u_solved = np.asarray(u_solved_gpu.to_numpy())
    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())
    nut_out = np.asarray(nut_out_gpu.to_numpy())
    k_out = np.asarray(k_out_gpu.to_numpy())
    omega_out = np.asarray(omega_out_gpu.to_numpy())
    materialization_wall_seconds = time.perf_counter() - materialization_started

    UEqn = gpu_matrix_output(
        "UEqn",
        "U",
        u_diag,
        u_matrix_source,
        u_np,
        kernel="gpu_rans_momentum_assembly",
        upper=u_upper,
        upper_kernel="gpu_rans_momentum_diffusion_coefficients",
        lower=u_lower,
        lower_kernel="gpu_rans_momentum_convection_coefficients",
    )
    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_add_vector_field")
    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",
        p_diag,
        p_source,
        p_np,
        kernel="gpu_rans_pressure_assembly",
        upper=p_upper,
        upper_kernel="gpu_rans_pressure_laplacian_coefficients",
    )
    pEqn["source_terms"] = {
        "internal_face_kernel": "gpu_rans_pressure_source_from_flux(phiHbyA)",
        "boundary_face_kernel": "gpu_rans_pressure_source_from_boundary_flux(phiHbyA.boundary)",
        "boundary_laplacian_kernel": None,
        "solver_boundary_diag_kernel": "gpu_rans_pressure_mixed_solver_diag",
        "operator_sign_convention": "openfoam_negative_diag_positive_upper",
        "consistent_rAtU_factor": DEFAULT_SIMPLE_CONSISTENT_RATU_FACTOR,
        "boundary_face_count": n_pressure_boundary_faces,
        "mixed_boundary_face_count": n_pressure_mixed_faces,
    }
    stages = [
        gpu_stage_result(
            "momentum_transport_predict",
            {"case_path": str(case), "solver_name": "quadrants_cuda_rans_solver"},
            kernels=[],
        ),
        gpu_stage_result(
            "assemble_momentum_terms",
            {"terms": [{"name": "gpu_momentum_ddt_diagonal", "kernel_entrypoint": "gpu_rans_momentum_assembly"}, {"name": "gpu_momentum_laminar_turbulent_diffusion", "kernel_entrypoint": "gpu_rans_momentum_diffusion_coefficients_from_face_coeff", "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_entrypoints": ["gpu_rans_momentum_diag_from_offdiag_by_cell", "gpu_rans_momentum_offdiag_abs_accumulate_by_cell", "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_from_face_coeff", "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_diag_from_offdiag_by_cell", "gpu_rans_momentum_offdiag_abs_accumulate_by_cell", "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_entrypoints": ["gpu_rans_momentum_diag_from_offdiag_by_cell", "gpu_rans_momentum_offdiag_abs_accumulate_by_cell", "gpu_rans_momentum_equation_relaxation"]}}, kernels=["gpu_rans_momentum_assembly", "gpu_rans_momentum_diffusion_coefficients_from_face_coeff", "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_diag_from_offdiag_by_cell", "gpu_rans_momentum_offdiag_abs_accumulate_by_cell", "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_add_vector_field"),
                "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_extract_vector_component", "gpu_scatter_vector_component", "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_dilu_forward_level_vector", "gpu_dilu_backward_level_vector", "gpu_bicgstab_update_intermediate_vector_preconditioned", "gpu_bicgstab_update_solution_residual_vector_preconditioned"],
            changed_fields=["U"],
        ),
        gpu_stage_result(
            "compute_pressure_inputs",
            {
                "rAU": gpu_array_output("rAU", rAU, kernel="gpu_rans_pressure_inputs"),
                "rAtU": gpu_array_output("rAtU", rAtU, kernel="gpu_rans_consistent_rAtU"),
                "HbyA": gpu_array_output("HbyA", HbyA, kernel="gpu_rans_momentum_hbyA_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",
                    "kernels": [
                        "gpu_rans_pressure_inputs",
                        "gpu_copy_scalar",
                        "gpu_rans_momentum_h1_finish",
                        "gpu_rans_consistent_rAtU",
                        "gpu_rans_momentum_hbyA_by_cell",
                        "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=OpenFOAM-weighted fvc::flux(HbyA)+OpenFOAM-weighted interpolate(rAtU-rAU)*snGrad(p)*magSf",
                },
            },
            kernels=[
                "gpu_rans_pressure_inputs",
                "gpu_copy_scalar",
                "gpu_rans_momentum_h1_finish",
                "gpu_rans_consistent_rAtU",
                "gpu_rans_momentum_hbyA_by_cell",
                "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_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", "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_only_scalar"),
                "residual": gpu_array_output("p_residual", p_residual, kernel="gpu_pcg_update_solution_residual_only_scalar"),
                "phi": gpu_array_output("phi", phi_solved, kernel="gpu_rans_pressure_flux_correction"),
            },
            kernels=["gpu_rans_negate_scalar_field", "gpu_rans_pressure_mixed_solver_diag", "gpu_ldu_matvec_scalar_symmetric_diag", "gpu_ldu_matvec_scalar_symmetric_face_accumulate", "gpu_pcg_initialize_residual_scalar", "gpu_dic_forward_level_scalar", "gpu_dic_backward_level_scalar", "gpu_cg_dot_scalar", "gpu_pcg_update_solution_residual_only_scalar", "gpu_pcg_update_direction_scalar"],
            changed_fields=["p", "phi"],
        ),
        gpu_stage_result(
            "update_phi_from_pEqn_flux",
            {"phi": gpu_array_output("phi", phi_solved, kernel="gpu_rans_pressure_flux_correction")},
            kernels=["gpu_rans_pressure_flux_correction"],
            changed_fields=["phi"],
        ),
        gpu_stage_result(
            "correct_velocity_pressure_flux",
            {
                "U": gpu_array_output("U", u_final, kernel="gpu_rans_pressure_velocity_correction"),
                "p": gpu_array_output("p", p_final, kernel="gpu_rans_final_correction"),
                "pressure_reference": {"mode": "initial_plus_gpu_correction", "kernel_entrypoint": "gpu_rans_final_correction"},
                "velocity_correction": {"gradient": "internal_face_pressure_jump", "kernel_entrypoint": "gpu_rans_pressure_velocity_correction"},
                "phi": gpu_array_output("phi", phi_solved, kernel="gpu_rans_pressure_flux_correction"),
            },
            kernels=["gpu_rans_final_correction", "gpu_rans_pressure_velocity_correction", "gpu_rans_pressure_flux_correction"],
            changed_fields=["U", "p", "phi"],
        ),
        gpu_stage_result(
            "momentum_transport_correct",
            {
                "U": gpu_array_output("U", u_final, kernel="gpu_rans_pressure_velocity_correction"),
                "p": gpu_array_output("p", p_final, kernel="gpu_rans_final_correction"),
                "phi": gpu_array_output("phi", phi_solved, kernel="gpu_rans_pressure_flux_correction"),
                "nut": gpu_array_output("nut", nut_out, kernel="gpu_rans_turbulence_update"),
                "k": gpu_array_output("k", k_out, kernel="gpu_rans_turbulence_update"),
                "omega": gpu_array_output("omega", omega_out, kernel="gpu_rans_omega_wall_update"),
                "omega_wall_function": {
                    "kernel_entrypoint": "gpu_rans_omega_wall_update",
                    "wall_cell_count": n_omega_wall_cells,
                    "wall_seed_face_count": omega_wall_seed_count,
                    "near_wall_layer_count": omega_wall_layer_count,
                    "wall_distance": "face-normal projection propagated across internal faces",
                    "formula": "max(omega, 6*nu/(beta1*y^2))",
                    "beta1": DEFAULT_OMEGA_WALL_BETA1,
                },
            },
            kernels=["gpu_rans_turbulence_update", "gpu_rans_omega_wall_update"],
            changed_fields=["nut", "k", "omega"],
        ),
    ]
    field_objects = {
        "U": gpu_field_from_source(fields["U"], u_final),
        "p": gpu_field_from_source(fields["p"], p_final),
        "phi": gpu_field_from_source(fields["phi"], phi_solved),
        "nut": gpu_field_from_source(fields["nut"], nut_out),
        "k": gpu_field_from_source(fields["k"], k_out),
        "omega": gpu_field_from_source(fields["omega"], omega_out),
    }
    expected_kernels = sorted({kernel for kernels in GPU_STAGE_KERNEL_ENTRYPOINTS.values() for kernel in kernels})
    gpu_solver_wall_seconds = time.perf_counter() - gpu_solver_started
    profiler_evidence = quadrants_kernel_evidence(expected_kernels)
    timing = {
        "schema_version": 1,
        "clock": "time.perf_counter",
        "gpu_solver_wall_seconds": round(gpu_solver_wall_seconds, 6),
        "input_transfer_wall_seconds": round(input_transfer_wall_seconds, 6),
        "kernel_graph_wall_seconds": round(kernel_graph_wall_seconds, 6),
        "result_materialization_wall_seconds": round(materialization_wall_seconds, 6),
        "device_time_ms_total": profiler_evidence.get("device_time_ms_total"),
        "profile_record_count": profiler_evidence.get("profile_record_count"),
        "timed_scope": "quadrants_cuda_stage_graph_with_input_transfer_and_result_materialization",
    }
    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 {
        "status": "executed",
        "case": case,
        "backend": {
            "selected": backend.get("selected"),
            "provider": backend.get("provider"),
            "device": backend.get("device"),
            "used_cpu_fallback": False,
        },
        "execution_path": "quadrants_cuda_gpu_rans_stage_graph",
        "graph": [stage.name for stage in stages],
        "stages": [transform_summary(stage) for stage in stages],
        "stage_objects": stages,
        "matrices": {"UEqn": UEqn, "pEqn": pEqn},
        "fields": {
            "U": gpu_array_output("U", u_final, kernel="gpu_rans_pressure_velocity_correction"),
            "p": gpu_array_output("p", p_final, kernel="gpu_rans_final_correction"),
            "phi": gpu_array_output("phi", phi_solved, kernel="gpu_rans_pressure_flux_correction"),
            "nut": gpu_array_output("nut", nut_out, kernel="gpu_rans_turbulence_update"),
            "k": gpu_array_output("k", k_out, kernel="gpu_rans_turbulence_update"),
            "omega": gpu_array_output("omega", omega_out, kernel="gpu_rans_omega_wall_update"),
        },
        "transfers": {
            "inputs": [u_transfer, p_transfer, phi_transfer, nut_transfer, k_transfer, omega_transfer, owner_transfer, neighbour_transfer, offdiag_cell_offsets_transfer, offdiag_cell_faces_transfer, offdiag_cell_sides_transfer, sf_transfer, face_centres_transfer, cell_centres_transfer, cell_volumes_transfer, mag_sf_transfer, face_weights_transfer, momentum_diffusion_coeff_transfer, momentum_internal_diag_transfer, momentum_internal_h1_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_pressure_source_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_u_values_transfer, pressure_mixed_face_centres_transfer, pressure_mixed_area_vectors_transfer, pressure_mixed_area_magnitudes_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, u_dilu_reciprocal_diag_transfer, p_preconditioner_diag_transfer, p_dic_reciprocal_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,
                p_diag_alloc,
                p_source_alloc,
                p_upper_alloc,
                u_solved_alloc,
                u_shadow_alloc,
                u_direction_alloc,
                u_operator_direction_alloc,
                u_intermediate_alloc,
                u_operator_intermediate_alloc,
                u_rr_alloc,
                u_rho_alloc,
                u_denominator_alloc,
                u_omega_numerator_alloc,
                u_omega_denominator_alloc,
                u_residual_alloc,
                p_solve_diag_alloc,
                p_solve_source_alloc,
                p_solve_upper_alloc,
                p_solved_alloc,
                p_work_alloc,
                p_residual_alloc,
                p_direction_alloc,
                p_operator_alloc,
                p_rr_alloc,
                p_denominator_alloc,
                phi_solved_alloc,
                u_final_alloc,
                p_final_alloc,
                nut_alloc,
                k_alloc,
                omega_alloc,
            ],
        },
        "profiler": profiler_evidence,
        "timing": timing,
        "diagnostic_artifacts": diagnostic_artifacts,
        "parity_integration": {
            "status": "ready",
            "modes": ["run_one", "split"],
        },
    }



def gpu_solver_input_blocker_summary(gpu_inputs: Mapping[str, Any]) -> dict[str, Any]:
    transfer = gpu_inputs.get("transfer") if isinstance(gpu_inputs.get("transfer"), Mapping) else {}
    required_fields = gpu_inputs.get("required_fields") if isinstance(gpu_inputs.get("required_fields"), Mapping) else {}
    return {
        "status": gpu_inputs.get("status"),
        "array_count": transfer.get("array_count"),
        "transferred_count": transfer.get("transferred_count"),
        "source_nbytes": transfer.get("source_nbytes"),
        "gpu_nbytes": transfer.get("gpu_nbytes"),
        "required_fields": {
            name: {
                "status": value.get("status") if isinstance(value, Mapping) else None,
                "internal_path": value.get("internal_path") if isinstance(value, Mapping) else None,
            }
            for name, value in required_fields.items()
        },
        "matrix_data": gpu_inputs.get("matrix_data"),
    }

def probe_quadrants_cuda_runtime(gpu_present: bool) -> dict[str, Any]:
    try:
        import quadrants as qd
    except Exception as exc:
        raise BackendExecutionError(
            "select_backend",
            "GPU runtime unavailable: failed to import Quadrants",
            details={
                "requested": "gpu",
                "failure_kind": GPU_RUNTIME_UNAVAILABLE,
                "gpu_device_present": gpu_present,
                "provider": GPU_BACKEND_PROVIDER,
                "cause": {"type": type(exc).__name__, "message": str(exc)},
            },
        ) from exc

    try:
        qd.init(arch=qd.cuda, kernel_profiler=True)
    except Exception as exc:
        raise BackendExecutionError(
            "select_backend",
            "GPU runtime unavailable: Quadrants could not initialize CUDA",
            details={
                "requested": "gpu",
                "failure_kind": GPU_RUNTIME_UNAVAILABLE,
                "gpu_device_present": gpu_present,
                "provider": GPU_BACKEND_PROVIDER,
                "device": "cuda",
                "used_cpu_fallback": False,
                "cause": {"type": type(exc).__name__, "message": str(exc)},
            },
        ) from exc

    identity = nvidia_device_identity()
    return {
        "framework": "quadrants",
        "provider": GPU_BACKEND_PROVIDER,
        "device": "cuda",
        "device_kind": "cuda",
        "arch_requested": "cuda",
        "arch_selected": "cuda",
        "gpu_device_present": gpu_present,
        "device_name": identity["device_name"],
        "device_uuid": identity["device_uuid"],
    }


def select_gpu_solver_backend(gpu_present: bool) -> dict[str, Any]:
    runtime = probe_quadrants_cuda_runtime(gpu_present)
    stage_contract = gpu_stage_contracts()
    stage_capabilities = [f"{GPU_STAGE_CAPABILITY_PREFIX}{name}" for name in stage_contract]
    return {
        "requested": "gpu",
        "selected": "gpu",
        "available_backends": ["cpu", "gpu"],
        "used_cpu_fallback": False,
        "cpu_fallback": {
            "allowed": False,
            "rejected_providers": ["foam_stepper_cpu", "openfoam", "host", "cpu"],
        },
        "execution_owner": "gpu_solver_backend_contract",
        "full_solver_guard": FULL_GPU_RANS_GUARD,
        "primitive_evidence": {
            "name": GPU_PRIMITIVE_NAME,
            "path": GPU_PRIMITIVE_PROOF,
            "role": "optional primitive component only",
            "counts_as_full_gpu_rans_solver": False,
        },
        "capabilities": [
            "gpu_runtime:cuda",
            "no_cpu_fallback",
            "gpu_solver_stage_kernels",
            "failure_diagnostics:runtime_unavailable",
            "failure_diagnostics:missing_kernels",
            "failure_diagnostics:unsupported_stage",
            "failure_diagnostics:numerical_mismatch",
            *stage_capabilities,
        ],
        "verifier_integration": {
            "status": "ready",
            "modes": ["run_one", "split"],
        },
        "solver_stage_contract": stage_contract,
        "failure_taxonomy": gpu_backend_failure_taxonomy(),
        **runtime,
    }


def gpu_solver_contract_blocker(backend: Mapping[str, Any], case: Path) -> BackendExecutionError | None:
    stage_contract = backend.get("solver_stage_contract")
    if not isinstance(stage_contract, Mapping):
        return BackendExecutionError(
            "execute_gpu_solver_contract",
            "selected GPU backend has no solver-stage contract",
            details={
                "failure_kind": GPU_STAGE_UNSUPPORTED,
                "backend": backend,
                "case": case,
                "used_cpu_fallback": False,
            },
        )

    unsupported = [
        name
        for name, contract in stage_contract.items()
        if isinstance(contract, Mapping) and contract.get("supported") is False
    ]
    if unsupported:
        return BackendExecutionError(
            "execute_gpu_solver_contract",
            "selected GPU backend does not support required solver stages",
            details={
                "failure_kind": GPU_STAGE_UNSUPPORTED,
                "unsupported_solver_stages": unsupported,
                "backend": backend,
                "case": case,
                "used_cpu_fallback": False,
            },
        )

    missing_kernels = [
        name
        for name, contract in stage_contract.items()
        if isinstance(contract, Mapping) and not contract.get("kernel_entrypoints")
    ]
    if missing_kernels:
        return BackendExecutionError(
            "execute_gpu_solver_contract",
            "selected GPU backend has no registered kernels for required RANS solver stages",
            details={
                "failure_kind": GPU_KERNELS_MISSING,
                "missing_kernel_stages": missing_kernels,
                "backend": backend,
                "case": case,
                "used_cpu_fallback": False,
            },
        )

    integration = backend.get("verifier_integration") if isinstance(backend.get("verifier_integration"), Mapping) else {}
    if integration.get("status") != "ready":
        return BackendExecutionError(
            "execute_gpu_solver_contract",
            "GPU solver stage kernels are registered but not yet wired into verifier parity modes",
            details={
                "failure_kind": GPU_STAGE_UNSUPPORTED,
                "integration_status": integration.get("status"),
                "implemented_kernel_stages": sorted(stage_contract),
                "backend": backend,
                "case": case,
                "used_cpu_fallback": False,
            },
        )

    return None


def gpu_solver_contract_failure(backend: Mapping[str, Any], case: Path) -> BackendExecutionError:
    blocker = gpu_solver_contract_blocker(backend, case)
    if blocker is not None:
        return blocker

    return BackendExecutionError(
        "execute_gpu_solver_contract",
        "selected GPU backend has solver-stage kernels but no solver executor is registered",
        details={
            "failure_kind": GPU_STAGE_UNSUPPORTED,
            "backend": backend,
            "case": case,
            "used_cpu_fallback": False,
        },
    )


def select_execution_backend(requested: str) -> dict[str, Any]:
    if requested not in BACKEND_CHOICES:
        raise BackendExecutionError(
            "select_backend",
            f"unknown backend {requested!r}",
            details={"requested": requested, "choices": list(BACKEND_CHOICES)},
        )

    available = ["cpu"]
    gpu_present = gpu_device_present()
    if requested == "gpu":
        return select_gpu_solver_backend(gpu_present)

    return {
        "requested": requested,
        "selected": "cpu",
        "device": "host",
        "provider": "foam_stepper_cpu",
        "available_backends": available,
        "gpu_device_present": gpu_present,
        "capabilities": [
            "foam_stepper_python_bridge",
            "openfoam_case_loader",
            "oracle_comparison",
        ],
        "used_cpu_fallback": False,
    }


def run_backend_iteration(stepper: Any, backend: Mapping[str, Any], case: Path) -> dict[str, Any]:
    selected = backend.get("selected")
    if selected == "gpu":
        gpu_run = run_gpu_solver_stage_smoke(stepper, backend, case)
        return {
            "backend": dict(backend),
            "case": case,
            "execution_path": "quadrants_cuda_gpu_rans_run_one",
            "result": GpuTransformResult(
                name="run_one_pimple_iteration",
                phase="gpu_solver",
                source=GpuSourceLocation(file=__file__, function="run_backend_iteration"),
                inputs={"case": case},
                outputs={
                    "fields": gpu_run["field_objects"],
                    "graph": gpu_run["stage_objects"],
                    "gpu_solver_stages": gpu_solver_stage_report(gpu_run),
                },
                changed_fields=list(REQUIRED_FIELDS),
                metadata={
                    "backend": "gpu",
                    "execution_path": "quadrants_cuda_gpu_rans_run_one",
                    "used_cpu_fallback": False,
                },
            ),
            "fields": gpu_run["field_objects"],
            "gpu_solver": gpu_solver_stage_report(gpu_run),
        }
    if selected != "cpu":
        raise BackendExecutionError(
            "execute_backend",
            f"backend {selected!r} is not executable",
            details={"backend": backend, "case": case},
        )
    result = stepper.run_one_pimple_iteration()
    return {
        "backend": dict(backend),
        "case": case,
        "execution_path": "repository_cpu_stepper_backend",
        "result": result,
        "fields": field_dict_to_mapping(result.outputs["fields"]),
    }


def run_gpu_split_iteration(foam: Any, stepper: Any, backend: Mapping[str, Any], case: Path, *, 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(
        "split",
        stages,
        evidence={
            "split_step_execution": True,
            "gpu_equivalent": True,
            "backend": backend,
            "execution_path": "quadrants_cuda_gpu_rans_split",
            "turbulence_fields": visible_turbulence_fields(fields),
        },
    )
    solver_state = export_solver_state_checked(foam, stepper, "split_gpu", fields)
    return {
        "fields": fields,
        "state": solver_state,
        "state_summary": foam.describe_solver_state(solver_state),
        "stages": stages,
        "graph": list(gpu_run["graph"]),
        "momentum_terms": ["gpu_identity_mass_momentum"],
        "UEqn": gpu_run["matrices"]["UEqn"],
        "pEqn": gpu_run["matrices"]["pEqn"],
        "matrix_states": gpu_run["matrices"],
        "observability": observability,
        "execution_path": "quadrants_cuda_gpu_rans_split",
        "backend": dict(backend),
        "gpu_solver": gpu_solver_stage_report(gpu_run),
        "diagnostic_artifacts": gpu_run.get("diagnostic_artifacts"),
    }

__all__ = [
    "BackendExecutionError",
    "GPU_INPUT_SCHEMA_VERSION",
    "GPU_NUMERICAL_MISMATCH",
    "prepare_gpu_solver_inputs",
    "gpu_solver_stage_report",
    "gpu_solver_contract_blocker",
    "gpu_solver_contract_failure",
    "gpu_solver_input_blocker_summary",
    "run_backend_iteration",
    "run_gpu_solver_stage_smoke",
    "run_gpu_split_iteration",
    "select_execution_backend",
]