From 8a3099b418cff8256aed69d6b74b20b1b80bc69f Mon Sep 17 00:00:00 2001 From: Tyagi Date: Wed, 26 Aug 2026 08:43:13 +0530 Subject: [PATCH 1/2] 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 537a79fe..c3bb6c0b 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 00000000..9083d3f1 --- /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 Date: Wed, 26 Aug 2026 08:43:19 +0530 Subject: [PATCH 2/2] 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 a0d9602b..f7489d2c 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 00000000..07512bcc --- /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 fb007f0b..e65ce571 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),