Page 1 of 1
Cpodes vs RK integrator for numerically stiff problems
Posted: Wed Aug 21, 2013 7:02 pm
by tlin067
Hi,
I have been running some forward dynamics simulations using HuntCrossleyForce contact elements. This contact has resulted in my problem being numerically stiff and I have used the Cpodes integrator to solve the equations of motion.
Code: Select all
SimTK::CPodesIntegrator integrator(osimModel.getMultibodySystem(),CPodes::BDF,CPodes::Newton);
integrator.setAccuracy(0.000001);
I am struggling to get the solution to converge but the model appears to be close to convergence with an accuracy of 1e-6. This simulation takes ~5min (compared to ~3 hours when using the rk merson integrator). What I have noticed is that the rk methods seems to converge to a solution with a significantly less accuracy (~1e-3) and I am wondering why this is.
I have also been unable to increase the accuracy further for the CPodes integrator. When I increase the accuracy to 1e-7, I get the following error;
Code: Select all
[CPODES WARNING] CPode Internal t = 10.7813 and h = 1.88757e-16 are such that t + h = t on the next step. The solver will continue anyway.
I was hoping you could comment on these different integrator accuracies and any possible characteristics of my model which will contribute to these integration errors?
Thanks alot,
Tom
Re: Cpodes vs RK integrator for numerically stiff problems
Posted: Wed Aug 21, 2013 7:26 pm
by sherm
Hi, Tom.
I don't know what you mean by "converge" -- please clarify what you're trying to do.
Your system seems unusually stiff if CPodes is making that big of a difference. I suspect there is something wrong with the model. Some random possibilities: near-massless or -inertialess bodies, excessive stiffness or damping parameters for contact or joint stops, very large friction coefficients or small transition velocity, use of Ball joints with coordinate limits or actuators, some kind of bug in a custom component. We also found a bug recently in the HuntCrossleyForce where large dissipation values could cause the the contact force to go negative causing an instability. Smaller dissipation (like 0.1-ish) is safe.
I would be inclined to systematically study the model one element at a time to see what is causing the performance problems.
Regards,
Sherm
Re: Cpodes vs RK integrator for numerically stiff problems
Posted: Wed Aug 21, 2013 7:42 pm
by tlin067
Hi Sherm,
Sorry, I am using the solved kinematics to compare against some measurements in an objective function for an optimisation routine. By convergence I mean that the integrator accuracy is sufficient to ensure a smooth objective function. i.e. I keep increasing the accuracy until a consistent solution is reached.
The force parameters I am currently using are;
Code: Select all
<stiffness>30000</stiffness>
<dissipation>5</dissipation>
<static_friction>0.95</static_friction>
<dynamic_friction>0.7</dynamic_friction>
<viscous_friction>0</viscous_friction>
I think I got these numbers from an existing OpenSim example and my dissipation is clearly high (compared to 0.1 you suggested). My model is of the head and neck and the impact occurs during impact of the head with the torso. Unfortunately, the forces are dominated by the mass of the head meaning that each vertebrae has relatively little mass/inertia which could be contributing to the problem as I clearly see better performance when these values are increased. Do you think this is a deal breaker or is there some way to get around this?
Thanks,
Tom
Re: Cpodes vs RK integrator for numerically stiff problems
Posted: Wed Aug 21, 2013 10:16 pm
by sherm
Hi, Tom. I think it is OK if you can integrate it with CPodes. You could try lower dissipation and larger inertias for the vertebrae (assuming that doesn't matter for your investigation) to reduce stiffness but it will be hard to get RK to run as fast as CPodes here (check the contact transition velocity also -- larger is faster but slips more).
We made a change in the final release version of OpenSim 3.1 that might be helpful for you. The Optimizer was formerly computing the gradient without knowledge of the underlying accuracy to which the objective function is being computed. There is now an optional parameter where you can provide an accuracy estimate -- in your case that will cause the optimizer to take larger numerical differencing steps which calculating the gradient, which will produce a better derivative.
You can see the extra parameter and some documentation for it here:
Optimizer::useNumericalJacobian()
This might allow you to get away with a looser accuracy for the objective and still get reasonable convergence. It's still somewhat of a black art though to squeeze the best derivative out of an approximate objective; give it a try and let us know if it helps.
Regards,
Sherm
Re: Cpodes vs RK integrator for numerically stiff problems
Posted: Wed Aug 21, 2013 11:32 pm
by tlin067
Hi Sherm,
That's a great addition to the optimizer classes. I hit a bit of a bottleneck with that and resorted to wrapping the forward dynamics with the matlab optimization methods. One last question, do you know of anywhere I can get good estimates for the contact parameters? It would be great if there were some models where they have validated the choice of contact parameters. I understand you may not be the person to ask but thought I would try my luck anyway.
Thanks for all your help,
Tom
Re: Cpodes vs RK integrator for numerically stiff problems
Posted: Thu Aug 22, 2013 8:37 am
by sherm
Hi, Tom. Sorry, I don't know where you might find validated contact parameters, but I think you might be able to get some suggestions from the forum if you repost that question with a different topic title.
Regards,
Sherm