Add generic spherical-shell geoid and self-gravity postprocessing

[gthyagi]

Summary

Add generic spherical-harmonic geoid and self-gravity postprocessing for two-boundary spherical shells, with an optional internal load.

The reusable mathematics now lives in uw.postprocessing.geoid. Zhong et al. (2008) supplies the validation case and exact benchmark parameters, but is not encoded in public class or function names.

This revision deliberately reuses UW3's existing rotated-free-slip stress recovery. It does not add another topography implementation, automatic BC dispatch, CBF fallback, constrained-multiplier path, Stokes.geoid() facade, or semi-analytical Stokes solver.

API

Pure coefficient operations, independent of a Stokes discretisation:

uw.postprocessing.geoid.spherical_shell_geoid_response(...)
uw.postprocessing.geoid.spherical_shell_self_gravity_response(...)

Convenience extraction from a completed rotated-free-slip 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 calls stokes.boundary_normal_traction(boundary, mass="auto") and projects surface/CMB responses onto the unnormalised axisymmetric P_l^0 harmonic. Other harmonic orders or topography-recovery methods can call the pure coefficient functions directly.

Inputs and conventions

• Surface, CMB, and internal-load coefficients must use one consistent spherical-harmonic normalisation.
• The radial potential kernel depends on harmonic degree but not order.
• The internal load is optional; a nonzero coefficient requires a radius strictly inside the shell.
• Shell radii and response coefficients are nondimensional.
• Density contrasts are signed and use SI units.
• Model-specific density contrasts, planet radius, and gravity are required explicitly when self-gravity is enabled. There are no hidden Zhong/Earth defaults.
• The universal gravitational constant has a CODATA default and can be overridden with a paper's rounded value.

Mathematics

One radial Green-function operator and optional load vector are shared by both paths:

N = G h + n_load
(I - Q G) h_self_gravity = h + Q n_load
N_self_gravity = G h_self_gravity + n_load

Here h contains surface/CMB topography coefficients, N contains surface/CMB geoid coefficients, and Q contains the self-gravity density factors.

This is the geoid/self-gravity algebra applied after topography is known. A Zhong propagator-matrix reference solver would solve the spherical-harmonic radial Stokes equations and belongs separately in uw.analytic.

Validation

• serial Level 2 geoid suite: 12 passed;
• four-rank focused MPI regression: 1 passed on every rank;
• complete Level 1 suite: 1476 passed, 49 skipped, 2 xfailed;
• pure no-load, internal-load, self-gravity, and validation-error coverage;
• end-to-end rotated Stokes comparison against Zhong Table 2;
• production isoviscous serial 1/4 and eight-rank 1/8 runs reproduce prior coefficients to printed precision;
• production layered eight-rank 1/8 run reproduces its prior coefficients.

Scope

This PR changes only the postprocessing package, focused tests, and concise developer documentation. The independent 3D boundary-flux memory-scaling correction remains #593. Mesh-files-only generation remains separately scoped in #592.

[gthyagi]

@lmoresi The benchmark integration is now updated to this narrowed API. Could you please double-check the Appendix A sign/normalisation conventions and whether this rotated-only adapter is the right UW3 abstraction? If an existing implementation is preferable, this focused commit can be discarded without affecting the separate memory fix in #593.

Production validation from benchmark commit 99fc132:

Case	Resolution / ranks	Surface topo	CMB topo	Surface geoid	CMB geoid	Surface topo SG	CMB topo SG	Surface geoid SG	CMB geoid SG
Isoviscous	1/4, serial	0.40907	0.77403	0.02383	0.03183	0.48410	0.93222	0.04173	0.05377
Isoviscous	1/8, 8 ranks	0.40599	0.79339	0.02357	0.03377	0.48070	0.95994	0.04156	0.05661
Layered	1/8, 8 ranks	0.43175	0.42516	0.02198	-0.00518	0.49313	0.41879	0.03414	-0.00216

The isoviscous run is within 3.82% for self-gravity topography and 7.36% for self-gravity geoid at 1/8. The layered surface/geoid response remains under-resolved at 1/8, while CMB velocity and self-gravity topography errors are +0.55% and -3.81%. Stokes/post-processing timings were 10.08/1.27 s (isoviscous) and 58.67/1.26 s (layered). Focused tests pass: 8 serial geoid tests and the dedicated adapter test on 4 MPI ranks.