|
setrun.py.html |
 |
|
Source file: setrun.py
|
|
Directory: /Users/rjl/clawpack_src/clawpack_master/geoclaw/examples/bouss/radial_flat
|
|
Converted: Sat Mar 23 2024 at 11:13:12
using clawcode2html
|
|
This documentation file will
not reflect any later changes in the source file.
|
"""
Module to set up run time parameters for Clawpack.
The values set in the function setrun are then written out to data files
that will be read in by the Fortran code.
"""
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')
# used to create ruled rectangle:
from clawpack.amrclaw import region_tools
#------------------------------
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] = 0. # xlower
clawdata.upper[0] = 5000. # xupper
clawdata.lower[1] = 0. # ylower
clawdata.upper[1] = 5000. # yupper
# Number of grid cells:
clawdata.num_cells[0] = 100 # mx
clawdata.num_cells[1] = 100 # my
# ---------------
# Size of system:
# ---------------
# Number of equations in the system:
clawdata.num_eqn = 5
# Number of auxiliary variables in the aux array (initialized in setaux)
clawdata.num_aux = 1
# Index of aux array corresponding to capacity function, if there is one:
clawdata.capa_index = 0
# -------------
# 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.chk00225' # 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 ntimes frames at equally spaced times up to tfinal:
# Can specify num_output_times = 0 for no output
clawdata.num_output_times = 20
clawdata.tfinal = 200.
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 = [251.] + list(np.arange(4.5, 17.1, 0.5)*60)
#clawdata.output_times = [251.] + list(np.arange(10,35,2)*60)
#list(np.arange(35,48.5,0.5)*60)
#list(np.arange(35,60.5,0.5)*60)
clawdata.output_t0 = True # output at initial (or restart) time?
elif clawdata.output_style == 3:
# Output every iout timesteps with a total of ntot time steps:
clawdata.output_step_interval = 1
clawdata.total_steps = 10
clawdata.output_t0 = True
#clawdata.output_t0 = False
clawdata.output_format = 'binary' # '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 = .001
clawdata.dt_initial = .03
# 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.limiter = [0,0,0]
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 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] = 'wall' # at xlower
clawdata.bc_upper[0] = 'extrap' # at xupper
clawdata.bc_lower[1] = 'wall' # 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]
# --------------
# Checkpointing:
# --------------
# Specify when checkpoint files should be created that can be
# used to restart a computation.
clawdata.checkpt_style = 1
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 = 50
# ---------------
# 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 = 2
# List of refinement ratios at each level (length at least amr_level_max-1)
amrdata.refinement_ratios_x = [4,2]
amrdata.refinement_ratios_y = [4,2]
amrdata.refinement_ratios_t = [4,2]
# 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']
# Flag using refinement routine flag2refine rather than richardson error
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 = 2
# width of buffer zone around flagged points:
# (typically the same as regrid_interval so waves don't escape):
amrdata.regrid_buffer_width = 3
# 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
# ---------------
# AMR flagregions:
# ---------------
flagregions = rundata.flagregiondata.flagregions # initialized to []
# now append as many flagregions as desired to this list:
from clawpack.amrclaw.data import FlagRegion
# The entire domain restricted to level 1 for illustration:
# Note that this is a rectangle specified in the new way:
# (other regions below will force/allow more refinement)
flagregion = FlagRegion(num_dim=2)
flagregion.name = 'Region_domain'
flagregion.minlevel = 1
flagregion.maxlevel = 1
flagregion.t1 = 0.
flagregion.t2 = 1e9
flagregion.spatial_region_type = 1 # Rectangle
flagregion.spatial_region = [clawdata.lower[0], clawdata.upper[0], \
clawdata.lower[1], clawdata.upper[1]]
flagregions.append(flagregion)
# A ruled rectangle covering abs(x-y) < 1000:
flagregion = FlagRegion(num_dim=2)
flagregion.name = 'Region_diagonal'
flagregion.minlevel = 1
flagregion.maxlevel = 2
flagregion.t1 = 0.
flagregion.t2 = 1e9
flagregion.spatial_region_type = 2 # Ruled Rectangle
flagregion.spatial_region_file = \
os.path.abspath('RuledRectangle_Diagonal.data')
flagregions.append(flagregion)
# code to make RuledRectangle_Diagonal.data:
rr = region_tools.RuledRectangle()
rr.method = 1 # piecewiselinear edges between s values
rr.ixy = 'x' # so s refers to x, lower & upper are limits in y
rr.s = np.array([0, 5000.])
rr.lower = np.array([-1000, 4000.])
rr.upper = np.array([1000., 6000.])
rr.write('RuledRectangle_Diagonal.data')
# Region near the origin where level 3 is allowed:
flagregion = FlagRegion(num_dim=2)
flagregion.name = 'Region_domain'
flagregion.minlevel = 1
flagregion.maxlevel = 3
flagregion.t1 = 0.
flagregion.t2 = 1e9
flagregion.spatial_region_type = 1 # Rectangle
flagregion.spatial_region = [0,1000,0,1000]
flagregions.append(flagregion)
# A ruled rectangle covering abs(x-y) < 1000 for x,y < 2000:
# Only relevant if amr_levels_max >= 3 to capture trailing waves better
flagregion = FlagRegion(num_dim=2)
flagregion.name = 'Region_diagonal2'
flagregion.minlevel = 1
flagregion.maxlevel = 3
flagregion.t1 = 0.
flagregion.t2 = 1e9
flagregion.spatial_region_type = 2 # Ruled Rectangle
flagregion.spatial_region_file = \
os.path.abspath('RuledRectangle_Diagonal2.data')
#flagregions.append(flagregion)
# code to make RuledRectangle_Diagonal2.data:
rr = region_tools.RuledRectangle()
rr.method = 1 # piecewiselinear edges between s values
rr.ixy = 'x' # so s refers to x, lower & upper are limits in y
rr.s = np.array([0, 1000., 2000.])
rr.lower = np.array([0, 0, 2000.])
rr.upper = np.array([1000., 2000., 2000.])
rr.write('RuledRectangle_Diagonal2.data')
# ----- 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 = 1
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 = 1.e-3
geo_data.friction_forcing = True
geo_data.manning_coefficient = 0.025
geo_data.friction_depth = 100.0
# Refinement settings
refinement_data = rundata.refinement_data
refinement_data.variable_dt_refinement_ratios = True
refinement_data.wave_tolerance = 0.02
# == settopo.data values ==
topofiles = rundata.topo_data.topofiles
# for topography, append lines of the form
# [topotype, fname]
topofiles.append([1, 'flat100.tt1'])
# == setdtopo.data values ==
dtopo_data = rundata.dtopo_data
# for moving topography, append lines of the form : (<= 1 allowed for now!)
# [topotype, fname]
# == setqinit.data values ==
rundata.qinit_data.qinit_type = 0
rundata.qinit_data.qinitfiles = []
qinitfiles = rundata.qinit_data.qinitfiles
# for qinit perturbations, append lines of the form: (<= 1 allowed for now!)
# [minlev, maxlev, fname]
# To use Boussinesq solver, add bouss_data parameters here
# Also make sure to use the correct Makefile pointing to bouss version
# and set clawdata.num_eqn = 5
from clawpack.geoclaw.data import BoussData
rundata.add_data(BoussData(),'bouss_data')
rundata.bouss_data.bouss_equations = 2 # 0=SWE, 1=MS, 2=SGN
rundata.bouss_data.bouss_min_level = 1 # coarsest level to apply bouss
rundata.bouss_data.bouss_max_level = 10 # finest level to apply bouss
rundata.bouss_data.bouss_min_depth = 1. # depth to switch to SWE
rundata.bouss_data.bouss_solver = 3 # 1=GMRES, 2=Pardiso, 3=PETSc
rundata.bouss_data.bouss_tstart = 0. # time to switch from SWE
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)