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/docs/developer/subsystems/constitutive-models.md b/docs/developer/subsystems/constitutive-models.md index 56170f2d3..b62c5ac29 100644 --- a/docs/developer/subsystems/constitutive-models.md +++ b/docs/developer/subsystems/constitutive-models.md @@ -97,6 +97,52 @@ 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). 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 | +| `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 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 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 0b44198e2..f16971949 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 @@ -1926,7 +1927,73 @@ def _logm_sym2(m): [k * c, (l1 + l2) / 2 + k * (b - mean)]]) -class ViscoElasticPlasticFlowModel(ViscousFlowModel): +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' + 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 = _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) + 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 _nondimensional_float(dt) if dt.sym is not sympy.oo else None + except (TypeError, ValueError): + return None + + +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) + _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. @@ -2119,6 +2186,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): @@ -2239,10 +2308,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) @@ -2376,38 +2445,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 - 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): - tau_eff = 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_eff) + 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() @@ -2518,8 +2569,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 @@ -3684,7 +3735,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` @@ -3850,6 +3901,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): @@ -3890,10 +3943,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( @@ -3935,7 +3995,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 @@ -4002,34 +4064,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 - 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): - tau_eff = 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_eff) + 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. """ @@ -4038,9 +4082,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. @@ -4165,8 +4207,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 @@ -4314,7 +4356,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) @@ -4653,8 +4695,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 @@ -4704,27 +4746,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 - 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._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. @@ -4886,9 +4912,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 @@ -4984,6 +5010,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__( @@ -5014,10 +5058,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 @@ -5036,6 +5076,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. @@ -5063,6 +5136,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") @@ -5080,8 +5158,105 @@ 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): + """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): + """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""" @@ -5091,36 +5266,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) + return self._level_set_weighted(lambda m: m.flux) - 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) - ) - - 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/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/ddt.py b/src/underworld3/systems/ddt.py index aca015454..5c39887fe 100644 --- a/src/underworld3/systems/ddt.py +++ b/src/underworld3/systems/ddt.py @@ -537,28 +537,30 @@ def _update_am_values(coeffs, effective_order, theta=0.5): def _create_exp_coefficients(instance_id): - """Create UWexpression objects for the ETD-2 coefficients ``[α, φ]``. + """Create UWexpression objects for the ETD coefficients ``[α, φ]`` of a + viscoelastic constitutive model. - 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. + 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/τ) (DDt instance {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/τ) (DDt instance {instance_id})", + description=f"Exp integrator φ = (1-α)/(Δt/τ) (constitutive model {instance_id})", _unique_name_generation=True, ) 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 @@ -573,11 +575,30 @@ 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. 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 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 - 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) + 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: alpha, phi = 0.0, 0.0 # viscous limit @@ -643,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 @@ -750,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 @@ -766,21 +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: - 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. @@ -792,28 +805,6 @@ def _adjoint_support(self): """ return (False, f"{type(self).__name__} declares no linearisation") - def update_exp_coefficients(self, dt, tau_eff): - 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) - - @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. @@ -1533,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. @@ -1815,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. @@ -3129,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 @@ -4280,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 @@ -4737,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. @@ -5160,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 @@ -5704,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 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 new file mode 100644 index 000000000..e0c2b2a65 --- /dev/null +++ b/tests/test_1053_ve_etd_variable_maxwell_time.py @@ -0,0 +1,129 @@ +"""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, 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 +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.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) + + +@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.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 new file mode 100644 index 000000000..a11dd05ab --- /dev/null +++ b/tests/test_1054_multimaterial_viscoelastic.py @@ -0,0 +1,366 @@ +"""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); 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 +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 +SMEARED_VISCOUS = 2.5e-2 +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): + """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]) + _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 ------------------------------------------------- + + +@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. + 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) + 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 + + +@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