Return the whole boundary traction from the constrained solver (#607)

docs/advanced/curved-boundary-conditions.md

@@ -71,29 +71,48 @@ stokes.add_nitsche_bc(0.0, "Fault", direction=fault_normal, gamma=10)
 at 1e4 gives 0.15%).
 
 
-### 2. Penalty Free-Slip (Simple but Fragile)
-
-Use the mesh-derived normals directly with a penalty parameter:
+### 2. Penalty Free-Slip (Simple, and Only With the Node Normal)
 
 ```python
-Gamma = mesh.Gamma
+n = mesh.boundary_normal("Boundary")     # measure-weighted node normal
 penalty = 10000
-stokes.add_natural_bc(penalty * Gamma.dot(v.sym) * Gamma, "Boundary")
+stokes.add_natural_bc(penalty * n.dot(v.sym) * n, "Boundary")
 ```
 
+**Use `mesh.boundary_normal(boundary)`, not `mesh.Gamma`.** A penalty written
+against the per-facet normal asks a node shared by two facets to satisfy two
+different constraints, which on a two-component velocity leaves nothing: push
+the coefficient up and the boundary freezes. Measured on an annulus with an
+exact solution (Kramer et al. 2021), coefficient `1e6`, cell 0.15 → 0.035:
+
+| normal | leak `u·n` | velocity error | surface stress error |
+|---|---|---|---|
+| `mesh.Gamma` (facet) | 1e-5 | 0.60, flat under refinement | 0.21 → 0.26, growing |
+| `mesh.boundary_normal` (node) | 3e-5 | 1.0e-2 → 4.9e-4 | 2.4e-2 → 1.4e-3 |
+
+The facet-normal row does not converge, and the leak cannot see it: at `1e3` that
+penalty leaks 3e-2 and gets the surface stress right to 2e-3, while at `1e8` it
+leaks 1e-7 and is 26% wrong. This is the classical over-constraint that the
+*consistent* normal was introduced to avoid (Engelman, Sani & Gresho 1982), and
+`mesh.boundary_normal` is that normal.
+
 **When to use:**
 - Quick prototyping where high accuracy isn't critical
 - When Nitsche is not yet available for your solver type
 
 **Limitations:**
 - Penalty must be tuned: too small → loose constraint, too large → ill-conditioning
 - On spherical shells, penalty can become unstable at moderate resolution
-- ~25-30% error on elliptical boundaries when using raw facet normals
+- Check the answer, not just the leak: a constraint that is satisfied is not
+  evidence that the solution is right
 
 
-### 3. Projected Normals (For Curved Boundaries with Penalty)
+### 3. Projected Normals (Superseded by `mesh.boundary_normal`)
 
-Project `mesh.Gamma` onto a continuous mesh variable, which interpolates and smooths the normals:
+`mesh.boundary_normal(boundary)` assembles the measure-weighted node normal
+directly, tracks mesh deformation, and is correct in parallel, so the recipe
+below is kept for reference rather than recommended. Project `mesh.Gamma` onto a
+continuous mesh variable, which interpolates and smooths the normals:
 
 ```python
 import sympy

docs/developer/CHANGELOG.md

@@ -6,6 +6,40 @@ This log tracks significant development work at a conceptual level, suitable for
 
 ## 2026 Q3 (July – September)
 
+### The Multiplier Was Not the Whole Traction (August 2026)
+
+**`Stokes_Constrained.topography()` now returns the traction the boundary is
+actually held with**, and a new `traction()` exposes it directly. The momentum
+row carries `λ + r(n·u − g)`, so the bare multiplier is short by the
+augmented-Lagrangian share — `r` times the discrete constraint residual. With the
+viscosity-weighted default `r = 1e4·μ(x)` that share is a few per cent of the
+surface traction on a uniform-viscosity annulus and most of it across a `1e6`
+viscosity step, where `λ` alone reads a tenth of the exact SolCx topography and
+is anti-correlated with it. `multiplier()` still returns `λ` and now says what it
+is not.
+
+The defect survived because the validation scored a **correlation** (0.9999)
+between the multiplier and the recovered normal stress. A correlation is
+scale-free and cannot see a systematic amplitude deficit, which is precisely what
+a missing share of the load is. The new guard,
+`tests/test_1063_constrained_traction.py`, scores a relative `l2` against the
+exact SolCx surface topography and carries the bare multiplier as its negative
+control.
+
+The corrected quantity is the consistent boundary flux: at convergence
+`M_Γ(λ + r(n·u − g))` balances the volume residual restricted to the boundary,
+which is the CBF nodal load (Zhong, Gurnis & Hulbert 1993). So the multiplier
+route and the rotated constraint's `boundary_normal_traction` are the same
+computation, and they agree to 3–5% — inside each route's own error against the
+exact answer.
+
+Documentation: `docs/advanced/curved-boundary-conditions.md` now writes the
+penalty free-slip recipe against `mesh.boundary_normal` rather than `mesh.Gamma`.
+A penalty against the per-facet normal over-constrains the shared nodes and does
+not converge — measured on an annulus at coefficient `1e6`, the velocity error
+stays at 0.60 and the surface-stress error grows from 0.21 to 0.26 as the mesh is
+refined, while the leak reads 1e-5 throughout. (underworld3#607, #608, #614)
+
 ### The Free Surface Reaches the Spherical Shell (July 2026)
 
 **`uw.systems.FreeSurface` now runs in 3D on a spherical shell** — the same

docs/developer/design/CONSTRAINED_FREESLIP_MULTIPLIER.md

@@ -2,8 +2,24 @@
 
 **Status**: shipped as `uw.systems.Stokes_Constrained` (serial). The constraint
 is enforced by a multiplier carried **inside** the saddle point (one coupled
-solve); the converged boundary multiplier is the normal traction = dynamic
-topography. An earlier augmented-Lagrangian **outer-loop** variant was removed in
+solve); the converged boundary traction is `λ + r(n·u − g)`, which is the dynamic
+topography, and it is returned by `traction()` / `topography()`.
+
+> **Correction, 2026-08-19.** This document said in several places that `λ` alone
+> is the normal traction and that the augmentation `r` is a pure speed knob. Both
+> are true only in the exact limit. Discretely the constraint row is satisfied to
+> the solver's tolerance and `r` multiplies that residual back into the traction,
+> so `λ` is short by `r(n·u − g)`. With the viscosity-weighted default
+> `r = 1e4·μ(x)` the omitted share is a few per cent of the surface traction on a
+> uniform-viscosity annulus, and across SolCx's `1e6` viscosity step it is most of
+> it — `λ` alone reads a tenth of the exact topography and is *anti-correlated*
+> with it. `traction()` and `topography()` now return the sum; `multiplier()`
+> still returns `λ`. See underworld3#607 and
+> `tests/test_1063_constrained_traction.py`.
+>
+> Why the existing validation did not catch it: it scored `corr(λ, −n·σ·n) ≈
+> 0.9999`. A correlation is scale-free and cannot see a systematic amplitude
+> deficit, which is exactly what a missing share of the load is. An earlier augmented-Lagrangian **outer-loop** variant was removed in
 favour of this in-saddle formulation (it is straightforward to reproduce in
 Python if needed). Validated against the exact SolCx analytic solution
 (`tests/test_1062_constrained_solcx.py`).
@@ -19,9 +35,9 @@ is ~100× too weak at Ra=1e6); too strong and the system ill-conditions and the
 Stokes solve diverges in line search.
 
 This feature enforces `u·n = g` on a curved boundary with a **true Lagrange
-multiplier** `λ` instead of a penalty. Because the converged multiplier *is* the
-normal traction holding the boundary, it is simultaneously a direct estimate of
-**dynamic surface topography**, `h = λ / (Δρ g)`. The equilibrium `λ` is also the
+multiplier** `λ` instead of a penalty. Because the converged boundary term *is*
+the normal traction holding the boundary, it is simultaneously a direct estimate
+of **dynamic surface topography**, `h = (λ + r(n·u − g)) / (Δρ g)`. The equilibrium `λ` is also the
 target end-state toward which a free surface can be integrated over a time
 interval (connecting to the ETD free-surface work on
 `feature/exp-integrator-freesurface`).
@@ -104,9 +120,10 @@ topo = stokes.topography("Upper", buoyancy_scale=delta_rho_g)   # h = lambda/(dr
 ```
 
 `solve()` does **one coupled solve** — no outer iteration or constraint tuning.
-The augmentation defaults to `1e4·μ(x)` (local-viscosity-weighted); accuracy is
-independent of it (the λ-row carries the exact constraint), so no per-problem
-tuning is needed.
+The augmentation defaults to `1e4·μ(x)` (local-viscosity-weighted). The
+CONSTRAINT is independent of it (the λ-row carries the exact constraint), so no
+per-problem tuning is needed, and so is the traction read through `traction()` /
+`topography()`. The bare `λ` is NOT: it is short by `r(n·u − g)`.
 
 Key design points:
 
@@ -139,15 +156,27 @@ Two regression tests cover the shipped solver:
   **exact analytic** solution: velocity `rel ≈ 8.7e-6` (== the Dirichlet
   baseline), constraint `RMS(u·n) ≈ 1.6e-10`.
 
-The consistent-boundary-flux identity `λ = −n·σ·n|_Γ` is the independent
-cross-check: the multiplier's boundary trace equals the recovered normal Cauchy
-stress (negative sign = the reaction traction holding the boundary), confirming
-`λ` is the dynamic topography signal.
+The consistent-boundary-flux identity is the independent cross-check, and it is
+an identity for the WHOLE boundary term rather than for `λ`:
+
+```
+M_Γ (λ + r(n·u − g)) = −(A·u − b)|_Γ        (the CBF nodal load)
+```
+
+so `λ + r(n·u − g)` is the CBF traction de-smeared by the boundary mass — the
+same computation `boundary_normal_traction` performs on a rotated constraint,
+arrived at by carrying the traction as an unknown instead of reading it out of
+the residual. Measured on SolCx at a `1e6` contrast, the corrected multiplier
+gives relative `l2` 0.047 against the exact topography where `λ` alone gives 1.04.
+
+A third test, `tests/test_1063_constrained_traction.py`, guards this with the
+bare multiplier as its negative control.
 
 ## The augmentation parameter `r`: true-work trade-off
 
-`r` is a *speed* knob, not an *accuracy* knob — this is the key advantage over a
-pure penalty, and it carries over to the in-saddle solver (accuracy is
+`r` is a *speed* knob for the CONSTRAINT, not an accuracy knob — this is the key
+advantage over a pure penalty, and it carries over to the in-saddle solver
+(accuracy is
 `r`-independent; `r` only sets the iteration count). The sweep table below is from
 the **historical outer-loop** variant (its "outer iterations" have no analogue in
 the one-shot coupled solve), but the shape and the conclusion stand. For the

src/underworld3/cython/petsc_generic_snes_solvers.pyx

@@ -8110,8 +8110,10 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
             # u-row residual:  fn_f = h·n  +  r(n·u − g)·n
             # The r-term is the augmented-Lagrangian penalty: it adds a uu
             # boundary stiffness r·(n⊗n) that conditions the Schur complement
-            # but does NOT bias the multiplier (the h-row stays the exact
-            # constraint, so h still converges to the true normal traction).
+            # but does not change what the constraint ENFORCES (the h-row stays
+            # the exact constraint). It does change what h IS: the traction is
+            # h + r(n.u - g), and only the sum is r-independent. Stokes_Constrained
+            # .traction() / .topography() return that sum; .multiplier() returns h.
             fn_f = sympy.Matrix(
                 [(hsym + r_sym * (u_dot_n - g_sym)) * n_row[i] for i in range(dim)]
             ).as_immutable()