GitHub

Toro Riemann Problem Suite

This example runs the five standard Euler Riemann problems from Toro (2009), Table 4.1. Each test exercises different aspects of the Riemann solver: rarefactions, contacts, shocks, near-vacuum, and strong discontinuities.

Reference

  • Toro, E.F. (2009). Riemann Solvers and Numerical Methods for Fluid Dynamics, 3rd ed., Springer. Chapter 4, Table 4.1.
using FiniteVolumeMethod
using StaticArrays
using CairoMakie

gamma = 1.4
eos = IdealGasEOS(gamma)
law = EulerEquations{1}(eos)
N = 400
400

Exact Riemann Solver (star-state values from Toro Table 4.1)

toro_tests = [
    (
        name = "Sod", rhoL = 1.0, uL = 0.0, pL = 1.0,
        rhoR = 0.125, uR = 0.0, pR = 0.1,
        x0 = 0.5, t = 0.2, p_star = 0.30313, u_star = 0.92745,
    ),
    (
        name = "123 (two rarefactions)", rhoL = 1.0, uL = -2.0, pL = 0.4,
        rhoR = 1.0, uR = 2.0, pR = 0.4,
        x0 = 0.5, t = 0.15, p_star = 0.00189, u_star = 0.0,
    ),
    (
        name = "Woodward-Colella blast", rhoL = 1.0, uL = 0.0, pL = 1000.0,
        rhoR = 1.0, uR = 0.0, pR = 0.01,
        x0 = 0.5, t = 0.012, p_star = 460.894, u_star = 19.5975,
    ),
    (
        name = "Two-shock collision", rhoL = 5.99924, uL = 19.5975, pL = 460.894,
        rhoR = 5.99242, uR = -6.19633, pR = 46.095,
        x0 = 0.4, t = 0.035, p_star = 1691.64, u_star = 8.68975,
    ),
    (
        name = "Stationary contact", rhoL = 1.0, uL = -19.5975, pL = 1000.0,
        rhoR = 1.0, uR = -19.5975, pR = 0.01,
        x0 = 0.8, t = 0.012, p_star = 460.894, u_star = 0.0,
    ),
]
5-element Vector{@NamedTuple{name::String, rhoL::Float64, uL::Float64, pL::Float64, rhoR::Float64, uR::Float64, pR::Float64, x0::Float64, t::Float64, p_star::Float64, u_star::Float64}}:
 (name = "Sod", rhoL = 1.0, uL = 0.0, pL = 1.0, rhoR = 0.125, uR = 0.0, pR = 0.1, x0 = 0.5, t = 0.2, p_star = 0.30313, u_star = 0.92745)
 (name = "123 (two rarefactions)", rhoL = 1.0, uL = -2.0, pL = 0.4, rhoR = 1.0, uR = 2.0, pR = 0.4, x0 = 0.5, t = 0.15, p_star = 0.00189, u_star = 0.0)
 (name = "Woodward-Colella blast", rhoL = 1.0, uL = 0.0, pL = 1000.0, rhoR = 1.0, uR = 0.0, pR = 0.01, x0 = 0.5, t = 0.012, p_star = 460.894, u_star = 19.5975)
 (name = "Two-shock collision", rhoL = 5.99924, uL = 19.5975, pL = 460.894, rhoR = 5.99242, uR = -6.19633, pR = 46.095, x0 = 0.4, t = 0.035, p_star = 1691.64, u_star = 8.68975)
 (name = "Stationary contact", rhoL = 1.0, uL = -19.5975, pL = 1000.0, rhoR = 1.0, uR = -19.5975, pR = 0.01, x0 = 0.8, t = 0.012, p_star = 460.894, u_star = 0.0)

Run All Five Tests

results = map(toro_tests) do tt
    wL = SVector(tt.rhoL, tt.uL, tt.pL)
    wR = SVector(tt.rhoR, tt.uR, tt.pR)
    ic(x) = x < tt.x0 ? wL : wR
    mesh = StructuredMesh1D(0.0, 1.0, N)
    prob = HyperbolicProblem(
        law, mesh, HLLCSolver(), CellCenteredMUSCL(MinmodLimiter()),
        DirichletHyperbolicBC(wL), DirichletHyperbolicBC(wR), ic;
        final_time = tt.t, cfl = 0.4,
    )
    x, U, t_end = solve_hyperbolic(prob)
    W = [conserved_to_primitive(law, U[i]) for i in eachindex(U)]
    (x = x, W = W, t = t_end, ref = tt)
end
5-element Vector{@NamedTuple{x::Vector{Float64}, W::Vector{StaticArraysCore.SVector{3, Float64}}, t::Float64, ref::@NamedTuple{name::String, rhoL::Float64, uL::Float64, pL::Float64, rhoR::Float64, uR::Float64, pR::Float64, x0::Float64, t::Float64, p_star::Float64, u_star::Float64}}}:
 (x = [0.00125, 0.00375, 0.00625, 0.00875, 0.01125, 0.01375, 0.01625, 0.01875, 0.02125, 0.02375  …  0.9762500000000001, 0.97875, 0.9812500000000001, 0.98375, 0.9862500000000001, 0.98875, 0.9912500000000001, 0.99375, 0.99625, 0.99875], W = [[1.0, 0.0, 1.0], [1.0, 0.0, 1.0], [1.0, 0.0, 1.0], [1.0, 0.0, 1.0], [1.0, 0.0, 1.0], [1.0, 0.0, 1.0], [1.0, 0.0, 1.0], [1.0, 0.0, 1.0], [1.0, 0.0, 1.0], [1.0, 0.0, 1.0]  …  [0.125, 0.0, 0.09999999999999998], [0.125, 0.0, 0.09999999999999998], [0.125, 0.0, 0.09999999999999998], [0.125, 0.0, 0.09999999999999998], [0.125, 0.0, 0.09999999999999998], [0.125, 0.0, 0.09999999999999998], [0.125, 0.0, 0.09999999999999998], [0.125, 0.0, 0.09999999999999998], [0.125, -7.467461762869597e-18, 0.09999999999999998], [0.125, -1.1251576643496256e-16, 0.09999999999999998]], t = 0.2, ref = (name = "Sod", rhoL = 1.0, uL = 0.0, pL = 1.0, rhoR = 0.125, uR = 0.0, pR = 0.1, x0 = 0.5, t = 0.2, p_star = 0.30313, u_star = 0.92745))
 (x = [0.00125, 0.00375, 0.00625, 0.00875, 0.01125, 0.01375, 0.01625, 0.01875, 0.02125, 0.02375  …  0.9762500000000001, 0.97875, 0.9812500000000001, 0.98375, 0.9862500000000001, 0.98875, 0.9912500000000001, 0.99375, 0.99625, 0.99875], W = [[0.9999999924470393, -1.9999999943478923, 0.3999999957703477], [0.9999999872205626, -1.999999990436756, 0.3999999928435213], [0.999999974961184, -1.9999999812626748, 0.3999999859782702], [0.9999999535229788, -1.9999999652197893, 0.3999999739728739], [0.9999999134681388, -1.9999999352454871, 0.3999999515421658], [0.9999998409735973, -1.9999998809955424, 0.399999910945226], [0.9999997102851319, -1.9999997831972498, 0.3999998377597041], [0.9999994775224721, -1.9999996090136, 0.39999970741266616], [0.9999990671115189, -1.9999993018901059, 0.3999994775826886], [0.9999983513408246, -1.9999987662561443, 0.39999907675157065]  …  [0.9999983513408257, 1.9999987662561427, 0.39999907675157115], [0.9999990671115142, 1.9999993018901012, 0.39999947758268967], [0.9999994775224644, 1.9999996090135994, 0.399999707412665], [0.9999997102851291, 1.9999997831972465, 0.39999983775970505], [0.9999998409735926, 1.9999998809955402, 0.399999910945228], [0.9999999134681334, 1.9999999352454871, 0.39999995154216583], [0.9999999535229729, 1.999999965219789, 0.39999997397287457], [0.9999999749611831, 1.9999999812626743, 0.3999999859782706], [0.9999999872205655, 1.9999999904367511, 0.39999999284352433], [0.9999999924470407, 1.999999994347888, 0.39999999577034995]], t = 0.15, ref = (name = "123 (two rarefactions)", rhoL = 1.0, uL = -2.0, pL = 0.4, rhoR = 1.0, uR = 2.0, pR = 0.4, x0 = 0.5, t = 0.15, p_star = 0.00189, u_star = 0.0))
 (x = [0.00125, 0.00375, 0.00625, 0.00875, 0.01125, 0.01375, 0.01625, 0.01875, 0.02125, 0.02375  …  0.9762500000000001, 0.97875, 0.9812500000000001, 0.98375, 0.9862500000000001, 0.98875, 0.9912500000000001, 0.99375, 0.99625, 0.99875], W = [[0.9999816646048082, 0.0006860496934482244, 999.974330540822], [0.9999727953352358, 0.0010179149659642246, 999.9619137088513], [0.9999558700285448, 0.0016512216421790368, 999.9382186304701], [0.9999321781279811, 0.0025377303057470334, 999.9050508004717], [0.9998953037627609, 0.003917538446349742, 999.8534286055984], [0.9998415555338904, 0.005928823592682054, 999.7781853817597], [0.9997627088899702, 0.008879462361115432, 999.66780945607], [0.9996498567386656, 0.013102995672045281, 999.5098363052141], [0.9994903150317591, 0.0190745574860844, 999.2865187585468], [0.9992688823940747, 0.02736393266237036, 998.9765944777629]  …  [1.0, 0.0, 0.01], [1.0, 0.0, 0.01], [1.0, 0.0, 0.01], [1.0, 0.0, 0.01], [1.0, 0.0, 0.01], [1.0, 0.0, 0.01], [1.0, 0.0, 0.01], [1.0, 0.0, 0.01], [1.0, 0.0, 0.01], [1.0, 0.0, 0.01]], t = 0.012, ref = (name = "Woodward-Colella blast", rhoL = 1.0, uL = 0.0, pL = 1000.0, rhoR = 1.0, uR = 0.0, pR = 0.01, x0 = 0.5, t = 0.012, p_star = 460.894, u_star = 19.5975))
 (x = [0.00125, 0.00375, 0.00625, 0.00875, 0.01125, 0.01375, 0.01625, 0.01875, 0.02125, 0.02375  …  0.9762500000000001, 0.97875, 0.9812500000000001, 0.98375, 0.9862500000000001, 0.98875, 0.9912500000000001, 0.99375, 0.99625, 0.99875], W = [[5.99924, 19.597499999999997, 460.8940000000003], [5.99924, 19.597499999999993, 460.8940000000004], [5.99924, 19.597499999999993, 460.8940000000004], [5.99924, 19.597499999999993, 460.8940000000004], [5.99924, 19.597499999999993, 460.8940000000004], [5.99924, 19.597499999999993, 460.8940000000004], [5.99924, 19.597499999999993, 460.8940000000004], [5.99924, 19.597499999999993, 460.8940000000004], [5.99924, 19.597499999999993, 460.8940000000004], [5.99924, 19.597499999999993, 460.8940000000004]  …  [5.992419999999999, -6.1963300000000014, 46.09499999999997], [5.992419999999999, -6.1963300000000014, 46.09499999999997], [5.992419999999999, -6.1963300000000014, 46.09499999999997], [5.992419999999999, -6.1963300000000014, 46.09499999999997], [5.992419999999999, -6.1963300000000014, 46.09499999999997], [5.992419999999999, -6.1963300000000014, 46.09499999999997], [5.992419999999999, -6.1963300000000014, 46.09499999999997], [5.992419999999999, -6.1963300000000014, 46.09499999999997], [5.992419999999999, -6.1963300000000014, 46.09499999999997], [5.992419999999999, -6.1963300000000014, 46.09499999999997]], t = 0.035, ref = (name = "Two-shock collision", rhoL = 5.99924, uL = 19.5975, pL = 460.894, rhoR = 5.99242, uR = -6.19633, pR = 46.095, x0 = 0.4, t = 0.035, p_star = 1691.64, u_star = 8.68975))
 (x = [0.00125, 0.00375, 0.00625, 0.00875, 0.01125, 0.01375, 0.01625, 0.01875, 0.02125, 0.02375  …  0.9762500000000001, 0.97875, 0.9812500000000001, 0.98375, 0.9862500000000001, 0.98875, 0.9912500000000001, 0.99375, 0.99625, 0.99875], W = [[0.9999999997665103, -19.597499991263707, 999.9999996731176], [0.9999999996248178, -19.597499985962088, 999.999999474749], [0.9999999993150854, -19.597499974372916, 999.9999990411216], [0.9999999988109625, -19.597499955510333, 999.9999983353481], [0.9999999979199174, -19.597499922170588, 999.9999970878882], [0.9999999964008268, -19.59749986533145, 999.9999949611617], [0.9999999938039246, -19.597499768164244, 999.9999913254974], [0.9999999894086053, -19.597499603706407, 999.9999851720486], [0.9999999820137644, -19.597499327016834, 999.9999748192723], [0.9999999696644861, -19.597498864949152, 999.9999575302833]  …  [1.0, -19.5975, 0.010000000000002271], [1.0, -19.5975, 0.010000000000002271], [1.0, -19.5975, 0.010000000000002271], [1.0, -19.5975, 0.010000000000002271], [1.0, -19.5975, 0.010000000000002271], [1.0, -19.5975, 0.010000000000002271], [1.0, -19.5975, 0.010000000000002271], [1.0, -19.5975, 0.010000000000002271], [1.0, -19.5975, 0.010000000000002271], [1.0, -19.5975, 0.010000000000002271]], t = 0.012, ref = (name = "Stationary contact", rhoL = 1.0, uL = -19.5975, pL = 1000.0, rhoR = 1.0, uR = -19.5975, pR = 0.01, x0 = 0.8, t = 0.012, p_star = 460.894, u_star = 0.0))

Visualisation — All Five Tests

fig = Figure(fontsize = 18, size = (1500, 900))
for (idx, r) in enumerate(results)
    row = (idx - 1) ÷ 3 + 1
    col = (idx - 1) % 3 + 1
    ax = Axis(
        fig[row, col], xlabel = "x", ylabel = L"\rho",
        title = "Test $(idx): $(r.ref.name)"
    )
    rho = [r.W[i][1] for i in eachindex(r.W)]
    lines!(ax, r.x, rho, color = :blue, linewidth = 1.5)
end
resize_to_layout!(fig)
fig
Example block output

Quantitative Verification

For each test, find the numerical pressure in the star region (the plateau between the waves) and compare to the exact value. The contact travels at speed u_star from the initial discontinuity x0, so we sample near x0 + u_star * t.

function measure_star_pressure(W, x, x_contact)
    # Sample cells near the contact location
    i_centre = argmin(abs.(x .- x_contact))
    # Average pressure over a few cells near the contact
    i_lo = max(1, i_centre - 5)
    i_hi = min(length(x), i_centre + 5)
    return sum(W[i][3] for i in i_lo:i_hi) / (i_hi - i_lo + 1)
end

star_pressures = [
    measure_star_pressure(r.W, r.x, r.ref.x0 + r.ref.u_star * r.ref.t) for r in results
]
exact_pressures = [r.ref.p_star for r in results]
rel_errors = [abs(star_pressures[i] - exact_pressures[i]) / max(abs(exact_pressures[i]), 1.0e-10) for i in eachindex(results)]
5-element Vector{Float64}:
 2.5912862641772868e-6
 1.0415515139224878
 2.0382073412819276e-5
 1.1103203473432732e-5
 1.1684046851510949e-5

Test Assertions

Star-region pressure should match published values within tolerance. Test 2 (near-vacuum, p_star ≈ 0.002) uses absolute tolerance since the relative error is dominated by numerics at near-zero pressure.

for i in eachindex(results)
end

Just the code

An uncommented version of this example is given below. You can view the source code for this file here.

using FiniteVolumeMethod
using StaticArrays
using CairoMakie

gamma = 1.4
eos = IdealGasEOS(gamma)
law = EulerEquations{1}(eos)
N = 400

toro_tests = [
    (
        name = "Sod", rhoL = 1.0, uL = 0.0, pL = 1.0,
        rhoR = 0.125, uR = 0.0, pR = 0.1,
        x0 = 0.5, t = 0.2, p_star = 0.30313, u_star = 0.92745,
    ),
    (
        name = "123 (two rarefactions)", rhoL = 1.0, uL = -2.0, pL = 0.4,
        rhoR = 1.0, uR = 2.0, pR = 0.4,
        x0 = 0.5, t = 0.15, p_star = 0.00189, u_star = 0.0,
    ),
    (
        name = "Woodward-Colella blast", rhoL = 1.0, uL = 0.0, pL = 1000.0,
        rhoR = 1.0, uR = 0.0, pR = 0.01,
        x0 = 0.5, t = 0.012, p_star = 460.894, u_star = 19.5975,
    ),
    (
        name = "Two-shock collision", rhoL = 5.99924, uL = 19.5975, pL = 460.894,
        rhoR = 5.99242, uR = -6.19633, pR = 46.095,
        x0 = 0.4, t = 0.035, p_star = 1691.64, u_star = 8.68975,
    ),
    (
        name = "Stationary contact", rhoL = 1.0, uL = -19.5975, pL = 1000.0,
        rhoR = 1.0, uR = -19.5975, pR = 0.01,
        x0 = 0.8, t = 0.012, p_star = 460.894, u_star = 0.0,
    ),
]

results = map(toro_tests) do tt
    wL = SVector(tt.rhoL, tt.uL, tt.pL)
    wR = SVector(tt.rhoR, tt.uR, tt.pR)
    ic(x) = x < tt.x0 ? wL : wR
    mesh = StructuredMesh1D(0.0, 1.0, N)
    prob = HyperbolicProblem(
        law, mesh, HLLCSolver(), CellCenteredMUSCL(MinmodLimiter()),
        DirichletHyperbolicBC(wL), DirichletHyperbolicBC(wR), ic;
        final_time = tt.t, cfl = 0.4,
    )
    x, U, t_end = solve_hyperbolic(prob)
    W = [conserved_to_primitive(law, U[i]) for i in eachindex(U)]
    (x = x, W = W, t = t_end, ref = tt)
end

fig = Figure(fontsize = 18, size = (1500, 900))
for (idx, r) in enumerate(results)
    row = (idx - 1) ÷ 3 + 1
    col = (idx - 1) % 3 + 1
    ax = Axis(
        fig[row, col], xlabel = "x", ylabel = L"\rho",
        title = "Test $(idx): $(r.ref.name)"
    )
    rho = [r.W[i][1] for i in eachindex(r.W)]
    lines!(ax, r.x, rho, color = :blue, linewidth = 1.5)
end
resize_to_layout!(fig)
fig

function measure_star_pressure(W, x, x_contact)
    # Sample cells near the contact location
    i_centre = argmin(abs.(x .- x_contact))
    # Average pressure over a few cells near the contact
    i_lo = max(1, i_centre - 5)
    i_hi = min(length(x), i_centre + 5)
    return sum(W[i][3] for i in i_lo:i_hi) / (i_hi - i_lo + 1)
end

star_pressures = [
    measure_star_pressure(r.W, r.x, r.ref.x0 + r.ref.u_star * r.ref.t) for r in results
]
exact_pressures = [r.ref.p_star for r in results]
rel_errors = [abs(star_pressures[i] - exact_pressures[i]) / max(abs(exact_pressures[i]), 1.0e-10) for i in eachindex(results)]

for i in eachindex(results)
end

This page was generated using Literate.jl.