Skip to content

yield_anchor: which side of exact Min the soft-min yield law sits on - #475

Open
lmoresi wants to merge 30 commits into
developmentfrom
feature/regime-diagram
Open

yield_anchor: which side of exact Min the soft-min yield law sits on#475
lmoresi wants to merge 30 commits into
developmentfrom
feature/regime-diagram

Conversation

@lmoresi

@lmoresi lmoresi commented Jul 29, 2026

Copy link
Copy Markdown
Member

What this is

A soft-min rounds the yield corner by δ, and rounding a corner moves the whole curve — so
each soft-min family has one free constant. Until now it was fixed implicitly, and
differently in each family, with the consequence buried in the algebra: both families
ended up sitting BELOW exact Min at and above the yield point.

For a homotopy that is the wrong direction. The entry problem is supposed to be an easier
version of the sharp problem; a law that sits below it is a weaker problem, which is a
different thing. The sqrt family undershoots for overstress ratios $f < 2$ — a third
under-stressed at nominal yield — and the power mean degenerates to the p-norm
$(\eta_{ve}^{-s} + \eta_{pl}^{-s})^{-1/s}$, whose viscosity falls away toward zero as δ
grows.

yield_anchor names the choice and makes it explicit, for both families:

yield_anchor pinned exactly at sqrt family power mean
"onset" (default, historical) $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; $\ge$ Min

The power-mean case is the neat one: it needs no offset at all. Averaging the two
terms — the /2 — is exactly what makes it a generalised mean rather than a norm, and
a generalised mean returns the common value when its arguments are equal. So
$\eta = \eta_{ve} = \eta_{pl}$ at $f = 1$ falls out of the algebra.

The cost, stated plainly

"yield" stiffens unyielded material: bounded by a factor 2 for the sqrt family, and by
$2^{\delta}$ for the power mean. At the power mean's own entry δ₀ = 1 that is a factor 2
on the entire background. Both decay to 1 as δ → 0, so the endpoint is unchanged, but the
entry problem is not merely "the sharp problem with a rounded corner" — the background is
genuinely stiffer along the path. That bound is also why the two families' entry δ are not
the same number, which is the second half of this change.

δ is not one parameter

For the power mean the sharpness is $s = 1/(\delta + 0.001)$, so δ ≤ 1 and δ = 1 is the
harmonic mean. For the sqrt family δ is a percentage stress deviation and a generous entry
is O(10). A sqrt-family δ = 16 handed to the power mean asks for a law a factor 2^16 away
from Min. The entry now comes from the family (_YIELD_HOMOTOPY_DELTA0) rather than
from a caller carrying one number across both, and yield_continuation grows smoother=
and anchor= so a march can pick both — refusing them, rather than silently ignoring
them, when handed a ready-made control.

What is NOT changed

The default stays "onset". Flipping it changes results for every existing
yield_mode="softmin" user, and that is a maintainer decision rather than a side effect of
naming the knob. Left open deliberately — see the question below.

Tests

34 passed. test_1059_yield_anchor.py is new and parametrised over both families; each
test that pins the new behaviour also asserts the OLD anchor fails it, so none of them
can pass against a no-op implementation. test_1055 now sets the anchor explicitly where
it exercises the below-Min side, rather than relying on the default.

Open question for the maintainer

Should the default flip to "yield"? The argument for: it is the correct direction for a
homotopy, and it is what any new user opting into softmin for continuation will want.
The argument against: it changes converged results (within tolerance, but not identically)
for anyone already using softmin as a plain regularisation rather than as a march, and
for the power mean it stiffens the background by up to 2^δ.

Underworld development team with AI support from Claude Code

lmoresi added 30 commits July 27, 2026 17:20
…ta constant

Add the power-mean soft-min yield smoother and make the soft-min softness delta a
constants[] atom (runtime-rampable with no JIT recompile) on the ViscousFlowModel base
class, so the visco-plastic subclasses inherit one implementation. This is the
generalisable, scalable substrate for a yield homotopy toward the sharp Min surface.

Base class (ViscousFlowModel):
- _combine_yield(eta_ve, eta_pl): the shared viscous/plastic combination, keyed on
  yield_mode -- "min" (exact hard Min), "harmonic", or "softmin" (the delta soft-min in
  the family chosen by yield_smoother). Development's tri-modal yield_mode semantics are
  preserved exactly (unlike the source branch, which collapsed "softmin" into "min").
- _get_yield_softness / _get_yield_offset: delta and the onset offset held as constants[]
  UWexpression atoms (the offset defined symbolically in terms of delta -- one symbol in
  the stress tensor, so no tensor blow-up -- so one delta update repacks both). Ramping
  delta via the atom + solver._update_constants() forces no recompile.
- yield_smoother property/setter: "sqrt" (default; overshoots tau_y in the transition) or
  "powermean" (undershoots tau_y -- eta_eff <= Min always -- overflow-safe on geodynamic
  viscosity ranges). Selecting "powermean" from delta=0 bumps delta to 1 (s=1 harmonic
  mean), since delta=0 is the singular Min limit.

Models (Stage 1 -- isotropic, test-covered):
- ViscoPlasticFlowModel (Drucker-Prager): defaults to yield_mode="min" == today's exact
  hard Min (zero behaviour change), and GAINS yield_mode / yield_softness knobs so it can
  opt into the softmin / power-mean homotopy -- the model the hard-case DP homotopy needs.
- ViscoElasticPlasticFlowModel (VEP): its inline sqrt soft-min routed through _combine_yield;
  delta moves float -> atom (value-identical); yield_softness setter syncs the atom.

Zero behaviour change at stock settings: _combine_yield at default settings is
bit-identical to the previous inline law (delta merely moves from a baked float to a
constants[] symbol evaluated to the same number). Power-mean and runtime-delta are opt-in.
The distrusted in-solve enable_yield_homotopy ramp driver is intentionally NOT part of
this PR. TI-VEP / TI-VEP-Split (4 anisotropic sites) are a Stage-2 follow-up.

Tests: tests/test_1055_yield_smoother.py -- delta=0 exact Min; power-mean undershoot +
overflow-safety + ->Min; smoother validation/default; runtime delta ramp no-recompile;
DP opt-in. level_1/tier_a gate green (402 passed, 0 failed). VEP _combine_yield verified
bit-identical to the prior float law at defaults.

Underworld development team with AI support from Claude Code
…th floors

The Drucker-Prager viscosity guarded its yield floor with
`if yield_stress_min.sym != 0:`, but the unset default yield_stress_min is a
-inf-valued Parameter (a UWexpression atom), which passes `!= 0`. So every yield
became `Max(<-inf parameter>, yield_stress)`. That -inf parameter cannot be
resolved by sympy's fuzzy `is_ge`, so canonicalising the outer Max for the
consistent-Newton tangent recursed without bound (is_ge -> _n2 ->
(a-b).evalf(2) -> is_ge ...). The guard now uses the model's own -inf "unset"
sentinel, so the spurious wrapper is not applied; a yield_stress_min of exactly
0 is now honoured as a real floor.

Also, because a hard Max floor reintroduces a non-differentiable kink under the
consistent tangent, route both DP floors (shear_viscosity_min, yield_stress_min)
through `_apply_floor`: exact hard Max in yield_mode="min" (unchanged), and a
delta-rounded smooth maximum in the smooth modes so the whole effective viscosity
stays differentiable (the same delta that regularises the yield transition).

Add public uw.maths.smooth_max / smooth_min primitives (regularised like
delta_function) - the Newton-safe way to write a physical rounded tension cap,
smooth_max(C + sin(phi)*p, 0, eps), in place of sympy.Max(sigma, 0).

Tests (test_1055): a consistent-Newton Jacobian build on a pressure-dependent
Max(C+sin(phi)*p, 0) yield (regression for the recursion), and a smooth_max/
smooth_min contract test.

Underworld development team with AI support from Claude Code
…tion (reference impl)

The trusted homotopy for hard viscoplastic (Drucker-Prager) yield: hold delta
constant for a full nonlinear solve to tolerance, warm-start the next (smaller)
delta from the converged state, march down and settle at the smallest feasible
delta. Extracted from the proven Spiegelman harness (convergence.py::
run_continuation). delta lives in constants[], so each step is a recompile-free
PetscDSSetConstants update.

This is the reference implementation the design doc
(nonlinear-solver-homotopy-warmstart) folds into stokes.solve(homotopy=True). It
orchestrates solves only; the caller configures the tangent/preconditioner/
smoother (see the driver's Notes and the design doc's trap list).

NOT the in-SNES delta-ramp, which is proven non-viable (see the design doc).

Underworld development team with AI support from Claude Code
…r yield homotopy

PROPOSED design for a dedicated, benchmarked implementation session. Captures:
- why setup is the real defect (fiddly config = API regression), with the trap list;
- the evidence: multi-solve delta-continuation converges the hard case, the in-SNES
  delta-ramp does NOT (proven on the proven config: DIVERGED_LINEAR_SOLVE, ~2h) and
  the mechanism why (sharpen delta only from a converged, well-conditioned iterate);
- the three-layer design: automatic Picard-on-cold warm-start + has_solution property
  + zero_init_guess auto-detect; constitutive-model-advertised single-parameter
  homotopy behind solve(homotopy=True) with residual-guided auto-descent; a
  non-symmetry-safe multigrid smoother default under the consistent tangent;
- staged implementation plan and open verification points (free-surface chain, mover
  reset, single-Picard-step-lands-in-basin, benchmark the default changes).

Underworld development team with AI support from Claude Code
… Picard warm-up (Layer 1a)

Layer 1a of the nonlinear-solver automatic warm-start / single-parameter yield
homotopy design (docs/developer/design/nonlinear-solver-homotopy-warmstart.md).

SolverBaseClass.has_solution — a public, read-only warm-start status flag. Set
True only after a solve whose SNES converged (reason > 0); reset to False
initially, after a diverged solve, and on a structural rebuild (the
is_setup=False invalidation used by remesh / adapt / mesh-mover). It survives
coefficient-only changes (viscosity, yield softness delta, BC values, time step),
which route through _update_constants / a direct function rewire, so parameter
continuation and time-stepping warm-start correctly. Recorded by
_record_convergence_status() at the end of every solve() — scalar, vector,
multi-component, Stokes, and the rotated free-slip path (from its info dict,
since that path runs its own KSP loop rather than self.snes).

Stokes cold-start warm-up — a single Picard step (reusing the existing picard=1
nrichardson machinery) is taken automatically on a cold (zero_init_guess) solve
under the consistent-Newton tangent, the opt-in nonlinear-yield regime that
needs a defect-correction step into the Newton basin. The default (frozen)
tangent path is bit-identical and no default is flipped — the zero_init_guess
auto-detect flip is the separate, benchmarked Layer 1b.

Additive and low-risk. Regression test test_0201 covers the has_solution
lifecycle (fresh / converged / structural-reset / coefficient-survival /
mesh-deform) and that the cold consistent-Newton warm-up converges.

Underworld development team with AI support from Claude Code
…rk Layer 1a landed

Seed a `nonlinear-solver` skill capturing the working recipe for hard
viscoplastic Stokes convergence — automatic warm-start plus the MULTI-SOLVE
δ-continuation (not an in-SNES ramp), the consistent-Newton tangent, and the
non-symmetry-safe smoother — together with the CONFIG TRAP LIST (the setup
mistakes that each produce a different failure a few steps in) and the
diagnostic tell (get_snes_diagnostics linear_iterations ≈ 1 = no linear work).
It records the proven-dead in-SNES δ-ramp as a do-not, superseding the
enable_yield_homotopy advice in the plasticity-solvers skill for the hard case,
and cross-links plasticity-solvers for the yield-law maths and tangent-per-model.

Mark Layer 1a landed in the design doc with the as-built scoping decisions
(warm-up gated on the consistent-Newton tangent; rotated path records via its
info dict; has_solution resets in the is_setup setter).

Underworld development team with AI support from Claude Code
…ost V-cycle (Layer 3)

Layer 3 of the nonlinear-solver design. The geometric-FMG option bundle now uses a
non-symmetry-safe Krylov smoother sized for a DEEP hierarchy, rather than the
stationary richardson+sor that preceded it:

  mg_levels_ksp_type      richardson -> gmres
  mg_levels_pc_type       sor        (unchanged)
  mg_levels_ksp_max_it    4          (unchanged)
  mg_levels_ksp_norm_type (new) none  -> fixed-cost V-cycle

Measured on the Spiegelman notch (Drucker-Prager, eta contrast 1e26, V=10) with a
NESTED hierarchy built by uniform refinement of an ultra-coarse 492-cell base, so
the multigrid actually has depth. Per-V-cycle contraction rho of the velocity block
at xi=0.01, at the SAME four smoother iterations:

  3 levels:  richardson 0.722   gmres 0.686   (5% better)
  4 levels:  richardson 0.746   gmres 0.560   (25% better)

The gmres margin GROWS with depth, because a deeper cycle applies the smoother on
more coarse operators, each non-symmetric under the consistent-Newton tangent. A
two-level cycle is a coarse-grid correction rather than a V-cycle and is not worth
special-casing, so 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 (rho^(1/4) = 0.87) beats
gmres/8 (0.478^(1/8) = 0.91).

norm_type=none makes the smoother run exactly max_it iterations with no residual
norm and no early exit, so every V-cycle costs the same. The resulting cycle is
non-stationary; the velocity block is already fgmres (flexible), which is what
tolerates that, so no other change is needed.

Not fixed by this: at the hard small-xi corner the inner solve fails under every
smoother tried (richardson outright, gmres with rho > 1 at 4 levels). That is
operator conditioning, and the lever for it is the delta/xi continuation (Layer 2),
not the smoother.

Verified: test_1014 (contract updated to the new bundle + the fixed-cost assertion),
FMG checkpoint hierarchy, Stokes cartesian/spherical/nullspace/rotated-freeslip and
the solver smoke suite (34 passed), plus ptest_0004 FMG hierarchy at np=2.

Underworld development team with AI support from Claude Code
…pth method

The Layer 3 design section predicted a one-line smoother swap gated on
consistent_jacobian. The benchmark changed the shape of the fix, so record what was
actually measured and what landed:

- the gmres-over-richardson margin GROWS with multigrid depth (5% at 3 levels, 25% at
  4), because a deeper cycle applies the smoother on more non-symmetric coarse
  operators — so a comparison run on a 2-level hierarchy is misleading;
- therefore the default is unconditional (no level-count or tangent gate), since a
  two-level cycle is a coarse-grid correction rather than a V-cycle and is not worth
  special-casing;
- four smoother iterations, not more, on a per-unit-work basis;
- the small-xi corner is NOT fixed by any smoother (richardson fails outright, gmres
  gives rho > 1 at 4 levels) — that is operator conditioning and belongs to the
  delta/xi continuation.

The skill gains a "Multigrid depth" section with the method: build depth from an
ultra-coarse NESTED base (492 cells + uniform refinement) rather than a fine mesh,
never a non-nested hierarchy, and probe the fieldsplit_velocity_ sub-KSP because
solve_report records the Newton rate and is blind to the smoother. Its trap-list row
for the non-symmetric smoother is now "this is the default" rather than a manual fix.

Underworld development team with AI support from Claude Code
…py=True) (Layer 2)

Layer 2 of the nonlinear-solver design. A hard Drucker-Prager problem does not
converge from a cold start on the sharp yield surface; the route that works is a
multi-solve continuation in the yield softness delta. This makes that the API rather
than a hand-built driver.

The constitutive model ADVERTISES the homotopy, exactly as it already advertises
stress-history support:

  supports_yield_homotopy   False on Constitutive_Model, True on the viscoplastic
                            and (TI-)VEP models.
  _yield_homotopy_control() switches the model into its smooth delta-parameterised
                            mode (softmin + power-mean) and returns a
                            YieldHomotopyControl: a model-owned delta setter, the
                            tangent to pair with it, and the delta atom.

The model owns what the parameter MEANS: the isotropic models update delta as a
constants[] atom (no recompile) while TI-VEP rebuilds, and the elastic models ask for
the frozen/Picard tangent because the consistent yield tangent across the elastic
stress-history block is indefinite and fails the linear solve outright.

stokes.solve(homotopy=True, homotopy_options=...) runs the continuation and returns
the march summary. Each delta-step re-enters solve() normally, so nothing recurses;
a VEP model's timestep is forwarded to the inner solves. Refusing an unsupported
model is explicit (there is no yield surface to sharpen).

The march itself is now residual-guided: a step that converges well inside its
iteration budget earns a bigger next step, a step that nearly exhausts it earns a
gentler one, and a failed delta is reverted and retried more gently before the march
settles. Measured on a small viscoplastic box it walks delta from 1.0 to 7.8e-4 in
four solves from a COLD start.

Cold-start safety: _yield_homotopy_control floors the strain-rate invariant with the
vanishing constant, because tau_y/(2 edot_II) is 0/0 at v=0 and the first residual is
otherwise NaN. Carries a TODO(BUG) — the hazard is not homotopy-specific, any cold
solve of a viscoplastic model with a finite yield stress hits it, and the floor
arguably belongs in the model.

Verified: new test_1057 (advertisement, control contract, tangent per model, refusal,
and an end-to-end march), plus yield-smoother, VEP stability, VE Stokes, has_solution
and multigrid regressions (34 passed).

Underworld development team with AI support from Claude Code
…NaN correction

Record what Layer 2 actually became: solve(homotopy=True) is the entry point, the
model-owned delta setter is what keeps delta from being reset by a later
_get_yield_softness(), the march reverts and retries a failed delta before settling,
and a VEP model's timestep is forwarded to the inner solves.

Correct a claim in the design and in the skill's trap list: starting delta LARGE does
not, on its own, make a cold plastic start safe. The singularity is tau_y/(2 edot_II)
inside the plastic viscosity, which is 0/0 at v=0 and NaNs before the soft-min ever
combines it with the viscous branch. The strain-rate invariant has to be floored.

Underworld development team with AI support from Claude Code
…yielded points)

A cold (v=0) viscoplastic solve died with DIVERGED_FNORM_NAN in yield_smoother=
"powermean", while the hard-Min and sqrt soft-min forms converged. The cause is not a
division by zero in the model: at edot_II = 0 the plastic viscosity tau_y/(2 edot_II)
is +inf, and infinity propagates CORRECTLY through Min and through the sqrt form to
give the viscous branch. It is specific to the power-mean algebra, which computed its
harmonic-mean prefactor as

    N = eta_ve * eta_pl / (eta_ve + eta_pl)     -> inf/inf = NaN

Written instead as the algebraically identical

    N = eta_ve / (1 + f),   f = eta_ve/eta_pl

the infinite-eta_pl limit is finite and correct (N -> eta_ve). No strain-rate floor,
no change to any converged result — the yield-smoother accuracy tests are unchanged.

This is NOT only a cold-start concern: eta_pl is infinite at every rigid (unyielded)
point, which is a physically real state in a viscoplastic flow, so the old form could
poison a converged solve wherever a plug formed. A cold start is just the case where
every point is unyielded at once.

Supersedes the strain-rate floor added with the Layer 2 homotopy control, which was
compensating for this bug rather than fixing it; that floor and its TODO(BUG) are
removed.

Regression: test_1057 parametrises a cold viscoplastic solve over all three yield
modes; the powermean case failed (reason -4) before this fix.

Underworld development team with AI support from Claude Code
…t a strain-rate floor

The design note and the skill trap list both said a cold plastic start needs the
strain-rate invariant floored. Measurement says otherwise: hard Min and the sqrt
soft-min cold-start cleanly because eta_pl = +inf propagates correctly to the viscous
branch. Only a harmonic mean written as a product over a sum breaks (inf/inf), and the
fix is to write it as eta_ve/(1+f). Steer readers away from the floor, which hides the
problem, and note that the singularity is a property of rigid/unyielded points rather
than of cold starts.

Underworld development team with AI support from Claude Code
… vs warm (Layer 1b)

Layer 1b, the benchmarked default change. solve(zero_init_guess=...) was True by
default, so a solve() in a loop silently re-zeroed the initial guess every iteration
and threw away a perfectly good warm start. It is now tri-state:

  None (default)  auto-detect: cold when has_solution is False, warm when True
  True            force a fresh start, discarding any existing solution
  False           insist on warming from the current field values

The onus flips from "flag when you have a solution" to "flag when you want to throw
one away". Detection is safe by construction: guessing cold when a solution was
available costs one iteration from a good starting point, while the harmful direction
— warming off stale field data — cannot happen, because has_solution is cleared both
by a structural rebuild (remesh / adapt / mesh-mover) and by a diverged solve.

Design gates checked before flipping:
- Free-surface chain-of-solves (the pattern flagged as most at risk): its repeated
  bare solve() is a LINEAR Poisson mesh-displacement solve with constant diffusivity,
  so a warm start cannot change the answer it converges to.
- Mover/adapt reset: already covered by has_solution's is_setup hook and its
  regression test.

Measured consequence, and the reason this is a real (if small) behaviour change:
warm and cold agree only to the CONVERGENCE TOLERANCE, not bitwise, because an
iterative solve stops anywhere inside the tolerance ball and a warm start enters it
from a different direction. The difference scales with the tolerance —
2.4e-5 at tol=1e-4, 7.4e-10 at 1e-8, 2.0e-14 at 1e-12 — so a test with a threshold
tighter than its own solver tolerance could shift.

Verified: core solvers, Stokes cartesian, multigrid, deform/rebuild, solver smoke
(50 passed) and the repeated-solve suites where this bites hardest — advection-
diffusion cartesian and annulus, VEP stability, VE Stokes, yield smoother, transient
Darcy (21 passed, 1 pre-existing xpass on a test its own header marks fragile).

Underworld development team with AI support from Claude Code
…claims

Mark the design IMPLEMENTED and Layer 1b DONE, with the measured caveat that warm and
cold agree to the convergence tolerance rather than bitwise.

Remove the remaining places where the skill still said a cold plastic start NaNs
unless large delta or a strain-rate floor saves it. Neither is true: eta_pl is simply
infinite at zero strain rate and the soft-min carries that to the viscous branch. The
trap list already carries the accurate version.

Underworld development team with AI support from Claude Code
…ic smooth_max

Maintainer ruling (2026-07-26): a negative yield stress is meaningless — tau_y is
compared against the second invariant of the stress — so yield_stress_min defaults to
0 rather than the -oo "unset" sentinel. The +/-oo defaults elsewhere in Parameters
exist so sympy can cancel an unused term away; that trick is wrong for this one.

This closes two defects found by adversarial review:

* A pressure-dependent Drucker-Prager yield C + sin(phi)*p goes negative in tension.
  Hard Min merely hid the resulting negative eta_pl; the power-mean the homotopy
  switches to raised a negative base to a non-integer power and returned NaN, so
  solve(homotopy=True) failed on the first residual for exactly the class of model
  the feature targets.
* _apply_floor's rounding scale is relative (delta * floor) and collapsed at a zero
  floor, so honouring yield_stress_min = 0 gave a zero-width smooth_max — the hard
  kink with a 0/0 derivative.

The floor is applied through uw.maths.smooth_max everywhere, including yield_mode=
"min" where eps = 0 makes it exactly Max. That matters for more than smoothness:
smooth_max is pure arithmetic on its operands, so the floor stays the symbolic
Parameter the JIT routes through constants[] and sympy is never asked for an ordering
test. sympy.Max(<UWexpression>, x) cannot resolve that comparison and recurses until
the stack dies — the fragility noted when the DP sentinel guard was fixed, which the
old -oo default had merely hidden by skipping the floor. Turning the floor on by
default made it fire on every yielding model, so it had to be fixed rather than
dodged.

Guards now test the DISABLING sentinel (!= -sympy.oo) instead of != 0, which would
have skipped the floor at the new default. That also fixes the VEP and TI-VEP models,
where the != 0 guard silently ignored a user-set yield_stress_min = 0.

Left open, marked TODO(DESIGN): a genuinely rounded tension cap (Griffith /
parabolic) needs an ABSOLUTE stress scale, which _apply_floor's signature cannot
supply. At floor = 0 the cutoff is therefore still a corner, and eta_pl is exactly
zero in tension.

Regression: test_1057 asserts the zero default and that a yield stress driven
negative by tension leaves the viscosity finite in both hard-Min and power-mean modes.

Underworld development team with AI support from Claude Code
…motopy march

Review report added at docs/reviews/2026-07/nonlinear-solver-homotopy-adversarial-review.md
(four independent adversarial reviewers plus an author pass; every finding re-verified
against source). Codes below refer to it.

C1 - iteration budgets are real again. entry_maxit/step_maxit were pushed through
petsc_options, which SNES_Stokes_SaddlePt.solve overwrites with its hardcoded
snes_max_it=50 via setFromOptions immediately before solving, so both documented
knobs were silent no-ops and the residual-guided step logic was comparing the real
iteration count against a budget that never applied. They now go through the same
_difficulty_probe / _difficulty_max_it hook estimate_difficulty uses, which is applied
after setFromOptions and restored in a finally.

C2 / M2 - a reverted iterate is a CONVERGED iterate, so the march now says so. The
failed inner solve cleared has_solution, so the retry's zero_init_guess resolved cold
and re-solved a sharper delta from zero, defeating the revert entirely; and a march
that settled through the failure exit handed back a good solution with
has_solution=False, so the caller's next solve threw it away.

C3 - refuse the march on a stress-history (VEP) solver. Each delta-step re-enters the
full time-integration tail, so an N-step march advanced the elastic stress history by
N*dt for one requested dt. Failing loudly beats silently corrupting the history; the
interaction needs designing before VEP can use this.

M1 - zero_init_guess is resolved AFTER the _force_setup block in all four solve
bodies. Resolving first meant solve(_force_setup=True) warm-started off the
has_solution flag that the very same call was about to clear.

M3 - the entry solve can warm-start again: has_solution is read before the rebuild
that clears it, instead of one line after.

M4 - the march is wrapped in try/finally, so an exception no longer leaves the
tangent, the iteration budget and the model's delta in march state.

M5 - consistent_jacobian is restored to the caller's value.

m1, m3, m5, m6 - pprint no longer prints a stray rank argument (the retired
pprint_old signature); the retry log names the delta that actually failed; delta0 is
validated; and max_steps caps a march that keeps easing off.

Tests strengthened where the review showed they would pass with the feature deleted:
the march must now descend (steps > 1, settled < delta0, reached_dmin) rather than
merely run once; _force_setup cold-start is driven through the public solve(); and
new cases cover the tangent restore and the VEP refusal.

Still open from the review: M6 (power-mean cannot reach exact Min - docs overstate
it), M7 (TI-VEP ignores yield_smoother), M9 (Newton tangent at zero strain rate),
the remaining minors, and parallel coverage.

Underworld development team with AI support from Claude Code
M6 - MEASURED AND MOOT, as the maintainer predicted. The power-mean's sharpness
s = 1/(delta + 0.001) saturates at 1000, so it is not exactly Min at delta = 0. What
that costs in practice: on a 45%-yielded box (a case a cold hard-Min solve cannot
solve at all), the exact hard-Min residual at the settled power-mean solution is 6e-10
of the initial residual, an order of magnitude INSIDE the 1e-8 solver tolerance -- and
a hard-Min solve warm-started from it converges in 0 nonlinear iterations, i.e. PETSc
already considers the exact problem solved. The docstring now states the limit
honestly along with the measurement instead of claiming exactness.

Two things worth recording from getting there. Sharpening the law does NOT help: with
the floor at 1e-6 (s up to 1e6) the achievable residual got WORSE, because the law
becomes too stiff to solve before it becomes more accurate. And an earlier reading of
this measurement said the opposite -- it normalised by the last warm solve's initial
residual, which is itself tiny, instead of the physical ||F(v=0)||.

M7 - transverse-isotropic VEP now uses the same yield envelope as the isotropic
models. Four sites reimplemented the harmonic / soft-min / hard-Min choice inline,
hardcoding the sqrt family and reading delta as a baked float, so yield_smoother was
silently ignored there and delta could not be ramped without a recompile. They all
call _combine_yield now, and the model's yield_softness setter keeps the constants[]
atom in step.

M9 - the design requirement is met rather than approximately met. A viscoplastic
Jacobian is NaN at zero strain rate (the residual survives, its derivative does not),
so the warm/cold machinery must make that state unreachable. It covered an explicit
cold start but not a nominally WARM solve whose solution had never been written --
the design's stated "secondary signal", never implemented. _solution_is_trivially_zero
supplies it, using the collective vector norm so every rank decides alike. The
"continuation" tangent needs no guard: it opens on a Picard stage by construction.

The M9 test is verified non-vacuous: with the previous guard restored, the
zero_init_guess=False case dies inside PETSc's KSPSetUpOnBlocks on the NaN operator.

Also proven for the first time, and the reason these were reachable at all: the
homotopy had never been exercised on a genuinely yielding problem. Earlier cases
(including my own tests) sat at 0% yielding -- a uniform body force in a box is just
hydrostatic. On a properly sheared box, solve(homotopy=True) converges at 45% yielding
where a cold hard-Min solve gives DIVERGED_MAX_IT, which is the feature doing its job.

Underworld development team with AI support from Claude Code
m2 - a diverged LINEAR rotated free-slip solve no longer records success. The linear
helper returns only ksp_reason (a linear solve can fail; _warn_if_ksp_diverged exists
for exactly that), so the .get("converged", True) fallback claimed convergence and the
next solve would warm-start off the bad iterate. It now falls back to ksp_reason > 0,
which is what the sibling _capture_rotated_report already did -- the two diagnostics
had been able to disagree.

m7 - solve(homotopy=True) forwards the per-solve arguments it can honour (timestep,
evalf, order, debug, debug_name, divergence_retries) to every step of the march
instead of dropping them, and REFUSES the two that contradict it: the march decides
cold-vs-warm per step and manages its own warm-up, so zero_init_guess and picard now
raise rather than being silently ignored.

m8 - resolved by C3: the base dispatch's missing solve_kwargs only mattered for a
stress-history model on a non-overriding subclass, and the march refuses those now.

m9 - YieldHomotopyControl and SolveReport are exported from underworld3.systems.
Both are named in public docstrings (the control= parameter, solver.solve_report) and
were deep-import-only, which Charter §6 forbids.

m11 - the Darcy and TransientDarcy solve docstrings described the pre-flip default.

m12 - test_1055_solve_report -> test_1058 (1055 collided with the yield smoother
tests; 1056 was already taken twice, 1058 was genuinely free).

PARALLEL COVERAGE (Charter §11) - ptest_0202 exercises all four layers at np>1, which
nothing did before: has_solution and the tri-state resolution must agree across ranks
or the solve collectives they gate would split (a hang, not a wrong answer); the
homotopy march's branches all derive from collective SNES state; and the new
gmres/sor fixed-cost smoother has to work with the redundant-LU coarse solve, which
test_1014's serial option-string assertions cannot show. Green at np=2 and np=4, with
the SAME settled delta and step count at both -- the march is partition-independent.

Full level_1/tier_a gate: 480 passed.

Underworld development team with AI support from Claude Code
…n the homotopy capability

strainrate_inv_II_min was DECLARED on ViscoPlasticFlowModel but never used -- setting
it had no effect at all, while the elastic models (ViscoElasticPlastic,
TransverseIsotropicVEP) have always applied it. It now enters the plastic viscosity the
same way: eta_pl = tau_y / (2 (edot_II + edot_min)), capping eta_pl at
tau_y/(2 edot_min) and so bounding the viscosity CONTRAST -- the knob the xi-style
regularisation needs, and a prerequisite for testing whether xi can open cases the
delta-homotopy alone cannot reach. Default 0 = off, so the unregularised law is
unchanged.

Two CAPABILITY tests, so the homotopy cannot regress at the user level into "runs
without error but no longer rescues anything":

  * test_homotopy_rescues_a_solve_the_cold_start_cannot_do asserts BOTH halves on one
    problem -- the cold sharp solve must fail AND solve(homotopy=True) must converge,
    with a guard that the fixture is genuinely yielding (>20%) so the comparison cannot
    go vacuous.
  * test_rate_regularisation_is_wired_into_the_plastic_viscosity would have caught the
    ignored parameter above.

The fixtures use a 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 -- which is why every earlier viscoplastic fixture in this file sits at
0% yielding and why the homotopy went so long without being exercised on a case that
actually yields.

Underworld development team with AI support from Claude Code
A working backup of constitutive_models.py was swept into 39be0ab. It is a
duplicate of a source file, so it would shadow real hits in greps and pattern
scans forever.

Underworld development team with AI support from Claude Code
A hard viscoplastic Stokes solve does not fail, it grinds -- for hours, inside a
single Newton step. Two measurements say why nothing outside PETSc can stop it:
an iteration cap never fires because the grind happens within one outer Krylov
iteration, so the counter it watches does not advance; and a Python signal.alarm
never fires because Python runs signal handlers only between bytecodes while
control sits in PETSc's C code (a 90s alarm measured unfired at 10 minutes).

So the deadline lives in PETSc's own convergence tests. solver.guard(
wall_per_step=...) restarts a budget at each outer KSP solve -- once per Newton
step -- and checks it inside the fieldsplit sub-KSP tests, where iterations are
frequent enough for seconds of granularity. A guard exit is a report, not an
error: solve_report.deadline_expired distinguishes it from a real divergence.

The reason returned matters and is easy to get wrong. PETSc's KSPCheckSolve
deliberately exempts KSP_DIVERGED_ITS from marking the preconditioner failed, so
a guard returning DIVERGED_MAX_IT from a sub-KSP does nothing at all -- measured,
the outer solve absorbed 36 truncated velocity solves and reported CONVERGED.
DIVERGED_BREAKDOWN marks the sub-PC failed and surfaces as DIVERGED_LINEAR_SOLVE.

In parallel the expiry decision is reduced over the solver's communicator: wall
clocks are not synchronised, and a convergence test returning different verdicts
on different ranks deadlocks the next collective. ptest_0203 provokes this with
deliberately skewed per-rank clocks -- it passes at np=2 and np=4, and deadlocks
when the reduction is removed (confirmed at a 90s limit).

Also adds the work axis the guard needs. solve_report.ksp_its counts OUTER
Krylov iterations, which Eisenstat-Walker collapses to about one per Newton step;
measured on a small box, outer 1 against 285 multigrid cycles in the velocity
block. solve_report.sub records iterations and applications per fieldsplit block,
collected by monitors so it cannot influence a solve. Cost of leaving it on: 624
callbacks in a 280ms solve, no measurable change in wall time.

Known and documented rather than hidden: on a solver's first solve the fieldsplit
blocks do not exist until PCSetUp runs part-way through, so one preconditioner
application is uncovered and that solve's sub-counts are a lower bound
(SubSolveReport.complete says so). guard() raises on the rotated free-slip path,
which solves outside self.snes and cannot be reached.

Underworld development team with AI support from Claude Code
CI caught what the level_1/tier_a gate deselected: this tier_b test compares
solutions at rtol 1e-12 and 1e-10, but its diffusivity depends on T, so the
solve is nonlinear and its answer is only pinned to the SNES tolerance. With
zero_init_guess defaulting to warm, the second solve stops at a different point
in the same tolerance ball and the comparison fails by ~2e-5.

The test is about whether a constant is read from its constants[] slot or baked
as a C literal. Holding the initial guess fixed makes it measure that, and makes
it independent of the warm-start default rather than accidentally dependent on
one.

Underworld development team with AI support from Claude Code
…guard

Two independent reviewers were given the change and told to break it. Four of
the findings were real defects in the new code and are fixed here; two more were
claims in the docs that measurement did not support.

MAJOR - the deadline was disabled for the first preconditioner application of
EVERY solve, not just the first. begin_solve() left the clock at infinity and
only outer-KSP iteration 0 armed it -- but PETSc's default for GMRES is LEFT
preconditioning, and KSPInitialResidual applies the preconditioner BEFORE
iteration 0. Verified: for pc_side LEFT the sub-block monitors fire before the
outer monitor. So one full velocity solve plus one pressure solve ran outside
the budget on every solve. The clock now starts in begin_solve().

MAJOR - the premise that forced a hand-rolled convergence test was false.
petsc4py DOES expose composition: KSP.addConvergenceTest(prepend=True) runs a
test in FRONT of the native one, and returning ITERATING hands the decision back
to KSPConvergedDefault with its options-configured context intact. Verified: a
prepended always-ITERATING test gives identical iteration count and reason to an
uninstrumented solve. _default_convergence is deleted, and with it four silently
dropped KSPConvergedDefault behaviours (-ksp_converged_maxits, the
non-zero-initial-guess residual convention, the norm-type early return,
-ksp_min_it) and a bug reporting an infinite residual as CONVERGED_RTOL. It also
makes unguard() exact: the test stays installed and goes inert, which IS the
uninstrumented behaviour, rather than being swapped for a fresh default context.

MAJOR - the instrumentation held petsc4py references to a destroyed SNES's KSP,
PC and multigrid hierarchy until the next solve, so a rebuilding solver carried
two hierarchies at once. That is the leak BUGFIX(#157) was written to close. It
now releases them before the destroy.

MAJOR - guard()'s rotated-free-slip refusal was order-dependent: arming first and
adding the BC afterwards armed cleanly and then sat inert on a path that never
reaches the instrumentation. Re-checked at solve time.

MINOR - preonly and richardson outer KSPs are now skipped: both change what they
DO when a monitor is attached (gated on ksp->numbermonitors in PETSc), and
neither iterates, so instrumenting them was all cost and no information.

MEASURED, NOT FIXED - a reviewer argued a retry or a continuation stage would be
handed a fresh budget after an expiry. PETSc prevents this on its own: once the
sub-preconditioner is marked failed the next snes.solve bails before iterating,
so removing the expiry latch changed the cost by 34 ticks against 32. The latch
is kept because that protection is implicit, but it is belt-and-braces and no
test claims to prove otherwise.

Tests reworked against the second reviewer's demonstration that several of them
passed with the feature broken. Now covered: cost parity (not just answer parity)
between a guarded and an unguarded solve; the per-Newton-step meaning of
wall_per_step, on a nonlinear model, since one linear solve cannot distinguish
per-step from per-solve; the first-solve lower-bound contract; arming before the
first solve; both call orders of the rotated refusal. The suite is tier_b, not
tier_a: Charter 8 reserves tier A for tests that have been through use in anger.
The parallel test is registered in mpi_runner.sh -- it was unreachable by any
runner before, so the hazard the design calls severe had no automated coverage.

Underworld development team with AI support from Claude Code
Underworld development team with AI support from Claude Code
… PETSc

Every positive SNES reason was shifted by one slot, and the table claimed a
code 1 that PETSc does not define:

  reason 2  PETSc CONVERGED_FNORM_ABS       reported as CONVERGED_FNORM_RELATIVE
  reason 3  PETSc CONVERGED_FNORM_RELATIVE  reported as CONVERGED_SNORM_RELATIVE
  reason 4  PETSc CONVERGED_SNORM_RELATIVE  reported as CONVERGED_ITS
  reason 5  PETSc CONVERGED_ITS             reported as UNKNOWN_5

The direction matters. A solve that stopped on the STEP norm -- the weakest
criterion, and what a stalled viscoplastic solve reports -- was labelled
CONVERGED_ITS, while a genuine residual convergence was labelled
CONVERGED_SNORM_RELATIVE. Reading a difficulty report is how a continuation
driver decides whether a parameter station is reachable, so a mislabelled
convergence is the kind of error that quietly puts false rescues on a regime
map. Found while checking why a hard plastic station was reporting
CONVERGED_ITS.

Also adds DIVERGED_OBJECTIVE_DOMAIN (-13) and DIVERGED_OBJECTIVE_NANORINF
(-14), which were unmapped and surfaced as UNKNOWN_n, and corrects -4 to its
PETSc name DIVERGED_FUNCTION_NANORINF.

The table is hand-written on purpose -- the module must import without the
Cython extension -- so it can drift, and the existing test pinned it to ITSELF
(reason_string(2) == 'CONVERGED_FNORM_RELATIVE') rather than to the enum, which
is how this survived. The KSP table next to it WAS checked against petsc4py and
was correct. test_snes_reason_table_matches_petsc now checks both directions:
every label matches the enum, and every reason PETSc can return is mapped.

Underworld development team with AI support from Claude Code
…in sits on

Smoothing a corner by delta moves the whole curve, so each soft-min family has one
free constant. Until now it was fixed implicitly, and differently per family, with
the consequence buried in the algebra: both families ended up sitting BELOW exact
Min at and above the yield point. For a homotopy that is the wrong direction — an
entry problem weaker than the sharp problem is not an easier version of it, it is a
different and softer problem.

`yield_anchor` names the choice and applies it to both families:

  "onset" (default, the historical law) pins the unyielded limit f -> 0. The sqrt
    family then undershoots for f < 2 — a third under-stressed at nominal yield —
    and the power mean degenerates to the p-norm, whose viscosity falls away
    toward zero as delta grows.

  "yield" pins f = 1, so tau/tau_y is exactly 1 at the yield point for every delta
    and the curve is >= exact Min everywhere. The sqrt family reaches it with the
    offset delta/2; the power mean needs no offset at all — averaging its two terms
    is what makes it a generalised MEAN rather than a p-NORM, and a generalised
    mean returns the common value when its arguments are equal.

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
delta is O(1) where the sqrt family's is O(10) — the two deltas are not the same
parameter, so the entry now comes from the family (`_YIELD_HOMOTOPY_DELTA0`) rather
than from a caller carrying one number across both.

`yield_continuation` gains `smoother=` and `anchor=` so a march can pick both, and
refuses them when handed a ready-made control rather than silently ignoring them.

The default stays "onset": flipping it changes results for every existing
yield_mode="softmin" user, and that is a maintainer decision, not a side effect of
naming the knob. Tests that exercise the old side now set the anchor explicitly
instead of relying on the default, and test_1059 asserts the OLD anchor FAILS the
new behaviour so it cannot pass against a no-op.

Underworld development team with AI support from Claude Code
Copilot AI review requested due to automatic review settings July 29, 2026 08:01

Copilot AI left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

@lmoresi

lmoresi commented Jul 29, 2026

Copy link
Copy Markdown
Member Author

Adversarial review

Read against the diff, not the description. Four findings, none blocking; the two worth acting on before merge are R1 and R4.

R1 — yield_anchor rebuilds the δ atom, so an already-built YieldHomotopyControl holds a stale handle.
The setter nulls both _yield_offset_expr and _yield_softness_expr and calls _reset(). That is correct for the value_get_yield_softness rebuilds from self._yield_softness, which the setter does not touch — but it creates a new expression object. control.delta is captured at construction, so

control = cm._yield_homotopy_control()
cm.yield_anchor = "yield"          # rebuilds the atom
control.delta                      # now detached from the model

leaves the diagnostics handle pointing at an orphan. set_delta is safe (it goes through the yield_softness property, which re-fetches), and _yield_homotopy_control sets the anchor before it reads delta=, so the driver path is correct — and yield_continuation already refuses anchor= alongside a ready-made control, which closes the obvious footgun. What is left is the hand-rolled path. Cheapest honest fix is a line in the YieldHomotopyControl.delta docstring saying it is a snapshot; the better one is to make it a property that re-reads.

R2 — the power mean's 2^δ background stiffening deserves a cross-reference, not just a bound.
The description and docstrings both state it, and the algebra checks out (as eta_pl → ∞, softmin → a^-s/2, so η → η_ve · 2^(δ+0.001)). But at the family's own entry δ0 = 1 that is a factor 2 on the entire unyielded domain — which makes the entry problem a materially different problem, not a rounded version of the target. That is very likely relevant to #473's first defect (the cold power-mean entry at delta0 = 1 failing at all four notch regime points). Worth a See also between them so the next person does not diagnose it twice.

R3 — _yield_offset_for_anchor reads self._yield_softness_expr directly.
It is only ever called from inside _get_yield_softness, on the line after that attribute is assigned, so it cannot currently see None. But the coupling is implicit and the function name does not suggest it. Taking δ as an argument would make it total.

R4 — no test that the anchor survives a δ update.
This is the coverage gap that matches the risk R1 describes. The setter rebuilds the offset atom from the anchor's formula, and the whole design rests on that atom continuing to track δ symbolically afterwards. Every test in test_1059 sets the anchor and then evaluates; none does

set anchor → set δ → check the pin → set δ again → check the pin again

which is exactly the sequence a march performs. A no-op rebuild bug (offset frozen at its construction-time δ) would pass the current suite and fail every homotopy. Cheap to add and it guards the one mechanism that is genuinely subtle here.

Checked and clean:

  • The /2 does what it claims: at a = b, ((2a^-s)/2)^(-1/s) = a exactly.
  • Sides are right in both directions — the generalised mean M_{-s} ≥ min, the p-norm ≤ min.
  • The sqrt family's δ/2 offset gives g → 1/2 as δ → ∞, hence the stated factor-2 bound.
  • Ordering in _yield_homotopy_control (smoother → anchor → read δ) is correct; reversing the first two would lose the power mean's δ bump.
  • test_1055's new explicit anchor="onset" lines are the right move — they stop those assertions from silently inverting if the default is ever flipped, which is precisely the open maintainer question.

@lmoresi

lmoresi commented Jul 30, 2026

Copy link
Copy Markdown
Member Author

Heads-up from the solver configuration audit: unpushed commits, and a collision with #471

Two coordination points found while auditing the solver configuration surface against development (report: docs/reviews/2026-07/solver-configuration-reachability-audit.md on branch docs/solver-config-audit).

1. Three commits are local only

feature/regime-diagram carries three commits that are not on origin:

954ff4f3 fix(solver): make Eisenstat-Walker switchable
5ad68b2f docs(solvers): the Stokes fieldsplit is two nested Krylov loops
f1555b3a fix(solver): make the Stokes solver settings reachable (#477, D-22/D2)

So this PR does not show the #477 fix. Auditing development or reading the PR list rediscovers a solved problem — which is exactly what happened here before we checked the local branch. Worth pushing so the fix is visible to anyone else working this surface.

Verified what it covers: snes_max_it and the two fieldsplit_*_ksp_rtol keys. It does not cover the MG bundle keys (fieldsplit_velocity_mg_levels_ksp_type / _max_it / mg_coarse_pc_type), which are clobbered at build time by _apply_preconditioner_options rather than at solve time — a different mechanism. That half is fixed in #471 (commit 5a390efe).

2. This PR and #471 both edit the same stale-key list, and #471 subsumes it

This PR patches the hand-written GAMG stale-key clear in _apply_preconditioner_options:

-                        "mg_levels_pc_type", "mg_coarse_pc_type",
+                        "mg_levels_pc_type", "mg_levels_ksp_norm_type",
+                        "mg_coarse_pc_type",

#471 deletes that list entirely and derives it — every key a sibling bundle sets that this one does not — and the derived GAMG stale set already contains mg_levels_ksp_norm_type. So the fix here is correct and already covered by the derivation.

They will conflict textually. Whichever lands second should take #471's derivation rather than merging the two lists, otherwise the hand-written list comes back and so does the class of bug that produced this patch.

3. Adjacent findings routed here

Filed as #483: Stokes and Stokes_Constrained disagree about which options are reachable, in opposite directions (ksp_rtol settable on one and not the other; the fieldsplit rtols the reverse) — plus snes_rtol, which _reassert_outer_tolerances() deliberately re-writes, so it is owned by design but still looks like a settable knob. Both belong on this branch, next to the fix they are adjacent to.

Underworld development team with AI support from Claude Code

lmoresi added a commit that referenced this pull request Jul 30, 2026
…483/#484

F-1 is fixed in #471 (5a390ef). F-3/F-4 filed as #483, F-6/F-7/F-8 as #484, F-5
raised on #475. Adds the two lessons from implementing F-1: the ownership record
must be keyed by the globally-qualified option name (an unqualified record made the
bundle silently stop applying on the custom-P route), and the defaults arm is what
caught that — the second ownership regression in two sessions caught by the arm
that was supposed to be unchanged.

Underworld development team with AI support from Claude Code
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants