Burgers_approximate.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# An approximate solver for Burgers' equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a first example, we return to the inviscid Burgers' equation that we studied in [Burgers](Burgers.ipynb):  \n",
    "\\begin{align} \\label{burgers}\n",
    "q_t + \\left(\\frac{1}{2}q^2\\right)_x & = 0.\n",
    "\\end{align}  \n",
    "Although it is easy to solve the Riemann problem for (\\ref{burgers}) exactly, it is nevertheless interesting to consider approximate solvers because a numerical scheme does not make use of the full Riemann solution.  Furthermore, Burgers' equation provides a simple setting in which some of the intricacies of more complicated approximate solvers can be introduced.  Recall that we are interested in **approximate solutions that consist entirely of traveling discontinuities**.\n",
    "\n",
    "Recall from [Approximate_solvers](Approximate_solvers.ipynb) that we use $Q^*$ to denote the state along $x/t = 0$ in the approximate Riemann solution, and $F^*$ for the corresponding interface flux that might be used in a finite volume method. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Shock wave solutions\n",
    "\n",
    "Recall that the exact Riemann solution for (\\ref{burgers}) consists of a single shock or rarefaction wave.  We have a shock if $q_\\ell>q_r$ and we have a rarefaction if $q_\\ell < q_r$.  In the case of a shock wave, we can simply use the exact solution as our approximation.  We have a single wave of strength ${\\cal W} = q_r-q_\\ell$ traveling at speed $s=(q_r+q_\\ell)/2$ (since there is only one wave we drop the subscript on ${\\cal W_1}$ and $s_1$).  \n",
    "\n",
    "In this case, the numerical flux is $F^*=f(q_\\ell)$ if $s\\geq 0$ and $F^*=f(q_r)$ if $s\\leq 0$.  In the special case $s=0$ we have a stationary shock, and it must be that $f(q_\\ell)=f(q_r)$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rarefaction wave solutions\n",
    "\n",
    "As discussed in [Approximate_solvers](Approximate_solvers.ipynb), for numerical purposes it is convenient to approximate a rarefaction wave by a traveling discontinuity.  For Burgers' equation this may seem unnecessary, but for more complicated solvers for systems of equations it will be essential.\n",
    "\n",
    "We will approximate the effect of the rarefaction wave by a fictitious shock\n",
    "$${\\cal W} = q_r-q_\\ell,$$  \n",
    "whose speed is given by the Rankine-Hugoniot jump condition:  \n",
    "$$s = \\frac{f(q_r)-f(q_\\ell)}{q_r-q_\\ell} = \\frac{q_r + q_\\ell}{2}.$$  \n",
    "Recall that this is indeed a weak solution of the Riemann problem.  This fictitious shock is not entropy-satisfying, but as long as it's traveling entirely to the left or entirely to the right, the effect on the numerical solution will be the same as if we used a (entropy-satisfying) rarefaction wave.  This must be true since both solutions are conservative and so averaging the wave over the single cell it affects must give the same value.  The numerical flux is again $F^*=f(q_\\ell)$ if $s\\geq 0$ and $F^*=f(q_r)$ if $s \\leq 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because this is a scalar equation with convex flux, both the Roe and HLL approaches will simplify to what we have already described.  But we briefly discuss them here to further illustrate the main ideas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A Roe solver"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we consider a linearized solver, in which we replace our nonlinear hyperbolic system with a linearization about some intermediate state $\\hat{q}$.  For Burgers' equation, the quasilinear form is $q_t + q q_x = 0$ and the linearization gives  \n",
    "$$q_t + \\hat{q}q_x = 0.$$  \n",
    "This is simply the advection equation with velocity $\\hat{q}$.  The solution of the Riemann problem for this equation consists of a wave ${\\cal W} = q_r - q_\\ell$ traveling at speed $s = \\hat{q}$.  It remains only to determine the state $\\hat{q}$, which is also the wave speed.\n",
    "\n",
    "The defining feature of a Roe linearization is that it should give the exact solution in case the states $(q_\\ell, q_r)$ lie on a single Hugoniot locus -- i.e., when the solution is a single shock.  We can achieve this by choosing  \n",
    "$$\\hat{q} = \\frac{q_r + q_\\ell}{2}$$  \n",
    "so that the Rankine-Hugoniot condition is satisfied.  This is equivalent to using the approximate solver described already in the sections above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Examples\n",
    "Below we show solutions for three examples; the first involves a shock, the second a (non-transonic) rarefaction, and the third a transonic rarefaction.  In the first row we plot the exact solution in terms of the waves in the $x$-$t$ plane; in the second row we plot numerical solutions obtained with PyClaw by using a simple first-order method combined with the Riemann solver.  \n",
    "\n",
    "To examine the Python code for this chapter, see:\n",
    "\n",
    " - [exact_solvers/burgers.py](exact_solvers/burgers.py) ...\n",
    "   [on github.](https://github.com/clawpack/riemann_book/blob/FA16/exact_solvers/burgers.py)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "%config InlineBackend.figure_format = 'svg'\n",
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import interact\n",
    "from ipywidgets import widgets\n",
    "from ipywidgets import IntSlider\n",
    "from exact_solvers import burgers\n",
    "from utils import riemann_tools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def setup(q_l, q_r, N=500, efix=False):\n",
    "\n",
    "    from clawpack import pyclaw\n",
    "    from clawpack import riemann\n",
    "    \n",
    "    rs = riemann.burgers_1D_py.burgers_1D\n",
    "    solver = pyclaw.ClawSolver1D(rs)        \n",
    "    solver.order = 1\n",
    "    solver.kernel_language = 'Python'\n",
    "    \n",
    "    solver.bc_lower[0]=pyclaw.BC.extrap\n",
    "    solver.bc_upper[0]=pyclaw.BC.extrap\n",
    "\n",
    "    x = pyclaw.Dimension(-1.0,1.0,N,name='x')\n",
    "    domain = pyclaw.Domain([x])\n",
    "    state = pyclaw.State(domain,1)\n",
    "    \n",
    "    state.problem_data['efix'] = efix\n",
    "\n",
    "    xc = state.grid.p_centers[0]\n",
    "\n",
    "    state.q[0 ,:] = (xc<=0)*q_l + (xc>0)*q_r\n",
    "\n",
    "    claw = pyclaw.Controller()\n",
    "    claw.tfinal = 0.5\n",
    "    claw.solution = pyclaw.Solution(state,domain)\n",
    "    claw.solver = solver\n",
    "    claw.keep_copy = True\n",
    "    claw.verbosity=0\n",
    "\n",
    "    return claw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 200  # number of grid cells to use\n",
    "\n",
    "shock = setup(2.,1.,N=N)\n",
    "shock.run()\n",
    "shocksol = burgers.exact_riemann_solution(2.,1.)\n",
    "\n",
    "raref = setup(1.,2.,N=N)\n",
    "raref.run()\n",
    "rarefsol = burgers.exact_riemann_solution(1.,2.)\n",
    "\n",
    "transonic = setup(-1.,2.,N=N)\n",
    "transonic.run()\n",
    "transonicsol = burgers.exact_riemann_solution(-1.,2.)\n",
    "\n",
    "def plot_frame(i):\n",
    "    fig, axes = plt.subplots(2,3,figsize=(8,4))\n",
    "    for ax in axes[0]:\n",
    "        ax.set_xlim((-1,1)); ax.set_ylim((-1.1,2.1))\n",
    "    \n",
    "    sxc = shock.frames[0].grid.x.centers\n",
    "    rxc = raref.frames[0].grid.x.centers\n",
    "    txc = transonic.frames[0].grid.x.centers\n",
    "    \n",
    "    t = i/10.\n",
    "    axes[1][0].plot(sxc, shock.frames[i].q[0,:],'-k',lw=2)\n",
    "    axes[1][0].set_title('Shock(t=' + str(t) + ')')\n",
    "    axes[1][1].plot(rxc, raref.frames[i].q[0,:],'-k',lw=2)\n",
    "    axes[1][1].set_title('Rarefaction(t=' + str(t) + ')')\n",
    "    axes[1][2].plot(txc, transonic.frames[i].q[0,:],'-k',lw=2)\n",
    "    axes[1][2].set_title('Transonic rarefaction(t=' + str(t) + ')')\n",
    "    \n",
    "    riemann_tools.plot_waves(*shocksol,ax=axes[0][0],t=t)\n",
    "    riemann_tools.plot_waves(*rarefsol,ax=axes[0][1],t=t)\n",
    "    riemann_tools.plot_waves(*transonicsol,ax=axes[0][2],t=t)\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "    \n",
    "interact(plot_frame, i=IntSlider(min=0, max=10, value=5, description='Frame'));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solutions obtained for the shock wave and for the first rarefaction wave are good approximations of the true solution.  In the case of the transonic rarefaction, however, we see that part of the rarefaction has been replaced by an entropy-violating shock.  Refining the grid in this numerical solution would not help; the method converges to a weak solution that does not satisfy the entropy condition.  *In the live notebook you can try adjusting $N$ in the cell above to investigate this.*\n",
    "\n",
    "At the end of this chapter we will show how to apply an *entropy fix* so that the solver gives a good approximation also in the transonic case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Two-wave solvers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For Burgers' equation, the exact Riemann solution consists of a single wave.  It is thus natural to modify the HLL approach by assuming that one of the waves vanishes, and denote the speed of the other wave simply by $s$.  Then the conservation condition discussed in [Approximate_solvers](Approximate_solvers.ipynb#Two-wave-solvers) reduces to  \n",
    "$$f(q_r) - f(q_\\ell) = s (q_r - q_\\ell),$$  \n",
    "which is just the Rankine-Hugoniot condition and again leads to the speed discussed above.  Since the solution involves only one wave, that wave must carry the entire jump $q_r - q_\\ell$, so this solver is entirely equivalent to that already described.\n",
    "\n",
    "It is also possible to follow the full HLL approach, taking\n",
    "\\begin{align*}\n",
    "s_1 & = \\min f'(q) = \\min(q_\\ell, q_r) \\\\\n",
    "s_2 & = \\max f'(q) = \\max(q_\\ell, q_r).\n",
    "\\end{align*}\n",
    "Since $f(q) = q^2/2$, it is easy to verify that regardless of the values of $q_\\ell$ and $q_r$, this leads to  \n",
    "$$q_m = \\frac{q_r + q_\\ell}{2},$$  \n",
    "so that each of the two waves carries half of the jump."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Transonic rarefactions\n",
    "\n",
    "In the approaches above, the solution was approximated by a single wave traveling either to the left or right.  For this scalar problem, this \"approximation\" is, in fact, an exact weak solution of the Riemann problem.  As discussed already, we do not typically need to worry about the fact that it may be entropy-violating, since the effect on the numerical solution (after averaging) is identical to that of the entropy-satisfying solution.\n",
    "\n",
    "However, if $q_\\ell < 0 < q_r$, then the true solution is a *transonic rarefaction*, in which part of the wave travels to the left and part travels to the right.  In this case, the true Riemann solution would lead to changes to both the left and right adjacent cells, whereas an approximate solution with a single wave will only modify one or the other.  This leads to an incorrect numerical solution (on the macroscopic level), as seen in the PyClaw example above.  It is therefore necessary to modify the Riemann solver, imposing what is known as an *entropy fix* in the transonic case.\n",
    "\n",
    "Specifically, we use a solution consisting of two waves, each of which captures the net effect of the corresponding rarefaction wave that appears in the exact solution:\n",
    "\n",
    "\\begin{align}\n",
    "{\\cal W}_1 & = q_m - q_\\ell, & s_1 = \\frac{q_\\ell + q_m}{2} \\\\\n",
    "{\\cal W}_2 & = q_r - q_m, & s_2 = \\frac{q_m + q_r}{2}.\n",
    "\\end{align}\n",
    "\n",
    "For Burgers' equation, the value $q_s=0$ is the *sonic point* satisfying $f'(q_s)=0$. A transonic rarefaction wave takes the value $Q^*=q_s$ along $x/t = 0$ and so it makes sense to choose $q_m = 0$.  The formulas above then simplify to\n",
    "\n",
    "\\begin{align}\n",
    "{\\cal W}_1 & = - q_\\ell, & s_1 = q_\\ell / 2 < 0 \\\\\n",
    "{\\cal W}_2 & = q_r, & s_2 = q_r / 2 > 0.\n",
    "\\end{align}\n",
    "\n",
    "Note that this can also be viewed as an HLL solver, although not with the usual choice of wave speeds.  Choosing the waves speeds $s_1=q_\\ell/2$ and $s_2=q_r/2$ and then solving for $q_m$ by requiring conservation gives $q_m=0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now repeat the PyClaw simulation performed above, but with the entropy fix applied."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 200  # number of grid cells to use\n",
    "shock_efix = setup(2.,1.,N=N,efix=True)\n",
    "shock_efix.run()\n",
    "shocksol = burgers.exact_riemann_solution(2.,1.)\n",
    "\n",
    "raref_efix = setup(1.,2.,N=N,efix=True)\n",
    "raref_efix.run()\n",
    "rarefsol = burgers.exact_riemann_solution(1.,2.)\n",
    "\n",
    "transonic_efix = setup(-1.,2.,N=N,efix=True)\n",
    "transonic_efix.run()\n",
    "transonicsol = burgers.exact_riemann_solution(-1.,2.)\n",
    "\n",
    "def plot_frame(i):\n",
    "    fig, axes = plt.subplots(2,3,figsize=(8,4))\n",
    "    for ax in axes[0]:\n",
    "        ax.set_xlim((-1,1)); ax.set_ylim((-1.1,2.1))\n",
    "        \n",
    "    sxc = shock_efix.frames[0].grid.x.centers\n",
    "    rxc = raref_efix.frames[0].grid.x.centers\n",
    "    txc = transonic_efix.frames[0].grid.x.centers\n",
    "    \n",
    "    t = i/10.\n",
    "    axes[1][0].plot(sxc, shock_efix.frames[i].q[0,:],'-k',lw=2)\n",
    "    axes[1][0].set_title('Shock(t=' + str(t) + ')')\n",
    "    axes[1][1].plot(rxc, raref_efix.frames[i].q[0,:],'-k',lw=2)\n",
    "    axes[1][1].set_title('Rarefaction(t=' + str(t) + ')')\n",
    "    axes[1][2].plot(txc, transonic_efix.frames[i].q[0,:],'-k',lw=2)\n",
    "    axes[1][2].set_title('Transonic rarefaction(t=' + str(t) + ')')\n",
    "    \n",
    "    riemann_tools.plot_waves(*shocksol,ax=axes[0][0],t=t)\n",
    "    riemann_tools.plot_waves(*rarefsol,ax=axes[0][1],t=t)\n",
    "    riemann_tools.plot_waves(*transonicsol,ax=axes[0][2],t=t)\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "    \n",
    "interact(plot_frame, i=IntSlider(min=0, max=10, value=5, description='Frame'));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The entropy fix has no effect on the first two solutions, since it is applied only in the case of a transonic rarefaction.  The third solution is greatly improved over the previous result.  There is still a small *entropy glitch* visible at the origin, but the numerical solution will converge to the correct weak solution as the grid is refined.  *In the live notebook you can try adjusting $N$ in the cell above to investigate this.*"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.6.5"
  },
  "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": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}