The automatic cold-start 'Picard' warm-up (and solve(picard=N)) runs nrichardson with no linear solve — the design specifies a frozen-coefficient Picard step

[lmoresi]

What

The Layer 1 design (docs/developer/design/nonlinear-solver-homotopy-warmstart.md) specifies the automatic cold-start warm-up as:

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.

A Picard step is a linear Stokes solve (velocity block, Schur complement) with the viscosity frozen at the current state. From a cold start the strain rate is ~0, so the yield branch is inactive and that frozen operator is the viscous one: the Picard step is the viscous solve. That is the whole point of it.

What the code does (SNES_Stokes.solve, petsc_generic_snes_solvers.pyx):

if picard != 0:
    self._push_snes_max_it(abs(picard))
    self.snes.setType("nrichardson")
    self.snes.setFromOptions()
    self.snes.solve(None, gvec)

SNESNRICHARDSON with no nonlinear preconditioner (nothing in the solver attaches an NPC) is x <- x - lambda F(x) — a residual step with no operator inverse and no linear solve. It cannot produce the viscous solution; each unknown moves only by its own residual.

The same path serves solve(picard=N), so the user-facing Picard option does not perform Picard iterations either.

The code comment on that block contradicts itself, which is how this survived review:

• line ~11208: "Automatic cold-start warm-up (Layer 1): a single Picard (frozen-coefficient) step moves a cold guess into the Newton basin — it is defect-correction iteration 1"
• line ~11240: "What this sweep IS: one nrichardson step, x <- x - lambda F(x), with no linear solve and no frozen tangent. It is not a Picard step and it is nearly inert (1-12% residual reduction on a linear Stokes, measured)"

Why it matters

• The documented cold-start behaviour of every nonlinear Stokes solve is not what runs. Cold solves currently get an inert residual step followed by Newton from an essentially-zero state, instead of Newton from the viscous solution.
• solve(picard=N), the documented rescue for stalled viscoplastic solves (see the plasticity-solvers skill), is not a Picard rescue.
• It confounds measurements. On the Spiegelman notch (ref 1, sqrt soft-min δ=0.4, tol 1e-10), a cold solve from zero and a warm solve from an explicit viscous seed ended 25× apart in residual after 80 Newton iterations (2.3e-3 vs 5.8e-2), and it took an extra round of bisection to establish that the cold path's "Picard step" was not the viscous solve it is documented to be.

Suggested fix

Make the warm-up a real Picard step: one Newton iteration with the frozen (Picard) tangent. The machinery exists — consistent_jacobian="continuation" at α = 0 is documented as bit-identical to the Picard tangent and is constants[]-routed, so it needs no recompile. Alternatively, attach an NPC to the nrichardson SNES so it becomes preconditioned defect correction.

Per the solver-stability rule this wants benchmarking before it lands: the cold-start path runs in every nonlinear solve, and the linear-problem short cut (ksponly) must stay untouched.

Related: #778 (solve_history loses the Picard stage of a continuation solve).