Simulating shallow water waves over a corrugated bottom with PyClaw

Author: David I. Ketcheson
License: CC-BY

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
    plotitem.add_colorbar = True
    
    # 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
    plotitem.add_colorbar = True
    

    # 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
    plotitem.add_colorbar = True
    
    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])
ax1 = fig.add_subplot(111)

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