Free surface in 3D: spherical shell support, dimension-general surface machinery, and the 2D compatibility-floor fix

docs/developer/CHANGELOG.md

@@ -6,6 +6,39 @@ This log tracks significant development work at a conceptual level, suitable for
 
 ## 2026 Q3 (July – September)
 
+### The Free Surface Reaches the Spherical Shell (July 2026)
+
+**`uw.systems.FreeSurface` now runs in 3D on a spherical shell** — the same
+exponential three-number integrator, held-lid σ_nn recovery and strong
+material-boundary datum, with the surface machinery made dimension-general
+rather than ported piecewise:
+
+- The datum gauge (mean removal) is an FE trace-mass reduction over
+  owned boundary facets — no ordered ring, no gather; the same code is the
+  2D line gauge and the 3D area gauge. On the way it resolved a real 2D
+  defect: the deforming-ring strong-datum solves used to stall at a ~2e-3
+  residual floor, which turned out to be three stacked causes (arc-length vs
+  FE trace weights; the datum's *directed* mean flux through the deformed
+  facet normals, now stripped with the same FE surface integral the residual
+  uses; and the constant-pressure gauge mode, which the inner solver projects
+  and the outer loop therefore now measures in the quotient space). With all
+  three closed, every step of the power-law acceptance run converges.
+- σ_nn on a 3D P2 boundary is recovered by **P1 projection** (edge-midpoint
+  loads folded exactly onto vertices, sound P1 lumped triangle mass) — chosen
+  over the consistent P2 mass because its vertex-integral checkerboard sits
+  exactly at the vertices the P1 topography field consumes.
+- The two genuinely 2D features (ring Taubin filter, tangential transport)
+  are refused explicitly in 3D; everything else is shared code.
+
+First 3D evidence (spherical Y20 topographic relaxation, constant-density
+shell): exponential decay at an O(1) shell correction below the half-space
+Cathles rate, in the physically correct direction, with the equilibrium
+modal bias falling 16% → 2% of the initial amplitude over one resolution
+step (the known discrete recovery defect, resolution-convergent). The
+detailed benchmarking — analytic shell-rate comparison, convergence study,
+low-Ra spherical convection, 3D parallel — is deliberately left to the
+review pass.
+
 ### One Owner for the Geometric-Multigrid Option Bundle (July 2026)
 
 **The PETSc option bundle that configures a Stokes velocity block's multigrid

src/underworld3/cython/petsc_generic_snes_solvers.pyx

@@ -6053,10 +6053,13 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         :meth:`solve`.
 
         ``mass="auto"`` (default) uses lumped recovery for 2D traces and 3D P1
-        triangles, and the required consistent surface-mass solve for 3D P2 triangles.
+        triangles, and the consistent surface-mass solve for 3D P2 triangles.
         Explicit ``"lumped"`` and ``"consistent"`` choices remain available where
-        mathematically valid. Three-dimensional recovery currently supports triangular
-        P1/P2 traces only.
+        mathematically valid, and ``"p1"`` selects P1-PROJECTED recovery on a 3D
+        P2 trace (edge-midpoint loads folded onto vertices, lumped P1 triangle
+        mass — sound where the consistent P2 path carries the vertex-integral
+        checkerboard; the FreeSurface default in 3D). Three-dimensional recovery
+        currently supports triangular P1/P2 traces only.
 
         .. warning::
            On CURVED boundaries, P2 vertex values of :math:`\sigma_{nn}` converge

src/underworld3/utilities/boundary_flux.py

@@ -165,8 +165,9 @@ def _desmear(solver, boundary, xs, R, mass, remove_mean, partial_reaction=True):
     csec = dm.getCoordinateSection()
     cvec = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, dim)
     v0, v1 = dm.getDepthStratum(0)
-    if mass not in ("auto", "lumped", "consistent"):
-        raise ValueError("mass must be 'auto', 'lumped', or 'consistent'.")
+    if mass not in ("auto", "lumped", "consistent", "p1"):
+        raise ValueError("mass must be 'auto', 'lumped', 'consistent', or 'p1' "
+                         "(P1-projected recovery on a 3D P2 trace).")
     if dim == 3:
         lsec = dm.getLocalSection()
         ncomp = lsec.getFieldComponents(0)
