Lesson_05_Approximate_Solvers.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Philosopher:400,700,400italic,700italic' rel='stylesheet' type='text/css'>\n",
       "\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: #888;\n",
       "    color: #f6f6f6;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff;\n",
       "    color: #333;\n",
       "    width: 1200px;\n",
       "    margin: auto;\n",
       "    padding-left: 1em;\n",
       "    padding-right: 1em;\n",
       "    padding-top: 2ex;\n",
       "    overflow-x: auto;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "margin-top:0.8em;\n",
       "}\n",
       "\n",
       "div.cell { /* set cell width to about 80 chars */\n",
       "    width: 900px;\n",
       "}\n",
       "\n",
       "div.cell.border-box-sizing.code_cell.running { \n",
       "    /* draw border around running cells */\n",
       "    border: 3px solid #111;\n",
       "}\n",
       "\n",
       "\n",
       "h1 {\n",
       "    font-family: 'Philosopher', sans-serif;\n",
       "}\n",
       "h2 {\n",
       "    font-family: 'Philosopher', serif;\n",
       "}\n",
       "h3{\n",
       "    font-family: 'Philosopher', serif;\n",
       "    font-style: 'italic';\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "}\n",
       "h4{\n",
       "    font-family: 'Philosopher', serif;\n",
       "}\n",
       "h5 {\n",
       "    font-family: 'Alegreya Sans', sans-serif;\n",
       "}\t\n",
       "h6 {\n",
       "    font-family: 'PT Mono', sans-serif;\n",
       "}\t\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Arvo' sans-serif;\n",
       "    line-height: 130%;\n",
       "    font-size: 115%;\n",
       "    width:800px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       ".CodeMirror{\n",
       "        font-family: \"PT Mono\";\n",
       "                    font-size: 100%;\n",
       "            }\n",
       "\n",
       ".text_cell_render h1 {\n",
       "    font-weight: 400;\n",
       "    font-size: 64pt;\n",
       "    line-height: 100%;\n",
       "    color: rgb(12,85,97);\n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-weight: 700;\n",
       "    font-size: 24pt;\n",
       "    line-height: 100%;\n",
       "    color: rgb(171,165,131);\n",
       "    margin-bottom: 0.1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-style: italic;\n",
       "    color: rgb(95,92,72);\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {\n",
       "    font-weight: 300;\n",
       "    font-size: 16pt;\n",
       "    color: grey;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .1em;\n",
       "    margin-top: 0.1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 {\n",
       "    font-weight: 300;\n",
       "    font-size: 10pt;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "div.cell.code_cell {  /* Group area containing both code and output */\n",
       "background-color: rgba(171,165,131,0.3); \n",
       "border-radius: 10px; /* rounded borders */\n",
       "padding: 1em;\n",
       "}\n",
       "\n",
       "</style>\n",
       "\n",
       "<script type=\"text/x-mathjax-config\">\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"],\n",
       "                           Macros: {N : \"\\\\mathbb{N}\"}\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {margin: \"40em 40em 40em 40em\",line-height : '200%', border: '1px solid ##CC0000' } },\n",
       "                    linebreaks: {'automatic': 'true'}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Approximate solutions to the Riemann Problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solutions in practice"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solutions to the Riemann problem are mainly used in two contexts:\n",
    "\n",
    "1. As reference solutions against which a numerical method is benchmarked, or\n",
    "2. As part of a numerical method, such as a high resolution shock capturing method, where the flux between two numerical cells is required."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the first case, accuracy is paramount and the complete solution (all wave speeds, and all intermediate states) is required. In the second case only one thing is required: the flux ${\\bf f}^*$ between the cells, which is the flux on the characteristic line $\\xi = x / t = 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this second case, the numerical method will have to repeatedly solve the Riemann problem. In a general problem, the solution may be needed tens of times *per cell, per timestep*, leading to millions (or more!) solutions in a simulation. The speed of the solution is then extremely important, and approximate solutions are often used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Roe-type solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The most obvious simplification is to reduce the nonlinear problem\n",
    "\n",
    "\\begin{equation}\n",
    "  \\partial_t {\\bf q} + \\partial_x {\\bf f}({\\bf q}) = {\\bf 0}\n",
    "\\end{equation}\n",
    "\n",
    "to the *linear* problem\n",
    "\n",
    "\\begin{equation}\n",
    "  \\partial_t {\\bf q} + A \\partial_x {\\bf q} = {\\bf 0},\n",
    "\\end{equation}\n",
    "\n",
    "where $A$ is a *constant* matrix that approximates the Jacobian $\\partial {\\bf f} / \\partial {\\bf q}$. We can then solve the linear problem exactly (e.g. by diagonalising the matrix and solving the resulting uncoupled advection equations), to find\n",
    "\n",
    "\\begin{align}\n",
    "  {\\bf q}(x, t) & = {\\bf q}_l + \\sum_{p: \\lambda^{(p)} < \\tfrac{x}{t}} \\left\\{ {\\bf l}^{(p)} \\cdot \\left( {\\bf q}_r - {\\bf q}_l \\right) \\right\\} {\\bf r}^{(p)}, \\\\\n",
    "   & = {\\bf q}_r - \\sum_{p: \\lambda^{(p)} > \\tfrac{x}{t}} \\left\\{ {\\bf l}^{(p)} \\cdot \\left( {\\bf q}_r - {\\bf q}_l \\right) \\right\\} {\\bf r}^{(p)}, \\\\\n",
    "   & = \\frac{1}{2} \\left( {\\bf q}_l + {\\bf q}_r \\right) + \\sum_{p: \\lambda^{(p)} < \\tfrac{x}{t}} \\left\\{ {\\bf l}^{(p)} \\cdot \\left( {\\bf q}_r - {\\bf q}_l \\right) \\right\\} {\\bf r}^{(p)} - \\sum_{p: \\lambda^{(p)} > \\tfrac{x}{t}} \\left\\{ {\\bf l}^{(p)} \\cdot \\left( {\\bf q}_r - {\\bf q}_l \\right) \\right\\} {\\bf r}^{(p)}.\n",
    "\\end{align}\n",
    "\n",
    "where $\\lambda^{(p)}, {\\bf r}^{(p)},$ and ${\\bf l}^{(p)}$ are the eigenvalues and the (right and left respectively) eigenvectors of $A$, ordered such that $\\lambda^{(1)} \\le \\dots \\le \\lambda^{(N)}$ as usual. All three solutions are equivalent; the last is typically used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given this complete solution, it is easily evaluated along $x = 0$, and the flux calculated from the result."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An even greater shortcut can be found by noting that we are approximating ${\\bf f} = A {\\bf q}$. Therefore the standard form is to write\n",
    "\n",
    "\\begin{equation}\n",
    "  {\\bf f}^* = \\frac{1}{2} \\left( {\\bf f}_l + {\\bf f}_r \\right) + \\sum_{p} \\left| \\lambda^{(p)} \\right| \\left\\{ {\\bf l}^{(p)} \\cdot \\left( {\\bf q}_r - {\\bf q}_l \\right) \\right\\} {\\bf r}^{(p)},\n",
    "\\end{equation}\n",
    "\n",
    "where now we are summing over all eigenvalues and eigenvectors. It should be noted that ${\\bf f}^* \\ne {\\bf f}({\\bf q}^*)$ in general, as the calculation of ${\\bf f}^*$ relied on an approximation to the flux."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to complete this specification of the solver, we only need to say how $A$ is defined. Roe gave the suggestion that\n",
    "\n",
    "\\begin{equation}\n",
    "  A = A({\\bf q}_{\\textrm{Roe}}) = \\left. \\frac{\\partial {\\bf f}}{\\partial {\\bf q}} \\right|_{{\\bf q}_{\\textrm{Roe}}},\n",
    "\\end{equation}\n",
    "\n",
    "where the *Roe average* ${\\bf q}_{\\textrm{Roe}}$ satisfies\n",
    "\n",
    "1. $A({\\bf q}_{\\textrm{Roe}}) \\left( {\\bf q}_r - {\\bf q}_l \\right) = {\\bf f}_r - {\\bf f}_l$,\n",
    "2. $A({\\bf q}_{\\textrm{Roe}})$ is diagonalizable with real eigenvalues, and\n",
    "3. $A({\\bf q}_{\\textrm{Roe}}) \\to \\partial {\\bf f} / \\partial {\\bf q}$ smoothly as ${\\bf q}_{\\textrm{Roe}} \\to {\\bf q}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is *possible* to construct the Roe average for many systems (such as the Euler equations, and the relativistic Euler equations). However, a simple arithmetic average is often nearly as good - in the sense that the algorithm will fail only slightly more often than the algorithm with the full Roe average!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The problem with Roe type solvers is that it approximates all waves as discontinuities. This leads to inaccuracies near rarefactions, and these can be catastrophically bad when the rarefaction fan crosses $\\xi = 0$ (a *sonic rarefaction*). It is possible to detect when these problems will occur (e.g. by looking at when $\\lambda^{(p)}$ changes sign between the left and right states) and change the approximation at this point, often known as an *entropy fix*. More systematic and complex methods that extend the Roe approach whilst avoiding this problem include the *Marquina* solver."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## HLL-type solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An alternative type of method simplifies the wave structure even more, by simplifying the number of waves. HLL (for Harten, Lax and van Leer) type solutions assume that\n",
    "\n",
    "1. there are two waves, both discontinuities, separating a constant central state in the solution, and\n",
    "2. the waves propagate at the (known) speeds $\\xi_{(\\pm)}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From these assumptions, and the Rankine-Hugoniot conditions, we have the two equations\n",
    "\n",
    "\\begin{align}\n",
    "  \\xi_{(-)} \\left[ {\\bf q}_m - {\\bf q}_l \\right] & = {\\bf f}_m - {\\bf f}_l, \\\\\n",
    "  \\xi_{(+)} \\left[ {\\bf q}_r - {\\bf q}_m \\right] & = {\\bf f}_r - {\\bf f}_m.\n",
    "\\end{align}\n",
    "\n",
    "These are immediately solved to give\n",
    "\n",
    "\\begin{align}\n",
    "  {\\bf q}_m & = \\frac{\\xi_{(+)} {\\bf q}_r - \\xi_{(-)} {\\bf q}_l - {\\bf f}_r + {\\bf f}_l}{\\xi_{(+)} - \\xi_{(-)}}, \\\\\n",
    "  {\\bf f}_m & = \\frac{\\hat{\\xi}_{(+)} {\\bf f}_l - \\hat{\\xi}_{(-)} {\\bf f}_r + \\hat{\\xi}_{(+)} \\hat{\\xi}_{(-)} \\left( {\\bf q}_r - {\\bf q}_r \\right)}{\\hat{\\xi}_{(+)} - \\hat{\\xi}_{(-)}},\n",
    "\\end{align}\n",
    "\n",
    "where\n",
    "\n",
    "\\begin{equation}\n",
    "  \\hat{\\xi}_{(-)} = \\min(0, \\xi_{(-)}), \\qquad \\hat{\\xi}_{(+)} = \\max(0, \\xi_{(+)}).\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again it should be noted that, in general, ${\\bf f}_m \\ne {\\bf f}({\\bf q}_m)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We still need some way to compute the wave speeds $\\xi_{(\\pm)}$. The simplest method is to make them as large as possible, compatible with stability. This means (via the CFL condition) setting\n",
    "\n",
    "\\begin{equation}\n",
    "  -\\xi_{(-)} = \\xi_{(+)} = \\frac{\\Delta x}{\\Delta t}\n",
    "\\end{equation}\n",
    "\n",
    "which implies that (as the central state is now guaranteed to include the origin, as the waves have different signs)\n",
    "\n",
    "\\begin{equation}\n",
    "  {\\bf f}^* = \\frac{1}{2} \\left( {\\bf f}_l + {\\bf f}_r + \\frac{\\Delta x}{\\Delta t} \\left[ {\\bf q}_l - {\\bf q}_r \\right] \\right).\n",
    "\\end{equation}\n",
    "\n",
    "This is the *Lax-Friedrichs* flux, as [used in HyperPython](https://github.com/ketch/HyperPython). We can also easily see how the *local* Lax-Friedrichs method, [used in lesson 3 of HyperPython](http://nbviewer.ipython.org/github/ketch/HyperPython/blob/master/Lesson_03_High-resolution_methods.ipynb), comes about: simply choose\n",
    "\n",
    "\\begin{equation}\n",
    "  -\\xi_{(-)} = \\xi_{(+)} = \\alpha = \\min \\left( \\left| \\lambda \\left( {\\bf q}_l \\right) \\right|, \\left| \\lambda \\left( {\\bf q}_r \\right) \\right| \\right)\n",
    "\\end{equation}\n",
    "\n",
    "to get\n",
    "\n",
    "\\begin{equation}\n",
    "  {\\bf f}^* = \\frac{1}{2} \\left( {\\bf f}_l + {\\bf f}_r + \\alpha \\left[ {\\bf q}_l - {\\bf q}_r \\right] \\right).\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "HLL type methods are straightforward to use but typically do not capture linear waves, such as the contact wave in the Euler equations, well. Extending the HLL method by including more waves is possible (see the *HLLC* method in Toro's book as an example), but rapidly increases the complexity of the solver."
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}