Grad-div penalty breaks POINTWISE traction recovery (28% on dynamic topography); facet-integrated is fine

[lmoresi]

The grad-div penalty corrupts the recovery of pointwise boundary traction from reaction loads, while leaving the facet-integrated value correct. This is the one thing blocking Stokes.DEFAULT_PENALTY from being turned on (PR #632, where everything else is in place).

Measured, test_1018_rotated_freeslip::test_rotated_freeslip_spherical3d_reaction_topography, spherical shell, rotated free-slip, penalty = 10:

quantity	expected	measured
surface (facet-integrated)	0.41920	0.4192	pass
surface_vertices (pointwise)	0.41920	0.30211	-28%
cmb (facet-integrated)	0.77060	0.7706	pass

At penalty = 0 all of them pass. That test's own docstring is "3D reaction loads must be divided by boundary mass to recover pointwise stress", which is exactly the step that stops working.

Not a p_mech bookkeeping error. The first hypothesis was that sigma_nn was being built from the Lagrange multiplier without the documented p_mech = p - lambda*mu*(div u) correction. It is not: boundary_normal_traction recovers sigma_nn from the Cartesian nodal reaction r_c = A.u - b, which already carries the penalty because the penalty is in the operator. Moving the term into the total stress (PR #632) left the measured value byte-identical at 0.30210742214244796.

Working hypothesis. lambda*mu*(div u) is non-zero cell-by-cell for P2-P0 — the discrete divergence is only weakly zero, against the pressure space — so it averages out over a facet integral but not at a vertex. The de-smearing from reaction loads to pointwise stress divides by a boundary mass that does not account for the augmented contribution.

Why it matters. Dynamic topography is the main product of this machinery, and a 28% error in it is not something to ship silently behind a default. It is also worth knowing independently of the penalty: if the pointwise recovery is sensitive to a term the integral is not, that is a property of the recovery, not of the penalty.

Suggested first step. Compare surface_vertices against surface_midpoints as lambda is swept 0, 1, 10, 100 — the test already computes both, and whether they diverge together or separately should distinguish a mass-weighting error from a sampling one.

Context: the penalty is measured to be a strong lever on the Schur count under FMG (Schur iterations per application 59 -> 18, total velocity work 546 -> 270, 21% faster) and the accuracy cost elsewhere is a consistent perturbation — same convergence rate, gap shrinking 1.102 -> 1.087 -> 1.066 over three refinements. So this is the last thing standing between that and a default.

Underworld development team with AI support from Claude Code

[lmoresi]

It is sampling, not mass weighting. Midpoints are essentially immune; vertices collapse. Sweep of Stokes.DEFAULT_PENALTY through the existing helper (_spherical3d_reaction_topography, Zhong l=2, spherical shell, cellSize 0.25):

lambda	surf (integrated)	surf (vertices)	surf (midpoints)	cmb (integrated)	cmb (vertices)	cmb (midpoints)
0	0.41151	0.38723	0.41966	0.79375	0.83904	0.77845
1	0.40825	0.37650	0.41890	0.79242	0.83470	0.77814
10	0.38787	0.30211	0.41664	0.80087	0.86581	0.77893
100	0.33128	0.08793	0.41292	0.87294	1.14231	0.78193

Reference: surface 0.41920, CMB 0.77060.

Midpoints barely move: surface 0.41966 -> 0.41292 (-1.6%) and CMB 0.77845 -> 0.78193 (+0.4%) across a hundredfold change in lambda. They are also the most accurate recovery at every lambda, including lambda=0.

Vertices diverge hard, and with opposite signs at the two boundaries: surface -77%, CMB +36%. The sign flip tracks the outward normal (+r at the surface, -r at the CMB), which is what an unaccounted radial term lambda*mu*(div u) n_hat in the nodal reaction would do.

So the earlier hypothesis — that _desmear divides by a boundary mass which ignores the augmented contribution — is ruled out. A mass-weighting error would move vertices and midpoints together, since both are divided by the same primitive. It does not.

What distinguishes the two is the geometry of the node, not the weighting. On P2 velocity the boundary carries vertices and edge midpoints; a vertex is shared by many incident cells and an edge midpoint by few. The grad-div term is volumetric, so the reaction at a node accumulates it from every incident cell — many at a vertex, two or three at a midpoint. That is the shape of the discrepancy and the next thing to test.

Worth recording separately: the integrated value is not the best recovery here. At lambda=0 the midpoint value (0.41966) is closer to 0.41920 than the facet-integrated one (0.41151). The test's tolerances (rtol 0.10 on integrated, 0.12 on sampled) hide that.

Also noted: this configuration is Taylor-Hood, P2-P1 continuous — not P2-P0. So the effect is not confined to a discontinuous pressure space.

Script: ~/+Simulations/topography_penalty_633/vertices_vs_midpoints.py.

Underworld development team with AI support from Claude Code

[lmoresi]

2-D does not reproduce the effect at all, so it cannot test the hypothesis. Reporting this as a null result rather than reading the small numbers as support.

Unit box, rotated free-slip on the top wall, P2-P1, smooth body force, RMS normal traction over the top boundary nodes split by node type:

element	lambda	rms(vertices)	rms(other)	n_vert	n_other
simplex	0	0.20451	0.23909	13	12
simplex	10	0.19924	0.24219	13	12
quad	0	0.21060	0.23612	13	12
quad	10	0.19816	0.24311	13	12

Ratios from lambda 0 -> 10: simplex vertices x0.974, quad vertices x0.941. Against the spherical 3-D case, where the same lambda takes vertices to x0.78 and lambda=100 to x0.23.

So the phenomenon is simply absent here — nothing moves by more than 6%, and if anything the quad vertices move slightly more than the simplex ones, which is the opposite of what the incident-cell hypothesis predicts. With 13 vertices against 12 edge nodes the sample is also far too small to read a 3% difference. The incident-cell hypothesis is untested, not supported.

That leaves the discriminating variable as whatever differs between the two configurations: three dimensions, the curved boundary, the spherical rotated constraint, or the internal-boundary loading.

Next step, and it is available: a 2-D annulus with rotated free-slip. That is curved and uses the rotated constraint but is 2-D, so it separates "3-D" from "curvature" — and unlike P3 or hexes it is supported by _desmear. If the annulus reproduces the vertex collapse, curvature (or the rotated constraint on a curved wall) is the variable and dimension is not; if it does not, the effect is specific to 3-D.

Also confirmed: the cross-element validation asked for cannot be run in 3-D. _desmear refuses anything else:

boundary_flux.py:331  NotImplementedError: 3D boundary-flux recovery supports P1 or P2 triangular traces.
boundary_flux.py:327  NotImplementedError: 3D boundary-flux recovery currently requires triangular facets.

P3/P2 and hexes both stop there. That is worth its own attention independently of the penalty: in 3-D, dynamic topography has exactly one supported discretisation.

Scripts: ~/+Simulations/topography_penalty_633/ (vertices_vs_midpoints.py, which_nodes_break.py, two_d_quads.py).

Underworld development team with AI support from Claude Code

[lmoresi]

Resolved: it is the de-smearing operator, not the penalty

The 3-D P2 triangular surface mass has vertex rows that sum to exactly zero
(boundary_flux.py:51 — row sums [0, 0, 0, 60, 60, 60]). Those rows annihilate a
constant, so solving M sigma = R amplifies whatever component of the nodal load lies
near that direction — without bound, at vertices only, and in a way refinement does not
remove. The penalty supplies such a component; it does not create the defect.

_desmear already half-knows this. It refuses mass="lumped" on a 3-D P2 trace for
precisely this reason (boundary_flux.py:405), and boundary_normal_traction's
docstring warns that consistent "carries the vertex-integral checkerboard on P2
triangles (#404 hold)". That warning turns out to describe the failure exactly.

Same solves, same reactions, two recovery operators

Judged against the analytic Zhong coefficient 0.41920, fixed before the run:

mass	penalty	all	vertices	midpoints
consistent	0	0.41151	0.38723	0.41966
consistent	10	0.38787	0.30211	0.41664
consistent	100	0.33128	0.08793	0.41292
p1	0	0.40390	0.40989	0.40189
p1	10	0.40140	0.40731	0.39943
p1	100	0.39671	0.40232	0.39483

mass="p1" holds to 0.6% at penalty=10 and 1.8% at penalty=100.

Note the penalty=0 row: consistent vertices are already 7.6% below the analytic
value with no penalty at all, where p1 is 2.2% below. The rtol=0.12 in
test_rotated_freeslip_spherical3d_reaction_topography hid that. The penalty did not
introduce this error, it made it large enough to fail a test.

How it was isolated

A 2-D annulus — curved, rotated free-slip on both arcs, cos(2 theta) internal load —
reproduces the sign pattern exactly (outer vertices down, outer midpoints up, inner
vertices up) but ~30x weaker. Matching its resolution and de-smearing mass to the
shell's closes most of that gap: 12.5% at penalty=10, cellSize=0.25, consistent
mass, against the shell's 22%. So the earlier "2-D does not reproduce" reading was
comparing a much finer mesh recovered with a different mass.

The two then part company under refinement, and that is the discriminating measurement:

	cellSize
annulus, vertex deficit	0.1245 (h=0.25)	0.0538 (0.15)	0.0091 (0.075)	0.0046 (0.0375)
shell, vertex deficit	0.2842 (h=0.30)	0.2550 (0.25)	0.2651 (0.20)	0.2793 (0.16)

The annulus converges away ~O(h^2). The shell is flat over a 3.2x node-count range. Not
a consistency error, and not curvature — the annulus is curved and recovers. The one
structural difference that survives is the mass:

2-D P2 LINE mass,     vertex row  [4, 2, -1] * h/30   row sum   5h/30 > 0
3-D P2 TRIANGLE mass, vertex row  [6,-1,-1, 0,-4, 0]  row sum   EXACTLY 0

Consequences

Production dynamic topography was never affected. mass="p1" is already the
FreeSurface default in 3-D. The exposure is mass="auto", which selects consistent
for a 3-D P2 trace — what a direct boundary_normal_traction caller gets.

Stokes.DEFAULT_PENALTY is no longer blocked by this issue. It still ships at 0.0;
the case for flipping it is the solver one in #625 and can now be argued on its own
merits.

Open, and a judgement call rather than a measurement: whether mass="auto" should
select "p1" instead of "consistent" for a 3-D P2 trace (boundary_flux.py:402, one
line). The evidence supports it — p1 is more accurate at every penalty tested,
including zero — but it moves shipped numbers and narrows the default from a
pointwise-exact P2 recovery to a P1-projected one. If it is taken, test_1018's
rtol=0.12 should tighten with it.

Also filed: #637 — 3-D recovery accepts only P1/P2 triangular traces, so there is no
second discretisation to cross-check any of this against. That blocked the P3/hex arm
and is why the discrimination had to be done in 2-D.

Scripts, logs and the full chain of reasoning: ~/+Simulations/topography_penalty_633/
(README.md).