Picard warm-up is a frozen-tangent step, not an nrichardson sweep (#791)

[lmoresi]

Fixes #791.

The defect

The Layer 1 design specifies the cold-start warm-up and solve(picard=N) as Picard (frozen-coefficient) iterations: linear Stokes solves, including the Schur complement, with the viscosity held at the current state. The standard SNES_Stokes.solve path ran SNES nrichardson with no nonlinear preconditioner instead: x <- x - lambda F(x), a residual step with no linear solve. It is not a Picard step and is nearly inert. The code comment contradicted itself on this.

The change

The standard path now follows the semantics the rotated free-slip path already had (utilities/rotated_bc.py):

consistent_jacobian	solve(picard=N)
False (default)	Every iteration already uses the frozen tangent, so this is a no-op.
"continuation"	The alpha=0 stage runs max(N, natural stage-1 count) iterations, both measured against the original residual. Per-stage counts are recorded in _continuation_stages, which is reset on every solve.
True	Raises NotImplementedError on a nonlinear residual, because no frozen tangent is compiled. A no-op on a linear residual.

The automatic cold-start warm-up is removed rather than repaired, since it was never a Picard step. When the rest state is exactly zero (homogeneous essential BCs, no stress history), Newton's first iteration is the Picard step. A boundary-driven or stress-history problem does not have that property, so for those the documented route to a Picard entry is "continuation" + picard=N. The design doc, both solver skills and test_0201 are updated to say so, and the dead phase="picard" branch of _warn_on_divergence is removed.

Behaviour note. Under the default frozen tangent, solve(picard=N) is now a no-op; previously it ran N nrichardson sweeps. Five example scripts use it this way (Ex_SurfaceIntegrals_BCs.py, Ex_Convection_Cylinder-FS.py, notched_beam_cascade.py, Ex_VP_Spiegelman_Benchmark.py). None of them raise.

Tests

• tests/test_1068_picard_warmup_is_a_frozen_tangent_step.py: 7 tests with hard baselines. Four of them fail on unpatched development.
  • With the frozen tangent, picard=3 matches picard=0 exactly.
  • With continuation, picard below the natural stage count gives exactly the natural count, and above it gives exactly N. This guards the case where PETSc's rtol restarts with each snes.solve call.
  • _continuation_stages does not go stale between solves.
  • With pure Newton, picard raises on a nonlinear residual and is a no-op on a linear one.
  • From an exactly-zero rest state, the first Newton step matches the first Picard step to 1e-6. The test documents that this does not extend to boundary-driven problems.
• The targeted solver suites (homotopy, cold start, warm start, rotated free-slip, Jacobian layout) pass 69/69.
• tier_a: 1,700 passed, 0 failed, with test_0053_hang_watchdog excluded. That test segfaults inside CPython's faulthandler, depending on accumulated full-suite state, with no UW code on the stack. It passes alone and when run with exactly its preceding test set. Filed separately as test_0053 hang_watchdog: SIGSEGV inside CPython faulthandler dump_traceback during full tier_a runs (state-dependent) #793. The first commit of this branch passed full tier_a (1,714/1,714).

Benchmark (solver-stability rule)

Spiegelman notch, ref 1, sqrt soft-min, cold start under Newton, cap 250 iterations. Before and after are built from the same base commit (1b4f3b5). KSP iterations are the cost measure.

case	before (nrichardson)	after
delta=0.4, tol 1e-6	converged, 80 nl, 1,287 KSP	converged, 97 nl, 1,057 KSP (-18%)
delta=0.1, tol 1e-6	linear-solve stall at 0.325, 10,046 KSP	identical stall, 10,054 KSP
delta=0.4, tol 1e-8	stall at 1.9e-3, 11,143 KSP	stall at 2.2e-3, 5,961 KSP (-47%)

Removing the sweep never costs Krylov effort.

Adversarial review

An independent review found one blocker: with "continuation", picard=N gave N iterations plus a full stage 1 restarted from scratch (measured 3+20 frozen iterations where the natural count was 10), because PETSc measures rtol from each call's own starting residual. It also found that the original test could not catch this, that the "from rest, Newton = Picard" claim was false for boundary-driven problems (first steps 45% apart on a sheared box), and that the nonlinear-solver skill and test_0201 still described the removed warm-up, plus smaller items. All are fixed in this branch.

🤖 Generated with Claude Code

https://claude.ai/code/session_01RdcatxUnjHx5nSP6DeGXGb

[Copilot]

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