Testing and Comparing Results

This page collects ways to confirm that you have not broken anything when you change MISO. The first half covers running and reading MISO’s unit tests (doctest and ctest); the second half covers how to compare demo output and a rule of thumb for float rounding error. It assumes the build procedure of Chapter 1.

What to compare
Running the unit tests
Reading a test
Comparing demo output
Rule of thumb for float rounding error
Recording the comparison conditions

What to compare

The degree of agreement you expect depends on what you compare.

Comparison	What it checks	Expected agreement
The result of a component versus the expected value from a formula or hand calculation	Indices, coefficients, and boundary handling	Exact for integers. About the rounding error for real numbers
Two runs with the same executable and the same settings	That the result does not change from run to run	Bit-for-bit identical
Before and after a change (a change that should not alter the computation)	That the rewrite did not change the computation	Bit-for-bit identical, or about the rounding error
Runs with a different number or direction of MPI decomposition	Domain decomposition and communication	Bit-for-bit identical as a whole array, or about the rounding error
The CPU and GPU versions, and differences in compiler or optimization	The effect of differences in the order of operations or in instructions	About the rounding error (it grows with time evolution)
Numerical and exact solutions	The discretization (grid and time step) error	That it becomes smaller as the grid is refined

In a MISO PR, describe the effect on numerical results as one of “bit-for-bit identical”, “about the rounding error”, “discretization error”, or “physically changed” (MISO’s dev/notes/rules_ai.md). The methods on this page give you the material to judge which one applies.

Compute the expected value from a formula or a separate simple computation, without using the function you want to check as it is. If you build the expected value with the same function, the function’s mistake also enters the expected value, and the comparison cannot find it.

Running the unit tests

Running them all together (ctest)

MISO’s unit tests are in miso/tests/, and they are built together when you build MISO by itself (ENABLE_TESTING is ON by default). The build folder may be outside the MISO folder.

cmake -S /path/to/MISO/miso -B build_miso -DMISO_USE_CUDA=OFF
cmake --build build_miso --parallel 4
cd build_miso
OMPI_MCA_btl_vader_single_copy_mechanism=none \
OMPI_MCA_rmaps_base_oversubscribe=yes \
ctest --output-on-failure

      Start  1: test_array1d_cpp
 1/15 Test  #1: test_array1d_cpp .................   Passed    0.01 sec
...
15/15 Test #15: test_parallel_grid_cpp ...........   Passed    2.08 sec

100% tests passed, 0 tests failed out of 15

The number and names of the tests change with the MISO version. --output-on-failure shows only the output of failed tests. It does not show the output of passing tests.

An environment variable written in the form NAME=value command is in effect only during that command. The two environment variables change the behavior of Open MPI.

Environment variable	Meaning
`OMPI_MCA_btl_vader_single_copy_mechanism=none`	Avoids a memory-copy mechanism for communication between processes on the same machine, which may be unavailable in containers and the like
`OMPI_MCA_rmaps_base_oversubscribe=yes`	Allows launching more processes than there are CPU cores. This lets you run the 2-process tests even on a machine with few cores

It may work on your own machine without them, but MISO’s procedure (dev/notes/usage_ai.md) and GitHub CI run in this form. Multi-process tests are also launched through miso/shared/mpiexec_with_env.sh, which also loads the settings in miso/shared/openmpi_config.env.

Selecting tests

Command	What it does
`ctest -N`	Lists the registered tests without running them
`ctest -R array`	Runs only the tests whose names contain `array` (a regular expression)
`ctest -R execution -V`	Runs the selected tests, showing the run command and all of the output

The name of a test registered with ctest is the source file name followed by _cpp (_cu for the GPU version). miso/tests/serial/CMakeLists.txt and miso/tests/parallel/CMakeLists.txt collect the files matching test_*.cpp and register them automatically. If you put a new test file in serial/, it is registered the next time you run the configuration.

Running a test executable directly

Each test is an independent executable, located in build_miso/tests/serial/ and build_miso/tests/parallel/. If you run it directly, you can narrow down what runs further with doctest options.

./tests/serial/test_execution_cpp --list-test-cases
./tests/serial/test_execution_cpp --test-case="Test reduce 1D CPU" --success

Option	Meaning
`--list-test-cases`	Lists the names of the test cases contained in the file
`--test-case=name`	Runs only the test cases whose names match. You can use `` (example: `"Test reduce"`)
`--success`	Also shows the successful checks
`--no-colors`	Suppresses the control characters for text color. Use it when you save the log to a file

The tests in parallel/ run with multiple MPI processes. To launch one in the same form as ctest, do the following.

/path/to/MISO/miso/shared/mpiexec_with_env.sh -n 2 ./tests/parallel/test_parallel_grid_cpp

Each process prints a summary of its own results, so the summary appears twice. Check the summaries of all the processes.

The GPU tests (*.cu) are built only with -DMISO_USE_CUDA=ON. See Chapter 9.

Reading a test

MISO’s unit tests are written with a testing library called doctest. doctest consists of a single header, and FetchContent in miso/tests/CMakeLists.txt fetches it from GitHub at configuration time.

Let us read the first test in miso/tests/serial/test_execution.cpp (excerpt).

#define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN
#include <doctest/doctest.h>

#include <miso/array1d.hpp>
#include <miso/array3d.hpp>
#include <miso/execution.hpp>

using namespace miso;

TEST_CASE("Test for_each 1D CPU" * doctest::test_suite("execution")) {
  Range1D range{2, 4};
  Array1D<int, backend::Host> arr(5);

  auto view = arr.view();
  for_each(backend::Host{}, range, [&](int i) { view[i] = i * i; });

  // This also works, but users should use the view interface as API.
  // for_each<backend::Host>(range, [&](int i) { arr[i] = i * i; });

  for (int i = 2; i < 4; ++i)
    REQUIRE(view[i] == i * i);
}

Syntax	Meaning
`#define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN`	This tells doctest to create the `main` function. You do not write `main` in the test file
`TEST_CASE("name" * doctest::test_suite("suite"))`	Defines one test case. What is inside `{ }` is the processing executed by that test. `test_suite` is the name of the suite that groups test cases
`REQUIRE(condition)`	If the condition does not hold, it counts as a failure and the rest of that test case is not executed
`CHECK(condition)`	If the condition does not hold, it is recorded as a failure, and execution continues

This test runs a for_each that puts i * i only into indices 2 and 3, and checks those two values. It is an integer computation, so it requires an exact match with ==. for_each, Range1D, Array1D, and lambda expressions are explained in Chapter 7.

Comparing real numbers: `doctest::Approx`

In real-number computations, even mathematically equivalent expressions can differ in the last digits depending on the order of computation (a later section). For this reason, to compare real numbers, use doctest::Approx, not ==, to check that they are “close enough”.

The following is an example that checks, in two ways, the result of adding 0.1 ten times in float. The doctest header is in _deps/doctest-src/ in the folder where you built MISO’s tests.

#define DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN
#include <doctest/doctest.h>

TEST_CASE("Sum of 0.1f" * doctest::test_suite("example")) {
  float s = 0.0f;
  for (int i = 0; i < 10; ++i) {
    s += 0.1f;  // add 0.1 ten times
  }
  CHECK(s == doctest::Approx(1.0));  // passes if close to within about 1e-5 relative
  CHECK(s == 1.0f);                  // passes only if bitwise identical
}

g++ -std=c++17 -I/path/to/build_miso/_deps/doctest-src test_sum.cpp -o test_sum
./test_sum --no-colors

test_sum.cpp:10: ERROR: CHECK( s == 1.0f ) is NOT correct!
  values: CHECK( 1 == 1 )

===============================================================================
[doctest] test cases: 1 | 0 passed | 1 failed | 0 skipped
[doctest] assertions: 2 | 1 passed | 1 failed |
[doctest] Status: FAILURE!

The Approx check succeeds, and the == check fails. s is 1.00000012, which differs from 1 only in the last digit. In the failure message values: CHECK( 1 == 1 ), both sides look the same, but this is because the display has few digits. When a comparison of real numbers fails, increase the number of digits yourself and print the values to check. If even one test fails, the executable returns a nonzero exit code, and ctest shows that test as Failed.

a == doctest::Approx(b) succeeds when the following condition holds (the implementation in doctest 2.4.12).

\[ |a - b| < \varepsilon \, \bigl(s + \max(|a|, |b|)\bigr) \]

You specify \(\varepsilon\) with .epsilon(value) and \(s\) with .scale(value). By default, \(\varepsilon\) is 100 times the machine epsilon of float (about \(1.2\times10^{-5}\)) and \(s = 0\), so the comparison uses the relative error. With .scale(1.0), for values near 0 the comparison uses \(\varepsilon\) as the tolerance of the absolute error. This is because, at a point where the expected value is 0, using a relative error makes any small difference a failure.

An example of an integration test in MISO

miso/tests/parallel/test_hd1d_shock_tube.cpp is a test that actually evolves a 1D shock tube (Sod’s problem) in time with MISO. It makes two kinds of comparison (excerpt).

  CHECK(ro_diff_sum / ro_sum < 5.e-2);
  CHECK(vx_diff_sum / vx_sum < 5.e-2);

  constexpr double single_precision_epsilon = 1e-5;
  for (std::size_t i = 0; i < model_x.grid.i_total; ++i) {
    CHECK(model_x.mhd.qq.ro[i] == doctest::Approx(model_y.mhd.qq.ro[i])
                                      .epsilon(single_precision_epsilon)
                                      .scale(1.0));

The first is a comparison with the exact solution. It checks that the sum of the absolute differences divided by the sum of the values is less than 5%. The difference from the exact solution includes the discretization error, so the tolerance is fairly large.

The second is a comparison between computations of the same problem placed along the x, y, and z directions. Mathematically these are the same computation, so it requires agreement to about the rounding error (\(\varepsilon = 10^{-5}\)).

Comparing demo output

To check whether a change to MISO has altered the computed results, run the demo before and after the change and compare the output.

Runs for comparison

The output destination of a demo, io.save_dir (data/ by default), is a path relative to the folder containing config.yaml. So if you create a separate folder for each run, copy config.yaml into it, and specify that file with --config, you can keep the outputs separate. This is an example of running twice, assuming demo/mhd2d_vortex has been built.

mkdir -p "$HOME/miso-work/compare/run_a" "$HOME/miso-work/compare/run_b"
cd "$HOME/miso-work/compare"
cp /path/to/MISO/demo/mhd2d_vortex/config.yaml run_a/
cp /path/to/MISO/demo/mhd2d_vortex/config.yaml run_b/
set -a
source /path/to/MISO/demo/shared/openmpi_config.env
set +a
mpiexec -np 2 /path/to/MISO/demo/mhd2d_vortex/build/mhd2d_vortex --config=run_a/config.yaml
mpiexec -np 2 /path/to/MISO/demo/mhd2d_vortex/build/mhd2d_vortex --config=run_b/config.yaml

From set -a to set +a, you make the environment variables written in openmpi_config.env take effect in this terminal (app_run.sh does the same thing).

The demo’s config.yaml has io.continue: true, so if results already exist in the output destination, the run continues from them. Do each run for comparison in an empty folder.

When comparing before and after a change, use the executables built from the respective versions of MISO. If you compare within the same MISO folder, either copy the pre-change executable under another name, or use separate build folders.

Bit-for-bit comparison (`cmp`, `diff -rq`)

Check whether two files differ by even one byte with cmp.

cmp run_a/data/mhd/mhd.00000080.00000000.bin run_b/data/mhd/mhd.00000080.00000000.bin

If they are the same it prints nothing; if they differ it prints the position of the first differing byte. You can compare whole folders with diff -rq.

diff -rq run_a/data run_b/data

If nothing is printed, all the outputs are bit-for-bit identical.

Bit-for-bit comparison can be used only when the layout of the output is the same. If you change the number or direction of the MPI decomposition, the range covered by each rank’s file (mhd.<n_output>.<rank>.bin) changes, so the files do not match even when the values inside are the same. In that case, assemble them into a whole array as follows and then compare.

Comparing as a whole array (pymiso)

pymiso.Data reads the file of each rank and assembles them into a whole array (the ghost cells are removed). The following script takes the last output numbers of the two outputs and, for each variable, prints the maximum difference divided by a representative value (the maximum absolute value of that variable). Some points have an expected value of 0, so we divide by the representative value, not use a pointwise relative error.

import sys

import numpy as np

import pymiso

dir_a, dir_b = sys.argv[1], sys.argv[2]
a = pymiso.Data(data_dir=dir_a)
b = pymiso.Data(data_dir=dir_b)
n = min(a.n_output, b.n_output)  # last output present in both
a.load(n)
b.load(n)
print(f"n_output = {n}, time = {a.time.time:.6g}, {b.time.time:.6g}")

for var in ["ro", "vx", "vy", "vz", "bx", "by", "bz", "ei"]:
    qa = getattr(a, var).astype(np.float64)
    qb = getattr(b, var).astype(np.float64)
    scale = np.max(np.abs(qa))  # divide by a representative value (the maximum magnitude)
    if scale == 0.0:
        print(f"{var}: variable is 0 in both")
        continue
    err = np.max(np.abs(qa - qb)) / scale
    print(f"{var}: max|a-b|/max|a| = {err:.3e}")

Save it as compare.py and run it in the virtual environment where you installed pymiso in Chapter 1.

python compare.py run_a/data run_b/data

a.time is an object that handles the time file, and the time value is a.time.time. Convert the differences to float64 before computing. If you subtract in float32, rounding also enters the difference itself. vz and bz are all 0 in this 2D problem, so they cannot be divided by the representative value.

Examples of comparison

These are the results of comparing under different conditions with the default settings of demo/mhd2d_vortex (128×128, 354 steps up to time 0.8, Real is float) (Ubuntu 24.04, GCC 13.3, Open MPI 4.1).

What was compared	`diff -rq`	`compare.py` (time 0.8)
Two runs with the same executable	No difference	All 0
Executables built as Release and Debug	No difference	All 0
Changing the MPI decomposition from 2 in the x direction to 2 in the y direction	All `mhd.*.bin` differ (and `mpi/coords.csv` too)	All 0
Executable built with `-DCMAKE_CXX_FLAGS=-march=native`	`mhd.*.bin` other than the initial values (output number 0), and the files of some times, differ	\(1.7\times10^{-6}\) to \(5.7\times10^{-6}\) depending on the variable

The computation with a changed decomposition differed as files, but was bit-for-bit identical as a whole array. -march=native tells the compiler to use instructions specific to the CPU in use (such as FMA, which merges multiply and add). Because the rounding of operations changes, the result changes by about the rounding error. The time step determined by the CFL condition also changes in its last digits, and some output times differed in the 7th digit. The change of the difference with time in this comparison was as follows.

Output number	Time	Difference in `ro` (the value from `compare.py`)
1	0.01	\(5.9\times10^{-7}\)
10	0.10	\(9.3\times10^{-7}\)
40	0.40	\(1.6\times10^{-6}\)
80	0.80	\(1.7\times10^{-6}\)

Whether Release and Debug agree depends on the compiler and the machine. The result above is one example.

Rule of thumb for float rounding error

MISO’s demos and tests compute with using Real = float; (single precision). Real numbers are represented with a finite number of digits, so each computation rounds the result to a representable value. The difference between 1 and the next larger representable number is called the machine epsilon.

Type	Machine epsilon	Approximate significant decimal digits
`float` (single precision)	\(2^{-23} \approx 1.19\times10^{-7}\)	7 digits
`double` (double precision)	\(2^{-52} \approx 2.22\times10^{-16}\)	16 digits

The rounding error of one arithmetic operation is at most half the machine epsilon in relative terms. The result of adding 0.1f ten times in the example above is 1.00000012, and its difference from 1, \(1.19\times10^{-7}\), is one machine epsilon of float.

The following is a rule of thumb for judging a difference.

Relative difference	How to judge it (for `float`)
0	Bit-for-bit identical. The order of computation has not changed either
About \(10^{-7}\) to \(10^{-5}\)	About the rounding error. It arises from differences in the order of operations or in instructions (the order of addition, FMA, the implementation of math functions, the difference between CPU and GPU)
Larger than that	The content of the computation may have changed. Check the input, the grid, the time step, the boundary conditions, and the expressions you changed

This is a rule of thumb, not a fixed tolerance. Rounding errors accumulate at each time-evolution step, and in a flow where small differences grow, such as turbulence, the difference can become large after a long time. In such a comparison, check that the difference is about the rounding error at an early time, and also look at how the difference grows with time. In an expression with cancellation (the loss of significant digits in a subtraction of nearby values), the difference can become larger than the machine epsilon.

Recording the comparison conditions

So that you can redo the same comparison later, record the following.

The MISO version (the output of git log -1 --oneline) and the content of the change
The compiler and its version, the build type, and whether it is the CPU or the GPU version
The demo and config.yaml you used, and the number of MPI processes and the decomposition
The output numbers and variables compared, the comparison method (cmp, compare.py, and so on), and the values obtained

When you open a PR for MISO, write the commands you ran and the numbers you compared in the body of the PR (MISO’s dev/notes/rules_ai.md).