#!/usr/bin/env python3
"""Verifier harness for the migrated AirfRANS OpenFOAM stepper case."""

from __future__ import annotations

import argparse
import dataclasses
import hashlib
import json
import math
import os
import shutil
import subprocess
import sys
import time
import traceback
from collections.abc import Iterable, Mapping
from pathlib import Path
from typing import Any


def _drop_ambient_pythonpath() -> None:
    pythonpath = os.environ.pop("PYTHONPATH", "")
    for entry in pythonpath.split(os.pathsep):
        if not entry:
            continue
        while entry in sys.path:
            sys.path.remove(entry)


_drop_ambient_pythonpath()

import numpy as np

from openfoam_env import apply_openfoam_env, openfoam_env
from prepare_airfrans_stepper_case import (
    DEFAULT_SOURCE,
    REQUIRED_COMPARISON_FIELDS,
    load_case_manifest,
    prepare_case,
    write_case_manifest,
)

ROOT = Path(__file__).resolve().parents[1]
DEFAULT_WORK = ROOT / "tmp/airfrans_stepper_verify"
PRIMARY_FIELDS = ("U", "p", "phi")
TURBULENCE_FIELDS = ("nut", "k", "omega")
REQUIRED_FIELDS = REQUIRED_COMPARISON_FIELDS
REQUIRED_EXECUTION_MODES = ("run_one", "split")
BACKEND_CHOICES = ("auto", "cpu", "gpu")
FULL_GPU_RANS_GUARD = "scripts/verify_gpu_rans_solver.sh"
GPU_PRIMITIVE_PROOF = "scripts/verify_gpu_algorithm.sh"
GPU_PRIMITIVE_NAME = "cell_flux_imbalance"
AIRFRANS_VERIFIER_SCOPE = (
    "AirfRANS/OpenFOAM stepper parity verifier; --backend gpu is a GPU-owned "
    "RANS solver contract with no CPU fallback or primitive-only acceptance"
)

LOADING_FAILURE = "loading_or_parsing_failure"
ORACLE_FAILURE = "openfoam_oracle_failure"
STEPPER_FAILURE = "stepper_execution_failure"
COMPARISON_FAILURE = "numerical_comparison_failure"
BACKEND_FAILURE = "backend_execution_failure"
INTERNAL_FAILURE = "internal_harness_failure"
REQUIRED_EVIDENCE_FAILURE = "required_evidence_failure"

EXIT_CODES = {
    LOADING_FAILURE: 2,
    ORACLE_FAILURE: 3,
    STEPPER_FAILURE: 4,
    COMPARISON_FAILURE: 5,
    BACKEND_FAILURE: 6,
    REQUIRED_EVIDENCE_FAILURE: 8,
    INTERNAL_FAILURE: 7,
}

STAGE_OBSERVABILITY_GROUPS = (
    {
        "name": "momentum_assembly",
        "split_stages": ("assemble_momentum_terms", "assemble_UEqn"),
        "run_one_stages": ("assemble_UEqn",),
        "split_outputs": {
            "assemble_momentum_terms": ("terms",),
            "assemble_UEqn": ("UEqn",),
        },
        "run_one_outputs": {
            "assemble_UEqn": ("UEqn",),
        },
    },
    {
        "name": "pressure_assembly",
        "split_stages": ("compute_pressure_inputs", "assemble_pEqn"),
        "run_one_stages": ("compute_pressure_inputs", "assemble_pEqn"),
        "split_outputs": {
            "compute_pressure_inputs": ("HbyA", "phiHbyA", "rAU", "rAtU"),
            "assemble_pEqn": ("pEqn",),
        },
        "run_one_outputs": {
            "compute_pressure_inputs": ("HbyA", "phiHbyA", "rAU", "rAtU"),
            "assemble_pEqn": ("pEqn",),
        },
    },
    {
        "name": "linear_solve_results",
        "split_stages": ("solve_UEqn", "solve_pEqn"),
        "run_one_stages": ("solve_UEqn", "solve_pEqn"),
        "split_outputs": {
            "solve_UEqn": ("performance", "field_after"),
            "solve_pEqn": ("performance", "p", "phi"),
        },
        "run_one_outputs": {
            "solve_UEqn": ("performance", "field_after"),
            "solve_pEqn": ("performance", "p", "phi"),
        },
    },
    {
        "name": "final_correction",
        "split_stages": ("correct_velocity_pressure_flux",),
        "run_one_stages": ("update_phi_from_pEqn_flux", "correct_velocity_pressure_flux"),
        "split_outputs": {
            "correct_velocity_pressure_flux": ("U", "p", "phi"),
        },
        "run_one_outputs": {
            "update_phi_from_pEqn_flux": ("phi",),
            "correct_velocity_pressure_flux": ("U", "p", "phi"),
        },
    },
    {
        "name": "turbulence_updates",
        "split_stages": ("momentum_transport_predict", "momentum_transport_correct"),
        "run_one_stages": ("momentum_transport_predict", "momentum_transport_correct"),
        "split_outputs": {
            "momentum_transport_predict": ("case_path", "solver_name"),
            "momentum_transport_correct": ("U", "p", "phi", "nut", "k", "omega"),
        },
        "run_one_outputs": {
            "momentum_transport_predict": ("case_path", "solver_name"),
            "momentum_transport_correct": ("U", "p", "phi", "nut", "k", "omega"),
        },
    },
)

STAGE_GROUP_ORDER = {group["name"]: index for index, group in enumerate(STAGE_OBSERVABILITY_GROUPS)}
MODE_COMPARISON_ORDER = {"run_one": 0, "split": 1}

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

INTERMEDIATE_ARTIFACT_FAMILIES = (
    {
        "name": "identity",
        "stage_groups": (),
        "reference_paths": ("case_preparation.prepared_case_identity", "mesh_identity.oracle"),
        "candidate_paths": ("mesh_identity.run_one", "mesh_identity.split"),
        "comparison_paths": ("modes.run_one.mesh_comparison", "modes.split.mesh_comparison"),
        "required": True,
    },
    {
        "name": "pressure_inputs",
        "stage_groups": ("pressure_assembly",),
        "reference_paths": ("modes.reference_split.stages.compute_pressure_inputs",),
        "candidate_paths": ("modes.split.stages.compute_pressure_inputs",),
        "comparison_paths": ("artifact_comparisons.pressure_inputs",),
        "required": True,
    },
    {
        "name": "matrix_operator",
        "stage_groups": ("momentum_assembly", "pressure_assembly"),
        "reference_paths": ("modes.reference_split.momentum_terms", "modes.reference_split.UEqn", "modes.reference_split.stages.relax_UEqn.outputs.UEqn", "modes.reference_split.pEqn", "state_exports.reference_split_matrices"),
        "candidate_paths": ("modes.split.momentum_terms", "modes.split.UEqn", "modes.split.pEqn", "state_exports.split_matrices"),
        "comparison_paths": ("artifact_comparisons.matrix_operator",),
        "required": True,
    },
    {
        "name": "solver",
        "stage_groups": ("linear_solve_results",),
        "reference_paths": ("modes.reference_split.stages.solve_UEqn", "modes.reference_split.stages.solve_pEqn"),
        "candidate_paths": ("modes.split.stages.solve_UEqn", "modes.split.stages.solve_pEqn"),
        "comparison_paths": ("artifact_comparisons.solver",),
        "required": True,
    },
    {
        "name": "correction",
        "stage_groups": ("final_correction",),
        "reference_paths": ("modes.oracle.fields.U", "modes.oracle.fields.p"),
        "candidate_paths": ("modes.split.stages.correct_velocity_pressure_flux", "modes.split.fields.U", "modes.split.fields.p"),
        "comparison_paths": ("modes.split.comparisons.U", "modes.split.comparisons.p"),
        "required": True,
    },
    {
        "name": "flux",
        "stage_groups": ("pressure_assembly", "final_correction"),
        "reference_paths": ("modes.oracle.fields.phi",),
        "candidate_paths": ("modes.split.stages.compute_pressure_inputs", "modes.split.fields.phi"),
        "comparison_paths": ("modes.split.comparisons.phi",),
        "required": True,
    },
    {
        "name": "turbulence",
        "stage_groups": ("turbulence_updates",),
        "reference_paths": ("modes.oracle.fields.nut", "modes.oracle.fields.k", "modes.oracle.fields.omega"),
        "candidate_paths": (
            "modes.split.stages.momentum_transport_predict",
            "modes.split.stages.momentum_transport_correct",
            "modes.split.fields.nut",
            "modes.split.fields.k",
            "modes.split.fields.omega",
        ),
        "comparison_paths": ("modes.split.comparisons.nut", "modes.split.comparisons.k", "modes.split.comparisons.omega"),
        "required": True,
    },
)

ARTIFACT_FAMILY_ORDER = {family["name"]: index for index, family in enumerate(INTERMEDIATE_ARTIFACT_FAMILIES)}

GPU_RUNTIME_UNAVAILABLE = "gpu_runtime_unavailable"
GPU_NUMERICAL_MISMATCH = "gpu_numerical_mismatch"
GPU_INPUT_SCHEMA_VERSION = 1

DIFFERENTIABILITY_VARIABLES = {
    "fields": REQUIRED_FIELDS,
    "mesh_state": ("points", "faces", "owner", "neighbour", "cells", "V", "C", "Cf", "Sf", "magSf", "patches"),
    "case_parameters": ("Uinf", "nu", "rhoInf", "angle_of_attack", "liftDir", "dragDir", "solver_tolerances"),
    "losses": ("per_field_linf_error", "per_field_l2_error", "oracle_parity_allclose"),
}

NON_DIFFERENTIABLE_BOUNDARIES = (
    {
        "name": "openfoam_case_io",
        "category": "io",
        "reason": "OpenFOAM dictionaries, mesh files, and field files are parsed and written as discrete external data.",
    },
    {
        "name": "mesh_topology_and_patch_addressing",
        "category": "topology",
        "reason": "Face/cell connectivity, owner/neighbour addressing, and boundary patch membership are discrete indices.",
    },
    {
        "name": "boundary_condition_selection",
        "category": "discrete_boundary_choice",
        "reason": "Patch types, coupled/constraint flags, and boundary condition dictionaries select code paths.",
    },
    {
        "name": "simple_pimple_and_linear_solver_control",
        "category": "convergence_control",
        "reason": "Iteration counts, stopping criteria, and solver convergence branches are discrete control flow.",
    },
    {
        "name": "turbulence_bounding_and_limiters",
        "category": "limiter",
        "reason": "k/omega/nut bounding and model limiters introduce clipping or branch-dependent updates.",
    },
    {
        "name": "openfoam_cpp_backend_calls",
        "category": "unsupported_solver_operation",
        "reason": "The current registered backend executes OpenFOAM C++ operations and exposes no AD, JVP, or VJP API.",
    },
)

CUSTOM_GRADIENT_REQUIREMENTS = (
    {
        "operation": "sparse_linear_solves",
        "stage_groups": ("linear_solve_results",),
        "reason": "Implicit sparse solves require an adjoint or custom VJP rather than differentiating solver iterations blindly.",
    },
    {
        "operation": "pressure_velocity_correction",
        "stage_groups": ("pressure_assembly", "final_correction"),
        "reason": "The coupled pressure/flux/velocity correction needs a consistent custom gradient for the assembled operators.",
    },
    {
        "operation": "turbulence_model_correctors_and_limiters",
        "stage_groups": ("turbulence_updates",),
        "reason": "Model correctors, wall functions, and bounding need explicit subgradient or smoothing choices.",
    },
    {
        "operation": "mesh_and_boundary_discrete_choices",
        "stage_groups": (),
        "reason": "Topology and patch-type changes are not differentiable; only fixed-topology numeric values can be checked.",
    },
)


class HarnessError(Exception):
    """Categorized verifier failure that can be serialized into the report."""

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

    def to_dict(self) -> dict[str, Any]:
        out: dict[str, Any] = {
            "category": self.category,
            "step": self.step,
            "message": str(self),
            "details": self.details,
        }
        if self.cause is not None:
            out["cause"] = {
                "type": type(self.cause).__name__,
                "message": str(self.cause),
            }
        return json_ready(out)

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

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

    def to_harness_error(self) -> HarnessError:
        return HarnessError(BACKEND_FAILURE, self.step, str(self), details=self.details, cause=self)




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

    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 report_path_value(report: Mapping[str, Any], path: str) -> Any:
    cursor: Any = report
    for part in path.split("."):
        if isinstance(cursor, Mapping):
            cursor = cursor.get(part)
            continue
        if isinstance(cursor, list):
            matched = None
            for item in cursor:
                if isinstance(item, Mapping) and item.get("name") == part:
                    matched = item
                    break
            cursor = matched
            continue
        return None
    return cursor


def report_path_exists(report: Mapping[str, Any], path: str) -> bool:
    value = report_path_value(report, path)
    if value is None:
        return False
    if isinstance(value, Mapping) or isinstance(value, list):
        return bool(value)
    return True



def array_shape(array: Any | None) -> list[int] | None:
    if array is None:
        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 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", {})
        backend = evidence.get("backend", {}) if isinstance(evidence.get("backend"), Mapping) else {}
        if mode == "run_one" and backend.get("selected") != "gpu":
            output_requirements = {}
        output_checks = {}
        for stage_name, required_outputs in output_requirements.items():
            stage = stage_by_name.get(stage_name)
            if stage is None:
                continue
            outputs = stage.get("outputs", {})
            output_keys = set(outputs) if isinstance(outputs, Mapping) else set()
            missing_outputs = [name for name in required_outputs if name not in output_keys]
            output_checks[stage_name] = {
                "required": list(required_outputs),
                "available": sorted(output_keys),
                "missing": missing_outputs,
            }
            if missing_outputs:
                failures.append(
                    {
                        "path": f"stage_observability.{mode}.{group_name}.{stage_name}.outputs",
                        "message": "missing inspectable stage output",
                        "expected": list(required_outputs),
                        "actual": sorted(output_keys),
                    }
                )
        if missing_stages:
            failures.append(
                {
                    "path": f"stage_observability.{mode}.{group_name}.stages",
                    "message": "missing required solver stage",
                    "expected": list(required_stages),
                    "actual": list(stage_by_name),
                }
            )
        groups.append(
            {
                "name": group_name,
                "required_stages": list(required_stages),
                "observed": not missing_stages and not any(check["missing"] for check in output_checks.values()),
                "missing_stages": missing_stages,
                "output_checks": output_checks,
            }
        )
    if failures:
        raise HarnessError(
            STEPPER_FAILURE,
            f"{mode}_stage_observability",
            f"{mode} solver-stage observability is incomplete",
            details={"failures": failures},
        )
    return {
        "mode": mode,
        "graph": [stage["name"] for stage in stages],
        "groups": groups,
        "evidence": json_ready(evidence),
    }


def selected_field_summaries(fields: Mapping[str, Any]) -> dict[str, Any]:
    return {name: field_summary(fields[name]) for name in REQUIRED_FIELDS if name in fields}


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 mesh_patch_table(mesh: Any) -> list[dict[str, Any]]:
    patches = []
    for patch in getattr(mesh, "boundary", []):
        patches.append(
            {
                "index": int(getattr(patch, "index", 0)),
                "name": getattr(patch, "name", ""),
                "type": getattr(patch, "type", ""),
                "start": int(getattr(patch, "start", 0)),
                "size": int(getattr(patch, "size", 0)),
                "coupled": bool(getattr(patch, "coupled", False)),
                "constraint": bool(getattr(patch, "constraint", False)),
            }
        )
    return patches


def mesh_topology_digest(mesh: Any) -> str:
    digest = hashlib.sha256()
    for name in ("n_points", "n_faces", "n_internal_faces", "n_cells"):
        update_hash_text(digest, f"{name}={int(getattr(mesh, name))}")
    update_hash_array(digest, "faces.offsets", mesh.faces.offsets)
    update_hash_array(digest, "faces.values", mesh.faces.values)
    update_hash_array(digest, "cells.offsets", mesh.cells.offsets)
    update_hash_array(digest, "cells.values", mesh.cells.values)
    update_hash_array(digest, "owner", mesh.owner)
    update_hash_array(digest, "neighbour", mesh.neighbour)
    ldu = getattr(mesh, "ldu", {})
    if isinstance(ldu, Mapping):
        if "lower_addr" in ldu:
            update_hash_array(digest, "ldu.lower_addr", ldu["lower_addr"])
        if "upper_addr" in ldu:
            update_hash_array(digest, "ldu.upper_addr", ldu["upper_addr"])
        for entry in ldu.get("patch_addr", []):
            update_hash_text(digest, f"ldu.patch_index={entry.get('patch_index')}")
            update_hash_array(digest, "ldu.patch_addr", entry.get("addr", []))
    for patch in mesh_patch_table(mesh):
        update_hash_text(digest, json.dumps(patch, sort_keys=True))
    return digest.hexdigest()


def mesh_geometry_digest(mesh: Any) -> str:
    digest = hashlib.sha256()
    for name in ("points", "V", "C", "Cf", "Sf", "magSf"):
        update_hash_array(digest, name, getattr(mesh, name))
    return digest.hexdigest()


