{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Finite volume discretizations with approximate Riemann solvers" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the previous few chapters we introduced several approximate Riemann solvers for different nonlinear conservation laws and showed some comparisons of the approximate solution to a single Riemann problem with the true solution, where there often appear to be significant differences. However, their intended use is as a building block for finite volume methods such as Godunov's method or high-resolution extensions. \n", "\n", "How do these solvers impact the solution accuracy when used within a finite volume discretization? To investigate, we will use them within [PyClaw](http://www.clawpack.org/pyclaw/) to solve several test problems.\n", "\n", "In particular, for both the shallow water equations and the Euler equations, we will apply the numerical method to solve dam break and shock tube problems, which are themselves Riemann problems and so we can compute the exact solution, but once discretized require the solution of a different set of Riemann problems at each cell interface every time step. \n", "\n", "For the Euler equations we will also consider the Woodward-Colella blast wave problem. The initial data consists of two Riemann problems, with resulting shock waves of differing strengths. These shock waves later interact with each other.\n", " \n", "In this chapter we include extensive sections of code in the notebook. This is meant to more easily allow the reader to use these as templates for setting up other problems in PyClaw. For the code in this chapter we use approximate Riemann solvers from PyClaw that can be found in these files:\n", "\n", "- [euler_1D_py.py,](https://github.com/clawpack/riemann/blob/FA16/riemann/euler_1D_py.py)\n", "- [shallow_roe_with_efix_1D.](https://github.com/clawpack/riemann/blob/FA16/riemann/shallow_1D_py.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 numpy as np\n", "from exact_solvers import euler\n", "from clawpack import riemann\n", "import matplotlib.pyplot as plt\n", "from ipywidgets import interact\n", "from ipywidgets import widgets" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Shallow water equations\n", "\n", "We compare results obtained with using the high-resolution wave propagation method implemented in PyClaw (with limiters to avoid nonphysical oscillations). We compare the results obtained when combined with two different approximate Riemann solvers: the Roe solver and HLLE." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from exact_solvers import shallow_water as sw\n", "from clawpack import riemann\n", "from clawpack.riemann.shallow_roe_with_efix_1D_constants import depth, momentum, num_eqn\n", "\n", "def setup(riemann_solver='roe',N=20,IC='dam-break'):\n", "\n", " from clawpack import pyclaw\n", "\n", " if riemann_solver.lower() == 'roe':\n", " rs = riemann.shallow_roe_with_efix_1D\n", " elif riemann_solver.lower() == 'hlle':\n", " rs = riemann.shallow_hlle_1D\n", " \n", " solver = pyclaw.ClawSolver1D(rs)\n", " \n", " solver.bc_lower[0] = pyclaw.BC.extrap\n", " solver.bc_upper[0] = pyclaw.BC.extrap\n", "\n", " xlower = -5.0\n", " xupper = 5.0\n", " x = pyclaw.Dimension(xlower,xupper,N,name='x')\n", " domain = pyclaw.Domain(x)\n", " state = pyclaw.State(domain,num_eqn)\n", "\n", " # Gravitational constant\n", " state.problem_data['grav'] = 1.0\n", " state.problem_data['dry_tolerance'] = 1e-3\n", " state.problem_data['sea_level'] = 0.0\n", " \n", " xc = state.grid.x.centers\n", "\n", " x0=0.\n", "\n", " hl = 10.\n", " ul = 0.\n", " hr = 0.5\n", " ur = 0.\n", " state.q[depth,:] = hl * (xc <= x0) + hr * (xc > x0)\n", " state.q[momentum,:] = hl*ul * (xc <= x0) + hr*ur * (xc > x0)\n", "\n", " claw = pyclaw.Controller()\n", " claw.keep_copy = True\n", " claw.output_format = None\n", " claw.tfinal = 1.0\n", " claw.solution = pyclaw.Solution(state,domain)\n", " claw.solver = solver\n", "\n", " return claw" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "N = 40 # number of grid cells to use\n", "q_l = [10,0]\n", "q_r = [0.5,0]\n", "\n", "# Roe solution\n", "roe_sw = setup(riemann_solver='roe',N=N)\n", "roe_sw.verbosity = 0\n", "status = roe_sw.run()\n", "xc_sw = roe_sw.grid.x.centers\n", "\n", "# HLLE solution\n", "hlle_sw = setup(riemann_solver='hlle',N=N)\n", "hlle_sw.verbosity = 0\n", "status = hlle_sw.run()\n", "\n", "# Exact solution\n", "xc_exact_sw = np.linspace(-5,5,2000)\n", "states_sw, speeds_sw, reval_sw, wave_types_sw = \\\n", " sw.exact_riemann_solution(q_l, q_r)\n", "def plot_frame(i):\n", " t = roe_sw.frames[i].t+1.e-13\n", " fig, ax = plt.subplots(2,1, sharex=True, figsize=(8,5))\n", " variablenames = [\"Depth\", \"Momentum\"]\n", " variables = [depth, momentum]\n", " ylims = [[-1,12], [-5,12]]\n", " plt.subplots_adjust(hspace=0)\n", " ax[0].title.set_text('Solutions at t={:.2f}'.format(t))\n", " ax[0].set_xlim((-5,5))\n", " ax[1].set(xlabel = 'x')\n", " for j, variable in enumerate(variables):\n", " ax[j].set_ylim(ylims[j])\n", " ax[j].plot(xc_exact_sw,reval_sw(xc_exact_sw/(t+1.e-16))[j],'-k',lw=1)\n", " ax[j].plot(xc_sw, hlle_sw.frames[i].q[variable,:],'-ob',lw=2,markersize=4)\n", " ax[j].plot(xc_sw, roe_sw.frames[i].q[variable,:],'-or',lw=0.5,markersize=3)\n", " ax[j].legend(['Exact','HLLE','Roe'],loc='best')\n", " ax[j].set(ylabel=variablenames[j])\n", " plt.show()\n", " \n", "interact(plot_frame, i=widgets.IntSlider(min=0,max=10,value=5));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Clearly, both solvers give very similar results for this problem. They converge to the same solution, as you can verify by increasing the value of $N$ above in the live notebook. You can also try varying the left and right states." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Euler equations\n", "\n", "We perform a similar test for the Euler equations with shock tube initial data." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Sod shock tube problem" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We first consider the classic shocktube problem proposed by Sod and already discussed in [Euler](Euler.ipynb). This is a particular Riemann problem in which the initial velocity is zero on both sides of a discontinuity in pressure and/or density of a gas, and so the exact Riemann solver for the Euler equations would provide the exact solution for all time, consisting of a right-going shock, a left-going rarefaction, and an intermediate contact discontinuity.\n", "\n", "In the numerical experiments done in this notebook, we use this initial data for a more general finite volume method that could be used to approximate the solution for any initial data. In the first time step there is a single cell interface with nontrivial Riemann data, but as the solution evolves on the grid the Riemann problems that arise in subsequent time steps are very different from the single problem we started with. Depending on the accuracy of the numerical method, the resolution of the grid, and the choice of approximate Riemann solver to use at each grid cell every time step, the numerical solution may deviate significantly from the exact solution to the original shocktube problem. This makes a good initial test problem for numerical methods because the exact solution can be computed for comparison purposes, and because it clearly shows whether the method introduces oscillations around discontinuities and/or smears them out. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from clawpack.riemann.euler_with_efix_1D_constants \\\n", " import density, momentum, energy, num_eqn\n", "\n", "def shocktube(q_l, q_r, N=50, riemann_solver='HLL', \n", " 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", " solver.limiters = pyclaw.limiters.tvd.MC\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.) + \\\n", " 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": [ "N = 80 # number of grid cells to use\n", "\n", "prim_l = [1.,0.,1.]\n", "prim_r = [1./8,0.,1./10]\n", "q_l = euler.conservative_to_primitive(*prim_l)\n", "q_r = euler.conservative_to_primitive(*prim_r)\n", "\n", "# Roe-based solution\n", "roe_st = shocktube(q_l,q_r,N=N,riemann_solver='Roe')\n", "roe_st.run()\n", "xc_st = roe_st.solution.state.grid.p_centers[0]\n", "\n", "# HLL-based solution\n", "hll_st = shocktube(q_l,q_r,N=N,riemann_solver='HLL')\n", "hll_st.run()\n", "\n", "# Exact solution\n", "xc_exact_st = np.linspace(-1,1,2000)\n", "states_st, speeds_st, reval_st, wave_types_st = euler.exact_riemann_solution(prim_l, prim_r)\n", "\n", "def plot_frame(i):\n", " t = roe_st.frames[i].t\n", " fig, ax = plt.subplots(3,1, sharex=True, figsize=(8,6))\n", " variablenames = [\"Density\", \"Momentum\", \"Energy\"]\n", " variables = [density, momentum, energy]\n", " ylims = [[0,1.1], [-0.05,0.35], [0,1.1]]\n", " plt.subplots_adjust(hspace=0)\n", " ax[0].title.set_text('Solutions at t={:.2f}'.format(t))\n", " ax[0].set_xlim((-1,1))\n", " ax[2].set(xlabel = 'x')\n", " for j, variable in enumerate(variables):\n", " ax[j].set_ylim(ylims[j])\n", " ax[j].plot(xc_exact_st,reval_st(xc_exact_st/(t+1.e-16))[j],'-k',lw=1)\n", " ax[j].plot(xc_st,hll_st.frames[i].q[variable,:],'-ob',lw=2,markersize=4)\n", " ax[j].plot(xc_st,roe_st.frames[i].q[variable,:],'-or',lw=0.5,markersize=3)\n", " ax[j].legend(['Exact','HLL','Roe'],loc='best')\n", " ax[j].set(ylabel=variablenames[j])\n", " plt.show()\n", " \n", "interact(plot_frame, i=widgets.IntSlider(min=0, max=10, value=5, 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. In the live notebook you can refine the grid and observe the convergence." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### High-order WENO + Runge-Kutta\n", "\n", "Next we look at the difference between the HLL and Roe solution when these solvers are employed within a higher-order method of lines discretization using fifth-order WENO and a 4th-order Runge-Kutta scheme." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "N = 80 # number of grid cells to use\n", "\n", "prim_l = [1.,0.,1.]\n", "prim_r = [1./8,0.,1./10]\n", "q_l = euler.conservative_to_primitive(*prim_l)\n", "q_r = euler.conservative_to_primitive(*prim_r)\n", "\n", "roe_weno = shocktube(q_l,q_r,N=N,riemann_solver='Roe',solver_type='sharpclaw')\n", "roe_weno.run()\n", "hll_weno = shocktube(q_l,q_r,N=N,riemann_solver='HLL',solver_type='sharpclaw')\n", "hll_weno.run()\n", "\n", "xc_weno = roe_weno.solution.state.grid.p_centers[0]\n", "\n", "# Exact solution\n", "xc_exact_weno = np.linspace(-1,1,2000)\n", "states_weno, speeds_weno, reval_weno, wave_types_weno = euler.exact_riemann_solution(prim_l, prim_r)\n", "\n", "def plot_frame(i):\n", " t = roe_weno.frames[i].t\n", " fig, ax = plt.subplots(3,1, sharex=True, figsize=(8,6))\n", " variablenames = [\"Density\", \"Momentum\", \"Energy\"]\n", " variables = [density, momentum, energy]\n", " ylims = [[0,1.1], [-0.05,0.35], [0,1.1]]\n", " plt.subplots_adjust(hspace=0)\n", " ax[0].title.set_text('Solutions at t={:.2f}'.format(t))\n", " ax[0].set_xlim((-1,1))\n", " ax[2].set(xlabel = 'x')\n", " for j, variable in enumerate(variables):\n", " ax[j].set_ylim(ylims[j])\n", " ax[j].plot(xc_exact_weno,reval_weno(xc_exact_weno/(t+1.e-16))[j],'-k',lw=1)\n", " ax[j].plot(xc_weno,hll_weno.frames[i].q[variable,:],'-ob',lw=2,markersize=4)\n", " ax[j].plot(xc_weno,roe_weno.frames[i].q[variable,:],'-or',lw=0.5,markersize=3)\n", " ax[j].legend(['Exact','HLL','Roe'],loc='best')\n", " ax[j].set(ylabel=variablenames[j])\n", " plt.show()\n", " \n", "interact(plot_frame, i=widgets.IntSlider(min=0, max=10, value=5, description='Frame'));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With higher-order discretizations, the difference in solutions due to using different Riemann solvers is less significant. This is partly because these high-order schemes use more accurate values as inputs to the Riemann problem, so that in smooth regions the jump between most cells is very small." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Woodward-Colella blast wave" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we consider the Woodward-Colella blast wave problem, which is discussed for example in (LeVeque 2002). Here the initial velocity is zero and the density is one everywhere. The pressure is\n", "\\begin{align}\n", " p_0(x) = \\begin{cases} 1000 & 0 \\le x \\le 0.1 \\\\\n", " 0.01 & 0.1 \\le x \\le 0.9 \\\\\n", " 100 & 0.9 \\le x \\le 1\n", " \\end{cases}\n", "\\end{align}\n", "The boundaries at $x=0$ and $x=1$ are solid walls. The solution involves a Riemann problem at $x=0.1$ and another at $x=0.9$. Later, the waves resulting from these Riemann problems interact with each other." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from clawpack.riemann.euler_with_efix_1D_constants \\\n", " import density, momentum, energy, num_eqn\n", "\n", "def blastwave(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", " kernel_language = 'Fortran'\n", " rs = riemann.euler_with_efix_1D\n", " elif riemann_solver == 'HLL':\n", " kernel_language = 'Python'\n", " rs = riemann.euler_1D_py.euler_hll_1D\n", "\n", " if solver_type == 'classic':\n", " solver = pyclaw.ClawSolver1D(rs)\n", " solver.limiters = pyclaw.limiters.tvd.MC\n", " else:\n", " solver = pyclaw.SharpClawSolver1D(rs)\n", "\n", " solver.kernel_language = kernel_language\n", " \n", " solver.bc_lower[0]=pyclaw.BC.wall\n", " solver.bc_upper[0]=pyclaw.BC.wall\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", " pressure = (xc<0.1)*1.e3 + (0.1<=xc)*(xc<0.9)*1.e-2 + (0.9<=xc)*1.e2\n", " \n", " state.q[density ,:] = 1.\n", " state.q[momentum,:] = 0.\n", " state.q[energy ,:] = pressure / (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 = 30\n", " claw.keep_copy = True\n", " claw.verbosity=0\n", "\n", " return claw" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "N = 400 # number of grid cells to use\n", "\n", "roe_bw = blastwave(N=N, riemann_solver='Roe')\n", "roe_bw.run()\n", "hll_bw = blastwave(N=N, riemann_solver='HLL')\n", "hll_bw.run()\n", "fine_bw = blastwave(N=4000,riemann_solver='Roe')\n", "fine_bw.run();\n", "xc_bw = roe_bw.solution.state.grid.p_centers[0]\n", "xc_fine_bw = fine_bw.solution.state.grid.p_centers[0]\n", "\n", "def plot_frame(i):\n", " t = roe_bw.frames[i].t\n", " fig, ax = plt.subplots(3,1, sharex=True, figsize=(8,6))\n", " variablenames = [\"Density\", \"Momentum\", \"Energy\"]\n", " variables = [density, momentum, energy]\n", " ylims = [[-1,10], [-50,120], [-300,3000]]\n", " plt.subplots_adjust(hspace=0)\n", " ax[0].title.set_text('Solutions at t={:.3f}'.format(t))\n", " ax[0].set_xlim((0,1))\n", " ax[2].set(xlabel = 'x')\n", " for j, variable in enumerate(variables):\n", " ax[j].set_ylim(ylims[j])\n", " ax[j].plot(xc_fine_bw,fine_bw.frames[i].q[variable,:],'-k',lw=1)\n", " ax[j].plot(xc_bw,hll_bw.frames[i].q[variable,:],'-b',lw=2)\n", " ax[j].plot(xc_bw,roe_bw.frames[i].q[variable,:],'--r',lw=2)\n", " ax[j].legend(['Fine','HLL','Roe'],loc='best')\n", " ax[j].set(ylabel=variablenames[j])\n", " plt.show()\n", " \n", "interact(plot_frame, i=widgets.IntSlider(min=0, max=30, value=15, description='Frame'));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here no exact solution is available, so we compare with a solution computed on a finer grid.\n", "Again the solutions are fairly similar, though the HLL solution is a bit more smeared.\n", "\n", "One should not conclude from these tests that, for instance, the Roe solver is always *better* than the HLL solver. Many factors besides accuracy should be considered, including cost and robustness. As we have seen the HLL solver is more robust in the presence of near-vacuum states." ] } ], "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 }