Fix 3D boundary reaction flux recovery (#404)

* Fix 3D boundary reaction flux recovery

Recover pointwise traction from 3D rotated free-slip reactions by assembling a lumped triangular boundary mass and dividing nodal reactions by that mass. Support P1 traces directly and P2 traces through four-triangle facet subdivision, with global coordinate-based assembly for partition independence and a surface-mass-weighted mean gauge.

Add serial and MPI regressions using the Zhong et al. l=2 internal-load shell. The corrected surface and CMB topography coefficients reproduce the reference response and remain invariant across MPI partitions. Retain focused 2D traction and flux coverage.

* Fix pointwise P2 surface reaction recovery

Replace the invalid four-subtriangle P2 lumping with the exact six-node triangular surface mass matrix and a sparse consistent solve. Add mass=auto so existing 2D and 3D P1 traces retain lumped recovery while 3D P2 selects the mathematically required consistent operator; reject explicit P2 lumping and unsupported 3D traces.

Add pointwise constant-flux P1/P2 tetrahedral-box regressions in serial and MPI, including two- and four-rank partitions. Extend the Zhong l=2 validation to compare all nodes, vertices, and edge midpoints independently, removing the prior all-node cancellation blind spot.

Document trace-order behavior, non-triangular limitations, analytic-normal partition guarantees, and the corrected 3D CBF behavior.

* Parameterize Zhong topography diagnostic resolution

Allow the serial and MPI Zhong l=2 rotated free-slip diagnostics to accept an optional shell cell size while retaining 0.25 as the regression-test default. Extend the standalone MPI diagnostic output with the selected cell size so manual convergence results are self-describing. This supports serial/MPI comparisons at 1/4, 1/8, and 1/16 resolution without adding costly high-resolution cases to routine CI.

Author: gthyagi

docs/developer/CHANGELOG.md

@@ -95,6 +95,13 @@ component exactly — correct on curved, tilted, and deformed boundaries (#293).
 - A general **consistent boundary flux (CBF) primitive** recovers boundary
   fluxes for any solver — surface heat flux / Nusselt number for scalar
   diffusion, boundary traction σ·n for Stokes (#294).
+- Three-dimensional CBF recovery now assembles the exact triangular trace mass:
+  P1 supports lumped or consistent recovery, while P2 uses the required
+  consistent six-node surface-mass solve. The default `mass="auto"` selects the
+  valid method; explicit P2 lumping and non-triangular 3D traces raise instead
+  of returning a non-pointwise reaction scaling. Strict MPI invariance of a
+  vector normal projection requires an analytic normal; geometric facet-normal
+  seam sensitivity is unchanged (#404).
 - Recorded as the preferred free-slip BC in the project guidance (#300);
   conda PETSc floor raised to ≥ 3.25 for FMG/rotation API consistency (#304).
 

src/underworld3/cython/petsc_generic_snes_solvers.pyx

@@ -2304,7 +2304,7 @@ class SolverBaseClass(uw_object):
             self.dm.restoreLocalVec(xlocal)
             self.dm.restoreGlobalVec(gvec)
 
-    def boundary_flux(self, boundary, mass="lumped", remove_mean=False, normal=None):
+    def boundary_flux(self, boundary, mass="auto", remove_mean=False, normal=None):
         r"""Consistent boundary flux on ``boundary``, recovered from the essential-BC
         reaction of the last solve (the Consistent Boundary Flux method).
 
@@ -2314,15 +2314,21 @@ 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 the ``"lumped"`` (diagonal, monotone —
-        no overshoot at a flux jump) or ``"consistent"`` boundary mass. ``remove_mean``
-        subtracts the boundary mean — leave ``False`` for a physical flux (the mean is
-        the Nusselt number); ``True`` gives a gauge-free field (e.g. dynamic topography).
-        Parallel-safe and partition-independent."""
+        ``mass`` de-smears the nodal reaction with ``"lumped"`` or ``"consistent"``
+        boundary mass. ``"auto"`` (default) selects lumped recovery for 2D traces and
+        3D P1 triangles, and the required consistent solve for 3D P2 triangles.
+        ``remove_mean`` subtracts the boundary mean — leave ``False`` for a physical
+        flux (the mean is the Nusselt number); ``True`` gives a gauge-free field.
+
+        Three-dimensional recovery supports triangular P1/P2 traces; quadrilateral
+        traces raise explicitly. Reaction and mass assembly are partition-independent.
+        For vector fluxes, supply an analytic ``normal`` when strict partition
+        independence of the normal projection is required; geometric facet-normal
+        averaging at partition seams has a small pre-existing partition sensitivity."""
         from underworld3.utilities.boundary_flux import boundary_flux as _bf
         return _bf(self, boundary, mass=mass, remove_mean=remove_mean, normal=normal)
 
-    def boundary_flux_field(self, boundary, field, mass="lumped",
+    def boundary_flux_field(self, boundary, field, mass="auto",
                             remove_mean=False, scale=1.0, normal=None):
         r"""Write the consistent boundary flux (see :meth:`boundary_flux`) onto a scalar
         MeshVariable ``field`` at the boundary nodes (interior untouched), multiplied by
@@ -5374,26 +5380,26 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
             J1.destroy(); J2.destroy()
         return rel > tol
 
-    def boundary_normal_traction(self, boundary, mass="lumped"):
+    def boundary_normal_traction(self, boundary, mass="auto"):
         r"""Return the boundary normal traction :math:`\sigma_{nn}` on a
         rotated-free-slip ``boundary`` as the constraint reaction from the last
         solve — the smooth, bounded quantity used for dynamic topography
         (:math:`h_\infty=-(\sigma_{nn}-\overline{\sigma_{nn}})/\rho g`). Requires a
         prior :meth:`add_rotated_freeslip_bc` on ``boundary`` and a completed
         :meth:`solve`.
 
-        ``mass`` chooses the boundary-mass de-smear of the nodal reaction:
-        ``"lumped"`` (default) is monotone — it cannot overshoot where the traction
-        jumps (e.g. across a viscosity contrast), so it is the safe choice for driving
-        a free surface; ``"consistent"`` uses the full P2 line mass (marginally sharper
-        on smooth tractions, but overshoots at discontinuities)."""
+        ``mass="auto"`` (default) uses lumped recovery for 2D traces and 3D P1
+        triangles, and the required 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."""
         if self._rotated_freeslip_info is None:
             raise RuntimeError(
                 "boundary_normal_traction requires a completed rotated-free-slip solve.")
         from underworld3.utilities.rotated_bc import boundary_normal_traction as _bnt
         return _bnt(self, boundary, self._rotated_freeslip_info, mass=mass)
 
-    def dynamic_topography(self, boundary, field, buoyancy_scale=1.0, mass="lumped"):
+    def dynamic_topography(self, boundary, field, buoyancy_scale=1.0, mass="auto"):
         r"""Write the dynamic topography
         :math:`h = -(\sigma_{nn}-\overline{\sigma_{nn}})/(\Delta\rho\,g)` on a
         rotated-free-slip ``boundary`` onto a scalar MeshVariable ``field``, from the
@@ -5403,9 +5409,9 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         pass it here after each :meth:`solve`; its boundary nodes are filled and the
         interior left untouched.
 
-        ``buoyancy_scale`` is :math:`\Delta\rho\,g` (traction → length). ``mass`` selects
-        the recovery de-smear (``"lumped"`` default is monotone — no overshoot at a
-        stress jump — and is the safe choice for a free surface). Requires a prior
+        ``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
         :meth:`add_rotated_freeslip_bc` on ``boundary`` and a completed :meth:`solve`."""
         if self._rotated_freeslip_info is None:
             raise RuntimeError(

src/underworld3/utilities/boundary_flux.py

@@ -15,15 +15,40 @@
 assembled here in ``_desmear`` by SUMMING each rank's partial contribution by coordinate
 (the same rock-solid gather used for the boundary mass — no hand-rolled global assembly).
 
-``mass="lumped"`` (default) uses the diagonal boundary mass: being an M-matrix it cannot
-overshoot where the flux jumps (no Gibbs wiggle) and is a purely local division.
+``mass="auto"`` (default) uses a diagonal lumped mass where the trace basis admits
+positive row sums (the 2D P2 line trace and the 3D P1 triangle trace), and the consistent
+mass otherwise. A 3D P2 triangle has exactly zero row sum at every vertex, so its
+pointwise recovery requires the consistent surface-mass solve.
 ``remove_mean=False`` (default) keeps the physical mean flux (the Nusselt number);
 set ``remove_mean=True`` for a gauge-free field (e.g. dynamic topography).
 """
 import numpy as np
 from mpi4py import MPI
 
 
+# M_e = (area / 12) * _P1_TRIANGLE_MASS.
+_P1_TRIANGLE_MASS = np.array(
+    (
+        (2.0, 1.0, 1.0),
+        (1.0, 2.0, 1.0),
+        (1.0, 1.0, 2.0),
+    )
+)
+
+# M_e = (area / 180) * _P2_TRIANGLE_MASS. Node order: vertices
+# (0, 1, 2), then edge nodes (01, 12, 20).
+_P2_TRIANGLE_MASS = np.array(
+    (
+        (6.0, -1.0, -1.0, 0.0, -4.0, 0.0),
+        (-1.0, 6.0, -1.0, 0.0, 0.0, -4.0),
+        (-1.0, -1.0, 6.0, -4.0, 0.0, 0.0),
+        (0.0, 0.0, -4.0, 32.0, 16.0, 16.0),
+        (-4.0, 0.0, 0.0, 16.0, 32.0, 16.0),
+        (0.0, -4.0, 0.0, 16.0, 16.0, 32.0),
+    )
+)
+
+
 def _key(c, dim):
     return tuple(round(float(t), 9) for t in np.asarray(c).ravel()[:dim])
 
@@ -140,12 +165,160 @@ 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 dim == 3:
+        lsec = dm.getLocalSection()
+        ncomp = lsec.getFieldComponents(0)
+        f0, f1 = dm.getHeightStratum(1)
+        e0, e1 = dm.getDepthStratum(1)
+
+        def coord(q):
+            return _point_coord(dm, dim, cvec, csec, v0, v1, q)
+
+        nodeR = {_key(x, dim): float(r) for x, r in zip(xs, R)}
+        sis = _boundary_stratum_is(dm, solver.mesh, boundary)
+        facets = [] if not (sis and sis.getSize() > 0) else [
+            int(q) for q in sis.getIndices() if f0 <= int(q) < f1
+        ]
+        local_elements = []
+
+        for facet in facets:
+            closure = [int(q) for q in dm.getTransitiveClosure(facet)[0]]
+            vertices = [q for q in closure if v0 <= q < v1]
+            edges = [q for q in closure if e0 <= q < e1]
+            if len(vertices) != 3:
+                raise NotImplementedError(
+                    "3D boundary-flux recovery currently requires triangular facets."
+                )
+            if lsec.getFieldDof(facet, 0) > 0:
+                raise NotImplementedError(
+                    "3D boundary-flux recovery supports P1 or P2 triangular traces."
+                )
+
+            vertex_coords = [np.asarray(coord(q), dtype=float) for q in vertices]
+            vertex_keys = [_key(value, dim) for value in vertex_coords]
+            edge_midpoints = {}
+            for edge in edges:
+                edge_dof = lsec.getFieldDof(edge, 0)
+                if edge_dof <= 0:
+                    continue
+                if edge_dof != ncomp:
+                    raise NotImplementedError(
+                        "3D boundary-flux recovery supports P1 or P2 triangular traces."
+                    )
+                edge_vertices = [
+                    int(q)
+                    for q in dm.getTransitiveClosure(edge)[0]
+                    if v0 <= int(q) < v1
+                ]
+                if len(edge_vertices) == 2:
+                    edge_key = frozenset(_key(coord(q), dim) for q in edge_vertices)
+                    edge_midpoints[edge_key] = _key(coord(edge), dim)
+
+            a, b, c = vertex_coords
+            area = 0.5 * float(np.linalg.norm(np.cross(b - a, c - a)))
+            if not edge_midpoints:
+                local_elements.append((1, tuple(vertex_keys), area))
+                continue
+            if len(edge_midpoints) != 3:
+                raise NotImplementedError(
+                    "3D boundary-flux recovery supports P1 or complete P2 triangular traces."
+                )
+
+            m01 = edge_midpoints[frozenset((vertex_keys[0], vertex_keys[1]))]
+            m12 = edge_midpoints[frozenset((vertex_keys[1], vertex_keys[2]))]
+            m20 = edge_midpoints[frozenset((vertex_keys[2], vertex_keys[0]))]
+            local_elements.append(
+                (
+                    2,
+                    tuple(vertex_keys) + (m01, m12, m20),
+                    area,
+                )
+            )
+
+        R_by = {}
+        for rank_values in comm.allgather(nodeR):
+            for key, value in rank_values.items():
+                R_by[key] = (
+                    R_by.get(key, 0.0) + value if partial_reaction else value
+                )
+
+        elements = {}
+        for rank_elements in comm.allgather(local_elements):
+            for order, nodes, area in rank_elements:
+                elements[(order, tuple(sorted(nodes)))] = (order, nodes, area)
+
+        orders = {order for order, _nodes, _area in elements.values()}
+        if len(orders) != 1:
+            raise RuntimeError(
+                f"Expected one trace order on boundary {boundary!r}, found {sorted(orders)}."
+            )
+        order = orders.pop()
+        if mass == "auto":
+            mass = "consistent" 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 vertices; "
+                "use mass='consistent' for pointwise boundary-flux recovery."
+            )
+
+        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)
+        if mass == "lumped":
+            boundary_mass = np.zeros(len(keys), dtype=float)
+            for _order, nodes, area in elements.values():
+                for key in nodes:
+                    boundary_mass[global_index[key]] += area / 3.0
+            missing = np.flatnonzero(boundary_mass <= 0.0)
+            if missing.size:
+                raise RuntimeError(
+                    f"Boundary mass is zero at {missing.size} nodes on {boundary!r}."
+                )
+            flux = reaction / boundary_mass
+        elif mass == "consistent":
+            from scipy.sparse import coo_matrix
+            from scipy.sparse.linalg import spsolve
+
+            rows = []
+            cols = []
+            values = []
+            reference_mass = (
+                _P1_TRIANGLE_MASS if order == 1 else _P2_TRIANGLE_MASS
+            )
+            mass_scale = 12.0 if order == 1 else 180.0
+            for _order, nodes, area in elements.values():
+                indices = [global_index[key] for key in nodes]
+                element_mass = (area / mass_scale) * reference_mass
+                for i, row in enumerate(indices):
+                    for j, col in enumerate(indices):
+                        rows.append(row)
+                        cols.append(col)
+                        values.append(element_mass[i, j])
+            surface_mass = coo_matrix(
+                (values, (rows, cols)), shape=(len(keys), len(keys))
+            ).tocsr()
+            surface_mass.sum_duplicates()
+            flux = np.asarray(spsolve(surface_mass, reaction), dtype=float)
+            boundary_mass = np.asarray(
+                surface_mass @ np.ones(len(keys), dtype=float)
+            )
+            if not np.all(np.isfinite(flux)):
+                raise RuntimeError(
+                    f"Consistent boundary-mass solve failed on boundary {boundary!r}."
+                )
+        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])
+
     if dim != 2:
-        # no line-mass geometry yet in 3D → global-mean lumped fallback
-        tot = comm.allreduce(float(np.sum(R)), op=MPI.SUM)
-        cnt = comm.allreduce(int(len(R)), op=MPI.SUM)
-        m = tot / max(cnt, 1)
-        return np.asarray(R) - (m if remove_mean else 0.0)
+        raise NotImplementedError(
+            f"Boundary-flux recovery is not implemented for mesh dimension {dim}."
+        )
+    if mass == "auto":
+        mass = "lumped"
 
     e0, e1 = dm.getDepthStratum(1)
     def vcoord(q): return cvec[csec.getOffset(q) // dim]
@@ -198,7 +371,7 @@ def vcoord(q): return cvec[csec.getOffset(q) // dim]
     return np.array([sig[gi[_key(x, dim)]] for x in xs])
 
 
-def boundary_flux(solver, boundary, mass="lumped", remove_mean=False, normal=None):
+def boundary_flux(solver, boundary, mass="auto", remove_mean=False, normal=None):
     """See ``SolverBaseClass.boundary_flux``. Returns ``(xs, flux)`` for this rank's
     boundary nodes; scalar solver → normal flux, vector solver → traction (or its normal
     component if ``normal`` is given)."""
@@ -253,7 +426,7 @@ def write_boundary_scalar_field(solver, field, value_by_key, dim):
     return field
 
 
-def boundary_flux_field(solver, boundary, field, mass="lumped",
+def boundary_flux_field(solver, boundary, field, mass="auto",
                         remove_mean=False, scale=1.0, normal=None):
     r"""See ``SolverBaseClass.boundary_flux_field`` (the documented entry point;
     this free function is its implementation and shares its name). Writes

src/underworld3/utilities/rotated_bc.py

@@ -1042,7 +1042,7 @@ def _mode_satisfies_constraints(solver, Q, normal_rows, tg, tol=1e-8):
     return viol < tol
 
 
-def boundary_normal_traction(solver, boundary, solve_result, mass="lumped"):
+def boundary_normal_traction(solver, boundary, solve_result, mass="auto"):
     """Boundary normal traction σ_nn on `boundary` from the constraint reaction of the
     last rotated-free-slip solve (``solve_result`` is the dict returned by
     ``solve_rotated_freeslip`` / ``solve_rotated_freeslip_nonlinear`` — see their
@@ -1055,22 +1055,19 @@ def boundary_normal_traction(solver, boundary, solve_result, mass="lumped"):
     rotated frame's normal row) is corner-correct — at a node shared with another
     rotated-free-slip boundary the rotated frame's first row is a mix of both walls'
     normals, but n̂·r_c is the true normal traction for this boundary. The pointwise
-    σ_nn is the boundary-mass de-smear of R (2D).
+    σ_nn is the boundary-mass de-smear of R.
 
     ``mass`` selects the de-smear:
-      * ``"lumped"`` (default) — the diagonal (row-sum) boundary mass. Being an M-matrix
-        it CANNOT overshoot at a stress discontinuity (no Gibbs wiggle where the traction
-        jumps, e.g. across a viscosity contrast), it is a purely local division (no global
-        mass solve → trivially parallel), and it is marginally more accurate than the
-        consistent mass on SolCx. Recommended for driving a free surface, where an
-        overshoot at a sharp feature injects a spurious surface-velocity pulse.
-      * ``"consistent"`` — the full consistent P2 line mass. Marginally sharper on smooth
-        tractions but overshoots at discontinuities.
+      * ``"auto"`` (default) — lumped for 2D traces and 3D P1 triangles, consistent for
+        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.
 
     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
     allgather of the boundary elements, so every rank produces the same de-smear and the
-    mean-removal gauge is global.
+    mean-removal gauge is global. In 3D, only triangular P1/P2 traces are supported.
     """
     dm = solver.dm
     dim = solver.mesh.dim
@@ -1119,12 +1116,12 @@ def boundary_normal_traction(solver, boundary, solve_result, mass="lumped"):
 
 
 def dynamic_topography_field(solver, boundary, solve_result, field,
-                             buoyancy_scale=1.0, mass="lumped"):
+                             buoyancy_scale=1.0, mass="auto"):
     """Populate a scalar MeshVariable ``field`` with the dynamic topography
     :math:`h = -(\\sigma_{nn}-\\overline{\\sigma_{nn}})/(\\Delta\\rho\\,g)` on ``boundary``,
-    recovered from the rotated-free-slip constraint reaction (lumped by default —
-    monotone, no Gibbs overshoot at a stress jump). Interior nodes are left untouched.
-    Returns ``field``.
+    recovered from the rotated-free-slip constraint reaction. ``mass="auto"`` uses
+    lumped recovery where valid and the consistent surface mass for 3D P2 triangles.
+    Interior nodes are left untouched. Returns ``field``.
 
     This is the hand-off to the free-surface machinery: the 3-number topography
     integrator drives node motion from a surface field, so σ_nn is written onto the