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

@@ -278,6 +278,12 @@ component exactly — correct on curved, tilted, and deformed boundaries (#293).
   existing boundary traction onto an axisymmetric harmonic; the pure functions
   also accept coefficients recovered by other methods and an optional internal
   load.
+- Rotated free slip now exposes
+  `Stokes.boundary_normal_traction_integral(boundary, fn)` for a distributed
+  weak contraction of the assembled normal reaction. The spherical-shell geoid
+  adapter accepts `projection="reaction"` to use this fitted quantity without
+  pointwise P2 recovery or a rank-zero surface triangulation; the existing
+  `projection="centroid"` behavior 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
@@ -137,10 +151,11 @@ Published reference solvers, such as the Zhong et al. propagator-matrix method,
 belong in `uw.analytic`; their computed topography coefficients can be passed to
 the pure post-processing functions above.
 
-The rotated harmonic projector gathers boundary samples to rank zero and
-reconstructs their spherical triangulation. A future boundary-reaction
-functional could replace this step with a direct distributed finite-element
-projection without changing the coefficient API.
+`Stokes.boundary_normal_traction_integral(boundary, fn)` is useful beyond the
+geoid adapter whenever only an integrated or fitted normal-traction diagnostic
+is required. Pointwise consumers should continue to call
+`boundary_normal_traction()` or `dynamic_topography()` and follow their curved
+P2 midpoint/field guidance.
 
 ## See also
 

src/underworld3/checkpoint/disk_snapshot.py

@@ -461,6 +461,15 @@ def read_snapshot(model, path: str) -> None:
     import h5py
 
     md = read_snapshot_metadata(path)
+    write_size = int(md.get("mpi_ranks_at_write", 1))
+    if write_size != int(uw.mpi.size):
+        raise ValueError(
+            f"snapshot at {path} was written on {write_size} MPI rank(s); "
+            f"this run uses {uw.mpi.size}. Exact disk restart requires the "
+            "same rank count; use mesh.write_timestep/read_timestep for "
+            "coordinate-remapped field transfer."
+        )
+
     bulk_dir = _bulk_dir_for(path)
     if not os.path.isdir(bulk_dir):
         raise FileNotFoundError(

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()
                 sgvec = gvec.getSubVector(_iss[0])
                 _subdms[0].localToGlobal(self.Unknowns.u.vec, sgvec)
@@ -3209,6 +3219,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
 
 
@@ -6502,6 +6524,25 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         return _dtf(self, boundary, self._rotated_freeslip_info, field,
                     buoyancy_scale=buoyancy_scale, mass=mass)
 
+    def boundary_normal_traction_integral(self, boundary, fn, remove_mean=True):
+        r"""Return the boundary integral of ``sigma_nn * fn`` directly from the
+        rotated-free-slip constraint reaction.
+
+        With ``remove_mean=True`` (default), the constant normal-traction gauge
+        is removed before projection.  Unlike pointwise
+        :meth:`boundary_normal_traction`, this weak projection does not recover
+        nodal traction values or gather a global surface mesh.  It is therefore
+        suitable for harmonic/integral diagnostics on curved P2 boundaries,
+        whose recovered vertex values converge slowly (issue #414).
+        """
+        if self._rotated_freeslip_info is None:
+            raise RuntimeError(
+                "boundary_normal_traction_integral requires a completed "
+                "rotated-free-slip solve.")
+        from underworld3.utilities.rotated_bc import boundary_normal_traction_integral as _bnti
+        return _bnti(self, boundary, self._rotated_freeslip_info, fn,
+                     remove_mean=remove_mean)
+
     def add_nitsche_bc(self, conds=None, boundary=None, direction=None, normal=None,
                        gamma=10.0, theta=1, mask=None, local_h=True, g=None):
         r"""Add Nitsche weak enforcement of a velocity constraint along a direction.

src/underworld3/discretisation/discretisation_mesh.py

@@ -4836,9 +4836,10 @@ def write_timestep(
         - ``create_xdmf=True`` writes ParaView/XDMF output. Variable files also
           receive ``/vertex_fields`` or ``/cell_fields`` compatibility groups,
           and rank 0 writes the companion ``.xdmf`` file.
-        - ``petsc_reload=True`` writes PETSc DMPlex section/vector metadata into
-          the same per-variable HDF5 files. These files can then be loaded with
-          ``MeshVariable.read_checkpoint()`` for PETSc-native same-mesh reload.
+        - ``petsc_reload=True`` writes PETSc DMPlex metadata and an owned
+          global-vector payload into the same per-variable HDF5 files. These
+          files can then be loaded with ``MeshVariable.read_checkpoint()`` for
+          exact same-layout reload.
 
         Common choices are:
 
@@ -5015,7 +5016,7 @@ def _write_petsc_reload_variable(self, viewer, var):
             subdm.destroy()
 
     def _write_petsc_reload_file(self, checkpoint_file, variables, mode="w"):
-        """Write PETSc DMPlex section/vector reload metadata."""
+        """Write DMPlex metadata and exact same-layout global vectors."""
 
         old_dm_name = self.dm.getName()
         self.dm.setName("uw_mesh")
@@ -5037,6 +5038,23 @@ def _write_petsc_reload_file(self, checkpoint_file, variables, mode="w"):
             if old_dm_name is not None:
                 self.dm.setName(old_dm_name)
 
+        viewer = PETSc.ViewerHDF5().create(
+            checkpoint_file, "a", comm=PETSc.COMM_WORLD
+        )
+        try:
+            viewer.pushGroup("/uw_checkpoint")
+            for var in variables:
+                var._sync_lvec_to_gvec()
+                checkpoint_vec = PETSc.Vec().createWithArray(
+                    var._gvec.array_r, comm=PETSc.COMM_WORLD
+                )
+                checkpoint_vec.setName(var.clean_name)
+                viewer(checkpoint_vec)
+                checkpoint_vec.destroy()
+            viewer.popGroup()
+        finally:
+            viewer.destroy()
+
     @timing.routine_timer_decorator
     def write_checkpoint(
         self,