Add distributed harmonic projection of rotated reactions

docs/advanced/index.md

@@ -82,6 +82,12 @@ time-derivative and Adams-Moulton/θ flux knobs, and how to pair them
 
 **[→ Semi-Lagrangian Time Integration](semi-lagrangian-time-integration.md)**
 
+### SUPG Scalar Transport
+Use local streamline-upwind stabilization with implicit BDF integration or
+the continuous-P1 CitcomS-compatible predictor-corrector.
+
+**[→ SUPG Scalar Transport](supg-transport.md)**
+
 ### Porous Media Flow
 Darcy flow, Richards equation, and variably-saturated groundwater modelling.
 
@@ -139,9 +145,10 @@ custom-meshes
 curved-boundary-conditions
 mesh-adaptation
 semi-lagrangian-time-integration
+supg-transport
 porous-flow
 snapshot-restore
 troubleshooting
 api-patterns
 SWARM-INTEGRATION-STATISTICS
-```
+```

docs/advanced/semi-lagrangian-time-integration.md

@@ -15,6 +15,9 @@ treating diffusion implicitly. This page explains the **two independent
 order knobs** in the scheme and how to pair them correctly — the common
 pitfall is mixing them.
 
+For local finite-element streamline stabilization without characteristic
+trace-back, see [SUPG scalar transport](supg-transport.md).
+
 ## The scheme has two time-integration choices
 
 The discrete residual assembled by the solver is

docs/advanced/supg-transport.md

