PyClaw is a solver for hyperbolic PDEs, based on Clawpack. You can read more about PyClaw in this paper (free version here.
In this notebook, we explore some basic PyClaw functionality. Before running the notebook, you should install Clawpack. The quick way is to just
pip install clawpack
%matplotlib inline
import matplotlib
from clawpack import pyclaw
from clawpack import riemann
Populating the interactive namespace from numpy and matplotlib
To solve a problem, we'll need to create the following:
Let's start by creating a controller and specifying the simulation end time:
claw = pyclaw.Controller()
claw.tfinal = 1.0
claw.keep_copy = True # Keep solution data in memory for plotting
claw.output_format = None # Don't write solution data to file
claw.num_output_times = 50 # Write 50 output frames
Like many solvers for nonlinear hyperbolic PDEs, PyClaw uses Riemann solvers. By specifying a Riemann solver, we will specify the system of PDEs that we want to solve.
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 Fortran, wrapped with f2py.
Note that this won't work if you're viewing the notebook online as HTML; you need to actually be running it.
riemann.
We'll solve the one-dimensional acoustics equations: $$\begin{align} p_t + K u_x & = 0 \\ u_t + \frac{1}{\rho} p_x & = 0. \end{align}$$ Here $p, u$ are the pressure and velocity as functions of $x,t$, while $\rho, K$ are constants representing the density and bulk modulus of the material transmitting the waves. We'll specify these constants later.
We can do this using the first solver in the list. Notice that the solver we create here belongs to the controller that we created above.
riemann_solver = riemann.acoustics_1D
claw.solver = pyclaw.ClawSolver1D(riemann_solver)
We also need to specify boundary conditions. We'll use periodic BCs, so that waves that go off one side of the domain come back in at the other:
claw.solver.all_bcs = pyclaw.BC.periodic
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:
domain = pyclaw.Domain( (0.,), (1.,), (100,))
Next we create a solution object that belongs to the controller and extends over the domain we specified:
claw.solution = pyclaw.Solution(claw.solver.num_eqn,domain)
The initial data is specified in an array named $q$. The pressure is contained in q[0,:]
and the velocity in q[1,:]
. We'll specify a wavepacket for the pressure and zero velocity.
x=domain.grid.x.centers
bet=100; gam=5; x0=0.75
claw.solution.q[0,:] = np.exp(-bet * (x-x0)**2) * np.cos(gam * (x - x0))
claw.solution.q[1,:] = 0.
plt.plot(x, claw.solution.q[0,:],'-o')
[<matplotlib.lines.Line2D at 0x1093166d0>]
The Riemann solver we've chosen requires some physical parameters to be specified. Press 'Tab' in the box below and you'll see what they are.
riemann_solver.cparam.
Two of these parameters are $\rho$ and $K$ in the equations above. The other two are the impedance $Z = \sqrt{\rho K}$ and sound speed $c = \sqrt{K/\rho}$. We specify these parameters in a dictionary that belongs to the solution object:
import numpy as np
density = 1.0
bulk_modulus = 1.0
impedance = np.sqrt(density*bulk_modulus)
sound_speed = np.sqrt(density/bulk_modulus)
claw.solution.state.problem_data = {
'rho' : density,
'bulk': bulk_modulus,
'zz' : np.sqrt(density*bulk_modulus),
'cc' : np.sqrt(bulk_modulus/density)
}
Finally, let's run the simulation.
status = claw.run()
2014-01-09 14:50:12,569 INFO CLAW: Solution 0 computed for time t=0.000000 2014-01-09 14:50:12,572 INFO CLAW: Solution 1 computed for time t=0.020000 2014-01-09 14:50:12,574 INFO CLAW: Solution 2 computed for time t=0.040000 2014-01-09 14:50:12,576 INFO CLAW: Solution 3 computed for time t=0.060000 2014-01-09 14:50:12,578 INFO CLAW: Solution 4 computed for time t=0.080000 2014-01-09 14:50:12,579 INFO CLAW: Solution 5 computed for time t=0.100000 2014-01-09 14:50:12,581 INFO CLAW: Solution 6 computed for time t=0.120000 2014-01-09 14:50:12,583 INFO CLAW: Solution 7 computed for time t=0.140000 2014-01-09 14:50:12,585 INFO CLAW: Solution 8 computed for time t=0.160000 2014-01-09 14:50:12,587 INFO CLAW: Solution 9 computed for time t=0.180000 2014-01-09 14:50:12,589 INFO CLAW: Solution 10 computed for time t=0.200000 2014-01-09 14:50:12,591 INFO CLAW: Solution 11 computed for time t=0.220000 2014-01-09 14:50:12,593 INFO CLAW: Solution 12 computed for time t=0.240000 2014-01-09 14:50:12,595 INFO CLAW: Solution 13 computed for time t=0.260000 2014-01-09 14:50:12,597 INFO CLAW: Solution 14 computed for time t=0.280000 2014-01-09 14:50:12,599 INFO CLAW: Solution 15 computed for time t=0.300000 2014-01-09 14:50:12,601 INFO CLAW: Solution 16 computed for time t=0.320000 2014-01-09 14:50:12,604 INFO CLAW: Solution 17 computed for time t=0.340000 2014-01-09 14:50:12,606 INFO CLAW: Solution 18 computed for time t=0.360000 2014-01-09 14:50:12,608 INFO CLAW: Solution 19 computed for time t=0.380000 2014-01-09 14:50:12,610 INFO CLAW: Solution 20 computed for time t=0.400000 2014-01-09 14:50:12,612 INFO CLAW: Solution 21 computed for time t=0.420000 2014-01-09 14:50:12,614 INFO CLAW: Solution 22 computed for time t=0.440000 2014-01-09 14:50:12,616 INFO CLAW: Solution 23 computed for time t=0.460000 2014-01-09 14:50:12,618 INFO CLAW: Solution 24 computed for time t=0.480000 2014-01-09 14:50:12,620 INFO CLAW: Solution 25 computed for time t=0.500000 2014-01-09 14:50:12,621 INFO CLAW: Solution 26 computed for time t=0.520000 2014-01-09 14:50:12,623 INFO CLAW: Solution 27 computed for time t=0.540000 2014-01-09 14:50:12,625 INFO CLAW: Solution 28 computed for time t=0.560000 2014-01-09 14:50:12,627 INFO CLAW: Solution 29 computed for time t=0.580000 2014-01-09 14:50:12,629 INFO CLAW: Solution 30 computed for time t=0.600000 2014-01-09 14:50:12,632 INFO CLAW: Solution 31 computed for time t=0.620000 2014-01-09 14:50:12,634 INFO CLAW: Solution 32 computed for time t=0.640000 2014-01-09 14:50:12,636 INFO CLAW: Solution 33 computed for time t=0.660000 2014-01-09 14:50:12,638 INFO CLAW: Solution 34 computed for time t=0.680000 2014-01-09 14:50:12,640 INFO CLAW: Solution 35 computed for time t=0.700000 2014-01-09 14:50:12,642 INFO CLAW: Solution 36 computed for time t=0.720000 2014-01-09 14:50:12,644 INFO CLAW: Solution 37 computed for time t=0.740000 2014-01-09 14:50:12,646 INFO CLAW: Solution 38 computed for time t=0.760000 2014-01-09 14:50:12,648 INFO CLAW: Solution 39 computed for time t=0.780000 2014-01-09 14:50:12,650 INFO CLAW: Solution 40 computed for time t=0.800000 2014-01-09 14:50:12,652 INFO CLAW: Solution 41 computed for time t=0.820000 2014-01-09 14:50:12,654 INFO CLAW: Solution 42 computed for time t=0.840000 2014-01-09 14:50:12,656 INFO CLAW: Solution 43 computed for time t=0.860000 2014-01-09 14:50:12,658 INFO CLAW: Solution 44 computed for time t=0.880000 2014-01-09 14:50:12,660 INFO CLAW: Solution 45 computed for time t=0.900000 2014-01-09 14:50:12,662 INFO CLAW: Solution 46 computed for time t=0.920000 2014-01-09 14:50:12,664 INFO CLAW: Solution 47 computed for time t=0.940000 2014-01-09 14:50:12,666 INFO CLAW: Solution 48 computed for time t=0.960000 2014-01-09 14:50:12,668 INFO CLAW: Solution 49 computed for time t=0.980000 2014-01-09 14:50:12,670 INFO CLAW: Solution 50 computed for time t=1.000000
Now we'll plot the results, which are contained in claw.frames[:]
. It's simple to plot a single frame with matplotlib:
pressure = claw.frames[50].q[0,:]
plt.plot(x,pressure,'-o')
[<matplotlib.lines.Line2D at 0x10935a0d0>]
To examine the evolution more thoroughly, it's nice to see all the frames in sequence. We can do this as follows.
from matplotlib import animation
import matplotlib.pyplot as plt
from clawpack.visclaw.JSAnimation import IPython_display
import numpy as np
fig = plt.figure()
ax = plt.axes(xlim=(0, 1), ylim=(-0.2, 1.2))
frame = claw.frames[0]
pressure = frame.q[0,:]
line, = ax.plot([], [], lw=2)
def fplot(frame_number):
frame = claw.frames[frame_number]
pressure = frame.q[0,:]
line.set_data(x,pressure)
return line,
animation.FuncAnimation(fig, fplot, frames=len(claw.frames), interval=30)
That's it! Here are some things you might try for fun:
pyclaw.BC.
+[Tab] to get a list of boundary conditions available).SharpClawSolver1D
instead of a ClawSolver1D
. How does this affect the final solution? You can read more about the methods in SharpClaw in this paper.