diff --git a/docs/developer/CHANGELOG.md b/docs/developer/CHANGELOG.md index 53c0126f..d578ca50 100644 --- a/docs/developer/CHANGELOG.md +++ b/docs/developer/CHANGELOG.md @@ -6,6 +6,76 @@ 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 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 +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. + +- 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 + 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.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 + 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) **`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 3b3da609..fef4e13b 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"``, + ``"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. @@ -6455,13 +6456,14 @@ 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 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 + currently supports triangular P1/P2 traces only (#637). .. warning:: On CURVED boundaries, P2 vertex values of :math:`\sigma_{nn}` converge @@ -6486,8 +6488,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 + 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/systems/solvers.py b/src/underworld3/systems/solvers.py index 8043975b..e8b82d84 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -2183,16 +2183,25 @@ 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. + #: **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 diff --git a/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py index c0b8e5b2..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) @@ -399,7 +438,34 @@ 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]). 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 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 " @@ -436,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) @@ -485,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 @@ -519,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 47a06683..bb8c1cac 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -2084,16 +2084,27 @@ 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, + 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. 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 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 @@ -2140,7 +2151,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 the consistent surface mass 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_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..e137463c 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,64 @@ 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) + # 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") + 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 midpoint-reconstructed recovery + 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 = ( _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.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) def test_rotated_freeslip_annulus_zero_leakage():