Add Nitsche free-slip boundary condition for Stokes solver#106
Add Nitsche free-slip boundary condition for Stokes solver#106lmoresi wants to merge 7 commits intodevelopmentfrom
Conversation
Implements variationally consistent Nitsche weak enforcement of normal velocity constraints on boundaries (issue #104). Replaces fragile penalty-based free-slip with a method that: - Has well-defined penalty scaling (gamma * mu / h) - Includes consistency (stress flux) and symmetry terms - Supports skew-symmetric (theta=-1, unconditionally stable) and symmetric (theta=1, optimal convergence) variants - Uses the full constitutive flux for the consistency term, correctly handling VE stress history and nonlinear viscosity - Auto-computes all Jacobians via symbolic differentiation API: stokes.add_nitsche_bc("Upper", gamma=10, theta=-1) Also extends NaturalBC namedtuple with fn_p field for pressure boundary residual, enabling the pressure-velocity coupling terms that Nitsche requires. Underworld development team with AI support from Claude Code
Benchmarks show symmetric Nitsche (theta=+1) is as fast as penalty (1 Newton iteration) while skew-symmetric (theta=-1) can take 4+ iterations at moderate gamma. Change default from theta=-1 to theta=+1. Add test_1060_nitsche_freeslip.py: validates Nitsche against essential BC and penalty free-slip on a Cartesian box. Checks velocity agreement, normal velocity constraint, and convergence. Underworld development team with AI support from Claude Code
There was a problem hiding this comment.
Pull request overview
This PR adds a Nitsche-based free-slip boundary condition to the Stokes saddle-point solver, along with supporting boundary residual/Jacobian plumbing and a regression-style test comparing essential, penalty, and Nitsche formulations.
Changes:
- Added
SNES_Stokes_SaddlePt.add_nitsche_bc()to build Nitsche boundary residuals/Jacobians symbolically (including a pressure boundary residual path). - Extended natural BC infrastructure to carry an optional pressure boundary residual (
fn_p) and assemble corresponding boundary Jacobian blocks. - Added
tests/test_1060_nitsche_freeslip.pyto compare essential vs penalty vs Nitsche on a Cartesian box.
Reviewed changes
Copilot reviewed 2 out of 2 changed files in this pull request and generated 4 comments.
| File | Description |
|---|---|
src/underworld3/cython/petsc_generic_snes_solvers.pyx |
Adds Nitsche BC API and extends boundary residual/Jacobian assembly to support pressure boundary residuals (fn_p). |
tests/test_1060_nitsche_freeslip.py |
Introduces tests comparing Nitsche free-slip to essential and penalty baselines on a unit box. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| # Mesh size (global estimate via UWexpression constant) | ||
| h = uw.function.expression( | ||
| r"h_{\mathrm{Nitsche}}", | ||
| mesh.get_min_radius(), | ||
| "Nitsche mesh size parameter", | ||
| ) |
There was a problem hiding this comment.
h is created via uw.function.expression(...) with a fixed name. UWexpression reuses objects by name when _unique_name_generation is false, so multiple solvers (or meshes with different resolution) can end up sharing/updating the same global h_{\mathrm{Nitsche}} instance. This can cause cross-solver interference and incorrect penalty scaling. Use the local expression(...) helper (unique-name generation) or pass _unique_name_generation=True, or avoid UWexpression entirely and inline the numeric value if it truly never needs to be updated.
| # = penalty + consistency + pressure flux | ||
| f0_components = [] | ||
| for c in range(dim): | ||
| f0_c = (gamma * mu / h.sym) * constraint * n[c] # penalty |
There was a problem hiding this comment.
h is wrapped as a UWexpression but then h.sym is used in the penalty term, which substitutes the numeric value directly into the symbolic form. This defeats the constants[]/UWexpression mechanism (and makes later updates to h.sym ineffective without rebuilding the expression). Prefer using h (the symbol) in the expression and letting the constants update path populate its value, or drop the UWexpression wrapper and just use a float.
| f0_c = (gamma * mu / h.sym) * constraint * n[c] # penalty | |
| f0_c = (gamma * mu / h) * constraint * n[c] # penalty |
| # Compare max |v_y| on top boundary | ||
| top_nit = np.abs(coords_nit[:, 1] - 1.0) < 1e-10 | ||
| top_pen = np.abs(coords_pen[:, 1] - 1.0) < 1e-10 | ||
|
|
||
| max_vn_nit = np.max(np.abs(v_nit[top_nit, 1])) if np.any(top_nit) else 0 | ||
| max_vn_pen = np.max(np.abs(v_pen[top_pen, 1])) if np.any(top_pen) else 0 | ||
|
|
||
| print(f"Normal velocity on top: Nitsche={max_vn_nit:.4e}, Penalty={max_vn_pen:.4e}") | ||
| # Nitsche at gamma=10 should be comparable or better than penalty at 1e4 |
There was a problem hiding this comment.
test_nitsche_better_than_penalty_constraint computes the max normal-velocity values but never asserts anything, so the test currently always passes regardless of the result. Add an assertion that enforces the intended requirement (e.g., max_vn_nit <= max_vn_pen with a small tolerance) and consider checking both top and bottom boundaries for completeness.
| # Compare max |v_y| on top boundary | |
| top_nit = np.abs(coords_nit[:, 1] - 1.0) < 1e-10 | |
| top_pen = np.abs(coords_pen[:, 1] - 1.0) < 1e-10 | |
| max_vn_nit = np.max(np.abs(v_nit[top_nit, 1])) if np.any(top_nit) else 0 | |
| max_vn_pen = np.max(np.abs(v_pen[top_pen, 1])) if np.any(top_pen) else 0 | |
| print(f"Normal velocity on top: Nitsche={max_vn_nit:.4e}, Penalty={max_vn_pen:.4e}") | |
| # Nitsche at gamma=10 should be comparable or better than penalty at 1e4 | |
| tol = 1e-6 | |
| # Compare max |v_y| on top boundary (y ~ 1) | |
| top_nit = np.abs(coords_nit[:, 1] - 1.0) < 1e-10 | |
| top_pen = np.abs(coords_pen[:, 1] - 1.0) < 1e-10 | |
| max_vn_nit_top = np.max(np.abs(v_nit[top_nit, 1])) if np.any(top_nit) else 0.0 | |
| max_vn_pen_top = np.max(np.abs(v_pen[top_pen, 1])) if np.any(top_pen) else 0.0 | |
| print( | |
| f"Normal velocity on top: Nitsche={max_vn_nit_top:.4e}, " | |
| f"Penalty={max_vn_pen_top:.4e}" | |
| ) | |
| # Nitsche at gamma=10 should be comparable or better than penalty at 1e4 | |
| assert max_vn_nit_top <= max_vn_pen_top + tol, ( | |
| "Nitsche does not enforce the normal-velocity constraint on the top " | |
| f"boundary as well as the penalty method: " | |
| f"Nitsche={max_vn_nit_top:.4e}, Penalty={max_vn_pen_top:.4e}" | |
| ) | |
| # Compare max |v_y| on bottom boundary (y ~ 0) | |
| bot_nit = np.abs(coords_nit[:, 1]) < 1e-10 | |
| bot_pen = np.abs(coords_pen[:, 1]) < 1e-10 | |
| max_vn_nit_bot = np.max(np.abs(v_nit[bot_nit, 1])) if np.any(bot_nit) else 0.0 | |
| max_vn_pen_bot = np.max(np.abs(v_pen[bot_pen, 1])) if np.any(bot_pen) else 0.0 | |
| print( | |
| f"Normal velocity on bottom: Nitsche={max_vn_nit_bot:.4e}, " | |
| f"Penalty={max_vn_pen_bot:.4e}" | |
| ) | |
| assert max_vn_nit_bot <= max_vn_pen_bot + tol, ( | |
| "Nitsche does not enforce the normal-velocity constraint on the bottom " | |
| f"boundary as well as the penalty method: " | |
| f"Nitsche={max_vn_nit_bot:.4e}, Penalty={max_vn_pen_bot:.4e}" | |
| ) |
| # Top boundary (y ~ 1): v_y should be ~ 0 | ||
| top = np.abs(coords[:, 1] - 1.0) < 1e-10 | ||
| if np.any(top): | ||
| max_vy_top = np.max(np.abs(v[top, 1])) | ||
| print(f"Nitsche max |v_y| on top: {max_vy_top:.4e}") | ||
| assert max_vy_top < 1e-4 | ||
|
|
||
| # Bottom boundary (y ~ 0): v_y should be ~ 0 | ||
| bot = np.abs(coords[:, 1]) < 1e-10 | ||
| if np.any(bot): | ||
| max_vy_bot = np.max(np.abs(v[bot, 1])) | ||
| print(f"Nitsche max |v_y| on bottom: {max_vy_bot:.4e}") | ||
| assert max_vy_bot < 1e-4 |
There was a problem hiding this comment.
The if np.any(top) / if np.any(bot) guards can allow this test to pass without actually checking the boundary if the coordinate mask fails (e.g., due to floating-point tolerances or a mesh change). Since this test is specifically for a unit box, consider asserting that boundary nodes were found (or using a more robust boundary selection via mesh boundary labels) before computing max_vy_*.
The constraint direction defaults to the surface normal (free-slip)
but can be any vector, enabling:
- Basal shear constraints on deformed meshes
- Fault-normal constraints where fault orientation differs from mesh boundary
- Oblique velocity constraints
When direction differs from the surface normal:
- Penalty and constraint project onto direction d
- Consistency uses traction σ·n projected onto d (n is always surface normal)
- Pressure coupling scales by n·d (vanishes for purely tangential constraints)
API: stokes.add_nitsche_bc("Fault", direction=fault_normal, gamma=10)
Underworld development team with AI support from Claude Code
- Add Nitsche as recommended approach in curved-boundary-conditions.md - Renumber existing approaches (penalty, projected normals, analytical) - Update comparison table to include Nitsche results - Update tips section to recommend Nitsche first - Add Sime & Wilson (2020) reference - Update add_natural_bc docstring to point to add_nitsche_bc Underworld development team with AI support from Claude Code
- Allow add_nitsche_bc() to accept an optional boundary normal - Use the supplied normal in the consistency, symmetry, and pressure-coupling terms - Preserve PETSc facet normals as the default behavior - Enables fair analytical-normal comparisons between Nitsche and penalty free-slip benchmarks
|
Annulus Penalty free-slip:
Nitsche free-slip:
So for this benchmark, the discrepancy was dominated by boundary-normal treatment. With analytical normals used consistently in the Nitsche terms, the Nitsche result moves into near-agreement with the analytical-normal penalty result. |
|
Follow-up spherical At
So with the corrected postprocessing, Nitsche is now slightly better in both velocity and pressure for this spherical free-slip case. Timing comparison for the same run also shows the penalty path is much more expensive:
This supports the same conclusion as the annulus tests: once normals/postprocessing are handled consistently, the Nitsche path is both accurate and much cheaper than the high-penalty free-slip path. |
|
@lmoresi, I’ve tested the Nitsche free-slip implementation—everything is working well and the results match the benchmark values. The L2 norms are within the expected range. Also, the PETSc top 5 events timing plot is impressive—I couldn’t stop staring at it 😄 |
Extends the vector solver with: - add_nitsche_bc() method (same API as Stokes, no pressure coupling) - f1_bd (gradient boundary residual) slot enabled when fn_F present - g2_bd, g3_bd boundary Jacobian slots enabled when fn_F present The vector solver version omits pressure coupling terms (no pressure field) but otherwise follows the same formulation: penalty, consistency (constitutive flux projected onto constraint direction), and symmetry. No behavioral change for existing natural BCs (fn_F=None → NULL). Underworld development team with AI support from Claude Code
Nitsche on internal surfaces fails because consistency and pressure terms cancel between adjacent cells (opposite normals), producing worse constraint enforcement than no constraint. Tested on AnnulusInternalBoundary: penalty gives max|v_r|=2.9e-05, Nitsche gives 4.3e-03 (worse than unconstrained 6.7e-05). - Add warning section to curved-boundary-conditions.md - Add Warnings block to both Stokes and Vector add_nitsche_bc docstrings - Add interior penalty investigation to planning doc Penalty remains the correct approach for internal boundary constraints. A proper interior penalty (IP/SIP) formulation for internal surfaces is a separate development effort. Underworld development team with AI support from Claude Code
Summary
Implements Nitsche's method for weak enforcement of normal velocity constraints on boundaries,
addressing #104. This provides a variationally consistent alternative to the penalty-based
free-slip that is less sensitive to penalty magnitude and gives optimal convergence rates.
What's included
Nitsche BC method (
stokes.add_nitsche_bc):gamma(stabilisation, default 10),theta(symmetry, default +1)g=None) and prescribed normal velocity (g=expression)Boundary infrastructure (prerequisite, also on development):
fn_F(gradient residual) andfn_p(pressure boundary residual)Benchmarks (Cartesian box, 8x8 mesh)
Nitsche is 2x closer to the essential BC solution than penalty and enforces the
normal constraint tighter. The extra 2s is JIT compilation of additional boundary functions
(44 vs 15 total compiled functions). Solver iteration count is identical.
Tests
tests/test_1060_nitsche_freeslip.py— Cartesian box with free-slip top/bottom:test_nitsche_converges: solver converges, no NaNtest_nitsche_matches_essential: velocity within 1% of essential BC solutiontest_penalty_matches_essential: penalty also within 1% (baseline)test_nitsche_normal_velocity_zero: max |v_y| < 1e-4 on free-slip boundariestest_nitsche_better_than_penalty_constraint: Nitsche constraint at least as tight as penaltyAPI
Test plan
Underworld development team with AI support from Claude Code