Advertisement

The Barnes-Hut Galaxy Simulator

An Introduction to solving large scale N-Body Problems

The N-Body Problem

A region of space contains the number of N bodies. Each body has its own potential and a resulting force field. This could be charges causing an electrical field (Coulomb's law) or planets in space (Law of Gravity). The potential of the bodies can be summed up yielding a combined potential which depends on the location of each of the bodies and its physical properties (i.e. mass or charge). According the Newtons first law the bodies themself will experience an acceleration caused by the field. As a result there is a force field that depends on the locations of the n-nodies whilst the location of the body depends on the force field. Mathematically this problem can be described with a differential equation of second order.

When solving such problems one usually starts with a known initial distribution of the bodies and known velocities and accelerations. From the initial conditions the system can be advanced in time by using numerical computations and a suitable integration scheme. Typical problems of this kind are for instance the computation of galaxy evolution in astrophysics or the behaviour of molecules in biology. Usually the equations of motion are relatively simple. For instance for simulating a galaxy they can be formulated using the Law of gravity:

Die Bewegungsgleichungen sind in aller Regel recht einfach aufgestellt. Wenn man beispielsweise eine Galaxie simulieren möchte, dann lautet die Gleichung um die Kraft auf ein Teilchen der Gesamtmenge zu berechnen:

\begin{equation} m_i \ddot{\vec{r_i}} = G \sum_{j=1;j \ne i}^N \Big( m_i m_j \frac{ \vec{r_i} - \vec{r_j}} {|\vec{r_i} - \vec{r_j}|^3} \Big) \end{equation}

with:

Lets assume Fij is the force acting between particles i and j. The total number of force calculations needed to compute the state of the system is N*(N-1) and according to Newtons third law every force has an opposite forces equal to itself: Fij = -Fji. The total number of force calculations can then be reduced to:

\begin{equation} {1 \over 2 } N (N-1) \approx O(N^2) \end{equation}

The problem is of order O(N2). If the number of particles double the number of calculations quadrupels. If the number of particles is increased by factor ten the number of calculations increases by a factor of 100 and so on... From this simple relation it is clear that computing the N-Body problem for large numbers of particles will quickly become very costly in terms of computational power. A more effective algorithm that scales better with increasing number of particles is needed.

The Barnes-Hut Algorithm

The Barnes-Hut Algorithm describes an effective methood for solving n-body problems. It was originally published in 1986 by Josh Barnes and Piet Hut [1]. Instead of directly summing up all forces it is using a tree based approximation scheme which reduces the computational complexity of the problem from O(N2) to O(N log N).

The Barnes-Hut algorithm reduces the number of force calculations by grouping particles. The basic idea behind the algorithm is that the force which a particle group excerts on a single particle can be approximated by the force of a pseudo particle located at the groups center of mass. For instance the force which the andromeda galaxy excerts on the milky way can be approximated by a point mass located at the centre of the andromeda galaxy. There is no need to integrate over all stars in the andromeda galaxy provided the distance between the two galaxies is large enough. This approximation is valid as long as the distance from the point group to the particle is large and the radius of the group is small in relation to the distance between the group and the particle. The ratio of group distance (d) and group radius (r) is called the Multipole-Acceptance-Criterion (MAC):

$$\theta = {d \over r }$$

If theta drops below a certain threshold the quality of the approximation starts to deteriorate resulting in larger errors. Barnes and hut use a value of θ=1. The general rule is: The smaller θ, the better the simulation results. The Barnes-Hut algorithm applies this basic principle by subdividing the space of the simulation domain recursively (see Image 2). This is described more indepth in the next chapter.

Image 1: The force a particle group excerts on a single particle at distance r can be approximated by the force of a point mass provided that the distance between the group and the particle is large in comparison with the group size d.
Image 2: The gravitational effect excerted by the star cluster and the single star B on star A can be approximated by a point mass. On a smaller scale the gravitational effect of the cluster on star B can also be approximated by a point mass since the distance of B from the cluster is much large compared to the size of the cluster (r2 >> d2).

Implementation

The following section will explain the algorithm. All explanations and examples are simplified and deal with a two dimensional domain. Extending the concept to a third dimension is pretty straightforward but not the scope of this article.

The tree structure

Lets start with a random distribution of points in a two dimensional domain. The first step of the Barnes-Hut algorithm is sorting the particles into a hierarchical tree structure. This structure is the so called quadtree (or octree in 3D). The algorithm starts by subdividing the domain into quadrants. Based on their location I will refer to these quadrants by the following names: NW, NO, SW, SO. In order to add a particle to the tree its quadrant is determined. If there is another particle in the same quadrant the quadrant is subdivided again. This is repeated until each particle is located in its own quadrant (however small this quadrant may be). Image 3 shows the result of this process for a random distribution of particles.

Image 3: right: Spacial distributiong of particles in their quadrants; left: The tree structure of the distribution.

The creation of the quadtree is a recursive process. It can be described by the following pseudocode:

  
Function MainApp::BuildTree
{
  ResetTree

  for all particles
    rootNode->InsertToNode(particle)
  end for
}

Function TreeNode::InsertToNode(newParticle)
{
  if number of particles in this node > 1
  {
    quad = GetQuadrant(newParticle);
    
    if subnode(quad) does not exist
      create subnode(quad)

    subnode(quad)->InsertToNode(newParticle)
  }
  else if number of particles in this node == 1
  {
    quad = GetQuadrant(existingParticle);
    if subnode(quad) does not exist
      create subnode(quad)
    subnode(quad)->InsertToNode(existingParticle)

    quad = GetQuadrant(newParticle);
    if subnode(quad) does not exist
      create subnode(quad)
    subnode(quad)->InsertToNode(newParticle);
  }
  else if node is empty
  {
    store newParticle in node as existingParticle
  }

  Increase number of particles
}

The tree has to be recreated for each step of the simulation. The computational cost of recreating the quadtree depends on the distribution of the particles (see also [2]). The cost of adding a particle to the tree is proportional to its distance from the root node. Distributions with many densly packed particles require more operations because the tree must be subdivided often to place all particles in their own subdivided quadrant.

Computing the mass distribution

Up until now the tree only contains the spacial distribution of the particles. The next step deals with calculating the mass distribution of the quadtree (or. octtree). It sonsists of computing the total mass of all particles contained in the child nodes of each tree node as well as their center of mass. This data must be computed for each tree node. The calculation can be implemented recursively by scanning all child nodes of each treenode. The following pseudo code should illustrate the algorithm:

 
Function TreeNode::ComputeMassDistribution()
{
  if number of particles in this node equals 1
  {
    CenterOfMass = Particle.Position
    Mass = Particle.Mass  
  }
  else
  {
    // Compute the center of mass based on the masses of all child quadrants
    // and the center of mass as the center of mass of the child quadrants 
    // weightes with their mass
    for all child quadrants that have particles in them
    {
      Quadrant.ComputeMassDistribution
      Mass += Quadrant.Mass
      CenterOfMass +=  Quadrant.Mass * Quadrant.CenterOfMass
    }
    CenterOfMass /= Mass
  }
}

This procedure is of order O(N log N) [2].

Force calculation

You probably guessed it by now the last step of the algorithm is computing the forces excerted on the particles. This is what makes the Barnes-Hut algorithm faster than the direct summation of all mutual forces. Do you remember the images 1 and 2 ? This is the basic principle we are going to implement here. In order to computer the force excerted on a single particle the quadtree must again be iterated recursively. For each node θ is computed:

$$\theta = {d \over r }$$ where:

If θ lies below a certain threshold the force can be approximated by inserting the quadrant nodes mass and its center of mass in the law of gravity. The child nodes don't have to be summed up separately. A reasonable relation is θ < 1. If θ is larger than the threshold the quadrants effect can't be approximated and all of its child nodes have to be tested again. The iteration stops only after all nodes have been tested. The worst case scenario is having to test all nodes without finding any node that meets the MAC. In such a case the result is similar to summing up all mutual particle forces (θ=0). The iteration depth can be finetuned by adjusting θ. Animation 1 illustrates the influence of θ on the number of force computations.

Abhängigkeit der Anzahl der Kraftberechnungen von der Wahl von theta
Animation 1: The animation shows a distribution of 5000 particles. The quadrants that are shown are the ones that are used for calculating the force excerted on the particle at the origin of the coordinate system. The higher θ is, the fewer the number of nodes that are necessary for the force calculation (and the larger the error).

Finally the pseude code for the force calculation:

 
Function MainApp::CalcForce
  for all particles
    force = RootNode.CalculateForceFromTree(particle)
  end for
end 

Function force = TreeNode::CalculateForce(targetParticle)
  force = 0
  if number of particle equals 1
    force = Gravitational force between targetParticle and particle
  else
    r = distance from nodes center of mass to targetParticle
    d = height of the node
    if (d/r < θ)
      force = Gravitational force between targetParticle and node using the nodes center of mass and the total mass of the node
    else
      for all child nodes n
        force += n.CalculateForce(particle)
      end for
    end if
  end
end

Colliding galaxies

By combining the Barnes-Hut-Algorithm with a suitable integration scheme (i.e. Adams-Bashforth) one can relatively easy compute a two dimensional galaxy collision:


You might also like: