Consistent Jacobian tangent (unwrap-before-differentiate fix) — opt-in, default-off#258
Open
lmoresi wants to merge 5 commits into
Open
Consistent Jacobian tangent (unwrap-before-differentiate fix) — opt-in, default-off#258lmoresi wants to merge 5 commits into
lmoresi wants to merge 5 commits into
Conversation
Gated behind consistent_jacobian (default False = bit-identical Picard). - symbolic_keep_constants unwrap mode (constants[]-safe, drift-guarded) - consistent_jacobian: False | True | 'continuation' (alpha-blend via constants[]) - model-owned flux_jacobian hook; Nitsche bd_F1 wired Validated: level_1 tier_a 225/225, units 64/64, constant-visc bit-identical. Underworld development team with AI support from Claude Code
Creating the continuation alpha UWexpression in every solver __init__ bumped the global unique-name counter, shifting JIT cache keys enough to flip two known-flaky VEP variable-dt yield-lock tests in full-suite runs. Construct alpha lazily (continuation mode only) so the default Picard path creates no extra expression. Verified: level_2 tier_a forked failure set now IDENTICAL to pristine origin/development (4 pre-existing failures). Underworld development team with AI support from Claude Code
…n-unwrap-to-constants # Conflicts: # src/underworld3/cython/petsc_generic_snes_solvers.pyx
…owned) ViscoElasticPlasticFlowModel.flux_jacobian returns the flux with Min->harmonic (1/(1/eta_ve+1/eta_pl)) and Max smoothed, for the Jacobian source ONLY (exact Min residual preserved). Pure symbolic substitution on the live flux — no state mutation (fixes the earlier err-77 from a yield_mode toggle hack). Runs clean. Guard _newton_flux on consistent_jacobian so the default (Picard) path never evaluates flux_jacobian — keeps default assembly allocation-free / bit-identical (test_1010 6/6). NOTE: smooth-Jacobian + hard-Min residual is an inconsistent tangent and does NOT improve convergence on the BDF-2 VEP loading test (8/15 vs Picard 3/15); only full harmonic (consistent smooth residual+Jacobian) converges (0/15). The hook is correct, clean, and available; VEP convergence remains a separate issue. Underworld development team with AI support from Claude Code
The smooth-Jacobian-with-Min-residual tangent is inconsistent (consistent with the harmonic problem, not Min) and converges worse than Picard on hard-yield VEP — so the VEP-specific harmonic flux_jacobian override is deferred to the yield-law / delta-homotopy follow-up PR. The generic Constitutive_Model. flux_jacobian hook (default None) + the _newton_flux guard remain. Add docs/developer/design/jacobian-consistent-tangent.md: the bug, the opt-in/default-off fix, non-regression evidence, the tangent hierarchy, and the delta-homotopy successor work. Underworld development team with AI support from Claude Code
Contributor
There was a problem hiding this comment.
Pull request overview
This PR introduces an opt-in “consistent Jacobian tangent” path for nonlinear/viscoplastic SNES solves by unwrapping UWexpression atoms before symbolic differentiation (while still preserving truly-constant atoms for the PETSc constants[] mechanism). It also adds a continuation mode to blend Picard→Newton without triggering JIT recompiles, plus a constitutive-model hook to provide a custom Jacobian-only flux.
Changes:
- Add
symbolic_keep_constantsunwrapping mode to expand non-constantUWexpressionatoms while preserving truly-constant atoms as the same symbols forconstants[]. - Add
solver.consistent_jacobianwith modesFalse(default Picard),True(Newton), and"continuation"(alpha-blended Picard→Newton). - Add
Constitutive_Model.flux_jacobianhook (defaultNone) so models can supply a Jacobian-only surrogate flux.
Reviewed changes
Copilot reviewed 4 out of 4 changed files in this pull request and generated 9 comments.
| File | Description |
|---|---|
src/underworld3/function/expressions.py |
Adds symbolic_keep_constants unwrapping mode used to expose coefficient dependence during Jacobian differentiation without breaking constants[]. |
src/underworld3/cython/petsc_generic_snes_solvers.pyx |
Implements Jacobian source selection (Picard/Newton/continuation), uses unwrap-before-differentiate for Jacobian assembly, and adds continuation solve control flow. |
src/underworld3/constitutive_models.py |
Introduces flux_jacobian optional hook for Jacobian-only tangent substitution and documents intended usage. |
docs/developer/design/jacobian-consistent-tangent.md |
New design note explaining the unwrap-before-differentiate bug, the opt-in fix, and the tangent/continuation rationale. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
Comment on lines
+45
to
+47
| if isinstance(expr, sympy.NDimArray): | ||
| return sympy.Array([f(e) for e in expr], expr.shape) | ||
| return f(expr) # scalar expression |
| This is a no-op for constant-viscosity problems (eta has no grad-v | ||
| dependence), so those Jacobians stay bit-identical. | ||
|
|
||
| See ``docs/developer/design/jacobian-unwrap-constants-bug.md``. |
| # the model's own smooth law (constitutive_model.flux_jacobian) when it | ||
| # provides one; otherwise the exact unwrapped flux. | ||
| # | ||
| # See docs/developer/design/jacobian-unwrap-constants-bug.md. |
| # The residual fns above (self._u_F0/_u_F1/_p_F0) are left untouched — | ||
| # getext() unwraps those itself. For constant-viscosity problems this | ||
| # is a no-op (eta has no grad-v dependence) so the Jacobian is | ||
| # bit-identical. See docs/developer/design/jacobian-unwrap-constants-bug.md |
Comment on lines
+1069
to
+1071
| # Restore a clean Picard tangent for any subsequent solve (next step). | ||
| self._set_newton_alpha(0.0) | ||
|
|
Comment on lines
+1024
to
1029
| if self.consistent_jacobian == "continuation": | ||
| self._continuation_solve(gvec, verbose=verbose) | ||
| else: | ||
| self.snes.solve(None, gvec) | ||
| if divergence_retries <= 0: | ||
| return |
Comment on lines
+580
to
+582
| Returns ``None`` by default, meaning the solver differentiates the | ||
| exact :attr:`flux` (the Newton fix unwraps it first; a generic Min/Max | ||
| kink-smoothing fallback then rounds any remaining yield kink). |
Comment on lines
+1718
to
+1719
| # (ramp the softmin softness δ→0), not a smooth tangent. See the design doc | ||
| # docs/developer/design/jacobian-unwrap-constants-bug.md. The generic |
Comment on lines
+2486
to
+2489
| # JACOBIAN: unwrap (keep constants) + smooth Min/Max kinks so the | ||
| # derivative sees the field-dependence of any nonlinear coefficient | ||
| # (full Newton) while the residual keeps the exact form. No-op for | ||
| # constant coefficients -> bit-identical. See _jacobian_source. |
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.
Consistent Jacobian tangent for nonlinear (viscoplastic) solves — opt-in, default-off
The bug
SNES Jacobian assembly differentiated the residual flux
F1while the effectiveviscosity was still a wrapped
UWexpressionatom, so∂η/∂(grad v)was silently droppedfrom every Jacobian. UW3 viscoplastic Stokes was therefore running an accidental Picard /
defect-correction tangent, not full Newton — the unwrap happened after the derivative
instead of before it. Constant-viscosity problems were unaffected, which is why it stayed
hidden behind the "≈20 Picard iterations is intrinsic" folklore.
The fix (default-off, bit-identical by default)
symbolic_keep_constantsunwrap mode — expandsUWexpressionatoms down to (but notincluding) truly-constant atoms (η₀, τ_y stay symbolic for the
constants[]mechanism).The keep-symbolic predicate is the same
_is_truly_constantused byconstants[]extraction, so they cannot drift (drift-guard test).
solver.consistent_jacobian:False(default, frozen/Picard, bit-identical) /True(full Newton) /"continuation"(Picard→Newton via anα-blend routed throughconstants[], so switching costs no JIT recompile andα=0is bit-identical).Constitutive_Model.flux_jacobianhook (defaultNone) for a model to supply a smoothtangent law;
_newton_fluxis guarded so the default path never evaluates it.The residual never goes through the new path, so a converged solution always satisfies the
exact constitutive law.
Non-regression evidence
test_1010, SolCxtest_1015,test_0610, asymmetric-Jacobian guard, units (64) pass;level_1 tier_a225/225.--forked)level_2 tier_afailure set is identical to pristinedevelopment— the pre-existing reds (test_1012gmsh crash + 3test_1052VEP) are notintroduced here.
Scope / why a smooth Jacobian is not the VEP fix
For a hard-
Minresidual, a smooth Jacobian is the consistent tangent of a different(harmonic) problem — it diverges more than Picard. The robust route for hard-yield VEP is
problem-space homotopy (ramp the softmin softness δ→0; δ=0 is identically
Min), which isa separate follow-up PR that reuses this PR's
constants[]-ramp machinery. Seedocs/developer/design/jacobian-consistent-tangent.mdfor the full tangent hierarchy and thesuccessor plan.
Underworld development team with AI support from Claude Code