03. 1D Advection and a Periodic Boundary

The equation to solve

\[::\]

  • Time integration: explicit Euler method

  • Spatial discretization: first-order upwind difference. Since a > 0, use the difference with the left neighbor.

\[u_i^{n+1} = u_i^n - \frac{a \Delta t}{\Delta x}\left(u_i^n - u_{i-1}^n\right)\]

Choose the time step so that the Courant number \(C = a \Delta t / \Delta x\) is greater than 0 and at most 1. Here we use \(C = 0.37\).

Grid and ghost cells

Place one ghost cell on each side of the NX physical cells. The upwind difference reads only the left neighbor u[i-1], so one point on the left would suffice, but we place them symmetrically on both sides, as MISO does. With a periodic boundary, at every step, fill the ghost cells with the values of the physical cells on the opposite side, and then update.

index   0  |  1  ...  NX  |  NX+1
    ghost   physical cells  ghost

What to implement

  1. Allocate arrays u and u_new of length NX + 2, each as an Array1D<double, Backend>.

  2. Take the initial values of the physical cells to be a periodic Gaussian. Let \(L = 1\), \(x_0 = 0.5\), and \(\sigma = 0.05\), and let the coordinate of physical cell i be x = (i - 1) * L / NX. Wrap \(d = x - x_0\) by the period \(L\) into \(-L/2 \le d \le L/2\), and then compute \(\exp(-d^2 / 2\sigma^2)\). For the wrapping, use util::fmod(d, L) (the remainder of d divided by L; the result lies between -L and L). If the remainder is larger than L/2, subtract L; if it is smaller than -L/2, add L.

  3. This initial-value formula is also used in the comparison with the exact solution, so make it a function. Because it is called from inside a lambda, mark it with __host__ __device__, as in __host__ __device__ inline double gaussian(double x) (in a CPU build this specifier does nothing).

  4. Before each step, copy the last physical cell u[NX] into the left ghost cell u[0], and the first physical cell u[1] into the right ghost cell u[NX + 1]. Write this step too as a for_each over Range1D{0, 1} (a range executed only once).

  5. In a for_each over the physical-cell range Range1D{1, NX + 1}, read the pre-update u and write the updated values to u_new.

  6. At the end of the step, swap the roles of the two arrays with std::swap(u, u_new). For std::swap, include <utility>.

Array1D cannot be copied, but it can be moved (Chapter 4). std::swap uses moves to exchange the internal pointers, so no elements are copied. After the swap, the variable named u holds the memory that u_new held before. A view keeps the pointer it had when you obtained it, so after the swap, obtain the view again with u.view(). If you write the ghost-cell fill and the update as functions that take the arrays as arguments and obtain the views inside the functions, you will not forget to obtain them again.

You never write into the pre-update array midway, so each cell reads its neighbors’ values at the same time level (“do not depend on the execution order” in Chapter 7).

Checking the result

First compute just one step with NX = 256, and compare with the values obtained by a host-side for loop from the same input and ghost cells. Taking the initial amplitude 1 as the reference, check that the maximum absolute error is at most 1e-14.

Next, advance in time, and compare with the exact solution, which is the initial distribution translated by the distance \(a t\). Obtain the exact solution by replacing \(x\) in the initial-value formula with \(x - a t\) and using the same periodic wrapping. Compute the relative L2 error over the physical cells only, as follows.

\[E = \sqrt{\frac{\sum_i (u_i-u_{\mathrm{exact},i})^2}{\sum_i u_{\mathrm{exact},i}^2}}\]

You can compute the sums in the numerator and denominator with the reduce from the previous assignment (call reduce twice). Increase NX to 256, 512, and 1024, keep \(C\) and the end time the same, and examine how the error changes. For example, with 200, 400, and 800 steps, you advance to the same time. This method is first-order accurate, so on a sufficiently fine grid, halving the grid spacing roughly halves the error.

Once you have checked the one-step comparison and the time evolution on the CPU, move on to the MPI of Chapter 8. You will learn how to prepare to run the same computation on the GPU in Chapter 9.

Running on the GPU

With CUDA as well, do the one-step comparison and the time evolution with the same inputs. gaussian, which is called from inside a lambda, needs __host__ __device__. Because of FMA and similar effects, rounding may differ between the CPU and the GPU. Do not require an exact match across the board; check with the tolerance above and with how the error changes as the grid is refined.