Stokes_Constrained: consolidated state of the open work (supersedes #244, #248, #254, #495)

[lmoresi]

Consolidating the Stokes_Constrained backlog. Six issues (#244, #248, #254, #495, #607,
#614) had accumulated on one subsystem, mostly from a single benchmark campaign, with
overlapping symptoms and no single place showing what is actually still broken. This is
that place. #244, #248, #254 and #495 are closed into it; #607 and #614 stay open because
each has distinct, actionable content (and #607 has a PR).

Status re-measured 2026-08-26, not copied forward.

Fixed — verified, not assumed

The [p, λ] gauge is partition-consistent now. #254's item 1: the combined
constant-(pressure, multiplier) gauge had a partition-dependent level — mean pressure
3.333 serial against 2.983 at np=8, a 10.5% difference that corrupted topography, since
topography is read from the multiplier. There is now a dedicated regression,
test_1063_constrained_freeslip_parallel.py::test_constrained_raw_gauge_partition_independent,
and it passes.

The 2-D partition-dependence is gone. #495 reported
test_constrained_freeslip_partition_independent failing at np=2 (iso −3.4%) and np=4
(ti +2.8%). Re-run today at np=1, 2 and 4: all pass. The question #495 raised —
whether that assertion should be a cross-build invariant rather than an absolute golden —
is worth keeping in mind if it recurs, but there is nothing failing to act on.

The Zhong velocity catastrophe (#248) was a SphericalShellInternalBoundary
Lower/Internal label bug, fixed by #242.

Still open

1. ~0.4% velocity partition-dependence in 3-D constraint/augmentation assembly.
#254's item 2, and the one to be careful about: it is 3-D-specific, and #254 states
plainly that 2-D test_1063 was bit-identical to 1e-9 even when this was first reported.
So today's np=1/2/4 passes above do not clear it — they are the 2-D case, which was
never the symptom. With the gauge pinned, velocity L2 / |Ut| / |Ub| still differed
0.36 / 0.37 / 0.41% serial against np=8, both sides fully converged, which points at the
assembled operator or RHS rather than a nullspace. Locus: the add_constraint_bc /
augmented-Lagrangian (augmentation_base=1e4) boundary assembly at inter-rank seams.
Needs a 3-D parallel regression; there is no test covering this today.

2. Interior-multiplier knockout is lossy, and the docstring says it is lossless.
#254's item 3. _constrain_interior_multipliers_in_section (gated by
_reduce_interior_multiplier, default True) pins interior multiplier DOFs. Disabling it
changes the serial answer by 2% (L2, |Ut|) to 5% (|Ub|), so the pinned DOFs are not
inert. The DOF count is partition-identical (5170 at np=1/2/4/8), so it is not gross loss.
Either the boundary-trace closure drops DOFs that matter, or "inert interior multiplier"
is the wrong model — decide which, then fix or re-document.

3. Monolithic direct solve is known-bad. A monolithic MUMPS LU on the constrained
system gives wrong physics (|Ut| = 4e-3 against a validated 1e-2) and segfaults at
np=8. The solver already warns that it is a serial diagnostic only. Capture as
known-bad in a test, or fix.

4. Slow in parallel, and not because of Krylov. #244: the 8-rank runtime is much worse
than the equivalent Nitsche free-slip solve while the constrained solve converges in two
KSP iterations, so the cost is setup, assembly or field handling for the block-constrained
system. Reproduced on SphericalShellInternalBoundary (ri 0.55, rint 0.775, ro 1.0,
cellSize 0.125, P2/P1, np=8). Never profiled to a line.

5. Topography from the constrained path is not validated against Zhong. The remaining
half of #248. The velocity half is fixed; the diagnostic half — constrained free slip
provides multiplier/topography diagnostics, Nitsche matches the Table 2 velocities but its
raw boundary sigma_rr is not enough for topography recovery — was never closed out.
Note this now interacts with #607: if multiplier() is missing the augmented-Lagrangian
share, a topography comparison built on it is measuring the wrong quantity.

Kept open separately

• Stokes_Constrained: multiplier() is not the whole boundary traction — the augmented-Lagrangian share is missing #607 — multiplier() is not the whole boundary traction; the momentum row is
  (h + r(n·u − g))·n, so the traction is h + r(n·u − g) and only h is returned.
  PR Return the whole boundary traction from the constrained solver (#607) #617 addresses it. Item 5 above depends on this landing.
• boundary_flux() returns ~1e12 on Stokes_Constrained #614 — boundary_flux() returns ~1e12 on Stokes_Constrained while the solve itself
  is fine (velocity error 8.8e-06, recovered n·σ·n good to 0.075). Suspected
  _assemble_volume_reaction() picking up multiplier rows. Distinct primitive bug; not
  re-verified today.

Method note, carried from #254 because it is worth not relearning

Point evaluation (uw.function.evaluate) is unreliable for serial-versus-parallel field
comparison here — it gave 134%-of-|v| pointwise noise while the volume L2 differed 0.4%.
Use partition-independent integrals only. (See also #641, where evaluate returns a
negative value for a provably non-negative expression.)

Ruled out and not worth redoing: velocity rigid-rotation nullspace (angular momentum
reproducible to ~1e-5 serial against np=8); iterative under-convergence (serial and np=8
each bit-identical at tol 1e-7 against 1e-11 — both converge fully, to different answers).

Reproduction: ~/+Simulations/repro_248_internal_load.py, args
<method> <write|check> <solver> <tol> <knockout on|off> <gauge none|pin>; amr-dev env,
parallel needs --with-mpi.

Suggested order

1. Stokes_Constrained: multiplier() is not the whole boundary traction — the augmented-Lagrangian share is missing #607 / PR Return the whole boundary traction from the constrained solver (#607) #617 — it is in flight and item 5 depends on it.
2. A 3-D parallel regression for item 1. There is no test today, so the defect could
   silently resolve or worsen and nobody would know.
3. boundary_flux() returns ~1e12 on Stokes_Constrained #614 — self-contained and likely a small fix.
4. Items 2 and 3 — decide the model, then fix or document.
5. Stokes_Constrained is slow in parallel despite low KSP iterations #244 — profile before optimising.