Programming Lab #3 (CSE 260, Winter 2014)
Many body simulation with MPI
Due Friday 3/14/14 at 10PM

Changelog

The Assignment

Your assignment is to parallelize, with MPI, the many body particle simulation discussed in class. Your starting point is serial code. You are also provided with a multithreaded executable that will enable you to assess performance on up to 8 cores (1 node of bang).

Conduct scaling studies on up to 64 cores of the Bang cluster, if in groups of 3, on up to 256 cores of the Carver sytem located at NERSC. (Due to queuing delays, 256 cores may be the practical limit). The provided code can also benefit from single processor optimizations. Try to make the code go as fast as you can. In the spirit of friendly competition, post your performance results to the A3 Moodle forum.

For instructions about how to compile and run MPI jobs on Bang, refer to the document Using MPI on Bang.

For instructions about how to compile and run MPI jobs on Carver refer to the document Using MPI on Carver (this document will be available in about a week)

The assignment includes a report, and the requirements for this report are described below. The report is worth 70% of your grade for this assignment.

The provided starter code is available via svn, as described on the following web page:

http://www-cse.ucsd.edu/classes/wi14/cse260-a/HW/A3/StarterCode.html

Extra credit opportunity

To encourage early progress, we are providing an opportunity for you to get extra credit by submitting a progress report. We will add up to 10% to your grade for the assignment if you can document, via a 1 page progress report, that you are getting parallel speedups on up to 8 cores. This is conditioned on our running your code and verifying the speedup you reported to us. Your progress report is due on Thursday Feb 27, midway through the assignment.

The accelerated algorithm

The algorithm used in this assignment takes advantage of knowledge about the underlying physical problem. In particular, because the force induced on a particle drops off rapidly with distance (see plot below), we may safely ignore forces at a distance greater than the cutoff, which is set to 0.01 in our simulations. This strategy significantly reduces the amount of computation.

The code sorts the particles, based on their coordinates, into a geographic data structure, a 2d array called a binning mesh. This mesh divides the simulation region into a uniform array of bins, which we’ll call world, containing all the particles residing in the corresponding region of space.

The force between two particles with the same y coordinate (i.e. on a horizontal line) is a function of their distance d:

Each bin contains a list of particles contained within the physical region of space covered by the bin. We organize the computation (apply_force) around this mesh and need only look to bins whose i and j coordinate differs by at most one. While the algorithm is still O(n²) in the worst case, where n is the number of particles, we have reduced the number of computations significantly since particles interact with nearest neighbor bins only (and including the bin itself). Put differently, given fixed time, we can solve much larger problem than if we didn’t use the mesh and had to compute distances between all pairs of particles.

The code allows you to specify M, the size of the binning mesh (binning meshes are square). As M shrinks, so does the the number of nearest neighbors to be checked. However, M must not drop below the cutoff value of 0.01 units or the computational will miss some interacting particles. In practice, there may be a large lower bound on M, since as M increases the bins grow as O(M²), and the a cost of handling the ever more numerous bins may cancel out some of the savings of reducing the number of distance checks.

Note that, as you increase the n parameter, the size of the computational region, as measured in floating point valued spatial units, grows as √n. This ensures that the density of particles remains constant as n increases (roughly speaking M/√n = constant.) Thus, you will need to increase M in order to optimize the running time, subject to the above constraint.

The provided code

To get you started, we have included some basic calls to MPI in the provided serial code. This code is enabled when you assert the mpi=1 make command line flag. Additional code is included to show you how to restrict output to one MPI process, for example, when printing information to standard out. This code is intended to illustrate the concepts and is by no means complete. You’ll have to make similar changes elsewhere.

We have also provided a batch script file, batch_serial.sh, to make serial runs, to run your code without MPI, and two strong scaling scripts, with MPI. Be sure to use the correct batch file for the capability you require; for example, the parallel batch job won't run correctly if you haven't enabled MPI when building yoru code. It will deadlock, tying up resources for the specified duration of the job.

Familiarize yourself with the various modules. The program defines two main classes to store the binning structure. The first class, called World, defines the 2D array, and is found in World.{cpp,h}. Each element of the array is an object of class Bin, containing a C++ vector of particles.

Class World defines the Simulate() member function which controls the overall simulation via two other member functions, World::move_particles() and World::apply_forces, that update the binning structure. In turn, these functions apply identically named member functions (but of class Bin) to each bin, i.e. Bin::move_particles() and Bin::apply_forces.

Your changes will be made to five (5) files:

World.cpp     World.h     Bin.cpp     Bin.h     nbody.cpp

If you need to support scalable buffering (groups of 3, described below), then you may also modify the plotting code. If not, you will not need to change that code.

To ensure that your program is correctly graded, do not rely on any changes to any other files besides the ones in the above list. We will substitute fresh copies of the remaining files when we grade your assignment. Also, don't change the signature of functions defined in particle.h as our test code makes callbacks to them.

To build your code, use the provided Makefile, which produces a target called nbody. We will build your code using your Makefile (with mpi=1 asserted, to enable MPI). You may modify the Makefile if you wish, however, do not change the name of make target as our tester depends on it.

Command line options

The code has various capabilities, which are controlled via command line flags explained in cmdLine.cpp. For example, the following invocation will run with a problem of size n=16384 particles for t=1000 timesteps on 32 processors with a process geometry of 8 x 4.

mpirun -np 4 ./nbody -n 16384 -t 1000 -px 8 -py 4 -M 40 -p 10 -s 1352723313 -o 16k.out

The particles will be organized in a 40 × 40 mesh binning structure. The results will be plotted every 10 timesteps. A random number seed is specified via the -s flag. (See below in the discussion of Testing your code.) We have asked to output the positions of the particles to a file called "16.out" This output may come in handy when debugging your code, in conjunction with summary output mentioned below.

Multicore executable

A multicore executable has also been provided, to be used in making performance comparisons on one node. You cannot use the same executable on both Bang and Carver; the executable on Bang is called nbody-bang, on Carver, nbody-carver. (The Carver executable will be made available soon)

To vary the number of threads, use the -nt flag, e.g. the following invocation will run on 4 threads, where X is one of bang or carver:

./nbody-X -M 80 -n 100000 -s 1359856190 -t 50 -nt 4

Testing your code

The provided code reports 2 accuracy measures that you can use to determine if your code is correct. These values reflect information about all the particles in the simulation and hence are summary errors. They may not catch small errors involving individual particles; for that you’ll need to output the particles’s positions.

Compare the two reported values against the single core run, obtained with the following command line parameters:

mpirun -np 1 ./nbody -M 40 -s 1359856190 -px 1 -py 1
...
(u,v) Max = ( 2.27641850e+00, 2.84668447e+00 )
(u,v) L2 = ( 5.81159812e-01, 5.81967427e-01 )

These values should agree to the full printing precision (9 significant figures) independent of the number of cores for smaller values of mesh size M. However, as you grow the -M parameter, you'll note changes in the last 3 digits. Such changes are due to roundoff errors that are a side effect of how the particles are organized in the chaining mesh. We will accept these errors, which should never affect the higher order 6 digits. Larger errors indicate an error. Note that, if you run for sufficiently long times (-t option) and/or with more particles, you may note larger errors, even in an otherwise "correct" code. These errors are an artifact of the computation so explore this regime with care as it is not guaranteed to be correct. For the default -t value chosen by the code, these errors should not arise.

Since the initial particle distribution is random, each run will use a different particle distribution, making it impossible to test for correctness unless we have a way of reproducing a given random number sequence. To this end, the appliation provides -s flag to start the random number sequence from a known place. For each combination of input variables such as -t and -n and processor geometry, run 3 times with the same seed and report the minimum time.

We will test your code's correctness and performance on Bang, and if you used it, on Carver. We will test against the public input described here, plus other inputs that will not be made public. Be sure to test against various combinations of the -px, -py, -n and -M command line parameters.

Like other parts of the application, the error reporting code also needs to be parallelized. The existing code will work on one core only as it assumes all the particles are contained within one process. Consult the code in particlesHelper.cpp which provides entries that you can use to compute the norm correctly on mulitple cores, without having to collect all the particles into one MPI process.

Plotting capability

As mentioned previously, the code includes with a plotting capability to observe the simulation dynamics which can be helpful when debugging. Plotting code should run on only one process and you will need to implement to move the particle data accordingly so that it can be plotted. When collecting performance data, do not not enable plotting (-p).

The plotter, implemented in Plotting.cpp, interfaces with gnuplot, which has been installed on bang (and is freeware). Only one process can control the interaction (via unix pipes) with gnuplot. So, when plotting on more than 1 core, you will need to communicate all the data to a distinguished processor in order to plot. The techniques will be discussed in class and entail collecting the solution onto one processor using an MPI Gather operation (see links to documentation below), taking care not to plot ghost cells. However, this solution isn't scalable, since all the particle data must be stored in one buffer residing in a single process. If working in a group of 3, you are required to implement a scalable solution, that requires a buffer of the size O(N/P), where N is the number of particles and P is the number of processes. In this case, you may modify the code in Plotting.cpp

Remote visual display can be slow, so either work from a campus machine, or increase the frequency to at least 20 (-p 20). The summary statistics may be more useful than plotting in some cases. If a display window does not appear on your screen, then your ssh display settings are not set correctly. If you are connecting to the remote machine from a Linux platform using ssh, use the -X option:

% ssh -X yourlogin@ccom-bang-login.ucsd.edu

If you are connecting from a Windows platform, and are using a 'secure shell client' AND 'exceed' (or xwin32, or cygwin) to log into a remote machine you must ensure that 'X11 tunneling' is set. To do this, you must have installed Cygwin /X with X11.

Due to an interaction between OpenMPI and Unix Pipes, we need to specify an additional command line to mpirun:

mpirun -np 4 --mca mpi_warn_on_fork 0 nbody -M 40 -s 1359856190 -px 1 -py 4

If you don't use the -mca command line argument, then you will get a large warning message beginning with the following text: "An MPI process has executed an operation involving a call to the "fork()" system call to create a child process." We will ignore this error but if you become aware of any problems, post to the A3 forum.

Alternatively you may set the following shell environment variable before running mpirun (e.g. via your shell startup scripts):

  export OMPI_MCA_mpi_warn_on_fork="0"

OR create a user-level configuration file $HOME/.openmpi/mca-params.conf containing the following single line:

  mpi_warn_on_fork = 0

How to Take timings

When you are measuring running times for your performance studies, do not enable plotting (-p) nor output the result to a file using the -o option, as doing so could affect your observations. Adjust the iteration count to keep the running times to within about 10-15 seconds.

Take timings using the provided timer code, which measures wall clock time. With MPI enabled, the timer uses MPI_Wtime() else gettimeofday(), as in A1. (Instructions for how to use the timer to take timings in MPI are described in the document Using MPI on Bang.)

On Bang, use an interactive qlogin node to debug your code and tune performance. (Interactive nodes are limited to 8 cores) Use the MPI batch queues to collect timings. Instructions for submitting jobs via the batch queue are described in

http://cseweb.ucsd.edu/classes/fa13/cse160-a/Bang/runningJobs.html

In performing your scaling studies, you need to take account of the fact that, the treatment and generation of random number sequences requires special care with MPI.

For example, when you require reproducible results while when debugging, you need to have a distinguished process generate all the random numbers required to produce the initial particle distribution, generate the particles, and distribute the particles to the other processes accordingly (i.e. with MPI_Scatter). This way, you can get reproducible results with the use of the -s flag. Note that if you let each process generate its own random distribution, and you specify the seed, the distribution will be identical across all processes and this will generate numerical errors.

The above method will work in your strong scaling studies. However, this not true for weak scaling studies, since the number of particles is varying with the number of cores. In this case you must average several runs for each different value of N, or any other parameter that affects performance (such as M). You should not use the -s flag. (There is a small likelihood that different processors will generate the same random number sequence, since they are seeded with the system clock. One could use a more ambitious random number generator, like SPRNG, see www.sprng.org but we will accept the small likelyhood.). This averaging process will smooth out noise due to different particle distributions. It is suggested that you take the average of 10 runs.

So, if you want to know if using a different value of the -M parameter affects performance, then for each distinct value of -M, you must perform 10 runs, take the average, and compare.

In sum.

• For weak scaling: you must take the averaging method.
• For strong scaling: the averaging method is preferred, but you may take the scattering approach This applies to both groups of 2 and 3 and it is not necessary for groups of 3 to initialize data in a scalable manner.
• To generate particles from the root, you will need to write a new init_particles() method, since you aren't allowed to change the existing code in particlesHelper.cpp. Just put the code in nbody.cpp and give the routine a different name so we can distinguish from the provided code.

MPI

For documentation and tutorials about MPI consult the following web page

http://www-cse.ucsd.edu/users/baden/Doc/mpi.html

Among this material we suggest the following

1. A tutorial published by LLNL with some worked examples.
2. Man pages for the MPI calls may be accessed with the man command and are also available On-line.

Assignment Checklist

1. Implement the MPI version of the provided code.
   1. Manage data partitioning so that the particles are distributed over all the processes according to the specified processor geometry. Correctly handle the case where the number of cores does not evenly divide the mesh size (-nx parameter). Implement 1D partitioning first, then 2D.
   2. Add point to point communication to manage ghost cells and to move particles that change owners.
   3. Parallelize the Plotting code. Groups of 2 can gather the data into a distinguised process and need not modify the plotting code. Groups of 3 must implement the scalable solution described in the Plotting Capability.
   4. Parallelize the code to compute the error checksums.
2. Conduct a strong scaling study, using the default parameters defined in the provided code. Be sure to use the same random seed on all your timing runs. If you have time, try different random seeds. You can let the program pick the random seed. In this case, you should run with at least 10 different seeds and average the timings obtained with each seed. Note that, for a given seed, you need to run at least 3 times, and take the minimum. Thus, if using K random seeds, you will make at least 3*K runs.
   1. Vary the number of cores from 1 to 8, in powers of two. Verify that results are independent of the number of cores.
   2. Determine the optimal value of -M on 1 core to within a resolution of 8, i.e. 24, 32, 40 etc. There is no need to explore this parameter at a finer resolution. Once you have established the value of -M on 1 core use it in your strong scaling study. Plot the parallel speedups and explain your results. You be able to improve performance further by tuning the -M parameter for each processor geometry. However, it is not necessary to conduct a two parameter study of -M and the number of threads over all the values of -M that you tested on one core.

3. Conduct a weak scaling study using the default parameters defined in the provided code. Run on 1 to 8 nodes on Bang (8 to 64 cores), maintaining a fixed workload per node. Groups of three should also conduct a strong weak scaling study on Carver, running on at least 256 cores (start at 32 cores.) It may be difficult to run on more than 256 cores, due to queue delays. Groups of 2 may run on Carver if they wish, but this is optional.
4. Report and comment on performance. The assignment includes a report, which should present and discuss the results of your experiments and your optimizations. The requirements for your report are described here, along with instructions on what to submit:

   http://www-cse.ucsd.edu/classes/wi14/cse260-a/HW/A3/WhatToTurnin.html
   

5. Groups of 3 have a somewhat more ambitious assignment.
   1. Do a weak scaling study on up to 256 cores of Carver.
   2. Implement scalable buffering for the plotting capability, as described in the section Plotting Capability.
   3. Explore the tradeoffs in hierarchical parallelism: MPI+threads. Each process will spawn multiple threads, which will carry out the computation using shared memory (processes will still communicate via messages). Compare performance with the pure MPI solution and offer insight in your report.
      Threading can be treated using OpenMP, pthreads, or C++ threads. MPI+threads programming will be discussed during the 2nd week of the assignment (there are also some changes needed in the batch submission file).

Afternotes

• In the weak scaling study, you'll need to optimize M for each value of N used. In generally, the optimal M value will grow with M.
• In the weak scaling study(ies), you'll need to choose the optimal geometry for each value of P. This is important when using more than 1 node.

Scott B. Baden [Wed Feb 12 12:12:07 PST 2014]