def mesh_identity(mesh: Any) -> dict[str, Any]:
    return {
        "n_points": int(mesh.n_points),
        "n_faces": int(mesh.n_faces),
        "n_internal_faces": int(mesh.n_internal_faces),
        "n_cells": int(mesh.n_cells),
        "patches": mesh_patch_table(mesh),
        "topology_sha256": mesh_topology_digest(mesh),
        "geometry_sha256": mesh_geometry_digest(mesh),
    }


def compare_mesh_identity(mode: str, actual: Mapping[str, Any], expected: Mapping[str, Any]) -> tuple[dict[str, Any], list[dict[str, Any]]]:
    keys = ("n_points", "n_faces", "n_internal_faces", "n_cells", "topology_sha256", "geometry_sha256", "patches")
    differences = {
        key: {"actual": actual.get(key), "expected": expected.get(key)}
        for key in keys
        if actual.get(key) != expected.get(key)
    }
    report = {"matches_oracle": not differences, "differences": differences}
    if not differences:
        return report, []
    return report, [{"mode": mode, "kind": "mesh_identity", "differences": differences}]


def field_comparison_attribution(
    mode: str,
    field: str,
    reason: str,
    stage_graph: Iterable[str] | None,
) -> dict[str, Any]:
    observed_stages = set(stage_graph or ())
    candidates = []
    for group_name in FIELD_COMPARISON_ATTRIBUTION.get(field, ()):
        group = next(item for item in STAGE_OBSERVABILITY_GROUPS if item["name"] == group_name)
        stages = list(group.get(f"{mode}_stages", ()))
        present_stages = [stage for stage in stages if not observed_stages or stage in observed_stages]
        candidates.append(
            {
                "group": group_name,
                "stages": present_stages,
                "evidence_path": f"stage_observability.{mode}.{group_name}",
            }
        )

    likely = candidates[0] if candidates else {"group": "unknown", "stages": [], "evidence_path": None}
    out = {
        "field": field,
        "reason": reason,
        "likely_stage_group": likely["group"],
        "likely_stages": likely["stages"],
        "candidate_stage_groups": candidates,
        "evidence_path": likely["evidence_path"],
        "basis": "field-to-stage dependency map for the selected solver graph",
    }
    if reason in {"missing_field", "shape_mismatch"}:
        out["precondition_failure"] = "field_export_or_state_mapping"
        out["field_evidence_path"] = f"modes.{mode}.fields.{field}"
        out["basis"] = "field is absent or has the wrong topology; inspect export/state mapping first, then the mapped solver stages"
    return out


def safe_array_item(array: Any, index: int) -> Any:
    arr = np.asarray(array)
    if index < 0 or index >= arr.shape[0]:
        return None
    value = arr[index]
    if isinstance(value, np.ndarray):
        return value.tolist()
    if isinstance(value, np.generic):
        return value.item()
    return json_ready(value)


def cell_neighbourhood_context(mesh: Any, cell_index: int) -> dict[str, Any]:
    owner = np.asarray(getattr(mesh, "owner", []))
    neighbour = np.asarray(getattr(mesh, "neighbour", []))
    owner_faces = np.flatnonzero(owner == cell_index)
    neighbour_faces = np.flatnonzero(neighbour == cell_index)
    touching_faces = sorted({int(face) for face in np.concatenate((owner_faces, neighbour_faces))})[:12]
    boundary_patches = []
    for patch in getattr(mesh, "boundary", []):
        face_cells = np.asarray(getattr(patch, "face_cells", []))
        hits = np.flatnonzero(face_cells == cell_index)
        if hits.size:
            boundary_patches.append(
                {
                    "patch_index": int(getattr(patch, "index", 0)),
                    "patch_name": getattr(patch, "name", ""),
                    "patch_type": getattr(patch, "type", ""),
                    "patch_local_entries": [int(item) for item in hits[:8]],
                    "face_indices": [int(np.asarray(getattr(patch, "face_indices", []))[item]) for item in hits[:8]],
                }
            )
    return {
        "entity": "cell",
        "cell_index": cell_index,
        "volume": safe_array_item(getattr(mesh, "V", []), cell_index),
        "centre": safe_array_item(getattr(mesh, "C", []), cell_index),
        "touching_internal_faces": touching_faces,
        "touching_internal_face_count": int(owner_faces.size + neighbour_faces.size),
        "boundary_patches": boundary_patches,
    }


def face_neighbourhood_context(mesh: Any, face_index: int) -> dict[str, Any]:
    owner = np.asarray(getattr(mesh, "owner", []))
    neighbour = np.asarray(getattr(mesh, "neighbour", []))
    context = {
        "entity": "face",
        "face_index": face_index,
        "owner": int(owner[face_index]) if 0 <= face_index < owner.shape[0] else None,
        "neighbour": int(neighbour[face_index]) if 0 <= face_index < neighbour.shape[0] else None,
        "centre": safe_array_item(getattr(mesh, "Cf", []), face_index),
        "area_vector": safe_array_item(getattr(mesh, "Sf", []), face_index),
        "area_magnitude": safe_array_item(getattr(mesh, "magSf", []), face_index),
        "patch": None,
    }
    for patch in getattr(mesh, "boundary", []):
        face_indices = np.asarray(getattr(patch, "face_indices", []))
        hits = np.flatnonzero(face_indices == face_index)
        if hits.size:
            context["patch"] = {
                "patch_index": int(getattr(patch, "index", 0)),
                "patch_name": getattr(patch, "name", ""),
                "patch_type": getattr(patch, "type", ""),
                "patch_local_entry": int(hits[0]),
                "face_cell": safe_array_item(getattr(patch, "face_cells", []), int(hits[0])),
            }
            break
    return context


def local_entity_context(mesh: Any | None, location: Mapping[str, Any] | None, field_name: str) -> dict[str, Any] | None:
    if mesh is None or not location:
        return None
    entity_kind = location.get("entity_kind")
    entity_index = location.get("entity_index")
    if not isinstance(entity_index, int):
        return None
    if entity_kind == "cell":
        context = cell_neighbourhood_context(mesh, entity_index)
    elif entity_kind in {"face", "internal_face"}:
        context = face_neighbourhood_context(mesh, entity_index)
    else:
        context = {"entity": entity_kind, "entity_index": entity_index}
    context["field"] = field_name
    context["array_index"] = location.get("array_index")
    context["component_index"] = location.get("component_index")
    return context


NUMERIC_ARTIFACT_SCHEMA_VERSION = 1
NUMERIC_ARTIFACT_TOP_N = 5


def numeric_artifact_array(value: Any) -> np.ndarray | None:
    if value is None:
        return None
    arr = np.asarray(value)
    if arr.dtype == object or not np.issubdtype(arr.dtype, np.number):
        return None
    return np.ascontiguousarray(arr)


def add_numeric_artifact(arrays: dict[str, np.ndarray], name: str, value: Any) -> None:
    arr = numeric_artifact_array(value)
    if arr is not None:
        arrays[name] = arr


def add_field_numeric_artifact(arrays: dict[str, np.ndarray], name: str, field: Any) -> None:
    if field is None:
        return
    if isinstance(field, Mapping):
        add_numeric_artifact(arrays, name, field.get("internal"))
    elif hasattr(field, "internal"):
        add_numeric_artifact(arrays, name, getattr(field, "internal"))


def _value_from_mapping_or_attr(value: Any, name: str) -> Any:
    if isinstance(value, Mapping):
        return value.get(name)
    return getattr(value, name, None)


def add_matrix_numeric_artifacts(arrays: dict[str, np.ndarray], prefix: str, matrix: Any) -> None:
    if matrix is None:
        return
    for component in ("diag", "upper", "lower", "source"):
        value = _value_from_mapping_or_attr(matrix, component)
        if value is not None:
            add_numeric_artifact(arrays, f"{prefix}.{component}", value)
    psi = _value_from_mapping_or_attr(matrix, "psi")
    if psi is not None:
        add_field_numeric_artifact(arrays, f"{prefix}.psi", psi)


def add_momentum_term_source_artifacts(arrays: dict[str, np.ndarray], momentum_terms: Any) -> None:
    term_outputs = stage_output(momentum_terms, "terms")
    if term_outputs is None:
        return
    source_names = {
        "assemble_momentum_ddt": "matrix_terms.UEqn.ddt.source",
        "assemble_momentum_div_phi_U": "matrix_terms.UEqn.div.source",
        "assemble_momentum_divDevSigma": "matrix_terms.UEqn.divDevSigma.source",
        "assemble_momentum_sources": "matrix_terms.UEqn.fvModels.source",
    }
    for term in term_outputs:
        term_name = _value_from_mapping_or_attr(term, "name")
        if term_name in source_names:
            matrix = _value_from_mapping_or_attr(term, "matrix")
            add_numeric_artifact(arrays, source_names[term_name], _value_from_mapping_or_attr(matrix, "source"))
            if term_name == "assemble_momentum_divDevSigma":
                diagnostics = _value_from_mapping_or_attr(term, "diagnostics")
                add_numeric_artifact(arrays, "matrix_terms.UEqn.divDevSigma.wall_dev_tau.source", _value_from_mapping_or_attr(diagnostics, "wall_dev_tau_source"))
        elif term_name == "assemble_momentum_MRF_DDt":
            add_field_numeric_artifact(arrays, "matrix_terms.UEqn.MRF_DDt.field", _value_from_mapping_or_attr(term, "field"))


def write_numeric_artifact_file(path: Path, arrays: Mapping[str, np.ndarray], *, role: str) -> dict[str, Any]:
    materialized = {name: np.ascontiguousarray(value) for name, value in arrays.items()}
    path.parent.mkdir(parents=True, exist_ok=True)
    np.savez_compressed(path, **materialized)
    return {
        "schema_version": NUMERIC_ARTIFACT_SCHEMA_VERSION,
        "role": role,
        "path": str(path),
        "format": "npz",
        "artifact_count": len(materialized),
        "artifacts": {name: array_stats(value) for name, value in materialized.items()},
    }


def stage_output(stage: Any, name: str) -> Any:
    outputs = getattr(stage, "outputs", None)
    if isinstance(outputs, Mapping):
        return outputs.get(name)
    return None


def write_split_numeric_artifacts(
    path: Path | None,
    *,
    role: str,
    momentum_terms: Any | None,
    assemble_UEqn: Any,
    relax_UEqn: Any,
    solve_UEqn: Any,
    compute_pressure_inputs: Any,
    assemble_pEqn: Any,
    solve_pEqn: Any,
    correct_velocity_pressure_flux: Any,
    momentum_transport_correct: Any,
    fields: Mapping[str, Any],
) -> dict[str, Any] | None:
    if path is None:
        return None

    arrays: dict[str, np.ndarray] = {}
    add_momentum_term_source_artifacts(arrays, momentum_terms)
    unrelaxed_UEqn = stage_output(assemble_UEqn, "UEqn")
    if unrelaxed_UEqn is not None:
        diag = getattr(unrelaxed_UEqn, "diag", None)
        if diag is not None:
            add_numeric_artifact(arrays, "matrix_operator.UEqn.unrelaxed_diag", diag)
        source = getattr(unrelaxed_UEqn, "source", None)
        if source is not None:
            add_numeric_artifact(arrays, "matrix_operator.UEqn.unrelaxed_source", source)
    add_matrix_numeric_artifacts(arrays, "matrix_operator.UEqn", stage_output(relax_UEqn, "UEqn"))
    add_matrix_numeric_artifacts(arrays, "matrix_operator.pEqn", stage_output(assemble_pEqn, "pEqn"))

    for name in ("HbyA", "phiHbyA", "rAU", "rAtU"):
        add_field_numeric_artifact(arrays, f"pressure_inputs.{name}", stage_output(compute_pressure_inputs, name))

    add_field_numeric_artifact(arrays, "solver.solve_UEqn.field_before", stage_output(solve_UEqn, "field_before"))
    add_field_numeric_artifact(arrays, "solver.solve_UEqn.field_after", stage_output(solve_UEqn, "field_after"))
    add_matrix_numeric_artifacts(arrays, "solver.solve_UEqn.matrix_before", stage_output(solve_UEqn, "matrix_before"))
    add_numeric_artifact(arrays, "solver.solve_UEqn.solve_diag", stage_output(solve_UEqn, "solve_diag"))
    add_numeric_artifact(arrays, "solver.solve_UEqn.solve_source", stage_output(solve_UEqn, "solve_source"))
    add_numeric_artifact(arrays, "solver.solve_UEqn.initial_residual", stage_output(solve_UEqn, "initial_residual"))
    add_numeric_artifact(arrays, "solver.solve_UEqn.preconditioned_residual", stage_output(solve_UEqn, "preconditioned_residual"))
    add_numeric_artifact(arrays, "solver.solve_UEqn.level_preconditioned_residual", stage_output(solve_UEqn, "level_preconditioned_residual"))
    add_numeric_artifact(arrays, "solver.solve_UEqn.operator_preconditioned_direction", stage_output(solve_UEqn, "operator_preconditioned_direction"))
    add_numeric_artifact(arrays, "solver.solve_UEqn.first_iteration.rho", stage_output(solve_UEqn, "first_iteration_rho"))
    add_numeric_artifact(arrays, "solver.solve_UEqn.first_iteration.denominator", stage_output(solve_UEqn, "first_iteration_denominator"))
    add_numeric_artifact(arrays, "solver.solve_UEqn.first_iteration.alpha", stage_output(solve_UEqn, "first_iteration_alpha"))
    add_numeric_artifact(arrays, "solver.solve_UEqn.first_iteration.intermediate_residual", stage_output(solve_UEqn, "first_iteration_intermediate_residual"))
    add_numeric_artifact(arrays, "solver.solve_UEqn.first_iteration.second_preconditioned_residual", stage_output(solve_UEqn, "first_iteration_second_preconditioned_residual"))
    add_numeric_artifact(arrays, "solver.solve_UEqn.first_iteration.operator_second_preconditioned_residual", stage_output(solve_UEqn, "first_iteration_operator_second_preconditioned_residual"))
    add_numeric_artifact(arrays, "solver.solve_UEqn.first_iteration.omega_numerator", stage_output(solve_UEqn, "first_iteration_omega_numerator"))
    add_numeric_artifact(arrays, "solver.solve_UEqn.first_iteration.omega_denominator", stage_output(solve_UEqn, "first_iteration_omega_denominator"))
    add_numeric_artifact(arrays, "solver.solve_UEqn.first_iteration.omega", stage_output(solve_UEqn, "first_iteration_omega"))
    add_numeric_artifact(arrays, "solver.solve_UEqn.first_iteration.residual_after_omega", stage_output(solve_UEqn, "first_iteration_residual_after_omega"))
    add_numeric_artifact(arrays, "solver.solve_UEqn.first_iteration.solution_after_omega", stage_output(solve_UEqn, "first_iteration_solution_after_omega"))
    add_field_numeric_artifact(arrays, "solver.solve_pEqn.field_before", stage_output(solve_pEqn, "field_before"))
    add_field_numeric_artifact(arrays, "solver.solve_pEqn.p", stage_output(solve_pEqn, "p"))
    add_field_numeric_artifact(arrays, "solver.solve_pEqn.phi", stage_output(solve_pEqn, "phi"))
    add_matrix_numeric_artifacts(arrays, "solver.solve_pEqn.matrix_before", stage_output(solve_pEqn, "matrix_before"))

    for name in ("U", "p", "phi"):
        add_field_numeric_artifact(arrays, f"final_correction.{name}", stage_output(correct_velocity_pressure_flux, name))
    for name in ("nut", "k", "omega"):
        add_field_numeric_artifact(arrays, f"turbulence.{name}", stage_output(momentum_transport_correct, name))
    for name in REQUIRED_FIELDS:
        add_field_numeric_artifact(arrays, f"fields.{name}", fields.get(name))

    return write_numeric_artifact_file(path, arrays, role=role)


def load_numeric_artifact_file(path: Path | None) -> dict[str, np.ndarray] | None:
    if path is None or not path.exists():
        return None
    try:
        with np.load(path, allow_pickle=False) as data:
            return {name: np.asarray(data[name]) for name in data.files}
    except Exception:
        return None


def diagnostic_artifact_path(report: Mapping[str, Any], role: str) -> Path | None:
    value = report_path_value(report, f"diagnostic_artifacts.{role}.path")
    if value is None:
        return None
    return Path(str(value))

def build_differential_trace(report: Mapping[str, Any], *, mesh_context: Any | None = None) -> dict[str, Any]:
    from foam_stepper.differential_trace import build_differential_trace_report, load_npz_artifacts

    reference_path = diagnostic_artifact_path(report, "reference_split")
    candidate_path = diagnostic_artifact_path(report, "split")
    reference_artifacts = load_npz_artifacts(reference_path)
    candidate_artifacts = load_npz_artifacts(candidate_path)
    tolerances = report.get("tolerances", {}) if isinstance(report.get("tolerances"), Mapping) else {}
    rtol = float(tolerances.get("rtol", 0.0) or 0.0)
    atol = float(tolerances.get("atol", 0.0) or 0.0)

    def context_for(name: str, index: tuple[int, ...], shape: tuple[int, ...]) -> dict[str, Any] | None:
        location = artifact_location(name, index, shape, mesh_context)
        return local_entity_context(mesh_context, location, name)

    substitution = report.get("differential_trace", {}).get("substitution", {}) if isinstance(report.get("differential_trace"), Mapping) else {}
    return build_differential_trace_report(
        reference_artifacts=reference_artifacts,
        candidate_artifacts=candidate_artifacts,
        reference_path=reference_path,
        candidate_path=candidate_path,
        rtol=rtol,
        atol=atol,
        top_n=NUMERIC_ARTIFACT_TOP_N,
        location_context=context_for,
        substitution=substitution,
    )


def load_trace_substitutions(report: Mapping[str, Any], checkpoints: Iterable[str]) -> tuple[dict[str, np.ndarray], dict[str, Any]]:
    from foam_stepper.differential_trace import extract_reference_substitutions, load_npz_artifacts

    reference_path = diagnostic_artifact_path(report, "reference_split")
    reference_artifacts = load_npz_artifacts(reference_path)
    if reference_artifacts is None:
        requested = list(dict.fromkeys(str(item) for item in checkpoints))
        return {}, {
            "schema_version": 1,
            "status": "failed" if requested else "not_requested",
            "requested": requested,
            "loaded": [],
            "missing_reference": requested,
            "unsupported": [],
            "source_artifact": str(reference_path) if reference_path is not None else None,
        }
    substitutions, manifest = extract_reference_substitutions(reference_artifacts, checkpoints)
    manifest["source_artifact"] = str(reference_path) if reference_path is not None else None
    return substitutions, manifest


def artifact_entity_kind(mesh: Any | None, shape: tuple[int, ...]) -> str | None:
    if mesh is None or not shape:
        return None
    first_dim = int(shape[0])
    if first_dim == int(getattr(mesh, "n_cells", -1)):
        return "cell"
    if first_dim == int(getattr(mesh, "n_internal_faces", -1)):
        return "internal_face"
    return None


def artifact_location(name: str, index: tuple[int, ...], shape: tuple[int, ...], mesh: Any | None) -> dict[str, Any]:
    entity_kind = artifact_entity_kind(mesh, shape) or "array"
    entity_index = int(index[0]) if index else None
    component_index = list(index[1:]) if len(index) > 1 else None
    return {
        "array_index": list(index),
        "entity_kind": entity_kind,
        "entity_index": entity_index,
        "component_index": component_index,
        "artifact": name,
    }


def top_numeric_differences(
    name: str,
    actual: np.ndarray,
    expected: np.ndarray,
    diff: np.ndarray,
    *,
    mesh_context: Any | None,
    limit: int = NUMERIC_ARTIFACT_TOP_N,
) -> list[dict[str, Any]]:
    if diff.size == 0:
        return []
    flat = diff.reshape(-1)
    finite = np.isfinite(flat)
    if np.any(finite):
        finite_indices = np.flatnonzero(finite)
        finite_values = flat[finite_indices]
        if finite_values.size > limit:
            local = np.argpartition(finite_values, -limit)[-limit:]
            candidate_indices = finite_indices[local]
        else:
            candidate_indices = finite_indices
        ordered = sorted(candidate_indices, key=lambda item: (-float(flat[int(item)]), int(item)))
    else:
        ordered = [int(item) for item in np.flatnonzero(~finite)[:limit]]

    out = []
    for flat_index in ordered[:limit]:
        index = tuple(int(item) for item in np.unravel_index(int(flat_index), diff.shape))
        location = artifact_location(name, index, tuple(diff.shape), mesh_context)
        out.append(
            {
                "location": location,
                "abs_difference": finite_float(diff[index]),
                "actual_value": value_at(actual, index),
                "reference_value": value_at(expected, index),
                "local_context": local_entity_context(mesh_context, location, name),
            }
        )
    return out


def compare_numeric_artifact_arrays(
    name: str,
    key: str,
    reference: np.ndarray | None,
    candidate: np.ndarray | None,
    *,
    rtol: float,
    atol: float,
    mesh_context: Any | None,
) -> dict[str, Any]:
    out: dict[str, Any] = {
        "name": name,
        "artifact_key": key,
        "reference_available": reference is not None,
        "candidate_available": candidate is not None,
        "comparison": "numpy.allclose plus absolute-difference diagnostics",
        "rtol": rtol,
        "atol": atol,
    }
    if reference is None or candidate is None:
        return {**out, "shape_matches": False, "allclose": False, "reason": "missing_numeric_artifact"}

    ref = np.asarray(reference)
    actual = np.asarray(candidate)
    shape_matches = ref.shape == actual.shape
    dtype_matches = ref.dtype == actual.dtype
    out.update(
        {
            "reference_shape": array_shape(ref),
            "candidate_shape": array_shape(actual),
            "shape_matches": shape_matches,
            "reference_dtype": str(ref.dtype),
            "candidate_dtype": str(actual.dtype),
            "dtype_matches": dtype_matches,
        }
    )
    if not shape_matches:
        return {**out, "allclose": False, "reason": "shape_mismatch"}
    if ref.size == 0:
        return {
            **out,
            "allclose": True,
            "reason": None,
            "max_abs": 0.0,
            "mean_abs": 0.0,
            "rms_abs": 0.0,
            "largest_difference": None,
            "top_differences": [],
        }

    diff = np.abs(actual - ref)
    finite = np.isfinite(diff)
    finite_diff = diff[finite]
    if finite_diff.size:
        max_flat = int(np.argmax(np.where(finite, diff, -np.inf)))
        max_abs = finite_float(diff.reshape(-1)[max_flat])
        mean_abs = finite_float(np.mean(finite_diff))
        rms_abs = finite_float(np.sqrt(np.mean(np.square(finite_diff))))
    else:
        max_flat = int(np.flatnonzero(~finite.reshape(-1))[0])
        max_abs = None
        mean_abs = None
        rms_abs = None

    max_index = tuple(int(item) for item in np.unravel_index(max_flat, diff.shape))
    location = artifact_location(name, max_index, tuple(diff.shape), mesh_context)
    expected_at_max = value_at(ref, max_index)
    tolerance_at_max = atol + rtol * abs(float(expected_at_max)) if isinstance(expected_at_max, (int, float)) else None
    allclose = bool(np.allclose(actual, ref, rtol=rtol, atol=atol, equal_nan=False))
    out.update(
        {
            "allclose": allclose,
            "reason": None if allclose else "value_mismatch",
            "max_abs": max_abs,
            "mean_abs": mean_abs,
            "rms_abs": rms_abs,
            "nonfinite_error_count": int(diff.size - np.count_nonzero(finite)),
            "largest_difference": {
                "location": location,
                "actual_value": value_at(actual, max_index),
                "reference_value": expected_at_max,
                "actual_entity_value": entity_value_at(actual, location.get("entity_index")),
                "reference_entity_value": entity_value_at(ref, location.get("entity_index")),
                "tolerance_at_max": finite_float(tolerance_at_max),
                "local_context": local_entity_context(mesh_context, location, name),
            },
            "top_differences": top_numeric_differences(name, actual, ref, diff, mesh_context=mesh_context),
        }
    )
    return out

def field_compare_report(
    name: str,
    actual_field: Any,
    expected_field: Any,
    *,
    rtol: float,
    atol: float,
    mesh_context: Any | None = None,
) -> tuple[dict[str, Any], dict[str, Any] | None]:
    actual = np.asarray(actual_field.internal)
    expected = np.asarray(expected_field.internal)
    actual_entity_kind = getattr(actual_field, "entity_kind", "")
    expected_entity_kind = getattr(expected_field, "entity_kind", "")
    report: dict[str, Any] = {
        "field": name,
        "actual_shape": array_shape(actual),
        "expected_shape": array_shape(expected),
        "shape_matches": actual.shape == expected.shape,
        "actual_entity_kind": actual_entity_kind,
        "expected_entity_kind": expected_entity_kind,
        "entity_kind_matches": actual_entity_kind == expected_entity_kind,
        "rtol": rtol,
        "atol": atol,
        "comparison": "numpy.allclose(equal_nan=False)",
    }

    if actual.shape != expected.shape:
        report.update({"allclose": False, "reason": "shape_mismatch"})
        return report, {"field": name, "reason": "shape_mismatch", **report}

    if actual.size == 0:
        report.update(
            {
                "allclose": True,
                "max_abs": 0.0,
                "mean_abs": 0.0,
                "rms_abs": 0.0,
                "location": None,
                "actual_at_max": None,
                "expected_at_max": None,
                "top_differences": [],
            }
        )
        return report, None

    diff = np.abs(actual - expected)
    finite = np.isfinite(diff)
    if not np.all(finite):
        flat_index = int(np.flatnonzero(~finite)[0])
        max_abs: float | None = None
    else:
        flat_index = int(np.argmax(diff))
        max_abs = finite_float(diff.reshape(-1)[flat_index])

    max_index = tuple(int(item) for item in np.unravel_index(flat_index, diff.shape))
    entity_index = max_index[0] if max_index else None
    component_index = list(max_index[1:]) if len(max_index) > 1 else None
    actual_at_max = value_at(actual, max_index)
    expected_at_max = value_at(expected, max_index)
    tolerance_at_max = None
    if isinstance(expected_at_max, (int, float)):
        tolerance_at_max = atol + rtol * abs(float(expected_at_max))

    finite_diff = diff[finite]
    allclose = bool(np.allclose(actual, expected, rtol=rtol, atol=atol, equal_nan=False))
    location = {
        "array_index": list(max_index),
        "entity_kind": getattr(actual_field, "entity_kind", ""),
        "entity_index": entity_index,
        "component_index": component_index,
    }
    report.update(
        {
            "allclose": allclose,
            "max_abs": max_abs,
            "mean_abs": finite_float(np.mean(finite_diff)) if finite_diff.size else None,
            "rms_abs": finite_float(np.sqrt(np.mean(np.square(finite_diff)))) if finite_diff.size else None,
            "location": location,
            "actual_at_max": actual_at_max,
            "expected_at_max": expected_at_max,
            "actual_entity_at_max": entity_value_at(actual, entity_index),
            "expected_entity_at_max": entity_value_at(expected, entity_index),
            "local_context": local_entity_context(mesh_context, location, name),
            "tolerance_at_max": finite_float(tolerance_at_max),
            "nonfinite_error_count": int(diff.size - np.count_nonzero(finite)),
            "top_differences": top_numeric_differences(name, actual, expected, diff, mesh_context=mesh_context),
        }
    )
    if allclose:
        return report, None
    report["reason"] = "value_mismatch"
    return report, {"field": name, "reason": "value_mismatch", **report}


def compare_fields(
    mode: str,
    actual: Mapping[str, Any],
    expected: Mapping[str, Any],
    *,
    rtol: float,
    atol: float,
    stage_graph: Iterable[str] | None = None,
    mesh_context: Any | None = None,
) -> tuple[dict[str, Any], list[dict[str, Any]]]:
    report: dict[str, Any] = {}
    mismatches: list[dict[str, Any]] = []
    for name in REQUIRED_FIELDS:
        if name not in actual or name not in expected:
            missing = {
                "field": name,
                "reason": "missing_field",
                "missing_actual": name not in actual,
                "missing_expected": name not in expected,
                "actual_fields": sorted(actual),
                "expected_fields": sorted(expected),
            }
            attribution = field_comparison_attribution(mode, name, "missing_field", stage_graph)
            report[name] = {"field": name, "allclose": False, "attribution": attribution, **missing}
            mismatches.append({"mode": mode, "attribution": attribution, **missing})
            continue
        field_report, mismatch = field_compare_report(name, actual[name], expected[name], rtol=rtol, atol=atol, mesh_context=mesh_context)
        report[name] = field_report
        if mismatch is not None:
            attribution = field_comparison_attribution(mode, name, mismatch["reason"], stage_graph)
            field_report["attribution"] = attribution
            mismatch["attribution"] = attribution
            mismatches.append({"mode": mode, **mismatch})
    return report, mismatches


def mismatch_sort_key(mismatch: Mapping[str, Any]) -> tuple[int, int, int]:
    attribution = mismatch.get("attribution") if isinstance(mismatch.get("attribution"), Mapping) else {}
    stage_group = attribution.get("likely_stage_group")
    field = str(mismatch.get("field", ""))
    try:
        field_index = REQUIRED_FIELDS.index(field)
    except ValueError:
        field_index = len(REQUIRED_FIELDS)
    return (
        MODE_COMPARISON_ORDER.get(str(mismatch.get("mode", "")), len(MODE_COMPARISON_ORDER)),
        STAGE_GROUP_ORDER.get(str(stage_group), len(STAGE_GROUP_ORDER)),
        field_index,
    )


def field_failure_context(mismatch: Mapping[str, Any]) -> dict[str, Any]:
    return {
        "field": mismatch.get("field"),
        "reason": mismatch.get("reason"),
        "shape": {
            "actual": mismatch.get("actual_shape"),
            "expected": mismatch.get("expected_shape"),
            "matches": mismatch.get("shape_matches"),
        },
        "entity_kind": {
            "actual": mismatch.get("actual_entity_kind"),
            "expected": mismatch.get("expected_entity_kind"),
            "matches": mismatch.get("entity_kind_matches"),
        },
        "tolerance": {
            "rtol": mismatch.get("rtol"),
            "atol": mismatch.get("atol"),
            "at_largest_difference": mismatch.get("tolerance_at_max"),
        },
        "error": {
            "max_abs": mismatch.get("max_abs"),
            "mean_abs": mismatch.get("mean_abs"),
            "rms_abs": mismatch.get("rms_abs"),
            "nonfinite_error_count": mismatch.get("nonfinite_error_count"),
        },
        "largest_difference": {
            "location": mismatch.get("location"),
            "actual_value": mismatch.get("actual_at_max"),
            "expected_value": mismatch.get("expected_at_max"),
            "actual_entity_value": mismatch.get("actual_entity_at_max"),
            "expected_entity_value": mismatch.get("expected_entity_at_max"),
        },
        "local_context": mismatch.get("local_context"),
        "top_differences": mismatch.get("top_differences"),
    }

def stage_group_observability(report: Mapping[str, Any], mode: Any, group_name: Any) -> dict[str, Any] | None:
    observability = report.get("stage_observability", {}).get(mode, {})
    groups = observability.get("groups", []) if isinstance(observability, Mapping) else []
    for group in groups:
        if isinstance(group, Mapping) and group.get("name") == group_name:
            return dict(group)
    return None


def finest_supported_target(mismatch: Mapping[str, Any], attribution: Mapping[str, Any], report: Mapping[str, Any]) -> dict[str, Any]:
    mode = mismatch.get("mode")
    stage_group = attribution.get("likely_stage_group", "unknown")
    stages = list(attribution.get("likely_stages") or [])
    observability = stage_group_observability(report, mode, stage_group)
    output_targets: list[dict[str, Any]] = []
    if mismatch.get("kind") == "mesh_identity":
        return {
            "kind": "identity_artifact_mismatch",
            "single_next_target": f"{mode}.mesh_identity",
            "finest_granularity": "mesh identity digest/patch table",
            "artifact_family": "mesh_identity",
            "evidence_path": f"modes.{mode}.mesh_comparison",
            "differences": mismatch.get("differences"),
            "why": "mesh identity failed before field or solver-stage comparisons, so all downstream numerical evidence is untrustworthy",
        }
    if observability is not None:
        output_checks = observability.get("output_checks", {})
        if isinstance(output_checks, Mapping):
            for stage_name in stages:
                check = output_checks.get(stage_name)
                if not isinstance(check, Mapping):
                    continue
                required = list(check.get("required") or [])
                output_targets.append(
                    {
                        "stage": stage_name,
                        "artifact_outputs": required,
                        "available_outputs": list(check.get("available") or []),
                        "missing_outputs": list(check.get("missing") or []),
                        "evidence_path": f"stage_observability.{mode}.{stage_group}.{stage_name}.outputs",
                    }
                )

    if mismatch.get("reason") in {"missing_field", "shape_mismatch"}:
        return {
            "kind": "field_export_precondition",
            "single_next_target": f"{mode}.{mismatch.get('field')} field export/state mapping",
            "finest_granularity": "field availability/shape",
            "evidence_path": attribution.get("field_evidence_path"),
            "why": "the required comparison field is absent or has incompatible topology, so solver-stage attribution is not yet trustworthy",
        }

    artifact_family = "solver" if stage_group == "linear_solve_results" else "matrix_operator" if stage_group in {"momentum_assembly", "pressure_assembly"} else None
    artifact_comparison = report.get("artifact_comparisons", {}).get(artifact_family) if artifact_family else None
    if isinstance(artifact_comparison, Mapping) and artifact_comparison.get("allclose") is False:
        failed_checks = artifact_comparison.get("failed_checks", [])
        first_failed = failed_checks[0] if failed_checks and isinstance(failed_checks[0], Mapping) else {}
        return {
            "kind": "failed_artifact_comparison",
            "single_next_target": f"artifact_comparisons.{artifact_family}.{first_failed.get('name')}",
            "finest_granularity": "reference-vs-candidate numeric artifact comparison",
            "stage": first_failed.get("name"),
            "artifact_family": artifact_family,
            "evidence_path": f"artifact_comparisons.{artifact_family}",
            "failed_check": first_failed,
            "why": "a direct reference-vs-candidate artifact comparison now exists and failed; fix this artifact before downstream field mismatches",
        }

    if output_targets:
        first_stage = output_targets[0]["stage"]
        return {
            "kind": "missing_artifact_comparison",
            "single_next_target": f"{mode}.{stage_group}.{first_stage} artifact comparison evidence",
            "finest_granularity": "stage output family",
            "stage": first_stage,
            "artifact_family": stage_group,
            "artifact_outputs": output_targets[0]["artifact_outputs"],
            "evidence_path": output_targets[0]["evidence_path"],
            "available_stage_outputs": output_targets,
            "why": "final field mismatch is the first numerical failure; no finer reference-vs-candidate artifact comparison is recorded yet, so the missing artifact comparison is the blocking evidence gap",
        }

    return {
        "kind": "missing_stage_granularity",
        "single_next_target": f"{mode}.{stage_group} observability/evidence",
        "finest_granularity": "stage group",
        "stage_names": stages,
        "evidence_path": attribution.get("evidence_path"),
        "why": "the harness cannot name a finer artifact from the current evidence; add or inspect finer stage artifact evidence before chasing downstream fields",
    }



def build_failure_diagnostics(mismatches: list[dict[str, Any]], report: Mapping[str, Any]) -> dict[str, Any]:
    if not mismatches:
        return {
            "status": "passed",
            "first_blocking_point": None,
            "downstream_unreliable": [],
            "basis": "all required field and mesh comparisons passed",
        }

    ordered = sorted(mismatches, key=mismatch_sort_key)
    first = ordered[0]
    attribution = first.get("attribution") if isinstance(first.get("attribution"), Mapping) else {}
    first_stage_group = attribution.get("likely_stage_group", "unknown")
    first_mode = first.get("mode")
    first_stage_index = STAGE_GROUP_ORDER.get(str(first_stage_group), len(STAGE_GROUP_ORDER))
    downstream = []
    for mismatch in ordered[1:]:
        mismatch_attr = mismatch.get("attribution") if isinstance(mismatch.get("attribution"), Mapping) else {}
        mismatch_stage = mismatch_attr.get("likely_stage_group", "unknown")
        if mismatch.get("mode") != first_mode or STAGE_GROUP_ORDER.get(str(mismatch_stage), len(STAGE_GROUP_ORDER)) >= first_stage_index:
            downstream.append(
                {
                    "mode": mismatch.get("mode"),
                    "field": mismatch.get("field"),
                    "reason": mismatch.get("reason"),
                    "stage_group": mismatch_stage,
                    "why_unreliable": "blocked by the first failing comparison; fix earlier evidence before interpreting this mismatch",
                }
            )

    observability_path = attribution.get("evidence_path")
    observed_group = None
    if observability_path:
        cursor: Any = report
        for part in str(observability_path).split("."):
            cursor = cursor.get(part) if isinstance(cursor, Mapping) else None
        if isinstance(cursor, Mapping):
            observed_group = cursor

    first_target = finest_supported_target(first, attribution, report)
    return {
        "status": "failed",
        "basis": "ordered by execution mode and solver-stage dependency before field declaration order; the first target is the finest granularity supported by recorded evidence",
        "first_blocking_point": {
            "mode": first_mode,
            "category": COMPARISON_FAILURE,
            "affected_artifact_family": "mesh_identity" if first.get("kind") == "mesh_identity" else first_target.get("artifact_family", "field_comparison"),
            "stage_group": first_stage_group,
            "stage_names": attribution.get("likely_stages", []),
            "evidence_path": first_target.get("evidence_path") or observability_path,
            "first_target": first_target,
            "field_context": field_failure_context(first),
            "stage_observability": observed_group,
        },
        "downstream_unreliable": downstream,
        "all_mismatch_count": len(mismatches),
    }


def run_openfoam_command(cmd: list[str], *, log_path: Path, timeout: int) -> dict[str, Any]:
    started = time.monotonic()
    env = openfoam_env()
    log_path.parent.mkdir(parents=True, exist_ok=True)
    try:
        with log_path.open("w") as log:
            subprocess.run(cmd, cwd=ROOT, env=env, check=True, stdout=log, stderr=subprocess.STDOUT, timeout=timeout)
    except subprocess.TimeoutExpired as exc:
        raise HarnessError(
            ORACLE_FAILURE,
            Path(cmd[0]).name,
            f"OpenFOAM command timed out after {timeout}s: {' '.join(cmd)}",
            details={"cmd": cmd, "timeout_seconds": timeout, "log_path": log_path},
            cause=exc,
        ) from exc
    except subprocess.CalledProcessError as exc:
        raise HarnessError(
            ORACLE_FAILURE,
            Path(cmd[0]).name,
            f"OpenFOAM command failed with exit code {exc.returncode}: {' '.join(cmd)}",
            details={"cmd": cmd, "returncode": exc.returncode, "log_path": log_path},
            cause=exc,
        ) from exc
    return {
        "cmd": cmd,
        "returncode": 0,
        "log_path": log_path,
        "duration_seconds": round(time.monotonic() - started, 6),
        "timeout_seconds": timeout,
    }

def diagnostic_artifact_root(args: argparse.Namespace) -> Path:
    return args.diagnostic_artifacts if args.diagnostic_artifacts is not None else args.work / "diagnostic_artifacts"


def run_reference_split_diagnostic(args: argparse.Namespace) -> tuple[dict[str, Any], dict[str, Any] | None]:
    """Run the CPU split reference in a separate verifier process for GPU artifact diagnostics."""
    work = args.work / "reference_split_diagnostic"
    report_path = work / "report.json"
    reference_artifact_root = diagnostic_artifact_root(args) / "reference_split"
    reference_artifact_path = reference_artifact_root / "split.npz"
    cmd = [
        sys.executable,
        str(Path(__file__).resolve()),
        "--backend",
        "cpu",
        "--source",
        str(args.source),
        "--work",
        str(work),
        "--report",
        str(report_path),
        "--rtol",
        str(args.rtol),
        "--atol",
        str(args.atol),
        "--diagnostic-artifacts",
        str(reference_artifact_root),
    ]
    started = time.monotonic()
    env = dict(os.environ)
    env.pop("PYTHONPATH", None)
    completed = subprocess.run(cmd, cwd=ROOT, env=env, stdout=subprocess.PIPE, stderr=subprocess.STDOUT, text=True, timeout=900)
    summary = {
        "cmd": cmd,
        "exit_code": completed.returncode,
        "duration_seconds": round(time.monotonic() - started, 6),
        "report": report_path,
        "diagnostic_artifacts": reference_artifact_path,
        "stdout_tail": completed.stdout[-4000:],
    }
    if not report_path.exists():
        return summary, None
    try:
        return summary, json.loads(report_path.read_text())
    except Exception as exc:
        summary["report_parse_error"] = {"type": type(exc).__name__, "message": str(exc)}
        return summary, None


def patch_start_from_latest(case: Path) -> None:
    path = case / "system/controlDict"
    text = path.read_text()
    old = "startFrom       startTime;"
    new = "startFrom       latestTime;"
    if old not in text:
        raise AssertionError(f"{path}: missing {old!r}")
    path.write_text(text.replace(old, new, 1))


def import_foam() -> Any:
    apply_openfoam_env()
    import foam_stepper as foam

    return foam


def clone_prepared_case(prepared: Path, dest: Path, source: Path, metadata: Any) -> dict[str, Any]:
    if dest.exists():
        shutil.rmtree(dest)
    shutil.copytree(prepared, dest, copy_function=shutil.copy2)
    return write_case_manifest(source, dest, metadata)


def prepared_case_digest(manifest: Mapping[str, Any]) -> str:
    prepared = manifest.get("prepared_case", {})
    if not isinstance(prepared, Mapping):
        return ""
    return str(prepared.get("files_sha256", ""))


def prepared_file_count(manifest: Mapping[str, Any]) -> int:
    prepared = manifest.get("prepared_case", {})
    files = prepared.get("files", []) if isinstance(prepared, Mapping) else []
    return len(files) if isinstance(files, list) else 0


def compare_prepared_case_manifests(manifests: Mapping[str, Mapping[str, Any]]) -> dict[str, dict[str, Any]]:
    baseline = prepared_case_digest(manifests["prepared"])
    return {
        role: {
            "matches_prepared": prepared_case_digest(manifest) == baseline,
            "files_sha256": prepared_case_digest(manifest),
            "file_count": prepared_file_count(manifest),
        }
        for role, manifest in manifests.items()
    }


def prepare_work_cases(source: Path, work: Path, *, include_split: bool) -> dict[str, Any]:
    if work.exists():
        shutil.rmtree(work)
    work.mkdir(parents=True)
    prepared_case = work / "prepared_case"
    oracle = work / "oracle_case"
    run_one = work / "run_one_case"
    split = work / "split_case" if include_split else None
    reference_split = work / "reference_split_case" if include_split else None

    prepared_meta = prepare_case(source, prepared_case, end_time=1)
    manifests: dict[str, Any] = {"prepared": load_case_manifest(prepared_case)}
    manifests["oracle"] = clone_prepared_case(prepared_case, oracle, source, prepared_meta)
    manifests["run_one"] = clone_prepared_case(prepared_case, run_one, source, prepared_meta)
    if split is not None:
        manifests["split"] = clone_prepared_case(prepared_case, split, source, prepared_meta)
        manifests["reference_split"] = clone_prepared_case(prepared_case, reference_split, source, prepared_meta)

    return {
        "prepared_case": prepared_case,
        "oracle_case": oracle,
        "run_one_case": run_one,
        "split_case": split,
        "reference_split_case": reference_split,
        "metadata": {
            "prepared": prepared_meta,
            "oracle": prepared_meta,
            "run_one": prepared_meta,
            "split": prepared_meta if split is not None else None,
            "reference_split": prepared_meta if reference_split is not None else None,
        },
        "manifests": manifests,
        "prepared_case_identity": compare_prepared_case_manifests(manifests),
    }


def _gpu_backend_module() -> Any:
    try:
        from foam_stepper.gpu import backend as gpu_backend
    except Exception as exc:
        raise BackendExecutionError(
            "select_backend",
            "GPU runtime unavailable: failed to import reusable foam_stepper.gpu backend",
            details={
                "requested": "gpu",
                "failure_kind": GPU_RUNTIME_UNAVAILABLE,
                "provider": "quadrants_cuda_rans_solver",
                "used_cpu_fallback": False,
                "cause": {"type": type(exc).__name__, "message": str(exc)},
            },
        ) from exc
    return gpu_backend


def _wrap_gpu_backend_error(exc: BaseException) -> BackendExecutionError:
    return BackendExecutionError(
        str(getattr(exc, "step", "execute_gpu_solver_contract")),
        str(exc),
        details=dict(getattr(exc, "details", {}) or {}),
    )


def _call_gpu_backend(function_name: str, *args: Any, **kwargs: Any) -> Any:
    gpu_backend = _gpu_backend_module()
    try:
        return getattr(gpu_backend, function_name)(*args, **kwargs)
    except gpu_backend.BackendExecutionError as exc:
        raise _wrap_gpu_backend_error(exc) from exc


def prepare_gpu_solver_inputs(foam: Any, stepper: Any, case: Path, prepared: Mapping[str, Any], backend: Mapping[str, Any]) -> dict[str, Any]:
    return _call_gpu_backend("prepare_gpu_solver_inputs", foam, stepper, case, prepared, backend)


def run_gpu_solver_stage_smoke(stepper: Any, backend: Mapping[str, Any], case: Path, *, diagnostic_artifact_path: Path | None = None, trace_substitutions: Mapping[str, np.ndarray] | None = None) -> dict[str, Any]:
    return _call_gpu_backend("run_gpu_solver_stage_smoke", stepper, backend, case, diagnostic_artifact_path=diagnostic_artifact_path, trace_substitutions=trace_substitutions)


def gpu_solver_stage_report(gpu_run: Mapping[str, Any]) -> dict[str, Any]:
    return _call_gpu_backend("gpu_solver_stage_report", gpu_run)


def gpu_solver_input_blocker_summary(gpu_inputs: Mapping[str, Any]) -> dict[str, Any]:
    return _call_gpu_backend("gpu_solver_input_blocker_summary", gpu_inputs)


def gpu_solver_contract_blocker(backend: Mapping[str, Any], case: Path) -> BackendExecutionError | None:
    gpu_backend = _gpu_backend_module()
    try:
        blocker = gpu_backend.gpu_solver_contract_blocker(backend, case)
    except gpu_backend.BackendExecutionError as exc:
        raise _wrap_gpu_backend_error(exc) from exc
    if blocker is None:
        return None
    return _wrap_gpu_backend_error(blocker)


def record_gpu_contract_failure(report: dict[str, Any], blocker: BackendExecutionError, gpu_inputs: Mapping[str, Any] | None = None) -> None:
    if gpu_inputs is not None:
        blocker.details["gpu_solver_inputs"] = gpu_solver_input_blocker_summary(gpu_inputs)
    report["backend"]["contract_preflight"] = json_ready(
        {
            "status": "failed",
            "step": blocker.step,
            "failure_kind": blocker.details.get("failure_kind"),
            "checked_before_oracle": True,
            "gpu_inputs_prepared": gpu_inputs is not None,
            "reason": "GPU requests must fail as incomplete before any CPU/OpenFOAM solver path can be substituted as the backend.",
        }
    )

    raise blocker.to_harness_error() from blocker


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

    if requested == "gpu":
        return _call_gpu_backend("select_execution_backend", requested)

    return {
        "requested": requested,
        "selected": "cpu",
        "device": "host",
        "provider": "foam_stepper_cpu",
        "available_backends": ["cpu"],
        "gpu_device_present": any(Path(path).exists() for path in ("/dev/nvidia0", "/dev/dri/renderD128")),
        "capabilities": [
            "foam_stepper_python_bridge",
            "openfoam_case_loader",
            "oracle_comparison",
        ],
        "used_cpu_fallback": False,
    }


def run_backend_iteration(stepper: Any, backend: Mapping[str, Any], case: Path) -> dict[str, Any]:
    selected = backend.get("selected")
    if selected == "gpu":
        return _call_gpu_backend("run_backend_iteration", stepper, backend, case)
    if selected != "cpu":
        raise BackendExecutionError(
            "execute_backend",
            f"backend {selected!r} is not executable",
            details={"backend": backend, "case": case},
        )
    result = stepper.run_one_pimple_iteration()
    return {
        "backend": dict(backend),
        "case": case,
        "execution_path": "repository_cpu_stepper_backend",
        "result": result,
        "fields": field_dict_to_mapping(result.outputs["fields"]),
    }


def run_gpu_split_iteration(foam: Any, stepper: Any, backend: Mapping[str, Any], case: Path, *, diagnostic_artifact_path: Path | None = None, trace_substitutions: Mapping[str, np.ndarray] | None = None) -> dict[str, Any]:
    return _call_gpu_backend("run_gpu_split_iteration", foam, stepper, backend, case, diagnostic_artifact_path=diagnostic_artifact_path, trace_substitutions=trace_substitutions)


def make_stepper(foam: Any, case: Path, label: str) -> Any:
    try:
        return foam.Case(case).make_stepper()
    except Exception as exc:
        raise HarnessError(
            LOADING_FAILURE,
            f"load_{label}_case",
            f"failed to load {label} case with foam_stepper: {case}",
            details={"case": case},
            cause=exc,
        ) from exc


def read_mesh_identity(stepper: Any, label: str) -> dict[str, Any]:
    try:
        return mesh_identity(stepper.mesh())
    except Exception as exc:
        raise HarnessError(
            LOADING_FAILURE,
            f"read_{label}_mesh",
            f"failed to read {label} mesh identity",
            details={"case": getattr(stepper, "case_path", "")},
            cause=exc,
        ) from exc


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

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


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


def export_matrix_state_checked(foam: Any, matrix: Any, solver_state: Mapping[str, Any], label: str) -> dict[str, Any]:
    try:
        return foam.export_matrix_state(matrix, mesh_state=solver_state)
    except Exception as exc:
        raise HarnessError(
            STEPPER_FAILURE,
            f"export_{label}_matrix_state",
            f"failed to export explicit {label} matrix state",
            details=validation_details(exc),
            cause=exc,
        ) from exc


def checked_step(stages: list[dict[str, Any]], name: str, fn: Any) -> Any:
    try:
        result = fn()
    except Exception as exc:
        raise HarnessError(
            STEPPER_FAILURE,
            name,
            f"split-step stage failed: {name}",
            cause=exc,
        ) from exc
    stages.append(transform_summary(result))
    return result


def run_split_iteration(foam: Any, stepper: Any, *, diagnostic_artifact_path: Path | None = None) -> dict[str, Any]:
    stages: list[dict[str, Any]] = []
    checked_step(stages, "pre_solve", stepper.pre_solve)
    checked_step(stages, "advance_time", stepper.advance_time)
    begin = checked_step(stages, "begin_pimple_iteration", stepper.begin_pimple_iteration)
    if begin.outputs.get("active") is not True:
        raise HarnessError(
            STEPPER_FAILURE,
            "begin_pimple_iteration",
            "split-step PIMPLE iteration was inactive",
            details={"outputs": summarize_value(begin.outputs)},
        )

    checked_step(stages, "fv_models_correct", stepper.fv_models_correct)
    checked_step(stages, "pre_predictor", stepper.pre_predictor)
    checked_step(stages, "momentum_transport_predict", stepper.momentum_transport_predictor)
    terms = checked_step(stages, "assemble_momentum_terms", stepper.assemble_momentum_terms)
    UEqn = checked_step(stages, "assemble_UEqn", stepper.assemble_momentum_matrix)
    relax_UEqn = checked_step(stages, "relax_UEqn", stepper.relax_matrix)
    checked_step(stages, "constrain_UEqn", stepper.constrain_matrix)
    solve_UEqn = checked_step(stages, "solve_UEqn", stepper.solve_momentum)
    compute_pressure_inputs = checked_step(stages, "compute_pressure_inputs", stepper.compute_pressure_inputs)
    pEqn = checked_step(stages, "assemble_pEqn", stepper.assemble_pressure_matrix)
    solve_pEqn = checked_step(stages, "solve_pEqn", stepper.solve_pressure)
    correct_velocity_pressure_flux = checked_step(stages, "correct_velocity_pressure_flux", stepper.correct_velocity_pressure_flux)
    momentum_transport_correct = checked_step(stages, "momentum_transport_correct", stepper.momentum_transport_corrector)
    checked_step(stages, "end_pimple_iteration", stepper.end_pimple_iteration)
    checked_step(stages, "post_solve", lambda: stepper.post_solve(write=False))

    try:
        fields = field_dict_to_mapping(stepper.fields())
    except Exception as exc:
        raise HarnessError(
            STEPPER_FAILURE,
            "split_fields",
            "failed to read split-step fields after execution",
            cause=exc,
        ) from exc

    solver_state = export_solver_state_checked(foam, stepper, "split", fields)

    momentum_terms = [term.get("name", "") for term in terms.outputs.get("terms", [])]
    UEqn_matrix = UEqn.outputs["UEqn"]
    pEqn_matrix = pEqn.outputs["pEqn"]
    matrix_states = {
        "UEqn": foam.describe_matrix_state(export_matrix_state_checked(foam, UEqn_matrix, solver_state, "UEqn")),
        "pEqn": foam.describe_matrix_state(export_matrix_state_checked(foam, pEqn_matrix, solver_state, "pEqn")),
    }
    diagnostic_artifacts = write_split_numeric_artifacts(
        diagnostic_artifact_path,
        role="split",
        momentum_terms=terms,
        assemble_UEqn=UEqn,
        relax_UEqn=relax_UEqn,
        solve_UEqn=solve_UEqn,
        compute_pressure_inputs=compute_pressure_inputs,
        assemble_pEqn=pEqn,
        solve_pEqn=solve_pEqn,
        correct_velocity_pressure_flux=correct_velocity_pressure_flux,
        momentum_transport_correct=momentum_transport_correct,
        fields=fields,
    )
    observability = stage_observability_report(
        "split",
        stages,
        evidence={
            "split_step_execution": True,
            "momentum_terms": momentum_terms,
            "matrix_states": sorted(matrix_states),
            "turbulence_fields": visible_turbulence_fields(fields),
        },
    )
    return {
        "fields": fields,
        "state": solver_state,
        "state_summary": foam.describe_solver_state(solver_state),
        "stages": stages,
        "graph": [stage["name"] for stage in stages],
        "momentum_terms": momentum_terms,
        "UEqn": matrix_summary(UEqn_matrix),
        "pEqn": matrix_summary(pEqn_matrix),
        "matrix_states": matrix_states,
        "observability": observability,
        "diagnostic_artifacts": diagnostic_artifacts,
    }


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


def differentiability_report(backend: Mapping[str, Any], *, case_name: str) -> dict[str, Any]:
    capabilities = tuple(backend.get("capabilities", ()))
    autodiff_available = any(capability in capabilities for capability in ("autodiff", "jvp", "vjp", "gradient"))
    status = "differentiable" if autodiff_available else "not_differentiable"
    return {
        "status": status,
        "backend": {
            "requested": backend.get("requested"),
            "selected": backend.get("selected"),
            "provider": backend.get("provider"),
            "capabilities": list(capabilities),
        },
        "variables": DIFFERENTIABILITY_VARIABLES,
        "claims": [
            {
                "name": "current_migrated_solver_path",
                "differentiable": autodiff_available,
                "variables": DIFFERENTIABILITY_VARIABLES,
                "applies_to": {
                    "modes": ["run_one", "split"],
                    "fields": list(REQUIRED_FIELDS),
                    "parameters": list(DIFFERENTIABILITY_VARIABLES["case_parameters"]),
                    "losses": list(DIFFERENTIABILITY_VARIABLES["losses"]),
                },
                "evidence": {
                    "backend_capabilities": list(capabilities),
                    "reason": "no registered backend capability exposes AD/JVP/VJP for solver updates"
                    if not autodiff_available
                    else "backend advertises differentiability capabilities",
                },
            }
        ],
        "non_differentiable": list(NON_DIFFERENTIABLE_BOUNDARIES),
        "custom_gradient_required": list(CUSTOM_GRADIENT_REQUIREMENTS),
        "checks": [
            {
                "name": "solver_sensitivity_against_finite_difference",
                "status": "not_run",
                "fixed_problem_setup": {
                    "case": case_name,
                    "fields": list(REQUIRED_FIELDS),
                    "mesh_topology": "fixed prepared AirfRANS mesh",
                    "backend": backend.get("selected"),
                },
                "computed_sensitivity": None,
                "independent_numerical_check": None,
                "reason": "no differentiable backend or custom gradient is currently available, so no computed sensitivity is claimed",
            }
        ],
    }


def command_timing_summary(command: Mapping[str, Any]) -> dict[str, Any]:
    cmd = list(command.get("cmd") or ())
    executable = Path(str(cmd[0])).name if cmd else None
    log_path = str(command.get("log_path", ""))
    role = "openfoam_command"
    if executable == "checkMesh":
        role = "oracle_mesh_check"
    elif executable == "foamRun" and "oracle" in log_path:
        role = "oracle_solver"
    return {
        "role": role,
        "executable": executable,
        "cmd": cmd,
        "duration_seconds": command.get("duration_seconds"),
        "timeout_seconds": command.get("timeout_seconds"),
        "log_path": command.get("log_path"),
        "returncode": command.get("returncode"),
    }


def gpu_timing_summary(source: Mapping[str, Any], *, label: str) -> dict[str, Any]:
    profiler = source.get("profiler", {}) if isinstance(source.get("profiler"), Mapping) else {}
    timing = source.get("timing", {}) if isinstance(source.get("timing"), Mapping) else {}
    return {
        "label": label,
        "execution_path": source.get("execution_path"),
        "timing": timing,
        "profiler": {
            "available": profiler.get("available"),
            "profile_record_count": profiler.get("profile_record_count"),
            "device_time_ms_total": profiler.get("device_time_ms_total"),
            "expected_kernel_entrypoints": profiler.get("expected_kernel_entrypoints"),
            "generated_kernel_names": profiler.get("generated_kernel_names"),
        },
    }


def build_timing_evidence(report: Mapping[str, Any]) -> dict[str, Any]:
    commands = [command_timing_summary(command) for command in report.get("commands", [])]
    oracle_solver = next((command for command in commands if command.get("role") == "oracle_solver"), None)
    gpu_sources: dict[str, Any] = {}
    gpu_stage_smoke = report.get("gpu_solver_stages")
    if isinstance(gpu_stage_smoke, Mapping) and gpu_stage_smoke.get("status") == "executed":
        gpu_sources["stage_preflight"] = gpu_timing_summary(gpu_stage_smoke, label="stage_preflight")
    for mode in ("run_one", "split"):
        mode_gpu = report.get("modes", {}).get(mode, {}).get("gpu_solver", {})
        if isinstance(mode_gpu, Mapping) and mode_gpu.get("status") == "executed":
            gpu_sources[mode] = gpu_timing_summary(mode_gpu, label=mode)

    gpu_solver_timing_ok = any(
        isinstance(source.get("timing"), Mapping) and source["timing"].get("gpu_solver_wall_seconds") is not None
        for source in gpu_sources.values()
    )
    gpu_profiler_ok = any(
        isinstance(source.get("profiler"), Mapping) and source["profiler"].get("available") is True and source["profiler"].get("profile_record_count", 0) > 0
        for source in gpu_sources.values()
    )
    openfoam_timing_ok = oracle_solver is not None and oracle_solver.get("duration_seconds") is not None
    return {
        "schema_version": 1,
        "status": "reported" if commands or gpu_sources else "not_reported",
        "openfoam_reference": {
            "oracle_solver_duration_seconds": None if oracle_solver is None else oracle_solver.get("duration_seconds"),
            "commands": commands,
        },
        "gpu_solver": gpu_sources,
        "criteria": {
            "openfoam_reference_timing": openfoam_timing_ok,
            "gpu_solver_wall_timing": gpu_solver_timing_ok,
            "gpu_profiler_timing": gpu_profiler_ok,
        },
    }

def path_presence(report: Mapping[str, Any], paths: Iterable[str]) -> tuple[list[str], list[str]]:
    available = []
    missing = []
    for path in paths:
        (available if report_path_exists(report, path) else missing).append(path)
    return available, missing

def artifact_stats(value: Any) -> Mapping[str, Any] | None:
    if not isinstance(value, Mapping):
        return None
    if "sha256" in value and "shape" in value:
        return value
    stats = value.get("stats")
    if isinstance(stats, Mapping) and "sha256" in stats and "shape" in stats:
        return stats
    internal = value.get("internal")
    if isinstance(internal, Mapping) and "sha256" in internal and "shape" in internal:
        return internal
    return None


def compare_artifact_stats(report: Mapping[str, Any], *, name: str, reference_path: str, candidate_path: str) -> dict[str, Any]:
    reference = artifact_stats(report_path_value(report, reference_path))
    candidate = artifact_stats(report_path_value(report, candidate_path))
    out: dict[str, Any] = {
        "name": name,
        "reference_path": reference_path,
        "candidate_path": candidate_path,
        "reference_available": reference is not None,
        "candidate_available": candidate is not None,
    }
    if reference is None or candidate is None:
        return {**out, "shape_matches": False, "allclose": False, "reason": "missing_artifact_stats"}
    shape_matches = reference.get("shape") == candidate.get("shape")
    dtype_matches = reference.get("dtype") == candidate.get("dtype")
    sha_matches = reference.get("sha256") == candidate.get("sha256")
    return {
        **out,
        "shape_matches": shape_matches,
        "dtype_matches": dtype_matches,
        "allclose": bool(shape_matches and dtype_matches and sha_matches),
        "reason": None if shape_matches and dtype_matches and sha_matches else "artifact_hash_mismatch",
        "reference": {key: reference.get(key) for key in ("shape", "dtype", "size", "sha256", "min", "max", "mean")},
        "candidate": {key: candidate.get(key) for key in ("shape", "dtype", "size", "sha256", "min", "max", "mean")},
    }


def summarize_artifact_comparison(name: str, checks: list[dict[str, Any]]) -> dict[str, Any]:
    failures = [check for check in checks if check.get("allclose") is not True]
    return {
        "name": name,
        "shape_matches": all(check.get("shape_matches") is True for check in checks),
        "allclose": not failures,
        "reason": None if not failures else failures[0].get("reason"),
        "checks": checks,
        "failed_checks": failures,
    }


def build_artifact_comparisons(report: Mapping[str, Any], *, mesh_context: Any | None = None) -> dict[str, Any]:
    pressure_input_pairs = [
        ("HbyA", "pressure_inputs.HbyA", "modes.reference_split.stages.compute_pressure_inputs.outputs.HbyA.internal", "modes.split.stages.compute_pressure_inputs.outputs.HbyA"),
        ("phiHbyA", "pressure_inputs.phiHbyA", "modes.reference_split.stages.compute_pressure_inputs.outputs.phiHbyA.internal", "modes.split.stages.compute_pressure_inputs.outputs.phiHbyA"),
        ("rAU", "pressure_inputs.rAU", "modes.reference_split.stages.compute_pressure_inputs.outputs.rAU.internal", "modes.split.stages.compute_pressure_inputs.outputs.rAU"),
        ("rAtU", "pressure_inputs.rAtU", "modes.reference_split.stages.compute_pressure_inputs.outputs.rAtU.internal", "modes.split.stages.compute_pressure_inputs.outputs.rAtU"),
    ]
    matrix_pairs = [
        ("UEqn.diag", "matrix_operator.UEqn.diag", "modes.reference_split.stages.relax_UEqn.outputs.UEqn.diag", "modes.split.UEqn.diag"),
        ("UEqn.upper", "matrix_operator.UEqn.upper", "modes.reference_split.stages.relax_UEqn.outputs.UEqn.upper", "modes.split.UEqn.upper"),
        ("UEqn.lower", "matrix_operator.UEqn.lower", "modes.reference_split.stages.relax_UEqn.outputs.UEqn.lower", "modes.split.UEqn.lower"),
        ("UEqn.source", "matrix_operator.UEqn.source", "modes.reference_split.stages.relax_UEqn.outputs.UEqn.source", "modes.split.UEqn.source"),
        ("UEqn.psi", "matrix_operator.UEqn.psi", "modes.reference_split.stages.relax_UEqn.outputs.UEqn.psi", "modes.split.UEqn.psi"),
        ("pEqn.diag", "matrix_operator.pEqn.diag", "modes.reference_split.pEqn.diag", "modes.split.pEqn.diag"),
        ("pEqn.upper", "matrix_operator.pEqn.upper", "modes.reference_split.pEqn.upper", "modes.split.pEqn.upper"),
        ("pEqn.source", "matrix_operator.pEqn.source", "modes.reference_split.pEqn.source", "modes.split.pEqn.source"),
        ("pEqn.psi", "matrix_operator.pEqn.psi", "modes.reference_split.pEqn.psi", "modes.split.pEqn.psi"),
    ]
    solver_pairs = [
        ("solve_UEqn.field_after", "solver.solve_UEqn.field_after", "modes.reference_split.stages.solve_UEqn.outputs.field_after", "modes.split.stages.solve_UEqn.outputs.field_after"),
        ("solve_pEqn.p", "solver.solve_pEqn.p", "modes.reference_split.stages.solve_pEqn.outputs.p", "modes.split.stages.solve_pEqn.outputs.p"),
        ("solve_pEqn.phi", "solver.solve_pEqn.phi", "modes.reference_split.stages.solve_pEqn.outputs.phi", "modes.split.stages.solve_pEqn.outputs.phi"),
    ]

    reference_artifact_path = diagnostic_artifact_path(report, "reference_split")
    candidate_artifact_path = diagnostic_artifact_path(report, "split")
    reference_artifacts = load_numeric_artifact_file(reference_artifact_path)
    candidate_artifacts = load_numeric_artifact_file(candidate_artifact_path)
    tolerances = report.get("tolerances", {}) if isinstance(report.get("tolerances"), Mapping) else {}
    rtol = float(tolerances.get("rtol", 0.0) or 0.0)
    atol = float(tolerances.get("atol", 0.0) or 0.0)
    backend_selected = report.get("backend", {}).get("selected") if isinstance(report.get("backend"), Mapping) else None


    def compare_pair(name: str, key: str, reference_path: str, candidate_path: str) -> dict[str, Any]:
        if backend_selected == "cpu" and reference_artifacts is None:
            if candidate_artifacts is not None:
                return compare_numeric_artifact_arrays(
                    name,
                    key,
                    candidate_artifacts.get(key),
                    candidate_artifacts.get(key),
                    rtol=rtol,
                    atol=atol,
                    mesh_context=mesh_context,
                )
            return compare_artifact_stats(report, name=name, reference_path=candidate_path, candidate_path=candidate_path)
        if reference_artifacts is not None and candidate_artifacts is not None:
            return compare_numeric_artifact_arrays(
                name,
                key,
                reference_artifacts.get(key),
                candidate_artifacts.get(key),
                rtol=rtol,
                atol=atol,
                mesh_context=mesh_context,
            )
        return compare_artifact_stats(report, name=name, reference_path=reference_path, candidate_path=candidate_path)

    pressure_input_checks = [compare_pair(name, key, reference, candidate) for name, key, reference, candidate in pressure_input_pairs]
    matrix_checks = [compare_pair(name, key, reference, candidate) for name, key, reference, candidate in matrix_pairs]
    solver_checks = [compare_pair(name, key, reference, candidate) for name, key, reference, candidate in solver_pairs]
    return {
        "schema_version": 2,
        "numeric_artifacts": {
            "reference_path": str(reference_artifact_path) if reference_artifact_path is not None else None,
            "candidate_path": str(candidate_artifact_path) if candidate_artifact_path is not None else None,
            "reference_loaded": reference_artifacts is not None,
            "candidate_loaded": candidate_artifacts is not None,
            "comparison_basis": "npz_numeric_artifacts" if reference_artifacts is not None and candidate_artifacts is not None else "candidate_numeric_self_check" if backend_selected == "cpu" and candidate_artifacts is not None else "summary_hash_stats",
        },
        "pressure_inputs": summarize_artifact_comparison("pressure_inputs", pressure_input_checks),
        "matrix_operator": summarize_artifact_comparison("matrix_operator", matrix_checks),
        "solver": summarize_artifact_comparison("solver", solver_checks),
    }


def comparison_paths_failed(report: Mapping[str, Any], paths: Iterable[str]) -> list[dict[str, Any]]:
    failed = []
    for path in paths:
        value = report_path_value(report, path)
        if isinstance(value, Mapping):
            if value.get("matches_oracle") is False or value.get("allclose") is False or value.get("shape_matches") is False:
                failed_checks = value.get("failed_checks", [])
                first_failed = failed_checks[0] if failed_checks and isinstance(failed_checks[0], Mapping) else value
                failed.append(
                    {
                        "path": path,
                        "status": "failed",
                        "reason": first_failed.get("reason") or value.get("reason"),
                        "first_failed_check": first_failed.get("name"),
                        "max_abs": first_failed.get("max_abs"),
                        "mean_abs": first_failed.get("mean_abs"),
                        "rms_abs": first_failed.get("rms_abs"),
                    }
                )
    return failed


def family_downstream_of_first(family: Mapping[str, Any], first_blocker: Mapping[str, Any]) -> bool:
    first_group = first_blocker.get("stage_group")
    if not first_group or first_group == "unknown":
        return False
    first_index = STAGE_GROUP_ORDER.get(str(first_group))
    family_groups = family.get("stage_groups", ())
    family_indexes = [STAGE_GROUP_ORDER[group] for group in family_groups if group in STAGE_GROUP_ORDER]
    return first_index is not None and bool(family_indexes) and min(family_indexes) > first_index


def intermediate_artifact_family_report(report: Mapping[str, Any]) -> dict[str, Any]:
    diagnostics = report.get("failure_diagnostics", {}) if isinstance(report.get("failure_diagnostics"), Mapping) else {}
    backend_selected = report.get("backend", {}).get("selected")
    first_blocker = diagnostics.get("first_blocking_point", {}) if isinstance(diagnostics.get("first_blocking_point"), Mapping) else {}
    families = []
    for family in INTERMEDIATE_ARTIFACT_FAMILIES:
        reference_available, reference_missing = path_presence(report, family["reference_paths"])
        candidate_available, candidate_missing = path_presence(report, family["candidate_paths"])
        comparison_available, comparison_missing = path_presence(report, family["comparison_paths"])
        comparison_failures = comparison_paths_failed(report, comparison_available)
        contract_gaps = []
        if backend_selected != "cpu":
            if not family["reference_paths"]:
                contract_gaps.append("reference_artifact_evidence")
            if not family["comparison_paths"]:
                contract_gaps.append("direct_reference_candidate_comparison")
        downstream = family_downstream_of_first(family, first_blocker)
        if downstream:
            status = "downstream_unreliable"
        elif comparison_failures:
            status = "failed"
        elif backend_selected == "cpu" and not candidate_missing and candidate_available and not comparison_failures:
            status = "passed"
        elif reference_missing or candidate_missing or comparison_missing or contract_gaps:
            status = "missing"
        elif comparison_available:
            status = "passed"
        else:
            status = "available"
        families.append(
            {
                "name": family["name"],
                "required": family["required"],
                "status": status,
                "stage_groups": list(family["stage_groups"]),
                "reference_evidence": {"available": reference_available, "missing": reference_missing},
                "candidate_evidence": {"available": candidate_available, "missing": candidate_missing},
                "comparison_evidence": {"available": comparison_available, "missing": comparison_missing, "failures": comparison_failures},
                "missing_core_evidence": contract_gaps,
                "comparison_basis": "cpu_openfoam_artifact_self_check" if backend_selected == "cpu" and status == "passed" and not comparison_available else "direct_reference_candidate_comparison",
                "blocked_by": first_blocker.get("first_target") if downstream else None,
                "why": (
                    "later artifact family is downstream of the first blocker"
                    if downstream
                    else "direct reference-vs-candidate comparison failed"
                    if comparison_failures
                    else "required reference, candidate, or direct comparison evidence is missing"
                    if status == "missing"
                    else "CPU/OpenFOAM artifact evidence is present on the trusted execution path"
                    if backend_selected == "cpu" and status == "passed" and not comparison_available
                    else "required direct comparison evidence is present and passed"
                    if status == "passed"
                    else "artifact evidence is present but has no direct comparison contract"
                ),
            }
        )
    first_actionable = next((family for family in families if family["status"] in {"failed", "missing"}), None)
    status = "complete" if all(family["status"] == "passed" for family in families if family["required"]) else "incomplete"
    return {
        "schema_version": 1,
        "status": status,
        "families": families,
        "first_actionable_family": first_actionable,
        "status_legend": {
            "passed": "required reference, candidate, and direct comparison evidence exists and passed",
            "failed": "direct comparison evidence exists and failed",
            "missing": "required reference, candidate, or comparison evidence is absent",
            "downstream_unreliable": "earlier blocker makes this family unsuitable for diagnosis",
            "available": "artifacts are present but not directly compared",
        },
    }

def first_divergence_summary(report: Mapping[str, Any]) -> dict[str, Any]:
    diagnostics = report.get("failure_diagnostics", {}) if isinstance(report.get("failure_diagnostics"), Mapping) else {}
    first_blocker = diagnostics.get("first_blocking_point") if isinstance(diagnostics.get("first_blocking_point"), Mapping) else None
    artifacts = report.get("intermediate_artifacts", {}) if isinstance(report.get("intermediate_artifacts"), Mapping) else {}
    families = artifacts.get("families", []) if isinstance(artifacts.get("families"), list) else []
    family_statuses = {
        family.get("name"): family.get("status")
        for family in families
        if isinstance(family, Mapping)
    }
    downstream_ignore = [
        {
            "mode": item.get("mode"),
            "field": item.get("field"),
            "artifact_family": item.get("stage_group"),
            "reason": item.get("why_unreliable"),
        }
        for item in diagnostics.get("downstream_unreliable", [])
        if isinstance(item, Mapping)
    ]
    trace = report.get("differential_trace", {}) if isinstance(report.get("differential_trace"), Mapping) else {}
    trace_first = trace.get("first_divergence") if isinstance(trace.get("first_divergence"), Mapping) else None
    if trace_first is not None:
        metadata = trace_first.get("metadata", {}) if isinstance(trace_first.get("metadata"), Mapping) else {}
        largest = trace_first.get("largest_difference", {}) if isinstance(trace_first.get("largest_difference"), Mapping) else {}
        return {
            "schema_version": 1,
            "status": "failed",
            "category": COMPARISON_FAILURE,
            "first_target": trace_first.get("name"),
            "target_kind": "differential_trace_checkpoint",
            "evidence_path": "differential_trace.first_divergence",
            "artifact_family": metadata.get("family"),
            "mode": "split",
            "stage_group": metadata.get("family"),
            "stage_names": [metadata.get("lifecycle_phase")],
            "field": metadata.get("field") or metadata.get("field_or_matrix"),
            "reason": trace_first.get("reason"),
            "local_context": largest.get("local_context"),
            "missing_evidence": None,
            "trace_lifecycle_phase": metadata.get("lifecycle_phase"),
            "trace_max_abs": trace_first.get("max_abs"),
            "trace_mean_abs": trace_first.get("mean_abs"),
            "trace_rms_abs": trace_first.get("rms_abs"),
            "trace_substitution_supported": metadata.get("substitution_supported"),
            "intermediate_artifact_statuses": family_statuses,
            "downstream_symptoms_to_ignore": downstream_ignore,
        }
    if first_blocker is not None:
        field_context = first_blocker.get("field_context", {}) if isinstance(first_blocker.get("field_context"), Mapping) else {}
        first_target = first_blocker.get("first_target", {}) if isinstance(first_blocker.get("first_target"), Mapping) else {}
        return {
            "schema_version": 1,
            "status": "failed",
            "category": first_blocker.get("category"),
            "first_target": first_target.get("single_next_target"),
            "target_kind": first_target.get("kind"),
            "evidence_path": first_target.get("evidence_path") or first_blocker.get("evidence_path"),
            "artifact_family": first_blocker.get("affected_artifact_family") or first_target.get("artifact_family"),
            "mode": first_blocker.get("mode"),
            "stage_group": first_blocker.get("stage_group"),
            "stage_names": first_blocker.get("stage_names"),
            "field": field_context.get("field"),
            "reason": field_context.get("reason") or first_target.get("why"),
            "local_context": field_context.get("local_context"),
            "missing_evidence": first_target.get("missing_core_evidence") or first_target.get("artifact_outputs"),
            "intermediate_artifact_statuses": family_statuses,
            "downstream_symptoms_to_ignore": downstream_ignore,
        }

    first_family = artifacts.get("first_actionable_family") if isinstance(artifacts.get("first_actionable_family"), Mapping) else None
    if first_family is not None:
        missing_evidence = {
            "reference": first_family.get("reference_evidence", {}).get("missing"),
            "candidate": first_family.get("candidate_evidence", {}).get("missing"),
            "comparison": first_family.get("comparison_evidence", {}).get("missing"),
            "core": first_family.get("missing_core_evidence"),
        }
        return {
            "schema_version": 1,
            "status": first_family.get("status") or "incomplete",
            "category": REQUIRED_EVIDENCE_FAILURE,
            "first_target": f"{first_family.get('name')} artifact evidence",
            "target_kind": "intermediate_artifact_family",
            "evidence_path": "intermediate_artifacts.first_actionable_family",
            "artifact_family": first_family.get("name"),
            "mode": None,
            "stage_group": (first_family.get("stage_groups") or [None])[0],
            "stage_names": first_family.get("stage_groups"),
            "field": None,
            "reason": first_family.get("why"),
            "local_context": None,
            "missing_evidence": missing_evidence,
            "intermediate_artifact_statuses": family_statuses,
            "downstream_symptoms_to_ignore": [
                {
                    "artifact_family": family.get("name"),
                    "reason": "ignore until the first missing or failed artifact family is closed",
                }
                for family in families
                if isinstance(family, Mapping) and family.get("status") == "downstream_unreliable"
            ],
        }

    return {
        "schema_version": 1,
        "status": "passed",
        "category": None,
        "first_target": None,
        "target_kind": None,
        "evidence_path": None,
        "artifact_family": None,
        "mode": None,
        "stage_group": None,
        "stage_names": [],
        "field": None,
        "reason": "all required comparisons and artifact families passed",
        "local_context": None,
        "missing_evidence": None,
        "intermediate_artifact_statuses": family_statuses,
        "downstream_symptoms_to_ignore": [],
    }



def enabled_mode_names(report: Mapping[str, Any]) -> list[str]:
    return [mode for mode in ("run_one", "split") if report.get("modes", {}).get(mode, {}).get("enabled", False)]


def comparison_regression_summary(comparisons: Mapping[str, Any]) -> dict[str, Any]:
    out: dict[str, Any] = {}
    for field in REQUIRED_FIELDS:
        data = comparisons.get(field, {})
        out[field] = {
            "actual_shape": data.get("actual_shape"),
            "expected_shape": data.get("expected_shape"),
            "shape_matches": data.get("shape_matches"),
            "allclose": data.get("allclose"),
            "max_abs": data.get("max_abs"),
            "mean_abs": data.get("mean_abs"),
            "tolerance_at_max": data.get("tolerance_at_max"),
            "location": data.get("location"),
        }
    return out


def backend_trust_evidence(report: Mapping[str, Any]) -> dict[str, Any]:
    backend = report.get("backend", {}) if isinstance(report.get("backend"), Mapping) else {}
    selected = backend.get("selected")
    requested = backend.get("requested")
    modes = {
        name: mode
        for name, mode in report.get("modes", {}).items()
        if name in {"run_one", "split"} and isinstance(mode, Mapping) and mode.get("enabled") is True
    }
    primitive = backend.get("primitive_evidence") if isinstance(backend.get("primitive_evidence"), Mapping) else {}
    capabilities = set(backend.get("capabilities", ())) if isinstance(backend.get("capabilities", ()), Iterable) else set()

    mode_gpu_ownership = {}
    for name, mode in modes.items():
        mode_backend = mode.get("backend") if isinstance(mode.get("backend"), Mapping) else {}
        gpu_solver = mode.get("gpu_solver") if isinstance(mode.get("gpu_solver"), Mapping) else {}
        gpu_solver_backend = gpu_solver.get("backend") if isinstance(gpu_solver.get("backend"), Mapping) else {}
        execution_path = str(mode.get("execution_path") or "").lower()
        mode_gpu_ownership[name] = {
            "backend_selected_gpu": mode_backend.get("selected") == "gpu",
            "execution_path_gpu_owned": "gpu" in execution_path and "cpu" not in execution_path,
            "gpu_solver_executed": gpu_solver.get("status") == "executed",
            "gpu_solver_backend_selected_gpu": gpu_solver_backend.get("selected") == "gpu",
            "used_cpu_fallback": mode_backend.get("used_cpu_fallback") is True or gpu_solver_backend.get("used_cpu_fallback") is True,
        }

    gpu_requested = requested == "gpu" or selected == "gpu"
    gpu_claim_checks = {
        "selected_gpu_backend": selected == "gpu",
        "provider_not_cpu": "cpu" not in str(backend.get("provider") or "").lower()
        and str(backend.get("provider") or "").lower() not in {"openfoam", "foam_stepper_cpu"},
        "device_not_host": str(backend.get("device") or "").lower() not in {"", "host", "cpu"},
        "no_cpu_fallback_capability": "no_cpu_fallback" in capabilities,
        "used_cpu_fallback_false": backend.get("used_cpu_fallback") is False,
        "cpu_fallback_disallowed": (backend.get("cpu_fallback") if isinstance(backend.get("cpu_fallback"), Mapping) else {}).get("allowed") is False,
        "primitive_not_acceptance": primitive.get("counts_as_full_gpu_rans_solver") is False,
        "full_solver_guard_named": bool(backend.get("full_solver_guard")),
        "modes_gpu_owned": bool(mode_gpu_ownership) and all(
            item["backend_selected_gpu"]
            and item["execution_path_gpu_owned"]
            and item["gpu_solver_executed"]
            and item["gpu_solver_backend_selected_gpu"]
            and not item["used_cpu_fallback"]
            for item in mode_gpu_ownership.values()
        ),
    }
    cpu_checks = {
        "selected_cpu_backend": selected == "cpu",
        "used_cpu_fallback_false": backend.get("used_cpu_fallback") is False,
    }
    checks = gpu_claim_checks if gpu_requested else cpu_checks
    return {
        "status": "trusted" if all(checks.values()) else "untrusted",
        "requested": requested,
        "selected": selected,
        "gpu_claimed": gpu_requested,
        "checks": checks,
        "mode_gpu_ownership": mode_gpu_ownership,
        "policy": "GPU-backed success requires GPU-owned run_one and split executions, no CPU fallback, and primitive evidence marked non-acceptance.",
    }


def build_verifier_evidence(report: Mapping[str, Any]) -> dict[str, Any]:
    modes = REQUIRED_EXECUTION_MODES
    prepared_identity = report.get("case_preparation", {}).get("prepared_case_identity", {})
    prepared_ok = bool(prepared_identity) and all(identity.get("matches_prepared") for identity in prepared_identity.values())

    mode_enabled = {
        mode: report.get("modes", {}).get(mode, {}).get("enabled") is True
        for mode in REQUIRED_EXECUTION_MODES
    }
    mesh_ok_by_mode = {
        mode: bool(report.get("modes", {}).get(mode, {}).get("mesh_comparison", {}).get("matches_oracle"))
        for mode in modes
    }
    comparison_ok_by_mode = {}
    observability_ok_by_mode = {}
    regression_modes: dict[str, Any] = {}
    for mode in modes:
        mode_report = report.get("modes", {}).get(mode, {})
        comparisons = mode_report.get("comparisons", {})
        comparison_ok_by_mode[mode] = all(
            comparisons.get(field, {}).get("shape_matches") is True and comparisons.get(field, {}).get("allclose") is True
            for field in REQUIRED_FIELDS
        )
        observability = report.get("stage_observability", {}).get(mode, {})
        groups = observability.get("groups", [])
        observability_ok_by_mode[mode] = bool(groups) and all(group.get("observed") is True for group in groups)
        regression_modes[mode] = {
            "enabled": mode_enabled[mode],
            "case_role": mode,
            "mesh_matches_oracle": mesh_ok_by_mode[mode],
            "observed_stage_groups": [group.get("name") for group in groups if group.get("observed")],
            "fields": comparison_regression_summary(comparisons),
        }

    timing = report.get("timing", {}) if isinstance(report.get("timing"), Mapping) else {}
    timing_criteria = timing.get("criteria", {}) if isinstance(timing.get("criteria"), Mapping) else {}
    if report.get("backend", {}).get("selected") == "gpu":
        timing_ok = (
            timing_criteria.get("openfoam_reference_timing") is True
            and timing_criteria.get("gpu_solver_wall_timing") is True
            and timing_criteria.get("gpu_profiler_timing") is True
        )
    else:
        timing_ok = timing.get("status") in {"reported", "not_reported"}
    backend_trust = backend_trust_evidence(report)
    intermediate_artifacts = report.get("intermediate_artifacts", {}) if isinstance(report.get("intermediate_artifacts"), Mapping) else {}



    criteria = {
        "required_execution_modes": all(mode_enabled.values()),
        "prepared_case_identity": prepared_ok,
        "mesh_identity": all(mesh_ok_by_mode.values()),
        "field_comparisons": all(comparison_ok_by_mode.values()),
        "stage_observability": all(observability_ok_by_mode.values()),
        "backend_selected": report.get("backend", {}).get("selected") is not None and backend_trust.get("status") == "trusted",
        "timing_reported": timing_ok,
        "intermediate_artifacts_complete": intermediate_artifacts.get("status") == "complete",
        "intermediate_artifacts_reported": bool(intermediate_artifacts.get("families")),
    }
    passed = all(criteria.values())
    oracle_mesh = report.get("mesh_identity", {}).get("oracle", {})
    prepared_digest = prepared_identity.get("prepared", {}).get("files_sha256")
    return {
        "passed": passed,
        "status_basis": "passed all verifier evidence criteria" if passed else "one or more verifier evidence criteria failed",
        "criteria": criteria,
        "mode_criteria": {
            "mesh_identity": mesh_ok_by_mode,
            "field_comparisons": comparison_ok_by_mode,
            "stage_observability": observability_ok_by_mode,
            "backend_trust": backend_trust,
            "intermediate_artifacts": {
                "status": intermediate_artifacts.get("status"),
                "families": {
                    family.get("name"): family.get("status")
                    for family in intermediate_artifacts.get("families", [])
                    if isinstance(family, Mapping)
                },
                "first_actionable_family": (
                    intermediate_artifacts.get("first_actionable_family", {}).get("name")
                    if isinstance(intermediate_artifacts.get("first_actionable_family"), Mapping)
                    else None
                ),
            },
        },
        "entrypoint": report.get("harness_entrypoint", {}),
        "regression_summary": {
            "schema_version": report.get("harness", {}).get("schema_version"),
            "source_case": report.get("problem_boundary", {}).get("airfrans_simulation"),
            "prepared_case_files_sha256": prepared_digest,
            "backend_selected": report.get("backend", {}).get("selected"),
            "backend_trust_status": backend_trust.get("status"),
            "differentiability_status": report.get("differentiability", {}).get("status"),
            "mesh_topology_sha256": oracle_mesh.get("topology_sha256"),
            "first_divergence_summary": report.get("first_divergence_summary"),
            "differential_trace_status": report.get("differential_trace", {}).get("status") if isinstance(report.get("differential_trace"), Mapping) else None,
            "differential_trace_first_divergence": report.get("differential_trace", {}).get("first_divergence") if isinstance(report.get("differential_trace"), Mapping) else None,
            "mesh_geometry_sha256": oracle_mesh.get("geometry_sha256"),
            "modes": regression_modes,
        },
    }



def base_report(args: argparse.Namespace) -> dict[str, Any]:
    report_path = args.report if args.report is not None else args.work / "verifier_report.json"
    return {
        "harness": {
            "name": "airfrans_stepper_verifier",
            "spec": "VERIFIER_HARNESS_SPEC.md",
            "schema_version": 1,
            "scope": AIRFRANS_VERIFIER_SCOPE,
        },
        "status": "running",
        "failure": None,
        "paths": {
            "root": ROOT,
            "source": args.source,
            "work": args.work,
            "report": report_path,
        },
        "harness_entrypoint": {
            "command": "uv run python scripts/verify_airfrans_stepper.py --work <work-dir> --report <report.json>",
            "arguments": {
                "source": args.source,
                "work": args.work,
                "report": report_path,
                "backend": args.backend,
                "skip_split": args.skip_split,
                "rtol": args.rtol,
                "atol": args.atol,
                "trace_substitute": list(args.trace_substitute),
                "trace_substituted_artifact": args.trace_substituted_artifact,
            },
            "prepares_case": True,
            "runs_oracle": True,
            "runs_migrated_path": True,
            "writes_report": True,
            "scope": AIRFRANS_VERIFIER_SCOPE,
        },
        "source_policy": "raw AirfRANS source is read-only; all solver runs use prepared work-directory copies",
        "tolerances": {
            "rtol": args.rtol,
            "atol": args.atol,
            "policy": "CPU stepper parity with the repository OpenFOAM oracle must pass np.allclose for every required field.",
        },
        "comparison_policy": {
            "function": "numpy.allclose(actual, expected, rtol, atol, equal_nan=False)",
            "scope": "internal arrays for U, p, phi, nut, k, and omega",
            "backend_tolerances": {
                "cpu": {
                    "rtol": args.rtol,
                    "atol": args.atol,
                    "reason": "strict CPU parity against the repository OpenFOAM oracle",
                },
                "gpu": {
                    "rtol": args.rtol,
                    "atol": args.atol,
                    "reason": "full GPU RANS backend must own solver execution and satisfy OpenFOAM oracle parity; primitive GPU diagnostics are not acceptance",
                    "availability_checked_during_backend_selection": True,
                    "full_solver_guard": FULL_GPU_RANS_GUARD,
                    "primitive_evidence": {
                        "name": GPU_PRIMITIVE_NAME,
                        "path": GPU_PRIMITIVE_PROOF,
                        "counts_as_full_gpu_rans_solver": False,
                    },
                },
            },
        },
        "backend": {
            "requested": args.backend,
            "selected": None,
            "failure_category": BACKEND_FAILURE,
            "execution_scope": AIRFRANS_VERIFIER_SCOPE,
        },
        "required_fields": {
            "primary": list(PRIMARY_FIELDS),
            "turbulence": list(TURBULENCE_FIELDS),
            "all": list(REQUIRED_FIELDS),
        },
        "problem_boundary": {
            "airfrans_simulation": args.source.name,
            "selected_source_case": args.source,
            "prepared_case_policy": "prepare one reproducible OpenFOAM v14 case and clone it into isolated oracle/run_one/split execution cases",
            "oracle_output": {
                "case_role": "oracle",
                "producer": "foamRun -solver incompressibleFluid -noFunctionObjects",
                "time": "1",
                "fields": list(REQUIRED_FIELDS),
            },
            "repository_outputs": {
                "run_one": "migrated run_one path selected by the requested backend; gpu requests must use the GPU solver contract and reject CPU fallback",
                "split": "split solver stages selected by the requested backend",
                "time": "1",
                "fields": list(REQUIRED_FIELDS),
            },
            "mesh_identity_evidence": [
                "n_points",
                "n_faces",
                "n_internal_faces",
                "n_cells",
                "patches",
                "topology_sha256",
                "geometry_sha256",
            ],
        },
        "commands": [],
        "case_preparation": {},
        "mesh_identity": {},
        "modes": {
            "run_one": {"enabled": True},
            "split": {"enabled": not args.skip_split},
        },
        "tracked_turbulence_fields": {},
        "comparison_mismatches": [],
        "comparison_attribution": [],
        "failure_diagnostics": {
            "status": "not_evaluated",
            "first_blocking_point": None,
            "downstream_unreliable": [],
        },
        "verifier_evidence": {
            "passed": False,
            "status_basis": "not_evaluated",
        },
        "timing": {
            "status": "not_evaluated",
        },
        "differentiability": {
            "status": "not_evaluated",
            "reason": "backend has not been selected yet",
        },
        "state_exports": {},
        "diagnostic_artifacts": {
            "schema_version": NUMERIC_ARTIFACT_SCHEMA_VERSION,
            "status": "not_evaluated",
            "root": diagnostic_artifact_root(args),
        },
        "differential_trace": {
            "schema_version": 1,
            "status": "not_evaluated",
            "comparison_basis": "reference_split_vs_split_npz_checkpoints",
            "substitution": {
                "requested": list(args.trace_substitute),
                "status": "not_requested" if not args.trace_substitute else "pending_reference_artifact",
                "cpu_fallback_allowed": False,
            },
        },
        "gpu_solver_inputs": {
            "schema_version": GPU_INPUT_SCHEMA_VERSION,
            "status": "not_requested",
        },
        "stage_observability": {
            "contract": json_ready(STAGE_OBSERVABILITY_GROUPS),
        },
    }


def run_harness(args: argparse.Namespace, report: dict[str, Any]) -> None:
    try:
        prepared = prepare_work_cases(args.source, args.work, include_split=not args.skip_split)
    except Exception as exc:
        raise HarnessError(
            LOADING_FAILURE,
            "prepare_cases",
            "failed to prepare reproducible AirfRANS work cases",
            details={"source": args.source, "work": args.work},
            cause=exc,
        ) from exc

    oracle_case = prepared["oracle_case"]
    run_one_case = prepared["run_one_case"]
    split_case = prepared["split_case"]
    report["case_preparation"] = json_ready(prepared)

    identity_mismatches = {
        role: identity
        for role, identity in prepared["prepared_case_identity"].items()
        if not identity["matches_prepared"]
    }
    if identity_mismatches:
        raise HarnessError(
            LOADING_FAILURE,
            "prepare_case_identity",
            "prepared execution case copies differ from the reproducible baseline case",
            details={"prepared_case_identity": identity_mismatches},
        )


    try:
        backend = select_execution_backend(args.backend)
    except BackendExecutionError as exc:
        raise exc.to_harness_error() from exc
    report["backend"] = json_ready(backend)
    report["differentiability"] = json_ready(differentiability_report(backend, case_name=args.source.name))
    trace_substitutions: dict[str, np.ndarray] | None = None
    if backend.get("selected") == "gpu" and split_case is not None:
        reference_summary, reference_report = run_reference_split_diagnostic(args)
        report["reference_split_diagnostic"] = json_ready(reference_summary)
        if reference_report is not None:
            reference_split = reference_report.get("modes", {}).get("split", {})
            report["modes"]["reference_split"] = {
                **reference_split,
                "role": "separate_process_openfoam_split_reference_for_gpu_artifact_comparison",
                "source_report": reference_summary["report"],
            }
            reference_state_exports = reference_report.get("state_exports", {})
            if "split" in reference_state_exports:
                report["state_exports"]["reference_split"] = reference_state_exports["split"]
            if "split_matrices" in reference_state_exports:
                report["state_exports"]["reference_split_matrices"] = reference_state_exports["split_matrices"]
            reference_artifacts = reference_report.get("diagnostic_artifacts", {})
            if isinstance(reference_artifacts, Mapping) and isinstance(reference_artifacts.get("split"), Mapping):
                report["diagnostic_artifacts"]["reference_split"] = reference_artifacts["split"]
    if args.trace_substitute:
        trace_substitutions, trace_substitution_manifest = load_trace_substitutions(report, args.trace_substitute)
        report["differential_trace"]["substitution"] = {
            **report["differential_trace"].get("substitution", {}),
            **trace_substitution_manifest,
            "cpu_fallback_allowed": False,
            "mode": "reference_checkpoint_replay_before_downstream_gpu_execution",
        }
    if backend.get("selected") == "gpu":
        try:
            foam = import_foam()
        except Exception as exc:
            raise HarnessError(
                LOADING_FAILURE,
                "import_foam_stepper",
                "failed to import foam_stepper before GPU input preparation",
                cause=exc,
            ) from exc
        gpu_input_stepper = make_stepper(foam, run_one_case, "gpu_inputs")
        report["mesh_identity"]["gpu_inputs"] = read_mesh_identity(gpu_input_stepper, "gpu_inputs")
        try:
            gpu_inputs = prepare_gpu_solver_inputs(foam, gpu_input_stepper, run_one_case, prepared, backend)
        except BackendExecutionError as exc:
            raise exc.to_harness_error() from exc
        report["gpu_solver_inputs"] = json_ready(gpu_inputs)
        try:
            gpu_stage_smoke = run_gpu_solver_stage_smoke(gpu_input_stepper, backend, run_one_case)
        except BackendExecutionError as exc:
            raise exc.to_harness_error() from exc
        report["gpu_solver_stages"] = json_ready(gpu_solver_stage_report(gpu_stage_smoke))
        gpu_contract_blocker = gpu_solver_contract_blocker(backend, run_one_case)
        if gpu_contract_blocker is not None:
            record_gpu_contract_failure(report, gpu_contract_blocker, gpu_inputs)
    report["commands"].append(
        run_openfoam_command(["checkMesh", "-case", str(oracle_case), "-constant"], log_path=args.work / "checkMesh.log", timeout=120)
    )
    report["commands"].append(
        run_openfoam_command(
            ["foamRun", "-case", str(oracle_case), "-solver", "incompressibleFluid", "-noFunctionObjects"],
            log_path=args.work / "foamRun_oracle.log",
            timeout=600,
        )
    )

    try:
        patch_start_from_latest(oracle_case)
    except Exception as exc:
        raise HarnessError(
            LOADING_FAILURE,
            "select_oracle_latest_time",
            "failed to configure oracle case to load latestTime output",
            details={"case": oracle_case},
            cause=exc,
        ) from exc

    try:
        foam = import_foam()
    except Exception as exc:
        raise HarnessError(
            LOADING_FAILURE,
            "import_foam_stepper",
            "failed to import foam_stepper in the repository OpenFOAM environment",
            cause=exc,
        ) from exc

    oracle_stepper = make_stepper(foam, oracle_case, "oracle")
    oracle_mesh = read_mesh_identity(oracle_stepper, "oracle")
    oracle_fields = read_fields(oracle_stepper, "oracle")
    report["mesh_identity"]["oracle"] = oracle_mesh
    report["modes"]["oracle"] = {
        "case": oracle_case,
        "fields": selected_field_summaries(oracle_fields),
    }
    report["tracked_turbulence_fields"]["oracle"] = visible_turbulence_fields(oracle_fields)
    oracle_state = export_solver_state_checked(foam, oracle_stepper, "oracle", oracle_fields)
    report["state_exports"]["oracle"] = foam.describe_solver_state(oracle_state)



    run_one_stepper = make_stepper(foam, run_one_case, "run_one")
    run_one_mesh = read_mesh_identity(run_one_stepper, "run_one")
    report["mesh_identity"]["run_one"] = run_one_mesh
    mesh_report, mesh_mismatches = compare_mesh_identity("run_one", run_one_mesh, oracle_mesh)
    report["modes"]["run_one"]["mesh_comparison"] = mesh_report

    try:
        backend_run = run_backend_iteration(run_one_stepper, backend, run_one_case)
        run_one_result = backend_run["result"]
        run_one_fields = backend_run["fields"]
    except BackendExecutionError as exc:
        raise exc.to_harness_error() from exc
    except Exception as exc:
        raise HarnessError(
            STEPPER_FAILURE,
            "run_one_pimple_iteration",
            "full one-iteration backend execution failed",
            details={"case": run_one_case, "backend": backend},
            cause=exc,
        ) from exc

    run_one_stages = graph_stage_summaries(run_one_result)
    run_one_observability = stage_observability_report(
        "run_one",
        run_one_stages,
        evidence={
            "one_iteration_execution": True,
            "compared_fields": list(REQUIRED_FIELDS),
            "turbulence_fields": visible_turbulence_fields(run_one_fields),
            "backend": backend,
        },
    )
    report["stage_observability"]["run_one"] = run_one_observability

    run_one_comparisons, run_one_mismatches = compare_fields(
        "run_one",
        run_one_fields,
        oracle_fields,
        rtol=args.rtol,
        atol=args.atol,
        stage_graph=[stage["name"] for stage in run_one_stages],
        mesh_context=run_one_stepper.mesh(),
    )
    report["modes"]["run_one"].update(
        {
            "case": run_one_case,
            "execution_path": backend_run["execution_path"],
            "backend": backend_run["backend"],
            "stage": transform_summary(run_one_result),
            "graph": [stage["name"] for stage in run_one_stages],
            "fields": selected_field_summaries(run_one_fields),
            "comparisons": run_one_comparisons,
            "observability": run_one_observability,
        }
    )
    if backend.get("selected") == "gpu":
        report["modes"]["run_one"]["gpu_solver"] = json_ready(backend_run.get("gpu_solver", {}))
    report["tracked_turbulence_fields"]["run_one"] = visible_turbulence_fields(run_one_fields)
    run_one_state = export_solver_state_checked(foam, run_one_stepper, "run_one", run_one_fields)
    report["state_exports"]["run_one"] = foam.describe_solver_state(run_one_state)


    mismatches = mesh_mismatches + run_one_mismatches

    split_artifact_mesh = None

    if split_case is not None:
        split_stepper = make_stepper(foam, split_case, "split")
        split_mesh = read_mesh_identity(split_stepper, "split")
        report["mesh_identity"]["split"] = split_mesh
        split_mesh_report, split_mesh_mismatches = compare_mesh_identity("split", split_mesh, oracle_mesh)
        report["modes"]["split"]["mesh_comparison"] = split_mesh_report

        split_artifact_path = diagnostic_artifact_root(args) / "split.npz"
        split_artifact_mesh = split_stepper.mesh()
        if backend.get("selected") == "gpu":
            split_result = run_gpu_split_iteration(foam, split_stepper, backend, split_case, diagnostic_artifact_path=split_artifact_path, trace_substitutions=trace_substitutions)
        else:
            split_result = run_split_iteration(foam, split_stepper, diagnostic_artifact_path=split_artifact_path)
        report["stage_observability"]["split"] = split_result["observability"]
        split_fields = split_result["fields"]
        split_comparisons, split_mismatches = compare_fields(
            "split",
            split_fields,
            oracle_fields,
            rtol=args.rtol,
            atol=args.atol,
            stage_graph=split_result["graph"],
            mesh_context=split_stepper.mesh(),
        )
        report["modes"]["split"].update(
            {
                "case": split_case,
                "graph": split_result["graph"],
                "stages": split_result["stages"],
                "momentum_terms": split_result["momentum_terms"],
                "UEqn": split_result["UEqn"],
                "pEqn": split_result["pEqn"],
                "fields": selected_field_summaries(split_fields),
                "comparisons": split_comparisons,
                "observability": split_result["observability"],
                "execution_path": split_result.get("execution_path"),
                "backend": split_result.get("backend"),
                "gpu_solver": split_result.get("gpu_solver"),
            }
        )
        report["tracked_turbulence_fields"]["split"] = visible_turbulence_fields(split_fields)
        report["state_exports"]["split"] = split_result["state_summary"]
        report["state_exports"]["split_matrices"] = split_result["matrix_states"]
        if isinstance(split_result.get("diagnostic_artifacts"), Mapping):
            report["diagnostic_artifacts"]["split"] = split_result["diagnostic_artifacts"]
            report["diagnostic_artifacts"]["status"] = "available"
        mismatches.extend(split_mesh_mismatches)
        mismatches.extend(split_mismatches)

    comparison_attribution = [mismatch["attribution"] for mismatch in mismatches if "attribution" in mismatch]
    report["comparison_attribution"] = json_ready(comparison_attribution)
    report["comparison_mismatches"] = json_ready(mismatches)
    report["artifact_comparisons"] = json_ready(build_artifact_comparisons(report, mesh_context=split_artifact_mesh))
    report["differential_trace"] = build_differential_trace(report, mesh_context=split_artifact_mesh)
    if args.trace_substituted_artifact is not None and args.trace_substitute:
        from foam_stepper.differential_trace import load_npz_artifacts, materialize_substituted_artifacts

        reference_path = diagnostic_artifact_path(report, "reference_split")
        candidate_path = diagnostic_artifact_path(report, "split")
        reference_artifacts = load_npz_artifacts(reference_path)
        candidate_artifacts = load_npz_artifacts(candidate_path)
        if reference_artifacts is not None and candidate_artifacts is not None:
            report["differential_trace"]["substituted_artifact"] = materialize_substituted_artifacts(
                reference_artifacts=reference_artifacts,
                candidate_artifacts=candidate_artifacts,
                checkpoints=args.trace_substitute,
                output_path=args.trace_substituted_artifact,
            )
        else:
            report["differential_trace"]["substituted_artifact"] = {
                "schema_version": 1,
                "status": "failed",
                "requested": list(args.trace_substitute),
                "reason": "reference_or_candidate_artifact_missing",
                "reference_path": str(reference_path) if reference_path is not None else None,
                "candidate_path": str(candidate_path) if candidate_path is not None else None,
            }
    report["failure_diagnostics"] = json_ready(build_failure_diagnostics(mismatches, report))
    report["timing"] = json_ready(build_timing_evidence(report))
    report["intermediate_artifacts"] = json_ready(intermediate_artifact_family_report(report))
    report["first_divergence_summary"] = json_ready(first_divergence_summary(report))
    report["verifier_evidence"] = json_ready(build_verifier_evidence(report))
    if mismatches:
        raise HarnessError(
            COMPARISON_FAILURE,
            "compare_oracle_stepper_outputs",
            f"{len(mismatches)} verifier comparison mismatch(es) observed",
            details={
                "failure_kind": GPU_NUMERICAL_MISMATCH if backend.get("selected") == "gpu" else COMPARISON_FAILURE,
                "mismatches": mismatches,
                "likely_causes": comparison_attribution,
                "first_divergence_summary": report["first_divergence_summary"],
                "diagnostics": report["failure_diagnostics"],
                "differential_trace": report["differential_trace"],
            },
        )

    if report["verifier_evidence"].get("passed") is not True:
        raise HarnessError(
            REQUIRED_EVIDENCE_FAILURE,
            "verifier_evidence",
            "required verifier evidence is incomplete or failed",
            details={"verifier_evidence": report["verifier_evidence"], "first_divergence_summary": report["first_divergence_summary"]},
        )


def write_report(report: Mapping[str, Any], path: Path) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    path.write_text(json.dumps(json_ready(report), indent=2, sort_keys=True, allow_nan=False) + "\n")


def fmt_sci(value: Any) -> str:
    if value is None:
        return "none"
    try:
        return f"{float(value):.3e}"
    except (TypeError, ValueError):
        return str(value)


def format_location(location: Mapping[str, Any] | None) -> str:
    if not location:
        return "none"
    component = location.get("component_index")
    component_text = "" if component is None else f" component={component}"
    return f"{location.get('entity_kind')}[{location.get('entity_index')}]{component_text}"


def print_human_summary(report: Mapping[str, Any]) -> None:
    status = report.get("status")
    print(f"airfrans_stepper_verification {status}")
    print(f"source={report['paths']['source']}")
    print(f"work={report['paths']['work']}")
    print(f"report={report['paths']['report']}")
    entrypoint = report.get("harness_entrypoint") or {}
    print(f"harness_entrypoint={entrypoint.get('command')}")
    print(f"harness_scope={entrypoint.get('scope') or report.get('harness', {}).get('scope')}")
    print(f"rtol={report['tolerances']['rtol']} atol={report['tolerances']['atol']}")
    backend = report.get("backend", {})
    print(f"backend_requested={backend.get('requested')} backend_selected={backend.get('selected')}")
    differentiability = report.get("differentiability") or {}
    if differentiability:
        claims = differentiability.get("claims") or []
        checks = differentiability.get("checks") or []
        print(f"differentiability_status={differentiability.get('status')}")
        print(f"differentiability_claims={len(claims)} non_differentiable_boundaries={len(differentiability.get('non_differentiable') or [])}")
        print(f"differentiability_checks={[check.get('status') for check in checks]}")
    evidence = report.get("verifier_evidence") or {}
    if evidence:
        print(f"evidence_passed={evidence.get('passed')} evidence_status_basis={evidence.get('status_basis')}")
        print(f"evidence_criteria={evidence.get('criteria')}")
    artifacts = report.get("intermediate_artifacts") or {}
    if artifacts:
        families = artifacts.get("families") or []
        print(f"intermediate_artifacts_status={artifacts.get('status')}")
        print(f"intermediate_artifact_families={[(family.get('name'), family.get('status')) for family in families]}")
        first_family = artifacts.get("first_actionable_family") or {}
        print(f"intermediate_first_actionable_family={first_family.get('name')}")
    trace = report.get("differential_trace") or {}
    if trace:
        first_trace = trace.get("first_divergence") or {}
        substitution = trace.get("substitution") or {}
        largest = first_trace.get("largest_difference") or {}
        print(f"differential_trace_status={trace.get('status')} checkpoint_count={trace.get('checkpoint_count')} failed={trace.get('failed_count')} missing={trace.get('missing_count')}")
        print(f"differential_trace_first_checkpoint={first_trace.get('name')} max_abs={fmt_sci(first_trace.get('max_abs'))} location={format_location((largest.get('location') or {}) if isinstance(largest, Mapping) else {})}")
        print(f"differential_trace_substitution_status={substitution.get('status')} loaded={substitution.get('loaded')}")
    compact = report.get("first_divergence_summary") or {}
    if compact:
        print(f"first_divergence_summary={compact}")
        print(f"first_divergence_status={compact.get('status')}")
        print(f"first_divergence_target={compact.get('first_target')}")
        print(f"first_divergence_evidence_path={compact.get('evidence_path')}")
        print(f"first_divergence_artifact_family={compact.get('artifact_family')}")

    if status != "passed":
        failure = report.get("failure") or {}
        print(f"failure_category={failure.get('category')}")
        print(f"failure_step={failure.get('step')}")
        print(f"failure_message={failure.get('message')}")
        timing = report.get("timing") or {}
        if timing.get("status") == "reported":
            openfoam = timing.get("openfoam_reference") or {}
            gpu_solver = timing.get("gpu_solver") or {}
            run_one_timing = ((gpu_solver.get("run_one") or {}).get("timing") or {})
            run_one_profiler = ((gpu_solver.get("run_one") or {}).get("profiler") or {})
            print(f"openfoam_oracle_duration_seconds={openfoam.get('oracle_solver_duration_seconds')}")
            print(f"gpu_solver_run_one_wall_seconds={run_one_timing.get('gpu_solver_wall_seconds')}")
            print(f"gpu_solver_run_one_device_time_ms={run_one_profiler.get('device_time_ms_total')}")
        diagnostics = ((failure.get("details") or {}).get("diagnostics") or report.get("failure_diagnostics") or {})
        first_blocker = diagnostics.get("first_blocking_point") or {}
        if first_blocker:
            field_context = first_blocker.get("field_context") or {}
            error_context = field_context.get("error") or {}
            largest = field_context.get("largest_difference") or {}
            print(f"first_blocking_mode={first_blocker.get('mode')}")
            print(f"first_blocking_stage_group={first_blocker.get('stage_group')}")
            print(f"first_blocking_stages={first_blocker.get('stage_names')}")
            first_target = first_blocker.get("first_target") or {}
            print(f"first_debug_target={first_target.get('single_next_target')}")
            print(f"first_debug_target_kind={first_target.get('kind')}")
            print(f"first_debug_target_evidence_path={first_target.get('evidence_path')}")
            print(f"first_debug_target_why={first_target.get('why')}")
            print(f"first_blocking_field={field_context.get('field')}")
            print(f"first_blocking_reason={field_context.get('reason')}")
            print(f"first_blocking_max_abs={fmt_sci(error_context.get('max_abs'))}")
            print(f"first_blocking_location={format_location(largest.get('location'))}")
            local_context = field_context.get("local_context") or {}
            print(f"first_blocking_local_context={local_context}")
            print(f"downstream_unreliable_count={len(diagnostics.get('downstream_unreliable') or [])}")
        mismatches = (failure.get("details") or {}).get("mismatches") or report.get("comparison_mismatches") or []
        if mismatches:
            first = mismatches[0]
            attribution = first.get("attribution") or {}
            print(f"comparison_failure_mode={first.get('mode')}")
            print(f"comparison_failure_field={first.get('field')}")
            print(f"comparison_failure_reason={first.get('reason')}")
            print(f"comparison_failure_likely_stage_group={attribution.get('likely_stage_group')}")
            print(f"comparison_failure_likely_stages={attribution.get('likely_stages')}")
            print(f"comparison_failure_max_abs={fmt_sci(first.get('max_abs'))}")
            print(f"comparison_failure_location={format_location(first.get('location'))}")
        return

    oracle_mesh = report["mesh_identity"]["oracle"]
    patch_names = [patch["name"] for patch in oracle_mesh["patches"]]
    print(f"oracle_case={report['case_preparation']['oracle_case']}")
    prepared_identity = report["case_preparation"].get("prepared_case_identity", {})
    prepared_digest = prepared_identity.get("prepared", {}).get("files_sha256")
    print(f"prepared_case={report['case_preparation']['prepared_case']}")
    print(f"prepared_case_files_sha256={prepared_digest}")
    print(f"run_one_case={report['modes']['run_one']['case']}")
    if report["modes"].get("split", {}).get("enabled"):
        print(f"split_case={report['modes']['split']['case']}")
    print(f"mesh_cells={oracle_mesh['n_cells']}")
    print(f"mesh_internal_faces={oracle_mesh['n_internal_faces']}")
    print(f"mesh_patches={patch_names}")
    print(f"mesh_topology_sha256={oracle_mesh['topology_sha256']}")

    for mode in ("run_one", "split"):
        mode_report = report["modes"].get(mode, {})
        if not mode_report.get("enabled", False):
            continue
        for field_name in REQUIRED_FIELDS:
            data = mode_report.get("comparisons", {}).get(field_name)
            if not data:
                continue
            print(
                f"mode={mode} field={field_name} actual_shape={tuple(data['actual_shape'])} "
                f"expected_shape={tuple(data['expected_shape'])} max_abs={fmt_sci(data.get('max_abs'))} "
                f"location={format_location(data.get('location'))} allclose={data['allclose']}"
            )
        graph = mode_report.get("graph")
        if graph:
            print(f"{mode}_graph={graph}")
        observability = mode_report.get("observability")
        if observability:
            observed_groups = [group["name"] for group in observability.get("groups", []) if group.get("observed")]
            print(f"{mode}_observed_stage_groups={observed_groups}")

    split_report = report["modes"].get("split", {})
    if split_report.get("enabled") and "UEqn" in split_report and "pEqn" in split_report:
        print(f"split_momentum_terms={split_report.get('momentum_terms')}")
        print(f"split_UEqn_diag_shape={tuple(split_report['UEqn']['diag']['shape'])}")
        print(f"split_pEqn_diag_shape={tuple(split_report['pEqn']['diag']['shape'])}")

    print(f"visible_turbulence_fields={report['tracked_turbulence_fields']}")


def parse_args(argv: list[str] | None = None) -> argparse.Namespace:
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("--source", type=Path, default=DEFAULT_SOURCE)
    parser.add_argument("--work", type=Path, default=DEFAULT_WORK)
    parser.add_argument("--report", type=Path, default=None, help="JSON report path; defaults to WORK/verifier_report.json")
    parser.add_argument("--diagnostic-artifacts", type=Path, default=None, help="Directory for focused numeric split-stage NPZ artifacts; defaults to WORK/diagnostic_artifacts")
    parser.add_argument("--rtol", type=float, default=1e-8)
    parser.add_argument("--atol", type=float, default=1e-8)
    parser.add_argument("--backend", choices=BACKEND_CHOICES, default="auto", help="Backend request; gpu means the full RANS solver backend and currently fails until that backend exists")
    parser.add_argument("--skip-split", action="store_true", help="Skip the explicit split-step parity check for local debugging")
    parser.add_argument("--trace-substitute", action="append", default=[], help="Replace this GPU split checkpoint with the reference_split artifact before downstream replay; may be repeated")
    parser.add_argument("--trace-substituted-artifact", type=Path, default=None, help="Optional NPZ path for an artifact-only candidate copy with requested trace substitutions applied")
    return parser.parse_args(argv)


def main(argv: list[str] | None = None) -> int:
    args = parse_args(argv)
    args.source = args.source.resolve()
    args.work = args.work.resolve()
    if args.report is None:
        args.report = args.work / "verifier_report.json"
    else:
        args.report = args.report.resolve()
    if args.trace_substituted_artifact is not None:
        args.trace_substituted_artifact = args.trace_substituted_artifact.resolve()

    report = base_report(args)
    exit_code = 0
    try:
        run_harness(args, report)
        evidence = report.get("verifier_evidence")
        if not isinstance(evidence, Mapping) or evidence.get("passed") is not True:
            raise HarnessError(
                INTERNAL_FAILURE,
                "verifier_evidence",
                "verifier completed without passing evidence criteria",
                details={"verifier_evidence": evidence},
            )
        report["status"] = "passed"
    except HarnessError as exc:
        report["status"] = "failed"
        report["failure"] = exc.to_dict()
        exit_code = EXIT_CODES.get(exc.category, EXIT_CODES[INTERNAL_FAILURE])
    except Exception as exc:  # pragma: no cover - keeps CLI failures categorized in the report.
        report["status"] = "failed"
        report["failure"] = json_ready(
            {
                "category": INTERNAL_FAILURE,
                "step": "run_harness",
                "message": str(exc),
                "cause": {"type": type(exc).__name__, "message": str(exc)},
                "traceback": traceback.format_exc(),
            }
        )
        exit_code = EXIT_CODES[INTERNAL_FAILURE]

    try:
        write_report(report, args.report)
    except Exception as exc:
        report["status"] = "failed"
        report["failure"] = json_ready(
            {
                "category": INTERNAL_FAILURE,
                "step": "write_report",
                "message": f"failed to write verifier report: {args.report}",
                "cause": {"type": type(exc).__name__, "message": str(exc)},
            }
        )
        exit_code = EXIT_CODES[INTERNAL_FAILURE]

    print_human_summary(json_ready(report))
    return exit_code


if __name__ == "__main__":
    raise SystemExit(main())