# Diffusion

## The Concept

We are going to look at the MPI communication related to a simplified diffusion problem, namely the Laplace equation in 2D. We will use the five point stencil as our discretisation scheme and a Dirichlet boundary condition. The resulting system is solved by simple Jacobi iteration: start by initialising the boundary and guess a solution for the inner nodes. Then while the system has not converged, update all the inner node to be the average of the 4 neighbouring nodes.

The work can be divided into squares of size NxN that are divided among the processes. The variables lying at the boundary have to be communicated with the neighbouring nodes resulting in a communication cost of order N (four sides of size N). The computation cost per node is N^2. Given that communication is expensive we want the total communication cost to be negligible compared to the total computational cost which will require N to be big.

## The Code

The domain is discreticed into a 6 (width) by 8 (height) grid which is to be divided into 2 by 2 squares handled by different processes; 12 processes in total. Each process also has to store the closest neighbouring data making the matrix storage requirement (width+2)*(height+2)=16. The processes are set up with a communicator that uses a two dimensional Cartesian topology. The key point here is to allow the library implementers to do hardware specific optimizations for such a topology as well as hiding the task of keeping track of neighbours from the programmer. MPI_Cart_coords and MPI_Cart_shift can easily be used to find the coordinates and the rank of neighbouring processes. The program exchange data across process boundaries in four steps by moving data up, down, left and right. The output of the program is information about the ranks of the processes, their coordinates and neighbours as well as the state of the global matrix (which is collected by the master) before the exchange, after the vertical exchange and after the horizontal exchange.

##### Previous instructions

MPI_Init() and MPI_Finalize(); Are used to initialize and end the MPI program.

MPI_Comm_rank(); Is used to get the rank as usual, only this time using both MPI_COMM_WORLD and the cartesian_communicator.

MPI_Comm_size(); Is used to get the total number of processes.

MPI_Send() and MPI_Recv(); Are used to pass data from the slaves to the master

##### New instructions

MPI_Comm cartesian_communicator; Declares the Cartesian communicator.

MPI_Cart_create( MPI_COMM_WORLD, PROCESS_DIMENSIONS, divisions, periods, reorder, &cartesian_communicator ); Creates a Cartesian topology from MPI_COMM_WORLD. where PROCESS_DIMENSIONS (2) is the dimension, divisions is an array specifying the size of each dimension.  periods is a Boolean specifying for each dimension if it is periodic. reorder is a Boolean value that sets if it is OK to reorder the processes (e.g. to improve performance based their physical location.) cartesian_communicator is the handle to the new communicator.

MPI_Cart_coords( cartesian_communicator, my_cartesian_rank, PROCESS_DIMENSIONS, my_coordinates ); Places the Cartesian coordinates of my_cartesian_rank within the cartesian_communicator in my_coordinates, which is has the same size as the number of dimensions specified by PROCESS_DIMENSIONS.

MPI_Cart_shift( cartesian_communicator, SIDEWAYS, RIGHT, &left_neighbour, &right_neighbour ); Is used to find the left_neighbour and right_neighbour within the cartesian_communicator. A slightly different call using UPDOWN and UP finds the bottom_neighbour and top_neighbour.

MPI_Sendrecv( &matrix[ROWS - 2][0], COLUMNS, MPI_INT,  top_neighbour, 4004, &matrix[0][0], COLUMNS, MPI_INT, bottom_neighbour, 4004, cartesian_communicator, &status ); Is used to pass the boundary of each process to the process above; each process sends to the top_neighbour and receives from the bottom_neighbour. The data being sent is the first row the top_neighbour does not have which is the second row from the top of the sending process (&matrix[ROWS - 2][0]). The data is placed in the first row of the receiving process (&matrix[0][0]). COLUMNS number of MPI_INT is both sent and received using the tag 4004. Similar calls are then used to transfer data from top to bottom, right to left and, finally, left to right.

MPI_Datatype column_type; Declares the column_type.

MPI_Type_vector(ROWS, 1, COLUMNS, MPI_INT, &column_type ); Is used to create a datatype for the columns to be used with MPI_Sendrecv(). The arrays in C are stored in row-major order, meaning that the stride between the each element in a row is 1. When sending a row with MPI_Sendrecv() it is therefore sufficient to use MPI_INT as type, but when sending columns we need to take into account the COLUMNS stride between each element in a column. This is done by creating the column_type, that contains ROWS elements with COLUMNS stride. The old datatype was MPI_INT and there is only 1 element in each block.

MPI_Type_commit(&column_type); Is used to commit the column_type. This is necessary in order to use it in a communication.

## Compile & Run

If you have not already done so, obtain all the example code here.

Switch to the Intel compiler (optional, only necessary once in each terminal session)

Compile the program using

Submit the job to the queue

The output from the program execution is placed in the output folder

• No labels