CSC/ECE 506 Spring 2011/ch4a zz
Introduction
The N-body problem stated as follows: Select the position and velocity of n celestial bodies as states. Given the initial condition of N bodies, compute their states at arbitrary time T. Normally a three-dimensional space is considered for N-body problem. There is a simplified N-body problem called restricted N-body problem where the mass of some of the bodies is negligible. Several remarkable three-body simulation can be found in [1].
Many mathematicians have proofed that it is impossible to find a general solution for n-body problem analytically[2][3]. The system could become unstable very easily. However, the problem can be solved numerically. The most common approach is to iterate over a sequence of small time steps. Within each time step, the acceleration on a body is approximated by the transient acceleration in the previous time step. The transient acceleration on a single body can be directly computed by summing the gravity from each of the other N-1 bodies. While this method is conceptually simple and is the algorithm of choice for many applications, its O(N2)
The simulation of N-body system can be used from simulation of celestial bodies (gravitational interaction)to interactions of a set of particles (electromagnetic interaction).
Parallel N-body problem
Since each particle will interact with others by the force of gravity, the simulation of N-body system is computationally expensive for large numbers of N. There a O(N2) interactions to compute for every iteration. Furthermore, in order to have a accurate result, the discrete time step must relatively small. Thus, there has been a huge interest in faster parallel algorithm for N-body problem.
Early parallel algorithms
In 1985, Appel took the first step of decomposes the problem by introducing a tree structure[4]. In the next year, Barnes and Hut extend the tree-based force calculation with logarithmic growth of force terms per particle[5]. The partition is shown on the left. The empty block have been pruned, thus the traversal time have been reduced fromO(N2) to O(N logN).
For each time step: 1. Build the BH tree. 2. Compute centers-of-mass bottom-up 3. For each body Star a depth-first traversal of the tree; Truncating the search at internal nodes where the approximation is applicable; Update the contribution of the node to the acceleration of the body. 4. Update the velocity and position of each body.
Data-parallel Implementation
One obvious approach for the data-parallel implementation is to divide the interactions into different sets based on distance. Thus the force from 'far away' particles can be updated less frequently or even can be ignored.
By separating the blocks by the Hamming distance to the R block. Three interaction list could be created by the distance: The direct interaction list Dr; the far interaction list Fr and the intermediate interaction list Hr
The steps for each iteration using the given decomposition of the simulation are [6]:
- Compute the multipole expansion coefficients for all leaves in the tree decomposition.
- Compute the multipole expansion coefficients for all internal nodes in the tree with depth ≥ 2.
- Compute the local expansion coefficients for a region R by summing R’s parent’s local expansion (shifted from the parent’s center to the center of R) with the sum of all of the multipole expansions in R’s far list, Fr(converting the multipole expansions to local expansions and shifting to the center of R) for all regions with depth ≥ 2.
- For each body, b, in each leaf region, R, compute all the direct forces on b from all the bodies in the regions in R’s direct interaction list (Dr).
- For each body, b, in each leaf region, R, compute the far force on b by evaluating the local expansion for region R at b’s position.
- For each body, b, in each leaf region, R, compute the intermediate force by evaluating the multipole expansion at b’s position for each region in R’s intermediate interaction list (Hr).
- Sum the 3 components of the force and potential for each body.
- Apply the forces, updating the positions and velocities, and move the bodies to their proper regions as indicated by boundary crossing.
Another data-parallel implementation can be found in [8]. The drawback of this type of approach is also obvious: the error may be increased significantly due to the clustering process.
Message-passing Implementation
In message-passing architectures, a processor’s local memory is often used as a software-controlled cache for communicated data, in addition to holding the processor’s own assigned data partition; that is, data are replicated in the memories of those processors that need to use them.
Warren and Salmon introduced a parallel hashed Oct-tree N-body algorithm. One of the advantages of using hash table is that the non-local data may be accessed by requesting a key, which is a uniform addressing scheme and easy to be implemented on a message passing architecture. The pseudocode of tree traversal are shown below:
ListTraverse((*MAC)(hcell *))
{
copy root to walk_list;
while (!Empty(walk_list)) {
for (each item on walk_list) {
for (each daughter of item) {
if (MAC(daughter))
copy daughter to interact_list;
else
copy daughter to output_walk_list;
}
}
walk_list = output_walk_list;
}
}
There are two major challenges of using message-passing parallel to solve the N-body problem:
- the nonuniform and dynamic nature of the particle distribution, which implies that the partitions assigned to processors change with time to maintain load balancing and data locality;
- the need for long-range communication which is irregular as a result of the nonuniformity.
The typical message-passing strategy of allowing communicated data to accumulate in local memory and flushing them out at certain points in the program does not work well here like it does in regular predictable programs. Thus, the only convenient points at which to flush communicated data in the application program are the boundaries between computational phases. Since the amount of nonlocal data read from the tree during the force computation phase can be quite large, the memory overhead due to data replication can be large when replacement is managed in this typical message-passing style.
Besides the early achievement of implementing the N-body simulation on shard-memory machines, there is a tremendous interests of using CUDA to solve the N-body problem by modern GPUs [10][11].
__global__ void
calculate_forces(void *devX, void *devA)
{
extern __shared__ float4[] shPosition;
float4 *globalX = (float4 *)devX;
float4 *globalA = (float4 *)devA;
float4 myPosition;
int i, tile;
float3 acc = {0.0f, 0.0f, 0.0f};
int gtid = blockIdx.x * blockDim.x + threadIdx.x;
myPosition = globalX[gtid];
for (i = 0, tile = 0; i < N; i += p, tile++) {
int idx = tile * blockDim.x + threadIdx.x;
shPosition[threadIdx.x] = globalX[idx];
__syncthreads();
acc = tile_calculation(myPosition, acc);
__syncthreads();
}
// Save the result in global memory for the integration step.
float4 acc4 = {acc.x, acc.y, acc.z, 0.0f};
globalA[gtid] = acc4;
}
CUDA shows an impressive potential to execute extremely parallel applications such as N-body problem. Twenty years ago, the N-body simulation can only runs on super computers. Because of the development of CUDA makes the classical N-body problem become a term project for some graduate students (see 2008 Fall CS 525 Project2 n-body simulation in CUDA).
Comparison of Performance
The following table is a summary of sequential and parallel implementations of N-body simulation from [7].
| Author | Method and Error | Program Model | N | Times (s) | P | Efficiency | Machine |
| Salmon (1990) | BH, quadrupole | Message passing | ? | ? | ? | ? | Ncude |
| Warren_Salmon (1993) | BH,ϵ = 10-3 | Message passing | 8.78 M | 114 | 512 | 28% | Intel Delta |
| Warren_Salmon (1995) | BH,ϵ = 10-2 | Message passing | 2 M | 10.8 | 256 | ? | CM-5E |
| Liu-Bhatt (1994) | BH, quadrupole | Message passing | 10 M | 59 | 256 | 30% | CM-5 |
| Leathrum-Board(1992) | GR,p=8 | Shared memory | 1 M | 1520 | 32 | 20% | KSR-1 |
| Elliott-Board(1994) | GR,FFT,p=8 | Shared memory | 1 M | 1420 | 32 | 14% | KSR-1 |
| Schmidt-Lee(1991) | GR,p=8 | ? | 40,000 | 94 | 1 | 39% | CRAY Y-MP 8/864 |
| Zhao-Johnsson(1991) | Zhao, p=3 | Data parallel | 16,000 | 5 | 8 K | 12% | CM-2 |
| Hu-Johnsson | Anderson | Data parallel | 100 M | 180 | 256 | 27%-35% | CM-5E |
| Singh et al.(1993) | GR, 2-D, adaptive | Shared memory | ? | ? | ? | ? | DASH, KSR-1 |
All performance numbers are for uniform particle distributions. Methods used are for 3-D, unless other wise stated. "?" imply unavailable data.
References
[1] Collection of remarkable three-body motions
[4] A. Appel, "An Efficient Program for Many-Body Simulation", SIAM J. Scientific and Statistical Computing, vol. 6, 1985
[5] Barnes. Josh, Hut. Piet, (12/04/1986). "A hierarchical O(N log N) force-calculation algorithm". Nature (London) (0028-0836), 324 (6096), p. 446.
[6] A Data-Parallel Implementation of the Adaptive Fast Multipole Algorithm 1993
[7] Johnsson, SL (01/01/1996). "A data-parallel implementation of hierarchical N-body methods". The international journal of supercomputer applications and high performance computing (1078-3482), 10 (1), p. 3.
[8] Implementing N-body Algorithms Efficiently in Data-Parallel Languages (1996)
[9] A Parallel Hashed Oct-Tree N-Body Algorithm (1993)