5

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[]

srcf = FileNameJoin[{$OpenCLLinkPath, "SupportFiles", "NBody.cl"}]

NBody = OpenCLFunctionLoad[{srcf}, 
"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;
pos = OpenCLMemoryLoad[{{0.,0.,0.,100.},{1.,0.,0.,1.}}, "Float"];
vel = OpenCLMemoryLoad[{{0.,0.,0.,0.},{0.,0.,0.,0.}}, "Float"];
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.

shrx
  • 7,807
  • 2
  • 22
  • 55
  • I don't see any code that represents the force of gravity and its action. – m_goldberg Dec 08 '13 at 16:18
  • Shoudln't NBody.cl simulate this? – shrx Dec 08 '13 at 16:19
  • Maybe it does, but you don't show the code in that file. – m_goldberg Dec 08 '13 at 16:20
  • You can see the code here, I have provided this link in the question. – shrx Dec 08 '13 at 16:21
  • 1
    It looks to me that the force equations are applied in the function IntegrateSystem. Where do you call that? – m_goldberg Dec 08 '13 at 16:41
  • I don't know, I just copied the code from the example. But the only possiblility is that it is called within the NBody function. – shrx Dec 08 '13 at 16:45
  • 2
    Just copying the code and expecting it to work without any understanding usually fails. – m_goldberg Dec 08 '13 at 16:58
  • Yes, I agree. I am not familiar with CLink, but I searched for "n-body simulation in mathematica" and found that example. So I wanted to try the already existing code before doing everything from scratch. – shrx Dec 08 '13 at 17:00
  • How did you get srcf = FileNameJoin[{$OpenCLLinkPath, "SupportFiles", "NBody.cl"}] to work? It doesn't work for me. – m_goldberg Dec 08 '13 at 17:17
  • It just works after I execute Needs["OpenCLLink`"] (I have to omit the code format here due to code formatting issues with the backtick). Make sure OpenCLQ[] returns True. Also, maybe it's OS-dependent? I am on OS X 10.9. – shrx Dec 08 '13 at 17:20
  • let us continue this discussion in chat – m_goldberg Dec 08 '13 at 17:22
  • @m_goldberg I have done more research on this problem and I belive your previous remark about "just copying the code" does not apply anymore. I now understand what the code does but still can't identify the cause of the problem. – shrx Dec 08 '13 at 22:23
  • I had the same problem. But I found that this is probably because we changed the number of particles in the simulation. If you set numParticles = n*localSize, where n is an integer, the simulation works! localSize = 16 on my platform. – Xeelee Apr 23 '15 at 02:37

0 Answers0