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.

[lmoresi]

Adversarial review — spherical-shell cluster (#591, #592, #593, #594)

Reviewed as a group because they overlap: #594 contains #591, and all four land
in the spherical/postprocessing path. Findings below are labelled by the PR they
apply to; the cluster-level ones come first because a single-PR reading misses
them.

Cluster

C1. #594 is stacked on #591 — merging #594 merges #591. pr591 is an
ancestor of pr594, and postprocessing/geoid.py and
postprocessing/__init__.py are byte-identical between them (git diff pr591 pr594 -- ... is empty). #594 is #591 plus two commits (95b3470, fb9ce95).
Neither PR body says so. Merge #591 first and #594 becomes a two-commit
follow-up; merge #594 first and #591 closes as already-merged. What must not
happen is reviewing the geoid code twice as though two implementations were on
offer.

C2. tests/test_1029_ is used twice after the merge. #550 landed
tests/test_1029_analytic_faulted_medium.py yesterday; #594 adds
tests/test_1029_analytic_zhong2008.py. Both survive the merge — we resolved
the tree and they are both there. Two different subjects under one number
defeats the numbering. Zhong should move to 1030.

C3. Not a defect, recorded because the diff view says otherwise. GitHub
shows #594 REMOVING FaultedMedium from analytic/__init__.py — the import,
the __all__ entry and the _SOLUTIONS entry. That is the diff against the
merge base, not against development. We resolved the merge (git merge-tree --write-tree origin/development pr594): it is clean, and the resulting
__init__.py carries FaultedMedium and Zhong2008 side by side, with
barr_houseman_96.py intact. No rebase is needed on that account.

#594 — Zhong 2008

1. Zhong2008 and Zhong2008Response are in __all__ but not in
_SOLUTIONS. So uw.analytic.available() does not list them, describe()
does not know them, and neither the conformance suite nor
convention_audit.py measures them. That follows from what they are — neither
subclasses AnalyticSolution, and there is no sample_points, boundaries,
apply_boundary_conditions or set_fields — so the registry would reject them
anyway. It is the right call and it is invisible: everything else reachable as
uw.analytic.X is a solution carrying the validation contract, and these two are
not. Worth a sentence in the module docstring saying they are a response
calculator rather than a solution, so the next reader does not go looking for the
momentum-residual row.

#593 — boundary-flux recovery

2. The memory bound moves, it does not fall. Replacing comm.allgather with
gather-to-rank-zero removes the ×nranks replication, which is what #408 asks
for. Rank zero still holds the entire global boundary and solves the entire
sparse system alone, so the worst-case single-rank memory is unchanged and the
solve is now serial where it used to be redundantly parallel. On the mesh sizes
that motivated #408 that is a clear win; at the scale where rank zero itself is
the constraint it is not, and the docstring's "avoids replicating the global P2
surface mesh and sparse solve on every rank" is true without saying which of the
two limits now binds.

3. Every rank must reach the gather. The recovery is collective by
construction now. A rank owning no part of the boundary has nothing to
contribute, and it still has to enter gather and the scatter that follows it.
We have not run this at a rank count that starves a rank of boundary facets, and
the tests in the PR do not construct one. This is the same shape as the defect
fixed in #596, and it is worth a parallel test with a boundary that does not
reach every rank.

#592 — internal spherical shell

4. It edits test_0502_boundary_integrals.py, which is on the collection
ratchet. That file does module-level work at import and is exempted by name in
tests/conftest.py (#587). Editing it is fine; adding more module-level work to
it is not, because the exemption hides it. Worth checking the added lines land
inside test functions.

5. discretisation_mesh.py is touched by #596 as well. The regions are far
apart — #592 near line 100, #596 near 8500 — so they merge textually. No action;
recorded so the collision is not discovered at merge time.

What we did not do

None of the four was run. The findings above come from reading the diffs and
resolving the merges mechanically; C1, C2 and C3 are verified against the git
objects, and 1 to 5 are read from the source. Where a claim needed a measurement
we have said so rather than asserting it.

Underworld development team with AI support from Claude Code

[lmoresi]

Merging this ahead of #594, per C1 in the cluster review above. #594 is stacked on this branch — pr591 is an ancestor of pr594 and the geoid files are byte-identical between them — so merging #594 first would have merged this one silently and closed it as already-merged, with the geoid work recorded against the Zhong PR rather than its own.

Taking them in order keeps the two changes separable in the history.