The task of N bodies or how to blow up a galaxy without leaving the kitchen

Not so long ago I read the science fiction novel “The Three-Body Problem” by Liu Qixin . In it, some aliens had a problem - they did not know how, with sufficient accuracy for them, to calculate the trajectory of their home planet. Unlike us, they lived in a three-star system, and the “weather” on the planet depended strongly on their mutual arrangement - from the incinerating heat to the icy cold. And I decided to check whether we can solve such problems.

Physics of the phenomenon

To understand the problem, it is necessary to deal with the physics of the phenomenon. In the framework of the classical theory, the attractive force of two bodies is determined by Newton's law :

$$

Where $$ - the position of bodies in space, $$ - masses of bodies, $$- gravitational constant.
In the system of$$ Bodies on each of them will be affected by the force of attraction from the rest, which is expressed by the equation:

$$

Using Newton’s second law, we write the acceleration for each particle:

$$

Recalling that acceleration is the second time derivative of the coordinate, we obtain a second-order partial differential equation, which must be solved to obtain the trajectory of each body:

$$

It is important to note here that the complexity of computing a function $$ is equal to $$ and increases significantly with increasing number of interacting bodies.

Maths

The first and simplest method for solving differential equations is the Euler method , which is designed to solve equations of the form:

$$

Upon transition to the discrete region, we obtain:

$$

Where $$ Is the integration step, and $$- the number of integration steps. Thus, if we need to calculate the position of the bodies at a time$$then we should do $$integration steps. Here the first problem is visible - if$$large, then we need to take a large number of integration steps.

To apply the Euler method to our problem, it should be reduced to a first-order system. To do this, we introduce an additional variable - particle velocity:

$$

The second problem in solving systems of differential equations is the accuracy of the solution and its control. Accuracy can be improved in two ways: by decreasing the integration step and choosing a method with a higher order of accuracy. Both methods lead to increased computational complexity, but in different ways. For example, you can use the classical fourth-order Runge-Kutta method ; it requires four function calculations$$ at every step, but has an order of accuracy $$ (for comparison, the Euler method has an order of accuracy $$ and requires one calculation $$) The accuracy of the solution can also be controlled in several ways: compare with the analytical solution, solve using different methods or with different steps and compare the results, control third-party parameters and limitations that the solution must comply with.
Also, each of these methods has its drawbacks. Analytical decisions may be absent, or, even in most cases, completely absent. For example, for our task$$ tel analytical solution is only for $$, but even this is enough to test the accuracy of the methods. Solving a problem by two methods or with different steps increases the computation time, but this approach can be applied to almost any task. Not every task has limitations, but for ours they have them: at each step of integration we can control the implementation of conservation laws . This approach also increases the complexity of the calculation, but there is plenty to choose from; calculating the sum of the momenta or angular momenta of all particles has complexity of the order$$, while calculating the total energy of a system has order complexity $$

Fig 1.  An example of an elliptical trajectory.

Consider the simple case of a satellite moving in an elliptical orbit around the Earth. As the satellite approaches the Earth, its speed increases, and when moving away from the Earth it decreases, accordingly, the possibility arises of decreasing the integration step in time in the green part of the orbit, and increasing it in red and blue without changing the accuracy of the solution. Let's try to compare in more detail.

Table 1.   Investigated methods for solving differential equations

Marked	Order	Description
adams	1-5	Adams-Bashfort Method
euler	1	Euler Method
rk4	4	Classic Runge-Kutta Method
rkck	5	Kash-Karp Method
rkdp	5	Dorman-Prince Method
rkdverk	6	Werner Method 1) p. 181
rkf	7	Felberg Method 1) p. 180
rkgl	6	Implicit Gauss-Legendre Method
rklc	4	Implicit Lobatto Method
trapeze	2	Trapezoid method

To select the best method for our task, we will compare several known methods. To do this, we simulate the collision of two systems of bodies$$and measure the relative change in the total energy , momentum and its moment at the end of the simulation (maximum simulation time$$) In this case, we will vary the step and parameters of the integration methods and measure the number of function calls$$, respectively, those methods that with a smaller number of calls will lead to less loss, we will consider more acceptable.

