Add Zhong 2008 propagator-matrix analytic response

[gthyagi]

Dependency

This is a stacked PR on top of #591, which adds the generic
uw.postprocessing.geoid coefficient API. The Zhong oracle calls those pure
functions rather than carrying a second geoid/self-gravity implementation.

Until #591 merges, this comparison also shows its three commits. After #591
merges into development, the remaining diff is the two commits
95b34701 and fb9ce95b.

What this adds

uw.analytic.Zhong2008 is a mesh-independent semi-analytical reference for the
spherical-shell delta-load Stokes benchmarks in Zhong et al. (2008).

For one harmonic degree l, it:

• propagates the four-component Hager--O'Connell poloidal state in
  v = log(r) with an exact matrix exponential in each constant-viscosity
  layer;
• accepts arbitrary piecewise-constant radial viscosity profiles;
• applies the radial delta load as the traction-state jump
  (0, 0, -rint * load, 0);
• solves the two remaining CMB state coefficients from impermeable free slip
  (radial velocity = tangential traction = 0) at both shell boundaries;
• returns no-self-gravity topography/geoid, signed characteristic horizontal
  velocity, and horizontal velocity divergence at the surface and CMB;
• returns self-gravity-corrected topography/geoid through the generic API from
  Add generic spherical-shell geoid and self-gravity postprocessing #591.

This is intentionally a numeric response oracle rather than an
AnalyticSolution(mesh): the paper publishes boundary response coefficients,
not pointwise symbolic velocity/pressure fields. It is exported from
uw.analytic but is not registered by available(), whose contract requires a
mesh-constructible symbolic solution.

Radial system

With L = l(l + 1) and state

u = (y1, y2, r sigma_rr / eta0, r sigma_r_perp / eta0)^T,

the constant-layer system is

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     ]

and u(rb) = exp(A log(rb / ra)) u(ra).

API example

reference = uw.analytic.Zhong2008(
    harmonic_degree=2,
    radius_inner=0.55,
    radius_outer=1.0,
    internal_load_radius=0.775,
).response()

reference.surface_characteristic_velocity
reference.self_gravity.surface_topography

The layered Table 3 case is selected with:

viscosity_interfaces=(0.971875,)
viscosities=(1.0, 1.0e4)

Published-reference validation

The focused test checks every parenthesized propagator entry in Zhong et al.
(2008) Tables 2 and 3:

Cases	Depths	Degrees	Quantities per case	Published values
Isoviscous	0.25d, 0.5d, 0.75d	2, 5, 8, 15	8	96
10^4 viscosity lid	0.25d, 0.5d, 0.75d	2, 5, 8, 15	8	96

All 192 values agree within the precision printed in the paper. Additional
tests cover both free-slip boundary states, response linearity in load
amplitude, no-self-gravity recovery, input validation, and public API behavior.

Representative degree-2, mid-mantle responses:

Case	s	b	h_t	h_b	U_t	U_b
Isoviscous oracle	0.499849	0.931260	0.0448629	0.0546149	-0.0100641	0.0118606
Paper	0.4998	0.9313	0.04486	0.05461	-0.01006	0.01186
Layered oracle	0.610374	0.435376	0.0578926	0.00675446	-1.25811e-05	0.00555605
Paper	0.6104	0.4354	0.05789	0.006754	-1.258e-05	0.005556

Checks

• python -m pytest tests/test_1029_analytic_zhong2008.py -v:
  34 passed
• analytic contract/optional/oracle group: 66 passed, 15 skipped
• ./uw test Level 1: 1476 passed, 49 skipped, 2 xfailed
• documentation audit still exits on the repository-wide baseline of 138
  undocumented API items; it reports no new missing Zhong API item.

Review focus

Please check:

1. the four-state matrix and traction scaling against the intended
   Hager--O'Connell convention;
2. the sign and rint factor in the delta-load jump;
3. whether a mesh-independent response oracle should remain directly exported
   but outside the symbolic available() registry;
4. the constructor representation of radial viscosity interfaces.

References

• Hager & O'Connell (1981), JGR 86, 4843--4878,
  doi:10.1029/JB086iB06p04843.
• Zhong et al. (2008), G3 9, Q10017, doi:10.1029/2008GC002048.

[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 after #591 (cdcec34), which this was stacked on — the geoid work is recorded there rather than here.

Two things from the cluster review that are carried forward rather than fixed by this merge:

• tests/test_1029_ is now used twice (C2): test_1029_analytic_faulted_medium.py from Barr & Houseman faulted-medium analytic solution in uw.analytic #550 and test_1029_analytic_zhong2008.py from here. Both are present after the merge and both run; the number is what collides. Tracked for a renumber to 1030.
• Zhong2008 and Zhong2008Response are exported but not registered (finding 1). That is correct — neither subclasses AnalyticSolution, so uw.analytic.available(), describe(), the conformance suite and convention_audit.py all pass them by. It is also invisible from the namespace, where everything else carries the validation contract.

The FaultedMedium removal shown in this PR's diff is an artefact of diffing against the merge base; we resolved the merge before merging and confirmed the result carries FaultedMedium and Zhong2008 side by side.