Add generic spherical-shell geoid and self-gravity postprocessing

docs/developer/CHANGELOG.md

@@ -273,6 +273,11 @@ component exactly — correct on curved, tilted, and deformed boundaries (#293).
   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).
+- `uw.postprocessing.geoid` provides generic spherical-shell geoid and
+  self-gravity coefficient functions. Its rotated-Stokes adapter projects the
+  existing boundary traction onto an axisymmetric harmonic; the pure functions
+  also accept coefficients recovered by other methods and an optional internal
+  load.
 
 ### Generalized Geometric Multigrid via Custom Prolongation (July 2026)
 

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

@@ -74,6 +74,74 @@ Combine the two rules — project the τ components with `linear_solver()` for t
 cheap linear solve, then compose `σ_rr` analytically — for accurate boundary
 stress recovery that also scales.
 
+## 3. Spherical-shell geoid and self-gravity response
+
+Rotated free slip already exposes normal traction through
+`Stokes.boundary_normal_traction()`. The convenience adapter projects that
+recovered traction onto the unnormalised axisymmetric `P_l^0` harmonic and
+applies the spherical-shell geoid operator:
+
+```python
+stokes.solve()
+
+response = uw.postprocessing.geoid.spherical_shell_response_from_rotated_stokes(
+    stokes=stokes,
+    radius_inner=0.55,
+    radius_outer=1.0,
+    harmonic_degree=2,
+    internal_load_radius=0.775,
+    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,
+)
+```
+
+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.
+
+When surface and CMB topography coefficients are already available, call
+`uw.postprocessing.geoid.spherical_shell_geoid_response()` or
+`uw.postprocessing.geoid.spherical_shell_self_gravity_response()` directly.
+These functions are pure post-processing, work for any spherical-harmonic
+order with a consistent coefficient normalisation, and do not require a Stokes
+object. They support non-negative harmonic degrees, and the internal load is
+optional. The rotated-Stokes adapter requires degree one or greater because
+normal-traction recovery removes the degree-zero mean.
+
+The density contrasts, dimensional outer-radius scale, and gravity are required
+when self-gravity is enabled. They deliberately have no Earth- or
+benchmark-specific defaults. The universal gravitational constant defaults to
+the current CODATA value and can be overridden when reproducing a paper's
+rounded constant.
+
+The calculation is expressed as one linear operator:
+
+```text
+N = G h + n_load
+(I - Q G) h_self_gravity = h + Q n_load
+```
+
+where `h` contains surface and CMB topography, `N` contains their geoid
+responses, and `Q` contains the two self-gravity density factors. Sharing `G`
+and `n_load` between the no-self-gravity and self-gravity paths avoids separate
+scalar implementations of the same coefficients.
+
+This module does not compute a benchmark's semi-analytical Stokes solution.
+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.
+
 ## See also
 
 - Issues [#156] (projection solver settings), [#157] (projection memory),

src/underworld3/__init__.py

@@ -201,6 +201,7 @@ def view():
 import underworld3.maths
 import underworld3.swarm
 import underworld3.systems
+import underworld3.postprocessing
 import underworld3.maths
 import underworld3.utilities
 import underworld3.model

src/underworld3/postprocessing/__init__.py

@@ -0,0 +1,10 @@
+r"""
+Post-processing helpers for Underworld3 models.
+
+The post-processing package contains derived diagnostics that are useful across
+multiple model scripts but do not belong in a solver implementation.
+"""
+
+from . import geoid
+
+__all__ = ["geoid"]