CMSoft

OpenCL ODE solver

OpenCL Runge-Kutta ODE solver

Get sample source code.
Includes:
- The example in this section, CLRungeKutta46;
- Particles position integration with sources that take turns attracting and repelling the particles, as in the video below:

1. Overview

The OpenCLTemplate.DifferentialEquations.floatODE46 and OpenCLTemplate.DifferentialEquations.doubleODE46 are programs created to make it possible to use GPU acceleration to compute solutions to ordinary differential equations using a fourth order ODE solver with sixth order error prediction - correction [2].

Differential equations are mathematical equations which describe unknown functions of one or several variables that relates the values of the function itself and its derivatives of various orders. Differential equations are crucial to all exact sciences, such as engineering, physics, chemistry and even economics.

Examples include calculating trajectories and collision of particles in game engines, electron-proton interactions, gravitational calculations, dynamic modeling of deformable bodies and many more.

Differential equations are present in many areas of science and technology: whenever a deterministic relation involving some continuously varying quantities (modeled by functions) and their rates of change in space and/or time (expressed as derivatives) is known, finding numerical values often involves numerical solution to systems of differential equations.

There are packages to use GPUs to compute solutions to other important problems, such as solving linear systems and computing FFT. This work tries to cover the lack of an easy-to-use ordinary differential equation solver for scientific applications and games.

The algorithm can solve ODEs using double-precision in systems that support it to accelerate scientific computing and it can solve ODEs using single-precision to compute visually correct real-time ordinary differential equation solving. It is suitable for solving differential equations that involve large amounts of functions.

The implemented function is an example of a series of masses connected by springs with the first mass connected to the ground. All masses are 1 and the spring constants are 0.2.

2. Using the ODE46 solver

The description below describes the floatODE46 ordinary differential equation system solver.  However, if you are using the double-precision version, just read double everywhere the explanation below says float and you should be fine.

2.1 Create a new solver

To use OpenCLTemplate ODE solver, you first need to create a floatODE46 object using:

floatrk46 = new OpenCLTemplate.DifferentialEquations.floatODE46(0, 0.1f, InitState, floatdydx);

The first argument is the initial value of the independent variable, usually zero (initial conditions referred to x = 0).

The second argument is the step size to be used. 0.1 should be enough for most applications because using 0.1 step size the error is proportional to 1E-6.

InitState is the initial conditions vector. It also informs floatODE46 the order of the system, N [1], [2].

floatdydx is a delegate that calculates the derivatives array dydx[N] as a function of the independent parameter x[0] and the system state y[N].

2.2 Creating the derivatives delegate

The derivative delegate gives the user versatility to create their own OpenCL routine to calculate complex derivatives or even simply read from GPU memory, process derivatives and write the answer back to the OpenCL solver (although this would be less effective). The preferred method looks like this:

private void floatdydx(CLCalc.Program.Variable x, CLCalc.Program.Variable y, CLCalc.Program.Variable dydx)

{
CLCalc.Program.Variable[] args = new CLCalc.Program.Variable[] { x, y, dydx };
kernelfloatdydx.Execute(args, n);
}

kernelfloatdydx should calculate all the components of the derivatives to allow the runge-kutta driver to procceed with integration. In this particular case, the kernel simulates a spring-mass system and the OpenCL code is:

__kernel void floatdydx(__global float * x,
__global float * y,
__global float * dydx)

{
int i = get_global_id(0);
int n = get_global_size(0);

//y[2*i] - position
//y[2*i+1] - velocity

dydx[2*i] = y[2*i + 1];

if (i == 0) dydx[1] = - 0.2 * y[0];
else dydx[2*i + 1] = - 0.2 * (y[2*i] - y[2*i - 2]);

if (i < n - 1) dydx[2*i + 1] += - 0.2 * (y[2*i] - y[2*i+2]);
}

Please refer to [1] and [2] if you need more information about transforming ordinary differential equations into a system of ODEs. As another example, if you want to integrate the Lorentz differential equations, you could use the following kernel (global_size = 3):

__kernel void floatDerivs( __global       float * x,
__global       float * y,
__global       float * dydx)
{
int i = get_global_id(0);

//Lorentz
float sigma = 10;
float rho = 28;
float beta = 2.66666666666667;
if (i==0) dydx[0] = sigma*(y[1]-y[0]);
if (i==1) dydx[1] = y[0]*(rho-y[2])-y[1];
if (i==2) dydx[2] = y[0]*y[1]-beta*y[2];
}

2.3 Integrating the system

Once the differential equations that rule the system are defined, you can:

- Set the current system state using floatrk46.ResetState(0, InitState);, where the first argument is the value of the independent variable and the second is a float array that describes the system state at that point;

- Integrate the system to a final value using floatrk46.Integrate(value);, where value is the new value of the independent argument;

- Retrieve the current system state (independent and dependent variables values) using properties floatrk46.State and floatrk46.IndepVar;

- Retrieve the current error estimate using floatrk46.AbsError.

The system will integrate to the final value provided that the final value is a multiple of the step size. The step size can be modified by using the floatrk46.SetStep(newStep); command. The final value may not bbe exactly the desired value due to round-off errors. Corrections to variable stepsize and round-off correction should be made manually because there are too many choices and the programmer should assess their impacts on integration performance.

3. Conclusion

This work makes available a package to solve both single-precision and double-precision (if hardware supports) ordinary differential equation systems, preferrably with very large number of states. This implementation is an important contribution because, unlike linear system solving and FFT packages, there is no such implementation using the processing available from GPU cards.

Single-precision (float) calculations should be enough to deploy acceptable physics effects to games whereas scientists will mostly be interested in the double-precision version of the OpenCL RK46 driver.

4. References

[1] PRESS, W. H., FLANNERY, B. P., TEUKOLSKY, S. A., VETTERLING, W. T., Numerical Recipes in Pascal, Cambridge University Press, 1989

[2] STOER, J., e BULIRSCH, R. Introduction to Numerical Analysis, New York: SpringerVerlag, 1993.

 OpenCL ODE solver