README.md

Barnes-Hut Simulation

Use the following commands to make a fresh clone of your repository:

git clone -b barneshut git@gitlab.epfl.ch:lamp/student-repositories-s21/cs206-GASPAR.git cs206-barneshut

Useful links

• The API documentation of the Scala standard library
• The API documentation of the Java standard library

If you have issues with the IDE, try reimporting the build, if you still have problems, use compile in sbt instead.

Introduction

In this assignment, you will implement the parallel Barnes-Hut algorithm for N-body simulation. N-body simulation is a simulation of a system of N particles that interact with physical forces, such as gravity or electrostatic force. Given initial positions and velocities of all the particles (or bodies), the N-body simulation computes the new positions and velocities of the particles as the time progresses. It does so by dividing time into discrete short intervals, and computing the positions of the particles after each interval.

Before we study the Barnes-Hut algorithm for the N-body simulation problem, we will focus on a simpler algorithm -- the direct sum N-body algorithm. The direct sum algorithm consists of multiple iterations, each of which performs the following steps for each particle:

1. The particle position is updated according to its current velocity (delta is a short time period).

    x' = x + v_x * delta
    y' = y + v_y * delta

2. The net force on the particle is computed by adding the individual forces from all the other particles.

    F_x = F_1x + F_2x + F_3x + ... + F_Nx
    F_y = F_1y + F_2y + F_3y + ... + F_Ny

3. The particle velocity is updated according to the net force on that particle.

    v_x' = v_x + F_x / mass * delta
    v_y' = v_y + F_y / mass * delta

In this exercise, we will assume that the force between particles is the gravitational force from classical mechanics. Let's recall the formula for the gravitational force between two stellar bodies:

F = G * (m1 * m2) / distance^2

Above, F is the absolute value of the gravitational force, m1 and m2 are the masses of the two bodies, and r is the distance between them. G is the gravitational constant.

For each particle, the net force is computed by summing the components of individual forces from all other particles, as shown in the following figure:

The direct sum N-body algorithm is very simple, but also inefficient. Since we need to update N particles, and compute N - 1 force contributions for each of those particles, the overall complexity of an iteration step of this algorithm is O(N^2). As the number of particles grows larger, the direct sum N-body algorithm becomes prohibitively expensive.

The Barnes-Hut algorithm is an optimization of the direct sum N-body algorithm, and is based on the following observation:

If a cluster of bodies is sufficiently distant from a body A, the net force on A from that cluster can be approximated with one big body with the mass of all the bodies in the cluster, positioned at the center of mass of the cluster.

This is illustrated in the following figure:

To take advantage of this observation, the Barnes-Hut algorithm relies on a quadtree -- a data structure that divides the space into cells, and answers queries such as 'What is the total mass and the center of mass of all the particles in this cell?'. The following figure shows an example of a quadtree for 6 bodies:

Above, the total force from the bodies B, C, D and E on the body A can be approximated by a single body with mass equal to the sum of masses B, C, D and E, positioned at the center of mass of bodies B, C, D and E. The center of mass (massX, massY) is computed as follows:

mass = m_B + m_C + m_D + m_E
massX = (m_B * x_B + m_C * x_C + m_D * x_D + m_E * x_E) / mass
massY = (m_B * y_B + m_C * y_C + m_D * y_D + m_E * y_E) / mass

An iteration of the Barnes-Hut algorithm is composed of the following steps:

1. Construct the quadtree for the current arrangement of the bodies.
   1. Determine the boundaries, i.e. the square into which all bodies fit.
   2. Construct a quadtree that covers the boundaries and contains all the bodies.
2. Update the bodies -- for each body:
   1. Update the body position according to its current velocity.
   2. Using the quadtree, compute the net force on the body by adding the individual forces from all the other bodies.
   3. Update the velocity according to the net force on that body.

It turns out that, for most spatial distribution of bodies, the expected number of cells that contribute to the net force on a body is log n, so the overall complexity of the Barnes-Hut algorithm is O(n log n).

Now that we covered all the necessary theory, let's finally dig into the implementation! You will implement:

• a quadtree and its combiner data structure
• an operation that computes the total force on a body using the quadtree
• a simulation step of the Barnes-Hut algorithm

Since this assignment consists of multiple components, we will follow the principles of test-driven development and test each component separately, before moving on to the next component. That way, if anything goes wrong, we will more precisely know where the error is. It is always better to detect errors sooner, rather than later.

Data Structures

We will start by implementing the necessary data structures: the quadtree, the body data-type and the sector matrix. You will find the stubs in the package.scala file of the barneshut package.

Quadtree Data Structure

In this part of the assignment, we implement the quadtree data structure, denoted with the abstract data-type Quad. Every Quad represents a square cell of space, and can be one of the following node types:

• an Empty node, which represents an empty quadtree
• a Leaf node, which represents one or more bodies
• a Fork node, which divides a spatial cell into four quadrants

The definition of Quad is as follows:

sealed abstract class Quad {
  def massX: Float
  def massY: Float
  def mass: Float
  def centerX: Float
  def centerY: Float
  def size: Float
  def total: Int
  def insert(b: Body): Quad
}

Here, massX and massY represent the center of mass of the bodies in the respective cell, mass is the total mass of bodies in that cell, centerX and centerY are the coordinates of the center of the cell, size is the length of the side of the cell, and total is the total number of bodies in the cell.

Note that we consider the top left corner to be at coordinate (0, 0). We also consider the x axis to grow to the right and the y axis to the bottom.

The method insert creates a new quadtree which additionally contains the body b, and covers the same area in space as the original quadtree. Quadtree is an immutable data structure -- insert does not modify the existing Quad object. Note that Body has the following signature:

class Body(val mass: Float, val x: Float, val y: Float, val xspeed: Float, val yspeed: Float)

In this part of the exercise, you only need to know about body's position x and y.

Let's start by implementing the simplest Quad type -- the empty quadtree:

case class Empty(centerX: Float, centerY: Float, size: Float) extends Quad

The center and the size of the Empty quadtree are specified in its constructor. The Empty tree does not contain any bodies, so we specify that its center of mass is equal to its center.

Next, let's implement the Fork quadtree:

case class Fork(nw: Quad, ne: Quad, sw: Quad, se: Quad) extends Quad

This node is specified by four child quadtrees nw, ne, sw and se, in the northwest, northeast, southwest and southeast quadrant, respectively.

The northwest is located on the top left, northeast on the top right, southwest on the bottom left and southeast on the bottom right.