Solar system simulation (n-body physics) with OpenCL

I would like to simulate a n-body solar system in Mathematica. I have found the following example: http://www.wolfram.com/mathematica/new-in-8/cuda-and-opencl-support/simulate-thousands-of-particles-in-real-time.html [1]

I have simplified the code from this example to create a minimal system with one "star" and one "planet". I plan to expand the solar system later with more bodies.

My code is:

Needs["OpenCLLink"]
OpenCLQ[]

"nbody_sim", {{"Float", _, "Input"}, {"Float", _, "Input"}, _Integer, "Float",
"Float", {"Local", "Float"}, {"Float", _, "Output"}, {"Float", _, "Output"}}, 256]


Note that I have changed all references to "Float4" from the example to "Float", because otherwise Mathematica would complain with errors.

numParticles = 2;
deltaT = 0.01;
epsSqrt = 1.0;
newPos = OpenCLMemoryAllocate["Float", {numParticles, 4}];
newVel = OpenCLMemoryAllocate["Float", {numParticles, 4}];


There are several possible issues with the last block of code. First, I don't know what is the epsSqrt parameter, I leave it the same as it is in the example code. Second, the bodies' initial positions are stored in the pos variable. I assume from reading the code of NBody.cl found here* that the specification for each body should be {positionX,positionY,positionZ,mass}, but I could be wrong. Third, the velocity vectors are stored in the vel variable. Here the first three components are {velocityX,velocityY,velocityZ}, but I do not know what is the fourth component.

If I simulate my simple system with

Graphics3D[{
AbsolutePointSize[20],
Point[Dynamic[Refresh[
NBody[pos, vel, numParticles, deltaT, epsSqrt, 256*4, newPos, newVel, numParticles];
NBody[newPos, newVel, numParticles, deltaT, epsSqrt, 256*4, pos, vel, numParticles];
Take[#, 3] & /@ OpenCLMemoryGet[pos],
UpdateInterval -> 0]]]
},
Boxed -> True, RotationAction -> "Clip", PlotRange -> All
]


the two bodies just stay where they are. Because they have no initial velocities, they should attract eachother (the heavier body should pull the lighter body towards itself basically), but this does not happen. How can I fix the code to simulate proper forces in a solar system?

*Edit: The file provided from the link is not the same that Mathematica uses. The location of the internal NBody.cl file is "/Applications/Mathematica.app/SystemFiles/Links/OpenCLLink/SupportFiles/NBody.cl". I am not sure if the source is available online.

Update

I have found this article [2] that explains how the NBody.cl program works in general, and an actual implementation. I have also read the Mathematica documentation page OpenCLFunctionLoad. From these two sites I obtained the following information:

• The proper data type in [1] should be "Float[4]" instead of "Float4" (and not "Float", which I used first).
• epsSqrt is "the minimum interaction distance (squared) between particles." If two particles are closer than this, the simulation will not compute the interaction between them.
• The fourth parameter in the pos variable is indeed the particle's mass.
• The fourth parameter in the vel variable is not used and can be set to any value.

I have gone through the Mathematica's internal NBody.cl file line-by-line comparing the code with the one provided in [2] and I have found no reason for the program to not work as it should. I have also written the code to solve the N-body system the "normal" way and it works as it should, however it's too slow to use for longer simulations with more bodies than just two in my minimal example here.

I believe I have done everything I can to identify the cause of my problem but I have not be able to find a solution yet.

-
I don't see any code that represents the force of gravity and its action. –  m_goldberg Dec 8 '13 at 16:18
Shoudln't NBody.cl simulate this? –  shrx Dec 8 '13 at 16:19
Maybe it does, but you don't show the code in that file. –  m_goldberg Dec 8 '13 at 16:20
It looks to me that the force equations are applied in the function IntegrateSystem`. Where do you call that? –  m_goldberg Dec 8 '13 at 16:41
Just copying the code and expecting it to work without any understanding usually fails. –  m_goldberg Dec 8 '13 at 16:58