In this example, we solve the Euler equations with initial data consisting of a different state in each quadrant of the unit square. The following block of code will compute the solution at 40 timeframes and should take 1-2 minutes to run.
from clawpack import pyclaw
from clawpack import riemann
claw = pyclaw.Controller()
claw.tfinal = 0.6
claw.num_output_times = 40
riemann_solver = riemann.euler_4wave_2D
claw.solver = pyclaw.ClawSolver2D(riemann_solver)
claw.solver.all_bcs = pyclaw.BC.extrap
grid_size = (300,300)
domain = pyclaw.Domain( (0.,0.), (1.,1.), grid_size)
claw.solution = pyclaw.Solution(claw.solver.num_eqn,domain)
gam = 1.4
claw.solution.problem_data['gamma'] = gam
# Set initial data
q = claw.solution.q
xx,yy = domain.grid.p_centers
l = xx<0.5; r = xx>=0.5; b = yy<0.5; t = yy>=0.5
q[0,...] = 2.*l*t + 1.*l*b + 1.*r*t + 3.*r*b
q[1,...] = 0.75*t - 0.75*b
q[2,...] = 0.5*l - 0.5*r
q[3,...] = 0.5*q[0,...]*(q[1,...]**2+q[2,...]**2) + 1./(gam-1.)
claw.keep_copy = True # Keep solution data in memory for plotting
claw.output_format = None # Don't write solution data to file
status = claw.run()
2014-01-09 10:32:47,173 INFO CLAW: Solution 0 computed for time t=0.000000 2014-01-09 10:32:47,875 INFO CLAW: Solution 1 computed for time t=0.015000 2014-01-09 10:32:48,594 INFO CLAW: Solution 2 computed for time t=0.030000 2014-01-09 10:32:49,328 INFO CLAW: Solution 3 computed for time t=0.045000 2014-01-09 10:32:50,020 INFO CLAW: Solution 4 computed for time t=0.060000 2014-01-09 10:32:50,720 INFO CLAW: Solution 5 computed for time t=0.075000 2014-01-09 10:32:51,432 INFO CLAW: Solution 6 computed for time t=0.090000 2014-01-09 10:32:52,149 INFO CLAW: Solution 7 computed for time t=0.105000 2014-01-09 10:32:52,879 INFO CLAW: Solution 8 computed for time t=0.120000 2014-01-09 10:32:53,614 INFO CLAW: Solution 9 computed for time t=0.135000 2014-01-09 10:32:54,359 INFO CLAW: Solution 10 computed for time t=0.150000 2014-01-09 10:32:55,115 INFO CLAW: Solution 11 computed for time t=0.165000 2014-01-09 10:32:55,878 INFO CLAW: Solution 12 computed for time t=0.180000 2014-01-09 10:32:56,644 INFO CLAW: Solution 13 computed for time t=0.195000 2014-01-09 10:32:57,414 INFO CLAW: Solution 14 computed for time t=0.210000 2014-01-09 10:32:58,189 INFO CLAW: Solution 15 computed for time t=0.225000 2014-01-09 10:32:58,972 INFO CLAW: Solution 16 computed for time t=0.240000 2014-01-09 10:32:59,764 INFO CLAW: Solution 17 computed for time t=0.255000 2014-01-09 10:33:00,552 INFO CLAW: Solution 18 computed for time t=0.270000 2014-01-09 10:33:01,346 INFO CLAW: Solution 19 computed for time t=0.285000 2014-01-09 10:33:02,136 INFO CLAW: Solution 20 computed for time t=0.300000 2014-01-09 10:33:02,939 INFO CLAW: Solution 21 computed for time t=0.315000 2014-01-09 10:33:03,743 INFO CLAW: Solution 22 computed for time t=0.330000 2014-01-09 10:33:04,551 INFO CLAW: Solution 23 computed for time t=0.345000 2014-01-09 10:33:05,362 INFO CLAW: Solution 24 computed for time t=0.360000 2014-01-09 10:33:06,158 INFO CLAW: Solution 25 computed for time t=0.375000 2014-01-09 10:33:06,959 INFO CLAW: Solution 26 computed for time t=0.390000 2014-01-09 10:33:07,756 INFO CLAW: Solution 27 computed for time t=0.405000 2014-01-09 10:33:08,557 INFO CLAW: Solution 28 computed for time t=0.420000 2014-01-09 10:33:09,369 INFO CLAW: Solution 29 computed for time t=0.435000 2014-01-09 10:33:10,182 INFO CLAW: Solution 30 computed for time t=0.450000 2014-01-09 10:33:10,988 INFO CLAW: Solution 31 computed for time t=0.465000 2014-01-09 10:33:11,793 INFO CLAW: Solution 32 computed for time t=0.480000 2014-01-09 10:33:12,604 INFO CLAW: Solution 33 computed for time t=0.495000 2014-01-09 10:33:13,411 INFO CLAW: Solution 34 computed for time t=0.510000 2014-01-09 10:33:14,225 INFO CLAW: Solution 35 computed for time t=0.525000 2014-01-09 10:33:15,032 INFO CLAW: Solution 36 computed for time t=0.540000 2014-01-09 10:33:15,838 INFO CLAW: Solution 37 computed for time t=0.555000 2014-01-09 10:33:16,643 INFO CLAW: Solution 38 computed for time t=0.570000 2014-01-09 10:33:17,454 INFO CLAW: Solution 39 computed for time t=0.585000 2014-01-09 10:33:18,264 INFO CLAW: Solution 40 computed for time t=0.600000
%pylab inline
frame = claw.frames[40]
density = frame.q[0,:,:]
(vx,vy) = np.gradient(density)
vs = np.sqrt(vx**2 + vy**2)
x, y = frame.state.grid.c_centers
plt.pcolormesh(x, y, vs, cmap='RdBu')
plt.axis('image')
Populating the interactive namespace from numpy and matplotlib
WARNING: pylab import has clobbered these variables: ['gamma'] `%pylab --no-import-all` prevents importing * from pylab and numpy
(0.0016666666666666668, 0.99833333333333341, 0.0016666666666666668, 0.99833333333333341)
Next we will plot an animation of all 40 frames, using matplotlib and Jake Vanderplas' JSAnimation code.
from matplotlib import animation
import matplotlib.pyplot as plt
from clawpack.visclaw.JSAnimation import IPython_display
import numpy as np
fig = plt.figure(figsize=[4,4])
frame = claw.frames[0]
density = frame.q[0,:,:]
(vx,vy) = np.gradient(density)
vs = np.sqrt(vx**2 + vy**2)
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(vs.T, cmap='Greys', vmin=vs.min(), vmax=vs.max()/20,
extent=[x.min(), x.max(), y.min(), y.max()],
interpolation='nearest', origin='lower')
def fplot(frame_number):
frame = claw.frames[frame_number]
density = frame.q[0,:,:]
(vx,vy) = np.gradient(density)
vs = np.sqrt(vx**2 + vy**2)
im.set_data(vs.T)
return im,
animation.FuncAnimation(fig, fplot, frames=len(claw.frames), interval=20)
Plotly is a new web service that makes it possible (among many other things) to embed interactive plots in an IPython notebook.
import plotly
un='jackp'
k='11m2qbzob9'
py = plotly.plotly(username_or_email=un, key=k)
frame = claw.frames[40]
density = frame.q[0,:,:]
(vx,vy) = np.gradient(density)
vs = np.sqrt(vx**2 + vy**2)
# Note that it's necessary to transpose the data array because plotly's heatmap uses imshow.
py.iplot({'z': vs.T, 'type':'heatmap'})
2014-01-09 10:46:43,720 INFO CLAW: Starting new HTTPS connection (1): plot.ly
To plot with VisClaw, we must first define a setplot function:
def setplot(plotdata):
plotfigure = plotdata.new_plotfigure(name='Density', figno=0)
plotaxes = plotfigure.new_plotaxes()
plotaxes.title = 'Density'
plotaxes.scaled = True # so aspect ratio is 1
plotitem = plotaxes.new_plotitem(plot_type='2d_schlieren')
plotitem.plot_var = 0
return plotdata
Now we can plot a single frame as follows (most of this code will be unnecessary in the next version):
from clawpack.visclaw import data
from clawpack.visclaw import frametools
plotdata = data.ClawPlotData()
plotdata.setplot = claw.setplot
claw.plotdata = frametools.call_setplot(claw.setplot,plotdata)
frame = claw.load_frame(40)
claw.plot_frame(frame)
Executed setplot successfully
VisClaw also has an interactive plotting loop capability. However, that doesn't work well in the notebook, because no plot appears until you exit the loop, and then you only get the last plot.
claw.plot()
Executed setplot successfully Interactive plotting... Type ? at IPLOT prompt for list of commands Start at which frame [default=0] ? Plotting Frame 0 at t = 0.0 IPLOT > Plotting Frame 1 at t = 0.015 IPLOT > Plotting Frame 2 at t = 0.03 IPLOT > Plotting Frame 3 at t = 0.045 IPLOT > q quitting...