Stokes_Constrained: in-saddle multiplier free-slip + topography (supersedes #219)#224
Merged
Merged
Conversation
…graphy Enforce u.n = g on curved boundaries with a recoverable Lagrange multiplier instead of a tuned penalty. An augmented-Lagrangian (ALG2) outer loop carries both the existing penalty BC and the multiplier; the multiplier removes the penalty's accuracy bias, so a moderate, well-conditioned augmentation gives an exact constraint in 2-3 Stokes solves. The converged multiplier is the normal traction holding the boundary -- a direct dynamic-topography estimate (h = lambda / (drho g)), recoverable as a clean MeshVariable (interior exactly zero; boundary trace correlates with -n.sigma.n at 1.0000). Purely additive: the validated 2x2 saddle-point assembly and fieldsplit config are untouched. New class behind uw.systems.Stokes_Constrained. Hands-off solve(): relative constraint tolerance (RMS(u.n-g) < rtol*RMS|u|) and viscosity-weighted augmentation (1e3*mu(x)) by default -- no user tuning. On the SolCx benchmark this gives ~2e-4 accuracy in 3 outer iterations across viscosity contrast 1 -> 1e6. Validation: tests/test_1061_constrained_freeslip_annulus.py (5 tests). Design note and production roadmap in docs/developer/design/. Underworld development team with AI support from Claude Code
Confirmed at the PETSc level that a 3-field DM (u, p, h) with p and h grouped yields a 2-way u | [p,h] Schur split (no nested fieldsplit), de-risking the monolithic "arbitrary saddle-point constraints" direction. Captures the bounded assembly touch-points and remaining caveats for a future block-Schur PR. Underworld development team with AI support from Claude Code
- add_constraint_bc validates the boundary name and the 'UW_Boundaries' label, and rejects empty/None point sets (petsc_dm_find_labeled_points_local returns the vertex-0 sentinel np.array([0]) when the label is absent, which an `is not None` check would silently accept and corrupt the mask). - The outer loop no longer applies a final, never-solved multiplier update when the iteration cap is hit (u/p stay consistent with lambda) and now warns loudly (RuntimeWarning) on non-convergence. - Raise NotImplementedError when run on more than one MPI rank (the mask and the node-wise update are serial-only), instead of silently producing wrong results. Adds a test that add_constraint_bc rejects an unknown boundary name. Underworld development team with AI support from Claude Code
…+ topography Enforces u.n = g on a boundary via a Lagrange multiplier h carried INSIDE the saddle-point system (grouped u | [p,h] 2-way Schur split), in one coupled solve. The converged boundary trace h|_G = -n.sigma.n is dynamic topography directly, with no penalty/Nitsche parameter to tune. - Guarded generalisation of SNES_Stokes_SaddlePt: every multiplier block wrapped `if self._multipliers:` so the 2-field path stays bit-identical (M1 gate). - Augmented-Lagrangian conditioning of the constraint Schur (r = 1e3.mu, viscosity-weighted); combined (p,h) gauge nullspace for enclosed domains. - FMG-ready: field-index fieldsplit grouping (pc_fieldsplit_N_fields) + option mirror, so geometric multigrid works on the velocity block. - Boundary-only multiplier reduction (pin interior h) DEFAULT OFF: refined-safe via DMPlexLabelComplete closure but not yet lossless on the clone DM (degrades boundary trace ~5e-4); the lossless PetscSection-constraint path is next. - add_constraint_bc(boundary, g, normal, screening, augmentation, degree); multiplier()/topography() recovery API. - tests/test_1062: 11 cases (box + annulus; constraint / Dirichlet-match / topography / variable-viscosity to 1e6), all passing. Key finding: on hard constraints (SolCx, 4-wall free-slip + viscosity jump) the block does one mesh-independent solve where the outer-loop Uzawa needs 9-11 and explodes to 33x Dirichlet -- the two methods fail in opposite places (block: per-solve cost from the fat [p,h] Schur; outer-loop: outer-iteration count). Underworld development team with AI support from Claude Code
Constrain the interior (off-boundary) multiplier DOFs directly in the fine local PetscSection rather than via DMAddBoundary. This is both lossless (the constraint-boundary closure is correct because createClosureIndex has finalised the section) and refined/FMG-safe (no "after section creation" restriction), so the solved [p,h] Schur block carries only the boundary trace (~sqrt(ndof) DOFs). Collapses the multiplier-block DOF premium from ~3x to ~1.1x Dirichlet with no accuracy loss (box+lu reduced-vs-full velocity relL2 ~1e-9). Underworld development team with AI support from Claude Code
…iant The in-saddle (single coupled solve) constrained free-slip solver is now the production path, exported as uw.systems.Stokes_Constrained. It is validated against the exact SolCx analytic solution and matches a Dirichlet free-slip reference to discretisation error, with the constraint enforced to machine precision. - rename SNES_Stokes_BlockConstrained -> SNES_Stokes_Constrained - remove the augmented-Lagrangian outer-loop solver + helper (_ConstraintBC); the Uzawa scheme is easy to reproduce in Python if wanted - single public export Stokes_Constrained (drop Stokes_BlockConstrained) - raise default augmentation_base 1e3 -> 1e4 (accuracy is r-independent; larger sits in the low-iteration plateau, well below roundoff) - tests: retarget the constrained test, add test_1062_constrained_solcx.py (free-slip vs the SolCx analytic); outer-loop annulus test removed - docs: update CONSTRAINED_FREESLIP_MULTIPLIER.md status Underworld development team with AI support from Claude Code
Contributor
There was a problem hiding this comment.
Pull request overview
This PR promotes an in-saddle-point constrained free-slip Stokes solver as uw.systems.Stokes_Constrained, introducing coupled Lagrange-multiplier boundary constraints (including a boundary-only multiplier reduction) and adding regression tests and a design note to document/validate the approach.
Changes:
- Add
SNES_Stokes_Constrainedwithadd_constraint_bc(),multiplier(), andtopography()APIs for enforcingu·n=gvia in-saddle multipliers. - Extend the PETSc/Cython Stokes saddle-point assembly to support multiplier fields, boundary coupling, nullspace handling, and boundary-only multiplier DOF reduction.
- Add new level_2 tests validating constrained free-slip (box + annulus) and matching SolCx analytic; add/update design documentation.
Reviewed changes
Copilot reviewed 6 out of 6 changed files in this pull request and generated 10 comments.
Show a summary per file
| File | Description |
|---|---|
src/underworld3/systems/solvers.py |
Adds the public constrained Stokes solver wrapper and constraint BC registration/topography API. |
src/underworld3/cython/petsc_generic_snes_solvers.pyx |
Implements the core multiplier field registration, boundary residual/Jacobian coupling, nullspace, fieldsplit grouping, and section-based DOF reduction. |
src/underworld3/systems/__init__.py |
Exports SNES_Stokes_Constrained as uw.systems.Stokes_Constrained. |
tests/test_1061_constrained_freeslip.py |
Adds constrained free-slip regression tests (box + annulus) and multiplier/topography API checks. |
tests/test_1062_constrained_solcx.py |
Adds SolCx analytic validation test for constrained free-slip on all walls. |
docs/developer/design/CONSTRAINED_FREESLIP_MULTIPLIER.md |
Adds a design note for constrained free-slip and dynamic-topography interpretation (currently contains several outdated outer-loop descriptions). |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| via DEFAULT solver options (no per-script PETSc options); compared against | ||
| the penalty reference, with topography recovery h ~ -n.sigma.n. | ||
|
|
||
| Run: pixi run python -m pytest tests/test_1062_block_constrained_freeslip.py -v |
Comment on lines
+1958
to
+1963
| # Boolean over the multiplier's nodes: True on INTERIOR nodes (off the | ||
| # constraint boundary). The constant-on-interior h mode is a near-null | ||
| # mode (interior h appears only in the screening eM), so we hand it to | ||
| # the solver as a nullspace vector -- the gauge-fixing analogue of the | ||
| # constant-pressure nullspace. Set in add_constraint_bc. | ||
| self.interior_mask = None |
Comment on lines
+2156
to
+2174
| # Interior mask: True on multiplier nodes OFF the constraint boundary. | ||
| # Build a P1 marker that is 1 on the boundary vertices and sample it at | ||
| # the multiplier's nodes (boundary mid-edge nodes interpolate to ~1). | ||
| from underworld3.discretisation.discretisation_mesh import ( | ||
| petsc_dm_find_labeled_points_local, | ||
| ) | ||
| marker = uw.discretisation.MeshVariable( | ||
| f"_bmask_{self.instance_number}_{idx}", self.mesh, 1, degree=1, | ||
| ) | ||
| marker.data[:] = 0.0 | ||
| if self.mesh.dm.hasLabel("UW_Boundaries"): | ||
| pts = petsc_dm_find_labeled_points_local( | ||
| self.mesh.dm, "UW_Boundaries", | ||
| getattr(self.mesh.boundaries, boundary).value, sectionIndex=False, | ||
| ) | ||
| if pts is not None and len(pts) > 0: | ||
| marker.data[pts] = 1.0 | ||
| on_boundary = np.array(uw.function.evaluate(marker.sym, h.coords)).reshape(-1) > 0.5 | ||
| cbc.interior_mask = ~on_boundary |
Comment on lines
+39
to
+43
| `C` is **co-dimension-1**: it touches only velocity DOFs on Γ. A monolithic | ||
| third FE field would therefore either waste interior DOFs or need a boundary | ||
| trace space PETSc/DMPlex does not provide on the same DM. We instead solve the | ||
| boundary Schur complement `S_λ = C K⁻¹ Cᵀ` by an **outer loop**, leaving the | ||
| validated 2×2 Stokes assembly untouched. |
Comment on lines
+45
to
+61
| ### Augmented Lagrangian (ALG2) — the production algorithm | ||
|
|
||
| Each Stokes solve carries a natural-BC traction with **both** the multiplier and | ||
| a penalty augmentation, reusing the existing penalty machinery: | ||
|
|
||
| $$\mathbf{t} = \bigl[\lambda + r\,(\mathbf{u}\cdot\mathbf{n} - g)\bigr]\,\mathbf{n} | ||
| \quad\text{on } \Gamma,$$ | ||
|
|
||
| and the multiplier is updated with the same `r`: | ||
|
|
||
| $$\lambda \leftarrow \lambda + r\,(\mathbf{u}\cdot\mathbf{n} - g).$$ | ||
|
|
||
| The penalty term `r(u·n)n` preconditions *every* boundary mode uniformly, so the | ||
| outer loop converges in a handful of iterations; the multiplier removes the | ||
| penalty's accuracy bias, so a **moderate, well-conditioned** `r` gives both fast | ||
| convergence *and* an exact constraint — unlike a pure penalty, which must be made | ||
| large and fragile. |
Comment on lines
+104
to
+108
| The outer loop is hidden behind ``solve()``: the tolerance is **relative** to the | ||
| velocity scale (``RMS(u.n-g) < rtol·RMS|u|``, default ``rtol=1e-3``) so it is | ||
| problem-independent, and the augmentation defaults to ``1e3·μ(x)`` (local- | ||
| viscosity-weighted). On SolCx this gives ~2e-4 accuracy in 3 outer iterations at | ||
| viscosity contrast 1 → 10⁶ with no user tuning. |
Comment on lines
+240
to
+242
| - `src/underworld3/systems/solvers.py` — `SNES_Stokes_Constrained`, `_ConstraintBC`. | ||
| - `src/underworld3/systems/__init__.py` — `Stokes_Constrained` export. | ||
| - `tests/test_1061_constrained_freeslip_annulus.py` — validation. |
| blk.add_dirichlet_bc((0.0, 0.0), "Lower") | ||
| hb = blk.add_constraint_bc("Upper", g=0.0, normal=unit_r) | ||
| # Enclosed -> pressure/multiplier gauge; DEFAULT grouped-Schur solver config. | ||
| blk._petsc_use_pressure_nullspace = True |
| s.add_constraint_bc("Right", g=0.0, normal=sympy.Matrix([[ 1.0, 0.0]])) | ||
| s.add_constraint_bc("Bottom", g=0.0, normal=sympy.Matrix([[ 0.0, -1.0]])) | ||
| s.add_constraint_bc("Top", g=0.0, normal=sympy.Matrix([[ 0.0, 1.0]])) | ||
| s._petsc_use_pressure_nullspace = True |
Comment on lines
+234
to
+236
| **Recommendation:** the architecture is de-risked, but since the outer loop converges in | ||
| 2–3 iterations with no user-facing tuning, monolithic is scheduled work (a general | ||
| saddle-point-constraint capability), not an urgent replacement. |
- remove the unused _BlockConstraintBC.interior_mask and its construction (a P1 marker field + evaluate) — dead since the section-based reduction; avoidable setup cost - tests: use the public petsc_use_pressure_nullspace property (its setter resets cached nullspace state); fix the stale run-command path in test_1061 - docs/CONSTRAINED_FREESLIP_MULTIPLIER.md: rewrite the outer-loop sections to the shipped in-saddle formulation (no outer iteration; r is AL stabilisation only; augmentation_base 1e4; boundary-only reduction); fix the Files list and the option-table verdicts; mark the r-sweep table as historical outer-loop data Underworld development team with AI support from Claude Code
Member
Author
|
Addressed the Copilot review (commit 9863744):
12/12 constrained tests still green. |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Makes the in-saddle (single coupled solve) constrained free-slip solver the production
uw.systems.Stokes_Constrained, and removes the augmented-Lagrangian outer-loop variant (PR #219). The constraintu·n = gis enforced by a Lagrange multiplier carried inside the saddle point; the converged boundary multiplier is the normal traction = dynamic topography (multiplier()/topography()).This supersedes #219 — the outer-loop solver is closer to half-formed and fragile on hard constraints (it explodes on viscosity jumps), whereas the in-saddle solver is one mesh-independent solve and is now validated against an exact analytic solution.
What's here
SNES_Stokes_Constrainedis the in-saddle solver; single public exportuw.systems.Stokes_Constrained(the outer-loop class and_ConstraintBCare removed — the Uzawa scheme is easy to reproduce in Python if wanted).[p,h]block DOF premium from ~3× to ~1.1× Dirichlet with no accuracy loss.augmentation_baseraised 1e3 → 1e4 (accuracy is independent of the AL value; larger sits in the low-iteration plateau, well below roundoff).tests/test_1062_constrained_solcx.py): free-slip via four in-saddle multipliers matches the analytic to8.7e-6— indistinguishable from a Dirichlet free-slip reference — with the constraint enforced to machine precision (RMS(u·n) ≈ 1.6e-10).Merge note
Please squash-merge — the branch was rebased onto
developmentcarrying #219's commits underneath, so squash collapses them into one clean commit. The net file diff is the in-saddle solver only (no outer-loop).Depends on #223 (SolCx analytic), already merged.
Underworld development team with AI support from Claude Code