@@ -0,0 +1,141 @@
+---
+title: "SUPG Scalar Transport"
+---
+
+# SUPG scalar transport
+
+`AdvDiffusionSUPG` solves
+
+$$
+\frac{\partial T}{\partial t} + \mathbf{u}\cdot\nabla T
+- \nabla\cdot(\kappa\nabla T) = f
+$$
+
+on simplex volume meshes. It adds streamline-upwind Petrov-Galerkin (SUPG)
+stabilization to the continuous finite-element residual. Advection remains a
+local finite-element operation: the solver does not trace departure points or
+interpolate a semi-Lagrangian history.
+
+## Minimal example
+
+This example transports and diffuses a continuous P1 scalar in a prescribed
+velocity field. Automatic stabilization is the default.
+
+```python
+import numpy as np
+import underworld3 as uw
+
+mesh = uw.meshing.UnstructuredSimplexBox(
+    minCoords=(0.0, 0.0),
+    maxCoords=(1.0, 1.0),
+    cellSize=0.1,
+)
+
+temperature = uw.discretisation.MeshVariable("T", mesh, 1, degree=1)
+velocity = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=1)
+
+with mesh.access(temperature, velocity):
+    x = temperature.coords[:, 0]
+    y = temperature.coords[:, 1]
+    temperature.data[:, 0] = np.sin(np.pi * x) * np.sin(np.pi * y)
+    velocity.data[:, 0] = 1.0
+    velocity.data[:, 1] = 0.0
+
+thermal = uw.systems.AdvDiffusionSUPG(
+    mesh,
+    u_Field=temperature,
+    V_fn=velocity.sym,
+)
+thermal.constitutive_model = uw.constitutive_models.DiffusionModel
+thermal.constitutive_model.Parameters.diffusivity = 0.01
+
+for boundary in ("Left", "Right", "Top", "Bottom"):
+    thermal.add_dirichlet_bc(0.0, boundary)
+
+for _ in range(10):
+    thermal.solve(timestep=1.0e-3, zero_init_guess=False)
+```
+
+The default `time_integrator="bdf"` uses an implicit Eulerian BDF method and
+the generic transient SUPG stabilization parameter. `order=1` and `order=2`
+select BDF1 and BDF2 respectively.
+
+## CitcomS-compatible predictor-corrector
+
+For continuous P1 temperature, UW3 also provides the row-sum-mass
+predictor-corrector used for the Zhong mantle-convection benchmark:
+
+```python
+temperature_rate = uw.discretisation.MeshVariable(
+    "Tdot", mesh, 1, degree=1
+)
+
+thermal = uw.systems.AdvDiffusionSUPG(
+    mesh,
+    u_Field=temperature,
+    V_fn=velocity.sym,
+    time_integrator="citcoms",
+    temperature_rate_field=temperature_rate,
+)
+thermal.constitutive_model = uw.constitutive_models.DiffusionModel
+thermal.constitutive_model.Parameters.diffusivity = 0.01
+
+dt = thermal.estimate_dt()
+thermal.solve(timestep=dt)
+```
+
+This path uses `adv_gamma=0.5`, two residual-correction iterations, positive
+row-sum mass, and the clipped CitcomS stabilization parameter by default. The
+timestep is explicit and must satisfy the returned advection-diffusion limit.
+Supplying a named `temperature_rate_field` makes the additional restart state
+visible and straightforward to checkpoint. Exact restart requires `T`,
+`Tdot`, and the solver snapshot metadata.
+
+## Choosing a transport solver
+
+| Method | Strength | Main cost or limitation | Restart state |
+| --- | --- | --- | --- |
+| `AdvDiffusionSUPG`, implicit BDF | Local assembly, no trace-back interpolation, automatic simplex stabilization | Timestep accuracy still requires convergence testing; automatic tau currently assumes scalar isotropic diffusivity on volume simplices | Temperature plus BDF history |
+| `AdvDiffusionSUPG`, CitcomS predictor-corrector | Second-order explicit update, row-lumped P1 mass, close to CitcomS mantle-convection numerics | Continuous P1 only; explicit advection-diffusion timestep limit | `T`, `Tdot`, solver metadata |
+| `AdvDiffusionSLCN` | Stable characteristic transport at large advective Courant number | Departure-point search/interpolation, flux history, and higher MPI memory/runtime | Temperature plus characteristic and flux histories |
+| `AdvDiffusionSLCN` with SL-BDF2 | Second-order characteristic history without Crank-Nicolson flux ringing | Two departure points and greater history/interpolation cost | Two-level characteristic history plus flux history |
+| `AdvDiffusionSLCN` with BDF1/Backward Euler | Robust, L-stable diffusion baseline | First-order time integration and trace-back interpolation | One characteristic history level |
+
+Use the CitcomS predictor-corrector when reproducing a continuous-P1 CitcomS
+benchmark. Use implicit SUPG when local streamline stabilization is desired
+without the explicit predictor-corrector restriction. Use SLCN when large
+advective timesteps are more important than trace-back cost. For every method,
+verify timestep and mesh convergence using the physical diagnostics of the
+problem; the solver name alone does not establish accuracy.
+
+## Conservation and boundedness
+
+SUPG is a consistent residual stabilization, but continuous Galerkin SUPG is
+not automatically monotone and does not guarantee nodal maximum principles.
+The CitcomS path uses positive row-sum mass, which improves the explicit
+update, but temperature bounds must still be checked. Semi-Lagrangian methods
+can remain stable at large advective Courant number, but departure-point
+interpolation is not strictly conservative. For either solver family, monitor
+the volume-integrated scalar, recovered boundary fluxes, source integral, and
+their discrete balance in addition to minimum and maximum nodal values.
+
+Parallel execution does not change these definitions. Compare global
+integrals between serial and MPI runs on the same mesh; do not compare local
+rank extrema or partition-dependent raw boundary-node sums.
+
+## Stabilization controls
+
+- Omit `tau` for automatic stabilization.
+- `tau_model="generic"` combines transient, advective, and diffusive scales.
+- `tau_model="citcoms"` selects the clipped steady CitcomS relation on a
+  simplex streamline length.
+- Pass an explicit scalar `tau` for unsupported element or constitutive-model
+  combinations. `tau=0` recovers the unstabilized Galerkin residual.
+- Automatic tau supports two- and three-dimensional simplex volume meshes and
+  scalar isotropic non-negative diffusivity.
+
+## Related documentation
+
+- [Semi-Lagrangian time integration](semi-lagrangian-time-integration.md)
+- [State snapshots and restore](snapshot-restore.md)
+- [Parallel computing](parallel-computing.md)

