Add cylindrical annulus geoid postprocessing

[gthyagi]

Summary

• add the Simons cylindrical sheet Green function for positive Fourier modes
• assemble two-boundary annulus potential, geoid, and self-gravity coefficients
• add an MPI-safe adapter for completed rotated-free-slip Stokes solves
• document sign, density, Fourier, unit, and n=0 gauge conventions

The coefficient API is independent of the physical-space Poisson meshes and convergence drivers.

Kramer-Simons annulus validation

The fixed nondimensional case uses Fourier mode n=2, radii ri=1.22, rload=2.0, ro=2.22, G=1, unit boundary density contrasts, unit internal sheet coefficient, and unit reference gravities. Consequently geoid N and potential Phi have the same numerical coefficient. Errors below use abs(numerical-reference)/abs(reference).

Spectral response with varying Stokes h

Here h is the target annulus Stokes element size. Each row recovers dynamic topography from that Stokes solve and applies the cylindrical spectral operator. The reference is the fully analytical Kramer traction plus Simons gravity response.

h	Mantle cells	Outer N spectral	Outer N analytic	Outer rel. error	Inner N spectral	Inner N analytic	Inner rel. error
1/4	476	+1.840897156	-0.663601170	377.410173%	+0.729297687	-0.502494144	245.135559%
1/8	1,684	-0.396340543	-0.663601170	40.274285%	-0.373554174	-0.502494144	25.659995%
1/16	6,578	-0.645005977	-0.663601170	2.802164%	-0.494736091	-0.502494144	1.543909%
1/32	26,376	-0.662116990	-0.663601170	0.223656%	-0.502457968	-0.502494144	0.007199%

The coarse 1/4 and 1/8 errors are dominated by recovered Stokes traction/topography. Because boundary topography and the internal load substantially cancel, the under-resolved 1/4 topography also changes the total-geoid sign. The converged 1/32 result has 0.224% outer and 0.0072% inner error.

Physical Poisson response with varying gravity h

This comparison isolates the Poisson discretization: every row uses the same numerical topography coefficients from the 1/32 Stokes solution, while only the P2 gravity mesh spacing h changes. The scalable gravity domain has one radial buffer layer on each side of the mantle: rmin = ri - h and rmax = ro + h. Exact regular/decaying single-mode Robin continuation is imposed at rmin and rmax. The spectral columns use exactly the same topography and source coefficients as the corresponding Poisson solve.

h	Gravity cells	Potential DOFs	Radial extent	Outer Phi Poisson	Outer Phi spectral	Outer rel. error	Inner Phi Poisson	Inner Phi spectral	Inner rel. error
1/4	662	1,414	0.97000-2.47000	-0.660430832	-0.662116990	0.254662%	-0.491171097	-0.502457968	2.246331%
1/8	2,040	4,254	1.09500-2.34500	-0.661847328	-0.662116990	0.040727%	-0.499466553	-0.502457968	0.595356%
1/16	7,282	14,912	1.15750-2.28250	-0.662076818	-0.662116990	0.006067%	-0.501710188	-0.502457968	0.148824%
1/32	27,772	56,238	1.18875-2.25125	-0.662107296	-0.662116990	0.001464%	-0.502269626	-0.502457968	0.037484%

The inner-interface error is approximately second order under each halving of h. At h=1/32 the one-layer domain uses 27,772 cells, only 18.73% of the 148,288-cell full-disc reference, while giving smaller errors at both interfaces.

Poisson buffer-width and exterior-radius study

At fixed h=1/32, the gravity-domain padding was varied independently. The first three rows extend the physical annulus by h, 2h, or 3h on both radial sides; the last row is the legacy full disc from the origin to exterior radius 4.44. Every row uses the same P2 discretization, 1/32 Stokes topography, and spectral targets: outer Phi=-0.662116990 and inner Phi=-0.502457968.

Domain extension	Inner truncation	Exterior radius	Gravity cells	Outer Phi Poisson	Outer rel. error	Inner Phi Poisson	Inner rel. error
h	1.18875	2.25125	27,772	-0.662107296	0.001464%	-0.502269626	0.037484%
2h	1.15750	2.28250	29,740	-0.662107026	0.001505%	-0.502269395	0.037530%
3h	1.12625	2.31375	31,856	-0.662106778	0.001542%	-0.502268883	0.037632%
Full disc	0.00000	4.44000	148,288	-0.662102946	0.002121%	-0.502262510	0.038901%

