{
"cells": [
{
"cell_type": "markdown",
"id": "3dba1684",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# 2021-12-01 Finite Volume methods\n",
"\n",
"## Last time\n",
"\n",
"* Notes on unstructured meshing workflow\n",
"* Finite volume methods for hyperbolic conservation laws\n",
"* Riemann solvers for scalar equations\n",
" * Shocks and the Rankine-Hugoniot condition\n",
" * Rarefactions and entropy solutions\n",
"\n",
"## Today\n",
"\n",
"* Higher order methods\n",
"* Godunov's Theorem\n",
"* Slope reconstruction\n",
"* Slope limiting"
]
},
{
"cell_type": "code",
"execution_count": 84,
"id": "7fd7f3b0",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [
{
"data": {
"text/plain": [
"riemann_traffic (generic function with 1 method)"
]
},
"execution_count": 84,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using LinearAlgebra\n",
"using Plots\n",
"default(linewidth=4)\n",
"\n",
"struct RKTable\n",
" A::Matrix\n",
" b::Vector\n",
" c::Vector\n",
" function RKTable(A, b)\n",
" s = length(b)\n",
" A = reshape(A, s, s)\n",
" c = vec(sum(A, dims=2))\n",
" new(A, b, c)\n",
" end\n",
"end\n",
"\n",
"rk4 = RKTable([0 0 0 0; .5 0 0 0; 0 .5 0 0; 0 0 1 0], [1, 2, 2, 1] / 6)\n",
"\n",
"function ode_rk_explicit(f, u0; tfinal=1., h=0.1, table=rk4)\n",
" u = copy(u0)\n",
" t = 0.\n",
" n, s = length(u), length(table.c)\n",
" fY = zeros(n, s)\n",
" thist = [t]\n",
" uhist = [u0]\n",
" while t < tfinal\n",
" tnext = min(t+h, tfinal)\n",
" h = tnext - t\n",
" for i in 1:s\n",
" ti = t + h * table.c[i]\n",
" Yi = u + h * sum(fY[:,1:i-1] * table.A[i,1:i-1], dims=2)\n",
" fY[:,i] = f(ti, Yi)\n",
" end\n",
" u += h * fY * table.b\n",
" t = tnext\n",
" push!(thist, t)\n",
" push!(uhist, u)\n",
" end\n",
" thist, hcat(uhist...)\n",
"end\n",
"\n",
"function testfunc(x)\n",
" max(1 - 4*abs.(x+2/3),\n",
" abs.(x) .< .2,\n",
" (2*abs.(x-2/3) .< .5) * cospi(2*(x-2/3)).^2\n",
" )\n",
"end\n",
"\n",
"riemann_advection(uL, uR) = 1*uL # velocity is +1\n",
"\n",
"function fv_solve1(riemann, u_init, n, tfinal=1)\n",
" h = 2 / n\n",
" x = LinRange(-1+h/2, 1-h/2, n) # cell midpoints (centroids)\n",
" idxL = 1 .+ (n-1:2*n-2) .% n\n",
" idxR = 1 .+ (n+1:2*n) .% n\n",
" function rhs(t, u)\n",
" fluxL = riemann(u[idxL], u)\n",
" fluxR = riemann(u, u[idxR])\n",
" (fluxL - fluxR) / h\n",
" end\n",
" thist, uhist = ode_rk_explicit(\n",
" rhs, u_init.(x), h=h, tfinal=tfinal)\n",
" x, thist, uhist\n",
"end\n",
"\n",
"function riemann_burgers(uL, uR)\n",
" flux = zero(uL)\n",
" for i in 1:length(flux)\n",
" fL = flux_burgers(uL[i])\n",
" fR = flux_burgers(uR[i])\n",
" flux[i] = if uL[i] > uR[i] # shock\n",
" max(fL, fR)\n",
" elseif uL[i] > 0 # rarefaction all to the right\n",
" fL\n",
" elseif uR[i] < 0 # rarefaction all to the left\n",
" fR\n",
" else\n",
" 0\n",
" end\n",
" end\n",
" flux\n",
"end\n",
"\n",
"function riemann_traffic(uL, uR)\n",
" flux = zero(uL)\n",
" for i in 1:length(flux)\n",
" fL = flux_traffic(uL[i])\n",
" fR = flux_traffic(uR[i])\n",
" flux[i] = if uL[i] < uR[i] # shock\n",
" min(fL, fR)\n",
" elseif uL[i] < .5 # rarefaction all to the right\n",
" fL\n",
" elseif uR[i] > .5 # rarefaction all to the left\n",
" fR\n",
" else\n",
" flux_traffic(.5)\n",
" end\n",
" end\n",
" flux\n",
"end"
]
},
{
"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": "3ef3d115",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Slope Reconstruction"
]
},
{
"cell_type": "markdown",
"id": "17d89dc5",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"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\n",
"$$ \\tilde u_i(x) = u_i + g_i \\cdot (x - x_i) . $$\n",
"We would like this reconstruction to be monotone in the sense that\n",
"$$ \\Big(\\tilde u_i(x) - \\tilde u_j(x) \\Big) \\Big( u_i - u_j \\Big) \\ge 0 $$\n",
"for any $x$ on the interface between element $i$ and element $j$."
]
},
{
"cell_type": "markdown",
"id": "a7121dc1",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"## Question\n",
"Is the symmetric slope\n",
"$$ \\hat g_i = \\frac{u_{i+1} - u_{i-1}}{2 \\Delta x} $$\n",
"monotone?"
]
},
{
"cell_type": "markdown",
"id": "2a75e74c",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Slope limiting\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "ba2741ca",
"metadata": {
"cell_style": "split"
},
"source": [
"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)](https://www.nas.nasa.gov/assets/pdf/techreports/2005/nas-05-007.pdf).\n",
"\n",
"![Berger et al, figure 3b, Cell averages](../img/Berger2005-3b-CellAverages.png)"
]
},
{
"cell_type": "markdown",
"id": "e0fef82e",
"metadata": {
"cell_style": "split"
},
"source": [
"We will express a slope limiter in terms of the ratio\n",
"$$ r_i = \\frac{u_i - u_{i-1}}{u_{i+1} - u_{i-1}} $$\n",
"and use as a gradient,\n",
"$$ g_i = \\phi(r_i) \\frac{u_{i+1} - u_{i-1}}{2 \\Delta x} . $$\n",
"Our functions $\\phi$ will be zero unless $0 < r < 1$ and $\\phi(1/2) = 1$.\n",
"![Berger et al, figure 3a, Symmetric form of slope limiters](../img/Berger2005-3a-SymSlope.png)\n",
"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."
]
},
{
"cell_type": "markdown",
"id": "a140c8da",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Common limiters"
]
},
{
"cell_type": "code",
"execution_count": 46,
"id": "c4171796",
"metadata": {
"cell_style": "split"
},
"outputs": [],
"source": [
"limit_zero(r) = 0\n",
"limit_none(r) = 1\n",
"limit_minmod(r) = max(min(2*r, 2*(1-r)), 0)\n",
"limit_sin(r) = (0 < r && r < 1) * sinpi(r)\n",
"limit_vl(r) = max(4*r*(1-r), 0)\n",
"limit_bj(r) = max(0, min(1, 4*r, 4*(1-r)))\n",
"limiters = [limit_zero limit_none limit_minmod limit_sin limit_vl limit_bj];"
]
},
{
"cell_type": "code",
"execution_count": 45,
"id": "8c74fd71",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(limiters, label=limiters, xlims=(-.1, 1.1))"
]
},
{
"cell_type": "markdown",
"id": "70c37c4b",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# A slope-limited solver"
]
},
{
"cell_type": "code",
"execution_count": 78,
"id": "a0a6c118",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"fv_solve2 (generic function with 3 methods)"
]
},
"execution_count": 78,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function fv_solve2(riemann, u_init, n, tfinal=1, limit=limit_sin)\n",
" h = 2 / n\n",
" x = LinRange(-1+h/2, 1-h/2, n) # cell midpoints (centroids)\n",
" idxL = 1 .+ (n-1:2*n-2) .% n\n",
" idxR = 1 .+ (n+1:2*n) .% n\n",
" function rhs(t, u)\n",
" jump = u[idxR] - u[idxL]\n",
" r = (u - u[idxL]) ./ jump\n",
" r[isnan.(r)] .= 0\n",
" g = limit.(r) .* jump / 2h\n",
" fluxL = riemann(u[idxL] + g[idxL]*h/2, u - g*h/2)\n",
" fluxR = fluxL[idxR]\n",
" (fluxL - fluxR) / h\n",
" end\n",
" thist, uhist = ode_rk_explicit(\n",
" rhs, u_init.(x), h=h, tfinal=tfinal)\n",
" x, thist, uhist\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 83,
"id": "6bec13ec",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 83,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x, thist, uhist = fv_solve2(riemann_advection, testfunc, 100, .5,\n",
" limit_sin)\n",
"plot(x, uhist[:,1:10:end], legend=:none)"
]
}
],
"metadata": {
"celltoolbar": "Slideshow",
"kernelspec": {
"display_name": "Julia 1.6.2",
"language": "julia",
"name": "julia-1.6"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.6.4"
},
"rise": {
"enable_chalkboard": true
}
},
"nbformat": 4,
"nbformat_minor": 5
}