import matplotlib.pyplot as plt
import time
from IPython.display import clear_output
import numpy as np
Here is an example python script for solving the LWR traffic model:
$$q_t + u_\text{max}( q(1-q) )_x = 0. $$def traffic(use_petsc=0,solver_type='classic'):
from clawpack import riemann
if use_petsc:
import clawpack.petclaw as pyclaw
else:
from clawpack import pyclaw
#===========================================================================
# Set up solver and solver parameters
#===========================================================================
if solver_type=='sharpclaw':
solver = pyclaw.SharpClawSolver1D(riemann.traffic_1D)
else:
solver = pyclaw.ClawSolver1D(riemann.traffic_1D)
solver.num_waves = 1
solver.bc_lower[0] = pyclaw.BC.periodic
solver.bc_upper[0] = pyclaw.BC.periodic
#======================================
# Set up domain and initialize solution
#======================================
x = pyclaw.Dimension('x',-30.0,30.0,500)
domain = pyclaw.Domain(x)
num_eqn = 1
state = pyclaw.State(domain,num_eqn)
grid = state.grid
xc=grid.x.centers
# Stopped traffic at a light:
#state.q[0,:] = 0.75*(xc<0) + 0.1*(xc>0.)
# Formation of a traffic jam
state.q[0,:] = 0.25 + 0.25*np.exp(-0.01*xc**2)
state.problem_data['efix']=True
state.problem_data['umax']=1.
#===========================================================================
# Setup controller and controller parameters. Then solve the problem
#===========================================================================
claw = pyclaw.Controller()
claw.tfinal = 50.0
claw.solution = pyclaw.Solution(state,domain)
claw.solver = solver
claw.num_output_times = 25
claw.keep_copy = True
status = claw.run()
return claw
claw = traffic()
%pylab inline
from matplotlib import animation
import matplotlib.pyplot as plt
from clawpack.visclaw.JSAnimation import IPython_display
xc=claw.solution.state.grid.x.centers
fig = plt.figure(figsize=[8,4])
ax = plt.axes(xlim=(xc[0], xc[-1]), ylim=(0, 1.0))
line, = ax.plot([], [], lw=1)
def fplot(i):
frame = claw.frames[i]
q = frame.q[0,:]
line.set_data(xc, q)
return line,
animation.FuncAnimation(fig, fplot, frames=len(claw.frames))