derivs
functionA lot of people had stability problems with their cloth solvers, and had to spend significant amounts of time tuning the parameters, etc. to get decent results.
After going through the projects and grading them, I can now say
with some confidence why this is the case. The problem for most
people comes down to the derivs
function.
Figure 1: The basic midpoint method for integrating a simple ODE.
The solvers in the GSL do some very fancy stuff using
embedded pairs to estimate error and adaptively change the timestep.
They also use higher-order methods than a simple Euler solver. As
a result, they will maintain stability in the vast majority of
cases (the exception being times when implicit solvers are strictly
necessary, which should not be the case in this assignment).
However, these methods are
entirely dependent on having a correctly-behaving
derivs
function. Specifically, the derivs
function needs to evaluate the state derivatives (f/m and v)
at the positions and velocities that are passed in.
This is easiest to see visually. Suppose we are trying to integrate
some simple ODE in 2 dimensions using the standard midpoint method.
This looks something like figure 1 at right, where we take two evaluations
of the derivs
function at different points in space and time,
and then use both of these to extrapolate a position at the end of the
timestep. The analytical solution (an approximately quadratic curve, in
this case) is shown in red.
Figure 2:The midpoint method using incorrect evaluations of the derivatives.
Now, suppose we always evaluate our derivatives using the positions
at the previous timestep. That is, in our derivs
function,
we evaluate at the passed in t value but using the x from the last time
we updated the positions/velocities in our system which was after the
previous timestep. In this case, the picture looks more like figure 2.
Notice that the result of this is that we are essentially
taking one big Euler step, despite our extra function evaluations.
Furthermore, we have crippled the solver's ability to estimate error,
since it cannot use the changes in derivatives for this purpose.
The key bit to see here is the necessity of using the
passed-in y vector in your derivs
function to compute the
forces in the system. That is, if you have code that looks like this,
int derivs (double t, const double y[], double f[], void *params) { ClothModel* model = reinterpret_cast< ClothModel *>(params); // Notice how the model here can't possibly be using the // positions from y to evaluate forces, since it is only // getting the out vector f. model->f(f); }then there is probably something wrong.
To see this in code, look at figure 3. You can also look at my full sample solution here. This code could probably be cleaned and tightened up a bit, I was in a bit of a hurry when I wrote it, but hopefully it gives you some idea of how things work. A sample excutable is available for both Windows and Linux. (note that you won't be able to compile the above file as it only includes the core cloth code and not the associated user interface code).
int derivs (double t, const double y[], double f[], void *params) { ClothModel* model = reinterpret_cast< ClothModel *>(params); // First part of the array is x, second part is v. // First, copy the velocities from the passed in state to // the state derivatives. std::copy( y+model->positionsStateSize(), y+model->stateSize(), f ); // Notice that the forces are evaluated at the passed // in positions and velocities! model->f( t, // current time y, // current positions y+model->positionsStateSize(), // current velocities f+model->positionsStateSize() ); // ougoing forces // Divide forces by mass (dv/dt = F/m) for( unsigned int i = 0; i < model->velocitiesStateSize(); ++i ) f[ model->positionsStateSize() + i ] *= model->invMass( i/3 ); return GSL_SUCCESS; }
Figure 3: the (annotated) derivs
function from the sample solution.
Questions? Contact Christopher Twigg.