Rotated free-slip: the vertex normal is a bisector but assembly is facet-wise, so a kinked boundary with unequal facets leaves the pressure gauge unpinned

[lmoresi]

Rotated strong free-slip loses the pressure gauge as soon as the boundary is curved. On a straight-walled box the constant pressure is pinned to machine zero and the solve is reproducible; deform the top of the same box and the pressure level runs away to ~1e4 while the pressure variation is 6.6e-2, and the answer stops being unique: a change of two machine epsilons in the mesh coordinates moves the velocity by 13%.

This is not a parallel effect (serial, np=1), not the workspace cache (the solver below is built fresh, after the deform), and not a convergence failure — every solve reports converged with a relative residual around 3e-10 against a tolerance of 1e-9.

Measured

Same 10x10 StructuredQuadBox, same load, penalty = 0, tolerance = 1e-9, rotated free-slip on all four walls, petsc_use_pressure_nullspace = True. Each configuration solved twice, the second time with every coordinate multiplied by 1 + 2*eps (a relative change of 4.44e-16):

configuration	mean(p)	std(p)	rel. move in v under a 2-eps coordinate change
rotated free-slip, straight walls	6.50e-14	6.42e-2	5.9e-15
rotated free-slip, top deformed by 0.02 y sin(pi x)	+1.86e+04	6.64e-2	1.33e-01
essential free-slip, same deformed mesh	-1.06e-16	6.73e-2	4.8e-13

The essential free-slip control is the point: on the same deformed mesh, with the same nullspace flag, the constant is pinned to 1e-16 and the answer is stable to 5e-13. Only the rotated path loses it.

The pressure constant is not merely large, it is arbitrary — the two runs above differ by +1.69e4 in the constant, and the size of the jump does not track the size of the perturbation (4.4e-16, 2.2e-15, 1e-14, 1e-12 and 1e-9 coordinate changes all move the velocity by between 4e-2 and 2e-1). That is the signature of a direction the operator does not pin at all: its amplitude is set by round-off, not by the data.

Decomposing the difference between two such solves: 99.95% of it is off the rigid-body span, so it is not the rigid-rotation gauge (which _mode_satisfies_constraints correctly rejects on a box — measured violation 3.045e-01 against a threshold of 1e-8, a margin of seven orders of magnitude). It is the constant-pressure mode, which drags the velocity with it because penalty = 0 leaves the constraint rows to couple them.

Why it matters beyond a test

• Any rotated free-slip model on a curved, deformed or adapted boundary — annulus, spherical shell, free surface, adapted meshes — has a pressure field offset by an arbitrary constant of order 1e4, and a velocity carrying a round-off-determined component of order 1e-1.
• Dynamic topography and boundary normal traction are recovered from the constraint reaction of this solve (boundary_normal_traction, dynamic_topography_field), which is exactly the workload that runs on curved boundaries.
• It makes results irreproducible across platforms: two builds that assemble the same operator to the last bit agree exactly, and two that differ in the last bit disagree at 1e-1.

How it surfaced

tests/test_1018_rotated_freeslip.py::test_rotated_workspace_deform_invalidates (from #543) compares a post-deform solve against a fresh solver on the deformed mesh and asserts they agree to 1e-6. That assertion is only satisfiable when the two assemblies agree bitwise. On macOS/arm64 they do, and err is exactly 0.0 in 81 consecutive runs across two PETSc toolchains and nine PYTHONHASHSEED values. On CI they do not, and the test fails intermittently with err around 7e-2 to 1.2e-1 — the size of the unpinned component, not a drift.

The CI failure reports are consistent with this: both solves converge from an identical |r0| = 0.02785134866629431 but to different residuals (1.7517680045617482e-11 and 1.9700774924116398e-11), where locally both give 8.267546030330813e-12. Perturbed runs here reproduce that residual range exactly (1.71e-11 to 1.79e-11).

Not caused by the branch that found it

Checked directly: the probe below is self-contained and runs unchanged against the merge base. Built at e475246f and at bugfix/issue-551-locator, the output is identical to every digit, including mean(p) = +1.860217e+04 and the 1.325836e-01 move. The point-locator work in #556 is not involved (measured separately: during the deform, 8 location calls over 2270 points, 38 returning -1, and 0 answers changed by the new rejection radius).

Reproduce

import numpy as np, sympy, underworld3 as uw

def rotated_solve(tag, deform, coord_scale=1.0):
    mesh = uw.meshing.StructuredQuadBox(
        elementRes=(10, 10), minCoords=(0, 0), maxCoords=(1, 1), qdegree=3)
    c = mesh.X.coords.copy()
    if deform:
        c[:, 1] += 0.02 * c[:, 1] * np.sin(np.pi * c[:, 0])
    if deform or coord_scale != 1.0:
        mesh.deform(c * coord_scale)
    v = uw.discretisation.MeshVariable(f"vR{tag}", mesh, mesh.dim, degree=2)
    p = uw.discretisation.MeshVariable(f"pR{tag}", mesh, 1, degree=1, continuous=False)
    s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p)
    s.constitutive_model = uw.constitutive_models.ViscousFlowModel
    s.constitutive_model.Parameters.shear_viscosity_0 = 1.0
    x, y = mesh.X
    s.bodyforce = sympy.Matrix([[0.0, sympy.sin(sympy.pi * x) * sympy.cos(sympy.pi * y)]])
    s.penalty = 0.0
    s.tolerance = 1e-9
    for wall in ("Top", "Bottom", "Left", "Right"):
        s.add_rotated_freeslip_bc(0, wall)
    s.petsc_use_pressure_nullspace = True
    s.solve()
    return np.asarray(v.data).copy(), np.asarray(p.data).ravel().copy()

eps = np.finfo(float).eps
for label, deform in (("straight", False), ("deformed", True)):
    out = [rotated_solve(f"{label}{i}", deform, cs)
           for i, cs in enumerate((1.0, 1.0 + 2 * eps))]
    (v0, p0), (v1, p1) = out
    print(label, "mean(p) =", p0.mean(),
          "| 2-eps move in v =", np.linalg.norm(v1 - v0) / np.linalg.norm(v0))

Suggested next steps

1. Find out why petsc_use_pressure_nullspace stops binding once the boundary is curved — whether the constant-pressure vector is still in the operator's null space after the rotation/elimination, or whether the nullspace is attached to a matrix the rotated path replaces.
2. Until then, test_rotated_workspace_deform_invalidates is asserting a property the solver does not have on a deformed mesh. It should either pin the gauge before comparing (compare velocity with the unpinned component projected out, and pressure with the constant removed) or be marked xfail against this issue — a green from it currently means "these two assemblies happened to agree bitwise", not "the solver is right".

Related: #348 (rotated LU path, pressure datum pin partition-dependent at np>1) is the parallel cousin of this; #543 introduced the test that exposed it.

Underworld development team with AI support from Claude Code

[lmoresi]

Blast radius: much smaller than the headline. Priority "important", not "drop everything".

The question was whether the unpinned mode reaches users' results or is masked downstream. Both — the consumers do remove a mean, and the configurations people actually run are not affected at all. Measured, not reasoned from the code.

The consumers do remove the mean

boundary_normal_traction recovers σ_nn with remove_mean=True hardcoded (rotated_bc.py:1886), and dynamic_topography_field is defined on σ_nn − mean(σ_nn). So a pure constant never reaches a user. That alone would not settle it, because the unpinned mode is not a pure constant — it carries a pressure shape component (≈5% of the pressure variation) and a velocity component (13%) as well. So the recovery was measured directly.

The dynamic-topography workload is unaffected

Annulus, radiusOuter=1.0, radiusInner=0.55, cellSize=0.1, rotated free-slip on both boundaries with the analytic radial normal, a degree-4 buoyancy load, solved twice with every coordinate multiplied by 1 + 2*eps:

quantity	value	2-eps change
mean(p)	−2.73e-04 (std 4.08e-02)	−7.0e-15
velocity	—	4.46e-13 relative
pressure shape (constant removed)	—	6.2e-12 relative
recovered σ_nn on Upper (128 nodes, range ±0.11)	mean 2.4e-19 after removal	4.84e-13 relative, max pointwise 8.7e-14

The gauge machinery works there: rotation_gauge_removed=True, one mode offered and admitted, which is correct for a closed circular free-slip boundary. The pressure constant is pinned. The mean removal is not even load-bearing on the annulus — the underlying solution is already reproducible to 5e-13.

Which configurations actually lose the gauge

Same 10×10 box, same load, penalty=0, tolerance=1e-9, petsc_use_pressure_nullspace=True; each solved twice with a 2-eps coordinate change:

#	configuration	mean(p)	2-eps move in v
A	straight box, rotated free-slip ×4	6.50e-14	5.9e-15
B	deformed top, rotated free-slip ×4	+1.86e+04	1.33e-01
C	deformed top, rotated free-slip on the 3 straight walls, top traction-free	1.37e-01	2.2e-14
D	deformed top, essential free-slip ×4	−1.06e-16	4.8e-13
E	annulus, rotated free-slip on both curved boundaries	−2.73e-04	4.5e-13

Only B is broken. The trigger is narrower than "a curved boundary": it needs rotated free-slip imposed on a deformed boundary and on the straight walls that meet it, enclosing the domain. Curvature alone is fine (E). A deformed boundary is fine as long as it is not itself a rotated wall (C — which is the free-surface shape: the moving boundary carries a traction condition, not a rotated constraint). Essential free-slip on the deformed boundary is fine (D).

Which is why the benchmarks passed

This explains the Y20 spherical relaxation and Crameri results directly rather than leaving them as a puzzle:

• spherical / annulus relaxation is configuration E — closed curved boundaries, gauge admitted and removed, stable to 5e-13;
• the free-surface benchmarks are configuration C — the deforming boundary carries the surface traction, and the rotated free-slip lives on the straight side walls;
• every rotated free-slip test on an undeformed box is configuration A.

Nothing in the validated set is configuration B, so nothing in the validated set was wrong.

Revised severity

Important, not urgent. No shipped benchmark or documented workflow is in the affected configuration, and the two curved-boundary workloads that matter for dynamic topography are reproducible to 5e-13. What remains is real and should still be fixed:

• configuration B is a legitimate thing for a user to build — a fully enclosed rotated free-slip box whose boundary has been moved by deform or by adaptation — and it silently produces a pressure field offset by ~1e4 and a velocity carrying a 13% round-off-determined component, with converged=True and a relative residual of 3e-10;
• it is a reproducibility hazard: two builds that assemble the same operator to the last bit agree exactly, two that differ in the last bit disagree at 1e-1;
• the narrowing above is itself a strong hint at the mechanism — the gauge survives when the rotated boundary set is closed and smooth (E) or axis-aligned (A), and fails when a tilted-normal wall meets axis-aligned ones, which points at the corner/constraint interaction rather than at the nullspace attachment in general.

Reproducer for the table above is the same script as in the issue body with the four config variants; the annulus measurement uses solver.boundary_normal_traction("Upper") on the two perturbed solves and compares node-for-node after sorting by polar angle.

Underworld development team with AI support from Claude Code

[lmoresi]

The mechanism: not the corners, and not the curvature — a kinked rotated node with unequal facets

The narrowing in the previous comment pointed at "a tilted-normal wall meeting axis-aligned
ones", i.e. at the corner. We built the experiment that separates the two candidates and
the corner is exonerated. What breaks is the nullspace content: the exact
constant-pressure vector is no longer a null vector of the constrained rotated operator,
but it is still the vector UW3 attaches. The wrong subspace is projected out, the RHS keeps
a component the solve cannot remove, and the amplitude along the true null direction is set
by round-off.

The residual is generated at boundary vertices where two facets of different
orientation meet, because the node gets a single nodal normal (the bisector of its facet
normals) while the assembler integrates facet by facet. It is proportional to
sin(Δ/2) · (|f1| − |f2|) — kink angle times facet-length difference — which we confirm
three ways: it grows as amplitude³ (§4), equalising the facets at identical curvature
kills it one-for-one (§5), and re-weighting the nodal normal by facet length makes it
machine zero (§6). A uniform annulus is clean because its facets are equal, not because it
is smooth; skew it and the defect appears there too (§7).

Everything below is measured on the issue's own reproducer (10x10 StructuredQuadBox,
penalty=0, tolerance=1e-9, petsc_use_pressure_nullspace=True, rotated free-slip x4),
serial, at origin/development (ffbcf98).

---

1. The decisive experiment: taper the corner, keep the curvature

Configuration B deforms the top by 0.02*y*sin(pi x). That displacement vanishes at
x=0 and x=1 but its slope does not, so each top corner is a tilted-normal wall meeting an
axis-aligned one (3.54 deg of tilt at res=10).

Two tapered variants remove that, keeping the top curved in its interior:

• T — 0.02*y*sin(pi x)**2: zero value and zero slope at both side walls.
• T2 — 0.02*y*bump(x), a C-infinity bump supported on (0.2, 0.8): the top is
  exactly flat for x<0.2 and x>0.8, so the corner node's Top normal is exactly
  (0,+1) — bit-identical corner geometry to the straight box A.

config	corner Top normal	mean(p)	2-eps move in v
A straight	(+0.0000, +1.0000)	6.50e-14	5.90e-15
B tilted corner	(-0.0617, +0.9981)	+1.86e+04	1.33e-01
T tapered corner	(-0.0191, +0.9998)	+1.25e+04	2.88e-04
T2 flat near corner	(+0.0000, +1.0000)	+1.43e+03	3.66e-06

T2 has A's corner exactly and B's disease. The pressure gauge is lost by three to four
orders of magnitude in every deformed case. The corner geometry does not decide it.

(The 2-eps move in v column is not the right severity metric — see §4: it is the
round-off amplitude of an unpinned direction and is chaotic. The gauge loss, mean(p),
is the stable symptom, and it is present in all three deformed configurations.)

---

2. The null space, computed rather than assumed

We captured the operator the rotated path actually hands to the KSP —
Ahat = zeroRowsColumns(Q J Q^T, normal_rows, diag) — together with the PETSc.NullSpace
built by _rotated_nullspace, both in the rotated frame, and took a dense SVD
(N=1282). Instrumentation was a monkeypatch on the probe side; src/ is untouched.

In every configuration the operator has a genuine null space of dimension exactly 1,
and that null vector is (numerically) pure pressure — velocity content 4e-16 to 9e-14 of
its norm. So the dimension is right and rigid rotations are not involved. What is wrong
is the direction.

config	sv_min/smax	|A z|/|A| for z = attached constant-p	true null vector's pressure part, std/|mean|	principal angle(attached, true null)
A straight	1.86e-16	2.73e-18	9.14e-12	5.24e-10 deg (sin 9.14e-12)
B tilted corner	4.44e-16	5.90e-10	1.93e-07	1.11e-05 deg (sin 1.94e-07)
T tapered corner	4.44e-16	2.14e-09	5.90e-07	3.39e-05 deg (sin 5.91e-07)
T2 flat near corner	2.72e-15	1.89e-08	5.05e-06	2.90e-04 deg (sin 5.06e-06)

Read it as: on the straight box the attached vector is the null vector to machine
precision. On every deformed box it is not, and the true null vector's pressure part is no
longer constant. The misalignment gets larger, not smaller, as the corner is tapered —
T2, with corners identical to A's, is the worst of the four.

The consequence is on the right-hand side. PETSc removes the attached vector's component
from the RHS exactly (|b.z|/|b| = 0.0 measured, all configs), but the component along the
true left null vector — the one that actually has to go for A x = b to be solvable — is
left in:

config	|b . l_true| / |b| (unremoved inconsistency)
A straight	2.64e-14
B tilted corner	1.66e-09
T tapered corner	4.51e-09
T2 flat near corner	3.72e-08

That is the defect stated exactly: the wrong subspace is being projected out, so the
solve is left with an inconsistent RHS and an unpinned direction. Decomposing the
difference between two 2-eps-perturbed solves onto the true null vector gives a component
of 1.0000 (to four decimals) in B, T and T2 — the run-to-run difference is the null
vector, nothing else.

(One column in that decomposition probe, "off both", was an artifact: the attached vector
and the true null vector are nearly parallel, so subtracting both removes the same
component twice. Ignore it; the "along n_true" column is the measurement.)

---

3. Corner census, and where the residual actually lives

For every boundary velocity node we counted the boundaries claiming it, the normals they
supply, the rank of the accumulated normal set, and the rows build_rotation constrains.
The corner handling is identical in all four configurations: 80 boundary velocity
nodes, 4 of them claimed by two boundaries, each with rank 2 and therefore
rows_constrained = 2/2 — a corner is pinned to v = 0. Nothing is overwritten, dropped,
or composed differently between A and B. There is no corner census finding.

The sharp version is where r = A z lives (z = the attached constant-pressure vector).
|A z| is entirely on velocity rows (pressure-row share 0.0000) and, in every deformed
case, 100% of |A z|^2 sits on non-corner Top boundary nodes:

config	corner nodes	other boundary nodes	interior
A straight	0.0000	0.0207	0.9793 (all of it at 1e-18 — noise)
B tilted corner	0.0000	1.0000	0.0000
T tapered corner	0.0000	1.0000	0.0000
T2 flat near corner	0.0000	1.0000	0.0000

Caveat, stated so it is not over-read: a corner node has both its rows in normal_rows, so
zeroRowsColumns guarantees zero residual there by construction. The corner's exoneration
does not rest on that row — it rests on T2, whose corner geometry is bit-identical to A's
and which is nonetheless the worst case in the table.

The residual concentrates on boundary vertex nodes where two facets of different
orientation meet — never on the P2 edge-midpoint nodes, which lie inside a single facet.
That is the mechanism in one line:

A boundary vertex gets one nodal normal (the normalised sum of its adjacent facet
normals — a bisector, not an area weighting), so its tangential row is orthogonal to the
average normal but to neither facet's actual normal. The constant-pressure column then
contributes ∮ φ_i (t̂_i · n) ds = sin(Δ/2) * (|f1| - |f2|) / 6 to that free tangential
row, where Δ is the kink angle at the node and |f1|, |f2| the two facet lengths. On a
flat wall Δ = 0. On an edge-midpoint node Δ = 0. On a curved wall with unequal facets it
is nonzero, and the constant pressure stops being a null vector.

The predicted scaling sin(Δ/2) * (|f1| - |f2|) is confirmed by §4 and, node by node,
by T2's residual distribution (largest at x = 0.2 and 0.8 where the flat part meets the
bump, exactly zero at the apex x = 0.5 where the two facets are equal by symmetry — the
ratio 0.366/0.134 between the x=0.2 and x=0.4 nodes matches the formula's 2.94 against
a measured 2.74).

---

4. Amplitude sweep: smooth growth from zero, no threshold

Deformation amplitude swept 0 -> 0.05 at res=10, both shapes:

shape	amp	corner tilt	max kink	|A z|/|A|	sin(angle)	|b.l_true|/|b|	mean(p)	2-eps dv
B	0.0000	0.00 deg	0.00 deg	2.73e-18	1.49e-08	2.64e-14	6.50e-14	5.90e-15
B	0.0005	0.09	0.03	9.37e-15	2.98e-08	6.48e-14	2.45e-09	4.92e-15
B	0.0020	0.35	0.11	5.99e-13	0.0	1.63e-12	3.93e-08	7.54e-15
B	0.0050	0.89	0.28	9.34e-12	0.0	2.55e-11	2.01e+02	6.47e-05
B	0.0100	1.77	0.56	7.44e-11	2.11e-08	2.06e-10	2.55e+05	2.63e-01
B	0.0200	3.54	1.12	5.90e-10	1.95e-07	1.67e-09	1.86e+04	1.33e-01
B	0.0500	8.78	2.80	8.92e-09	2.99e-06	2.65e-08	3.05e+03	6.96e-09
T	0.0005	0.03	0.05	3.39e-14	0.0	7.80e-14	3.91e-09	6.19e-15
T	0.0020	0.11	0.22	2.17e-12	0.0	4.49e-12	6.29e-08	6.76e-15
T	0.0050	0.27	0.55	3.38e-11	0.0	7.05e-11	4.90e+02	1.16e-03
T	0.0100	0.55	1.09	2.70e-10	7.60e-08	5.64e-10	5.44e+04	9.79e-02
T	0.0200	1.09	2.19	2.14e-09	5.92e-07	4.51e-09	1.25e+04	2.88e-04
T	0.0500	2.73	5.45	3.24e-08	9.06e-06	6.90e-08	8.35e+02	1.22e-07

Two things to read off:

1. |A z| grows smoothly from zero as amplitude^3 — 4x amplitude gives 64x, 2x gives
   8x, at every step of both sweeps. No threshold, no switch. That is exactly the
   sin(Δ/2) * (|f1|-|f2|) law: the kink angle Δ is O(amp), the facet-length difference is
   O(amp^2). A discrete misapplication that switches on at some corner angle would not
   produce a clean cubic. The |b.l_true| inconsistency tracks it one-for-one.
2. mean(p) and the 2-eps velocity move are chaotic, not monotone — B at amp 0.05
   moves by 7e-9 while B at amp 0.01 moves by 2.6e-1. They are the round-off amplitude of
   an unpinned direction, so they carry the presence of the defect but not its size. Quote
   |A z|/|A| or |b.l_true|/|b| when measuring a fix, not the velocity move.

