diff --git a/tsunami/radial-ocean/awr2011_5.x/Makefile b/tsunami/radial-ocean/awr2011_5.x/Makefile new file mode 100644 index 00000000..f2ca5bc6 --- /dev/null +++ b/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 + diff --git a/tsunami/radial-ocean/awr2011_5.x/Makefile_220 b/tsunami/radial-ocean/awr2011_5.x/Makefile_220 new file mode 100644 index 00000000..daf4de8a --- /dev/null +++ b/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 + diff --git a/tsunami/radial-ocean/awr2011_5.x/Makefile_260 b/tsunami/radial-ocean/awr2011_5.x/Makefile_260 new file mode 100644 index 00000000..481cca81 --- /dev/null +++ b/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 + diff --git a/tsunami/radial-ocean/awr2011_5.x/README.rst b/tsunami/radial-ocean/awr2011_5.x/README.rst new file mode 100644 index 00000000..ff6f6ba3 --- /dev/null +++ b/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 `_ + +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 diff --git a/tsunami/radial-ocean/awr2011_5.x/compare_gauges.py b/tsunami/radial-ocean/awr2011_5.x/compare_gauges.py new file mode 100644 index 00000000..f256e0a9 --- /dev/null +++ b/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") diff --git a/tsunami/radial-ocean/awr2011_5.x/interp.py b/tsunami/radial-ocean/awr2011_5.x/interp.py new file mode 100644 index 00000000..1f8c60a4 --- /dev/null +++ b/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[-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 + diff --git a/tsunami/radial-ocean/awr2011_5.x/make_all.sh b/tsunami/radial-ocean/awr2011_5.x/make_all.sh new file mode 100644 index 00000000..0d0c63b8 --- /dev/null +++ b/tsunami/radial-ocean/awr2011_5.x/make_all.sh @@ -0,0 +1,18 @@ + +make .exe + +make data -f Makefile_220 +make topo -f Makefile_220 +make output -f Makefile_220 +make plots -f Makefile_220 + +make data -f Makefile_260 +make topo -f Makefile_260 +make output -f Makefile_260 +make plots -f Makefile_260 + +python compare_gauges.py + +# To make other plots in paper: +python plot_ocean.py +python plot_xsec.py diff --git a/tsunami/radial-ocean/awr2011_5.x/maketopo.py b/tsunami/radial-ocean/awr2011_5.x/maketopo.py new file mode 100644 index 00000000..77c768ed --- /dev/null +++ b/tsunami/radial-ocean/awr2011_5.x/maketopo.py @@ -0,0 +1,128 @@ + +""" +Module to create topo and qinit data files for this example. +""" + +from __future__ import print_function +from clawpack.geoclaw import topotools +from clawpack.geoclaw.data import Rearth # radius of earth +from interp import pwcubic +from clawpack.clawutil.data import ClawData +from numpy import * + +from mapper import latlong, gcdist + +probdata = ClawData() +probdata.read('setprob.data', force=True) +theta_island = probdata.theta_island +print("theta_island = ",theta_island) + +(xisland,yisland) = latlong(1600.e3, theta_island, 40., Rearth) + +def maketopo(): + """ + Output topography file for the entire domain + and near-shore in one location. + """ + nxpoints=400 + nypoints=400 + xlower=-20.e0 + xupper= 20.e0 + ylower= 20.e0 + yupper= 60.e0 + outfile= "ocean.topotype2" + topotools.topo2writer(outfile,topo,xlower,xupper,ylower,yupper,nxpoints,nypoints) + + nxpoints=200 + nypoints=200 + xlower= xisland - 1. + xupper= xisland + 1. + ylower= yisland - 1. + yupper= yisland + 1. + outfile= "island.topotype2" + topotools.topo2writer(outfile,topo,xlower,xupper,ylower,yupper,nxpoints,nypoints) + +def makeqinit(): + """ + Create qinit data file + """ + nxpoints=100 + nypoints=100 + xlower=-20.e0 + xupper= 20.e0 + ylower= 20.e0 + yupper= 60.e0 + outfile= "hump.xyz" + topotools.topo1writer(outfile,qinit,xlower,xupper,ylower,yupper,nxpoints,nypoints) + +def shelf1(r): + """ + Ocean followed by continental slope, continental shelf, and beach. + The ocean is flat, the slope is cubic, and the shelf and beach are linear. + """ + + rad1 = 1500e3 # beginning of continental slope + rad2 = rad1 + 60e3 # end of slope, start of shelf + rad3 = rad2 + 80e3 # end of shelf, start of beach + rad4 = rad3 + 5e3 # radius of shoreline where z=0 + z1 = -4000. # depth from r=0 to rad1 (ocean) + z2 = -100. # depth at r=rad2 (start of shelf) + z3 = -100. # depth at r=rad3 (start of beach) + z4 = 0. # depth at shoreline + xi = array([0., rad1, rad2, rad3, rad4]) + zl = array([z1, z1, z2, z3, z4]) + zr = zl # continuous! + slope_of_shelf = (z3 - z2) / (rad3 - rad2) + slope_of_beach = (z4 - z3) / (rad4 - rad3) + print("Slope of shelf = ",slope_of_shelf) + print("Slope of beach = ",slope_of_beach) + slopel = array([0., 0., slope_of_shelf, slope_of_shelf, slope_of_beach]) + sloper = array([0., 0., slope_of_shelf, slope_of_beach, slope_of_beach]) + z = pwcubic(xi, zl, zr, slopel, sloper, r) + return z + +def island1(r): + """ + Island created using piecewise cubic with radius rad with zero slope at peak and + base. Raises topo by ztop at the peak. + """ + rad = 30e3 + ztop = 120. # 20m above sealevel on the 100m deep shelf. + xi = array([0., rad]) + zl = array([ztop, 0.]) + zr = zl + slopel = array([0., 0.]) + sloper = slopel + z = pwcubic(xi, zl, zr, slopel, sloper, r) + return z + +def topo(x,y): + """ + x = longitude, y = latitude in degrees + """ + import numpy as np + x0 = 0.0 + y0 = 40.0 + d = gcdist(x0,y0,x,y,Rearth) + z = shelf1(d) + d = gcdist(xisland,yisland,x,y,Rearth) + z = z + island1(d) + z = where(z>200, 200, z) # cut off topo at 200 m above sea level + return z + + +def qinit(x,y): + """ + Gaussian hump: + """ + from numpy import where + x0 = 0.0 + y0 = 40.0 + d = gcdist(x0,y0,x,y,Rearth) + ze = -d**2/2.e9 + z = where(ze>-100., 20.e0*exp(ze), 0.) + return z + +if __name__=='__main__': + maketopo() + makeqinit() diff --git a/tsunami/radial-ocean/awr2011_5.x/mapper.py b/tsunami/radial-ocean/awr2011_5.x/mapper.py new file mode 100644 index 00000000..8bd01831 --- /dev/null +++ b/tsunami/radial-ocean/awr2011_5.x/mapper.py @@ -0,0 +1,99 @@ + +from __future__ import print_function +from pylab import * + +def latlong(d,theta,phi,Rearth): + """ + Take a point at distance d from the North pole with longitude theta + and return longitude and latitude of point after rotating pole down + to latitude phi. + Applying this function to a set of points on a circle of radius d + will result in points that are equi-distant (distance d on the earth) + from the point at latitude phi, longitude 0. + + Used to construct a radially symmetric ocean on the earth and place + gauges, islands, etc. + """ + + # Convert phi to radians: + theta = theta * pi/180. + phi = phi * pi/180. + + alpha = pi/2. - d/Rearth # latitude of original point + # x,y,z coordinates of this point on the earth: + x1 = cos(alpha)*cos(theta) + y1 = cos(alpha)*sin(theta) + z1 = sin(alpha) + + # rotate so centered at latitude phi: + x2 = x1*sin(phi) + z1*cos(phi) + y2 = y1 + z2 = -x1*cos(phi) + z1*sin(phi) + + # compute longitude xhat and latitude yhat: + xhat = -arctan(y2/x2) + yhat = arcsin(z2) + + # convert to degrees: + xhat = xhat * 180./pi + yhat = yhat * 180./pi + + return xhat,yhat + +def plot_ocean_and_shelf(d1=1580e3, d2=1645e3, phi=40.): + + theta = linspace(0, 360., 200) + d = d1*ones(theta.shape) + xhat, yhat = latlong(d, theta, phi) + + clf() + plot(xhat,yhat,'b') + hold(True) + + d = d2*ones(theta.shape) + xhat, yhat = latlong(d, theta, phi) + + plot(xhat,yhat,'r') + legend(['Continental shelf', 'Shoreline'],loc='lower left') + + (xi1,yi1) = latlong(1600.e3,220.,40.) + plot([xi1],[yi1],'ko') + (xi2,yi2) = latlong(1600.e3,260.,40.) + plot([xi2],[yi2],'ko') + axis('scaled') + xlim([-20, 20]) + ylim([15, 60]) + +#========================================================== +def gcdist(x1,y1,x2,y2,Rearth,units='degrees'): + """ + Compute the great circle distance on the earth between points + (x1,y1) and (x2,y2), where: + x = longitude, y = latitude + + Taken from Clawpack 4.x. + """ + from numpy import pi,sin,cos,arccos,arcsin,sqrt + if units=='degrees': + # convert to radians: + x1 = x1 * pi/180. + y1 = y1 * pi/180. + x2 = x2 * pi/180. + y2 = y2 * pi/180. + elif units != 'radians': + raise Exception("unrecognized units") + + dx = x1 - x2 + dy = y1 - y2 + + # angle subtended by two points, using Haversine formula: + dsigma = 2. * arcsin(sqrt(sin(0.5*dy)**2 + cos(y1)*cos(y2)*sin(0.5*dx)**2)) + + # alternative formula that may have more rounding error: + #dsigma2 = arccos(sin(y1)*sin(y2)+ cos(y1)*cos(y2)*cos(dx)) + #print("max diff in dsigma: ", abs(dsigma-dsigma2).max()) + + d = Rearth * dsigma + return d + + diff --git a/tsunami/radial-ocean/awr2011_5.x/plot_ocean.py b/tsunami/radial-ocean/awr2011_5.x/plot_ocean.py new file mode 100644 index 00000000..b8baa660 --- /dev/null +++ b/tsunami/radial-ocean/awr2011_5.x/plot_ocean.py @@ -0,0 +1,90 @@ + +from __future__ import print_function +from pylab import * + +from clawpack.geoclaw import topotools +from clawpack.geoclaw.data import Rearth # radius of earth +from clawpack.clawutil.data import ClawData +from mapper import latlong +import maketopo + +probdata = ClawData() +probdata.read('setprob.data', force=True) +theta_island = probdata.theta_island +print("theta_island = ",theta_island) + + +#close(1) +#figure(1, figsize=(10,6)) +#axes([.1,.1,.5,.8]) +clf() +s = linspace(0, 360., 1001) +xshore,yshore = latlong(1645.e3, s, 40., Rearth) +plot(xshore,yshore,'k',linewidth=2) + +xshelf,yshelf = latlong(1560.e3, s, 40., Rearth) +plot(xshelf,yshelf,'k--',linewidth=2) + +theta_island = 220. +(xi,yi) = latlong(1600.e3, theta_island, 40., Rearth) +xb = [xi-1, xi+1, xi+1, xi-1, xi-1] +yb = [yi-1, yi-1, yi+1, yi+1, yi-1] +plot(xb,yb,'b',linewidth=1.5) +text(5,49,'Test 1', fontsize=20) + +theta_island = 260. +(xi,yi) = latlong(1600.e3, theta_island, 40., Rearth) +xb = [xi-1, xi+1, xi+1, xi-1, xi-1] +yb = [yi-1, yi-1, yi+1, yi+1, yi-1] +plot(xb,yb,'b',linewidth=1.5) +text(9,40,'Test 2', fontsize=20) + +title('(a) Radial ocean',fontsize=20) + +x = linspace(-5,5,51) +y = linspace(35,45,51) +X,Y = meshgrid(x,y) +Z = maketopo.qinit(X,Y) +contour(X,Y,Z,[2],linewidth=2,colors='k') + +axis('scaled') +xlim([-21,21]) +ylim([20,60]) +xticks(fontsize='15') +yticks([25,35,45,55],fontsize='15') +xlabel('Longitude',fontsize=15) +ylabel('Latitude',fontsize=15) + +savefig('ocean.png') +print("Created ocean.png") + +if 0: + axes([.65,.6,.25,.25]) + plot(xshore,yshore,'k',linewidth=2) + plot(xshelf,yshelf,'k--',linewidth=2) + maketopo.theta_island = 220. + (xisland,yisland) = latlong(1600.e3, maketopo.theta_island, 40., Rearth) + x = linspace(xisland-1, xisland+1, 200) + y = linspace(yisland-1, yisland+1, 200) + X,Y = meshgrid(x,y) + Z = maketopo.topo(X,Y) + contour(X,Y,Z,[0,7,15]) + axis('scaled') + xlim([xisland-0.5, xisland+0.5]) + ylim([yisland-0.5, yisland+0.5]) + + + axes([.65,.2,.25,.25]) + plot(xshore,yshore,'k',linewidth=2) + plot(xshelf,yshelf,'k--',linewidth=2) + maketopo.theta_island = 260. + (xisland,yisland) = latlong(1600.e3, maketopo.theta_island, 40., Rearth) + x = linspace(xisland-1, xisland+1, 200) + y = linspace(yisland-1, yisland+1, 200) + X,Y = meshgrid(x,y) + Z = maketopo.topo(X,Y) + contour(X,Y,Z,[0,7,15]) + axis('scaled') + xlim([xisland-0.5, xisland+0.5]) + ylim([yisland-0.5, yisland+0.5]) + diff --git a/tsunami/radial-ocean/awr2011_5.x/plot_xsec.py b/tsunami/radial-ocean/awr2011_5.x/plot_xsec.py new file mode 100644 index 00000000..f953a7e4 --- /dev/null +++ b/tsunami/radial-ocean/awr2011_5.x/plot_xsec.py @@ -0,0 +1,74 @@ + +from __future__ import print_function +from clawpack.geoclaw import topotools +from clawpack.geoclaw.data import Rearth # radius of earth +from interp import pwcubic +from clawpack.clawutil.data import ClawData +import numpy as np +from matplotlib import pylab as plt + +from mapper import latlong +import maketopo + +probdata = ClawData() +probdata.read('setprob.data', force=True) +theta_island = probdata.theta_island +print("theta_island = ",theta_island) + + +(xisland,yisland) = latlong(1600.e3, theta_island, 40., Rearth) + + +def ocean(): + dmax = 1650.e3 + dp = np.linspace(0., dmax, 3301) + zp = maketopo.shelf1(dp) + plt.clf() + plt.plot(dp*1e-3,zp,'k-',linewidth=3) + plt.fill_between(dp*1e-3,zp,0.,where=(zp<0.), color=[0,0.5,1]) + plt.xlim([0.,1650]) + plt.ylim([-4500.,500]) + plt.title("Topography as function of radius",fontsize=18) + plt.xlabel("kilometers from center",fontsize=15) + plt.ylabel("meters",fontsize=15) + plt.xticks(fontsize=15) + plt.yticks(fontsize=15) + plt.savefig('topo1.png') + #plt.savefig('topo1.tif') + print("Cross-section plot in topo1.png") + + plt.xlim([1500,1650]) + plt.ylim([-200.,20]) + plt.title("Topography of shelf and beach",fontsize=18) + plt.xlabel("kilometers from center",fontsize=15) + plt.ylabel("meters",fontsize=15) + plt.xticks(fontsize=15) + plt.yticks(fontsize=15) + plt.savefig('topo2.png') + #plt.savefig('topo2.tif') + print("Zoom near shore in topo2.png") + + +def with_island(): + dmax = 1650.e3 + dp = np.linspace(1500e3, dmax, 2000) + zp = maketopo.shelf1(dp) + maketopo.island1(abs(dp-1600e3)) + + plt.clf() + plt.plot(dp*1e-3,zp,'k-',linewidth=3) + plt.fill_between(dp*1e-3,zp,0.,where=(zp<0.), color=[0,0.5,1]) + + plt.xlim([1500,1650]) + plt.ylim([-200.,30]) + plt.title("Cross-section through island center",fontsize=18) + plt.xlabel("kilometers from center",fontsize=15) + plt.ylabel("meters",fontsize=15) + plt.xticks(fontsize=15) + plt.yticks(fontsize=15) + plt.savefig('topo3.png') + #plt.savefig('topo3.tif') + print("Cross-section through island in topo3.png") + +if __name__=="__main__": + ocean() + with_island() diff --git a/tsunami/radial-ocean/awr2011_5.x/setplot.py b/tsunami/radial-ocean/awr2011_5.x/setplot.py new file mode 100644 index 00000000..e90f45ff --- /dev/null +++ b/tsunami/radial-ocean/awr2011_5.x/setplot.py @@ -0,0 +1,240 @@ + +""" +Set up the plot figures, axes, and items to be done for each frame. + +This module is imported by the plotting routines and then the +function setplot is called to set the plot parameters. + +""" + +from clawpack.geoclaw import topotools +from clawpack.geoclaw.data import Rearth # radius of earth +from clawpack.clawutil.data import ClawData + +from mapper import latlong, gcdist + +probdata = ClawData() +probdata.read('setprob.data', force=True) +theta_island = probdata.theta_island + + +#-------------------------- +def setplot(plotdata): +#-------------------------- + + """ + Specify what is to be plotted at each frame. + Input: plotdata, an instance of clawpack.visclaw.data.ClawPlotData. + Output: a modified version of plotdata. + + """ + + + from clawpack.visclaw import colormaps, geoplot + from numpy import linspace + + plotdata.clearfigures() # clear any old figures,axes,items data + plotdata.format = 'binary' + + + # To plot gauge locations on pcolor or contour plot, use this as + # an afteraxis function: + + def addgauges(current_data): + from clawpack.visclaw import gaugetools + gaugetools.plot_gauge_locations(current_data.plotdata, \ + gaugenos='all', format_string='ko', add_labels=True) + + + #----------------------------------------- + # Figure for pcolor plot + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='pcolor', figno=0) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.title = 'Surface' + plotaxes.scaled = True + #plotaxes.afteraxes = addgauges + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + #plotitem.plot_var = geoplot.surface + plotitem.plot_var = geoplot.surface_or_depth + plotitem.pcolor_cmap = geoplot.tsunami_colormap + plotitem.pcolor_cmin = -1.0 + plotitem.pcolor_cmax = 1.0 + plotitem.add_colorbar = True + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 1 + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.land + plotitem.pcolor_cmap = colormaps.all_white # to print better in B&W + #plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 100.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [1,0] + plotitem.patchedges_show = 1 + plotaxes.xlimits = [-20,20] + plotaxes.ylimits = [20,60] + + + # add contour lines of bathy if desired: + plotitem = plotaxes.new_plotitem(plot_type='2d_contour') + plotitem.show = False + plotitem.plot_var = geoplot.topo + plotitem.contour_levels = linspace(-1000,-1000,1) + plotitem.amr_contour_colors = ['g'] # color on each level + plotitem.kwargs = {'linestyles':'dashed','linewidths':2} + plotitem.amr_contour_show = [1,0,0] + plotitem.celledges_show = 0 + plotitem.patchedges_show = 0 + + + #----------------------------------------- + # Figure for zoom + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Zoom1', figno=7) + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('zoom on island') + plotaxes.title = 'Surface elevation' + plotaxes.scaled = True + xisland,yisland = latlong(1600e3, theta_island, 40., Rearth) + plotaxes.xlimits = [xisland-0.6, xisland+0.6] + plotaxes.ylimits = [yisland-0.6, yisland+0.6] + def bigfont(current_data): + import pylab + t = current_data.t + addgauges(current_data) + pylab.title("Surface at t = %8.1f" % t, fontsize=20) + pylab.xticks(fontsize=15) + pylab.yticks(fontsize=15) + plotaxes.afteraxes = bigfont + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.show = True + #plotitem.plot_var = geoplot.surface + plotitem.plot_var = geoplot.surface_or_depth + plotitem.pcolor_cmap = geoplot.tsunami_colormap + plotitem.pcolor_cmin = -1.0 + plotitem.pcolor_cmax = 1.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0] + plotitem.patchedges_show = 1 + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.show = True + plotitem.plot_var = geoplot.land + plotitem.pcolor_cmap = colormaps.all_white # to print better in B&W + #plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 100.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [1,0,0,0,0] + plotitem.patchedges_show = 1 + + # contour lines: + plotitem = plotaxes.new_plotitem(plot_type='2d_contour') + plotitem.show = True + plotitem.plot_var = geoplot.surface + plotitem.contour_levels = [-0.8, -0.4, 0.4, 0.8] + plotitem.amr_contour_colors = ['k'] # color on each level + plotitem.kwargs = {'linewidths':2} + plotitem.amr_contour_show = [0,0,0,1,1] + plotitem.celledges_show = 0 + plotitem.patchedges_show = 0 + + # add contour lines of bathy if desired: + plotitem = plotaxes.new_plotitem(plot_type='2d_contour') + plotitem.show = False + plotitem.plot_var = geoplot.topo + plotitem.contour_levels = linspace(-1000,-1000,1) + plotitem.amr_contour_colors = ['g'] # color on each level + plotitem.kwargs = {'linestyles':'dashed','linewidths':2} + plotitem.amr_contour_show = [0,0,1] + plotitem.celledges_show = 0 + plotitem.patchedges_show = 0 + + #----------------------------------------- + # Figures for gauges + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Surface & topo', figno=300, \ + type='each_gauge') + plotfigure.clf_each_gauge = True + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = [6000, 14000] + plotaxes.ylimits = [-2.0, 3.0] + plotaxes.title = 'Surface' + + # Plot surface as blue curve: + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.plot_var = 3 + plotitem.plotstyle = 'b-' + + # Plot topo as green curve: + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + + def gaugetopo(current_data): + q = current_data.q + h = q[0,:] + eta = q[3,:] + topo = eta - h + return topo + + plotitem.plot_var = gaugetopo + plotitem.plotstyle = 'g-' + def add_zeroline(current_data): + from pylab import plot, legend + t = current_data.t + #legend(('surface','topography'),loc='lower left') + plot(t, 0*t, 'k') + + plotaxes.afteraxes = add_zeroline + + + #----------------------------------------- + # Figure for grids alone + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='grids', figno=2) + plotfigure.show = False + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = [0,1] + plotaxes.ylimits = [0,1] + plotaxes.title = 'grids' + plotaxes.scaled = True + + # Set up for item on these axes: + plotitem = plotaxes.new_plotitem(plot_type='2d_patch') + plotitem.amr_patch_bgcolor = ['#ffeeee', '#eeeeff', '#eeffee'] + plotitem.amr_celledges_show = [1,1,0] + plotitem.amr_patchedges_show = [1] + + #----------------------------------------- + + # Parameters used only when creating html and/or latex hardcopy + # e.g., via clawpack.visclaw.frametools.printframes: + + plotdata.printfigs = True # print figures + plotdata.print_format = 'png' # file format + plotdata.print_framenos = 'all' # list of frames to print + plotdata.print_gaugenos = 'all' # list of gauges to print + plotdata.print_fignos = 'all' # list of figures to print + plotdata.html = True # create html files of plots? + plotdata.html_homelink = '../README.html' # pointer for top of index + plotdata.latex = True # create latex file of plots? + plotdata.latex_figsperline = 2 # layout of plots + plotdata.latex_framesperline = 1 # layout of plots + plotdata.latex_makepdf = False # also run pdflatex? + + return plotdata + + diff --git a/tsunami/radial-ocean/awr2011_5.x/setrun.py b/tsunami/radial-ocean/awr2011_5.x/setrun.py new file mode 100644 index 00000000..ca0f91d3 --- /dev/null +++ b/tsunami/radial-ocean/awr2011_5.x/setrun.py @@ -0,0 +1,462 @@ +""" +Module to set up run time parameters for Clawpack -- AMRClaw code. + +The values set in the function setrun are then written out to data files +that will be read in by the Fortran code. + +""" + +from __future__ import print_function +import os +import numpy as np + +from clawpack.geoclaw.data import Rearth # radius of earth +from mapper import latlong + + +#------------------------------ +def setrun(claw_pkg='geoclaw'): +#------------------------------ + + """ + Define the parameters used for running Clawpack. + + INPUT: + claw_pkg expected to be "geoclaw" for this setrun. + + OUTPUT: + rundata - object of class ClawRunData + + """ + + from clawpack.clawutil import data + + + assert claw_pkg.lower() == 'geoclaw', "Expected claw_pkg = 'geoclaw'" + + num_dim = 2 + rundata = data.ClawRunData(claw_pkg, num_dim) + + + #------------------------------------------------------------------ + # Problem-specific parameters to be written to setprob.data: + #------------------------------------------------------------------ + theta_island = 220. + probdata = rundata.new_UserData(name='probdata',fname='setprob.data') + probdata.add_param('theta_island', theta_island, 'angle to island') + + + #------------------------------------------------------------------ + # GeoClaw specific parameters: + #------------------------------------------------------------------ + rundata = setgeo(rundata) + + #------------------------------------------------------------------ + # Standard Clawpack parameters to be written to claw.data: + #------------------------------------------------------------------ + + clawdata = rundata.clawdata # initialized when rundata instantiated + + + # Set single grid parameters first. + # See below for AMR parameters. + + + # --------------- + # Spatial domain: + # --------------- + + # Number of space dimensions: + clawdata.num_dim = num_dim + + # Lower and upper edge of computational domain: + clawdata.lower[0] = -20.0 # xlower + clawdata.upper[0] = 20.0 # xupper + clawdata.lower[1] = 20.0 # ylower + clawdata.upper[1] = 60.0 # yupper + + # Number of grid cells: + clawdata.num_cells[0] = 40 # mx + clawdata.num_cells[1] = 40 # my + + + # --------------- + # Size of system: + # --------------- + + # Number of equations in the system: + clawdata.num_eqn = 3 + + # Number of auxiliary variables in the aux array (initialized in setaux) + clawdata.num_aux = 3 + + # Index of aux array corresponding to capacity function, if there is one: + clawdata.capa_index = 2 + + + # ------------- + # Initial time: + # ------------- + + clawdata.t0 = 0.0 + + + # Restart from checkpoint file of a previous run? + # Note: If restarting, you must also change the Makefile to set: + # RESTART = True + # If restarting, t0 above should be from original run, and the + # restart_file 'fort.chkNNNNN' specified below should be in + # the OUTDIR indicated in Makefile. + + clawdata.restart = False # True to restart from prior results + clawdata.restart_file = 'fort.chk00006' # File to use for restart data + + + # ------------- + # Output times: + #-------------- + + # Specify at what times the results should be written to fort.q files. + # Note that the time integration stops after the final output time. + + clawdata.output_style = 1 + + if clawdata.output_style==1: + # Output ntimes frames at equally spaced times up to tfinal: + # Can specify num_output_times = 0 for no output + clawdata.num_output_times = 14 + clawdata.tfinal = 14000.0 + clawdata.output_t0 = True # output at initial (or restart) time? + + elif clawdata.output_style == 2: + # Specify a list or numpy array of output times: + # Include t0 if you want output at the initial time. + clawdata.output_times = [0., 0.1] + + elif clawdata.output_style == 3: + # Output every step_interval timesteps over total_steps timesteps: + clawdata.output_step_interval = 2 + clawdata.total_steps = 4 + clawdata.output_t0 = True # output at initial (or restart) time? + + + clawdata.output_format = 'binary' # 'ascii', 'binary' + + clawdata.output_q_components = 'all' # could be list such as [True,True] + clawdata.output_aux_components = 'none' # could be list + clawdata.output_aux_onlyonce = True # output aux arrays only at t0 + + + # --------------------------------------------------- + # Verbosity of messages to screen during integration: + # --------------------------------------------------- + + # The current t, dt, and cfl will be printed every time step + # at AMR levels <= verbosity. Set verbosity = 0 for no printing. + # (E.g. verbosity == 2 means print only on levels 1 and 2.) + clawdata.verbosity = 1 + + + + # -------------- + # Time stepping: + # -------------- + + # if dt_variable==True: variable time steps used based on cfl_desired, + # if dt_variable==Falseixed time steps dt = dt_initial always used. + clawdata.dt_variable = True + + # Initial time step for variable dt. + # (If dt_variable==0 then dt=dt_initial for all steps) + clawdata.dt_initial = 16.0 + + # Max time step to be allowed if variable dt used: + clawdata.dt_max = 1e+99 + + # Desired Courant number if variable dt used + clawdata.cfl_desired = 0.75 + # max Courant number to allow without retaking step with a smaller dt: + clawdata.cfl_max = 1.0 + + # Maximum number of time steps to allow between output times: + clawdata.steps_max = 5000 + + + # ------------------ + # Method to be used: + # ------------------ + + # Order of accuracy: 1 => Godunov, 2 => Lax-Wendroff plus limiters + clawdata.order = 2 + + # Use dimensional splitting? (not yet available for AMR) + clawdata.dimensional_split = 'unsplit' + + # For unsplit method, transverse_waves can be + # 0 or 'none' ==> donor cell (only normal solver used) + # 1 or 'increment' ==> corner transport of waves + # 2 or 'all' ==> corner transport of 2nd order corrections too + clawdata.transverse_waves = 2 + + + # Number of waves in the Riemann solution: + clawdata.num_waves = 3 + + # List of limiters to use for each wave family: + # Required: len(limiter) == num_waves + # Some options: + # 0 or 'none' ==> no limiter (Lax-Wendroff) + # 1 or 'minmod' ==> minmod + # 2 or 'superbee' ==> superbee + # 3 or 'vanleer' ==> van Leer + # 4 or 'mc' ==> MC limiter + clawdata.limiter = ['vanleer', 'vanleer', 'vanleer'] + + clawdata.use_fwaves = True # True ==> use f-wave version of algorithms + + # Source terms splitting: + # src_split == 0 or 'none' ==> no source term (src routine never called) + # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used, + # src_split == 2 or 'strang' ==> Strang (2nd order) splitting used, not recommended. + clawdata.source_split = 1 + + + # -------------------- + # Boundary conditions: + # -------------------- + + # Number of ghost cells (usually 2) + clawdata.num_ghost = 2 + + # Choice of BCs at xlower and xupper: + # 0 or 'user' => user specified (must modify bcNamr.f to use this option) + # 1 or 'extrap' => extrapolation (non-reflecting outflow) + # 2 or 'periodic' => periodic (must specify this at both boundaries) + # 3 or 'wall' => solid wall for systems where q(2) is normal velocity + + clawdata.bc_lower[0] = 'extrap' # at xlower + clawdata.bc_upper[0] = 'extrap' # at xupper + + clawdata.bc_lower[1] = 'extrap' # at ylower + clawdata.bc_upper[1] = 'extrap' # at yupper + + + # --------------- + # Gauges: + # --------------- + gauges = rundata.gaugedata.gauges + # for gauges append lines of the form [gaugeno, x, y, t1, t2] + + gaugeno = 0 + for d in np.linspace(1550e3,1650e3,11): + gaugeno = gaugeno+1 + x,y = latlong(d, theta_island, 40., Rearth) + gauges.append([gaugeno, x, y, 0., 1e10]) + + + # -------------- + # Checkpointing: + # -------------- + + # Specify when checkpoint files should be created that can be + # used to restart a computation. + + clawdata.checkpt_style = 0 + + if clawdata.checkpt_style == 0: + # Do not checkpoint at all + pass + + elif clawdata.checkpt_style == 1: + # Checkpoint only at tfinal. + pass + + elif clawdata.checkpt_style == 2: + # Specify a list of checkpoint times. + clawdata.checkpt_times = [0.1,0.15] + + elif clawdata.checkpt_style == 3: + # Checkpoint every checkpt_interval timesteps (on Level 1) + # and at the final time. + clawdata.checkpt_interval = 5 + + + + # --------------- + # AMR parameters: (written to amr.data) + # --------------- + amrdata = rundata.amrdata + + # max number of refinement levels: + amrdata.amr_levels_max = 5 + + # List of refinement ratios at each level (length at least amr_level_max-1) + amrdata.refinement_ratios_x = [4, 4, 4, 4] + amrdata.refinement_ratios_y = [4, 4, 4, 4] + amrdata.refinement_ratios_t = [4, 4, 4, 1] + + + # Specify type of each aux variable in amrdata.auxtype. + # This must be a list of length num_aux, each element of which is one of: + # 'center', 'capacity', 'xleft', or 'yleft' (see documentation). + amrdata.aux_type = ['center', 'capacity', 'yleft'] + + + # Flag for refinement based on Richardson error estimater: + amrdata.flag_richardson = False # use Richardson? + amrdata.flag_richardson_tol = 1.0 # Richardson tolerance + + # Flag for refinement using routine flag2refine: + amrdata.flag2refine = True # use this? + amrdata.flag2refine_tol = 0.5 # tolerance used in this routine + # Note: in geoclaw the refinement tolerance is set as wave_tolerance below + # and flag2refine_tol is unused! + + # steps to take on each level L between regriddings of level L+1: + amrdata.regrid_interval = 3 + + # width of buffer zone around flagged points: + # (typically the same as regrid_interval so waves don't escape): + amrdata.regrid_buffer_width = 2 + + # clustering alg. cutoff for (# flagged pts) / (total # of cells refined) + # (closer to 1.0 => more small grids may be needed to cover flagged cells) + amrdata.clustering_cutoff = 0.7 + + # print info about each regridding up to this level: + amrdata.verbosity_regrid = 0 + + + # --------------- + # Regions: + # --------------- + regions = rundata.regiondata.regions + # to specify regions of refinement append lines of the form + # [minlevel,maxlevel,t1,t2,x1,x2,y1,y2] + + # allow 2 levels over portion of ocean in correct direction: + # and 3 levels up to time 7000. + if theta_island == 190: + regions.append([1, 2, 0.e0, 1.e10, -5., 5., 35., 55.]) + regions.append([1, 3, 0.e0, 7000., -5., 5., 35., 55.]) + if theta_island == 220: + regions.append([1, 2, 0.e0, 7000., -5., 20., 35., 55.]) + regions.append([1, 3, 0.e0, 7000., -5., 20., 35., 55.]) + regions.append([1, 2, 7000., 9000., 10., 20., 45., 52.]) + if theta_island == 260: + regions.append([1, 2, 0.e0, 1.e10, -5., 20., 35., 45.]) + regions.append([1, 3, 0.e0, 7000., -5., 20., 35., 45.]) + + # Force refinement near the island as the wave approaches: + + (xisland,yisland) = latlong(1600.e3, theta_island, 40., Rearth) + x1 = xisland - 1. + x2 = xisland + 1. + y1 = yisland - 1. + y2 = yisland + 1. + regions.append([4, 4, 7000., 1.e10, x1,x2,y1,y2]) + + x1 = xisland - 0.2 + x2 = xisland + 0.2 + y1 = yisland - 0.2 + y2 = yisland + 0.2 + regions.append([4, 5, 8000., 1.e10, x1,x2,y1,y2]) + + + + # ----- For developers ----- + # Toggle debugging print statements: + amrdata.dprint = False # print domain flags + amrdata.eprint = False # print err est flags + amrdata.edebug = False # even more err est flags + amrdata.gprint = False # grid bisection/clustering + amrdata.nprint = False # proper nesting output + amrdata.pprint = False # proj. of tagged points + amrdata.rprint = False # print regridding summary + amrdata.sprint = False # space/memory output + amrdata.tprint = False # time step reporting each level + amrdata.uprint = False # update/upbnd reporting + + return rundata + + # end of function setrun + # ---------------------- + + +#------------------- +def setgeo(rundata): +#------------------- + """ + Set GeoClaw specific runtime parameters. + """ + + try: + geo_data = rundata.geo_data + except: + print("*** Error, this rundata has no geo_data attribute") + raise AttributeError("Missing geo_data attribute") + + # == Physics == + geo_data.gravity = 9.81 + geo_data.coordinate_system = 2 + geo_data.earth_radius = 6367500.0 + + # == Forcing Options + geo_data.coriolis_forcing = False + + # == Algorithm and Initial Conditions == + geo_data.sea_level = 0.0 + geo_data.dry_tolerance = 0.001 + geo_data.friction_forcing = True + geo_data.manning_coefficient = 0.025 + geo_data.friction_depth = 1000000.0 + + # Refinement settings + refinement_data = rundata.refinement_data + refinement_data.variable_dt_refinement_ratios = False + refinement_data.wave_tolerance = 0.01 + refinement_data.deep_depth = 100.0 + refinement_data.max_level_deep = 3 + + # == settopo.data values == + rundata.topo_data.topofiles = [[2, 1, 1, 0.0, 10000000000.0, 'ocean.topotype2'], [2, 3, 4, 7.0, 10000000000.0, 'island.topotype2']] + topofiles = rundata.topo_data.topofiles + # for topography, append lines of the form + # [topotype, minlevel, maxlevel, t1, t2, fname] + + # == setdtopo.data values == + rundata.dtopo_data.dtopofiles = [] + dtopofiles = rundata.dtopo_data.dtopofiles + # for moving topography, append lines of the form : + # [topotype, minlevel,maxlevel,fname] + + + # == setqinit.data values == + rundata.qinit_data.qinit_type = 4 + rundata.qinit_data.qinitfiles = [[1, 2, 'hump.xyz']] + qinitfiles = rundata.qinit_data.qinitfiles + # for qinit perturbations, append lines of the form: (<= 1 allowed for now!) + # [minlev, maxlev, fname] + + # == fixedgrids.data values == + rundata.fixed_grid_data.fixedgrids = [] + fixedgrids = rundata.fixed_grid_data.fixedgrids + # for fixed grids append lines of the form + # [t1,t2,noutput,x1,x2,y1,y2,xpoints,ypoints,\ + # ioutarrivaltimes,ioutsurfacemax] + + # == fgmax.data values == + fgmax_files = rundata.fgmax_data.fgmax_files + # for fixed grids append to this list names of any fgmax input files + + + return rundata + # end of function setgeo + # ---------------------- + +if __name__ == '__main__': + # Set up run-time parameters and write all data files. + import sys + rundata = setrun(*sys.argv[1:]) + rundata.write() + diff --git a/tsunami/radial-ocean/awr2011_5.x/setrun_220.py b/tsunami/radial-ocean/awr2011_5.x/setrun_220.py new file mode 100644 index 00000000..e7af2b55 --- /dev/null +++ b/tsunami/radial-ocean/awr2011_5.x/setrun_220.py @@ -0,0 +1,454 @@ +""" +Module to set up run time parameters for Clawpack -- AMRClaw code. + +The values set in the function setrun are then written out to data files +that will be read in by the Fortran code. + +""" + +from __future__ import print_function +import os +import numpy as np + +from clawpack.geoclaw.data import Rearth # radius of earth +from mapper import latlong + + +#------------------------------ +def setrun(claw_pkg='geoclaw'): +#------------------------------ + + """ + Define the parameters used for running Clawpack. + + INPUT: + claw_pkg expected to be "geoclaw" for this setrun. + + OUTPUT: + rundata - object of class ClawRunData + + """ + + from clawpack.clawutil import data + + + assert claw_pkg.lower() == 'geoclaw', "Expected claw_pkg = 'geoclaw'" + + num_dim = 2 + rundata = data.ClawRunData(claw_pkg, num_dim) + + + #------------------------------------------------------------------ + # Problem-specific parameters to be written to setprob.data: + #------------------------------------------------------------------ + theta_island = 220. + probdata = rundata.new_UserData(name='probdata',fname='setprob.data') + probdata.add_param('theta_island', theta_island, 'angle to island') + + + #------------------------------------------------------------------ + # GeoClaw specific parameters: + #------------------------------------------------------------------ + rundata = setgeo(rundata) + + #------------------------------------------------------------------ + # Standard Clawpack parameters to be written to claw.data: + #------------------------------------------------------------------ + + clawdata = rundata.clawdata # initialized when rundata instantiated + + + # Set single grid parameters first. + # See below for AMR parameters. + + + # --------------- + # Spatial domain: + # --------------- + + # Number of space dimensions: + clawdata.num_dim = num_dim + + # Lower and upper edge of computational domain: + clawdata.lower[0] = -20.0 # xlower + clawdata.upper[0] = 20.0 # xupper + clawdata.lower[1] = 20.0 # ylower + clawdata.upper[1] = 60.0 # yupper + + # Number of grid cells: + clawdata.num_cells[0] = 40 # mx + clawdata.num_cells[1] = 40 # my + + + # --------------- + # Size of system: + # --------------- + + # Number of equations in the system: + clawdata.num_eqn = 3 + + # Number of auxiliary variables in the aux array (initialized in setaux) + clawdata.num_aux = 3 + + # Index of aux array corresponding to capacity function, if there is one: + clawdata.capa_index = 2 + + + # ------------- + # Initial time: + # ------------- + + clawdata.t0 = 0.0 + + + # Restart from checkpoint file of a previous run? + # Note: If restarting, you must also change the Makefile to set: + # RESTART = True + # If restarting, t0 above should be from original run, and the + # restart_file 'fort.chkNNNNN' specified below should be in + # the OUTDIR indicated in Makefile. + + clawdata.restart = False # True to restart from prior results + clawdata.restart_file = 'fort.chk00006' # File to use for restart data + + + # ------------- + # Output times: + #-------------- + + # Specify at what times the results should be written to fort.q files. + # Note that the time integration stops after the final output time. + + clawdata.output_style = 1 + + if clawdata.output_style==1: + # Output ntimes frames at equally spaced times up to tfinal: + # Can specify num_output_times = 0 for no output + clawdata.num_output_times = 14 + clawdata.tfinal = 14000.0 + clawdata.output_t0 = True # output at initial (or restart) time? + + elif clawdata.output_style == 2: + # Specify a list or numpy array of output times: + # Include t0 if you want output at the initial time. + clawdata.output_times = [0., 0.1] + + elif clawdata.output_style == 3: + # Output every step_interval timesteps over total_steps timesteps: + clawdata.output_step_interval = 2 + clawdata.total_steps = 4 + clawdata.output_t0 = True # output at initial (or restart) time? + + + clawdata.output_format = 'binary' # 'ascii', 'binary' + + clawdata.output_q_components = 'all' # could be list such as [True,True] + clawdata.output_aux_components = 'none' # could be list + clawdata.output_aux_onlyonce = True # output aux arrays only at t0 + + + # --------------------------------------------------- + # Verbosity of messages to screen during integration: + # --------------------------------------------------- + + # The current t, dt, and cfl will be printed every time step + # at AMR levels <= verbosity. Set verbosity = 0 for no printing. + # (E.g. verbosity == 2 means print only on levels 1 and 2.) + clawdata.verbosity = 1 + + + + # -------------- + # Time stepping: + # -------------- + + # if dt_variable==True: variable time steps used based on cfl_desired, + # if dt_variable==Falseixed time steps dt = dt_initial always used. + clawdata.dt_variable = True + + # Initial time step for variable dt. + # (If dt_variable==0 then dt=dt_initial for all steps) + clawdata.dt_initial = 16.0 + + # Max time step to be allowed if variable dt used: + clawdata.dt_max = 1e+99 + + # Desired Courant number if variable dt used + clawdata.cfl_desired = 0.75 + # max Courant number to allow without retaking step with a smaller dt: + clawdata.cfl_max = 1.0 + + # Maximum number of time steps to allow between output times: + clawdata.steps_max = 5000 + + + # ------------------ + # Method to be used: + # ------------------ + + # Order of accuracy: 1 => Godunov, 2 => Lax-Wendroff plus limiters + clawdata.order = 2 + + # Use dimensional splitting? (not yet available for AMR) + clawdata.dimensional_split = 'unsplit' + + # For unsplit method, transverse_waves can be + # 0 or 'none' ==> donor cell (only normal solver used) + # 1 or 'increment' ==> corner transport of waves + # 2 or 'all' ==> corner transport of 2nd order corrections too + clawdata.transverse_waves = 2 + + + # Number of waves in the Riemann solution: + clawdata.num_waves = 3 + + # List of limiters to use for each wave family: + # Required: len(limiter) == num_waves + # Some options: + # 0 or 'none' ==> no limiter (Lax-Wendroff) + # 1 or 'minmod' ==> minmod + # 2 or 'superbee' ==> superbee + # 3 or 'vanleer' ==> van Leer + # 4 or 'mc' ==> MC limiter + clawdata.limiter = ['vanleer', 'vanleer', 'vanleer'] + + clawdata.use_fwaves = True # True ==> use f-wave version of algorithms + + # Source terms splitting: + # src_split == 0 or 'none' ==> no source term (src routine never called) + # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used, + # src_split == 2 or 'strang' ==> Strang (2nd order) splitting used, not recommended. + clawdata.source_split = 1 + + + # -------------------- + # Boundary conditions: + # -------------------- + + # Number of ghost cells (usually 2) + clawdata.num_ghost = 2 + + # Choice of BCs at xlower and xupper: + # 0 or 'user' => user specified (must modify bcNamr.f to use this option) + # 1 or 'extrap' => extrapolation (non-reflecting outflow) + # 2 or 'periodic' => periodic (must specify this at both boundaries) + # 3 or 'wall' => solid wall for systems where q(2) is normal velocity + + clawdata.bc_lower[0] = 'extrap' # at xlower + clawdata.bc_upper[0] = 'extrap' # at xupper + + clawdata.bc_lower[1] = 'extrap' # at ylower + clawdata.bc_upper[1] = 'extrap' # at yupper + + + # --------------- + # Gauges: + # --------------- + gauges = rundata.gaugedata.gauges + # for gauges append lines of the form [gaugeno, x, y, t1, t2] + + gaugeno = 0 + for d in [1570e3, 1590e3, 1610e3, 1630e3]: + gaugeno = gaugeno+1 + x,y = latlong(d, theta_island, 40., Rearth) + gauges.append([gaugeno, x, y, 0., 1e10]) + + + # -------------- + # Checkpointing: + # -------------- + + # Specify when checkpoint files should be created that can be + # used to restart a computation. + + clawdata.checkpt_style = 0 + + if clawdata.checkpt_style == 0: + # Do not checkpoint at all + pass + + elif clawdata.checkpt_style == 1: + # Checkpoint only at tfinal. + pass + + elif clawdata.checkpt_style == 2: + # Specify a list of checkpoint times. + clawdata.checkpt_times = [0.1,0.15] + + elif clawdata.checkpt_style == 3: + # Checkpoint every checkpt_interval timesteps (on Level 1) + # and at the final time. + clawdata.checkpt_interval = 5 + + + + # --------------- + # AMR parameters: (written to amr.data) + # --------------- + amrdata = rundata.amrdata + + # max number of refinement levels: + amrdata.amr_levels_max = 5 + + # List of refinement ratios at each level (length at least amr_level_max-1) + amrdata.refinement_ratios_x = [4, 4, 4, 4] + amrdata.refinement_ratios_y = [4, 4, 4, 4] + amrdata.refinement_ratios_t = [4, 4, 4, 1] + + + # Specify type of each aux variable in amrdata.auxtype. + # This must be a list of length num_aux, each element of which is one of: + # 'center', 'capacity', 'xleft', or 'yleft' (see documentation). + amrdata.aux_type = ['center', 'capacity', 'yleft'] + + + # Flag for refinement based on Richardson error estimater: + amrdata.flag_richardson = False # use Richardson? + amrdata.flag_richardson_tol = 1.0 # Richardson tolerance + + # Flag for refinement using routine flag2refine: + amrdata.flag2refine = True # use this? + amrdata.flag2refine_tol = 0.5 # tolerance used in this routine + # Note: in geoclaw the refinement tolerance is set as wave_tolerance below + # and flag2refine_tol is unused! + + # steps to take on each level L between regriddings of level L+1: + amrdata.regrid_interval = 3 + + # width of buffer zone around flagged points: + # (typically the same as regrid_interval so waves don't escape): + amrdata.regrid_buffer_width = 2 + + # clustering alg. cutoff for (# flagged pts) / (total # of cells refined) + # (closer to 1.0 => more small grids may be needed to cover flagged cells) + amrdata.clustering_cutoff = 0.7 + + # print info about each regridding up to this level: + amrdata.verbosity_regrid = 0 + + + # --------------- + # Regions: + # --------------- + regions = rundata.regiondata.regions + # to specify regions of refinement append lines of the form + # [minlevel,maxlevel,t1,t2,x1,x2,y1,y2] + + regions.append([1, 3, 0., 5000., -5., 20., 35., 55.]) + regions.append([1, 2, 5000., 6900., -5., 20., 35., 55.]) + regions.append([1, 3, 5000., 6900., 10., 20., 45., 55.]) + regions.append([1, 2, 6900., 9000., 10., 20., 47., 52.]) + + # Force refinement near the island as the wave approaches: + + (xisland,yisland) = latlong(1600.e3, theta_island, 40., Rearth) + x1 = xisland - 1. + x2 = xisland + 1. + y1 = yisland - 1. + y2 = yisland + 1. + regions.append([4, 4, 7000., 1.e10, x1,x2,y1,y2]) + + x1 = xisland - 0.2 + x2 = xisland + 0.2 + y1 = yisland - 0.2 + y2 = yisland + 0.2 + regions.append([4, 5, 8000., 1.e10, x1,x2,y1,y2]) + + + + # ----- For developers ----- + # Toggle debugging print statements: + amrdata.dprint = False # print domain flags + amrdata.eprint = False # print err est flags + amrdata.edebug = False # even more err est flags + amrdata.gprint = False # grid bisection/clustering + amrdata.nprint = False # proper nesting output + amrdata.pprint = False # proj. of tagged points + amrdata.rprint = False # print regridding summary + amrdata.sprint = False # space/memory output + amrdata.tprint = False # time step reporting each level + amrdata.uprint = False # update/upbnd reporting + + return rundata + + # end of function setrun + # ---------------------- + + +#------------------- +def setgeo(rundata): +#------------------- + """ + Set GeoClaw specific runtime parameters. + """ + + try: + geo_data = rundata.geo_data + except: + print("*** Error, this rundata has no geo_data attribute") + raise AttributeError("Missing geo_data attribute") + + # == Physics == + geo_data.gravity = 9.81 + geo_data.coordinate_system = 2 + geo_data.earth_radius = 6367500.0 + + # == Forcing Options + geo_data.coriolis_forcing = False + + # == Algorithm and Initial Conditions == + geo_data.sea_level = 0.0 + geo_data.dry_tolerance = 0.001 + geo_data.friction_forcing = True + geo_data.manning_coefficient = 0.025 + geo_data.friction_depth = 1000000.0 + + # Refinement settings + refinement_data = rundata.refinement_data + refinement_data.variable_dt_refinement_ratios = False + refinement_data.wave_tolerance = 0.01 + refinement_data.deep_depth = 100.0 + refinement_data.max_level_deep = 3 + + # == settopo.data values == + rundata.topo_data.topofiles = [[2, 1, 1, 0.0, 10000000000.0, 'ocean.topotype2'], [2, 3, 4, 7.0, 10000000000.0, 'island.topotype2']] + topofiles = rundata.topo_data.topofiles + # for topography, append lines of the form + # [topotype, minlevel, maxlevel, t1, t2, fname] + + # == setdtopo.data values == + rundata.dtopo_data.dtopofiles = [] + dtopofiles = rundata.dtopo_data.dtopofiles + # for moving topography, append lines of the form : + # [topotype, minlevel,maxlevel,fname] + + + # == setqinit.data values == + rundata.qinit_data.qinit_type = 4 + rundata.qinit_data.qinitfiles = [[1, 2, 'hump.xyz']] + qinitfiles = rundata.qinit_data.qinitfiles + # for qinit perturbations, append lines of the form: (<= 1 allowed for now!) + # [minlev, maxlev, fname] + + # == fixedgrids.data values == + rundata.fixed_grid_data.fixedgrids = [] + fixedgrids = rundata.fixed_grid_data.fixedgrids + # for fixed grids append lines of the form + # [t1,t2,noutput,x1,x2,y1,y2,xpoints,ypoints,\ + # ioutarrivaltimes,ioutsurfacemax] + + # == fgmax.data values == + fgmax_files = rundata.fgmax_data.fgmax_files + # for fixed grids append to this list names of any fgmax input files + + + return rundata + # end of function setgeo + # ---------------------- + +if __name__ == '__main__': + # Set up run-time parameters and write all data files. + import sys + rundata = setrun(*sys.argv[1:]) + rundata.write() + diff --git a/tsunami/radial-ocean/awr2011_5.x/setrun_260.py b/tsunami/radial-ocean/awr2011_5.x/setrun_260.py new file mode 100644 index 00000000..418f2e62 --- /dev/null +++ b/tsunami/radial-ocean/awr2011_5.x/setrun_260.py @@ -0,0 +1,454 @@ +""" +Module to set up run time parameters for Clawpack -- AMRClaw code. + +The values set in the function setrun are then written out to data files +that will be read in by the Fortran code. + +""" + +from __future__ import print_function +import os +import numpy as np + +from clawpack.geoclaw.data import Rearth # radius of earth +from mapper import latlong + + +#------------------------------ +def setrun(claw_pkg='geoclaw'): +#------------------------------ + + """ + Define the parameters used for running Clawpack. + + INPUT: + claw_pkg expected to be "geoclaw" for this setrun. + + OUTPUT: + rundata - object of class ClawRunData + + """ + + from clawpack.clawutil import data + + + assert claw_pkg.lower() == 'geoclaw', "Expected claw_pkg = 'geoclaw'" + + num_dim = 2 + rundata = data.ClawRunData(claw_pkg, num_dim) + + + #------------------------------------------------------------------ + # Problem-specific parameters to be written to setprob.data: + #------------------------------------------------------------------ + theta_island = 260. + probdata = rundata.new_UserData(name='probdata',fname='setprob.data') + probdata.add_param('theta_island', theta_island, 'angle to island') + + + #------------------------------------------------------------------ + # GeoClaw specific parameters: + #------------------------------------------------------------------ + rundata = setgeo(rundata) + + #------------------------------------------------------------------ + # Standard Clawpack parameters to be written to claw.data: + #------------------------------------------------------------------ + + clawdata = rundata.clawdata # initialized when rundata instantiated + + + # Set single grid parameters first. + # See below for AMR parameters. + + + # --------------- + # Spatial domain: + # --------------- + + # Number of space dimensions: + clawdata.num_dim = num_dim + + # Lower and upper edge of computational domain: + clawdata.lower[0] = -20.0 # xlower + clawdata.upper[0] = 20.0 # xupper + clawdata.lower[1] = 20.0 # ylower + clawdata.upper[1] = 60.0 # yupper + + # Number of grid cells: + clawdata.num_cells[0] = 40 # mx + clawdata.num_cells[1] = 40 # my + + + # --------------- + # Size of system: + # --------------- + + # Number of equations in the system: + clawdata.num_eqn = 3 + + # Number of auxiliary variables in the aux array (initialized in setaux) + clawdata.num_aux = 3 + + # Index of aux array corresponding to capacity function, if there is one: + clawdata.capa_index = 2 + + + # ------------- + # Initial time: + # ------------- + + clawdata.t0 = 0.0 + + + # Restart from checkpoint file of a previous run? + # Note: If restarting, you must also change the Makefile to set: + # RESTART = True + # If restarting, t0 above should be from original run, and the + # restart_file 'fort.chkNNNNN' specified below should be in + # the OUTDIR indicated in Makefile. + + clawdata.restart = False # True to restart from prior results + clawdata.restart_file = 'fort.chk00006' # File to use for restart data + + + # ------------- + # Output times: + #-------------- + + # Specify at what times the results should be written to fort.q files. + # Note that the time integration stops after the final output time. + + clawdata.output_style = 1 + + if clawdata.output_style==1: + # Output ntimes frames at equally spaced times up to tfinal: + # Can specify num_output_times = 0 for no output + clawdata.num_output_times = 14 + clawdata.tfinal = 14000.0 + clawdata.output_t0 = True # output at initial (or restart) time? + + elif clawdata.output_style == 2: + # Specify a list or numpy array of output times: + # Include t0 if you want output at the initial time. + clawdata.output_times = [0., 0.1] + + elif clawdata.output_style == 3: + # Output every step_interval timesteps over total_steps timesteps: + clawdata.output_step_interval = 2 + clawdata.total_steps = 4 + clawdata.output_t0 = True # output at initial (or restart) time? + + + clawdata.output_format = 'binary' # 'ascii', 'binary' + + clawdata.output_q_components = 'all' # could be list such as [True,True] + clawdata.output_aux_components = 'none' # could be list + clawdata.output_aux_onlyonce = True # output aux arrays only at t0 + + + # --------------------------------------------------- + # Verbosity of messages to screen during integration: + # --------------------------------------------------- + + # The current t, dt, and cfl will be printed every time step + # at AMR levels <= verbosity. Set verbosity = 0 for no printing. + # (E.g. verbosity == 2 means print only on levels 1 and 2.) + clawdata.verbosity = 1 + + + + # -------------- + # Time stepping: + # -------------- + + # if dt_variable==True: variable time steps used based on cfl_desired, + # if dt_variable==Falseixed time steps dt = dt_initial always used. + clawdata.dt_variable = True + + # Initial time step for variable dt. + # (If dt_variable==0 then dt=dt_initial for all steps) + clawdata.dt_initial = 16.0 + + # Max time step to be allowed if variable dt used: + clawdata.dt_max = 1e+99 + + # Desired Courant number if variable dt used + clawdata.cfl_desired = 0.75 + # max Courant number to allow without retaking step with a smaller dt: + clawdata.cfl_max = 1.0 + + # Maximum number of time steps to allow between output times: + clawdata.steps_max = 5000 + + + # ------------------ + # Method to be used: + # ------------------ + + # Order of accuracy: 1 => Godunov, 2 => Lax-Wendroff plus limiters + clawdata.order = 2 + + # Use dimensional splitting? (not yet available for AMR) + clawdata.dimensional_split = 'unsplit' + + # For unsplit method, transverse_waves can be + # 0 or 'none' ==> donor cell (only normal solver used) + # 1 or 'increment' ==> corner transport of waves + # 2 or 'all' ==> corner transport of 2nd order corrections too + clawdata.transverse_waves = 2 + + + # Number of waves in the Riemann solution: + clawdata.num_waves = 3 + + # List of limiters to use for each wave family: + # Required: len(limiter) == num_waves + # Some options: + # 0 or 'none' ==> no limiter (Lax-Wendroff) + # 1 or 'minmod' ==> minmod + # 2 or 'superbee' ==> superbee + # 3 or 'vanleer' ==> van Leer + # 4 or 'mc' ==> MC limiter + clawdata.limiter = ['vanleer', 'vanleer', 'vanleer'] + + clawdata.use_fwaves = True # True ==> use f-wave version of algorithms + + # Source terms splitting: + # src_split == 0 or 'none' ==> no source term (src routine never called) + # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used, + # src_split == 2 or 'strang' ==> Strang (2nd order) splitting used, not recommended. + clawdata.source_split = 1 + + + # -------------------- + # Boundary conditions: + # -------------------- + + # Number of ghost cells (usually 2) + clawdata.num_ghost = 2 + + # Choice of BCs at xlower and xupper: + # 0 or 'user' => user specified (must modify bcNamr.f to use this option) + # 1 or 'extrap' => extrapolation (non-reflecting outflow) + # 2 or 'periodic' => periodic (must specify this at both boundaries) + # 3 or 'wall' => solid wall for systems where q(2) is normal velocity + + clawdata.bc_lower[0] = 'extrap' # at xlower + clawdata.bc_upper[0] = 'extrap' # at xupper + + clawdata.bc_lower[1] = 'extrap' # at ylower + clawdata.bc_upper[1] = 'extrap' # at yupper + + + # --------------- + # Gauges: + # --------------- + gauges = rundata.gaugedata.gauges + # for gauges append lines of the form [gaugeno, x, y, t1, t2] + + gaugeno = 0 + for d in [1570e3, 1590e3, 1610e3, 1630e3]: + gaugeno = gaugeno+1 + x,y = latlong(d, theta_island, 40., Rearth) + gauges.append([gaugeno, x, y, 0., 1e10]) + + + # -------------- + # Checkpointing: + # -------------- + + # Specify when checkpoint files should be created that can be + # used to restart a computation. + + clawdata.checkpt_style = 0 + + if clawdata.checkpt_style == 0: + # Do not checkpoint at all + pass + + elif clawdata.checkpt_style == 1: + # Checkpoint only at tfinal. + pass + + elif clawdata.checkpt_style == 2: + # Specify a list of checkpoint times. + clawdata.checkpt_times = [0.1,0.15] + + elif clawdata.checkpt_style == 3: + # Checkpoint every checkpt_interval timesteps (on Level 1) + # and at the final time. + clawdata.checkpt_interval = 5 + + + + # --------------- + # AMR parameters: (written to amr.data) + # --------------- + amrdata = rundata.amrdata + + # max number of refinement levels: + amrdata.amr_levels_max = 5 + + # List of refinement ratios at each level (length at least amr_level_max-1) + amrdata.refinement_ratios_x = [4, 4, 4, 4] + amrdata.refinement_ratios_y = [4, 4, 4, 4] + amrdata.refinement_ratios_t = [4, 4, 4, 1] + + + # Specify type of each aux variable in amrdata.auxtype. + # This must be a list of length num_aux, each element of which is one of: + # 'center', 'capacity', 'xleft', or 'yleft' (see documentation). + amrdata.aux_type = ['center', 'capacity', 'yleft'] + + + # Flag for refinement based on Richardson error estimater: + amrdata.flag_richardson = False # use Richardson? + amrdata.flag_richardson_tol = 1.0 # Richardson tolerance + + # Flag for refinement using routine flag2refine: + amrdata.flag2refine = True # use this? + amrdata.flag2refine_tol = 0.5 # tolerance used in this routine + # Note: in geoclaw the refinement tolerance is set as wave_tolerance below + # and flag2refine_tol is unused! + + # steps to take on each level L between regriddings of level L+1: + amrdata.regrid_interval = 3 + + # width of buffer zone around flagged points: + # (typically the same as regrid_interval so waves don't escape): + amrdata.regrid_buffer_width = 2 + + # clustering alg. cutoff for (# flagged pts) / (total # of cells refined) + # (closer to 1.0 => more small grids may be needed to cover flagged cells) + amrdata.clustering_cutoff = 0.7 + + # print info about each regridding up to this level: + amrdata.verbosity_regrid = 0 + + + # --------------- + # Regions: + # --------------- + regions = rundata.regiondata.regions + # to specify regions of refinement append lines of the form + # [minlevel,maxlevel,t1,t2,x1,x2,y1,y2] + + regions.append([1, 3, 0., 5000., -5., 20., 35., 45.]) + regions.append([1, 2, 5000., 6900., 10., 20., 35., 45.]) + regions.append([1, 3, 5000., 6900., 12., 20., 39., 43.]) + + + # Force refinement near the island as the wave approaches: + + (xisland,yisland) = latlong(1600.e3, theta_island, 40., Rearth) + x1 = xisland - 1. + x2 = xisland + 1. + y1 = yisland - 1. + y2 = yisland + 1. + regions.append([4, 4, 7000., 1.e10, x1,x2,y1,y2]) + + x1 = xisland - 0.2 + x2 = xisland + 0.2 + y1 = yisland - 0.2 + y2 = yisland + 0.2 + regions.append([4, 5, 8000., 1.e10, x1,x2,y1,y2]) + + + + # ----- For developers ----- + # Toggle debugging print statements: + amrdata.dprint = False # print domain flags + amrdata.eprint = False # print err est flags + amrdata.edebug = False # even more err est flags + amrdata.gprint = False # grid bisection/clustering + amrdata.nprint = False # proper nesting output + amrdata.pprint = False # proj. of tagged points + amrdata.rprint = False # print regridding summary + amrdata.sprint = False # space/memory output + amrdata.tprint = False # time step reporting each level + amrdata.uprint = False # update/upbnd reporting + + return rundata + + # end of function setrun + # ---------------------- + + +#------------------- +def setgeo(rundata): +#------------------- + """ + Set GeoClaw specific runtime parameters. + """ + + try: + geo_data = rundata.geo_data + except: + print("*** Error, this rundata has no geo_data attribute") + raise AttributeError("Missing geo_data attribute") + + # == Physics == + geo_data.gravity = 9.81 + geo_data.coordinate_system = 2 + geo_data.earth_radius = 6367500.0 + + # == Forcing Options + geo_data.coriolis_forcing = False + + # == Algorithm and Initial Conditions == + geo_data.sea_level = 0.0 + geo_data.dry_tolerance = 0.001 + geo_data.friction_forcing = True + geo_data.manning_coefficient = 0.025 + geo_data.friction_depth = 1000000.0 + + # Refinement settings + refinement_data = rundata.refinement_data + refinement_data.variable_dt_refinement_ratios = False + refinement_data.wave_tolerance = 0.01 + refinement_data.deep_depth = 100.0 + refinement_data.max_level_deep = 3 + + # == settopo.data values == + rundata.topo_data.topofiles = [[2, 1, 1, 0.0, 10000000000.0, 'ocean.topotype2'], [2, 3, 4, 7.0, 10000000000.0, 'island.topotype2']] + topofiles = rundata.topo_data.topofiles + # for topography, append lines of the form + # [topotype, minlevel, maxlevel, t1, t2, fname] + + # == setdtopo.data values == + rundata.dtopo_data.dtopofiles = [] + dtopofiles = rundata.dtopo_data.dtopofiles + # for moving topography, append lines of the form : + # [topotype, minlevel,maxlevel,fname] + + + # == setqinit.data values == + rundata.qinit_data.qinit_type = 4 + rundata.qinit_data.qinitfiles = [[1, 2, 'hump.xyz']] + qinitfiles = rundata.qinit_data.qinitfiles + # for qinit perturbations, append lines of the form: (<= 1 allowed for now!) + # [minlev, maxlev, fname] + + # == fixedgrids.data values == + rundata.fixed_grid_data.fixedgrids = [] + fixedgrids = rundata.fixed_grid_data.fixedgrids + # for fixed grids append lines of the form + # [t1,t2,noutput,x1,x2,y1,y2,xpoints,ypoints,\ + # ioutarrivaltimes,ioutsurfacemax] + + # == fgmax.data values == + fgmax_files = rundata.fgmax_data.fgmax_files + # for fixed grids append to this list names of any fgmax input files + + + return rundata + # end of function setgeo + # ---------------------- + +if __name__ == '__main__': + # Set up run-time parameters and write all data files. + import sys + rundata = setrun(*sys.argv[1:]) + rundata.write() +