a)	b)
Figure 2.  Relative change in energy a), momentum b), at the end of system simulation$$ bodies by various methods depending on the number of function calculations $$ double precision

From the graphs of Figure 2 it can be seen that the best ratio of the amount of function calculation $$and the relative change in the energy of the system of bodies in the fifth-order Adams and Dorman-Prince methods. It is also seen that for all methods with an increase in the number of computations$$the relative change in the momentum of the system increases. For a relative change in energy, this is also noticeable, but only for a few methods that could reach the threshold$$. This effect can no longer be attributed to the error of methods, but to the error of calculations, and a further increase in accuracy is possible only together with an increase in the accuracy of calculations with floating point.

a)	b)
Fig. 3.  Relative change in energy a), momentum b), at the end of system simulation$$ bodies by various methods depending on the number of function calculations $$ with quadruple precision (__float128)

Figures 3a and 3b show that the use of calculations with quadruple accuracy can reduce the relative energy losses up to $$, but you need to understand that the computation time is increased by two orders of magnitude compared to double precision. As in the case of double-precision calculations, the best ratio of accuracy and the number of function calculations$$possess methods of fifth-order Adams and Dorman-Prince.

The Dorman-Prince and Werner methods belong to the class of nested methods and can simultaneously calculate two solutions with a high and low order of accuracy (for the Dorman-Prince method, orders 5 and 4, and for the Werner method, orders 6 and 5). If these two solutions are very different, then we can break the current integration step into smaller ones. That allows us to dynamically change the integration step and reduce it only in those areas where it is required.

Let us compare the fifth-order Dorman-Prince, Werner, and Adams methods in more detail, over a longer simulation interval of our system ($$)

a) The  relative change in energy, momentum and its moment in the modeling process

b)  Number of function calculations$$ on the simulation interval $$

Fig. 4.  Dependences of the relative change in the physical parameters of the system and the number of function calculations$$ from time to time modeling with the Dorman-Prince variable-step method $$

a) The  relative change in energy, momentum and its moment in the modeling process

b)  Number of function calculations$$ on the simulation interval $$

Fig. 5.  Dependences of the relative change in the physical parameters of the system and the number of function calculations$$ versus time modeling by the Werner method with a variable step $$

Fig. 6.  Relative change in energy, momentum and its moment in the process of modeling by the fifth-order Adams-Bashfort method with a step$$

Fig 7.  Dependences of the number of function calculations$$ for fifth-order Adams, Dorman-Prince and Werner methods from simulation time

It can be seen that at the initial stage of the evolution of our system ($$) all three methods show similar characteristics, but at later stages some events occur in the system, as a result of which errors in the main parameters of the system (total energy, momentum and its momentum) jump sharply. But the Dorman-Prince and Werner methods cope with these changes much better due to the ability to reduce the integration step in “complex” sections, as a result of which the number of function calculations increases$$as seen in Figures 4b and 5b, but the total number of calculations $$ nested methods are smaller than the Adams-Bashfort method, as can be seen in Figure 7.


The video demonstrates that up to the point in time $$the motion is relatively calm, and after a collision of the centers of the “galaxies” occurs, which leads to a sharp change in the trajectories and a strong increase in the velocities of some particles. Moreover, to maintain the accuracy of the solution, it is necessary to reduce the integration step. Nested methods can do this automatically; the graphs show that in some parts of the evolution of the system, the integration step was reduced by almost two orders of magnitude with$$ before $$. When using the Adams method and such a step on the entire interval of the evolution of the system, we would not have waited for a solution.

Total

For the solution, it is better to use nested methods that allow you to dynamically control the integration step, and reduce it only on "complex" sections of the trajectory.

Do not chase the highest order methods. Even when using the data type 'double', they do not reach their potential capabilities, using the data types with greater accuracy greatly increases the time it takes to solve the problem.

CPU implementation

Now that the choice of a method for solving equations is defined, let us try to figure out the calculation of the interaction force for each particle. We get a double cycle for all particles:


The gravitational forces for each body are calculated independently, and to use all the processor cores, it is enough to write the OpenMP directive before the first cycle:


Because Since each body interacts with each, to reduce the number of processor interactions with RAM and improve cache usage, we have the ability to load part of the data into the cache and use it repeatedly:


Further optimization consists in taking out the contents of the function of calculating the force in the main cycle and eliminating division and multiplication by body mass m1 [b1]. Besides the fact that we have reduced the calculations a little, the compiler will be able to apply vector instructions of the SSE and AVX processor on such an expanded cycle.


Table 2.   Dependence of the calculation time (in seconds) of the function$$ on the number of interacting bodies $$ for various CPU implementations

$$	2048	4096	8192	16384	32768
simple	0.0425	0.1651	0.6594	2.65	10.52
openmp	0.0078	0.0260	0.1079	0.417	1.655
openmp + block + optimization	0.0037	0.0128	0.0495	0.194	0.774

It is clearly seen that the version with OpenMP support is accelerated by six times, exactly in terms of the number of cores, and the optimized version is faster by a little more than two times. So, with optimization, you should not rely only on concurrency. Interestingly, in single-stream calculations, the processor worked at a frequency of 3.6 GHz, in the parallel version (openmp) it dropped the frequency to 3.4 GHz, and in the parallel and optimized (openmp + block + optimization) it dropped to 3.3 GHz, but this it did not stop her from working 13.6 times faster. It is also seen that an increase in the computation time with an increase in the size of the problem is quadratic, and a further increase$$ makes the task unsolvable in a reasonable time.

GPU implementation

But there is a desire to make calculations even faster. There are several directions available for acceleration: GPU computation, function approximation$$. First, for the GPU computing, I chose OpenCL technology. For more convenient work, the CLHPP library was used . The main advantage of OpenCL is that the code can be run on the processor and on the GPU, which simplifies writing and debugging, as well as expanding the list of hardware to run. The Oclgrind tool helps in debugging , which in runtime shows incorrect OpenCL API calls and memory access problems.

Opencl

To get started with OpenCL, you need to get a list of available platforms. The most common platforms are AMD, Intel, and NVidia.


Next, after choosing a platform, you need to select the computing device that this platform represents:


And at the end of the preparatory phase, it is necessary to create a context and queues within which memory will be allocated and calculations will be performed. For example, a context that combines all the computing devices of a selected platform is created as follows:


To download the source code to a computing device, it must be compiled, the cl :: Program class is intended for this.


To describe the parameters of a function (kernel) that is executed on a computing device, there is a cl template: make_kernel.


Further, everything is simple: we declare a variable with the type of the kernel, pass the compiled program and the name of the computational kernel into it, we can start the kernel almost like a normal function.


The computational core for OpenCL itself is very similar to the 'openmp + block + optimization' option for the CPU, only in contrast to the CPU version, the first cycle is controlled using OpenCL (the cycle range is determined by the global_range variable from the kernel launch code, and the current iteration number is available from the kernel using the function get_global_id (0)). First, part of the body data is loaded from local memory, then processed. Local memory is common to all threads in the group, so the download occurs once for the group, and is processed by each thread in the group and, since local memory is much faster than global, calculations are much faster.

Cuda

The implementation for the NVidia CUDA platform is a bit simpler than OpenCL, we don’t need to create the device context and manage the execution queue (at least until we want to make a multi-GPU implementation). As in the case of OpenCL, we need to allocate memory on the GPU, copy our data into it, and then we can start the computing core.

You can read more about working with CUDA here .


Unlike OpenCL, CUDA does not specify the full range of iterations (in the OpenCL implementation it is global_range), but sets the grid size and block sizes in the grid, accordingly, the calculation of the current body number changes slightly, otherwise the kernel is very similar to OpenCL, with the exception of other names synchronization and specifier functions for shared memory. Another useful distinguishing feature of CUDA is that we can specify the required size of shared memory when the kernel starts. As in the OpenCL implementation, at the beginning of each iteration block we copy part of the data into shared memory and then work with this memory from all the threads of the block.


