diff --git a/docs/developer/design/jacobian-consistent-tangent.md b/docs/developer/design/jacobian-consistent-tangent.md index 343c30dff..22f148986 100644 --- a/docs/developer/design/jacobian-consistent-tangent.md +++ b/docs/developer/design/jacobian-consistent-tangent.md @@ -2,6 +2,11 @@ **Status**: Implemented — PR #258 (2026-07-02): opt-in `solver.consistent_jacobian`, default off (Picard tangent unchanged). +**Since #794** the consistent tangent is the default (`consistent_jacobian = True`), and a +cold start takes one automatic Picard step first; `False` is the frozen tangent as an +explicit choice. The defaults named below are those of #258. See the `nonlinear-solver` +guide. + **Since #823 (tier 2)** the Newton source is no longer the `symbolic_keep_constants` expansion described below. `_jacobian_unwrap` replaces each non-constant atom by a *node* of the shared graph (`_jit_graph`), whose partial derivatives SymPy's chain rule diff --git a/docs/developer/design/jit-shared-graph-codegen.md b/docs/developer/design/jit-shared-graph-codegen.md index dd1cdac44..530db697a 100644 --- a/docs/developer/design/jit-shared-graph-codegen.md +++ b/docs/developer/design/jit-shared-graph-codegen.md @@ -225,7 +225,10 @@ $(\dot\varepsilon_{II}/\dot\varepsilon_0)^{1/n-1}$, is not merged, and both rout a finite Newton tangent at a state of rest, equal entry for entry. A power law on the bare named invariant, $\eta = \dot\varepsilon_{II}^{\,1/n-1}$, is merged into $g^{(1/n-1)/2}$: the tree's Newton flux is NaN at rest and the graph's is finite, because -$\dot\varepsilon_{II}$ stays a node whose body carries the guard (`test_0024`). +$\dot\varepsilon_{II}$ stays a node whose body carries the guard (`test_0024`). Since #841 +both guards cover every fractional power (a non-integer exponent on a base with free +symbols), so the merged power is guarded on both routes and both tangents are finite at +rest; `test_0024` keeps the old tree as its reference. ### Each lowering is read back by `getext()` from its nodes @@ -608,7 +611,7 @@ multiplied by a name that holds its reciprocal. | a node class from an earlier compile reused | a changed `.sym` is ignored | the compile serial in each class's identity; a test that changes a body and rebuilds in one process without clearing SymPy's cache | | source that depends on the hash seed, a creation counter or a re-declaration | every parallel run repairs; a new process or a re-run notebook cell misses the cache | canonical emission; `test_0105` extended with a Newton fixture, a preamble of extra objects and a re-declared law | | a cancellation across a name lost | NaN at a state where the tree gave a finite limit, in an unguarded residual or Picard kernel | every residual and Picard kernel compared with the tree at a state of rest on every fixture | -| `test_0022` (tier 1) asserts the expanded, guarded tree of `_jacobian_unwrap`, and checks a block for `Derivative` with `has()`, which does not see node bodies | one test fails by construction, the other cannot fail | the guard test now holds `_jit_graph.guard_half_integer_powers` to SymPy's `replace` on eight laws; the tests of deleted helpers went with them | +| `test_0022` (tier 1) asserts the expanded, guarded tree of `_jacobian_unwrap`, and checks a block for `Derivative` with `has()`, which does not see node bodies | one test fails by construction, the other cannot fail | the guard test now holds `_jit_graph.guard_fractional_powers` (named `guard_half_integer_powers` until #841) to SymPy's `replace` on eight laws; the tests of deleted helpers went with them | | an entry that is zero today becomes a non-zero expression | assembly of an entry that evaluates to zero | numbers and leaves inlined; zero patterns compared on every fixture | | guard placement changes cold-start behaviour | NaN at a state of rest | `test_0024`: the guarded lowering against a guarded tree, at a random state and at rest, for a viscoplastic law and a power law on a named invariant; `test_1067` (cold Newton start) | | `getext()` or another walker expands nodes back into the tree | the cost returns, silently | `getext()` lowers atoms itself; `test_0024` emits a twelve-layer law whose tree doubles with each layer and bounds its C | diff --git a/docs/developer/design/nonlinear-solver-homotopy-warmstart.md b/docs/developer/design/nonlinear-solver-homotopy-warmstart.md index 27d58e0ee..a012bd86d 100644 --- a/docs/developer/design/nonlinear-solver-homotopy-warmstart.md +++ b/docs/developer/design/nonlinear-solver-homotopy-warmstart.md @@ -66,6 +66,35 @@ every nonlinear solver; layers 2–3 build on it. ### Layer 1 — automatic warm-start (all nonlinear solvers) +```{warning} +**Implementation status (#791, 2026-09-26).** 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 implemented as specified: + +- `consistent_jacobian=True`: `solve(picard=N)` takes exactly N Picard steps, then Newton, and a + cold start takes one automatically (`picard=-1` disables it). The Jacobian is compiled in the + blended form J_p + alpha (J_n - J_p); alpha = 0 is the Picard step, alpha = 1 Newton, no + recompile to switch. +- `"continuation"`: ONE solve with alpha keyed on the residual. alpha leaves 0 at + F = 0.1 F_ref and reaches 1 at F = 5e-4 F_ref, linear in log F, never decreasing, with + F_ref = ||F(u=0)|| (one residual evaluation). `l2` line search for that solve; stop at + `tolerance` x F_ref. `picard=N` holds alpha = 0 for the first N iterations. A warm start + enters the ramp at the alpha its residual calls for — keyed on the solve's own starting + residual instead, a warm start re-ran the whole Picard phase (763 vs 10 iterations on the + hard notch). Measured on the notch and a power-law variant, this converges cases where pure + Newton and the former two-stage continuation both fail (eta_bg 5e24; power law n = 5); it + still stalls where the blended iteration oscillates at intermediate alpha (refinement 2; + power law n = 3 with yield). The rotated-free-slip path keeps its own two-phase loop. +- `False`: every iteration is already a Picard step. + +Measured on a boundary-driven sheared box: without the warm-up, pure Newton stalls on step size +(CONVERGED_SNORM_RELATIVE at ||F|| = 4.5e-5, whatever the tolerance); with it, Newton converges on +the residual to machine precision in fewer iterations (`tests/test_1068_picard_warmup_is_a_frozen_tangent_step.py`). +Measurements quoted elsewhere for "opening Picard steps" before 2026-09-26 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 @@ -168,9 +197,15 @@ domain in the softened region. The cost of `"yield"` is stiffened unyielded mate bounded by 2 for the sqrt family and by $2^{\delta}$ for the power mean; that bound is why the power mean's entry δ is O(1) while the sqrt family's is O(10). -The default is `"onset"` because it is the historical law, not because it is the better -one. Changing it alters results for every existing `yield_mode="softmin"` user, so it is -a separate maintainer decision. +The default is `"onset"`, and since 2026-10-05 it is also the default for the default +`"softmin"` law. `"yield"` was tried as that default and refused in review (PR #794): +for the sqrt family the stiffening is not confined to the corner — every unyielded cell +is $1/(1-\delta/2)$ stiffer (5 % at δ = 0.1, so a model that never yields changes its +answer) — and above yield the stress falls back to $\tau_y$ from $1.033\,\tau_y$, a +negative yielded tangent (−0.029 of the viscous slope). Under `"onset"` both limits are +exact and the yielded tangent is $\ge -10^{-4}$ of the viscous one. The approach from +above that `"yield"` was meant to provide is the declared `plastic_rate_strengthening` +ladder (`solve(homotopy="rate_strengthening")`). ### Layer 3 — smoother as a consistent-Newton consequence diff --git a/docs/developer/guides/nonlinear-solver.md b/docs/developer/guides/nonlinear-solver.md index 1674d0bf7..f275274df 100644 --- a/docs/developer/guides/nonlinear-solver.md +++ b/docs/developer/guides/nonlinear-solver.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 warm-up (`solve(picard=N)` gives N real Picard steps then Newton, and one is automatic 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`. families: [Stokes, ViscoPlasticFlowModel] kind: recipe --- @@ -22,9 +22,11 @@ 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, and take a few **Picard steps** into the Newton + basin: `solve(picard=N)` — N successive-substitution steps (linear Stokes solves with + the viscosity frozen), then Newton. A cold start takes one automatically. From a + *warm* iterate, take **no** Picard step (it wastes the good quadratic start). + ⚠️ #791: before 2026-09-26 the "Picard step" was an nrichardson residual sweep. 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). @@ -47,7 +49,7 @@ 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 its Picard warm-up (`solve(picard=N)`) 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. @@ -79,7 +81,8 @@ case felt like whack-a-mole. Check them first. | Trap | Symptom | Fix | |---|---|---| -| **Perfect plasticity's consistent tangent is SINGULAR along the flow**: on the hard-`Min` plastic branch η = τ_y/2ε̇_II, so 2η + 2η′ε̇_II = 0 — the velocity block is symmetric but semi-definite in every yielded cell. (An earlier version of this row blamed *asymmetry*; that is wrong for any η(ε̇_II) law — the rank-one term η′ ε̇⊗ε̇/ε̇_II is symmetric. Pressure-dependent yield adds a non-symmetric v–p coupling, not a non-symmetric velocity block. Corrected 2026-08-26, maintainer review.) | benign while yielded cells are few (the viscous neighbours regularise); with a large yielded fraction the velocity sub-solve caps out and Newton stalls at ~1e-3, no failure reason | give the plastic branch a positive tangent: a small δ soft-min (`yield_mode="softmin"`, powermean, `yield_anchor="yield"`), a rounded viscosity floor, or rate-strengthening ξ; Picard converges regardless (full 2η stiffness) but is linear-rate. The FMG bundle's `gmres`+`sor` smoother is Newton-safe either way | +| **A pressure-dependent yield stress makes the Newton Schur complement INDEFINITE far from the solution**: the consistent tangent differentiates η by p, so J_vp is no longer −Bᵀ (Spiegelman notch at η_bg 5e24: 41% off, 40 of 1968 Schur eigenvalues change sign). In the basin of convergence the full tangent is benign again (no wrong-sign eigenvalues, no negative pivots). | fieldsplit Schur diverges at ~0.99 per iteration even with the velocity block solved exactly; `DIVERGED_LINEAR_SOLVE` on the first Newton steps of a cold Drucker-Prager solve | `newton_pressure_coupling=False` (the Newton velocity block with the frozen pressure column, a partial tangent with the standard saddle-point structure) on the approach, or `consistent_jacobian="continuation"`; the default keeps the full tangent, which is right once inside the basin | +| **Perfect plasticity's consistent tangent is SINGULAR along the flow**: on the hard-`Min` plastic branch η = τ_y/2ε̇_II, so 2η + 2η′ε̇_II = 0 — the velocity block is symmetric but semi-definite in every yielded cell. (An earlier version of this row blamed *asymmetry*; that is wrong for any η(ε̇_II) law — the rank-one term η′ ε̇⊗ε̇/ε̇_II is symmetric. Pressure-dependent yield adds a non-symmetric v–p coupling, not a non-symmetric velocity block. Corrected 2026-08-26, maintainer review.) | benign while yielded cells are few (the viscous neighbours regularise); with a large yielded fraction the velocity sub-solve caps out and Newton stalls at ~1e-3, no failure reason | give the plastic branch a positive tangent: the declared `Parameters.plastic_rate_strengthening` (η_reg adds 2η_reg along the flow in every yielded cell; its scale is the model's), a rounded viscosity floor (clips only the fastest cells). The default soft-min (sqrt, `"onset"`, δ = 0.1) does NOT do it — its yielded tangent is ≥ −1e-4 of the viscous one, still effectively singular; the sqrt law under `yield_anchor="yield"` is negative (−0.029). Picard converges regardless (full 2η stiffness) but is linear-rate. The FMG bundle's `gmres`+`sor` smoother is Newton-safe either way | | `preconditioner="fmg"` (vs explicit `pc_type=mg` + manual mg opts) | outer KSP "converges" in **1 iteration** → no real Newton correction → stall → `DIVERGED_LINE_SEARCH` | use explicit `pc_type=mg` with the smoother opts above; bound the outer KSP (`ksp_max_it`~80) so a hostile step fails fast | | Cold plastic start `v=0`, or any rigid/unyielded point | `DIVERGED_FNORM_NAN` at iteration 0 | **Not** a div/0: `ε̇=0` gives `η_pl=+inf`, which `Min` and the sqrt soft-min carry correctly to the viscous branch. Only a soft-min form that computes `η_ve·η_pl/(η_ve+η_pl)` breaks (`inf/inf`). Fixed in the power-mean; if you hand-roll a blend, write the harmonic mean as `η_ve/(1+η_ve/η_pl)`. **Do not reach for a strain-rate floor** — it hides this rather than fixing it | | LU velocity block with all-Dirichlet-ish BC | pressure nullspace singular | attach the Stokes nullspace / avoid a bare LU there | @@ -100,23 +103,71 @@ 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. +`consistent_jacobian=True` (the consistent Newton tangent) is the DEFAULT for every +solver class since 2026-10-05. On a **cold** solve ONE Picard step is taken +automatically (a real one — #791 replaced the former nrichardson sweep), in Stokes and +in the Scalar, Vector and MultiComponent solvers alike. No automatic step is taken for a +`ksponly` (linear) solve, for a flux whose coefficient does not depend on the unknowns +(the two tangents coincide), or on the rotated free-slip and fault-contact paths, which +run their own loops; an explicit `picard=N` on the rotated path stages N Picard +iterations then Newton on the same blended Jacobian. A warm start takes none. ```python -stokes.consistent_jacobian = True -stokes.solve() # cold → one automatic Picard step, then Newton +stokes.solve() # cold -> one Picard step, then Newton +stokes.solve(picard=3) # the recipe: three Picard steps, then Newton +stokes.solve(picard=-1) # no warm-up if stokes.has_solution: ... ``` +The warm-up runs inside the solve's own contracts: update callbacks, the +`guard(wall_per_step=...)` deadline, an `estimate_difficulty` cap (the steps count +against it), and the report. `solve_report.nl_its`, `ksp_its` and the history include +the warm-up, the reduction is measured from the cold residual, and +`solve_report.config["tangent"]["picard_warmup_its"]` says how many steps were taken. +The Newton stage that follows stops at `tolerance`·‖F(u_cold)‖, where one uninterrupted +solve would. + +Changing `consistent_jacobian` or `newton_pressure_coupling` after a solve recompiles the +Jacobian on the next one; the report always describes the tangent that ran. + +Under `False` every iteration is already a Picard step. Under `"continuation"` the whole +Picard→Newton passage is ONE solve: α is set each iteration from the residual relative to +‖F(u=0)‖ (0 until 0.1, 1 from 5e-4, linear in log F between, never decreasing), with an +`l2` line search and a stop at `tolerance`·‖F(u=0)‖; `picard=N` holds α = 0 for the first N +iterations. Tunables: `continuation_alpha_start`, `continuation_newton_at`, +`continuation_linesearch`. This is the robust choice for a hard cold start AND for warm +continuation steps (δ, η, BC changes) — a warm start enters the ramp near α = 1. + +### What the solve records (`solve_report.config`) + +Every solve records the configuration that actually ran, read from the live PETSc +objects at the call: SNES type and tolerances, line search, outer KSP (with its Krylov +`restart`) and the fieldsplit sub-solvers, the tangent (`consistent_jacobian`, +`newton_pressure_coupling`, α where the blended Jacobian is live), and for Stokes the +pressure `gauge` (nullspace requested / attached, whether the rheology reads p, pressure +Dirichlet BCs, natural boundaries) and the yield-law `regularisation` (mode, smoother, +anchor, δ, `plastic_rate_strengthening` and its ladder multiplier). The step transcript's +solve event carries `tangent`, `gauge` and `regularisation`, so a run can be debugged from +the record without a live process. + +Two warnings come from the same facts, once per solve: a mean-pressure nullspace with a +viscosity that reads p (the yield law then sees a gauge-fixed pressure), and a +pressure-dependent yield stress under the consistent tangent with no stated +`plastic_rate_strengthening` (see `plasticity-solvers`). + +The Krylov restart is 100 by default (outer KSP of every solver, and the Stokes velocity +split), not PETSc's 30: a yielded layer gives the Newton velocity block a near-null +mechanism mode that FGMRES(30) cannot resolve at all (rate 1.000) and FGMRES(100) does. +A user's own `ksp_gmres_restart` is left alone. + --- ## 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 (corrected by #791):** `has_solution` + cold consistent-Newton + Picard warm-up, now a real Picard step (blended kernel at α=0); `solve(picard=N)` for + N steps (tests `test_0201`, `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. @@ -129,7 +180,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 + Picard warm-up and grid sequencing have failed. --- diff --git a/docs/developer/guides/plasticity-solvers.md b/docs/developer/guides/plasticity-solvers.md index eab59b7e7..b32328b7d 100644 --- a/docs/developer/guides/plasticity-solvers.md +++ b/docs/developer/guides/plasticity-solvers.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) warmed by real Picard steps via solve(picard=N) (one is automatic on a cold start), 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`. families: [Stokes, ViscoPlasticFlowModel, ViscoElasticPlasticFlowModel] kind: recipe --- @@ -18,7 +18,8 @@ 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() # cold: ONE Picard step automatically, then Newton +stokes.solve(picard=3) # the recipe: a few Picard steps, then Newton ``` --- @@ -28,20 +29,25 @@ 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. **"A few Picard steps, then Newton" is one flag: `solve(picard=N)`.** Under + `consistent_jacobian=True` it takes exactly N Picard steps — successive substitution, + a linear Stokes solve with the viscosity frozen at the current state — then Newton. + A cold start takes ONE automatically; `picard=-1` switches it off. It matters most for + BOUNDARY-DRIVEN problems: measured on a sheared box, pure Newton from rest stalls on + step size (`CONVERGED_SNORM_RELATIVE` at ||F|| = 4.5e-5 — reported as "converged"), + while one or three Picard steps first give true residual convergence in fewer + iterations. (Mechanism: the Jacobian is compiled in the blended form + J_p + alpha (J_n - J_p); alpha = 0 is the Picard step, alpha = 1 is Newton, and + switching needs no recompile.) ⚠️ #791: the former warm-up was an nrichardson + residual sweep, not a Picard step; measurements quoted for "opening Picard steps" + before 2026-09-26 were of that sweep and are 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: `solve(picard=N)` (N Picard steps, then Newton), or + `consistent_jacobian="continuation"` (one solve, α-blend keyed on the residual + relative to ‖F(u=0)‖, α a `constants[]` atom, no recompile), or `consistent_jacobian=False` for Picard throughout. Rescue on *stagnation*, not only on a failure reason — the failure-only trigger measured byte-identical to doing nothing at the cliff. @@ -55,7 +61,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. --- @@ -84,14 +90,16 @@ rescue, and grid sequencing are. ## Which tangent for which model (measured) -`solver.consistent_jacobian` takes `False` | `True` | `"continuation"`: +`solver.consistent_jacobian` takes `False` | `True` | `"continuation"`; the default is +`True` for every solver class since 2026-10-05, so a model that wants the frozen tangent +must ask for it: | Model | Use | Why | |-------|-----|-----| -| `ViscoPlasticFlowModel` (non-elastic) | **`True`** (Newton) | Quadratic near the solution; the automatic Picard entry handles the cold start. | -| `ViscoElasticPlasticFlowModel` (VEP) | **`False`** (Picard) | The consistent yield tangent over the elastic stress-history block makes the Jacobian **indefinite → `DIVERGED_LINEAR_SOLVE`**. Picard is contractive. | +| `ViscoPlasticFlowModel` (non-elastic) | **`True`** (Newton) | Quadratic near the solution; a cold start takes one Picard step automatically — use `solve(picard=N)` for more. | +| `ViscoElasticPlasticFlowModel` (VEP) | **`False`** (Picard), set explicitly | The consistent yield tangent over the elastic stress-history block makes the Jacobian **indefinite → `DIVERGED_LINEAR_SOLVE`**. Picard is contractive. (The review of PR #794 found `test_1052`'s loading-through-yield converging under the Newton default with the smooth law, nl = 1 per step; the measurement below predates the smooth default.) | | `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. | +| Any, far from the solution | **`"continuation"`** | One solve: α (0 Picard → 1 Newton) set each iteration from F/‖F(u=0)‖ — leaves 0 at 0.1, reaches 1 at 5e-4 — with an `l2` line search (`bt` stops the Picard end). Warm starts enter at the α their residual calls for. Converged the hard notch (η_bg 5e24, also at δ 0.1) and a power law n = 5 where pure Newton fails; 3–9x fewer KSP than cold Newton where both converge. Known stall: the blended iteration oscillates at α ≈ 0.4–0.6 (notch at refinement 2; power law n = 3 with yield). | > Measured: VEP loading-through-yield — Picard converges (σ locks at τ_y), > Newton diverges every step (`DIVERGED_LINEAR_SOLVE`). @@ -155,8 +163,11 @@ If you want a *rounded* yield law at all (as physics or as a formulation choice) the substrate is three model properties; δ is a `constants[]` atom, so changing it never recompiles: -- **`yield_mode`**: `"min"` (default — exact hard `Min`), `"softmin"` (the - δ-parameterised family below), `"harmonic"` (a **distinct physical model**, a +- **`yield_mode`**: `"softmin"`, the δ-parameterised family below (DEFAULT since + 2026-10-05: `sqrt` smoother, `"onset"` anchor, δ = 0.1 — the viscous branch and + the plastic limit exact, the yield point itself 4.5 % under-stressed, the yielded + tangent ≥ −1e-4 of the viscous one); `"min"` (exact hard `Min`, an explicit choice, + the default until 2026-10-05); `"harmonic"` (a **distinct physical model**, a parallel blend — not an approximation to `Min`). - **`yield_smoother`**: `"sqrt"` or `"powermean"`. **δ is NOT the same parameter in the two families**: the power mean's sharpness is `s = 1/(δ + 0.001)`, so @@ -165,11 +176,15 @@ never recompiles: δ = 0 lands within 0.07 % of `Min` — an order of magnitude inside a 1e-8 solver tolerance. - **`yield_anchor`**: which point is pinned to the exact law — the SIDE of `Min` - belongs to the anchor, not the family. `"onset"` (default, historical) is exact - on the unyielded branch but sits BELOW `Min` at and above yield — a *weaker* - problem than the sharp one. `"yield"` pins τ/τ_y = 1 exactly and sits on-or-above - `Min` everywhere; the cost is stiffer unyielded material (bounded ×2 sqrt, - ×2^δ powermean, both → 1 as δ → 0). + belongs to the anchor, not the family. `"onset"` (default) is exact on the + unyielded branch but sits BELOW `Min` at and above yield — a *weaker* problem than + the sharp one, by 4.5 % at the yield point at δ = 0.1 and by a third at δ = 64. + `"yield"` pins τ/τ_y = 1 exactly and sits on-or-above `Min` everywhere; the cost is + stiffer unyielded material EVERYWHERE below yield, not just near it (×1/(1 − δ/2) + for sqrt — 5 % at δ = 0.1; ×2^δ for powermean), and for sqrt a stress that falls back + to τ_y from 1.033 τ_y above yield, i.e. a NEGATIVE yielded tangent (−0.029 of the + viscous slope at δ = 0.1). A δ-march entry, not a single-solve law; it was refused + as the default in review (PR #794). **If you march δ toward the sharp law, the only sound discipline is multi-solve:** hold δ constant for a full solve to tolerance, warm-start the next smaller δ, @@ -184,8 +199,43 @@ measured, it never has. --- +## Declared regularisation: `plastic_rate_strengthening` + +A yielded cell's consistent tangent has no modulus along its own strain-rate direction +(perfect plasticity), so a yielded layer is a near-null mechanism mode of the velocity +block (Spiegelman notch: Rayleigh quotient 1.8e-7). The model can state a +regularisation for it: + +```python +cm.Parameters.plastic_rate_strengthening = 1.0e19 # eta_reg, units of viscosity +``` + +The law becomes `eta_pl -> tau_y / (2 edot_II) + m * eta_reg` on ViscoPlastic, VEP and +TI-VEP. It floors the plastic viscosity AND adds `2 m eta_reg` to the tangent along the +flow in every yielded cell. Its scale is the model's: there is no universal +dimensionless value (a 1 % overstress on the notch is 20 x sigma_y in a convective lid), +which is why the default is 0 (no term) and a pressure-dependent yield stress under the +consistent tangent without it warns. `m` (`rate_strengthening_scale`, default 1) is a +`constants[]` multiplier, so changing it is a repack, not a recompile. + +`stokes.solve(homotopy="rate_strengthening")` ladders `m` from the viscous limit down to +exactly 1 (ViscoPlastic only — the elastic models would advance their stress history on +every rung): one solve at an entry multiplier `1e3 x max(eta_ve / eta_reg)` at which +nothing yields, then the largest `m` at which a cell still yields there, read from the +fields, then rungs down to `m = 1`, each warm-started, failed rungs retried from the last +converged one. A ladder that stops above 1 leaves the model at `m = 1`, the fields at +the last converged rung, and warns. The summary carries `scale0`, `rungs`, +`reached_one`, `converged`. Measured on the notch: the equivalent ladder keeps the +velocity block healthy on every rung and reaches a converged state no single solve +could. + ## Floors +- With `plastic_rate_strengthening` stated, the viscosity floor is redundant on the + plastic branch wherever `m * eta_reg >= shear_viscosity_min`; the floor remains a + safeguard for the creep branch. `yield_stress_min` (the POSITIVE residual strength in + tension) is a statement of the yield law and stays. + - **`shear_viscosity_min`** (default `-oo` = off) is applied through `uw.maths.smooth_max`, but the default rounding scale is zero under `yield_mode="min"` and `δ·|floor|` under the smooth modes — so it **vanishes as @@ -208,8 +258,8 @@ 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 | -| 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 | +| Newton stalls with no divergence reason | admissible uselessness — steps accepted, residual flat | revert to best iterate, Picard block (`solve(picard=N)` or `"continuation"`); consider grid sequencing | +| Converges but σ sits **below** τ_y | a fixed δ>0 soft-min under the default `"onset"` anchor is a WEAKER law (4.5 % at the yield point at δ = 0.1) | that is the modelling choice you made — a smaller δ, or `yield_mode="min"` for the exact surface (`yield_anchor="yield"` reaches τ_y but stiffens every unyielded cell and softens above yield) | | 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/reviews/2026-10/nonlinear-defaults-794-adversarial-review.md b/docs/reviews/2026-10/nonlinear-defaults-794-adversarial-review.md new file mode 100644 index 000000000..b653fd1ac --- /dev/null +++ b/docs/reviews/2026-10/nonlinear-defaults-794-adversarial-review.md @@ -0,0 +1,122 @@ +# 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 a8aa68c0): 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 8010d335. + +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 (8010d335, 6535935b) 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 8010d335) + +- **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. + +## Round 3 — the merge with development (2c6bc2f9: #823 tier 2, #835) + +The textual conflicts were five hunks in three files: `_jacobian_source` takes the +solver's JIT route and keeps the warm-up blend, the `consistent_jacobian` docstring keeps +this branch's defaults with tier 2's description of the Newton source, `Stokes.solve` +keeps `_skip_stress_history` and the untyped `homotopy`, and `test_0104` keeps the frozen +tangent and the route. The review then attacked the seams between this branch's tangent +machinery and the two JIT routes. + +- **blocker, fixed (#841, 4cdc00dd).** 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)`, which the + half-integer guard missed, so `J(0)` was NaN and the blended warm-up kernel inherited it + at alpha = 0. The guard now covers every non-integer exponent on a base with free + symbols, on both routes. `test_1067` raised 77 on both routes before the change. +- **should-fix, fixed (bf4c32ff).** The rate-strengthening probe projection compiled on + the process default route, so `jit_route = "expanded"` did not rule the graph out of the + ladder's entry multiplier. It now takes its solver's route. +- **should-fix, fixed (bf4c32ff).** Docstrings: `True` returns the blend once a cold start + has planned its Picard step; alpha = 0 is bit-identical to Picard only while the Newton + blocks are finite. The consistent-tangent design note said "default off" (3eaa99be). +- **open.** `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}}`. The rate-strengthening ladder ends at + the stated problem; the march does not. +- **open.** On the DP box the automatic warm-up doubles the iteration count (19–20 against + 9 with `picard=-1`): after the Picard step the backtracking 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. +- **open, #827.** The PETSc banner from `_check_velocity_preconditioner` also prints after + a rewire followed by a warm solve that converges at iteration 0 (a tangent toggle, a + route switch, the tangent restored after either homotopy). Triggers added to #827. + +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, and 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. On the DP box with the default law, the default tangent takes 20 +iterations on both routes and the solutions agree to 2.2e-16 (v) and 7e-16 (p); `J(0)` is +finite and equal on both routes. `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. + +## Open (not fixed on this branch) + +- **#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. diff --git a/docs/reviews/README.md b/docs/reviews/README.md index e47f7cfe8..4ed31e180 100644 --- a/docs/reviews/README.md +++ b/docs/reviews/README.md @@ -262,6 +262,7 @@ docs/reviews/ | Review | Component | Tests | Pass Rate | Date | Status | |--------|-----------|-------|-----------|------|--------| +| Nonlinear defaults and η_reg ([2026-10](2026-10/nonlinear-defaults-794-adversarial-review.md), PR #794) | Solvers/Constitutive | adversarial, 2 rounds | Passing | 2026-10-06 | Fixes on branch | | JIT setup cost ([2026-10](2026-10/jit-setup-cost-823-adversarial-review.md), #823) | JIT/Solvers | adversarial, 3 rounds | Passing | 2026-10-07 | Fixes on branch | | Units System (2026-02) | Units/Scaling | 60+ | Passing | 2026-02-01 | Approved (2026-03-10) | | Data Access (2026-02) | Array/Math Interface | 75+ | Passing | 2026-02-01 | Approved (2026-03-10) | diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 94b8b3595..0b44198e2 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -54,6 +54,39 @@ from underworld3.function import expression as public_expression + +def _parameter_is_zero(par): + """True when a constitutive parameter is exactly zero: ``0``, ``0.0``, a units + quantity of zero, or an expression that unwraps to zero. ``None`` counts as zero. + + The one test every "is this term stated?" decision uses. ``sympy.Float(0.0) == 0`` + is False (SymPy >= 1.13), so ``== 0`` / ``!= 0`` on a parameter's ``.sym`` gives + different answers for ``0`` and ``0.0``; ``is_zero`` does not. + """ + if par is None: + return True + from underworld3.function.expressions import unwrap_expression + try: + value = unwrap_expression(par, mode="nondimensional") + value = getattr(value, "magnitude", value) + return sympy.sympify(value).is_zero is True + except (TypeError, ValueError, AttributeError, sympy.SympifyError): + # an expression we cannot reduce to a number is a stated (non-zero) term + return False + + +def _plastic_rate_strengthening_parameter(): + """The declared plastic rate strengthening, shared by every yielding model's + Parameters (ViscoPlastic, ViscoElasticPlastic, TransverseIsotropicVEP).""" + import underworld3.utilities._api_tools as api_tools + return api_tools.Parameter( + R"{\eta_{\mathrm{reg}}}", + lambda inner_self: 0, + "Plastic rate strengthening: eta_pl -> tau_y/(2 edot_II) + m eta_reg (a stated " + "regularisation; floors the plastic viscosity AND its tangent modulus; 0 = no term)", + units="Pa*s", + ) + expression = lambda *x, **X: public_expression(*x, _unique_name_generation=True, **X) @@ -1091,6 +1124,84 @@ def _yield_offset_for_anchor(self): return delta / 2 return (-1 + sympy.sqrt(1 + delta**2)) / 2 + def _with_rate_strengthening(self, viscosity_yield, parameters): + """``viscosity_yield + m * eta_reg`` when the model states a plastic rate + strengthening, else ``viscosity_yield`` unchanged (no term compiled).""" + eta_reg = getattr(parameters, "plastic_rate_strengthening", None) + if _parameter_is_zero(eta_reg): + return viscosity_yield + return viscosity_yield + self._get_rate_strengthening_scale() * eta_reg + + def _get_rate_strengthening_scale(self): + r"""The homotopy multiplier on ``plastic_rate_strengthening`` as a ``constants[]`` atom. + + The law is ``eta_pl -> tau_y/(2 edot_II) + m * eta_reg`` with ``m = 1`` the stated + problem. A ladder in ``m`` from a large value (effectively viscous) down to 1 is an + approach from above that ends at the problem as stated; ``m`` is a rampable + constant, so each rung is a repack and not a recompile. The + regularisation is a DECLARED parameter with its scale stated per model — there is + no universal dimensionless value. + """ + value = getattr(self, "_rate_strengthening_scale", 1.0) + if getattr(self, "_rate_strengthening_scale_expr", None) is None: + self._rate_strengthening_scale_expr = expression( + R"{m_{\mathrm{reg}}}", + sympy.Float(value), + "Homotopy multiplier on the plastic rate strengthening (1 = the stated problem)", + ) + else: + self._rate_strengthening_scale_expr.sym = sympy.Float(value) + return self._rate_strengthening_scale_expr + + @property + def rate_strengthening_scale(self): + """Multiplier ``m`` on ``Parameters.plastic_rate_strengthening``; 1 = as stated.""" + return getattr(self, "_rate_strengthening_scale", 1.0) + + @rate_strengthening_scale.setter + def rate_strengthening_scale(self, value): + # No _reset(): m is a constants[] atom and the next solve repacks it. A reset + # would force a JIT rewire on every rung (measured 5.2 s against 0.23 s per rung + # on a 0.25 box). Before the atom exists the value is read when it is created. + self._rate_strengthening_scale = float(value) + if getattr(self, "_rate_strengthening_scale_expr", None) is not None: + self._rate_strengthening_scale_expr.sym = sympy.Float(self._rate_strengthening_scale) + + def _rate_strengthening_control(self): + """How the solver ladders the declared rate strengthening (``m`` from the viscous + limit down to 1). ``None`` when the model states no ``plastic_rate_strengthening``.""" + from underworld3.systems.yield_continuation import RateStrengtheningControl + eta_reg = getattr(self.Parameters, "plastic_rate_strengthening", None) + if _parameter_is_zero(eta_reg): + return None + + def set_scale(value): + self.rate_strengthening_scale = value + + return RateStrengtheningControl( + set_scale=set_scale, + tangent=self._yield_homotopy_tangent, + scale=self._get_rate_strengthening_scale(), + parameter=eta_reg, + probe_terms=self._rate_strengthening_probe_terms(), + ) + + def _rate_strengthening_probe_terms(self): + r"""Cell-wise terms the rate-strengthening ladder reads to place its rungs, or + ``None`` when this model does not provide them (the ladder then needs an + explicit ``scale0``). + + Returns a dict of expressions: ``viscous`` and ``eta_reg`` (their ratio, at any + state, sets an entry ``m`` at which nothing yields), and ``onset_num`` / + ``onset_den``, :math:`2\dot\varepsilon_{II}\eta_{ve} - \tau_y` and + :math:`2\dot\varepsilon_{II}\eta_{reg}`, whose ratio is the ``m`` at which a cell + just stops yielding at the current fields. Kept as a numerator/denominator pair + so a rigid cell (zero strain rate) or a material with no regularisation (zero + ``eta_reg``) gives finite values that the ladder masks, instead of an infinite + or undefined quotient inside the kernel. + """ + return None + def _get_yield_offset(self): """The offset constant atom (lazily created alongside δ).""" if getattr(self, "_yield_offset_expr", None) is None: @@ -1114,9 +1225,15 @@ def yield_anchor(self): :math:`f = \eta_{ve} / \eta_{pl}` for the overstress ratio (so :math:`f = 1` is the yield point and :math:`f < 1` is unyielded): - ``"onset"`` (default, historical) + ``"onset"`` (default) Pins the unyielded limit :math:`f \to 0`, where :math:`\eta = \eta_{ve}` - exactly. The curve then sits BELOW the exact law at and above the yield + exactly, so a material that never yields gives the viscous answer. At the + default δ = 0.1 the sqrt law under-stresses the yield point by 4.5% + (:math:`\tau/\tau_y = 0.955` at :math:`f = 1`), recovers the plastic limit + :math:`\tau \to \tau_y` from within 0.04%, and its tangent above yield is + at least :math:`-10^{-4}` of the viscous one — perfect plasticity to within + rounding, made positive by any ``plastic_rate_strengthening``. For a large + homotopy δ the curve sits well BELOW the exact law at and above the yield point — the sqrt family undershoots for :math:`f < 2` (at :math:`\delta = 64`, :math:`\eta/\eta_{ve} = 0.80` at :math:`f = 0.5` where the exact law gives 1.0, and the yield point itself is a third @@ -1138,7 +1255,13 @@ def yield_anchor(self): arguments are equal. The cost is stiffened unyielded material, bounded by a factor 2 for the sqrt family and by :math:`2^{\delta}` for the power mean, both decaying to 1 as :math:`\delta \to 0`. The power-mean bound is why its - entry δ is O(1) and not O(10). + entry δ is O(1) and not O(10). The stiffening is not confined to the corner: + the sqrt law is :math:`1/(1 - \delta/2)` times the viscous one for EVERY + :math:`f \ll 1` (5% at δ = 0.1), and above yield its stress falls back to + :math:`\tau_y` from :math:`1.033\,\tau_y` at :math:`f = 1.2`, so the + yielded tangent is NEGATIVE (−0.029 of the viscous slope at δ = 0.1); the + power mean's plastic limit is :math:`2^{\delta}\tau_y`, not :math:`\tau_y`. + Use it for a δ march's entry, not as a single-solve law. """ return getattr(self, "_yield_anchor", "onset") @@ -1475,10 +1598,16 @@ class ViscoPlasticFlowModel(ViscousFlowModel): .. math:: - \eta_\mathrm{eff} = \min\left(\eta_0, \frac{\tau_y}{2\dot{\varepsilon}_{II}}\right) + \eta_\mathrm{eff} = \operatorname{softmin}_\delta\left(\eta_0, + \frac{\tau_y}{2\dot{\varepsilon}_{II}} + m\,\eta_\mathrm{reg}\right) and :math:`\tau_y` is the yield stress and :math:`\dot{\varepsilon}_{II}` - is the second invariant of the strain rate. + is the second invariant of the strain rate. The combination is the δ-smoothed + minimum by default (sqrt family, ``"onset"`` anchor, δ = 0.1: see + :attr:`yield_mode`); ``yield_mode="min"`` gives the exact ``min``. + :math:`\eta_\mathrm{reg}` is the declared ``plastic_rate_strengthening`` + (default 0, no term) and :math:`m` its homotopy multiplier + (``rate_strengthening_scale``, 1 = as stated). Parameters ---------- @@ -1520,12 +1649,20 @@ def __init__(self, unknowns, material_name: str = None): ) # Yield-combination mode (see _combine_yield on the base class). Default - # "min" = the exact hard Min(η_0, η_yield) this model has always used, so the - # default behaviour is unchanged. Opt into "softmin" (+ yield_smoother / - # yield_softness) for the δ-parameterised smooth-min homotopy. - self._yield_mode = "min" - self._yield_softness = 0.0 # δ; 0 ⇒ exact Min + # DEFAULT = the smooth yield law (since 2026-10-05): the exact hard Min is no + # longer the default. The kink broke the consistent tangent at the yield corner + # (and sympy's fuzzy comparisons on UWexpressions), and every hard case in the + # Spiegelman campaign had to set these by hand. The "onset" anchor keeps the + # viscous branch exact (a material that never yields is unchanged) and the + # plastic limit exact; delta = 0.1 under-stresses the yield point itself by + # 4.5% and leaves the yielded tangent at >= -1e-4 of the viscous one. (The + # "yield" anchor was tried as the default and refused in review: 5% stiffer + # unyielded material everywhere and a negative yielded tangent.) yield_mode = + # "min" with yield_softness = 0 is still available, as an explicit choice. + self._yield_mode = "softmin" + self._yield_softness = 0.1 # δ (sqrt family: width of the rounding in f) self._yield_smoother = "sqrt" # smooth-min family: "sqrt" | "powermean" + self._yield_anchor = "onset" # exact viscous branch and plastic limit self._yield_softness_expr = None # constants[] δ atom (created lazily) self._yield_offset_expr = None # onset-offset atom (created lazily) self._yield_sharpness_expr = None # power-mean sharpness atom (created lazily) @@ -1579,6 +1716,7 @@ class _Parameters(_ParameterBase, _ViscousParameterAlias): "Strain rate invariant minimum value", units="1/s", ) + plastic_rate_strengthening = _plastic_rate_strengthening_parameter() def __init__(inner_self, _owning_model): inner_self._owning_model = _owning_model @@ -1589,9 +1727,11 @@ def viscosity(self): r"""Effective viscosity with plastic yielding. .. math:: - \eta_{\mathrm{eff}} = \min\left(\eta_0, \frac{\tau_y}{2\dot{\varepsilon}_{II}}\right) + \eta_{\mathrm{eff}} = \operatorname{softmin}_\delta\left(\eta_0, + \frac{\tau_y}{2\dot{\varepsilon}_{II}} + m\,\eta_\mathrm{reg}\right) - where :math:`\dot{\varepsilon}_{II}` is the second invariant of strain rate. + where :math:`\dot{\varepsilon}_{II}` is the second invariant of strain rate + (see the class docstring for the combination and the regularisation). """ inner_self = self.Parameters # detect if values we need are defined or are placeholder symbols @@ -1632,11 +1772,16 @@ def viscosity(self): ) else: viscosity_yield = yield_stress / (2 * self._strainrate_inv_II) - - # Combine the viscous and plastic (yield) viscosities. The default - # yield_mode="min" gives the exact hard Min(η_0, η_yield); yield_mode="softmin" - # opts into the δ-parameterised smooth-min (sqrt or powermean family) for a - # scalable homotopy toward the sharp yield surface. + # The declared rate strengthening: eta_pl -> eta_pl + m * eta_reg. It floors the + # plastic viscosity AND the tangent modulus along the strain-rate direction in + # every yielded element (a static viscosity floor clips only the fastest cells and + # leaves the tangent at zero elsewhere — the mechanism mode of a yielded layer). + # Scale stated by the model; m is the homotopy multiplier, 1 = as stated. + viscosity_yield = self._with_rate_strengthening(viscosity_yield, inner_self) + + # Combine the viscous and plastic (yield) viscosities: the δ-parameterised + # smooth-min by default (see yield_mode / yield_smoother / yield_anchor); + # yield_mode="min" gives the exact hard Min(η_0, η_yield). effective_viscosity = self._combine_yield( inner_self.shear_viscosity_0, viscosity_yield ) @@ -1662,15 +1807,35 @@ def supports_yield_homotopy(self): can march it. See :meth:`_yield_homotopy_control`.""" return True + def _rate_strengthening_probe_terms(self): + # See ViscousFlowModel._rate_strengthening_probe_terms: the same floors as + # the viscosity above, so the ladder reads the law the solve runs. + inner_self = self.Parameters + tau = inner_self.yield_stress + if inner_self.yield_stress_min.sym != -sympy.oo: + tau = self._apply_floor(tau, inner_self.yield_stress_min) + edot = self._strainrate_inv_II + if not _parameter_is_zero(inner_self.strainrate_inv_II_min): + edot = edot + inner_self.strainrate_inv_II_min + eta_reg = inner_self.plastic_rate_strengthening + return dict( + viscous=inner_self.shear_viscosity_0, + eta_reg=eta_reg, + onset_num=2 * edot * inner_self.shear_viscosity_0 - tau, + onset_den=2 * edot * eta_reg, + ) + @property def yield_mode(self): r"""How the viscous and plastic (yield) viscosities are combined. - - ``"min"`` (default): exact hard ``Min(η_0, η_yield)`` — the sharp yield - surface this model has always used. + - ``"softmin"`` (default since 2026-10-05): the δ-parameterised smooth-min + (family set by ``yield_smoother``, pinned by ``yield_anchor``, rounded by + ``yield_softness``; default sqrt, "onset", δ = 0.1) — smooth enough for the + consistent tangent, and the scale of a homotopy toward the sharp surface. + - ``"min"``: exact hard ``Min(η_0, η_yield)`` — the sharp yield surface this + model used until 2026-10-05; its kink breaks the Newton tangent at yield. - ``"harmonic"``: ``1/(1/η_0 + 1/η_yield)`` — a smooth blend. - - ``"softmin"``: the δ-parameterised smooth-min (family set by - ``yield_smoother``), for a scalable homotopy toward the sharp surface. """ return self._yield_mode @@ -1918,6 +2083,7 @@ def __init__(self, unknowns, order=1, integrator: str = "bdf", self._check_order_supported(order) self._order = order self._yield_mode = "softmin" # "min", "harmonic", "smooth", or "softmin" + self._yield_anchor = "onset" # exact viscous branch and plastic limit self._yield_softness = 0.1 # δ parameter for "softmin" mode self._yield_smoother = "sqrt" # smooth-min family: "sqrt" | "powermean" self._yield_softness_expr = None # constants[] δ atom (created lazily) @@ -1940,7 +2106,7 @@ def __init__(self, unknowns, order=1, integrator: str = "bdf", # Persistent container for the yield-limited VEP effective viscosity # (the ViscoPlasticFlowModel._plastic_eff_viscosity pattern): the # combined coefficient is stored INSIDE this atom so the c-tensor - # bakes a wrapped atom and the default (Picard) tangent stays frozen; + # bakes a wrapped atom and the frozen (Picard, consistent_jacobian=False) tangent stays frozen; # only _jacobian_unwrap (Newton) sees the strain-rate dependence of # the yield law. Same freezing contract as the TI classes (#457/#493). self._vep_eff_viscosity = expression( @@ -2033,6 +2199,7 @@ def dt_elastic(inner_self, value): "Strain rate invariant minimum value", units="1/s", ) + plastic_rate_strengthening = _plastic_rate_strengthening_parameter() def __init__( inner_self, @@ -2632,7 +2799,7 @@ def viscosity(self): # Store the combined coefficient INSIDE the persistent container (the # ViscoPlasticFlowModel._plastic_eff_viscosity pattern) so the # c-tensor bakes ONE wrapped atom: sympy.diff treats it as a constant - # and the default (Picard) tangent stays frozen; only the Newton path + # and the frozen (Picard) tangent stays frozen; only the Newton path # (_jacobian_unwrap) sees the yield law's strain-rate dependence. # (Same freezing contract as the TI classes — issue #457 / PR #493.) self._vep_eff_viscosity._sym = effective_viscosity @@ -2687,6 +2854,7 @@ def _plastic_effective_viscosity(self): ) else: viscosity_yield = yield_stress / (2 * strainrate_inv_II) + viscosity_yield = self._with_rate_strengthening(viscosity_yield, parameters) return viscosity_yield @@ -2823,8 +2991,9 @@ def yield_mode(self): ``η_ve / g(f)`` where ``g(f) ≈ max(1, f)`` with smoothing parameter δ (``yield_softness``, default 0.1). Approaches exact Min as δ → 0; smooth derivatives at the kink. - Recommended default: gets within ~2 % of the true yield - surface while avoiding the SNES kink penalties of ``"min"``. + Recommended default: with the "onset" anchor it under-stresses the + yield point itself by 4.5 % and recovers the plastic limit to + 0.04 %, while avoiding the SNES kink penalties of ``"min"``. ``"harmonic"``: parallel blending — ``1/(1/η_ve + 1/η_pl)``. Smooth but undershoots τ_y for soft materials. ``"min"``: sharp cutoff — ``Min(η_ve, η_pl)``. @@ -3141,7 +3310,7 @@ def flux(inner_self, value: sympy.Matrix): Stores the raw sympy expression directly. Previously each component was wrapped in a UWexpression via ``validate_parameters`` — and the - default (Picard) tangent differentiates the flux WITHOUT unwrapping + frozen (Picard, consistent_jacobian=False) tangent differentiates the flux WITHOUT unwrapping (frozen-coefficient semantics; unwrap-before-differentiate runs only on the consistent-Newton path). An opaque wrapper holding the ENTIRE flux therefore has zero derivative w.r.t. the unknown and its @@ -3453,7 +3622,7 @@ def _build_c_tensor(self): # Bake the WRAPPED Parameter atoms into the tensor, exactly like the # isotropic ViscousFlowModel path. sympy.diff treats a UWexpression - # atom as a constant, so the default (Picard) tangent stays genuinely + # atom as a constant, so the frozen (Picard) tangent stays genuinely # frozen even when eta_0/eta_1 depend on strain rate; only the Newton # path (_jacobian_unwrap) exposes that dependence. Baking `.sym` # (unwrapped contents) here silently un-froze the TI Picard tangent — @@ -3748,6 +3917,7 @@ def dt_elastic(inner_self, value): lambda inner_self: 0, "Strain rate invariant minimum value", units="1/s", ) + plastic_rate_strengthening = _plastic_rate_strengthening_parameter() def __init__(inner_self, _owning_model): inner_self._owning_model = _owning_model @@ -4112,6 +4282,7 @@ def _plastic_effective_viscosity(self): ) else: viscosity_yield = tau_y / (2 * gamma_dot_abs) + viscosity_yield = self._with_rate_strengthening(viscosity_yield, parameters) return viscosity_yield diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 2e397319d..d128145d6 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -54,6 +54,70 @@ from underworld3.function._function import diff_wrt_field, derive_by_array_wrt_f from underworld3.function.expressions import unwrap_expression as _unwrap_expression +_DEFAULT_KRYLOV_RESTART = 100 + + +def _warn_from_caller(message, category=RuntimeWarning): + """Warn, attributed to the first frame outside the underworld3 package, so the + warning names the user's ``solve()`` call rather than an internal frame.""" + import os + import sys + import warnings + if sys.version_info >= (3, 12): + warnings.warn(message, category, + skip_file_prefixes=(os.path.dirname(uw.__file__),)) + else: + warnings.warn(message, category, stacklevel=3) + + +def _reaches(expr, targets): + """True when ``expr`` reaches any of ``targets``, following each UWexpression into + its contents. Every node is visited once. + + Unwrapping an expression and then searching it rebuilds the shared expression + graph as a tree, which is exponential in the nesting. Measured on the Spiegelman + notch (Drucker-Prager yield, residual strength, viscosity floor, smooth yield law): + the unwrap-and-differentiate linearity test ran for hours, and the unwrap-and-search + pressure test for minutes. This walk is linear in the size of the graph. + """ + from underworld3.function.expressions import UWexpression + targets = set(targets) + if not targets: + return False + seen = set() + stack = [expr] + while stack: + e = stack.pop() + key = id(e) + if key in seen: + continue + seen.add(key) + if isinstance(e, (sympy.MatrixBase, sympy.NDimArray)): + stack.extend(e) + continue + if not isinstance(e, sympy.Basic): + continue + if e in targets: + return True + if isinstance(e, UWexpression): + inner = getattr(e, "_sym", None) + if isinstance(inner, sympy.Basic) or isinstance(inner, (sympy.MatrixBase, sympy.NDimArray)): + stack.append(inner) + continue + stack.extend(e.args) + return False + + +def _check_homotopy_value(homotopy): + """``solve(homotopy=...)`` takes False, True or "rate_strengthening"; anything else + (a misspelt string is truthy) would silently run the δ march.""" + if homotopy is False or homotopy is None or homotopy is True or homotopy == "rate_strengthening": + return + raise ValueError( + f"homotopy must be False, True (the yield-softness march) or " + f"'rate_strengthening' (the declared regularisation ladder), got {homotopy!r}") + + def _public_names(cls): """The names ``uw.systems`` exports a solver class under.""" try: @@ -84,7 +148,7 @@ def _jacobian_unwrap(expr, route=None): No-op for constant-viscosity problems (eta has no grad-v dependence). The source is DIFFERENTIATION-SAFE: in every node body and at the top level, - each half-integer power whose base has free symbols gets ``+1e-36`` in its base. + each fractional power whose base has free symbols gets ``+1e-36`` in its base. Differentiating a bare invariant :math:`\dot\varepsilon_{II} = \sqrt{g}` produces :math:`\partial\sqrt{g}/\partial L = \dot\varepsilon/(2\dot\varepsilon_{II})` @@ -96,7 +160,10 @@ def _jacobian_unwrap(expr, route=None): kernel inherits it even at alpha = 0 because IEEE 0*NaN = NaN). The guard makes the derivative exactly zero at the singular point and perturbs it by under one part in 1e24 at any resolvable strain rate. - The RESIDUAL is never routed through here, and the default (Picard) + SymPy merges a power of the bare invariant into one that is not a + half-integer (``sqrt(g)**(1/3)`` is ``g**(1/6)``), so the guard covers + every non-integer exponent (#841). + The RESIDUAL is never routed through here, and the frozen (Picard) tangent never calls this function. See ``docs/developer/design/jit-shared-graph-codegen.md`` and @@ -108,21 +175,23 @@ def _jacobian_unwrap(expr, route=None): from underworld3.utilities import _jit_graph graph = _jit_graph.KernelGraph() - f = lambda e: _jit_graph.guard_half_integer_powers(graph.lower(e, guarded=True)) + f = lambda e: _jit_graph.guard_fractional_powers(graph.lower(e, guarded=True)) if isinstance(expr, sympy.MatrixBase): return expr.applyfunc(f) if isinstance(expr, sympy.NDimArray): return sympy.Array([f(e) for e in expr], expr.shape) return f(expr) # scalar expression - # the expanded route: development's body, unchanged + # the expanded route: the JIT before #823 tier 2, with the guard of #841 + from underworld3.utilities._jit_graph import is_fractional_power + eps2 = sympy.Float(1.0e-36) def _guard_sqrts(e): - # every HALF-INTEGER power: +1/2 (the invariant itself), -1/2 - # (its reciprocal in eta_pl = tau_y/(2 edot_II)), -3/2 (their - # derivatives), ... — all singular in value or derivative at a - # zero-argument state. + # every FRACTIONAL power (is_fractional_power): +1/2 (the invariant + # itself), -1/2 (its reciprocal in eta_pl = tau_y/(2 edot_II)), -3/2 + # (their derivatives), and the powers SymPy merges them into — all + # singular in value or derivative at a zero-argument state. # The same bottom-up rebuild as `e.replace(query, value)`, memoised on node # identity: the unwrapped flux repeats its shared sub-expressions as the same # object, and `replace` walked every occurrence (measured 16 s of a 108 s @@ -144,8 +213,7 @@ def _jacobian_unwrap(expr, route=None): if any(out == a and a != b for a, b in zip(args, new_args)): memo[id(n)] = (n, out) return out - if (out.is_Pow and out.exp.is_Rational and out.exp.q == 2 - and out.args[0].free_symbols): + if is_fractional_power(out): out = sympy.Pow(out.args[0] + eps2, out.exp) memo[id(n)] = (n, out) return out @@ -184,16 +252,69 @@ class SolverBaseClass(uw_object): # Jacobian tangent selection — validated property, see the # consistent_jacobian docstring below for the mode semantics. - self.consistent_jacobian = False + # DEFAULT = the consistent (Newton) tangent (since 2026-10-05): every Jacobian is + # assembled anyway; Picard is only the ONE automatic frozen-tangent step a cold + # start takes before Newton (see solve()). False (Picard throughout) and + # "continuation" (the residual-keyed blend) remain explicit choices. + self.consistent_jacobian = True # Picard->Newton continuation parameter (constants[]-routed so it can be # ramped at solve time without a JIT recompile). 0 = Picard, 1 = Newton. # Created LAZILY (see _get_newton_alpha) only when continuation is used, # so the default path constructs no extra UWexpression and therefore # cannot perturb global symbol-naming / JIT-cache state. self._newton_alpha = None - # Switch threshold: ramp alpha toward 1 once the relative residual falls - # below this (the basin-of-attraction heuristic for Newton). + # True once a Picard warm-up has been requested under the consistent tangent: + # the Jacobian is then compiled in the blended form J_p + alpha (J_n - J_p), + # alpha = 0 for the Picard steps, alpha = 1 for Newton. Sticky, so a solver + # that has warmed up once does not recompile back and forth (#791). + self._picard_blend = False + # Per-solve warm-up records (reset by _plan_picard_warmup); initialised here so + # a solver whose path never plans one (rotated free-slip) still has them. + self._picard_stages = None + self._picard_warmup_record = None + self._pending_picard = 0 + # Switch threshold of the rotated-free-slip two-phase path (rotated_bc.py): + # Picard until the relative residual falls below this, then Newton. self.newton_switch_rtol = 1.0e-2 + # consistent_jacobian = "continuation" (#791): alpha is keyed on the residual, + # alpha = clip(log10(F_ref / F) - log10(1 / start)) / (log10(1 / newton_at) + # - log10(1 / start)), never decreasing, + # with F_ref = ||F(u = 0)||, the problem's own residual scale (so a warm start + # enters the ramp where its residual puts it, not at alpha = 0). alpha leaves 0 + # at F = start * F_ref and reaches 1 (Newton) at F = newton_at * F_ref: the + # consistent tangent's linear solve fails from partly-yielded states above + # ~1e-3 F_ref on the Spiegelman notch. `l2` rather than `bt`: backtracking on + # ||F|| stops the frozen-tangent end, whose direction need not descend ||F||. + self.continuation_alpha_start = 1.0e-1 + self.continuation_newton_at = 5.0e-4 + self.continuation_linesearch = "l2" + # Back-off: when F climbs above `continuation_backoff` x the lowest F of the solve, + # alpha falls to the value F calls for (and the minimum resets). Without it, alpha + # cannot leave a high value once reached: measured on a power law n = 5, alpha + # reached 0.95, F then grew 40x and the solve ran out of iterations. None disables. + self.continuation_backoff = 10.0 + # Optional: switch the line search to `bt` (monotone ||F||) once alpha reaches this + # value; below it the frozen-tangent direction need not descend ||F||, so `bt` + # would stop it. None keeps `continuation_linesearch` throughout. (Experimental.) + self.continuation_bt_from = None + # When a continuation solve diverges after leaving alpha = 0 (Newton's linear + # solve failing, or a line-search failure), re-enter it warm up to + # `continuation_restarts` times, holding alpha = 0 (Picard) for the first + # `continuation_fallback_picard` iterations of each re-entry. + self.continuation_restarts = 3 + self.continuation_fallback_picard = 10 + # Whether the consistent (Newton) tangent differentiates the viscosity with respect + # to PRESSURE. With a pressure-dependent yield stress (Drucker-Prager, sigma_y = + # C + sin(phi) p) that derivative puts a d(2 eta eps)/dp term in the velocity- + # pressure block, so J_vp is no longer -B^T and the Schur complement is no longer + # negative definite: measured on the Spiegelman notch at eta_bg 5e24 the Newton + # J_vp was 41% away from -B^T, 40 of 1968 Schur eigenvalues changed sign, and the + # fieldsplit Schur preconditioner diverged at 0.99 per iteration with the velocity + # block solved EXACTLY, while the same operator with J_vp := -B^T converged. + # False keeps the Newton velocity block and takes the pressure column from the + # exact (frozen-viscosity) flux: a partial Newton tangent with the standard + # saddle-point structure. The residual is untouched either way. + self._newton_pressure_coupling = True # Fine-grained rebuild flags backing the is_setup property. See the # is_setup docstring and _build() for how these are consumed. @@ -302,6 +423,15 @@ class SolverBaseClass(uw_object): # the `tolerance` and `strategy` setters defeated an earlier attempt (#477). # Internal writers therefore go through _push_managed_option(). self._managed_pc_options = {} + # Krylov restart 100, not PETSc's 30, for every solver's outer KSP (Stokes adds + # its velocity split). The cost is memory (two vectors per iteration for FGMRES) + # and a general problem rarely runs past 30 anyway; the gain is that a near-null + # mode of the operator — a yielded layer's mechanism mode on the Spiegelman notch + # — stalls FGMRES(30) completely (rate 1.000) and is resolved by FGMRES(100) + # (188 its). Set once here, as a managed option, so a later `tolerance` or + # `strategy` assignment does not overwrite a user's own value. Recorded in + # solve_report.config["ksp"]["restart"]. + self._push_managed_option("ksp_gmres_restart", _DEFAULT_KRYLOV_RESTART) # Has _apply_preconditioner_options actually made the resolution decision # yet? Until it has, the options database holds only this solver's __init__ # defaults, which is NOT what the next solve will run — reporting them as @@ -452,25 +582,36 @@ class SolverBaseClass(uw_object): dispatch; the residual is never affected, so the converged solution always satisfies the exact constitutive law. - ``False`` (default) + ``False`` Differentiate the residual flux *as wrapped* — the effective - viscosity is frozen, giving a Picard / defect-correction tangent. - Globally robust; load-bearing for the tuned hard-yield viscoplastic - paths. - ``True`` + viscosity is frozen, giving a Picard / defect-correction tangent at + every iteration. The default until 2026-10-05; now an explicit choice. + ``True`` (default) Differentiate through the named quantities of the flux, so the tangent captures :math:`\partial\eta/\partial(\nabla v)` (full - Newton). Fast near the solution; its yield kink can stall the line - search far from it. + Newton). A COLD start takes one automatic Picard step first + (``solve(picard=...)`` adjusts it); a warm start takes none. With the + smooth yield law that is now the model default the yield kink is gone; + far from the solution a pressure-dependent yield still needs + ``newton_pressure_coupling`` / ``"continuation"`` (see the + ``nonlinear-solver`` guide). ``"continuation"`` Picard :math:`\rightarrow` Newton. Blend :math:`J(\alpha) = J_{\mathrm{picard}} + \alpha\,(J_{\mathrm{newton}} - J_{\mathrm{picard}})` with :math:`\alpha` a ``constants[]`` - parameter ramped :math:`0 \rightarrow 1` by a SNES monitor as the - residual drops. Picard locates the basin, Newton gives quadratic - convergence inside it (cf. Spiegelman et al. 2016; ASPECT - defect-correction-then-Newton). :math:`\alpha = 0` is bit-identical - to Picard, so no recompile is needed to switch. + parameter set at the start of every iteration from the residual: + :math:`\alpha` leaves 0 at :math:`F = s\,F_{\mathrm{ref}}` and reaches 1 + at :math:`F = n\,F_{\mathrm{ref}}`, linear in :math:`\log F`, never + decreasing, with :math:`F_{\mathrm{ref}} = \lVert F(u=0)\rVert`, + :math:`s` = ``continuation_alpha_start`` (0.1) and :math:`n` = + ``continuation_newton_at`` (5e-4). One SNES solve, ``l2`` line search + (``continuation_linesearch``), stopping at ``tolerance`` times + :math:`F_{\mathrm{ref}}`; a warm start enters the ramp at the + :math:`\alpha` its residual calls for. Picard locates the basin, Newton + converges inside it (cf. Spiegelman et al. 2016). :math:`\alpha = 0` is + bit-identical to Picard while the Newton blocks are finite (the guard of + :func:`_jacobian_unwrap` keeps them finite at rest), so no recompile is + needed to switch. The Newton flux for a model whose flux has a non-smooth yield kink is the model's own smooth law (``constitutive_model.flux_jacobian``) when @@ -492,13 +633,32 @@ class SolverBaseClass(uw_object): @consistent_jacobian.setter def consistent_jacobian(self, mode): if not mode: - self._consistent_jacobian = False - elif mode is True or mode == "continuation": - self._consistent_jacobian = mode - else: + mode = False + elif not (mode is True or mode == "continuation"): raise ValueError( f"consistent_jacobian must be False, True or 'continuation'; " f"got {mode!r}") + previous = getattr(self, "_consistent_jacobian", None) + self._consistent_jacobian = mode + if previous is not None and previous != mode: + # The tangent is compiled into the Jacobian kernel: a change after the + # first build must recompile it, or the next solve runs the old tangent + # while the solve report records the new one. + self._needs_function_rewire = True + + @property + def newton_pressure_coupling(self): + """Whether the consistent tangent differentiates the viscosity with respect to + the pressure (see the constructor comment). A change recompiles the Jacobian.""" + return getattr(self, "_newton_pressure_coupling", True) + + @newton_pressure_coupling.setter + def newton_pressure_coupling(self, value): + value = bool(value) + previous = getattr(self, "_newton_pressure_coupling", None) + self._newton_pressure_coupling = value + if previous is not None and previous != value: + self._needs_function_rewire = True @property def jit_route(self): @@ -545,18 +705,21 @@ class SolverBaseClass(uw_object): """Prepare a residual flux for Jacobian differentiation. ``expr`` is the exact (Picard) flux; ``newton_expr`` is the consistent - (Newton) flux to use — when None it is the unwrapped ``expr`` (the - derivative then captures d(eta)/d(grad v)). The residual itself is never + (Newton) flux to use — when None it is ``_jacobian_unwrap(expr)`` on the + solver's JIT route (the derivative then captures d(eta)/d(grad v)). The residual itself is never passed through here, so the converged solution always satisfies the exact constitutive law. Returns, by ``consistent_jacobian`` mode: * False -> ``expr`` (frozen viscosity; bit-identical Picard tangent) - * True -> ``newton_expr`` (full consistent Newton tangent) + * True -> ``newton_expr`` (full consistent Newton tangent); once a + cold start has planned its Picard step (``_picard_blend``), the + blend below, held at alpha = 1 except during that step * "continuation" -> ``expr + alpha*(newton_expr - expr)`` with alpha a constants[] parameter ramped 0->1 at solve time. Differentiating this gives ``G_picard + alpha*(G_newton - G_picard)`` because alpha - is constant w.r.t. the unknowns. alpha=0 is bit-identical to Picard. + is constant w.r.t. the unknowns. alpha=0 is bit-identical to Picard + while ``newton_expr``'s blocks are finite (0 * NaN is NaN). No-op for constant viscosity (``newton_expr`` == ``expr``). """ @@ -565,7 +728,7 @@ class SolverBaseClass(uw_object): return expr if newton_expr is None: newton_expr = _jacobian_unwrap(expr, route=self._jit_route_in_use()) - if mode == "continuation": + if mode == "continuation" or (mode is True and getattr(self, "_picard_blend", False)): a = self._get_newton_alpha() if isinstance(expr, sympy.MatrixBase): ne = newton_expr if isinstance(newton_expr, sympy.MatrixBase) \ @@ -592,7 +755,7 @@ class SolverBaseClass(uw_object): model whose flux has a non-smooth yield kink provide a physically motivated smooth tangent without changing the residual. """ - # Default (Picard) path never needs the Newton flux — short-circuit so + # The frozen (Picard) path never needs the Newton flux — short-circuit so # the (potentially expensive) model.flux_jacobian is not even evaluated, # keeping the default assembly allocation-free and bit-identical. if not self.consistent_jacobian: @@ -615,12 +778,14 @@ class SolverBaseClass(uw_object): def _get_newton_alpha(self): """Lazily construct the Picard->Newton continuation parameter. - Built on first use (continuation mode only) so the default Picard path - creates no extra UWexpression and leaves global symbol-naming / JIT-cache - state byte-identical to historical behaviour. + Built on first use (blended Jacobian only) so the plain Picard path + creates no extra UWexpression. One per solver: expressions are persistent + containers keyed by NAME, so a shared ``\alpha_{N}`` let one solver's + warm-up or ramp move another solver's tangent (measured: a second solver's + lazy creation reset the first one's alpha from 0.7 to 0.0). """ if self._newton_alpha is None: - self._newton_alpha = uw.function.expression( + self._newton_alpha = expression( r"\alpha_{N}", sympy.Float(0.0), "Picard->Newton continuation fraction (0=Picard, 1=Newton)", ) @@ -631,13 +796,15 @@ class SolverBaseClass(uw_object): alpha is routed through ``constants[]`` (it appears as a constant atom in the blended Jacobian), so this updates the live tangent WITHOUT a JIT - recompile. No-op outside ``consistent_jacobian == "continuation"`` (alpha - is then absent from the constants manifest). + recompile. Without a blended Jacobian (pure Newton never warmed up, or the + frozen tangent) alpha is absent from the constants manifest and has no effect. """ self._get_newton_alpha().sym = sympy.Float(value) try: self._update_constants() except Exception: + # Sanctioned: before the first build there is no DS to push constants to; + # the value is read when the constants are first packed. pass def set_custom_mg(self, coarse_meshes, kind="barycentric", verbose=False): @@ -1231,6 +1398,7 @@ class SolverBaseClass(uw_object): self._needs_dm_rebuild = False self._needs_bc_reregister = False self._needs_function_rewire = False + self._form_facts = None # forms (re)built: cached symbolic facts are stale else: self._needs_dm_rebuild = True self._needs_bc_reregister = True @@ -1264,6 +1432,153 @@ class SolverBaseClass(uw_object): """ return self._has_solution + def _plan_picard_warmup(self, picard, zero_init_guess, own_loop=False, staged_loop=False): + """Decide this solve's Picard warm-up (Layer 1). Every solver class calls this + from ``solve()`` before ``_build``. + + Arms ``self._pending_picard``, the number of frozen-tangent steps + :meth:`_snes_solve_with_retries` takes before Newton (inside the solve's + update hooks, wall-clock guard, iteration cap and report), and resets the + per-solve records. By :attr:`consistent_jacobian`: + + * ``True``: ``picard`` Picard steps, then Newton. A cold start takes ONE + automatically, except for a ``ksponly`` solve, a flux whose coefficient does + not depend on the unknowns, and a solver that runs its own loop + (``own_loop``: rotated free-slip, fault contact). ``picard < 0`` forbids it. + The first warm-up compiles the Jacobian in the blended form + ``J_p + alpha (J_n - J_p)`` (one rewire; sticky). + * ``"continuation"``: ``picard`` holds alpha = 0 for the first ``picard`` + iterations of the ramp. + * ``False``: every iteration already is a Picard step. + + ``staged_loop``: the solver's own loop stages Picard then Newton itself + (rotated free-slip); an explicit ``picard > 0`` under ``True`` then only + compiles the blended Jacobian for it to stage with. + """ + picard = int(picard or 0) + self._picard_stages = None + self._picard_warmup_record = None + self._pending_picard = 0 + self._warned_this_solve = set() + self._continuation_min_picard = 0 + self._continuation_stages = None # never report a previous solve's stages + mode = self.consistent_jacobian + if mode == "continuation": + self._continuation_min_picard = max(picard, 0) + return + if mode is not True: + return + if own_loop: + if staged_loop and picard > 0 and not self._picard_blend: + self._picard_blend = True + self._needs_function_rewire = True + return + # Cheapest test first: warm solves (every time step after the first) stop at + # zero_init_guess / the norm and never reach the symbolic linearity test. + if (picard == 0 and (zero_init_guess or self._solution_is_trivially_zero()) + and not self._snes_is_linear() and not self._flux_is_linear_in_unknowns()): + picard = 1 + if picard <= 0: + return + self._pending_picard = picard + if not self._picard_blend: + self._picard_blend = True + self._needs_function_rewire = True # recompile J in the blended form + + def _form_fact(self, key, compute): + """A symbolic fact about the solver's forms, cached until they can have changed. + + Recomputed whenever a rewire is pending or the solver is not set up (a + constitutive-parameter change sets ``_needs_function_rewire``; a new model or + new terms clear ``is_setup``), and dropped wherever ``_build`` consumes a + rewire or completes a setup — a warm solve's plan never reads the cache, so + without the drop the read after the build saw the facts from before the change + (measured: a viscous model switched to Drucker-Prager recorded + ``pressure_in_rheology=False`` and fired no warning). The facts cached here (flux linearity, pressure in the rheology) cost a + symbolic unwrap-and-differentiate — measured 3.7 s per call on a 16x8 VEP box — + which a time step must not pay on every warm solve. + """ + facts = getattr(self, "_form_facts", None) + if (facts is None or getattr(self, "_needs_function_rewire", False) + or not getattr(self, "is_setup", False)): + facts = self._form_facts = {} + if key not in facts: + facts[key] = compute() + return facts[key] + + def _flux_is_linear_in_unknowns(self): + """True when the residual flux's coefficient does not depend on the unknowns + (cached: see :meth:`_form_fact`). + + The frozen and consistent tangents then coincide, so the automatic cold-start + Picard step would be a second linear solve of the same operator for nothing + (and, through the blended kernel it compiles, a different code path — measured: + it diverged a linear constrained Stokes solve that Newton and Picard each solve + in one iteration). Test, on the WRAPPED flux: the coefficient is + unknown-dependent if any named sub-expression in it (the viscosity, the yield + stress) reaches an unknown — velocity, gradient or pressure symbol — searched + by :func:`_reaches` rather than unwrapped (see there for why), or if the flux + differentiated by a gradient symbol still carries an unknown at top level. A + spatially varying coefficient passes as linear, as it should. Scope: the flux F1 by the + gradient only — a nonlinear F0 (source) term, or F1 = grad u + b(u), reads as + linear here; that is what the guard needs (the warm-up changes the TANGENT of + the flux) but not a general linearity test. + """ + return self._form_fact("flux_linear", self._compute_flux_linearity) + + def _compute_flux_linearity(self): + try: + F1 = getattr(self, "F1", None) + if F1 is None: + return False + unknowns = set() + u = getattr(self, "u", None) + if u is not None: + unknowns.update(sympy.Matrix(u.sym)) + L = getattr(getattr(self, "Unknowns", None), "L", None) + if L is not None: + unknowns.update(sympy.Matrix(L)) + p = getattr(self, "p", None) + if p is not None: + unknowns.update(sympy.Matrix(p.sym)) + grads = list(sympy.Matrix(L)) if L is not None else [] + if not unknowns or not grads: + return False + # The WRAPPED flux: its named sub-expressions (the viscosity, the yield + # stress) are atoms here, so this is small. (1) Any named sub-expression that + # reaches an unknown makes the coefficient unknown-dependent; (2) the top + # level, differentiated by each gradient symbol, must carry no unknown. + from underworld3.function.expressions import UWexpression + flux = sympy.Array(F1.sym) + named = [a for a in flux.atoms(sympy.Symbol) if isinstance(a, UWexpression)] + if any(_reaches(a, unknowns) for a in named): + return False + for g in grads: + d = sympy.Array(diff_wrt_field(flux, g)) + if any(sympy.sympify(e).has(*unknowns) for e in sympy.flatten(d.tolist())): + return False + return True + except Exception: + # Sanctioned: a flux we cannot analyse is treated as nonlinear, so the cold + # start keeps the documented Layer-1 Picard step. + return False + + def _snes_is_linear(self): + """True when the configured SNES type is ``ksponly`` (a linear solve). + + A linear solve takes no automatic Picard step: the frozen and consistent + tangents coincide, and under Eisenstat-Walker the extra step only loosens the + first linear tolerance (measured on the spherical-shell Nitsche response: 12% + error with the step, 2% without — PR #716's finding, re-applied here). + """ + try: + opts = self.petsc_options + return opts.hasName("snes_type") and opts.getString("snes_type") == "ksponly" + except Exception: + # Sanctioned: an options database we cannot read means "not ksponly" — the + # nonlinear default, under which the automatic Picard step is the rule + return False + def _solution_is_trivially_zero(self): """True when the solution field is still identically zero. @@ -1276,13 +1591,18 @@ class SolverBaseClass(uw_object): the viscous one) but its derivative is not, and assembling the consistent tangent produces NaN. - Uses the PETSc vector norm, which is collective and therefore rank-uniform, - so every rank reaches the same decision. + Collective and rank-uniform: ``u.vec`` is the variable's LOCAL (sequential) + vector, so its norm differs between ranks, and the verdict decides a rewire and + an extra ``snes.solve`` — a rank-dependent answer hangs the job (measured at + np=2: one rank took the warm-up, the other did not). Each rank tests its own + entries and the answer is reduced over all of them. """ u = getattr(self, "u", None) if u is None: return False - return float(u.vec.norm()) == 0.0 + from mpi4py import MPI + local_nonzero = int(bool(np.any(np.asarray(u.vec.array) != 0.0))) + return not bool(uw.mpi.comm.allreduce(local_nonzero, op=MPI.MAX)) def _resolve_zero_init_guess(self, zero_init_guess): """Resolve the tri-state ``zero_init_guess`` argument of ``solve()``. @@ -1865,6 +2185,121 @@ class SolverBaseClass(uw_object): def solve_history(self, value): raise AttributeError("solve_history is read-only (set by solve()).") + def _snapshot_numerics(self): + """Record the numerical configuration in force for the NEXT ``snes.solve`` (#806). + + Read from the live PETSc objects, not from ``petsc_options``: a path that switches + a setting for the duration of one solve (the continuation ramp's line search and + atol, the difficulty probe's max_it) restores it before the report is captured, + so only a snapshot taken here, at the call, shows what actually ran. Stored for + :meth:`_capture_solve_report`. Read-only: nothing is set up or changed. + """ + snes = self.snes + cfg = {} + + def _safe(fn, default=None): + # Sanctioned swallow: a getter PETSc refuses on this object/state (e.g. a PC + # not yet set up) leaves the field at its honest "unavailable" value. + try: + return fn() + except Exception: + return default + + rtol, atol, stol, max_it = snes.getTolerances() + cfg["snes"] = dict(type=snes.getType(), rtol=rtol, atol=atol, stol=stol, + max_it=max_it, + max_funcs=_safe(snes.getMaxFunctionEvaluations), + ew=bool(_safe(snes.getUseEW, False))) + if cfg["snes"]["ew"]: + cfg["snes"]["ew_params"] = _safe(snes.getParamsEW) + # ⚠️ SNESGetNPC CREATES a nonlinear preconditioner when there is none, and the solve + # then runs it (measured: DIVERGED_INNER at iteration 0). Ask first. + if _safe(snes.hasNPC, False): + cfg["snes"]["npc"] = _safe(lambda: snes.getNPC().getType()) + # SNESGetLineSearch also creates one when absent; only line-search SNES types + # own one to report. + ls = (_safe(snes.getLineSearch) + if cfg["snes"]["type"] in ("newtonls", "nrichardson", "ncg", "qn", + "ngmres", "anderson", "ms") + else None) + if ls is not None: + cfg["linesearch"] = dict(type=_safe(ls.getType), order=_safe(ls.getOrder), + tolerances=_safe(ls.getTolerances)) + + def _ksp(k): + kr, ka, kd, km = k.getTolerances() + pc = k.getPC() + d = dict(prefix=_safe(k.getOptionsPrefix), type=_safe(k.getType), rtol=kr, + atol=ka, divtol=kd, max_it=km, pc=_safe(pc.getType)) + # The Krylov restart decides whether a near-null mode of the operator can be + # resolved at all (measured on the Spiegelman notch: a yielded layer gives the + # Newton velocity block a mechanism mode; FGMRES(30) stalls on it with NO + # progress, FGMRES(100) converges). petsc4py has no getter; read the option + # PETSc read, with PETSc's default. + if d["type"] in ("gmres", "fgmres", "lgmres", "dgmres", "pgmres", "pipefgmres"): + pfx = d["prefix"] or "" + d["restart"] = _safe(lambda: PETSc.Options().getInt(pfx + "ksp_gmres_restart", 30), 30) + return d, pc + + ksp = _safe(snes.getKSP) + if ksp is not None: + cfg["ksp"], pc = _ksp(ksp) + # The fieldsplit sub-solvers are read AFTER the solve (_capture_solve_report): + # before the first solve the PC is not set up, and PCFieldSplitGetSubKSP then + # errors (and prints, even when caught). Nothing switches them per solve, so the + # post-solve read is what ran. + + tangent = dict(consistent_jacobian=self.consistent_jacobian, + newton_pressure_coupling=getattr(self, "newton_pressure_coupling", True)) + # alpha only where it is live: the blended Jacobian (a warm-up compiled it, or + # the continuation ramp). Its value before the solve is the one the solve + # starts from; the ramp's end value is added by _capture_solve_report. + a = getattr(self, "_newton_alpha", None) + if a is not None and (self.consistent_jacobian == "continuation" + or (self.consistent_jacobian is True and self._picard_blend)): + tangent["alpha"] = _safe(lambda: float(a.sym)) + if self.consistent_jacobian == "continuation": + tangent["continuation_alpha_start"] = self.continuation_alpha_start + tangent["continuation_newton_at"] = self.continuation_newton_at + tangent["continuation_linesearch"] = self.continuation_linesearch + tangent["continuation_backoff"] = self.continuation_backoff + tangent["continuation_bt_from"] = self.continuation_bt_from + cfg["tangent"] = tangent + # Model-level facts a solver class knows about its own formulation (the pressure + # gauge for Stokes). Recorded WITH the numerics so a transcript can say what + # convention a solve ran under without introspecting a live process. + try: + self._snapshot_model_facts(cfg) + except Exception: + # Sanctioned: formulation facts are a record, never a reason to stop a + # solve; a class whose facts cannot be read records none. + pass + notes = cfg.pop("_warnings", ()) + self._numerics_at_solve = cfg + self._emit_solve_notes(notes) + return cfg + + def _emit_solve_notes(self, notes): + """Emit formulation warnings noted by :meth:`_snapshot_model_facts`, once per + solve (retries and re-entries snapshot again), attributed to the caller.""" + warned = getattr(self, "_warned_this_solve", None) + if warned is None: + warned = self._warned_this_solve = set() + for msg in notes: + if msg not in warned: + warned.add(msg) + _warn_from_caller(msg) + + def _snapshot_model_facts(self, cfg): + """Add formulation facts to the solve configuration snapshot (``solve_report.config``). + + Base: nothing. Stokes adds ``cfg["gauge"]``: whether a pressure nullspace was + requested and is attached, whether the rheology reads the pressure, pressure + Dirichlet BCs and natural (open) boundaries — the facts that decide whether the + pressure level a pressure-dependent yield law sees is physical or a gauge choice. + """ + return + def _capture_solve_report(self, *, bounded=False): """Record a SolveReport for the just-completed solve into ``self._solve_report`` and append it to ``self._solve_history``. Reads SNES state only (side-effect free). @@ -1902,6 +2337,15 @@ class SolverBaseClass(uw_object): hist = tuple(float(h) for h in self.snes.getConvergenceHistory()[0]) except Exception: hist = () + # A Picard warm-up ran as its own snes.solve just before this one: count it, and + # start the history (so fnorm0 and the reduction) at the cold residual. + warm = getattr(self, "_picard_warmup_record", None) + self._picard_warmup_record = None + if warm: + nl_its += warm["nl_its"] + ksp_its += warm["ksp_its"] + if warm["history"]: + hist = warm["history"] + hist[1:] fnorm0 = hist[0] if hist else None reduction = (fnorm / fnorm0) if (fnorm0 not in (None, 0.0)) else None rho = contraction(hist) @@ -1913,6 +2357,34 @@ class SolverBaseClass(uw_object): deadline_expired = (instrumentation is not None and instrumentation.deadline_expired) + config = getattr(self, "_numerics_at_solve", None) or {} + if warm and "tangent" in config: + config["tangent"]["picard_warmup_its"] = warm["nl_its"] + if "alpha" in config.get("tangent", {}): + try: + config["tangent"]["alpha_end"] = float(self._newton_alpha.sym) + except (AttributeError, TypeError): + # Sanctioned: no alpha container (never built); the start value stands + pass + _kcfg = config.get("ksp") + # Only once a linear solve has run: before that the PC is not set up, and + # PCFieldSplitGetSubKSP errors — printing to stderr even when caught (measured on + # a zero-iteration build solve). + if (_kcfg is not None and _kcfg.get("pc") == "fieldsplit" and "splits" not in _kcfg + and ksp_its > 0): + try: + _splits = [] + for _sk in self.snes.getKSP().getPC().getFieldSplitSubKSP(): + _r, _a, _d, _m = _sk.getTolerances() + _splits.append(dict(prefix=_sk.getOptionsPrefix(), type=_sk.getType(), + rtol=_r, atol=_a, divtol=_d, max_it=_m, + pc=_sk.getPC().getType())) + _kcfg["splits"] = _splits + except Exception: + # Sanctioned swallow: a solve that failed before PC set-up has no + # sub-solvers to report; the field stays absent. + pass + report = SolveReport( reason=reason, reason_str=reason_string(reason), @@ -1928,7 +2400,9 @@ class SolverBaseClass(uw_object): bounded=bool(bounded), sub=sub, deadline_expired=deadline_expired, + config=config, ) + self._numerics_at_solve = None # never report a previous solve's config self._solve_report = report self._solve_history.append(report) self._record_solve_outcome(report) @@ -1967,11 +2441,23 @@ class SolverBaseClass(uw_object): rnorm0 = info.get("rnorm0") fnorm0 = float(rnorm0) if rnorm0 else None reduction = (fnorm / fnorm0) if (fnorm0 and fnorm == fnorm) else None + # The rotated loop runs its own Krylov solves, so there is no SNES snapshot to + # read; record the formulation facts the transcript needs (tangent, gauge, + # regularisation) from the solver itself. + config = dict(tangent=dict(consistent_jacobian=self.consistent_jacobian, + newton_pressure_coupling=self.newton_pressure_coupling, + path="rotated_freeslip")) + try: + self._snapshot_model_facts(config) + except Exception: + # Sanctioned: formulation facts are a record; none is recorded on failure + pass + self._emit_solve_notes(config.pop("_warnings", ())) report = SolveReport( reason=reason, reason_str=ksp_reason_string(reason), converged=converged, nl_its=nl_its, ksp_its=ksp_its, fnorm=fnorm, fnorm0=fnorm0, reduction=reduction, rho=None, fev=None, history=(), - bounded=False, + bounded=False, config=config, ) self._solve_report = report self._solve_history.append(report) @@ -2088,8 +2574,8 @@ class SolverBaseClass(uw_object): ----- Not available with rotated free-slip BCs: that path solves outside ``self.snes``, so the cap and anchor cannot reach it. Under - ``consistent_jacobian="continuation"`` the cap applies to EACH stage (Picard - then Newton), so one probe may run up to twice ``max_nl_its``. + ``consistent_jacobian="continuation"`` the cap applies to the one ramped + solve (Picard and Newton iterations together). """ if getattr(self, "_rotated_freeslip_bcs", None) \ or getattr(self, "_fault_contact_faults", None): @@ -2344,9 +2830,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() @@ -2354,10 +2839,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: @@ -2431,18 +2912,26 @@ class SolverBaseClass(uw_object): # start/stop/restart terminates at the same point an uninterrupted solve would. # Tolerances are saved and restored so nothing leaks into later solves. _probe = bool(self._difficulty_probe) + _warm = int(getattr(self, "_pending_picard", 0) or 0) + self._pending_picard = 0 _saved_tol = None - if _probe: + if _probe or _warm: _saved_tol = self.snes.getTolerances() + if _probe: if self._difficulty_max_it is not None: self.snes.setTolerances(max_it=int(self._difficulty_max_it)) if self._resume_abs_target is not None: self.snes.setTolerances(atol=float(self._resume_abs_target)) try: + if _warm: + # The planned Picard warm-up (_plan_picard_warmup), inside this solve's + # hooks, guard, cap and report rather than before them. + self._picard_warmup(gvec, _warm, count_against_cap=_probe, verbose=verbose) if self.consistent_jacobian == "continuation": self._continuation_solve(gvec, verbose=verbose) else: + self._snapshot_numerics() self.snes.solve(None, gvec) if divergence_retries > 0: for _r in range(divergence_retries): @@ -2455,7 +2944,13 @@ class SolverBaseClass(uw_object): f"warm-start retry {_r + 1}/{divergence_retries}", flush=True, ) - self.snes.solve(None, gvec) + if self.consistent_jacobian == "continuation": + # a warm re-entry: the residual-keyed ramp resumes at the alpha + # the current residual calls for + self._continuation_solve(gvec, verbose=verbose) + else: + self._snapshot_numerics() + self.snes.solve(None, gvec) finally: if _saved_tol is not None: self.snes.setTolerances(rtol=_saved_tol[0], atol=_saved_tol[1], @@ -2467,32 +2962,204 @@ class SolverBaseClass(uw_object): # produced under an estimate_difficulty cap. self._capture_solve_report(bounded=_probe) - def _continuation_solve(self, gvec, verbose=False): - """Picard -> Newton continuation via the constants[]-routed alpha. - - Stage 1 solves with the frozen (Picard) tangent (alpha=0) to a loose - tolerance to enter Newton's basin of attraction; stage 2 ramps to the - consistent (Newton) tangent (alpha=1) and warm-starts to the requested - tolerance. alpha is toggled through ``constants[]`` so neither stage - triggers a JIT recompile (cf. Spiegelman et al. 2016; ASPECT). + def _picard_warmup(self, gvec, n, count_against_cap=False, verbose=False): + """Take ``n`` Picard steps — Newton iterations with the FROZEN tangent — in place, + and set the tolerances of the Newton solve that follows. + + Called by :meth:`_snes_solve_with_retries` (which saves and restores the SNES + tolerances around the whole solve). Requires the blended Jacobian + (``_picard_blend``): alpha = 0 is the frozen (Picard) operator, i.e. a linear + solve with the coefficient held at the current state; alpha is a + ``constants[]`` atom, so switching costs no recompile. Runs exactly ``n`` steps + (no relative stop; an absolute stop at the solve's ``atol`` is honoured), then + restores alpha = 1. The Newton solve is then anchored to the COLD residual: + ``atol >= rtol * ||F(u_cold)||``, so the solve stops where one uninterrupted + solve would. Without the anchor PETSc measures ``rtol`` from the post-warm-up + residual and asks Newton for a further ``rtol`` below an already-converged + state (measured: SNORM stagnation on a near-linear viscoplastic box that + converges in 2 iterations without the warm-up). With ``count_against_cap`` + (an ``estimate_difficulty`` probe) the steps come out of the probe's + iteration cap. Recorded in ``_picard_stages`` and, for the solve report, + ``_picard_warmup_record``. """ rtol, atol, stol, max_it = self.snes.getTolerances() + n = int(n) + if count_against_cap and max_it: + n = min(n, int(max_it)) + its = 0 + hist = () + try: + self._set_newton_alpha(0.0) + self.snes.setTolerances(rtol=0.0, atol=atol, stol=stol, max_it=n) + self.snes.solve(None, gvec) + its = int(self.snes.getIterationNumber()) + reason = int(self.snes.getConvergedReason()) + try: + hist = tuple(float(h) for h in self.snes.getConvergenceHistory()[0]) + except Exception: + # Sanctioned: no history recorded; the report keeps the Newton one + hist = () + # Kept for the solve report, which reads only the SNES state of the Newton + # solve that follows: without this the warm-up's iterations, linear solves + # and initial residual vanish from it (measured: a two-step solve reported + # nl_its = 1 and a reduction relative to the post-warm-up residual). + self._picard_warmup_record = dict( + nl_its=its, ksp_its=int(self.snes.getLinearSolveIterations()), history=hist) + if reason < 0 and reason != int(PETSc.SNES.ConvergedReason.DIVERGED_MAX_IT): + from underworld3.systems.solve_report import reason_string + inst = getattr(self, "_instrumentation", None) + if inst is not None and getattr(inst, "deadline_expired", False): + _warn_from_caller( + f"[{self.name}] the wall-clock deadline expired during the Picard " + f"warm-up ({reason_string(reason)}); the solve stops there.") + else: + _warn_from_caller( + f"[{self.name}] the Picard warm-up step ended {reason_string(reason)} " + f"({reason}); Newton continues from its iterate.") + if verbose: + uw.pprint(f"Picard warm-up: {its} frozen-tangent step(s), reason={reason}") + finally: + self._set_newton_alpha(1.0) + new_atol = atol + if hist and rtol > 0.0: + new_atol = max(atol, rtol * hist[0]) + new_max_it = max(int(max_it) - its, 0) if (count_against_cap and max_it) else max_it + self.snes.setTolerances(rtol=rtol, atol=new_atol, stol=stol, max_it=new_max_it) + self._picard_stages = dict(picard_iterations=its) + + def _continuation_alpha(self, F, F_ref): + """alpha for residual norm ``F`` on the continuation ramp (before monotonicity).""" + if F <= 0.0: + return 1.0 + L = np.log10(F_ref / F) + L0 = np.log10(1.0 / self.continuation_alpha_start) + L1 = np.log10(1.0 / self.continuation_newton_at) + if L1 <= L0: + # start == newton_at: a SWITCH, Picard above the gate and Newton below it, + # no blended tangent in between + return 1.0 if L >= L1 else 0.0 + return float(min(1.0, max(0.0, (L - L0) / (L1 - L0)))) - # 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) - if verbose and uw.mpi.rank == 0: - print(f"continuation Picard: reason={self.snes.getConvergedReason()} " - f"it={self.snes.getIterationNumber()}", flush=True) + def _continuation_solve(self, gvec, verbose=False): + """Picard -> Newton continuation in ONE SNES solve, alpha keyed on the residual. + + The Jacobian is the blended form J_p + alpha (J_n - J_p) (alpha a + ``constants[]`` atom, so no recompile). At the start of every iteration alpha + is set from the current residual F relative to F_ref = ||F(u = 0)|| + (:meth:`_continuation_alpha`, never decreasing within a solve); an explicit + ``solve(picard=N)`` holds alpha = 0 for the first N iterations. The solve + stops on the ABSOLUTE target ``tolerance * F_ref`` (or tighter, if already + set), the point a cold solve stops at, so a warm start is not asked for a + further ``tolerance`` reduction of an already small residual. The line search + is ``continuation_linesearch`` for this solve (None leaves it alone). + + Measured (Spiegelman notch, power-law and DP cases; #791): converges cases + where both pure Newton and the old two-stage continuation fail, and makes + warm starts cheaper than cold ones. Known limitation: at intermediate alpha + the iteration can oscillate and stall (notch at refinement 2; power law n=3 + with yield). + """ + rtol, atol, stol, max_it = self.snes.getTolerances() + ls = self.snes.getLineSearch() + ls_type = ls.getType() + + # F_ref: the residual at u = 0 — one residual evaluation, whatever the start. + _x = gvec.duplicate() + _x.zeroEntries() + _f = gvec.duplicate() + self.snes.computeFunction(_x, _f) + F_ref = _f.norm() + _x.destroy() + _f.destroy() + _fs = gvec.duplicate() + self.snes.computeFunction(gvec, _fs) + F_start = _fs.norm() + _fs.destroy() + # u = 0 can be (near) the solution, e.g. a problem driven only by a weak body + # force from a good guess: never key on a scale below the starting residual. + F_ref = max(F_ref, F_start) + + hold = dict(n=int(getattr(self, "_continuation_min_picard", 0) or 0)) + backoff = self.continuation_backoff + state = dict(alpha=0.0, frozen=0, blended=0, newton=0, backoffs=0, F_min=None, + restarts=0, nl=0, ksp=0) + user_hook = bool(self._snes_update_callbacks) + + def _ramp(snes, k): + F = snes.getFunctionNorm() + a_key = self._continuation_alpha(F, F_ref) + if state["F_min"] is None or F < state["F_min"]: + state["F_min"] = F + if k < hold["n"]: + a = 0.0 + elif backoff and F > backoff * state["F_min"] and a_key < state["alpha"]: + # the iteration is moving away from the solution: follow F back down + a = a_key + state["F_min"] = F + state["backoffs"] += 1 + else: + a = max(state["alpha"], a_key) + state["alpha"] = a + self._set_newton_alpha(a) + bt_from = self.continuation_bt_from + if bt_from is not None: + want = "bt" if a >= bt_from else (self.continuation_linesearch or ls_type) + if ls.getType() != want: + ls.setType(want) + if a <= 0.0: + state["frozen"] += 1 + elif a >= 1.0: + state["newton"] += 1 + else: + state["blended"] += 1 + if user_hook: + self._dispatch_snes_update(snes, k) - # Stage 2 — Newton (alpha=1), requested tolerance, warm-started. - 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) - if verbose and uw.mpi.rank == 0: - print(f"continuation Newton: reason={self.snes.getConvergedReason()} " - f"it={self.snes.getIterationNumber()}", flush=True) + try: + if self.continuation_linesearch: + ls.setType(str(self.continuation_linesearch)) + self.snes.setTolerances(rtol=rtol, atol=max(atol, self.tolerance * F_ref), + stol=stol, max_it=max_it) + self.snes.setUpdate(_ramp) + self._snapshot_numerics() + self.snes.solve(None, gvec) + state["nl"] += int(self.snes.getIterationNumber()) + state["ksp"] += int(self.snes.getLinearSolveIterations()) + # Fall back to Picard when the solve diverged after leaving alpha = 0. + # Only failures a Picard block can repair: Newton's linear solve failing, the + # line search finding no step, a non-finite residual. NOT DIVERGED_MAX_IT — + # restarting an exhausted budget just spends it again (measured: 4 x 800). + _R = PETSc.SNES.ConvergedReason + _recoverable = (_R.DIVERGED_LINEAR_SOLVE, _R.DIVERGED_LINE_SEARCH, + _R.DIVERGED_FUNCTION_NANORINF, _R.DIVERGED_INNER) + while (self.snes.getConvergedReason() in _recoverable and state["alpha"] > 0.0 + and state["restarts"] < int(self.continuation_restarts or 0)): + state["restarts"] += 1 + if verbose: + uw.pprint(f"continuation: {self.snes.getConvergedReason()} at alpha=" + f"{state['alpha']:.3f}; Picard fallback " + f"{state['restarts']}/{self.continuation_restarts}") + hold["n"] = int(self.continuation_fallback_picard or 0) + state["alpha"] = 0.0 + state["F_min"] = None + self.snes.solve(None, gvec) + state["nl"] += int(self.snes.getIterationNumber()) + state["ksp"] += int(self.snes.getLinearSolveIterations()) + finally: + self.snes.setUpdate(None) + self._attach_snes_update_hook() + ls.setType(ls_type) + self.snes.setTolerances(rtol=rtol, atol=atol, stol=stol, max_it=max_it) + + # Iterations by tangent (the hook fires once at the start of each iteration). + self._continuation_stages = dict( + frozen_iterations=state["frozen"], blended_iterations=state["blended"], + newton_iterations=state["newton"], final_alpha=state["alpha"], + backoffs=state["backoffs"], restarts=state["restarts"], + total_nl_iterations=state["nl"], total_ksp_iterations=state["ksp"], F_ref=F_ref) + if verbose: + uw.pprint(f"continuation: reason={self.snes.getConvergedReason()} " + f"it={self.snes.getIterationNumber()} {self._continuation_stages}") # Restore a clean Picard tangent for any subsequent solve (next step). self._set_newton_alpha(0.0) @@ -2538,6 +3205,7 @@ class SolverBaseClass(uw_object): self._current_jit_cache_key == self._last_jit_cache_key: self._update_constants() self.is_setup = True + self._form_facts = None # forms (re)built: cached symbolic facts are stale if hasattr(self, "constitutive_model") and \ self.constitutive_model is not None and \ hasattr(self.constitutive_model, "_solver_is_setup"): @@ -2558,6 +3226,7 @@ class SolverBaseClass(uw_object): self._setup_solver(verbose, _rewire_only=True) self.is_setup = True + self._form_facts = None # forms (re)built: cached symbolic facts are stale # _last_jit_cache_key must track "which bundle the solver is # currently wired for", not "what the last full-build was". @@ -2661,6 +3330,7 @@ class SolverBaseClass(uw_object): self._setup_solver(verbose) self.is_setup = True + self._form_facts = None # forms (re)built: cached symbolic facts are stale # Record cache key after full build — used by the fast path # on subsequent _build() calls to detect constants-only changes. @@ -4500,6 +5170,7 @@ class SNES_Scalar(SolverBaseClass): UW_DMPlexSetSNESLocalFEM(cdm.dm, PETSC_FALSE, NULL) self.is_setup = True + self._form_facts = None # forms (re)built: cached symbolic facts are stale self.constitutive_model._solver_is_setup = True @timing.routine_timer_decorator @@ -4510,7 +5181,8 @@ class SNES_Scalar(SolverBaseClass): debug: bool=False, debug_name: str=None, time=None, - divergence_retries: int=0, ): + divergence_retries: int=0, + picard: int = 0, ): """ Solve the system of equations. @@ -4547,6 +5219,11 @@ class SNES_Scalar(SolverBaseClass): If SNES reports DIVERGED after the solve, re-call it with warm start up to this many times. A single retry rescues most VEP yield-surface kink divergences. 0 preserves legacy behaviour. + picard : int, default=0 + Picard (frozen-tangent) steps before Newton under the default + ``consistent_jacobian=True``. A cold start of a nonlinear flux takes + ONE automatically; ``picard=N`` asks for N, ``picard=-1`` for none. + Keyword only in practice: it follows every other argument. Returns ------- @@ -4592,7 +5269,13 @@ class SNES_Scalar(SolverBaseClass): # and resolving earlier would warm-start off the flag it just cleared. zero_init_guess = self._resolve_zero_init_guess(zero_init_guess) + # Layer 1: the Picard warm-up for this solve, decided before _build (a warm-up + # under the consistent tangent needs the blended Jacobian compiled in). + self._plan_picard_warmup(picard, zero_init_guess, + own_loop=bool(getattr(self, "_rotated_freeslip_bcs", None))) self._build(verbose, debug, debug_name) + if self._picard_blend and self.consistent_jacobian is True: + self._set_newton_alpha(1.0) # Newton unless the warm-up says otherwise # Set time on the DM so petsc_t is available in pointwise functions cdef DM _time_dm @@ -5557,6 +6240,7 @@ class SNES_Vector(SolverBaseClass): self.is_setup = True + self._form_facts = None # forms (re)built: cached symbolic facts are stale self.constitutive_model._solver_is_setup = True @@ -5568,6 +6252,7 @@ class SNES_Vector(SolverBaseClass): debug=False, debug_name=None, divergence_retries: int=0, + picard: int = 0, ): """ Solve the vector field system of equations. @@ -5597,6 +6282,11 @@ class SNES_Vector(SolverBaseClass): divergence_retries : int, default=0 If SNES reports DIVERGED after the solve, re-call it with warm start up to this many times. 0 preserves legacy behaviour. + picard : int, default=0 + Picard (frozen-tangent) steps before Newton under the default + ``consistent_jacobian=True``. A cold start of a nonlinear flux takes + ONE automatically; ``picard=N`` asks for N, ``picard=-1`` for none. + Keyword only in practice: it follows every other argument. Returns ------- @@ -5625,7 +6315,13 @@ class SNES_Vector(SolverBaseClass): # and resolving earlier would warm-start off the flag it just cleared. zero_init_guess = self._resolve_zero_init_guess(zero_init_guess) + # Layer 1: the Picard warm-up for this solve, decided before _build (a warm-up + # under the consistent tangent needs the blended Jacobian compiled in). + self._plan_picard_warmup(picard, zero_init_guess, + own_loop=bool(getattr(self, "_rotated_freeslip_bcs", None))) self._build(verbose, debug, debug_name) + if self._picard_blend and self.consistent_jacobian is True: + self._set_newton_alpha(1.0) # Newton unless the warm-up says otherwise gvec = self.dm.getGlobalVec() @@ -6273,6 +6969,7 @@ class SNES_MultiComponent(SolverBaseClass): UW_DMPlexSetSNESLocalFEM(cdm.dm, PETSC_FALSE, NULL) self.is_setup = True + self._form_facts = None # forms (re)built: cached symbolic facts are stale self.constitutive_model._solver_is_setup = True @timing.routine_timer_decorator @@ -6283,6 +6980,7 @@ class SNES_MultiComponent(SolverBaseClass): debug=False, debug_name=None, divergence_retries: int=0, + picard: int = 0, ): """Solve the multi-component SNES problem. @@ -6294,6 +6992,11 @@ class SNES_MultiComponent(SolverBaseClass): divergence_retries : int, default=0 If SNES reports DIVERGED after the solve, re-call it with warm start up to this many times. 0 preserves legacy behaviour. + picard : int, default=0 + Picard (frozen-tangent) steps before Newton under the default + ``consistent_jacobian=True``. A cold start of a nonlinear flux takes + ONE automatically; ``picard=N`` asks for N, ``picard=-1`` for none. + Keyword only in practice: it follows every other argument. """ @@ -6309,7 +7012,13 @@ class SNES_MultiComponent(SolverBaseClass): # and resolving earlier would warm-start off the flag it just cleared. zero_init_guess = self._resolve_zero_init_guess(zero_init_guess) + # Layer 1: the Picard warm-up for this solve, decided before _build (a warm-up + # under the consistent tangent needs the blended Jacobian compiled in). + self._plan_picard_warmup(picard, zero_init_guess, + own_loop=bool(getattr(self, "_rotated_freeslip_bcs", None))) self._build(verbose, debug, debug_name) + if self._picard_blend and self.consistent_jacobian is True: + self._set_newton_alpha(1.0) # Newton unless the warm-up says otherwise gvec = self.dm.getGlobalVec() @@ -6585,6 +7294,9 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self._tolerance = 1.0e-4 self._strategy = "default" + # the velocity split's Krylov restart, beside the outer one set in the base + # constructor (see there): the yielded layer's mechanism mode lives in this block + self._push_managed_option("fieldsplit_velocity_ksp_gmres_restart", _DEFAULT_KRYLOV_RESTART) # Owned-option latch state (see _resolve_owned_option): the value THIS # solver last pushed per key, and any user value latched per key. @@ -7910,6 +8622,115 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): # changes (e.g. mesh deformation) — see _build_velocity_rotation_nullspace. self._velocity_rotation_nullspace = None + def _pressure_in_rheology(self): + """True when the constitutive model's viscosity depends on the solver's pressure. + + Drucker-Prager (``tau_y = C + sin(phi) p``) is the usual case. It matters because + the pressure LEVEL is then physics, not a gauge: a mean-pressure nullspace projection + changes the yield stress, and the Newton tangent gains a d(stress)/dp column. + """ + return self._form_fact("pressure_in_rheology", + lambda: self._reads_pressure(lambda cm: cm.viscosity)) + + def _yield_reads_pressure(self): + """True when the YIELD STRESS depends on the pressure (Drucker-Prager), as opposed + to a pressure-dependent viscous branch. ``None`` when it cannot be determined.""" + return self._form_fact( + "yield_reads_pressure", + lambda: self._reads_pressure( + lambda cm: getattr(getattr(cm, "Parameters", None), "yield_stress", None))) + + def _reads_pressure(self, pick): + cm = getattr(self, "constitutive_model", None) + psym = getattr(getattr(self, "p", None), "sym", None) + if cm is None or psym is None: + return False + try: + expr = pick(cm) + if expr is None: + return False + return _reaches(expr, {psym[0]}) + except Exception: + # Sanctioned: an expression we cannot unwrap is recorded as unknown (None), + # and neither pressure warning fires on an unknown. + return None + + def _snapshot_model_facts(self, cfg): + """``cfg["gauge"]``: the pressure-level convention this Stokes solve runs under.""" + attached = False + try: + # Before the first solve SNESGetJacobian hands back NULL matrices; a method + # call on one is a segfault, not an exception. Check the handle first. + if self.snes is not None: + jac = self.snes.getJacobian() + J = jac[0] if jac else None + if J is not None and J.handle != 0: + ns = J.getNullSpace() + attached = bool(ns is not None and ns.handle != 0) + except Exception: + # Sanctioned: a Jacobian PETSc will not describe yet records "not attached" + attached = False + try: + pbcs = sorted({bc.boundary for bc in self._pressure_dirichlet_bcs()}) + except Exception: + # Sanctioned: no pressure field decomposition yet; recorded as none + pbcs = [] + cm = getattr(self, "constitutive_model", None) + reg = {} + for name in ("plastic_rate_strengthening", "shear_viscosity_min", "yield_stress_min", "strainrate_inv_II_min"): + par = getattr(getattr(cm, "Parameters", None), name, None) + if par is not None: + try: + reg[name] = str(par.sym) + except Exception: + # Sanctioned: a value with no printable form is recorded as None + reg[name] = None + for name in ("yield_mode", "yield_smoother", "yield_anchor", "yield_softness", "rate_strengthening_scale"): + if hasattr(cm, name): + reg[name] = getattr(cm, name) + cfg["regularisation"] = reg + in_rheology = self._pressure_in_rheology() + # Warnings are NOTED here and emitted by _snapshot_numerics after the snapshot is + # stored (so a warnings-as-errors caller cannot leave the config half-written), + # once per solve. + notes = [] + if self.consistent_jacobian is not False and self._yield_reads_pressure() is True: + eta_reg = getattr(getattr(cm, "Parameters", None), "plastic_rate_strengthening", None) + from underworld3.constitutive_models import _parameter_is_zero + if eta_reg is not None and _parameter_is_zero(eta_reg): + notes.append( + f"[{self.name}] a pressure-dependent yield law with no stated " + f"plastic_rate_strengthening under the consistent tangent: the yielded " + f"tangent has no modulus along the strain-rate direction, which is the " + f"near-singular velocity block measured on the Spiegelman notch. State " + f"the regularisation (Parameters.plastic_rate_strengthening, units of " + f"viscosity) — its scale is the model's, there is no universal default " + f"— or accept the risk. Recorded in solve_report.config['regularisation'].") + if self._petsc_use_pressure_nullspace and in_rheology is True: + # The projection fixes mean(p) = 0 on every Krylov solve. With a pressure- + # dependent viscosity (Drucker-Prager) the yield stress then follows that + # gauge choice, and moves with it as the yielded region changes between + # iterations. Legitimate only when the physical level is supplied separately + # (lithostatic pressure written into the yield law, p the dynamic part). + notes.append( + f"[{self.name}] petsc_use_pressure_nullspace=True with a viscosity that " + f"depends on the pressure: the nullspace projection imposes mean(p) = 0, so " + f"the yield stress / viscosity reads a GAUGE-fixed pressure, not a physical " + f"one. Either the yield law must add the reference (lithostatic) pressure " + f"explicitly, or the level must come from the boundary conditions (a free " + f"surface or traction boundary, in which case turn the nullspace flag off; " + f"a pressure Dirichlet BC for a closed box). Recorded in " + f"solve_report.config['gauge'].") + if notes: + cfg["_warnings"] = notes + cfg["gauge"] = dict( + pressure_nullspace_requested=bool(self._petsc_use_pressure_nullspace), + pressure_nullspace_attached=attached, + pressure_in_rheology=in_rheology, + pressure_dirichlet_bcs=pbcs, + natural_bcs=sorted({bc.boundary for bc in (getattr(self, "natural_bcs", None) or [])}), + ) + def _pressure_dirichlet_bcs(self): """Return essential boundary conditions applied to the pressure field.""" @@ -8330,7 +9151,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): # a no-op (eta has no grad-v dependence). See # docs/developer/design/jacobian-consistent-tangent.md # - # (see consistent_jacobian / _jacobian_source: default Picard, bit- + # (see consistent_jacobian / _jacobian_source: False = Picard, bit- # identical; True -> Newton; "continuation" -> alpha-blended.) F0_jac = self._jacobian_source(F0) PF0_jac = self._jacobian_source(PF0) @@ -8437,8 +9258,13 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): G1[fc, dg] = dF0_dGp[fc] # up_G2[fc, df] = dF1[fc, df] / dp + # newton_pressure_coupling=False differentiates the exact flux here: its viscosity + # is an opaque UWexpression atom, so d/dp sees only the -p I term (= Picard's + # pressure column, -B^T) while the uu block above keeps the Newton tangent. + F1_for_up = F1_for_jac if self.newton_pressure_coupling \ + else sympy.Array(F1).reshape(dim, dim) G2 = sympy.zeros(dim, dim) - dF1_dp = diff_wrt_field(F1_for_jac, p_scalar) + dF1_dp = diff_wrt_field(F1_for_up, p_scalar) for fc in range(dim): for df in range(dim): G2[fc, df] = dF1_dp[fc, df] @@ -8446,7 +9272,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): # up_G3[fc*dim + df, dg] = dF1[fc, df] / d(dp/dx_dg) G3 = sympy.zeros(dim * dim, dim) for dg in range(dim): - dF1_dGp = diff_wrt_field(F1_for_jac, Gp[0, dg]) + dF1_dGp = diff_wrt_field(F1_for_up, Gp[0, dg]) for fc in range(dim): for df in range(dim): G3[fc * dim + df, dg] = dF1_dGp[fc, df] @@ -9478,6 +10304,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self._setup_region_ds() self.is_setup = True + self._form_facts = None # forms (re)built: cached symbolic facts are stale self.constitutive_model._solver_is_setup = True def set_active_region(self, region_label_name, region_label_value): @@ -9949,7 +10776,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): _force_setup: bool =False, time=None, divergence_retries: int = 0, - homotopy: bool = False, + homotopy=False, homotopy_options: dict = None, ): """ Solve the Stokes system for velocity and pressure. @@ -9970,9 +10797,23 @@ 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. + Picard warm-up: iterations with the FROZEN (Picard) tangent — linear + Stokes solves with the viscosity held at the current state — before + the consistent tangent takes over. Meaning by + :attr:`consistent_jacobian`: + + * ``True``: exactly ``picard`` Picard steps, then Newton. A cold start + takes ONE Picard step automatically; pass ``picard=-1`` to switch + that off. (The Jacobian is compiled in a blended Picard/Newton form + the first time this is used; switching between the two needs no + recompile.) + * ``"continuation"``: alpha is held at 0 (Picard) for the first + ``picard`` iterations; after that the residual-keyed ramp sets it + (see :attr:`consistent_jacobian`). + * ``False``: every iteration already is a Picard step, so this is + satisfied by the solve itself. + + Negative values mean "no warm-up". verbose : bool, default=False Print solver progress and timing information. debug : bool, default=False @@ -9990,21 +10831,24 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): If the final SNES solve reports DIVERGED, re-call it with warm start up to this many times. A single retry rescues most VEP yield-surface kink divergences. 0 preserves legacy behaviour. - homotopy : bool, default=False - Solve a yielding (viscoplastic) model by a **yield homotopy** instead of - a single solve: the constitutive model is put in its smooth, - δ-parameterised yield mode and δ is marched down to the sharp yield - surface as a sequence of warm-started solves. This is the robust route - for a hard Drucker-Prager problem, which does not converge from a cold - start on the sharp surface. Requires a model with - ``supports_yield_homotopy`` (raises otherwise). Returns the march - summary rather than ``None``. + homotopy : bool or "rate_strengthening", default=False + Solve a yielding (viscoplastic) model as a sequence of warm-started + solves instead of one. ``True``: the **yield homotopy** — the model is + put in its smooth, δ-parameterised yield mode and δ is marched down + toward the sharp yield surface; requires ``supports_yield_homotopy``. + ``"rate_strengthening"``: the ladder on the model's declared + ``Parameters.plastic_rate_strengthening`` (η_reg) — its multiplier + ``m`` runs from the viscous limit, found from the fields, down to + exactly 1, the problem as stated + (:func:`~underworld3.systems.yield_continuation.rate_strengthening_continuation`). + Any other value raises. Returns the march summary rather than ``None``. homotopy_options : dict, optional March settings passed to - :func:`~underworld3.systems.yield_continuation.yield_continuation` — - ``smoother``, ``anchor``, ``delta0``, ``down``, ``dmin``, - ``entry_maxit``, ``step_maxit``, ``retries``. All are defaulted; tuning - them is optional. + :func:`~underworld3.systems.yield_continuation.yield_continuation` + (``smoother``, ``anchor``, ``delta0``, ``down``, ``dmin``, + ``entry_maxit``, ``step_maxit``, ``retries``) or, for + ``"rate_strengthening"``, to ``rate_strengthening_continuation`` + (``scale0``, ``down``, ``rung_tangent``, ...). All are defaulted. ``smoother`` picks the soft-min family — ``"powermean"`` (default) or ``"sqrt"``. Which gives the better cold entry is problem-dependent, so @@ -10032,7 +10876,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 Picard steps, then Newton + >>> stokes.consistent_jacobian = True >>> stokes.solve(picard=3) >>> # Time-stepping with previous solution @@ -10048,9 +10893,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 -------- @@ -10068,9 +10914,17 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): # nothing before it is refused. self._reject_mixed_constraint_mechanisms("solve") + _check_homotopy_value(homotopy) if homotopy: - # The march runs a SEQUENCE of ordinary solves at successively sharper - # yield surfaces; each one re-enters this method with homotopy=False. + # The march runs a SEQUENCE of ordinary solves; each one re-enters this + # method with homotopy=False. True = the delta soft-min march; + # "rate_strengthening" = the ladder on the declared plastic rate + # strengthening (m from the viscous limit down to 1). + if homotopy == "rate_strengthening": + from underworld3.systems.yield_continuation import rate_strengthening_continuation + options = dict(homotopy_options or {}) + options.setdefault("verbose", verbose) + return rate_strengthening_continuation(self, **options) return self._solve_yield_homotopy(homotopy_options, verbose=verbose) if _force_setup: @@ -10085,7 +10939,20 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): # and resolving earlier would warm-start off the flag it just cleared. zero_init_guess = self._resolve_zero_init_guess(zero_init_guess) + # ---- Picard warm-up policy (#791) — decided BEFORE _build, because under the + # consistent tangent a Picard step needs the frozen operator compiled in. + # Rotated free-slip and fault contact run their own loop (utilities/rotated_bc.py) + # with their own picard semantics: no automatic step there. + self._plan_picard_warmup( + picard, zero_init_guess, + own_loop=bool(getattr(self, "_rotated_freeslip_bcs", None) + or getattr(self, "_fault_contact_faults", None)), + staged_loop=bool(getattr(self, "_rotated_freeslip_bcs", None))) + picard = max(int(picard or 0), 0) # picard < 0 = "no warm-up" from here on + self._build(verbose, debug, debug_name) + if self._picard_blend and self.consistent_jacobian is True: + self._set_newton_alpha(1.0) # Newton unless a warm-up says otherwise # Set time on the DM so petsc_t is available in pointwise functions. # Non-dimensionalise if the scaling system is active. @@ -10189,46 +11056,40 @@ 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 (successive substitution + # A(u_k) u_{k+1} = f). From a cold start it IS the viscous solve. This used to + # run SNES `nrichardson` with no nonlinear preconditioner — x <- x - lambda F(x), + # a residual step with no linear solve, nearly inert, NOT a Picard step. + # Semantics now, by tangent (policy decided before _build, above): # - # 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 - - if verbose and uw.mpi.rank == 0: - print(f"SNES solve - picard = {picard}", flush=True) - - # Picard solves if requested + # False — every iteration already is a Picard step, so `picard` + # is satisfied by the solve itself. + # "continuation" — alpha is held at 0 for the first `picard` iterations of + # the residual-keyed ramp in _continuation_solve. + # True — `picard` Picard steps, THEN Newton (_picard_warmup): the Jacobian is + # compiled in the blended form J_p + alpha (J_n - J_p) (alpha = 0 frozen, + # alpha = 1 Newton; a constants[] atom, no recompile to switch). A cold + # start takes ONE such step automatically (the Layer 1 design); + # picard = -1 switches it off. This matters most for boundary-driven + # problems, where Newton's first step from rest is NOT a Picard step + # (the Dirichlet values yield the driven layer at once). + # + # 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 fractional-power guard in _jacobian_unwrap, not by any + # warm-up. + # + # The policy (and the per-solve resets of _continuation_min_picard, + # _continuation_stages, _picard_stages) was applied by _plan_picard_warmup + # above. consistent_jacobian False: every iteration is already a Picard step. + # True: `_snes_solve_with_retries` takes the planned frozen-tangent steps first. - 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") + if verbose: + uw.pprint(f"SNES solve - picard = {picard} " + f"(tangent: {self.consistent_jacobian!r})") # The standard Newton solve, run whether or not the optional Picard # warmup above was taken: restore the configured SNES type and diff --git a/src/underworld3/model.py b/src/underworld3/model.py index b87841839..fb0a68bd7 100644 --- a/src/underworld3/model.py +++ b/src/underworld3/model.py @@ -1705,6 +1705,19 @@ def _record_solve_outcome(self, part: str, report) -> None: event["deadline_expired"] = True if getattr(report, "bounded", False): event["bounded"] = True + # Formulation facts the solve ran under (solve_report.config, #806): the + # tangent in force and, for Stokes, the pressure gauge and the yield-law + # regularisation (yield mode, softness, anchor, the declared rate + # strengthening and its ladder multiplier). These are what a reader + # needs to debug the MODEL from the transcript alone — which pressure a + # pressure-dependent yield law saw, whether a mean was projected out, + # which tangent and which law produced the step, one ladder rung from + # the next — without a live process to introspect. + config = getattr(report, "config", None) or {} + for key in ("tangent", "gauge", "regularisation"): + facts = config.get(key) + if facts: + event[key] = {k: v for k, v in facts.items() if not str(k).startswith("_")} return def _record_warning(self, message, category, filename, lineno) -> None: diff --git a/src/underworld3/systems/__init__.py b/src/underworld3/systems/__init__.py index 5ba37b93a..4d8522ef2 100644 --- a/src/underworld3/systems/__init__.py +++ b/src/underworld3/systems/__init__.py @@ -101,5 +101,6 @@ from .ddt import EulerianSUPG as EulerianSUPG_DDt # δ-continuation driver for hard viscoplastic (Drucker–Prager) yield -from .yield_continuation import yield_continuation, YieldHomotopyControl +from .yield_continuation import (yield_continuation, YieldHomotopyControl, + rate_strengthening_continuation, RateStrengtheningControl) from .solve_report import SolveReport diff --git a/src/underworld3/systems/advection_diffusion_eulerian.py b/src/underworld3/systems/advection_diffusion_eulerian.py index 8f7424fe4..813544d75 100644 --- a/src/underworld3/systems/advection_diffusion_eulerian.py +++ b/src/underworld3/systems/advection_diffusion_eulerian.py @@ -381,6 +381,10 @@ def _set_linear_solver(self, multigrid: bool): for key in self._SCHWARZ_OPTIONS: opts.delValue(key) self._push_managed_option("ksp_type", "gmres") + # the Schwarz route's restart 200 is gone with its keys; the multigrid + # route takes the solvers' default (see SolverBaseClass), not PETSc's 30 + from underworld3.cython.generic_solvers import _DEFAULT_KRYLOV_RESTART + self._push_managed_option("ksp_gmres_restart", _DEFAULT_KRYLOV_RESTART) for key, value in multigrid_options.gamg_bundle().settings.items(): self._push_managed_option(key, value) else: diff --git a/src/underworld3/systems/solve_report.py b/src/underworld3/systems/solve_report.py index f3861b856..85be887a8 100644 --- a/src/underworld3/systems/solve_report.py +++ b/src/underworld3/systems/solve_report.py @@ -135,6 +135,13 @@ class SolveReport: ``True`` when a wall-clock guard armed with ``solver.guard(...)`` cut this solve short. Distinguishes "the budget ran out" from a genuine divergence: both report ``DIVERGED_LINEAR_SOLVE``. + config + The numerical configuration the solve RAN with, snapshotted from the live PETSc + objects immediately before ``snes.solve`` (#806): ``snes`` (type, rtol, atol, stol, + max_it, Eisenstat-Walker), ``linesearch`` (type, order, tolerances), ``ksp`` (type, + tolerances, pc, and each fieldsplit sub-solver once set up), and ``tangent`` + (``consistent_jacobian``, alpha, continuation parameters). Two solves that should + be the same computation can be compared field by field. """ reason: int @@ -151,6 +158,12 @@ class SolveReport: bounded: bool = False sub: Mapping[str, "SubSolveReport"] = field(default_factory=dict) deadline_expired: bool = False + # The numerical configuration in force during the solve, read from the live PETSc + # objects immediately before ``snes.solve`` (#806): SNES/line-search/KSP types and + # tolerances, Eisenstat-Walker, fieldsplit sub-solvers, tangent mode and alpha. + # Captured AT the solve because some paths switch settings for its duration and + # restore them before the report is taken. Empty when no snapshot was taken. + config: Mapping[str, object] = field(default_factory=dict) def __str__(self) -> str: rho = f"{self.rho:.3f}" if self.rho is not None else "n/a" diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index aa0d6a5df..81511de3f 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -1921,7 +1921,7 @@ def solve( picard: int = 0, divergence_retries: int = 0, _skip_stress_history: bool = False, - homotopy: bool = False, + homotopy=False, homotopy_options: dict = None, ): """Solve the Stokes system, with optional viscoelastic stress history. @@ -1967,19 +1967,26 @@ def solve( kinks) to step off a bad Newton iterate. ``0`` preserves legacy behaviour (divergence is terminal). Typical useful value is 1. Only applies in the VE/VEP branch (``DFDt is not None``). - homotopy : bool, default=False - Solve a yielding model by marching the **yield homotopy** — the model's - δ-parameterised yield law is sharpened toward the exact ``Min`` over a - sequence of warm-started solves — instead of attempting the sharp surface - in one go. Requires ``constitutive_model.supports_yield_homotopy``. - Returns the march summary instead of ``None``. + homotopy : bool or "rate_strengthening", default=False + Solve a yielding model as a sequence of warm-started solves instead of + one. ``True`` marches the **yield homotopy** — the model's + δ-parameterised yield law is sharpened toward the exact ``Min`` — + and requires ``constitutive_model.supports_yield_homotopy``. + ``"rate_strengthening"`` ladders the model's declared + ``Parameters.plastic_rate_strengthening`` from the viscous limit down + to the problem as stated (``m = 1``). Any other value raises. Returns + the march summary instead of ``None``. homotopy_options : dict, optional March settings for :func:`~underworld3.systems.yield_continuation.yield_continuation` (``delta0``, ``down``, ``dmin``, ``entry_maxit``, ``step_maxit``, - ``retries``). All defaulted. + ``retries``) or + :func:`~underworld3.systems.yield_continuation.rate_strengthening_continuation`. + All defaulted. """ + from underworld3.cython.generic_solvers import _check_homotopy_value + _check_homotopy_value(homotopy) if homotopy: # Each δ-step re-enters this method with homotopy=False, so per-solve # arguments are forwarded to EVERY step of the march rather than dropped. @@ -2004,6 +2011,12 @@ def solve( "picard. Use homotopy_options={'entry_maxit': ...} to size the " "first solve." ) + if homotopy == "rate_strengthening": + from underworld3.systems.yield_continuation import rate_strengthening_continuation + options = dict(homotopy_options or {}) + options.setdefault("verbose", verbose) + options.setdefault("solve_kwargs", inner) + return rate_strengthening_continuation(self, **options) return self._solve_yield_homotopy( homotopy_options, verbose=verbose, solve_kwargs=inner ) diff --git a/src/underworld3/systems/yield_continuation.py b/src/underworld3/systems/yield_continuation.py index b6feeb10a..d612e0e8b 100644 --- a/src/underworld3/systems/yield_continuation.py +++ b/src/underworld3/systems/yield_continuation.py @@ -54,6 +54,301 @@ class YieldHomotopyControl: delta0: float = 1.0 +@dataclass(frozen=True) +class RateStrengtheningControl: + """How to ladder one constitutive model's declared plastic rate strengthening. + + Built by ``constitutive_model._rate_strengthening_control()``; ``None`` when the model + states no ``plastic_rate_strengthening``. The model owns the meaning (``eta_pl -> + tau_y/(2 edot_II) + m eta_reg``); the solver only moves ``m``, from the viscous limit + down to 1, and never past 1 — the end of the ladder is the problem as stated. + + Attributes + ---------- + set_scale + Sets the multiplier ``m`` (a ``constants[]`` atom: no recompile). + tangent + The ``consistent_jacobian`` value the model pairs with its yield law. + scale, parameter + The ``m`` atom and the ``eta_reg`` parameter, for diagnostics. + probe_terms + Cell-wise expressions the ladder reads to place its first rungs (see + ``ViscousFlowModel._rate_strengthening_probe_terms``), or ``None`` when the + model does not provide them; the ladder then needs an explicit ``scale0``. + """ + set_scale: Callable[[float], None] + tangent: Any + scale: Optional[Any] = None + parameter: Optional[Any] = None + probe_terms: Optional[dict] = None + + +def rate_strengthening_continuation( + solver, + control=None, + scale0=None, + scale_big=None, + entry_factor=1.0e3, + down=0.5, + rung_tangent="continuation", + entry_maxit=60, + step_maxit=80, + retries=2, + max_steps=40, + solve_kwargs=None, + verbose=True, +): + r"""Ladder the declared plastic rate strengthening from the viscous limit to the + stated problem (``m`` from ``scale0`` down to 1), each rung a full solve warm-started + from the previous one. + + Measured on the Spiegelman notch (2026-10-04): the equivalent ``xi`` ladder keeps the + velocity block healthy on every rung (~14 multigrid cycles per Krylov iteration) and + reaches a converged state on a case no single solve could; the blend + (``consistent_jacobian="continuation"``) on every rung is 3-4x cheaper than Newton + with a blend fallback, which is why it is the rung tangent here. + + ViscoPlasticFlowModel only: the elastic models advance a stress history on every + solve, so a ladder of solves has to be driven around one history update by hand. + + Parameters + ---------- + solver + A configured Stokes solver whose model states ``plastic_rate_strengthening``. + control : RateStrengtheningControl, optional + Defaults to ``solver.constitutive_model._rate_strengthening_control()``. + scale0 : float, optional + First rung. Default: the viscous limit found from the fields — one solve at the + entry multiplier (``scale_big``), at which nothing yields, then ``0.95 *`` the + largest ``m`` at which a cell still yields there. If that is ``<= 1`` the stated + problem does not yield beyond its onset and the ladder is one solve. + scale_big : float, optional + The entry multiplier. Default ``entry_factor * max(eta_ve / eta_reg)`` over the + domain, so that ``m * eta_reg`` exceeds the viscous branch everywhere whatever the + model's own ratio of the two scales. + down : float + Multiplicative step in ``m`` per rung (adapted: a cheap rung widens it, a + laboured one narrows it); the last rung is clamped to exactly 1. + rung_tangent + ``consistent_jacobian`` for the rungs where the model's tangent is Newton; + ``None`` keeps the model's own. + entry_maxit, step_maxit, retries, max_steps, solve_kwargs, verbose + As for :func:`yield_continuation`. + + Returns + ------- + dict + ``scale0``, ``rungs`` (list of ``(m, reason, its)``), ``reached_one``, + ``converged``, ``reason``, ``message``. A ladder that stops above ``m = 1`` leaves + the model at ``m = 1`` (as stated), the fields at the last converged rung, and + warns. + """ + import math + import warnings + + if not (0.0 < down < 1.0): + raise ValueError(f"down must satisfy 0 < down < 1, got {down}") + if getattr(getattr(solver, "Unknowns", None), "DFDt", None) is not None: + raise NotImplementedError( + "the rate-strengthening ladder runs several solves and would advance a stress " + "history each time; drive it yourself around one history update") + cm = getattr(solver, "constitutive_model", None) + if control is None: + control = cm._rate_strengthening_control() if cm is not None and hasattr(cm, "_rate_strengthening_control") else None + if control is None: + raise TypeError( + "rate_strengthening_continuation needs a constitutive model that STATES a " + "plastic_rate_strengthening (Parameters.plastic_rate_strengthening > 0); the " + "scale is the model's, not a solver default.") + if scale0 is None and control.probe_terms is None: + raise NotImplementedError( + f"{type(cm).__name__} does not provide the terms the ladder reads to find the " + f"viscous limit; pass scale0 explicitly") + + import numpy as np + + saved_probe = solver._difficulty_probe + saved_probe_maxit = solver._difficulty_max_it + saved_resume = solver._resume_abs_target + saved_tangent = solver.consistent_jacobian + saved_scale = cm.rate_strengthening_scale + saved_has_solution = bool(solver.has_solution) + u, p = solver.Unknowns.u, solver.Unknowns.p + # The state the ladder started from: restored if no rung converges, so a failed + # ladder never leaves a non-solution in the fields marked as one. + u_start, p_start = np.array(u.array), np.array(p.array) + if rung_tangent is not None: + solver.consistent_jacobian = rung_tangent if control.tangent is True else control.tangent + rungs = [] + state = dict(reason=0, reached_one=False, message=None, scale0=scale0) + + def _solve_rung(m, first, budget): + control.set_scale(m) + if first: + solver.is_setup = False + else: + solver._update_constants() + solver._difficulty_probe = True + solver._difficulty_max_it = budget + solver._resume_abs_target = None + solver.solve(zero_init_guess=not solver.has_solution, **dict(solve_kwargs or {})) + reason = int(solver.snes.getConvergedReason()) + rungs.append((float(m), reason, int(solver.snes.getIterationNumber()))) + return reason + + def _ladder(): + # returns with state["message"] set when it stops above m = 1 + nonlocal scale_big + if state["scale0"] is None: + terms = control.probe_terms + if scale_big is None: + visc, reg = _cell_values(solver, (terms["viscous"], terms["eta_reg"])) + scale_big = entry_factor * _masked_max_ratio(visc, reg) + if not (math.isfinite(float(scale_big)) and float(scale_big) > 0.0): + state["message"] = f"the entry multiplier is not a positive finite number (m = {scale_big})" + return + # the viscous limit, found from the fields rather than from the problem + state["reason"] = _solve_rung(float(scale_big), True, entry_maxit) + if state["reason"] <= 0: + state["message"] = (f"the viscous-limit solve (m={float(scale_big):g}) failed " + f"(reason={state['reason']})") + return + num, den = _cell_values(solver, (terms["onset_num"], terms["onset_den"])) + m_star = _masked_max_ratio(num, den) + if math.isnan(m_star) or m_star == math.inf: + state["message"] = f"the viscous limit could not be read from the fields (m* = {m_star})" + return + state["scale0"] = 0.95 * m_star if m_star > 1.0 else 1.0 + if verbose: + uw.pprint(f" [rate-strengthening] viscous limit m*={m_star:.4g}; " + f"ladder starts at m={state['scale0']:.4g}") + m = max(float(state["scale0"]), 1.0) + step = float(down) + u_good, p_good = None, None + if rungs and rungs[-1][1] > 0: + u_good, p_good = np.array(u.array), np.array(p.array) + failures = 0 + first = not rungs + while len(rungs) < max_steps: + budget = entry_maxit if first else step_maxit + reason = state["reason"] = _solve_rung(m, first, budget) + nit = rungs[-1][2] + first = False + if reason > 0: + u_good, p_good = np.array(u.array), np.array(p.array) + if verbose: + uw.pprint(f" [rate-strengthening] m={m:<10.4g} its={nit:3d} -> converged") + if m <= 1.0: + state["reached_one"] = True + return + if nit <= max(2, budget // 5): + step = max(step * step, 0.05) + elif nit >= 0.8 * budget: + step = min(step ** 0.5, 0.95) + m = max(m * step, 1.0) + continue + # a failed rung: back to the last converged state (or the start) + back_u, back_p = (u_good, p_good) if u_good is not None else (u_start, p_start) + with uw.synchronised_array_update("rate_strengthening revert"): + u.array[...] = back_u + p.array[...] = back_p + solver._record_convergence_status( + converged=True if u_good is not None else saved_has_solution) + failures += 1 + if failures > retries: + state["message"] = f"m={m:.4g} failed (reason={reason}) after {retries} retries" + return + step = min(step ** 0.5, 0.95) + # retry from the LAST converged rung — the state just restored + last_ok = next((r[0] for r in reversed(rungs) if r[1] > 0), m / step) + m = max(last_ok * step, 1.0) + if verbose: + uw.pprint(f" [rate-strengthening] failed (reason={reason}); retrying at m={m:.4g}") + state["message"] = (f"max_steps={max_steps} rungs used, last converged m=" + f"{next((r[0] for r in reversed(rungs) if r[1] > 0), float('nan')):.4g}") + + try: + _ladder() + finally: + solver._difficulty_probe = saved_probe + solver._difficulty_max_it = saved_probe_maxit + solver._resume_abs_target = saved_resume + solver.consistent_jacobian = saved_tangent + if not state["reached_one"]: + # leave the model as it was stated, never at an intermediate m; the next + # solve repacks the constant, so no repack is attempted here (it would + # mask the exception that brought us to this block) + control.set_scale(saved_scale) + if not state["reached_one"]: + if not any(r[1] > 0 for r in rungs): + # nothing converged: the fields go back to where the ladder started + with uw.synchronised_array_update("rate_strengthening restore"): + u.array[...] = u_start + p.array[...] = p_start + solver._record_convergence_status(converged=saved_has_solution) + message = state["message"] or "the ladder stopped above m = 1" + state["message"] = message + if verbose: + uw.pprint(f" [rate-strengthening] {message}; the model is back at m={saved_scale:g}") + warnings.warn( + f"[{solver.name}] rate-strengthening ladder did not reach the stated problem: " + f"{message}. The fields are the last converged rung's (or the starting state if " + f"none converged); the model is back at m={saved_scale:g}.", + RuntimeWarning, stacklevel=2) + return dict(scale0=state["scale0"], rungs=rungs, reached_one=state["reached_one"], + converged=bool(state["reached_one"] and state["reason"] > 0), + reason=state["reason"], message=state["message"]) + + +def _cell_values(solver, exprs): + """Cell averages (P0 projections) of ``exprs`` on this rank's cells. + + One P0 variable and one projector per solver, created on first use and reused, so + every rank creates the same objects in the same order. Raises if a projection + fails rather than returning the previous contents of the variable. + """ + import numpy as np + probe = getattr(solver, "_rate_strengthening_probe", None) + if probe is None: + mesh = solver.mesh + var = uw.discretisation.MeshVariable( + f"rs_probe_{solver.instance_number}", mesh, 1, degree=0, continuous=False) + proj = uw.systems.Projection(mesh, var) + proj.petsc_options.delValue("ksp_monitor") + probe = solver._rate_strengthening_probe = (var, proj) + var, proj = probe + # the probe compiles on its solver's JIT route, so "expanded" rules the graph out + # of the entry multiplier as well as of the solve + proj.jit_route = solver.jit_route + values = [] + for expr in exprs: + proj.uw_function = expr + proj.solve() + if proj.snes.getConvergedReason() <= 0: + raise RuntimeError( + f"rate-strengthening probe projection failed " + f"(reason={proj.snes.getConvergedReason()}) for {expr}") + values.append(np.asarray(var.array).reshape(-1).copy()) + return values + + +def _masked_max_ratio(num, den): + """Global ``max(num / den)`` over cells with ``den > 0`` (collective). + + ``-inf`` when no cell on any rank qualifies; NaN if any qualifying value is NaN. + """ + import numpy as np + from mpi4py import MPI + num, den = np.asarray(num, dtype=float), np.asarray(den, dtype=float) + ok = den > 0.0 + ratio = num[ok] / den[ok] + # NaN is reduced as its own flag: MPI.MAX compares with >, which drops a NaN + if int(uw.mpi.comm.allreduce(int(np.isnan(ratio).any()), op=MPI.MAX)): + return float("nan") + local = float(np.max(ratio)) if ratio.size else -np.inf + return float(uw.mpi.comm.allreduce(local, op=MPI.MAX)) + + def yield_continuation( solver, control=None, @@ -107,9 +402,12 @@ def yield_continuation( ``ViscousFlowModel.yield_anchor``. Applies to BOTH families. Default (``None``) leaves the model's own choice. ``"yield"`` keeps the stress exact at the yield point for every δ and puts the law on or above the exact one throughout; - ``"onset"`` (the historical default) is exact on the unyielded branch but sits + ``"onset"`` (the model default) is exact on the unyielded branch but sits BELOW the exact law at and above the yield point, which makes the entry problem - WEAKER than the sharp problem it is supposed to lead to. + WEAKER than the sharp problem it is supposed to lead to. The march therefore + runs under ``"onset"`` unless ``anchor="yield"`` is passed; for an approach + from above, pass it (or use the declared ``plastic_rate_strengthening`` + ladder, ``solve(homotopy="rate_strengthening")``). delta0 : float, optional Starting (smooth) δ. Defaults to the family's own entry (``control.delta0``), because δ is NOT the same parameter in the two families — see diff --git a/src/underworld3/utilities/_jit_graph.py b/src/underworld3/utilities/_jit_graph.py index 585246894..f2fd1f612 100644 --- a/src/underworld3/utilities/_jit_graph.py +++ b/src/underworld3/utilities/_jit_graph.py @@ -31,11 +31,22 @@ _nodes_made = False -def guard_half_integer_powers(e): - r"""``e`` with :math:`10^{-36}` added to the base of every half-integer power - whose base has free symbols: the sqrt guard of ``_jacobian_unwrap``, the same rule - and the same rebuild. Node applications are not entered (a guarded node's body was - guarded when it was made).""" +def is_fractional_power(n): + """True for a power with a non-integer exponent whose base has free symbols: its + value or its derivative can be singular where the base is zero, as the strain-rate + invariant is at rest. Half-integer powers are the invariant and its reciprocal; + SymPy merges a power of the bare invariant into others, ``sqrt(g)**(1/3)`` into + ``g**(1/6)`` (#841).""" + return n.is_Pow and n.exp.is_integer is not True and bool(n.args[0].free_symbols) + + +def guard_fractional_powers(e): + r"""``e`` with :math:`10^{-36}` added to the base of every fractional power + (:func:`is_fractional_power`): the guard of ``_jacobian_unwrap``, the same rule + and the same rebuild on both routes. In floating point the added term vanishes + against any base above about :math:`10^{-20}`, so only a state at rest sees it. + Node applications are not entered (a guarded node's body was guarded when it was + made).""" memo = {} def guard(n): @@ -53,8 +64,7 @@ def guard(n): if any(out == a and a != b for a, b in zip(args, new_args)): memo[id(n)] = (n, out) return out - if (out.is_Pow and out.exp.is_Rational and out.exp.q == 2 - and out.args[0].free_symbols): + if is_fractional_power(out): out = sympy.Pow(out.args[0] + _EPS2, out.exp) memo[id(n)] = (n, out) return out @@ -332,7 +342,7 @@ def node_of(self, atom, guarded): try: body = self.lower(atom.sym, guarded) if guarded: - body = guard_half_integer_powers(body) + body = guard_fractional_powers(body) # a matrix-valued atom cannot be one scalar temporary: it is expanded # in place, as the tree route expands it out = self.make_node(body, label=atom.name) if _is_scalar(body) else body diff --git a/src/underworld3/utilities/custom_mg.py b/src/underworld3/utilities/custom_mg.py index 4c57b1ec0..23ac11eb0 100644 --- a/src/underworld3/utilities/custom_mg.py +++ b/src/underworld3/utilities/custom_mg.py @@ -1059,6 +1059,13 @@ def _configure_pcmg(pc, Ps, coarse="redundant", smoother="robust", owned=None, # ``PCSetUp``, after the write. if ksp is not None: ksp.setType(opts.getString(prefix + "ksp_type", ksp.getType())) + # KSPSetType re-creates the Krylov context with ITS defaults, so a restart the + # options database carries (the solvers' default 100, or the user's own) is + # lost here: re-apply it to the live object (measured: a single-field FMG + # solve ran FGMRES(30) while its report said 100). + if ksp.getType() in ("gmres", "fgmres", "lgmres", "dgmres", "pgmres", "pipefgmres") \ + and opts.hasName(prefix + "ksp_gmres_restart"): + ksp.setGMRESRestart(opts.getInt(prefix + "ksp_gmres_restart")) pc.setType("mg") pc.setMGLevels(nlev) pc.setMGType(PETSc.PC.MGType.FULL) diff --git a/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 81ccd4571..ec14fa256 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -983,13 +983,14 @@ def solve_rotated_freeslip(solver, boundaries, remove_rotation_gauge=True, the loose ``solver.newton_switch_rtol`` (and at least ``picard`` iterations if given) to enter Newton's basin; phase 2 sets α=1 (consistent Newton tangent) and drives to the requested tolerance. α is restored to 0 afterwards. - * ``consistent_jacobian is True`` (pure Newton) with ``picard > 0`` — a Picard - warmup needs the frozen tangent, which the pure-Newton compile does not carry, so - this **raises** ``NotImplementedError`` pointing to ``"continuation"`` rather than - silently ignoring the request. - * ``consistent_jacobian is False`` (default, frozen/Picard) — the whole solve is - already the frozen tangent, so ``picard`` is inherently satisfied (matches the - standard path, whose post-warmup Newton phase also uses the frozen tangent here). + * ``consistent_jacobian is True`` (the default, consistent Newton) with + ``picard > 0`` — the same two-phase solve on the blended Jacobian that + ``solve()`` compiles for an explicit warm-up. A cold start takes NO automatic + Picard step on this path (pure Newton from rest converges here). Called directly, + without that compile, a nonlinear residual **raises** ``NotImplementedError`` + rather than silently ignoring the request. + * ``consistent_jacobian is False`` (frozen/Picard) — the whole solve is already the + frozen tangent, so ``picard`` is inherently satisfied. Returns ------- @@ -1051,10 +1052,15 @@ def solve_rotated_freeslip(solver, boundaries, remove_rotation_gauge=True, rtol = float(solver.tolerance) # Resolve the tangent / warmup policy (see the docstring). ``continuation`` runs a - # staged α=0 → α=1 solve; pure Newton + picard>0 is unsupported (no frozen tangent - # to warm up with) and errors loudly; frozen (default) already satisfies a warmup. + # staged α=0 → α=1 solve. Under the consistent tangent (the default) an explicit + # picard>0 runs the same staging on the blended Jacobian that solve() compiled for + # it (_plan_picard_warmup, staged_loop); the frozen tangent already satisfies a + # warmup. No automatic cold-start step is taken here: from rest, pure Newton on + # this path converges (measured 8-9 its against 27-29 frozen on a VP box). mode = getattr(solver, "consistent_jacobian", False) - continuation = (mode == "continuation") + continuation = (mode == "continuation" + or (mode is True and bool(picard) and picard > 0 + and getattr(solver, "_picard_blend", False))) if picard and not continuation and mode is True: # The ONLY place linearity must be known in advance (hence the lazy probe — # two trial assemblies, paid only by this contradictory configuration): a diff --git a/tests/parallel/test_1081_cold_start_verdict_is_rank_uniform.py b/tests/parallel/test_1081_cold_start_verdict_is_rank_uniform.py new file mode 100644 index 000000000..50d905dfd --- /dev/null +++ b/tests/parallel/test_1081_cold_start_verdict_is_rank_uniform.py @@ -0,0 +1,45 @@ +#!/usr/bin/env python3 +"""The cold-start verdict that arms the automatic Picard warm-up must be the same on +every rank. + +`_solution_is_trivially_zero()` decides whether a solve that was told to warm-start is +really a cold start, and the verdict arms a JIT rewire and an extra `snes.solve` — both +collective. It read `u.vec.norm()`, and `u.vec` is the variable's LOCAL vector, so the +norm was per rank. Measured at np=2 (review of PR #794): a nonlinear Poisson solve with +an initial guess non-zero in one corner only — rank 0 saw a zero norm and took the +warm-up, rank 1 did not; rank 1 failed with MPI_ERR_BUFFER in the Jacobian assembly and +rank 0 spun until killed. + +The verdicts are gathered and compared on every rank BEFORE the solve, so a regression +fails the test on all ranks together instead of hanging the job. +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + + +def test_trivially_zero_is_rank_uniform_and_the_solve_completes(): + uw.reset_default_model() + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.1) + u = uw.discretisation.MeshVariable("u1081", mesh, 1, degree=1) + pois = uw.systems.Poisson(mesh, u_Field=u) + pois.constitutive_model = uw.constitutive_models.DiffusionModel + pois.constitutive_model.Parameters.diffusivity = 1.0 + u.sym[0] ** 2 + pois.f = 1.0 + for b in ("Left", "Right", "Top", "Bottom"): + pois.add_dirichlet_bc(0.0, b) + coords = u.coords + corner = (coords[:, 0] < 0.15) & (coords[:, 1] < 0.15) + with uw.synchronised_array_update(): + u.array[:, 0, 0] = np.where(corner, 0.1, 0.0) + + verdicts = uw.mpi.comm.allgather(bool(pois._solution_is_trivially_zero())) + assert verdicts == [False] * uw.mpi.size, verdicts + + pois.solve(zero_init_guess=False) + assert pois.snes.getConvergedReason() > 0 diff --git a/tests/test_0022_unwrap_memoised_matches_fixed_point.py b/tests/test_0022_unwrap_memoised_matches_fixed_point.py index ef02864cb..064999001 100644 --- a/tests/test_0022_unwrap_memoised_matches_fixed_point.py +++ b/tests/test_0022_unwrap_memoised_matches_fixed_point.py @@ -103,14 +103,14 @@ def test_memoised_unwrap_is_the_fixed_point(laws, mode): def test_the_jacobian_sqrt_guard_matches_replace(laws): """The expanded JIT route's Newton source: its guard, a memoised rebuild, against - SymPy's own replace.""" + SymPy's own replace with the same rule (every fractional power, #841).""" from underworld3.cython.generic_solvers import _jacobian_unwrap eps2 = sympy.Float(1.0e-36) def old_guard(e): return e.replace( - lambda n: (n.is_Pow and n.exp.is_Rational and n.exp.q == 2 + lambda n: (n.is_Pow and n.exp.is_integer is not True and n.args[0].free_symbols), lambda n: sympy.Pow(n.args[0] + eps2, n.exp)) @@ -153,21 +153,21 @@ def test_the_shared_xreplace_matches_xreplace(laws): def test_the_graph_guard_matches_replace(laws): - """The guard every Newton node body gets (``_jit_graph.guard_half_integer_powers``), + """The guard every Newton node body gets (``_jit_graph.guard_fractional_powers``), a memoised rebuild, against SymPy's own ``replace`` with the same rule.""" - from underworld3.utilities._jit_graph import guard_half_integer_powers + from underworld3.utilities._jit_graph import guard_fractional_powers eps2 = sympy.Float(1.0e-36) def old_guard(e): return e.replace( - lambda n: (n.is_Pow and n.exp.is_Rational and n.exp.q == 2 + lambda n: (n.is_Pow and n.exp.is_integer is not True and n.args[0].free_symbols), lambda n: sympy.Pow(n.args[0] + eps2, n.exp)) for name, expr in laws: tree = _unwrap_each(expr, lambda e: _fixed_point_unwrap(e, "symbolic_keep_constants")) - new = _unwrap_each(tree, guard_half_integer_powers) + new = _unwrap_each(tree, guard_fractional_powers) old = _unwrap_each(tree, old_guard) assert sympy.srepr(new) == sympy.srepr(old), name diff --git a/tests/test_0024_jit_graph_lowering.py b/tests/test_0024_jit_graph_lowering.py index 8da9edac3..e289a95dd 100644 --- a/tests/test_0024_jit_graph_lowering.py +++ b/tests/test_0024_jit_graph_lowering.py @@ -148,9 +148,10 @@ def test_a_changed_body_is_lowered_afresh_without_clearing_the_cache(box): def _guarded_tree(e): - """The expanded-tree Newton source the graph replaced: every non-constant atom - expanded, then 1e-36 added to the base of every half-integer power with free - symbols, by SymPy's own replace.""" + """The expanded-tree Newton source the graph replaced, with the guard it had then: + every non-constant atom expanded, then 1e-36 added to the base of every + half-integer power with free symbols, by SymPy's own replace. (Both routes have + guarded every fractional power since #841.)""" eps2 = sympy.Float(1.0e-36) tree = ex.unwrap_expression(e, mode="symbolic_keep_constants") return tree.replace( diff --git a/tests/test_0104_constant_slot_still_constant.py b/tests/test_0104_constant_slot_still_constant.py index 8a0f1d86e..d1e2df2fe 100644 --- a/tests/test_0104_constant_slot_still_constant.py +++ b/tests/test_0104_constant_slot_still_constant.py @@ -42,6 +42,12 @@ def _build(initial, route="graph"): poisson.add_dirichlet_bc(0.0, "Top") poisson.add_dirichlet_bc(0.0, "Bottom") poisson.petsc_options.delValue("ksp_monitor") + # Frozen tangent on purpose: under the consistent default (2026-10-05) the first solve + # whose flux has become field-dependent takes the automatic Picard warm-up, and + # compiling the blended Jacobian for it rewires the kernel — which repacks the slot + # correctly and so never reaches the staleness check this file pins. The check is + # about the constants manifest, not the tangent; hold the tangent still. + poisson.consistent_jacobian = False poisson.jit_route = route return uw, poisson, T, m diff --git a/tests/test_0201_solver_has_solution_warmstart.py b/tests/test_0201_solver_has_solution_warmstart.py index 8ce6cc907..77d0ff142 100644 --- a/tests/test_0201_solver_has_solution_warmstart.py +++ b/tests/test_0201_solver_has_solution_warmstart.py @@ -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; a LINEAR flux takes no automatic Picard step (a nonlinear one takes one). (Before #791 that step was an + nrichardson residual sweep; 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,10 @@ 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 converges and sets has_solution. Since + 2026-10-05 the automatic Picard warm-up is taken only when the flux's coefficient + depends on the unknowns; this viscosity is constant, so no warm-up runs (the frozen + and consistent tangents coincide). The nonlinear case is pinned in test_1068.""" mesh = uw.meshing.StructuredQuadBox( elementRes=(8, 8), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0) ) @@ -202,6 +205,8 @@ 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, linear flux: Newton directly, no Picard step + assert stokes._picard_stages is None, ( + f"a linear flux took an automatic Picard step: {stokes._picard_stages}") assert stokes.snes.getConvergedReason() > 0 assert stokes.has_solution is True diff --git a/tests/test_0641_wave_c_api_shims.py b/tests/test_0641_wave_c_api_shims.py index b143ef893..821fec267 100644 --- a/tests/test_0641_wave_c_api_shims.py +++ b/tests/test_0641_wave_c_api_shims.py @@ -250,9 +250,9 @@ def test_invalid_values_raise(self, mesh, stokes): with pytest.raises(ValueError, match="consistent_jacobian"): stokes.consistent_jacobian = value - def test_default_is_false(self, mesh): + def test_default_is_true(self, mesh): solver = uw.systems.Poisson(mesh) - assert solver.consistent_jacobian is False + assert solver.consistent_jacobian is True # --------------------------------------------------------------------------- diff --git a/tests/test_1018_rotated_freeslip.py b/tests/test_1018_rotated_freeslip.py index e137463ce..6a7f56c30 100644 --- a/tests/test_1018_rotated_freeslip.py +++ b/tests/test_1018_rotated_freeslip.py @@ -918,7 +918,8 @@ def _powerlaw_stokes(mesh, prefix, amp=2.0, nexp=3.0, cj=None): """A genuinely NONLINEAR Stokes: power-law viscosity eta = eps_II^(1/n - 1) (smooth, so Newton/Picard iterates robustly), driven by a horizontally-varying vertical body force so there is real shear. ``cj`` sets ``consistent_jacobian`` - (None=default frozen/Picard, True=consistent Newton, "continuation"=staged).""" + (None = the solver default, consistent Newton since 2026-10-05; False = frozen/ + Picard; True = consistent Newton; "continuation" = staged).""" x, y = mesh.X v = uw.discretisation.MeshVariable(prefix + "v", mesh, mesh.dim, degree=2, continuous=True) p = uw.discretisation.MeshVariable(prefix + "p", mesh, 1, degree=1, continuous=False) @@ -962,12 +963,14 @@ def plaw_box_ref(): def test_rotated_freeslip_nonlinear_matches_essential(plaw_box_ref): - """Default (frozen/Picard) tangent through the rotated path: genuinely iterates - and converges to the native essential nonlinear free-slip answer (both impose - v_n=0 — identical discrete problem), with machine-zero wall-normal flow on every - wall. Exercises the rotated constraint INSIDE the nonlinear iteration.""" + """Frozen (Picard) tangent through the rotated path: genuinely iterates and + converges to the native essential nonlinear free-slip answer (both impose v_n=0 — + identical discrete problem), with machine-zero wall-normal flow on every wall. + Exercises the rotated constraint INSIDE the nonlinear iteration. Pinned to + consistent_jacobian=False: the default is Newton since 2026-10-05, and this is the + rotated path's only frozen-tangent nonlinear test.""" mesh, vE, _ = plaw_box_ref - s, vR, pR = _powerlaw_stokes(mesh, "nlP") # default (Picard) tangent + s, vR, pR = _powerlaw_stokes(mesh, "nlP", cj=False) # frozen (Picard) tangent for wall in ("Top", "Bottom", "Left", "Right"): s.add_rotated_freeslip_bc(0, wall) s.solve() @@ -1028,17 +1031,23 @@ def test_rotated_freeslip_continuation_tangent(plaw_box_ref): "picard did not extend the continuation Picard phase") -def test_rotated_freeslip_picard_newton_unsupported_raises(plaw_box_ref): - """picard>0 with the pure consistent-Newton tangent has no frozen warmup tangent - to form, so the rotated path raises a clear NotImplementedError pointing to - 'continuation' rather than silently ignoring it. Also confirms the SNES_Stokes - wrapper forwards picard to the plain-Stokes solve at all.""" - mesh, _, _ = plaw_box_ref - s, v, p = _powerlaw_stokes(mesh, "prN", cj=True) +def test_rotated_freeslip_picard_under_newton_stages_picard_then_newton(plaw_box_ref): + """picard>0 under the consistent-Newton tangent (the default since 2026-10-05) stages + at least N Picard iterations, then Newton, on the blended Jacobian solve() compiles + for it, and converges to the essential-BC answer. It used to RAISE + NotImplementedError, which made solve(picard=N) — the guides' recipe — fail on + every nonlinear rotated free-slip model under the new default. Also confirms the + SNES_Stokes wrapper forwards picard to the plain-Stokes solve at all (without it + there is no Picard phase to switch from).""" + mesh, vE, _ = plaw_box_ref + s, vR, pR = _powerlaw_stokes(mesh, "prN", cj=True) for wall in ("Top", "Bottom", "Left", "Right"): s.add_rotated_freeslip_bc(0, wall) - with pytest.raises(NotImplementedError, match="continuation"): - s.solve(picard=3) + s.solve(picard=3) + info = s._rotated_freeslip_info + assert info["continuation_switched"], "no Picard phase: picard was not staged" + assert info["nonlinear_iterations"] >= 3, info["nonlinear_iterations"] + assert np.linalg.norm(vR.data - vE) / np.linalg.norm(vE) < 1e-6 def test_rotated_freeslip_nonlinear_warm_start(plaw_box_ref): diff --git a/tests/test_1055_yield_smoother.py b/tests/test_1055_yield_smoother.py index 580b29df4..b17290cb8 100644 --- a/tests/test_1055_yield_smoother.py +++ b/tests/test_1055_yield_smoother.py @@ -157,8 +157,11 @@ def test_yield_softness_runtime_no_recompile(): @pytest.mark.level_1 @pytest.mark.tier_a def test_dp_model_smoother_optin(): - """The non-elastic Drucker–Prager model defaults to exact hard Min and can opt - into the δ soft-min / power-mean homotopy (its new capability).""" + """The non-elastic Drucker–Prager model defaults to the SMOOTH yield law (since + 2026-10-05: softmin, sqrt, anchor "onset", delta 0.1 — the viscous branch and the + plastic limit exact, the yield point itself 4.5% under-stressed); the exact hard + Min and the power-mean homotopy are explicit choices. Hard baselines: the closed + form eta_ve / g(f), g = 1 + (f - 1 + sqrt((f-1)^2 + d^2))/2 - (sqrt(1 + d^2) - 1)/2.""" mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) v = uw.discretisation.MeshVariable("Udp", mesh, mesh.dim, degree=2) p = uw.discretisation.MeshVariable("Pdp", mesh, 1, degree=1) @@ -166,9 +169,24 @@ def test_dp_model_smoother_optin(): cm = uw.constitutive_models.ViscoPlasticFlowModel(s.Unknowns) s.constitutive_model = cm - # default: exact hard Min - assert cm.yield_mode == "min" + # default: the smooth law pinned at the unyielded limit + assert cm.yield_mode == "softmin" assert cm.yield_smoother == "sqrt" + assert cm.yield_anchor == "onset" + assert cm.yield_softness == 0.1 + for eta_ve, eta_pl, expected in [ + (1.0, 1.0e6, 0.9999999975185925), # unyielded: the viscous answer + (1.0, 2.0, 0.9975488283821262), # just below yield: 0.25% soft + (1.0, 1.0, 0.9546482702586252), # the yield point: 4.5% under-stressed + (2.0, 1.0, 1.0), # f = 2: exact + (1.0, 1.0e-3, 0.00100000249128476), # deep yield: the plastic limit + ]: + comb = cm._combine_yield(sympy.Float(eta_ve), sympy.Float(eta_pl)) + val = float(unwrap_expression(comb, mode="nondimensional")) + assert abs(val - expected) <= 1.0e-12 * expected, (eta_ve, eta_pl, val, expected) + + # the exact hard Min remains an explicit choice + cm.yield_mode = "min" for eta_ve, eta_pl in [(1.0, 2.0), (2.0, 1.0), (0.3, 5.0)]: comb = cm._combine_yield(sympy.Float(eta_ve), sympy.Float(eta_pl)) val = float(unwrap_expression(comb, mode="nondimensional")) @@ -176,7 +194,14 @@ def test_dp_model_smoother_optin(): # opt into the power-mean homotopy: undershoots Min under the onset anchor cm.yield_mode = "softmin" - cm.yield_anchor = "onset" # the side under test; not the default + cm.yield_anchor = "onset" # the side under test (also the default) + # δ carries across a family switch: from the default 0.1 the power mean gets + # sharpness s = 1/(0.1 + 0.001) ~ 9.9, not its entry value 1.0 (set that explicitly) + cm.yield_smoother = "powermean" + assert cm.yield_softness == 0.1 + cm.yield_smoother = "sqrt" + # only the singular δ = 0 is bumped, to the power mean's entry value + cm.yield_softness = 0.0 cm.yield_smoother = "powermean" # bumps δ→1 assert cm.yield_softness == 1.0 for eta_ve, eta_pl in [(1.0, 2.0), (5.0, 0.3), (1e25, 1e21)]: diff --git a/tests/test_1057_yield_homotopy_solve.py b/tests/test_1057_yield_homotopy_solve.py index 8bd74fd64..a1da222db 100644 --- a/tests/test_1057_yield_homotopy_solve.py +++ b/tests/test_1057_yield_homotopy_solve.py @@ -262,19 +262,25 @@ def _yielding_box(tag, tau_y, cellSize=0.2): @pytest.mark.level_2 def test_homotopy_rescues_a_solve_the_cold_start_cannot_do(): - """THE user-level guarantee: ``solve(homotopy=True)`` converges on a genuinely - yielding problem where a direct cold solve of the sharp law does not. - - This is a CAPABILITY test, deliberately asserting both halves on the same problem, - so the feature cannot silently regress into "runs without error but no longer - rescues anything". At tau_y = 0.30 (~45% of the domain yielding) the cold hard-Min - solve gives DIVERGED_MAX_IT; the march settles near 1e-4 and converges. + """``solve(homotopy=True)`` converges on a genuinely yielding problem where a direct + cold solve of the sharp law under the FROZEN (Picard) tangent does not. + + A CAPABILITY test, asserting both halves on the same problem, so the feature cannot + silently regress into "runs without error but no longer rescues anything". At + tau_y = 0.30 (~45% of the domain yielding) the cold hard-Min Picard solve gives + DIVERGED_MAX_IT; the march settles near 1e-4 and converges. Since 2026-10-05 the + DEFAULT cold solve (one Picard step, then Newton) converges on this box too, so the + rescue is measured against the Picard-throughout solve, not the default one. """ import numpy as np - # (a) the direct cold solve of the sharp law FAILS + # (a) the direct cold solve of the sharp law under the FROZEN tangent FAILS. Pinned to + # consistent_jacobian=False: with the Newton default (2026-10-05, one automatic cold + # Picard step then Newton) this cold hard-Min solve CONVERGES (reason 3), so the + # baseline the homotopy is measured against is the Picard-throughout solve. mesh, cold, cm_cold = _yielding_box("c", 0.30) cm_cold.yield_mode = "min" + cold.consistent_jacobian = False cold.solve() cold_reason = int(cold.snes.getConvergedReason()) diff --git a/tests/test_1059_yield_anchor.py b/tests/test_1059_yield_anchor.py index 7747d12d2..37d3cd0bc 100644 --- a/tests/test_1059_yield_anchor.py +++ b/tests/test_1059_yield_anchor.py @@ -114,6 +114,8 @@ def test_yield_anchor_removes_the_power_mean_degenerate_basin(): @pytest.mark.level_1 @pytest.mark.tier_a def test_anchor_rejects_nonsense_and_defaults_to_onset(): + # "onset" stays the default with the smooth default law (2026-10-05): "yield" made + # every unyielded cell 1/(1 - delta/2) stiffer and the yielded tangent negative c = _model() assert c.yield_anchor == "onset" with pytest.raises(ValueError): diff --git a/tests/test_1066_stokes_jacobian_layout.py b/tests/test_1066_stokes_jacobian_layout.py index 976444caa..95ea75d62 100644 --- a/tests/test_1066_stokes_jacobian_layout.py +++ b/tests/test_1066_stokes_jacobian_layout.py @@ -15,7 +15,7 @@ 2. ``TransverseIsotropicFlowModel._build_c_tensor`` baked the UNWRAPPED ``.sym`` contents of its viscosity Parameters into the c-tensor, exposing the strain-rate dependence to ``sympy.diff`` and silently un-freezing the - default (Picard) tangent. Tests 2/3 pin the freezing contract: coefficients + frozen (Picard) tangent. Tests 2/3 pin the freezing contract: coefficients hide inside UWexpression atoms under Picard; only the Newton unwrap sees through them. @@ -166,7 +166,7 @@ def test_stokes_uu_g3_matches_fd_oracle(): @pytest.mark.level_1 @pytest.mark.tier_a def test_ti_picard_tangent_is_frozen(): - """Under the default (Picard) tangent the TI viscosity coefficients must + """Under the frozen (Picard) tangent the TI viscosity coefficients must stay hidden inside UWexpression atoms: no raw velocity-gradient (Derivative) atoms may appear in the assembled uu_G3. The Newton build of the same problem MUST expose them (positive control that the dependence @@ -188,7 +188,7 @@ def test_ti_picard_tangent_is_frozen(): ] assert not picard_live, ( "Picard uu_G3 depends on the velocity gradient — the TI c-tensor " - "build has un-frozen the default tangent (issue #457): first hits " + "build has un-frozen the Picard tangent (issue #457): first hits " f"{picard_live[:4]}" ) @@ -371,7 +371,7 @@ def test_ti_vep_viscosity_property_reports_yield_limited_eta1(): def test_vep_c_tensor_coefficients_are_frozen(): """Isotropic VEP with yield active: same freezing contract as the TI classes (the .sym tidy deferred on PR #493) — the c-tensor coefficients - must live inside UWexpression atoms so the default (Picard) tangent stays + must live inside UWexpression atoms so the frozen (Picard) tangent stays frozen; only the Newton unwrap sees the yield law's strain-rate dependence.""" import itertools diff --git a/tests/test_1067_newton_cold_start.py b/tests/test_1067_newton_cold_start.py index 3d9165dd0..849842fc8 100644 --- a/tests/test_1067_newton_cold_start.py +++ b/tests/test_1067_newton_cold_start.py @@ -8,9 +8,10 @@ (error 77) on the standard path, the rotated free-slip path, and the split-node fault path alike. The alpha-blended continuation kernel inherited it even in its Picard phase (IEEE 0*NaN = NaN). The guard in -``_jacobian_unwrap`` regularises exactly the half-integer-power family -(+1e-36 under the root), leaving the residual and the default Picard -tangent bit-identical. +``_jacobian_unwrap`` regularised the half-integer-power family (+1e-36 +under the root), leaving the residual and the frozen (Picard) tangent +bit-identical. Since #841 it covers every fractional power: SymPy merges a +power of the bare invariant into one that is not a half-integer. """ import numpy as np import pytest @@ -53,3 +54,53 @@ def test_cold_consistent_tangent_is_finite_and_solves(tangent): U.set(0.0) snes.computeJacobian(U, J) assert np.isfinite(J.norm()), "cold-state Jacobian contains NaN/Inf" + + +def _bounded_thinning_box(route, tangent): + """A bounded shear-thinning law written on the bare strain-rate invariant, + eta = 1 / (1 + edot**(1/3)). The residual is finite at rest. SymPy merges + sqrt(g)**(1/3) into g**(1/6), a power the half-integer guard did not cover, so + the Newton tangent was 0 * Inf = NaN at rest, on both JIT routes; under the + default tangent the cold start's Picard step runs the blended kernel at + alpha = 0, and 0 * NaN is NaN there too (#841).""" + import sympy + + tag = f"{route[0]}{'n' if tangent is True else 'p'}" + mesh = uw.meshing.UnstructuredSimplexBox(cellSize=0.25) + v = uw.discretisation.MeshVariable(f"V1067{tag}", mesh, 2, degree=2) + p = uw.discretisation.MeshVariable(f"P1067{tag}", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + edot = stokes.Unknowns.Einv2 + stokes.constitutive_model.Parameters.shear_viscosity_0 = uw.expression( + rf"\eta_{{1067{tag}}}", 1 / (1 + edot ** sympy.Rational(1, 3)), + "bounded shear thinning") + stokes.add_essential_bc((1.0, 0.0), "Top") + stokes.add_essential_bc((0.0, 0.0), "Bottom") + stokes.petsc_options.delValue("ksp_monitor") + stokes.tolerance = 1.0e-8 + stokes.jit_route = route + stokes.consistent_jacobian = tangent + return stokes, v + + +@pytest.mark.parametrize("route", ["graph", "expanded"]) +def test_a_fractional_power_of_the_invariant_solves_cold_under_newton(route): + """Newton from a cold start reaches the solution the frozen tangent reaches, and + its Jacobian at rest is finite.""" + frozen, v_frozen = _bounded_thinning_box(route, False) + frozen.solve() + assert frozen.snes.getConvergedReason() > 0 + + newton, v_newton = _bounded_thinning_box(route, True) + newton.solve() # cold: raised error 77 before the guard + assert newton.snes.getConvergedReason() > 0 + + a, b = np.asarray(v_newton.data), np.asarray(v_frozen.data) + assert np.abs(a - b).max() <= 1.0e-6 * np.abs(b).max() + + J = newton.snes.getJacobian()[0].copy() + U = newton.dm.getGlobalVec() + U.set(0.0) + newton.snes.computeJacobian(U, J) + assert np.isfinite(J.norm()), "the Newton Jacobian at rest is not finite" 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..940329ea7 --- /dev/null +++ b/tests/test_1068_picard_warmup_is_a_frozen_tangent_step.py @@ -0,0 +1,467 @@ +"""A Picard step is a Newton iteration with the FROZEN tangent — not a residual sweep (#791). + +The recipe for hard viscoplastic solves is "a few Picard steps, then Newton", and that must be one +flag: ``solve(picard=N)``. A Picard step is successive substitution, A(u_k) u_{k+1} = f — a linear +Stokes solve (velocity block and Schur complement) with the viscosity frozen at the current state. +The standard path used to run 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 contract, per tangent mode: + + * ``consistent_jacobian=True`` — ``picard=N`` gives exactly N Picard steps, then Newton; a cold + start takes ONE automatically (``picard=-1`` switches it off). Held against the definition: + after ``picard=1`` with the Newton stage capped at zero iterations, the state must equal ONE + frozen-tangent iteration — on a BOUNDARY-DRIVEN problem, where Newton's first step from rest is + NOT a Picard step (measured 45% apart), so the test cannot pass by coincidence. + * ``"continuation"`` — ONE solve with alpha keyed on the residual relative to F_ref = ||F(u=0)||; + ``picard=N`` holds alpha = 0 for the first N iterations, so the frozen-tangent count is + max(N, natural). A warm start enters the ramp at the alpha its residual calls for. + * ``False`` — every iteration is already a Picard step, so ``picard`` changes nothing. +""" +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 + # These fixtures pin the frozen-tangent semantics on the EXACT hard Min, the law whose + # kink makes the warm-up necessary. Since 2026-10-05 the model default is the smooth + # law (softmin, sqrt, anchor "onset", delta 0.1), under which pure Newton converges here. + cm.yield_mode = "min" + 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 _driven_box(tag, tangent, tau_y=0.30): + """Sheared by a moving top wall: boundary-driven, so the rest state is NOT zero — the + Dirichlet values yield the driven layer at once, and the consistent tangent at rest + differs from the frozen one.""" + uw.reset_default_model() + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.25) + v = uw.discretisation.MeshVariable("Vdb" + tag, mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("Pdb" + tag, mesh, 1, degree=1, continuous=True) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + s.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel + s.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + s.constitutive_model.Parameters.yield_stress = tau_y + s.constitutive_model.yield_mode = "min" # exact hard Min, see _yielding_box + s.add_essential_bc((1.0, 0.0), "Top") + s.add_essential_bc((0.0, 0.0), "Bottom") + s.petsc_use_pressure_nullspace = True + s.petsc_options.delValue("ksp_monitor") + s.consistent_jacobian = tangent + s.tolerance = 1.0e-6 + s.petsc_options["snes_max_it"] = 300 + return s, v + + +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 holds alpha = 0 for the first N iterations of the ramp, so the frozen-tangent + count is max(N, natural). Hard baselines: + + * N below the natural count -> exactly the natural count (picard adds nothing); + * N above it -> exactly N. + + And all three land on the same solution.""" + 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_newton_warmup_step_is_exactly_a_picard_step(): + """THE DEFINITION. Under consistent_jacobian=True, picard=1 followed by a Newton stage + capped at zero iterations must leave exactly the state of ONE frozen-tangent iteration. + Boundary-driven, so an nrichardson sweep (old code) or no warm-up (first fix) both differ.""" + sN, vN = _driven_box("n", True) + sN.petsc_options["snes_max_it"] = 0 # Newton stage does nothing + sN.solve(zero_init_guess=True, picard=1) + assert sN._picard_stages == dict(picard_iterations=1), sN._picard_stages + + sP, vP = _driven_box("p", False) + sP.petsc_options["snes_max_it"] = 1 # one frozen-tangent iteration + sP.solve(zero_init_guess=True) + assert sP.solve_report.nl_its == 1 + + a, b = np.array(vN.data), np.array(vP.data) + rel = np.linalg.norm(a - b) / np.linalg.norm(b) + assert rel < 1.0e-8, ( + f"the Newton-path Picard step differs from a frozen-tangent step by {rel:.3e}") + + +FNORM = ("CONVERGED_FNORM_ABS", "CONVERGED_FNORM_RELATIVE") + + +def test_newton_warmup_takes_n_steps_then_converges_to_the_same_answer(): + """picard=N gives exactly N Picard steps, then Newton converges on the RESIDUAL, to the + same solution as a different warm-up count. ⚠️ Require FNORM convergence: without a + warm-up, pure Newton on this boundary-driven box stops on step stagnation — + CONVERGED_SNORM_RELATIVE at ||F|| = 4.5e-5 whatever the tolerance — which PETSc reports + as "converged". That run is exactly what the warm-up exists to prevent, and it must not + be used as a reference.""" + runs = {} + for n in (1, 3): + s, v = _driven_box("r%d" % n, True) + s.solve(zero_init_guess=True, picard=n) + assert s._picard_stages == dict(picard_iterations=n), s._picard_stages + assert str(s.solve_report.reason_str) in FNORM, ( + f"picard={n}: {s.solve_report.reason_str} at ||F|| = {s.solve_report.fnorm:.2e}") + runs[n] = np.array(v.data, copy=True) + assert np.linalg.norm(runs[3] - runs[1]) / np.linalg.norm(runs[1]) < 1.0e-6 + + +def test_pure_newton_without_warmup_stagnates_on_a_driven_box(): + """Why the automatic warm-up exists, pinned: switch it off (picard=-1) and pure Newton + on the boundary-driven box stalls on step size instead of converging on the residual. + If this ever starts converging on FNORM, the fixture no longer demonstrates the need + and the warm-up tests above lose their contrast — re-tune tau_y.""" + s, _ = _driven_box("nw", True) + s.solve(zero_init_guess=True, picard=-1) + assert s._picard_stages is None + assert str(s.solve_report.reason_str) not in FNORM, ( + f"pure Newton converged on the residual ({s.solve_report.reason_str}); the fixture " + "no longer shows why the warm-up is needed") + + +def test_cold_newton_takes_one_automatic_picard_step_and_warm_takes_none(): + """Layer 1: a cold start under the consistent tangent takes ONE Picard step + automatically; a warm re-solve takes none.""" + s, v = _driven_box("a", True) + s.solve(zero_init_guess=True) + assert s._picard_stages == dict(picard_iterations=1), s._picard_stages + assert str(s.solve_report.reason_str).startswith("CONVERGED") + s.solve(zero_init_guess=False) + assert s._picard_stages is None, "a warm re-solve took a Picard warm-up" + + +def test_continuation_alpha_is_linear_in_log_residual(): + """The ramp itself: alpha = 0 at F = start * F_ref, 1 at F = newton_at * F_ref, linear in + log F between, clipped outside.""" + s, _, _ = _yielding_box("ra", "continuation") + a0, a1 = s.continuation_alpha_start, s.continuation_newton_at + assert (a0, a1) == (1.0e-1, 5.0e-4) + assert s._continuation_alpha(1.0, 1.0) == 0.0 + assert s._continuation_alpha(a0, 1.0) == 0.0 + assert abs(s._continuation_alpha(a1, 1.0) - 1.0) < 1.0e-12 + assert abs(s._continuation_alpha(np.sqrt(a0 * a1), 1.0) - 0.5) < 1.0e-12 + assert s._continuation_alpha(1.0e-9, 1.0) == 1.0 + assert s._continuation_alpha(0.0, 1.0) == 1.0 + + +def test_continuation_warm_start_enters_the_ramp_at_its_residual(): + """#791: alpha is keyed on F_ref = ||F(u=0)||, not on the solve's own starting residual. + Keyed on its own F0, a warm start re-ran the whole Picard phase (measured on the notch: + 763 vs 10 iterations). So a warm start near the solution must take NO frozen-tangent + iteration, converge to the cold answer, and a converged state must take none at all.""" + s, v, _ = _yielding_box("wc", "continuation") + s.solve(zero_init_guess=True) + cold = s.solve_report + ref = np.array(v.data, copy=True) + assert str(cold.reason_str).startswith("CONVERGED"), cold.reason_str + st = s._continuation_stages + assert st["frozen_iterations"] >= 1 and st["final_alpha"] == 1.0, st + F_ref = st["F_ref"] + + s.solve(zero_init_guess=False) # re-solve the converged state + assert s.solve_report.nl_its == 0, s.solve_report.nl_its + + with uw.synchronised_array_update(): + v.data[...] = 1.02 * ref # 2% off the solution + s.solve(zero_init_guess=False) + st = s._continuation_stages + assert st["F_ref"] == pytest.approx(F_ref, rel=1.0e-12), "F_ref must not depend on the start" + assert st["frozen_iterations"] == 0, ( + f"a warm start 2% from the solution took {st['frozen_iterations']} frozen-tangent " + "iterations: the ramp is keyed on the solve's own residual again") + assert str(s.solve_report.reason_str).startswith("CONVERGED"), s.solve_report.reason_str + assert s.solve_report.nl_its < cold.nl_its + assert np.linalg.norm(np.array(v.data) - ref) / np.linalg.norm(ref) < 1.0e-5 + + +def test_continuation_restores_line_search_and_user_update_hook(): + """The ramp borrows the SNES update hook and the line search for its own solve. Both must + be handed back: user callbacks still fire every iteration inside the ramp, and the line + search is the configured one again afterwards.""" + s, _, _ = _yielding_box("hk", "continuation") + calls = [] + s.add_update_callback(lambda solver, it: calls.append(it)) + s.solve(zero_init_guess=True) + n = s.solve_report.nl_its + assert str(s.solve_report.reason_str).startswith("CONVERGED") + inside = [c for c in calls if c >= 0] + assert len(inside) == n, f"user callback fired {len(inside)} times in {n} iterations" + # nothing configures the line search on this fixture, so it is PETSc's newtonls default + assert s.snes.getLineSearch().getType() == "bt", ( + f"line search after the ramp is {s.snes.getLineSearch().getType()!r}: the " + "continuation solve's l2 leaked out") + + +def test_solve_report_config_is_what_the_solve_ran_with(): + """#806: the report carries the configuration IN FORCE during the solve, snapshotted at + the call. The continuation ramp switches the line search to l2 and the atol to + tolerance * F_ref for its own solve and restores both after, so the report must say l2 + while the live SNES says bt again. And the snapshot must not change the solver: reading + the nonlinear preconditioner with SNESGetNPC creates one (measured: DIVERGED_INNER).""" + s, _, _ = _yielding_box("cf", "continuation") + s.solve(zero_init_guess=True) + r = s.solve_report + assert str(r.reason_str).startswith("CONVERGED"), r.reason_str + cfg = r.config + assert cfg["linesearch"]["type"] == "l2", cfg["linesearch"] + assert s.snes.getLineSearch().getType() == "bt" + F_ref = s._continuation_stages["F_ref"] + assert cfg["snes"]["atol"] == pytest.approx(s.tolerance * F_ref, rel=1.0e-12) + assert cfg["tangent"]["consistent_jacobian"] == "continuation" + assert cfg["ksp"]["pc"] == "fieldsplit" and len(cfg["ksp"]["splits"]) == 2 + assert "npc" not in cfg["snes"] and not s.snes.hasNPC(), ( + "the config snapshot attached a nonlinear preconditioner") + + s2, _, _ = _yielding_box("cf2", False) + s2.solve(zero_init_guess=True) + assert s2.solve_report.config["tangent"]["consistent_jacobian"] is False + assert s2.solve_report.config["linesearch"]["type"] == "bt" + + +def _pressure_yield_box(tag, pressure_coupling): + """Lid-driven, open sides (no pressure nullspace), yield stress that depends on PRESSURE + so the consistent tangent has a d(eta)/dp term in the velocity-pressure block.""" + uw.reset_default_model() + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.25) + v = uw.discretisation.MeshVariable("Vpy" + tag, mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("Ppy" + tag, mesh, 1, degree=1, continuous=True) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + s.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel + s.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + # positive for every p, with a non-zero derivative in p (a smoothed |p|) + s.constitutive_model.Parameters.yield_stress = 0.3 + 0.3 * sympy.sqrt(p.sym[0] ** 2 + 1.0e-2) + s.add_essential_bc((1.0, 0.0), "Top") + s.add_essential_bc((0.0, 0.0), "Bottom") + s.petsc_options.delValue("ksp_monitor") + s.consistent_jacobian = True + s.newton_pressure_coupling = pressure_coupling + s.tolerance = 1.0e-6 + s.petsc_options["snes_max_it"] = 300 + return s, v, p + + +def _blocks(s): + """(A, J_vp, J_pv) of the assembled Jacobian at the solver's current state.""" + snes = s.snes + J, P, _ = snes.getJacobian() + snes.computeJacobian(snes.getSolution(), J, P) + keys = list(s._subdict) + isv, isp = s._subdict[keys[0]][0], s._subdict[keys[1]][0] + return J.createSubMatrix(isv, isv), J.createSubMatrix(isv, isp), J.createSubMatrix(isp, isv) + + +def test_partial_newton_keeps_the_velocity_block_and_restores_the_pressure_column(): + """newton_pressure_coupling=False: uu block identical to Newton's, J_vp = -B^T exactly + (the Picard column), where the full Newton tangent's J_vp is NOT -B^T.""" + sN, vN, pN = _pressure_yield_box("N", True) + sP, vP, pP = _pressure_yield_box("P", False) + # the same state in both: one Newton step from rest, copied across + sN.petsc_options["snes_max_it"] = 1 + sN.solve(zero_init_guess=True, picard=-1) + with uw.synchronised_array_update(): + vP.data[...] = vN.data; pP.data[...] = pN.data + sP.petsc_options["snes_max_it"] = 0 + sP.solve(zero_init_guess=False, picard=-1) + A_N, Jvp_N, Jpv_N = _blocks(sN) + A_P, Jvp_P, Jpv_P = _blocks(sP) + + dA = A_N.copy(); dA.axpy(-1.0, A_P) + assert dA.norm() / A_N.norm() < 1.0e-12, "partial Newton changed the velocity block" + BT = Jpv_P.copy(); BT.transpose() + dP = Jvp_P.copy(); dP.axpy(1.0, BT) + assert dP.norm() / Jpv_P.norm() < 1.0e-12, "partial Newton's pressure column is not -B^T" + BT_N = Jpv_N.copy(); BT_N.transpose() + dN = Jvp_N.copy(); dN.axpy(1.0, BT_N) + assert dN.norm() / Jpv_N.norm() > 1.0e-3, \ + "fixture has no pressure coupling in the Newton tangent; the test proves nothing" + assert sP.solve_report.config["tangent"]["newton_pressure_coupling"] is False + assert sN.solve_report.config["tangent"]["newton_pressure_coupling"] is True + + +def test_partial_newton_converges_to_the_frozen_tangent_answer(): + """The residual is untouched, so the partial tangent reaches the same solution as the + frozen (Picard) one.""" + sF, vF, pF = _pressure_yield_box("F", True) + sF.consistent_jacobian = False + sF.solve(zero_init_guess=True) + assert sF.solve_report.converged, str(sF.solve_report) + sP, vP, pP = _pressure_yield_box("Q", False) + sP.solve(zero_init_guess=True) + assert sP.solve_report.converged, str(sP.solve_report) + dv = np.abs(np.asarray(vP.data) - np.asarray(vF.data)).max() + assert dv < 1.0e-4 * np.abs(np.asarray(vF.data)).max(), dv + + +# ---- the warm-up inside the solve's contracts (review of PR #794) --------------------- + +def test_the_newton_stage_is_anchored_to_the_cold_residual(): + """After the warm-up, Newton stops where an uninterrupted solve would: at + tolerance x ||F(u_cold)||, not a further tolerance below the post-warm-up residual. + Without the anchor a secretly-linear yielding box (tau_y far above any stress), + which converges in 2 Newton iterations without the warm-up, stagnated on step size + (SNORM_RELATIVE) after it. The report counts the warm-up and measures the + reduction from the cold residual.""" + s, _ = _driven_box("anc", True, tau_y=1.0e4) + s.solve(zero_init_guess=True) + assert s._picard_stages == dict(picard_iterations=1), s._picard_stages + r = s.solve_report + assert str(r.reason_str) in FNORM, r.reason_str + assert r.config["tangent"]["picard_warmup_its"] == 1 + # The one Picard step solves this secretly-linear problem, so the anchored Newton + # stage has nothing to do: the solve is that ONE step, counted in the report (the + # report used to show only the Newton solve's iterations). Before the anchor, Newton + # took 5 more iterations and stopped on SNORM. + assert r.nl_its == 1, r.nl_its + assert r.reduction <= s.tolerance * (1.0 + 1.0e-9), r.reduction + + +def test_force_setup_is_not_read_as_picard(): + """`picard` follows every other argument of the Scalar/Vector/MultiComponent solve(): + five subclasses pass `_force_setup` positionally, and with `picard` second it landed + there — a forced Picard warm-up and a blended-kernel compile on a LINEAR + advection-diffusion solve.""" + uw.reset_default_model() + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(-1.0, -1.0), maxCoords=(1.0, 1.0), cellSize=0.5) + x, y = mesh.X + T = uw.discretisation.MeshVariable("T_fs", mesh, 1, degree=2) + adv = uw.systems.AdvDiffusion(mesh, T, sympy.Matrix([[-y, x]])) + for b in ("Left", "Right", "Top", "Bottom"): + adv.add_dirichlet_bc(0.0, b) + adv.solve(timestep=0.01, _force_setup=True) + assert adv._picard_stages is None, adv._picard_stages + assert adv._picard_blend is False + + +def test_changing_the_tangent_recompiles_it(): + """The tangent is compiled into the Jacobian, so changing consistent_jacobian or + newton_pressure_coupling after a solve marks a rewire; without it the next solve ran + the old tangent while the report recorded the new one. alpha is per solver: a shared + container let one solver's warm-up move another solver's tangent.""" + s, _ = _driven_box("tg", True) + s.solve(zero_init_guess=True) + assert s._needs_function_rewire is False + s.consistent_jacobian = False + assert s._needs_function_rewire is True + s.solve(zero_init_guess=False) + assert s._needs_function_rewire is False + s.newton_pressure_coupling = False + assert s._needs_function_rewire is True + s2, _ = _driven_box("tg2", True) + assert s._get_newton_alpha() is not s2._get_newton_alpha() + + +def test_a_flux_made_nonlinear_after_a_solve_takes_the_cold_warmup(): + """Linear viscosity first (no warm-up: the tangents coincide), then a viscosity that + depends on the velocity, solved warm and then cold: the cold solve must take the + automatic Picard step. The cached linearity verdict survived the rebuild and the + step was skipped on a nonlinear model.""" + uw.reset_default_model() + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.25) + v = uw.discretisation.MeshVariable("Vfl", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("Pfl", mesh, 1, degree=1, continuous=True) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + s.constitutive_model = uw.constitutive_models.ViscousFlowModel + s.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + s.add_essential_bc((1.0, 0.0), "Top") + s.add_essential_bc((0.0, 0.0), "Bottom") + s.petsc_use_pressure_nullspace = True + s.petsc_options.delValue("ksp_monitor") + s.tolerance = 1.0e-6 + s.solve(zero_init_guess=True) + assert s._picard_stages is None + s.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + 5.0 * v.sym.dot(v.sym) + s.solve(zero_init_guess=False) + assert s._picard_stages is None # warm: no warm-up by design + s.solve(zero_init_guess=True) + assert s._picard_stages == dict(picard_iterations=1), s._picard_stages diff --git a/tests/test_1069_pressure_gauge_is_recorded.py b/tests/test_1069_pressure_gauge_is_recorded.py new file mode 100644 index 000000000..989da001d --- /dev/null +++ b/tests/test_1069_pressure_gauge_is_recorded.py @@ -0,0 +1,182 @@ +"""The pressure gauge a Stokes solve runs under is recorded, and a gauge-fixed pressure +feeding a pressure-dependent viscosity is warned about. + +Why: with ``petsc_use_pressure_nullspace=True`` the projection imposes mean(p) = 0 on every +Krylov solve. A Drucker-Prager yield stress ``C + sin(phi) p`` then reads a gauge choice, not +a physical pressure, and the choice moves with the yielded region between iterations. On an +open domain (free surface, traction boundary) there is no nullspace and the flag must be +off; on a closed box the physical level has to be supplied separately. Either way the +transcript must say which it was, without a live process to introspect. +""" +import warnings + +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + + +def _box(tag, closed, pressure_dependent, nullspace): + uw.reset_default_model() + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.25) + v = uw.discretisation.MeshVariable("Vg" + tag, mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("Pg" + tag, mesh, 1, degree=1, continuous=True) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + s.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel + s.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + if pressure_dependent: + s.constitutive_model.Parameters.yield_stress = 0.3 + 0.3 * sympy.sqrt(p.sym[0] ** 2 + 1.0e-2) + # a stated regularisation, so the only warning these tests can see is the gauge one + s.constitutive_model.Parameters.plastic_rate_strengthening = 0.05 + else: + s.constitutive_model.Parameters.yield_stress = 0.3 + s.add_essential_bc((1.0, 0.0), "Top") + s.add_essential_bc((0.0, 0.0), "Bottom") + if closed: + s.add_essential_bc((0.0, sympy.oo), "Left") + s.add_essential_bc((0.0, sympy.oo), "Right") + s.petsc_use_pressure_nullspace = nullspace + s.petsc_options.delValue("ksp_monitor") + s.tolerance = 1.0e-5 + s.petsc_options["snes_max_it"] = 60 + return s + + +def test_gauge_facts_are_in_the_solve_config_open_domain(): + """Open sides: no nullspace attached, the rheology's pressure dependence is recorded.""" + s = _box("open", closed=False, pressure_dependent=True, nullspace=False) + with warnings.catch_warnings(): + warnings.simplefilter("error", RuntimeWarning) + try: + s.solve() + except RuntimeWarning as w: # only the gauge warning is a failure here + assert "GAUGE" not in str(w), str(w) + g = s.solve_report.config["gauge"] + assert g["pressure_nullspace_requested"] is False + assert g["pressure_nullspace_attached"] is False + assert g["pressure_in_rheology"] is True + assert g["pressure_dirichlet_bcs"] == [] + assert s.solve_report.config["ksp"]["restart"] == 100 + + +def test_closed_box_with_pressure_dependent_viscosity_and_nullspace_warns(): + """Closed box + nullspace + p-dependent yield: the projection fixes the gauge the yield + law reads. Warned, and recorded as attached.""" + s = _box("closed", closed=True, pressure_dependent=True, nullspace=True) + with pytest.warns(RuntimeWarning, match="GAUGE-fixed pressure"): + s.solve() + g = s.solve_report.config["gauge"] + assert g["pressure_nullspace_requested"] is True + assert g["pressure_nullspace_attached"] is True + assert g["pressure_in_rheology"] is True + + +def test_closed_box_with_constant_yield_does_not_warn(): + s = _box("const", closed=True, pressure_dependent=False, nullspace=True) + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + s.solve() + assert not [w for w in caught if "GAUGE" in str(w.message)] + g = s.solve_report.config["gauge"] + assert g["pressure_in_rheology"] is False + assert g["pressure_nullspace_attached"] is True + + +def test_transcript_event_carries_the_gauge_and_the_tangent(): + """The step transcript's solve event records the gauge and tangent facts, and the gauge + warning itself lands in the step as a warning.""" + s = _box("tr", closed=True, pressure_dependent=True, nullspace=True) + model = uw.get_default_model() + with warnings.catch_warnings(): + warnings.simplefilter("always") + with model.step(0.1, label="gauge step"): + s.solve() + entry = model.transcript[-1] + solves = [e for e in entry.events if e["kind"] == "solve"] + assert solves, [e["kind"] for e in entry.events] + ev = solves[-1] + assert ev["gauge"]["pressure_nullspace_attached"] is True + assert ev["gauge"]["pressure_in_rheology"] is True + assert "consistent_jacobian" in ev["tangent"] + assert "newton_pressure_coupling" in ev["tangent"] + recorded = [e for e in entry.events + if e["kind"] == "warning" and "GAUGE" in str(e.get("message", ""))] + assert recorded, "the gauge warning was not recorded in the step" + + +def test_krylov_restart_default_is_in_force_and_leaves_a_user_value_alone(): + """Restart 100 is set at construction for the outer KSP (every class) and the Stokes + velocity split. It was first written inside the `tolerance` / `strategy` setters, so + a Stokes solver kept FGMRES(30) unless `strategy` was assigned, and a `tolerance` + assignment overwrote a user's own restart.""" + from petsc4py import PETSc + opts = PETSc.Options() + s = _box("rs", closed=False, pressure_dependent=False, nullspace=False) + pfx = s.petsc_options_prefix + assert opts.getInt(pfx + "ksp_gmres_restart") == 100 + assert opts.getInt(pfx + "fieldsplit_velocity_ksp_gmres_restart") == 100 + s.petsc_options["ksp_gmres_restart"] = 30 + s.tolerance = 1.0e-7 + assert opts.getInt(pfx + "ksp_gmres_restart") == 30 + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.5) + T = uw.discretisation.MeshVariable("Trs", mesh, 1, degree=1) + pois = uw.systems.Poisson(mesh, u_Field=T) + assert opts.getInt(pois.petsc_options_prefix + "ksp_gmres_restart") == 100 + + +def _live_restart(ksp, tmp_path, tag): + """The restart the LIVE KSP runs with (petsc4py has no getter): from its view.""" + import re + from petsc4py import PETSc + f = tmp_path / f"ksp_{tag}.txt" + viewer = PETSc.Viewer().createASCII(str(f)) + ksp.view(viewer) + viewer.destroy() + m = re.search(r"restart=(\d+)", f.read_text()) + return int(m.group(1)) if m else None + + +def test_krylov_restart_reaches_the_live_solver(tmp_path): + """The option database is not the solver: a KSPSetType after setFromOptions (the + single-field custom-P FMG route) re-created the Krylov context at PETSc's 30 while + the database and the report said 100. Checked on the live objects after a solve: + default Stokes (outer and velocity split) and an FMG Poisson.""" + s = _box("lv", closed=False, pressure_dependent=False, nullspace=False) + s.solve() + assert _live_restart(s.snes.getKSP(), tmp_path, "stokes") == 100 + vel = s.snes.getKSP().getPC().getFieldSplitSubKSP()[0] + assert _live_restart(vel, tmp_path, "vel") == 100 + uw.reset_default_model() + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.5, refinement=2) + T = uw.discretisation.MeshVariable("Tlv", mesh, 1, degree=1) + pois = uw.systems.Poisson(mesh, u_Field=T) + pois.constitutive_model = uw.constitutive_models.DiffusionModel + pois.constitutive_model.Parameters.diffusivity = 1.0 + pois.f = 1.0 + pois.add_dirichlet_bc(0.0, "Bottom") + pois.preconditioner = "fmg" + pois.solve() + assert _live_restart(pois.snes.getKSP(), tmp_path, "pois_fmg") == 100 + + +def test_switching_on_a_pressure_dependent_yield_is_seen_by_the_record(): + """Solve with a constant yield stress, then switch on Drucker-Prager and solve WARM: + the record and the gauge warning must see the new rheology. The cached symbolic + facts survived the rebuild (a warm solve's plan never reads them, and the read after + the build hit the old dict): measured pressure_in_rheology=False and no warning, + while the kernel had changed.""" + s = _box("sw", closed=True, pressure_dependent=False, nullspace=True) + s.solve() + assert s.solve_report.config["gauge"]["pressure_in_rheology"] is False + cm = s.constitutive_model + cm.Parameters.yield_stress = 0.3 + 0.3 * sympy.sqrt(s.p.sym[0] ** 2 + 1.0e-2) + cm.Parameters.plastic_rate_strengthening = 0.05 + with pytest.warns(RuntimeWarning, match="GAUGE-fixed pressure"): + s.solve(zero_init_guess=False) + assert s.solve_report.config["gauge"]["pressure_in_rheology"] is True diff --git a/tests/test_1071_plastic_rate_strengthening.py b/tests/test_1071_plastic_rate_strengthening.py new file mode 100644 index 000000000..2d02eb79e --- /dev/null +++ b/tests/test_1071_plastic_rate_strengthening.py @@ -0,0 +1,133 @@ +"""The declared plastic rate strengthening eta_reg: eta_pl -> tau_y/(2 edot_II) + m eta_reg. + +Why: a yielded element's consistent tangent has no modulus along its own strain-rate +direction, so a yielded layer is a mechanism the velocity multigrid cannot solve +(Spiegelman notch, 2026-10-03: Rayleigh quotient 1.8e-7). A static viscosity floor clips +only the fastest cells; eta_reg floors the plastic viscosity AND the tangent everywhere. +Its scale is the model's (no universal dimensionless value), it is laddered as a +multiplier m from the viscous limit down to 1, and the solve records it. +""" +import warnings + +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + + +def _box(tag, eta_reg): + uw.reset_default_model() + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.25) + v = uw.discretisation.MeshVariable("Vr" + tag, mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("Pr" + 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 = 0.3 + 0.3 * sympy.sqrt(p.sym[0] ** 2 + 1.0e-2) + if eta_reg is not None: # None = never stated; 0.0 = stated as zero + cm.Parameters.plastic_rate_strengthening = eta_reg + s.add_essential_bc((1.0, 0.0), "Top") + s.add_essential_bc((0.0, 0.0), "Bottom") + s.petsc_options.delValue("ksp_monitor") + s.tolerance = 1.0e-6 + s.petsc_options["snes_max_it"] = 300 + return s, v, p, cm + + +def test_the_term_floors_the_plastic_viscosity_and_is_recorded(): + """At a strain rate far above yield the effective viscosity sits at m*eta_reg (the + plastic branch tau_y/2edot -> 0), and the solve config says what was stated.""" + s, v, p, cm = _box("a", 0.05) + from underworld3.function.expressions import unwrap_expression + eta = unwrap_expression(cm.viscosity, mode="symbolic_keep_constants") + # substitute a huge strain rate and p = 0 by evaluating the unwrapped symbolic law + L = s.Unknowns.L + subs = {L[i, j]: (1.0e6 if (i, j) == (0, 1) else 0.0) for i in range(2) for j in range(2)} + subs[p.sym[0]] = 0.0 + val = float(unwrap_expression(sympy.sympify(eta).subs(subs), mode="nondimensional")) + assert abs(val - 0.05) < 2.0e-3, val # eta_reg, with the soft-min's rounding + with warnings.catch_warnings(): + warnings.simplefilter("error", RuntimeWarning) + try: + s.solve() + except RuntimeWarning as w: + assert "plastic_rate_strengthening" not in str(w), str(w) + reg = s.solve_report.config["regularisation"] + assert float(str(reg["plastic_rate_strengthening"]).split()[0]) == 0.05, reg + assert reg["rate_strengthening_scale"] == 1.0 + + +@pytest.mark.parametrize("stated", [None, 0, 0.0], ids=["unstated", "int0", "float0"]) +def test_pressure_dependent_yield_without_the_term_warns(stated): + """Unstated and stated-as-zero are the same model: no term, and the warning. (0.0 was + missed once: sympy.Float(0.0) == 0 is False, and `!= 0` compiled a zero term.)""" + s, v, p, cm = _box("b" + str(stated), stated) + assert cm._rate_strengthening_control() is None + with pytest.warns(RuntimeWarning, match="no stated plastic_rate_strengthening"): + s.solve() + recorded = str(s.solve_report.config["regularisation"]["plastic_rate_strengthening"]) + assert float(recorded.split()[0]) == 0.0, recorded # "0 [pascal * second]" unstated + + +def test_the_ladder_ends_at_the_stated_problem(): + """solve(homotopy="rate_strengthening") finds the viscous limit from the fields, + ladders m down to exactly 1, and the result is the stated problem's solution.""" + s, v, p, cm = _box("c", 0.05) + report = s.solve(homotopy="rate_strengthening", homotopy_options=dict(verbose=False)) + assert report["reached_one"] and report["converged"], report + rungs = report["rungs"] + # the entry: 1e3 x max(eta_ve / eta_reg) = 1e3 x 1 / 0.05, at which nothing yields + assert rungs[0][0] == pytest.approx(2.0e4, rel=1.0e-9), rungs + # the ladder RAN: the viscous limit read from the fields is well above 1, and the + # rungs descend from it to exactly 1 (a skipped ladder is [entry, 1.0]) + assert report["scale0"] > 5.0, report + assert len(rungs) >= 4, rungs + ms = [m for m, reason, its in rungs[1:]] + assert all(a > b for a, b in zip(ms, ms[1:])), ms + assert rungs[-1][0] == 1.0 + assert cm.rate_strengthening_scale == 1.0 + v_ladder = np.asarray(v.data).copy() + s2, v2, p2, cm2 = _box("d", 0.05) + s2.solve() + assert s2.solve_report.converged + dv = np.abs(v_ladder - np.asarray(v2.data)).max() + assert dv < 1.0e-4 * np.abs(np.asarray(v2.data)).max(), dv + + +def test_unknown_homotopy_value_is_refused(): + """A misspelt homotopy string is truthy; it must not silently run the delta march.""" + s, v, p, cm = _box("e", 0.05) + with pytest.raises(ValueError, match="rate_strengthening"): + s.solve(homotopy="rate-strengthening") + + +def test_the_term_is_in_the_viscoelastic_plastic_law(): + """VEP carries the same declared term; the ladder refuses VEP (it would advance the + stress history on every rung) with NotImplementedError, not an AttributeError.""" + uw.reset_default_model() + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.25) + v = uw.discretisation.MeshVariable("Vve", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("Pve", mesh, 1, degree=1, continuous=True) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + s.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel(s.Unknowns, order=1) + cm = s.constitutive_model + cm.Parameters.shear_viscosity_0 = 1 + cm.Parameters.shear_modulus = 1 + cm.Parameters.dt_elastic = sympy.sympify(1) / 10 + cm.Parameters.yield_stress = 0.3 + from underworld3.function.expressions import unwrap_expression + def _has_m(expr): + e = unwrap_expression(expr, mode="symbolic_keep_constants") + m = cm._get_rate_strengthening_scale() + return m in sympy.sympify(e).atoms(sympy.Symbol) + assert not _has_m(cm.viscosity) + cm.Parameters.plastic_rate_strengthening = 0.05 + assert _has_m(cm.viscosity) + with pytest.raises(NotImplementedError, match="stress history"): + s.solve(homotopy="rate_strengthening", homotopy_options=dict(verbose=False))