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.41% 73.62% +0.20%
==========================================
Files 269 269
Lines 40337 40335 -2
Branches 6755 6744 -11
==========================================
+ Hits 29612 29695 +83
+ Misses 10482 10397 -85
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
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.
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...
There was a problem hiding this comment.
Originally, this class bundled the parameters for active stress and active strain.
The active stress configuration was moved to ActiveStress and its derived classes, so I removed it from here.
The active strain configuration is left here, but it seems it is inaccessible to the user, that is there's no way to set cemModelType::aStrain to true. The same goes for cemModelType::cpld.
There's a few places in the code where these two flags are read in if statements, and the code they guard is mostly unreachable because of this.
I think we can do one of two things:
- we keep everything as it currently is, waiting for the day when support for active strain will be restored;
- we remove the class and all the dead code its flags are guarding, and if that code is needed when active strain will be reimplemented we rely on
githistory to retrieve it.
I prefer option 2, because the resulting code is less misleading (e.g. does not suggest to the reader/new developer that active strain is supported when it is not).
There was a problem hiding this comment.
I also noticed that cemModelType::cpld also guards branches in cep_3d and cep_2d that modified the diffusion tensor to account for deformation. Those branches are currently dead code, since cpld can never be set to true, so I'm still in favor of dropping them.
Also, in light of this paper and my previous experience, those changes to the conduction tensor are not particularly significant and can generally be safely neglected (unlike e.g. stretch-activated currents). Still, should we want to bring them back, I think this should be done in a cleaner way by externalizing the evaluation of the deformation gradient (rather than implementing it again in cep_3d and cep_2d), so the current code would need to be reworked anyways.
There was a problem hiding this comment.
@michelebucelli Thanks! This makes sense. I'm in favor of removing the cem class and the evaluation of the deformation gradients. @aabrown100-git would you be okay with this? I think you added comments related to the active strain formulation in mat_models.cpp.
There was a problem hiding this comment.
I have tentatively removed the class in bf75a63, pending Aaron's opinion on this.
I restored SACs since that was a one-line change, but I removed the dead code for active strain and stretch-dependence of the conductivity tensor.
| 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.
There was a problem hiding this comment.
I agree with this. @michelebucelli and I discussed adding a general evaluator class for evaluating different quantities at integration points. I think it could be nice to be able to separate struct3d into kinematics (Nx, vx, F), passive stress, active stress, and viscosity, and prestress, since all of these contribute linearly to S.
As for recomputing Hff, Hss, and Hnn, is there any reason we need to store the fiber vectors at all or can we just store the structural tensors immediately after the fiber files are read? For interpolation of fiber direcitons, it makes sense to use the structure tensors anyway, and if we do need to access the vector itself for any reason, we can just compute the primary eigenvector of Hff, etc.
There was a problem hiding this comment.
I hadn't thought too much about this, but it's probably a good idea.
I think that long term the function compute_pk2cc should be retired and replaced by something like passive_model->compute_stress(F, Hff, Hss, Hnn, ...), where passive_model is an abstract constitutive model class (see also #178). This way, passive and active stress would be clearly separated, and all kinematics- and geometry-related quantities (F, Hff etc.) would be computed upstream (e.g. by struct_*) and then forwarded to the evaluators for the different contributions.
I think we can already start moving in that direction by splitting active stress out of this function. However, this relates to multiple issues/PRs currently open (at least this one and #640, and the evaluator class follow up I was discussing with @dseyler).
I think the evaluator class PR that will follow will be a good place to make this change.
|
@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. |
Good point, I will see if I can add something for this to the contraction model unit tests. |
…nvariantsAndDerivatives
The flags the class bundled were never set to true, so the class was guarding dead code everywhere. This dead code related to active strain (currently unsupported), changing conductivity in cep based on deformation, and stretch-activated currents. I have only restored support for stretch-activated currents (currently not covered by any automatic test).
…ins during nonlinear iterations
|
@javijv4 I've added a check in the unit test that compares the output of On top of this, the electromechanics tests also exercise |
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