Table 3.   The dependence of the calculation time (in seconds) function$$ on the number of interacting bodies $$ for various GPU implementations and better CPU implementations

$$	4096	8192	16384	32768	65536	131072
openmp + block + optimization	0.0128	0.0495	0.194	0.774	---	---
OpenCL + half NVidia K80	0.004	0.008	0.026	0.134	0.322	1.18
CUDA + half NVidia K80	0.004	0.008	0.0245	0.115	0.291	1.13

In general, the result is disappointing, only 5-6 times faster than the CPU implementation. Even if we conduct calculations on the entire K80, we will get acceleration up to 12 times, but since Since the complexity of the task is quadratic, then in a reasonable time we will be able to process not 32768 interacting bodies, but 131072, which is only 4 times more.

Function Approximation $$

If you look closely at the function, which is set by the attractive force of two bodies, it is clear that it decreases quadratically with distance. Therefore, we can accurately calculate the force of interaction between close bodies, and approximately between distant ones. One well-known approach
is the treecode algorithm proposed by D. Barnes and P. Hat. Octo- tree is being built in the simulated spacecontaining in its leaves the coordinates and masses of the modeled bodies. The parent nodes contain the center of mass, the total mass of the child nodes and the radius of the sphere described around the bodies of the child nodes. The root of the tree contains the center of mass of all bodies, their total mass and the radius of the sphere described around them. When calculating the interaction force, the distance to the root of the tree is first considered if the ratio of the distance to the node to its radius is greater than some constant$$, then the root is considered one body with coordinates equal to the coordinates of the center of mass of the bodies included in it, and a mass equal to the sum of the masses of the daughter nodes, if the ratio is less than or equal to $$, then the procedure is recursively repeated for each child node. If a leaf of a tree is reached, then the interaction force is considered in the usual way. Thus, if one body is far removed from the compact group of other bodies, this group is presented to it as one body, and the interaction force is calculated not up to each body, but only up to one body. Due to this, the complexity of the algorithm decreases with$$ before $$at the cost of some loss of accuracy.

In my approach, I decided not to use the octree tree , but the kd tree, because it is easier to use and has lower storage overhead compared to the octree.

Let's go back to the implementation on the CPU. A kd-tree node can be represented in the form of a class containing pointers to the left and right descendants and information about the coordinates and mass:


With this method of storing the tree, we have two possible options for traversing the tree: either use explicit recursion, or use the stack ourselves. I settled on the second option.


As in the case of the “exact” CPU implementation, the force calculation function is called for each body. The cycle across all bodies can be easily parallelized using OpenMP directives.

But neighboring loop iterations in this case will refer to completely different parts of the tree, which does not allow efficient use of the processor cache. To overcome this problem, all bodies need to be traversed not in the original order, but in the order in which the bodies are located in the leaves of the kd-tree, then neighboring iterations will occur for bodies that are close in space, and will traverse the tree along almost the same paths.


This implementation has another problem - there is no universal way to parallelize such a tree traversal. But we can completely change the way the tree is stored in memory - we can store all the nodes in one linear array and completely abandon the storage of pointers to descendants, by analogy with the construction of a binary heap . When starting the numbering of nodes with$$if the node is in cell number $$, then his left child is in the cell $$, right child in the cell $$, and the parent in the cell $$. The right node corresponding to the left with the number$$will have a number $$. The array itself will have a length$$, and all leaf nodes will be located in the last $$cells, moreover, closely spaced nodes will be located in close cells of the array. The tree traversal function will change a bit, but the framework remains the same:


But this is not all the possibilities that storing the nodes in the array opens up to us - we can refuse the stack when crawling. To do this, in the code branch in which we go to the children of the node, we add the function to calculate the next node ($$), and in the branch in which we calculate the interaction force, we add the calculation of the next node with the current subtree skipped ($$)

