2. Variables, Functions, and Control Flow

In this chapter, you will learn to write short sequential computations in C++ using variables, operators, conditional branches, loops, and functions. Next, you will learn how to split code into files and namespaces, input and output in the standard library, strings and variable-length arrays, exceptions, and the basics of the preprocessor. In each section, a small illustrative example is followed by real MISO code that uses the same syntax.

References and pointers are covered in Chapter 3, classes and structs in Chapter 4, and templates in Chapter 6. When these appear in an excerpt from MISO, the text points out which parts you can skip.

1. Program structure

The following is a complete illustrative program. As in Chapter 1, you can build it with g++ -std=c++17 hello.cpp -o hello and run it with ./hello.

#include <iostream>

int main() {
  double dt = 0.1;  // time step
  int n_step = 3;   // number of steps
  std::cout << "time = " << dt * n_step << std::endl;
}

time = 0.3

A program starts from a function called main and executes the statements inside its braces { ... } from top to bottom. This corresponds to the part between program and end program in Fortran.

#include <iostream> makes std::cout, which is used for output to the screen, available (Section 6).
A comment runs from // to the end of the line. It corresponds to ! in Fortran. You can also enclose multiple lines with /* ... */.
End each statement with ;. Unlike in Fortran, a newline does not end a statement, so you can write one statement across multiple lines.
std::cout << value prints the value to the screen. std::endl is a newline (Section 8).
Uppercase and lowercase letters are distinguished. dt and DT are different names.

When execution reaches the final } of main, the program terminates normally. Unless stated otherwise, the short code examples below show the part that goes inside a function such as main. Put the #include lines at the top of the file.

A real MISO example — the main function in 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();
}

int argc, char **argv in the parentheses receives the arguments given on the command line at run time (such as --config=...). How to read char ** is explained in Chapter 3. In the body, the program sets up the environment such as MPI (env), reads the configuration file (config), creates the model that represents the problem, and advances the computation with model.run(). Env, Config, and Model are classes (Chapter 4); for now it is enough to read each as “create one variable of that type.” The main of every MISO demo is about this short; the substance of the computation lives in the headers and in Model.

2. Variables and types

2.1 Type, name, and initial value

The declaration corresponding to Fortran’s integer :: nx = 48 is written as follows in C++.

int nx = 48;
double dt = 0.1;
bool io_enabled = true;

Read int nx = 48; from left to right as “type,” “variable name,” and “initial value.” Telling the compiler a variable’s name and type is called a declaration, and giving it its first value when it is created is called initialization. The later nx = 64; is an assignment to the existing variable. In C++, unlike in Fortran, you do not have to gather declarations at the top of a program unit; you can declare a variable where you start using it. There is nothing corresponding to implicit none; using an undeclared name is always a compile error.

Type	Values it handles	Usual Fortran counterpart	Declaration example
`int`	Integer (32 bits in this course’s environment)	`integer`	`int n_step = 0;`
`long long`	Integer with a wider range	`integer(8)`	`long long n_cells = 1000000;`
`float`	Single-precision real (32 bits)	`real` / `real(4)`	`float x = 0.5f;`
`double`	Double-precision real (64 bits)	`real(8)`, `double precision`	`double time = 0.0;`
`bool`	Boolean `true` / `false`	`logical`	`bool periodic = false;`
`char`	A single character	`character`	`char axis = 'x';`
`std::string`	String (`#include <string>`)	`character(len=*)`	`std::string name = "mhd";`

Numeric literals also have types. 48 is an int, 0.1 is a double, and adding f, as in 0.1f, makes it a float. Note that this is the reverse of Fortran, where 0.1 is single precision and 0.1d0 is double precision. 'x' is a single character (char) and "x" is a string, so the kind of quotation mark also matters.

The value of a variable declared in a function with just double value; is undefined. It is not guaranteed to be 0, so either give it an initial value or always assign to it before reading it.

2.2 Ways to initialize

You can give the initial value with =, parentheses, or braces.

int a = 48;
int b(48);   // parentheses also give the initial value 48
int c{48};   // braces also give the initial value 48
int zero{};  // empty braces initialize to 0

Braces turn conversions that would change the value into errors.

int n = 2.9;   // compiles; the fractional part is dropped, so n is 2
int m{2.9};    // error: narrowing conversion of '2.8999999999999999e+0' from 'double' to 'int'

All three forms appear in the MISO source. It is enough to read each as “put a value in when creating it.” For classes, parentheses and braces can mean different things; this is covered in Chapter 4.

2.3 const and constexpr

const specifies that a variable will not be modified after initialization. It is close to Fortran’s parameter, but it can also be initialized with a value determined at run time.

const int nx = 48;  // nx can no longer be modified
int ny = nx;        // a separate variable holding a copy of nx
ny = 72;            // only ny changes; nx stays 48
// nx = 64;         // error: cannot assign to a const variable

const int nx and int const nx mean the same thing.

constexpr declares a constant whose value is determined at compile time. Its initial value must also be an expression that can be evaluated at compile time.

constexpr int n_dim = 3;
constexpr int n_faces = 2 * n_dim;  // fixed to 6 at compile time

const means “not modified,” and constexpr means “determined at compile time (and therefore not modified).” constexpr is closer in meaning to Fortran’s parameter.

A real MISO example — the following line appears inside IdealEOS::gas_pressure in miso/include/miso/eos.hpp.

const Real gm_ = gm;  // device lambdas must not capture `this`

Real is an alias for the real-number type (Section 2.6), and gm is the ratio of specific heats. The program copies gm, the ratio of specific heats read at run time, into a local variable gm_ that is never modified. Why the copy is necessary is explained in Chapter 7 and Chapter 9.

A real example of constexpr is in parse_config_filepath in Section 4.5 (constexpr std::string_view key = "--config=";).

2.4 auto

auto asks the compiler to determine the type from the initial value. The type does not change while the program runs.

auto n_step = 100;       // int
auto dt = 0.1;           // double
auto x = 0.1f;           // float
const auto limit = 100;  // const int

You cannot write auto dt; without an initial value. auto is often used when a type name is long or when you receive the return value of a function as it is.

A real MISO example — a line in main in demo/mhd2d_vortex/src/main.cpp.

auto config_path = parse_config_filepath(argc, argv);

The return type of parse_config_filepath is std::optional<std::string>, a standard library type that means “either there is one string or there is nothing.” The code receives it with auto instead of writing the long type name. The config_path.value_or("./config.yaml") on the next line means “use the value if there is one, otherwise use "./config.yaml".” In other words, if you do not specify --config, the program reads config.yaml in the current directory.

2.5 Integers, floating-point numbers, and type conversion

As in Fortran, division of integers gives an integer with the fractional part truncated. The remainder is obtained with % (Fortran’s mod).

int a = 7;
int b = 2;
std::cout << a / b << "\n";                       // 3
std::cout << a % b << "\n";                       // 1
std::cout << 1 / 2 << "\n";                       // 0
std::cout << 1.0 / 2 << "\n";                     // 0.5
std::cout << static_cast<double>(a) / b << "\n";  // 3.5

In an expression that mixes integers and reals, the integer is converted to a real before the calculation (1.0 / 2 is 1.0 / 2.0). If you mix float and double, the float is converted to double and the calculation is done in double. This automatic conversion is called implicit conversion.

To change a type explicitly, write static_cast<type>(value). It corresponds to Fortran’s dble(a) or int(x). Writing the target type inside <> uses the same notation as templates (Chapter 6), but here it is enough to read it as “change the type.” The older forms double(a) and (double)a also exist, and MISO often uses the form Real(0.5) (next section).

float and double differ in the number of digits they can represent.

float x = 0.1f;
double y = 0.1;
std::cout.precision(17);           // print with 17 digits
std::cout << x << " " << y << "\n";

0.10000000149011612 0.10000000000000001

A float has about 7 significant digits, and a double has about 16.

A real MISO example — the constructor of Config in miso/include/miso/config.hpp contains the following line.

int yaml_str_length = static_cast<int>(yaml_str.length());

The string length yaml_str.length() is returned as a non-negative integer type (std::size_t). It is converted to int with static_cast<int> so that it can be passed to an MPI function.

2.6 Type aliases and Real

using new_name = type; gives a type another name. It is not a variable declaration.

using Real = float;  // Real is an alias for float
Real dt = 0.1f;      // same as float dt = 0.1f;

MISO demos define Real at the beginning of the file. The following is an excerpt from demo/mhd2d_vortex/src/main.cpp.

using Real = float;

All physical variables are declared as Real, so changing this one line to using Real = double; makes the computation run in double precision. This is the same idea as defining integer, parameter :: rk = kind(1.0d0) in Fortran and writing real(rk). (Inside the MISO library itself, Real is received as a template argument; that mechanism is explained in Chapter 6.)

When you use Real, MISO’s convention is to write constants so that they match Real.

using Real = float;
Real a = 1.0f / 3.0f;
Real h1 = a * 0.5;        // converts a to double, computes in double, then converts the result back to float
Real h2 = a * Real(0.5);  // computes in float
auto h3 = a * 0.5;        // auto, so h3 is double

0.5 is a double, so a * 0.5 is computed in double. It goes back to float when stored in h1, so the value is almost the same, but an extra conversion is inserted, which slows down double-precision computation on GPUs. If you receive the result with auto, as in h3, it becomes double even though you meant float.

A real MISO example — an excerpt from flux_core in miso/include/miso/limiter.hpp.

Real pp = (Real(0.5) + util::copysign(Real(0.5), ra)) *
          util::max2(Real(0.0), Real(1.0) + fh * (ra - Real(1.0)));

Real(0.5) is “the value 0.5 as type Real.” Every term is a Real, so the computation stays in Real. util::copysign and util::max2 are small MISO functions that, respectively, copy a sign onto a value and return the larger of two values. (template <typename Real> at the head of the functions in this file is explained in Chapter 6, and __host__ __device__ in Chapter 9. When reading, you can think of Real as float.)

2.7 Enumerations

With enum class, you can create a type that represents one of a fixed set of choices.

enum class Axis { X, Y, Z };  // type Axis with 3 choices
Axis direction = Axis::Y;     // store one of them in a variable
if (direction == Axis::Y) {
  // handle the y direction
}

Writing Axis::Y is clearer than representing a direction with the integers 0, 1, and 2, and mixing them up causes a compile error.

A real MISO example — miso/include/miso/types.hpp contains a series of enumerations that represent directions and boundary sides (excerpt).

enum class Direction { X, Y, Z };

enum class Side { INNER, OUTER };

Direction::X represents the x direction, and Side::INNER represents the boundary on the inner side of the domain (the side with smaller coordinates). They are used in functions for boundary conditions and similar code (Section 4.6).

3. Operators

Arithmetic operators are almost the same as in Fortran. The difference is that there is no exponentiation operator. You cannot write x**2, so write x * x or use a function (in MISO, util::pow2(x)).

Meaning	C++	Fortran
Four arithmetic operations	`a + b`, `a - b`, `a * b`, `a / b`	Same
Integer remainder	`a % b`	`mod(a, b)`
Exponentiation	None (`x * x`, `std::pow(x, y)`)	`x**2`
Accumulate	`a += b` (`a = a + b`). `-=`, `*=`, and `/=` work the same way	`a = a + b`
Increment and decrement by 1	`++i`, `--i`	`i = i + 1`
Equal and not equal	`a == b`, `a != b`	`a == b`, `a /= b`
Less than and greater than	`a < b`, `a <= b`, `a > b`, `a >= b`	Same
And, or, and not	`a && b`, `a \|\| b`, `!a`	`.and.`, `.or.`, `.not.`

There are three points to watch for.

= is assignment, and == is comparison. if (i = 3) is not a comparison; it assigns 3 to i.
/= means “divide and assign.” Fortran’s “not equal” is != in C++.
&& and || evaluate from the left, and if the left side alone determines the result, they do not evaluate the right side (short-circuit evaluation). In if (n > 0 && total / n > 1.0), the division is not performed when n is 0. Fortran’s .and. has no such guarantee.

Write a range check as 0 <= i && i < n. The mathematical form 0 <= i < n does not work as intended.

There is also the conditional operator condition ? value_if_true : value_if_false, which chooses between two values.

double dt = 0.2;
double limited = (dt > 0.1) ? 0.1 : dt;  // the condition is true, so 0.1

A real MISO example — is_root in miso/include/miso/env.hpp returns the result of a comparison (a bool) as it is.

inline bool is_root() { return rank() == 0; }

rank() is a function that returns the MPI rank number, and returns true for the process with rank 0 (for how to write functions, see Section 5; for inline, see Section 6).

clamp in miso/include/miso/utility.hpp nests conditional operators to keep a value v between lo and hi inclusive.

return v < lo ? lo : (v > hi ? hi : v);

If v < lo, the result is lo; otherwise it is “hi if v > hi, otherwise v.”

parse_config_filepath in miso/include/miso/config.hpp uses short-circuit evaluation.

if (arg.size() >= key.size() && arg.compare(0, key.size(), key) == 0) {

The condition on the left confirms that the argument is at least as long as --config=, and then the right side checks whether the beginning matches --config=. If the left side is false, the right side is not executed.

4. Conditional branches and loops

4.1 if and else

double dt = 0.2;
if (dt > 0.1) {   // 0.2 > 0.1, so this branch runs
  dt = 0.1;
} else {
  dt = 0.05;      // not executed in this example
}
// dt is 0.1

This corresponds to Fortran’s if (...) then ... else ... end if. Put the condition in parentheses and enclose the statements to execute in braces { ... }. The range enclosed by braces is called a block. You do not write then or end if. To add a condition, use else if (...). You can omit the braces when there is only one statement, but this course includes them so that the extent is clear.

A real MISO example — an excerpt from the constructor of Time in miso/include/miso/time.hpp.

if (!(tend > 0)) {
  throw std::runtime_error("time.tend must be positive.");
}

If time.tend (the end time) in config.yaml is not positive, the program reports an error and stops (throw is covered in Section 9). The code writes !(tend > 0) instead of tend <= 0 so that it also catches the case where the value is NaN (not a number). Every comparison with NaN is false, so tend <= 0 would miss a NaN.

4.2 for and the range of indices

int total = 0;
for (int i = 0; i < 3; ++i) {  // runs for i = 0, 1, 2
  total += i;                  // total goes 0 -> 1 -> 3
}
// total is 3. i, declared inside the for, cannot be used here

The parentheses consist of three parts separated by ;.

Part	Example	When it is executed
Initialization	`int i = 0`	Once, when entering the loop
Condition to continue	`i < 3`	Every time, before executing the block
Update	`++i`	Every time, after finishing the block

In Fortran, this is written as follows.

total = 0
do i = 0, 2
  total = total + i
end do

Fortran’s do i = 0, 2 includes the end point 2. C++’s i < 3 does not include 3. In C++, the indices of an array with n elements run from 0 to n - 1, and the basic loop is written for (int i = 0; i < n; ++i). If you simply turn Fortran’s do i = 1, n into i <= n, you are off by one and read outside the array.

You can choose the start value, condition, and increment freely.

How the `for` begins	`i` used in the block
`for (int i = 1; i <= 3; ++i)`	1, 2, 3
`for (int i = 0; i < 6; i += 2)`	0, 2, 4
`for (int i = 3; i >= 1; --i)`	3, 2, 1

The condition to continue is evaluated every time. If you change n inside the loop, the result of the i < n test changes too. Unlike Fortran’s do, the number of iterations is not fixed at the start. Also, if you write ; after the parentheses, as in for (...);, only that ; is repeated, and the block that follows is outside the loop.

::

C++ indices start at 0 because it inherits the way C arrays work. For an array a, a[i] means “the element i places from the start,” and the first element is 0 places from the start (pointer arithmetic in Chapter 3).

A short note by Dijkstra is often cited for the convention of writing a range as “including the lower bound and excluding the upper bound,” 0 <= i < n. He lists as advantages that in this form the number of elements equals the difference of the bounds, n - 0, that an empty range can be written naturally as 0 <= i < 0, and that, if indices start from 0, “the index of an element equals the number of elements before it.” Source: E. W. Dijkstra, Why numbering should start at zero (EWD831, 1982).

MISO also holds ranges in the form [begin, end). For example, the range of a parallel loop, Range1D{0, n}, represents 0 <= i < n (Chapter 7).

4.3 Nested loops

int visits = 0;
for (int k = 0; k < 2; ++k) {
  for (int j = 0; j < 3; ++j) {
    for (int i = 0; i < 4; ++i) {  // the innermost i changes fastest
      ++visits;
    }
  }
}
// visits is 2 x 3 x 4 = 24

This is the same as nested do loops in Fortran. Simply writing an ordinary for does not run it in parallel. You will learn how to write parallel computations in Chapter 7.

A real MISO example — InitialCondition::apply in demo/mhd2d_vortex/src/main.cpp is a triple loop that sets initial values at all grid points. We excerpt part of it and add explanatory comments where code is omitted.

for (int k = 0; k < grid.k_total; ++k) {
  for (int j = 0; j < grid.j_total; ++j) {
    for (int i = 0; i < grid.i_total; ++i) {
      qq.ro(i, j, k) = 1.0;
      qq.vx(i, j, k) = -v0 * util::sin(2.0 * pi<Real> * grid.y[j]);
      qq.vy(i, j, k) = +v0 * util::sin(2.0 * pi<Real> * grid.x[i]);
      // (same for the other variables)
    }
  }
}

grid.i_total is the number of grid points in the x direction (including ghost cells), and ghost cells are margin cells added at the edges of the domain (Chapter 8). i runs from 0 to grid.i_total - 1. qq.ro(i, j, k) is an element of the density array, and grid.x[i] is the x coordinate. The mechanism that lets you write indices in parentheses, as in qq.ro(i, j, k), is explained in Chapter 4, and pi<Real> (the constant for pi expressed as type Real) in Chapter 6. For now, it is enough to read this as “loop over the indices in the three directions starting from 0 and put a value at each point.”

4.4 while

Use while to repeat based on a condition rather than a count. It corresponds to Fortran’s do while.

double residual = 1.0;
int iterations = 0;
while (residual > 0.1) {  // check the condition before every iteration
  residual *= 0.5;
  ++iterations;
}
// residual is 0.0625, iterations is 4

If you do not update the variables that affect the condition inside the block, the loop never ends. There is also do { ... } while (condition);, which executes once first and then checks the condition.

A real MISO example — ModelBase::run in miso/include/miso/mhd_model_base.hpp, the loop that advances the time evolution.

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

Until the current time time.time reaches the end time time.tend, the loop repeats one-step updates d.update() and, when needed, output. d is a name that refers to the model for each problem; the mechanism is explained in Chapter 5.

4.5 break, continue, and return

C++	Action	Fortran
`break;`	Finish the current loop and go to the statement right after it	`exit`
`continue;`	Skip the rest of this iteration and go to the next one	`cycle`
`return;`, `return value;`	Finish the function itself and return to the caller	`return`

int total = 0;
for (int i = 0; i < 6; ++i) {
  if (i == 2) {
    continue;  // skip the addition and go to ++i
  }
  if (i == 5) {
    break;     // exit the loop
  }
  total += i;  // adds 0, 1, 3, 4
}
// total is 8

C++’s continue is different from Fortran’s continue statement (a statement that does nothing). It corresponds to Fortran’s cycle. If you write break inside the inner loop of nested loops, only the inner loop is exited.

A real MISO example — the whole of parse_config_filepath in miso/include/miso/config.hpp (for inline, see Section 6; for std::optional, Section 2.4; for throw, Section 9).

inline std::optional<std::string> parse_config_filepath(int argc, char **argv) {
  for (int i = 1; i < argc; ++i) {
    std::string_view arg(argv[i]);

    if (arg == "--config") {
      if (i + 1 >= argc)
        throw std::runtime_error("--config requires filepath");
      return std::string(argv[i + 1]);
    }

    constexpr std::string_view key = "--config=";
    if (arg.size() >= key.size() && arg.compare(0, key.size(), key) == 0) {
      return std::string(arg.substr(key.size()));
    }
  }
  return std::nullopt;
}

argv[i] is the string of the i-th command-line argument (you will read the type of argv in Chapter 3). std::string_view (<string_view>) is a type for reading a string without copying it; you can compare it with == and use size() (the length), compare (comparing a part), and substr (extracting a part). argv[0] is the name of the program itself, so the loop starts at i = 1 and runs up to argc - 1. For the form --config path, it returns the next argument argv[i + 1]; for the form --config=path, it returns what follows =; and in either case it finishes the function at that point. If the loop ends without finding anything, the final return std::nullopt; returns “no value.” Leaving the whole function with return instead of break out of the loop is commonly used in searches.

4.6 switch

You can also use switch to branch on an integer or enumeration value. It is close to Fortran’s select case.

enum class Axis { X, Y, Z };
Axis axis = Axis::Y;
int n = 0;
switch (axis) {
case Axis::X:
  n = 64;
  break;
case Axis::Y:  // execution starts here
  n = 32;
  break;       // leave the switch
case Axis::Z:
  n = 1;
  break;
}
// n is 32

Execution does not stop at the boundary between case labels. If you forget break, the statements of the next case also run. Write the handling for values that match no case under default:.

A real MISO example — symmetric_index in miso/include/miso/boundary_condition.hpp computes, on the inner and outer sides of the boundary, the index of a ghost cell and the index of the source cell whose value is copied to it.

switch (side) {
case Side::INNER:
  i_ghst = i;
  i_trgt = 2 * i_margin - i - 1;
  break;
case Side::OUTER:
  i_ghst = i_total - i_margin + i;
  i_trgt = i_total - i_margin - i - 1;
  break;
}

side is the Side type from Section 2.7. If the ghost-cell width i_margin is 2, then on the inner boundary (i_ghst, i_trgt) is (0, 3) and (1, 2) for i = 0, 1, giving pairs of positions mirrored across the boundary. The mechanism that returns the resulting i_ghst and i_trgt to the caller (reference parameters) is covered in Chapter 3.

Exercise: loop conditions

For the following loop, give the value of i each time the block is entered and the final total. Also find total when the condition is changed to i < 5.

int total = 0;
for (int i = 1; i <= 5; i += 2) {
  if (i == 3) {
    continue;
  }
  total += i;
}

::

i takes the values 1, 3, and 5. At 3 the addition is skipped, so total is 1 + 5 = 6. Even with continue, the update i += 2 is executed. If you change the condition to i < 5, i takes only 1 and 3, and total is 1.

5. Functions

5.1 Declaration, definition, and call

In C++, both computations that return a value and computations corresponding to Fortran subroutines are written as functions.

double square(double x);   // declaration

double square(double x) {  // definition
  return x * x;
}

Part	Meaning
The leading `double`	The type of the returned value (the return type)
`square`	The function name
`(double x)`	The types of the parameters and the names used inside the function
The `;` at the end of a declaration	Indicates that the body is not written here
The `{ ... }` of a definition	The statements to execute
`return`	Return a value to the caller and end the function

The caller writes double area = square(2.0);, and area becomes 4.0. A function must be declared before the place where it is called. If the definition comes later, write only the declaration beforehand.

A function that returns no value has the return type void. It corresponds to Fortran’s subroutine. Call it as write_output(); without call. The parentheses cannot be omitted even when there are no arguments.

You can also give a parameter a default value (a default argument).

double scale(double x, double factor = 1.0) {
  return x * factor;
}
// scale(2.0) is 2.0, scale(2.0, 3.0) is 6.0

A real MISO example — rank in miso/include/miso/env.hpp is a function that takes no arguments and returns an int.

inline int rank() {
  if (!is_initialized()) {
    throw std::runtime_error(
        "mpi::Env must be instantiated at the beginning of the program.");
  }
  return world::rank;
}

If MPI has not been set up, it throws an error; otherwise it returns the stored rank number world::rank (world:: is covered in Section 7).

5.2 Pass by value

If you write a parameter as double x, a copy of the caller’s value is passed to the function. Even if you modify x inside the function, the caller’s variable does not change.

void set_copy(double value) {
  value = 5.0;  // only the copy inside the function changes
}

double a = 1.0;
set_copy(a);    // a is still 1.0 after the call

In Fortran, arguments are normally treated as the caller’s variable itself, and modifying one inside a subroutine also changes the caller’s variable. A C++ parameter such as double value does not behave this way. Pass by reference, which you use when you want to modify the caller’s variable, is explained in Chapter 3. Passing a copy of a value in this way is called pass by value. For a small value such as a single number, receiving it by value is the basic approach.

A real MISO example — zfill in miso/include/miso/utility.hpp is a function that turns an integer into a zero-padded string.

/// @brief zero-fill integer to string
/// @param num target integer number
/// @param width width of the string
/// @return zero-filled string representation of the integer
inline std::string zfill(int num, int width) {
  std::ostringstream oss;
  oss << std::setw(width) << std::setfill('0') << num;
  return oss.str();
}

It receives two int values by value and returns a std::string. zfill(3, 8) is "00000003". The body writes num into oss, a stream that writes into a string, with width width and zero padding, and returns the resulting string (Section 8). This function is used to make the numbers in the MHD output file name mhd.00000003.00000000.bin (output number and rank number) and in the time file names (miso/include/miso/mhd_checkpoint.hpp, time.hpp). A comment beginning with /// above a function is the notation used by a tool called Doxygen to generate API documentation.

5.3 Overloading

In C++, you can define multiple functions with the same name if their parameter types or counts differ. This is called overloading.

int square(int x) { return x * x; }
double square(double x) { return x * x; }

int n = square(3);       // calls the int version; n is 9
double y = square(1.5);  // calls the double version; y is 2.25

Which one is called is determined at compile time from the types of the arguments you pass. Functions that differ only in return type cannot be distinguished. In Fortran, you declare a generic name with interface to group individual functions, but in C++ you simply define them side by side with the same name.

A real MISO example — the namespace miso::util in miso/include/miso/utility.hpp has math functions for float and for double side by side. The source has only the single line MISO_MATH_UTILS(sqrt, sqrtf), but this is a macro (Section 10), and when you build for the CPU it is replaced by the following two function definitions.

// result of expanding MISO_MATH_UTILS(sqrt, sqrtf) (line breaks adjusted for clarity)
inline float sqrt(float x) { return std::sqrt(x); }
inline double sqrt(double x) { return std::sqrt(x); }

The float version and the double version are defined under the same name sqrt. If you write util::sqrt(x), the version for the type of the argument is called, whether Real is float or double. sin, cos, exp, log, fabs, and others are defined in the same way.

You can check this with the following program (specify the location of the MISO headers with -I).

#include <iostream>

#include <miso/constants.hpp>
#include <miso/utility.hpp>

int main() {
  using Real = float;
  std::cout << miso::util::zfill(3, 8) << "\n";
  Real r = miso::util::sqrt(Real(2));    // float version
  double d = miso::util::sqrt(2.0);      // double version
  std::cout.precision(10);
  std::cout << r << " " << d << "\n";
  std::cout << miso::pi<Real> << " " << miso::pi<double> << "\n";
}

g++ -std=c++17 -I/path/to/MISO/miso/include miso_util.cpp -o miso_util
./miso_util

00000003
1.414213538 1.414213562
3.141592741 3.141592654

The result of the float version differs from the double version at about the 7th digit. miso::pi<Real> is the value of pi in miso/include/miso/constants.hpp, and the notation that specifies a type in <> is explained in Chapter 6.

Exercise: type conversion and overloading

With using Real = float;, answer in which type the right-hand side of each of the following lines is computed and which version of square is called. square refers to the two versions in Section 5.3.

Real a = 2.0f;
auto b = a * 0.5;
auto c = a * Real(0.5);
auto d = square(3) / 2;
auto e = square(3.0) / 2;

::

b: 0.5 is a double, so a is converted to double and the calculation is done in double. b is the double value 1.0.
c: This is a calculation between float values, and c is the float value 1.0.
d: The int version of square is called and gives 9. 9 / 2 is integer division, so d is the int value 4.
e: The double version is called and gives 9.0. 9.0 / 2 is the double value 4.5.

If you receive the result with auto, as with b, it becomes double even though you meant float. This is why MISO writes constants as Real(0.5).

6. Splitting files and headers

When a program gets long, you split the functions into separate files. In C++, you write declarations and definitions that other files will use in a header (.hpp) and #include it on the side that uses them.

physics.hpp   header with the function
main.cpp      includes "physics.hpp" and uses the function

#include "physics.hpp" pastes the contents of the file as they are at that position before compilation. Fortran’s use makes the names in an already-compiled module available, so the mechanism is different. By convention, headers of your own project are enclosed in "...", and headers of the standard library or installed libraries are enclosed in <...>. The compiler also searches for headers in the directories you specify with -I (Chapter 1).

#pragma once: If the same header is #included many times within one .cpp (this happens when a header includes other headers), the same definitions are duplicated and cause an error. If you write #pragma once at the top of a header, the second and later inclusions are ignored. Every MISO header has this line.

inline: If you write a function definition (the body) in a header and #include that header from two .cpp files, you get two definitions of the same function, which causes a link error. In the next example, both main.cpp and a.cpp #include physics.hpp.

// physics.hpp
#pragma once

double kinetic_energy(double ro, double v) {
  return 0.5 * ro * v * v;
}

g++ -std=c++17 main.cpp a.cpp -o odr

/usr/bin/ld: ...: in function `kinetic_energy(double, double)':
a.cpp:(.text+0x0): multiple definition of `kinetic_energy(double, double)'; ...main.cpp:(.text+0x0): first defined here

#pragma once only prevents duplication within a single .cpp; it cannot prevent duplicate definitions that arise in separate .cpp files. Adding inline to a function means “the same definition may appear in multiple .cpp files,” which makes it possible to link.

inline double kinetic_energy(double ro, double v) {
  return 0.5 * ro * v * v;
}

The name inline comes from the optimization of “expanding the body at the call site,” but its role in current C++ is that it lets you “place a definition in a header.” Whether the function is actually expanded is decided by the compiler.

MISO is a header-only library. All function and class definitions are in the headers (miso/include/miso/*.hpp), and there is no .cpp to compile separately as a library. For this reason, MISO’s convention is to add inline to free functions (functions outside a class) that are not templates. This is why inline is attached to zfill in Section 5.2 and rank in Section 5.1. (Templates and functions defined inside a class may appear in multiple .cpp files even without inline. They are covered in Chapter 4 and Chapter 6.)

A real MISO example — miso/include/miso/core.hpp is a header that does nothing but #include other headers together. Here is an excerpt of the beginning.

#pragma once

// Umbrella header for the common layer shared by MHD and RT.
// Include this together with a physics header, e.g.
//   #include <miso/core.hpp>
//   #include <miso/mhd.hpp>  // or <miso/rt.hpp>

// fundamentals
#include "backend.hpp"
#include "constants.hpp"
#include "types.hpp"
#include "utility.hpp"

The demo main.cpp loads all the MISO headers it needs with just two lines, #include <miso/core.hpp> and #include <miso/mhd.hpp>. Headers inside MISO include each other with "..." because they are in the same directory, and the demo includes them with <miso/...> because it searches the location given by -I.

7. Namespaces

Names are divided into groups so that functions with the same name in different libraries can be told apart. Such a group is called a namespace.

namespace eos {
double gamma = 5.0 / 3.0;
}

namespace grid {
int gamma = 3;  // same name, but a different variable from eos::gamma
}

int main() {
  std::cout << eos::gamma << " " << grid::gamma << "\n";  // 1.66667 3
}

To use a name from outside, write namespace::name. :: means “inside this namespace.” All names in the standard library are inside the namespace std, which is why we write std::cout and std::string. Namespaces can be nested, and you follow them from the outside in, as in a::b::name. They play a role similar to distinguishing names in Fortran modules, but they are unrelated to files, and one namespace can be written across multiple files.

There are two ways to write names more briefly.

After using namespace eos;, you can omit eos:: and write gamma. This is similar to Fortran’s use eos, but it makes name collisions more likely, so by convention it is not used in headers.
Writing namespace fs = std::filesystem; gives a namespace a short alias.

A real MISO example — all MISO names are in the namespace miso, which is further divided by area.

Namespace	Examples of contents	File
`miso`	`Config`, `Env`, `Array3D`, `for_each`, `Direction`	`config.hpp`, `env.hpp`, `array3d.hpp`, etc.
`miso::util`	`zfill`, `sqrt`, `pow2`, `clamp`	`utility.hpp`
`miso::mpi`	`rank()`, `size()`, `is_root()`	`env.hpp`, `mpi_util.hpp`
`miso::mhd`	`ModelBase`, `Fields`, `Integrator`	`mhd_*.hpp`
`miso::eos`	`IdealEOS`	`eos.hpp`
`miso::backend`	`Host`, `CUDA`	`backend.hpp`

At the beginning of miso/include/miso/env.hpp, namespaces are nested three levels deep and hold the variable that stores MPI information.

namespace miso {

/// @brief Global environment for MPI
namespace mpi {

// global (multi-node) MPI info
namespace world {
inline bool is_initialized = false;
inline MPI_Comm comm = MPI_COMM_NULL;
inline int rank = -1;
inline int size = -1;
}  // namespace world

From outside, this is the variable named miso::mpi::world::rank, the value returned by the function rank() in Section 5.1. inline on a variable has the same meaning as on a function: a variable placed in a header is treated as a single variable across multiple .cpp files (since C++17). The // namespace world after } is a comment showing which namespace ends there. The function miso::mpi::rank() and the variable miso::mpi::world::rank are distinguished because their namespaces differ.

miso/include/miso/utility.hpp gives an alias with namespace fs = std::filesystem; and uses it as in fs::create_directories(...). The demo main.cpp writes using namespace miso; at the beginning, so you can write miso::Config as Config and miso::mhd::ModelBase as mhd::ModelBase.

8. Standard library: input/output, strings, and vector

The collection of components that come with C++ is called the standard library. Here we introduce three that are often used in numerical code.

#include <iostream>
#include <string>
#include <vector>

int main() {
  std::string name = "mhd";
  std::string file = name + ".00000003.bin";  // concatenate with +
  std::cout << file << " (" << file.size() << " chars)\n";

  std::vector<double> x(4);  // 4 doubles, all 0.0
  for (int i = 0; i < static_cast<int>(x.size()); ++i) {
    x[i] = 0.5 * i;           // indices run from 0 to 3
  }
  x.push_back(2.0);           // append to the end, making 5
  std::cout << x.size() << " " << x[4] << "\n";
}

mhd.00000003.bin (16 chars)
5 2

std::cout (<iostream>): sends values in order with << to print them. A newline is written as "\n" or std::endl. The number of digits and the width of the output are specified with std::setw, std::setprecision, and so on from <iomanip>. Writing to a file uses std::ofstream (<fstream>), which is used in the same way, and writing to a string uses std::ostringstream (<sstream>).
std::string (<string>): a string whose length can change. You concatenate with + and get the length with size(). Unlike Fortran’s fixed-length character, no trailing spaces are added. To convert a number to a string, use std::to_string(value).
std::vector<double> (<vector>): a one-dimensional array whose number of elements can be decided at run time. It is close to a one-dimensional allocatable array in Fortran: you allocate it with x(4), and it is released automatically when the variable’s lifetime ends. Indices are written in square brackets as x[i] and run from 0 to x.size() - 1. The part <double> specifies the element type (this notation is covered in Chapter 6).

A function that you call by adding . after a variable name, such as x.size(), is called a member function (Chapter 4). x.size() returns a non-negative integer type, so when comparing it with the int variable i, the code aligns the types with static_cast<int>.

MISO’s array computations use MISO’s own Array1D to Array4D, not std::vector, because arrays that can be placed in GPU memory are needed (Chapter 7). std::vector is used for one-dimensional data that is used only on the CPU side, such as grid coordinates.

A real MISO example — Grid in miso/include/miso/grid.hpp holds coordinates in a std::vector (excerpt).

/// @brief coordinate in x direction
std::vector<Real> x;
/// @brief coordinate in y direction
std::vector<Real> y;

Time::log in miso/include/miso/time.hpp prints the time and step count to the screen (excerpt of the body).

if (mpi::is_root()) {
  std::cout << std::fixed << std::setprecision(2) << "time = " << std::setw(6)
            << time << ";  n_step = " << std::setw(8) << n_step
            << ";  n_output = " << std::setw(8) << n_output << std::endl;
}

Only MPI process 0 prints. std::fixed << std::setprecision(2) specifies “print in fixed-point with 2 digits after the decimal point,” and std::setw(6) specifies “right-align the next value in a width of 6,” which corresponds to the Fortran format (f6.2). In Time::time_filepath in the same file, the output file name is built by string concatenation, as in "time." + util::zfill(...) + ".txt".

9. Error handling with exceptions

When a state arises in which processing cannot continue, such as a file that cannot be opened or an invalid setting value, you report it by throwing an exception.

#include <stdexcept>
#include <string>

double checked_dt(double dt) {
  if (!(dt > 0.0)) {
    throw std::runtime_error("dt must be positive: " + std::to_string(dt));
  }
  return dt;
}

When throw is executed, the function ends there, and control goes back through the callers in order. std::runtime_error (<stdexcept>) is a standard exception that holds a string describing the error. If the caller uses try and catch, it can catch the exception and continue processing.

try {
  checked_dt(0.1);
  checked_dt(-1.0);  // an exception is thrown here
  std::cout << "not printed\n";
} catch (const std::runtime_error &e) {  // const and & are explained in Chapter 3
  std::cout << "error: " << e.what() << "\n";
}

error: dt must be positive: -1.000000

An exception that no catch handles terminates the program abnormally. In many environments, the string from what() is displayed at that time. This is used much like Fortran’s stop "message" or error stop, but the difference is that it leaves a way for the caller to catch the exception and recover.

A similar facility is assert(condition) (<cassert>). It stops the program if the condition is false, but in a Release build it is removed and checks nothing. In MISO, use assert to check for program mistakes (such as an index out of range), and use exceptions for errors that can occur at run time, such as a file that cannot be opened.

A real MISO example — the constructor of Config in miso/include/miso/config.hpp throws an exception if the configuration file does not exist (excerpt).

if (!fs::exists(load_filepath)) {
  throw std::runtime_error("Config file not found: " + load_filepath);
}

MISO demos do not catch exceptions, so if you specify a nonexistent configuration file, the program stops as follows (the beginning of the output when we built and tried mhd2d_vortex for the CPU).

./mhd2d_vortex --config=no_such.yaml

terminate called after throwing an instance of 'std::runtime_error'
  what():  Config file not found: no_such.yaml

The what(): line shows the string passed to throw. Reading this line tells you the cause. The Time constructor in Section 4.1 and rank() in Section 5.1 have the same form.

10. The preprocessor

A line that begins with # is an instruction processed by the preprocessor before compilation. #include is one of them. The preprocessor does not understand C++ syntax; it only replaces text and keeps or removes lines.

Directive	Meaning
`#include <file>`	Paste the contents of the file at that position
`#define NAME value`	Replace later occurrences of `NAME` with `value` (a macro)
`#ifdef NAME` 〜 `#else` 〜 `#endif`	Keep the first half if `NAME` is defined, and the second half if not
`#ifndef NAME`	Keep it if `NAME` is not defined

In Fortran as well, you sometimes give a file the extension .F90 and use #ifdef. This is the same thing.

#include <iostream>

#define N_GHOST 2

#ifdef USE_DOUBLE
using Real = double;
#else
using Real = float;
#endif

int main() {
  std::cout << "N_GHOST = " << N_GHOST << ", sizeof(Real) = " << sizeof(Real) << "\n";
}

g++ -std=c++17 pp.cpp -o pp && ./pp                 # N_GHOST = 2, sizeof(Real) = 4
g++ -std=c++17 -DUSE_DOUBLE pp.cpp -o pp && ./pp    # N_GHOST = 2, sizeof(Real) = 8

-DUSE_DOUBLE is a compiler option that says “start with USE_DOUBLE defined.” sizeof(type) is the size of that type in bytes. For constants, using constexpr from Section 2.3 is safer than a macro because it has a type. Use macros for switching between builds and for repetition that cannot be written with types.

A real MISO example — demo/mhd2d_vortex/src/main.cpp switches the hardware used for the computation with USE_CUDA.

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

When you run CMake with -DUSE_CUDA=ON, target_compile_definitions(${app_name} PRIVATE USE_CUDA) in the demo’s CMakeLists.txt defines USE_CUDA, and Backend becomes backend::CUDA, which represents the GPU. Otherwise it is backend::Host, which represents the CPU. How Backend is used is explained in Chapter 6 and Chapter 7.

miso/include/miso/cuda_compat.hpp switches on __CUDACC__, which the CUDA compiler (nvcc) defines.

#ifdef __CUDACC__

#include <cuda_runtime.h>

#else  // __CUDACC__

#ifndef __device__
#define __device__
#endif

#ifndef __host__
#define __host__
#endif

When compiling for the GPU, it loads the CUDA headers; when compiling for the CPU (g++ only), it replaces __host__ and __device__ with nothing. This way, a function marked with the GPU markers __host__ __device__ can be compiled by g++ as it is. The meaning of these markers is explained in Chapter 9.

You can also write macros that take arguments. MISO_MATH_UTILS from Section 5.3 is defined as follows in miso/include/miso/utility.hpp (excerpt of the CPU part).

#define MISO_MATH_UTILS(name, namef)                                             \
  __host__ __device__ inline float name(float x) { return std::name(x); }        \
  __host__ __device__ inline double name(double x) { return std::name(x); }

The \ at the end of a line means “continues on the next line,” so the macro definition is treated as a single line. Writing MISO_MATH_UTILS(sqrt, sqrtf) replaces name with sqrt and yields two function definitions, one for float and one for double. (The sqrtf in namef is used as the function name for float in the GPU definition.) The macro saves the effort of writing many pairs of functions of the same form.

11. Exercise: reading the beginning of the demo `main.cpp`

The beginning of demo/mhd2d_vortex/src/main.cpp is as follows. Explain what each line (or each group of lines) does, using the terms from this chapter.

#include <miso/core.hpp>
#include <miso/mhd.hpp>

using namespace miso;

using Real = float;

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

::

#include <miso/core.hpp>, #include <miso/mhd.hpp>: These paste in the common part of MISO and the MHD header (Section 6). MISO is header-only, so these alone make all of MISO’s definitions available. Because of the <...>, the compiler searches the location given by -I (in CMake, this is set when you link MISO::MISO).
using namespace miso;: From here on, you can write names in the namespace miso without miso:: (Section 7). backend::CUDA is really miso::backend::CUDA.
using Real = float;: This makes the real-number type Real single precision (Section 2.6). If you change this line to double, the computation runs in double precision.
#ifdef USE_CUDA to #endif: If the macro USE_CUDA is defined (-DUSE_CUDA=ON in CMake), Backend becomes an alias of backend::CUDA; otherwise it becomes an alias of backend::Host (Section 10). Only one of the two lines is passed to the compiler.

The two types Real and Backend are decided at the beginning, and the rest of the code (Model and so on) is written using only these two names. This is a structure in which precision and hardware can each be switched by changing one line.

For details of the syntax, see the C++17 draft standard N4659: if and switch, for, while and do, break, continue, and return, && and ||, and the conditional operator. For unfamiliar terms, see the Glossary; for how Fortran constructs correspond to C++, see the Fortran Correspondence Table.