sin(angle) reads 0.0 at small amplitudes only because the misalignment falls below the
resolution of a dense SVD whose smallest singular value is already at machine level; the
|A z| column is the sensitive instrument there.

At amp 0.01 the rotated solver also printed
WARNING: nonlinear rotated free-slip did NOT converge in 7 iterations (rel |F| = 2.94e-09)
— the inconsistent RHS floors the residual, which is the same signature as the CI reports
in the issue body (converged-but-to-different-residuals).

Resolution: it converges away, so it is a consistency error

At amp = 0.02, shape B:

res	max kink	|A z|/|A|	|b.l_true|/|b|	mean(p)	2-eps dv
5	2.08 deg	6.16e-09	6.37e-09	1.06e+04	1.08e-10
10	1.12	5.90e-10	1.67e-09	1.86e+04	1.33e-01
16	0.70	1.16e-10	6.57e-10	3.07e+00	2.52e-06

|A z| falls roughly as h^3.3 under refinement. This is a consistency error of the
discrete rotated constraint, not an O(1) misapplication — but it never reaches zero, and
because it is an unpinned direction rather than a bounded error, a 1e-10 inconsistency
still produces an O(1e4) pressure offset. Refining the mesh does not make the answer
reproducible; it only changes the size of the seed.

---

5. Confirming the law: equal facets clean it up, curvature and all

The law says the two facet contributions at a vertex cancel exactly when the facets are
the same length, however curved the boundary is. Test: keep configuration B's top
surface exactly, re-space its top vertices to equal arc length. Same curve, same kink
angles, equal facets. Control: the same non-uniform column spacing with a flat top
(unequal facets, no kink).

config	top facet-length spread	max kink	|A z|/|A|	mean(p)	2-eps dv
B_ref (uniform x)	1.858e-03	1.1217 deg	5.902e-10	1.86e+04	1.33e-01
B_equal_arc	1.537e-05	1.1227 deg	4.897e-12	4.06e+00	1.19e-07
A_nonuniform_x (flat top)	1.843e-03	0.0000 deg	2.985e-18	6.51e-14	3.81e-15

Equalising the facets at identical curvature drops |A z| by 120.6x while the
facet-length spread drops by 120.9x — one for one, exactly as r ∝ sin(Δ/2)·(|f1|-|f2|)
predicts. And unequal facets alone, with no kink, are completely harmless (2.99e-18,
machine zero). It takes both the kink and the length inequality.

That is why configuration E passes: a uniform annulus has equal facets everywhere, so
the two terms cancel at every node however curved the boundary is. Curvature was never the
trigger; non-uniformity of a curved rotated wall is.

---

6. The fix the law implies, measured

The residual at a boundary vertex row (2-D, P2 velocity, ∫_f φ_vertex = |f|/6) is

r_i = (1/6) [ |f1| (t̂_i · n1) + |f2| (t̂_i · n2) ]

which vanishes identically iff the nodal normal is the length-weighted (area-weighted
in 3-D) average n_i ∝ |f1| n1 + |f2| n2. _boundary_velocity_nodes currently normalises
each facet normal to unit length before accumulating (rotated_bc.py:171-172, 187), so
the nodal normal is the bisector n1 + n2 — right only when |f1| = |f2|.

Monkeypatching that one weighting on the probe side (src untouched: multiply ne by the
computeCellGeometryFVM volume before accumulating):

config	|A z|/|A| now	with length weighting	mean(p) now	with weighting	2-eps dv now	with weighting
A straight	2.73e-18	2.73e-18	6.50e-14	6.50e-14	5.90e-15	5.90e-15
B tilted corner	5.90e-10	2.65e-18	1.86e+04	4.05e-06	1.33e-01	7.08e-15
T tapered corner	2.14e-09	2.80e-18	1.25e+04	6.48e-06	2.88e-04	5.88e-15
T2 flat near corner	1.89e-08	2.75e-18	1.43e+03	1.01e-05	3.66e-06	7.51e-15

The constant-pressure vector becomes a null vector to machine precision, the gauge is
pinned, and the 2-eps reproducibility becomes machine level — while the straight box is
bit-identical either way (all facets equal, so the weighting is a no-op there).

Caveats before anyone treats this as the patch:

• the algebra above is 2-D and P2-specific; 3-D triangular faces have their own
  ∫ φ weights and need the same derivation and measurement;
• it changes what the imposed normal means at a node (length-weighted rather than
  bisector), so boundary_normal_traction / dynamic_topography recovery must be
  re-validated against the existing benchmarks;
• it only touches the geometric normal path. A boundary given an analytic normal
  (normal=X/|X|, which is what the annulus and spherical shells use) bypasses this code
  entirely and is not fixed by it — see §7.

---

7. The annulus is not clean for the reason we assumed

Since a uniform annulus is clean because its facets are equal, the obvious question is
what happens when a curved rotated boundary is not uniform — which is what adaptation
and grading produce. Test: keep both annulus boundaries exactly circular (so the
analytic radial normal stays exact) and move the nodes to non-uniform angles,
θ → θ + δ sin 3θ. Only the facet lengths change.

The annulus admits a rigid rotation as well as the constant pressure, so two vectors are
attached and the null space is genuinely 2-D; principal angles are between subspaces of
dimension 2.

config	sv_1/smax	sv_2/smax	|A·(const-p)|/|A|	principal angles(attached, true null)	mean(p)	2-eps dv
E uniform	5.78e-17	4.43e-16	5.69e-14	2.30e-07, 9.97e-08 deg	-4.98e-04	6.34e-12
E skew 0.01	4.44e-16	1.68e-12	3.01e-07	1.61e-02, 3.05e-08 deg	-4.98e-04	2.28e-15
E skew 0.04	4.44e-16	2.52e-11	1.13e-06	6.50e-02, 7.43e-08 deg	-4.96e-04	1.49e-14

The same defect appears on the annulus as soon as it stops being uniform, and much
larger than in the box: the constant-pressure vector's non-nullity goes from 5.7e-14 to
1.1e-6 (2e7x) and the subspace misalignment from 2.3e-07 to 6.5e-02 degrees (3e5x). The
analytic normal does not protect it — it is exactly tangent to the true circle, but the
assembler integrates over the straight facets.

It nevertheless solves cleanly. The reason is visible in the singular values: the skew
also lifts the second singular value out of the machine-zero range (4.4e-16 → 2.5e-11),
so the mode is pinned — weakly, but enough that the solution component along it stays
O(10)·(RHS component 3.6e-10) instead of blowing up. In box configuration B the operator
stays exactly singular (sv_min = 4.4e-16) while the unremoved RHS component is 1.7e-9,
and 1.7e-9 / 4.4e-16 ≈ 4e6 — which, times |b| = 0.028, is the observed mean(p) of ~1e5.
That arithmetic reproduces the symptom to an order of magnitude and is the real severity
criterion:

the observable failure needs the operator to stay exactly singular in a direction the
attached nullspace does not span. The misalignment is the defect; exact singularity is
what converts it into a 1e4 pressure offset instead of a rounding error.

So the previous comment's "annulus is fine" should be narrowed: a uniform annulus is
fine, a skewed one carries the defect but does not (at these skews) express it. We did not
find a curved-boundary case where it does surface, so this is a flag on the blast-radius
assessment rather than a contradiction of it. The next experiment, if that matters, is a
genuinely adapted / strongly graded annulus or spherical shell, checking |A z|/|A| and
sv_min together rather than mean(p) alone.

---

8. Which mechanism the evidence supports, and what a fix must change

Not the constraint application at corners. The corner rows are built and eliminated
identically in every configuration (rank 2, both rows constrained, nothing dropped or
overwritten), and configuration T2 — whose corners are bit-identical to the clean straight
box — is the worst case measured.

Not the curvature of a rotated wall either. A curved wall with equal facets is clean
(§5, and the uniform annulus); a wall with unequal facets and no kink is clean (§5). It
takes both.

It is the nullspace attachment, in content rather than dimension. The constrained
rotated operator's null space has the right dimension and is nearly constant pressure;
UW3 attaches exactly constant pressure. Where every rotated facet is equal-length or
axis-aligned those coincide to 1e-18. Otherwise they differ by sin(angle) 2e-7 to 5e-6
(box) or up to 1e-3 (skewed annulus), PETSc projects out the wrong direction, the RHS keeps
an inconsistent component of 1e-9 to 4e-8, and the amplitude along the true null direction
is set by round-off. Whether that surfaces as a 1e4 pressure offset depends on the second
condition in §7 — the operator must remain exactly singular in that direction — which is
why the same defect is invisible on the skewed annulus and catastrophic in box B.

A fix has to make those two agree. Two families:

• Fix the operator — the real repair. Build the nodal frame so the tangential row is
  orthogonal to the facet-weighted normal the assembler actually integrates against,
  rather than to the bisector. §6 measures this: length/area-weighting the accumulation in
  _boundary_velocity_nodes restores |A z|/|A| to 2.7e-18 in every deformed
  configuration, pins mean(p) at 1e-5 or better, and leaves the straight box bit-identical.
  It is a one-line weighting change, but it needs the 3-D derivation, and it does not
  touch the analytic-normal path (§7).
• Fix the nullspace — the guard. Attach the operator's actual null vector instead of
  the analytic constant (one inverse iteration against the assembled Ahat starting from
  the constant-pressure vector), or pin the gauge with a datum row rather than a nullspace
  whenever the rotated boundary set is not flat. This also covers the analytic-normal path
  and any future rotated frame whose tangential rows are not exactly divergence-orthogonal.

Both are worth having: the first removes the consistency error at its source, the second
makes the solve safe even when some frame is imperfect. And whatever _rotated_nullspace
attaches, it is worth asserting that it is a null vector — |A z| / |A| is one matvec
and would have turned this into an immediate, loud failure instead of a silent 1e4 offset.

Measure any candidate on |A z|/|A| and |b.l_true|/|b| together with sv_min, not
on mean(p) or the velocity move — §4 and §7 show those two are chaotic and can report a
"fix" that is luck (or hide a real defect that happens not to be expressed).

---

Reproduction

Probes are self-contained and run unchanged against origin/development (ffbcf98), serial:

probe	section
probe1_taper.py	§1 tapered corner vs A and B
probe2_nullspace.py	§2 dense SVD, principal angles
probe3_decompose.py	§2 RHS inconsistency, difference decomposition
probe4_sweep.py	§4 amplitude and resolution sweeps
probe5_census.py	§3 corner census and residual map
probe6_equalarc.py	§5 equal-arc-length top
probe7_areaweight.py	§6 length-weighted nodal normal
probe8b_annulus.py	§7 skewed annulus

All instrumentation is a monkeypatch inside the probe; src/ is unmodified. One caution
for anyone re-running: Mat.convert("dense") on the captured operator corrupts enough
state to break the GAMG setup that follows — extract with getValuesCSR() instead.

Underworld development team with AI support from Claude Code

[lmoresi]

Fixed by PR #561, merged to development as bdb5671.

The nodal normal is now weighted by the facet measure the assembly actually integrates over — derived from the assembly (the free tangential row collects Σ_f (∫_f φ_i ds)(t̂_i·n_f), which vanishes only when the nodal normal is parallel to Σ_f (∫_f φ_i ds) n_f; for simplicial P1/P2 that integral is the measure times a facet-independent constant), not from the probe. |A z|/|A| collapses to ~6.6e-20 in every configuration, serial and parallel.

The parallel half turned out to matter more than the original defect: the weighted sum was rank-local, and the uniform annulus at np=4 was the worst case of all (1.40e-05 → 6.72e-20) — the configuration everyone treats as clean. The seam defect does not need a kink, it needs a partition. Three further defects were found and fixed underneath it: a rank-local outward test that let seam facets cancel, ranks owning nodes whose labelled facets all lived on neighbours, and a missing getSupportSize == 1 guard that would have let internal-boundary facets cancel.

fault_contact carried a verbatim copy and was more exposed than a wall — a constant pressure cancels in the mean rows but doubles in the jump rows, and the jump-tangential row is the slip row, so the bisector produced pressure-driven spurious slip scaling with lithostatic pressure, invisible to a leak diagnostic that contracts with the same normal that defined the constraint.

Convention change recorded: inner-boundary σ_nn now carries the sign the docstrings always claimed (outward = out of the domain). No golden value moved — every curved-boundary golden already passes an explicit normal=, and the CMB ones already pass -unit_r. Crameri Case-2 is unchanged to ten significant figures.

Not exact for non-affine 3-D hexes (leading-order, still 3e8–2e9 better than the bisector) — stated in the docstring. Two stale copies of the old rule carry TODO(BUG): boundary_flux._node_normals (dead) and FreeSurface._normal_direction (unreachable today).