Add reaction-based spherical-shell geoid projection

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/postprocessing/geoid.py

@@ -892,11 +892,36 @@ def _rotated_topography_coefficient(
     harmonic_degree: int,
     buoyancy_scale: float,
     response_sign: float,
+    projection: str,
 ) -> float:
     buoyancy_scale = float(buoyancy_scale)
     if not np.isfinite(buoyancy_scale) or buoyancy_scale == 0.0:
         raise ValueError("Boundary buoyancy scales must be finite and nonzero.")
 
+    if projection == "reaction":
+        import sympy
+        from underworld3.maths import BdIntegral
+
+        theta = stokes.mesh.CoordinateSystem.xR[1]
+        harmonic = sympy.assoc_legendre(harmonic_degree, 0, sympy.cos(theta))
+        traction_integral = stokes.boundary_normal_traction_integral(
+            boundary,
+            harmonic,
+            remove_mean=True,
+        )
+        # The reaction is the load functional assembled on the faceted FE
+        # boundary.  Use the matching discrete inner product for the fitted
+        # coefficient; an analytical spherical norm would mix geometries and
+        # introduce a chord-area bias, especially on the smaller CMB.
+        harmonic_norm = float(
+            BdIntegral(stokes.mesh, fn=harmonic**2, boundary=boundary).evaluate()
+        )
+        return float(
+            -response_sign * traction_integral / (buoyancy_scale * harmonic_norm)
+        )
+    if projection != "centroid":
+        raise ValueError("projection must be 'centroid' or 'reaction'.")
+
     coords, sigma_nn = stokes.boundary_normal_traction(boundary, mass="auto")
     local_rows = np.column_stack(
         (
@@ -958,16 +983,18 @@ def spherical_shell_response_from_rotated_stokes(
     planet_radius: float | None = None,
     gravity: float | None = None,
     gravitational_constant: float = 6.67430e-11,
+    projection: str = "centroid",
 ) -> SphericalShellResponse:
     r"""Compute spherical-shell response from a rotated-free-slip Stokes solve.
 
-    Normal traction recovery is delegated to the existing
-    :meth:`Stokes.boundary_normal_traction` implementation. This adapter only
-    projects the two boundary responses onto the unnormalised
-    axisymmetric :math:`P_l^0` harmonic.  Use the pure coefficient functions
-    directly for other harmonic orders or topography-recovery methods. Density
-    contrasts, planet radius, and gravity are required when
-    ``include_self_gravity`` is true.
+    ``projection="centroid"`` (default) recovers pointwise traction and fits it
+    over a triangulation of the boundary samples. ``projection="reaction"``
+    contracts the assembled nodal reaction directly with the harmonic test
+    function. The latter is a distributed weak/integral quantity that avoids
+    consuming slowly converging P2 vertex values on curved boundaries (issue
+    #414). Use the pure coefficient functions directly for other harmonic
+    orders or topography-recovery methods. Density contrasts, planet radius,
+    and gravity are required when ``include_self_gravity`` is true.
     """
 
     ri, ro, degree = _validate_geometry(
@@ -977,6 +1004,8 @@ def spherical_shell_response_from_rotated_stokes(
     )
     if not isinstance(include_self_gravity, bool):
         raise TypeError("include_self_gravity must be True or False.")
+    if projection not in ("centroid", "reaction"):
+        raise ValueError("projection must be 'centroid' or 'reaction'.")
     if degree == 0:
         raise ValueError(
             "The rotated-Stokes adapter requires harmonic_degree >= 1 because "
@@ -1005,6 +1034,7 @@ def spherical_shell_response_from_rotated_stokes(
         degree,
         surface_buoyancy_scale,
         1.0,
+        projection,
     )
     cmb_topography = _rotated_topography_coefficient(
         stokes,
@@ -1013,6 +1043,7 @@ def spherical_shell_response_from_rotated_stokes(
         degree,
         cmb_buoyancy_scale,
         -1.0,
+        projection,
     )
     geoid = spherical_shell_geoid_response(
         radius_inner=ri,

tests/parallel/test_1071_spherical_shell_geoid_parallel.py

@@ -56,6 +56,7 @@ def test_rotated_spherical_shell_geoid_matches_serial_reference():
         planet_radius=6370000.0,
         gravity=9.8,
         gravitational_constant=6.67e-11,
+        projection="reaction",
     )
     values = np.array(
         [
@@ -71,7 +72,7 @@ def test_rotated_spherical_shell_geoid_matches_serial_reference():
     )
 
     for rank_values in uw.mpi.comm.allgather(values):
-        assert np.array_equal(rank_values, values)
+        np.testing.assert_allclose(rank_values, values, rtol=1.0e-13, atol=1.0e-14)
 
     zhong_table_2 = np.array(
         [0.41920, 0.77060, 0.02579, 0.03206, 0.49980, 0.93130, 0.04486, 0.05461]

tests/test_1070_postprocessing_geoid.py

@@ -182,6 +182,17 @@ def test_rotated_adapter_requires_explicit_self_gravity_parameters():
         )
 
 
+def test_rotated_adapter_rejects_unknown_projection():
+    with pytest.raises(ValueError, match="projection must be"):
+        uw.postprocessing.geoid.spherical_shell_response_from_rotated_stokes(
+            stokes=object(),
+            radius_inner=0.55,
+            radius_outer=1.0,
+            harmonic_degree=2,
+            projection="nodal",
+        )
+
+
 def test_rotated_stokes_adapter_matches_zhong_table_2():
     radius_inner = 0.55
     radius_outer = 1.0
@@ -229,6 +240,21 @@ def test_rotated_stokes_adapter_matches_zhong_table_2():
         gravity=9.8,
         gravitational_constant=6.67e-11,
     )
+    reaction_response = uw.postprocessing.geoid.spherical_shell_response_from_rotated_stokes(
+        stokes=stokes,
+        radius_inner=radius_inner,
+        radius_outer=radius_outer,
+        harmonic_degree=2,
+        internal_load_radius=rint,
+        internal_load_coefficient=1.0,
+        include_self_gravity=True,
+        surface_density_contrast=3300.0,
+        cmb_density_contrast=5400.0,
+        planet_radius=6370000.0,
+        gravity=9.8,
+        gravitational_constant=6.67e-11,
+        projection="reaction",
+    )
 
     assert np.isclose(response.surface_topography, 0.41920, rtol=0.10)
     assert np.isclose(response.cmb_topography, 0.77060, rtol=0.10)
@@ -238,3 +264,7 @@ def test_rotated_stokes_adapter_matches_zhong_table_2():
     assert np.isclose(response.self_gravity.cmb_topography, 0.93130, rtol=0.10)
     assert np.isclose(response.self_gravity.surface_geoid, 0.04486, rtol=0.10)
     assert np.isclose(response.self_gravity.cmb_geoid, 0.05461, rtol=0.10)
+    assert np.isclose(reaction_response.surface_topography, 0.41920, rtol=0.03)
+    assert np.isclose(reaction_response.cmb_topography, 0.77060, rtol=0.03)
+    assert np.isclose(reaction_response.surface_geoid, 0.02579, rtol=0.03)
+    assert np.isclose(reaction_response.cmb_geoid, 0.03206, rtol=0.03)