underworldcode · lmoresi · Sep 26, 2026
@@ -1,6 +1,6 @@
 ---
 name: nonlinear-solver
-description: How to make a hard nonlinear Stokes solve (Drucker-Prager / yield-stress viscoplastic) CONVERGE reliably in Underworld3 the way the working recipe actually does it — automatic warm-start (one Picard step on a cold start) plus a MULTI-SOLVE δ-continuation (constant δ per solve, warm-start the next, sharper δ), the consistent-Newton tangent, and a non-symmetry-safe multigrid smoother. Reach for THIS when a viscoplastic solve stalls / diverges and you are about to hand-tune PETSc options, ramp δ, or "just add a monitor". It carries the CONFIG TRAP LIST — the setup mistakes that each produce a different failure a few steps in — and the one thing you must NOT do (ramp δ inside a single SNES solve). For the yield-law maths and which tangent per model, see `plasticity-solvers`.
+description: How to make a hard nonlinear Stokes solve (Drucker-Prager / yield-stress viscoplastic) CONVERGE reliably in Underworld3 the way the working recipe actually does it — a Picard entry where one is genuinely needed (`consistent_jacobian="continuation"` + `picard=N` — there is NO automatic one, #791) plus a MULTI-SOLVE δ-continuation (constant δ per solve, warm-start the next, sharper δ), the consistent-Newton tangent, and a non-symmetry-safe multigrid smoother. Reach for THIS when a viscoplastic solve stalls / diverges and you are about to hand-tune PETSc options, ramp δ, or "just add a monitor". It carries the CONFIG TRAP LIST — the setup mistakes that each produce a different failure a few steps in — and the one thing you must NOT do (ramp δ inside a single SNES solve). For the yield-law maths and which tangent per model, see `plasticity-solvers`.
 ---
 
 # nonlinear-solver
@@ -20,9 +20,14 @@ Yield-law maths, tangent-per-model, quadratic-convergence check: `plasticity-sol
 ## The recipe (what actually converges)
 
 1. **Warm start.** Start the continuation at **large δ**, where the yield surface is
-   smooth and the problem is easy, and take **one Picard step** into the Newton
-   basin. One Picard step is defect-correction iteration 1 — contractive, cheap. From
-   a *warm* iterate, take **no** Picard step (it wastes the good quadratic start).
+   smooth and the problem is easy. If a Picard entry is needed, it must be asked for:
+   `consistent_jacobian="continuation"` with `solve(picard=N)` gives N frozen-tangent
+   (defect-correction) iterations before Newton. From a *warm* iterate, take **no**
+   Picard step (it wastes the good quadratic start). ⚠️ #791: the former "automatic
+   Picard step" on a cold start was an nrichardson residual sweep, not a Picard step,
+   and has been removed. With homogeneous BCs and no stress history, Newton's first
+   step from rest coincides with the Picard step anyway; a boundary-driven or
+   stress-history problem gets no such entry (measured 45% apart on a sheared box).
    A cold `v=0` start is safe on its own terms: `ε̇=0` makes `η_pl` infinite, which
    the soft-min carries to the viscous branch (see the trap list for the one form
    that must be written carefully).
@@ -45,7 +50,8 @@ Yield-law maths, tangent-per-model, quadratic-convergence check: `plasticity-sol
    and the surviving evidence), and the driver's documented cold-start
    guarantee does not currently hold (issue #473: entry can fail on a
    pressure-dependent yield, and the step control is effectively one-shot).
-   Newton + the automatic Picard entry handles the standard cases without it.
+   Newton (with `"continuation"` + `picard=N` where a Picard entry is needed) handles
+   the standard cases without it.
 
 3. **Consistent-Newton tangent** for non-elastic DP (`consistent_jacobian=True`);
    **Picard** for elastic VEP — see `plasticity-solvers` for the per-model table.
@@ -98,23 +104,26 @@ mesh-mover — the `is_setup=False` hook); kept through coefficient changes (vis
 δ, BC values, time step). A **diverged** solve leaves it `False`, so the next solve
 auto-cold-starts rather than warming off a corrupted iterate.
 
