Fix instability of force-length coupling in electromechanics, and uniform formulation of active stress tensor - #650
Conversation
![claude[bot]](https://avatars.githubusercontent.com/in/1236702?s=60&v=4){kind=small}
There was a problem hiding this comment.
Copilot review overview
🟡 Changes recommended
Unresolved critical restart-state and directional-validation issues, along with moderate solver coupling inconsistencies.
Get a fresh assessment by requesting another Copilot review.
Review effort: Lite
Findings: 2
Open (3)
What changed in this PR
This pull request standardizes active-stress formulation and improves electromechanical coupling stability through quadrature-point evaluation, stretch derivatives, and optional implicit state coupling.
Changes:
- Unifies active-stress evaluation across material models.
- Adds active-tension derivatives and quadrature-point evaluators.
- Updates solver integration, tangents, postprocessing, tests, and reference outputs.
| File | Reviewed change / status |
|---|---|
tests/unitTests/material_model_tests/test_material_common.h |
Updates material-test helpers for the new active-stress API. |
tests/cases/ustruct/LV_Guccione_active/result_001.vtu |
Regenerated regression reference output. |
tests/cases/struct/tensile_adventitia_Guccione_active/result_002.vtu |
Regenerated regression reference output. |
tests/cases/electromechanics/slab/result_Regazzoni_001.vtu |
Regenerated regression reference output. |
Code/Source/solver/ustruct.h |
Integrates active stress into incompressible assembly interfaces. |
Code/Source/solver/ustruct.cpp |
Passes evaluated active stress into incompressible assembly. |
Code/Source/solver/sv_struct.h |
Integrates active stress into structural assembly interfaces. |
Code/Source/solver/sv_struct.cpp |
Passes evaluated active stress into structural assembly. |
Code/Source/solver/post.cpp |
Uses quadrature-point active stress during postprocessing. |
Code/Source/solver/Parameters.h |
Declares implicit state-coupling configuration. |
Code/Source/solver/Parameters.cpp |
Implements implicit state-coupling configuration. |
Code/Source/solver/mat_models.h |
Updates constitutive active-stress APIs. |
Code/Source/solver/mat_models.cpp |
Applies the unified active stress and tangent. Critical (2 votes): sheet-normal stress with one fiber can access fl.col(1) without validating the number of fiber directions. Nit (1 vote): add nonzero active-tension finite-difference tangent coverage. |
Code/Source/solver/Integrator.h |
Declares stretch and active-state coupling helpers. |
Code/Source/solver/Integrator.cpp |
Integrates stretch computation and implicit state updates. Moderate (2 votes): stretch rate is computed before the displacement predictor, disabling Regazzoni force–velocity feedback. Moderate (1 vote): final Newton correction can leave mechanics and active state inconsistent. |
Code/Source/solver/fsi.cpp |
Updates FSI structural assembly interfaces. |
Code/Source/solver/ActiveStressUniformUnsteady.h |
Implements the zero direct stretch derivative. |
Code/Source/solver/ActiveStressUniformSteady.h |
Implements the zero direct stretch derivative. |
Code/Source/solver/ActiveStressRegazzoni.h |
Declares Regazzoni stretch-derivative support. |
Code/Source/solver/ActiveStressRegazzoni.cpp |
Implements Regazzoni stretch-derivative support. |
Code/Source/solver/ActiveStressNashPanfilov.h |
Implements the zero direct stretch derivative. |
Code/Source/solver/ActiveStress.h |
Defines active-tension data, evaluator, and lifecycle APIs. |
Code/Source/solver/ActiveStress.cpp |
Implements state gathering and tension evaluation. Critical (1 vote): restart and bin-to-VTK paths do not restore model state arrays, so freshly initialized state may be used instead of checkpoint state. |
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## main #650 +/- ##
==========================================
+ Coverage 73.08% 73.15% +0.06%
==========================================
Files 268 268
Lines 40270 40336 +66
Branches 6738 6752 +14
==========================================
+ Hits 29433 29507 +74
+ Misses 10594 10586 -8
Partials 243 243 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
…ance_active_stress, to simplify interface of update_active_stress
| // @todo[michelebucelli] Tfa is unused in this call, and it should | ||
| // probably be removed from the function signature. Active stress is | ||
| // added below in any case. | ||
| CANNModel.computeInvariantsAndDerivatives<nsd>(C, fl, nfd, J2d, J4d, Ci, Idm, Tfa, N1, psi, Inv, dInv, ddInv); |
There was a problem hiding this comment.
Is there anyone more familiar with this material model that can confirm removing Tfa from the arguments here is safe? The code would be equivalent, I'm not sure if there's some design reason or planned development that would justify keeping it.
There was a problem hiding this comment.
@msbazzi and @divya-adil should weigh in here. they are actively working on the CANN model implementation for more than one paper - so yes it is important and we need to keep it. They can explain the structure.
There was a problem hiding this comment.
@michelebucelli, you are correct that Tfa is not currently used in computeInvariantsAndDerivatives, so in the present CANN implementation, removing it from this particular call would not change the results.
The reason I would like to preserve the active-tension capability is because of where we are taking the CANN model next. We have been developing the CANN framework as a more flexible constitutive representation for vascular tissue, particularly for vein mechanics, rather than treating it only as another passive material model. The next step is to use the CANN-based constitutive description within a growth and remodeling framework.
My current plan is to build on the Fortran-based constrained-mixture/FSG implementation we discussed in #2, while making the changes needed to make it clearer and more general within svMultiPhysics. In that setting, the passive constitutive response and active smooth-muscle tension will ultimately need to coexist as part of the evolving wall mechanics.
That said, I’m also open to suggestions on the cleanest way to handle this. If Tfa does not need to be passed into computeInvariantsAndDerivatives yet, we can refactor it while keeping the active-tension pathway for the planned G&R implementation.
There was a problem hiding this comment.
@msbazzi thanks for explaining!
I think that, for the purpose of this PR, the question to address is whether the CANN-based stress tensor evaluation pipeline returns (now or in the planned future) the passive stress tensor only (option A) or the sum of passive and active stress tensor (option B).
Under option A, keeping Tfa as an input argument would mean that the passive stress can depend on the active tension, which is maybe counterintuitive but plausible in principle, and the code of this PR can be left as it is. Tfa can maybe be removed in the interest of code hygiene, but it is somewhat of a minor issue.
On the other hand, under option B the code of this PR is incorrect, because the lines after the switch add active stress regardless of the passive constitutive model, and the CANN-based model would end up having the active stress added to it twice (perhaps with two different formulations even).
I am not familiar with the CANN-based model, so there may be something I'm missing, but I think option A would be the cleanest, both conceptually and in the code, because it neatly follows the additive decomposition of stress which is behind the active stress formalism, and keeps passive and active constitutive models separate, which is a good step towards #178. However, if option B is necessary, then I'd have to rework this PR a little bit to accomodate for that.
I also note that the CANN model currently evaluates the stress as the derivative of an energy function (through CANNModel.evaluate), but active stress in general cannot be expressed in that way (I think), which makes me lean even more towards option A.
There was a problem hiding this comment.
@michelebucelli, currently (and I would say likely for the near future), CANN-based models are used only to describe passive stress, while active stress will come from a different formulation, given that I agree with you that A is the best and cleanest option.
There was a problem hiding this comment.
I see, thanks! I think the current code for this PR is mostly fine then.
I propose we remove the argument Tfa from this call, since it is currently unused and it might be a bit counterintuitive (as in, it might suggest option B to the reader, like it did for me), and we can bring it back if necessary in the future.
kko27
left a comment
kko27
left a comment
There was a problem hiding this comment.
@michelebucelli Thanks for making these changes! This definitely simplifies many of the function signatures. I mostly have clarification questions and comments related to future PRs.
| std::shared_ptr<IonicModel> ionic_model; | ||
| }; | ||
|
|
||
| /// @brief Cardiac electromechanics model type |
There was a problem hiding this comment.
@michelebucelli do you know why we have this cemModelType class? It looks like legacy code for stretch activated currents...
| /// @brief Compute active tension along fibers, sheets and sheet normals at | ||
| /// every mesh node. | ||
| // | ||
| void active_tension(const ComMod &com_mod, const int iEq, const mshType &lM, |
There was a problem hiding this comment.
@michelebucelli. Will the eventual goal be to make active_tension, fib_stretch, and fib_stretch_rate purely post-processing functions (just for writing to output files)? I see active_tension is already just for post-processing, and fib_stretch_rate is a bug if it is used in the integrator (as discussed previously).
| // Compute body forces | ||
| set_body_forces(); | ||
|
|
||
| // Implicit state coupling of the active stress: re-advance the state of the |
There was a problem hiding this comment.
@michelebucelli Just wanted to check if this is intended behavior. So if one of my domains has implicit_state_coupling set to true, then all of my domains get treated implicitly
| // If we need to compute fiber stretch, we iterate through all meshes, compute | ||
| // the stretch for each mesh, and then copy the mesh-local resulting vector | ||
| // into the global vector. | ||
| if (have_active_stress || need_fiber_stretch) { |
There was a problem hiding this comment.
We may at some point want to figure out how much computation time is spent from calling fib_stretch and fib_stretch_rate, especially for larger meshes. Maybe for implicit coupling, we don't always have to call this at every iteration
There was a problem hiding this comment.
@michelebucelli Does this mean that by default, all active stress models will treat the direct dependence on active tension implicitly? (i.e., there is no need to specify a xml parameter for implicit coupling, unless I want the states to implicitly depend on the fiber stretch).
| double Ja; | ||
| mat_models::compute_pk2cc(com_mod, cep_mod, dmn, F, nFn, fN, ya_g_f, ya_g_s, | ||
| ya_g_n, S, Dm, Ja); | ||
| mat_models::compute_pk2cc(com_mod, cep_mod, dmn, F, nFn, fN, Ta, S, Dm, Ja); |
There was a problem hiding this comment.
Could we move the active stress contribution outside mat_models::compute_pk2cc? I think this would be beneficial because it would isolate the use of ' mat_models ' to only passive constitutive relationships, and because then adding the active component would follow a similar logic to adding the viscous terms (S = S + Svis), which would make the different contributions to PK2 very clear in the struct_* functions.
That being said, I am guessing the reason not to do this is to avoid recomputing the structural tensors Hff, Hss, Hnn. If the loss in time is too important, maybe we should leave it as is and revisit this once the struct/ustruct physics are refactored.
|
@michelebucelli, another question: this PR adds derivatives to the active stress models. Should we add tests for these terms as well? I think we can do this in another PR, but it might be good to include this in future plans. |
Fixes #635; fixes #641.
Current situation
Release Notes
The PR addresses both issues above. Active tension, in particular, is now evaluated implicitly, and appropriate terms are added to the
(u)structmatrices to account for the coupling. A detailed list of changes follows.Modified
compute_pk2ccfrommat_models.cppso that all material models use the same formulation for the active stress tensor.Added the purely virtual function$T_\text{act}$ with respect to fiber stretch $\lambda$ . The function was implemented in concrete models, with all models returning
ActiveStress::compute_active_tension_derivative_localto compute the partial derivative of the active tension0.0except forActiveStressRegazzoni.Added the helper function$\lambda$ , bundled in an object of the new type
ActiveStress::compute_tension, to return active tensions along principal directions and their derivatives with respect toActiveStress::ActiveTension, introduced to be easily passed to downstream functions.The previous code computed active tension at integration points by
In this framework, computing the tangent contribution of active stress would have required computing the derivative of$T_\text{act}$ wrt displacement, which would have meant differentiating through the $L^2$ projection. This would be a nonlocal, and thus nontrivial, operation.
The new code evaluates$\lambda$ and $T_\text{act}$ directly at integration points. The nodal evaluation of $T_\text{act}$ is kept for the sake of output.
Added helper class
ActiveStress::evaluatorto compute active tension and its derivative at integration points.Modified
compute_pk2ccfrommat_models.cppto also account for the derivative of active tension when computing the tensorDm. This also required modifyingcc_to_voigt_eigento relax its symmetry assumption.By default, the active stress model state is still advanced explicitly (i.e. implicit coupling is only used to evaluate the active tension function). Optionally, it can be made implicit by setting the XML parameter
Active_stress/Implicit_state_couplingto true. This coupling does not contribute to the tangent, and is thus a fixed-point coupling between the two models. I implemented this experimentally, but opted to keep it in view of exploring its effect on force-velocity coupling stability (see Fiber stretch rate always evaluated to zero for active stress coupling #639).Refactoring of
Integratorto accomodate for the implicit coupling.Updated all downstream consumers of active stress to accomodate for changed interfaces.
Documentation
All modified/new classes or functions were given Doxygen documentation.
Testing
A few tests (
electromechanics/slab/Regazzoni,struct/tensile_adventitia_Guccione_active,ustruct/LV_Guccione_active) failed due to the change to the active stress formulation. I have regenerated their reference solution, as per the discussion in #635.Besides that, all automatic tests are passing.
Additional context
Code of Conduct & Contributing Guidelines