# 并行N-body模拟

## 算法

### 加速度计算

According to Newton's law of universal gravitation, every body attracts every other body by a force pointing along the line intersecting both centers. The force is proportional to the product of the two masses and inversely proportional to the square of the distance between them:

$\begin{eqnarray} F = G \frac{M m}{r^2} \end{eqnarray}$

By Newton's Second Law, the acceleration is expressed by $$F = m a$$.

Combing the two laws, we can easily deduce an direct-sum algorithm to calculate the accelerations of these bodies, which would be $$O(n^2)$$. However, the brute-force simulation has low data dependency and is easy to be reformed to a parallelized algorithm.

Barnes-Hut simulation decreases the time complexity to $$O(n \log n)$$ which exploting a quad-tree (in the 2-dimensional case) and approximate a distant cluster of bodies to a point. Each node in the quad-tree represents a subspace. The root node represents the whole space and its four children represent the four quadrants of the space. The resultant force attracted by a distant cluster of bodies is treated as being attracted by the center of mass.

To determine if the cluster is sufficiently far away, I compute the quotient $$s/d$$, where $$s$$ is the size of the subspace and $$d$$ the distance between the body and the center of mass. If $$s/d < \theta$$, the cluster is considered distant.

### 碰撞检测

Each body along with the imaginary position when it reaches the destination after the time slice, is enclosed in a axis-aligned minimum bounding rectangle. It is easy to decuce that if two body collides in the time slice, then their bounding rectangles intersect. The property can be used to accelerate the performance of collision detecting.

Several existing approaches deal with the intersection detection problem straightforwardly.

• R-tree and its variants. R-trees are tree data structures used for spatial access methods and the most immediate solution to the problem.
• Binary space partitioning tree. They are developed in the context of computer graphics are are also suitable for the problem.
• 2-dimensional segment tree. First notice that two rectangles overlaps iff their projections on coordinate axes overlap. We need to construct a 2-dimensional segment tree with $$O(n \log^2 n)$$ space which supports $$O(n \log^2 n)$$ query.
• Spatial hashing. Subdivide the plane into a set of voxels (a term used in image processing, the smallest distinguishable box-shaped part) and thus forms a mesh. Each voxel contains a bucket of bodies. Put each body into the buckets corresponsing its occupied voxels. Before insertng into a voxel, each existing body in the cell is a candidate for potential collision.

Other approaches including using a sweep line algorithm to reduce the problem to the 1-dimensional case (overlapping interval search problem).

• Interval tree. It is a tree data structure holding intervals. Space complexity $$\Theta(n)$$, construction time complexity $$O(n \log n)$$, query time complexity $$O(\log n + k)$$.
• Segment tree. Space/construction time complexity $$O(n \log n)$$, query time complexity $$O(\log n + k)$$.

Approaches exploiting the property of temporal coherence:

• TPR-tree and a few ramification. A variant of R-tree taking advantages of spatio-temporal updates.

I'll expatiate three data structures implemented in my project.

#### Binary space partitioning tree

The implementation is in the file =src/BSPTree.cc=.

The construction process is rather simple and straightforward. I start with the whole scene and a list of bounding boxes, recursively divide the scene into two until the number of bodies in the subspace is less than a threshold (=NODE_SIZE_THRESHOLD=) or there is no balanced partition scheme. Each separating line is parallel to X, Y or Z axis.

For a range query, we examine whether the cube in question intersects with the left-subspace. If so, recursively do a query in the left-subspace. The same process also applies to the right-subspace.

#### Sweep line algorithm + Interval tree

The implementation is in the file =src/IntervalTree.cc=.

Well, in this project I also use a much easier but also efficient approach combining a sweep line algorithm and an interval tree besides a BSP tree implementation.

• Each bounding rectangle has an upper and a lower edges. They mark the timestamp between which the rectangle is valid.
• Sort these edges by their y coordinates and process these edges in ascending order. For an upper edge, search all the valid rectangles and find those x-axis spans (intervals) which overlaps the span of the current rectangle, and then insert the x-axis span into a data structure. For these overlapping intervals, check whether the corresponding rectangle is marked valid and if so, do a precise detection. For a lower edge, delete the x-axis span from the data structure.

The unamed data structure should support following operations efficiently: insert an interval, query all intervals overlapping given interval (overlapping interval search problem).

The query can be fulfilled by an interval tree. Construct an interval tree with all intervals but do not really store these intervals. This makes the constructed tree balanced.

Suppose $$n$$ is the number of bodies. The interval tree can be constructed in $$\Theta(n \log n)$$ time while each insertion runs in $$O(\log n)$$ time and each query runs in $$O(\log n + k)$$ time where $$k$$ is the number of overlapping intervals.

#### Sweep line algorithm + Segment tree

Ditto on the reduction part. Create a sorted list of unique X coordinates extracted from those bounding boxes: $$p_0, p_1, p_2, \ldots, p_{m}$$. These coordinates form $$m$$ segments: $$[p_0,p_1), [p_1,p_2), \ldots, [p_{m-1},p_m)$$. Each leaf node represents a segment while each internal node represents several segments.

### 碰撞模拟

The most accurate approach is to deduce the nearest collision event, simulate the event, and then find another. However, this is incapable of handling hundreds of bodies, especially when there are some high-speed and light bodies. They will frequently change directions and affect other bodies.

A comparatively coarse one is to calculate $$n(n-1)/2$$ pair-wise collision events and simulate them one after another. It is coarse in the sense that it may not handle multiple collisions taking place on one single body correctly.

The two mentioned ways are neither efficient nor parallelizable. One simple mend making the latter parallelizable is to update velocities after all collision events are processed, but incurs the inaccuracy that the updated velocity of the other body in a collision won't be reflected immediately.

I used a non-parallelizable, asymptotically optimal sequential approach that is a bit accurate then the above one.

By the way, to make the calculation of the collision point more accurate, I'll use a techinique stated below:

Suppose $$\vec{p_0}, \vec{p_0}$$ are positions of the two bodies and $$\vec{\Delta p}, \vec{\Delta q}$$ their translational vectors. Then we can parametrize their positions:

$\begin{eqnarray} \vec{p}(t) = \vec{p_0} + t \vec{\Delta p} \\ \vec{q}(t) = \vec{q_0} + t \vec{\Delta q} \end{eqnarray}$

We will try to solve the equation $$|\vec{p}(t)-\vec{q}(t)|=R+r$$ where $$R,r$$ are radii of the two bodies. This is a quadratic equation and we should care whether the smaller root $$t_0$$ satisfies $$0 \le t_0 \le 1$$ or the two bodies overlap. If so, teleports the two bodies to their collision moment which is $$now + t_0 \Delta t$$ and carry out an elastic collision.

Since collision only imparts force along the line of collision, only the velocities that are along the line of collision should be changed, thus degrades to 1-dimensional case.

## 截图

src/n-body -f text -d 3 -n 20 --fps 20 -t 2 | python scripts/a.py