## Last time#

• Variable coefficients

• Conservative/divergence form vs non-divergence forms

• Verification with discontinuities

## Today#

• The (cell) Péclet number and oscillations

• Godunov’s Theorem

using Plots
default(linewidth=3)
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 my_spy(A)
cmax = norm(vec(A), Inf)
s = max(1, ceil(120 / size(A, 1)))
spy(A, marker=(:square, s), c=:diverging_rainbow_bgymr_45_85_c67_n256, clims=(-cmax, cmax))
end

function poisson_nondivergence(n, kappa, forcing; leftbc=(0, zero), rightbc=(0, zero), spoints=3)
x = LinRange(-1, 1, n)
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])
x, L, rhs
end

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

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 my_spy(A)
cmax = norm(vec(A), Inf)
s = max(1, ceil(120 / size(A, 1)))
spy(A, marker=(:square, s), c=:diverging_rainbow_bgymr_45_85_c67_n256, clims=(-cmax, cmax))
end

my_spy (generic function with 1 method)


# Solving with the smoothed step $$\kappa$$#

manufactured(x) = tanh(2x)
d_manufactured(x) = 2*cosh(2x)^-2
mdd_manufactured(x) = 8 * tanh(2x) / cosh(2x)^2
kappa_tanh(x, epsilon=.02)= .55 + .45 * tanh(x / epsilon)
d_kappa_tanh(x, epsilon=.02) = .45/epsilon * cosh(x/epsilon)^-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(30, kappa_tanh,
forcing_manufactured_kappa_tanh,
leftbc=(0, manufactured), rightbc=(0, manufactured));

u = L \ rhs
plot(x, u, marker=:circle, legend=:bottomright, title="Error $(norm(u - manufactured.(x), Inf))") plot!([manufactured flux_manufactured_kappa_tanh forcing_manufactured_kappa_tanh kappa_tanh], label=["u_manufactured" "flux" "forcing" "kappa"], ylim=(-3, 3))  # MMS Convergence: nondivergence vs conservative# 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\")  # Outlook# • Manufactured solutions can be misleading for verification in the presence of rapid coefficient variation. • Forcing terms may be poorly resolved on the grid • Non-conservative methods are scary • Can look great during verification on smooth problems, then break catastrophically • Conservative methods are generally robust, but now we have choices to make • What if $$\kappa(x)$$ is gridded data? How do we evaluate at staggered points? • What does this mean for boundary conditions, $$\kappa u_x = g$$ • What do staggered points mean on unstructured/adaptive meshes? • Which terms should we evaluate where? # 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], legend=:none)  ## What boundary conditions can be specified on the steady-state ($$t\to\infty$$) problem?# ## What is the steady state solution?# # Advection-diffusion with Chebyshev# $\underbrace{u_t}_{\text{ignore today}} + \underbrace{(wu)_x}_{\text{advection}} - \underbrace{u_{xx}}_{\text{diffusion}} = \underbrace{f}_{\text{source}}$ 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)  wind = 200 x, L, rhs = advdiff_cheb(40, wind, one) @show minimum(diff(x)) plot(x, L \ rhs, legend=:none, marker=:auto)  minimum(diff(x)) = 0.0032426918657900705  # Advection-diffusion with conservative FD# $\big(wu - \kappa u_x\big)_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} \big(wu - \kappa u_x\big)_x &= f & \mathrm{Pe}_h &= \frac{h \lvert w \rvert}{\kappa} \end{align} n = 60; h = 2/n kappa = 1 wind = 8000 x, L, rhs = advdiff_conservative(n, x -> kappa, wind, one) @show minimum(diff(x)) plot(x, L \ rhs, marker=:auto, legend=:none, title="\Pe_h = $(wind*h/kappa)\$")

minimum(diff(x)) = 0.03389830508474567


# 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 = 100000
x, L, rhs = advdiff_upwind(n, x -> kappa, wind, one)
@show minimum(diff(x))
plot(x, L \ rhs, marker=:auto, legend=:none, title="\$Pe_h =$(wind*h/kappa)\\$")

minimum(diff(x)) = 0.06896551724137923


# 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, x) \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.

# Next week#

• Sparse linear algebra (necessary for more than 1 spatial dimension)

• Nonlinear solvers

• Time integration for stiff and non-stiff systems