2023-12-04 Reconstruction#

Last time#

  • Notes on unstructured meshing workflow

  • Finite volume methods for hyperbolic conservation laws

  • Riemann solvers for scalar equations

    • Shocks and the Rankine-Hugoniot condition

    • Rarefactions and entropy solutions

Today#

  • Godunov’s Theorem

  • Slope reconstruction and limiting

  • Hyperbolic systems

    • Rankine-Hugoniot and Riemann Invariants

    • Basic strategy for exact Riemann solvers

using LinearAlgebra
using Plots
default(linewidth=3)

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)

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
    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
    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

function riemann_traffic(uL, uR)
    flux = zero(uL)
    for i in 1:length(flux)
        fL = flux_traffic(uL[i])
        fR = flux_traffic(uR[i])
        flux[i] = if uL[i] < uR[i] # shock
            min(fL, fR)
        elseif uL[i] < .5 # rarefaction all to the right
            fL
        elseif uR[i] > .5 # rarefaction all to the left
            fR
        else
            flux_traffic(.5)
        end
    end
    flux 
end

riemann_traffic (generic function with 1 method)

Godunov methods (first order accurate)#

init_func(x) = testfunc(x) - .5
x, thist, uhist = fv_solve1(riemann_burgers, init_func, 200, 1)
plot(x, uhist[:,1:20:end], legend=:none)

Burgers#

  • flux \(u^2/2\) has speed \(u\)

  • negative values make sense

  • satisfies a maximum principle

Traffic#

  • flux \(u - u^2\) has speed \(1 - 2u\)

  • state must satisfy \(u \in [0, 1]\)

Godunov’s Theorem (1954)#

Linear numerical methods

\[ \dot u_i = \sum_j a_{ij} u_j \]
for solving time-dependent PDE can be at most first order accurate if they are monotone.

For our purposes, monotonicity is equivalent to positivity preservation,

\[ \min_x u(0, x) \ge 0 \quad \Longrightarrow \quad \min_x u(t, 0) \ge 0 .\]

Discontinuities#

A numerical method for representing a discontinuous function on a stationary grid can be no better than first order accurate in the \(L^1\) norm,

\[ \lVert u - u^* \rVert_{L^1} = \int \lvert u - u^* \rvert . \]
If we merely sample a discontinuous function, say
\[\begin{split} u(x) = \begin{cases} 0, & x \le a \\ 1, & x > a \end{cases} \end{split}\]
onto a grid with element size \(\Delta x\) then we will have errors of order 1 on an interval proportional to \(\Delta x\).

In light of these two observations, we may still ask for numerical methods that are more than first order accurate for smooth solutions, but those methods must be nonlinear.

Slope Reconstruction#

One method for constructing higher order methods is to use the state in neighboring elements to perform a conservative reconstruction of a piecewise polynomial, then compute numerical fluxes by solving Riemann problems at the interfaces. If \(x_i\) is the center of cell \(i\) and \(g_i\) is the reconstructed gradient inside cell \(i\), our reconstructed solution is

\[ \tilde u_i(x) = u_i + g_i \cdot (x - x_i) . \]
We would like this reconstruction to be monotone in the sense that
\[ \Big(\tilde u_i(x) - \tilde u_j(x) \Big) \Big( u_i - u_j \Big) \ge 0 \]
for any \(x\) on the interface between element \(i\) and element \(j\).

Question#

Is the symmetric slope

\[ \hat g_i = \frac{u_{i+1} - u_{i-1}}{2 \Delta x} \]
monotone?

Slope limiting#

We will determine gradients by “limiting” the above slope using a nonlinear function that reduces to 1 when the solution is smooth. There are many ways to express limiters and our discussion here roughly follows Berger, Aftosmis, and Murman (2005).

Berger et al, figure 3b, Cell averages

We will express a slope limiter in terms of the ratio

\[ r_i = \frac{u_i - u_{i-1}}{u_{i+1} - u_{i-1}} \]
and use as a gradient,
\[ g_i = \phi(r_i) \frac{u_{i+1} - u_{i-1}}{2 \Delta x} . \]
Our functions \(\phi\) will be zero unless \(0 < r < 1\) and \(\phi(1/2) = 1\). Berger et al, figure 3a, Symmetric form of slope limiters All of these limiters are second order accurate and TVD; those that fall below minmod are not second order accurate and those that are above Barth-Jesperson are not second order accurate, not TVD, or produce artifacts.

Common limiters#

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];

plot(limiters, label=limiters, xlims=(-.1, 1.1))

A slope-limited solver#

function fv_solve2(riemann, u_init, n, tfinal=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
    function rhs(t, u)
        jump = u[idxR] - u[idxL]
        r = (u - u[idxL]) ./ jump
        r[isnan.(r)] .= 0
        g = limit.(r) .* jump / 2h
        fluxL = riemann(u[idxL] + g[idxL]*h/2, u - g*h/2)
        fluxR = fluxL[idxR]
        (fluxL - fluxR) / h
    end
    thist, uhist = ode_rk_explicit(
        rhs, u_init.(x), h=h, tfinal=tfinal)
    x, thist, uhist
end

fv_solve2 (generic function with 3 methods)

x, thist, uhist = fv_solve2(riemann_advection, testfunc, 200, 10,
    limit_sin)
plot(x, uhist[:,1:400:end], legend=:none)

Hyperbolic systems#

\[ U_t + f(U)_x = 0 \]

Isentropic gas dynamics#

(64)#\[\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

  • Equation of state \(p(\rho) = C \rho^\gamma\) with \(\gamma = 1.4\) (typical air).

  • “isothermal” gas dynamics: \(p(\rho) = c^2 \rho\), wave speed \(c\).

  • Compute as \( \rho u^2 = \frac{(\rho u)^2}{\rho} .\)

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 and eigenvectors,
\[ f'(U) = W \Lambda W^{-1} \]
as

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

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)\) is the corresponding eigenvalue of the flux Jacobian \(f'(U)\). The inequality in the entropy condition is strict if the wave is genuinely nonlinear. For isothermal gas dynamics, the Rankine-Hugoniot condition is

(66)#\[\begin{align} (\rho_* - \rho_L) s &= \rho_* u_* - \rho_L u_L \\ (\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}\]

Solving the first equation for \(s\) and substituting into the second, we compute \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} and will need to use the entropy condition to learn which sign to take.

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 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 \(\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

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

Generalized Riemann invariants#

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

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

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].\]

Basic algorithm#

  1. Find \(\rho_*\)

  2. 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

First a miracle happens#

In general we will need to use a Newton method to solve for the state in the star region.