Created
March 29, 2017 17:27
-
-
Save ketch/cc7b0d015c81cdcab660b9e4bf96a8e3 to your computer and use it in GitHub Desktop.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "# The shallow water equations" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "In this chapter we study a model for shallow water waves in one dimension:\n", | |
| "\n", | |
| "\\begin{align}\n", | |
| " h_t + (hu)_x & = 0 \\label{SW_mass} \\\\\n", | |
| " (hu)_t + \\left(hu^2 + \\frac{1}{2}gh^2\\right)_x & = 0. \\label{SW_mom}\n", | |
| "\\end{align}\n", | |
| "\n", | |
| "Here $h$ is the depth, $u$ is the velocity, and $g$ is a constant representing the force of gravity. These equations are a relatively simple but surprisingly effective model for water waves; they are bassed on several important physical assumptions; for instance, it ignores viscosity and compressibility, and assumes the pressure is hydrostatic. Nevertheless, it is a surprisingly effective model for many applications.\n", | |
| "\n", | |
| "Previously we have looked at a *linear hyperbolic system* (acoustics) and *scalar hyperbolic equations* (Burgers', traffic flow). The shallow water system is our first example of a *nonlinear hyperbolic system*; solutions of the Riemann problem will consist of multiple waves, each of which may be a shock or rarefaction.\n", | |
| "\n", | |
| "The discussion here largely follows that of (Leveque, FVM, Chapter 13)." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "## Hyperbolic structure\n", | |
| "We can write (\\ref{SW_mass})-(\\ref{SW_mom}) in the canonical form $q_t + f(q)_x = 0$ if we define\n", | |
| "\n", | |
| "\\begin{align}\n", | |
| "q & = \\begin{pmatrix} h \\\\ hu \\end{pmatrix} & f & = \\begin{pmatrix} hu \\\\ hu^2 + \\frac{1}{2}gh^2 \\end{pmatrix}.\n", | |
| "\\end{align}\n", | |
| "\n", | |
| "In terms of the conserved quantities, the flux is\n", | |
| "\n", | |
| "\\begin{align}\n", | |
| "f(q) & = \\begin{pmatrix} q_2 \\\\ q_2^2/q_1 + \\frac{1}{2}g q_1^2 \\end{pmatrix}\n", | |
| "\\end{align}\n", | |
| "\n", | |
| "Thus the flux jacobian is\n", | |
| "\\begin{align}\n", | |
| "f'(q) & = \\begin{pmatrix} 0 & 1 \\\\ -(q_2/q_1)^2 + g q_1 & 2 q_2/q_1 \\end{pmatrix} \n", | |
| " = \\begin{pmatrix} 0 & 1 \\\\ -u^2 + g h & 2 u \\end{pmatrix}.\n", | |
| "\\end{align}\n", | |
| "\n", | |
| "Its eigenvalues are\n", | |
| "\\begin{align}\n", | |
| " \\lambda_1 & = u - \\sqrt{gh} & \\lambda_2 & = u + \\sqrt{gh},\n", | |
| "\\end{align}\n", | |
| "with corresponding eigenvectors\n", | |
| "\\begin{align} \\label{SW:fjac-evecs}\n", | |
| " r_1 & = \\begin{bmatrix} 1 \\\\ u-\\sqrt{gh} \\end{bmatrix} &\n", | |
| " r_2 & = \\begin{bmatrix} 1 \\\\ u+\\sqrt{gh} \\end{bmatrix}\n", | |
| "\\end{align}\n", | |
| "\n", | |
| "Notice that -- unlike for acoustics, but similar to the LWR traffic model -- the eigenvalues and eigenvectors depend on $q$. Because of this, the waves appearing in the Riemann problem move at different speeds and may be either jump discontinuities (shocks) or smoothly-varying regions (rarefactions)." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "## The Riemann problem\n", | |
| "Consider the Riemann problem with left and right states\n", | |
| "\\begin{align*}\n", | |
| "q_l & = \\begin{bmatrix} h_l \\\\ h_r u_l \\end{bmatrix} &\n", | |
| "q_r & = \\begin{bmatrix} h_r \\\\ h_r u_r \\end{bmatrix}. \n", | |
| "\\end{align*}\n", | |
| "Typically the Riemann solution will consist of two waves, one related to each of the eigenvectors in \\eqref{SW:fjac-evecs}. Each wave may be a shock or rarefaction. There will be an intermediate state $q_m$ between them. Here we illustrate a typical situation, in which the 1-wave happens to be a rarefaction and the 2-wave is a shock:\n", | |
| "\n", | |
| "\n", | |
| "\n", | |
| "In the sketch above we have one wave going in each direction, but since the wave speeds depend on $q$ and can have either sign, it is possible to have both waves going left, or both going right.\n", | |
| "\n", | |
| "To solve the Riemann problem, we must find $q_m$. To do so we must find a state that can be connected to $q_l$ by a 1-shock or 1-rarefaction and to $q_r$ by a 2-shock or 2-rarefaction. We must also ensure that the resulting waves satisfy the entropy condition." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "# Shock waves\n", | |
| "Let us now consider a shock wave connecting a known state $q_*=(h_*, h_* u_*)$ to an unknown state $q=(h,hu)$. In the context of the Riemann problem $q_*$ will be the left or right state and $q$ will be the middle state.\n", | |
| "The Rankine-Hugoniot jump conditions $s(q_r - q_l) = f(q_r) - f(q_l)$ for a shock wave moving\n", | |
| "at speed $s$ read\n", | |
| "\n", | |
| "\\begin{align}\n", | |
| "s (h_* - h) & = h_* u_* - h u \\\\\n", | |
| "s (h_* u_* - h u) & = h_* u_*^2 - h u^2 + \\frac{g}{2}(h_*^2 - h^2).\n", | |
| "\\end{align}\n", | |
| "\n", | |
| "It is convenient to define $\\alpha = h - h_*$. Then we can combine the two equations above to show that the momenta are related by\n", | |
| "\n", | |
| "\\begin{align} \\label{SW:hugoniot-locus}\n", | |
| " h u & = h_* u_* + \\alpha \\left[u_* \\pm \\sqrt{gh_* \\left(1+\\frac{\\alpha}{h_*}\\right)\\left(1+\\frac{\\alpha}{2h_*}\\right)}\\right]\n", | |
| "\\end{align}\n", | |
| "\n", | |
| "for a shock moving with speed\n", | |
| "\n", | |
| "\\begin{align} \\label{SW:shock-speed}\n", | |
| " s & = \\frac{h_* u_* - h u}{h_* - h} = u_* \\mp \\sqrt{gh_* \\left(1+\\frac{\\alpha}{h_*}\\right)\\left(1+\\frac{\\alpha}{2h_*}\\right)} = u_* \\mp \\sqrt{gh \\frac{h_* + h}{2h_*} }\n", | |
| "\\end{align}\n", | |
| "\n", | |
| "Choosing the plus sign in \\eqref{SW:hugoniot-locus} (which yields the minus sign in \\eqref{SW:shock-speed}) gives a 1-shock. Choosing the minus sign in \\eqref{SW:hugoniot-locus} (which yields the plus sign in \\eqref{SW:shock-speed}) gives a 2-shock. Notice from the last expression in \\eqref{SW:shock-speed} that as $h \\to h^*$ the shock speed approaches the corresponding characteristic speed, while the ratio of the jumps in momentum and depth approaches that suggested by the eigenvectors \\eqref{SW:fjac-evecs}.\n", | |
| "\n", | |
| "We can now consider the set of all possible states $(h, h u)$ given by \\eqref{SW:hugoniot-locus}. This curve in the $h - hu$ plane is known as the hugoniot locus; there is a family of hugoniot loci for 1-shocks and another for 2-shocks. Below we plot some members of these two families of curves." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 1, | |
| "metadata": { | |
| "collapsed": false, | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "%matplotlib inline\n", | |
| "import numpy as np\n", | |
| "import matplotlib.pyplot as plt\n", | |
| "from exact_solvers import shallow_water\n", | |
| "from collections import namedtuple\n", | |
| "from utils import riemann_tools\n", | |
| "from ipywidgets import interact, widgets\n", | |
| "Primitive_State = namedtuple('State', shallow_water.primitive_variables)\n", | |
| "plt.style.use('seaborn-talk')\n", | |
| "g = 1." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 2, | |
| "metadata": { | |
| "collapsed": false, | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "application/vnd.jupyter.widget-view+json": { | |
| "model_id": "96855aa295344b178ce009a3df64991a" | |
| } | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "interact(shallow_water.plot_hugoniot_loci,y_axis=widgets.fixed('hu'));" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "The plot above shows the Hugoniot loci in the $h$-$hu$ plane, the natural phase plane in terms of the two conserved quantities. Note that they all approach the origin as $h \\rightarrow 0$. Alternatively, we can plot these same curves in the $h$-$u$ plane:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 3, | |
| "metadata": { | |
| "collapsed": false, | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "application/vnd.jupyter.widget-view+json": { | |
| "model_id": "2d06aeffe658459495b81531d855fdc6" | |
| } | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "interact(shallow_water.plot_hugoniot_loci,y_axis=widgets.fixed('u'));" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "Note that in this plane, the curves in each family are simply vertical translates of one another, and all curves asymptote to $\\pm \\infty$ as $h \\rightarrow 0$." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "## All-shock Riemann solutions\n", | |
| "If we knew *a priori* that both waves in the Riemann solution were shock waves, we could find $q_m$ by finding the intersection of the 1-Hugoniot locus passing through $q_l$ with the 2-Hugoniot locus passing through $q_r$. The Riemann solution would then consist simply of two discontinuities propagating at the speeds given by \\eqref{SW:shock-speed}.\n", | |
| "\n", | |
| "Here is an example for which this is the correct solution." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 11, | |
| "metadata": { | |
| "collapsed": false, | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "def c1(q, x, g=1.):\n", | |
| " \"Characteristic speed for shallow water 1-waves.\"\n", | |
| " h = q[0]\n", | |
| " u = q[1]/q[0]\n", | |
| " return u - np.sqrt(g*h)\n", | |
| "\n", | |
| "def c2(q, x, g=1.):\n", | |
| " \"Characteristic speed for shallow water 2-waves.\"\n", | |
| " h = q[0]\n", | |
| " u = q[1]/q[0]\n", | |
| " return u + np.sqrt(g*h)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 12, | |
| "metadata": { | |
| "collapsed": false, | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "def make_plot_function(q_l,q_r,g=1.,force_waves=None):\n", | |
| " states, speeds, reval, wave_types = shallow_water.exact_riemann_solution(q_l,q_r,g,force_waves=force_waves)\n", | |
| " def plot_function(t,plot_1_chars=False,plot_2_chars=False):\n", | |
| " ax = riemann_tools.plot_riemann(states,speeds,reval,wave_types,t=t,t_pointer=0,\n", | |
| " extra_axes=True,conserved_variables=shallow_water.conserved_variables);\n", | |
| " if plot_1_chars:\n", | |
| " riemann_tools.plot_characteristics(reval,c1,ax[0])\n", | |
| " if plot_2_chars:\n", | |
| " riemann_tools.plot_characteristics(reval,c2,ax[0])\n", | |
| " shallow_water.phase_plane_plot(q_l,q_r,g,ax=ax[3],force_waves=force_waves)\n", | |
| " plt.show()\n", | |
| " return plot_function\n", | |
| "\n", | |
| "def plot_riemann_SW(q_l,q_r,g=1.,force_waves=None):\n", | |
| " plot_function = make_plot_function(q_l,q_r,g,force_waves,)\n", | |
| " interact(plot_function, t=widgets.FloatSlider(value=0.,min=0,max=.9))" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 13, | |
| "metadata": { | |
| "collapsed": false, | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "application/vnd.jupyter.widget-view+json": { | |
| "model_id": "c9fb77acfb36406590485990d60f58e7" | |
| } | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "q_l = Primitive_State(Depth = 1.,\n", | |
| " Velocity = 1.)\n", | |
| "q_r = Primitive_State(Depth = 1.,\n", | |
| " Velocity = -1.)\n", | |
| "\n", | |
| "plot_riemann_SW(q_l,q_r)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "Check the boxes to plot the characteristic families, and notice that the 1-characteristics impinge on the 1-shock; the 2-characteristics impinge on the 2-shock. Thus these shocks satisfy the entropy condition." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "### Dam-break problem: all-shock solution\n", | |
| "Here's another example, which you might think of as modeling a dam that suddenly breaks. The water is initially at rest on both sides, but much deeper on the left." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 14, | |
| "metadata": { | |
| "collapsed": false, | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "application/vnd.jupyter.widget-view+json": { | |
| "model_id": "3449326df5094fe28d63e9244a2d7d0a" | |
| } | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "q_l = Primitive_State(Depth = 4.,\n", | |
| " Velocity = 0.)\n", | |
| "q_r = Primitive_State(Depth = 1.,\n", | |
| " Velocity = 0.)\n", | |
| "\n", | |
| "plot_riemann_SW(q_l,q_r,force_waves='shock')" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "Check the box to plot the 1-characteristics. Notice that they don't impinge on the 1-shock; instead, the characteristics to either side of the shock are diverging away from it. This shock does not satisfy the entropy condition and should be replaced with a rarefaction. The corresponding Hugoniot locus is plotted with a dashed line to remind us that it is unphysical.\n", | |
| "\n", | |
| "Now check the box to plot the 2-characteristics. Notice that these *do* impinge on the 2-shock; this is an entropy-satisfying shock, and the Hugoniot locus is plotted as a solid line." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "## The entropy condition\n", | |
| "To be physically correct, each of these waves must satisfy the Lax entropy condition. Notice that the shock speed \\eqref{SW:shock-speed} can be written as\n", | |
| "\\begin{align}\n", | |
| " s & = u_* \\pm \\sqrt{gh \\frac{h_* + h}{2h_*} }\n", | |
| "\\end{align}\n", | |
| "Thus the 1-shock entropy condition reads\n", | |
| "\\begin{align*}\n", | |
| "u_l - \\sqrt{gh_l} > u_l - \\sqrt{gh_l \\frac{h_l+h_m}{2h_l}} = u_m - \\sqrt{gh_m \\frac{h_m+h_l}{2h_m}} > u_m - \\sqrt{gh_m}.\n", | |
| "\\end{align*}\n", | |
| "The second expression for the shock speed is obtained by noticing\n", | |
| "that the shock speed is invariant under transposition of the left and right states.\n", | |
| "\n", | |
| "The first and second inequalities both simplify to\n", | |
| "\\begin{align} \\label{SW:1-entropy}\n", | |
| " h_m > h_l.\n", | |
| "\\end{align}\n", | |
| "Similarly, it can be shown that the entropy condition for a 2-shock connecting $q_m$ and $q_r$ is satisfied if and only if\n", | |
| "\\begin{align} \\label{SW:2-entropy}\n", | |
| " h_m > h_r.\n", | |
| "\\end{align}\n", | |
| "If \\eqref{SW:1-entropy} or \\eqref{SW:2-entropy} is violated, then the\n", | |
| "corresponding wave must be a rarefaction rather than a shock." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "# Rarefaction waves\n", | |
| "\n", | |
| "## Integral curves\n", | |
| "In the LWR traffic flow model, we saw that a rarefaction wave\n", | |
| "consisted of a smoothly-varying density. For the shallow water\n", | |
| "system, both the depth $h$ and the momentum $hu$ will vary smoothly\n", | |
| "through a rarefaction wave. A given rarefaction wave is associated\n", | |
| "with just one characteristic family, and the variation within the\n", | |
| "wave is at all points proportional to the corresponding eigenvector\n", | |
| "$r_p$. In other words, all values $\\tilde{q}=(h,hu)$ in the rarefaction\n", | |
| "lie on a curve defined by\n", | |
| "\n", | |
| "\\begin{align} \\label{SW:gen_ode}\n", | |
| " \\tilde{q}'(\\xi) & = r_p(\\tilde{q}(\\xi)).\n", | |
| "\\end{align}\n", | |
| "Here $\\tilde{q}(\\xi)$ is a parameterization of the solution.\n", | |
| "For the shallow water system, \\eqref{SW:gen_ode} reads\n", | |
| "\\begin{align}\n", | |
| " h'(\\xi) = q_1'(\\xi) & = 1 \\label{SW:dh} \\\\\n", | |
| " (hu)'(\\xi) = q_2'(\\xi) & = u \\pm \\sqrt{gh} = \\tilde{q}_2/\\tilde{q}_1 - \\sqrt{g\\tilde{q}_1}, \\label{SW:dhu}\n", | |
| "\\end{align}\n", | |
| "where we take the minus sign for 1-waves and the plus sign for 2-waves.\n", | |
| "We can think of \\eqref{SW:dh}-\\eqref{SW:dhu} as an initial value ODE;\n", | |
| "fixing the value of $\\tilde{q}$ at one point in the rarefaction wave\n", | |
| "determines the whole solution of \\eqref{SW:dh}-\\eqref{SW:dhu}. We\n", | |
| "refer to each of these solutions as an *integral curve*." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "For the shallow water system, these equations can be integrated exactly. If we fix one point $(h_*, h_*u_*)$ on the curve, the whole integral curve is defined by\n", | |
| "\n", | |
| "\\begin{align}\n", | |
| " hu & = hu_* \\pm 2h(\\sqrt{gh_*} - \\sqrt{gh}).\n", | |
| "\\end{align}\n", | |
| "\n", | |
| "where now the plus sign is for 1-waves and the minus sign for 2-waves.\n", | |
| "This can be rewritten as\n", | |
| "\n", | |
| "\\begin{align} \\label{SW:riemann-invariant}\n", | |
| "u \\pm 2 \\sqrt{gh} = u_* \\pm 2 \\sqrt{gh_*}.\n", | |
| "\\end{align}\n", | |
| "From \\eqref{SW:riemann-invariant} we see that the value $u + 2 \\sqrt{gh}$\n", | |
| "is constant along any 1-integral curve, while the value $u-2\\sqrt{gh}$ is constant along any 2-integral curve. We refer to these quantities as *Riemann invariants* for the shallow water system:\n", | |
| "\n", | |
| "\\begin{align}\n", | |
| "w_1(q) & = u + 2 \\sqrt{gh} \\\\\n", | |
| "w_2(q) & = u - 2 \\sqrt{gh}.\n", | |
| "\\end{align}\n", | |
| "\n", | |
| "In other words, the trajectories plotted above are just the level curves of $w_1$ and $w_2$. \n", | |
| "This allows us to easily plot the integral curves as a contour plot of $w_1$ or $w_2$:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 10, | |
| "metadata": { | |
| "collapsed": false, | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "application/vnd.jupyter.widget-view+json": { | |
| "model_id": "8d914320b8134f7ab2f0661dba09f38e" | |
| } | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "def plot_int_curves(plot_1=True,plot_2=False):\n", | |
| " N = 400\n", | |
| " h, hu = np.meshgrid(np.linspace(0.01,3,N),np.linspace(-3,3,N))\n", | |
| " g = 1.\n", | |
| " u = hu/h\n", | |
| " w1 = u + 2 * np.sqrt(g*h)\n", | |
| " w2 = u - 2 * np.sqrt(g*h)\n", | |
| " if plot_1:\n", | |
| " clines = np.linspace(-4,6,21)\n", | |
| " plt.contour(h,hu,w1,clines,colors='coral',linestyles='solid')\n", | |
| " if plot_2:\n", | |
| " clines = np.linspace(-6,4,21)\n", | |
| " plt.contour(h,hu,w2,clines,colors='maroon',linestyles='solid')\n", | |
| "\n", | |
| " plt.axis((0,3,-3,3))\n", | |
| " plt.xlabel('depth h')\n", | |
| " plt.ylabel('momentum hu')\n", | |
| " plt.title('Integral curves')\n", | |
| " plt.show()\n", | |
| "\n", | |
| "interact(plot_int_curves,\n", | |
| " plot_1=widgets.Checkbox(description='1-wave curves',value=True),\n", | |
| " plot_2=widgets.Checkbox(description='2-wave curves'));" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 9, | |
| "metadata": { | |
| "collapsed": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "widgets.Checkbox?" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "We can also plot the integral curves in the $h$--$u$ plane:" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 4, | |
| "metadata": { | |
| "collapsed": false, | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "application/vnd.jupyter.widget-view+json": { | |
| "model_id": "fe5aff9549ac41e898563d134ce56203" | |
| } | |
| }, | |
| "metadata": {}, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "def plot_int_curves_u(plot_1=True,plot_2=False):\n", | |
| " N = 100\n", | |
| " h, u = np.meshgrid(np.linspace(0.01,3,N),np.linspace(-3,3,N))\n", | |
| " g = 1.\n", | |
| " w1 = u + 2 * np.sqrt(g*h)\n", | |
| " w2 = u - 2 * np.sqrt(g*h)\n", | |
| " if plot_1:\n", | |
| " clines = np.linspace(-4,6,21)\n", | |
| " plt.contour(h,u,w1,clines,colors='coral',linestyles='solid')\n", | |
| " if plot_2:\n", | |
| " clines = np.linspace(-6,4,21)\n", | |
| " plt.contour(h,u,w2,clines,colors='maroon',linestyles='solid')\n", | |
| "\n", | |
| " plt.axis((0,3,-3,3))\n", | |
| " plt.xlabel('depth h')\n", | |
| " plt.ylabel('velocity u')\n", | |
| " plt.title('Integral curves')\n", | |
| " plt.show()\n", | |
| "\n", | |
| "interact(plot_int_curves_u);" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "Note that in this plane the integral curves of each family are simply vertical translates of one another due to the form of the functions $w_1$ and $w_2$. Unlike the Hugoniot loci, the integral curves do not asymptote to $\\pm\\infty$ as $h \\rightarrow 0$ and instead approach a finite value. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "## The structure of centered rarefaction waves\n", | |
| "Rarefactions appearing in the Riemann solution are also similarity solutions; that is they are constant along any ray $\\xi=x/t$. For a $p$-rarefaction connecting states $q_\\ell$ and $q_r$, the rarefaction thus has the form\n", | |
| "\n", | |
| "\\begin{align}\n", | |
| " q(x,t) & = \\begin{cases} q_\\ell & \\text{if } x/t \\le \\lambda_p(q_\\ell) \\\\\n", | |
| " \\tilde{q}(x/t) & \\text{if } \\lambda_p(q_\\ell) \\le \\lambda_p(q_r) \\\\\n", | |
| " q_r & \\text{if } x/t \\ge \\lambda_p(q_r). \\end{cases}\n", | |
| "\\end{align}\n", | |
| "\n", | |
| "It remains to be determined how $\\tilde{q}$ varies as a function of $x/t$. This can be determined by noting that each solution value $\\tilde{q}$ propagates at the characteristic speed $\\lambda_p(\\tilde{q})$. Letting $\\xi=x/t$, this gives us the condition\n", | |
| "\n", | |
| "\\begin{align} \\label{SW:char-sim}\n", | |
| " \\lambda_p(\\tilde{q}(\\xi)) = \\xi.\n", | |
| "\\end{align}\n", | |
| "\n", | |
| "By combining \\eqref{SW:char-sim} with the Riemann invariant condition, we can find $\\tilde{h}$ and $\\tilde{u}$ within the rarefaction. For concreteness, let us consider a 1-rarefaction connecting $q_l$ and $q_m$. Then \\eqref{SW:char-sim} reads $\\tilde{u} - \\sqrt{g\\tilde{h}} = \\xi$, which simplifies to\n", | |
| "\n", | |
| "\\begin{align} \\label{SW:char-sim1}\n", | |
| "\\tilde{h} = \\frac{(\\tilde{u}-\\xi)^2}{g}.\n", | |
| "\\end{align}\n", | |
| "\n", | |
| "Meanwhile, the fact that $w_1(\\tilde{q})$ is constant means that\n", | |
| "\n", | |
| "\\begin{align} \\label{SW:w1-const}\n", | |
| " \\tilde{u} + 2 \\sqrt{g\\tilde{h}} = u_\\ell + 2 \\sqrt{gh_\\ell}.\n", | |
| "\\end{align}\n", | |
| "Combining \\eqref{SW:char-sim1} with \\eqref{SW:w1-const}, we find\n", | |
| "\\begin{align}\n", | |
| " \\tilde{h} & = \\frac{(u_\\ell + 2\\sqrt{gh_\\ell} - \\xi)^2}{9g} = \\frac{(w_1(q_l) - \\xi)^2}{9g} \\\\\n", | |
| " \\tilde{u} & = \\frac{1}{3}w_1(q_l) + \\frac{2}{3}\\xi.\n", | |
| "\\end{align}\n", | |
| "We see that $h$ varies quadratically through the rarefaction, while $u$ varies linearly." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "## Connecting states\n", | |
| "Suppose we know that the middle state $q_m$ is connected to the left and right states by rarefactions. Then we can find $q_m$ by finding the intersection of the corresponding integral curves. Here is an example for which this is the physically correct solution:\n", | |
| "\n", | |
| "\\begin{align}\n", | |
| " q_l & = \\begin{bmatrix} 1 \\\\ -1 \\end{bmatrix} &\n", | |
| " q_r & = \\begin{bmatrix} 1 \\\\ 1 \\end{bmatrix}\n", | |
| "\\end{align}" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 5, | |
| "metadata": { | |
| "collapsed": false, | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "outputs": [ | |
| { | |
| "ename": "NameError", | |
| "evalue": "name 'plot_riemann_SW' is not defined", | |
| "output_type": "error", | |
| "traceback": [ | |
| "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", | |
| "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", | |
| "\u001b[0;32m<ipython-input-5-08551a9817ad>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 4\u001b[0m Velocity = 1)\n\u001b[1;32m 5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m \u001b[0mplot_riemann_SW\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mq_l\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mq_r\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", | |
| "\u001b[0;31mNameError\u001b[0m: name 'plot_riemann_SW' is not defined" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "q_l = Primitive_State(Depth = 1.,\n", | |
| " Velocity = -1)\n", | |
| "q_r = Primitive_State(Depth = 1.,\n", | |
| " Velocity = 1)\n", | |
| "\n", | |
| "plot_riemann_SW(q_l,q_r)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "*Note that the velocities can be changed to $2, -2$ and the cell above still works with resulting $h_m \\approx 0$, but it breaks for larger outflow velocities. Say something here about how the $h$--$u$ phase plane plots of integral curves show the dry state solution?*" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "### Dam-break problem: all-rarefaction solution\n", | |
| "Next let us revisit the problem above where we found that the all-shock solution was incorrect. This time we'll try to find an all-rarefaction solution." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "metadata": { | |
| "collapsed": false, | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "q_l = Primitive_State(Depth = 4.,\n", | |
| " Velocity = 0.)\n", | |
| "q_r = Primitive_State(Depth = 1.,\n", | |
| " Velocity = 0.)\n", | |
| "\n", | |
| "plot_riemann_SW(q_l,q_r,force_waves='raref')" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "Plot the 2-characteristics and notice that they cross each other within the 2-rarefaction. This rarefaction is not physical and should be replaced with a shock; the corresponding integral curve is hence shown as a dashed line." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "# The Riemann solution\n", | |
| "To determine $q_m$, we need to know whether each of the two waves is a shock or rarefaction, so that we can use the appopriate Hugoniot locus or integral curve. But to determine the nature of each wave, we need to check \\eqref{SW:1-entropy} or \\eqref{SW:2-entropy}, which requires knowledge of $h_m$. To deal with this is we define a piecewise function $\\phi_l(h)$ that gives:\n", | |
| "\n", | |
| "- for $h>h_l$, the momentum connected to $q_l$ by a 1-shock;\n", | |
| "- for $h<h_l$, the momentum connected to $q_l$ by a 1-integral curve.\n", | |
| " \n", | |
| "Thus $\\phi$ is partly a 1-Hugoniot locus and partly a 1-integral curve; it gives the value of each precisely in the regime where they are physically correct. We can similarly define a function $\\phi_r(h)$ related to the 2-Hugoniot locus and 2-integral curve. The middle state depth $h_m$ is the value for which $\\phi_l(h) = \\phi_r(h)$." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "### Dam-break problem: correct solution\n", | |
| "Below we plot (at last) the correct solution to the dam-break problem, involving a 1-rarefaction and a 2-shock. Plot the characteristics and notice how each family behaves correctly with respect the corresponding wave in the Riemann solution." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "metadata": { | |
| "collapsed": false, | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "q_l = Primitive_State(Depth = 4.,\n", | |
| " Velocity = 0.)\n", | |
| "q_r = Primitive_State(Depth = 1.,\n", | |
| " Velocity = 0.)\n", | |
| "\n", | |
| "plot_riemann_SW(q_l,q_r)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "source": [ | |
| "## Interactive phase plane\n", | |
| "In the phase plane image below, the hugoniot loci are plotted in red and the integral curves in blue. The portion of each curve that satisfies the entropy condition is plotted as a solid line, while the entropy-violating parts are plotted as dashed lines. The curves are plotted for both the left and right states. Thus the function $\\phi_l$ just described consists of the solid curve (blue and red parts) passing through the left state $q_\\ell$; the function $\\phi_r$ is the solid curve (blue and red parts) passing through $q_r$. The intersection of these curves is labeled $q_m$. [Click here to go to the interactive phase plane](phase_plane_shallow_water_small.html). While you're there, be sure to click and drag the dark circles indicating the values of $q_\\ell$ and $q_r$ in the top left plot. The time can also be adjusted by dragging the circle in the top right plot." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "metadata": { | |
| "collapsed": false, | |
| "deletable": true, | |
| "editable": true | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "# Perhaps this is unnecessary\n", | |
| "def plot_exact_riemann_solution(h_l=3.,u_l=0.,h_r=1.,u_r=0.,t=0.2): \n", | |
| " q_l = shallow_water.primitive_to_conservative(h_l,u_l)\n", | |
| " q_r = shallow_water.primitive_to_conservative(h_r,u_r)\n", | |
| " \n", | |
| " x = np.linspace(-1.,1.,1000)\n", | |
| " states, speeds, reval, wave_types = shallow_water.exact_riemann_solution(q_l, q_r, grav=g)\n", | |
| " q = reval(x/t)\n", | |
| " primitive = shallow_water.conservative_to_primitive(q[0],q[1])\n", | |
| " \n", | |
| " fig = plt.figure(figsize=(18,6))\n", | |
| " names = shallow_water.primitive_variables\n", | |
| " axes = [0]*len(q_l)\n", | |
| " for i in range(len(q_l)):\n", | |
| " axes[i] = fig.add_subplot(1,len(q_l),i+1)\n", | |
| " q = primitive[i]\n", | |
| " plt.plot(x,q,linewidth=3)\n", | |
| " plt.title(names[i])\n", | |
| " qmax = max(q)\n", | |
| " qmin = min(q)\n", | |
| " qdiff = qmax - qmin\n", | |
| " axes[i].set_ylim((qmin-0.1*qdiff,qmax+0.1*qdiff))\n", | |
| " axes[i].set_xlim(-1,1)\n", | |
| " plt.show()\n", | |
| "\n", | |
| "interact(plot_exact_riemann_solution,\n", | |
| " h_l=widgets.FloatSlider(min=1.,max=10.,step=0.1,value=3.,description=r'$h_l$'),\n", | |
| " u_l=widgets.FloatSlider(min=-5.,max=5.,step=0.1,value=0.,description=r'$u_l$'),\n", | |
| " h_r=widgets.FloatSlider(min=1.,max=10.,step=0.1,value=1.,description=r'$h_r$'),\n", | |
| " u_r=widgets.FloatSlider(min=-5.,max=5.,step=0.1,value=0.,description=r'$u_r$'),\n", | |
| " t=(0.1,1.,0.1));" | |
| ] | |
| } | |
| ], | |
| "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" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 2 | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment