Add Zhong 2008 propagator-matrix analytic response

docs/api/analytic.md

@@ -17,6 +17,18 @@ whichever one you use.
    :show-inheritance:
 ```
 
+## Zhong spherical-shell response oracle
+
+`Zhong2008` is a mesh-independent propagator-matrix oracle. It returns the
+boundary response coefficients used by the Zhong et al. (2008) benchmark; it
+does not implement the symbolic mesh-field contract above.
+
+```{eval-rst}
+.. automodule:: underworld3.analytic.zhong2008
+   :members:
+   :show-inheritance:
+```
+
 ## See also
 
 - {doc}`solvers` — the solvers these solutions validate.

docs/developer/CHANGELOG.md

@@ -273,6 +273,16 @@ 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.
+- `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
+  response printed in Zhong Tables 2 and 3; geoid and self-gravity are delegated
+  to the generic postprocessing functions above.
 
 ### Generalized Geometric Multigrid via Custom Prolongation (July 2026)
 

docs/developer/subsystems/analytic-solutions.md

@@ -818,6 +818,98 @@ Two class attributes carry this:
 | `symbolic` | fields are SymPy on `mesh.X`. False means the residual gates cannot be applied at all |
 | `requires` | name of an optional package, or `None` |
 
+## Numeric response oracles
+
+Not every published analytical benchmark is a pointwise field on a mesh. The
+Zhong et al. (2008) spherical-shell benchmark publishes scalar response
+functions computed by a semi-analytic radial propagator. It belongs in
+`uw.analytic`, but it does **not** subclass `AnalyticSolution` and is not listed
+by `available()` because it cannot satisfy that class's symbolic mesh-field
+contract.
+
+```python
+reference = uw.analytic.Zhong2008(
+    harmonic_degree=2,
+    radius_inner=0.55,
+    radius_outer=1.0,
+    internal_load_radius=0.775,
+)
+response = reference.response()
+
+print(response.surface_characteristic_velocity)
+print(response.self_gravity.surface_topography)
+```
+
+### Radial state and propagator
+
+For one spherical-harmonic degree `l`, define `L = l(l + 1)` and use
+`v = log(r)`. The four-component poloidal state follows Hager and O'Connell
+(1981):
+
+```text
+u = (y1, y2, r sigma_rr / eta0, r sigma_r_perp / eta0)^T
+```
+
+`y1` is radial velocity and `y2` is the characteristic horizontal velocity.
+Within one constant-viscosity layer, with dimensionless viscosity `eta`,
+
+```text
+du/dv = A u
+
+    [ -2          L          0   0     ]
+    [ -1          1          0   1/eta ]
+A = [ 12 eta     -6 L eta    1   L     ]
+    [ -6 eta  2(2L-1) eta   -1  -2     ]
+```
+
+The exact layer transfer is `exp(A log(r_b/r_a))`. Transfers are multiplied in
+increasing-radius order. Because the state stores physical tractions relative
+to one fixed reference viscosity, all four components remain continuous at a
+viscosity interface.
+
+The outward radial delta load of coefficient `a` at `rint` gives the jump
+
+```text
+u(rint+) - u(rint-) = (0, 0, -rint a, 0)^T.
+```
+
+Impermeable free slip requires `u[0] = u[3] = 0` at the CMB and surface. The
+two remaining CMB state components are obtained from a dense two-by-two solve.
+
+### Recovering published quantities
+
+With CMB and surface states `ub` and `us`, respectively:
+
+```text
+surface topography without self-gravity = -us[2] / radius_outer
+CMB topography without self-gravity     =  ub[2] / radius_inner
+surface characteristic velocity         =  us[1]
+CMB characteristic velocity             =  ub[1]
+surface horizontal divergence            = -L us[1] / radius_outer
+CMB horizontal divergence                = -L ub[1] / radius_inner
+```
+
+The topography coefficients then pass through the generic
+`uw.postprocessing.geoid.spherical_shell_geoid_response()` and
+`spherical_shell_self_gravity_response()` functions. The propagator module does
+not carry a second geoid implementation.
+
+`Zhong2008Response.surface_topography` and `.surface_geoid` are the
+no-self-gravity quantities. The corresponding Table 2 or Table 3 quantities
+are under `.self_gravity`. Velocity is unchanged by that postprocessing.
+
+### Validation
+
+`tests/test_1029_analytic_zhong2008.py` checks every parenthesized analytical
+entry in Zhong et al. (2008) Tables 2 and 3: three load depths, four harmonic
+degrees, isoviscous and `10^4`-lid cases, and eight response quantities per
+case. The 192 comparisons agree within the precision printed in the paper. The
+test also checks both free-slip boundary states, load linearity, no-self-gravity
+recovery, and invalid input handling.
+
+This is stronger than copying table values into benchmark scripts: the values
+are independently recomputed from the governing radial system for every case.
+
 ## Adding a new solution
 
 ### The format to supply a new solution in

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/analytic/__init__.py

@@ -37,7 +37,21 @@
 from .kramer import CylindricalStokes
 from .richards import GardnerSteady, GardnerTransient
 from .transport import AdvectedFront, ErfcDiffusion, Poisson1D, TwoLayerDarcy
-from .velic import SolA, SolB, SolC, SolCx, SolDA, SolDB2d, SolDB3d, SolH, SolKx, SolKz, SolM, SolNL
+from .velic import (
+    SolA,
+    SolB,
+    SolC,
+    SolCx,
+    SolDA,
+    SolDB2d,
+    SolDB3d,
+    SolH,
+    SolKx,
+    SolKz,
+    SolM,
+    SolNL,
+)
+from .zhong2008 import Zhong2008, Zhong2008Response
 
 __all__ = [
     "AnalyticSolution",
@@ -64,6 +78,8 @@
     "SolM",
     "SolNL",
     "TwoLayerDarcy",
+    "Zhong2008",
+    "Zhong2008Response",
     "available",
     "describe",
     "is_available",