Merge pull request #35 from harpolea/swe

add a shallow water solver

Author: zingale

README.md

@@ -11,7 +11,7 @@ finite-volume framework.  The code is mainly written in python and is
 designed with simplicity in mind.  The algorithms are written to
 encourage experimentation and allow for self-learning of these code
 methods.
- 
+
 The latest version of pyro is always available at:
 
 https://github.com/zingale/pyro2
@@ -154,6 +154,8 @@ pyro provides the following solvers (all in 2-d):
     variable-coefficient Poisson equation (which inherits from the
     constant-coefficient solver).
 
+  - `swe`: a solver for the shallow water equations.
+
 
 ## Working with data
 
@@ -175,6 +177,12 @@ with their data.
 
   - `analysis/`
 
+      * `dam_compare.py`: this takes an output file from the
+        shallow water dam break problem and plots a slice through the domain
+        together with the analytic solution (calculated in the script).
+
+         usage: `./dam_compare.py file`
+
       * `gauss_diffusion_compare.py`: this is for the diffusion solver's
         Gaussian diffusion problem. It takes a sequence of output
         files as arguments, computes the angle-average, and the plots

analysis/dam_compare.py

@@ -0,0 +1,183 @@
+#!/usr/bin/env python3
+
+from __future__ import print_function
+
+import numpy as np
+from scipy.optimize import brentq
+import sys
+import os
+import matplotlib.pyplot as plt
+from util import msg, runparams, io
+
+usage = """
+      compare the output for a dam problem with the exact solution contained
+      in dam-exact.out.
+
+      usage: ./dam_compare.py file
+"""
+
+
+def abort(string):
+    print(string)
+    sys.exit(2)
+
+
+if not len(sys.argv) == 2:
+    print(usage)
+    sys.exit(2)
+
+try:
+    file1 = sys.argv[1]
+except IndexError:
+    print(usage)
+    sys.exit(2)
+
+sim = io.read(file1)
+myd = sim.cc_data
+myg = myd.grid
+
+# time of file
+t = myd.t
+if myg.nx > myg.ny:
+    # x-problem
+    xmin = myg.xmin
+    xmax = myg.xmax
+    param_file = "inputs.dam.x"
+else:
+    # y-problem
+    xmin = myg.ymin
+    xmax = myg.ymax
+    param_file = "inputs.dam.y"
+
+
+height = myd.get_var("height")
+xmom = myd.get_var("x-momentum")
+ymom = myd.get_var("y-momentum")
+
+# get the 1-d profile from the simulation data -- assume that whichever
+# coordinate is the longer one is the direction of the problem
+
+# parameter defaults
+rp = runparams.RuntimeParameters()
+rp.load_params("../_defaults")
+rp.load_params("../swe/_defaults")
+rp.load_params("../swe/problems/_dam.defaults")
+
+# now read in the inputs file
+if not os.path.isfile(param_file):
+    # check if the param file lives in the solver's problems directory
+    param_file = "../swe/problems/" + param_file
+    if not os.path.isfile(param_file):
+        msg.fail("ERROR: inputs file does not exist")
+
+rp.load_params(param_file, no_new=1)
+
+if myg.nx > myg.ny:
+    # x-problem
+    x = myg.x[myg.ilo:myg.ihi+1]
+    jj = myg.ny//2
+
+    h = height[myg.ilo:myg.ihi+1, jj]
+
+    u = xmom[myg.ilo:myg.ihi+1, jj]/h
+    ut = ymom[myg.ilo:myg.ihi+1, jj]/h
+
+else:
+
+    # y-problem
+    x = myg.y[myg.jlo:myg.jhi+1]
+    ii = myg.nx//2
+
+    h = height[ii, myg.jlo:myg.jhi+1]
+
+    u = ymom[ii, myg.jlo:myg.jhi+1]/h
+    ut = xmom[ii, myg.jlo:myg.jhi+1]/h
+
+print(myg)
+
+x_exact = x
+h_exact = np.zeros_like(x)
+u_exact = np.zeros_like(x)
+
+# find h0, h1
+h1 = rp.get_param("dam.h_left")
+h0 = rp.get_param("dam.h_right")
+
+
+def find_h2(h2):
+    return (h2/h1)**3 - 9*(h2/h1)**2*(h0/h1) + \
+        16*(h2/h1)**1.5*(h0/h1) - (h2/h1)*(h0/h1)*(h0/h1+8) + \
+        (h0/h1)**3
+
+
+h2 = brentq(find_h2, min(h0, h1), max(h0, h1))
+
+# calculate sound speeds
+g = rp.get_param("swe.grav")
+
+c0 = np.sqrt(g*h0)
+c1 = np.sqrt(g*h1)
+c2 = np.sqrt(g*h2)
+
+u2 = 2 * (c1 - c2)
+
+# shock speed
+xi = c0 * np.sqrt(1/8 * ((2*(c2/c0)**2 + 1)**2 - 1))
+
+xctr = 0.5*(xmin + xmax)
+
+# h0
+idx = x >= xctr + xi*t
+h_exact[idx] = h0
+u_exact[idx] = 0
+
+# h1
+idx = x <= xctr - c1*t
+h_exact[idx] = h1
+u_exact[idx] = 0
+
+# h2
+idx = ((x >= xctr + (u2-c2)*t) & (x < xctr + xi*t))
+h_exact[idx] = h2
+u_exact[idx] = u2
+
+# h3
+idx = ((x >= xctr - c1*t) & (x < xctr + (u2-c2)*t))
+c3 = 1/3 * (2*c1 - (x-xctr)/t)
+h_exact[idx] = c3[idx]**2 / g
+u_exact[idx] = 2 * (c1-c3[idx])
+
+# plot
+fig, axes = plt.subplots(nrows=2, ncols=1, num=1)
+
+plt.rc("font", size=10)
+
+ax = axes.flat[0]
+
+ax.plot(x_exact, h_exact, label='Exact')
+ax.scatter(x, h, marker="x", s=7, color="r", label='Pyro')
+
+ax.set_ylabel(r"$h$")
+ax.set_xlim(0, 1.0)
+ax.set_ylim(0, 1.1)
+
+ax = axes.flat[1]
+
+ax.plot(x_exact, u_exact)
+ax.scatter(x, u, marker="x", s=7, color="r")
+
+ax.set_ylabel(r"$u$")
+ax.set_xlim(0, 1.0)
+
+if (myg.nx > myg.ny):
+    ax.set_xlabel(r"x")
+else:
+    ax.set_xlabel(r"y")
+
+lgd = axes.flat[0].legend()
+
+plt.subplots_adjust(hspace=0.25)
+
+fig.set_size_inches(8.0, 8.0)
+
+plt.savefig("dam_compare.png", bbox_inches="tight")

docs/source/analysis.rst

@@ -27,6 +27,12 @@ with their data.
     the finer data and then computes the norm of the difference.  This
     is used to test the convergence of solvers.
 
+  * ``dam_compare.py``: this takes an output file from the
+    shallow water dam break problem and plots a slice through the domain
+    together with the analytic solution (calculated in the script).
+
+    usage: ``./dam_compare.py file``
+
   * ``gauss_diffusion_compare.py``: this is for the diffusion solver's
     Gaussian diffusion problem. It takes a sequence of output files as
     arguments, computes the angle-average, and the plots the resulting
@@ -65,6 +71,3 @@ with their data.
     the analytic solution (read in from ``sod-exact.out``).
 
     usage: ``./sod_compare.py file``
-
-
-

docs/source/design.rst

@@ -85,6 +85,11 @@ The overall structure is:
   * ``problems/``: The problem setups for when the multigrid solver is used in a stand-alone fashion.
   * ``tests/``: Reference multigrid solver solutions (from when the multigrid solver is used stand-alone) for regression testing.
 
+* ``swe/``: The shallow water solver.
+
+  * ``problems/``: The problem setups for the shallow water solver.
+  * ``tests/``: Reference shallow water output for regression testing.
+
 * ``util/``: Various service modules used by the pyro routines,
   including runtime parameters, I/O, profiling, and pretty output
   modes.

docs/source/index.rst

@@ -37,6 +37,7 @@ pyro: a python hydro code
    diffusion_basics
    incompressible_basics
    lowmach_basics
+   swe_basics
 
 .. toctree::
    :maxdepth: 1

docs/source/modules.rst

@@ -23,5 +23,6 @@ pyro2
    plot
    pyro
    simulation_null
+   swe
    test
    util

docs/source/swe.problems.rst

@@ -0,0 +1,58 @@
+swe\.problems package
+==============================
+
+.. automodule:: swe.problems
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+Submodules
+----------
+
+swe\.problems\.acoustic\_pulse module
+----------------------------------------------
+
+.. automodule:: swe.problems.acoustic_pulse
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+swe\.problems\.advect module
+-------------------------------------
+
+.. automodule:: swe.problems.advect
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+swe\.problems\.dam module
+----------------------------------
+
+.. automodule:: swe.problems.dam
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+swe\.problems\.kh module
+---------------------------------
+
+.. automodule:: swe.problems.kh
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+swe\.problems\.quad module
+-----------------------------------
+
+.. automodule:: swe.problems.quad
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+swe\.problems\.test module
+-----------------------------------
+
+.. automodule:: swe.problems.test
+    :members:
+    :undoc-members:
+    :show-inheritance:

docs/source/swe.rst

@@ -0,0 +1,41 @@
+swe package
+====================
+
+.. automodule:: swe
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+Subpackages
+-----------
+
+.. toctree::
+
+    swe.problems
+
+Submodules
+----------
+
+swe\.derives module
+----------------------------
+
+.. automodule:: swe.derives
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+swe\.simulation module
+-------------------------------
+
+.. automodule:: swe.simulation
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+swe\.unsplit\_fluxes module
+------------------------------------
+
+.. automodule:: swe.unsplit_fluxes
+    :members:
+    :undoc-members:
+    :show-inheritance: