![[Mesh]](./FDMesh.png)
![[5-pt stencil]](./5-point.png)
| Date | Description |
| 20-Nov-13 | Original posting |
| 03-Dec-13 | Clarification about Initialization |
In this assignment you'll implement the heart electrophysiology simulator discussed in class using the Message Passing Interface, MPI, and run on the Bang cluster. Since implementing parallel codes under message passing can be more time consuming than under threads, give yourself sufficient time to complete the assignment.
You've been provided with serial code, written in C++, and an executable reference implementation to assess performance goals. (There are two, one for each course platform. Code compiled on Bang won't run on Ieng6, and vice versa).
The procedure for obtaining the starter code — and a working reference executable — is similar to the previous assignments, and can be found by visiting the following web page:
http://www-cse.ucsd.edu/classes/fa13/cse160-a/HW/A4/StarterCode.html
The assignment includes a report, and the requirements for this report are described below.
For instructions about how to compile and run MPI jobs on Bang, refer to the document Using MPI on Bang.
Simulations play an important role in science, medicine, and engineering. A cardiac electrophysiology simulator like the one we are using in this assignment can be used for clinical diagnostic and therapeutic purposes. To reduce computational costs, simplifying assumptions are often made to make the simulations feasible. In our clinical example, time may be of the essence in determining an appropriate life-saving treatment.
Cell simulators entail solving a coupled set of a equations: a system of Ordinary Differential Equations (ODEs) together with a Partial Differential Equation (PDE). If you are interested in learning more about cardiac modeling peruse www.labheart.org or thevirtualheart.org or visit the web page of the UCSD Cardiac Mechanics Research Group.
Our simulator models the propagation of electrical signals in the heart, and it incorporates a cell model describing the kinetics of the membrane of a single cell. There can be different cell models, of varying degrees of complexity. We will use a model known as the Aliev-Panfilov model, that maintains 2 state variables, and solves one PDE (The PDE couples multiple cells into a system). This simple model can account for complex behavior such as how spiral waves break up and form elaborate patterns. Spiral waves can lead to life threatening situations such as ventricular fibrillation.
Our simulator models electrical signals in an idealized biological system in which voltages vary at discrete points in time, called timesteps, at the discrete positions of a mesh. The simulator is an iterative procedure that updates the simulation state variables (on the mesh) many times over a sequence of discrete timesteps. A simulation will run for some amount of simulated time, which is a real number. (For our purposes, the units of time are unimportant.) The timesteps subdivide the simulation time line that starts at t=0 and ends at some specified time t=T. The unit of subdivision is called Δt, or the time discretization. It is important to keep in mind that when we use the term "time," we need to specify whether we are talking about real-valued simulation time, or the discrete timestep. We will use the fully qualified name to avoid ambiguity. By convention, we number the discrete timesteps with integers. Thus, if we say that we are at timestep k of the simulation, then the current simulation time is kΔt.
To discretize space, we'll use a uniformly spaced two-dimensional mesh, shown in the left side of the figure (Irregular meshes are also possible, but are more difficult to parallelize).
|
|
At each time step, the simulator updates the voltage according to nearest neighboring positions in space and time. This is done first in space, and next in time. Nearest neighbors in space are defined on the 4 Manhattan directions only: north, south, east and west, as shown on the right side of the above figure.
In general, the finer the mesh — and hence the more numerous the points — the more accurate the solution, but at an increased computational cost. In a similar fashion, the smaller the timestep, the more accurate the solution, but at a higher computational cost. Thus, if h is the mesh spacing, Δt is the time discretization, then the running time is proportional to (Δt×h2)-1. However, to simplify our performance studies, we will run our simulations for a given number of iterations rather than a given amount of simulated time.
The simulator maintains two state variables characterizing the cell electrophysiology model: excitation (E) and recovery (R). Since we are using the method of finite differences to solve the problem, we discretize the variables E and R by considering the values only at a regularly spaced set of discrete points. We represent the variables as 2D arrays: E[][] and R[][], respectively. The PDE solver also maintains the voltage at the previous timestep, Eprev.
Here is the formula for solving the PDE, where the superscripts t and t-1 refer to the current and previous timestep, respectively, and the constant α is defined in the simulator (see apf.cpp and solve.cpp):
Et[i,j] = Et-1[i,j] + α*(Et-1[i+1,j] + Et-1[i-1,j] + Et-1[i,j+1] + Et-1[i,j-1] - 4*Et-1[i,j])
Here are the formulas for solving the ODE:
E = E - Δt*(k*E * (E-a) * (E-1)+E * R); R = R + Δt*(ε+ μ1*R / ( E+μ2)) * (-R-k*E * (E-b-1));
The variables k, a, b, ε, μ1 and μ2 are constant floating point scalars, and are defined in the simulator. References to E and R correspond to whole arrays. Expressions involving whole arrays (that is, array names without subscripts) are pointwise operations on corresponding array elements.
We interpret a pointwise whole array operation as follows. If A, B and C are 2D arrays of size N × N, then A = B + C is equivalent to the following loop nest:
for (i = 0; i < N; i++)
for (j = 0; j < N; j++)
A[i][j] = B[i][j] + C[i][j];
When an array is combined with a constant, e.g. (E-a), the constant is extended to cover each element of the array, i.e. we compute the E[i][j] - a for all i and j.
To organize our work, we observe a number of bookkeeping conventions. We divide each axis into N-1 subintervals using N points positioned uniformly along the line at positions x=i×h, i=0, 1, 2, ..., N-1. We also pad the boundaries with a cell on each side, in order to properly handle updates next to a boundary. (This padding will simplify the code and you'll see comments about it in solve.cpp.) This is shown in below, for 1 axis of the mesh drawn above, where the mesh spacing h = 1/N. The points labeled 'o' correspond to the physical boundary, the points labeled 'P' are padding.
i=0 i=1 i=2 i=N-1 i=N i=N+1
P------o------|------| ...... |------|------|------o------P
^ ^
| |
x=0.0 x=h x=2h x=1.0-h x=1.0
Consider the case where we set N=64. We allocate an array of size 67 × 67, but we solve the problem on a central 65 × 65 submesh. Thus, our loops range from 1 to N+1, that is, 1 to 65. Points 0 and 66 correspond to buffer points.
It follows that when we subdivide the data to partition it for parallel execution, we are actually subdividing N+1 rather than N.
The provided serial code includes some basic calls to MPI to get you started, which are enabled when you assert the mpi=1 make command line flag. We have also provided a batch script file, batch_serial.sh, to make serial runs, to build your code without enabling MPI compilation, and batch_parallel.sh, to enable MPI compilation. 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. It will deadlock, tying up resources for the specified duration of the job.
Most (or all) of your changes will made to two files: apf.cpp, the driver program, and solve.cpp, containing the solve() function. The solver is computationally intensive, so the focus of your parallelization and performance programming effort will be on this routine.
To build your code, use the provided Makefile, as it sets up the appropriate flags for the compiler via the "arch" file located in $(PUB)/Arch/arch.gnu-4.7_c++11.generic. The Makefile produces a target called apf. We will build your code from the provided Makefile (with mpi=1 asserted) so don't depend on changes to this file. In particular, don't change the name of make target nor add additional code files. Moreover, to ensure that your program is correctly graded, do not change any of the following files:
cmdLine.cpp Report.cpp Timer.cpp util.cpp apf.h
The code and Makefile have been set up to enable you to compile with or without MPI, which may come in handy when debugging. The default is to compile without MPI. To compile with MPI, specify following make command line argument:
mpi=1
This setting defines the _MPI_ macro which causes certain code to be conditionally compiled. To take advantage of this feature in any code you write, follow the convention used in the provided starter code.
The simulator has various options, controlled via command line arguments (parsed in in cmdLine.cpp), with default values set in apf.cpp:
| -n <int> | Number of mesh points in the x and y dimensions |
| -i <int> | Duration of simulation in iterations |
| -s <int> | Print statistics, L2 and L∞ norms of the
excitation variable at specified regular intervals (in timesteps).
This feature is useful when debugging. |
| -p <int> | Plot the excitation variable, at the specified regular intervals (in timesteps). Useful when debugging. |
| -x <int> -y <int> |
x and y-axis of the the processor geometry (Used only for your MPI implementation) |
| -k | Disable communication (Used only for your MPI implementation) |
Thus, the following command line will conduct the simulation on a 150 × 150 box, run for 400 iterations, plot the excitation variable every 20 iterations, and print out summary values of the excitation variable every 10 iterations.
./apf -n 150 -i 400 -s 10 -p 20
If we leave off the -s and -p flags, then the output is as follows. The last line of output, beginning with the string "@ " outputs performance data along with relevant parameters, that is suitable for automated extraction from a sequence of runs.
[Simulation begins] Thu Oct 31 17:09:56 2013
dt= 0.0636729, # iters = 400
m x n = 151 x 151
processor geometry: 1 x 1
Compiled with MPI Enabled
End at iteration: 400
[Simulation completes] Thu Oct 31 17:09:56 2013
End at iteration 400
Max norm: 9.928177e-01, L2norm: 6.128517e-01
Running Time: 0.220473 sec.
GFlop rate: 1.17
M x N, [px x py] [#iter] {T_p, Gflops} Linf, L2
@ 151 151 [1 1] [ 400] {0.2205 1.174} 9.92818e-01 6.12852e-01
The plotting capability, controlled via the -p flag, uses gnuplot, and may be useful in debugging. Do not enable plotting when you are measuring running times for your performance studies. The plotter capability only works on a single core and is discussed in more detail later in this writeup.
You'll notice that the provided code uses a different naming convention for the bounds of the computational box. In particular, while we use the single letter N in this document, the code uses different bounds, for the x and y axes, in terms of variables m and n. This convention enables the code to accommodate non-square local problem sizes that arise when the data has been partitioned for parallel execution. The computational boxes have been padded by one cell in each dimension. This change is needed to handle boundary conditions, which are described in solve()
Also note the conventions used to index the 2D solution arrays: the indices are swapped, with the i index being the most rapidly varying index, that is, indexing along the last dimension of the array. Thus, array elements in the same row have different i indices and elements i-1, i and i+1 are adjacent in memory.
Parallelization under message passing is actually cleaner than under threads, since all data dependencies must be met through explicit data motion where both sender and receiver are aware of one another. Thus, MPI programs don't need critical sections, and only infrequently use barrier synchronization. The MPI message passing calls, including collectives, are usually sufficient to meet synchronization needs. However, the benefit of having to specify data motion in such detail comes at the cost of complicating your code. For these reasons, MPI can have a steep learning curve, but once code has been written, it can be easier to debug. For example, if your program hangs, you know that there is a missing Receive or send call, or somehow the calls don't match up, e.g. the calls used different message tags. These errors are usually much easier to spot than missing or incorrect synchronization.
If you run the single processor code with MPI compilation enabled, the program will run the identical computation on every processors. To parallelize the code, you'll have to decompose (partition) the 2D computational arrays and manage ghost cells as discussed in class. We advise that you begin with 1 dimensional (strip partitioning) along the slowest varying (first) dimension. Under these conditions, all ghost cell data are contiguous, and no intermediate buffering is needed. Once you have gotten one-dimensional partitioning to work, next implement two-dimensional partitioning. To do this, you'll need to buffer non-contiguous data prior to sending it, and unpack incoming data into the ghost cells since the originate from non-contiguous locations in the sending process.
You'll also need to add collective communication to properly compute the global L2 and L∞ values. These values need be printed from a single process only. You may use whichever process you like, though traditionally process zero is used. To this end, it makes sense to define a variable called ROOT, which you can use to designated the distinguished process that handles terminal I/O. It also comes in handy when in collective communication.
Initialization
Initialization requires special treatment, and should be parallelized differently from the solver.
E_prev and R should be initialized in global coordinates, namely, that you set the left half of E_prev (and the lower half of R), consistently with how you set up the array up on 1 core. It is not correct to use the local bounds to determine where to set values to 1. Rather, you should loop over your local bounds, and convert i,j to global indices (in the sense of the 1 core run). If you haven't already, implement the -p plotting option to help you see what is going on. Compare the output of the 1 core run with the multicore run and, if things are implemented correctly, they will look identical. But not if you used local coordinates to determine where to set values to 1.0.
Try to make your code go as fast as you can. In the spirit of friendly competition, post your results to the A4 Moodle Forum.
Your code will be tested for correctness and performance against various combinations of the N, niter and the number of processors and processor geometries.
Test your code by comparing against the provided implementation. You can catch obvious errors using the graphical output (see below), but a more precise test is to compare two summary quantities reported by the simulator: the L2 norm and the L∞ norm, the maximum value (magnitude) of the excitation variable. The former value is obtained by summing the squares of all solution values on the mesh, dividing by the number of points in the mesh (called normalization), and taking the square root of the result. For the latter measure, we take the "absolute max," that is the maximum of the absolute value of the solution, taken over all points. While there may be slight variations in the last digits of the reported quantities, variations in the first few digits indicate an error.
Verify that results are independent of the number of cores and the processor geometry, by examining the summary quantities. For example, for the command line options
-n 50 -i 1000
should result in the simulator reporting the following summary quantities at the end of the simulation:
Max norm: 9.866452e-01, L2norm: 5.919411e-01
To facilitate testing, the code outputs certain information to a file called Log.txt: the number of threads, their geometry, the running time, and the summary measures. Any previous copy of Log.txt will be overwritten, so if you are doing automated testing, be sure to rename the files between invocations.
As mentioned previously, the code includes with a plotting capability to observe the simulation dynamics which can be helpful when debugging. The plotter (Plotting.cpp) interfaces with gnuplot, which has been installedon bang (and is freeware). Note that two entries are provided, one to the provided multidimensional array implementation, and another to support flattened arrays.
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.
The simulations will produce a spiral pattern that evolves with time. However, if you set the number of iterations (-i) to a sufficiently large value, the organized spiral pattern will break up. This indicates non-physical behavior, and the length of time that the simulator produces physically meaningful results depends on the mesh size N.
Be sure that, while collecting performance data, that you have not enabled plotting (-p) or statistics reporting (-s).
Note that the plotter works on 1 core only, so if you want to use the plotting capability on more than 1 core, you will first need to gather all the data onto one processor before plotting. The techniques will be discussed in class an entail collecting the solution onto one processor using an MPI Gather operation (see links to documentation below), taking care not to plot ghost cells.
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 apf -n 100 -t 400
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 A4 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
Adjust the iteration count to keep the running times to within about 10-15 seconds. Use the mpi.q batch queue on Bang to collect timings. Use the a qlogin node to debug your code and to tune performance. (Interactive nodes are limited to 8 cores)
To take timings, use MPI's timer function MPI_Wtime() to measure wall clock time. MPI_Wtime() reports time in seconds and its man-page is available using the man command. To measure time, take the difference in times sampled by two successive calls to MPI_Wtime():
double t0 = MPI_Wtime();
// some code....
double t1 = MPI_Wtime();
double time = (t1-t0);
Recall that the running time of the application is that of the last processor to finish. Don't time individual iterations, but rather, time an entire run, but excluding initialization, calls to MPI_Init() and MPI_Finalize() and any other activities that are not representative of a run with many iterations.
Although Bang provides exclusive access to allocated processors, the Infiniband interconnect is a shared resource. Thus, if you run on one node, you'll observe no contention from other jobs, but on more than one node To reduce the effects of switch contention run the program at least 3 times (for each input size and number of cores, or any variable that affects the running time) until you see a consistent timing emerge, such that a majority of the timings agree to within 10%. If you observe greater variation, then try doubling the number of runs to see if you can obtain a better consensus. Exclude outlier timings, but note them in your report. Report the minimal running time rather than the average, as previously discussed in class.
The instructions for turning in this assignment are the same as for the previous assignment. For remaining instructions on what and how to turnin, visit the following web page
http://www-cse.ucsd.edu/classes/fa13/cse160-a/HW/A4/WhatToTurnin.html
Your grade will depend on 3 factors: performance, correctness and the depth and your explanations of observed performance and how the design of your code influenced performance. In determining your grade for the assignment, explanations and insight will be emphasized more heavily than performance itself, provided you were able to get your code to scale.
Ieng6 should be used for development only. Rather than using ieng6.ucsd.edu, use one of the 15 named servers ieng6-2NN where NN as an integer between 40 and 54 inclusive, i.e. ieng6-240 through ieng6-254. In an effort to evenly distribute the class over the servers, choose a server based on the following formula:
YourServer = Your AID MOD 12
That is, the remainder after dividing your AID number by 12. While this strategy isn't perfect, it will help. Don't use more than 4 MPI processes as there is likely no advantage to doing so.
For documentation and tutorials about MPI see this web page at: http://www-cse.ucsd.edu/users/baden/Doc/mpi.html. Among this material we suggest the following.
[Mon Nov 18 07:51:53 PST 2013]