diff --git a/illapel2015_chile/Makefile b/illapel2015_chile/Makefile new file mode 100644 index 0000000..2106f06 --- /dev/null +++ b/illapel2015_chile/Makefile @@ -0,0 +1,68 @@ +# 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/illapel2015_chile/README.md b/illapel2015_chile/README.md new file mode 100644 index 0000000..b0f837b --- /dev/null +++ b/illapel2015_chile/README.md @@ -0,0 +1,80 @@ +# 2015 Illapel Earthquake Report + +This example contains code for running a simulation of the tsunami generated by the 2015 Illapel earthquake on the Chilean coast. + +## The Earthquake + +*Sources: USGS (https://earthquake.usgs.gov/earthquakes/eventpage/us20003k7a/executive), NOAA (https://www.ngdc.noaa.gov/hazard/16sep2015.html)* + +On September 16, 2015, an earthquake with a magnitude of 8.3 occurred on the boundary between the Nazca and South America tectonic plates, north of the rupture zone from the 2010 Chile earthquake. It produced significant aftershocks, some of which exceeded magnitude 7.0, and it was reported to have caused 15 deaths and 14 injuries. + +In this example, we model the tsunami resulting from the earthquake. We validate our model using data from two tide gauges in the region. + +## Topography/Bathymetry Data + +*Topography data was sourced from the ETOPO 10-arc-minute model (Esri ASCII format) here: http://depts.washington.edu/clawpack/geoclaw/topo/etopo/etopo10min120W60W60S0S.asc* + +In the example, [maketopo.py](maketopo.py) downloads the data directly from source into the scratch directory within GeoClaw. + +```python +# maketopo.py, lines 26-27 +url = 'http://depts.washington.edu/clawpack/geoclaw/topo/etopo/' + topo_fname +clawpack.clawutil.data.get_remote_file(url, output_dir=scratch_dir, file_name=topo_fname, verbose=True) +``` + +If you want to run a simulation with more refined topography, simply download the file into the scratch directory (or edit the download URL) and edit the appropriate filenames within [maketopo.py](maketopo.py) and [setrun.py](setrun.py). If the geographical bounds or topography file format are different, note that these parameters must be changed in [setrun.py](setrun.py). + +## Moving Topography Data + +*Moving topography (dtopo) data was sourced from the USGS finite fault model for the earthquake here: https://earthquake.usgs.gov/product/finite-fault/us20003k7a/us/1539809967421/basic_inversion.param* + +As with topography data, moving topography (dtopo) data is directly downloaded. + +```python +# maketopo.py, lines 57-59 +param_fname = 'basic_inversion.param' +csv_fname = 'chile2015.csv' +get_csv('https://earthquake.usgs.gov/product/finite-fault/us20003k7a/us/1539809967421/' + param_fname, param_fname, csv_fname) +``` +In order to produce a CSV file of the subfaults readable by the `clawpack.geoclaw.dtopotools` module, the data file is reformatted in [get_dtopo_csv.py](get_dtopo_csv.py). Therefore, to use different finite fault data, edit this file alongside the download URL in [maketopo.py](maketopo.py) (and potentially the file format in [setrun.py](setrun.py)). + +## GeoClaw Results + +The simulation was set to run 5 hours from the origin time of the earthquake, with a maximum of 6 adaptive mesh refinement layers. Areas of high refinement were manually implemented around tide gauges. Surface plots and tide gauge plots are generated by [setplot.py](setplot.py), and they show predicted tide residuals ranging from -0.075 to +0.125. Overall, these values match NOAA tide data relatively well. + +## Validation Results + +*Tide data was downloaded directly from the NOAA database for the earthquake here: https://www.ngdc.noaa.gov/hazard/dart/2015chile.html* + +We test our simulation against tide residual values at two DART stations: gauge no. 32412 (SW of Lima, Peru) and gauge no. 32402 (W of Caldera, Chile). The simulation matched the residual patterns at the closer gauge (32402) relatively well but those at farther gauge (32412) relatively poorly. + +For additional testing, one can add more gauges for plotting. To do so, the download must be specified in [maketopo.py](maketopo.py) and the plotting in [setplot.py](setplot.py). + +## Conclusion + +Topography and moving topography data for the 2015 Illapel earthquake were obtained from the internet and used by GeoClaw to run a simulation of the regional sea surface. Tide residuals were compared with observed NOAA data and demonstrated good local predictiveness but poor faraway predictiveness. Better topography data and adaptive mesh refinement customization may yield improved results. More tide gauges may provide a better picture of the simulation's predictive scope. + +## Installation + +Make sure to have the following dependencies downloaded in order to run this example: +- [Clawpack](https://www.clawpack.org/installing.html) +- [pandas](https://pandas.pydata.org/) + +## Usage + +Make sure you have GNU Make installed, and run the following terminal commands inside of the `illapel2015_chile` directory: + +``` +make clobber +make .data +make xgeoclaw +make .output +make .plots +``` + +To optimize simulation speed, consider [setting OpenMP compiler flags for multiprocessing](https://www.clawpack.org/openmp.html). + +--- + +For questions, please contact Ashwin Padaki at ap4025@columbia.edu. + diff --git a/illapel2015_chile/format_gauge.py b/illapel2015_chile/format_gauge.py new file mode 100644 index 0000000..ecaa770 --- /dev/null +++ b/illapel2015_chile/format_gauge.py @@ -0,0 +1,14 @@ +"""Script for formatting NOAA gauge data for the Chile 2015 earthquake""" +def format_gauges_chile2015(gaugeinfo): + import pandas as pd + for gaugeno, offset in gaugeinfo: + df = pd.read_csv('dart{}.txt'.format(gaugeno), delim_whitespace=True) + df.columns = ['time', 'year', 'month', 'day', 'hours', 'minutes', 'seconds', 'raw_obs', 'fitted_tidal', 'residual', 'blank'] + # only store time and residual data + df = df[['time', 'residual']] + # earthquake origin time (from USGS): 2015-09-16 22:54:32.860 (UTC) + origin_time = 259 + (22 * 60**2 + 54*60 + 32.860)/(24 * 60**2) + df.time = df.time.map(lambda time : 24*60**2 * (time-origin_time)) + # option to include manual offsets to correct for equilibria that are not at sea level + df.residual = df.residual.map(lambda res : res - offset) + df.to_csv('{}_notide.txt'.format(gaugeno), index=False, header=False, sep=' ') \ No newline at end of file diff --git a/illapel2015_chile/get_dtopo_csv.py b/illapel2015_chile/get_dtopo_csv.py new file mode 100644 index 0000000..53d59f2 --- /dev/null +++ b/illapel2015_chile/get_dtopo_csv.py @@ -0,0 +1,22 @@ +"""Script for formatting USGS finite fault data into a CSV file readable by GeoClaw's dtopotools module""" +def get_csv(url, param_fname, csv_fname): + import clawpack.clawutil.data + import os + try: + CLAW = os.environ['CLAW'] + except: + raise Exception("*** Must first set CLAW enviornment variable") + scratch_dir = os.path.join(CLAW, 'geoclaw', 'scratch') + clawpack.clawutil.data.get_remote_file(url, output_dir=scratch_dir, file_name=param_fname, verbose=True) + + import pandas as pd + df = pd.read_csv(os.path.join(scratch_dir,param_fname), skiprows=9, delim_whitespace=True) + df.columns = ['latitude', 'longitude', 'depth', 'slip', 'rake', 'strike', 'dip', 'rupture_time', 'rise_time', 'fall_time', 'mo'] + # length/width can be automated, but for convenience, they're manually written here + df['length'] = 12.00 + df['width'] = 8.80 + # rigidity = (moment seismicity) / (slip * area) + # cm * km * km = 10^4 m^3 + df['mu'] = df['mo'] / (10**4 * df['length'] * df['width'] * df['slip']) + df.drop(['fall_time', 'mo'], axis=1, inplace=True) + df.to_csv(csv_fname, index=False) \ No newline at end of file diff --git a/illapel2015_chile/maketopo.py b/illapel2015_chile/maketopo.py new file mode 100644 index 0000000..5efed5d --- /dev/null +++ b/illapel2015_chile/maketopo.py @@ -0,0 +1,153 @@ +""" +Script for downloading gauge and topography data, and generating differential topography data +Call functions with makeplots==True to create plots of topo, slip, and dtopo. +""" + +from __future__ import absolute_import +from __future__ import print_function +import os + +import clawpack.clawutil.data + +try: + CLAW = os.environ['CLAW'] +except: + raise Exception("*** Must first set CLAW enviornment variable") + +scratch_dir = os.path.join(CLAW, 'geoclaw', 'scratch') + +def get_topo(makeplots=False): + """ + Retrieve the topo file from the GeoClaw repository. + """ + from clawpack.geoclaw import topotools + # download 10 arc-minute topography data in region around Chile + topo_fname = 'etopo10min120W60W60S0S.asc' + url = 'http://depts.washington.edu/clawpack/geoclaw/topo/etopo/' + topo_fname + clawpack.clawutil.data.get_remote_file(url, output_dir=scratch_dir, file_name=topo_fname, verbose=True) + + if makeplots: + from matplotlib import pyplot as plt + topo = topotools.Topography(os.path.join(scratch_dir,topo_fname), topo_type=2) + topo.plot() + fname = os.path.splitext(topo_fname)[0] + '.png' + plt.savefig(fname) + print("Created ",fname) + +# generate side-by-side fault slip and seafloor deformation contour plots +def plot_subfaults_dZ(t, fig, fault, dtopo, xlower, xupper, ylower, yupper, xylim, dz_max): + fig.clf() + ax1 = fig.add_subplot(121) + ax2 = fig.add_subplot(122) + fault.plot_subfaults(axes=ax1, slip_color=True, + slip_time=t, xylim=xylim) + dtopo.plot_dZ_colors(axes=ax2, t=t, cmax_dZ=dz_max) + ax2.set_xlim(xlower,xupper) + ax2.set_ylim(ylower,yupper) + return fig + +# create dtopo file from finite fault data +def make_dtopo(makeplots=False): + from clawpack.geoclaw import dtopotools + from matplotlib import pyplot as plt + import numpy + from get_dtopo_csv import get_csv + + # specify the URL for the USGS finite fault data file and properly format into a CSV file + param_fname = 'basic_inversion.param' + csv_fname = 'chile2015.csv' + get_csv('https://earthquake.usgs.gov/product/finite-fault/us20003k7a/us/1539809967421/' + param_fname, param_fname, csv_fname) + + dtopo_fname = os.path.join(scratch_dir, "dtopo_usgs150916.tt3") + + # manually specify the bounds for the fault model (leave unchanged) + xlower = -73 + xupper = -70.5 + ylower = -33 + yupper = -29.5 + xylim = [xlower,xupper,ylower,yupper] + + input_units = { + 'depth' : 'km', + 'slip' : 'cm', + 'length' : 'km', + 'width' : 'km', + 'mu' : 'Pa' + } + + fault = dtopotools.CSVFault() + fault.read(csv_fname,input_units=input_units, coordinate_specification='noaa sift') + fault.rupture_type = 'dynamic' + + print ("%s subfaults read in " % len(fault.subfaults)) + + if os.path.exists(dtopo_fname): + print("*** Not regenerating dtopo file (already exists): %s" \ + % dtopo_fname) + else: + print("Using Okada model to create dtopo file") + points_per_degree = 60 + mx = int((xupper - xlower)*points_per_degree + 1) + my = int((yupper - ylower)*points_per_degree + 1) + x = numpy.linspace(xlower, xupper, mx) + y = numpy.linspace(ylower, yupper, my) + + # list of the times (in seconds) when the dtopo data is captured + # here, `times` consists of the interval [0,200] divided evenly into 20 markers (0 and 200 included) + # note: ensure that the maximum value in `times` exceeds the maximum rupture time of any subfault + times = numpy.linspace(0,200,20) + + # create and write dtopo file + fault.create_dtopography(x,y,times) + dtopo = fault.dtopo + dtopo.write(dtopo_fname, dtopo_type=3) + + import clawpack.visclaw.JSAnimation.JSAnimation_frametools as J + plotdir = '_fault_slip_plots' + J.make_plotdir(plotdir, clobber=True) + fig = plt.figure(figsize=(12,5)) + dzmax = abs(dtopo.dZ).max() + + # generate fault slip and deformation plots + for k,t in enumerate(times): + plot_subfaults_dZ(t,fig, fault, dtopo, xlower, xupper, ylower, yupper, xylim, dzmax) + J.save_frame(k, plotdir = plotdir, verbose=True) + + if makeplots: + if fault.dtopo is None: + # read in the pre-existing file: + print("Reading in dtopo file...") + dtopo = dtopotools.DTopography() + dtopo.read(dtopo_fname, dtopo_type=3) + x = dtopo.x + y = dtopo.y + plt.figure(figsize=(12,7)) + ax1 = plt.subplot(121) + ax2 = plt.subplot(122) + fault.plot_subfaults(axes=ax1,slip_color=True) + ax1.set_xlim(x.min(),x.max()) + ax1.set_ylim(y.min(),y.max()) + dtopo.plot_dZ_colors(1.,axes=ax2) + fname = os.path.splitext(os.path.split(dtopo_fname)[-1])[0] + '.png' + plt.savefig(fname) + print("Created ",fname) + +# grabs gauge data from NOAA and formats it +def get_gauge(): + # find more gauge numbers here: https://www.ngdc.noaa.gov/hazard/dart/2015chile.html + gaugeinfo = [] + # gaugeinfo.append([gauge_no, offset]) where `offset` is the nonzero equilibrium in the data + gaugeinfo.append([32412, -0.03]) + gaugeinfo.append([32402, 0.00]) + + for gaugeno, offset in gaugeinfo: + url = 'https://www.ngdc.noaa.gov/hazard/data/DART/20150916_chile/dart{}_20150915to20150921_meter.txt'.format(gaugeno) + clawpack.clawutil.data.get_remote_file(url, output_dir='.', file_name='dart{}.txt'.format(gaugeno), verbose=True) + + from format_gauge import format_gauges_chile2015 + format_gauges_chile2015(gaugeinfo) + +if __name__=='__main__': + get_topo(False) + make_dtopo(False) + get_gauge() \ No newline at end of file diff --git a/illapel2015_chile/setplot.py b/illapel2015_chile/setplot.py new file mode 100644 index 0000000..c530863 --- /dev/null +++ b/illapel2015_chile/setplot.py @@ -0,0 +1,212 @@ + +""" +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 __future__ import absolute_import +from __future__ import print_function +import numpy as np + +from six.moves import range + +try: + TG32412 = np.loadtxt('32412_notide.txt') + TG32402 = np.loadtxt('32402_notide.txt') + # add a NumPy reader for each tide gauge downloaded in maketopo.py +except: + print("*** Could not load DART data file") + +#-------------------------- +def setplot(plotdata=None): +#-------------------------- + + """ + Specify what is to be plotted at each frame. + Input: plotdata, an instance of pyclaw.plotters.data.ClawPlotData. + Output: a modified version of plotdata. + + """ + + from clawpack.visclaw import colormaps, geoplot + from numpy import linspace + + if plotdata is None: + from clawpack.visclaw.data import ClawPlotData + plotdata = ClawPlotData() + + + plotdata.clearfigures() # clear any old figures,axes,items data + plotdata.format = 'ascii' # 'ascii' or 'binary' to match setrun.py + + + # 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 surface + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Surface', figno=0) + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.title = 'Surface' + plotaxes.scaled = True + + def fixup(current_data): + import pylab + addgauges(current_data) + t = current_data.t + t = t / 3600. # hours + pylab.title('Surface at %4.2f hours' % t, fontsize=20) + pylab.xticks(fontsize=15) + pylab.yticks(fontsize=15) + plotaxes.afteraxes = fixup + + # 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 = -0.05 + plotitem.pcolor_cmax = 0.05 + 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 = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 100.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [1,1,0] + plotitem.patchedges_show = 1 + plotaxes.xlimits = [-120,-60] + plotaxes.ylimits = [-60,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(-3000,-3000,1) + plotitem.amr_contour_colors = ['y'] # color on each level + plotitem.kwargs = {'linestyles':'solid','linewidths':2} + plotitem.amr_contour_show = [1,0,0] + plotitem.celledges_show = 0 + plotitem.patchedges_show = 0 + + + #----------------------------------------- + # Figures for gauges + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Surface at gauges', figno=300, \ + type='each_gauge') + plotfigure.clf_each_gauge = True + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = 'auto' + plotaxes.ylimits = 'auto' + 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') + plotitem.show = False + + 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, xticks, floor, axis, xlabel + t = current_data.t + gaugeno = current_data.gaugeno + + # add conditional blocks of this form for each tide gauge, with the desired axis bounds + if gaugeno == 32412: + try: + plot(TG32412[:,0], TG32412[:,1], 'r') + legend(['GeoClaw','Obs'],loc='lower right') + except: pass + axis((0,t.max(),-0.1,0.1)) + + if gaugeno == 32402: + try: + plot(TG32402[:,0], TG32402[:,1], 'r') + legend(['GeoClaw','Obs'],loc='lower right') + except: pass + axis((0,t.max(),-0.2,0.2)) + + plot(t, 0*t, 'k') + n = int(floor(t.max()/3600.) + 2) + xticks([3600*i for i in range(n)], ['%i' % i for i in range(n)]) + xlabel('time (hours)') + + plotaxes.afteraxes = add_zeroline + + + #----------------------------------------- + # Plots of timing (CPU and wall time): + + def make_timing_plots(plotdata): + from clawpack.visclaw import plot_timing_stats + import os,sys + try: + timing_plotdir = plotdata.plotdir + '/_timing_figures' + os.system('mkdir -p %s' % timing_plotdir) + # adjust units for plots based on problem: + units = {'comptime':'seconds', 'simtime':'hours', + 'cell':'millions'} + plot_timing_stats.make_plots(outdir=plotdata.outdir, + make_pngs=True, + plotdir=timing_plotdir, + units=units) + except: + print('*** Error making timing plots') + + otherfigure = plotdata.new_otherfigure(name='timing plots', + fname='_timing_figures/timing.html') + otherfigure.makefig = make_timing_plots + + + #----------------------------------------- + + # Parameters used only when creating html and/or latex hardcopy + # e.g., via pyclaw.plotters.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? + plotdata.parallel = True # make multiple frame png's at once + + return plotdata \ No newline at end of file diff --git a/illapel2015_chile/setrun.py b/illapel2015_chile/setrun.py new file mode 100644 index 0000000..d39c785 --- /dev/null +++ b/illapel2015_chile/setrun.py @@ -0,0 +1,454 @@ +""" +Module to set up run time parameters for Clawpack. +Left relatively unchanged from file in chile2010 example, outside of filenames and user parameters + +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 absolute_import +from __future__ import print_function +import os +import numpy as np + +try: + CLAW = os.environ['CLAW'] +except: + raise Exception("*** Must first set CLAW enviornment variable") + +# Scratch directory for storing topo and dtopo files: +scratch_dir = os.path.join(CLAW, 'geoclaw', 'scratch') + + +#------------------------------ +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: + #------------------------------------------------------------------ + + #probdata = rundata.new_UserData(name='probdata',fname='setprob.data') + + + #------------------------------------------------------------------ + # GeoClaw specific parameters: + #------------------------------------------------------------------ + rundata = setgeo(rundata) + + #------------------------------------------------------------------ + # Standard Clawpack parameters to be written to claw.data: + # (or to amr2ez.data for AMR) + #------------------------------------------------------------------ + 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] = -120.0 # west longitude + clawdata.upper[0] = -60.0 # east longitude + + clawdata.lower[1] = -60.0 # south latitude + clawdata.upper[1] = 0.0 # north latitude + + + + # Number of grid cells: Coarsest grid + clawdata.num_cells[0] = 30 + clawdata.num_cells[1] = 30 + + # --------------- + # 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? + # 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.chk00096' # 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. + # The solution at initial time t0 is always written in addition. + + clawdata.output_style = 1 + + if clawdata.output_style==1: + # Output nout frames at equally spaced times up to tfinal: + clawdata.num_output_times = 20 + clawdata.tfinal = 5*3600. + clawdata.output_t0 = True # output at initial (or restart) time? + + elif clawdata.output_style == 2: + # Specify a list of output times. + clawdata.output_times = [0.5, 1.0] + + elif clawdata.output_style == 3: + # Output every iout timesteps with a total of ntot time steps: + clawdata.output_step_interval = 1 + clawdata.total_steps = 3 + clawdata.output_t0 = True + + + clawdata.output_format = 'ascii' # 'ascii' or 'binary' + + clawdata.output_q_components = 'all' # need all + clawdata.output_aux_components = 'none' # eta=h+B is in q + clawdata.output_aux_onlyonce = False # output aux arrays each frame + + + + # --------------------------------------------------- + # 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==1: variable time steps used based on cfl_desired, + # if dt_variable==0: fixed time steps dt = dt_initial will always be 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 = 0.2 + + # Max time step to be allowed if variable dt used: + clawdata.dt_max = 1e+99 + + # Desired Courant number if variable dt used, and max to allow without + # retaking step with a smaller dt: + clawdata.cfl_desired = 0.75 + 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 'mc' ==> MC limiter + # 4 or 'vanleer' ==> van Leer + clawdata.limiter = ['mc', 'mc', 'mc'] + + 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 = 'godunov' + + + # -------------------- + # Boundary conditions: + # -------------------- + + # Number of ghost cells (usually 2) + clawdata.num_ghost = 2 + + # Choice of BCs at xlower and xupper: + # 0 => user specified (must modify bcN.f to use this option) + # 1 => extrapolation (non-reflecting outflow) + # 2 => periodic (must specify this at both boundaries) + # 3 => solid wall for systems where q(2) is normal velocity + + clawdata.bc_lower[0] = 'extrap' + clawdata.bc_upper[0] = 'extrap' + + clawdata.bc_lower[1] = 'extrap' + clawdata.bc_upper[1] = 'extrap' + + + + # -------------- + # 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 np.abs(clawdata.checkpt_style) == 1: + # Checkpoint only at tfinal. + pass + + elif np.abs(clawdata.checkpt_style) == 2: + # Specify a list of checkpoint times. + clawdata.checkpt_times = [0.1,0.15] + + elif np.abs(clawdata.checkpt_style) == 3: + # Checkpoint every checkpt_interval timesteps (on Level 1) + # and at the final time. + clawdata.checkpt_interval = 5 + + + # --------------- + # AMR parameters: + # --------------- + amrdata = rundata.amrdata + + # maximum size of patches in each direction (matters in parallel): + amrdata.max1d = 60 + + # max number of refinement levels: + amrdata.amr_levels_max = 6 + + # List of refinement ratios at each level (length at least mxnest-1) + # note: for the topography data used in this example, model accuracy plateaus at level 7 or 8 + amrdata.refinement_ratios_x = [2,5,5,5,5] + amrdata.refinement_ratios_y = [2,5,5,5,5] + amrdata.refinement_ratios_t = [2,5,5,5,5] + + # Specify type of each aux variable in amrdata.auxtype. + # This must be a list of length maux, each element of which is one of: + # 'center', 'capacity', 'xleft', or 'yleft' (see documentation). + + amrdata.aux_type = ['center','capacity','yleft'] + + + # Flag using refinement routine flag2refine rather than richardson error + amrdata.flag_richardson = False # use Richardson? + amrdata.flag_richardson_tol = 0.002 # Richardson tolerance + amrdata.flag2refine = True + + # 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.700000 + + # print info about each regridding up to this level: + amrdata.verbosity_regrid = 0 + + # ----- 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 = True # time step reporting each level + amrdata.uprint = False # update/upbnd reporting + + # More AMR parameters can be set -- see the defaults in pyclaw/data.py + + # --------------- + # Regions: + # --------------- + rundata.regiondata.regions = [] + # to specify regions of refinement append lines of the form + # # [minlevel,maxlevel,t1,t2,x1,x2,y1,y2] + # rundata.regiondata.regions.append([3, 3, 0., 10000., -77,-67,-40,-30]) + # rundata.regiondata.regions.append([3, 3, 8000., 26000., -90,-80,-30,-15]) + + # --------------- + # Gauges: + # --------------- + rundata.gaugedata.gauges = [] + # for gauges append lines of the form [gaugeno, x, y, t1, t2] + # note: make sure the gauge numbers here match with those in setplot.py and maketopo.py + rundata.gaugedata.gauges.append([32412, -86.392, -17.975, 0., 1.e10]) + rundata.gaugedata.gauges.append([32402, -73.983, -26.743, 0., 1.e10]) + + # Manually-specified regions around gauges for higher refinement + # AMR level 5-6 reserved for within 1 latitude/longitude around gauges + # MAR level 1-4 for the remaining region + rundata.regiondata.regions.append([1, 4, 0., 1e9, clawdata.lower[0], clawdata.upper[0], clawdata.lower[1], clawdata.upper[1]]) + from clawpack.amrclaw.data import FlagRegion + + for gaugeno,x1,x2,y1,y2 in [[32412,-86.392-1,-86.392+1,-17.975-1,-17.975+1], + [32402,-73.983-1,-73.983+1,-26.743-1,-26.743+1]]: + flagregion = FlagRegion(num_dim=2) + flagregion.name = 'gauge{}'.format(gaugeno) + flagregion.minlevel = 5 + flagregion.maxlevel = 6 + flagregion.t1 = 0. + flagregion.t2 = 1e10 + flagregion.spatial_region_type = 1 + flagregion.spatial_region = [x1,x2,y1,y2] + rundata.flagregiondata.flagregions.append(flagregion) + + + return rundata + # end of function setrun + # ---------------------- + + +#------------------- +def setgeo(rundata): +#------------------- + """ + Set GeoClaw specific runtime parameters. + For documentation see .... + """ + + 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 = 6367.5e3 + + # == Forcing Options + geo_data.coriolis_forcing = False + + # == Algorithm and Initial Conditions == + geo_data.sea_level = 0.0 + geo_data.dry_tolerance = 1.e-3 + geo_data.friction_forcing = True + geo_data.manning_coefficient =.025 + geo_data.friction_depth = 1e6 + + # Refinement settings + refinement_data = rundata.refinement_data + refinement_data.variable_dt_refinement_ratios = True + refinement_data.wave_tolerance = 1.e-1 + + # == settopo.data values == + topo_data = rundata.topo_data + # for topography, append lines of the form + # [topotype, fname] + topo_path = os.path.join(scratch_dir, 'etopo10min120W60W60S0S.asc') + + topo_data.topofiles.append([2, topo_path]) + + # == setdtopo.data values == + dtopo_data = rundata.dtopo_data + # for moving topography, append lines of the form : (<= 1 allowed for now!) + # [topotype, fname] + dtopo_path = os.path.join(scratch_dir, 'dtopo_usgs150916.tt3') + dtopo_data.dtopofiles.append([3,dtopo_path]) + dtopo_data.dt_max_dtopo = 0.2 + + + # == setqinit.data values == + rundata.qinit_data.qinit_type = 0 + rundata.qinit_data.qinitfiles = [] + # for qinit perturbations, append lines of the form: (<= 1 allowed for now!) + # [fname] + + # == setfixedgrids.data values == + fixed_grids = rundata.fixed_grid_data + # for fixed grids append lines of the form + # [t1,t2,noutput,x1,x2,y1,y2,xpoints,ypoints,\ + # ioutarrivaltimes,ioutsurfacemax] + + return rundata + # end of function setgeo + # ---------------------- + + + +if __name__ == '__main__': + # Set up run-time parameters and write all data files. + import sys + from clawpack.geoclaw import kmltools + + rundata = setrun(*sys.argv[1:]) + rundata.write() + + kmltools.make_input_data_kmls(rundata)