18) Hyperbolic Systems

18) Hyperbolic Systems#

Last time:

Godunov’s Theorem
Slope reconstruction and limiting
Common limiters in the FV literature
Introduction on hyperbolic systems
Example: Isentropic gas dynamics

Today:

Hyperbolic systems
1.1 Riemann problem for systems: shocks
1.2 Admissible shocks
Generalized Riemann invariants
Exact Riemann solver for isothermal gas dynamics
Approximate Riemann solvers
4.1 HLL (Harten, Lax, and van Leer)

using LinearAlgebra
using Plots
default(linewidth=3)

# utility functions seen in class

struct RKTable
    A::Matrix
    b::Vector
    c::Vector
    function RKTable(A, b)
        s = length(b)
        A = reshape(A, s, s)
        c = vec(sum(A, dims=2))
        new(A, b, c)
    end
end

rk4 = RKTable([0 0 0 0; .5 0 0 0; 0 .5 0 0; 0 0 1 0], [1, 2, 2, 1] / 6)

# explicit RK solver
function ode_rk_explicit(f, u0; tfinal=1., h=0.1, table=rk4)
    u = copy(u0)
    t = 0.
    n, s = length(u), length(table.c)
    fY = zeros(n, s)
    thist = [t]
    uhist = [u0]
    while t < tfinal
        tnext = min(t+h, tfinal)
        h = tnext - t
        for i in 1:s
            ti = t + h * table.c[i]
            Yi = u + h * sum(fY[:,1:i-1] * table.A[i,1:i-1], dims=2)
            fY[:,i] = f(ti, Yi)
        end
        u += h * fY * table.b
        t = tnext
        push!(thist, t)
        push!(uhist, u)
    end
    thist, hcat(uhist...)
end

function testfunc(x)
    max(1 - 4*abs.(x+2/3),
        abs.(x) .< .2,
        (2*abs.(x-2/3) .< .5) * cospi(2*(x-2/3)).^2
    )
end

flux_advection(u) = u
flux_burgers(u) = u^2/2
flux_traffic(u) = u * (1 - u)

riemann_advection(uL, uR) = 1*uL # velocity is +1

function fv_solve1(riemann, u_init, n, tfinal=1)
    h = 2 / n # 1D domain [-1,1]
    x = LinRange(-1+h/2, 1-h/2, n) # cell midpoints (centroids)
    idxL = 1 .+ (n-1:2*n-2) .% n # left interface indices
    idxR = 1 .+ (n+1:2*n) .% n # right interface indices
    function rhs(t, u)
        fluxL = riemann(u[idxL], u)
        fluxR = riemann(u, u[idxR])
        (fluxL - fluxR) / h
    end
    thist, uhist = ode_rk_explicit(rhs, u_init.(x), h=h, tfinal=tfinal)
    x, thist, uhist
end

function riemann_burgers(uL, uR)
    flux = zero(uL)
    for i in 1:length(flux)
        fL = flux_burgers(uL[i])
        fR = flux_burgers(uR[i])
        flux[i] = if uL[i] > uR[i] # shock
            max(fL, fR)
        elseif uL[i] > 0 # rarefaction all to the right
            fL
        elseif uR[i] < 0 # rarefaction all to the left
            fR
        else
            0
        end
    end
    flux
end

limit_zero(r) = 0
limit_none(r) = 1
limit_minmod(r) = max(min(2*r, 2*(1-r)), 0)
limit_sin(r) = (0 < r && r < 1) * sinpi(r)
limit_vl(r) = max(4*r*(1-r), 0)
limit_bj(r) = max(0, min(1, 4*r, 4*(1-r)))
limiters = [limit_zero limit_none limit_minmod limit_sin limit_vl limit_bj];

1. Hyperbolic systems#

Systems of the form

\[ \mathbf U_t + \nabla \cdot F(\mathbf U) = 0 \]

It would be too good if we could safely reconstruct each component of the flux in the same way.

But the system has waves traveling at different speeds. So we need to be more careful with reconstructions.

As in the scalar case, elementary waves may propagate as shock waves, or as rarefaction fans, except now they may occur in each of the different characteristic fields.

Isentropic gas dynamics#

(41)#\[\begin{align} U &= \begin{bmatrix} \rho \\ \rho u \end{bmatrix} \, ,& F(U) &= \begin{bmatrix} \rho u \\ \rho u^2 + p \end{bmatrix} \end{align}\]

Variable	meaning
\(\rho\)	density
\(u\)	velocity
\(\rho u\)	momentum
\(p\)	pressure

Isentrtropic equation of state (pressure depends only on density) \(p(\rho) = C \rho^\gamma\) with \(\gamma = 1.4\) (typical air).
“isothermal” gas dynamics: \(p(\rho) = c^2 \rho\), wave speed \(c\).
Compute the needed \(\rho u^2\) as \( \rho u^2 = \frac{(\rho u)^2}{\rho} .\)

The flux Jacobian is

\[ F'(U) = \frac{\partial F}{\partial U} \]

Smooth wave structure#

For perturbations of a constant state, systems of equations produce multiple waves with speeds equal to the eigenvalues of the flux Jacobian,

\[\begin{split} F'(U) = \begin{bmatrix} 0 & 1 \\ -u^2 + c^2 & 2 u \end{bmatrix}. \end{split}\]

We can compute the eigenvalues \(\Lambda \) and eigenvector \(W\) matrices,

\[ F'(U) = W \Lambda W^{-1} \]

as

(42)#\[\begin{align} W &= \begin{bmatrix} 1 & 1 \\ u-c & u+c \end{bmatrix} & \Lambda &= \begin{bmatrix} u-c & \\ & u+c \end{bmatrix} . \end{align}\]

1.1 Riemann problem for systems: shocks#

Given states \(U_L\) and \(U_R\), we will see two waves with some state \(U^*\) in between.

There could be two shocks, two rarefactions, or one of each.
The type of wave will determine the condition that must be satisfied to connect \(U_L\) to \(U_*\) and \(U_*\) to \(U_R\).

Left-moving wave#

If there is a shock between \(U_L\) and \(U_*\), the Rankine-Hugoniot condition

\[ S \Delta U = \Delta F \]

will be satisfied along with the entropy condition

\[ \lambda_i(U_L) \ge S \ge \lambda_i(U_*) \]

where \(\lambda_i(U)\) are the corresponding eigenvalues of the flux Jacobian \(F'(U)\). These are the characteristic speeds.

This condition states:

The shock speed lies between the characteristic speeds on each side.

The inequality in the entropy condition is strict if the wave is genuinely nonlinear (more on this in a second). Meaning, \( \lambda_i(U_L) \ge S \ge \lambda_i(U_*) \).

For isentropic/isothermal gas dynamics,

(43)#\[\begin{align} U &= \begin{bmatrix} \rho \\ \rho u \end{bmatrix} \, ,& F(U) &= \begin{bmatrix} \rho u \\ \rho u^2 + c^2 \rho \end{bmatrix} \end{align}\]

If a shock connects \(U_L = (\rho_L, \rho_L u_L)\) to \(U_* = (\rho_* , \rho_* u_*)\), with shock speed \(S\), and

\[\begin{align*} S = \frac{\Delta F}{\Delta U} \end{align*}\]

in this case we have

(44)#\[\begin{align} S \left( \begin{bmatrix} \rho_* \\ \rho_* u_* \end{bmatrix} - \begin{bmatrix} \rho_L \\ \rho_L u_L \end{bmatrix} \right) &= \begin{bmatrix} \rho_* u_* \\ \rho_* u_*^2 + c^2 \rho_* \end{bmatrix} - \begin{bmatrix} \rho_L u_L \\ \rho_L u_L^2 + c^2 \rho_L \end{bmatrix} \end{align}\]

the Rankine-Hugoniot condition is, component-wise:

(45)#\[\begin{align} \textrm{(i) }(\rho_* - \rho_L) S &= \rho_* u_* - \rho_L u_L \\ \textrm{(ii) }(\rho_* u_* - \rho_L u_L) S &= \Big(\rho_* u_*^2 + c^2 \rho_* \Big) - \Big( \rho_L u_L^2 + c^2 \rho_L \Big) \end{align}\]

Nonlinear Systems#

For a system

\[ \mathbf U_t + \nabla \cdot F(\mathbf U) = 0 \]

We have again, the Jacobian

\[ F'(U) = \frac{\partial F}{\partial U} \]

\(\lambda_i(U)\) are the eigenvalues (wave speeds)
\(w_i(U)\) are the corresponding eigenvectors

A field (wave family) \(i\) is called genuinely nonlinear if

\[ \nabla \lambda_i(U) \cdot w_i(U) \neq 0, \forall U \]

and we say that the \(i\)-th characteristic field is linearly degenerate if

\[ \nabla \lambda_i(U) \cdot w_i(U) = 0 \]

Observations#

These conditions generalize the scalar case \(u_t + f(u)_x = 0 \).
The genuinely nonlinear case, in the scalar case, reduces to the scalar convexity condition.
A genuinely nonlinear field admits elementary waves that are either shock waves (converging characteristics) or expansion fans (diverging characteristics)
While a linearly degenerate field admits simple waves called contact discontinuities, across which the characteristics remain parallel as in the linear case. In fact, in this case wave speed stays constant along the wave - no steepening occurs.
Again, for the genuinely nonlinear case, the entropy condition is strict \( \lambda_i(U_L) \ge S \ge \lambda_i(U_*) \).

Rankine–Hugoniot for the isentropic gas dynamics example#

Going back to the isentropic/isothermal gas dynamics example, we solve the first equation (i) for \(S\) and substituting into the second (ii):

\[\begin{split} \begin{split} \frac{(\rho_* u_* - \rho_L u_L)^2}{\rho_* - \rho_L} &= \rho_* u_*^2 - \rho_L u_L^2 + c^2 (\rho_* - \rho_L) \\ \rho_*^2 u_*^2 - 2 \rho_* \rho_L u_* u_L + \rho_L^2 u_L^2 &= \rho_* (\rho_* - \rho_L) u_*^2 - \rho_L (\rho_* - \rho_L) u_L^2 + c^2 (\rho_* - \rho_L)^2 \\ \rho_* \rho_L \Big( u_*^2 - 2 u_* u_L + u_L^2 \Big) &= c^2 (\rho_* - \rho_L)^2 \\ (u_* - u_L)^2 &= c^2 \frac{(\rho_* - \rho_L)^2}{\rho_* \rho_L} \\ u_* - u_L &= \pm c \frac{\rho_* - \rho_L}{\sqrt{\rho_* \rho_L}} \end{split} \end{split}\]

By doing so, we have found the Hugoniot locus, the algebraic relation that admissible shock-connected states \(U_*\) must satisfy.

Note that algebraically, there are two possible solutions (\(\pm\) sign) and we will need to use the entropy condition to learn which sign to take. One corresponds to a physical shock, the other corresponds to a nonphysical solution (e.g. expansion shock). The Rankine–Hugoniot condition(s) alone does not distinguish them. That’s why we also need the entropy condition.

Admissible shocks#

We need to choose the sign

\[u_* - u_L = \pm c \frac{\rho_* - \rho_L}{\sqrt{\rho_*\rho_L}}\]

The entropy condition requires that

\[ \lambda(U_L) = u_L - c \ge S \ge u_* - c = \lambda(U_*) . \]

Meanwhile, the first Rankine-Hugoniot condition (i) and \(\rho > 0\) provides

\[ \rho_L c \le \rho_L (u_L - S) = \rho_* (u_* - S) \le \rho_* c \]

which is to say that \(\rho_* \ge \rho_L\) and consequently we take the negative sign,

\[ u_* - u_L = - c\frac{\rho_* - \rho_L}{\sqrt{\rho_* \rho_L}} .\]

Right-moving wave#

The same process for the right wave with characteristic speed \(\lambda(U) = u + c\) yields a shock when \(\rho_* \ge \rho_R\), in which case the velocity jump is

\[ u_R - u_* = - c \frac{\rho_* - \rho_R}{\sqrt{\rho_* \rho_R}} . \]

Rarefactions#

A rarefaction occurs when the entropy condition is not satisfied, that is

\[ \lambda_i(U_L) < \lambda_i(U_R) \]

2. Generalized Riemann invariants#

For a system, with eigenvalues \(\lambda_i(U)\) and corresponding eigenvectors \(w_i(U)\), a Riemann invariant for the \(i\)-th field is a scalar function \(r(U)\) such that

\[ \nabla r(U) \cdot w_i(U) = 0 \]

(Derivation of this condition is beyond the scope of this class.)

This is a first-order linear PDE, whose solution gives invariants. \(r\) does not change along the \(i\)-th characteristic direction (hence, the invariant name).

Why are they useful?

They let you reduce a system to something like a scalar problem along each wave family.

In particular:

Across a rarefaction wave, Riemann invariants are constant
They give you explicit relations between states

We will use:

\[ \frac{d U^1}{w_i^1} = \frac{d U^2}{w_i^2} = \dotsb \]

Isentropic/isothermal gas dynamics#

\[ \frac{d \rho}{1} = \frac{d (\rho u)}{u - c} \]

across the wave \(\lambda = u-c\). We can rearrange to

\[ (u-c) d \rho = \rho d u + u d \rho \]

or

\[ d u + \frac{c}{\rho} d \rho = 0 . \]

Integration yields

\[ u + \int \frac{c}{\rho} d \rho = u + c \log \rho = \text{constant} . \]

So if we have a rarefaction between \(U_L\) and \(U_*\) then

\[ u_L + c \log \rho_L = u_* + c \log \rho_* . \]

Similarly, for the wave \(\lambda = u + c\), a rarefaction between \(U_*\) and \(U_R\) results in

\[ u_* - c \log \rho_* = u_R - c \log \rho_R . \]

If the rarefaction is sonic, we have \(u_0 - c = 0\) for a left rarefaction (or \(u_0 + c =0\) for the right) and can compute \(\rho_0\) from

\[ u_* -c \log \rho_* = u_0 - c \log \rho_0,\]

via

\[ \rho_0 = \exp\big[(u_0 - u_*)/c + \log \rho_* \big].\]

First a miracle happens: in this case it is exactly integrable
In general we will need to use a Newton method to solve for the state in the star region.

Basic algorithm#

Find \(\rho_*\)
Use entropy condition (for shock speeds)

If it’s a shock: find \(u_*\) using Rankine-Hugoniot
If it’s a rarefaction: use generalized Riemann invariant

3. Exact Riemann solver for isothermal gas dynamics#

function flux_isogas(U, c=1)
    rho = U[1]
    u = U[2] / rho
    [U[2], U[2]*u + c^2*rho]
end

function ujump_isogas(rho_L, rho_R, cw)
    if sign(rho_L - rho_R) == sign(cw) # shock
        sign(cw) * (rho_R - rho_L) / sqrt(rho_L*rho_R)
    else # rarefaction
        cw * (log(rho_R) - log(rho_L))
    end
end

function dujump_isogas(rho_L, drho_L, rho_R, drho_R, cw)
    if sign(rho_L - rho_R) == sign(cw) # shock
        sign(cw) * ((drho_R - drho_L) / sqrt(rho_L*rho_R)
            - .5*(rho_R - rho_L) * (rho_L*rho_R)^(-3/2) 
        * (drho_L * rho_R + rho_L * drho_R))
    else
        cw * (drho_R / rho_R - drho_L / rho_L)
    end
end

dujump_isogas (generic function with 1 method)

function riemann_isogas(UL, UR, maxit=20; show=false)
    rho_L, u_L = UL[1], UL[2]/UL[1]
    rho_R, u_R = UR[1], UR[2]/UR[1]
    rho = .5 * (rho_L + rho_R) # initial guess
    U_star = zero.(UL)
    for i in 1:maxit
        f = (ujump_isogas(UL[1], rho, -1)
            + ujump_isogas(rho, UR[1], 1)
            - (UR[2]/UR[1] - UL[2]/UL[1]))
        if norm(f) < 1e-10
            u = u_L + ujump_isogas(UL[1], rho, -1)
            U_star[:] = [rho, rho*u]
            break
        end
        J = (dujump_isogas(rho_L, 0, rho, 1, -1)
            + dujump_isogas(rho, 1, rho_R, 0, 1))
        delta_rho = -f / J
        rho += delta_rho
        ## Line search not needed in practice
        #while rho + delta_rho <= 0
        #    delta_rho /= 2 # line search to prevent negative rho
        #end
    end
    U0 = resolve_isogas(UL, UR, U_star)
    if show; @show U0, U0[2]/U0[1]; end
    flux_isogas(U0)
end

riemann_isogas (generic function with 2 methods)

Resolving waves for isothermal gas dynamics#

function resolve_isogas(UL, UR, U_star; c=1)
    rho_L, u_L = UL[1], UL[2]/UL[1]
    rho_R, u_R = UR[1], UR[2]/UR[1]
    rho, u = U_star[1], U_star[2] / U_star[1]
    if ((u_L - c < 0 < u - c) || 
        (u + c < 0 < u_R + c))
        # inside left (right) sonic rarefaction
        u0 = sign(u) * c
        rho0 = exp((u0-u)/c + log(rho))
        @show "sonic"
        [rho0, rho0 * u0]
    elseif ((rho_L >= rho && 0 <= u_L - c) ||
            (rho_L < rho && 0 < rho*u - UL[2]))
        # left rarefaction or shock is supersonic
        UL
    elseif ((rho_R >= rho && u_R + c <= 0) ||
            (rho_R < rho && rho*u - UR[2] < 0))
        # right rarefaction or shock is supersonic
        UR
    else # sample star region
        U_star
    end
end

resolve_isogas (generic function with 1 method)

Test the Riemann solver#

isogas(rho, u) = [rho, rho*u]
riemann_isogas(isogas(1, .1), isogas(.9, .1), show=true)

(U0, U0[2] / U0[1]) = ([0.9486804089487428, 0.14484765853290457], 0.15268330321421308)

2-element Vector{Float64}:
 0.14484765853290457
 0.9707962279163911

function initial_isogas(x)
    f(s) = isogas(1 + 2 * exp(-4s ^ 2), 0.1)
    vcat(f.(x)'...)
end

x = LinRange(-1, 1, 100)
U = initial_isogas(x)
plot(x, [U (U[:, 2] ./ U[:, 1])], 
  label=["density" "momentum" "velocity"])

Solver#

function fv_solve2system(riemann, initfunc, n, tfinal=1; dt_scale=1, limit=limit_sin)
    h = 2 / n
    x = LinRange(-1+h/2, 1-h/2, n) # cell midpoints (centroids)
    idxL = 1 .+ (n-1:2*n-2) .% n
    idxR = 1 .+ (n+1:2*n) .% n
    U0 = initfunc(x)
    n, k = size(U0)
    function rhs(t, U)
        U = reshape(U, n, k)
        jump = U[idxR, :] - U[idxL, :]
        r = (U - U[idxL, :]) ./ jump
        r[.~isfinite.(r)] .= 0
        g = limit.(r) .* jump / 2h
        flux = zero.(U)
        for i in 1:n
            UL = U[idxL[i],:] + g[idxL[i],:] * h/2
            UR = U[i,:] - g[i,:] * h/2
            flux[i, :] = riemann(UL, UR)
        end
        vec(flux - flux[idxR, :]) / h
    end
    thist, uhist = ode_rk_explicit(
        rhs, vec(U0), h=h*dt_scale, tfinal=tfinal)
    x, thist, reshape(uhist, n, k, length(thist))
end

fv_solve2system (generic function with 2 methods)

x, t_hist, U_hist = fv_solve2system(riemann_isogas,
initial_isogas, 200, 1.5, dt_scale=0.8)
rho_hist = U_hist[:, 1, :]
u_hist = U_hist[:, 2, :] ./ rho_hist
plot(x, rho_hist[:, 1:20:end], 
label=round.(t_hist[1:20:end]', digits=2)) # check both components

We can see the initial smooth solutions develop shock and rarefaction structure.

Visualize solutions#

function initial_isogas_2(x)
    f(s) = isogas(1 + 2 * (abs(s) < 0.25), 0.0)
    vcat(f.(x)'...)
end
x, t_hist, U_hist = fv_solve2system(riemann_isogas,
initial_isogas_2, 200, 1, dt_scale=0.5, limit=limit_sin)
step = length(t_hist) ÷ 5; rho_hist = U_hist[:, 1, :]
plot(x, rho_hist[:, 1:step:end], label=round.(t_hist[1:step:end]', digits=2), title="density")

u_hist = U_hist[:, 2, :] ./ rho_hist
plot(x, u_hist[:, 1:step:end], label=round.(t_hist[1:step:end]', digits=2), title="velocity")

4. Approximate Riemann solvers#

Exact Riemann solvers are

complicated to implement
fragile in the sense that small changes to the physics, such as in the equation of state \(p(\rho)\), can require changing many conditionals
the need to solve for \(\rho^*\) using a Newton method and then evaluate each case of the wave structure can be expensive.
An exact Riemann solver has never been implemented for some equations.

4.1 HLL (Harten, Lax, and van Leer)#

Assume two shocks with speeds \(S_L\) and \(S_R\).

These speeds will be estimated and must be at least as fast as the fastest left- and right-traveling waves.

If the wave speeds \(S_L\) and \(S_R\) are given, we have the Rankine-Hugoniot conditions across both shocks,

(46)#\[\begin{align} S_L (U_* - U_L) &= f(U_*) - f(U_L) \\ S_R (U_R - U_*) &= f(U_R) - f(U_*) . \end{align}\]

Adding these together gives

\[S_R U_R - (S_R - S_L) U_* - S_L U_L = f(U_R) - f(U_L)\]

which can be solved for \(U_*\) as

\[ U_* = \frac{f(U_L) - f(U_R) - S_L U_L + S_R U_R}{S_R - S_L} . \]

We can save an extra evaluation of the flux by using the Rankine-Hugoniot conditions

\[ f(U_*) = \frac{S_R f(U_L) - S_L f(U_R) + S_L S_R (U_R - U_L)}{S_R - S_L} . \]

Isentropic/isothermal gas dynamics example#

(47)#\[\begin{align} S_L &= \min( u_L - c, u_R - c ) \\ S_R &= \max(u_L + c, u_R + c) \end{align}\]

Rusanov method#

Special case \(S_L = - S_R\), in which case the wave structure is always subsonic and the flux is simply

\[f(U_*) = \frac 1 2 \Big( f(U_L) + f(U_R) \Big) - \frac S 2 \Big( U_R - U_L \Big) . \]

Observations on HLL solvers#