To skip the current subtree, we need to go down to the root (direction$$), while we are in the right descendant, as soon as we reach the left, we go to the corresponding right subtree (direction $$), if we get to the root, then the tree traversal is completed.



Figure 8.   Skipping a subtree$$.

To go to the next node, if possible, go to the left child (direction$$), and if there is no descendant, then go to the next node 'from below' using the function $$.



Figure 9.   Transitions to the next node$$.

It might seem that we replaced the recursion with a loop$$ in function $$, and such a replacement does not give anything, but let's see how to determine if a node with a number $$right descendant. This can be done simply by checking its odd number (the right child is in the cell$$), for this it is enough to calculate $$. Those. we divide the number$$two if the least significant bit is set to one. But this can be done without a loop, in many processors there is a find first set instruction that returns the position of the first set bit, using it we turn the loop into four processor instructions.


After that, we can exclude the stack from the tree traversal function and replace it with a pair $$, after that the function looks even simpler:


In total, we got six combinations of tree traversal and force calculation. Compare these approaches in terms of computation time and quality. We take as a measure of quality the relative change in the total energy of the system after 100 iterations. As a model system, we take two interacting “galaxies” consisting of$$bodies each.

Table 4.   Combinations of tree traversal method and force calculation

Tree traversal / force calculation type	Tree with stack	'Pile' with a stack	'Heap' without a stack
Iterations by body number	cycle + tree	cycle + heap	cycle + heapstackless
Leaf bypass	nestedtree + tree	nestedtree + heap	nestedtree + heapstackless

a)  Function calculation time dependence$$ from the ratio of the critical distance to the tree node to its radius ($$)	b)  The dependence of the relative change in energy in the system on the ratio of the critical distance to the tree node to its radius ($$)
Figure 10.  Function calculation results$$ in various ways on the CPU (number of bodies $$)

It can be seen that the implementation of 'nestedtree + tree' is hopelessly behind in speed, because it lacks concurrency. Implementations with the location of tree nodes in the array and indexing as in a binary heap are in the lead. The relative change in energy is negligible for all variants with$$. Also in fig. 10a shows that for ($$) all variants of function calculation $$ significantly overtake in speed the most optimized version of the exact calculation ('openmp + block + optimization'), with a further increase $$ implementations with a tree lose the exact version.

GPU tree traversal

I tried to bypass the tree on the GPU using both OpenCL technology and CUDA. The option of storing nodes in the form of a tree was immediately discarded, and only the options were left with storing the tree in an array with indexing as in a binary heap. In general, the implementation of the computing core is not much different from the CPU version.


In the first version, tree traversal began in the order of numbering the bodies in the original array, so neighboring threads bypassed completely different parts of the tree, which negatively affected the performance of the GPU memory cache. Therefore, in the second embodiment, a traversal was applied, starting from the top of the tree, in this case, neighboring threads begin to traverse the tree along the same path, because neighboring tree tops are nearby and in space. It is also important that we chose numbering in the array of tree nodes not from scratch, but from one, in this case the leaves of the tree are stored in the second half of the array, and with the number of bodies equal to the power of two, we will have equal access to memory by index tn1.


When traversing the order of numbering of tree nodes, we got a performance increase. But this option can also be improved. The global memory in which the tree nodes are currently located is optimized for shared access., i.e. threads of one group should read words located in one memory block. In our case, the tree traversal starts along the same paths, and we request the same data with all the threads of the group, but the further we go deeper into the tree, the paths of neighboring threads diverge more and more, and we have to request different data, which reduces the performance of the subsystem several times memory. But the nodes of each subtree belonging to the same level are located in relatively close memory cells. Those. when traversing the rest of the tree, the neighboring threads of the computing core do not access the same nodes of the tree, but closely located in memory. To optimize this memory access, texture memory can be used. But there is one snag. At the moment, textures do not support working with double precision data (we want to accurately calculate). But in CUDA there is a function__hiloint2double , which collects a double-precision number from two integers.