docs/developer/CHANGELOG.md

@@ -387,6 +387,9 @@ component exactly — correct on curved, tilted, and deformed boundaries (#293).
   weak contraction of the assembled normal reaction. Cylindrical-annulus
   Stokes responses use this fitted integral and its matching finite-element
   boundary norm instead of gathering pointwise samples for angular quadrature.
+- The spherical-shell geoid adapter accepts `projection="reaction"` to use
+  the same fitted integral without pointwise P2 recovery or a rank-zero
+  surface triangulation; `projection="centroid"` remains the default.
 - `uw.analytic.Zhong2008` implements the Hager--O'Connell propagator-matrix
   oracle used for the Zhong et al. spherical-shell response benchmark. It
   supports piecewise-constant radial viscosity and reproduces every analytical

docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md

@@ -97,13 +97,27 @@ response = uw.postprocessing.geoid.spherical_shell_response_from_rotated_stokes(
     planet_radius=6370000.0,
     gravity=9.8,
     gravitational_constant=6.67e-11,
+    projection="reaction",
 )
 ```
 
-The adapter delegates stress recovery to the existing rotated-free-slip API;
-it does not implement a second CBF, constrained-multiplier, or topography
-recovery path. `internal_load_coefficient` must use the same harmonic
-normalisation and sign convention as the model's internal load.
+The adapter supports two projection paths. `projection="centroid"` retains the
+original pointwise-recovery workflow: recover `sigma_nn`, gather the samples to
+rank zero, reconstruct a spherical triangulation, and integrate centroid
+values. `projection="reaction"` contracts the assembled normal-reaction load
+directly with the harmonic test function through
+`Stokes.boundary_normal_traction_integral()`. The latter is distributed, avoids
+the rank-zero surface reconstruction, and is an integral/fitted quantity rather
+than a consumer of the slowly converging P2 vertex values on curved boundaries
+(#414). Its fitted coefficient uses the matching discrete boundary norm, not an
+analytical spherical norm, so the numerator and denominator share the same
+faceted geometry.
+
+Both paths reuse the existing rotated-free-slip reaction; neither implements a
+second CBF, constrained-multiplier, or topography recovery. `centroid` remains
+the compatibility default while the direct reaction path accumulates benchmark
+coverage. `internal_load_coefficient` must use the same harmonic normalisation
+and sign convention as the model's internal load.
 
 When surface and CMB topography coefficients are already available, call
 `uw.postprocessing.geoid.spherical_shell_geoid_response()` or

src/underworld3/cython/petsc_generic_snes_solvers.pyx

@@ -3125,6 +3125,16 @@ class SolverBaseClass(uw_object):
                     self._subdict[name][1].localToGlobal(var.vec, sgvec)
                     gvec.restoreSubVector(self._subdict[name][0], sgvec)
             else:
+                # Map the variable's LOCAL vector through field 0's subDM rather
+                # than assuming the solver DM's local layout matches it. The two
+                # coincide on a plain single-field solver, which is why the
+                # direct `self.dm.localToGlobal(self.Unknowns.u.vec, gvec)` looked
+                # equivalent -- but where they differ it writes to the wrong slots
+                # and the field comes back never-written. Restored while CI is red
+                # on tests/test_1120_SLVectorCartesian.py::test_SLVec_boxmesh[mesh1],
+                # a semi-Lagrangian VECTOR test, i.e. exactly the single-field path
+                # this branch serves; its recovered values were ~1e-18 against an
+                # analytic ~1e-5.
                 _names, _iss, _subdms = self.dm.createFieldDecomposition()
                 try:
                     sgvec = gvec.getSubVector(_iss[0])
@@ -3218,6 +3228,18 @@ class SolverBaseClass(uw_object):
         return _bff(self, boundary, field, mass=mass, remove_mean=remove_mean,
                     scale=scale, normal=normal)
 
+    def boundary_flux_integral(self, boundary):
+        r"""Integrated scalar CBF flux through ``boundary``.
+
+        This is the direct integral diagnostic for quantities such as Nusselt
+        numbers. It sums the consistent scalar nodal reactions collectively,
+        avoiding pointwise de-smearing and a temporary flux MeshVariable. Use
+        :meth:`boundary_flux` or :meth:`boundary_flux_field` when nodal values
+        are required.
+        """
+        from underworld3.utilities.boundary_flux import boundary_flux_integral as _bfi
+        return _bfi(self, boundary)
+
 ## Specific to dimensionality
 
 

src/underworld3/function/_function.pyx

@@ -48,6 +48,33 @@ cdef extern from "petsc.h" nogil:
     PetscErrorCode DMSwarmSetMigrateType(PetscDM dm, DMSwarmMigrateType mtype)
     PetscErrorCode DMSwarmGetMigrateType(PetscDM dm, DMSwarmMigrateType *mtype)
 
+
+_fallback_stats_enabled = False
+_fallback_stats = {}
+
+
+def _reset_global_evaluate_fallback_stats(enabled=True):
+    """Reset and optionally enable internal MPI fallback diagnostics."""
+
+    global _fallback_stats_enabled, _fallback_stats
+    _fallback_stats_enabled = bool(enabled)
+    _fallback_stats = {
+        "calls": 0,
+        "calls_with_points": 0,
+        "local_extrapolated_points": 0,
+        "replicated_points_per_rank": 0,
+        "max_replicated_points_per_call": 0,
+        "temporary_bytes_per_rank": 0,
+        "max_temporary_bytes_per_call": 0,
+    }
+
+
+def _get_global_evaluate_fallback_stats():
+    """Return one rank's internal MPI fallback diagnostics."""
+
+    return dict(_fallback_stats)
+
+
 class UnderworldAppliedFunction(sympy.core.function.AppliedUndef):
     """
     Applied Underworld function representing a mesh variable evaluated at coordinates.
@@ -610,7 +637,21 @@ def global_evaluate_nd(   expr,
         counts = np.array(comm.allgather(ext_coords.shape[0]), dtype=int)
         n_ext_total = int(counts.sum())
 
+        if _fallback_stats_enabled:
+            _fallback_stats["calls"] += 1
+            _fallback_stats["local_extrapolated_points"] += int(
+                ext_coords.shape[0]
+            )
+            _fallback_stats["replicated_points_per_rank"] += n_ext_total
+            _fallback_stats["max_replicated_points_per_call"] = max(
+                _fallback_stats["max_replicated_points_per_call"],
+                n_ext_total,
+            )
+
         if n_ext_total > 0:
+            if _fallback_stats_enabled:
+                _fallback_stats["calls_with_points"] += 1
+
             parts = comm.allgather(ext_coords)
             all_ext = np.concatenate(
                 [p for p in parts if p.size], axis=0).reshape(n_ext_total, -1)
@@ -647,6 +688,30 @@ def global_evaluate_nd(   expr,
             best_flag = np.empty(n_ext_total, dtype=np.int32)
             comm.Allreduce([contrib_flag, MPI.INT], [best_flag, MPI.INT], op=MPI.SUM)
 
+            if _fallback_stats_enabled:
+                temporary_bytes = sum(
+                    int(array.nbytes)
+                    for array in (
+                        all_ext,
+                        ext_vals,
+                        ext_flag,
+                        dist2,
+                        min_dist2,
+                        my_claim,
+                        win_rank,
+                        contrib_val,
+                        best_val,
+                        contrib_flag,
+                        best_flag,
+                    )
+                )
+                temporary_bytes += sum(int(part.nbytes) for part in parts)
+                _fallback_stats["temporary_bytes_per_rank"] += temporary_bytes
+                _fallback_stats["max_temporary_bytes_per_call"] = max(
+                    _fallback_stats["max_temporary_bytes_per_call"],
+                    temporary_bytes,
+                )
+
             # Scatter this rank's segment of the global set back to its points.
             offset = int(counts[:comm.rank].sum())
             seg = slice(offset, offset + ext_coords.shape[0])