.ipynb .pdf

repository open issue

Binder Live Code

Contents

2021-09-15 Advection and stability

Last time

• Variable coefficients

• Conservative/divergence form vs non-divergence forms

• Verification with discontinuities

Today

• Advection and boundary layers

• The (cell) Péclet number and oscillations

using Plots
using LinearAlgebra

function vander(x, k=nothing)
    if k === nothing
        k = length(x)
    end
    V = ones(length(x), k)
    for j = 2:k
        V[:, j] = V[:, j-1] .* x
    end
    V
end

function fdstencil(source, target, k)
    "kth derivative stencil from source to target"
    x = source .- target
    V = vander(x)
    rhs = zero(x)'
    rhs[k+1] = factorial(k)
    rhs / V
end


function poisson(x, spoints, forcing; left=(0, zero), right=(0, zero))
    n = length(x)
    L = zeros(n, n)
    rhs = forcing.(x)
    for i in 2:n-1
        jleft = min(max(1, i-spoints÷2), n-spoints+1)
        js = jleft : jleft + spoints - 1
        L[i, js] = -fdstencil(x[js], x[i], 2)
    end
    L[1,1:spoints] = fdstencil(x[1:spoints], x[1], left[1])
    L[n,n-spoints+1:n] = fdstencil(x[n-spoints+1:n], x[n], right[1])
    rhs[1] = left[2](x[1])
    rhs[n] = right[2](x[n])
    L, rhs
end

CosRange(a, b, n) = (a + b)/2 .+ (b - a)/2 * cos.(LinRange(-pi, 0, n))

function vander_chebyshev(x, n=nothing)
    if isnothing(n)
        n = length(x) # Square by default
    end
    m = length(x)
    T = ones(m, n)
    if n > 1
        T[:, 2] = x
    end
    for k in 3:n
        T[:, k] = 2 * x .* T[:,k-1] - T[:, k-2]
    end
    T
end

function chebdiff(x, n=nothing)
    T = vander_chebyshev(x, n)
    m, n = size(T)
    dT = zero(T)
    dT[:,2:3] = [one.(x) 4*x]
    for j in 3:n-1
        dT[:,j+1] = j * (2 * T[:,j] + dT[:,j-1] / (j-2))
    end
    ddT = zero(T)
    ddT[:,3] .= 4
    for j in 3:n-1
        ddT[:,j+1] = j * (2 * dT[:,j] + ddT[:,j-1] / (j-2))
    end
    T, dT, ddT
end

function poisson_cheb(n, rhsfunc, leftbc=(0, zero), rightbc=(0, zero))
    x = CosRange(-1, 1, n)
    T, dT, ddT = chebdiff(x)
    dT /= T
    ddT /= T
    T /= T
    L = -ddT
    rhs = rhsfunc.(x)
    for (index, deriv, func) in
            [(1, leftbc...), (n, rightbc...)]
        L[index,:] = (T, dT)[deriv+1][index,:]
        rhs[index] = func(x[index])
    end
    x, L, rhs
end

┌ Info: Precompiling Plots [91a5bcdd-55d7-5caf-9e0b-520d859cae80]
└ @ Base loading.jl:1342

poisson_cheb (generic function with 3 methods)

Variable coefficients

• Heat conduction: steel, brick, wood, foam

• Electrical conductivity: copper, rubber, air

• Elasticity: steel, rubber, concrete, rock

• Linearization of nonlinear materials

  • ketchup, glacier ice, rocks (mantle/lithosphere)

kappa_step(x) = .1 + .9 * (x > 0)
kappa_tanh(x, epsilon=.1) = .55 + .45 * tanh(x / epsilon)
plot([kappa_step, kappa_tanh], xlims=(-1, 1), ylims=(0, 1), label="κ")

┌ Info: Precompiling GR_jll [d2c73de3-f751-5644-a686-071e5b155ba9]
└ @ Base loading.jl:1342

A naive finite difference solver

Conservative (divergence) form	Non-divergence form
\(-(\kappa u_x)_x = 0\)	\(-\kappa u_{xx} - \kappa_x u_x = 0\)

function poisson_nondivergence(x, spoints, kappa, forcing; leftbc=(0, zero), rightbc=(0, zero))
    n = length(x)
    L = zeros(n, n)
    rhs = forcing.(x)
    kappax = kappa.(x)
    for i in 2:n-1
        jleft = min(max(1, i-spoints÷2), n-spoints+1)
        js = jleft : jleft + spoints - 1
        kappa_x = fdstencil(x[js], x[i], 1) * kappax[js]
        L[i, js] = -fdstencil(x[js], x[i], 2) .* kappax[i] - fdstencil(x[js], x[i], 1) * kappa_x
    end
    L[1,1:spoints] = fdstencil(x[1:spoints], x[1], leftbc[1])
    if leftbc[1] == 1
        L[1, :] *= kappax[1]
    end
    L[n,n-spoints+1:n] = fdstencil(x[n-spoints+1:n], x[n], rightbc[1])
    if rightbc[1] == 1
        L[n, :] *= kappax[n]
    end
    rhs[1] = leftbc[2](x[1])
    rhs[n] = rightbc[2](x[n])
    L, rhs
