From 6b2611dd03903aa20a19f9443c69df736be250a2 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Tue, 13 Oct 2026 11:17:00 -0700 Subject: [PATCH 1/8] Find the checkpoint vector wherever the writer put it (test_1079) test_1079's `_checkpoint_fields` helper asserted the field's saved vector is at `uw_checkpoint/`. #756 moved the reload payload under /restart/petsc, where the same vector sits one group deeper as `.../dms//vecs//`, so the assertion fails on every current file and development has been red since #708 merged. Nothing was broken. By the time that line runs, both the write and the resume phases have already returned 0 and the saved metadata has already compared equal -- the restart works. The test is marked `xfail(strict=True, raises=SymbolicDiskContinuationMismatch)`, so only that one exception counts as the expected failure and an AssertionError from the layout check surfaces as a genuine failure, which is why a cosmetic relocation reads as a broken restart. The helper now accepts either spelling, so a file written by either version stays readable to the test. This is the same assumption that regressed #640 in the reader (fixed in #847): the payload is found by structure, not by one hard-coded path. Proven on one build: with the single-spelling assertion restored the test fails at line 34; with this change it reaches its intended strict xfail. Checkpoint and snapshot suite: 122 passed, 1 xfailed. How this reached development: #708 added this test and its CI was green, but that run predates #756 by a day. #708 then merged twelve hours after #756 without a re-run, so the layout change and the assertion never met in CI until afterwards. Underworld development team with AI support from Claude Code Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01E87Q7KrpapxeQiLD1RiNXv --- tests/test_1079_diffusion_disk_restart.py | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/tests/test_1079_diffusion_disk_restart.py b/tests/test_1079_diffusion_disk_restart.py index 2bce0e560..aa2afb1e7 100644 --- a/tests/test_1079_diffusion_disk_restart.py +++ b/tests/test_1079_diffusion_disk_restart.py @@ -31,7 +31,19 @@ def _checkpoint_fields(wrapper): data.visititems(lambda key, obj: arrays.update({key: obj[()]}) if isinstance(obj, h5py.Dataset) else None) assert arrays, f"Missing checkpoint datasets for {name}" - assert f"uw_checkpoint/{name}" in arrays + # The field's saved vector, wherever the writer put it. It sat + # at uw_checkpoint/ until #756 moved the payload under + # /restart/petsc, where the same vector is one group deeper as + # .../dms//vecs//. Pinning one spelling made + # this test fail on a relocation that broke nothing -- both + # phases had already restarted and the metadata already + # compared equal by the time this line ran -- and it is the + # same assumption that regressed #640 in the reader (#847). + assert any( + key == f"uw_checkpoint/{name}" + or key.endswith(f"dms/{name}/vecs/{name}/{name}") + for key in arrays + ), sorted(arrays) fields[mesh_name, name] = arrays assert fields return fields From 84f2075c6b05ee555cecc1d6bf4c2ebccab46758 Mon Sep 17 00:00:00 2001 From: Juan Carlos Graciosa Date: Wed, 14 Oct 2026 11:04:58 +1100 Subject: [PATCH 2/8] ci: build a version-tagged container on the release tag MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit docker-image.yaml matched branch pushes only, so v3.2.0 produced no version-tagged laptop image and :v3.2.0 was added by hand afterwards. `tags: ['v*']` makes a release tag build like any other push; the semver entries already in metadata-action then emit 3.2.0 and 3.2. The `release: published` trigger in both workflows has never fired — release.yml creates the Release with GITHUB_TOKEN, which raises no workflow events. The comment in hpc-uw3-image.yml said the opposite. Co-Authored-By: Claude Opus 5 (1M context) --- .github/workflows/docker-image.yaml | 4 ++++ .github/workflows/hpc-uw3-image.yml | 6 +++--- 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/.github/workflows/docker-image.yaml b/.github/workflows/docker-image.yaml index ad3a3d300..2c6c9ad8a 100644 --- a/.github/workflows/docker-image.yaml +++ b/.github/workflows/docker-image.yaml @@ -9,6 +9,10 @@ name: "GHCR: Command-Line Container" on: push: branches: [main, development] + # release.yml creates the Release with GITHUB_TOKEN, which raises no + # workflow events, so `release:` below never fires. This is the path + # that produces a version-tagged image. + tags: ['v*'] paths: # Only rebuild when these files change - 'container/Containerfile' diff --git a/.github/workflows/hpc-uw3-image.yml b/.github/workflows/hpc-uw3-image.yml index fbac9eef0..9f9116329 100644 --- a/.github/workflows/hpc-uw3-image.yml +++ b/.github/workflows/hpc-uw3-image.yml @@ -21,9 +21,9 @@ on: - 'pyproject.toml' - 'setup.py' - 'src/**' - # Releases are the guaranteed path for a tagged image: the `tags:` filter above - # is combined with `paths:`, so a tag whose diff touches nothing in that list - # would not build on its own. + # Never fires in practice: release.yml creates the Release with + # GITHUB_TOKEN, which raises no workflow events. The `tags:` filter + # above is what builds the tagged image. Kept for a hand-published Release. release: types: [published] workflow_dispatch: From a7775fea86f421bd863e9f0859727a2db36f96f0 Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sun, 18 Oct 2026 16:29:32 +1100 Subject: [PATCH 3/8] 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 4/8] 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 5/8] 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 6/8] 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 7/8] 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 d570a8efca11e2c5a74c8f2bcb1ce850e4e1a3fa Mon Sep 17 00:00:00 2001 From: Louis Moresi Date: Sun, 18 Oct 2026 22:29:04 +1100 Subject: [PATCH 8/8] 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