From 8a3099b418cff8256aed69d6b74b20b1b80bc69f Mon Sep 17 00:00:00 2001
From: Tyagi <gthyagi@gmail.com>
Date: Wed, 26 Aug 2026 08:43:13 +0530
Subject: [PATCH 1/3] Add cylindrical annulus geoid response API

Implement the Simons cylindrical sheet Green's function for positive Fourier modes and assemble two-boundary potential, geoid, and self-gravity responses.

Add a rotated-free-slip Stokes adapter that projects wall reactions on rank zero and broadcasts deterministic coefficients. Cover kernel branches, Poisson jumps, mode validation, superposition, feedback solves, public exports, and partitioned MPI sampling.

The coefficient API remains independent of the mantle-convection benchmark drivers and physical-space Poisson meshes.
---
 src/underworld3/postprocessing/geoid.py       | 948 +++++++++++++++++-
 ...t_1071_postprocessing_cylindrical_geoid.py | 339 +++++++
 2 files changed, 1277 insertions(+), 10 deletions(-)
 create mode 100644 tests/test_1071_postprocessing_cylindrical_geoid.py

diff --git a/src/underworld3/postprocessing/geoid.py b/src/underworld3/postprocessing/geoid.py
index 537a79fee..c3bb6c0b9 100644
--- a/src/underworld3/postprocessing/geoid.py
+++ b/src/underworld3/postprocessing/geoid.py
@@ -1,16 +1,21 @@
-r"""Spherical-harmonic geoid and self-gravity response functions.
+r"""Spherical and cylindrical geoid and self-gravity response functions.
 
 The pure coefficient functions in this module are independent of a particular
-Stokes discretisation.  They combine surface, CMB, and optional internal-load
-coefficients through the radial Green's function for one spherical-harmonic
-degree.  A separate convenience adapter obtains the two topography coefficients
-from a completed rotated-free-slip Stokes solve.
+Stokes discretisation. They combine boundary and optional internal-load
+coefficients through the appropriate radial Green's function. Separate
+convenience adapters obtain topography coefficients from completed
+rotated-free-slip Stokes solves.
+
+The cylindrical sheet kernel follows Simons (1996), Appendix B, with its
+normalisation fixed directly by potential continuity and the radial-derivative
+jump condition.
 """
 
 from __future__ import annotations
 
 from dataclasses import dataclass
 from numbers import Integral
+from typing import Any
 
 from mpi4py import MPI
 import numpy as np
@@ -20,9 +25,19 @@
     "GeoidResponse",
     "SelfGravityResponse",
     "SphericalShellResponse",
+    "CylindricalGravityResponse",
+    "CylindricalSelfGravityResponse",
+    "CylindricalAnnulusResponse",
     "spherical_shell_geoid_response",
     "spherical_shell_self_gravity_response",
     "spherical_shell_response_from_rotated_stokes",
+    "cylindrical_sheet_potential_coefficient",
+    "cylindrical_sheet_radial_derivative_coefficient",
+    "cylindrical_annulus_potential_operator",
+    "cylindrical_annulus_geoid_response",
+    "cylindrical_annulus_self_gravity_response",
+    "cylindrical_cosine_boundary_coefficient",
+    "cylindrical_annulus_response_from_rotated_stokes",
 ]
 
 
@@ -57,6 +72,48 @@ class SphericalShellResponse:
     self_gravity: SelfGravityResponse | None = None
 
 
+@dataclass(frozen=True)
+class CylindricalGravityResponse:
+    """Potential and geoid coefficients at the outer and inner boundaries."""
+
+    outer_potential: float
+    inner_potential: float
+    outer_geoid: float
+    inner_geoid: float
+
+
+@dataclass(frozen=True)
+class CylindricalSelfGravityResponse:
+    """Self-gravity-corrected cylindrical response coefficients."""
+
+    outer_topography: float
+    inner_topography: float
+    outer_potential: float
+    inner_potential: float
+    outer_geoid: float
+    inner_geoid: float
+    q_outer: float
+    q_inner: float
+    matrix_residual_norm: float
+
+
+@dataclass(frozen=True)
+class CylindricalAnnulusResponse:
+    """Rotated-Stokes topography and cylindrical gravity coefficients."""
+
+    outer_reaction: float
+    inner_reaction: float
+    outer_reaction_mean: float
+    inner_reaction_mean: float
+    outer_topography: float
+    inner_topography: float
+    outer_potential: float
+    inner_potential: float
+    outer_geoid: float
+    inner_geoid: float
+    self_gravity: CylindricalSelfGravityResponse | None = None
+
+
 def _validate_geometry(
     radius_inner: float,
     radius_outer: float,
@@ -68,7 +125,9 @@ def _validate_geometry(
     except (TypeError, ValueError) as error:
         raise TypeError("The shell radii must be real numbers.") from error
     if not np.all(np.isfinite((ri, ro))) or not 0.0 < ri < ro:
-        raise ValueError("Expected finite radii ordered as 0 < radius_inner < radius_outer.")
+        raise ValueError(
+            "Expected finite radii ordered as 0 < radius_inner < radius_outer."
+        )
     if isinstance(harmonic_degree, bool) or not isinstance(harmonic_degree, Integral):
         raise TypeError("harmonic_degree must be an integer.")
     degree = int(harmonic_degree)
@@ -100,9 +159,13 @@ def _spherical_shell_geoid_operator(
         try:
             rint = float(internal_load_radius)
         except (TypeError, ValueError) as error:
-            raise TypeError("internal_load_radius must be a real number or None.") from error
+            raise TypeError(
+                "internal_load_radius must be a real number or None."
+            ) from error
         if not np.isfinite(rint) or not ri < rint < ro:
-            raise ValueError("internal_load_radius must lie strictly between the shell radii.")
+            raise ValueError(
+                "internal_load_radius must lie strictly between the shell radii."
+            )
     elif internal_load_coefficient != 0.0:
         raise ValueError(
             "internal_load_radius is required when internal_load_coefficient is nonzero."
@@ -203,8 +266,12 @@ def spherical_shell_self_gravity_response(
         [surface_topography_coefficient, cmb_topography_coefficient],
         dtype=float,
     )
-    density_contrasts = np.array([surface_density_contrast, cmb_density_contrast], dtype=float)
-    physical_constants = np.array([planet_radius, gravity, gravitational_constant], dtype=float)
+    density_contrasts = np.array(
+        [surface_density_contrast, cmb_density_contrast], dtype=float
+    )
+    physical_constants = np.array(
+        [planet_radius, gravity, gravitational_constant], dtype=float
+    )
     if not np.all(np.isfinite(topography)):
         raise ValueError("The topography coefficients must be finite.")
     if not np.all(np.isfinite(density_contrasts)):
@@ -231,6 +298,556 @@ def spherical_shell_self_gravity_response(
     )
 
 
+def _finite_float(value: Any, name: str) -> float:
+    try:
+        result = float(value)
+    except (TypeError, ValueError) as error:
+        raise TypeError(f"{name} must be a real number.") from error
+    if not np.isfinite(result):
+        raise ValueError(f"{name} must be finite.")
+    return result
+
+
+def _positive_float(value: Any, name: str) -> float:
+    result = _finite_float(value, name)
+    if result <= 0.0:
+        raise ValueError(f"{name} must be positive.")
+    return result
+
+
+def _validate_cylindrical_mode(wavenumber: int) -> int:
+    if isinstance(wavenumber, bool) or not isinstance(wavenumber, Integral):
+        raise TypeError("wavenumber must be an integer.")
+    mode = int(wavenumber)
+    if mode < 0:
+        raise ValueError("wavenumber must be non-negative.")
+    if mode == 0:
+        raise ValueError(
+            "The n=0 cylindrical mode has a logarithmic radial solution and "
+            "requires a potential gauge."
+        )
+    return mode
+
+
+def _validate_cylindrical_annulus(
+    radius_inner: float,
+    radius_outer: float,
+    wavenumber: int,
+) -> tuple[float, float, int]:
+    radius_inner = _positive_float(radius_inner, "radius_inner")
+    radius_outer = _positive_float(radius_outer, "radius_outer")
+    if radius_inner >= radius_outer:
+        raise ValueError("Expected radius_inner < radius_outer.")
+    return (
+        radius_inner,
+        radius_outer,
+        _validate_cylindrical_mode(wavenumber),
+    )
+
+
+def cylindrical_sheet_potential_coefficient(
+    *,
+    source_radius: float,
+    target_radius: float,
+    wavenumber: int,
+    surface_density_coefficient: float,
+    gravitational_constant: float = 1.0,
+) -> float:
+    r"""Return one cylindrical mass sheet's potential coefficient.
+
+    The result multiplies the unnormalised real Fourier basis
+    :math:`\cos(n\theta)`. For :math:`n\geq 1`, a sheet at radius :math:`r_s`
+    has
+
+    .. math::
+
+       \Phi_n(r_s) = \frac{2\pi G r_s\sigma_n}{n},
+
+    with radial factors :math:`(r/r_s)^n` inside and :math:`(r_s/r)^n`
+    outside. Positive density gives positive potential under the convention
+    :math:`\nabla^2\Phi=-4\pi G\rho` and :math:`\mathbf{g}=\nabla\Phi`.
+
+    Radii, density, and ``gravitational_constant`` may be dimensional or
+    nondimensional, but they must use one mutually consistent unit system.
+    The axisymmetric ``n=0`` mode is intentionally excluded because its
+    logarithmic exterior branch requires a separate potential gauge.
+
+    Parameters
+    ----------
+    source_radius, target_radius : real
+        Positive source-sheet and evaluation radii.
+    wavenumber : int
+        Positive azimuthal Fourier wavenumber.
+    surface_density_coefficient : real
+        Sheet-density coefficient multiplying :math:`\cos(n\theta)`.
+    gravitational_constant : real, default=1
+        Positive gravitational constant in the selected unit system.
+
+    Returns
+    -------
+    float
+        Potential coefficient at ``target_radius``.
+    """
+
+    source_radius = _positive_float(source_radius, "source_radius")
+    target_radius = _positive_float(target_radius, "target_radius")
+    mode = _validate_cylindrical_mode(wavenumber)
+    density = _finite_float(
+        surface_density_coefficient,
+        "surface_density_coefficient",
+    )
+    gravity_constant = _positive_float(
+        gravitational_constant,
+        "gravitational_constant",
+    )
+
+    radial_factor = (
+        min(source_radius, target_radius) / max(source_radius, target_radius)
+    ) ** mode
+    source_amplitude = 2.0 * np.pi * gravity_constant * source_radius * density / mode
+    return float(source_amplitude * radial_factor)
+
+
+def cylindrical_sheet_radial_derivative_coefficient(
+    *,
+    source_radius: float,
+    target_radius: float,
+    wavenumber: int,
+    surface_density_coefficient: float,
+    gravitational_constant: float = 1.0,
+    source_side: str | None = None,
+) -> float:
+    r"""Return :math:`d\Phi_n/dr` on either side of a cylindrical sheet.
+
+    At ``target_radius == source_radius``, ``source_side`` must be
+    ``"inside"`` or ``"outside"`` because the derivative is discontinuous.
+    The returned branches satisfy
+    :math:`[d\Phi_n/dr]_{outside-inside}=-4\pi G\sigma_n`.
+
+    Parameters
+    ----------
+    source_radius, target_radius : real
+        Positive source-sheet and evaluation radii.
+    wavenumber : int
+        Positive azimuthal Fourier wavenumber.
+    surface_density_coefficient : real
+        Sheet-density coefficient multiplying :math:`\cos(n\theta)`.
+    gravitational_constant : real, default=1
+        Positive gravitational constant in the selected unit system.
+    source_side : {"inside", "outside"}, optional
+        Radial branch when evaluating exactly on the sheet.
+
+    Returns
+    -------
+    float
+        Radial derivative coefficient at ``target_radius``.
+    """
+
+    source_radius = _positive_float(source_radius, "source_radius")
+    target_radius = _positive_float(target_radius, "target_radius")
+    mode = _validate_cylindrical_mode(wavenumber)
+    potential = cylindrical_sheet_potential_coefficient(
+        source_radius=source_radius,
+        target_radius=target_radius,
+        wavenumber=mode,
+        surface_density_coefficient=surface_density_coefficient,
+        gravitational_constant=gravitational_constant,
+    )
+
+    if target_radius == source_radius:
+        if source_side not in ("inside", "outside"):
+            raise ValueError(
+                "source_side must be 'inside' or 'outside' at the sheet radius."
+            )
+        branch_sign = 1.0 if source_side == "inside" else -1.0
+    else:
+        if source_side is not None:
+            raise ValueError(
+                "source_side is only valid when target_radius equals source_radius."
+            )
+        branch_sign = 1.0 if target_radius < source_radius else -1.0
+    return float(branch_sign * mode * potential / target_radius)
+
+
+def cylindrical_annulus_potential_operator(
+    *,
+    radius_inner: float,
+    radius_outer: float,
+    wavenumber: int,
+    outer_density_contrast: float,
+    inner_density_contrast: float,
+    gravitational_constant: float = 1.0,
+) -> np.ndarray:
+    r"""Return the two-boundary operator :math:`\Phi=G_n h`.
+
+    Rows are target boundaries ``[outer, inner]`` and columns are topographic
+    sheet sources ``[outer, inner]``. Density contrasts are signed as density
+    on the smaller-radius side minus density on the larger-radius side. Thus
+    positive outward topography creates sheet density
+    :math:`\Delta\rho\,h`.
+
+    Parameters
+    ----------
+    radius_inner, radius_outer : real
+        Positive annulus radii ordered from inner to outer.
+    wavenumber : int
+        Positive azimuthal Fourier wavenumber.
+    outer_density_contrast, inner_density_contrast : real
+        Signed density contrasts at the two boundaries.
+    gravitational_constant : real, default=1
+        Positive gravitational constant in the selected unit system.
+
+    Returns
+    -------
+    numpy.ndarray
+        Two-by-two potential operator with targets in rows and sources in
+        columns, both ordered ``[outer, inner]``.
+    """
+
+    radius_inner, radius_outer, mode = _validate_cylindrical_annulus(
+        radius_inner,
+        radius_outer,
+        wavenumber,
+    )
+    density_contrasts = np.array(
+        [
+            _finite_float(
+                outer_density_contrast,
+                "outer_density_contrast",
+            ),
+            _finite_float(
+                inner_density_contrast,
+                "inner_density_contrast",
+            ),
+        ]
+    )
+    gravity_constant = _positive_float(
+        gravitational_constant,
+        "gravitational_constant",
+    )
+    source_radii = (radius_outer, radius_inner)
+    target_radii = (radius_outer, radius_inner)
+
+    operator = np.empty((2, 2), dtype=float)
+    for row, target_radius in enumerate(target_radii):
+        for column, (source_radius, density) in enumerate(
+            zip(source_radii, density_contrasts)
+        ):
+            operator[row, column] = cylindrical_sheet_potential_coefficient(
+                source_radius=source_radius,
+                target_radius=target_radius,
+                wavenumber=mode,
+                surface_density_coefficient=density,
+                gravitational_constant=gravity_constant,
+            )
+    return operator
+
+
+def _cylindrical_internal_load_vector(
+    *,
+    radius_inner: float,
+    radius_outer: float,
+    wavenumber: int,
+    internal_load_radius: float | None,
+    internal_surface_density_coefficient: float,
+    gravitational_constant: float,
+) -> np.ndarray:
+    """Return boundary potential coefficients from an internal mass sheet."""
+
+    load_density = _finite_float(
+        internal_surface_density_coefficient,
+        "internal_surface_density_coefficient",
+    )
+    load = np.zeros(2, dtype=float)
+    if internal_load_radius is None:
+        if load_density != 0.0:
+            raise ValueError(
+                "internal_load_radius is required for a nonzero internal load."
+            )
+        return load
+
+    internal_load_radius = _positive_float(
+        internal_load_radius,
+        "internal_load_radius",
+    )
+    if not radius_inner < internal_load_radius < radius_outer:
+        raise ValueError("internal_load_radius must lie strictly inside the annulus.")
+    for index, target_radius in enumerate((radius_outer, radius_inner)):
+        load[index] = cylindrical_sheet_potential_coefficient(
+            source_radius=internal_load_radius,
+            target_radius=target_radius,
+            wavenumber=wavenumber,
+            surface_density_coefficient=load_density,
+            gravitational_constant=gravitational_constant,
+        )
+    return load
+
+
+def cylindrical_annulus_geoid_response(
+    *,
+    radius_inner: float,
+    radius_outer: float,
+    wavenumber: int,
+    outer_topography_coefficient: float,
+    inner_topography_coefficient: float,
+    outer_density_contrast: float,
+    inner_density_contrast: float,
+    outer_reference_gravity: float,
+    inner_reference_gravity: float,
+    internal_load_radius: float | None = None,
+    internal_surface_density_coefficient: float = 0.0,
+    gravitational_constant: float = 1.0,
+) -> CylindricalGravityResponse:
+    r"""Assemble annulus potential and geoid coefficients for one mode.
+
+    Potential and topography retain their physical signs. Geoid is defined as
+    :math:`N_n=\Phi_n/g_{reference}` independently at the outer and inner
+    boundaries. The optional internal source is a cylindrical sheet-density
+    coefficient in the same Fourier normalisation and unit system.
+
+    Parameters
+    ----------
+    radius_inner, radius_outer : real
+        Positive annulus radii ordered from inner to outer.
+    wavenumber : int
+        Positive azimuthal Fourier wavenumber.
+    outer_topography_coefficient, inner_topography_coefficient : real
+        Signed boundary topography coefficients.
+    outer_density_contrast, inner_density_contrast : real
+        Signed density contrasts at the two boundaries.
+    outer_reference_gravity, inner_reference_gravity : real
+        Positive gravity magnitudes used to convert potential to geoid.
+    internal_load_radius : real, optional
+        Radius of an internal sheet, strictly inside the annulus.
+    internal_surface_density_coefficient : real, default=0
+        Density coefficient of the optional internal sheet.
+    gravitational_constant : real, default=1
+        Positive gravitational constant in the selected unit system.
+
+    Returns
+    -------
+    CylindricalGravityResponse
+        Outer and inner potential and geoid coefficients.
+    """
+
+    radius_inner, radius_outer, mode = _validate_cylindrical_annulus(
+        radius_inner,
+        radius_outer,
+        wavenumber,
+    )
+    topography = np.array(
+        [
+            _finite_float(
+                outer_topography_coefficient,
+                "outer_topography_coefficient",
+            ),
+            _finite_float(
+                inner_topography_coefficient,
+                "inner_topography_coefficient",
+            ),
+        ],
+        dtype=float,
+    )
+    reference_gravity = np.array(
+        [
+            _positive_float(
+                outer_reference_gravity,
+                "outer_reference_gravity",
+            ),
+            _positive_float(
+                inner_reference_gravity,
+                "inner_reference_gravity",
+            ),
+        ],
+        dtype=float,
+    )
+    gravity_constant = _positive_float(
+        gravitational_constant,
+        "gravitational_constant",
+    )
+    operator = cylindrical_annulus_potential_operator(
+        radius_inner=radius_inner,
+        radius_outer=radius_outer,
+        wavenumber=mode,
+        outer_density_contrast=outer_density_contrast,
+        inner_density_contrast=inner_density_contrast,
+        gravitational_constant=gravity_constant,
+    )
+    load = _cylindrical_internal_load_vector(
+        radius_inner=radius_inner,
+        radius_outer=radius_outer,
+        wavenumber=mode,
+        internal_load_radius=internal_load_radius,
+        internal_surface_density_coefficient=(internal_surface_density_coefficient),
+        gravitational_constant=gravity_constant,
+    )
+
+    potential = operator @ topography + load
+    geoid = potential / reference_gravity
+    return CylindricalGravityResponse(
+        outer_potential=float(potential[0]),
+        inner_potential=float(potential[1]),
+        outer_geoid=float(geoid[0]),
+        inner_geoid=float(geoid[1]),
+    )
+
+
+def cylindrical_annulus_self_gravity_response(
+    *,
+    radius_inner: float,
+    radius_outer: float,
+    wavenumber: int,
+    outer_topography_coefficient: float,
+    inner_topography_coefficient: float,
+    outer_density_contrast: float,
+    inner_density_contrast: float,
+    outer_reference_gravity: float,
+    inner_reference_gravity: float,
+    internal_load_radius: float | None = None,
+    internal_surface_density_coefficient: float = 0.0,
+    gravitational_constant: float = 1.0,
+    outer_feedback_factor: float | None = None,
+    inner_feedback_factor: float | None = None,
+) -> CylindricalSelfGravityResponse:
+    r"""Return the two-boundary cylindrical self-gravity correction.
+
+    Holding the hydrodynamic traction fixed gives
+    :math:`h_{sg}=h+Q\Phi_{sg}`. For positive reference-gravity magnitudes the
+    default factors are :math:`Q=diag(1/g_o,1/g_i)`. Explicit factors may be
+    supplied to test a signed convention or disable either feedback row. The
+    solved equation is
+
+    .. math::
+
+       (I-QG_n)h_{sg}=h+Q\phi_{load}.
+
+    Parameters
+    ----------
+    radius_inner, radius_outer : real
+        Positive annulus radii ordered from inner to outer.
+    wavenumber : int
+        Positive azimuthal Fourier wavenumber.
+    outer_topography_coefficient, inner_topography_coefficient : real
+        Hydrodynamic topography coefficients before self-gravity feedback.
+    outer_density_contrast, inner_density_contrast : real
+        Signed density contrasts at the two boundaries.
+    outer_reference_gravity, inner_reference_gravity : real
+        Positive reference-gravity magnitudes.
+    internal_load_radius : real, optional
+        Radius of an internal sheet, strictly inside the annulus.
+    internal_surface_density_coefficient : real, default=0
+        Density coefficient of the optional internal sheet.
+    gravitational_constant : real, default=1
+        Positive gravitational constant in the selected unit system.
+    outer_feedback_factor, inner_feedback_factor : real, optional
+        Explicit diagonal entries of :math:`Q`; defaults are reciprocal
+        reference-gravity magnitudes.
+
+    Returns
+    -------
+    CylindricalSelfGravityResponse
+        Corrected topography, potential, geoid, feedback, and residual values.
+    """
+
+    radius_inner, radius_outer, mode = _validate_cylindrical_annulus(
+        radius_inner,
+        radius_outer,
+        wavenumber,
+    )
+    reference_gravity = np.array(
+        [
+            _positive_float(
+                outer_reference_gravity,
+                "outer_reference_gravity",
+            ),
+            _positive_float(
+                inner_reference_gravity,
+                "inner_reference_gravity",
+            ),
+        ],
+        dtype=float,
+    )
+    topography = np.array(
+        [
+            _finite_float(
+                outer_topography_coefficient,
+                "outer_topography_coefficient",
+            ),
+            _finite_float(
+                inner_topography_coefficient,
+                "inner_topography_coefficient",
+            ),
+        ],
+        dtype=float,
+    )
+    gravity_constant = _positive_float(
+        gravitational_constant,
+        "gravitational_constant",
+    )
+    operator = cylindrical_annulus_potential_operator(
+        radius_inner=radius_inner,
+        radius_outer=radius_outer,
+        wavenumber=mode,
+        outer_density_contrast=outer_density_contrast,
+        inner_density_contrast=inner_density_contrast,
+        gravitational_constant=gravity_constant,
+    )
+    load = _cylindrical_internal_load_vector(
+        radius_inner=radius_inner,
+        radius_outer=radius_outer,
+        wavenumber=mode,
+        internal_load_radius=internal_load_radius,
+        internal_surface_density_coefficient=(internal_surface_density_coefficient),
+        gravitational_constant=gravity_constant,
+    )
+    feedback = np.array(
+        [
+            (
+                1.0 / reference_gravity[0]
+                if outer_feedback_factor is None
+                else _finite_float(
+                    outer_feedback_factor,
+                    "outer_feedback_factor",
+                )
+            ),
+            (
+                1.0 / reference_gravity[1]
+                if inner_feedback_factor is None
+                else _finite_float(
+                    inner_feedback_factor,
+                    "inner_feedback_factor",
+                )
+            ),
+        ],
+        dtype=float,
+    )
+    q_matrix = np.diag(feedback)
+    system_matrix = np.eye(2) - q_matrix @ operator
+    right_hand_side = topography + q_matrix @ load
+    try:
+        corrected_topography = np.linalg.solve(
+            system_matrix,
+            right_hand_side,
+        )
+    except np.linalg.LinAlgError as error:
+        raise ValueError("The self-gravity feedback matrix is singular.") from error
+    corrected_potential = operator @ corrected_topography + load
+    corrected_geoid = corrected_potential / reference_gravity
+    residual = system_matrix @ corrected_topography - right_hand_side
+
+    return CylindricalSelfGravityResponse(
+        outer_topography=float(corrected_topography[0]),
+        inner_topography=float(corrected_topography[1]),
+        outer_potential=float(corrected_potential[0]),
+        inner_potential=float(corrected_potential[1]),
+        outer_geoid=float(corrected_geoid[0]),
+        inner_geoid=float(corrected_geoid[1]),
+        q_outer=float(feedback[0]),
+        q_inner=float(feedback[1]),
+        matrix_residual_norm=float(np.linalg.norm(residual)),
+    )
+
+
 def _spherical_triangle_area(a, b, c, radius: float) -> float:
     determinant = abs(float(np.dot(a, np.cross(b, c))))
     denominator = float(1.0 + np.dot(a, b) + np.dot(b, c) + np.dot(c, a))
@@ -293,6 +910,8 @@ def _rotated_topography_coefficient(
     root_error = None
     if MPI.COMM_WORLD.rank == 0:
         try:
+            if gathered_rows is None:
+                raise RuntimeError("MPI gather returned no root payload.")
             nonempty_rows = [rows for rows in gathered_rows if rows.size]
             if not nonempty_rows:
                 raise RuntimeError(f"No samples found on boundary {boundary!r}.")
@@ -406,6 +1025,10 @@ def spherical_shell_response_from_rotated_stokes(
     )
     self_gravity = None
     if include_self_gravity:
+        assert surface_density_contrast is not None
+        assert cmb_density_contrast is not None
+        assert planet_radius is not None
+        assert gravity is not None
         self_gravity = spherical_shell_self_gravity_response(
             radius_inner=ri,
             radius_outer=ro,
@@ -428,3 +1051,308 @@ def spherical_shell_response_from_rotated_stokes(
         cmb_geoid=geoid.cmb_geoid,
         self_gravity=self_gravity,
     )
+
+
+def _trapezoidal_integral(values: np.ndarray, coordinates: np.ndarray) -> float:
+    """Integrate samples without requiring NumPy's version-specific helpers."""
+
+    widths = np.diff(coordinates)
+    averages = 0.5 * (values[:-1] + values[1:])
+    return float(np.sum(widths * averages))
+
+
+def cylindrical_cosine_boundary_coefficient(
+    coords,
+    values,
+    wavenumber: int,
+) -> tuple[float, float]:
+    r"""Project circular samples onto :math:`\cos(n\theta)` and the mean.
+
+    ``coords`` must contain at least three two-dimensional Cartesian boundary
+    points. The samples may begin at any angle and need not include a duplicate
+    endpoint; this function sorts and closes the periodic interval.
+
+    Parameters
+    ----------
+    coords : array-like, shape (n, 2)
+        Cartesian circular-boundary coordinates.
+    values : array-like, shape (n,)
+        Scalar values at ``coords``.
+    wavenumber : int
+        Positive azimuthal Fourier wavenumber.
+
+    Returns
+    -------
+    coefficient, mean : tuple of float
+        Cosine-mode coefficient and degree-zero mean.
+    """
+
+    mode = _validate_cylindrical_mode(wavenumber)
+    coords = np.asarray(coords, dtype=float)
+    values = np.asarray(values, dtype=float).reshape(-1)
+    if coords.ndim != 2 or coords.shape[1] != 2:
+        raise ValueError("coords must have shape (n, 2).")
+    if coords.shape[0] != values.size:
+        raise ValueError("coords and values must contain the same number of samples.")
+    if coords.shape[0] < 3:
+        raise ValueError("At least three circular-boundary samples are required.")
+    if not np.all(np.isfinite(coords)) or not np.all(np.isfinite(values)):
+        raise ValueError("Boundary coordinates and values must be finite.")
+
+    theta = np.mod(np.arctan2(coords[:, 1], coords[:, 0]), 2.0 * np.pi)
+    order = np.argsort(theta)
+    theta = theta[order]
+    values = values[order]
+    theta = np.append(theta, theta[0] + 2.0 * np.pi)
+    values = np.append(values, values[0])
+
+    coefficient = (
+        _trapezoidal_integral(
+            values * np.cos(mode * theta),
+            theta,
+        )
+        / np.pi
+    )
+    mean = _trapezoidal_integral(values, theta) / (2.0 * np.pi)
+    return float(coefficient), float(mean)
+
+
+def _merge_unique_boundary_samples(gathered_rows, dimension: int):
+    """Merge duplicate partition-boundary samples on the root rank."""
+
+    nonempty_rows = [rows for rows in gathered_rows if rows.size]
+    if not nonempty_rows:
+        raise RuntimeError("No boundary samples were recovered.")
+
+    merged = {}
+    for row in np.vstack(nonempty_rows):
+        key = tuple(np.round(row[:dimension], 12))
+        if key not in merged:
+            merged[key] = [row[:dimension].copy(), 0.0, 0]
+        merged[key][1] += float(row[dimension])
+        merged[key][2] += 1
+    coords = np.array([entry[0] for entry in merged.values()])
+    values = np.array(
+        [entry[1] / entry[2] for entry in merged.values()],
+        dtype=float,
+    )
+    return coords, values
+
+
+def _rotated_cylindrical_boundary_response(
+    *,
+    stokes,
+    boundary: str,
+    wavenumber: int,
+    buoyancy_scale: float,
+) -> tuple[float, float, float]:
+    """Recover one circular boundary's reaction and topography coefficients."""
+
+    mode = _validate_cylindrical_mode(wavenumber)
+    buoyancy_scale = _finite_float(buoyancy_scale, "buoyancy_scale")
+    if buoyancy_scale == 0.0:
+        raise ValueError("Boundary buoyancy scales must be nonzero.")
+    if getattr(stokes.mesh, "dim", None) != 2:
+        raise ValueError("The cylindrical adapter requires a two-dimensional mesh.")
+
+    local_error = None
+    try:
+        local_coords, local_reaction = stokes.boundary_normal_traction(
+            boundary,
+            mass="auto",
+        )
+        local_coords = np.asarray(local_coords, dtype=float)
+        local_reaction = np.asarray(local_reaction, dtype=float).reshape(-1)
+        if local_coords.ndim != 2 or local_coords.shape[1] != 2:
+            raise ValueError("Boundary coordinates must have shape (n, 2).")
+        if local_coords.shape[0] != local_reaction.size:
+            raise ValueError(
+                "Boundary coordinates and reactions must have equal lengths."
+            )
+        local_rows = np.column_stack((local_coords, local_reaction))
+    except Exception as error:
+        local_rows = np.empty((0, 3), dtype=float)
+        local_error = f"{type(error).__name__}: {error}"
+
+    rank_errors = MPI.COMM_WORLD.allgather(local_error)
+    rank_errors = [error for error in rank_errors if error is not None]
+    if rank_errors:
+        raise RuntimeError(
+            "Cylindrical boundary traction recovery failed: " + rank_errors[0]
+        )
+
+    gathered_rows = MPI.COMM_WORLD.gather(local_rows, root=0)
+    result = None
+    root_error = None
+    if MPI.COMM_WORLD.rank == 0:
+        try:
+            coords, reaction = _merge_unique_boundary_samples(
+                gathered_rows,
+                dimension=2,
+            )
+            reaction_coefficient, reaction_mean = (
+                cylindrical_cosine_boundary_coefficient(
+                    coords,
+                    reaction,
+                    mode,
+                )
+            )
+            result = (
+                reaction_coefficient,
+                reaction_mean,
+                -reaction_coefficient / buoyancy_scale,
+            )
+        except Exception as error:
+            root_error = f"{type(error).__name__}: {error}"
+
+    root_error = MPI.COMM_WORLD.bcast(root_error, root=0)
+    if root_error is not None:
+        raise RuntimeError("Cylindrical Fourier projection failed: " + root_error)
+    global_result = MPI.COMM_WORLD.bcast(result, root=0)
+    if global_result is None:
+        raise RuntimeError("MPI broadcast returned no cylindrical response.")
+    return (
+        float(global_result[0]),
+        float(global_result[1]),
+        float(global_result[2]),
+    )
+
+
+def cylindrical_annulus_response_from_rotated_stokes(
+    *,
+    stokes,
+    radius_inner: float,
+    radius_outer: float,
+    wavenumber: int,
+    outer_density_contrast: float,
+    inner_density_contrast: float,
+    outer_reference_gravity: float,
+    inner_reference_gravity: float,
+    internal_load_radius: float | None = None,
+    internal_surface_density_coefficient: float = 0.0,
+    outer_boundary: str = "Upper",
+    inner_boundary: str = "Lower",
+    outer_buoyancy_scale: float = 1.0,
+    inner_buoyancy_scale: float = -1.0,
+    gravitational_constant: float = 1.0,
+    include_self_gravity: bool = False,
+    outer_feedback_factor: float | None = None,
+    inner_feedback_factor: float | None = None,
+) -> CylindricalAnnulusResponse:
+    r"""Compute a cylindrical response from rotated-free-slip wall reactions.
+
+    :meth:`Stokes.boundary_normal_traction` returns the wall reaction. This
+    adapter uses
+
+    .. math::
+
+       h=-reaction_{nn}/signed\_buoyancy\_scale
+
+    and projects each circular boundary onto the unnormalised real basis
+    :math:`\cos(n\theta)`. Boundary samples are gathered only to MPI rank zero;
+    the resulting coefficients are broadcast to all ranks.
+
+    Parameters
+    ----------
+    stokes : underworld3.systems.Stokes
+        Completed two-dimensional rotated-free-slip Stokes solve.
+    radius_inner, radius_outer : float
+        Positive annulus radii ordered from inner to outer.
+    wavenumber : int
+        Positive azimuthal Fourier wavenumber.
+    outer_density_contrast, inner_density_contrast : float
+        Signed density contrasts at the two boundaries.
+    outer_reference_gravity, inner_reference_gravity : float
+        Positive gravity magnitudes used to convert potential to geoid.
+    internal_load_radius : float, optional
+        Radius of an internal sheet, strictly inside the annulus.
+    internal_surface_density_coefficient : float, default=0
+        Density coefficient of the optional internal sheet.
+    outer_boundary, inner_boundary : str
+        Mesh boundary labels used for reaction recovery.
+    outer_buoyancy_scale, inner_buoyancy_scale : float
+        Signed scales converting wall reaction to dynamic topography.
+    gravitational_constant : float, default=1
+        Positive gravitational constant in the selected unit system.
+    include_self_gravity : bool, default=False
+        Return the self-gravity-corrected response when true.
+    outer_feedback_factor, inner_feedback_factor : float, optional
+        Explicit self-gravity feedback factors.
+
+    Returns
+    -------
+    CylindricalAnnulusResponse
+        Recovered reactions, topography, gravity response, and optional
+        self-gravity correction, identical on every MPI rank.
+    """
+
+    radius_inner, radius_outer, mode = _validate_cylindrical_annulus(
+        radius_inner,
+        radius_outer,
+        wavenumber,
+    )
+    if not isinstance(include_self_gravity, bool):
+        raise TypeError("include_self_gravity must be True or False.")
+
+    outer_reaction, outer_mean, outer_topography = (
+        _rotated_cylindrical_boundary_response(
+            stokes=stokes,
+            boundary=outer_boundary,
+            wavenumber=mode,
+            buoyancy_scale=outer_buoyancy_scale,
+        )
+    )
+    inner_reaction, inner_mean, inner_topography = (
+        _rotated_cylindrical_boundary_response(
+            stokes=stokes,
+            boundary=inner_boundary,
+            wavenumber=mode,
+            buoyancy_scale=inner_buoyancy_scale,
+        )
+    )
+    gravity = cylindrical_annulus_geoid_response(
+        radius_inner=radius_inner,
+        radius_outer=radius_outer,
+        wavenumber=mode,
+        outer_topography_coefficient=outer_topography,
+        inner_topography_coefficient=inner_topography,
+        outer_density_contrast=outer_density_contrast,
+        inner_density_contrast=inner_density_contrast,
+        outer_reference_gravity=outer_reference_gravity,
+        inner_reference_gravity=inner_reference_gravity,
+        internal_load_radius=internal_load_radius,
+        internal_surface_density_coefficient=internal_surface_density_coefficient,
+        gravitational_constant=gravitational_constant,
+    )
+    self_gravity = None
+    if include_self_gravity:
+        self_gravity = cylindrical_annulus_self_gravity_response(
+            radius_inner=radius_inner,
+            radius_outer=radius_outer,
+            wavenumber=mode,
+            outer_topography_coefficient=outer_topography,
+            inner_topography_coefficient=inner_topography,
+            outer_density_contrast=outer_density_contrast,
+            inner_density_contrast=inner_density_contrast,
+            outer_reference_gravity=outer_reference_gravity,
+            inner_reference_gravity=inner_reference_gravity,
+            internal_load_radius=internal_load_radius,
+            internal_surface_density_coefficient=internal_surface_density_coefficient,
+            gravitational_constant=gravitational_constant,
+            outer_feedback_factor=outer_feedback_factor,
+            inner_feedback_factor=inner_feedback_factor,
+        )
+
+    return CylindricalAnnulusResponse(
+        outer_reaction=outer_reaction,
+        inner_reaction=inner_reaction,
+        outer_reaction_mean=outer_mean,
+        inner_reaction_mean=inner_mean,
+        outer_topography=outer_topography,
+        inner_topography=inner_topography,
+        outer_potential=gravity.outer_potential,
+        inner_potential=gravity.inner_potential,
+        outer_geoid=gravity.outer_geoid,
+        inner_geoid=gravity.inner_geoid,
+        self_gravity=self_gravity,
+    )
diff --git a/tests/test_1071_postprocessing_cylindrical_geoid.py b/tests/test_1071_postprocessing_cylindrical_geoid.py
new file mode 100644
index 000000000..9083d3f16
--- /dev/null
+++ b/tests/test_1071_postprocessing_cylindrical_geoid.py
@@ -0,0 +1,339 @@
+from types import SimpleNamespace
+
+from mpi4py import MPI
+import numpy as np
+import pytest
+import underworld3 as uw
+
+
+pytestmark = [pytest.mark.level_2, pytest.mark.tier_c]
+
+RADIUS_INNER = 1.22
+RADIUS_INTERNAL = 2.0
+RADIUS_OUTER = 2.22
+GRAVITATIONAL_CONSTANT = 3.7
+
+
+def _circle_samples(radius, coefficient, mean, wavenumber, count=256):
+    theta = np.linspace(0.0, 2.0 * np.pi, count, endpoint=False)
+    coords = radius * np.column_stack((np.cos(theta), np.sin(theta)))
+    values = coefficient * np.cos(wavenumber * theta) + mean
+    return coords, values
+
+
+class PartitionedFakeRotatedStokes:
+    """Return a disjoint subset of each synthetic boundary on every rank."""
+
+    def __init__(self, boundary_data):
+        self.mesh = SimpleNamespace(dim=2)
+        self.boundary_data = boundary_data
+
+    def boundary_normal_traction(self, boundary, mass="auto"):
+        assert mass == "auto"
+        coords, values = self.boundary_data[boundary]
+        indices = np.arange(coords.shape[0])[MPI.COMM_WORLD.rank :: MPI.COMM_WORLD.size]
+        return coords[indices], values[indices]
+
+
+def _adapter_kwargs(wavenumber=2):
+    return {
+        "radius_inner": RADIUS_INNER,
+        "radius_outer": RADIUS_OUTER,
+        "wavenumber": wavenumber,
+        "outer_density_contrast": 0.06,
+        "inner_density_contrast": 0.09,
+        "outer_reference_gravity": 1.7,
+        "inner_reference_gravity": 2.4,
+        "internal_load_radius": RADIUS_INTERNAL,
+        "internal_surface_density_coefficient": 0.027,
+        "outer_buoyancy_scale": 1.0,
+        "inner_buoyancy_scale": -1.0,
+        "gravitational_constant": 0.1,
+    }
+
+
+def _pure_gravity_kwargs(adapter_kwargs):
+    return {
+        key: value
+        for key, value in adapter_kwargs.items()
+        if key not in ("outer_buoyancy_scale", "inner_buoyancy_scale")
+    }
+
+
+def test_cylindrical_geoid_api_is_public():
+    expected = (
+        "CylindricalGravityResponse",
+        "CylindricalSelfGravityResponse",
+        "CylindricalAnnulusResponse",
+        "cylindrical_sheet_potential_coefficient",
+        "cylindrical_sheet_radial_derivative_coefficient",
+        "cylindrical_annulus_potential_operator",
+        "cylindrical_annulus_geoid_response",
+        "cylindrical_annulus_self_gravity_response",
+        "cylindrical_cosine_boundary_coefficient",
+        "cylindrical_annulus_response_from_rotated_stokes",
+    )
+    for name in expected:
+        assert name in uw.postprocessing.geoid.__all__
+        assert hasattr(uw.postprocessing.geoid, name)
+
+
+@pytest.mark.parametrize("wavenumber", [1, 2, 3, 4, 8])
+def test_sheet_kernel_supports_positive_fourier_modes(wavenumber):
+    density = -0.73
+    source = uw.postprocessing.geoid.cylindrical_sheet_potential_coefficient(
+        source_radius=RADIUS_INTERNAL,
+        target_radius=RADIUS_INTERNAL,
+        wavenumber=wavenumber,
+        surface_density_coefficient=density,
+        gravitational_constant=GRAVITATIONAL_CONSTANT,
+    )
+    expected = (
+        2.0 * np.pi * GRAVITATIONAL_CONSTANT * RADIUS_INTERNAL * density / wavenumber
+    )
+    assert source == pytest.approx(expected)
+
+    inner_radius = 0.4 * RADIUS_INTERNAL
+    outer_radius = 3.0 * RADIUS_INTERNAL
+    inner = uw.postprocessing.geoid.cylindrical_sheet_potential_coefficient(
+        source_radius=RADIUS_INTERNAL,
+        target_radius=inner_radius,
+        wavenumber=wavenumber,
+        surface_density_coefficient=density,
+        gravitational_constant=GRAVITATIONAL_CONSTANT,
+    )
+    outer = uw.postprocessing.geoid.cylindrical_sheet_potential_coefficient(
+        source_radius=RADIUS_INTERNAL,
+        target_radius=outer_radius,
+        wavenumber=wavenumber,
+        surface_density_coefficient=density,
+        gravitational_constant=GRAVITATIONAL_CONSTANT,
+    )
+    assert inner == pytest.approx(
+        source * (inner_radius / RADIUS_INTERNAL) ** wavenumber
+    )
+    assert outer == pytest.approx(
+        source * (RADIUS_INTERNAL / outer_radius) ** wavenumber
+    )
+
+
+@pytest.mark.parametrize("wavenumber", [1, 2, 3, 4, 8])
+def test_sheet_derivative_has_poisson_jump(wavenumber):
+    density = 0.41
+    common = {
+        "source_radius": RADIUS_INTERNAL,
+        "target_radius": RADIUS_INTERNAL,
+        "wavenumber": wavenumber,
+        "surface_density_coefficient": density,
+        "gravitational_constant": GRAVITATIONAL_CONSTANT,
+    }
+    derivative_inside = (
+        uw.postprocessing.geoid.cylindrical_sheet_radial_derivative_coefficient(
+            source_side="inside",
+            **common,
+        )
+    )
+    derivative_outside = (
+        uw.postprocessing.geoid.cylindrical_sheet_radial_derivative_coefficient(
+            source_side="outside",
+            **common,
+        )
+    )
+    expected_jump = -4.0 * np.pi * GRAVITATIONAL_CONSTANT * density
+    assert derivative_outside - derivative_inside == pytest.approx(expected_jump)
+
+
+def test_axisymmetric_mode_requires_separate_logarithmic_solution():
+    with pytest.raises(ValueError, match="logarithmic radial solution"):
+        uw.postprocessing.geoid.cylindrical_sheet_potential_coefficient(
+            source_radius=2.0,
+            target_radius=2.2,
+            wavenumber=0,
+            surface_density_coefficient=1.0,
+        )
+
+
+def test_negative_mode_is_rejected():
+    with pytest.raises(ValueError, match="non-negative"):
+        uw.postprocessing.geoid.cylindrical_sheet_potential_coefficient(
+            source_radius=2.0,
+            target_radius=2.2,
+            wavenumber=-1,
+            surface_density_coefficient=1.0,
+        )
+
+
+@pytest.mark.parametrize("wavenumber", [True, 2.0])
+def test_noninteger_modes_are_rejected(wavenumber):
+    with pytest.raises(TypeError, match="must be an integer"):
+        uw.postprocessing.geoid.cylindrical_sheet_potential_coefficient(
+            source_radius=2.0,
+            target_radius=2.2,
+            wavenumber=wavenumber,
+            surface_density_coefficient=1.0,
+        )
+
+
+def test_annulus_geoid_response_superposes_boundaries_and_internal_load():
+    mode = 3
+    outer_topography = -0.77
+    inner_topography = -0.32
+    outer_density = 0.6
+    inner_density = -0.9
+    internal_density = 0.27
+    outer_gravity = 1.7
+    inner_gravity = 2.4
+    operator = uw.postprocessing.geoid.cylindrical_annulus_potential_operator(
+        radius_inner=RADIUS_INNER,
+        radius_outer=RADIUS_OUTER,
+        wavenumber=mode,
+        outer_density_contrast=outer_density,
+        inner_density_contrast=inner_density,
+        gravitational_constant=GRAVITATIONAL_CONSTANT,
+    )
+    expected_potential = operator @ np.array([outer_topography, inner_topography])
+    for index, target_radius in enumerate((RADIUS_OUTER, RADIUS_INNER)):
+        expected_potential[
+            index
+        ] += uw.postprocessing.geoid.cylindrical_sheet_potential_coefficient(
+            source_radius=RADIUS_INTERNAL,
+            target_radius=target_radius,
+            wavenumber=mode,
+            surface_density_coefficient=internal_density,
+            gravitational_constant=GRAVITATIONAL_CONSTANT,
+        )
+
+    response = uw.postprocessing.geoid.cylindrical_annulus_geoid_response(
+        radius_inner=RADIUS_INNER,
+        radius_outer=RADIUS_OUTER,
+        wavenumber=mode,
+        outer_topography_coefficient=outer_topography,
+        inner_topography_coefficient=inner_topography,
+        outer_density_contrast=outer_density,
+        inner_density_contrast=inner_density,
+        outer_reference_gravity=outer_gravity,
+        inner_reference_gravity=inner_gravity,
+        internal_load_radius=RADIUS_INTERNAL,
+        internal_surface_density_coefficient=internal_density,
+        gravitational_constant=GRAVITATIONAL_CONSTANT,
+    )
+    np.testing.assert_allclose(
+        [response.outer_potential, response.inner_potential],
+        expected_potential,
+    )
+    np.testing.assert_allclose(
+        [response.outer_geoid, response.inner_geoid],
+        expected_potential / np.array([outer_gravity, inner_gravity]),
+    )
+
+
+def test_cylindrical_self_gravity_satisfies_matrix_equation():
+    kwargs = _adapter_kwargs(wavenumber=4)
+    pure_kwargs = _pure_gravity_kwargs(kwargs)
+    original_topography = np.array([-0.77, -0.32])
+    response = uw.postprocessing.geoid.cylindrical_annulus_self_gravity_response(
+        **pure_kwargs,
+        outer_topography_coefficient=original_topography[0],
+        inner_topography_coefficient=original_topography[1],
+    )
+    operator = uw.postprocessing.geoid.cylindrical_annulus_potential_operator(
+        radius_inner=kwargs["radius_inner"],
+        radius_outer=kwargs["radius_outer"],
+        wavenumber=kwargs["wavenumber"],
+        outer_density_contrast=kwargs["outer_density_contrast"],
+        inner_density_contrast=kwargs["inner_density_contrast"],
+        gravitational_constant=kwargs["gravitational_constant"],
+    )
+    load = np.array(
+        [
+            uw.postprocessing.geoid.cylindrical_sheet_potential_coefficient(
+                source_radius=kwargs["internal_load_radius"],
+                target_radius=target_radius,
+                wavenumber=kwargs["wavenumber"],
+                surface_density_coefficient=kwargs[
+                    "internal_surface_density_coefficient"
+                ],
+                gravitational_constant=kwargs["gravitational_constant"],
+            )
+            for target_radius in (RADIUS_OUTER, RADIUS_INNER)
+        ]
+    )
+    q_matrix = np.diag([response.q_outer, response.q_inner])
+    corrected_topography = np.array(
+        [response.outer_topography, response.inner_topography]
+    )
+    residual = (
+        (np.eye(2) - q_matrix @ operator) @ corrected_topography
+        - original_topography
+        - q_matrix @ load
+    )
+    np.testing.assert_allclose(residual, 0.0, atol=2.0e-16)
+    assert response.matrix_residual_norm < 2.0e-16
+
+
+def test_cylindrical_projection_recovers_mode_and_mean():
+    coords, values = _circle_samples(2.22, -0.73, 0.12, 8)
+    coefficient, mean = uw.postprocessing.geoid.cylindrical_cosine_boundary_coefficient(
+        coords,
+        values,
+        8,
+    )
+    assert coefficient == pytest.approx(-0.73, abs=2.0e-15)
+    assert mean == pytest.approx(0.12, abs=2.0e-15)
+
+
+def test_partitioned_rotated_adapter_matches_pure_response():
+    kwargs = _adapter_kwargs(wavenumber=3)
+    outer_reaction = 0.8
+    inner_reaction = -0.4
+    stokes = PartitionedFakeRotatedStokes(
+        {
+            "Upper": _circle_samples(
+                RADIUS_OUTER,
+                outer_reaction,
+                0.03,
+                kwargs["wavenumber"],
+            ),
+            "Lower": _circle_samples(
+                RADIUS_INNER,
+                inner_reaction,
+                -0.02,
+                kwargs["wavenumber"],
+            ),
+        }
+    )
+    response = uw.postprocessing.geoid.cylindrical_annulus_response_from_rotated_stokes(
+        stokes=stokes,
+        include_self_gravity=True,
+        **kwargs,
+    )
+    expected_outer_topography = -outer_reaction / kwargs["outer_buoyancy_scale"]
+    expected_inner_topography = -inner_reaction / kwargs["inner_buoyancy_scale"]
+    pure_kwargs = _pure_gravity_kwargs(kwargs)
+    expected = uw.postprocessing.geoid.cylindrical_annulus_geoid_response(
+        **pure_kwargs,
+        outer_topography_coefficient=expected_outer_topography,
+        inner_topography_coefficient=expected_inner_topography,
+    )
+
+    assert response.outer_reaction == pytest.approx(outer_reaction)
+    assert response.inner_reaction == pytest.approx(inner_reaction)
+    np.testing.assert_allclose(
+        [response.outer_topography, response.inner_topography],
+        [expected_outer_topography, expected_inner_topography],
+    )
+    np.testing.assert_allclose(
+        [
+            response.outer_potential,
+            response.inner_potential,
+            response.outer_geoid,
+            response.inner_geoid,
+        ],
+        [
+            expected.outer_potential,
+            expected.inner_potential,
+            expected.outer_geoid,
+            expected.inner_geoid,
+        ],
+    )
+    assert response.self_gravity is not None

