diff --git a/docs/developer/design/nonlinear-solver-homotopy-warmstart.md b/docs/developer/design/nonlinear-solver-homotopy-warmstart.md index 0ce8dda9..27d58e0e 100644 --- a/docs/developer/design/nonlinear-solver-homotopy-warmstart.md +++ b/docs/developer/design/nonlinear-solver-homotopy-warmstart.md @@ -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 diff --git a/docs/developer/subsystems/solvers.md b/docs/developer/subsystems/solvers.md index e03c9309..167bd935 100644 --- a/docs/developer/subsystems/solvers.md +++ b/docs/developer/subsystems/solvers.md @@ -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 diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 4be2d14c..d7db6e05 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -958,25 +958,86 @@ def _get_yield_softness(self): 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). + # The offset fixes WHERE the smoothed curve is pinned to the exact law. + # 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). See `yield_anchor`. 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", + self._yield_offset_for_anchor(), + "Yield soft-min offset; tracks δ so the chosen anchor is exact", ) else: self._yield_softness_expr.sym = sympy.Float(delta_value) return self._yield_softness_expr + def _yield_offset_for_anchor(self): + r"""The sqrt soft-min offset that pins the chosen anchor, symbolic in δ.""" + delta = self._yield_softness_expr + if getattr(self, "_yield_anchor", "onset") == "yield": + return delta / 2 + return (-1 + sympy.sqrt(1 + delta**2)) / 2 + def _get_yield_offset(self): - """The onset-offset constant atom (lazily created alongside δ).""" + """The offset constant atom (lazily created alongside δ).""" if getattr(self, "_yield_offset_expr", None) is None: self._get_yield_softness() return self._yield_offset_expr + @property + def yield_anchor(self): + r"""Which point of the ``"softmin"`` yield law is pinned to the exact ``Min``. + + A soft-min rounds the yield corner by δ, and rounding a corner moves the whole + curve — so the family has one free constant, fixed by choosing WHICH point the + smoothed law must reproduce exactly. Both choices apply to both smoother + families, and both recover ``Min`` as :math:`\delta \to 0`; they differ only in + which SIDE of the exact law the curve sits on in between. Writing + :math:`f = \eta_{ve} / \eta_{pl}` for the overstress ratio (so :math:`f = 1` is + the yield point and :math:`f < 1` is unyielded): + + ``"onset"`` (default, historical) + Pins the unyielded limit :math:`f \to 0`, where :math:`\eta = \eta_{ve}` + exactly. The curve then sits BELOW the exact law at and above the yield + point — the sqrt family undershoots for :math:`f < 2` (at + :math:`\delta = 64`, :math:`\eta/\eta_{ve} = 0.80` at :math:`f = 0.5` where + the exact law gives 1.0, and the yield point itself is a third + under-stressed, :math:`\tau/\tau_y = 0.67`), and the power mean degenerates + to the p-norm :math:`(\eta_{ve}^{-s} + \eta_{pl}^{-s})^{-1/s}`, which falls + away to :math:`\eta \to 0` as δ grows. So neither family approaches the + yield surface from above under this anchor. That matters whenever the + problem's overstress ratio is O(1), because then the smoothed problem is + WEAKER than the sharp one it is supposed to be an easy version of. + + ``"yield"`` + Pins :math:`f = 1`, so :math:`\tau/\tau_y = 1` exactly at the yield point + for every δ and the curve is :math:`\ge` the exact law everywhere — a + genuine approach from above, which is the safe direction for a homotopy + entry. The sqrt family gets there with the offset :math:`\delta/2`; the + power mean needs no offset at all, since averaging the two terms, + :math:`\eta = ((\eta_{ve}^{-s} + \eta_{pl}^{-s})/2)^{-1/s}`, makes it a + generalised mean, and a generalised mean returns the common value when its + arguments are equal. The cost is stiffened unyielded material, bounded by a + factor 2 for the sqrt family and by :math:`2^{\delta}` for the power mean, + both decaying to 1 as :math:`\delta \to 0`. The power-mean bound is why its + entry δ is O(1) and not O(10). + """ + return getattr(self, "_yield_anchor", "onset") + + @yield_anchor.setter + def yield_anchor(self, value): + if value not in ("onset", "yield"): + raise ValueError( + f"yield_anchor must be 'onset' or 'yield', got {value!r}") + self._yield_anchor = value + # The sqrt family's offset atom holds the anchor's FORMULA, so it must be + # rebuilt, not merely revalued — and _get_yield_softness only builds it + # alongside δ. (The power mean carries the anchor in _combine_yield itself, so + # for that family the _reset is the whole of the update.) + self._yield_offset_expr = None + self._yield_softness_expr = None + self._reset() + 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``: @@ -984,14 +1045,18 @@ def _combine_yield(self, eta_ve, eta_pl): - ``"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). + ``self.yield_smoother`` (``"sqrt"`` or ``"powermean"``) and on the side of + exact ``Min`` chosen by ``self.yield_anchor``. Which side is the property + that matters for a homotopy, and it belongs to the ANCHOR, not the family: + ``"onset"`` puts both families below ``Min`` near the yield point, + ``"yield"`` puts both on or above it. 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). + lands within a factor :math:`2^{\pm(\delta + 0.001)}` of ``Min``, i.e. 0.07 % + at ``δ = 0``. 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 @@ -1008,7 +1073,7 @@ def _combine_yield(self, eta_ve, eta_pl): 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 + # Soft-min of order -s 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)) @@ -1021,7 +1086,14 @@ def _combine_yield(self, eta_ve, eta_pl): # 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) + softmin = a ** (-s) + b ** (-s) + if self.yield_anchor == "yield": + # Averaging the two terms is what makes this a power MEAN rather than + # a p-NORM, and it is the whole of the "yield" anchor for this family: + # a generalised mean returns the common value when its arguments are + # equal, so η = η_ve = η_pl exactly at f = 1 — no offset atom needed. + softmin = softmin / 2 + return N * softmin ** (-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() @@ -1037,19 +1109,30 @@ def yield_smoother(self): 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. + ``g = 1 + ½(f−1+√((f−1)²+δ²)) − offset``, ``f = η_ve/η_pl``. Exact + ``Min`` at ``δ=0``. Its smoothing authority saturates: as ``δ → ∞``, + ``η → 2 η_ve/(f+2)``, so the most it can move a point is a factor + ``2f/(f+2)`` — which tends to 1 as ``f → 1``. On a problem whose + overstress ratio is O(1) this family therefore has very little room to + build an easier problem, however large δ is. + - ``"powermean"``: the soft-min of order ``−s``, + ``η_eff = (η_ve^(−s) + η_pl^(−s))^(−1/s)`` with ``s = 1/(δ + 0.001)`` + (``s→∞`` ⇒ ``Min``). Unlike the sqrt family it does not saturate — it + keeps moving with δ — which is what makes it usable where the sqrt family + has no range. + + ``δ`` is NOT the same parameter in the two families and their schedules do + not transfer: for the power mean ``s = 1/(δ + 0.001)``, so ``δ = 1`` is + already the harmonic mean and useful entries are ``δ ≤ 1``; for the sqrt + family δ is a percentage stress deviation and a generous entry is O(10). + Take the entry from ``_yield_homotopy_control``, which asks the family. + + Which SIDE of ``Min`` a family sits on is set by ``yield_anchor``, not by the + family — see that property. 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*. + (``s≈1``, the parameter-free harmonic mean) since ``δ=0`` is the sharp + hard-``Min`` limit it only *approaches*. """ return getattr(self, "_yield_smoother", "sqrt") @@ -1074,7 +1157,44 @@ def yield_smoother(self, value): # 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): + @property + def viscosity_min_rounding(self): + r"""Rounding scale :math:`\epsilon` for the VISCOSITY floor, in viscosity units. + + A floor is a corner, and at :math:`\epsilon = 0` the smooth-max expression is + algebraically exactly ``Max`` — non-differentiable, which the consistent-Newton + tangent cannot cope with (it dies at ``nl=0, DIVERGED_LINEAR_SOLVE``). Without + this property the only rounding available came from ``yield_mode``: zero under + ``"min"``, and :math:`\delta\,|floor|` under the smooth modes, so the smoothing + **vanished as δ → 0** — exactly where a yield homotopy lands. + + The floor is a different corner from the yield surface and needs a scale of its + own. Set this to a small fraction of the floor (a few per cent is usually ample) + and the cutoff becomes usable with ``consistent_jacobian=True`` at any δ, + including δ = 0. + + ``None`` (the default) keeps the historical ``yield_mode``-derived behaviour, so + nothing changes for existing models. + + Notes + ----- + This matters for more than differentiability. A viscosity floor bounds how weak + yielded material can become, and therefore bounds the viscosity CONTRAST the + solution can develop — which is to say it bounds how localised the solution can + be. A floor that is relaxed toward zero is a continuation from a diffuse, + uniquely-solvable viscous problem toward the localised (near rigid-plastic) + limit, tracking one branch as it sharpens. That is a solution-SELECTION + homotopy, not a smoothness fix, and it is why the floor needs to be usable + independently of δ. + """ + return getattr(self, "_viscosity_min_rounding", None) + + @viscosity_min_rounding.setter + def viscosity_min_rounding(self, value): + self._viscosity_min_rounding = value + self._reset() + + def _apply_floor(self, value, floor, rounding=None): r"""Impose the lower bound :math:`value \ge floor`. In ``yield_mode="min"`` this is the exact hard ``sympy.Max(value, floor)`` — @@ -1101,8 +1221,9 @@ def _apply_floor(self, value, floor): # 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 + if rounding is None: + 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 @@ -1112,7 +1233,15 @@ def _apply_floor(self, value, floor): # solve fails outright (DIVERGED_LINEAR_SOLVE). _yield_homotopy_tangent = True - def _yield_homotopy_control(self): + #: Entry δ per soft-min family. These are NOT interchangeable numbers: for the + #: power mean the sharpness is s = 1/(δ + 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 δ handed to the power mean asks for a law a factor 2^δ + #: away from Min — 65000 at δ=16 — which is why the entry has to come from the + #: family, not the caller. + _YIELD_HOMOTOPY_DELTA0 = {"powermean": 1.0, "sqrt": 16.0} + + def _yield_homotopy_control(self, smoother=None, anchor=None): """Put this model in its smooth (δ-parameterised) yield mode and describe how to march it. @@ -1122,22 +1251,52 @@ def _yield_homotopy_control(self): ``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. + Defaults to the power-mean 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. + + **Which family suits a given problem is not settled by theory** — it is set by + the problem's overstress ratio :math:`f = \\eta_{ve}/\\eta_{pl}`, which is cheap + to measure on the viscous seed before choosing. The sqrt family's smoothing + saturates at a factor :math:`2f/(f+2)`, so on a problem where f is O(1) it can + barely build an easier problem at all, however large δ is; the power mean does + not saturate and has range there. Where f is large the position reverses. + + The side the regularised law sits on is set by ``anchor=``, not by the family. + An entry problem should be no WEAKER than the sharp problem it precedes, so + ``anchor="yield"`` is the safe direction for both families — see + :attr:`yield_anchor`. + + Parameters + ---------- + smoother : {"powermean", "sqrt"}, optional + Soft-min family to march. ``None`` keeps the default (power mean). + anchor : {"onset", "yield"}, optional + Which point the smoothed law reproduces exactly. ``None`` leaves the + model's current :attr:`yield_anchor` alone. Returns ------- YieldHomotopyControl ``set_delta`` (model-owned setter for δ), ``tangent`` (the - ``consistent_jacobian`` value to use), and ``delta`` (the ``constants[]`` - atom itself, for diagnostics). + ``consistent_jacobian`` value to use), ``delta`` (the ``constants[]`` + atom itself, for diagnostics) and ``delta0`` (the family's entry δ). """ from underworld3.systems.yield_continuation import YieldHomotopyControl + if smoother is None: + smoother = "powermean" + if smoother not in self._YIELD_HOMOTOPY_DELTA0: + raise ValueError( + f"smoother must be one of {sorted(self._YIELD_HOMOTOPY_DELTA0)}, " + f"got {smoother!r}" + ) + self.yield_mode = "softmin" - self.yield_smoother = "powermean" + self.yield_smoother = smoother + if anchor is not None: + self.yield_anchor = anchor # 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 @@ -1155,6 +1314,7 @@ def set_delta(value): set_delta=set_delta, tangent=self._yield_homotopy_tangent, delta=self._get_yield_softness(), + delta0=self._YIELD_HOMOTOPY_DELTA0[smoother], ) @@ -1347,7 +1507,8 @@ def viscosity(self): if inner_self.shear_viscosity_min.sym != -sympy.oo: self._plastic_eff_viscosity._sym = self._apply_floor( - effective_viscosity, inner_self.shear_viscosity_min + effective_viscosity, inner_self.shear_viscosity_min, + rounding=self.viscosity_min_rounding, ) else: @@ -2024,6 +2185,15 @@ def viscosity(self): # BDF-2 Jacobian. Those modes are already smooth and bounded. if inner_self.shear_viscosity_min.sym != -sympy.oo: + rounding = self.viscosity_min_rounding + if rounding is not None: + # An explicit rounding scale makes the cutoff differentiable, which is + # what the skip below existed to avoid needing: there is no outer Max to + # nest, so the floor can be honoured in the smooth yield modes too. + return self._apply_floor( + effective_viscosity, inner_self.shear_viscosity_min, + rounding=rounding, + ) if self.is_viscoplastic and self._yield_mode in ("harmonic", "softmin"): return effective_viscosity else: diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index a5c36ee8..fb3e7f82 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -6389,6 +6389,15 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): """ return self._tolerance + #: Inner-solve tolerance margins, as a fraction of the outer tolerance. The inner + #: solves are deliberately inexact (Citcom; Moresi & Solomatov 1995) — which is why + #: the sub-blocks are flexible Krylov methods — but the inexactness must stay well + #: BELOW the tolerance demanded of the outer solve. These factors are that margin. + #: Their existence is principled; their size is inherited convention, so they are + #: DEFAULTS a user may override rather than values this property owns outright. + _INNER_RTOL_MARGIN = {"fieldsplit_pressure_ksp_rtol": 0.1, + "fieldsplit_velocity_ksp_rtol": 0.033} + @tolerance.setter def tolerance(self, value): self._tolerance = value @@ -6397,8 +6406,65 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): self.petsc_options["snes_ksp_ew_version"] = 3 self.petsc_options["ksp_atol"] = self._tolerance * 1.0e-6 - self.petsc_options["fieldsplit_pressure_ksp_rtol"] = self._tolerance * 0.1 # rule of thumb - self.petsc_options["fieldsplit_velocity_ksp_rtol"] = self._tolerance * 0.033 + + # Setting the tolerance re-derives the inner margins from it — that is this + # property's job, and a user who changes the tolerance expects it. What must NOT + # happen is `solve()` re-deriving them on every call: it did (pyx `solve()` + # round-trips `self.tolerance` immediately before `setFromOptions()`), which + # overwrote any user value between it being set and PETSc reading it and made + # both documented options silently unreachable (#477). `solve()` now re-asserts + # only the outer keys, via `_reassert_outer_tolerances`. + for key, margin in self._INNER_RTOL_MARGIN.items(): + self.petsc_options[key] = self._tolerance * margin + + def _resolve_snes_max_it(self, default): + """The nonlinear iteration cap for this solve, honouring a user-set + ``snes_max_it``. + + ``solve()`` pushes this option before every solve, so it has to be able to tell + its OWN previous push from a value the user set — otherwise the first solve makes + the option permanently unreachable. Call once, before this solve pushes anything. + """ + pushed = getattr(self, "_snes_max_it_pushed", None) + user = getattr(self, "_snes_max_it_user", None) + if self.petsc_options.hasName("snes_max_it"): + try: + current = int(self.petsc_options.getInt("snes_max_it")) + except Exception: + return default + if pushed is None or current != pushed: + # Never pushed by us, or the user has moved it since. Latch: from here on + # the option is theirs, because the next solve will read back OUR push of + # THEIR value and would otherwise mistake it for our own default. + self._snes_max_it_user = current + return current + if user is not None: + return user + return default + + def _push_snes_max_it(self, value): + """Push ``snes_max_it`` and remember what was pushed, so a later solve can tell + this solver's own value apart from a user override.""" + value = int(value) + self.petsc_options.setValue("snes_max_it", value) + self._snes_max_it_pushed = value + + def _reassert_outer_tolerances(self): + """Re-push the OUTER tolerance keys before a solve, leaving the inner margins be. + + `solve()` may have changed `snes_max_it` and the SNES type for a Picard warm-up, + so the outer settings are re-asserted before the real solve. The sub-block rtols + are deliberately excluded: they belong to whoever set them last, which may be the + user (#477). Overwriting them here is what made them unsettable.""" + self.petsc_options["snes_rtol"] = self._tolerance + self.petsc_options["ksp_atol"] = self._tolerance * 1.0e-6 + # The Eisenstat-Walker flags are NOT re-asserted here. solve() never changes + # them, so they stay in the options DB from the `tolerance` setter and + # setFromOptions picks them up regardless — while re-asserting would make + # `snes_ksp_ew` impossible to switch OFF, which is the same defect as #477 on a + # knob that matters: EW re-picks the outer KSP rtol every Newton step and + # OVERRIDES ksp_rtol, so "how hard is the linear solve actually being asked to + # work" is not answerable without being able to disable it. @property @@ -8741,8 +8807,14 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): 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. + ``smoother``, ``delta0``, ``down``, ``dmin``, ``entry_maxit``, + ``step_maxit``, ``retries``. All are defaulted; tuning them is optional. + ``smoother`` picks the soft-min family — ``"powermean"`` (default, + approaches the yield surface from below) or ``"sqrt"`` (from above). Which + gives the better cold entry is problem-dependent, so it is worth trying both + when a march will not start. Leave ``delta0`` unset unless you have a + reason: each family supplies its own entry, and the two δ are not the same + parameter. Returns ------- @@ -8863,15 +8935,12 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): UW_DMSetTime(_time_dm_stokes.dm, t_nd) # Keep a record of these set-up parameters - tolerance = self.tolerance snes_type = self.snes.getType() - # NOTE: hardcoded iteration cap. Unlike tolerance/snes_type above (read - # from current state), this value is pushed to petsc_options before the - # final solve below, silently clobbering any user-set snes_max_it on - # every solve. Respecting the user's setting is a behaviour change - # deferred to the benchmarked D2 sub-wave — see - # docs/reviews/2026-07/REMEDIATION-WORKLIST.md rows D-22/D2 (ruling D18). - snes_max_it = 50 + # Resolved ONCE, before any of this solve's own pushes, so the resolver compares + # against the previous solve's push rather than this one's. (Was a hardcoded 50 + # pushed unconditionally, which made `snes_max_it` unreachable — worklist rows + # D-22/D2, ruling D18; same failure shape as the sub-block rtols in #477.) + snes_max_it = self._resolve_snes_max_it(50) self.mesh.update_lvec() self.dm.setAuxiliaryVec(self.mesh.lvec, None) @@ -8888,7 +8957,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): print(f"SNES pre-solve - non-zero initial guess", flush=True) - self.petsc_options.setValue("snes_max_it", 0) + self._push_snes_max_it(0) self.snes.setType("nrichardson") self.snes.setFromOptions() # PETSc may rebuild operator state after setFromOptions(), so reattach @@ -8931,8 +9000,8 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): # Picard solves if requested if picard != 0: - self.petsc_options.setValue("snes_max_it", abs(picard)) - self.tolerance = tolerance + self._push_snes_max_it(abs(picard)) + self._reassert_outer_tolerances() self.snes.atol = self.atol self.snes.setType("nrichardson") self.snes.setFromOptions() @@ -8944,9 +9013,9 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): # warmup above was taken: restore the configured SNES type and # tolerances, then solve. self.snes.setType(snes_type) - self.tolerance = tolerance + self._reassert_outer_tolerances() self.snes.atol = self.atol - self.petsc_options.setValue("snes_max_it", snes_max_it) + self._push_snes_max_it(snes_max_it) self.snes.setFromOptions() self._attach_stokes_nullspace() # Custom geometric-MG prolongation on the velocity block (if registered diff --git a/src/underworld3/systems/yield_continuation.py b/src/underworld3/systems/yield_continuation.py index ff162619..b6feeb10 100644 --- a/src/underworld3/systems/yield_continuation.py +++ b/src/underworld3/systems/yield_continuation.py @@ -40,17 +40,26 @@ class YieldHomotopyControl: yield, the frozen/Picard tangent for the elastic (VEP) models. delta The δ ``constants[]`` atom itself, for diagnostics. + delta0 + The entry δ this family should START from. Model-owned because δ does NOT mean + the same thing in the two soft-min families: for the power mean the sharpness is + ``s = 1/(δ + 0.001)``, so δ ≤ 1 and δ = 1 IS the harmonic mean, while for the + sqrt family δ is a percentage stress deviation and a generous entry is O(10). + Carrying one number across both families is a category error. """ set_delta: Callable[[float], None] tangent: Any delta: Optional[Any] = None + delta0: float = 1.0 def yield_continuation( solver, control=None, - delta0=1.0, + smoother=None, + anchor=None, + delta0=None, down=0.5, dmin=1.0e-3, entry_maxit=30, @@ -74,10 +83,9 @@ def yield_continuation( 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. + Starting δ is deliberately **large**, where the power-mean soft-min stays 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 ---------- @@ -86,8 +94,26 @@ def yield_continuation( 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. + smoother : {"powermean", "sqrt"}, optional + Which soft-min FAMILY to march. Default (``None``) leaves the model's own + choice, which is the power mean. The families differ in how much room δ gives + them: the sqrt law's smoothing saturates at a factor :math:`2f/(f+2)` in the + overstress ratio :math:`f = \eta_{ve}/\eta_{pl}`, so where f is O(1) it can + barely build an easier problem however large δ is, while the power mean does not + saturate. Where f is large the position reverses. f is cheap to measure on the + viscous seed, so choose the family from it rather than by reputation. + anchor : {"onset", "yield"}, optional + Which point of the smoothed law is pinned to the exact one — see + ``ViscousFlowModel.yield_anchor``. Applies to BOTH families. Default (``None``) + leaves the model's own choice. ``"yield"`` keeps the stress exact at the yield + point for every δ and puts the law on or above the exact one throughout; + ``"onset"`` (the historical default) is exact on the unyielded branch but sits + BELOW the exact law at and above the yield point, which makes the entry problem + WEAKER than the sharp problem it is supposed to lead to. + delta0 : float, optional + Starting (smooth) δ. Defaults to the family's own entry (``control.delta0``), + because δ is NOT the same parameter in the two families — see + :class:`YieldHomotopyControl`. Do not carry one number across both. down : float Nominal multiplicative step, :math:`0 < down < 1` (δ ← δ·down per success). Adapted during the march as described above. @@ -118,8 +144,6 @@ def yield_continuation( """ 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) @@ -129,7 +153,19 @@ def yield_continuation( "homotopy (supports_yield_homotopy). Got " f"{type(cm).__name__ if cm is not None else None}." ) - control = cm._yield_homotopy_control() + control = cm._yield_homotopy_control(smoother=smoother, anchor=anchor) + elif smoother is not None or anchor is not None: + raise ValueError( + "smoother=/anchor= configure the soft-min when this driver builds the " + "control; they cannot override a control that was passed in ready-made. " + "Set them on the model before building the control instead." + ) + + # Each family has its own natural entry, and they are not interchangeable numbers. + if delta0 is None: + delta0 = control.delta0 + if not delta0 > 0.0: + raise ValueError(f"delta0 must be positive, got {delta0}") # 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. diff --git a/tests/test_1055_yield_smoother.py b/tests/test_1055_yield_smoother.py index b1383a80..580b29df 100644 --- a/tests/test_1055_yield_smoother.py +++ b/tests/test_1055_yield_smoother.py @@ -78,8 +78,12 @@ def test_delta0_is_exact_min(): @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/δ→∞).""" + """Power-mean smooth-min UNDER THE ONSET ANCHOR: ≤ Min everywhere (no over-yield), + overflow-safe on geodynamic ranges, and → Min as δ→0 (s=1/δ→∞). + + Which side of Min the law sits on belongs to ``yield_anchor``, not to the family, so + the anchor is set explicitly here: under ``"yield"`` this same family is ≥ Min by + construction and every assertion below would (correctly) invert.""" 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) @@ -88,6 +92,7 @@ def test_powermean_smoother_undershoots_min(): s.constitutive_model = cm cm.yield_mode = "softmin" # dev: soft-min family lives under "softmin" + cm.yield_anchor = "onset" # the law under test; not left to the default 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" @@ -169,8 +174,9 @@ def test_dp_model_smoother_optin(): 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 + # opt into the power-mean homotopy: undershoots Min under the onset anchor cm.yield_mode = "softmin" + cm.yield_anchor = "onset" # the side under test; not the default cm.yield_smoother = "powermean" # bumps δ→1 assert cm.yield_softness == 1.0 for eta_ve, eta_pl in [(1.0, 2.0), (5.0, 0.3), (1e25, 1e21)]: diff --git a/tests/test_1059_yield_anchor.py b/tests/test_1059_yield_anchor.py new file mode 100644 index 00000000..16765d85 --- /dev/null +++ b/tests/test_1059_yield_anchor.py @@ -0,0 +1,116 @@ +"""Which point of the soft-min yield law is pinned to exact Min (`yield_anchor`). + +A soft-min rounds the yield corner by delta, and rounding a corner moves the whole +curve -- so the family has one free constant, fixed by choosing which point the +smoothed law must reproduce exactly. That choice decides which SIDE of the exact law +the curve sits on, and for a homotopy the side is the property that matters: an entry +problem that is WEAKER than the sharp problem is not an easier version of it. + +The historical anchor pins the unyielded limit f -> 0. Both families then sit BELOW +exact Min at and above the yield point -- the sqrt family undershoots for f < 2 (a +third under-stressed at nominal yield), and the power mean degenerates to the p-norm, +whose viscosity falls away towards zero as delta grows. + +Anchoring at the yield point f = 1 instead puts both families on or above exact Min +everywhere. The sqrt family gets there with the offset delta/2; the power mean needs no +offset, because averaging its two terms makes it a generalised mean and a generalised +mean returns the common value when its arguments are equal. + +The tests that pin the new behaviour also assert the OLD anchor FAILS it, so they +cannot pass against a no-op implementation. +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw +from underworld3.function.expressions import unwrap + + +def _model(smoother="sqrt"): + m = uw.meshing.UnstructuredSimplexBox(minCoords=(0,0), maxCoords=(1,1), cellSize=0.5, qdegree=2) + v = uw.discretisation.MeshVariable('U', m, m.dim, degree=2) + p = uw.discretisation.MeshVariable('P', m, 1, degree=1) + s = uw.systems.Stokes(m, velocityField=v, pressureField=p) + s.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel + c = s.constitutive_model; c.yield_mode='softmin'; c.yield_smoother=smoother + return c + +def eta(c, f, d): + """eta_eff/eta_ve at overstress f (eta_ve=1 => eta_pl=1/f).""" + c.yield_softness = d + e = c._combine_yield(sympy.Float(1.0), sympy.Float(1.0/f) if f>0 else sympy.oo) + return float(sympy.N(unwrap(e, keep_constants=False))) + +FS = (0.25, 0.5, 1.0, 2.0, 4.0, 10.0) + +# delta is NOT the same parameter in the two families: for the power mean +# s = 1/(delta + 0.001), so delta = 1 is already the harmonic mean and the law sits a +# factor 2^delta from Min. Each family is therefore probed over its own useful range. +DELTAS = {"sqrt": (1.0, 16.0, 64.0), "powermean": (0.1, 0.25, 1.0)} + +@pytest.mark.level_1 +@pytest.mark.tier_a +def test_delta_zero_is_exact_min_for_both_anchors(): + """Only the sqrt family is exactly Min at delta=0; the power mean's sharpness + saturates at s=1000, so it is checked against its own 2^0.001 bound.""" + c = _model("sqrt") + for anchor in ("onset", "yield"): + c.yield_anchor = anchor + for f in FS: + assert abs(eta(c, f, 0.0) - 1.0/max(1.0, f)) < 1e-12, (anchor, f) + c = _model("powermean") + for anchor in ("onset", "yield"): + c.yield_anchor = anchor + for f in FS: + # The bound is attained, not approached: 2^0.001 - 1 = 6.93e-4 exactly. + assert abs(eta(c, f, 0.0)/(1.0/max(1.0, f)) - 1.0) < 1e-3, (anchor, f) + +@pytest.mark.level_1 +@pytest.mark.tier_a +@pytest.mark.parametrize("smoother", ("sqrt", "powermean")) +def test_yield_anchor_keeps_onset_stress_exact(smoother): + """tau/tau_y = f*eta/eta_ve = 1 at f=1 for every delta. The onset anchor does NOT.""" + c = _model(smoother); c.yield_anchor = "yield" + for d in DELTAS[smoother]: + assert abs(eta(c, 1.0, d) - 1.0) < 1e-10, d + # The defect being fixed: sqrt is a third under-stressed at nominal yield, and the + # p-norm the power mean falls back to is HALF the correct viscosity at delta=1. + c.yield_anchor = "onset" + assert abs(eta(c, 1.0, DELTAS[smoother][-1]) - 1.0) > 0.3 + +@pytest.mark.level_1 +@pytest.mark.tier_a +@pytest.mark.parametrize("smoother", ("sqrt", "powermean")) +def test_yield_anchor_never_falls_below_exact_min(smoother): + """A genuine approach from ABOVE. The onset anchor undershoots.""" + c = _model(smoother); c.yield_anchor = "yield" + for d in DELTAS[smoother]: + for f in FS: + assert eta(c, f, d) >= 1.0/max(1.0, f) - 1e-12, (d, f) + c.yield_anchor = "onset" + assert eta(c, 0.5, DELTAS[smoother][-1]) < 1.0 - 1e-3 + +@pytest.mark.level_1 +@pytest.mark.tier_a +def test_yield_anchor_removes_the_power_mean_degenerate_basin(): + """The p-norm's viscosity collapses towards zero as delta grows -- the 'spurious + fully-yielded fixed point'. It is an artefact of the missing 1/2: the power mean is + bounded below by exact Min for every delta, and tends to the geometric mean.""" + c = _model("powermean") + c.yield_anchor = "onset" + collapse = [eta(c, 10.0, d) for d in (1.0, 4.0, 16.0)] + assert collapse[-1] < 0.01 * 0.1 # 0.1 is exact Min at f=10 + assert collapse[0] > collapse[1] > collapse[2] # monotone collapse + c.yield_anchor = "yield" + bounded = [eta(c, 10.0, d) for d in (1.0, 4.0, 16.0)] + assert all(v >= 0.1 for v in bounded) + assert max(bounded) < np.sqrt(1.0 * 0.1) * 1.001 # geometric mean is the ceiling + +@pytest.mark.level_1 +@pytest.mark.tier_a +def test_anchor_rejects_nonsense_and_defaults_to_onset(): + c = _model() + assert c.yield_anchor == "onset" + with pytest.raises(ValueError): + c.yield_anchor = "middle"