ETD Maxwell integrator: a viscosity that varies in space keeps its elastic response

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
 

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