THE FAST MULTIPOLE METHOD (C++)
UNIVERSITY PROJECT
The fast multipole method is an "approximation algorithm" to solve n-body problems. N-body problems are all about the interactions between a large set of bodies. If each body interacts with every other body, calculating these interactions can be very computationally expensive.
The FMM is a method that is complex, but hopefully I can explain my implementation in a simple way. I implemented this algorithm as part of my n-body solution application that I made for my dissertation.
THE FORCE OF GRAVITY
In our large set of particles, the force acting on any given particle i, due to any other particle j, is given in vector form as:
​
F = G*mi*mj / r^3
​
where G is the gravitational constant,
mi is the mass of body i,
mj is the mass of body j,
and r is the vector from i to j.
​
This gives us the value F, which is the force vector acting on the body i. So calculating the force vector for every particle is going to really expensive, we'll have to run this function for each pair, giving a time complexity of O(N^2). This is where the approximation comes in.
OCTREE
Similar to the Barnes-Hut algorithm we are going to create an octree: The octree contains a single root node, and is broken down into 8 children nodes, these nodes are again broken down further. In a Barnes-Hut simulation we will break each node down until there is 1 or 0 bodies in a node. However in a FMM we don't break down based on particle count, we'll just break the node down to a certain layer.
​
If we break the simulation down 3 times, we will have 3 different layers: Our first layer has the single root node, our second layer has the 8 children, and the third layer contains 64 nodes - the 8 children of our previous 8 children.
​
As you can see there is a basic node count relation here. The number of nodes on any given layer is 8 to the power of the layer depth.
​
Cell count = 8^L
​
where L is the depth.
CELL INTERACTIONS
The cells on each layer are going to interact with each other. For every layer a cell will approximate the forces produced by far away cells. For each cell we will create an interaction list, this list will contain all the other cells that a given cell should approximate its forces with. This means when it comes to actually running the algorithm, the algorithm will loop over every single cell and approximate the forces due to the cells in its interaction list. For any given cell i, another cell j is considered "far away" if it follows 3 rules:
-
Cell j is on the same layer as cell i.
-
The parent of cell j is a neighbour of the parent of cell i.
-
Cell i and cell j are not neighbours.
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
​
[1]
In this example, light grey cells are cells in B's interaction list,
dark grey cells are cells that have been approximated by B's parent,
white cells are neighbours and are not approximated.
​
​
So how do we build the interaction list for a cell? Since we are not using a dynamic octree we can build this list for each cell before the algorithm actually starts running, and this also means it will never change so we can just calculate it once.
​
To build this list you can do a nested for loop for i and j and check if all 3 rules are being followed. I did it in a slightly more efficient way. I perform a breadth-first search of the octree. This search algorithm looks through a tree structure layer by layer; This means that if I wanted to extract a list of cells on any given layer, I can stop the search when the queue has a length of 8^L.
​
Stopping this search at each layer, and doing a nested "for loop" for i and j to check for rules 1 and 2, builds the interaction list for each cell slightly faster.
THE MULTIPOLES
In this context: A multipole is a way of accurately expressing the centre position, and the centre of mass of a cell. The multipole is expressed by two components: the inner expansion, and the outer expansion. The multipole is expressed as:
[2]
​
The components containing the values of j are outer expansion expressions, and the inner expansion are the expressions containing i. Therefore the full multipole would be the dot product of these expansions.
​
Each frame, we perform the upward sweep, this is where we start at each cell on the bottom level and calculate the inner expansion of a cell. For the cells on other layers: their multipoles can be calculated using their children. This is why it is called the upward sweep, as we calculate the ones at the bottom and then work our way up.
​
After this upward sweep, we do our downward sweep, where we start at the root node and work our way down, calculating the approximate force acting on a cell due to its interaction list. At the very bottom level we must also do a brute force calculation on the cells that are neighbours!
RESULTS
I calculated the time for the algorithm to complete at 100 bodies, incrementing by 100 until 1,000; Then incrementing by 1,000 until at 10,000 bodies; Then incrementing by 10,000 until at 100,000 bodies.
The accuracy of the algorithm was determined by also calculating force using the regular force calculation and comparing the results. These results created an accuracy error value using this equation:
In this context, "potential" just means force.
You may notice that the error seems to decrease over time. This is not a feature of the algorithm but just a side effect of the simulation space. For fast multipole methods, the algorithm has optimal performance when the bodies are evenly distributed [1]. Since we are adding more and more bodies over time, the simulation space is becoming more dense, and thus more direct calculations are being done with neighbours at the lowest level; Since more direct calculations are being done, the accuracy increases, and the error decreases.
REFERENCES
[1] Ying, L. (2012) “A pedestrian introduction to fast multipole methods.”, Sci China Math, 55(5), pp. 1043-1051. doi: 10.1007/s11425-012-4392-0
​
[2] Clementi, N. (2014) “FMM Tutorial” Available at: https://github.com/barbagroup/FMM_tutorial (Accessed 27 March 2021)