Overview

For this project we will we will be writing a program to solve Laplace's equation.

This equation can be used to describe heat transfer. For example, suppose you have a piece of sheet metal with heat applied to certain points around the edge. The heat will gradually spread through the sheet metal until it reaches a point where Laplace's equation is satisfied. The image below shows an example of a surface satisfying Laplace's equation where heat has been applied to the left and right edge.

Source: https://nibot-lab.livejournal.com/48933.html

As you can see in the image, the sheet has been divided into a grid of cells, and each cell has been assigned a temperature. This discretization makes it possible to solve the Laplace equation using the Jacobi method. This is an iterative method where you repeatedly replace the value of a cell with the average of its neighbors.

However, when performing this method in your program, you should not actually make any changes to your matrix because they would affect future calculations. Instead, you should create a new version of your matrix that contains the updated values. This is reflected below by adding a superscript to describe each new version of the matrix.

You will repeat this calculation until you reach a fixed point—that is, a point where \(m^k\) = \(m^{k+1}\).

Implementation

We will represent our sheet metal as an \(M \times N\) array, where \(M\) and \(N\) are selected by the user. We will assume there is a layer of “air” around the sheet metal which does not change. Thus, you should actually represent your sheet metal as an \((M + 2) \times (N + 2)\) array so that there is an extra row and column at the beginning and the end. (This has an added advantage of eliminating edge cases: every cell in the sheet has all four neighbors.)

To avoid confusion, we will always express cells in this sheet as row-column pairs, not X-Y coordinates. Thus, the cell \((3, 7)\) refers to row \(3\), column \(7\). Note that row \(0\) is the top air layer, row \(M + 1\) is the bottom air layer, column \(0\) is the left air layer, and column \(N + 1\) is the right air layer.

We will define our heat sources as rectangular regions in our sheet. We will express rectangular regions as two cells: the upper left corner and the lower right corner. For example, the region \(\{(10, 10), (20, 50)\}\) would be all cells with row \(r\) and column \(c\) such that \(10 \le r \le 20\) and \(10 \le c \le 50\).

To simplify things, we will assume that all heat sources are lines contained in the air around the sheet. If our sheet is an \(M \times N\) grid, then our heat sources will always be rectangles \(\{(r_1, c_1), (r_2, c_2)\}\) such that at least one of the following is true:

• \(r_1 = r_2 = 0\)
• \(r_1 = r_2 = M + 1\)
• \(c_1 = c_2 = 0\)
• \(c_1 = c_2 = N + 1\)

By placing our heat sources in the air region (which does not change) we ensure that the heat sources remain constant.

We will represent temperature as a number between 0 and 1, with 1 being the highest and 0 being the lowest.

Requirements

• The main process should input the following:
  • The number \(M\) of rows in the sheet
  • The number \(N\) of columns in the sheet
  • The number \(h\) of heat sources
  • For each heat source:
    • The row, then column of the upper-left corner
    • The row, then column of the lower-right corner
    • The temperature as a number between 0 and 1, inclusive
• Use the Jacobi method to solve Laplace's equation. To do this, repeatedly update each cell with the average of its neighbors as described above. You should do this until an update results in no significant changes. We will consider a significant change any change \(\ge 0.001\). Thus, you should stop at the first \(m^k\) such that \(\forall_{i, j} | m^{k}_{i, j} - m^{k-1}_{i, j} | \lt 0.001\).
• Output the entire sheet, including the air layer. Output one row per line, separating columns with spaces. You may have a trailing space at the end of your lines.

Parallelization

You must parallelize your solution using MPI. Here is the recommended approach:

• Divide the sheet into even regions by dividing up the rows. You may assume that \(M\) is divisible by the number of processes. Each process will actually need to store \(M + 2\) rows: its \(M\) rows plus the row above and the row beneath.

  For example, suppose \(M = 8\) and we have \(4\) processes. That means we actually have \(M + 2 = 10\) rows counting the air layers. Each process needs \(M / 4 = 2\) rows as follows:

  Process 0 0 1 2 3

  Process 1 2 3 4 5

  Process 2 4 5 6 7

  Process 3 6 7 8 9

  The rows in parentheses will not be edited by that process; they will only be used to update the rows that process is responsible for.

• After each update, the processes will need to communicate the updated rows to their neighbors. For example, process 0 will need to send the updated row 2 to process 1, and process 1 will need to send the updated row 3 to process 0.
• Each process will need to notice if its rows have had significant changes. Then your processes will need to communicate to decide whether any process has had significant changes. (Hint: consider MPI_Allreduce.) If any process had significant changes, then all processes need to continue for another round.

Design

You are welcome to use my matrix class if you wish.

This is a significantly larger MPI program than we have written so far. There are a couple of implications of that fact.

1. You need to have checkpoints. Do not attempt to write the whole thing and then begin debugging. There will be way too many places for bugs to hide. Instead, try to get one small piece working, then verify that it works. For example:
   • Have the main process input and distribute the input. Have each process print out the input to make sure it was received correctly.
   • Have each process print its initial chunk of the sheet to make sure it is correct.
   • Have each process make one iteration of the algorithm and print its chunk to make sure it is correct.
   • Have processes exchange rows and print their chunks to make sure it is correct.
   • etc.

   When printing chunks, I recommend using barriers to make sure output comes out synchronized and not in a jumbled mess. Here is the output function I used:

   void print_all_chunks(int tag) {
   	for (int i = 0; i < num_procs; ++i) {
   		if (my_rank == i) {
   			cout << '[' << tag << "] Process " << my_rank << " has chunk:\n" << my_chunk << endl;
   		}
   		MPI_Barrier(MPI_COMM_WORLD);
   	}
   }
   

   It is recommended that you use debugging macros to make it easy to include/exclude debugging output, like so:

   #ifdef IS_DEBUG
   		if (tag % 1000 == 0)
   			print_all_chunks(tag);
   		++tag;
   #endif
   

   Then if you want to see the debug output you can simply add the following to the top of your program:

   #define IS_DEBUG
   

   Begin by testing your program with small test cases. Here is a good one for 1, 2, or 4 processes:

   4 4
   1
   0 0 0 5 1

2. You need to organize your program to keep it manageable. Your subroutines should be no more than one screenful and should be easy to read and digest. Here's the organization I used for my program.

struct laplace {
	/*----------DATA----------*/

	// Basic parameters
	int num_procs;
	int my_rank;
	int num_rows;
	int num_cols;
	int num_rows_per_proc;

	// The "extended" values include the extra rows and columns
	int ext_rows() { return num_rows + 2; }
	int ext_cols() { return num_cols + 2; }
	int ext_rows_per_proc() { return num_rows_per_proc + 2; }

	// The entire sheet (main only)
	matrix<double> sheet;

	// Each process has their own chunk of the sheet
	matrix<double> my_chunk;

	// Used to calculate next version of chunk without messing up my_chunk
	// When done, you can std::swap(my_chunk, my_scratch)
	matrix<double> my_scratch;


	/*----------FUNCTIONS----------*/

	// Initialize MPI and initialize variables
	void init(int argc, char* argv[]);

	// Main process gets input from user
	void get_input();

	// Main process shares input with others
	void share_input();

	// Outer edges (row 0 and M + 1) only need to be shared once (they don't change)
	void share_outer_edges();

	// Inner edges need to be shared after every update
	void share_inner_edges();

	// Replace each cell with the average of its neigbhors
	int advance();

	// Iterate until a fixed point is reached
	void solve();

	// Each process sends its results back to main process
	void collect_results();

	// Main process outputs final results
	void print_results();

	// Debugging functions
	void print_all_chunks(int tag);
	void log(const std::string& message) {
		std::cout << my_rank << ": " << message << endl;
	}
};