Correct usage of the derivs function

A 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);
}

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;
}

Questions? Contact Christopher Twigg.