From c2cff1a2568a92dd2ecce785c9923ff197249d16 Mon Sep 17 00:00:00 2001
From: Tyagi <gthyagi@gmail.com>
Date: Wed, 26 Aug 2026 08:43:19 +0530
Subject: [PATCH 2/3] Document cylindrical annulus gravity postprocessing

Expose the postprocessing API in the documentation index and describe the Fourier convention, sheet kernel, density and sign conventions, self-gravity matrix, unit contract, n=0 limitation, and MPI adapter workflow.

Keep the Kramer-Simons convergence and physical-space Poisson comparisons in the mantle-convection benchmark repository.
---
 docs/api/index.md                             |  4 +
 docs/api/postprocessing.md                    | 19 +++++
 ...ry-stress-and-projection-postprocessing.md | 82 +++++++++++++++++++
 3 files changed, 105 insertions(+)
 create mode 100644 docs/api/postprocessing.md

diff --git a/docs/api/index.md b/docs/api/index.md
index a0d9602bf..f7489d2cd 100644
--- a/docs/api/index.md
+++ b/docs/api/index.md
@@ -30,6 +30,7 @@ utilities
 visualisation
 adaptivity
 analytic
+postprocessing
 ```
 
 ## Quick Links
@@ -54,6 +55,9 @@ analytic
 ### Validation
 - **{doc}`analytic`** - Exact solutions for benchmarking and convergence testing
 
+### Post-processing
+- **{doc}`postprocessing`** - Boundary-response, geoid, and self-gravity coefficients
+
 ### Infrastructure
 - **{doc}`model`** - Model management and configuration
 - **{doc}`utilities`** - I/O, mesh import, and helper functions
diff --git a/docs/api/postprocessing.md b/docs/api/postprocessing.md
new file mode 100644
index 000000000..07512bcc0
--- /dev/null
+++ b/docs/api/postprocessing.md
@@ -0,0 +1,19 @@
+# Post-processing
+
+```{eval-rst}
+.. automodule:: underworld3.postprocessing
+   :members:
+   :show-inheritance:
+```
+
+## Geoid and self-gravity responses
+
+The geoid module provides coefficient-only spherical-shell and
+cylindrical-annulus gravity operators plus adapters for completed
+rotated-free-slip Stokes solves.
+
+```{eval-rst}
+.. automodule:: underworld3.postprocessing.geoid
+   :members:
+   :show-inheritance:
+```
diff --git a/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md b/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md
index fb007f0ba..e65ce5715 100644
--- a/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md
+++ b/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md
@@ -142,6 +142,88 @@ reconstructs their spherical triangulation. A future boundary-reaction
 functional could replace this step with a direct distributed finite-element
 projection without changing the coefficient API.
 
+## 4. Cylindrical-annulus gravity and geoid response
+
+The cylindrical API uses the unnormalised real Fourier basis
+`cos(n theta)`. For a sheet-density coefficient `sigma_n` at radius `r_s`,
+the convention is
+
+```text
+laplacian(Phi) = -4*pi*G*rho
+gravity = grad(Phi)
+Phi_n(r_s) = 2*pi*G*r_s*sigma_n/n
+```
+
+The coefficient varies as `(r/r_s)^n` inside the sheet and `(r_s/r)^n`
+outside it. These branches are regular toward the axis and decay at infinity.
+They apply for integer modes `n >= 1`. The axisymmetric `n=0` solution is
+logarithmic and requires an explicit potential gauge, so this API rejects it.
+
+When topography coefficients are already available, use the pure operator:
+
+```python
+response = uw.postprocessing.geoid.cylindrical_annulus_geoid_response(
+    radius_inner=1.22,
+    radius_outer=2.22,
+    wavenumber=2,
+    outer_topography_coefficient=-0.77,
+    inner_topography_coefficient=-0.32,
+    outer_density_contrast=0.06,
+    inner_density_contrast=0.09,
+    outer_reference_gravity=1.7,
+    inner_reference_gravity=2.4,
+    internal_load_radius=2.0,
+    internal_surface_density_coefficient=0.027,
+    gravitational_constant=0.1,
+)
+```
+
+Potential and topography keep their physical signs; geoid is returned as
+`Phi_n/g_reference` independently at both boundaries. Radii, topography,
+sheet density, gravity, and the gravitational constant may be dimensional or
+nondimensional, but every input must use one consistent unit system. Density
+contrast is defined as the smaller-radius density minus the larger-radius
+density, so positive outward topography creates sheet density
+`Delta_rho*h`.
+
+Self-gravity solves the two-boundary coefficient equation
+
+```text
+(I - Q G_n) h_self_gravity = h + Q phi_load
+Q = diag(1/g_outer, 1/g_inner)
+```
+
+with `cylindrical_annulus_self_gravity_response()`. Explicit feedback factors
+can disable either row or represent another signed convention.
+
+For a completed two-dimensional rotated-free-slip Stokes solve, the adapter
+recovers the wall reactions, defines
+`h=-reaction_nn/signed_buoyancy_scale`, and performs the Fourier projection:
+
+```python
+response = (
+    uw.postprocessing.geoid.cylindrical_annulus_response_from_rotated_stokes(
+        stokes=stokes,
+        radius_inner=1.22,
+        radius_outer=2.22,
+        wavenumber=2,
+        outer_density_contrast=0.06,
+        inner_density_contrast=0.09,
+        outer_reference_gravity=1.7,
+        inner_reference_gravity=2.4,
+        outer_buoyancy_scale=1.0,
+        inner_buoyancy_scale=-1.0,
+        include_self_gravity=True,
+    )
+)
+```
+
+Only boundary samples are gathered to rank zero; the projected coefficients
+are broadcast to all ranks. The coefficient kernel follows Simons (1996),
+Appendix B. Complete Kramer--Simons finite-element convergence and
+physical-space Poisson comparisons remain in the separate mantle-convection
+benchmark repository.
+
 ## See also
 
 - Issues [#156] (projection solver settings), [#157] (projection memory),

From b30a02bc73f957ccf694cf1ecb7895a340e0c143 Mon Sep 17 00:00:00 2001
From: Tyagi <gthyagi@gmail.com>
Date: Thu, 3 Sep 2026 16:14:55 +1000
Subject: [PATCH 3/3] Project annulus reactions with distributed boundary
 integrals

Expose Stokes.boundary_normal_traction_integral() as a weak contraction of the assembled rotated-free-slip reaction. Count only owned reaction degrees of freedom, reduce on the mesh communicator, accept numeric or symbolic scalar tests, and remove the constant traction gauge through boundary integrals.

Use this primitive for cylindrical-annulus Fourier projection and normalize with the matching BdIntegral harmonic norm. This removes pointwise traction recovery, angular sorting, COMM_WORLD collection, and the rank-zero gather from the finite-element adapter while retaining the sampled-circle helper for external data.

Add serial and MPI coverage with a real annulus Stokes solve, document the fitted projection, and renumber the cylindrical geoid test from 1071 to 1073 to avoid the development-branch collision.
---
 docs/developer/CHANGELOG.md                   |   5 +
 ...ry-stress-and-projection-postprocessing.md |  31 ++--
 .../cython/petsc_generic_snes_solvers.pyx     |  19 +++
 src/underworld3/postprocessing/geoid.py       | 134 ++++++---------
 src/underworld3/utilities/rotated_bc.py       |  83 ++++++++++
 ..._1073_postprocessing_cylindrical_geoid.py} | 152 +++++++++++++++---
 6 files changed, 309 insertions(+), 115 deletions(-)
 rename tests/{test_1071_postprocessing_cylindrical_geoid.py => test_1073_postprocessing_cylindrical_geoid.py} (70%)

diff --git a/docs/developer/CHANGELOG.md b/docs/developer/CHANGELOG.md
index 53c0126f8..d92498873 100644
--- a/docs/developer/CHANGELOG.md
+++ b/docs/developer/CHANGELOG.md
@@ -278,6 +278,11 @@ component exactly — correct on curved, tilted, and deformed boundaries (#293).
   existing boundary traction onto an axisymmetric harmonic; the pure functions
   also accept coefficients recovered by other methods and an optional internal
   load.
+- Rotated free slip exposes
+  `Stokes.boundary_normal_traction_integral(boundary, fn)` for a distributed
+  weak contraction of the assembled normal reaction. Cylindrical-annulus
+  Stokes responses use this fitted integral and its matching finite-element
+  boundary norm instead of gathering pointwise samples for angular quadrature.
 - `uw.analytic.Zhong2008` implements the Hager--O'Connell propagator-matrix
   oracle used for the Zhong et al. spherical-shell response benchmark. It
   supports piecewise-constant radial viscosity and reproduces every analytical
diff --git a/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md b/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md
index e65ce5715..6aef79807 100644
--- a/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md
+++ b/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md
@@ -137,10 +137,13 @@ Published reference solvers, such as the Zhong et al. propagator-matrix method,
 belong in `uw.analytic`; their computed topography coefficients can be passed to
 the pure post-processing functions above.
 
-The rotated harmonic projector gathers boundary samples to rank zero and
-reconstructs their spherical triangulation. A future boundary-reaction
-functional could replace this step with a direct distributed finite-element
-projection without changing the coefficient API.
+`Stokes.boundary_normal_traction_integral(boundary, fn)` contracts an assembled
+normal-reaction load directly with a scalar test function over owned degrees of
+freedom, followed by reductions on the mesh communicator. It is useful whenever
+only an integrated or fitted traction diagnostic is required. It avoids
+pointwise recovery and global boundary reconstruction; consumers that need a
+nodal field should continue to use `boundary_normal_traction()` or
+`dynamic_topography()` and follow their curved-P2 guidance.
 
 ## 4. Cylindrical-annulus gravity and geoid response
 
@@ -218,11 +221,21 @@ response = (
 )
 ```
 