end

poisson_nondivergence (generic function with 1 method)

Discretizing in conservative form

Conservative (divergence) form	Non-divergence form
\(-(\kappa u_x)_x = 0\)	\(-\kappa u_{xx} - \kappa_x u_x = 0\)

function poisson_conservative(n, kappa, forcing; leftbc=(0, zero), rightbc=(0, zero))
    x = LinRange(-1, 1, n)
    xstag = (x[1:end-1] + x[2:end]) / 2
    L = zeros(n, n)
    rhs = forcing.(x)
    kappa_stag = kappa.(xstag)
    for i in 2:n-1
        flux_L = kappa_stag[i-1] * fdstencil(x[i-1:i], xstag[i-1], 1)
        flux_R = kappa_stag[i] * fdstencil(x[i:i+1], xstag[i], 1)
        js = i-1:i+1
        weights = -fdstencil(xstag[i-1:i], x[i], 1)
        L[i, i-1:i+1] = weights[1] * [flux_L..., 0] + weights[2] * [0, flux_R...]
    end
    if leftbc[1] == 0
        L[1, 1] = 1
        rhs[1] = leftbc[2](x[1])
        rhs[2:end] -= L[2:end, 1] * rhs[1]
        L[2:end, 1] .= 0
    end
    if rightbc[1] == 0
        L[end,end] = 1
        rhs[end] = rightbc[2](x[end])
        rhs[1:end-1] -= L[1:end-1,end] * rhs[end]
        L[1:end-1,end] .= 0
    end
    x, L, rhs
end

poisson_conservative (generic function with 1 method)

Compare conservative vs non-divergence forms

forcing = zero # one
kappa_tanh(x, epsilon=.01) = .55 + .45 * tanh(x / epsilon)
x, L, rhs = poisson_conservative(20, kappa_tanh,
    forcing, leftbc=(0, zero), rightbc=(0, one))
u = L \ rhs
plot(x, u, marker=:auto, legend=:bottomright)

x = LinRange(-1, 1, 20)
L, rhs = poisson_nondivergence(x, 3, kappa_tanh,
    forcing, leftbc=(0, zero), rightbc=(0, one))
u = L \ rhs
plot(x, u, marker=:auto, legend=:bottomright)

Continuity of flux

forcing = zero
x, L, rhs = poisson_conservative(20, kappa_step,
    forcing, leftbc=(0, zero), rightbc=(0, one))
u = L \ rhs
plot(x, u, marker=:auto, legend=:bottomright)

xstag = (x[1:end-1] + x[2:end]) ./ 2
du = (u[1:end-1] - u[2:end]) ./ diff(x)
plot(xstag, [du .* kappa_step.(xstag)], marker=:auto, ylims=[-1, 1])

Manufactured solutions with discontinuous coefficients

• We need to be able to evaluate derivatives of the flux \(-\kappa u_x\).

• A physically-realizable solution would have continuous flux, but we we’d have to be making a physical solution to have that in verification.

• Idea: replace the discontinuous function with a continuous one with a rapid transition.

kappa_tanh(x, epsilon=.1) = .55 + .45 * tanh(x / epsilon)
d_kappa_tanh(x, epsilon=.1) = .45/epsilon * cosh(x/epsilon)^-2
plot([kappa_tanh])

Solving with the smoothed step \(\kappa\)

kappa_tanh(x, epsilon=.01) = .55 + .45 * tanh(x / epsilon)
d_kappa_tanh(x, epsilon=.01) = .45/epsilon * cosh(x/epsilon)^-2
manufactured(x) = tanh(2*x)
d_manufactured(x) = 2*cosh(2*x).^-2
flux_manufactured_kappa_tanh(x) = kappa_tanh(x) * d_manufactured(x)
function forcing_manufactured_kappa_tanh(x)
    8 * tanh(2x) / cosh(2x)^2 * kappa_tanh(x) -
    d_kappa_tanh(x) * d_manufactured(x)
end
x, L, rhs = poisson_conservative(200, kappa_tanh,
    forcing_manufactured_kappa_tanh,
    leftbc=(0, manufactured), rightbc=(0, manufactured))
u = L \ rhs
plot(x, u, marker=:auto, legend=:bottomright, title="Error $(norm(u - manufactured.(x), Inf))")
plot!([manufactured flux_manufactured_kappa_tanh forcing_manufactured_kappa_tanh kappa_tanh])

Convergence

function poisson_error(n, spoints=3)
    x, L, rhs = poisson_conservative(n, kappa_tanh,
        forcing_manufactured_kappa_tanh,
        leftbc=(0, manufactured), rightbc=(0, manufactured))
    u = L \ rhs
    norm(u - manufactured.(x), Inf)
end
ns = 2 .^ (3:10)
plot(ns, poisson_error.(ns, 3), marker=:auto, xscale=:log10, yscale=:log10)
plot!(n -> n^-2, label="\$1/n^2\$")

Advection

Advection represents transport in a “wind” \(w\). The time-dependent model is

\[ u_t + (w u)_x = 0 \]

for which \(u(t,x) = u(0, x - tw)\).

u0(x) = @. tanh(2*x)
usolution(t,x; w=1) = @. u0(x - t*w)
x = LinRange(-1, 1, 50)
plot(x, [usolution(t, x) for t in 0:.2:2])

What boundary conditions can be specified on the steady-state (\(t\to\infty\)) problem?

What is the steady state solution?

Advection-diffusion with Chebyshev

\[ -u_{xx} + (wu)_x = f \]

function advdiff_cheb(n, wind, rhsfunc, leftbc=(0, zero), rightbc=(0, zero))
    x = CosRange(-1, 1, n)
    T, dT, ddT = chebdiff(x)
    dT /= T
    ddT /= T
    T /= T
    L = -ddT + wind * dT
    rhs = rhsfunc.(x)
    for (index, deriv, func) in
            [(1, leftbc...), (n, rightbc...)]
        L[index,:] = (T, dT)[deriv+1][index,:]
        rhs[index] = func(x[index])
    end
    x, L, rhs
end

advdiff_cheb (generic function with 3 methods)

x, L, rhs = advdiff_cheb(25, 100, one)
@show minimum(diff(x))
plot(x, L \ rhs, legend=:none, marker=:auto)

minimum(diff(x)) = 0.008555138626189618

Advection-diffusion with conservative FD

\[(- \kappa u_x + wu)_x = f\]

function advdiff_conservative(n, kappa, wind, forcing)
    x = LinRange(-1, 1, n)
    xstag = (x[1:end-1] + x[2:end]) / 2
    L = zeros(n, n)
    rhs = forcing.(x)
    kappa_stag = kappa.(xstag)
    for i in 2:n-1
        flux_L = -kappa_stag[i-1] * fdstencil(x[i-1:i], xstag[i-1], 1) + wind * [.5 .5]
        flux_R = -kappa_stag[i] * fdstencil(x[i:i+1], xstag[i], 1) + wind * [.5 .5]
        js = i-1:i+1
        weights = fdstencil(xstag[i-1:i], x[i], 1)
        L[i, i-1:i+1] = weights[1] * [flux_L..., 0] + weights[2] * [0, flux_R...]
    end
    L[1, 1] = 1
    rhs[1] = 0
    L[2:end, 1] .= 0
    L[end,end] = 1
    rhs[end] = 0
    L[1:end-1,end] .= 0
    x, L, rhs
end

advdiff_conservative (generic function with 1 method)

Experiments, The (cell) Péclet number

(10)\[\begin{align} (- \kappa u_x + wu)_x &= f & \mathrm{Pe}_h &= \frac{h \lvert w \rvert}{\kappa} \end{align}\]

n = 30; h = 2/n
kappa = 1
wind = 10
x, L, rhs = advdiff_conservative(n, x -> kappa, wind, one)
@show minimum(diff(x))
plot(x, L \ rhs, marker=:auto, legend=:none, title="ratio $(wind*h/kappa)")

minimum(diff(x)) = 0.06896551724137923

Upwinded discretization

Idea: incoming advective flux should depend only on upwind value, outgoing should depend only on my value.

function advdiff_upwind(n, kappa, wind, forcing)
    x = LinRange(-1, 1, n)
    xstag = (x[1:end-1] + x[2:end]) / 2
    L = zeros(n, n)
    rhs = forcing.(x)
    kappa_stag = kappa.(xstag)
    for i in 2:n-1
        flux_L = -kappa_stag[i-1] * fdstencil(x[i-1:i], xstag[i-1], 1) +
            wind * (wind > 0 ? [1 0] : [0 1])
        flux_R = -kappa_stag[i] * fdstencil(x[i:i+1], xstag[i], 1) +
            wind * (wind > 0 ? [1 0] : [0 1])
        js = i-1:i+1
        weights = fdstencil(xstag[i-1:i], x[i], 1)
        L[i, i-1:i+1] = weights[1] * [flux_L..., 0] + weights[2] * [0, flux_R...]
    end
    L[1, 1] = 1
    rhs[1] = 0
    L[2:end, 1] .= 0
    L[end,end] = 1
    rhs[end] = 0
    L[1:end-1,end] .= 0
    x, L, rhs
end

advdiff_upwind (generic function with 1 method)

Try it for robustness

n = 30; h = 2/n
kappa = 1
wind = 10
x, L, rhs = advdiff_upwind(n, x -> kappa, wind, one)
@show minimum(diff(x))
plot(x, L \ rhs, marker=:auto, legend=:none, title="ratio $(wind*h/kappa)")

minimum(diff(x)) = 0.06896551724137923

previous

2021-09-13 Variable coefficients

next

2021-09-17 Sparsity, solvers, and nonlinear intro