Radial ocean with island test problem

tsunami/radial-ocean/awr2011_5.x/Makefile

@@ -0,0 +1,69 @@
+
+# Makefile for Clawpack code in this directory.
+# This version only sets the local files and frequently changed
+# options, and then includes the standard makefile pointed to by CLAWMAKE.
+CLAWMAKE = $(CLAW)/clawutil/src/Makefile.common
+
+# See the above file for details and a list of make options, or type
+#   make .help
+# at the unix prompt.
+
+
+# Adjust these variables if desired:
+# ----------------------------------
+
+CLAW_PKG = geoclaw                  # Clawpack package to use
+EXE = xgeoclaw                 # Executable to create
+SETRUN_FILE = setrun.py        # File containing function to make data
+OUTDIR = _output               # Directory for output
+SETPLOT_FILE = setplot.py      # File containing function to set plots
+PLOTDIR = _plots               # Directory for plots
+
+# Environment variable FC should be set to fortran compiler, e.g. gfortran
+
+# Compiler flags can be specified here or set as an environment variable
+FFLAGS ?= 
+
+# ---------------------------------
+# package sources for this program:
+# ---------------------------------
+
+GEOLIB = $(CLAW)/geoclaw/src/2d/shallow
+include $(GEOLIB)/Makefile.geoclaw
+
+# ---------------------------------------
+# package sources specifically to exclude
+# (i.e. if a custom replacement source 
+#  under a different name is provided)
+# ---------------------------------------
+
+EXCLUDE_MODULES = \
+
+EXCLUDE_SOURCES = \
+
+# ----------------------------------------
+# List of custom sources for this program:
+# ----------------------------------------
+
+
+MODULES = \
+
+SOURCES = \
+  $(CLAW)/riemann/src/rpn2_geoclaw.f \
+  $(CLAW)/riemann/src/rpt2_geoclaw.f \
+  $(CLAW)/riemann/src/geoclaw_riemann_utils.f \
+
+#-------------------------------------------------------------------
+# Include Makefile containing standard definitions and make options:
+include $(CLAWMAKE)
+
+# Construct the topography data
+.PHONY: topo all
+topo:
+	python maketopo.py
+
+all: 
+	$(MAKE) topo
+	$(MAKE) .plots
+	$(MAKE) .htmls
+

tsunami/radial-ocean/awr2011_5.x/Makefile_220

@@ -0,0 +1,69 @@
+
+# Makefile for Clawpack code in this directory.
+# This version only sets the local files and frequently changed
+# options, and then includes the standard makefile pointed to by CLAWMAKE.
+CLAWMAKE = $(CLAW)/clawutil/src/Makefile.common
+
+# See the above file for details and a list of make options, or type
+#   make .help
+# at the unix prompt.
+
+
+# Adjust these variables if desired:
+# ----------------------------------
+
+CLAW_PKG = geoclaw                  # Clawpack package to use
+EXE = xgeoclaw                 # Executable to create
+SETRUN_FILE = setrun_220.py    # File containing function to make data
+OUTDIR = _output_220           # Directory for output
+SETPLOT_FILE = setplot.py      # File containing function to set plots
+PLOTDIR = _plots_220           # Directory for plots
+
+# Environment variable FC should be set to fortran compiler, e.g. gfortran
+
+# Compiler flags can be specified here or set as an environment variable
+FFLAGS ?= 
+
+# ---------------------------------
+# package sources for this program:
+# ---------------------------------
+
+GEOLIB = $(CLAW)/geoclaw/src/2d/shallow
+include $(GEOLIB)/Makefile.geoclaw
+
+# ---------------------------------------
+# package sources specifically to exclude
+# (i.e. if a custom replacement source 
+#  under a different name is provided)
+# ---------------------------------------
+
+EXCLUDE_MODULES = \
+
+EXCLUDE_SOURCES = \
+
+# ----------------------------------------
+# List of custom sources for this program:
+# ----------------------------------------
+
+
+MODULES = \
+
+SOURCES = \
+  $(CLAW)/riemann/src/rpn2_geoclaw.f \
+  $(CLAW)/riemann/src/rpt2_geoclaw.f \
+  $(CLAW)/riemann/src/geoclaw_riemann_utils.f \
+
+#-------------------------------------------------------------------
+# Include Makefile containing standard definitions and make options:
+include $(CLAWMAKE)
+
+# Construct the topography data
+.PHONY: topo all
+topo:
+	python maketopo.py
+
+all: 
+	$(MAKE) topo
+	$(MAKE) .plots
+	$(MAKE) .htmls
+

tsunami/radial-ocean/awr2011_5.x/Makefile_260

@@ -0,0 +1,69 @@
+
+# Makefile for Clawpack code in this directory.
+# This version only sets the local files and frequently changed
+# options, and then includes the standard makefile pointed to by CLAWMAKE.
+CLAWMAKE = $(CLAW)/clawutil/src/Makefile.common
+
+# See the above file for details and a list of make options, or type
+#   make .help
+# at the unix prompt.
+
+
+# Adjust these variables if desired:
+# ----------------------------------
+
+CLAW_PKG = geoclaw                  # Clawpack package to use
+EXE = xgeoclaw                 # Executable to create
+SETRUN_FILE = setrun_260.py    # File containing function to make data
+OUTDIR = _output_260           # Directory for output
+SETPLOT_FILE = setplot.py      # File containing function to set plots
+PLOTDIR = _plots_260           # Directory for plots
+
+# Environment variable FC should be set to fortran compiler, e.g. gfortran
+
+# Compiler flags can be specified here or set as an environment variable
+FFLAGS ?= 
+
+# ---------------------------------
+# package sources for this program:
+# ---------------------------------
+
+GEOLIB = $(CLAW)/geoclaw/src/2d/shallow
+include $(GEOLIB)/Makefile.geoclaw
+
+# ---------------------------------------
+# package sources specifically to exclude
+# (i.e. if a custom replacement source 
+#  under a different name is provided)
+# ---------------------------------------
+
+EXCLUDE_MODULES = \
+
+EXCLUDE_SOURCES = \
+
+# ----------------------------------------
+# List of custom sources for this program:
+# ----------------------------------------
+
+
+MODULES = \
+
+SOURCES = \
+  $(CLAW)/riemann/src/rpn2_geoclaw.f \
+  $(CLAW)/riemann/src/rpt2_geoclaw.f \
+  $(CLAW)/riemann/src/geoclaw_riemann_utils.f \
+
+#-------------------------------------------------------------------
+# Include Makefile containing standard definitions and make options:
+include $(CLAWMAKE)
+
+# Construct the topography data
+.PHONY: topo all
+topo:
+	python maketopo.py
+
+all: 
+	$(MAKE) topo
+	$(MAKE) .plots
+	$(MAKE) .htmls
+

tsunami/radial-ocean/awr2011_5.x/README.rst

@@ -0,0 +1,34 @@
+
+.. _apps_tsunami_radial-ocean_awr2011_5.x:
+
+Radial ocean test case from AWR 2011 paper
+==========================================
+
+
+These codes were developed to accompany the paper:
+
+
+The GeoClaw software for depth-averaged flows with adaptive refinement, 
+by M.J. Berger, D.L. George, R.J. LeVeque, and K.M.  Mandli.  
+Advances in Water Resources 34 (2011), pp. 1195-1206.
+`link <http://www.amath.washington/edu/~rjl/pubs/awr10>`_
+
+Radially symmetric ocean on a lat-long grid, so it doesn't look circular in the
+computational coordinates.
+
+The profile of the ocean is set in function topo of `maketopo.py`.
+You must first set the desired value `theta_island` in `setrun.py` and do::
+
+    make data
+    make topo
+    make output
+
+To run both test problems from the paper and produce some plots::
+
+    source make_all.sh
+
+To make comparison plots of gauges from the two runs:
+
+    python compare_gauges.py
+
+This creates gauges1-2.png and gauges3-4.png

tsunami/radial-ocean/awr2011_5.x/compare_gauges.py

@@ -0,0 +1,64 @@
+
+import pylab
+from clawpack.visclaw.data import ClawPlotData
+
+outdir1 = '_output_260'
+outdir2 = '_output_220'
+gaugenos = [1,2]
+
+plotdata = ClawPlotData()
+
+pylab.figure(301)
+pylab.clf()
+
+plotdata.outdir = outdir1
+for gaugeno in gaugenos:
+    gs = plotdata.getgauge(gaugeno)
+    pylab.plot(gs.t, gs.q[3,:], 'b',linewidth=2)
+
+
+plotdata.outdir = outdir2
+for gaugeno in gaugenos:
+    gs = plotdata.getgauge(gaugeno)
+    pylab.plot(gs.t, gs.q[3,:], 'r--',linewidth=2)
+
+pylab.xlim([7500,14000])
+pylab.ylim([-2.5,4])
+pylab.xticks(fontsize=15)
+pylab.yticks(fontsize=15)
+
+pylab.annotate('Gauge 1',[8200,1.9],[7700,2.5],arrowprops={'width':1,'color':'k'})
+pylab.annotate('Gauge 2',[9300,2.7] ,[9700,3.1],arrowprops={'width':1,'color':'k'})
+
+pylab.savefig('gauges1-2.png')
+
+print("Created gauges1-2.png")
+
+#===============================
+
+gaugenos = [3,4]
+pylab.figure(302)
+pylab.clf()
+
+plotdata.outdir = outdir1
+for gaugeno in gaugenos:
+    gs = plotdata.getgauge(gaugeno)
+    pylab.plot(gs.t, gs.q[3,:], 'b',linewidth=2)
+
+
+plotdata.outdir = outdir2
+for gaugeno in gaugenos:
+    gs = plotdata.getgauge(gaugeno)
+    pylab.plot(gs.t, gs.q[3,:], 'r--',linewidth=2)
+
+pylab.xlim([7500,14000])
+pylab.ylim([-2.5,4])
+pylab.xticks(fontsize=15)
+pylab.yticks(fontsize=15)
+
+pylab.annotate('Gauge 3',[12500,2.6],[13000,3.1],arrowprops={'width':1,'color':'k'})
+pylab.annotate('Gauge 4',[11550,1.0] ,[11200,2.0],arrowprops={'width':1,'color':'k'})
+
+pylab.savefig('gauges3-4.png')
+
+print("Created gauges3-4.png")

tsunami/radial-ocean/awr2011_5.x/interp.py

@@ -0,0 +1,83 @@
+"""
+Tools for interpolating data.
+
+Includes:
+    pwcubic:        piecewise cubic interpolation.
+    pwlinear:       piecewise linear interpolation.
+"""
+
+
+def pwcubic(xi, zl, zr, slopel, sloper, x):
+    """
+    Construct a piecewise cubic function and evaluate at the points in x.
+    The interpolation data consists of points xi, and values
+        zl[i] = value of z in limit as x approaches xi[i] from the left,
+        zr[i] = value of z in limit as x approaches xi[i] from the right,
+        slopel[i] = value of z'(x) in limit as x approaches xi[i] from the left,
+        sloper[i] = value of z'(x) in limit as x approaches xi[i] from the right.
+    To the left of xi[0] the linear function based on 
+       zl[0] and slopel[0] is used.
+    To the right of xi[-1] the linear function based on 
+       zr[-1] and sloper[-1] is used.
+
+    In the interior intervals a cubic us used that interpolates 
+    the 4 values at the ends of the interval.  
+
+    Note that the function will be linear in the i'th interval if
+        sloper[i]   = (zl[i+1] - zr[i]) / (xi[i+1] - xi[i])
+        slopel[i+1] = (zl[i+1] - zr[i]) / (xi[i+1] - xi[i])
+
+    A piecewise linear interpolant can be created by setting this everywhere:
+        slopel[1:]  = (zl[1:]-zr[:-1]) / (xi[1:] - xi[:-1])
+        sloper[:-1] = (zl[1:]-zr[:-1]) / (xi[1:] - xi[:-1])
+    This can be done automatically with the pwlinear function defined below.
+
+    Note that the pwcubic function is continuous if zl == zr and is C^1 if
+    in addition slopel == sloper.
+    """
+
+    from numpy import where, zeros
+
+    dx = xi[1:] - xi[:-1]
+
+    s = (zl[1:] - zr[:-1]) / dx
+    c2 = (s - sloper[:-1]) / dx
+    c3 = (slopel[1:] - 2.*s + sloper[:-1]) / (dx**2)
+
+    # set to linear function for x<xi[0] or x>= xi[-1]
+    z = where(x < xi[0],  zl[0] + (x-xi[0]) * slopel[0], 
+                          zr[-1] + (x-xi[-1]) * sloper[-1])
+
+    # replace by appropriate cubic in intervals:
+    for i in range(len(xi)-1):
+        cubic = zr[i] + sloper[i]*(x - xi[i]) \
+                    + (x-xi[i])**2 * (c2[i] + c3[i]*(x-xi[i+1]))
+        z = where((x >= xi[i]) & (x < xi[i+1]), cubic, z)
+
+    return z
+
+
+def pwlinear(xi, zl, zr, x, extrap=0):
+    """
+    Construct a piecewise linear function and evaluate at the points in x.
+    The interpolation data consists of points xi, and values
+        zl[i] = value of z in limit as x approaches xi[i] from the left,
+        zr[i] = value of z in limit as x approaches xi[i] from the right,
+    Extrapolate outside the interval by:
+        if extrap==0:  slope = 0 outside of intervals
+        if extrap==1:  slope taken from neighboring interval.
+    Note that the pwlinear function is continuous if zl == zr.
+    """
+
+    from numpy import zeros
+    slopel = zeros(len(zl))
+    sloper = zeros(len(zr))
+    dx = xi[1:] - xi[:-1]
+    slopel[1:]  = (zl[1:]-zr[:-1]) / dx
+    sloper[:-1] = slopel[1:]
+    if extrap==1:
+        slopel[0] = sloper[0]
+        sloper[-1] = slopel[-1]
+    z = pwcubic(xi, zl, zr, slopel, sloper, x)
+    return z
+