On Distributed Gravitational N-Body Simulations

Alexander Brandt
Department of Computer Science

University of Western Ontario

London, Canada

Email<EMAIL_ADDRESS>
The $N$-body problem is a classic problem involving
a system of $N$ discrete bodies mutually interacting
in a dynamical system.  At any moment in time there
are ${N(N-1)}/{2}$ such interactions occurring.
This scaling as $N^2$ leads to computational difficulties
where simulations range from tens of thousands of bodies
to many millions.
Approximation algorithms, such as the famous Barnes-Hut algorithm,
simplify the number of interactions to scale as $N\log{N}$.
Even still, this improvement in complexity is insufficient to achieve
the desired performance for very large simulations on
computing clusters with many nodes and many cores.
In this work we explore a variety of algorithmic techniques
for distributed and parallel variations on the Barnes-Hut algorithm
to improve parallelism and reduce inter-process communication requirements.
Our MPI implementation of distributed gravitational $N$-body
simulation is evaluated on a cluster of 10 nodes,
each with two 6-core CPUs, to test the effectiveness and
scalability of the aforementioned techniques.

Simulating the collision of two globular clusters, each consisting of 1000 stars.

§ INTRODUCTION

$N$-body problems involve a large collection of problems
wherein a set of $N$ bodies, particles, or points, interact under
mutually influential forces together in a dynamical system.
Often, the dynamical system is governed by a potential field,
for example, electrical, magnetic, or gravitational.
The result is a system of second order differential equations
where the acceleration of each body is computed from forces or potential
created by the other $N-1$ bodies. Then, velocities and positions are updated
from this acceleration.
Targeted applications include fluid dynamics,
molecular dynamics, quantum mechanics, astrophysics, and cosmology.
Historically, the $N$-body problem refers specifically to the problem where bodies are
governed by gravity and, moreover, will be the focus of this work.
Whether the physical scale of each body is very small or very large, these $N$-body problems
share two common issues: the lack of analytic solutions, and scalability.

As the number of bodies in a system increases, the number
of differential equations describing the dynamical system grows proportionally.
Due to the mutual interaction of those bodies,
the number of interactions, and the number of terms
in the differential equations, grows as $N^2$.
It is therefore prohibitively challenging to determine exact solutions
to these systems of equations as the value of $N$ increases.
For example, in quantum mechanics, and when restricting the number of
dimensions to only 2, exact solutions are only known up to the order of $N=10$ <cit.>.
For bodies moving under the force of gravity, general solutions are known only for $N=2$,
while for $N=3$ solutions exist for only very specific initial conditions <cit.>.
Simulations are therefore required to study systems with any practically meaningful number of bodies.
The efficiency of simulations with such a practical number of bodies is thus highly important.

In astronomy and gravitational $N$-body systems there are several different
simulation scales of practical importance:

* celestial mechanics, where a few dozen bodies—planets, stars—typically orbit a single body
of much higher mass than the others, e.g., a solar system;

* stellar dynamics, where a collection of roughly equal mass stars interact under their mutual gravity,
e.g., $10^4$–$10^6$ stars bound in a globular cluster;

* galaxy cluster evolution, where each body represents an entire galaxy and $10^2$–$10^3$ such galaxies are gravitationally-bound in a so-called galaxy cluster; and

* cosmological simulation, where $10^6$–$10^{11}$ bodies simulate the evolving structure of the universe.

These latter two, however, must also account for dark matter and relativistic behaviour
(see <cit.> and references therein).
The former two categories are dominated rather by Newtonian dynamics,
posing less of a mathematical challenge and more of a computational one.
Cosmological simulation also requires extreme-scale computing (petascale computing and beyond).
For these reasons, we focus on stellar dynamics throughout this report.

Where the number of bodies in a simulation increases from hundreds to thousands, or millions, or even billions,
the ability to simulate that number quickly becomes infeasible.
While, technically, simulating a system with $N^2$ interactions
is a polynomial-time algorithm, these $N^2$ interactions must
be computed at every time-step of the simulation.
The result is a $\mathcal{O}(TN^2)$ running time, for $T$
simulation steps. $T$ is relatively unbounded as scientists
wish to examine the long-term evolution of these systems.
As either $T$ or $N$ grows in orders of magnitude,
direct simulation of these $N^2$ interactions
becomes impossible.

Research into approximation algorithms
to reduce the running time complexity of such simulations,
while maintaining accurate results, has been ongoing for decades.
Several different methods appeared in the late 1980s
which brought the complexity of computing interactions down to $\mathcal{O}(NlogN)$.
The Particle-Mesh method <cit.> represents
particles[Particles, points, bodies are all terms used interchangeably to mean
the same thing: discrete positions in space. The differences only come from different application domains
and bear little difference in interpretation or understanding. Throughout this text we interchange the terms to keep the writing varied and interesting.] as a grid or mesh and then solves the potential on this density-mesh.
This method is well-suited where separations between particles are large, but
is insufficient where particle density is high.
Hierarchical methods, or treecodes, on the other hand, look to balance
approximating long-range interactions with using the direct method for close-range interactions.
The Fast Multipole Method <cit.> and
the Barnes-Hut algorithm <cit.>
are two examples of hierarchical methods. We explore both in more detail in
Section <ref>.

As computer hardware began to evolve in the 1990s,
so did support for parallel execution of $N$-body simulations
on multiprocessors and distributed systems.
Two different pioneering works appeared
simultaneously at Supercomputing '93.
The costzones approach by Singh et al. <cit.>, and the
hashed octree method of Warren and Salmon <cit.>.
Both works attempt to parallelize hierarchical methods through
adaptive domain decomposition, dynamic load-balancing across processors,
and parallel interaction computation.
The costzones approach was initially developed for a (distributed) shared memory system, while
the hashed octree method was developed explicitly for distributed systems.
On modern multi-core multiprocessor systems, the shared memory costzones approach
is simple and efficient on a single compute node.
Yet, as simulation sizes and durations grow, explicit distributed
computing is required to achieve higher performance without limiting simulations
to a single node. Nonetheless, the costzones approach can be adapted to
distributed computing, as we will explore.

Throughout this report we present a step-by-step guide to the implementation
of gravitational $N$-body simulation and, in particular, algorithmic
techniques which bolster a distributed and parallel
implementation. We collect and coherently present the work of many papers: <cit.>, <cit.>,
<cit.>,  <cit.>, <cit.>, <cit.>, and <cit.>;
with plenty of otherwise missing details coming from <cit.>.
The parallelizations and optimization techniques of the costzones and hashed octree methods
are presented incrementally to build up to a robust distributed algorithm (see Section <ref>).
The required details are made explicit through ample figures and algorithmic listings.
Further, our implementation, written in C/C++ and making using of the OpenMPI <cit.>
implementation of Message Passing Interface (MPI), is freely and openly available at
While no algorithm or technique here is cutting edge, nor at the level of modern super-computing,
it is nonetheless foundational and can be applied to small local area network clusters.

The remainder of this report is organized as follows.
We begin in Section <ref> by describing the mathematical model
underpinning gravitational simulations.
Section <ref> is further background describing the basis
of hierarchical methods through their use of spatial tree representations
(quadtrees, octrees) and force approximations.
Section <ref> explores the algorithmic techniques of
the costzones and hashed octree methods to parallelize treecodes.
These techniques are applied incrementally to produce continually more complex code
but with continually increasing scalability, resulting in
8 total algorithm configurations.
Finally, Section <ref>
evaluates the implementation of these algorithms
and presents our experimental data., while Section <ref>

§ GRAVITATIONAL $N$-BODY SIMULATION

This section reviews the mathematical background necessary to perform
an $N$-body simulation governed by the force of gravity.
While this information is not strictly necessary to understand the
distributed computing and algorithmic techniques discussed later in Section <ref>,
it does concretely define how “interactions” between bodies are computed.
Moreover, this section presents the mathematical reasoning, particularly gravitational potential,
which allows for the development of the approximation algorithms described in Section <ref>.

We begin in Section <ref> describing the direct interactions between
particles using Newton's law of universal gravitation. Then, Section <ref>
describes gravitational potential and the gravitational field. Through superposition of
gravitational potentials, and a multipole expansion of that superposition, the
acceleration due to the collective gravity of many bodies can be accurately approximated.
Section <ref> explains how this acceleration can be used to
update the dynamical system over time (i.e. perform integration).
Finally, Section <ref> explains how $N$-body simulations
can be scaled and standardized for consistent comparison across simulation scales.

This section is largely based on a comprehensive text by Sverre Aarseth;
see <cit.> and references therein for further details.

Throughout this section, scalars are given by Latin characters, and vectors are three-dimensional
and given by bold Latin characters.
The vector norm $||\cdot||$ is the Euclidean norm.

§.§ Gravitational Force

Newton's law of universal gravitation describes the attractive forces
between any two particles.
For two point masses with mass $m_i$ and $m_j$, respectively, separated
by a distance $r$, we know the magnitude of their mutual attractive gravitational force is given by

\begin{align}
F = \frac{G{m_i}{m_j}}{r^2}   \label{eqn:gforcesimple}
\end{align}

where $G$ is the gravitational constant.

For brevity, let us describe a particle whose mass is $m_i$ as the mass $m_i$.
Less mass $m_i$ be at position $\bm{r}_i$ and mass $m_j$ be at position $\bm{r}_j$.
Describing the force exerted on mass $m_i$ by mass $m_j$ can then be given by:

\begin{align}
\bm{F}_{ij} = \frac{G{m_i}{m_j}}{\lVert \bm{r}_i - \bm{r}_j \rVert^2} \cdot \frac{-(\bm{r}_i - \bm{r}_j)}{\lVert \bm{r}_i - \bm{r}_j \rVert}. \label{eqn:gforce}
\end{align}

Here, $\lVert \bm{r}_i - \bm{r}_j \rVert = r$ and $\nicefrac{-(\bm{r}_i - \bm{r}_j)}{\lVert \bm{r}_i - \bm{r}_j \rVert}$ is a unit vector pointing from mass $m_i$ to mass $m_j$. Using  Newton's second law,
the acceleration $\bm{a}_{ij}$ induced on mass $m_i$ by the force
$\bm{F}_{ij}$ is given by $\bm{F}_{ij} = m_i\bm{a}_{ij}$.

Using Newtonian dynamics
combined with this law of gravity yields the equations of motions for a system of bodies moving
under gravitational forces.
For a system of $N$ bodies there exists $3N$ second-order differential equations
describing the motion of the bodies. Using vector notation reduces the number of
equations to $N$ and improves readability.
Then, the position of the body with index $i$ evolves as:

\begin{align}
\ddot{\bm{r}}_i = \bm{a}_i &= \sum_{\substack{j=1\\j\neq i}}^{N} {\bm{a}_{ij}} \label{eqn:motion1} \\ &= \frac{1}{m_i}\sum_{\substack{j=1\\j\neq i}}^{N} {\bm{F}_{ij}} \\ &= -G \sum_{\substack{j=1\\j\neq i}}^{N} \frac{m_j(\bm{r}_i - \bm{r}_j)}{\lVert\bm{r}_i - \bm{r}_j\rVert^3}.  \label{eqn:motion}
\end{align}

Therefore, determining the interaction between two bodies $i$ and $j$
in an $N$-body simulation is simply the computation of
this $\bm{F}_{ij}$ or, equivalently, the acceleration $\bm{a}_{ij}$.
This defines a so-called particle-particle interaction.
Summing (<ref>) over all other bodies, as given by (<ref>), and dividing through by $m_i$,
yields the acceleration of body $i$, $\bm{a}_i$, as in (<ref>).

Due to the finite and numerical nature of the simulation,
it is useful to introduce a so-called
softening factor to a particle-particle interaction.
Including a small value $\epsilon$
in computing particle separations helps reduce the effects of close-encounters.
One simply replaces  $\lVert \bm{r}_i - \bm{r}_j \rVert$ with
$(\lVert \bm{r}_i - \bm{r}_j \rVert^2 + \epsilon^2)^{1/2}$.
This softening avoids numerical errors associated with singularities
and prevents odd behaviour arising from discretizing time steps during close encounters <cit.>.

However, we know that to compute $\bm{F}_{ij}$ for every pair of points
quickly becomes computationally infeasible.
One key observation, which is employed by hierarchical $N$-body methods,
is the approximation of the gravitational force exerted by many discrete bodies at once.
The fundamental observation is that for some subset of bodies $J$,
whose centre of mass is $\bm{r}_J$ and total mass is $m_J$, the force acting on particle $i$
from this ensemble of points is:

\begin{align}
\bm{F}_{iJ} = \sum_{j \in J} \frac{-G{m_i}{m_j}(\bm{r}_i - \bm{r}_j)}{\lVert \bm{r}_i - \bm{r}_j \rVert^3} \approx
\frac{-G{m_i}{m_J}(\bm{r}_i - \bm{r}_J)}{\lVert \bm{r}_i - \bm{r}_J \rVert^3} \label{eqn:approxforce}.
\end{align}

Recall computing the centre of mass is simply $\bm{r}_J = \sum_{j \in J}m_j\bm{r}_j$.
This approximation approaches equality as the separation distance $\lVert\bm{r}_i - \bm{r}_J\rVert$ grows larger.
However, this approximation can be further refined by using gravitational potential.

§.§ Gravitational Field, Gravitational Potential, and Multipole Expansion

Expanding from Newton's original idea of an attractive force between point masses,
a better model is that of the gravitational field.
The gravitational field of a point mass $M$ is a vector field describing the
force of gravity per unit mass felt by a body at any point in space
in the field. The gravitational field is typically denoted by $\bm{g}$ and,
for body with mass $m$ at position $\bm{r}$, is defined as:

\begin{align}
\bm{g}(\bm{r}) = \frac{\bm{F}(\bm{r})}{m} = -\frac{GM}{r^3}\bm{r}.
\end{align}

Notice that force per unit mass is precisely acceleration as from Newton's second law.
Hence, the acceleration due to gravity of mass $m$ by mass $M$ is precisely $\bm{g}$.

The gravitational potential $\Phi$ relates
energy, work, and force. It can be defined
as the gravitational potential energy per unit mass at a particular location.
Equivalently, it can be defined as the amount of
work per unit mass done against gravity
to move a mass from infinity to that location.
In the frame of reference of the point mass $M$
inducing a gravitational field (i.e. putting the point mass at the origin),
the gravitational potential at position $\bm{r}$
at a distance $r$ from the point mass is:

\begin{align}
\Phi(\bm{r}) = \frac{1}{m}\int_{\infty}^{r} \bm{F}\cdot d\bm{r} = \int_{\infty}^{r} \frac{GM}{r^2} dr = -\frac{GM}{r}.
\label{eqn:workpotential}
\end{align}

Gravitational potential is related to the gravitational field by the gradient.
Since gravity is an attractive force, gravitational potential
is negative while the gravitational field is positive.
The acceleration due to gravity of an object at position
$\bm{r}$ is then given as the negative gradient of the potential:

\begin{align}
\bm{a} = \bm{g}(\bm{r}) = -\nabla\Phi(\bm{r}) = -\frac{GM}{r^3}\bm{r}. \label{eqn:accelpot}
\end{align}

This implies that we can use gravitational potential
for the force calculations and interactions in our $N$-body system.
The useful consequence of this is that gravitational potentials
can be combined via superposition. Given a subset of points $J$, as before,
and a point $\bm{r}$ external to the sphere enclosing the points in $J$, the resulting potential is:

\begin{align}
\Phi(\bm{r}) &= \sum_{j \in J} -\frac{Gm_j}{\lVert\bm{r} - \bm{r}_j\rVert} \label{eqn:exactpotentialsum} \\
&\approx -\frac{Gm_J}{\lVert\bm{r} - \bm{r}_J\rVert}. \label{eqn:approxpotentialsum}
\end{align}

Much like forces, we can approximate the superposition of
potentials using their total mass $m_J$ and their centre of mass $\bm{r}_J$ as seen in (<ref>).
However, an even better approximation can be obtained using a multipole expansion.

For simplicity of notation, let us take $G=1$ henceforth.
A multipole expansion is a series expansion of a function
depending on angles. Often, the set of basis functions
for this series expansion are the spherical harmonics,
a set of orthogonal functions defined on the surface of a sphere.
In Cartesian coordinates these functions are the
Legendre polynomials. Following <cit.>,
let us denote the 	Legendre  polynomial of order (degree) $k$ as $P_k$.
The multipole expansion of (<ref>)
is then given by:

\begin{align}
\Phi(\bm{r}) = -\frac{1}{r}\sum_{j\in J}m_j \sum_{k=0}^n \left(\frac{r_j}{r}\right)^k P_k(\cos\theta_j), \label{eqn:multipole}
\end{align}

with $r_j = \lVert\bm{r}_j\rVert$ and $\cos\theta = {\bm{r}\cdot\bm{r}_j}/{rr_j}$.
The direct summation of (<ref>) is approached with
higher order $n$ in the series (<ref>).

In this series expansion the first few terms are called the monopole, dipole, and quadrupole terms.
Including even higher-ordered terms in the approximation yields diminishing
returns since their contribution to the approximation decreases with higher order
meanwhile their computational cost increases.
Often only up to quadrupole terms are included, yet more
recent works have included up to octopole terms <cit.>.

Using the compact notation of <cit.>,
the multipole expansion of $\Phi$ to order 2 is:

\begin{align}
\Phi_2(\bm{r}) = -\frac{M_0}{r} - \frac{D_ax_a}{r^3} - \frac{Q_{ab}x_ax_b}{2r^5}, \label{eqn:truncmultipole}
\end{align}

where $\bm{r} = (x_1, x_2, x_3)$ and the summation over $\{1,2,3\}$ for indices $a$ and $b$
is implied. There are therefore 3 dipole coefficients $D_a$ and 9 quadrupole coefficients $Q_{ab}$.
The multipole coefficients hide the summation over the
bodies in $J$ as:

\begin{equation*}
\begin{aligned}
M_0 &= \sum_{j \in J}m_j, \\
Q_{aa} &= \sum_{j \in J}m_j\left(3x_{a,j}^2 - r_j^2 \right),
\end{aligned}
\qquad\quad
\begin{aligned}
D_a &= \sum_{j \in J}m_jx_{a,j}, \\
Q_{ab} &= \sum_{j \in J}m_j\left(3x_{a,j}x_{b,j}\right), \ \ \text{if } (a \neq b),
\end{aligned}
\end{equation*}

where $x_{a,j}$ is the $a$th coordinate of $\bm{r}_j$.

Notice that the monopole term  is precisely the sum of the masses in $J$
as in our initial approximation.
Next, notice that taking the frame of reference to be the centre of mass of
the points in $J$ causes the dipole term to vanish ($D_a = 0$) by definition of $\bm{r}_J$.
Via the relation $\bm{g}(\bm{r}) = -\nabla\Phi(\bm{r})$,
the multipole expansion (<ref>) can be used in place of
the approximation (<ref>) for a set of masses $J$.
Therefore, we need only to add the quadrupole terms
to obtain an approximation which is better by two orders.

§.§ Energy of the System and Standard Units

In order to understand the accuracy of our simulation,
and in particular our approximation method,
a quantitative measure of accuracy or stability is very useful.
In a gravitational $N$-body simulation, this value is the total system's energy.
This is natural following the law of conservation of energy.
Since our simulation, by construction, simulates an isolated system, the total
energy should be conserved. The total energy of an $N$-body
system is given by the kinetic energy and the potential energy
of the $N$ particles. For masses $m_i$, positions $\bm{r}_i$, and velocities $\bm{v}_i$,
we have:

\begin{align}
E_{kin} = \sum_{i=1}^N \frac{1}{2}m_i\lVert \bm{v}_i\rVert^2, \qquad
E_{pot} = \sum_{i=1}^N \sum_{\substack{j=1\\j\neq i}}^N -\frac{Gm_im_j}{\lVert\bm{r}_i - \bm{r}_j\rVert}, \qquad E_{tot} = E_{kin} + E_{pot}.
\end{align}

By computing $E_{tot}$ at the beginning and the end of a simulation (or even throughout),
it is possible to obtain a measure of how well the simulation performed. A perfect simulation
would perfectly conserve energy. Of course this is not possible due to floating point errors,
truncation errors in the numerical integrator (see Section <ref>),
and the approximations introduced in the approximation algorithm itself.

While looking at conservation of energy for any one system is a good measure of its
stability, it is not completely adequate.
Different systems, with different initial conditions,
can have vastly different values for their total energy.
We would like a way to measure accuracy agnostic to particular initial conditions or
simulation parameters like $N$ or total mass.
The few practitioners of $N$-body simulation in the 1980s then decided on the so-called
standard units for gravitational $N$-body simulations,
see <cit.> and <cit.>.

The standard units define a scaling whereby all simulations
have an initial total energy $E_0 = -0.25$.
In the standard units, the total mass of the system is scaled to 1, the gravitational
constant is 1, and the system's coordinates are shifted to be centred at its centre of mass.
This latter fact makes the system's centre of mass the origin and the
net velocity $\vec{0}$.
To scale a system to standard units it is first shifted to its centre of mass frame
and then its kinetic energy and potential energy calculated.
From the virial theorem, a system in equilibrium
has kinetic energy equal to a negative half of its potential energy.
One can scale velocities and positions to independently and respectively
scale $E_{kin}$ and $E_{pot}$ to obtain $E_{kin} = \nicefrac{-1}{2}E_{pot}$.
Finally, the positions and velocities are scaled together
so that $E_0 = -0.25$;
see Algorithm 7.2 in <cit.>.

§.§ Simulating Dynamical Systems with Leapfrog Integration

Up to this point we have focused on how to calculate the interactions between bodies
in our simulation, and thus each body's acceleration.
Now we must use this acceleration to update the positions and velocities
of the bodies in our system so that it may evolve over time.

The simplest method is the basic Euler method.
Velocities are updated from accelerations, and positions are updated from velocities.
Let $\bm{r}_i^{(k)}$ and $\bm{v}_i^{(k)}$
be the position and velocity, respectively, of the $i$th particle at the $k$th time-step.
The forward Euler method uses values of derivatives at the current step
to approximate values at the next step:

\begin{align}
\bm{r}^{(k+1)} = \bm{r}^{(k)} + \bm{v}^{(k)}\Delta t, \\
\bm{v}^{(k+1)} = \bm{v}^{(k)} + \bm{a}^{(k)}\Delta t.
\end{align}

While simple, the forward Euler method is only a first-order method,
meaning that the (local) error introduced at each time-step
is proportional to $\mathcal{O}(\Delta t^2)$ <cit.>.
A higher-order method would see that the error reduces more quickly
with smaller time-steps.
For example, a common fourth-order method is the Hermite scheme <cit.>.
This scheme is a more complex generalization
of the classic and simpler scheme known as leapfrog, which is sufficient
in many cases.

The leapfrog method is as simple and as computationally expensive as the Euler method,
but is a second-order method, thus providing better accuracy.
Leapfrog specifically solves second-order differential equations of the form
$\ddot{\bm{r}} = \bm{a(\bm{r})}$ or, equivalently, $\dot{\bm{v}} =  \bm{a}(\bm{r}), \dot{\bm{r}} = \bm{v}$.
The method gets its name from the way that the two differential equations
are solved at interleaved time-steps, so that each equation jumps over the other in time.
The classic formulation of leapfrog shows this jumping structure:

\begin{align}
\bm{v}^{(k + \frac{1}{2})} &= \bm{v}^{(k - \frac{1}{2})} + \bm{a}^{(k)}\Delta t, \label{eqn:classicleapv} \\
\bm{r}^{(k + 1)} &= \bm{r}^{(k)} + \bm{v}^{(k + \frac{1}{2})}\Delta t. \label{eqn:classicleapr}
\end{align}

One of the key advantages to this scheme is its simplicity and time-symmetry.
That is, it provides the same results when run forward in time as
when run backward in time. In the context of gravitational simulations,
this time-symmetry prevents systematic build-up in error in the total energy of
the system over time <cit.>. Why this is a favourable condition
is discussed next in Section <ref>.

With a simple rearrangement and translation of time-steps, these equations
can be written at integer time-steps.
In particular, we can arrive at the
“Kick-Drift-Kick” (KDK) version[The KDK leapfrog is also known as the velocity Verlet method.] of leapfrog:

\begin{align}
\bm{v}^{(k+\frac{1}{2})} &= \bm{v}^{(k)} + \bm{a}^{(k)}\frac{\Delta t}{2}, \label{eqn:leapkick1}\\
\bm{r}^{(k+1)} &= \bm{r}^{(k)} + \bm{v}^{(k+\frac{1}{2})}\Delta t, \label{eqn:leapdrift}\\
\bm{v}^{(k+1)} &= \bm{v}^{(k + \frac{1}{2})} + \bm{a}^{(k+1)}\frac{\Delta t}{2}. \label{eqn:leapkick2}
\end{align}

The idea behind KDK leapfrog is to use acceleration from the previous step ($k$) to update velocity
a little bit, to time-step $k + \nicefrac{1}{2}$; this is the first kick.
Then, position is allowed to update by drifting along the trajectory induced by the updated velocity.
Then, a second kick occurs to update the velocity based on the acceleration at the current time step ($k+1$).

This variation of leapfrog allows for variable time-steps. That is,
changing the value of $\Delta t$ over the course of a simulation to give
more precision during volatile points of the simulation (i.e. close encounters).
While the KDK leapfrog is employed in our implementation,
we currently only use a fixed time-step.

§ HIERARCHICAL $N$-BODY METHODS

In $N$-body simulations, the dynamical system describing the
bodies includes $N$ differential equations.
By the mutual interaction of these $N$ bodies,
the differential equation describing each body includes $N-1$ terms
for the $N-1$ interactions of that body with all others.
The details of this dynamical system were described in Section <ref>.

From the tools established in the previous section—force or potential calculation,
acceleration calculation, and discretized updates to velocity and position
via leapfrog integration—we have everything required to perform an $N$-body
Directly computing the $N(N-1)$ interaction terms (or $\nicefrac{1}{2}N(N-1)$,
if taking advantage of symmetries), of the dynamical is infeasible for large-scale simulations.
Nonetheless, Algorithm <ref> depicts this simple scheme.
For each discrete time step in the simulation, the acceleration
is computed for each body (Lines 2–5) and then each body's position and velocity are updated (Lines 6–7).

Direct $N$-Body Simulation($P$, $\Delta t$, $t_{end}$)The general scheme loops over small time intervals, first computing the forces on and the acceleration of each body for that interval, then updating each body's position. $P$ is the list of bodies and $t_{end}$ and $\Delta t$ are simulation parameters for the total simulation time and time-step, respectively.
[1]
$t$ from 0 to $t_{end}$ by $\Delta t$
each particle $p$ in $P$ Force Calculation
Set the acceleration of $p$ to $\vec{0}$.
each particle $q$ in $P \setminus \{p\}$

Add the force contribution of $q$ acting on $p$ to the acceleration of $p$.
each particle $p$ in $P$ Position Update
Update position and velocity of $p$ using its current acceleration.

Computing these direct particle-particle interactions scales as $N^2$,
and is computationally infeasible in practice.
Much research has been dedicated to reducing the number
of interactions computed while maintaining accuracy in the simulation.
As we have seen in Sections <ref> and <ref>,
a collection of bodies, or a mass distribution, can be approximated
by its total mass and centre of mass. The multipole expansion
can also be used for even better approximation.
Moreover, recall that the approximation (<ref>) in Section <ref>
approaches equality with larger distances from this centre of mass.
Hierarchical methods (treecodes) take advantage of this approximation
to reduce the number of interaction terms computed in the force calculations of each time step
in the simulation.
These methods only modify the force calculation loop of the simulation,
leaving the position update, i.e. numerical integration, unchanged.

In the late 1980s two such hierarchical methods arose to reduce the computational efforts
required in the force calculation step of $N$-body simulations.
The first method to appear was by Barnes and Hut <cit.>
and eventually named the Barnes-Hut algorithm.
The second method of Greengard and Rokhlin is the Fast Multipole Method (FMM) <cit.>.
Both of these approximation algorithms—approximate in the way force acting on a body is
calculated, not by numerical or truncation errors—are based on a hierarchical
tree representation of the geometric space containing the bodies; hence, treecodes.
In two-dimensions these methods employ quadtrees, while in three-dimensions octrees are used.

The key idea in treecodes is to describe an ensemble of bodies, a mass distribution,
by a single representative pseudo-body in place of its
constituents. Then, a far-enough away test body can approximate its interactions
with each of the constituent bodies by a single interaction with the pseudo-body.
These are so-called long-range interactions.
Applying this technique hierarchically, the further away a target point is from
a group of points, the more of those points can be included in an ensemble
to further reduce the number of interactions to be computed for the target point.
Whenever two points are considered to be too close, and thus the approximation too inaccurate,
computations fall back to the direct method of computing particle-particle interactions.
Hierarchical methods thus drastically reduce the total number of interactions needed to be computed,
down to $\mathcal{O}(NlogN)$.

The Barnes-Hut and FMM methods are both hierarchical methods.
They differ in two places: ($i$) how an interaction is determined to be “long-range”,
and ($ii$) how long-range interactions are specifically computed.
However, let us begin with the commonalities.

Both methods begin by constructing a hierarchical tree representing the entire geometric space
which bounds the particles being simulated.
The root of this tree is precisely the bounding box of all the particles.
The root is then split in half in each of its dimensions. In two dimensions,
a quadtree has each node split into 4 children, the area of each being $\nicefrac{1}{4}$
of its parent. In three dimensions, an octree splits each of its nodes into 8 children,
the volume of each being $\nicefrac{1}{8}$ of its parent.
This splitting continues recursively until each leaf node contains exactly one particle.
A quadtree is shown in Fig. <ref> while an octree is shown in Fig. <ref>.

[step=2cm] (4,4) grid +(4,4)
[step=2cm] (0,4) grid +(4,4)
[step=1cm] (2,4) grid +(2,2)
[step=1cm] (4,4) grid +(2,2)
[step=2cm] (4,0) grid + (4,4)
[step=0.5cm] (4, 2) grid +(1,1)
[step=0.5cm] (4, 3) grid +(1,1)
[step=1cm] (5, 2) grid +(1,1)
[step=0.5cm] (2, 4) grid +(1,1);
[step=4cm, semithick, line cap=round] (0,0) grid (8,8);

[fill] (0.88, 1.77) circle (0.06cm);
[fill] (0.43, 4.65) circle (0.06cm);
[fill] (0.77, 7.23) circle (0.06cm);
[fill] (2.11, 4.13) circle (0.06cm);
[fill] (2.41, 4.81) circle (0.06cm);
[fill] (2.75, 4.25) circle (0.06cm);
[fill] (2.70, 4.60) circle (0.06cm);

[fill] (3.60, 4.44) circle (0.06cm);
[fill] (3.44, 5.77) circle (0.06cm);

[fill] (4.20, 2.78) circle (0.06cm);
[fill] (4.65, 2.42) circle (0.06cm);
[fill] (4.33, 3.27) circle (0.06cm);
[fill] (4.26, 3.88) circle (0.06cm);
[fill] (5.18, 2.77) circle (0.06cm);
[fill] (5.80, 3.13) circle (0.06cm);

[fill] (7.22, 0.8) circle (0.06cm);
[fill] (6.77, 3.2) circle (0.06cm);

[fill] (4.64, 4.76) circle (0.06cm);
[fill] (5.72, 5.31) circle (0.06cm);
[fill] (6.22, 4.81) circle (0.06cm);
[fill] (6.4, 6.81) circle (0.06cm);

The particle distribution of 21 points in two dimensions contained in a quadtree.

The particle distribtuion of 20 points in three dimensions contained in an octree.

With each particle assigned to a leaf node, each method now beings parsing the tree in a bottom-up fashion.
Starting from the leaves, for each level in the tree, a single representative particle
is computed for the ensemble of particles in the child (and grandchild, etc.) nodes below it.
To compute this representative particle,
one can simply place it at the the centre of mass of the child particles and
set its mass equal to the sum of masses of its children.
This is the original formulation of the Barnes-Hut method <cit.>.
However, in the FMM method and later variations of the Barnes-Hut method <cit.>,
the representative particle is instead described by a multipole expansion
of the gravitational potential of the ensemble of points. Recall
that the centre of mass formulation is equivalent to only the monopole term
in the multipole expansion, and including more terms is only required for increased
accuracy in approximation; see Section <ref>.

When including terms in the multipole expansion beyond the centre of mass, this representative particle,
or pseudo-body, at each level is no longer really a single particle
but rather a formula expressing the mass distribution of its children.
For a better description we then refer to such a pseudo-body as simply a cell, or node,
since there is a one-to-one correspondence between these
and the nodes or cells of a tree structure.[By both
node and cell we refer to the discrete data units of a tree structure.
Cell is more common in the context of treecodes, while node is more common for generic data structures.]
In treecodes, the tree structure thus fulfills two aspects:

* a spatial decomposition of the bodies being simulated; and

* succinctly encoding hierarchical groupings of bodies and the approximations of those mass distributions.

With the tree structure fully constructed,

it is now used for force calculations for each body.
Recall that the key behaviour of treecodes is to approximate the interaction between
a body and an ensemble of bodies.
Hierarchical ensembles are precisely what are encoded as the internal cells of the tree.

To accurately approximate the interaction between a body and an internal cell,
they must be “far enough away” or well-separated.
This notion depends on the particular method.

* In the FMM method, a cell $d$ is statically determined to be well-separated from another cell $c$ if the separation is greater
than the side length of $c$. Note that this definition is not symmetric.

* In the Barnes-Hut method, a parameter $\theta$, called the opening angle or multipole acceptance criterion (MAC), controls whether a cell should be considered well-separated. If the ratio between a cell's side length $\ell$ and the separation $d$ between a point and the cell's centre of mass is less than $\theta$ (i.e. $\nicefrac{\ell}{d} < \theta$), then they are well-separated. This is depicted in Fig. <ref>.

angles, quotes, decorations.pathreplacing

[step=2cm] (4,4) grid +(4,4)
[step=2cm] (0,4) grid +(4,4)
[step=1cm] (2,4) grid +(2,2)
[step=1cm] (4,4) grid +(2,2)
[step=2cm] (4,0) grid + (4,4)
[step=0.5cm] (4, 2) grid +(1,1)
[step=0.5cm] (4, 3) grid +(1,1)
[step=1cm] (5, 2) grid +(1,1)
[step=0.5cm] (2, 4) grid +(1,1);
[step=4cm, semithick, line cap=round] (0,0) grid (8,8);

[fill] (0.88, 1.77) circle (0.06cm);
[fill] (0.43, 4.65) circle (0.06cm);
[fill] (0.77, 7.23) circle (0.06cm);
[fill] (2.11, 4.13) circle (0.06cm);
[fill] (2.41, 4.81) circle (0.06cm);
[fill] (2.75, 4.25) circle (0.06cm);
[fill] (2.70, 4.60) circle (0.06cm);

[fill] (3.60, 4.44) circle (0.06cm);
[fill] (3.44, 5.77) circle (0.06cm);

[fill] (4.20, 2.78) circle (0.06cm);
[fill] (4.65, 2.42) circle (0.06cm);
[fill] (4.33, 3.27) circle (0.06cm);
[fill] (4.26, 3.88) circle (0.06cm);
[fill] (5.18, 2.77) circle (0.06cm);
[fill] (5.80, 3.13) circle (0.06cm);

[fill] (7.22, 0.8) circle (0.06cm);
[fill] (6.77, 3.2) circle (0.06cm);

[fill] (4.64, 4.76) circle (0.06cm);
[fill] (5.72, 5.31) circle (0.06cm);
[fill] (6.22, 4.81) circle (0.06cm);
[fill] (6.4, 6.81) circle (0.06cm);

(particle) at (0.77, 7.23);
(leftcorner) at (3,4);
(rightcorner) at (2,5);
[draw=red, thick] (leftcorner) rectangle (rightcorner);

[blue, |-|, thick] (2.48, 4.4425) – node[above=0.09, xshift=0.9, black] $\,d_1$ (particle);
[decorate,decoration=brace, amplitude=3pt, raise=2pt] (2,4) to node[black,midway,left= 0.12cm] $\ell_1$ ++(0, 1);
[fill=yellow] (0.77, 7.23) circle (0.06cm);

[red, thick] (4,0) rectangle (8,4);
[decorate, decoration=brace, amplitude=5pt, raise=2pt, aspect=0.75] (4,0) to node[black, pos=0.75, left=0.15cm] $\ell_2$ (4,4);
[|-|, thick, blue] (0.88, 1.77) – node[above, black, pos=0.4] $d_2$ (5.322, 2.661);
[fill=purple] (0.88, 1.77) circle (0.06cm);

[blue, thick, |-|] (5.72, 5.31) – (6.22, 4.81);
[fill=green] (5.72, 5.31) circle (0.06cm);
[fill=green] (6.22, 4.81) circle (0.06cm);

Determining particle-particle and particle-cell interactions in the Barnes-Hut method. The yellow point is well-separated from the cell with side length $\ell_1$ if $\nicefrac{\ell_1}{d_1} < \theta$. The purple point is well-separated from the large cell with side length $\ell_2$ if $\nicefrac{\ell_2}{d_2} < \theta$. The non-particle endpoints of the two line segments are at the centre of mass of the respective ensemble of points. The two green points are not well-separated and their close-range interaction must be computed directly.

These two definitions also imply the second major difference between the two methods.
While both methods directly compute the so-called particle-particle interactions when two particles are close together,
the approximation of long-range interactions differ. In the Barnes-Hut method, all long-range interactions
are computed as particle-cell interactions. That is, a cell contributes to the net force acting on a particle.
In contrast, in FMM, cell-cell interactions are possible, where all bodies in an ensemble experience the same exerted force by
another ensemble. This additional type of interaction allows the FMM method to reach an asymptotic complexity of
$\mathcal{O}(N)$. Despite this attractive computational complexity, this method is more mathematically difficult
and more challenging to implement <cit.>.
Further, particularly in three-dimensions, the constants discarded by the big-$\mathcal{O}$ notation in fact
dominate running time for practical values of $N$ <cit.>.
For these reasons, we focus on the Barnes-Hut algorithm.
We refer the reader to <cit.>
for further information on the FMM method.

In the Barnes-Hut method, the forces acting on each body are
computed via a top-down tree traversal which can be implemented recursively, starting at the root node.
Let $\ell$ and $d$ be defined as before for the current node.
If the multipole approximation is acceptable, $\nicefrac{\ell}{d} < \theta$
(i.e. the particle is far enough away from the cell),
then a particle-cell interaction is computed and the traversal along
this branch of the tree is complete.
Otherwise, the cell is “opened” and the traversal continues down
each of the current node's children.
If a leaf node is encountered, a typical particle-particle interaction is computed.

The entire Barnes-Hut method is detailed in Algorithm <ref>.
It beings by building the tree which describes the current positions of bodies (Lines 2–7).
Notice that the tree construction here proceeds differently as previously mentioned.
Rather than loading all particles into the tree and then recursively dividing cells until the leaves
contain exactly one particle, the tree is constructed one particle at a time,
splitting cells and pushing particles down the tree until each particle is alone in its leaf node.
Once the tree is constructed, it is parsed bottom-up to compute the internal cell multipole moments (Lines 8–9).
Next, forces acting on each particle are computed via a tree traversal (Lines 13–22), where
particle-particle interactions or particle-cell interactions are computed as appropriate. This is repeated
for every body in the simulation.
Finally, just as in the direct method, velocities and positions are updated (Lines 23–24).

We conclude the discussion of the classic Barnes-Hut method with a final remark.
The essential benefit of the Barnes-Hut method is
that to traverse the tree and compute the forces on a single particle in this hierarchical manner
takes time $\mathcal{O}(\log N)$. Thus, to compute forces acting on all $N$ particles is $\mathcal{O}(N\log N)$.
First, notice that constructing the tree also takes time $\mathcal{O}(N\log N)$ since $N$ bodies
are inserted, on average, to a depth of $\log N$ in the tree.
This depth arises from the fact that the average size of a leaf node is on the order of
the average inter-particle spacing. The best case is a uniform distribution of particles
which creates a balanced tree with $\log_8 N$ levels.
Then, traversing the tree to compute the forces acting on a particle is also
on the order of $\mathcal{O}(\log N)$. In the worse case, the target particle is “close”
to every other particle and no approximations can be made, thus
computing $N$ interactions for that target particle. However, it is not possible
for every particle to be close to every particle. Indeed,
for branches of the tree not including the target particle,
it is much more common for those branches to be approximated with a single
particle-cell interaction. Since well-separateness is controlled
by the parameter $\theta$, the amount of work spent in
the force calculation phase is dependent on $\theta$.
It is found that the force calculation for a single particle
typically scales as $\nicefrac{1}{\theta^2}\log N$ <cit.>;
a smaller value of $\theta$ requires more cells to be opened and thus more interactions to be computed.
Typical values of $\theta$ are on the order of $\nicefrac{1}{2}$ <cit.>.

Barnes-Hut Simultation Loop Body($P$, $\theta$)Given a list of particles $P$ and a MAC criterion $\theta$, compute forces and update positions for the current time-step. This is repeated for each time-step in the simulation.
[1]
Build Tree:
Create a tree $T$ with a single root node that is the bounding box of all particles.
each particle $p$ in $P$ Spatial decomposition of bodies
Insert $p$ into the leaf node $L$ of $T$ which encloses the point's coordinates.
$L$ now holds two particles
Split $L$ in half in each of its dimensions to create children, $L$ is now an internal node.
Insert each particle of $L$ into the appropriate child.
Repeat until each particle is in a separate leaf node.

each node in $T$ from leaves to root Compute mass distribution approximations
Compute sum of masses of contained points, centre of mass, and multipole values.

Compute Forces:
each particle $p$ in $P$
Set the acceleration of p to $\vec{0}$.
$Node$ $\gets$ the root node of $T$

Compute Cell Force:
each child cell $C$ of $Node$
$C$ contains 1 particle $C$ is a leaf
Add to $p$'s acceleration the contribution from $C$'s particle. particle-particle
$\ell$ $\gets$ $C$'s side length
$d$ $\gets$ distance between $C$'s centre of mass (c.o.m.) and $p$'s position.
${\ell} / {d} < \theta$
Add to $p$'s acceleration the contribution from $C$. particle-cell
Compute Cell Force with $Node \gets C$ recursive call

Integration step (update velocities and positions):
each particle $p$ in $P$
Update $p$'s position and velocity using its acceleration.

§ HIGH-PERFORMANCE AND PARALLEL COMPUTING

In order to better understand the design and implementation of a parallel
Barnes-Hut method, we must discuss aspects of high-performance computing
and parallel computing. In this section we give a brief overview of some core concepts
including multithreading, data locality, parallel task granularity and load balancing,
and irregular parallelism. For further details
readers should see <cit.> and <cit.>.

§.§.§ Multithreading

The multithreading programming and execution model is a way for programs to achieve thread-level parallelism (TLP).
TLP is achieved by creating multiple independent and concurrent threads of execution,
each usually executing an independent task within the program.
True concurrency, although, is only achieved where the computer hardware executing a program supports
multithreading through multiple physical threads. Typically, one physical thread is
one core (i.e. processor) within a multi-core processor or a multiprocessor.
With simultaneous multithreading (SMT), each physical core can support several physical threads (often only 2).
However, simultaneous multithreading does not provide the same level of parallelism as distinct cores,
and brings additional overheads. SMT with two physical threads on a single core can
give a parallel speedup between 0.9 (a slowdown) and 1.6 for various applications <cit.>.
Overall, SMT can give a slight boost to parallel speedup over the same number of cores without SMT.
This parallel speedup is defined as the ratio between a completely serial execution and the parallel execution.
Ideal speedup, or linear speedup, is obtained when the parallel speedup equals the number of physical threads.

§.§.§ Data Locality

On modern computer architectures with a cache memory hierarchy,
a key consideration for program performance is data locality.
This refers to the way a program accesses and traverses data
during the program's execution. Two types of locality should be considered:

* temporal locality refers to when data which was recently accessed is accessed again, and

* spatial locality refers to when data is accessed that is close to recently accessed data.

Here, data being close refers to their memory positions, e.g. adjacent elements in an array.
Good data locality avoids expensive cache misses where data must be retrieved
from main memory as opposed to the processor's local cache of data.
The processor-memory gap <cit.>,
tells us that modern processors are exponentially faster than modern random-access
memory, and thus data locality is important for performance.
Key considerations to obtain data locality include accessing an array in a linear order (spatial locality),
and accessing only a small amount of the overall data at one time (temporal locality).
Data locality is made further important in the case of multithreading. Independent threads
should access independent sections of memory.
Otherwise overhead is introduced where coherency of data
must be maintained by constantly sharing data between each core's local cache.

§.§.§ Task Granularity, Load Balancing

An important consideration when implementing thread-level parallelism is task granularity and load balancing.
Task granularity refers to the amount of work to be executed by a thread. Since creating or spawning
a thread takes time, as does scheduling threads for execution and copying data into each core's local cache,
the amount of work to perform must be substantial
to warrant those parallelism overheads.

To achieve optimal parallel speedup, the amount of work performed by each thread should also be equal.
Otherwise, the overall program performance is limited by the thread performing the most amount of work.
The act of spreading work evenly between available threads is known as load balancing.

§.§.§ Irregular Parallelism

Lastly, the parallel execution of a program can be categorized into one of regular parallelism or irregular parallelism <cit.>.
Regular parallelism describes situations where tasks to be executed in parallel are similar and predictable.
Each task performs a known amount of work and load balancing is easy. In contrast, dissimilar tasks
with unpredictable amounts of work lead to irregular parallelism.
Problems with irregular parallelism are difficult to load balance.
Moreover, the ability of a program to scale—to increase parallel speedup given more processors—
is further challenging.

In treecodes, the number of interactions to be computed for each body changes at
every time-step; exhibiting irregular parallelism.
In contrast, updating positions and velocities requires
a consistent and known amount of work per body, yielding regular parallelism.

§ PARALLEL IMPLEMENTATION OF $N$-BODY SIMULATIONS AND TREECODES

To gain further performance in $N$-body simulations,
and hierarchical methods for such simulations,
parallel execution is a natural choice.
Toward that goal we will examine two pioneering works:
the costzones approach of Singh et al. <cit.>, and the hashed octree method
of Warren and Salmon <cit.>.
The costzones approach was first developed for a shared-memory architecture, meanwhile
the hashed octree method was explicitly developed for a distributed system.
Nonetheless, all the techniques we will discuss can be applied to either a shared-memory system
or a distributed system thanks to the message passing model being agnostic to the underlying
communication protocol.

We begin this section with a high-level overview of the
designs and parallelization techniques of the costzones and hashed octree methods.
The details follow in the below subsections.

Just as with any parallel algorithm, both methods are concerned with:

* minimizing synchronization—where processes must wait for all others before proceeding;

* load-balancing—ensuring each process has an equal amount of work, thus reducing wait times at synchronization points;

* reducing communication overheads—the amount data and number of messages to be shared between processes; and

* reducing span—the maximum amount of work done by any one process.

Reducing span means parallelizing as many parts of the algorithm as possible.
Both the costzones and hashed octree algorithms attempt to parallelize the Barnes-Hut algorithm by:

* parallelizing the building of the octree;

* parallelizing the force calculation on each body, i.e. parallelizing the tree traversal; and

* parallelizing the integration step.

An obvious and simple approach to parallelizing an $N$-body simulation is to simply
assign $\nicefrac{N}{s}$ bodies to each of the $s$ parallel processes.
This would be sufficient for the integration step, where each body's
position and velocity is easily updated from its own acceleration.
However, due to the mutual interaction between all bodies,
the computation of that acceleration requires a global view of the entire simulation
domain, and thus a global view of the octree.
Both methods have a solution to this global view requirement.

In the costzones method, each process, in parallel, builds a so-called local octree
for its assigned bodies. Then, these local octrees are merged into a global octree which is then
shared with all processes. Finally, with each process having a global view of the simulation domain,
each process can proceed in parallel to compute the forces on each of its assigned bodies.

In the hashed octree method, again, each process builds its local octree in parallel
for its assigned bodies. This method, however, explicitly avoids creating a single global tree.
Instead, during the tree traversal and force calculations, when a local particle
needs information about the position or multipole values of a non-local particle or cell,
that information is requested dynamically from the owning process. An optional improvement
is to use asynchronous communication and a latency-hiding traversal to avoid
stalling computations while waiting for communications (see Section <ref>).
To facilitate consistent addressing of octree cells stored on different processes,
the octree is implemented as a hashed octree or linear octree <cit.>;
hence the naming of this method.

Both methods are in fact very similar in design. They differ only in how a process
receives octree information for non-local bodies.
The costzones approach uses upfront communication to construct and share the global octree,
and then proceeds with tree traversals in parallel, without any communication.
On the other hand, the global tree is not strictly needed to calculate
the force on any one particle due the particle-cell interactions of the Barnes-Hut method.
The hashed octree method thus avoids computing the entire global tree at the beginning
and instead uses extra communication during the tree traversal to obtain only
the essential data for its assigned bodies.
In principle, this latter approach is the most scalable, where simulation sizes may
grow so large that the entire domain cannot be stored in the memory of a single process.

Nonetheless, both techniques share many similarities to reduce span and improve load-balancing.
Firstly, a spatial decomposition of the
simulation domain is performed. Each process is then assigned a partition of the spatial domain,
and the bodies it contains, rather than statically assigning $\nicefrac{N}{s}$ bodies to each process.
Having the bodies assigned to one process be spatially local greatly improves the
tree merging of the costzones method
and reduces the number of communications required during the 	tree traversal
step of the hashed octree method. This latter fact is obvious
considering when particle-particle interactions are computed in a Barnes-Hut algorithm.
Secondly, dynamic load-balancing is achieved by re-partitioning the spatial domain at every time step.
Particles in regions dense with others require more work during force computation
due to the higher number of particle-particle interactions. Those spatial partitions
should be smaller and contain fewer particles so that each process
computes roughly the same number of interactions.
Thirdly, the spatial partitioning of bodies is implemented by linearizing three-dimensional
space via a space-filling curve, sorting the list of bodies based on this linearization,
and then partitioning the sorted list.
This sort can be parallelized to further reduce span.

Clearly, there are many different opportunities for parallelization in
the Barnes-Hut algorithm and many different algorithmic techniques which can be implemented
to improve this parallelization through improved load-balancing, etc.
We now detail each parallelization and technique, one at a time, in the following subsections.
These techniques are added incrementally to our implementation
to step-by-step improve the parallelization and scalability of the algorithm;
see Section <ref> for that experimentation.

We begin in Section <ref> by describing a simple parallelization
of the direct method and the parallelization of the integration step.
This serves as a point of comparison for the remaining methods to show
that the additional work of creating and using the octree is worthwhile.
Section <ref> presents a first parallelization
of the Barnes-Hut method where parallel force calculation via
tree traversal is added to the previously described parallel integration step.
Next, building the octree in parallel and the octree merging process is detailed in Section <ref>.
The method of spatial decomposition is then detailed in Section <ref>.
Spatial decomposition is further enhanced by dynamic load-balancing in Section <ref>
and by distributed sorting in Section <ref>.
Section <ref> introduces the hashed octree and the techniques required to
build local octrees which are globally consistent.
Finally, we add latency-hiding asynchronous communication to the hashed octree tree traversal
in Section <ref>.
Therefore, the algorithm described in Section <ref> describes the “final”
version of our costzones-based method, meanwhile the algorithm described in Section <ref>
describes the “final” version of our hashed octree method.

In the follow sections we use common terminology from parallel processing
and distributed computing, e.g., gather, scatter, and broadcast. For definitions
and further details see <cit.>.

§.§ Naive Parallelization of Integration and Particle-Particle Interactions

We have already seen the serial direct method as
Algorithm <ref> in Section <ref>.
We now formulate this algorithm as a distributed algorithm,
as shown in Algorithm <ref>.
We begin by partitioning the bodies evenly across the $s$ processes.
No care is taken for spatial decomposition because all $\mathcal{O}(N^2)$ interactions
will be computed.
Next, we compute the mutual forces on the local bodies by the other local bodies (Lines 6–8).
Each process sets a temporary list of particles to its own list of local particles.
Then, $s - 1$ rounds of message passing occur (Lines 9–14). Assuming a ring topology,
each process sends its current temporary list of particles to the process on the right
and receives a new list from the process on the left. The force
from each particle in this temporary list is then added to each particle
in the local list. Then, the next round begins. After $s-1$ rounds process $id-1$ has
received process $id$'s particles and all other particles. This concludes
the force calculation in this time-step.
Each process then updates the positions and velocities of each local particle (Lines 15–16).
In this algorithm and all algorithms which follow, we assume the processes are ranked $0,\ldots,s-1$.

Parallel Direct $N$-Body Simulation($id$, $s$, $\Delta t$, $t_{end}$, $N$, $P_{all}$)$P_{all}$ is the list of size $N$ of bodies to simulate from 0 to $t_{end}$ going by $\Delta t$ time-steps. $id$ is the executing process's id, $s$ is the number of processes. Assume $s \mid N$.
[1]
$N' \gets N / s$
$P \gets$ Scatter $P_{all}$, sending to process $id$ $P_{all}[id\cdot N', \ldots, id\cdot N' + N' - 1]$ Get Local Bodies
$t$ from 0 to $t_{end}$ by $\Delta t$
$P' \gets P$
for each particle $p$ in $P$ set acceleration to $\vec{0}$.
each particle $p$ in $P$ Force Calculation with Locals
each particle $q$ in $P \setminus \{p\}$

Add the force contribution of $q$ acting on $p$ to the acceleration of $p$.
$i$ from 0 to $s-1$
Send $P'$ to process $id + 1 \mod s$
$P' \gets $ particles received from process $id - 1 \mod s$
each particle $p$ in $P$ Force Calculation with Remotes
each particle $q$ in $P'$

Add the force contribution of $q$ acting on $p$ to the acceleration of $p$.
each particle $p$ in $P$ Position Update
Update position and velocity of $p$ using its current acceleration.

In this algorithm, each process computes only $N'(N'-1) + N'(N')(s-1) = \nicefrac{1}{s}N(N-1)$ interactions.
Further, each process integrates for only $\nicefrac{N}{s}$ particles.
The span and load-balance of this algorithm are thus optimal, since the total work is $N(N-1)$ interactions
plus $N$ integrations. Nonetheless, the number of interactions still scales as $\mathcal{O}(N^2)$.
We thus look to parallelize the Barnes-Hut method with its $\mathcal{O}(N\log N)$ number of interactions.

§.§ Parallel Tree Traversal

Having seen our first distributed algorithm, we are now ready to modify
the Barnes-Hut algorithm to be executed in parallel.
As seen at the beginning of Section <ref>,
there are three aspects to parallelize: $(a)$ the building of the octree,
$(b)$ the force calculation, and $(c)$ the integration step.
Parallelizing the integration is embarrassingly parallel and
is as simple as assigning a subset of the bodies to each process as
was just seen in Algorithm <ref>.
The next simplest procedure is to parallelize the force calculation.
In the Barnes-Hut method, if the octree encompassing all bodies in the simulation has already been constructed,
then parallel force calculation is also embarrassingly parallel as
each process simply traverses the octree for each of its local bodies.

Let us begin with a short review of the implementation of an octree.
A simple recursive implementation based on pointers is a direct generalization
of a rooted binary tree based on linked lists;
see, e.g., <cit.>.
Each cell in an octree is represented by its data or payload and 8 pointers
to its children. Since the octree is deeply connected with the spatial domain,
each node is also affixed with its spatial centre and its
size—the side length of the cell.
The payload of each node, for the purpose of $N$-body simulation,
also includes the position and mass of the body
it contains, if a leaf node, or the center of mass and multipole coefficients,
if an internal node; see Section <ref>.
Fig. <ref> shows the encoding
of an octree with three spatial divisions as a rooted tree of

[circle, semithick, draw, scale=1.2] (root) at (9,4) ;
in 5,...,12
[circle, draw, scale=1.2, blue] (-1) at (,2) ;
in 6,...,13
[circle, draw,scale=1.2, green!80!black] (-2) at (.5,0) ;
in 5,...,12
[->, >=Stealth] (root.south) – (-1);
in 6,...,13
[->, >=Stealth] (10-1.south) – (-2);

(O) at (0,0,0);
(A) at (0,,0);
(B) at (0,,);
(C) at (0,0,);
(D) at (,0,0);
(E) at (,,0);
(F) at (,,);
(G) at (,0,);

[black, semithick, line cap=round] (O) – (C) – (G) – (D) – cycle;
[black, semithick, line cap=round] (O) – (A) – (E) – (D) – cycle;
[black, semithick, line cap=round] (O) – (A) – (B) – (C) – cycle;
[black, semithick, line cap=round] (D) – (E) – (F) – (G) – cycle;
[black, semithick, line cap=round] (C) – (B) – (F) – (G) – cycle;
[black, semithick, line cap=round] (A) – (B) – (F) – (E) – cycle;

[blue] (2, 0, 0) – (2, 4, 0);
[blue] (0, 2, 0) – (4, 2, 0);

[blue] (4, 2, 0) – (4, 2, 4);
[blue] (4, 0, 2) – (4, 4, 2);

[blue] (2, 0, 0) – (2, 0, 4);
[blue] (0, 0, 2) – (4, 0, 2);

[blue] (0, 2, 0) – (0, 2, 4);
[blue] (0, 0, 2) – (0, 4, 2);

[blue] (2, 4, 0) – (2, 4, 4);
[blue] (0, 4, 2) – (4, 4, 2);

[blue] (2, 0, 4) – (2, 4, 4);
[blue] (0, 2, 4) – (4, 2, 4);

[blue] (0, 2, 2) – (4, 2, 2);
[blue] (2, 2, 0) – (2, 2, 4);
[blue] (2, 0, 2) – (2, 4, 2);

[fill=blue!50,fill opacity=0.08, draw opacity=0] (0,4,0) – (4,4,0) – (4,0,0) – (4,0,4) – (0,0,4) – (0,4,4) – (0,4,0);
[fill=blue!50,fill opacity=0.08,draw opacity=0] (0,2,0) – (4,2,0) – (4,2,4) – (0,2,4) – cycle;
[fill=blue!50,fill opacity=0.08,draw opacity=0] (0,0,2) – (4,0,2) – (4,4,2) – (0,4,2) – cycle;
[fill=blue!50,fill opacity=0.08,draw opacity=0] (2,0,0) – (2,4,0) – (2,4,4) – (2,0,4) – cycle;

[blue,fill=green!50, opacity=0.1] (2,2,2) – (2,4,2) – (2,4,4) – (2,2,4) – cycle;
[blue,fill=green!50,opacity=0.1] (2,2,2) – (4,2,2) – (4,4,2) – (2,4,2) – cycle;
[blue,fill=green!50,opacity=0.1] (2,2,2) – (4,2,2) – (4,2,4) – (2,2,4) – cycle;

[green!80!black] (3, 2, 2) – (3, 4, 2);
[green!80!black] (2, 3, 2) – (4, 3, 2);

[green!80!black] (4, 3, 2) – (4, 3, 4);
[green!80!black] (4, 2, 3) – (4, 4, 3);

[green!80!black] (3, 2, 2) – (3, 2, 4);
[green!80!black] (2, 2, 3) – (4, 2, 3);

[green!80!black] (2, 3, 2) – (2, 3, 4);
[green!80!black] (2, 2, 3) – (2, 4, 3);

[green!80!black] (3, 4, 2) – (3, 4, 4);
[green!80!black] (2, 4, 3) – (4, 4, 3);

[green!80!black] (3, 2, 4) – (3, 4, 4);
[green!80!black] (2, 3, 4) – (4, 3, 4);

[green!80!black] (2, 3, 3) – (4, 3, 3);
[green!80!black] (3, 3, 2) – (3, 3, 4);
[green!80!black] (3, 2, 3) – (3, 4, 3);

[fill=green!50,fill opacity=0.2,draw opacity=0] (2,3,2) – (4,3,2) – (4,3,4) – (2,3,4) – cycle;
[fill=green!50,fill opacity=0.2,draw opacity=0] (2,2,3) – (4,2,3) – (4,4,3) – (2,4,3) – cycle;
[fill=green!50,fill opacity=0.2,draw opacity=0] (3,2,2) – (3,4,2) – (3,4,4) – (3,2,4) – cycle;

[decorate,decoration=brace, amplitude=8pt, raise=20pt] (0,0,4) to node[black,midway,left=0.83] $\ell_1$ (0, 4, 4);
[decorate,decoration=brace, amplitude=8pt, raise=2pt, blue] (0,2,4) to node[black,midway,left=0.22,blue] $\ell_2$ (0, 4, 4);
[decorate,decoration=brace, amplitude=4pt, raise=2pt, mirror,green!80!black] (4,3,2) to node[black,midway,right=0.15,green!80!black] $\ell_3$ (4, 4, 2);

An octree with three levels of spatial division, and its corresponding encoding as
a rooted tree. The root is black, with side length $\ell_1$; this side length
is the extents of the entire simulation domain. The next division is blue, with side length $\ell_2$,
and the next green, with side length $\ell_3$.

Algorithm <ref> shows the Barnes-Hut method with simple parallel tree traversal.
It begins by assigning $\nicefrac{N}{s}$ bodies to each process. Then, the root process (id 0),
gathers all bodies (since their positions were updated in the previous simulation step).
The root process simply builds the global tree from all bodies at once. This global tree
is then broadcast to all other processes (Lines 5–7).
With each process having the global octree, they can proceed in parallel to
compute the forces on their local bodies with respect to the entire octree (Lines 8-20).
Integration steps then proceed in parallel as before (Lines 21-22).

Parallel Barnes-Hut V1($id$, $s$, $\Delta t$, $t_{end}$, $N$, $P_{all}$)$P_{all}$ is the list of size $N$ of bodies to simulate from 0 to $t_{end}$ going by $\Delta t$ time-steps. $id$ is the executing process's id, $s$ is the number of processes. Assume $s \mid N$.
[1]
$N' \gets N / s$
$P \gets$ Scatter $P_{all}$, sending to process $id$ $P_{all}[id\cdot N', \ldots, id\cdot N' + N' - 1]$ Get Local Bodies
$t$ from 0 to $t_{end}$ by $\Delta t$
$P_{all} \gets$ Gather $P$ to process 0.
$id = 0$
$T \gets$ build octree from $P_{all}$
Broadcast $T$ from process 0 to all others.

each particle $p$ in $P$ Compute Forces
Set the acceleration of $p$ to $\vec{0}$.
$Node$ $\gets$ the root node of $T$

$\quad\ $Compute Cell Force:
each child cell $C$ of $Node$
$C$ contains 1 particle $C$ is a leaf
Add to $p$'s acceleration the contribution from $C$'s particle. particle-particle
$\ell$ $\gets$ $C$'s side length
$d$ $\gets$ distance between $C$'s c.o.m. and $p$'s position.
${\ell} / {d} < \theta$
Add to $p$'s acceleration the contribution from $C$. particle-cell
Compute Cell Force with $Node \gets C$ recursive call
each particle $p$ in $P$ Position Update
Update position and velocity of $p$ using its current acceleration.

Notice that for this algorithm to fit within the message passing model for
distributed algorithms, we must have a way to serialize
the octree into a message. Since its encoding is based on pointers and linked lists,
this is not immediately obvious.
One solution, which we implement,
is to add to each node an 8-bit mask whose bits encode
whether any of the 8 children actually exist.
Then, serializing the octree is as simple as serializing each node
and concatenating them in the order of pre-order tree traversal.
Deseralizing is also made simple by following a pre-order traversal
in combination with the deserialized bitmask.

In this simple parallel Barnes-Hut algorithm, the span
is not very good since the root process must construct the entire octree
in time $\mathcal{O}(N\log N)$. The force calculation is still done
in parallel, but now taking time $\mathcal{O}(\nicefrac{N}{s}\log N)$ thanks to the
octree and particle-cell interactions. Updating the particle positions
and velocities remains at $\nicefrac{N}{s}$ integrations.

§.§ Parallel Octree Construction and Merging

Parallelizing the octree construction is a non-trivial task.
Naive attempts, such as all processes simultaneously inserting their
assigned bodies into a single global octree, requires an excessive amount of
synchronization. This technique was explored, and abandoned, in <cit.>.
Their improved variation involves each process creating
its own local octree, from its local bodies, and then merging
trees together.
One fundamental aspect of this merge operation is that each local tree
must cover the same spatial domain. That is, their root nodes
must have the same size. Then, all cells in all of the local trees
have the same sizes. This allows entire subtrees to be inserted into
a merged tree where there is an otherwise empty child node.
Algorithm <ref> presents this merge algorithm;
Lines 5–6 and Lines 11-12 show this subtree replacement.

Octree Merge($T_1$, $T_2$)Merges the two octrees rooted at nodes $T_1$ and $T_2$ by moving data from $T_2$ into $T_1$. This assumes both $T_1$ and $T_2$ have at least one level of children.
[1]
$i$ from 0 to 7
$T_1$ and $T_2$ both have an $i$th child
$(C_1,C_2) \gets$  (the $i$th child of $T_1$, the $i$th child of $T_2$)
$C_1$ is a leaf node
Insert the particle of $C_1$ into $C_2$ see Build Tree of Algorithm <ref>
Replace $C_1$ in $T_1$ with $C_2$
$C_2$ is a leaf node
Insert the particle of $C_2$ into $C_1$
Octree Merge($C_1$, $C_2$) Recursive Call
$T_2$ has an $i$th child $T_1$ does not have an $i$th child
Set the $i$th child of $T_1$ to $C_2$

In the original octree merge of <cit.>,
each process's local tree was merged one at a time into the global tree.
This was feasible in a shared-memory model.
In a message passing model, or in any more generic parallel scheme,
this would require excessive communication and synchronization.
Instead, we perform merges of local trees pairwise, in a reduction (see <cit.>),
until a single tree contains all $N$ bodies. This single tree is then broadcast to all other processes.
Note that while this is phrased in terms of a Map-Reduce pattern, the reduction in fact
does not “reduce” the size of data, but moreso performs a gather of the individual local trees.
Algorithm <ref> shows this scheme, where the parallel tree building occurs in Lines 4–6,
and the remaining lines are the same as Algorithm <ref>.
Notice that the reduce step hides many serializations, inter-process messages, and deserializations.

Parallel Barnes-Hut V2($id$, $s$, $\Delta t$, $t_{end}$, $N$, $P_{all}$)$P_{all}$ is the list of size $N$ of bodies to simulate from 0 to $t_{end}$ going by $\Delta t$ time-steps. $id$ is the executing process's id, $s$ is the number of processes. Assume $s \mid N$.
[1]
$N' \gets N / s$
$P \gets$ Scatter $P_{all}$, sending to process $id$ $P_{all}[id\cdot N', \ldots, id\cdot N' + N' - 1]$ Get Local Bodies		$t$ from 0 to $t_{end}$ by $\Delta t$
$T' \gets$ build octree from $P$ Parallel Tree Build
$T \gets $ Reduce, via Octree Merge, $T'$ to process 0.
Broadcast $T$ from process 0 to all others.

each particle $p$ in $P$ Compute Forces; unchanged from Algorithm <ref>
each particle $p$ in $P$ Position Update
Update position and velocity of $p$ using its current acceleration.

While the above algorithm parallelizes the tree building step,
and building each local octree now only incurs time $\mathcal{O}(\nicefrac{N}{s}\log(\nicefrac{N}{s}))$,
the span of the reduce still requires $\log(s)$ merges to occur. In the worst case,
where each process's local bodies are uniformly distributed across
the entire spatial domain, merging octrees is a laborious process
where both trees must be traversed all the way down to their leaves
(see the recursive call in Algorithm <ref>),
and then leaf nodes inserted into the other subtree
(see Build Tree step of Algorithm <ref>).
A single merge of two trees with $n$ bodies each thus requires time
$\mathcal{O}(n\log n)$, resulting in a
span of the parallel octree build of $\mathcal{O}(\log s (N\log N))$.
This does not exhibit strong scaling—where the span
would reduce with an increased number of processes—since the
dominant term is still $N\log N$.
A better parallel solution would result in a span
with a dominant term of $\mathcal{O}(\nicefrac{N}{s}\log(\nicefrac{N}{s}))$

§.§ Spatial Decomposition via Morton's Ordering

A key consideration in both the
costzones and hashed octree methods is the spatial decomposition of bodies
rather than static assignment of bodies to processes.
Both methods are based on the idea that each process
is responsible for a spatial subdomain of the simulation, and the bodies
contained therein, rather than a set of bodies
spread across the entire domain.
In the costzones method this allows for dynamic load-balancing into
the namesake costzones (see Section <ref>)
and improves the octree merge, as we will soon see.
In the hashed octree method, the spatial decomposition minimizes
inter-process communication (see Section <ref>).

However, one key challenge for $N$-body simulation
is its dynamicity and irregularity. It would be wholly
insufficient to simply divide the spatial domain
into regular pieces, and assign one piece to each process.
As bodies move throughout the spatial domain during the simulation,
some areas would have more bodies, some would have less, and
some may not contain any bodies at all.
This load-balancing problem is further exacerbated by the formation of
highly clustered areas as the gravitational system evolves <cit.>.

A first attempt then at spatial decomposition involves partitioning
the spatial domain so that each process's domain contains
the same number of bodies. This will naturally create different
shapes and sizes for each spatial partition.
First, a linear ordering is given to the leaf nodes of the octree
or, equivalently, the bodies themselves.
This is done through the use of a space-filling curve—Peano-Hilbert
curve in the case of costzones, and Morton ordering in the case of hashed octree.
This linear ordering on the bodies is then partitioned evenly and each partition assigned
to a process.

The more simple Morton ordering, or Z-ordering, involves
transforming a multi-dimensional coordinate tuple into a single integer,
which can then be sorted in a typical integer order.
This single integer spatial key is created by interleaving
the bits of each coordinate's value.
For example, a 2-dimension ($x,y$) coordinate $(3,5) = (0b\textcolor{blue}{0011},0b\textcolor{green!60!black}{0101})$
can be interleaved into the 8-bit spatial key $0b(\textcolor{blue}{0}\textcolor{green!60!black}{0})(\textcolor{blue}{0}\textcolor{green!60!black}{1})(\textcolor{blue}{1}\textcolor{green!60!black}{0})(\textcolor{blue}{1}\textcolor{green!60!black}{1}) = 27$.
Notice this interleave orders the bits such that those from the x-coordinate
are more significant, thus we say this ordering has $x > y$. A different orientation would be
produced with $y > x$, but it would still create a space-filling curve.

With such a bit interleave, in $d$ dimensions, the spatial key
requires $d\cdot b$ bits to encode spatial coordinates
which range over the integers $[0, 2^b)$.
Moreover, spatial coordinates are often encoded using floating point
values, not integers. To remedy this, we perform several mappings
to obtain integer spatial coordinates in the range $[0, 2^{21}-1]$;
$2^{21}$ is chosen so that the three-dimensional spatial key of $3\cdot21=63$ bits can be
stored in a single 64-bit machine word.
First, we must know the extents of the spatial domain. Let us
assume the position of each body falls in the range $(-R,R)$, for some integer $R$,
in each dimension.
We first translate each coordinate of each body to the range $(0,2R)$.
Then, we map the range $(0,2R)$ to $(0,2^{21})$. Notice
this mapping causes a loss of precision, with at most $2^{21}$
discrete values in each dimension. This is sufficient for simulations
with up to several million bodies. In larger simulations,
a 128-bit key is more than sufficient.

[step=1, thin] (0,0) grid (8,8);
[step=2,thick] (0,0) grid (8,8);
[step=4,very thick, line cap=round] (0,0) grid (8,8);

[green!60!black, semithick, rounded corners=5] (0.5,0.5) – (0.5,1.5) – (1.5,0.5) – (1.5,1.5) – (0.5,2.5) – (0.5,3.5) – (1.5,2.5) – (1.5,3.5) – (2.5,0.5) – (2.5,1.5) – (3.5,0.5) – (3.5,1.5) – (2.5,2.5) – (2.5,3.5) – (3.5,2.5) – (3.5,3.5)
– (0.5,4.5) – (0.5,5.5) – (1.5,4.5) – (1.5,5.5) – (0.5,6.5) – (0.5,7.5) – (1.5,6.5) – (1.5,7.5) – (2.5,4.5) – (2.5,5.5) – (3.5,4.5) – (3.5,5.5) – (2.5,6.5) – (2.5,7.5) – (3.5,6.5) – (3.5,7.5) –
(4.5,0.5) – (4.5,1.5) – (5.5,0.5) – (5.5,1.5) – (4.5,2.5) – (4.5,3.5) – (5.5,2.5) – (5.5,3.5) – (6.5,0.5) – (6.5,1.5) – (7.5,0.5) – (7.5,1.5) – (6.5,2.5) – (6.5,3.5) – (7.5,2.5) – (7.5,3.5) –
(4.5,4.5) – (4.5,5.5) – (5.5,4.5) – (5.5,5.5) – (4.5,6.5) – (4.5,7.5) – (5.5,6.5) – (5.5,7.5) – (6.5,4.5) – (6.5,5.5) – (7.5,4.5) – (7.5,5.5) – (6.5,6.5) – (6.5,7.5) – (7.5,6.5) – (7.5,7.5);

[->, >=Stealth] (-0.4, -0.4) – (-0.4, 3);
[font=] at (-0.25, 2.6) $y$;
[->, >=Stealth] (-0.4, -0.4) – (3, -0.4);
[font=] at (2.6, -0.25) $x$;

The Morton ordering of a quadtree with three levels of spatial division. This ordering treats the $x$ dimension as more significant than the $y$ dimension in the bit interleave. $y > x$ would instead create a mirrored “Z” shape.

[fill=blue!80, fill opacity=0.1, draw opacity=0] (0.1,0.1) – (0.1, 7.9) – (1.9, 7.9) – (1.9, 4.9) – (2.9, 4.9) – (2.9,3.9) – (3.9,3.9) – (3.9,0.1) – cycle;
[fill=orange!90!yellow, fill opacity=0.1, draw opacity=0] (3.1,5.1) – (2.1,5.1) – (2.1, 7.9) – (3.1,7.9) – (3.9,7.9) – (3.9,4.1) – (3.1,4.1) – (3.1, 5.1) – cycle;
[fill=orange!90!yellow, fill opacity=0.1, draw opacity=0] (4.1, 3.9) – (4.9, 3.9) – (4.9, 2.9) – (5.9, 2.9) – (5.9, 2.1) – (4.1,2.1) – cycle;

[fill=red, fill opacity=0.1, draw opacity=0] (4.1,1.9) – (4.1, 0.1) – (7.9,0.1) – (7.9, 7.9) – (4.1, 7.9) – (4.1, 4.1) – (5.1, 4.1) – (5.1, 3.1) – (6.1, 3.1) – (6.1, 1.9) – cycle;

[step=2cm, thin] (4,4) grid +(4,4)
[step=2cm, thin] (0,4) grid +(4,4)
[step=1cm, thin] (2,4) grid +(2,2)
[step=1cm, thin] (4,4) grid +(2,2)
[step=2cm, thin] (4,0) grid + (4,4)
[step=0.5cm, thin] (4, 2) grid +(1,1)
[step=0.5cm, thin] (4, 3) grid +(1,1)
[step=1cm, thin] (5, 2) grid +(1,1)
[step=0.5cm, thin] (2, 4) grid +(1,1);
[step=4cm, very thick, line cap=round] (0,0) grid (8,8);

[green!40!black,semithick, rounded corners=2](0.88, 1.77) – (0.43, 4.65) – (0.77, 7.23) – (2.11, 4.1) – (2.41, 4.81)– (2.75, 4.25)– (2.70, 4.60) – (3.60, 4.44) – (3.44, 5.77)

– (4.20, 2.78) – (4.65, 2.42) – (4.33, 3.27) – (4.26, 3.88)– (5.18, 2.77) – (5.80, 3.13) 	– (7.22, 0.8) – (6.77, 3.2) – (4.64, 4.76) – (5.72, 5.31)  – (6.22, 4.81) – (6.4, 6.81);

[fill] (0.88, 1.77) circle (0.06cm);
[fill] (0.43, 4.65) circle (0.06cm);
[fill] (0.77, 7.23) circle (0.06cm);
[fill] (2.11, 4.13) circle (0.06cm);
[fill] (2.41, 4.81) circle (0.06cm);
[fill] (2.75, 4.25) circle (0.06cm);
[fill] (2.70, 4.60) circle (0.06cm);

[fill] (3.60, 4.44) circle (0.06cm);
[fill] (3.44, 5.77) circle (0.06cm);

[fill] (4.20, 2.78) circle (0.06cm);
[fill] (4.65, 2.42) circle (0.06cm);
[fill] (4.33, 3.27) circle (0.06cm);
[fill] (4.26, 3.88) circle (0.06cm);
[fill] (5.18, 2.77) circle (0.06cm);
[fill] (5.80, 3.13) circle (0.06cm);

[fill] (7.22, 0.8) circle (0.06cm);
[fill] (6.77, 3.2) circle (0.06cm);

[fill] (4.64, 4.76) circle (0.06cm);
[fill] (5.72, 5.31) circle (0.06cm);
[fill] (6.22, 4.81) circle (0.06cm);
[fill] (6.4, 6.81) circle (0.06cm);

[->, >=Stealth,opacity=0] (-0.4, -0.4) – (-0.4, 3);
[font=,opacity=0] at (-0.25, 2.6) $x$;
[->, >=Stealth,opacity=0] (-0.4, -0.4) – (3, -0.4);
[font=,opacity=0] at (2.6, -0.25) $y$;

The Morton ordering of 21 points in two dimensions. The spatial domain is partitioned into three pieces according to the cell colors. Notice that the second domain contains a spatial discontinuity.

The Morton ordering on the cells of a quad tree and on a collection of bodies in 2 dimensions.

Fig. <ref> presents the Morton ordering, for $x > y$, for
a quadtree with three levels of spatial division. Fig. <ref>
then shows this Morton ordering applied to a collection of bodies in two-dimensional space;
the line connects the bodies in order. In this same figure, a spatial decomposition
into 3 partitions is shown where the list of bodies is separated into 3 sections, each
containing 7 bodies. Notice that
a partition may span a spatial discontinuity. This is a result of the
Morton ordering; a solution would be to use the more complicated Peano-Hilbert ordering.

Algorithm <ref> presents the parallel Barnes-Hut algorithm which makes use of
spatial decomposition. At each simulation step, all bodies
must be gathered to the root process so that
their spatial keys can be computed, the bodies sorted,
and then the bodies scattered again to each process (Lines 4–8).
The sort followed by the scatter implies the spatial decomposition, where
the scatter evenly distributes the sorted bodies to each process.
The algorithm proceeds then as Algorithm <ref>, constructing
and sharing a global octree by a reduce and broadcast.

Parallel Barnes-Hut V3($id$, $s$, $\Delta t$, $t_{end}$, $N$, $P_{all}$)$P_{all}$ is the list of size $N$ of bodies to simulate from 0 to $t_{end}$ going by $\Delta t$ time-steps. $id$ is the executing process's id, $s$ is the number of processes. Assume $s \mid N$.
[1]
$N' \gets N / s$
$P \gets$ Scatter $P_{all}$, sending to process $id$ $P_{all}[id\cdot N', \ldots, id\cdot N' + N' - 1]$ Get Local Bodies
$t$ from 0 to $t_{end}$ by $\Delta t$
$P_{all} \gets$ Gather $P$ to process 0. Spatial Decomposition
id = 0
$keys \gets $ ConstructSpatialKeys($P_{all}$).
Sort $P_{all}$ by $keys$.
$P \gets$ Scatter $P_{all}$, sending to process $id$ $P_{all}[id\cdot N', \ldots, id\cdot N' + N' - 1]$
$T' \gets$ build octree from $P$ Parallel Tree Build
$T \gets $ Reduce, via Octree Merge, $T'$ to process 0.
Broadcast $T$ from process 0 to all others.

each particle $p$ in $P$ Compute Forces; unchanged from Algorithm <ref>
each particle $p$ in $P$ Position Update
Update position and velocity of $p$ using its current acceleration.

One crucial side effect of the spatial decomposition is its effect
on the octree merge. Recall that, in the worst case, each
merge requires $\mathcal{O}(n\log n)$ steps to merge two
trees which each have $n$ uniformly distributed bodies.
As noted in <cit.>, since
each process now contains a distinct spatial partition,
the leaf nodes in each local octree do not overlap.
Each local octree thus only shares a few high-level nodes
with other local octrees. A merge
between two octrees then occurs in time $\mathcal{O}(\log n)$,
since only one traversal is needed to find where to place
entire subtrees.
The entire reduce then has a span of $\mathcal{O}(\log N)$, meanwhile
the span of the local tree construction remains as $\mathcal{O}(\nicefrac{N}{s}\log(\nicefrac{N}{s}))$
The dominant term is the latter, and thus the span
of the entire parallel tree construction becomes $\mathcal{O}(\nicefrac{N}{s}\log(\nicefrac{N}{s}))$.

Notice further that, for a sufficiently small time-step $\Delta t$,
the bodies will only move very slightly between time-steps.
Thus, a sorted list of bodies will, after the integration step,
still be mostly sorted.
We thus implement an insertion sort to sort the particles and keys.
Sorting therefore only requires a few swaps and time $\mathcal{O}(N)$.
Yet, it still requires high amounts of communication, and sub-optimal scalability,
where the particles are gathered to the root process for sorting.
We tackle this problem in Section <ref>.
However, we first explore a more crucial part of the algorithm: load-balancing
the force calculations.

§.§ Dynamic Load-Balancing

The dominant operation in an $N$-body simulation is the force calculation.
In the direct method it takes time $\mathcal{O}(N^2)$ and in a
hierarchical method $\mathcal{O}(N \log N)$. However,
not every interaction is created equal when computing forces acting on a particle.
A particle in a dense region will have more particle-particle interactions,
thus approaching an $N$ number of interactions.
A very distant particle in a sparse region may have only particle-cell interactions,
thus requiring only $\mathcal{O}(1)$
interactions, depending on the value of the MAC; see Section <ref>.
Moreover, the big-$\mathcal{O}$ notation can be misleading.
The force calculation step is by far the most dominant part of the $N$-body
simulation—despite it and tree building both taking $\mathcal{O}(N\log N)$
time in serial—taking between 70% and 90% of the time of each simulation time-step <cit.>.

Therefore, our spatial decomposition scheme, which assigns
an equal $\nicefrac{N}{s}$ number of bodies to each of the
$s$ processes, is likely insufficient.

The resolution is that of the aptly named costzones,
to determine zones—spatial partitions—of equal cost.
The exact same operation is performed in the hashed octree method,
but is not named nor emphasized.
The method of determining costzones is simple.
Begin with a work estimate $w_i$ for each body as the number of
interactions computed for that body in the previous time-step
(in the very first time-step, assume $w_i=N$).
Then, each partition (costzone) should be responsible for
$w = \left(\sum_{i=0}^N w_i\right) / s$ amount of work.
Determining costzones proceeds by iterating over
the spatially-linear list of bodies
(the list of bodies after sorting by the Morton ordering)
and partitioning it into sets whose work estimate is at least $w$.
Algorithm <ref> shows this iteration.

Costzones($W$, $s$)$W$ is a list of size $N$ of work estimates
such that $W[i]$ is the work estimate of the $i$th body in the Morton (or other) ordering; $s$ is the number of processes.Returns a list $Z$ of size $s+1$, returning the index of the first body in process $k$'s domain partition as $Z[k]$.
[1]
$w_{tot} \gets 0$
$i$ from 0 to $N-1$
$w_{tot} \gets w_{tot} + W[i]$
$w_{targ} \gets w_{tot} / s$
$k \gets 0$; $w \gets 0$
$j \gets 0$; $Z \gets [\ ]$
$i$ from 0 to $N-1$
$w \gets w + W[i]$
$w \geq w_{targ}$
$Z[k] \gets j$; $j \gets i+1$
$k \gets k + 1$; $w \gets 0$
$Z[s-1] \gets j$ Commit final partition
$Z[s] \gets N$ Used to compute end of partition for $id$ $s-1$

Algorithm <ref> adds
the Costzones load balancing
to Algorithm <ref>.
It begins by gathering the work estimates of each particle
to the root process (Line 6).
This list of work estimates $W_{all}$ is then sorted
along with $P_{all}$ as parallel arrays (Lines 8–9).
Costzones are computed using $W_{all}$ and then
$P_{all}$ is scattered to all processes
based on the costzones partitioning (Lines 10–11).
Lastly, the tree traversal is modified
to accumulate the number of interactions of
each particle in the array $W$ (Lines 16, 22, 28).
A subtle optimization for the work estimate is to
weight the work of a particle-cell interaction
twice as much as a particle-particle interaction.
This is because of the multipole approximation incurred by a
particle-cell interaction; computing acceleration
from a particle-cell interaction involves twice as many floating point
operations as a simple particle-particle interaction (see Section <ref>).

Parallel Barnes-Hut V4($id$, $s$, $\Delta t$, $t_{end}$, $N$, $P_{all}$)$P_{all}$ is the list of size $N$ of bodies to simulate from 0 to $t_{end}$ going by $\Delta t$ time-steps. $id$ is the executing process's id, $s$ is the number of processes. Assume $s \mid N$.
[1]
$N' \gets N / s$
$P \gets$ Scatter $P_{all}$, sending to process $id$ $P_{all}[id\cdot N', \ldots, id\cdot N' + N' - 1]$ Get Local Bodies
$W \gets [N,\dots,N]$ $W$ of size $N'$
$t$ from 0 to $t_{end}$ by $\Delta t$
$P_{all} \gets$ Gather $P$ to process 0. Spatial Decomposition
$W_{all} \gets$ Gather $W$ to process 0.
id = 0
$keys \gets $ ConstructSpatialKeys($P_{all}$).
Sort $P_{all}$ and $W_{all}$ by $keys$.
$Z \gets$ Costzones($W_{all}$, $s$)
$P \gets$ Scatter $P_{all}$ sending to process $id$ $P_{all}[Z[id], \ldots, Z[id+1]-1]$
$T' \gets$ build octree from $P$ Parallel Tree Build
$T \gets $ Reduce, via Octree Merge, $T'$ to process 0.
Broadcast $T$ from process 0 to all others.
$i$ from 0 to $|P|-1$ Compute Forces
$p \gets P[i]$; $W[i] \gets 0$
Set the acceleration of $p$ to $\vec{0}$.
$Node$ $\gets$ the root node of $T$

$\quad\ $Compute Cell Force:
each child cell $C$ of $Node$
$C$ contains 1 particle $C$ is a leaf
Add to $p$'s acceleration the contribution from $C$'s particle. particle-particle
$W[i] \gets W[i] + 1$
$\ell$ $\gets$ $C$'s side length
$d$ $\gets$ distance between $C$'s c.o.m. and $p$'s position.
${\ell} / {d} < \theta$
Add to $p$'s acceleration the contribution from $C$. particle-cell
$W[i] \gets W[i] + 2$
Compute Cell Force with $Node \gets C$ recursive call
each particle $p$ in $P$ Position Update
Update position and velocity of $p$ using its current acceleration.

The load balancing of costzones more evenly distributes the amount of work
performed by each process during the force calculation step.
This comes at the expense of an uneven number of bodies
assigned to each process for the octree construction
and integration steps. Experiments suggest the number of
bodies assigned to any one process is roughly twice as many as another.
In long-term evolution, clustering may lead to the maximum number of bodies
on one process being 10 times as many as the one with the fewest <cit.>.
Nonetheless, since the force calculation is such a dominant part of the algorithm,
load-balancing the force calculation results in better overall performance
despite the loss of load-balancing in the other steps.

§.§ Parallel and Distributed Sorting

Our next step towards improved scalability is to avoid
the gather of bodies and work estimates to the root process.
This gather occurs so that the spatial decomposition via sorting,
and load-balancing via costzones, can occur.
This gather also incurs high amounts of communication cost,
particularly as the number of processes increases.
We propose here a solution which provides a distributed algorithm for both
the sorting and load-balancing steps so as to avoid excess communication.
We note that methods detailed in this section did not appear
in the original costzones method <cit.>, but are
a result of applying this method to a distributed system.

We begin with the sorting algorithm. One could implement
a generalized and highly scalable distributed sort, e.g.
methods based on sample sort; see <cit.>.
However, we have already made the observation that
bodies remain mostly sorted between time-steps.
Hence we implement a modification of bubble sort to a distributed system.
Simply stated, elements out of place
bubble up and are exchanged with the process's right neighbour, or sink down and
are exchanged with the process's left neighbour.

Assume a linear topology of the processes.
Neighbour processes exchange fixed amounts of data,
both neighbours sort this exchanged data,
and then each keeps their respective halves.
For example, process with $id$ 0 would send its
$k$ largest elements to $id$ 1, and $id$ 1 would send
its $k$ smallest elements. These $2k$ elements are locally sorted
by both processes and $id$ 0 keeps the smallest $k$, while $id$ 1 keeps the largest $k$.
Finally, the new data is sorted into the local data collection on each process
to update which elements are now the $k$ largest or $k$ smallest.
These exchanges repeat between all sets of neighbours until an iteration
occurs with no changes on any process.

This distributed sort is shown in Algorithm <ref>.
Local sorting is implemented
by insertion sort.
Again, since bodies are already mostly sorted,
this bubble sort generalization should
behave as in the best case, resulting
in only $\mathcal{O}(1)$ rounds of
exchanges and $\mathcal{O}(\nicefrac{N}{s})$ total work.

With the addition of a distributed sort, we would like to continue to avoid gathering the bodies and work estimates
to the root process for the computation of the costzones.
We thus propose a modified costzones approach which load-balances
without a gather, requiring only exchanges between neighbouring processes.
This is a simple greedy algorithm, shown in Algorithm <ref>.
Each process begins by determining the work target via a reduce and
broadcast (Lines 1–6). Then, iterating from process 0 to process $s-2$, each process
“pushes” (Lines 14–16) or “pulls”  (Lines 17–19) data to or from its right neighbour
to achieve it's own load balance.
Process $s-1$ is left with whatever remains.
This redistribution and data transfer must be careful to
not change the ordering. Notice the ordering of the concatenations
upon receiving data (Lines 19 and 24).

Having modified the spatial decomposition and load-balancing
to be in a distributed manner,
Algorithm <ref> presents
this finalized costzones method, with reduced communications
and still parallelized tree construction, tree traversal, and integration steps.
This presents the “final” version of a costzones-based Barnes-Hut simulation.
The next improvement to make on top of Algorithm <ref>
is to avoid the octree merge and broadcast. Indeed, for the highest scalability
on the largest simulations, we should avoid having the entire global tree
present on one process. This is precisely the intent of the hashed-octree method.

SE[DOWHILE]DodoWhile[1] #1

Parallel Sort($V$, $id$, $s$, $k=64$)$V$ is a list of size $n$ of values to sort, in place. $id$ is the executing process's id, $s$ is the number of processes. $k$ controls the number of values to exchange between neighbours at each step. The result is any element of $V$ on process $id$ is less than any element of $V$ on process $id+1$, and $V$ is locally sorted. This assumes $n$ > $2k$. An experimentally determined default is $k=64$.
[1]
Sort $V$.
$Left \gets V[0,...,k-1] $ concat $V[0,...,k-1]$
$Right \gets V[n-k,...,n-1] $ concat $V[n-k,...,n-1]$
$id - 1 > 0$
Send $Left[0,...,k-1]$ to $id - 1$.
$Left[0,...,k-1] \gets k$ values from $id - 1$.
Sort $Left$.
$V[0,...,k-1] \gets Left[k,...,2k-1]$ Keep upper half
$id + 1 < s$
Send $Right[k,...,2k-1]$ to $id + 1$.
$Right[k,...,2k-1] \gets k$ values from $id + 1$.
Sort $Right$.
$V[n-k,...,n-1] \gets Right[0,...,k-1]$ Keep lower half
Sort $V$.
$c' \gets V$ modified during sort?
$c \gets $Reduce, via logical OR, $c'$ to process 0.
$c \gets $Broadcast $c$ from process 0 to all others.
$c \neq 0$

ParallelCostzones($P$, $W$, $id$, $s$)$P$ and $W$ are lists of size $n$ of particles and work estimates such that $W[i]$ is the work estimate for $P[i]$. $id$ is the executing process's id, $s$ is the number of processes. Distributes work between neighbours to achieve costzones-like load-balancing.
[1]
$w_{tot} \gets 0$
$i$ from 0 to $N-1$
$w_{tot} \gets w_{tot} + W[i]$
$w_{tot} \gets $Reduce, via addition, $w_{tot}$ to process 0.
$w_{tot} \gets $Broadcast $w_{tot}$ from process 0 to all others.
$w_{targ} \gets w_{tot} / s$
$k$ from 0 to $s-2$
$id = k$
$w \gets 0$; $j \gets 0$
$i$ from 0 to $n-1$
$w \gets w + W[i]$
$w \geq w_{targ}$ and $j = 0$ $j \gets i$ Send to process $id + 1$ $(w - w_{targ})$
$w > w_{targ}$ Push excess work to right neighbour
$P,P' \gets P[0,...,j], P[j+1,...,n]$; $W,W' \gets W[0,...,j], W[j+1,...,n]$
Send to process $id + 1$ $(P',W')$
$w < w_{targ}$ Pull work from right neighbour
$P',W' \gets $ Receive data from process $id + 1$
$P \gets P $ concat $P'$; $W \gets W $ concat $W'$
$id = k + 1$
$w' \gets $ Receive work difference from $id - 1$
$w' > 0$ Pull excess work from left neighbour
$P',W' \gets $ Receive data from process $id - 1$
$P \gets P' $ concat $P$; $W \gets W' $ concat $W$
$w' < 0$ Push work to left neighbour
$w \gets 0$; $j \gets n$
$i$ from 0 to $n-1$
$w \gets w + W[i]$
$w \geq |w'|$ $j \gets i$; break; $P',P \gets P[0,...,j], P[j+1,...,n]$; $W',W \gets W[0,...,j], W[j+1,...,n]$
Send to process $id - 1$ $(P',W')$
$n \gets |P|$
Synchronize all processes Ensure process $k$ has finished before continuing

Parallel Barnes-Hut V5($id$, $s$, $\Delta t$, $t_{end}$, $N$, $P_{all}$)$P_{all}$ is the list of size $N$ of bodies to simulate from 0 to $t_{end}$ going by $\Delta t$ time-steps. $id$ is the executing process's id, $s$ is the number of processes. Assume $s \mid N$.
[1]
$N' \gets N / s$
$P \gets$ Scatter $P_{all}$, sending to process $id$ $P_{all}[id\cdot N', \ldots, id\cdot N' + N' - 1]$ Get Local Bodies
$W \gets [N,\dots,N]$ $W$ of size $N'$
$t$ from 0 to $t_{end}$ by $\Delta t$
$keys \gets $ ConstructSpatialKeys($P$)
ParallelSort($\{P,W,keys\}, id, s$)Spatial Decomp.; sorting particle-work-key tuples by key
ParallelCostzones($P, W, id, s$) Load balance
$T' \gets$ build octree from $P$ Parallel Tree Build
$T \gets $ Reduce, via Octree Merge, $T'$ to process 0.
Broadcast $T$ from process 0 to all others.
$i$ from 0 to $|P|-1$ Compute Forces; unchanged from Algorithm <ref>
each particle $p$ in $P$ Position Update
Update position and velocity of $p$ using its current acceleration.

§.§ A Distributed Octree Implementation

The hashed octree method is concerned with keeping the
scalability of the simulation high and avoiding
excess communication caused by broadcasting the global octree.
Rather, this method sees that
each process dynamically requests the data it needs
from other processes during the tree traversal of its local tree.

Such a distributed data structure requires a consistent addressing
scheme such that an octree node on one process
can be translated into a node on a different process.
Moreover, a data structure which allows direct access to a particular node
would be beneficial. Recall the basic implementation of
an octree based on pointers and linked lists from Section <ref>.
Accessing this structure involves traversing pointers and moving from root to leaf—
accessing a leaf thus requires time $\mathcal{O}(\log N)$.
Both aspects do not fulfill the needs of a distributed data structure.
However, the linear octree or hashed octree (originally presented in <cit.>) does serve this purpose.

In a hashed octree, each cell is given a unique key. This key
directly serves as the index—modulo the table size—into a hash table.
The hash table is a standard implementation using chaining and linked lists
for conflict resolution; see, e.g., <cit.>
or <cit.>.
A child cell's key is derived from its parent's key plus
its linear order in the parent's list of children.
These keys will be represented in octal notation for simplicity.
The root has key 1, its 8 children have keys: 10, 11, 12, 13, 14, 15, 16, 17.
This procedure continues recursively; for example the children of node 12 are:
120, 121, 122, 123, 124, 125, 126, 127.
Notice that this scheme allows for determining the key of a parent
or child of a particular node through simple bit operations.
The parent's key is obtained by a 3-bit right shift;
the child's key is obtained by a 3-bit left shift followed
by a bit-wise or with the child's index (0–7).

Any ordering may be used when assigning the indices 0–7 to the children
of a particular node, as long as the same ordering is used at each level.
A natural ordering in our context is of course the
Morton ordering of Section <ref>.
This Morton ordering of a cell's position to a spatial key
serves the purpose of a traditional hash function,
distributing cells of the octree evenly across the hash table.
The hash function is then just a mod function over the hash table's size.
Moreover, the spatial keys computed for each body using
the Morton ordering directly correspond
to the key of the octree node—or the node's child, grandchild, etc.—in which they are contained.
That is to say, the body's key corresponds to a leaf node's key
if all levels of spatial division are used.
It is often not needed to actually divide
the octree into that many levels.

Fig. <ref> shows the key assignment via Morton's ordering
for up to three levels
of an octree and an encoding of that octree as a hash table.
The nodes of the pointer-based octree can be repurposed
to be nodes of the hash table linked list with a simple transformation:
$(i)$ the list of pointers to children is replaced by a single pointer
to the next cell for collision resolution via chaining; and $(ii)$
the cell's key is included in the payload.

[->,>=Stealth[width=1mm,length=1mm]] (-1.25, 0, 3) – (-1.25, 0, 3.75);
[->,>=Stealth[width=1mm,length=1mm]] (-1.25, 0, 3) – (-1.25, 0.55, 3);
[->,>=Stealth[width=1mm,length=1mm]] (-1.25, 0, 3) – (-0.75, 0, 3);
[font=] at (-1.1, 0.45, 3) $y$;
[font=] at (-1.1, 0, 3.65) $z$;
[font=] at (-0.91, 0.09, 3) $x$;

(O) at (0,0,0);
(A) at (0,,0);
(B) at (0,,);
(C) at (0,0,);
(D) at (,0,0);
(E) at (,,0);
(F) at (,,);
(G) at (,0,);

[black, semithick, line cap=round] (O) – (C) – (G) – (D) – cycle;

[black, semithick, line cap=round] (A) – (B) – (C) – (O);
[black, semithick, opacity=0.2] (A) – (O);
[black, semithick, line cap=round] (D) – (E) – (F) – (G) – cycle;
[black, semithick, line cap=round] (C) – (B) – (F) – (G) – cycle;
[black, semithick, line cap=round] (A) – (B) – (F) – (E) – cycle;

[blue] (2, 0, 0) – (2, 4, 0);
[blue] (0, 2, 0) – (4, 2, 0);

[blue] (4, 2, 0) – (4, 2, 4);
[blue] (4, 0, 2) – (4, 4, 2);

[blue] (2, 0, 0) – (2, 0, 4);
[blue] (0, 0, 2) – (4, 0, 2);

[blue] (0, 2, 0) – (0, 2, 4);
[blue] (0, 0, 2) – (0, 4, 2);

[blue] (2, 4, 0) – (2, 4, 4);
[blue] (0, 4, 2) – (4, 4, 2);

[blue] (2, 0, 4) – (2, 4, 4);
[blue] (0, 2, 4) – (4, 2, 4);

[blue] (0, 2, 2) – (4, 2, 2);
[blue] (2, 2, 0) – (2, 2, 4);
[blue] (2, 0, 2) – (2, 4, 2);

[fill=blue!50,fill opacity=0.08, draw opacity=0] (0,4,0) – (4,4,0) – (4,0,0) – (4,0,4) – (0,0,4) – (0,4,4) – (0,4,0);
[fill=blue!50,fill opacity=0.08,draw opacity=0] (0,2,0) – (4,2,0) – (4,2,4) – (0,2,4) – cycle;
[fill=blue!50,fill opacity=0.08,draw opacity=0] (0,0,2) – (4,0,2) – (4,4,2) – (0,4,2) – cycle;
[fill=blue!50,fill opacity=0.08,draw opacity=0] (2,0,0) – (2,4,0) – (2,4,4) – (2,0,4) – cycle;

[blue,fill=green!50, opacity=0.1] (2,2,2) – (2,4,2) – (2,4,4) – (2,2,4) – cycle;
[blue,fill=green!50,opacity=0.1] (2,2,2) – (4,2,2) – (4,4,2) – (2,4,2) – cycle;
[blue,fill=green!50,opacity=0.1] (2,2,2) – (4,2,2) – (4,2,4) – (2,2,4) – cycle;

[fill=green!50,fill opacity=0.2,draw opacity=0] (2,3,2) – (4,3,2) – (4,3,4) – (2,3,4) – cycle;
[fill=green!50,fill opacity=0.2,draw opacity=0] (2,2,3) – (4,2,3) – (4,4,3) – (2,4,3) – cycle;
[fill=green!50,fill opacity=0.2,draw opacity=0] (3,2,2) – (3,4,2) – (3,4,4) – (3,2,4) – cycle;

[green!80!black] (3, 2, 2) – (3, 4, 2);
[green!80!black] (2, 3, 2) – (4, 3, 2);

[green!80!black] (4, 3, 2) – (4, 3, 4);
[green!80!black] (4, 2, 3) – (4, 4, 3);

[green!80!black] (3, 2, 2) – (3, 2, 4);
[green!80!black] (2, 2, 3) – (4, 2, 3);

[green!80!black] (2, 3, 2) – (2, 3, 4);
[green!80!black] (2, 2, 3) – (2, 4, 3);

[green!80!black] (3, 4, 2) – (3, 4, 4);
[green!80!black] (2, 4, 3) – (4, 4, 3);

[green!80!black] (3, 2, 4) – (3, 4, 4);
[green!80!black] (2, 3, 4) – (4, 3, 4);

[green!80!black] (2, 3, 3) – (4, 3, 3);
[green!80!black] (3, 3, 2) – (3, 3, 4);
[green!80!black] (3, 2, 3) – (3, 4, 3);

[] at (0.7,0.2,1.8) 10;
[] at (0.7,0.2,3.8) 11;
[] at (0.7,4.2,1.8) 12;
[] at (0.7,4.2,3.8) 13;
[] at (2.7,0.2,1.8) 14;
[] at (2.7,0.2,3.8) 15;
[] at (2.7,4.2,1.8) 16;

[font=] at (2.4, 2.1, 2.9) 170;
[font=] at (2.4, 2.1, 3.9) 171;
[font=] at (2.42, 4.1, 2.9) 172;
[font=] at (2.42, 4.1, 3.9) 173;
[font=] at (3.37, 2.1, 2.9) 174;
[font=] at (3.4, 2.1, 3.9) 175;
[font=] at (3.42,4.1, 2.9) 176;
[font=] at (3.42, 4.1, 3.9) 177;

[count=from 0] in 0,...,-7
[font=] at (-0.35,-0.5) ;
[thick] (0,0) – (1.25,0) – (1.25,-8) – (0,-8) – cycle;
in 0,...,-7
(0,) – (1.25,);
in 0,...,-7
(0.4,-0.3) rectangle (0.8,-0.7);
in 0,...,-7 [->, >=Stealth] (0.625,-0.5)   – (2.0,-0.5);
in 0,...,-7 [green!60!black] (2.0,-0.1)  rectangle (3.75,-0.9);
[count=from 0] in 0,...,-7 [font=] at (2.45,-0.5) 17;
in 0,...,-7 (3.0,-0.1) – (3.0,-0.9);
in 0,...,-7
(3.15,-0.3) rectangle (3.55,-0.7);
in 0,...,-7 [->, >=Stealth] (3.35,-0.5)   – (4.5,-0.5);

in 0,...,-7 [blue!50!black] (4.5,-0.1)  rectangle (6.25,-0.9);
[count=from 0] in 0,...,-7 [font=] at (4.99,-0.5) 1;
in 0,...,-7 (5.5,-0.1) – (5.5,-0.9);

in 0,...,-7
(5.65,-0.3) rectangle (6.05,-0.7);

in 0,...,-7 (5.85,-0.5)   – (6.8,-0.5);

at (6.8,-0.5) $\times$;
in -2,...,-7 at (6.8,-0.5) $\times$;

[->,>=Stealth] (5.85,-1.5) – (7.0,-1.5);
[black] (7.0,-1-0.1)  rectangle (8.75,-1-0.9);
[font=] at (7.45,-1-0.5) 1;
(8.0,-1-0.1) – (8.0,-1-0.9);
(8.15,-1-0.3) rectangle (8.55,-1-0.7);
(8.35,-1.5) – (9.4,-1.5);
at (9.4,-1.5) $\times$;

An octree with cells labelled by their hashed octree key, written in octal. The root node's key of 1 is not shown. Keys are created
from a Morton ordering with $x > y > z$. That octree's encoding as a hash table—an array of linked lists—with
8 entries is shown to the right. Note that the payload of each cell is not shown.

A hashed octree implementation allows for consistent indexing of nodes
in the tree, rather than reliance on pointer traversal. Further,
it allows access to any node in the tree in time $\mathcal{O}(1)$.
This scheme allows the $N$-body simulation
to avoid the creation of a global tree
and for each process to dynamically and efficiently request remote data
from other processes as needed during the tree traversal.

The last requirement is for each process to share the boundary
of its spatial domain so that each process knows from where
to request data during during the tree traversal.
This sharing proceeds in three steps.

* Assume a linear topology of processes, just as in distributed sorting and costzones. In the first step, each process shares the first body in its spatial domain with its left neighbour and the last body in its spatial domain with its right neighbour, if such a neighbour exists. The “neighbour bodies” (or single body for process 0 and $s-1$) are inserted into the local octree, splitting as needed so that each leaf node contains exactly one body. Due to the spatial ordering of bodies, this insertion will cause the required splitting of leaf nodes so that the union of all leaf nodes across the $s$ local octrees is exactly the leaf nodes that would be created if a global octree were to be constructed. That is to say, all leaf nodes (excluding the leaves containing neighbour bodes) are distinct across all processes and a consistent spatial partitioning has been created.

* In the second step, each process computes and broadcasts its so-called branch nodes. These are the set of nodes of the octree which are closest to the root, but which do not contain a neighbour body in their domain.
They represent the entire spatial domain of a process at the coarsest level possible. Branch nodes can be found easily by a depth-first traversal of the the local octree: if a cell's domain does not contain a neighbour body, but the cell's parent's domain does, then that cell is a branch node.

* In the third step, a process receives sets of branch nodes as they are broadcast. For each set, those branch nodes are inserted into the local octree and marked with their source process. This source id will allow for dynamically requesting data during tree traversal. If a branch node is received for which its parent is not in the local octree, the parent, grandparent, etc. nodes are “filled” in to the local octree as well.

Upon inserting all the branch nodes, each process has a local octree which
is globally and mutually consistent with the other local trees.
Fig. <ref> shows this neighbour exchange
and resulting branch nodes for three processes.
After sharing and inserting these branch nodes,
each local tree contains the required information
to request nodes below that of the received branch nodes.
Such nodes must be requested during tree traversal when nodes
fail the multiple acceptance criterion and must be “opened” (see Section <ref>).

DistributeBranchNodes($P$, $T$, $id$, $s$)$P$ is the list of local bodies and $T$ is the octree containing those bodies. $id$ is the executing process's id, $s$ is the number of processes.
[1]
$id - 1 \geq 0$
Send $P[0]$ to process $id -1$
$p \gets $ Receive neighbour body from $id -1$
Insert $p$ into $T$, splitting as required.
$id + 1 < s$
Send $P[n]$ to process $id + 1$
$p \gets $ Receive neighbour body from $id + 1$
Insert $p$ into $T$, splitting as required.
$B \gets$ Branch nodes of $T$, the coarsest noes which do not contain a neighbour body.
$i$ from 0 to $s-1$
Broadcast $B$ from process $i$ to all others.
$id \neq i$
$B' \gets$ Receive branch nodes from $i$
each node $b$ in $B'$ insert $b$ into $T$ and fill in parents

[step=2cm, thin] (0,4) grid +(4,4)
[step=1cm, thin] (2,4) grid +(2,2)

[step=0.5cm, thin] (2, 4) grid +(1,1);
[step=4cm, very thick, line cap=round] (0,0) grid (8,8);

[fill] (0.88, 1.77) circle (0.06cm);
[fill] (0.43, 4.65) circle (0.06cm);
[fill] (0.77, 7.23) circle (0.06cm);
[fill] (2.11, 4.13) circle (0.06cm);
[fill] (2.41, 4.81) circle (0.06cm);
[fill] (2.75, 4.25) circle (0.06cm);
[fill] (2.70, 4.60) circle (0.06cm);

[fill=blue!80, fill opacity=0.1] (0.1,0.1) – (0.1,3.9) – (3.9,3.9) – (3.9, 0.1) – cycle;
[fill=blue!80, fill opacity=0.1] (0.1,4.1) – (1.9,4.1) – (1.9,5.9) – (0.1,5.9) – cycle;
[fill=blue!80, fill opacity=0.1] (0.1,6.1) – (1.9,6.1) – (1.9,7.9) – (0.1,7.9) – cycle;
[fill=blue!80, fill opacity=0.1] (2.1,4.1) – (2.9,4.1) – (2.9,4.9) – (2.1,4.9) – cycle;

[step=2cm, thin] (0,4) grid +(4,4)
[step=1cm, thin] (2,4) grid +(2,2)

[step=0.5cm, thin] (2, 4) grid +(1,1);
[step=4cm, very thick, line cap=round] (0,0) grid (8,8);

[fill] (0.88, 1.77) circle (0.06cm);
[fill] (0.43, 4.65) circle (0.06cm);
[fill] (0.77, 7.23) circle (0.06cm);
[fill] (2.14, 4.14) circle (0.06cm);
[fill] (2.41, 4.81) circle (0.06cm);
[fill] (2.75, 4.25) circle (0.06cm);
[fill] (2.70, 4.60) circle (0.06cm);

[fill=orange!30!yellow,draw=orange!30!yellow!40!black] (3.60, 4.44) circle (0.07cm);

[step=2cm, thin] (0,4) grid +(4,4)
[step=1cm, thin] (2,4) grid +(2,2)

[step=2cm, thin] (4,0) grid + (4,4)
[step=0.5cm, thin] (4, 2) grid +(1,1)
[step=0.5cm, thin] (4, 3) grid +(1,1)
[step=1cm, thin] (5, 2) grid +(1,1);

[step=4cm, very thick, line cap=round] (0,0) grid (8,8);

[fill] (3.60, 4.44) circle (0.06cm);
[fill] (3.44, 5.77) circle (0.06cm);

[fill] (4.20, 2.78) circle (0.06cm);
[fill] (4.65, 2.42) circle (0.06cm);
[fill] (4.33, 3.27) circle (0.06cm);
[fill] (4.26, 3.88) circle (0.06cm);
[fill] (5.18, 2.77) circle (0.06cm);

[fill=orange!90!yellow, fill opacity=0.1] (3.1,5.1) – (3.9,5.1) – (3.9,5.9) – (3.1,5.9) – cycle;
[fill=orange!90!yellow, fill opacity=0.1] (3.1,4.1) – (3.9,4.1) – (3.9,4.9) – (3.1,4.9) – cycle;
[fill=orange!90!yellow, fill opacity=0.1] (4.1,3.1) – (4.9,3.1) – (4.9,3.92) – (4.1,3.92) – cycle;
[fill=orange!90!yellow, fill opacity=0.1] (4.1,2.1) – (4.9,2.1) – (4.9,2.9) – (4.1,2.9) – cycle;
[fill=orange!90!yellow, fill opacity=0.1] (5.1,2.1) – (5.9,2.1) – (5.9,2.9) – (5.1,2.9) – cycle;

[step=2cm, thin] (0,4) grid +(4,4)
[step=1cm, thin] (2,4) grid +(2,2)

[step=2cm, thin] (4,0) grid + (4,4)
[step=0.5cm, thin] (4, 2) grid +(1,1)
[step=0.5cm, thin] (4, 3) grid +(1,1)
[step=1cm, thin] (5, 2) grid +(1,1);

[step=4cm, very thick, line cap=round] (0,0) grid (8,8);

[fill=blue!60,draw=blue!80] (2.70, 4.60) circle (0.07cm);

[fill] (3.60, 4.44) circle (0.06cm);
[fill] (3.44, 5.77) circle (0.06cm);

[fill] (4.20, 2.78) circle (0.06cm);
[fill] (4.65, 2.42) circle (0.06cm);
[fill] (4.33, 3.27) circle (0.06cm);
[fill] (4.26, 3.85) circle (0.06cm);
[fill] (5.18, 2.77) circle (0.07cm);
[fill=red!80!blue,fill opacity=0.5,draw=red!20!black,very thin] (5.80, 3.13) circle (0.07cm);

[step=2cm, thin] (4,4) grid +(4,4)

[step=1cm, thin] (4,4) grid +(2,2)
[step=2cm, thin] (4,0) grid + (4,4);

[step=4cm, very thick, line cap=round] (0,0) grid (8,8);

[fill] (5.80, 3.13) circle (0.06cm);

[fill] (7.22, 0.8) circle (0.06cm);
[fill] (6.77, 3.2) circle (0.06cm);

[fill] (4.64, 4.76) circle (0.06cm);
[fill] (5.72, 5.31) circle (0.06cm);
[fill] (6.22, 4.81) circle (0.06cm);
[fill] (6.4, 6.81) circle (0.06cm);

[fill=red, fill opacity=0.1] (4.1, 4.1) – (4.1,7.9) – (7.9, 7.9) – (7.9,4.1) – cycle;
[fill=red, fill opacity=0.1] (6.1, 0.1) – (7.9,0.1) – (7.9, 1.9) – (6.1,1.9) – cycle;
[fill=red, fill opacity=0.1] (6.1, 2.1) – (7.9,2.1) – (7.9, 3.9) – (6.1,3.9) – cycle;
[fill=red, fill opacity=0.1] (5.07, 3.07) – (5.93,3.07) – (5.93, 3.9) – (5.07,3.9) – cycle;

[step=2cm, thin] (4,4) grid +(4,4)

[step=1cm, thin] (4,4) grid +(2,2)
[step=2cm, thin] (4,0) grid + (4,4)

[step=1cm, thin] (4, 2) grid +(2,2);

[step=4cm, very thick, line cap=round] (0,0) grid (8,8);

[fill=orange!30!yellow,draw=orange!30!yellow!40!black] (5.18, 2.77) circle (0.07cm);
[fill] (5.80, 3.13) circle (0.06cm);

[fill] (7.22, 0.8) circle (0.06cm);
[fill] (6.77, 3.2) circle (0.06cm);

[fill] (4.64, 4.76) circle (0.06cm);
[fill] (5.72, 5.31) circle (0.06cm);
[fill] (6.22, 4.81) circle (0.06cm);
[fill] (6.4, 6.81) circle (0.06cm);

Three local quadtrees, (a), (b), (c), before sharing neighbour bodies, and their respective updated quadtrees (d), (e), (f), after sharing neighbour bodies. The neighbour bodies are colored corresponding to their owning process. Notice in (f) that the addition of the neighbour body causes a leaf node to split, creating consistent leaf nodes
between (d) and (f). The branch nodes of each process are colored in (b), (d), (f). The branch nodes collectively cover every body in the simulation domain.

With the hashed octree data structure and the required
branch nodes fully specified, we can now take advantage of this
distributed data structure to remove the final reduce and broadcast
(apart from a few branch nodes) of our parallel Barnes-Hut algorithm.
The span of the simulation is now minimized.
Algorithm <ref> presents this modified Barnes-Hut method
which–implicitly–makes use of a hashed octree encoding
of the octree.
First, the reduce-broadcast of the global tree is replaced
with a call to DistributeBranchNodes (Line 9).
Second, the tree traversal is slightly modified to
check if a node's children are local or not.
If not, they are requested from their source process, as determined
by their containing branch node (Lines 14–16).
Finally, all processes synchronize before moving to the integration step
to ensure all requests from all processes have been fulfilled (Lines 19–20).

In this configuration the requesting process stalls until its request has been fulfilled.
Therefore, before a request is actually sent, the process fulfills any
and all requests made to it. This avoids the situation
where two processes are waiting on each other to fulfill a request.
In the next section we will improve this
to avoid stalling and instead perform asynchronous communication.

Parallel Barnes-Hut V6($id$, $s$, $\Delta t$, $t_{end}$, $N$, $P_{all}$)$P_{all}$ is the list of size $N$ of bodies to simulate from 0 to $t_{end}$ going by $\Delta t$ time-steps. $id$ is the executing process's id, $s$ is the number of processes. Assume $s \mid N$.
[1]
$N' \gets N / s$
$P \gets$ Scatter $P_{all}$, sending to process $id$ $P_{all}[id\cdot N', \ldots, id\cdot N' + N' - 1]$ Get Local Bodies
$W \gets [N,\dots,N]$ $W$ of size $N'$
$t$ from 0 to $t_{end}$ by $\Delta t$
$keys \gets $ ConstructSpatialKeys($P$)
ParallelSort($\{P,W,keys\}, id, s$)Spatial Decomp.; sorting particle-work-key tuples by key
ParallelCostzones($P, W, id, s$) Load balance
$T \gets$ build hashed octree from $P$ Parallel Tree Build
DistributeBranchNodes($P$, $T$, $id$, $s$)
$i$ from 0 to $|P|-1$ Compute Forces
$p \gets P[i]$; $W[i] \gets 0$
Set the acceleration of $p$ to $\vec{0}$.
$Node$ $\gets$ the root node of $T$

$\quad\ $Compute Cell Force:
children of $Node$ are not stored locally
Fulfill any waiting requests for local children from other processes avoid deadlock
Request and receive the children of $Node$ from its owning process
each child cell $C$ of $Node$ Tree traversal unchanged from Algorithm <ref>
not all processes done force computation
Fulfill any waiting requests for local children from other processes
each particle $p$ in $P$ Position Update
Update position and velocity of $p$ using its current acceleration.

§.§ Latency-Hiding Tree Traversal

The hashed octree method
avoids each process from needing
the global octree for force computation.
Instead, branch nodes are used as entry points
to remote trees, where any child node
can be dynamically requested from the source process
during the tree traversal.
In Algorithm <ref>
that dynamic request is made simple through synchronous
and blocking communication.
We now suggest a latency-hiding scheme
which uses asynchronous communication to avoid communication overheads.

The concept is simple:
during the tree traversal for a particle,
maintain a walk list—a queue of nodes which have yet to be traversed—and a
defer list—the queue of nodes whose children have been requested from other processes.
This strategy also transforms the tree traversal
from a recursive method to an iterative one.
Algorithm <ref> shows this iterative tree traversal
to compute the acceleration of one particular body.
The walk list begins containing the root node.
For each iteration of the traversal,
the children of the active node are iterated over,
testing for multipole acceptance and thus a particle-cell
interaction. If the particle-cell interaction is not acceptable,
the child is added to the walk list.
If, at any point, a node is encountered whose
child cells are not stored locally,
then an asynchronous request is sent to the source process
for the children, and the node is enqueued to the defer list.
The traversal then skips to the next iteration.
If, at any point, the walk list becomes empty,
then a node is dequeued from the defer list.
If the child data of that dequeued node has not yet arrived, the traversal waits synchronously for it.
The traversal then continues as normal.

AsyncHashedOctreeTraversal($p$, $T$)$p$ is the body
for which to compute its acceleration, $T$ is the local octree. Returns the work estimate for $p$ for the next time step.
[1]
$w \gets 0$; Set the acceleration of $p$ to $\vec{0}$.
$walk \gets [\text{the root node of}\ T]$; $defer \gets [\ ]$
$walk$ is not empty or $defer$ is not empty
$walk$ is not empty
$Node \gets$ dequeue $walk$
$Node \gets$ dequeue $defer$
Fulfill any waiting requests for local children from other processes
Wait for receipt of children of $Node$ from its owning process
Children of $Node$ are not stored locally
Fulfill any waiting requests for local children from other processes
Send synchronous request for children of $Node$ to its owning process
enqueue $Node$ to $defer$
each child cell $C$ of $Node$
$C$ contains 1 particle $C$ is a leaf
Add to $p$'s acceleration the contribution from $C$'s particle. particle-particle
$w \gets w + 1$
$\ell$ $\gets$ $C$'s side length
$d$ $\gets$ distance between $C$'s c.o.m. and $p$'s position.
${\ell} / {d} < \theta$
Add to $p$'s acceleration the contribution from $C$. particle-cell
$w \gets w + 2$
enqueue $C$ to $walk$

The individual processes operating in this asynchronous
request-and-reply scheme must act cooperatively to achieve best results.
To ensure no request is left waiting for too long, other processes
should frequently check for and fulfill incoming requests for data.
Such a fulfillment occurs every time a process wishes to make a request itself,
thus minimizing access to the communication subprograms and improving locality of operation.

The asynchronous tree traversal easily slots into
the previous Barnes-Hut algorithm based on hashed octrees (Algorithm <ref>).
The modified algorithm is presented in Algorithm <ref>.
Notice that each force calculation step still must end with a synchronization
so that all requests are fulfilled before any process continues
to the integration step (Lines 12–13).
This algorithm is the final version of our many parallel Barnes-Hut methods.

Parallel Barnes-Hut V7($id$, $s$, $\Delta t$, $t_{end}$, $N$, $P_{all}$)$P_{all}$ is the list of size $N$ of bodies to simulate from 0 to $t_{end}$ going by $\Delta t$ time-steps. $id$ is the executing process's id, $s$ is the number of processes. Assume $s \mid N$.
[1]
$N' \gets N / s$
$P \gets$ Scatter $P_{all}$, sending to process $id$ $P_{all}[id\cdot N', \ldots, id\cdot N' + N' - 1]$ Get Local Bodies
$W \gets [N,\dots,N]$ $W$ of size $N'$
$t$ from 0 to $t_{end}$ by $\Delta t$
$keys \gets $ ConstructSpatialKeys($P$)
ParallelSort($\{P,W,keys\}, id, s$)Spatial Decomp.; sorting particle-work-key tuples by key
ParallelCostzones($P, W, id, s$) Load balance
$T \gets$ build hashed octree from $P$ Parallel Tree Build
DistributeBranchNodes($P$, $T$, $id$, $s$)
$i$ from 0 to $|P|-1$ Compute Forces
$W[i] \gets$ AsyncHashedOctreeTraversal($P[i]$, $T$)
not all processes done force computation
Fulfill any waiting requests for local children from other processes
each particle $p$ in $P$ Position Update
Update position and velocity of $p$ using its current acceleration.

§ IMPLEMENTATION

§.§ Data Representation

§.§ Spatial Decomposition and Morton Ordering

§.§ Octrees and Hashed Octrees

§.§ Parallel Sort

The implementation of our treecode is written from scratch in
the C/C++ language, making use of no external libraries.
Nearly everything is written in C, except for the
the use of C++11 threads. Our implementation
closely follows all of the preceding discussion,
including the Barnes-Hut method
described in Section <ref> and Algorithm <ref>,
the leapfrog integrator of Section <ref>,
and the parallel aspects described in Section <ref>.
We take this section to describe implementation-level details
which differ slightly in practice from their conceptual description
or from their original description in <cit.> and <cit.>.
In particular, we detail here our data structures,
the use of Morton ordering to compute costzones,
our particular tree merging technique,
and lastly, our thread implementation using the C++11 built-in thread support library.
Finally, we conclude with an overview of our OpenGL visualization tool.

Our main data structure is a set of parallel arrays (i.e. structure of arrays)
maintaining the masses, positions, velocities, and accelerations of the particles.
For a particle of index $i$, its values are represented
in positions $3i$, $3i + 1$, and $3i + 2$ for the x, y, and z components,
respectively, for each of the latter three arrays,
while mass is a simple array of scalar values.
This data layout allows our code to make use of spatial data locality
throughout the simulation. In particular, it also
allows the use of each type of value independently
for even better locality (cf. an array of structures).
For example, when building the tree only positions are needed,
or when computing forces velocity is not needed.
Using this structure for the integration time-step is simple and follows
exactly from the equations.

The only other data structure for our simulation is the tree itself.
In our case it is an octree in three dimensions.
Our octree is a simple recursive data structure with
each cell storing several values:

* its centre in physical space,

* its side length,

* its centre of mass,

* the total mass of its children,

* the total number of particles contained in its children,

* the quadrupole coefficients (See Equation (<ref>)), and

* a list of its 8 children if it is an internal cell.

To construct the octree in parallel we first determine the costzones.
This begins with an estimate of work for each particle.
We augment our structure of arrays with one more array holding the work estimate of each particle. This array is automatically
updated during the force calculation of the previous step to count
the number of interactions
and serve as the estimate of work for the current step.[In fact, it counts not only the number of interactions but an integer value representing a work estimate. Since particle-cell interactions use quadrupole moments that particle-particle interactions do not, a particle-cell interaction is counted twice to make up for this work difference.]

Next, we must obtain a linear ordering of the particles.
Directly using their indices as an order will not account for physical locality
as the simulation evolves. We thus employ the Morton ordering.
In code this is referred to a spatial key.
The Morton ordering is a bit-wise interleaving of multi-dimensional
integer coordinates to yield a single (one-dimensional) integer which acts as key
to give the Morton ordering for the coordinates
creating the key.
On a 64-bit machine, it is useful for this one-dimensional value
to be a single machine word. Thus each key is 64 bits and, in three dimensions,
each coordinate is a 21-bit number.
We begin by determining the bounding box of all points in the system.
Let $R$ be the smallest value such that
all points in the system are contained in the cube with extents $[-R,R]$.
We then map each particle's floating-point coordinates from the range $[-R,R]$ to $[0, 2R]$
to $[0, 2^{21})$ and finally to an integer by taking its ceiling.
Given these three integers, it is now simple to produce the spatial key
using Morton order bit-interleave (see <cit.> for details)
or any other space-filling curve.
Each array in our structure of arrays is then sorted based on the ordering
induced by the spatial keys.

Finally, the costzones for this time-step are determined.
One simply iterates through the now sorted
array of work estimates, aggregating contiguous partitions
so that each partition holds $\nicefrac{1}{p}$ of the total work, given $p$ processors.

Not only does this sorting yield an ordering of the particles in physical space,
but, since the structure of arrays is sorted by the spatial keys,
spatial data locality is obtained within each costzone.
As we know, improved data locality yields improved performance.
With the costzones computed, we can now build each processor's
local octree in parallel.

Building each local tree simply follows the first step of Algorithm <ref>.
Then, a global tree must be computed by merging each local tree.
In contrast with the original work in <cit.>,
we avoid any use of mutual exclusion to a global data structure through
a restructuring of the merge.
Whereas all local trees were merged into a single global one,
we merge trees pair-wise, in parallel, until only one remains, naming that
the global one.
The local tree construction and merge thus follows the celebrated map-reduce
parallel pattern <cit.>.

Lastly, we have also added an optimization where
trees can be constructed in-place, reusing the existing nodes and structure
of the trees from the previous step. In order to employ this, the bounding box
of all the points must be scaled up slightly so that the size of the root node
does not change at every step. Only when the root node changes size
do the trees need to be truly created from scratch again.
With the octree built and costzones computed, it is now
straightforward to compute the forces on each particle in parallel,
using one thread per costzone.

§.§ C++ Thread Support for Parallelism

For each step in the $N$-body simulation there are several
distinct parallel steps. There is the map-reduce
creation of the global octree, the force computation, and the integration step.
To spawn and join threads for each of these steps would lead
to too much overhead.
Instead, we employ a thread pool—a collection of long-running
threads which wait to be given a task, execute that task,
and then return to the pool. This avoids
the overhead of repeatedly spawning threads.
Often, threads in a pool execute some statically-determined task.
However, our threads must be receptive to each of parallel steps.
We thus implement functor-executing threads.
A functor is simply an object encoding a function call.
Thus, these threads may execute any function that can be
encapsulated as a functor.
To facilitate the object passing of functors to threads,
we have implemented an asynchronous producer-consumer queue.
Nonetheless, these are all standard constructs in modern multithreaded programming,
so we exclude further details.

Our threads, functors, and asynchronous queues are all built
using standard features of C++11. In particular,
(which wraps  on Linux),
for functor objects,
and  for sleeping and
waking threads waiting for a new functor.

§.§ Three Dimensional Visualization

For both debugging and the natural beauty of it,
we have implemented a visualization tool for
our three-dimensional $N$-body simulation
using OpenGL. This visualization tool
has already been shown in Figure <ref>