A quick introduction to PyClaw¶
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
If you are viewing this notebook online, you can download it and run it from here.
In [2]:
%matplotlib inline
import matplotlib
from clawpack import pyclaw
from clawpack import riemann
Populating the interactive namespace from numpy and matplotlib
Setting up a problem¶
To solve a problem, we’ll need to create the following:
- A controller, which handles the running, output, and can be used for plotting (you don’t absolutely need a controller, but it makes life simpler)
- 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.
- A domain over which to solve the problem
- A solution, where we will provide the initial data. After running, the solution will contain – you guessed it! – the solution.
Let’s start by creating a controller and specifying the simulation end time:
In [3]:
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
Riemann solvers¶
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.
In [ ]:
riemann.
We’ll solve the one-dimensional acoustics equations:
\(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.
In [4]:
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:
In [5]:
claw.solver.all_bcs = pyclaw.BC.periodic
The problem domain¶
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:
In [6]:
domain = pyclaw.Domain( (0.,), (1.,), (100,))
The initial solution¶
Next we create a solution object that belongs to the controller and extends over the domain we specified:
In [7]:
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.
In [9]:
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')
Out[9]:
[<matplotlib.lines.Line2D at 0x1093166d0>]
Problem-specific parameters¶
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.
In [ ]:
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:
In [10]:
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.
In [11]:
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
Plotting¶
Now we’ll plot the results, which are contained in claw.frames[:].
It’s simple to plot a single frame with matplotlib:
In [12]:
pressure = claw.frames[50].q[0,:]
plt.plot(x,pressure,'-o')
Out[12]:
[<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.
In [13]:
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)
Out[13]:
That’s it! Here are some things you might try for fun:
- Change the boundary conditions to reflecting or outflow (hint: type
pyclaw.BC.+[Tab] to get a list of boundary conditions available). - Change the grid to use a larger or smaller number of grid cells. How does this affect the final solution?
- Use higher-order methods by instantiating a
SharpClawSolver1Dinstead of aClawSolver1D. How does this affect the final solution? You can read more about the methods in SharpClaw in this paper.
In [ ]:
