{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Shallow water waves breaking on a beach\n", "\n", "In this example we examine the water waves breaking on a beach. We solve the shallow water equations:\n", "\n", "\\begin{align}\n", "(h)_t + (hu)_x & = 0 \\\\\n", "(hu)_t + (hu^2 + \\frac{1}{2}gh^2)_x & = -ghb_x\n", "\\end{align}\n", "\n", "Here $h$ is the water depth, $u$ is the water velocity, $g$ is a constant representing the force of gravity, and $b$ is the height of the bottom (referred to as \"bathymetry\").\n", "This problem is surprisingly challenging, for two reasons. First, if the surface height $h+b$ and the velocity $u$ are constant in space, then they should remain so for all time. But if $b$ is not constant, then this requires that terms on the left and right side of the momentum equation cancel exactly. Numerical discretizations that achieve this are said to be \"well-balanced\". Second, on the dry beach the depth $h$ is zero, and the location of the wet-dry interface moves as waves approach. Even the smallest numerical errors can lead to a negative depth near that interface.\n", "\n", "For this example we make use of [a special Riemann solver developed by David George](http://www.bu.edu/pasi-tsunami/files/2012/11/George2008.pdf) that is well-balanced and handles dry states without generating negative depths. Although the problem we're interested in is 1D, the solver is written for 2D flows. We therefore use a very narrow 2D domain." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def wave_maker_bc(state,dim,t,qbc,auxbc,num_ghost):\n", " \"Generate waves at left boundary as if there were a moving wall there.\"\n", " if dim.on_lower_boundary:\n", " qbc[0,:num_ghost,:]=qbc[0,num_ghost,:] \n", " t=state.t;\n", " amp = state.problem_data['amp'];\n", " if t<=state.problem_data['t1']: \n", " vwall = amp*(np.sin(t*np.pi/1.5))\n", " else: \n", " vwall=0.\n", " for ibc in range(num_ghost-1):\n", " qbc[1,num_ghost-ibc-1,:] = 2*vwall - qbc[1,num_ghost+ibc,:]\n", "\n", "def qinit(state):\n", " \"Gaussian surface perturbation\"\n", " b = state.aux[0,:,:] # Bathymetry\n", " X,Y = state.grid.p_centers\n", " xleft = X.min()\n", " surface = ambient_surface_height + pulse_amplitude * np.exp(-(X-(xleft+2.))**2/pulse_width)\n", " state.q[0,:,:] = np.maximum(0.0, surface - b)\n", " state.q[1,:,:] = 0.0\n", " state.q[2,:,:] = 0.0\n", " \n", "def bathymetry(x):\n", " \"Flat bottom for x<3; then a steep slope to x=5, followed by a gentle slope.\"\n", " return (x>3)*(x<5)*((x-3)*0.4) + (x>=5)*(0.8+(x-5)/40.)\n", " \n", "def setup(num_cells=500,tfinal=30,solver_type='classic',num_output_times=150):\n", "\n", " from clawpack import riemann\n", " from clawpack import pyclaw\n", "\n", " if solver_type == 'classic':\n", " solver = pyclaw.ClawSolver2D(riemann.sw_aug_2D)\n", " solver.dimensional_split=True\n", " solver.limiters = pyclaw.limiters.tvd.minmod\n", " solver.cfl_max = 0.45\n", " solver.cfl_desired = 0.4\n", " elif solver_type == 'sharpclaw':\n", " solver = pyclaw.SharpClawSolver2D(riemann.sw_aug_2D)\n", "\n", " solver.bc_lower[0] = pyclaw.BC.custom \n", " solver.user_bc_lower = wave_maker_bc\n", " solver.bc_upper[0] = pyclaw.BC.extrap\n", " solver.bc_lower[1] = pyclaw.BC.periodic\n", " solver.bc_upper[1] = pyclaw.BC.periodic\n", "\n", " solver.aux_bc_lower[0] = pyclaw.BC.extrap\n", " solver.aux_bc_upper[0] = pyclaw.BC.extrap\n", " solver.aux_bc_lower[1] = pyclaw.BC.periodic\n", " solver.aux_bc_upper[1] = pyclaw.BC.periodic\n", "\n", " solver.fwave = True\n", "\n", " # Domain:\n", " xlower, xupper = -15.0, 15.0\n", " ylower, yupper = -0.5, 0.5\n", "\n", " mx = num_cells\n", " my = 2\n", "\n", " x = pyclaw.Dimension(xlower,xupper,mx,name='x')\n", " y = pyclaw.Dimension(ylower,yupper,my,name='y')\n", " domain = pyclaw.Domain([x,y])\n", "\n", " num_aux = 1\n", " state = pyclaw.State(domain,solver.num_eqn,num_aux)\n", " state.aux[:,:,:] = bathymetry(state.p_centers[0])\n", "\n", " state.problem_data['grav'] = 10.0 # Gravitational force\n", " state.problem_data['t1'] = 50.0 # Stop generating waves after this time\n", " state.problem_data['amp'] = 0.1 # Amplitude of incoming waves\n", " qinit(state)\n", "\n", " #===========================================================================\n", " # Set up controller and controller parameters\n", " #===========================================================================\n", " claw = pyclaw.Controller()\n", " claw.tfinal = tfinal\n", " claw.solution = pyclaw.Solution(state,domain)\n", " claw.solver = solver\n", " claw.num_output_times = num_output_times\n", " claw.keep_copy = True\n", " claw.output_format = None\n", "\n", " return claw" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we set the key parameters and call the problem setup function. Running the code should take no more than a couple of minutes. With only 250 cells, there is significant numerical dissipation; if you wish you can run it with more cells to get a more accurate solution." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# These parameters are used in qinit.\n", "ambient_surface_height = 1.0\n", "pulse_amplitude = 0.0 # Use this to add an initial Gaussian wave\n", "pulse_width = 1.0\n", "\n", "claw = setup(num_cells=250,tfinal=30.0)\n", "claw.verbosity = 0 # Use this to suppress output during the run\n", "claw.run()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Plotting the results\n", "\n", "The code below is used to plot individual frames of the solution." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from matplotlib import animation\n", "from IPython.display import HTML\n", "\n", "def plot_waves(claw,ylim=(0,1.2),save_plots=False):\n", " fig = plt.figure(figsize=[12,4])\n", " ax1 = fig.add_subplot(111)\n", " fills = []\n", " frame = claw.frames[0]\n", " b = frame.aux[0,:,:]\n", " h = frame.q[0,:,:]\n", " surface = np.maximum(b,h+b)\n", "\n", " x, y = frame.state.grid.p_centers \n", " slice = 1\n", " #line, = ax1.plot(x[:,0],surface[:,slice],'-k',linewidth=3)\n", " fill = ax1.fill_between(x[:,0],b[:,slice],surface[:,slice],facecolor='blue')\n", " fill2 = ax1.fill_between(x[:,0],0*b[:,slice],b[:,slice],facecolor='brown')\n", " fills = [fill,fill2]\n", " ax1.set_xlim(-15,15)\n", " if ylim: ax1.set_ylim(ylim)\n", "\n", " def fplot(frame_number):\n", " fills[-2].remove()\n", " fills[-1].remove()\n", " frame = claw.frames[frame_number]\n", " b = frame.aux[0,:,:]\n", " h = frame.q[0,:,:]\n", " surface = np.maximum(b,h+b)\n", " #line.set_data(x[:,0],surface[:,slice])\n", " fill = ax1.fill_between(x[:,0],b[:,slice],surface[:,slice],facecolor='blue',\n", " where=b[:,slice]