Programming Lab 3

Cardiac Electrophysiology Simulation

Changelog

Contents

• Cardiac Electrophysiology Simulation
• The provided starter code
• Testing your code
• How to take timings
• Assignment Checklist
• Turning in your assignment
• Grading
• Frequently asked questions

Part 0. The Assignment

In this assignment you'll develop a multithreaded implementation of the heart electrophysiology simulator discussed in class and conduct experiments to assess performance. The assignment includes a brief report, and the requirements for your report are described at the end of this writeup.

Your grade will based on 4 factors: performance attained by your implementation; correctness; well-displayed results; descriptive comments in the code and in your repository.

In the spirit of friendly competition, post your performance results to the HW3 Piazza folder.

☞ Due Friday2/5/16 @ 9PM

‣ ‣ ‣ Required solution files: solve.cpp helper.cpp

This means that these files, containing your solution to the assignment, must be turned in by the deadline.

You may develop your code anywhere, but you must ensure that it compiles and runs correctly on Bang because that’s where we’ll be testing it. (If you'd like to develop your software on another machine, like a laptop, see these instructions in the FAQ section below.

In doing this assignment, it's your responsibility to understand the course rules for integrity of scholarship.

Part 1. Getting Started with GitLab

All assignments will be turned in using UCSD's GitLab. The turnin procedure is described below in the section "Turning in your assignment". You must follow these steps to get the starter code and to get set up to do assignment 3.

Refer to the link below for instructions on how to obtain the starter code to get started on your assignment

A3 Git Instructions

Part 2. Cardiac Electrophysiology Simulation

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 and the Simula Computational Cardiac Modeling Department

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 into timesteps. The time line starts at t=0 and ends at some specified maximum 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 or discrete simulation time. We will use the fully qualified name e.g. "discrete simulation time" to avoid ambiguity. By convention, we number the discrete timesteps on the simulation time line 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

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 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×h²)^-1. However, to simplify our performance studies, we will run our simulations for a given number of timesteps (iterations) rather than a given amount of simulated time.

The simulator algorithm

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. (our notation will use a pseudo C syntax, E[i,j] is synonymous with E[i][j]) The PDE solver also maintains the voltage at the previous timestep, E_prev.

Here is the formula for solving the PDE (see solve.cpp), where the superscripts t and t-1 refer to the current and previous timestep, respectively (α is defined in apf.cpp):

