diff --git a/docs/api/analytic.md b/docs/api/analytic.md index 77698aad..e4ff784e 100644 --- a/docs/api/analytic.md +++ b/docs/api/analytic.md @@ -17,6 +17,18 @@ whichever one you use. :show-inheritance: ``` +## Zhong spherical-shell response oracle + +`Zhong2008` is a mesh-independent propagator-matrix oracle. It returns the +boundary response coefficients used by the Zhong et al. (2008) benchmark; it +does not implement the symbolic mesh-field contract above. + +```{eval-rst} +.. automodule:: underworld3.analytic.zhong2008 + :members: + :show-inheritance: +``` + ## See also - {doc}`solvers` — the solvers these solutions validate. diff --git a/docs/developer/CHANGELOG.md b/docs/developer/CHANGELOG.md index 0b6f1607..53c0126f 100644 --- a/docs/developer/CHANGELOG.md +++ b/docs/developer/CHANGELOG.md @@ -273,6 +273,16 @@ component exactly — correct on curved, tilted, and deformed boundaries (#293). seam sensitivity is unchanged (#404). - Recorded as the preferred free-slip BC in the project guidance (#300); conda PETSc floor raised to ≥ 3.25 for FMG/rotation API consistency (#304). +- `uw.postprocessing.geoid` provides generic spherical-shell geoid and + self-gravity coefficient functions. Its rotated-Stokes adapter projects the + existing boundary traction onto an axisymmetric harmonic; the pure functions + also accept coefficients recovered by other methods and an optional internal + load. +- `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 + response printed in Zhong Tables 2 and 3; geoid and self-gravity are delegated + to the generic postprocessing functions above. ### Generalized Geometric Multigrid via Custom Prolongation (July 2026) diff --git a/docs/developer/subsystems/analytic-solutions.md b/docs/developer/subsystems/analytic-solutions.md index 9f7f079b..71536e6c 100644 --- a/docs/developer/subsystems/analytic-solutions.md +++ b/docs/developer/subsystems/analytic-solutions.md @@ -818,6 +818,98 @@ Two class attributes carry this: | `symbolic` | fields are SymPy on `mesh.X`. False means the residual gates cannot be applied at all | | `requires` | name of an optional package, or `None` | +## Numeric response oracles + +Not every published analytical benchmark is a pointwise field on a mesh. The +Zhong et al. (2008) spherical-shell benchmark publishes scalar response +functions computed by a semi-analytic radial propagator. It belongs in +`uw.analytic`, but it does **not** subclass `AnalyticSolution` and is not listed +by `available()` because it cannot satisfy that class's symbolic mesh-field +contract. + +```python +reference = uw.analytic.Zhong2008( + harmonic_degree=2, + radius_inner=0.55, + radius_outer=1.0, + internal_load_radius=0.775, +) +response = reference.response() + +print(response.surface_characteristic_velocity) +print(response.self_gravity.surface_topography) +``` + +### Radial state and propagator + +For one spherical-harmonic degree `l`, define `L = l(l + 1)` and use +`v = log(r)`. The four-component poloidal state follows Hager and O'Connell +(1981): + +```text +u = (y1, y2, r sigma_rr / eta0, r sigma_r_perp / eta0)^T +``` + +`y1` is radial velocity and `y2` is the characteristic horizontal velocity. +Within one constant-viscosity layer, with dimensionless viscosity `eta`, + +```text +du/dv = A u + + [ -2 L 0 0 ] + [ -1 1 0 1/eta ] +A = [ 12 eta -6 L eta 1 L ] + [ -6 eta 2(2L-1) eta -1 -2 ] +``` + +The exact layer transfer is `exp(A log(r_b/r_a))`. Transfers are multiplied in +increasing-radius order. Because the state stores physical tractions relative +to one fixed reference viscosity, all four components remain continuous at a +viscosity interface. + +The outward radial delta load of coefficient `a` at `rint` gives the jump + +```text +u(rint+) - u(rint-) = (0, 0, -rint a, 0)^T. +``` + +Impermeable free slip requires `u[0] = u[3] = 0` at the CMB and surface. The +two remaining CMB state components are obtained from a dense two-by-two solve. + +### Recovering published quantities + +With CMB and surface states `ub` and `us`, respectively: + +```text +surface topography without self-gravity = -us[2] / radius_outer +CMB topography without self-gravity = ub[2] / radius_inner +surface characteristic velocity = us[1] +CMB characteristic velocity = ub[1] +surface horizontal divergence = -L us[1] / radius_outer +CMB horizontal divergence = -L ub[1] / radius_inner +``` + +The topography coefficients then pass through the generic +`uw.postprocessing.geoid.spherical_shell_geoid_response()` and +`spherical_shell_self_gravity_response()` functions. The propagator module does +not carry a second geoid implementation. + +`Zhong2008Response.surface_topography` and `.surface_geoid` are the +no-self-gravity quantities. The corresponding Table 2 or Table 3 quantities +are under `.self_gravity`. Velocity is unchanged by that postprocessing. + +### Validation + +`tests/test_1029_analytic_zhong2008.py` checks every parenthesized analytical +entry in Zhong et al. (2008) Tables 2 and 3: three load depths, four harmonic +degrees, isoviscous and `10^4`-lid cases, and eight response quantities per +case. The 192 comparisons agree within the precision printed in the paper. The +test also checks both free-slip boundary states, load linearity, no-self-gravity +recovery, and invalid input handling. + +This is stronger than copying table values into benchmark scripts: the values +are independently recomputed from the governing radial system for every case. + ## Adding a new solution ### The format to supply a new solution in diff --git a/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md b/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md index e23e2e2f..fb007f0b 100644 --- a/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md +++ b/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md @@ -74,6 +74,74 @@ Combine the two rules — project the τ components with `linear_solver()` for t cheap linear solve, then compose `σ_rr` analytically — for accurate boundary stress recovery that also scales. +## 3. Spherical-shell geoid and self-gravity response + +Rotated free slip already exposes normal traction through +`Stokes.boundary_normal_traction()`. The convenience adapter projects that +recovered traction onto the unnormalised axisymmetric `P_l^0` harmonic and +applies the spherical-shell geoid operator: + +```python +stokes.solve() + +response = uw.postprocessing.geoid.spherical_shell_response_from_rotated_stokes( + stokes=stokes, + radius_inner=0.55, + radius_outer=1.0, + harmonic_degree=2, + internal_load_radius=0.775, + internal_load_coefficient=1.0, + include_self_gravity=True, + surface_density_contrast=3300.0, + cmb_density_contrast=5400.0, + planet_radius=6370000.0, + gravity=9.8, + gravitational_constant=6.67e-11, +) +``` + +The adapter delegates stress recovery to the existing rotated-free-slip API; +it does not implement a second CBF, constrained-multiplier, or topography +recovery path. `internal_load_coefficient` must use the same harmonic +normalisation and sign convention as the model's internal load. + +When surface and CMB topography coefficients are already available, call +`uw.postprocessing.geoid.spherical_shell_geoid_response()` or +`uw.postprocessing.geoid.spherical_shell_self_gravity_response()` directly. +These functions are pure post-processing, work for any spherical-harmonic +order with a consistent coefficient normalisation, and do not require a Stokes +object. They support non-negative harmonic degrees, and the internal load is +optional. The rotated-Stokes adapter requires degree one or greater because +normal-traction recovery removes the degree-zero mean. + +The density contrasts, dimensional outer-radius scale, and gravity are required +when self-gravity is enabled. They deliberately have no Earth- or +benchmark-specific defaults. The universal gravitational constant defaults to +the current CODATA value and can be overridden when reproducing a paper's +rounded constant. + +The calculation is expressed as one linear operator: + +```text +N = G h + n_load +(I - Q G) h_self_gravity = h + Q n_load +``` + +where `h` contains surface and CMB topography, `N` contains their geoid +responses, and `Q` contains the two self-gravity density factors. Sharing `G` +and `n_load` between the no-self-gravity and self-gravity paths avoids separate +scalar implementations of the same coefficients. + +This module does not compute a benchmark's semi-analytical Stokes solution. +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. + ## See also - Issues [#156] (projection solver settings), [#157] (projection memory), diff --git a/src/underworld3/__init__.py b/src/underworld3/__init__.py index 818ff233..450e2368 100644 --- a/src/underworld3/__init__.py +++ b/src/underworld3/__init__.py @@ -201,6 +201,7 @@ def view(): import underworld3.maths import underworld3.swarm import underworld3.systems +import underworld3.postprocessing import underworld3.maths import underworld3.utilities import underworld3.model diff --git a/src/underworld3/analytic/__init__.py b/src/underworld3/analytic/__init__.py index 6d7690c9..95de55bd 100644 --- a/src/underworld3/analytic/__init__.py +++ b/src/underworld3/analytic/__init__.py @@ -37,7 +37,21 @@ from .kramer import CylindricalStokes from .richards import GardnerSteady, GardnerTransient from .transport import AdvectedFront, ErfcDiffusion, Poisson1D, TwoLayerDarcy -from .velic import SolA, SolB, SolC, SolCx, SolDA, SolDB2d, SolDB3d, SolH, SolKx, SolKz, SolM, SolNL +from .velic import ( + SolA, + SolB, + SolC, + SolCx, + SolDA, + SolDB2d, + SolDB3d, + SolH, + SolKx, + SolKz, + SolM, + SolNL, +) +from .zhong2008 import Zhong2008, Zhong2008Response __all__ = [ "AnalyticSolution", @@ -64,6 +78,8 @@ "SolM", "SolNL", "TwoLayerDarcy", + "Zhong2008", + "Zhong2008Response", "available", "describe", "is_available", diff --git a/src/underworld3/analytic/zhong2008.py b/src/underworld3/analytic/zhong2008.py new file mode 100644 index 00000000..e5e29d1c --- /dev/null +++ b/src/underworld3/analytic/zhong2008.py @@ -0,0 +1,346 @@ +r"""Semi-analytic spherical-shell Stokes responses from Zhong et al. (2008). + +The benchmark applies one radial delta-function load in a spherical harmonic +and asks for the surface and CMB response. Following Hager and O'Connell +(1981), the poloidal Stokes equations are written as a four-component first +order system in :math:`v=\ln r` and propagated exactly through each +constant-viscosity layer with a matrix exponential. + +This is a numeric reference oracle, not a finite-element solve and not a SymPy +field on a mesh. It returns the scalar response coefficients published in +Tables 2--4 of Zhong et al. (2008). The no-self-gravity topography is recovered +from radial stress first; geoid and self-gravity feedback are then delegated to +the generic functions in :mod:`underworld3.postprocessing.geoid`. + +References +---------- +Hager, B.H. & O'Connell, R.J. (1981). A simple global model of plate dynamics +and mantle convection. *Journal of Geophysical Research* 86, 4843--4878. +doi:10.1029/JB086iB06p04843 + +Zhong, S. et al. (2008). A benchmark study on mantle convection in a 3-D +spherical shell using CitcomS. *Geochemistry, Geophysics, Geosystems* 9, +Q10017. doi:10.1029/2008GC002048 +""" + +from __future__ import annotations + +from dataclasses import dataclass +from numbers import Integral + +import numpy as np +from scipy.linalg import expm + +from underworld3.postprocessing.geoid import ( + SelfGravityResponse, + spherical_shell_geoid_response, + spherical_shell_self_gravity_response, +) + + +__all__ = ["Zhong2008", "Zhong2008Response"] + + +@dataclass(frozen=True) +class Zhong2008Response: + r"""Boundary response coefficients for one harmonic delta load. + + ``surface_topography``, ``cmb_topography``, and the two corresponding geoid + coefficients exclude self-gravity. The values including self-gravity are + grouped in ``self_gravity``. This mirrors the separation in Zhong et al. + (2008): self-gravity changes the topography/geoid interpretation but not the + Stokes velocity. + """ + + surface_topography: float + cmb_topography: float + surface_geoid: float + cmb_geoid: float + surface_characteristic_velocity: float + cmb_characteristic_velocity: float + surface_velocity_divergence: float + cmb_velocity_divergence: float + self_gravity: SelfGravityResponse + + +class Zhong2008: + r"""Propagator-matrix reference solution for Zhong et al. (2008). + + Parameters + ---------- + harmonic_degree : int + Spherical-harmonic degree :math:`l`. The radial response is independent + of harmonic order :math:`m`. + radius_inner, radius_outer : float + Nondimensional CMB and surface radii. + internal_load_radius : float + Radius of the radial delta-function load. + internal_load_coefficient : float + Coefficient of the outward load. The Zhong benchmark uses one. + viscosity_interfaces : sequence of float + Strictly increasing radii at which viscosity changes. Omit for an + isoviscous shell. + viscosities : sequence of float + Positive dimensionless viscosities from the innermost to outermost + layer. Its length must be one greater than ``viscosity_interfaces``. + surface_density_contrast, cmb_density_contrast : float + Signed density contrasts used only for self-gravity postprocessing, in + kg/m^3. Defaults reproduce Zhong et al. (2008). + planet_radius : float + Dimensional planet radius in metres, used only for self-gravity. + gravity : float + Surface gravitational acceleration in m/s^2. + gravitational_constant : float + Gravitational constant in SI units. The default is the value used in + Zhong et al. (2008), not the current CODATA value. + + Notes + ----- + The propagated state is + + .. math:: + + \mathbf u=(y_1,y_2,r\sigma_{rr}/\eta_0,r\sigma_{r\perp}/\eta_0)^T, + + where :math:`y_1` is radial velocity and :math:`y_2` is the characteristic + horizontal velocity. Impermeable free slip is therefore ``u[0] = u[3] = + 0`` at both boundaries. Across the load at :math:`r=r_0`, only the radial + traction state jumps: + + .. math:: + + \mathbf u(r_0^+)-\mathbf u(r_0^-) + =(0,0,-r_0 A,0)^T. + + Examples + -------- + The isoviscous, degree-two, mid-mantle case from Table 2: + + >>> result = Zhong2008(harmonic_degree=2, internal_load_radius=0.775).response() + >>> round(result.self_gravity.surface_topography, 4) + 0.4998 + + The layered case from Table 3 has a :math:`10^4` viscosity lid: + + >>> result = Zhong2008( + ... harmonic_degree=2, + ... internal_load_radius=0.775, + ... viscosity_interfaces=(0.971875,), + ... viscosities=(1.0, 1.0e4), + ... ).response() + >>> round(result.self_gravity.surface_geoid, 5) + 0.05789 + """ + + reference = ( + "Zhong et al. (2008), Geochem. Geophys. Geosyst. 9, Q10017, " + "doi:10.1029/2008GC002048; propagator formulation after Hager and " + "O'Connell (1981), doi:10.1029/JB086iB06p04843." + ) + + def __init__( + self, + *, + harmonic_degree: int = 2, + radius_inner: float = 0.55, + radius_outer: float = 1.0, + internal_load_radius: float = 0.775, + internal_load_coefficient: float = 1.0, + viscosity_interfaces=(), + viscosities=(1.0,), + surface_density_contrast: float = 3300.0, + cmb_density_contrast: float = 5400.0, + planet_radius: float = 6_370_000.0, + gravity: float = 9.8, + gravitational_constant: float = 6.67e-11, + ): + if isinstance(harmonic_degree, bool) or not isinstance( + harmonic_degree, Integral + ): + raise TypeError("harmonic_degree must be an integer.") + if harmonic_degree < 1: + raise ValueError("harmonic_degree must be positive.") + + ri = float(radius_inner) + ro = float(radius_outer) + rint = float(internal_load_radius) + if not np.all(np.isfinite((ri, ro, rint))) or not 0.0 < ri < rint < ro: + raise ValueError( + "Expected finite radii ordered as 0 < radius_inner < " + "internal_load_radius < radius_outer." + ) + + interfaces = np.asarray(tuple(viscosity_interfaces), dtype=float) + layer_viscosities = np.asarray(tuple(viscosities), dtype=float) + if interfaces.ndim != 1 or layer_viscosities.ndim != 1: + raise ValueError( + "Viscosity interfaces and viscosities must be one-dimensional." + ) + if len(layer_viscosities) != len(interfaces) + 1: + raise ValueError( + "viscosities must contain one value more than viscosity_interfaces." + ) + if ( + not np.all(np.isfinite(interfaces)) + or np.any(interfaces <= ri) + or np.any(interfaces >= ro) + or np.any(np.diff(interfaces) <= 0.0) + ): + raise ValueError( + "viscosity_interfaces must be finite, strictly increasing, and " + "strictly inside the shell." + ) + if ( + not len(layer_viscosities) + or not np.all(np.isfinite(layer_viscosities)) + or np.any(layer_viscosities <= 0.0) + ): + raise ValueError("viscosities must be finite and positive.") + + load = float(internal_load_coefficient) + physical = np.asarray( + ( + surface_density_contrast, + cmb_density_contrast, + planet_radius, + gravity, + gravitational_constant, + ), + dtype=float, + ) + if not np.isfinite(load): + raise ValueError("internal_load_coefficient must be finite.") + if not np.all(np.isfinite(physical[:2])): + raise ValueError("Density contrasts must be finite.") + if not np.all(np.isfinite(physical[2:])) or np.any(physical[2:] <= 0.0): + raise ValueError("Physical constants must be finite and positive.") + + self.harmonic_degree = int(harmonic_degree) + self.radius_inner = ri + self.radius_outer = ro + self.internal_load_radius = rint + self.internal_load_coefficient = load + self.viscosity_interfaces = tuple(float(value) for value in interfaces) + self.viscosities = tuple(float(value) for value in layer_viscosities) + self.surface_density_contrast = float(surface_density_contrast) + self.cmb_density_contrast = float(cmb_density_contrast) + self.planet_radius = float(planet_radius) + self.gravity = float(gravity) + self.gravitational_constant = float(gravitational_constant) + + def _system_matrix(self, viscosity: float) -> np.ndarray: + r"""Hager--O'Connell four-state poloidal matrix for one layer.""" + + degree = self.harmonic_degree + angular_eigenvalue = degree * (degree + 1) + eta = float(viscosity) + return np.array( + [ + [-2.0, angular_eigenvalue, 0.0, 0.0], + [-1.0, 1.0, 0.0, 1.0 / eta], + [12.0 * eta, -6.0 * angular_eigenvalue * eta, 1.0, angular_eigenvalue], + [-6.0 * eta, 2.0 * (2.0 * angular_eigenvalue - 1.0) * eta, -1.0, -2.0], + ], + dtype=float, + ) + + def _propagate( + self, state: np.ndarray, radius_a: float, radius_b: float + ) -> np.ndarray: + """Propagate one state upward, applying the load at its exact radius.""" + + breakpoints = sorted( + { + self.radius_inner, + self.radius_outer, + self.internal_load_radius, + *self.viscosity_interfaces, + float(radius_a), + float(radius_b), + } + ) + points = [point for point in breakpoints if radius_a <= point <= radius_b] + propagated = np.asarray(state, dtype=float).copy() + + for lower, upper in zip(points[:-1], points[1:]): + midpoint = 0.5 * (lower + upper) + layer = int(np.searchsorted(self.viscosity_interfaces, midpoint)) + matrix = self._system_matrix(self.viscosities[layer]) + propagated = expm(matrix * np.log(upper / lower)) @ propagated + if upper == self.internal_load_radius: + propagated[2] -= ( + self.internal_load_radius * self.internal_load_coefficient + ) + + return propagated + + def _boundary_states(self) -> tuple[np.ndarray, np.ndarray]: + """Solve the two free-slip boundary constraints for the boundary states.""" + + basis = np.zeros((4, 2), dtype=float) + basis[1, 0] = 1.0 + basis[2, 1] = 1.0 + offset = self._propagate( + np.zeros(4, dtype=float), self.radius_inner, self.radius_outer + ) + transferred_basis = np.column_stack( + [ + self._propagate(basis[:, column], self.radius_inner, self.radius_outer) + - offset + for column in range(2) + ] + ) + unknowns = np.linalg.solve(transferred_basis[[0, 3], :], -offset[[0, 3]]) + cmb_state = basis @ unknowns + surface_state = self._propagate(cmb_state, self.radius_inner, self.radius_outer) + return cmb_state, surface_state + + def response(self) -> Zhong2008Response: + r"""Compute no-self-gravity and self-gravity response coefficients.""" + + cmb_state, surface_state = self._boundary_states() + degree = self.harmonic_degree + angular_eigenvalue = degree * (degree + 1) + + surface_topography = -surface_state[2] / self.radius_outer + cmb_topography = cmb_state[2] / self.radius_inner + surface_velocity = surface_state[1] + cmb_velocity = cmb_state[1] + surface_divergence = -angular_eigenvalue * surface_velocity / self.radius_outer + cmb_divergence = -angular_eigenvalue * cmb_velocity / self.radius_inner + + geoid = spherical_shell_geoid_response( + radius_inner=self.radius_inner, + radius_outer=self.radius_outer, + harmonic_degree=degree, + surface_topography_coefficient=surface_topography, + cmb_topography_coefficient=cmb_topography, + internal_load_radius=self.internal_load_radius, + internal_load_coefficient=self.internal_load_coefficient, + ) + self_gravity = spherical_shell_self_gravity_response( + radius_inner=self.radius_inner, + radius_outer=self.radius_outer, + harmonic_degree=degree, + surface_topography_coefficient=surface_topography, + cmb_topography_coefficient=cmb_topography, + internal_load_radius=self.internal_load_radius, + internal_load_coefficient=self.internal_load_coefficient, + surface_density_contrast=self.surface_density_contrast, + cmb_density_contrast=self.cmb_density_contrast, + planet_radius=self.planet_radius, + gravity=self.gravity, + gravitational_constant=self.gravitational_constant, + ) + + return Zhong2008Response( + surface_topography=float(surface_topography), + cmb_topography=float(cmb_topography), + surface_geoid=geoid.surface_geoid, + cmb_geoid=geoid.cmb_geoid, + surface_characteristic_velocity=float(surface_velocity), + cmb_characteristic_velocity=float(cmb_velocity), + surface_velocity_divergence=float(surface_divergence), + cmb_velocity_divergence=float(cmb_divergence), + self_gravity=self_gravity, + ) diff --git a/src/underworld3/postprocessing/__init__.py b/src/underworld3/postprocessing/__init__.py new file mode 100644 index 00000000..4e4325a1 --- /dev/null +++ b/src/underworld3/postprocessing/__init__.py @@ -0,0 +1,10 @@ +r""" +Post-processing helpers for Underworld3 models. + +The post-processing package contains derived diagnostics that are useful across +multiple model scripts but do not belong in a solver implementation. +""" + +from . import geoid + +__all__ = ["geoid"] diff --git a/src/underworld3/postprocessing/geoid.py b/src/underworld3/postprocessing/geoid.py new file mode 100644 index 00000000..537a79fe --- /dev/null +++ b/src/underworld3/postprocessing/geoid.py @@ -0,0 +1,430 @@ +r"""Spherical-harmonic 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. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from numbers import Integral + +from mpi4py import MPI +import numpy as np + + +__all__ = [ + "GeoidResponse", + "SelfGravityResponse", + "SphericalShellResponse", + "spherical_shell_geoid_response", + "spherical_shell_self_gravity_response", + "spherical_shell_response_from_rotated_stokes", +] + + +@dataclass(frozen=True) +class GeoidResponse: + """Surface and CMB geoid coefficients without self-gravity feedback.""" + + surface_geoid: float + cmb_geoid: float + + +@dataclass(frozen=True) +class SelfGravityResponse: + """Self-gravity-corrected topography and geoid coefficients.""" + + surface_topography: float + cmb_topography: float + surface_geoid: float + cmb_geoid: float + q_surface: float + q_cmb: float + + +@dataclass(frozen=True) +class SphericalShellResponse: + """Spherical-shell topography and corresponding geoid coefficients.""" + + surface_topography: float + cmb_topography: float + surface_geoid: float + cmb_geoid: float + self_gravity: SelfGravityResponse | None = None + + +def _validate_geometry( + radius_inner: float, + radius_outer: float, + harmonic_degree: int, +) -> tuple[float, float, int]: + try: + ri = float(radius_inner) + ro = float(radius_outer) + 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.") + if isinstance(harmonic_degree, bool) or not isinstance(harmonic_degree, Integral): + raise TypeError("harmonic_degree must be an integer.") + degree = int(harmonic_degree) + if degree < 0: + raise ValueError("harmonic_degree must be non-negative.") + return ri, ro, degree + + +def _spherical_shell_geoid_operator( + radius_inner: float, + radius_outer: float, + harmonic_degree: int, + internal_load_radius: float | None, + internal_load_coefficient: float, +) -> tuple[np.ndarray, np.ndarray]: + """Return the boundary-topography operator and fixed-load vector.""" + + ri, ro, degree = _validate_geometry( + radius_inner, + radius_outer, + harmonic_degree, + ) + internal_load_coefficient = float(internal_load_coefficient) + if not np.isfinite(internal_load_coefficient): + raise ValueError("internal_load_coefficient must be finite.") + + rint = None + if internal_load_radius is not None: + 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 + if not np.isfinite(rint) or not ri < rint < ro: + 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." + ) + + denominator = float(2 * degree + 1) + operator = np.array( + [ + [ + ro / denominator, + ri * (ri / ro) ** (degree + 1) / denominator, + ], + [ + ro * (ri / ro) ** degree / denominator, + ri / denominator, + ], + ], + dtype=float, + ) + load = np.zeros(2, dtype=float) + if rint is not None and internal_load_coefficient != 0.0: + load = -internal_load_coefficient * np.array( + [ + rint * (rint / ro) ** (degree + 1) / denominator, + rint * (ri / rint) ** degree / denominator, + ], + dtype=float, + ) + return operator, load + + +def spherical_shell_geoid_response( + *, + radius_inner: float, + radius_outer: float, + harmonic_degree: int, + surface_topography_coefficient: float, + cmb_topography_coefficient: float, + internal_load_radius: float | None = None, + internal_load_coefficient: float = 0.0, +) -> GeoidResponse: + r"""Return spherical-shell geoid coefficients without self-gravity. + + All source coefficients must use one consistent spherical-harmonic + normalisation. The radial potential kernel depends on degree but not order, + so coefficients for any order :math:`m` may be used. A positive optional + internal-load coefficient follows the outward radial-load convention and + enters the potential with the opposite sign to boundary topography. + """ + + operator, load = _spherical_shell_geoid_operator( + radius_inner, + radius_outer, + harmonic_degree, + internal_load_radius, + internal_load_coefficient, + ) + topography = np.array( + [surface_topography_coefficient, cmb_topography_coefficient], + dtype=float, + ) + if not np.all(np.isfinite(topography)): + raise ValueError("The topography coefficients must be finite.") + geoid = operator @ topography + load + return GeoidResponse(float(geoid[0]), float(geoid[1])) + + +def spherical_shell_self_gravity_response( + *, + radius_inner: float, + radius_outer: float, + harmonic_degree: int, + surface_topography_coefficient: float, + cmb_topography_coefficient: float, + surface_density_contrast: float, + cmb_density_contrast: float, + planet_radius: float, + gravity: float, + internal_load_radius: float | None = None, + internal_load_coefficient: float = 0.0, + gravitational_constant: float = 6.67430e-11, +) -> SelfGravityResponse: + r"""Return self-gravity-corrected spherical-shell response coefficients. + + Density contrasts are signed and in kg/m³, ``planet_radius`` is in metres, + ``gravity`` is in m/s², and ``gravitational_constant`` is in SI units. + Shell radii and response coefficients remain nondimensional. + """ + + operator, load = _spherical_shell_geoid_operator( + radius_inner, + radius_outer, + harmonic_degree, + internal_load_radius, + internal_load_coefficient, + ) + topography = np.array( + [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) + if not np.all(np.isfinite(topography)): + raise ValueError("The topography coefficients must be finite.") + if not np.all(np.isfinite(density_contrasts)): + raise ValueError("Density contrasts must be finite.") + if not np.all(np.isfinite(physical_constants)) or np.any(physical_constants <= 0.0): + raise ValueError("Physical constants must be finite and positive.") + + gravity_scale = 4.0 * np.pi * gravitational_constant * planet_radius / gravity + q = gravity_scale * density_contrasts + q_matrix = np.diag(q) + corrected_topography = np.linalg.solve( + np.eye(2) - q_matrix @ operator, + topography + q_matrix @ load, + ) + corrected_geoid = operator @ corrected_topography + load + + return SelfGravityResponse( + surface_topography=float(corrected_topography[0]), + cmb_topography=float(corrected_topography[1]), + surface_geoid=float(corrected_geoid[0]), + cmb_geoid=float(corrected_geoid[1]), + q_surface=float(q[0]), + q_cmb=float(q[1]), + ) + + +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)) + return 2.0 * np.arctan2(determinant, denominator) * radius**2 + + +def _project_spherical_harmonic_samples( + coords: np.ndarray, + values: np.ndarray, + radius: float, + harmonic_degree: int, + response_sign: float, +) -> float: + from scipy.spatial import ConvexHull + + unit_coords = coords / np.linalg.norm(coords, axis=1)[:, None] + hull = ConvexHull(unit_coords, qhull_options="QJ") + projected_integral = 0.0 + harmonic = np.polynomial.legendre.Legendre.basis(harmonic_degree) + + for simplex in hull.simplices: + area = _spherical_triangle_area( + unit_coords[simplex[0]], + unit_coords[simplex[1]], + unit_coords[simplex[2]], + radius, + ) + centroid = unit_coords[simplex].mean(axis=0) + centroid /= np.linalg.norm(centroid) + projected_integral += ( + area * float(values[simplex].mean()) * harmonic(np.clip(centroid[2], -1, 1)) + ) + + harmonic_norm = 4.0 * np.pi / (2 * harmonic_degree + 1) + return float(response_sign * projected_integral / (radius**2 * harmonic_norm)) + + +def _rotated_topography_coefficient( + stokes, + boundary: str, + radius: float, + harmonic_degree: int, + buoyancy_scale: float, + response_sign: float, +) -> float: + buoyancy_scale = float(buoyancy_scale) + if not np.isfinite(buoyancy_scale) or buoyancy_scale == 0.0: + raise ValueError("Boundary buoyancy scales must be finite and nonzero.") + + coords, sigma_nn = stokes.boundary_normal_traction(boundary, mass="auto") + local_rows = np.column_stack( + ( + np.asarray(coords, dtype=float), + -np.asarray(sigma_nn, dtype=float) / buoyancy_scale, + ) + ) + gathered_rows = MPI.COMM_WORLD.gather(local_rows, root=0) + + coefficient = None + root_error = None + if MPI.COMM_WORLD.rank == 0: + try: + 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}.") + merged = {} + for row in np.vstack(nonempty_rows): + key = tuple(np.round(row[:3], 12)) + if key not in merged: + merged[key] = [row[:3].copy(), 0.0, 0] + merged[key][1] += float(row[3]) + merged[key][2] += 1 + global_coords = np.array([item[0] for item in merged.values()]) + global_values = np.array([item[1] / item[2] for item in merged.values()]) + coefficient = _project_spherical_harmonic_samples( + global_coords, + global_values, + radius, + harmonic_degree, + response_sign, + ) + 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(f"Spherical-harmonic projection failed: {root_error}") + return float(MPI.COMM_WORLD.bcast(coefficient, root=0)) + + +def spherical_shell_response_from_rotated_stokes( + *, + stokes, + radius_inner: float, + radius_outer: float, + harmonic_degree: int, + internal_load_radius: float | None = None, + internal_load_coefficient: float = 0.0, + surface_boundary: str = "Upper", + cmb_boundary: str = "Lower", + surface_buoyancy_scale: float = 1.0, + cmb_buoyancy_scale: float = 1.0, + include_self_gravity: bool = False, + surface_density_contrast: float | None = None, + cmb_density_contrast: float | None = None, + planet_radius: float | None = None, + gravity: float | None = None, + gravitational_constant: float = 6.67430e-11, +) -> SphericalShellResponse: + r"""Compute spherical-shell response from a rotated-free-slip Stokes solve. + + Normal traction recovery is delegated to the existing + :meth:`Stokes.boundary_normal_traction` implementation. This adapter only + projects the two boundary responses onto the unnormalised + axisymmetric :math:`P_l^0` harmonic. Use the pure coefficient functions + directly for other harmonic orders or topography-recovery methods. Density + contrasts, planet radius, and gravity are required when + ``include_self_gravity`` is true. + """ + + ri, ro, degree = _validate_geometry( + radius_inner, + radius_outer, + harmonic_degree, + ) + if not isinstance(include_self_gravity, bool): + raise TypeError("include_self_gravity must be True or False.") + if degree == 0: + raise ValueError( + "The rotated-Stokes adapter requires harmonic_degree >= 1 because " + "boundary_normal_traction() removes the degree-zero mean." + ) + if include_self_gravity: + missing = [ + name + for name, value in ( + ("surface_density_contrast", surface_density_contrast), + ("cmb_density_contrast", cmb_density_contrast), + ("planet_radius", planet_radius), + ("gravity", gravity), + ) + if value is None + ] + if missing: + raise ValueError( + "Self-gravity requires explicit values for " + ", ".join(missing) + "." + ) + + surface_topography = _rotated_topography_coefficient( + stokes, + surface_boundary, + ro, + degree, + surface_buoyancy_scale, + 1.0, + ) + cmb_topography = _rotated_topography_coefficient( + stokes, + cmb_boundary, + ri, + degree, + cmb_buoyancy_scale, + -1.0, + ) + geoid = spherical_shell_geoid_response( + radius_inner=ri, + radius_outer=ro, + harmonic_degree=degree, + surface_topography_coefficient=surface_topography, + cmb_topography_coefficient=cmb_topography, + internal_load_radius=internal_load_radius, + internal_load_coefficient=internal_load_coefficient, + ) + self_gravity = None + if include_self_gravity: + self_gravity = spherical_shell_self_gravity_response( + radius_inner=ri, + radius_outer=ro, + harmonic_degree=degree, + surface_topography_coefficient=surface_topography, + cmb_topography_coefficient=cmb_topography, + internal_load_radius=internal_load_radius, + internal_load_coefficient=internal_load_coefficient, + surface_density_contrast=float(surface_density_contrast), + cmb_density_contrast=float(cmb_density_contrast), + planet_radius=float(planet_radius), + gravity=float(gravity), + gravitational_constant=gravitational_constant, + ) + + return SphericalShellResponse( + surface_topography=surface_topography, + cmb_topography=cmb_topography, + surface_geoid=geoid.surface_geoid, + cmb_geoid=geoid.cmb_geoid, + self_gravity=self_gravity, + ) diff --git a/tests/parallel/test_1071_spherical_shell_geoid_parallel.py b/tests/parallel/test_1071_spherical_shell_geoid_parallel.py new file mode 100644 index 00000000..189d6121 --- /dev/null +++ b/tests/parallel/test_1071_spherical_shell_geoid_parallel.py @@ -0,0 +1,79 @@ +import numpy as np +import pytest +import sympy +import underworld3 as uw + + +pytestmark = pytest.mark.level_2 + + +def test_rotated_spherical_shell_geoid_matches_serial_reference(): + if uw.mpi.size == 1: + pytest.skip("Run this regression with at least two MPI ranks.") + + radius_inner = 0.55 + radius_outer = 1.0 + rint = 0.775 + mesh = uw.meshing.SphericalShellInternalBoundary( + radiusOuter=radius_outer, + radiusInternal=rint, + radiusInner=radius_inner, + cellSize=0.25, + qdegree=2, + degree=1, + ) + velocity = uw.discretisation.MeshVariable("U_geoid_mpi", mesh, mesh.dim, degree=2) + pressure = uw.discretisation.MeshVariable("P_geoid_mpi", 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 + + theta = mesh.CoordinateSystem.xR[1] + unit_r = mesh.CoordinateSystem.unit_e_0 + load = sympy.assoc_legendre(2, 0, sympy.cos(theta)) * unit_r + stokes.add_natural_bc(load, "Internal") + stokes.add_rotated_freeslip_bc(0, "Upper", normal=unit_r) + stokes.add_rotated_freeslip_bc(0, "Lower", normal=-unit_r) + stokes.petsc_use_nullspace = True + stokes.petsc_options["snes_type"] = "ksponly" + stokes.tolerance = 1.0e-5 + stokes.solve() + + response = uw.postprocessing.geoid.spherical_shell_response_from_rotated_stokes( + stokes=stokes, + radius_inner=radius_inner, + radius_outer=radius_outer, + harmonic_degree=2, + internal_load_radius=rint, + internal_load_coefficient=1.0, + include_self_gravity=True, + surface_density_contrast=3300.0, + cmb_density_contrast=5400.0, + planet_radius=6370000.0, + gravity=9.8, + gravitational_constant=6.67e-11, + ) + values = np.array( + [ + response.surface_topography, + response.cmb_topography, + response.surface_geoid, + response.cmb_geoid, + response.self_gravity.surface_topography, + response.self_gravity.cmb_topography, + response.self_gravity.surface_geoid, + response.self_gravity.cmb_geoid, + ] + ) + + for rank_values in uw.mpi.comm.allgather(values): + assert np.array_equal(rank_values, values) + + zhong_table_2 = np.array( + [0.41920, 0.77060, 0.02579, 0.03206, 0.49980, 0.93130, 0.04486, 0.05461] + ) + assert np.allclose(values, zhong_table_2, rtol=0.10, atol=1.0e-12) diff --git a/tests/test_1029_analytic_zhong2008.py b/tests/test_1029_analytic_zhong2008.py new file mode 100644 index 00000000..859f2b8a --- /dev/null +++ b/tests/test_1029_analytic_zhong2008.py @@ -0,0 +1,454 @@ +r"""Propagator-matrix responses from Zhong et al. (2008). + +The values in parentheses in Tables 2 and 3 are the semi-analytic propagator +solutions, not the CitcomS finite-element results. This test pins every +published propagator row so a sign, stress scaling, layer-order, or self-gravity +regression cannot hide behind agreement with one selected case. + +Run: pixi run python -m pytest tests/test_1029_analytic_zhong2008.py -v +""" + +import numpy as np +import pytest +import underworld3 as uw + + +pytestmark = pytest.mark.level_2 + + +# depth fraction, degree, surface topography, CMB topography, surface geoid, +# CMB geoid, surface characteristic velocity, CMB characteristic velocity, +# surface velocity divergence, CMB velocity divergence. These are the +# parenthesized propagator values in Zhong et al. (2008), Tables 2 and 3. +ISOVISCOUS_TABLE_2 = [ + ( + 0.25, + 2, + 0.7716, + 0.5646, + 0.04058, + 0.04062, + -9.949e-3, + 7.709e-3, + 5.969e-2, + -8.410e-2, + ), + ( + 0.25, + 5, + 0.7562, + 0.3894, + 0.02986, + 0.01556, + -4.546e-3, + 2.134e-3, + 1.364e-1, + -1.164e-1, + ), + ( + 0.25, + 8, + 0.6458, + 0.1629, + 0.02018, + 0.004454, + -2.128e-3, + 4.480e-4, + 1.532e-1, + -5.865e-2, + ), + ( + 0.25, + 15, + 0.3905, + 0.01094, + 0.008355, + 0.0001739, + -5.074e-4, + 1.044e-5, + 1.218e-1, + -4.557e-3, + ), + ( + 0.50, + 2, + 0.4998, + 0.9313, + 0.04486, + 0.05461, + -1.006e-2, + 1.186e-2, + 6.038e-2, + -1.294e-1, + ), + ( + 0.50, + 5, + 0.4238, + 0.7235, + 0.02427, + 0.02543, + -3.593e-3, + 3.733e-3, + 1.078e-1, + -2.036e-1, + ), + ( + 0.50, + 8, + 0.2679, + 0.3794, + 0.01122, + 0.009474, + -1.170e-3, + 9.890e-4, + 8.427e-2, + -1.295e-1, + ), + ( + 0.50, + 15, + 0.06905, + 0.05554, + 0.001804, + 0.0008399, + -1.091e-4, + 5.095e-5, + 2.618e-2, + -2.223e-2, + ), + ( + 0.75, + 2, + 0.2323, + 1.090, + 0.02788, + 0.04267, + -5.486e-3, + 1.046e-2, + 3.291e-2, + -1.141e-1, + ), + ( + 0.75, + 5, + 0.1664, + 1.008, + 0.01143, + 0.02741, + -1.600e-3, + 4.201e-3, + 4.801e-2, + -2.292e-1, + ), + ( + 0.75, + 8, + 0.07634, + 0.7471, + 0.003644, + 0.01542, + -3.707e-4, + 1.637e-3, + 2.669e-2, + -2.143e-1, + ), + ( + 0.75, + 15, + 0.007289, + 0.3006, + 0.0002061, + 0.004023, + -1.238e-5, + 2.453e-4, + 2.973e-3, + -1.071e-1, + ), +] + +LAYERED_TABLE_3 = [ + ( + 0.25, + 2, + 0.8798, + 0.07981, + 0.05333, + -0.006165, + -1.231e-5, + 1.543e-3, + 7.384e-5, + -1.683e-2, + ), + ( + 0.25, + 5, + 0.7976, + 0.1254, + 0.03326, + 0.002543, + -1.964e-6, + 7.093e-4, + 5.893e-5, + -3.869e-2, + ), + ( + 0.25, + 8, + 0.6913, + 0.07955, + 0.02284, + 0.001778, + -6.509e-7, + 2.191e-4, + 4.687e-5, + -2.869e-2, + ), + ( + 0.25, + 15, + 0.4501, + 0.008396, + 0.01028, + 0.0001290, + -1.103e-7, + 7.999e-6, + 2.647e-5, + -3.490e-3, + ), + ( + 0.50, + 2, + 0.6104, + 0.4354, + 0.05789, + 0.006754, + -1.258e-5, + 5.556e-3, + 7.549e-5, + -6.061e-2, + ), + ( + 0.50, + 5, + 0.4567, + 0.5110, + 0.02696, + 0.01496, + -1.577e-6, + 2.587e-3, + 4.730e-5, + -1.411e-1, + ), + ( + 0.50, + 8, + 0.2928, + 0.3323, + 0.01268, + 0.007963, + -3.654e-7, + 8.598e-4, + 2.631e-5, + -1.126e-1, + ), + ( + 0.50, + 15, + 0.08166, + 0.05497, + 0.002211, + 0.0008298, + -2.436e-8, + 5.040e-5, + 5.846e-6, + -2.199e-2, + ), + ( + 0.75, + 2, + 0.2927, + 0.8192, + 0.03500, + 0.01650, + -6.879e-6, + 7.009e-3, + 4.127e-5, + -7.646e-2, + ), + ( + 0.75, + 5, + 0.1811, + 0.9128, + 0.01263, + 0.02272, + -7.054e-7, + 3.688e-3, + 2.116e-5, + -2.012e-1, + ), + ( + 0.75, + 8, + 0.08421, + 0.7321, + 0.004105, + 0.01493, + -1.164e-7, + 1.596e-3, + 8.384e-6, + -2.089e-1, + ), + ( + 0.75, + 15, + 0.008713, + 0.3006, + 0.0002520, + 0.004022, + -2.789e-9, + 2.453e-4, + 6.693e-7, + -1.070e-1, + ), +] + + +def _computed_row(depth_fraction, harmonic_degree, *, layered): + radius_inner = 0.55 + radius_outer = 1.0 + rint = radius_outer - depth_fraction * (radius_outer - radius_inner) + options = {} + if layered: + # Four of the benchmark's 64 radial elements form the high-viscosity lid. + lid_base = radius_outer - 4.0 * (radius_outer - radius_inner) / 64.0 + options = { + "viscosity_interfaces": (lid_base,), + "viscosities": (1.0, 1.0e4), + } + + response = uw.analytic.Zhong2008( + harmonic_degree=harmonic_degree, + radius_inner=radius_inner, + radius_outer=radius_outer, + internal_load_radius=rint, + **options, + ).response() + return np.array( + [ + response.self_gravity.surface_topography, + response.self_gravity.cmb_topography, + response.self_gravity.surface_geoid, + response.self_gravity.cmb_geoid, + response.surface_characteristic_velocity, + response.cmb_characteristic_velocity, + response.surface_velocity_divergence, + response.cmb_velocity_divergence, + ] + ) + + +@pytest.mark.parametrize("published", ISOVISCOUS_TABLE_2) +def test_isoviscous_responses_match_every_analytic_value_in_table_2(published): + depth_fraction, harmonic_degree, *expected = published + computed = _computed_row(depth_fraction, harmonic_degree, layered=False) + + np.testing.assert_allclose(computed, expected, rtol=5.0e-4, atol=1.0e-12) + + +@pytest.mark.parametrize("published", LAYERED_TABLE_3) +def test_layered_responses_match_every_analytic_value_in_table_3(published): + depth_fraction, harmonic_degree, *expected = published + computed = _computed_row(depth_fraction, harmonic_degree, layered=True) + + np.testing.assert_allclose(computed, expected, rtol=5.0e-4, atol=1.0e-12) + + +def test_no_self_gravity_outputs_match_the_mid_mantle_reference(): + response = uw.analytic.Zhong2008().response() + + assert np.isclose(response.surface_topography, 0.4191904157) + assert np.isclose(response.cmb_topography, 0.7705825500) + assert np.isclose(response.surface_geoid, 0.0257906289) + assert np.isclose(response.cmb_geoid, 0.0320605845) + + +def test_boundary_states_satisfy_impermeable_free_slip(): + solution = uw.analytic.Zhong2008( + harmonic_degree=8, + internal_load_radius=0.6625, + viscosity_interfaces=(0.8, 0.95), + viscosities=(1.0, 30.0, 0.2), + ) + cmb_state, surface_state = solution._boundary_states() + + np.testing.assert_allclose(cmb_state[[0, 3]], 0.0, atol=1.0e-12) + np.testing.assert_allclose(surface_state[[0, 3]], 0.0, atol=1.0e-10) + + +def test_response_is_linear_in_load_amplitude(): + unit = uw.analytic.Zhong2008(internal_load_coefficient=1.0).response() + doubled = uw.analytic.Zhong2008(internal_load_coefficient=2.0).response() + + unit_values = np.array( + [ + unit.surface_topography, + unit.cmb_topography, + unit.surface_geoid, + unit.cmb_geoid, + unit.surface_characteristic_velocity, + unit.cmb_characteristic_velocity, + unit.self_gravity.surface_topography, + unit.self_gravity.cmb_topography, + unit.self_gravity.surface_geoid, + unit.self_gravity.cmb_geoid, + ] + ) + doubled_values = np.array( + [ + doubled.surface_topography, + doubled.cmb_topography, + doubled.surface_geoid, + doubled.cmb_geoid, + doubled.surface_characteristic_velocity, + doubled.cmb_characteristic_velocity, + doubled.self_gravity.surface_topography, + doubled.self_gravity.cmb_topography, + doubled.self_gravity.surface_geoid, + doubled.self_gravity.cmb_geoid, + ] + ) + + np.testing.assert_allclose(doubled_values, 2.0 * unit_values, rtol=1.0e-12) + + +@pytest.mark.parametrize( + "options, error, message", + [ + ({"harmonic_degree": 2.0}, TypeError, "must be an integer"), + ({"harmonic_degree": 0}, ValueError, "must be positive"), + ({"internal_load_radius": 1.1}, ValueError, "ordered as"), + ( + {"viscosity_interfaces": (0.8,), "viscosities": (1.0,)}, + ValueError, + "one value more", + ), + ( + {"viscosity_interfaces": (0.9, 0.8), "viscosities": (1.0, 2.0, 3.0)}, + ValueError, + "strictly increasing", + ), + ({"viscosities": (-1.0,)}, ValueError, "finite and positive"), + ], +) +def test_invalid_parameters_are_rejected(options, error, message): + with pytest.raises(error, match=message): + uw.analytic.Zhong2008(**options) + + +def test_zhong_oracle_is_public_but_not_a_symbolic_mesh_solution(): + assert "Zhong2008" in uw.analytic.__all__ + assert "Zhong2008Response" in uw.analytic.__all__ + assert "Zhong2008" not in uw.analytic.available() + assert "Zhong et al. (2008)" in uw.analytic.Zhong2008.reference diff --git a/tests/test_1070_postprocessing_geoid.py b/tests/test_1070_postprocessing_geoid.py new file mode 100644 index 00000000..dce646a8 --- /dev/null +++ b/tests/test_1070_postprocessing_geoid.py @@ -0,0 +1,240 @@ +import math + +import numpy as np +import pytest +import sympy +import underworld3 as uw + + +pytestmark = pytest.mark.level_2 + + +def test_spherical_shell_postprocessing_is_public(): + assert hasattr(uw, "postprocessing") + assert uw.postprocessing.__all__ == ["geoid"] + assert hasattr(uw.postprocessing.geoid, "spherical_shell_geoid_response") + assert hasattr(uw.postprocessing.geoid, "spherical_shell_self_gravity_response") + assert hasattr(uw.postprocessing.geoid, "spherical_shell_response_from_rotated_stokes") + + +def test_spherical_shell_geoid_response_matches_direct_formula_and_load(): + radius_inner = 0.55 + radius_outer = 1.0 + rint = 0.775 + degree = 2 + internal_load_coefficient = 1.75 + surface_topography = 0.3918264884592169 + cmb_topography = 0.2653834107104151 + + response = uw.postprocessing.geoid.spherical_shell_geoid_response( + radius_inner=radius_inner, + radius_outer=radius_outer, + harmonic_degree=degree, + surface_topography_coefficient=surface_topography, + cmb_topography_coefficient=cmb_topography, + internal_load_radius=rint, + internal_load_coefficient=internal_load_coefficient, + ) + + denominator = 2 * degree + 1 + surface_expected = ( + radius_inner * (radius_inner / radius_outer) ** (degree + 1) * cmb_topography + + radius_outer * surface_topography + - internal_load_coefficient * rint * (rint / radius_outer) ** (degree + 1) + ) / denominator + cmb_expected = ( + radius_inner * cmb_topography + + radius_outer * (radius_inner / radius_outer) ** degree * surface_topography + - internal_load_coefficient * rint * (radius_inner / rint) ** degree + ) / denominator + + assert math.isclose(response.surface_geoid, surface_expected) + assert math.isclose(response.cmb_geoid, cmb_expected) + + +def test_spherical_shell_geoid_response_without_internal_load(): + response = uw.postprocessing.geoid.spherical_shell_geoid_response( + radius_inner=0.55, + radius_outer=1.0, + harmonic_degree=3, + surface_topography_coefficient=0.4, + cmb_topography_coefficient=-0.2, + ) + + denominator = 7.0 + expected = np.array( + [ + (0.4 + 0.55 * 0.55**4 * -0.2) / denominator, + (0.55**3 * 0.4 + 0.55 * -0.2) / denominator, + ] + ) + assert np.allclose([response.surface_geoid, response.cmb_geoid], expected) + + +def test_spherical_shell_self_gravity_response_matches_direct_solve(): + surface_topography = 0.3918264884592169 + cmb_topography = 0.2653834107104151 + response = uw.postprocessing.geoid.spherical_shell_self_gravity_response( + radius_inner=0.55, + radius_outer=1.0, + harmonic_degree=2, + surface_topography_coefficient=surface_topography, + cmb_topography_coefficient=cmb_topography, + internal_load_radius=0.775, + internal_load_coefficient=1.0, + surface_density_contrast=3300.0, + cmb_density_contrast=5400.0, + planet_radius=6370000.0, + gravity=9.8, + gravitational_constant=6.67e-11, + ) + + q_surface = 4.0 * np.pi * 6.67e-11 * 6370000.0 * 3300.0 / 9.8 + q_cmb = 4.0 * np.pi * 6.67e-11 * 6370000.0 * 5400.0 / 9.8 + denominator = 5.0 + operator = np.array( + [ + [1.0 / denominator, 0.55 * 0.55**3 / denominator], + [0.55**2 / denominator, 0.55 / denominator], + ] + ) + load = np.array( + [ + -0.775 * 0.775**3 / denominator, + -0.775 * (0.55 / 0.775) ** 2 / denominator, + ] + ) + q_matrix = np.diag([q_surface, q_cmb]) + expected_topography = np.linalg.solve( + np.eye(2) - q_matrix @ operator, + np.array([surface_topography, cmb_topography]) + q_matrix @ load, + ) + expected_geoid = operator @ expected_topography + load + + assert math.isclose(response.q_surface, q_surface) + assert math.isclose(response.q_cmb, q_cmb) + assert np.allclose( + [response.surface_topography, response.cmb_topography], + expected_topography, + ) + assert np.allclose( + [response.surface_geoid, response.cmb_geoid], + expected_geoid, + ) + + +def test_spherical_shell_geoid_supports_degree_zero(): + response = uw.postprocessing.geoid.spherical_shell_geoid_response( + radius_inner=0.55, + radius_outer=1.0, + harmonic_degree=0, + surface_topography_coefficient=0.4, + cmb_topography_coefficient=-0.2, + ) + + assert math.isclose(response.surface_geoid, 0.4 - 0.55**2 * 0.2) + assert math.isclose(response.cmb_geoid, 0.4 - 0.55 * 0.2) + + +@pytest.mark.parametrize("harmonic_degree", [-1, 2.0, True]) +def test_spherical_shell_geoid_rejects_invalid_harmonic_degree(harmonic_degree): + error = TypeError if isinstance(harmonic_degree, (float, bool)) else ValueError + with pytest.raises(error): + uw.postprocessing.geoid.spherical_shell_geoid_response( + radius_inner=0.55, + radius_outer=1.0, + harmonic_degree=harmonic_degree, + surface_topography_coefficient=0.4, + cmb_topography_coefficient=0.7, + ) + + +def test_rotated_adapter_rejects_degree_zero(): + with pytest.raises(ValueError, match="requires harmonic_degree >= 1"): + uw.postprocessing.geoid.spherical_shell_response_from_rotated_stokes( + stokes=object(), + radius_inner=0.55, + radius_outer=1.0, + harmonic_degree=0, + ) + + +def test_internal_load_requires_a_radius(): + with pytest.raises(ValueError, match="internal_load_radius is required"): + uw.postprocessing.geoid.spherical_shell_geoid_response( + radius_inner=0.55, + radius_outer=1.0, + harmonic_degree=2, + surface_topography_coefficient=0.4, + cmb_topography_coefficient=0.7, + internal_load_coefficient=1.0, + ) + + +def test_rotated_adapter_requires_explicit_self_gravity_parameters(): + with pytest.raises(ValueError, match="Self-gravity requires explicit values"): + uw.postprocessing.geoid.spherical_shell_response_from_rotated_stokes( + stokes=object(), + radius_inner=0.55, + radius_outer=1.0, + harmonic_degree=2, + include_self_gravity=True, + ) + + +def test_rotated_stokes_adapter_matches_zhong_table_2(): + radius_inner = 0.55 + radius_outer = 1.0 + rint = 0.775 + mesh = uw.meshing.SphericalShellInternalBoundary( + radiusOuter=radius_outer, + radiusInternal=rint, + radiusInner=radius_inner, + cellSize=0.25, + qdegree=2, + degree=1, + ) + velocity = uw.discretisation.MeshVariable("U_geoid", mesh, mesh.dim, degree=2) + pressure = uw.discretisation.MeshVariable("P_geoid", 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 + + theta = mesh.CoordinateSystem.xR[1] + unit_r = mesh.CoordinateSystem.unit_e_0 + load = sympy.assoc_legendre(2, 0, sympy.cos(theta)) * unit_r + stokes.add_natural_bc(load, "Internal") + stokes.add_rotated_freeslip_bc(0, "Upper", normal=unit_r) + stokes.add_rotated_freeslip_bc(0, "Lower", normal=-unit_r) + stokes.petsc_use_nullspace = True + stokes.petsc_options["snes_type"] = "ksponly" + stokes.tolerance = 1.0e-5 + stokes.solve() + + response = uw.postprocessing.geoid.spherical_shell_response_from_rotated_stokes( + stokes=stokes, + radius_inner=radius_inner, + radius_outer=radius_outer, + harmonic_degree=2, + internal_load_radius=rint, + internal_load_coefficient=1.0, + include_self_gravity=True, + surface_density_contrast=3300.0, + cmb_density_contrast=5400.0, + planet_radius=6370000.0, + gravity=9.8, + gravitational_constant=6.67e-11, + ) + + assert np.isclose(response.surface_topography, 0.41920, rtol=0.10) + assert np.isclose(response.cmb_topography, 0.77060, rtol=0.10) + assert np.isclose(response.surface_geoid, 0.02579, rtol=0.10) + assert np.isclose(response.cmb_geoid, 0.03206, rtol=0.10) + assert np.isclose(response.self_gravity.surface_topography, 0.49980, rtol=0.10) + assert np.isclose(response.self_gravity.cmb_topography, 0.93130, rtol=0.10) + assert np.isclose(response.self_gravity.surface_geoid, 0.04486, rtol=0.10) + assert np.isclose(response.self_gravity.cmb_geoid, 0.05461, rtol=0.10)