Fix Stokes velocity FE dual-space setup for high-order simplex solves - #84
Merged
lmoresi merged 4 commits intoMar 19, 2026
Conversation
The Stokes velocity FE fix (previous commit) is mirrored here for the mesh coordinate projection FE in discretisation_mesh.py, which had the same omission. This only affects higher-order (curved) coordinate meshes on simplices — P1 meshes work by accident with PETSc defaults. Added comments to both sites explaining why these options are required. Underworld development team with AI support from Claude Code
lmoresi
approved these changes
Mar 19, 2026
lmoresi
left a comment
Member
There was a problem hiding this comment.
Reviewed the fix and audited all PETSc.FE().createDefault() call sites across the codebase.
Stokes velocity fix (Tyagi's original commit): Correct and minimal. The Stokes velocity FE was the only solver field missing the dual-space options — Scalar solver, Vector solver, and the Stokes pressure FE all had them already.
Added in follow-up commit:
- Comment explaining why these options matter (node placement on simplices, P2 works by accident, P3+ breaks)
- Same fix applied to the mesh coordinate projection FE in
discretisation_mesh.py— same omission, currently harmless for P1 coordinate meshes but would break higher-order (curved) coordinate meshes on simplices
All other FE creation sites (MeshVariable, gradient evaluation, field projection) already set the dual-space options correctly.
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
This fixes a Stokes FE setup inconsistency that breaks higher-order triangular velocity solves.
In
SNES_Stokes_SaddlePt._setup_discretisation(), the private velocityPetscFEwas created with the requested polynomial degree, but without the matching PETSc dual-space options that UW3 already uses for mesh variables and the other solver setup paths.This PR adds the missing velocity FE options before
PETSc.FE().createDefault(...):private_<prefix>_u_petscdualspace_lagrange_continuity = self.u.continuousprivate_<prefix>_u_petscdualspace_lagrange_node_endpoints = FalseReproducer
A minimal reproducer script is included in this PR:
examples/stokes_box_mms_simplex.pyRun it directly from the repo root, for example:
The script uses a simple manufactured Stokes solution on the unit square with triangular elements:
psi = x^2 y^2u = (2 x^2 y, -2 x y^2)p = x^2 - 1/3This is a clean test because the exact velocity is cubic and the exact pressure is quadratic, so
P3/P2should represent the solution essentially exactly on a straight box.Error norms before and after
Using the same box MMS setup at
h = 0.125:Before this fix:
P2/P1: velocity1.5068686519735996e-04, pressure3.9062499996873525e-03P3/P2: velocity1.2277197041312973e-02, pressure2.0668703128596483e-10After this fix:
P2/P1: velocity1.5068686519735996e-04, pressure3.9062499996873525e-03P3/P2: velocity3.6635897127143294e-12, pressure1.796728310502028e-10So the fix leaves
P2/P1unchanged and restoresP3/P2to the expected near-machine-precision result on this exact polynomial test.Notes
This PR is intentionally minimal:
It does not include the separate imported-mesh / XDMF work that came up during debugging.