01. Serial Version: Solving Heat Conduction on the CPU

You write the program heat1d_serial, which solves the problem in the common setup on a single-process CPU. The equations and how to check them are explained in Chapter 10.

What you use

Item

Header

miso::Array1D<Real> (a host array)

miso/array1d.hpp

miso::for_each、miso::reduce、miso::Range1D、MISO_LAMBDA、miso::backend::Host

miso/execution.hpp

miso::pi<Real>

miso/constants.hpp

miso::util::sin、util::exp、util::fabs

miso/utility.hpp

Include it with the miso/ prefix, as in #include <miso/array1d.hpp>.

What you build

Arrays

  • Add the ghost cells at both ends to the number of physical cells \(N\), and allocate Array1D<Real> u(N + 2) and the destination of the update, u_new(N + 2).

  • Physical cell \(j\) (\(j = 0, \dots, N-1\)) is stored at array index j + 1. The range of the physical cells is Range1D{1, N + 1}.

Steps

  1. Initial condition: With for_each, put \(\sin(\pi x_j)\) into the physical cells.

  2. Time integration: Repeat the following \(n_\mathrm{steps}\) times.

    1. Put the boundary-condition values into the ghost cells: u[0] = -u[1], u[N+1] = -u[N]. So that you can move this to the GPU later, write it inside for_each(backend::Host{}, Range1D{0, 1}, ...).

    2. With for_each, compute the following expression for the physical cells. Read u and write to u_new.

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

      Write 2 as Real(2). Here \(r = \kappa\Delta t/\Delta x^2\).

    3. Swap the arrays with std::swap(u, u_new);.

  3. Error: With reduce, compute the maximum difference from the discrete solution and the maximum difference from the analytical solution.

  4. Output: Print the 3 lines described in the common setup, and write the result to result_N<N>_serial.txt (result_N64_serial.txt for \(N = 64\)).

Points to watch

  • Inside lambdas, use a view (auto uv = u.view();), not the Array1D.

  • After std::swap, the view points to the data from before the swap. Obtain the view again every time inside the step loop.

  • Write the maximum operation of reduce as MISO_LAMBDA(Real a, Real b) { return a > b ? a : b; }. The initial value is Real(0) (the absolute value of a difference is at least 0, so Real(0) serves in place of the identity element of the maximum).

  • Compute the coefficient of the discrete solution, \(G^{n_\mathrm{steps}}\), only once on the host as std::pow(G, n_steps) (<cmath>), and capture it in the lambda.

What to check

  1. Initial condition: Run the time integration for 0 steps (or write the output right after the initial condition), and check that \(u_j = \sin(\pi x_j)\).

  2. One step: Advance just one step instead of \(n_\mathrm{steps}\), and check that the maximum difference from the discrete solution \(G\sin(\pi x_j)\) is about \(10^{-15}\).

  3. Default input (\(N = 64\)): Check that you get the following values.

    N = 64, steps = 1024, r = 0.4000
max error (discrete) = 1.249e-14
max error (exact)    = 1.034113e-04

    max error (discrete) is a rounding error, so it changes with the compiler and the way you write the code. It is fine as long as it is well below \(10^{-12}\). max error (exact) should agree to the 4th or 5th digit.

  4. Order of convergence: Run with \(N = 32, 64, 128, 256, 512\), tabulate max error (exact), and compute the ratios of neighboring values and the order of convergence \(p = \log_2(E_N/E_{2N})\).

    \(N\)

    steps

    max error (exact)

    32

    256

    4.136835e-04

    64

    1024

    1.034113e-04

    128

    4096

    2.585224e-05

    256

    16384

    6.463023e-06

    512

    65536

    1.615753e-06

    Every \(p\) is close to 2.000. max error (discrete) is at most \(10^{-11}\) for every \(N\).

  5. Stability condition (optional): Run ./heat1d_serial 64 0.6 and check that the error becomes very large (the computation breaks down). Also check that ./heat1d_serial 64 0.5 does not break down.

If the difference from the discrete solution is large, check in order the ghost-cell values, the range of the physical cells (the upper limit of Range1D is excluded), and the expression for \(x_j\) (the \(+1/2\) and the index offset). If it agrees with the discrete solution but the difference from the analytical solution differs from the table, check \(\Delta t\), \(n_\mathrm{steps}\), and \(t_\mathrm{end}\).

Criteria for moving on

Once you obtain the values in 3 and 4 above, move on to the 02. MPI version. Keep the files you wrote, such as result_N64_serial.txt, because you will use them in the comparison with the MPI version.