Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
35 commits
Select commit Hold shift + click to select a range
7775c5f
feat(constitutive): power-mean yield soft-min smoother + rampable del…
lmoresi Jul 22, 2026
f8df9ff
fix(constitutive): DP consistent-Newton recursion on Max yield + smoo…
lmoresi Jul 22, 2026
691e0bb
feat(systems): yield_continuation driver — multi-solve delta-continua…
lmoresi Jul 24, 2026
753b042
docs(design): nonlinear-solver automatic warm-start + single-paramete…
lmoresi Jul 24, 2026
940d743
feat(solver): has_solution warm-start status + cold consistent-Newton…
lmoresi Jul 24, 2026
c13e69a
docs(solver): nonlinear-solver skill (recipe + config trap list) + ma…
lmoresi Jul 24, 2026
9789630
feat(solver): default the FMG velocity smoother to gmres+sor, fixed-c…
lmoresi Jul 25, 2026
528055b
docs(solver): record the measured Layer 3 result and the multigrid-de…
lmoresi Jul 25, 2026
6ddafd6
feat(solver): model-advertised yield homotopy and stokes.solve(homoto…
lmoresi Jul 25, 2026
7a97505
docs(solver): Layer 2 landed — one-call homotopy, and the cold-start …
lmoresi Jul 25, 2026
39572ef
fix(constitutive): power-mean yield NaN at zero strain rate (rigid/un…
lmoresi Jul 26, 2026
a5e8372
docs(solver): correct the cold-start diagnosis — inf-safe algebra, no…
lmoresi Jul 26, 2026
476a662
feat(solver): zero_init_guess becomes tri-state and auto-detects cold…
lmoresi Jul 26, 2026
7a27a11
docs(solver): all four layers landed; drop the last stale cold-start …
lmoresi Jul 26, 2026
69e82a2
fix(constitutive): yield_stress_min defaults to 0, floored via symbol…
lmoresi Jul 26, 2026
ca26077
fix(solver): close the critical adversarial-review findings in the ho…
lmoresi Jul 26, 2026
a6db75d
fix(constitutive/solver): close review findings M6, M7, M9
lmoresi Jul 27, 2026
ef2765b
fix(solver): close the remaining review minors and add parallel coverage
lmoresi Jul 27, 2026
eec04ce
feat(constitutive): apply the rate regularisation on the DP model; pi…
lmoresi Jul 27, 2026
3d2d27b
chore: drop the constitutive_models.py.bak copy committed by accident
lmoresi Jul 27, 2026
da5a37b
feat(solver): wall-clock guard and sub-solve work gauge
lmoresi Jul 27, 2026
1121368
Merge remote-tracking branch 'origin/feature/nonlinear-warmstart-homo…
lmoresi Jul 27, 2026
7e90eba
Merge branch 'feature/solver-wallclock-guard' into feature/regime-dia…
lmoresi Jul 27, 2026
384fafd
fix(test): pin the JIT-constant ramp test to a cold initial guess
lmoresi Jul 27, 2026
2206c00
fix(solver): close the adversarial-review findings on the wall-clock …
lmoresi Jul 27, 2026
928e7fe
docs(review): adversarial review of the wall-clock guard
lmoresi Jul 27, 2026
968a3f1
Merge branch 'feature/solver-wallclock-guard' into feature/regime-dia…
lmoresi Jul 27, 2026
a552b55
Merge remote-tracking branch 'origin/feature/nonlinear-warmstart-homo…
lmoresi Jul 27, 2026
c105994
fix(solve-report): SNES converged-reason table was off by one against…
lmoresi Jul 27, 2026
cfa10de
feat(constitutive): yield_anchor — which side of exact Min the soft-m…
lmoresi Jul 29, 2026
f1555b3
fix(solver): make the Stokes solver settings reachable (#477, D-22/D2)
lmoresi Jul 29, 2026
5ad68b2
docs(solvers): the Stokes fieldsplit is two nested Krylov loops
lmoresi Jul 29, 2026
954ff4f
fix(solver): make Eisenstat-Walker switchable
lmoresi Jul 29, 2026
5d5a59e
viscosity_min_rounding: give the floor its own rounding scale, decoup…
lmoresi Jul 30, 2026
ea99e45
Merge remote-tracking branch 'origin/development' into feature/regime…
lmoresi Jul 30, 2026
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 26 additions & 0 deletions docs/developer/design/nonlinear-solver-homotopy-warmstart.md
Original file line number Diff line number Diff line change
Expand Up @@ -146,6 +146,32 @@ multi-parameter is not, and is unnecessary here. The rate-strengthening ξ term
*non-homotopic* regularisation, not a homotopy — it belongs in a user-domain loop
*around* `solve()`, not inside it.

#### Which side of `Min` the entry problem sits on (`yield_anchor`)

Smoothing a corner by δ moves the whole curve, so each soft-min family has one free
constant. It is fixed by choosing which point the smoothed law must reproduce exactly,
and that choice — not the choice of family — decides whether the regularised problem is
*stiffer* or *weaker* than the sharp one. Writing $f = \eta_{ve}/\eta_{pl}$ for the
overstress ratio ($f = 1$ is the yield point):

| `yield_anchor` | exact at | sqrt family | power mean |
|---|---|---|---|
| `"onset"` (default) | $f \to 0$, the unyielded branch | undershoots for $f < 2$ | the p-norm; $\eta \to 0$ as δ grows |
| `"yield"` | $f = 1$, the yield point | offset $\delta/2$; $\ge$ `Min` | the true power mean (the $1/2$); $\ge$ `Min` |

**A homotopy wants `"yield"`.** An entry problem that sits *below* the sharp law is
weaker than the problem it is supposed to lead to, which is the opposite of an easier
version of it. Under `"onset"` both families are weaker near the yield point, and that
bites hardest exactly where the overstress ratio is O(1) — since the rate term floors
the plastic viscosity ($\eta_{pl} \ge \xi$, so $f \le 1/\xi$), a held ξ puts the whole
domain in the softened region. The cost of `"yield"` is stiffened unyielded material,
bounded by 2 for the sqrt family and by $2^{\delta}$ for the power mean; that bound is
why the power mean's entry δ is O(1) while the sqrt family's is O(10).

The default is `"onset"` because it is the historical law, not because it is the better
one. Changing it alters results for every existing `yield_mode="softmin"` user, so it is
a separate maintainer decision.

### Layer 3 — smoother as a consistent-Newton consequence

Turning on the consistent tangent adds the `∂η/∂(grad v)` term, which makes the
Expand Down
106 changes: 106 additions & 0 deletions docs/developer/subsystems/solvers.md
Original file line number Diff line number Diff line change
Expand Up @@ -178,6 +178,112 @@ def validate_unknowns_sharing(multi_material_model):
- ⚠️ Preconditioner *selection* is partly covered — see "Choosing the Krylov method for a fieldsplit sub-solve" and "The multigrid option bundle, and its one owner" below; Schur preconditioner choice is still undocumented
- Could benefit from optimization examples

## The Stokes fieldsplit: two nested Krylov loops

The Stokes saddle-point solver is configured as

```
pc_type = fieldsplit
pc_fieldsplit_type = schur
pc_fieldsplit_schur_fact_type = full
pc_fieldsplit_schur_precondition = a11
```

which is easy to read as one solver and is in fact **two nested Krylov loops with a third
inside them**. Knowing which loop does the work is the difference between tuning the
solver and permuting options.

For the saddle system

$$
\begin{bmatrix} A & B^{T} \\ B & 0 \end{bmatrix}
\begin{bmatrix} u \\ p \end{bmatrix} =
\begin{bmatrix} f \\ g \end{bmatrix}
$$

`fieldsplit` is a **preconditioner for a Krylov method on the full coupled system**. With
`schur_fact_type = full`, one application of that preconditioner is the block
factorisation, which needs $A^{-1}$ twice and $\hat{S}^{-1}$ once, where
$S = -B A^{-1} B^{T}$ is the Schur complement and $\hat S$ is its preconditioner — with
`schur_precondition = a11`, the $1/\mu$-weighted pressure mass matrix that
`saddle_preconditioner` supplies.

The three solves:

| loop | iterates on | cost of one iteration |
|---|---|---|
| **outer KSP** | the full coupled residual | one preconditioner application |
| **pressure sub-KSP** | the Schur system $S p = r$ | a MatMult by $S$ — i.e. **a velocity solve** |
| **velocity sub-KSP** | $A u = r$, preconditioned by FMG or GAMG | one multigrid cycle per Krylov step |

### Two designs, and the same work in different loops

The pressure block's `ksp_type` chooses between two genuinely different algorithms.

**A pressure Krylov solve gives you Citcom.** The Schur system is actually solved, so the
block preconditioner is close to the exact inverse and the outer KSP converges in one or
two iterations — a formality wrapped around a Uzawa solve. This is the classical
arrangement (Moresi & Solomatov 1995): an outer *pressure* iteration with an inexact
multigrid velocity solve inside it. UW3's inner-tolerance margins (`0.033` velocity,
`0.1` pressure, applied by the `tolerance` setter) were designed for this shape.

**`fieldsplit_pressure_ksp_type = preonly` gives you a block preconditioner.** The Schur
system is never solved; $\hat S^{-1}$ is a single mass-matrix application, the fieldsplit
becomes a cheap approximate preconditioner, and the **outer** Krylov does all the coupling
work. This is the Elman–Silvester–Wathen approach, standard in the finite-element
literature. It is not a degenerate configuration — but it is not Citcom, and the margins
above are not aimed at it.

Neither is universally better, and UW3 does not take a position: the choice is
problem-dependent and worth measuring on the problem at hand. What matters is knowing
which one is in force, because the diagnostics differ.

```{important} The two choices interact with the outer restart
Under `preonly`, *every* coupling iteration is an **outer** iteration — and the outer loop
is where GMRES's restart lives. A problem needing more outer iterations than
`ksp_gmres_restart` (PETSc default 30) discards its Krylov space and stagnates, exactly
where a deep residual is being ground out. Under a pressure Krylov the outer count stays
at one or two, so the restart never bites there; it moves into the pressure KSP, which can
be given its own.

A `preonly` configuration therefore wants a **flexible outer method and a generous
restart**. Note also that the outer preconditioner *varies between applications* whenever
a sub-block is itself a Krylov solve — which a non-flexible GMRES assumes away.
```

```{tip} Report the outer iteration count
`solve_report.sub["velocity"].its` is the total multigrid cycle count — a **cost** measure.
A restart stagnation looks like a large **outer** count (`solve_report.ksp_its`) with a
poor residual, and is invisible in the velocity total. When diagnosing a Stokes solve that
grinds, print both, plus `stokes.snes.getKSP().getConvergedReason()`.
```

### Inner tolerances are deliberately inexact, and bounded

The inner solves are inexact on purpose — you do not solve the velocity block exactly in
order to apply the Schur complement. Two consequences hold together:

1. inexact inner solves perturb the outer search directions, so **the outer Krylov must be
flexible**; and
2. the inexactness is **bounded** — the inner solves must still converge well below the
tolerance demanded of the outer solve. The `0.033` and `0.1` factors are that margin.

"Flexible or not" and "how tight is the inner tolerance" are two halves of one decision,
not independent axes: flexibility buys tolerance of inexactness, and the margin bounds how
inexact. `preonly` is the degenerate end — no tolerance, hence no margin — which is why it
belongs to a design where the *outer* loop is doing the converging.

Both `fieldsplit_velocity_ksp_rtol` and `fieldsplit_pressure_ksp_rtol` are settable, and
setting `tolerance` re-derives them from it — so set `tolerance` first, then any override.

```{warning} `preonly` under the Schur is a different matter
`preonly` on the **pressure** block is a design choice, as described above. `preonly` on
the **velocity** block, under `schur_fact_type = full`, is a defect: PCFieldSplit applies
$A^{-1}$ *through* the velocity sub-KSP when forming the Schur action, so an inexact
velocity solve hands the pressure Krylov a different operator than the one its
preconditioner was built for.
```

## Choosing the Krylov method for a fieldsplit sub-solve

This choice gets re-argued periodically. The confusion is that it looks like three
Expand Down
Loading
Loading