Tension cutoff at yield_stress_min = 0 gives exactly zero viscosity and a singular operator (DIVERGED_LINEAR_SOLVE)

[lmoresi]

What we see

With gravity off, or anywhere the dynamic pressure goes negative, a Drucker-Prager yield
stress tau_y = C + sin(phi) p goes negative in tension. yield_stress_min defaults to
0, so tau_y is clamped to exactly zero, and then

eta_pl = tau_y / (2 edot_II) = 0

A zero-viscosity cell is a singular operator, not a weak one, and the tangent through it
is undefined. The linear solve then fails with DIVERGED_LINEAR_SOLVE.

The code already describes this. From _apply_floor in constitutive_models.py:

TODO(DESIGN): the rounding scale is RELATIVE (delta * floor), so it collapses for a
tension cutoff at floor = 0 — now the default for yield_stress_min. That leaves a hard
corner: the floored yield stress is exactly 0 in tension, hence
eta_pl = tau_y/(2 edot_II) is exactly 0 and the tangent through the soft-min is
undefined there. A properly rounded cap (Griffith / parabolic) needs an ABSOLUTE stress
scale, which this signature cannot supply. Maintainer decision pending (2026-07-26).

This issue is to record that the corner is reached in practice, and what it costs.

Measured

Pure-shear box, no gravity, cohesion softening C0 (1-D)^n. Once the band organises the
dynamic pressure changes sign (min p goes +2.13 to -0.38 over one damage step), and the
cohesion has collapsed with damage — C0 (1-D)^n is 1e-4 C0 at n = 4, D = 0.9, and 1e-8
at n = 8. tau_y reaches the zero floor, and every run in the affected regime ends in
DIVERGED_LINEAR_SOLVE rather than a nonlinear failure:

softening exponent	outcome
n = 4	runs to damage saturation
n = 8, four separate configurations	all DIVERGED_LINEAR_SOLVE, survival ordered by how far each delays the corner

Harder softening reaches the corner sooner, which is the signature: it is the cohesion
collapsing, not the mesh or the regime.

Two things a user cannot currently do

1. Set a residual strength in physical units — that is yield_stress_min is not scaled like yield_stress: a dimensional floor changes the answer, a non-dimensional one does not #758.
2. Get a rounded tension cap at all. The rounding scale is relative to the floor, so at
   floor = 0 there is no scale and the corner stays hard even in yield_mode="softmin".

Suggestion

A non-zero default residual strength would fix both the singular operator and the missing
rounding scale, since the rounding is relative to the floor. That is also what most codes
do — rock retains residual strength, and a tension cutoff at exactly zero cohesion is not
physical. The value has to come from somewhere explicit rather than being invented, so it
may belong as a required parameter rather than a default.

Related: #758 (the same parameter, its scaling).

Underworld development team with AI support from Claude Code