Constrained free-slip via recoverable Lagrange multiplier (dynamic topography)

[lmoresi]

Summary

Adds uw.systems.Stokes_Constrained (SNES_Stokes_Constrained): enforce u·n = g
on curved boundaries with a recoverable Lagrange multiplier instead of a tuned
penalty. The converged multiplier is the normal traction holding the boundary — a
direct dynamic-topography estimate h = λ/(Δρ g).

Motivated by the fragility of penalty/Nitsche free-slip on annulus/spherical
boundaries (the magnitude must be tuned against forcing strength and viscosity; too
weak gives spurious radial throughflow, too strong diverges).

How it works

An augmented-Lagrangian (ALG2) outer loop reusing the existing penalty BC: each
Stokes solve carries t = [λ + r(u·n − g)]·n, and λ ← λ + r(u·n − g). The penalty
term preconditions every boundary mode (fast convergence); the multiplier removes the
penalty's accuracy bias (exact constraint from a moderate, well-conditioned r).

Purely additive — the validated 2×2 saddle-point assembly and fieldsplit
configuration are untouched (honours "solver stability is paramount"). The outer loop
lives inside solve().

Hands-off solve() (no tuning)

• Relative tolerance: RMS(u·n − g) < rtol·RMS|u| (default rtol=1e-3) — problem-independent.
• Viscosity-weighted augmentation: default r = 10³·μ(x), so high-viscosity
  boundary regions are constrained as well as low-viscosity ones (a flat r is not).

stokes=uw.systems.Stokes_Constrained(mesh, velocityField=v, pressureField=p)
stokes.constitutive_model=uw.constitutive_models.ViscousFlowModelstokes.constitutive_model.Parameters.shear_viscosity_0=mulam=stokes.add_constraint_bc("Upper", g=0.0) # free-slip, no penalty coefficientstokes.solve()
topo=stokes.topography("Upper", buoyancy_scale=delta_rho_g)

Validation

• Annulus (tests/test_1061_constrained_freeslip_annulus.py, 5 tests, level_2/tier_b):
  RMS(u·n)=6.5e-5, velocity within 0.3% of the penalty reference, 2 outer iterations,
  interior λ exactly 0, boundary corr(λ, −n·σ·n) = 0.9999.
• SolCx (variable viscosity): ~2e-4 vs Dirichlet reference in 3 outer iterations across
  contrast 1 → 10⁶, no manual tuning.

Trade-offs and roadmap

The augmentation r is a speed knob, not an accuracy knob (accuracy is
r-independent — the AL property): too small just costs iterations, too large stiffens
the inner solve, the answer is never wrong. True-work is U-shaped with an efficient
basin r ∈ [300, 10⁴].

Deferred (see docs/developer/design/CONSTRAINED_FREESLIP_MULTIPLIER.md):
adaptive r₀ for very coarse meshes; parallel (mask + nodal update are serial);
rotation-null-space removal for both-boundaries free-slip; and the monolithic
P'=[p,h] fieldsplit ("arbitrary constraints in the saddle point") as a separate
feasibility spike.

Design note: docs/developer/design/CONSTRAINED_FREESLIP_MULTIPLIER.md.

Underworld development team with AI support from Claude Code

[Copilot]

Pull request overview

Adds a new Stokes solver variant that enforces curved-boundary normal-velocity constraints (u·n=g) using an augmented-Lagrangian outer loop with a recoverable Lagrange multiplier, enabling direct access to the converged boundary traction (dynamic-topography proxy) without penalty tuning.

Changes:

• Introduces SNES_Stokes_Constrained (exported as uw.systems.Stokes_Constrained) with add_constraint_bc(), multiplier(), and topography() APIs.
• Implements the augmented-Lagrangian outer iteration inside solve() and a boundary-only multiplier update strategy via a nodal mask.
• Adds annulus validation tests and a developer design note documenting formulation, trade-offs, and validation results.

Reviewed changes

Copilot reviewed 4 out of 4 changed files in this pull request and generated 3 comments.

File	Description
src/underworld3/systems/solvers.py	Adds constrained Stokes solver, constraint BC registration, multiplier/topography accessors, and augmented-Lagrangian outer loop.
src/underworld3/systems/__init__.py	Exports the new solver as Stokes_Constrained.
tests/test_1061_constrained_freeslip_annulus.py	Adds regression/validation coverage comparing constrained vs penalty free-slip on an annulus and checks multiplier/topography behavior.
docs/developer/design/CONSTRAINED_FREESLIP_MULTIPLIER.md	Design/validation documentation for the constrained free-slip multiplier approach and roadmap.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[lmoresi]

Superseded by #224. The constrained free-slip work is replaced by the in-saddle (single coupled solve)Stokes_Constrained solver, which is more production-ready: one mesh-independent solve (no fragile outer loop), validated against the exact SolCx analytic solution, with a lossless boundary-only multiplier reduction. The augmented-Lagrangian outer-loop variant introduced here is removed in #224 (it's easy to reproduce in Python if wanted). Closing in favour of #224.