In this case, two implementations were made, in one each element of the tree (x, y, z, tree_crit_r2) was requested independently, and in the second implementation these requests were combined. The request for the mass of the node occurs much less frequently, only if the condition r2> tree_crit_r2 [curr] is fulfilled , so it makes no sense to combine this request with the others. Another useful feature of the CUDA framework is the ability to control the ratio of the sizes of the L1 cache and the size of shared memory ( cudaFuncSetCacheConfig ). In case of tree traversal, we do not use shared memory, so we can increase the L1 cache to the detriment of it.


An analysis of the program in the nvprof profiler showed that even when using texture memory to store a tree, there is still a very high load on global memory.

Indeed, in CUDA, all kernel memory that is addressed to 'calculated' addresses is stored in global memory, and accordingly, the stack that is needed to traverse the tree is located in global memory and eats away a significant part of the memory chip bandwidth, because the stack has each run thread, and there are a lot of threads.

But, fortunately, we already know how to traverse a tree without using a stack. Complementing the previous computational core with the functions of computing the next tree node, we get a new kernel, moreover, more compact.


The performance of the cores running on the GPU is highly dependent on the size of the blocks into which we divide the task. Depends on this size how many registers, local memory and other resources will be available for each computing thread. You also need to keep in mind that while waiting for memory access in one thread, another thread can perform calculations on the shader processor, thus, with a sufficient number of simultaneously executed threads on one processor, the memory access time will be hidden behind the calculations. Therefore, before comparing the performance of our cores, we need to calculate the optimal block size for each of them. Let's make a comparison on the half available to us from NVidia K80.

Table 5.   The dependence of the calculation time (in seconds) function$$ от размера блока для различных GPU реализаций при количестве тел $$ и $$

Размер блока/ядро	8	16	32	64	128	256	512	1024
opencl+dense	5.77	2.84	1.46	1.13	1.15	1.14	1.14	1.13
cuda+dense	5.44	2.55	1.27	0.96	0.97	0.99	0.99	-
opencl+heap+cycle	5.88	5.65	5.74	5.96	5.37	5.38	5.35	5.38
opencl+heap+nested	4.54	3.68	3.98	5.25	5.46	5.41	5.48	5.31
cuda+heap+nested	3.62	2.81	2.68	4.26	4.84	4.75	4.8	4.67
cuda+heap+nested+tex	2.6	1.51	0.912	0.7	1.85	1.75	1.69	1.61
cuda+heap+nested+tex+stackless	2.3	1.29	0.773	0.5	0.51	0.52	0.52	0.52

Таблица 6.  Зависимость времени вычисления (в секундах) функции $$ от размера блока для различных GPU реализаций при количестве тел $$ и $$

Размер блока/ядро	8	16	32	64	128	256	512	1024
opencl+dense	366	179	89.9	69.3	70.3	69.1	68.9	68.0
cuda+dense	346	162	79.6	60.8	60.8	60.7	59.6	-
opencl+heap+cycle	16.2	18.2	20.1	21.2	21.2	21.3	21.2	21.1
opencl+heap+nested	10.5	7.63	6.38	8.23	9.95	9.89	9.65	9.66
cuda+heap+nested	8.67	6.44	5.39	5.93	8.65	8.61	8.41	8.27
cuda+heap+nested+tex	6.38	3.57	2.13	1.44	3.56	3.46	3.30	3.29
cuda+heap+nested+tex+stackless	5.78	3.19	1.83	1.21	1.11	1.10	1.11	1.13

