I have a 3-body simulation which must run millions of times.

As far as I know, the GPU shines when it gets to preform simple operations on huge matrices/arrays. Currently I'm debugging and running my simulation one at a time, in which case there is basically no need for the GPU, since 3 objects really isn't a lot. I'm also not using CPU parallelization for the same reason.

I'm trying to figure out the best way to use the GPU's strengths for calculating many 3-body simulations in concert. Figuring out how to use the GPU once all the acceleration vectors have been calculated seems pretty simple: with most algorithms I've seen, and specifically the ones that I use like RK4, accelerations are added to velocities (times dt or some other factor), and velocities to positions, so this just involves scalar multiplication and matrix addition, which I'd hope a GPU would be good at.

My main problem is figuring out how to use the GPU well for calculating the accelerations of many (mostly unrelated) bodies at once.

Some ideas I thought about:

I don't think using a sparse matrix would be very efficient, but it is true that:

$$\begin{bmatrix} 0 & _1{\vec{R}_{12}} & _1{\vec{R}_{13}} &0&0&0&0& \dots\\\ _1{\vec{R}_{21}} & 0 & _1{\vec{R}_{23}} &0&0&0&0& \dots \\\ _1{\vec{R}_{31}} & _1{\vec{R}_{32}} & 0 &0&0&0&0& \dots\\\ 0&0&0& 0 & _2{\vec{R}_{12}} & _2{\vec{R}_{13}} &0& \dots\\\ 0&0&0& _2{\vec{R}_{21}} & 0 & _2{\vec{R}_{23}} &0& \dots \\\ 0&0&0&_2{\vec{R}_{31}} & _2{\vec{R}_{32}} & 0 &0& \dots \\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots &\vdots & \ddots\end{bmatrix}\cdot \begin{pmatrix} _1 m_1\\\ _1 m_2\\\ _1 m_3\\\ _2 m_1\\\ _2 m_2\\\ _2 m_3\\\ \vdots\end{pmatrix} = \begin{pmatrix} _1 \vec{a}_1\\\ _1 \vec{a}_2\\\ _1 \vec{a}_3\\\ _2 \vec{a}_1\\\ _2 \vec{a}_2\\\ _2 \vec{a}_3\\\ \vdots\end{pmatrix}$$

Where the left index is the 3-body system number, and $\vec{R}_{ij} = \frac{\hat{r}_{ij}}{\lVert\vec{r}_{ij}\rVert ^2}$. Notably, this matrix is skew-symmetric.

An unrelated idea is to see that 3-body problem acceleration is somewhat reminiscent of a cross product: $$G\cdot \begin{pmatrix} m_1 \\\ m_2 \\\ m_3 \end{pmatrix} \times \begin{pmatrix} \vec{R}_{23} \\\ \vec{R}_{31} \\\ \vec{R}_{12} \end{pmatrix} = \begin{pmatrix} \vec{a}_1 \\\ \vec{a}_2 \\\ \vec{a}_3 \end{pmatrix}$$

As far as I know, GPUs need to be very good at cross products, to calculate normals to surfaces for graphics. This SO answer seems to be what I need in this case. I don't yet know enough about the subject to recognize if this is a viable option, though. Also, since $\vec{R}$s are vectors, I'd either need to cross product $\vec{m}$ with the $x$ and $y$ components separately, or "cheat" with complex numbers ($\hat{y}=i$), which I hear are possible in CUDA.

Another avenue is to somehow use matrices to calculate many squared distances (or normalized directions) between all the bodies at once, but I haven't been able to come up with how.

  • $\begingroup$ You don't need matrices, you can implement this matrix-free. Just feed an array with your data and have 1 thread work on one problem. $\endgroup$
    – lightxbulb
    Feb 6 at 8:57

1 Answer 1


CUDA has an example repository with codes for different use cases. There are samples for dealing with matrix operations like the approach you took above. One of the advanced use cases listed is the simulation of gravitational n-body systems. That might be a good start. Here are links:

List of code samples

n-Body simulation via CUDA

  • 1
    $\begingroup$ It worries me to know that my use case is considered advanced... Also, because their simulation is of N-bodies, perhaps some of their optimizations don't apply. This is still incredibly useful though, so thank you. $\endgroup$
    – Remeraze
    Jan 30 at 18:36
  • 1
    $\begingroup$ Maybe you could using batching instead (docs.nvidia.com/cuda/cublas/…)? $\endgroup$
    – IPribec
    Feb 5 at 19:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.