Consistent Jacobian tangent (unwrap-before-differentiate fix) — opt-in, default-off by lmoresi · Pull Request #258 · underworldcode/underworld3

lmoresi · 2026-06-19T10:35:34Z

Consistent Jacobian tangent for nonlinear (viscoplastic) solves — opt-in, default-off

The bug

SNES Jacobian assembly differentiated the residual flux F1 while the effective
viscosity was still a wrapped UWexpression atom, so ∂η/∂(grad v) was silently dropped
from 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_constants unwrap mode — expands UWexpression atoms down to (but not
including) truly-constant atoms (η₀, τ_y stay symbolic for the constants[] mechanism).
The keep-symbolic predicate is the same _is_truly_constant used by constants[]
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 through
constants[], so switching costs no JIT recompile and α=0 is bit-identical).
Constitutive_Model.flux_jacobian hook (default None) for a model to supply a smooth
tangent law; _newton_flux is 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

Constant-viscosity Jacobian symbolically bit-identical (0/N blocks change).
test_1010, SolCx test_1015, test_0610, asymmetric-Jacobian guard, units (64) pass;
level_1 tier_a 225/225.
Crash-isolated (--forked) level_2 tier_a failure set is identical to pristine
development — the pre-existing reds (test_1012 gmsh crash + 3 test_1052 VEP) are not
introduced here.

Scope / why a smooth Jacobian is not the VEP fix

For a hard-Min residual, 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 is
a separate follow-up PR that reuses this PR's constants[]-ramp machinery. See
docs/developer/design/jacobian-consistent-tangent.md for the full tangent hierarchy and the
successor plan.

Underworld development team with AI support from Claude Code

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

Copilot

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_constants unwrapping mode to expand non-constant UWexpression atoms while preserving truly-constant atoms as the same symbols for constants[].
Add solver.consistent_jacobian with modes False (default Picard), True (Newton), and "continuation" (alpha-blended Picard→Newton).
Add Constitutive_Model.flux_jacobian hook (default None) 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.

+    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


+        # Restore a clean Picard tangent for any subsequent solve (next step).
+        self._set_newton_alpha(0.0)
+


+        if self.consistent_jacobian == "continuation":
+            self._continuation_solve(gvec, verbose=verbose)
+        else:
+            self.snes.solve(None, gvec)
        if divergence_retries <= 0:
            return


+        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).


+    # (ramp the softmin softness δ→0), not a smooth tangent. See the design doc
+    # docs/developer/design/jacobian-unwrap-constants-bug.md. The generic


+        # 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.


…n-unwrap-to-constants

bknight1 · 2026-06-25T14:47:01Z

@lmoresi I've had a look at this PR and it's vastly improved the convergence of the VP models I've been working on, both when using continuous and True. The only issue is compilation time with multimaterials, which can be >5 mins, with 'continuous' taking longer than True. Another issue is the compilation appears to occur twice after updating the constitutive model (one to check the hash is different then a rebuild the jacobians).

I have made some modifications to fix these issues, happy to discuss committing these changes.

…n-unwrap-to-constants # Conflicts: # src/underworld3/cython/petsc_generic_snes_solvers.pyx

lmoresi · 2026-07-02T08:41:53Z

Reviewed + rebased onto current development (resolved the petsc_generic_snes_solvers.pyx conflict — kept both this PR's Jacobian methods and development's set_custom_mg; builds clean, Stokes 6/6 passes).

Safe to merge: the consistent tangent is opt-in and default-off (consistent_jacobian = False), so the default assembly path is bit-identical to prior behaviour (green CI confirms no regression). The opt-in Newton benefit and the deferred harmonic-Min flux_jacobian are follow-up work under #239. Merging via admin (team-authored).

Underworld development team with AI support from Claude Code

…og, value-first call-site sweep (WE-01..03,05,06,08,09,10) (#338) * docs(WE-01): adopt the one-governing-doc-per-topic authority map Repoint CLAUDE.md's Data Access 'Authoritative Reference' from the stale UW3_Style_and_Patterns_Guide.md to subsystems/data-access.md (the guide it crowned teaches patterns the code deprecates at runtime — DOC-04), and record the Style Charter §10 authority table in docs/developer/index.md as the master authority index. The Charter is added to the Getting Started toctree (removes a baseline 'not included in any toctree' warning). Finding: DOC-04 (docs/reviews/2026-07/DOCS-STANDARDS-COHERENCE.md). Underworld development team with AI support from Claude Code * docs(WE-02): de-drift the Style Guide's four stale normative sections Rewrites the sections DOC-01 verified as contradicting the settled standards: - Docstring format: the 'Markdown Docstrings for pdoc/pdoc3' section is replaced by the NumPy/Sphinx RST standard (worked example with :math: and Parameters/Returns/Examples/Notes; conversion tracked in docs/plans/docstring-conversion-plan.md), per Style Charter section 6. - Doc file format: Quarto .qmd prescription (zero .qmd files exist in the repo) replaced by MyST .md/Sphinx guidance matching CLAUDE.md; migration table row updated. - Data access examples: 'Preferred' coordinate examples now use the real, runnable API — mesh.X.coords (read), mesh.deform() (coordinate changes), and the swarm.coords getter/setter for particle positions. The previous 'Preferred' example swarm.data += displacement raises AttributeError (getter-only property — SWARM-13 evidence); mesh.data warns at runtime. The private-attribute migration advice (swarm._particle_coordinates, mesh._deform_mesh presented as the NEW pattern) is deleted. - Front matter: the 21-line Quarto YAML header is replaced by a minimal MyST title block, and the guide now states that the UW3 Style Charter is the normative contract and wins on conflict. All replacement examples verified against current source: Swarm.coords setter (swarm.py), Mesh.deform (discretisation_mesh.py:3133), uw.synchronised_array_update / NDArray_With_Callback.delay_callbacks_global. Findings: DOC-01, SWARM-13 (style-guide part). Underworld development team with AI support from Claude Code * docs(WE-03): regenerate the docstring review queue; add the sweep to the release checklist The queue (last generated 2026-01-13, cdf5bb2) misrepresented the codebase both ways: it flagged now-complete items (solve, SNES_Scalar) as missing and contained zero entries for the June 2026 API (DOC-02). Regenerated over src/underworld3/**/*.py + **/*.pyx at the current tip. Two bugs in scripts/docstring_sweep.py's regex-based Cython parser made the regenerated queue lie about .pyx docstrings and are fixed as part of making the regeneration meaningful: - the indent group '(\s*)' with re.MULTILINE consumed preceding blank lines, shifting the computed definition line so the docstring search started ON the def/class line and always missed; - the docstring search started at the definition line rather than after the (possibly multi-line) signature, so long signatures hid their docstrings; - raw-string docstrings (r""", the norm in the solver .pyx) were not recognised. DOC-02 cross-validation on the regenerated queue now passes: solve / SNES_Scalar in the solver pyx are no longer flagged 'none'; the queue contains the June API (add_nitsche_bc, add_rotated_freeslip_bc, boundary_flux, set_custom_fmg, consistent_jacobian: 13 mentions) and flags the DOC-05 targets (Swarm.advection x2, read_timestep, write_proxy) as undocumented. Also adds the sweep to the quarterly release checklist (guides/release-process.md) so the queue cannot go stale unnoticed again. Findings: DOC-02 (docs/reviews/2026-07/DOCS-STANDARDS-COHERENCE.md). Underworld development team with AI support from Claude Code * docs(WE-05): backfill the changelog for May - early July 2026; add the changelog sweep to the release checklist The changelog (the quarterly CIG/stakeholder record) ended in April 2026 while ~117 first-parent commits landed May through early July (DOC-03). Backfilled at the existing conceptual granularity — 14 grouped entries, grouped by subsystem rather than by PR, matching the established format (### Title (Month Year), bold lead sentence, hyphen bullets, inline PR references): - New '2026 Q3 (July - September)' section: the July 2026 quality campaign (#309-#313, #317, #322-#326, #329, #334 as grouped entries), rotated strong free-slip / boundary traction / dynamic topography (#293, #294, #298, #306), generalized geometric multigrid via custom prolongation (#290, #297), consistent Jacobian tangent (#258), swarm correctness (#216, #313, #323, #329), numpy 2 support (#301, #305). - Extended '2026 Q2' section with the May-June entries: mesh adaptation movers (#190, #209, #213, #228, #259, #264, #266), moving-mesh field transfer / deform() (#246, #249, #251), semi-Lagrangian accuracy controls (#164, #183, #185-#189, #208, #220), snapshot/checkpoint toolkit (#146, #195, #196, #198), Stokes_Constrained (#224, #229, #240, #265), local-h Nitsche + boundary-slip surfaces (#225, #241, #275), units interoperability (#277, #278, #283, #284), memory/evaluation/solver infrastructure (#161, #177-#179, #181, #182, #222, #237, #250, ...). Every entry is backed by a merged commit on development (verified against git log --first-parent aed517f..3184a40). Also adds a quarterly-changelog sweep step beside the docstring sweep in the release checklist (guides/release-process.md) per DOC-03's proposed fix. Findings: DOC-03 (docs/reviews/2026-07/DOCS-STANDARDS-COHERENCE.md). Underworld development team with AI support from Claude Code * docs(WE-06): status headers on the unmarked design docs (per-doc git verification) Adds one-to-three-line Status markers to the 13 design docs that lacked one, following the directory's existing conventions (**Status**: line under the title; status: key inside existing YAML frontmatter for the three frontmatter-only docs), and corrects the stale 'Design Phase' marker on MATHEMATICAL_MIXIN_DESIGN.md (the mixin ships in utilities/mathematical_mixin.py). Every stamp was verified against git history (git log --follow dates) and the current source tree before writing: - Implemented: jacobian-consistent-tangent (PR #258, c63cd70), fmg-checkpoint-hierarchy (3cd73cd), petsc-dmplex-checkpoint-reload-plan (PR #146, write_timestep(petsc_reload=True) in tree), fault-refinement-simplification (smooth_mesh_interior / metric_density_from_gradient / fault_comb_metric all in tree), MATHEMATICAL_MIXIN_DESIGN. - Current reference/contract: mesh-adaptation-formulation, ND_UNITS_BOUNDARY_CONTRACT (PR #278, e0ece9a). - Investigation records (preserved via PR #245, 34a9dd4; production geometric-MG is custom prolongation, PR #290): snesfas-feasibility, snesfas-vanka-feasibility-study. - Design notes / prototypes with honest gaps: in_memory_checkpoint_design (not implemented, per its own trailing Status section), submesh-solver-architecture (extract_region/extract_surface exist; coarsened_companion does not). - Historical: ARCHITECTURE_ANALYSIS (persistence.py layout superseded), COORDINATE_MIGRATION_GUIDE (transition shipped), WHY_UNITS_NOT_DIMENSIONALITY (decision record). The audit's ~16 estimate over-counted: re-derived at this tip, 13 docs were unmarked plus one marked-but-stale (DOC-07). Findings: DOC-07 (docs/reviews/2026-07/DOCS-STANDARDS-COHERENCE.md). Underworld development team with AI support from Claude Code * docs(WE-08): convert units.py public docstrings Google -> NumPy style Docstring-only conversion of the 18 public module-level functions that carried Google-style Args:/Returns:/Raises:/Examples: labels (check_units_consistency, get_dimensionality, get_units, non_dimensionalise, show_nondimensional_form, simplify_units, create_quantity, convert_units, to_base_units, to_reduced_units, to_compact, get_scaling_coefficients, set_scaling_coefficients, validate_expression_units, assert_dimensionality, validate_coordinates_dimensionality, enforce_units_consistency, require_units_if_active, convert_angle_to_degrees) to the NumPy/Sphinx standard (Style Charter section 6). dimensionalise was already NumPy style; one-line docstrings and private helpers are untouched. No code, signature, or behaviour changes (verified: every diff hunk is inside a docstring; ast.parse clean). Finding: API-12 (docs/reviews/2026-07/API-CONSISTENCY-REVIEW.md). Underworld development team with AI support from Claude Code * docs(WE-09): sweep call sites of the newer BC methods to value-first (conds, boundary, ...) order Wave C (#334) made the ORIGINAL value-first order canonical for add_nitsche_bc / add_rotated_freeslip_bc / add_constraint_bc (maintainer decisions D2/D3; Style Charter section 6) with deprecation shims for the legacy boundary-first and g= spellings. This sweep updates every call site of those THREE methods to the canonical order so nothing in the repository exercises the shims — 74 sites total: - tests/: 63 call sites across 12 files (test_1017, test_1018, test_1060, test_1061, test_1062, test_1064, test_1065 x2 serial; parallel test_1017, test_1062, test_1063, test_1064). tests/test_0641_wave_c_api_shims.py is deliberately untouched — its legacy-order calls ARE the deprecation contract. - docs/: 7 sites (curved-boundary-conditions.md x4, CONSTRAINED_FREESLIP_MULTIPLIER.md call + signature line, examples/submesh_investigation/test_region_ds_nitsche.py). - .claude/skills/: 3 sites (adapt-on-top-faults x2, free-surface-convection x1). - CLAUDE.md: 1 signature reference (free-slip BC preference section). The ~1,370 legacy-trio (add_dirichlet_bc/add_natural_bc/add_essential_bc) sites already conform and are untouched per the D2 decision. The audit review documents under docs/reviews/2026-07/ record the pre-decision state as evidence and are not swept. Discovered while verifying the swept tests run warning-free: the Wave C zero-datum guard in add_rotated_freeslip_bc rejects FLOAT zero (sympy.sympify(0.0) != 0 is structurally True), so the canonical add_rotated_freeslip_bc(0.0, boundary) raises NotImplementedError while conds=0 works. Filed as issue #336 with a TODO(BUG) marker at the guard (comment-only src touch); the swept call sites use the working integer form add_rotated_freeslip_bc(0, boundary). No fix applied here (Charter section 9 scope discipline). Findings: API-01/API-02 sweep (WE-09, REMEDIATION-WORKLIST.md). Underworld development team with AI support from Claude Code

lmoresi added 5 commits June 18, 2026 16:55

Merge remote-tracking branch 'origin/development' into bugfix/jacobia…

d7e0cd4

…n-unwrap-to-constants # Conflicts: # src/underworld3/cython/petsc_generic_snes_solvers.pyx

Copilot AI review requested due to automatic review settings June 19, 2026 10:35

Copilot started reviewing on behalf of lmoresi June 19, 2026 10:36 View session

Copilot AI reviewed Jun 19, 2026

View reviewed changes

Merge remote-tracking branch 'origin/development' into bugfix/jacobia…

e657dd4

…n-unwrap-to-constants

lmoresi mentioned this pull request Jun 24, 2026

Nitsche free-slip BC: inconsistent Jacobian for anisotropic (transverse-isotropic) constitutive models #239

Open

Merge remote-tracking branch 'origin/development' into bugfix/jacobia…

6c0a54f

…n-unwrap-to-constants # Conflicts: # src/underworld3/cython/petsc_generic_snes_solvers.pyx

lmoresi merged commit c63cd70 into development Jul 2, 2026
1 check passed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Consistent Jacobian tangent (unwrap-before-differentiate fix) — opt-in, default-off#258

Consistent Jacobian tangent (unwrap-before-differentiate fix) — opt-in, default-off#258
lmoresi merged 7 commits into
developmentfrom
bugfix/jacobian-unwrap-to-constants

lmoresi commented Jun 19, 2026

Uh oh!

Copilot AI left a comment

Uh oh!

bknight1 commented Jun 25, 2026

Uh oh!

lmoresi commented Jul 2, 2026

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

		# Restore a clean Picard tangent for any subsequent solve (next step).
		self._set_newton_alpha(0.0)

		# (ramp the softmin softness δ→0), not a smooth tangent. See the design doc
		# docs/developer/design/jacobian-unwrap-constants-bug.md. The generic

Uh oh!

Conversation

lmoresi commented Jun 19, 2026

Consistent Jacobian tangent for nonlinear (viscoplastic) solves — opt-in, default-off

The bug

The fix (default-off, bit-identical by default)

Non-regression evidence

Scope / why a smooth Jacobian is not the VEP fix

Uh oh!

Copilot AI left a comment

Choose a reason for hiding this comment

Pull request overview

Reviewed changes

Uh oh!

bknight1 commented Jun 25, 2026

Uh oh!

lmoresi commented Jul 2, 2026

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants