'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?

Coffman

Could you also publish the Initial module, so this compiles? – Imponderabilia

github.com/idontgetoutmuch/Leapfrog/blob/master/Initial.hs – Coffman

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)

Chaussure

Note that NR codes are actually copyrighted, so you can't really use them in public codes. – Spotter

Thank you but this doesn't seem to address my question which is about the performance of the repa library for this particular problem. As an aside RK4 suffers from the same problem as the explicit and implicit Euler methods which is that energy is not conserved; I could have added "softening" to address the problem of planets getting too close but for the roughly circular orbits considered in the article this is unnecessary. – Coffman

Kyle, I am not implying that the code in Numerical Recipes be used without permission. Many examples shown in the book are based on non-copyrighted numerical computation methods, some even predating development of digital computers. It's a good place to start if one is unfamiliar with the various techniques. – Chaussure