diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index e2cf145..498fd4e 100644
--- a/.github/workflows/ci.yml
+++ b/.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
 
diff --git a/euler/init/closed_box.hpp b/euler/init/closed_box.hpp
index 5cc5d69..7cd6437 100644
--- a/euler/init/closed_box.hpp
+++ b/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
diff --git a/euler/init/double_mach.hpp b/euler/init/double_mach.hpp
index c0b5625..0428131 100644
--- a/euler/init/double_mach.hpp
+++ b/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
diff --git a/euler/init/double_rarefaction.hpp b/euler/init/double_rarefaction.hpp
index 85a89e7..592e7d2 100644
--- a/euler/init/double_rarefaction.hpp
+++ b/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
diff --git a/euler/init/free_stream.hpp b/euler/init/free_stream.hpp
index 4977830..92b998c 100644
--- a/euler/init/free_stream.hpp
+++ b/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
diff --git a/euler/init/isentropic_vortex.hpp b/euler/init/isentropic_vortex.hpp
index e0d7eae..9decc8e 100644
--- a/euler/init/isentropic_vortex.hpp
+++ b/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.
diff --git a/euler/init/kelvin_helmholtz.hpp b/euler/init/kelvin_helmholtz.hpp
index af9433c..255b141 100644
--- a/euler/init/kelvin_helmholtz.hpp
+++ b/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;
diff --git a/euler/init/riemann_2d.hpp b/euler/init/riemann_2d.hpp
index 35f0948..9317f9c 100644
--- a/euler/init/riemann_2d.hpp
+++ b/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
diff --git a/euler/init/sedov_blast.hpp b/euler/init/sedov_blast.hpp
index 4486545..7870f23 100644
--- a/euler/init/sedov_blast.hpp
+++ b/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,14 +57,16 @@ 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()
     {
@@ -47,11 +74,11 @@ namespace test_case::sedov_blast
 
         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
         {
diff --git a/euler/init/sod.hpp b/euler/init/sod.hpp
index ce23ffc..3082300 100644
--- a/euler/init/sod.hpp
+++ b/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
diff --git a/python/exact_riemann.py b/python/exact_riemann.py
new file mode 100644
index 0000000..3eda40a
--- /dev/null
+++ b/python/exact_riemann.py
@@ -0,0 +1,187 @@
+# Copyright 2025 the samurai team
+# SPDX-License-Identifier:  BSD-3-Clause
+"""
+Exact solution of the one-dimensional Riemann problem for an ideal gas.
+
+The shock tube cases are checked against this. A convergence order says how
+fast an error decays; this says where the shock, the contact and the
+rarefaction are.
+
+The algorithm is the classical one 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
+
+chapter 4: the star pressure solves the pressure-function equation
+
+    f_L(p) + f_R(p) + (u_R - u_L) = 0,
+
+f_K being the shock (Rankine-Hugoniot) branch for p > p_K and the rarefaction
+(isentropic) branch for p <= p_K. It is solved by Newton iterations from the
+two-rarefaction guess, then the solution is sampled at a given x/t.
+
+Run the module to reproduce the star states of Toro's Table 4.3:
+
+    python exact_riemann.py
+"""
+
+import numpy as np
+
+__all__ = ["star_state", "sample", "solution"]
+
+
+def _f(p, rho_k, p_k, gamma):
+    """Toro eq. (4.6)-(4.7): pressure function of one side, and its derivative."""
+    a_k = np.sqrt(gamma * p_k / rho_k)
+
+    if p > p_k:  # shock
+        A = 2.0 / ((gamma + 1.0) * rho_k)
+        B = (gamma - 1.0) / (gamma + 1.0) * p_k
+        root = np.sqrt(A / (B + p))
+        return (p - p_k) * root, root * (1.0 - 0.5 * (p - p_k) / (B + p))
+
+    # rarefaction
+    ratio = p / p_k
+    power = (gamma - 1.0) / (2.0 * gamma)
+    return (
+        2.0 * a_k / (gamma - 1.0) * (ratio**power - 1.0),
+        1.0 / (rho_k * a_k) * ratio ** (-(gamma + 1.0) / (2.0 * gamma)),
+    )
+
+
+def star_state(left, right, gamma=1.4, tol=1e-14, maxiter=100):
+    """Pressure and velocity of the star region between the two states.
+
+    `left` and `right` are (rho, u, p) triples. Raises if the two states have
+    already produced a vacuum, which is the case the iteration cannot solve
+    (Toro sec. 4.3.3, pressure positivity condition).
+    """
+    rho_l, u_l, p_l = left
+    rho_r, u_r, p_r = right
+    a_l = np.sqrt(gamma * p_l / rho_l)
+    a_r = np.sqrt(gamma * p_r / rho_r)
+
+    du = u_r - u_l
+    if 2.0 / (gamma - 1.0) * (a_l + a_r) <= du:
+        raise ValueError("vacuum generated by the initial data; no star state")
+
+    # two-rarefaction guess (Toro eq. 4.46), always positive
+    power = (gamma - 1.0) / (2.0 * gamma)
+    p = (
+        (a_l + a_r - 0.5 * (gamma - 1.0) * du)
+        / (a_l / p_l**power + a_r / p_r**power)
+    ) ** (1.0 / power)
+
+    for _ in range(maxiter):
+        f_l, df_l = _f(p, rho_l, p_l, gamma)
+        f_r, df_r = _f(p, rho_r, p_r, gamma)
+        step = (f_l + f_r + du) / (df_l + df_r)
+        p_new = max(p - step, tol)
+        if abs(p_new - p) / (0.5 * (p_new + p)) < tol:
+            p = p_new
+            break
+        p = p_new
+    else:
+        raise RuntimeError("the star pressure iteration did not converge")
+
+    f_l, _ = _f(p, rho_l, p_l, gamma)
+    f_r, _ = _f(p, rho_r, p_r, gamma)
+    return p, 0.5 * (u_l + u_r + f_r - f_l)
+
+
+def _sample_side(s, state, p_star, u_star, gamma, mirror):
+    """Sample the side of the contact the ray s falls on.
+
+    Written once for both sides. `mirror` is +1 on the left and -1 on the right:
+    negating x and u turns the right state into a left one, so the left-running
+    wave formulas of Toro sec. 4.5 serve both, which is all his two symmetric
+    sampling branches say.
+    """
+    rho_k, u_k, p_k = state
+    a_k = np.sqrt(gamma * p_k / rho_k)
+
+    s, u_k, u_star = mirror * s, mirror * u_k, mirror * u_star
+    ratio = p_star / p_k
+
+    if p_star > p_k:  # shock, Toro eq. (4.52)-(4.53)
+        beta = (gamma - 1.0) / (gamma + 1.0)
+        speed = u_k - a_k * np.sqrt((gamma + 1.0) / (2.0 * gamma) * ratio + (gamma - 1.0) / (2.0 * gamma))
+        if s <= speed:
+            return rho_k, mirror * u_k, p_k
+        return rho_k * (ratio + beta) / (beta * ratio + 1.0), mirror * u_star, p_star
+
+    # rarefaction, Toro eq. (4.55)-(4.56)
+    a_star = a_k * ratio ** ((gamma - 1.0) / (2.0 * gamma))
+    head, tail = u_k - a_k, u_star - a_star
+    if s <= head:
+        return rho_k, mirror * u_k, p_k
+    if s >= tail:
+        return rho_k * ratio ** (1.0 / gamma), mirror * u_star, p_star
+
+    # inside the fan
+    u = 2.0 / (gamma + 1.0) * (a_k + 0.5 * (gamma - 1.0) * u_k + s)
+    a = 2.0 / (gamma + 1.0) * (a_k + 0.5 * (gamma - 1.0) * (u_k - s))
+    return (
+        rho_k * (a / a_k) ** (2.0 / (gamma - 1.0)),
+        mirror * u,
+        p_k * (a / a_k) ** (2.0 * gamma / (gamma - 1.0)),
+    )
+
+
+def sample(s, left, right, p_star, u_star, gamma=1.4):
+    """State (rho, u, p) on the ray x/t = s."""
+    if s <= u_star:
+        return _sample_side(s, left, p_star, u_star, gamma, +1)
+    return _sample_side(s, right, p_star, u_star, gamma, -1)
+
+
+def solution(x, t, left, right, x0=0.5, gamma=1.4):
+    """The exact solution at time t on the points x, as (rho, u, p) arrays."""
+    p_star, u_star = star_state(left, right, gamma)
+    s = (np.asarray(x, dtype=float) - x0) / t
+    out = np.array([sample(si, left, right, p_star, u_star, gamma) for si in s])
+    return out[:, 0], out[:, 1], out[:, 2]
+
+
+# ---------------------------------------------------------------------------
+# Self test
+# ---------------------------------------------------------------------------
+# The five shock tubes of Toro's Table 4.1, with the star states of his Table
+# 4.3. Reproducing them is what makes this module usable as a reference. The
+# values are kept as the book prints them, and the comparison only asks for the
+# digits it actually printed: one unit of the last one, or six significant
+# digits, whichever is looser. Test 5 needs the second bound, its initial data
+# being the star state of tests 3 and 4 as the book rounded it before reusing
+# it.
+TORO_TESTS = {
+    1: ((1.0, 0.0, 1.0), (0.125, 0.0, 0.1), "0.30313", "0.92745"),
+    2: ((1.0, -2.0, 0.4), (1.0, 2.0, 0.4), "0.00189", "0.00000"),
+    3: ((1.0, 0.0, 1000.0), (1.0, 0.0, 0.01), "460.894", "19.5975"),
+    4: ((1.0, 0.0, 0.01), (1.0, 0.0, 100.0), "46.0950", "-6.19633"),
+    5: ((5.99924, 19.5975, 460.894), (5.99242, -6.19633, 46.0950), "1691.64", "8.68975"),
+}
+
+
+def _printed_tolerance(text):
+    """How far from a published value a value may sit and still match it."""
+    decimals = len(text.split(".")[1]) if "." in text else 0
+    return max(10.0**-decimals, 5e-6 * abs(float(text)))
+
+
+def _self_test():
+    print(f"{'test':>4}  {'p* (Toro)':>10} {'p* (here)':>14}  {'u* (Toro)':>10} {'u* (here)':>14}")
+    ok = True
+    for number, (left, right, p_ref, u_ref) in TORO_TESTS.items():
+        p_star, u_star = star_state(left, right)
+        print(f"{number:>4}  {p_ref:>10} {p_star:>14.6g}  {u_ref:>10} {u_star:>14.6g}")
+        ok &= abs(p_star - float(p_ref)) <= _printed_tolerance(p_ref)
+        ok &= abs(u_star - float(u_ref)) <= _printed_tolerance(u_ref)
+    print("\nagrees with Table 4.3 to the printed digits" if ok else "\nDOES NOT agree with Table 4.3")
+    return ok
+
+
+if __name__ == "__main__":
+    import sys
+
+    sys.exit(0 if _self_test() else 1)
diff --git a/python/sedov_exact.py b/python/sedov_exact.py
new file mode 100644
index 0000000..fda2d93
--- /dev/null
+++ b/python/sedov_exact.py
@@ -0,0 +1,154 @@
+# Copyright 2025 the samurai team
+# SPDX-License-Identifier:  BSD-3-Clause
+"""
+Exact self-similar solution of the Sedov blast wave, in one, two and three
+dimensions.
+
+This is what the Sedov case is checked against. A point-like release of energy E
+in an ambient gas at rest drives a shock whose radius obeys
+
+    r_s(t) = ( E t^2 / (alpha rho_0) )^(1/(nu+2)),
+
+nu = 1, 2, 3 being planar, cylindrical and spherical geometry. The constant
+alpha depends on nu and gamma alone, and is computed here rather than quoted:
+the similarity equations are integrated inward from the shock, where the strong
+shock jump conditions give the starting values, and alpha follows from the
+energy integral over the profile.
+
+    L.I. Sedov, "Similarity and Dimensional Methods in Mechanics", Academic
+    Press, New York, 1959, chapter IV.
+    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 planar constant is computed with the energy counted on both sides of the
+plane, which is how euler/init/sedov_blast.hpp deposits it. Kamm's planar
+standard case counts one side, so his planar energy corresponds to half of
+this one; the cylindrical and spherical cases have no such convention.
+
+Run the module to print the constants and check the spherical one against the
+value published for gamma = 1.4:
+
+    python sedov_exact.py
+"""
+
+import numpy as np
+
+__all__ = ["alpha", "shock_radius", "profile"]
+
+# The similarity integration is deterministic and takes about half a second, so
+# the constants are memoised rather than tabulated: no number in this file is
+# copied from a paper.
+_ALPHA_CACHE = {}
+
+
+def _derivatives(lam, y, nu, gamma):
+    """Right-hand side of the similarity equations, in the variables
+
+        rho = rho_0 G(lam),  u = rdot_s V(lam),  p = rho_0 rdot_s^2 P(lam),
+
+    lam = r / r_s. Continuity, momentum and the entropy condition, with the
+    time dependence eliminated through r_s ~ t^(2/(nu+2)), leave three coupled
+    ordinary equations in lam.
+    """
+    G, V, P = y
+    D = V - lam
+    Vp = (nu - gamma * (nu - 1) * V / lam - nu * V * D * G / (2.0 * P)) / (gamma - D * D * G / P)
+    Gp = -(G * Vp + (nu - 1) * G * V / lam) / D
+    Pp = G * (nu * V / 2.0 - D * Vp)
+    return np.array([Gp, Vp, Pp])
+
+
+def _integrate(nu, gamma, steps, lam_min):
+    """Integrate from the shock at lam = 1 down to lam_min, by RK4.
+
+    Returns lam increasing, with G, V, P alongside.
+    """
+    y = np.array(
+        [
+            (gamma + 1.0) / (gamma - 1.0),  # strong shock jump conditions
+            2.0 / (gamma + 1.0),
+            2.0 / (gamma + 1.0),
+        ]
+    )
+
+    h = -(1.0 - lam_min) / steps
+    lam = 1.0
+    lams, ys = [lam], [y.copy()]
+
+    for _ in range(steps):
+        k1 = _derivatives(lam, y, nu, gamma)
+        k2 = _derivatives(lam + 0.5 * h, y + 0.5 * h * k1, nu, gamma)
+        k3 = _derivatives(lam + 0.5 * h, y + 0.5 * h * k2, nu, gamma)
+        k4 = _derivatives(lam + h, y + h * k3, nu, gamma)
+        y = y + h / 6.0 * (k1 + 2.0 * k2 + 2.0 * k3 + k4)
+        lam += h
+        if not np.all(np.isfinite(y)) or y[2] <= 0.0:
+            break
+        lams.append(lam)
+        ys.append(y.copy())
+
+    return np.array(lams)[::-1], np.array(ys)[::-1]
+
+
+def alpha(dim, gamma=1.4, steps=20000, lam_min=1e-6):
+    """Dimensionless energy of the similarity solution, for this geometry.
+
+    The energy of the solution, in the similarity variables, is
+
+        E = sigma_nu rho_0 r_s^nu rdot_s^2 Int, Int = int_0^1 (G V^2 / 2 +
+            P / (gamma - 1)) lam^(nu-1) dlam,
+
+    and alpha = E t^2 / (rho_0 r_s^(nu+2)) follows from rdot_s = 2 r_s /
+    ((nu+2) t). sigma_nu is 2, 2 pi, 4 pi: the planar solution is counted on
+    both sides of the plane.
+    """
+    key = (dim, gamma, steps, lam_min)
+    if key not in _ALPHA_CACHE:
+        lam, y = _integrate(dim, gamma, steps, lam_min)
+        G, V, P = y[:, 0], y[:, 1], y[:, 2]
+        integral = np.trapezoid((0.5 * G * V * V + P / (gamma - 1.0)) * lam ** (dim - 1), lam)
+        sigma = {1: 2.0, 2: 2.0 * np.pi, 3: 4.0 * np.pi}[dim]
+        _ALPHA_CACHE[key] = 4.0 * sigma * integral / (dim + 2.0) ** 2
+    return _ALPHA_CACHE[key]
+
+
+def shock_radius(energy, t, dim, gamma=1.4, rho_ambient=1.0):
+    """Where the shock stands at time t."""
+    return (energy * t * t / (alpha(dim, gamma) * rho_ambient)) ** (1.0 / (dim + 2.0))
+
+
+def profile(energy, t, dim, gamma=1.4, rho_ambient=1.0, p_ambient=0.0):
+    """The exact solution at time t, as (r, rho, u, p) sampled on the profile.
+
+    The points are those the integration produced, which cluster where the
+    solution varies; the last one is the shock. Outside it the gas is still
+    ambient, which the caller usually knows already.
+    """
+    lam, y = _integrate(dim, gamma, 20000, 1e-6)
+    r_s = shock_radius(energy, t, dim, gamma, rho_ambient)
+    speed = 2.0 * r_s / ((dim + 2.0) * t)
+    return (
+        lam * r_s,
+        rho_ambient * y[:, 0],
+        speed * y[:, 1],
+        rho_ambient * speed * speed * y[:, 2] + p_ambient,
+    )
+
+
+if __name__ == "__main__":
+    import sys
+
+    # Published for gamma = 1.4: the spherical blast reaches r = 1 at t = 1 for
+    # E = 0.851072, so alpha is that same number. It is the one constant of the
+    # three that several codes quote, hence the one worth checking against.
+    SPHERICAL = 0.851072
+
+    print(f"gamma = 1.4\n{'geometry':>12} {'alpha':>10}   E for r_s = 1 at t = 1")
+    for dim, name in ((1, "planar"), (2, "cylindrical"), (3, "spherical")):
+        a = alpha(dim)
+        print(f"{name:>12} {a:>10.6f}   {a:.6f}")
+
+    deviation = abs(alpha(3) - SPHERICAL)
+    print(f"\nspherical against the published {SPHERICAL}: {deviation:.1e}")
+    sys.exit(0 if deviation < 5e-7 else 1)
diff --git a/tests/reference/euler_1d_sedov_blast_hll.npz b/tests/reference/euler_1d_sedov_blast_hll.npz
index a520356..2e4cbf4 100644
Binary files a/tests/reference/euler_1d_sedov_blast_hll.npz and b/tests/reference/euler_1d_sedov_blast_hll.npz differ
diff --git a/tests/reference/euler_1d_sedov_blast_hllc.npz b/tests/reference/euler_1d_sedov_blast_hllc.npz
index 65c6c73..0490b02 100644
Binary files a/tests/reference/euler_1d_sedov_blast_hllc.npz and b/tests/reference/euler_1d_sedov_blast_hllc.npz differ
diff --git a/tests/reference/euler_1d_sedov_blast_rusanov.npz b/tests/reference/euler_1d_sedov_blast_rusanov.npz
index df75d03..bd694cf 100644
Binary files a/tests/reference/euler_1d_sedov_blast_rusanov.npz and b/tests/reference/euler_1d_sedov_blast_rusanov.npz differ
diff --git a/tests/reference/euler_2d_riemann2d_config12_hll.npz b/tests/reference/euler_2d_riemann2d_config12_hll.npz
index 308baff..f512902 100644
Binary files a/tests/reference/euler_2d_riemann2d_config12_hll.npz and b/tests/reference/euler_2d_riemann2d_config12_hll.npz differ
diff --git a/tests/reference/euler_2d_riemann2d_config12_hllc.npz b/tests/reference/euler_2d_riemann2d_config12_hllc.npz
index d30f604..cff9930 100644
Binary files a/tests/reference/euler_2d_riemann2d_config12_hllc.npz and b/tests/reference/euler_2d_riemann2d_config12_hllc.npz differ
diff --git a/tests/reference/euler_2d_riemann2d_config12_rusanov.npz b/tests/reference/euler_2d_riemann2d_config12_rusanov.npz
index ac6330a..a988481 100644
Binary files a/tests/reference/euler_2d_riemann2d_config12_rusanov.npz and b/tests/reference/euler_2d_riemann2d_config12_rusanov.npz differ
diff --git a/tests/reference/euler_2d_riemann2d_config3_hll.npz b/tests/reference/euler_2d_riemann2d_config3_hll.npz
index 5ba167f..05dc737 100644
Binary files a/tests/reference/euler_2d_riemann2d_config3_hll.npz and b/tests/reference/euler_2d_riemann2d_config3_hll.npz differ
diff --git a/tests/reference/euler_2d_riemann2d_config3_hllc.npz b/tests/reference/euler_2d_riemann2d_config3_hllc.npz
index 06c6023..4f38dba 100644
Binary files a/tests/reference/euler_2d_riemann2d_config3_hllc.npz and b/tests/reference/euler_2d_riemann2d_config3_hllc.npz differ
diff --git a/tests/reference/euler_2d_riemann2d_config3_rusanov.npz b/tests/reference/euler_2d_riemann2d_config3_rusanov.npz
index 01f6115..a3f6bbc 100644
Binary files a/tests/reference/euler_2d_riemann2d_config3_rusanov.npz and b/tests/reference/euler_2d_riemann2d_config3_rusanov.npz differ
diff --git a/tests/reference/euler_2d_riemann2d_config4_hll.npz b/tests/reference/euler_2d_riemann2d_config4_hll.npz
index 9a5b57c..bf56302 100644
Binary files a/tests/reference/euler_2d_riemann2d_config4_hll.npz and b/tests/reference/euler_2d_riemann2d_config4_hll.npz differ
diff --git a/tests/reference/euler_2d_riemann2d_config4_hllc.npz b/tests/reference/euler_2d_riemann2d_config4_hllc.npz
index 18c1fc1..495e86f 100644
Binary files a/tests/reference/euler_2d_riemann2d_config4_hllc.npz and b/tests/reference/euler_2d_riemann2d_config4_hllc.npz differ
diff --git a/tests/reference/euler_2d_riemann2d_config4_rusanov.npz b/tests/reference/euler_2d_riemann2d_config4_rusanov.npz
index 1f1de4b..4d5176f 100644
Binary files a/tests/reference/euler_2d_riemann2d_config4_rusanov.npz and b/tests/reference/euler_2d_riemann2d_config4_rusanov.npz differ
diff --git a/tests/reference/euler_2d_sedov_blast_hll.npz b/tests/reference/euler_2d_sedov_blast_hll.npz
index 6e3f49c..21c8840 100644
Binary files a/tests/reference/euler_2d_sedov_blast_hll.npz and b/tests/reference/euler_2d_sedov_blast_hll.npz differ
diff --git a/tests/reference/euler_2d_sedov_blast_hllc.npz b/tests/reference/euler_2d_sedov_blast_hllc.npz
index ad86e9a..f5ae8b9 100644
Binary files a/tests/reference/euler_2d_sedov_blast_hllc.npz and b/tests/reference/euler_2d_sedov_blast_hllc.npz differ
diff --git a/tests/reference/euler_2d_sedov_blast_rusanov.npz b/tests/reference/euler_2d_sedov_blast_rusanov.npz
index 89686d7..89f1b1c 100644
Binary files a/tests/reference/euler_2d_sedov_blast_rusanov.npz and b/tests/reference/euler_2d_sedov_blast_rusanov.npz differ
diff --git a/tests/test_validation.py b/tests/test_validation.py
index debc3a7..dfa5635 100644
--- a/tests/test_validation.py
+++ b/tests/test_validation.py
@@ -17,7 +17,7 @@
 import numpy as np
 import pytest
 
-from util import ROOT, level_count, read, run_case
+from util import ROOT, level_count, read, run_case, sedov_blast_energy
 
 sys.path.insert(0, str(ROOT / "python"))
 from error_analysis import errors, exact_vortex  # noqa: E402
@@ -91,3 +91,118 @@ def test_adaptation_costs_no_accuracy(tmp_path):
     assert adapted_l1 < 1.5 * uniform_l1, (
         f"adapted L1 {adapted_l1:.3e} against uniform {uniform_l1:.3e}"
     )
+
+
+# ---------------------------------------------------------------------------
+# The published test cases, against the solution their paper gives
+# ---------------------------------------------------------------------------
+# These check that the cases reproduce the papers they were taken from.
+# test_regression.py checks something else, that they still compute what they
+# computed yesterday. The thresholds are what a first-order scheme
+# reaches at these resolutions, with room to spare; they should be tightened
+# when the MUSCL-Hancock update lands.
+from exact_riemann import solution as exact_riemann_solution  # noqa: E402
+from exact_riemann import star_state  # noqa: E402
+from sedov_exact import shock_radius  # noqa: E402
+
+
+def test_double_rarefaction_matches_toro_test_2(tmp_path):
+    """The 123 problem against the exact solution of Toro's Table 4.1, test 2."""
+    level, tf = 12, 0.15
+    left, right = (1.0, -2.0, 0.4), (1.0, 2.0, 0.4)
+
+    out, stem = run_case("euler_1d", tmp_path, "double_rarefaction",
+                         min_level=level, max_level=level, Tf=tf)
+    centers, volume, fields = read(out / stem)
+    x = centers[:, 0]
+    rho, u, p = exact_riemann_solution(x, tf, left, right, x0=0.5)
+
+    def l1(computed, exact):
+        return float(np.sum(np.abs(computed - exact) * volume) / np.sum(volume))
+
+    assert l1(fields["rho"], rho) < 1e-2
+    assert l1(fields["velocity"][:, 0], u) < 3e-2
+    assert l1(fields["pressure"], p) < 1e-2
+
+    # the near-vacuum is what the case is for: a first-order scheme sits above
+    # the exact star pressure, and must not sit far above it
+    p_star, _ = star_state(left, right)
+    assert p_star < fields["pressure"].min() < 3.0 * p_star
+
+
+def test_sod_matches_its_exact_solution(tmp_path):
+    """Sod's tube, rotated 45 degrees, against the exact solution.
+
+    Two things at once: the waves must be in the right place, and the solution
+    must stay one-dimensional along the diagonal. The transverse velocity is the
+    isotropy measure the rotation was introduced for.
+    """
+    level, tf = 8, 0.2
+    left, right = (1.0, 0.0, 1.0), (0.125, 0.0, 0.1)
+
+    out, stem = run_case("euler_2d", tmp_path, "sod",
+                         min_level=level, max_level=level, Tf=tf)
+    centers, volume, fields = read(out / stem)
+    x, y = centers[:, 0], centers[:, 1]
+
+    # coordinate across the interface, and the velocity split along and across it
+    band = np.abs(x - y) < 0.3  # away from the corners, where outflow is not exact
+    xi = (x + y - 1.0) / np.sqrt(2.0)
+    along = (fields["velocity"][:, 0] + fields["velocity"][:, 1]) / np.sqrt(2.0)
+    across = (fields["velocity"][:, 0] - fields["velocity"][:, 1]) / np.sqrt(2.0)
+
+    rho, u, p = exact_riemann_solution(xi[band], tf, left, right, x0=0.0)
+    weight = volume[band]
+
+    def l1(computed, exact):
+        return float(np.sum(np.abs(computed - exact) * weight) / np.sum(weight))
+
+    assert l1(fields["rho"][band], rho) < 2e-2
+    assert l1(along[band], u) < 2e-2
+    assert l1(fields["pressure"][band], p) < 2e-2
+    assert np.abs(across[band]).max() < 1e-2
+
+
+@pytest.mark.parametrize("binary,dim,level,tf", [("euler_2d", 2, 9, 0.6), ("euler_1d", 1, 12, 0.6)])
+def test_sedov_shock_sits_where_the_similarity_solution_puts_it(binary, dim, level, tf, tmp_path):
+    """The blast energy is only meaningful through the shock radius it produces.
+
+    Measured as the outermost radius at which the density is still above halfway
+    to its peak, which is where a smeared shock front has its middle. A wrong
+    blast energy shows up here and nowhere else: every other check of this case
+    is a symmetry or a positivity, and both survive any energy at all.
+    """
+    out, stem = run_case(binary, tmp_path, "sedov_blast",
+                         min_level=level, max_level=level, Tf=tf)
+    centers, _, fields = read(out / stem)
+
+    r = np.linalg.norm(centers[:, :dim], axis=1)
+    rho = fields["rho"]
+    front = r[rho > 0.5 * (rho.max() + 1.0)].max()
+
+    expected = shock_radius(sedov_blast_energy(dim), tf, dim)
+    assert abs(front / expected - 1.0) < 0.1, f"shock at {front:.4f}, similarity solution at {expected:.4f}"
+
+
+@pytest.mark.parametrize("case", ["riemann2d_config3", "riemann2d_config4", "riemann2d_config12"])
+def test_riemann_2d_keeps_the_symmetry_of_its_configuration(case, tmp_path):
+    """Lax & Liu configurations 3, 4 and 12 are symmetric about the diagonal.
+
+    Their initial data is invariant under (x,y,u,v) -> (y,x,v,u), so the solution
+    is too, and a uniform mesh carries that symmetry exactly. It is the cheapest
+    check that the states were copied correctly: a single mistyped digit in one
+    quadrant breaks it, while leaving a picture that still looks plausible.
+    """
+    level, tf = 7, 0.2
+    out, stem = run_case("euler_2d", tmp_path, case,
+                         min_level=level, max_level=level, Tf=tf)
+    centers, _, fields = read(out / stem)
+
+    n = int(round(np.sqrt(centers.shape[0])))
+    order = np.lexsort((centers[:, 0], centers[:, 1]))
+    rho = fields["rho"][order].reshape(n, n)
+    vx = fields["velocity"][order, 0].reshape(n, n)
+    vy = fields["velocity"][order, 1].reshape(n, n)
+
+    assert np.abs(rho - rho.T).max() < 1e-12
+    assert np.abs(vx - vy.T).max() < 1e-12
diff --git a/tests/util.py b/tests/util.py
index 18fc6ea..b879fea 100644
--- a/tests/util.py
+++ b/tests/util.py
@@ -10,6 +10,7 @@
 """
 
 import os
+import re
 import subprocess
 from pathlib import Path
 
@@ -120,6 +121,19 @@ def finest_cell_size(volume, dim):
     return volume.min() ** (1.0 / dim)
 
 
+def sedov_blast_energy(dim):
+    """The blast energy euler/init/sedov_blast.hpp compiles in, read from it.
+
+    Restating the three numbers here would make a second source of truth that
+    eventually stops agreeing with the first.
+    """
+    text = (ROOT / "euler" / "init" / "sedov_blast.hpp").read_text()
+    body = text.split("constexpr double blast_energy()", 1)[1]
+    values = re.findall(r"return\s+([0-9.eE+-]+);", body)[:3]
+    assert len(values) == 3, f"expected one blast energy per dimension, read {values}"
+    return float(values[dim - 1])
+
+
 def level_count(volume):
     """Number of distinct cell sizes: 1 on a uniform mesh, more once adapted.