E^t[i,j] = E^t-1[i,j] + α*(E^t-1[i+1,j] + E^t-1[i-1,j] + E^t-1[i,j+1] + E^t-1[i,j-1] - 4*E^t-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*(ε+ μ₁*R / ( E+μ₂)) * (-R-k*E * (E-b-1));

The variables k, a, b, ε, μ₁ and μ₂ 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. We interpret a pointwise 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, where A, B and C have the same size, N:

for (i = 0; i < N; i++)
    for (j = 0; j < N; j++)
        A[i][j] = B[i][j] + C[i][j];

Pointwise operations operations are parallelizable, as there is no possibility of any dependencies. They extend to mixed array expressions involving scalars. When whole array is combined with a scalar, e.g. (E-a), the scalar 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 the following 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. In order to properly handle updates next to a boundary, we pad the boundaries with a cell on each side. This padding will simplify the code (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-1). 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 66 × 66, but we solve the problem on a central 64 × 64 submesh. Thus, our loops range from 1 to N, that is, 1 to 64. Points 0 and 65 correspond to padding points.

We've provided you with starter code that implements the Aliev-Panfilov solver as well as an executable so that you may compare your results with our reference. The file is found in $PUB/HW/A3/apf-ref on Bang. It may be possible to improve on the performance of the reference, since it doesn't include any significant optimizations to conserve memory locality.

All of your changes will made to just two files:

solve.cpp       helper.cpp

Solve.cpp contains the solver, helper.cpp provides various helper functions such as storage allocation. (You may therefore change the storage allocation scheme if you wish.) The solver is computationally intensive, and so the focus of your parallelization and performance programming effort will be on this routine.

As with the previous assignments, build your code using the provided Makefile, which sets the appropriate flags for the compiler via the "arch" file located in $(PUB)/Arch/arch.gnu-c++11.generic. The Makefile produces a target called apf. Do not change the name of this target. While you may modify any files not listed above, remember that we will substitute fresh copies when we grade your assignment, so don't really on any custom changes as they will be lost during the grading process.

We have also provided one batch script, batch_job.sh. This script runs the serial provided code, and you'll need to modify it to run your multithreaded code.

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 separate bounds for the x and y axes: the variables m and n. We need two variables so handle the rectangular subdomains that result from partitioning the data over multiple threads. You'll noticed that array indices are one-dimensional. We've "flattened" the logical two-dimensional arrays so that they are subscripted with just 1 index.

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. Be careful to use the correct convention when writing flattened subscript expressions.

The simulator has various options, controlled via command line arguments (parsed in in cmdLine.cpp), with default values set in apf.cpp:

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, and you can use this marker to extract performance data from an output file recording many runs.

$ ./apf -n 150 -i 400

[Simulation begins] Sat Feb  6 22:50:45 2016

dt= 0.0639292, # iters = 400
m x n = 150 x 150
# threads: 1

[Simulation ends] Sat Feb  6 22:50:45 2016

End at iteration 400
Max norm: 9.92837e-01, L2norm: 6.12580e-01
Running Time: 0.216509 sec. [1.16 GFlop/sec]

   M x N    NT    #iter    T_p     Gflops        Linf, L2
@ 150 150    1      400  0.2165    1.164    9.92837e-01 6.12580e-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 Part 4.

Your code will be tested for correctness and performance against various combinations of the N, niter and the number of threads. Test your code by comparing against the provided implementation. While you can catch some obvious errors using the graphical output (see below), a more precise test is to compare two summary quantities reported by the simulator: the L₂ norm and the L_∞ norm, the maximum value (magnitude) of the excitation variable, which are included in the 1 line of performance output. The first 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 second measurement, we take the "absolute max:" the maximum of the absolute value of each of the elements in the solution array. While there may be slight variations in the last digits of the reported quantities, variations in the first few digits indicate an error.

For example, you may find that after adding a new performance optimization that your code indeed ran faster but now the L2 errors differ in the 3rd decimal place, while the Max errors are the same. There is an error in your code. The L2 norm is giving a clearer picture of what is going on here; the maximum difference between the two solutions might not change, but the sum of the squares of their differences is another matter. Note, however, that the reverse situation could occur. The max norm may be the more useful measure. In this case your code is probably introducing small errors in many elements of the solution mesh, and a small number of larger errors. These errors aren't visible in the L2 norm, which takes an average. But the errors are exposed by the max norm. We give you both measures as they can each catch different kinds of errors. So be sure to check that each matches.

Plotting capability

The code includes a plotting capability to enable you to visualize obvious errors (The summary statistics will be more useful than plotting for catching subtle errors). 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. If you see this behavior occur for short lengths of time, or other non-physical behavior, you've located an error.

The plotter (Plotting.cpp) interfaces with gnuplot, as in Assignment 1. Remote visual display can be slow, so either work from a campus machine, or increase the frequency to at least 20 (-p 20). 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.

Part 5. How to take timings

As in the previous assignments, use either an interactive node (via qlogin) or the non-interactive batch queue to collect timings.

Take timings using the provided timer code, which measures wall clock time via gettimeofday(), as in the previous assignments. Adjust the iteration count so that the running time never drops below 1 or 2 seconds.

When you are measuring running times for your performance studies, do not enable plotting (-p) nor output statistics using the -s option, as doing so could affect your observations.

Part 6. Assignment Checklist

1. Parallelize your code using C++11 threads. Partition the data in one direction only, corresponding to the 2D index that varies the slowest. (Look at how the provided code indexes the arrays to work this out.)
   1. While you need support square global meshes only, your code must support rectangular local meshes that each thread will work with.
   2. Your code must correctly handle the case when the number of threads geometry doesn't divide the mesh evenly, including odd values of the number of threads and N. In particular, it must divide the mesh as evenly as possible using partitionings having the following property:
      

      The amount of work assigned to the most heavily and most lightly-loaded process along a given dimension must differ by not more than one row of the solution array.

   3. If you find small errors (4th or 5th decimal place), the error is small. We advise that you hold off on fixing this subtle error, which may be an off-by-one indexing coding bug. Instead, assess performance and carry out and needed performance optimizations (Step 2), culminating with a scalability study (Step 3). Once you've done this assessment, and if you have the time, try to locate the source of the error. If not, write up your results and take some time time to answer the following 3 questions, which you should add as addendum to your scaling study plot:
   1. Where is the error coming from? Make your best guess by examining the code. Show us code fragments as needed to clarify
   2. How severe is it? Try successively doubling the number of iterations with the -i flag, say for 2000, 4000 and 8000 iterations and see how much the errors grows.
   3. How much of an increase in performance did you get at a cost of introducing this error? Was it worth it? If you got say only 10% performance improvement, probably not worth it. But if you went much faster, could be worth it.
2. Optimize your code to improve performance. To help you assess your performance, we've provided a reference binary executable. You may be able to beat our reference.
   1. Your implementation should scale strongly: for a sufficiently large problem, performance should be proportional to the number of cores. We suggest you take this in steps as follows
      1. Determine if your code can scale on 2 cores. Set N=400 and # iterations (niter)=500; the provided code will complete in about 2 seconds on one core. On two cores, the running time should be about 1/2 that of the single core run. If you are observing this behavior, move to the next step.
      2. Next, set niter=2000. Conduct a strong scaling study: observe the running time as you successively double the number of cores, starting at p=1 core and ending at 8 cores, while keeping N fixed. Your code should run in under 10 seconds on 1 core of Bang and run in about 1 second on 8 cores.
      3. Compare the single processor performance of your multithreaded implementation against the performance of the provided code to asses any overheads introduced in parallelizing the code. These overheads should be negligible; if not, investigate the source of the overhead and eliminate it.
      4. Optimize the inner loops to improve cache locality. Do not use SSE - this will come in assignment #4.
   2. For larger N--a few thousand--the code won't scale as well as it did for smaller N. What is the largest value of N that will enable your code to scale perfectly to 8 cores? After completing the strong scaling study for n=400, try larger values of n, successively doubling N to 800, 1600 and so on. You will need to decrease the number of iterations since the running time grows as the square of N. Keep your runs on 8 cores to about 1 second.
3. Report and comment on your results. The general requirements for the document are described in §7.2.
   1. Present the result of your scaling study both as a plot and as a table. Also provide a plot showing the speedup on 8 cores as a function of N, starting with N=400. Double N until you begin to notice that the speedup drops below 7. What is the largest N that will enable you to obtain a speedup of at least 7?
   2. If you have negative results, document them (see item (1.c) of the checklist) as we way be able to award partial credit if you can offer insight into program errors or performance bottlenecks that prevented you from collecting all your data.

Part 7. Turning in your assignment

7.1 Submission via GitLab

This procedure is the same as the previous assignment and is reproduced here:

Git Lab Submission

1. A PDF document that presents all requested performance data along with any addendum explaining negative results. Use the same format as in the previous assignments, e.g. A3-report-baden-powell.pdf. This report is worth 20% of your grade.

   Provide all data in both tabular and graphical form. Plot your data carefully, plotting only what's necessary to get your points across, and choosing log and linear plots carefully. Be sure to provide a caption for each table or plot that describes the experiment you conducted, e.g. "strong scaling."

   We encourage you to report negative results which are also valid especially valued if you can document them, with the underlying causes, or best guesses if the causes are not clear. So if results weren't as you expected, let us know why and we may be able to give you some partial credit. Include this discussion on the PDF containing your results.

2. Source code. Even though we will be substituting fresh copies of all but one source module, include all code files your worked with, in case we need to unravel a problem with your turnin. This could happen because you made changes to one of the files that we replaced. Also include your Makefile.
3. Team evaluation. Your turnin must must also contain a team evaluation, the file teameval.txt. This document must be filled in completely. If you don't turn in it in, your assignment will be considered late. (Minimum 10% lateness penalty).

There are 25 possible points on the assignment and your grade will depend on 4 factors: performance attained by your code, correctness, well-displayed results (including an explanation if you have negative results) and comments. You must not only put comments inside your code, but every time you commit to your repository. There should be a steady stream of comments that enable us to retrace your development steps. Your performance score will be based on the ability to scale strongly.

We will grade your code on Bang using our own Makefile with the same compilation flags as in the provided in starter code. To get full credit, design, style and comments (including gitlab commit comments) count.

1. Q. How can I use the statistics output (L2 and Linfo) to help me debug my code?

   A. Your first step after partitioning the data structures and setting up the initial conditions is to make sure you get the statistics reporting to work, so that you will have diagnostic tests enabling you to compare your code with the provided reference. Having done this, you can test that you've set up the initial conditions correctly, and that your solver runs correctly on the local partitions.

2. Q. What performance are you looking for?

   A. We are looking for perfect scaling with the public test. Follow the guidelines in the check list. For larger N (a few thousand), performance won't scale as nearly well. At some point you may reach a plateau where adding more processors doesn't help much, or not at all.

3. Q. We think we have subtle errors in the solution, that the visual display can't help us with.

   A. To see what the solution, modify the provided routine printMat() [helper.cpp] so it prints out larger matrices. Redirect to a file. Or, if you like, shut off the test in apf.cpp so that the code will run when N <= 25 and use a small matrix like N=8 or 9. This output is easily readable at the terminal. (Remember to turn back on that check when done since the computation will fail when N <= 25.)

4. Q. How I can get to develop the code on my laptop? Do I need to modify the Makefile?

   A. Copy over the "arch" file that the Makefile includes, and modify the path in the makefile to point at the place where you put the arch file. Make sure that you have g++ version 4.8.4. Earlier versions might also work (but not early than 4.7, which supports the C++11 standard) but not guaranteed. Before you turn in your code, be sure that it runs correctly on Bang and meets the performance goals on that machine If you are using an IDE you might need to add pthreads lib, as C++11 threads relies on pthreads. In some environments, this information may be needed.

5. Q. How can I get assembler output?.

   A. Add asm=1 to the make command line as in

   make asm=1 solve.o
   

   You'll get a file called solve.cpp.lst, containing the assembler and the source code, with the assembler interleaved among the lines of source. This file is easier to read than the raw assembler output, but having said that, the optimizer makes life "interesting." Try turning the optimizer off to (-O0) the first time you try this out. To get the raw assembler output just do "make solve.s" For further information this blog: https://panthema.net/2013/0124-GCC-Output-Assembler-Code/ or this web site, http://gcc.godbolt.org/ that provides a graphical interface.

6. Q. We have errors in the 3rd decimal place of the reported L2 or Linf values.

   A. Here are a few tips.

   1. Make sure that the code that computes L2 and Linf values is correct.
   2. Choose N such that N+2 is divisible by the # of threads. This will eliminate certain possible sources of errors. (Recall that the simulator allocates an array of size N+2)
   3. Make liberal use of barriers. Don't be concerned about performance, unless you are debugging a performance optimization. You can remove those barrsier one by one to improve your understanding of of whether or not your code handles synchronization properly
   4. Use the plot -p option. While the plotter code itself is single threaded, the point of using this option is to help you determine that your multithreaded code is correct.
   5. Make sure to properly handle synchronization on code that must run on 1 thread, such as the plotter.
   6. Use the -s option with a value of 1 to see if the error begins immediately or whether it takes time to appear. If the error appears immediately, it's more likely you don't have a simple "off by one error", then if the error takes several iterations to appear.
7. Q. For N=512, we have obtained good spedups (greater than 4) on 8 cores. But our speedup on 2 cores is still a little lower than 2. We have tried iterations = 500, 1000, and this is a steady result. Is this okay, or do we need to realize at least a speedup of 2 on 2 threads?

   A. Compare with the reference and satisfy yourself that, for various values of N, you are doing about as well. We interpret "about as well" as "at least 90% of the performance" Remember that we do have some tests that won't be made public, so try a few different values of N, and be sure to do the study to determine the threshold for realizing scalability, such that performance is within 90% of the reference. It isn't about speedup alone; you need to get within 10% of the running time of the reference. Matching the speedup realized by the reference is not sufficient: your code could actually run slower (though it probably won't)

[Mon Feb 8 21:24:27 PST 2016]