shear_viscosity_min: the floor's rounding scale is tied to yield_mode (and the VEP path silently skips it)

[lmoresi]

What

Constitutive_Model._apply_floor (src/underworld3/constitutive_models.py:1160) picks the
rounding scale for a viscosity/yield-stress floor from yield_mode:

rounding = 0 if getattr(self, "_yield_mode", "min") == "min" \
    else self._get_yield_softness() * floor
return uw.maths.smooth_max(value, floor, rounding)

Two consequences, and they are the two halves of the same design gap — the floor has no
smoothing parameter of its own.

1. Under yield_mode="min" the floor is a hard sympy.Max. That is correct for the
frozen/Picard tangent, but consistent_jacobian=True differentiates through the kink and
the solve dies immediately. Measured on the Spiegelman notch, isolated with a
floor-on/off × Newton/frozen 2×2 — the only failing cell is (floor on, Newton), at
nl=0, DIVERGED_LINEAR_SOLVE. So a user who sets a physically sensible viscosity floor
silently loses the consistent tangent.

2. Under yield_mode="softmin" the rounding is δ·floor, so it vanishes as δ→0.
The floor is smooth at the start of a yield homotopy and hard at the end — i.e. the
smoothing disappears exactly where the march is trying to land. δ is the yield corner's
regularisation parameter; the floor is a different corner and should not inherit it.

Why it matters now

A viscosity floor that is relaxed toward zero is one of the classical continuation
routes for viscoplastic Stokes (it bounds the plastic viscosity from below and is then
removed — structurally the same idea as the Perzyna-style rate term ξ). That route cannot
be measured against the soft-min homotopies at all while (1) forces it onto the Picard
tangent and (2) couples its smoothing to the parameter being marched.

Suggested shape of a fix

Give the floor its own rounding scale rather than deriving one from yield_mode:

• an explicit smoothing parameter on the floor (absolute, or relative to the floor with
  its own coefficient), defaulting to 0 so today's min behaviour is unchanged;
• _apply_floor uses that parameter regardless of yield_mode.

There is a related TODO(DESIGN) already in the same function (lines 1180-1186) about the
relative rounding collapsing for a tension cutoff at floor = 0 — the two want settling
together, since both are "this corner has no length scale of its own".

Underworld development team with AI support from Claude Code

[lmoresi]

Correction to the description above

I mis-stated the mechanism when filing this. The title and body say _apply_floor "emits a
hard sympy.Max" — that is wrong for the ViscoPlasticFlowModel path, and it matters,
because it points at the wrong fix.

_apply_floor already routes through uw.maths.smooth_max:

def smooth_max(a, b, epsilon):
    return (a + b + sympy.sqrt((a - b) ** 2 + epsilon**2)) / 2

That is pure arithmetic. sympy is never asked for an ordering test, so the
UWexpression-recursion problem I invoked does not arise — and notch_model.py has been
using exactly this construction for the tension cutoff on σ_y all along.

The actual defect is only the choice of ε, which is what the rest of the issue
describes correctly:

• rounding = 0 under yield_mode="min" makes the expression algebraically exactly
  Max — a genuine kink, which is what breaks the consistent tangent. Not a sympy.Max
  object, but the same non-differentiability.
• rounding = delta * floor under the smooth modes vanishes as δ → 0, i.e. exactly where a
  yield homotopy lands.

One place the original wording was right, and worse than I said:
ViscoElasticPlasticFlowModel (constitutive_models.py ~2148) uses a literal
sympy.Max(effective_viscosity, shear_viscosity_min) and silently returns the unfloored
viscosity when is_viscoplastic and yield_mode in ("harmonic", "softmin") — so on that
path a user-set floor is quietly ignored in the smooth yield modes.

What has been done

On feature/regime-diagram: a viscosity_min_rounding property gives the floor a rounding
scale of its own, independent of yield_mode and of δ. None (default) keeps the historical
behaviour, so nothing changes for existing models. With an explicit rounding there is no
outer Max to nest, so the VEP path can honour the floor in the smooth modes too and the
silent skip is bypassed.

Still open, and why this issue should stay open

I have not reproduced the original failure, so the fix is unverified. The first attempt
used a floor of 0.02 (non-dimensional) which turned out not to bind: measured, the
domain minimum of η_eff is 1.000e-02 at every ηV tested and that is the weak substrate,
while yielded material in the strong layer only reaches 6.5e-02. eps=0 and eps>0 gave
identical iteration counts and the eps=0 arm that was supposed to fail converged in 7.

With a material-scoped floor that does bind, there is now a reproducer:

floor | strong min | weak min  | nl | |F|/|F0|   reason
0     | 6.5229e-02 | 1.000e-02 | 30 | 5.891e-01  DIVERGED_MAX_IT
0.3   | 3.0090e-01 | 1.000e-02 |  0 | 1.000e+00  DIVERGED_LINEAR_SOLVE   <-- this
0.6   | 6.0117e-01 | 1.000e-02 | 30 | 3.979e-02  DIVERGED_MAX_IT

floor = 0.3 fails at nl=0 with a stronger floor (0.6) working either side of it, under
consistent Newton with the rounding set to 5% of the floor. A non-monotonic failure in the
middle of the range is not "the floor is too strong", so either the rounding is not doing
its job or the Piecewise floor is doing something to the operator at the material
boundary. That is the next thing to diagnose here.

Underworld development team with AI support from Claude Code

[lmoresi]

Addressed by PR #475, merged to development 2026-08-02 (manual close — auto-close doesn't fire on development-targeted PRs). viscosity_min_rounding decouples the floor's rounding scale from yield_mode, and the VEP path honours the floor whenever it is set (constitutive_models.py ~2187-2196). Default is None = historical behavior, preserved deliberately — the defect as described no longer applies once the knob is set. The unvalidated setter is tracked in #490. Reopen if the DEFAULT behavior should change rather than being opt-in.