Rebuild P2 boundary-recovery vertices from the midpoints, and exonerate the grad-div penalty (#633)

docs/developer/CHANGELOG.md

@@ -6,6 +6,76 @@ This log tracks significant development work at a conceptual level, suitable for
 
 ## 2026 Q3 (July – September)
 
+### A Singular Recovery Mass, Mistaken for a Penalty Defect (August 2026)
+
+**The grad-div penalty default stays off**, but the reason it was held off turned out to
+be a defect somewhere else entirely (#633) — so the objection that had blocked it is gone,
+and a different one took its place.
+
+With the penalty at 10, the spherical dynamic topography recovered from the rotated
+free-slip reaction dropped 28% at *vertices* while the facet-integrated value stayed
+correct. The natural reading — that grad-div augmentation corrupts the de-smearing from
+reaction loads to pointwise stress — was wrong.
+
+The de-smearing mass for a 3-D **P2 triangular** trace has vertex rows that sum to
+**exactly zero**. Those rows annihilate a constant, so solving `M σ = R` amplifies any
+perturbation of the nodal load at vertices by O(1) — and, being an instability rather
+than a discretisation error, does so independently of mesh resolution. The recovery was
+already 7.6% low with no penalty at all; the penalty only made it large enough to fail a
+test whose 12% tolerance had been hiding it.
+
+The discrimination needed a case that was curved but not 3-D. A 2-D annulus reproduces
+the signature exactly and then parts company under refinement: its error falls ~O(h²)
+while the shell's stays flat at ~0.28 over a 3.2× node-count range. The 2-D P2 **line**
+mass has positive vertex row sums, which is why 2-D never showed the defect and why
+dimension, curvature and the rotated constraint were all red herrings.
+
+- The zero-mean P2 vertex basis was **already known and documented in #414**, which
+  recorded the same drift-away-under-refinement we re-measured here. Its mechanism is
+  the sharper one and is adopted: because the vertex basis has zero surface mean, the
+  vertex reaction carries essentially only the O(h) facet-normal/geometry error, and the
+  consistent solve *faithfully reconstructs that error* — it is not amplifying noise.
+  What #633 adds is the separation from the grad-div penalty (which was blamed for it)
+  and the fix below, which is #414's own unactioned recommendation (2).
+- So `mass="auto"` stops asking. On a 3-D P2 trace it now takes the **consistent** solve,
+  keeps its superconvergent midpoints, and **reconstructs the vertices from them**: the
+  three midpoints of a facet determine a unique linear function, so a vertex reads its
+  two adjacent midpoints and subtracts the opposite one, averaged over incident facets.
+  Worst-node error against the analytic coefficient, over cellSize 0.25 → 0.11:
+
+  | | 0.25 | 0.20 | 0.16 | 0.13 | 0.11 |
+  |---|---|---|---|---|---|
+  | surface, P1-projected | 0.041 | 0.026 | 0.018 | 0.013 | 0.008 |
+  | surface, reconstructed | 0.016 | 0.012 | 0.012 | 0.003 | 0.004 |
+  | CMB, P1-projected | 0.116 | 0.067 | 0.047 | 0.030 | 0.025 |
+  | CMB, reconstructed | 0.094 | 0.058 | 0.043 | 0.024 | 0.015 |
+
+  Better at every resolution on both boundaries, by 1.8x to 4.9x, and converging. The
+  simpler P1-projected recovery stays available as `mass="p1"` — it is sound, it just
+  discards the good data along with the bad.
+- `FreeSurface` already used the P1-projected recovery in 3-D, so production dynamic
+  topography was never affected. The exposure was `mass="auto"`.
+- The spherical topography test is refined (cellSize 0.25 → 0.13) and its tolerances
+  tightened from 0.10/0.12 to 0.01/0.05, set from measured discretisation error with
+  ~2× headroom, and now assert every node class rather than the aggregate — the failure
+  was confined to one class and an aggregate assertion passed straight through it.
+- Filed #637: 3-D recovery accepts only P1/P2 triangular traces, so dynamic topography
+  has exactly one supported discretisation there and cannot be cross-validated. That
+  blocked the P3/hex arm of this investigation.
+- `Stokes.DEFAULT_PENALTY` was flipped to 10 on the #625 evidence and then **reverted**.
+  Three tier-A/B tests fail at 10 and pass at 0, and the same three run **8.3x slower**
+  (9.27 s to 76.92 s, warm cache both ways). The #625 win needs an FMG hierarchy; without
+  `refinement>=1` the velocity block falls back to GAMG, which is where grad-div
+  augmentation drives the solve into its iteration cap. The default path is the one
+  without a hierarchy, so the default serves it; set `penalty=10` explicitly where FMG
+  is available.
+- Two of those three failures are not penalty defects. The Nitsche free-slip leak
+  (1.234e-4 against a 1e-4 bound) is augmentation perturbing a *weakly* imposed
+  constraint — a strong rotated constraint is untouched. The swarm one exposed #641:
+  `evaluate` returns −0.4976 for `sqrt((E**2).trace()/2)` at in-domain points near the
+  lid-corner singularity, at `penalty=0` as well; the penalty merely moved an accumulated
+  total across zero.
+
 ### The Free Surface Reaches the Spherical Shell (July 2026)
 
 **`uw.systems.FreeSurface` now runs in 3D on a spherical shell** — the same

docs/developer/subsystems/rotated-freeslip.md

@@ -324,6 +324,27 @@ the shared boundary-mass machinery in `utilities/boundary_flux.py`;
 `dynamic_topography_field` writes `h = −(σ_nn − σ̄_nn)/(Δρ g)` onto a surface
 field for the free-surface integrator.
 
+### Which de-smearing mass (3-D)
+
+`mass="auto"` uses the **P1-projected** recovery on a 3-D P2 trace. The full
+("consistent") P2 **triangle** mass has vertex rows summing to exactly zero, so
+it is singular on constants along those rows and `M⁻¹` amplifies any
+perturbation of the nodal load at *vertices* by O(1) — independently of
+resolution. Measured on the Zhong l=2 shell (#633): consistent-mass vertex
+values are 7.6% low with no grad-div penalty, 28% low at `penalty=10` and 79%
+low at 100, and the error is flat across a 3.2× refinement. The 2-D P2 **line**
+mass has positive vertex row sums and is unaffected — which is why this appears
+only in 3-D, and not because of dimension or curvature as such.
+
+This is a trade rather than a strict improvement. Consistent-mass *midpoint*
+values are superconvergent (0.1–1.5% at every penalty) and the P1 projection
+gives that up, because under it the midpoints are the P1 interpolant of the
+vertices. `"p1"` is the default because its error converges under refinement
+(worst node 0.170 → 0.049 over cellSize 0.30 → 0.16) where the consistent
+vertex error does not. **If you sample only midpoints, pass
+`mass="consistent"`.** See also #404 (the vertex-integral checkerboard) and
+#637 (only P1/P2 triangular traces are supported in 3-D at all).
+
 ## Tests
 
 `tests/test_1018_rotated_freeslip.py` (serial: essential-equivalence, FMG,

src/underworld3/cython/petsc_generic_snes_solvers.pyx

@@ -3168,11 +3168,12 @@ class SolverBaseClass(uw_object):
         number); for a **vector** solver the traction :math:`\sigma\cdot\hat n` (pass
         ``normal`` to get the scalar normal component :math:`\hat n\cdot\sigma\cdot\hat n`).
 
-        ``mass`` de-smears the nodal reaction with ``"lumped"`` or ``"consistent"``
-        boundary mass. ``"auto"`` (default) selects lumped recovery for 2D P1/P2
-        traces and 3D P1 triangles, and the consistent solve for 3D P2 triangles and
-        2D traces of degree >= 3 (where row-sum lumping is respectively invalid and
-        only O(h) pointwise).
+        ``mass`` de-smears the nodal reaction with ``"lumped"``, ``"consistent"``,
+        ``"p1"`` or ``"midpoint"`` boundary mass. ``"auto"`` (default) selects lumped
+        recovery for 2D P1/P2 traces and 3D P1 triangles, MIDPOINT-RECONSTRUCTED
+        recovery for 3D P2 triangles (row-sum lumping is invalid there and the
+        consistent solve amplifies at vertices — #633), and the consistent solve for 2D
+        traces of degree >= 3 (where lumping is only O(h) pointwise).
         ``remove_mean`` subtracts the boundary mean — leave ``False`` for a physical
         flux (the mean is the Nusselt number); ``True`` gives a gauge-free field.
 
@@ -6455,13 +6456,14 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         :meth:`solve`.
 
         ``mass="auto"`` (default) uses lumped recovery for 2D traces and 3D P1
-        triangles, and the consistent surface-mass solve for 3D P2 triangles.
-        Explicit ``"lumped"`` and ``"consistent"`` choices remain available where
-        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.
+        triangles, and MIDPOINT-RECONSTRUCTED recovery for 3D P2 triangles (the
+        consistent solve, keeping its superconvergent midpoints, with vertices rebuilt
+        from them). ``"p1"`` selects the simpler P1-projected recovery; explicit
+        ``"lumped"`` and ``"consistent"`` remain available where mathematically valid;
+        ``"consistent"`` is pointwise-exact on a P2 trace in exact arithmetic but
+        its zero vertex row sums amplify any load perturbation at VERTICES by O(1),
+        independently of h (#404, measured in #633). Three-dimensional recovery
+        currently supports triangular P1/P2 traces only (#637).
 
         .. warning::
            On CURVED boundaries, P2 vertex values of :math:`\sigma_{nn}` converge
@@ -6486,8 +6488,8 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         interior left untouched.
 
         ``buoyancy_scale`` is :math:`\Delta\rho\,g` (traction → length).
-        ``mass="auto"`` selects lumped recovery where valid and the consistent
-        surface-mass solve for 3D P2 triangles. Requires a prior
+        ``mass="auto"`` selects lumped recovery where valid and
+        midpoint-reconstructed recovery for 3D P2 triangles. Requires a prior
         :meth:`add_rotated_freeslip_bc` on ``boundary`` and a completed :meth:`solve`.
 
         .. warning::

src/underworld3/systems/solvers.py

@@ -2183,16 +2183,25 @@ def penalty(self, value):
     #: not inside P0, so the term does not vanish at the discrete solution --
     #: but it converges away. ``penalty = 0`` restores the unaugmented operator.
     #:
-    #: **Held at 0 pending the pointwise-traction question.** At 10 the
-    #: spherical dynamic topography recovered from the rotated free-slip
-    #: reaction drops 0.4192 -> 0.3021, 28%, on the *vertex-sampled* value while
-    #: the facet-integrated value stays correct (``test_1018``). So augmentation
-    #: corrupts the de-smearing from reaction loads to pointwise stress —
-    #: presumably because lambda*mu*(div u) is non-zero cell-by-cell for P2-P0
-    #: and averages out over a facet integral but not at a vertex. Dynamic
-    #: topography is the main product of that machinery, so the value stays 0
-    #: until that is resolved; everything needed for the change is in place and
-    #: it is this constant.
+    #: **Held at 0, but no longer because of #633.** That issue turned out to be
+    #: the recovery operator, not the penalty: the 3-D P2 TRIANGLE de-smearing
+    #: mass has vertex rows summing to exactly zero, and ``mass="auto"`` now
+    #: picks the monotone P1-projected recovery instead. The topography
+    #: objection is gone.
+    #:
+    #: It stays at 0 because of the GAMG cost above, which is far worse than
+    #: #625 recorded. Flipping it to 10 was tried and reverted: three tier-A/B
+    #: tests failed (``test_0112`` swarm, ``test_1060`` Nitsche free-slip leak
+    #: 1.234e-4 against a 1e-4 bound, ``test_1061`` topography correlation
+    #: 0.979 against 0.99) and the same three ran **8.3x slower** -- 9.27 s to
+    #: 76.92 s, warm JIT cache both ways. The 21% FMG win needs a hierarchy;
+    #: without ``refinement>=1`` the velocity block falls back to GAMG, which
+    #: is the default path and the one that pays. A default has to serve the
+    #: default.
+    #:
+    #: Set ``penalty = 10`` explicitly on a solver that has an FMG hierarchy.
+    #: Note it also perturbs a Nitsche free-slip constraint, whose leak bound
+    #: is not slack enough to absorb it.
     DEFAULT_PENALTY = 0.0
 
 

src/underworld3/utilities/boundary_flux.py

@@ -60,6 +60,44 @@
 )
 
 
+def _rebuild_vertices_from_midpoints(elements, global_index, flux):
+    """Replace the VERTEX entries of a 3D P2 recovery with values reconstructed from
+    its edge-midpoint entries, averaged over the facets incident on each vertex.
+
+    The P2 triangle vertex basis function has zero mean (``_P2_TRIANGLE_MASS`` vertex
+    row sums are exactly 0), so the vertex DOF of an L2 recovery carries no mass and is
+    fixed by cancellation — O(1) error, independent of h. The midpoints carry all of it
+    and are superconvergent, so the vertices are better obtained from them than asked
+    for directly.
+
+    On one triangle the three midpoint values determine a unique linear function.
+    Because ``m01 = (v0+v1)/2`` and so on, that function takes the value
+    ``m01 + m20 - m12`` at ``v0`` — each vertex reads its two ADJACENT midpoints and
+    subtracts the OPPOSITE one. The subtraction amplifies noise up to 3x, which is
+    tolerable only because the midpoints are superconvergent; averaging over the
+    (typically ~6) incident facets damps it further.
+
+    Mean removal is unaffected: it weights by ``M·1``, which is zero at exactly these
+    vertices, so the gauge is set by the midpoints either way — and the reconstruction
+    is linear, so removing the mean before or after gives the same answer.
+    """
+    total = np.zeros_like(flux)
+    count = np.zeros(len(flux), dtype=np.int64)
+    for _order, nodes, _area in elements.values():
+        v0, v1, v2 = (global_index[k] for k in nodes[:3])
+        m01, m12, m20 = (global_index[k] for k in nodes[3:])
+        for vertex, near_a, near_b, far in (
+                (v0, m01, m20, m12),      # v0 lies on edges 01 and 20
+                (v1, m01, m12, m20),      # v1 lies on edges 01 and 12
+                (v2, m12, m20, m01)):     # v2 lies on edges 12 and 20
+            total[vertex] += flux[near_a] + flux[near_b] - flux[far]
+            count[vertex] += 1
+    rebuilt = np.array(flux, dtype=float, copy=True)
+    touched = count > 0
+    rebuilt[touched] = total[touched] / count[touched]
+    return rebuilt
+
+
 def _key(c, dim):
     return tuple(round(float(t), 9) for t in np.asarray(c).ravel()[:dim])
 
@@ -300,9 +338,10 @@ 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", "p1"):
-        raise ValueError("mass must be 'auto', 'lumped', 'consistent', or 'p1' "
-                         "(P1-projected recovery on a 3D P2 trace).")
+    if mass not in ("auto", "lumped", "consistent", "p1", "midpoint"):
+        raise ValueError("mass must be 'auto', 'lumped', 'consistent', 'p1' "
+                         "(P1-projected) or 'midpoint' (midpoint-reconstructed) "
+                         "— the last two apply to a 3D P2 trace.")
     if dim == 3:
         lsec = dm.getLocalSection()
         ncomp = lsec.getFieldComponents(0)
@@ -399,7 +438,34 @@ def coord(q):
                     )
                 order = orders.pop()
                 if mass == "auto":
-                    mass = "consistent" if order == 2 else "lumped"
+                    # A P2 TRIANGLE mass has vertex rows that sum to EXACTLY zero
+                    # (_P2_TRIANGLE_MASS, row sums [0,0,0,60,60,60]). Row sum i of a
+                    # mass matrix IS the integral of basis function i, because the
+                    # basis is a partition of unity:
+                    #     sum_j INT(phi_i phi_j) = INT(phi_i sum_j phi_j) = INT(phi_i)
+                    # so this says INT(phi_vertex) = 0 — a property of the P2 triangle
+                    # element, not of this code. A vertex DOF carries no mass, so an
+                    # L2 recovery fixes it by cancellation and amplifies whatever noise
+                    # is present: O(1), and independently of h. Measured on the Zhong
+                    # l=2 shell (#633), the consistent vertex error over cellSize
+                    # 0.25 -> 0.11 runs 0.076 / 0.074 / 0.070 / 0.119 / 0.093 — flat
+                    # and erratic, never converging.
+                    #
+                    # The midpoints carry all of the mass (row sums 60) and are
+                    # superconvergent (to 5e-4). So do not ask the recovery for vertex
+                    # values at all: keep the midpoints and RECONSTRUCT the vertices
+                    # from them ('midpoint'). Worst-node error against the analytic
+                    # coefficient beats the P1-projected recovery at every resolution
+                    # measured, on both boundaries, by 1.8x to 4.9x:
+                    #     h      0.25   0.20   0.16   0.13   0.11
+                    #     p1     .041   .026   .018   .013   .008   (Upper)
+                    #     mid    .016   .012   .012   .003   .004
+                    #     p1     .116   .067   .047   .030   .025   (Lower)
+                    #     mid    .094   .058   .043   .024   .015
+                    # 'p1' remains available and is sound; it simply discards the good
+                    # data along with the bad. The 2-D P2 LINE mass has vertex row sums
+                    # of 5 and needs none of this.
+                    mass = "midpoint" if order == 2 else "lumped"
                 if order == 2 and mass == "lumped":
                     raise ValueError(
                         "A 3D P2 triangular trace has zero row-sum mass at its "
@@ -436,6 +502,13 @@ def coord(q):
                         order = 1
                         mass = "lumped"
 
+                # 'midpoint' IS the consistent solve, plus a vertex reconstruction
+                # from its midpoint values afterwards.
+                reconstruct = mass == "midpoint"
+                if reconstruct:
+                    mass = "consistent" if order == 2 else "lumped"
+                    reconstruct = order == 2   # a P1 trace has no midpoints to use
+
                 keys = sorted(R_by)
                 global_index = {key: i for i, key in enumerate(keys)}
                 reaction = np.array([R_by[key] for key in keys], dtype=float)
@@ -485,6 +558,9 @@ def coord(q):
                         raise RuntimeError(
                             "Consistent boundary-mass solve failed on boundary " f"{boundary!r}."
                         )
+                    if reconstruct:
+                        flux = _rebuild_vertices_from_midpoints(
+                            elements, global_index, flux)
                 if remove_mean:
                     mean = float(np.dot(flux, boundary_mass) / np.sum(boundary_mass))
                     flux -= mean
@@ -519,8 +595,11 @@ def value_at(key):
         raise NotImplementedError(
             f"Boundary-flux recovery is not implemented for mesh dimension {dim}."
         )
-    if mass == "p1":
-        mass = "lumped"    # P1 consumers read vertex values; vertex lumping is sound
+    if mass in ("p1", "midpoint"):
+        # Both exist to work around the 3D P2 TRIANGLE's zero-mean vertex basis. A 2D
+        # P2 LINE trace has vertex row sums of 5, so its vertex values are sound as
+        # recovered and neither workaround is needed here.
+        mass = "lumped"
 
     lsec = dm.getLocalSection()
     ncomp = lsec.getFieldComponents(0)