Merge pull request #466 from mandli/add-sill-example

Add shallow water sill example

examples/__init__.py

@@ -12,6 +12,7 @@
 from kpp import kpp
 from psystem_2d import psystem_2d
 from shallow_1d import dam_break
+from shallow_1d import sill
 from shallow_2d import radial_dam_break
 from shallow_sphere import Rossby_wave
 from stegoton_1d import stegoton

examples/shallow_1d/sill.py

@@ -0,0 +1,128 @@
+#!/usr/bin/env python
+# encoding: utf-8
+
+r"""
+Shallow water flow
+==================
+
+Solve the one-dimensional shallow water equations including bathymetry:
+
+.. math::
+    h_t + (hu)_x & = 0 \\
+    (hu)_t + (hu^2 + \frac{1}{2}gh^2)_x & = -g h b_x.
+
+Here h is the depth, u is the velocity, g is the gravitational constant, and b
+the bathymetry.  
+"""
+
+import numpy
+from clawpack import riemann
+
+def setup(use_petsc=False, outdir='./_output', solver_type='classic'):
+
+    if use_petsc:
+        import clawpack.petclaw as pyclaw
+    else:
+        from clawpack import pyclaw
+
+    solver = pyclaw.ClawSolver1D(riemann.shallow_1D_py.shallow_fwave_1d)
+    solver.limiters = pyclaw.limiters.tvd.vanleer
+    solver.kernel_language = "Python"
+    solver.fwave = True
+    solver.num_waves = 2
+    solver.num_eqn = 2
+    solver.bc_lower[0] = pyclaw.BC.extrap
+    solver.bc_upper[0] = pyclaw.BC.extrap
+    solver.aux_bc_lower[0] = pyclaw.BC.extrap
+    solver.aux_bc_upper[0] = pyclaw.BC.extrap
+
+    xlower = -1.0
+    xupper = 1.0
+    x = pyclaw.Dimension('x', xlower, xupper, 500)
+    domain = pyclaw.Domain(x)
+    state = pyclaw.State(domain, 2, 1)
+
+    # Gravitational constant
+    state.problem_data['grav'] = 9.8
+    state.problem_data['sea_level'] = 0.0
+
+    xc = state.grid.x.centers
+    state.aux[0, :] = 0.8 * numpy.exp(-xc**2 / 0.2**2) - 1.0
+    state.q[0, :] = 0.1 * numpy.exp(-(xc + 0.4)**2 / 0.2**2) - state.aux[0, :]
+    state.q[1, :] = 0.0
+
+    claw = pyclaw.Controller()
+    claw.keep_copy = True
+    claw.tfinal = 1.0
+    claw.solution = pyclaw.Solution(state, domain)
+    claw.solver = solver
+    claw.outdir = outdir
+    claw.setplot = setplot
+    claw.write_aux_init = 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.
+    """ 
+    plotdata.clearfigures()  # clear any old figures,axes,items data
+
+    # Plot variables
+    def bathy(current_data):
+        return current_data.aux[0, :]
+
+    def eta(current_data):
+        return current_data.q[0, :] + bathy(current_data)
+
+    def velocity(current_data):
+        return current_data.q[1, :] / current_data.q[0, :]
+
+    rgb_converter = lambda triple: [float(rgb) / 255.0 for rgb in triple]
+
+    # Figure for depth
+    plotfigure = plotdata.new_plotfigure(name='Depth', figno=0)
+
+    # Axes for water depth
+    plotaxes = plotfigure.new_plotaxes()
+    plotaxes.xlimits = [-1.0, 1.0]
+    plotaxes.ylimits = [-1.1, 0.2]
+    plotaxes.title = 'Water Depth'
+    plotaxes.axescmd = 'subplot(211)'
+
+    plotitem = plotaxes.new_plotitem(plot_type='1d_fill_between')
+    plotitem.plot_var = eta
+    plotitem.plot_var2 = bathy
+    plotitem.color = rgb_converter((67,183,219))
+
+    plotitem = plotaxes.new_plotitem(plot_type='1d_plot')
+    plotitem.plot_var = bathy
+    plotitem.color = 'k'
+
+    plotitem = plotaxes.new_plotitem(plot_type='1d_plot')
+    plotitem.plot_var = eta
+    plotitem.color = 'k'
+
+    # Axes for velocity
+    plotaxes = plotfigure.new_plotaxes()
+    plotaxes.axescmd = 'subplot(212)'
+    plotaxes.xlimits = [-1.0, 1.0]
+    plotaxes.ylimits = [-0.5, 0.5]
+    plotaxes.title = 'Velocity'
+
+    plotitem = plotaxes.new_plotitem(plot_type='1d_plot')
+    plotitem.plot_var = velocity
+    plotitem.color = 'b'
+    plotitem.kwargs = {'linewidth':3}
+
+    return plotdata
+
+
+if __name__=="__main__":
+    from clawpack.pyclaw.util import run_app_from_main
+    output = run_app_from_main(setup,setplot)