{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison of full discretizations using approximate Riemann solvers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How do these solvers impact the solution accuracy when used within a finite volume discretization?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide"
    ]
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from exact_solvers import Euler\n",
    "from clawpack import riemann\n",
    "from utils import riemann_tools\n",
    "import matplotlib.pyplot as plt\n",
    "from collections import namedtuple\n",
    "from ipywidgets import interact\n",
    "from ipywidgets import widgets\n",
    "Primitive_State = namedtuple('State', Euler.primitive_variables)\n",
    "gamma = 1.4\n",
    "problem_data = {}\n",
    "problem_data['gamma'] = gamma\n",
    "problem_data['gamma1'] = gamma - 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from clawpack.riemann.euler_with_efix_1D_constants import density, momentum, energy, num_eqn\n",
    "\n",
    "def shocktube(q_l, q_r, N=50, riemann_solver='HLL', solver_type='classic'):\n",
    "\n",
    "    from clawpack import pyclaw\n",
    "    from clawpack import riemann\n",
    "\n",
    "    if riemann_solver == 'Roe':\n",
    "        rs = riemann.euler_1D_py.euler_roe_1D\n",
    "    elif riemann_solver == 'HLL':\n",
    "        rs = riemann.euler_1D_py.euler_hll_1D\n",
    "\n",
    "    if solver_type == 'classic':\n",
    "        solver = pyclaw.ClawSolver1D(rs)        \n",
    "    else:\n",
    "        solver = pyclaw.SharpClawSolver1D(rs)\n",
    "\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,num_eqn)\n",
    "\n",
    "    gamma = 1.4\n",
    "    state.problem_data['gamma']= gamma\n",
    "    state.problem_data['gamma1']= gamma-1.\n",
    "\n",
    "    state.problem_data['efix'] = False\n",
    "\n",
    "    xc = state.grid.p_centers[0]\n",
    "    \n",
    "    velocity = (xc<=0)*q_l[1] + (xc>0)*q_r[1]\n",
    "    pressure = (xc<=0)*q_l[2] + (xc>0)*q_r[2]\n",
    "\n",
    "    state.q[density ,:] = (xc<=0)*q_l[0] + (xc>0)*q_r[0]\n",
    "    state.q[momentum,:] = velocity * state.q[density,:]\n",
    "    state.q[energy  ,:] = pressure/(gamma - 1.) + 0.5 * state.q[density,:] * velocity**2\n",
    "\n",
    "    claw = pyclaw.Controller()\n",
    "    claw.tfinal = 0.5\n",
    "    claw.solution = pyclaw.Solution(state,domain)\n",
    "    claw.solver = solver\n",
    "    claw.num_output_times = 10\n",
    "    claw.keep_copy = True\n",
    "    claw.verbosity=0\n",
    "\n",
    "    return claw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from clawpack.riemann.euler_with_efix_1D_constants import density, momentum, energy, num_eqn\n",
    "\n",
    "def blastwave(q_l, q_r, N=400, riemann_solver='HLL', solver_type='classic'):\n",
    "\n",
    "    from clawpack import pyclaw\n",
    "    from clawpack import riemann\n",
    "\n",
    "    if riemann_solver == 'Roe':\n",
    "        rs = riemann.euler_1D_py.euler_roe_1D\n",
    "    elif riemann_solver == 'HLL':\n",
    "        rs = riemann.euler_1D_py.euler_hll_1D\n",
    "\n",
    "    if solver_type == 'classic':\n",
    "        solver = pyclaw.ClawSolver1D(rs)        \n",
    "    else:\n",
    "        solver = pyclaw.SharpClawSolver1D(rs)\n",
    "        solver.time_integrator = 'SSP33'\n",
    "        solver.cfl_max = 0.65\n",
    "        solver.cfl_desired = 0.6\n",
    "        solver.lim_type = 1\n",
    "        #solver.char_decomp = 2\n",
    "\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(0.0,1.0,N,name='x')\n",
    "    domain = pyclaw.Domain([x])\n",
    "    state = pyclaw.State(domain,num_eqn)\n",
    "\n",
    "    gamma = 1.4\n",
    "    state.problem_data['gamma']= gamma\n",
    "    state.problem_data['gamma1']= gamma-1.\n",
    "\n",
    "    state.problem_data['efix'] = False\n",
    "\n",
    "    xc = state.grid.p_centers[0]\n",
    "\n",
    "    state.q[density ,:] = 1.\n",
    "    state.q[momentum,:] = 0.\n",
    "    state.q[energy  ,:] = ( (xc<0.1)*1.e3 + (0.1<=xc)*(xc<0.9)*1.e-2 + (0.9<=xc)*1.e2 ) / (gamma - 1.)\n",
    "\n",
    "    claw = pyclaw.Controller()\n",
    "    claw.tfinal = 0.038\n",
    "    claw.solution = pyclaw.Solution(state,domain)\n",
    "    claw.solver = solver\n",
    "    claw.num_output_times = 10\n",
    "    claw.keep_copy = True\n",
    "    claw.verbosity=0\n",
    "\n",
    "    return claw"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Second-order modified Lax-Wendroff algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_l = Euler.conservative_to_primitive(1.,0.,1.)\n",
    "q_r = Euler.conservative_to_primitive(1./8,0.,1./10)\n",
    "roe = shocktube(q_l,q_r,N=50,riemann_solver='Roe')\n",
    "roe.run()\n",
    "hll = shocktube(q_l,q_r,N=50,riemann_solver='HLL')\n",
    "hll.run()\n",
    "fine = shocktube(q_l,q_r,N=1000,riemann_solver='Roe')\n",
    "fine.run();\n",
    "xc = roe.solution.state.grid.p_centers[0]\n",
    "xc_fine = fine.solution.state.grid.p_centers[0]\n",
    "\n",
    "def plot_frame(i):\n",
    "    fig, ax = plt.subplots(figsize=(12,6))\n",
    "    ax.set_xlim((-1,1)); ax.set_ylim((0,1.1))\n",
    "    ax.plot(xc_fine,fine.frames[i].q[density,:],'-k',lw=1)\n",
    "    ax.plot(xc,hll.frames[i].q[density,:],'-ob',lw=2)\n",
    "    ax.plot(xc,roe.frames[i].q[density,:],'-or',lw=2)\n",
    "    plt.legend(['Fine','HLL','Roe'],loc='best')\n",
    "    plt.show()\n",
    "    \n",
    "interact(plot_frame, i=widgets.IntSlider(min=0, max=10, description='Frame'));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you might expect, the HLL solver smears the middle wave (contact discontinuity) significantly more than the Roe solver does. Perhaps surprisingly, it captures the shock just as accurately as the Roe solver does."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_l = Euler.conservative_to_primitive(1.,0.,1.)\n",
    "q_r = Euler.conservative_to_primitive(1./8,0.,1./10)\n",
    "roe = blastwave(q_l,q_r,riemann_solver='Roe')\n",
    "roe.run()\n",
    "hll = blastwave(q_l,q_r,riemann_solver='HLL')\n",
    "hll.run()\n",
    "fine = blastwave(q_l,q_r,N=2000,riemann_solver='Roe')\n",
    "fine.run();\n",
    "xc = roe.solution.state.grid.p_centers[0]\n",
    "xc_fine = fine.solution.state.grid.p_centers[0]\n",
    "\n",
    "def plot_frame(i):\n",
    "    fig, ax = plt.subplots(figsize=(12,6))\n",
    "    ax.set_xlim((0.,1)); ax.set_ylim((0,20))\n",
    "    ax.plot(xc_fine,fine.frames[i].q[density,:],'-k',lw=1)\n",
    "    ax.plot(xc,hll.frames[i].q[density,:],'-ob',lw=2)\n",
    "    ax.plot(xc,roe.frames[i].q[density,:],'-or',lw=2)\n",
    "    plt.legend(['Fine','HLL','Roe'],loc='best')\n",
    "    plt.show()\n",
    "    \n",
    "interact(plot_frame, i=widgets.IntSlider(min=0, max=10, description='Frame'));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## High-order WENO + Runge-Kutta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q_l = Euler.conservative_to_primitive(1.,0.,1.)\n",
    "q_r = Euler.conservative_to_primitive(1./8,0.,1./10)\n",
    "roe2 = shocktube(q_l,q_r,N=50,riemann_solver='Roe',solver_type='sharpclaw')\n",
    "roe2.run()\n",
    "hll2 = shocktube(q_l,q_r,N=50,riemann_solver='HLL',solver_type='sharpclaw')\n",
    "hll2.run()\n",
    "fine2 = shocktube(q_l,q_r,N=1000,riemann_solver='Roe',solver_type='sharpclaw')\n",
    "fine2.run();\n",
    "xc = roe2.solution.state.grid.p_centers[0]\n",
    "xc_fine2 = fine2.solution.state.grid.p_centers[0]\n",
    "\n",
    "def plot_frame(i):\n",
    "    fig, ax = plt.subplots(figsize=(12,6))\n",
    "    ax.set_xlim((-1,1)); ax.set_ylim((0,1.1))\n",
    "    ax.plot(xc_fine2,fine2.frames[i].q[density,:],'-k',lw=1)\n",
    "    ax.plot(xc,hll2.frames[i].q[density,:],'-ob',lw=2)\n",
    "    ax.plot(xc,roe2.frames[i].q[density,:],'-or',lw=2)\n",
    "    plt.legend(['Fine','HLL','Roe'],loc='best')\n",
    "    plt.show()\n",
    "    \n",
    "interact(plot_frame, i=widgets.IntSlider(min=0, max=10, description='Frame'));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With higher-order discretizations, the difference in the Riemann solvers is less significant."
   ]
  }
 ],
 "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.13"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}