# Simulating shallow water waves over a corrugated bottom with PyClaw¶

Author: David I. Ketcheson

What happens when shallow water waves run over a periodic bottom with alternating channels and ridges? Something interesting, as it turns out. Let's investigate.

To run this notebook, you will need Clawpack version >= 5.4.1.

In [1]:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt


## Problem setup¶

The code below sets this problem up. A few things to notice:

• The channels and ridges run along the $x$-direction
• We start with an initial perturbation in the surface height that is homogeneous in $y$, so it's clear that the initial propagation will be primarily in the $x$-direction, too.
• Since the bathymetry and initial condition are periodic in $y$, we save time by taking just one period as the problem domain and imposing periodic boundary conditions in $y$.
• We also save time by switching to periodic boundary boundary conditions in $x$ once the initial disturbance has moved toward the middle of the domain. This allows us to continue propagating the waves for a long time without using an extremely large domain (in $x$).
In [6]:
def qinit(state):
"Gaussian surface perturbation"
x0=0.
y0=0.

b = state.aux[0,:,:] # Bathymetry

X,Y = state.grid.p_centers
surface = ambient_surface_height+pulse_amplitude*np.exp(-X**2/pulse_width)
state.q[0,:,:] = surface - b
state.q[1,:,:] = 0.
state.q[2,:,:] = 0.

def bathymetry(y):
b_A = 0.
b_B = 0.5
return (y>0)*b_A + (y<=0)*b_B

def switch_to_periodic_BCs(solver,state):
from clawpack import pyclaw
#Change to periodic BCs after initial pulse
if state.t>periodic_bc_time and solver.bc_lower[0]==pyclaw.BC.wall:
solver.bc_lower[0]=pyclaw.BC.periodic
solver.bc_upper[0]=pyclaw.BC.periodic
solver.aux_bc_lower[0]=pyclaw.BC.periodic
solver.aux_bc_upper[0]=pyclaw.BC.periodic

def setup(cells_per_period=20,tfinal=300,solver_type='classic'):

from clawpack import riemann
from clawpack import pyclaw

if solver_type == 'classic':
solver = pyclaw.ClawSolver2D(riemann.sw_aug_2D)
solver.limiters = pyclaw.limiters.tvd.minmod
solver.dimensional_split=False
solver.cfl_max     = 0.45
solver.cfl_desired = 0.4
elif solver_type == 'sharpclaw':
solver = pyclaw.SharpClawSolver2D(riemann.sw_aug_2D)

solver.bc_lower[0] = pyclaw.BC.wall
solver.bc_upper[0] = pyclaw.BC.extrap
solver.bc_lower[1] = pyclaw.BC.periodic
solver.bc_upper[1] = pyclaw.BC.periodic

solver.aux_bc_lower[0] = pyclaw.BC.wall
solver.aux_bc_upper[0] = pyclaw.BC.extrap
solver.aux_bc_lower[1] = pyclaw.BC.periodic
solver.aux_bc_upper[1] = pyclaw.BC.periodic

solver.fwave = True
solver.before_step = switch_to_periodic_BCs

# Domain:
xlower =   0.;  xupper =  100.
ylower = -0.5;  yupper =  0.5

mx = (xupper-xlower)*cells_per_period
my = (yupper-ylower)*cells_per_period

x = pyclaw.Dimension(xlower,xupper,mx,name='x')
y = pyclaw.Dimension(ylower,yupper,my,name='y')
domain = pyclaw.Domain([x,y])

num_aux = 1
state = pyclaw.State(domain,solver.num_eqn,num_aux)
state.aux[:,:,:] = bathymetry(state.p_centers[1])

grav = 1.0 # Parameter (global auxiliary variable)
state.problem_data['grav'] = grav

qinit(state)

#===========================================================================
# Set up controller and controller parameters
#===========================================================================
claw = pyclaw.Controller()
claw.tfinal = 300.
claw.solution = pyclaw.Solution(state,domain)
claw.solver = solver
claw.num_output_times = 200
claw.keep_copy = True

return claw


Next we set the key parameters and call the problem setup function. You'll probably want to run the first setup line below, which will run quickly. But in order to fully observe the interesting behavior that appears, you can come back and run the more highly resolved, longer version (just uncomment the last line below).

In [7]:
periodic_bc_time        = 50.
ambient_surface_height  = 0.75
pulse_amplitude         = 0.05
pulse_width             = 10.

#This will run in a couple of minutes:
#n = 10; T = 50

#This will run in about 20 minutes:
n = 20; T = 300
claw = setup(cells_per_period=n,tfinal=T)

In [ ]:
claw.run()


# Plotting the results¶

The code below is used to plot individual frames of the solution.

In [9]:
def surface_height(current_data):
h = current_data.q[0,:,:]
b = bathymetry(current_data.y)
return h+b

#--------------------------
def setplot(plotdata):
#--------------------------
"""
Specify what is to be plotted at each frame.
Input:  plotdata, an instance of visclaw.data.ClawPlotData.
Output: a modified version of plotdata.
"""
from clawpack.visclaw import colormaps

plotdata.clearfigures()  # clear any old figures,axes,items data

# Figure for q[0]
plotfigure = plotdata.new_plotfigure(name='Water height', figno=0)

# Set up for axes in this figure:
plotaxes = plotfigure.new_plotaxes()
plotaxes.title = 'Water height'
plotaxes.scaled = False

# Set up for item on these axes:
plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor')
plotitem.plot_var = surface_height
plotitem.pcolor_cmap = colormaps.red_yellow_blue
plotitem.pcolor_cmin = 0.5
plotitem.pcolor_cmax = 1.5

# Figure for q[1]
plotfigure = plotdata.new_plotfigure(name='Momentum in x direction', figno=2)

# Set up for axes in this figure:
plotaxes = plotfigure.new_plotaxes()
plotaxes.title = 'Momentum in x direction'
plotaxes.scaled = False

# Set up for item on these axes:
plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor')
plotitem.plot_var = 1
plotitem.pcolor_cmap = colormaps.yellow_red_blue

# Figure for q[2]
plotfigure = plotdata.new_plotfigure(name='Momentum in y direction', figno=3)

# Set up for axes in this figure:
plotaxes = plotfigure.new_plotaxes()
plotaxes.title = 'Momentum in y direction'
plotaxes.scaled = False

# Set up for item on these axes:
plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor')
plotitem.plot_var = 2
plotitem.pcolor_cmap = colormaps.yellow_red_blue

return plotdata

claw.setplot = setplot


## Animation of surface height¶

The next two blocks of code are used to set up an animation of the results.

In [11]:
from matplotlib import animation
from IPython.display import HTML

fig = plt.figure(figsize=[12,4])

frame = claw.frames[0]
b = frame.aux[0,:,:]
h = frame.q[0,:,:]
surface = h+b

x, y = frame.state.grid.p_centers

ax2 = ax1.twinx()
slice = int(round(n/4))
line, = ax2.plot(x[:,0],surface[:,slice],'-k',linewidth=3)

im = ax1.imshow(surface.T, cmap='Blues',aspect=35,
extent=[x.min(), x.max(), y.min(), y.max()],
#extent=[x.min(), x.max(), y.min(), y.max()],
interpolation='nearest', origin='lower')

def fplot(frame_number):
frame = claw.frames[frame_number]
b = frame.aux[0,:,:]
h = frame.q[0,:,:]
surface = h+b
im.set_data(surface.T)
line.set_data(x[:,0],surface[:,slice])
return im,

anim = animation.FuncAnimation(fig, fplot, frames=len(claw.frames), interval=20, repeat=False)
plt.close()
HTML(anim.to_jshtml())

2018-05-31 10:01:59,485 INFO CLAW: Animation.save using <class 'matplotlib.animation.HTMLWriter'>

Out[11]:

Once Loop Reflect