diff --git a/examples/storm-surge/synthetic/.gitignore b/examples/storm-surge/synthetic/.gitignore deleted file mode 100644 index 7c3cf68c4..000000000 --- a/examples/storm-surge/synthetic/.gitignore +++ /dev/null @@ -1,2 +0,0 @@ -fake.storm -topo.tt2 \ No newline at end of file diff --git a/examples/storm-surge/synthetic/Makefile b/examples/storm-surge/synthetic/Makefile deleted file mode 100644 index 110c8f7a3..000000000 --- a/examples/storm-surge/synthetic/Makefile +++ /dev/null @@ -1,59 +0,0 @@ - -# 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 -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: -# ---------------------------------------- - -RIEMANN = $(CLAW)/riemann/src - -MODULES = \ - -SOURCES = \ - $(RIEMANN)/rpn2_geoclaw.f \ - $(RIEMANN)/rpt2_geoclaw.f \ - $(RIEMANN)/geoclaw_riemann_utils.f \ - -#------------------------------------------------------------------- -# Include Makefile containing standard definitions and make options: -include $(CLAWMAKE) - -### DO NOT remove this line - make depends on it ### diff --git a/examples/storm-surge/synthetic/README.rst b/examples/storm-surge/synthetic/README.rst deleted file mode 100644 index dbdf650c8..000000000 --- a/examples/storm-surge/synthetic/README.rst +++ /dev/null @@ -1,13 +0,0 @@ - -.. _geoclaw_examples_storm_surge_square_basin: - -Basic Example of Storm Surge from Hurricane Ike -=============================================== - -This example contains the data and setup for running a storm surge simulation -for a square(ish) basin in cartesian coordinates. The example can be run via:: - - make all - -This is a good example to use if you are interested in idealized storm surge -as the storm is created in the *setrun.py* rather than downloaded. \ No newline at end of file diff --git a/examples/storm-surge/synthetic/setplot.py b/examples/storm-surge/synthetic/setplot.py deleted file mode 100644 index 89292549f..000000000 --- a/examples/storm-surge/synthetic/setplot.py +++ /dev/null @@ -1,230 +0,0 @@ - -from __future__ import absolute_import -from __future__ import print_function - -import os - -import numpy -import matplotlib.pyplot as plt -import datetime - -import clawpack.visclaw.colormaps as colormap -import clawpack.visclaw.gaugetools as gaugetools -import clawpack.clawutil.data as clawutil -import clawpack.amrclaw.data as amrclaw -import clawpack.geoclaw.data as geodata - - -import clawpack.geoclaw.surge.plot as surgeplot - -try: - from setplotfg import setplotfg -except: - setplotfg = None - - -def setplot(plotdata=None): - """""" - - if plotdata is None: - from clawpack.visclaw.data import ClawPlotData - plotdata = ClawPlotData() - - # clear any old figures,axes,items data - plotdata.clearfigures() - plotdata.format = 'binary' - - # Load data from output - clawdata = clawutil.ClawInputData(2) - clawdata.read(os.path.join(plotdata.outdir, 'claw.data')) - physics = geodata.GeoClawData() - physics.read(os.path.join(plotdata.outdir, 'geoclaw.data')) - surge_data = geodata.SurgeData() - surge_data.read(os.path.join(plotdata.outdir, 'surge.data')) - friction_data = geodata.FrictionData() - friction_data.read(os.path.join(plotdata.outdir, 'friction.data')) - - # Set indices for plotting - surgeplot.friction_field = 1 - surgeplot.wind_field = 2 - surgeplot.pressure_field = 4 - - # Load storm track - track = surgeplot.track_data(os.path.join(plotdata.outdir, 'fort.track')) - - # Set afteraxes function - def surge_afteraxes(cd): - surgeplot.surge_afteraxes(cd, track, plot_direction=False, - kwargs={"markersize": 4}) - - # Color limits - surface_limits = [-5.0, 5.0] - speed_limits = [0.0, 3.0] - wind_limits = [0, 64] - pressure_limits = [935, 1013] - friction_bounds = [0.01, 0.04] - - def friction_after_axes(cd): - plt.title(r"Manning's $n$ Coefficient") - - # ========================================================================== - # Plot specifications - # ========================================================================== - regions = {"domain": {"xlimits": (clawdata.lower[0], clawdata.upper[0]), - "ylimits": (clawdata.lower[1], clawdata.upper[1]), - "figsize": (6.4, 4.8)}} - - for (name, region_dict) in regions.items(): - - # Surface Figure - plotfigure = plotdata.new_plotfigure(name="Surface - %s" % name) - plotfigure.kwargs = {"figsize": region_dict['figsize']} - plotaxes = plotfigure.new_plotaxes() - plotaxes.title = "Surface" - plotaxes.xlimits = region_dict["xlimits"] - plotaxes.ylimits = region_dict["ylimits"] - plotaxes.afteraxes = surge_afteraxes - - surgeplot.add_surface_elevation(plotaxes, bounds=surface_limits) - surgeplot.add_land(plotaxes, bounds=[0.0, 20.0]) - plotaxes.plotitem_dict['surface'].amr_patchedges_show = [0] * 10 - plotaxes.plotitem_dict['land'].amr_patchedges_show = [0] * 10 - - # Speed Figure - plotfigure = plotdata.new_plotfigure(name="Currents - %s" % name) - plotfigure.kwargs = {"figsize": region_dict['figsize']} - plotaxes = plotfigure.new_plotaxes() - plotaxes.title = "Currents" - plotaxes.xlimits = region_dict["xlimits"] - plotaxes.ylimits = region_dict["ylimits"] - plotaxes.afteraxes = surge_afteraxes - - surgeplot.add_speed(plotaxes, bounds=speed_limits) - surgeplot.add_land(plotaxes, bounds=[0.0, 20.0]) - plotaxes.plotitem_dict['speed'].amr_patchedges_show = [0] * 10 - plotaxes.plotitem_dict['land'].amr_patchedges_show = [0] * 10 - - # - # Friction field - # - plotfigure = plotdata.new_plotfigure(name='Friction') - plotfigure.show = friction_data.variable_friction and True - - plotaxes = plotfigure.new_plotaxes() - plotaxes.xlimits = regions['domain']['xlimits'] - plotaxes.ylimits = regions['domain']['ylimits'] - # plotaxes.title = "Manning's N Coefficient" - plotaxes.afteraxes = friction_after_axes - plotaxes.scaled = True - - surgeplot.add_friction(plotaxes, bounds=friction_bounds, shrink=0.9) - plotaxes.plotitem_dict['friction'].amr_patchedges_show = [0] * 10 - plotaxes.plotitem_dict['friction'].colorbar_label = "$n$" - - # - # Hurricane Forcing fields - # - # Pressure field - plotfigure = plotdata.new_plotfigure(name='Pressure') - plotfigure.show = surge_data.pressure_forcing and True - - plotaxes = plotfigure.new_plotaxes() - plotaxes.xlimits = regions['domain']['xlimits'] - plotaxes.ylimits = regions['domain']['ylimits'] - plotaxes.title = "Pressure Field" - plotaxes.afteraxes = surge_afteraxes - plotaxes.scaled = True - surgeplot.add_pressure(plotaxes, bounds=pressure_limits) - surgeplot.add_land(plotaxes, bounds=[0.0, 20.0]) - - # Wind field - plotfigure = plotdata.new_plotfigure(name='Wind Speed') - plotfigure.show = surge_data.wind_forcing and True - - plotaxes = plotfigure.new_plotaxes() - plotaxes.xlimits = regions['domain']['xlimits'] - plotaxes.ylimits = regions['domain']['ylimits'] - plotaxes.title = "Wind Field" - plotaxes.afteraxes = surge_afteraxes - plotaxes.scaled = True - surgeplot.add_wind(plotaxes, bounds=wind_limits) - surgeplot.add_land(plotaxes, bounds=[0.0, 20.0]) - - # ======================================================================== - # Figures for gauges - # ======================================================================== - plotfigure = plotdata.new_plotfigure(name='Gauge Surfaces', figno=300, - type='each_gauge') - plotfigure.show = True - plotfigure.clf_each_gauge = True - - # Set up for axes in this figure: - plotaxes = plotfigure.new_plotaxes() - plotaxes.xlimits = [-2, 1] - # plotaxes.xlabel = "Days from landfall" - # plotaxes.ylabel = "Surface (m)" - plotaxes.ylimits = [-1, 5] - plotaxes.title = 'Surface' - - def gauge_afteraxes(cd): - - axes = plt.gca() - surgeplot.plot_landfall_gauge(cd.gaugesoln, axes) - - # Fix up plot - in particular fix time labels - axes.set_title('Station %s' % cd.gaugeno) - axes.set_xlabel('Days relative to landfall') - axes.set_ylabel('Surface (m)') - axes.set_xlim([-2, 1]) - axes.set_ylim([-1, 5]) - axes.set_xticks([-2, -1, 0, 1]) - axes.set_xticklabels([r"$-2$", r"$-1$", r"$0$", r"$1$"]) - axes.grid(True) - plotaxes.afteraxes = gauge_afteraxes - - # Plot surface as blue curve: - plotitem = plotaxes.new_plotitem(plot_type='1d_plot') - # plotitem.plot_var = 3 - # plotitem.plotstyle = 'b-' - - # - # Gauge Location Plot - # - def gauge_location_afteraxes(cd): - plt.subplots_adjust(left=0.12, bottom=0.06, right=0.97, top=0.97) - surge_afteraxes(cd) - gaugetools.plot_gauge_locations(cd.plotdata, gaugenos='all', - format_string='ko', add_labels=True) - - plotfigure = plotdata.new_plotfigure(name="Gauge Locations") - plotfigure.show = True - - # Set up for axes in this figure: - plotaxes = plotfigure.new_plotaxes() - plotaxes.title = 'Gauge Locations' - plotaxes.scaled = True - plotaxes.xlimits = [-95.5, -94] - plotaxes.ylimits = [29.0, 30.0] - plotaxes.afteraxes = gauge_location_afteraxes - surgeplot.add_surface_elevation(plotaxes, bounds=surface_limits) - surgeplot.add_land(plotaxes, bounds=[0.0, 20.0]) - plotaxes.plotitem_dict['surface'].amr_patchedges_show = [0] * 10 - plotaxes.plotitem_dict['land'].amr_patchedges_show = [0] * 10 - - # ----------------------------------------- - # 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.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 # parallel plotting - - return plotdata diff --git a/examples/storm-surge/synthetic/setrun.py b/examples/storm-surge/synthetic/setrun.py deleted file mode 100644 index 0b2a16832..000000000 --- a/examples/storm-surge/synthetic/setrun.py +++ /dev/null @@ -1,475 +0,0 @@ -#!/usr/bin/env python -# encoding: utf-8 -from __future__ import absolute_import -from __future__ import print_function - -#ensure# encoding: utf-8 -""" -Setup a run for a simple hurricane in a square ocean basin - -""" - -import os -import datetime - -import numpy - -import clawpack.geoclaw.data -import clawpack.geoclaw.topotools as tt -import clawpack.geoclaw.units as units -from clawpack.geoclaw.surge.storm import Storm - -# Initial ramp up time for hurricane -RAMP_UP_TIME = 12 * 60**2 - -def days2seconds(days): - return days * 24 * 60**2 - -# ------------------------------ -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) - - # ------------------------------------------------------------------ - # 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] = -200e3 # west boundary (m) - clawdata.upper[0] = 500e3 # east boundary (m) - - clawdata.lower[1] = -300e3 # south boundary (m) - clawdata.upper[1] = 300e3 # north boundary (m) - - # Number of grid cells: - clawdata.num_cells[0] = 70 - clawdata.num_cells[1] = 60 - - # --------------- - # Size of system: - # --------------- - - # Number of equations in the system: - clawdata.num_eqn = 3 - - # Number of auxiliary variables in the aux array (initialized in setaux) - # First three are from shallow GeoClaw, fourth is friction and last 3 are - # storm fields - clawdata.num_aux = 1 + 1 + 3 - - # Index of aux array corresponding to capacity function, if there is one: - clawdata.capa_index = 0 - - # ------------- - # Initial time: - # ------------- - clawdata.t0 = -RAMP_UP_TIME - - # 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.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. - # 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.tfinal = days2seconds(3) - recurrence = 4 - clawdata.num_output_times = int((clawdata.tfinal - clawdata.t0) * - recurrence / (60**2 * 24)) - - 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 = 1 - clawdata.output_t0 = True - - clawdata.output_format = 'binary' # 'ascii' or 'binary' - clawdata.output_q_components = 'all' # could be list such as [True,True] - clawdata.output_aux_components = 'all' - clawdata.output_aux_onlyonce = False # 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==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.016 - - # 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 = 50000 - - # ------------------ - # 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' - - # 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 - - # max number of refinement levels: - amrdata.amr_levels_max = 3 - - # List of refinement ratios at each level (length at least mxnest-1) - amrdata.refinement_ratios_x = [2,6] - amrdata.refinement_ratios_y = [2,6] - amrdata.refinement_ratios_t = [2,6] - - # 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', 'center', 'center', 'center', 'center'] - - # Flag using refinement routine flag2refine rather than richardson error - amrdata.flag_richardson = False # use Richardson? - 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 = False # time step reporting each level - amrdata.uprint = False # update/upbnd reporting - - # More AMR parameters can be set -- see the defaults in pyclaw/data.py - - # == setregions.data values == - regions = rundata.regiondata.regions - # to specify regions of refinement append lines of the form - # [minlevel,maxlevel,t1,t2,x1,x2,y1,y2] - - # Force the gauges to also record the wind and pressure fields - # rundata.gaugedata.aux_out_fields = [4, 5, 6] - - # ------------------------------------------------------------------ - # GeoClaw specific parameters: - # ------------------------------------------------------------------ - rundata = setgeo(rundata) - - return rundata - # end of function setrun - # ---------------------- - - -# ------------------- -def setgeo(rundata): - """ - Set GeoClaw specific runtime parameters. - For documentation see .... - """ - - geo_data = rundata.geo_data - - # == Physics == - geo_data.gravity = 9.81 - geo_data.coordinate_system = 1 - geo_data.earth_radius = 6367.5e3 - geo_data.rho = 1025.0 - geo_data.rho_air = 1.15 - geo_data.ambient_pressure = 101.3e3 - - # == Forcing Options - geo_data.coriolis_forcing = True - geo_data.theta_0 = 25.0 # Beta-plane approximation center - geo_data.friction_forcing = True - geo_data.friction_depth = 1e10 - - # == Algorithm and Initial Conditions == - # Due to seasonal swelling of gulf we set sea level higher - geo_data.sea_level = 0.0 - geo_data.dry_tolerance = 1.e-2 - - # Refinement Criteria - refine_data = rundata.refinement_data - refine_data.wave_tolerance = 1.0 - refine_data.speed_tolerance = [1.0, 2.0, 3.0, 4.0] - refine_data.variable_dt_refinement_ratios = True - - # == settopo.data values == - topo_data = rundata.topo_data - topo_data.topofiles = [] - # for topography, append lines of the form - # [topotype, fname] - topo_data.topofiles.append([2, 'topo.tt2']) - - # == setfixedgrids.data values == - rundata.fixed_grid_data.fixedgrids = [] - # for fixed grids append lines of the form - # [t1,t2,noutput,x1,x2,y1,y2,xpoints,ypoints,\ - # ioutarrivaltimes,ioutsurfacemax] - - # ================ - # Set Surge Data - # ================ - data = rundata.surge_data - - # Source term controls - data.wind_forcing = True - data.drag_law = 1 - data.pressure_forcing = True - - data.wind_index = 2 - data.pressure_index = 4 - - data.display_landfall_time = True - - # AMR parameters, m/s and m respectively - data.wind_refine = [20.0, 40.0, 60.0] - data.R_refine = [60.0e3, 40e3, 20e3] - - # Storm parameters - Parameterized storm (Holland 1980) - data.storm_specification_type = 'holland80' # (type 1) - data.storm_file = os.path.expandvars(os.path.join(os.getcwd(), - 'fake.storm')) - - # Contruct storm - forward_velocity = units.convert(20, 'km/h', 'm/s') - theta = 0.0 * numpy.pi / 180.0 # degrees from horizontal to radians - - my_storm = Storm() - - # Take seconds and time period of 30 minutes and turn them into datatimes - t_ref = datetime.datetime.now() - t_sec = numpy.arange(-RAMP_UP_TIME, rundata.clawdata.tfinal, 30.0 * 60.0) - my_storm.t = [t_ref + datetime.timedelta(seconds=t) for t in t_sec] - - ramp_func = lambda t: (t + (2 * RAMP_UP_TIME)) * (t < 0) / (2 * RAMP_UP_TIME) \ - + numpy.ones(t_sec.shape) * (t >= 0) - - my_storm.time_offset = t_ref - my_storm.eye_location = numpy.empty((t_sec.shape[0], 2)) - my_storm.eye_location[:, 0] = forward_velocity * t_sec * numpy.cos(theta) - my_storm.eye_location[:, 1] = forward_velocity * t_sec * numpy.sin(theta) - my_storm.max_wind_speed = units.convert(56, 'knots', 'm/s') * ramp_func(t_sec) - my_storm.central_pressure = units.convert(1024, "mbar", "Pa") \ - - (units.convert(1024, "mbar", "Pa") \ - - units.convert(950, 'mbar', 'Pa')) \ - * ramp_func(t_sec) - my_storm.max_wind_radius = [units.convert(8, 'km', 'm')] * t_sec.shape[0] - my_storm.storm_radius = [units.convert(100, 'km', 'm')] * t_sec.shape[0] - - my_storm.write(data.storm_file, file_format="geoclaw") - - - # ======================= - # Set Variable Friction - # ======================= - data = rundata.friction_data - - # Variable friction - data.variable_friction = False - data.friction_index = 1 - - return rundata - # end of function setgeo - # ---------------------- - - -def write_topo_file(run_data, out_file, **kwargs): - """Create simple topography""" - - # Make topography - topo = tt.Topography() - topo.x = numpy.linspace(run_data.clawdata.lower[0], - run_data.clawdata.upper[0], - run_data.clawdata.num_cells[0] + 8) - topo.y = numpy.linspace(run_data.clawdata.lower[1], - run_data.clawdata.upper[1], - run_data.clawdata.num_cells[1] + 8) - - # Create bathymetry profile - beach_slope = 0.05 - basin_depth = -3000 - shelf_depth = -200 - x0 = 350e3 - x1 = 450e3 - x2 = 480e3 - topo_profile = [(run_data.clawdata.lower[0], basin_depth), - (x0, basin_depth), (x1, shelf_depth), (x2, shelf_depth), - (run_data.clawdata.upper[0], - beach_slope * (run_data.clawdata.upper[0] - x2) - + shelf_depth)] - topo.topo_func = tt.create_topo_func(topo_profile) - topo.write(out_file) - - return topo - - -if __name__ == '__main__': - # Set up run-time parameters and write all data files. - import sys - if len(sys.argv) == 2: - rundata = setrun(sys.argv[1]) - else: - rundata = setrun() - - rundata.write() - - topo = write_topo_file(rundata, 'topo.tt2')