feat: some progress on gpu implementation, more granularity on solver numerical comparison

.gitignore

@@ -4,6 +4,14 @@
 python/src/foam_stepper.egg-info/
 tmp/
 
+# local planning/spec context injected by the loop harness
+CURRENT_STATUS.md
+MATCHING_AIRFRANS_SIMULATION.md
+PROJECT_GROUNDING.md
+PYTHON_STEPPER_SPEC.md
+STAGE_ORACLE_SPEC.md
+VERIFIER_HARNESS_SPEC.md
+
 # foreign
 OpenFOAM-14/
 ThirdParty-14/

python/src/foam_stepper/gpu/constants.py

@@ -19,9 +19,10 @@ GPU_STAGE_CAPABILITY_PREFIX = "solver_stage_contract:"
 GPU_INPUT_SCHEMA_VERSION = 1
 DEFAULT_LAMINAR_NU = 1.5e-5
 DEFAULT_MOMENTUM_RELAXATION_ALPHA = 0.9
-DEFAULT_MOMENTUM_PBICGSTAB_ITERATIONS = 50
+DEFAULT_MOMENTUM_PBICGSTAB_ITERATIONS = 100
 DEFAULT_MOMENTUM_PBICGSTAB_RESIDUAL_TOLERANCE_SQUARED = 1.0e-16
-DEFAULT_PRESSURE_CG_ITERATIONS = 300
+DEFAULT_PRESSURE_CG_ITERATIONS = 2000
+DEFAULT_PRESSURE_NON_ORTHOGONAL_CORRECTORS = 3
 DEFAULT_SIMPLE_CONSISTENT_RATU_FACTOR = 10.0
 DEFAULT_OMEGA_WALL_BETA1 = 0.075
 
@@ -34,13 +35,20 @@ STAGE_OBSERVABILITY_GROUPS = (
             "assemble_momentum_terms": ("terms",),
             "assemble_UEqn": ("UEqn",),
         },
+        "run_one_outputs": {
+            "assemble_UEqn": ("UEqn",),
+        },
     },
     {
         "name": "pressure_assembly",
         "split_stages": ("compute_pressure_inputs", "assemble_pEqn"),
         "run_one_stages": ("compute_pressure_inputs", "assemble_pEqn"),
         "split_outputs": {
-            "compute_pressure_inputs": ("HbyA", "phiHbyA", "rAU"),
+            "compute_pressure_inputs": ("HbyA", "phiHbyA", "rAU", "rAtU"),
+            "assemble_pEqn": ("pEqn",),
+        },
+        "run_one_outputs": {
+            "compute_pressure_inputs": ("HbyA", "phiHbyA", "rAU", "rAtU"),
             "assemble_pEqn": ("pEqn",),
         },
     },
@@ -52,6 +60,10 @@ STAGE_OBSERVABILITY_GROUPS = (
             "solve_UEqn": ("performance", "field_after"),
             "solve_pEqn": ("performance", "p", "phi"),
         },
+        "run_one_outputs": {
+            "solve_UEqn": ("performance", "field_after"),
+            "solve_pEqn": ("performance", "p", "phi"),
+        },
     },
     {
         "name": "final_correction",
@@ -60,6 +72,10 @@ STAGE_OBSERVABILITY_GROUPS = (
         "split_outputs": {
             "correct_velocity_pressure_flux": ("U", "p", "phi"),
         },
+        "run_one_outputs": {
+            "update_phi_from_pEqn_flux": ("phi",),
+            "correct_velocity_pressure_flux": ("U", "p", "phi"),
+        },
     },
     {
         "name": "turbulence_updates",
@@ -69,13 +85,17 @@ STAGE_OBSERVABILITY_GROUPS = (
             "momentum_transport_predict": ("case_path", "solver_name"),
             "momentum_transport_correct": ("U", "p", "phi", "nut", "k", "omega"),
         },
+        "run_one_outputs": {
+            "momentum_transport_predict": ("case_path", "solver_name"),
+            "momentum_transport_correct": ("U", "p", "phi", "nut", "k", "omega"),
+        },
     },
 )
 
 GPU_STAGE_KERNEL_ENTRYPOINTS = {
-    "momentum_assembly": ["gpu_rans_momentum_assembly", "gpu_rans_momentum_diffusion_coefficients", "gpu_rans_momentum_wall_diffusion_coefficients", "gpu_rans_momentum_convection_coefficients", "gpu_rans_momentum_equation_relaxation"],
+    "momentum_assembly": ["gpu_rans_momentum_assembly", "gpu_rans_momentum_diffusion_coefficients", "gpu_rans_momentum_wall_diffusion_coefficients", "gpu_rans_momentum_convection_coefficients", "gpu_rans_momentum_bounded_convection_sp_internal", "gpu_rans_momentum_bounded_convection_sp_boundary", "gpu_rans_momentum_convection_boundary_coefficients", "gpu_rans_momentum_convection_boundary_source", "gpu_rans_momentum_boundary_internal_diag", "gpu_rans_momentum_boundary_relaxation_coefficients", "gpu_rans_add_scalar_field", "gpu_rans_zero_tensor_field", "gpu_rans_momentum_gauss_grad_u_internal", "gpu_rans_momentum_gauss_grad_u_boundary", "gpu_rans_momentum_gauss_grad_u_finish", "gpu_rans_momentum_linear_upwind_source", "gpu_rans_zero_scalar_field", "gpu_rans_momentum_offdiag_abs_accumulate", "gpu_rans_momentum_equation_relaxation", "gpu_rans_momentum_pressure_gradient_source", "gpu_rans_momentum_pressure_boundary_source"],
-    "pressure_assembly": ["gpu_rans_pressure_inputs", "gpu_rans_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"],
+    "pressure_assembly": ["gpu_rans_pressure_inputs", "gpu_rans_momentum_h1_face_accumulate", "gpu_rans_momentum_h1_finish", "gpu_rans_consistent_rAtU", "gpu_rans_momentum_hbyA_source", "gpu_rans_momentum_hbyA_face_accumulate", "gpu_rans_momentum_hbyA_finish", "gpu_rans_surface_flux_from_cells", "gpu_rans_consistent_phiHbyA_correction", "gpu_rans_pressure_assembly", "gpu_rans_pressure_laplacian_coefficients", "gpu_rans_pressure_mixed_boundary_laplacian", "gpu_rans_pressure_source_from_flux", "gpu_rans_pressure_source_from_boundary_flux"],
-    "linear_solve_results": ["gpu_ldu_matvec_vector_asymmetric_diag", "gpu_ldu_matvec_vector_asymmetric_face_accumulate", "gpu_bicgstab_initialize_vector", "gpu_bicgstab_dot_vector", "gpu_bicgstab_update_direction_vector", "gpu_bicgstab_precondition_vector", "gpu_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"],
+    "linear_solve_results": ["gpu_ldu_matvec_vector_asymmetric_diag", "gpu_ldu_matvec_vector_asymmetric_face_accumulate", "gpu_bicgstab_initialize_vector", "gpu_bicgstab_dot_vector", "gpu_bicgstab_update_direction_vector", "gpu_bicgstab_precondition_vector", "gpu_dilu_apply_vector_asymmetric_faces", "gpu_bicgstab_update_intermediate_vector_preconditioned", "gpu_bicgstab_update_solution_residual_vector_preconditioned", "gpu_vector_residual_squared", "gpu_rans_negate_scalar_field", "gpu_ldu_matvec_scalar_symmetric_diag", "gpu_ldu_matvec_scalar_symmetric_face_accumulate", "gpu_pcg_initialize_scalar", "gpu_cg_dot_scalar", "gpu_pcg_update_solution_residual_scalar", "gpu_pcg_update_direction_scalar"],
     "final_correction": ["gpu_rans_pressure_flux_correction", "gpu_rans_final_correction", "gpu_rans_pressure_velocity_correction"],
     "turbulence_updates": ["gpu_rans_turbulence_update", "gpu_rans_omega_wall_update"],
 }

python/src/foam_stepper/gpu/kernels.py

@@ -10,14 +10,16 @@ from .constants import GPU_STAGE_KERNEL_ENTRYPOINTS
 def gpu_rans_momentum_assembly(
     n_cells: int,
     u_internal: qd.types.NDArray[qd.f64, 2],
+    cell_volumes: qd.types.NDArray[qd.f64, 1],
     diag: qd.types.NDArray[qd.f64, 1],
     source: qd.types.NDArray[qd.f64, 2],
 ) -> None:
     for cell in range(n_cells):
-        diag[cell] = 1.0
+        # fvSchemes uses steadyState ddt; fvm::ddt(U) contributes no diagonal.
-        source[cell, 0] = u_internal[cell, 0]
+        diag[cell] = 0.0
-        source[cell, 1] = u_internal[cell, 1]
+        source[cell, 0] = 0.0
-        source[cell, 2] = u_internal[cell, 2]
+        source[cell, 1] = 0.0
+        source[cell, 2] = 0.0
 
 @qd.kernel
 def gpu_rans_momentum_diffusion_coefficients(
@@ -64,7 +66,9 @@ def gpu_rans_momentum_wall_diffusion_coefficients(
     face_centres: qd.types.NDArray[qd.f64, 2],
     face_area_vectors: qd.types.NDArray[qd.f64, 2],
     face_area_magnitudes: qd.types.NDArray[qd.f64, 1],
-    diag: qd.types.NDArray[qd.f64, 1],
+    boundary_relax_add: qd.types.NDArray[qd.f64, 1],
+    boundary_relax_subtract: qd.types.NDArray[qd.f64, 1],
+    boundary_diag: qd.types.NDArray[qd.f64, 1],
     source: qd.types.NDArray[qd.f64, 2],
 ) -> None:
     for boundary_face in range(n_boundary_faces):
@@ -79,11 +83,15 @@ def gpu_rans_momentum_wall_diffusion_coefficients(
             projected_delta = 1.0e-300
         mag_sf = face_area_magnitudes[boundary_face]
         coeff = (laminar_nu + nut_boundary[boundary_face]) * mag_sf * mag_sf / projected_delta
-        qd.atomic_add(diag[cell], coeff)
+        qd.atomic_add(boundary_relax_add[cell], coeff)
+        qd.atomic_add(boundary_relax_subtract[cell], coeff)
+        qd.atomic_add(boundary_diag[cell], coeff)
         qd.atomic_add(source[cell, 0], coeff * boundary_values[boundary_face, 0])
         qd.atomic_add(source[cell, 1], coeff * boundary_values[boundary_face, 1])
         qd.atomic_add(source[cell, 2], coeff * boundary_values[boundary_face, 2])
 
+
+
 @qd.kernel
 def gpu_rans_momentum_convection_coefficients(
     n_internal_faces: int,
@@ -106,25 +114,321 @@ def gpu_rans_momentum_convection_coefficients(
             qd.atomic_add(upper[face], flux)
 
 
+@qd.kernel
+def gpu_rans_momentum_bounded_convection_sp_internal(
+    n_internal_faces: int,
+    owner: qd.types.NDArray[qd.i32, 1],
+    neighbour: qd.types.NDArray[qd.i32, 1],
+    phi: qd.types.NDArray[qd.f64, 1],
+    diag: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for face in range(n_internal_faces):
+        owner_cell = owner[face]
+        neighbour_cell = neighbour[face]
+        flux = phi[face]
+        qd.atomic_add(diag[owner_cell], -flux)
+        qd.atomic_add(diag[neighbour_cell], flux)
+
+
+@qd.kernel
+def gpu_rans_momentum_bounded_convection_sp_boundary(
+    n_boundary_faces: int,
+    face_cells: qd.types.NDArray[qd.i32, 1],
+    phi_boundary: qd.types.NDArray[qd.f64, 1],
+    diag: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for face in range(n_boundary_faces):
+        cell = face_cells[face]
+        qd.atomic_add(diag[cell], -phi_boundary[face])
+
+
+@qd.kernel
+def gpu_rans_momentum_convection_boundary_coefficients(
+    n_boundary_faces: int,
+    face_cells: qd.types.NDArray[qd.i32, 1],
+    phi_boundary: qd.types.NDArray[qd.f64, 1],
+    value_internal_coeffs: qd.types.NDArray[qd.f64, 2],
+    value_boundary_coeffs: qd.types.NDArray[qd.f64, 2],
+    diag: qd.types.NDArray[qd.f64, 1],
+    source: qd.types.NDArray[qd.f64, 2],
+) -> None:
+    for face in range(n_boundary_faces):
+        cell = face_cells[face]
+        flux = phi_boundary[face]
+        internal_average = (
+            value_internal_coeffs[face, 0]
+            + value_internal_coeffs[face, 1]
+            + value_internal_coeffs[face, 2]
+        ) / 3.0
+        qd.atomic_add(diag[cell], flux * internal_average)
+        qd.atomic_add(source[cell, 0], -flux * value_boundary_coeffs[face, 0])
+        qd.atomic_add(source[cell, 1], -flux * value_boundary_coeffs[face, 1])
+        qd.atomic_add(source[cell, 2], -flux * value_boundary_coeffs[face, 2])
+
+
+@qd.kernel
+def gpu_rans_momentum_convection_boundary_source(
+    n_boundary_faces: int,
+    face_cells: qd.types.NDArray[qd.i32, 1],
+    boundary_source: qd.types.NDArray[qd.f64, 2],
+    source: qd.types.NDArray[qd.f64, 2],
+) -> None:
+    for face in range(n_boundary_faces):
+        cell = face_cells[face]
+        qd.atomic_add(source[cell, 0], boundary_source[face, 0])
+        qd.atomic_add(source[cell, 1], boundary_source[face, 1])
+        qd.atomic_add(source[cell, 2], boundary_source[face, 2])
+
+
+@qd.kernel
+def gpu_rans_momentum_boundary_internal_diag(
+    n_boundary_faces: int,
+    face_cells: qd.types.NDArray[qd.i32, 1],
+    internal_coeffs: qd.types.NDArray[qd.f64, 2],
+    diag: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for face in range(n_boundary_faces):
+        cell = face_cells[face]
+        coeff = (
+            internal_coeffs[face, 0]
+            + internal_coeffs[face, 1]
+            + internal_coeffs[face, 2]
+        ) / 3.0
+        qd.atomic_add(diag[cell], coeff)
+
+
+@qd.kernel
+def gpu_rans_momentum_boundary_relaxation_coefficients(
+    n_boundary_faces: int,
+    face_cells: qd.types.NDArray[qd.i32, 1],
+    internal_coeffs: qd.types.NDArray[qd.f64, 2],
+    boundary_relax_add: qd.types.NDArray[qd.f64, 1],
+    boundary_relax_subtract: qd.types.NDArray[qd.f64, 1],
+    boundary_diag: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for face in range(n_boundary_faces):
+        cell = face_cells[face]
+        c0 = internal_coeffs[face, 0]
+        c1 = internal_coeffs[face, 1]
+        c2 = internal_coeffs[face, 2]
+        a0 = c0
+        if a0 < 0.0:
+            a0 = -a0
+        a1 = c1
+        if a1 < 0.0:
+            a1 = -a1
+        a2 = c2
+        if a2 < 0.0:
+            a2 = -a2
+        max_abs = a0
+        if a1 > max_abs:
+            max_abs = a1
+        if a2 > max_abs:
+            max_abs = a2
+        min_coeff = c0
+        if c1 < min_coeff:
+            min_coeff = c1
+        if c2 < min_coeff:
+            min_coeff = c2
+        qd.atomic_add(boundary_relax_add[cell], max_abs)
+        qd.atomic_add(boundary_relax_subtract[cell], min_coeff)
+        qd.atomic_add(boundary_diag[cell], (c0 + c1 + c2) / 3.0)
+
+
+@qd.kernel
+def gpu_rans_add_scalar_field(
+    n_cells: int,
+    addend: qd.types.NDArray[qd.f64, 1],
+    field: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for cell in range(n_cells):
+        field[cell] += addend[cell]
+
+
+@qd.kernel
+def gpu_rans_zero_tensor_field(
+    n_cells: int,
+    tensor: qd.types.NDArray[qd.f64, 3],
+) -> None:
+    for cell in range(n_cells):
+        for component in range(3):
+            for direction in range(3):
+                tensor[cell, component, direction] = 0.0
+
+
+@qd.kernel
+def gpu_rans_momentum_gauss_grad_u_internal(
+    n_internal_faces: int,
+    owner: qd.types.NDArray[qd.i32, 1],
+    neighbour: qd.types.NDArray[qd.i32, 1],
+    u_internal: qd.types.NDArray[qd.f64, 2],
+    sf: qd.types.NDArray[qd.f64, 2],
+    grad_u: qd.types.NDArray[qd.f64, 3],
+) -> None:
+    for face in range(n_internal_faces):
+        owner_cell = owner[face]
+        neighbour_cell = neighbour[face]
+        for component in range(3):
+            face_value = 0.5 * (u_internal[owner_cell, component] + u_internal[neighbour_cell, component])
+            for direction in range(3):
+                flux_value = face_value * sf[face, direction]
+                qd.atomic_add(grad_u[owner_cell, component, direction], flux_value)
+                qd.atomic_add(grad_u[neighbour_cell, component, direction], -flux_value)
+
+
+@qd.kernel
+def gpu_rans_momentum_gauss_grad_u_boundary(
+    n_boundary_faces: int,
+    face_cells: qd.types.NDArray[qd.i32, 1],
+    u_boundary: qd.types.NDArray[qd.f64, 2],
+    sf_boundary: qd.types.NDArray[qd.f64, 2],
+    grad_u: qd.types.NDArray[qd.f64, 3],
+) -> None:
+    for face in range(n_boundary_faces):
+        cell = face_cells[face]
+        for component in range(3):
+            face_value = u_boundary[face, component]
+            for direction in range(3):
+                qd.atomic_add(grad_u[cell, component, direction], face_value * sf_boundary[face, direction])
+
+
+@qd.kernel
+def gpu_rans_momentum_gauss_grad_u_finish(
+    n_cells: int,
+    cell_volumes: qd.types.NDArray[qd.f64, 1],
+    grad_u: qd.types.NDArray[qd.f64, 3],
+) -> None:
+    for cell in range(n_cells):
+        inv_volume = 1.0 / cell_volumes[cell]
+        for component in range(3):
+            for direction in range(3):
+                grad_u[cell, component, direction] *= inv_volume
+
+
+@qd.kernel
+def gpu_rans_momentum_linear_upwind_source(
+    n_internal_faces: int,
+    owner: qd.types.NDArray[qd.i32, 1],
+    neighbour: qd.types.NDArray[qd.i32, 1],
+    phi: qd.types.NDArray[qd.f64, 1],
+    cell_centres: qd.types.NDArray[qd.f64, 2],
+    face_centres: qd.types.NDArray[qd.f64, 2],
+    grad_u: qd.types.NDArray[qd.f64, 3],
+    source: qd.types.NDArray[qd.f64, 2],
+) -> None:
+    for face in range(n_internal_faces):
+        owner_cell = owner[face]
+        neighbour_cell = neighbour[face]
+        flux = phi[face]
+        upwind_cell = owner_cell
+        if flux < 0.0:
+            upwind_cell = neighbour_cell
+        dx0 = face_centres[face, 0] - cell_centres[upwind_cell, 0]
+        dx1 = face_centres[face, 1] - cell_centres[upwind_cell, 1]
+        dx2 = face_centres[face, 2] - cell_centres[upwind_cell, 2]
+        for component in range(3):
+            correction = (
+                dx0 * grad_u[upwind_cell, component, 0]
+                + dx1 * grad_u[upwind_cell, component, 1]
+                + dx2 * grad_u[upwind_cell, component, 2]
+            )
+            flux_correction = flux * correction
+            qd.atomic_add(source[owner_cell, component], -flux_correction)
+            qd.atomic_add(source[neighbour_cell, component], flux_correction)
+
+
+@qd.kernel
+def gpu_rans_zero_scalar_field(
+    n_cells: int,
+    field: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for cell in range(n_cells):
+        field[cell] = 0.0
+
+
+@qd.kernel
+def gpu_rans_momentum_offdiag_abs_accumulate(
+    n_internal_faces: int,
+    owner: qd.types.NDArray[qd.i32, 1],
+    neighbour: qd.types.NDArray[qd.i32, 1],
+    upper: qd.types.NDArray[qd.f64, 1],
+    lower: qd.types.NDArray[qd.f64, 1],
+    offdiag_sum: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for face in range(n_internal_faces):
+        upper_abs = upper[face]
+        if upper_abs < 0.0:
+            upper_abs = -upper_abs
+        lower_abs = lower[face]
+        if lower_abs < 0.0:
+            lower_abs = -lower_abs
+        qd.atomic_add(offdiag_sum[owner[face]], upper_abs)
+        qd.atomic_add(offdiag_sum[neighbour[face]], lower_abs)
+
 
 @qd.kernel
 def gpu_rans_momentum_equation_relaxation(
     n_cells: int,
     u_internal: qd.types.NDArray[qd.f64, 2],
     alpha: float,
+    offdiag_sum: qd.types.NDArray[qd.f64, 1],
+    boundary_relax_add: qd.types.NDArray[qd.f64, 1],
+    boundary_relax_subtract: qd.types.NDArray[qd.f64, 1],
     diag: qd.types.NDArray[qd.f64, 1],
     source: qd.types.NDArray[qd.f64, 2],
 ) -> None:
     for cell in range(n_cells):
-        old_diag = diag[cell]
+        raw_diag = diag[cell]
-        relaxed_diag = old_diag / alpha
+        dominant_diag = raw_diag + boundary_relax_add[cell]
-        source_scale = relaxed_diag - old_diag
+        if dominant_diag < 0.0:
+            dominant_diag = -dominant_diag
+        if offdiag_sum[cell] > dominant_diag:
+            dominant_diag = offdiag_sum[cell]
+        relaxed_diag = dominant_diag / alpha - boundary_relax_subtract[cell]
+        source_scale = relaxed_diag - raw_diag
         diag[cell] = relaxed_diag
         source[cell, 0] += source_scale * u_internal[cell, 0]
         source[cell, 1] += source_scale * u_internal[cell, 1]
         source[cell, 2] += source_scale * u_internal[cell, 2]
 
 
+@qd.kernel
+def gpu_rans_momentum_pressure_gradient_source(
+    n_internal_faces: int,
+    owner: qd.types.NDArray[qd.i32, 1],
+    neighbour: qd.types.NDArray[qd.i32, 1],
+    p_internal: qd.types.NDArray[qd.f64, 1],
+    sf: qd.types.NDArray[qd.f64, 2],
+    source: qd.types.NDArray[qd.f64, 2],
+) -> None:
+    for face in range(n_internal_faces):
+        owner_cell = owner[face]
+        neighbour_cell = neighbour[face]
+        p_face = 0.5 * (p_internal[owner_cell] + p_internal[neighbour_cell])
+        qd.atomic_add(source[owner_cell, 0], -p_face * sf[face, 0])
+        qd.atomic_add(source[owner_cell, 1], -p_face * sf[face, 1])
+        qd.atomic_add(source[owner_cell, 2], -p_face * sf[face, 2])
+        qd.atomic_add(source[neighbour_cell, 0], p_face * sf[face, 0])
+        qd.atomic_add(source[neighbour_cell, 1], p_face * sf[face, 1])
+        qd.atomic_add(source[neighbour_cell, 2], p_face * sf[face, 2])
+
+
+@qd.kernel
+def gpu_rans_momentum_pressure_boundary_source(
+    n_boundary_faces: int,
+    face_cells: qd.types.NDArray[qd.i32, 1],
+    p_boundary: qd.types.NDArray[qd.f64, 1],
+    sf_boundary: qd.types.NDArray[qd.f64, 2],
+    source: qd.types.NDArray[qd.f64, 2],
+) -> None:
+    for face in range(n_boundary_faces):
+        cell = face_cells[face]
+        p_face = p_boundary[face]
+        qd.atomic_add(source[cell, 0], -p_face * sf_boundary[face, 0])
+        qd.atomic_add(source[cell, 1], -p_face * sf_boundary[face, 1])
+        qd.atomic_add(source[cell, 2], -p_face * sf_boundary[face, 2])
+
+
 @qd.kernel
 def gpu_rans_momentum_hbyA_source(
     n_cells: int,
@@ -173,18 +477,7 @@ def gpu_rans_momentum_hbyA_finish(
         HbyA[cell, 2] *= inv_diag
 
 
-@qd.kernel
+
-def gpu_rans_momentum_hbyA_fixed_value_boundary(
-    n_boundary_faces: int,
-    face_cells: qd.types.NDArray[qd.i32, 1],
-    boundary_values: qd.types.NDArray[qd.f64, 2],
-    HbyA: qd.types.NDArray[qd.f64, 2],
-) -> None:
-    for boundary_face in range(n_boundary_faces):
-        cell = face_cells[boundary_face]
-        HbyA[cell, 0] = boundary_values[boundary_face, 0]
-        HbyA[cell, 1] = boundary_values[boundary_face, 1]
-        HbyA[cell, 2] = boundary_values[boundary_face, 2]
 
 
 @qd.kernel
@@ -193,20 +486,53 @@ def gpu_rans_pressure_inputs(
     u_diag: qd.types.NDArray[qd.f64, 1],
     cell_volumes: qd.types.NDArray[qd.f64, 1],
     rAU: qd.types.NDArray[qd.f64, 1],
+    H1: qd.types.NDArray[qd.f64, 1],
 ) -> None:
     for cell in range(n_cells):
         rAU[cell] = cell_volumes[cell] / u_diag[cell]
+        H1[cell] = 0.0
+
+
+@qd.kernel
+def gpu_rans_momentum_h1_face_accumulate(
+    n_internal_faces: int,
+    owner: qd.types.NDArray[qd.i32, 1],
+    neighbour: qd.types.NDArray[qd.i32, 1],
+    upper: qd.types.NDArray[qd.f64, 1],
+    lower: qd.types.NDArray[qd.f64, 1],
+    H1: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for face in range(n_internal_faces):
+        owner_cell = owner[face]
+        neighbour_cell = neighbour[face]
+        qd.atomic_add(H1[owner_cell], -upper[face])
+        qd.atomic_add(H1[neighbour_cell], -lower[face])
+
+
+@qd.kernel
+def gpu_rans_momentum_h1_finish(
+    n_cells: int,
+    cell_volumes: qd.types.NDArray[qd.f64, 1],
+    H1: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for cell in range(n_cells):
+        H1[cell] /= cell_volumes[cell]
 
 
 @qd.kernel
 def gpu_rans_consistent_rAtU(
     n_cells: int,
     rAU: qd.types.NDArray[qd.f64, 1],
-    ratu_factor: float,
+    H1: qd.types.NDArray[qd.f64, 1],
     rAtU: qd.types.NDArray[qd.f64, 1],
 ) -> None:
     for cell in range(n_cells):
-        rAtU[cell] = ratu_factor * rAU[cell]
+        inv_rAU = 1.0 / rAU[cell]
+        floor = 0.1 * inv_rAU
+        denominator = inv_rAU - H1[cell]
+        if denominator < floor:
+            denominator = floor
+        rAtU[cell] = 1.0 / denominator
 
 
 @qd.kernel
@@ -228,6 +554,39 @@ def gpu_rans_surface_flux_from_cells(
         )
 
 
+
+@qd.kernel
+def gpu_rans_consistent_phiHbyA_correction(
+    n_internal_faces: int,
+    owner: qd.types.NDArray[qd.i32, 1],
+    neighbour: qd.types.NDArray[qd.i32, 1],
+    rAU: qd.types.NDArray[qd.f64, 1],
+    rAtU: qd.types.NDArray[qd.f64, 1],
+    p_internal: qd.types.NDArray[qd.f64, 1],
+    cell_centres: qd.types.NDArray[qd.f64, 2],
+    sf: qd.types.NDArray[qd.f64, 2],
+    mag_sf: qd.types.NDArray[qd.f64, 1],
+    phiHbyA: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for face in range(n_internal_faces):
+        owner_cell = owner[face]
+        neighbour_cell = neighbour[face]
+        dx0 = cell_centres[neighbour_cell, 0] - cell_centres[owner_cell, 0]
+        dx1 = cell_centres[neighbour_cell, 1] - cell_centres[owner_cell, 1]
+        dx2 = cell_centres[neighbour_cell, 2] - cell_centres[owner_cell, 2]
+        projected_delta = dx0 * sf[face, 0] + dx1 * sf[face, 1] + dx2 * sf[face, 2]
+        if projected_delta < 0.0:
+            projected_delta = -projected_delta
+        if projected_delta < 1.0e-300:
+            projected_delta = 1.0e-300
+        interpolated_delta = 0.5 * (
+            (rAtU[owner_cell] - rAU[owner_cell])
+            + (rAtU[neighbour_cell] - rAU[neighbour_cell])
+        )
+        pressure_jump = p_internal[neighbour_cell] - p_internal[owner_cell]
+        phiHbyA[face] += interpolated_delta * pressure_jump * mag_sf[face] * mag_sf[face] / projected_delta
+
+
 @qd.kernel
 def gpu_rans_pressure_assembly(
     n_cells: int,
@@ -311,6 +670,16 @@ def gpu_rans_pressure_mixed_boundary_laplacian(
         qd.atomic_add(source[cell], -coeff * boundary_values[face])
 
 
+@qd.kernel
+def gpu_rans_negate_scalar_field(
+    n_values: int,
+    source: qd.types.NDArray[qd.f64, 1],
+    out: qd.types.NDArray[qd.f64, 1],
+) -> None:
+    for index in range(n_values):
+        out[index] = -source[index]
+
+
 @qd.kernel
 def gpu_rans_face_flux_copy(
     n_internal_faces: int,
@@ -424,17 +793,34 @@ __all__ = [
     "gpu_rans_momentum_assembly",
     "gpu_rans_momentum_diffusion_coefficients",
     "gpu_rans_momentum_convection_coefficients",
+    "gpu_rans_momentum_bounded_convection_sp_internal",
+    "gpu_rans_momentum_bounded_convection_sp_boundary",
+    "gpu_rans_momentum_convection_boundary_coefficients",
+    "gpu_rans_momentum_convection_boundary_source",
+    "gpu_rans_momentum_boundary_internal_diag",
+    "gpu_rans_momentum_boundary_relaxation_coefficients",
+    "gpu_rans_add_scalar_field",
+    "gpu_rans_zero_tensor_field",
+    "gpu_rans_momentum_gauss_grad_u_internal",
+    "gpu_rans_momentum_gauss_grad_u_boundary",
+    "gpu_rans_momentum_gauss_grad_u_finish",
+    "gpu_rans_momentum_linear_upwind_source",
     "gpu_rans_momentum_equation_relaxation",
+    "gpu_rans_momentum_pressure_gradient_source",
     "gpu_rans_momentum_hbyA_source",
     "gpu_rans_momentum_hbyA_face_accumulate",
     "gpu_rans_momentum_hbyA_finish",
     "gpu_rans_pressure_inputs",
+    "gpu_rans_momentum_h1_face_accumulate",
+    "gpu_rans_momentum_h1_finish",
     "gpu_rans_consistent_rAtU",
+    "gpu_rans_consistent_phiHbyA_correction",
     "gpu_rans_pressure_assembly",
     "gpu_rans_pressure_laplacian_coefficients",
     "gpu_rans_pressure_source_from_flux",
     "gpu_rans_pressure_source_from_boundary_flux",
     "gpu_rans_pressure_mixed_boundary_laplacian",
+    "gpu_rans_negate_scalar_field",
     "gpu_rans_face_flux_copy",
     "gpu_rans_surface_flux_from_cells",
     "gpu_rans_pressure_flux_correction",

python/src/foam_stepper/gpu/linear_solve.py

@@ -21,6 +21,20 @@ class GpuLduCsr:
     metadata: Mapping[str, Any]
 
 
+@dataclasses.dataclass(frozen=True)
+class GpuLduLevelSchedule:
+    """Device level schedule for triangular sweeps over OpenFOAM LDU faces."""
+
+    level_offsets: Any
+    level_cells: Any
+    incoming_offsets: Any
+    incoming_faces: Any
+    outgoing_offsets: Any
+    outgoing_faces: Any
+    n_levels: int
+    metadata: Mapping[str, Any]
+
+
 @qd.kernel
 def gpu_ldu_jacobi_scalar(
     n_cells: int,
@@ -332,16 +346,18 @@ def gpu_ldu_pcg_scalar_symmetric_faces(
     operator_work: Any,
     residual_squared: Any,
     denominator: Any,
+    preconditioner_diag: Any | None = None,
     *,
     iterations: int,
     residual_tolerance_squared: float = 0.0,
 ) -> dict[str, Any]:
-    """Run diagonal-preconditioned GPU CG for a symmetric per-face LDU matrix."""
+    """Run GPU CG for a symmetric per-face LDU matrix with a diagonal preconditioner."""
 
+    precond_diag = diag if preconditioner_diag is None else preconditioner_diag
     gpu_ldu_matvec_scalar_symmetric_faces(n_cells, n_internal_faces, owner, neighbour, upper, diag, initial, operator_work)
     gpu_zero_scalar_accumulator(residual_squared)
     gpu_zero_scalar_accumulator(denominator)
-    gpu_pcg_initialize_scalar(n_cells, source, operator_work, diag, initial, out, residual, preconditioned_residual, direction, residual_squared, denominator)
+    gpu_pcg_initialize_scalar(n_cells, source, operator_work, precond_diag, initial, out, residual, preconditioned_residual, direction, residual_squared, denominator)
     qd.sync()
     rr_value = float(np.asarray(residual_squared.to_numpy())[0])
     rho_value = float(np.asarray(denominator.to_numpy())[0])
@@ -362,7 +378,7 @@ def gpu_ldu_pcg_scalar_symmetric_faces(
         alpha = rho_value / denominator_value
         gpu_zero_scalar_accumulator(residual_squared)
         gpu_zero_scalar_accumulator(denominator)
-        gpu_pcg_update_solution_residual_scalar(n_cells, alpha, diag, out, direction, residual, operator_work, preconditioned_residual, residual_squared, denominator)
+        gpu_pcg_update_solution_residual_scalar(n_cells, alpha, precond_diag, out, direction, residual, operator_work, preconditioned_residual, residual_squared, denominator)
         qd.sync()
         next_rr_value = float(np.asarray(residual_squared.to_numpy())[0])
         next_rho_value = float(np.asarray(denominator.to_numpy())[0])
@@ -382,7 +398,7 @@ def gpu_ldu_pcg_scalar_symmetric_faces(
         "final_preconditioned_dot": rho_value,
         "converged": rr_value <= residual_tolerance_squared,
         "residual_tolerance_squared": residual_tolerance_squared,
-        "preconditioner": "diagonal_jacobi",
+        "preconditioner": "diagonal_jacobi" if preconditioner_diag is None else "dic_reciprocal_diagonal",
     }
 def gpu_ldu_jacobi_scalar_symmetric_faces(
     n_cells: int,
@@ -552,6 +568,8 @@ def gpu_bicgstab_dot_vector(
         )
 
 
+
+
 @qd.kernel
 def gpu_bicgstab_update_direction_vector(
     n_cells: int,
@@ -642,6 +660,148 @@ def gpu_bicgstab_precondition_vector(
             out[cell, 2] = source[cell, 2]
 
 
+@qd.kernel
+def gpu_dilu_apply_vector_asymmetric_faces(
+    n_cells: int,
+    n_internal_faces: int,
+    owner: qd.types.NDArray[qd.i32, 1],
+    neighbour: qd.types.NDArray[qd.i32, 1],
+    losort: qd.types.NDArray[qd.i32, 1],
+    upper: qd.types.NDArray[qd.f64, 1],
+    lower: qd.types.NDArray[qd.f64, 1],
+    reciprocal_diag: qd.types.NDArray[qd.f64, 1],
+    source: qd.types.NDArray[qd.f64, 2],
+    out: qd.types.NDArray[qd.f64, 2],
+) -> None:
+    """Apply OpenFOAM DILU forward/back substitution to one vector residual."""
+
+    for worker in range(1):
+        for cell in range(n_cells):
+            scale = reciprocal_diag[cell]
+            out[cell, 0] = scale * source[cell, 0]
+            out[cell, 1] = scale * source[cell, 1]
+            out[cell, 2] = scale * source[cell, 2]
+
+        for sorted_index in range(n_internal_faces):
+            face = losort[sorted_index]
+            owner_cell = owner[face]
+            neighbour_cell = neighbour[face]
+            scale = reciprocal_diag[neighbour_cell] * lower[face]
+            out[neighbour_cell, 0] -= scale * out[owner_cell, 0]
+            out[neighbour_cell, 1] -= scale * out[owner_cell, 1]
+            out[neighbour_cell, 2] -= scale * out[owner_cell, 2]
+
+        for reverse_index in range(n_internal_faces):
+            face = n_internal_faces - 1 - reverse_index
+            owner_cell = owner[face]
+            neighbour_cell = neighbour[face]
+            scale = reciprocal_diag[owner_cell] * upper[face]
+            out[owner_cell, 0] -= scale * out[neighbour_cell, 0]
+            out[owner_cell, 1] -= scale * out[neighbour_cell, 1]
+            out[owner_cell, 2] -= scale * out[neighbour_cell, 2]
+
+@qd.kernel
+def gpu_dilu_forward_level_vector(
+    level: int,
+    level_offsets: qd.types.NDArray[qd.i32, 1],
+    level_cells: qd.types.NDArray[qd.i32, 1],
+    owner: qd.types.NDArray[qd.i32, 1],
+    incoming_offsets: qd.types.NDArray[qd.i32, 1],
+    incoming_faces: qd.types.NDArray[qd.i32, 1],
+    lower: qd.types.NDArray[qd.f64, 1],
+    reciprocal_diag: qd.types.NDArray[qd.f64, 1],
+    source: qd.types.NDArray[qd.f64, 2],
+    out: qd.types.NDArray[qd.f64, 2],
+) -> None:
+    for index in range(level_offsets[level], level_offsets[level + 1]):
+        cell = level_cells[index]
+        scale = reciprocal_diag[cell]
+        value_x = scale * source[cell, 0]
+        value_y = scale * source[cell, 1]
+        value_z = scale * source[cell, 2]
+        for face_slot in range(incoming_offsets[cell], incoming_offsets[cell + 1]):
+            face = incoming_faces[face_slot]
+            owner_cell = owner[face]
+            coeff = scale * lower[face]
+            value_x -= coeff * out[owner_cell, 0]
+            value_y -= coeff * out[owner_cell, 1]
+            value_z -= coeff * out[owner_cell, 2]
+        out[cell, 0] = value_x
+        out[cell, 1] = value_y
+        out[cell, 2] = value_z
+
+
+@qd.kernel
+def gpu_dilu_backward_level_vector(
+    level: int,
+    level_offsets: qd.types.NDArray[qd.i32, 1],
+    level_cells: qd.types.NDArray[qd.i32, 1],
+    neighbour: qd.types.NDArray[qd.i32, 1],
+    outgoing_offsets: qd.types.NDArray[qd.i32, 1],
+    outgoing_faces: qd.types.NDArray[qd.i32, 1],
+    upper: qd.types.NDArray[qd.f64, 1],
+    reciprocal_diag: qd.types.NDArray[qd.f64, 1],
+    out: qd.types.NDArray[qd.f64, 2],
+) -> None:
+    for index in range(level_offsets[level], level_offsets[level + 1]):
+        cell = level_cells[index]
+        scale = reciprocal_diag[cell]
+        value_x = out[cell, 0]
+        value_y = out[cell, 1]
+        value_z = out[cell, 2]
+        begin = outgoing_offsets[cell]
+        end = outgoing_offsets[cell + 1]
+        for reverse_slot in range(end - begin):
+            face = outgoing_faces[end - 1 - reverse_slot]
+            neighbour_cell = neighbour[face]
+            coeff = scale * upper[face]
+            value_x -= coeff * out[neighbour_cell, 0]
+            value_y -= coeff * out[neighbour_cell, 1]
+            value_z -= coeff * out[neighbour_cell, 2]
+        out[cell, 0] = value_x
+        out[cell, 1] = value_y
+        out[cell, 2] = value_z
+
+
+
+def gpu_dilu_apply_vector_asymmetric_levels(
+    schedule: GpuLduLevelSchedule,
+    owner: Any,
+    neighbour: Any,
+    upper: Any,
+    lower: Any,
+    reciprocal_diag: Any,
+    source: Any,
+    out: Any,
+) -> None:
+    """Apply OpenFOAM DILU using parallel cell work within each dependency level."""
+
+    for level in range(schedule.n_levels):
+        gpu_dilu_forward_level_vector(
+            level,
+            schedule.level_offsets,
+            schedule.level_cells,
+            owner,
+            schedule.incoming_offsets,
+            schedule.incoming_faces,
+            lower,
+            reciprocal_diag,
+            source,
+            out,
+        )
+    for reverse_level in range(schedule.n_levels):
+        level = schedule.n_levels - 1 - reverse_level
+        gpu_dilu_backward_level_vector(
+            level,
+            schedule.level_offsets,
+            schedule.level_cells,
+            neighbour,
+            schedule.outgoing_offsets,
+            schedule.outgoing_faces,
+            upper,
+            reciprocal_diag,
+            out,
+        )
 @qd.kernel
 def gpu_bicgstab_update_intermediate_vector_preconditioned(
     n_cells: int,
@@ -1032,11 +1192,12 @@ def gpu_ldu_pbicgstab_vector_asymmetric_faces(
     denominator: Any,
     omega_numerator: Any,
     omega_denominator: Any,
+    preconditioner_diag: Any | None = None,
     *,
     iterations: int,
     residual_tolerance_squared: float = 0.0,
 ) -> dict[str, Any]:
-    """Run diagonal-preconditioned GPU BiCGStab for an asymmetric vector LDU matrix."""
+    """Run GPU BiCGStab for an asymmetric vector LDU matrix with a diagonal preconditioner."""
 
     gpu_ldu_matvec_vector_asymmetric_faces(
         n_cells,
@@ -1068,6 +1229,7 @@ def gpu_ldu_pbicgstab_vector_asymmetric_faces(
     alpha = 1.0
     omega = 1.0
     performed_iterations = 0
+    precond_diag = diag if preconditioner_diag is None else preconditioner_diag
 
     for _ in range(iterations):
         if residual_value <= residual_tolerance_squared:
@@ -1088,8 +1250,19 @@ def gpu_ldu_pbicgstab_vector_asymmetric_faces(
             beta = (rho_new / rho_old) * (alpha / omega)
 
         gpu_bicgstab_update_direction_vector(n_cells, beta, omega, residual, direction, operator_direction)
-        gpu_bicgstab_precondition_vector(n_cells, diag, direction, intermediate)
+        gpu_bicgstab_precondition_vector(n_cells, precond_diag, direction, intermediate)
-        gpu_ldu_matvec_vector_asymmetric_faces(n_cells, n_internal_faces, owner, neighbour, upper, lower, diag, intermediate, operator_direction)
+        gpu_ldu_matvec_vector_asymmetric_faces(
+            n_cells,
+            n_internal_faces,
+            owner,
+            neighbour,
+            upper,
+            lower,
+            diag,
+            intermediate,
+            operator_direction,
+        )
+
         gpu_zero_scalar_accumulator(denominator)
         gpu_bicgstab_dot_vector(n_cells, shadow, operator_direction, denominator)
         qd.sync()
@@ -1116,7 +1289,7 @@ def gpu_ldu_pbicgstab_vector_asymmetric_faces(
             residual_value = intermediate_residual
             break
 
-        gpu_bicgstab_precondition_vector(n_cells, diag, intermediate, residual)
+        gpu_bicgstab_precondition_vector(n_cells, precond_diag, intermediate, residual)
         gpu_ldu_matvec_vector_asymmetric_faces(n_cells, n_internal_faces, owner, neighbour, upper, lower, diag, residual, operator_intermediate)
         gpu_zero_scalar_accumulator(omega_numerator)
         gpu_zero_scalar_accumulator(omega_denominator)
@@ -1149,7 +1322,7 @@ def gpu_ldu_pbicgstab_vector_asymmetric_faces(
         "final_residual_squared": residual_value,
         "converged": residual_value <= residual_tolerance_squared,
         "residual_tolerance_squared": residual_tolerance_squared,
-        "preconditioner": "diagonal_jacobi",
+        "preconditioner": "diagonal_jacobi" if preconditioner_diag is None else "dilu_reciprocal_diagonal",
     }
 
 def gpu_ldu_jacobi_vector_symmetric_faces(
@@ -1228,6 +1401,28 @@ def build_ldu_csr(
     return offsets, columns, coefficients
 
 
+def build_losort_addr(n_cells: int, neighbour: np.ndarray) -> np.ndarray:
+    """Build OpenFOAM lduAddressing::losortAddr from upper/neighbour cells."""
+
+    neighbour_i32 = np.asarray(neighbour, dtype=np.int32).reshape(-1)
+    if neighbour_i32.size == 0:
+        return np.zeros(0, dtype=np.int32)
+    if int(neighbour_i32.min()) < 0 or int(neighbour_i32.max()) >= n_cells:
+        raise ValueError("neighbour addresses exceed cell range")
+
+    counts = np.bincount(neighbour_i32, minlength=n_cells).astype(np.int32, copy=False)
+    offsets = np.empty(n_cells + 1, dtype=np.int32)
+    offsets[0] = 0
+    np.cumsum(counts, out=offsets[1:])
+    losort = np.empty(int(offsets[-1]), dtype=np.int32)
+    cursor = offsets[:-1].copy()
+    for face, neighbour_cell in enumerate(neighbour_i32):
+        slot = int(cursor[int(neighbour_cell)])
+        losort[slot] = int(face)
+        cursor[int(neighbour_cell)] += 1
+    return np.ascontiguousarray(losort, dtype=np.int32)
+
+
 def _gpu_i32(values: np.ndarray) -> Any:
     host = np.ascontiguousarray(values.astype(np.int32, copy=False))
     alloc_shape = host.shape if host.size else (1,)
@@ -1274,6 +1469,7 @@ def empty_ldu_csr_gpu(n_cells: int, name: str) -> GpuLduCsr:
 __all__ = [
     "GpuLduCsr",
     "build_ldu_csr",
+    "build_losort_addr",
     "empty_ldu_csr_gpu",
     "gpu_bicgstab_dot_vector",
     "gpu_bicgstab_initialize_vector",
@@ -1285,6 +1481,7 @@ __all__ = [
     "gpu_bicgstab_update_solution_residual_vector_preconditioned",
     "gpu_copy_scalar",
     "gpu_copy_vector",
+    "gpu_dilu_apply_vector_asymmetric_faces",
     "gpu_ldu_bicgstab_vector_asymmetric_faces",
     "gpu_ldu_pbicgstab_vector_asymmetric_faces",
     "gpu_ldu_cg_scalar_symmetric_faces",
@@ -1306,5 +1503,6 @@ __all__ = [
     "gpu_scalar_jacobi_finish",
     "gpu_vector_jacobi_finish",
     "gpu_vector_residual_squared",
+    "gpu_zero_scalar_accumulator",
     "ldu_csr_to_gpu",
 ]

scripts/update_loop_diagnostic_context.py

@@ -0,0 +1,430 @@
+#!/usr/bin/env python3
+"""Generate concise loop context from GPU RANS verifier reports."""
+
+from __future__ import annotations
+
+import argparse
+import json
+import math
+from datetime import datetime, timezone
+from pathlib import Path
+from typing import Any, Mapping
+
+ROOT = Path(__file__).resolve().parents[1]
+DEFAULT_CONTEXT = ROOT / ".loop/diagnostic-context.md"
+DEFAULT_BASELINE = ROOT / ".loop/diagnostic-baseline.json"
+REPORT_GLOBS = (
+    "tmp/**/verifier_report.json",
+    "tmp/**/report.json",
+)
+TMP_REPORT_GLOBS = (
+    "*gpu*/report.json",
+    "*gpu*/verifier_report.json",
+    "worker_gpu_rans_solver*/report.json",
+    "judge_gpu_rans_solver*/report.json",
+    "*gpu*rans*/report.json",
+    "*gpu*rans*/verifier_report.json",
+)
+
+
+def load_json(path: Path) -> Any | None:
+    try:
+        return json.loads(path.read_text())
+    except Exception:
+        return None
+
+
+def is_verifier_report(value: Any) -> bool:
+    if not isinstance(value, Mapping):
+        return False
+    return any(key in value for key in ("verifier_evidence", "first_divergence_summary", "artifact_comparisons"))
+
+
+def is_gpu_report(value: Mapping[str, Any], path: Path) -> bool:
+    backend = value.get("backend") if isinstance(value.get("backend"), Mapping) else {}
+    requested = str(backend.get("requested") or "").lower()
+    selected = str(backend.get("selected") or "").lower()
+    return requested == "gpu" or selected == "gpu" or "gpu" in path.parent.name.lower()
+
+
+def discover_latest_report(root: Path) -> tuple[Path | None, Mapping[str, Any] | None]:
+    candidates: list[Path] = []
+    for pattern in REPORT_GLOBS:
+        candidates.extend(root.glob(pattern))
+    tmp_root = Path("/tmp")
+    if tmp_root.exists():
+        for pattern in TMP_REPORT_GLOBS:
+            candidates.extend(tmp_root.glob(pattern))
+    unique = sorted({path.resolve() for path in candidates if path.is_file()}, key=lambda path: path.stat().st_mtime, reverse=True)
+    reports: list[tuple[Path, Mapping[str, Any]]] = []
+    for path in unique:
+        data = load_json(path)
+        if is_verifier_report(data):
+            reports.append((path, data))  # type: ignore[arg-type]
+    if not reports:
+        return None, None
+    gpu_reports = [(path, data) for path, data in reports if is_gpu_report(data, path)]
+    return (gpu_reports or reports)[0]
+
+
+def get_path(value: Mapping[str, Any], path: str) -> Any:
+    cursor: Any = value
+    for part in path.split("."):
+        if not isinstance(cursor, Mapping):
+            return None
+        cursor = cursor.get(part)
+    return cursor
+
+
+def finite_number(value: Any) -> float | None:
+    try:
+        number = float(value)
+    except (TypeError, ValueError):
+        return None
+    return number if math.isfinite(number) else None
+
+
+def fmt(value: Any) -> str:
+    number = finite_number(value)
+    if number is None:
+        if value is True:
+            return "true"
+        if value is False:
+            return "false"
+        if value is None:
+            return "-"
+        return str(value)
+    if number == 0.0:
+        return "0"
+    if abs(number) >= 1e4 or abs(number) < 1e-3:
+        return f"{number:.3e}"
+    return f"{number:.6g}"
+
+
+def status_word(value: Any) -> str:
+    if value is True:
+        return "passed"
+    if value is False:
+        return "failed"
+    return str(value or "unknown")
+
+
+def first_nested_key(value: Any, target: str, path: str = "") -> tuple[str, Mapping[str, Any]] | None:
+    if isinstance(value, Mapping):
+        for key, item in value.items():
+            next_path = f"{path}.{key}" if path else str(key)
+            if key == target and isinstance(item, Mapping):
+                return next_path, item
+            found = first_nested_key(item, target, next_path)
+            if found is not None:
+                return found
+    elif isinstance(value, list):
+        for index, item in enumerate(value):
+            found = first_nested_key(item, target, f"{path}[{index}]")
+            if found is not None:
+                return found
+    return None
+
+
+def artifact_checks(report: Mapping[str, Any]) -> list[dict[str, Any]]:
+    out: list[dict[str, Any]] = []
+    comparisons = report.get("artifact_comparisons", {}) if isinstance(report.get("artifact_comparisons"), Mapping) else {}
+    for family in ("pressure_inputs", "matrix_operator", "solver"):
+        family_report = comparisons.get(family) if isinstance(comparisons.get(family), Mapping) else {}
+        for check in family_report.get("checks", []) if isinstance(family_report.get("checks"), list) else []:
+            if not isinstance(check, Mapping):
+                continue
+            largest = check.get("largest_difference") if isinstance(check.get("largest_difference"), Mapping) else {}
+            location = largest.get("location") if isinstance(largest.get("location"), Mapping) else {}
+            out.append(
+                {
+                    "family": family,
+                    "name": check.get("name"),
+                    "allclose": check.get("allclose"),
+                    "reason": check.get("reason"),
+                    "max_abs": check.get("max_abs"),
+                    "mean_abs": check.get("mean_abs"),
+                    "rms_abs": check.get("rms_abs"),
+                    "location": location,
+                }
+            )
+    out.sort(key=lambda item: (item.get("allclose") is True, item["family"], str(item.get("name"))))
+    return out
+
+
+def field_checks(report: Mapping[str, Any]) -> list[dict[str, Any]]:
+    out: list[dict[str, Any]] = []
+    modes = report.get("modes", {}) if isinstance(report.get("modes"), Mapping) else {}
+    for mode_name in ("run_one", "split"):
+        mode = modes.get(mode_name) if isinstance(modes.get(mode_name), Mapping) else {}
+        comparisons = mode.get("comparisons", {}) if isinstance(mode.get("comparisons"), Mapping) else {}
+        for field, data in comparisons.items():
+            if not isinstance(data, Mapping):
+                continue
+            out.append(
+                {
+                    "mode": mode_name,
+                    "field": field,
+                    "allclose": data.get("allclose"),
+                    "max_abs": data.get("max_abs"),
+                    "mean_abs": data.get("mean_abs"),
+                    "rms_abs": data.get("rms_abs"),
+                    "location": data.get("location"),
+                }
+            )
+    out.sort(key=lambda item: (item.get("allclose") is True, item["mode"], str(item.get("field"))))
+    return out
+
+
+def extract_metrics(report: Mapping[str, Any]) -> dict[str, Any]:
+    metrics: dict[str, Any] = {
+        "status.passed": report.get("status") == "passed",
+        "verifier.passed": get_path(report, "verifier_evidence.passed") is True,
+    }
+    first = report.get("first_divergence_summary") if isinstance(report.get("first_divergence_summary"), Mapping) else {}
+    if first:
+        metrics["first.target"] = first.get("first_target")
+        metrics["first.family"] = first.get("artifact_family")
+        metrics["first.stage_group"] = first.get("stage_group")
+    for check in artifact_checks(report):
+        base = f"artifact.{check['family']}.{check['name']}"
+        metrics[f"{base}.allclose"] = check.get("allclose") is True
+        for key in ("max_abs", "mean_abs", "rms_abs"):
+            number = finite_number(check.get(key))
+            if number is not None:
+                metrics[f"{base}.{key}"] = number
+    for check in field_checks(report):
+        base = f"field.{check['mode']}.{check['field']}"
+        metrics[f"{base}.allclose"] = check.get("allclose") is True
+        for key in ("max_abs", "mean_abs", "rms_abs"):
+            number = finite_number(check.get(key))
+            if number is not None:
+                metrics[f"{base}.{key}"] = number
+    preconditioner = first_nested_key(report, "preconditioner_diagnostic")
+    if preconditioner is not None:
+        _, data = preconditioner
+        for key in (
+            "gpu_dilu_vs_reference_delta_l2",
+            "diagonal_vs_reference_delta_l2",
+        ):
+            number = finite_number(data.get(key))
+            if number is not None:
+                metrics[f"preconditioner.{key}"] = number
+        for path, metric_name in (
+            ("gpu_dilu_vs_openfoam_reference.max_abs", "preconditioner.gpu_dilu_vs_reference.max_abs"),
+            ("gpu_dilu_vs_openfoam_reference.rms_abs", "preconditioner.gpu_dilu_vs_reference.rms_abs"),
+            ("diagonal_vs_openfoam_reference.max_abs", "preconditioner.diagonal_vs_reference.max_abs"),
+            ("diagonal_vs_openfoam_reference.rms_abs", "preconditioner.diagonal_vs_reference.rms_abs"),
+        ):
+            number = finite_number(get_path(data, path))
+            if number is not None:
+                metrics[metric_name] = number
+    return metrics
+
+
+def classify_delta(current: Mapping[str, Any], baseline: Mapping[str, Any] | None) -> list[dict[str, Any]]:
+    if not baseline:
+        return [{"metric": name, "status": "newly_available", "current": value, "baseline": None} for name, value in sorted(current.items())]
+    out: list[dict[str, Any]] = []
+    previous = baseline.get("metrics", {}) if isinstance(baseline.get("metrics"), Mapping) else {}
+    for name, value in sorted(current.items()):
+        old = previous.get(name)
+        status = "newly_available"
+        if old is not None:
+            if isinstance(value, bool) and isinstance(old, bool):
+                if value == old:
+                    status = "unchanged"
+                elif value and not old:
+                    status = "improved"
+                else:
+                    status = "regressed"
+            else:
+                now_num = finite_number(value)
+                old_num = finite_number(old)
+                if now_num is not None and old_num is not None:
+                    tolerance = max(1e-15, abs(old_num) * 1e-9)
+                    if abs(now_num - old_num) <= tolerance:
+                        status = "unchanged"
+                    elif now_num < old_num:
+                        status = "improved"
+                    else:
+                        status = "regressed"
+                else:
+                    status = "unchanged" if value == old else "changed"
+        out.append({"metric": name, "status": status, "current": value, "baseline": old})
+    return out
+
+
+def metric_priority(item: Mapping[str, Any]) -> tuple[int, str]:
+    status_order = {"regressed": 0, "improved": 1, "newly_available": 2, "changed": 3, "unchanged": 4}
+    return status_order.get(str(item.get("status")), 9), str(item.get("metric"))
+
+
+def render_location(location: Any) -> str:
+    if not isinstance(location, Mapping):
+        return "-"
+    entity = location.get("entity_kind") or "array"
+    index = location.get("entity_index")
+    component = location.get("component_index")
+    if component is None:
+        return f"{entity}[{index}]"
+    return f"{entity}[{index}] component={component}"
+
+
+def render_context(report_path: Path | None, report: Mapping[str, Any] | None, baseline: Mapping[str, Any] | None) -> str:
+    generated = datetime.now(timezone.utc).isoformat()
+    if report is None or report_path is None:
+        return "\n".join(
+            [
+                "# GPU RANS Loop Diagnostic Context",
+                "",
+                f"Generated: {generated}",
+                "",
+                "No verifier report was found under repo tmp/ or /tmp GPU RANS work directories.",
+                "Next action: run `scripts/verify_gpu_rans_solver.sh --work <work> --report <work>/report.json` or the focused verifier, then rerun this hook.",
+                "",
+            ]
+        )
+
+    first = report.get("first_divergence_summary") if isinstance(report.get("first_divergence_summary"), Mapping) else {}
+    artifacts = report.get("intermediate_artifacts", {}) if isinstance(report.get("intermediate_artifacts"), Mapping) else {}
+    families = artifacts.get("families", []) if isinstance(artifacts.get("families"), list) else []
+    metrics = extract_metrics(report)
+    deltas = classify_delta(metrics, baseline)
+    checks = artifact_checks(report)
+    fields = field_checks(report)
+    preconditioner = first_nested_key(report, "preconditioner_diagnostic")
+    solver_trace = first_nested_key(report, "linear_solver_trace")
+
+    lines = [
+        "# GPU RANS Loop Diagnostic Context",
+        "",
+        f"Generated: {generated}",
+        f"Latest report: `{report_path}`",
+        f"Report status: `{report.get('status')}`",
+        f"Verifier evidence passed: `{get_path(report, 'verifier_evidence.passed')}`",
+        "",
+        "## First divergence",
+        "",
+        f"- Target: `{first.get('first_target') if first else None}`",
+        f"- Artifact family: `{first.get('artifact_family') if first else None}`",
+        f"- Stage group: `{first.get('stage_group') if first else None}`",
+        f"- Evidence path: `{first.get('evidence_path') if first else None}`",
+        f"- Field/reason: `{first.get('field') if first else None}` / `{first.get('reason') if first else None}`",
+        "",
+        "## Solver phase evidence",
+        "",
+        "| Family | Status | Why |",
+        "|---|---:|---|",
+    ]
+    for family in families:
+        if not isinstance(family, Mapping):
+            continue
+        lines.append(f"| {family.get('name')} | {family.get('status')} | {family.get('why')} |")
+
+    lines.extend(["", "## Numeric artifact comparisons", "", "| Family | Check | Status | max_abs | mean_abs | rms_abs | Location |", "|---|---|---:|---:|---:|---:|---|"])
+    for check in checks[:16]:
+        lines.append(
+            "| {family} | {name} | {status} | {max_abs} | {mean_abs} | {rms_abs} | {location} |".format(
+                family=check.get("family"),
+                name=check.get("name"),
+                status=status_word(check.get("allclose")),
+                max_abs=fmt(check.get("max_abs")),
+                mean_abs=fmt(check.get("mean_abs")),
+                rms_abs=fmt(check.get("rms_abs")),
+                location=render_location(check.get("location")),
+            )
+        )
+
+    lines.extend(["", "## Field comparison symptoms", "", "| Mode | Field | Status | max_abs | mean_abs | rms_abs | Location |", "|---|---|---:|---:|---:|---:|---|"])
+    for check in fields[:12]:
+        lines.append(
+            "| {mode} | {field} | {status} | {max_abs} | {mean_abs} | {rms_abs} | {location} |".format(
+                mode=check.get("mode"),
+                field=check.get("field"),
+                status=status_word(check.get("allclose")),
+                max_abs=fmt(check.get("max_abs")),
+                mean_abs=fmt(check.get("mean_abs")),
+                rms_abs=fmt(check.get("rms_abs")),
+                location=render_location(check.get("location")),
+            )
+        )
+
+    lines.extend(["", "## Linear solver and preconditioner trace", ""])
+    if preconditioner is None:
+        lines.append("No preconditioner diagnostic artifact found in the latest report.")
+    else:
+        path, data = preconditioner
+        lines.extend(
+            [
+                f"Preconditioner evidence path: `{path}`",
+                f"- GPU DILU vs OpenFOAM reference max_abs: `{fmt(get_path(data, 'gpu_dilu_vs_openfoam_reference.max_abs'))}`",
+                f"- GPU DILU vs OpenFOAM reference rms_abs: `{fmt(get_path(data, 'gpu_dilu_vs_openfoam_reference.rms_abs'))}`",
+                f"- Diagonal/current vs OpenFOAM reference max_abs: `{fmt(get_path(data, 'diagonal_vs_openfoam_reference.max_abs'))}`",
+                f"- Residual entering preconditioner recorded: `{data.get('residual_entering_preconditioner') is not None}`",
+            ]
+        )
+    if solver_trace is not None:
+        path, data = solver_trace
+        lines.append(f"Solver trace path: `{path}`")
+        for point in data.get("trace_points", []) if isinstance(data.get("trace_points"), list) else []:
+            if isinstance(point, Mapping):
+                lines.append(f"- `{point.get('name')}`: {point.get('step')}")
+
+    lines.extend(["", "## Delta versus retained baseline", "", "| Metric | Status | Current | Baseline |", "|---|---:|---:|---:|"])
+    for item in sorted(deltas, key=metric_priority)[:24]:
+        lines.append(f"| `{item.get('metric')}` | {item.get('status')} | {fmt(item.get('current'))} | {fmt(item.get('baseline'))} |")
+
+    lines.extend(
+        [
+            "",
+            "## Next target hint",
+            "",
+            f"Focus first on `{first.get('first_target') if first else 'unknown'}`. Treat downstream field symptoms as unreliable until that artifact or missing evidence closes.",
+            "",
+        ]
+    )
+    return "\n".join(lines)
+
+
+def parse_args(argv: list[str] | None = None) -> argparse.Namespace:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--root", type=Path, default=ROOT)
+    parser.add_argument("--report", type=Path, default=None, help="Explicit report path; otherwise discover newest verifier report")
+    parser.add_argument("--context", type=Path, default=DEFAULT_CONTEXT)
+    parser.add_argument("--baseline", type=Path, default=DEFAULT_BASELINE)
+    return parser.parse_args(argv)
+
+
+def main(argv: list[str] | None = None) -> int:
+    args = parse_args(argv)
+    root = args.root.resolve()
+    if args.report is not None:
+        report_path = args.report.resolve()
+        loaded = load_json(report_path)
+        report = loaded if is_verifier_report(loaded) else None
+    else:
+        report_path, report = discover_latest_report(root)
+
+    baseline = load_json(args.baseline) if args.baseline.exists() else None
+    baseline_mapping = baseline if isinstance(baseline, Mapping) else None
+    args.context.parent.mkdir(parents=True, exist_ok=True)
+    args.context.write_text(render_context(report_path, report, baseline_mapping), encoding="utf-8")
+
+    if report is not None and report_path is not None:
+        current = {
+            "updated_at": datetime.now(timezone.utc).isoformat(),
+            "report": str(report_path),
+            "metrics": extract_metrics(report),
+            "first_divergence_summary": report.get("first_divergence_summary"),
+        }
+        args.baseline.parent.mkdir(parents=True, exist_ok=True)
+        args.baseline.write_text(json.dumps(current, indent=2, sort_keys=True) + "\n", encoding="utf-8")
+        print(f"updated diagnostic context: {args.context} from {report_path}")
+    else:
+        print(f"updated diagnostic context without report: {args.context}")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())

scripts/verify_gpu_rans_solver.sh

@@ -96,17 +96,23 @@ def require(condition: bool, message: str) -> None:
 
 backend = report.get("backend") if isinstance(report.get("backend"), dict) else {}
 verifier_evidence = report.get("verifier_evidence") if isinstance(report.get("verifier_evidence"), dict) else {}
+backend_trust = ((verifier_evidence.get("mode_criteria") or {}).get("backend_trust") or {}) if isinstance(verifier_evidence.get("mode_criteria"), dict) else {}
 modes = report.get("modes") if isinstance(report.get("modes"), dict) else {}
 provider = str(backend.get("provider") or "").lower()
 device = str(backend.get("device") or "").lower()
-
+capabilities = set(backend.get("capabilities") or [])
+primitive = backend.get("primitive_evidence") if isinstance(backend.get("primitive_evidence"), dict) else {}
 require(verify_status == 0, f"underlying AirfRANS verifier exited {verify_status}")
 require(report.get("status") == "passed", f"report.status is {report.get('status')!r}, expected 'passed'")
 require(backend.get("requested") == "gpu", f"backend.requested is {backend.get('requested')!r}, expected 'gpu'")
 require(backend.get("selected") == "gpu", f"backend.selected is {backend.get('selected')!r}, expected 'gpu'")
 require("cpu" not in provider and provider not in {"foam_stepper_cpu", "openfoam"}, f"backend.provider looks CPU-backed: {backend.get('provider')!r}")
 require(device not in {"host", "cpu"}, f"backend.device looks CPU-backed: {backend.get('device')!r}")
-require(backend.get("counts_as_gpu_algorithm_progress") is not False, "backend explicitly says it does not count as GPU progress")
+require(backend.get("used_cpu_fallback") is False, f"backend.used_cpu_fallback is {backend.get('used_cpu_fallback')!r}, expected false")
+require((backend.get("cpu_fallback") or {}).get("allowed") is False, "backend.cpu_fallback.allowed is not false")
+require("no_cpu_fallback" in capabilities, "backend capabilities do not include no_cpu_fallback")
+require(primitive.get("counts_as_full_gpu_rans_solver") is False, "primitive GPU evidence is not explicitly marked non-acceptance")
+require(backend_trust.get("status") == "trusted", f"backend trust status is {backend_trust.get('status')!r}, expected 'trusted'")
 require(verifier_evidence.get("passed") is True, "verifier_evidence.passed is not true")
 
 for mode_name in required_modes:
@@ -116,6 +122,12 @@ for mode_name in required_modes:
     require(mode_backend.get("selected") == "gpu", f"modes.{mode_name}.backend.selected is not 'gpu'")
     execution_path = str(mode.get("execution_path") or "").lower()
     require("gpu" in execution_path, f"modes.{mode_name}.execution_path does not identify a GPU path: {mode.get('execution_path')!r}")
+    require(mode_backend.get("used_cpu_fallback") is False, f"modes.{mode_name}.backend.used_cpu_fallback is not false")
+    gpu_solver = mode.get("gpu_solver") if isinstance(mode.get("gpu_solver"), dict) else {}
+    gpu_solver_backend = gpu_solver.get("backend") if isinstance(gpu_solver.get("backend"), dict) else {}
+    require(gpu_solver.get("status") == "executed", f"modes.{mode_name}.gpu_solver.status is not executed")
+    require(gpu_solver_backend.get("selected") == "gpu", f"modes.{mode_name}.gpu_solver.backend.selected is not 'gpu'")
+    require(gpu_solver_backend.get("used_cpu_fallback") is False, f"modes.{mode_name}.gpu_solver.backend.used_cpu_fallback is not false")
     comparisons = mode.get("comparisons") if isinstance(mode.get("comparisons"), dict) else {}
     for field in required_fields:
         comparison = comparisons.get(field) if isinstance(comparisons.get(field), dict) else {}