yield_stress_min is not scaled like yield_stress: a dimensional floor changes the answer, a non-dimensional one does not

[lmoresi]

What we see

ViscoPlasticFlowModel.Parameters.yield_stress_min is not reconciled with the same
scaling as yield_stress. A floor set well below the yield stress, which must be inert,
changes the effective viscosity when given dimensionally, and is correctly inert when the
same physical value is given as a plain non-dimensional float.

Pure-shear box, reference quantities set (length 10.4 km, velocity 44.496 mm/yr,
viscosity 5.8333e19 Pa s), so the stress scale is 7.9086e6 Pa. Floor = 0.35 MPa, against a
median yield stress of 29.1 MPa — 83x below it.

yield_stress_min	median tau_y	min eta
unset (library default)	2.9143e+07	0.6673
plain float 0.0443 (non-dimensional)	2.9143e+07	0.6673
uw.quantity(0.350, "MPa") — same physical value	2.9154e+07	0.7455

A third spelling, 0.02 * C0 where C0 is the dimensional cohesion expression, gives a
third answer: yielding stops entirely, min eta exactly 1.0000.

So three ways of writing one physical stress give three different results, and only the
non-dimensional one behaves as a floor that far below the yield stress should.

Why it matters

The failure is silent and it looks like physics. In our case the model simply stopped
yielding, which reads as a modelling choice rather than a units fault. Anyone setting a
residual strength or a tension cutoff in physical units is affected, and the natural
spelling is the dimensional one.

Where it comes from

constitutive_models.py, in viscosity:

if inner_self.yield_stress_min.sym != -sympy.oo:
    yield_stress = self._apply_floor(inner_self.yield_stress, inner_self.yield_stress_min)

_apply_floor is smooth_max(value, floor, rounding), pure arithmetic on the two
operands. The two operands do not appear to arrive in the same system.

Reproduction

Full script: ~/+Simulations/gradient_damage/repro/yield_stress_min_units.py. The shape
of it:

model.set_reference_quantities(length=..., velocity=..., viscosity=...)
stress_scale = eta_bg * V_ms / H_m

cm.Parameters.yield_stress = smooth_max(C0 * seed + sin(phi) * p, 0, 0)
cm.yield_mode = "min"
cm.Parameters.yield_stress_min = <one of the three spellings above>
# viscous seed against an unreachable yield stress, then the plastic solve

Sample cm.viscosity at the velocity nodes and compare across the three.

What we would expect

Either the dimensional forms are scaled like yield_stress, or yield_stress_min refuses
a dimensional value rather than accepting it and mis-scaling it.

Underworld development team with AI support from Claude Code