# How to leverage the GPU for parallel 3-body problem computations

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.

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.

• 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. Feb 6 at 8:57