From ffdbfd5fcee92953d49179912baa9002ec707821 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 22 Jul 2026 11:50:00 +1000 Subject: [PATCH] feat(constitutive): power-mean yield soft-min smoother + rampable delta constant Add the power-mean soft-min yield smoother and make the soft-min softness delta a constants[] atom (runtime-rampable with no JIT recompile) on the ViscousFlowModel base class, so the visco-plastic subclasses inherit one implementation. This is the generalisable, scalable substrate for a yield homotopy toward the sharp Min surface. Base class (ViscousFlowModel): - _combine_yield(eta_ve, eta_pl): the shared viscous/plastic combination, keyed on yield_mode -- "min" (exact hard Min), "harmonic", or "softmin" (the delta soft-min in the family chosen by yield_smoother). Development's tri-modal yield_mode semantics are preserved exactly (unlike the source branch, which collapsed "softmin" into "min"). - _get_yield_softness / _get_yield_offset: delta and the onset offset held as constants[] UWexpression atoms (the offset defined symbolically in terms of delta -- one symbol in the stress tensor, so no tensor blow-up -- so one delta update repacks both). Ramping delta via the atom + solver._update_constants() forces no recompile. - yield_smoother property/setter: "sqrt" (default; overshoots tau_y in the transition) or "powermean" (undershoots tau_y -- eta_eff <= Min always -- overflow-safe on geodynamic viscosity ranges). Selecting "powermean" from delta=0 bumps delta to 1 (s=1 harmonic mean), since delta=0 is the singular Min limit. Models (Stage 1 -- isotropic, test-covered): - ViscoPlasticFlowModel (Drucker-Prager): defaults to yield_mode="min" == today's exact hard Min (zero behaviour change), and GAINS yield_mode / yield_softness knobs so it can opt into the softmin / power-mean homotopy -- the model the hard-case DP homotopy needs. - ViscoElasticPlasticFlowModel (VEP): its inline sqrt soft-min routed through _combine_yield; delta moves float -> atom (value-identical); yield_softness setter syncs the atom. Zero behaviour change at stock settings: _combine_yield at default settings is bit-identical to the previous inline law (delta merely moves from a baked float to a constants[] symbol evaluated to the same number). Power-mean and runtime-delta are opt-in. The distrusted in-solve enable_yield_homotopy ramp driver is intentionally NOT part of this PR. TI-VEP / TI-VEP-Split (4 anisotropic sites) are a Stage-2 follow-up. Tests: tests/test_1055_yield_smoother.py -- delta=0 exact Min; power-mean undershoot + overflow-safety + ->Min; smoother validation/default; runtime delta ramp no-recompile; DP opt-in. level_1/tier_a gate green (402 passed, 0 failed). VEP _combine_yield verified bit-identical to the prior float law at defaults. Underworld development team with AI support from Claude Code --- src/underworld3/constitutive_models.py | 216 +++++++++++++++++++++---- tests/test_1055_yield_smoother.py | 184 +++++++++++++++++++++ 2 files changed, 373 insertions(+), 27 deletions(-) create mode 100644 tests/test_1055_yield_smoother.py diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 2944ea0c..bb8a70ce 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -899,6 +899,125 @@ def _object_viewer(self): ) ) + # --- Yield soft-min smoother (shared by the visco-plastic subclasses) ----------- + # The δ soft-min regularisation and the smooth-min FAMILY selection live on the + # base class so every yielding model inherits one implementation. δ is held as a + # constants[] UWexpression atom (not a baked float) so a homotopy can ramp it at + # runtime via PetscDSSetConstants with no JIT recompile. + + def _get_yield_softness(self): + r"""The soft-min regularisation δ as a ``constants[]`` UWexpression atom. + + Created lazily and kept in sync with ``self._yield_softness`` (the + configured numeric value). Storing δ as a UWexpression — rather than + baking the float into the compiled flux — lets a yield homotopy ramp + δ at runtime via ``PetscDSSetConstants`` with no JIT recompile. + ``δ = 0`` makes the sqrt law identically ``Min``. + """ + delta_value = getattr(self, "_yield_softness", 0.0) + if getattr(self, "_yield_softness_expr", None) is None: + self._yield_softness_expr = expression( + R"{\updelta_{y}}", + 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). + 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", + ) + else: + self._yield_softness_expr.sym = sympy.Float(delta_value) + return self._yield_softness_expr + + def _get_yield_offset(self): + """The onset-offset constant atom (lazily created alongside δ).""" + if getattr(self, "_yield_offset_expr", None) is None: + self._get_yield_softness() + return self._yield_offset_expr + + 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``: + + - ``"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). + Both approach exact ``Min`` as ``δ → 0``. + + Behaviour at the stock settings (``yield_mode`` per model default, + ``yield_smoother="sqrt"``) is identical to the previous inline law — δ merely + moves from a baked float to a ``constants[]`` atom (same value). + """ + mode = getattr(self, "_yield_mode", "softmin") + if mode == "harmonic": + return 1 / (1 / eta_ve + 1 / eta_pl) + if mode == "min": + return sympy.Min(eta_ve, eta_pl) + + # "softmin": δ-parameterised smooth-min family. + smoother = getattr(self, "_yield_smoother", "sqrt") + 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 + # 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)) + a = 1 + f + b = 1 + 1 / f + N = eta_ve * eta_pl / (eta_ve + eta_pl) + return N * (a ** (-s) + b ** (-s)) ** (-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() + g = 1 + (f - 1 + sympy.sqrt((f - 1) ** 2 + delta**2)) / 2 - offset + return eta_ve / g + + @property + def yield_smoother(self): + r"""Which smooth-min FAMILY regularises the ``"softmin"`` yield mode. + + Both families use the same softness parameter ``δ`` (``yield_softness``) + and approach exact ``Min`` as ``δ → 0``, but differ in how they round + 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. + + 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*. + """ + return getattr(self, "_yield_smoother", "sqrt") + + @yield_smoother.setter + def yield_smoother(self, value): + if value not in ("sqrt", "powermean"): + raise ValueError( + 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() + ## NOTE - retrofit VEP into here @@ -963,6 +1082,16 @@ def __init__(self, unknowns, material_name: str = None): "Effective viscosity (plastic)", ) + # Yield-combination mode (see _combine_yield on the base class). Default + # "min" = the exact hard Min(η_0, η_yield) this model has always used, so the + # default behaviour is unchanged. Opt into "softmin" (+ yield_smoother / + # yield_softness) for the δ-parameterised smooth-min homotopy. + self._yield_mode = "min" + self._yield_softness = 0.0 # δ; 0 ⇒ exact Min + self._yield_smoother = "sqrt" # smooth-min family: "sqrt" | "powermean" + self._yield_softness_expr = None # constants[] δ atom (created lazily) + self._yield_offset_expr = None # onset-offset atom (created lazily) + class _Parameters(_ParameterBase, _ViscousParameterAlias): """Any material properties that are defined by a constitutive relationship are collected in the parameters which can then be defined/accessed by name in @@ -1044,15 +1173,13 @@ def viscosity(self): viscosity_yield = yield_stress / (2 * self._strainrate_inv_II) - ## Question is, will sympy reliably differentiate something - ## with so many Max / Min statements. The smooth version would - ## be a reasonable alternative: - - # effective_viscosity = sympy.sympify( - # 1 / (1 / inner_self.shear_viscosity_0 + 1 / viscosity_yield), - # ) - - effective_viscosity = sympy.Min(inner_self.shear_viscosity_0, viscosity_yield) + # Combine the viscous and plastic (yield) viscosities. The default + # yield_mode="min" gives the exact hard Min(η_0, η_yield); yield_mode="softmin" + # opts into the δ-parameterised smooth-min (sqrt or powermean family) for a + # scalable homotopy toward the sharp yield surface. + effective_viscosity = self._combine_yield( + inner_self.shear_viscosity_0, viscosity_yield + ) # If we want to apply limits to the viscosity but see caveat above # Keep this as an sub-expression for clarity @@ -1068,6 +1195,44 @@ def viscosity(self): # Returns an expression that has a different description return self._plastic_eff_viscosity + @property + def yield_mode(self): + r"""How the viscous and plastic (yield) viscosities are combined. + + - ``"min"`` (default): exact hard ``Min(η_0, η_yield)`` — the sharp yield + surface this model has always used. + - ``"harmonic"``: ``1/(1/η_0 + 1/η_yield)`` — a smooth blend. + - ``"softmin"``: the δ-parameterised smooth-min (family set by + ``yield_smoother``), for a scalable homotopy toward the sharp surface. + """ + return self._yield_mode + + @yield_mode.setter + def yield_mode(self, value): + if value not in ("min", "harmonic", "softmin"): + raise ValueError( + f"yield_mode must be 'min', 'harmonic', or 'softmin', got '{value}'" + ) + self._yield_mode = value + self._reset() + + @property + def yield_softness(self): + r"""Soft-min regularisation δ for ``yield_mode="softmin"`` (0 ⇒ exact Min). + + δ is held as a ``constants[]`` atom, so ramping it at runtime (this setter, + or ``cm._get_yield_softness().sym = ...`` + ``solver._update_constants()``) + does not trigger a JIT recompile — the basis of a scalable yield homotopy. + """ + return self._yield_softness + + @yield_softness.setter + def yield_softness(self, value): + self._yield_softness = float(value) + if getattr(self, "_yield_softness_expr", None) is not None: + self._yield_softness_expr.sym = sympy.Float(self._yield_softness) + self._reset() + def plastic_correction(self) -> float: r"""Scaling factor to reduce stress to yield surface. @@ -1217,6 +1382,9 @@ def __init__(self, unknowns, order=1, integrator: str = "bdf", self._order = order self._yield_mode = "softmin" # "min", "harmonic", "smooth", or "softmin" self._yield_softness = 0.1 # δ parameter for "softmin" mode + self._yield_smoother = "sqrt" # smooth-min family: "sqrt" | "powermean" + self._yield_softness_expr = None # constants[] δ atom (created lazily) + self._yield_offset_expr = None # onset-offset atom (created lazily) # Timestep — set by the solver before each solve(). Not a user parameter. # Initialised to oo (viscous limit). The solver overwrites this with the @@ -1676,23 +1844,13 @@ def viscosity(self): if self.is_viscoplastic: vp_effective_viscosity = self._plastic_effective_viscosity - if self._yield_mode == "harmonic": - effective_viscosity = 1 / (1 / effective_viscosity + 1 / vp_effective_viscosity) - elif self._yield_mode == "softmin": - # Smooth approximation to Min(η_ve, η_pl): - # η_eff = η_ve / g(f) - # g(f) = 1 + softplus(f-1) - softplus(-1) ≈ max(1, f) - # where softplus(x) = (x + √(x² + δ²))/2 and f = η_ve/η_pl. - # Corrected so g(0) = 1 exactly (no spurious yield below onset). - # Approaches exact Min as δ→0. No Min/Max in expression. - delta = self._yield_softness - f = effective_viscosity / vp_effective_viscosity - import math # float offset avoids sympy expression blowup in tensor - offset = (-1 + math.sqrt(1 + delta**2)) / 2 - g = 1 + (f - 1 + sympy.sqrt((f - 1)**2 + delta**2)) / 2 - offset - effective_viscosity = effective_viscosity / g - else: - effective_viscosity = sympy.Min(effective_viscosity, vp_effective_viscosity) + # Combine η_ve with the plastic viscosity per yield_mode (harmonic / exact + # Min / δ-soft-min). The soft-min softness δ now lives in a constants[] atom + # (see _combine_yield / _get_yield_softness) — value-identical to the former + # inline float law at the same δ, but runtime-rampable with no recompile. + effective_viscosity = self._combine_yield( + effective_viscosity, vp_effective_viscosity + ) # Apply viscosity floor — but skip for smooth-blend yield modes # where the outer Max creates a nested Min/Max that breaks the @@ -1984,7 +2142,11 @@ def yield_softness(self): @yield_softness.setter def yield_softness(self, value): - self._yield_softness = value + self._yield_softness = float(value) + # Keep the constants[] δ atom in sync (created lazily on first use) so a + # numeric δ assignment is reflected without a JIT recompile. + if getattr(self, "_yield_softness_expr", None) is not None: + self._yield_softness_expr.sym = sympy.Float(self._yield_softness) self._reset() @property diff --git a/tests/test_1055_yield_smoother.py b/tests/test_1055_yield_smoother.py new file mode 100644 index 00000000..9e0640f1 --- /dev/null +++ b/tests/test_1055_yield_smoother.py @@ -0,0 +1,184 @@ +"""δ-parameterised yield soft-min smoother (sqrt + power-mean families). + +Pins the smoother contract on development (the ``_combine_yield`` helper + the +runtime-rampable δ ``constants[]`` atom, ported from the yield-homotopy work, with +development's tri-modal ``yield_mode`` semantics preserved): + +1. ``test_delta0_is_exact_min`` — the sqrt soft-min at δ=0 equals ``Min(η_ve, η_pl)`` + to machine precision (g(f) = max(1, f)). +2. ``test_powermean_smoother_undershoots_min`` — the power-mean family is ≤ Min + everywhere (never over-yields), overflow-safe on geodynamic ranges, and → Min as δ→0. +3. ``test_yield_smoother_validation`` — only the two known families; default "sqrt". +4. ``test_yield_softness_runtime_no_recompile`` — ramping δ through its constants[] atom + does NOT change the solver's JIT cache key. +5. ``test_dp_model_smoother_optin`` — the non-elastic Drucker–Prager model defaults to + exact hard Min and can opt into the softmin/power-mean homotopy. + +NOTE (vs the branch): development keeps ``yield_mode`` tri-modal — ``"min"`` is the exact +hard ``Min`` (not the soft-min), so the smoother tests select ``yield_mode="softmin"``. +The ``enable_yield_homotopy`` in-solve ramp driver is intentionally NOT part of this PR. +""" + +import numpy as np +import pytest +import sympy + +import underworld3 as uw +from underworld3.function import expression +from underworld3.function.expressions import unwrap_expression + +ETA = 1.0 +MU = 1.0 +TAU_Y = 0.5 +V0 = 0.5 +T_R = ETA / MU + + +def _build_vep(label, order=2): + mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 8), minCoords=(-1.0, -0.5), maxCoords=(1.0, 0.5), + ) + v = uw.discretisation.MeshVariable(f"U_{label}", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable(f"P_{label}", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + cm = uw.constitutive_models.ViscoElasticPlasticFlowModel(stokes.Unknowns, order=order) + stokes.constitutive_model = cm + cm.Parameters.shear_viscosity_0 = ETA + cm.Parameters.shear_modulus = MU + cm.Parameters.yield_stress = TAU_Y + cm.Parameters.strainrate_inv_II_min = 1.0e-6 + V_top = expression(rf"V_{{{label}}}", sympy.Float(V0), "Top V") + stokes.add_dirichlet_bc((V_top, 0.0), "Top") + stokes.add_dirichlet_bc((-V_top, 0.0), "Bottom") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Left") + stokes.add_dirichlet_bc((sympy.oo, 0.0), "Right") + stokes.tolerance = 1.0e-6 + stokes.petsc_options["snes_force_iteration"] = True + return mesh, stokes, cm, V_top + + +@pytest.mark.level_1 +@pytest.mark.tier_a +def test_delta0_is_exact_min(): + """δ=0 sqrt soft-min law equals Min(η_ve, η_pl) to machine precision.""" + mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) + v = uw.discretisation.MeshVariable("Um", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("Pm", mesh, 1, degree=1) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + cm = uw.constitutive_models.ViscoElasticPlasticFlowModel(s.Unknowns, order=2) + s.constitutive_model = cm + cm.yield_mode = "softmin" # the δ-family path (dev "min" = hard Min) + cm.yield_softness = 0.0 # δ=0 ⇒ g(f)=max(1,f) ⇒ exact Min + for eta_ve, eta_pl in [(1.0, 2.0), (2.0, 1.0), (1.0, 1.0), (0.3, 5.0), (5.0, 0.3)]: + comb = cm._combine_yield(sympy.Float(eta_ve), sympy.Float(eta_pl)) + val = float(unwrap_expression(comb, mode="nondimensional")) + assert abs(val - min(eta_ve, eta_pl)) < 1.0e-12 + + +@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/δ→∞).""" + 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) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + cm = uw.constitutive_models.ViscoElasticPlasticFlowModel(s.Unknowns, order=2) + s.constitutive_model = cm + + cm.yield_mode = "softmin" # dev: soft-min family lives under "softmin" + 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" + assert cm.yield_softness == 1.0 + + # Includes a geodynamic-range pair (1e25, 1e21) that overflows the naive + # η^(−s) form above s≈40 but is finite in the harmonic-normalised form. + pairs = [(1.0, 2.0), (2.0, 1.0), (1.0, 1.0), (0.3, 5.0), (5.0, 0.3), (1e25, 1e21)] + for delta in (1.0, 0.5, 0.1, 0.02): + cm.yield_softness = delta + for eta_ve, eta_pl in pairs: + comb = cm._combine_yield(sympy.Float(eta_ve), sympy.Float(eta_pl)) + val = float(unwrap_expression(comb, mode="nondimensional")) + assert np.isfinite(val), f"overflow δ={delta} ({eta_ve},{eta_pl})" + assert val <= min(eta_ve, eta_pl) * (1 + 1e-9), \ + f"power-mean over-yields δ={delta}: {val} > {min(eta_ve, eta_pl)}" + + # smallest δ (largest s) is close to exact Min. + cm.yield_softness = 0.02 + for eta_ve, eta_pl in pairs: + comb = cm._combine_yield(sympy.Float(eta_ve), sympy.Float(eta_pl)) + val = float(unwrap_expression(comb, mode="nondimensional")) + m = min(eta_ve, eta_pl) + assert abs(val - m) <= 0.05 * m, f"not near Min at δ=0.02: {val} vs {m}" + + +@pytest.mark.level_1 +@pytest.mark.tier_a +def test_yield_smoother_validation(): + """yield_smoother only accepts the two known families; default is 'sqrt'.""" + mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) + v = uw.discretisation.MeshVariable("Usm", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("Psm", mesh, 1, degree=1) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + cm = uw.constitutive_models.ViscoElasticPlasticFlowModel(s.Unknowns, order=2) + s.constitutive_model = cm + assert cm.yield_smoother == "sqrt" # default family unchanged + with pytest.raises(ValueError): + cm.yield_smoother = "not_a_family" + + +@pytest.mark.level_2 +@pytest.mark.tier_a +def test_yield_softness_runtime_no_recompile(): + """Ramping δ through its constants[] atom must not trigger a JIT recompile.""" + _, stokes, cm, V_top = _build_vep("norecompile") + cm.yield_softness = 0.3 # δ>0 so the soft-min atom is exercised + dt = 0.20 * T_R + cm.Parameters.dt_elastic = dt + stokes.solve(zero_init_guess=False, timestep=dt, divergence_retries=2) + key0 = stokes._current_jit_cache_key + assert key0 is not None + + cm._get_yield_softness() + for d in (0.2, 0.1, 0.0): + cm._yield_softness_expr.sym = sympy.Float(d) + stokes._update_constants() + stokes.solve(zero_init_guess=False, timestep=dt, divergence_retries=2) + assert stokes._current_jit_cache_key == key0, "δ ramp forced a JIT recompile" + + +@pytest.mark.level_1 +@pytest.mark.tier_a +def test_dp_model_smoother_optin(): + """The non-elastic Drucker–Prager model defaults to exact hard Min and can opt + into the δ soft-min / power-mean homotopy (its new capability).""" + mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) + v = uw.discretisation.MeshVariable("Udp", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("Pdp", mesh, 1, degree=1) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + cm = uw.constitutive_models.ViscoPlasticFlowModel(s.Unknowns) + s.constitutive_model = cm + + # default: exact hard Min + assert cm.yield_mode == "min" + assert cm.yield_smoother == "sqrt" + for eta_ve, eta_pl in [(1.0, 2.0), (2.0, 1.0), (0.3, 5.0)]: + comb = cm._combine_yield(sympy.Float(eta_ve), sympy.Float(eta_pl)) + 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 + cm.yield_mode = "softmin" + 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)]: + comb = cm._combine_yield(sympy.Float(eta_ve), sympy.Float(eta_pl)) + val = float(unwrap_expression(comb, mode="nondimensional")) + assert np.isfinite(val) + assert val <= min(eta_ve, eta_pl) * (1 + 1e-9) + + # bad mode rejected + with pytest.raises(ValueError): + cm.yield_mode = "not_a_mode"