11. Reading MISO’s Time Evolution

The heat conduction calculation you wrote in Chapter 10 had these parts: determining the time step, the finite difference of the right-hand side, time integration, boundary handling, and output. MISO’s MHD calculation is built from the same parts. In this chapter, starting from the demo’s main(), we trace the contents of one step and the output through the MISO source, and relate the equations and finite-difference formulas to the code. Finally, we identify the places you rewrite when you create your own problem.

Follow the paths in the text from the folder where you cloned MISO. Even if you have not finished the GPU chapter (Chapter 9), you can read on with the CPU source alone. For the CUDA version, we show only how it corresponds to the CPU version.

Contents

1. Basic equations and MISO’s variables

1.1 State variables

MISO holds the following 9 variables at each grid point. Fields in miso/include/miso/mhd_fields.hpp bundles them as members of type Array3D (Chapter 4).

Member name

Symbol

Meaning

ro

\(\rho\)

Density

vx, vy, vz

\(v_x, v_y, v_z\)

Velocity

bx, by, bz

\(B_x, B_y, B_z\)

Magnetic field

ei

\(e\)

Internal energy per unit mass

ph

\(\phi\)

Auxiliary variable for reducing the divergence of the magnetic field (the potential of divB cleaning)

The conserved quantities are the momentum \(\rho\boldsymbol{v}\), the total energy, and so on, but the arrays hold the velocity \(\boldsymbol{v}\) and the internal energy per unit mass \(e\). In the time integration, we first form the conserved quantities and advance them, and at the end divide by \(\rho\) to recover the velocity and \(e\) (Section 8).

1.2 Basic equations

MISO solves the ideal MHD equations in Gaussian units (the form in which the magnetic terms carry a factor \(1/4\pi\)). \(\partial_j\) is the partial derivative in the \(x, y, z\) directions, and when the same index \(j\) appears twice, it is summed over \(j = x, y, z\). \(S_\rho\) and the like are source terms (external forces and heating), which the user provides (Section 12).

\[\begin{split} \begin{aligned} \frac{\partial \rho}{\partial t} &= -\partial_j(\rho v_j) + S_\rho, \\ \frac{\partial (\rho v_i)}{\partial t} &= -\partial_j(\rho v_i v_j) - \partial_i p - \partial_i\left(\frac{B^2}{8\pi}\right) + \frac{1}{4\pi}\partial_j(B_i B_j) + S_{v_i}, \\ \frac{\partial B_i}{\partial t} &= \partial_j(v_i B_j - v_j B_i) - \partial_i \phi + S_{B_i}, \\ \frac{\partial \phi}{\partial t} &= -c_h^2\, \partial_j B_j - \frac{\phi}{\tau}, \\ \frac{\partial E}{\partial t} &= -\partial_j\left[\left(\rho e + p + \frac{1}{2}\rho v^2 + \frac{B^2}{4\pi}\right) v_j - \frac{(\boldsymbol{v}\cdot\boldsymbol{B})}{4\pi} B_j\right] + S_E, \\ E &= \rho e + \frac{1}{2}\rho v^2 + \frac{B^2}{8\pi}. \end{aligned} \end{split}\]

The term \(-\partial_i\phi\) in the magnetic field equation and the equation for \(\phi\) form a mechanism that carries the error of \(\nabla\cdot\boldsymbol{B}\) generated in the numerical calculation away as a wave (with speed \(c_h\)) and damps it over a time \(\tau\). This method is called hyperbolic divergence cleaning (Dedner et al. 2002, J. Comput. Phys. 175, 645).

1.3 Equation of state

The pressure \(p\) and the sound speed \(c_s\) are obtained from the equation of state (EOS). The ideal-gas EOS that MISO provides is eos::IdealEOS in miso/include/miso/eos.hpp. With the ratio of specific heats \(\gamma\) (eos.gm in the configuration),

\[ p = (\gamma - 1)\rho e, \qquad c_s = \sqrt{\gamma(\gamma - 1) e} = \sqrt{\gamma p/\rho} \]

are computed by the two member functions gas_pressure and sound_speed. Both process all grid points with for_each (Chapter 7).

// miso/include/miso/eos.hpp: IdealEOS::gas_pressure (excerpt)
template <typename Backend>
void gas_pressure(Backend btag, mhd::FieldsView<const Real> qq,
                  Array3DView<Real> pr) const {
  const Real gm_ = gm;  // device lambdas must not capture `this`
  Range1D range{0, qq.size()};
  for_each(
      btag, range,
      MISO_LAMBDA(int i) { pr[i] = (gm_ - Real(1)) * qq.ro[i] * qq.ei[i]; });
}

The member variable gm is copied to the local variable gm_ before being used in the lambda in order to avoid the this capture explained in Chapter 9.

2. main() and Model

The main() in demo/mhd2d_vortex/src/main.cpp consists of just the following 4 lines.

// 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();
}

Line

Role

In Fortran

Env env(argc, argv);

Initializes MPI. At the end of main(), when the lifetime of env ends, the destructor calls MPI_Finalize (Chapter 8)

call MPI_Init(ierr) and the final call MPI_Finalize(ierr)

parse_config_filepath(argc, argv)

Looks for --config=<path> or --config <path> on the command line. If it is absent, returns an empty std::optional

Code that reads arguments with get_command_argument

Config config(...)

Reads the configuration file and fills in default values (Section 3). value_or("./config.yaml") means to use ./config.yaml if none is specified

Code that reads a namelist

Model model(config);

Creates the grid, arrays, and time-integration components together (Section 4)

The calls to each *_init and allocate

model.run();

Runs the time loop (Section 4)

The do while (time < tend) time loop

Model is a struct that you write for each demo.

// demo/mhd2d_vortex/src/main.cpp
struct Model : public mhd::ModelBase<Model, Real, Backend> {
  eos::IdealEOS<Real> eos;
  InitialCondition ic;
  mhd::EmptyBoundaryCondition<Real> bc;
  mhd::EmptySourceTerm<Real> src;

  Model(Config &config) : ModelBase(config), eos(config), ic(eos), bc(), src() {}
};

Model inherits from mhd::ModelBase and passes itself (Model) as the template argument. This is the CRTP you learned about in Chapter 5. The derived class holds the 4 components that change from problem to problem.

Member

Role

Type in this example

eos

Equation of state

eos::IdealEOS<Real> (ideal gas)

ic

Initial condition

InitialCondition (defined in the same file)

bc

Boundary condition

mhd::EmptyBoundaryCondition<Real> (does nothing; used when all directions are periodic boundaries)

src

Source term

mhd::EmptySourceTerm<Real> (no source term)

Real and Backend are decided at the top of the file by using Real = float; and #ifdef USE_CUDA. Real = float gives single precision, and the calculation runs on backend::CUDA (GPU) if USE_CUDA is defined and on backend::Host (CPU) otherwise (Chapter 2, Chapter 6).

3. Configuration file and default values

3.1 Sections of config.yaml

demo/mhd2d_vortex/config.yaml has the following sections.

Section

Contents

Where it is read

io

Whether to write output, the output destination, and whether to restart (continue)

config.hpp, time.hpp, mhd_checkpoint.hpp

time

End time tend and output interval dt_output

Time in time.hpp

grid

Number of grid points i_size and so on (the global number of grid points), ghost cell width margin, domain x_min and so on

Grid in grid.hpp

mpi

Number of processes in each direction, x_procs and so on

mpi::Shape in mpi_util.hpp

domain

Directions made periodic, periodic

mpi::Shape in mpi_util.hpp

mhd

CFL number cfl_number and artificial viscosity parameters

mhd_integrator_host.hpp, mhd_artificial_viscosity_host.hpp

eos

Ratio of specific heats gm

IdealEOS in eos.hpp

To read a value, specify the section name and the key in turn with [], as in config["time"]["tend"].as<Real>(), and extract it by choosing the type with as<Type>(). The <> in as<Real>() is a template argument (Chapter 6).

3.2 Default values and merge_yaml

So that you do not have to write every key each time, MISO keeps default values as a YAML string in miso/include/miso/config_defaults.hpp.

// miso/include/miso/config_defaults.hpp(excerpt; some comments omitted)
inline constexpr const char *config_default_yaml = R"yaml(
io:
  enabled: true
  continue: true
  save_dir: data/
  ...
mhd:
    # CFL number for time-stepping
    cfl_number: 0.5
    artificial_viscosity:
        ep: 1.0
        fh: 1.0
        ...
    )yaml";

R"yaml( ... )yaml" is a raw string literal. It turns everything from R"delimiter( to )delimiter", including newlines and ", into a single string as it is. It is used to embed YAML directly in the C++ source.

The constructor of Config in miso/include/miso/config.hpp does the following in order.

  1. Only rank 0 reads the file (YAML::LoadFile), converts it to a string, and distributes it to all ranks with MPI_Bcast.

  2. All ranks read the default-value string (YAML::Load(config_default_yaml)) and add it to the user’s settings with merge_yaml.

  3. Interpret io.save_dir as a path relative to the folder containing the configuration file, and create the output destination.

// miso/include/miso/config.hpp
inline void merge_yaml(YAML::Node &target, const YAML::Node &defaults) {
  if (!defaults.IsMap()) {
    return;
  }

  if (!target || !target.IsMap()) {
    target = YAML::Node(YAML::NodeType::Map);
  }

  for (const auto &it : defaults) {
    const std::string key = it.first.as<std::string>();
    const YAML::Node &default_value = it.second;

    if (default_value.IsMap()) {
      YAML::Node child = target[key];
      merge_yaml(child, default_value);  // fill nested sections recursively
      target[key] = child;
    } else {
      if (!target[key]) {               // fill only keys the user did not write
        target[key] = default_value;
      }
    }
  }
}

(The two comments in merge_yaml were added for explanation.) Values written by the user are not overwritten, and only the keys that were not written are filled with default values. Sections inside sections (such as mhd.artificial_viscosity) are also filled in key by key. The completed result is saved as config.yaml in the output destination (Section 11). You can check the values actually used in this saved file.

4. Constructing ModelBase and run()

4.1 Construction order of the members

ModelBase in miso/include/miso/mhd_model_base.hpp has the following members.

// miso/include/miso/mhd_model_base.hpp(excerpt)
template <class Derived, class Real, class Backend> class ModelBase {
public:
  Config &config;
  mpi::Shape mpi_shape;
  Time<Real> time;
  Grid<Real, backend::Host> grid;

  ExecContext<Real, Backend> exec_ctx;
  MHD<Real, Backend> mhd;

  explicit ModelBase(Config &cfg)
      : config(cfg), mpi_shape(cfg), time(cfg), grid(cfg, mpi_shape),
        exec_ctx(mpi_shape, grid), mhd(cfg, grid, exec_ctx) {}
  ...

Members are constructed in the order of declaration (Chapter 4). Later members use members that were constructed earlier.

Order

Member

What is created

1

config

Reference to the configuration

2

mpi_shape

The Cartesian process layout created by MPI_Cart_create, and the neighboring ranks in each direction (Chapter 8)

3

time

Time, step count, and output number

4

grid

The grid this rank is responsible for (coordinates including ghost cells, and grid spacing). Created on the CPU

5

exec_ctx

Extra information for each compute location (on the CPU, only a reference to mpi_shape; on the GPU, it also holds the kernel launch configuration)

6

mhd

The state variables qq, the time-integration Integrator, and the output Checkpoint

The type MHD of mhd is in miso/include/miso/mhd_solver.hpp.

// miso/include/miso/mhd_solver.hpp(excerpt)
template <typename Real, typename Backend> struct MHD {
  Grid<Real, Backend> grid;
  Fields<Real, Backend> qq;
  ExecContext<Real, Backend> &exec_ctx;
  Integrator<Real, Backend> integrator;
  Checkpoint<Real> checkpoint;
  ...
  template <typename EOS, typename BoundaryCondition, typename Source>
  void update(Real dt, const EOS &eos, const BoundaryCondition &bc,
              const Source &src) {
    integrator.update(dt, eos, bc, src, qq);
  }
  ...
};

MHD holds the state variables qq and the computation components, and acts as the front end that hands the work to Integrator. Grid<Real, Backend> is a copy of the CPU grid grid_h placed on the compute location (Backend).

The members of the derived class Model (eos, ic, bc, src) are constructed after the construction of the base class ModelBase is finished. For this reason, the constructor of ModelBase does not use them, and run() receives them through derived().

4.2 run(): the time loop

// miso/include/miso/mhd_model_base.hpp(excerpt)
void run() {
  // Ensure to call methods in the derived class
  auto &d = derived();

  // Assert that d has eos, ic, bc, and src as member variables.
  static_assert(has_eos<Derived>::value,
                "Derived must have member variable 'eos'.");
  ...
  mhd.apply_initial_condition(d.ic, d.bc);
  if (config["io"]["continue"].template as<bool>() &&
      fs::exists(time.n_output_filepath())) {
    load_state();
  }

  save_metadata();
  save_if_needed();

  while (time.time < time.tend) {
    d.update();
    save_if_needed();
  }
}

Processing

Contents

derived()

Treats itself as the derived class Model with static_cast<Derived &>(*this) (Chapter 5)

static_assert(has_eos<Derived>::value, ...)

Causes a compile-time error if Model lacks members such as eos (Chapter 6)

mhd.apply_initial_condition(d.ic, d.bc)

Fills a temporary array on the CPU with initial values via ic.apply, copies it to qq on the compute location, and performs the boundary condition and halo exchange once

load_state()

If io.continue is true and the previous output remains, restarts from the latest output (Section 11)

save_metadata()

Writes out the configuration (with defaults filled in), the grid, and the MPI layout

save_if_needed()

Writes out the state if the output time has been reached

d.update()

Advances one step

The .template in config["io"]["continue"].template as<bool>() is the syntax required inside a template when calling a member template of a dependent type. You can read it as meaning the same as as<bool>().

d.update() calls update() if Model defines it, and otherwise calls the base class ModelBase::update(). The Model in demo/mhd3d_magnetosphere/src/main.cpp defines its own update() and changes how the time step is determined.

4.3 update(): the entry point of one step

// miso/include/miso/mhd_model_base.hpp(excerpt)
void update() {
  auto &d = derived();
  using Source = decltype(d.src);
  ...
  auto dt = mhd.cfl(d.eos);
  if constexpr (has_dt_limit<Source>::value) {
    // dt_limit() returns the limit of this rank: take the global minimum.
    ...
    Real dt_src = static_cast<Real>(d.src.dt_limit());
    ...
    MPI_Allreduce(&dt_src, &dt_src_g, 1, mpi::data_type<Real>(), MPI_MIN,
                  mpi::comm());
    ...
    dt = std::min(dt, dt_src_g);
  }
  mhd.update(dt, d.eos, d.bc, d.src);
  time.update(dt);
}

  1. Determine the time step from the CFL condition with mhd.cfl(d.eos) (Section 6).

  2. If the source term has dt_limit(), also apply that limit (the minimum over all ranks). If it does not, the inside of the if constexpr is not compiled.

  3. Advance the state one step with mhd.update (Section 5).

  4. Advance the time and step count with time.update(dt).

5. Flow of one step

Integrator::update, called from MHD::update, is the body of one step. Integrator is declared in miso/include/miso/mhd_integrator.hpp, and the CPU and GPU versions are specialized in separate files (Chapter 6).

File

Contents

miso/include/miso/mhd_integrator.hpp

The declaration of Integrator, source_ro and the like that extract source terms, and the fourth-order central difference space_centered_4th

miso/include/miso/mhd_integrator_host.hpp

Integrator<Real, backend::Host> (CPU version)

miso/include/miso/mhd_integrator_cuda.cuh

Integrator<Real, backend::CUDA> (GPU version; included only when compiling with nvcc)

The CPU version of update is 3 lines.

// miso/include/miso/mhd_integrator_host.hpp
template <typename EOS, typename BoundaryCondition, typename Source>
void update(Real dt, const EOS &eos, const BoundaryCondition &bc,
            const Source &src, Fields<Real, backend::Host> &qq) {
  divb_parameters_set(dt);
  runge_kutta_4step(dt, eos, bc, src, qq);
  apply_artificial_viscosity(dt, eos, bc, qq);
}

The whole step, arranged in the order of function calls, is as follows.

ModelBase::update
├─ MHD::cfl → Integrator::cfl                      time step dt (section 6)
├─ MHD::update → Integrator::update
│   ├─ divb_parameters_set(dt)                     set c_h and τ
│   ├─ runge_kutta_4step                           4-stage time integration (section 7)
│   │   ├─ update_sc4(dt/4) → boundary conditions and halo exchange
│   │   ├─ update_sc4(dt/3) → boundary conditions and halo exchange
│   │   ├─ update_sc4(dt/2) → boundary conditions and halo exchange
│   │   └─ update_sc4(dt)   → boundary conditions and halo exchange
│   └─ apply_artificial_viscosity                   artificial viscosity (section 10)
│       ├─ x direction → boundary conditions and halo exchange
│       ├─ y direction → boundary conditions and halo exchange
│       └─ z direction → boundary conditions and halo exchange
└─ Time::update                                    advance the time

divb_parameters_set determines the divergence-cleaning speed \(c_h\) and the damping time \(\tau\) to match the time step.

// miso/include/miso/mhd_integrator_host.hpp
void divb_parameters_set(Real dt) {
  ch_divb = 0.8 * cfl_number * grid.min_dxyz / dt;
  ch_divb_square = ch_divb * ch_divb;
  tau_divb = 2.0 * dt;
}

Since \(c_h \Delta t = 0.8\, C \min(\Delta x, \Delta y, \Delta z)\) (\(C\) is the CFL number), the waves of \(\phi\) also propagate at a speed that satisfies the CFL condition.

The CPU and GPU versions have the same function names and call order. The difference is whether the per-grid-point computation is written as a for statement or as a CUDA kernel.

Processing

CPU version (mhd_integrator_host.hpp)

GPU version (mhd_integrator_cuda.cuh)

Preparing \(B^2\), \(h_t\), and \(\boldsymbol{v}\cdot\boldsymbol{B}\)

The first triple loop of update_sc4

bb_ht_vb_kernel

Update of each variable

The second triple loop of update_sc4

update_ro_kernel, update_vx_kernel to update_ei_kernel

CFL condition

Take the minimum in a triple loop, then MPI_Allreduce

reduce(backend::CUDA{}, ...) and MPI_Allreduce

runge_kutta_4step, apply_boundary_condition

Same order

Same order

Artificial viscosity

impl_host::ArtificialViscosity (mhd_artificial_viscosity_host.hpp)

impl_cuda::ArtificialViscosity (mhd_artificial_viscosity_cuda.cuh)

The CPU version of update_sc4 is written with an ordinary triple for statement, not for_each. With i as the outer loop and k as the inner loop, it reads memory along the layout of Array3D, in which the last index is contiguous (Chapter 7).

6. Time step from the CFL condition

The time step is chosen so that, at every grid point, the fastest wave does not cross the grid spacing in one step.

\[ \Delta t = C \min_{i,j,k} \frac{\min(\Delta x_i, \Delta y_j, \Delta z_k)}{c_s + |\boldsymbol{v}| + c_A}, \qquad c_A = \frac{|\boldsymbol{B}|}{\sqrt{4\pi\rho}} \]

\(C\) is mhd.cfl_number in the configuration, and \(c_A\) is the Alfvén speed.

// miso/include/miso/mhd_integrator_host.hpp: Integrator::cfl (excerpt)
Real dt = 1.e10;
eos.sound_speed(backend::Host{}, qq.const_view(), cs.view());
for (int i = grid.i_margin; i < grid.i_total - grid.i_margin; ++i) {
  for (int j = grid.j_margin; j < grid.j_total - grid.j_margin; ++j) {
    for (int k = grid.k_margin; k < grid.k_total - grid.k_margin; ++k) {
      ...
      Real total_vel = cs(i, j, k) + vv + ca;
      dt = util::min2(dt, cfl_number
        * util::min3<Real>(grid.dx[i], grid.dy[j], grid.dz[k])/total_vel);
    }
  }
}
Real dt_global;
MPI_Allreduce(&dt, &dt_global, 1, mpi::data_type<Real>(), MPI_MIN,
              mpi::comm());
return dt_global;

The loop runs from i_margin up to just before i_total - i_margin, that is, over the grid points this rank is responsible for, excluding ghost cells. After each rank finds its minimum, MPI_MIN in MPI_Allreduce makes the global minimum common to all ranks (Chapter 8). This is because, unless all ranks advance with the same \(\Delta t\), the times of the values exchanged with neighbors would be inconsistent.

7. Four-stage Runge–Kutta method

7.1 runge_kutta_4step

Write the equation as \(\partial q/\partial t = R(q)\). \(q\) is the conserved quantities (\(\rho\), \(\rho\boldsymbol{v}\), \(\boldsymbol{B}\), \(\phi\), \(E\)), and \(R\) is the right-hand side (spatial derivatives and source terms). MISO builds \(q^{n+1}\) from \(q^n\) in the following 4 stages.

\[\begin{split} \begin{aligned} q^{(1)} &= q^n + \frac{\Delta t}{4} R(q^n), \\ q^{(2)} &= q^n + \frac{\Delta t}{3} R(q^{(1)}), \\ q^{(3)} &= q^n + \frac{\Delta t}{2} R(q^{(2)}), \\ q^{n+1} &= q^n + \Delta t\, R(q^{(3)}). \end{aligned} \end{split}\]

Each stage has the form “add the right-hand side, evaluated with the value from the previous stage, to the original value \(q^n\).” The only arrays needed are three sets: the original value, the value of the previous stage, and the result.

// miso/include/miso/mhd_integrator_host.hpp
template <typename EOS, typename BoundaryCondition, typename Source>
void runge_kutta_4step(const Real dt, const EOS &eos,
                       const BoundaryCondition &bc, const Source &src,
                       Fields<Real, backend::Host> &qq) {
  update_sc4(dt / 4.0, eos, src, qq, qq, qq_rslt);
  qq_argm.copy_from(qq_rslt);
  apply_boundary_condition(bc, qq_argm);

  update_sc4(dt / 3.0, eos, src, qq, qq_argm, qq_rslt);
  qq_argm.copy_from(qq_rslt);
  apply_boundary_condition(bc, qq_argm);

  update_sc4(dt / 2.0, eos, src, qq, qq_argm, qq_rslt);
  qq_argm.copy_from(qq_rslt);
  apply_boundary_condition(bc, qq_argm);

  update_sc4(dt, eos, src, qq, qq_argm, qq_rslt);
  qq.copy_from(qq_rslt);
  apply_boundary_condition(bc, qq);
}

The arguments of update_sc4(stage dt, eos, src, original value, value for evaluating the right-hand side, result) correspond to the formulas as follows.

Argument

Array

Symbol in the formula

qq_orgn (the 3rd Fields)

qq

\(q^n\) (the value at the start of the step; unchanged during the 4 stages)

qq_argm (the 4th)

qq for the 1st stage, qq_argm for the 2nd stage onward

The value from the previous stage that is fed into \(R(\cdot)\)

qq_rslt (the 5th)

qq_rslt

The result of that stage

The result of each stage is written to qq_rslt, copied to qq_argm with copy_from, and then the boundary condition is applied. update_sc4 updates only the grid points excluding the ghost cells, so the ghost cells must be refilled before the right-hand side of the next stage is computed. This is the same order as in the heat conduction calculation of Chapter 10, where the boundary values were prepared before each stage of the time integration. Only in the last stage is the result put back into qq, and this becomes \(q^{n+1}\).

7.2 Checking the stage coefficients

For the linear equation \(dq/dt = \lambda q\), we have \(R(q) = \lambda q\), so substituting in order with \(h = \lambda\Delta t\) gives

\[ q^{n+1} = \left(1 + h + \frac{h^2}{2} + \frac{h^3}{6} + \frac{h^4}{24}\right) q^n \]

which agrees with the Taylor expansion of the exact solution \(e^{h} q^n\) up to the fourth-order term. The following short program solves \(dq/dt = -q\) up to \(t = 1\) with the same 4 stages and prints the difference from the exact solution \(e^{-1}\).

// rk4_linear.cpp: solve dq/dt = -q with the same 4-stage coefficients as MISO (illustrative example)
#include <cmath>
#include <iostream>

double solve(int n_step) {
  const double dt = 1.0 / n_step;
  double q = 1.0;                       // value at time t (qq in MISO)
  for (int n = 0; n < n_step; ++n) {
    double q_argm = q;                  // value used to evaluate the right-hand side (qq_argm in MISO)
    q_argm = q + dt / 4.0 * (-q_argm);  // stage 1
    q_argm = q + dt / 3.0 * (-q_argm);  // stage 2
    q_argm = q + dt / 2.0 * (-q_argm);  // stage 3
    q = q + dt * (-q_argm);             // stage 4 updates q
  }
  return q;
}

int main() {
  for (int n : {10, 20, 40}) {
    std::cout << "n_step = " << n
              << "  error = " << std::abs(solve(n) - std::exp(-1.0)) << "\n";
  }
}

g++ -std=c++17 -O2 -o rk4_linear rk4_linear.cpp
./rk4_linear

n_step = 10  error = 3.33241e-07
n_step = 20  error = 1.99761e-08
n_step = 40  error = 1.22274e-09

Doubling the number of steps reduces the error to about \(1/16 = 1/2^4\), which shows that the method is fourth-order accurate.

::

This form of multistage method is known from its use by Jameson, Schmidt, and Turkel (1981, AIAA Paper 81-1259) for compressible flow calculations. Its advantage is that few arrays need to be kept at each stage.

The coefficients \(1/4, 1/3, 1/2, 1\) are chosen so that the method is fourth-order accurate for linear equations. For a nonlinear equation \(dq/dt = f(q)\), expanding and comparing the coefficient of \(f''f^2\) in the \(\Delta t^3\) term gives \(1/8\) for this method and \(1/6\) for the exact solution, which do not agree. The formal accuracy for nonlinear equations is therefore second order. On the other hand, the linear stability region is the same as that of the classical four-stage Runge–Kutta method. The “stability condition of RK4, \(\lambda\Delta t < 2.78\)” in the comment of miso/include/miso/mhd_resistivity.hpp also comes from this stability region.

8. Fourth-order central difference and the right-hand side

8.1 space_centered_4th

Spatial derivatives are approximated with a fourth-order central difference. On a uniformly spaced grid,

\[ \left.\frac{\partial f}{\partial x}\right|_i \approx \frac{-f_{i+2} + 8 f_{i+1} - 8 f_{i-1} + f_{i-2}}{12\, \Delta x} \]

The error is proportional to \(\Delta x^4\). Because it uses 2 points on each side, the ghost cell width (grid.margin) must be at least 2.

// miso/include/miso/mhd_integrator.hpp
template <typename Real>
__host__ __device__ inline Real
space_centered_4th(Array3DView<const Real> qq, Real dxyzi, int i, int j, int k,
                   int is, int js, int ks) {
  // clang-format off
  return (
    -     qq(i + 2*is, j + 2*js, k + 2*ks)
    + 8.0*qq(i +   is, j +   js, k +   ks)
    - 8.0*qq(i -   is, j -   js, k -   ks)
    +     qq(i - 2*is, j - 2*js, k - 2*ks)
  )*inv12<Real>*dxyzi;
  // clang-format on
};

Argument

Meaning

qq

The array to differentiate (a const view)

dxyzi

The reciprocal of the grid spacing, \(1/\Delta x\) (grid.dxi[i] and so on)

i, j, k

The grid point at which the derivative is evaluated

is, js, ks

The direction of differentiation. For the \(x\) direction, (grid.is, 0, 0)

Because the direction is selected with is, js, ks, one function can handle the derivatives in all 3 directions. grid.is is 1 if the number of grid points in the \(x\) direction is 2 or more, and 0 if it is 1. In a direction with one grid point, the index does not shift, and the derivative in that direction becomes 0 (Exercise 4). inv12<Real> is a variable template that represents \(1/12\) in type Real (miso/include/miso/constants.hpp, Chapter 6).

There are also functions with the same name that differentiate the product of two or three arrays (overloads, Chapter 2). space_centered_4th(c_ro, c_vx, dxi, ...) computes \(\partial_x(\rho v_x)\), and space_centered_4th(c_ro, c_vx, c_vy, dyi, ...) computes \(\partial_y(\rho v_x v_y)\). The product is formed at each grid point before the difference is taken.

8.2 Preparation in update_sc4

update_sc4 first prepares in arrays the quantities that are used many times on the right-hand side.

// miso/include/miso/mhd_integrator_host.hpp: update_sc4 (excerpt; formulas omitted, explanatory comments added)
eos.gas_pressure(backend::Host{}, qq_argm.const_view(), pr.view());

for (int i = 0; i < grid.i_total; ++i) {
  for (int j = 0; j < grid.j_total; ++j) {
    for (int k = 0; k < grid.k_total; ++k) {
      bb(i, j, k) = ...;  // B^2
      ht(i, j, k) = ...;  // ρe + p + B^2/4π + ρv^2/2
      vb(i, j, k) = ...;  // v·B
    }
  }
}

Array

Formula

Where it is used

pr

\(p = (\gamma-1)\rho e\)

Pressure gradient in the equation of motion

bb

\(B^2\)

Magnetic pressure \(B^2/8\pi\) in the equation of motion

ht

\(\rho e + p + \frac{1}{2}\rho v^2 + \frac{B^2}{4\pi}\)

Energy flux

vb

\(\boldsymbol{v}\cdot\boldsymbol{B}\)

Energy flux

Because differences are taken, these are computed at all grid points including ghost cells. All of them are made from qq_argm (the value used to evaluate the right-hand side).

Next, read-only views are created with const_view() and given short names (c_ro, c_vx, … c_ph).

8.3 Correspondence between the equations and the code

The second triple loop updates the grid points excluding the ghost cells. We write the time step of a stage as \(\Delta t'\) (the dt argument; for example \(\Delta t/4\)) and space_centered_4th as \(D_x, D_y, D_z\).

Continuity equation

// miso/include/miso/mhd_integrator_host.hpp: update_sc4 (excerpt)
qq_rslt.ro(i, j, k) = qq_orgn.ro(i, j, k) + dt * (
-space_centered_4th(c_ro, c_vx, dxi, i, j, k, grid.is, 0, 0)
-space_centered_4th(c_ro, c_vy, dyi, i, j, k, 0, grid.js, 0)
-space_centered_4th(c_ro, c_vz, dzi, i, j, k, 0, 0, grid.ks)
+source_ro(src, c_qq, i, j, k)
);
\[ \rho^{\text{rslt}} = \rho^{n} + \Delta t' \left[ -D_x(\rho v_x) - D_y(\rho v_y) - D_z(\rho v_z) + S_\rho \right] \]

The \(\rho v_x\) and so on on the right-hand side are values from qq_argm, and the \(\rho^n\) added to it is the value from qq_orgn. This matches the form of the formula in Section 7.

Equation of motion (\(x\) component)

// miso/include/miso/mhd_integrator_host.hpp: update_sc4 (excerpt)
qq_rslt.vx(i, j, k) = (
    qq_orgn.ro(i, j, k) * qq_orgn.vx(i, j, k) + dt * (
    -space_centered_4th(c_ro, c_vx, c_vx, dxi, i, j, k, grid.is, 0, 0)
    -space_centered_4th(c_ro, c_vx, c_vy, dyi, i, j, k, 0, grid.js, 0)
    -space_centered_4th(c_ro, c_vx, c_vz, dzi, i, j, k, 0, 0, grid.ks)
    -space_centered_4th(c_pr, dxi, i, j, k, grid.is, 0, 0)
    -space_centered_4th(c_bb, dxi, i, j, k, grid.is, 0, 0)*pii8<Real>
    +pii4<Real>*(
        +space_centered_4th(c_bx, c_bx, dxi, i, j, k, grid.is, 0, 0)
        +space_centered_4th(c_bx, c_by, dyi, i, j, k, 0, grid.js, 0)
        +space_centered_4th(c_bx, c_bz, dzi, i, j, k, 0, 0, grid.ks)
        )
    +source_vx(src, c_qq, i, j, k)
    )
)/qq_rslt.ro(i, j, k);
\[ v_x^{\text{rslt}} = \frac{1}{\rho^{\text{rslt}}}\left\{ \rho^n v_x^n + \Delta t' \left[ -D_j(\rho v_x v_j) - D_x p - \frac{1}{8\pi} D_x(B^2) + \frac{1}{4\pi} D_j(B_x B_j) + S_{v_x} \right]\right\} \]

The conserved quantity \(\rho v_x\) is advanced and then divided by the new density \(\rho^{\text{rslt}}\) to recover the velocity. For this reason, the continuity equation (qq_rslt.ro) must be computed first. pii4<Real> and pii8<Real> are \(1/4\pi\) and \(1/8\pi\) (miso/include/miso/constants.hpp). The \(y, z\) components have the same form.

Induction equation (\(x\) component)

// miso/include/miso/mhd_integrator_host.hpp: update_sc4 (excerpt)
qq_rslt.bx(i, j, k) = qq_orgn.bx(i, j, k) + dt * (
    -space_centered_4th(c_vy, c_bx, dyi, i, j, k, 0, grid.js, 0)
    -space_centered_4th(c_vz, c_bx, dzi, i, j, k, 0, 0, grid.ks)
    +space_centered_4th(c_vx, c_by, dyi, i, j, k, 0, grid.js, 0)
    +space_centered_4th(c_vx, c_bz, dzi, i, j, k, 0, 0, grid.ks)
    -space_centered_4th(c_ph, dxi, i, j, k, grid.is, 0, 0)
    +source_bx(src, c_qq, i, j, k)
);
\[ B_x^{\text{rslt}} = B_x^n + \Delta t' \left[ D_y(v_x B_y - v_y B_x) + D_z(v_x B_z - v_z B_x) - D_x\phi + S_{B_x} \right] \]

This is the induction equation of Section 1.2 with \(i = x\) (the \(j = x\) term is not written because \(v_x B_x - v_x B_x = 0\)). The last term, \(-D_x\phi\), is the divergence-cleaning term.

The divergence-cleaning \(\phi\)

The formula for qq_rslt.ph takes the following form, using ch_divb_square (\(c_h^2\)) and util::exp(-dt/tau_divb).

\[ \phi^{\text{rslt}} = \left[ \phi^n - \Delta t'\, c_h^2 \left( D_x B_x + D_y B_y + D_z B_z \right) \right] e^{-\Delta t'/\tau} \]

The term \(-c_h^2\nabla\cdot\boldsymbol{B}\) is advanced with the finite difference, and the damping term \(-\phi/\tau\) is handled by multiplying by the exact-solution factor \(e^{-\Delta t'/\tau}\).

Energy equation

// miso/include/miso/mhd_integrator_host.hpp: update_sc4 (excerpt; parts elided with ...)
// Et: total energy per unit volume
// ei: is the internal energy per unit mass
const Real Et =
    + qq_orgn.ro(i,j,k)*qq_orgn.ei(i,j,k)
    + 0.5*qq_orgn.ro(i, j, k)*( ... )   // v^2 of qq_orgn
    + pii8<Real>*( ... );               // B^2 of qq_orgn

qq_rslt.ei(i, j, k) = ( Et + dt * (
    -space_centered_4th(c_ht, c_vx, dxi, i, j, k, grid.is, 0, 0)
    -space_centered_4th(c_ht, c_vy, dyi, i, j, k, 0, grid.js, 0)
    -space_centered_4th(c_ht, c_vz, dzi, i, j, k, 0, 0, grid.ks)
    + pii4<Real>* (
        +space_centered_4th(c_vb, c_bx, dxi, i, j, k, grid.is, 0, 0)
        +space_centered_4th(c_vb, c_by, dyi, i, j, k, 0, grid.js, 0)
        +space_centered_4th(c_vb, c_bz, dzi, i, j, k, 0, 0, grid.ks)
        )
    )
    -0.5*qq_rslt.ro(i, j, k)*( ... )   // v^2 of qq_rslt (subtract the new kinetic energy)
    - pii8<Real>*( ... )               // B^2 of qq_rslt (subtract the new magnetic energy)
    + dt * source_ei(src, c_qq, i, j, k)
)/qq_rslt.ro(i, j, k);
\[ e^{\text{rslt}} = \frac{1}{\rho^{\text{rslt}}}\left\{ E^n + \Delta t' \left[ -D_j(h_t v_j) + \frac{1}{4\pi} D_j\big((\boldsymbol{v}\cdot\boldsymbol{B}) B_j\big) + S_E \right] - \frac{1}{2}\rho^{\text{rslt}} (v^{\text{rslt}})^2 - \frac{(B^{\text{rslt}})^2}{8\pi} \right\} \]

The total energy \(E\) is advanced, and the kinetic and magnetic energies of the new \(\rho\), \(\boldsymbol{v}\), and \(\boldsymbol{B}\) computed in the same stage are subtracted to recover the internal energy. Because it uses the new velocity and magnetic field, the energy equation is placed at the end of the loop. The source term source_ei is the right-hand side \(S_E\) of the total energy equation (per unit volume and unit time).

8.4 Extracting the source terms

source_ro(src, c_qq, i, j, k) and the like are functions generated by the macro MISO_MHD_DEFINE_SOURCE_TERM in miso/include/miso/mhd_integrator.hpp.

The following shows the body of the function that the macro generates for vy, with the macro argument name replaced by vy (expanded for explanation).

// miso/include/miso/mhd_integrator.hpp: MISO_MHD_DEFINE_SOURCE_TERM expanded for vy (excerpt)
if constexpr (impl_source::has_vy<Source, Real>::value) {
  return src.vy(qq, i, j, k);
} else {
  return Real(0);
}

If the source term type Source has a member function such as vy(FieldsView<const Real>, int, int, int) const, it is called; if not, 0 is returned (Chapter 6). mhd::EmptySourceTerm has nothing, so all 8 source terms are 0. The decision is made at compile time, so terms that are not defined add no computation. The source terms are called at each stage with c_qq (qq_argm), the value used to evaluate the right-hand side.

9. Boundary conditions and halo exchange

9.1 apply_boundary_condition

The following function is called after each stage of RK4 and after each direction of the artificial viscosity.

// miso/include/miso/mhd_integrator_host.hpp
template <typename BoundaryCondition>
void apply_boundary_condition(const BoundaryCondition &bc,
                              Fields<Real, backend::Host> &qq) {
  bc.apply(qq.view(), grid.const_view());
  halo_exchanger.apply(qq);
}

  1. Apply the physical boundary condition given by the user with bc.apply.

  2. Exchange ghost cell values with the neighboring ranks with halo_exchanger.apply (halo exchange, Chapter 8).

The argument form of bc.apply is a convention that the bc of Model must follow.

void apply(mhd::FieldsView<Real> qq, GridView<const Real> grid) const;

9.2 Halo exchange and periodic boundaries

Halo exchange is in HaloExchanger in miso/include/miso/mhd_halo_exchange.hpp. The CPU version of apply does the following for each direction, in the order \(x\) → \(y\) → \(z\).

  1. Pack the range to send (margin layers from the edge) of the 9 variables into a send Array4D.

  2. Send and receive on the positive and negative sides at the same time with MPI_Isend and MPI_Irecv, and wait with MPI_Waitall.

  3. Write the received values into the ghost cells.

Because the directions are exchanged one at a time, the range sent in the \(y\) direction includes the ghost cells received just before in the \(x\) direction. As a result, the corner ghost cells also receive correct values.

The neighboring ranks are x_procs_pos and the like, which mpi::Shape obtains with MPI_Cart_shift. In a direction made periodic with domain.periodic, the neighbor of the rank at one end of the domain is the rank at the opposite end. If the number of processes in that direction is 1, the neighbor is the rank itself. In other words, the halo exchange takes care of the periodic boundary. For a problem in which all directions are periodic, you can use mhd::EmptyBoundaryCondition, which does nothing, as bc.

At the end of a direction that is not periodic, the neighbor is MPI_PROC_NULL and no exchange takes place. Filling those ghost cells is the role of bc.apply.

9.3 Example of a physical boundary

The BoundaryCondition in demo/hd2d_rayleigh_taylor/src/main.cpp places symmetric boundaries at both ends in the \(y\) direction.

// demo/hd2d_rayleigh_taylor/src/main.cpp: BoundaryCondition::apply (excerpt)
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);
  bc::symmetric(btag, qq.vx, grid, Sign::Pos, Direction::Y, Side::INNER);
  bc::symmetric(btag, qq.vy, grid, Sign::Neg, Direction::Y, Side::INNER);
  ...
}

is_physical_boundary (miso/include/miso/boundary_condition.hpp) checks whether this rank has no neighboring rank on the negative side in the \(y\) direction (Side::INNER), that is, whether it is responsible for the edge of the domain. When the domain is decomposed with MPI, some ranks are not responsible for an edge, so this check is necessary.

symmetric copies the value at the mirror-image position across the boundary into the ghost cell. With Sign::Neg, the sign is reversed. Only the velocity vy perpendicular to the boundary is reversed, which makes the flow through the boundary 0.

10. Artificial viscosity

With the central difference alone, numerical oscillations occur near sharp changes such as shock waves. After RK4, MISO applies artificial viscosity (artificial diffusion) once, selectively smoothing only the neighborhood of sharp changes. The implementation is impl_host::ArtificialViscosity in miso/include/miso/mhd_artificial_viscosity_host.hpp, and the method is based on Rempel (2014, ApJ 789, 132) (see the comment in config_defaults.hpp).

10.1 Update in each direction

// miso/include/miso/mhd_integrator_host.hpp(excerpt)
template <typename EOS, typename BoundaryCondition>
void apply_artificial_viscosity(const Real dt, const EOS &eos,
                                const BoundaryCondition &bc,
                                Fields<Real, backend::Host> &qq) {
  artdiff.characteristic_velocity_eval(qq, eos);

  // x direction
  artdiff.update(qq, qq_rslt, Direction::X, dt);
  qq.copy_from(qq_rslt);
  apply_boundary_condition(bc, qq);

  ...  // same 3 lines for y and z (Direction::Y, Direction::Z)
}

First, compute the characteristic speed \(c = f_{cs} c_s + f_{vv}|\boldsymbol{v}| + f_{ca} c_A\) at all grid points (the coefficients are mhd.artificial_viscosity.cs_fac and so on in the configuration). Then update one direction at a time in the order \(x, y, z\), applying the boundary condition and halo exchange each time.

10.2 Diffusive flux

The update in one direction is the difference of the fluxes \(F_{i+1/2}\) at the cell boundaries \(i+1/2\).

\[ q_i^{\text{new}} = q_i - \frac{\Delta t}{\Delta x}\left(F_{i+1/2} - F_{i-1/2}\right) \]

The targets are \(\rho\), \(\rho\boldsymbol{v}\), \(\boldsymbol{B}\), \(\phi\), and \(\rho e\). The flux is built by two functions in miso/include/miso/limiter.hpp.

Function

Computation

dqq_eval(q_{i-1}, q_i, q_{i+1}, ep)

The slope \(\Delta q_i\). If \(\varepsilon(q_i - q_{i-1})\), \(\varepsilon(q_{i+1} - q_i)\), and \((q_{i+1} - q_{i-1})/2\) have the same sign, it is the one with the smallest absolute value; otherwise it is 0 (generalized minmod limiter). \(\varepsilon\) is ep in the configuration

flux_core(q_i, q_{i+1}, Δq_i, Δq_{i+1}, c, fh)

Builds the values on the left and right of the cell boundary, \(q_L = q_i + \Delta q_i/2\) and \(q_R = q_{i+1} - \Delta q_{i+1}/2\), and returns \(F = -\frac{1}{2} c\, \Phi\, (q_R - q_L)\)

\(\Phi\) is defined using \(r = \min\left(1, (q_R - q_L)/(q_{i+1} - q_i)\right)\): it is \(\max(0, 1 + h(r - 1))\) if \(r \ge 0\) and 0 if \(r < 0\). \(h\) is fh in the configuration. Where the variation is smooth, \(q_R \approx q_L\), so the flux is almost 0. Diffusion acts only where there is a sharp change and \(q_R - q_L\) remains large.

For the energy, the total energy is updated with the flux of \(\rho e\), and then the new kinetic and magnetic energies are subtracted to recover the internal energy. The kinetic and magnetic energy lost to artificial viscosity is transferred to the internal energy.

11. Output and restart

11.1 Output times

// miso/include/miso/mhd_model_base.hpp
void save_if_needed() {
  if (time.time >= time.dt_output * time.n_output) {
    save_state();
    time.log();
    time.n_output++;
  }
}

Output is written at the first step at which the time reaches dt_output * n_output or more. The time step is not adjusted to match the output time, so the output time is not exactly a multiple of dt_output (Exercise 6). time.log() prints lines like the following on rank 0 (an example running demo/mhd2d_vortex with tend: 0.03 on 2 processes).

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
time =   0.03;  n_step =       13;  n_output =        3

11.2 Files that are written

All of them are written under io.save_dir (by default, data/ in the same folder as the configuration file).

File

Function that writes it

Contents

config.yaml

Config::save

The configuration with defaults filled in

grid.bin

Grid::save

Coordinates of the global grid

mpi/coords.csv

mpi::Shape::save

MPI coordinates of each rank

time/time.NNNNNNNN.txt

Time::save (miso/include/miso/time.hpp)

Time, output number, and step count, one per line

time/n_output.txt

Time::save

The latest output number

mhd/mhd.<n_output>.<rank>.bin

Checkpoint::save (miso/include/miso/mhd_checkpoint.hpp)

The 9 variables of each rank

Checkpoint::save copies qq on the compute location to an array on the CPU (copy_from), writes the element size (4 for float) at the beginning as a std::uint32_t, and then writes out the arrays, including ghost cells, as they are, in the order ro, vx, vy, vz, bx, by, bz, ei, ph. The layout is the in-memory layout of Array3D (the last index k is contiguous). In the example above, the grid of each rank is \((64 + 2\cdot 2) \times (128 + 2\cdot 2) \times 1\), so the file size is \(4 + 68 \times 132 \times 9 \times 4 = 323140\) bytes. In Python, read it with pymiso.Data (Chapter 1).

11.3 Restart

If io.continue is true (the default is also true) and time/n_output.txt remains, run() creates the initial condition, then calls load_state() and continues the calculation from the latest output. Time::load reads the time file on rank 0 and distributes it with MPI_Bcast, and in Checkpoint::load each rank reads its own mhd.<n_output>.<rank>.bin. The design uses the output as is as the data for restarting (a checkpoint). To recompute from the beginning, delete data/ or set io.continue: false.

12. Where to make changes when creating your own problem

In most cases, you create a new problem by copying a demo’s main.cpp and rewriting the 4 components of Model and config.yaml. You do not need to modify anything inside miso/include/miso/.

What you want to change

Where to write it

Convention (the form in which it is called)

Demo to refer to

Initial condition

Type of ic

void apply(mhd::FieldsView<Real> qq, GridView<const Real> grid) const. Fills in values at all grid points including ghost cells. Called only once, on arrays on the CPU

demo/mhd2d_vortex

Boundary condition

Type of bc

void apply(mhd::FieldsView<Real> qq, GridView<const Real> grid) const. Fills the ghost cells at the ends of directions that are not periodic. Called many times per step, on arrays on the compute location

demo/hd2d_rayleigh_taylor

Source term

Type of src

Define only the terms you need, in the form __host__ __device__ Real vy(mhd::FieldsView<const Real> qq, int i, int j, int k) const (ro, vx to vz, bx to bz, ei). If there is an upper limit on the time step, define Real dt_limit() const

demo/hd2d_rayleigh_taylor, demo/mhd2d_reconnection

Equation of state

Type of eos

Provide gas_pressure(btag, qq, pr) and sound_speed(btag, qq, cs) with the same arguments as IdealEOS

miso/include/miso/eos.hpp

How the time step is determined

Model::update()

Override it in the same form as ModelBase::update()

demo/mhd3d_magnetosphere

Problem parameters

A custom section in config.yaml

Read them in the constructor, as in config["rayleigh_taylor"]["g_grav"].as<Real>()

demo/hd2d_rayleigh_taylor

Points to note are as follows.

  • In the GPU version, boundary conditions and source terms are called on arrays on the GPU. Attach __host__ __device__ to the member functions of the source term and use only values that the GPU can read (numeric members and views of arrays on the compute location) (Chapter 9).

  • If you write a for statement directly in a boundary condition, it works only on the CPU. Use functions that take the compute-location tag btag, such as boundary_condition::symmetric, or use for_each.

  • If you make all directions periodic, set all of domain.periodic in config.yaml to true, and you can use mhd::EmptyBoundaryCondition<Real> as bc.

How to write each component, with examples, is summarized in Customization in the MISO documentation.

13. Exercises: checking by reading the source

Exercise 1: How many times boundary conditions are called in one step

When advancing one step with the settings of demo/mhd2d_vortex, read miso/include/miso/mhd_integrator_host.hpp and count how many times bc.apply is called. Also count the number of calls made when creating the initial condition, and the number of times MPI_Allreduce is called in one step.

::

It is called 7 times in one step.

  • In the 4 stages of runge_kutta_4step, apply_boundary_condition is called once in each (4 times).

  • In apply_artificial_viscosity, it is called once in each of the \(x, y, z\) directions (3 times).

apply_boundary_condition calls halo_exchanger.apply after bc.apply, so the halo exchange also happens 7 times in one step. Even in a 2D calculation, the call for the artificial viscosity in the \(z\) direction is not skipped, so the count is the same.

When creating the initial condition, it is called once at the end of MHD::apply_initial_condition.

MPI_Allreduce is called once in Integrator::cfl. If the source term has dt_limit(), it is called once more in ModelBase::update.

You can also verify the counts by using, as bc, a boundary condition that counts the number of calls inside apply, and calling update() directly. Since apply is a const member function, it cannot modify member variables inside. Put the counting variable outside the boundary condition struct (for example, a variable declared outside any function in main.cpp). Counted this way, you get 1 for the initial condition, 7 for one step, and 14 for two steps.

Exercise 2: Filling in default values

mhd.cfl_number is not written in demo/hd2d_rayleigh_taylor/config.yaml. Read config_defaults.hpp and merge_yaml and answer what the CFL number and the artificial viscosity parameter ep are in this calculation. Also explain in which file, after the run, you can check the values that were actually used.

::

The CFL number is the default value, 0.5, and ep is 2.0, as written by the user.

merge_yaml descends into the mhd section of the defaults, and since cfl_number is not in the user’s settings, it adds 0.5. artificial_viscosity is a section (a YAML map), so it descends into it recursively and fills in key by key. ep is written by the user, so it is not overwritten. mhd.resistivity.cfl_number, which the user did not write, is also filled in with 0.5.

The completed result remains in data/config.yaml, which save_metadata() writes to the output destination. In the settings of demo/hd2d_rayleigh_taylor, io.save_dir is the default data/, so the file is demo/hd2d_rayleigh_taylor/data/config.yaml.

Exercise 3: Gravity source term

The SourceTerm in demo/hd2d_rayleigh_taylor/src/main.cpp provides a uniform gravity (magnitude \(g\)) in the \(-y\) direction. Explain which term of which formula in update_sc4 changes when this source term is used. Also, explain from the energy equation in Section 8.3 what must be returned as the source term for ei in order to conserve the total energy.

::

SourceTerm defines vx, vy, vz, and ei, and the only one that returns a nonzero value is \(-g\rho\) for vy.

// demo/hd2d_rayleigh_taylor/src/main.cpp
__host__ __device__ inline Real vy(mhd::FieldsView<const Real> qq, int i, int j,
                                   int k) const noexcept {
  return -g_grav * qq.ro(i, j, k);
}

What changes is the +source_vy(src, c_qq, i, j, k) term in the equation of motion in the \(y\) direction.

\[ \rho^{\text{rslt}} v_y^{\text{rslt}} = \rho^n v_y^n + \Delta t' \left[ \cdots - g\rho \right] \]

\(\rho\) is the value in c_qq, that is, qq_argm of that stage, and it is evaluated in each of the 4 RK4 stages. ro, bx, by, and bz are not defined, so they become 0 through if constexpr.

Regarding the energy: MISO’s energy equation is the equation for the total energy \(E\), and the internal energy is recovered with \(e = (E - \rho v^2/2 - B^2/8\pi)/\rho\). If the work that gravity does on the fluid, \(\rho\boldsymbol{v}\cdot\boldsymbol{g} = -g\rho v_y\), is not included in \(S_E\), the kinetic energy gained from gravity is subtracted without having been added to \(E\), and the internal energy decreases. To conserve the total energy, return \(-g\rho v_y\) in ei.

// example of the source term for ei (illustrative)
__host__ __device__ inline Real ei(mhd::FieldsView<const Real> qq, int i, int j,
                                   int k) const noexcept {
  return -g_grav * qq.ro(i, j, k) * qq.vy(i, j, k);
}

The ei in demo/hd2d_rayleigh_taylor returns 0, so in this demo the work done by gravity is deducted from the internal energy side.

When adding gravity, besides the source term you also need to consider making the initial condition satisfy hydrostatic equilibrium (in this demo, \(p = p_0 - g\rho y\)) and the boundary conditions at the ends in the \(y\) direction.

Exercise 4: Differences in the z direction in a 2D calculation

demo/mhd2d_vortex is a 2D calculation with k_size: 1. Read AxisGrid in miso/include/miso/grid.hpp and explain what happens in this case to the \(z\)-direction difference in update_sc4, the CFL condition, the halo exchange, and the \(z\)-direction update of the artificial viscosity.

::

In a direction where the number of grid points size is 1, AxisGrid sets stride = 0, margin = 0, and total = 1, sets the grid spacing ds to 1.e30, and sets its reciprocal dsi to 0. grid.ks is the stride, so it is 0; grid.k_margin is 0; and grid.k_total is 1.

  • Difference: In space_centered_4th(..., 0, 0, grid.ks), ks = 0, so all 4 points are the same grid point k, and the combination \((-1 + 8 - 8 + 1)\) gives 0. In addition, dzi is also 0. The derivative terms in the \(z\) direction become 0.

  • CFL condition: In min3(dx[i], dy[j], dz[k]), dz is 1.e30, so the \(z\) direction does not restrict the time step.

  • Halo exchange: The \(z\) direction of HaloExchanger is not performed because grid.k_size > 1 is false (the GPU version decides with grid.k_total > 1).

  • Artificial viscosity: For Direction::Z, ks = 0, so the left and right values are all the same, and the flux difference \(F_{i+1/2} - F_{i-1/2}\) is 0. The dzi it is multiplied by is also 0, so the values do not change.

In every case the loop range is just k = 0, and no out-of-range array access occurs.

Exercise 5: Checking the accuracy of the fourth-order central difference (advanced)

Write a program that uses MISO, compute the derivative of \(f(x) = \sin x\) with mhd::space_centered_4th, and find the maximum difference from \(\cos x\). Divide \(0 \le x < 2\pi\) into \(n\) cells and place the values at the cell centers \(x_i = (i + 1/2)\Delta x\). Check how the error changes for \(n = 16, 32, 64\). For the build, use the “minimal CMake configuration for using MISO from your own program” in Chapter 1.

::

What to build is as follows.

  1. Write #include <miso/core.hpp> and #include <miso/mhd.hpp> (space_centered_4th is included from mhd.hpp).

  2. Create an Array3D<double> with 2 ghost cells on each side (of size n + 4, 1, 1), and put \(\sin x\) at the position \(x = (i - 2 + 1/2)\Delta x\) of index i. If you fill the ghost cells with the same formula, the result is the same as a periodic boundary.

  3. For i = 2 to n + 1, find the maximum of the difference between mhd::space_centered_4th(f.const_view(), 1.0 / dx, i, 0, 0, 1, 0, 0) and \(\cos x\).

space_centered_4th does not use MPI, so this program works even without creating an Env.

An example of the result is as follows.

n = 16  error = 0.000763339
n = 32  error = 4.90805e-05
n = 64  error = 3.08927e-06
ratio 16/32 = 15.5528
ratio 32/64 = 15.8874

Doubling the number of grid points reduces the error to about \(1/16 = 1/2^4\), which shows that the scheme is fourth-order accurate. The leading term of the error of the difference formula is \(-\frac{\Delta x^4}{30} f^{(5)}\), and for \(f = \sin x\) its maximum is \(\Delta x^4/30\). For \(n = 16\), \(\Delta x = \pi/8\), so \((\pi/8)^4/30 \approx 7.9\times 10^{-4}\), which agrees with the result.

Exercise 6: Output times

When you run demo/mhd2d_vortex with dt_output: 0.01, the first line of data/time/time.00000001.txt is a value slightly off from 0.01, such as 0.0100069. Read ModelBase::save_if_needed and ModelBase::run and explain why this offset occurs.

::

run() calls save_if_needed() after every step and writes output at the first step at which time.time >= dt_output * n_output holds. The time step \(\Delta t\) is determined only by the CFL condition and is not shortened to land exactly on the output time. As a result, the time that is written is that of the step that first exceeds 0.01, which is later than 0.01 by at most \(\Delta t\).

The exact value of the output time is written on the first line of time/time.NNNNNNNN.txt. When you use the time in analysis, use the value in this file (the value that pymiso reads), not a multiple of dt_output.

Contents · Previous chapter · Fortran Correspondence Table