A difficult situation, but, unlike the CPU version of the tree traversal, it is clear that each optimization step brings tangible results. The implementation of 'opencl + heap + cycle' is almost 6 times slower than an exact solution with full function calculation$$. An implementation of 'opencl + heap + nested', in which a tree traversal starts from neighboring nodes, 1.4 times faster than the previous one, because Cache memory is used better. In the implementation of 'cuda + heap + nested', the L1 cache was increased to the detriment of shared memory, which increased the speed by 1.4 times, although it is possible that in the cuda implementation the compilation kernel is more optimally compiled (in the versions of 'opencl + dense' and 'cuda + dense 'cores are identical, and cuda's version performance is ~ 1.2 times higher). The most significant increase in computational speed (3.8 times) is achieved when the tree is located in the texture memory and the queries to the elements of the tree node are combined. Implementation with tree traversal without using the 'cuda + heap + nested + tex + stackless' stack is 1.4 times faster as well. in it, the entire bandwidth of the memory bus is used only to access data about tree nodes and is not spent on the stack. Thus, with$$ It was possible to achieve acceleration twice as compared with the full calculation of the function $$. But$$ an excessively large value of the parameter, on the graph of the relative change in energy in the system from the ratio of the critical distance to the tree node to its radius for the CPU implementation, it can be seen that lower values ​​can be used $$without a visible loss of accuracy. Let's try to vary$$ with the optimal block sizes that we determined in the previous step.

a) $$	b) $$
Fig. 11.  Dependence of the function calculation time.$$ from the ratio of the critical distance to the tree node to its radius ($$) for various GPU tree walk implementations

It can be seen that for small $$ all methods of function calculation $$go to close values, determined by the time of building the kd-tree and preparing the data for the GPU. Moreover, the time of building a tree makes a significant contribution to the total time to$$, then this time can be neglected. It is interesting to note that when$$ performance improves again upon reaching $$most likely this is due to the fact that all the threads of the GPU go around the same tree vertices at the same time, which improves the use of the L1 cache and completely eliminates the 'branch divergence' . It is also seen that the best performance for an implementation without a stack (cuda + heap + nested + tex + stackless), it outperforms the version with the stack at about$$times. Other implementations are several times slower. To consolidate the result, we will calculate the time on a GPU with a newer architecture.


When using the GeForce GTX 1080 Ti for calculation, the difference between tree walks with a stack and without a stack is up to two times, provided that we neglect the time it takes to build the tree. This fact pushes us to ensure that in the future porting to the GPU and building a tree. A comparison of tables 5-7 shows that there is no single optimal block size for a different number of bodies, and even more so for different GPU architectures, moreover, the calculation time difference can reach several times, even if you do not take boundary values ​​into account. Thus, before lengthy calculations, it makes sense to determine the optimal block size for each configuration.

The main thing that we have achieved is the ability to model pairwise interaction of a little over a million bodies ($$) for a reasonable time on one, not the newest GPU. On newer GPUs (Tesla V100), apparently, it will be possible to process already about two million interacting bodies in a reasonable time, because their performance is about four times higher than the half of the Tesla K80. Although this number is also incomparable with the number of stars even in dwarf galaxies, such as the Small Magellanic Cloud , but, in my opinion, is impressive.

Conclusion

The use of embedded methods for solving differential equations makes it possible to solve similar problems with a high level of accuracy at a relatively small time cost, and the application of algorithms for approximating the function of calculating the pairwise attraction forces allows us to process colossal volumes of interacting bodies. Thus, unlike the aliens from the “Three-body Tasks”, we are able to solve the problem$$ bodies, and for a small number of bodies and for entire star clusters, although at the cost of some loss of accuracy.

Visualization

For those who have read to the end, I will give a few more videos with visualization of the evolution of the systems of bodies.

Simulation of a collision of two galaxies. The total number of bodies is 60 thousand.


Modeling the evolution of a galaxy consisting of a million stars. In the center is a body with a mass of 99% of the total. Single particles are almost indistinguishable. Already more like a wave in a drop of liquid. Colorized according to speed module. Low speed - blue, medium - green, high - red. It is seen that in the center the speed is higher, and gradually decreases to the edges, and the lowest in the equatorial plane.


A more 'dynamic' simulation of the evolution of a million-star galaxy. In the center is a body with a mass of 10% of the total. The central body affects the rest significantly less. First, the “stars” fly apart, then gather back into several clusters, and at the end again form one large cluster.

In the process of modeling, about 5% of the “stars” left the initial region irrevocably.


At the 10th second, it very much resembles the existing cartwheel galaxy .

The project code can be found on the github .