Rebuild P2 boundary-recovery vertices from the midpoints, and exonerate the grad-div penalty (#633) - #649
Rebuild P2 boundary-recovery vertices from the midpoints, and exonerate the grad-div penalty (#633)#649lmoresi wants to merge 6 commits into
Conversation
…es (#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
) 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
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
) 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
…licitly 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
Adversarial review of this PRPosting the case against it, since nobody else has the context to and it should not merge on the strength of its own description. The headline mechanism was not discovered here, and the first version of this branch overclaimed it. #414 had the zero-mean P2 vertex basis, the drift-away-under-refinement, and the recommendation that became the fix. The commit history shows the correction ( The reconstruction inherits its safety from the midpoints.
The tolerance tightening overrides an explicit request in #414 not to tighten these goldens. The justification — that The CMB is the weak arm throughout. Reconstructed worst-node error is 0.015 at cellSize 0.11 against 0.004 at the surface, and The penalty revert rests on a small sample. The 8.3x figure is three tests — the three that failed, so selected for penalty sensitivity. It was measured back-to-back with a warm cache and I believe it, but it is not a suite-wide benchmark, and the earlier attempt at one was contaminated and discarded rather than repeated. What I am confident in: the zero row sums (exact, arithmetic); the h-independence of the consistent vertex error over five resolutions on the shell against convergence on the annulus; the negative control, which fails the new tolerance with the old recovery; and tier A+B at 1963 passed / 0 failed. Underworld development team with AI support from Claude Code |
Closes #633.
The penalty was innocent
With
stokes.penalty = 10, spherical dynamic topography recovered from the rotatedfree-slip reaction was 28% low at vertices while the facet-integrated value stayed
correct. Grad-div augmentation was the obvious suspect. It is not responsible.
The de-smearing mass for a 3-D P2 triangular trace has vertex rows summing to exactly
zero. Row sum i of a mass matrix is
∫φᵢ, sinceΣⱼ∫φᵢφⱼ = ∫φᵢ·Σⱼφⱼ = ∫φᵢfor apartition of unity — so
[0,0,0,60,60,60]says the P2 triangle vertex basis function haszero surface mean.
That was already known: #414 documented it and gave the sharper mechanism, which this
PR adopts rather than re-deriving. Because the vertex basis has zero 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 is new here is
the separation from the penalty, and the fix — which is #414's own action item (2),
"consider having
dynamic_topographyon curved boundaries return a midpoint-weighted orfitted field", previously unactioned.
How it was isolated
A 2-D annulus — curved, rotated free-slip,
cos(2θ)internal load — reproduces the signpattern exactly but ~30x weaker. Matching its resolution and mass to the shell's closes
most of that gap (12.5% vs 22% at penalty=10, cellSize 0.25). The two then part company
under refinement, and that is the discriminating measurement:
The annulus converges ~O(h²); the shell is flat over a 3.2x node-count range. So it is
neither a consistency error nor curvature — the annulus is curved and recovers. The 2-D P2
line mass has vertex row sums of 5, which is why 2-D never showed it.
The fix
mass="auto"on a 3-D P2 trace no longer asks for vertex values at all. It takes theconsistent 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:
p1p1Better at every resolution on both boundaries, by 1.8x to 4.9x.
mass="p1"remainsavailable and unchanged.
It is not only topography: on a flat boundary with a varying flux (Poisson heat
flux, manufactured harmonic solution) the reconstruction is ~3x better at vertices than
consistent— 0.169/0.124/0.084 against 0.459/0.358/0.265 over cellSize 0.45 → 0.22. Notethis sharpens #414's "exact on flat boundaries" claim, which holds for a uniform flux
(what
test_1019checks) but not a varying one.The penalty default was flipped and reverted
#633 no longer blocks
Stokes.DEFAULT_PENALTY, so it was set to 10 on the #625 evidenceand then reverted in the same branch. 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 JIT cache both ways. The
#625 win needs an FMG hierarchy; without
refinement>=1the velocity block falls back toGAMG, 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.
Two of those three failures are not penalty defects and are recorded as such:
test_1060(Nitsche free-slip leak 1.234e-4 against 1e-4) — augmentation perturbing aweakly imposed constraint. A strong rotated constraint is untouched.
test_0112(swarm) — actually evaluate() returns large negative values for a non-negative expression (sqrt invariant) at in-domain points near a boundary singularity #641:evaluatereturns −0.4976 forsqrt((E**2).trace()/2)at in-domain points near the lid-corner singularity, atpenalty=0as well. The penalty only moved an accumulated total across zero.Tests
test_1018refined cellSize 0.25 → 0.13 and tolerances tightened 0.10/0.12 → 0.01/0.05,set from measured error with ~2-4x headroom. Every node class is asserted, not the
aggregate — the failure was confined to one class and an aggregate assertion passed
straight through it.
bias. That caveat is discharged, not ignored:
autono longer consumes the biasedvertex reactions, so the quantity under test changed. The test says so and cites Curved-boundary P2 VERTEX flux/traction values converge only slowly — document and steer consumers to midpoint/fitted quantities #414.
consistentgives 0.31525 at the surface vertices andfails the new tolerance; the reconstruction gives 0.41604 and passes.
test_p2_triangle_mass_has_zero_vertex_row_sums— a no-solve guard on the identity, sothe reason stays attached to the choice.
development.Scope and what is left
This changes
mass="auto", which is what a direct caller ofboundary_flux/boundary_normal_tractiongets. It does not touch the free-surface path:free_surface.py:146already translates tomass="p1"in 3-D and consumes a P1 field.Deliberately not addressed here:
supported discretisation in 3-D and none of this can be cross-validated against a second
element. That is the real remaining gap.
boundary_normal_traction_integral, which contracts the reaction with the test functionand never forms a pointwise field. This PR improves the pointwise default; it does not
compete with that.
Measurements and scripts:
~/+Simulations/topography_penalty_633/.Underworld development team with AI support from Claude Code