-Only boundary samples are gathered to rank zero; the projected coefficients
-are broadcast to all ranks. The coefficient kernel follows Simons (1996),
-Appendix B. Complete Kramer--Simons finite-element convergence and
-physical-space Poisson comparisons remain in the separate mantle-convection
-benchmark repository.
+The Stokes adapter contracts the assembled reaction directly with
+`cos(n theta)` through `boundary_normal_traction_integral()`. The coefficient is
+normalised by the matching `BdIntegral` of `cos(n theta)**2`, so its numerator
+and denominator use the same faceted finite-element boundary geometry. The
+calculation counts owned reaction degrees of freedom and reduces on the mesh
+communicator; it does not recover pointwise traction, sort samples by angle, or
+gather a global boundary on rank zero. This construction remains valid on a
+deformed boundary, where the fitted coefficient is weighted by the actual
+finite-element boundary measure.
+
+`cylindrical_cosine_boundary_coefficient()` remains available for external
+sampled data on a mathematical circle; it is not used for finite-element
+reaction loads. The coefficient kernel follows Simons (1996), Appendix B.
+Complete Kramer--Simons finite-element convergence and physical-space Poisson
+comparisons remain in the separate mantle-convection benchmark repository.
 
 ## See also
 
diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
index 3b3da609f..712879606 100644
--- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx
+++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
@@ -6502,6 +6502,25 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         return _dtf(self, boundary, self._rotated_freeslip_info, field,
                     buoyancy_scale=buoyancy_scale, mass=mass)
 
+    def boundary_normal_traction_integral(self, boundary, fn, remove_mean=True):
+        r"""Return the boundary integral of ``sigma_nn * fn`` directly from the
+        rotated-free-slip constraint reaction.
+
+        With ``remove_mean=True`` (default), the constant normal-traction gauge
+        is removed before projection. Unlike pointwise
+        :meth:`boundary_normal_traction`, this weak projection does not recover
+        nodal traction values or gather a global boundary mesh. It is therefore
+        suitable for harmonic and integral diagnostics on curved P2 boundaries,
+        whose recovered vertex values converge slowly (issue #414).
+        """
+        if self._rotated_freeslip_info is None:
+            raise RuntimeError(
+                "boundary_normal_traction_integral requires a completed "
+                "rotated-free-slip solve.")
+        from underworld3.utilities.rotated_bc import boundary_normal_traction_integral as _bnti
+        return _bnti(self, boundary, self._rotated_freeslip_info, fn,
+                     remove_mean=remove_mean)
+
     def add_nitsche_bc(self, conds=None, boundary=None, direction=None, normal=None,
                        gamma=10.0, theta=1, mask=None, local_h=True, g=None):
         r"""Add Nitsche weak enforcement of a velocity constraint along a direction.
diff --git a/src/underworld3/postprocessing/geoid.py b/src/underworld3/postprocessing/geoid.py
index c3bb6c0b9..3a2e18eb6 100644
--- a/src/underworld3/postprocessing/geoid.py
+++ b/src/underworld3/postprocessing/geoid.py
@@ -1066,11 +1066,15 @@ def cylindrical_cosine_boundary_coefficient(
     values,
     wavenumber: int,
 ) -> tuple[float, float]:
-    r"""Project circular samples onto :math:`\cos(n\theta)` and the mean.
+    r"""Project sampled circular data onto :math:`\cos(n\theta)` and the mean.
 
     ``coords`` must contain at least three two-dimensional Cartesian boundary
     points. The samples may begin at any angle and need not include a duplicate
-    endpoint; this function sorts and closes the periodic interval.
+    endpoint; this function sorts and closes the periodic interval. This helper
+    is intended for external sampled data. Finite-element reaction loads use
+    :meth:`Stokes.boundary_normal_traction_integral` instead, so their projection
+    follows the actual boundary facets without a global gather or reconstructed
+    angular ordering.
 
     Parameters
     ----------
@@ -1117,28 +1121,6 @@ def cylindrical_cosine_boundary_coefficient(
     return float(coefficient), float(mean)
 
 
-def _merge_unique_boundary_samples(gathered_rows, dimension: int):
-    """Merge duplicate partition-boundary samples on the root rank."""
-
-    nonempty_rows = [rows for rows in gathered_rows if rows.size]
-    if not nonempty_rows:
-        raise RuntimeError("No boundary samples were recovered.")
-
-    merged = {}
-    for row in np.vstack(nonempty_rows):
-        key = tuple(np.round(row[:dimension], 12))
-        if key not in merged:
-            merged[key] = [row[:dimension].copy(), 0.0, 0]
-        merged[key][1] += float(row[dimension])
-        merged[key][2] += 1
-    coords = np.array([entry[0] for entry in merged.values()])
-    values = np.array(
-        [entry[1] / entry[2] for entry in merged.values()],
-        dtype=float,
-    )
-    return coords, values
-
-
 def _rotated_cylindrical_boundary_response(
     *,
     stokes,
@@ -1146,7 +1128,7 @@ def _rotated_cylindrical_boundary_response(
     wavenumber: int,
     buoyancy_scale: float,
 ) -> tuple[float, float, float]:
-    """Recover one circular boundary's reaction and topography coefficients."""
+    """Project one boundary's assembled reaction and topography coefficients."""
 
     mode = _validate_cylindrical_mode(wavenumber)
     buoyancy_scale = _finite_float(buoyancy_scale, "buoyancy_scale")
@@ -1155,66 +1137,48 @@ def _rotated_cylindrical_boundary_response(
     if getattr(stokes.mesh, "dim", None) != 2:
         raise ValueError("The cylindrical adapter requires a two-dimensional mesh.")
 
-    local_error = None
-    try:
-        local_coords, local_reaction = stokes.boundary_normal_traction(
+    import sympy
+    from underworld3.maths import BdIntegral
+
+    theta = stokes.mesh.CoordinateSystem.xR[1]
+    harmonic = sympy.cos(mode * theta)
+    reaction_integral = float(
+        stokes.boundary_normal_traction_integral(
             boundary,
-            mass="auto",
+            harmonic,
+            remove_mean=True,
         )
-        local_coords = np.asarray(local_coords, dtype=float)
-        local_reaction = np.asarray(local_reaction, dtype=float).reshape(-1)
-        if local_coords.ndim != 2 or local_coords.shape[1] != 2:
-            raise ValueError("Boundary coordinates must have shape (n, 2).")
-        if local_coords.shape[0] != local_reaction.size:
-            raise ValueError(
-                "Boundary coordinates and reactions must have equal lengths."
-            )
-        local_rows = np.column_stack((local_coords, local_reaction))
-    except Exception as error:
-        local_rows = np.empty((0, 3), dtype=float)
-        local_error = f"{type(error).__name__}: {error}"
-
-    rank_errors = MPI.COMM_WORLD.allgather(local_error)
-    rank_errors = [error for error in rank_errors if error is not None]
-    if rank_errors:
+    )
+    reaction_total = float(
+        stokes.boundary_normal_traction_integral(
+            boundary,
+            1.0,
+            remove_mean=False,
+        )
+    )
+    boundary_measure = float(
+        BdIntegral(stokes.mesh, fn=1.0, boundary=boundary).evaluate()
+    )
+    harmonic_norm = float(
+        BdIntegral(stokes.mesh, fn=harmonic**2, boundary=boundary).evaluate()
+    )
+    if not np.all(
+        np.isfinite(
+            (reaction_integral, reaction_total, boundary_measure, harmonic_norm)
+        )
+    ):
+        raise RuntimeError("Cylindrical reaction projection produced non-finite data.")
+    if boundary_measure <= 0.0 or harmonic_norm <= 0.0:
         raise RuntimeError(
-            "Cylindrical boundary traction recovery failed: " + rank_errors[0]
+            "Cylindrical reaction projection requires positive boundary integrals."
         )
 
-    gathered_rows = MPI.COMM_WORLD.gather(local_rows, root=0)
-    result = None
-    root_error = None
-    if MPI.COMM_WORLD.rank == 0:
-        try:
-            coords, reaction = _merge_unique_boundary_samples(
-                gathered_rows,
-                dimension=2,
-            )
-            reaction_coefficient, reaction_mean = (
-                cylindrical_cosine_boundary_coefficient(
-                    coords,
-                    reaction,
-                    mode,
-                )
-            )
-            result = (
-                reaction_coefficient,
-                reaction_mean,
-                -reaction_coefficient / buoyancy_scale,
-            )
-        except Exception as error:
-            root_error = f"{type(error).__name__}: {error}"
-
-    root_error = MPI.COMM_WORLD.bcast(root_error, root=0)
-    if root_error is not None:
-        raise RuntimeError("Cylindrical Fourier projection failed: " + root_error)
-    global_result = MPI.COMM_WORLD.bcast(result, root=0)
-    if global_result is None:
-        raise RuntimeError("MPI broadcast returned no cylindrical response.")
+    reaction_coefficient = reaction_integral / harmonic_norm
+    reaction_mean = reaction_total / boundary_measure
     return (
-        float(global_result[0]),
-        float(global_result[1]),
-        float(global_result[2]),
+        float(reaction_coefficient),
+        float(reaction_mean),
+        float(-reaction_coefficient / buoyancy_scale),
     )
 
 
@@ -1241,16 +1205,18 @@ def cylindrical_annulus_response_from_rotated_stokes(
 ) -> CylindricalAnnulusResponse:
     r"""Compute a cylindrical response from rotated-free-slip wall reactions.
 
-    :meth:`Stokes.boundary_normal_traction` returns the wall reaction. This
-    adapter uses
+    :meth:`Stokes.boundary_normal_traction_integral` contracts the assembled
+    wall reaction directly with a boundary test function. This adapter uses
 
     .. math::
 
        h=-reaction_{nn}/signed\_buoyancy\_scale
 
-    and projects each circular boundary onto the unnormalised real basis
-    :math:`\cos(n\theta)`. Boundary samples are gathered only to MPI rank zero;
-    the resulting coefficients are broadcast to all ranks.
+    and projects each boundary onto the unnormalised real basis
+    :math:`\cos(n\theta)`. The numerator and harmonic norm use the same finite-
+    element boundary measure. Owned reaction degrees of freedom are reduced on
+    the mesh communicator; no pointwise traction recovery, angular sorting, or
+    rank-zero boundary gather is performed.
 
     Parameters
     ----------
diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py
index 47a06683a..3f7832567 100644
--- a/src/underworld3/utilities/rotated_bc.py
+++ b/src/underworld3/utilities/rotated_bc.py
@@ -2135,6 +2135,89 @@ def boundary_normal_traction(solver, boundary, solve_result, mass="auto"):
                         remove_mean=True, partial_reaction=False)
 
 
+def boundary_normal_traction_integral(solver, boundary, solve_result, fn,
+                                      remove_mean=True):
+    r"""Return ``integral((sigma_nn - mean) * fn, boundary)`` directly from the
+    assembled rotated-constraint reaction.
+
+    This is the weak/integral counterpart of :func:`boundary_normal_traction`.
+    It contracts the nodal reaction with ``fn`` at the velocity interpolation
+    nodes before any pointwise boundary-mass recovery. On curved P2 boundaries
+    this is a fitted quantity and therefore does not consume the slowly
+    converging recovered vertex values described in issue #414.
+
+    Each assembled reaction degree of freedom is counted on its owning rank,
+    followed by an MPI sum on the mesh communicator. The operation does not
+    gather boundary topology or recovered values onto rank zero.
+    ``remove_mean=True`` removes the constant-traction gauge using boundary
+    integrals of ``fn`` and one.
+    """
+    if not isinstance(remove_mean, (bool, np.bool_)):
+        raise TypeError("remove_mean must be True or False.")
+
+    import underworld3 as uw
+
+    fn = sympy.sympify(fn)
+    dm = solver.dm
+    comm = dm.comm.tompi4py()
+    dim = solver.mesh.dim
+    rc = solve_result["reaction"]
+    rstart, rend = rc.getOwnershipRange()
+    rcl = dm.getLocalVec()
+    dm.globalToLocal(rc, rcl)
+
+    try:
+        rca = np.asarray(rcl.getArray())
+        lsec = dm.getLocalSection()
+        l2g = dm.getLGMap()
+        csec = dm.getCoordinateSection()
+        cvec = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, dim)
+        v0, v1 = dm.getDepthStratum(0)
+        normal = dict(_boundary_spec(s) for s in solve_result["boundaries"]).get(
+            boundary
+        )
+        nodes = _boundary_velocity_nodes(solver, boundary, normal=normal)
+
+        owned_coords = []
+        owned_reactions = []
+        for q, nrm in nodes:
+            lo = lsec.getFieldOffset(q, _VELOCITY_FIELD)
+            global_row = int(l2g.apply([lo])[0])
+            if not rstart <= global_row < rend:
+                continue
+            owned_coords.append(_point_coord(dm, dim, cvec, csec, v0, v1, q))
+            # sigma_nn load = -n.r_c, matching boundary_normal_traction().
+            owned_reactions.append(-float(np.dot(nrm, rca[lo:lo + dim])))
+    finally:
+        dm.restoreLocalVec(rcl)
+
+    if owned_coords:
+        coords = np.ascontiguousarray(owned_coords, dtype=float)
+        weights = np.asarray(uw.function.evaluate(fn, coords), dtype=float).reshape(-1)
+        if weights.size == 1 and len(owned_reactions) != 1:
+            weights = np.full(len(owned_reactions), float(weights[0]))
+        if weights.size != len(owned_reactions):
+            raise ValueError("fn must evaluate to one scalar per boundary node.")
+        local_weighted = float(np.dot(owned_reactions, weights))
+        local_total = float(np.sum(owned_reactions))
+    else:
+        local_weighted = 0.0
+        local_total = 0.0
+
+    weighted = float(comm.allreduce(local_weighted))
+    if not remove_mean:
+        return weighted
+
+    total = float(comm.allreduce(local_total))
+    area = float(uw.maths.BdIntegral(solver.mesh, fn=1.0, boundary=boundary).evaluate())
+    if not np.isfinite(area) or area <= 0.0:
+        raise RuntimeError(f"Boundary {boundary!r} has non-positive area {area}.")
+    fn_integral = float(
+        uw.maths.BdIntegral(solver.mesh, fn=fn, boundary=boundary).evaluate()
+    )
+    return weighted - (total / area) * fn_integral
+
+
 def dynamic_topography_field(solver, boundary, solve_result, field,
                              buoyancy_scale=1.0, mass="auto"):
     """Populate a scalar MeshVariable ``field`` with the dynamic topography
diff --git a/tests/test_1071_postprocessing_cylindrical_geoid.py b/tests/test_1073_postprocessing_cylindrical_geoid.py
similarity index 70%
rename from tests/test_1071_postprocessing_cylindrical_geoid.py
rename to tests/test_1073_postprocessing_cylindrical_geoid.py
index 9083d3f16..8643ae34c 100644
--- a/tests/test_1071_postprocessing_cylindrical_geoid.py
+++ b/tests/test_1073_postprocessing_cylindrical_geoid.py
@@ -1,8 +1,8 @@
 from types import SimpleNamespace
 
-from mpi4py import MPI
 import numpy as np
 import pytest
+import sympy
 import underworld3 as uw
 
 
@@ -21,18 +21,44 @@ def _circle_samples(radius, coefficient, mean, wavenumber, count=256):
     return coords, values
 
 
-class PartitionedFakeRotatedStokes:
-    """Return a disjoint subset of each synthetic boundary on every rank."""
+class IntegralFakeRotatedStokes:
+    """Expose fitted reaction functionals while rejecting pointwise recovery."""
 
     def __init__(self, boundary_data):
-        self.mesh = SimpleNamespace(dim=2)
+        radius, theta = sympy.symbols("r theta", real=True)
+        self.mesh = SimpleNamespace(
+            dim=2,
+            CoordinateSystem=SimpleNamespace(xR=(radius, theta)),
+            boundary_data=boundary_data,
+        )
         self.boundary_data = boundary_data
+        self.integral_calls = []
 
     def boundary_normal_traction(self, boundary, mass="auto"):
-        assert mass == "auto"
-        coords, values = self.boundary_data[boundary]
-        indices = np.arange(coords.shape[0])[MPI.COMM_WORLD.rank :: MPI.COMM_WORLD.size]
-        return coords[indices], values[indices]
+        raise AssertionError("The finite-element adapter must not recover point values.")
+
+    def boundary_normal_traction_integral(self, boundary, fn, remove_mean=True):
+        radius, coefficient, mean = self.boundary_data[boundary]
+        self.integral_calls.append((boundary, bool(remove_mean)))
+        if remove_mean:
+            return coefficient * np.pi * radius
+        assert float(sympy.sympify(fn)) == 1.0
+        return mean * 2.0 * np.pi * radius
+
+
+class FakeBoundaryIntegral:
+    """Exact circular measure used to isolate the adapter's projection route."""
+
+    def __init__(self, mesh, fn, boundary):
+        self.mesh = mesh
+        self.fn = sympy.sympify(fn)
+        self.boundary = boundary
+
+    def evaluate(self):
+        radius = self.mesh.boundary_data[self.boundary][0]
+        if self.fn.is_number and float(self.fn) == 1.0:
+            return 2.0 * np.pi * radius
+        return np.pi * radius
 
 
 def _adapter_kwargs(wavenumber=2):
@@ -76,6 +102,89 @@ def test_cylindrical_geoid_api_is_public():
     for name in expected:
         assert name in uw.postprocessing.geoid.__all__
         assert hasattr(uw.postprocessing.geoid, name)
+    assert hasattr(uw.systems.Stokes, "boundary_normal_traction_integral")
+
+
+def test_reaction_integral_requires_boolean_mean_removal():
+    from underworld3.utilities.rotated_bc import boundary_normal_traction_integral
+
+    with pytest.raises(TypeError, match="remove_mean must be True or False"):
+        boundary_normal_traction_integral(
+            solver=None,
+            boundary="Upper",
+            solve_result=None,
+            fn=1.0,
+            remove_mean="yes",
+        )
+
+
+def test_distributed_reaction_integral_on_finite_element_annulus():
+    """Exercise the complete adapter on an assembled rotated-free-slip reaction."""
+
+    radius_inner = 0.5
+    radius_outer = 1.0
+    mode = 4
+    mesh = uw.meshing.Annulus(
+        radiusInner=radius_inner,
+        radiusOuter=radius_outer,
+        cellSize=0.15,
+        qdegree=3,
+    )
+    x, y = mesh.X
+    radius = sympy.sqrt(x**2 + y**2)
+    theta = sympy.atan2(y, x)
+    velocity = uw.discretisation.MeshVariable(
+        "V_geoid_integral", mesh, mesh.dim, degree=2, continuous=True
+    )
+    pressure = uw.discretisation.MeshVariable(
+        "P_geoid_integral", mesh, 1, degree=1, continuous=True
+    )
+    stokes = uw.systems.Stokes(
+        mesh,
+        velocityField=velocity,
+        pressureField=pressure,
+    )
+    stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
+    stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0
+    radial_force = (
+        sympy.cos(mode * theta)
+        * (radius - radius_inner)
+        * (radius_outer - radius)
+        * 40.0
+    )
+    stokes.bodyforce = sympy.Matrix(
+        [[x * radial_force / radius, y * radial_force / radius]]
+    )
+    normal = sympy.Matrix([[x / radius, y / radius]])
+    stokes.add_rotated_freeslip_bc(0, "Lower", normal=normal)
+    stokes.add_rotated_freeslip_bc(0, "Upper", normal=normal)
+    stokes.petsc_use_pressure_nullspace = True
+    stokes.petsc_options["snes_type"] = "ksponly"
+    stokes.solve()
+
+    response = uw.postprocessing.geoid.cylindrical_annulus_response_from_rotated_stokes(
+        stokes=stokes,
+        radius_inner=radius_inner,
+        radius_outer=radius_outer,
+        wavenumber=mode,
+        outer_density_contrast=1.0,
+        inner_density_contrast=1.0,
+        outer_reference_gravity=1.0,
+        inner_reference_gravity=1.0,
+    )
+
+    assert response.outer_reaction == pytest.approx(0.4044747, rel=2.0e-4)
+    assert response.inner_reaction == pytest.approx(0.5158471, rel=2.0e-4)
+    assert np.all(
+        np.isfinite(
+            [
+                response.outer_reaction_mean,
+                response.inner_reaction_mean,
+                response.outer_geoid,
+                response.inner_geoid,
+            ]
+        )
+    )
 
 
 @pytest.mark.parametrize("wavenumber", [1, 2, 3, 4, 8])
@@ -282,26 +391,17 @@ def test_cylindrical_projection_recovers_mode_and_mean():
     assert mean == pytest.approx(0.12, abs=2.0e-15)
 
 
-def test_partitioned_rotated_adapter_matches_pure_response():
+def test_rotated_adapter_uses_distributed_reaction_integral(monkeypatch):
     kwargs = _adapter_kwargs(wavenumber=3)
     outer_reaction = 0.8
     inner_reaction = -0.4
-    stokes = PartitionedFakeRotatedStokes(
+    stokes = IntegralFakeRotatedStokes(
         {
-            "Upper": _circle_samples(
-                RADIUS_OUTER,
-                outer_reaction,
-                0.03,
-                kwargs["wavenumber"],
-            ),
-            "Lower": _circle_samples(
-                RADIUS_INNER,
-                inner_reaction,
-                -0.02,
-                kwargs["wavenumber"],
-            ),
+            "Upper": (RADIUS_OUTER, outer_reaction, 0.03),
+            "Lower": (RADIUS_INNER, inner_reaction, -0.02),
         }
     )
+    monkeypatch.setattr(uw.maths, "BdIntegral", FakeBoundaryIntegral)
     response = uw.postprocessing.geoid.cylindrical_annulus_response_from_rotated_stokes(
         stokes=stokes,
         include_self_gravity=True,
@@ -318,6 +418,14 @@ def test_partitioned_rotated_adapter_matches_pure_response():
 
     assert response.outer_reaction == pytest.approx(outer_reaction)
     assert response.inner_reaction == pytest.approx(inner_reaction)
+    assert response.outer_reaction_mean == pytest.approx(0.03)
+    assert response.inner_reaction_mean == pytest.approx(-0.02)
+    assert stokes.integral_calls == [
+        ("Upper", True),
+        ("Upper", False),
+        ("Lower", True),
+        ("Lower", False),
+    ]
     np.testing.assert_allclose(
         [response.outer_topography, response.inner_topography],
         [expected_outer_topography, expected_inner_topography],