@@ -260,8 +261,44 @@ def coord(q):
         if order == 2 and mass == "lumped":
             raise ValueError(
                 "A 3D P2 triangular trace has zero row-sum mass at its vertices; "
-                "use mass='consistent' for pointwise boundary-flux recovery."
+                "use mass='consistent' (pointwise, carries the vertex-integral "
+                "checkerboard risk) or mass='p1' (P1-projected, monotone — the "
+                "choice for driving a P1 surface field) for boundary-flux recovery."
             )
+        mid_owners = {}
+        if mass == "p1":
+            if order != 2:
+                mass = "lumped"                    # P1 trace: p1 IS lumped
+            else:
+                # P1-PROJECTED recovery on a P2 trace: the consistent P2 path has
+                # the ∫φ_vertex = 0 vertex checkerboard (#404 hold), while the P1
+                # trace is sound — and a P1 surface field only consumes vertex
+                # values anyway. Fold each edge-midpoint load onto its two edge
+                # vertices (φ^{P1}(edge-mid) = 1/2 exactly, P1 ⊂ P2 — the load
+                # transfer is the interpolation transpose, so the total load is
+                # conserved), then de-smear with the P1 lumped triangle mass.
+                # Midpoint outputs are read back as the P1 interpolant (vertex
+                # average).
+                new_elements = {}
+                for _order, nodes, area in elements.values():
+                    vk = nodes[:3]
+                    m01, m12, m20 = nodes[3:]
+                    mid_owners[m01] = (vk[0], vk[1])
+                    mid_owners[m12] = (vk[1], vk[2])
+                    mid_owners[m20] = (vk[2], vk[0])
+                    new_elements[(1, tuple(sorted(vk)))] = (1, vk, area)
+                folded = {}
+                for key, value in R_by.items():
+                    if key in mid_owners:
+                        va, vb = mid_owners[key]
+                        folded[va] = folded.get(va, 0.0) + 0.5 * value
+                        folded[vb] = folded.get(vb, 0.0) + 0.5 * value
+                    else:
+                        folded[key] = folded.get(key, 0.0) + value
+                R_by = folded
+                elements = new_elements
+                order = 1
+                mass = "lumped"
 
         keys = sorted(R_by)
         global_index = {key: i for i, key in enumerate(keys)}
@@ -311,14 +348,23 @@ def coord(q):
         if remove_mean:
             mean = float(np.dot(flux, boundary_mass) / np.sum(boundary_mass))
             flux -= mean
-        return np.array([flux[global_index[_key(x, dim)]] for x in xs])
+
+        def value_at(x):
+            key = _key(x, dim)
+            if key in global_index:
+                return flux[global_index[key]]
+            # P1-projected mode: a P2 edge midpoint reads the P1 interpolant
+            va, vb = mid_owners[key]
+            return 0.5 * (flux[global_index[va]] + flux[global_index[vb]])
+
+        return np.array([value_at(x) for x in xs])
 
     if dim != 2:
         raise NotImplementedError(
             f"Boundary-flux recovery is not implemented for mesh dimension {dim}."
         )
-    if mass == "auto":
-        mass = "lumped"
+    if mass in ("auto", "p1"):
+        mass = "lumped"                    # 2D: the lumped line mass is sound
 
     e0, e1 = dm.getDepthStratum(1)
     def vcoord(q): return cvec[csec.getOffset(q) // dim]

src/underworld3/utilities/rotated_bc.py

@@ -657,13 +657,37 @@ def solve_rotated_freeslip(solver, boundaries, remove_rotation_gauge=True,
     vel_its_last = []
     pres_its_last = []
 
+    # With the constant-pressure nullspace active, the outer residual is measured
+    # in the pressure-gauge QUOTIENT space: the inner KSP projects the constant
+    # mode out of every increment (it is the gauge of an enclosed incompressible
+    # domain), so the loop cannot reduce that component and must not measure it.
+    # On a DEFORMED faceted boundary the component is not exactly zero — free
+    # tangential DOFs carry a small net flux through the node-vs-facet normal
+    # mismatch, an irreducible discrete incompatibility (measured: an unprojected
+    # outer norm floors at rel ~2e-3, 100% pressure rows, and the line search
+    # stalls against the constant offset). This mirrors PETSc's own projected
+    # residual for singular systems. The Cartesian reaction stash is UNPROJECTED
+    # (σ_nn reads velocity rows only).
+    use_pnull = bool(getattr(solver, "_petsc_use_pressure_nullspace", False))
+    # L2 norm of the LAST projected-out gauge component (|mean|·√N over the
+    # pressure rows) — kept observable so an INCOMPATIBLE datum cannot hide
+    # behind the projection (see the guard after the Newton loop).
+    pnull_gauge = [0.0]
+
     def rotated_residual(uvec, keep_cartesian=False):
         snes.computeFunction(uvec, Fc)
         if keep_cartesian:
             Fc.copy(reaction)                # stash the Cartesian reaction for σ_nn
         Fh = Fc.duplicate()
         Q.mult(Fc, Fh)
         _zero_rows_local(Fh, normal_rows)
+        if use_pnull:
+            sp = Fh.getSubVector(pres_is)
+            n_p = sp.getSize()
+            mean = sp.sum() / max(n_p, 1)
+            sp.shift(-mean)                  # project out the pressure-gauge mode
+            Fh.restoreSubVector(pres_is, sp)
+            pnull_gauge[0] = abs(mean) * (max(n_p, 1) ** 0.5)
         return Fh
 
     # Convergence reference: max(initial residual, REST-STATE residual ‖F̂(0)‖).
@@ -831,6 +855,25 @@ def rotated_residual(uvec, keep_cartesian=False):
                    f"in {newton_its} iterations (rel |F̂| = {rel:.2e} of the reference "
                    f"residual); the fields hold the last (unconverged) iterate.")
 
+    # The quotient projection makes the outer test blind to the pressure-gauge
+    # component. A COMPATIBLE datum leaves that component at discretisation
+    # level; an INCOMPATIBLE one (net wall-normal flux into an enclosed
+    # incompressible domain) parks an O(forcing) constant there — and the loop
+    # above would report convergence while mass conservation is violated.
+    # Surface the component: readable state always, loud when it dominates.
+    if use_pnull:
+        solver._rotated_pressure_gauge_residual = pnull_gauge[0]
+        if converged and pnull_gauge[0] > 10.0 * max(rnorm, atol):
+            mpi.pprint(f"[rotated_bc] WARNING: converged in the pressure-gauge "
+                       f"QUOTIENT space, but the projected-out gauge component "
+                       f"(|F̂_p·1|/√N = {pnull_gauge[0]:.2e}) dominates the "
+                       f"converged residual ({rnorm:.2e}). The wall-normal datum "
+                       f"carries a net boundary flux this enclosed incompressible "
+                       f"domain cannot absorb — the velocity field violates mass "
+                       f"conservation at that level. Check the datum's surface "
+                       f"integral (solver._rotated_pressure_gauge_residual holds "
+                       f"this number).")
+
     Fc.destroy()                     # residual output buffer (reaction persists in the result dict)
     _destroy_rotated_ksp_ctx(ctx)            # KSP/PC + the owned Schur pmat
     if Ahat is not None:
@@ -1294,7 +1337,12 @@ def boundary_normal_traction(solver, boundary, solve_result, mass="auto"):
         3D P2 triangles.
       * ``"lumped"`` — the diagonal row-sum mass. It is monotone for supported traces,
         but invalid for 3D P2 triangles because their vertex row sums are exactly zero.
-      * ``"consistent"`` — the full trace mass. Required for pointwise 3D P2 recovery.
+      * ``"consistent"`` — the full trace mass. Pointwise-exact 3D P2 recovery, but
+        carries the vertex-integral checkerboard on P2 triangles (#404 hold).
+      * ``"p1"`` — P1-PROJECTED recovery on a 3D P2 trace (edge-midpoint loads folded
+        onto vertices, lumped P1 triangle mass). Sound where the consistent P2 path
+        checkerboards; the FreeSurface default in 3D. On a P1 trace, identical to
+        ``"lumped"``.
 
     Parallel-safe: r_c is scattered to a local vector (ghosts included) and read by LOCAL
     section offset; the boundary mass is assembled globally by a coordinate-keyed

tests/test_1070_free_surface_plume.py

@@ -311,8 +311,19 @@ def test_freesurface_ring_quadrature_is_exact_in_parallel():
     fs = uw.systems.FreeSurface(stokes, "Upper", buoyancy_scale=1.0, normal=rhat)
 
     weights = fs._ring_weights()
-    assert abs(float(weights.sum()) - 2.0 * np.pi * r_out) < 1.0e-9, \
-        f"ring weights sum to {weights.sum():.10f}, not the circumference (seam double-count?)"
+    # The gauge must match the FE boundary quadrature EXACTLY: the datum's
+    # flux-free condition lives in the discrete space, so the weight total is the
+    # FE measure of the (polygonal) boundary — NOT the ideal-circle arc length,
+    # which differs at O(h^2) and was the strong-datum compatibility floor
+    # (block-split evidence: the whole stalled residual sat in the pressure rows).
+    # A seam double-count breaks this equality at the first shared node.
+    fe_len = float(uw.maths.BdIntegral(mesh=mesh, fn=sympy.S.One,
+                                       boundary="Upper").evaluate())
+    assert abs(float(weights.sum()) - fe_len) < 1.0e-9, \
+        f"ring weights sum to {weights.sum():.10f}, FE boundary measure is {fe_len:.10f} " \
+        "(seam double-count?)"
+    # sanity: the polygonal measure approximates the circle at O(h^2)
+    assert abs(fe_len - 2.0 * np.pi * r_out) < 1.0e-2
 
     coords = fs._ring_coords
     theta = np.arctan2(coords[:, 1], coords[:, 0])

tests/test_1072_free_surface_spherical.py

@@ -0,0 +1,59 @@
+"""3D FreeSurface on a spherical shell: the end-to-end loop must run and produce
+physically sensible topography (guards the dimension-general surface machinery:
+owned-facet trace-mass gauge, P1-projected sigma_nn recovery, radial deform).
+
+The quantitative benchmarking (analytic Y_lm shell rate, convergence of the
+h_inf modal bias, 3D parallel) is the review-team's scope — this test pins the
+CAPABILITY: construction, one solve/advance cycle, finite mean-free h_inf, and
+the explicit refusal of the 2D-only features.
+"""
+import numpy as np
+import pytest
+import sympy
+import underworld3 as uw
+
+pytestmark = [pytest.mark.level_2, pytest.mark.tier_b]
+
+
+def _shell_stokes(cell=0.35):
+    mesh = uw.meshing.SphericalShell(radiusOuter=1.0, radiusInner=0.547,
+                                     cellSize=cell, qdegree=3)
+    x, y, z = mesh.X
+    r = sympy.sqrt(x ** 2 + y ** 2 + z ** 2)
+    rhat = sympy.Matrix([[x / r, y / r, z / r]])
+    stokes = uw.systems.Stokes(mesh)
+    stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
+    stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0
+    blob = sympy.exp(-(((x - 0.75) ** 2 + y ** 2 + z ** 2) / 0.05))
+    stokes.bodyforce = 50.0 * blob * rhat.T
+    stokes.add_essential_bc((0.0, 0.0, 0.0), "Lower")
+    stokes.tolerance = 1.0e-5
+    return mesh, stokes, rhat
+
+
+def test_freesurface_spherical_shell_end_to_end():
+    """Construction + one full solve/advance on the shell; h_inf finite and
+    mean-free; the surface responds toward equilibrium (|h| grows from flat
+    under the one-sided load and stays bounded by |h_inf|)."""
+    mesh, stokes, rhat = _shell_stokes()
+    fs = uw.systems.FreeSurface(stokes, "Upper", buoyancy_scale=50.0, normal=rhat)
+    fs.solve()
+    h_inf = np.asarray(fs._h_inf)
+    assert np.isfinite(h_inf).all(), "3D h_inf recovery produced non-finite values"
+    assert abs(fs._surface_mean(h_inf)) < 1.0e-8 * (np.abs(h_inf).max() + 1e-30), \
+        "h_inf datum is not mean-free under the trace-mass gauge"
+    assert np.abs(h_inf).max() > 1.0e-4, "no topographic response to the load"
+    fs.advance(fs.estimate_dt(advect_scale=10.0))
+    shape = fs._current_shape()
+    assert np.isfinite(shape).all()
+    assert 0.0 < np.abs(shape).max() <= 1.5 * np.abs(h_inf).max(), \
+        "surface did not move toward (or overshot) equilibrium"
+
+
+def test_freesurface_spherical_refuses_2d_only_features():
+    """The 2D-only features fail loudly at construction in 3D, not silently."""
+    mesh, stokes, rhat = _shell_stokes()
+    with pytest.raises(NotImplementedError, match="tangential"):
+        uw.systems.FreeSurface(stokes, "Upper", normal=rhat, tangent_advect="shape")
+    with pytest.raises(NotImplementedError, match="filter"):
+        uw.systems.FreeSurface(stokes, "Upper", normal=rhat, surface_filter=10)