GitHub

Limiter Comparison

MUSCL reconstruction achieves second-order accuracy by extrapolating cell-centred values to face centres. The slope limiter controls how aggressively this extrapolation is performed near discontinuities. A more aggressive limiter resolves smooth features better but can introduce oscillations near shocks; a more dissipative limiter is robust but smears sharp features.

This tutorial compares all five available limiters on the Sod shock tube to visualise these trade-offs.

Problem Setup

We use the standard 1D Sod shock tube with $N = 200$ cells:

using FiniteVolumeMethod
using StaticArrays

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

wL = SVector(1.0, 0.0, 1.0)
wR = SVector(0.125, 0.0, 0.1)

N = 200
mesh = StructuredMesh1D(0.0, 1.0, N)
bc_left = DirichletHyperbolicBC(wL)
bc_right = DirichletHyperbolicBC(wR)
ic(x) = x < 0.5 ? wL : wR
ic (generic function with 1 method)

Running with Each Limiter

We solve the same problem with all five limiters: Minmod (most dissipative), Superbee (least dissipative), Van Leer, Koren, and Ospre.

limiters = [
    ("Minmod", MinmodLimiter()),
    ("Superbee", SuperbeeLimiter()),
    ("Van Leer", VanLeerLimiter()),
    ("Koren", KorenLimiter()),
    ("Ospre", OspreLimiter()),
]

results = map(limiters) do (name, lim)
    prob = HyperbolicProblem(
        law, mesh, HLLCSolver(), CellCenteredMUSCL(lim),
        bc_left, bc_right, ic;
        final_time = 0.2, cfl = 0.5
    )
    x, U, t = solve_hyperbolic(prob)
    rho = [conserved_to_primitive(law, U[i])[1] for i in eachindex(U)]
    (name = name, x = x, rho = rho)
end
5-element Vector{@NamedTuple{name::String, x::Vector{Float64}, rho::Vector{Float64}}}:
 (name = "Minmod", x = [0.0025, 0.0075, 0.0125, 0.0175, 0.0225, 0.0275, 0.0325, 0.0375, 0.0425, 0.0475  …  0.9525, 0.9575, 0.9625, 0.9675, 0.9725, 0.9775, 0.9825, 0.9875, 0.9925, 0.9975], rho = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125])
 (name = "Superbee", x = [0.0025, 0.0075, 0.0125, 0.0175, 0.0225, 0.0275, 0.0325, 0.0375, 0.0425, 0.0475  …  0.9525, 0.9575, 0.9625, 0.9675, 0.9725, 0.9775, 0.9825, 0.9875, 0.9925, 0.9975], rho = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125])
 (name = "Van Leer", x = [0.0025, 0.0075, 0.0125, 0.0175, 0.0225, 0.0275, 0.0325, 0.0375, 0.0425, 0.0475  …  0.9525, 0.9575, 0.9625, 0.9675, 0.9725, 0.9775, 0.9825, 0.9875, 0.9925, 0.9975], rho = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125])
 (name = "Koren", x = [0.0025, 0.0075, 0.0125, 0.0175, 0.0225, 0.0275, 0.0325, 0.0375, 0.0425, 0.0475  …  0.9525, 0.9575, 0.9625, 0.9675, 0.9725, 0.9775, 0.9825, 0.9875, 0.9925, 0.9975], rho = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125, 0.125])
 (name = "Ospre", x = [0.0025, 0.0075, 0.0125, 0.0175, 0.0225, 0.0275, 0.0325, 0.0375, 0.0425, 0.0475  …  0.9525, 0.9575, 0.9625, 0.9675, 0.9725, 0.9775, 0.9825, 0.9875, 0.9925, 0.9975], rho = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  0.12499999999999989, 0.12499999999999951, 0.12499999999999968, 0.1249999999999998, 0.12499999999999968, 0.12499999999999989, 0.12499999999999967, 0.1249999999999998, 0.12499999999999985, 0.12499999999999985])

We also include a first-order (no reconstruction) baseline:

prob_first = HyperbolicProblem(
    law, mesh, HLLCSolver(), NoReconstruction(),
    bc_left, bc_right, ic;
    final_time = 0.2, cfl = 0.5
)
x_first, U_first, _ = solve_hyperbolic(prob_first)
rho_first = [conserved_to_primitive(law, U_first[i])[1] for i in eachindex(U_first)]
200-element Vector{Float64}:
 0.999999999996201
 0.9999999999919227
 0.999999999982993
 0.9999999999645519
 0.9999999999268738
 0.999999999850709
 0.9999999996984046
 0.9999999993971627
 0.9999999988078995
 0.999999997668061
 ⋮
 0.12500000125531985
 0.1250000004779681
 0.12500000018148116
 0.12500000006870016
 0.12500000002592285
 0.12500000000974779
 0.12500000000365216
 0.1250000000013627
 0.12500000000050648

Exact Solution

function sod_exact(x, t; x0 = 0.5, gamma = 1.4)
    cL = sqrt(gamma * 1.0 / 1.0)
    P_star = 0.30313017805064707
    v_star = 0.92745262004895057
    rho_star_L = 0.42631942817849544
    rho_star_R = 0.26557371170530708
    c_star_L = sqrt(gamma * P_star / rho_star_L)
    S_shock = v_star + 1.0 / (gamma * rho_star_R) *
        (P_star - 0.1) / (v_star + 1.0e-30)
    xi = (x - x0) / t
    if xi < -cL
        return 1.0
    elseif xi < v_star - c_star_L
        v_fan = 2.0 / (gamma + 1) * (cL + xi)
        c_fan = cL - 0.5 * (gamma - 1) * v_fan
        return 1.0 * (c_fan / cL)^(2.0 / (gamma - 1))
    elseif xi < v_star
        return rho_star_L
    elseif xi < S_shock
        return rho_star_R
    else
        return 0.125
    end
end
sod_exact (generic function with 1 method)

Visualisation

using CairoMakie

x_exact = range(0.0, 1.0, length = 1000)
rho_exact = [sod_exact(xi, 0.2) for xi in x_exact]

colors = [:blue, :red, :green, :orange, :purple]

fig = Figure(fontsize = 20, size = (1200, 800))
ax1 = Axis(fig[1, 1], xlabel = "x", ylabel = L"\rho", title = "Full Domain")
lines!(ax1, x_exact, rho_exact, color = :black, linewidth = 2, label = "Exact")
scatter!(ax1, x_first, rho_first, color = :gray, markersize = 3, label = "1st order")
for (i, r) in enumerate(results)
    scatter!(ax1, r.x, r.rho, color = colors[i], markersize = 3, label = r.name)
end
axislegend(ax1, position = :cb)
Legend()

Zoom into the contact discontinuity region:

ax2 = Axis(
    fig[2, 1], xlabel = "x", ylabel = L"\rho",
    title = "Contact Discontinuity (zoom)"
)
lines!(ax2, x_exact, rho_exact, color = :black, linewidth = 2)
scatter!(ax2, x_first, rho_first, color = :gray, markersize = 5, label = "1st order")
for (i, r) in enumerate(results)
    scatter!(ax2, r.x, r.rho, color = colors[i], markersize = 5, label = r.name)
end
xlims!(ax2, 0.55, 0.85)
ylims!(ax2, 0.15, 0.5)
axislegend(ax2, position = :rt)

resize_to_layout!(fig)
fig
Example block output

Observations

  • Superbee resolves the contact discontinuity most sharply, but can produce slight overshoots near the shock.
  • Minmod is the most dissipative, producing the smoothest profiles.
  • Van Leer, Koren, and Ospre offer a middle ground between sharpness and robustness.
  • The first-order solution (no reconstruction) is very diffusive, illustrating why MUSCL reconstruction is important.
true

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

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

wL = SVector(1.0, 0.0, 1.0)
wR = SVector(0.125, 0.0, 0.1)

N = 200
mesh = StructuredMesh1D(0.0, 1.0, N)
bc_left = DirichletHyperbolicBC(wL)
bc_right = DirichletHyperbolicBC(wR)
ic(x) = x < 0.5 ? wL : wR

limiters = [
    ("Minmod", MinmodLimiter()),
    ("Superbee", SuperbeeLimiter()),
    ("Van Leer", VanLeerLimiter()),
    ("Koren", KorenLimiter()),
    ("Ospre", OspreLimiter()),
]

results = map(limiters) do (name, lim)
    prob = HyperbolicProblem(
        law, mesh, HLLCSolver(), CellCenteredMUSCL(lim),
        bc_left, bc_right, ic;
        final_time = 0.2, cfl = 0.5
    )
    x, U, t = solve_hyperbolic(prob)
    rho = [conserved_to_primitive(law, U[i])[1] for i in eachindex(U)]
    (name = name, x = x, rho = rho)
end

prob_first = HyperbolicProblem(
    law, mesh, HLLCSolver(), NoReconstruction(),
    bc_left, bc_right, ic;
    final_time = 0.2, cfl = 0.5
)
x_first, U_first, _ = solve_hyperbolic(prob_first)
rho_first = [conserved_to_primitive(law, U_first[i])[1] for i in eachindex(U_first)]

function sod_exact(x, t; x0 = 0.5, gamma = 1.4)
    cL = sqrt(gamma * 1.0 / 1.0)
    P_star = 0.30313017805064707
    v_star = 0.92745262004895057
    rho_star_L = 0.42631942817849544
    rho_star_R = 0.26557371170530708
    c_star_L = sqrt(gamma * P_star / rho_star_L)
    S_shock = v_star + 1.0 / (gamma * rho_star_R) *
        (P_star - 0.1) / (v_star + 1.0e-30)
    xi = (x - x0) / t
    if xi < -cL
        return 1.0
    elseif xi < v_star - c_star_L
        v_fan = 2.0 / (gamma + 1) * (cL + xi)
        c_fan = cL - 0.5 * (gamma - 1) * v_fan
        return 1.0 * (c_fan / cL)^(2.0 / (gamma - 1))
    elseif xi < v_star
        return rho_star_L
    elseif xi < S_shock
        return rho_star_R
    else
        return 0.125
    end
end

using CairoMakie

x_exact = range(0.0, 1.0, length = 1000)
rho_exact = [sod_exact(xi, 0.2) for xi in x_exact]

colors = [:blue, :red, :green, :orange, :purple]

fig = Figure(fontsize = 20, size = (1200, 800))
ax1 = Axis(fig[1, 1], xlabel = "x", ylabel = L"\rho", title = "Full Domain")
lines!(ax1, x_exact, rho_exact, color = :black, linewidth = 2, label = "Exact")
scatter!(ax1, x_first, rho_first, color = :gray, markersize = 3, label = "1st order")
for (i, r) in enumerate(results)
    scatter!(ax1, r.x, r.rho, color = colors[i], markersize = 3, label = r.name)
end
axislegend(ax1, position = :cb)

ax2 = Axis(
    fig[2, 1], xlabel = "x", ylabel = L"\rho",
    title = "Contact Discontinuity (zoom)"
)
lines!(ax2, x_exact, rho_exact, color = :black, linewidth = 2)
scatter!(ax2, x_first, rho_first, color = :gray, markersize = 5, label = "1st order")
for (i, r) in enumerate(results)
    scatter!(ax2, r.x, r.rho, color = colors[i], markersize = 5, label = r.name)
end
xlims!(ax2, 0.55, 0.85)
ylims!(ax2, 0.15, 0.5)
axislegend(ax2, position = :rt)

resize_to_layout!(fig)
fig

This page was generated using Literate.jl.