Fix the test case constants, and check every case against its paper

.github/workflows/ci.yml

@@ -16,6 +16,12 @@ concurrency:
 jobs:
   build-and-test:
     runs-on: ubuntu-latest
+    env:
+      # Some of the Azure runners expose a MANA RDMA device. UCX picks it up,
+      # fails to open it, and MPI_Init aborts, which fails every test before it
+      # runs. The suite drives one rank, so keeping UCX on shared memory and TCP
+      # costs nothing.
+      UCX_TLS: sm,self,tcp
     steps:
       - uses: actions/checkout@v4
 

euler/init/closed_box.hpp

@@ -22,6 +22,9 @@
 //  check that a flux is conservative, which is the first thing a reconstruction
 //  can quietly break, and it is the only case that exercises the Reflective
 //  boundary condition.
+//
+//  Like the free stream, this case has no publication behind it: what it checks
+//  is a conservation identity of the scheme, held to round-off.
 // =============================================================================
 
 namespace test_case::closed_box

euler/init/double_mach.hpp

@@ -19,6 +19,19 @@
 //  wall along the bottom boundary starting at x0, the shock entering upstream of
 //  it. The top boundary follows the analytic shock position, which is why the
 //  boundary condition needs the current time.
+//
+//  The set-up is the one of
+//
+//      P. Woodward, P. Colella, "The numerical simulation of two-dimensional
+//      fluid flow with strong shocks", J. Comput. Phys. 54 (1) (1984) 115-173,
+//      https://doi.org/10.1016/0021-9991(84)90142-6
+//
+//  section IVa: domain [0,4] x [0,1], wall from x0 = 1/6, final time 0.2.
+//
+//  The post-shock state below follows from the Rankine-Hugoniot conditions for
+//  a Mach 10 shock running into (rho, p) = (1.4, 1), a gas whose sound speed is
+//  therefore 1: density 1.4 (2.4 M^2)/(0.4 M^2 + 2) = 8, pressure
+//  (2.8 M^2 - 0.4)/2.4 = 116.5, speed (1 - 1.4/8) x 10 = 8.25.
 // =============================================================================
 
 namespace test_case::double_mach_reflection

euler/init/double_rarefaction.hpp

@@ -22,7 +22,15 @@
 //  it fails loudly rather than quietly. It is the initial state euler_1d used
 //  before the test cases were shared, kept here as its default.
 //
-//  Toro, "Riemann Solvers and Numerical Methods for Fluid Dynamics", test 2.
+//  The states are test 2 of Table 4.1 of
+//
+//      E.F. Toro, "Riemann Solvers and Numerical Methods for Fluid Dynamics.
+//      A Practical Introduction", 3rd ed., Springer, 2009,
+//      https://doi.org/10.1007/b79761
+//
+//  whose exact solution python/exact_riemann.py computes: p* = 0.00189,
+//  u* = 0, and a star density of 0.02185 on both sides. That pressure is four
+//  orders of magnitude below the initial one.
 // =============================================================================
 
 namespace test_case::double_rarefaction

euler/init/free_stream.hpp

@@ -23,6 +23,10 @@
 //
 //  The exact solution at any time is the same constant state, so the error is
 //  obtained by comparing to the initial uniform values.
+//
+//  No paper to compare this against: free-stream preservation is a property of
+//  the discretisation, so it is checked against machine precision rather than
+//  against a published figure.
 // =============================================================================
 
 namespace test_case::free_stream

euler/init/isentropic_vortex.hpp

@@ -23,10 +23,19 @@
 //  case the reference for measuring convergence rates; the moving structure also
 //  exercises the mesh adaptation (the refined region must track the vortex).
 //
-//  Notation and formulas follow:
+//  The case goes back to
+//      C.-W. Shu, "Essentially non-oscillatory and weighted essentially
+//      non-oscillatory schemes for hyperbolic conservation laws", in: Advanced
+//      Numerical Approximation of Nonlinear Hyperbolic Equations, Lecture Notes
+//      in Mathematics 1697, Springer, 1998, pp. 325-432,
+//      https://doi.org/10.1007/BFb0096355
+//
+//  but the many variants of it in circulation differ in their scaling, so the
+//  notation and formulas here follow the survey that sorts them out:
 //      S. C. Spiegel, H. T. Huynh, J. R. DeBonis, "A Survey of the Isentropic
 //      Euler Vortex Problem using High-Order Methods", AIAA Paper 2015-2444,
-//      NASA Glenn Research Center, 2015 (NTRS 20150018403).
+//      NASA Glenn Research Center, 2015 (NTRS 20150018403),
+//      https://doi.org/10.2514/6.2015-2444.
 //  Equation/section numbers below refer to that paper. The parameter set is the
 //  "Shu" row of Table 1, in the paper's (sound-speed) non-dimensionalization:
 //  rho_inf = 1, a_inf = 1, T_inf = 1, R_gas = 1, so p_inf = 1/gamma.

euler/init/kelvin_helmholtz.hpp

@@ -13,9 +13,42 @@
 #include "../variables.hpp"
 #include "registry.hpp"
 
-// Kelvin-Helmholtz instability: two horizontal layers in shear. A small
-// vertical-velocity perturbation localized at the two interfaces grows into
-// the characteristic rolled-up billows.
+// =============================================================================
+//  Kelvin-Helmholtz instability
+// -----------------------------------------------------------------------------
+//  Two horizontal layers in shear. A small vertical-velocity perturbation
+//  localized at the two interfaces grows into the characteristic rolled-up
+//  billows.
+//
+//  The layer states, the shear and the seeded perturbation follow the widely
+//  reproduced set-up of
+//
+//      V. Springel, "E pur si muove: Galilean-invariant cosmological
+//      hydrodynamical simulations on a moving mesh", Mon. Not. R. Astron. Soc.
+//      401 (2) (2010) 791-851,
+//      https://doi.org/10.1111/j.1365-2966.2009.15715.x
+//
+//  with two departures from it, which matter as soon as a figure here is
+//  compared to a published one:
+//
+//    - gamma is 1.4 here against 5/3 there, so that the case shares the
+//      equation of state of every other case in this repository;
+//    - the boundaries are outflow rather than periodic, periodicity not being
+//      wired through the mesh configuration yet.
+//
+//  The interface is a discontinuity, which makes the growth of the billows
+//  depend on the resolution rather than converge to one answer:
+//
+//      C.P. McNally, W. Lyra, J.-C. Passy, "A well-posed Kelvin-Helmholtz
+//      instability test and comparison", Astrophys. J. Suppl. Ser. 201 (2)
+//      (2012) 18, https://doi.org/10.1088/0067-0049/201/2/18
+//
+//  So the case is qualitative: it shows the scheme developing the instability
+//  and the adaptation following it. No number should be read off it. What earns
+//  it a place in the non-regression suite is being a moving, structured,
+//  adapted solution.
+// =============================================================================
+
 namespace test_case::kelvin_helmholtz
 {
     using field_t = config<2>::field_t;

euler/init/riemann_2d.hpp

@@ -18,15 +18,33 @@
 //  The nineteen admissible combinations are classified in
 //
 //      P.D. Lax, X.-D. Liu, "Solution of two-dimensional Riemann problems of gas
-//      dynamics by positive schemes", SIAM J. Sci. Comput. 19 (1998) 319-340.
+//      dynamics by positive schemes", SIAM J. Sci. Comput. 19 (2) (1998)
+//      319-340, https://doi.org/10.1137/S1064827595291819
+//
+//  and are tabulated, with this numbering, in
+//
+//      A. Kurganov, E. Tadmor, "Solution of two-dimensional Riemann problems
+//      for gas dynamics without Riemann problem solvers", Numer. Methods
+//      Partial Differ. Equ. 18 (5) (2002) 584-608,
+//      https://doi.org/10.1002/num.10025
+//
+//  The states registered below are those tables, verbatim; the article this
+//  repository reproduces (eq. 4) takes configuration 3 as its reference case
+//  for performance.
 //
 //  Quadrants are numbered counter-clockwise from the upper right:
 //      q[0] : x >= x0, y >= y0        q[1] : x <  x0, y >= y0
 //      q[2] : x <  x0, y <  y0        q[3] : x >= x0, y <  y0
 //
-//  TODO (lot 2) the interface should sit at x0 = y0 = 0.8, and configuration 3
-//  has p = 0.029 in quadrant 3, not 0.29. Both are corrected in their own commit
-//  so that this one stays behaviour preserving.
+//  The interfaces sit at x0 = y0 = 0.8 on [0,1]^2, as in Lax & Liu and in the
+//  article this repository reproduces: with t_f = 0.8 the waves then fill
+//  the domain without reaching its boundary, where the outflow condition
+//  would pollute them.
+//
+//  All three configurations below are invariant under the reflection about the
+//  diagonal (x,y,u,v) -> (y,x,v,u), so their solution must be too.
+//  tests/test_validation.py asserts that, and a single mistyped digit in one
+//  quadrant is enough to break it.
 // =============================================================================
 
 namespace test_case::riemann_2d
@@ -78,32 +96,32 @@ namespace test_case::riemann_2d
 
     // --- configuration 3 -----------------------------------------------------
     inline const Config config_3{
-        0.5,
-        0.5,
+        0.8,
+        0.8,
         {PrimState<2>{1.5, 1.5, {0., 0.}},
           PrimState<2>{0.5323, 0.3, {1.206, 0.}},
-          PrimState<2>{0.138, 0.29, {1.206, 1.206}},
+          PrimState<2>{0.138, 0.029, {1.206, 1.206}},
           PrimState<2>{0.5323, 0.3, {0., 1.206}}}
     };
 
     // --- configuration 4 -----------------------------------------------------
     inline const Config config_4{
-        0.5,
-        0.5,
+        0.8,
+        0.8,
         {PrimState<2>{1.1, 1.1, {0., 0.}},
           PrimState<2>{0.5065, 0.35, {0.8939, 0.}},
-          PrimState<2>{1.1, 1.1, {0.8939, 0.89396}},
-          PrimState<2>{0.5065, 0.35, {0., 0.89396}}}
+          PrimState<2>{1.1, 1.1, {0.8939, 0.8939}},
+          PrimState<2>{0.5065, 0.35, {0., 0.8939}}}
     };
 
     // --- configuration 12 ----------------------------------------------------
     inline const Config config_12{
-        0.5,
-        0.5,
-        {PrimState<2>{0.5197, 0.4, {0., 0.}},
-          PrimState<2>{1., 1., {-0.6259, 0.}},
-          PrimState<2>{0.8, 1., {-0.6259, -0.6259}},
-          PrimState<2>{1., 1., {0., -0.6259}}}
+        0.8,
+        0.8,
+        {PrimState<2>{0.5313, 0.4, {0., 0.}},
+          PrimState<2>{1., 1., {0.7276, 0.}},
+          PrimState<2>{0.8, 1., {0., 0.}},
+          PrimState<2>{1., 1., {0., 0.7276}}}
     };
 
     // One definition per configuration; `--riemann-config` will replace these

euler/init/sedov_blast.hpp

@@ -23,7 +23,32 @@
 //  up immediately as a broken symmetry.
 //
 //  The definition is dimension agnostic — only the volume over which the energy
-//  is spread changes — so the same case serves euler_2d and euler_3d.
+//  is spread changes — so the same case serves euler_1d, euler_2d and euler_3d.
+//
+//  Reference solution
+//  ------------------
+//      L.I. Sedov, "Similarity and Dimensional Methods in Mechanics", Academic
+//      Press, New York, 1959, chapter IV.
+//      J.R. Kamm, "Evaluation of the Sedov-von Neumann-Taylor Blast Wave
+//      Solution", report LA-UR-00-6055, Los Alamos National Laboratory, 2000.
+//      J.R. Kamm, F.X. Timmes, "On Efficient Generation of Numerically Robust
+//      Sedov Solutions", report LA-UR-07-2849, Los Alamos National Laboratory,
+//      2007.
+//
+//  The similarity solution puts the shock at
+//
+//      r_s(t) = ( E t^2 / (alpha rho_0) )^(1/(nu+2)),
+//
+//  nu = 1, 2, 3 being the geometry and alpha the dimensionless energy of the
+//  solution, a number that depends on nu and gamma alone. Integrating the
+//  similarity equations for gamma = 1.4 gives
+//
+//      alpha = 1.077485 (planar, energy counted on both sides of the plane)
+//              0.984074 (cylindrical)
+//              0.851072 (spherical)
+//
+//  python/sedov_exact.py computes these constants and the radius they give,
+//  and the validation suite measures the computed shock against it.
 // =============================================================================
 
 namespace test_case::sedov_blast
@@ -32,26 +57,28 @@ namespace test_case::sedov_blast
     inline constexpr double p_ambient   = 1e-5; // ambient pressure (very small)
     inline constexpr double r_blast     = 0.1;  // blast radius
 
-    // Blast energy. These are the values that put the shock at r = 1 at t = 1 for
-    // gamma = 1.4 in the Kamm & Timmes verification suite (planar 0.0673185,
-    // spherical 0.851072; the same spherical value is used by clawpack and by
-    // lanl/HARD).
+    // Blast energy, one per geometry: the standard cases of the Kamm report,
+    // which fix E so that the shock sits at r = 0.5, 0.75 and 1 at t = 1 for
+    // gamma = 1.4 and rho_0 = 1. The spherical value is the one clawpack and
+    // lanl/HARD use as well, and equals alpha, spherical being the geometry
+    // whose shock is normalised to r = 1.
     //
-    // NOTE the two-dimensional value below is NOT that suite's cylindrical one,
-    // which is 0.311357. It predates this file and is left as it was, since
-    // changing it changes results; worth settling separately.
+    // The planar value is twice Kamm's 0.0673185 because his planar blast
+    // releases its energy on one side of the plane, while the deposit below
+    // straddles the origin and feeds two shocks. Doubling it gives each shock
+    // his energy, and puts both where his solution puts its own, at 0.5.
     template <std::size_t dim>
     constexpr double blast_energy()
     {
         static_assert(dim >= 1 && dim <= 3, "no blast energy tabulated for this dimension");
 
         if constexpr (dim == 1)
         {
-            return 0.0673185;
+            return 0.134637; // = 2 x 0.0673185, see above
         }
         else if constexpr (dim == 2)
         {
-            return 0.244816;
+            return 0.311357;
         }
         else
         {

euler/init/sod.hpp

@@ -19,6 +19,22 @@
 //  that the solution is not aligned with the mesh. Any directional bias in the
 //  scheme or in the adaptation shows up as a distortion of what should stay a
 //  planar wave.
+//
+//  The states are those of
+//
+//      G.A. Sod, "A survey of several finite difference methods for systems of
+//      nonlinear hyperbolic conservation laws", J. Comput. Phys. 27 (1) (1978)
+//      1-31, https://doi.org/10.1016/0021-9991(78)90023-2
+//
+//  reprinted as test 1 of Table 4.1 of
+//
+//      E.F. Toro, "Riemann Solvers and Numerical Methods for Fluid Dynamics.
+//      A Practical Introduction", 3rd ed., Springer, 2009,
+//      https://doi.org/10.1007/b79761
+//
+//  whose exact solution python/exact_riemann.py computes: p* = 0.30313,
+//  u* = 0.92745. The 45 degree rotation is ours, not Sod's; the article this
+//  repository reproduces uses the axis-aligned tube.
 // =============================================================================
 
 namespace test_case::sod