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), 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