'Repa' performance for planetary simulation

Question

I have written a simulation of the outer planets of the solar system using the Euler symplectic method and implemented this a) using repa and b) using yarr.

yarr seems to perform about x30 quicker than repa.

Given this, I didn't even try to use parallelism. Is there any obvious performance problems in my repa code? The repository is at github. I can produce a cut-down repa-only version if this is helpful, but then you won't get the performance comparison against yarr.

Alternatively, how do I debug performance issues in repa?

Answer 1

Most Euler numeric integration methods suffer from cumulative round-off error that will eventually cause the simulation to "blow up". You may want to investigate advanced numerical integration methods, such as 4th-order Runge-Kutta or predictor-corrector.

Another place where n-body problem simulations become sticky is when two bodies get very close, such as a moon with a very eccentric orbit about its planet. If one uses fixed time increments for the simulation, the error during large changes of angular velocity can lead to division-by-zero errors or division by very small values that result in the simulation blowing up. Use of a variable delta-t that depends on the angular velocity can be beneficial.

These suggestions are based on running many such simulations as a project for an undergraduate physics course I took in 1973, while testing various numerical integration methods. Runge-Kutta and predictor corrector methods have been around since the dawn of digital computing and a number of books are available. See, e.g., Numerical Recipes: The Art of Scientific Computing by William H. Press, Brian P. Flannery, Saul A. Teukolsky and William T. Vetterling. (Cambridge University Press, 1989)