Samurai Euler Tutorial

This project demonstrates how to solve the Euler equations using the Samurai library with adaptive mesh refinement (AMR).

Prerequisites

You need a package manager like conda or mamba to install the dependencies.

Installation

1. Create the Conda environment:

   Use the provided environment.yml file to create the environment.

   mamba env create -f conda/environment.yml

   Or with conda:

   conda env create -f conda/environment.yml

2. Activate the environment:

   mamba activate samurai-euler-env

Building the Project

1. Create a build directory:

   mkdir build
   cd build

2. Configure the project with CMake:

   cmake .. -DCMAKE_BUILD_TYPE=Release

3. Build the executables:

   make

Running the Tests

Euler 2D

To run the 2D Euler simulation:

./euler_2d

This will generate output files (e.g., HDF5/XDMF) in the results directory (or current directory depending on configuration), which can be visualized using ParaView.

Other Executables

• euler_1d: 1D Euler simulation.
• euler_user_pred_1d: 1D Euler simulation with user-defined prediction.

Command Line Options

The euler_2d executable accepts several command-line arguments to control the simulation, multiresolution parameters, and output.

Simulation Parameters

Option	Description	Default
--cfl	The CFL number	0.4
--Ti	Initial time	0.0
--Tf	Final time	0.25
--scheme	Riemann solver (rusanov, hll, hllc)	hllc
--order	Order in space: 1 on cell averages, 2 on a MUSCL reconstruction	1
--slope-limiter	Slope limiter of the reconstruction (minmod, vanleer, moncen, none)	moncen
--time-integrator	Time stepping (auto, euler, ssprk2, strang)	auto
--test-case	Test case to run	double_mach_reflection
--gamma	Ratio of specific heats; overrides the value of the test case	(test case)
--restart-file	Path to a file to restart the simulation from	(empty)
--check-positivity	Check positivity of density and pressure at each iteration	off

Run ./euler_2d --help for the list of available test cases: it is read from the registry, so it always matches what the binary actually supports. euler_3d accepts the same options and exposes the cases that are defined in 3D.

Test case parameters

A case that comes in variants declares its own option, and --help lists them under Test case parameters. There are two today, both belonging to lax_liu:

Option	Description	Default
--riemann-config	Lax & Liu configuration, 1 to 19	3
--riemann-interface	Where the quadrants meet, in each direction	0.8

lax_liu is the four-quadrant Riemann problem in two dimensions and the eight octant one in three. In 3D only configuration 3 exists, the one the article extends, and --riemann-config there accepts nothing else.

The interface position is not part of the classification. Lax & Liu and Kurganov & Tadmor put the four quadrants of the unit square at its centre; the default here is 0.8, which is what the article uses for configuration 3 so that the waves fill the domain by t_f = 0.8 without reaching the boundary. It suits configurations 3 and 4 and misleads for the other seventeen, whose structure is born 0.2 from a corner and reaches it early: pass --riemann-interface 0.5 to compare those against their published figures.

The options of every case are declared, not only those of the selected one: which case runs is itself decided by the parse. Two cases must therefore not ask for the same option name, and naming the option after the case is what keeps them apart.

Order and time stepping

--scheme and --order are independent: the first says which Riemann solver settles an interface, the second what states it is handed. At order 1 those are the two cell averages; at order 2 they are values reconstructed on the face from a limited slope, on a stencil of four cells.

--time-integrator says how the step is taken.

Value	What it does
euler	One explicit step. At order 2 the flux also advances the face values half a step, which makes it MUSCL-Hancock.
ssprk2	Two stages averaged (Heun).
strang	One directional sweep at a time, X(dt/2) Y(dt) X(dt/2) in two dimensions.
auto	Explicit Euler at order 1, Strang at order 2.

The three are not interchangeable at order 2 in more than one dimension. The Hancock predictor advances a face value with the equations taken normal to that face, and the transverse terms it leaves out are of the same order as the ones it keeps, so unsplit it converges at 1.05 on the isentropic vortex. Under a directional sweep there is no transverse direction and the same predictor is exact, which is why Strang reaches 2.23 there, and SSP-RK2 2.13. In one dimension Strang is the single Hancock step, to the bit.

--slope-limiter none is unlimited and oscillates at a shock. It is there to measure an order on a smooth solution, where every limiter clips the extremum and costs a fraction of an order that says nothing about the scheme.

Multiresolution Parameters

Option	Description	Default
--min-level	Minimum level of the multiresolution	8
--max-level	Maximum level of the multiresolution	8

Output Parameters

Option	Description	Default
--path	Output directory path	results
--filename	Output file name prefix	<test-case>_<scheme>
--nfiles	Number of output files to generate	1
--metrics-file	Write the performance metrics of the run, as JSON, to this file	(none)

Example Usage

Run the simulation with a specific final time and output 10 files:

./euler_2d --Tf 0.5 --nfiles 10

Run with adaptive mesh refinement (levels 5 to 10):

./euler_2d --min-level 5 --max-level 10

Run a case with a different gas, writing somewhere else:

./euler_2d --test-case sedov_blast --gamma 1.6666667 --path out --filename sedov_g53

Run the reference case of the article, configuration 3 of Lax & Liu, to its final time on an adapted mesh:

./euler_2d --test-case lax_liu --riemann-config 3 --min-level 4 --max-level 10 --Tf 0.8 --order 2

Run configuration 5 as Kurganov & Tadmor publish it, quadrants meeting at the centre:

./euler_2d --test-case lax_liu --riemann-config 5 --riemann-interface 0.5 --Tf 0.3

Performance

The article this repository reproduces is first of all a performance paper, and its tables are made of three numbers. Every run reports them:

performance
  cells            17524 -> 151480 of 262144 uniform
  sparsity index   6.68% -> 57.79%
  cell updates     262193561 over 2672 time steps
  time to solution 102.92 s, 28.4% of it adapting the mesh
  throughput       2.55 Mcu/s

Metric	What it is
sparsity index	cells of the adapted mesh over the cells of the uniform mesh at --max-level, in percent. 100% is a mesh that never coarsened. Given at the initial and at the final time, as the article gives it: a Riemann problem fills its mesh up as the waves spread, and one number taken at one end would flatter or damn it.
cell updates, Mcu/s	one cell advanced by one time step is one cell update; the throughput is millions of those per second over the whole run. Counted once per cell per step whatever the integrator does inside, so that the 2*dim - 1 sweeps of Strang do not read as more work done.
time to solution	the time loop. The mesh adaptation is part of it and is reported apart; writing files is not, --nfiles being a choice of whoever runs the solver.

--metrics-file <name> writes the same numbers as JSON, which is what python/performance.py builds a table out of:

python python/performance.py --levels 6 7 8 9 --uniform

One run per resolution, on the reference case of the article — configuration 3 of Lax & Liu, to t_f = 0.8, second order — gives, on one core:

 l_min  l_max   resolution    mr-eps    Mcu/s  time (s)    AMR           cells ti/tf    sparsity ti/tf
     3      6         64^2   default      2.2      0.39  25.1%           1720 / 3784     42.0% / 92.4%
     3      7        128^2   default      2.4      2.49  29.6%          3916 / 13522     23.9% / 82.5%
     3      8        256^2   default      2.5     15.41  28.1%          8416 / 45394     12.8% / 69.3%
     3      9        512^2   default      2.5    102.92  28.4%        17524 / 151480      6.7% / 57.8%
     9      9        512^2   default      4.7    147.64   0.0%       262144 / 262144   100.0% / 100.0%

The cell counts are reproducible to the cell; the times are wall clock on one core and move by ten percent or so between runs, which is worth remembering before reading anything into a small difference.

The last row is the uniform mesh at the same resolution, which is the reference the article puts at the bottom of its own table. Reading the two bottom rows together is the whole point of the exercise: the adapted run does 2.7 times fewer cell updates and is 1.4 times faster, because it runs at a little more than half the throughput of the uniform one. Roughly a third of what is lost is the adaptation itself, at 28% of the time to solution; the rest is what an adapted mesh costs per cell — level interfaces, prediction, intervals that are shorter than a uniform row.

Comparing that with the table of the article takes some care, and the sparsity column is where it goes wrong most easily:

• Equivalent resolution is the only fair pairing. The article varies the number of cells per octree leaf at a fixed equivalent resolution of 4096²; samurai carries one cell per leaf, so that axis does not exist here and the table above sweeps the resolution instead. A sparsity index quoted without the max-level it was measured at compares nothing: what the adaptation keeps is a neighbourhood of the discontinuities, which are curves in a plane, so their share of the mesh falls as the resolution rises — 92%, 83%, 69%, 58% over the four rows above.

• The refinement criterion is not the same one, and this is the real difference between the two codes rather than a defect of either. The article refines on a Löhner criterion, a normalised second difference thresholded at r_refine = 0.4; samurai refines on the details of the multiresolution thresholded at --mr-eps. The multiresolution comes with an error estimate that the gradient criterion has not, and it keeps more cells for it. The threshold is the knob that trades the two against each other, and --mr-eps sweeps it:

python python/performance.py --levels 9 --mr-eps 1e-4 1e-3 1e-2

 l_min  l_max   resolution    mr-eps    Mcu/s  time (s)    AMR           cells ti/tf    sparsity ti/tf
     3      9        512^2     1e-04      2.3    115.04  28.5%        17524 / 151480      6.7% / 57.8%
     3      9        512^2     1e-03      1.8     69.19  38.9%         17524 / 72616      6.7% / 27.7%
     3      9        512^2     1e-02      1.5     58.86  43.6%         17524 / 42139      6.7% / 16.1%

Two decades of threshold take the final sparsity from 58% to 16%, which is the order of magnitude the article reports, and the time to solution from 115 s to 59 s. The initial mesh does not move at all: the details of a piecewise constant state are of order one at the discontinuities and far above every threshold in this range, so the three runs start from the same cells and part company as the solution develops structure. These rows say nothing about accuracy, and a large enough threshold makes any mesh sparse and any solution wrong; the error of an adapted run against a uniform one is what python/error_analysis.py measures, on the cases that have an exact solution.

• The AMR share is not measured over the same cadence. The article runs its AMR cycle once every 10 time steps and this solver adapts at every one, which is most of the distance between the 28% above and the few percent it reports on CPU at its nominal block size.

• Throughput is architecture, not method. The numbers of the article are measured on 72 ARM cores or on a Hopper GPU, against one core here, and its solver works on blocks of 16² cells where this one works on intervals. The column worth comparing is the sparsity index; the Mcu/s column is worth comparing against itself, between two runs of this solver.

Tests

The suite drives the built binaries as subprocesses, so it checks what a user actually runs, command line included.

ctest --test-dir build --output-on-failure       # the fast tests, about 30 s
ctest --test-dir build -L slow                   # the validation runs as well

It has three tiers, and they are not interchangeable.

test_invariants.py owns no reference file. It asserts properties that stay true when the numerics legitimately change: a uniform flow stays uniform, a closed box conserves mass and energy and a periodic one conserves momentum as well, density and pressure stay positive, the Sedov blast keeps its rotational symmetry, a restart reproduces the run. A better scheme cannot make these fail, and no amount of regenerating can make them pass.

test_regression.py compares whole fields against references under tests/reference. One test there compares no field against a reference but the references against each other: three identical files for the three Riemann solvers mean the entry exercises none of them, which happens when a run is too short to leave the initial state or too coarse to resolve it. Every case there runs on a uniform mesh, on purpose: on an adapted mesh a rounding difference of the order of 1e-16 near the multiresolution threshold flips a refinement decision, the mesh changes, and the comparison fails on another compiler without anything being wrong. Regenerate the references with pytest --generate-ref, and say in the commit message why they moved.

test_validation.py is marked slow and asserts on scalars rather than fields, which is what makes it usable on an adapted mesh. It measures the convergence order of the isentropic vortex against its exact solution, checks that adaptation reaches the same error as a uniform mesh with fewer cells, holds the two shock tubes to the exact Riemann solution of python/exact_riemann.py — L1 errors and the position of each wave — and holds the six Lax & Liu configurations that are symmetric about the diagonal to that symmetry, which a single mistyped digit in one quadrant breaks.

Adding a test case

A test case is a domain, an initial state, a set of boundary conditions and the gas it is written for. Add a header in euler/init/, expose definition<Field>() in a namespace of its own, close the file with

REGISTER_TEST_CASE(my_case, test_case::my_case, 2, 3)

and add one #include to euler/init/cases.hpp. There is no list to keep in sync: the trailing numbers are the dimensions the case is written for, and a case whose definition does not depend on the dimension (free_stream, sedov_blast, sod_x) passes all three and serves euler_1d, euler_2d and euler_3d. The boundary conditions the cases need — outflow, solid wall, an imposed state — are in euler/bc.hpp and work in any dimension; a case with no boundary at all sets periodic instead, as blast_periodic does.

A case that takes a parameter of its own gives the registry an options function, which declares the command line option and keeps the variable the parser writes into. lax_liu is the worked example: the parameter is read when a cell is initialised, which is after the parse, and the values the option accepts depend on the dimension the case was registered for.

Give the case the gas it was written for through its eos field: monofluid cases use EOS::ideal_gas(gamma). euler/eos.hpp also defines a stiffened gas for a future two-phase model; the solver is templated on the state law, so the monofluid path does not pay for the coefficients it never uses.