{
"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]<surface[:,slice])\n",
" fill2 = ax1.fill_between(x[:,0],0*b[:,slice],b[:,slice],facecolor='brown')\n",
" fills.append(fill)\n",
" fills.append(fill2)\n",
" if save_plots:\n",
" fname = 'frame'+str(frame_number).zfill(4)+'.png'\n",
" fig.savefig(fname) \n",
" return fill,\n",
"\n",
" anim = animation.FuncAnimation(fig, fplot, frames=len(claw.frames), interval=100, repeat=False)\n",
" plt.close()\n",
" return HTML(anim.to_jshtml())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plot full view"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plot_waves(claw,ylim=(0,1.5))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Zoomed-in view\n",
"The next animation is zoomed in in the vertical dimension, to better show the structure of the waves and the runup on the beach."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"plot_waves(claw,ylim=(0.9,1.25))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The characteristic speeds for the shallow water equations are $u \\pm \\sqrt{gh}$. Thus waves travel slower in shallower water.\n",
"\n",
"In deep water, we see that the waves travel almost without changing shape. As the waves reach the shallow sloping beach, two things happen. First, their wavelength decreases because their speed decreases. Second, they steepen because the difference in depth between the crest and trough becomes significant and the crest catches up to the trough.\n",
"\n",
"Unlike real water waves, our waves cannot break (overturn) because $h$ is required to be a single-valued function of $x$. Instead they form discontinuous shock waves."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.0"
},
"toc": {
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"toc_cell": false,
"toc_position": {},
"toc_section_display": "block",
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 4
}