5. Inheritance and CRTP

Every MISO demo’s main.cpp contains a definition of the following form (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() {}
};

The time-evolution loop, the output, and the MHD computation are handled by ModelBase in the MISO library, and the parts that differ from problem to problem (the equation of state eos, the initial condition ic, the boundary condition bc, and the source term src) are supplied by the demo’s Model. The goal of this chapter is for you to be able to read this first line. We learn, in order, the basics of inheritance, virtual functions for comparison, and CRTP, which MISO uses.

5.1 Basics of inheritance
5.2 Functions with the same name and calling base-class functions
5.3 Run-time switching with virtual functions (for comparison)
5.4 CRTP: passing the derived-class type to the base class
5.5 MISO’s ModelBase
5.6 The demo’s Model
5.7 Overriding update()
Exercises

In this chapter too, read the angle brackets <> as “a marker that specifies a type later” (Chapter 6).

5.1 Basics of inheritance

Base class and derived class

We prepare a class that manages the time and the step count, and use it to build a class that solves the ordinary differential equation \(du/dt = -u\).

class Simulation {
public:
  explicit Simulation(double tend) : tend(tend) {}

  double time = 0.0;  // current time
  double tend;        // end time
  int n_step = 0;     // number of steps advanced

protected:
  void advance(double dt) {  // advances the time and the step count
    time += dt;
    ++n_step;
  }
};

class Decay : public Simulation {
public:
  double u;  // member added in Decay

  Decay(double tend, double u0) : Simulation(tend), u(u0) {}

  void step(double dt) {
    u -= dt * u;   // advances du/dt = -u by one step with the Euler method
    advance(dt);   // can call the protected function of the base class
  }
};

class Decay : public Simulation is a declaration that “defines Decay after inheriting all the members of Simulation.” This relationship is called inheritance; the original Simulation is called the base class, and the newly created Decay is called the derived class. An object of Decay has both the Simulation part (time, tend, n_step) and the part added in Decay (u).

Decay d(1.0, 1.0);
d.step(0.125);
d.step(0.125);
std::cout << d.time << " " << d.n_step << " " << d.u << "\n";  // 0.25 2 0.765625

Simulation &s = d;  // refers to the Simulation part of d
std::cout << s.time << " " << s.tend << "\n";                  // 0.25 1
// s.u = 2.0;       // error: Simulation has no u

Inherited members can be used just like members of Decay, as in d.time. When you inherit with public, you can handle an object of the derived class through a reference (Simulation &) or pointer to the base class. Through s, only the members of Simulation are visible.

The public in the inheritance makes the public members of the base class remain public in the derived class. For a derived class defined with struct, the inheritance is public if you write nothing, but it is usual to write public explicitly, as in the MISO demos.

protected

advance is written after protected:. protected is an access level that cannot be used from outside but can be used from member functions of the class itself and of derived classes. You can call advance(dt) from inside Decay::step, but writing d.advance(0.1); in main() produces the error 'void Simulation::advance(double)' is protected within this context. private members cannot be used from the derived class either.

Calling the base-class constructor

Simulation(tend) is written at the beginning of the member initializer list of the Decay constructor. This specifies that the base-class part is initialized with the base-class constructor. Simulation has no constructor without parameters, so omitting this specification causes an error.

Initialization happens in the following order.

The base-class part (Simulation(tend); inside it, the members of Simulation are initialized in declaration order)
The members added in the derived class (in declaration order; here, u)
The body of the derived-class constructor

The base class is initialized first no matter where it is written in the initializer list. Cleanup happens in the reverse order: the derived-class destructor is called first, then the base-class destructor.

Simulation(double tend) : tend(tend) is a way of initializing the member tend with the parameter of the same name, tend (the same as gm(gm) of IdealEOS in Chapter 4).

A MISO example: inheritance and member initialization in Model

We read the Model at the beginning of the chapter using the terms of this section.

The base class of struct Model : public mhd::ModelBase<Model, Real, Backend> is mhd::ModelBase<Model, Real, Backend>. The reason Model itself appears inside the angle brackets is explained in Section 5.4.
Model inherits the members of ModelBase (config, mpi_shape, time, grid, exec_ctx, mhd) and adds eos, ic, bc, and src.
ModelBase(config) in the initializer list is the call to the base-class constructor. Inside Model, you can write the base class as ModelBase, without angle brackets.

The order of initialization tells us the following. The constructor of Model in demo/hd1d_shock_tube/src/main.cpp is as follows.

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

mpi_shape in bc(mpi_shape) is a member of the base class ModelBase. The base-class part is initialized first, so mpi_shape already exists when bc is created. Similarly, the eos passed to ic(config, eos) is a member declared before ic, so it exists before ic. As we read in Chapter 4, InitialCondition holds eos by reference, so this order matters.

5.2 Functions with the same name and calling base-class functions

We put the time-evolution loop run() and the default one-step update() in the base class. In the derived class, we want to define an update() with the same name and replace it with our own computation.

class Simulation {
public:
  explicit Simulation(double tend) : tend(tend) {}

  double time = 0.0;
  double tend;
  int n_step = 0;

  void update() { advance(0.25); }  // default single step: advances only the time

  void run() {
    while (time < tend) {
      update();  // note which update() is called
    }
  }

protected:
  void advance(double dt) {
    time += dt;
    ++n_step;
  }
};

class Decay : public Simulation {
public:
  double u;

  Decay(double tend, double u0) : Simulation(tend), u(u0) {}

  void update() {  // defines a function with the same name as in the base class
    const double dt = 0.125;
    u -= dt * u;
    advance(dt);
  }
};

Decay a(1.0, 1.0);
a.update();              // Decay::update() is called
a.Simulation::update();  // calls the base class update() by name
std::cout << a.time << " " << a.n_step << " " << a.u << "\n";  // 0.375 2 0.875

Decay b(1.0, 1.0);
b.run();                 // update() inside run() is Simulation::update()
std::cout << b.time << " " << b.n_step << " " << b.u << "\n";  // 1 4 1

The first result tells us the following.

a.update() calls Decay::update(). If a derived class defines a function with the same name, the derived-class function is used when it is called from an object of the derived class. The base-class function with the same name is hidden (name hiding).
You can still call the hidden base-class function by adding class_name::, as in Simulation::update(). From inside a member function of the derived class, you can use this to call the base-class processing first and then add more processing (Section 5.7).

The second result is different from what we expected. b.run() ends after 4 steps, and u stays at 1. Decay::update() has never been called.

The reason is that run() is a member function of Simulation. At the point where Simulation::run() is written, the only update() the compiler knows is Simulation::update(). The update(); inside run() means this->update(); (Chapter 4), and the type of this is Simulation *, so Simulation::update() is called. Simulation does not know that a class inherits from it, nor that such a class has an update().

To call a derived-class function (Decay::update()) from a common procedure written in the base class (run()), you need a different mechanism. C++ has two ways: virtual functions, which switch at run time, and CRTP, which decides at compile time. MISO uses CRTP. For comparison, we first explain virtual functions briefly.

5.3 Run-time switching with virtual functions (for comparison)

Change Simulation and Decay from the previous section in just the following three places.

class Simulation {
public:
  // ...
  virtual void update() { advance(0.25); }  // add virtual
  virtual ~Simulation() = default;          // add a virtual destructor
  // ... (run() and advance() are the same)
};

class Decay : public Simulation {
public:
  // ...
  void update() override {  // add override
    // ... (the body is the same)
  }
};

Then the result of b.run() becomes 1 8 0.343609, and Decay::update() is called 8 times (\(0.875^8 \approx 0.3436\)).

A member function with virtual is called a virtual function. When you call a virtual function, the function that matches the actual type of the object (here, Decay) is selected at run time, rather than the declared type of the pointer or reference (here, Simulation *). Redefining a virtual function in a derived class is called overriding, and the trailing override makes the compiler check that a corresponding virtual function exists in the base class. We also add virtual to the destructor in case a derived-class object is destroyed through the base class.

For the run-time selection, an object of a class with virtual functions gets a hidden pointer to “a table of which function to call,” and the table is consulted at every call (this is how many compilers implement it). Thanks to this mechanism, you can, for example, read a configuration file and then decide at run time which derived-class object to create.

MISO does not use virtual functions. virtual does not appear in the sources under miso/include/miso/ or demo/. What MISO needs is to “build parts that are fixed at compile time for each demo into a common procedure,” and it does not need run-time switching. Moreover, the way MISO’s parts are used has points that do not fit well with virtual functions.

ModelBase uses not only the member functions of the derived class but also its member variables eos, ic, bc, and src. Their types differ from demo to demo (IdealEOS<Real>, InitialCondition, EmptyBoundaryCondition<Real>, and so on) and they have no common base class. To write this with virtual functions, you would need to prepare a common base class for each part.
Source terms and similar parts are used from inside functions that run on the GPU (Chapter 9), and for that purpose objects created on the CPU side are copied to the GPU. In CUDA, copying an object of a class with virtual functions between the CPU and the GPU (including passing it as an argument of a GPU function) is undefined behavior (Polymorphic Classes in the CUDA Programming Guide).

CRTP in the next section resolves these issues by “deciding types at compile time.”

5.4 CRTP: passing the derived-class type to the base class

Mechanism

The problem in Section 5.2 was that the base class Simulation does not know the derived class Decay. In that case, if we tell the base class the derived-class type, the base class can call the derived-class functions. To pass a type in later, we use templates (Chapter 6).

template <class Derived>  // Derived is the type of the derived class, specified later
class Simulation {
public:
  explicit Simulation(double tend) : tend(tend) {}

  double time = 0.0;
  double tend;
  int n_step = 0;

  void update() { advance(0.25); }  // default single step

  void run() {
    while (time < tend) {
      derived().update();  // calls the derived class's update()
    }
  }

protected:
  void advance(double dt) {
    time += dt;
    ++n_step;
  }

  // treats *this as a reference to the derived class
  Derived &derived() { return static_cast<Derived &>(*this); }
};

class Decay : public Simulation<Decay> {  // passes its own type to the base
public:
  double u;

  Decay(double tend, double u0) : Simulation(tend), u(u0) {}

  void update() {
    const double dt = 0.125;
    u -= dt * u;
    advance(dt);
  }
};

class Clock : public Simulation<Clock> {  // does not define update()
public:
  explicit Clock(double tend) : Simulation(tend) {}
};

Decay d(1.0, 1.0);
d.run();  // Decay::update() is called 8 times
std::cout << d.time << " " << d.n_step << " " << d.u << "\n";  // 1 8 0.343609

Clock c(1.0);
c.run();  // the base update() is called 4 times
std::cout << c.time << " " << c.n_step << "\n";                // 1 4

Without using virtual, we got the same result as with the virtual function. Let us read it in order.

template <class Derived> is a marker that Simulation receives one type later. Inside the definition, the received type is used under the name Derived.
In class Decay : public Simulation<Decay>, Decay itself is passed as Derived. This form, in which the derived class itself is written inside the angle brackets of the base class, is called CRTP (Curiously Recurring Template Pattern).
Inside Simulation<Decay>, Derived is Decay. derived() returns a Decay &, and derived().update() calls Decay::update().
Clock does not define update(). derived().update() looks for update() as a member of Clock, and if Clock does not have it, it finds the base class’s Simulation<Clock>::update(). If the derived class does not define it, the default processing is used; if it does, the derived-class processing is used.

static_cast<Derived &>(*this)

Let us read the contents of derived().

*this is the object that called run(). run() is a member function of Simulation<Decay>, so the type of *this is Simulation<Decay>.
However, the Simulation<Decay> part is actually a part of a Decay object (it is used so that only Decay inherits from it).
static_cast<Decay &>(*this) is a conversion that means “treat this Simulation<Decay> as a reference to the Decay that contains it.” static_cast<type>(expression) is the notation for converting types at compile time, the same as static_cast<size_t>(...) that we saw in Chapter 4.

This conversion works correctly only when the object really is a Decay. The compiler does not check at run time. In CRTP, you make this always hold by keeping the convention of passing your own type, as in class Decay : public Simulation<Decay>.

Comparison of virtual functions and CRTP

Point of comparison	Virtual functions	CRTP
When the function to call is decided	Run time	Compile time
How to choose the derived class	Can be chosen at run time (after reading the configuration)	Determined by the types written in the source
What the base class can use	Virtual functions declared in the base class	Member functions and member variables of the derived class (including their types)
Type of the base class	One common type, `Simulation`	A separate type for each derived class: `Simulation<Decay>`, `Simulation<Clock>`
Use in MISO	Not used	Used in `ModelBase`

In CRTP, the function to call is decided at compile time, so there is no cost of consulting a table, and the compiler can also expand the function at the call site (inlining). If the derived class lacks a needed member, you get a compile error before running. On the other hand, Simulation<Decay> and Simulation<Clock> are different types, so you cannot make “a single variable that can hold either one.” MISO demos solve one problem per executable, so this restriction is not an issue.

5.5 MISO’s ModelBase

The following is an excerpt of ModelBase in miso/include/miso/mhd_model_base.hpp. update() is read in Section 5.7. The functions for saving and loading are omitted.

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) {}

  // ... (update(), save_metadata(), save_state(), load_state(), save_if_needed())

  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'.");
    // ... (the same check for ic, bc, and src)

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

protected:
  Derived &derived() { return static_cast<Derived &>(*this); }
};

Template arguments and members

There are three items in the angle brackets. Derived is the CRTP derived class (the demo’s Model), Real is the real-number type, and Backend is where the computation runs (CPU or GPU).
The members are the configuration config (a reference), the MPI domain decomposition mpi_shape, the time time (the Time of Chapter 4), the grid grid, the information used to run computations exec_ctx, and the MHD variables and computation mhd.
The constructor’s initializer list is in declaration order, and it uses earlier members to initialize later ones (such as grid(cfg, mpi_shape)).

The flow of run()

run() is the procedure for the whole MISO computation.

auto &d = derived(); obtains a reference d to the derived class Model. auto & is a reference declaration whose type is determined from the right-hand side (auto in Chapter 2).
static_assert(...) checks at compile time that Model has the member variables eos, ic, bc, and src. The mechanism is explained in Chapter 6.
mhd.apply_initial_condition(d.ic, d.bc) sets the initial condition using the derived class’s members ic and bc. We read in Chapter 4 that ic.apply(...) is called inside.
If io.continue in the configuration is true and there is output from a previous run, the program loads it in order to continue the computation from there.
Until the time reaches the end time, it advances one step with d.update() and writes output when needed.

This has the same form as Simulation<Derived>::run() in Section 5.4. d.update() calls Model::update() if Model defines update(), and ModelBase::update() if not. The base class uses the derived class’s member variables, such as d.ic and d.bc, with their actual types; this is what is hard to write with virtual functions.

The template in .template as<bool>() inside the if (...) is a specification required when calling a function template inside a template. When reading, you may think of it as the same as .as<bool>() (Chapter 6).

If Model lacks a required member, you get a compile error. For example, if you forget to write src, static assertion failed: Derived must have member variable 'src'. appears in the g++ error output (another error saying that src does not exist also appears before it).

5.6 The demo’s Model

We now read again the Model and main() in demo/mhd2d_vortex/src/main.cpp at the beginning of the chapter.

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

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

Member	Type	Role
`eos`	`eos::IdealEOS<Real>`	The equation of state for an ideal gas. The ratio of specific heats is `eos.gm` in the configuration file
`ic`	`InitialCondition` (defined in the demo)	The initial condition of the vortex. `apply` writes values into the `FieldsView`
`bc`	`mhd::EmptyBoundaryCondition<Real>`	A boundary condition that does nothing. All directions are periodic boundaries, so MPI communication alone is enough
`src`	`mhd::EmptySourceTerm<Real>`	No source term

EmptyBoundaryCondition and EmptySourceTerm are structs with almost no contents, at the end of miso/include/miso/mhd_model_base.hpp. EmptyBoundaryCondition has only an apply that does nothing.

template <typename Real> struct EmptyBoundaryCondition {
  void apply(mhd::FieldsView<Real>, GridView<const Real>) const {}
};

The parameters are unnamed because they are not used. From ModelBase’s point of view, all it needs is that bc.apply(...) can be called; it does not care what the type of the boundary condition is. To create your own boundary condition, write a struct with an apply of the same form and change the type of the Model member bc to it (BoundaryCondition in demo/hd1d_shock_tube/src/main.cpp is an example).

The flow of main() is: set up MPI (Env, Chapter 4), read the configuration file (Config), create the Model, and call model.run(). Model has no run(), so model.run() calls the base class’s ModelBase::run(), and the parts of Model are used inside it.

5.7 Overriding update()

ModelBase’s update()

ModelBase::update() is the default processing for one step. We extract it from miso/include/miso/mhd_model_base.hpp, omitting the handling for when the source term has an upper limit on the time step.

  /// @brief Default implementation of one timestep update.
  /// @details The derived class may provide its own `update()` method.
  void update() {
    auto &d = derived();
    // ... (preparation and checks for detecting the source term type)
    auto dt = mhd.cfl(d.eos);
    // ... (if the source term has the time step upper limit dt_limit(), limits dt to at most that)
    mhd.update(dt, d.eos, d.bc, d.src);
    time.update(dt);
  }

mhd.cfl(d.eos) determines the time step dt from the CFL condition.
mhd.update(dt, d.eos, d.bc, d.src) advances the MHD variables by one step. The boundary condition and the source term use the derived class’s parts.
time.update(dt) advances the time and the step count (Time::update in Chapter 4).

As the original comment The derived class may provide its own update() method. says, a derived class can define its own update(). The d.update() in run() calls the derived class’s update() by the same mechanism as Decay in Section 5.4. This is not overriding of a virtual function but a combination of the name hiding of Section 5.2 and CRTP. In this course, we call defining a function with the same name in a derived class to replace the processing “overwriting.”

A MISO example: the magnetosphere demo

The Model in demo/mhd3d_magnetosphere/src/main.cpp overwrites update().

struct Model : public mhd::ModelBase<Model, Real, Backend> {
  eos::IdealEOS<Real> eos;
  InitialCondition ic;
  BoundaryCondition bc;
  SourceTerm src;
  TimeStep timestep;

  void update() {
    const auto dt = timestep.cfl(mhd.qq, mhd.grid, eos);
    mhd.update(dt, eos, bc, src);
    time.update(dt);
  }

  Model(Config &config)
      : ModelBase(config), eos(config), ic(config, eos),
        bc(config, mpi_shape, grid, eos, ic), src(config, mhd.grid),
        timestep(config, mhd.grid) {}
};

The only difference from ModelBase::update() is how the time step is determined. TimeStep::cfl (demo/mhd3d_magnetosphere/src/timestep.hpp) determines dt by excluding the region inside the planet (inside the radius magnetosphere.radius) from the CFL condition calculation. Because this is inside the derived class, you can use base-class members such as mhd and time, and your own members eos, bc, and src, as they are. You do not need to go through derived().

If Model defines run(), model.run() calls Model::run() (name hiding). However, no MISO demo overwrites run(). The time-evolution loop and the output procedure are left to ModelBase::run(), and only the contents of one step are overwritten when needed.

Exercises

Exercise 1: Initialization order of Model members

For the Model in demo/hd1d_shock_tube/src/main.cpp, explain the following.

The order in which the members of ModelBase (config, mpi_shape, time, grid, exec_ctx, mhd) and the members of Model (eos, ic, bc, src) are initialized when model is created.
The problem that can occur if you swap the order of the declarations of eos and ic in the declaration of Model.

::

The base-class part is initialized first, with config, mpi_shape, time, grid, exec_ctx, and mhd in declaration order. Then the members of Model are initialized in declaration order: eos, ic, bc, and src. Finally, the body of the Model constructor (empty) runs. The reason bc(mpi_shape) can be used is that mpi_shape is a member of the base class and has already been initialized.
ic would be initialized before eos. Then ic(config, eos) would receive the not-yet-constructed eos as a reference. The current constructor of InitialCondition only records eos as a reference and does not read its value, so in practice eos exists by the time apply is called. However, if you changed the constructor to read eos.gm, it would read a member that holds no value. The order in the initializer list (eos, ic) would not change, so this is an easily overlooked mistake. With -Wall, g++ warns about the difference in order (-Wreorder) and about using the uninitialized eos (-Wuninitialized).

Exercise 2: Adding processing to the default update()

Add an update() to the Model in demo/mhd2d_vortex/src/main.cpp so that, after running the default one step of ModelBase, it prints the step count and the time every 10 steps. Build it with the CMake setup from Chapter 1 for using MISO from your own program, run it with 1 MPI process, and check the result. For the configuration file, start from demo/mhd2d_vortex/config.yaml, set the number of grid points to 64 x 64, the mpi decomposition to 1 x 1 x 1, time.tend to 0.1, and io.enabled to false; then it finishes quickly and writes no files.

::

Add the following function to Model (it uses std::cout for printing, so also add #include <iostream>).

  void update() {
    ModelBase::update();  // default single step (computes dt, updates MHD, advances time)
    if (time.n_step % 10 == 0) {
      std::cout << "n_step = " << time.n_step << ", time = " << time.time << "\n";
    }
  }

ModelBase::update() is a way of calling a hidden base-class function by name, like a.Simulation::update() in Section 5.2. Inside Model, you can write the base class as ModelBase.
The d.update() in ModelBase::run() calls this Model::update().
Note that if you write update(); inside update() without the ModelBase::, it keeps calling Model::update() itself, and the program never stops (it finally terminates abnormally from running out of stack).

This is an example of a result from running on the CPU (float) with a 64 x 64 grid, time.tend: 0.1, and time.dt_output: 0.05.

n_step = 10, time = 0.0496095
n_step = 20, time = 0.0984416

The time step is determined at every step by the CFL condition, so the time values vary with the number of grid points and the settings.

Exercise 3: Which function CRTP calls

For Simulation<Derived>, Decay, and Clock in Section 5.4, explain which class’s update() is called by derived().update() in run() for each of the following three calls. In the third, assume that the definition of update() has been removed from Decay.

Decay d(1.0, 1.0); d.run();
Clock c(1.0); c.run();
(with Decay::update() removed) Decay d(1.0, 1.0); d.run();

::

It is Decay::update(). Inside Simulation<Decay>, Derived is Decay, so derived() returns a Decay &, and the update() that is a member of Decay is found.
It is Simulation<Clock>::update(). derived() returns a Clock &, but Clock does not define update(), so the base class’s update() is used.
It is Simulation<Decay>::update(), for the same reason as in 2. The result is that time is 1 and n_step is 4, and u does not change.

In every case, the function to call is decided at compile time. The same is true for MISO’s Model: if it defines update(), Model::update() is called; if not, ModelBase::update() is called.