Summary
SNES_Stokes(_SaddlePt).add_nitsche_bc (and the SNES_Vector version) produces a boundary Jacobian that is inconsistent with its residual when the constitutive model is anisotropic (e.g. TransverseIsotropicFlowModel). The residual is correct, so solves still reach the right answer — but Newton degrades to a slow defect-correction (tens of oscillating SNES steps instead of 1–2), because the linearisation it follows is wrong by a few percent. The bulk Jacobian, derive_by_array, and the JIT are all correct; the error is confined to the Nitsche boundary terms for anisotropic stress.
Reproduce
snes_test_jacobian on a minimal TI Stokes (tilted director, η₁ ≠ η₀), varying only the boundary condition / geometry:
| configuration |
‖J − J_fd‖/‖J‖ |
| box + TI |
3e-11 ✓ |
| annulus + Dirichlet + TI |
2e-9 ✓ |
| annulus + field-varying viscosity + TI |
3e-11 ✓ |
annulus + penalty free-slip (add_natural_bc) + TI |
6.7e-11 ✓ |
| annulus + Nitsche + isotropic |
2e-9 ✓ |
| annulus + Nitsche + TI |
0.069 (7%) ✗ |
The failure is precisely the Nitsche × anisotropy intersection. Amat and Pmat are equally wrong, and the ratio is invariant to pressure-nullspace on/off, monolithic-lu vs fieldsplit-Schur — so it is the tangent the BC contributes, not solver configuration.
A self-contained reproducer is scripts/jac_minimal_box.py (on the feature/fault-convection worktree); it sweeps GEOM/BC/VISC/NOWEAK and runs snes_test_jacobian. Happy to upstream a trimmed version as a regression test.
Mechanism (likely a small fix)
-
Consistency term → f0 gradient dependence. Nitsche injects the gradient-dependent traction σ(u)·n into the f0 (value) boundary residual, so the boundary Jacobian g1 = ∂f0/∂L is nonzero and is exercised for the first time anywhere — the bulk f0 is just body force (∂f0/∂L = 0), so this assembly path is never tested in the bulk. The boundary assembly (cython/petsc_generic_snes_solvers.pyx ~L5665–5686) reuses the bulk (0,2,1,3) index permutation/reshape, which is not correct for the anisotropic boundary traction.
-
Symmetry term uses scalar viscosity. The Nitsche symmetry term hardcodes mu = constitutive_model.viscosity (= shear_viscosity_0) rather than the full C : ε̇(v) response (add_nitsche_bc, vector ~L2864 / saddle ~L4557). Correct for isotropic, wrong for TI.
Impact
Latent for every anisotropic model under Nitsche — TI viscous, TI-VEP, anisotropic diffusion. For the fault-convection work it presented as "TI Stokes is intractable" (≈23 oscillating Newton steps); it is in fact a regression introduced when free-slip was switched from penalty (add_natural_bc) to add_nitsche_bc (penalty never exercises the bad path).
Workaround
Use penalty free-slip add_natural_bc(KFS · v·n̂ · n̂) instead of add_nitsche_bc for anisotropic rheologies — verified to give an exact Jacobian (6.7e-11).
Acceptance / regression gate
snes_test_jacobian on annulus + Nitsche + TI goes from 0.069 to ~1e-8, with isotropic + Nitsche unchanged. Add scripts/jac_minimal_box.py (trimmed) to the test suite.
Underworld development team with AI support from Claude Code
Summary
SNES_Stokes(_SaddlePt).add_nitsche_bc(and theSNES_Vectorversion) produces a boundary Jacobian that is inconsistent with its residual when the constitutive model is anisotropic (e.g.TransverseIsotropicFlowModel). The residual is correct, so solves still reach the right answer — but Newton degrades to a slow defect-correction (tens of oscillating SNES steps instead of 1–2), because the linearisation it follows is wrong by a few percent. The bulk Jacobian,derive_by_array, and the JIT are all correct; the error is confined to the Nitsche boundary terms for anisotropic stress.Reproduce
snes_test_jacobianon a minimal TI Stokes (tilted director,η₁ ≠ η₀), varying only the boundary condition / geometry:‖J − J_fd‖/‖J‖3e-11✓2e-9✓3e-11✓add_natural_bc) + TI6.7e-11✓2e-9✓0.069(7%) ✗The failure is precisely the Nitsche × anisotropy intersection.
AmatandPmatare equally wrong, and the ratio is invariant to pressure-nullspace on/off, monolithic-luvs fieldsplit-Schur — so it is the tangent the BC contributes, not solver configuration.A self-contained reproducer is
scripts/jac_minimal_box.py(on thefeature/fault-convectionworktree); it sweepsGEOM/BC/VISC/NOWEAKand runssnes_test_jacobian. Happy to upstream a trimmed version as a regression test.Mechanism (likely a small fix)
Consistency term →
f0gradient dependence. Nitsche injects the gradient-dependent tractionσ(u)·ninto thef0(value) boundary residual, so the boundary Jacobiang1 = ∂f0/∂Lis nonzero and is exercised for the first time anywhere — the bulkf0is just body force (∂f0/∂L = 0), so this assembly path is never tested in the bulk. The boundary assembly (cython/petsc_generic_snes_solvers.pyx~L5665–5686) reuses the bulk(0,2,1,3)index permutation/reshape, which is not correct for the anisotropic boundary traction.Symmetry term uses scalar viscosity. The Nitsche symmetry term hardcodes
mu = constitutive_model.viscosity(=shear_viscosity_0) rather than the fullC : ε̇(v)response (add_nitsche_bc, vector ~L2864 / saddle ~L4557). Correct for isotropic, wrong for TI.Impact
Latent for every anisotropic model under Nitsche — TI viscous, TI-VEP, anisotropic diffusion. For the fault-convection work it presented as "TI Stokes is intractable" (≈23 oscillating Newton steps); it is in fact a regression introduced when free-slip was switched from penalty (
add_natural_bc) toadd_nitsche_bc(penalty never exercises the bad path).Workaround
Use penalty free-slip
add_natural_bc(KFS · v·n̂ · n̂)instead ofadd_nitsche_bcfor anisotropic rheologies — verified to give an exact Jacobian (6.7e-11).Acceptance / regression gate
snes_test_jacobianon annulus + Nitsche + TI goes from0.069to~1e-8, with isotropic + Nitsche unchanged. Addscripts/jac_minimal_box.py(trimmed) to the test suite.Underworld development team with AI support from Claude Code