4. Classes and Structs

Building on the variables, functions, references, and pointers from Chapter 3, you will define your own types that bundle data together with the functions that operate on it. By the end of this chapter, you will be able to read what data the MISO array Array3D, the MHD variable bundle Fields, the time manager Time, the configuration reader Config, and the demo’s InitialCondition hold, how they are created, and how they are cleaned up.

4.1 Defining a Type and Creating Objects
4.2 Access Control: public and private
4.3 Constructors and Member Initializer Lists
explicit and implicit conversion
4.4 const Member Functions and this
4.5 Destructors and Cleanup (RAII)
4.6 Prohibiting Copies, and Moves
4.7 Operator Overloading and operator()
4.8 Owning Arrays and Non-Owning Views
4.9 Reading the demo’s InitialCondition
Exercises

The MISO excerpts in this chapter contain names with angle brackets <>, such as template <typename T> and Array3D<Real, Backend>. These are “marks that specify a type later.” For example, read Array3D<float> as “an Array3D whose elements are float,” and Time<Real> as “a Time that holds time as type Real.” The mechanism is explained in Chapter 6. In this chapter, read what is inside the angle brackets as “the type specification used by that class” and move on.

4.1 Defining a Type and Creating Objects

We start with an example that bundles the grid counts nx, ny, and nz into one. The following GridSize is a type for explanation.

struct GridSize {
  int nx = 1;
  int ny = 1;
  int nz = 1;

  int cell_count() const {
    return nx * ny * nz;  // multiplies the cell counts in the 3 directions and returns the total number of cells
  }
};

struct GridSize { ... }; defines a type named GridSize. This definition alone does not yet create a variable for any particular grid, just as writing the type name int does not create an integer variable. Variables written inside a type are called member variables, and functions are called member functions. The = 1 in int nx = 1; is the initial value of that member. The ; after the closing brace cannot be omitted.

Written as a Fortran derived type, the data part corresponds to the following.

type :: grid_size
  integer :: nx = 1, ny = 1, nz = 1
end type grid_size

In C++, besides data, you can also write the functions that use that data (here cell_count) inside the same type.

Here is an example that creates variables of the type and uses the members.

GridSize a{4, 3, 2};  // initializes a's nx, ny, nz to 4, 3, 2
GridSize b{8, 3, 2};  // b has its own nx, ny, nz, separate from a
a.nx = 5;             // only a.nx becomes 5. b.nx stays 8
std::cout << a.cell_count() << " " << b.cell_count() << "\n";  // 30 48

GridSize c;           // uses the initial value 1 written in the declaration. c.cell_count() is 1

a and b each have their own storage for nx, ny, and nz. Something that has a region for storing values according to its type is called an object, and an object created from a class type is also called an instance.

. is the symbol that specifies a member of an object, and corresponds to the % in Fortran’s a%nx. While a.cell_count() is being called, nx inside the function means a.nx. For b.cell_count(), it is b.nx. Which member variables a member function uses depends on which object it was called from. The const at the end of cell_count() specifies that “this function does not modify member variables.” It is explained in Section 4.4.

To use a member through a pointer, write -> instead of ..

GridSize *p = &a;  // p points to a
p->ny = 6;         // same as (*p).ny = 6. a.ny becomes 6

p->ny means “the ny of the object that p points to.” Use . with an object or a reference to it, and -> with a pointer.

In C++, a type that bundles data and functions is called a class. A class can be defined with either the struct or the class keyword. Something defined with struct is also called a structure, but in C++ it is a kind of class. The difference between the two is explained in the next section.

A MISO example: Time, which manages time

Time in miso/include/miso/time.hpp is a struct that bundles the time and the step count. We extract the member variable part (the original comments /// @brief ... have been replaced with short English comments).

template <typename Real> struct Time {
  Real time;                  // current time
  Real tend;                  // end time
  Real dt_output;             // time interval between outputs
  int n_step;                 // time step number
  int n_output;               // output number
  int n_output_digits;        // number of digits in output file names
  std::string time_save_dir;  // directory where the time file is written
  bool io_enabled;            // whether to do input/output
  // ... (member functions are covered from Section 4.3)
};

The template <typename Real> at the top is a mark that specifies the type Real for representing time later (Chapter 6). The demo uses using Real = float;, so time and tend are float.

The Time object is created as the member time of the class that manages MISO’s time evolution (ModelBase, which you will read in Chapter 5). ModelBase::run() in miso/include/miso/mhd_model_base.hpp contains the following line.

while (time.time < time.tend) {

The time on the left is the name of the object, and .time is its member variable. The spelling is the same, but the roles differ. Read it as “while the current time time of the object time is smaller than the end time tend.”

4.2 Access Control: public and private

In GridSize in the previous section, the code using it could freely read and write any member. Here we create a class that represents a one-dimensional uniform grid, and decide which members can be read and written from outside. The number of cells and the grid spacing are decided at creation, and afterward only reading them is allowed.

class Axis {
public:
  // divides the interval [x_min, x_max] into n equal cells
  explicit Axis(int n, double x_min = 0.0, double x_max = 1.0)
      : n_(n), x_min_(x_min), dx_((x_max - x_min) / n) {}

  int size() const { return n_; }      // only reads the number of cells
  double dx() const { return dx_; }    // only reads the grid spacing
  double x(int i) const {              // coordinate of the center of cell i (i starts from 0)
    return x_min_ + (i + 0.5) * dx_;
  }

private:
  int n_;
  double x_min_;
  double dx_;
};

Members written after public: can be used from the code that uses Axis, such as main(). n_, x_min_, and dx_, written after private:, cannot be read or written directly by the user and can be used only from inside the member functions of Axis. The part at the top beginning with Axis(...) is the processing when an object is created, and is explained in the next section. explicit is explained in the explicit part of Section 4.3.

On the user side, write the following.

Axis a(4);              // divides the interval [0, 1] into 4
Axis b(10, -1.0, 1.0);  // divides the interval [-1, 1] into 10
std::cout << a.size() << " " << a.dx() << " " << a.x(0) << "\n";  // 4 0.25 0.125
std::cout << b.size() << " " << b.dx() << " " << b.x(9) << "\n";  // 10 0.2 0.9
// a.n_ = 8;           // error: a private member cannot be assigned from outside

If you enable the last line, which is commented out, g++ gives the error 'int Axis::n_' is private within this context. Because no operation that changes the number of cells is made public, it cannot happen that only the cell count changes after creation and no longer matches the grid spacing. Bundling data and the functions that operate on it into a class and deciding which operations to expose to the outside in this way is called encapsulation.

The only difference between struct and class is the default when no access is written. A struct starts as public, and a class starts as private. The trailing _ in member names is a naming convention for distinguishing private member variables, not C++ grammar.

A MISO example: Array3D with private members, and structs with everything public

MISO uses the two styles according to their roles.

Array3D in miso/include/miso/array3d.hpp is a class, and keeps the array’s start address data_ and its size shape_ private. Users go through public functions such as extent(), size(), and a(i, j, k). This is so that users cannot arbitrarily rewrite the size of the allocated array.
Time, Config, Fields, and others are structs with all member variables public. They are designed as bundles of data that MISO’s internals and the demo use directly, as in time.n_step and qq.ro.

The Array3D excerpt is read from Section 4.5 onward.

4.3 Constructors and Member Initializer Lists

A special member function with the same name as the class, like Axis(...), is called a constructor. It is called automatically when an object is created and initializes the member variables. You do not write a return type (not even void).

explicit Axis(int n, double x_min = 0.0, double x_max = 1.0)
    : n_(n), x_min_(x_min), dx_((x_max - x_min) / n) {}

Read these two lines in the following parts.

Part	Meaning
`Axis(int n, double x_min = 0.0, double x_max = 1.0)`	A constructor that takes three arguments
`= 0.0`、`= 1.0`	A default argument. If the caller omits it, this value is used
From `:` to just before `{`	The member initializer list
`n_(n)`	Initializes the member `n_` with the value of the argument `n`
`dx_((x_max - x_min) / n)`	Initializes the member `dx_` with the value of the expression
`{}`	The body executed after initialization. It is empty in this example

Writing Axis a(4); uses the default values 0.0 and 1.0 for x_min and x_max. Default arguments can be omitted starting from the last one, and you cannot omit only a middle argument. Writing with braces, as in Axis a{4};, calls the same constructor. Note that writing empty parentheses, Axis a();, is interpreted as “the declaration of a function a that returns an Axis.”

Fortran has no procedure that is always executed when a type is created. In C++, an object of a class with a constructor is always created through a constructor. Axis has only a constructor that takes arguments, so you cannot create it with no arguments as Axis a; (it is an error). For GridSize, which has no constructor written, the compiler provides a constructor with no arguments, so GridSize c; worked.

The order in which members are initialized

Members are initialized not in the order written in the initializer list, but in the order the member variables are declared in the class. In Axis, n_, x_min_, and dx_ are declared in this order, and the initializer list is written in the same order. In the following BadAxis, dx_ is declared before n_, and n_ is used to initialize dx_.

class BadAxis {
public:
  BadAxis(int n, double x_min, double x_max)
      : n_(n), x_min_(x_min), dx_((x_max - x_min) / n_) {}

private:
  double dx_;  // declared before n_, so it is initialized first
  int n_;
  double x_min_;
};

At the time dx_ is initialized, n_ does not yet hold a value. With -Wall, g++ warns about the order mismatch (-Wreorder) and about using n_, which has no value (-Wuninitialized). Write the initializer list in the same order as the declarations, and check that order when an initialization uses another member. If you use the argument n instead of a member, as in Axis above, the result does not depend on the declaration order.

Members whose values can be set only in the initializer list

Members declared const and reference members cannot be assigned after they are created. Their values or referents are set in the initializer list. This is because, by the time execution enters the constructor body { ... }, the members’ initialization is already finished.

struct Scaled {
  const double factor;   // not changed after construction
  double &target;        // refers to an outside variable (Chapter 3)

  Scaled(double f, double &t) : factor(f), target(t) {}
  // Scaled(double f, double &t) { factor = f; target = t; }  // error
};

A MISO example: the constructor of Time

We read the constructor of Time in miso/include/miso/time.hpp (the latter half, which creates directories, is omitted).

  /// @brief Default constructor
  Time(const Config &config)
      : tend(config["time"]["tend"].as<Real>()),
        dt_output(config["time"]["dt_output"].as<Real>()),
        n_output_digits(config["io"]["n_output_digits"].as<int>()),
        io_enabled(config.yaml_obj["io"]["enabled"].as<bool>()) {
    if (!(tend > 0)) {
      throw std::runtime_error("time.tend must be positive.");
    }
    if (!(dt_output > 0)) {
      throw std::runtime_error("time.dt_output must be positive.");
    }

    initialize();
    // ... (decides time_save_dir and creates the directory)
  }

The argument const Config &config is a reference for reading the Config, which holds the contents of the configuration file, without copying it (Chapter 3).
config["time"]["tend"] is how to extract the value of tend: under time: in config.yaml. The mechanism of [] is explained in Section 4.7. .as<Real>() is a yaml-cpp function that reads the extracted value as a number of type Real.
The initializer list initializes four members in declaration order (tend, dt_output, n_output_digits, io_enabled).
The body checks the values, then calls initialize() to set time, n_step, and n_output to 0. If a value is invalid, it throws an exception and the object is not created.

The original comment says Default constructor, but in C++ terminology a “default constructor” is a constructor that can be called with no arguments, and this constructor takes a Config. The word in the comment does not match the C++ term, so we rely on the explanation in the text.

explicit and implicit conversion

We added explicit to the constructor of Axis. explicit specifies that “this constructor is not used for automatic conversion from another type.”

Because there is a default argument, an Axis can be created from a single integer. Without explicit, C++ also performs “conversion from an integer to an Axis” automatically. The following two lines are errors with explicit and compile without it.

double length(const Axis &axis) { return axis.size() * axis.dx(); }

// Axis c = 4;   // without explicit, Axis(4) is created automatically
// length(4);    // without explicit, a temporary Axis(4) is created and passed
length(Axis(8));  // explicitly creating it works even with explicit (returns 1)

A conversion in which the compiler changes the type for you, without you writing the type change, is called an implicit conversion. The role of explicit here is to prevent an unintended grid from being created when you pass an integer by mistake in place of the intended argument.

IdealEOS, the ideal-gas equation of state in miso/include/miso/eos.hpp, has two constructors, and both have explicit.

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

  /// @brief Construct the EOS from the specific heat ratio (gamma)
  explicit IdealEOS(Real gm) : gm(gm) {}
  // ...
};

The second one, explicit IdealEOS(Real gm), creates the equation of state from the ratio of specific heats \(\gamma\). Because of explicit, accidentally passing a bare number such as 1.4f to a function that takes an equation of state is an error. : gm(gm) is the notation that initializes the member gm with the argument of the same name, gm. The name outside the parentheses in the initializer list is the member, and the name inside the parentheses is the argument.

In this way, you can write multiple constructors with the same name and different argument types (the same as function overloading in Chapter 2). IdealEOS<Real> eos(config); calls the first, and IdealEOS<Real> eos(Real(1.4)); calls the second.

4.4 const Member Functions and this

The const at the end of int size() const specifies that “this function does not modify member variables.” Writing n_ = 8; inside size() is a compile error. Put a trailing const on functions that only read, and not on functions that modify.

This distinction determines which functions can be called from a const object or a const reference.

const Axis a(4);   // forbids changes after construction
int n = a.size();  // functions with a trailing const can be called

Even from an argument received as const Axis& (the length function above), only size() and dx(), which have a trailing const, can be called. If you forget the trailing const, even a function that only reads cannot be called from a const reference.

We compare the two kinds of functions in Time in miso/include/miso/time.hpp (excerpt).

  /// @brief update time parameters
  void update(const Real dt) {
    time += dt;  // modifies members, so no trailing const
    n_step++;
  };

  /// @brief Save time parameters to file
  void save() const {  // only reads members and writes to a file
    // ...
  }

The ; after the closing brace of update is unnecessary, but writing it is not an error.

this: a pointer to the object that made the call

Inside a member function, you could use a member variable just by writing its name. This is because C++ passes to the member function “which object it was called from” as a pointer named this. When you call a.size(), this inside size() points to a, and n_ is shorthand for this->n_.

int size() const { return this->n_; }  // same as return n_;

this is a pointer, so the object it points to is written *this. The time += dt; in Time::update is, to be exact, this->time += dt;. This property that “member variables are used through this” becomes important when you learn lambda expressions in Chapter 7 and the pitfalls on the GPU in Chapter 9.

4.5 Destructors and Cleanup (RAII)

Once you allocate an array, you must release it when you are done. The allocate and deallocate pair in Fortran can be left to a class in C++. The function that is called automatically when an object’s lifetime ends is called a destructor. Its name is ~ followed by the class name, and it has no arguments or return type.

The following is a class for explanation that allocates an array in the constructor and releases it in the destructor. It prints text at allocation and release to confirm when they are called (std::printf is a function that prints with a format, available with #include <cstdio>).

#include <cstdio>

class Buffer {
public:
  explicit Buffer(int n) : data_(new double[n]), n_(n) {
    std::printf("allocated (%d elements)\n", n_);
  }
  ~Buffer() {
    std::printf("releasing (%d elements)\n", n_);
    delete[] data_;  // pairs with new[] in the constructor
  }

private:
  double *data_;
  int n_;
};

int main() {
  std::puts("before block");
  {
    Buffer work(100);  // the constructor is called here
    std::puts("while in use");
  }  // the lifetime of work ends, and the destructor is called automatically
  std::puts("after block");
}

The output is as follows.

before block
allocated (100 elements)
while in use
releasing (100 elements)
after block

You do not need to write a line that calls the destructor yourself. It is called when execution reaches the end of the block, when you return in the middle of a function, and when you leave the function by an exception. When there are multiple objects in the same block, they are cleaned up in the reverse order of creation. If a member variable is of a class type, the member’s destructor is also called automatically after the body of the enclosing destructor finishes.

Acquiring a resource (memory, a file, the MPI environment, etc.) in the constructor and releasing it in the destructor, tying it to the object’s lifetime in this way, is called RAII (Resource Acquisition Is Initialization). It prevents forgetting to write the release, or skipping the release by an early return.

A MISO example: allocating and releasing Array3D

This is the constructor and destructor of the CPU Array3D in miso/include/miso/array3d.hpp (part of the size check is omitted).

template <typename T> class Array3D<T, backend::Host> {
private:
  T *data_ = nullptr;
  std::array<int, 3> shape_ = {-1, -1, -1};

public:
  Array3D(int nx0, int nx1, int nx2) : shape_{nx0, nx1, nx2} {
    assert(shape_[0] > 0 && shape_[1] > 0 && shape_[2] > 0);
    const std::size_t array_size = static_cast<size_t>(shape_[0]) *
                                   static_cast<size_t>(shape_[1]) *
                                   static_cast<size_t>(shape_[2]);
    // ... (checks that array_size fits in the range of int)
    data_ = new T[array_size];
  }

  ~Array3D() { delete[] data_; }
  // ...
};

T is the element type, and backend::Host is a specification that indicates “place in CPU memory” (Chapter 6).
std::array<int, 3> shape_ is a fixed-length array that holds the sizes in the three directions, and represents the not-yet-created state with -1.
The constructor allocates the elements with new T[array_size], and the destructor releases them with delete[]. This is the same form as Buffer above.
static_cast<size_t>(...) is the notation that converts an int to the large integer type size_t before the multiplication. This is because multiplying the grid counts in the three directions can exceed the range of int.

The GPU Array3D<T, backend::CUDA> uses functions that allocate and release GPU memory instead of new/delete[] (Chapter 9). From the user’s side, both are arrays that are “allocated when created and released when their lifetime ends.”

A MISO example: Env, which handles MPI initialization and finalization

The first line of the demo’s main() is as follows (demo/mhd2d_vortex/src/main.cpp).

int main(int argc, char **argv) {
  // Initialize MPI and CUDA environments
  Env env(argc, argv);
  // ...
}

miso::Env in miso/include/miso/env.hpp has mpi::Env as a member. The constructor of mpi::Env calls MPI_Init, and the destructor calls MPI_Finalize (excerpt).

struct Env {
  bool owns_mpi = false;
  // ...
  ~Env() noexcept {
    // ...
    int is_mpi_finalized = false;
    MPI_Finalized(&is_mpi_finalized);
    if (owns_mpi && !is_mpi_finalized) {
      MPI_Finalize();
    }
  }
};

env is created at the start of main() and its lifetime ends at the end of main(), so you cannot forget to write the MPI finalization. noexcept is a declaration that “this function does not throw exceptions to the outside.” The details of MPI are explained in Chapter 8.

4.6 Prohibiting Copies, and Moves

The problem with copying, and = delete

An int variable can be copied with int b = a;. A class object can also be copied if you specify nothing, in which case the member variables are copied one by one. If you copy the Buffer from the previous section, only the address held in the member data_ is copied. The array itself is not duplicated, and the two objects point to the same array. Then the two destructors would delete[] the same array, which is an error such as an abnormal program termination.

For a class that manages an array, prohibit copying. Writing = delete after the declaration of a special member function makes that function unusable.

class Buffer {
public:
  // ... (the constructor and destructor are the same as before)
  Buffer(const Buffer &) = delete;             // forbids creating by copy
  Buffer &operator=(const Buffer &) = delete;  // forbids copy assignment
};

Buffer(const Buffer &) on the first line is a constructor that creates an object from an object of the same type, and is called the copy constructor. It is called by Buffer b = a; or Buffer b(a);. operator= on the second line is the function called by b = a; after creation (names containing operator are explained in Section 4.7). With both prohibited, Buffer b = a; becomes the compile error use of deleted function 'Buffer::Buffer(const Buffer&)'. The mistake is found at compile time rather than at run time.

Move: transferring ownership of the contents

Even when copying is prohibited, there are cases where you want to “transfer the owner of an array to another object.” One example is returning an array created inside a function as the return value. This is called a move. We add a constructor for moving to Buffer.

  Buffer(Buffer &&other) noexcept      // creates by move
      : data_(other.data_), n_(other.n_) {
    other.data_ = nullptr;  // the source object gives up the array
    other.n_ = 0;
  }

Buffer && is a reference that receives “an object whose contents may be moved.” The move constructor takes the address of other’s array and sets other’s side to nullptr. The array stays single, and only the owner is swapped. delete[] on nullptr does nothing, so the destructor of the emptied other releases nothing.

The following usage example changes the destructor to print only when it holds an array (if (data_ != nullptr)) and adds size(). std::move is available with #include <utility>.

int main() {
  Buffer a(10);
  // Buffer b = a;          // error: copying is forbidden
  Buffer c = std::move(a);  // moves a's array to c. No new allocation
  std::printf("a.size() = %d, c.size() = %d\n", a.size(), c.size());
}

allocated (10 elements)
a.size() = 0, c.size() = 10
releasing (10 elements)

std::move(a) is a specification that treats a as “an object whose contents may be moved,” and does not move anything by itself. What moves is the move constructor. After the move, a is empty, so you do not use it as an array.

When you return it as a function’s return value, a move (or the elision of the transfer itself) is used even if you do not write std::move.

A MISO example: Array3D and Fields

At the end of Array3D in miso/include/miso/array3d.hpp, there is the same set as in Buffer (excerpt; the // comments are from the original code).

  // Prohibit copy semantics
  Array3D(const Array3D &) = delete;
  Array3D &operator=(const Array3D &) = delete;

  // Allow move semantics (allow uninitialized state)
  Array3D(Array3D &&other) noexcept : data_(other.data_), shape_(other.shape_) {
    other.data_ = nullptr;
    other.shape_ = std::array<int, 3>{-1, -1, -1};
  }

After a move, Array3D is in the “not created” state, with data_ equal to nullptr and the size equal to -1. The comment allow uninitialized state refers to this state. The move assignment operator=(Array3D &&other) is also defined, and it releases the array it held before taking over other’s array.

Fields in miso/include/miso/mhd_fields.hpp is a struct that holds the nine MHD variables (density ro, velocity vx, vy, vz, magnetic field bx, by, bz, internal energy ei, and the variable for suppressing divergence ph) as Array3Ds. This is an excerpt of the CPU definition.

template <typename Real> struct Fields<Real, backend::Host> {
  Array3D<Real, backend::Host> ro, vx, vy, vz, bx, by, bz, ei, ph;

  Fields(const int i_total, const int j_total, const int k_total)
      : ro(i_total, j_total, k_total), vx(i_total, j_total, k_total),
        // ... (same from vy to ph)
        ph(i_total, j_total, k_total) {}
  // ...
  // Prohibit copy semantics
  Fields(const Fields &) = delete;
  Fields &operator=(const Fields &) = delete;

  // Allow move semantics
  Fields(Fields &&) = default;
  Fields &operator=(Fields &&) = default;
};

Fields itself does not allocate arrays directly. The constructors of the member Array3Ds each allocate, and their destructors each release. For that reason, no destructor is written in Fields. = default specifies “use this function in the standard form the compiler provides.” For a move, it becomes the processing “move the members one by one,” and the moves of the nine Array3Ds are each called.

An example in the demo of a function returning an Array3D is make_eta in demo/mhd2d_reconnection/src/main.cpp. It creates the resistivity array and returns it, initializing the Model member eta (excerpt).

inline Array3D<Real, Backend> make_eta(const Grid<Real, backend::Host> &grid,
                                       const Parameters &par) {
  // ... (computes values on the CPU and copies them to the Backend array eta)
  return eta;
}

// part of the member initializer list of the Model constructor
      eta(make_eta(grid, par)),

Array3D cannot be copied, but the return value is passed by move (or by eliding the transfer), so this style works.

::

There is a guideline for the five functions: the destructor, the copy constructor and assignment, and the move constructor and assignment. If you define any one of them yourself, you either define the rest too or state them explicitly with = delete or = default. It is explained as the “rule of five” in C.21 of the C++ Core Guidelines. Array3D releases in the destructor, prohibits copying, and defines moves, so it decides all five as this guideline says. Conversely, a class that does not directly hold resources, such as Fields, can leave cleanup to the member classes and write almost nothing itself (C.20 of the same document).

4.7 Operator Overloading and operator()

In MISO, an element of a three-dimensional array is written a(i, j, k). It looks similar to Fortran’s a(i, j, k), but in C++, a(i, j, k) means “call the member function named operator() of the object a with arguments i, j, k.” Defining a class-specific meaning for operators such as +, [], and () is called operator overloading. The name of the function you define is operator followed by the operator’s symbol, and operator(), operator[], and the operator= from the previous section are examples.

We check it with a two-dimensional array for explanation.

class Array2 {
public:
  Array2(int n0, int n1) : data_(new double[n0 * n1]), n0_(n0), n1_(n1) {}
  ~Array2() { delete[] data_; }
  Array2(const Array2 &) = delete;
  Array2 &operator=(const Array2 &) = delete;

  // writing a(i, j) calls this function
  double &operator()(int i, int j) { return data_[i * n1_ + j]; }

  int extent(int dim) const { return dim == 0 ? n0_ : n1_; }

private:
  double *data_;
  int n0_, n1_;
};

Array2 a(2, 3);
for (int i = 0; i < a.extent(0); ++i) {
  for (int j = 0; j < a.extent(1); ++j) {
    a(i, j) = 10 * i + j;  // the return value is a reference, so an element can be assigned to
  }
}
std::cout << a(1, 2) << "\n";  // 12
double &e = a(0, 1);           // alias of element (0, 1)
e = -1.0;
std::cout << a(0, 1) << "\n";  // -1

The return type of operator(), double &, is important. If it were defined to return a value (double), a(i, j) would become a temporary value that is a copy of the element’s value, and a(i, j) = ... would be a compile error. Because it returns a reference, a(i, j) becomes another name for the array’s element itself, and it can be both read and written. dim == 0 ? n0_ : n1_ is an expression whose value is n0_ if dim == 0 and n1_ otherwise (the conditional operator).

data_[i * n1_ + j] is an expression that converts the two-dimensional index into a one-dimensional position. When j increases by 1, you move to the adjacent element, and when i increases by 1, you advance n1_ elements, so the last index is contiguous. In Fortran the first index is contiguous, so the ordering is reversed. The difference in ordering and the loop order are explained in detail in Chapter 7.

A MISO example: operator() of Array3D

These are operator() and operator[] of the CPU Array3D in miso/include/miso/array3d.hpp.

  /// @brief Return reference to the element at the given indices.
  T &operator()(int i0, int i1, int i2) const noexcept {
    assert(data_);
    assert(i0 >= 0 && i0 < shape_[0]);
    assert(i1 >= 0 && i1 < shape_[1]);
    assert(i2 >= 0 && i2 < shape_[2]);
    return data_[(i0 * shape_[1] + i1) * shape_[2] + i2];
  }

  /// @brief Return reference to the element at the given linear index.
  T &operator[](int idx) const noexcept {
    assert(data_);
    assert(idx >= 0 && idx < size());
    return data_[idx];
  }

(i0 * shape_[1] + i1) * shape_[2] + i2 is the Array2 expression extended to three dimensions, and the last index i2 is contiguous.
assert(...) is a check that stops the program if the condition does not hold. It stops if an index is out of range. However, assert works only in Debug builds and is removed in MISO’s default Release build (Chapter 1).
operator[] specifies an element by the one-dimensional position idx. You can write a[n]. It is used when processing all elements in the same way.

Config in miso/include/miso/config.hpp also has operator[], so you can write config["time"].

  /// @brief Accessor for YAML configuration object
  const YAML::Node operator[](const std::string &key) const {
    return yaml_obj[key];
  }

It extracts and returns the item with key key from the member yaml_obj (a yaml-cpp object). The returned YAML::Node also has [], so you can chain them as config["time"]["tend"]. Because of the trailing const, it can also be called from a configuration received as const Config & (such as the argument of the Time constructor).

4.8 Owning Arrays and Non-Owning Views

The idea of ownership and views

Array2 and Array3D allocate arrays and release them when their lifetime ends. Being responsible for cleaning up a resource in this way is called “owning.” Owning objects prohibit copying, so you pass them to functions by reference.

On the other hand, all a computation function needs is the array’s start address and size. So we prepare a lightweight class that has only the start address and size, and neither allocates nor releases, and pass that to functions. This is called a View.

// non-owning View: holds only the address of the first element and the size
struct Array2View {
  double *data;
  int n0, n1;
  double &operator()(int i, int j) const { return data[i * n1 + j]; }
};

// owning array: handles allocation and release
class Array2 {
public:
  // ... (the constructor, destructor, and copy prohibition are the same as before)
  Array2View view() { return Array2View{data_, n0_, n1_}; }
  // ...
};

// takes the View by value. Only the address and the size are copied
void fill(Array2View v, double value) {
  for (int i = 0; i < v.n0; ++i) {
    for (int j = 0; j < v.n1; ++j) {
      v(i, j) = value;  // modifies elements of the original Array2
    }
  }
}

Array2 a(2, 3);
fill(a.view(), 7.0);
Array2View v = a.view();
std::cout << v(1, 2) << " " << sizeof(Array2View) << "\n";  // 7 16

fill receives an Array2View by value, so v inside the function is a copy. What was copied is only the address and size, so v(i, j) refers to the element of the original a, and a is rewritten. sizeof(Array2View) is the size of this object (in bytes). In this check it is 16 bytes (an 8-byte address and two ints), and it does not depend on the number of elements in the array.

A View does not own the array, so you must not use it longer than its owner lives. If you use a View after the lifetime of the owning Array2 has ended, it points to freed memory (the same problem as the dangling reference in Chapter 3).

A MISO example: Array3DView

This is an excerpt of Array3DView in miso/include/miso/array3d.hpp.

/// @brief Lightweight non-owning view of 3D array data.
template <typename T> class Array3DView {
private:
  T *data_ = nullptr;
  int shape_[3] = {-1, -1, -1};

  // Always constructed from Array3D
  // to ensure that the array size is within the numeric limits of int.
  template <typename, typename> friend class Array3D;
  Array3DView(T *data, int nx0, int nx1, int nx2) noexcept
      : data_(data), shape_{nx0, nx1, nx2} {}

public:
  // ... (data(), extent(), size(), operator[])
  __host__ __device__ T &operator()(int i0, int i1, int i2) const noexcept {
    // ... (check with assert)
    return data_[(i0 * shape_[1] + i1) * shape_[2] + i2];
  }

  // Allow copy semantics
  __host__ __device__ Array3DView(const Array3DView &) noexcept = default;
  __host__ __device__ Array3DView &
  operator=(const Array3DView &) noexcept = default;
};

It holds only the start address data_ and the size shape_. It does not release in its destructor.
Copying is allowed with = default. This is in contrast to Array3D, which prohibits copying. A View is meant to be copied cheaply and passed to functions.
The constructor is private. friend class Array3D is a specification that “allows Array3D to use private members,” so only view() of Array3D can create a View. As the original comment says, this is so that Views can be created only from an Array3D whose size check is done.
__host__ __device__ is a mark that makes this function usable on the GPU as well. In CPU-only builds, miso/include/miso/cuda_compat.hpp defines it as an empty mark, so you can ignore it when reading (Chapter 9).

The Array3D side has three functions that return a View (excerpt).

  /// @brief Return a lightweight view of the array.
  Array3DView<T> view() noexcept { /* ... */ }

  /// @brief Return a constant lightweight view of the array.
  Array3DView<const T> view() const noexcept { /* ... */ }

  /// @brief Return a constant lightweight view of the array.
  Array3DView<const T> const_view() const noexcept { return view(); }

Array3DView<T> is a View whose elements can be modified, and Array3DView<const T> is a View that only reads elements. Because the element type has const, the reference returned by operator() is also const T &, and assigning to it is a compile error. The first and second have the same name, and differ only in the trailing const. The first is called from a non-const Array3D, and the second from a const Array3D (or a const reference). const_view() is a function for explicitly obtaining a read-only View even from a non-const array.

Checking by running: Array3D and Array3DView

We check with a short program that uses only the MISO headers. miso/include/miso/array3d.hpp does not use MPI or yaml-cpp, so you can build it just by specifying the include path.

#include <iostream>
#include <utility>

#include <miso/array3d.hpp>

using Real = float;

// takes a writable View by value and modifies elements of the original array
void fill(miso::Array3DView<Real> a) {
  for (int i = 0; i < a.extent(0); ++i) {
    for (int j = 0; j < a.extent(1); ++j) {
      for (int k = 0; k < a.extent(2); ++k) {
        a(i, j, k) = 100 * i + 10 * j + k;
      }
    }
  }
}

// takes a read-only View and returns the sum
Real total(miso::Array3DView<const Real> a) {
  Real s = 0;
  for (int n = 0; n < a.size(); ++n) {
    s += a[n];
    // a[n] = 0;  // error: cannot assign to an element of const Real
  }
  return s;
}

int main() {
  miso::Array3D<Real> a(2, 3, 4);  // allocates a 2 x 3 x 4 array
  fill(a.view());
  std::cout << a(1, 2, 3) << " " << a[23] << "\n";  // the same element
  std::cout << total(a.const_view()) << "\n";

  // miso::Array3D<Real> b = a;          // error: copying is forbidden
  miso::Array3D<Real> b = std::move(a);  // moves ownership of the array from a to b
  std::cout << a.extent(0) << " " << b.extent(0) << "\n";
}

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

123 123
1476
-1 2

If you omit the second type specification, as in miso::Array3D<Real>, it becomes backend::Host (CPU).
a(1, 2, 3) is at position (1 * 3 + 2) * 4 + 3 = 23 in the one-dimensional layout, so it is the same element as a[23].
After the move, a is in the “not created” state with size -1, and b holds a 2 x 3 x 4 array.

operator() of Array3D returns T & even though it has a trailing const. Therefore, even from an array received as const Array3D<Real> &, you can modify elements with a(i, j, k) = ... (it compiles if you try). The const on the Array3D side guarantees only that “the array’s size and owner do not change.” When you want to show in the type that a function only reads the array, pass an Array3DView<const Real> made with const_view().

A MISO example: Fields and FieldsView

In miso/include/miso/mhd_fields.hpp, Fields, which bundles the nine variables, also has a corresponding View, FieldsView (excerpt).

template <typename T>  // Use T to represent Real or const Real
struct FieldsView {
  Array3DView<T> ro, vx, vy, vz, bx, by, bz, ei, ph;
  // ... (explicit constructor taking 9 Views, extent(), size())
};

view() of Fields bundles the view()s of the nine Array3Ds to make a FieldsView<Real>, and const_view() makes a FieldsView<const Real>. As the original comment says, T is either Real (modifiable) or const Real (read-only).

The user of FieldsView is the demo’s initial condition. MHD::apply_initial_condition in miso/include/miso/mhd_solver.hpp creates a working Fields on the CPU and passes its View to the initial condition’s apply (excerpt).

    Fields<Real, backend::Host> qq_h(grid.i_total, grid.j_total, grid.k_total);
    ic.apply(qq_h.view(), grid_h.const_view());
    qq.copy_from(qq_h);

qq_h owns the arrays, and apply receives the View by value. When apply writes values through the View, the arrays of qq_h are rewritten. After that, copy_from copies them to qq, which is used for computation (when computing on the GPU, they are transferred to GPU memory here). grid_h.const_view() is a View GridView<const Real> that only reads the grid coordinates and so on.

4.9 Reading the demo’s InitialCondition

Using what you learned in this chapter, we read InitialCondition in demo/hd1d_shock_tube/src/main.cpp. It is a struct that reads the left and right states (density, pressure, velocity) of the shock tube problem from the configuration file and gives the initial condition (part of apply is omitted).

struct InitialCondition {
  eos::IdealEOS<Real> &eos;
  Real rol, prl, vvl;
  Real ror, prr, vvr;

  explicit InitialCondition(Config &config, eos::IdealEOS<Real> &eos)
      : eos(eos), rol(config["shock_tube"]["rol"].as<Real>()),
        prl(config["shock_tube"]["prl"].as<Real>()),
        vvl(config["shock_tube"]["vvl"].as<Real>()),
        ror(config["shock_tube"]["ror"].as<Real>()),
        prr(config["shock_tube"]["prr"].as<Real>()),
        vvr(config["shock_tube"]["vvr"].as<Real>()) {}

  // The signature must not be changed as it is called inside miso::mhd::MHD.
  void apply(mhd::FieldsView<Real> qq, GridView<const Real> grid) const {

    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) {
          // ... (sets velocity, magnetic field, etc. to 0 and selects the coordinates xyz)
          if (xyz < 0.5) {
            qq.ro(i, j, k) = rol;
            qq.ei(i, j, k) = prl / (eos.gm - 1.0) / qq.ro(i, j, k);
            qq.vx(i, j, k) = vvl;
          } else {
            // ... (right-hand state)
          }
        }
      }
    }
  }
};

Part	How to read it
`struct InitialCondition`	No access is written, so all members are public
`eos::IdealEOS<Real> &eos;`	A member that holds the equation of state by reference. It holds no copy and refers to the `eos` held by `Model`
`Real rol, prl, vvl;` and so on	Members that hold the left and right density, pressure, and velocity by value
`explicit InitialCondition(Config &config, eos::IdealEOS<Real> &eos)`	A constructor that receives the configuration and the equation of state
`: eos(eos), rol(...), ...`	The reference member `eos` can have its referent set only in the initializer list. The other members are initialized in declaration order with values from the configuration file
`void apply(...) const`	A function that writes the initial condition. The trailing `const` means it does not modify `InitialCondition`’s own members (such as `rol`)
`mhd::FieldsView<Real> qq`	Receives a writable View by value. `qq.ro(i, j, k) = rol;` rewrites the arrays of the caller’s `Fields`
`GridView<const Real> grid`	A read-only View of the grid. `grid.i_total` is the number of grid points including ghost cells, and `grid.x[i]` is the coordinate of the cell center
`eos.gm`	The ratio of specific heats of the equation of state being referenced

The trailing const of apply and writing through qq do not contradict each other. What const prohibits is changing the member variables of InitialCondition, and qq is a View received as an argument. The reason the form of apply’s arguments is fixed is that MISO’s internals (MHD::apply_initial_condition in the previous section) call it in this form. The original comment says the same.

This struct is created as the member ic of the demo’s Model, and is used at the start of MISO’s time evolution. How to read Model is explained in Chapter 5.

Exercises

Exercise 1: Where each member of Time gets its value

Read the constructor of Time in miso/include/miso/time.hpp and classify, for each of the eight member variables (time, tend, dt_output, n_step, n_output, n_output_digits, time_save_dir, io_enabled), whether it gets its value in (a) the initializer list or (b) the constructor body (including functions called from the body).

::

Member	Where it gets its value
`tend`, `dt_output`, `n_output_digits`, `io_enabled`	(a) The initializer list. It is initialized with a value from the configuration file
`time`, `n_step`, `n_output`	(b) `initialize()`, called from the body, assigns 0
`time_save_dir`	(b) The body assigns a path joining `config.save_dir` and `io.time_save_dir` of the configuration file

The members in the initializer list are in the same order as the declarations. The members assigned in the body have been created once in the default way before entering the body (Real and int have indeterminate values, and std::string is an empty string). tend and dt_output are used in the first check in the body, so their values must be set in the initializer list.

Exercise 2: The position of an Array3D element

For miso::Array3D<float> a(2, 3, 4);, find the following.

Which element from the start of the array is a(1, 0, 2) (the n of a[n])?
Answer how many elements apart each of these pairs is: a(0, 0, 1) and a(0, 0, 0), a(0, 1, 0) and a(0, 0, 0), and a(1, 0, 0) and a(0, 0, 0).

::

Substitute shape_ = {2, 3, 4} into the expression of operator(), (i0 * shape_[1] + i1) * shape_[2] + i2.

(1 * 3 + 0) * 4 + 2 = 14, so it is a[14].
When the last index increases by 1, the distance is 1 element; when the middle index increases by 1, it is 4 elements (shape_[2]); and when the first index increases by 1, it is 12 elements (shape_[1] * shape_[2]).

The last index i2 is contiguous. In Fortran’s a(2, 3, 4), the first index is contiguous, so it is the reverse. The relationship with loop order is covered in Chapter 7.

Exercise 3: The argument type of a function that receives an array

Suppose you write a function max_value that finds the maximum of the elements of an Array3D<Real>, and a function clear that sets all elements to 0. Among the candidates below, choose the suitable argument type for each function and explain why.

(a) Array3D<Real> a
(b) const Array3D<Real> &a
(c) Array3DView<Real> a
(d) Array3DView<const Real> a

::

(d) Array3DView<const Real> is suitable for max_value. The type shows that it only reads, and accidentally assigning inside the function becomes a compile error. The caller writes max_value(a.const_view()).
(c) Array3DView<Real> is suitable for clear. Even if you pass a View by value, what it points to is the original array, so the original elements can be set to 0. The caller writes clear(a.view()).

(a) tries to make a copy of Array3D, so it is a compile error because copying is prohibited. (b) compiles and does not copy. However, operator() of Array3D returns T & even with a trailing const, so elements can be modified even through const Array3D<Real> &. It is inferior to (d) in that the type cannot guarantee that it only reads. Also, inside parallel loops such as MISO’s for_each, Views are used (Chapter 7). Making the function’s argument a View matches that style.