Problem
Fiber stretch and stretch rate are currently computed for all active stress models, regardless of whether each models use them (Lines 466-530 in Integrator.cpp). Computing these quantities at each time step is expensive and accounts for the majority (~75%) of the predictor step and about 5% of the total runtime for a 500k element cardiac structural mechanics simulation.
In Integrator::Predictor():
// Determine if we need to compute fiber stretch and stretch rate, by going
// through all domains of all equations until we find one for which active
// stress is enabled.
bool need_fiber_stretch = false;
bool need_fiber_stretch_rate = false;
int fiber_stretch_eq_index = -1;
for (int iEq = 0; iEq < com_mod.nEq; ++iEq) {
const auto &eq = com_mod.eq[iEq];
if (supports_active_stress(eq.phys)) {
fiber_stretch_eq_index = iEq;
for (const auto &dmn : eq.dmn) {
if (dmn.active_stress != nullptr) {
need_fiber_stretch = true;
need_fiber_stretch_rate = true;
}
}
}
Less importantly, but another area for improvement, even when fiber stretch and stretch rate are required, the fiber stretch rate is computed by calling fiber stretch at the current and old time step, resulting in 3 total fib_stretch function calls per time step.
In post.cpp, 818-821:
fib_stretch(com_mod, iEq, lM, solutions.current.get_displacement(), res);
fib_stretch(com_mod, iEq, lM, solutions.old.get_displacement(), lambda_old);
res = (res - lambda_old) / dt;
This could be fairly easily shortened to two fib_stretch calls by reusing the current fib_stretch value from the fiber stretch computation itself. More complicated, but number of fib_stretch calls could be reduced to one if we cache the fiber stretch value across time steps.
Solution
The simplest way to implement this would be to add two new flags for each fiber model stating whether stretch or stretch rate are needed.
Under the ActiveStress class in active_stress.h:
virtual bool needs_fiber_stretch() const { return true; }
virtual bool needs_fiber_stretch_rate() const { return true; }
Then in Integrator.cpp (line 475):
if (supports_active_stress(eq.phys)) {
fiber_stretch_eq_index = iEq;
for (const auto &dmn : eq.dmn) {
if (dmn.active_stress != nullptr) {
have_active_stress = true;
need_fiber_stretch |= dmn.active_stress->needs_fiber_stretch();
need_fiber_stretch_rate |=
dmn.active_stress->needs_fiber_stretch_rate();
}
}
}
@kko27 and I discussed this and wondered whether it would be better to eventually encapsulate the fiber kinematics logic and computation into the active stress class or elsewhere so the predictor method doesn't become cluttered. @zasexton Do you have thoughts? Are you planning to refactor the Integrator?
Additional context
No response
Code of Conduct
Problem
Fiber stretch and stretch rate are currently computed for all active stress models, regardless of whether each models use them (Lines 466-530 in Integrator.cpp). Computing these quantities at each time step is expensive and accounts for the majority (~75%) of the predictor step and about 5% of the total runtime for a 500k element cardiac structural mechanics simulation.
In Integrator::Predictor():
Less importantly, but another area for improvement, even when fiber stretch and stretch rate are required, the fiber stretch rate is computed by calling fiber stretch at the current and old time step, resulting in 3 total fib_stretch function calls per time step.
In post.cpp, 818-821:
This could be fairly easily shortened to two fib_stretch calls by reusing the current fib_stretch value from the fiber stretch computation itself. More complicated, but number of fib_stretch calls could be reduced to one if we cache the fiber stretch value across time steps.
Solution
The simplest way to implement this would be to add two new flags for each fiber model stating whether stretch or stretch rate are needed.
Under the ActiveStress class in active_stress.h:
Then in Integrator.cpp (line 475):
@kko27 and I discussed this and wondered whether it would be better to eventually encapsulate the fiber kinematics logic and computation into the active stress class or elsewhere so the predictor method doesn't become cluttered. @zasexton Do you have thoughts? Are you planning to refactor the Integrator?
Additional context
No response
Code of Conduct