where the collision of two star clusters is being simulated.

The overall scene is rather simple. There is a camera whose
sight is fixed at the origin, but can rotate and zoom freely.
The origin is augmented with axes which act as a reference point
during movement. The x-y plane is given a translucent
grid to give both scale to the simulation and a reference point
in the otherwise black space. The most interesting detail here
is the implementation of the bodies in the space.

At each time-step in the simulation the positions of the bodies
are passed to the visualization tool for rendering.
Each is rendered as single point but using a number of shaders.
There are three rendering passes of importance.
The first renders a small two-dimensional Gaussian texture at each point.
This gives the appearance of limb-darkening—where the
central part of a star appears brighter than the edges.
The next pass slightly blurs each texture for a subtle haze and twinkle
around each point. Finally, lighting bloom is calculated
where each point acts as a light source, creating the impression of
ambient light in the inter-particle space.
Lastly, it may appear that the color of each body is random,
but it is only to a point. In fact, we interpolate the apparent
color of main sequence stars from red through orange and white to blue.
In equal-mass simulations, a random color assignment in this interpolated range
gives visual interest and the ability to visually track a body over time.
Otherwise, one could assign color based on the mapping between mass and temperature.

§ EXPERIMENTAL RESULTS AND DISCUSSION

Throughout the last section, we
presented and developed many algorithms
for parallel $N$-body simulation.
These algorithms saw various techniques applied incrementally
to improve scalability and load-balancing.
This lead to one parallel version of the direct method
and seven versions of the Barnes-Hut method.
We summarize those methods here, and rename them as follows.

Algorithm 0: The parallelization of the direct $N^2$ method.
Each process is assigned an equal number of bodies,
and the force calculations and integration steps
are computed in parallel. This method passes each process's local set
of bodies around in a ring formation to compute forces at each time-step.

Algorithm 1: Barnes-Hut parallelization version 1. Bodies are gathered to the root process, a global octree is constructed serially, and then broadcast to all processes. Force computation (tree traversals) and integration steps are performed in parallel.

Algorithm 2: Parallel octree construction is added to Algorithm 1
where each process builds a local octree from its local bodies, and then a global tree
is constructing through pair-wise merges. The global octree is then broadcast. Force computation and integration steps continue to be computed in parallel.

Algorithm 3: Spatial decomposition is added to Algorithm 2. Each process is now responsible for some partition of the spatial domain at each time-step rather than a fixed set of local bodies.
This domain decomposition follows a Morton ordering on the bodies to assign an equal number of bodies,
but an unequal amount of space, to each process.  The re-partitioning occurs every time step.

Algorithm 4: Load-balancing via costzones is added to Algorithm 4. Bodies are constantly re-distributed such that each process computes roughly the same number of interactions per time-step, rather than having the same number of bodies. The spatial partitioning is still applied.

Algorithm 5: Collective gather and scatter operations are removed from Algorithm 4
with the addition of a distributed sort for the spatial decomposition and a greedy distributed re-balance
into costzones. This is the final variation of our costzones-based method.

Algorithm 6: The final collective operations are removed from Algorithm 5: the octree merge and broadcast of the global octree. The basic octree implementation is reformed as a hashed octree implementation to allow for unique and consistent indexing into all local octrees. The global octree
thus remains distributed. During tree traversal,
when a process requires non-local data, it is requested dynamically from the owning remote process.

Algorithm 7: The tree traversal of Algorithm 6 is improved by asynchronous communication and a latency-hiding tree traversal. Rather than waiting for requested remote data to be returned, those requested nodes are put on a so-called defer list and the algorithm proceeds in the meantime by traversing
other branches.

These eight algorithms have been implemented
in C and C++ using the OpenMPI v2.1.1 implementation <cit.>.
We also implement a OpenGL 3.3 implementation of
a 3D visualization tool for the simulations. Screenshots of which are shown
in Fig. <ref> and Fig. <ref>.
Our source code is freely available at

To test these algorithms, we consider a standard test
involving the interaction and eventual collision of two globular clusters.
Each spherical cluster contains equal-mass stars and follows a Plummer model <cit.>.
In this model, the mass density of a core of some Plummer radius is highly dense
and then falls off as $r^{\frac{-5}{2}}$ outside the core.
These clusters are initially separated from each other in each of the spatial dimensions
and then brought together under their mutual gravity.
In a simulation with $N$ bodies, each cluster is given $\nicefrac{N}{2}$.
Fig. <ref> shows how this simulation evolves over time for $N = 4000$.

All experiments in this section
simulate this interaction of two Plummer clusters,
varying only in $N$, the number of bodies in the system,
and $p$, the number of parallel processes executing the problem.
All trials use the same random seed when constructing the clusters.
All simulations are run in the standard units (see Section <ref>)
and with a multipole acceptance criterion $\theta = 0.5$,
a time-step of 0.01, and 500 total time-steps.
Our experiments were performed on a local area network cluster of 10 nodes,
each with two Intel Xeon X5650 processors (6 cores each; 12 cores per node),
a 12x4GB DDR3 main memory at 1.33 GHz, and running Ubuntu 18.04.4 and GCC 7.5.0.
For tests where $p \leq 10$, each process is mapped to a distinct
node of the cluster; where $p > 10$, processes are distributed evenly across
nodes and across sockets on each node.

$t = 0.0$

$t = 2.0$

$t = 4.0$

$t = 6.0$

The interaction and evolution of two Plummer model clusters, each with $2000$ bodies.

We begin our experimentation with validating the correctness of
our implementation.
Recall from Section <ref> that a stable
simulation will—up to numerical errors—conserve energy.
Table <ref> shows the percent change in total energy of the system
at the end of the 500 simulation steps.
These changes represent a change in energy of less than 0.0005 in standard units.
The original implementation of Barnes and Hut conserved
energy to about 1% <cit.>, thus these results are satisfactory.
Note also that for any one algorithm the change in energy was identical for
every number of processes $p$—even up to machine precision—suggesting
that the parallelization itself is correct and does not modify the operation
of the simulation itself.

$N$        10000        20000        40000        60000        80000

Algorithm     $p$

0     1       0.1391     0.1527     0.1538        –        –

100     0.1391     0.1527     0.1538      0.1415      0.1520

1     1       0.1324     0.1497     0.1483     0.1415     0.1520

100     0.1324     0.1497     0.1483     0.1415     0.1520

2     1       0.1324     0.1497     0.1483     0.1415     0.1520

100     0.1324     0.1497     0.1483     0.1415     0.1520

3     1       0.1324     0.1497     0.1483     0.1415     0.1520

100     0.1324     0.1497     0.1483     0.1415     0.1520

4     1       0.1324     0.1497     0.1483     0.1415     0.1520

100     0.1324     0.1497     0.1483     0.1415     0.1520

5     1       0.1324     0.1497     0.1483     0.1415     0.1520

100     0.1324     0.1497     0.1483     0.1415     0.1520

6     1       0.1324     0.1497     0.1483     0.1415     0.1520

100     0.1324     0.1497     0.1483     0.1415     0.1520

7     1       0.1324     0.1497     0.1483     0.1415     0.1520

100     0.1324     0.1497     0.1483     0.1415     0.1520

Absolute value of percent change in total energy for various experimental trials. Change in total energy was identical for all number of processes $p$ for any one Algorithm.
Change in energy for the direct method (Algorithm 0) can be attributed to the truncation error in integration and floating point error. Change in energy for all others also includes the added error of force approximation.

Next, we consider a large test suite of our eight algorithms.
For each algorithm we run the previously described
two cluster scenario for $N \in \{10000,20000,40000,60000,80000\}$
and $p \in \{1,2,4,8,10,20,40,80,100\}$.
Fig. <ref> summarizes these results.
As expected, Algorithm 0 runs much slower than the other algorithms
with its $\mathcal{O}(N^2)$ complexity. However, its
optimal parallelization and minimal communication overhead
lead to near linear parallel speedup.
For $N=10000$, algorithms achieve moderate speedup
for small $p$ but then communication overheads quickly dominate.
This can be seen as the up then down shape of parallel speedup curves for Algorithms 4–7
and the flattening of the curve for the other algorithms.
For $N=80000$, the achieved speedup is much better for the more complex algorithms,
with Algorithm 7 achieving the highest speedup, as hoped for.
Nonetheless, speedup still quickly flattens after $p=40$. This
again indicates communications are dominating and a larger $N$ is needed
to take full advantage of the scalability of the algorithms. The algorithms exhibit weak scaling
but not strong scaling.

First column: execution time vs simulation size (number of bodies $N$) for $p=10$ and $p=40$, respectively. Second column: parallel speedup vs number of processes, for $N=10000$ and $N=80000$, respectively. No data for Algorithm 0 for $N=80000$. Third column: the same as the second column, but focusing on the range $p\in[0,10]$. Parallel speedup plots include the ideal speedup as a dashed line. Algorithm 0 achieves near linear speedup but an exorbitant running time. The more simple algorithms achieve better execution time, but worse parallel speedup, indicating that larger values of $N$ are needed for other algorithms to perform favourably.

Throughout Section <ref>
we presented incremental changes as an attempt
to continually improve the parallelization
of the Barnes-Hut algorithms. Despite this,
the data shows some surprising results, listed below.

* The most simple Barnes-Hut algorithm, Algorithm 1,
achieves the best performance in terms of execution time.
Further, even for $N=80000$, it achieves the third-best parallel speedup.
The scheme is simple and does not involve much additional work for computing spatial decomposition or load-balancing, just a single gather of bodies and broadcast of the serially-constructed tree.
However, the slopes of the execution times for $p=40$ suggests this would no longer hold
for slightly larger values of $N$.

* Algorithm 7 achieves superlinear speedup for a low to moderate number of processes. While the work law tells us that $pT_p \geq T_1$ (the time to execute a program in parallel times the number of processes should always be more than the serial execution time), notice that Algorithm 7 in fact
can perform less overall work when more processes are added. This is due to each process
computing only a local tree and the global tree never fully constructed.

* Algorithm 4 achieves very poor performance. The change from Algorithm 3
to Algorithm 4 is the simple addition of costzones load-balancing.
Computing costzones is a simple $\mathcal{O}(N)$ operation which requires only a single
iteration over the list of work estimates. This alone cannot explain the difference.
Rather, this experimentation suggests that the performance of a regular gather and
scatter—where each process sends and receives the exact same size of data—is able to achieve
much better performance than a irregular gather and scatter—where each process sends
and receives a different amount of data. The costzones load-balancing requires the latter.
In implementation this causes the  and  calls
to be replaced with  and , respectively. Some slowdowns are expected when moving to the more generalized irregular routine <cit.>,
but this experimentation suggests the performance difference in OpenMPI is quite drastic. It could perhaps be
mitigated with a newer version of OpenMPI. Regardless, moving to a distributed sort and costzones in Algorithm 5 removes these irregular gather and scatters and again improves performance.

Top row: the parellel speedup vs number $p$ of processes for Algorithms 1, 5, 7.
Middle row: the execution time vs number $N$ of bodies in the system for Algorithms 1, 5, 7. Note log-log scaling.
Bottom row: the parallel cost (number $p$ of processes times execution time) per body vs number $N$ of bodies in the system for Barnes-Hut versions 1, 5, and 7. Note log-log scaling.

To examine more closely the performance characteristics of the different
parallelization algorithms, Fig. <ref> presents
data for Algorithm 1, Algorithm 5, and Algorithm 7.
Recall that Algorithm 5 is the final variation based on costzones
and Algorithm 7 is the final variation based on hashed octree.
Fig. <ref> presents 3 plots for each algorithm:
($i$) the parallel speedup achieved for various sizes of $N$;
($ii$) the execution time for various numbers of processes $p$
(as a log-log plot); and ($iii$) the parallel cost—$p$ times
the execution time in parallel $T_p$—scaled by $N$ and for
various number of processes (as a log-log plot).

First, the parallel speedups for Algorithm 1 show that
increasing speedup has reached a limit as $N$ increases.
Parallel speedup does not continue to increase with a
per-process increase in problem size, thus Algorithm 1
has poor scalability.
Algorithm 5 improves in its speedup
with increased work per process, but overall speedup is still lacking.
Communication costs begin to dominate.
Algorithm 7 achieves great and continual improvements in speedup
for increasing problem size, indicating its superior scalability.

Second, the execution time for each algorithm for
various process sizes is plotted against the value of $N$.
Where these log-log plots show nearly-parallel lines
indicates that adding more processes scales down work proportionately.
The separation between lines shows continual improvement in
speedup with adding more processes; for Algorithm 1
the lines for $p=40,80,100$ are not well-separated again indicating
its poor scalability. For Algorithm 5 and Algorithm 7
the curvature of the lines for high values of $p$ and low values of $N$
indicate how communication and parallelization overheads dominate the running time.
Where the lines have flattened for larger values of $N$ indicates that the algorithm
has reaches a sufficient problem size where the simulation work rather dominates.
Importantly, we see that for $p \geq 80$, neither algorithm has really reached its
full potential with these smaller data sets.

Third, we present the value $\nicefrac{pT_p}{N}$ vs $N$ for various values of $p$.
This parallel cost indicates how the additional work for parallelization
and communication overheads changes with adding processes and more work per process.
The distance between the $p=1$ line and any other line indicates the
parallel overhead. For larger values of $N$ these lines should
converge but remain parallel as communication becomes
a fixed proportion of the overall running time.
For Algorithm 1 all lines have reached
the parallel stage for low to moderate values of $N$,
indicating they have reached there limit in scalability.
For Algorithm 5, trends have begun to flatten for high values of $N$,
meanwhile the vertical distance between, say, $p=100$ and $p=1$ is still quite large.
This indicates that a large proportion of Algorithm 5
is spend on parallel overheads.
For large values of $N$ we see that the lines for Algorithm 7 converge,
and even intersect. Where any line is below that of $p=1$ indicates superlinear speedup.
For $p \geq 40$ we see that the lines are still trending downward,
indicting that Algorithm 7 has not reached its “large-$N$ limit”,
and further speedup can be achieved on any fixed number of processes
by increasing the simulation size.

§ CONCLUSION

Throughout this report we have examined
the necessary mathematical background and
the classical formulations of hierarchical
methods for gravitational $N$-body simulations.
We have collected and presented the work of many papers: <cit.>, <cit.>,
<cit.>,  <cit.>, <cit.>, <cit.>, and <cit.>;
applying their algorithmic techniques to progressively
improve a parallel and distributed Barnes-Hut algorithm.
The implementation of these algorithms in C/C++ using MPI
was tested on a small compute cluster with great results.

Overall, the experimental data suggests two things.
First, that the costzones method and octree merge <cit.> used
in Algorithm 3, Algorith 4, and Algorithm 5
incur a high degree of parallel overhead.
This method was first developed for a shared-memory system
where communication and inter-process synchronization is significantly cheaper.
On distributed systems this method is less practical, although much more
simple and easy to implement.
Second, the scalability of Algorithm 7 is highly encouraging.
While, for raw execution time, the simplest parallelization scheme of Algorithm 1
won out, this trend would not continue for larger simulation sizes
or for a larger number of processes. The hashed octree method <cit.>
can provide superlinear speedup for large simulations and
provides excellent scaling and speedup for increasing amounts of work per process.
Further refinement and experimentation is needed to test
the limits of this most complex but most promising algorithm.

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

Cambridge University Press

titleGravitational N-body simulations: tools and algorithms

titleA hierarchical O(N log N) force-calculation algorithm
Computational Astrophysics,Many Body Problem,Numerical Integration,Stellar Motions,Algorithms,Hierarchies,Physics (General)

journaltitleSIAM journal on scientific and statistical computing
titleA fast adaptive multipole algorithm for particle simulations

titleA graduate introduction to numerical methods

MIT Press

titleIntroduction to Algorithms, 3rd Edition

Budapest, Hungary

booktitleProceedings, 11th European PVM/MPI Users' Group Meeting
titleOpen MPI: Goals, Concept, and Design of a Next Generation MPI Implementation

journaltitleCommun. ACM
titleAn Effective Way to Represent Quadtrees

titleComputer Simulation Methods

Academic Press

journaltitleJournal of Computational Physics
titleA fast algorithm for particle simulations

booktitleThe use of supercomputers in stellar dynamics
titleStandardised units and time scales

journaltitleComputer Physics Communications
titleHierarchical N-body methods
N-Body Problem:Numerical Methods,Numerical Methods:N-Body Problem

crc Press

titleComputer simulation using particles

EDP Sciences

journaltitleAstronomy & Astrophysics
titleSEREN–a new SPH code for star and planet formation simulations-Algorithms and tests

titleThe Art of Computational Science

Institute of Physics Publishing

journaltitleThe Astrophysical Journal
titleBuilding a better leapfrog

titleStructured parallel programming: patterns for efficient computation

journaltitleInt. J. High Perform. Comput. Appl.
titleFast Parallel Tree Codes for Gravitational and Fluid Dynamical N-Body Problems

journaltitleACM Trans. Comput. Syst.
titleImplications of Hierarchical N-Body Methods for Multiprocessor Architectures

booktitleProceedings of the 1993 ACM/IEEE conference on Supercomputing
titleA parallel adaptive fast multipole method

journaltitleJournal of Parallel and Distributed Computing
titleLoad Balancing and Data locality in Adaptive Hierarchical N-Body Methods: Barnes-Hut, Fast Multipole, and Rasiosity

booktitle24th IEEE International Symposium on Parallel and Distributed Processing, IPDPS 2010, Proceedings
titleHighly scalable parallel sorting

World Scientific Publishing Company

titleBeautiful models: 70 years of exactly solved quantum many-body problems

journaltitleIEEE Trans. Parallel Distributed Syst.
titleSelf-Consistent MPI Performance Guidelines

booktitleProceedings of the 1993 ACM/IEEE Conference on Supercomputing
seriesSupercomputing ’93
titleA Parallel Hashed Oct-Tree N-Body Algorithm

journaltitleScientific Programming
title2HOT: an improved parallel hashed oct-tree n-body algorithm for cosmological simulation