10. Heat Conduction Calculation

In this chapter, you combine the arrays, for_each and reduce, MPI, and GPUs you have learned so far into a 1D heat conduction calculation. The chapter explains how to solve the equation with finite differences and how to check that the result is correct. You write the program yourself, following the specification of the Practice heat1d.

You build the program in the following three stages.

1. Serial version: solve the problem on the CPU in a single process, using Array1D and the for_each and reduce of backend::Host.

2. MPI version: divide the domain and exchange ghost cell values with the neighboring ranks using MPI_Sendrecv.

3. GPU version: replace the arrays with Array1D<Real, backend::CUDA> and the loops with the for_each of backend::CUDA.

Stages 1 and 2 can be completed on the CPU of an Ubuntu or Mac machine. For stage 3, if you have no GPU, you can build the same source for the CPU and still check that the structure is correct.

Contents

• The heat equation and its analytical solution

• Grid and finite differences

• Stability condition of the explicit method

• Boundary conditions and ghost cells

• How to compare the results with expected values

• Building the serial version

• MPI version: domain decomposition and ghost cell exchange

• GPU version: moving arrays and loops to the GPU

• Correspondence with MISO itself

The heat equation and its analytical solution

The way the temperature \(u(x, t)\) diffuses and becomes smooth is described by the following heat equation.

\[\frac{\partial u}{\partial t} = \kappa\,\frac{\partial^2 u}{\partial x^2}\]

\(\kappa\) is a positive constant (the thermal diffusivity). In this chapter, we solve the following problem.

Item	Setting
Domain	\(0 < x < 1\)
Boundary condition	\(u = 0\) at both ends (fixed temperature)
Initial condition	\(u(x, 0) = \sin(\pi x)\)
Constant	\(\kappa = 1\)
End time	\(t_\mathrm{end} = 0.1\)

This problem has a solution that can be written as a formula (the analytical solution).

\[u(x, t) = e^{-\pi^2 \kappa t}\sin(\pi x)\]

Substituting \(u = e^{-\pi^2\kappa t}\sin(\pi x)\) into the equation makes both the left and right sides equal to \(-\pi^2\kappa u\), and it also satisfies \(u = 0\) at both ends. The shape stays \(\sin(\pi x)\), and only the amplitude decays exponentially. The amplitude at \(t = 0.1\) is \(e^{-\pi^2 \times 0.1} \approx 0.3727\).

Because an analytical solution exists, you can measure numerically how accurate the computed values are. For the MHD problems solved with MISO, such a solution usually does not exist. We therefore first check that the parts of the program (finite differences, boundaries, communication) are correct on a simple problem that has an analytical solution.

Grid and finite differences

Cells and grid spacing

We divide the domain \([0, 1]\) into \(N\) equal intervals (cells). The cell width (the grid spacing) is \(\Delta x = 1/N\). Values are stored at the cell centers. The center of the \(j\)-th cell (\(j = 0, 1, \dots, N-1\)) is at the following position.

\[x_j = \left(j + \tfrac12\right)\Delta x\]

The MISO grid also stores values at cell centers. The boundaries (\(x = 0\) and \(x = 1\)) lie on the faces of the end cells.

Central difference for the second derivative

We approximate the second derivative \(\partial^2 u/\partial x^2\) using differences between neighboring values. Letting \(u_j\) be the value in cell \(j\), we use the following central difference.

\[\left.\frac{\partial^2 u}{\partial x^2}\right|_{x_j} \approx \frac{u_{j-1} - 2u_j + u_{j+1}}{\Delta x^2}\]

A Taylor expansion shows that the error of this approximation is proportional to \(\Delta x^2\) (second-order accuracy). Because it uses one neighbor on each side, the end cells (\(j = 0\) and \(j = N-1\)) each need one value from outside the domain. We deal with this in Boundary conditions and ghost cells.

Explicit method in time

We discretize time as \(t^n = n\,\Delta t\) and write the value at time \(t^n\) as \(u_j^n\). Approximating the time derivative with the forward difference \((u_j^{n+1} - u_j^n)/\Delta t\) gives the following update formula.

\[u_j^{n+1} = u_j^n + r\left(u_{j-1}^n - 2u_j^n + u_{j+1}^n\right),\qquad r = \frac{\kappa\,\Delta t}{\Delta x^2}\]

The right-hand side can be computed from the values at time \(n\) alone, so the new values can be obtained one after another. Such a method is called an explicit method, and this approximation in time is called the forward Euler method. \(r\) is a dimensionless number called the diffusion number, and it appears in the stability condition in the next section.

In Fortran, one step looks like the following (an example for illustration).

do j = 1, n
  u_new(j) = u(j) + r * (u(j-1) - 2.0d0 * u(j) + u(j+1))
end do
u(1:n) = u_new(1:n)

The u read on the right-hand side and the u_new written on the left-hand side are separate arrays. If you overwrite the same array, then when you compute u(j+1) after updating u(j), you read the already updated u(j). The cells are then no longer all computed from values at the same time, and the calculation differs from the update formula above. Separating the array you read from the array you write is necessary in both serial and parallel code. It matters even more in parallel execution, because the order in which cells are updated is not determined.

Stability condition of the explicit method

In an explicit method, if \(\Delta t\) is too large, the values grow while oscillating and the calculation breaks down. Let us find how large it can be.

Consider the shortest wave, that is, the distribution \(u_j = (-1)^j\) whose sign alternates from cell to cell. Since \(u_{j-1} - 2u_j + u_{j+1} = -4u_j\), one step gives the following.

\[u_j^{n+1} = (1 - 4r)\,u_j^n\]

This component is multiplied by \((1 - 4r)\) at every step, so we need \(|1 - 4r| \le 1\), that is,

\[r = \frac{\kappa\,\Delta t}{\Delta x^2} \le \frac12\]

Otherwise this component keeps growing. This is the stability condition for the forward Euler method combined with the central difference. Components of other wavelengths do not grow as fast as this one, so this condition makes the whole solution stable. The factor applied in one step is called the amplification factor. We derive it for a general wavelength in How to compare the results with expected values.

If you halve the grid spacing, the upper limit of \(\Delta t\) becomes one quarter. To reach the same end time, the number of steps becomes 4 times larger and the number of cells 2 times larger, so the amount of computation becomes 8 times larger.

In the Practice, we set the upper limit of \(r\) to \(r_\mathrm{max} = 0.4\) and choose \(\Delta t\) so that the run ends exactly at the end time.

\[n_\mathrm{steps} = \left\lceil \frac{t_\mathrm{end}}{r_\mathrm{max}\,\Delta x^2/\kappa} \right\rceil,\qquad \Delta t = \frac{t_\mathrm{end}}{n_\mathrm{steps}}\]

\(\lceil\ \rceil\) means rounding up (std::ceil in C++). This keeps \(r \le r_\mathrm{max}\).

Exercise 10.1 For \(N = 64\) and \(r_\mathrm{max} = 0.4\), find \(n_\mathrm{steps}\) and \(r\). Also, explain from the amplification factor above what happens if you set \(r_\mathrm{max} = 0.6\).

::

Since \(\Delta x = 1/64\) and \(\Delta x^2 = 1/4096\), we have \(r_\mathrm{max}\Delta x^2/\kappa = 0.4/4096\). Then \(n_\mathrm{steps} = \lceil 0.1 \times 4096 / 0.4 \rceil = 1024\), \(\Delta t = 0.1/1024\), and \(r = 0.4\).

With \(r_\mathrm{max} = 0.6\), we get \(n_\mathrm{steps} = 683\) and \(r \approx 0.5997\), and the amplification factor of the sign-alternating component is \(1 - 4r \approx -1.4\). The initial condition \(\sin(\pi x)\) contains almost none of this component, but it is always mixed in as rounding error (about \(10^{-16}\)). It is multiplied by 1.4 at every step, and after 683 steps it has grown by a factor of \(1.4^{683} \approx 10^{100}\), so the calculation breaks down. In fact, if you give r = 0.6 to the serial version of the Practice, the maximum error becomes about \(10^{82}\).

Boundary conditions and ghost cells

Ghost cells

The difference at the end cell \(j = 0\) needs \(u_{-1}\), and the difference at \(j = N-1\) needs \(u_N\). So we add one extra cell at each end of the array and fill it with a value determined by the boundary condition. These extra cells are called ghost cells. The MISO grid follows the same idea and has i_margin ghost cells (Chapter 8).

The array indices are assigned as follows. C++ arrays start at 0, so the physical cell \(j\) is stored at the array index i = j + 1.

array index i :  0  |  1   2   ...   N  |  N+1
               ghost|  N physical cells | ghost
cell index j :      |  0   1   ...  N-1  |

If the ghost cell values are filled in beforehand, the update of every physical cell can be written with the same formula,

u_new[i] = u[i] + r * (u[i-1] - 2 u[i] + u[i+1])     (i = 1, ..., N)

so you no longer need to write a separate formula for the end cells.

The boundary condition \(u = 0\)

The boundary \(x = 0\) lies exactly midway between the ghost cell (center \(x = -\Delta x/2\)) and the first physical cell (center \(x = \Delta x/2\)). To make the straight line connecting the two values equal 0 at the midpoint, put the sign-reversed value into the ghost cell.

\[u_\mathrm{ghost,\,left} = -u_0,\qquad u_\mathrm{ghost,\,right} = -u_{N-1}\]

In terms of array indices, this is u[0] = -u[1] and u[N+1] = -u[N]. The ghost cell values must be refilled every time the physical cell values change. Fill the ghost cells before the finite difference of each step.

Exercise 10.2 For the initial condition \(\sin(\pi x)\), putting \(-u_0\) into the ghost cell gives exactly the value of \(\sin(\pi x)\) evaluated at \(x = -\Delta x/2\). Verify this. Verify the same for the right end.

::

Left end: \(\sin(\pi \cdot (-\Delta x/2)) = -\sin(\pi\Delta x/2) = -u_0\) (\(\sin\) is an odd function).

Right end: the center of the ghost cell is \(x = 1 + \Delta x/2\). Since \(\sin(\pi + \theta) = -\sin\theta\) and \(\sin(\pi - \theta) = \sin\theta\), we have \(\sin(\pi(1 + \Delta x/2)) = -\sin(\pi\Delta x/2) = -\sin(\pi(1 - \Delta x/2)) = -u_{N-1}\).

In other words, \(\sin(\pi x)\) can be expressed by the same formula including the ghost cells. In the next section, we use this to derive the result of the finite difference as a formula.

How to compare the results with expected values

We check the correctness of the computed values with the following three comparisons. Each has a different purpose, so the criterion for agreement also differs.

Comparison	What it checks	Expected difference
Comparison with the discrete solution	Whether the finite difference, boundary, and time integration are computed exactly as the formulas say	About the rounding error
Comparison with the analytical solution and the order of convergence	Whether the finite-difference approximation has the correct accuracy	About 1/4 when the number of grid points is doubled
Comparison across different numbers of ranks and compute locations	Whether the domain decomposition, communication, and the move to the GPU are correct	0 (bitwise identical) or about the rounding error

Comparison with the discrete solution

Applying the central difference to \(\sin(\pi x_j)\) and using the trigonometric addition formula \(\sin(a - b) + \sin(a + b) = 2\sin a\cos b\), we get

\[\sin(\pi x_{j-1}) - 2\sin(\pi x_j) + \sin(\pi x_{j+1}) = -4\sin^2\!\left(\frac{\pi\Delta x}{2}\right)\sin(\pi x_j)\]

so the shape stays \(\sin(\pi x_j)\) and only a constant factor is applied. As Exercise 10.2 shows, this relation also holds at the boundaries, so one step of the update formula multiplies the values in all cells by the same factor

\[G = 1 - 4r\sin^2\!\left(\frac{\pi\Delta x}{2}\right)\]

This is the amplification factor for the \(\sin(\pi x)\) component. The value after \(n\) steps can be computed as follows.

\[u_j^n = G^n \sin(\pi x_j)\]

We call this the discrete solution. If the finite difference and time integration are computed exactly as the formulas say, the only difference between the numerical solution and the discrete solution is rounding error. In double precision (double) with a few thousand steps, the maximum difference over all cells is about \(10^{-12}\) or less. If there is a mistake, such as an index shifted by one or a forgotten ghost cell fill, this difference becomes much larger than the rounding error.

Comparison with the analytical solution and the order of convergence

The difference from the analytical solution includes the approximation error of the finite difference (the truncation error). We compute the maximum difference over all cells,

\[E_N = \max_j \left| u_j^{n_\mathrm{steps}} - e^{-\pi^2\kappa t_\mathrm{end}}\sin(\pi x_j)\right|\]

doubling \(N\) each time. The scheme is second-order accurate in space and first-order accurate in time, but since we set \(\Delta t \propto \Delta x^2\), the error in time is also proportional to \(\Delta x^2\). Therefore, doubling \(N\) should reduce the error to about 1/4. From the ratio of the errors on two grids, we obtain the order of convergence \(p\) with the following formula.

\[p = \log_2\frac{E_N}{E_{2N}}\]

If \(p \approx 2\), the finite-difference approximation has the expected accuracy. The values under the Practice conditions (\(r_\mathrm{max} = 0.4\)) are as follows.

\(N\)	\(n_\mathrm{steps}\)	\(E_N\)	\(E_{N/2}/E_N\)	\(p\)
32	256	\(4.1368 \times 10^{-4}\)	—	—
64	1024	\(1.0341 \times 10^{-4}\)	4.0004	2.0001
128	4096	\(2.5852 \times 10^{-5}\)	4.0001	2.0000
256	16384	\(6.4630 \times 10^{-6}\)	4.0000	2.0000
512	65536	\(1.6158 \times 10^{-6}\)	4.0000	2.0000

The error is much larger than the rounding error, so we do not require agreement at the rounding-error level as in the comparison with the discrete solution. What we check is whether the error decreases at a fixed rate as the grid is refined.

Exercise 10.3 Compute the order of convergence \(p\) from \(E_N\) for \(N = 32\) and \(N = 64\) in the table. Also, estimate from the table roughly what \(N\) is needed to bring the error down to \(10^{-6}\) or less.

::

\(E_{32}/E_{64} = 4.1368\times10^{-4} / 1.0341\times10^{-4} \approx 4.0004\), and \(p = \log_2 4.0004 \approx 2.000\).

Since \(E_N \propto N^{-2}\), going from \(E_{512} = 1.6\times10^{-6}\) to \(10^{-6}\) requires multiplying \(N\) by \(\sqrt{1.6} \approx 1.27\), that is, \(N \approx 650\). The number of steps then grows in proportion to \(N^2\), and the amount of computation grows in proportion to \(N^3\).

::

Expanding the difference in amplitude between the discrete solution \(G^n\sin(\pi x_j)\) and the analytical solution in powers of \(\Delta x\), the leading term is as follows (for \(\kappa = 1\)).

\[E_N \approx t_\mathrm{end}\,\pi^4\,\Delta x^2\left(\frac{r}{2} - \frac{1}{12}\right)e^{-\pi^2 t_\mathrm{end}}\]

The term \(-1/12\) comes from the error of the central difference in space, and \(r/2\) from the error of the forward Euler method. Substituting \(r = 0.4\) and \(N = 32\) gives \(4.136\times10^{-4}\), which agrees well with the table above.

The two terms have opposite signs, so with \(r = 1/6\) the leading terms cancel and the error becomes proportional to \(\Delta x^4\). If you give r = 0.16666666666666667 to the serial version of the Practice, the error drops to about \(1.6\times10^{-8}\) at \(N = 64\) (because \(n_\mathrm{steps}\) is rounded up, \(r\) deviates slightly from 1/6, and you do not get a clean fourth-order behavior). This cancellation is specific to this combination of equation and difference scheme, and you cannot expect it for general problems.

Comparison across different numbers of ranks and compute locations

In the MPI version, the global grid, the initial condition, \(\Delta t\), and the update formula are the same regardless of the number of ranks. The update of each cell is determined only by that cell and its left and right neighbors, and the order of the computation (the order of the additions) does not change either. Therefore, as long as the correct values reach the ghost cells, the results are bitwise identical even if you change the number of ranks. They also match the serial version.

To check agreement, compare files in which the values of all cells are written with 17 digits. If the files differ by even one byte, cmp shows the first place where they differ. If nothing is printed, they match.

cmp result_N64_serial.txt result_N64_P4.txt

Comparing only the maximum error can miss a discrepancy. Compare the values of all cells.

When comparing the GPU version with the CPU version, bitwise agreement is not guaranteed. The GPU compiler may use an instruction that computes a multiplication and an addition with a single rounding (FMA, fused multiply-add), and the implementations of mathematical functions such as sin also differ from those on the CPU. The last digit may differ, so look at the maximum difference over all cells and regard the results as matching if it is about \(10^{-12}\) or less relative to the size of the values. You can compute the maximum difference with awk, for example, as follows (in both files, column 1 is \(x\) and column 2 is \(u\)).

paste result_N64_P1.txt result_N64_P1_gpu.txt \
  | awk '{d = $2 - $4; if (d < 0) d = -d; if (d > m) m = d} END {print m}'

How to compare computed results is also summarized in the reference material “Testing and Comparing Results”.

Building the serial version

In Practice 01, you write the formulas so far using MISO’s Array1D, for_each, and reduce. The basics of arrays and loops were explained in Chapter 7. Here we explain the points to watch for in this Practice.

Value type

Use using Real = double; to compute in double precision. The MISO demos run with using Real = float; (single precision), but the rounding error of single precision is about \(10^{-7}\), which is too large to resolve the comparison with the discrete solution (about \(10^{-12}\)) or the error at \(N = 512\) (about \(10^{-6}\)).

Arrays and views

Allocate Array1D<Real> u(N + 2) for \(N\) physical cells and 2 ghost cells. Allocate the destination u_new separately with the same size.

Inside the lambda passed to for_each, use a view obtained with view(), not the Array1D itself. MISO_LAMBDA captures variables by value ([=]), so using an Array1D directly tries to copy an Array1D, whose copying is prohibited, and causes a compile error.

error: use of deleted function 'miso::Array1D<T, miso::backend::Host>::Array1D(const miso::Array1D<T, miso::backend::Host>&) [with T = double]'

A view is a lightweight object that holds only the location of the data (a pointer) and its size, so copying it does not duplicate the contents of the array (Chapter 4).

Flow of one step

Proceed through one step in the following order.

1. Fill the ghost cells with the boundary condition values (u[0] = -u[1], u[N+1] = -u[N]).

2. For the physical cells (indices 1 to N), read u and write to u_new (for_each, with the range Range1D{1, N + 1}).

3. Swap u and u_new.

In the serial version, writing the two ghost cell values directly as u[0] = ... works. However, when you later move to the GPU, you cannot write to an array on the GPU by index from the host. If you write them from the start inside a for_each over Range1D{0, 1} (a loop that runs only once), you can carry the code over to the GPU version as it is.

For the swap in step 3, you can use std::swap(u, u_new); (<utility>). Array1D supports move (reassigning the pointers to the contents), so the contents of the two arrays are exchanged without copying elements. A view still points to the data from before the swap, so after the swap, obtain the views again in the next step. If you call view() every time inside the step loop, this mistake does not happen.

Finding the maximum error with reduce

Find the maximum error over all cells with reduce. Pass the initial value Real(0), a lambda that returns the value of each cell, and a lambda that returns the larger of two values. The error is 0 or greater, so Real(0) serves as the identity element of the maximum (Chapter 7).

auto max_op = MISO_LAMBDA(Real a, Real b) { return a > b ? a : b; };

In the lambda that returns the value of each cell, compute \(x_j\) and the analytical solution (or the discrete solution), and return the absolute value of the difference. If you use miso::util::sin, miso::util::exp, and miso::util::fabs (miso/include/miso/utility.hpp) for the mathematical functions, you can write the same code for the GPU version. The constant pi is miso::pi<Real> (miso/include/miso/constants.hpp).

As a real example from MISO, look again at the Host versions of for_each and reduce in miso/include/miso/execution.hpp. Both simply call the lambda inside a plain for statement.

/// @brief 1D parallel (concurrent) for loop in the Host backend.
template <class F> inline void for_each(backend::Host, Range1D range, F f) {
  for (int i = range.begin; i < range.end; ++i) {
    f(i);
  }
}

/// @brief 1D parallel reduction in the Host backend.
template <class T, class F, class Op>
inline T reduce(backend::Host, Range1D range, T init, F map, Op op) {
  T acc = init;
  for (int i = range.begin; i < range.end; ++i)
    acc = op(acc, map(i));
  return acc;
}

Checking part by part

If you write everything at once and then run it, it is hard to tell where the mistake is. Proceed in the following order, checking as you go.

Stage	What it checks
Initial condition	Write out the physical cell values and confirm that they equal \(\sin(\pi x_j)\)
One step	The maximum difference from the discrete solution \(G\sin(\pi x_j)\) is about the rounding error
\(n_\mathrm{steps}\) steps	The difference from the discrete solution \(G^{n}\sin(\pi x_j)\) is about \(10^{-12}\) or less, and the difference from the analytical solution agrees with the table above
Varying \(N\)	The order of convergence is close to 2

MPI version: domain decomposition and ghost cell exchange

In Practice 02, you split the serial version across \(P\) ranks. The basics of MPI (ranks, MPI_Sendrecv, MPI_Allreduce) were explained in Chapter 8. In this Practice, you do not use MISO’s mpi::Shape; you write the 1D decomposition and exchange yourself.

Grid on each rank

Divide the total number of cells \(N\) equally among \(P\) ranks, so that each rank has \(n = N/P\) physical cells (we only treat the case where \(N\) is divisible by \(P\)). The array on rank \(p\) has n + 2 elements, and the physical cell with local index i (1 to n) corresponds to the global number \(j = p\,n + (i - 1)\).

overall (N = 8, P = 2): j = 0 1 2 3 | 4 5 6 7
rank 0 array          :  [g] 0 1 2 3 [g]      g is a ghost cell
rank 1 array          :           [g] 4 5 6 7 [g]

Compute \(x_j\) for the initial condition from the global number \(j\) using the same formula as the serial version, \(x_j = (j + 1/2)\Delta x\). Also determine \(\Delta x\), \(\Delta t\), and \(n_\mathrm{steps}\) from the global \(N\), and do not change them with the number of ranks.

Exchanging ghost cells with neighboring ranks

The left neighbor of rank \(p\) is rank \(p - 1\) and its right neighbor is rank \(p + 1\). There is no neighbor to the left of the left end of the domain (rank 0) or to the right of the right end (rank \(P - 1\)). There, use MPI_PROC_NULL in place of the neighbor’s rank number. A send or receive with MPI_PROC_NULL as the partner does nothing and returns immediately, and the values in the receive location do not change. This way, the end ranks can make the same calls as the other ranks.

Before the finite difference of each step, fill the left and right ghost cells with the following two MPI_Sendrecv calls.

Call	Value sent → partner	Receive location ← partner
First (tag 0)	First physical cell u[1] → left neighbor	Right ghost cell u[n+1] ← right neighbor
Second (tag 1)	Last physical cell u[n] → right neighbor	Left ghost cell u[0] ← left neighbor

In the first call, all ranks together “send to the left and receive from the right,” so the value that rank \(p\) sends to the left goes into the right ghost cell of rank \(p - 1\). MPI_Sendrecv performs the send and the receive as a pair, so it does not happen that all ranks start only sending at the same time and wait for each other (deadlock).

After the exchange, only the ranks that hold an end of the domain fill in the physical boundary value (u[0] = -u[1] or u[n+1] = -u[n]). You can tell whether a rank holds the left end by checking whether its left neighbor is MPI_PROC_NULL. The distinction of using the value received from the neighbor at an internal boundary (between ranks) and the boundary condition value at a physical boundary is the same in MISO itself.

The address of an element of Array1D<Real> (a host array) can be taken as &u[1], so you can pass it directly to MPI_Sendrecv.

Aggregating errors and gathering results

Each rank finds the maximum error over its own physical cells with reduce. The global maximum error is the maximum of these. Obtain it by calling MPI_Allreduce with MPI_MAX.

When writing to a file, gather the physical cells of all ranks onto rank 0 with MPI_Gather. If every rank sends the same number of elements (\(n\)), they are arranged in the receive array on rank 0 in rank order, which is the order of the global number \(j\).

// args: send buffer, send count, type, receive buffer, receive count per rank, type, source rank, communicator
MPI_Gather(&u[1], n, MPI_DOUBLE, u_all.data(), n, MPI_DOUBLE, 0, comm);

For the receive location u_all, only rank 0 needs to prepare \(N\) elements (it is not used on the other ranks). Note that MPI_Abort(comm, 1), which the Practice uses when \(N\) is not divisible by \(P\), is a function that terminates all processes in that communicator.

Exercise 10.4 In the MPI version, the global maximum error matches the serial version bitwise regardless of the number of ranks. On the other hand, if each rank computes the sum of the values over all cells (\(\sum_j u_j\)) and you add them up with MPI_Allreduce (MPI_SUM), the last digit may change with the number of ranks. Explain why this difference occurs.

::

The value in each cell itself is bitwise identical regardless of the number of ranks (How to compare the results with expected values).

For the maximum, the same value is chosen in whatever order you compare, and no rounding occurs in the comparison. The result is therefore the same regardless of the order of aggregation.

For the sum, floating-point addition is not associative (the rounding of \((a + b) + c\) and \(a + (b + c)\) differs), so the result changes with the order of addition. The serial version adds from \(u_0\) in order, whereas the MPI version first forms the partial sum on each rank and then adds them, so the order of addition changes when the number of ranks changes. The difference is about the size of the rounding error. When you use the sum in a comparison across different numbers of ranks, do not require bitwise agreement; allow a difference at the rounding-error level.

GPU version: moving arrays and loops to the GPU

In Practice 03, you move the arrays and loops of the MPI version to the GPU. Host and device memory, MISO_LAMBDA, and capturing this were explained in Chapter 9. Even without a GPU, you can build the same source for the CPU and check that it gives the same result as the MPI version.

Choosing the compute location by type

As in the MISO demos, decide the type Backend, which represents the compute location, in one place.

#ifdef __CUDACC__
using Backend = backend::CUDA;
#else
using Backend = backend::Host;
#endif

__CUDACC__ is a macro that is defined only when the source is compiled with the CUDA compiler (nvcc). Write the array as Array1D<Real, Backend>, the loop as for_each(Backend{}, ...), and the reduction as reduce(Backend{}, ...). If you build with g++, Backend becomes backend::Host and the program performs the same calculation as the MPI version. MISO_LAMBDA, which you use to write lambdas, is switched in the same way in miso/include/miso/execution.hpp ([=] __host__ __device__ for CUDA and [=] for the CPU).

Reading and writing arrays on the GPU

The contents of Array1D<Real, backend::CUDA> are in GPU memory. Unlike the host array, it has no operator[], so writing u[i] in host code causes a compile error. There are two ways to read and write values on the GPU.

• Read and write through a view inside the lambda of for_each (the lambda runs on the GPU).

• Copy the whole array to or from a host Array1D<Real> with copy_from.

Array1D<T, backend::CUDA>::copy_from in miso/include/miso/array1d.hpp copies from a host array to the GPU with cudaMemcpy.

  /// @brief Copy data from another Array1D in host memory.
  void copy_from(const Array1D<T, backend::Host> &other) {
    assert(other.data() && data());
    assert(other.size() == size());
    MISO_CUDA_CHECK(cudaMemcpy(data(), other.data(), sizeof(T) * size(),
                               cudaMemcpyHostToDevice));
  }

For the opposite direction (from the GPU to the host), pass the GPU array to copy_from of the host array. The initial condition, boundary conditions, update, and error computation are all done on the GPU inside for_each and reduce. The only things copied to the host are the values used for the ghost cell exchange and the final result to be written out.

Ghost cell exchange

Whether you can pass a GPU memory address to MPI_Sendrecv depends on whether the MPI supports GPUs (CUDA-aware MPI, Chapter 9). In this Practice, so that it works with any MPI, we first copy the values to be sent and received to host arrays.

1. Pack the two values to send (the first and last physical cells) into a small array send (2 elements) on the GPU with for_each.

2. Copy send to the host array h_send with copy_from.

3. Perform the same two MPI_Sendrecv calls as in the MPI version, using h_send and h_recv (host, 2 elements).

4. Copy h_recv to the GPU array recv with copy_from.

5. With for_each, put the values of recv (the boundary condition value on a side with no neighbor) into the ghost cells.

A copy_from (cudaMemcpy) with no stream specified starts copying after the GPU work launched before it has finished, and does not return until the copy is complete. Therefore, the values packed in step 1 cannot change before they are copied in step 2.

Do not capture this

If you gather the exchange procedure into a member function of a struct, you will want to use member variables (such as the number of physical cells n) inside the lambda. However, if you use the member variable n inside a lambda, it is actually read as this->n, and this (a pointer to the struct in host memory) is captured. The GPU cannot read what that pointer points to, so the program fails at run time. The error is reported not at the lambda but at a later CUDA call (for example, the next copy_from), as follows.

CUDA Error: an illegal memory access was encountered .../miso/include/miso/array1d.hpp ...

Copy member variables into local variables before the lambda and then use those. MISO itself does the same. The following is an excerpt from Integrator<Real, backend::CUDA>::cfl in miso/include/miso/mhd_integrator_cuda.cuh.

    // device lambdas must not capture `this`
    const Real cfl_number_ = cfl_number;

Exercise 10.5 If you build and run the member function of the following struct in the GPU version, it stops with illegal memory access (an example for illustration). Identify the line that causes it and explain how to fix it.

struct Boundary {
  int n;  // number of physical cells
  void apply(Array1D<Real, Backend> &u) {
    auto uv = u.view();
    for_each(Backend{}, Range1D{0, 1}, MISO_LAMBDA(int) {
      uv[0] = -uv[1];
      uv[n + 1] = -uv[n];
    });
  }
};

::

It is the line uv[n + 1] = -uv[n];. Since n is a member variable, it is read as this->n inside the lambda, and this is captured. The lambda running on the GPU tries to read the struct in host memory and fails.

uv is a local variable (a view), so it is captured by value and is not a problem.

Put const int n_ = n; before the lambda and use n_ inside the lambda.

    const int n_ = n;
    for_each(Backend{}, Range1D{0, 1}, MISO_LAMBDA(int) {
      uv[0] = -uv[1];
      uv[n_ + 1] = -uv[n_];
    });

When built for the CPU, the memory that this points to can be read directly, so you do not notice this mistake. When you rewrite the code for the GPU, check whether any member variables appear inside lambdas.

Comparing the CPU and the GPU

Run the GPU version with the same \(N\) and number of ranks as the MPI version, and compare the values of all cells. As described in How to compare the results with expected values, the last digit may differ, so regard the results as matching if the maximum difference is about \(10^{-12}\) or less. Within the GPU version, a comparison across different numbers of ranks (one rank versus several ranks) is bitwise identical, as in the MPI version.

To build the GPU version, give CMake -DMISO_USE_CUDA=ON so that the source is compiled as CUDA. The details of the configuration are in the Practice specification. If you run into trouble with the GPU, see the reference material “GPU Troubleshooting”.

Correspondence with MISO itself

The processing you wrote in this chapter is found in MISO itself in the following places.

Processing in this chapter	MISO itself
Arrays u and u_new	mhd::Fields, which bundles Array3D arrays (miso/include/miso/mhd_fields.hpp)
Ghost cell exchange in each step	mhd::HaloExchanger::apply (miso/include/miso/mhd_halo_exchange.hpp)
MPI_PROC_NULL when there is no neighbor	Check whether x_procs_pos and the like of mpi::Shape are MPI_PROC_NULL before sending and receiving
Physical boundary values	Boundary conditions (miso/include/miso/boundary_condition.hpp and the boundary conditions of each demo)
Aggregation with reduce and MPI_Allreduce	Integrator::cfl (finds \(\Delta t\) with reduce in the CUDA version and with a for statement in the Host version, and takes the global minimum with MPI_MIN)
Exchange via the host	The GPU version of MISO passes GPU memory addresses directly to MPI (using CUDA-aware MPI)

MISO’s HaloExchanger::apply uses MPI_Isend, MPI_Irecv, and MPI_Waitall instead of MPI_Sendrecv, and sends 9 variables in 3 directions together. The flow of packing the values to send into a buffer and returning the received values to the ghost cells is the same as in the GPU version in this chapter.

In Chapter 11, we read MISO’s time evolution in order starting from main(), and relate the processing you wrote here to the classes and functions of MISO itself.