From f900ac5c9256cfe36877b0851ba7297cd548cd09 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sat, 17 Oct 2026 10:25:33 +1100 Subject: [PATCH 01/11] Maxwell ETD integrator: a viscosity that varies in space keeps its elastic response MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The exponential coefficients alpha = exp(-dt/tau) and phi were formed from float(viscosity)/float(modulus). A viscosity written as an expression of position (a damage band, a temperature-dependent law) raised there and the Maxwell time fell silently to None, which the coefficient update reads as the viscous limit: the model lost its elastic response wherever the viscosity was not a plain number (measured on the box fault model: a first-step end velocity of 0.007, the viscous value, against 0.061 with a constant). The coefficients are now formed symbolically when the Maxwell time is an expression, in the shared helper and at both constitutive-model callers; a constant viscosity takes the numeric path as before. Test: tests/test_1093 — a uniform-valued expression must match the constant it equals, and a weak band must not make the block viscous. Both shown failing before the fix. Underworld development team with AI support from Claude Code --- src/underworld3/constitutive_models.py | 9 ++- src/underworld3/systems/ddt.py | 21 ++++++- .../test_1093_ve_etd_expression_viscosity.py | 62 +++++++++++++++++++ 3 files changed, 88 insertions(+), 4 deletions(-) create mode 100644 tests/test_1093_ve_etd_expression_viscosity.py diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 0b44198e2..4110144a9 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -2394,7 +2394,11 @@ def _update_history_coefficients(self): mu_val = float(params.shear_modulus.sym) tau_eff = eta_val / mu_val if mu_val > 0 else sympy.oo except (TypeError, ValueError): - tau_eff = None + # a viscosity (or modulus) that varies in space: the Maxwell + # time is a field and the ETD coefficients are formed + # symbolically (ddt._update_exp_values); None here meant the + # viscous limit and lost the elastic response silently + tau_eff = sympy.sympify(params.shear_viscosity_0.sym) / sympy.sympify(params.shear_modulus.sym) try: dt_val = ( float(params.dt_elastic.sym) @@ -4014,7 +4018,8 @@ def _update_etd_coefficients(self): mu_val = float(params.shear_modulus.sym) tau_eff = eta_val / mu_val if mu_val > 0 else sympy.oo except (TypeError, ValueError): - tau_eff = None + # spatially varying: symbolic coefficients (see the isotropic model) + tau_eff = sympy.sympify(params.shear_viscosity_1.sym) / sympy.sympify(params.shear_modulus.sym) try: dt_val = ( float(params.dt_elastic.sym) diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index aca015454..bb7bb2299 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -573,11 +573,28 @@ def _update_exp_values(coeffs, dt, tau_eff): Two-element list ``[α, φ]`` to update. dt : float or None Current timestep. - tau_eff : float or None + tau_eff : float, sympy expression or None Maxwell relaxation time (η_eff/μ). When None or non-positive, - defaults to the viscous limit. + defaults to the viscous limit. A sympy expression (a viscosity that + varies in space) gives symbolic coefficients, evaluated pointwise. """ dt_f = _as_float(dt) + if (isinstance(tau_eff, sympy.Basic) and not tau_eff.is_number + and dt_f is not None and dt_f > 0.0): + # A spatially varying Maxwell time (a viscosity written as an expression + # of position or of a field): the coefficients are FIELDS, formed + # symbolically and evaluated by the kernels where they are used. Taking + # the float here failed and fell to the viscous limit, which silently + # removed the elastic response wherever the viscosity was not a plain + # number (measured: a damage band made the whole block viscous). + # tau <= 0 anywhere is read as the viscous limit there (x -> inf). + tau_pos = sympy.Max(tau_eff, sympy.Float(1e-300)) + x = sympy.Float(dt_f) / tau_pos + alpha = sympy.exp(-x) + phi = (1 - alpha) / x + coeffs[0].sym = alpha + coeffs[1].sym = phi + return tau_f = _as_float(tau_eff) if dt_f is None or tau_f is None or tau_f <= 0.0 or dt_f <= 0.0: alpha, phi = 0.0, 0.0 # viscous limit diff --git a/tests/test_1093_ve_etd_expression_viscosity.py b/tests/test_1093_ve_etd_expression_viscosity.py new file mode 100644 index 000000000..08da7ce76 --- /dev/null +++ b/tests/test_1093_ve_etd_expression_viscosity.py @@ -0,0 +1,62 @@ +"""The ETD Maxwell integrator with a spatially varying viscosity. + +The exponential (ETD) coefficients alpha = exp(-dt/tau) and phi are formed +from the Maxwell time tau = eta/mu. When the viscosity is a sympy EXPRESSION +of position (a damage band, a temperature-dependent viscosity) rather than a +number, the coefficient update could not take its float, fell back to the +viscous limit (alpha = phi = 0) and the model silently lost its elastic +response: a block under a step traction moved at the steady viscous rate +instead of taking the elastic jump (measured on the box fault study, +2026-10-17: first-step end velocity 7e-3 against 6e-2). + +The check: a viscosity written as an expression whose VALUE is uniform must +give the same first step as the constant it equals; and a genuinely varying +viscosity must not collapse onto the viscous limit. +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + + +def _first_step_end_velocity(eta_expr, dt=0.125, tag="a"): + mesh = uw.meshing.UnstructuredSimplexBox(minCoords=(0.0, 0.0), maxCoords=(3.0, 1.0), + cellSize=0.25, regular=False, qdegree=3) + x, y = mesh.X + v = uw.discretisation.MeshVariable(f"v_{tag}", mesh, 2, degree=2) + p = uw.discretisation.MeshVariable(f"p_{tag}", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel( + stokes.Unknowns, order=1, integrator="etd", objective_rate="none") + stokes.constitutive_model.Parameters.shear_modulus = 1.0 + stokes.constitutive_model.Parameters.shear_viscosity_0 = eta_expr(x, y) + stokes.bodyforce = sympy.Matrix([0.0, 0.0]) + stokes.tolerance = 1e-8 + stokes.add_dirichlet_bc((0.0, 0.0), "Left") + stokes.add_dirichlet_bc((None, 0.0), "Bottom") + stokes.add_natural_bc(sympy.Matrix([[-0.01, 0.0]]), "Right") # tension on the free end + stokes.solve(timestep=dt, zero_init_guess=True) + pts = np.array([[2.95, 0.5]]) + return float(np.asarray(uw.function.evaluate(v.sym, pts)).reshape(-1, 2)[0, 0]) + + +def test_uniform_expression_matches_the_constant(): + u_const = _first_step_end_velocity(lambda x, y: 1.0, tag="c") + # numerically uniform, but a sympy expression of position: float() of it fails + u_expr = _first_step_end_velocity(lambda x, y: 1.0 + 1.0e-12 * x, tag="e") + assert abs(u_expr - u_const) < 1e-6 * abs(u_const), (u_const, u_expr) + + +def test_varying_viscosity_keeps_the_elastic_response(): + # a weak band in the middle third: the response must stay near the elastic + # jump of the constant case (the band is small), not fall to the viscous rate, + # which is an order of magnitude smaller at this dt (tau = 1 period, dt = 1/8) + u_const = _first_step_end_velocity(lambda x, y: 1.0, tag="c2") + u_visc = _first_step_end_velocity(lambda x, y: 1.0, dt=1.0e6, tag="v") # dt >> tau: viscous limit + u_band = _first_step_end_velocity( + lambda x, y: 1.0 - 0.99 * sympy.exp(-((x - 1.5) / 0.1) ** 2) * sympy.exp(-((y - 0.5) / 0.1) ** 2), tag="b") + assert u_visc < 0.3 * u_const # the two limits are distinct here + assert u_band > 0.7 * u_const, (u_const, u_visc, u_band) # the band did not make it viscous From 73709b1f9a5cd7a33e2c5f3fc088ffde7050026a Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sun, 18 Oct 2026 10:59:18 +1100 Subject: [PATCH 02/11] ETD symbolic coefficients: floor the Maxwell time with the declared regulariser The symbolic path floored tau with a literal 1e-300; the Charter's regulariser is uw.maths.functions.vanishing. The comment records where phi = (1 - alpha)/x loses digits (x below ~1e-12) and that order 1 does not use phi. Underworld development team with AI support from Claude Code --- src/underworld3/systems/ddt.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index bb7bb2299..e7802435b 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -587,8 +587,12 @@ def _update_exp_values(coeffs, dt, tau_eff): # the float here failed and fell to the viscous limit, which silently # removed the elastic response wherever the viscosity was not a plain # number (measured: a damage band made the whole block viscous). - # tau <= 0 anywhere is read as the viscous limit there (x -> inf). - tau_pos = sympy.Max(tau_eff, sympy.Float(1e-300)) + # tau <= 0 anywhere is read as the viscous limit there (x -> inf); the + # floor is the declared vanishing regulariser, not a literal. + # At order 1 phi is replaced by alpha; at order 2 phi = (1 - alpha)/x + # loses digits only for x below ~1e-12, which a field of Maxwell times + # on the scale of the step does not reach. + tau_pos = sympy.Max(tau_eff, uw.maths.functions.vanishing) x = sympy.Float(dt_f) / tau_pos alpha = sympy.exp(-x) phi = (1 - alpha) / x From 60f19fe18fa2671f6eedf77f9e81b205a467fa0d Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sun, 18 Oct 2026 12:22:08 +1100 Subject: [PATCH 03/11] ETD Maxwell time: one helper, the step as its own constant, units honoured Review of the ported fix (four findings): - The symbolic coefficients took the step as a number, frozen into the compiled kernel: a change of step left the kernels on an earlier step's alpha (Stokes 1-8 %, Navier-Stokes 77 % on a step change). The step now enters as the dt_elastic atom, so a new value reaches the kernels as a constant. - The Maxwell time is formed from the parameters' own atoms, so the units system scales it: a dimensional modulus with a varying viscosity raised SympifyError (alpha now evaluates to exp(-0.125) in a units model). - One helper, _etd_maxwell_time, for all three callers; the experimental TransverseIsotropicVEPSplitFlowModel still had the silent viscous fallback. - The test moves to test_1053 (test_109* is outside every CI batch) and covers order 1 and 2, a changing step (fails when the step is frozen as a number), the transversely isotropic model, and a weak band against the uniform block and the band's own viscous limit; all at 1e-10 against the constant. All five fail on a build without the fix. Underworld development team with AI support from Claude Code --- src/underworld3/constitutive_models.py | 95 +++++++--------- src/underworld3/systems/ddt.py | 38 ++++--- .../test_1053_ve_etd_variable_maxwell_time.py | 101 ++++++++++++++++++ .../test_1093_ve_etd_expression_viscosity.py | 62 ----------- 4 files changed, 155 insertions(+), 141 deletions(-) create mode 100644 tests/test_1053_ve_etd_variable_maxwell_time.py delete mode 100644 tests/test_1093_ve_etd_expression_viscosity.py diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 4110144a9..eea59ddf0 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -1926,6 +1926,30 @@ def _logm_sym2(m): [k * c, (l1 + l2) / 2 + k * (b - mean)]]) +def _etd_maxwell_time(eta, mu): + """The Maxwell time tau = eta/mu the ETD coefficients are formed from: a plain + number when both parameters are numbers, otherwise the ratio of the parameters' + own expression atoms -- a field the kernels evaluate pointwise, scaled by the + units system like any other parameter.""" + if mu.sym is sympy.oo: + return sympy.oo + try: + eta_val = float(eta.sym) + mu_val = float(mu.sym) + except (TypeError, ValueError): + return eta / mu + return eta_val / mu_val if mu_val > 0 else sympy.oo + + +def _etd_timestep(dt): + """The elastic timestep as a plain number, or None when it is unset (infinite) + or cannot be read as one.""" + try: + return float(dt.sym) if dt.sym is not sympy.oo else None + except (TypeError, ValueError): + return None + + class ViscoElasticPlasticFlowModel(ViscousFlowModel): r""" Viscoelastic-plastic flow constitutive model. @@ -2386,28 +2410,10 @@ def _update_history_coefficients(self): if self.Unknowns.DFDt is None: return params = self.Parameters - if params.shear_modulus.sym is sympy.oo: - tau_eff = sympy.oo - else: - try: - eta_val = float(params.shear_viscosity_0.sym) - mu_val = float(params.shear_modulus.sym) - tau_eff = eta_val / mu_val if mu_val > 0 else sympy.oo - except (TypeError, ValueError): - # a viscosity (or modulus) that varies in space: the Maxwell - # time is a field and the ETD coefficients are formed - # symbolically (ddt._update_exp_values); None here meant the - # viscous limit and lost the elastic response silently - tau_eff = sympy.sympify(params.shear_viscosity_0.sym) / sympy.sympify(params.shear_modulus.sym) - try: - dt_val = ( - float(params.dt_elastic.sym) - if params.dt_elastic.sym is not sympy.oo - else None - ) - except (TypeError, ValueError): - dt_val = None - self.Unknowns.DFDt.update_exp_coefficients(dt_val, tau_eff) + self.Unknowns.DFDt.update_exp_coefficients( + _etd_timestep(params.dt_elastic), + _etd_maxwell_time(params.shear_viscosity_0, params.shear_modulus), + dt_symbol=params.dt_elastic) if self._order == 1: # ETD-1 reduction: φ = α makes the (φ-α)·ε̇* term zero # AND turns (1-φ)·ε̇ into (1-α)·ε̇. @@ -4010,25 +4016,10 @@ def _update_etd_coefficients(self): if self.Unknowns.DFDt is None: return params = self.Parameters - if params.shear_modulus.sym is sympy.oo: - tau_eff = sympy.oo - else: - try: - eta_val = float(params.shear_viscosity_1.sym) - mu_val = float(params.shear_modulus.sym) - tau_eff = eta_val / mu_val if mu_val > 0 else sympy.oo - except (TypeError, ValueError): - # spatially varying: symbolic coefficients (see the isotropic model) - tau_eff = sympy.sympify(params.shear_viscosity_1.sym) / sympy.sympify(params.shear_modulus.sym) - try: - dt_val = ( - float(params.dt_elastic.sym) - if params.dt_elastic.sym is not sympy.oo - else None - ) - except (TypeError, ValueError): - dt_val = None - self.Unknowns.DFDt.update_exp_coefficients(dt_val, tau_eff) + self.Unknowns.DFDt.update_exp_coefficients( + _etd_timestep(params.dt_elastic), + _etd_maxwell_time(params.shear_viscosity_1, params.shear_modulus), + dt_symbol=params.dt_elastic) def _update_history_coefficients(self): r"""Pre-solve hook — dispatches BDF, ETD (order 1/2), or hybrid. @@ -4712,24 +4703,10 @@ def _update_history_coefficients(self): if self._integrator != "etd" or self.Unknowns.DFDt is None: return super()._update_history_coefficients() params = self.Parameters - if params.shear_modulus.sym is sympy.oo: - tau_perp = sympy.oo - else: - try: - eta_val = float(params.shear_viscosity_0.sym) - mu_val = float(params.shear_modulus.sym) - tau_perp = eta_val / mu_val if mu_val > 0 else sympy.oo - except (TypeError, ValueError): - tau_perp = None - try: - dt_val = ( - float(params.dt_elastic.sym) - if params.dt_elastic.sym is not sympy.oo - else None - ) - except (TypeError, ValueError): - dt_val = None - self.Unknowns.DFDt.update_exp_coefficients(dt_val, tau_perp) + self.Unknowns.DFDt.update_exp_coefficients( + _etd_timestep(params.dt_elastic), + _etd_maxwell_time(params.shear_viscosity_0, params.shear_modulus), + dt_symbol=params.dt_elastic) def _eta_par_eff(self): """Yield-clipped ``η_∥_eff`` — same softmin/min/harmonic as parent. diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index e7802435b..23b7475c6 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -558,7 +558,7 @@ def _create_exp_coefficients(instance_id): return [alpha, phi] -def _update_exp_values(coeffs, dt, tau_eff): +def _update_exp_values(coeffs, dt, tau_eff, dt_symbol=None): r"""Update exponential-integrator coefficient values for current state. Computes :math:`\alpha = e^{-\Delta t/\tau_\mathrm{eff}}` and @@ -575,29 +575,27 @@ def _update_exp_values(coeffs, dt, tau_eff): Current timestep. tau_eff : float, sympy expression or None Maxwell relaxation time (η_eff/μ). When None or non-positive, - defaults to the viscous limit. A sympy expression (a viscosity that - varies in space) gives symbolic coefficients, evaluated pointwise. + defaults to the viscous limit. An expression (a viscosity that varies + in space) gives symbolic coefficients, evaluated pointwise. + dt_symbol : expression, optional + The timestep's own atom, used in the symbolic coefficients so that a + new step value reaches the kernels as a constant. """ dt_f = _as_float(dt) if (isinstance(tau_eff, sympy.Basic) and not tau_eff.is_number and dt_f is not None and dt_f > 0.0): - # A spatially varying Maxwell time (a viscosity written as an expression - # of position or of a field): the coefficients are FIELDS, formed - # symbolically and evaluated by the kernels where they are used. Taking - # the float here failed and fell to the viscous limit, which silently - # removed the elastic response wherever the viscosity was not a plain - # number (measured: a damage band made the whole block viscous). - # tau <= 0 anywhere is read as the viscous limit there (x -> inf); the - # floor is the declared vanishing regulariser, not a literal. - # At order 1 phi is replaced by alpha; at order 2 phi = (1 - alpha)/x - # loses digits only for x below ~1e-12, which a field of Maxwell times - # on the scale of the step does not reach. - tau_pos = sympy.Max(tau_eff, uw.maths.functions.vanishing) - x = sympy.Float(dt_f) / tau_pos + # A Maxwell time that varies in space: alpha and phi are fields the + # kernels evaluate pointwise. The step enters as its own constant atom + # (dt_symbol), so a change of step updates a constant rather than a + # number frozen into the compiled kernel. tau <= 0 is the viscous limit + # there (x -> inf); the floor is the declared vanishing regulariser. + # phi = (1 - alpha)/x, used at order 2 only, loses digits for x below + # about 1e-12. + dt_s = dt_symbol if dt_symbol is not None else sympy.Float(dt_f) + x = dt_s / sympy.Max(tau_eff, uw.maths.functions.vanishing) alpha = sympy.exp(-x) - phi = (1 - alpha) / x coeffs[0].sym = alpha - coeffs[1].sym = phi + coeffs[1].sym = (1 - alpha) / x return tau_f = _as_float(tau_eff) if dt_f is None or tau_f is None or tau_f <= 0.0 or dt_f <= 0.0: @@ -813,7 +811,7 @@ def _adjoint_support(self): """ return (False, f"{type(self).__name__} declares no linearisation") - def update_exp_coefficients(self, dt, tau_eff): + def update_exp_coefficients(self, dt, tau_eff, dt_symbol=None): r"""Set the exponential (ETD) coefficients for this step. ``self._exp_coeffs[0].sym = α = exp(-Δt/τ_eff)`` and @@ -824,7 +822,7 @@ def update_exp_coefficients(self, dt, tau_eff): Every flavour that can carry a Maxwell stress allocates these coefficients (#739); the particle flavours do not. """ - _update_exp_values(self._exp_coeffs, dt, tau_eff) + _update_exp_values(self._exp_coeffs, dt, tau_eff, dt_symbol) @property def _exp_alpha(self): diff --git a/tests/test_1053_ve_etd_variable_maxwell_time.py b/tests/test_1053_ve_etd_variable_maxwell_time.py new file mode 100644 index 000000000..48ffe6156 --- /dev/null +++ b/tests/test_1053_ve_etd_variable_maxwell_time.py @@ -0,0 +1,101 @@ +"""The ETD Maxwell integrator with a viscosity that varies in space. + +The exponential (ETD) coefficients alpha = exp(-dt/tau) and phi = (1 - alpha)/x +come from the Maxwell time tau = eta/mu. A viscosity written as an expression of +position (a damage band, a temperature-dependent law) has no float, and the +coefficients used to fall to the viscous limit (alpha = phi = 0): the model lost +its elastic response with no warning (measured on a box fault model, 2026-10-17: +first-step end velocity 7e-3, the viscous value, against 6e-2). + +Checked here against the constant the expression equals (a numerically uniform +expression, 1 + 1e-12 x, has no float but the same value): the first step at +order 1 and 2; a run whose step changes, since the step must reach the kernels as +a constant and not as a number frozen into them; the transversely isotropic model. +And a weak band is more compliant than the uniform block but far from the viscous +limit of the same banded model. +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw + +pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] + +UNIFORM = lambda x, y: 1.0 + 1.0e-12 * x + + +@pytest.fixture(autouse=True) +def _reset(): + uw.reset_default_model() + yield + uw.reset_default_model() + + +def _end_velocities(eta_expr, dts, order=1, model="isotropic"): + """A block clamped on the left under a tension on its right end: the end + velocity after each step of ``dts``.""" + uw.reset_default_model() + mesh = uw.meshing.UnstructuredSimplexBox(minCoords=(0.0, 0.0), maxCoords=(3.0, 1.0), + cellSize=0.25, regular=False, qdegree=3) + x, y = mesh.X + v = uw.discretisation.MeshVariable("v", mesh, 2, degree=2) + p = uw.discretisation.MeshVariable("p", mesh, 1, degree=1) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + if model == "isotropic": + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel( + stokes.Unknowns, order=order, integrator="etd", objective_rate="none") + stokes.constitutive_model.Parameters.shear_viscosity_0 = eta_expr(x, y) + else: + stokes.constitutive_model = uw.constitutive_models.TransverseIsotropicVEPFlowModel( + stokes.Unknowns, order=order, integrator="etd") + stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + stokes.constitutive_model.Parameters.shear_viscosity_1 = eta_expr(x, y) + stokes.constitutive_model.Parameters.director = sympy.Matrix([0, 1]) + stokes.constitutive_model.Parameters.shear_modulus = 1.0 + stokes.bodyforce = sympy.Matrix([0.0, 0.0]) + stokes.tolerance = 1e-10 + stokes.add_dirichlet_bc((0.0, 0.0), "Left") + stokes.add_dirichlet_bc((None, 0.0), "Bottom") + stokes.add_natural_bc(sympy.Matrix([[-0.01, 0.0]]), "Right") + out = [] + for k, dt in enumerate(dts): + stokes.constitutive_model.Parameters.dt_elastic = dt + stokes.solve(timestep=dt, zero_init_guess=(k == 0)) + out.append(float(np.asarray(uw.function.evaluate(v.sym, np.array([[2.95, 0.5]]))).reshape(-1, 2)[0, 0])) + return np.array(out) + + +def _close(a, b): + return np.abs(a - b).max() <= 1e-10 * np.abs(b).max() + + +@pytest.mark.parametrize("order", [1, 2]) +def test_a_uniform_expression_is_the_constant(order): + dts = [0.125, 0.125, 0.125] + const = _end_velocities(lambda x, y: 1.0, dts, order) + expr = _end_velocities(UNIFORM, dts, order) + assert _close(expr, const), (const, expr) + + +def test_a_change_of_step_reaches_the_kernels(): + dts = [0.125, 0.5, 0.05, 1.0e6] + const = _end_velocities(lambda x, y: 1.0, dts) + expr = _end_velocities(UNIFORM, dts) + assert _close(expr, const), (const, expr) + + +def test_the_transversely_isotropic_model(): + dts = [0.125, 0.5] + const = _end_velocities(lambda x, y: 1.0, dts, model="ti") + expr = _end_velocities(UNIFORM, dts, model="ti") + assert _close(expr, const), (const, expr) + + +def test_a_weak_band_keeps_the_elastic_response(): + band = lambda x, y: 1.0 - 0.9 * sympy.exp(-((x - 1.5) / 0.1) ** 2) + u_const = _end_velocities(lambda x, y: 1.0, [0.125])[0] + u_band = _end_velocities(band, [0.125])[0] + u_band_viscous = _end_velocities(band, [1.0e6])[0] + assert u_band > u_const # a weak band only adds compliance + assert u_band > 3.0 * u_band_viscous, (u_const, u_band, u_band_viscous) diff --git a/tests/test_1093_ve_etd_expression_viscosity.py b/tests/test_1093_ve_etd_expression_viscosity.py deleted file mode 100644 index 08da7ce76..000000000 --- a/tests/test_1093_ve_etd_expression_viscosity.py +++ /dev/null @@ -1,62 +0,0 @@ -"""The ETD Maxwell integrator with a spatially varying viscosity. - -The exponential (ETD) coefficients alpha = exp(-dt/tau) and phi are formed -from the Maxwell time tau = eta/mu. When the viscosity is a sympy EXPRESSION -of position (a damage band, a temperature-dependent viscosity) rather than a -number, the coefficient update could not take its float, fell back to the -viscous limit (alpha = phi = 0) and the model silently lost its elastic -response: a block under a step traction moved at the steady viscous rate -instead of taking the elastic jump (measured on the box fault study, -2026-10-17: first-step end velocity 7e-3 against 6e-2). - -The check: a viscosity written as an expression whose VALUE is uniform must -give the same first step as the constant it equals; and a genuinely varying -viscosity must not collapse onto the viscous limit. -""" -import numpy as np -import pytest -import sympy - -import underworld3 as uw - -pytestmark = [pytest.mark.level_2, pytest.mark.tier_b] - - -def _first_step_end_velocity(eta_expr, dt=0.125, tag="a"): - mesh = uw.meshing.UnstructuredSimplexBox(minCoords=(0.0, 0.0), maxCoords=(3.0, 1.0), - cellSize=0.25, regular=False, qdegree=3) - x, y = mesh.X - v = uw.discretisation.MeshVariable(f"v_{tag}", mesh, 2, degree=2) - p = uw.discretisation.MeshVariable(f"p_{tag}", mesh, 1, degree=1) - stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) - stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel( - stokes.Unknowns, order=1, integrator="etd", objective_rate="none") - stokes.constitutive_model.Parameters.shear_modulus = 1.0 - stokes.constitutive_model.Parameters.shear_viscosity_0 = eta_expr(x, y) - stokes.bodyforce = sympy.Matrix([0.0, 0.0]) - stokes.tolerance = 1e-8 - stokes.add_dirichlet_bc((0.0, 0.0), "Left") - stokes.add_dirichlet_bc((None, 0.0), "Bottom") - stokes.add_natural_bc(sympy.Matrix([[-0.01, 0.0]]), "Right") # tension on the free end - stokes.solve(timestep=dt, zero_init_guess=True) - pts = np.array([[2.95, 0.5]]) - return float(np.asarray(uw.function.evaluate(v.sym, pts)).reshape(-1, 2)[0, 0]) - - -def test_uniform_expression_matches_the_constant(): - u_const = _first_step_end_velocity(lambda x, y: 1.0, tag="c") - # numerically uniform, but a sympy expression of position: float() of it fails - u_expr = _first_step_end_velocity(lambda x, y: 1.0 + 1.0e-12 * x, tag="e") - assert abs(u_expr - u_const) < 1e-6 * abs(u_const), (u_const, u_expr) - - -def test_varying_viscosity_keeps_the_elastic_response(): - # a weak band in the middle third: the response must stay near the elastic - # jump of the constant case (the band is small), not fall to the viscous rate, - # which is an order of magnitude smaller at this dt (tau = 1 period, dt = 1/8) - u_const = _first_step_end_velocity(lambda x, y: 1.0, tag="c2") - u_visc = _first_step_end_velocity(lambda x, y: 1.0, dt=1.0e6, tag="v") # dt >> tau: viscous limit - u_band = _first_step_end_velocity( - lambda x, y: 1.0 - 0.99 * sympy.exp(-((x - 1.5) / 0.1) ** 2) * sympy.exp(-((y - 0.5) / 0.1) ** 2), tag="b") - assert u_visc < 0.3 * u_const # the two limits are distinct here - assert u_band > 0.7 * u_const, (u_const, u_visc, u_band) # the band did not make it viscous From f4c46ac99a80e314261625331062a6d14d8340a5 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sun, 18 Oct 2026 12:37:22 +1100 Subject: [PATCH 04/11] ETD phi with expm1; a vanishing modulus is the elastic limit Second review round (no blockers): at order 2 phi = (1 - e^-x)/x lost its digits as x -> 0 and was 0/0 where the modulus vanishes. phi is now -expm1(-x)/x (printed as C expm1), and the symbolic Maxwell time floors the modulus with the declared regulariser, so mu -> 0 is the elastic limit as on the numeric path. Test docstring states intent; a redundant dt_elastic assignment goes. Underworld development team with AI support from Claude Code --- src/underworld3/constitutive_models.py | 3 ++- src/underworld3/systems/ddt.py | 10 +++++----- tests/test_1053_ve_etd_variable_maxwell_time.py | 9 ++++----- 3 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index eea59ddf0..566149aee 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -1937,7 +1937,8 @@ def _etd_maxwell_time(eta, mu): eta_val = float(eta.sym) mu_val = float(mu.sym) except (TypeError, ValueError): - return eta / mu + # a vanishing modulus is the elastic limit (tau -> inf), as below + return eta / sympy.Max(mu, uw.maths.functions.vanishing) return eta_val / mu_val if mu_val > 0 else sympy.oo diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index 23b7475c6..b50c215f0 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -589,13 +589,13 @@ def _update_exp_values(coeffs, dt, tau_eff, dt_symbol=None): # (dt_symbol), so a change of step updates a constant rather than a # number frozen into the compiled kernel. tau <= 0 is the viscous limit # there (x -> inf); the floor is the declared vanishing regulariser. - # phi = (1 - alpha)/x, used at order 2 only, loses digits for x below - # about 1e-12. + # phi = (1 - e^-x)/x is written with expm1, which keeps its digits as + # x -> 0 (the elastic limit) and gives 1/x as x -> inf. + from sympy.codegen.cfunctions import expm1 dt_s = dt_symbol if dt_symbol is not None else sympy.Float(dt_f) x = dt_s / sympy.Max(tau_eff, uw.maths.functions.vanishing) - alpha = sympy.exp(-x) - coeffs[0].sym = alpha - coeffs[1].sym = (1 - alpha) / x + coeffs[0].sym = sympy.exp(-x) + coeffs[1].sym = -expm1(-x) / x return tau_f = _as_float(tau_eff) if dt_f is None or tau_f is None or tau_f <= 0.0 or dt_f <= 0.0: diff --git a/tests/test_1053_ve_etd_variable_maxwell_time.py b/tests/test_1053_ve_etd_variable_maxwell_time.py index 48ffe6156..b00fd554c 100644 --- a/tests/test_1053_ve_etd_variable_maxwell_time.py +++ b/tests/test_1053_ve_etd_variable_maxwell_time.py @@ -2,10 +2,10 @@ The exponential (ETD) coefficients alpha = exp(-dt/tau) and phi = (1 - alpha)/x come from the Maxwell time tau = eta/mu. A viscosity written as an expression of -position (a damage band, a temperature-dependent law) has no float, and the -coefficients used to fall to the viscous limit (alpha = phi = 0): the model lost -its elastic response with no warning (measured on a box fault model, 2026-10-17: -first-step end velocity 7e-3, the viscous value, against 6e-2). +position (a damage band, a temperature-dependent law) has no float, so the +coefficients are formed symbolically and evaluated pointwise; reading such a +viscosity as the viscous limit (alpha = phi = 0) would remove the elastic response +everywhere. Checked here against the constant the expression equals (a numerically uniform expression, 1 + 1e-12 x, has no float but the same value): the first step at @@ -60,7 +60,6 @@ def _end_velocities(eta_expr, dts, order=1, model="isotropic"): stokes.add_natural_bc(sympy.Matrix([[-0.01, 0.0]]), "Right") out = [] for k, dt in enumerate(dts): - stokes.constitutive_model.Parameters.dt_elastic = dt stokes.solve(timestep=dt, zero_init_guess=(k == 0)) out.append(float(np.asarray(uw.function.evaluate(v.sym, np.array([[2.95, 0.5]]))).reshape(-1, 2)[0, 0])) return np.array(out) From a7775fea86f421bd863e9f0859727a2db36f96f0 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sun, 18 Oct 2026 16:29:32 +1100 Subject: [PATCH 05/11] ETD coefficients belong to the model, not to the shared stress history The exponential integrator's alpha = exp(-dt/tau) and phi lived once on the DFDt, written by each model's pre-solve hook from its own Maxwell time. Two models sharing one history (the constituents of a multi-material model) therefore both relaxed on whichever Maxwell time was written last: 60 % error in the other material. Each VEP / TI-VEP / TI-VEP-split model now owns its (alpha, phi) (_ExponentialCoefficients), created like the DDt's and refreshed with _update_exp_values from its own Maxwell time, keeping the dt_symbol path for field Maxwell times; the ETD-1 reduction phi = alpha is applied to the model's own pair. Nothing in the library reads the DDt's coefficients any more; they are kept for outside callers of update_exp_coefficients (TODO(DESIGN) at their allocation). Single-material results are bit-identical (velocity, history, sampled stress) for VEP BDF-1/2, ETD-1/2, ETD with a field viscosity, TI-VEP BDF and ETD, and the TI split model. Underworld development team with AI support from Claude Code --- src/underworld3/constitutive_models.py | 97 ++++++++++++++++---------- src/underworld3/systems/ddt.py | 17 +++-- 2 files changed, 73 insertions(+), 41 deletions(-) diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 566149aee..92d48bb10 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -49,6 +49,7 @@ from underworld3.discretisation import MeshVariable from underworld3.systems.ddt import SemiLagrangian as SemiLagrangian_DDt from underworld3.systems.ddt import _bdf_coefficients, _as_float +from underworld3.systems.ddt import _create_exp_coefficients, _update_exp_values from underworld3.function.quantities import UWQuantity from underworld3.systems.ddt import Lagrangian as Lagrangian_DDt @@ -1951,7 +1952,39 @@ def _etd_timestep(dt): return None -class ViscoElasticPlasticFlowModel(ViscousFlowModel): +class _ExponentialCoefficients: + r"""The exponential (ETD) integrator's coefficients, owned by the model: + :math:`\alpha = e^{-\Delta t/\tau}` and :math:`\varphi = (1-\alpha)\,\tau/\Delta t` + with :math:`\tau` the model's own Maxwell time. Models that share one stress + history (the constituents of a multi-material model) each relax on theirs; + a single pair on the history would carry whichever Maxwell time was written + last.""" + + def _allocate_exp_coefficients(self): + """Allocate ``[alpha, phi]`` at the viscous limit (both zero); called once + the model has its instance number.""" + self._exp_coeffs = _create_exp_coefficients( + self._global_instance_number, owner="constitutive model") + _update_exp_values(self._exp_coeffs, None, None) + + @property + def _exp_alpha(self): + r"""The ETD :math:`\alpha` coefficient UWexpression.""" + return self._exp_coeffs[0] + + @property + def _exp_phi(self): + r"""The ETD :math:`\varphi` coefficient UWexpression.""" + return self._exp_coeffs[1] + + def _update_exp_coefficients(self, maxwell_time): + """Refresh ``[alpha, phi]`` for this step from ``maxwell_time``; the step + enters as the ``dt_elastic`` atom, so a field Maxwell time stays symbolic.""" + dt = self.Parameters.dt_elastic + _update_exp_values(self._exp_coeffs, _etd_timestep(dt), maxwell_time, dt_symbol=dt) + + +class ViscoElasticPlasticFlowModel(_ExponentialCoefficients, ViscousFlowModel): r""" Viscoelastic-plastic flow constitutive model. @@ -2144,6 +2177,8 @@ def __init__(self, unknowns, order=1, integrator: str = "bdf", super().__init__(unknowns, material_name=material_name) + self._allocate_exp_coefficients() + return class _Parameters(_ParameterBase, _ViscousParameterAlias): @@ -2401,24 +2436,20 @@ def _update_history_coefficients(self): Dispatches on ``(self._integrator, self._order)``: - ``"bdf"`` (order 1 or 2): updates BDF c-coefficients via :py:meth:`_update_bdf_coefficients`. - - ``"etd"`` order=2 (Phase B ETD-2): updates α, φ on the DDt + - ``"etd"`` order=2 (Phase B ETD-2): updates the model's own α, φ from ``τ_VE = η/μ``; forcing-history slot active. - ``"etd"`` order=1 (ETD-1): updates α, φ as for ETD-2 then forces ``φ = α`` so the ``(φ-α)·ε̇*`` term zeros out — fully L-stable single-step, no forcing-history slot needed. """ if self._integrator == "etd": - if self.Unknowns.DFDt is None: - return params = self.Parameters - self.Unknowns.DFDt.update_exp_coefficients( - _etd_timestep(params.dt_elastic), - _etd_maxwell_time(params.shear_viscosity_0, params.shear_modulus), - dt_symbol=params.dt_elastic) + self._update_exp_coefficients( + _etd_maxwell_time(params.shear_viscosity_0, params.shear_modulus)) if self._order == 1: # ETD-1 reduction: φ = α makes the (φ-α)·ε̇* term zero # AND turns (1-φ)·ε̇ into (1-α)·ε̇. - self.Unknowns.DFDt._exp_phi.sym = self.Unknowns.DFDt._exp_alpha.sym + self._exp_phi.sym = self._exp_alpha.sym else: self._update_bdf_coefficients() @@ -2529,8 +2560,8 @@ def E_eff(self): # ETD-1 (order=1): φ = α (set in _update_history_coefficients), # so the (φ-α)·ε̇* term zeros out and (1-φ)·ε̇ → (1-α)·ε̇ — same # expression tree, no separate code path needed. - alpha = DDt._exp_alpha - phi = DDt._exp_phi + alpha = self._exp_alpha + phi = self._exp_phi sigma_star = self._carried_stress_sym(0) if DDt.forcing_star is not None: edot_star = DDt.forcing_star.sym @@ -3695,7 +3726,7 @@ def _build_c_tensor(self): -class TransverseIsotropicVEPFlowModel(TransverseIsotropicFlowModel): +class TransverseIsotropicVEPFlowModel(_ExponentialCoefficients, TransverseIsotropicFlowModel): r"""Transversely isotropic viscoelastic-plastic flow model for fault mechanics. Combines the anisotropic viscosity tensor from :class:`TransverseIsotropicFlowModel` @@ -3861,6 +3892,8 @@ def __init__(self, unknowns, order=1, integrator: str = "bdf", super().__init__(unknowns, material_name=material_name) + self._allocate_exp_coefficients() + return class _Parameters(_ParameterBase, _ViscousParameterAlias): @@ -4013,20 +4046,16 @@ def effective_order(self): return self._order def _update_etd_coefficients(self): - """Refresh DDt's (α, φ) UWexpressions from τ_eff = η_1/μ.""" - if self.Unknowns.DFDt is None: - return + """Refresh the model's (α, φ) UWexpressions from τ_eff = η_1/μ.""" params = self.Parameters - self.Unknowns.DFDt.update_exp_coefficients( - _etd_timestep(params.dt_elastic), - _etd_maxwell_time(params.shear_viscosity_1, params.shear_modulus), - dt_symbol=params.dt_elastic) + self._update_exp_coefficients( + _etd_maxwell_time(params.shear_viscosity_1, params.shear_modulus)) def _update_history_coefficients(self): r"""Pre-solve hook — dispatches BDF, ETD (order 1/2), or hybrid. - BDF: refresh ``_bdf_c0..c3``. ETD-2 (order=2): refresh ``α, φ`` - on the DDt from ``η_1/μ``. ETD-1 (order=1): same plus force + BDF: refresh ``_bdf_c0..c3``. ETD-2 (order=2): refresh the model's + ``α, φ`` from ``η_1/μ``. ETD-1 (order=1): same plus force ``φ = α`` so the ``(φ-α)·ε̇*`` history term vanishes. Hybrid: refresh both BDF and ETD-2 — flux uses both per-spatial weight. """ @@ -4035,9 +4064,7 @@ def _update_history_coefficients(self): if self._order == 1: # ETD-1 reduction: φ = α zeros (φ-α)·ε̇* and turns # (1-φ)·ε̇ into (1-α)·ε̇. - DDt = self.Unknowns.DFDt - if DDt is not None: - DDt._exp_phi.sym = DDt._exp_alpha.sym + self._exp_phi.sym = self._exp_alpha.sym elif self._integrator == "hybrid": # Hybrid uses BOTH integrators with spatial blend; update # both coefficient sets each step. @@ -4162,8 +4189,8 @@ def _e_eff_for(self, integrator_mode): DDt = self.Unknowns.DFDt if integrator_mode == "etd": - alpha = DDt._exp_alpha - phi = DDt._exp_phi + alpha = self._exp_alpha + phi = self._exp_phi sigma_star = DDt.psi_star[0].sym if DDt.forcing_star is not None: edot_star = DDt.forcing_star.sym @@ -4650,8 +4677,8 @@ class TransverseIsotropicVEPSplitFlowModel(TransverseIsotropicVEPFlowModel): relaxes each channel on its proper timescale; the analytical floor on ``σ_∥`` is then ``≲ τ_y`` by construction. - Implementation choice: ``α_⊥, φ_⊥`` come from the DDt's existing - scalar ``_exp_coeffs`` (matrix viscosity is fixed and spatially + Implementation choice: ``α_⊥, φ_⊥`` are the model's own scalar + ``_exp_coeffs`` (matrix viscosity is fixed and spatially uniform — a single per-step scalar is right). ``α_∥, φ_∥`` are inlined as sympy expressions of the yield-clipped ``η_∥_eff`` so the JIT evaluates them per quadrature point (spatial heterogeneity @@ -4701,13 +4728,11 @@ def _update_history_coefficients(self): ``(α_∥, φ_∥)`` are inlined per-quadrature in :py:meth:`stress`. Picks ``τ_⊥ = η_0/μ`` (raw matrix viscosity). """ - if self._integrator != "etd" or self.Unknowns.DFDt is None: + if self._integrator != "etd": return super()._update_history_coefficients() params = self.Parameters - self.Unknowns.DFDt.update_exp_coefficients( - _etd_timestep(params.dt_elastic), - _etd_maxwell_time(params.shear_viscosity_0, params.shear_modulus), - dt_symbol=params.dt_elastic) + self._update_exp_coefficients( + _etd_maxwell_time(params.shear_viscosity_0, params.shear_modulus)) def _eta_par_eff(self): """Yield-clipped ``η_∥_eff`` — same softmin/min/harmonic as parent. @@ -4869,9 +4894,9 @@ def stress(self): DDt = self.Unknowns.DFDt E = self.Unknowns.E - # Matrix branch: scalar (α_⊥, φ_⊥) from DDt - alpha_perp = DDt._exp_alpha - phi_perp = DDt._exp_phi + # Matrix branch: the model's scalar (α_⊥, φ_⊥) + alpha_perp = self._exp_alpha + phi_perp = self._exp_phi eta_perp = params.shear_viscosity_0 eta_perp_sym = eta_perp.sym diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index b50c215f0..ad742a224 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -536,23 +536,25 @@ def _update_am_values(coeffs, effective_order, theta=0.5): coeffs[i].sym = sympy.Integer(0) -def _create_exp_coefficients(instance_id): +def _create_exp_coefficients(instance_id, owner="DDt instance"): """Create UWexpression objects for the ETD-2 coefficients ``[α, φ]``. These are named to render as ``α_exp`` and ``φ_exp`` (with the - DDt instance id appended) so they remain visually distinct from - BDF/AM coefficients in symbolic output. + owner's instance id appended) so they remain visually distinct from + BDF/AM coefficients in symbolic output. ``owner`` names what holds them + in the description: a DDt, or a constitutive model, which carries its + own because the Maxwell time they come from is the model's. """ alpha = _UWexpression( rf"{{\alpha^{{\mathrm{{exp}}}}_{{[{instance_id}]}}}}", sym=0.0, - description=f"Exp integrator α = exp(-Δt/τ) (DDt instance {instance_id})", + description=f"Exp integrator α = exp(-Δt/τ) ({owner} {instance_id})", _unique_name_generation=True, ) phi = _UWexpression( rf"{{\varphi^{{\mathrm{{exp}}}}_{{[{instance_id}]}}}}", sym=0.0, - description=f"Exp integrator φ = (1-α)/(Δt/τ) (DDt instance {instance_id})", + description=f"Exp integrator φ = (1-α)/(Δt/τ) ({owner} {instance_id})", _unique_name_generation=True, ) return [alpha, phi] @@ -794,6 +796,11 @@ def _init_coefficient_expressions(self, order, theta, with_exp): self._bdf_coeffs = _create_coefficients(order, r"c^{\mathrm{BDF}}", self.instance_number) self._am_coeffs = _create_coefficients(order, r"a^{\mathrm{AM}}", self.instance_number) if with_exp: + # TODO(DESIGN): no constitutive model reads these any more: each + # viscoelastic model owns its (alpha, phi), formed from its own + # Maxwell time, so that constituents sharing this history each relax + # on theirs. Kept for update_exp_coefficients callers outside the + # library; remove with that method once none remain. self._exp_coeffs = _create_exp_coefficients(self.instance_number) # Initialise to order-1 / viscous values _update_bdf_values(self._bdf_coeffs, 1, None, []) From d242dfaf12d4c9798b9b0cfa3c258b43f1a345be Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sun, 18 Oct 2026 16:29:45 +1100 Subject: [PATCH 06/11] An infinite elastic timestep is the viscous limit in the VE effective viscosity Before the solver sets a step, dt_elastic is infinite. With a finite modulus, eta mu dt / (c0 eta + mu dt) then evaluates to inf/inf = nan: - ViscoElasticPlasticFlowModel.ve_effective_viscosity guarded mu = inf but not dt = inf; it now returns the viscosity for either. - TransverseIsotropicVEPFlowModel: the same guard in ve_effective_viscosity and in _eta_for_tensor. Its tensor is cached (_is_setup), and a model given its parameters before it meets a solver was built at dt = inf and kept the nan tensor once the step became finite; the dt_elastic setter now resets the model when the step crosses between infinite and finite. Regression tests are in tests/test_1054_multimaterial_viscoelastic.py (next commit). Single-material results with the parameters set after assignment are unchanged bit for bit. Underworld development team with AI support from Claude Code --- src/underworld3/constitutive_models.py | 25 +++++++++++++++++-------- 1 file changed, 17 insertions(+), 8 deletions(-) diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 92d48bb10..7e2c8c223 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -2299,10 +2299,10 @@ def __init__( @property def ve_effective_viscosity(inner_self): r"""Visco-elastic effective viscosity: :math:`\eta_{\mathrm{eff}} = \frac{\eta G \Delta t}{\eta + G \Delta t}`.""" - # the dt_elastic defaults to infinity, t_relax to zero, - # so this should be well behaved in the viscous limit - - if inner_self.shear_modulus == sympy.oo: + # An infinite modulus or an infinite step (the default before the + # solver sets one) is the viscous limit; written out, either is + # inf/inf or 0*inf. + if inner_self.shear_modulus == sympy.oo or inner_self.dt_elastic.sym is sympy.oo: return inner_self.shear_viscosity_0 # BDF-k effective viscosity: eta_eff = eta*mu*dt / (c0*eta + mu*dt) @@ -3934,10 +3934,17 @@ def dt_elastic(inner_self): @dt_elastic.setter def dt_elastic(inner_self, value): + model = inner_self._owning_model + was_viscous = model._dt.sym is sympy.oo if hasattr(value, 'sym'): - inner_self._owning_model._dt.sym = value.sym + model._dt.sym = value.sym else: - inner_self._owning_model._dt.sym = value + model._dt.sym = value + # The tensor is cached, and built from the raw viscosities while the + # step is infinite (the viscous limit) but from the step's atom once + # it is finite: crossing between the two is a different tensor. + if was_viscous != (model._dt.sym is sympy.oo): + model._reset() # Viscosity limits shear_viscosity_min = api_tools.Parameter( @@ -3979,7 +3986,9 @@ def __init__(inner_self, _owning_model): def ve_effective_viscosity(inner_self): r"""VE effective viscosity using η₁ (fault-plane viscosity).""" mu_val = inner_self.shear_modulus.sym if hasattr(inner_self.shear_modulus, 'sym') else inner_self.shear_modulus - if mu_val is sympy.oo: + # infinite modulus or infinite step: the viscous limit (see the + # isotropic model) + if mu_val is sympy.oo or inner_self.dt_elastic.sym is sympy.oo: return inner_self.shear_viscosity_1 eta = inner_self.shear_viscosity_1 @@ -4338,7 +4347,7 @@ def _eta_for_tensor(self, integrator_mode, apply_yield): dt_e = self.Parameters.dt_elastic c0 = self._bdf_c0 mu_val = mu.sym if hasattr(mu, 'sym') else mu - if mu_val is sympy.oo: + if mu_val is sympy.oo or dt_e.sym is sympy.oo: eta_0 = eta_0_raw else: eta_0 = eta_0_raw * mu * dt_e / (c0 * eta_0_raw + mu * dt_e) From 45b6d584d3ab901b22814c0bd6048603e54b3538 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sun, 18 Oct 2026 17:21:33 +1100 Subject: [PATCH 07/11] Multi-material model forwards the stress-history interface to its constituents MultiMaterialConstitutiveModel took none of the stress-history interface from its constituents: it reported requires_stress_history = False (so the solver made no history), had no Parameters.dt_elastic (the VE solve raised), and never called the constituents' history hooks, so a viscoelastic constituent kept dt_elastic = inf -- nan (BDF) or a purely viscous response (ETD), silently. The composite now forwards: requires_stress_history (any constituent); order (highest) and effective_order (highest: the solver rebuilds when it changes, so it must follow the ramp of the highest-order constituent); stress_history_ddt_kwargs (union; contradictory requests raise); Parameters.dt_elastic (set on every viscoelastic constituent); _update_history_coefficients / _update_history_post_solve; history_flux (level-set sum of each constituent's history_flux or flux, sharing one weighting helper with flux); Parameters._solver and _solver_is_setup, so a constituent's parameter change rebuilds the solver's functions. Mixed integrators and orders are allowed: the history stores the composite stress, which inside each material is its own stress, and each constituent reads it through its own scheme and coefficients (its own ETD alpha, phi since the previous commits). The log-conformation refusal moves into _validate_model_compatibility, which now states what combines. tests/test_1054_multimaterial_viscoelastic.py: pure shear, two materials, closed forms away from the interface -- identical constituents (5.6e-12), two Maxwell times BDF and ETD (5.6e-4, interface smearing of the projected history), BDF-1 beside ETD-1, BDF-2 beside BDF-1, VEP beside TI-VEP, yield in one constituent, a constituent's parameter change between steps; plus the single-material TI-VEP parameters-before-assignment and infinite-step regressions. 15 of 18 fail on the base build; all pass serial and on two ranks. Subsystem doc gains the combination rules. Underworld development team with AI support from Claude Code --- .../subsystems/constitutive-models.md | 40 +++ src/underworld3/constitutive_models.py | 193 +++++++++--- tests/test_1054_multimaterial_viscoelastic.py | 276 ++++++++++++++++++ 3 files changed, 476 insertions(+), 33 deletions(-) create mode 100644 tests/test_1054_multimaterial_viscoelastic.py diff --git a/docs/developer/subsystems/constitutive-models.md b/docs/developer/subsystems/constitutive-models.md index 56170f2d3..f62991e40 100644 --- a/docs/developer/subsystems/constitutive-models.md +++ b/docs/developer/subsystems/constitutive-models.md @@ -97,6 +97,46 @@ stress_history = self.Unknowns.DFDt.psi_star[0] # Correct $\psi^*[0]$ ``` ``` +### Viscoelastic constituents in a multi-material model + +The history stores the composite stress $\sigma^* = \sum_i \phi_i \sigma_i^*$, which +inside material $i$ is that material's own stress. Each constituent advances it with +its own law, integrator and coefficients, so the composite is the set of Maxwell laws, +each integrated where it applies. `MultiMaterialConstitutiveModel` forwards the +stress-history interface to its constituents: + +| Interface | Composite behaviour | +|---|---| +| `requires_stress_history` | true when any constituent integrates a history | +| `order` / `effective_order` | the highest over the viscoelastic constituents (the history is allocated at that order; a rebuild follows the start-up ramp of the highest-order one) | +| `stress_history_ddt_kwargs` | every slot any constituent asks for; contradictory requests raise | +| `Parameters.dt_elastic` | set on every viscoelastic constituent; read from the first | +| `_update_history_coefficients`, `_update_history_post_solve` | called on every constituent | +| `history_flux` | $\sum_i \phi_i h_i$, with $h_i$ the constituent's `history_flux` (its `flux` when it has none) | +| `Parameters._solver`, `_solver_is_setup` | passed to / gathered from every constituent, so a constituent's parameter change rebuilds the solver's functions | + +**What combines.** Any mix of `ViscousFlowModel`, `ViscoElasticPlasticFlowModel`, +`TransverseIsotropicVEPFlowModel` (and the split variant), with `integrator="bdf"` +or `"etd"`, first or second order, with or without yield. Mixing integrators is +consistent because the shared history is only the stress: a BDF constituent reads +$\sigma^*$ through its own BDF weights, an exponential one through its own +$\alpha = e^{-\Delta t/\tau_i}$, $\varphi$. Each viscoelastic model therefore owns its +`(alpha, phi)` coefficients (`_ExponentialCoefficients`), formed from its own Maxwell +time; one pair on the history would carry the Maxwell time of whichever constituent +wrote it last. A BDF-1 constituent beside a BDF-2 one reads a second-order history with +a zero second weight; an ETD-1 constituent ignores the forcing slot an ETD-2 one asks +for ($\varphi = \alpha$ removes that term). + +**Refused.** A constituent with `stress_history="log_conformation"` (the history +would no longer be the stress), and constituents whose history requests contradict +each other. + +**Accuracy.** The history is a continuous projection of a stress that jumps at the +interface, so it is smeared over about a cell there. Away from the interface each +material follows its own closed form: in pure shear, a 4:1 Maxwell-time contrast on +an 8x8 mesh is within 6e-4 of the peak stress at 0.35 from the interface +(`tests/test_1054_multimaterial_viscoelastic.py`). + ### Retargeting a model to another solver Because a constitutive model reads the solver's unknowns, its parameter expressions diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 7e2c8c223..c9e208f2a 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -5001,6 +5001,24 @@ class MultiMaterialConstitutiveModel(Constitutive_Model): - Solver owns Unknowns (including :math:`D\mathbf{F}/Dt` stress history) - All constituent models share solver's Unknowns - Composite flux becomes stress history for all materials + + Viscoelastic constituents + ------------------------- + The stress history stores the composite stress + :math:`\sigma^* = \sum_i \phi_i \sigma_i^*`, which inside material :math:`i` + is that material's own stress. Each constituent advances it with its own + law, its own integrator and its own coefficients (the BDF weights, or the + exponential :math:`\alpha, \varphi` from its own Maxwell time), so any mix of + viscous, viscoelastic(-plastic) and transversely isotropic constituents, + BDF or exponential, first or second order, is the same set of Maxwell laws + each integrated where it applies. The composite forwards the elastic + timestep and the history hooks to every constituent; the history is + allocated at the highest order any constituent asks for, with every slot + (a forcing history) any of them asks for. Refused: a constituent that + stores something other than the stress (``stress_history='log_conformation'``) + and constituents whose history requests contradict each other. + Across the interface the history is the projection of a stress that jumps, + so it is smeared over about a cell there. """ def __init__( @@ -5031,10 +5049,6 @@ def __init__( # utilities/_jitextension.py and tests/test_0103_jit_rampable_constants.py. # Validate compatibility before initialization self._validate_model_compatibility(constitutive_models) - if any(getattr(m, "_stress_history", "stress") != "stress" for m in constitutive_models): - raise NotImplementedError( - "the multi-material model stores its averaged flux as a stress; a constituent " - "with stress_history='log_conformation' cannot share that history") self._material_var = material_swarmVariable self._constitutive_models = constitutive_models @@ -5053,6 +5067,39 @@ def __init__( # Composite model doesn't have its own material_name - constituents do super().__init__(unknowns, material_name=None) + class _Parameters(_ParameterBase): + """The composite has no material parameters: each constituent carries its + own. It exposes the elastic timestep, which the solver sets once per step + and every viscoelastic constituent integrates over, and passes the + solver reference on to every constituent, so that a change to any + constituent's parameters rebuilds the solver's functions.""" + + def __init__(inner_self, _owning_model): + inner_self._owning_model = _owning_model + + @property + def dt_elastic(inner_self): + """The elastic timestep the viscoelastic constituents share.""" + constituents = inner_self._owning_model._viscoelastic_constituents() + if not constituents: + raise AttributeError("no constituent of this multi-material model is viscoelastic") + return constituents[0].Parameters.dt_elastic + + @dt_elastic.setter + def dt_elastic(inner_self, value): + for constituent in inner_self._owning_model._viscoelastic_constituents(): + constituent.Parameters.dt_elastic = value + + @property + def _solver(inner_self): + return inner_self._solver_reference + + @_solver.setter + def _solver(inner_self, solver): + inner_self._solver_reference = solver + for constituent in inner_self._owning_model._constitutive_models: + constituent.Parameters._solver = solver + def _setup_shared_unknowns(self, constitutive_models, unknowns): """ Ensure all constituent models share the solver's authoritative unknowns. @@ -5080,6 +5127,11 @@ def _validate_model_compatibility(self, models: list) -> bool: - Same spatial dimension (2D/3D consistency) - Compatible flux tensor shapes - All models properly initialized + - A stress history, if any, stores the stress: the averaged flux is a + stress, so a constituent that stores another quantity (the + log-conformation) cannot share it. Mixed integrators and orders are + allowed (see the class docstring); contradictory history requests are + refused when the history is created (:attr:`stress_history_ddt_kwargs`). """ if not models: raise ValueError("At least one constitutive model required") @@ -5097,8 +5149,100 @@ def _validate_model_compatibility(self, models: list) -> bool: if not hasattr(model, "Unknowns"): raise ValueError(f"Model {i} is not properly initialized") + if any(getattr(m, "_stress_history", "stress") != "stress" for m in models): + raise NotImplementedError( + "the multi-material model stores its averaged flux as a stress; a constituent " + "with stress_history='log_conformation' cannot share that history") + return True + # ----- the stress-history interface, forwarded to the constituents ----- + + def _viscoelastic_constituents(self): + """The constituents that integrate a stress history.""" + return [m for m in self._constitutive_models if m.requires_stress_history] + + @property + def requires_stress_history(self): + """True when any constituent integrates a stress history.""" + return bool(self._viscoelastic_constituents()) + + @property + def order(self): + """The history order: the highest any constituent asks for (0 when none + integrates a history).""" + return max((m.order for m in self._viscoelastic_constituents()), default=0) + + @order.setter + def order(self, value): + """Set the order of every viscoelastic constituent (the solver does this + when its history already fixes one).""" + for constituent in self._viscoelastic_constituents(): + constituent.order = value + + @property + def effective_order(self): + """The highest order any constituent is integrating at this step. The + solver rebuilds its functions when this changes, so it must follow the + start-up ramp of the highest-order constituent; the lowest would hide + it behind a first-order one.""" + return max((m.effective_order for m in self._viscoelastic_constituents()), default=0) + + @property + def stress_history_ddt_kwargs(self): + """Every history slot any constituent asks for; contradictory requests + (one constituent asking for a slot another refuses) are an error.""" + kwargs = {} + for constituent in self._constitutive_models: + for key, value in constituent.stress_history_ddt_kwargs.items(): + if key in kwargs and kwargs[key] != value: + raise ValueError( + f"constituents of the multi-material model ask the stress history for " + f"{key}={kwargs[key]!r} and {key}={value!r}; one history cannot be both") + kwargs[key] = value + return kwargs + + def _update_history_coefficients(self): + """Pre-solve: each constituent refreshes its own integrator coefficients.""" + for constituent in self._constitutive_models: + constituent._update_history_coefficients() + + def _update_history_post_solve(self): + """Post-solve: each constituent stores its own integrator state. The + forcing history (ETD-2) is the strain rate, the same for every + constituent, so a second second-order constituent rewrites the same + values.""" + for constituent in self._constitutive_models: + constituent._update_history_post_solve() + + @property + def _solver_is_setup(self): + """Set up only while every constituent is: a change to any constituent's + parameters must rebuild what the solver compiled from the composite.""" + return self._composite_is_setup and all( + m._solver_is_setup for m in self._constitutive_models) + + @_solver_is_setup.setter + def _solver_is_setup(self, value): + self._composite_is_setup = value + for constituent in self._constitutive_models: + constituent._solver_is_setup = value + + def _level_set_weighted(self, constituent_flux): + r""":math:`\sum_i \phi_i\, f_i` over the constituents, with + :math:`f_i` = ``constituent_flux(model_i)`` and the level sets + normalised to a partition of unity when asked for.""" + weights = [self._material_var.sym[i] for i in range(self._material_var.indices)] + if self._normalize_levelsets: + total_levelset = sum(weights) + weights = [w / total_levelset for w in weights] + + fluxes = [constituent_flux(m) for m in self._constitutive_models] + combined_flux = sympy.Matrix.zeros(*fluxes[0].shape) + for weight, flux in zip(weights, fluxes): + combined_flux += weight * flux + return combined_flux + @property def flux(self): r""" @@ -5108,36 +5252,19 @@ def flux(self): all constituent models (including elastic ones) will read via ``DFDt.psi_star[0]`` in the next time step. """ - # Get reference flux shape from first model - reference_flux = self._constitutive_models[0].flux - combined_flux = sympy.Matrix.zeros(*reference_flux.shape) - - if self._normalize_levelsets: - # Compute normalization factor to ensure partition of unity - total_levelset = sum( - self._material_var.sym[i] for i in range(self._material_var.indices) - ) + return self._level_set_weighted(lambda m: m.flux) - for i in range(self._material_var.indices): - # Get normalized level-set function for material i - material_fraction = self._material_var.sym[i] / total_levelset - - # Get flux contribution from constituent model i - model_flux = self._constitutive_models[i].flux - - # Add weighted contribution to composite flux - combined_flux += material_fraction * model_flux - else: - # Use level-sets directly (assuming they already maintain partition of unity) - for i in range(self._material_var.indices): - # Get flux contribution from constituent model i - model_flux = self._constitutive_models[i].flux - - # Add weighted contribution using level-set directly - combined_flux += self._material_var.sym[i] * model_flux - - # This combined_flux will become the stress history for ALL materials - return combined_flux + @property + def history_flux(self): + r"""The stress the history carries: the level-set average of each + constituent's own history stress (its flux without any solvent part, + see :attr:`ViscoElasticPlasticFlowModel.history_flux`), or its flux + when it has none.""" + def history_part(constituent): + part = getattr(constituent, "history_flux", None) + return constituent.flux if part is None else part + + return self._level_set_weighted(history_part) @property def K(self): diff --git a/tests/test_1054_multimaterial_viscoelastic.py b/tests/test_1054_multimaterial_viscoelastic.py new file mode 100644 index 000000000..a78a97846 --- /dev/null +++ b/tests/test_1054_multimaterial_viscoelastic.py @@ -0,0 +1,276 @@ +"""Viscoelastic constituents inside the multi-material constitutive model. + +A box in pure shear, :math:`v = (\\dot\\varepsilon x, -\\dot\\varepsilon y)`, prescribed +on every wall, so the strain rate is uniform and known and each material relaxes +on its own Maxwell time :math:`\\lambda = \\eta/\\mu` towards +:math:`2\\eta\\dot\\varepsilon`. Material A fills :math:`x < 1/2` and material B +the rest; the deviatoric :math:`\\sigma_{xx}` of the committed stress history is +sampled at :math:`x = 0.15` and :math:`x = 0.85`, away from the interface. Closed +forms per step :math:`n`: + +- ETD-1, exact for constant forcing: :math:`2\\eta\\dot\\varepsilon(1 - e^{-n\\Delta t/\\lambda})`. +- BDF-1: :math:`\\sigma_{n} = \\eta_e(2\\dot\\varepsilon + \\sigma_{n-1}/(\\mu\\Delta t))` + with :math:`\\eta_e = \\eta\\mu\\Delta t/(\\eta + \\mu\\Delta t)`. + +The history is a continuous projection of the composite stress, which smears +the stress jump at the interface into the sample points. Measured 2026-10-18 +(8x8 quads, P2/P1, six steps of 0.25), error over the run relative to the peak +stress, BDF / ETD: identical constituents 5.6e-12 / 5.7e-12 (the composite is +the single material); 4:1 Maxwell times 5.6e-4 / 5.2e-4; BDF-1 beside ETD-1 +6.0e-4, BDF-2 beside BDF-1 2.3e-4; yield in one +constituent 5.7e-4 / 6.0e-4; isotropic beside transversely isotropic 1.42e-3 / +1.40e-3 (the larger jump, 1.55 against 0.47, read relative to the smaller +stress). Tolerances are about three times these. + +Also here, single-material regressions found on the way: a transversely +isotropic model given its parameters before it meets a solver (its tensor was +built, and cached, with an infinite timestep) and the viscoelastic effective +viscosity at an infinite timestep (0 * inf). +""" +import numpy as np +import pytest +import sympy + +import underworld3 as uw +from underworld3.function.expressions import unwrap_expression + +pytestmark = pytest.mark.level_2 + +EDOT = 1.0 +MU = 1.0 +DT = 0.25 +NSTEPS = 6 +POINTS = np.array([[0.15, 0.5], [0.85, 0.5]]) +# Interface smearing of the projected history at the sample points (see the +# module docstring): about three times the largest measured error. +SMEARED = 2.0e-3 +SMEARED_TI = 4.0e-3 +EXACT = 1.0e-8 + + +@pytest.fixture(autouse=True) +def _reset(): + uw.reset_default_model() + yield + uw.reset_default_model() + + +def _pure_shear_box(tag): + mesh = uw.meshing.StructuredQuadBox(elementRes=(8, 8), qdegree=3) + v = uw.discretisation.MeshVariable(f"v{tag}", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable(f"p{tag}", mesh, 1, degree=1, continuous=False) + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p, p_continuous=False) + x, y = mesh.X + for wall in ("Top", "Bottom", "Left", "Right"): + stokes.add_dirichlet_bc((EDOT * x, -EDOT * y), wall) + stokes.tolerance = 1.0e-10 + return mesh, stokes + + +def _left_right_materials(mesh, tag): + swarm = uw.swarm.Swarm(mesh) + material = uw.swarm.IndexSwarmVariable(f"M{tag}", swarm, indices=2, proxy_degree=1) + swarm.populate(fill_param=3) + material.array[:, 0, 0] = (np.asarray(swarm.coords)[:, 0] > 0.5).astype(int) + return swarm, material + + +def _vep(stokes, eta, integrator, name=None, yield_stress=None, order=1): + constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel( + stokes.Unknowns, order=order, integrator=integrator, material_name=name) + constitutive_model.Parameters.shear_viscosity_0 = eta + constitutive_model.Parameters.shear_modulus = MU + if yield_stress is not None: + constitutive_model.Parameters.yield_stress = yield_stress + constitutive_model.yield_mode = "min" + return constitutive_model + + +def _ti_vep(stokes, integrator, name=None): + """Fault-plane viscosity 0.25 on a plane at 45 degrees: pure shear in x, y is + pure shear on that plane, so sigma_xx relaxes on eta_1 / mu.""" + constitutive_model = uw.constitutive_models.TransverseIsotropicVEPFlowModel( + stokes.Unknowns, order=1, integrator=integrator, material_name=name) + _ti_parameters(constitutive_model) + return constitutive_model + + +def _ti_parameters(constitutive_model): + constitutive_model.Parameters.shear_viscosity_0 = 1.0 + constitutive_model.Parameters.shear_viscosity_1 = 0.25 + constitutive_model.Parameters.shear_modulus = MU + constitutive_model.Parameters.director = sympy.Matrix([[1, 1]]) / sympy.sqrt(2) + + +def _history_sigma_xx(stokes, nsteps=NSTEPS): + """sigma_xx of the committed history at the two sample points, per step.""" + out = [] + for _ in range(nsteps): + stokes.solve(timestep=DT) + sigma = stokes.Unknowns.DFDt.psi_star[0].sym + out.append(np.asarray(uw.function.global_evaluate(sigma[0, 0], POINTS)).reshape(-1)) + return np.array(out) + + +def _closed_form(eta, integrator, nsteps=NSTEPS, yield_stress=np.inf, order=1): + """sigma_xx after each step; a yield stress caps the pure-shear invariant |sigma_xx|. + BDF-2 (weights 3/2, -2, 1/2) starts once the history holds two committed + levels, so its first two steps are BDF-1.""" + sigma, previous, out = 0.0, 0.0, [] + for n in range(nsteps): + if integrator == "etd": + alpha = np.exp(-DT * MU / eta) + new = alpha * sigma + (1.0 - alpha) * 2.0 * eta * EDOT + elif order == 2 and n >= 2: + eta_e = eta * MU * DT / (1.5 * eta + MU * DT) + new = eta_e * (2.0 * EDOT + (2.0 * sigma - 0.5 * previous) / (MU * DT)) + else: + eta_e = eta * MU * DT / (eta + MU * DT) + new = eta_e * (2.0 * EDOT + sigma / (MU * DT)) + previous, sigma = sigma, min(new, yield_stress) + out.append(sigma) + return np.array(out) + + +def _relative_error(sim, ref): + return np.abs(sim - ref).max() / np.abs(ref).max() + + +@pytest.mark.tier_a +@pytest.mark.parametrize("integrator", ["bdf", "etd"]) +def test_identical_constituents_are_the_single_material(integrator): + mesh, stokes = _pure_shear_box(f"i{integrator}") + swarm, material = _left_right_materials(mesh, f"i{integrator}") + stokes.constitutive_model = uw.MultiMaterialConstitutiveModel( + stokes.Unknowns, material, + [_vep(stokes, 1.0, integrator, "A"), _vep(stokes, 1.0, integrator, "B")]) + sim = _history_sigma_xx(stokes) + ref = _closed_form(1.0, integrator) + for side in (0, 1): + assert _relative_error(sim[:, side], ref) < EXACT, (side, sim[:, side], ref) + + +@pytest.mark.tier_a +@pytest.mark.parametrize("integrator", ["bdf", "etd"]) +def test_each_material_relaxes_on_its_own_maxwell_time(integrator): + mesh, stokes = _pure_shear_box(f"m{integrator}") + swarm, material = _left_right_materials(mesh, f"m{integrator}") + stokes.constitutive_model = uw.MultiMaterialConstitutiveModel( + stokes.Unknowns, material, + [_vep(stokes, 1.0, integrator, "A"), _vep(stokes, 4.0, integrator, "B")]) + sim = _history_sigma_xx(stokes) + for side, eta in ((0, 1.0), (1, 4.0)): + ref = _closed_form(eta, integrator) + assert _relative_error(sim[:, side], ref) < SMEARED, (side, sim[:, side], ref) + + +@pytest.mark.tier_b +@pytest.mark.parametrize("a, b", [(("bdf", 1), ("etd", 1)), (("bdf", 2), ("bdf", 1))], + ids=["bdf1-etd1", "bdf2-bdf1"]) +def test_mixed_integrators_and_orders(a, b): + """Constituents read the same stored stress through their own schemes.""" + tag = f"x{a[0]}{a[1]}{b[0]}{b[1]}" + mesh, stokes = _pure_shear_box(tag) + swarm, material = _left_right_materials(mesh, tag) + stokes.constitutive_model = uw.MultiMaterialConstitutiveModel( + stokes.Unknowns, material, + [_vep(stokes, 1.0, a[0], "A", order=a[1]), _vep(stokes, 4.0, b[0], "B", order=b[1])]) + sim = _history_sigma_xx(stokes) + for side, eta, (integrator, order) in ((0, 1.0, a), (1, 4.0, b)): + ref = _closed_form(eta, integrator, order=order) + assert _relative_error(sim[:, side], ref) < SMEARED, (side, sim[:, side], ref) + + +@pytest.mark.tier_b +@pytest.mark.parametrize("integrator", ["bdf", "etd"]) +def test_isotropic_and_transversely_isotropic_constituents(integrator): + mesh, stokes = _pure_shear_box(f"t{integrator}") + swarm, material = _left_right_materials(mesh, f"t{integrator}") + stokes.constitutive_model = uw.MultiMaterialConstitutiveModel( + stokes.Unknowns, material, + [_vep(stokes, 1.0, integrator, "A"), _ti_vep(stokes, integrator, "B")]) + sim = _history_sigma_xx(stokes) + for side, eta in ((0, 1.0), (1, 0.25)): + ref = _closed_form(eta, integrator) + assert _relative_error(sim[:, side], ref) < SMEARED_TI, (side, sim[:, side], ref) + + +@pytest.mark.tier_b +@pytest.mark.parametrize("integrator", ["bdf", "etd"]) +def test_yield_in_one_constituent(integrator): + tau_y = 0.6 + mesh, stokes = _pure_shear_box(f"y{integrator}") + swarm, material = _left_right_materials(mesh, f"y{integrator}") + stokes.constitutive_model = uw.MultiMaterialConstitutiveModel( + stokes.Unknowns, material, + [_vep(stokes, 1.0, integrator, "A", yield_stress=tau_y), + _vep(stokes, 1.0, integrator, "B")]) + sim = _history_sigma_xx(stokes) + yielding = _closed_form(1.0, integrator, yield_stress=tau_y) + elastic = _closed_form(1.0, integrator) + assert _relative_error(sim[:, 0], yielding) < SMEARED, (sim[:, 0], yielding) + assert _relative_error(sim[:, 1], elastic) < SMEARED, (sim[:, 1], elastic) + + +@pytest.mark.tier_b +@pytest.mark.parametrize("integrator", ["bdf", "etd"]) +def test_a_constituent_change_reaches_the_solver(integrator): + """Changing a constituent's parameter between steps changes the next step.""" + mesh, stokes = _pure_shear_box(f"c{integrator}") + swarm, material = _left_right_materials(mesh, f"c{integrator}") + a, b = _vep(stokes, 1.0, integrator, "A"), _vep(stokes, 1.0, integrator, "B") + stokes.constitutive_model = uw.MultiMaterialConstitutiveModel(stokes.Unknowns, material, [a, b]) + first = _history_sigma_xx(stokes, nsteps=1)[0] + b.Parameters.shear_viscosity_0 = 4.0 + second = _history_sigma_xx(stokes, nsteps=1)[0] + # B continues from A's first-step stress, now with eta = 4 + if integrator == "etd": + alpha = np.exp(-DT * MU / 4.0) + expected_b = alpha * first[1] + (1.0 - alpha) * 2.0 * 4.0 * EDOT + else: + eta_e = 4.0 * MU * DT / (4.0 + MU * DT) + expected_b = eta_e * (2.0 * EDOT + first[1] / (MU * DT)) + expected_a = _closed_form(1.0, integrator, nsteps=2)[1] + assert abs(second[0] - expected_a) < SMEARED * expected_a, (second, expected_a) + assert abs(second[1] - expected_b) < SMEARED * expected_b, (second, expected_b) + + +# --- single-material regressions ------------------------------------------------- + + +@pytest.mark.tier_a +@pytest.mark.parametrize("integrator", ["bdf", "etd"]) +@pytest.mark.parametrize("parameters_first", [True, False]) +def test_transversely_isotropic_parameters_before_the_solver(integrator, parameters_first): + """Setting the parameters before the model meets the solver (its tensor built + at an infinite timestep) gives the same relaxation as setting them after.""" + mesh, stokes = _pure_shear_box(f"p{integrator}{int(parameters_first)}") + constitutive_model = uw.constitutive_models.TransverseIsotropicVEPFlowModel( + stokes.Unknowns, order=1, integrator=integrator) + if parameters_first: + _ti_parameters(constitutive_model) + stokes.constitutive_model = constitutive_model + else: + stokes.constitutive_model = constitutive_model + _ti_parameters(constitutive_model) + sim = _history_sigma_xx(stokes, nsteps=3) + ref = _closed_form(0.25, integrator, nsteps=3) + for side in (0, 1): + assert _relative_error(sim[:, side], ref) < EXACT, (side, sim[:, side], ref) + + +@pytest.mark.tier_a +@pytest.mark.parametrize("model_class", ["isotropic", "transverse"]) +def test_infinite_timestep_is_the_viscous_limit(model_class): + """With no timestep (dt = inf) the viscoelastic effective viscosity is the + viscosity itself, whatever the modulus: eta mu dt / (eta + mu dt) -> eta.""" + mesh, stokes = _pure_shear_box(f"d{model_class}") + if model_class == "isotropic": + constitutive_model = _vep(stokes, 2.0, "bdf") + else: + constitutive_model = _ti_vep(stokes, "bdf") + assert constitutive_model.Parameters.dt_elastic.sym is sympy.oo + eta = constitutive_model.Parameters.ve_effective_viscosity + value = sympy.sympify(unwrap_expression(eta, mode="nondimensional")) + expected = 2 if model_class == "isotropic" else sympy.Rational(1, 4) + assert value.is_finite and sympy.nsimplify(value) == expected, value From a8dd765914de8a3331e5e016a4ddb7733f5bdb39 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sun, 18 Oct 2026 19:58:31 +1100 Subject: [PATCH 08/11] Remove the stress history's own ETD coefficients Since each viscoelastic model owns its (alpha, phi), nothing in src, tests, examples or notebooks read the DDt's _exp_coeffs, _exp_alpha, _exp_phi or update_exp_coefficients. Remove them, and the with_exp switch of _init_coefficient_expressions; _create_exp_coefficients now names its owner as the constitutive model. EXPONENTIAL_VE_INTEGRATOR.md no longer says the coefficients live on the DDt. (docs/docstrings/review_queue.md still lists them; it is generated by docstring_sweep.py.) Single-material results stay bit-identical to the base build (VEP BDF-1/2, ETD-1/2, ETD with a field viscosity, TI-VEP BDF/ETD, TI split). Underworld development team with AI support from Claude Code --- .../design/EXPONENTIAL_VE_INTEGRATOR.md | 2 +- src/underworld3/constitutive_models.py | 3 +- src/underworld3/systems/ddt.py | 79 ++++++------------- 3 files changed, 24 insertions(+), 60 deletions(-) diff --git a/docs/developer/design/EXPONENTIAL_VE_INTEGRATOR.md b/docs/developer/design/EXPONENTIAL_VE_INTEGRATOR.md index 4441e81d9..01cb984bd 100644 --- a/docs/developer/design/EXPONENTIAL_VE_INTEGRATOR.md +++ b/docs/developer/design/EXPONENTIAL_VE_INTEGRATOR.md @@ -262,7 +262,7 @@ The ETD-2 result at parity-with-BDF-1-production — 1.10-1.14·τ_y at the faul ### 5. Architectural collapse landed at the parameter level -The plan envisaged Phase B with sibling classes (`MaxwellExponentialFlowModel`, TI variant) and Phase D moving integrator state onto the DDt with a strategy parameter. The collapse landed at the **constructor parameter** level instead: `ViscoElasticPlasticFlowModel(unknowns, integrator='etd')` (and the same on TI-VEP). Coefficients still live where the existing infrastructure naturally wants them (`_bdf_c0..c3` on the model, `_exp_coeffs` on the DDt) — the dispatch in `E_eff`, `viscosity`, `_build_c_tensor`, and the uniform `_update_history_*` hooks all branch on `self._integrator`. Sibling classes survive as ~10-line aliases for backwards compatibility (commit `ae79664`). +The plan envisaged Phase B with sibling classes (`MaxwellExponentialFlowModel`, TI variant) and Phase D moving integrator state onto the DDt with a strategy parameter. The collapse landed at the **constructor parameter** level instead: `ViscoElasticPlasticFlowModel(unknowns, integrator='etd')` (and the same on TI-VEP). Both coefficient sets live on the constitutive model (`_bdf_c0..c3`, and `_exp_coeffs` from `_ExponentialCoefficients`). The exponential pair first lived on the DDt, which left one Maxwell time per history: the constituents of a multi-material model sharing it all relaxed on the last one written (2026-10, `tests/test_1054_multimaterial_viscoelastic.py`); each model now forms its own — the dispatch in `E_eff`, `viscosity`, `_build_c_tensor`, and the uniform `_update_history_*` hooks all branch on `self._integrator`. Sibling classes survive as ~10-line aliases for backwards compatibility (commit `ae79664`). ### 6. Unit-handling in any new array touchpoints needs explicit care diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index c9e208f2a..caca4040e 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -1963,8 +1963,7 @@ class _ExponentialCoefficients: def _allocate_exp_coefficients(self): """Allocate ``[alpha, phi]`` at the viscous limit (both zero); called once the model has its instance number.""" - self._exp_coeffs = _create_exp_coefficients( - self._global_instance_number, owner="constitutive model") + self._exp_coeffs = _create_exp_coefficients(self._global_instance_number) _update_exp_values(self._exp_coeffs, None, None) @property diff --git a/src/underworld3/systems/ddt.py b/src/underworld3/systems/ddt.py index ad742a224..5c39887fe 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -536,25 +536,25 @@ def _update_am_values(coeffs, effective_order, theta=0.5): coeffs[i].sym = sympy.Integer(0) -def _create_exp_coefficients(instance_id, owner="DDt instance"): - """Create UWexpression objects for the ETD-2 coefficients ``[α, φ]``. - - These are named to render as ``α_exp`` and ``φ_exp`` (with the - owner's instance id appended) so they remain visually distinct from - BDF/AM coefficients in symbolic output. ``owner`` names what holds them - in the description: a DDt, or a constitutive model, which carries its - own because the Maxwell time they come from is the model's. +def _create_exp_coefficients(instance_id): + """Create UWexpression objects for the ETD coefficients ``[α, φ]`` of a + viscoelastic constitutive model. + + These are named to render as ``α_exp`` and ``φ_exp`` (with the model's + instance id appended) so they remain visually distinct from BDF/AM + coefficients in symbolic output. The model owns them because the + Maxwell time they come from is the model's. """ alpha = _UWexpression( rf"{{\alpha^{{\mathrm{{exp}}}}_{{[{instance_id}]}}}}", sym=0.0, - description=f"Exp integrator α = exp(-Δt/τ) ({owner} {instance_id})", + description=f"Exp integrator α = exp(-Δt/τ) (constitutive model {instance_id})", _unique_name_generation=True, ) phi = _UWexpression( rf"{{\varphi^{{\mathrm{{exp}}}}_{{[{instance_id}]}}}}", sym=0.0, - description=f"Exp integrator φ = (1-α)/(Δt/τ) ({owner} {instance_id})", + description=f"Exp integrator φ = (1-α)/(Δt/τ) (constitutive model {instance_id})", _unique_name_generation=True, ) return [alpha, phi] @@ -664,8 +664,6 @@ class _DDtBase(uw_object): - **History symbols**: Symbolic stores raw sympy matrices in ``psi_star``; the storage-backed flavors store variables and contribute ``.sym`` (see :meth:`_history_syms`). - - **ETD-2 exp coefficients** exist on every flavour that can carry a - viscoelastic stress (``with_exp=True``: all but the particle flavours). """ @classmethod @@ -771,8 +769,8 @@ def delta_t(self): """ return self._delta_t - def _init_coefficient_expressions(self, order, theta, with_exp): - """Create BDF/AM (and optionally ETD-2 exp) coefficient UWexpressions. + def _init_coefficient_expressions(self, order, theta): + """Create the BDF/AM coefficient UWexpressions. The coefficients are routed through PetscDS constants[] (see the module-level coefficient helpers), so order ramp-up and variable @@ -787,26 +785,15 @@ def _init_coefficient_expressions(self, order, theta, with_exp): Adams-Moulton order-1 implicitness for the startup values — ``self.theta`` for flavors that expose it, 0.5 (Crank-Nicolson) for the Lagrangian flavors that don't. - with_exp : bool - Also create the ETD-2 ``[α, φ]`` coefficients used by - Maxwell-relaxation integration; values are pushed via - PetscDSSetConstants every step in ``update_exp_coefficients`` - (every flavour but the particle ones, #739). + + The exponential (ETD) coefficients are not here: each viscoelastic + constitutive model owns its own, formed from its own Maxwell time. """ self._bdf_coeffs = _create_coefficients(order, r"c^{\mathrm{BDF}}", self.instance_number) self._am_coeffs = _create_coefficients(order, r"a^{\mathrm{AM}}", self.instance_number) - if with_exp: - # TODO(DESIGN): no constitutive model reads these any more: each - # viscoelastic model owns its (alpha, phi), formed from its own - # Maxwell time, so that constituents sharing this history each relax - # on theirs. Kept for update_exp_coefficients callers outside the - # library; remove with that method once none remain. - self._exp_coeffs = _create_exp_coefficients(self.instance_number) # Initialise to order-1 / viscous values _update_bdf_values(self._bdf_coeffs, 1, None, []) _update_am_values(self._am_coeffs, 1, theta) - if with_exp: - _update_exp_values(self._exp_coeffs, None, None) def _adjoint_support(self): """Whether this history's shift admits a discrete adjoint. @@ -818,28 +805,6 @@ def _adjoint_support(self): """ return (False, f"{type(self).__name__} declares no linearisation") - def update_exp_coefficients(self, dt, tau_eff, dt_symbol=None): - r"""Set the exponential (ETD) coefficients for this step. - - ``self._exp_coeffs[0].sym = α = exp(-Δt/τ_eff)`` and - ``self._exp_coeffs[1].sym = φ = (1-α)/(Δt/τ_eff)``, with τ_eff the - Maxwell relaxation time :math:`\eta_\mathrm{eff}/\mu`. Called by the - constitutive model, which owns τ_eff, before each solve -- peer to the - BDF/AM coefficient updates that ``update_pre_solve`` makes itself. - Every flavour that can carry a Maxwell stress allocates these - coefficients (#739); the particle flavours do not. - """ - _update_exp_values(self._exp_coeffs, dt, tau_eff, dt_symbol) - - @property - def _exp_alpha(self): - """The ETD ``α`` coefficient UWexpression.""" - return self._exp_coeffs[0] - - @property - def _exp_phi(self): - """The ETD ``φ`` coefficient UWexpression.""" - return self._exp_coeffs[1] def _note_history_shift(self, dt, **detail): """Tell the model's open step that this history advanced. @@ -1559,7 +1524,7 @@ def __init__( # Create the history list: each element is a Matrix of shape _shape. self.psi_star = [sympy.zeros(*self._shape) for _ in range(order)] - self._init_coefficient_expressions(order, self.theta, with_exp=True) + self._init_coefficient_expressions(order, self.theta) # Register with the active default model as a Snapshottable # state-bearer. Safe if no model is active. @@ -1841,7 +1806,7 @@ def __init__( ) ) - self._init_coefficient_expressions(order, self.theta, with_exp=True) + self._init_coefficient_expressions(order, self.theta) # Register with the active default model as a Snapshottable # state-bearer. Safe if no model is active. @@ -3155,7 +3120,7 @@ def __init__( self.forcing_star.remesh_policy = RemeshPolicy.CARRY self.forcing_star._remesh_managed_by = self - self._init_coefficient_expressions(order, self.theta, with_exp=True) + self._init_coefficient_expressions(order, self.theta) # Phase-2 remesh redesign: mark every DDt-owned mesh variable as # CARRY + operator-managed so the generic per-variable REMAP pass @@ -4306,7 +4271,7 @@ def __init__( # No user-settable theta on the swarm-based Lagrangian flavors: # Crank-Nicolson (0.5) is their fixed Adams-Moulton value. - self._init_coefficient_expressions(order, 0.5, with_exp=False) + self._init_coefficient_expressions(order, 0.5) dudt_swarm.populate(fill_param) # The class owns this swarm and carries the stress on it, so it also @@ -4763,7 +4728,7 @@ def _initialise_before_first_move(): # No user-settable theta on the swarm-based Lagrangian flavors: # Crank-Nicolson (0.5) is their fixed Adams-Moulton value. - self._init_coefficient_expressions(order, 0.5, with_exp=False) + self._init_coefficient_expressions(order, 0.5) # Register with the active default model as a Snapshottable # state-bearer. Safe if no model is active. @@ -5186,7 +5151,7 @@ def __init__( # V_fn evaluated at the nodes at that time, so V_fn may be any # expression (variables, ramping constants, swarm proxies). self._n_v = max(order, 2) # velocity levels the segments read - self._init_coefficient_expressions(order, self.theta, with_exp=True) + self._init_coefficient_expressions(order, self.theta) self._register_with_default_model() # The forcing (strain-rate) history the second-order exponential # integrator reads. Unlike the nodal flavour, which re-evaluates the @@ -5730,7 +5695,7 @@ def __init__( degree=1, continuous=True, varsymbol=rf"{{ F^{{\mathrm{{nodal}}}}_{{ [{inst}] }} }}") self._flux_projection = None self._n_v = 2 - self._init_coefficient_expressions(1, self.theta, with_exp=True) + self._init_coefficient_expressions(1, self.theta) self._register_with_default_model() @property From 83e8b76bd5b204db305c72c4c191c03a72b23150 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sun, 18 Oct 2026 19:58:40 +1100 Subject: [PATCH 09/11] Multi-material order setter keeps each constituent's order The solver passes its history's order to a model it is given again (constitutive_model.order = self._order). The composite forwarded that to every viscoelastic constituent, so reassigning it, or a VE_Stokes(order=2) or user history, silently turned a BDF-1 constituent into BDF-2 (error against its BDF-1 closed form 3e-4 -> 9.2e-2). The setter now leaves the constituents alone and refuses a history of lower order than a constituent integrates at. Tests (test_1054): - reassigning the composite keeps orders (2, 1) and each closed form; fails with the previous forwarding setter. - the constituent-change test now makes a structural change (a yield stress given to B between steps); the previous viscosity change was a rampable constant that needs no rebuild. Fails with the rebuild forwarding disabled (B [0.976, 1.18] against tau_y = 0.5). - viscous beside VEP, tolerance from the measured 7.9e-3. - the TI parameters-before test names its three controls. constitutive-models.md: the smearing error grows with the stress jump against the local stress; measured cases stated, including three materials in thirds (4.1e-2). Underworld development team with AI support from Claude Code --- .../subsystems/constitutive-models.md | 16 ++-- src/underworld3/constitutive_models.py | 13 +++- tests/test_1054_multimaterial_viscoelastic.py | 75 +++++++++++++++---- 3 files changed, 79 insertions(+), 25 deletions(-) diff --git a/docs/developer/subsystems/constitutive-models.md b/docs/developer/subsystems/constitutive-models.md index f62991e40..c1a161aa2 100644 --- a/docs/developer/subsystems/constitutive-models.md +++ b/docs/developer/subsystems/constitutive-models.md @@ -108,7 +108,7 @@ stress-history interface to its constituents: | Interface | Composite behaviour | |---|---| | `requires_stress_history` | true when any constituent integrates a history | -| `order` / `effective_order` | the highest over the viscoelastic constituents (the history is allocated at that order; a rebuild follows the start-up ramp of the highest-order one) | +| `order` / `effective_order` | the highest over the viscoelastic constituents (the history is allocated at that order; a rebuild follows the start-up ramp of the highest-order one). Setting `order` (the solver does, when it already holds a history) leaves each constituent's own order alone and refuses a history of lower order than a constituent needs | | `stress_history_ddt_kwargs` | every slot any constituent asks for; contradictory requests raise | | `Parameters.dt_elastic` | set on every viscoelastic constituent; read from the first | | `_update_history_coefficients`, `_update_history_post_solve` | called on every constituent | @@ -132,10 +132,16 @@ would no longer be the stress), and constituents whose history requests contradi each other. **Accuracy.** The history is a continuous projection of a stress that jumps at the -interface, so it is smeared over about a cell there. Away from the interface each -material follows its own closed form: in pure shear, a 4:1 Maxwell-time contrast on -an 8x8 mesh is within 6e-4 of the peak stress at 0.35 from the interface -(`tests/test_1054_multimaterial_viscoelastic.py`). +interface, so it is smeared there, and the error reaching a point grows with the +stress jump measured against the local stress and falls with the distance in cells. +Measured in pure shear on an 8x8 mesh, at 0.35 from the interface (about 2.8 cells), +relative to the local peak stress: two VEP materials with a 4:1 Maxwell-time +contrast, 6e-4; VEP beside TI-VEP (stresses 1.55 and 0.47), 1.4e-3 in the TI +material; a viscous material ($\eta = 4$, stress 8) beside a VEP one (stress up to +1.55), 7.9e-3 (BDF) and 7.3e-3 (ETD) in the VEP material. Three materials in thirds +(about 1.3 cells from each interface in the middle band) reach 4.1e-2. Where the +stress jump is large against the stress of interest, resolve the interface or read +the stress further from it (`tests/test_1054_multimaterial_viscoelastic.py`). ### Retargeting a model to another solver diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index caca4040e..fd55802f2 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -5174,10 +5174,15 @@ def order(self): @order.setter def order(self, value): - """Set the order of every viscoelastic constituent (the solver does this - when its history already fixes one).""" - for constituent in self._viscoelastic_constituents(): - constituent.order = value + """The solver passes the order of a history it already holds. Each + constituent keeps its own order: a history of at least the highest of + them serves every one (a lower-order constituent reads its extra levels + with zero weight), and one of lower order cannot serve them all.""" + if value < self.order: + raise ValueError( + f"the stress history has order {value}, but a constituent of the " + f"multi-material model integrates at order {self.order}; create the " + f"history at that order or lower the constituent's order") @property def effective_order(self): diff --git a/tests/test_1054_multimaterial_viscoelastic.py b/tests/test_1054_multimaterial_viscoelastic.py index a78a97846..52b058279 100644 --- a/tests/test_1054_multimaterial_viscoelastic.py +++ b/tests/test_1054_multimaterial_viscoelastic.py @@ -20,7 +20,9 @@ 6.0e-4, BDF-2 beside BDF-1 2.3e-4; yield in one constituent 5.7e-4 / 6.0e-4; isotropic beside transversely isotropic 1.42e-3 / 1.40e-3 (the larger jump, 1.55 against 0.47, read relative to the smaller -stress). Tolerances are about three times these. +stress); viscous eta = 4 beside VEP, 7.9e-3 / 7.3e-3 in the VEP material. +The smearing grows with the stress jump measured against the local stress. +Tolerances are about three times these. Also here, single-material regressions found on the way: a transversely isotropic model given its parameters before it meets a solver (its tensor was @@ -45,6 +47,7 @@ # module docstring): about three times the largest measured error. SMEARED = 2.0e-3 SMEARED_TI = 4.0e-3 +SMEARED_VISCOUS = 2.5e-2 EXACT = 1.0e-8 @@ -215,24 +218,62 @@ def test_yield_in_one_constituent(integrator): @pytest.mark.tier_b @pytest.mark.parametrize("integrator", ["bdf", "etd"]) def test_a_constituent_change_reaches_the_solver(integrator): - """Changing a constituent's parameter between steps changes the next step.""" + """A yield stress given to B between steps caps B's stress at the next step. + + The change is structural (a new term in B's flux, not a new value of a + constant), so it reaches the solve only if the composite tells the solver + to rebuild what it compiled. Measured 2026-10-18: B after the change + [0.4996, 0.5067] (BDF), within 1.4e-2 / 1.5e-2 (BDF / ETD) of tau_y = 0.5 + (the stress of the unyielded A beside it smears in), A within 7.5e-4 / + 7.8e-4 of its closed form; without the rebuild B carries on to + [0.976, 1.18].""" + tau_y = 0.5 mesh, stokes = _pure_shear_box(f"c{integrator}") swarm, material = _left_right_materials(mesh, f"c{integrator}") a, b = _vep(stokes, 1.0, integrator, "A"), _vep(stokes, 1.0, integrator, "B") stokes.constitutive_model = uw.MultiMaterialConstitutiveModel(stokes.Unknowns, material, [a, b]) - first = _history_sigma_xx(stokes, nsteps=1)[0] - b.Parameters.shear_viscosity_0 = 4.0 - second = _history_sigma_xx(stokes, nsteps=1)[0] - # B continues from A's first-step stress, now with eta = 4 - if integrator == "etd": - alpha = np.exp(-DT * MU / 4.0) - expected_b = alpha * first[1] + (1.0 - alpha) * 2.0 * 4.0 * EDOT - else: - eta_e = 4.0 * MU * DT / (4.0 + MU * DT) - expected_b = eta_e * (2.0 * EDOT + first[1] / (MU * DT)) - expected_a = _closed_form(1.0, integrator, nsteps=2)[1] - assert abs(second[0] - expected_a) < SMEARED * expected_a, (second, expected_a) - assert abs(second[1] - expected_b) < SMEARED * expected_b, (second, expected_b) + _history_sigma_xx(stokes, nsteps=2) + b.Parameters.yield_stress = tau_y + b.yield_mode = "min" + after = _history_sigma_xx(stokes, nsteps=2) + assert np.abs(after[:, 1] - tau_y).max() < 0.05 * tau_y, after[:, 1] + assert _relative_error(after[:, 0], _closed_form(1.0, integrator, nsteps=4)[2:]) < SMEARED, after[:, 0] + + +@pytest.mark.tier_a +def test_reassigning_the_composite_keeps_each_constituent_order(): + """The solver passes its history's order to a model it is given again; a + BDF-1 constituent beside a BDF-2 one stays BDF-1.""" + mesh, stokes = _pure_shear_box("o") + swarm, material = _left_right_materials(mesh, "o") + a = _vep(stokes, 1.0, "bdf", "A", order=2) + b = _vep(stokes, 4.0, "bdf", "B", order=1) + composite = uw.MultiMaterialConstitutiveModel(stokes.Unknowns, material, [a, b]) + stokes.constitutive_model = composite + stokes.constitutive_model = composite + assert (a.order, b.order) == (2, 1) + sim = _history_sigma_xx(stokes) + for side, eta, order in ((0, 1.0, 2), (1, 4.0, 1)): + ref = _closed_form(eta, "bdf", order=order) + assert _relative_error(sim[:, side], ref) < SMEARED, (side, sim[:, side], ref) + + +@pytest.mark.tier_b +@pytest.mark.parametrize("integrator", ["bdf", "etd"]) +def test_viscous_beside_viscoelastic(integrator): + """A purely viscous constituent (eta = 4: stress 8 from the first step) + beside a VEP one. The stress jump is large against A's stress, so the + smearing is larger than for two VEP constituents. Measured 2026-10-18, A / + B, BDF 7.9e-3 / 8.1e-4, ETD 7.3e-3 / 8.0e-4.""" + mesh, stokes = _pure_shear_box(f"v{integrator}") + swarm, material = _left_right_materials(mesh, f"v{integrator}") + viscous = uw.constitutive_models.ViscousFlowModel(stokes.Unknowns, material_name="B") + viscous.Parameters.shear_viscosity_0 = 4.0 + stokes.constitutive_model = uw.MultiMaterialConstitutiveModel( + stokes.Unknowns, material, [_vep(stokes, 1.0, integrator, "A"), viscous]) + sim = _history_sigma_xx(stokes) + assert _relative_error(sim[:, 0], _closed_form(1.0, integrator)) < SMEARED_VISCOUS, sim[:, 0] + assert _relative_error(sim[:, 1], np.full(NSTEPS, 2.0 * 4.0 * EDOT)) < SMEARED, sim[:, 1] # --- single-material regressions ------------------------------------------------- @@ -243,7 +284,9 @@ def test_a_constituent_change_reaches_the_solver(integrator): @pytest.mark.parametrize("parameters_first", [True, False]) def test_transversely_isotropic_parameters_before_the_solver(integrator, parameters_first): """Setting the parameters before the model meets the solver (its tensor built - at an infinite timestep) gives the same relaxation as setting them after.""" + at an infinite timestep) gives the same relaxation as setting them after. + The regression is [True-bdf]; the other three are controls (setting after, + and the exponential tensor, which never read the step, passed before).""" mesh, stokes = _pure_shear_box(f"p{integrator}{int(parameters_first)}") constitutive_model = uw.constitutive_models.TransverseIsotropicVEPFlowModel( stokes.Unknowns, order=1, integrator=integrator) From 189aa8812ffa35ac20c6ec9c4656e40b82982208 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sun, 18 Oct 2026 22:19:13 +1100 Subject: [PATCH 10/11] ETD coefficients from the scaled parameters, not their magnitudes With units, float(par.sym) is the magnitude in the parameter's own units: 1e21 Pa s / 3.2e7 Pa is a Maxwell time in seconds, against a step of 100 kyr, so alpha = 1 - 3e-12 and the model was purely elastic. The helpers now read the non-dimensional value the kernels see (dt/tau = 0.1 in the new test, for the isotropic and transversely isotropic models). Pre-existing: the base read the same magnitudes. Underworld development team with AI support from Claude Code --- src/underworld3/constitutive_models.py | 16 ++++++++-- .../test_1053_ve_etd_variable_maxwell_time.py | 29 +++++++++++++++++++ 2 files changed, 42 insertions(+), 3 deletions(-) diff --git a/src/underworld3/constitutive_models.py b/src/underworld3/constitutive_models.py index 566149aee..633a5fbd9 100644 --- a/src/underworld3/constitutive_models.py +++ b/src/underworld3/constitutive_models.py @@ -1926,6 +1926,16 @@ def _logm_sym2(m): [k * c, (l1 + l2) / 2 + k * (b - mean)]]) +def _nondimensional_float(par): + """A parameter's value as the kernels see it: scaled by the units system, as a + float. ``float(par.sym)`` would be the magnitude in the parameter's own units + (a viscosity in Pa s against a step in kyr). Raises TypeError or ValueError + when the parameter is not a number.""" + from underworld3.function.expressions import unwrap_expression + value = unwrap_expression(par, mode="nondimensional") + return float(getattr(value, "magnitude", value)) + + def _etd_maxwell_time(eta, mu): """The Maxwell time tau = eta/mu the ETD coefficients are formed from: a plain number when both parameters are numbers, otherwise the ratio of the parameters' @@ -1934,8 +1944,8 @@ def _etd_maxwell_time(eta, mu): if mu.sym is sympy.oo: return sympy.oo try: - eta_val = float(eta.sym) - mu_val = float(mu.sym) + eta_val = _nondimensional_float(eta) + mu_val = _nondimensional_float(mu) except (TypeError, ValueError): # a vanishing modulus is the elastic limit (tau -> inf), as below return eta / sympy.Max(mu, uw.maths.functions.vanishing) @@ -1946,7 +1956,7 @@ def _etd_timestep(dt): """The elastic timestep as a plain number, or None when it is unset (infinite) or cannot be read as one.""" try: - return float(dt.sym) if dt.sym is not sympy.oo else None + return _nondimensional_float(dt) if dt.sym is not sympy.oo else None except (TypeError, ValueError): return None diff --git a/tests/test_1053_ve_etd_variable_maxwell_time.py b/tests/test_1053_ve_etd_variable_maxwell_time.py index b00fd554c..dd63ecf7a 100644 --- a/tests/test_1053_ve_etd_variable_maxwell_time.py +++ b/tests/test_1053_ve_etd_variable_maxwell_time.py @@ -98,3 +98,32 @@ def test_a_weak_band_keeps_the_elastic_response(): u_band_viscous = _end_velocities(band, [1.0e6])[0] assert u_band > u_const # a weak band only adds compliance assert u_band > 3.0 * u_band_viscous, (u_const, u_band, u_band_viscous) + + +@pytest.mark.parametrize("model", ["isotropic", "ti"]) +def test_dimensional_parameters_give_the_nondimensional_maxwell_time(model): + """eta, mu and dt stated with units: the coefficients come from their scaled + values (dt/tau = 0.1 here), not their magnitudes (1e21 Pa s / 3.2e7 Pa is a + Maxwell time in seconds against a step of 100 kyr).""" + uw.get_default_model().set_reference_quantities( + length=uw.quantity(1.0, "km"), viscosity=uw.quantity(1.0e21, "Pa*s"), + time=uw.quantity(1.0, "Myr")) + mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) + v = uw.discretisation.MeshVariable("U", mesh, 2, degree=2, units="km/Myr") + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1, units="Pa") + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + if model == "isotropic": + stokes.constitutive_model = uw.constitutive_models.ViscoElasticPlasticFlowModel( + stokes.Unknowns, order=1, integrator="etd") + stokes.constitutive_model.Parameters.shear_viscosity_0 = uw.quantity(1.0e21, "Pa*s") + else: + stokes.constitutive_model = uw.constitutive_models.TransverseIsotropicVEPFlowModel( + stokes.Unknowns, order=1, integrator="etd") + stokes.constitutive_model.Parameters.shear_viscosity_0 = uw.quantity(1.0e21, "Pa*s") + stokes.constitutive_model.Parameters.shear_viscosity_1 = uw.quantity(1.0e21, "Pa*s") + stokes.constitutive_model.Parameters.director = sympy.Matrix([0, 1]) + stokes.constitutive_model.Parameters.shear_modulus = uw.quantity(1.0e21 / 3.15576e13, "Pa") + stokes.constitutive_model.Parameters.dt_elastic = uw.quantity(100.0, "kyr") + stokes.constitutive_model._update_history_coefficients() + alpha = float(stokes.DFDt._exp_alpha.sym) + assert abs(alpha - np.exp(-0.1)) < 1.0e-6, alpha From d570a8efca11e2c5a74c8f2bcb1ce850e4e1a3fa Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sun, 18 Oct 2026 22:29:04 +1100 Subject: [PATCH 11/11] Multi-material review round 2: a refused model is not assigned; units - The solver sets an explicit history order on a model before assigning it, so a composite whose constituent needs a higher-order history is refused and the solver keeps the model it had (it was left holding the refused one). - With #851's scaled Maxwell time merged in, each constituent's ETD coefficient comes from its own scaled tau when eta, mu and dt carry units (test; the units test from #851 reads alpha from the model, where this branch keeps it). - Stale comments about the history's exponential coefficients removed; the guide states how a composite and a single model take the solver's order. Underworld development team with AI support from Claude Code --- .../subsystems/constitutive-models.md | 2 +- .../cython/petsc_generic_snes_solvers.pyx | 9 ++-- src/underworld3/systems/solvers.py | 12 ++--- .../test_1053_ve_etd_variable_maxwell_time.py | 2 +- tests/test_1054_multimaterial_viscoelastic.py | 47 +++++++++++++++++++ 5 files changed, 59 insertions(+), 13 deletions(-) diff --git a/docs/developer/subsystems/constitutive-models.md b/docs/developer/subsystems/constitutive-models.md index c1a161aa2..b62c5ac29 100644 --- a/docs/developer/subsystems/constitutive-models.md +++ b/docs/developer/subsystems/constitutive-models.md @@ -108,7 +108,7 @@ stress-history interface to its constituents: | Interface | Composite behaviour | |---|---| | `requires_stress_history` | true when any constituent integrates a history | -| `order` / `effective_order` | the highest over the viscoelastic constituents (the history is allocated at that order; a rebuild follows the start-up ramp of the highest-order one). Setting `order` (the solver does, when it already holds a history) leaves each constituent's own order alone and refuses a history of lower order than a constituent needs | +| `order` / `effective_order` | the highest over the viscoelastic constituents (the history is allocated at that order; a rebuild follows the start-up ramp of the highest-order one). Setting `order` (the solver does, when it already holds a history) leaves each constituent's own order alone and refuses a history of lower order than a constituent needs; the solver then keeps the model it had. A single viscoelastic model is different: the solver sets its order to the history's, so an order-1 model given to an order-2 history integrates at order 2 | | `stress_history_ddt_kwargs` | every slot any constituent asks for; contradictory requests raise | | `Parameters.dt_elastic` | set on every viscoelastic constituent; read from the first | | `_update_history_coefficients`, `_update_history_post_solve` | called on every constituent | diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index 5e3c23dda..274c8f349 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -4904,13 +4904,14 @@ class SolverBaseClass(uw_object): ### checking if it's an instance - it will need to be reset if isinstance(model_or_class, uw.constitutive_models.Constitutive_Model): - self._constitutive_model = model_or_class - self._constitutive_model.Unknowns = self.Unknowns - self._constitutive_model._solver_is_setup = False # Only override the constitutive model's order if the solver has # an explicit VE order (> 0). Otherwise preserve the model's default. + # Done first: a model that refuses the order is never assigned. if self._order > 0: - self._constitutive_model.order = self._order + model_or_class.order = self._order + self._constitutive_model = model_or_class + self._constitutive_model.Unknowns = self.Unknowns + self._constitutive_model._solver_is_setup = False # Establish bidirectional reference so parameter changes can propagate to solver self._constitutive_model.Parameters._solver = self diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 605b11814..8bc6790c6 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -1751,10 +1751,9 @@ def _stress_history_advance(self, timestep, verbose=False, evalf=False): self.DFDt.update_pre_solve(timestep, verbose=verbose, evalf=evalf, store_result=False) - # Uniform pre-solve coefficient hook: VEP delegates to - # _update_bdf_coefficients(); MaxwellExponentialFlowModel updates - # α, φ on the DDt via _update_exp_coefficients(). No isinstance - # checks at the solver layer. + # Uniform pre-solve coefficient hook: each model refreshes its own BDF + # or exponential (α, φ) coefficients. No isinstance checks at the + # solver layer. self.constitutive_model._update_history_coefficients() def _stress_history_post_solve(self, timestep, verbose=False, evalf=False): @@ -1847,9 +1846,8 @@ def _create_stress_history_ddt(self, order=2): f"{type(cm).__name__} asks its stress history for " f"{sorted(ddt_kwargs)}, which the particle Lagrangian flavour does " "not provide; use stress_transport='semi_lagrangian' for it.") - # Order 1 BDF only for now: the particle flavour has no exponential - # coefficients (it is built with_exp=False), and order 2 is not yet - # validated. Refuse cleanly rather than crash inside the first solve. + # Order 1 BDF only for now: the exponential integrator is not + # validated on the particle flavour, and neither is order 2. Refuse cleanly rather than crash inside the first solve. if getattr(cm, "_integrator", "bdf") != "bdf": raise NotImplementedError( "the particle Lagrangian stress history supports the BDF " diff --git a/tests/test_1053_ve_etd_variable_maxwell_time.py b/tests/test_1053_ve_etd_variable_maxwell_time.py index dd63ecf7a..e0c2b2a65 100644 --- a/tests/test_1053_ve_etd_variable_maxwell_time.py +++ b/tests/test_1053_ve_etd_variable_maxwell_time.py @@ -125,5 +125,5 @@ def test_dimensional_parameters_give_the_nondimensional_maxwell_time(model): stokes.constitutive_model.Parameters.shear_modulus = uw.quantity(1.0e21 / 3.15576e13, "Pa") stokes.constitutive_model.Parameters.dt_elastic = uw.quantity(100.0, "kyr") stokes.constitutive_model._update_history_coefficients() - alpha = float(stokes.DFDt._exp_alpha.sym) + alpha = float(stokes.constitutive_model._exp_alpha.sym) assert abs(alpha - np.exp(-0.1)) < 1.0e-6, alpha diff --git a/tests/test_1054_multimaterial_viscoelastic.py b/tests/test_1054_multimaterial_viscoelastic.py index 52b058279..a11dd05ab 100644 --- a/tests/test_1054_multimaterial_viscoelastic.py +++ b/tests/test_1054_multimaterial_viscoelastic.py @@ -317,3 +317,50 @@ def test_infinite_timestep_is_the_viscous_limit(model_class): value = sympy.sympify(unwrap_expression(eta, mode="nondimensional")) expected = 2 if model_class == "isotropic" else sympy.Rational(1, 4) assert value.is_finite and sympy.nsimplify(value) == expected, value + + +@pytest.mark.tier_a +def test_a_refused_composite_is_not_assigned(): + """A history of order 1 cannot serve a BDF-2 constituent: the solver refuses + the composite and keeps the model it had.""" + mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) + v = uw.discretisation.MeshVariable("vr", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("pr", mesh, 1, degree=1) + with pytest.warns(DeprecationWarning): + stokes = uw.systems.VE_Stokes(mesh, velocityField=v, pressureField=p, order=1) + before = stokes.constitutive_model + swarm, material = _left_right_materials(mesh, "r") + a = _vep(stokes, 1.0, "bdf", "A", order=2) + b = _vep(stokes, 4.0, "bdf", "B", order=1) + composite = uw.MultiMaterialConstitutiveModel(stokes.Unknowns, material, [a, b]) + with pytest.raises(ValueError, match="order"): + stokes.constitutive_model = composite + assert stokes.constitutive_model is before + + +@pytest.mark.tier_a +def test_dimensional_constituents_take_the_scaled_maxwell_time(): + """eta, mu and dt with units: each constituent's ETD coefficient comes from + its scaled Maxwell time (dt/tau = 0.1 and 0.025), not the parameters' + magnitudes in their own units.""" + uw.get_default_model().set_reference_quantities( + length=uw.quantity(1.0, "km"), viscosity=uw.quantity(1.0e21, "Pa*s"), + time=uw.quantity(1.0, "Myr")) + mesh = uw.meshing.StructuredQuadBox(elementRes=(4, 4)) + v = uw.discretisation.MeshVariable("vu", mesh, mesh.dim, degree=2, units="km/Myr") + p = uw.discretisation.MeshVariable("pu", mesh, 1, degree=1, units="Pa") + stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + swarm, material = _left_right_materials(mesh, "u") + constituents = [] + for name, eta in (("A", 1.0e21), ("B", 4.0e21)): + constituent = uw.constitutive_models.ViscoElasticPlasticFlowModel( + stokes.Unknowns, order=1, integrator="etd", material_name=name) + constituent.Parameters.shear_viscosity_0 = uw.quantity(eta, "Pa*s") + constituent.Parameters.shear_modulus = uw.quantity(1.0e21 / 3.15576e13, "Pa") + constituent.Parameters.dt_elastic = uw.quantity(100.0, "kyr") + constituents.append(constituent) + stokes.constitutive_model = uw.MultiMaterialConstitutiveModel( + stokes.Unknowns, material, constituents) + stokes.constitutive_model._update_history_coefficients() + alphas = [float(c._exp_alpha.sym) for c in constituents] + assert np.allclose(alphas, np.exp([-0.1, -0.025]), rtol=0, atol=1.0e-6), alphas