8. Domain Decomposition with MPI
In Fortran numerical codes, a widely used approach is to divide the computational domain into pieces, let each process compute its own piece, and pass only the boundary values with MPI. MISO takes the same approach to parallel computing. In this chapter, we first run the basic MPI calls in short test programs, and then read how MISO divides the domain, how it exchanges values with its neighbors, and how it synchronizes the time step.
As we saw in Chapter 7, the Host version of MISO’s for_each is just a plain for loop and does not use OpenMP. The only way to run MISO in parallel on a CPU is the MPI covered in this chapter. Everything in this chapter can be tried on a CPU alone. Communication when using a GPU is covered in Chapter 9.
Contents
1. Processes and Ranks
Running the same program in multiple processes
mpiexec -np 4 ./program starts the same program ./program simultaneously as four processes. Each process has its own variables and arrays. Even if one process modifies an array, the arrays of the other processes do not change. When a process needs the values of another process, it must explicitly send and receive them with MPI functions. This distinction holds even when the processes run on the same computer.
The group of processes that take part in communication is called a communicator. MPI_COMM_WORLD is the group of all processes started by mpiexec. The number, starting from 0, that each process receives within a communicator is its rank. With four processes, the ranks are 0, 1, 2, and 3.
Every process runs the same program, but each one checks its own rank and changes the region or the work it handles. This pattern is called SPMD (Single Program, Multiple Data).
Test program: printing the rank
The following program has each process print its own rank and the total number of processes. Save it as hello_mpi.cpp.
#include <iostream>
#include <mpi.h>
int main(int argc, char **argv) {
MPI_Init(&argc, &argv); // start using MPI (called by all ranks)
int rank = 0; // own rank number
int size = 0; // total number of ranks
MPI_Comm_rank(MPI_COMM_WORLD, &rank); // written to rank
MPI_Comm_size(MPI_COMM_WORLD, &size); // written to size
std::cout << "rank " << rank << " / size " << size << std::endl;
MPI_Finalize(); // finish using MPI (called by all ranks)
return 0;
}
#include <mpi.h> brings in the declarations of the MPI functions. MPI functions are called after MPI_Init and before MPI_Finalize. &rank is the address of a variable, which we learned about in Chapter 3. MPI_Comm_rank writes the rank number to that address.
From C++, we call MPI’s C-language functions. The differences from Fortran are as follows.
Fortran |
C++ (C functions) |
|
|---|---|---|
Initialization |
|
|
Rank |
|
|
Error code |
Returned in the last argument |
Returned as the function’s return value ( |
Type specification |
|
|
Fortran passes arguments by reference, so you pass rank as it is. For a C function to write into it, however, you pass the address &rank. You do not have to receive the returned error code if you do not use it.
Compile with the mpicxx of the OpenMPI we installed in Chapter 1, and launch with mpiexec. mpicxx is a command that calls g++ after adding the locations of the MPI headers and libraries.
mpicxx -std=c++17 -o hello_mpi hello_mpi.cpp
mpiexec -np 2 ./hello_mpi
mpiexec -np 4 ./hello_mpi
This is an example of the output when run with four processes.
rank 0 / size 4
rank 2 / size 4
rank 1 / size 4
rank 3 / size 4
The order of the lines may change from run to run, because each process writes to the screen independently. If you give -np a number larger than the number of cores on the computer, OpenMPI refuses to start by default. When you practice on a computer with few cores, start the run with the environment variable OMPI_MCA_rmaps_base_oversubscribe=yes (Section 10).
2. MISO’s Env: Starting and Ending MPI
The main of a MISO demo first creates an Env (demo/mhd2d_vortex/src/main.cpp).
int main(int argc, char **argv) {
// Initialize MPI and CUDA environments
Env env(argc, argv);
// Read configuration file
auto config_path = parse_config_filepath(argc, argv);
Config config(config_path.value_or("./config.yaml"));
// Run simulation
Model model(config);
model.run();
}
Neither MPI_Init nor MPI_Finalize appears here. Env takes care of both (miso/include/miso/env.hpp). The MPI part is miso::mpi::Env. Here is an excerpt. In the excerpts in this chapter, comments of the form // ...(description) were added by this course to mark omitted parts.
struct Env {
bool owns_mpi = false;
void setup(int *argc, char ***argv) {
int is_mpi_initialized = false;
MPI_Initialized(&is_mpi_initialized);
if (!is_mpi_initialized) {
owns_mpi = true;
MPI_Init(argc, argv);
}
world::comm = MPI_COMM_WORLD;
MPI_Comm_rank(world::comm, &world::rank);
MPI_Comm_size(world::comm, &world::size);
// ...(finds the ranks within the node; used in Chapter 9)
world::is_initialized = true;
}
Env(int &argc, char **&argv) { setup(&argc, &argv); }
// ...(there is also a constructor with no arguments)
~Env() noexcept {
// ...(frees the intra-node communicator)
int is_mpi_finalized = false;
MPI_Finalized(&is_mpi_finalized);
if (owns_mpi && !is_mpi_finalized) {
MPI_Finalize();
}
}
};
The constructor calls MPI_Init and the destructor calls MPI_Finalize. This is the RAII pattern we learned in Chapter 4. When the lifetime of env ends at the end of main, the destructor runs, so you cannot forget to write MPI_Finalize. MPI_Initialized checks whether MPI has already been initialized, and the finalization is done only when Env itself initialized MPI (owns_mpi). This is so that Env also works when MPI has been started elsewhere, as in tests.
The rank number is stored in a variable declared in the namespace mpi::world of the same file, like inline int rank = -1;. A variable placed outside functions lasts from the start to the end of the program (similar to a module variable in Fortran). inline ensures that the variable becomes a single one even when the file is included from several files (the same reason as for the inline functions in Chapter 2).
Elsewhere in MISO, these variables are not read directly. Instead, the following functions are used.
Function |
Value returned |
|---|---|
|
Own rank in |
|
Total number of ranks |
|
|
|
The global communicator ( |
Calling any of these functions before creating the Env throws the exception std::runtime_error. This reports at run time that the convention “create Env at the start of main” was not followed.
The miso::Env used in main bundles mpi::Env and, only in CUDA builds, cuda::Env (selection of the GPU to use). In a CPU build, only mpi::Env is active. cuda::Env is covered in Chapter 9.
3. The Arrangement of Ranks and Neighboring Ranks: mpi::Shape
The mpi section of config.yaml
In the configuration file, MISO specifies how many pieces the three-dimensional grid is divided into in each of the x, y, and z directions. This is the relevant part of demo/mhd2d_vortex/config.yaml.
grid:
i_size: 128
j_size: 128
k_size: 1
margin: 2
# ...
mpi:
x_procs: 2
y_procs: 1
z_procs: 1
domain:
periodic:
x: true
y: true
z: true
x_procs: 2 means dividing into two in the x direction. x_procs * y_procs * z_procs must equal the number of processes given to mpiexec -np. The mpi section has no default values (there is no mpi entry in miso/include/miso/config_defaults.hpp), so it is written in the configuration file of every demo. Among the other demos, demo/hd2d_rayleigh_taylor/config.yaml divides into two in the y direction (y_procs: 2), and the one-dimensional shock tube demos (such as demo/hd1d_shock_tube/config/config_x.yaml) are not divided (all values are 1).
MPI_Cart_create and MPI_Cart_shift
mpi::Shape reads this setting and builds the arrangement of ranks (miso/include/miso/mpi_util.hpp). Here is an excerpt of the main part of the constructor.
x_procs = config["mpi"]["x_procs"].as<int>();
y_procs = config["mpi"]["y_procs"].as<int>();
z_procs = config["mpi"]["z_procs"].as<int>();
int dims[ndims] = {x_procs, y_procs, z_procs};
int periods[ndims];
periods[0] = config["domain"]["periodic"]["x"].as<bool>() ? 1 : 0;
periods[1] = config["domain"]["periodic"]["y"].as<bool>() ? 1 : 0;
periods[2] = config["domain"]["periodic"]["z"].as<bool>() ? 1 : 0;
MPI_Cart_create(MPI_COMM_WORLD, ndims, dims, periods, 0, &cart_comm);
MPI_Comm_rank(cart_comm, &myrank);
MPI_Comm_size(cart_comm, &n_procs);
// ...(check the number of processes)
MPI_Cart_coords(cart_comm, myrank, ndims, coord);
MISO_MPI_CHECK(MPI_Cart_shift(cart_comm, 0, 1, &x_procs_neg, &x_procs_pos));
MISO_MPI_CHECK(MPI_Cart_shift(cart_comm, 1, 1, &y_procs_neg, &y_procs_pos));
MISO_MPI_CHECK(MPI_Cart_shift(cart_comm, 2, 1, &z_procs_neg, &z_procs_pos));
The roles of the MPI functions used are as follows.
MPI_Cart_create: creates a new communicatorcart_commin which all processes are arranged in a grid of the shapedims(here \(2 \times 1 \times 1\)). In the directions given inperiods, the two ends are connected (periodic boundary). The fifth argument0means “do not reorder the rank numbers”.MPI_Cart_coords: finds where your own rank is in the grid (coord[0],coord[1],coord[2]).MPI_Cart_shift: finds the neighboring rank shifted by one in the specified direction (0 is x). The one on the negative side goes intox_procs_neg, and the one on the positive side goes intox_procs_pos.
At the ends in a non-periodic direction there is no neighbor. In that case, MPI_PROC_NULL is stored as the neighboring rank. MPI_PROC_NULL is a special value meaning “no partner”, and sending to or receiving from it does nothing and finishes immediately. MISO uses this value to decide on the halo exchange and the physical boundary, which we will see in later sections.
MISO_MPI_CHECK is a macro that, if the return value (error code) of an MPI function is not success, prints a description of the error and stops all processes with MPI_Abort (check_error in the same file).
The “process count check” in the constructor prints a description and stops when x_procs * y_procs * z_procs differs from the actual number of processes. However, when we checked with OpenMPI 4.1, if the numbers differ, the earlier MPI_Cart_create (or the MPI_Comm_rank right after it) stops with an MPI error, and this description is not shown. If you get an error from MPI_Cart_create or MPI_Comm_rank, first check whether the mpi section matches the number given to -np.
The order of the rank numbers
MPI_Cart_create orders the rank numbers so that the last dimension (z) changes fastest. When we divide into two in the x direction and two in the y direction, the correspondence between ranks and positions is as follows.
y
^
coord_y=1 | rank 1 rank 3
coord_y=0 | rank 0 rank 2
+--------------------> x
coord_x=0 coord_x=1
This is the same rule as the layout of C++ arrays that we learned in Chapter 7 (the last index is contiguous). Note that it is the opposite of the layout of Fortran arrays. At the start of a run, MISO writes this correspondence to mpi/coords.csv under io.save_dir (Shape::save). When divided as \(2 \times 2\), the contents are as follows.
rank,x,y,z
0,0,0,0
1,0,1,0
2,1,0,0
3,1,1,0
pymiso reads this file and joins the per-rank outputs to build the global array.
4. The Global Grid and the Per-Rank Grid
Assigned cells and ghost cells
Let us think in one dimension. If we divide 16 cells in total between two ranks, each rank is responsible for 8 cells.
global index 0 1 2 3 4 5 6 7 | 8 9 10 11 12 13 14 15
owner rank 0 | rank 1
MISO’s spatial derivative is a fourth-order central difference, which uses the values from i-2 to i+2 for the derivative at cell i. At the cells at the edge of the assigned range, two values held by the neighboring rank are needed. So we prepare margin extra cells at each end of the array and put the values received from the neighbor there. These are called ghost cells or margins (“nori-shiro” in Japanese). The process of receiving values from the neighbor and filling them in is called halo exchange.
When margin: 2, the array of rank 1 in the x direction looks as follows.
rank 1 local index 0 1 | 2 3 4 5 6 7 8 9 | 10 11
role ghost | 8 owned cells | ghost
global index 6 7 | 8 9 10 11 12 13 14 15 | 16 17
The global indices 16 and 17 of the ghost cells at the right end correspond to 0 and 1 of rank 0 in the case of a periodic boundary. The ghost cells at a non-periodic end (a physical boundary) are not received from a neighbor but are filled by the boundary condition (Section 8).
The values of Grid on each rank
MISO’s Grid<Real, backend::Host> (miso/include/miso/grid.hpp) holds the grid of its own rank, not the global grid. This is an excerpt of the constructor.
const int margin = config["grid"]["margin"].as<int>();
const int i_size_g = config["grid"]["i_size"].as<int>();
// ...(same for j_size_g, k_size_g)
i_size = i_size_g / mpi_shape.x_procs;
j_size = j_size_g / mpi_shape.y_procs;
k_size = k_size_g / mpi_shape.z_procs;
// ...
const int i_offset = mpi_shape.coord[0] * i_size;
// ...
auto x_grid = AxisGrid<Real>(i_size, margin, i_offset, x_provider);
i_size in the configuration file is the total number of cells, whereas the member i_size of Grid is the number of cells that one rank is responsible for. They have the same name but different meanings, so the code distinguishes the global count by adding _g (for global). i_offset is the position, counted over the whole domain, of the first cell assigned to this rank.
The constructor of AxisGrid, the grid for one direction, determines the following values.
size = size_;
stride = size > 1 ? 1 : 0;
margin = margin_ * stride;
total = size + 2 * margin;
begin = margin;
end = margin + size;
offset = offset_;
// ...
for (int i = 0; i < total; ++i) {
const int ig = i - begin + offset;
s[i] = ap.loc(ig);
// ...
}
Grid holds these as i_total (the number of cells including ghost cells) and i_margin (the number of ghost cells). The local indices of the assigned cells are from i_margin up to but not including i_total - i_margin. This is why the loops in Integrator::cfl and others in Chapter 7 started from grid.i_margin.
In a direction where size is 1 (such as the z direction in a two-dimensional calculation), stride becomes 0 and margin also becomes 0. No ghost cells are created in that direction, and it is not divided.
The coordinate s[i] is computed after converting the local index i to the global index ig. The ig of a ghost cell can fall outside the assigned range (a negative value, or a value not less than the total number of cells), but the coordinate is obtained directly from the equation for a uniform grid.
If the total number of cells is not divisible by the number of divisions, check_divisibility throws an exception and the program stops. MISO handles only divisions in which all ranks are responsible for the same number of cells.
Test program: MISO’s Shape and Grid
Using MISO’s headers, we print the position of each rank, its neighboring ranks, and the grid values. We use the CMake configuration from Chapter 1 for “using MISO from your own program” (linking MISO::MISO with add_subdirectory) and build the following main.cpp.
#include <iostream>
#include <miso/core.hpp>
using namespace miso;
using Real = float;
int main(int argc, char **argv) {
Env env(argc, argv); // MPI_Init (MPI_Finalize in the destructor)
Config config("shape_check.yaml"); // read the configuration file
mpi::Shape mpi_shape(config); // rank layout and neighboring ranks
Grid<Real, backend::Host> grid(config, mpi_shape); // grid of this rank
namespace bc = miso::boundary_condition;
const bool x_inner = bc::is_physical_boundary(Direction::X, Side::INNER, mpi_shape);
const bool y_inner = bc::is_physical_boundary(Direction::Y, Side::INNER, mpi_shape);
for (int r = 0; r < mpi::size(); ++r) { // one rank at a time so the output does not interleave
if (r == mpi::rank()) {
std::cout << "rank " << mpi::rank() << ": coord = (" << mpi_shape.coord[0]
<< ", " << mpi_shape.coord[1] << ")"
<< ", x_neg/x_pos = " << mpi_shape.x_procs_neg << "/"
<< mpi_shape.x_procs_pos
<< ", y_neg/y_pos = " << mpi_shape.y_procs_neg << "/"
<< mpi_shape.y_procs_pos << "\n"
<< " i_size = " << grid.i_size
<< ", i_total = " << grid.i_total << ", x[i_margin] = "
<< grid.x[grid.i_margin] << ", y[j_margin] = "
<< grid.y[grid.j_margin] << "\n"
<< " physical boundary: x inner = " << x_inner
<< ", y inner = " << y_inner << std::endl;
}
MPI_Barrier(mpi::comm());
}
}
MPI_Barrier is a function that waits until all ranks reach that point. We use it here to make the ranks print in order starting from rank 0. is_physical_boundary is explained in Section 8.
Put the configuration file shape_check.yaml in the folder where you run the program. It is set to be periodic in the x direction and non-periodic in the y direction. io.enabled: false disables file output.
io:
enabled: false
grid:
i_size: 8
j_size: 8
k_size: 1
margin: 2
x_min: 0.0
x_max: 1.0
y_min: 0.0
y_max: 1.0
z_min: 0.0
z_max: 1.0
mpi:
x_procs: 2
y_procs: 2
z_procs: 1
domain:
periodic:
x: true
y: false
z: true
This is the output when run with four processes (assuming the executable was created under build).
mpiexec -np 4 ./build/shape_check
rank 0: coord = (0, 0), x_neg/x_pos = 2/2, y_neg/y_pos = -2/1
i_size = 4, i_total = 8, x[i_margin] = 0.0625, y[j_margin] = 0.0625
physical boundary: x inner = 0, y inner = 1
rank 1: coord = (0, 1), x_neg/x_pos = 3/3, y_neg/y_pos = 0/-2
i_size = 4, i_total = 8, x[i_margin] = 0.0625, y[j_margin] = 0.5625
physical boundary: x inner = 0, y inner = 0
rank 2: coord = (1, 0), x_neg/x_pos = 0/0, y_neg/y_pos = -2/3
i_size = 4, i_total = 8, x[i_margin] = 0.5625, y[j_margin] = 0.0625
physical boundary: x inner = 0, y inner = 1
rank 3: coord = (1, 1), x_neg/x_pos = 1/1, y_neg/y_pos = 2/-2
i_size = 4, i_total = 8, x[i_margin] = 0.5625, y[j_margin] = 0.5625
physical boundary: x inner = 0, y inner = 0
Let us summarize what we can read from it.
The positions of the ranks are as in the figure in Section 3 (rank 1 is at
(0, 1)).The x direction is periodic and divided into two, so the left and right neighbors are the same rank (both the left and right of rank 0 are rank 2).
The y direction is not periodic, so
y_negof ranks 0 and 2 at the bottom end andy_posof ranks 1 and 3 at the top end are-2.-2is the value ofMPI_PROC_NULLin this OpenMPI (the value differs between MPI implementations, so in code you must always compare withMPI_PROC_NULL).The 8 cells in total are divided into two, so
i_size = 4, and adding two ghost cells on each side givesi_total = 8.The coordinates of the center of the first assigned cell are \(0.5/8 = 0.0625\) and \(4.5/8 = 0.5625\), depending on the position of the rank.
If you run with mpiexec -np 2, it does not match x_procs * y_procs * z_procs = 4, so it stops with an error from MPI_Cart_create.
Exercise 1: Per-rank grid values
With demo/mhd2d_vortex/config.yaml (i_size: 128, margin: 2, x_min: 0.0, x_max: 1.0, x_procs: 2), find grid.i_size, grid.i_total, and i_offset for the rank with coord[0] = 1. Then, for the cells with local indices i = 2 and i = 0 on that rank, find the global index ig and the coordinate grid.x[i].
i_size = 128 / 2 = 64, i_total = 64 + 2 * 2 = 68, and i_offset = 1 * 64 = 64.
In AxisGrid, ig = i - begin + offset, and begin = margin = 2. For i = 2, ig = 2 - 2 + 64 = 64, which is the first assigned cell. Its coordinate is \((64 + 0.5) \times (1/128) = 0.50390625\). For i = 0, ig = 62, which is the left ghost cell. It corresponds to the second-to-last assigned cell of rank 0, and its coordinate is \((62 + 0.5)/128 = 0.48828125\).
5. Exchanging Values with Neighbors: MPI_Sendrecv
Matching sends and receives
In one-to-one communication, the sending side and the receiving side each specify the following.
Argument |
Meaning |
|---|---|
Location |
The address of the start of the values to send or of the place to receive into (such as |
Count |
The number of elements |
Type |
|
Partner |
The rank of the partner |
Tag |
An integer that distinguishes communications |
Communicator |
The group of processes within which the communication takes place |
A send and a receive are matched when the “partner” and “tag” specified by the receiving side agree with the rank and tag of the sending side. When you receive two kinds of values from the same partner, you can tell them apart by the tag.
Test program: exchange with neighbors under a periodic boundary
We build a case in which each rank is responsible for 4 cells in one dimension and has one ghost cell at each end. We put each cell’s global index into it as its value, and exchange with the left and right neighbors. Save it as sendrecv.cpp.
#include <iostream>
#include <vector>
#include <mpi.h>
int main(int argc, char **argv) {
MPI_Init(&argc, &argv);
int rank = 0, size = 0;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &size);
const int n = 4; // number of cells owned by each rank
std::vector<double> u(n + 2, -1.0); // u[0] and u[n+1] are ghost cells
for (int i = 1; i <= n; ++i) {
u[i] = rank * n + (i - 1); // value = global index
}
// left and right neighbors under periodic boundaries
const int left = (rank - 1 + size) % size;
const int right = (rank + 1) % size;
// (1) send u[n] to the right, receive into u[0] from the left
MPI_Sendrecv(&u[n], 1, MPI_DOUBLE, right, 100,
&u[0], 1, MPI_DOUBLE, left, 100,
MPI_COMM_WORLD, MPI_STATUS_IGNORE);
// (2) send u[1] to the left, receive into u[n+1] from the right
MPI_Sendrecv(&u[1], 1, MPI_DOUBLE, left, 200,
&u[n + 1], 1, MPI_DOUBLE, right, 200,
MPI_COMM_WORLD, MPI_STATUS_IGNORE);
std::cout << "rank " << rank << ": ghost left = " << u[0]
<< ", inner = " << u[1] << ".." << u[n]
<< ", ghost right = " << u[n + 1] << std::endl;
MPI_Finalize();
return 0;
}
MPI_Sendrecv does “send” and “receive” at the same time in a single call. The first five arguments are for the send, the next five are for the receive, and the last are the communicator and the place to store detailed information about the receive. MPI_STATUS_IGNORE means “detailed information is not needed”.
% size is the remainder of a division; it connects the left of rank 0 to the last rank and the right of the last rank to rank 0 (periodic boundary). We add size in (rank - 1 + size) because in C++ the remainder of a negative number is negative (-1 % 4 is -1).
In (1), all ranks simultaneously “send to the right and receive from the left”. The value that rank 1 sends to the right is received by rank 2 as “receive from the left”. They match because both use tag 100. (2) is the exchange in the opposite direction.
mpicxx -std=c++17 -o sendrecv sendrecv.cpp
mpiexec -np 2 ./sendrecv
mpiexec -np 4 ./sendrecv
This is the output with four processes (sorted in rank order).
rank 0: ghost left = 15, inner = 0..3, ghost right = 4
rank 1: ghost left = 3, inner = 4..7, ghost right = 8
rank 2: ghost left = 7, inner = 8..11, ghost right = 12
rank 3: ghost left = 11, inner = 12..15, ghost right = 0
The left ghost cell of each rank holds the value of the last cell of its left neighbor, and the right ghost cell holds the value of the first cell of its right neighbor. To the left of rank 0, the value 15 of the last rank arrived because of the periodic boundary.
A caution when writing send and receive separately
You can also write MPI_Send (send) and MPI_Recv (receive) separately. In that case, if all ranks call MPI_Send first, then when the values to send are large, everyone may end up in a state of “waiting until the partner receives” and be unable to proceed. This state is called deadlock. MPI_Sendrecv leaves sending and receiving together to MPI, so it avoids this kind of mutual waiting.
Exercise 2: Values that go into the ghost cells
When you run the sendrecv above with two processes, what values go into the left and right ghost cells of rank 0, and from which rank and with which tag does each arrive? Explain. Then predict the output when run with one process (mpiexec -np 1), and run it to check.
With two processes, rank 0 is responsible for global cells 0 to 3, and rank 1 for 4 to 7. Both the left neighbor and the right neighbor of rank 0 are rank 1. To the left ghost cell, the last cell 7, which rank 1 sent to the right in (1), arrives with tag 100. To the right ghost cell, the first cell 4, which rank 1 sent to the left in (2), arrives with tag 200. Even when the partner is the same rank, (1) and (2) are separate calls with different tags, so they are not mixed up.
rank 0: ghost left = 7, inner = 0..3, ghost right = 4
rank 1: ghost left = 3, inner = 4..7, ghost right = 0
With one process, both left and right are the rank itself (0), which sends to itself and receives from itself.
rank 0: ghost left = 3, inner = 0..3, ghost right = 0
Its own right-end value 3 goes into the left ghost cell, and its left-end value 0 goes into the right ghost cell. This is exactly a periodic boundary condition. As we will see in Section 7, MISO handles the periodic boundary in an undivided direction by this same “exchange with itself”.
6. Starting Communication and Waiting for It Later: MPI_Isend and MPI_Irecv
MISO’s halo exchange uses MPI_Isend, MPI_Irecv, and MPI_Waitall instead of MPI_Sendrecv. MPI_Isend and MPI_Irecv only start the communication and return without waiting for it to complete. This pattern is called non-blocking communication (the I in the name stands for “immediate”). On return, they write a handle for checking completion later into a variable of type MPI_Request. MPI_Waitall waits until the communications of all the handles passed to it have completed.
Rewriting (1) and (2) of the previous section with non-blocking communication gives the following.
MPI_Request reqs[4];
MPI_Isend(&u[n], 1, MPI_DOUBLE, right, 100, MPI_COMM_WORLD, &reqs[0]);
MPI_Irecv(&u[0], 1, MPI_DOUBLE, left, 100, MPI_COMM_WORLD, &reqs[1]);
MPI_Isend(&u[1], 1, MPI_DOUBLE, left, 200, MPI_COMM_WORLD, &reqs[2]);
MPI_Irecv(&u[n + 1], 1, MPI_DOUBLE, right, 200, MPI_COMM_WORLD, &reqs[3]);
// only started so far; do not read u[0], u[n+1] yet
MPI_Waitall(4, reqs, MPI_STATUSES_IGNORE);
Even if you replace the two MPI_Sendrecv calls in sendrecv.cpp with this, the output is the same.
We start four communications and then wait for them all together, so no mutual waiting caused by the order of sends and receives occurs. There are two rules to follow.
Do not modify the values being sent (
u[n],u[1]) until the send has completed.Do not read or write the places being received into (
u[0],u[n + 1]) until the receive has completed.
After MPI_Waitall, you can start the finite-difference computation using the values of the ghost cells.
7. MISO’s Halo Exchange
The overall structure of HaloExchanger
The halo exchange of the nine MHD variables is done by HaloExchanger in miso/include/miso/mhd_halo_exchange.hpp. Using the partial specialization we learned in Chapter 6, separate versions are written for backend::Host and backend::CUDA. Here we read the Host version. The CUDA version differs in that the packing of values into buffers is done by a GPU kernel. It is covered in Chapter 9.
The Host version has send and receive buffers for each direction and each orientation.
template <typename Real> struct HaloExchanger<Real, backend::Host> {
// MPI communication buffers
Array4D<Real> recv_x_pos, recv_x_neg;
// ...(same for y, z and send_*)
Grid<Real, backend::Host> &grid;
mpi::Shape &mpi_shape;
HaloExchanger(Grid<Real, backend::Host> &grid,
ExecContext<Real, backend::Host> &exec_ctx)
: recv_x_pos(grid.i_margin, grid.j_total, grid.k_total, n_fields),
// ...
The shape of the buffer in the x direction is (i_margin, j_total, k_total, n_fields). What is sent in the x direction is the i_margin layers of faces at the ends, packed into one array for all nine variables (n_fields = 9). Sending them together once reduces the number of communications compared with communicating nine times, once per variable.
Packing and sending the values
This is the part of apply for the x direction, positive orientation.
if (mpi_shape.x_procs_pos != MPI_PROC_NULL && grid.i_size > 1) {
for (int i = 0; i < grid.i_margin; ++i) {
for (int j = 0; j < grid.j_total; ++j) {
for (int k = 0; k < grid.k_total; ++k) {
for (int v = 0; v < 9; ++v) {
send_x_pos(i, j, k, v) =
(*vars[v])(grid.i_total - 2 * grid.i_margin + i, j, k);
}
}
}
}
MPI_Isend(send_x_pos.data(), send_x_pos.size(), mpi::data_type<Real>(),
mpi_shape.x_procs_pos, 100, mpi_shape.cart_comm,
&reqs[req_count++]);
MPI_Irecv(recv_x_pos.data(), recv_x_pos.size(), mpi::data_type<Real>(),
mpi_shape.x_procs_pos, 200, mpi_shape.cart_comm,
&reqs[req_count++]);
}
We explain how to read it step by step.
varsis an array of pointers to the nine variables (ro,vx, …,ph) created at the start ofapply.(*vars[v])(i, j, k)is the value at(i, j, k)of thev-th variable.The exchange is done only when there is a neighbor on the positive side (it is not
MPI_PROC_NULL) and the x direction has more than one cell (i_size > 1).What is packed is from
i_total - 2 * i_margintoi_total - i_margin - 1, that is, thei_marginlayers at the right end of the assigned range.send_x_pos.data()is the pointer to the start of the buffer, andsend_x_pos.size()is the number of elements.mpi::data_type<Real>()is a function template that returnsMPI_FLOATifRealisfloatandMPI_DOUBLEif it isdouble(miso/include/miso/mpi_util.hpp; the specialization example in Chapter 6).We send to the positive side with tag 100 and receive from the positive side with tag 200.
In the negative orientation, the range to pack is the i_margin layers starting at i_margin (the left end of the assigned range); we send to the negative side with tag 200 and receive from the negative side with tag 100. Tag 100 means “a value sent in the positive orientation”, and tag 200 means “a value sent in the negative orientation”. A value that a rank sends to the positive side with tag 100 is, from the receiver’s point of view, a value arriving from the negative side, so recv_x_neg receives it with tag 100. This is the same correspondence as (1) (send to the right with tag 100) and (2) (send to the left with tag 200) in the test program of the previous section.
req_count++ is a notation that uses the current req_count as the index and then increases it by 1. The handles of the started communications go into reqs[0], reqs[1], … in order.
Waiting for completion and putting the values back into the array
Once the communications in both orientations have been started, we wait for completion with MPI_Waitall and put the received values back into the ghost cells.
if (req_count > 0) {
MPI_Waitall(req_count, reqs, MPI_STATUSES_IGNORE);
}
req_count = 0;
// positive x-direction
if (mpi_shape.x_procs_pos != MPI_PROC_NULL && grid.i_size > 1) {
// ...(4-level loop over v, i, j, k)
(*vars[v])(grid.i_total - grid.i_margin + i, j, k) =
recv_x_pos(i, j, k, v);
The values received from the positive side go into the right-end ghost cells (i_total - i_margin to i_total - 1). The values received from the negative side go into the left-end ghost cells (0 to i_margin - 1).
The process of packing values into a buffer is called packing, and the process of putting them back from the buffer into the array is called unpacking. The names pack_x_send and unpack_x_recv in the CUDA version refer to these processes.
Why we exchange one direction at a time in the order x, y, z
apply completes the exchange in the x direction, then exchanges the y direction, and after that the z direction. The buffer in the y direction has the shape (i_total, j_margin, k_total, n_fields), and in the x direction it sends the whole range including the ghost cells (i_total). The ghost values filled by the x-direction exchange reach the diagonal neighboring ranks through the y-direction exchange, so correct values also go into the corners (cells that are ghosts in both x and y). The comment at the top of the code (Exchange one direction at a time ...) states this.
Periodic boundaries and undivided directions
As we saw in Section 3, in a periodic direction the neighbor of a rank at an end is the rank at the opposite end. If that direction is not divided (x_procs: 1), the neighbor is the rank itself, and the halo exchange becomes an exchange with itself (the same as the one-process case of Exercise 2). For this reason, MISO realizes the periodic boundary condition not with a boundary condition function but with the halo exchange. demo/mhd2d_vortex uses mhd::EmptyBoundaryCondition (a boundary condition that does nothing) because it is periodic in both x and y, and all the ghost cells are filled by the halo exchange.
The radiative transfer (RT) computation has a separate halo exchange (mpi_exchange_halo in miso/include/miso/rt_host.hpp). The idea is the same.
Exercise 3: The range to send and the range to receive
Consider a rank with i_total = 68 and i_margin = 2. Answer which cells (by local index i) this rank sends to the positive side in x, and into which cells (by index) the neighboring rank on the positive side puts the values it receives. Also explain, with the reason, which of recv_x_pos and recv_x_neg the receiving side uses with tag 100.
What is sent is i_total - 2 * i_margin + i (i = 0, 1), that is, i = 64, 65. These are the two layers at the right end of the assigned range (2 to 65).
From the point of view of the neighboring rank on the positive side, these values arrive from the negative side. It receives them with recv_x_neg and puts them into the left-end ghost cells i = 0, 1. Our 64 and 65 become the neighbor’s 0 and 1, so our last two cells line up immediately to the left of the neighbor’s first assigned cell (2).
The tag used when sending to the positive side is 100. The one that receives that value is, for the receiving side, communication with the negative side, so it is recv_x_neg that specifies tag 100. The sender says “send in the positive orientation” and the receiver says “receive from the negative side”; they are looking at the same communication from opposite standpoints.
8. Physical and Internal Boundaries
There are two kinds of ghost cells.
Internal boundary: a boundary that has a neighboring rank. The ghost cells are filled with the neighbor’s values (halo exchange). The ends of a periodic boundary also belong here.
Physical boundary: the edge of the computational domain itself, a non-periodic boundary. There is no neighboring rank (
MPI_PROC_NULL), so the ghost cells are filled by a boundary condition (symmetric, antisymmetric, and so on).
When the domain is divided, only the ranks at the ends face a physical boundary. If an inner rank overwrote the ghost cells with a boundary condition function, it would destroy values that should be received from the neighbor. So before applying a boundary condition, we check whether our own face is a physical boundary. This check is boundary_condition::is_physical_boundary (miso/include/miso/boundary_condition.hpp).
bool is_physical_boundary(const Direction direction, const Side side,
const mpi::Shape &mpi_shape) {
switch (direction) {
case Direction::X:
return (side == Side::INNER) ? (mpi_shape.x_procs_neg == MPI_PROC_NULL)
: (mpi_shape.x_procs_pos == MPI_PROC_NULL);
// ...(same for Y, Z)
}
return false;
}
Side::INNER is the side with smaller coordinates (the negative side), and Side::OUTER is the side with larger coordinates (the positive side). If the neighbor on that side is MPI_PROC_NULL, it is a physical boundary. condition ? A : B is an expression whose value is A if the condition is true and B if it is false.
The boundary conditions in demo/hd2d_kelvin_helmholtz/src/main.cpp, where the y direction is not periodic, are used as follows (excerpt).
void apply(mhd::FieldsView<Real> qq, GridView<const Real> grid) const {
namespace bc = miso::boundary_condition;
Backend btag{};
if (bc::is_physical_boundary(Direction::Y, Side::INNER, mpi_shape)) {
bc::symmetric(btag, qq.ro, grid, Sign::Pos, Direction::Y, Side::INNER);
// ...(same for the other 8 variables; the Side::OUTER side has the same form)
}
}
namespace bc = miso::boundary_condition; is a way to give a long namespace a short alias bc. The x direction is a periodic boundary, so no boundary condition is written for the x direction.
The order in which the boundary condition and the halo exchange are called is determined in apply_boundary_condition in miso/include/miso/mhd_integrator_host.hpp.
void apply_boundary_condition(const BoundaryCondition &bc,
Fields<Real, backend::Host> &qq) {
bc.apply(qq.view(), grid.const_view());
halo_exchanger.apply(qq);
}
First the ghost cells at the physical boundary are filled by the boundary condition, and then the ghost cells at the internal boundary are filled by the halo exchange. The values filled at the physical boundary in the y direction also reach the neighboring ranks through the halo exchange in the x direction (which sends the whole j_total). Where this function is called in the time evolution is read in Chapter 11.
In the output of the test program in Section 4, only ranks 0 and 2 at the bottom end in the y direction had y inner = 1. The x direction is a periodic boundary, so x inner = 0 on every rank.
9. Combining Values across All Ranks: MPI_Allreduce
Collective communication
Even if each rank finds its own candidate dt for the time step in its own region, the times would drift apart if each rank advanced with a different dt. We need to find the minimum of the candidates over all ranks and have everyone use it. Communication in which all processes of a communicator take part is called collective communication. MPI_Allreduce combines the values of all ranks with an operation (minimum, sum, and so on) and distributes the result to all ranks.
Test program: minimum and sum
Save it as allreduce.cpp.
#include <iostream>
#include <mpi.h>
int main(int argc, char **argv) {
MPI_Init(&argc, &argv);
int rank = 0, size = 0;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &size);
// each rank's time step candidate from its own region (dummy values here)
double dt_local = 0.1 / (rank + 1);
double dt_global = 0.0;
MPI_Allreduce(&dt_local, &dt_global, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD);
int one = 1;
int count = 0;
MPI_Allreduce(&one, &count, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
std::cout << "rank " << rank << ": dt_local = " << dt_local
<< ", dt_global = " << dt_global << ", count = " << count
<< std::endl;
MPI_Finalize();
return 0;
}
The arguments are, in order, the address of the value to send, the address to receive the result, the count, the type, the operation, and the communicator. If you change MPI_MIN to MPI_SUM, you get the sum. In the second call, we add up the 1 of each rank to count the number of ranks.
mpicxx -std=c++17 -o allreduce allreduce.cpp
mpiexec -np 4 ./allreduce
rank 0: dt_local = 0.1, dt_global = 0.025, count = 4
rank 1: dt_local = 0.05, dt_global = 0.025, count = 4
rank 2: dt_local = 0.0333333, dt_global = 0.025, count = 4
rank 3: dt_local = 0.025, dt_global = 0.025, count = 4
dt_local differs from rank to rank, but dt_global became the same minimum value on all ranks.
The dt in Integrator::cfl
Integrator::cfl, which we read in Chapter 7 (miso/include/miso/mhd_integrator_host.hpp), finds dt for the CFL condition using only the cells in its own assigned range, and at the end takes the minimum over all ranks.
Real dt_global;
MPI_Allreduce(&dt, &dt_global, 1, mpi::data_type<Real>(), MPI_MIN,
mpi::comm());
return dt_global;
The reduce of Chapter 7 combines the values of an array within one process, and MPI_Allreduce further combines the results of each process. It is easier to read if you think of it as finding the global value in two stages.
All ranks call in the same order
Collective communication must be called by all ranks of the communicator, in the same order. If only one rank does not call it, or calls a different collective first, the other ranks keep waiting for a partner who never comes, and stop.
ModelBase::update in miso/include/miso/mhd_model_base.hpp has a device to keep this convention. It handles the case where a source term has an upper limit dt_limit() on the time step.
Real dt_src = static_cast<Real>(d.src.dt_limit());
if (!(dt_src > Real(0))) {
dt_src = Real(-1);
}
Real dt_src_g;
MPI_Allreduce(&dt_src, &dt_src_g, 1, mpi::data_type<Real>(), MPI_MIN,
mpi::comm());
if (dt_src_g <= Real(0)) {
throw std::runtime_error(/* ...(error description)... */);
}
Even if dt_limit() returns a non-positive (invalid) value on some rank, that rank does not throw an exception immediately; it sets the value to -1 and takes part in the MPI_Allreduce. Taking the minimum makes every rank receive -1, so all ranks throw the exception at the same time. This prevents the situation in which only one rank stops first while the other ranks keep waiting in MPI_Allreduce.
Elsewhere in MISO, the constructor of Config (miso/include/miso/config.hpp) uses MPI_Bcast (distributing the contents of the configuration file read by rank 0 to all ranks), and Shape::save uses MPI_Gather (collecting the coord of all ranks on rank 0).
Exercise 4: Conditions for calling collective communication
Suppose you moved the MPI_Allreduce in Integrator::cfl inside if (mpi::is_root()) { ... }. Explain what happens when you run with two processes (you do not need to actually rewrite and run it).
Only rank 0 calls MPI_Allreduce, and rank 1 does not. Rank 0 keeps waiting for the value from rank 1 to arrive, and the computation does not proceed. Rank 1 returns the uninitialized dt_global and moves on, and then in the next halo exchange or collective communication it ends up waiting for a communication that does not correspond to anything on rank 0. In either case, the convention that all ranks call the same collective communication in the same order is broken, and the program either hangs or gives incorrect results.
What may be done on rank 0 alone is processing that involves no communication, such as writing to a file or displaying on the screen.
10. Running a MISO Demo with MPI
app_run.sh and the OpenMPI environment variables
The basics of building and running a demo were covered in Chapter 1. Looking at the execution part of demo/mhd2d_vortex/app_run.sh shows how to launch with MPI (excerpt).
# Number of processes
NUM_PROCS=2
# ...
# Load environment variables
set -a
source "${THIS_DIR}"/../shared/openmpi_config.env
set +a
# ...
mpiexec -np ${NUM_PROCS} "${APP_NAME}" --config="${CONFIG_PATH}"
NUM_PROCS=2 matches the product of x_procs: 2, y_procs: 1, and z_procs: 1 in config.yaml. If you change the division in the configuration file, change this number accordingly.
Variables loaded between set -a and set +a are passed on to mpiexec as environment variables. The contents of the file loaded, demo/shared/openmpi_config.env, are as follows.
# Environment variables for OpenMPI
# * Suppress errors when using MPI in containers
# Usage (shell): set -a && source ./openmpi_config.env && set +a
# Usage (docker): docker run --env-file openmpi_config.env
OMPI_MCA_pml=ob1
OMPI_MCA_btl=tcp,self,vader
OMPI_MCA_coll=^hcoll
OMPI_MCA_osc=^ucx
OMPI_MCA_hwloc_base_binding_policy=none
Environment variables starting with OMPI_MCA_ specify how to choose the internal components of OpenMPI (MCA: Modular Component Architecture). As the comments say, these settings avoid errors in environments, such as inside containers, where certain communication libraries cannot be used. pml and btl specify the layer and the paths for one-to-one communication (TCP, itself, and the shared memory vader within the same computer), coll specifies the components for collective communication, and osc those for one-sided communication; ^ means “exclude”. The last line specifies that processes are not pinned to particular CPU cores.
Depending on the environment, the following two may also be needed. This is why they are attached when running MISO’s unit tests with ctest.
OMPI_MCA_btl_vader_single_copy_mechanism=none: avoids warnings and errors in environments (such as containers) where the feature of copying memory directly between processes on the same computer cannot be used.OMPI_MCA_rmaps_base_oversubscribe=yes: allows starting more processes than the number of cores on the computer.
Both are OpenMPI-only settings and have no effect with other MPI implementations (such as MPICH). Detailed descriptions of each component are in the OpenMPI documentation (Open MPI 4.1 FAQ).
Launching by hand
You can also launch by hand without using app_run.sh. After building the demo, run the following in the MISO folder.
cd /path/to/MISO
set -a; source demo/shared/openmpi_config.env; set +a
mpiexec -np 2 demo/mhd2d_vortex/build/mhd2d_vortex --config=demo/mhd2d_vortex/config.yaml
--config= specifies the configuration file. The output destination io.save_dir (default data/) is a path relative to the folder containing the configuration file, so output is written to demo/mhd2d_vortex/data/. For a short trial, copy the configuration file to another folder, make time.tend small (for example 0.02), and pass that file to --config=. The output is then written to data/ in the folder you copied to. This is the screen display when run with tend: 0.02.
time = 0.00; n_step = 0; n_output = 0
time = 0.01; n_step = 4; n_output = 1
time = 0.02; n_step = 9; n_output = 2
Only rank 0 produces the display.
Output files and ranks
mhd/mhd.<n_output>.<rank>.bin under data/ (example: mhd.00000001.00000000.bin) holds the MHD variables of each rank. Each rank writes only its own array (including ghost cells). grid.bin is the global grid, written by rank 0 only. mpi/coords.csv is the correspondence between ranks and positions that we saw in Section 3; with x_procs: 2, two lines, 0,0,0,0 and 1,1,0,0, appear. pymiso.Data reads coords.csv and the file of each rank, removes the ghost cells, and assembles the global arrays (such as ro). Readers can analyze the data without being aware of how the domain was divided.
Comparing results with different divisions
Domain decomposition does not change the procedure of the computation, so the result should not change even if the way of dividing changes. This can also be used as a way to check your work when you have modified the MPI part. Suppose you set demo/mhd2d_vortex to tend: 0.02, set the mpi section to three cases, \(1 \times 1\) (-np 1), \(2 \times 1\) (-np 2), and \(2 \times 2\) (-np 4), and ran each in a separate folder. Use pymiso to examine the difference of the values at output number 2.
python3 - run_1x1/data run_2x1/data run_2x2/data <<'EOF'
import sys
import numpy as np
import pymiso
ref = pymiso.Data(sys.argv[1])
ref.load(2)
for d in sys.argv[2:]:
da = pymiso.Data(d)
da.load(2)
for v in ["ro", "vx", "bx", "ei"]:
print(d, v, float(np.max(np.abs(getattr(ref, v) - getattr(da, v)))))
EOF
run_1x1 and so on are the folders holding the respective configuration files. In the result we checked, the difference was 0.0 for every variable, and the results were identical bit for bit. This is because the halo exchange just copies values as they are, and the minimum by MPI_MIN also gives the same value regardless of the order of computation. If there is a step that takes a floating-point sum with MPI_SUM, the order of addition changes with the division, so differences of about the size of rounding error may appear.
Exercise 5: Running with a division in the y direction
Copy the configuration file of demo/mhd2d_vortex, change tend to 0.02 and the mpi section to x_procs: 1, y_procs: 2, and run with two processes. Check the contents of data/mpi/coords.csv, and compare with the result run with x_procs: 2, y_procs: 1 using the method above. Also compute the size (in bytes) of the per-rank output file mhd.00000000.00000000.bin and compare it with the actual size. The file begins with a 4-byte integer giving the element size, followed by the values of the 9 variables including ghost cells as Real = float (4 bytes).
coords.csv looks as follows. The ranks are lined up in the y direction, so the y coordinate of rank 1 is 1.
rank,x,y,z
0,0,0,0
1,0,1,0
The difference from the result with x_procs: 2 is 0.0 for every variable.
The grid of one rank has no division in the x direction, so i_total = 128 + 2 * 2 = 132; in the y direction j_total = 64 + 2 * 2 = 68; and in the z direction k_size = 1, so there are no ghost cells and k_total = 1.
This matches the size displayed by ls -l data/mhd. When x_procs: 2, i_total = 68 and j_total = 132, so the size is the same.
MPI is not the name of a particular piece of software but a standard that defines the names, arguments, and behavior of functions. The standard is defined by the MPI Forum, and the latest document is published on the MPI Forum web page. OpenMPI and MPICH are libraries that implement this standard. The meaning of the functions described in the standard (MPI_Sendrecv, MPI_Allreduce, and so on) is the same in every implementation, but the concrete value of MPI_PROC_NULL and run-time settings such as OMPI_MCA_ differ between implementations. The -2 displayed for MPI_PROC_NULL in this chapter’s output is the value in OpenMPI 4.1.
The functions used from C++ are the C-language ones. The MPI 2.x standard also had C++ functions (such as MPI::COMM_WORLD.Get_rank()), but they were removed in MPI 3.0. Current C++ programs call the C functions, as MISO does.
11. Summary
MPI runs the same program in multiple processes, and the processes are distinguished by rank. Values are sent and received explicitly between processes.
In MISO, the constructor and destructor of
Envtake care ofMPI_InitandMPI_Finalize. The rank is obtained withmpi::rank().The number of divisions is set in the
mpisection ofconfig.yaml, andmpi::Shapefinds the arrangement of ranks and the neighboring ranks withMPI_Cart_createandMPI_Cart_shift. If there is no neighbor, it isMPI_PROC_NULL.Gridis the grid of its own rank. It hasi_marginghost cells on both sides of the assigned cells, and the count including ghost cells isi_total.HaloExchangerpacks the nine variables into buffers, exchanges with the neighbors usingMPI_Isend/MPI_Irecv/MPI_Waitall, and fills the ghost cells one direction at a time in the order x, y, z. The periodic boundary is also handled by this exchange.Whether a boundary is physical is checked with
is_physical_boundary, and only the ghost cells at physical boundaries are filled by the boundary condition.The time step is synchronized to the minimum over all ranks with
MPI_Allreduce(MPI_MIN). All ranks call collective communication in the same order.
The halo exchange when using a GPU, and the mapping between ranks and GPUs, are covered in Chapter 9. If you proceed with the CPU only, in the heat conduction Practice of Chapter 10 you will write the one-dimensional halo exchange and the MPI_Allreduce of this chapter yourself.