{
"cells": [
{
"cell_type": "markdown",
"id": "3471b94b",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# 2022-09-09 Advection"
]
},
{
"cell_type": "markdown",
"id": "4159be0e",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"## Last time\n",
"\n",
"* Variable coefficients\n",
"* Conservative/divergence form vs non-divergence forms\n",
"* Verification with discontinuities\n",
"\n",
"## Today\n",
"\n",
"* Advection and boundary layers\n",
"* The (cell) PĂ©clet number and oscillations\n",
"* Godunov's Theorem"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "e82d841d",
"metadata": {
"hideOutput": true,
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [
{
"data": {
"text/plain": [
"poisson_conservative (generic function with 1 method)"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Plots\n",
"default(linewidth=3)\n",
"using LinearAlgebra\n",
"\n",
"function vander(x, k=nothing)\n",
" if k === nothing\n",
" k = length(x)\n",
" end\n",
" V = ones(length(x), k)\n",
" for j = 2:k\n",
" V[:, j] = V[:, j-1] .* x\n",
" end\n",
" V\n",
"end\n",
"\n",
"function fdstencil(source, target, k)\n",
" \"kth derivative stencil from source to target\"\n",
" x = source .- target\n",
" V = vander(x)\n",
" rhs = zero(x)'\n",
" rhs[k+1] = factorial(k)\n",
" rhs / V\n",
"end\n",
"\n",
"function my_spy(A)\n",
" cmax = norm(vec(A), Inf)\n",
" s = max(1, ceil(120 / size(A, 1)))\n",
" spy(A, marker=(:square, s), c=:diverging_rainbow_bgymr_45_85_c67_n256, clims=(-cmax, cmax))\n",
"end\n",
"\n",
"function poisson_nondivergence(n, kappa, forcing; leftbc=(0, zero), rightbc=(0, zero), spoints=3)\n",
" x = LinRange(-1, 1, n)\n",
" L = zeros(n, n)\n",
" rhs = forcing.(x)\n",
" kappax = kappa.(x)\n",
" for i in 2:n-1\n",
" jleft = min(max(1, i-spointsĂ·2), n-spoints+1)\n",
" js = jleft : jleft + spoints - 1\n",
" kappa_x = fdstencil(x[js], x[i], 1) * kappax[js]\n",
" L[i, js] = -fdstencil(x[js], x[i], 2) .* kappax[i] - fdstencil(x[js], x[i], 1) * kappa_x\n",
" end\n",
" L[1,1:spoints] = fdstencil(x[1:spoints], x[1], leftbc[1])\n",
" if leftbc[1] == 1\n",
" L[1, :] *= kappax[1]\n",
" end\n",
" L[n,n-spoints+1:n] = fdstencil(x[n-spoints+1:n], x[n], rightbc[1])\n",
" if rightbc[1] == 1\n",
" L[n, :] *= kappax[n]\n",
" end\n",
" rhs[1] = leftbc[2](x[1])\n",
" rhs[n] = rightbc[2](x[n])\n",
" x, L, rhs\n",
"end\n",
"\n",
"function poisson_conservative(n, kappa, forcing; leftbc=(0, zero), rightbc=(0, zero))\n",
" x = LinRange(-1, 1, n)\n",
" xstag = (x[1:end-1] + x[2:end]) / 2\n",
" L = zeros(n, n)\n",
" rhs = forcing.(x)\n",
" kappa_stag = kappa.(xstag)\n",
" for i in 2:n-1\n",
" flux_L = kappa_stag[i-1] * fdstencil(x[i-1:i], xstag[i-1], 1)\n",
" flux_R = kappa_stag[i] * fdstencil(x[i:i+1], xstag[i], 1)\n",
" js = i-1:i+1\n",
" weights = -fdstencil(xstag[i-1:i], x[i], 1)\n",
" L[i, i-1:i+1] = weights[1] * [flux_L..., 0] + weights[2] * [0, flux_R...]\n",
" end\n",
" if leftbc[1] == 0\n",
" L[1, 1] = 1\n",
" rhs[1] = leftbc[2](x[1])\n",
" rhs[2:end] -= L[2:end, 1] * rhs[1]\n",
" L[2:end, 1] .= 0\n",
" end\n",
" if rightbc[1] == 0\n",
" L[end,end] = 1\n",
" rhs[end] = rightbc[2](x[end])\n",
" rhs[1:end-1] -= L[1:end-1,end] * rhs[end]\n",
" L[1:end-1,end] .= 0\n",
" end\n",
" x, L, rhs\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "443c90e4",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Solving with the smoothed step $\\kappa$"
]
},
{
"cell_type": "code",
"execution_count": 69,
"id": "f6710a06",
"metadata": {
"cell_style": "split"
},
"outputs": [],
"source": [
"manufactured(x) = tanh(2x)\n",
"d_manufactured(x) = 2*cosh(2x)^-2\n",
"mdd_manufactured(x) = 8 * tanh(2x) / cosh(2x)^2\n",
"kappa_tanh(x, epsilon=.02)= .55 + .45 * tanh(x / epsilon)\n",
"d_kappa_tanh(x, epsilon=.02) = .45/epsilon * cosh(x/epsilon)^-2\n",
"flux_manufactured_kappa_tanh(x) = kappa_tanh(x) * d_manufactured(x)\n",
"function forcing_manufactured_kappa_tanh(x)\n",
" 8 * tanh(2x) / cosh(2x)^2 * kappa_tanh(x) -\n",
" d_kappa_tanh(x) * d_manufactured(x)\n",
"end\n",
"x, L, rhs = poisson_conservative(40, kappa_tanh,\n",
" forcing_manufactured_kappa_tanh,\n",
" leftbc=(0, manufactured), rightbc=(0, manufactured));"
]
},
{
"cell_type": "code",
"execution_count": 70,
"id": "b6f536fe",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 70,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"u = L \\ rhs\n",
"plot(x, u, marker=:circle, legend=:bottomright, title=\"Error $(norm(u - manufactured.(x), Inf))\")\n",
"plot!([manufactured flux_manufactured_kappa_tanh forcing_manufactured_kappa_tanh kappa_tanh],\n",
"label=[\"u_manufactured\" \"flux\" \"forcing\" \"kappa\"], ylim=(-5, 5))"
]
},
{
"cell_type": "markdown",
"id": "7b533f94",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# MMS Convergence: nondivergence vs conservative"
]
},
{
"cell_type": "code",
"execution_count": 72,
"id": "fc5d08d6",
"metadata": {
"slideshow": {
"slide_type": ""
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 72,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function poisson_error(n, spoints=3)\n",
" x, L, rhs = poisson_conservative(n, kappa_tanh,\n",
" forcing_manufactured_kappa_tanh,\n",
" leftbc=(0, manufactured), rightbc=(0, manufactured))\n",
" u = L \\ rhs\n",
" norm(u - manufactured.(x), Inf)\n",
"end\n",
"ns = 2 .^ (3:10)\n",
"plot(ns, poisson_error.(ns, 3), marker=:auto, xscale=:log10, yscale=:log10)\n",
"plot!(n -> n^-2, label=\"\\$1/n^2\\$\")"
]
},
{
"cell_type": "markdown",
"id": "7f09d92d",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Outlook\n",
"\n",
"* Manufactured solutions can be misleading for verification in the presence of rapid coefficient variation.\n",
" * Forcing terms may be poorly resolved on the grid\n",
"* Non-conservative methods are scary\n",
" * Can look great during verification on smooth problems, then break catastrophically\n",
"* Conservative methods are generally robust, but now we have choices to make\n",
" * What if $\\kappa(x)$ is gridded data? How do we evaluate at staggered points?\n",
" * What does this mean for boundary conditions, $\\kappa u_x = g$\n",
" * What do staggered points mean on unstructured/adaptive meshes?\n",
" * Which terms should we evaluate where?"
]
},
{
"cell_type": "markdown",
"id": "f2b911c6",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Advection\n",
"\n",
"Advection represents transport in a \"wind\" $w$. The time-dependent model is\n",
"$$ u_t + (w u)_x = 0 $$\n",
"for which $u(t,x) = u(0, x - tw)$."
]
},
{
"cell_type": "code",
"execution_count": 32,
"id": "2b09c440",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"u0(x) = @. tanh(2*x)\n",
"usolution(t,x; w=1) = @. u0(x - t*w)\n",
"x = LinRange(-1, 1, 50)\n",
"plot(x, [usolution(t, x) for t in 0:.2:2], legend=:none)"
]
},
{
"cell_type": "markdown",
"id": "c7a59c8d",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"## What boundary conditions can be specified on the steady-state ($t\\to\\infty$) problem?\n",
"\n",
"## What is the steady state solution?"
]
},
{
"cell_type": "markdown",
"id": "7233ddb3",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Advection-diffusion with Chebyshev\n",
"$$ \\underbrace{u_t}_{\\text{ignore today}} + \\underbrace{(wu)_x}_{\\text{advection}} - \\underbrace{u_{xx}}_{\\text{diffusion}} = \\underbrace{f}_{\\text{source}} $$"
]
},
{
"cell_type": "code",
"execution_count": 73,
"id": "bc035ef0",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"advdiff_cheb (generic function with 3 methods)"
]
},
"execution_count": 73,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function advdiff_cheb(n, wind, rhsfunc, leftbc=(0, zero), rightbc=(0, zero))\n",
" x = CosRange(-1, 1, n)\n",
" T, dT, ddT = chebdiff(x)\n",
" dT /= T\n",
" ddT /= T\n",
" T /= T\n",
" L = -ddT + wind * dT\n",
" rhs = rhsfunc.(x)\n",
" for (index, deriv, func) in\n",
" [(1, leftbc...), (n, rightbc...)]\n",
" L[index,:] = (T, dT)[deriv+1][index,:]\n",
" rhs[index] = func(x[index])\n",
" end\n",
" x, L, rhs\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 81,
"id": "0c4705b4",
"metadata": {
"cell_style": "split",
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"minimum(diff(x)) = 0.0032426918657900705\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 81,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"wind = 100\n",
"x, L, rhs = advdiff_cheb(40, wind, one)\n",
"@show minimum(diff(x))\n",
"plot(x, L \\ rhs, legend=:none, marker=:auto)"
]
},
{
"cell_type": "markdown",
"id": "b76caa3d",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Advection-diffusion with conservative FD\n",
"\n",
"$$\\big(wu - \\kappa u_x\\big)_x = f$$"
]
},
{
"cell_type": "code",
"execution_count": 82,
"id": "27f723a2",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": ""
}
},
"outputs": [
{
"data": {
"text/plain": [
"advdiff_conservative (generic function with 1 method)"
]
},
"execution_count": 82,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function advdiff_conservative(n, kappa, wind, forcing)\n",
" x = LinRange(-1, 1, n)\n",
" xstag = (x[1:end-1] + x[2:end]) / 2\n",
" L = zeros(n, n)\n",
" rhs = forcing.(x)\n",
" kappa_stag = kappa.(xstag)\n",
" for i in 2:n-1\n",
" flux_L = -kappa_stag[i-1] * fdstencil(x[i-1:i], xstag[i-1], 1) + wind * [.5 .5]\n",
" flux_R = -kappa_stag[i] * fdstencil(x[i:i+1], xstag[i], 1) + wind * [.5 .5]\n",
" js = i-1:i+1\n",
" weights = fdstencil(xstag[i-1:i], x[i], 1)\n",
" L[i, i-1:i+1] = weights[1] * [flux_L..., 0] + weights[2] * [0, flux_R...]\n",
" end\n",
" L[1, 1] = 1\n",
" rhs[1] = 0\n",
" L[2:end, 1] .= 0\n",
" L[end,end] = 1\n",
" rhs[end] = 0\n",
" L[1:end-1,end] .= 0\n",
" x, L, rhs\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "359f81d7",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Experiments, The (cell) [PĂ©clet number](https://en.wikipedia.org/wiki/P%C3%A9clet_number)\n",
"\n",
"\\begin{align}\n",
"\\big(wu - \\kappa u_x\\big)_x &= f & \\mathrm{Pe}_h &= \\frac{h \\lvert w \\rvert}{\\kappa}\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": 98,
"id": "4cb9b977",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"minimum(diff(x)) = 0.06896551724137923\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 98,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = 30; h = 2/n\n",
"kappa = 1\n",
"wind = 20\n",
"x, L, rhs = advdiff_conservative(n, x -> kappa, wind, one)\n",
"@show minimum(diff(x))\n",
"plot(x, L \\ rhs, marker=:auto, legend=:none, title=\"\\$Pe_h = $(wind*h/kappa)\\$\")"
]
},
{
"cell_type": "markdown",
"id": "862a1396",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Upwinded discretization\n",
"\n",
"Idea: incoming **advective flux** should depend only on upwind value, outgoing should depend only on my value."
]
},
{
"cell_type": "code",
"execution_count": 99,
"id": "12cc82ae",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"advdiff_upwind (generic function with 1 method)"
]
},
"execution_count": 99,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function advdiff_upwind(n, kappa, wind, forcing)\n",
" x = LinRange(-1, 1, n)\n",
" xstag = (x[1:end-1] + x[2:end]) / 2\n",
" L = zeros(n, n)\n",
" rhs = forcing.(x)\n",
" kappa_stag = kappa.(xstag)\n",
" for i in 2:n-1\n",
" flux_L = -kappa_stag[i-1] * fdstencil(x[i-1:i], xstag[i-1], 1) +\n",
" wind * (wind > 0 ? [1 0] : [0 1])\n",
" flux_R = -kappa_stag[i] * fdstencil(x[i:i+1], xstag[i], 1) +\n",
" wind * (wind > 0 ? [1 0] : [0 1])\n",
" js = i-1:i+1\n",
" weights = fdstencil(xstag[i-1:i], x[i], 1)\n",
" L[i, i-1:i+1] = weights[1] * [flux_L..., 0] + weights[2] * [0, flux_R...]\n",
" end\n",
" L[1, 1] = 1\n",
" rhs[1] = 0\n",
" L[2:end, 1] .= 0\n",
" L[end,end] = 1\n",
" rhs[end] = 0\n",
" L[1:end-1,end] .= 0\n",
" x, L, rhs\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "8ea1c290",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Try it for robustness"
]
},
{
"cell_type": "code",
"execution_count": 103,
"id": "a1b82c78",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"minimum(diff(x)) = 0.06896551724137923\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 103,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = 30; h = 2/n\n",
"kappa = 1\n",
"wind = 100000\n",
"x, L, rhs = advdiff_upwind(n, x -> kappa, wind, one)\n",
"@show minimum(diff(x))\n",
"plot(x, L \\ rhs, marker=:auto, legend=:none, title=\"\\$Pe_h = $(wind*h/kappa)\\$\")"
]
},
{
"cell_type": "markdown",
"id": "99d7b03f",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"## What is the order of accuracy of the upwinded discretization?"
]
},
{
"cell_type": "markdown",
"id": "389c0694",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Godunov's Theorem (1954)"
]
},
{
"cell_type": "markdown",
"id": "fae87e3c",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"Linear numerical methods\n",
"$$ \\dot u_i = \\sum_j a_{ij} u_j $$\n",
"for solving time-dependent PDE can be at most first order accurate if they are monotone.\n",
"\n",
"For our purposes, monotonicity is equivalent to positivity preservation,\n",
"$$ \\min_x u(0, x) \\ge 0 \\quad \\Longrightarrow \\quad \\min_x u(t, 0) \\ge 0 .$$"
]
},
{
"cell_type": "markdown",
"id": "01364912",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"## Discontinuities\n",
"\n",
"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,\n",
"$$ \\lVert u - u^* \\rVert_{L^1} = \\int \\lvert u - u^* \\rvert . $$\n",
"If we merely sample a discontinuous function, say\n",
"$$ u(x) = \\begin{cases} 0, & x \\le a \\\\ 1, & x > a \\end{cases} $$\n",
"onto a grid with element size $\\Delta x$ then we will have errors of order 1 on an interval proportional to $\\Delta x$."
]
},
{
"cell_type": "markdown",
"id": "4608bb6f",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": ""
}
},
"source": [
"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."
]
},
{
"cell_type": "markdown",
"id": "d1597f41",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Next week\n",
"\n",
"* Sparse linear algebra (necessary for more than 1 spatial dimension)\n",
"* Nonlinear solvers\n",
"* Time integration for stiff and non-stiff systems"
]
}
],
"metadata": {
"@webio": {
"lastCommId": null,
"lastKernelId": null
},
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Julia 1.7.2",
"language": "julia",
"name": "julia-1.7"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.8.1"
},
"rise": {
"enable_chalkboard": true
}
},
"nbformat": 4,
"nbformat_minor": 5
}