From fd0a3bccd8e9f073272ba97ec4bff57ef897d4a8 Mon Sep 17 00:00:00 2001 From: Tyagi Date: Sat, 5 Sep 2026 13:27:24 +1000 Subject: [PATCH] Add spherical-shell reaction projection using the upstream traction integral Extract only the outstanding spherical postprocessing changes from PR #646 onto upstream development at 9da04b72. Reuse boundary_normal_traction_integral already provided by #648; no solver, Cython, SLCN/SUPG, checkpoint, or evaluation changes are included. Add opt-in projection=reaction to contract the normal reaction with the axisymmetric harmonic and normalize by its matching finite-element boundary inner product. Preserve centroid as the compatibility default. Carry the projection-validation and Zhong serial/MPI regressions, including near-machine-precision rank comparisons, and update the postprocessing documentation and changelog. Post-extraction build and regression results are recorded separately in the benchmark repository. --- docs/developer/CHANGELOG.md | 3 ++ ...ry-stress-and-projection-postprocessing.md | 22 +++++++-- src/underworld3/postprocessing/geoid.py | 45 ++++++++++++++++--- ...est_1071_spherical_shell_geoid_parallel.py | 3 +- tests/test_1070_postprocessing_geoid.py | 30 +++++++++++++ 5 files changed, 91 insertions(+), 12 deletions(-) diff --git a/docs/developer/CHANGELOG.md b/docs/developer/CHANGELOG.md index d1af2c6f8..04289452b 100644 --- a/docs/developer/CHANGELOG.md +++ b/docs/developer/CHANGELOG.md @@ -387,6 +387,9 @@ component exactly — correct on curved, tilted, and deformed boundaries (#293). 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. +- The spherical-shell geoid adapter accepts `projection="reaction"` to use + the same fitted integral without pointwise P2 recovery or a rank-zero + surface triangulation; `projection="centroid"` remains the default. - `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 6aef79807..114e7e783 100644 --- a/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md +++ b/docs/developer/subsystems/boundary-stress-and-projection-postprocessing.md @@ -97,13 +97,27 @@ response = uw.postprocessing.geoid.spherical_shell_response_from_rotated_stokes( planet_radius=6370000.0, gravity=9.8, gravitational_constant=6.67e-11, + projection="reaction", ) ``` -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. +The adapter supports two projection paths. `projection="centroid"` retains the +original pointwise-recovery workflow: recover `sigma_nn`, gather the samples to +rank zero, reconstruct a spherical triangulation, and integrate centroid +values. `projection="reaction"` contracts the assembled normal-reaction load +directly with the harmonic test function through +`Stokes.boundary_normal_traction_integral()`. The latter is distributed, avoids +the rank-zero surface reconstruction, and is an integral/fitted quantity rather +than a consumer of the slowly converging P2 vertex values on curved boundaries +(#414). Its fitted coefficient uses the matching discrete boundary norm, not an +analytical spherical norm, so the numerator and denominator share the same +faceted geometry. + +Both paths reuse the existing rotated-free-slip reaction; neither implements a +second CBF, constrained-multiplier, or topography recovery. `centroid` remains +the compatibility default while the direct reaction path accumulates benchmark +coverage. `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 diff --git a/src/underworld3/postprocessing/geoid.py b/src/underworld3/postprocessing/geoid.py index 3a2e18eb6..8b4a3dd7f 100644 --- a/src/underworld3/postprocessing/geoid.py +++ b/src/underworld3/postprocessing/geoid.py @@ -892,11 +892,36 @@ def _rotated_topography_coefficient( harmonic_degree: int, buoyancy_scale: float, response_sign: float, + projection: str, ) -> 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.") + if projection == "reaction": + import sympy + from underworld3.maths import BdIntegral + + theta = stokes.mesh.CoordinateSystem.xR[1] + harmonic = sympy.assoc_legendre(harmonic_degree, 0, sympy.cos(theta)) + traction_integral = stokes.boundary_normal_traction_integral( + boundary, + harmonic, + remove_mean=True, + ) + # The reaction is the load functional assembled on the faceted FE + # boundary. Use the matching discrete inner product for the fitted + # coefficient; an analytical spherical norm would mix geometries and + # introduce a chord-area bias, especially on the smaller CMB. + harmonic_norm = float( + BdIntegral(stokes.mesh, fn=harmonic**2, boundary=boundary).evaluate() + ) + return float( + -response_sign * traction_integral / (buoyancy_scale * harmonic_norm) + ) + if projection != "centroid": + raise ValueError("projection must be 'centroid' or 'reaction'.") + coords, sigma_nn = stokes.boundary_normal_traction(boundary, mass="auto") local_rows = np.column_stack( ( @@ -958,16 +983,18 @@ def spherical_shell_response_from_rotated_stokes( planet_radius: float | None = None, gravity: float | None = None, gravitational_constant: float = 6.67430e-11, + projection: str = "centroid", ) -> 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. + ``projection="centroid"`` (default) recovers pointwise traction and fits it + over a triangulation of the boundary samples. ``projection="reaction"`` + contracts the assembled nodal reaction directly with the harmonic test + function. The latter is a distributed weak/integral quantity that avoids + consuming slowly converging P2 vertex values on curved boundaries (issue + #414). 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( @@ -977,6 +1004,8 @@ def spherical_shell_response_from_rotated_stokes( ) if not isinstance(include_self_gravity, bool): raise TypeError("include_self_gravity must be True or False.") + if projection not in ("centroid", "reaction"): + raise ValueError("projection must be 'centroid' or 'reaction'.") if degree == 0: raise ValueError( "The rotated-Stokes adapter requires harmonic_degree >= 1 because " @@ -1005,6 +1034,7 @@ def spherical_shell_response_from_rotated_stokes( degree, surface_buoyancy_scale, 1.0, + projection, ) cmb_topography = _rotated_topography_coefficient( stokes, @@ -1013,6 +1043,7 @@ def spherical_shell_response_from_rotated_stokes( degree, cmb_buoyancy_scale, -1.0, + projection, ) geoid = spherical_shell_geoid_response( radius_inner=ri, diff --git a/tests/parallel/test_1071_spherical_shell_geoid_parallel.py b/tests/parallel/test_1071_spherical_shell_geoid_parallel.py index 189d61217..9ad0b6e15 100644 --- a/tests/parallel/test_1071_spherical_shell_geoid_parallel.py +++ b/tests/parallel/test_1071_spherical_shell_geoid_parallel.py @@ -56,6 +56,7 @@ def test_rotated_spherical_shell_geoid_matches_serial_reference(): planet_radius=6370000.0, gravity=9.8, gravitational_constant=6.67e-11, + projection="reaction", ) values = np.array( [ @@ -71,7 +72,7 @@ def test_rotated_spherical_shell_geoid_matches_serial_reference(): ) for rank_values in uw.mpi.comm.allgather(values): - assert np.array_equal(rank_values, values) + np.testing.assert_allclose(rank_values, values, rtol=1.0e-13, atol=1.0e-14) zhong_table_2 = np.array( [0.41920, 0.77060, 0.02579, 0.03206, 0.49980, 0.93130, 0.04486, 0.05461] diff --git a/tests/test_1070_postprocessing_geoid.py b/tests/test_1070_postprocessing_geoid.py index dce646a86..784b9440c 100644 --- a/tests/test_1070_postprocessing_geoid.py +++ b/tests/test_1070_postprocessing_geoid.py @@ -182,6 +182,17 @@ def test_rotated_adapter_requires_explicit_self_gravity_parameters(): ) +def test_rotated_adapter_rejects_unknown_projection(): + with pytest.raises(ValueError, match="projection must be"): + uw.postprocessing.geoid.spherical_shell_response_from_rotated_stokes( + stokes=object(), + radius_inner=0.55, + radius_outer=1.0, + harmonic_degree=2, + projection="nodal", + ) + + def test_rotated_stokes_adapter_matches_zhong_table_2(): radius_inner = 0.55 radius_outer = 1.0 @@ -229,6 +240,21 @@ def test_rotated_stokes_adapter_matches_zhong_table_2(): gravity=9.8, gravitational_constant=6.67e-11, ) + reaction_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, + projection="reaction", + ) assert np.isclose(response.surface_topography, 0.41920, rtol=0.10) assert np.isclose(response.cmb_topography, 0.77060, rtol=0.10) @@ -238,3 +264,7 @@ def test_rotated_stokes_adapter_matches_zhong_table_2(): 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) + assert np.isclose(reaction_response.surface_topography, 0.41920, rtol=0.03) + assert np.isclose(reaction_response.cmb_topography, 0.77060, rtol=0.03) + assert np.isclose(reaction_response.surface_geoid, 0.02579, rtol=0.03) + assert np.isclose(reaction_response.cmb_geoid, 0.03206, rtol=0.03)