02. MPI Version: Dividing the Domain

You extend the program of the 01. serial version into heat1d_mpi, which shares the work among \(P\) ranks. The ideas of domain decomposition and ghost-cell exchange are explained in Chapter 10, and the MPI functions in Chapter 8.

What you use

In addition to what the serial version uses, you use the following.

Item

Header

miso::Env (MPI_Init and MPI_Finalize), miso::mpi::rank(), mpi::size(), mpi::comm()

miso/env.hpp

MPI_Sendrecv、MPI_Allreduce、MPI_Gather、MPI_Abort、MPI_PROC_NULL

mpi.h

If you create Env env(argc, argv); at the start of main, MPI is initialized, and the finalization is performed when you leave main.

What you build

The grid of each rank

  • Divide the total number of cells \(N\) by the number of ranks \(P\), and let the number of physical cells of each rank be \(n = N/P\). If \(N\) is not divisible by \(P\), rank 0 prints the reason and ends all ranks with MPI_Abort.

  • The arrays are Array1D<Real> u(n + 2) and u_new(n + 2).

  • Index i (1 to n) of the array on rank \(p\) corresponds to the global index \(j = p\,n + (i - 1)\). In the initial condition and the error computation, obtain \(x_j = (j + 1/2)\Delta x\) from this \(j\).

  • Obtain \(\Delta x\), \(n_\mathrm{steps}\), \(\Delta t\), and \(r\) from the global \(N\) with the same expressions as in the serial version.

Neighboring ranks

  • The left neighbor is rank - 1, and the right neighbor is rank + 1.

  • Set the left neighbor of rank 0 and the right neighbor of rank \(P - 1\) to MPI_PROC_NULL.

Each step

  1. Fill the ghost cells with the following two calls to MPI_Sendrecv. The type is MPI_DOUBLE, the number of elements is 1, and the communicator is mpi::comm().

    Call

    Value to send → partner

    Where to receive ← partner

    Tag

    First

    u[1] → left neighbor

    u[n+1] ← right neighbor

    0

    Second

    u[n] → right neighbor

    u[0] ← left neighbor

    1

    The receive status is not used, so pass MPI_STATUS_IGNORE as the last argument.

  2. A rank whose left neighbor is MPI_PROC_NULL sets u[0] = -u[1], and a rank whose right neighbor is MPI_PROC_NULL sets u[n+1] = -u[n].

  3. With the same expression as in the serial version, update the physical cells (Range1D{1, n + 1}) and swap u and u_new.

Error and output

  • On each rank, use reduce to compute the maximum error over its own physical cells, and turn it into the global maximum error with MPI_Allreduce (MPI_MAX).

  • With MPI_Gather, collect the \(n\) physical cells of each rank (starting from &u[1]) on rank 0. Rank 0 prepares a receive array (std::vector<Real>) of \(N\) elements.

  • Only rank 0 prints the same 3 lines as the serial version (the first line also shows \(P\), as in N = 64, P = 2, steps = 1024, r = 0.4000), and writes the result to result_N<N>_P<P>.txt (for example result_N64_P2.txt).

What to check

In the same folder where you ran the serial version, run the following.

./build/heat1d_serial 64
for P in 1 2 4; do
  mpiexec -np $P ./build/heat1d_mpi 64
  cmp result_N64_serial.txt result_N64_P$P.txt
done

  1. Printed output: For any of \(P = 1, 2, 4\), the two printed errors have the same values as in the serial version.

  2. Bit-for-bit agreement: Check that cmp prints nothing (the files do not differ by even one byte).

  3. Non-divisible decomposition: With mpiexec -np 3 ./build/heat1d_mpi 64, check that the program prints the reason and exits.

  4. Three ranks: Run \(P = 1, 2, 3, 4\) with \(N = 120\), and check that all agree with the serial result_N120_serial.txt. In this case, max error (exact) is 2.941414e-05.

If it does not agree with the serial version, check in the following order.

  • If it differs even for \(P = 1\), check \(x_j\) in the initial condition (how you compute the global index \(j\)) and how you compute \(\Delta t\).

  • If it agrees for \(P = 1\) but differs for \(P = 2\), compare the value to send, the place to receive, the partner, and the tag of MPI_Sendrecv with the table above.

  • If only the end cells differ, check the condition for putting in the physical-boundary values (which rank puts them in).

You can find the line where the files first differ with diff result_N64_serial.txt result_N64_P2.txt | head. From the line number you get the global index \(j\) (line number − 1), which narrows down which end of which rank has the wrong value.

Criteria for moving on

Once the results match the serial version bit for bit for several numbers of ranks, move on to the 03. GPU version. Even without a GPU, you can go as far as the CPU build of Assignment 3.