-On a **cold** (`zero_init_guess=True`) Stokes solve under the **consistent-Newton
-tangent**, a single Picard step is now taken automatically (reusing the existing
-`picard=1` machinery). The default (frozen) tangent path is bit-identical.
+There is **no automatic Picard step** on a cold start (#791 — the former one was an
+nrichardson residual sweep and has been removed). For a genuine Picard entry:
 
 ```python
-stokes.consistent_jacobian = True
-stokes.solve()                 # cold → one automatic Picard step, then Newton
+stokes.consistent_jacobian = "continuation"
+stokes.solve(picard=3)         # >= 3 frozen-tangent iterations, then Newton
 if stokes.has_solution:
     ...
 ```
 
+Under `consistent_jacobian=True`, `solve(picard=N)` raises on a nonlinear residual (no
+frozen tangent is compiled); under `False` every iteration is already a Picard step.
+
 ---
 
 ## Implementation status (this line of work)
 
-- **Layer 1a — DONE:** `has_solution` + cold consistent-Newton Picard warm-up
-  (`petsc_generic_snes_solvers.pyx`; test `test_0201`).
+- **Layer 1a — DONE (revised by #791):** `has_solution`. The cold-start "Picard
+  warm-up" was an nrichardson sweep and is removed; Picard entry is explicit via
+  `"continuation"` + `picard=N` (test `test_1068`).
 - **Layer 1b — DONE:** `zero_init_guess` is tri-state — `None` (default) auto-detects
   from `has_solution`, `True` forces fresh, `False` insists on warm. Note warm and cold
   agree only to the *convergence tolerance*, not bitwise.
@@ -127,7 +136,7 @@ if stokes.has_solution:
   continuation, returning the march summary. The doctrine that made this the
   recommended entry point was retracted (unit-scaling error — see
   `plasticity-solvers`), and its cold-start guarantee is broken (issue #473);
-  use it after Newton + Picard entry and grid sequencing have failed.
+  use it after Newton (with an explicit Picard entry) and grid sequencing have failed.
 
 ---
 

@@ -1,6 +1,6 @@
 ---
 name: plasticity-solvers
-description: How to get hard-Min viscoplastic / visco-elastic-plastic (VEP) Stokes solves to CONVERGE in Underworld3 — Newton with the automatic Picard entry (solver.consistent_jacobian), which tangent per model, grid sequencing for the hard cases, and the δ-soft-min substrate (yield_mode / yield_smoother / yield_anchor) as a modelling choice. Reach for THIS first when a Drucker-Prager / yield-stress Stokes solve stalls, diverges (DIVERGED_LINEAR_SOLVE / line-search fail), or grinds through ~20+ nonlinear iterations. Tells you which tangent to use per model, how to confirm you are actually running Newton, and the measured failure modes. For the solver-config trap list and multigrid, see `nonlinear-solver`.
+description: How to get hard-Min viscoplastic / visco-elastic-plastic (VEP) Stokes solves to CONVERGE in Underworld3 — Newton (solver.consistent_jacobian) with an explicit Picard entry via "continuation" + picard=N, which tangent per model, grid sequencing for the hard cases, and the δ-soft-min substrate (yield_mode / yield_smoother / yield_anchor) as a modelling choice. Reach for THIS first when a Drucker-Prager / yield-stress Stokes solve stalls, diverges (DIVERGED_LINEAR_SOLVE / line-search fail), or grinds through ~20+ nonlinear iterations. Tells you which tangent to use per model, how to confirm you are actually running Newton, and the measured failure modes. For the solver-config trap list and multigrid, see `nonlinear-solver`.
 ---
 
 # plasticity-solvers
@@ -16,7 +16,7 @@ and records what was retired.
 stokes.constitutive_model = cm           # ViscoPlastic / ViscoElasticPlastic / TI-VEP
 cm.Parameters.yield_stress = tau_y       # finite -> plasticity active
 stokes.consistent_jacobian = True        # Newton tangent (non-elastic DP; see table)
-stokes.solve()                           # cold start takes ONE Picard step automatically
+stokes.solve()                           # no automatic Picard entry (#791) — see item 1
 ```
 
 ---
@@ -26,20 +26,27 @@ stokes.solve()                           # cold start takes ONE Picard step auto
 Yielding viscoplasticity is `η_eff = Min(η_visc, η_yield)`,
 `η_yield = τ_y/(2·ε̇_II)`. The `Min` kink is what makes it hard.
 
-1. **Picard is an ENTRY requirement, not an accelerator.** On a cold start under
-   the consistent tangent, `solve()` takes one automatic Picard (frozen-tangent)
-   step and then runs Newton (fires only when `picard==0`,
-   `consistent_jacobian is True`, and the start is cold — from a warm iterate,
-   0 Picard really is 0). Do NOT front-load Picard where Newton works: measured,
-   opening with 5 / 25 Picard steps cost 12 / 30 total iterations against pure
-   Newton's 7.
+1. **There is no automatic Picard entry; ask for one where it matters.** When the
+   rest state is exactly zero (homogeneous essential BCs, no stress history) Newton's
+   first step from rest IS the Picard step (pinned by `test_1068`, < 1e-6). A
+   BOUNDARY-DRIVEN or STRESS-HISTORY problem does not have that property — the driven
+   layer yields immediately (measured: first steps 45% apart on a sheared box) — so use
+   `consistent_jacobian="continuation"` + `solve(picard=N)` for a real entry. ⚠️ #791:
+   the old "automatic Picard entry" was never a Picard step — it ran SNES
+   `nrichardson` with no linear solve, a nearly inert residual sweep — and it has
+   been removed. The earlier measurement "opening with 5 / 25 Picard steps cost
+   12 / 30 iterations against pure Newton's 7" was of those nrichardson sweeps, NOT
+   of Picard, and is withdrawn.
 
 2. **Newton-first; spend Picard only to rescue.** When Newton fails
    (`DIVERGED_LINE_SEARCH` / `DIVERGED_LINEAR_SOLVE`) **or stalls admissibly**
    (steps accepted, residual flat — no FAIL reason ever fires), revert to the best
-   iterate and buy a Picard block: `solve(picard=N)`, or
-   `consistent_jacobian="continuation"` (staged Picard→Newton α-blend; α is a
-   `constants[]` atom, no recompile). Rescue on *stagnation*, not only on a
+   iterate and buy a Picard block. A Picard block needs the FROZEN tangent, which
+   the pure-Newton compile does not carry: use `consistent_jacobian="continuation"`
+   with `solve(picard=N)` (at least N frozen-tangent iterations, then Newton; α is a
+   `constants[]` atom, no recompile), or `consistent_jacobian=False` for Picard
+   throughout. Under `consistent_jacobian=True`, `solve(picard=N)` RAISES on a
+   nonlinear residual rather than silently running something else (#791). Rescue on *stagnation*, not only on a
    failure reason — the failure-only trigger measured byte-identical to doing
    nothing at the cliff.
 
@@ -53,7 +60,7 @@ Yielding viscoplasticity is `η_eff = Min(η_visc, η_yield)`,
 
 4. **`solver.has_solution` / tri-state `zero_init_guess`** make warm-started
    campaigns safe: `None` (default) auto-detects; a diverged solve or a remesh
-   clears the flag so the next solve cold-starts (with its Picard entry) instead
+   clears the flag so the next solve cold-starts instead
    of warming off a corrupted iterate.
 
 ---
@@ -86,7 +93,7 @@ rescue, and grid sequencing are.
 
 | Model | Use | Why |
 |-------|-----|-----|
-| `ViscoPlasticFlowModel` (non-elastic) | **`True`** (Newton) | Quadratic near the solution; the automatic Picard entry handles the cold start. |
+| `ViscoPlasticFlowModel` (non-elastic) | **`True`** (Newton) | Quadratic near the solution; for a boundary-driven cold start use `"continuation"` + `picard=N` for a Picard entry. |
 | `ViscoElasticPlasticFlowModel` (VEP) | **`False`** (Picard) | The consistent yield tangent over the elastic stress-history block makes the Jacobian **indefinite → `DIVERGED_LINEAR_SOLVE`**. Picard is contractive. |
 | `TransverseIsotropicVEPFlowModel` (TI-VEP) | **`False`** (Picard) | Same as VEP (elastic). |
 | Any, far from the solution | **`"continuation"`** | Staged Picard→Newton; Picard locates the basin, Newton finishes. Beat pure Newton at every notch point measured — but its stage switch is a residual LEVEL and one-way, so it can overspend Picard on easy problems. |
@@ -182,7 +189,7 @@ measured, it never has.
 |---------|-------|-----|
 | `DIVERGED_LINEAR_SOLVE`, 0 iters, VEP | consistent Newton over the elastic block → indefinite | Picard (`consistent_jacobian=False`) |
 | `DIVERGED_LINEAR_SOLVE` at nl=0 with a viscosity floor set | δ→0 leaves the floor's `Max` corner exact | set `viscosity_min_rounding` |
-| Newton stalls with no divergence reason | admissible uselessness — steps accepted, residual flat | revert to best iterate, Picard block (`picard=N` / `"continuation"`); consider grid sequencing |
+| Newton stalls with no divergence reason | admissible uselessness — steps accepted, residual flat | revert to best iterate, Picard block (`"continuation"` + `picard=N`; under pure Newton `picard` raises); consider grid sequencing |
 | Converges but σ sits **below** τ_y | a fixed δ>0 soft-min under the default `"onset"` anchor is a WEAKER law | that is the modelling choice you made — use `yield_anchor="yield"`, or δ→0 / `yield_mode="min"` for the exact surface |
 | Linear (~20-iter) convergence | Picard tangent when you wanted Newton | `consistent_jacobian=True` on a non-elastic model (see "Confirm" above) |
 

@@ -66,6 +66,29 @@ every nonlinear solver; layers 2–3 build on it.
 
 ### Layer 1 — automatic warm-start (all nonlinear solvers)
 
+```{warning}
+**Implementation status (#791, 2026-09-25).** As originally landed, this layer did not do what
+it specifies. The warm-up and `solve(picard=N)` ran SNES `nrichardson` with no nonlinear
+preconditioner — `x <- x - lambda F(x)`, a residual step with no linear solve — which is not a
+Picard step and is nearly inert. It is now corrected, with the semantics the rotated free-slip
+path already had:
+
+- `consistent_jacobian=False`: every iteration uses the frozen tangent, so `picard` is
+  satisfied by the solve itself.
+- `"continuation"`: `picard=N` guarantees at least N frozen-tangent iterations (stage 1, α=0).
+- `True`: `picard>0` raises on a nonlinear residual (no frozen tangent is compiled).
+
+The **automatic** cold-start warm-up is removed rather than repaired — it was never a Picard step,
+so nothing is lost. When the rest state is exactly zero (homogeneous essential BCs, no stress
+history) Newton's first iteration from rest IS the Picard step (verified to < 1e-6 in
+`tests/test_1068_picard_warmup_is_a_frozen_tangent_step.py`). A boundary-driven problem does NOT
+have this property — the Dirichlet values yield the driven layer immediately (first steps measured
+45% apart on a sheared box) — nor does a stress-history model; for those, a genuine Picard entry is
+`consistent_jacobian="continuation"` with `solve(picard=N)`. The Rules and cold/warm text below
+describe the ORIGINAL design intent of an automatic entry, which is not what the code now does.
+Measurements quoted elsewhere for "opening Picard steps" were of the nrichardson sweep.
+```
+
 A single **Picard (frozen-coefficient) step is a general cold-start warm-up**. It is
 defect-correction iteration 1: contractive, moves a cold guess into the Newton
 basin. UW3 already exposes it (`solve(picard=N)`) and already has the Picard tangent