Temporal instability in electromechanics force-fiber-stretch feedback

Problem description

Coupled electromechanics simulations, in some settings, give rise to unphysical oscillations in the structural solution, that build up over time eventually leading up to the crash of the solver.

I attach some screenshots of this behavior below, for the slab test case and a ventricular simulation, both using ActiveStressRegazzoni. Both were taken at the last computed timestep, after which the linear solver for struct failed to converge.

The slab test case was slightly modified with respect to the test suite, to use active stress parameters from Fedele et al. (2023), and I have introduced a small perturbation to the fiber orientation in one point to break the symmetry (otherwise, the test would use an axis-aligned domain, with an axis-aligned structured grid). I plan to test this again on an unstructured tetrahedral grid without artificial perturbations to the fibers.

I think that @kko27 has observed a similar behavior in testing the implementation of the Land-Niederer model (#606).

Probable explanation

For realistic active stress models such as ActiveStressRegazzoni, the active tension $T_{act}$ depends on the fiber stretch $\lambda$, both directly through the expression of the active tension itself, $T_{act} = T_{act}(\mathbf{s}, \lambda)$, and indirectly through the dependence of the ODE system on $\lambda$ (so that the state vector $\mathbf{s}$ ends up depending on $\lambda$ too).

This dependence is treated explicitly. More precisely:

1. in Integrator::predictor, the active stress model is advanced by calling ActiveStress::advance_time_step passing the fiber stretch from previous time step;
2. Integrator::initiator and Integrator::corrector proceed to solve the struct equation keeping the active tension fixed.

I believe that the oscillations and blowup are due to the instability (or maybe conditional stability) of this explicit coupling.

Evidence supporting the explanation

1. In the slab test, reducing $\Delta t$ delays or removes the instability: with $\Delta t = 1 \text{ms}$, the solver blows up at $t = 193 \text{ms}$; with $\Delta t = 0.5 \text{ms}$, it blows up at $t = 230 \text{ms}$; with $\Delta t = 0.25 \text{ms}$, blows up at around $t = 340 \text{ms}$.
2. The oscillations in velocity have a period of 2 timesteps, regardless of timestep size, which is a hallmark of instability due to explicit timestepping.
3. Removing the feedback removes the issue: if I hardcode $\lambda = 1$ in the active stress model code, oscillations disappear.
4. Experimental (AI-generated) code implementing implicit coupling (the active stress model is updated within the nonlinear iterations in which struct is also solved) also removes the instabilities (at the price of much increased number of iterations and computational cost). I trust this code only up to a point, but I'll verify it, push it and update this description.

Possible solutions

1. Implicit partitioned coupling (that is, point 3 from above). It seems to increase computational cost substantially, because convergence requires many more (2-3x) iterations. I experimented adding Aitken acceleration but it didn't help much.
2. Semi-implicit discretization. The dependence of $T_{act}$ on $\lambda$ could be made implicit only in direct part (i.e. in the way $\lambda$ appears in the expression of $T_{act}$), but not in the indirect part (i.e. in the dependence of $\mathbf{s}$ on $\lambda$ through the ODE timestepping). The tangent for this term would be manageable, I think. I have not verified whether this solves the problem or not.

Open questions

1. The explicit coupling is not uncommon in electromechanics models (see e.g. Fedele et al. (2023)). It is not obvious to me why it should not work here. Some candidate (currently unverified) explanations:
   1. other solvers may use more dissipative timestepping methods, and numerical dissipation compensates for instability;
   2. other models may have further sources of damping (e.g. more viscous boundary conditions)