From 520b7d14aefbeea37c791441e6a18885ff8e15e5 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 25 Sep 2026 18:26:12 -0700 Subject: [PATCH] Picard warm-up is a frozen-tangent step, not an nrichardson sweep (#791) The Layer 1 design specifies the cold-start warm-up and solve(picard=N) as Picard (frozen-coefficient) iterations: linear Stokes solves with the viscosity held at the current state. The standard path instead ran SNES nrichardson with no nonlinear preconditioner (x <- x - lambda F(x)), a residual step with no linear solve, and not a Picard step. The standard path now has 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 makes the alpha=0 stage run max(N, natural stage-1 count) iterations, both measured against the original residual (stage 1 after the N-iteration block uses an absolute target, since PETSc's rtol restarts with each snes.solve call). Per-stage counts are recorded in _continuation_stages, reset every solve. - True: picard>0 raises on a nonlinear residual (no frozen tangent is compiled) and is a no-op on a linear one. The automatic cold-start warm-up is removed rather than repaired; 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 (test_1068, < 1e-6). A boundary-driven or stress-history problem does not have that property (first steps measured 45% apart on a sheared box); a genuine Picard entry there is "continuation" + picard=N. Docs, both solver skills and test_0201 are updated to say so; the dead phase="picard" branch of _warn_on_divergence is removed. Behaviour note: solve(picard=N) under the default frozen tangent is now a no-op (it previously ran N nrichardson sweeps). Five example scripts use it that way; none raise. Tests: test_1068 (7 tests, hard baselines). Targeted solver suites 69/69; tier_a 1714 passed, 0 failed (first revision). Notch benchmark, cold Newton, before/after on the same base commit: KSP effort -18% (delta=0.4, tol 1e-6, both converge) and -47% (tol 1e-8, both stall at the same residual); delta=0.1 stalls identically on both. Underworld development team with AI support from Claude Code --- .claude/skills/nonlinear-solver/SKILL.md | 35 ++-- .claude/skills/plasticity-solvers/SKILL.md | 37 ++-- .../nonlinear-solver-homotopy-warmstart.md | 23 +++ .../cython/petsc_generic_snes_solvers.pyx | 174 ++++++++++++------ ...test_0201_solver_has_solution_warmstart.py | 15 +- ..._picard_warmup_is_a_frozen_tangent_step.py | 172 +++++++++++++++++ 6 files changed, 368 insertions(+), 88 deletions(-) create mode 100644 tests/test_1068_picard_warmup_is_a_frozen_tangent_step.py diff --git a/.claude/skills/nonlinear-solver/SKILL.md b/.claude/skills/nonlinear-solver/SKILL.md index 297b3df85..61f0472f4 100644 --- a/.claude/skills/nonlinear-solver/SKILL.md +++ b/.claude/skills/nonlinear-solver/SKILL.md @@ -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. --- diff --git a/.claude/skills/plasticity-solvers/SKILL.md b/.claude/skills/plasticity-solvers/SKILL.md index d5f83bfd0..3c13b627a 100644 --- a/.claude/skills/plasticity-solvers/SKILL.md +++ b/.claude/skills/plasticity-solvers/SKILL.md @@ -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) | diff --git a/docs/developer/design/nonlinear-solver-homotopy-warmstart.md b/docs/developer/design/nonlinear-solver-homotopy-warmstart.md index 27d58e0ee..0c7b6f9ae 100644 --- a/docs/developer/design/nonlinear-solver-homotopy-warmstart.md +++ b/docs/developer/design/nonlinear-solver-homotopy-warmstart.md @@ -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 diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 3f7ec1426..53a438fe7 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -2189,9 +2189,8 @@ class SolverBaseClass(uw_object): Parameters ---------- phase : str - ``"picard"`` -- ignore ``DIVERGED_MAX_IT`` (truncation is - intentional during the Picard-to-Newton transition). - ``"solve"`` -- any divergence is a real problem. + Label for the report. (A ``"picard"`` phase that ignored + ``DIVERGED_MAX_IT`` served the nrichardson warm-up removed in #791.) """ reason = self.snes.getConvergedReason() @@ -2199,10 +2198,6 @@ class SolverBaseClass(uw_object): if reason >= 0: return # converged -- nothing to report - # Picard truncation (DIVERGED_MAX_IT) is expected and harmless - if phase == "picard" and reason == -5: - return - # Bounded difficulty probe: hitting the iteration cap (DIVERGED_MAX_IT) is # the intended stop, not a failure — the effort is reported, not warned. if self._difficulty_probe and reason == -5: @@ -2325,8 +2320,44 @@ class SolverBaseClass(uw_object): # Stage 1 — Picard (alpha=0), loose tolerance. self._set_newton_alpha(0.0) - self.snes.setTolerances(rtol=max(self.newton_switch_rtol, rtol)) - self.snes.solve(None, gvec) + # An explicit solve(picard=N) makes the frozen stage run max(N, the natural + # stage-1 count) iterations, both measured against the ORIGINAL residual — the + # rotated path's rule (`rnorm <= switch_rtol*r0 and iters >= picard`, #791). + # ⚠️ PETSc measures rtol from EACH snes.solve call's own starting residual, so + # a second relative stage 1 after the N-iteration block would restart the + # clock and ask for a further `switch` reduction from the already-reduced + # residual (measured: picard=3 gave 3 + 20 frozen iterations, not max(3, 10)). + # Hence stage 1 after the block runs on an ABSOLUTE target derived from the + # block's initial residual, and returns at once if already below it. + n_min = int(getattr(self, "_continuation_min_picard", 0) or 0) + switch = max(self.newton_switch_rtol, rtol) + frozen_its = 0 + try: + if n_min > 0: + # respect a tighter cap already in force (estimate_difficulty probe) + cap = n_min if not max_it else min(n_min, int(max_it)) + self.snes.setTolerances(rtol=0.0, max_it=cap) + self.snes.solve(None, gvec) # ends DIVERGED_MAX_IT by design + frozen_its += int(self.snes.getIterationNumber()) + F0 = None + try: + hist = self.snes.getConvergenceHistory()[0] + F0 = float(hist[0]) if len(hist) else None + except Exception: + F0 = None + if F0: + self.snes.setTolerances(rtol=0.0, atol=max(atol, switch * F0), + stol=stol, max_it=max_it) + else: # no history available: fall back to the relative stage + self.snes.setTolerances(rtol=switch, atol=atol, stol=stol, + max_it=max_it) + else: + self.snes.setTolerances(rtol=switch, atol=atol, stol=stol, + max_it=max_it) + self.snes.solve(None, gvec) + frozen_its += int(self.snes.getIterationNumber()) + finally: + self.snes.setTolerances(rtol=rtol, atol=atol, stol=stol, max_it=max_it) if verbose and uw.mpi.rank == 0: print(f"continuation Picard: reason={self.snes.getConvergedReason()} " f"it={self.snes.getIterationNumber()}", flush=True) @@ -2335,6 +2366,11 @@ class SolverBaseClass(uw_object): self._set_newton_alpha(1.0) self.snes.setTolerances(rtol=rtol, atol=atol, stol=stol, max_it=max_it) # restore self.snes.solve(None, gvec) + # Per-stage iteration counts. solve_report only sees the last snes.solve + # (#778), so record the frozen-tangent stage explicitly. + self._continuation_stages = dict( + frozen_iterations=frozen_its, + newton_iterations=int(self.snes.getIterationNumber())) if verbose and uw.mpi.rank == 0: print(f"continuation Newton: reason={self.snes.getConvergedReason()} " f"it={self.snes.getIterationNumber()}", flush=True) @@ -9889,9 +9925,24 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): off stale data because a remesh or a diverged solve clears ``has_solution``. picard : int, default=0 - Number of Picard iterations before switching to Newton. - Picard iterations use a simplified Jacobian and can help - convergence for strongly nonlinear problems. + Minimum number of Picard iterations — Newton iterations with the + FROZEN (Picard) tangent, i.e. linear Stokes solves with the viscosity + held at the current state — before the consistent tangent takes over. + Meaning by :attr:`consistent_jacobian`: + + * ``False`` (default): every iteration already uses the frozen + tangent, so this is satisfied by the solve itself. + * ``"continuation"``: stage 1 (alpha = 0) runs at least ``picard`` + frozen-tangent iterations before the Newton stage. + * ``True``: raises ``NotImplementedError`` for a nonlinear residual — + the pure-Newton compile has no frozen tangent. Use + ``"continuation"`` for a Picard entry. (When the rest state is + exactly zero — homogeneous essential BCs, no stress history — + Newton's first iteration from rest already is the Picard step; a + boundary-driven or stress-history problem does not have that + property.) + + Negative values are treated as 0. verbose : bool, default=False Print solver progress and timing information. debug : bool, default=False @@ -9951,7 +10002,8 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): >>> velocity = stokes.u.array[:, 0, :] >>> pressure = stokes.p.array[:, 0, 0] - >>> # Nonlinear solve with Picard warmup + >>> # Nonlinear solve: 3 frozen-tangent iterations, then Newton + >>> stokes.consistent_jacobian = "continuation" >>> stokes.solve(picard=3) >>> # Time-stepping with previous solution @@ -9967,9 +10019,10 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): ----- This is a **collective operation** - all MPI ranks must call it. - For nonlinear viscosity (e.g., power-law, viscoplastic), the solver - uses Newton iteration by default. The ``picard`` parameter can help - with initial convergence by using simpler Jacobian approximations. + For nonlinear viscosity (e.g., power-law, viscoplastic) the tangent is + chosen by :attr:`consistent_jacobian` (frozen/Picard, consistent Newton, + or staged Picard->Newton), and ``picard`` sets a minimum number of + frozen-tangent iterations where that is meaningful (see above). See Also -------- @@ -10108,46 +10161,61 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): else: self.atol = 0.0 - # 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, contractive and cheap. The default - # (frozen) tangent path is left bit-identical, and an explicit Picard count - # is honoured as given. + # ⚠️ #791. A Picard step is a Newton iteration with the FROZEN (Picard) + # tangent: a linear Stokes solve — velocity block and Schur complement — + # with the viscosity held at the current state. From a cold start that IS + # the viscous solve. This block previously ran SNES `nrichardson` with no + # nonlinear preconditioner, i.e. x <- x - lambda F(x): a residual step with + # no linear solve, nearly inert, and NOT a Picard step. The semantics below + # mirror the rotated free-slip path (utilities/rotated_bc.py), which already + # had them right: # - # This warm-up is a CONVERGENCE aid (move a cold guess toward the Newton - # basin), NOT the correctness mechanism for the zero-strain-rate state. - # Two measured facts (issue #507) retired the old "make the NaN state - # unreachable" framing: (1) one nrichardson sweep only propagates - # boundary data a single element layer, so on a BOUNDARY-DRIVEN problem - # the deep interior stays exactly zero after the warm-up (body-force - # problems fill F(0) everywhere, which is why the yield campaigns never - # saw it); (2) a rigidly-translating stuck region has edot = 0 at the - # CONVERGED solution — the state is physics, not a start-up artifact. - # Finiteness of the consistent tangent at edot = 0 is owned by the - # half-integer-power guard in _jacobian_unwrap (the derivative's - # removable-singularity limit, implemented). NOTE the "continuation" - # tangent is NOT protected by its alpha = 0 phase (the blended kernel - # still evaluates the Newton branch pointwise, and IEEE 0*NaN = NaN); - # with the guard in place both tangents are finite everywhere. See - # docs/developer/design/nonlinear-solver-homotopy-warmstart.md (Layer 1). - if (picard == 0 and self.consistent_jacobian is True - and (zero_init_guess or self._solution_is_trivially_zero())): - picard = 1 + # consistent_jacobian False (default) — the whole solve already uses the + # frozen tangent, so `picard` warm-up iterations are inherently the + # first iterations of the solve itself. Nothing extra to run. + # "continuation" — stage 1 of _continuation_solve runs at alpha = 0, the + # frozen tangent; `picard` sets a minimum number of those iterations. + # True — the pure-Newton compile carries no frozen tangent. An explicit + # picard > 0 on a nonlinear residual RAISES (as the rotated path does) + # rather than silently running something else. The AUTOMATIC cold-start + # warm-up is removed, not repaired: it was never a Picard step, so + # nothing is lost. When the rest state is exactly zero (homogeneous + # essential BCs, no stress history) the strain rate is zero, the yield + # branch is inactive and Newton's first iteration IS the Picard step. + # ⚠️ That does NOT hold for a boundary-driven problem — the Dirichlet + # values put edot != 0 in the driven layer (measured: first steps 45% + # apart on a sheared box) — nor for a stress-history model, whose + # effective strain rate includes sigma*/(2 mu dt) at u = 0. A genuine + # Picard entry under the consistent tangent needs "continuation". + # + # Still true from the retired warm-up (#507): a zero strain rate is not only + # a start-up state — a rigidly translating stuck region has edot = 0 at the + # CONVERGED solution — so finiteness of the consistent tangent at edot = 0 is + # owned by the half-integer-power guard in _jacobian_unwrap, not by any + # warm-up. + # + # picard < 0 is accepted and treated as 0 (explicit "no warm-up"). + if picard < 0: + picard = 0 + self._continuation_min_picard = 0 + self._continuation_stages = None # never report a previous solve's stages + if picard > 0: + if self.consistent_jacobian == "continuation": + self._continuation_min_picard = int(picard) + elif self.consistent_jacobian is True: + if self._residual_is_nonlinear(): + raise NotImplementedError( + f"solve(picard={picard}) with consistent_jacobian=True: a Picard " + "warm-up needs the frozen (Picard) tangent, which the pure-Newton " + "compile does not carry. Use consistent_jacobian='continuation' " + "(staged Picard then Newton), consistent_jacobian=False (Picard " + "throughout), or drop `picard` to run pure Newton — whose first " + "iteration from a cold start already is the Picard step.") + picard = 0 # linear residual: the frozen tangent IS the tangent if verbose and uw.mpi.rank == 0: - print(f"SNES solve - picard = {picard}", flush=True) - - # Picard solves if requested - - if picard != 0: - self._push_snes_max_it(abs(picard)) - self._reassert_outer_tolerances() - self.snes.atol = self.atol - self.snes.setType("nrichardson") - self.snes.setFromOptions() - self._attach_stokes_nullspace() - self.snes.solve(None, gvec) - self._warn_on_divergence(phase="picard") + print(f"SNES solve - picard = {picard} " + f"(tangent: {self.consistent_jacobian!r})", flush=True) # The standard Newton solve, run whether or not the optional Picard # warmup above was taken: restore the configured SNES type and diff --git a/tests/test_0201_solver_has_solution_warmstart.py b/tests/test_0201_solver_has_solution_warmstart.py index 8ce6cc907..39ddbf84b 100644 --- a/tests/test_0201_solver_has_solution_warmstart.py +++ b/tests/test_0201_solver_has_solution_warmstart.py @@ -1,5 +1,5 @@ """Layer 1a of the nonlinear-solver warm-start / homotopy design: -``solver.has_solution`` and the automatic cold-start Picard warm-up. +``solver.has_solution`` and cold starts under the consistent-Newton tangent. Contract under test (design: ``docs/developer/design/nonlinear-solver-homotopy-warmstart.md``, Layer 1): @@ -11,9 +11,10 @@ cold-starts rather than warming off a stale iterate. * It survives a coefficient-only change (a new viscosity value), so parameter continuation and time-stepping warm-start correctly. - * The cold-start Picard warm-up runs under the consistent-Newton tangent - without breaking convergence, and leaves the default (frozen) tangent path - untouched. + * A cold start under the consistent-Newton tangent converges and sets the + flag. (The former automatic "Picard warm-up" on this path was an nrichardson + residual sweep and has been removed — #791; the Picard contract itself is + tested in test_1068.) These are cheap serial checks on the public API — the hard-case δ-continuation that motivates the design is validated separately against the Spiegelman study. @@ -186,8 +187,8 @@ def test_repeated_default_solve_agrees_to_solver_tolerance(): def test_cold_warmstart_under_consistent_newton_converges(): - """A cold consistent-Newton Stokes solve exercises the automatic - single-Picard warm-up branch — it must converge and set has_solution.""" + """A cold consistent-Newton Stokes solve must converge and set has_solution. + (No warm-up runs on this path any more — #791.)""" mesh = uw.meshing.StructuredQuadBox( elementRes=(8, 8), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0) ) @@ -202,6 +203,6 @@ def test_cold_warmstart_under_consistent_newton_converges(): stokes.add_essential_bc((0.0, None), "Right") stokes.consistent_jacobian = True - stokes.solve() # cold (zero_init_guess default True) → warm-up branch taken + stokes.solve() # cold: Newton straight from rest assert stokes.snes.getConvergedReason() > 0 assert stokes.has_solution is True diff --git a/tests/test_1068_picard_warmup_is_a_frozen_tangent_step.py b/tests/test_1068_picard_warmup_is_a_frozen_tangent_step.py new file mode 100644 index 000000000..ece44e80f --- /dev/null +++ b/tests/test_1068_picard_warmup_is_a_frozen_tangent_step.py @@ -0,0 +1,172 @@ +"""A Picard step is a Newton iteration with the FROZEN tangent — not a residual sweep (#791). + +The Layer 1 design (docs/developer/design/nonlinear-solver-homotopy-warmstart.md) specifies the +cold-start warm-up and ``solve(picard=N)`` as Picard (frozen-coefficient) iterations: linear +Stokes solves — velocity block and Schur complement — with the viscosity held at the current +state. The standard path instead ran SNES ``nrichardson`` with no nonlinear preconditioner, +``x <- x - lambda F(x)``: a residual step with no linear solve, nearly inert, and not a Picard +step. The rotated free-slip path already had the right semantics; these tests hold the standard +path to the same contract, per tangent mode: + + * ``consistent_jacobian=False`` — every iteration uses the frozen tangent, so ``picard`` is + satisfied by the solve itself and must change NOTHING (bit-identical result); + * ``"continuation"`` — ``picard=N`` guarantees at least N frozen-tangent iterations; + * ``True`` — the pure-Newton compile has no frozen tangent: ``picard>0`` raises on a nonlinear + residual, and is a no-op on a linear one. + +It also pins where the removed automatic warm-up is genuinely unnecessary: when the rest state is +exactly zero (homogeneous BCs, no stress history) the consistent and frozen tangents coincide, so +Newton's first iteration IS the Picard step. That does NOT extend to boundary-driven problems. +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + + +def _yielding_box(tag, tangent, tau_y=0.30, cellSize=0.25): + """Sheared hard enough to yield. The body force MUST vary in x: a uniform one is + hydrostatic, nothing moves and the yield law never engages.""" + uw.reset_default_model() + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=cellSize) + x, y = mesh.X + v = uw.discretisation.MeshVariable("Vpw" + tag, mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("Ppw" + tag, mesh, 1, degree=1, continuous=True) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + s.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel + cm = s.constitutive_model + cm.Parameters.shear_viscosity_0 = 1.0 + cm.Parameters.yield_stress = tau_y + s.bodyforce = sympy.Matrix([[0.0, -2.0 * sympy.cos(sympy.pi * x)]]) + s.add_essential_bc((sympy.oo, 0.0), "Top") + s.add_essential_bc((sympy.oo, 0.0), "Bottom") + s.add_essential_bc((0.0, sympy.oo), "Left") + s.add_essential_bc((0.0, sympy.oo), "Right") + s.petsc_use_pressure_nullspace = True + s.petsc_options.delValue("ksp_monitor") + s.consistent_jacobian = tangent + s.tolerance = 1.0e-6 + # The frozen tangent converges LINEARLY: at 1e-8 it exhausted the default 50-iteration + # cap on this fixture. Give it room, so a converged baseline exists to compare against. + s.petsc_options["snes_max_it"] = 300 + return s, v, p + + +def test_frozen_tangent_picard_changes_nothing(): + """Under consistent_jacobian=False the whole solve already takes Picard steps, so + picard=3 must reproduce picard=0 exactly — same iterations, same answer. The old + nrichardson sweeps moved the starting state and so changed the path.""" + s0, v0, _ = _yielding_box("f0", False) + s0.solve(zero_init_guess=True, picard=0) + r0 = s0.solve_report + ref = np.array(v0.data, copy=True) + + s3, v3, _ = _yielding_box("f3", False) + s3.solve(zero_init_guess=True, picard=3) + r3 = s3.solve_report + + assert str(r0.reason_str).startswith("CONVERGED"), r0.reason_str + assert r3.nl_its == r0.nl_its, ( + f"picard=3 took {r3.nl_its} Newton iterations vs {r0.nl_its} for picard=0 under the " + "frozen tangent — something other than frozen-tangent iterations ran first " + "(the #791 nrichardson sweep)") + # tight tolerance rather than bit-equality: the two solvers carry different variable + # names, and bit-identity across builds would depend on JIT term ordering (#752) + assert np.linalg.norm(np.array(v3.data) - ref) <= 1.0e-12 * np.linalg.norm(ref), ( + "picard=3 changed the frozen-tangent solution: the warm-up is not a no-op") + + +def test_continuation_picard_gives_max_of_n_and_the_natural_stage(): + """picard=N makes the alpha=0 (frozen-tangent) stage run max(N, natural) iterations, + both measured against the ORIGINAL residual — the rotated path's rule. Hard baselines: + + * N below the natural count -> exactly the natural count (picard adds nothing); + * N above it -> exactly N. + + Guards the restart trap: PETSc's rtol is relative to EACH snes.solve call's own + starting residual, so re-running a relative stage 1 after the N-iteration block + asked for a further reduction and gave N + (another whole stage 1).""" + s0, v0, _ = _yielding_box("c0", "continuation") + s0.solve(zero_init_guess=True) + natural = s0._continuation_stages["frozen_iterations"] + ref = np.array(v0.data, copy=True) + assert str(s0.solve_report.reason_str).startswith("CONVERGED") + assert natural >= 3, f"fixture too easy to discriminate: natural stage = {natural}" + + low = max(1, natural // 3) + s1, v1, _ = _yielding_box("cl", "continuation") + s1.solve(zero_init_guess=True, picard=low) + assert s1._continuation_stages["frozen_iterations"] == natural, ( + f"picard={low} (< natural {natural}) gave " + f"{s1._continuation_stages['frozen_iterations']} frozen iterations; expected " + f"exactly {natural}. More means stage 1 restarted its relative clock.") + + high = natural + 5 + s2, v2, _ = _yielding_box("ch", "continuation") + s2.solve(zero_init_guess=True, picard=high) + assert s2._continuation_stages["frozen_iterations"] == high, ( + f"picard={high} (> natural {natural}) gave " + f"{s2._continuation_stages['frozen_iterations']} frozen iterations; expected " + f"exactly {high}.") + for s, v in ((s1, v1), (s2, v2)): + assert str(s.solve_report.reason_str).startswith("CONVERGED"), s.solve_report.reason_str + assert np.linalg.norm(np.array(v.data) - ref) / np.linalg.norm(ref) < 1.0e-5 + + +def test_continuation_stages_do_not_go_stale(): + """A non-continuation solve after a continuation one must not report the old stages.""" + s, _, _ = _yielding_box("cs", "continuation") + s.solve(zero_init_guess=True) + assert s._continuation_stages is not None + s.consistent_jacobian = False + s.solve(zero_init_guess=True) + assert s._continuation_stages is None + + +def test_pure_newton_picard_raises_on_a_nonlinear_residual(): + """The pure-Newton compile carries no frozen tangent; silently running something else + (as the nrichardson sweep did) is replaced by a clear refusal, matching the rotated path.""" + s, _, _ = _yielding_box("n", True) + with pytest.raises(NotImplementedError, match="continuation"): + s.solve(zero_init_guess=True, picard=2) + + +def test_pure_newton_picard_is_ignored_on_a_linear_residual(): + """On a linear residual the frozen tangent IS the tangent, so picard is meaningless: + it must neither raise nor change anything (same iterations as picard=0).""" + runs = {} + for n in (0, 2): + s, v, _ = _yielding_box("l%d" % n, True, tau_y=1.0e6) # never yields + s.solve(zero_init_guess=True, picard=n) + assert str(s.solve_report.reason_str).startswith("CONVERGED") + runs[n] = (s.solve_report.nl_its, np.array(v.data, copy=True)) + assert runs[2][0] == runs[0][0] + assert np.linalg.norm(runs[2][1] - runs[0][1]) <= 1.0e-12 * np.linalg.norm(runs[0][1]) + + +def test_newton_first_step_from_rest_is_the_picard_step(): + """The claim that replaced the automatic warm-up. From a state of rest the strain rate is + zero, the yield branch is inactive, and the consistent tangent coincides with the frozen + one — so ONE Newton iteration from rest equals ONE Picard iteration. + + ⚠️ SCOPE: this holds when the rest state is exactly zero — this fixture is body-force + driven with homogeneous essential BCs. It does NOT hold for a boundary-driven problem: + the Dirichlet values yield the driven layer at once (measured: first steps 45% apart on + a sheared box), nor for a stress-history model. There, a Picard entry needs + consistent_jacobian="continuation" + picard=N.""" + states = {} + for tangent in (False, True): + s, v, _ = _yielding_box("s%s" % tangent, tangent) + s.petsc_options["snes_max_it"] = 1 + s.solve(zero_init_guess=True) + assert s.solve_report.nl_its == 1 + states[tangent] = np.array(v.data, copy=True) + rel = (np.linalg.norm(states[True] - states[False]) + / np.linalg.norm(states[False])) + assert rel < 1.0e-6, ( + f"first Newton step from rest differs from the Picard step by {rel:.3e} — the " + "consistent tangent at rest is NOT the frozen one, and a real warm-up is needed")