stash: partial completion of RANS GPU-based implementation before better observability harness build

2026-07-26 10:47:36 +04:00 · 2026-07-26 10:47:36 +04:00 · 6ea17a78cd
commit 6ea17a78cd
parent 1bdbaa9572
17 changed files with 6406 additions and 64 deletions
--- a/.gitignore
+++ b/.gitignore
@ -2,6 +2,7 @@
 .loop/
 **/__pycache__
 python/src/foam_stepper.egg-info/
+tmp/

 # foreign
 OpenFOAM-14/
--- a/python/src/foam_stepper/init.py
+++ b/python/src/foam_stepper/init.py
@ -3,7 +3,7 @@
 from __future__ import annotations

 from pathlib import Path
-from typing import Any
+from typing import Any, Iterable

 from ._runtime import configure_openfoam_environment

@ -23,6 +23,15 @@ from .types import (
    TransformResult,
    _wrap_value,
 )
+from .state import (
+    SolverStateValidationError,
+    describe_matrix_state,
+    describe_solver_state,
+    export_matrix_state,
+    export_solver_state,
+    validate_matrix_state,
+    validate_solver_state,
+)

 OpenFoamError = _native.OpenFoamError

@ -89,6 +98,12 @@ class SimpleStepper:
    def state(self) -> dict[str, Any]:
        return dict(self._native.state())

+    def export_state(self, *, required_fields: Iterable[str] = ()) -> dict[str, Any]:
+        return export_solver_state(self.mesh(), self.fields(), required_fields=required_fields)
+
+    def state_summary(self, *, required_fields: Iterable[str] = ()) -> dict[str, Any]:
+        return describe_solver_state(self.export_state(required_fields=required_fields))
+
    def control_dict(self) -> str:
        return self._native.control_dict()

@ -201,5 +216,12 @@ __all__ = [
    "SolveResult",
    "SourceLocation",
    "TransformResult",
+    "SolverStateValidationError",
+    "describe_matrix_state",
+    "describe_solver_state",
+    "export_matrix_state",
+    "export_solver_state",
+    "validate_matrix_state",
+    "validate_solver_state",
    "version",
 ]
--- a/python/src/foam_stepper/gpu/init.py
+++ b/python/src/foam_stepper/gpu/init.py
@ -0,0 +1,5 @@
+"""Reusable GPU backend modules for foam_stepper."""
+
+from __future__ import annotations
+
+__all__ = ["backend", "constants", "kernels", "linear_solve"]
--- a/python/src/foam_stepper/gpu/backend.py
+++ b/python/src/foam_stepper/gpu/backend.py
--- a/python/src/foam_stepper/gpu/constants.py
+++ b/python/src/foam_stepper/gpu/constants.py
@ -0,0 +1,81 @@
+"""GPU RANS backend constants."""
+
+from __future__ import annotations
+
+FULL_GPU_RANS_GUARD = "scripts/verify_gpu_rans_solver.sh"
+GPU_PRIMITIVE_PROOF = "scripts/verify_gpu_algorithm.sh"
+GPU_PRIMITIVE_NAME = "cell_flux_imbalance"
+
+PRIMARY_FIELDS = ("U", "p", "phi")
+TURBULENCE_FIELDS = ("nut", "k", "omega")
+REQUIRED_FIELDS = (*PRIMARY_FIELDS, *TURBULENCE_FIELDS)
+
+GPU_BACKEND_PROVIDER = "quadrants_cuda_rans_solver"
+GPU_RUNTIME_UNAVAILABLE = "gpu_runtime_unavailable"
+GPU_KERNELS_MISSING = "missing_gpu_solver_kernels"
+GPU_STAGE_UNSUPPORTED = "unsupported_gpu_solver_stage"
+GPU_NUMERICAL_MISMATCH = "gpu_numerical_mismatch"
+GPU_STAGE_CAPABILITY_PREFIX = "solver_stage_contract:"
+GPU_INPUT_SCHEMA_VERSION = 1
+DEFAULT_LAMINAR_NU = 1.5e-5
+DEFAULT_MOMENTUM_RELAXATION_ALPHA = 0.9
+DEFAULT_MOMENTUM_PBICGSTAB_ITERATIONS = 50
+DEFAULT_MOMENTUM_PBICGSTAB_RESIDUAL_TOLERANCE_SQUARED = 1.0e-16
+DEFAULT_PRESSURE_CG_ITERATIONS = 300
+DEFAULT_SIMPLE_CONSISTENT_RATU_FACTOR = 10.0
+DEFAULT_OMEGA_WALL_BETA1 = 0.075
+
+STAGE_OBSERVABILITY_GROUPS = (
+    {
+        "name": "momentum_assembly",
+        "split_stages": ("assemble_momentum_terms", "assemble_UEqn"),
+        "run_one_stages": ("assemble_UEqn",),
+        "split_outputs": {
+            "assemble_momentum_terms": ("terms",),
+            "assemble_UEqn": ("UEqn",),
+        },
+    },
+    {
+        "name": "pressure_assembly",
+        "split_stages": ("compute_pressure_inputs", "assemble_pEqn"),
+        "run_one_stages": ("compute_pressure_inputs", "assemble_pEqn"),
+        "split_outputs": {
+            "compute_pressure_inputs": ("HbyA", "phiHbyA", "rAU"),
+            "assemble_pEqn": ("pEqn",),
+        },
+    },
+    {
+        "name": "linear_solve_results",
+        "split_stages": ("solve_UEqn", "solve_pEqn"),
+        "run_one_stages": ("solve_UEqn", "solve_pEqn"),
+        "split_outputs": {
+            "solve_UEqn": ("performance", "field_after"),
+            "solve_pEqn": ("performance", "p", "phi"),
+        },
+    },
+    {
+        "name": "final_correction",
+        "split_stages": ("correct_velocity_pressure_flux",),
+        "run_one_stages": ("update_phi_from_pEqn_flux", "correct_velocity_pressure_flux"),
+        "split_outputs": {
+            "correct_velocity_pressure_flux": ("U", "p", "phi"),
+        },
+    },
+    {
+        "name": "turbulence_updates",
+        "split_stages": ("momentum_transport_predict", "momentum_transport_correct"),
+        "run_one_stages": ("momentum_transport_predict", "momentum_transport_correct"),
+        "split_outputs": {
+            "momentum_transport_predict": ("case_path", "solver_name"),
+            "momentum_transport_correct": ("U", "p", "phi", "nut", "k", "omega"),
+        },
+    },
+)
+
+GPU_STAGE_KERNEL_ENTRYPOINTS = {
+    "momentum_assembly": ["gpu_rans_momentum_assembly", "gpu_rans_momentum_diffusion_coefficients", "gpu_rans_momentum_wall_diffusion_coefficients", "gpu_rans_momentum_convection_coefficients", "gpu_rans_momentum_equation_relaxation"],
+    "pressure_assembly": ["gpu_rans_pressure_inputs", "gpu_rans_consistent_rAtU", "gpu_rans_momentum_hbyA_source", "gpu_rans_momentum_hbyA_face_accumulate", "gpu_rans_momentum_hbyA_finish", "gpu_rans_momentum_hbyA_fixed_value_boundary", "gpu_rans_surface_flux_from_cells", "gpu_rans_pressure_assembly", "gpu_rans_pressure_laplacian_coefficients", "gpu_rans_pressure_mixed_boundary_laplacian", "gpu_rans_pressure_source_from_flux", "gpu_rans_pressure_source_from_boundary_flux"],
+    "linear_solve_results": ["gpu_ldu_matvec_vector_asymmetric_diag", "gpu_ldu_matvec_vector_asymmetric_face_accumulate", "gpu_bicgstab_initialize_vector", "gpu_bicgstab_dot_vector", "gpu_bicgstab_update_direction_vector", "gpu_bicgstab_precondition_vector", "gpu_bicgstab_update_intermediate_vector_preconditioned", "gpu_bicgstab_update_solution_residual_vector_preconditioned", "gpu_vector_residual_squared", "gpu_ldu_matvec_scalar_symmetric_diag", "gpu_ldu_matvec_scalar_symmetric_face_accumulate", "gpu_pcg_initialize_scalar", "gpu_cg_dot_scalar", "gpu_pcg_update_solution_residual_scalar", "gpu_pcg_update_direction_scalar"],
+    "final_correction": ["gpu_rans_pressure_flux_correction", "gpu_rans_final_correction", "gpu_rans_pressure_velocity_correction"],
+    "turbulence_updates": ["gpu_rans_turbulence_update", "gpu_rans_omega_wall_update"],
+}
--- a/python/src/foam_stepper/gpu/kernels.py
+++ b/python/src/foam_stepper/gpu/kernels.py
@ -0,0 +1,445 @@
+"""Quadrants CUDA kernels for the GPU RANS stage graph."""
+
+from __future__ import annotations
+
+import quadrants as qd
+
+from .constants import GPU_STAGE_KERNEL_ENTRYPOINTS
+
+@qd.kernel
+def gpu_rans_momentum_assembly(
+    n_cells: int,
+    u_internal: qd.types.NDArray[qd.f64, 2],
+    diag: qd.types.NDArray[qd.f64, 1],
+    source: qd.types.NDArray[qd.f64, 2],
+) -> None:
+    for cell in range(n_cells):
+        diag[cell] = 1.0
+        source[cell, 0] = u_internal[cell, 0]
+        source[cell, 1] = u_internal[cell, 1]
+        source[cell, 2] = u_internal[cell, 2]
+
+@qd.kernel
+def gpu_rans_momentum_diffusion_coefficients(
+    n_internal_faces: int,
+    owner: qd.types.NDArray[qd.i32, 1],
+    neighbour: qd.types.NDArray[qd.i32, 1],
+    nut_internal: qd.types.NDArray[qd.f64, 1],
+    laminar_nu: float,
+    cell_centres: qd.types.NDArray[qd.f64, 2],
+    sf: qd.types.NDArray[qd.f64, 2],
+    mag_sf: qd.types.NDArray[qd.f64, 1],
+    diag: qd.types.NDArray[qd.f64, 1],
+    upper: qd.types.NDArray[qd.f64, 1],
+    lower: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for face in range(n_internal_faces):
+        owner_cell = owner[face]
+        neighbour_cell = neighbour[face]
+        dx0 = cell_centres[neighbour_cell, 0] - cell_centres[owner_cell, 0]
+        dx1 = cell_centres[neighbour_cell, 1] - cell_centres[owner_cell, 1]
+        dx2 = cell_centres[neighbour_cell, 2] - cell_centres[owner_cell, 2]
+        projected_delta = dx0 * sf[face, 0] + dx1 * sf[face, 1] + dx2 * sf[face, 2]
+        if projected_delta < 0.0:
+            projected_delta = -projected_delta
+        if projected_delta < 1.0e-300:
+            projected_delta = 1.0e-300
+        effective_nu = laminar_nu + 0.5 * (nut_internal[owner_cell] + nut_internal[neighbour_cell])
+        coeff = effective_nu * mag_sf[face] * mag_sf[face] / projected_delta
+        upper[face] = -coeff
+        lower[face] = -coeff
+        qd.atomic_add(diag[owner_cell], coeff)
+        qd.atomic_add(diag[neighbour_cell], coeff)
+
+
+
+@qd.kernel
+def gpu_rans_momentum_wall_diffusion_coefficients(
+    n_boundary_faces: int,
+    face_cells: qd.types.NDArray[qd.i32, 1],
+    boundary_values: qd.types.NDArray[qd.f64, 2],
+    nut_boundary: qd.types.NDArray[qd.f64, 1],
+    laminar_nu: float,
+    cell_centres: qd.types.NDArray[qd.f64, 2],
+    face_centres: qd.types.NDArray[qd.f64, 2],
+    face_area_vectors: qd.types.NDArray[qd.f64, 2],
+    face_area_magnitudes: qd.types.NDArray[qd.f64, 1],
+    diag: qd.types.NDArray[qd.f64, 1],
+    source: qd.types.NDArray[qd.f64, 2],
+) -> None:
+    for boundary_face in range(n_boundary_faces):
+        cell = face_cells[boundary_face]
+        dx0 = face_centres[boundary_face, 0] - cell_centres[cell, 0]
+        dx1 = face_centres[boundary_face, 1] - cell_centres[cell, 1]
+        dx2 = face_centres[boundary_face, 2] - cell_centres[cell, 2]
+        projected_delta = dx0 * face_area_vectors[boundary_face, 0] + dx1 * face_area_vectors[boundary_face, 1] + dx2 * face_area_vectors[boundary_face, 2]
+        if projected_delta < 0.0:
+            projected_delta = -projected_delta
+        if projected_delta < 1.0e-300:
+            projected_delta = 1.0e-300
+        mag_sf = face_area_magnitudes[boundary_face]
+        coeff = (laminar_nu + nut_boundary[boundary_face]) * mag_sf * mag_sf / projected_delta
+        qd.atomic_add(diag[cell], coeff)
+        qd.atomic_add(source[cell, 0], coeff * boundary_values[boundary_face, 0])
+        qd.atomic_add(source[cell, 1], coeff * boundary_values[boundary_face, 1])
+        qd.atomic_add(source[cell, 2], coeff * boundary_values[boundary_face, 2])
+
+@qd.kernel
+def gpu_rans_momentum_convection_coefficients(
+    n_internal_faces: int,
+    owner: qd.types.NDArray[qd.i32, 1],
+    neighbour: qd.types.NDArray[qd.i32, 1],
+    phi: qd.types.NDArray[qd.f64, 1],
+    diag: qd.types.NDArray[qd.f64, 1],
+    upper: qd.types.NDArray[qd.f64, 1],
+    lower: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for face in range(n_internal_faces):
+        owner_cell = owner[face]
+        neighbour_cell = neighbour[face]
+        flux = phi[face]
+        if flux >= 0.0:
+            qd.atomic_add(diag[owner_cell], flux)
+            qd.atomic_add(lower[face], -flux)
+        else:
+            qd.atomic_add(diag[neighbour_cell], -flux)
+            qd.atomic_add(upper[face], flux)
+
+
+
+@qd.kernel
+def gpu_rans_momentum_equation_relaxation(
+    n_cells: int,
+    u_internal: qd.types.NDArray[qd.f64, 2],
+    alpha: float,
+    diag: qd.types.NDArray[qd.f64, 1],
+    source: qd.types.NDArray[qd.f64, 2],
+) -> None:
+    for cell in range(n_cells):
+        old_diag = diag[cell]
+        relaxed_diag = old_diag / alpha
+        source_scale = relaxed_diag - old_diag
+        diag[cell] = relaxed_diag
+        source[cell, 0] += source_scale * u_internal[cell, 0]
+        source[cell, 1] += source_scale * u_internal[cell, 1]
+        source[cell, 2] += source_scale * u_internal[cell, 2]
+
+
+@qd.kernel
+def gpu_rans_momentum_hbyA_source(
+    n_cells: int,
+    source: qd.types.NDArray[qd.f64, 2],
+    HbyA: qd.types.NDArray[qd.f64, 2],
+) -> None:
+    for cell in range(n_cells):
+        HbyA[cell, 0] = source[cell, 0]
+        HbyA[cell, 1] = source[cell, 1]
+        HbyA[cell, 2] = source[cell, 2]
+
+
+@qd.kernel
+def gpu_rans_momentum_hbyA_face_accumulate(
+    n_internal_faces: int,
+    owner: qd.types.NDArray[qd.i32, 1],
+    neighbour: qd.types.NDArray[qd.i32, 1],
+    upper: qd.types.NDArray[qd.f64, 1],
+    lower: qd.types.NDArray[qd.f64, 1],
+    u_internal: qd.types.NDArray[qd.f64, 2],
+    HbyA: qd.types.NDArray[qd.f64, 2],
+) -> None:
+    for face in range(n_internal_faces):
+        owner_cell = owner[face]
+        neighbour_cell = neighbour[face]
+        upper_coeff = upper[face]
+        lower_coeff = lower[face]
+        qd.atomic_add(HbyA[owner_cell, 0], -upper_coeff * u_internal[neighbour_cell, 0])
+        qd.atomic_add(HbyA[owner_cell, 1], -upper_coeff * u_internal[neighbour_cell, 1])
+        qd.atomic_add(HbyA[owner_cell, 2], -upper_coeff * u_internal[neighbour_cell, 2])
+        qd.atomic_add(HbyA[neighbour_cell, 0], -lower_coeff * u_internal[owner_cell, 0])
+        qd.atomic_add(HbyA[neighbour_cell, 1], -lower_coeff * u_internal[owner_cell, 1])
+        qd.atomic_add(HbyA[neighbour_cell, 2], -lower_coeff * u_internal[owner_cell, 2])
+
+
+@qd.kernel
+def gpu_rans_momentum_hbyA_finish(
+    n_cells: int,
+    diag: qd.types.NDArray[qd.f64, 1],
+    HbyA: qd.types.NDArray[qd.f64, 2],
+) -> None:
+    for cell in range(n_cells):
+        inv_diag = 1.0 / diag[cell]
+        HbyA[cell, 0] *= inv_diag
+        HbyA[cell, 1] *= inv_diag
+        HbyA[cell, 2] *= inv_diag
+
+
+@qd.kernel
+def gpu_rans_momentum_hbyA_fixed_value_boundary(
+    n_boundary_faces: int,
+    face_cells: qd.types.NDArray[qd.i32, 1],
+    boundary_values: qd.types.NDArray[qd.f64, 2],
+    HbyA: qd.types.NDArray[qd.f64, 2],
+) -> None:
+    for boundary_face in range(n_boundary_faces):
+        cell = face_cells[boundary_face]
+        HbyA[cell, 0] = boundary_values[boundary_face, 0]
+        HbyA[cell, 1] = boundary_values[boundary_face, 1]
+        HbyA[cell, 2] = boundary_values[boundary_face, 2]
+
+
+@qd.kernel
+def gpu_rans_pressure_inputs(
+    n_cells: int,
+    u_diag: qd.types.NDArray[qd.f64, 1],
+    cell_volumes: qd.types.NDArray[qd.f64, 1],
+    rAU: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for cell in range(n_cells):
+        rAU[cell] = cell_volumes[cell] / u_diag[cell]
+
+
+@qd.kernel
+def gpu_rans_consistent_rAtU(
+    n_cells: int,
+    rAU: qd.types.NDArray[qd.f64, 1],
+    ratu_factor: float,
+    rAtU: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for cell in range(n_cells):
+        rAtU[cell] = ratu_factor * rAU[cell]
+
+
+@qd.kernel
+def gpu_rans_surface_flux_from_cells(
+    n_internal_faces: int,
+    owner: qd.types.NDArray[qd.i32, 1],
+    neighbour: qd.types.NDArray[qd.i32, 1],
+    cell_vector: qd.types.NDArray[qd.f64, 2],
+    sf: qd.types.NDArray[qd.f64, 2],
+    out: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for face in range(n_internal_faces):
+        owner_cell = owner[face]
+        neighbour_cell = neighbour[face]
+        out[face] = 0.5 * (
+            (cell_vector[owner_cell, 0] + cell_vector[neighbour_cell, 0]) * sf[face, 0]
+            + (cell_vector[owner_cell, 1] + cell_vector[neighbour_cell, 1]) * sf[face, 1]
+            + (cell_vector[owner_cell, 2] + cell_vector[neighbour_cell, 2]) * sf[face, 2]
+        )
+
+
+@qd.kernel
+def gpu_rans_pressure_assembly(
+    n_cells: int,
+    p_internal: qd.types.NDArray[qd.f64, 1],
+    diag: qd.types.NDArray[qd.f64, 1],
+    source: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for cell in range(n_cells):
+        diag[cell] = 0.0
+        source[cell] = 0.0
+
+
+@qd.kernel
+def gpu_rans_pressure_laplacian_coefficients(
+    n_internal_faces: int,
+    owner: qd.types.NDArray[qd.i32, 1],
+    neighbour: qd.types.NDArray[qd.i32, 1],
+    rAtU: qd.types.NDArray[qd.f64, 1],
+    cell_centres: qd.types.NDArray[qd.f64, 2],
+    sf: qd.types.NDArray[qd.f64, 2],
+    mag_sf: qd.types.NDArray[qd.f64, 1],
+    diag: qd.types.NDArray[qd.f64, 1],
+    upper: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for face in range(n_internal_faces):
+        owner_cell = owner[face]
+        neighbour_cell = neighbour[face]
+        dx0 = cell_centres[neighbour_cell, 0] - cell_centres[owner_cell, 0]
+        dx1 = cell_centres[neighbour_cell, 1] - cell_centres[owner_cell, 1]
+        dx2 = cell_centres[neighbour_cell, 2] - cell_centres[owner_cell, 2]
+        projected_delta = dx0 * sf[face, 0] + dx1 * sf[face, 1] + dx2 * sf[face, 2]
+        if projected_delta < 0.0:
+            projected_delta = -projected_delta
+        if projected_delta < 1.0e-300:
+            projected_delta = 1.0e-300
+        coeff = 0.5 * (rAtU[owner_cell] + rAtU[neighbour_cell]) * mag_sf[face] * mag_sf[face] / projected_delta
+        upper[face] = coeff
+        qd.atomic_add(diag[owner_cell], -coeff)
+        qd.atomic_add(diag[neighbour_cell], -coeff)
+
+
+@qd.kernel
+def gpu_rans_pressure_source_from_flux(
+    n_internal_faces: int,
+    owner: qd.types.NDArray[qd.i32, 1],
+    neighbour: qd.types.NDArray[qd.i32, 1],
+    phi: qd.types.NDArray[qd.f64, 1],
+    source: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for face in range(n_internal_faces):
+        flux = phi[face]
+        qd.atomic_add(source[owner[face]], flux)
+        qd.atomic_add(source[neighbour[face]], -flux)
+
+
+@qd.kernel
+def gpu_rans_pressure_source_from_boundary_flux(
+    n_boundary_faces: int,
+    face_cells: qd.types.NDArray[qd.i32, 1],
+    phi_boundary: qd.types.NDArray[qd.f64, 1],
+    source: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for face in range(n_boundary_faces):
+        qd.atomic_add(source[face_cells[face]], phi_boundary[face])
+
+
+@qd.kernel
+def gpu_rans_pressure_mixed_boundary_laplacian(
+    n_boundary_faces: int,
+    face_cells: qd.types.NDArray[qd.i32, 1],
+    boundary_scales: qd.types.NDArray[qd.f64, 1],
+    boundary_values: qd.types.NDArray[qd.f64, 1],
+    rAtU: qd.types.NDArray[qd.f64, 1],
+    diag: qd.types.NDArray[qd.f64, 1],
+    source: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for face in range(n_boundary_faces):
+        cell = face_cells[face]
+        coeff = rAtU[cell] * boundary_scales[face]
+        qd.atomic_add(diag[cell], -coeff)
+        qd.atomic_add(source[cell], -coeff * boundary_values[face])
+
+
+@qd.kernel
+def gpu_rans_face_flux_copy(
+    n_internal_faces: int,
+    source: qd.types.NDArray[qd.f64, 1],
+    out: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for face in range(n_internal_faces):
+        out[face] = source[face]
+
+
+@qd.kernel
+def gpu_rans_pressure_flux_correction(
+    n_internal_faces: int,
+    owner: qd.types.NDArray[qd.i32, 1],
+    neighbour: qd.types.NDArray[qd.i32, 1],
+    p_solved: qd.types.NDArray[qd.f64, 1],
+    upper: qd.types.NDArray[qd.f64, 1],
+    phiHbyA: qd.types.NDArray[qd.f64, 1],
+    phi_out: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for face in range(n_internal_faces):
+        pressure_jump = p_solved[neighbour[face]] - p_solved[owner[face]]
+        phi_out[face] = phiHbyA[face] - upper[face] * pressure_jump
+
+
+@qd.kernel
+def gpu_rans_final_correction(
+    n_cells: int,
+    HbyA: qd.types.NDArray[qd.f64, 2],
+    p_solved: qd.types.NDArray[qd.f64, 1],
+    p_internal: qd.types.NDArray[qd.f64, 1],
+    u_out: qd.types.NDArray[qd.f64, 2],
+    p_out: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for cell in range(n_cells):
+        u_out[cell, 0] = HbyA[cell, 0]
+        u_out[cell, 1] = HbyA[cell, 1]
+        u_out[cell, 2] = HbyA[cell, 2]
+        p_out[cell] = p_internal[cell] + p_solved[cell]
+
+
+@qd.kernel
+def gpu_rans_pressure_velocity_correction(
+    n_internal_faces: int,
+    owner: qd.types.NDArray[qd.i32, 1],
+    neighbour: qd.types.NDArray[qd.i32, 1],
+    rAU: qd.types.NDArray[qd.f64, 1],
+    p_solved: qd.types.NDArray[qd.f64, 1],
+    cell_centres: qd.types.NDArray[qd.f64, 2],
+    sf: qd.types.NDArray[qd.f64, 2],
+    u_out: qd.types.NDArray[qd.f64, 2],
+) -> None:
+    for face in range(n_internal_faces):
+        owner_cell = owner[face]
+        neighbour_cell = neighbour[face]
+        dx0 = cell_centres[neighbour_cell, 0] - cell_centres[owner_cell, 0]
+        dx1 = cell_centres[neighbour_cell, 1] - cell_centres[owner_cell, 1]
+        dx2 = cell_centres[neighbour_cell, 2] - cell_centres[owner_cell, 2]
+        projected_delta = dx0 * sf[face, 0] + dx1 * sf[face, 1] + dx2 * sf[face, 2]
+        if projected_delta < 0.0:
+            projected_delta = -projected_delta
+        if projected_delta < 1.0e-300:
+            projected_delta = 1.0e-300
+        gradient_scale = (p_solved[neighbour_cell] - p_solved[owner_cell]) / projected_delta
+        owner_scale = rAU[owner_cell] * gradient_scale
+        neighbour_scale = rAU[neighbour_cell] * gradient_scale
+        u_out[owner_cell, 0] -= owner_scale * sf[face, 0]
+        u_out[owner_cell, 1] -= owner_scale * sf[face, 1]
+        u_out[owner_cell, 2] -= owner_scale * sf[face, 2]
+        u_out[neighbour_cell, 0] -= neighbour_scale * sf[face, 0]
+        u_out[neighbour_cell, 1] -= neighbour_scale * sf[face, 1]
+        u_out[neighbour_cell, 2] -= neighbour_scale * sf[face, 2]
+
+
+@qd.kernel
+def gpu_rans_turbulence_update(
+    n_cells: int,
+    nut_internal: qd.types.NDArray[qd.f64, 1],
+    k_internal: qd.types.NDArray[qd.f64, 1],
+    omega_internal: qd.types.NDArray[qd.f64, 1],
+    nut_out: qd.types.NDArray[qd.f64, 1],
+    k_out: qd.types.NDArray[qd.f64, 1],
+    omega_out: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for cell in range(n_cells):
+        nut_out[cell] = nut_internal[cell]
+        k_out[cell] = k_internal[cell]
+        omega_out[cell] = omega_internal[cell]
+
+
+@qd.kernel
+def gpu_rans_omega_wall_update(
+    n_wall_cells: int,
+    wall_cells: qd.types.NDArray[qd.i32, 1],
+    wall_distances: qd.types.NDArray[qd.f64, 1],
+    laminar_nu: float,
+    beta1: float,
+    omega_out: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for index in range(n_wall_cells):
+        distance = wall_distances[index]
+        if distance < 1.0e-300:
+            distance = 1.0e-300
+        omega_wall = 6.0 * laminar_nu / (beta1 * distance * distance)
+        cell = wall_cells[index]
+        if omega_wall > omega_out[cell]:
+            omega_out[cell] = omega_wall
+
+__all__ = [
+    "GPU_STAGE_KERNEL_ENTRYPOINTS",
+    "gpu_rans_momentum_assembly",
+    "gpu_rans_momentum_diffusion_coefficients",
+    "gpu_rans_momentum_convection_coefficients",
+    "gpu_rans_momentum_equation_relaxation",
+    "gpu_rans_momentum_hbyA_source",
+    "gpu_rans_momentum_hbyA_face_accumulate",
+    "gpu_rans_momentum_hbyA_finish",
+    "gpu_rans_pressure_inputs",
+    "gpu_rans_consistent_rAtU",
+    "gpu_rans_pressure_assembly",
+    "gpu_rans_pressure_laplacian_coefficients",
+    "gpu_rans_pressure_source_from_flux",
+    "gpu_rans_pressure_source_from_boundary_flux",
+    "gpu_rans_pressure_mixed_boundary_laplacian",
+    "gpu_rans_face_flux_copy",
+    "gpu_rans_surface_flux_from_cells",
+    "gpu_rans_pressure_flux_correction",
+    "gpu_rans_pressure_velocity_correction",
+    "gpu_rans_final_correction",
+    "gpu_rans_turbulence_update",
+    "gpu_rans_omega_wall_update",
+]
--- a/python/src/foam_stepper/gpu/linear_solve.py
+++ b/python/src/foam_stepper/gpu/linear_solve.py
--- a/python/src/foam_stepper/state.py
+++ b/python/src/foam_stepper/state.py
@ -0,0 +1,453 @@
+"""Explicit solver-state exports and validation diagnostics."""
+
+from __future__ import annotations
+
+from collections.abc import Mapping
+from typing import Any, Iterable
+
+import numpy as np
+
+
+DEFAULT_TURBULENCE_FIELDS = ("nut", "k", "omega")
+
+
+class SolverStateValidationError(ValueError):
+    """Localized validation failure for exported solver-state data."""
+
+    def __init__(self, path: str, message: str, *, expected: Any = None, actual: Any = None) -> None:
+        super().__init__(f"{path}: {message}")
+        self.path = path
+        self.message = message
+        self.expected = expected
+        self.actual = actual
+
+    def to_dict(self) -> dict[str, Any]:
+        out: dict[str, Any] = {"path": self.path, "message": self.message}
+        if self.expected is not None:
+            out["expected"] = _json_ready(self.expected)
+        if self.actual is not None:
+            out["actual"] = _json_ready(self.actual)
+        return out
+
+
+def export_solver_state(
+    mesh: Any,
+    fields: Mapping[str, Any] | Any,
+    *,
+    required_fields: Iterable[str] = (),
+    turbulence_fields: Iterable[str] = DEFAULT_TURBULENCE_FIELDS,
+) -> dict[str, Any]:
+    """Return mesh, field, boundary, and turbulence arrays as plain Python/NumPy data."""
+
+    field_map = _field_mapping(fields)
+    state = {
+        "schema_version": 1,
+        "mesh": _export_mesh(mesh),
+        "fields": {name: _export_field(field) for name, field in field_map.items()},
+        "turbulence": {name: _export_field(field_map[name]) for name in turbulence_fields if name in field_map},
+    }
+    validate_solver_state(state, required_fields=required_fields)
+    return state
+
+
+def export_matrix_state(matrix: Any, mesh_state: Mapping[str, Any] | None = None) -> dict[str, Any]:
+    """Return an assembled matrix and solver intermediates as explicit arrays."""
+
+    state = {
+        "name": getattr(matrix, "name", ""),
+        "field_name": getattr(matrix, "field_name", ""),
+        "value_rank": getattr(matrix, "value_rank", ""),
+        "dimensions": getattr(matrix, "dimensions", ""),
+        "flags": {
+            "has_diag": bool(getattr(matrix, "has_diag", False)),
+            "has_upper": bool(getattr(matrix, "has_upper", False)),
+            "has_lower": bool(getattr(matrix, "has_lower", False)),
+            "diagonal": bool(getattr(matrix, "diagonal", False)),
+            "symmetric": bool(getattr(matrix, "symmetric", False)),
+            "asymmetric": bool(getattr(matrix, "asymmetric", False)),
+        },
+        "diag": _array(getattr(matrix, "diag")),
+        "upper": _optional_array(getattr(matrix, "upper", None)),
+        "lower": _optional_array(getattr(matrix, "lower", None)),
+        "source": _array(getattr(matrix, "source")),
+        "psi": _export_field(getattr(matrix, "psi")),
+        "internal_coeffs": [_array(item) for item in getattr(matrix, "internal_coeffs", [])],
+        "boundary_coeffs": [_array(item) for item in getattr(matrix, "boundary_coeffs", [])],
+        "derived": {},
+    }
+    for name in ("A", "H", "H1", "flux", "face_flux_correction"):
+        value = getattr(matrix, name, None)
+        if callable(value):
+            value = value()
+        if value is not None:
+            state["derived"][name] = _export_field(value)
+    for name in ("residual", "D", "DD"):
+        value = getattr(matrix, name, None)
+        if callable(value):
+            value = value()
+        if value is not None:
+            state["derived"][name] = _array(value)
+    validate_matrix_state(state, mesh_state=mesh_state)
+    return state
+
+
+def validate_solver_state(
+    state: Mapping[str, Any],
+    *,
+    required_fields: Iterable[str] = (),
+    turbulence_fields: Iterable[str] = (),
+) -> None:
+    mesh = _mapping(state.get("mesh"), "mesh")
+    sizes = _mapping(mesh.get("sizes"), "mesh.sizes")
+    n_points = _positive_int(sizes.get("n_points"), "mesh.sizes.n_points", allow_zero=False)
+    n_faces = _positive_int(sizes.get("n_faces"), "mesh.sizes.n_faces", allow_zero=False)
+    n_internal_faces = _positive_int(sizes.get("n_internal_faces"), "mesh.sizes.n_internal_faces", allow_zero=True)
+    n_cells = _positive_int(sizes.get("n_cells"), "mesh.sizes.n_cells", allow_zero=False)
+    if n_internal_faces > n_faces:
+        raise SolverStateValidationError(
+            "mesh.sizes.n_internal_faces",
+            "cannot exceed n_faces",
+            expected={"max": n_faces},
+            actual=n_internal_faces,
+        )
+
+    connectivity = _mapping(mesh.get("connectivity"), "mesh.connectivity")
+    geometry = _mapping(mesh.get("geometry"), "mesh.geometry")
+    _require_array(geometry.get("points"), "mesh.geometry.points", shape=(n_points, 3), numeric=True)
+    _require_array(geometry.get("V"), "mesh.geometry.V", shape=(n_cells,), numeric=True)
+    _require_array(geometry.get("C"), "mesh.geometry.C", shape=(n_cells, 3), numeric=True)
+    _require_array(geometry.get("Cf"), "mesh.geometry.Cf", shape=(n_internal_faces, 3), numeric=True)
+    _require_array(geometry.get("Sf"), "mesh.geometry.Sf", shape=(n_internal_faces, 3), numeric=True)
+    _require_array(geometry.get("magSf"), "mesh.geometry.magSf", shape=(n_internal_faces,), numeric=True)
+
+    _validate_ragged(connectivity.get("faces"), "mesh.connectivity.faces", rows=n_faces)
+    _validate_ragged(connectivity.get("cells"), "mesh.connectivity.cells", rows=n_cells)
+    _require_array(connectivity.get("owner"), "mesh.connectivity.owner", shape=(n_internal_faces,), integer=True)
+    _require_array(connectivity.get("neighbour"), "mesh.connectivity.neighbour", shape=(n_internal_faces,), integer=True)
+
+    ldu = _mapping(connectivity.get("ldu"), "mesh.connectivity.ldu")
+    _require_array(ldu.get("lower_addr"), "mesh.connectivity.ldu.lower_addr", shape=(n_internal_faces,), integer=True)
+    _require_array(ldu.get("upper_addr"), "mesh.connectivity.ldu.upper_addr", shape=(n_internal_faces,), integer=True)
+
+    patches = _patches(mesh.get("patches"), n_faces=n_faces, n_cells=n_cells)
+    fields = _mapping(state.get("fields"), "fields")
+    for name in required_fields:
+        if name not in fields:
+            raise SolverStateValidationError("fields", "missing required field", expected=name, actual=sorted(fields))
+    for name in turbulence_fields:
+        if name not in fields:
+            raise SolverStateValidationError("turbulence", "missing turbulence field", expected=name, actual=sorted(fields))
+
+    for name, field in fields.items():
+        _validate_field(_mapping(field, f"fields.{name}"), f"fields.{name}", patches=patches, n_cells=n_cells, n_internal_faces=n_internal_faces)
+
+
+def validate_matrix_state(state: Mapping[str, Any], *, mesh_state: Mapping[str, Any] | None = None) -> None:
+    path = f"matrices.{state.get('field_name', '<unknown>')}"
+    rank = _rank_from_value_rank(str(state.get("value_rank", "")), path)
+    n_cells = None
+    n_internal_faces = None
+    patches: list[Mapping[str, Any]] = []
+    if mesh_state is not None:
+        mesh = _mapping(mesh_state.get("mesh", mesh_state), "mesh")
+        sizes = _mapping(mesh.get("sizes"), "mesh.sizes")
+        n_cells = _positive_int(sizes.get("n_cells"), "mesh.sizes.n_cells", allow_zero=False)
+        n_internal_faces = _positive_int(sizes.get("n_internal_faces"), "mesh.sizes.n_internal_faces", allow_zero=True)
+        patches = list(_mapping(mesh, "mesh").get("patches", []))
+
+    if n_cells is not None:
+        _require_array(state.get("diag"), f"{path}.diag", shape=(n_cells,), numeric=True)
+        _require_array(state.get("source"), f"{path}.source", shape=_ranked_shape(n_cells, rank), numeric=True)
+    else:
+        _require_array(state.get("diag"), f"{path}.diag", numeric=True)
+        _require_array(state.get("source"), f"{path}.source", numeric=True)
+
+    if n_internal_faces is not None:
+        _optional_valid_array(state.get("upper"), f"{path}.upper", shape=(n_internal_faces,), numeric=True)
+        _optional_valid_array(state.get("lower"), f"{path}.lower", shape=(n_internal_faces,), numeric=True)
+    else:
+        _optional_valid_array(state.get("upper"), f"{path}.upper", numeric=True)
+        _optional_valid_array(state.get("lower"), f"{path}.lower", numeric=True)
+
+    if patches:
+        _validate_coeff_list(state.get("internal_coeffs", []), f"{path}.internal_coeffs", patches=patches, rank=rank)
+        _validate_coeff_list(state.get("boundary_coeffs", []), f"{path}.boundary_coeffs", patches=patches, rank=rank)
+
+
+def describe_solver_state(state: Mapping[str, Any]) -> dict[str, Any]:
+    return _describe(state)
+
+
+def describe_matrix_state(state: Mapping[str, Any]) -> dict[str, Any]:
+    return _describe(state)
+
+
+def _export_mesh(mesh: Any) -> dict[str, Any]:
+    return {
+        "sizes": {
+            "n_points": int(getattr(mesh, "n_points")),
+            "n_faces": int(getattr(mesh, "n_faces")),
+            "n_internal_faces": int(getattr(mesh, "n_internal_faces")),
+            "n_cells": int(getattr(mesh, "n_cells")),
+        },
+        "connectivity": {
+            "faces": _export_ragged(getattr(mesh, "faces")),
+            "cells": _export_ragged(getattr(mesh, "cells")),
+            "owner": _array(getattr(mesh, "owner")),
+            "neighbour": _array(getattr(mesh, "neighbour")),
+            "ldu": _export_ldu(getattr(mesh, "ldu", {})),
+        },
+        "geometry": {
+            "points": _array(getattr(mesh, "points")),
+            "V": _array(getattr(mesh, "V")),
+            "C": _array(getattr(mesh, "C")),
+            "Cf": _array(getattr(mesh, "Cf")),
+            "Sf": _array(getattr(mesh, "Sf")),
+            "magSf": _array(getattr(mesh, "magSf")),
+        },
+        "patches": [_export_patch(patch) for patch in getattr(mesh, "boundary")],
+    }
+
+
+def _export_ragged(value: Any) -> dict[str, np.ndarray]:
+    return {"offsets": _array(getattr(value, "offsets")), "values": _array(getattr(value, "values"))}
+
+
+def _export_ldu(ldu: Mapping[str, Any]) -> dict[str, Any]:
+    return {
+        "lower_addr": _array(ldu.get("lower_addr")),
+        "upper_addr": _array(ldu.get("upper_addr")),
+        "patch_addr": [
+            {"patch_index": int(entry.get("patch_index")), "addr": _array(entry.get("addr"))}
+            for entry in ldu.get("patch_addr", [])
+        ],
+    }
+
+
+def _export_patch(patch: Any) -> dict[str, Any]:
+    return {
+        "name": getattr(patch, "name"),
+        "type": getattr(patch, "type"),
+        "index": int(getattr(patch, "index")),
+        "start": int(getattr(patch, "start")),
+        "size": int(getattr(patch, "size")),
+        "coupled": bool(getattr(patch, "coupled")),
+        "constraint": bool(getattr(patch, "constraint")),
+        "face_cells": _array(getattr(patch, "face_cells")),
+        "face_indices": _array(getattr(patch, "face_indices")),
+        "Cf": _array(getattr(patch, "Cf")),
+        "Sf": _array(getattr(patch, "Sf")),
+        "magSf": _array(getattr(patch, "magSf")),
+    }
+
+
+def _export_field(field: Any) -> dict[str, Any]:
+    return {
+        "name": getattr(field, "name"),
+        "kind": getattr(field, "kind"),
+        "rank": _rank_from_kind(getattr(field, "kind", "")),
+        "dimensions": getattr(field, "dimensions"),
+        "entity_kind": getattr(field, "entity_kind"),
+        "entity_count": int(getattr(field, "entity_count")),
+        "internal": _array(getattr(field, "internal")),
+        "boundary": {name: _export_patch_field(patch) for name, patch in getattr(field, "boundary").items()},
+    }
+
+
+def _export_patch_field(patch: Any) -> dict[str, Any]:
+    return {
+        "name": getattr(patch, "name"),
+        "type": getattr(patch, "type"),
+        "values": _array(getattr(patch, "values")),
+        "fixes_value": bool(getattr(patch, "fixes_value")),
+        "assignable": bool(getattr(patch, "assignable")),
+        "coupled": bool(getattr(patch, "coupled")),
+        "updated": bool(getattr(patch, "updated")),
+        "patch_internal": _optional_array(getattr(patch, "patch_internal", None)),
+        "value_internal_coeffs": _optional_array(getattr(patch, "value_internal_coeffs", None)),
+        "value_boundary_coeffs": _optional_array(getattr(patch, "value_boundary_coeffs", None)),
+        "gradient_internal_coeffs": _optional_array(getattr(patch, "gradient_internal_coeffs", None)),
+        "gradient_boundary_coeffs": _optional_array(getattr(patch, "gradient_boundary_coeffs", None)),
+    }
+
+
+def _field_mapping(fields: Mapping[str, Any] | Any) -> Mapping[str, Any]:
+    if isinstance(fields, Mapping):
+        return fields
+    if hasattr(fields, "fields") and isinstance(fields.fields, Mapping):
+        return fields.fields
+    raise SolverStateValidationError("fields", "expected field mapping", actual=type(fields).__name__)
+
+
+def _array(value: Any) -> np.ndarray:
+    if value is None:
+        raise SolverStateValidationError("array", "missing required array")
+    return np.asarray(value)
+
+
+def _optional_array(value: Any) -> np.ndarray | None:
+    return None if value is None else np.asarray(value)
+
+
+def _optional_valid_array(value: Any, path: str, *, shape: tuple[int, ...] | None = None, numeric: bool = False) -> None:
+    if value is None:
+        return
+    _require_array(value, path, shape=shape, numeric=numeric)
+
+
+def _require_array(
+    value: Any,
+    path: str,
+    *,
+    shape: tuple[int, ...] | None = None,
+    numeric: bool = False,
+    integer: bool = False,
+) -> np.ndarray:
+    if value is None:
+        raise SolverStateValidationError(path, "missing required array")
+    arr = np.asarray(value)
+    if shape is not None and tuple(arr.shape) != shape:
+        raise SolverStateValidationError(path, "shape mismatch", expected=list(shape), actual=list(arr.shape))
+    if numeric and not np.issubdtype(arr.dtype, np.number):
+        raise SolverStateValidationError(path, "expected numeric dtype", actual=str(arr.dtype))
+    if integer and not np.issubdtype(arr.dtype, np.integer):
+        raise SolverStateValidationError(path, "expected integer dtype", actual=str(arr.dtype))
+    return arr
+
+
+def _validate_ragged(value: Any, path: str, *, rows: int) -> None:
+    mapping = _mapping(value, path)
+    offsets = _require_array(mapping.get("offsets"), f"{path}.offsets", shape=(rows + 1,), integer=True)
+    _require_array(mapping.get("values"), f"{path}.values", integer=True)
+    if offsets.size and int(offsets[0]) != 0:
+        raise SolverStateValidationError(f"{path}.offsets", "first offset must be zero", expected=0, actual=int(offsets[0]))
+    if offsets.size > 1 and bool(np.any(offsets[1:] < offsets[:-1])):
+        raise SolverStateValidationError(f"{path}.offsets", "offsets must be monotonically nondecreasing")
+
+
+def _patches(value: Any, *, n_faces: int, n_cells: int) -> list[Mapping[str, Any]]:
+    if not isinstance(value, list):
+        raise SolverStateValidationError("mesh.patches", "expected patch list", actual=type(value).__name__)
+    patches = []
+    names: set[str] = set()
+    for index, patch_value in enumerate(value):
+        path = f"mesh.patches[{index}]"
+        patch = _mapping(patch_value, path)
+        name = str(patch.get("name", ""))
+        if not name:
+            raise SolverStateValidationError(f"{path}.name", "missing patch name")
+        if name in names:
+            raise SolverStateValidationError(f"{path}.name", "duplicate patch name", actual=name)
+        names.add(name)
+        size = _positive_int(patch.get("size"), f"{path}.size", allow_zero=True)
+        start = _positive_int(patch.get("start"), f"{path}.start", allow_zero=True)
+        if start + size > n_faces:
+            raise SolverStateValidationError(f"{path}.start", "patch faces exceed mesh face count", expected={"n_faces": n_faces}, actual={"start": start, "size": size})
+        _require_array(patch.get("face_cells"), f"{path}.face_cells", shape=(size,), integer=True)
+        _require_array(patch.get("face_indices"), f"{path}.face_indices", shape=(size,), integer=True)
+        _require_array(patch.get("Cf"), f"{path}.Cf", shape=(size, 3), numeric=True)
+        _require_array(patch.get("Sf"), f"{path}.Sf", shape=(size, 3), numeric=True)
+        _require_array(patch.get("magSf"), f"{path}.magSf", shape=(size,), numeric=True)
+        face_cells = np.asarray(patch.get("face_cells"))
+        if face_cells.size and (int(face_cells.min()) < 0 or int(face_cells.max()) >= n_cells):
+            raise SolverStateValidationError(f"{path}.face_cells", "face cell index out of range", expected={"min": 0, "max_exclusive": n_cells})
+        patches.append(patch)
+    return patches
+
+
+def _validate_field(field: Mapping[str, Any], path: str, *, patches: list[Mapping[str, Any]], n_cells: int, n_internal_faces: int) -> None:
+    rank = str(field.get("rank") or _rank_from_kind(str(field.get("kind", ""))))
+    entity_kind = str(field.get("entity_kind", ""))
+    entity_count = _positive_int(field.get("entity_count"), f"{path}.entity_count", allow_zero=True)
+    expected_count = n_internal_faces if entity_kind == "internal_face" else n_cells
+    if entity_count != expected_count:
+        raise SolverStateValidationError(f"{path}.entity_count", "topology entity count mismatch", expected=expected_count, actual=entity_count)
+    _require_array(field.get("internal"), f"{path}.internal", shape=_ranked_shape(expected_count, rank), numeric=True)
+    boundary = _mapping(field.get("boundary"), f"{path}.boundary")
+    patch_names = {str(patch["name"]) for patch in patches}
+    if set(boundary) != patch_names:
+        raise SolverStateValidationError(f"{path}.boundary", "patch set mismatch", expected=sorted(patch_names), actual=sorted(boundary))
+    patches_by_name = {str(patch["name"]): patch for patch in patches}
+    for patch_name, patch_field_value in boundary.items():
+        patch_field_path = f"{path}.boundary.{patch_name}"
+        patch_field = _mapping(patch_field_value, patch_field_path)
+        patch_size = int(patches_by_name[str(patch_name)]["size"])
+        expected_shape = _ranked_shape(patch_size, rank)
+        _require_array(patch_field.get("values"), f"{patch_field_path}.values", shape=expected_shape, numeric=True)
+        _optional_valid_array(patch_field.get("patch_internal"), f"{patch_field_path}.patch_internal", shape=expected_shape, numeric=True)
+        for coeff_name in ("value_internal_coeffs", "value_boundary_coeffs", "gradient_internal_coeffs", "gradient_boundary_coeffs"):
+            _optional_valid_array(patch_field.get(coeff_name), f"{patch_field_path}.{coeff_name}", numeric=True)
+
+
+def _validate_coeff_list(value: Any, path: str, *, patches: list[Mapping[str, Any]], rank: str) -> None:
+    if not isinstance(value, list):
+        raise SolverStateValidationError(path, "expected coefficient list", actual=type(value).__name__)
+    if len(value) != len(patches):
+        raise SolverStateValidationError(path, "patch coefficient count mismatch", expected=len(patches), actual=len(value))
+    for index, coeffs in enumerate(value):
+        patch_size = int(patches[index]["size"])
+        _require_array(coeffs, f"{path}[{index}]", shape=_ranked_shape(patch_size, rank), numeric=True)
+
+
+def _mapping(value: Any, path: str) -> Mapping[str, Any]:
+    if not isinstance(value, Mapping):
+        raise SolverStateValidationError(path, "expected mapping", actual=type(value).__name__)
+    return value
+
+
+def _positive_int(value: Any, path: str, *, allow_zero: bool) -> int:
+    try:
+        number = int(value)
+    except (TypeError, ValueError) as exc:
+        raise SolverStateValidationError(path, "expected integer", actual=value) from exc
+    if number < 0 or (number == 0 and not allow_zero):
+        raise SolverStateValidationError(path, "invalid count", actual=number)
+    return number
+
+
+def _rank_from_kind(kind: str) -> str:
+    if kind.endswith("Vector"):
+        return "vector"
+    if kind.endswith("Scalar"):
+        return "scalar"
+    return "scalar"
+
+
+def _rank_from_value_rank(value_rank: str, path: str) -> str:
+    if value_rank in {"scalar", "vector"}:
+        return value_rank
+    raise SolverStateValidationError(f"{path}.value_rank", "expected scalar or vector rank", actual=value_rank)
+
+
+def _ranked_shape(size: int, rank: str) -> tuple[int, ...]:
+    if rank == "scalar":
+        return (size,)
+    if rank == "vector":
+        return (size, 3)
+    raise SolverStateValidationError("rank", "expected scalar or vector rank", actual=rank)
+
+
+def _describe(value: Any) -> Any:
+    if isinstance(value, np.ndarray):
+        return {"shape": [int(dim) for dim in value.shape], "dtype": str(value.dtype), "size": int(value.size)}
+    if isinstance(value, Mapping):
+        return {str(key): _describe(item) for key, item in value.items()}
+    if isinstance(value, list):
+        return [_describe(item) for item in value]
+    if isinstance(value, tuple):
+        return [_describe(item) for item in value]
+    return _json_ready(value)
+
+
+def _json_ready(value: Any) -> Any:
+    if isinstance(value, np.generic):
+        return value.item()
+    if isinstance(value, np.ndarray):
+        return value.tolist()
+    if isinstance(value, Mapping):
+        return {str(key): _json_ready(item) for key, item in value.items()}
+    if isinstance(value, (list, tuple)):
+        return [_json_ready(item) for item in value]
+    if isinstance(value, (str, int, float, bool)) or value is None:
+        return value
+    return repr(value)
--- a/scripts/build_openfoam_airfrans_subset.sh
+++ b/scripts/build_openfoam_airfrans_subset.sh
@ -11,20 +11,59 @@ THIRDPARTY_REPO_URL="${THIRDPARTY_REPO_URL:-https://github.com/OpenFOAM/ThirdPar
 AIRFRANS_REPO_URL="${AIRFRANS_REPO_URL:-https://github.com/Extrality/AirfRANS.git}"
 JOBS="${JOBS:-$(nproc)}"

-clone_if_missing() {
-    local url="$1"
-    local dir="$2"
+validate_dependency_tree() {
+    local dir="$1"
+    shift

-    if [[ -d "${ROOT_DIR}/${dir}/.git" ]]; then
-        printf 'using existing %s\n' "${dir}"
-    else
-        git clone --depth 1 "${url}" "${ROOT_DIR}/${dir}"
+    local root="${ROOT_DIR}/${dir}"
+    if [[ ! -d "${root}" ]]; then
+        printf 'dependency %s exists but is not a directory: %s\n' "${dir}" "${root}" >&2
+        return 1
+    fi
+
+    local -a missing=()
+    local expected
+    for expected in "$@"; do
+        if [[ ! -e "${root}/${expected}" ]]; then
+            missing+=("${expected}")
+        fi
+    done
+
+    if ((${#missing[@]})); then
+        printf 'dependency %s exists but is incomplete; missing expected path(s):' "${dir}" >&2
+        printf ' %s' "${missing[@]}" >&2
+        printf '\nMove or remove %s, or replace it with a complete checkout before rerunning.\n' "${root}" >&2
+        return 1
    fi
 }

-clone_if_missing "${OPENFOAM_REPO_URL}" OpenFOAM-14
-clone_if_missing "${THIRDPARTY_REPO_URL}" ThirdParty-14
-clone_if_missing "${AIRFRANS_REPO_URL}" airfrans
+clone_if_missing() {
+    local url="$1"
+    local dir="$2"
+    shift 2
+
+    local root="${ROOT_DIR}/${dir}"
+    if [[ -e "${root}" || -L "${root}" ]]; then
+        validate_dependency_tree "${dir}" "$@"
+        printf 'using existing %s\n' "${dir}"
+    else
+        git clone --depth 1 "${url}" "${root}"
+        validate_dependency_tree "${dir}" "$@"
+    fi
+}
+
+clone_if_missing "${OPENFOAM_REPO_URL}" OpenFOAM-14 \
+    etc/bashrc \
+    wmake/wmake \
+    src/OpenFOAM/Make/files \
+    applications/solvers/foamRun/Make/files \
+    tutorials/incompressibleFluid/venturiTube
+clone_if_missing "${THIRDPARTY_REPO_URL}" ThirdParty-14 \
+    Allwmake \
+    etc/tools
+clone_if_missing "${AIRFRANS_REPO_URL}" airfrans \
+    README.md \
+    dataset.py

 # Dummy MPI keeps this subset serial and avoids requiring mpicc/OpenMPI for p1.
 # Use SYSTEMOPENMPI instead only after installing a system MPI development package.
--- a/scripts/guard_gpu_rans_solver_when_done.sh
+++ b/scripts/guard_gpu_rans_solver_when_done.sh
@ -0,0 +1,17 @@
+#!/usr/bin/env bash
+set -euo pipefail
+
+ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd -P)"
+NOTES_FILE="${ROOT_DIR}/.loop/notes.md"
+STATUS_LINE=""
+
+if [[ -f "${NOTES_FILE}" ]]; then
+    IFS= read -r STATUS_LINE < "${NOTES_FILE}" || STATUS_LINE=""
+fi
+
+if [[ "${STATUS_LINE}" != "STATUS: DONE" ]]; then
+    printf 'Skipping full GPU RANS guard because .loop/notes.md is not STATUS: DONE (%s)\n' "${STATUS_LINE:-missing}"
+    exit 0
+fi
+
+exec "${ROOT_DIR}/scripts/verify_gpu_rans_solver.sh" "$@"
--- a/scripts/prepare_airfrans_stepper_case.py
+++ b/scripts/prepare_airfrans_stepper_case.py
@ -9,6 +9,7 @@ OpenFOAM/stepper iteration: mesh, initial fields, fvSchemes, and fvSolution.
 from __future__ import annotations

 import argparse
+import hashlib
 import json
 import math
 import re
@ -21,6 +22,51 @@ RAW_ROOT = ROOT.parent / "airfrans/data/raw/OF_dataset"
 DEFAULT_SIMULATION = "airFoil2D_SST_93.213_3.79_0.418_0.0_9.665"
 DEFAULT_SOURCE = RAW_ROOT / DEFAULT_SIMULATION
 DEFAULT_DEST = ROOT / "tmp/airfrans_stepper_case" / f"{DEFAULT_SIMULATION}_v14"
+METADATA_FILENAME = "airfrans_case_metadata.json"
+MANIFEST_FILENAME = "airfrans_case_manifest.json"
+REQUIRED_COMPARISON_FIELDS = ("U", "p", "phi", "nut", "k", "omega")
+REQUIRED_SOURCE_FILES = (
+    "system/controlDict",
+    "system/fvSchemes",
+    "system/fvSolution",
+    "system/blockMeshDict",
+    "0/U",
+    "0/p",
+    "0/nut",
+    "0/k",
+    "0/omega",
+    "constant/transportProperties",
+    "constant/turbulenceProperties",
+    "constant/polyMesh/boundary",
+    "constant/polyMesh/points.gz",
+    "constant/polyMesh/faces.gz",
+    "constant/polyMesh/owner.gz",
+    "constant/polyMesh/neighbour.gz",
+)
+SOURCE_PARAMETER_FILES = (
+    "system/controlDict",
+)
+COPIED_SOURCE_FILES = (
+    "system/fvSchemes",
+    "system/fvSolution",
+    "system/blockMeshDict",
+    "0/U",
+    "0/p",
+    "0/nut",
+    "0/k",
+    "0/omega",
+    "constant/transportProperties",
+    "constant/turbulenceProperties",
+)
+TOPOLOGY_EVIDENCE_KEYS = (
+    "n_points",
+    "n_faces",
+    "n_internal_faces",
+    "n_cells",
+    "patches",
+    "topology_sha256",
+    "geometry_sha256",
+)

 FLOAT_RE = re.compile(r"[-+]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][-+]?\d+)?")

@ -198,6 +244,109 @@ def metadata_from_source(src: Path, migrated_end_time: int) -> AirfransCaseMetad
        migrated_end_time=migrated_end_time,
    )

+def sha256_file(path: Path) -> str:
+    digest = hashlib.sha256()
+    with path.open("rb") as handle:
+        for chunk in iter(lambda: handle.read(1024 * 1024), b""):
+            digest.update(chunk)
+    return digest.hexdigest()
+
+
+def file_record(root: Path, relative: str) -> dict[str, object]:
+    path = root / relative
+    stat = path.stat()
+    return {
+        "path": relative,
+        "size": stat.st_size,
+        "sha256": sha256_file(path),
+    }
+
+
+def case_file_inventory(root: Path, *, exclude: tuple[str, ...] = ()) -> list[dict[str, object]]:
+    excluded = set(exclude)
+    records = []
+    for path in sorted(root.rglob("*")):
+        if not path.is_file():
+            continue
+        relative = path.relative_to(root).as_posix()
+        if relative in excluded:
+            continue
+        records.append(file_record(root, relative))
+    return records
+
+def source_input_inventory(src: Path) -> list[dict[str, object]]:
+    seen: set[str] = set()
+    records: list[dict[str, object]] = []
+
+    def add(relative: str) -> None:
+        if relative in seen:
+            return
+        seen.add(relative)
+        records.append(file_record(src, relative))
+
+    for relative in SOURCE_PARAMETER_FILES + COPIED_SOURCE_FILES:
+        add(relative)
+    for path in sorted((src / "constant/polyMesh").rglob("*")):
+        if path.is_file():
+            add(path.relative_to(src).as_posix())
+    return records
+
+
+
+def inventory_digest(records: list[dict[str, object]]) -> str:
+    digest = hashlib.sha256()
+    for record in records:
+        payload = json.dumps(record, sort_keys=True, separators=(",", ":")).encode("utf-8")
+        digest.update(len(payload).to_bytes(8, "little"))
+        digest.update(payload)
+    return digest.hexdigest()
+
+
+def build_case_manifest(src: Path, dst: Path, meta: AirfransCaseMetadata) -> dict[str, object]:
+    source_records = source_input_inventory(src)
+    prepared_records = case_file_inventory(dst, exclude=(MANIFEST_FILENAME,))
+    return {
+        "schema_version": 1,
+        "simulation": meta.simulation,
+        "source_case": {
+            "path": str(src),
+            "read_only": True,
+            "required_files": source_records,
+            "required_files_sha256": inventory_digest(source_records),
+        },
+        "prepared_case": {
+            "path": str(dst),
+            "openfoam_version": 14,
+            "migrated_end_time": meta.migrated_end_time,
+            "files": prepared_records,
+            "files_sha256": inventory_digest(prepared_records),
+        },
+        "comparison_contract": {
+            "fields": list(REQUIRED_COMPARISON_FIELDS),
+            "oracle_output": {
+                "producer": "foamRun -solver incompressibleFluid -noFunctionObjects",
+                "time": str(meta.migrated_end_time),
+            },
+            "repository_outputs": {
+                "run_one": "foam_stepper run_one_pimple_iteration",
+                "split": "foam_stepper split solver stages",
+                "time": str(meta.migrated_end_time),
+            },
+            "mesh_identity_evidence": list(TOPOLOGY_EVIDENCE_KEYS),
+        },
+        "metadata": asdict(meta),
+    }
+
+
+def write_case_manifest(src: Path, dst: Path, meta: AirfransCaseMetadata) -> dict[str, object]:
+    manifest = build_case_manifest(src, dst, meta)
+    (dst / MANIFEST_FILENAME).write_text(json.dumps(manifest, indent=2, sort_keys=True) + "\n")
+    return manifest
+
+
+def load_case_manifest(case: Path) -> dict[str, object]:
+    return json.loads((case / MANIFEST_FILENAME).read_text())
+

 def prepare_case(src: Path, dst: Path, *, end_time: int = 1) -> AirfransCaseMetadata:
    src = src.resolve()
@ -205,24 +354,7 @@ def prepare_case(src: Path, dst: Path, *, end_time: int = 1) -> AirfransCaseMeta
    if not src.exists():
        raise FileNotFoundError(src)

-    required = [
-        "system/fvSchemes",
-        "system/fvSolution",
-        "system/blockMeshDict",
-        "0/U",
-        "0/p",
-        "0/nut",
-        "0/k",
-        "0/omega",
-        "constant/transportProperties",
-        "constant/turbulenceProperties",
-        "constant/polyMesh/boundary",
-        "constant/polyMesh/points.gz",
-        "constant/polyMesh/faces.gz",
-        "constant/polyMesh/owner.gz",
-        "constant/polyMesh/neighbour.gz",
-    ]
-    missing = [relative for relative in required if not (src / relative).exists()]
+    missing = [relative for relative in REQUIRED_SOURCE_FILES if not (src / relative).exists()]
    if missing:
        raise FileNotFoundError(f"missing required AirfRANS files: {missing}")

@ -231,18 +363,7 @@ def prepare_case(src: Path, dst: Path, *, end_time: int = 1) -> AirfransCaseMeta
    dst.mkdir(parents=True)

    copy_required_tree(src, dst, "constant/polyMesh")
-    for relative in [
-        "system/fvSchemes",
-        "system/fvSolution",
-        "system/blockMeshDict",
-        "0/U",
-        "0/p",
-        "0/nut",
-        "0/k",
-        "0/omega",
-        "constant/transportProperties",
-        "constant/turbulenceProperties",
-    ]:
+    for relative in COPIED_SOURCE_FILES:
        copy_required_file(src, dst, relative)

    meta = metadata_from_source(src, end_time)
@ -254,7 +375,8 @@ def prepare_case(src: Path, dst: Path, *, end_time: int = 1) -> AirfransCaseMeta
    shutil.move(dst / "constant/transportProperties", dst / "constant/transportProperties.v2112")
    shutil.move(dst / "constant/turbulenceProperties", dst / "constant/turbulenceProperties.v2112")

-    (dst / "airfrans_case_metadata.json").write_text(json.dumps(asdict(meta), indent=2, sort_keys=True) + "\n")
+    (dst / METADATA_FILENAME).write_text(json.dumps(asdict(meta), indent=2, sort_keys=True) + "\n")
+    write_case_manifest(src, dst, meta)
    return meta


@ -268,6 +390,7 @@ def main() -> None:
    meta = prepare_case(args.source, args.dest, end_time=args.end_time)
    print(json.dumps(asdict(meta), indent=2, sort_keys=True))
    print(f"prepared={args.dest}")
+    print(f"manifest={args.dest / MANIFEST_FILENAME}")


 if __name__ == "__main__":
--- a/scripts/verify_airfrans_stepper.py
+++ b/scripts/verify_airfrans_stepper.py
--- a/scripts/verify_gpu_algorithm.py
+++ b/scripts/verify_gpu_algorithm.py
@ -0,0 +1,445 @@
+#!/usr/bin/env python3
+"""Run a real finite-volume GPU primitive over exported OpenFOAM arrays."""
+
+from __future__ import annotations
+
+import argparse
+import os
+import shutil
+import sys
+import time
+from dataclasses import dataclass
+from pathlib import Path
+from typing import Any, Mapping
+
+
+def _drop_ambient_pythonpath() -> None:
+    pythonpath = os.environ.pop("PYTHONPATH", "")
+    if not pythonpath:
+        return
+    for entry in pythonpath.split(os.pathsep):
+        if entry and entry in sys.path:
+            sys.path.remove(entry)
+
+
+_drop_ambient_pythonpath()
+
+import numpy as np
+import quadrants as qd
+from quadrants.profiler.kernel_profiler import get_default_kernel_profiler
+
+from verify_gpu_step_timing import json_ready, nvidia_device_identity, prepare_case, time_openfoam_step, write_report
+
+ROOT = Path(__file__).resolve().parents[1]
+DEFAULT_WORK = ROOT / "tmp/gpu_algorithm_check"
+DEFAULT_REPORT = DEFAULT_WORK / "report.json"
+ALGORITHM_NAME = "cell_flux_imbalance"
+
+
+@dataclass(frozen=True)
+class FluxInputs:
+    owner: np.ndarray
+    neighbour: np.ndarray
+    phi: np.ndarray
+    n_cells: int
+    n_internal_faces: int
+    source_summary: dict[str, Any]
+
+
+class VerificationFailure(Exception):
+    """Verifier failure with JSON-reportable details."""
+
+    def __init__(self, message: str, *, details: Mapping[str, Any] | None = None):
+        super().__init__(message)
+        self.details = dict(details or {})
+
+
+@qd.kernel
+def zero_cell_flux_imbalance(n_cells: int, out: qd.types.NDArray[qd.f64, 1]) -> None:
+    for cell in range(n_cells):
+        out[cell] = 0.0
+
+
+@qd.kernel
+def cell_flux_imbalance(
+    n_internal_faces: int,
+    owner: qd.types.NDArray[qd.i32, 1],
+    neighbour: qd.types.NDArray[qd.i32, 1],
+    phi_internal: qd.types.NDArray[qd.f64, 1],
+    out: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for face in range(n_internal_faces):
+        flux = phi_internal[face]
+        qd.atomic_add(out[owner[face]], flux)
+        qd.atomic_add(out[neighbour[face]], -flux)
+
+
+def _require_mapping(value: Any, path: str) -> Mapping[str, Any]:
+    if not isinstance(value, Mapping):
+        raise ValueError(f"{path}: expected mapping, got {type(value).__name__}")
+    return value
+
+
+def _contiguous_1d_array(value: Any, path: str) -> np.ndarray:
+    array = np.asarray(value)
+    if array.ndim != 1:
+        raise ValueError(f"{path}: expected 1-D array, got shape {array.shape}")
+    return np.ascontiguousarray(array)
+
+
+def load_openfoam_flux_inputs(case: Path) -> FluxInputs:
+    import foam_stepper as foam
+
+    state = foam.Case(case).make_stepper().export_state(required_fields=("phi",))
+    mesh = _require_mapping(state.get("mesh"), "mesh")
+    sizes = _require_mapping(mesh.get("sizes"), "mesh.sizes")
+    connectivity = _require_mapping(mesh.get("connectivity"), "mesh.connectivity")
+    fields = _require_mapping(state.get("fields"), "fields")
+    phi_field = _require_mapping(fields.get("phi"), "fields.phi")
+
+    try:
+        n_cells = int(sizes["n_cells"])
+        n_internal_faces = int(sizes["n_internal_faces"])
+    except KeyError as exc:
+        raise ValueError(f"mesh.sizes.{exc.args[0]}: missing required size") from exc
+
+    owner_raw = _contiguous_1d_array(connectivity.get("owner"), "mesh.connectivity.owner")
+    neighbour_raw = _contiguous_1d_array(connectivity.get("neighbour"), "mesh.connectivity.neighbour")
+    phi_raw = _contiguous_1d_array(phi_field.get("internal"), "fields.phi.internal")
+
+    expected_shape = (n_internal_faces,)
+    for path, array in (
+        ("mesh.connectivity.owner", owner_raw),
+        ("mesh.connectivity.neighbour", neighbour_raw),
+        ("fields.phi.internal", phi_raw),
+    ):
+        if array.shape != expected_shape:
+            raise ValueError(f"{path}: expected shape {expected_shape}, got {array.shape}")
+
+    if not np.issubdtype(owner_raw.dtype, np.integer):
+        raise ValueError(f"mesh.connectivity.owner: expected integer dtype, got {owner_raw.dtype}")
+    if not np.issubdtype(neighbour_raw.dtype, np.integer):
+        raise ValueError(f"mesh.connectivity.neighbour: expected integer dtype, got {neighbour_raw.dtype}")
+    if not np.issubdtype(phi_raw.dtype, np.floating):
+        raise ValueError(f"fields.phi.internal: expected floating dtype, got {phi_raw.dtype}")
+
+    owner_min = int(owner_raw.min(initial=0)) if owner_raw.size else 0
+    owner_max = int(owner_raw.max(initial=0)) if owner_raw.size else -1
+    neighbour_min = int(neighbour_raw.min(initial=0)) if neighbour_raw.size else 0
+    neighbour_max = int(neighbour_raw.max(initial=0)) if neighbour_raw.size else -1
+    if owner_min < 0 or neighbour_min < 0 or owner_max >= n_cells or neighbour_max >= n_cells:
+        raise ValueError(
+            "mesh.connectivity owner/neighbour indices out of cell range: "
+            f"owner=[{owner_min}, {owner_max}], neighbour=[{neighbour_min}, {neighbour_max}], n_cells={n_cells}"
+        )
+    if owner_max > np.iinfo(np.int32).max or neighbour_max > np.iinfo(np.int32).max:
+        raise ValueError("mesh.connectivity owner/neighbour exceed int32 GPU index range")
+
+    owner = np.ascontiguousarray(owner_raw.astype(np.int32, copy=False))
+    neighbour = np.ascontiguousarray(neighbour_raw.astype(np.int32, copy=False))
+    phi = np.ascontiguousarray(phi_raw.astype(np.float64, copy=False))
+
+    return FluxInputs(
+        owner=owner,
+        neighbour=neighbour,
+        phi=phi,
+        n_cells=n_cells,
+        n_internal_faces=n_internal_faces,
+        source_summary={
+            "case": str(case),
+            "mesh": {
+                "n_cells": n_cells,
+                "n_internal_faces": n_internal_faces,
+                "owner_shape": list(owner_raw.shape),
+                "owner_dtype": str(owner_raw.dtype),
+                "neighbour_shape": list(neighbour_raw.shape),
+                "neighbour_dtype": str(neighbour_raw.dtype),
+            },
+            "fields": {
+                "phi": {
+                    "entity_kind": phi_field.get("entity_kind"),
+                    "entity_count": int(phi_field.get("entity_count", -1)),
+                    "internal_shape": list(phi_raw.shape),
+                    "internal_dtype": str(phi_raw.dtype),
+                }
+            },
+            "gpu_input_dtypes": {
+                "owner": str(owner.dtype),
+                "neighbour": str(neighbour.dtype),
+                "phi_internal": str(phi.dtype),
+            },
+        },
+    )
+
+
+def run_gpu_flux_imbalance(inputs: FluxInputs, *, repeats: int) -> tuple[np.ndarray, dict[str, Any]]:
+    if repeats < 1:
+        raise ValueError("repeats must be >= 1")
+
+    qd.init(arch=qd.cuda, kernel_profiler=True)
+    owner_gpu = qd.ndarray(qd.i32, shape=inputs.owner.shape)
+    neighbour_gpu = qd.ndarray(qd.i32, shape=inputs.neighbour.shape)
+    phi_gpu = qd.ndarray(qd.f64, shape=inputs.phi.shape)
+    out_gpu = qd.ndarray(qd.f64, shape=(inputs.n_cells,))
+    owner_gpu.from_numpy(inputs.owner)
+    neighbour_gpu.from_numpy(inputs.neighbour)
+    phi_gpu.from_numpy(inputs.phi)
+
+    zero_cell_flux_imbalance(inputs.n_cells, out_gpu)
+    cell_flux_imbalance(inputs.n_internal_faces, owner_gpu, neighbour_gpu, phi_gpu, out_gpu)
+    qd.sync()
+    qd.profiler.clear_kernel_profiler_info()
+
+    start = time.perf_counter()
+    for _ in range(repeats):
+        zero_cell_flux_imbalance(inputs.n_cells, out_gpu)
+        cell_flux_imbalance(inputs.n_internal_faces, owner_gpu, neighbour_gpu, phi_gpu, out_gpu)
+    qd.sync()
+    wall_ms = (time.perf_counter() - start) * 1000.0
+
+    profiler = get_default_kernel_profiler()
+    profiler._update_records()
+    records = list(profiler._traced_records)
+    generated_kernel_names = sorted({str(record.name) for record in records})
+    if not any(ALGORITHM_NAME in name for name in generated_kernel_names):
+        raise AssertionError(
+            f"Quadrants CUDA profiler recorded no {ALGORITHM_NAME!r} kernel; recorded {generated_kernel_names}"
+        )
+    device_time_ms_total = float(sum(record.kernel_time for record in records))
+    if device_time_ms_total <= 0.0:
+        raise AssertionError(f"Quadrants CUDA profiler recorded nonpositive device time: {device_time_ms_total}")
+
+    actual = np.asarray(out_gpu.to_numpy())
+    identity = nvidia_device_identity()
+    evidence = {
+        "backend_requested": "gpu",
+        "backend_selected": "gpu",
+        "framework": "quadrants",
+        "arch_requested": "cuda",
+        "arch_selected": "cuda",
+        "device_kind": "cuda",
+        "device_name": identity["device_name"],
+        "device_uuid": identity["device_uuid"],
+        "kernel_names": [ALGORITHM_NAME],
+        "generated_kernel_names": generated_kernel_names,
+        "requested_kernel_calls": repeats,
+        "profile_record_count": len(records),
+        "used_cpu_fallback": False,
+        "synchronized_before_timing": True,
+        "synchronized_after_timing": True,
+        "wall_ms_total": wall_ms,
+        "wall_ms_per_step": wall_ms / repeats,
+        "device_time_ms_total": device_time_ms_total,
+        "device_time_ms_per_step": device_time_ms_total / repeats,
+        "device_time_ms_min_record": float(min(record.kernel_time for record in records)),
+        "device_time_ms_max_record": float(max(record.kernel_time for record in records)),
+    }
+    return actual, evidence
+
+
+def compute_cpu_flux_imbalance_reference(inputs: FluxInputs) -> np.ndarray:
+    reference = np.zeros(inputs.n_cells, dtype=np.float64)
+    np.add.at(reference, inputs.owner, inputs.phi)
+    np.add.at(reference, inputs.neighbour, -inputs.phi)
+    return reference
+
+
+def compare_flux_imbalance(actual: np.ndarray, expected: np.ndarray, *, rtol: float, atol: float) -> dict[str, Any]:
+    actual_array = np.asarray(actual)
+    expected_array = np.asarray(expected)
+    if actual_array.shape != expected_array.shape:
+        report = {
+            "allclose": False,
+            "reason": "shape_mismatch",
+            "actual_shape": list(actual_array.shape),
+            "expected_shape": list(expected_array.shape),
+            "actual_dtype": str(actual_array.dtype),
+            "expected_dtype": str(expected_array.dtype),
+            "rtol": rtol,
+            "atol": atol,
+        }
+        raise VerificationFailure("cell_flux_imbalance shape mismatch", details={"comparison": report})
+    if not np.issubdtype(actual_array.dtype, np.floating):
+        report = {
+            "allclose": False,
+            "reason": "actual_dtype_not_floating",
+            "shape": list(actual_array.shape),
+            "actual_dtype": str(actual_array.dtype),
+            "expected_dtype": str(expected_array.dtype),
+            "rtol": rtol,
+            "atol": atol,
+        }
+        raise VerificationFailure("GPU output dtype is not floating", details={"comparison": report})
+    if not np.issubdtype(expected_array.dtype, np.floating):
+        report = {
+            "allclose": False,
+            "reason": "expected_dtype_not_floating",
+            "shape": list(actual_array.shape),
+            "actual_dtype": str(actual_array.dtype),
+            "expected_dtype": str(expected_array.dtype),
+            "rtol": rtol,
+            "atol": atol,
+        }
+        raise VerificationFailure("CPU reference dtype is not floating", details={"comparison": report})
+
+    actual64 = actual_array.astype(np.float64, copy=False)
+    expected64 = expected_array.astype(np.float64, copy=False)
+    abs_diff = np.abs(actual64 - expected64)
+    max_abs = float(np.max(abs_diff)) if abs_diff.size else 0.0
+    denominator = np.maximum(np.abs(expected64), atol)
+    rel_diff = np.divide(abs_diff, denominator, out=np.zeros_like(abs_diff), where=denominator > 0.0)
+    max_rel = float(np.max(rel_diff)) if rel_diff.size else 0.0
+    largest_index = None
+    if abs_diff.size:
+        largest_index = [int(index) for index in np.unravel_index(np.argmax(abs_diff), abs_diff.shape)]
+    allclose = bool(np.allclose(actual64, expected64, rtol=rtol, atol=atol))
+    report = {
+        "allclose": allclose,
+        "shape": list(actual_array.shape),
+        "dtype": str(actual_array.dtype),
+        "actual_shape": list(actual_array.shape),
+        "expected_shape": list(expected_array.shape),
+        "actual_dtype": str(actual_array.dtype),
+        "expected_dtype": str(expected_array.dtype),
+        "rtol": rtol,
+        "atol": atol,
+        "max_abs_error": max_abs,
+        "max_rel_error": max_rel,
+        "largest_difference_index": largest_index,
+    }
+    if not allclose:
+        raise VerificationFailure("cell_flux_imbalance numerical mismatch", details={"comparison": report})
+    return report
+
+
+def run_check(args: argparse.Namespace) -> dict[str, Any]:
+    if args.work.exists():
+        shutil.rmtree(args.work)
+    args.work.mkdir(parents=True)
+
+    openfoam_case = args.work / "openfoam_step_case"
+    case = args.work / "gpu_algorithm_case"
+    prepare_case(openfoam_case)
+    prepare_case(case)
+    openfoam_step = time_openfoam_step(openfoam_case)
+    inputs = load_openfoam_flux_inputs(case)
+    cpu_reference_start = time.perf_counter()
+    expected = compute_cpu_flux_imbalance_reference(inputs)
+    cpu_reference_wall_ms = (time.perf_counter() - cpu_reference_start) * 1000.0
+    actual, gpu_evidence = run_gpu_flux_imbalance(inputs, repeats=args.repeats)
+    comparison = compare_flux_imbalance(actual, expected, rtol=args.rtol, atol=args.atol)
+    conservation_sum = float(np.sum(actual, dtype=np.float64)) if actual.size else 0.0
+    gpu_wall_per_step = gpu_evidence["wall_ms_per_step"]
+    cpu_reference_speedup = cpu_reference_wall_ms / gpu_wall_per_step if gpu_wall_per_step > 0.0 else float("inf")
+    openfoam_step_speedup = openfoam_step["wall_ms"] / gpu_wall_per_step if gpu_wall_per_step > 0.0 else float("inf")
+
+    return {
+        "status": "passed",
+        "backend_requested": "gpu",
+        "backend_selected": "gpu",
+        "device_kind": gpu_evidence["device_kind"],
+        "device_name": gpu_evidence["device_name"],
+        "device_uuid": gpu_evidence["device_uuid"],
+        "used_cpu_fallback": False,
+        "algorithm": {
+            "name": ALGORITHM_NAME,
+            "kind": "internal_face_finite_volume_primitive",
+            "formula": "cell_flux_imbalance[cell] = sum(owner phi_internal) - sum(neighbour phi_internal)",
+            "inputs": [
+                "mesh.connectivity.owner",
+                "mesh.connectivity.neighbour",
+                "mesh.sizes.n_internal_faces",
+                "fields.phi.internal",
+            ],
+            "cpu_reference_complete": True,
+        },
+        "openfoam_step": openfoam_step,
+        "cpu_reference": {
+            "operation": ALGORITHM_NAME,
+            "implementation": "numpy.add.at owner(+phi) and neighbour(-phi)",
+            "output": "cell_flux_imbalance",
+            "output_shape": list(expected.shape),
+            "output_dtype": str(expected.dtype),
+            "wall_ms": cpu_reference_wall_ms,
+        },
+        "openfoam_inputs": inputs.source_summary,
+        "gpu_algorithm": {
+            "operation": ALGORITHM_NAME,
+            "output": "cell_flux_imbalance",
+            "output_shape": list(actual.shape),
+            "output_dtype": str(actual.dtype),
+            "repeats": args.repeats,
+            "evidence": gpu_evidence,
+        },
+        "comparison": comparison,
+        "timing": {
+            "openfoam_step_wall_ms": openfoam_step["wall_ms"],
+            "openfoam_operation": openfoam_step["operation"],
+            "cpu_numpy_reference_wall_ms": cpu_reference_wall_ms,
+            "gpu_algorithm_wall_ms_total": gpu_evidence["wall_ms_total"],
+            "gpu_algorithm_wall_ms_per_step": gpu_evidence["wall_ms_per_step"],
+            "gpu_algorithm_device_ms_total": gpu_evidence["device_time_ms_total"],
+            "gpu_algorithm_device_ms_per_step": gpu_evidence["device_time_ms_per_step"],
+            "gpu_repeats": args.repeats,
+            "speedup_vs_openfoam_step_wall": openfoam_step_speedup,
+            "speedup_vs_cpu_numpy_reference_wall": cpu_reference_speedup,
+        },
+        "conservation_check": {
+            "description": "owner and neighbour scatter signs should make the global internal-face imbalance sum cancel",
+            "sum": conservation_sum,
+            "abs_sum": abs(conservation_sum),
+        },
+        "work": args.work,
+    }
+
+
+def parse_args(argv: list[str] | None = None) -> argparse.Namespace:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--work", type=Path, default=DEFAULT_WORK)
+    parser.add_argument("--report", type=Path, default=None)
+    parser.add_argument("--repeats", type=int, default=20)
+    parser.add_argument("--rtol", type=float, default=1e-10)
+    parser.add_argument("--atol", type=float, default=1e-12)
+    return parser.parse_args(argv)
+
+
+def main(argv: list[str] | None = None) -> int:
+    args = parse_args(argv)
+    report_path = args.report if args.report is not None else args.work / "report.json"
+    try:
+        report = run_check(args)
+    except Exception as exc:
+        failure = {
+            "status": "failed",
+            "backend_requested": "gpu",
+            "backend_selected": None,
+            "used_cpu_fallback": False,
+            "algorithm": {"name": ALGORITHM_NAME},
+            "work": args.work,
+            "failure": {"type": type(exc).__name__, "message": str(exc)},
+        }
+        details = getattr(exc, "details", None)
+        if details:
+            failure["failure"]["details"] = details
+        write_report(failure, report_path)
+        print(f"gpu algorithm verification failed: {exc}", file=sys.stderr)
+        print(f"report={report_path}", file=sys.stderr)
+        return 1
+
+    write_report(report, report_path)
+    evidence = report["gpu_algorithm"]["evidence"]
+    print("gpu algorithm verification passed")
+    print(f"report={report_path}")
+    print(f"algorithm={report['algorithm']['name']}")
+    print(f"backend_selected={report['backend_selected']}")
+    print(f"profile_record_count={evidence['profile_record_count']}")
+    print(f"gpu_wall_ms_per_step={evidence['wall_ms_per_step']:.6f}")
+    print(f"gpu_device_ms_per_step={evidence['device_time_ms_per_step']:.6f}")
+    print(f"openfoam_step_wall_ms={report['timing']['openfoam_step_wall_ms']:.3f}")
+    print(f"cpu_numpy_reference_wall_ms={report['timing']['cpu_numpy_reference_wall_ms']:.6f}")
+    print(f"output_shape={report['gpu_algorithm']['output_shape']}")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
--- a/scripts/verify_gpu_algorithm.sh
+++ b/scripts/verify_gpu_algorithm.sh
@ -0,0 +1,36 @@
+#!/usr/bin/env bash
+set -euo pipefail
+
+ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd -P)"
+JOBS="${JOBS:-$(nproc)}"
+PYTHON_BIN="${PYTHON_BIN:-${ROOT_DIR}/.venv/bin/python}"
+
+if [[ ! -x "${PYTHON_BIN}" ]]; then
+    uv sync --dev
+fi
+
+if ! PYTHONPATH= "${PYTHON_BIN}" -c 'import pybind11, numpy, quadrants' >/dev/null 2>&1; then
+    uv sync --dev
+fi
+
+"${ROOT_DIR}/scripts/build_openfoam_airfrans_subset.sh" >/dev/null
+
+set +u
+source "${ROOT_DIR}/OpenFOAM-14/etc/bashrc" \
+    WM_MPLIB=Dummy \
+    ParaView_TYPE=none \
+    SCOTCH_TYPE=none \
+    ZOLTAN_TYPE=none
+set -u
+unset FOAM_SIGFPE
+
+(
+    cd "${ROOT_DIR}/OpenFOAM-14"
+    wmake -j "${JOBS}" libso src/mesh/blockMesh >/dev/null
+    wmake -j "${JOBS}" applications/utilities/mesh/generation/blockMesh >/dev/null
+    wmake -j "${JOBS}" applications/utilities/mesh/manipulation/createZones >/dev/null
+)
+
+"${ROOT_DIR}/scripts/build_python_stepper.sh" >/dev/null
+
+PYTHONPATH= "${PYTHON_BIN}" "${ROOT_DIR}/scripts/verify_gpu_algorithm.py" "$@"
--- a/scripts/verify_gpu_rans_solver.sh
+++ b/scripts/verify_gpu_rans_solver.sh
@ -0,0 +1,136 @@
+#!/usr/bin/env bash
+set -euo pipefail
+
+ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd -P)"
+JOBS="${JOBS:-$(nproc)}"
+PYTHON_BIN="${PYTHON_BIN:-${ROOT_DIR}/.venv/bin/python}"
+REPORT=""
+PASSTHRU=()
+
+while [[ $# -gt 0 ]]; do
+    case "$1" in
+        --report)
+            if [[ $# -lt 2 ]]; then
+                echo "--report requires a path" >&2
+                exit 2
+            fi
+            REPORT="$2"
+            PASSTHRU+=("$1" "$2")
+            shift 2
+            ;;
+        --report=*)
+            REPORT="${1#--report=}"
+            PASSTHRU+=("$1")
+            shift
+            ;;
+        --backend|--backend=*)
+            echo "verify_gpu_rans_solver.sh always requests --backend gpu; do not pass --backend" >&2
+            exit 2
+            ;;
+        *)
+            PASSTHRU+=("$1")
+            shift
+            ;;
+    esac
+done
+
+if [[ -z "${REPORT}" ]]; then
+    echo "--report is required so the full GPU RANS guard can inspect verifier evidence" >&2
+    exit 2
+fi
+
+if [[ ! -x "${PYTHON_BIN}" ]]; then
+    uv sync --dev
+fi
+
+if ! PYTHONPATH= "${PYTHON_BIN}" -c 'import pybind11, numpy, quadrants' >/dev/null 2>&1; then
+    uv sync --dev
+fi
+
+"${ROOT_DIR}/scripts/build_openfoam_airfrans_subset.sh" >/dev/null
+
+set +u
+source "${ROOT_DIR}/OpenFOAM-14/etc/bashrc" \
+    WM_MPLIB=Dummy \
+    ParaView_TYPE=none \
+    SCOTCH_TYPE=none \
+    ZOLTAN_TYPE=none
+set -u
+unset FOAM_SIGFPE
+
+"${ROOT_DIR}/scripts/build_python_stepper.sh" >/dev/null
+
+set +e
+PYTHONPATH= "${PYTHON_BIN}" "${ROOT_DIR}/scripts/verify_airfrans_stepper.py" --backend gpu "${PASSTHRU[@]}"
+VERIFY_STATUS=$?
+set -e
+
+PYTHONPATH= "${PYTHON_BIN}" - "${REPORT}" "${VERIFY_STATUS}" <<'PY'
+from __future__ import annotations
+
+import json
+import sys
+from pathlib import Path
+
+report_path = Path(sys.argv[1])
+verify_status = int(sys.argv[2])
+required_fields = ("U", "p", "phi", "nut", "k", "omega")
+required_modes = ("run_one", "split")
+failures: list[str] = []
+
+if not report_path.exists():
+    failures.append(f"report missing: {report_path}")
+    report = {}
+else:
+    try:
+        report = json.loads(report_path.read_text())
+    except Exception as exc:  # pragma: no cover - shell guard diagnostic
+        failures.append(f"report is not valid JSON: {exc}")
+        report = {}
+
+
+def require(condition: bool, message: str) -> None:
+    if not condition:
+        failures.append(message)
+
+
+backend = report.get("backend") if isinstance(report.get("backend"), dict) else {}
+verifier_evidence = report.get("verifier_evidence") if isinstance(report.get("verifier_evidence"), dict) else {}
+modes = report.get("modes") if isinstance(report.get("modes"), dict) else {}
+provider = str(backend.get("provider") or "").lower()
+device = str(backend.get("device") or "").lower()
+
+require(verify_status == 0, f"underlying AirfRANS verifier exited {verify_status}")
+require(report.get("status") == "passed", f"report.status is {report.get('status')!r}, expected 'passed'")
+require(backend.get("requested") == "gpu", f"backend.requested is {backend.get('requested')!r}, expected 'gpu'")
+require(backend.get("selected") == "gpu", f"backend.selected is {backend.get('selected')!r}, expected 'gpu'")
+require("cpu" not in provider and provider not in {"foam_stepper_cpu", "openfoam"}, f"backend.provider looks CPU-backed: {backend.get('provider')!r}")
+require(device not in {"host", "cpu"}, f"backend.device looks CPU-backed: {backend.get('device')!r}")
+require(backend.get("counts_as_gpu_algorithm_progress") is not False, "backend explicitly says it does not count as GPU progress")
+require(verifier_evidence.get("passed") is True, "verifier_evidence.passed is not true")
+
+for mode_name in required_modes:
+    mode = modes.get(mode_name) if isinstance(modes.get(mode_name), dict) else {}
+    require(mode.get("enabled") is True, f"modes.{mode_name}.enabled is not true")
+    mode_backend = mode.get("backend") if isinstance(mode.get("backend"), dict) else backend
+    require(mode_backend.get("selected") == "gpu", f"modes.{mode_name}.backend.selected is not 'gpu'")
+    execution_path = str(mode.get("execution_path") or "").lower()
+    require("gpu" in execution_path, f"modes.{mode_name}.execution_path does not identify a GPU path: {mode.get('execution_path')!r}")
+    comparisons = mode.get("comparisons") if isinstance(mode.get("comparisons"), dict) else {}
+    for field in required_fields:
+        comparison = comparisons.get(field) if isinstance(comparisons.get(field), dict) else {}
+        require(comparison.get("allclose") is True, f"modes.{mode_name}.comparisons.{field}.allclose is not true")
+        require("actual_shape" in comparison or "shape" in comparison, f"modes.{mode_name}.comparisons.{field} has no actual shape evidence")
+        require("expected_shape" in comparison or "shape" in comparison, f"modes.{mode_name}.comparisons.{field} has no expected shape evidence")
+        require("max_abs" in comparison or "max_abs_error" in comparison, f"modes.{mode_name}.comparisons.{field} has no max error evidence")
+
+if failures:
+    print("full GPU RANS solver guard failed:", file=sys.stderr)
+    for failure in failures:
+        print(f"- {failure}", file=sys.stderr)
+    print(f"report={report_path}", file=sys.stderr)
+    sys.exit(1)
+
+print("full GPU RANS solver guard passed")
+print(f"report={report_path}")
+PY
--- a/scripts/verify_gpu_step_timing.py
+++ b/scripts/verify_gpu_step_timing.py
@ -0,0 +1,333 @@
+#!/usr/bin/env python3
+"""Compare OpenFOAM step timing with a proven GPU kernel over OpenFOAM-derived data.
+
+This is deliberately not a RANS-solver speedup claim. The GPU step here is the
+smallest honest executable GPU operation we can verify today: a Quadrants CUDA
+kernel consuming the OpenFOAM velocity field and producing per-cell squared
+speed. The report records OpenFOAM one-iteration timing next to the GPU kernel
+timing and fails unless the GPU kernel actually ran on the CUDA backend.
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+import os
+import shutil
+import subprocess
+import sys
+import time
+from pathlib import Path
+from typing import Any, Mapping
+
+
+def _drop_ambient_pythonpath() -> None:
+    pythonpath = os.environ.pop("PYTHONPATH", "")
+    if not pythonpath:
+        return
+    for entry in pythonpath.split(os.pathsep):
+        if entry and entry in sys.path:
+            sys.path.remove(entry)
+
+
+_drop_ambient_pythonpath()
+
+import numpy as np
+import quadrants as qd
+from quadrants.profiler.kernel_profiler import get_default_kernel_profiler
+
+from openfoam_env import apply_openfoam_env, openfoam_env
+
+ROOT = Path(__file__).resolve().parents[1]
+TUTORIAL = ROOT / "OpenFOAM-14/tutorials/incompressibleFluid/venturiTube"
+DEFAULT_WORK = ROOT / "tmp/gpu_step_timing_check"
+DEFAULT_REPORT = DEFAULT_WORK / "report.json"
+KERNEL_NAME = "squared_speed"
+HARNESS_SCOPE = (
+    "squared_speed is only a CUDA proof/timing harness over OpenFOAM-derived "
+    "velocity data; it is not the target finite-volume or RANS GPU algorithm"
+)
+
+
+@qd.kernel
+def squared_speed(n: int, u: qd.types.NDArray[qd.f32, 2], out: qd.types.NDArray[qd.f32, 1]) -> None:
+    for i in range(n):
+        ux = u[i, 0]
+        uy = u[i, 1]
+        uz = u[i, 2]
+        out[i] = ux * ux + uy * uy + uz * uz
+
+
+def json_ready(value: Any) -> Any:
+    if isinstance(value, Mapping):
+        return {str(key): json_ready(item) for key, item in value.items()}
+    if isinstance(value, (list, tuple)):
+        return [json_ready(item) for item in value]
+    if isinstance(value, np.generic):
+        return value.item()
+    if isinstance(value, np.ndarray):
+        return value.tolist()
+    if isinstance(value, Path):
+        return str(value)
+    return value
+
+
+def patch_control_dict(case: Path) -> None:
+    path = case / "system/controlDict"
+    text = path.read_text()
+    replacements = {
+        "startFrom       startTime;": "startFrom       startTime;",
+        "endTime         1000;": "endTime         1;",
+        "writeInterval   50;": "writeInterval   1;",
+    }
+    for old, new in replacements.items():
+        text = text.replace(old, new)
+    path.write_text(text)
+
+
+def run_openfoam_command(cmd: list[str]) -> None:
+    subprocess.run(cmd, cwd=ROOT, check=True, stdout=subprocess.DEVNULL, stderr=subprocess.STDOUT, env=openfoam_env())
+
+
+def prepare_case(dst: Path) -> None:
+    if dst.exists():
+        shutil.rmtree(dst)
+    shutil.copytree(TUTORIAL, dst, ignore=shutil.ignore_patterns("processor*", "postProcessing", "*.log"))
+    for orig in (dst / "0").glob("*.orig"):
+        shutil.copyfile(orig, orig.with_suffix(""))
+    patch_control_dict(dst)
+    run_openfoam_command(["blockMesh", "-case", str(dst)])
+    run_openfoam_command(["createZones", "-case", str(dst)])
+
+
+def import_foam() -> Any:
+    apply_openfoam_env()
+    import foam_stepper as foam
+
+    return foam
+
+
+def time_openfoam_step(case: Path) -> dict[str, Any]:
+    foam = import_foam()
+    stepper = foam.Case(case).make_stepper()
+    start = time.perf_counter()
+    result = stepper.run_one_pimple_iteration()
+    wall_ms = (time.perf_counter() - start) * 1000.0
+    fields = result.outputs["fields"]
+    return {
+        "operation": "foam_stepper.run_one_pimple_iteration",
+        "backend": "OpenFOAM C++ CPU stepper",
+        "wall_ms": wall_ms,
+        "output_shapes": {
+            "U": list(np.asarray(fields["U"].internal).shape),
+            "p": list(np.asarray(fields["p"].internal).shape),
+            "phi": list(np.asarray(fields["phi"].internal).shape),
+        },
+    }
+
+
+def load_openfoam_velocity(case: Path) -> np.ndarray:
+    foam = import_foam()
+    fields = foam.Case(case).make_stepper().fields()
+    velocity = np.asarray(fields.U.internal, dtype=np.float32)
+    if velocity.ndim != 2 or velocity.shape[1] != 3:
+        raise AssertionError(f"expected vector U field with shape (cells, 3), got {velocity.shape}")
+    return np.ascontiguousarray(velocity)
+
+
+def nvidia_device_identity() -> dict[str, str | None]:
+    try:
+        completed = subprocess.run(
+            ["nvidia-smi", "--query-gpu=name,uuid", "--format=csv,noheader,nounits"],
+            check=True,
+            stdout=subprocess.PIPE,
+            stderr=subprocess.DEVNULL,
+            text=True,
+        )
+    except (FileNotFoundError, subprocess.CalledProcessError):
+        return {"device_name": None, "device_uuid": None}
+
+    first = completed.stdout.strip().splitlines()[0] if completed.stdout.strip() else ""
+    if not first:
+        return {"device_name": None, "device_uuid": None}
+    parts = [part.strip() for part in first.split(",", 1)]
+    return {"device_name": parts[0], "device_uuid": parts[1] if len(parts) > 1 else None}
+
+
+def compare_arrays(actual: np.ndarray, expected: np.ndarray, *, rtol: float, atol: float) -> dict[str, Any]:
+    if actual.shape != expected.shape:
+        raise AssertionError(f"GPU output shape mismatch: {actual.shape} != {expected.shape}")
+    abs_diff = np.abs(actual - expected)
+    max_abs = float(np.max(abs_diff)) if abs_diff.size else 0.0
+    max_rel = float(np.max(abs_diff / np.maximum(np.abs(expected), atol))) if abs_diff.size else 0.0
+    largest_index = None
+    if abs_diff.size:
+        largest_index = [int(index) for index in np.unravel_index(np.argmax(abs_diff), abs_diff.shape)]
+    allclose = bool(np.allclose(actual, expected, rtol=rtol, atol=atol))
+    report = {
+        "allclose": allclose,
+        "shape": list(actual.shape),
+        "dtype": str(actual.dtype),
+        "rtol": rtol,
+        "atol": atol,
+        "max_abs_error": max_abs,
+        "max_rel_error": max_rel,
+        "largest_difference_index": largest_index,
+    }
+    if not allclose:
+        raise AssertionError(f"GPU output mismatch: {json.dumps(report, sort_keys=True)}")
+    return report
+
+
+def run_gpu_velocity_step(velocity: np.ndarray, *, repeats: int) -> tuple[np.ndarray, dict[str, Any]]:
+    if repeats < 1:
+        raise ValueError("repeats must be >= 1")
+
+    qd.init(arch=qd.cuda, kernel_profiler=True)
+    cells = int(velocity.shape[0])
+    u_gpu = qd.ndarray(qd.f32, shape=velocity.shape)
+    out_gpu = qd.ndarray(qd.f32, shape=(cells,))
+    u_gpu.from_numpy(velocity)
+    squared_speed(cells, u_gpu, out_gpu)
+    qd.sync()
+    qd.profiler.clear_kernel_profiler_info()
+
+    start = time.perf_counter()
+    for _ in range(repeats):
+        squared_speed(cells, u_gpu, out_gpu)
+    qd.sync()
+    wall_ms = (time.perf_counter() - start) * 1000.0
+
+    profiler = get_default_kernel_profiler()
+    profiler._update_records()
+    records = list(profiler._traced_records)
+    generated_kernel_names = sorted({str(record.name) for record in records})
+    if not any(KERNEL_NAME in name for name in generated_kernel_names):
+        raise AssertionError(f"Quadrants CUDA profiler recorded no {KERNEL_NAME!r} kernel; recorded {generated_kernel_names}")
+    device_time_ms_total = float(sum(record.kernel_time for record in records))
+    if device_time_ms_total <= 0.0:
+        raise AssertionError(f"Quadrants CUDA profiler recorded nonpositive device time: {device_time_ms_total}")
+    actual = out_gpu.to_numpy()
+
+    identity = nvidia_device_identity()
+    evidence = {
+        "backend_requested": "gpu",
+        "backend_selected": "gpu",
+        "framework": "quadrants",
+        "arch_requested": "cuda",
+        "arch_selected": "cuda",
+        "device_kind": "cuda",
+        "device_name": identity["device_name"],
+        "device_uuid": identity["device_uuid"],
+        "kernel_names": [KERNEL_NAME],
+        "generated_kernel_names": generated_kernel_names,
+        "requested_kernel_calls": repeats,
+        "profile_record_count": len(records),
+        "used_cpu_fallback": False,
+        "synchronized_before_timing": True,
+        "synchronized_after_timing": True,
+        "wall_ms_total": wall_ms,
+        "wall_ms_per_step": wall_ms / repeats,
+        "device_time_ms_total": device_time_ms_total,
+        "device_time_ms_per_step": device_time_ms_total / repeats,
+        "device_time_ms_min_record": float(min(record.kernel_time for record in records)),
+        "device_time_ms_max_record": float(max(record.kernel_time for record in records)),
+    }
+    return np.asarray(actual), evidence
+
+
+def write_report(report: Mapping[str, Any], path: Path) -> None:
+    path.parent.mkdir(parents=True, exist_ok=True)
+    path.write_text(json.dumps(json_ready(report), indent=2, sort_keys=True, allow_nan=False) + "\n")
+
+
+def run_check(args: argparse.Namespace) -> dict[str, Any]:
+    if args.work.exists():
+        shutil.rmtree(args.work)
+    args.work.mkdir(parents=True)
+
+    openfoam_case = args.work / "openfoam_step_case"
+    gpu_case = args.work / "gpu_input_case"
+    prepare_case(openfoam_case)
+    prepare_case(gpu_case)
+
+    openfoam_step = time_openfoam_step(openfoam_case)
+    velocity = load_openfoam_velocity(gpu_case)
+    expected = np.einsum("ij,ij->i", velocity, velocity).astype(np.float32, copy=False)
+    actual, gpu_evidence = run_gpu_velocity_step(velocity, repeats=args.repeats)
+    comparison = compare_arrays(actual, expected, rtol=args.rtol, atol=args.atol)
+
+    speedup_vs_openfoam_wall = openfoam_step["wall_ms"] / gpu_evidence["wall_ms_per_step"] if gpu_evidence["wall_ms_per_step"] > 0 else float("inf")
+    report = {
+        "status": "passed",
+        "scope": HARNESS_SCOPE,
+        "algorithm": {
+            "name": KERNEL_NAME,
+            "role": "proof_timing_harness",
+            "target_algorithm_complete": False,
+        },
+        "work": args.work,
+        "openfoam_step": openfoam_step,
+        "gpu_step": {
+            "operation": KERNEL_NAME,
+            "input_field": "U",
+            "input_shape": list(velocity.shape),
+            "output": "squared_speed_per_cell",
+            "output_shape": list(actual.shape),
+            "output_dtype": str(actual.dtype),
+            "repeats": args.repeats,
+            "evidence": gpu_evidence,
+        },
+        "comparison": comparison,
+        "timing": {
+            "openfoam_wall_ms": openfoam_step["wall_ms"],
+            "gpu_wall_ms_total": gpu_evidence["wall_ms_total"],
+            "gpu_wall_ms_per_step": gpu_evidence["wall_ms_per_step"],
+            "gpu_device_ms_per_step": gpu_evidence["device_time_ms_per_step"],
+            "speedup_vs_openfoam_wall": speedup_vs_openfoam_wall,
+        },
+    }
+    return report
+
+
+def parse_args(argv: list[str] | None = None) -> argparse.Namespace:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--work", type=Path, default=DEFAULT_WORK)
+    parser.add_argument("--report", type=Path, default=None)
+    parser.add_argument("--repeats", type=int, default=20)
+    parser.add_argument("--rtol", type=float, default=5e-6)
+    parser.add_argument("--atol", type=float, default=1e-6)
+    return parser.parse_args(argv)
+
+
+def main(argv: list[str] | None = None) -> int:
+    args = parse_args(argv)
+    report_path = args.report if args.report is not None else args.work / "report.json"
+    try:
+        report = run_check(args)
+    except Exception as exc:
+        failure = {
+            "status": "failed",
+            "scope": HARNESS_SCOPE,
+            "work": args.work,
+            "failure": {"type": type(exc).__name__, "message": str(exc)},
+        }
+        write_report(failure, report_path)
+        print(f"gpu step timing verification failed: {exc}", file=sys.stderr)
+        print(f"report={report_path}", file=sys.stderr)
+        return 1
+
+    write_report(report, report_path)
+    print("gpu step timing verification passed")
+    print(f"report={report_path}")
+    print(f"openfoam_wall_ms={report['timing']['openfoam_wall_ms']:.3f}")
+    print(f"gpu_wall_ms_per_step={report['timing']['gpu_wall_ms_per_step']:.6f}")
+    print(f"gpu_device_ms_per_step={report['timing']['gpu_device_ms_per_step']:.6f}")
+    print(f"speedup_vs_openfoam_wall={report['timing']['speedup_vs_openfoam_wall']:.3f}")
+    print(f"backend_selected={report['gpu_step']['evidence']['backend_selected']}")
+    print(f"profile_record_count={report['gpu_step']['evidence']['profile_record_count']}")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())
--- a/scripts/verify_gpu_step_timing.sh
+++ b/scripts/verify_gpu_step_timing.sh
@ -0,0 +1,36 @@
+#!/usr/bin/env bash
+set -euo pipefail
+
+ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd -P)"
+JOBS="${JOBS:-$(nproc)}"
+PYTHON_BIN="${PYTHON_BIN:-${ROOT_DIR}/.venv/bin/python}"
+
+if [[ ! -x "${PYTHON_BIN}" ]]; then
+    uv sync --dev
+fi
+
+if ! PYTHONPATH= "${PYTHON_BIN}" -c 'import pybind11, numpy, quadrants' >/dev/null 2>&1; then
+    uv sync --dev
+fi
+
+"${ROOT_DIR}/scripts/build_openfoam_airfrans_subset.sh" >/dev/null
+
+set +u
+source "${ROOT_DIR}/OpenFOAM-14/etc/bashrc" \
+    WM_MPLIB=Dummy \
+    ParaView_TYPE=none \
+    SCOTCH_TYPE=none \
+    ZOLTAN_TYPE=none
+set -u
+unset FOAM_SIGFPE
+
+(
+    cd "${ROOT_DIR}/OpenFOAM-14"
+    wmake -j "${JOBS}" libso src/mesh/blockMesh >/dev/null
+    wmake -j "${JOBS}" applications/utilities/mesh/generation/blockMesh >/dev/null
+    wmake -j "${JOBS}" applications/utilities/mesh/manipulation/createZones >/dev/null
+)
+
+"${ROOT_DIR}/scripts/build_python_stepper.sh" >/dev/null
+
+PYTHONPATH= "${PYTHON_BIN}" "${ROOT_DIR}/scripts/verify_gpu_step_timing.py" "$@"