One h of padding is sufficient: increasing to 2h or 3h changes the interface coefficients only slightly, while the h mesh is also more accurate than the full 0-4.44 disc. At h=1/32 it uses 18.73% of the full-disc cells.

Validation

• 89 relevant serial Kramer, Zhong, spherical-geoid, and cylindrical-geoid tests passed
• 20 tests passed per rank in a two-rank MPI run
• all 12 outputs matched the frozen Kramer-Simons reference at the stored precision
• Kramer Stokes, Poisson geoid, and physical self-gravity passed their standard serial 1/4 and eight-rank 1/8 runs
• the 1/32 extension study and all three retained convergence checkers passed
• Black, mypy, deprecated-pattern, and diff checks passed
• Sphinx documentation build succeeded
• repository Level-1 run reached 1636 passed, 35 skipped, and 2 xfailed before an unrelated timeout-reliability test left its divergent subprocess running

Scope and limitations

• Fourier modes n >= 1 are supported
• n = 0 is rejected explicitly because its logarithmic solution requires a potential gauge
• tests begin at tier C while the API receives review

[lmoresi]

Reviewed. The physics and the validation are solid; my concern is architectural, and it is cheap to address now and expensive later.

What is good

The Kramer–Simons convergence campaign is real evidence, not a smoke test: 0.224% outer and 0.0072% inner at h=1/32, with clean second-order behaviour on the Poisson side under each halving. The input validation in cylindrical_cosine_boundary_coefficient is careful — shape, finiteness, minimum sample count all checked with specific messages. The buffer-width study justifying the one-layer gravity domain (27,772 cells against a 148,288-cell full-disc reference, with smaller errors at both interfaces) is a genuinely useful result on its own.

The concern: this adds a second bespoke projection path

cylindrical_cosine_boundary_coefficient gathers boundary samples to rank zero, sorts them by arctan2, and integrates with a hand-rolled trapezoid. That is the 2-D analogue of _project_spherical_harmonic_samples, which #647 argues should be retired: gather the boundary, re-derive an ordering the mesh already knows, integrate with bespoke quadrature instead of the FE boundary integral.

Three consequences worth weighing:

1. It assumes the boundary is a circle. Sorting by θ and trapezoid-integrating is exact only for samples on a circle. A deformed annulus — a free surface, the case this machinery ultimately serves — breaks it, in the same way convexity limits the spherical path.
2. Chord versus arc. Integrating in θ over a faceted boundary is the Curved-boundary P2 VERTEX flux/traction values converge only slowly — document and steer consumers to midpoint/fitted quantities #414/CBF sigma_nn recovery on deformed/faceted boundaries: ~2.4% amplitude deficit (chord-vs-arc boundary mass suspected) #431 error family again, applied to an FE reaction.
3. It gathers to rank zero (MPI.COMM_WORLD.gather), which is the scaling objection in boundary_flux 3D assembly gathers the full global boundary onto every rank (memory scaling) #408, and uses COMM_WORLD rather than the mesh communicator.

The alternative is already in your own #646

boundary_normal_traction_integral contracts the assembled reaction with the test function over owned DOFs plus MPI reductions — no gather, no ordering, no pointwise recovery. For a coefficient that is strictly the better construction, and it is the one I argued in #647 should replace the spherical path. uw.maths.BdIntegral is the other UW-native option where a nodal field genuinely exists.

Routing this projection through one of those is a change to one function, not a rewrite. Doing it before merge avoids maintaining two bespoke paths and then migrating both.

Practical collisions if both land

• Add distributed harmonic projection of rotated reactions #646 and this PR both modify src/underworld3/postprocessing/geoid.py — expect a conflict.
• Both add a test_1071_* file (tests/parallel/…spherical_shell_geoid_parallel vs tests/…postprocessing_cylindrical_geoid). Different directories, same number; worth renumbering one.

Verdict

Not blocking on correctness — CI is green and the validation stands. But I would rather see the projection go through #646's integral (or BdIntegral) than merge a second gather-and-trapezoid path that #647 then has to unwind twice.

Underworld development team with AI support from Claude Code