{
"cells": [
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from IPython.core.display import HTML\n",
"css_file = './custom.css'\n",
"HTML(open(css_file, \"r\").read())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###### Content provided under a Creative Commons Attribution license, CC-BY 4.0; code under MIT License. (c)2014 [David I. Ketcheson](http://davidketcheson.info)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"##### version 0.1 - May 2014"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# High-resolution methods"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\newcommand{Dx}{\\Delta x}\n",
"\\newcommand{Dt}{\\Delta t}\n",
"\\newcommand{imh}{{i-1/2}}\n",
"\\newcommand{iph}{{i+1/2}}\n",
"\\newcommand{\\DQ}{\\Delta Q}\n",
"\\DeclareMathOperator{\\sgn}{sgn}\n",
"$$\n",
"The methods we have used so far (the *upwind method* and the *Lax-Friedrichs method*) are both dissipative. This means that over time they artificially smear out the solution -- especially shocks. Furthermore, both of these methods are only *first order accurate*, meaning that if we reduce the values of $\\Dt$ and $\\Dx$ by a factor of two, the overall error decreases (only) by a factor of two. We can do better."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Reducing dissipation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The first step in improving our accuracy is to use a more accurate representation of $q(x)$ at each step. Instead of assuming that $q$ is piecewise-constant, we will now approximate it by a piecewise-linear function:\n",
"\n",
"$$q(x) = Q^n_i + \\sigma^n_i (x-x_i).$$\n",
"\n",
"Here $\\sigma_i$ represents the slope in cell $i$. The most obvious choice to ensure that this results in a second-order accurate approximation is to take the centered approximation\n",
"\n",
"$$\\sigma^n_i = \\frac{Q^n_{i+1} - Q^n_{i-1}}{2\\Dx}.$$\n",
"\n",
"We use this to obtain values at the cell interfaces:\n",
"\n",
"\\begin{align}\n",
"q^+_\\imh & = Q_i - \\sigma_i \\frac{\\Dx}{2} \\\\\n",
"q^-_\\iph & = Q_i + \\sigma_i \\frac{\\Dx}{2}.\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"![](./figures/linear_reconstruction.png)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We'll use these interface values to approximate the flux, based on the **Lax-Friedrichs flux**:\n",
"\n",
"$$F_\\imh = \\frac{1}{2} \\left( f(q^-_\\imh) + f(q^+_\\imh) - \\frac{\\Dx}{\\Dt} (q^+_\\imh - q^-_\\imh)\\right)$$\n",
"\n",
"This provides second-order accuracy in space. We also need to make the method second-order accurate in time. We can do so by using a second-order Runge--Kutta method. Then the full method is\n",
"\n",
"\\begin{align}\n",
"Q^*_i & = Q^n_i - \\frac{\\Dt}{\\Dx}\\left( F^n_\\iph - F^n_\\imh \\right) \\\\\n",
"Q^{n+1}_i & = \\frac{1}{2} Q^n_i + \\frac{1}{2}\\left( Q^*_i - \\frac{\\Dt}{\\Dx}\\left( F^*_\\iph - F^*_\\imh \\right) \\right) \\\\\n",
"\\end{align}"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"import sys\n",
"sys.path.append('./util')\n",
"from ianimate import ianimate\n",
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"
\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def f(q):\n",
" return q*(1.0-q)\n",
"\n",
"m = 100 # number of points\n",
"dx = 1./m # Size of 1 grid cell\n",
"x = np.arange(-3*dx/2, 1.+5*dx/2, dx)\n",
"\n",
"t = 0. # Initial time\n",
"T = 0.5 # Final time\n",
"dt = 0.9 * dx # Time step\n",
"\n",
"Q = 0.9*np.exp(-100*(x-0.5)**2)\n",
"Qnew = np.zeros(m+4)\n",
"Qstar = np.zeros(m+4)\n",
"sigma = np.zeros(m+4)\n",
"\n",
"F = np.zeros(m+4)\n",
"QQ = [Q]\n",
"\n",
"while t < T:\n",
" \n",
" sigma[1:-1] = (Q[2:] - Q[:-2])/(2.0*dx)\n",
" qplus = Q[1:-1] - sigma[1:-1] * dx/2.0 # q^+_{i-1/2}\n",
" qminus = Q[:-2] + sigma[:-2] * dx/2.0 # q^-_{i-1/2}\n",
" F[1:-1] = 0.5*(f(qplus)+f(qminus) - dx/dt*(qplus-qminus) )# F_{i-1/2}\n",
" \n",
" Qstar[2:-2] = Q[2:-2] - dt/dx*(F[3:-1]-F[2:-2])\n",
" Qstar[0:2] = Qstar[2]\n",
" Qstar[-2:] = Qstar[-3]\n",
" \n",
" sigma[1:-1] = (Qstar[2:] - Qstar[:-2])/(2.0*dx)\n",
" qplus = Qstar[1:-1] - sigma[1:-1] * dx/2.0 # q^+_{i-1/2}\n",
" qminus = Qstar[:-2] + sigma[:-2] * dx/2.0 # q^-_{i-1/2}\n",
" F[1:-1] = 0.5*(f(qplus)+f(qminus) - dx/dt*(qplus-qminus) )# F_{i-1/2}\n",
" \n",
" Qnew[2:-2] = 0.5*Q[2:-2] + 0.5*(Qstar[2:-2] - dt/dx*(F[3:-1]-F[2:-2]))\n",
" \n",
" Q = Qnew.copy()\n",
" Q[0:2] = Q[2]\n",
" Q[-2:] = Q[-3]\n",
" t = t + dt\n",
" QQ.append(Q)\n",
" \n",
"ianimate(x,QQ)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The shock wave is much sharper now, but we have a new problem. Do you see the little dip behind the shock? If you look closely, you'll see that the solution is actually negative there! Obviously, a negative density of cars makes no sense. What's more, the solution shouldn't dip there -- and it shouldn't have that funny bump just in front of the shock either. What to do?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Slope limiting"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The spurious oscillations in our solution are not a particular feature of the method we've chosen. In fact, *any* second-order (or higher) method that computes $q^\\pm_\\imh$ as a linear function of the cell averages will have oscillations (this is a famous result known as *Godunov's Theorem*).\n",
"\n",
"We can get around this difficulty by choosing the slope $\\sigma$ as a *nonlinear* function of the cell averages. In particular, to avoid oscillations we can choose the smaller of the two one-sided slopes. Let $\\DQ_\\imh = Q_i - Q_{i-1}$. Then we use the slope\n",
"\n",
"\\begin{align}\n",
"\\sigma_i & = \\text{minmod}(\\DQ_\\imh,\\DQ_\\iph)/\\Dx \\\\\n",
"& = \\begin{cases} \\min(\\DQ_\\imh, \\DQ_\\iph)/\\Dx & \\text{ if } \\DQ_\\imh, \\DQ_\\iph > 0 \\\\\n",
"\\max(\\DQ_\\imh, \\DQ_\\iph)/\\Dx & \\text{ if } \\DQ_\\imh, \\DQ_\\iph < 0 \\\\\n",
"0 & \\text{ if } \\DQ_\\imh\\cdot \\DQ_\\iph < 0.\n",
"\\end{cases}\n",
"\\end{align}\n",
"\n",
"This choice of slope is known as the minimum-modulus, or *minmod* slope."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Local Lax-Friedrichs flux"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Lax-Friedrichs flux ensures that our solution is stable, but it does so by adding a lot of dissipation everywhere. In fact, we could get away with using less dissipation over most of the domain. A variant that does this is called the **local Lax-Friedrichs flux**. It is little more accurate than the Lax-Friedrichs flux because it uses the local characteristic speeds to determine how much dissipation is needed at each interface. It is\n",
"\n",
"$$F_\\imh = \\frac{1}{2} \\left( f(q^-_\\imh) + f(q^+_\\imh) - \\alpha_\\imh\\frac{\\Dx}{\\Dt} (q^+_\\imh - q^-_\\imh)\\right)$$\n",
"\n",
"where\n",
"\n",
"$$\\alpha_\\imh = \\min(\\left|f'(q^-_\\imh)\\right|,\\left|f'(q^+_\\imh)\\right|).$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exercise"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the cell below, copy and modify the second-order method above to use the minmod slope and local Lax-Friedrichs flux. \n",
"\n",
"*Hint 1*: You will want to use the functions `np.minimum` and `np.maximum`, which compare two arrays elementwise (not `np.min`, which finds the minimum of a single array).\n",
"\n",
"*Hint 2*: to avoid using a loop for the slope computation, note that the minmod function can be written as\n",
"$$\n",
"\\text{minmod}(a,b) = \\frac{\\sgn(a)+\\sgn(b)}{2} \\min(|a|,|b|).\n",
"$$\n",
"The signum function is implemented as `np.sign()`. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As you can see (after completing the exercise), this method keeps the shock fairly sharp and avoids the creation of negative solution values. This method falls into the class of [MUSCL]() schemes and is proven to avoid oscillations."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## High-order methods"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Just like the method we implemented above, most methods that are more than first-order accurate consist of three components:\n",
"1. **Reconstruction**: A method for computing interface values $q^\\pm_\\imh$ from cell averages $Q_i$. This involves some kind of limiting in order to avoid oscillations. Higher-order reconstruction is often done using [weighted essentially non-oscillatory (WENO)](http://www.dam.brown.edu/scicomp/media/report_files/BrownSC-2006-21.ps.gz) methods.\n",
"2. **Numerical flux**: An approximation of the flux, computed based on the interface values $q^\\pm_\\imh$. The Lax-Friedrichs flux above is one of the simplest. Much more accurate fluxes can be computed using an exact or approximate [Riemann solver](http://en.wikipedia.org/wiki/Riemann_solver).\n",
"3. **Time integrator**: In order to get high-order accuracy in time, usually a [Runge-Kutta method](http://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods) is used. [Strong stability preserving methods](http://www.davidketcheson.info/assets/papers/sspreview.pdf) are particularly popular."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It's possible to devise methods of very high order by increasing the order of the polynomial reconstruction and of the time integrator."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.9"
}
},
"nbformat": 4,
"nbformat_minor": 0
}