{
"cells": [
{
"cell_type": "markdown",
"id": "4159be0e",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# 2021-09-13 Variable coefficients\n",
"\n",
"## Last time\n",
"* Conditioning and stability of polynomial interpolation\n",
"* Runge Phenomenon\n",
"* Chebyshev differentiation\n",
"* Chebyshev-based solver for the Poisson problem\n",
"\n",
"## Today\n",
"* Variable coefficients\n",
"* Conservative/divergence form vs non-divergence forms\n",
"* Verification with discontinuities"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "009bd630",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [
{
"data": {
"text/plain": [
"chebdiff (generic function with 2 methods)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"using Plots\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",
"\n",
"function poisson(x, spoints, forcing; left=(0, zero), right=(0, zero))\n",
" n = length(x)\n",
" L = zeros(n, n)\n",
" rhs = forcing.(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",
" L[i, js] = -fdstencil(x[js], x[i], 2)\n",
" end\n",
" L[1,1:spoints] = fdstencil(x[1:spoints], x[1], left[1])\n",
" L[n,n-spoints+1:n] = fdstencil(x[n-spoints+1:n], x[n], right[1])\n",
" rhs[1] = left[2](x[1])\n",
" rhs[n] = right[2](x[n])\n",
" L, rhs\n",
"end\n",
"\n",
"CosRange(a, b, n) = (a + b)/2 .+ (b - a)/2 * cos.(LinRange(-pi, 0, n))\n",
"\n",
"function vander_chebyshev(x, n=nothing)\n",
" if isnothing(n)\n",
" n = length(x) # Square by default\n",
" end\n",
" m = length(x)\n",
" T = ones(m, n)\n",
" if n > 1\n",
" T[:, 2] = x\n",
" end\n",
" for k in 3:n\n",
" T[:, k] = 2 * x .* T[:,k-1] - T[:, k-2]\n",
" end\n",
" T\n",
"end\n",
"\n",
"function chebdiff(x, n=nothing)\n",
" T = vander_chebyshev(x, n)\n",
" m, n = size(T)\n",
" dT = zero(T)\n",
" dT[:,2:3] = [one.(x) 4*x]\n",
" for j in 3:n-1\n",
" dT[:,j+1] = j * (2 * T[:,j] + dT[:,j-1] / (j-2))\n",
" end\n",
" ddT = zero(T)\n",
" ddT[:,3] .= 4\n",
" for j in 3:n-1\n",
" ddT[:,j+1] = j * (2 * dT[:,j] + ddT[:,j-1] / (j-2))\n",
" end\n",
" T, dT, ddT\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "b3462747",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Lagrange interpolating polynomials on `CosRange`\n",
"\n",
"\"nodal\" --- \"modal\""
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "f51dc7a2",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"n = 5\n",
"x = CosRange(-1, 1, n)\n",
"xx = LinRange(-1, 1, 100)\n",
"Tx = vander_chebyshev(x)\n",
"Txx, dTxx, ddTxx = chebdiff(xx, n)\n",
"plot(xx, Txx / Tx)\n",
"scatter!(x, [one.(x) zero.(x)])"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "8becaa50",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot(xx, Txx)"
]
},
{
"cell_type": "markdown",
"id": "f0e71078",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Solving a BVP on Chebyshev nodes"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "e759ef71",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"text/plain": [
"poisson_cheb (generic function with 3 methods)"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function poisson_cheb(n, 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\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": 7,
"id": "64fb16ff",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"manufactured(x) = tanh(2x)\n",
"d_manufactured(x) = 2*cosh(2x)^-2\n",
"mdd_manufactured(x) = 8 * tanh(2x) / cosh(2x)^2\n",
"x, A, rhs = poisson_cheb(12, mdd_manufactured,\n",
" (0, manufactured), (1, d_manufactured))\n",
"plot(x, A \\ rhs, marker=:auto)\n",
"plot!(manufactured, legend=:bottomright)"
]
},
{
"cell_type": "markdown",
"id": "16543cbb",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# \"spectral\" (exponential) convergence"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "a35af886",
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function poisson_error(n)\n",
" x, A, rhs = poisson_cheb(n, mdd_manufactured, (0, manufactured), (1, d_manufactured))\n",
" u = A \\ rhs\n",
" norm(u - manufactured.(x), Inf)\n",
"end\n",
"\n",
"ns = 3:20\n",
"plot(ns, abs.(poisson_error.(ns)), marker=:auto, yscale=:log10)\n",
"ps = [1 2 3]\n",
"plot!([n -> n^-p for p in ps], label=map(p -> \"\\$n^{-$p}\\$\", ps))"
]
},
{
"cell_type": "markdown",
"id": "e1c32aa4",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Variable coefficients"
]
},
{
"cell_type": "markdown",
"id": "21a362be",
"metadata": {
"cell_style": "split",
"slideshow": {
"slide_type": ""
}
},
"source": [
"* Heat conduction: steel, brick, wood, foam\n",
"* Electrical conductivity: copper, rubber, air\n",
"* Elasticity: steel, rubber, concrete, rock\n",
"* Linearization of nonlinear materials\n",
" * ketchup, glacier ice, rocks (mantle/lithosphere)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "1286ddd0",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"kappa_step(x) = .1 + .9 * (x > 0)\n",
"kappa_smooth(x) = .55 + .45 * sin(pi*x/2)\n",
"plot([kappa_step, kappa_smooth], xlims=(-1, 1), ylims=(0, 1), label=\"κ\")"
]
},
{
"cell_type": "markdown",
"id": "c0f1667d",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"\\begin{align}\n",
"-(\\kappa u_x)_x &= 0 & u(-1) &= 0 & \\kappa u_x(1) &= 1\n",
"\\end{align}\n",
"\n",
"* What physical scenario could this represent?\n",
"* Qualitatively, what would a solution look like?"
]
},
{
"cell_type": "markdown",
"id": "bd4bb0f7",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# A naive finite difference solver\n",
"\n",
"| Conservative (divergence) form | Non-divergence form |\n",
"| ------------------------------ | ------------------- |\n",
"| $-(\\kappa u_x)_x = 0$ | $-\\kappa u_{xx} - \\kappa_x u_x = 0$ |"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "81c405d9",
"metadata": {
"cell_style": "center",
"slideshow": {
"slide_type": ""
}
},
"outputs": [
{
"data": {
"text/plain": [
"poisson_nondivergence (generic function with 1 method)"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function poisson_nondivergence(x, spoints, kappa, forcing; leftbc=(0, zero), rightbc=(0, zero))\n",
" n = length(x)\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",
" L, rhs\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "c6fa1017",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Try it"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "b95036fe",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = LinRange(-1, 1, 30)\n",
"L, rhs = poisson_nondivergence(x, 3, kappa_smooth, zero, rightbc=(1, one))\n",
"u = L \\ rhs\n",
"plot(x, u)\n",
"plot!(x -> 5*kappa_smooth(x))"
]
},
{
"cell_type": "markdown",
"id": "3fa28122",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Manufactured solutions for variable coefficients"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "d71b825f",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"manufactured(x) = tanh(2x)\n",
"d_manufactured(x) = 2*cosh(2x)^-2\n",
"d_kappa_smooth(x) = .45*pi/2 * cos(pi*x/2)\n",
"flux_manufactured_kappa_smooth(x) = kappa_smooth(x) * d_manufactured(x)\n",
"function forcing_manufactured_kappa_smooth(x)\n",
" 8 * tanh(2x) / cosh(2x)^2 * kappa_smooth(x) -\n",
" d_kappa_smooth(x) * d_manufactured(x)\n",
"end\n",
"x = LinRange(-1, 1, 20)\n",
"L, rhs = poisson_nondivergence(x, 3, kappa_smooth,\n",
" forcing_manufactured_kappa_smooth,\n",
" leftbc=(0, manufactured), rightbc=(1, flux_manufactured_kappa_smooth))\n",
"u = L \\ rhs\n",
"plot(x, u, marker=:auto, legend=:bottomright)\n",
"\n",
"plot!([manufactured flux_manufactured_kappa_smooth forcing_manufactured_kappa_smooth kappa_smooth])"
]
},
{
"cell_type": "markdown",
"id": "b5449e60",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Convergence"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "2c0e80c9",
"metadata": {
"slideshow": {
"slide_type": ""
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function poisson_error(n, spoints=3)\n",
" x = LinRange(-1, 1, n)\n",
" L, rhs = poisson_nondivergence(x, spoints, kappa_smooth,\n",
" forcing_manufactured_kappa_smooth,\n",
" leftbc=(0, manufactured), rightbc=(1, flux_manufactured_kappa_smooth))\n",
" u = L \\ rhs\n",
" norm(u - manufactured.(x), Inf)\n",
"end\n",
"ns = 2 .^ (3:10)\n",
"plot(ns, poisson_error.(ns, 5), marker=:auto, xscale=:log10, yscale=:log10)\n",
"plot!([n -> n^-2, n -> n^-4], label=[\"\\$1/n^2\\$\" \"\\$1/n^4\\$\"])"
]
},
{
"cell_type": "markdown",
"id": "31bcee39",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# ✅ Verified!"
]
},
{
"cell_type": "markdown",
"id": "9d90b570",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Let's try with the discontinuous coefficients"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "08fc1386",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = LinRange(-1, 1, 20)\n",
"L, rhs = poisson_nondivergence(x, 3, kappa_step,\n",
" zero,\n",
" leftbc=(0, zero), rightbc=(0, one))\n",
"u = L \\ rhs\n",
"plot(x, u, marker=:auto, legend=:bottomright)"
]
},
{
"cell_type": "markdown",
"id": "f7eef608",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Discretizing in conservative form\n",
"\n",
"| Conservative (divergence) form | Non-divergence form |\n",
"| ------------------------------ | ------------------- |\n",
"| $-(\\kappa u_x)_x = 0$ | $-\\kappa u_{xx} - \\kappa_x u_x = 0$ |"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "4dd13b8a",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"poisson_conservative (generic function with 1 method)"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"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": "2cc4b9a8",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Compare conservative vs non-divergence forms"
]
},
{
"cell_type": "code",
"execution_count": 41,
"id": "29f4a24f",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"forcing = zero # one\n",
"x, L, rhs = poisson_conservative(20, kappa_step,\n",
" forcing, leftbc=(0, zero), rightbc=(0, one))\n",
"u = L \\ rhs\n",
"plot(x, u, marker=:auto, legend=:bottomright)"
]
},
{
"cell_type": "code",
"execution_count": 42,
"id": "9f8fecb0",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = LinRange(-1, 1, 200)\n",
"L, rhs = poisson_nondivergence(x, 3, kappa_step,\n",
" forcing, leftbc=(0, zero), rightbc=(0, one))\n",
"u = L \\ rhs\n",
"plot(x, u, marker=:auto, legend=:bottomright)"
]
},
{
"cell_type": "markdown",
"id": "f5e28364",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Continuity of flux"
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "d22860ce",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"forcing = zero\n",
"x, L, rhs = poisson_conservative(20, kappa_step,\n",
" forcing, leftbc=(0, zero), rightbc=(0, one))\n",
"u = L \\ rhs\n",
"plot(x, u, marker=:auto, legend=:bottomright)"
]
},
{
"cell_type": "code",
"execution_count": 29,
"id": "7d9e5788",
"metadata": {
"cell_style": "split"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xstag = (x[1:end-1] + x[2:end]) ./ 2\n",
"du = (u[1:end-1] - u[2:end]) ./ diff(x)\n",
"plot(xstag, [du .* kappa_step.(xstag)], marker=:auto, ylims=[-1, 1])"
]
},
{
"cell_type": "markdown",
"id": "89df9254",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Manufactured solutions with discontinuous coefficients\n",
"\n",
"* We need to be able to evaluate derivatives of the flux $-\\kappa u_x$.\n",
"* A physically-realizable solution would have continuous flux, but we we'd have to be making a physical solution to have that in verification.\n",
"* Idea: replace the discontinuous function with a continuous one with a rapid transition."
]
},
{
"cell_type": "code",
"execution_count": 44,
"id": "dd4d0e83",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 44,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"kappa_tanh(x, epsilon=.1) = .55 + .45 * tanh(x / epsilon)\n",
"d_kappa_tanh(x, epsilon=.1) = .45/epsilon * cosh(x/epsilon)^-2\n",
"plot([kappa_tanh])"
]
},
{
"cell_type": "markdown",
"id": "443c90e4",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Solving with the smoothed step $\\kappa$"
]
},
{
"cell_type": "code",
"execution_count": 35,
"id": "b6f536fe",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"