Point-to-point communications
The equation governing heat diffusion (or thermal conduction) is known as the heat equation:
where is a positive coefficient called the thermal diffusivity, which depends on the material properties.
In one spatial dimension, the heat equation is given by:
Discretizing the Heat Equation
For numerical simulations, the heat equation is often discretized using finite differences. For example, using the explicit finite difference method, the time derivative can be approximated as:
and the second spatial derivative as:
Substituting these approximations into the heat equation gives a scheme to calculate from :
How to make it parallel with MPI
In parallelizing the heat diffusion equation, the computational domain is divided into subdomains, each handled by a separate process. This division allows each process to handle a smaller subset of the domain, reducing the computational load per process and enabling larger simulations to run more efficiently.
Let’s assume a one-dimensional domain with points divided into subdomains, where is the number of processes. Each process is responsible for updating the temperature values for a certain range of -values within the domain. For example:
- Process will handle points from to .
- Process will handle points from to .
- And so on, until the last process handles to .
For example, if we have points and processes, each process will handle points. We will divide the domain as follow:
- Process will handle points in range .
- Process will handle points in range .
- Process will handle points in range .
- Process will handle points in range .
Boundary Communication
In each time step, updating the temperature at a given point requires information from neighboring points. Specifically, points near the boundaries of each process’s subdomain need the latest temperature values from the adjacent processes. This is where boundary communication becomes necessary:
- Exchange of Boundary Values: Each process needs to send its boundary data (the temperature at its boundary points) to its neighboring processes. For instance:
- Process must send the temperature at its rightmost point to Process and receive the leftmost boundary value from Process .
- Similarly, it sends its leftmost boundary to Process i-1 and receives the rightmost boundary value from Process .
- Using
MPI_SendandMPI_Recv: Communication can be implemented withMPI_SendandMPI_Recv:
- Non-blocking communication (e.g., using
MPI_IsendandMPI_Irecv) can allow computation to overlap with communication, potentially improving performance. - Boundary update ensures that each process has the necessary information before updating the temperature values at boundary points.
C code
Here the MPI C code:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <mpi.h>
// Initialize temperature
void initialize(double *temperature, int local_n, int rank, int size) {
// Initialize the temperature array (e.g., set an initial heat distribution)
for (int i = 0; i < local_n; i++) {
temperature[i] = 0.;
}
// Boundary conditions: high temperature at left edge
if (rank == 0) temperature[0] = 1.0;
}
// Update temperature using the heat equation
void update_temperature(double *temp, double *new_temp, int local_n, double alpha, double dx, double dt, int rank, int size) {
for (int i = 1; i < local_n-1; i++) {
new_temp[i] = temp[i] + alpha*(dt / (dx * dx)) * (temp[i + 1] - 2 * temp[i] + temp[i - 1]);
}
}
// Write data to a file for Gnuplot visualization
void write_results(int rank, int size, int local_n, double *temp, int step, int NX) {
char filename[256];
sprintf(filename, "heat_output_step_%d.dat", step);
FILE *fp = fopen(filename, rank == 0 ? "w" : "a");
//Write header
if (rank == 0) fprintf(fp, "# Step %d\n", step);
//Buffers used to send and recv data during REDUCE (used to print data)
double *sendbuff = (double *)malloc(NX * sizeof(double));
double *recvbuff = (double *)malloc(NX * sizeof(double));
for (int i = 0; i < NX; i++){
sendbuff[i]=0;
recvbuff[i]=0;
}
for (int i = 1; i < local_n-1; i++) sendbuff[(local_n-2)*rank+(i-1)] = temp[i]; //We remove from local_n the boundaries (that are 2: left and right)
//All processors reduce on root their local values
MPI_Reduce(sendbuff, recvbuff, NX, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
//Only root prints
if(rank==0){
for (int i = 0; i < NX; i++) fprintf(fp, "%f\n", recvbuff[i]);
}
fclose(fp);
free(sendbuff);
free(recvbuff);
}
int main(int argc, char *argv[]) {
double dx, dt = 0.001, alpha = 0.01;
int NX = 100; //Number of grid points
int NT = 1e5; //Number of time steps
int time_print = 1e3; //Print frequency
int ln, rn; //Left and Right Neighbours ranks
int rank, size;
int local_n;
double *temp, *new_temp;
MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &size);
ln = rank-1;
rn = rank+1;
if(rank==0) ln = rank;
if(rank==size-1) rn = rank;
dx = 1./((double)NX);
local_n = NX / size + 2; //We add 2 beacuse, for each process, we add the "boundaries" at 0 and Nx/size+1
temp = (double *)malloc(local_n * sizeof(double));
new_temp = (double *)malloc(local_n * sizeof(double));
// Initialize the temperature field
initialize(temp, local_n, rank, size);
// Time-stepping loop
for (int t = 0; t < NT; t++) {
// Exchange boundary points
if(rank!=0) MPI_Send(&temp[1], 1, MPI_DOUBLE, ln, 0, MPI_COMM_WORLD);
if(rank!=size-1) MPI_Send(&temp[local_n - 2], 1, MPI_DOUBLE, rn, 1, MPI_COMM_WORLD);
if(rank!=size-1) MPI_Recv(&temp[local_n - 1], 1, MPI_DOUBLE, rn, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
if(rank!=0) MPI_Recv(&temp[0], 1, MPI_DOUBLE, ln, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
// Update temperature field
update_temperature(temp, new_temp, local_n, alpha, dx, dt, rank, size);
// Swap pointers
double *swap = temp;
temp = new_temp;
new_temp = swap;
// Write results for Gnuplot
//if(t%time_print==0) write_results(rank, size, local_n, temp, t, NX);
}
free(temp);
free(new_temp);
MPI_Finalize();
return 0;
}
MPI_Send
The MPI_Send function is a core part of MPI that sends a message from one process to another within an MPI communicator. It performs a standard-mode blocking send (below there is a section explaining blocking vs. non-blocking operations), meaning it only returns once the send buffer is safe to modify, which may occur when the message is delivered. Here’s a breakdown of its syntax in both C and Fortran.
The signature in C is:
int MPI_Send(const void *buf, int count, MPI_Datatype datatype, int dest, int tag, MPI_Comm comm);
buf: Pointer to the data being sent.count: Number of elements in the send buffer.datatype: Data type of each element (e.g.,MPI_INT,MPI_DOUBLE).dest: Rank of the destination process within the communicator.tag: Message tag to identify the message type, used in pairing sends and receives.comm: Communicator (usuallyMPI_COMM_WORLD).
The signature in Fortran is:
MPI_Send(buf, count, datatype, dest, tag, comm, ierror) <type>, dimension(*) :: buf integer :: count, datatype, dest, tag, comm, ierror
buf: Array containing data to be sent.count: Number of elements inbuf.datatype: MPI data type of each element (e.g.,MPI_INTEGER,MPI_DOUBLE_PRECISION).dest: Destination process rank.tag: Message tag identifier.comm: Communicator.ierror: Error code.
MPI_Recv
The MPI_Recv function is the counterpart to MPI_Send, and it allows a process to receive a message from another process in a blocking fashion. This means that the receiving process will wait (block) until the requested message is received, and it can then proceed with execution. The message is stored in a buffer that is provided by the receiving process.
The matching between the send and receive operations is done using the message's tag and destination rank. The tag allows the program to distinguish between different types of messages.
The signature in C is:
int MPI_Recv(void *buf, int count, MPI_Datatype datatype, int source, int tag, MPI_Comm comm, MPI_Status *status);
buf: Pointer to the buffer where the received data will be stored.count: The maximum number of elements to receive (this should match the number of elements sent by the sender).datatype: Data type of each element (e.g.,MPI_INT,MPI_DOUBLE).source: The rank of the sending process, orMPI_ANY_SOURCEif the message can come from any process.tag: The message tag used to identify the message, orMPI_ANY_TAGto accept messages with any tag.comm: Communicator (usuallyMPI_COMM_WORLD).status: A pointer to anMPI_Statusstructure, which is used to store information about the received message, such as the actualsourceandtagof the message.
The signature in Fortran is:
MPI_Recv(buf, count, datatype, source, tag, comm, status, ierror) <type>, dimension(*) :: buf integer :: count, datatype, source, tag, comm, ierror integer, dimension(MPI_STATUS_SIZE) :: status
buf: The array where the received data will be stored.count: The number of elements to be received.datatype: MPI data type of each element (e.g.,MPI_INTEGER,MPI_DOUBLE_PRECISION).source: The rank of the sending process, orMPI_ANY_SOURCEif the message can come from any process.tag: The message tag used to identify the message, orMPI_ANY_TAGto accept messages with any tag.comm: Communicator (usuallyMPI_COMM_WORLD).status: An array to hold the message status information (e.g.,source,tag).ierror: Error code.
Blocking vs. Non-Blocking calls
In MPI, a blocking call refers to a function that causes the calling process to wait until the operation completes before it continues execution. This means the process will block, or be paused, until it receives a response or the data it is waiting for has been transmitted or received.
For example, when a process calls MPI_Send to send a message, it will block until the message is sent (or at least successfully enqueued for transmission, depending on the specific call).
Similarly, MPI_Recv is blocking: the receiver will block until the message is received, meaning it waits for the message to arrive.
Blocking calls can simplify program flow because the program doesn't need to handle asynchronous communication. However, blocking calls can also lead to inefficiencies, especially in parallel applications where multiple processes might be waiting on each other.
In contrast, non-blocking calls like MPI_Isend and MPI_Irecv allow processes to continue execution while the communication is being handled in the background, and completion must be explicitly checked (using MPI_Wait, for example).
Blocking calls are generally easier to program but can lead to performance bottlenecks, especially when processes have to wait for each other to send or receive messages.
Example
Here is an example in which the code gets blocked due to incorrect use of MPI_Send and MPI_Recv. This happens when the sender and receiver are not properly matched, or if there is a mismatch between the expected number of elements or tags.
Suppose we have two processes. One process is sending data, and the other process is supposed to receive it. The issue occurs when the receiver is not correctly set up to receive the message, such as by using the wrong tag or not providing enough space in the receive buffer.
#include <stdio.h>
#include <mpi.h>
int main(int argc, char *argv[]) {
int rank, size;
int send_data = 123; // Data to send
int recv_data; // Buffer to receive data
MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &size);
if (size < 2) {
printf("This example requires at least 2 processes.\n");
MPI_Finalize();
return 1;
}
if (rank == 0) {
// Process 0 sends data with tag 100
printf("Process %d sending data: %d\n", rank, send_data);
MPI_Send(&send_data, 1, MPI_INT, 1, 100, MPI_COMM_WORLD);
}
else if (rank == 1) {
// Process 1 expects data with tag 200 (incorrect tag)
printf("Process %d waiting for data...\n", rank);
MPI_Recv(&recv_data, 1, MPI_INT, 0, 200, MPI_COMM_WORLD, MPI_STATUS_IGNORE); // Incorrect tag here
printf("Process %d received data: %d\n", rank, recv_data);
}
MPI_Finalize();
return 0;
}
To fix it, just change the tag to match the send and receive message.
Visualisation
You might use the following Python script to create an animation showing the time evolution of the temperature field:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
# Parameters
NX = 100 # Number of spatial points
NT = 99900 # Total number of time steps
every = 1000 # Time step interval for displaying the frame
# Initialize x positions and dummy y positions for a "bar" look
x = np.linspace(0, NX - 1, NX)
y = np.ones(NX) # Fixed y-position for visualization
# Initialize the figure
fig, ax = plt.subplots(figsize=(10, 2))
colorbar = plt.colorbar(plt.cm.ScalarMappable(cmap="coolwarm"), ax=ax, orientation='horizontal', label='Temperature')
plt.xlim([0, NX])
plt.ylim([0.995, 1.005])
ax.set_title("1D Heat Diffusion")
# Define update function to load data from file and update plot
def update(frame):
# Construct the filename for the current time step
file = f'heat_output_step_{frame}.dat'
# Load temperature data
try:
T = np.loadtxt(file, skiprows=1) # Load temperature data from the file
# Clear previous scatter and reinitialize it with the new color data
ax.collections.clear() # Clear old scatter plot objects
scat = ax.scatter(x, y, c=T, cmap="coolwarm", marker='s', s=100) # Reinitialize color data
ax.set_title(f"1D Heat Diffusion - Step {frame}")
except Exception as e:
print(f"Error loading or processing file {file}: {e}")
return scat,
# Create the animation
frames = range(0, NT, every) # Use only every nth frame
ani = FuncAnimation(fig, update, frames=frames, blit=True)
plt.close(fig)
# Display the animation in the notebook
HTML(ani.to_jshtml())