Add Regazzoni active stress model

[samibismar]

Current situation

This PR adds the RegazzoniActiveStress model, implementing the RDQ20-MF formulation described by Regazzoni, Dedè, and Quarteroni (2020), using the ActiveStress framework introduced in PR #578.

The implementation includes:

• the 20-state RDQ20-MF system, consisting of 16 regulatory-unit probabilities and 4 crossbridge moments;
• model-specific XML parameter parsing and MPI distribution;
• registration with the ActiveStress model factory;
• calcium-, sarcomere-length-, and velocity-dependent state evolution;
• active-tension computation;
• and an electromechanics regression test.

The local active-stress interface is extended to provide fiber stretch to each model. RDQ20-MF first computes the paper tension

$$ \widetilde T_{\mathrm{act}} = a_{\mathrm{XB}} \left(\mu_P^1+\mu_N^1\right)\phi(SL), $$

then returns

$$ T_{\mathrm{code}} = \frac{\widetilde T_{\mathrm{act}}}{\lambda_f}. $$

This normalization is required because svMultiPhysics assembles the returned coefficient as a second Piola–Kirchhoff stress contribution:

$$ \mathbf S_{\mathrm{act}} = T_{\mathrm{code}}\ \mathbf f_0\otimes\mathbf f_0. $$

After conversion to first Piola–Kirchhoff stress, this gives

$$\mathbf P_{\mathrm{act}} = \widetilde T_{\mathrm{act}}\ \mathbf f\otimes\mathbf f_0, $$

consistent with the RDQ20 paper convention.

Related to #594.

Release Notes

• Add the RDQ20-MF active-stress model.
• Add RDQ20-MF XML parameters and MPI parameter distribution.
• Extend the local active-stress interface to receive fiber stretch.
• Normalize the RDQ20-MF returned coefficient for the existing PK2 mechanics assembly.
• Add the slab_RDQ20MF electromechanics regression case.

Documentation

The RDQ20-MF class documentation describes:

• the model states and parameters
• the distinction between the paper tension $\widetilde T_{\mathrm{act}}$ and the returned PK2 coefficient $T_{\mathrm{code}}$
• the fiber-stretch normalization
• and the exception raised for zero or near-zero fiber stretch

The slab_RDQ20MF README documents:

• the coupled electromechanics setup
• the svMultiPhysics fiber/sheet/sheet-normal directional weighting
• the meaning of the stored Active_tension_* fields
• and the reference-validation procedure

Testing

• The project build completed successfully.
• The slab_RDQ20MF pytest case passed for all three configured process-count variants.
• The RDQ20-MF state evolution and paper tension were validated against the authors' C++ reference implementation at commit 26f05df.
• After normalization, the reconstructed paper tension

$$ \widetilde{T}_{\mathrm{act}} = \lambda_f T_{\mathrm{code}} $$

was compared node-by-node against the same reference implementation and agreed to approximately machine precision.

• Additional checks at $\lambda_f = 0.8$, $1.0$, and $1.2$ confirmed

$$ \lambda_f T_{\mathrm{code}} = \widetilde{T}_{\mathrm{act}} $$

to floating-point precision.

• A local full-suite run completed with 206 tests passing and 27 skipped. The skipped tests require optional PETSc or Trilinos build configurations.
• Twelve sv0D coupling tests could not run because the local svZeroDSolver shared library was unavailable. The same twelve tests fail with the same missing-library error on main in the same environment.

References

Regazzoni, F., Dede', L., & Quarteroni, A.
“Biophysically detailed mathematical models of multiscale cardiac active mechanics.”
PLOS Computational Biology, 2020.
https://doi.org/10.1371/journal.pcbi.1008294

Code of Conduct & Contributing Guidelines

• I agree to follow the Code of Conduct and Contributing Guidelines.

[claude]

Claude Code Review

This pull request is from a fork — automated review is disabled. A repository maintainer can comment @claude review to run a one-time review.

[codecov]

Codecov Report

❌ Patch coverage is 96.87500% with 7 lines in your changes missing coverage. Please review.
✅ Project coverage is 72.73%. Comparing base (b5fbed0) to head (577f7a6).
⚠️ Report is 1 commits behind head on main.

Files with missing lines	Patch %	Lines
Code/Source/solver/active_stress_regazzoni.cpp	97.39%	5 Missing ⚠️
Code/Source/solver/active_stress_regazzoni.h	96.55%	1 Missing ⚠️
Code/Source/solver/active_stress_uniform_steady.h	0.00%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #596      +/-   ##
==========================================
+ Coverage   72.54%   72.73%   +0.18%     
==========================================
  Files         252      255       +3     
  Lines       39032    39290     +258     
  Branches     6678     6731      +53     
==========================================
+ Hits        28317    28576     +259     
+ Misses      10480    10472       -8     
- Partials      235      242       +7     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[ktbolt]

@samibismar May I suggest using a more meaningful name than RDQ20MF ? Maybe just MeanField perhaps qualified by what it actually represents rather than authors and publication date ?

[samibismar]

@ktbolt Would MeanFieldActiveStress be an appropriate class name, with MeanField as the factory label used in the XML input, for example <Model>MeanField</Model>?

I would retain “RDQ20-MF” in the documentation and citations to identify the specific Regazzoni–Dedè–Quarteroni formulation being implemented.

[ktbolt]

@samibismar I think MeanFieldActiveStress is fine and just document it what the implementation is based on.

[michelebucelli]

@samibismar @ktbolt Sorry to comment on this this late, but: I'm not sure that MeanFieldActiveStress is a descriptive name for the model. Mean-field based models are more a general category than a description of this specific model, and the fact that the model is based on the mean-field formalism is, I think, a relatively secondary aspect of it.

As far as I know, the standard way of referring to these cellular models is by the name of the author (and sometimes year, if the same authors have published multiple models), since (I think) it is often impossible to come up with a unique descriptive name. This is the same as with ionic models (and indeed we have the classes AlievPanfilov, BuenoOrovio and so on).

If you feel the acronym RDQ is a bit too obscure, I wouldn't mind using Regazzoni instead (or RegazzoniActiveStress).

[samibismar]

@michelebucelli That makes sense, thanks. I agree that MeanFieldActiveStress is probably too broad if mean-field models are a general category rather than a unique description of this formulation.

I think RegazzoniActiveStress for the C++ class, with Regazzoni as the model/XML label, would fit the existing naming pattern while remaining clear in the active-stress code. I would continue to identify it as the RDQ20-MF formulation of Regazzoni, Dedè, and Quarteroni (2020) in the documentation and validation notes.

@ktbolt, would that naming approach work for you as well?

[michelebucelli]

Thanks @samibismar! Please see my comments below.

[ktbolt]

@samibismar Using RegazzoniActiveStress for the C++ class is fine and does follow the naming convention for these sort of things.

[michelebucelli]

@samibismar As we discussed, this model needs some kind of stabilization to enable force-strain-rate feedback. I think this can and should be addressed in a follow-up PR, but in the meantime I suggest the following:

1. to add a @todo comment in the class documentation to point this out, and
2. to add a boolean parameter (read from XML) to switch off force-strain-rate feedback (setting the parameter to true would replace $d\lambda / dt$ with 0 in the equations).

Having done this, I think this is good to merge!

[samibismar]

@michelebucelli I added the @todo for the follow-up stabilization work and the Disable_force_strain_rate_feedback parameter. When enabled, it replaces $\frac{d\lambda}{dt}$ with 0 in the XB update; the default is false so the existing behavior is preserved.

[javijv4]

Hey @samibismar, this looks good to me! Just left a small comment regarding the test. Not super important though.

[javijv4]

@michelebucelli, @samibismar, I don't think this is an issue for this PR, but I was playing around with the test and saw noticeably different results when using 3 or 4 processors when I let it run longer (below is timestep 50).
4 proc

3 proc

I checked using Nash Panfilov, and the same occurs. Then I tried changing GMRES to CG (just because, in my experience, CG works better for EP), and the differences became much smaller. Here are tables quantifying the differences.
=== GMRES: 3 vs 4 ranks at t=50 (max abs diff) ===

model	dVm	dTact	dCa	dDisp
Regazzoni	5.0823e+00	3.1990e-04	1.0133e-04	8.5070e-04
NashPanfilov	5.0823e+00	1.0954e-02	1.0133e-04	3.4885e-02

=== CG: 3 vs 4 ranks at t=50 (max abs diff) ===

model	dVm	dTact	dCa	dDisp
Regazzoni	2.4833e-06	1.5796e-11	8.5102e-11	1.8047e-11
NashPanfilov	2.4833e-06	1.5706e-09	8.5102e-11	1.9004e-09

I do not really know what is causing this. I can open an independent issue if you think that is a good idea.

[ktbolt]

@javijv4 Aren't the matrices produced by CEP symmetric and positive definite ? If so then you would use CG and not GMRES; same for the struct equation.

[michelebucelli]

Thanks @javijv4 for running these experiments!

I agree with @ktbolt that the CEP matrices are symmetric and positive definite, and so CG should be the best choice to solve them. It might be a good idea to switch the two electromechanics tests to use CG by default for electrophysiology.

However, for nonlinear hyperelasticity the matrix is not symmetric in general, I think (I'm not sure whether this holds fror the generalized $\alpha$ method, but the stiffness matrix at least is generally not symmetric). Therefore, struct should still require GMRES.

Having said that, GMRES should still work, and should work consistently when varying the number of cores. So I agree with @javijv4 that it is quite suspicious that the results change so much between 3 and 4 processes. I think this might be another hint that something, somewhere is very sensitive to the order of operations. I have encountered other clues to this in #584 (fast math optimization broke all the tests) and #577 (using x += ... in place of x = x + ... broke some tests). When prompted about this, Claude indicated that this might be related to the GMRES implementation, but I haven't investigated this claim so take it with a very big grain of salt.

@javijv4 if you want to open an issue about this I'll be happy to add more details about the similar problems I've found, and we can investigate further. I suspect this to be independent of the changes in this PR, so I suggest merging this PR regardless.

[javijv4]

Sounds good! I'll open a new issue.