## Barnes-Hut: An Implementation and Study

### Tree walk and unbounded errors

The Barnes Hut force calculation stage consists of taking each body and walking the tree performing the interaction calculation with the body and other bodies or nodes. The interact subroutine looks like:
```void Interact(body *b1,vector b2pos,double b2mass)
{
double maga; /* magnitude of acceleration */
double r; /* distance from b1 to b2 */
double r2; /* r*r */
vector b1b2; /* vector from b1 to b2 (b2.pos - b1.pos) */
vector da; /* the change that will be applied to the acceleration vector (delta a) */

/* b1 is the body that will be moved, b2 is the body causing the change in acceleration */

/* Given that F = (G m1 m2) / (R^2), and F = m1 a, we can calculate
|a| = G m2 / R^2, and then use that magnitude change to adjust the acceleration
along the direction of the interaction */

SUBV(b1b2,b2pos,b1->pos);
DOTVP(r2,b1b2,b1b2);
r2 += constants.epssq; /* Weaken closeness of the two objects */
r = sqrt(r2);
maga = G * b2mass / (r2*r); /* Account for b1b2 not being a unit vector
in maga */
MULVS(da,b1b2,maga); /* make da the right length */
ADDV(b1->acc,b1->acc,da); /* Accumulate the acceleration */
}
```
This function will correctly calculate the interaction of any two bodies. The tree-walk optimization pre-computes the center of mass of a group of bodies and calculates the interaction of a single body with that entire group all in a single step. The criteria for whether or not to recurse down the oct-tree or perform the interaction with the entire group looks like:
```void TW_DownWalk(node *global root,body *particle,double dsq)
{
SUBV(dr,root->pos,particle->pos);
DOTVP(drsq,dr,dr);
if (constants.tolsq * drsq < dsq) {
/* need to open node, i.e. recurse down the tree */
} else {
/* perform the interaction */
}
}
```
The important piece to note is that the Interaction function uses the mass of the second body in the calculation to determine the amount of acceleration, but the opening criteria (whether or not to recurse) does not take into account the accumulated mass at the node, but only the tolerance. If the simulation is sufficiently well behaved to keep large clumps from occurring, then having a single tolerance is sufficient. However, since most of these simulations move from randomly distributed particles to clustered particles, it is important to take the mass into account. Our implementation does not do so, and as a result has the same problem with unbounded error as many other implementations. An area for further work would be to either use the opening criteria in Warren and Salmon's paper, or work out another criteria for opening the node. Further reduction in the amount of work for the Barnes-Hut algorithm could be found by performing node-node interactions when the two nodes are sufficiently far apart as to minimize the error.

-----

maintained by <hodes@cs.berkeley.edu>