from IPython.core.display import HTML
css_file = './custom.css'
HTML(open(css_file, "r").read())
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:
$$q_t + (a(x,t) q)_x = 0.$$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.
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:
$$q_t + f(q)_x = S(q).$$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.
PyClaw has built-in functions for the following kinds of boundary conditions:
For other types of boundary conditions, you can provide a function that fills the ghost cells at the boundary of the domain.
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. 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:
$$ \rho_t + (\rho u)_z + (\rho v)_r = - \frac{\rho v}{r} \\ (\rho u)_t + (\rho u^2 + p)_z + (\rho u v)_r = - \frac{\rho u v}{r} \\ (\rho v)_t + (\rho u v)_z + (\rho v^2 + p)_r = -\frac{\rho v^2}{r} \\ (\rho E)_t + ((E+p)u)_z + ((E+p)v)_r = - \frac{(E+p)v}{r}. $$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.
gamma = 1.4
gamma1 = gamma - 1.
def step_Euler_radial(solver,state,dt):
"""
Geometric source terms for Euler equations with radial symmetry.
Integrated using a 2-stage, 2nd-order Runge-Kutta method.
This is a Clawpack-style source term routine.
"""
dt2 = dt/2.
ndim = 2
aux=state.aux
q = state.q
rad = aux[0,:,:]
rho = q[0,:,:]
u = q[1,:,:]/rho
v = q[2,:,:]/rho
press = gamma1 * (q[3,:,:] - 0.5*rho*(u**2 + v**2))
qstar = np.empty(q.shape)
qstar[0,:,:] = q[0,:,:] - dt2*(ndim-1)/rad * q[2,:,:]
qstar[1,:,:] = q[1,:,:] - dt2*(ndim-1)/rad * rho*u*v
qstar[2,:,:] = q[2,:,:] - dt2*(ndim-1)/rad * rho*v*v
qstar[3,:,:] = q[3,:,:] - dt2*(ndim-1)/rad * v * (q[3,:,:] + press)
rho = qstar[0,:,:]
u = qstar[1,:,:]/rho
v = qstar[2,:,:]/rho
press = gamma1 * (qstar[3,:,:] - 0.5*rho*(u**2 + v**2))
q[0,:,:] = q[0,:,:] - dt*(ndim-1)/rad * qstar[2,:,:]
q[1,:,:] = q[1,:,:] - dt*(ndim-1)/rad * rho*u*v
q[2,:,:] = q[2,:,:] - dt*(ndim-1)/rad * rho*v*v
q[3,:,:] = q[3,:,:] - dt*(ndim-1)/rad * v * (qstar[3,:,:] + press)
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.
def qinit(state,z0=0.5,r0=0.,bubble_radius=0.2,rhoin=0.1,pinf=5.):
grid = state.grid
# Background state
rhoout = 1.
pout = 1.
# Bubble state
pin = 1.
# Shock state
zshock = 0.2
rinf = (gamma1 + pinf*(gamma+1.))/ ((gamma+1.) + gamma1*pinf)
vinf = 1./np.sqrt(gamma) * (pinf - 1.) / np.sqrt(0.5*((gamma+1.)/gamma) * pinf+0.5*gamma1/gamma)
einf = 0.5*rinf*vinf**2 + pinf/gamma1
z = grid.z.centers
r = grid.r.centers
R,Z = np.meshgrid(r,z)
rad = np.sqrt((Z-z0)**2 + (R-r0)**2)
state.q[0,:,:] = rinf*(Z<zshock) + rhoin*(rad<=bubble_radius) + rhoout*(rad>bubble_radius)*(Z>=zshock)
state.q[1,:,:] = rinf*vinf*(Z<zshock)
state.q[2,:,:] = 0.
state.q[3,:,:] = einf*(Z<zshock) + (pin*(rad<=bubble_radius) + pout*(rad>bubble_radius)*(Z>=zshock))/gamma1
state.q[4,:,:] = 1.*(rad<=bubble_radius)
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.
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.
def shockbc(state,dim,t,qbc,num_ghost):
"""
Incoming shock at left boundary.
"""
p = 5.
rho = (gamma1 + p*(gamma+1.))/ ((gamma+1.) + gamma1*p)
v = 1./np.sqrt(gamma) * (p - 1.) / np.sqrt(0.5*((gamma+1.)/gamma) * p+0.5*gamma1/gamma)
e = 0.5*rho*v**2 + p/gamma1
for i in xrange(num_ghost):
qbc[0,i,...] = rho
qbc[1,i,...] = rho*v
qbc[2,i,...] = 0.
qbc[3,i,...] = e
qbc[4,i,...] = 0.
Now the main code to set up the problem:
from clawpack import riemann
from clawpack import pyclaw
import numpy as np
solver = pyclaw.ClawSolver2D(riemann.euler_5wave_2D)
solver.step_source=step_Euler_radial
# Boundary conditions
solver.bc_lower[0]=pyclaw.BC.custom
solver.bc_upper[0]=pyclaw.BC.extrap
solver.bc_lower[1]=pyclaw.BC.wall
solver.bc_upper[1]=pyclaw.BC.extrap
solver.user_bc_lower=shockbc
# We must also set boundary conditions for the auxiliary variable (r)
# But in this case the values don't matter
solver.aux_bc_lower[0]=pyclaw.BC.extrap
solver.aux_bc_upper[0]=pyclaw.BC.extrap
solver.aux_bc_lower[1]=pyclaw.BC.extrap
solver.aux_bc_upper[1]=pyclaw.BC.extrap
# Initialize domain
mx=320; my=80
z = pyclaw.Dimension('z',0.0,2.0,mx)
r = pyclaw.Dimension('r',0.0,0.5,my)
domain = pyclaw.Domain([z,r])
num_aux=1
state = pyclaw.State(domain,solver.num_eqn,num_aux)
state.problem_data['gamma']= gamma
state.problem_data['gamma1']= gamma1
qinit(state)
rc = state.grid.r.centers
for j,rcoord in enumerate(rc):
state.aux[0,:,j] = rcoord
solver.cfl_max = 0.5
solver.cfl_desired = 0.45
solver.source_split = 1 # Use first-order splitting for the source term
claw = pyclaw.Controller()
claw.keep_copy = True
claw.output_format = None
claw.tfinal = 0.6
claw.solution = pyclaw.Solution(state,domain)
claw.solver = solver
claw.num_output_times = 60
claw.run()
from matplotlib import animation
import matplotlib.pyplot as plt
from clawpack.visclaw.JSAnimation import IPython_display
import numpy as np
fig = plt.figure(figsize=[8,4])
frame = claw.frames[0]
rho = frame.q[0,:,:]
x, y = frame.state.grid.c_centers
# This essentially does a pcolor plot, but it returns the appropriate object
# for use in animation. See http://matplotlib.org/examples/pylab_examples/pcolor_demo.html.
# Note that it's necessary to transpose the data array because of the way imshow works.
im = plt.imshow(rho.T, cmap='Greys', vmin=rho.min(), vmax=rho.max(),
extent=[x.min(), x.max(), y.min(), y.max()],
interpolation='nearest', origin='lower')
def fplot(frame_number):
frame = claw.frames[frame_number]
rho = frame.q[0,:,:]
im.set_data(rho.T)
return im,
animation.FuncAnimation(fig, fplot, frames=len(claw.frames), interval=20)
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.
Clawpack includes Riemann solvers for the following systems of equations:
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 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.
I think the simplest way to learn more about using PyClaw is through examples. Here are a few more:
before_step
function.Many additional examples are included in the PyClaw examples module.
from clawpack.pyclaw import examples
examples.