Add spherical-shell topography and geoid postprocessing

[gthyagi]

Summary

This PR proposes a reusable UW3 workflow for spherical-shell dynamic topography, geoid, and Zhong-style self-gravity response calculations. It also fixes the MPI memory scaling of the consistent P2 boundary-flux recovery used by rotated free-slip topography.

The implementation is intentionally split into two commits:

1. Add spherical-shell topography and geoid postprocessing
2. Avoid replicated 3D boundary-flux recovery

Proposed workflow

uw.postprocessing.spherical_shell_dynamic_response(...) computes the surface and CMB response in three stages:

1. Recover topography coefficients using the available free-slip formulation:
   • constrained multiplier when both boundaries have multipliers;
   • rotated constraint reaction when both boundaries use rotated strong free slip;
   • CBF volume-residual recovery as the fallback.
2. Project the boundary response onto the requested unnormalised P_l^0 harmonic, with the outer/CMB sign conventions handled internally.
3. Apply the Zhong et al. (2008) Appendix A two-boundary operator to obtain no-self-gravity geoid coefficients and, optionally, solve the coupled 2 x 2 self-gravity correction.

Stokes.geoid(...) is a convenience facade. It validates spherical geometry, radii, harmonic degree, and solve convergence, then delegates all mathematics to uw.postprocessing.

The intended user code is therefore:

response=stokes.geoid(
radius_inner=0.55,
radius_outer=1.0,
radius_internal=0.775,
harmonic_degree=5,
self_gravity=True,
)

Users do not need to select separate constrained, rotated, or CBF topography routines.

Parallel memory fix

The previous 3D consistent-P2 boundary recovery used MPI.allgather, after which every rank retained the global Python boundary mesh, assembled the same sparse surface-mass matrix, and called spsolve independently.

The new path:

• gathers partial reactions, trace elements, and requested coordinate keys to rank zero;
• assembles and solves the same global boundary-mass system once;
• scatters only each rank's requested recovered values;
• preallocates COO arrays and propagates root-side failures collectively.

The recovery mathematics, mean-removal gauge, and local return values are unchanged.

Zhong benchmark evidence

The production benchmark is bench_010_zhong2008_isoviscous_response.py.

For cellsize=1/64, l=5, depth 0.5, 192 MPI ranks:

Quantity	Zhong	UW3	Error
Surface velocity	0.003593	0.00359301	+0.00%
CMB velocity	0.003733	0.00373243	-0.02%
Surface topography, self-gravity	0.4238	0.421141	-0.63%
CMB topography, self-gravity	0.7235	0.730447	+0.96%
Surface geoid, self-gravity	0.02427	0.0240307	-0.99%
CMB geoid, self-gravity	0.02543	0.0257665	+1.32%

Every coefficient improves from 1/32 to 1/64. The fixed 1/64 run completed in 7:29 with 301.77 GB peak memory; the old path was killed during stokes.geoid() after exhausting 576 GB. Detailed measurements are in the high-resolution memory report.

Validation

• tests/test_1070_postprocessing_geoid.py: 5 passed.
• tests/test_1019_boundary_flux.py -k 3d: 3 passed.
• Four-rank tests/parallel/test_1065_boundary_flux_parallel.py -k 3d: 2 passed on every rank.
• Four-rank spherical rotated-topography comparison reached the platform-specific Gmsh reference-anchor skip after partition independence was asserted.
• Gadi 1/32 results match the pre-change coefficients to floating-point roundoff.
• Gadi 1/64, 192-rank production run completed with finite metrics and complete HDF5 output.

Review request

@lmoresi, could you please double-check the geoid computation workflow, especially:

• whether this functionality belongs in uw.postprocessing with a Stokes.geoid() facade;
• the automatic constrained/rotated/CBF topography-source dispatch;
• the surface/CMB harmonic projection and sign conventions;
• the Zhong Appendix A operator and self-gravity correction;
• whether the API should return response coefficients or a more general boundary field representation.

If UW3 already has, or should use, a better implementation or abstraction, please recommend that path. I am happy to discard or replace this implementation rather than merge duplicate or incorrectly placed machinery.

Scope

This PR intentionally excludes the unrelated mesh-files-only work already under review in #581 and the parallel non-finite interpolation fallback. It is limited to spherical topography/geoid postprocessing and the boundary-flux memory correction required by that workflow.

[gthyagi]

Superseded by #591 after renaming the fork branch from codex/spherical-geoid-postprocessing to feature/spherical-geoid-postprocessing. The commit and PR description are unchanged.