Lesson_04_Fluid_dynamics.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:70f483d0cdc9b19b25cb4d7712103b68b352d934bdc7ec6ef6d2c6b8e126a633"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from IPython.core.display import HTML\n",
      "css_file = './custom.css'\n",
      "HTML(open(css_file, \"r\").read())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "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=Philosopher:400,700,400italic,700italic' rel='stylesheet' type='text/css'>\n",
        "\n",
        "<style>\n",
        "\n",
        "#notebook_panel { /* main background */\n",
        "    background: #888;\n",
        "    color: #f6f6f6;\n",
        "}\n",
        "\n",
        "div.cell { /* set cell width to about 80 chars */\n",
        "    width: 800px;\n",
        "}\n",
        "\n",
        "div #notebook { /* centre the content */\n",
        "    background: #fff; /* white background for content */\n",
        "    width: 1000px;\n",
        "    margin: auto;\n",
        "    padding-left: 1em;\n",
        "}\n",
        "\n",
        "#notebook li { /* More space between bullet points */\n",
        "margin-top:0.8em;\n",
        "}\n",
        "\n",
        "/* draw border around running cells */\n",
        "div.cell.border-box-sizing.code_cell.running { \n",
        "    border: 3px solid #111;\n",
        "}\n",
        "\n",
        "/* Put a solid color box around each cell and its output, visually linking them together */\n",
        "div.cell.code_cell {\n",
        "    background-color: rgba(171,165,131,1.0); \n",
        "    border-radius: 10px; /* rounded borders */\n",
        "    padding: 1em;\n",
        "    margin-top: 1em;\n",
        "}\n",
        "\n",
        "div.text_cell_render{\n",
        "    font-family: 'Arvo' sans-serif;\n",
        "    line-height: 130%;\n",
        "    font-size: 115%;\n",
        "    width:700px;\n",
        "    margin-left:auto;\n",
        "    margin-right:auto;\n",
        "}\n",
        "\n",
        "\n",
        "/* Formatting for header cells */\n",
        ".text_cell_render h1 {\n",
        "    font-family: 'Philosopher', sans-serif;\n",
        "    font-weight: 400;\n",
        "    font-size: 40pt;\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-family: 'Philosopher', serif;\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-family: 'Philosopher', serif;\n",
        "    margin-top:12px;\n",
        "    margin-bottom: 3px;\n",
        "    font-style: italic;\n",
        "    color: rgb(95,92,72);\n",
        "}\n",
        "\n",
        ".text_cell_render h4 {\n",
        "    font-family: 'Philosopher', serif;\n",
        "}\n",
        "\n",
        ".text_cell_render h5 {\n",
        "    font-family: 'Alegreya Sans', sans-serif;\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-family: 'PT Mono', sans-serif;\n",
        "    font-weight: 300;\n",
        "    font-size: 10pt;\n",
        "    color: grey;\n",
        "    margin-bottom: 1px;\n",
        "    margin-top: 1px;\n",
        "}\n",
        "\n",
        ".CodeMirror{\n",
        "    font-family: \"PT Mono\";\n",
        "    font-size: 100%;\n",
        "}\n",
        "\n",
        "</style>\n"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 1,
       "text": [
        "<IPython.core.display.HTML at 0x1027a4710>"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "heading",
     "level": 6,
     "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": "heading",
     "level": 5,
     "metadata": {},
     "source": [
      "version 0.1 - April 2014"
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Fluid dynamics"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In this lesson we will look at the system of hyperbolic PDEs that governs the motions of fluids in the absence of viscosity.  These consist of conservation laws for **mass, momentum**, and **energy**.  Together, they are referred to as the **compressible Euler equations**, or simply the Euler equations."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Mass conservation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We will use $\\rho(x,t)$ to denote the fluid density and $u(x,t)$ for its velocity.  Then the equation for conservation of mass is just the **continuity equation** we discussed in [Lesson 1](Lesson_01_Advection.ipynb):\n",
      "\n",
      "$$\\rho_t + (\\rho u)_x = 0.$$"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Momentum conservation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The momentum is given by the product of density and velocity, i.e. $\\rho u$.  The momentum flux has two components.  First, the momentum is transported in the same way that the density is; this flux is given by the momentum times the density; i.e. $\\rho u^2$.\n",
      "\n",
      "To understand the second term in the momentum flux, we must realize that a fluid is made up of many tiny molecules.  The density and velocity we are modeling are average values over some small region of space.  The individual molecules in that region are not all moving with exactly velocity $u$; that's just their average.  Each molecule also has some additional random velocity component.  These random velocities are what accounts for the **pressure** of the fluid, which we'll denote by $p$.  These velocity components also lead to a net flux of momentum.  Thus the momentum conservation equation is\n",
      "\n",
      "$$(\\rho u)_t + (\\rho u^2 + p)_x = 0.$$"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Energy conservation"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The energy has two components: internal energy $\\rho e$ and kinetic energy $\\rho u^2/2$:\n",
      "\n",
      "$$E = \\rho e + \\frac{1}{2}\\rho u^2.$$\n",
      "\n",
      "Like the momentum flux, the energy flux involves both bulk transport ($Eu$) and transport due to pressure ($pu$):\n",
      "\n",
      "$$E_t + (u(E+p))_x = 0.$$"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Equation of state"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "You may have noticed that we have 4 unknowns (density, momentum, energy, and pressure) but only 3 conservation laws.  We need one more relation to close the system.  That relation, known as the equation of state, expresses how the pressure is related to the other quantities.  We'll focus on the case of an ideal gas, for which\n",
      "\n",
      "$$p = \\rho e (\\gamma-1).$$\n",
      "\n",
      "Here $\\gamma$ is the ratio of specific heats, which for air is approximately 1.4."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "The Euler equations"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can write the three conservation laws as a single system $q_t + f(q)_x = 0$ by defining\n",
      "\\begin{align}\n",
      "q & = \\begin{pmatrix} \\rho \\\\ \\rho u \\\\ E\\end{pmatrix}, & \n",
      "f(q) & = \\begin{pmatrix} \\rho u \\\\ \\rho u^2 + p \\\\ u(E+p)\\end{pmatrix}.\n",
      "\\end{align}"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In three dimensions, the equations are similar.  We have two additional velocity components $v, w$, and their corresponding fluxes.  Additionally, we have to account for fluxes in the $y$ and $z$ directions.  We can write the full system as\n",
      "\n",
      "$$ q_t + f(q)_x + g(q)_y + h(q)_z = 0$$\n",
      "\n",
      "with\n",
      "\n",
      "\\begin{align}\n",
      "q & = \\begin{pmatrix} \\rho \\\\ \\rho u \\\\ \\rho v \\\\ \\rho w \\\\ E\\end{pmatrix}, &\n",
      "f(q) & = \\begin{pmatrix} \\rho u \\\\ \\rho u^2 + p \\\\ \\rho u v \\\\ \\rho u w \\\\ u(E+p)\\end{pmatrix} &\n",
      "g(q) & = \\begin{pmatrix} \\rho v \\\\ \\rho uv \\\\ \\rho v^2 + p \\\\ \\rho v w \\\\ v(E+p)\\end{pmatrix} &\n",
      "h(q) & = \\begin{pmatrix} \\rho w \\\\ \\rho uw \\\\ \\rho vw \\\\ \\rho w^2 + p \\\\ w(E+p)\\end{pmatrix}.\n",
      "\\end{align}"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Solving the Euler equations"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "These equations can be solved in a manner similar to what we used for advection and traffic flow.  As you might guess, computing the flux gets significantly more complicated since we now have 3 (or 5) equations and more complicated flux expressions."
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "PyClaw"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Implementing a solver for the Euler equations from scratch would be a lot of fun, but to save some time we'll use a  package called [PyClaw](http://clawpack.github.io/doc/pyclaw/), which is part of the [Clawpack](http://clawpack.github.io/) software (Clawpack stands for Conservation LAWs PACKage).  PyClaw allows us to quickly and easily set up and solve problems modeled by hyperbolic PDEs."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now let's get started.  First, import the parts of Clawpack that we'll use:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from clawpack import pyclaw\n",
      "from clawpack import riemann"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Setting up a problem\n",
      "To solve a problem, we'll need to create the following:\n",
      "\n",
      "- A **domain** over which to solve the problem\n",
      "- A **solution**, where we will provide the initial data.  After running, the solution will contain -- you guessed it! -- the solution.\n",
      "- A **solver**, which is responsible for actually evolving the solution in time.  Here we'll need to specify the equations to be solved and the boundary conditions.\n",
      "- A **controller**, which handles the running, output, and can be used for plotting\n",
      "\n",
      "This might sound complicated at first, but stick with me."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's start by creating a controller and specifying the simulation end time:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "claw = pyclaw.Controller()\n",
      "claw.tfinal = 0.1                    # Simulation end time\n",
      "\n",
      "claw.keep_copy = True          # Keep solution data in memory for plotting\n",
      "claw.output_format = None    # Don't write solution data to file\n",
      "claw.num_output_times = 50  # Write 50 output frames"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Riemann solvers\n",
      "\n",
      "The method used to compute the flux between each pair of cells is referred to as a *Riemann solver*.  By specifying a Riemann solver, we will specify the system of PDEs that we want to solve.  So far we have only used very simple approximate Riemann solvers.  Clawpack includes much more sophisticated Riemann solvers for many hyperbolic systems.\n",
      "\n",
      "Place your cursor at the end of the line in the box below and hit 'Tab' (for autocompletion).  You'll see a dropdown list of all the Riemann solvers currently available in PyClaw.  The ones with 'py' at the end of the name are written in pure Python; the others are written in Fortran and wrapped with f2py."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "riemann."
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We'll start with a simple 1D problem, using the Riemann solver `riemann.euler_with_efix_1D`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "riemann_solver = riemann.euler_with_efix_1D\n",
      "claw.solver = pyclaw.ClawSolver1D(riemann_solver)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We also need to specify boundary conditions.  We'll use extrapolation boundary conditions, so that waves simply pass out of the domain:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "claw.solver.all_bcs = pyclaw.BC.extrap"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### The problem domain\n",
      "Next we need to specify the domain and the grid.  We'll solve on the unit line $[0,1]$ using 100 grid cells.  Note that each argument to the Domain constructor must be a tuple:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "domain = pyclaw.Domain( (0.,), (1.,), (100,))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### The initial solution\n",
      "Next we create a solution object that belongs to the controller and extends over the domain we specified:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "claw.solution = pyclaw.Solution(claw.solver.num_eqn,domain)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The initial data is specified in an array named `solution.q`.  The density is contained in `q[0,:]`, the momentum in `q[1,:]`, and the energy in `q[2,:]`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x=domain.grid.x.centers # grid cell centers\n",
      "gam = 1.4 # ratio of specific heats\n",
      "\n",
      "rho_left = 1.0; rho_right = 0.125\n",
      "p_left = 1.0; p_right = 0.1\n",
      "\n",
      "claw.solution.q[0,:] = (x<0.5)*rho_left + (x>=0.5)*rho_right\n",
      "claw.solution.q[1,:] = 0.\n",
      "claw.solution.q[2,:] = ((x<0.5)*p_left + (x>=0.5)*p_right)/(gam-1.0)\n",
      "\n",
      "%matplotlib inline\n",
      "import matplotlib.pyplot as plt\n",
      "plt.plot(x, claw.solution.q[0,:],'-o')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This problem is known as the **Sod shock-tube**.  It amounts to setting up a tube with a thin separator between a high-pressure, high-density region and a low-pressure, low-density region, then suddenly removing the separator."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Next we need to specify the value of $\\gamma$, the ratio of specific heats."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "problem_data = claw.solution.problem_data\n",
      "problem_data['gamma'] = 1.4\n",
      "problem_data['gamma1'] = 0.4"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Finally, let's run the simulation."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "claw.run()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Plotting\n",
      "Now we'll plot the results, which are contained in a list called `claw.frames`.  It's simple to plot a single frame with matplotlib:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from matplotlib import animation\n",
      "import matplotlib.pyplot as plt\n",
      "from clawpack.visclaw.JSAnimation import IPython_display\n",
      "import numpy as np\n",
      "\n",
      "fig = plt.figure()\n",
      "ax = plt.axes(xlim=(0, 1), ylim=(-0.2, 1.2))\n",
      "\n",
      "frame = claw.frames[0]\n",
      "pressure = frame.q[0,:]\n",
      "line, = ax.plot([], [], 'o-', lw=2)\n",
      "\n",
      "def fplot(frame_number):\n",
      "    frame = claw.frames[frame_number]\n",
      "    pressure = frame.q[0,:]\n",
      "    line.set_data(x,pressure)\n",
      "    return line,\n",
      "\n",
      "animation.FuncAnimation(fig, fplot, frames=len(claw.frames), interval=30)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Waves"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In the solution, 3 waves are visible:\n",
      "1. A **shock wave** moving rapidly to the right as the low-density fluid is compressed.\n",
      "2. A **rarefaction** wave moving to the left as the high-density fluid expands.\n",
      "3. A **contact discontinuity** moving more slowly to the right.  This discontinuity in the density separates the region containing fluid that started in the high-pressure region and fluid that started in the low-pressure region."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In fact, the solution of any Riemann problem consists of some combination of these three types of waves.  In the Euler equations, one of the waves is always a contact discontinuity, but each of the other two waves may be a shock or a rarefaction, depending on the left and right states."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Putting it all together"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For convenience, all of the code from the cells above to set up and run the shocktube problem is pasted together below.  Play around with the code.  You might:\n",
      "- Increase the number of grid points to see what the solution converges to.  Notice that the code still runs pretty fast even for larger grids.  This is because the bottom layer of code in PyClaw is compiled Fortran, not Python.\n",
      "- Change the initial left and right states, or set up a completely different initial condition.  See if you can generate a solution with two shock waves, or two rarefaction waves (some physical intuition is helpful here).\n",
      "- Change the ratio of specific heats\n",
      "- Make the boundaries periodic, so that there is a second shock wave moving left from $x=1$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import matplotlib.pyplot as plt\n",
      "from clawpack import pyclaw\n",
      "from clawpack import riemann\n",
      "\n",
      "claw = pyclaw.Controller()\n",
      "claw.tfinal = 0.1\n",
      "\n",
      "claw.keep_copy = True       # Keep solution data in memory for plotting\n",
      "claw.output_format = None   # Don't write solution data to file\n",
      "claw.num_output_times = 50  # Write 50 output frames\n",
      "\n",
      "riemann_solver = riemann.euler_with_efix_1D\n",
      "claw.solver = pyclaw.ClawSolver1D(riemann_solver)\n",
      "\n",
      "claw.solver.all_bcs = pyclaw.BC.extrap\n",
      "\n",
      "domain = pyclaw.Domain( (0.,), (1.,), (100,))\n",
      "x=domain.grid.x.centers # grid cell centers\n",
      "\n",
      "claw.solution = pyclaw.Solution(claw.solver.num_eqn,domain)\n",
      "\n",
      "gam = 1.4 # ratio of specific heats\n",
      "claw.solution.problem_data['gamma'] = gam\n",
      "claw.solution.problem_data['gamma1'] = gam-1.0\n",
      "\n",
      "rho_left = 1.0; rho_right = 0.125\n",
      "p_left = 1.0; p_right = 0.1\n",
      "\n",
      "claw.solution.q[0,:] = (x<0.5)*rho_left + (x>=0.5)*rho_right\n",
      "claw.solution.q[1,:] = 0.\n",
      "claw.solution.q[2,:] = ((x<0.5)*p_left + (x>=0.5)*p_right)/(gam-1.0)\n",
      "\n",
      "status = claw.run()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}