Simbody is useful for internal coordinate and coarse grained molecule modeling, large scale mechanical models like skeletons, and anything else that can be modeled as bodies interconnected by joints, acted upon by forces, and restricted by constraints.
-
Tom Lintern
- Posts: 61
- Joined: Mon Dec 13, 2010 5:43 pm
Post
by Tom Lintern » Wed Aug 21, 2013 7:02 pm
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
-
Michael Sherman
- Posts: 806
- Joined: Fri Apr 01, 2005 6:05 pm
Post
by Michael Sherman » Wed Aug 21, 2013 7:26 pm
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
-
Tom Lintern
- Posts: 61
- Joined: Mon Dec 13, 2010 5:43 pm
Post
by Tom Lintern » Wed Aug 21, 2013 7:42 pm
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
-
Michael Sherman
- Posts: 806
- Joined: Fri Apr 01, 2005 6:05 pm
Post
by Michael Sherman » Wed Aug 21, 2013 10:16 pm
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
-
Tom Lintern
- Posts: 61
- Joined: Mon Dec 13, 2010 5:43 pm
Post
by Tom Lintern » Wed Aug 21, 2013 11:32 pm
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
-
Michael Sherman
- Posts: 806
- Joined: Fri Apr 01, 2005 6:05 pm
Post
by Michael Sherman » Thu Aug 22, 2013 8:37 am
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