From 225a1c78e49454ca77d6948c67112b7aff8c6a40 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 25 Aug 2026 13:53:07 +1000 Subject: [PATCH 1/5] Recover pointwise traction with a mass that is not singular at vertices (#633) The grad-div penalty was blamed for a 28% loss in vertex-sampled spherical dynamic topography. It was not responsible. The de-smearing mass for a 3-D P2 TRIANGULAR trace has vertex rows that sum to exactly zero (row sums [0,0,0,60,60,60]). Those rows annihilate a constant, so solving M sigma = R amplifies any perturbation of the nodal load at vertices by O(1) -- and, being an instability rather than a discretisation error, independently of h. The recovery was already 7.6% low with no penalty at all; the penalty only made it large enough to fail a test whose 12% tolerance had been hiding it. Isolating it needed a case that was curved but not 3-D. A 2-D annulus reproduces the signature exactly and then parts company under refinement: its error falls ~O(h^2) while the shell's stays flat at ~0.28 over a 3.2x node-count range. The 2-D P2 LINE mass has vertex row sums of 5, which is why 2-D never showed this. Dimension, curvature and the rotated constraint were all red herrings. mass="auto" now selects the P1-projected recovery on a 3-D P2 trace. This is a trade, not a free win: "consistent" midpoint values are superconvergent (0.1-1.5% at every penalty) and the P1 projection gives that up. It is the default because its error CONVERGES -- worst node 0.170 -> 0.049 over cellSize 0.30 -> 0.16 -- where the consistent vertex error does not. A default should not be unstable. "consistent" remains right for a caller sampling only midpoints, and the docstrings now say which is which. FreeSurface already used the P1-projected recovery in 3-D, so production dynamic topography was never affected; the exposure was mass="auto". With that closed, DEFAULT_PENALTY goes to 10.0 on the #625 evidence (59 -> 18 Schur iterations per application under FMG, 21% faster). test_1018's spherical topography test is refined (cellSize 0.25 -> 0.13, 4.0s -> 14.7s) and its tolerances tightened from 0.10/0.12 to 0.03/0.05, set from measured discretisation error with ~2x headroom. It now asserts every node class rather than the aggregate: the failure was confined to one class and an aggregate assertion passed straight through it. Negative control at that resolution -- "consistent" gives 0.31525 at vertices and fails, "p1" gives 0.41604 and passes. test_p2_triangle_mass_has_zero_vertex_row_sums guards the identity itself, so the reason stays attached to the choice. Measurements: ~/+Simulations/topography_penalty_633/README.md Also filed #637 (3-D recovery supports only P1/P2 triangular traces, so this cannot be cross-validated against a second discretisation). 25 passed in test_1018. Underworld development team with AI support from Claude Code --- docs/developer/CHANGELOG.md | 39 +++++++++++++++ docs/developer/subsystems/rotated-freeslip.md | 21 ++++++++ .../cython/petsc_generic_snes_solvers.pyx | 29 +++++------ src/underworld3/systems/solvers.py | 20 ++++---- src/underworld3/utilities/boundary_flux.py | 22 ++++++++- src/underworld3/utilities/rotated_bc.py | 19 ++++---- tests/test_0206_automatic_penalty_pairing.py | 7 +-- tests/test_1018_rotated_freeslip.py | 48 +++++++++++++++---- 8 files changed, 160 insertions(+), 45 deletions(-) diff --git a/docs/developer/CHANGELOG.md b/docs/developer/CHANGELOG.md index 53c0126f..f87254ea 100644 --- a/docs/developer/CHANGELOG.md +++ b/docs/developer/CHANGELOG.md @@ -6,6 +6,45 @@ This log tracks significant development work at a conceptual level, suitable for ## 2026 Q3 (July – September) +### A Singular Recovery Mass, Mistaken for a Penalty Defect (August 2026) + +**The grad-div penalty default is now on (`Stokes.DEFAULT_PENALTY = 10.0`)**, and the +reason it was held off turned out to be a defect somewhere else entirely (#633). + +With the penalty at 10, the spherical dynamic topography recovered from the rotated +free-slip reaction dropped 28% at *vertices* while the facet-integrated value stayed +correct. The natural reading — that grad-div augmentation corrupts the de-smearing from +reaction loads to pointwise stress — was wrong. + +The de-smearing mass for a 3-D **P2 triangular** trace has vertex rows that sum to +**exactly zero**. Those rows annihilate a constant, so solving `M σ = R` amplifies any +perturbation of the nodal load at vertices by O(1) — and, being an instability rather +than a discretisation error, does so independently of mesh resolution. The recovery was +already 7.6% low with no penalty at all; the penalty only made it large enough to fail a +test whose 12% tolerance had been hiding it. + +The discrimination needed a case that was curved but not 3-D. A 2-D annulus reproduces +the signature exactly and then parts company under refinement: its error falls ~O(h²) +while the shell's stays flat at ~0.28 over a 3.2× node-count range. The 2-D P2 **line** +mass has positive vertex row sums, which is why 2-D never showed the defect and why +dimension, curvature and the rotated constraint were all red herrings. + +- `mass="auto"` now selects the monotone **P1-projected** recovery on a 3-D P2 trace. + This is a trade, not a free win: `"consistent"` midpoint values are superconvergent + (0.1–1.5% at every penalty) and the P1 projection gives that up. It is the default + because its error *converges* (worst node 0.170 → 0.049 over cellSize 0.30 → 0.16) + where the consistent vertex error does not. `"consistent"` remains the right choice + for a caller sampling only midpoints, and the docstrings now say which is which. +- `FreeSurface` already used the P1-projected recovery in 3-D, so production dynamic + topography was never affected. The exposure was `mass="auto"`. +- The spherical topography test is refined (cellSize 0.25 → 0.13) and its tolerances + tightened from 0.10/0.12 to 0.03/0.05, set from measured discretisation error with + ~2× headroom, and now assert every node class rather than the aggregate — the failure + was confined to one class and an aggregate assertion passed straight through it. +- Filed #637: 3-D recovery accepts only P1/P2 triangular traces, so dynamic topography + has exactly one supported discretisation there and cannot be cross-validated. That + blocked the P3/hex arm of this investigation. + ### The Free Surface Reaches the Spherical Shell (July 2026) **`uw.systems.FreeSurface` now runs in 3D on a spherical shell** — the same diff --git a/docs/developer/subsystems/rotated-freeslip.md b/docs/developer/subsystems/rotated-freeslip.md index 83606e1c..668947b7 100644 --- a/docs/developer/subsystems/rotated-freeslip.md +++ b/docs/developer/subsystems/rotated-freeslip.md @@ -324,6 +324,27 @@ the shared boundary-mass machinery in `utilities/boundary_flux.py`; `dynamic_topography_field` writes `h = −(σ_nn − σ̄_nn)/(Δρ g)` onto a surface field for the free-surface integrator. +### Which de-smearing mass (3-D) + +`mass="auto"` uses the **P1-projected** recovery on a 3-D P2 trace. The full +("consistent") P2 **triangle** mass has vertex rows summing to exactly zero, so +it is singular on constants along those rows and `M⁻¹` amplifies any +perturbation of the nodal load at *vertices* by O(1) — independently of +resolution. Measured on the Zhong l=2 shell (#633): consistent-mass vertex +values are 7.6% low with no grad-div penalty, 28% low at `penalty=10` and 79% +low at 100, and the error is flat across a 3.2× refinement. The 2-D P2 **line** +mass has positive vertex row sums and is unaffected — which is why this appears +only in 3-D, and not because of dimension or curvature as such. + +This is a trade rather than a strict improvement. Consistent-mass *midpoint* +values are superconvergent (0.1–1.5% at every penalty) and the P1 projection +gives that up, because under it the midpoints are the P1 interpolant of the +vertices. `"p1"` is the default because its error converges under refinement +(worst node 0.170 → 0.049 over cellSize 0.30 → 0.16) where the consistent +vertex error does not. **If you sample only midpoints, pass +`mass="consistent"`.** See also #404 (the vertex-integral checkerboard) and +#637 (only P1/P2 triangular traces are supported in 3-D at all). + ## Tests `tests/test_1018_rotated_freeslip.py` (serial: essential-equivalence, FMG, diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index c42bbb47..1b9892de 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -3168,11 +3168,12 @@ class SolverBaseClass(uw_object): number); for a **vector** solver the traction :math:`\sigma\cdot\hat n` (pass ``normal`` to get the scalar normal component :math:`\hat n\cdot\sigma\cdot\hat n`). - ``mass`` de-smears the nodal reaction with ``"lumped"`` or ``"consistent"`` - boundary mass. ``"auto"`` (default) selects lumped recovery for 2D P1/P2 - traces and 3D P1 triangles, and the consistent solve for 3D P2 triangles and - 2D traces of degree >= 3 (where row-sum lumping is respectively invalid and - only O(h) pointwise). + ``mass`` de-smears the nodal reaction with ``"lumped"``, ``"consistent"`` or + ``"p1"`` boundary mass. ``"auto"`` (default) selects lumped recovery for 2D + P1/P2 traces and 3D P1 triangles, P1-PROJECTED recovery for 3D P2 triangles + (row-sum lumping is invalid there and the consistent solve amplifies at + vertices — #633), and the consistent solve for 2D traces of degree >= 3 + (where lumping is only O(h) pointwise). ``remove_mean`` subtracts the boundary mean — leave ``False`` for a physical flux (the mean is the Nusselt number); ``True`` gives a gauge-free field. @@ -6453,13 +6454,13 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): :meth:`solve`. ``mass="auto"`` (default) uses lumped recovery for 2D traces and 3D P1 - triangles, and the consistent surface-mass solve for 3D P2 triangles. - Explicit ``"lumped"`` and ``"consistent"`` choices remain available where - mathematically valid, and ``"p1"`` selects P1-PROJECTED recovery on a 3D - P2 trace (edge-midpoint loads folded onto vertices, lumped P1 triangle - mass — sound where the consistent P2 path carries the vertex-integral - checkerboard; the FreeSurface default in 3D). Three-dimensional recovery - currently supports triangular P1/P2 traces only. + triangles, and P1-PROJECTED recovery for 3D P2 triangles (edge-midpoint + loads folded onto vertices, lumped P1 triangle mass). Explicit ``"lumped"`` + and ``"consistent"`` choices remain available where mathematically valid; + ``"consistent"`` is pointwise-exact on a P2 trace in exact arithmetic but + its zero vertex row sums amplify any load perturbation at VERTICES by O(1), + independently of h (#404, measured in #633). Three-dimensional recovery + currently supports triangular P1/P2 traces only (#637). .. warning:: On CURVED boundaries, P2 vertex values of :math:`\sigma_{nn}` converge @@ -6484,8 +6485,8 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): interior left untouched. ``buoyancy_scale`` is :math:`\Delta\rho\,g` (traction → length). - ``mass="auto"`` selects lumped recovery where valid and the consistent - surface-mass solve for 3D P2 triangles. Requires a prior + ``mass="auto"`` selects lumped recovery where valid and P1-projected + recovery for 3D P2 triangles. Requires a prior :meth:`add_rotated_freeslip_bc` on ``boundary`` and a completed :meth:`solve`. .. warning:: diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 8043975b..b752b9f5 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -2183,17 +2183,15 @@ def penalty(self, value): #: not inside P0, so the term does not vanish at the discrete solution -- #: but it converges away. ``penalty = 0`` restores the unaugmented operator. #: - #: **Held at 0 pending the pointwise-traction question.** At 10 the - #: spherical dynamic topography recovered from the rotated free-slip - #: reaction drops 0.4192 -> 0.3021, 28%, on the *vertex-sampled* value while - #: the facet-integrated value stays correct (``test_1018``). So augmentation - #: corrupts the de-smearing from reaction loads to pointwise stress — - #: presumably because lambda*mu*(div u) is non-zero cell-by-cell for P2-P0 - #: and averages out over a facet integral but not at a vertex. Dynamic - #: topography is the main product of that machinery, so the value stays 0 - #: until that is resolved; everything needed for the change is in place and - #: it is this constant. - DEFAULT_PENALTY = 0.0 + #: This was held at 0 while #633 was open — at 10 the spherical dynamic + #: topography recovered from the rotated free-slip reaction dropped + #: 0.4192 -> 0.3021 on the *vertex-sampled* value. That turned out not to be + #: the penalty: the 3-D P2 TRIANGLE de-smearing mass has vertex rows summing + #: to exactly zero, so ``mass="consistent"`` amplified the perturbation at + #: vertices by O(1) regardless of h, and was already 7.6% low at + #: ``penalty = 0``. ``mass="auto"`` now picks the monotone P1-projected + #: recovery there, which holds to 0.6% at 10 and 1.8% at 100. + DEFAULT_PENALTY = 10.0 def _apply_automatic_penalty(self): diff --git a/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py index c0b8e5b2..00db4b98 100644 --- a/src/underworld3/utilities/boundary_flux.py +++ b/src/underworld3/utilities/boundary_flux.py @@ -399,7 +399,27 @@ def coord(q): ) order = orders.pop() if mass == "auto": - mass = "consistent" if order == 2 else "lumped" + # A P2 TRIANGLE mass has vertex rows that sum to EXACTLY zero + # (_P2_TRIANGLE_MASS, row sums [0,0,0,60,60,60]), so those rows + # annihilate a constant and M^-1 amplifies any component of the + # load near that direction — at vertices only, by O(1), and + # independently of h. Measured on the Zhong l=2 shell (#633): + # 'consistent' recovers 0.38723 at penalty=0 and 0.08793 at + # penalty=100 against an analytic 0.41920, where 'p1' holds + # 0.40989 / 0.40232. + # + # The choice is NOT that p1 is more accurate — it is not. It + # trades: p1's midpoints are the P1 interpolant of its vertices, + # where 'consistent' midpoints are superconvergent (0.1-1.5% at + # every penalty). What decides it is that p1's error CONVERGES + # (worst node 0.170 -> 0.049 over cellSize 0.30 -> 0.16) while + # the consistent vertex error is FLAT at ~0.30 — an instability, + # not a discretisation error. A default should not be unstable. + # 'consistent' stays available, and is the right call for a + # caller who samples ONLY midpoints (#404 checkerboard aside). + # The 2-D P2 LINE mass has vertex row sums of 5 and needs none + # of this. + mass = "p1" if order == 2 else "lumped" if order == 2 and mass == "lumped": raise ValueError( "A 3D P2 triangular trace has zero row-sum mass at its " diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 8400a887..8f392d84 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -2045,16 +2045,19 @@ def boundary_normal_traction(solver, boundary, solve_result, mass="auto"): σ_nn is the boundary-mass de-smear of R. ``mass`` selects the de-smear: - * ``"auto"`` (default) — lumped for 2D traces and 3D P1 triangles, consistent for - 3D P2 triangles. + * ``"auto"`` (default) — lumped for 2D traces and 3D P1 triangles, P1-PROJECTED + for 3D P2 triangles. * ``"lumped"`` — the diagonal row-sum mass. It is monotone for supported traces, but invalid for 3D P2 triangles because their vertex row sums are exactly zero. - * ``"consistent"`` — the full trace mass. Pointwise-exact 3D P2 recovery, but - carries the vertex-integral checkerboard on P2 triangles (#404 hold). + * ``"consistent"`` — the full trace mass. Its 3D P2 MIDPOINT values are + superconvergent (0.1–1.5% on the Zhong l=2 shell at every penalty), and that is + the reason to choose it. Its VERTEX values are not usable: the zero vertex row + sums make M singular on constants there, so M⁻¹ amplifies any perturbation of + the load by O(1) independently of h (7.6% low at ``penalty=0``, 28% at 10, + 79% at 100) — the vertex-integral checkerboard (#404), measured in #633. * ``"p1"`` — P1-PROJECTED recovery on a 3D P2 trace (edge-midpoint loads folded - onto vertices, lumped P1 triangle mass). Sound where the consistent P2 path - checkerboards; the FreeSurface default in 3D. On a P1 trace, identical to - ``"lumped"``. + onto vertices, lumped P1 triangle mass). Monotone, and the ``"auto"`` choice on + a 3D P2 trace since #633. On a P1 trace, identical to ``"lumped"``. Parallel-safe: r_c is scattered to a local vector (ghosts included) and read by LOCAL section offset. In 3D, coordinate-keyed reactions and boundary elements are gathered @@ -2101,7 +2104,7 @@ def dynamic_topography_field(solver, boundary, solve_result, field, """Populate a scalar MeshVariable ``field`` with the dynamic topography :math:`h = -(\\sigma_{nn}-\\overline{\\sigma_{nn}})/(\\Delta\\rho\\,g)` on ``boundary``, recovered from the rotated-free-slip constraint reaction. ``mass="auto"`` uses - lumped recovery where valid and the consistent surface mass for 3D P2 triangles. + lumped recovery where valid and P1-projected recovery for 3D P2 triangles. Interior nodes are left untouched. Returns ``field``. This is the hand-off to the free-surface machinery: the 3-number topography diff --git a/tests/test_0206_automatic_penalty_pairing.py b/tests/test_0206_automatic_penalty_pairing.py index e53fca6e..55b78812 100644 --- a/tests/test_0206_automatic_penalty_pairing.py +++ b/tests/test_0206_automatic_penalty_pairing.py @@ -70,9 +70,10 @@ def test_the_penalty_default_does_not_depend_on_the_preconditioner(monkeypatch): the solver means the preconditioner changes the answer, and `test_1017`, `test_0835` and `test_0836` all assert it does not. - The default is patched to a non-zero value for the duration, so the test - still means something while `DEFAULT_PENALTY` is held at 0 — otherwise it - would pass by comparing 0 to 0 however the value were chosen. + The default is patched to a value that is not the shipped one, so the test + cannot pass by accident if `DEFAULT_PENALTY` is ever changed to whatever is + asserted here — and could not pass by comparing 0 to 0 back when it shipped + at 0. """ monkeypatch.setattr(uw.systems.Stokes, "DEFAULT_PENALTY", 7.0) diff --git a/tests/test_1018_rotated_freeslip.py b/tests/test_1018_rotated_freeslip.py index 5b122c84..bc84c3a7 100644 --- a/tests/test_1018_rotated_freeslip.py +++ b/tests/test_1018_rotated_freeslip.py @@ -611,7 +611,7 @@ def test_rotated_freeslip_spherical_shell_3d(): assert rotfrac < 1e-8, f"rotation mode {k} gauge {rotfrac:.2e} not removed" -def _spherical3d_reaction_topography(cell_size=0.25): +def _spherical3d_reaction_topography(cell_size=0.13): """Zhong l=2 topography recovered directly from rotated constraint reactions.""" radius_inner = 0.55 @@ -689,19 +689,51 @@ def fit(mask): ) +def test_p2_triangle_mass_has_zero_vertex_row_sums(): + """WHY mass='auto' is P1-projected on a 3D P2 trace, and not the consistent solve. + + The P2 TRIANGLE mass has vertex rows summing to exactly zero, so it is singular + on constants along those rows and M^-1 amplifies any perturbation of the nodal + load at VERTICES without bound -- by O(1), and independently of h (#633: the + vertex error is flat at ~0.28 over a 3.2x node-count range, while the same + measurement on a 2-D annulus converges away ~O(h^2)). The 2-D P2 LINE mass has + positive vertex row sums and needs none of this, which is why 2-D never showed + the defect. Guarding the identity here keeps the reason attached to the choice. + """ + from underworld3.utilities.boundary_flux import _P2_TRIANGLE_MASS + + row_sums = _P2_TRIANGLE_MASS.sum(axis=1) + assert np.allclose(row_sums[:3], 0.0), ( + f"P2 triangle vertex row sums {row_sums[:3]} are no longer zero — the " + "reason mass='auto' avoids the consistent solve has changed") + assert np.all(row_sums[3:] > 0.0), "P2 triangle midpoint rows should be positive" + + @pytest.mark.level_2 def test_rotated_freeslip_spherical3d_reaction_topography(): - """3D reaction loads must be divided by boundary mass to recover pointwise stress.""" + """3D reaction loads must be divided by boundary mass to recover pointwise stress. + + Every node class is asserted, because the failure this guards is confined to ONE + of them: under the consistent P2 surface mass the VERTEX values were 7.6% low with + no penalty at all and 28% low at penalty=10, while the midpoints stayed within 1% + and the facet-integrated value stayed correct (#633). An assertion on the + aggregate alone passes straight through that. + + Tolerances are the measured discretisation error at this resolution with ~2x + headroom, not round numbers: at cellSize=0.13 the P1-projected recovery lands + within 1.3% on the surface and 3.1% at the CMB. They were 0.10/0.12 at + cellSize=0.25, loose enough to pass a recovery that was 7.6% wrong. + """ surface, surface_vertices, surface_midpoints, cmb, cmb_vertices, cmb_midpoints = ( _spherical3d_reaction_topography() ) - assert np.isclose(surface, 0.41920, rtol=0.10) - assert np.isclose(cmb, 0.77060, rtol=0.10) - assert np.isclose(surface_vertices, 0.41920, rtol=0.12) - assert np.isclose(surface_midpoints, 0.41920, rtol=0.12) - assert np.isclose(cmb_vertices, 0.77060, rtol=0.12) - assert np.isclose(cmb_midpoints, 0.77060, rtol=0.12) + assert np.isclose(surface, 0.41920, rtol=0.03) + assert np.isclose(surface_vertices, 0.41920, rtol=0.03) + assert np.isclose(surface_midpoints, 0.41920, rtol=0.03) + assert np.isclose(cmb, 0.77060, rtol=0.05) + assert np.isclose(cmb_vertices, 0.77060, rtol=0.05) + assert np.isclose(cmb_midpoints, 0.77060, rtol=0.05) def test_rotated_freeslip_annulus_zero_leakage(): From cd592411e905a6462feae366634708606ee560d6 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 25 Aug 2026 16:21:15 +1000 Subject: [PATCH 2/5] Keep the penalty default at 0: the FMG win is not the default path (#625) Flipping Stokes.DEFAULT_PENALTY to 10 was justified by #625 (59 -> 18 Schur iterations per application under FMG, 21% faster) once #633 stopped blocking it. Measured against the tier A/B suite, it does not pay. Three tests fail at 10 and pass at 0, deterministically: test_0112 swarm accumulated strain goes negative test_1060 Nitsche free-slip leak 1.234e-4 against a 1e-4 bound test_1061 topography/normal-traction correlation 0.979 against 0.99 and the same three run 8.3x slower -- 9.27s to 76.92s, warm JIT cache both ways, back to back. The mechanism is the one UW3 already warns about: the FMG win needs a mesh hierarchy, and without refinement>=1 the velocity block falls back to GAMG, where grad-div augmentation is exactly what drives the solve into its iteration cap. #625 recorded that cost as 33%; on this evidence it is much worse. The default path is the one without a hierarchy, so the default has to serve it. Set penalty=10 explicitly on a solver that has FMG. Two of the three failures are worth separating from the performance case, because neither is a penalty defect: test_1060 is predictable rather than surprising -- augmentation perturbs a WEAKLY imposed constraint, and a Nitsche leak bound has no slack to absorb it. A strong rotated constraint is untouched. Noted in the docstring. test_0112 is not about the penalty at all. It accumulates 0.5*evaluate(Einv2, swarm.data) and asserts the total is non-negative, but evaluate returns -0.4976 at in-domain points against the left wall below the lid corner -- for sqrt((E**2).trace()/2), whose range is non-negative. That happens at penalty=0 too; the penalty only perturbs the field enough for the running total to cross zero. Filed as #641, with the extrapolation explanation ruled out: the 815 deliberately out-of-domain samples all came back positive. The #633 recovery fix is unaffected and stays. It is what actually closed that issue, and it is independent of this constant: at penalty=0 the tightened test_1018 tolerances hold with 1.6x-4.1x headroom, and the recovered values barely move between penalty 0 and 10, which is the point of the P1-projected recovery. Underworld development team with AI support from Claude Code --- src/underworld3/systems/solvers.py | 29 ++++++++++++++++++++--------- 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index b752b9f5..e8b82d84 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -2183,15 +2183,26 @@ def penalty(self, value): #: not inside P0, so the term does not vanish at the discrete solution -- #: but it converges away. ``penalty = 0`` restores the unaugmented operator. #: - #: This was held at 0 while #633 was open — at 10 the spherical dynamic - #: topography recovered from the rotated free-slip reaction dropped - #: 0.4192 -> 0.3021 on the *vertex-sampled* value. That turned out not to be - #: the penalty: the 3-D P2 TRIANGLE de-smearing mass has vertex rows summing - #: to exactly zero, so ``mass="consistent"`` amplified the perturbation at - #: vertices by O(1) regardless of h, and was already 7.6% low at - #: ``penalty = 0``. ``mass="auto"`` now picks the monotone P1-projected - #: recovery there, which holds to 0.6% at 10 and 1.8% at 100. - DEFAULT_PENALTY = 10.0 + #: **Held at 0, but no longer because of #633.** That issue turned out to be + #: the recovery operator, not the penalty: the 3-D P2 TRIANGLE de-smearing + #: mass has vertex rows summing to exactly zero, and ``mass="auto"`` now + #: picks the monotone P1-projected recovery instead. The topography + #: objection is gone. + #: + #: It stays at 0 because of the GAMG cost above, which is far worse than + #: #625 recorded. Flipping it to 10 was tried and reverted: three tier-A/B + #: tests failed (``test_0112`` swarm, ``test_1060`` Nitsche free-slip leak + #: 1.234e-4 against a 1e-4 bound, ``test_1061`` topography correlation + #: 0.979 against 0.99) and the same three ran **8.3x slower** -- 9.27 s to + #: 76.92 s, warm JIT cache both ways. The 21% FMG win needs a hierarchy; + #: without ``refinement>=1`` the velocity block falls back to GAMG, which + #: is the default path and the one that pays. A default has to serve the + #: default. + #: + #: Set ``penalty = 10`` explicitly on a solver that has an FMG hierarchy. + #: Note it also perturbs a Nitsche free-slip constraint, whose leak bound + #: is not slack enough to absorb it. + DEFAULT_PENALTY = 0.0 def _apply_automatic_penalty(self): From 7f1cf7a56a46183e62396025a0b986b4e4973d1a Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 25 Aug 2026 16:22:01 +1000 Subject: [PATCH 3/5] Correct the changelog: the penalty default stayed off The entry was written while the flip was in place and announced it as landed. It was reverted the same day on the GAMG evidence (8.3x slower on the three tests it broke), so the changelog claimed the opposite of what shipped. The #633 recovery fix it sits beside is unchanged. Underworld development team with AI support from Claude Code --- docs/developer/CHANGELOG.md | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/docs/developer/CHANGELOG.md b/docs/developer/CHANGELOG.md index f87254ea..30764903 100644 --- a/docs/developer/CHANGELOG.md +++ b/docs/developer/CHANGELOG.md @@ -8,8 +8,9 @@ This log tracks significant development work at a conceptual level, suitable for ### A Singular Recovery Mass, Mistaken for a Penalty Defect (August 2026) -**The grad-div penalty default is now on (`Stokes.DEFAULT_PENALTY = 10.0`)**, and the -reason it was held off turned out to be a defect somewhere else entirely (#633). +**The grad-div penalty default stays off**, but the reason it was held off turned out to +be a defect somewhere else entirely (#633) — so the objection that had blocked it is gone, +and a different one took its place. With the penalty at 10, the spherical dynamic topography recovered from the rotated free-slip reaction dropped 28% at *vertices* while the facet-integrated value stayed @@ -44,6 +45,19 @@ dimension, curvature and the rotated constraint were all red herrings. - Filed #637: 3-D recovery accepts only P1/P2 triangular traces, so dynamic topography has exactly one supported discretisation there and cannot be cross-validated. That blocked the P3/hex arm of this investigation. +- `Stokes.DEFAULT_PENALTY` was flipped to 10 on the #625 evidence and then **reverted**. + Three tier-A/B tests fail at 10 and pass at 0, and the same three run **8.3x slower** + (9.27 s to 76.92 s, warm cache both ways). The #625 win needs an FMG hierarchy; without + `refinement>=1` the velocity block falls back to GAMG, which is where grad-div + augmentation drives the solve into its iteration cap. The default path is the one + without a hierarchy, so the default serves it; set `penalty=10` explicitly where FMG + is available. +- Two of those three failures are not penalty defects. The Nitsche free-slip leak + (1.234e-4 against a 1e-4 bound) is augmentation perturbing a *weakly* imposed + constraint — a strong rotated constraint is untouched. The swarm one exposed #641: + `evaluate` returns −0.4976 for `sqrt((E**2).trace()/2)` at in-domain points near the + lid-corner singularity, at `penalty=0` as well; the penalty merely moved an accumulated + total across zero. ### The Free Surface Reaches the Spherical Shell (July 2026) From 424771ab6dd8f29fa040f1224675c4cc1e9494c0 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 25 Aug 2026 20:34:40 +1000 Subject: [PATCH 4/5] Rebuild the vertices from the midpoints instead of discarding both (#633) The zero row sums are a statement about the ELEMENT, not the code. Row sum i of a mass matrix is INT(phi_i), because the basis is a partition of unity: sum_j INT(phi_i phi_j) = INT(phi_i sum_j phi_j) = INT(phi_i) So [0,0,0,60,60,60] says the P2 triangle vertex basis function has ZERO MEAN. Its DOF carries no mass and an L2 recovery fixes it by cancellation. Reliable pointwise topography at those vertices was never something we failed to implement -- it is not available from an L2 recovery at all. mass="auto" therefore stops asking for it. On a 3-D P2 trace it now takes the consistent solve, keeps its superconvergent midpoints, and reconstructs the vertices from them: the three midpoints of a facet determine a unique linear function, so a vertex reads its two adjacent midpoints and subtracts the opposite one, averaged over incident facets. Worst-node error against the analytic Zhong coefficients, cellSize 0.25 -> 0.11: surface p1 .041 .026 .018 .013 .008 surface mid .016 .012 .012 .003 .004 cmb p1 .116 .067 .047 .030 .025 cmb mid .094 .058 .043 .024 .015 Better at every resolution on both boundaries, by 1.8x to 4.9x, and converging. The same ladder re-confirms why the consistent vertices cannot be used directly: .076 / .074 / .070 / .119 / .093 -- flat and erratic, no better at h=0.11 than at h=0.25. The previous commit made auto SAFE by switching to the P1-projected recovery. That was sound but it discarded the superconvergent midpoints along with the unusable vertices; this keeps them. "p1" remains available and unchanged. test_1018's surface tolerances tighten 0.03 -> 0.01 on the improvement (worst 0.26% at cellSize=0.13, ~4x headroom); the CMB stays at 0.05 (worst 2.4%). 25 passed. Caveat recorded in the helper: the reconstruction amplifies noise up to 3x by construction, which is tolerable only BECAUSE the midpoints are superconvergent. It would not transfer to an element whose midpoints are not. Measurements: ~/+Simulations/topography_penalty_633/midpoint_reconstruction.py Underworld development team with AI support from Claude Code --- docs/developer/CHANGELOG.md | 30 +++-- .../cython/petsc_generic_snes_solvers.pyx | 23 ++-- src/underworld3/utilities/boundary_flux.py | 107 ++++++++++++++---- src/underworld3/utilities/rotated_bc.py | 19 +++- tests/test_1018_rotated_freeslip.py | 15 ++- 5 files changed, 142 insertions(+), 52 deletions(-) diff --git a/docs/developer/CHANGELOG.md b/docs/developer/CHANGELOG.md index 30764903..e22e839f 100644 --- a/docs/developer/CHANGELOG.md +++ b/docs/developer/CHANGELOG.md @@ -30,16 +30,32 @@ while the shell's stays flat at ~0.28 over a 3.2× node-count range. The 2-D P2 mass has positive vertex row sums, which is why 2-D never showed the defect and why dimension, curvature and the rotated constraint were all red herrings. -- `mass="auto"` now selects the monotone **P1-projected** recovery on a 3-D P2 trace. - This is a trade, not a free win: `"consistent"` midpoint values are superconvergent - (0.1–1.5% at every penalty) and the P1 projection gives that up. It is the default - because its error *converges* (worst node 0.170 → 0.049 over cellSize 0.30 → 0.16) - where the consistent vertex error does not. `"consistent"` remains the right choice - for a caller sampling only midpoints, and the docstrings now say which is which. +- The row-sum result is a statement about the **element**, not the code. Row sum *i* of + a mass matrix is `∫φᵢ`, since `Σⱼ∫φᵢφⱼ = ∫φᵢ·Σⱼφⱼ = ∫φᵢ` for a partition of unity. So + `[0,0,0,60,60,60]` says **the P2 triangle vertex basis function has zero mean**. Its + DOF carries no mass; an L2 recovery fixes it by cancellation. Reliable pointwise + topography at those vertices is not something we failed to implement — it is not + available from an L2 recovery at all. +- So `mass="auto"` stops asking. On a 3-D P2 trace it now takes the **consistent** solve, + keeps its superconvergent midpoints, and **reconstructs the vertices from them**: the + three midpoints of a facet determine a unique linear function, so a vertex reads its + two adjacent midpoints and subtracts the opposite one, averaged over incident facets. + Worst-node error against the analytic coefficient, over cellSize 0.25 → 0.11: + + | | 0.25 | 0.20 | 0.16 | 0.13 | 0.11 | + |---|---|---|---|---|---| + | surface, P1-projected | 0.041 | 0.026 | 0.018 | 0.013 | 0.008 | + | surface, reconstructed | 0.016 | 0.012 | 0.012 | 0.003 | 0.004 | + | CMB, P1-projected | 0.116 | 0.067 | 0.047 | 0.030 | 0.025 | + | CMB, reconstructed | 0.094 | 0.058 | 0.043 | 0.024 | 0.015 | + + Better at every resolution on both boundaries, by 1.8x to 4.9x, and converging. The + simpler P1-projected recovery stays available as `mass="p1"` — it is sound, it just + discards the good data along with the bad. - `FreeSurface` already used the P1-projected recovery in 3-D, so production dynamic topography was never affected. The exposure was `mass="auto"`. - The spherical topography test is refined (cellSize 0.25 → 0.13) and its tolerances - tightened from 0.10/0.12 to 0.03/0.05, set from measured discretisation error with + tightened from 0.10/0.12 to 0.01/0.05, set from measured discretisation error with ~2× headroom, and now assert every node class rather than the aggregate — the failure was confined to one class and an aggregate assertion passed straight through it. - Filed #637: 3-D recovery accepts only P1/P2 triangular traces, so dynamic topography diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 1b9892de..b011e26d 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -3168,12 +3168,12 @@ class SolverBaseClass(uw_object): number); for a **vector** solver the traction :math:`\sigma\cdot\hat n` (pass ``normal`` to get the scalar normal component :math:`\hat n\cdot\sigma\cdot\hat n`). - ``mass`` de-smears the nodal reaction with ``"lumped"``, ``"consistent"`` or - ``"p1"`` boundary mass. ``"auto"`` (default) selects lumped recovery for 2D - P1/P2 traces and 3D P1 triangles, P1-PROJECTED recovery for 3D P2 triangles - (row-sum lumping is invalid there and the consistent solve amplifies at - vertices — #633), and the consistent solve for 2D traces of degree >= 3 - (where lumping is only O(h) pointwise). + ``mass`` de-smears the nodal reaction with ``"lumped"``, ``"consistent"``, + ``"p1"`` or ``"midpoint"`` boundary mass. ``"auto"`` (default) selects lumped + recovery for 2D P1/P2 traces and 3D P1 triangles, MIDPOINT-RECONSTRUCTED + recovery for 3D P2 triangles (row-sum lumping is invalid there and the + consistent solve amplifies at vertices — #633), and the consistent solve for 2D + traces of degree >= 3 (where lumping is only O(h) pointwise). ``remove_mean`` subtracts the boundary mean — leave ``False`` for a physical flux (the mean is the Nusselt number); ``True`` gives a gauge-free field. @@ -6454,9 +6454,10 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): :meth:`solve`. ``mass="auto"`` (default) uses lumped recovery for 2D traces and 3D P1 - triangles, and P1-PROJECTED recovery for 3D P2 triangles (edge-midpoint - loads folded onto vertices, lumped P1 triangle mass). Explicit ``"lumped"`` - and ``"consistent"`` choices remain available where mathematically valid; + triangles, and MIDPOINT-RECONSTRUCTED recovery for 3D P2 triangles (the + consistent solve, keeping its superconvergent midpoints, with vertices rebuilt + from them). ``"p1"`` selects the simpler P1-projected recovery; explicit + ``"lumped"`` and ``"consistent"`` remain available where mathematically valid; ``"consistent"`` is pointwise-exact on a P2 trace in exact arithmetic but its zero vertex row sums amplify any load perturbation at VERTICES by O(1), independently of h (#404, measured in #633). Three-dimensional recovery @@ -6485,8 +6486,8 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): interior left untouched. ``buoyancy_scale`` is :math:`\Delta\rho\,g` (traction → length). - ``mass="auto"`` selects lumped recovery where valid and P1-projected - recovery for 3D P2 triangles. Requires a prior + ``mass="auto"`` selects lumped recovery where valid and + midpoint-reconstructed recovery for 3D P2 triangles. Requires a prior :meth:`add_rotated_freeslip_bc` on ``boundary`` and a completed :meth:`solve`. .. warning:: diff --git a/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py index 00db4b98..ba99c55c 100644 --- a/src/underworld3/utilities/boundary_flux.py +++ b/src/underworld3/utilities/boundary_flux.py @@ -60,6 +60,44 @@ ) +def _rebuild_vertices_from_midpoints(elements, global_index, flux): + """Replace the VERTEX entries of a 3D P2 recovery with values reconstructed from + its edge-midpoint entries, averaged over the facets incident on each vertex. + + The P2 triangle vertex basis function has zero mean (``_P2_TRIANGLE_MASS`` vertex + row sums are exactly 0), so the vertex DOF of an L2 recovery carries no mass and is + fixed by cancellation — O(1) error, independent of h. The midpoints carry all of it + and are superconvergent, so the vertices are better obtained from them than asked + for directly. + + On one triangle the three midpoint values determine a unique linear function. + Because ``m01 = (v0+v1)/2`` and so on, that function takes the value + ``m01 + m20 - m12`` at ``v0`` — each vertex reads its two ADJACENT midpoints and + subtracts the OPPOSITE one. The subtraction amplifies noise up to 3x, which is + tolerable only because the midpoints are superconvergent; averaging over the + (typically ~6) incident facets damps it further. + + Mean removal is unaffected: it weights by ``M·1``, which is zero at exactly these + vertices, so the gauge is set by the midpoints either way — and the reconstruction + is linear, so removing the mean before or after gives the same answer. + """ + total = np.zeros_like(flux) + count = np.zeros(len(flux), dtype=np.int64) + for _order, nodes, _area in elements.values(): + v0, v1, v2 = (global_index[k] for k in nodes[:3]) + m01, m12, m20 = (global_index[k] for k in nodes[3:]) + for vertex, near_a, near_b, far in ( + (v0, m01, m20, m12), # v0 lies on edges 01 and 20 + (v1, m01, m12, m20), # v1 lies on edges 01 and 12 + (v2, m12, m20, m01)): # v2 lies on edges 12 and 20 + total[vertex] += flux[near_a] + flux[near_b] - flux[far] + count[vertex] += 1 + rebuilt = np.array(flux, dtype=float, copy=True) + touched = count > 0 + rebuilt[touched] = total[touched] / count[touched] + return rebuilt + + def _key(c, dim): return tuple(round(float(t), 9) for t in np.asarray(c).ravel()[:dim]) @@ -300,9 +338,10 @@ def _desmear(solver, boundary, xs, R, mass, remove_mean, partial_reaction=True, csec = dm.getCoordinateSection() cvec = np.asarray(dm.getCoordinatesLocal().array).reshape(-1, dim) v0, v1 = dm.getDepthStratum(0) - if mass not in ("auto", "lumped", "consistent", "p1"): - raise ValueError("mass must be 'auto', 'lumped', 'consistent', or 'p1' " - "(P1-projected recovery on a 3D P2 trace).") + if mass not in ("auto", "lumped", "consistent", "p1", "midpoint"): + raise ValueError("mass must be 'auto', 'lumped', 'consistent', 'p1' " + "(P1-projected) or 'midpoint' (midpoint-reconstructed) " + "— the last two apply to a 3D P2 trace.") if dim == 3: lsec = dm.getLocalSection() ncomp = lsec.getFieldComponents(0) @@ -400,26 +439,33 @@ def coord(q): order = orders.pop() if mass == "auto": # A P2 TRIANGLE mass has vertex rows that sum to EXACTLY zero - # (_P2_TRIANGLE_MASS, row sums [0,0,0,60,60,60]), so those rows - # annihilate a constant and M^-1 amplifies any component of the - # load near that direction — at vertices only, by O(1), and - # independently of h. Measured on the Zhong l=2 shell (#633): - # 'consistent' recovers 0.38723 at penalty=0 and 0.08793 at - # penalty=100 against an analytic 0.41920, where 'p1' holds - # 0.40989 / 0.40232. + # (_P2_TRIANGLE_MASS, row sums [0,0,0,60,60,60]). Row sum i of a + # mass matrix IS the integral of basis function i, because the + # basis is a partition of unity: + # sum_j INT(phi_i phi_j) = INT(phi_i sum_j phi_j) = INT(phi_i) + # so this says INT(phi_vertex) = 0 — a property of the P2 triangle + # element, not of this code. A vertex DOF carries no mass, so an + # L2 recovery fixes it by cancellation and amplifies whatever noise + # is present: O(1), and independently of h. Measured on the Zhong + # l=2 shell (#633), the consistent vertex error over cellSize + # 0.25 -> 0.11 runs 0.076 / 0.074 / 0.070 / 0.119 / 0.093 — flat + # and erratic, never converging. # - # The choice is NOT that p1 is more accurate — it is not. It - # trades: p1's midpoints are the P1 interpolant of its vertices, - # where 'consistent' midpoints are superconvergent (0.1-1.5% at - # every penalty). What decides it is that p1's error CONVERGES - # (worst node 0.170 -> 0.049 over cellSize 0.30 -> 0.16) while - # the consistent vertex error is FLAT at ~0.30 — an instability, - # not a discretisation error. A default should not be unstable. - # 'consistent' stays available, and is the right call for a - # caller who samples ONLY midpoints (#404 checkerboard aside). - # The 2-D P2 LINE mass has vertex row sums of 5 and needs none - # of this. - mass = "p1" if order == 2 else "lumped" + # The midpoints carry all of the mass (row sums 60) and are + # superconvergent (to 5e-4). So do not ask the recovery for vertex + # values at all: keep the midpoints and RECONSTRUCT the vertices + # from them ('midpoint'). Worst-node error against the analytic + # coefficient beats the P1-projected recovery at every resolution + # measured, on both boundaries, by 1.8x to 4.9x: + # h 0.25 0.20 0.16 0.13 0.11 + # p1 .041 .026 .018 .013 .008 (Upper) + # mid .016 .012 .012 .003 .004 + # p1 .116 .067 .047 .030 .025 (Lower) + # mid .094 .058 .043 .024 .015 + # 'p1' remains available and is sound; it simply discards the good + # data along with the bad. The 2-D P2 LINE mass has vertex row sums + # of 5 and needs none of this. + mass = "midpoint" if order == 2 else "lumped" if order == 2 and mass == "lumped": raise ValueError( "A 3D P2 triangular trace has zero row-sum mass at its " @@ -456,6 +502,13 @@ def coord(q): order = 1 mass = "lumped" + # 'midpoint' IS the consistent solve, plus a vertex reconstruction + # from its midpoint values afterwards. + reconstruct = mass == "midpoint" + if reconstruct: + mass = "consistent" if order == 2 else "lumped" + reconstruct = order == 2 # a P1 trace has no midpoints to use + keys = sorted(R_by) global_index = {key: i for i, key in enumerate(keys)} reaction = np.array([R_by[key] for key in keys], dtype=float) @@ -505,6 +558,9 @@ def coord(q): raise RuntimeError( "Consistent boundary-mass solve failed on boundary " f"{boundary!r}." ) + if reconstruct: + flux = _rebuild_vertices_from_midpoints( + elements, global_index, flux) if remove_mean: mean = float(np.dot(flux, boundary_mass) / np.sum(boundary_mass)) flux -= mean @@ -539,8 +595,11 @@ def value_at(key): raise NotImplementedError( f"Boundary-flux recovery is not implemented for mesh dimension {dim}." ) - if mass == "p1": - mass = "lumped" # P1 consumers read vertex values; vertex lumping is sound + if mass in ("p1", "midpoint"): + # Both exist to work around the 3D P2 TRIANGLE's zero-mean vertex basis. A 2D + # P2 LINE trace has vertex row sums of 5, so its vertex values are sound as + # recovered and neither workaround is needed here. + mass = "lumped" lsec = dm.getLocalSection() ncomp = lsec.getFieldComponents(0) diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 8f392d84..d9f84af0 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -2045,8 +2045,8 @@ def boundary_normal_traction(solver, boundary, solve_result, mass="auto"): σ_nn is the boundary-mass de-smear of R. ``mass`` selects the de-smear: - * ``"auto"`` (default) — lumped for 2D traces and 3D P1 triangles, P1-PROJECTED - for 3D P2 triangles. + * ``"auto"`` (default) — lumped for 2D traces and 3D P1 triangles, + MIDPOINT-RECONSTRUCTED for 3D P2 triangles. * ``"lumped"`` — the diagonal row-sum mass. It is monotone for supported traces, but invalid for 3D P2 triangles because their vertex row sums are exactly zero. * ``"consistent"`` — the full trace mass. Its 3D P2 MIDPOINT values are @@ -2056,8 +2056,16 @@ def boundary_normal_traction(solver, boundary, solve_result, mass="auto"): the load by O(1) independently of h (7.6% low at ``penalty=0``, 28% at 10, 79% at 100) — the vertex-integral checkerboard (#404), measured in #633. * ``"p1"`` — P1-PROJECTED recovery on a 3D P2 trace (edge-midpoint loads folded - onto vertices, lumped P1 triangle mass). Monotone, and the ``"auto"`` choice on - a 3D P2 trace since #633. On a P1 trace, identical to ``"lumped"``. + onto vertices, lumped P1 triangle mass). Monotone and sound, but it discards the + superconvergent midpoints along with the unusable vertices. On a P1 trace, + identical to ``"lumped"``. + * ``"midpoint"`` — the ``"auto"`` choice on a 3D P2 trace (#633). The consistent + solve, keeping its superconvergent midpoint values, with the VERTEX values + reconstructed from them: on each facet the three midpoints determine a unique + linear function, so a vertex reads its two adjacent midpoints and subtracts the + opposite one, averaged over incident facets. Beats ``"p1"`` on worst-node error + at every resolution measured (h 0.25 → 0.11), by 1.8x to 4.9x. On a P1 trace, + identical to ``"lumped"``. Parallel-safe: r_c is scattered to a local vector (ghosts included) and read by LOCAL section offset. In 3D, coordinate-keyed reactions and boundary elements are gathered @@ -2104,7 +2112,8 @@ def dynamic_topography_field(solver, boundary, solve_result, field, """Populate a scalar MeshVariable ``field`` with the dynamic topography :math:`h = -(\\sigma_{nn}-\\overline{\\sigma_{nn}})/(\\Delta\\rho\\,g)` on ``boundary``, recovered from the rotated-free-slip constraint reaction. ``mass="auto"`` uses - lumped recovery where valid and P1-projected recovery for 3D P2 triangles. + lumped recovery where valid and midpoint-reconstructed recovery for 3D P2 + triangles. Interior nodes are left untouched. Returns ``field``. This is the hand-off to the free-surface machinery: the 3-number topography diff --git a/tests/test_1018_rotated_freeslip.py b/tests/test_1018_rotated_freeslip.py index bc84c3a7..166f8d5c 100644 --- a/tests/test_1018_rotated_freeslip.py +++ b/tests/test_1018_rotated_freeslip.py @@ -703,6 +703,11 @@ def test_p2_triangle_mass_has_zero_vertex_row_sums(): from underworld3.utilities.boundary_flux import _P2_TRIANGLE_MASS row_sums = _P2_TRIANGLE_MASS.sum(axis=1) + # Row sum i IS the integral of basis function i: sum_j INT(phi_i phi_j) = + # INT(phi_i sum_j phi_j) = INT(phi_i), the basis being a partition of unity. + # So this asserts INT(phi_vertex) = 0 — the P2 triangle vertex basis has ZERO + # MEAN, which is why its DOF carries no mass and why 'auto' reconstructs the + # vertices from the midpoints instead of recovering them. assert np.allclose(row_sums[:3], 0.0), ( f"P2 triangle vertex row sums {row_sums[:3]} are no longer zero — the " "reason mass='auto' avoids the consistent solve has changed") @@ -720,17 +725,17 @@ def test_rotated_freeslip_spherical3d_reaction_topography(): aggregate alone passes straight through that. Tolerances are the measured discretisation error at this resolution with ~2x - headroom, not round numbers: at cellSize=0.13 the P1-projected recovery lands - within 1.3% on the surface and 3.1% at the CMB. They were 0.10/0.12 at + headroom, not round numbers: at cellSize=0.13 the midpoint-reconstructed recovery + lands within 0.26% on the surface and 2.4% at the CMB. They were 0.10/0.12 at cellSize=0.25, loose enough to pass a recovery that was 7.6% wrong. """ surface, surface_vertices, surface_midpoints, cmb, cmb_vertices, cmb_midpoints = ( _spherical3d_reaction_topography() ) - assert np.isclose(surface, 0.41920, rtol=0.03) - assert np.isclose(surface_vertices, 0.41920, rtol=0.03) - assert np.isclose(surface_midpoints, 0.41920, rtol=0.03) + assert np.isclose(surface, 0.41920, rtol=0.01) + assert np.isclose(surface_vertices, 0.41920, rtol=0.01) + assert np.isclose(surface_midpoints, 0.41920, rtol=0.01) assert np.isclose(cmb, 0.77060, rtol=0.05) assert np.isclose(cmb_vertices, 0.77060, rtol=0.05) assert np.isclose(cmb_midpoints, 0.77060, rtol=0.05) From e8dc357ce8ba6c5c980a3171d3a6490078ad7f23 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 26 Aug 2026 12:52:19 +1000 Subject: [PATCH 5/5] Credit #414 for the mechanism, and discharge its tolerance caveat explicitly The zero-mean P2 vertex basis was not a new finding: #414 documented it, and with a better mechanism than the one first written here. Because the vertex basis has zero surface mean, the vertex reaction carries essentially only the O(h) facet-normal/geometry error and the consistent solve faithfully reconstructs it -- it is not amplifying noise, as the earlier wording had it. #414 also asked that the test_1018 vertex goldens NOT be tightened, since at rtol 0.10-0.12 they were absorbing that bias. The previous commit tightened them. That is defensible only because the quantity under test changed -- mass='auto' no longer consumes the biased vertex reactions at all -- so the test now says so, and cites #414 rather than silently overriding it. The reconstruction is #414's own action item (2), 'consider having dynamic_topography on curved boundaries return a midpoint-weighted or fitted field', which had not been actioned. Underworld development team with AI support from Claude Code --- docs/developer/CHANGELOG.md | 13 +++++++------ tests/test_1018_rotated_freeslip.py | 12 ++++++++++-- 2 files changed, 17 insertions(+), 8 deletions(-) diff --git a/docs/developer/CHANGELOG.md b/docs/developer/CHANGELOG.md index e22e839f..d578ca50 100644 --- a/docs/developer/CHANGELOG.md +++ b/docs/developer/CHANGELOG.md @@ -30,12 +30,13 @@ while the shell's stays flat at ~0.28 over a 3.2× node-count range. The 2-D P2 mass has positive vertex row sums, which is why 2-D never showed the defect and why dimension, curvature and the rotated constraint were all red herrings. -- The row-sum result is a statement about the **element**, not the code. Row sum *i* of - a mass matrix is `∫φᵢ`, since `Σⱼ∫φᵢφⱼ = ∫φᵢ·Σⱼφⱼ = ∫φᵢ` for a partition of unity. So - `[0,0,0,60,60,60]` says **the P2 triangle vertex basis function has zero mean**. Its - DOF carries no mass; an L2 recovery fixes it by cancellation. Reliable pointwise - topography at those vertices is not something we failed to implement — it is not - available from an L2 recovery at all. +- The zero-mean P2 vertex basis was **already known and documented in #414**, which + recorded the same drift-away-under-refinement we re-measured here. Its mechanism is + the sharper one and is adopted: because the vertex basis has zero surface mean, the + vertex reaction carries essentially only the O(h) facet-normal/geometry error, and the + consistent solve *faithfully reconstructs that error* — it is not amplifying noise. + What #633 adds is the separation from the grad-div penalty (which was blamed for it) + and the fix below, which is #414's own unactioned recommendation (2). - So `mass="auto"` stops asking. On a 3-D P2 trace it now takes the **consistent** solve, keeps its superconvergent midpoints, and **reconstructs the vertices from them**: the three midpoints of a facet determine a unique linear function, so a vertex reads its diff --git a/tests/test_1018_rotated_freeslip.py b/tests/test_1018_rotated_freeslip.py index 166f8d5c..e137463c 100644 --- a/tests/test_1018_rotated_freeslip.py +++ b/tests/test_1018_rotated_freeslip.py @@ -726,8 +726,16 @@ def test_rotated_freeslip_spherical3d_reaction_topography(): Tolerances are the measured discretisation error at this resolution with ~2x headroom, not round numbers: at cellSize=0.13 the midpoint-reconstructed recovery - lands within 0.26% on the surface and 2.4% at the CMB. They were 0.10/0.12 at - cellSize=0.25, loose enough to pass a recovery that was 7.6% wrong. + lands within 0.26% on the surface and 2.4% at the CMB. + + #414 asked that these goldens NOT be tightened, because at rtol 0.10-0.12 they were + absorbing the curved-boundary vertex bias and tightening them would have been + tightening against the wrong reference. That caveat is discharged rather than + ignored: `mass="auto"` no longer consumes the biased vertex reactions at all. It + reconstructs vertices from the superconvergent midpoints, which is #414's own + action item (2) — "consider having dynamic_topography on curved boundaries return a + midpoint-weighted or fitted field". The quantity under test changed; the tolerance + follows it. """ surface, surface_vertices, surface_midpoints, cmb, cmb_vertices, cmb_midpoints = (