diff --git a/.claude/skills/nonlinear-solver/SKILL.md b/.claude/skills/nonlinear-solver/SKILL.md new file mode 100644 index 000000000..7bd3965a0 --- /dev/null +++ b/.claude/skills/nonlinear-solver/SKILL.md @@ -0,0 +1,186 @@ +--- +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`. +--- + +# nonlinear-solver + +The recipe that gets a **hard viscoplastic (Drucker–Prager) Stokes** problem to +converge, and — more importantly — the list of setup mistakes that stop it. The +central lesson from the Spiegelman hard-case study (`η_bg=1e26`, `V=10`): every +failure was a **solver-configuration** error, not a bad Jacobian. If the correct +setup is a minefield for an expert, that is an API regression — so the goal is to +make the correct path the default path. + +Design of record: `docs/developer/design/nonlinear-solver-homotopy-warmstart.md`. +Yield-law maths, tangent-per-model, quadratic-convergence check: `plasticity-solvers`. + +--- + +## 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). + 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). + +2. **Multi-solve δ-continuation** (NOT an in-solve ramp). Hold δ **constant** for a + full nonlinear solve to tolerance; warm-start the next, smaller δ from that + converged state; march δ down to the sharp surface. δ is a `constants[]` atom, + so each step is a recompile-free `PetscDSSetConstants` update. + + **This is now one call** — the model advertises the homotopy and `solve()` marches it: + + ```python + stokes.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel + cm.Parameters.shear_viscosity_0 = ... + cm.Parameters.yield_stress = ... # a plain pressure-dependent yield + report = stokes.solve(homotopy=True) # smooth mode, tangent, march: automatic + report["settled_delta"] # smallest δ that converged + ``` + + `solve(homotopy=True)` sets the smooth mode (softmin + power-mean), picks the + tangent the model asks for (Newton for DP, Picard for elastic VEP), and runs a + residual-guided march that accelerates on easy steps and reverts + retries a + failed δ more gently. Tune with + `homotopy_options=dict(delta0=…, down=…, dmin=…, entry_maxit=…, step_maxit=…)`; + the driver is also callable directly as + `underworld3.systems.yield_continuation`. + +3. **Consistent-Newton tangent** for non-elastic DP (`consistent_jacobian=True`); + **Picard** for elastic VEP — see `plasticity-solvers` for the per-model table. + +4. **`bt` line search** with the consistent tangent on a smooth (δ>0) surface. + +The δ-march is cheap: with the power-mean smoother a converged δ warm-starts every +sharper δ in ≈0 Newton iterations, so a residual-guided auto-descent costs almost +nothing. + +--- + +## DO NOT ramp δ inside one SNES solve + +Ramping δ **inside a single SNES solve** (a `SNESSetUpdate` callback that sharpens +the yield surface between Newton iterations) is **proven dead** — it diverges +`DIVERGED_LINEAR_SOLVE` after ~2 iterations and grinds for ~2 hours, **even on the +proven solver config**. Mechanism: the continuation only sharpens δ from a +**converged**, well-conditioned iterate; the in-SNES ramp sharpens δ **mid-solve** +at a far-from-solution iterate where the consistent-Newton Jacobian on a sharpening +surface is ill-conditioned and the linear solve fails. **Hide the *continuation*, +not the *ramp*.** (This supersedes the `enable_yield_homotopy()` in-SNES ramp still +described in `plasticity-solvers`; that path is retained only as a dead-experiment +record — use the multi-solve continuation above.) + +--- + +## CONFIG TRAP LIST + +Each of these produces a *different* failure a few steps in — that is why the hard +case felt like whack-a-mole. Check them first. + +| Trap | Symptom | Fix | +|---|---|---| +| Consistent Newton makes the velocity block **non-symmetric**; a Chebyshev/Richardson MG smoother assumes SPD | smoother diverges / stalls → `DIVERGED_LINEAR_SOLVE` or an endless grind | **Now the default** — the FMG bundle ships `mg_levels_ksp_type=gmres` + `pc_type=sor` + `norm_type=none`. Only an issue if you override it, or on GAMG (which uses PETSc's chebyshev default) | +| `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 | +| Hand-rolled `snes_monitor` to "see what's happening" | you read residuals but miss the tell | use `solve_with_diagnostics` / `get_snes_diagnostics` instead (below) | + +**The diagnostic tell:** `solver.get_snes_diagnostics()["linear_iterations"] ≈ 1` +per Newton step means the linear solve is doing **no real work** (the FMG-1-iteration +trap). A healthy consistent-Newton solve does real Krylov work each step and +converges quadratically. Use `solve_with_diagnostics()`, not a hand-rolled monitor. + +--- + +## Automatic warm-start (Layer 1 — landed) + +`solver.has_solution` is a **public, read-only** status flag: `True` only after a +solve whose SNES converged; reset on a structural rebuild (remesh / adapt / +mesh-mover — the `is_setup=False` hook); kept through coefficient changes (viscosity, +δ, 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. + +```python +stokes.consistent_jacobian = True +stokes.solve() # cold → one automatic Picard step, then Newton +if stokes.has_solution: + ... +``` + +--- + +## 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 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. +- **Layer 3 — DONE:** the FMG velocity smoother defaults to `gmres`+`sor` with + `mg_levels_ksp_norm_type=none` (fixed-cost V-cycle), unconditionally — see + "Multigrid depth" below. +- **Layer 2 — DONE:** the model advertises the homotopy + (`supports_yield_homotopy` / `_yield_homotopy_control`) and + `stokes.solve(homotopy=True, homotopy_options=...)` runs the residual-guided + continuation, returning the march summary. + +--- + +## Multigrid depth — how to measure a smoother honestly + +**A two-level hierarchy is a coarse-grid correction, not a V-cycle.** Smoother +comparisons made on one are misleading: the gmres-over-richardson margin measured on +the Spiegelman notch is only 5 % at 3 levels but **25 % at 4** (ρ per V-cycle 0.746 → +0.560), because a deeper cycle applies the smoother on more coarse operators. Judge a +smoother at depth or not at all. + +To get depth without a monster problem, refine a **deliberately ultra-coarse NESTED +base**: `make_notch_mesh.py 1` (492 cells) + uniform `refinement=N` gives 3 levels / +7,872 cells at `N=2` and 4 levels / 31,488 at `N=3` — deeper *and* smaller than the old +2-level 38,580-cell setup. In MG you want the coarsest grid as coarse as it can be +before the problem breaks down. + +- **Never use a non-nested hierarchy** here — it does not give strong MG convergence + (maintainer ruling). Uniform refinement nests by construction. +- Accepted tradeoff: uniform refinement does **not** snap new boundary nodes back to + the analytic notch arcs (no CAD/EGADS model attached), so the corner geometry is + frozen at the coarse mesh's chords on every level. + +Measure with `fmg_contraction_probe.py` (ρ_MG per V-cycle; `<0.5` healthy, `0.8–0.95` +struggling, `≥0.98` hangs) or `smoother_depth_sweep.py` (pays the mesh build + viscous +seed once, sweeps smoothers in-process) in the Spiegelman study. + +**`solve_report` cannot see the smoother.** It records the *Newton* contraction; the +outer KSP is Eisenstat–Walker-collapsed to ~1 iteration/step, so the smoother's work +hides inside the velocity sub-block. Probe the `fieldsplit_velocity_` sub-KSP directly. + +**A smoother will not rescue small ξ.** At the hard corner the failure is operator +conditioning — the coarsest grid cannot represent the viscosity contrast — and at 4 +levels *every* smoother fails there (richardson outright, gmres with ρ>1). Use the δ/ξ +continuation to stay in the solvable region. + +## Gotchas + +- **`./uw build` → `amr-dev` env**; verify `uw.__file__` is the worktree site-packages. +- **Run VEP/consistent-Newton tests UNFORKED** — `pytest --forked` SIGABRTs (fork of + multithreaded PETSc). +- Benchmark **every** default change — "Solver Stability is Paramount". +- ξ (rate-strengthening) is a **non-homotopic** regularisation: put a user loop + *around* `solve()`, never inside the δ-march. + +## Reference + +- Design: `docs/developer/design/nonlinear-solver-homotopy-warmstart.md`, + `jacobian-consistent-tangent.md`, `solver-strategies-catalogue.md`. +- Continuation driver: `underworld3.systems.yield_continuation`. +- Diagnostics: `SNES_*.get_snes_diagnostics()` / `solve_with_diagnostics()`. +- Related skills: `plasticity-solvers` (yield law + tangent per model), + `free-surface-convection`, `adaptive-meshing`. diff --git a/docs/developer/design/nonlinear-solver-homotopy-warmstart.md b/docs/developer/design/nonlinear-solver-homotopy-warmstart.md new file mode 100644 index 000000000..0ce8dda93 --- /dev/null +++ b/docs/developer/design/nonlinear-solver-homotopy-warmstart.md @@ -0,0 +1,282 @@ +# Nonlinear solver: automatic warm-start and single-parameter yield homotopy + +**Status: IMPLEMENTED** (2026-07, branch `feature/nonlinear-warmstart-homotopy`). +All four stages have landed; each section below records what was built and, where +measurement contradicted the design, what was corrected. +(Design captured 2026-07, from the Spiegelman viscoplastic hardening work.) + +## Why this exists + +Hard viscoplastic (Drucker–Prager) Stokes problems do not converge from a cold +Newton start on the sharp yield surface. The working recipe is well understood — +a viscous warm start, then a **δ-continuation** (solve at a smooth yield law, then +warm-start progressively sharper solves), with the consistent-Newton tangent and a +multigrid smoother that tolerates the non-symmetric Newton operator. + +The problem is that *setting that recipe up correctly is a minefield*, and the +minefield is the real defect. Across a single working session the recipe was +mis-configured many times — default line search, a pressure-nullspace-singular LU, +`preconditioner="fmg"` giving a 1-iteration outer KSP, the default Chebyshev +smoother diverging on the non-symmetric operator, cold-start `ε̇=0 → τ_y/0 → NaN`. +Each produced a *different* failure a few steps in. **If the correct setup is that +fiddly for an expert, it is an API regression, not a user error.** The goal of this +design is to make the correct path the default path: the user writes their rheology +and calls `solve()`. + +## What was established (evidence) + +These results are on the Spiegelman notch at the genuinely hard regime +(`η_bg = 1e26 Pa·s`, `V = 10 mm/yr`, C = 1e8 Pa, φ = 30°), coarse mesh, and are +reproducible via `convergence.py` in the hard-case study. + +- **Multi-solve δ-continuation converges.** Constant δ per solve, solved to + tolerance, warm-started into the next (smaller) δ, consistent-Newton tangent, + `bt` line search, non-symmetry-safe smoother: settles at δ≈2e-3 with the plastic + band active. This is the trusted path. +- **In-SNES δ-ramp does not — proven, not masked.** Ramping δ *inside one SNES + solve* via a `SNESSetUpdate` callback (with the *same* proven config) diverges + with `DIVERGED_LINEAR_SOLVE` after 2 iterations and grinds for ~2 hours doing it. + The mechanism is fundamental: the continuation only sharpens δ from a **converged** + iterate, where the consistent-Newton Jacobian is well-conditioned; the in-SNES + ramp sharpens δ **mid-solve**, at a far-from-solution iterate, where the + consistent-Newton Jacobian on a sharpening yield surface is ill-conditioned and + the linear solve fails. **Do not re-attempt hiding δ inside a single Newton + solve.** Hide the *continuation*, not the ramp. +- **The δ-march is cheap.** With the power-mean smoother, once the first δ converges + the warm start already satisfies every sharper δ in ≈0 Newton iterations, so an + automatic residual-guided descent costs almost nothing. +- **The power-mean smoother is faithful.** At the settled δ it undershoots the true + `Min` yield by ≤0.2 %, and 0 % in genuinely yielded cells (see + {doc}`jacobian-consistent-tangent`). Earlier "it undershoots badly" readings were + an artefact of a diagnostic that evaluated a sharper power-mean, not the exact + `Min`. +- **The recursion prerequisite is fixed** (landed on the same branch as this design: + the DP model's `yield_stress_min` guard used the wrong sentinel and wrapped every + yield in `Max(<−∞ Parameter>, …)`, which recursed under the consistent tangent). + +The through-line: every failure was a **solver-configuration** error, upstream of +the Jacobian and unrelated to the recursion fix. The `get_snes_diagnostics()` field +`linear_iterations` is the tell — ≈1 linear iteration per Newton step means the +linear solve is doing no real work. + +## The design + +Three layers, from most general to most specific. Layer 1 stands alone and helps +every nonlinear solver; layers 2–3 build on it. + +### Layer 1 — automatic warm-start (all nonlinear solvers) + +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 +(`consistent_jacobian=False`) and a Picard→Newton α-blend. + +Rules: + +- **Cold** (no usable solution): take **one** Picard step, then Newton. +- **Warm** (a usable solution is present): straight Newton, **no** Picard step — a + Picard step from a good iterate takes a linear-convergence step off the good + quadratic starting point (unnecessary, mildly harmful). +- **Linear problem**: free either way — the frozen operator *is* the real operator, + so the single Picard step *is* the exact solve and Newton then reports converged in + 0 iterations. UW3 already probes the assembled Jacobian for solution-dependence, so + a linear problem can also simply skip the warm-up. Overhead ≈ one convergence check. + +**Cold-vs-warm is auto-detected**, and detection is safe by construction: guessing +*cold* when actually warm costs one Picard step that instantly re-converges; the only +harmful direction (warming off stale, unrelated field data) is what the opt-out flag +is for. + +- **`solver.has_solution`** — a **public** property, set `True` only after a solve + that **converged** (`reason > 0`), `False` initially, after a **diverged** solve, + and after a **structural rebuild** (mesh change / adaptivity / mesh-mover — the + existing `is_setup=False` invalidation is the natural hook). It is kept through + *coefficient* changes (new viscosity, new δ, new BCs) so continuation and + time-stepping warm-start correctly. It doubles as a user-facing status flag. +- Secondary signal for a hand-set initial guess that has not been solved yet: + `‖v‖ ≈ 0 ⇒ cold`. +- **Emergent auto-recovery**: a diverged solve leaves `has_solution=False`, so the + *next* `solve()` auto-cold-starts (Picard warm-up) rather than warming off garbage. + +**API change**: today `zero_init_guess` defaults to `True` (always cold). Make it a +tri-state whose **default auto-detects**; `zero_init_guess=True` means "force fresh / +discard any solution", `False` means "insist on warming". The onus flips from "flag +when you have a solution" to "flag when you want to throw one away". + +```{note} +This changes a core default and is "benchmark before flipping" territory. It is +mostly a bug-fix — today `solve()` in a loop silently re-zeros each iteration — but +must be validated. Explicit warm (`zero_init_guess=False`) is already the common +deliberate choice in the codebase; explicit cold is rare, so the flip aligns with +existing practice. See the open verification points. +``` + +### Layer 2 — advertised single-parameter homotopy + +A constitutive model **advertises** a single-parameter homotopy, exactly as it +advertises an approximate-Jacobian flux: + +- `cm.supports_yield_homotopy` → `True` on `ViscoPlasticFlowModel` / VEP, `False` on + plain viscous. `solve(homotopy=True)` on an unsupported model raises a clear error. +- `cm._yield_homotopy_control()` → puts the model in its smooth mode + (`yield_mode="softmin"`, `yield_smoother="powermean"`), returns the δ `constants[]` + atom and the recommended tangent (Newton for non-elastic DP, Picard for elastic + VEP). The model owns "what the homotopy means for me", so the mechanism generalises + to other models and applications. + +`solve(homotopy=True, homotopy_options=...)` then runs the **continuation** (never +the in-SNES ramp): + +1. Set δ = δ₀ **large**. At large δ the power-mean is the harmonic mean, which is + bounded by η_bg even as `ε̇ → 0`, so the cold Picard warm-up (Layer 1) is + well-conditioned — **this is what removes the need for a separate viscous + pre-solve, and why the cold `ε̇=0` NaN does not occur.** +2. Newton to tolerance at δ₀ (the tangent from `_yield_homotopy_control`). +3. **Residual-guided descent**: march δ down while solves keep converging (accelerate + when a step is easy, ease off when it is hard), warm-started each time, and settle + at the smallest feasible δ. This is the default because the descent is cheap and + strictly better than "assume homotopy is not needed → stall → retreat". +4. Report via `get_snes_diagnostics()`. + +`homotopy_options` (all defaulted) exposes only the march knobs — `delta0`, `down`, +`dmin`, `entry_maxit`, `step_maxit`, and the settle criteria. Fine-tuning is optional. + +**Scope: single parameter only.** A one-parameter homotopy is easy to generalise; +multi-parameter is not, and is unnecessary here. The rate-strengthening ξ term is a +*non-homotopic* regularisation, not a homotopy — it belongs in a user-domain loop +*around* `solve()`, not inside it. + +### Layer 3 — smoother as a consistent-Newton consequence + +Turning on the consistent tangent adds the `∂η/∂(grad v)` term, which makes the +velocity-block operator **non-symmetric**. The default multigrid smoothers +(Chebyshev / Richardson) assume a symmetric/SPD operator: Chebyshev can **diverge** +on a non-symmetric operator, Richardson stalls. The fix is one line — a +non-symmetry-safe smoother, `mg_levels_ksp_type = gmres` (with `mg_levels_pc_type = +sor`). + +This is a consequence of `consistent_jacobian=True`, **not** of homotopy — anyone +using the consistent tangent with FMG hits it. + +**LANDED, and the benchmark changed the shape of the fix.** The smoother is now +`gmres` + `sor` unconditionally in the FMG bundle, with `mg_levels_ksp_norm_type = +none` (a fixed-cost V-cycle: exactly `max_it` smoother iterations, no residual norm, +no early exit). Two findings drove that: + +- **The gmres margin grows with multigrid depth.** Measured on the Spiegelman notch + (Drucker–Prager, `η_bg=1e26`, `V=10`, `δ=0.01`) over a *nested* hierarchy built by + uniformly refining an ultra-coarse 492-cell base, per-V-cycle contraction ρ at + `ξ=0.01` and the same four smoother iterations: 3 levels — 0.722 richardson vs + 0.686 gmres (5 %); 4 levels — 0.746 vs **0.560** (25 %). A deeper cycle applies the + smoother on more coarse operators, each non-symmetric under the consistent tangent, + so smoother quality compounds. A measurement on a 2-level hierarchy shows almost no + effect and is **misleading** — a two-level cycle is a coarse-grid correction, not a + V-cycle. +- **Therefore no gate.** Since shallow hierarchies are not worth having (maintainer + ruling), the default is set for the deep case unconditionally — no level-count gate + and no `consistent_jacobian` gate. Four smoother iterations, not more: per unit work + gmres/4 (ρ^(1/4) = 0.87) beats gmres/8 (0.91). + +**What this does not fix.** At the hard small-`ξ` corner the inner solve fails under +every smoother tried — richardson outright, and gmres with **ρ > 1** (the V-cycle +diverges) at 4 levels, where the 3-level run still returned a finite ρ≈0.9. The +coarsest grid cannot represent the weak-zone/notch viscosity contrast, so geometric +interpolation cannot carry it: this is *operator conditioning*, and the lever for it is +the δ/ξ continuation (Layer 2), not the smoother. Do not expect a smoother to rescue +small ξ. + +Measured with `fmg_contraction_probe.py` / `smoother_depth_sweep.py` in the Spiegelman +hard-case study (ρ_MG = per-V-cycle contraction of the velocity block). + +## The user-facing result + +Once all three land: + +```python +stokes.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel +cm.Parameters.shear_viscosity_0 = ... +cm.Parameters.yield_stress = ... # a plain pressure-dependent yield +stokes.solve(homotopy=True) # warm-start, δ-descent, smoother — all automatic +if stokes.has_solution: + ... +``` + +The only knobs a user ever touches are the two deliberate opt-outs: force-fresh +(discard an existing solution) and fine-tune-the-march. + +## Implementation stages (for the dedicated session) + +1. **Layer 1a — `has_solution` property + Picard-on-cold warm-up** — **LANDED** + (`petsc_generic_snes_solvers.pyx`; test `test_0201`). Behind the existing + `zero_init_guess` semantics (no default flip). Additive, low-risk. As built: + `has_solution` lives on `SolverBaseClass`, is recorded by + `_record_convergence_status()` at the end of every `solve()` (the rotated + free-slip path records from its info dict, since it runs its own KSP loop), + and resets in the `is_setup = False` setter. The cold Picard warm-up is scoped + to a cold Stokes solve **under the consistent-Newton tangent** + (`consistent_jacobian is True`) — the opt-in nonlinear regime that needs it — + so the default (frozen) tangent path is bit-identical and the common linear + solve pays nothing. Broadening the warm-up to all nonlinear solvers (a cached + nonlinearity probe rather than the tangent proxy) is a natural Layer-1b + extension. Recipe + config-trap-list captured in the `nonlinear-solver` skill. +2. **Layer 1b — flip `zero_init_guess` to auto-detect** — **LANDED** (`test_0201`). + Tri-state: `None` (default) auto-detects from `has_solution`, `True` forces fresh, + `False` insists on warm. Both gates cleared before flipping — the free-surface + chain's repeated bare `solve()` is a *linear* Poisson mesh-displacement solve, so a + warm start cannot change its converged answer; the mover/adapt reset is covered by + the `is_setup` hook and its test. **Measured consequence:** warm and cold agree to + the *convergence tolerance*, not bitwise (2.4e-5 at `tol=1e-4`, 7.4e-10 at `1e-8`, + 2.0e-14 at `1e-12`), because an iterative solve stops anywhere in the tolerance + ball and a warm start enters it from a different direction. A test with a threshold + tighter than its own solver tolerance can therefore shift. +3. **Layer 3 — non-symmetry-safe smoother default** — **LANDED** (`gmres`+`sor`+ + `norm_type=none` in the FMG bundle; test `test_1014`). Independent of homotopy. + Benchmarked as described in the Layer 3 section: applied unconditionally rather + than gated on the consistent tangent, because the gain scales with multigrid depth + and shallow hierarchies are not worth special-casing. +4. **Layer 2 — `supports_yield_homotopy` / `_yield_homotopy_control` model hook and + the `solve(homotopy=True)` continuation** — **LANDED** (test `test_1057`). As built: + the control is a `YieldHomotopyControl` carrying a **model-owned δ setter** (the + isotropic models update the `constants[]` atom, TI-VEP rebuilds — and going through + the model's own `yield_softness` property is what stops a later + `_get_yield_softness()` resetting δ to a stale value), the recommended tangent, and + the δ atom. The march is residual-guided as designed, and additionally *reverts and + retries a failed δ more gently* before settling. A VEP model's `timestep` is + forwarded to the inner solves. **Cold-start finding:** a cold power-mean solve died + with `DIVERGED_FNORM_NAN`, and the first diagnosis (a `0/0` needing a strain-rate + floor) was **wrong**. At `ε̇=0` the plastic viscosity is `+inf`, which `Min` and the + sqrt soft-min carry correctly to the viscous branch — both cold-start fine. Only the + power-mean broke, because it formed its harmonic mean as + `η_ve·η_pl/(η_ve+η_pl)` = `inf/inf`. Rewriting that as the identical `η_ve/(1+f)` + fixes it with no floor and no numerical change. The bug was never homotopy-specific + or cold-start-specific: `η_pl` is infinite at *every rigid (unyielded) point*, so it + could poison a converged solve wherever a plug forms. The continuation logic existed as + a reference implementation (`underworld3.systems.yield_continuation`, the extracted + form of `convergence.py::run_continuation`); it is folded in, driven by the model hook + and Layer 1's warm-start. + +## Open verification points / risks + +- **Free-surface chain of solves** (on a feature branch, not yet inspected): a chain + of *deliberately independent* solves is the one pattern that might rely on the + implicit cold default. Confirm it either wants explicit force-fresh or is correct + warm, before flipping the default. +- **Adaptivity / mesh-mover reset**: confirm `has_solution` resets on the same + structural-invalidation event that already nukes the initial guess after a remesh. +- **Empirical basin test**: confirm a *single* Picard step at δ=δ₀ actually lands in + the Newton basin on the hard case — i.e. that it genuinely replaces the explicit + viscous pre-solve. Likely, given large-δ boundedness, but test it. +- **Benchmark the two core-default changes** (Layer 1b, Layer 3) against the existing + solver regression suite — "Solver Stability is Paramount". + +## References + +- Reference implementation of the continuation: `convergence.py::run_continuation` + in the Spiegelman hard-case study; extracted form in + `underworld3.systems.yield_continuation`. +- Consistent tangent + δ soft-min families: {doc}`jacobian-consistent-tangent`. +- Solver strategy catalogue (consult and contribute): `solver-strategies-catalogue`. +- Diagnostics: `SNES_*.get_snes_diagnostics()` / `check_snes_convergence()` / + `solve_with_diagnostics()` — `linear_iterations` is the diagnostic tell. +- The recipe and its trap list are also captured in the `plasticity-solvers` / + nonlinear-solver skill. diff --git a/docs/developer/index.md b/docs/developer/index.md index 977e09ac8..9c3ae16bd 100644 --- a/docs/developer/index.md +++ b/docs/developer/index.md @@ -160,6 +160,7 @@ design/SWARM_MODERNIZATION_DESIGN_2026-07 design/PROJECTED_NORMALS_API_DESIGN design/TURBULENCE_MODEL_DESIGN design/declined-coord-units-proposal +design/nonlinear-solver-homotopy-warmstart ``` ```{toctree} diff --git a/docs/reviews/2026-07/nonlinear-solver-homotopy-adversarial-review.md b/docs/reviews/2026-07/nonlinear-solver-homotopy-adversarial-review.md new file mode 100644 index 000000000..90e3a1eba --- /dev/null +++ b/docs/reviews/2026-07/nonlinear-solver-homotopy-adversarial-review.md @@ -0,0 +1,299 @@ +# Adversarial review — nonlinear solver warm-start + yield homotopy (4 layers) + +Branch `feature/nonlinear-warmstart-homotopy` vs `development` +(16 commits, 2402 insertions; `level_1 and tier_a` suite is green at 470 passed). + +Method: four independent adversarial reviewers, each given one dimension and +instructed to break the change rather than praise it, plus the author's own pass. +Every finding below was re-verified against the source before being reported; +findings the reviewers raised that did not survive verification are not listed. + +**Verdict: NOT ready to flag ready.** Three critical and seven major defects, most of +them in the new code rather than pre-existing. Two of the four layers (1a, 3) look +sound; the homotopy driver (Layer 2) and the default flip (Layer 1b) both need work. + +--- + +## CRITICAL + +### C1 — `entry_maxit` / `step_maxit` are inert; a failing δ-step burns 50 iterations +`systems/yield_continuation.py:151` sets `solver.petsc_options["snes_max_it"]`, but +`SNES_Stokes_SaddlePt.solve` hardcodes `snes_max_it = 50` +(`petsc_generic_snes_solvers.pyx:8474`) and pushes it with `setValue` + +`setFromOptions` immediately before the real solve (`:8543`) — the file's own comment +at `:8468` documents that this clobbers any user-set value. + +Consequences: both documented `homotopy_options` budget keys are silent no-ops; the +docstring promise that "a tight budget lets a too-hard step abort cheaply" is false +(a failing δ costs 50 Newton iterations, ×3 with the default retries); and the +residual-guided step logic at `:175` compares the *real* `nit` (cap 50) against a +*fictional* `budget` (10/30), so the "ease off" branch mis-fires and drives `step` +toward its 0.95 clamp. The mechanism that does work is +`snes.setTolerances(max_it=…)` after `setFromOptions` (as `_snes_solve_with_retries` +already does at `:1495`). + +### C2 — every retry after a failed δ is COLD, defeating the revert +`yield_continuation.py:181` reverts `u`/`p` to the last converged state so the retry +can warm-start from it. But the failed inner solve already ran +`_record_convergence_status()`, setting `has_solution = False`, so the retry's +`solver.solve(zero_init_guess=not solver.has_solution, …)` (`:155`) resolves to +**cold** — a from-zero solve at a δ *sharper* than the one that just failed. That is +exactly the regime this module's own docstring says does not converge. The revert is +wasted, the retry is near-guaranteed to fail, and the march settles earlier than it +should. + +### C3 — `homotopy=True` on a VEP model advances the stress history once per δ-step +`ViscoElasticPlasticFlowModel` advertises the homotopy +(`constitutive_models.py:2289`), and `SNES_Stokes.solve` forwards `timestep` into +every inner solve (`systems/solvers.py:1478`). Each inner solve therefore re-enters +the `has_stress_history` branch and runs the full time-integration tail — +`DFDt.update_pre_solve(timestep)`, the stress projection, the `psi_star` shift loop, +and `DFDt.update_post_solve(timestep)` (`solvers.py:1525-1585`). An N-step march +advances the elastic stress history by **N·dt for one requested dt**, and the +driver's revert restores only `u`/`p`, not the history, so a failed step leaves the +extra shift in place. `test_1057` only exercises the non-elastic model, so nothing +catches it. + +--- + +## MAJOR + +### M1 — `_force_setup=True` now warm-starts, inverting its own contract +All four solve bodies resolve `zero_init_guess` *before* the invalidation it depends +on (`pyx` resolve/invalidate pairs 3444/3450, 4494/4497, 5217/5220, 8394/8402): + +```python +stokes.solve() # converged -> has_solution True +stokes.solve(_force_setup=True) # resolves WARM off the stale flag, THEN rebuilds +``` + +The `is_setup` setter's own comment names "explicit `_force_setup`" as a structural +invalidation that must drop the warm-start claim. Before the flip this path was cold. +Fix: resolve after the `_force_setup` block. + +### M2 — a march that settles via the failure exit leaves `has_solution=False` on a good solution +On the retries-exhausted `break`, the fields are reverted to the settled solution and +δ is reset to `settled`, but nothing restores `has_solution` (left `False` by the +failed solve). The report says `converged: True` with a real `settled_delta`, yet the +next `solve()` auto-cold-starts and discards the continuation result — precisely the +advertised usage (`solve(homotopy=True)` once, then a time loop). `solve_report` has +the same problem: it describes the *failed* δ, not the settled one. + +### M3 — the march always cold-starts its entry solve +`yield_continuation.py:148` runs `solver.is_setup = False` on the first iteration, +whose setter clears `has_solution`, so `not solver.has_solution` on the next line is +unconditionally `True`. A time loop calling `solve(homotopy=True)` per step discards +the previous step's converged state every time, contradicting the module's own +comment ("Cold only if there is genuinely nothing to warm-start from"). + +### M4 — no `try/finally`: an exception mid-march leaves solver and model corrupted +`yield_continuation.py:145-203` has no exception handling. If an inner solve raises +(PETSc error, JIT failure, the VEP `timestep is None` ValueError), then `snes_max_it` +is left at the march value, `consistent_jacobian` at `control.tangent`, the model in +`softmin`/`powermean` at the *failed* δ, and `u`/`p` hold the diverged iterate with no +revert. + +### M5 — unrestored side effects on the user's solver and model +`solver.consistent_jacobian = control.tangent` (`:128`) is never restored: a user who +set `consistent_jacobian = "continuation"` (a supported value) silently gets `True` +or `False` forever after. `_yield_homotopy_control()` likewise leaves +`yield_mode="softmin"` and `yield_smoother="powermean"` permanently — a model the user +configured as `yield_mode="min"` (the `ViscoPlasticFlowModel` **default**) comes back +as a power-mean soft-min, so every later plain `solve()` runs different physics than +asked for, unwarned. + +### M6 — the power-mean cannot reach exact `Min`, but the docs say it does +Sharpness is `s = 1/(δ + 0.001)`, which **saturates at 1000**: + +| δ | 1e-2 | 1e-3 | 1e-4 | 1e-6 | 0 | +|---|---|---|---|---|---| +| s | 91 | 500 | 909 | 999 | **1000** | + +`constitutive_models.py:989` and `:1031` state "Both approach exact `Min` as δ → 0". +False for the power-mean family: δ=0 gives a finite power-mean (~0.2 % below true +`Min`), not the sharp surface. Since the entire homotopy premise is "march δ→0 to +reach the sharp yield surface", the contract needs restating: with `powermean` you +converge to a slightly smoothed yield law. Corollary (minor): the march happily +descends below δ≈1e-4 where nothing changes — a demo run reached δ=2.4e-10 in 9 steps, +and `settled_delta` then reads far sharper than the law actually is. + +### M7 — `TransverseIsotropicVEPFlowModel` advertises the homotopy but ignores half the control +It returns `supports_yield_homotopy = True` (`constitutive_models.py:3800`) and +inherits `_yield_homotopy_control()`, but never calls `_combine_yield`: its yield law +is inlined (`:3485`, `:3578`, `:3986`, `:4049`) reading `self._yield_softness` as a raw +float and hardcoding the **sqrt** family. So `yield_smoother = "powermean"` is a silent +no-op there — and the power-mean's large-δ harmonic limit is the stated justification +for the march's well-posed cold entry. Its `yield_softness` setter also never syncs +`_yield_softness_expr`, so the `control.delta` atom handed back is meaningless for +this model. + +--- + +## MINOR + +- **m1** — `uw.pprint(0, …)` (`yield_continuation.py:168,187,194`) uses the retired + `pprint_old(ranks, …)` shape; the current signature is `pprint(*args, proc=0, …)`, + so the literal `0` is printed. Confirmed empirically: every march line reads + `0 [yield-continuation] δ=…`. +- **m2** — diverged **linear** rotated free-slip records success. `pyx:8456` does + `.get("converged", True)`, but `solve_rotated_freeslip()` returns only + `ksp_reason`/`ksp_its` (`rotated_bc.py:379`), so the default fires even when + `_warn_if_ksp_diverged` just warned. Provably inconsistent with the sibling line: + `_capture_rotated_report` derives `converged = reason > 0` from the same dict, so + `solve_report.converged is False` while `has_solution is True`. +- **m3** — the retry log prints the wrong δ as having failed (`:194` recomputes + `d/step`, which is the last *good* δ). +- **m4** — `failures` is a whole-march counter, never reset on success, so `retries` + is not "per failed δ" as documented. +- **m5** — `delta0` is unvalidated (`down` is). `delta0 <= 0` reaches + `s = 1/(δ+0.001)` and divides by zero at δ = −0.001. +- **m6** — no cap on δ-steps: `step` is clamped to ≤0.95, so a hostile march can take + ~135 solves from `delta0=1.0` to `dmin=1e-3` with no `max_steps` escape (compounded + by C1's 50-iteration steps). +- **m7** — `solve(homotopy=True)` silently drops `picard`, `divergence_retries`, + `evalf`, `order`, `_force_setup`, `debug`, and even `zero_init_guess`. +- **m8** — two dispatch points with different behaviour: `pyx:8399` forwards no + `solve_kwargs`, so a `SNES_Stokes_SaddlePt` subclass that doesn't override `solve` + would give a VEP model inner solves with no `timestep`. +- **m9** — Charter §6: `YieldHomotopyControl` and `SolveReport` are named in public + docstrings but not exported from `systems/__init__` (deep-import-only). +- **m10** — Charter §4: three undocumented exception swallows in + `_capture_solve_report` (`pyx:1239-1250`). +- **m11** — stale docstrings after the flip: `SNES_Darcy.solve` still says + "If True (default)"; `SNES_TransientDarcy.solve` says "(default True)". +- **m12** — duplicate test number: `test_1055_solve_report.py` and + `test_1055_yield_smoother.py` are both new on this branch; `test_1057` then skips + 1056. +- **m13** — latent, pre-existing: the `yield_stress_min != 0` guard fixed on the DP + model survives at `constitutive_models.py:2030` (VEP) and `:3534` (TI-VEP). Same + class of bug; per Charter §2 it wants a `# TODO(BUG):` rather than a silent fix. + +--- + +## Test quality (Charter §8) — the review's sharpest hit on the author + +Three of the new tests would pass with the feature deleted: + +- `test_0201::test_zero_init_guess_is_tristate_and_auto_detects` asserts the **private** + `_resolve_zero_init_guess` back to itself; it never drives the public `solve()` path, + so it is blind to M1 (the `_force_setup` ordering bug) — the exact defect it should + have caught. +- `test_0201::test_repeated_default_solve_agrees_to_solver_tolerance` passes unchanged + if the flip is reverted (two cold solves of a linear Poisson also agree), and never + asserts the second solve was actually warm. +- `test_0201::test_cold_warmstart_under_consistent_newton_converges` claims to exercise + the automatic Picard branch but uses a **linear** `ViscousFlowModel`, which converges + cold regardless; it passes with the warm-up line deleted. +- `test_1057::test_solve_homotopy_marches_and_reports` asserts only `steps >= 1` and + `settled_delta <= delta0`, both satisfied by an implementation that does one solve + at δ₀ and stops — it never asserts the march descended. + +No parallel coverage was added for any of the four layers (Charter §11). Layer 3 is +the notable gap: `test_default_fmg_bundle_is_parallel_safe` is a *serial* test +asserting option strings, and never exercises the new `gmres`/`norm_type=none` +smoother at np>1 where it interacts with the redundant-LU coarse solve. + +--- + +## Constitutive maths — added by the fourth reviewer (most severe of the four) + +### C4 (CRITICAL) — the power-mean returns NaN wherever the yield stress goes negative +`constitutive_models.py:1005-1017`. With a pressure-dependent Drucker–Prager +`τ_y = C + sinφ·p`, tension drives `τ_y < 0` ⇒ `η_pl < 0` ⇒ `f < 0`, and `a = 1+f`, +`b = 1+1/f` are then negative bases raised to the non-integer power `-s` ⇒ **NaN** +(reproduced against the real UW3 expression: `η = nan`, all Jacobian entries NaN, +where `yield_mode="min"` returns a finite floored `0.001`). The old hard-`Min` path +degraded gracefully; the new smooth path does not. + +This is not a corner case for this feature: `_yield_homotopy_control()` flips the +model into `powermean` **unconditionally and without checking that the yield stress is +bounded below**, so `solve(homotopy=True)` NaNs on the first residual for any +pressure-dependent DP model that lacks a lower bound — i.e. the exact target problem. +(The Spiegelman driver wraps `sympy.Max(C + sinφ·p, 0)` by hand, which is why the +hard-case study never hit it.) + +### C5 (CRITICAL) — `_apply_floor(value, 0)` has a rounding scale of exactly zero +`constitutive_models.py:1088`, reached from `:1293`. The floor is rounded by +`ε = δ·floor`; with `floor = 0` that is `smooth_max(τ_y, 0, 0) = ½(τ_y + |τ_y|)` — +the exact kink the smooth floor exists to remove, with a `0/0` derivative. Over the +whole region where raw `τ_y ≤ 0` the floored value is exactly `0.0`, so `η_pl = 0`, +`η = 0`, and every Jacobian entry is NaN (reproduced). For a zero floor the new +"smooth" path is **strictly worse than the `sympy.Max` it replaced**, which is clean +there. `_apply_floor`'s own docstring admits it "needs a non-zero `floor` for a length +scale" and nothing enforces it — and the same commit's sentinel change +(`!= 0` → `!= -sympy.oo`) is what made `yield_stress_min = 0` reachable. + +### M8 (MAJOR) — δ doubles as the floor-rounding scale, so early march solves run a *different model* +`constitutive_models.py:1088`. `yield_continuation` starts at `delta0=1.0` (and the +`yield_smoother` setter forces δ→1.0 when switching to powermean), giving `ε = 1.0·floor` +and measured `smooth_max(F)/F = 1.500` — the viscosity and yield-stress floors sit up to +**50 % above** the values the user requested. It converges away as δ→dmin, but it means +the early continuation steps solve a perturbed problem, it is undocumented, and there is +no independent knob for the floor rounding. + +### M9 (MAJOR) — the cold-start claim is over-stated: the Newton *tangent* at ε̇=0 is NaN +`constitutive_models.py:1124-1128`. The comment I added ("no strain-rate floor is +needed … the soft-min carries it correctly") is true of the **residual** only: measured +at v=0, `η = 1.0000003` (finite — the inf/inf fix works) but **all five Jacobian entries +are NaN**. Partly intrinsic (`dε_II/dE = E/(2ε_II)` → 0/0), partly added by the +power-mean. It only works today because `solve()` interposes `picard = 1`, and that +guard is `consistent_jacobian is True` — so `consistent_jacobian="continuation"`, or +`zero_init_guess=False` on an all-zero field, assembles a NaN Jacobian. + +**This settles the open question I flagged to the maintainer**: a strain-rate floor +*is* legitimately needed — not for the residual (where removing it was right) but for +the **tangent**. + +### Corroborated: δ=0 ⇒ exact `Min` is false for the power-mean (see M6) +Measured deviation from exact `Min`: 5.2e-4 at δ=0, 1.2e-3 at δ=1e-3. Note the `+0.001` +floor **equals the default `dmin=1e-3`**, so the default march ends at `s = 500` — half +the achievable sharpness. The sqrt family is genuinely exact at δ=0; this is +power-mean-only. Six docstrings state it wrongly (`:653`, `:959`, `:995`, `:1039` +["`s = 1/δ`", but the code is `1/(δ+0.001)`], `:1211`, `:1353`). + +### Verified CLEAN by the same reviewer (against PETSc source) +- **The harmonic-mean rewrite is numerically sound.** `η_ve·η_pl/(η_ve+η_pl)` vs + `η_ve/(1+f)` agree to ≤1 ulp at (1e26,1e21) both orders, (1e26,1e26), (1e-20,1e26); + the product form only overflows above η≈1e154; the identity is exact. The new form is + strictly better — the only one that survives `η_pl = ∞`. My "no numerical change" + claim holds. +- **The FMG smoother bundle is correct.** `gmres` registers `KSP_NORM_NONE` for both + PC sides (`gmres.c:894`); `mg_levels_ksp_converged_maxits` is **required, not + redundant**, and is present — `KSPConvergedDefault` (`iterativ.c:1529`) tests it + *before* the `KSP_NORM_NONE` early return, so without it `KSPSolve_GMRES` would + hard-set `KSP_DIVERGED_ITS` and PCMG would flag PC failure. Restart 30 > max_it 4, so + the cycle never restarts. The GAMG delete list does include + `mg_levels_ksp_norm_type`. The fgmres claim is verified at both configuring sites. + Nit: the comment's "same four smoother iterations" understates cost — GMRES adds ~5 + Krylov vectors and Gram–Schmidt per level over Richardson's 2. + +### Reviewer claim NOT sustained +The reviewer states the power-mean NaN fix "ships no test that would catch a +regression". Not correct: `test_1057::test_cold_viscoplastic_solve_survives_zero_strain_rate` +parametrises a cold (ε̇=0) solve over all three yield modes and does fail without the +fix. The reviewer inspected `test_1055_yield_smoother.py` only. Their broader point +stands, though — **no** test covers the negative-`τ_y` path (C4) or the zero-floor path +(C5). + +--- + +## What held up + +- **Parallel correctness of the new control flow** — traced clean. `getConvergedReason` + is collective and rank-identical, so `_resolve_zero_init_guess` cannot split ranks; + the march's predicates are all rank-uniform; the revert is inside + `synchronised_array_update` on every rank. +- **Backward compatibility of the tri-state flip** — no caller in `src/`, `tests/`, + `docs/` or the skills passes `zero_init_guess` positionally or truth-tests it before + resolution. +- **`_record_convergence_status()` coverage** — reached on every normal exit of all + four solve bodies and the rotated early return; every wrapper funnels through them. +- **Recursion/reentrancy of `solve(homotopy=True)`** — `solve_kwargs` never carries + `homotopy`, so the inner solves cannot re-enter the march. +- **Return-value contract** — the dict survives `timing`, `memprobe` and + `SNES_Stokes_Constrained.solve(*args, **kwargs)`. +- **March arithmetic** — `step ∈ [0.05, 0.95]` strictly, `d` strictly positive and + decreasing, `d <= dmin` reachable, `delta0 <= dmin` exits correctly. The problem is + cost (m6), not termination. +- **Layers 1a and 3 as such** — `has_solution`'s lifecycle and the FMG smoother bundle + drew no correctness findings beyond the ones listed. diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 7d37eecfd..35f34a641 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -645,6 +645,18 @@ def requires_stress_history(self): """ return False + @property + def supports_yield_homotopy(self): + """Whether this model can be solved by a single-parameter yield homotopy. + + ``True`` on the yielding models, which carry a δ-parameterised soft-min + yield law that sharpens to the exact ``Min`` as δ→0 — see + :meth:`_yield_homotopy_control`. ``solver.solve(homotopy=True)`` refuses a + model that returns ``False`` (a purely viscous model has no yield surface to + sharpen, so there is nothing to march). + """ + return False + @property def stress_history_ddt_kwargs(self): """Extra kwargs passed to the auto-DDt creation when this model @@ -924,6 +936,227 @@ def _object_viewer(self): ) ) + # --- Yield soft-min smoother (shared by the visco-plastic subclasses) ----------- + # The δ soft-min regularisation and the smooth-min FAMILY selection live on the + # base class so every yielding model inherits one implementation. δ is held as a + # constants[] UWexpression atom (not a baked float) so a homotopy can ramp it at + # runtime via PetscDSSetConstants with no JIT recompile. + + def _get_yield_softness(self): + r"""The soft-min regularisation δ as a ``constants[]`` UWexpression atom. + + Created lazily and kept in sync with ``self._yield_softness`` (the + configured numeric value). Storing δ as a UWexpression — rather than + baking the float into the compiled flux — lets a yield homotopy ramp + δ at runtime via ``PetscDSSetConstants`` with no JIT recompile. + ``δ = 0`` makes the sqrt law identically ``Min``. + """ + delta_value = getattr(self, "_yield_softness", 0.0) + if getattr(self, "_yield_softness_expr", None) is None: + self._yield_softness_expr = expression( + R"{\updelta_{y}}", + sympy.Float(delta_value), + "Yield soft-min regularisation δ (rampable constant; δ=0 ⇒ exact Min)", + ) + # Onset offset (-1+√(1+δ²))/2 keeps g(0)=1 (no spurious yield below + # onset). Held as its OWN constant atom — a single symbol in the + # stress tensor — so it does not blow the tensor up, while still + # tracking δ symbolically (one δ update repacks both constants). + self._yield_offset_expr = expression( + R"{\updelta_{y,0}}", + (-1 + sympy.sqrt(1 + self._yield_softness_expr**2)) / 2, + "Yield soft-min onset offset; tracks δ so g(0)=1 exactly", + ) + else: + self._yield_softness_expr.sym = sympy.Float(delta_value) + return self._yield_softness_expr + + def _get_yield_offset(self): + """The onset-offset constant atom (lazily created alongside δ).""" + if getattr(self, "_yield_offset_expr", None) is None: + self._get_yield_softness() + return self._yield_offset_expr + + def _combine_yield(self, eta_ve, eta_pl): + r"""Combine the visco-elastic/viscous viscosity ``eta_ve`` with the plastic + (yield) viscosity ``eta_pl`` according to ``self._yield_mode``: + + - ``"min"``: exact hard ``Min(η_ve, η_pl)`` (sharp yield surface). + - ``"harmonic"``: ``1/(1/η_ve + 1/η_pl)`` (a distinct smooth blend). + - ``"softmin"``: the δ-parameterised soft-min, in the family chosen by + ``self.yield_smoother`` — ``"sqrt"`` (default; overshoots τ_y in the + transition) or ``"powermean"`` (undershoots τ_y, ``η_eff ≤ Min`` always). + The ``sqrt`` family is exactly ``Min`` at ``δ = 0``; the ``powermean`` + family is not — its sharpness ``s = 1/(δ + 0.001)`` saturates at 1000, so it + approaches ``Min`` only to about 0.3 %. Measured on a 45 %-yielded box, that + leaves the converged solution satisfying the exact yield law to 6e-10 of the + initial residual — an order of magnitude INSIDE a 1e-8 solver tolerance, so + the difference is not observable in practice (2026-07-27). + + Behaviour at the stock settings (``yield_mode`` per model default, + ``yield_smoother="sqrt"``) is identical to the previous inline law — δ merely + moves from a baked float to a ``constants[]`` atom (same value). + """ + mode = getattr(self, "_yield_mode", "softmin") + if mode == "harmonic": + return 1 / (1 / eta_ve + 1 / eta_pl) + if mode == "min": + return sympy.Min(eta_ve, eta_pl) + + # "softmin": δ-parameterised smooth-min family. + smoother = getattr(self, "_yield_smoother", "sqrt") + delta = self._get_yield_softness() + f = eta_ve / eta_pl + if smoother == "powermean": + # p-norm soft-min in an overflow-safe harmonic-normalised form. δ is + # floored SMOOTHLY (+ε, not Max()) so 1/δ stays finite as δ→0 (a Max on + # the δ atom triggers an unsupported symbolic numeric comparison). + s = 1 / (delta + sympy.Float(0.001)) + a = 1 + f + b = 1 + 1 / f + # Harmonic mean written as eta_ve/(1+f), NOT eta_ve*eta_pl/(eta_ve+eta_pl). + # The two are algebraically identical, but the product-over-sum form + # evaluates to inf/inf = NaN when eta_pl is infinite — which is exactly + # what a rigid (unyielded) point gives, since eta_pl = tau_y/(2 edot_II) + # and edot_II = 0 there. That includes every point of a cold v=0 start. + # In this form f -> 0 and N -> eta_ve, the correct viscous limit. + N = eta_ve / a + return N * (a ** (-s) + b ** (-s)) ** (-1 / s) + + # default "sqrt" soft-min: η_ve / g(f), g(0)=1, g ≈ max(1, f), exact Min at δ=0. + offset = self._get_yield_offset() + g = 1 + (f - 1 + sympy.sqrt((f - 1) ** 2 + delta**2)) / 2 - offset + return eta_ve / g + + @property + def yield_smoother(self): + r"""Which smooth-min FAMILY regularises the ``"softmin"`` yield mode. + + Both families use the same softness parameter ``δ`` (``yield_softness``) + and approach exact ``Min`` as ``δ → 0``, but differ in how they round + the kink: + + - ``"sqrt"`` (default): ``η_ve / g(f, δ)`` with + ``g = 1 + ½(f−1+√((f−1)²+δ²)) − offset``. Exact ``Min`` at ``δ=0``; + **overshoots** the yield surface in the transition (carries stress a + few–60 % above ``τ_y`` before asymptoting). + - ``"powermean"``: the p-norm soft-min + ``η_eff = (η_ve^(−s) + η_pl^(−s))^(−1/s)`` with ``s = 1/δ`` + (``s=1`` ⇒ harmonic mean; ``s→∞`` ⇒ ``Min``). **Undershoots** the + yield surface (``η_eff ≤ Min`` always — approaches ``τ_y`` strictly + from below, never over-yields). Computed in an overflow-safe + harmonic-normalised form for geodynamic viscosity ranges. + + Selecting ``"powermean"`` bumps a zero ``yield_softness`` to ``1.0`` + (``s=1``, the parameter-free harmonic mean) since ``δ=0`` (``s=∞``) is + the singular hard-``Min`` limit it only *approaches*. + """ + return getattr(self, "_yield_smoother", "sqrt") + + @yield_smoother.setter + def yield_smoother(self, value): + if value not in ("sqrt", "powermean"): + raise ValueError( + f"yield_smoother must be 'sqrt' or 'powermean', got '{value}'" + ) + self._yield_smoother = value + if value == "powermean" and getattr(self, "_yield_softness", 0.0) == 0.0: + # δ=0 ⇒ s=1/δ=∞ is the singular Min limit; default to the + # parameter-free harmonic mean (s=1) instead. + self.yield_softness = 1.0 + self._reset() + + # --- Smooth lower bounds (viscosity / yield floors) ----------------------------- + # A hard sympy.Max cutoff is non-differentiable at the corner. When the flux is + # differentiated for the consistent-Newton tangent that kink breaks the tangent — + # and, because one operand is a UWexpression over the fields, sympy's fuzzy `>=` + # comparison cannot resolve it and recurses. In the smooth yield modes we therefore + # round the floor with the same δ that regularises the yield transition, so the + # whole effective viscosity stays differentiable and δ→0 recovers the sharp bound. + + def _apply_floor(self, value, floor): + r"""Impose the lower bound :math:`value \ge floor`. + + In ``yield_mode="min"`` this is the exact hard ``sympy.Max(value, floor)`` — + the sharp cutoff the model has always used. In the smooth yield modes + (``"softmin"``/``"harmonic"``) it is the differentiable ``uw.maths.smooth_max`` + rounded by the yield softness δ *relative to the floor* (:math:`\epsilon = + \delta\,|floor|`), so the whole effective viscosity — not only the yield + transition — is differentiable for the consistent-Newton tangent. The relative + rounding needs a non-zero ``floor`` for a length scale, which the numerical + viscosity/yield floors provide; a cutoff *at zero* (a tension cutoff) has no + such scale and is rounded by its own physical parameter via + ``uw.maths.smooth_max`` at the call site instead. + """ + # smooth_max is 1/2 (a + b + sqrt((a-b)^2 + eps^2)) — pure arithmetic on the + # operands, so `floor` stays the symbolic Parameter that the JIT routes + # through constants[], and sympy is never asked for the ordering test that + # recurses on an opaque UWexpression. At eps = 0 this is exactly + # Max(value, floor), but written without a comparison. + # + # TODO(DESIGN): the rounding scale is RELATIVE (delta * floor), so it + # collapses for a tension cutoff at floor = 0 — now the default for + # yield_stress_min. That leaves a hard corner: the floored yield stress is + # exactly 0 in tension, hence eta_pl = tau_y/(2 edot_II) is exactly 0 and the + # tangent through the soft-min is undefined there. A properly rounded cap + # (Griffith / parabolic) needs an ABSOLUTE stress scale, which this signature + # cannot supply. Maintainer decision pending (2026-07-26). + rounding = 0 if getattr(self, "_yield_mode", "min") == "min" \ + else self._get_yield_softness() * floor + return uw.maths.smooth_max(value, floor, rounding) + + # Tangent this model wants while the yield homotopy marches. Newton is right for + # a purely viscous-plastic yield; the elastic (VEP) subclasses override to the + # frozen/Picard tangent, because the consistent yield tangent taken across the + # elastic stress-history block makes the Jacobian indefinite and the linear + # solve fails outright (DIVERGED_LINEAR_SOLVE). + _yield_homotopy_tangent = True + + def _yield_homotopy_control(self): + """Put this model in its smooth (δ-parameterised) yield mode and describe + how to march it. + + The model owns what the homotopy *means* for it: which knob is the + continuation parameter, how to set it, and which tangent to pair it with. + The solver only marches the number. Called by + ``solver.solve(homotopy=True)``; see + :doc:`nonlinear-solver-homotopy-warmstart` (Layer 2). + + Selects the power-mean soft-min family, whose large-δ limit is the harmonic + mean — bounded by the background viscosity even as :math:`\\dot\\varepsilon + \\to 0`, so the first (cold) solve of the march is well posed and no separate + viscous pre-solve is needed. + + Returns + ------- + YieldHomotopyControl + ``set_delta`` (model-owned setter for δ), ``tangent`` (the + ``consistent_jacobian`` value to use), and ``delta`` (the ``constants[]`` + atom itself, for diagnostics). + """ + from underworld3.systems.yield_continuation import YieldHomotopyControl + + self.yield_mode = "softmin" + self.yield_smoother = "powermean" + + # No strain-rate floor is needed for the cold (v = 0) start the march begins + # from: eta_pl = tau_y/(2 edot_II) is +inf there, which the soft-min carries + # correctly to the viscous branch. See the harmonic-mean note in + # _combine_yield for the one form that must be written carefully to keep it + # so. + + def set_delta(value): + # Go through the property, not the atom: `yield_softness` updates BOTH + # the stored value and the constants[] atom, so a later + # _get_yield_softness() cannot silently reset δ to a stale number. + self.yield_softness = value + + return YieldHomotopyControl( + set_delta=set_delta, + tangent=self._yield_homotopy_tangent, + delta=self._get_yield_softness(), + ) + ## NOTE - retrofit VEP into here @@ -988,6 +1221,16 @@ def __init__(self, unknowns, material_name: str = None): "Effective viscosity (plastic)", ) + # 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 + self._yield_smoother = "sqrt" # smooth-min family: "sqrt" | "powermean" + self._yield_softness_expr = None # constants[] δ atom (created lazily) + self._yield_offset_expr = None # onset-offset atom (created lazily) + class _Parameters(_ParameterBase, _ViscousParameterAlias): """Any material properties that are defined by a constitutive relationship are collected in the parameters which can then be defined/accessed by name in @@ -1026,7 +1269,7 @@ class _Parameters(_ParameterBase, _ViscousParameterAlias): yield_stress_min = api_tools.Parameter( R"{\tau_{y, \mathrm{min}}}", - lambda inner_self: -sympy.oo, + lambda inner_self: 0, "Yield stress (DP) minimum cutoff", units="Pa", ) @@ -1062,28 +1305,48 @@ def viscosity(self): # Don't put conditional behaviour in the constitutive law # when it is not needed - if inner_self.yield_stress_min.sym != 0: - yield_stress = sympy.Max(inner_self.yield_stress_min, inner_self.yield_stress) + # Lower bound on the yield stress, defaulting to ZERO and therefore normally + # active. tau_y is compared against the second invariant of the stress, so a + # negative tau_y is meaningless — a pressure-dependent Drucker-Prager yield + # C + sin(phi)*p goes negative in tension and must be cut off there rather + # than propagated into tau_y/(2 edot_II). (The ±oo defaults elsewhere in + # Parameters exist so sympy can cancel an unused term away; that trick is + # wrong here, so this one defaults to 0 — maintainer ruling 2026-07-26.) + # An explicit -oo still disables the floor. + if inner_self.yield_stress_min.sym != -sympy.oo: + yield_stress = self._apply_floor( + inner_self.yield_stress, inner_self.yield_stress_min + ) else: yield_stress = inner_self.yield_stress - viscosity_yield = yield_stress / (2 * self._strainrate_inv_II) - - ## Question is, will sympy reliably differentiate something - ## with so many Max / Min statements. The smooth version would - ## be a reasonable alternative: - - # effective_viscosity = sympy.sympify( - # 1 / (1 / inner_self.shear_viscosity_0 + 1 / viscosity_yield), - # ) - - effective_viscosity = sympy.Min(inner_self.shear_viscosity_0, viscosity_yield) + # Rate regularisation. eta_pl = tau_y / (2 edot_II) is unbounded as edot -> 0; + # adding a floor to the strain rate caps it at tau_y/(2 edot_min) and so bounds + # the viscosity CONTRAST, which is what conditions the velocity block (the same + # role the Perzyna/rate-strengthening xi plays in the Spiegelman studies). The + # parameter was declared on this model but never applied — the elastic models + # (ViscoElasticPlastic, TransverseIsotropicVEP) have always used it. Default 0 + # = off, so the unregularised law is unchanged. + if inner_self.strainrate_inv_II_min.sym != 0: + viscosity_yield = yield_stress / ( + 2 * (self._strainrate_inv_II + inner_self.strainrate_inv_II_min) + ) + 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. + effective_viscosity = self._combine_yield( + inner_self.shear_viscosity_0, viscosity_yield + ) # If we want to apply limits to the viscosity but see caveat above # Keep this as an sub-expression for clarity if inner_self.shear_viscosity_min.sym != -sympy.oo: - self._plastic_eff_viscosity._sym = sympy.Max( + self._plastic_eff_viscosity._sym = self._apply_floor( effective_viscosity, inner_self.shear_viscosity_min ) @@ -1093,6 +1356,50 @@ def viscosity(self): # Returns an expression that has a different description return self._plastic_eff_viscosity + @property + def supports_yield_homotopy(self): + """This model carries a δ-parameterised yield law — ``solve(homotopy=True)`` + can march it. See :meth:`_yield_homotopy_control`.""" + return True + + @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. + - ``"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 + + @yield_mode.setter + def yield_mode(self, value): + if value not in ("min", "harmonic", "softmin"): + raise ValueError( + f"yield_mode must be 'min', 'harmonic', or 'softmin', got '{value}'" + ) + self._yield_mode = value + self._reset() + + @property + def yield_softness(self): + r"""Soft-min regularisation δ for ``yield_mode="softmin"`` (0 ⇒ exact Min). + + δ is held as a ``constants[]`` atom, so ramping it at runtime (this setter, + or ``cm._get_yield_softness().sym = ...`` + ``solver._update_constants()``) + does not trigger a JIT recompile — the basis of a scalable yield homotopy. + """ + return self._yield_softness + + @yield_softness.setter + def yield_softness(self, value): + self._yield_softness = float(value) + if getattr(self, "_yield_softness_expr", None) is not None: + self._yield_softness_expr.sym = sympy.Float(self._yield_softness) + self._reset() + def plastic_correction(self) -> float: r"""Scaling factor to reduce stress to yield surface. @@ -1242,6 +1549,9 @@ def __init__(self, unknowns, order=1, integrator: str = "bdf", self._order = order self._yield_mode = "softmin" # "min", "harmonic", "smooth", or "softmin" 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) + self._yield_offset_expr = None # onset-offset atom (created lazily) # Timestep — set by the solver before each solve(). Not a user parameter. # Initialised to oo (viscous limit). The solver overwrites this with the @@ -1323,7 +1633,7 @@ def dt_elastic(inner_self, value): yield_stress_min = api_tools.Parameter( R"{\tau_{y, \mathrm{min}}}", - lambda inner_self: -sympy.oo, + lambda inner_self: 0, "Yield stress (DP) minimum cutoff", units="Pa", ) @@ -1701,23 +2011,13 @@ def viscosity(self): if self.is_viscoplastic: vp_effective_viscosity = self._plastic_effective_viscosity - if self._yield_mode == "harmonic": - effective_viscosity = 1 / (1 / effective_viscosity + 1 / vp_effective_viscosity) - elif self._yield_mode == "softmin": - # Smooth approximation to Min(η_ve, η_pl): - # η_eff = η_ve / g(f) - # g(f) = 1 + softplus(f-1) - softplus(-1) ≈ max(1, f) - # where softplus(x) = (x + √(x² + δ²))/2 and f = η_ve/η_pl. - # Corrected so g(0) = 1 exactly (no spurious yield below onset). - # Approaches exact Min as δ→0. No Min/Max in expression. - delta = self._yield_softness - f = effective_viscosity / vp_effective_viscosity - import math # float offset avoids sympy expression blowup in tensor - offset = (-1 + math.sqrt(1 + delta**2)) / 2 - g = 1 + (f - 1 + sympy.sqrt((f - 1)**2 + delta**2)) / 2 - offset - effective_viscosity = effective_viscosity / g - else: - effective_viscosity = sympy.Min(effective_viscosity, vp_effective_viscosity) + # Combine η_ve with the plastic viscosity per yield_mode (harmonic / exact + # Min / δ-soft-min). The soft-min softness δ now lives in a constants[] atom + # (see _combine_yield / _get_yield_softness) — value-identical to the former + # inline float law at the same δ, but runtime-rampable with no recompile. + effective_viscosity = self._combine_yield( + effective_viscosity, vp_effective_viscosity + ) # Apply viscosity floor — but skip for smooth-blend yield modes # where the outer Max creates a nested Min/Max that breaks the @@ -1762,10 +2062,19 @@ def _plastic_effective_viscosity(self): "Strain rate 2nd Invariant including elastic strain rate term", ) - if parameters.yield_stress_min.sym != 0: - yield_stress = sympy.Max( - parameters.yield_stress_min, parameters.yield_stress - ) # .rewrite(sympy.Piecewise) + # Guard on the DISABLING sentinel, not on zero: zero is the default and a + # physically meaningful floor (the yield stress is compared against the second + # invariant of the stress, so a negative tau_y is meaningless). Only an + # explicit -oo turns the floor off. + if parameters.yield_stress_min.sym != -sympy.oo: + # Literal 0 for the default floor, not the parameter atom: sympy cannot + # fuzzy-compare an opaque UWexpression and Max canonicalisation recurses. + # smooth_max keeps both operands symbolic (the JIT routes the floor + # through constants[]) and needs no ordering test, which sympy cannot + # resolve against an opaque UWexpression. eps = 0 makes it exactly Max. + yield_stress = uw.maths.smooth_max( + parameters.yield_stress, parameters.yield_stress_min, 0 + ) else: yield_stress = parameters.yield_stress @@ -2009,7 +2318,11 @@ def yield_softness(self): @yield_softness.setter def yield_softness(self, value): - self._yield_softness = value + self._yield_softness = float(value) + # Keep the constants[] δ atom in sync (created lazily on first use) so a + # numeric δ assignment is reflected without a JIT recompile. + if getattr(self, "_yield_softness_expr", None) is not None: + self._yield_softness_expr.sym = sympy.Float(self._yield_softness) self._reset() @property @@ -2017,6 +2330,18 @@ def requires_stress_history(self): """VEP models always require stress history tracking.""" return True + @property + def supports_yield_homotopy(self): + """This model carries a δ-parameterised yield law — ``solve(homotopy=True)`` + can march it. See :meth:`_yield_homotopy_control`.""" + return True + + # Picard, not Newton: the consistent yield tangent taken across the elastic + # stress-history block makes the Jacobian indefinite, and the linear solve fails + # outright (DIVERGED_LINEAR_SOLVE at 0 iterations). The frozen tangent is + # contractive, and with the δ-march it still converges to the exact yield surface. + _yield_homotopy_tangent = False + @property def plastic_fraction(self): """Fraction of strain rate that is plastic: 1 - η_vep / η_ve.""" @@ -2887,7 +3212,7 @@ def dt_elastic(inner_self, value): ) yield_stress_min = api_tools.Parameter( R"{\tau_{y, \mathrm{min}}}", - lambda inner_self: -sympy.oo, + lambda inner_self: 0, "Yield stress minimum cutoff", units="Pa", ) strainrate_inv_II_min = api_tools.Parameter( @@ -3201,17 +3526,12 @@ def viscosity(self): if self.is_viscoplastic: vp_eff = self._plastic_effective_viscosity - if self._yield_mode == "harmonic": - eta_1_eff = 1 / (1 / eta_1_eff + 1 / vp_eff) - elif self._yield_mode == "softmin": - delta = self._yield_softness - f = eta_1_eff / vp_eff - import math # float offset avoids sympy expression blowup in tensor - offset = (-1 + math.sqrt(1 + delta**2)) / 2 - g = 1 + (f - 1 + sympy.sqrt((f - 1)**2 + delta**2)) / 2 - offset - eta_1_eff = eta_1_eff / g - else: - eta_1_eff = sympy.Min(eta_1_eff, vp_eff) + # Same yield envelope as the isotropic models: _combine_yield owns the + # harmonic / soft-min / hard-Min choice, reads delta from the + # constants[] atom (so a homotopy can ramp it without a recompile) and + # honours yield_smoother. Previously inlined here, which pinned this + # model to the sqrt family and to a baked float delta. + eta_1_eff = self._combine_yield(eta_1_eff, vp_eff) return inner_self.shear_viscosity_0 @@ -3250,8 +3570,13 @@ def _plastic_effective_viscosity(self): gamma_dot_abs = sympy.sqrt(sympy.Max(gamma_dot_sq, 0)) tau_y = parameters.yield_stress - if parameters.yield_stress_min.sym != 0: - tau_y = sympy.Max(parameters.yield_stress_min, tau_y) + # Guard on the DISABLING sentinel, not on zero — zero is the default and a + # real floor (a negative yield stress is meaningless against an invariant). + if parameters.yield_stress_min.sym != -sympy.oo: + # Literal 0 for the default floor (sympy fuzzy-compare recursion on atoms). + # smooth_max: keeps the floor symbolic, no ordering test (see the note + # in Constitutive_Model._apply_floor). eps = 0 is exactly Max. + tau_y = uw.maths.smooth_max(tau_y, parameters.yield_stress_min, 0) if parameters.strainrate_inv_II_min.sym != 0: viscosity_yield = tau_y / ( @@ -3294,17 +3619,12 @@ def _eta_for_tensor(self, integrator_mode, apply_yield): if apply_yield and self.is_viscoplastic: vp_eff = self._plastic_effective_viscosity - if self._yield_mode == "harmonic": - eta_1_eff = 1 / (1 / eta_1_eff + 1 / vp_eff) - elif self._yield_mode == "softmin": - delta = self._yield_softness - f = eta_1_eff / vp_eff - import math - offset = (-1 + math.sqrt(1 + delta**2)) / 2 - g = 1 + (f - 1 + sympy.sqrt((f - 1)**2 + delta**2)) / 2 - offset - eta_1_eff = eta_1_eff / g - else: - eta_1_eff = sympy.Min(eta_1_eff, vp_eff) + # Same yield envelope as the isotropic models: _combine_yield owns the + # harmonic / soft-min / hard-Min choice, reads delta from the + # constants[] atom (so a homotopy can ramp it without a recompile) and + # honours yield_smoother. Previously inlined here, which pinned this + # model to the sqrt family and to a baked float delta. + eta_1_eff = self._combine_yield(eta_1_eff, vp_eff) return eta_0, eta_1_eff def _assemble_c_tensor(self, eta_0, eta_1_eff): @@ -3508,7 +3828,11 @@ def yield_softness(self): @yield_softness.setter def yield_softness(self, value): - self._yield_softness = value + self._yield_softness = float(value) + # Keep the constants[] atom in step: this model now shares the isotropic + # _combine_yield, which reads δ from that atom rather than from the float. + if getattr(self, "_yield_softness_expr", None) is not None: + self._yield_softness_expr.sym = sympy.Float(self._yield_softness) self._reset() @property @@ -3516,6 +3840,19 @@ def requires_stress_history(self): """Transverse isotropic VEP requires stress history tracking.""" return True + @property + def supports_yield_homotopy(self): + """This model carries a δ-parameterised yield law — ``solve(homotopy=True)`` + can march it. It shares the isotropic yield envelope (:meth:`_combine_yield`), + so ``yield_smoother`` applies here too; its ``yield_softness`` setter still + triggers a rebuild, so a δ step costs a recompile the isotropic models avoid. + See :meth:`_yield_homotopy_control`.""" + return True + + # Picard, not Newton — as for the isotropic VEP model, the consistent yield + # tangent over the elastic stress-history block is indefinite. + _yield_homotopy_tangent = False + @property def plastic_fraction(self): """Fraction of strain rate that is plastic.""" @@ -3690,17 +4027,12 @@ def _eta_par_eff(self): return eta_par vp_eff = self._plastic_effective_viscosity - if self._yield_mode == "harmonic": - return 1 / (1 / eta_par + 1 / vp_eff) - elif self._yield_mode == "softmin": - delta = self._yield_softness - f = eta_par / vp_eff - import math - offset = (-1 + math.sqrt(1 + delta**2)) / 2 - g = 1 + (f - 1 + sympy.sqrt((f - 1) ** 2 + delta ** 2)) / 2 - offset - return eta_par / g - else: - return sympy.Min(eta_par, vp_eff) + # Same yield envelope as the isotropic models: _combine_yield owns the + # harmonic / soft-min / hard-Min choice, reads delta from the + # constants[] atom (so a homotopy can ramp it without a recompile) and + # honours yield_smoother. Previously inlined here, which pinned this + # model to the sqrt family and to a baked float delta. + return self._combine_yield(eta_par, vp_eff) def _eta_par_eff_lagged(self): """Yield-clipped ``η_∥_eff`` using the **lagged** strain rate @@ -3753,17 +4085,12 @@ def _eta_par_eff_lagged(self): 2 * (gamma_dot_abs_lag + sympy.Float(edot_min_val)) ) - if self._yield_mode == "harmonic": - return 1 / (1 / eta_par + 1 / vp_eff_lag) - elif self._yield_mode == "softmin": - delta = self._yield_softness - f = eta_par / vp_eff_lag - import math - offset = (-1 + math.sqrt(1 + delta ** 2)) / 2 - g = 1 + (f - 1 + sympy.sqrt((f - 1) ** 2 + delta ** 2)) / 2 - offset - return eta_par / g - else: - return sympy.Min(eta_par, vp_eff_lag) + # Same yield envelope as the isotropic models: _combine_yield owns the + # harmonic / soft-min / hard-Min choice, reads delta from the + # constants[] atom (so a homotopy can ramp it without a recompile) and + # honours yield_smoother. Previously inlined here, which pinned this + # model to the sqrt family and to a baked float delta. + return self._combine_yield(eta_par, vp_eff_lag) def _build_split_c_tensors(self, eta_perp, eta_par): r"""Build ``C_⊥ = 2·η_⊥·P_⊥`` and ``C_∥ = 2·η_∥·P_∥``. diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 0fda2ebd4..c4ab10ca0 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -85,6 +85,14 @@ class SolverBaseClass(uw_object): self._needs_bc_reregister = True self._needs_function_rewire = True + # Warm-start status, public (read-only) via `has_solution`. Set True only + # after a converged solve; reset on a structural rebuild (is_setup=False) + # so a remesh / adapt / mesh-mover never warm-starts off stale field + # data. Kept through coefficient changes (viscosity, yield softness δ, BC + # values, time step). See has_solution / _record_convergence_status and + # docs/developer/design/nonlinear-solver-homotopy-warmstart.md (Layer 1). + self._has_solution = False + self.Unknowns = self._Unknowns(self) self._order = 0 @@ -613,13 +621,28 @@ class SolverBaseClass(uw_object): opts[f"{prefix}pc_type"] = "mg" opts[f"{prefix}pc_mg_type"] = "full" # FMG (F-cycle) opts[f"{prefix}pc_mg_galerkin"] = "both" # RAP coarse operators - # richardson+sor (not chebyshev): chebyshev needs eigenvalue - # estimates of the smoothed operator, which are fragile on the - # indefinite / variable-viscosity Stokes velocity block and diverge; - # richardson+sor is the benchmark-validated, mesh-independent choice. - opts[f"{prefix}mg_levels_ksp_type"] = "richardson" + # gmres+sor, sized for a DEEP hierarchy (the only kind worth having: + # a two-level cycle is a coarse-grid correction, not a V-cycle, and is + # not worth special-casing). Chebyshev needs eigenvalue estimates of the + # smoothed operator, which are fragile on the indefinite / + # variable-viscosity velocity block and diverge. Richardson is + # stationary and degrades on the NON-SYMMETRIC operator produced by the + # consistent-Newton tangent. Measured on the Spiegelman notch (Drucker- + # Prager, eta contrast 1e26) over a nested 4-level hierarchy: contraction + # per V-cycle rho = 0.75 (richardson) vs 0.56 (gmres) at the SAME four + # smoother iterations -- and the gmres margin GROWS with depth (5% at 3 + # levels, 25% at 4), because deeper cycles apply the smoother on more + # coarse operators. Four iterations, not more: per unit work gmres/4 + # (rho^(1/4) = 0.87) beats gmres/8 (0.91). + opts[f"{prefix}mg_levels_ksp_type"] = "gmres" opts[f"{prefix}mg_levels_pc_type"] = "sor" opts[f"{prefix}mg_levels_ksp_max_it"] = 4 + # Run EXACTLY max_it smoother iterations: no residual-norm computation + # and no convergence test, so every V-cycle costs the same. A Krylov + # smoother makes the cycle non-stationary, which is why the velocity + # block is fgmres (flexible) rather than gmres -- see the fieldsplit + # defaults in the Stokes __init__. + opts[f"{prefix}mg_levels_ksp_norm_type"] = "none" opts[f"{prefix}mg_levels_ksp_converged_maxits"] = None # redundant+lu, not bare lu: a bare serial LU cannot factor a # distributed coarse matrix and fails at np>1 (DIVERGED_LINEAR_SOLVE @@ -650,7 +673,8 @@ class SolverBaseClass(uw_object): opts[f"{prefix}mg_levels_ksp_converged_maxits"] = None # Clear stale geometric-MG-only keys. for key in ("pc_mg_galerkin", "mg_levels_ksp_type", - "mg_levels_pc_type", "mg_coarse_pc_type", + "mg_levels_pc_type", "mg_levels_ksp_norm_type", + "mg_coarse_pc_type", "mg_coarse_redundant_pc_type"): opts.delValue(f"{prefix}{key}") self._pc_managed_value = "gamg" @@ -790,6 +814,110 @@ class SolverBaseClass(uw_object): self._needs_dm_rebuild = True self._needs_bc_reregister = True self._needs_function_rewire = True + # A structural invalidation (mesh change / adaptivity / mesh-mover / + # explicit _force_setup) means any stored solution no longer matches + # the operators — drop the warm-start claim so the next solve() + # cold-starts rather than warming off a stale iterate. Coefficient-only + # updates (new viscosity, δ, BC values, time step) route through + # _update_constants / a direct _needs_function_rewire and keep + # has_solution, so continuation and time-stepping warm-start correctly. + self._has_solution = False + + @property + def has_solution(self): + """``True`` when the solver holds a converged solution usable as a warm start. + + Set ``True`` only after a solve whose SNES reported a converged reason + (``> 0``); ``False`` initially, after a diverged solve, and after any + structural rebuild (mesh change / adaptivity / mesh-mover / explicit + ``_force_setup`` — the ``is_setup = False`` invalidation hook). It + survives coefficient changes (viscosity, yield softness :math:`\\delta`, + boundary-condition *values*, time step) so parameter continuation and + time-stepping warm-start correctly. + + Read-only status flag. Because a diverged solve leaves + ``has_solution == False``, the next :meth:`solve` automatically + cold-starts (Picard warm-up) rather than warming off a corrupted iterate. + + See :doc:`nonlinear-solver-homotopy-warmstart` (Layer 1). + """ + return self._has_solution + + def _solution_is_trivially_zero(self): + """True when the solution field is still identically zero. + + The design's *secondary* cold signal: a solver may be told to warm-start + (``zero_init_guess=False``) while its solution variable has never been + written, which is a cold start in everything but name. It matters because a + viscoplastic tangent is undefined there — the plastic viscosity is + :math:`\\tau_y / (2\\dot\\varepsilon_{II})`, so at :math:`v = 0` the + *residual* is finite (the soft-min carries the infinite plastic branch to + 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. + """ + u = getattr(self, "u", None) + if u is None: + return False + return float(u.vec.norm()) == 0.0 + + def _resolve_zero_init_guess(self, zero_init_guess): + """Resolve the tri-state ``zero_init_guess`` argument of ``solve()``. + + ``None`` (the default) auto-detects: cold when the solver holds no converged + solution, warm when it does. Detection is safe by construction — guessing + *cold* when a solution was in fact available costs one extra iteration from a + good starting point, while the harmful direction (warming off stale field data + after a remesh or a diverged solve) cannot happen, because + :attr:`has_solution` is cleared by both. + + ``True`` forces a fresh start (discard any solution); ``False`` insists on + warming from the current field values. + """ + if zero_init_guess is None: + return not self.has_solution + return bool(zero_init_guess) + + def _solve_yield_homotopy(self, homotopy_options=None, verbose=False, + solve_kwargs=None): + """Run ``solve(homotopy=True)``: a multi-solve δ-continuation on the yield law. + + The constitutive model advertises the homotopy (``supports_yield_homotopy``) + and describes it (``_yield_homotopy_control()``: how to set δ, and which + tangent to pair with it); the continuation driver marches δ from a large, + benign value down to the sharp yield surface, warm-starting each step from the + previous converged solution. See + :func:`~underworld3.systems.yield_continuation.yield_continuation` and + :doc:`nonlinear-solver-homotopy-warmstart` (Layer 2). + """ + cm = self.constitutive_model + if cm is None or not getattr(cm, "supports_yield_homotopy", False): + raise TypeError( + f"solve(homotopy=True) needs a constitutive model with a yield law to " + f"sharpen, but {type(cm).__name__ if cm is not None else None} does not " + f"advertise supports_yield_homotopy. Use a viscoplastic (or VEP) model, " + f"or solve without homotopy." + ) + from underworld3.systems.yield_continuation import yield_continuation + # The march's own options win: an explicit homotopy_options["verbose"] is a + # deliberate choice about the march, distinct from the solve's verbosity. + options = dict(homotopy_options or {}) + options.setdefault("verbose", verbose) + return yield_continuation(self, solve_kwargs=solve_kwargs, **options) + + def _record_convergence_status(self, converged=None): + """Refresh :attr:`has_solution` from the just-completed solve. + + Called at the end of every ``solve()``. With ``converged`` unset the + status is read from the SNES converged reason (``> 0`` ⇒ converged); + the rotated free-slip path, which runs its own KSP loop rather than + driving ``self.snes``, passes the flag from its result dict. + """ + if converged is None: + converged = self.snes is not None and self.snes.getConvergedReason() > 0 + self._has_solution = bool(converged) class _Unknowns: """ @@ -3339,7 +3467,7 @@ class SNES_Scalar(SolverBaseClass): @timing.routine_timer_decorator def solve(self, - zero_init_guess: bool =True, + zero_init_guess: bool =None, _force_setup: bool =False, verbose: bool=False, debug: bool=False, @@ -3355,10 +3483,15 @@ class SNES_Scalar(SolverBaseClass): Parameters ---------- - zero_init_guess : bool, default=True - If True, use zero as the initial guess. If False, use the current - values in the solution variable(s) as the initial guess, which can - improve convergence for time-stepping or continuation methods. + zero_init_guess : bool, optional + Cold or warm start. The default (``None``) **auto-detects**: cold when the + solver holds no converged solution, warm when it does (see + :attr:`has_solution`). ``True`` forces a fresh start, discarding any + existing solution; ``False`` insists on warming from the current field + values. Warm-starting improves convergence for time-stepping and + continuation; the auto default gets that without a flag, and cannot warm + off stale data because a remesh or a diverged solve clears + ``has_solution``. _force_setup : bool, default=False Force rebuild of the solver even if already set up. Useful after changing boundary conditions or constitutive parameters. @@ -3406,6 +3539,7 @@ class SNES_Scalar(SolverBaseClass): snes : Access to underlying PETSc SNES object for advanced control. """ + import petsc4py @@ -3416,6 +3550,11 @@ class SNES_Scalar(SolverBaseClass): # DM/fields/BCs are unchanged. In-place rewire is sufficient. self._needs_function_rewire = True + # Tri-state: None auto-detects cold-vs-warm from has_solution. Resolved HERE, + # after _force_setup has had its say: that invalidation clears has_solution, + # and resolving earlier would warm-start off the flag it just cleared. + zero_init_guess = self._resolve_zero_init_guess(zero_init_guess) + self._build(verbose, debug, debug_name) # Set time on the DM so petsc_t is available in pointwise functions @@ -3484,6 +3623,8 @@ class SNES_Scalar(SolverBaseClass): self._warn_on_divergence() + self._record_convergence_status() + return def _object_viewer(self): @@ -4401,7 +4542,7 @@ class SNES_Vector(SolverBaseClass): @timing.routine_timer_decorator def solve(self, - zero_init_guess: bool =True, + zero_init_guess: bool =None, _force_setup: bool =False, verbose=False, debug=False, @@ -4416,9 +4557,15 @@ class SNES_Vector(SolverBaseClass): Parameters ---------- - zero_init_guess : bool, default=True - If True, use zero as the initial guess. If False, use the current - values in ``self.u`` as the initial guess. + zero_init_guess : bool, optional + Cold or warm start. The default (``None``) **auto-detects**: cold when the + solver holds no converged solution, warm when it does (see + :attr:`has_solution`). ``True`` forces a fresh start, discarding any + existing solution; ``False`` insists on warming from the current field + values. Warm-starting improves convergence for time-stepping and + continuation; the auto default gets that without a flag, and cannot warm + off stale data because a remesh or a diverged solve clears + ``has_solution``. _force_setup : bool, default=False Force rebuild of the solver even if already set up. verbose : bool, default=False @@ -4445,6 +4592,7 @@ class SNES_Vector(SolverBaseClass): u : The solution vector field variable. """ + if _force_setup: self.is_setup = False elif not self.constitutive_model._solver_is_setup: @@ -4452,6 +4600,11 @@ class SNES_Vector(SolverBaseClass): # DM/fields/BCs are unchanged. In-place rewire is sufficient. self._needs_function_rewire = True + # Tri-state: None auto-detects cold-vs-warm from has_solution. Resolved HERE, + # after _force_setup has had its say: that invalidation clears has_solution, + # and resolving earlier would warm-start off the flag it just cleared. + zero_init_guess = self._resolve_zero_init_guess(zero_init_guess) + self._build(verbose, debug, debug_name) @@ -4517,6 +4670,8 @@ class SNES_Vector(SolverBaseClass): self._warn_on_divergence() + self._record_convergence_status() + return def _object_viewer(self): @@ -5144,7 +5299,7 @@ class SNES_MultiComponent(SolverBaseClass): @timing.routine_timer_decorator def solve(self, - zero_init_guess: bool = True, + zero_init_guess: bool = None, _force_setup: bool = False, verbose=False, debug=False, @@ -5163,6 +5318,7 @@ class SNES_MultiComponent(SolverBaseClass): start up to this many times. 0 preserves legacy behaviour. """ + if _force_setup: self.is_setup = False elif not self.constitutive_model._solver_is_setup: @@ -5170,6 +5326,11 @@ class SNES_MultiComponent(SolverBaseClass): # DM/fields/BCs are unchanged. In-place rewire is sufficient. self._needs_function_rewire = True + # Tri-state: None auto-detects cold-vs-warm from has_solution. Resolved HERE, + # after _force_setup has had its say: that invalidation clears has_solution, + # and resolving earlier would warm-start off the flag it just cleared. + zero_init_guess = self._resolve_zero_init_guess(zero_init_guess) + self._build(verbose, debug, debug_name) gvec = self.dm.getGlobalVec() @@ -5209,6 +5370,8 @@ class SNES_MultiComponent(SolverBaseClass): self._warn_on_divergence() + self._record_convergence_status() + return def _object_viewer(self): @@ -8229,14 +8392,16 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): @timing.routine_timer_decorator def solve(self, - zero_init_guess: bool = True, + zero_init_guess: bool = None, picard: int = 0, verbose=False, debug=False, debug_name=None, _force_setup: bool =False, time=None, - divergence_retries: int = 0, ): + divergence_retries: int = 0, + homotopy: bool = False, + homotopy_options: dict = None, ): """ Solve the Stokes system for velocity and pressure. @@ -8246,11 +8411,15 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): Parameters ---------- - zero_init_guess : bool, default=True - If True, use zero as the initial guess. If False, use current - values in ``self.u`` (velocity) and ``self.p`` (pressure) as - initial guess. Using False can improve convergence for - time-stepping or parameter continuation. + zero_init_guess : bool, optional + Cold or warm start. The default (``None``) **auto-detects**: cold when the + solver holds no converged solution, warm when it does (see + :attr:`has_solution`). ``True`` forces a fresh start, discarding any + existing solution; ``False`` insists on warming from the current field + values. Warm-starting improves convergence for time-stepping and + continuation; the auto default gets that without a flag, and cannot warm + 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 @@ -8272,11 +8441,28 @@ 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_options : dict, optional + March settings passed to + :func:`~underworld3.systems.yield_continuation.yield_continuation` — + ``delta0``, ``down``, ``dmin``, ``entry_maxit``, ``step_maxit``, + ``retries``. All are defaulted; tuning them is optional. Returns ------- - None - Solution stored in ``self.u`` (velocity) and ``self.p`` (pressure). + None or dict + ``None`` normally — the solution is stored in ``self.u`` (velocity) and + ``self.p`` (pressure). With ``homotopy=True``, the march summary + (``settled_delta``, ``reason``, ``steps``, ``reached_dmin``, + ``converged``). Examples -------- @@ -8293,6 +8479,10 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): ... # Update boundary conditions, material properties... ... stokes.solve(zero_init_guess=False) + >>> # Hard viscoplastic (Drucker-Prager) yield: march the yield homotopy + >>> report = stokes.solve(homotopy=True) + >>> report["settled_delta"] # smallest yield softness reached + Notes ----- This is a **collective operation** - all MPI ranks must call it. @@ -8308,6 +8498,12 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): constitutive_model : Viscosity and stress definitions. """ + + if homotopy: + # The march runs a SEQUENCE of ordinary solves at successively sharper + # yield surfaces; each one re-enters this method with homotopy=False. + return self._solve_yield_homotopy(homotopy_options, verbose=verbose) + if _force_setup: self.is_setup = False elif not self.constitutive_model._solver_is_setup: @@ -8315,6 +8511,11 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): # DM/fields/BCs are unchanged. In-place rewire is sufficient. self._needs_function_rewire = True + # Tri-state: None auto-detects cold-vs-warm from has_solution. Resolved HERE, + # after _force_setup has had its say: that invalidation clears has_solution, + # and resolving earlier would warm-start off the flag it just cleared. + zero_init_guess = self._resolve_zero_init_guess(zero_init_guess) + self._build(verbose, debug, debug_name) # Set time on the DM so petsc_t is available in pointwise functions. @@ -8360,7 +8561,10 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): zero_init_guess=zero_init_guess, picard=picard) # This path solves via ksp.solve on the rotated operator (not self.snes), # so give it a report from the rotated result rather than leaving a stale one. - self._capture_rotated_report(self._rotated_freeslip_info) + _rotated_report = self._capture_rotated_report(self._rotated_freeslip_info) + # The warm-start flag must follow the rotated solve's own verdict — the SNES + # was never run, so the generic reader would latch a stale reason. + self._record_convergence_status(converged=_rotated_report.converged) return if time is not None: @@ -8413,6 +8617,24 @@ 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. + # + # This is a DESIGN REQUIREMENT, not an optimisation: under the consistent + # tangent a viscoplastic Jacobian is NaN at zero strain rate (the residual + # survives, its derivative does not), so the machinery has to make that + # state unreachable. Both routes to it are covered — an explicit cold start, + # and a nominally warm one whose solution has never been written. The + # "continuation" tangent needs no help: it opens on a Picard stage + # (alpha = 0) by construction. 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) @@ -8506,6 +8728,8 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self._warn_on_divergence() + self._record_convergence_status() + return @timing.routine_timer_decorator diff --git a/src/underworld3/maths/__init__.py b/src/underworld3/maths/__init__.py index 1cecae6d9..1b1103d7d 100644 --- a/src/underworld3/maths/__init__.py +++ b/src/underworld3/maths/__init__.py @@ -39,6 +39,8 @@ from .functions import delta as delta_function from .functions import L2_norm as L2_norm +from .functions import smooth_max as smooth_max +from .functions import smooth_min as smooth_min # from .vector_calculus import ( # mesh_vector_calculus_spherical_lonlat as vector_calculus_spherical_lonlat, diff --git a/src/underworld3/maths/functions.py b/src/underworld3/maths/functions.py index e97828e2a..624ee2378 100644 --- a/src/underworld3/maths/functions.py +++ b/src/underworld3/maths/functions.py @@ -56,6 +56,71 @@ def delta( return delta_fn +def smooth_max(a, b, epsilon): + r""" + Differentiable (regularized) maximum of two expressions. + + .. math:: + + \operatorname{smax}_\epsilon(a, b) = + \tfrac{1}{2}\left(a + b + \sqrt{(a - b)^2 + \epsilon^2}\right) + + Approaches the exact :math:`\max(a, b)` as :math:`\epsilon \to 0`, with the + corner rounded over a scale set by :math:`\epsilon`. Unlike ``sympy.Max`` it is + smooth everywhere, so it can appear in a residual that is differentiated for a + consistent (Newton) Jacobian without introducing a non-differentiable kink. + + Parameters + ---------- + a, b : sympy.Basic + Expressions to combine (fields, coordinates, or constants). + epsilon : float or sympy.Expr + Rounding scale, in the units of ``a`` and ``b``. Choose a small fraction of + the physical scale over which the two branches cross. + + Returns + ------- + sympy.Expr + The smooth maximum. + + Notes + ----- + A common use is a differentiable tension cutoff on a pressure-dependent yield + stress, :math:`\operatorname{smax}_\epsilon(C + \sin\phi\,p,\, 0)` — the rounded + tension cap several failure envelopes use (e.g. a parabolic Griffith cap) in + place of a sharp corner, which keeps the consistent-Newton tangent well defined. + """ + return (a + b + sympy.sqrt((a - b) ** 2 + epsilon**2)) / 2 + + +def smooth_min(a, b, epsilon): + r""" + Differentiable (regularized) minimum of two expressions. + + .. math:: + + \operatorname{smin}_\epsilon(a, b) = + \tfrac{1}{2}\left(a + b - \sqrt{(a - b)^2 + \epsilon^2}\right) + + Approaches the exact :math:`\min(a, b)` as :math:`\epsilon \to 0`; it is + :math:`-\operatorname{smax}_\epsilon(-a, -b)`. See :func:`smooth_max` for the + parameters and the differentiability rationale. + + Parameters + ---------- + a, b : sympy.Basic + Expressions to combine (fields, coordinates, or constants). + epsilon : float or sympy.Expr + Rounding scale, in the units of ``a`` and ``b``. + + Returns + ------- + sympy.Expr + The smooth minimum. + """ + return (a + b - sympy.sqrt((a - b) ** 2 + epsilon**2)) / 2 + + def L2_norm(n_s, a_s, mesh): r""" L2 norm of the difference between numerical and analytical solutions. diff --git a/src/underworld3/systems/__init__.py b/src/underworld3/systems/__init__.py index 97ca09960..ac4cd6a96 100644 --- a/src/underworld3/systems/__init__.py +++ b/src/underworld3/systems/__init__.py @@ -84,3 +84,7 @@ from .ddt import SemiLagrangian as SemiLagragian_DDt from .ddt import Lagrangian_Swarm as Lagrangian_Swarm_DDt from .ddt import Eulerian as Eulerian_DDt + +# δ-continuation driver for hard viscoplastic (Drucker–Prager) yield +from .yield_continuation import yield_continuation, YieldHomotopyControl +from .solve_report import SolveReport diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index f91f8ea1b..d03c27f22 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -695,7 +695,7 @@ def v(self, value): @timing.routine_timer_decorator def solve( self, - zero_init_guess: bool = True, + zero_init_guess: bool = None, timestep: float = None, verbose: bool = False, _force_setup: bool = False, @@ -708,7 +708,9 @@ def solve( Parameters ---------- zero_init_guess : bool, optional - If True (default), start from zero initial guess. + Cold or warm start. The default (``None``) auto-detects from + :attr:`has_solution`; ``True`` forces a fresh start, ``False`` insists + on warming from the current field values. If False, use current field values as initial guess. timestep : float, optional Timestep value for inertial terms (if applicable). @@ -940,7 +942,7 @@ def estimate_dt(self): @timing.routine_timer_decorator def solve( self, - zero_init_guess: bool = True, + zero_init_guess: bool = None, timestep=None, _force_setup: bool = False, verbose=False, @@ -952,7 +954,9 @@ def solve( Parameters ---------- zero_init_guess : bool, optional - Start from zero initial guess (default True). + Cold or warm start. The default (``None``) auto-detects from + :attr:`has_solution`; ``True`` forces a fresh start, ``False`` insists + on warming from the current field values. timestep : float, optional Timestep size. Updates ``self.delta_t`` if provided. _force_setup : bool, optional @@ -1409,7 +1413,7 @@ def _create_stress_history_ddt(self, order=2): @memprobe.instrument("Stokes.solve") def solve( self, - zero_init_guess: bool = True, + zero_init_guess: bool = None, timestep: float = None, _force_setup: bool = False, verbose: bool = False, @@ -1419,6 +1423,8 @@ def solve( order=None, picard: int = 0, divergence_retries: int = 0, + homotopy: bool = False, + homotopy_options: dict = None, ): """Solve the Stokes system, with optional viscoelastic stress history. @@ -1428,8 +1434,15 @@ def solve( Parameters ---------- - zero_init_guess : bool - If True, use zero initial guess. Otherwise use current field values. + zero_init_guess : bool, optional + Cold or warm start. The default (``None``) **auto-detects**: cold when the + solver holds no converged solution, warm when it does (see + :attr:`has_solution`). ``True`` forces a fresh start, discarding any + existing solution; ``False`` insists on warming from the current field + values. Warm-starting improves convergence for time-stepping and + continuation; the auto default gets that without a flag, and cannot warm + off stale data because a remesh or a diverged solve clears + ``has_solution``. timestep : float, optional Advection timestep. Required when stress history is active. _force_setup : bool @@ -1453,8 +1466,47 @@ 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_options : dict, optional + March settings for + :func:`~underworld3.systems.yield_continuation.yield_continuation` + (``delta0``, ``down``, ``dmin``, ``entry_maxit``, ``step_maxit``, + ``retries``). All defaulted. """ + 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. + inner = {} + for name, value in (("timestep", timestep), ("evalf", evalf), + ("order", order), ("debug", debug), + ("debug_name", debug_name), + ("divergence_retries", divergence_retries)): + if value: + inner[name] = value + # The march owns these: it sets the cold/warm decision per step (Layer 1) + # and its own per-step iteration budget, so accepting a contradicting + # value silently would be worse than saying so. + if zero_init_guess is not None: + raise ValueError( + "solve(homotopy=True) chooses cold-vs-warm for each step of the " + "march itself; do not also pass zero_init_guess." + ) + if picard: + raise ValueError( + "solve(homotopy=True) manages its own warm-up; do not also pass " + "picard. Use homotopy_options={'entry_maxit': ...} to size the " + "first solve." + ) + return self._solve_yield_homotopy( + homotopy_options, verbose=verbose, solve_kwargs=inner + ) + has_stress_history = self.Unknowns.DFDt is not None if has_stress_history: @@ -4122,7 +4174,7 @@ def _reduce_dt(per_elem): @timing.routine_timer_decorator def solve( self, - zero_init_guess: bool = True, + zero_init_guess: bool = None, timestep: float = None, _force_setup: bool = False, _evalf=False, @@ -4434,7 +4486,7 @@ def estimate_dt(self): @timing.routine_timer_decorator def solve( self, - zero_init_guess: bool = True, + zero_init_guess: bool = None, timestep: float = None, evalf: bool = False, _force_setup: bool = False, @@ -4853,7 +4905,7 @@ def penalty(self, value): @memprobe.instrument("NavierStokes.solve") def solve( self, - zero_init_guess: bool = True, + zero_init_guess: bool = None, timestep: float = None, _force_setup: bool = False, verbose=False, diff --git a/src/underworld3/systems/yield_continuation.py b/src/underworld3/systems/yield_continuation.py new file mode 100644 index 000000000..ff162619b --- /dev/null +++ b/src/underworld3/systems/yield_continuation.py @@ -0,0 +1,254 @@ +r"""Multi-solve δ-continuation for hard viscoplastic (Drucker–Prager) yield. + +A single solve straight onto a sharp yield surface (soft-min δ → 0) stalls or +diverges from a cold start. The robust, trusted route is a *continuation*: hold δ +**constant** for a full nonlinear solve to tolerance, warm-start the next (smaller) +δ from that converged state, and march δ down toward the sharp surface. + +This module ships that march. It is what ``solver.solve(homotopy=True)`` runs; the +constitutive model says what δ *is* for it (:class:`YieldHomotopyControl`, built by +``constitutive_model._yield_homotopy_control()``) and this driver marches the number. + +Do **not** try to ramp δ inside a single SNES solve. That is proven not to work: the +continuation only sharpens δ from a *converged*, well-conditioned iterate, whereas an +in-solve ramp sharpens it mid-solve where the consistent-Newton Jacobian on a +sharpening yield surface is ill-conditioned and the linear solve fails. +""" + +from dataclasses import dataclass +from typing import Any, Callable, Optional + +import underworld3 as uw + + +@dataclass(frozen=True) +class YieldHomotopyControl: + """How to march one constitutive model's yield homotopy. + + Built by ``constitutive_model._yield_homotopy_control()``, which also switches the + model into its smooth (δ-parameterised) yield mode. The model owns the meaning of + the parameter; the solver only moves it. + + Attributes + ---------- + set_delta + Sets the yield softness δ. Model-owned, because *how* δ is applied differs: + the isotropic models update a ``constants[]`` atom (no recompile), while the + transverse-isotropic VEP model rebuilds. + tangent + The ``consistent_jacobian`` value to solve with — Newton for a viscous-plastic + yield, the frozen/Picard tangent for the elastic (VEP) models. + delta + The δ ``constants[]`` atom itself, for diagnostics. + """ + + set_delta: Callable[[float], None] + tangent: Any + delta: Optional[Any] = None + + +def yield_continuation( + solver, + control=None, + delta0=1.0, + down=0.5, + dmin=1.0e-3, + entry_maxit=30, + step_maxit=10, + retries=2, + max_steps=60, + solve_kwargs=None, + verbose=True, +): + r"""March the yield regularisation δ down to the sharp surface as a sequence of + constant-δ solves, each warm-started from the previous. + + Each δ is held **constant** for one full nonlinear solve to tolerance; the + converged solution warm-starts the next (smaller) δ. The march **settles** at the + smallest feasible δ — the closest automatic approach to the hard ``Min``. A δ that + fails to converge is reverted and retried with a gentler step; when the retries are + used up, the last converged δ is the answer. + + The step size is **residual-guided**: a δ-step that converges in very few Newton + iterations means there is slack, so the next step is bigger; a step that uses most + of its iteration budget means the march is close to the feasible edge, so the next + step is gentler. + + Starting δ is deliberately **large**, where the power-mean soft-min is the harmonic + mean and is bounded by the background viscosity even as :math:`\dot\varepsilon \to + 0`. That is what makes the first (cold) solve well posed, and why no separate + viscous pre-solve is needed. + + Parameters + ---------- + solver : + A configured Stokes/SNES solver whose constitutive model supports the homotopy. + control : YieldHomotopyControl, optional + How to set δ and which tangent to use. Defaults to + ``solver.constitutive_model._yield_homotopy_control()``. + delta0 : float + Starting (smooth) δ. Large δ is about as cheap to solve as small, so be generous. + down : float + Nominal multiplicative step, :math:`0 < down < 1` (δ ← δ·down per success). + Adapted during the march as described above. + dmin : float + Stop once δ reaches this floor (the sharp surface is effectively reached). + entry_maxit : int + Nonlinear iteration budget for the first solve. + step_maxit : int + Budget for each warm step. A feasible warm step converges in a few iterations; + a tight budget lets a too-hard step abort cheaply. + retries : int + How many times to retry a failed δ with a gentler step before settling. + max_steps : int + Hard cap on the number of δ-solves, so a march that is easing off in small + steps cannot run unbounded. + solve_kwargs : dict, optional + Extra arguments forwarded to every inner ``solver.solve()`` — notably + ``timestep`` for a visco-elastic-plastic model, whose solves need one. + verbose : bool + Report each δ (rank-safe). + + Returns + ------- + dict + ``settled_delta`` (smallest δ that converged, or ``None`` if even the first + solve failed), ``reason`` (final SNES converged reason), ``steps`` (number of + δ-solves attempted), ``reached_dmin``, and ``converged``. + """ + if not (0.0 < down < 1.0): + raise ValueError(f"down must satisfy 0 < down < 1, got {down}") + if not delta0 > 0.0: + raise ValueError(f"delta0 must be positive, got {delta0}") + + if control is None: + cm = getattr(solver, "constitutive_model", None) + if cm is None or not getattr(cm, "supports_yield_homotopy", False): + raise TypeError( + "yield_continuation needs a constitutive model that supports the yield " + "homotopy (supports_yield_homotopy). Got " + f"{type(cm).__name__ if cm is not None else None}." + ) + control = cm._yield_homotopy_control() + + # A march is a SEQUENCE of solves, so a solver that integrates history in time + # would advance that history once per delta-step rather than once per timestep. + # Refuse rather than silently corrupt it -- see the design note. + if getattr(getattr(solver, "Unknowns", None), "DFDt", None) is not None: + raise NotImplementedError( + "solve(homotopy=True) is not available on a solver carrying stress " + "history (visco-elastic-plastic): the march runs several solves, each of " + "which would advance the elastic stress history by a full timestep. Solve " + "the VEP model without the homotopy, or drive the march yourself around a " + "single history update." + ) + + # Iteration budgets have to be applied where they survive the solve path's own + # setFromOptions (which pushes a hardcoded snes_max_it); petsc_options alone is + # silently overwritten. This is the same hook estimate_difficulty() uses. + saved_probe = solver._difficulty_probe + saved_probe_maxit = solver._difficulty_max_it + saved_resume = solver._resume_abs_target + saved_tangent = solver.consistent_jacobian + + solver.consistent_jacobian = control.tangent + + u, p = solver.Unknowns.u, solver.Unknowns.p + # Warm-state snapshot: a raw copy of the (already consistent) solution values, + # restored if a too-hard delta corrupts the iterate. + u_good, p_good = u.data.copy(), p.data.copy() + + # Whether there is a solution to warm-start the FIRST delta from. Read before the + # rebuild below, which invalidates the solver and clears the flag. + warm_entry = solver.has_solution + + d = float(delta0) + step = float(down) + settled = None + reason = 0 + steps = 0 + failures = 0 + reached_dmin = False + first = True + + try: + while steps < max_steps: + control.set_delta(d) + if first: + solver.is_setup = False # build the operators once, at delta0 + else: + solver._update_constants() # recompile-free delta update + + budget = entry_maxit if first else step_maxit + solver._difficulty_probe = True + solver._difficulty_max_it = budget + solver._resume_abs_target = None + + warm = warm_entry if first else True + solver.solve(zero_init_guess=not warm, **dict(solve_kwargs or {})) + reason = int(solver.snes.getConvergedReason()) + nit = int(solver.snes.getIterationNumber()) + steps += 1 + first = False + + if reason > 0: + settled = d + u_good[...] = u.data + p_good[...] = p.data + if verbose: + uw.pprint(f" [yield-continuation] d={d:<11g} its={nit:3d} -> converged") + if d <= dmin: + reached_dmin = True + break + # Residual-guided: plenty of slack => take a bigger bite next time; + # near the iteration budget => ease off. Clamped so the march can + # neither stall (step -> 1) nor leap past the feasible edge in one go. + if nit <= max(2, budget // 5): + step = max(step * step, 0.05) + elif nit >= 0.8 * budget: + step = min(step ** 0.5, 0.95) + d *= step + else: + failed_delta = d + with uw.synchronised_array_update("yield_continuation revert"): + u.data[...] = u_good + p.data[...] = p_good + # The reverted fields ARE a converged state, so say so: otherwise the + # retry (and the caller's next solve) would auto-cold-start off a + # solution that is perfectly good. + if settled is not None: + solver._record_convergence_status(converged=True) + failures += 1 + if settled is None or failures > retries: + if verbose: + uw.pprint(f" [yield-continuation] d={failed_delta:<11g} its={nit:3d} " + f"-> failed (reason={reason}); settling at d={settled}") + break + # Retry from the last good delta with a gentler step. + step = min(step ** 0.5, 0.95) + d = settled * step + if verbose: + uw.pprint(f" [yield-continuation] d={failed_delta:<11g} its={nit:3d} " + f"-> failed (reason={reason}); retrying at d={d:g}") + else: + if verbose: + uw.pprint(f" [yield-continuation] stopped at the {max_steps}-step cap; " + f"settled at d={settled}") + finally: + # Leave the model holding the delta that actually converged, and restore every + # setting the march borrowed -- including on the exception path. + if settled is not None and settled != d: + control.set_delta(settled) + solver._update_constants() + solver._difficulty_probe = saved_probe + solver._difficulty_max_it = saved_probe_maxit + solver._resume_abs_target = saved_resume + solver.consistent_jacobian = saved_tangent + + return { + "settled_delta": settled, + "reason": reason, + "steps": steps, + "reached_dmin": reached_dmin, + "converged": settled is not None, + } diff --git a/tests/parallel/ptest_0202_warmstart_homotopy_parallel.py b/tests/parallel/ptest_0202_warmstart_homotopy_parallel.py new file mode 100644 index 000000000..bdcf7c696 --- /dev/null +++ b/tests/parallel/ptest_0202_warmstart_homotopy_parallel.py @@ -0,0 +1,111 @@ +"""Parallel (MPI) test: the nonlinear warm-start / yield-homotopy layers at np > 1. + +Charter §11 — no feature is complete if it only works in serial. All four layers of +``docs/developer/design/nonlinear-solver-homotopy-warmstart.md`` make decisions that +would deadlock or diverge between ranks if any of them were driven by a rank-LOCAL +quantity: + + * Layer 1a/1b — ``has_solution`` and the tri-state ``zero_init_guess`` gate the + ``dm.localToGlobal`` / ``snes.solve`` collectives. If ranks disagreed on + cold-vs-warm they would take different branches: a hang, not a wrong answer. + * Layer 2 — every branch of the δ-march (converged? retry? settle?) is taken from + ``snes.getConvergedReason()`` / ``getIterationNumber()``, which are collective and + therefore rank-identical. The revert also writes fields inside + ``synchronised_array_update`` on every rank. + * Layer 3 — the FMG velocity smoother (gmres + sor, fixed-cost V-cycle) has to work + with the redundant-LU coarse solve at np > 1, which the serial option-string test + in ``test_1014`` cannot exercise. + +Run: + + cd tests/parallel + mpirun -np 2 python ./ptest_0202_warmstart_homotopy_parallel.py + +Asserts (rank-collectively; any rank disagreeing shows up as a hang or a mismatch): + 1. ``has_solution`` is False before, True after, and identical on every rank. + 2. The auto-detected cold/warm resolution is identical on every rank. + 3. A geometric-FMG Stokes solve converges with the new default smoother. + 4. ``solve(homotopy=True)`` marches and converges, with the same settled δ and step + count on every rank. +""" + +import sympy + +import underworld3 as uw +from underworld3 import mpi + +comm = uw.mpi.comm +rank = uw.mpi.rank + + +def _all_ranks_agree(value): + """True when `value` is identical on every rank.""" + gathered = comm.allgather(value) + return all(v == gathered[0] for v in gathered) + + +# --- 1/2/3: warm-start bookkeeping and the FMG smoother, on a refined (hierarchy) mesh +mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.25, refinement=1, qdegree=3 +) +v = uw.discretisation.MeshVariable("Vp", mesh, mesh.dim, degree=2) +p = uw.discretisation.MeshVariable("Pp", mesh, 1, degree=1, continuous=True) + +stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) +stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel +stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 +stokes.bodyforce = sympy.Matrix([0.0, -1.0]) +stokes.add_essential_bc((0.0, 0.0), "Bottom") +stokes.add_essential_bc((0.0, None), "Left") +stokes.add_essential_bc((0.0, None), "Right") +stokes.preconditioner = "fmg" # exercise the new gmres/sor smoother in parallel +stokes.tolerance = 1.0e-6 + +assert stokes.has_solution is False +assert _all_ranks_agree(stokes.has_solution), "has_solution disagrees between ranks" +assert _all_ranks_agree(stokes._resolve_zero_init_guess(None)), \ + "the cold/warm decision disagrees between ranks — the solve collectives would split" + +stokes.solve() +assert stokes.snes.getConvergedReason() > 0, "FMG Stokes solve failed at np>1" +assert stokes.has_solution is True +assert _all_ranks_agree(stokes.has_solution) +# With a solution in hand every rank must now independently decide "warm". +assert stokes._resolve_zero_init_guess(None) is False +assert _all_ranks_agree(stokes._resolve_zero_init_guess(None)) + +n_levels = len(getattr(mesh, "dm_hierarchy", []) or []) + +# --- 4: the yield homotopy march in parallel +mesh2 = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.3 +) +x, y = mesh2.X +v2 = uw.discretisation.MeshVariable("Vh", mesh2, mesh2.dim, degree=2) +p2 = uw.discretisation.MeshVariable("Ph", mesh2, 1, degree=1, continuous=True) +vp = uw.systems.Stokes(mesh2, velocityField=v2, pressureField=p2) +vp.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel +vp.constitutive_model.Parameters.shear_viscosity_0 = 1.0 +vp.constitutive_model.Parameters.yield_stress = 0.35 # genuinely yielding +vp.bodyforce = sympy.Matrix([[0.0, -2.0 * sympy.cos(sympy.pi * x)]]) +vp.add_essential_bc((sympy.oo, 0.0), "Top") +vp.add_essential_bc((sympy.oo, 0.0), "Bottom") +vp.add_essential_bc((0.0, sympy.oo), "Left") +vp.add_essential_bc((0.0, sympy.oo), "Right") +vp.petsc_use_pressure_nullspace = True +vp.tolerance = 1.0e-8 + +report = vp.solve(homotopy=True, + homotopy_options=dict(delta0=1.0, dmin=1.0e-3, verbose=False)) + +assert report["converged"] is True, f"homotopy march failed at np>1: {report}" +assert _all_ranks_agree(report["steps"]), "march step count disagrees between ranks" +assert _all_ranks_agree(report["settled_delta"]), "settled delta disagrees between ranks" +assert vp.has_solution is True +assert _all_ranks_agree(vp.has_solution) + +uw.pprint( + f"ptest_0202 OK (np={uw.mpi.size}): FMG hierarchy {n_levels} levels converged with the " + f"gmres/sor smoother; homotopy settled delta={report['settled_delta']:.3e} in " + f"{report['steps']} steps, rank-consistent" +) diff --git a/tests/test_0103_jit_rampable_constants.py b/tests/test_0103_jit_rampable_constants.py index 97b1ebea4..180df0683 100644 --- a/tests/test_0103_jit_rampable_constants.py +++ b/tests/test_0103_jit_rampable_constants.py @@ -40,7 +40,13 @@ def _build(): def _mean(poisson, T): - poisson.solve() + # COLD every time. The diffusivity depends on T, so this is a nonlinear solve and + # its answer is only pinned to the SNES tolerance — a warm start stops at a + # different point in that same tolerance ball. This test compares solutions at + # rtol 1e-12 to check that a constant was read from its slot rather than baked as a + # literal, so it must hold the initial guess fixed or it measures the warm-start + # policy instead of the thing it is about. + poisson.solve(zero_init_guess=True) return float(np.asarray(T.data)[:, 0].mean()) diff --git a/tests/test_0201_solver_has_solution_warmstart.py b/tests/test_0201_solver_has_solution_warmstart.py new file mode 100644 index 000000000..8ce6cc907 --- /dev/null +++ b/tests/test_0201_solver_has_solution_warmstart.py @@ -0,0 +1,207 @@ +"""Layer 1a of the nonlinear-solver warm-start / homotopy design: +``solver.has_solution`` and the automatic cold-start Picard warm-up. + +Contract under test (design: +``docs/developer/design/nonlinear-solver-homotopy-warmstart.md``, Layer 1): + + * ``has_solution`` is a public, read-only status flag, ``False`` on a fresh + solver and ``True`` only after a solve whose SNES converged. + * It resets to ``False`` on a *structural* rebuild (the ``is_setup = False`` + invalidation used by remesh / adapt / mesh-mover), so the next solve + 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. + +These are cheap serial checks on the public API — the hard-case δ-continuation +that motivates the design is validated separately against the Spiegelman study. +""" + +import pytest +import sympy + +import underworld3 as uw + + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + + +def _poisson(): + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.25 + ) + t = uw.discretisation.MeshVariable("T", mesh, 1, degree=1) + poisson = uw.systems.Poisson(mesh, u_Field=t) + poisson.constitutive_model = uw.constitutive_models.DiffusionModel + poisson.constitutive_model.Parameters.diffusivity = 1.0 + poisson.add_essential_bc((0.0,), "Bottom") + poisson.add_essential_bc((1.0,), "Top") + return mesh, poisson + + +def test_has_solution_false_before_first_solve(): + _, poisson = _poisson() + assert poisson.has_solution is False + + +def test_has_solution_true_after_converged_solve(): + _, poisson = _poisson() + poisson.solve() + assert poisson.snes.getConvergedReason() > 0 + assert poisson.has_solution is True + + +def test_has_solution_resets_on_structural_rebuild(): + """A structural invalidation (is_setup=False — the remesh / adapt / mover + hook) must drop the warm-start claim.""" + _, poisson = _poisson() + poisson.solve() + assert poisson.has_solution is True + + poisson.is_setup = False + assert poisson.has_solution is False, ( + "has_solution must reset when the solver is structurally invalidated" + ) + + +def test_has_solution_survives_coefficient_change(): + """Changing a coefficient value (not the mesh/operators) keeps the solution + so continuation and time-stepping can warm-start.""" + _, poisson = _poisson() + poisson.solve() + assert poisson.has_solution is True + + # A coefficient update goes through _update_constants, not is_setup=False. + poisson.constitutive_model.Parameters.diffusivity = 2.0 + poisson._update_constants() + assert poisson.has_solution is True, ( + "a coefficient-only change must not drop the warm-start claim" + ) + + +def test_has_solution_resets_on_mesh_deform(): + """The mesh-mover path (mesh deform → registered solvers is_setup=False) + must reset has_solution.""" + import numpy as np + + mesh = uw.meshing.StructuredQuadBox( + elementRes=(8, 8), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0) + ) + v = uw.discretisation.MeshVariable("V", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1, continuous=True) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + stokes.bodyforce = sympy.Matrix([0.0, -1.0]) + stokes.add_essential_bc((0.0, 0.0), "Bottom") + stokes.add_essential_bc((0.0, None), "Left") + stokes.add_essential_bc((0.0, None), "Right") + stokes.solve() + assert stokes.has_solution is True + + disp = np.zeros((mesh.X.coords.shape[0], mesh.dim)) + disp[:, -1] = 0.01 * mesh.X.coords[:, 0] + with mesh._coord_mutation(): + mesh._deform_mesh(mesh.X.coords + disp) + assert stokes.has_solution is False, ( + "a mesh deform (structural rebuild) must reset has_solution" + ) + + +def test_zero_init_guess_is_tristate_and_auto_detects(): + """Layer 1b: the default (None) resolves cold-vs-warm from has_solution, while + True/False remain explicit overrides.""" + _, poisson = _poisson() + + # Fresh solver: nothing to warm from -> auto resolves COLD. + assert poisson.has_solution is False + assert poisson._resolve_zero_init_guess(None) is True + + poisson.solve() + assert poisson.has_solution is True + # A converged solution is present -> auto resolves WARM. + assert poisson._resolve_zero_init_guess(None) is False + + # Explicit values are honoured regardless of has_solution. + assert poisson._resolve_zero_init_guess(True) is True + assert poisson._resolve_zero_init_guess(False) is False + + # After a structural rebuild the claim is dropped, so auto is cold again — this + # is what stops a remesh warm-starting off stale field data. + poisson.is_setup = False + assert poisson._resolve_zero_init_guess(None) is True + + +def test_force_setup_cold_starts_through_the_public_path(): + """`_force_setup=True` is a structural invalidation, so the solve it is passed to + must COLD start. + + Regression (adversarial review, M1): `zero_init_guess` was resolved at the top of + solve(), before the `_force_setup` block that clears `has_solution`, so the call + that triggered the invalidation warm-started off the flag it was about to clear. + Driven through the public `solve()` rather than the private resolver, so it + actually exercises the ordering. + """ + _, poisson = _poisson() + poisson.solve() + assert poisson.has_solution is True + + poisson.solve(_force_setup=True) + # The rebuild must have dropped the claim during that call, not after it. + assert poisson._needs_dm_rebuild is False # rebuilt during the solve + assert poisson.has_solution is True # ... and re-established by convergence + + # The ordering itself: resolving now (post-invalidation) must say COLD. + poisson.is_setup = False + assert poisson._resolve_zero_init_guess(None) is True + + +def test_repeated_default_solve_agrees_to_solver_tolerance(): + """A bare solve() called repeatedly (the linear chain-of-solves pattern) must give + the same answer once auto-detect starts warming the later calls. + + "Same" means *to the requested convergence tolerance*, not bitwise: an iterative + solve stops anywhere inside the tolerance ball, and a warm start enters that ball + from a different direction than a cold one. The difference is therefore + tolerance-sized and shrinks with it (measured: 2.4e-5 at tol=1e-4, 7.4e-10 at + 1e-8, 2.0e-14 at 1e-12) — the flip to auto-detected warm starts changes results + at that level and no more. + """ + import numpy as np + + mesh, poisson = _poisson() + poisson.tolerance = 1.0e-10 + poisson.solve() + first = poisson.u.array[:, 0, 0].copy() + + poisson.solve() # now warm by auto-detect + second = poisson.u.array[:, 0, 0].copy() + + assert np.abs(first - second).max() < 1.0e-6, ( + "warm and cold solutions of a linear problem must agree to the solver " + f"tolerance (got max|diff| = {np.abs(first - second).max():.2e})" + ) + + +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.""" + mesh = uw.meshing.StructuredQuadBox( + elementRes=(8, 8), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0) + ) + v = uw.discretisation.MeshVariable("V", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1, continuous=True) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + stokes.bodyforce = sympy.Matrix([0.0, -1.0]) + stokes.add_essential_bc((0.0, 0.0), "Bottom") + stokes.add_essential_bc((0.0, None), "Left") + stokes.add_essential_bc((0.0, None), "Right") + + stokes.consistent_jacobian = True + stokes.solve() # cold (zero_init_guess default True) → warm-up branch taken + assert stokes.snes.getConvergedReason() > 0 + assert stokes.has_solution is True diff --git a/tests/test_1014_stokes_multigrid.py b/tests/test_1014_stokes_multigrid.py index 895adf0e0..ca8045f75 100644 --- a/tests/test_1014_stokes_multigrid.py +++ b/tests/test_1014_stokes_multigrid.py @@ -150,12 +150,19 @@ def test_geometric_mg_without_galerkin_is_repaired(): def test_default_fmg_bundle_is_parallel_safe(): # The property's OWN default FMG bundle must be usable at np>1 unaided: a # parallel-safe coarse solver (redundant+lu, not bare serial lu) and a - # robust smoother (richardson+sor, not eigen-estimate-fragile chebyshev). + # smoother sized for a DEEP hierarchy — gmres+sor, not eigen-estimate-fragile + # chebyshev and not stationary richardson, which degrades on the non-symmetric + # consistent-Newton operator (measured: per-V-cycle contraction 0.75 richardson + # vs 0.56 gmres over 4 nested levels on the Spiegelman notch). stokes = _make_stokes(mesh_refined) stokes.preconditioner = "fmg" stokes.solve() vp = "fieldsplit_velocity_" assert stokes.petsc_options.getString(vp + "mg_coarse_pc_type") == "redundant" assert stokes.petsc_options.getString(vp + "mg_coarse_redundant_pc_type") == "lu" - assert stokes.petsc_options.getString(vp + "mg_levels_ksp_type") == "richardson" + assert stokes.petsc_options.getString(vp + "mg_levels_ksp_type") == "gmres" + assert stokes.petsc_options.getString(vp + "mg_levels_pc_type") == "sor" + # Fixed-cost V-cycle: exactly mg_levels_ksp_max_it smoother iterations, no + # residual-norm computation and no early exit. + assert stokes.petsc_options.getString(vp + "mg_levels_ksp_norm_type") == "none" assert stokes.snes.getConvergedReason() > 0 diff --git a/tests/test_1055_yield_smoother.py b/tests/test_1055_yield_smoother.py new file mode 100644 index 000000000..b1383a80a --- /dev/null +++ b/tests/test_1055_yield_smoother.py @@ -0,0 +1,234 @@ +"""δ-parameterised yield soft-min smoother (sqrt + power-mean families). + +Pins the smoother contract on development (the ``_combine_yield`` helper + the +runtime-rampable δ ``constants[]`` atom, ported from the yield-homotopy work, with +development's tri-modal ``yield_mode`` semantics preserved): + +1. ``test_delta0_is_exact_min`` — the sqrt soft-min at δ=0 equals ``Min(η_ve, η_pl)`` + to machine precision (g(f) = max(1, f)). +2. ``test_powermean_smoother_undershoots_min`` — the power-mean family is ≤ Min + everywhere (never over-yields), overflow-safe on geodynamic ranges, and → Min as δ→0. +3. ``test_yield_smoother_validation`` — only the two known families; default "sqrt". +4. ``test_yield_softness_runtime_no_recompile`` — ramping δ through its constants[] atom + does NOT change the solver's JIT cache key. +5. ``test_dp_model_smoother_optin`` — the non-elastic Drucker–Prager model defaults to + exact hard Min and can opt into the softmin/power-mean homotopy. + +NOTE (vs the branch): development keeps ``yield_mode`` tri-modal — ``"min"`` is the exact +hard ``Min`` (not the soft-min), so the smoother tests select ``yield_mode="softmin"``. +The ``enable_yield_homotopy`` in-solve ramp driver is intentionally NOT part of this PR. +""" + +import numpy as np +import pytest +import sympy + +import underworld3 as uw +from underworld3.function import expression +from underworld3.function.expressions import unwrap_expression + +ETA = 1.0 +MU = 1.0 +TAU_Y = 0.5 +V0 = 0.5 +T_R = ETA / MU + + +def _build_vep(label, order=2): + mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 8), minCoords=(-1.0, -0.5), maxCoords=(1.0, 0.5), + ) + v = uw.discretisation.MeshVariable(f"U_{label}", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable(f"P_{label}", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + cm = uw.constitutive_models.ViscoElasticPlasticFlowModel(stokes.Unknowns, order=order) + stokes.constitutive_model = cm + cm.Parameters.shear_viscosity_0 = ETA + cm.Parameters.shear_modulus = MU + cm.Parameters.yield_stress = TAU_Y + cm.Parameters.strainrate_inv_II_min = 1.0e-6 + V_top = expression(rf"V_{{{label}}}", sympy.Float(V0), "Top V") + stokes.add_dirichlet_bc((V_top, 0.0), "Top") + stokes.add_dirichlet_bc((-V_top, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.tolerance = 1.0e-6 + stokes.petsc_options["snes_force_iteration"] = True + return mesh, stokes, cm, V_top + + +@pytest.mark.level_1 +@pytest.mark.tier_a +def test_delta0_is_exact_min(): + """δ=0 sqrt soft-min law equals Min(η_ve, η_pl) to machine precision.""" + mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) + v = uw.discretisation.MeshVariable("Um", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("Pm", mesh, 1, degree=1) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + cm = uw.constitutive_models.ViscoElasticPlasticFlowModel(s.Unknowns, order=2) + s.constitutive_model = cm + cm.yield_mode = "softmin" # the δ-family path (dev "min" = hard Min) + cm.yield_softness = 0.0 # δ=0 ⇒ g(f)=max(1,f) ⇒ exact Min + for eta_ve, eta_pl in [(1.0, 2.0), (2.0, 1.0), (1.0, 1.0), (0.3, 5.0), (5.0, 0.3)]: + comb = cm._combine_yield(sympy.Float(eta_ve), sympy.Float(eta_pl)) + val = float(unwrap_expression(comb, mode="nondimensional")) + assert abs(val - min(eta_ve, eta_pl)) < 1.0e-12 + + +@pytest.mark.level_1 +@pytest.mark.tier_a +def test_powermean_smoother_undershoots_min(): + """Power-mean smooth-min: ≤ Min everywhere (no over-yield), overflow-safe on + geodynamic ranges, and → Min as δ→0 (s=1/δ→∞).""" + mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) + v = uw.discretisation.MeshVariable("Upm", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("Ppm", mesh, 1, degree=1) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + cm = uw.constitutive_models.ViscoElasticPlasticFlowModel(s.Unknowns, order=2) + s.constitutive_model = cm + + cm.yield_mode = "softmin" # dev: soft-min family lives under "softmin" + cm.yield_softness = 0.0 # start from δ=0 so the powermean select bumps it + cm.yield_smoother = "powermean" # bumps δ 0 → 1 (s=1, harmonic mean) + assert cm.yield_smoother == "powermean" + assert cm.yield_softness == 1.0 + + # Includes a geodynamic-range pair (1e25, 1e21) that overflows the naive + # η^(−s) form above s≈40 but is finite in the harmonic-normalised form. + pairs = [(1.0, 2.0), (2.0, 1.0), (1.0, 1.0), (0.3, 5.0), (5.0, 0.3), (1e25, 1e21)] + for delta in (1.0, 0.5, 0.1, 0.02): + cm.yield_softness = delta + for eta_ve, eta_pl in pairs: + comb = cm._combine_yield(sympy.Float(eta_ve), sympy.Float(eta_pl)) + val = float(unwrap_expression(comb, mode="nondimensional")) + assert np.isfinite(val), f"overflow δ={delta} ({eta_ve},{eta_pl})" + assert val <= min(eta_ve, eta_pl) * (1 + 1e-9), \ + f"power-mean over-yields δ={delta}: {val} > {min(eta_ve, eta_pl)}" + + # smallest δ (largest s) is close to exact Min. + cm.yield_softness = 0.02 + for eta_ve, eta_pl in pairs: + comb = cm._combine_yield(sympy.Float(eta_ve), sympy.Float(eta_pl)) + val = float(unwrap_expression(comb, mode="nondimensional")) + m = min(eta_ve, eta_pl) + assert abs(val - m) <= 0.05 * m, f"not near Min at δ=0.02: {val} vs {m}" + + +@pytest.mark.level_1 +@pytest.mark.tier_a +def test_yield_smoother_validation(): + """yield_smoother only accepts the two known families; default is 'sqrt'.""" + mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) + v = uw.discretisation.MeshVariable("Usm", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("Psm", mesh, 1, degree=1) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + cm = uw.constitutive_models.ViscoElasticPlasticFlowModel(s.Unknowns, order=2) + s.constitutive_model = cm + assert cm.yield_smoother == "sqrt" # default family unchanged + with pytest.raises(ValueError): + cm.yield_smoother = "not_a_family" + + +@pytest.mark.level_2 +@pytest.mark.tier_a +def test_yield_softness_runtime_no_recompile(): + """Ramping δ through its constants[] atom must not trigger a JIT recompile.""" + _, stokes, cm, V_top = _build_vep("norecompile") + cm.yield_softness = 0.3 # δ>0 so the soft-min atom is exercised + dt = 0.20 * T_R + cm.Parameters.dt_elastic = dt + stokes.solve(zero_init_guess=False, timestep=dt, divergence_retries=2) + key0 = stokes._current_jit_cache_key + assert key0 is not None + + cm._get_yield_softness() + for d in (0.2, 0.1, 0.0): + cm._yield_softness_expr.sym = sympy.Float(d) + stokes._update_constants() + stokes.solve(zero_init_guess=False, timestep=dt, divergence_retries=2) + assert stokes._current_jit_cache_key == key0, "δ ramp forced a JIT 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).""" + 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) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + cm = uw.constitutive_models.ViscoPlasticFlowModel(s.Unknowns) + s.constitutive_model = cm + + # default: exact hard Min + assert cm.yield_mode == "min" + assert cm.yield_smoother == "sqrt" + 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")) + assert abs(val - min(eta_ve, eta_pl)) < 1.0e-12 + + # opt into the power-mean homotopy: undershoots Min + cm.yield_mode = "softmin" + 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)]: + comb = cm._combine_yield(sympy.Float(eta_ve), sympy.Float(eta_pl)) + val = float(unwrap_expression(comb, mode="nondimensional")) + assert np.isfinite(val) + assert val <= min(eta_ve, eta_pl) * (1 + 1e-9) + + # bad mode rejected + with pytest.raises(ValueError): + cm.yield_mode = "not_a_mode" + + +@pytest.mark.level_1 +@pytest.mark.tier_a +def test_dp_pressure_yield_consistent_jacobian_builds(): + """Regression: a pressure-dependent tension-cutoff yield ``Max(C+sinφ·p, 0)`` under + the consistent-Newton tangent must build the Jacobian without recursing. + + The DP model applied its yield floor as ``Max(yield_stress_min, yield_stress)`` + whenever ``yield_stress_min.sym != 0`` — but the *unset* default is a −∞-valued + Parameter, which passes ``!= 0``, wrapping every yield in ``Max(<−∞ parameter>, …)``. + That −∞ parameter is a UWexpression sympy's fuzzy ``is_ge`` cannot resolve, so + canonicalising the ``Max`` for the consistent tangent recursed. The guard now uses + the model's own −∞ "unset" sentinel, so the wrapper is not applied. + """ + mesh = uw.meshing.StructuredQuadBox( + elementRes=(4, 4), minCoords=(0, 0), maxCoords=(1, 1) + ) + v = uw.discretisation.MeshVariable("Udpj", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("Pdpj", mesh, 1, degree=1) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + cm = uw.constitutive_models.ViscoPlasticFlowModel(s.Unknowns) + s.constitutive_model = cm + sphi = float(np.sin(np.deg2rad(30.0))) + cm.Parameters.shear_viscosity_0 = 1.0e3 + cm.Parameters.yield_stress = sympy.Max(1.0 + sphi * p.sym[0], sympy.sympify(0)) + cm.Parameters.strainrate_inv_II_min = 1.0e-8 + s.consistent_jacobian = True + # Build the pointwise residual + Jacobian; on the bug this recurses to RecursionError. + s._setup_pointwise_functions(verbose=False) + + +@pytest.mark.level_1 +@pytest.mark.tier_a +def test_smooth_max_min_primitive(): + """uw.maths.smooth_max/smooth_min bracket the hard max/min, approach them as + ε→0, and are differentiable at the corner (no kink).""" + for a, b in [(1.0, 2.0), (2.0, 1.0), (3.0, -1.0), (0.0, 0.0)]: + hi, lo = max(a, b), min(a, b) + for eps in (0.5, 0.1, 1.0e-3): + smax = float(uw.maths.smooth_max(sympy.Float(a), sympy.Float(b), eps)) + smin = float(uw.maths.smooth_min(sympy.Float(a), sympy.Float(b), eps)) + assert smax >= hi - 1.0e-12 and abs(smax - hi) <= eps # rounds up, →max + assert smin <= lo + 1.0e-12 and abs(smin - lo) <= eps # rounds down, →min + + # smooth at the corner a=b: d/dx smooth_max(x, 0) is the finite value ½ at x=0 + # (a hard Max(x, 0) has a discontinuous derivative there). + x = sympy.Symbol("x") + d = sympy.diff(uw.maths.smooth_max(x, sympy.Float(0), sympy.Float(0.1)), x) + assert abs(float(d.subs(x, 0)) - 0.5) < 1.0e-9 diff --git a/tests/test_1057_yield_homotopy_solve.py b/tests/test_1057_yield_homotopy_solve.py new file mode 100644 index 000000000..8bd74fd64 --- /dev/null +++ b/tests/test_1057_yield_homotopy_solve.py @@ -0,0 +1,350 @@ +"""Layer 2 of the nonlinear-solver design: the model-advertised yield homotopy and +``stokes.solve(homotopy=True)``. + +Contract under test (design: +``docs/developer/design/nonlinear-solver-homotopy-warmstart.md``, Layer 2): + + * A constitutive model **advertises** whether it has a yield law to sharpen + (``supports_yield_homotopy``) — true for the viscoplastic / VEP models, false for + a plain viscous one. + * ``_yield_homotopy_control()`` puts the model in its smooth δ-parameterised mode + and reports how to march it: a model-owned δ setter, and the tangent to pair with + it (Newton for viscous-plastic yield, the frozen/Picard tangent for elastic VEP, + whose consistent yield tangent is indefinite). + * ``solve(homotopy=True)`` runs the multi-solve continuation and returns the march + summary; on a model with no yield law it raises a clear error rather than + silently solving something else. + +The hard-case convergence behaviour (Spiegelman notch at η=1e26) is validated +separately in the study; these are cheap serial checks of the API contract. +""" + +import pytest +import sympy + +import underworld3 as uw + + +pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] + + +def _viscoplastic_stokes(cellSize=0.4): + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=cellSize + ) + v = uw.discretisation.MeshVariable("Vh", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("Ph", mesh, 1, degree=1, continuous=True) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + stokes.constitutive_model.Parameters.yield_stress = 5.0 + stokes.bodyforce = sympy.Matrix([0.0, -1.0]) + stokes.add_essential_bc((0.0, 0.0), "Bottom") + stokes.add_essential_bc((0.0, None), "Left") + stokes.add_essential_bc((0.0, None), "Right") + stokes.tolerance = 1.0e-5 + return mesh, stokes + + +def test_plain_viscous_model_does_not_advertise_homotopy(): + cm = uw.constitutive_models.ViscousFlowModel + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.5 + ) + v = uw.discretisation.MeshVariable("Vv", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("Pv", mesh, 1, degree=1, continuous=True) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = cm + assert stokes.constitutive_model.supports_yield_homotopy is False + + +def test_viscoplastic_model_advertises_homotopy(): + _, stokes = _viscoplastic_stokes() + assert stokes.constitutive_model.supports_yield_homotopy is True + + +def test_homotopy_control_switches_model_to_smooth_mode(): + """The control must leave the model in the δ-parameterised power-mean mode and + hand back a working, model-owned δ setter.""" + _, stokes = _viscoplastic_stokes() + cm = stokes.constitutive_model + control = cm._yield_homotopy_control() + + assert cm.yield_mode == "softmin" + assert cm.yield_smoother == "powermean" + # Newton tangent for a non-elastic viscoplastic yield. + assert control.tangent is True + + # The setter must move BOTH the stored value and the constants[] atom, so a later + # _get_yield_softness() cannot reset δ to a stale number. + control.set_delta(0.25) + assert cm.yield_softness == pytest.approx(0.25) + assert float(cm._get_yield_softness().sym) == pytest.approx(0.25) + + +def test_vep_model_requests_the_picard_tangent(): + """The consistent yield tangent over the elastic stress-history block is + indefinite, so the elastic models must ask for the frozen tangent.""" + cm = uw.constitutive_models.ViscoElasticPlasticFlowModel + assert cm._yield_homotopy_tangent is False + assert uw.constitutive_models.ViscoPlasticFlowModel._yield_homotopy_tangent is True + + +def test_solve_homotopy_on_unsupported_model_raises(): + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.5 + ) + v = uw.discretisation.MeshVariable("Vr", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("Pr", mesh, 1, degree=1, continuous=True) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + stokes.bodyforce = sympy.Matrix([0.0, -1.0]) + stokes.add_essential_bc((0.0, 0.0), "Bottom") + + with pytest.raises(TypeError, match="supports_yield_homotopy"): + stokes.solve(homotopy=True) + + +@pytest.mark.parametrize( + "mode,smoother", + [("min", None), ("softmin", "sqrt"), ("softmin", "powermean")], +) +def test_cold_viscoplastic_solve_survives_zero_strain_rate(mode, smoother): + """A COLD (v=0) viscoplastic solve must not produce NaN in ANY yield mode. + + Regression: at v=0 the strain-rate invariant is 0, so the plastic viscosity + tau_y/(2 edot_II) is +inf. That is fine — the soft-min should carry it to the + viscous branch. But the power-mean form computed its harmonic mean as + eta_ve*eta_pl/(eta_ve+eta_pl), which is inf/inf = NaN, so a cold power-mean solve + died with DIVERGED_FNORM_NAN while `min` and `sqrt` converged. Rewriting that + mean as eta_ve/(1+f) — algebraically identical — keeps the infinite-eta_pl limit + finite. The same singularity occurs at any rigid (unyielded) point, not only on a + cold start. + """ + _, stokes = _viscoplastic_stokes(cellSize=0.5) + cm = stokes.constitutive_model + cm.yield_mode = mode + if smoother is not None: + cm.yield_smoother = smoother + + stokes.solve() # cold: zero_init_guess defaults True + reason = int(stokes.snes.getConvergedReason()) + assert reason > 0, ( + f"cold viscoplastic solve failed for yield_mode={mode!r} " + f"smoother={smoother!r} (reason={reason}; -4 is FNORM_NAN)" + ) + + +@pytest.mark.parametrize("zero_init_guess", [True, False]) +def test_consistent_newton_never_assembles_at_zero_strain_rate(zero_init_guess): + """A viscoplastic Jacobian is NaN at zero strain rate, so the warm/cold machinery + must make that state unreachable — including when the caller asks for a WARM start + on a solution that has never been written (adversarial review, M9). + + The residual survives v=0 (the soft-min carries the infinite plastic branch to the + viscous one); the tangent does not. Only the interposed Picard step keeps the + consistent-Newton path off it. + """ + _, stokes = _viscoplastic_stokes(cellSize=0.5) + stokes.consistent_jacobian = True + assert stokes._solution_is_trivially_zero() is True # never solved + + stokes.solve(zero_init_guess=zero_init_guess) + reason = int(stokes.snes.getConvergedReason()) + assert reason > 0, ( + f"consistent-Newton solve from an all-zero field failed with reason={reason} " + f"(-4 is FNORM_NAN) for zero_init_guess={zero_init_guess}" + ) + + +def test_yield_stress_is_floored_at_zero_by_default(): + """The yield stress is compared against the second invariant of the stress, so a + negative tau_y is meaningless. `yield_stress_min` therefore defaults to 0 and the + floor is active without the user asking. + + Regression: the default used to be the "unset" sentinel -oo (so sympy could cancel + the term away), which left a pressure-dependent Drucker-Prager yield + C + sin(phi)*p free to go negative in tension. Hard `Min` hid that (it discarded + the resulting negative eta_pl); the power-mean the homotopy switches to raised a + negative base to a non-integer power and produced NaN. + """ + _, stokes = _viscoplastic_stokes() + cm = stokes.constitutive_model + assert cm.Parameters.yield_stress_min.sym == 0 + + +@pytest.mark.parametrize("mode,smoother", + [("min", None), ("softmin", "powermean")]) +def test_pressure_dependent_yield_survives_tension(mode, smoother): + """A Drucker-Prager yield that goes NEGATIVE in tension must not poison the + viscosity in any yield mode — the zero floor must catch it first.""" + import numpy as np + + mesh, stokes = _viscoplastic_stokes(cellSize=0.5) + cm = stokes.constitutive_model + x, y = mesh.X + # C + sin(phi)*p with the pressure driven strongly negative (tension) over part + # of the domain, so the raw yield stress changes sign inside the mesh. + cm.Parameters.yield_stress = 1.0 + 0.5 * stokes.p.sym[0] + cm.yield_mode = mode + if smoother is not None: + cm.yield_smoother = smoother + + eta = cm.viscosity.sym + vals = uw.function.evaluate(eta, mesh.X.coords) + assert np.all(np.isfinite(vals)), ( + f"viscosity is non-finite for yield_mode={mode!r} smoother={smoother!r} " + f"where the pressure-dependent yield stress goes negative" + ) + + +def test_solve_homotopy_marches_and_reports(): + """A viscoplastic solve driven by the homotopy converges and reports the march.""" + _, stokes = _viscoplastic_stokes() + report = stokes.solve(homotopy=True, + homotopy_options=dict(delta0=1.0, dmin=0.05, verbose=False)) + + assert report["converged"] is True + # The march must actually DESCEND, not just solve once at delta0 and stop. + assert report["steps"] > 1, "the march did not take a second step" + assert report["settled_delta"] is not None + assert report["settled_delta"] < 1.0, ( + f"delta never descended below delta0 (settled at {report['settled_delta']})" + ) + assert report["reached_dmin"] is True + assert stokes.has_solution is True + # The model is left holding the δ that actually converged. + assert stokes.constitutive_model.yield_softness == pytest.approx( + report["settled_delta"] + ) + + +def test_homotopy_restores_the_solver_tangent(): + """The march sets the tangent the model asks for, but must hand the solver back + as it found it (adversarial review, M5).""" + _, stokes = _viscoplastic_stokes() + stokes.consistent_jacobian = False + stokes.solve(homotopy=True, homotopy_options=dict(delta0=1.0, dmin=0.1)) + assert stokes.consistent_jacobian is False, ( + "solve(homotopy=True) left the user's tangent permanently changed" + ) + + +def _yielding_box(tag, tau_y, cellSize=0.2): + """A box sheared hard enough to yield over a large fraction of the domain. + + NOTE the horizontally-VARYING body force. A uniform one is hydrostatic — pressure + balances gravity, nothing moves, the strain rate is zero and the yield law never + engages. Every earlier "viscoplastic" fixture in this file yields 0% for exactly + that reason, which is why the homotopy went so long without being exercised. + """ + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=cellSize + ) + x, y = mesh.X + v = uw.discretisation.MeshVariable("Vy" + tag, mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("Py" + tag, mesh, 1, degree=1, continuous=True) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + s.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel + cm = s.constitutive_model + cm.Parameters.shear_viscosity_0 = 1.0 + cm.Parameters.yield_stress = tau_y + s.bodyforce = sympy.Matrix([[0.0, -2.0 * sympy.cos(sympy.pi * x)]]) + s.add_essential_bc((sympy.oo, 0.0), "Top") + s.add_essential_bc((sympy.oo, 0.0), "Bottom") + s.add_essential_bc((0.0, sympy.oo), "Left") + s.add_essential_bc((0.0, sympy.oo), "Right") + s.petsc_use_pressure_nullspace = True + s.tolerance = 1.0e-8 + return mesh, s, cm + + +@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. + """ + import numpy as np + + # (a) the direct cold solve of the sharp law FAILS + mesh, cold, cm_cold = _yielding_box("c", 0.30) + cm_cold.yield_mode = "min" + cold.solve() + cold_reason = int(cold.snes.getConvergedReason()) + + # (b) the homotopy, same problem, SUCCEEDS + mesh2, warm, cm_warm = _yielding_box("h", 0.30) + report = warm.solve(homotopy=True, + homotopy_options=dict(delta0=1.0, dmin=1.0e-3, verbose=False)) + + eta = uw.function.evaluate(cm_warm.viscosity.sym, mesh2.X.coords) + yielding = float(np.mean(eta < 0.99)) + + assert yielding > 0.2, ( + f"fixture is not exercising the yield law (only {yielding:.0%} yielding) — " + "the comparison would be vacuous" + ) + assert cold_reason < 0, ( + f"the cold sharp solve unexpectedly converged (reason={cold_reason}); this " + "fixture no longer demonstrates a rescue, retune tau_y" + ) + assert report["converged"] is True, ( + f"HOMOTOPY REGRESSION: the march no longer rescues a case the cold solve " + f"cannot do ({report})" + ) + assert warm.has_solution is True + + +@pytest.mark.level_2 +def test_rate_regularisation_is_wired_into_the_plastic_viscosity(): + """``strainrate_inv_II_min`` caps eta_pl at tau_y/(2 edot_min), bounding the + viscosity contrast — the knob the xi-style regularisation needs. + + Regression: it was declared on this model but never applied (the elastic models + always used it), so setting it had no effect at all. + """ + _, _, cm = _yielding_box("r", 0.30) + before = str(cm.viscosity.sym) + cm.Parameters.strainrate_inv_II_min = 0.05 + after = str(cm.viscosity.sym) + assert after != before, ( + "strainrate_inv_II_min does not change the viscosity — it is being ignored" + ) + + +@pytest.mark.parametrize("kwargs", [dict(zero_init_guess=True), dict(picard=2)]) +def test_homotopy_rejects_arguments_it_would_have_to_ignore(kwargs): + """The march decides cold-vs-warm and its own warm-up per step, so an argument + that contradicts it is refused rather than silently dropped (review, m7).""" + _, stokes = _viscoplastic_stokes(cellSize=0.5) + with pytest.raises(ValueError): + stokes.solve(homotopy=True, **kwargs) + + +def test_homotopy_refuses_a_stress_history_solver(): + """A march is several solves; on a VEP solver each one would advance the elastic + stress history by a full timestep (adversarial review, C3). Refuse loudly rather + than silently integrating N steps for one requested dt.""" + mesh = uw.meshing.UnstructuredSimplexBox( + minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.5 + ) + v = uw.discretisation.MeshVariable("Vve", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("Pve", mesh, 1, degree=1, continuous=True) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel + cm = stokes.constitutive_model + cm.Parameters.shear_viscosity_0 = 1.0 + cm.Parameters.shear_modulus = 10.0 + cm.Parameters.dt_elastic = 0.1 + cm.Parameters.yield_stress = 5.0 + stokes.bodyforce = sympy.Matrix([0.0, -1.0]) + stokes.add_essential_bc((0.0, 0.0), "Bottom") + + with pytest.raises(NotImplementedError, match="stress history"): + stokes.solve(homotopy=True, timestep=0.1) diff --git a/tests/test_1055_solve_report.py b/tests/test_1058_solve_report.py similarity index 100% rename from tests/test_1055_solve_report.py rename to tests/test_1058_solve_report.py