Guide_to_limiters.ipynb

{
 "cells": [
  {
   "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": [
    "# An illustrated guide to limiters\n",
    "## Or: how to interpolate non-smooth data without creating wiggles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many interesting wave phenomena -- like fluid dynamics, lasers, and water waves -- are described by nonlinear hyperbolic partial differential equations.  The solutions of these problems are discontinuous.  So-called **limiters** (sometimes referred to as *slope limiters* or *flux limiters* are one of the key ingredients in approximating these discontinuous solutions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of contents\n",
    "- [Motivation: interpolation and wiggles](#Interpolation-and-wiggles)\n",
    "- [The simplest limiter: Minmod](#The-simplest-limiter:-Minmod)\n",
    "- [Other TVD limiters](#TVD-limiters)\n",
    "- [WENO](#Higher-order-interpolation:-WENO)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "matplotlib.rcParams.update({'font.size': 22})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mpld3             # Skip this cell if you don't have mpld3 installed\n",
    "mpld3.enable_notebook()  # or just go and do it now: pip install mpld3\n",
    "                         # This allows you to zoom and pan the plots."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpolation and wiggles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose you're given a set of data samples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 5\n",
    "x=np.arange(-k+1,k)\n",
    "y=np.sin(x/2.)+1.\n",
    "width = 12\n",
    "size = (width,4)\n",
    "plt.figure(figsize=size)\n",
    "plt.plot(x,y,'or',markersize=10,alpha=0.5)\n",
    "plt.axis( (-k, k, -0.1, 2.1) );"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now what you really want to know is, what is the state of the system at the points halfway between your samples?  And to figure that out, you need to guess at what's going on in the times in-between those samples.  The simplest approximation would be to assume that the system just jumps from one value to the next somewhere in-between:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def piecewise_constant_interp(x,y,xx):\n",
    "    \"From samples (x,y) generate piecewise constant function sampled at points xx.\"\n",
    "    diff = np.abs(x.reshape(1,-1) - xx.reshape(-1,1))  # Here we use numpy broadcasting.\n",
    "    closest = np.argmin(diff,axis=1)\n",
    "    return y[closest]\n",
    "\n",
    "xx = np.linspace(-k+1,k-1,1000)\n",
    "yy = piecewise_constant_interp(x,y,xx)\n",
    "plt.figure(figsize=size)\n",
    "plt.plot(xx,yy,'-k',lw=2)\n",
    "plt.plot(x,y,'or',markersize=10,alpha=0.5)\n",
    "plt.axis( (-k, k, -0.1, 2.1) );\n",
    "plt.title('Piecewise-constant approximation',fontsize=20);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this set of data, you don't really believe that's what's happening, do you?  But our goal is to deal with systems that exhibit non-smooth (possibly discontinuous) behavior, so we need to at least admit the possibility of sudden jumps.  That's why we won't simply \"connect the dots\" to get a continuous approximation.\n",
    "\n",
    "Instead, we can try to approximate the slope around each of our sample points.  The simplest way to do so is using finite differences.  If we let $\\sigma_i$ denote our approximation of the slope at $x_i$, then three common approximations are:\n",
    "\n",
    "- Forward difference: $\\sigma_i = \\frac{y_{i+1}-y_i}{x_{i+1}-x_i}$\n",
    "\n",
    "\n",
    "- Backward difference: $\\sigma_i = \\frac{y_i - y_{i-1}}{x_i - x_{i-1}}$\n",
    "\n",
    "\n",
    "- Centered difference: $\\sigma_i = \\frac{y_{i+1}-y_{i-1}}{x_{i+1}-x_{i-1}}$\n",
    "\n",
    "Here's what each of these approximations looks like for our data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def piecewise_linear_interp(x,y,xx, fd='centered'):\n",
    "    \"From samples (x,y) generate piecewise-linear function sampled at points xx using finite difference slopes.\"\n",
    "    diff = np.abs(x.reshape(1,-1) - xx.reshape(-1,1))\n",
    "    closest = np.argmin(diff,axis=1)\n",
    "    \n",
    "    sigma = np.zeros_like(y)\n",
    "    if fd == 'centered':\n",
    "        sigma[1:-1] = (y[2:]-y[:-2])/(x[2:]-x[:-2])\n",
    "    elif fd == 'forward':\n",
    "        sigma[:-1] = (y[1:]-y[:-1])/(x[1:]-x[:-1])\n",
    "    elif fd == 'backward':\n",
    "        sigma[1:] = (y[1:]-y[:-1])/(x[1:]-x[:-1])\n",
    "    return y[closest] + sigma[closest]*(xx-x[closest])\n",
    "\n",
    "def compare_fd(x,y,xx, axis=(-4, 4, -0.1, 2.1)):\n",
    "    fig, ax = plt.subplots(3,1,figsize=(width,8))\n",
    "    for i, fd in enumerate( ('centered','forward','backward') ):\n",
    "        yy = piecewise_linear_interp(x,y,xx,fd=fd)\n",
    "        ax[i].plot(xx,yy,'-k',lw=2)\n",
    "        ax[i].plot(x,y,'or',markersize=10,alpha=0.5)\n",
    "        ax[i].axis( axis );\n",
    "        ax[i].text(.5,.9,fd,\n",
    "            horizontalalignment='center',\n",
    "            transform=ax[i].transAxes,fontsize=20)\n",
    "        \n",
    "compare_fd(x,y,xx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I've used $\\sigma=0$ for the points at the edges where we don't have enough data to compute the appropriate slope."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The problem: overshoots!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking closely, you can see that each of these approximations adds little jumps (called \"overshoots\") in some region where the data itself was monotone.  Worse still, each of them generates negative values, whereas the original values were non-negative!  If our data represent concentrations or probabilities, then we have no way to make sense of negative values.\n",
    "\n",
    "Things look even worse if we take data samples from a function that is in fact discontinuous:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.sin(x/2.)+1. + 2.*(x>0)\n",
    "compare_fd(x,y,xx,axis=(-4,4,-0.5,4.8))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now all three approaches have large, obvious overshoots.  This becomes even more problematic when actually solving a hyperbolic PDE; see [Lesson 3 of my HyperPython course](http://nbviewer.ipython.org/github/ketch/HyperPython/blob/master/Lesson_03_High-resolution_methods.ipynb) for details.\n",
    "\n",
    "Is there a better way?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The simplest limiter: Minmod"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'd like to avoid those overshoots and ensure that monotone regions of the data give monotone interpolations.  We can do that by choosing the slope $\\sigma_i$ small enough that the interpolant near $x_i$ stays bounded between the neighboring averages $(y_{i-1}+y_i)/2$ and $(y_i+y_{i+1})/2$.  There's an easy way to do that: just compute the forward and backward differences (like we did above), and then use *whichever is smaller* in absolute value.  If $y_i$ is an extremum, then to avoid increasing the overall range of the data we always choose $\\sigma_i=0$.\n",
    "\n",
    "Here's what that looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pw_minmod(x,y,xx):\n",
    "    \"From samples (x,y) generate piecewise-linear function sampled at points xx using Minmod slopes.\"\n",
    "    diff = np.abs(x.reshape(1,-1) - xx.reshape(-1,1))\n",
    "    closest = np.argmin(diff,axis=1)\n",
    "    \n",
    "    forward  = np.zeros_like(y)\n",
    "    backward = np.zeros_like(y)\n",
    "    sigma    = np.zeros_like(y)\n",
    "    \n",
    "    forward[:-1] = (y[1:]-y[:-1])/(x[1:]-x[:-1])\n",
    "    backward[1:] = (y[1:]-y[:-1])/(x[1:]-x[:-1])\n",
    "\n",
    "    sigma = (np.sign(forward)+np.sign(backward))/2. * np.minimum(np.abs(forward),np.abs(backward))\n",
    "    \n",
    "    return y[closest] + sigma[closest]*(xx-x[closest])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "yy = pw_minmod(x,y,xx)\n",
    "plt.figure(figsize=size)\n",
    "plt.plot(xx,yy,'-k',lw=2)\n",
    "plt.plot(x,y,'or',markersize=10,alpha=0.5)\n",
    "plt.axis( (-4,4,-0.5,4.8) );\n",
    "plt.title('Minmod approximation',fontsize=20);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's apply minmod to a monotone sequence of values, to illustrate the average-boundedness property:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.patches import Rectangle\n",
    "\n",
    "y = np.exp(x/3.)\n",
    "yy = pw_minmod(x,y,xx)\n",
    "plt.figure(figsize=(width,6))\n",
    "plt.plot(xx,yy,'-k',lw=2)\n",
    "plt.plot(x,y,'or',markersize=10,alpha=0.5)\n",
    "plt.axis( (-4,4,-0.1,4.1) );\n",
    "plt.title('minmod approximation',fontsize=20);\n",
    "for i in range(len(y)-1):\n",
    "    if 1<=i<len(y):\n",
    "        x_avgs = [(x[i]+x[i-1])/2.,(x[i]+x[i+1])/2.]\n",
    "        y_avgs = [(y[i]+y[i-1])/2.,(y[i]+y[i+1])/2.]\n",
    "        currentAxis = plt.gca()\n",
    "        currentAxis.add_patch(Rectangle((x_avgs[0], y_avgs[0]), \n",
    "                                         x_avgs[1]-x_avgs[0], y_avgs[1]-y_avgs[0], \n",
    "                                         facecolor=\"grey\",alpha=0.5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The grey regions show the average of each value and its neighbors; the Minmod interpolant is guaranteed to stay within these bounds, so that the right edge of one interpolated region doesn't \"overshoot\" the leftmost value of the next."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Total variation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In fact, the minmod limiter guarantees an important mathematical property: it does not increase the total variation (TV) of the data.  For discrete data, TV is defined as\n",
    "\n",
    "$$ TV(y) = \\sum_i |y_i - y_{i-1}|.$$\n",
    "\n",
    "Essentially, TV is a measure of how much a function wiggles!  This *total variation diminishing* (TVD) property of minmod is nice, because the true solutions of scalar hyperbolic PDEs have the same property.  Because of this, in the 1970s-80s, researchers developed many other limiters that enforce the same property; they are referred to as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TVD limiters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just for fun, let's look at 3 more such limiters.  If we define $F_i, B_i$ as forward- and backward-difference slope approximations at point $i$, and  $\\theta_i$ as the ratio of the two,  then each of these limiters can be described in terms of a function $\\phi(\\theta)$.  The slope to be used is given by\n",
    "\n",
    "$$\\sigma_i = \\phi(\\theta_i) B_i.$$\n",
    "\n",
    "The limiter names and functions are:\n",
    "\n",
    "- Monotonized centered difference (**MC**): $$\\phi(\\theta) = \\max(0,\\min((1+\\theta)/2,2,\\theta))$$\n",
    "\n",
    "- **van Leer**: $$\\phi(\\theta) = \\frac{\\theta + |\\theta|}{1+ |\\theta|}$$\n",
    "\n",
    "- **Superbee**: $$\\phi(\\theta) = \\max(0,\\min(1,2\\theta),\\min(2,\\theta))$$\n",
    "\n",
    "Here's what each of them looks like for our data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def phi(theta,limiter):\n",
    "    if limiter == 'minmod':\n",
    "        phi = (1+np.sign(theta))/2. * np.minimum(1,theta)\n",
    "    elif limiter == 'vanleer':\n",
    "        phi = (theta + np.abs(theta))/(1+np.abs(theta))\n",
    "    elif limiter == 'MC':\n",
    "        phi = np.maximum(0,np.minimum( (1.+theta)/2., np.minimum(2.,theta)))\n",
    "    elif limiter == 'superbee':\n",
    "        phi = np.maximum(0,np.maximum(np.minimum(1.,2*theta),np.minimum(2.,theta)))\n",
    "    return phi\n",
    "\n",
    "def pw_limited(x,y,xx,limiter='minmod'):\n",
    "    \"From samples (x,y) generate piecewise-linear function sampled at points xx using Minmod slopes.\"\n",
    "    diff = np.abs(x.reshape(1,-1) - xx.reshape(-1,1))\n",
    "    closest = np.argmin(diff,axis=1)\n",
    "    \n",
    "    forward  = np.zeros_like(y)\n",
    "    backward = np.zeros_like(y)\n",
    "    theta    = np.zeros_like(y)\n",
    "    \n",
    "    forward[:-1] = (y[1:]-y[:-1])/(x[1:]-x[:-1])\n",
    "    backward[1:] = (y[1:]-y[:-1])/(x[1:]-x[:-1])\n",
    "    theta[1:-1] = forward[1:-1]/backward[1:-1]\n",
    "    \n",
    "    sigma = phi(theta,limiter) * backward\n",
    "        \n",
    "    return y[closest] + sigma[closest]*(xx-x[closest])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.sin(x/2.)+1. + 2.*(x>0)\n",
    "fig, ax = plt.subplots(4,1,figsize=(width,10))\n",
    "for i, limiter in enumerate( ('minmod', 'vanleer','superbee','MC') ):\n",
    "    yy = pw_limited(x,y,xx,limiter=limiter)\n",
    "    ax[i].plot(xx,yy,'-k',lw=2)\n",
    "    ax[i].plot(x,y,'or',markersize=10,alpha=0.5)\n",
    "    ax[i].axis( (-4,4,-0.1,4.4) );\n",
    "    ax[i].text(.8,.2,limiter,\n",
    "        horizontalalignment='center',\n",
    "        transform=ax[i].transAxes,fontsize=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare these with [the finite difference approximations above](#The-problem:-overshoots!).\n",
    "\n",
    "If you look closely (or zoom in) you'll notice that -- except for minmod -- all the limiters *do* produce some overshoot near the discontinuity.  What gives?  Well, these limiters are used within a larger algorithm for solving hyperbolic PDEs, and it turns out that if the overshoots are small enough, they'll go away in a full step of the algorithm.  These limiters produce \"small enough\" overshoots so that no oscillations appear in the PDE solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interactive comparison\n",
    "In each region, these limiter take three data points and give back a linear interpolant.  It's illuminating to compare their behavior on a single set of 3 points.  Note that the interactive plot below doesn't work on nbviewer; you'll need to download and run the notebook yourself."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipywidgets import interact, FloatSlider, RadioButtons\n",
    "from IPython.display import display\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xx = np.linspace(-0.5,0.5)\n",
    "\n",
    "def compare_limiters(y1,y3):\n",
    "    fig, ax = plt.subplots(figsize=(width,4))\n",
    "    x = np.array((-1.,0.,1.))\n",
    "    y = np.array((y1,0.,y3))\n",
    "    ax.set_xlim(-1.1,1.1)\n",
    "    ax.set_ylim(-1.1,1.1)\n",
    "    \n",
    "    if y1 == 0:\n",
    "        theta = y3\n",
    "    else:\n",
    "        theta = y3/(-y1)\n",
    "    forward_slope = y3\n",
    "    backward_slope = -y1\n",
    "    plt.fill_between(xx,xx*forward_slope,xx*backward_slope,color='k',alpha=0.2,zorder=0)\n",
    "    for limiter in ('minmod', 'vanleer','superbee','MC'):\n",
    "        sigma = phi(np.array(theta),limiter)*(-y1)\n",
    "        ax.plot(xx,sigma*xx,alpha=0.5,lw=2)\n",
    "    ax.legend( ('','minmod', 'vanleer','superbee','MC'), loc='best' )\n",
    "    ax.plot(x,y,'ok',markersize=15,alpha=0.5)\n",
    "\n",
    "    #return fig\n",
    "\n",
    "interact(compare_limiters,y1=FloatSlider(min=-1., max=1., step=0.1, value=-0.3,description='$y_{i-1}$',labelcolor='k'),#,orientation='vertical'),\n",
    "                          y3=FloatSlider(min=-1., max=1., step=0.1, value=0.8,description='$y_{i+1}$'));#,orientation='vertical'));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The shaded region in the plot above shows the range of slopes that would give at least 2nd-order accuracy.  Play with the sliders and answer the following questions:\n",
    "\n",
    "- Which limiter usually chooses the flattest approximation?  Does it always?\n",
    "- Which limiter usually chooses the steepest approximation?  Does it always?\n",
    "- In which situations do all the limiters give the same slope?  Why?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Higher-order interpolation: WENO\n",
    "If we want to get higher-order accuracy (in smooth regions), then we have to give up the TVD property -- at least, in the sense defined above.  The most common approach for higher order non-oscillatory piecewise interpolation is known as weighted essentially non-oscillatory (WENO) interpolation.\n",
    "\n",
    "WENO is a very effective technique for interpolating or reconstructing functions that contain discontinuities without introducing oscillations.  We'll focus on 5th-order WENO interpolation, which is the most commonly used.\n",
    "\n",
    "Let's generate some function values to interpolate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note: uses PyWENO v. 0.11.2\n",
    "import sympy\n",
    "from pyweno import symbolic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mpld3             # Skip this cell if you don't have mpld3 installed\n",
    "mpld3.enable_notebook()  # or just go and do it now: pip install mpld3\n",
    "                         # This allows you to zoom and pan the plots.\n",
    "matplotlib.rcParams.update({'font.size': 18})\n",
    "colors = 'brg'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weno_order = 5  # must be odd\n",
    "k = (weno_order+1)/2\n",
    "\n",
    "size = (width,4)\n",
    "plt.figure(figsize=size)\n",
    "x=np.arange(-k+1,k)\n",
    "y=np.random.rand(len(x))\n",
    "#y = np.array((1.,1.,1.,0.,0.))\n",
    "plt.plot(x,y,'ok')\n",
    "plt.axis((-(k-.5),k-.5,-0.5,2.1));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def stencil_interpolant(x,y,n,offset):\n",
    "    \"\"\"Return the polynomial interpolant (of degree n-1) \n",
    "       through the points (x_j,y_j) for offset <= j <= offset+n-1.\n",
    "    \"\"\"\n",
    "    return np.poly1d(np.polyfit(x[offset:offset+n],y[offset:offset+n],n-1))\n",
    "\n",
    "def plot_interpolants(x,y,interpolants,axis=None,color='kbrg'):\n",
    "    if axis is None:\n",
    "        fig, axis = plt.subplots(figsize=size)\n",
    "    xc = np.linspace(-0.5,0.5)\n",
    "    xx = np.linspace(-(k-1),k-1)\n",
    "    for i, interpolant in enumerate(interpolants):\n",
    "        axis.plot(xx,interpolant(xx),'-'+color[i])\n",
    "        axis.plot(xc,interpolant(xc),'-'+color[i],linewidth=5,alpha=0.5)\n",
    "    axis.plot(x,y,'ok')\n",
    "    axis.axis((-(k-.5),k-.5,-0.5,2.1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ordinary polynomial interpolation yields an oscillatory polynomial that also exceeds the bounds of the data:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For application to hyperbolic conservation laws, our main interest is in getting values of the function at the half-integer points (interfaces).  Let's suppose we're trying to interpolate around $x=0$, at $x=\\pm 1/2$.  Instead of using all 5 points, we could just use three points, which might give us a less oscillatory interpolant, at least in that interval.  Using the 5 points we're given, there are three natural choices of interpolation stencil: the leftmost three, the middle three, or the rightmost three.  Let's see what each of these quadratic interpolants looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_opt = stencil_interpolant(x,y,5,0)\n",
    "plot_interpolants(x,y,[p_opt])\n",
    "plt.title('Quartic interpolant');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(3,1,figsize=(width,10))\n",
    "names = ['left','right','center']\n",
    "p = []\n",
    "for i in range(int(k)):\n",
    "    p.append(stencil_interpolant(x,y,k,i))\n",
    "    plot_interpolants(x,y,[p[i]],axis=ax[i],color=[colors[i]])\n",
    "    ax[i].set_title(names[i]+' interpolant')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are all three quadratic interpolants together with the quartic interpolant for comparison:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_interpolants(x,y,p+[p_opt],color='brgk')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The quadratic interpolants look less oscillatory, but they're also less accurate.  The WENO idea is to use the high-order interpolant (with all 5 points) if the data is smooth, but to use one of the lower-order interpolants (or a combination of them) if the data is not smooth.  This is achieved by computing point values of the interpolant as weighted averages of the point values of the candidate polynomials, e.g.\n",
    "\n",
    "$$y_{x_{i-1/2}} = w_{1,-1/2} p_\\text{left}(x_{i-1/2}) + w_{2,-1/2} p_\\text{center}(x_{i-1/2}) + w_{3,-1/2} p_\\text{right}(x_{i-1/2}).$$\n",
    "\n",
    "Of course, there is some particular set of weights that gives the quartic interpolant:\n",
    "\n",
    "$$y_{x_{i-1/2}} = \\gamma_{1,-1/2} p_\\text{left}(x_{i-1/2}) + \\gamma_{2,-1/2} p_\\text{center}(x_{i-1/2}) + \\gamma_{3,-1/2} p_\\text{right}(x_{i-1/2}).$$\n",
    "\n",
    "We will want to have $w_{j,-1/2} \\approx \\gamma_{j,-1/2}$ for smooth data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_opt_weights(k,xi):\n",
    "    \"\"\"\n",
    "    Get the optimal weights (gamma) at points xi.\n",
    "    \"\"\"\n",
    "    if not hasattr(xi,'__iter__'): xi = [xi]\n",
    "    opt_weights = symbolic.optimal_weights(k,xi)\n",
    "    gamma = {}\n",
    "    for i, xi_val in enumerate(xi):\n",
    "        gamma[xi_val] = np.empty(k)\n",
    "        for j in range(k):\n",
    "            gamma[xi_val][j] = opt_weights[0][(i,j)]\n",
    "\n",
    "    return gamma\n",
    "\n",
    "gamma = compute_opt_weights(k,(-1,0.5,1))\n",
    "\n",
    "print \"$\\gamma_{j,-1/2}$:\",  gamma[-1]\n",
    "print \"$\\gamma_{j,+1/2}$:\",  gamma[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How does one determine if a polynomial is non-oscillatory?  There are several ways proposed in the literature, but the original and most widely used is the weighted Sobolev norm:\n",
    "\n",
    "$$\\beta = \\sum_{l=1}^k \\Delta x^{2l-1} \\int_{x_{i-1/2}}^{x_{i+1/2}} \\left(\\frac{d^l}{dx^l}p(x)\\right)^2 dx.$$\n",
    "\n",
    "Put simply, $\\beta$ is a scaled sum of the square $L^2$ norms of all the derivatives of the polynomial over the interval where it will be used.  The scaling is chosen to make the \"smoothness indicator\" $\\beta$ independent of the choice of $\\Delta x$ (note that $\\Delta x = 1$ in our example data).\n",
    "\n",
    "As each of the interpolants above is a linear function of the values $y_i$, the smoothness indicators are quadratic functions of the $y_i$ and can be expressed in the generic form\n",
    "\n",
    "$$\\beta = \\sum_{m=-2}^{2} \\sum_{n=-2}^{m}\n",
    "    C_{m,n} y_{i-k+m} y_{i-k+n}$$\n",
    "\n",
    "Of course, the coefficients $C_{m,n}$ will be different for each of the candidate polynomials $p_\\text{left},p_\\text{center},p_\\text{right}$.  We can use the Python package PyWeno to automatically compute these coefficients and then apply them to our data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_smoothness_indicators(y,k):\n",
    "    C = symbolic.jiang_shu_smoothness_coefficients(k)\n",
    "    beta = np.zeros((k,1))\n",
    "    for m in range(k):\n",
    "        for n in range(m+1):\n",
    "            for r in range(len(beta)):\n",
    "                beta[r] = beta[r] + C[(r,n,m)] * y[r+m] * y[r+n]\n",
    "    return beta\n",
    "\n",
    "beta = compute_smoothness_indicators(y,k)\n",
    "print beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we use these smoothness indicators to determine a weighting for the candidate polynomials.  Observe that a large smoothness indicator means a polynomial has large derivatives, so we will want to give it less weight (perhaps they should be called non-smoothness indicators). \n",
    "\n",
    "$$\\tilde{w}_j = \\frac{\\gamma_j}{(\\epsilon + \\beta_j)^2}$$\n",
    "\n",
    "Here $\\epsilon$ is a small number used to avoid division by zero.  We also normalize the weights so that they sum to unity:\n",
    "\n",
    "$$w_j = \\frac{\\tilde{w}_j}{\\sum_j\\tilde{w}_j}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_weights(gamma, beta, epsilon=1.e-6):\n",
    "    k = len(beta)\n",
    "    w = np.empty(k)\n",
    "    for j in range(k):\n",
    "        w[j] = gamma[j]/(epsilon+beta[j])**2\n",
    "        \n",
    "    wsum = np.sum(w)\n",
    "    return w/wsum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = {}\n",
    "for xi in (-1,1):\n",
    "    w = compute_weights(gamma[xi],beta)\n",
    "    q[xi] = w[0]*p[0](xi/2.) + w[1]*p[1](xi/2.) + w[2]*p[2](xi/2.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the final reconstructed values given by WENO (indicated by the large grey circles):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_interpolants(x,y,p+[p_opt],color=['b','r','g','k'])\n",
    "plt.hold(True)\n",
    "plt.plot(-0.5,q[-1],'ok',alpha=0.3,markersize=15)\n",
    "plt.plot(0.5,q[1],'ok',alpha=0.3,markersize=15)\n",
    "\n",
    "plt.axis((-(k-.5),k-.5,-0.5,2.1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's some code to plot everything for some given $(x,y)$ values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "styles = { 'left' : 'b', 'center' : 'r', 'right' : 'g'}\n",
    "size = (16,4); fs = 20\n",
    "def WENO_visualization(x,y,xi=(-1,1)):\n",
    "    \"\"\"\n",
    "        (x,y): data to interpolate\n",
    "        xi: points at which to evaluate interpolant (w.r.t. (-1,1) reference interval)\n",
    "    \"\"\"\n",
    "    xx = np.linspace(np.min(x),np.max(x))\n",
    "    color=['b','r','g']\n",
    "    plt.figure(figsize=size)\n",
    "    plt.hold(True)\n",
    "    ax1 = plt.subplot2grid((1,8), (0,0), colspan=6)\n",
    "    ax2 = plt.subplot2grid((1,8), (0,6))\n",
    "    ax3 = plt.subplot2grid((1,8), (0,7))\n",
    "\n",
    "    K = len(y)\n",
    "    k = (K+1)/2\n",
    "    assert len(x)==K\n",
    "    p_opt=np.poly1d(np.polyfit(x,y,K-1))\n",
    "    p = {}\n",
    "    for name, offset in zip(('left','right','center'),range(k)):\n",
    "        p[name] = stencil_interpolant(x,y,k,offset)\n",
    "        \n",
    "    gamma = compute_opt_weights(k,xi)\n",
    "    beta = compute_smoothness_indicators(y,k)\n",
    "    \n",
    "    w = {}; q = {}\n",
    "    for loc in xi:\n",
    "        w[loc] = compute_weights(gamma[loc],beta)\n",
    "        q[loc] = w[loc][0]*p['left'](loc/2.) \\\n",
    "               + w[loc][1]*p['center'](loc/2.) \\\n",
    "               + w[loc][2]*p['right'](loc/2.)\n",
    "    ax2.bar(range(3),w[-1],color=color,align='center'); \n",
    "    ax2.set_title(r'$w_{i-1/2}$',fontsize=fs)\n",
    "    ax3.bar(range(3),w[1],color=color,align='center')\n",
    "    ax3.set_title(r'$w_{i+1/2}$',fontsize=fs)\n",
    "    for ax in (ax2,ax3):\n",
    "        ax.set_xticks(range(3)); \n",
    "        ax.set_xticklabels(('left','center','right'))\n",
    "        ax.set_ylim(0,1); ax.set_yticks((0,1))\n",
    "\n",
    "    for name, interpolant in p.iteritems():\n",
    "        ax1.plot(xx,interpolant(xx),styles[name])\n",
    "        xc = np.linspace(-0.5,0.5)\n",
    "        ax1.plot(xc,interpolant(xc),styles[name],linewidth=5)\n",
    "        ax1.plot(x,y,'ok')\n",
    "    ax1.hold(True)\n",
    "    ax1.plot(xx,p_opt(xx),'-k',x,y,'ok',linewidth=2)\n",
    "    for loc in xi:\n",
    "        ax1.plot(loc/2.,q[loc],'ok', alpha=0.3,markersize=15)\n",
    "        ax1.plot(loc/2.,q[loc],'ok',alpha=0.3,markersize=15)\n",
    "        \n",
    "    ax1.axis((-(k-0.8),k-0.8,-0.5,2.1));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "y=np.random.rand(len(x))\n",
    "WENO_visualization(x,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The bar charts on the right show the relative weight given to each of the quadratic interpolants when computing the left and right interpolated values.\n",
    "Try running the box above a few times, or insert your own $y$ values.  What happens if you use a step function for $y$?  Let's see:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y=np.array( (1,1,1,0,0) )\n",
    "WENO_visualization(x,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a perfect step function, WENO picks the flat interpolant, just like any TVD limiter would!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison of several limiters for advection\n",
    "In practice, all of these limiters are used as components of numerical solvers for hyperbolic PDEs.  The simplest hyperbolic PDE is the advection equation:\n",
    "\n",
    "$$ q_t + a q_x = 0.$$\n",
    "\n",
    "The solution $q$ simply translates at velocity $a$; if you're not familiar with this, take a look at my [HyperPython lesson on advection](http://nbviewer.ipython.org/github/ketch/HyperPython/blob/master/Lesson_01_Advection.ipynb) and then the [lesson on high-resolution methods](http://nbviewer.ipython.org/github/ketch/HyperPython/blob/master/Lesson_03_High-resolution_methods.ipynb).\n",
    "\n",
    "The cells below solve the advection equation using several of the limiters we've discussed.  To run this part, you need to [install PyClaw](http://www.clawpack.org/installing.html) and Visclaw, which can be most easily accomplished via\n",
    "\n",
    "    pip install clawpack"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from clawpack import pyclaw\n",
    "from clawpack import riemann\n",
    "import matplotlib\n",
    "from matplotlib import animation\n",
    "from clawpack.visclaw.JSAnimation import IPython_display\n",
    "\n",
    "def setup(scheme='minmod',cfl_max=0.9,IC='gauss_square',mx=100):\n",
    "    if 'weno' in scheme:\n",
    "        solver = pyclaw.SharpClawSolver1D(riemann.advection_1D)\n",
    "    else:\n",
    "        solver = pyclaw.ClawSolver1D(riemann.advection_1D)\n",
    "\n",
    "    solver.bc_lower[0] = pyclaw.BC.periodic\n",
    "    solver.bc_upper[0] = pyclaw.BC.periodic\n",
    "    \n",
    "    if scheme in ('minmod','superbee','MC','vanleer'):\n",
    "        solver.limiters = getattr(pyclaw.limiters.tvd,scheme)\n",
    "    #elif scheme == 'CT':\n",
    "        #solver.limiters = pyclaw.limiters.tvd.cada_torrilhon_limiter\n",
    "    elif scheme == 'Lax-Wendroff':\n",
    "        solver.limiters = 0\n",
    "    elif scheme == 'first-order':\n",
    "        solver.order = 1\n",
    "    elif 'weno' in scheme:\n",
    "        solver.weno_order = int(scheme[4:]) #weno5, weno7, ...\n",
    "    else:\n",
    "        raise Exception('Unrecognized limiter')\n",
    "\n",
    "    solver.cfl_max = cfl_max\n",
    "    solver.cfl_desired = cfl_max*0.9\n",
    "\n",
    "    x = pyclaw.Dimension(0.0,1.0,mx)\n",
    "    domain = pyclaw.Domain(x)\n",
    "    num_eqn = 1\n",
    "    state = pyclaw.State(domain,num_eqn)\n",
    "    state.problem_data['u']=1.\n",
    "\n",
    "    grid = state.grid\n",
    "    xc = grid.x.centers\n",
    "    if IC=='gauss_square':\n",
    "        beta=200.; x0=0.3\n",
    "        state.q[0,:] = np.exp(-beta * (xc-x0)**2) + (xc>0.6)*(xc<0.8)\n",
    "    elif IC=='wavepacket':\n",
    "        beta=100.; x0=0.5\n",
    "        state.q[0,:] = np.exp(-beta * (xc-x0)**2) * np.sin(80.*xc)\n",
    "    else:\n",
    "        raise Exception('Unrecognized initial condition.')\n",
    "\n",
    "    claw = pyclaw.Controller()\n",
    "    claw.solution = pyclaw.Solution(state,domain)\n",
    "    claw.solver = solver\n",
    "    claw.keep_copy = True\n",
    "    claw.output_format = None\n",
    "\n",
    "    claw.tfinal =10.0\n",
    "    return claw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This cell may take a few seconds to run\n",
    "results = []\n",
    "schemes = ('first-order','Lax-Wendroff','minmod','superbee','MC','vanleer','weno5','weno7','weno9')\n",
    "for scheme in schemes:\n",
    "    claw = setup(scheme=scheme)\n",
    "    claw.verbosity = 0\n",
    "    claw.run()\n",
    "    results.append(claw.frames)\n",
    "    \n",
    "def animate(results,ymin=-0.1):\n",
    "    fig = plt.figure(figsize=(width,8))\n",
    "\n",
    "    N = len(results)\n",
    "    n = int(np.ceil(np.sqrt(N)))\n",
    "    axes = []\n",
    "    gs1 = matplotlib.gridspec.GridSpec(n, n)\n",
    "    gs1.update(wspace=0.,hspace=0.)\n",
    "\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            k = n*i + j\n",
    "            if k<N:\n",
    "                axes.append(plt.subplot(gs1[i,j]));\n",
    "                if i<n-1:\n",
    "                    axes[-1].xaxis.set_ticklabels(())\n",
    "                if j>0:\n",
    "                    axes[-1].yaxis.set_ticklabels(())\n",
    "\n",
    "    lines = [0]*len(schemes)\n",
    "    for i in range(len(lines)):\n",
    "        lines[i], = axes[i].plot([], [], lw=2)\n",
    "\n",
    "    xc = results[0][0].p_centers[0]\n",
    "\n",
    "    for i,ax in enumerate(axes):\n",
    "        ax.set_xlim(0,1); ax.set_ylim(ymin,1.3)\n",
    "        #ax.grid()\n",
    "        ax.set_title(schemes[i], x = 0.5, y=0.85 )\n",
    "        ax.plot(xc,results[i][0].q[0,:],color='k',alpha=0.3)\n",
    "    \n",
    "\n",
    "    def fplot(frame_number):\n",
    "        fig.suptitle('Solution after %s cycles' % frame_number, fontsize=20)\n",
    "        for i, line in enumerate(lines):\n",
    "            line.set_data(xc,results[i][frame_number].q[0,:])\n",
    "        return lines,\n",
    "\n",
    "    return matplotlib.animation.FuncAnimation(fig, fplot, frames=len(claw.frames), interval=30)\n",
    "\n",
    "animate(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the plot above, the solution advects across the full domain once between each frame of the animation (the boundary is periodic).  By stepping through the animation, you can see how each limiter modifies the shape of the solution over time.  The Lax-Wendroff method is based on a centered-difference approximation with no limiting; notice that it creates oscillations and is also less accurate than the limiter-based methods.  For the advection equation, oscillations and overshoots are not a serious problem, but in the context of fluid dynamics or water wave simulations, they can be catastrophic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "3.9",
   "language": "python",
   "name": "3.9"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.10"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}