6. Templates

Contents

Introduction

MISO’s computations run in single precision (float) or double precision (double), on a CPU or a GPU. The demo’s main.cpp chooses among them with the following two lines (excerpt from demo/mhd2d_vortex/src/main.cpp).

using Real = float;

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

The library side (miso/include/miso/) does not have separate code for float and double, or for CPU and GPU. Types are left as “names to be specified later,” and when the user provides Real and Backend, the compiler creates the concrete functions and classes. This mechanism is called a template.

In Chapter 5, we used only the act of passing types into the angle brackets of ModelBase<Model, Real, Backend>. This chapter starts with how to write templates for functions, classes, and variables, then reads how MISO selects CPU or GPU processing according to Backend, and how it checks at compile time whether a model has the required members. You do not need to be able to write the last Section 6.7 yourself; being able to read MISO’s source is enough.

6.1 Function Templates

The expression for computing a square is x * x for both integers and real numbers. In the overloading of Chapter 2, we defined separate functions for int and double. If the form of the processing is the same regardless of type, you can replace the type part with a placeholder name T and combine them into one definition.

#include <iostream>

template <typename T>  // T is a type name chosen at the point of use
T square(T x) {        // same T for the argument and the return value
  return x * x;
}

int main() {
  double a = square<double>(1.5);  // T is specified as double: 2.25
  int b = square<int>(3);          // T is specified as int: 9
  auto c = square(2.5f);           // T is deduced as float from the argument 2.5f: 6.25
  std::cout << a << " " << b << " " << c << "\n";
}

2.25 9 6.25

A function defined this way is called a function template.

  • template <typename T> is a declaration that “in the following definition, T is used as the name of a type.” Writing class instead of typename means the same. Even if you write class, T is not limited to class types.

  • The double written in the angle brackets of square<double> is called a template argument. (1.5) is a value argument, and has a different role.

  • If you omit the angle brackets, as in square(2.5f), the compiler determines T from the argument’s type. This is called template argument deduction.

When square<double> is used, the compiler creates a function with T replaced by double. This is called instantiation. The type is not swapped at run time; instead, a function for each type is prepared at compile time. square<int> and square<double> are separate functions with similar names.

Deduction requires that the type be determined uniquely. The following max_of uses the same T for both arguments, so passing a mix of int and double is a compile error.

template <typename T> T max_of(T a, T b) { return a > b ? a : b; }

double x = max_of(1, 2.5);          // error: T is ambiguous (int or double)
double y = max_of<double>(1, 2.5);  // with T = double given explicitly, 1 is converted to 1.0

g++ gives an explanation like the following.

error: no matching function for call to 'max_of(int, double)'
note:   deduced conflicting types for parameter 'T' ('int' and 'double')

Put the definition of a template in the header. The compiler needs to see the function body at instantiation, so for a template you write not only the declaration but also the body in a place visible to the user. MISO is header-only (a library of only .hpp files) because most of it is templates. In Chapter 2, we explained that ordinary functions placed in a header need inline. For function templates, the same definition is allowed to appear in multiple .cpp files, so no link problem occurs without inline. In MISO, some templates also have inline. Adding it does not change the meaning.

A MISO example: fourth-order central difference

space_centered_4th in miso/include/miso/mhd_integrator.hpp is a function template in which Real is decided later. We show only the declaration part (the body is read in Chapter 11).

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

  • __host__ __device__ is a specification that makes this function callable from both the CPU and the GPU (Chapter 9). In CPU-only builds, it is replaced with an empty string.

  • The caller does not write <Real>, as in space_centered_4th(qq_view, dxi, i, j, k, 1, 0, 0). If the type of qq_view is Array3DView<const float>, Real is deduced to be float. A Real inside the angle brackets of another template, as in Array3DView<const Real>, also serves as a clue for deduction.

6.2 Class Templates

Classes (structs) also have a mechanism for specifying types later. The following is a class template that bundles two values of the same type and a member function that adds them.

#include <iostream>

template <typename T>
struct ValuePair {
  T first;   // same type T for both members
  T second;
  T sum() const { return first + second; }
};

int main() {
  ValuePair<double> reals{1.5, 2.5};  // instance with two double members
  ValuePair<int> integers{2, 3};      // instance with two int members
  std::cout << reals.sum() << " " << integers.sum() << "\n";  // 4 5
}

ValuePair<double> as a whole is the name of one type. ValuePair<int> is a different type, and they cannot be assigned to each other. Both are created from the same definition, but ValuePair alone is not a type (C++17 also has a mechanism that deduces what is inside the angle brackets from the initial values, but this course does not use it).

The angle brackets of std::vector<double> and Array3D<Real> in Chapter 4 are also template arguments of this kind.

Default values for template arguments

As with function arguments, template arguments can have default values. The declaration in miso/include/miso/array3d.hpp is as follows.

/// @brief 3D Array in general execution/memory space.
template <typename T, typename Backend = backend::Host> class Array3D;

T is the element type, and Backend is where the array is placed. Writing Array3D<float> means Array3D<float, backend::Host>, and it is placed in CPU memory. To place it in GPU memory, write the second argument: Array3D<float, backend::CUDA>. Fields (mhd_fields.hpp), Grid (grid.hpp), and others also set the default of Backend to backend::Host in the same way.

A MISO example: classes that take Real

Many MISO classes are class templates that take Real. Reduced to just the member variables and the constructor, the ideal-gas equation of state in miso/include/miso/eos.hpp looks like this.

template <typename Real> struct IdealEOS {
  /// @brief Ratio of specific heats (gamma)
  Real gm;

  /// @brief Construct the EOS from the configuration file
  explicit IdealEOS(const Config &config) : gm(config["eos"]["gm"].as<Real>()) {}
  // (rest omitted)
};

The demo declares eos::IdealEOS<Real> eos;. If Real is float, gm becomes a float member. The as<Real> in config["eos"]["gm"].as<Real>() is also a call of a function template that extracts a value from the configuration file as type Real (a yaml-cpp feature).

ModelBase from Chapter 5 is a class template with three template arguments.

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

(Excerpt from miso/include/miso/mhd_model_base.hpp)

The received Real and Backend are passed as they are to the template arguments of Time<Real> and MHD<Real, Backend>. Once the demo decides Real and Backend, they propagate in turn to the member classes, and the types of all of MISO are determined. Here only grid is backend::Host regardless of Backend, because the grid coordinates are first prepared on the CPU side.

6.3 Variable Templates

Sometimes you want to use constants such as \(\pi\) and \(1/12\) at the precision of Real. If you mix a double constant into a float computation, the whole computation is done in double and then converted back to float, which slows things down or produces warnings about type mismatch. To prepare a constant for each type, use templates for variables as well. This is called a variable template.

#include <cstdio>

template <typename T>
inline constexpr T third = static_cast<T>(1.0 / 3.0);  // constant: 1/3 converted to T

int main() {
  float a = third<float>;    // 1/3 rounded to float
  double b = third<double>;  // 1/3 as double
  std::printf("%.17g\n%.17g\n", a, b);
}

0.3333333432674408
0.33333333333333331

third<float> and third<double> are separate constants, each holding a value at the precision of its type.

  • constexpr indicates a constant whose value is determined at compile time (Chapter 2).

  • static_cast<T>(...) converts 1.0 / 3.0, computed in double, to T.

  • inline indicates that this definition is placed in a header and may be included from multiple .cpp files (it can be applied to variables from C++17).

A MISO example: constants.hpp

This is the whole of miso/include/miso/constants.hpp.

#pragma once

#include <cmath>

namespace miso {

template <typename Real>
inline constexpr Real inv12 = static_cast<Real>(1.0 / 12.0);

template <typename Real> inline constexpr Real pi = static_cast<Real>(M_PI);

template <typename Real>
inline constexpr Real pii4 = static_cast<Real>(1.0 / (4.0 * M_PI));  // 1 / 4π

template <typename Real>
inline constexpr Real pii8 = static_cast<Real>(1.0 / (8.0 * M_PI));  // 1 / 8π

}  // namespace miso

The user writes pi<Real> and pii4<Real>. For example, the initial condition in demo/mhd2d_vortex/src/main.cpp has util::sin(2.0 * pi<Real> * grid.y[j]), and the differences in mhd_integrator.hpp multiply by inv12<Real>. If you print miso::pii4<float> and miso::pii4<double>, the precision differs as follows.

0.07957746833562851
0.079577471545947673

M_PI is not a constant of the C++ standard, but a macro that <cmath> provides in many environments. C++20 has a variable template with the same role, std::numbers::pi_v<T>, but MISO is C++17, so it defines its own.

6.4 Selecting Overloads with Tag Types

Consider a situation where you want to select CPU or GPU processing at compile time while keeping the same call form. In the overloading of Chapter 2, the function called was determined by the argument types. So we prepare empty types that only express “which one to compute with” and pass them as arguments.

#include <iostream>

struct Cpu {};  // type with no members (tag type)
struct Gpu {};

void run(Cpu, int n) { std::cout << "Cpu version: n = " << n << "\n"; }
void run(Gpu, int n) { std::cout << "Gpu version: n = " << n << "\n"; }

template <typename Backend>
void compute(Backend b, int n) {
  run(b, n);  // the type of b (Cpu or Gpu) selects which run is called
}

int main() {
  compute(Cpu{}, 10);  // Backend is deduced as Cpu
  compute(Gpu{}, 20);  // Backend is deduced as Gpu
}

Cpu version: n = 10
Gpu version: n = 20

  • Cpu{} is an expression that creates one value of type Cpu with no initial value. It has no contents, so the value carries no meaning, and only the type is information.

  • As in void run(Cpu, int n), you can omit the name of an unused argument. It is an argument only for distinguishing types.

  • Which run the run(b, n) inside compute calls is determined at compile time. It is not branching with if at run time.

An empty type used as a selection marker in this way is called a tag type, and the style of selecting an overload by a tag-type argument is called tag dispatch.

A MISO example: backend.hpp and for_each

MISO’s tag types are the two in miso/include/miso/backend.hpp.

namespace miso {
namespace backend {

/// @brief Tag type representing CPU (host) execution backend.
/// @details Used for compile-time backend dispatch and specialization.
struct Host {};

/// @brief Tag type representing CUDA (GPU device) execution backend.
/// @details Used for compile-time backend dispatch and specialization.
struct CUDA {};

}  // namespace backend
}  // namespace miso

The parallel loop for_each (miso/include/miso/execution.hpp) selects CPU or GPU by the tag type of its first argument. We omit the body and show only the form of the declarations.

/// @brief 1D parallel (concurrent) for loop in the Host backend.
template <class F> inline void for_each(backend::Host, Range1D range, F f) {
  // CPU: plain for loop (Chapter 7)
}

#ifdef __CUDACC__
/// @brief 1D parallel (concurrent) for loop in the CUDA backend.
template <class F> inline void for_each(backend::CUDA, Range1D range, F f) {
  // GPU: launch a kernel (Chapter 9)
}
#endif

(They are replaced with explanatory comments. The versions for the three-dimensional Range3D are lined up in the same form.)

  • The first argument is an unnamed tag type. It has the same form as run(Cpu, int) above.

  • F is the type of the lambda expression (an unnamed function) passed as the loop body. Lambda expressions are explained in Chapter 7.

  • The CUDA version is inside #ifdef __CUDACC__ and is defined only when building with the CUDA compiler (nvcc). If you pass backend::CUDA{} in a CPU-only build, there is no matching for_each, so it is a compile error.

An example that receives a tag and passes it down is the member function of IdealEOS in miso/include/miso/eos.hpp (the body is read in Chapter 7).

  template <typename Backend>
  void gas_pressure(Backend btag, mhd::FieldsView<const Real> qq,
                    Array3DView<Real> pr) const {
    // (omitted)
    for_each(btag, range, /* lambda expression */);
  }

gas_pressure is a member function that is made a template inside a class template (a member function template). The CPU Integrator (miso/include/miso/mhd_integrator_host.hpp) calls eos.gas_pressure(backend::Host{}, qq_argm.const_view(), pr.view());, and the GPU one passes backend::CUDA{}. Passing backend::Host{} makes Backend deduced as backend::Host, and the CPU version of the inner for_each is selected as well. The equation-of-state side needs to be written only once, without distinguishing CPU and GPU.

A tag value can be created not only as backend::Host{} but also as Backend{} from the type name Backend. ReduceWorkspace in miso/include/miso/mhd_resistivity.hpp has the call miso::reduce(Backend{}, range, init, f, op).

6.5 Template Specialization

Function templates and class templates use the same definition for every type. If you want to use a different definition only for a specific type, you add a definition for that type. This is called specialization.

6.5.1 Explicit specialization

A definition in which all template arguments are decided is called an explicit specialization (also called a full specialization). It is written with an empty angle bracket, template <>.

#include <iostream>

template <typename T> const char *type_name();  // declaration only; no definition
template <> const char *type_name<float>() { return "float"; }    // when T = float
template <> const char *type_name<double>() { return "double"; }  // when T = double

int main() {
  std::cout << type_name<float>() << " " << type_name<double>() << "\n";  // float double
}

There are no arguments, so deduction is impossible, and you call it by specifying the type explicitly, as in type_name<float>(). Calling type_name<int>() gives an error at link time because there is no definition with a body.

In MISO, miso/include/miso/mpi_util.hpp uses this form to return the MPI data type corresponding to Real.

/// @brief MPI Datatype corresponding to Real type
template <typename Real> MPI_Datatype data_type();
template <> inline MPI_Datatype data_type<float>() { return MPI_FLOAT; }
template <> inline MPI_Datatype data_type<double>() { return MPI_DOUBLE; }

If you write mpi::data_type<Real>(), MPI_FLOAT is returned when Real is float (Chapter 8). An explicit specialization is no longer a template but an ordinary function, so it needs inline to be placed in a header.

6.5.2 Partial specialization

For class templates, you can also write a definition that fixes only some of the template arguments. This is called partial specialization. This is an example where the element type T is left free and the contents change for each placement Backend.

#include <iostream>

struct Cpu {};
struct Gpu {};

template <typename T, typename Backend> struct Storage;  // declaration only (primary template)

template <typename T> struct Storage<T, Cpu> {  // definition for Backend = Cpu
  void describe() const { std::cout << "Placed in CPU memory\n"; }
};

template <typename T> struct Storage<T, Gpu> {  // definition for Backend = Gpu
  void describe() const { std::cout << "Placed in GPU memory\n"; }
};

int main() {
  Storage<float, Cpu> a;   // T = float, uses the Cpu definition
  Storage<double, Gpu> b;  // T = double, uses the Gpu definition
  a.describe();
  b.describe();
}

Placed in CPU memory
Placed in GPU memory

  • The first template <typename T, typename Backend> struct Storage; is called the primary template. Here it is only a declaration, with no body.

  • template <typename T> struct Storage<T, Cpu> is the definition for when “T is free and the second is Cpu.” In the angle brackets after template, write only the T that is not yet decided.

  • When you write Storage<float, Cpu>, the compiler looks for a matching partial specialization and uses it. If none matches, it uses the primary template. In this example, the primary template has no body, so Storage<float, int> is a compile error.

For each partial specialization, you can freely change both member variables and member functions. It is normal to keep the names and arguments the same so that the user can call member functions of the same name.

Function templates have no partial specialization. That is why the CPU and GPU versions of for_each are written separately with the tag-type overloads of Section 6.4 rather than with specialization.

A MISO example: Array3D and Integrator

After the declaration seen in Section 6.2, miso/include/miso/array3d.hpp defines a partial specialization for each Backend. We show them with the bodies omitted.

/// @brief 3D Array in general execution/memory space.
template <typename T, typename Backend = backend::Host> class Array3D;

/// @brief 3D Array in host memory.
template <typename T> class Array3D<T, backend::Host> {
  // allocate CPU memory with new T[...] (Chapter 4)
};

#ifdef __CUDACC__
/// @brief 3D Array in CUDA device memory.
template <typename T> class Array3D<T, backend::CUDA> {
  // allocate GPU memory with cudaMalloc (Chapter 9)
};
#endif  // __CUDACC__

(They are replaced with explanatory comments.)

This is the same form as Storage above. Array3D<float, backend::Host> and Array3D<float, backend::CUDA> both have member functions of the same names, such as view(), extent(), and copy_from(), but differ in how memory is allocated and released. Array1D, Array2D, Array4D, Fields, and Grid are built the same way.

Integrator, which handles time integration, also has only a declaration as its primary template (miso/include/miso/mhd_integrator.hpp).

template <typename Real, typename Backend> struct Integrator;

The CPU definition is in miso/include/miso/mhd_integrator_host.hpp, and the GPU one is in miso/include/miso/mhd_integrator_cuda.cuh.

template <typename Real> struct Integrator<Real, backend::Host> {
  // (omitted)
};

template <typename Real> struct Integrator<Real, backend::CUDA> {
  // (omitted)
};

At the end of mhd_integrator.hpp, the CPU one is always included, and the GPU one only when __CUDACC__ is defined.

#include "mhd_integrator_host.hpp"

#ifdef __CUDACC__
#include "mhd_integrator_cuda.cuh"
#endif

MHD<Real, Backend> (miso/include/miso/mhd_solver.hpp), which uses Integrator, has Integrator<Real, Backend> integrator; as a member and just calls integrator.update(...) and integrator.cfl(...). Once Backend is decided, which definition is used is decided too. The two definitions are aligned so that both have update, cfl, and apply_boundary_condition with the same names and argument lists, and only the GPU one has members such as the CUDA kernel shape (cuda::KernelShape3D) and the workspace for reduction.

6.6 if constexpr and static_assert

if constexpr

Template arguments can be not only types but also bool and int values. The following is a function that specifies, by a template argument, whether to flip the sign.

template <bool Negate>
double change_sign(double value) {
  if constexpr (Negate) {  // Negate is known at compile time
    return -value;
  } else {
    return value;
  }
}

change_sign<true>(3.0) returns −3, and change_sign<false>(3.0) returns 3. if constexpr is a branch used when the condition is determined at compile time. At instantiation, the side not chosen is not built into the function. Only return -value; remains in change_sign<true>. The ordinary if of Chapter 2 can take a condition determined at run time, but both sides are compiled. The condition of if constexpr can only be an expression whose value is determined at compile time.

To change the processing according to the type, use std::is_same_v<A, B> from <type_traits>. It is a compile-time constant that is true if A and B are the same type and false otherwise.

#include <type_traits>

template <typename T>
const char *precision_name() {
  if constexpr (std::is_same_v<T, float>) {
    return "single precision";
  } else {
    return "double precision or higher";
  }
}

precision_name<float>() returns "single precision", and precision_name<double>() returns "double precision or higher".

That the side not chosen is not built in has another implication. Processing that cannot be written for some template arguments (such as calling a member the type does not have) does not cause an error if you put it on the side that is not chosen. The MISO example in Section 6.7 uses this property.

static_assert

When you want to put a condition on a template argument, use static_assert (a compile-time check).

#include <type_traits>

template <typename T>
T half(T x) {
  static_assert(std::is_floating_point_v<T>, "T must be float or double");
  return x / 2;
}

static_assert(condition, "message"); makes it a compile error and prints the message if the condition is false. half(3.0) compiles as is, but half(3) makes T an int, so it stops with the following error.

error: static assertion failed: T must be float or double

std::is_floating_point_v<T> is a constant that is true if T is a floating-point type such as float or double. Unlike the run-time assert, static_assert does not disappear in Release builds. It can find incorrect usage before the program runs.

You can check the details of the grammar in the C++17 draft standard N4659: template parameters, if constexpr, and type properties (type traits).

A MISO example: branching in the reduction kernel

In the GPU reduction (reduce) in miso/include/miso/execution.hpp, how the indices are extracted changes depending on whether the range is one-dimensional or three-dimensional. Range is a template argument, and takes either Range1D, which represents a one-dimensional range, or Range3D, which represents a three-dimensional range (Chapter 7).

    if constexpr (std::is_same_v<Range, Range1D>) {
      const int i = range.unravel(idx);
      v = map(i);
    } else {
      Index3D p = range.unravel(idx);
      v = map(p.i0, p.i1, p.i2);
    }

unravel of Range1D returns an int, and unravel of Range3D returns an Index3D (a triple of indices). Written with an ordinary if, Index3D p = range.unravel(idx); would also be compiled for Range1D, and the types would not match, causing an error. With if constexpr, only the side that fits the type is built in.

6.7 Mechanisms for Examining the Properties of Types

Users of MISO prepare the members eos, ic, bc, and src in a model that inherits from ModelBase (Chapter 5). In the source term src, you write only the terms you need (for example, vy for gravity), and terms you did not write are treated as 0. To achieve this, you must check at compile time “whether this type has a member function vy,” and branch with if constexpr based on the result. This section explains the mechanism just enough for you to read it.

6.7.1 Structs that hold a type’s property as a value

<type_traits> has two structs that represent true or false.

std::true_type::value   // true
std::false_type::value  // false

value is a constant belonging to the class (static constexpr bool). It can be read as type_name::value without creating an instance. If a struct inherits from one of these two, that struct’s ::value also becomes true or false.

std::is_same_v<A, B> in Section 6.6 is a short way of writing std::is_same<A, B>::value. std::is_same<A, B> is built to inherit std::true_type if A and B are the same and std::false_type otherwise. A class template that returns the result of examining a type through ::value in this way is called a type trait.

6.7.2 Checking for members with std::void_t

We build the type trait has_gm, which checks “whether type T has a member gm,” from two definitions.

#include <iostream>
#include <type_traits>

// (1) primary template: false if nothing else matches
template <class T, class = void>
struct has_gm : std::false_type {};

// (2) partial specialization: matches only when T::gm is valid; true
template <class T>
struct has_gm<T, std::void_t<decltype(&T::gm)>> : std::true_type {};

struct Gas {
  double gm;  // has a member gm
};
struct Dust {
  double size;  // has no gm
};

int main() {
  std::cout << has_gm<Gas>::value << " " << has_gm<Dust>::value << "\n";  // 1 0
}

We read the parts in order.

  • class = void is an unnamed second template argument whose default is void. The user writes only one, has_gm<Gas>, and the second becomes void.

  • &T::gm is an expression meaning “point to the member gm of T.” If T has no gm, this expression cannot be written.

  • decltype(expression) finds the type of the expression without evaluating it. Here only whether the expression can be written matters.

  • std::void_t<...> is a type that becomes void regardless of its contents, as long as what is inside the angle brackets can be written correctly.

For has_gm<Gas>, the compiler reasons as follows.

  1. has_gm<Gas>, with the default filled in, is has_gm<Gas, void>.

  2. Applying T = Gas to partial specialization (2), std::void_t<decltype(&Gas::gm)> becomes void, and (2) becomes the definition for has_gm<Gas, void>.

  3. 1 and 2 match, so (2) is used and value is true.

For has_gm<Dust>, &Dust::gm cannot be written in 2. In this case, the compiler does not report an error and merely removes (2) from the candidates. The remaining (1) is used, and value is false. The C++ rule that when applying template arguments produces a form that cannot be written, the candidate is simply removed without an error is called SFINAE (Substitution Failure Is Not An Error). The rule is in template argument deduction in the C++17 draft standard.

To check whether a member function can be called with a given argument list, write a call expression instead of &T::name.

#include <utility>  // std::declval

template <class T, class = void> struct has_dt_limit : std::false_type {};
template <class T>
struct has_dt_limit<T, std::void_t<decltype(std::declval<const T &>().dt_limit())>>
    : std::true_type {};

std::declval<const T &>() is a hypothetical expression meaning “if there were a value of type const T &,” and can be used only inside decltype. It lets you check whether t.dt_limit() can be written without creating an instance of T. Because it is called through const T &, the result is false unless dt_limit is a const member function.

struct A { double dt_limit() const { return 0.1; } };  // has_dt_limit<A>::value is true
struct B { double dt_limit() { return 0.1; } };        // false: not const
struct C {};                                           // false: no dt_limit

A MISO example: has_eos and has_dt_limit

The beginning of miso/include/miso/mhd_model_base.hpp has type traits in the same form as above (excerpt).

// SFINAE to check if the derived class has required member variables.
template <class T, class = void> struct has_eos : std::false_type {};
template <class T>
struct has_eos<T, std::void_t<decltype(&T::eos)>> : std::true_type {};

// Optional time step limit of the source term: Real dt_limit() const
template <class T, class = void> struct has_dt_limit : std::false_type {};
template <class T>
struct has_dt_limit<T,
                    std::void_t<decltype(std::declval<const T &>().dt_limit())>>
    : std::true_type {};
template <class T, class = void> struct has_member_dt_limit : std::false_type {};
template <class T>
struct has_member_dt_limit<T, std::void_t<decltype(&T::dt_limit)>>
    : std::true_type {};

has_ic, has_bc, and has_src have the same form as has_eos. ModelBase::run() checks them with static_assert (excerpt).

    static_assert(has_eos<Derived>::value,
                  "Derived must have member variable 'eos'.");

If you forget to write eos in the model, compilation stops with this message before the program runs.

ModelBase::update() uses the source term’s upper limit on the time step, dt_limit(), only when the source term has it (excerpt; the MPI part is omitted).

  void update() {
    auto &d = derived();
    using Source = decltype(d.src);
    static_assert(!has_member_dt_limit<Source>::value ||
                      has_dt_limit<Source>::value,
                  "Source::dt_limit must be callable as dt_limit() const");
    auto dt = mhd.cfl(d.eos);
    if constexpr (has_dt_limit<Source>::value) {
      Real dt_src = static_cast<Real>(d.src.dt_limit());
      // (omitted: compute dt_src_g, the minimum over all ranks)
      dt = std::min(dt, dt_src_g);
    }
    mhd.update(dt, d.eos, d.bc, d.src);
    time.update(dt);
  }

  • using Source = decltype(d.src); gives the type of the member src the alias Source.

  • d.src.dt_limit() inside if constexpr is an expression that cannot be written for a source term without dt_limit. As seen in Section 6.6, it is not built in when the condition is false, so there is no error.

  • The condition of the static_assert is “there is no name dt_limit, or it can be called as dt_limit() const” (! is negation and || is “or”). If you wrote dt_limit but forgot const, checking with has_dt_limit alone would give false and it would be silently ignored. By checking separately for the presence of the name (has_member_dt_limit) and for callability (has_dt_limit), the code turns such a mistake into a compile error.

6.7.3 Combining macros and templates (source terms)

There are eight source terms, ro, vx, vy, vz, bx, by, bz, and ei, and each can be omitted. If you wrote the two kinds of type traits above and the function that fills in 0 for all eight, the same definition would be repeated eight times, differing only in the name. In miso/include/miso/mhd_integrator.hpp, this is written only once using a preprocessor macro (Chapter 2).

#define MISO_MHD_DEFINE_SOURCE_TERM(name)                                      \
  namespace impl_source {                                                      \
  template <class Source, class = void>                                        \
  struct has_member_##name : std::false_type {};                               \
  template <class Source>                                                      \
  struct has_member_##name<Source, std::void_t<decltype(&Source::name)>>       \
      : std::true_type {};                                                     \
  template <class Source, class Real, class = void>                            \
  struct has_##name : std::false_type {};                                      \
  template <class Source, class Real>                                          \
  struct has_##name<Source, Real,                                              \
                    std::void_t<decltype(std::declval<const Source &>().name(  \
                        std::declval<FieldsView<const Real>>(), 0, 0, 0))>>    \
      : std::true_type {};                                                     \
  }                                                                            \
  /** @brief Source term of `name` (0 if the Source does not define it) */     \
  template <typename Real, typename Source>                                    \
  __host__ __device__ inline Real source_##name(                               \
      const Source &src, const FieldsView<const Real> &qq, int i, int j,       \
      int k) {                                                                 \
    static_assert(!impl_source::has_member_##name<Source>::value ||            \
                      impl_source::has_##name<Source, Real>::value,            \
                  "Source::" #name " must be callable as " #name               \
                  "(FieldsView<const Real>, int, int, int) const");            \
    if constexpr (impl_source::has_##name<Source, Real>::value) {              \
      return src.name(qq, i, j, k);                                            \
    } else {                                                                   \
      return Real(0);                                                          \
    }                                                                          \
  }
MISO_MHD_DEFINE_SOURCE_TERM(ro)
MISO_MHD_DEFINE_SOURCE_TERM(vx)
// (same for vy, vz, bx, by, bz, ei)
#undef MISO_MHD_DEFINE_SOURCE_TERM

We read the macro symbols first.

Symbol

Meaning

\ at the end of a line

Continues the macro definition onto the next line. A macro is originally written on one line

name

The macro’s argument. For MISO_MHD_DEFINE_SOURCE_TERM(vx), every name is replaced with vx

has_##name

## joins what is before and after into one name. If name is vx, it becomes has_vx

#name

# turns the argument into the string "vx"

#undef

Cancels the definition of a macro that is no longer needed

"Source::" #name " must be ..." becomes strings placed side by side, like "Source::" "vx" " must be ...". In C++, adjacent strings are joined into one, so the message becomes "Source::vx must be callable as vx(...) const".

Expanding MISO_MHD_DEFINE_SOURCE_TERM(vx) by hand and formatting it for readability gives the following (it is expanded for explanation, and is not written in this form in the source).

namespace impl_source {
// has a member named vx?
template <class Source, class = void>
struct has_member_vx : std::false_type {};
template <class Source>
struct has_member_vx<Source, std::void_t<decltype(&Source::vx)>>
    : std::true_type {};
// can src.vx(qq, i, j, k) be called? (src is const, qq is FieldsView<const Real>)
template <class Source, class Real, class = void>
struct has_vx : std::false_type {};
template <class Source, class Real>
struct has_vx<Source, Real,
              std::void_t<decltype(std::declval<const Source &>().vx(
                  std::declval<FieldsView<const Real>>(), 0, 0, 0))>>
    : std::true_type {};
}  // namespace impl_source

/** @brief Source term of `vx` (0 if the Source does not define it) */
template <typename Real, typename Source>
__host__ __device__ inline Real source_vx(const Source &src,
                                          const FieldsView<const Real> &qq,
                                          int i, int j, int k) {
  static_assert(!impl_source::has_member_vx<Source>::value ||
                    impl_source::has_vx<Source, Real>::value,
                "Source::vx must be callable as vx(FieldsView<const Real>, int, int, int) const");
  if constexpr (impl_source::has_vx<Source, Real>::value) {
    return src.vx(qq, i, j, k);  // call vx if it exists
  } else {
    return Real(0);              // otherwise 0
  }
}

It is the has_member_dt_limit and has_dt_limit of Section 6.7.2, and the static_assert and if constexpr of update(), combined for vx. has_vx also takes Real in order to determine the argument type FieldsView<const Real>.

In the time integration (update_sc4 in miso/include/miso/mhd_integrator_host.hpp), source_vx(src, c_qq, i, j, k) is added to the right-hand side of the equation of motion. <Real> is not written because Real is deduced from the type of c_qq, FieldsView<const Real> (Section 6.1). If src has no vx, source_vx becomes a function that merely returns Real(0), and no extra computation remains.

The source term in demo/hd2d_rayleigh_taylor/src/main.cpp defines four terms: vx, vy, vz, and ei (only vy is excerpted).

struct SourceTerm {
  Real g_grav;
  // (omitted)

  /// External force: y-direction
  __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);
  }
  // (rest omitted)
};

If you write vy without const, has_member_vy is true and has_vy is false, and compilation stops with the following message.

error: static assertion failed: Source::vy must be callable as vy(FieldsView<const Real>, int, int, int) const
::

In C++20, you can write “whether this expression can be written” directly with a requires expression. The same check as has_dt_limit can be written as follows (-std=c++20 is required).

template <class T>
concept HasDtLimit = requires(const T &t) { t.dt_limit(); };

HasDtLimit<T> becomes true or false, and can be used as in if constexpr (HasDtLimit<Source>). The method of Section 6.7.2 that combines std::void_t with partial specialization is a technique for doing the same thing up to C++17. MISO is written in C++17, so it uses the std::void_t method. The rules are in constraints and requires expressions in the C++20 draft standard N4861.

6.7.4 .template inside templates

ModelBase::run() has the following line (miso/include/miso/mhd_model_base.hpp).

    if (config["io"]["continue"].template as<bool>() &&

In IdealEOS in Section 6.2, we wrote config["eos"]["gm"].as<Real>(), but here there is a template between . and as. This is a mark that tells the compiler “the name that follows is a function template.” It is needed in cases like the following.

#include <iostream>

struct Reader {
  template <class U> U get() const { return U(2.5); }  // member function template
};

template <class T>
double read_twice(const T &r) {
  // T, the type of r, is a template argument, so whether get is a template is unknown until instantiation
  return r.template get<double>() * 2;
}

int main() {
  Reader r;
  std::cout << r.get<int>() << " " << read_twice(r) << "\n";  // 2 5
}

The type of r in main is determined to be Reader, so you can write r.get<int>() as it is. The type of r in read_twice is T, and what goes into T is not known until instantiation. In this case, the compiler reads get as a non-template name and tries to interpret r.get < double > () as comparisons, which causes an error. With .template, it reads get<double> as a template call. typename, which tells the compiler that something is a type name, as in typename T::value_type, is used for the same reason.

The config in ModelBase has type Config &, whose type is determined regardless of the template arguments, so this line compiles even if .template is omitted. Adding it where it is not needed does not change the meaning (the same style also appears in miso/include/miso/mhd_resistivity.hpp). When reading, it is enough to recognize it as “a mark for a template call,” and the compiler tells you with an error where you need to write it. The rules are in name resolution in the C++17 draft standard.

Exercises: Specifying and Selecting Templates

  1. For max_of in Section 6.1, give the return type and value of max_of(2, 7) and max_of<double>(1, 2.5).

  2. Write a variable template two_pi<T> that represents \(2\pi\), and print two_pi<float> and two_pi<double> with %.17g.

  3. Add struct Fpga {}; to Cpu, Gpu, run, and compute of Section 6.4, check what happens when you call compute(Fpga{}, 1), and explain why.

  4. For the following types, answer which definition in which file under miso/include/miso/ is used.

    • Array3D<float>

    • Integrator<double, backend::Host>

    • Array3D<float, backend::CUDA> in a CPU-only build

  5. Following has_gm in Section 6.7.2, write has_size<T>, which checks whether there is a member size. Include miso/include/miso/execution.hpp, print has_size<miso::Range1D>::value and has_size<miso::Index3D>::value, and compare the results with the definitions in execution.hpp.

  6. Read SourceTerm in demo/hd2d_rayleigh_taylor/src/main.cpp, and, among the eight functions source_ro through source_ei, separate those that call a member of SourceTerm from those that return Real(0). Also answer whether the inside of if constexpr (has_dt_limit<Source>::value) in ModelBase::update() is built in.

::

1. For max_of(2, 7), both arguments are int, so T is deduced to be int and it returns the int 7. For max_of<double>(1, 2.5), T is specified as double, so 1 is converted to 1.0 and it returns the double 2.5. Even for a call that would be an error because T cannot be determined without a specification, if you specify it explicitly, ordinary type conversion is used.

2. Example:

#include <cstdio>

template <typename T>
inline constexpr T two_pi = static_cast<T>(2.0 * 3.14159265358979323846);

int main() { std::printf("%.17g\n%.17g\n", two_pi<float>, two_pi<double>); }

6.2831854820251465
6.2831853071795862

float has about 7 significant digits, so the values differ from the 8th digit. Following constants.hpp, you may also use M_PI and write static_cast<T>(2.0 * M_PI).

3. It is a compile error. In g++, no matching function for call to 'run(Fpga&, int&)' is shown. Instantiating compute<Fpga> makes b in run(b, n) have type Fpga, but no run that takes an Fpga is defined. Selecting overloads by tag types is done at compile time, so a missing corresponding function is found before running. The same applies when for_each(backend::CUDA{}, ...) is called in a CPU-only build of MISO.

4.

  • Array3D<float> becomes Array3D<float, backend::Host> by default, and the partial specialization template <typename T> class Array3D<T, backend::Host> in array3d.hpp is used.

  • For Integrator<double, backend::Host>, template <typename Real> struct Integrator<Real, backend::Host> in mhd_integrator_host.hpp is used (Real is double).

  • In a CPU-only build, the partial specialization of Array3D<T, backend::CUDA> is inside #ifdef __CUDACC__, so it is not defined. The primary template is only a declaration, so writing Array3D<float, backend::CUDA> b(2, 2, 2); gives the g++ compile error has initializer but incomplete type (a type without a definition).

5. Example:

#include <iostream>
#include <type_traits>
#include <miso/execution.hpp>

template <class T, class = void> struct has_size : std::false_type {};
template <class T>
struct has_size<T, std::void_t<decltype(&T::size)>> : std::true_type {};

int main() {
  std::cout << has_size<miso::Range1D>::value << " "
            << has_size<miso::Index3D>::value << "\n";
}

Compile with the location of MISO’s headers specified by -I/path/to/MISO/miso/include. The output is 1 0. Range1D in execution.hpp has the member function size(), and Index3D has only int i0, i1, i2;, which matches the result. &T::size can be used for both member variables and member functions.

6. SourceTerm defines vx, vy, vz, and ei, all in the form (mhd::FieldsView<const Real>, int, int, int) const.

Function

Processing used

source_vx, source_vy, source_vz, source_ei

Calls src.vx(qq, i, j, k) and so on

source_ro, source_bx, source_by, source_bz

Returns Real(0)

vx, vz, and ei just return 0, but they are defined, so they are called (the result is the same if the definitions are omitted). SourceTerm has no dt_limit, so has_dt_limit<SourceTerm>::value is false. The inside of if constexpr is not built in, and the time step is determined only by the CFL condition (mhd.cfl(d.eos)).