{ "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", "
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "\n", " \n", " \n", " Once \n", " Loop \n", " Reflect \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 }