Yield edge cases: units in sums (#758), a stated residual strength (#759), no silent smoother bump (#782)

docs/advanced/benchmarks/_ti_vep_bdf2_isolation.py

@@ -66,6 +66,7 @@ def build(label, *, spatial_yield, yield_mode):
     cm.Parameters.shear_viscosity_1 = ETA_1
     cm.Parameters.shear_modulus = MU
     cm.Parameters.yield_stress = ty
+    cm.Parameters.yield_stress_min = -sympy.oo  # the fault ty = 1/weakness reads a field but stays >= TAU_Y > 0
     cm.Parameters.director = sympy.Matrix([0.0, 1.0])  # θ=0 throughout
     cm.Parameters.strainrate_inv_II_min = 1.0e-6
     cm._yield_mode = yield_mode

docs/advanced/benchmarks/bench_ti_vep_harmonic.py

@@ -130,6 +130,7 @@ def build_ti_stokes(label, theta_deg, tau_y, bdf_order):
     cm.Parameters.shear_viscosity_1 = ETA_1
     cm.Parameters.shear_modulus = MU
     cm.Parameters.yield_stress = tau_y_field
+    cm.Parameters.yield_stress_min = -sympy.oo  # reads the distance field, but 1/weakness >= tau_y > 0 by construction
     cm.Parameters.director = director
     cm.Parameters.shear_viscosity_min = ETA_0 * 1.0e-3
     cm.Parameters.strainrate_inv_II_min = 1.0e-6

docs/advanced/fault-networks.md

@@ -91,6 +91,7 @@ yield. Read the next paragraph before using it as glue.
 ```python
 stokes.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel
 stokes.constitutive_model.Parameters.yield_stress = net.band_yield(4.0)
+stokes.constitutive_model.Parameters.yield_stress_min = -sympy.oo  # a field, but tau_y or tau_far: positive by construction
 stokes.consistent_jacobian = True
 ```
 

docs/advanced/listric-fault-extension.md

@@ -69,6 +69,7 @@ net.prepare(h=CELL, ligament=1.0)
 net.mesh = base_mesh().add_fault(net.prepared)   # same base as the ribbons
 ...
 stokes.constitutive_model.Parameters.yield_stress = net.damage_yield(v, dial=0.05)
+stokes.constitutive_model.Parameters.yield_stress_min = -sympy.oo  # dial*(1 + 2 edot_II) > 0 by construction
 net.apply_contact(stokes, conds=ETA_F)
 ```
 

docs/developer/guides/nonlinear-solver.md

@@ -84,6 +84,7 @@ case felt like whack-a-mole. Check them first.
 | **A pressure-dependent yield stress makes the Newton Schur complement INDEFINITE far from the solution**: the consistent tangent differentiates η by p, so J_vp is no longer −Bᵀ (Spiegelman notch at η_bg 5e24: 41% off, 40 of 1968 Schur eigenvalues change sign). In the basin of convergence the full tangent is benign again (no wrong-sign eigenvalues, no negative pivots). | fieldsplit Schur diverges at ~0.99 per iteration even with the velocity block solved exactly; `DIVERGED_LINEAR_SOLVE` on the first Newton steps of a cold Drucker-Prager solve | `newton_pressure_coupling=False` (the Newton velocity block with the frozen pressure column, a partial tangent with the standard saddle-point structure) on the approach, or `consistent_jacobian="continuation"`; the default keeps the full tangent, which is right once inside the basin |
 | **Perfect plasticity's consistent tangent is SINGULAR along the flow**: on the hard-`Min` plastic branch η = τ_y/2ε̇_II, so 2η + 2η′ε̇_II = 0 — the velocity block is symmetric but semi-definite in every yielded cell. (An earlier version of this row blamed *asymmetry*; that is wrong for any η(ε̇_II) law — the rank-one term η′ ε̇⊗ε̇/ε̇_II is symmetric. Pressure-dependent yield adds a non-symmetric v–p coupling, not a non-symmetric velocity block. Corrected 2026-08-26, maintainer review.) | benign while yielded cells are few (the viscous neighbours regularise); with a large yielded fraction the velocity sub-solve caps out and Newton stalls at ~1e-3, no failure reason | give the plastic branch a positive tangent: the declared `Parameters.plastic_rate_strengthening` (η_reg adds 2η_reg along the flow in every yielded cell; its scale is the model's), a rounded viscosity floor (clips only the fastest cells). The default soft-min (sqrt, `"onset"`, δ = 0.1) does NOT do it — its yielded tangent is ≥ −1e-4 of the viscous one, still effectively singular; the sqrt law under `yield_anchor="yield"` is negative (−0.029). Picard converges regardless (full 2η stiffness) but is linear-rate. The FMG bundle's `gmres`+`sor` smoother is Newton-safe either way |
 | `preconditioner="fmg"` (vs explicit `pc_type=mg` + manual mg opts) | outer KSP "converges" in **1 iteration** → no real Newton correction → stall → `DIVERGED_LINE_SEARCH` | use explicit `pc_type=mg` with the smoother opts above; bound the outer KSP (`ksp_max_it`~80) so a hostile step fails fast |
+| A yield stress that reads a field (pressure, damage) left at the zero `yield_stress_min` | `ValueError: ... the yield stress reads ... residual strength` at setup (the cut-off would make η_pl = 0, a singular operator, #759) | state the residual strength (`yield_stress_min > 0`, stress units), or `-sympy.oo` if the yield stress is positive by construction |
 | Cold plastic start `v=0`, or any rigid/unyielded point | `DIVERGED_FNORM_NAN` at iteration 0 | **Not** a div/0: `ε̇=0` gives `η_pl=+inf`, which `Min` and the sqrt soft-min carry correctly to the viscous branch. Only a soft-min form that computes `η_ve·η_pl/(η_ve+η_pl)` breaks (`inf/inf`). Fixed in the power-mean; if you hand-roll a blend, write the harmonic mean as `η_ve/(1+η_ve/η_pl)`. **Do not reach for a strain-rate floor** — it hides this rather than fixing it |
 | LU velocity block with all-Dirichlet-ish BC | pressure nullspace singular | attach the Stokes nullspace / avoid a bare LU there |
 | Hand-rolled `snes_monitor` to "see what's happening" | you read residuals but miss the tell | use `solve_with_diagnostics` / `get_snes_diagnostics` instead (below) |

docs/developer/guides/plasticity-solvers.md

@@ -17,6 +17,7 @@ and records what was retired.
 ```python
 stokes.constitutive_model = cm           # ViscoPlastic / ViscoElasticPlastic / TI-VEP
 cm.Parameters.yield_stress = tau_y       # finite -> plasticity active
+cm.Parameters.yield_stress_min = tau_r   # residual strength: REQUIRED if tau_y reads a field (-sympy.oo if positive by construction)
 stokes.consistent_jacobian = True        # Newton tangent (non-elastic DP; see table)
 stokes.solve()                           # cold: ONE Picard step automatically, then Newton
 stokes.solve(picard=3)                   # the recipe: a few Picard steps, then Newton
@@ -235,6 +236,11 @@ could.
   plastic branch wherever `m * eta_reg >= shear_viscosity_min`; the floor remains a
   safeguard for the creep branch. `yield_stress_min` (the POSITIVE residual strength in
   tension) is a statement of the yield law and stays.
+- **`yield_stress_min` must be stated when the yield stress reads a field** (a
+  Drucker-Prager `C + sin(phi) p`, a cohesion weakening with damage). At the zero
+  default a yield stress cut off in tension gives `eta_pl = 0`, a singular operator
+  (`DIVERGED_LINEAR_SOLVE`, #759); the model refuses it. State the residual strength
+  in stress units, or `-sympy.oo` for a yield stress positive by construction.
 
 - **`shear_viscosity_min`** (default `-oo` = off) is applied through
   `uw.maths.smooth_max`, but the default rounding scale is zero under
@@ -258,6 +264,7 @@ could.
 |---------|-------|-----|
 | `DIVERGED_LINEAR_SOLVE`, 0 iters, VEP | consistent Newton over the elastic block → indefinite | Picard (`consistent_jacobian=False`) |
 | `DIVERGED_LINEAR_SOLVE` at nl=0 with a viscosity floor set | δ→0 leaves the floor's `Max` corner exact | set `viscosity_min_rounding` |
+| `ValueError: the yield stress depends on a field ... residual strength` | a field-dependent yield stress with the zero default floor (#759) | set `yield_stress_min` to the residual strength, or `-sympy.oo` if the yield stress cannot reach zero |
 | Newton stalls with no divergence reason | admissible uselessness — steps accepted, residual flat | revert to best iterate, Picard block (`solve(picard=N)` or `"continuation"`); consider grid sequencing |
 | Converges but σ sits **below** τ_y | a fixed δ>0 soft-min under the default `"onset"` anchor is a WEAKER law (4.5 % at the yield point at δ = 0.1) | that is the modelling choice you made — a smaller δ, or `yield_mode="min"` for the exact surface (`yield_anchor="yield"` reaches τ_y but stiffens every unyielded cell and softens above yield) |
 | Linear (~20-iter) convergence | Picard tangent when you wanted Newton | `consistent_jacobian=True` on a non-elastic model (see "Confirm" above) |

docs/examples/WIP/VEP_Fault_Investigation.py

@@ -87,6 +87,7 @@
 cm.Parameters.shear_viscosity_0 = ETA
 cm.Parameters.shear_modulus = MU
 cm.Parameters.yield_stress = tau_y_field
+cm.Parameters.yield_stress_min = -sympy.oo  # reads the distance field, but 1/weakness >= TAU_Y_FAULT > 0 by construction
 # yield_mode="softmin" is the default (smooth approximation to Min)
 cm.Parameters.shear_viscosity_min = ETA * 1.0e-2
 cm.Parameters.strainrate_inv_II_min = 1.0e-5

docs/examples/solid_mechanics/advanced/Ex_Shear_Band_Plasticity_SS.py

@@ -182,7 +182,7 @@ class boundaries(Enum):
 stokes.constitutive_model = uw.constitutive_models.ViscoPlasticFlowModel
 stokes.constitutive_model.Parameters.shear_viscosity_0 = 1
 stokes.constitutive_model.Parameters.yield_stress = mu * p_soln.sym[0] + C
-stokes.constitutive_model.Parameters.yield_min = 1000
+stokes.constitutive_model.Parameters.yield_stress_min = 1000  # residual strength (= C): the Drucker-Prager yield reads p
 stokes.penalty = 1
 
 stokes.petsc_options["ksp_monitor"] = None

docs/examples/solid_mechanics/advanced/Ex_Sheared_Layer_Test.py

@@ -427,6 +427,7 @@ class boundaries(Enum):
 C = 2 * (y * 2) ** 2 + (C0 - Cinf) * (1 - sympy.tanh((strain.sym[0] / eps_ref - 1) * scale)) / 2 + Cinf
 
 stokes.constitutive_model.Parameters.yield_stress = C + mu * p_soln.sym[0]
+stokes.constitutive_model.Parameters.yield_stress_min = 0.1 * Cinf  # residual strength: 10% of the softened cohesion (yield reads p and strain)
 stokes.constitutive_model.Parameters.edot_II_fn = stokes.Unknowns.Einv2
 stokes.constitutive_model.Parameters.min_viscosity = 0.1
 stokes.saddle_preconditioner = 1 / stokes.constitutive_model.viscosity

docs/examples/solid_mechanics/intermediate/Ex_TI_VEP_Angled_Fault.py

@@ -165,6 +165,7 @@
 cm.Parameters.shear_viscosity_1 = ETA_1
 cm.Parameters.shear_modulus = MU
 cm.Parameters.yield_stress = tau_y_field
+cm.Parameters.yield_stress_min = -sympy.oo  # reads the distance field, but 1/weakness >= TAU_Y > 0 by construction
 cm.Parameters.director = director
 cm.Parameters.shear_viscosity_min = ETA_0 * 1.0e-3
 cm.Parameters.strainrate_inv_II_min = 1.0e-6

src/underworld3/constitutive_models.py

@@ -75,6 +75,72 @@ def _parameter_is_zero(par):
         return False
 
 
+def _fields_read_by(par):
+    """The names of the mesh or swarm variables a parameter reads (the pressure, a
+    damage field, ...): a value not known until the solve, which can change sign."""
+    from underworld3.function.expressions import unwrap_expression
+    try:
+        value = unwrap_expression(par, mode="nondimensional")
+    except (TypeError, ValueError, AttributeError):
+        # Sanctioned: a parameter that cannot be unwrapped is not inspected, and the
+        # residual-strength check then trusts it (fails open, as before #759)
+        return []
+    names = set()
+    for f in sympy.sympify(value).atoms(sympy.Function):
+        if hasattr(f, "meshvar"):
+            variable = f.meshvar() if callable(f.meshvar) else f.meshvar
+            names.add(getattr(variable, "name", None) or str(f.func))
+    return sorted(names)
+
+
+def _numeric_value(par):
+    """A parameter's non-dimensional value as a float, or None when it is not a
+    number (an expression of a field)."""
+    from underworld3.function.expressions import unwrap_expression
+    try:
+        value = unwrap_expression(par, mode="nondimensional")
+        return float(getattr(value, "magnitude", value))
+    except (TypeError, ValueError, AttributeError):
+        # Sanctioned: an expression has no single value; None says so
+        return None
+
+
+def _floor_is_off(floor):
+    """True when the yield-stress floor is disabled: -oo, bare or wrapped in an
+    expression (a wrapped -oo applied as a floor evaluates to NaN)."""
+    if floor.sym == -sympy.oo:
+        return True
+    value = _numeric_value(floor)
+    return value is not None and value == -float("inf")
+
+
+def _require_residual_strength(constitutive_model, parameters):
+    """Refuse a yield stress that reads a field (Drucker-Prager through the pressure,
+    a cohesion that weakens with damage) with a floor that is not positive: where it
+    reaches the floor the plastic viscosity tau_y / (2 edot_II) is zero or negative
+    and the operator singular (#759). The residual strength must be stated, > 0 in
+    stress units; an explicit -oo removes the floor for a yield stress known to stay
+    positive. A floor that is itself a field expression is trusted as stated."""
+    floor = parameters.yield_stress_min
+    if _floor_is_off(floor):
+        return
+    value = _numeric_value(floor)
+    if value is not None and value > 0.0:
+        return
+    if value is None and not _parameter_is_zero(floor):
+        return                                     # a floor that is itself a field
+    fields = _fields_read_by(parameters.yield_stress)
+    if fields:
+        material = getattr(constitutive_model, "_material_name", None)
+        where = type(constitutive_model).__name__ + (f" ({material})" if material else "")
+        raise ValueError(
+            f"{where}: the yield stress reads {', '.join(fields)} and can fall to the "
+            f"floor yield_stress_min = {floor.sym}, where the plastic viscosity is zero or "
+            "negative and the operator singular. Set Parameters.yield_stress_min to the "
+            "residual strength (> 0, in stress units), or to -sympy.oo if the yield "
+            "stress is known to stay positive.")
+
+
 def _plastic_rate_strengthening_parameter():
     """The declared plastic rate strengthening, shared by every yielding model's
     Parameters (ViscoPlastic, ViscoElasticPlastic, TransverseIsotropicVEP)."""
@@ -1379,9 +1445,11 @@ def yield_smoother(self):
         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`` is the sharp
-        hard-``Min`` limit it only *approaches*.
+        Selecting a family never changes ``yield_softness``: the power mean runs
+        at the δ the model states, ``s = 1/(δ + 0.001)`` -- ``δ = 0`` is
+        ``s = 1000``, close to the exact ``Min``; ``δ = 1`` is the harmonic mean,
+        a distinct, strongly smoothed law (#782: selecting it used to bump a zero
+        δ to 1 silently).
         """
         return getattr(self, "_yield_smoother", "sqrt")
 
@@ -1392,10 +1460,6 @@ def yield_smoother(self, value):
                 f"yield_smoother must be 'sqrt' or 'powermean', got '{value}'"
             )
         self._yield_smoother = value
-        if value == "powermean" and getattr(self, "_yield_softness", 0.0) == 0.0:
-            # δ=0 ⇒ s=1/δ=∞ is the singular Min limit; default to the
-            # parameter-free harmonic mean (s=1) instead.
-            self.yield_softness = 1.0
         self._reset()
 
     # --- Smooth lower bounds (viscosity / yield floors) -----------------------------
@@ -1706,7 +1770,7 @@ class _Parameters(_ParameterBase, _ViscousParameterAlias):
         yield_stress_min = api_tools.Parameter(
             R"{\tau_{y, \mathrm{min}}}",
             lambda inner_self: 0,
-            "Yield stress (DP) minimum cutoff",
+            "Residual strength: the floor of the yield stress (stated, > 0, when the yield stress reads a field)",
             units="Pa",
         )
 
@@ -1751,8 +1815,10 @@ def viscosity(self):
         # than propagated into tau_y/(2 edot_II). (The ±oo defaults elsewhere in
         # Parameters exist so sympy can cancel an unused term away; that trick is
         # wrong here, so this one defaults to 0 — maintainer ruling 2026-07-26.)
-        # An explicit -oo still disables the floor.
-        if inner_self.yield_stress_min.sym != -sympy.oo:
+        # An explicit -oo still disables the floor. A yield stress that reads a field
+        # needs a stated, positive floor (_require_residual_strength).
+        _require_residual_strength(self, inner_self)
+        if not _floor_is_off(inner_self.yield_stress_min):
             yield_stress = self._apply_floor(
                 inner_self.yield_stress, inner_self.yield_stress_min
             )
@@ -1812,7 +1878,7 @@ def _rate_strengthening_probe_terms(self):
         # the viscosity above, so the ladder reads the law the solve runs.
         inner_self = self.Parameters
         tau = inner_self.yield_stress
-        if inner_self.yield_stress_min.sym != -sympy.oo:
+        if not _floor_is_off(inner_self.yield_stress_min):
             tau = self._apply_floor(tau, inner_self.yield_stress_min)
         edot = self._strainrate_inv_II
         if not _parameter_is_zero(inner_self.strainrate_inv_II_min):
@@ -2189,7 +2255,7 @@ def dt_elastic(inner_self, value):
         yield_stress_min = api_tools.Parameter(
             R"{\tau_{y, \mathrm{min}}}",
             lambda inner_self: 0,
-            "Yield stress (DP) minimum cutoff",
+            "Residual strength: the floor of the yield stress (stated, > 0, when the yield stress reads a field)",
             units="Pa",
         )
 
@@ -2835,8 +2901,10 @@ def _plastic_effective_viscosity(self):
         # Guard on the DISABLING sentinel, not on zero: zero is the default and a
         # physically meaningful floor (the yield stress is compared against the second
         # invariant of the stress, so a negative tau_y is meaningless). Only an
-        # explicit -oo turns the floor off.
-        if parameters.yield_stress_min.sym != -sympy.oo:
+        # explicit -oo turns the floor off. A yield stress that reads a field needs a
+        # stated, positive floor (_require_residual_strength).
+        _require_residual_strength(self, parameters)
+        if not _floor_is_off(parameters.yield_stress_min):
             # Literal 0 for the default floor, not the parameter atom: sympy cannot
             # fuzzy-compare an opaque UWexpression and Max canonicalisation recurses.
             # smooth_max keeps both operands symbolic (the JIT routes the floor
@@ -3910,7 +3978,7 @@ def dt_elastic(inner_self, value):
         yield_stress_min = api_tools.Parameter(
             R"{\tau_{y, \mathrm{min}}}",
             lambda inner_self: 0,
-            "Yield stress minimum cutoff", units="Pa",
+            "Residual strength: the floor of the yield stress (stated, > 0, when the yield stress reads a field)", units="Pa",
         )
         strainrate_inv_II_min = api_tools.Parameter(
             R"{\dot\varepsilon_{II,\mathrm{min}}}",
@@ -4270,7 +4338,9 @@ def _plastic_effective_viscosity(self):
         tau_y = parameters.yield_stress
         # Guard on the DISABLING sentinel, not on zero — zero is the default and a
         # real floor (a negative yield stress is meaningless against an invariant).
-        if parameters.yield_stress_min.sym != -sympy.oo:
+        # A yield stress that reads a field needs a stated, positive floor.
+        _require_residual_strength(self, parameters)
+        if not _floor_is_off(parameters.yield_stress_min):
             # Literal 0 for the default floor (sympy fuzzy-compare recursion on atoms).
             # smooth_max: keeps the floor symbolic, no ordering test (see the note
             # in Constitutive_Model._apply_floor). eps = 0 is exactly Max.
@@ -4794,6 +4864,8 @@ def _eta_par_eff_lagged(self):
             except (TypeError, ValueError):
                 edot_min_val = 1.0e-6
 
+        # TODO(BUG): the lagged estimate uses the unfloored yield stress, so it ignores
+        # yield_stress_min and the residual-strength rule (#759); see _require_residual_strength.
         tau_y_sym = params.yield_stress.sym
         vp_eff_lag = tau_y_sym / (
             2 * (gamma_dot_abs_lag + sympy.Float(edot_min_val))