# 2-dimensional shallow water equations¶

## 2D shallow water: radial dam break¶

Solve the 2D shallow water equations:

$\begin{split}h_t + (hu)_x + (hv)_y = 0 \\ (hu)_t + (hu^2 + \frac{1}{2}gh^2)_x + (huv)_y = 0 \\ (hv)_t + (huv)_x + (hv^2 + \frac{1}{2}gh^2)_y = 0.\end{split}$

The initial condition is a circular area with high depth surrounded by lower-depth water. The top and right boundary conditions reflect, while the bottom and left boundaries are outflow.

## Source:¶

#!/usr/bin/env python
# encoding: utf-8
r"""
2D shallow water: radial dam break
==================================

Solve the 2D shallow water equations:

.. math::
h_t + (hu)_x + (hv)_y = 0 \\
(hu)_t + (hu^2 + \frac{1}{2}gh^2)_x + (huv)_y = 0 \\
(hv)_t + (huv)_x + (hv^2 + \frac{1}{2}gh^2)_y = 0.

The initial condition is a circular area with high depth surrounded by lower-depth water.
The top and right boundary conditions reflect, while the bottom and left boundaries
are outflow.
"""

from __future__ import absolute_import
import numpy as np
from clawpack import riemann
from clawpack.riemann.shallow_roe_with_efix_2D_constants import depth, x_momentum, y_momentum, num_eqn

x0=0.
y0=0.
X, Y = state.p_centers
r = np.sqrt((X-x0)**2 + (Y-y0)**2)

state.q[x_momentum,:,:] = 0.
state.q[y_momentum,:,:] = 0.

def setup(use_petsc=False,outdir='./_output',solver_type='classic'):
if use_petsc:
import clawpack.petclaw as pyclaw
else:
from clawpack import pyclaw

if solver_type == 'classic':
solver = pyclaw.ClawSolver2D(riemann.shallow_roe_with_efix_2D)
solver.limiters = pyclaw.limiters.tvd.MC
solver.dimensional_split=1
elif solver_type == 'sharpclaw':
solver = pyclaw.SharpClawSolver2D(riemann.shallow_roe_with_efix_2D)

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

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

state = pyclaw.State(domain,num_eqn)

# Gravitational constant
state.problem_data['grav'] = 1.0

qinit(state)

claw = pyclaw.Controller()
claw.tfinal = 2.5
claw.solution = pyclaw.Solution(state,domain)
claw.solver = solver
claw.outdir = outdir
claw.num_output_times = 10
claw.setplot = setplot
claw.keep_copy = True

return claw

#--------------------------
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 depth
plotfigure = plotdata.new_plotfigure(name='Water height', figno=0)

# Set up for axes in this figure:
plotaxes = plotfigure.new_plotaxes()
plotaxes.xlimits = [-2.5, 2.5]
plotaxes.ylimits = [-2.5, 2.5]
plotaxes.title = 'Water height'
plotaxes.scaled = True

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

# Scatter plot of depth
plotfigure = plotdata.new_plotfigure(name='Scatter plot of h', figno=1)

# Set up for axes in this figure:
plotaxes = plotfigure.new_plotaxes()
plotaxes.xlimits = [0., 2.5]
plotaxes.ylimits = [0., 2.1]
plotaxes.title = 'Scatter plot of h'

# Set up for item on these axes:
plotitem = plotaxes.new_plotitem(plot_type='1d_from_2d_data')
plotitem.plot_var = depth
from numpy import sqrt
x = current_data.x
y = current_data.y
r = sqrt(x**2 + y**2)
q = current_data.q[depth,:,:]
return r,q
plotitem.plotstyle = 'o'

# Figure for x-momentum
plotfigure = plotdata.new_plotfigure(name='Momentum in x direction', figno=2)

# Set up for axes in this figure:
plotaxes = plotfigure.new_plotaxes()
plotaxes.xlimits = [-2.5, 2.5]
plotaxes.ylimits = [-2.5, 2.5]
plotaxes.title = 'Momentum in x direction'
plotaxes.scaled = True

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

# Figure for y-momentum
plotfigure = plotdata.new_plotfigure(name='Momentum in y direction', figno=3)

# Set up for axes in this figure:
plotaxes = plotfigure.new_plotaxes()
plotaxes.xlimits = [-2.5, 2.5]
plotaxes.ylimits = [-2.5, 2.5]
plotaxes.title = 'Momentum in y direction'
plotaxes.scaled = True

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