Skip to content

Muscle FiberVelocityInfo and MuscleDynamicsInfo refactor #14

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 13 commits into
base: draft-merge-target
Choose a base branch
from

Conversation

pepbos
Copy link
Collaborator

@pepbos pepbos commented Dec 8, 2023

This PR combines Muscle::FiberVelocityInfo and Muscle::MuscleDynamicsInfo into a single cached variable, and removes member fields that are cheap to compute.

Considering that fetching memory is slower than multiplying a few numbers, the following fields are cheap to compute, and are removed:

  • FiberVelocityInfo::fiberVelocityAlongTendon,
  • FiberVelocityInfo::normFiberVelocity,
  • FiberVelocityInfo::pennationAngularVelocity,
  • FiberVelocityInfo::normTendonVelocity,
  • MuscleDynamicsInfo::fiberForceAlongTendon,
  • MuscleDynamicsInfo::normFiberForce,
  • MuscleDynamicsInfo::normTendonForce,
  • MuscleDynamicsInfo::fiberStiffnessAlongTendon,
  • MuscleDynamicsInfo::muscleStiffness,
  • MuscleDynamicsInfo::tendonPower,
  • MuscleDynamicsInfo::musclePower,

This leaves two fields in FiberVelocityInfo.
When computing the fiber velocity, with the exeption of Millard2012AccelerationMuscle, all muscles need to compute the forces as well. The forces can thus be seen as a byproduct of computing the velocities,
Furthermore, the FiberVelocityInfo and MuscleDynamicsInfo cached variables are invalidated at the same stage (Stage::Velocity).
Finally, it turns out that distributing the Muscle information over 3 cached variables is slower than over 2.

Since computing the forces is actually part of computing the muscle velocities, I argue for combining MuscleDynamicsInfo with FiberVelocityInfo.

Performance

Measuring the wall clock time for several simulations, shows around ~10% performance benefit:

main this PR
CMC Running Model 26.0 [secs] 21.4 [secs] (-22%)
Rajagopal forward integration 6.5 [secs] 5.9 [secs] (-10%)

Brief summary of changes

  • Moved all fields from MuscleDynamicsInfo to FiberVelocityInfo,
  • Removed Muscle::calcMuscleDynamicsInfo(), Muscle::getMuscleDynamicsInfo and the cached variable.
  • Added cached variable to Millard2012AccelerationMuscle for caching the muscle-specific curve values.
  • Added cached variable FiberVelocityInfoCache containing the MuscleLengthInfo and FiberVelocityInfo in a single struct
  • Moved fiberForceLengthMultiplier and fiberActiveForceLengthMultiplier from MuscleLengthInfo to FiberVelocityInfo (since they are used to compute forces/velocities and not lengths).
  • Added calc methods on FiberVelocityInfoCache variable for the removed fields mentioned above.

Some notes:

  • I think a more fitting name for the combined cached variable would be MuscleForceInfo, but adding that to this PR will create a very large diff, so that would be a follow up PR.
  • TODO: Thelen output changed, still debugging.
  • TODO: DeGrooteMuscle unit tests are failing on small numerical differences (<1e-13).

Testing I've completed

  • Rajagopal simulations showed no change in output.
  • Thelen simulation showing no change in output.
  • Unit tests passing

Looking for feedback on

@aseth1 would be great if you have time to take a look at Muscle.h. I moved the curve computations from MuscleLengthInfo to FiberVelocityInfo because they appeared a byproduct of computing the fiber force. But I might have misunderstood the invalidating behavior of CMC.
@adamkewley I left the calc methods in Muscle.h in the header e.g. Muscle.h::685

CHANGELOG.md (choose one)

  • no need to update because...
  • updated.

@pepbos pepbos requested review from adamkewley and aseth1 December 8, 2023 08:07
@@ -303,70 +287,61 @@ void DeGrooteFregly2016Muscle::calcFiberVelocityInfoHelper(
FiberVelocityInfo& fvi, const SimTK::Real& normTendonForce,
const SimTK::Real& normTendonForceDerivative) const {

// Multipliers.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What are Multipliers?

fvi.fiberActiveForceLengthMultiplier =
calcActiveForceLengthMultiplier(mli.normFiberLength);

double normFiberVelocity = SimTK::NaN;
if (isTendonDynamicsExplicit && !ignoreTendonCompliance) {
const auto& normFiberForce = normTendonForce / mli.cosPennationAngle;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Reference to temporary?

fvi.fiberVelocity =
fvi.fiberVelocityAlongTendon * mli.cosPennationAngle;
fvi.normFiberVelocity =
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The original code updates normFiberVelocity, whereas the new version NaNs it?

fvi.normTendonVelocity = fvi.tendonVelocity / get_tendon_slack_length();
const double tendonVelocity =
muscleTendonVelocity - fiberVelocityAlongTendon;
const double normTendonVelocity = tendonVelocity / get_tendon_slack_length();
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These local variables seem to have no purpose?

calcForceVelocityInverseCurve(fvi.fiberForceVelocityMultiplier);
fvi.fiberVelocity = fvi.normFiberVelocity *
fvi.fiberVelocity = normFiberVelocity *
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

These are the only side-effect-full parts of this block scope?

double fpe = fpeCurve.calcValue(mli.normFiberLength);

// Get multipliers specific to this muscle.
const ForceMultipliersCV multipliers = getForceMultipliers(s);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copies


void calcFiberVelocityInfo(const SimTK::State& s,
const MuscleLengthInfo& mli,
FiberVelocityInfo& fvi) const override;
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You plan on making this override-able?

fvi.fiberActiveForceLengthMultiplier =
get_active_force_length_curve().calcValue(arg);
fvi.fiberPassiveForceLengthMultiplier =
SimTK::clamp(0, get_passive_force_length_curve().calcValue(arg), 10);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is this arbitrarily clamped between 0 and 10?

@@ -458,8 +461,7 @@ void Thelen2003Muscle::calcFiberVelocityInfo(const SimTK::State& s,
//1. Get fiber/tendon kinematic information

//clamp activation to a legal range
double a = getActivationModel().clampActivation(getStateVariableValue(s,
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This one is clamped?

MuscleDynamicsInfo():
activation(SimTK::NaN),
FiberVelocityInfo():
fiberVelocity(SimTK::NaN),
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Use member initializers instead?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants