Notebook

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]:

The animation above shows the water surface height perturbation (in blue). The line plot is a slice through the domain, about one-fourth of the way up from the bottom.

If you ran the low-res version, then you see that the initial disturbance starts to develop multiple peaks. After that, numerical dissipation takes over and destroys the physical behavior. So go back and run the high-res version now. You'll want to do something else while you wait for it to finish.

Animation of y-velocity¶

What causes the solitary wave formation observed above? Remember that the wave speed for shallow water waves is $u\pm \sqrt{gh}$ , where $h$ is the distance from the bottom to the surface. So the parts of the wave propagating over the deep channels move faster than those passing over the ridges. The result? Diffraction. The code below plots the $y$ -component of velocity, revealing the transverse motions due to diffraction. This microscopic diffraction has a macroscopic effect like that of dispersion, and is responsible for the solitary wave formation.

These papers give a more detailed mathematical explanation:

In [12]:

from matplotlib import animation
from clawpack.visclaw.JSAnimation import IPython_display

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

frame = claw.frames[0]
v = frame.q[2,:,:]

x, y = frame.state.grid.p_centers    

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

im = ax1.imshow(v.T, cmap='RdBu',aspect=35,vmin=-0.005, vmax=0.005,
           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]
    v = frame.q[2,:,:]
    im.set_data(v.T)
    line.set_data(x[:,0],v[:,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:02:44,619 INFO CLAW: Animation.save using <class 'matplotlib.animation.HTMLWriter'>

Out[12]: