{
"metadata": {
"name": "",
"signature": "sha256:a6012fdd078ee5922bad6e7a84c31a42690e4803260228bcf9ac2e37d98e9742"
},
"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": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n"
],
"metadata": {},
"output_type": "pyout",
"prompt_number": 1,
"text": [
""
]
}
],
"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": [
"Solving problems with PyClaw"
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Variable coefficients"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Many hyperbolic problems involve coefficients that vary in space and/or time. The simplest example is our old friend the advection equation, but where the velocity may vary:\n",
"\n",
"$$q_t + (a(x,t) q)_x = 0.$$\n",
"\n",
"Such problems can be solved with PyClaw by storing the variable coefficients the `state.aux` array. For coefficients that vary in time, you can provide a function that updates the `aux` array at each time step."
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Non-hyperbolic terms"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Many physical models include both hyperbolic terms that represent transport and other terms that may represent reactions, diffusion, external forces, etc. It is common to refer to such terms as source terms, denoted by $S(q)$ and to the resulting equation as a *balance law*:\n",
"\n",
"$$q_t + f(q)_x = S(q).$$\n",
"\n",
"Such problems can be solved with PyClaw, but you must provide a routine that integrates the source term over one time step. Alternatively, the source terms may be incorporated into the Riemann solver by using an *f-wave solver*; this is particularly useful in situations where $f(q)_x$ and $S(q)$ are nearly equal."
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Boundary conditions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"PyClaw has built-in functions for the following kinds of boundary conditions:\n",
"- Periodic\n",
"- Non-reflecting (zero-order extrapolation)\n",
"- Reflecting\n",
"\n",
"For other types of boundary conditions, you can provide a function that fills the ghost cells at the boundary of the domain."
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Example: Shock-bubble interaction"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's look at an example that includes variable coefficients, source terms, and a custom boundary condition. This example appeared in a [2013 paper in SISC](http://epubs.siam.org/doi/abs/10.1137/110856976). It models the interaction of a shock wave with a low-density bubble. The problem is 3-dimensional, but we'll assume the shock is planar and the bubble is a sphere, which allows us to reduce the problem to two dimensions using cylindrical symmetry. The resulting equations are:\n",
"\n",
"$$\n",
"\\rho_t + (\\rho u)_z + (\\rho v)_r = - \\frac{\\rho v}{r} \\\\\n",
"(\\rho u)_t + (\\rho u^2 + p)_z + (\\rho u v)_r = - \\frac{\\rho u v}{r} \\\\\n",
"(\\rho v)_t + (\\rho u v)_z + (\\rho v^2 + p)_r = -\\frac{\\rho v^2}{r} \\\\\n",
"(\\rho E)_t + ((E+p)u)_z + ((E+p)v)_r = - \\frac{(E+p)v}{r}.\n",
"$$\n",
"\n",
"The left hand side of these equations is identical to the 2D Euler equations in cartesian geometry; the right hand side terms arise because of the difference between $\\nabla$ in cartesian and cylindrical coordinates. We refer to those as *geometric source terms*. Here is a function that integrates just those source terms using a second-order Runge-Kutta method."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"gamma = 1.4\n",
"gamma1 = gamma - 1.\n",
"\n",
"def step_Euler_radial(solver,state,dt):\n",
" \"\"\"\n",
" Geometric source terms for Euler equations with radial symmetry.\n",
" Integrated using a 2-stage, 2nd-order Runge-Kutta method.\n",
" This is a Clawpack-style source term routine.\n",
" \"\"\"\n",
" \n",
" dt2 = dt/2.\n",
" ndim = 2\n",
"\n",
" aux=state.aux\n",
" q = state.q\n",
"\n",
" rad = aux[0,:,:]\n",
"\n",
" rho = q[0,:,:]\n",
" u = q[1,:,:]/rho\n",
" v = q[2,:,:]/rho\n",
" press = gamma1 * (q[3,:,:] - 0.5*rho*(u**2 + v**2))\n",
"\n",
" qstar = np.empty(q.shape)\n",
"\n",
" qstar[0,:,:] = q[0,:,:] - dt2*(ndim-1)/rad * q[2,:,:]\n",
" qstar[1,:,:] = q[1,:,:] - dt2*(ndim-1)/rad * rho*u*v\n",
" qstar[2,:,:] = q[2,:,:] - dt2*(ndim-1)/rad * rho*v*v\n",
" qstar[3,:,:] = q[3,:,:] - dt2*(ndim-1)/rad * v * (q[3,:,:] + press)\n",
"\n",
" rho = qstar[0,:,:]\n",
" u = qstar[1,:,:]/rho\n",
" v = qstar[2,:,:]/rho\n",
" press = gamma1 * (qstar[3,:,:] - 0.5*rho*(u**2 + v**2))\n",
"\n",
" q[0,:,:] = q[0,:,:] - dt*(ndim-1)/rad * qstar[2,:,:]\n",
" q[1,:,:] = q[1,:,:] - dt*(ndim-1)/rad * rho*u*v\n",
" q[2,:,:] = q[2,:,:] - dt*(ndim-1)/rad * rho*v*v\n",
" q[3,:,:] = q[3,:,:] - dt*(ndim-1)/rad * v * (qstar[3,:,:] + press)"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "heading",
"level": 3,
"metadata": {},
"source": [
"Initial condition"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here is the code to set up the initial condition. Note that we have simplified it a bit: cells that lie on the boundary of the bubble should properly be initialized with a weighted average of states, but here we just initialize based on whether the cell center is inside the bubble or not. For the full code that does weighted averaging, see the corresponding example in the PyClaw examples directory."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def qinit(state,z0=0.5,r0=0.,bubble_radius=0.2,rhoin=0.1,pinf=5.):\n",
" grid = state.grid\n",
"\n",
" # Background state\n",
" rhoout = 1.\n",
" pout = 1.\n",
" \n",
" # Bubble state\n",
" pin = 1.\n",
"\n",
" # Shock state\n",
" zshock = 0.2\n",
" rinf = (gamma1 + pinf*(gamma+1.))/ ((gamma+1.) + gamma1*pinf)\n",
" vinf = 1./np.sqrt(gamma) * (pinf - 1.) / np.sqrt(0.5*((gamma+1.)/gamma) * pinf+0.5*gamma1/gamma)\n",
" einf = 0.5*rinf*vinf**2 + pinf/gamma1\n",
" \n",
" z = grid.z.centers\n",
" r = grid.r.centers\n",
" R,Z = np.meshgrid(r,z)\n",
" rad = np.sqrt((Z-z0)**2 + (R-r0)**2)\n",
"\n",
" state.q[0,:,:] = rinf*(Zbubble_radius)*(Z>=zshock)\n",
" state.q[1,:,:] = rinf*vinf*(Zbubble_radius)*(Z>=zshock))/gamma1\n",
" state.q[4,:,:] = 1.*(rad<=bubble_radius)"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "heading",
"level": 3,
"metadata": {},
"source": [
"Boundary conditions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Symmetry dictates that we use a reflecting boundary condition at the bottom boundary, $r=0$. We will impose outflow boundary conditions at the top and right boundaries.\n",
"\n",
"Next we will write the function for the boundary condition at the left of the domain, where a shock is entering. This function operates on an array called `qbc` (for *q with boundary conditions*) and must fill in the values of the ghost cells at the left boundary of the grid."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def shockbc(state,dim,t,qbc,num_ghost):\n",
" \"\"\"\n",
" Incoming shock at left boundary.\n",
" \"\"\"\n",
" p = 5.\n",
" rho = (gamma1 + p*(gamma+1.))/ ((gamma+1.) + gamma1*p)\n",
" v = 1./np.sqrt(gamma) * (p - 1.) / np.sqrt(0.5*((gamma+1.)/gamma) * p+0.5*gamma1/gamma)\n",
" e = 0.5*rho*v**2 + p/gamma1\n",
"\n",
" for i in xrange(num_ghost):\n",
" qbc[0,i,...] = rho\n",
" qbc[1,i,...] = rho*v\n",
" qbc[2,i,...] = 0.\n",
" qbc[3,i,...] = e\n",
" qbc[4,i,...] = 0."
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now the main code to set up the problem:"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from clawpack import riemann\n",
"from clawpack import pyclaw\n",
"import numpy as np\n",
"\n",
"solver = pyclaw.ClawSolver2D(riemann.euler_5wave_2D)\n",
"solver.step_source=step_Euler_radial\n",
"\n",
"# Boundary conditions\n",
"solver.bc_lower[0]=pyclaw.BC.custom\n",
"solver.bc_upper[0]=pyclaw.BC.extrap\n",
"solver.bc_lower[1]=pyclaw.BC.wall\n",
"solver.bc_upper[1]=pyclaw.BC.extrap\n",
"solver.user_bc_lower=shockbc\n",
"\n",
"# We must also set boundary conditions for the auxiliary variable (r)\n",
"# But in this case the values don't matter\n",
"solver.aux_bc_lower[0]=pyclaw.BC.extrap\n",
"solver.aux_bc_upper[0]=pyclaw.BC.extrap\n",
"solver.aux_bc_lower[1]=pyclaw.BC.extrap\n",
"solver.aux_bc_upper[1]=pyclaw.BC.extrap\n",
"\n",
"# Initialize domain\n",
"mx=320; my=80\n",
"z = pyclaw.Dimension('z',0.0,2.0,mx)\n",
"r = pyclaw.Dimension('r',0.0,0.5,my)\n",
"domain = pyclaw.Domain([z,r])\n",
"num_aux=1\n",
"state = pyclaw.State(domain,solver.num_eqn,num_aux)\n",
"\n",
"state.problem_data['gamma']= gamma\n",
"state.problem_data['gamma1']= gamma1\n",
"\n",
"qinit(state)\n",
"\n",
"rc = state.grid.r.centers\n",
"for j,rcoord in enumerate(rc):\n",
" state.aux[0,:,j] = rcoord\n",
"\n",
"solver.cfl_max = 0.5\n",
"solver.cfl_desired = 0.45\n",
"solver.source_split = 1 # Use first-order splitting for the source term\n",
"\n",
"claw = pyclaw.Controller()\n",
"claw.keep_copy = True\n",
"claw.output_format = None\n",
"claw.tfinal = 0.6\n",
"claw.solution = pyclaw.Solution(state,domain)\n",
"claw.solver = solver\n",
"claw.num_output_times = 60"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"claw.run()"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"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(figsize=[8,4])\n",
"\n",
"frame = claw.frames[0]\n",
"rho = frame.q[0,:,:]\n",
"\n",
"x, y = frame.state.grid.c_centers \n",
"\n",
"# This essentially does a pcolor plot, but it returns the appropriate object\n",
"# for use in animation. See http://matplotlib.org/examples/pylab_examples/pcolor_demo.html.\n",
"# Note that it's necessary to transpose the data array because of the way imshow works.\n",
"im = plt.imshow(rho.T, cmap='Greys', vmin=rho.min(), vmax=rho.max(),\n",
" extent=[x.min(), x.max(), y.min(), y.max()],\n",
" interpolation='nearest', origin='lower')\n",
"\n",
"def fplot(frame_number):\n",
" frame = claw.frames[frame_number]\n",
" rho = frame.q[0,:,:]\n",
" im.set_data(rho.T)\n",
" return im,"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"animation.FuncAnimation(fig, fplot, frames=len(claw.frames), interval=20)"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Mapped grids"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So far we have solved problems on uniform cartesian grids. For many problems, the natural domain is not cartesian; in other problems, there may be internal interfaces not aligned with a cartesian grid. PyClaw can be used on any domain for which there exists a mapping to a cartesian domain. Such mappings are surprisingly flexible, and can be used to solve in disks and spheres, as well as many more complicated geometries. Solving on a mapped grid requires the use of a Riemann solver that takes the geometry into account."
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Other systems of equations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Clawpack includes Riemann solvers for the following systems of equations:\n",
"- Acoustics in 1, 2, and 3 dimensions\n",
"- Constant and variable-coefficient advection in 1 and 2 dimensions\n",
"- Burgers equation\n",
"- The Euler equations of compressible fluid dynamics in 1, 2, and 3 dimensions\n",
"- The reactive Euler equations in 1 dimension\n",
"- Nonlinear elasticity in 1 dimension\n",
"- Linear elasticity in 2 dimensions\n",
"- The shallow water equations, in 1 and 2 dimensions, including variable bathymetry\n",
"- Traffic flow in 1 dimension\n",
"\n",
"Additional solvers are being developed by users all the time; if you don't see the system that interests you in this list, it's worth asking on the [Clawpack Users Google group](claw-users@groups.google.com) to see if someone has developed a solver already. You can also write your own Riemann solver, but that is beyond the scope of this course."
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"More examples"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"I think the simplest way to learn more about using PyClaw is through examples. Here are a few more:\n",
"\n",
"- [Acoustics](http://nbviewer.ipython.org/8332861) - An introductory example similar to Lesson 4 of this course, but using the acoustics equations.\n",
"- [Quadrants](http://nbviewer.ipython.org/8333043) - Another 2-dimensional Euler equations example.\n",
"- [Traffic](http://nbviewer.ipython.org/348ea7a9ae8b4a558c5b) - Solving the LWR traffic model with PyClaw.\n",
"- [Stegotons](http://nbviewer.ipython.org/8554686) - Illustration of a curious class of solitary waves that arise in periodic nonlinear wave problems (using a nonlinear elasticity model). This example includes a custom boundary condition and a demonstration of how to use a `before_step` function.\n",
"- [Wave runup on a beach](http://nbviewer.ipython.org/github/ketch/shallow_water_periodic_bathymetry/blob/master/pyclaw/beach.ipynb) - using the 1D shallow water equations\n",
"- [Shallow water solitary waves](http://nbviewer.ipython.org/9250942) - in 2D, with periodic bathymetry"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Many additional examples are included in the PyClaw examples module."
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from clawpack.pyclaw import examples"
],
"language": "python",
"metadata": {},
"outputs": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"examples."
],
"language": "python",
"metadata": {},
"outputs": []
}
],
"metadata": {}
}
]
}