From c2b0c180bbf3f95c6d8f9e7a9b97048839052045 Mon Sep 17 00:00:00 2001
From: Tyagi <gthyagi@gmail.com>
Date: Tue, 16 Jun 2026 07:39:53 +0530
Subject: [PATCH 01/11] Fix spherical shell internal boundary labels

Rework SphericalShellInternalBoundary mesh construction so the shell volume is retained, the internal spherical surface is embedded, and duplicate internal geometry is removed before meshing.

The previous nested-sphere fragment path could leave the Lower boundary unusable for BdIntegral; the benchmark probe observed lower_area=0.0 even though the Stokes solve converged.

Add a boundary-integral regression test that checks nonzero, close-to-analytic Lower, Internal, and Upper surface areas.

Validation: ./uw build passed. py_compile passed for src/underworld3/meshing/spherical.py. pytest -q tests/test_0502_boundary_integrals.py::test_bd_integral_spherical_internal_boundary_areas passed. 005_internal_boundary_delta_probe.py passed with -uw_mesh_source uw3_builtin for serial Nitsche, serial constrained, 8-rank Nitsche, and 8-rank constrained runs.
---
 src/underworld3/meshing/spherical.py  | 151 +++++++++++++-------------
 tests/test_0502_boundary_integrals.py |  46 ++++++++
 2 files changed, 122 insertions(+), 75 deletions(-)

diff --git a/src/underworld3/meshing/spherical.py b/src/underworld3/meshing/spherical.py
index fb03e1b8d..4e7288d95 100644
--- a/src/underworld3/meshing/spherical.py
+++ b/src/underworld3/meshing/spherical.py
@@ -580,96 +580,97 @@ class regions(Enum):
         gmsh.option.setNumber("General.Verbosity", gmsh_verbosity)
         gmsh.model.add("SphereShell_with_Internal_Surface")
 
-        # Create three concentric spheres and use OCC fragment to split
-        # into two non-overlapping shell volumes sharing the internal surface
-        ball_outer = gmsh.model.occ.addSphere(0, 0, 0, radiusOuter)
-        ball_internal = gmsh.model.occ.addSphere(0, 0, 0, radiusInternal)
-        ball_inner = gmsh.model.occ.addSphere(0, 0, 0, radiusInner)
-
-        # Fragment creates non-overlapping pieces from the boolean intersection
-        out_dimtags, out_map = gmsh.model.occ.fragment(
-            [(3, ball_outer)],
-            [(3, ball_internal), (3, ball_inner)],
+        # Create the spherical shell volume.
+        outer = gmsh.model.occ.addSphere(0.0, 0.0, 0.0, radiusOuter)
+        inner = gmsh.model.occ.addSphere(0.0, 0.0, 0.0, radiusInner)
+        gmsh.model.occ.cut(
+            [(3, outer)],
+            [(3, inner)],
+            removeObject=True,
+            removeTool=True,
+        )
+
+        # Create an internal shell only to obtain a clean spherical surface at
+        # radiusInternal. That surface is embedded into the shell volume below;
+        # the duplicate volume and duplicate lower surface are removed before
+        # meshing.
+        internal = gmsh.model.occ.addSphere(0.0, 0.0, 0.0, radiusInternal)
+        inner_copy = gmsh.model.occ.addSphere(0.0, 0.0, 0.0, radiusInner)
+        gmsh.model.occ.cut(
+            [(3, internal)],
+            [(3, inner_copy)],
+            removeObject=True,
+            removeTool=True,
         )
 
         gmsh.model.occ.synchronize()
         gmsh.option.setNumber("Mesh.CharacteristicLengthMax", cellSize)
 
-        # Identify volumes and surfaces by bounding box
-        # For a sphere, bbox diagonal = sqrt(3) * radius
-        volumes = gmsh.model.getEntities(3)
-        surfaces = gmsh.model.getEntities(2)
-
         def bbox_radius(dimtag):
-            """Estimate the sphere radius from a bounding box diagonal."""
+            """Estimate a concentric sphere radius from the bounding box."""
             bb = gmsh.model.get_bounding_box(dimtag[0], dimtag[1])
             return np.sqrt(bb[3]**2 + bb[4]**2 + bb[5]**2) / np.sqrt(3.0)
 
-        inner_vols = []
-        outer_vols = []
-        solid_ball_vols = []  # r < radiusInner — to be removed
-
-        for vol in volumes:
-            r_est = bbox_radius(vol)
-            if np.isclose(r_est, radiusInner, atol=cellSize):
-                solid_ball_vols.append(vol)
-            elif np.isclose(r_est, radiusInternal, atol=cellSize):
-                inner_vols.append(vol)
-            elif np.isclose(r_est, radiusOuter, atol=cellSize):
-                outer_vols.append(vol)
-
-        # Remove the solid inner ball (r < radiusInner)
-        if solid_ball_vols:
-            gmsh.model.occ.remove(solid_ball_vols, recursive=True)
-            gmsh.model.occ.synchronize()
-
-        # Re-query after removal
         volumes = gmsh.model.getEntities(3)
-        surfaces = gmsh.model.getEntities(2)
+        shell_vols = [vol for vol in volumes if np.isclose(bbox_radius(vol), radiusOuter, atol=cellSize * 0.5)]
+        duplicate_vols = [vol for vol in volumes if np.isclose(bbox_radius(vol), radiusInternal, atol=cellSize * 0.5)]
+
+        if len(shell_vols) != 1:
+            raise RuntimeError(
+                "Could not identify the spherical-shell volume while building "
+                "SphericalShellInternalBoundary."
+            )
 
-        # Classify surfaces by bounding box radius
-        for surface in surfaces:
+        shell_vol = shell_vols[0]
+        shell_boundary = {
+            dimtag[1]
+            for dimtag in gmsh.model.getBoundary([shell_vol], oriented=False, recursive=False)
+            if dimtag[0] == 2
+        }
+
+        lower_surface_tags = []
+        internal_surface_tags = []
+        upper_surface_tags = []
+        duplicate_lower_tags = []
+
+        for surface in gmsh.model.getEntities(2):
+            surface_tag = surface[1]
             r_est = bbox_radius(surface)
             if np.isclose(r_est, radiusInner, atol=cellSize * 0.5):
-                gmsh.model.addPhysicalGroup(
-                    surface[0], [surface[1]],
-                    boundaries.Lower.value, name=boundaries.Lower.name,
-                )
+                if surface_tag in shell_boundary:
+                    lower_surface_tags.append(surface_tag)
+                else:
+                    duplicate_lower_tags.append(surface_tag)
             elif np.isclose(r_est, radiusOuter, atol=cellSize * 0.5):
-                gmsh.model.addPhysicalGroup(
-                    surface[0], [surface[1]],
-                    boundaries.Upper.value, name=boundaries.Upper.name,
-                )
+                upper_surface_tags.append(surface_tag)
             elif np.isclose(r_est, radiusInternal, atol=cellSize * 0.5):
-                gmsh.model.addPhysicalGroup(
-                    surface[0], [surface[1]],
-                    boundaries.Internal.value, name=boundaries.Internal.name,
-                )
-
-        # Classify remaining volumes into Inner and Outer
-        inner_vol_tags = [v[1] for v in inner_vols if v not in solid_ball_vols]
-        outer_vol_tags = [v[1] for v in outer_vols]
-        # Re-classify from current volumes in case tags changed after removal
-        inner_vol_tags = []
-        outer_vol_tags = []
-        for vol in volumes:
-            r_est = bbox_radius(vol)
-            if r_est < radiusInternal + cellSize * 0.5:
-                inner_vol_tags.append(vol[1])
-            else:
-                outer_vol_tags.append(vol[1])
-
-        # Region physical groups
-        if inner_vol_tags:
-            gmsh.model.addPhysicalGroup(3, inner_vol_tags,
-                                        regions.Inner.value, name=regions.Inner.name)
-        if outer_vol_tags:
-            gmsh.model.addPhysicalGroup(3, outer_vol_tags,
-                                        regions.Outer.value, name=regions.Outer.name)
-
-        # Combined elements group
-        all_vol_tags = inner_vol_tags + outer_vol_tags
-        gmsh.model.addPhysicalGroup(3, all_vol_tags, 99999, "Elements")
+                internal_surface_tags.append(surface_tag)
+
+        if not lower_surface_tags or not upper_surface_tags or not internal_surface_tags:
+            raise RuntimeError(
+                "Could not identify Lower, Internal, and Upper spherical surfaces "
+                "while building SphericalShellInternalBoundary."
+            )
+
+        gmsh.model.mesh.embed(2, internal_surface_tags, shell_vol[0], shell_vol[1])
+
+        remove_dimtags = duplicate_vols + [(2, tag) for tag in duplicate_lower_tags]
+        if remove_dimtags:
+            gmsh.model.remove_entities(remove_dimtags, recursive=False)
+            gmsh.model.occ.remove(remove_dimtags, recursive=False)
+            gmsh.model.occ.synchronize()
+
+        gmsh.model.addPhysicalGroup(
+            2, lower_surface_tags, boundaries.Lower.value, name=boundaries.Lower.name
+        )
+        gmsh.model.addPhysicalGroup(
+            2, internal_surface_tags, boundaries.Internal.value, name=boundaries.Internal.name
+        )
+        gmsh.model.addPhysicalGroup(
+            2, upper_surface_tags, boundaries.Upper.value, name=boundaries.Upper.name
+        )
+
+        gmsh.model.addPhysicalGroup(shell_vol[0], [shell_vol[1]], 99999, "Elements")
 
         gmsh.model.mesh.generate(3)
         gmsh.write(uw_filename)
diff --git a/tests/test_0502_boundary_integrals.py b/tests/test_0502_boundary_integrals.py
index 0cd9e18d9..1721db7db 100644
--- a/tests/test_0502_boundary_integrals.py
+++ b/tests/test_0502_boundary_integrals.py
@@ -293,6 +293,52 @@ def test_bd_integral_annulus_internal_normal_tangential():
     assert abs(value) < 0.05, f"Expected ~0, got {value}"
 
 
+# --- Spherical shell internal boundary tests ---
+
+from underworld3.meshing import SphericalShellInternalBoundary
+
+_R_SHELL_INNER = 0.55
+_R_SHELL_INTERNAL = 0.775
+_R_SHELL_OUTER = 1.0
+_mesh_spherical_internal = None
+
+
+def _get_spherical_internal_mesh():
+    global _mesh_spherical_internal
+    if _mesh_spherical_internal is None:
+        _mesh_spherical_internal = SphericalShellInternalBoundary(
+            radiusOuter=_R_SHELL_OUTER,
+            radiusInternal=_R_SHELL_INTERNAL,
+            radiusInner=_R_SHELL_INNER,
+            cellSize=0.25,
+            degree=1,
+            qdegree=2,
+        )
+        uw.discretisation.MeshVariable(
+            "T_spherical_internal", _mesh_spherical_internal, 1, degree=1
+        )
+    return _mesh_spherical_internal
+
+
+def test_bd_integral_spherical_internal_boundary_areas():
+    """SphericalShellInternalBoundary preserves Lower/Internal/Upper labels."""
+
+    mesh_spherical = _get_spherical_internal_mesh()
+    expected_areas = {
+        "Lower": 4.0 * np.pi * _R_SHELL_INNER**2,
+        "Internal": 4.0 * np.pi * _R_SHELL_INTERNAL**2,
+        "Upper": 4.0 * np.pi * _R_SHELL_OUTER**2,
+    }
+
+    for boundary, expected in expected_areas.items():
+        value = uw.maths.BdIntegral(mesh_spherical, fn=1.0, boundary=boundary).evaluate()
+        relative_error = abs(value - expected) / expected
+        assert relative_error < 0.06, (
+            f"{boundary} area should be close to {expected:.4f}; "
+            f"got {value:.4f} (relative error {relative_error:.3f})"
+        )
+
+
 def _build_spherical_shell_for_integrals():
     from underworld3.meshing import SphericalShell
 

From 55566d980d13466b2c591fc8d3e0426f3294262e Mon Sep 17 00:00:00 2001
From: Tyagi <gthyagi@gmail.com>
Date: Tue, 16 Jun 2026 08:03:39 +0530
Subject: [PATCH 02/11] Polish spherical shell internal boundary validation

Add explicit radius ordering validation for SphericalShellInternalBoundary and wrap the entity-selection comprehensions introduced by the internal-boundary label fix.

Validation: py_compile passed for src/underworld3/meshing/spherical.py and test_bd_integral_spherical_internal_boundary_areas passed.
---
 src/underworld3/meshing/spherical.py | 18 +++++++++++++++---
 1 file changed, 15 insertions(+), 3 deletions(-)

diff --git a/src/underworld3/meshing/spherical.py b/src/underworld3/meshing/spherical.py
index 4e7288d95..2ea99fa9d 100644
--- a/src/underworld3/meshing/spherical.py
+++ b/src/underworld3/meshing/spherical.py
@@ -571,9 +571,13 @@ class regions(Enum):
     else:
         uw_filename = filename
 
-    # Check if r_i is greater than 0
     if radiusInner <= 0:
         raise ValueError("The inner radius must be greater than 0.")
+    if not radiusInner < radiusInternal < radiusOuter:
+        raise ValueError(
+            "SphericalShellInternalBoundary requires "
+            "radiusInner < radiusInternal < radiusOuter."
+        )
 
     if uw.mpi.rank == 0:
         gmsh.initialize()
@@ -612,8 +616,16 @@ def bbox_radius(dimtag):
             return np.sqrt(bb[3]**2 + bb[4]**2 + bb[5]**2) / np.sqrt(3.0)
 
         volumes = gmsh.model.getEntities(3)
-        shell_vols = [vol for vol in volumes if np.isclose(bbox_radius(vol), radiusOuter, atol=cellSize * 0.5)]
-        duplicate_vols = [vol for vol in volumes if np.isclose(bbox_radius(vol), radiusInternal, atol=cellSize * 0.5)]
+        shell_vols = [
+            vol
+            for vol in volumes
+            if np.isclose(bbox_radius(vol), radiusOuter, atol=cellSize * 0.5)
+        ]
+        duplicate_vols = [
+            vol
+            for vol in volumes
+            if np.isclose(bbox_radius(vol), radiusInternal, atol=cellSize * 0.5)
+        ]
 
         if len(shell_vols) != 1:
             raise RuntimeError(

From 7f9235172750f803e9f8d132fb51d86f34b5a814 Mon Sep 17 00:00:00 2001
From: Tyagi <gthyagi@gmail.com>
Date: Wed, 17 Jun 2026 15:25:09 +0530
Subject: [PATCH 03/11] Add constrained spherical shell response regression

Record the Zhong-style SphericalShellInternalBoundary constrained free-slip behaviour in UW3 tests.

The new regression separates the constrained weak-form check from the practical solver failure: monolithic LU constrained matches monolithic Nitsche, while the fast grouped pressure-plus-multiplier Schur path remains a strict expected failure against the Zhong velocity response.

This documents that the remaining fix belongs in the practical grouped Schur/preconditioner path, not in the basic 3-D vector-scalar multiplier assembly.
---
 ...64_constrained_spherical_shell_response.py | 184 ++++++++++++++++++
 1 file changed, 184 insertions(+)
 create mode 100644 tests/test_1064_constrained_spherical_shell_response.py

diff --git a/tests/test_1064_constrained_spherical_shell_response.py b/tests/test_1064_constrained_spherical_shell_response.py
new file mode 100644
index 000000000..3071047d3
--- /dev/null
+++ b/tests/test_1064_constrained_spherical_shell_response.py
@@ -0,0 +1,184 @@
+"""3-D spherical-shell constrained free-slip response regression.
+
+This test records two facts exposed by the Zhong et al. (2008)-style benchmark:
+
+* with monolithic LU, ``Stokes_Constrained`` and Nitsche free slip solve the same
+  internal-boundary Stokes response to within a tight tolerance;
+* the practical fast grouped-Schur constrained path currently does not reproduce
+  the Zhong velocity response and remains an expected failure.
+
+Run:
+    pixi run -e amr-dev pytest -q tests/test_1064_constrained_spherical_shell_response.py
+"""
+
+from functools import cache
+
+import numpy as np
+import pytest
+import sympy
+
+import underworld3 as uw
+
+pytestmark = [pytest.mark.level_3, pytest.mark.slow, pytest.mark.tier_c]
+
+RADIUS_INNER = 0.55
+RADIUS_INTERNAL = 0.775
+RADIUS_OUTER = 1.0
+CELL_SIZE = 1.0 / 8.0
+HARMONIC_DEGREE = 2
+NITSCHE_GAMMA = 10.0
+
+ZHONG_SURFACE_VELOCITY = 1.006e-2
+ZHONG_CMB_VELOCITY = 1.186e-2
+
+
+@cache
+def solve_response(method, solver_mode):
+    mesh = uw.meshing.SphericalShellInternalBoundary(
+        radiusOuter=RADIUS_OUTER,
+        radiusInternal=RADIUS_INTERNAL,
+        radiusInner=RADIUS_INNER,
+        cellSize=CELL_SIZE,
+        qdegree=2,
+        degree=1,
+    )
+
+    velocity = uw.discretisation.MeshVariable(
+        f"U_{method}_{solver_mode}",
+        mesh,
+        mesh.dim,
+        degree=2,
+        vtype=uw.VarType.VECTOR,
+    )
+    pressure = uw.discretisation.MeshVariable(
+        f"P_{method}_{solver_mode}",
+        mesh,
+        1,
+        degree=1,
+        continuous=True,
+    )
+
+    theta = mesh.CoordinateSystem.xR[1]
+    unit_r = mesh.CoordinateSystem.unit_e_0
+    y_l0 = sympy.assoc_legendre(HARMONIC_DEGREE, 0, sympy.cos(theta))
+    harmonic_norm = 4.0 * np.pi / (2 * HARMONIC_DEGREE + 1)
+
+    if method == "constrained":
+        stokes = uw.systems.Stokes_Constrained(
+            mesh,
+            velocityField=velocity,
+            pressureField=pressure,
+        )
+    elif method == "nitsche":
+        stokes = uw.systems.Stokes(
+            mesh,
+            velocityField=velocity,
+            pressureField=pressure,
+        )
+    else:
+        raise ValueError(f"Unknown response method {method!r}")
+
+    stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
+    stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0
+    stokes.bodyforce = sympy.Matrix([0.0, 0.0, 0.0])
+    stokes.add_natural_bc(y_l0 * unit_r, mesh.boundaries.Internal.name)
+
+    if method == "nitsche":
+        stokes.add_nitsche_bc("Upper", normal=unit_r, gamma=NITSCHE_GAMMA)
+        stokes.add_nitsche_bc("Lower", normal=-unit_r, gamma=NITSCHE_GAMMA)
+    else:
+        stokes.add_constraint_bc(
+            "Upper",
+            g=0.0,
+            normal=unit_r,
+            augmentation_base=1.0e4,
+            degree=2,
+        )
+        stokes.add_constraint_bc(
+            "Lower",
+            g=0.0,
+            normal=-unit_r,
+            augmentation_base=1.0e4,
+            degree=2,
+        )
+
+    stokes.petsc_use_nullspace = True
+    stokes.tolerance = 1.0e-7
+    stokes.petsc_options["snes_type"] = "ksponly"
+
+    if solver_mode == "monolithic":
+        stokes.petsc_options["ksp_type"] = "preonly"
+        stokes.petsc_options["pc_type"] = "lu"
+        stokes.petsc_options["pc_factor_mat_solver_type"] = "mumps"
+        stokes.petsc_options["pc_use_amat"] = None
+    elif solver_mode == "fast_schur":
+        if method != "constrained":
+            raise ValueError("fast_schur mode is only defined for constrained runs")
+        stokes.petsc_options["pc_fieldsplit_schur_precondition"] = "selfp"
+        stokes.petsc_options["fieldsplit_1_ksp_type"] = "preonly"
+        stokes.petsc_options["fieldsplit_1_pc_type"] = "gasm"
+    else:
+        raise ValueError(f"Unknown solver mode {solver_mode!r}")
+
+    stokes.solve()
+
+    horizontal_v2 = velocity.sym.dot(velocity.sym) - velocity.sym.dot(unit_r) ** 2
+
+    surface_velocity = np.sqrt(
+        uw.maths.BdIntegral(mesh, horizontal_v2, boundary="Upper").evaluate()
+        / (
+            RADIUS_OUTER**2
+            * HARMONIC_DEGREE
+            * (HARMONIC_DEGREE + 1)
+            * harmonic_norm
+        )
+    )
+    cmb_velocity = np.sqrt(
+        uw.maths.BdIntegral(mesh, horizontal_v2, boundary="Lower").evaluate()
+        / (
+            RADIUS_INNER**2
+            * HARMONIC_DEGREE
+            * (HARMONIC_DEGREE + 1)
+            * harmonic_norm
+        )
+    )
+
+    return (
+        float(surface_velocity),
+        float(cmb_velocity),
+        int(stokes.snes.getConvergedReason()),
+    )
+
+
+def test_monolithic_constrained_matches_monolithic_nitsche_response():
+    nitsche_surface, nitsche_cmb, nitsche_reason = solve_response(
+        "nitsche",
+        "monolithic",
+    )
+    constrained_surface, constrained_cmb, constrained_reason = solve_response(
+        "constrained",
+        "monolithic",
+    )
+
+    assert nitsche_reason > 0
+    assert constrained_reason > 0
+    assert abs(constrained_surface - nitsche_surface) / nitsche_surface < 0.01
+    assert abs(constrained_cmb - nitsche_cmb) / nitsche_cmb < 0.01
+
+
+@pytest.mark.xfail(
+    reason=(
+        "Known fast grouped-Schur constrained response failure for the "
+        "3-D SphericalShellInternalBoundary Zhong-style load."
+    ),
+    strict=True,
+)
+def test_fast_schur_constrained_matches_zhong_velocity_response():
+    surface_velocity, cmb_velocity, snes_reason = solve_response(
+        "constrained",
+        "fast_schur",
+    )
+
+    assert snes_reason > 0
+    assert abs(surface_velocity - ZHONG_SURFACE_VELOCITY) / ZHONG_SURFACE_VELOCITY < 0.05
+    assert abs(cmb_velocity - ZHONG_CMB_VELOCITY) / ZHONG_CMB_VELOCITY < 0.05

From c086517ba556fae94cc43c5faecb468273d938ac Mon Sep 17 00:00:00 2001
From: Tyagi <gthyagi@gmail.com>
Date: Thu, 18 Jun 2026 00:01:19 +0530
Subject: [PATCH 04/11] Add exact constrained fieldsplit regression

Record that the constrained weak form matches monolithic Nitsche while PETSc field-split constrained solves remain incorrect even with exact LU sub-solves on the velocity and grouped pressure-multiplier blocks.

Validation: pixi run -e amr-dev pytest -q tests/test_1064_constrained_spherical_shell_response.py -> 1 passed, 2 xfailed in 137.51s.
---
 ...64_constrained_spherical_shell_response.py | 33 +++++++++++++++++++
 1 file changed, 33 insertions(+)

diff --git a/tests/test_1064_constrained_spherical_shell_response.py b/tests/test_1064_constrained_spherical_shell_response.py
index 3071047d3..36040743b 100644
--- a/tests/test_1064_constrained_spherical_shell_response.py
+++ b/tests/test_1064_constrained_spherical_shell_response.py
@@ -111,6 +111,16 @@ def solve_response(method, solver_mode):
         stokes.petsc_options["pc_type"] = "lu"
         stokes.petsc_options["pc_factor_mat_solver_type"] = "mumps"
         stokes.petsc_options["pc_use_amat"] = None
+    elif solver_mode == "fieldsplit_exact":
+        if method != "constrained":
+            raise ValueError(
+                "fieldsplit_exact mode is only defined for constrained runs"
+            )
+        stokes.petsc_options["pc_fieldsplit_schur_precondition"] = "selfp"
+        stokes.petsc_options["fieldsplit_velocity_ksp_type"] = "preonly"
+        stokes.petsc_options["fieldsplit_velocity_pc_type"] = "lu"
+        stokes.petsc_options["fieldsplit_1_ksp_type"] = "preonly"
+        stokes.petsc_options["fieldsplit_1_pc_type"] = "lu"
     elif solver_mode == "fast_schur":
         if method != "constrained":
             raise ValueError("fast_schur mode is only defined for constrained runs")
@@ -166,6 +176,29 @@ def test_monolithic_constrained_matches_monolithic_nitsche_response():
     assert abs(constrained_cmb - nitsche_cmb) / nitsche_cmb < 0.01
 
 
+@pytest.mark.xfail(
+    reason=(
+        "Known constrained field-split algebra failure: exact LU sub-solves in "
+        "the velocity | [p,h] split still do not reproduce monolithic LU."
+    ),
+    strict=True,
+)
+def test_exact_fieldsplit_constrained_matches_monolithic_nitsche_response():
+    nitsche_surface, nitsche_cmb, nitsche_reason = solve_response(
+        "nitsche",
+        "monolithic",
+    )
+    constrained_surface, constrained_cmb, constrained_reason = solve_response(
+        "constrained",
+        "fieldsplit_exact",
+    )
+
+    assert nitsche_reason > 0
+    assert constrained_reason > 0
+    assert abs(constrained_surface - nitsche_surface) / nitsche_surface < 0.01
+    assert abs(constrained_cmb - nitsche_cmb) / nitsche_cmb < 0.01
+
+
 @pytest.mark.xfail(
     reason=(
         "Known fast grouped-Schur constrained response failure for the "

From 2ed383fe55e8dd8a6b087acfd644d4045011963b Mon Sep 17 00:00:00 2001
From: Tyagi <gthyagi@gmail.com>
Date: Thu, 18 Jun 2026 00:24:01 +0530
Subject: [PATCH 05/11] Correct constrained response regression reference

Use the validated Nitsche/default field-split response as the Zhong reference and relabel the direct-LU Nitsche/constrained comparison as a diagnostic path rather than a monolithic benchmark reference.

Validation: FI_PROVIDER=tcp pixi run -e amr-dev pytest -q tests/test_1064_constrained_spherical_shell_response.py -> 2 passed, 2 xfailed in 174.07s.
---
 ...64_constrained_spherical_shell_response.py | 33 ++++++++++++++-----
 1 file changed, 25 insertions(+), 8 deletions(-)

diff --git a/tests/test_1064_constrained_spherical_shell_response.py b/tests/test_1064_constrained_spherical_shell_response.py
index 36040743b..56acda5ad 100644
--- a/tests/test_1064_constrained_spherical_shell_response.py
+++ b/tests/test_1064_constrained_spherical_shell_response.py
@@ -1,9 +1,12 @@
 """3-D spherical-shell constrained free-slip response regression.
 
-This test records two facts exposed by the Zhong et al. (2008)-style benchmark:
+This test records facts exposed by the Zhong et al. (2008)-style benchmark:
 
-* with monolithic LU, ``Stokes_Constrained`` and Nitsche free slip solve the same
-  internal-boundary Stokes response to within a tight tolerance;
+* the validated Nitsche/default field-split path reproduces the Zhong velocity
+  scale for this low-resolution response case;
+* a direct-LU diagnostic path gives matching Nitsche and constrained responses,
+  but it does not reproduce the validated Nitsche/default response and should
+  not be treated as the benchmark reference;
 * the practical fast grouped-Schur constrained path currently does not reproduce
   the Zhong velocity response and remains an expected failure.
 
@@ -111,6 +114,8 @@ def solve_response(method, solver_mode):
         stokes.petsc_options["pc_type"] = "lu"
         stokes.petsc_options["pc_factor_mat_solver_type"] = "mumps"
         stokes.petsc_options["pc_use_amat"] = None
+    elif solver_mode == "default":
+        pass
     elif solver_mode == "fieldsplit_exact":
         if method != "constrained":
             raise ValueError(
@@ -160,7 +165,18 @@ def solve_response(method, solver_mode):
     )
 
 
-def test_monolithic_constrained_matches_monolithic_nitsche_response():
+def test_default_nitsche_matches_zhong_velocity_response():
+    surface_velocity, cmb_velocity, snes_reason = solve_response(
+        "nitsche",
+        "default",
+    )
+
+    assert snes_reason > 0
+    assert abs(surface_velocity - ZHONG_SURFACE_VELOCITY) / ZHONG_SURFACE_VELOCITY < 0.05
+    assert abs(cmb_velocity - ZHONG_CMB_VELOCITY) / ZHONG_CMB_VELOCITY < 0.05
+
+
+def test_direct_lu_diagnostic_constrained_matches_direct_lu_diagnostic_nitsche():
     nitsche_surface, nitsche_cmb, nitsche_reason = solve_response(
         "nitsche",
         "monolithic",
@@ -178,15 +194,16 @@ def test_monolithic_constrained_matches_monolithic_nitsche_response():
 
 @pytest.mark.xfail(
     reason=(
-        "Known constrained field-split algebra failure: exact LU sub-solves in "
-        "the velocity | [p,h] split still do not reproduce monolithic LU."
+        "Known constrained field-split failure: LU sub-solves in the "
+        "velocity | [p,h] preconditioner still do not reproduce the validated "
+        "Nitsche/default velocity response."
     ),
     strict=True,
 )
-def test_exact_fieldsplit_constrained_matches_monolithic_nitsche_response():
+def test_lu_subsolve_fieldsplit_constrained_matches_default_nitsche_response():
     nitsche_surface, nitsche_cmb, nitsche_reason = solve_response(
         "nitsche",
-        "monolithic",
+        "default",
     )
     constrained_surface, constrained_cmb, constrained_reason = solve_response(
         "constrained",

From 815e0368e494df09a64c72c7ed94142260726fa9 Mon Sep 17 00:00:00 2001
From: Tyagi <gthyagi@gmail.com>
Date: Thu, 18 Jun 2026 09:24:17 +0530
Subject: [PATCH 06/11] Add Stokes FEM residual diagnostic hook

Expose a low-level Stokes saddle-point helper for diagnostic residual recovery.

The helper assembles the volume-only PETSc FEM residual into local field layouts, with an optional cell-index path through DMPlexComputeResidualByKey for boundary-strip probes. This is intended for CBF/topography debugging and does not alter the normal Stokes solve path.

Validated by rebuilding the mantle-convection worktree with ./uw build and rerunning the Zhong 010 serial and 8-rank benchmark paths from the benchmark repository.
---
 src/underworld3/cython/petsc_extras.pxi       |   9 ++
 .../cython/petsc_generic_snes_solvers.pyx     | 121 ++++++++++++++++++
 2 files changed, 130 insertions(+)

diff --git a/src/underworld3/cython/petsc_extras.pxi b/src/underworld3/cython/petsc_extras.pxi
index 23785c99c..941207e5e 100644
--- a/src/underworld3/cython/petsc_extras.pxi
+++ b/src/underworld3/cython/petsc_extras.pxi
@@ -31,6 +31,13 @@ cdef CHKERRQ(PetscErrorCode ierr):
     cdef int interr = <int>ierr
     if ierr != 0: raise RuntimeError(f"PETSc error code '{interr}' was encountered.\nhttps://www.mcs.anl.gov/petsc/petsc-current/include/petscerror.h.html")
 
+cdef extern from "petscdstypes.h":
+    ctypedef struct PetscFormKey:
+        PetscDMLabel label
+        PetscInt value
+        PetscInt field
+        PetscInt part
+
 cdef extern from "petsc_compat.h":
 
     PetscErrorCode PetscDSAddBoundary_UW( PetscDM, DMBoundaryConditionType, const char[], const char[] , PetscInt, PetscInt, PetscInt *, void (*)(), void (*)(), PetscInt, const PetscInt *, void *)
@@ -52,6 +59,8 @@ cdef extern from "petsc_compat.h":
 cdef extern from "petsc.h" nogil:
     PetscErrorCode PetscDSSetConstants(PetscDS, PetscInt, const PetscScalar[])
     PetscErrorCode DMPlexSNESComputeBoundaryFEM( PetscDM, void *, void *)
+    PetscErrorCode DMPlexSNESComputeResidualFEM( PetscDM, PetscVec, PetscVec, void *)
+    PetscErrorCode DMPlexComputeResidualByKey( PetscDM, PetscFormKey, PetscIS, PetscReal, PetscVec, PetscVec, PetscReal, PetscVec, void *)
     # PetscErrorCode DMPlexSetSNESLocalFEM( PetscDM, void *, void *, void *)
     # PetscErrorCode DMPlexSetSNESLocalFEM( PetscDM, PetscBool, void *)
     PetscErrorCode DMPlexComputeGeometryFVM( PetscDM dm, PetscVec *cellgeom, PetscVec *facegeom)
diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
index d64b3c122..9eae9ee9d 100644
--- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx
+++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
@@ -6582,6 +6582,127 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         if uw.mpi.rank == 0 and self.verbose:
             print(f"Region DS: inactive region '{label_name}' gets trivial DS", flush=True)
 
+    def compute_volume_residual_fields(self, time=None, verbose=False, cell_indices=None, residual_field_id=None):
+        """Return the volume-only FEM residual in each solver field's local layout.
+
+        This is a low-level diagnostic hook for post-processing derived
+        boundary quantities such as consistent-boundary-flux traction. By
+        default it calls PETSc's ``DMPlexSNESComputeResidualFEM`` directly. If
+        ``cell_indices`` is supplied, it instead calls
+        ``DMPlexComputeResidualByKey`` on those local cells and the requested
+        test field. Boundary residuals registered through natural, Nitsche, or
+        multiplier boundary terms are not included. The returned arrays are
+        local to each rank and have the same flat layout as the corresponding
+        MeshVariable PETSc vector.
+        """
+        cdef DM _time_dm_residual
+        cdef DM dm
+        cdef Vec xvec
+        cdef Vec fvec
+        cdef PetscFormKey key
+        cdef IS ccell_is
+
+        self._build(verbose, False, None)
+
+        if time is not None:
+            if hasattr(time, 'magnitude') or hasattr(time, '_pint_qty'):
+                t_nd = float(uw.non_dimensionalise(time))
+            else:
+                t_nd = float(time)
+            _time_dm_residual = self.dm
+            UW_DMSetTime(_time_dm_residual.dm, t_nd)
+
+        self.mesh.update_lvec()
+        self.dm.setAuxiliaryVec(self.mesh.lvec, None)
+        self._update_constants()
+
+        gvec = self.dm.getGlobalVec()
+        xlocal = self.dm.getLocalVec()
+        flocal = self.dm.getLocalVec()
+        gvec.setArray(0.0)
+        xlocal.setArray(0.0)
+        flocal.setArray(0.0)
+
+        try:
+            for name, var in self.fields.items():
+                sgvec = gvec.getSubVector(self._subdict[name][0])
+                subdm = self._subdict[name][1]
+                subdm.localToGlobal(var.vec, sgvec)
+                gvec.restoreSubVector(self._subdict[name][0], sgvec)
+
+            self.dm.globalToLocal(gvec, xlocal)
+
+            dm = self.dm
+            xvec = xlocal
+            fvec = flocal
+            if cell_indices is None:
+                CHKERRQ(DMPlexSNESComputeResidualFEM(dm.dm, xvec.vec, fvec.vec, NULL))
+            else:
+                if residual_field_id is None:
+                    residual_field_id = 0
+                cell_is = PETSc.IS().createGeneral(
+                    list(cell_indices), comm=PETSc.COMM_SELF
+                )
+                try:
+                    ccell_is = cell_is
+                    key.label = NULL
+                    key.value = 0
+                    key.field = <PetscInt>residual_field_id
+                    key.part = 0
+                    CHKERRQ(DMPlexComputeResidualByKey(
+                        dm.dm, key, ccell_is.iset, <PetscReal>-1.7976931348623157e308,
+                        xvec.vec, NULL, 0.0, fvec.vec, NULL,
+                    ))
+                finally:
+                    cell_is.destroy()
+
+            local_section = self.dm.getLocalSection()
+            pStart, pEnd = local_section.getChart()
+            out = {}
+
+            for name, var in self.fields.items():
+                field_id = getattr(var, "_solver_field_id", None)
+                if field_id is None:
+                    field_id = getattr(var, "field_id", None)
+                if field_id is None:
+                    continue
+
+                is_field = None
+                created_is_field = False
+                if name == "velocity" and getattr(self, "_velocity_is", None) is not None:
+                    is_field = self._velocity_is
+                elif name == "pressure" and getattr(self, "_pressure_is", None) is not None:
+                    is_field = self._pressure_is
+                elif getattr(self, "_multiplier_is", None) is not None and name in self._multiplier_is:
+                    is_field = self._multiplier_is[name]
+                else:
+                    indices = []
+                    for point in range(pStart, pEnd):
+                        dof = local_section.getFieldDof(point, field_id)
+                        if dof > 0:
+                            offset = local_section.getFieldOffset(point, field_id)
+                            for i in range(dof):
+                                indices.append(offset + i)
+
+                    is_field = PETSc.IS().createGeneral(indices, comm=PETSc.COMM_SELF)
+                    created_is_field = True
+
+                try:
+                    subvec = flocal.getSubVector(is_field)
+                    try:
+                        out[name] = np.array(subvec.array, copy=True)
+                    finally:
+                        flocal.restoreSubVector(is_field, subvec)
+                finally:
+                    if created_is_field:
+                        is_field.destroy()
+
+            return out
+        finally:
+            self.dm.restoreLocalVec(flocal)
+            self.dm.restoreLocalVec(xlocal)
+            self.dm.restoreGlobalVec(gvec)
+
     @timing.routine_timer_decorator
     def solve(self,
               zero_init_guess: bool = True,

From eec6d0c7e1c2d1967f828fad492d5f95b5cb12c3 Mon Sep 17 00:00:00 2001
From: Tyagi <gthyagi@gmail.com>
Date: Thu, 18 Jun 2026 09:46:48 +0530
Subject: [PATCH 07/11] Add boundary weak residual diagnostic hook

Expose a small PETSc compatibility helper to fetch the boundary-specific PetscWeakForm registered by DMAddBoundary/PetscDSAddBoundary_UW.

Add SNES_Stokes_SaddlePt.compute_boundary_residual_fields() so benchmark diagnostics can assemble the registered Nitsche boundary residual for a named boundary and field through DMPlexComputeBdResidualSingle.

This fixes the earlier zero boundary-residual diagnostic, which used the global PetscDS weak form instead of the boundary weak form that actually stores UW3 Nitsche terms.
---
 src/underworld3/cython/petsc_compat.h         |   6 +
 src/underworld3/cython/petsc_extras.pxi       |   5 +
 .../cython/petsc_generic_snes_solvers.pyx     | 125 ++++++++++++++++++
 3 files changed, 136 insertions(+)

diff --git a/src/underworld3/cython/petsc_compat.h b/src/underworld3/cython/petsc_compat.h
index 2fe0c649e..1e74600f8 100644
--- a/src/underworld3/cython/petsc_compat.h
+++ b/src/underworld3/cython/petsc_compat.h
@@ -45,6 +45,12 @@ PetscErrorCode DMSetAuxiliaryVec_UW(DM dm, DMLabel label, PetscInt value, PetscI
     return DMSetAuxiliaryVec(dm, label, value, part, aux);
 }
 
+PetscErrorCode UW_PetscDSGetBoundaryWeakForm(PetscDS ds, PetscInt bd, PetscWeakForm *wf)
+{
+    PetscCall(PetscDSGetBoundary(ds, bd, wf, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
+    return PETSC_SUCCESS;
+}
+
 // copy paste function signitures from $PETSC_DIR/include/petscds.h - would be nice to automate this.
 #define UW_SIG_F0 PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]
 #define UW_SIG_G0 PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]
diff --git a/src/underworld3/cython/petsc_extras.pxi b/src/underworld3/cython/petsc_extras.pxi
index 941207e5e..029efeb81 100644
--- a/src/underworld3/cython/petsc_extras.pxi
+++ b/src/underworld3/cython/petsc_extras.pxi
@@ -32,6 +32,7 @@ cdef CHKERRQ(PetscErrorCode ierr):
     if ierr != 0: raise RuntimeError(f"PETSc error code '{interr}' was encountered.\nhttps://www.mcs.anl.gov/petsc/petsc-current/include/petscerror.h.html")
 
 cdef extern from "petscdstypes.h":
+    ctypedef void *PetscWeakForm "PetscWeakForm"
     ctypedef struct PetscFormKey:
         PetscDMLabel label
         PetscInt value
@@ -42,6 +43,7 @@ cdef extern from "petsc_compat.h":
 
     PetscErrorCode PetscDSAddBoundary_UW( PetscDM, DMBoundaryConditionType, const char[], const char[] , PetscInt, PetscInt, PetscInt *, void (*)(), void (*)(), PetscInt, const PetscInt *, void *)
     PetscErrorCode DMSetAuxiliaryVec_UW(PetscDM, PetscDMLabel, PetscInt, PetscInt, PetscVec)
+    PetscErrorCode UW_PetscDSGetBoundaryWeakForm(PetscDS, PetscInt, PetscWeakForm *)
     # PetscErrorCode UW_PetscDSSetBdResidual(PetscDS, PetscDMLabel, PetscInt, PetscInt, PetscInt, PetscInt, void*, PetscInt, void*)
 
     PetscErrorCode UW_PetscDSSetBdResidual(PetscDS, PetscDMLabel, PetscInt, PetscInt, PetscInt, PetscInt, void*, void*)
@@ -61,6 +63,7 @@ cdef extern from "petsc.h" nogil:
     PetscErrorCode DMPlexSNESComputeBoundaryFEM( PetscDM, void *, void *)
     PetscErrorCode DMPlexSNESComputeResidualFEM( PetscDM, PetscVec, PetscVec, void *)
     PetscErrorCode DMPlexComputeResidualByKey( PetscDM, PetscFormKey, PetscIS, PetscReal, PetscVec, PetscVec, PetscReal, PetscVec, void *)
+    PetscErrorCode DMPlexComputeBdResidualSingle( PetscDM, PetscWeakForm, PetscFormKey, PetscVec, PetscVec, PetscReal, PetscVec )
     # PetscErrorCode DMPlexSetSNESLocalFEM( PetscDM, void *, void *, void *)
     # PetscErrorCode DMPlexSetSNESLocalFEM( PetscDM, PetscBool, void *)
     PetscErrorCode DMPlexComputeGeometryFVM( PetscDM dm, PetscVec *cellgeom, PetscVec *facegeom)
@@ -71,6 +74,7 @@ cdef extern from "petsc.h" nogil:
     PetscErrorCode PetscDSSetJacobian( PetscDS, PetscInt, PetscInt, PetscDSJacobianFn, PetscDSJacobianFn, PetscDSJacobianFn, PetscDSJacobianFn)
     PetscErrorCode PetscDSSetJacobianPreconditioner( PetscDS, PetscInt, PetscInt, PetscDSJacobianFn, PetscDSJacobianFn, PetscDSJacobianFn, PetscDSJacobianFn)
     PetscErrorCode PetscDSSetResidual( PetscDS, PetscInt, PetscDSResidualFn, PetscDSResidualFn )
+    PetscErrorCode PetscDSGetWeakForm( PetscDS, PetscWeakForm * )
     
     PetscErrorCode PetscDSSetBdJacobian( PetscDS, PetscInt, PetscInt, PetscDSBdJacobianFn, PetscDSBdJacobianFn, PetscDSBdJacobianFn, PetscDSBdJacobianFn)
     PetscErrorCode PetscDSSetBdJacobianPreconditioner( PetscDS, PetscInt, PetscInt, PetscDSBdJacobianFn, PetscDSBdJacobianFn, PetscDSBdJacobianFn, PetscDSBdJacobianFn)
@@ -92,6 +96,7 @@ cdef extern from "petsc.h" nogil:
     PetscErrorCode DMGetRegionDS(PetscDM dm, PetscDMLabel label, PetscIS *fields, PetscDS *ds, PetscDS *dsIn)
     PetscErrorCode DMGetRegionNumDS(PetscDM dm, PetscInt num, PetscDMLabel *label, PetscIS *fields, PetscDS *ds, PetscDS *dsIn)
     PetscErrorCode DMSetRegionNumDS(PetscDM dm, PetscInt num, PetscDMLabel label, PetscIS fields, PetscDS ds, PetscDS dsIn)
+    PetscErrorCode DMGetDS(PetscDM dm, PetscDS *ds)
     PetscErrorCode DMGetNumDS(PetscDM dm, PetscInt *num)
     PetscErrorCode DMGetCellDS(PetscDM dm, PetscInt point, PetscDS *ds, PetscDS *dsIn)
     PetscErrorCode PetscDSSetCoordinateDimension(PetscDS ds, PetscInt dim)
diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
index 9eae9ee9d..f0b832e2f 100644
--- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx
+++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
@@ -6703,6 +6703,131 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
             self.dm.restoreLocalVec(xlocal)
             self.dm.restoreGlobalVec(gvec)
 
+    def compute_boundary_residual_fields(self, boundary, time=None, verbose=False, residual_field_id=0):
+        """Return the registered FEM boundary residual for one named boundary.
+
+        This is a low-level diagnostic hook for weak-boundary-condition
+        debugging. It assembles PETSc's boundary residual terms registered on
+        ``boundary`` through ``DMPlexComputeBdResidualSingle``. For Nitsche
+        free slip, this includes the full registered weak boundary residual,
+        not only the scalar penalty term. The returned arrays are local to each
+        rank and have the same flat layout as the corresponding MeshVariable
+        PETSc vector.
+        """
+        cdef DM _time_dm_boundary_residual
+        cdef DM dm
+        cdef Vec xvec
+        cdef Vec fvec
+        cdef PetscFormKey key
+        cdef PetscDS ds
+        cdef PetscWeakForm wf
+        cdef DMLabel c_label
+
+        self._build(verbose, False, None)
+
+        boundary_bc = None
+        for bc in self.natural_bcs:
+            if bc.boundary == boundary and bc.f_id == residual_field_id:
+                boundary_bc = bc
+                break
+        if boundary_bc is None:
+            raise ValueError(
+                f"No natural/Nitsche boundary residual is registered for "
+                f"boundary '{boundary}' and field {residual_field_id}."
+            )
+
+        if time is not None:
+            if hasattr(time, 'magnitude') or hasattr(time, '_pint_qty'):
+                t_nd = float(uw.non_dimensionalise(time))
+            else:
+                t_nd = float(time)
+            _time_dm_boundary_residual = self.dm
+            UW_DMSetTime(_time_dm_boundary_residual.dm, t_nd)
+
+        self.mesh.update_lvec()
+        self.dm.setAuxiliaryVec(self.mesh.lvec, None)
+        self._update_constants()
+
+        gvec = self.dm.getGlobalVec()
+        xlocal = self.dm.getLocalVec()
+        flocal = self.dm.getLocalVec()
+        gvec.setArray(0.0)
+        xlocal.setArray(0.0)
+        flocal.setArray(0.0)
+
+        try:
+            for name, var in self.fields.items():
+                sgvec = gvec.getSubVector(self._subdict[name][0])
+                subdm = self._subdict[name][1]
+                subdm.localToGlobal(var.vec, sgvec)
+                gvec.restoreSubVector(self._subdict[name][0], sgvec)
+
+            self.dm.globalToLocal(gvec, xlocal)
+
+            dm = self.dm
+            xvec = xlocal
+            fvec = flocal
+            CHKERRQ(DMGetDS(dm.dm, &ds))
+            CHKERRQ(UW_PetscDSGetBoundaryWeakForm(
+                ds, <PetscInt>boundary_bc.PETScID, &wf,
+            ))
+
+            c_label = self.dm.getLabel("UW_Boundaries")
+            key.label = c_label.dmlabel
+            key.value = <PetscInt>boundary_bc.boundary_label_val
+            key.field = <PetscInt>residual_field_id
+            key.part = 0
+            CHKERRQ(DMPlexComputeBdResidualSingle(
+                dm.dm, wf, key, xvec.vec, NULL, 0.0, fvec.vec,
+            ))
+
+            local_section = self.dm.getLocalSection()
+            pStart, pEnd = local_section.getChart()
+            out = {}
+
+            for name, var in self.fields.items():
+                field_id = getattr(var, "_solver_field_id", None)
+                if field_id is None:
+                    field_id = getattr(var, "field_id", None)
+                if field_id is None:
+                    continue
+
+                is_field = None
+                created_is_field = False
+                if name == "velocity" and getattr(self, "_velocity_is", None) is not None:
+                    is_field = self._velocity_is
+                elif name == "pressure" and getattr(self, "_pressure_is", None) is not None:
+                    is_field = self._pressure_is
+                elif getattr(self, "_multiplier_is", None) is not None and name in self._multiplier_is:
+                    is_field = self._multiplier_is[name]
+                else:
+                    indices = []
+                    for point in range(pStart, pEnd):
+                        dof = local_section.getFieldDof(point, field_id)
+                        if dof > 0:
+                            offset = local_section.getFieldOffset(point, field_id)
+                            for i in range(dof):
+                                indices.append(offset + i)
+
+                    is_field = PETSc.IS().createGeneral(indices, comm=PETSc.COMM_SELF)
+                    created_is_field = True
+
+                try:
+                    subvec = flocal.getSubVector(is_field)
+                    try:
+                        out[name] = np.array(subvec.array, copy=True)
+                    finally:
+                        flocal.restoreSubVector(is_field, subvec)
+                finally:
+                    if created_is_field:
+                        is_field.destroy()
+
+            return out
+        finally:
+            self.dm.restoreLocalVec(flocal)
+            self.dm.restoreLocalVec(xlocal)
+            self.dm.restoreGlobalVec(gvec)
+
     @timing.routine_timer_decorator
     def solve(self,
               zero_init_guess: bool = True,

From 734666a0d118e90df5b02a8b5960f1047914cff1 Mon Sep 17 00:00:00 2001
From: Tyagi <gthyagi@gmail.com>
Date: Thu, 18 Jun 2026 10:06:09 +0530
Subject: [PATCH 08/11] Clarify Stokes residual diagnostic semantics

Document that compute_volume_residual_fields uses PETSc's keyed residual path, which appends registered boundary residuals in the implicit Stokes mode. This prevents the helper from being treated as a true volume-only CBF recovery path.
---
 src/underworld3/cython/petsc_generic_snes_solvers.pyx | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
index f0b832e2f..89b63d6b3 100644
--- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx
+++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
@@ -6583,17 +6583,18 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
             print(f"Region DS: inactive region '{label_name}' gets trivial DS", flush=True)
 
     def compute_volume_residual_fields(self, time=None, verbose=False, cell_indices=None, residual_field_id=None):
-        """Return the volume-only FEM residual in each solver field's local layout.
+        """Return PETSc FEM residual fields in each solver field's local layout.
 
         This is a low-level diagnostic hook for post-processing derived
         boundary quantities such as consistent-boundary-flux traction. By
         default it calls PETSc's ``DMPlexSNESComputeResidualFEM`` directly. If
         ``cell_indices`` is supplied, it instead calls
         ``DMPlexComputeResidualByKey`` on those local cells and the requested
-        test field. Boundary residuals registered through natural, Nitsche, or
-        multiplier boundary terms are not included. The returned arrays are
-        local to each rank and have the same flat layout as the corresponding
-        MeshVariable PETSc vector.
+        test field. PETSc's keyed residual path also appends registered
+        boundary residuals, so this is a total residual diagnostic, not a
+        volume-only CBF recovery. The returned arrays are local to each rank
+        and have the same flat layout as the corresponding MeshVariable PETSc
+        vector.
         """
         cdef DM _time_dm_residual
         cdef DM dm

From 798dc33eb1edfda373788079a93e4e64c703e9de Mon Sep 17 00:00:00 2001
From: Tyagi <gthyagi@gmail.com>
Date: Thu, 18 Jun 2026 10:13:44 +0530
Subject: [PATCH 09/11] Add cloned-DM volume residual diagnostic

Add a UW PETSc compatibility helper that assembles DMPlexComputeResidualByKey on a cloned DM with copied sections, fields, equations, constants, and auxiliary vector, but without registered boundary objects. This gives selected-cell volume residuals without mutating the live solver PetscDS.

Update compute_volume_residual_fields() so selected-cell calls use the volume-only cloned-DM path by default, while include_boundary_terms=True preserves PETSc's original keyed residual behavior for cancellation/debug checks.
---
 src/underworld3/cython/petsc_compat.h         | 27 +++++++++++++
 src/underworld3/cython/petsc_extras.pxi       |  1 +
 .../cython/petsc_generic_snes_solvers.pyx     | 39 +++++++++++++------
 3 files changed, 56 insertions(+), 11 deletions(-)

diff --git a/src/underworld3/cython/petsc_compat.h b/src/underworld3/cython/petsc_compat.h
index 1e74600f8..70683fc46 100644
--- a/src/underworld3/cython/petsc_compat.h
+++ b/src/underworld3/cython/petsc_compat.h
@@ -51,6 +51,33 @@ PetscErrorCode UW_PetscDSGetBoundaryWeakForm(PetscDS ds, PetscInt bd, PetscWeakF
     return PETSC_SUCCESS;
 }
 
+PetscErrorCode UW_DMPlexComputeResidualByKeyVolumeOnly(DM dm, PetscFormKey key, IS cellIS, PetscReal time, Vec locX, Vec locX_t, PetscReal t, Vec locF, PetscCtx ctx)
+{
+    DM vdm = NULL;
+    PetscSection section = NULL;
+    PetscSection global_section = NULL;
+    PetscDS ds = NULL;
+    PetscDS vds = NULL;
+    Vec aux = NULL;
+
+    PetscCall(DMClone(dm, &vdm));
+    PetscCall(DMGetLocalSection(dm, &section));
+    PetscCall(DMSetLocalSection(vdm, section));
+    PetscCall(DMGetGlobalSection(dm, &global_section));
+    PetscCall(DMSetGlobalSection(vdm, global_section));
+    PetscCall(DMCopyFields(dm, PETSC_DETERMINE, PETSC_DETERMINE, vdm));
+    PetscCall(DMCreateDS(vdm));
+    PetscCall(DMGetDS(dm, &ds));
+    PetscCall(DMGetDS(vdm, &vds));
+    PetscCall(PetscDSCopyConstants(ds, vds));
+    PetscCall(PetscDSCopyEquations(ds, vds));
+    PetscCall(DMGetAuxiliaryVec(dm, key.label, key.value, key.part, &aux));
+    if (aux) PetscCall(DMSetAuxiliaryVec(vdm, key.label, key.value, key.part, aux));
+    PetscCall(DMPlexComputeResidualByKey(vdm, key, cellIS, time, locX, locX_t, t, locF, ctx));
+    PetscCall(DMDestroy(&vdm));
+    return PETSC_SUCCESS;
+}
+
 // copy paste function signitures from $PETSC_DIR/include/petscds.h - would be nice to automate this.
 #define UW_SIG_F0 PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]
 #define UW_SIG_G0 PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, PetscReal, const PetscReal[], const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]
diff --git a/src/underworld3/cython/petsc_extras.pxi b/src/underworld3/cython/petsc_extras.pxi
index 029efeb81..f864e0d8c 100644
--- a/src/underworld3/cython/petsc_extras.pxi
+++ b/src/underworld3/cython/petsc_extras.pxi
@@ -44,6 +44,7 @@ cdef extern from "petsc_compat.h":
     PetscErrorCode PetscDSAddBoundary_UW( PetscDM, DMBoundaryConditionType, const char[], const char[] , PetscInt, PetscInt, PetscInt *, void (*)(), void (*)(), PetscInt, const PetscInt *, void *)
     PetscErrorCode DMSetAuxiliaryVec_UW(PetscDM, PetscDMLabel, PetscInt, PetscInt, PetscVec)
     PetscErrorCode UW_PetscDSGetBoundaryWeakForm(PetscDS, PetscInt, PetscWeakForm *)
+    PetscErrorCode UW_DMPlexComputeResidualByKeyVolumeOnly( PetscDM, PetscFormKey, PetscIS, PetscReal, PetscVec, PetscVec, PetscReal, PetscVec, void *)
     # PetscErrorCode UW_PetscDSSetBdResidual(PetscDS, PetscDMLabel, PetscInt, PetscInt, PetscInt, PetscInt, void*, PetscInt, void*)
 
     PetscErrorCode UW_PetscDSSetBdResidual(PetscDS, PetscDMLabel, PetscInt, PetscInt, PetscInt, PetscInt, void*, void*)
diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
index 89b63d6b3..fc530d73f 100644
--- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx
+++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
@@ -6582,19 +6582,27 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         if uw.mpi.rank == 0 and self.verbose:
             print(f"Region DS: inactive region '{label_name}' gets trivial DS", flush=True)
 
-    def compute_volume_residual_fields(self, time=None, verbose=False, cell_indices=None, residual_field_id=None):
+    def compute_volume_residual_fields(
+        self,
+        time=None,
+        verbose=False,
+        cell_indices=None,
+        residual_field_id=None,
+        include_boundary_terms=False,
+    ):
         """Return PETSc FEM residual fields in each solver field's local layout.
 
         This is a low-level diagnostic hook for post-processing derived
         boundary quantities such as consistent-boundary-flux traction. By
         default it calls PETSc's ``DMPlexSNESComputeResidualFEM`` directly. If
         ``cell_indices`` is supplied, it instead calls
-        ``DMPlexComputeResidualByKey`` on those local cells and the requested
-        test field. PETSc's keyed residual path also appends registered
-        boundary residuals, so this is a total residual diagnostic, not a
-        volume-only CBF recovery. The returned arrays are local to each rank
-        and have the same flat layout as the corresponding MeshVariable PETSc
-        vector.
+        a UW wrapper around ``DMPlexComputeResidualByKey`` on a cloned DM with
+        a copied ``PetscDS`` that has no registered boundary objects, so the
+        selected-cell path returns volume terms only. Set
+        ``include_boundary_terms=True`` to call PETSc's original keyed
+        residual behavior, which appends registered boundary residuals. The
+        returned arrays are local to each rank and have the same flat layout as
+        the corresponding MeshVariable PETSc vector.
         """
         cdef DM _time_dm_residual
         cdef DM dm
@@ -6602,6 +6610,8 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         cdef Vec fvec
         cdef PetscFormKey key
         cdef IS ccell_is
+        cdef PetscReal residual_time = 0.0
+        cdef PetscReal implicit_form_time = <PetscReal>-1.7976931348623157e308
 
         self._build(verbose, False, None)
 
@@ -6612,6 +6622,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
                 t_nd = float(time)
             _time_dm_residual = self.dm
             UW_DMSetTime(_time_dm_residual.dm, t_nd)
+            residual_time = <PetscReal>t_nd
 
         self.mesh.update_lvec()
         self.dm.setAuxiliaryVec(self.mesh.lvec, None)
@@ -6650,10 +6661,16 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
                     key.value = 0
                     key.field = <PetscInt>residual_field_id
                     key.part = 0
-                    CHKERRQ(DMPlexComputeResidualByKey(
-                        dm.dm, key, ccell_is.iset, <PetscReal>-1.7976931348623157e308,
-                        xvec.vec, NULL, 0.0, fvec.vec, NULL,
-                    ))
+                    if include_boundary_terms:
+                        CHKERRQ(DMPlexComputeResidualByKey(
+                            dm.dm, key, ccell_is.iset, implicit_form_time,
+                            xvec.vec, NULL, residual_time, fvec.vec, NULL,
+                        ))
+                    else:
+                        CHKERRQ(UW_DMPlexComputeResidualByKeyVolumeOnly(
+                            dm.dm, key, ccell_is.iset, implicit_form_time,
+                            xvec.vec, NULL, residual_time, fvec.vec, NULL,
+                        ))
                 finally:
                     cell_is.destroy()
 

From ff700141a8cae36e4d4d5923b02d150d8e1c06d4 Mon Sep 17 00:00:00 2001
From: Tyagi <gthyagi@gmail.com>
Date: Thu, 18 Jun 2026 13:07:02 +0530
Subject: [PATCH 10/11] Fix residual diagnostic context signature

Use void * for the UW_DMPlexComputeResidualByKeyVolumeOnly context argument so the shared PETSc compatibility header builds in all Cython extension modules. The Cython declaration already exposes this wrapper with a void * context.

Validation: ./uw build passed. pytest -q tests/test_0502_boundary_integrals.py::test_bd_integral_spherical_internal_boundary_areas passed. pytest -q tests/test_1064_constrained_spherical_shell_response.py reported 2 passed and 2 expected xfails.
---
 src/underworld3/cython/petsc_compat.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/underworld3/cython/petsc_compat.h b/src/underworld3/cython/petsc_compat.h
index 70683fc46..74075e043 100644
--- a/src/underworld3/cython/petsc_compat.h
+++ b/src/underworld3/cython/petsc_compat.h
@@ -51,7 +51,7 @@ PetscErrorCode UW_PetscDSGetBoundaryWeakForm(PetscDS ds, PetscInt bd, PetscWeakF
     return PETSC_SUCCESS;
 }
 
-PetscErrorCode UW_DMPlexComputeResidualByKeyVolumeOnly(DM dm, PetscFormKey key, IS cellIS, PetscReal time, Vec locX, Vec locX_t, PetscReal t, Vec locF, PetscCtx ctx)
+PetscErrorCode UW_DMPlexComputeResidualByKeyVolumeOnly(DM dm, PetscFormKey key, IS cellIS, PetscReal time, Vec locX, Vec locX_t, PetscReal t, Vec locF, void *ctx)
 {
     DM vdm = NULL;
     PetscSection section = NULL;

From b4472b6aa2c846b34535ee3d9514ec843135eec3 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Thu, 18 Jun 2026 22:13:36 +1000
Subject: [PATCH 11/11] SphericalShellInternalBoundary: stop advertising
 non-existent Inner/Outer regions
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The single-volume cut+embed design (this PR) produces one shell volume with the
internal sphere embedded as a conformal surface — it has no Inner/Outer region
sub-volumes. But mesh.regions was still set to the Inner/Outer enum, so
extract_region('Inner') failed with an opaque 'Label not found on DM'. Leave
mesh.regions as None (drop the phantom enum) so region extraction reports a
clean 'no regions defined' and the API doesn't advertise a capability the DM
can't back. Use the 'Internal' boundary label for the internal interface.

Also re-tier the new 3D-mesh-gen area test to level_2/tier_b (it inherits the
module-level level_1/tier_a but builds a full 3D gmsh+embed mesh and the
generator is not yet production-soaked). Addresses the PR #242 review.

Underworld development team with AI support from Claude Code
---
 src/underworld3/meshing/spherical.py  | 14 ++++++++++----
 tests/test_0502_boundary_integrals.py |  9 ++++++++-
 2 files changed, 18 insertions(+), 5 deletions(-)

diff --git a/src/underworld3/meshing/spherical.py b/src/underworld3/meshing/spherical.py
index 2ea99fa9d..82cbf5371 100644
--- a/src/underworld3/meshing/spherical.py
+++ b/src/underworld3/meshing/spherical.py
@@ -557,9 +557,13 @@ class boundaries(Enum):
         Internal = 12
         Upper = 13
 
-    class regions(Enum):
-        Inner = 101
-        Outer = 102
+    # NOTE: this generator builds a SINGLE shell volume [radiusInner, radiusOuter]
+    # with the radiusInternal sphere *embedded* as a conformal internal surface
+    # (the `Internal` boundary). Unlike the old occ.fragment approach it does NOT
+    # split the volume into Inner/Outer region sub-volumes, so no Inner/Outer
+    # region physical groups exist and `mesh.regions` is left as None — region
+    # extraction is intentionally unsupported here (use the `Internal` boundary
+    # label for the internal interface). See PR #242.
 
     import gmsh
 
@@ -741,7 +745,9 @@ def spherical_mesh_refinement_callback(dm):
 
     # boundary_normals deprecated — use mesh.Gamma_P1 for boundary normals
 
-    new_mesh.regions = regions
+    # Single-volume embed design (see note above): no Inner/Outer region groups
+    # exist on the DM, so leave mesh.regions as None rather than advertising
+    # labels that extract_region() cannot resolve.
 
     # Full spherical shell with internal boundary: 3 rigid rotation modes
     x, y, z = new_mesh.X
diff --git a/tests/test_0502_boundary_integrals.py b/tests/test_0502_boundary_integrals.py
index 1721db7db..9080fd9ad 100644
--- a/tests/test_0502_boundary_integrals.py
+++ b/tests/test_0502_boundary_integrals.py
@@ -320,8 +320,15 @@ def _get_spherical_internal_mesh():
     return _mesh_spherical_internal
 
 
+@pytest.mark.level_2
+@pytest.mark.tier_b
 def test_bd_integral_spherical_internal_boundary_areas():
-    """SphericalShellInternalBoundary preserves Lower/Internal/Upper labels."""
+    """SphericalShellInternalBoundary preserves Lower/Internal/Upper labels.
+
+    Overrides the module-level level_1/tier_a marks: this builds a full 3D
+    gmsh+embed mesh (not a seconds-scale level_1 op), and the embed generator
+    is not yet production-soaked for tier_a. See PR #242 review.
+    """
 
     mesh_spherical = _get_spherical_internal_mesh()
     expected_areas = {