It's amazing how far our numerical methods have taken us. We're doing the most basic thing (forward Euler) which has the worst behavior in terms of numerical stability and yet for our purposes it has been fine, more or less. Let's look at the limitations and what our options are for the future.
First off, we need to distinguish between two problems. On the advection equation, where the result should be a smooth symmetric peak:
Issue 1
The result is distorted.
Changing the timestep doesn't help with this issue. With timestep = 0.0002, at t = 1048576 0.0002 = 209.7152:
With timestep = 0.0001, at t = 2097152 0.0001 = 209.7152:
Using double or Runge-Kutta 4 (see below) doesn't help either. This isn't a precision issue.
But using larger stencils does help with this issue. Larger stencils can have higher orders of accuracy. With our new stencil engine we can support more stencils pretty easily, and can make their size an option if we want to. It slows down the simulation of course.
Using a 4th-order stencil for the x-gradient instead of a 2nd-order one, timestep = 0.0002, at t = 1048576 * 0.0002 = 209.7152:
2nd-order stencil for x-gradient: [ -1, 0, 1 ] / 2
4th-order stencil for x-gradient: [ 1, -8, 0, 8, -1 ] / 12
Issue 2
With different settings we see a different problem:
The issue here is numerical stability. The errors soon explode.
Taking smaller timesteps helps. With timestep = 0.05 we can get to t = 7.0 until instability. With timestep = 0.025 we can get to t = 15.0. With timestep = 0.0125 we can get to t = 27.5.
Using double helps. With timestep = 0.05 we can get to t = 15.0 until instability.
Higher-order stencils don't help. In fact they make this issue worse. With 4th-order stencil for x-gradient, timestep = 0.05 we can only get to t = 4.0 before instability.
Using higher-order explicit methods helps. Runge-Kutta 4 is the most popular of these - see advection_RungeKutta4.vti for an implementation of explicit RK4. With timestep = 0.2 we can get further than forward-Euler does with timestep = 0.01:
RK4 with timestep = 0.2, after 136 steps (t = 0.2 * 136 = 27.2):
Forward-Euler with timestep = 0.01, after 2720 steps (t = 0.01 * 2720 = 27.2):
Using implicit methods would help. These would require a completely different implementation though, maybe not even OpenCL. And while they help with stability, they have other problems, and aren't a silver bullet.
Thoughts:
I'm coming to the conclusion that we need to do two things, for different reasons:
1) Make an easy-to-use option to switch to using higher-order stencils, for formula rules. The formula will remain unchanged and laplacian_a still means the same thing but would have higher accuracy. This would be useful for the researcher who is writing a paper and the reviewer asks them whether their results might not be the result of simulation inaccuracy. (Not an idle worry, see Ready's citations!) They can also reduce the timestep but as we've seen above, this doesn't guarantee that the result is correct.
2) Implement RK4 for formula rules. The motivation here is to find out if this allows us to run faster. RK4 is more expensive to compute but as long as it isn't 20 times more expensive then it is likely to be a win. I'm hoping that we can do this and still support all the formula rules. (We would need to forbid certain things, like access to get_global_id().)
It's amazing how far our numerical methods have taken us. We're doing the most basic thing (forward Euler) which has the worst behavior in terms of numerical stability and yet for our purposes it has been fine, more or less. Let's look at the limitations and what our options are for the future.
First off, we need to distinguish between two problems. On the advection equation, where the result should be a smooth symmetric peak:
Issue 1
Changing the timestep doesn't help with this issue. With timestep = 0.0002, at t = 1048576 0.0002 = 209.7152:
With timestep = 0.0001, at t = 2097152 0.0001 = 209.7152:
Using double or Runge-Kutta 4 (see below) doesn't help either. This isn't a precision issue.
But using larger stencils does help with this issue. Larger stencils can have higher orders of accuracy. With our new stencil engine we can support more stencils pretty easily, and can make their size an option if we want to. It slows down the simulation of course.
Using a 4th-order stencil for the x-gradient instead of a 2nd-order one, timestep = 0.0002, at t = 1048576 * 0.0002 = 209.7152:
2nd-order stencil for x-gradient: [ -1, 0, 1 ] / 2
4th-order stencil for x-gradient: [ 1, -8, 0, 8, -1 ] / 12
Issue 2
With different settings we see a different problem:
The issue here is numerical stability. The errors soon explode.
Taking smaller timesteps helps. With timestep = 0.05 we can get to t = 7.0 until instability. With timestep = 0.025 we can get to t = 15.0. With timestep = 0.0125 we can get to t = 27.5.
Using
doublehelps. With timestep = 0.05 we can get to t = 15.0 until instability.Higher-order stencils don't help. In fact they make this issue worse. With 4th-order stencil for x-gradient, timestep = 0.05 we can only get to t = 4.0 before instability.
Using higher-order explicit methods helps. Runge-Kutta 4 is the most popular of these - see advection_RungeKutta4.vti for an implementation of explicit RK4. With timestep = 0.2 we can get further than forward-Euler does with timestep = 0.01:
RK4 with timestep = 0.2, after 136 steps (t = 0.2 * 136 = 27.2):
Forward-Euler with timestep = 0.01, after 2720 steps (t = 0.01 * 2720 = 27.2):
Using implicit methods would help. These would require a completely different implementation though, maybe not even OpenCL. And while they help with stability, they have other problems, and aren't a silver bullet.
Thoughts:
I'm coming to the conclusion that we need to do two things, for different reasons:
1) Make an easy-to-use option to switch to using higher-order stencils, for formula rules. The formula will remain unchanged and laplacian_a still means the same thing but would have higher accuracy. This would be useful for the researcher who is writing a paper and the reviewer asks them whether their results might not be the result of simulation inaccuracy. (Not an idle worry, see Ready's citations!) They can also reduce the timestep but as we've seen above, this doesn't guarantee that the result is correct.
2) Implement RK4 for formula rules. The motivation here is to find out if this allows us to run faster. RK4 is more expensive to compute but as long as it isn't 20 times more expensive then it is likely to be a win. I'm hoping that we can do this and still support all the formula rules. (We would need to forbid certain things, like access to
get_global_id().)