CS 160 Homework Assignment #4
Cardiac Cell Electrophysiology Simulation with Message Passing

Changelog

Introduction

In this assignment you'll parallelize, using MPI, a simulator used to study the electrophysiology of the heart. A working serial version of the code has been provided as the gzipped tarball sim.tgz. Try to make the code go as fast as you can; thus, you should not only optimize message passing performance, but also single processor performance.

Cardiac Electrophysiology Simulation

Simulations play an important role in science, medicine, and engineering. For example, a cardiac electrophysiology simulator can be used for clinical diagnostic and therapeutic purposes. To reduce computational costs, simplifying assumptions are often made to make the simulations practical. 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 (called ODEs for short) together with Partial Differential Equations (called PDEs for sort). We will not be focusing on the underlying numerical issues, but 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 a the membrane of a single cell. The PDE couples multiple cells into a system. There can be different cell models, with varying degrees of complexity. We will use a model known as the Aliev-Panfilov model, that maintains 2 state variables, and solves one PDE. 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 system in which voltages vary at discrete points in time, called timesteps on discrete positions of a mesh of points. In our case we’ll use a uniformly spaced mesh (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.

In general, the finer the mesh (and hence the more numerous the points) the more accurate the solution, but for 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, and Δt is the discrete timestep, then the running time is proportional to (Δt×h²)^-1.

The simulator keeps track of two state variables that characterize the electrophysiology we are simulating. Both are represented as 2D arrays. The first variable, called the excitation, is stored in the E[][] array. The second variable, called the recovery variable, is stored in the R[][] array. Lastly, we store E_prev, the voltage at the previous timestep. We need this to advance the voltage over time.  

Discretization (3/4/10)

Since we are using the method of finite differences to solve the problem, we discretize the variables E and E_prev by considering the values only at a regularly spaced set of discrete points. For example, if we look at the decomposition along one axis, then we are divide the solution interval into N-1 subintervals using N points positioned uniformly along the line at positions x=i×h, i=0, 1, 2, ..., N-1. This is shown in the following figure, where the interval spacing h = 1/N and is known as the mesh spacing

i=0   i=1   i=2                      i=N-1
|-----|-----|------| ..... |-----|-----|
^                                      ^
|                                      |
x=0   x=h  x=2h                 x=1-h  x=1

Consider the case where we set n=51. We allocate an array of size 54 × 54, but we solve the problem on a central 52 × 52 submesh. Thus, our loops range from 1 to n+1, that is, 1 to 52. Points 0 and 53 correspond to buffer points.

It follows that when we subdivide the data to partition it, we are actually are subdividing n+1 rather than n. Thus, we choose a value of n such that n+1 is divisible by the number of threads. It's fine if you work with 2, 4 and 8 threads only.

Serial Code

We are providing you with a working serial simulator that uses the Aliev-Panfilov cell model. The simulator includes a plotting capability (using gnuplot) which you can use to debug your code, and also to observe the simulation dynamics. The plot frequency is hardwired into the code but you might want to adjust it. Your timings results will be taken when the plotting is disabled. The simulator has various options as follows.

With the arguments interpreted as follows (The default parameters are described in cmdLine.C):

Thus, the command line

will simulate on a 50 × 50 box, and run to 1000 units of simulated time, plotting the evolving solution every 100 units of simulated time and periodically printing out simulation statistics.

Additional arguments have also been provided for the parallel version of the code: -x and -y for specifying the processor geometry, and -k for shutting off communication. These options will be described below.

NOTE: If a display window does not appear on your screen, then your SSH display settings are not set correctly. If you are using a Linux based machine and are connecting to a Linux machine using SSH, use the -X option:

If you are using a Windows based machine, and are using a 'secure shell client' AND 'exceed' (or xwin32, or cygwin) to log into Linux machine, you must ensure that 'X11 tunneling' is set.

The Assignment

In this assignment you'll make modifications to the provided code. In particular, you'll parallelize it with MPI. You'll also conduct strong scaling studies. With a fixed problem size N run on 1, 2, 4 , and 8 cores. What speedups do you observe? Discuss your observations and explore the underlying causes of what you’ve observed. If you see perfect speedup, shrink the problem size until you see differently, and explain your results. If you see poor speedups, grow the problem size until the speedups improve. Also perform the following study on 8 cores. Determine the smallest value of N for which you obtain linear speedups.

Three command line arguments have been provided for analyzing performance. The first two, -x, -y will enable you to adjust the processor geometry. You should use these arguments to partition the data with either horizontally (-x 1 -y < p >) and vertically (-x < p > -y 1) oriented partitions. Measure the performance on 4 processors using both horizontally and vertically oriented partitions and note the difference in performance. To further explore the cause of the performance difference, report the communication time.

It can be difficult to measure communication times accurately, owing to the short time durations involved. Therefore, we will use an indirect measurement technique. We do this by disabling communication of ghost cells (using the -k option) and subtracting the resultant time from that obtained from a run that includes communication.

Note that, when plotting data, you may need to collect all the data on the root processor before plotting. Though this implementation strategy isn’t scalable, we won’t be running on enough processors for this issue to be a concern.

Be sure that, when collecting performance data, that you have not enabled plotting (-p) or statistics reporting (-s) when collecting performance data.

The use of 2D arrays may or may not be the most efficient way to implement the simulator. As in the image morphing assignment, you may need to reorganize the data structures. More generally, try to make the code go as fast as you can.

How to test your code

./sim -n 100 -t 2000 > a
dt=7.648953e-02, T=2.000000e+03
End at time: 2000, iteration 26148
Max:  9.760256e-01, L2norm:  6.351231e-01
Wallclock time: 14.606487 sec.

Performance Measurement

To take timings, use MPI's timer function MPI_Wtime() to measure wall clock time (like a stopwatch). Your timings will reflect contention from other jobs. In a true dedicated run, you would be the only user on the machine, and contention would be minimal; running times would be shorter and reproducible. Since it is unrealistic to exclude other users, you’ll use a different technique, that will filter out variations. Run the program at least 3 times (for each input size and number of processors, or any variable that affects the running time) until you see a consistent timing emerge, such that a majority agree within 10%. If you observe a lot more variation try doubling the number of runs to see if you can obtain a better consensus. Exclude outlying timings, but note them in your report. Report the minimal running time rather than the average, as previously discussed in class and in the text.

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);

The man-page for MPI_Wtime() is found on line. Note, however, that all times are reported in seconds, and there is no need to divide by MPI_Wtick() to convert to seconds. (You can obtain the resolution of MPI_Wtime using MPI_Wtick())

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.

Performance Guidelines

Building and Running MPI Programs

mpirun -np 4 ring -t 10000 -s 16

To establish that your environment has been set up correctly, copy, compile and run the parallel "hello world" program. This program prints "Hello World'' from each process along with the process ID. It also reports the total number of processes in the run. The hello world program is found in $(pub)/examples/MPI/Basic. To run the program use mpirun as follows:

mpirun -np 2 hello

Here is some sample output:

# processes: 2
Hello world from node 0
Hello world from node 1

MPI Documentation

Here is a web page devoted to MPI documentation Man pages for the MPI calls may be accessed with the man command. More extensive documentation can be found at at http://www-cse.ucsd.edu/users/baden/Doc/mpi.html. You can obtain man pages online for MPI calls used in the example programs described here:

Things you should turn in

Document your work in a well-written report which discusses your findings carefully. This part of your assignment will be worth 50% the grade. Your report should present a clear evaluation of the design of your code, including bottlenecks of the implementation, and describe any special coding or tuned parameters. Provide pseudo code listings that illustrates how you parallelized the code.

Your report should cite any written works or software as appropriate in the text of the report. Citations should include any software you used that was written by someone else, or that was written by you for purposes other than this class.

Transmit the report electronically along with your source code and your report. by the deadline. Be sure and delete all .o files and executables before creating the tar file. Transmit a single archive file, e.g. .zip, .tgz, containing your writeup, and source code, plus any output files.

Email the archive file as an attachment to the following address:

The name of your tar file should start with your login ID, e.g. baden.tgz. Your archive file should create a directory with the name of your login ID. The directory should contain two subdirectories

src:     All source code including your makefile, so we can build the software.
report: Your report including any separate image files, screen shots, or other output.