{
"cells": [
{
"cell_type": "markdown",
"id": "3471b94b",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# 2022-11-30 Reconstruction"
]
},
{
"cell_type": "markdown",
"id": "4159be0e",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"## 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",
"* Godunov's Theorem\n",
"* Slope reconstruction and limiting\n",
"* Hyperbolic systems\n",
" * Rankine-Hugoniot and Riemann Invariants\n",
" * Exact Riemann solvers\n",
" * Approximate Riemann solvers"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "7fd7f3b0",
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [
{
"data": {
"text/plain": [
"riemann_traffic (generic function with 1 method)"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using LinearAlgebra\n",
"using Plots\n",
"default(linewidth=3)\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",
"flux_advection(u) = u\n",
"flux_burgers(u) = u^2/2\n",
"flux_traffic(u) = u * (1 - u)\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(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": "d6067fcd",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Godunov methods (first order accurate)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "ff0c7b0d",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"init_func(x) = testfunc(x)\n",
"x, thist, uhist = fv_solve1(riemann_traffic, init_func, 100, 2)\n",
"plot(x, uhist[:,1:5:end], legend=:none)"
]
},
{
"cell_type": "markdown",
"id": "2184d13f",
"metadata": {
"cell_style": "split"
},
"source": [
"## Burgers\n",
"\n",
"* flux $u^2/2$ has speed $u$\n",
"* negative values make sense\n",
"* satisfies a maximum principle\n",
"\n",
"## Traffic\n",
"\n",
"* flux $u - u^2$ has speed $1 - 2u$\n",
"* state must satisfy $u \\in [0, 1]$"
]
},
{
"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": 3,
"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": 4,
"id": "8c74fd71",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 4,
"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": 14,
"id": "a0a6c118",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"fv_solve2 (generic function with 3 methods)"
]
},
"execution_count": 14,
"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": 20,
"id": "6bec13ec",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 20,
"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)"
]
},
{
"cell_type": "markdown",
"id": "25d011b3",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Hyperbolic systems\n",
"\n",
"$$ U_t + f(U)_x = 0 $$"
]
},
{
"cell_type": "markdown",
"id": "598defce",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"## Isentropic gas dynamics\n",
"\n",
"\\begin{align} U &= \\begin{bmatrix} \\rho \\\\ \\rho u \\end{bmatrix} & f(U) &= \\begin{bmatrix} \\rho u \\\\ \\rho u^2 + p \\end{bmatrix} \\end{align}\n",
"\n",
"| Variable | meaning |\n",
"|----------|--------|\n",
"| $\\rho$ | density |\n",
"| $u$ | velocity |\n",
"| $\\rho u$ | momentum |\n",
"| $p$ | pressure |\n",
"\n",
"\n",
"* Equation of state $p(\\rho) = C \\rho^\\gamma$ with $\\gamma = 1.4$ (typical air). \n",
"* \"isothermal\" gas dynamics: $p(\\rho) = c^2 \\rho$, wave speed $c$.\n",
"* Compute as $ \\rho u^2 = \\frac{(\\rho u)^2}{\\rho} .$"
]
},
{
"cell_type": "markdown",
"id": "35a1e510",
"metadata": {
"cell_style": "split"
},
"source": [
"## Smooth wave structure\n",
"\n",
"For perturbations of a constant state, systems of equations produce multiple waves with speeds equal to the eigenvalues of the flux Jacobian,\n",
"$$ f'(U) = \\begin{bmatrix} 0 & 1 \\\\ -u^2 + c^2 & 2 u \\end{bmatrix}. $$\n",
"We can compute the eigenvalues and eigenvectors,\n",
"$$ f'(U) = W \\Lambda W^{-1} $$\n",
"as\n",
"\\begin{align} \n",
"W &= \\begin{bmatrix} 1 & 1 \\\\ u-c & u+c \\end{bmatrix} &\n",
"\\Lambda &= \\begin{bmatrix} u-c & \\\\ & u+c \\end{bmatrix} .\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"id": "1d945225",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Riemann problem for systems: shocks\n",
"\n",
"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.\n",
"The type of wave will determine the condition that must be satisfied to connect $U_L$ to $U_*$ and $U_*$ to $U_R$."
]
},
{
"cell_type": "markdown",
"id": "15492bf5",
"metadata": {
"cell_style": "split"
},
"source": [
"## Left-moving wave\n",
"\n",
"If there is a shock between $U_L$ and $U_*$, the Rankine-Hugoniot condition\n",
"$$ s \\Delta U = \\Delta f $$\n",
"will be satisfied along with the entropy condition\n",
"$$ \\lambda_i(U_L) \\ge s \\ge \\lambda_i(U_*) $$\n",
"where $\\lambda_i(U)$ is the corresponding eigenvalue of the flux Jacobian $f'(U)$.\n",
"The inequality in the entropy condition is strict if the wave is genuinely nonlinear.\n",
"For isothermal gas dynamics, the Rankine-Hugoniot condition is\n",
"\\begin{align}\n",
" (\\rho_* - \\rho_L) s &= \\rho_* u_* - \\rho_L u_L \\\\\n",
" (\\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)\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"id": "927703d5",
"metadata": {
"cell_style": "split"
},
"source": [
"Solving the first equation for $s$ and substituting into the second, we compute\n",
"\\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) \\\\\n",
" \\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 \\\\\n",
" \\rho_* \\rho_L \\Big( u_*^2 - 2 u_* u_L + u_L^2 \\Big) &= c^2 (\\rho_* - \\rho_L)^2 \\\\\n",
" (u_* - u_L)^2 &= c^2 \\frac{(\\rho_* - \\rho_L)^2}{\\rho_* \\rho_L} \\\\\n",
" u_* - u_L &= \\pm c \\frac{\\rho_* - \\rho_L}{\\sqrt{\\rho_* \\rho_L}}\n",
"\\end{split}\n",
"and will need to use the entropy condition to learn which sign to take."
]
},
{
"cell_type": "markdown",
"id": "56dbb58e",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Admissible shocks"
]
},
{
"cell_type": "markdown",
"id": "7df0f6a7",
"metadata": {
"cell_style": "split"
},
"source": [
"We need to choose the sign\n",
"$$u_* - u_L = \\pm c \\frac{\\rho_* - \\rho_L}{\\sqrt{\\rho_*\\rho_L}}$$"
]
},
{
"cell_type": "markdown",
"id": "1591c7d4",
"metadata": {
"cell_style": "split"
},
"source": [
"The entropy condition requires that\n",
"$$ \\lambda(U_L) = u_L - c \\ge s \\ge u_* - c = \\lambda(U_*) . $$\n",
"Meanwhile, the first Rankine-Hugoniot condition and $\\rho > 0$ provides\n",
"$$ \\rho_L c \\le \\rho_L (u_L - s) = \\rho_* (u_* - s) \\le \\rho_* c $$\n",
"which is to say that $\\rho_* \\ge \\rho_L$ and consequently we take the negative sign,\n",
"$$ u_* - u_L = - c\\frac{\\rho_* - \\rho_L}{\\sqrt{\\rho_* \\rho_L}} .$$\n",
"## Right-moving wave\n",
"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\n",
"$$ u_R - u_* = - c \\frac{\\rho_* - \\rho_R}{\\sqrt{\\rho_* \\rho_R}} . $$"
]
},
{
"cell_type": "markdown",
"id": "82495d67",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Rarefactions\n",
"\n",
"A rarefaction occurs when\n",
"$$ \\lambda_i(U_L) < \\lambda_i(U_R) $$"
]
},
{
"cell_type": "markdown",
"id": "4417d692",
"metadata": {
"cell_style": "split"
},
"source": [
"## Generalized Riemann invariants\n",
"\n",
"$$ \\frac{d U^1}{W_i^1} = \\frac{d U^2}{W_i^2} = \\dotsb $$\n",
"\n",
"(Derivation of this condition is beyond the scope of this class.)\n",
"\n",
"### Isothermal gas dynamics\n",
"$$ \\frac{d \\rho}{1} = \\frac{d (\\rho u)}{u - c} $$\n",
"across the wave $\\lambda = u-c$.\n",
"We can rearrange to\n",
"$$ (u-c) d \\rho = \\rho d u + u d \\rho $$\n",
"or\n",
"$$ d u + \\frac{c}{\\rho} d \\rho = 0 . $$"
]
},
{
"cell_type": "markdown",
"id": "71a05127",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"Integration yields\n",
"$$ u + \\int \\frac{c}{\\rho} d \\rho = u + c \\log \\rho = \\text{constant} . $$\n",
"So if we have a rarefaction between $U_L$ and $U_*$ then\n",
"$$ u_L + c \\log \\rho_L = u_* + c \\log \\rho_* . $$\n",
"Similarly, for the wave $\\lambda = u + c$, a rarefaction between $U_*$ and $U_R$ results in\n",
"$$ u_* - c \\log \\rho_* = u_R - c \\log \\rho_R . $$\n",
"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\n",
"$$ u_* -c \\log \\rho_* = u_0 - c \\log \\rho_0,$$\n",
"via\n",
"$$ \\rho_0 = \\exp\\big[(u_0 - u_*)/c + \\log \\rho_* \\big].$$"
]
},
{
"cell_type": "markdown",
"id": "67b8903d",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Basic algorithm"
]
},
{
"cell_type": "markdown",
"id": "3717562c",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"1. Find $\\rho_*$\n",
"2. Use entropy condition (for shock speeds)\n",
" * If it's a shock: find $u_*$ using Rankine-Hugoniot\n",
" * If it's a rarefaction: use generalized Riemann invariant"
]
},
{
"cell_type": "markdown",
"id": "2cdf121d",
"metadata": {
"cell_style": "split"
},
"source": [
"![](../img/Toro-ShockOrRarefaction-2.19.png)"
]
},
{
"cell_type": "markdown",
"id": "a4f2a40f",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"## First a miracle happens\n",
"In general we will need to use a Newton method to solve for the state in the star region."
]
},
{
"cell_type": "markdown",
"id": "f4091dea",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Exact Riemann solver for isothermal gas dynamics"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "49ef36eb",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"dujump_isogas (generic function with 2 methods)"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function flux_isogas(U, c=1)\n",
" rho = U[1]\n",
" u = U[2] / rho\n",
" [U[1], U[1]*u + c^2*rho]\n",
"end\n",
"\n",
"function ujump_isogas(rho_L, rho_R, c=1)\n",
" if rho_R > rho_L # shock\n",
" (rho_L - rho_R) / sqrt(rho_L*rho_R)\n",
" else # rarefaction\n",
" c * (log(rho_L) - log(rho_R))\n",
" end\n",
"end\n",
"\n",
"function dujump_isogas(rho_L, drho_L, rho_R, drho_R, c=1)\n",
" if rho_R > rho_L # shock\n",
" ((drho_L - drho_R) / sqrt(rho_L*rho_R)\n",
" - .5*(rho_L - rho_R) * (rho_L*rho_R)^(-3/2) * (drho_L * rho_R + rho_L * drho_R))\n",
" else\n",
" c * (drho_L / rho_L - drho_R / rho_R)\n",
" end\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "97985958",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"riemann_isogas (generic function with 2 methods)"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function riemann_isogas(UL, UR, maxit=20)\n",
" rho_L, u_L = UL[1], UL[2]/UL[1]\n",
" rho_R, u_R = UR[1], UR[2]/UR[1]\n",
" rho = .5 * (rho_L + rho_R) # initial guess\n",
" for i in 1:maxit\n",
" f = (ujump_isogas(UL[1], rho)\n",
" + ujump_isogas(rho, UR[1])\n",
" - (UR[2]/UR[1] - UL[2]/UL[1]))\n",
" if norm(resid) < 1e-10\n",
" u = u_L + ujump_isogas(UL[1], rho)\n",
" break\n",
" end\n",
" J = (dujump_isogas(UL[1], 0, rho, 1)\n",
" + dujump_isogas(rho, 1, uR[1], 0))\n",
" delta_rho = -f / J\n",
" while min(rho + delta_rho <= 0)\n",
" delta_rho /= 2 # line search to prevent negative rho\n",
" end\n",
" end\n",
" U0 = resolve_isogas(UL, UR, rho)\n",
" flux_isogas(U0)\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "78ef5f14",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Resolving waves for isothermal gas dynamics"
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "d4e2f3b8",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"resolve_isogas (generic function with 1 method)"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function resolve_isogas(UL, UR, U_star)\n",
" rho_L, u_L = UL[1], UL[2]/UL[1]\n",
" rho_R, u_R = UR[1], UR[2]/UR[1]\n",
" rho, u = U_star[1], U_star[2] / U_star[1]\n",
" if u_L - c < 0 < u - c || u + c < 0 < u_R + c\n",
" # inside left (right) sonic rarefaction\n",
" u0 = sign(u) * c\n",
" rho0 = exp((u0-u)/c + log(rho))\n",
" return [rho0, rho0 * u0]\n",
" elseif ((rho_L >= rho && 0 <= u_L - c) ||\n",
" (rho_L < rho && 0 < rho*u - UL[2]))\n",
" return UL\n",
" elseif ((rho_R >= rho && u_R + c <= 0) ||\n",
" (rho_R < rho && UR[2] - rho*u > 0))\n",
" return UR\n",
" end\n",
" [rho, rho*u]\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 42,
"id": "6d864ecb",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function initial_isogas(x)\n",
" [1 .+ 2 * exp.(-4x .^ 2) 0*x]\n",
"end\n",
"\n",
"x = LinRange(-1, 1, 100)\n",
"U = initial_isogas(x)\n",
"plot(x, U)"
]
},
{
"cell_type": "markdown",
"id": "e73d27c1",
"metadata": {},
"source": [
"# Solver"
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "cd777176",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [
{
"data": {
"text/plain": [
"fv_solve2system (generic function with 3 methods)"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function fv_solve2system(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",
" U0 = u_init(x)\n",
" m, n = size(U0)\n",
" @show size(U0)\n",
" function rhs(t, U)\n",
" U = reshape(U, m, n)\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": 31,
"id": "fde9d804",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"size(U0) = (100,)\n"
]
},
{
"ename": "LoadError",
"evalue": "MethodError: no method matching +(::Vector{Float64}, ::Float64)\nFor element-wise addition, use broadcasting with dot syntax: array .+ scalar\n\u001b[0mClosest candidates are:\n\u001b[0m +(::Any, ::Any, \u001b[91m::Any\u001b[39m, \u001b[91m::Any...\u001b[39m) at operators.jl:591\n\u001b[0m +(\u001b[91m::T\u001b[39m, ::T) where T<:Union{Float16, Float32, Float64} at float.jl:383\n\u001b[0m +(::AbstractArray, \u001b[91m::StaticArraysCore.StaticArray\u001b[39m) at ~/.julia/packages/StaticArrays/6QFsp/src/linalg.jl:13\n\u001b[0m ...",
"output_type": "error",
"traceback": [
"MethodError: no method matching +(::Vector{Float64}, ::Float64)\nFor element-wise addition, use broadcasting with dot syntax: array .+ scalar\n\u001b[0mClosest candidates are:\n\u001b[0m +(::Any, ::Any, \u001b[91m::Any\u001b[39m, \u001b[91m::Any...\u001b[39m) at operators.jl:591\n\u001b[0m +(\u001b[91m::T\u001b[39m, ::T) where T<:Union{Float16, Float32, Float64} at float.jl:383\n\u001b[0m +(::AbstractArray, \u001b[91m::StaticArraysCore.StaticArray\u001b[39m) at ~/.julia/packages/StaticArrays/6QFsp/src/linalg.jl:13\n\u001b[0m ...",
"",
"Stacktrace:",
" [1] _broadcast_getindex_evalf",
" @ ./broadcast.jl:670 [inlined]",
" [2] _broadcast_getindex",
" @ ./broadcast.jl:643 [inlined]",
" [3] getindex",
" @ ./broadcast.jl:597 [inlined]",
" [4] copy",
" @ ./broadcast.jl:899 [inlined]",
" [5] materialize",
" @ ./broadcast.jl:860 [inlined]",
" [6] broadcast_preserving_zero_d",
" @ ./broadcast.jl:849 [inlined]",
" [7] +(A::Vector{Vector{Float64}}, Bs::Vector{Float64})",
" @ Base ./arraymath.jl:16",
" [8] ode_rk_explicit(f::var\"#rhs#6\"{typeof(riemann_isogas), Int64, typeof(limit_sin), Vector{Int64}, Vector{Int64}, Float64}, u0::Vector{Vector{Float64}}; tfinal::Float64, h::Float64, table::RKTable)",
" @ Main ./In[2]:31",
" [9] fv_solve2system(riemann::typeof(riemann_isogas), u_init::Function, n::Int64, tfinal::Float64, limit::Function)",
" @ Main ./In[30]:18",
" [10] fv_solve2system(riemann::Function, u_init::Function, n::Int64, tfinal::Float64)",
" @ Main ./In[30]:2",
" [11] top-level scope",
" @ In[31]:1",
" [12] eval",
" @ ./boot.jl:368 [inlined]",
" [13] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)",
" @ Base ./loading.jl:1428"
]
}
],
"source": [
"x, t_hist, U_hist = fv_solve2system(riemann_isogas, initial_isogas, 100, 0.1)\n",
"size(U_hist)"
]
},
{
"cell_type": "markdown",
"id": "522627ec",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Approximate Riemann solvers"
]
},
{
"cell_type": "markdown",
"id": "dd50fb6f",
"metadata": {
"cell_style": "split"
},
"source": [
"Exact Riemann solvers are\n",
"* complicated to implement\n",
"* fragile in the sense that small changes to the physics, such as in the equation of state $p(\\rho)$, can require changing many conditionals\n",
"* the need to solve for $\\rho^*$ using a Newton method and then evaluate each case of the wave structure can be expensive.\n",
"* An exact Riemann solver has never been implemented for some equations.\n",
"\n",
"## HLL (Harten, Lax, and van Leer)\n",
"\n",
"Assume two shocks with speeds $s_L$ and $s_R$. These speeds will be estimated and must be at least as fast as the fastest left- and right-traveling waves. If the wave speeds $s_L$ and $s_R$ are given, we have the Rankine-Hugoniot conditions across both shocks,"
]
},
{
"cell_type": "markdown",
"id": "6e98da72",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"\\begin{align}\n",
" s_L (U_* - U_L) &= f(U_*) - f(U_L) \\\\\n",
" s_R (U_R - U_*) &= f(U_R) - f(U_*) .\n",
"\\end{align}\n",
"\n",
"Adding these together gives\n",
"$$s_R U_R - (s_R - s_L) U_* - s_L U_L = f(U_R) - f(U_L)$$\n",
"which can be solved for $U_*$ as\n",
"$$ U_* = \\frac{f(U_L) - f(U_R) - s_L U_L + s_R U_R}{s_R - s_L} . $$\n",
"We can save an extra evaluation of the flux by using the Rankine-Hugoniot conditions\n",
"$$ f(U_*) = \\frac{s_R f(U_L) - s_L f(U_R) + s_L s_R (U_R - U_L)}{s_R - s_L} . $$\n",
"### Isothermal gas dynamics\n",
"\\begin{align} s_L &= \\min( u_L - c, u_R - c ) \\\\\n",
"s_R &= \\max(u_L + c, u_R + c)\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"id": "8a906089",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": "skip"
}
},
"source": [
"### Rusanov method\n",
"Special case $s_L = - s_R$, in which case the wave structure is always subsonic and the flux is simply\n",
"$$f(U_*) = \\frac 1 2 \\Big( f(U_L) + f(U_R) \\Big) - \\frac s 2 \\Big( U_R - U_L \\Big) . $$"
]
},
{
"cell_type": "markdown",
"id": "498a9410",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Observations on HLL solvers\n",
"\n",
"$$ f(U_*) = \\frac{s_R f(U_L) - s_L f(U_R) + s_L s_R (U_R - U_L)}{s_R - s_L} . $$\n",
"\n",
"* The term involving $U_R-U_L$ represents diffusion and will cause entropy to decay (physical entropy is produced).\n",
"* If our Riemann problem produces shocks and we have correctly calculated the wave speeds, the HLL solver is exact and produces the minimum diffusion necessary for conservation.\n",
"* If the wave speed estimates are slower than reality, the method will be unstable due to CFL.\n",
"* If the wave speed estimates are faster than reality, the method will be more diffusive than an exact Riemann solver."
]
}
],
"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.3"
},
"rise": {
"enable_chalkboard": true
}
},
"nbformat": 4,
"nbformat_minor": 5
}