Merge pull request #7 from geoflows/radial_slide

Initial addition of radial slide example

examples/1d_slide/Makefile

@@ -0,0 +1,78 @@
+# 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 = dclaw                  # Clawpack package to use
+
+
+EXE = $(PWD)/xgeoclaw
+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
+FC = gfortran
+FFLAGS ?= -O  -gno-strict-dwarf -fbounds-check  -fopenmp \
+      -std=legacy -ffpe-trap='invalid,overflow,zero'
+
+# ---------------------------------
+# package sources for this program:
+# ---------------------------------
+
+# DIGLIB contains library code for D-Claw solvers,
+# AMRLIB and GEOLIB are standard libraries, defined in case you need to
+# exclude some modele or source file from one of these.
+DIGLIB = $(CLAW)/dclaw/src/2d/dig
+AMRLIB = $(CLAW)/amrclaw/src/2d
+GEOLIB = $(CLAW)/geoclaw/src/2d/shallow
+
+# See this Makefile for the list of library routines used:
+include $(DIGLIB)/Makefile.dclaw
+
+# ---------------------------------------
+# 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 = \
+
+#-------------------------------------------------------------------
+# Include Makefile containing standard definitions and make options:
+include $(CLAWMAKE)
+
+# Construct the topography data
+.PHONY: topo all
+topo:
+	$(CLAW_PYTHON) maketopo.py
+
+all: 
+	$(MAKE) topo
+	$(MAKE) .plots
+	$(MAKE) .htmls
+

examples/1d_slide/maketopo.py

@@ -0,0 +1,168 @@
+
+from pylab import *
+from clawpack.geoclaw import topotools
+from clawpack.visclaw import colormaps
+from clawpack.geoclaw.data import Rearth  # radius of earth
+from clawpack.clawutil.data import ClawData
+from scipy.interpolate import interp1d
+
+
+# # 1d symmetric landslide (and tsunami)
+
+def dist(x,y):
+    #d = np.sqrt(x**2 + y**2)  # radial
+    d = abs(x)  # 1d in x
+    #d = abs(y)  # 1d in y
+    return d
+
+m0 = 0.63  # mass fraction in landslide
+
+r0 = 0.     # left boundary (center of radially symmetric)
+r1 = 1100.  # top of landslide
+r2 = 1300.
+r3 = 1600.  # toe of landslide
+r4 = 2000.  # start of flat
+r5 = 5000.*sqrt(2)   # right boundary for radial coordinate
+
+B4 = B5 = 0 # flat region
+Bslope = 25 # degrees
+B3 = (r4 - r3) * sin(Bslope * pi/180)
+B1 = (r4 - r1) * sin(Bslope * pi/180)
+B2 = (r4 - r2) * sin(Bslope * pi/180)
+B0 = B1
+etaslope = 40 # degrees
+eta2 = B1
+
+rr = array([r0,r1,r2,r3,r4,r5])
+Br = array([B0,B1,B2,B3,B4,B5])
+etar = Br.copy()
+#etar[2] = eta2
+etar[:3] = 420.
+
+def make_plots():
+    figure(figsize=(10,3))
+    fill_between(rr,Br,etar,color='brown')
+    plot(rr,Br,'g')
+    grid(True)
+    #axis('equal')
+    plot([r1,r1],[0,500],'k--')
+    text(r1,-20,'r1\n%.0f' % r1,va='top',ha='center')
+    plot([r2,r2],[0,500],'k--')
+    text(r2,-20,'r2\n%.0f' % r2,va='top',ha='center')
+    plot([r3,r3],[0,500],'k--')
+    text(r3,-20,'r3\n%.0f' % r3,va='top',ha='center')
+    plot([r4,r4],[0,500],'k--')
+    text(r4,-20,'r4\n%.0f' % r4,va='top',ha='center')
+    axis([0,3000,-100,500])
+    fname = '1d_topo.png'
+    savefig(fname)
+    print('Created ',fname)
+
+    basal = topotools.Topography('basal_topo_1dx.tt3',3)
+    basal.plot()
+    title('Basal topo')
+    fname = 'basal_topo.png'
+    savefig(fname)
+    print('Created ',fname)
+
+    eta = topotools.Topography('surface_topo_1dx.tt3',3)
+    eta.plot()
+    title('Surface topo eta')
+    fname = 'surface_topo.png'
+    savefig(fname)
+    print('Created ',fname)
+
+    h = eta.Z - basal.Z
+    figure()
+    pcolormesh(eta.X,eta.Y,h,cmap=colormaps.white_red)
+    axis('equal')
+    colorbar()
+    title('Landslide depth')
+    fname = 'landslide_depth.png'
+    savefig(fname)
+    print('Created ',fname)
+
+
+landslide_depth = eta2 - B2
+print('Maximum landslide depth: %.2f m' % landslide_depth)
+
+#x1d, z1d = loadtxt('1d_radial/celledges.data',skiprows=1,unpack=True)
+B1d_func = interp1d(rr, Br, bounds_error=False, fill_value = Br[-1])
+
+def basal(x,y):
+    """
+    Cartesian: x,y in meters
+    """
+    import numpy as np
+    x0 = 0.
+    y0 = 0.
+    d = dist(x-x0, y-y0)
+    z = B1d_func(d)
+    return z
+
+eta1d_func = interp1d(rr, etar, bounds_error=False, fill_value = etar[-1])
+
+def eta(x,y):
+    """
+    Cartesian: x,y in meters
+    """
+    import numpy as np
+    x0 = 0.
+    y0 = 0.
+    d = dist(x-x0, y-y0)
+    z = eta1d_func(d)
+    return z
+
+
+m1d_func = interp1d(rr, etar, bounds_error=False, fill_value = 0.)
+
+def mfrac(x,y):
+    """
+    mass fraction
+    Cartesian: x,y in meters
+    """
+    import numpy as np
+    x0 = 0.
+    y0 = 0.
+    d = dist(x-x0, y-y0)
+    eta = eta1d_func(d)
+    B = B1d_func(d)
+    mfrac = where(eta-B > 0, m0, 0.)
+    return mfrac
+
+
+xylim2d = 4e3
+
+def maketopo():
+    """
+    Output topography file for the entire domain
+    """
+    nxpoints = 501
+    nypoints = 501
+    xlower= -xylim2d
+    xupper=  xylim2d
+    ylower= -xylim2d
+    yupper=  xylim2d
+    outfile= "basal_topo_1dx.tt3"
+    topotools.topo3writer(outfile,basal,xlower,xupper,ylower,yupper,nxpoints,nypoints)
+    outfile= "mass_frac_1dx.tt3"
+    topotools.topo3writer(outfile,mfrac,xlower,xupper,ylower,yupper,nxpoints,nypoints)
+
+def make_surface():
+    """
+    Output surface topography file for the entire domain
+    (Could be for smaller region)
+    """
+    nxpoints = 501
+    nypoints = 501
+    xlower= -xylim2d
+    xupper=  xylim2d
+    ylower= -xylim2d
+    yupper=  xylim2d
+    outfile= "surface_topo_1dx.tt3"
+    topotools.topo3writer(outfile,eta,xlower,xupper,ylower,yupper,nxpoints,nypoints)
+
+if __name__=='__main__':
+    maketopo()
+    make_surface()
+    make_plots()