Nonlinear solver defaults: Newton with one cold Picard step, smooth yield law, declared plastic regularisation (#791)

[lmoresi]

Fixes #791. Changes the nonlinear solver defaults and adds a declared plastic regularisation.

What changes for a user

	before	after
consistent_jacobian (every solver class)	False: Picard throughout	True: Newton, with ONE automatic Picard step on a cold start
ViscoPlasticFlowModel yield law	exact hard Min	smooth: softmin, sqrt family, anchor "onset", δ = 0.1
Krylov restart	PETSc's 30	100 (outer KSP of every class, Stokes velocity split)
solve(picard=N) under Newton	raised NotImplementedError	N real Picard steps, then Newton

The smooth default is exact on the viscous branch, so a model that never yields gives the viscous answer. It is exact in the plastic limit too. At the yield point it under-stresses by 4.5%, and its yielded tangent is at least −1e-4 of the viscous one. yield_mode="min" and consistent_jacobian=False remain as explicit choices.

The automatic Picard step is a real frozen-tangent step, not the former nrichardson sweep (#791). It is skipped:

• for a ksponly solve;
• for a flux whose coefficient does not depend on the unknowns;
• on the rotated free-slip and fault-contact paths, which run their own loops. An explicit picard=N on the rotated path stages Picard, then Newton.

The step runs inside the solve's update hooks, wall-clock guard, estimate_difficulty cap and report. The Newton stage that follows stops at tolerance·‖F(u_cold)‖.

New

• Parameters.plastic_rate_strengthening (η_reg) on ViscoPlastic, VEP and TI-VEP. The plastic viscosity becomes η_pl = τ_y/(2ε̇_II) + m·η_reg.

  • It floors the plastic viscosity AND the tangent along the flow. Without it, a yielded layer is a near-singular mechanism mode of the velocity block (notch: Rayleigh quotient 1.8e-7).
  • Its scale is the model's. There is no universal dimensionless value: a 1% overstress on the notch is 20 σ_y in a convective lid. The default is therefore 0 (no term).
  • A pressure-dependent yield stress under the consistent tangent without it warns.

• solve(homotopy="rate_strengthening") ladders m from the viscous limit, read from the fields, down to exactly 1 (ViscoPlastic only). Any other homotopy string now raises.

• newton_pressure_coupling=False: the Newton velocity block with the frozen pressure column. Far from the solution, the full tangent's d(η)/dp term makes the Schur complement indefinite (notch: J_vp 41% off −Bᵀ, 40/1968 eigenvalues change sign). The default stays True, which is benign in the basin of convergence.

• consistent_jacobian="continuation": one solve, with α keyed on the residual relative to ‖F(u=0)‖. It backs off when F climbs and re-enters warm, with Picard iterations, on divergence.

• solve_report.config records what actually ran (Solver introspection does not report the numerical configuration a solve ran with #806):

  • SNES, line search, KSP with its restart, and the fieldsplit sub-solvers;
  • the tangent (including warm-up steps and α);
  • for Stokes, the pressure gauge and the yield-law regularisation.

  The transcript's solve event carries the tangent, gauge and regularisation. A mean-pressure nullspace feeding a pressure-dependent viscosity warns once per solve.

Behaviour notes

• Iteration counts move for every nonlinear solve. Results move by the smooth law near yield (up to 4.5% at the yield point), never on the unyielded branch.
• solve_report.nl_its, ksp_its and history include the warm-up, and reduction is measured from the cold residual.
• picard is now the LAST argument of the Scalar/Vector/MultiComponent solve(). Five subclasses pass _force_setup positionally.
• Changing consistent_jacobian or newton_pressure_coupling after a solve recompiles the Jacobian on the next one.
• homotopy=True runs the δ march under the model's anchor ("onset" by default). Pass anchor="yield" for an approach from above.

Benchmark (solver-stability rule)

Spiegelman notch, refinement 1 (two multigrid levels: too few to test the multigrid, and silent on the resolution trend, which is the open measurement), cold start from zero, tolerance 1e-6, cap 250 iterations, same build and fixture. Only the changed defaults differ. KSP iterations are the cost measure.

η_bg	old defaults: exact Min, Picard, restart 30	new defaults: smooth (onset, δ 0.1), Newton + 1 Picard step, restart 100
1e24	DIVERGED_LINE_SEARCH, 43 its, 81 KSP, reduction 1.4e-2	converged, 144 its, 2,656 KSP
5e24	DIVERGED_LINE_SEARCH, 35 its, 125 KSP, reduction 9.0e-3	DIVERGED_MAX_IT, 251 its, 9,785 KSP, reduction 1.6e-3
5e24, cap 1000		DIVERGED_MAX_IT, 1,001 its, 50,544 KSP, reduction 1.54e-3

The hard case, on the new defaults, by the two routes the guides give:

5e24	result
consistent_jacobian="continuation", one solve	converged, 115 its, 768 KSP
η_reg = 1e-3 η_bg, homotopy="rate_strengthening"	11 rungs, m 1e6 → 1 (two rungs retried), converged at m = 1: 441 its, 2,912 KSP in total

The old defaults stop on a line-search failure at a reduction of about 1e-2 on both cases. The new defaults converge the 1e24 case in one solve. On 5e24 a single Newton solve reaches 1.6e-3 by iteration 140 and then stalls: monotone, but 1.604e-3 at 150 and 1.543e-3 at 1000 (0.99995 per iteration, so more iterations do not help), at about 39 KSP per iteration (13–15 in the basin of convergence), because far from the solution the full tangent's linear solve is the expensive part. That case needs the blend or the declared ladder, as plasticity-solvers and nonlinear-solver say.

Tests

• test_1068: the warm-up contract per tangent.
  • A warm-up step equals one frozen-tangent iteration on a boundary-driven box.
  • N steps converge to the same answer, and a warm solve takes none.
  • The Newton stage is anchored to the cold residual.
  • _force_setup is not read as picard.
  • A tangent change recompiles, and α is per solver.
  • A flux made nonlinear after a solve takes the cold step.
  • The partial tangent restores J_vp = −Bᵀ.
• test_1069: the gauge record and warning.
  • The Krylov restart is checked on the live KSP, after a solve.
  • Switching to Drucker-Prager after a solve is seen by the record.
• test_1071: η_reg floors the plastic viscosity and is recorded. Unstated, 0 and 0.0 are the same model and warn. The ladder runs and ends at m = 1. A misspelt homotopy raises. VEP carries the term and refuses the ladder.
• test_1055: the default law pinned to its closed form to 1e-12.
• tests/parallel/test_1081: the cold-start verdict is rank-uniform. The old per-rank norm disagreed at np=2 and np=3, and in a solve it hung the job.
• test_1018: rotated free-slip keeps frozen-tangent coverage and stages an explicit picard under Newton.
• Updated for the new defaults: 0104, 0201, 0641, 1057, 1059, 1066.
• Full serial suite (tests/test_*.py without tier C, plus the serially-run parallel guard; test_0053_hang_watchdog excluded, test_0053 hang_watchdog: SIGSEGV inside CPython faulthandler dump_traceback during full tier_a runs (state-dependent) #793) on 6535935: 2,998 passed, 0 failed. tests/parallel/test_1081 passes at np=2 and np=3. The style gate and the test-coverage check pass.

Adversarial review

Two rounds: three reviewers on the branch, then one on the fix commit. The findings, and the attacks that failed, are in the review comment below. All blockers are fixed in this branch.

Open, not fixed here:

• Automatic cold Picard warm-up diverges a LINEAR constrained Stokes solve (Stokes_Constrained + multiplier free-slip) #817 is misdiagnosed. The constrained Stokes solver fails on its second Newton iteration whatever the tangent; details are on the issue.
• Whether VEP should default to the frozen tangent: the guide still recommends it, and test_1052 converges under the new default.

Merged with development (#823 tier 2, #835), 2026-10-11

Development's JIT now lowers kernels onto a graph of named quantities by default, with the previous route kept as jit_route = "expanded". The merge (2c6bc2f) takes the solver's route in _jacobian_source and keeps this branch's warm-up blend; details are in the review record, round 3.

The merge review found one blocker this branch's default exposes, fixed here:

• A fractional power of the bare strain-rate invariant (1/(1 + Einv2**(1/3))) had a NaN Newton tangent at rest on both routes, so the cold start raised PETSc error 77. The guard in _jacobian_unwrap now covers every non-integer exponent (Closes Newton tangent is NaN at rest for a fractional (non-half) power of the strain-rate invariant: the Jacobian guard covers half-integer powers only #841). test_1067 gains the regression; it failed on both routes before the change.
• The rate-strengthening probe projection now compiles on its solver's route, and three docstrings that the merge made false are corrected.

Tests on 72f0bd8: the serial suite (test.sh batches, not tier C) passes on both JIT routes, 3,139 each, 0 failed; the parallel directory at np=4 passes apart from test_1063[ti], which is #826 and fails the same way on development.

Open from the merge review (for decision, not fixed here): solve(homotopy=True) ends the delta march at the power-mean law with delta = 9.8e-4 rather than the model's stated delta = 0.1; on the DP box the automatic warm-up takes 19-20 iterations against 9 with picard=-1.

Underworld development team with AI support from Claude Code

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.

[lmoresi]

Adversarial review — nonlinear solver defaults and the declared plastic regularisation (PR #794)

Branch bugfix/picard-warmup-791 vs development.

Method: two rounds.

• Round 1 (at a8aa68c): three independent reviewers, one each on the solver control flow, the constitutive models and η_reg, and the solve record, docs and tests.
• Round 2: one reviewer on the fix commit 8010d33.

Every reviewer was told to break the change and to verify each finding by probe. Probe scripts are in ~/+Simulations/review_794/. All findings below are fixed on the branch (8010d33, 6535935) unless listed as open.

Round 1 — blockers

• Cold-start verdict was per rank (_solution_is_trivially_zero read u.vec.norm(), a LOCAL vector).

  • At np=2, with the initial guess non-zero in one corner, rank 0 took the warm-up and rank 1 did not. Rank 1 failed with MPI_ERR_BUFFER; rank 0 spun until killed.
  • The verdict is now a reduced any-nonzero test. tests/parallel/test_1081 checks it; the old verdict disagrees at np=2 and np=3.

• Krylov restart 100 never took effect for Stokes. It was written inside the strategy setter, which construction never runs; the live KSP and the report both showed 30. It is now set once at construction as a managed option.

• picard was inserted second in the Scalar/Vector/MultiComponent solve(). Five subclasses pass _force_setup positionally, so _force_setup=True became picard=1: a warm-up and a blended compile on a linear AdvDiffusion solve. picard is now the last argument.

• solve(picard=N) raised on rotated free-slip under the new default. It now stages N Picard iterations, then Newton (power law: 17 its, error 1.2e-8 against the essential-BC answer).

• The default yield anchor "yield", sqrt family, δ 0.1, had two defects:

  • every unyielded cell was 1/(1 − δ/2) stiffer: η/η_ve = 1.0499 at f → 0, and a box that never yields moved 4.75%;
  • the yielded tangent was negative: τ/τ_y peaks at 1.033 at f = 1.2, and dτ/dε̇ = −0.029 at f = 1.5.

  The maintainer chose sqrt / "onset" / δ 0.1: η(0) = 1, τ(f=1) = 0.955 τ_y, tangent ≥ −1e-4.

• A YAML colon in the guide front matter dropped nonlinear-solver from uw.capabilities(); test_0030 failed 2 tests.

• An uncommented except: pass in yield_continuation.py failed the Charter S4 style gate.

Round 1 — should-fix

• Linearity test on every warm solve. The symbolic test ran before the cheap cold-start test: 3.7 s per warm step on a 16×8 VEP box, and 5.8 against 2.2 s/step overall. The cold test now runs first, and the facts are cached.
• The warm-up ran outside the solve's contracts: before the update hooks and the wall-clock guard, outside the estimate_difficulty cap (cap 1 ran 2 iterations), and missing from the report (fnorm0 was 4.04e-2 against a cold F of 0.120). It now runs inside _snes_solve_with_retries.
• Newton after the warm-up measured rtol from the post-warm-up residual, and stagnated on SNORM on near-linear boxes (2 its without the warm-up). It is now anchored at tolerance·‖F(u_cold)‖.
• Changing the tangent after a solve did not recompile, so the report described a tangent that was not running (Poisson False → True: 9 its, against 3 for a fresh solver). α was one process-wide container shared by every solver. Both setters now mark a rewire, and α is unique per solver.
• Ladder defects:
  • a retry restarted at the largest converged m (1e6), not the last good rung;
  • a zero or Piecewise η_reg gave zoo or stale zeros, and a ladder that did nothing reported converged=True;
  • the entry multiplier was absolute;
  • each rung recompiled (5.2 s against 0.23 s);
  • VEP failed with an AttributeError instead of a clean refusal;
  • the probe variable was named from id(solver), which differs between ranks.
• Zero test spelled three ways. Float(0.0) == 0 is False in SymPy 1.14, so plastic_rate_strengthening = 0.0 compiled a zero term and suppressed the warning.
• Warning hygiene. The gauge warning fired twice per solve from an internal frame. The Drucker-Prager warning fired under the frozen tangent and for a p-dependent η0 with a constant yield stress.
• The record had gaps. Rotated reports carried no config. The transcript dropped the regularisation. Guides and docstrings still described the old defaults, and the new controls were undocumented.
• Tests that could not fail. test_1069 asserted only that a "restart" key existed (it was 30). test_1071's string checks could not tell 0.05 from 0, and it never checked that the ladder ran. test_1018's rotated frozen-tangent coverage had silently become Newton.
• Duplicated code. The warm-up policy was pasted 4 times and had drifted: "continuation" ignored picard=N in three classes.

Round 2 (on 8010d33)

• blocker: the form-fact cache went stale across a rebuild. Viscous → Drucker-Prager, solved warm, recorded pressure_in_rheology=False with no warning. A flux made nonlinear skipped the cold step. The cache is now dropped wherever a rewire or a setup completes.
• blocker: test_1018 still expected the removed NotImplementedError. It now asserts the staging.
• should-fix: the rotated path discarded the warning notes.
• should-fix: KSPSetType on the single-field FMG route reset the live restart to 30 while the report said 100. The restart is re-applied, and test_1069 now reads it from the live KSP.
• should-fix: the ladder's early exits skipped its warning. A failed first rung marked non-solution fields as converged.
• nit: the duplicated homotopy check, the warm-up warning frame, the δ-march anchor docstring.

Attacks that failed

• Linearity guard: correct on 19 cases (Poisson k=1/k(T)/k(u), Stokes η=1/η(T)/η(|v|), VE, VP, DP, VEP, power law, units-aware viscous and VP, projections, AdvDiffusion, Navier-Stokes).
• Warm loop: no warm-up on steps 2+. Nonlinear Poisson matches Picard-throughout to 7.7e-11 (3 vs 9 its).
• Eisenstat-Walker: the warm-up iterate is identical with EW on and off.
• newton_pressure_coupling=False with a p-independent rheology: identical iterations, dv = 0, ‖J_T − J_F‖ = 0.
• Units: uw.quantity(1e20, "Pa*s") with η_ref 1e21 gives deep-yield η = 0.1005. The default 0 [Pa s] omits the term.
• Ladder: ends at exactly 1.0 on a clean run, after a mid-ladder failure, and after a failure at m = 1. It refuses NaN through a separate flag reduction, and is bounded by max_steps.
• Names: two VP models get distinct m_reg atoms. describe() names η_reg and m_reg.
• Record: JSON-safe. A transcript round-trips through uw.Transcript. The snapshot creates no NPC or line search and prints nothing.
• Probe cap: picard 2/5/3 with cap 3 gives 3 iterations. Tolerances are restored.
• Anchor: never loosens a solve (boundary-driven DP, VP, n = 0.5 and 3).
• Rank uniformity: the ladder runs at np=3. Every new collective is reached by every rank.
• Restart: reaches the live KSP for default, fast, robust, gamg and refinement=2 FMG Stokes, and for plain, gamg and fmg Poisson. A user value is kept. The composed solver keeps its 200.

Open (not fixed on this branch)

• Automatic cold Picard warm-up diverges a LINEAR constrained Stokes solve (Stokes_Constrained + multiplier free-slip) #817 is misdiagnosed. Stokes_Constrained fails on its second Newton iteration whatever the tangent:
  • with η(1 + 1e-6 ε̇_II): automatic warm-up DIVERGED_DTOL; picard=-1 DIVERGED_MAX_IT oscillating; consistent_jacobian=False identical;
  • the linearity guard only hides this for linear fluxes.
• VEP under the Newton default. The guide recommends the frozen tangent for VEP, but test_1052's loading through yield converges under the new default (nl = 1 per step).
• Retries in the record: with divergence retries, the report keeps the warm-up plus the last attempt (pre-existing behaviour for retries).
• mpi4py reductions: _solution_is_trivially_zero and _masked_max_ratio reduce through uw.mpi.comm (mpi4py), which the library already does elsewhere. A PETSc global-Vec norm would serve.

[github-actions]

✅ Test Suite: success

[lmoresi]

Adversarial review, round 3: the merge with development (2c6bc2f)

One reviewer attacked the seams between this branch's tangent machinery and development's two JIT routes, with probes on the installed build. The full record is in docs/reviews/2026-10/nonlinear-defaults-794-adversarial-review.md, round 3.

Fixed in this branch

• blocker (Newton tangent is NaN at rest for a fractional (non-half) power of the strain-rate invariant: the Jacobian guard covers half-integer powers only #841, 4cdc00d): a fractional power of the bare invariant, eta = 1/(1 + Einv2**(1/3)), raised PETSc error 77 on a cold solve under the default tangent, on both routes. SymPy merges sqrt(g)**(1/3) into g**(1/6), and the guard covered half-integer powers only, so J(0) was NaN and the warm-up blend inherited it at alpha = 0 (0 * NaN). The guard now covers every non-integer exponent on a base with free symbols, on both routes. In floating point the added 1e-36 vanishes against any base above about 1e-20, so only states at rest see it.
• should-fix (bf4c32f): the rate-strengthening probe projection compiled on the process default route; it now takes its solver's route.
• should-fix (bf4c32f, 3eaa99b): docstrings and the design note that the merge made false (the blend under True, alpha = 0 exact only while the Newton blocks are finite, "default off").

Open, for decision

• solve(homotopy=True) ends the delta march at the power-mean law with delta = 9.8e-4, not at the model's stated delta = 0.1 and sqrt smoother; the manifest carries {s_{y}} in place of {\updelta_{y,0}}.
• On the DP box the automatic warm-up doubles the iteration count (19-20 against 9 with picard=-1): after the Picard step the line search sits at lambda = 0.1 for 13 steps. plasticity-solvers.md reads a linear tail as "you are on the Picard tangent", which this case contradicts.
• Every viscoelastic Stokes solve prints a PETSc error banner from _check_velocity_preconditioner #827: the PETSc banner also prints after a rewire followed by a warm solve that converges at iteration 0 (triggers added to the issue).
• test_1063 [ti]: constrained free-slip topography 1.8% off its recorded anchor on the anchor's own mesh #826: test_1063[ti] at np=4, unchanged from development (37.7525479 against the 37.078 anchor).

Attacks that failed

• Every _jacobian_unwrap and getext call passes the solver's route. With the solver on "expanded", the Jacobian blocks hold no graph nodes and the C no uwt_ temporaries, for the default tangent, "continuation" and newton_pressure_coupling=False.
• \alpha_{N}, the yield softness, m_reg and \eta_reg are constants[] slots on both routes; ramping them leaves the cache key unchanged.
• A mid-run route switch keeps _picard_blend and recompiles the blend on the new route with the same answer. Form facts are structural (_reaches), so they cannot go stale with the route.
• DP box, default law and tangent: 20 iterations on both routes, solutions equal to 2.2e-16 (v) and 7e-16 (p), J(0) finite and equal.
• test_0104 under the default tangent: the ramp keeps its cache key and matches a fresh build to 1e-13, and the stale-slot check still raises on a warm solve, so the frozen-tangent pin hides nothing.

Tests on 72f0bd8: serial suite on both routes 3,139 passed, 0 failed (before the guard fix, the same suite failed only the two new test_1067 cases); tests/parallel at np=4, all pass except #826.

Underworld development team with AI support from Claude Code

[github-actions]

✅ Test Suite: success