diff --git a/docs/api/index.md b/docs/api/index.md index a0d9602bf..f7489d2cd 100644 --- a/docs/api/index.md +++ b/docs/api/index.md @@ -30,6 +30,7 @@ utilities visualisation adaptivity analytic +postprocessing ``` ## Quick Links @@ -54,6 +55,9 @@ analytic ### Validation - **{doc}`analytic`** - Exact solutions for benchmarking and convergence testing +### Post-processing +- **{doc}`postprocessing`** - Boundary-response, geoid, and self-gravity coefficients + ### Infrastructure - **{doc}`model`** - Model management and configuration - **{doc}`utilities`** - I/O, mesh import, and helper functions diff --git a/docs/api/postprocessing.md b/docs/api/postprocessing.md new file mode 100644 index 000000000..07512bcc0 --- /dev/null +++ b/docs/api/postprocessing.md @@ -0,0 +1,19 @@ +# Post-processing + +```{eval-rst} +.. automodule:: underworld3.postprocessing + :members: + :show-inheritance: +``` + +## Geoid and self-gravity responses + +The geoid module provides coefficient-only spherical-shell and +cylindrical-annulus gravity operators plus adapters for completed +rotated-free-slip Stokes solves. + +```{eval-rst} +.. automodule:: underworld3.postprocessing.geoid + :members: + :show-inheritance: +``` diff --git a/docs/developer/CHANGELOG.md b/docs/developer/CHANGELOG.md index 6bc812aa3..d1af2c6f8 100644 --- a/docs/developer/CHANGELOG.md +++ b/docs/developer/CHANGELOG.md @@ -382,6 +382,11 @@ component exactly — correct on curved, tilted, and deformed boundaries (#293). existing boundary traction onto an axisymmetric harmonic; the pure functions also accept coefficients recovered by other methods and an optional internal load. +- Rotated free slip exposes + `Stokes.boundary_normal_traction_integral(boundary, fn)` for a distributed + weak contraction of the assembled normal reaction. Cylindrical-annulus + Stokes responses use this fitted integral and its matching finite-element + boundary norm instead of gathering pointwise samples for angular quadrature. - `uw.analytic.Zhong2008` implements the Hager--O'Connell propagator-matrix oracle used for the Zhong et al. spherical-shell response benchmark. It supports piecewise-constant radial viscosity and reproduces every analytical diff --git a/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md b/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md index fb007f0ba..6aef79807 100644 --- a/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md +++ b/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md @@ -137,10 +137,105 @@ Published reference solvers, such as the Zhong et al. propagator-matrix method, belong in `uw.analytic`; their computed topography coefficients can be passed to the pure post-processing functions above. -The rotated harmonic projector gathers boundary samples to rank zero and -reconstructs their spherical triangulation. A future boundary-reaction -functional could replace this step with a direct distributed finite-element -projection without changing the coefficient API. +`Stokes.boundary_normal_traction_integral(boundary, fn)` contracts an assembled +normal-reaction load directly with a scalar test function over owned degrees of +freedom, followed by reductions on the mesh communicator. It is useful whenever +only an integrated or fitted traction diagnostic is required. It avoids +pointwise recovery and global boundary reconstruction; consumers that need a +nodal field should continue to use `boundary_normal_traction()` or +`dynamic_topography()` and follow their curved-P2 guidance. + +## 4. Cylindrical-annulus gravity and geoid response + +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, + ) +) +``` + +The Stokes adapter contracts the assembled reaction directly with +`cos(n theta)` through `boundary_normal_traction_integral()`. The coefficient is +normalised by the matching `BdIntegral` of `cos(n theta)**2`, so its numerator +and denominator use the same faceted finite-element boundary geometry. The +calculation counts owned reaction degrees of freedom and reduces on the mesh +communicator; it does not recover pointwise traction, sort samples by angle, or +gather a global boundary on rank zero. This construction remains valid on a +deformed boundary, where the fitted coefficient is weighted by the actual +finite-element boundary measure. + +`cylindrical_cosine_boundary_coefficient()` remains available for external +sampled data on a mathematical circle; it is not used for finite-element +reaction loads. The coefficient kernel follows Simons (1996), Appendix B. +Complete Kramer--Simons finite-element convergence and physical-space Poisson +comparisons remain in the separate mantle-convection benchmark repository. ## See also diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 2b7dc4e9c..d307a6677 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -6512,6 +6512,25 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): return _dtf(self, boundary, self._rotated_freeslip_info, field, buoyancy_scale=buoyancy_scale, mass=mass) + def boundary_normal_traction_integral(self, boundary, fn, remove_mean=True): + r"""Return the boundary integral of ``sigma_nn * fn`` directly from the + rotated-free-slip constraint reaction. + + With ``remove_mean=True`` (default), the constant normal-traction gauge + is removed before projection. Unlike pointwise + :meth:`boundary_normal_traction`, this weak projection does not recover + nodal traction values or gather a global boundary mesh. It is therefore + suitable for harmonic and integral diagnostics on curved P2 boundaries, + whose recovered vertex values converge slowly (issue #414). + """ + if self._rotated_freeslip_info is None: + raise RuntimeError( + "boundary_normal_traction_integral requires a completed " + "rotated-free-slip solve.") + from underworld3.utilities.rotated_bc import boundary_normal_traction_integral as _bnti + return _bnti(self, boundary, self._rotated_freeslip_info, fn, + remove_mean=remove_mean) + def add_nitsche_bc(self, conds=None, boundary=None, direction=None, normal=None, gamma=10.0, theta=1, mask=None, local_h=True, g=None): r"""Add Nitsche weak enforcement of a velocity constraint along a direction. diff --git a/src/underworld3/postprocessing/geoid.py b/src/underworld3/postprocessing/geoid.py index 537a79fee..3a2e18eb6 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,274 @@ 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 sampled circular data onto :math:`\cos(n\theta)` and the mean. + + ``coords`` must contain at least three two-dimensional Cartesian boundary + points. The samples may begin at any angle and need not include a duplicate + endpoint; this function sorts and closes the periodic interval. This helper + is intended for external sampled data. Finite-element reaction loads use + :meth:`Stokes.boundary_normal_traction_integral` instead, so their projection + follows the actual boundary facets without a global gather or reconstructed + angular ordering. + + Parameters + ---------- + 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 _rotated_cylindrical_boundary_response( + *, + stokes, + boundary: str, + wavenumber: int, + buoyancy_scale: float, +) -> tuple[float, float, float]: + """Project one boundary's assembled 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.") + + import sympy + from underworld3.maths import BdIntegral + + theta = stokes.mesh.CoordinateSystem.xR[1] + harmonic = sympy.cos(mode * theta) + reaction_integral = float( + stokes.boundary_normal_traction_integral( + boundary, + harmonic, + remove_mean=True, + ) + ) + reaction_total = float( + stokes.boundary_normal_traction_integral( + boundary, + 1.0, + remove_mean=False, + ) + ) + boundary_measure = float( + BdIntegral(stokes.mesh, fn=1.0, boundary=boundary).evaluate() + ) + harmonic_norm = float( + BdIntegral(stokes.mesh, fn=harmonic**2, boundary=boundary).evaluate() + ) + if not np.all( + np.isfinite( + (reaction_integral, reaction_total, boundary_measure, harmonic_norm) + ) + ): + raise RuntimeError("Cylindrical reaction projection produced non-finite data.") + if boundary_measure <= 0.0 or harmonic_norm <= 0.0: + raise RuntimeError( + "Cylindrical reaction projection requires positive boundary integrals." + ) + + reaction_coefficient = reaction_integral / harmonic_norm + reaction_mean = reaction_total / boundary_measure + return ( + float(reaction_coefficient), + float(reaction_mean), + float(-reaction_coefficient / buoyancy_scale), + ) + + +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_integral` contracts the assembled + wall reaction directly with a boundary test function. This adapter uses + + .. math:: + + h=-reaction_{nn}/signed\_buoyancy\_scale + + and projects each boundary onto the unnormalised real basis + :math:`\cos(n\theta)`. The numerator and harmonic norm use the same finite- + element boundary measure. Owned reaction degrees of freedom are reduced on + the mesh communicator; no pointwise traction recovery, angular sorting, or + rank-zero boundary gather is performed. + + Parameters + ---------- + 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/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 4b3ac812e..992d7b2b2 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -2198,6 +2198,89 @@ def boundary_normal_traction(solver, boundary, solve_result, mass="auto"): remove_mean=True, partial_reaction=False) +def boundary_normal_traction_integral(solver, boundary, solve_result, fn, + remove_mean=True): + r"""Return ``integral((sigma_nn - mean) * fn, boundary)`` directly from the + assembled rotated-constraint reaction. + + This is the weak/integral counterpart of :func:`boundary_normal_traction`. + It contracts the nodal reaction with ``fn`` at the velocity interpolation + nodes before any pointwise boundary-mass recovery. On curved P2 boundaries + this is a fitted quantity and therefore does not consume the slowly + converging recovered vertex values described in issue #414. + + Each assembled reaction degree of freedom is counted on its owning rank, + followed by an MPI sum on the mesh communicator. The operation does not + gather boundary topology or recovered values onto rank zero. + ``remove_mean=True`` removes the constant-traction gauge using boundary + integrals of ``fn`` and one. + """ + if not isinstance(remove_mean, (bool, np.bool_)): + raise TypeError("remove_mean must be True or False.") + + import underworld3 as uw + + fn = sympy.sympify(fn) + dm = solver.dm + comm = dm.comm.tompi4py() + dim = solver.mesh.dim + rc = solve_result["reaction"] + rstart, rend = rc.getOwnershipRange() + rcl = dm.getLocalVec() + dm.globalToLocal(rc, rcl) + + try: + rca = np.asarray(rcl.getArray()) + lsec = dm.getLocalSection() + l2g = dm.getLGMap() + csec = dm.getCoordinateSection() + cvec = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, dim) + v0, v1 = dm.getDepthStratum(0) + normal = dict(_boundary_spec(s) for s in solve_result["boundaries"]).get( + boundary + ) + nodes = _boundary_velocity_nodes(solver, boundary, normal=normal) + + owned_coords = [] + owned_reactions = [] + for q, nrm in nodes: + lo = lsec.getFieldOffset(q, _VELOCITY_FIELD) + global_row = int(l2g.apply([lo])[0]) + if not rstart <= global_row < rend: + continue + owned_coords.append(_point_coord(dm, dim, cvec, csec, v0, v1, q)) + # sigma_nn load = -n.r_c, matching boundary_normal_traction(). + owned_reactions.append(-float(np.dot(nrm, rca[lo:lo + dim]))) + finally: + dm.restoreLocalVec(rcl) + + if owned_coords: + coords = np.ascontiguousarray(owned_coords, dtype=float) + weights = np.asarray(uw.function.evaluate(fn, coords), dtype=float).reshape(-1) + if weights.size == 1 and len(owned_reactions) != 1: + weights = np.full(len(owned_reactions), float(weights[0])) + if weights.size != len(owned_reactions): + raise ValueError("fn must evaluate to one scalar per boundary node.") + local_weighted = float(np.dot(owned_reactions, weights)) + local_total = float(np.sum(owned_reactions)) + else: + local_weighted = 0.0 + local_total = 0.0 + + weighted = float(comm.allreduce(local_weighted)) + if not remove_mean: + return weighted + + total = float(comm.allreduce(local_total)) + area = float(uw.maths.BdIntegral(solver.mesh, fn=1.0, boundary=boundary).evaluate()) + if not np.isfinite(area) or area <= 0.0: + raise RuntimeError(f"Boundary {boundary!r} has non-positive area {area}.") + fn_integral = float( + uw.maths.BdIntegral(solver.mesh, fn=fn, boundary=boundary).evaluate() + ) + return weighted - (total / area) * fn_integral + + def dynamic_topography_field(solver, boundary, solve_result, field, buoyancy_scale=1.0, mass="auto"): """Populate a scalar MeshVariable ``field`` with the dynamic topography diff --git a/tests/test_1073_postprocessing_cylindrical_geoid.py b/tests/test_1073_postprocessing_cylindrical_geoid.py new file mode 100644 index 000000000..8643ae34c --- /dev/null +++ b/tests/test_1073_postprocessing_cylindrical_geoid.py @@ -0,0 +1,447 @@ +from types import SimpleNamespace + +import numpy as np +import pytest +import sympy +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 IntegralFakeRotatedStokes: + """Expose fitted reaction functionals while rejecting pointwise recovery.""" + + def __init__(self, boundary_data): + radius, theta = sympy.symbols("r theta", real=True) + self.mesh = SimpleNamespace( + dim=2, + CoordinateSystem=SimpleNamespace(xR=(radius, theta)), + boundary_data=boundary_data, + ) + self.boundary_data = boundary_data + self.integral_calls = [] + + def boundary_normal_traction(self, boundary, mass="auto"): + raise AssertionError("The finite-element adapter must not recover point values.") + + def boundary_normal_traction_integral(self, boundary, fn, remove_mean=True): + radius, coefficient, mean = self.boundary_data[boundary] + self.integral_calls.append((boundary, bool(remove_mean))) + if remove_mean: + return coefficient * np.pi * radius + assert float(sympy.sympify(fn)) == 1.0 + return mean * 2.0 * np.pi * radius + + +class FakeBoundaryIntegral: + """Exact circular measure used to isolate the adapter's projection route.""" + + def __init__(self, mesh, fn, boundary): + self.mesh = mesh + self.fn = sympy.sympify(fn) + self.boundary = boundary + + def evaluate(self): + radius = self.mesh.boundary_data[self.boundary][0] + if self.fn.is_number and float(self.fn) == 1.0: + return 2.0 * np.pi * radius + return np.pi * radius + + +def _adapter_kwargs(wavenumber=2): + 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) + assert hasattr(uw.systems.Stokes, "boundary_normal_traction_integral") + + +def test_reaction_integral_requires_boolean_mean_removal(): + from underworld3.utilities.rotated_bc import boundary_normal_traction_integral + + with pytest.raises(TypeError, match="remove_mean must be True or False"): + boundary_normal_traction_integral( + solver=None, + boundary="Upper", + solve_result=None, + fn=1.0, + remove_mean="yes", + ) + + +def test_distributed_reaction_integral_on_finite_element_annulus(): + """Exercise the complete adapter on an assembled rotated-free-slip reaction.""" + + radius_inner = 0.5 + radius_outer = 1.0 + mode = 4 + mesh = uw.meshing.Annulus( + radiusInner=radius_inner, + radiusOuter=radius_outer, + cellSize=0.15, + qdegree=3, + ) + x, y = mesh.X + radius = sympy.sqrt(x**2 + y**2) + theta = sympy.atan2(y, x) + velocity = uw.discretisation.MeshVariable( + "V_geoid_integral", mesh, mesh.dim, degree=2, continuous=True + ) + pressure = uw.discretisation.MeshVariable( + "P_geoid_integral", mesh, 1, degree=1, continuous=True + ) + stokes = uw.systems.Stokes( + mesh, + velocityField=velocity, + pressureField=pressure, + ) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + radial_force = ( + sympy.cos(mode * theta) + * (radius - radius_inner) + * (radius_outer - radius) + * 40.0 + ) + stokes.bodyforce = sympy.Matrix( + [[x * radial_force / radius, y * radial_force / radius]] + ) + normal = sympy.Matrix([[x / radius, y / radius]]) + stokes.add_rotated_freeslip_bc(0, "Lower", normal=normal) + stokes.add_rotated_freeslip_bc(0, "Upper", normal=normal) + stokes.petsc_use_pressure_nullspace = True + stokes.petsc_options["snes_type"] = "ksponly" + stokes.solve() + + response = uw.postprocessing.geoid.cylindrical_annulus_response_from_rotated_stokes( + stokes=stokes, + radius_inner=radius_inner, + radius_outer=radius_outer, + wavenumber=mode, + outer_density_contrast=1.0, + inner_density_contrast=1.0, + outer_reference_gravity=1.0, + inner_reference_gravity=1.0, + ) + + assert response.outer_reaction == pytest.approx(0.4044747, rel=2.0e-4) + assert response.inner_reaction == pytest.approx(0.5158471, rel=2.0e-4) + assert np.all( + np.isfinite( + [ + response.outer_reaction_mean, + response.inner_reaction_mean, + response.outer_geoid, + response.inner_geoid, + ] + ) + ) + + +@pytest.mark.parametrize("wavenumber", [1, 2, 3, 4, 8]) +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_rotated_adapter_uses_distributed_reaction_integral(monkeypatch): + kwargs = _adapter_kwargs(wavenumber=3) + outer_reaction = 0.8 + inner_reaction = -0.4 + stokes = IntegralFakeRotatedStokes( + { + "Upper": (RADIUS_OUTER, outer_reaction, 0.03), + "Lower": (RADIUS_INNER, inner_reaction, -0.02), + } + ) + monkeypatch.setattr(uw.maths, "BdIntegral", FakeBoundaryIntegral) + response = uw.postprocessing.geoid.cylindrical_annulus_response_from_rotated_stokes( + stokes=stokes, + include_self_gravity=True, + **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) + assert response.outer_reaction_mean == pytest.approx(0.03) + assert response.inner_reaction_mean == pytest.approx(-0.02) + assert stokes.integral_calls == [ + ("Upper", True), + ("Upper", False), + ("Lower", True), + ("Lower", False), + ] + np.testing.assert_allclose( + [response.outer_topography, response.inner_topography], + [expected_outer_topography, expected_inner_topography], + ) + 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