****************************************************************************
* *
* SHALLOW WATER EQUATION SOLVER: *
* Use Flux Calculation Approach of D.L. George
* *
****************************************************************************
This README explains how to use the AMR swe sample application
code in the SAMRAI library. The files included in this directory
are as follows:
main.C - application driver program
swe.[C,h] - C++ class containing routines specific to the
swe example and which are called
from SAMRAI AMR library routines
sweFort.h - header file defining external fortran routines
fortran - Fortran numerical routines for 2d swe
To understand the construction of this application code from pieces
specific to the swe model and more generic pieces in the SAMRAI
library consider the object construction pattern in the main.C file.
For more details on these files and the swe class consult the comments
in those files and the swe.h class header file.
COMPILE and RUN
---------------
Compile: make <prog name> where <prog name> = "main" for 2D case
Run:
main <input file> ... Serial run (w/o MPI)
mpirun -np <nprocs> main <input file> ... Parallel run (with MPI)
NOTE: <input file> is the FULL-PATH of the input file
There are sample input files in the subdirectory sample_input.
example:
main sample_input/ramp.input
mpirun -np 2 main sample_input/ramp.input
Restarted Run:
same format as regular run but add <restart dir> <restore num>
to end of command.
<restart dir> is directory where restart files are stored
(set in input file).
<restore num> is the identifier of which restart file to use
(generally the timestep).
example: (to start from timestep #5)
main sample_input/ramp.input restart_ramp 5
mpirun -np 2 main sample_input/ramp.input restart_ramp 5
NOTE: The ordinary use of our restart functionality requires that
a problem can only be restarted on on the same number of
processors on which the problem was originally run. To restart
on a different number of processors, use the
restart-redistribute tool in SAMRAI/tools/restart.
INPUT FILES:
------------
One example input files are included with the swe application:
sample_input/ramp.input - supercritical flow over a ramp
The swe application code is composed of various classes that implement
the algorithmic and numerical routines in the AMR solution process (see
comments in the main.C code for more discussion of the various objects
used to build the application). User-specified inputs control the
behavior of the objects. Those objects expecting input are:
Main program .... creates objects used in the simulation - controls output
of viz and restart information.
swe ........... defines variables and provides numerical kernels for
Euler problem, including data initialization, time
integration, boundary routines, etc.
CartesianGeometry ..... manages Cartesian grid coordinate system
on AMR patch hierarchy.
GriddingAlgorithm ..... constructs AMR patch hierarchy and
regrids levels during simulation.
StandardTagAndInitialize .... manages tagging and initialization during
regridding operation.
HyperbolicLevelIntegrator ... manages time integration and data management
on levels in AMR patch hierarchy.
TimeRefinementIntegrator .... orchestrates time integration and
regridding on levels in AMR patch hierarchy.
LoadBalancer .... manages mapping of patches to processors.
DataWriter .... specifies variables to write for visualization.
TimerManager ..... controls timers used that measure performance of
sections of code.
The input options for these objects are discussed below. The TYPE is
shown in (parenthesis) and the DEFAULT setting (value set if no input
is given) for the option is shown in [brackets]. Required inputs are
designated by [REQD].
"Main" input section
--------------------
Main {
// Log file where output is directed
// log_file_name -- (char string) name of log file ["swe.log"]
// log_all_nodes -- (bool) option to specify whether output on different
// processors should be directed to multiple log files
// (named "log_file_name.<processor id>") or to the
// single file. [FALSE]
log_file_name = "ramp.log"
log_all_nodes = TRUE
// Visualization dump parameters.
// viz_writer -- (char string [, char string]) which visualization
// data writer to use. Use "VisIt" for the VisIt Data
// Writer, or "Vizamrai" for the Vizamrai Data Writer,
/ or "VisIt", "Vizamrai" for both.
// viz_dump_interval -- (int) the interval in timesteps over coarsest
// level between writing of consecutive
// visualization dumps. Setting it to zero turns off
// writing of visualization files.
// viz_dump_dirname -- (char string) directory where viz files are
// writen (may include a path). VisIt requires a
// non-empty string. Vizamrai permits an empty string,
// in which case dumps go to directory where main
// exists. If both data writers are used, directory name
// must be nonempty, and main.C will append _VisIt,
// _Vizamrai to dirname.
// viz_dump_filename -- (char string) Nonempty file prefix required for
// vizamrai filenames. Not used by VisIt.
// visit_number_procs_per_file
// -- (int >= 1) for parallel runs, number of processors
// that share a single common dump file. If this parameter
// > number processors being used, all processors share
// single dump file. Not used by Vizamrai.
//
// For example, the following inputs:
viz_writer = "Vizamrai", "VisIt"
viz_dump_interval = 1
viz_dump_dirname = "viz_room-2d"
viz_dump_filename = "room-2d"
// would result in the following dump structure:
// main
// viz_room-2d_VisIt/
// visit_dump.00000/
// visit_dump.00001/
// ...
// dumps.visit
// viz_room-2d_Vizamrai/
// room-2d.00000.vis
// room-2d.00001.vis
// ...
// If only one data writer is selected, the
// _VisIt, _Vizamrai extension will be omitted.
// Restart dump parameters.
// restart_interval -- (int) the interval in timesteps over coarsest
// level between writing of consecutive restart
// files. Setting it to zero turns off writing
// of restart files. [0]
// restart_write_dirname -- (char string) directory where restart files
// are written [empty string]
// (exception thrown if not specified
// and restart_interval > 0)
//
restart_interval = 1
restart_write_dirname = "restart_room-2d"
// Run-time plotting arguments (only apply in 2d):
//
// X-windows plotting specification for each primitive variable.
// Options: NO_PLOTTING, BOXES_ONLY, CONTOURS_ONLY, SHADE_ONLY,
// CONTOUR_WITH_SHADE, BOTH_CONTOUR_AND_SHADE
// (see "euler_const.input" below for discussion of these options).
// plotting_interval -- (int) interval for which values are plotted
// in X window during run [1]
// plot2d_density -- (above options) specification of plotting for
// density [NO_PLOTTING]
// plot2d_pressure -- (above options) specification of plotting for
// pressure [NO_PLOTTING]
// plot2d_velocity -- (above options) specification of plotting for
// velocity [NO_PLOTTING]
plotting_interval = 5
plot2d_density = NO_PLOTTING
plot2d_pressure = SHADE_ONLY
plot2d_velocity = NO_PLOTTING, NO_PLOTTING
}
"Euler" input section
---------------------
Euler {
// Problem specification parameters
// gamma -- (double) ratio of specific heats [1.4]
// godunov_order -- (int) order of Godunov slopes (1, 2, or 4) [1]
// flux_method -- (one of APPROX_RIEM_SOLVE, EXACT_RIEM_SOLVE,
// see "euler_const.input" for more on these
// options) solver to be used in flux calculation
// across cell face [APPROX_RIEM_SOLVE]
// corner transport scheme -- (CORNER_TRANSPORT_1 -or-
// CORNER_TRANSPORT_2 - see "euler_const.input"
// for more on these options) [CORNER_TRANSPORT_1]
// data_problem -- (SPHERE -or- PIECEWISE_CONSTANT_[X,Y,Z])
// specification of the problem to be solved [REQD]
//
gamma = 1.4
godunov_order = 2
flux_method = APPROX_RIEM_SOLVE
corner_transport = CORNER_TRANSPORT_1
data_problem = SPHERE
// Initial data for SPHERE problem:
// radius -- (double) radius of sphere [REQD]
// center -- (double array) location of sphere center [REQD]
// density_inside -- (double) fluid density inside sphere [REQD]
// velocity_inside -- (double array) fluid velocity inside sphere
// [REQD]
// pressure_inside -- (double) pressure outside sphere [REQD]
// density_outside -- (double) fluid density outside sphere [REQD]
// velocity_outside -- (double array) fluid velocity outside sphere
// [REQD]
// pressure_outside -- (double) pressure outside sphere [REQD]
data_problem = SPHERE
Initial_data {
radius = 2.9
center = 5.5 , 5.5
density_inside = 8.0
velocity_inside = 0.0 , 0.0
pressure_inside = 40.0
density_outside = 1.0
velocity_outside = 0.0 , 0.0
pressure_outside = 1.0
}
// Initial data for PIECEWISE_CONSTANT_[X,Y,Z]:
// front_position -- (double array) points between intervals [REQD]
// interval_n {
// density -- (double) density on interval n [REQD]
// velocity -- (double array) velocity in interval n [REQD]
// pressure -- (double) pressure in interval n [REQD]
// }
//
// interval 0 interval 1 interval n
// |------------|------------------|--------...-------|------------|
// front_position[0] front_position[1] front_position[n-1]
//
data_problem = PIECEWISE_CONSTANT_X
Initial_data {
front_position = 0.0
interval_0 {
density = 1.0
velocity = 0.0 , 0.0
pressure = 1.0
}
interval_1 {
density = 0.125
velocity = 0.0 , 0.0
pressure = 0.1
}
}
// Refinement criteria
// Data for tagging cells to refine for gradient detection and
// Richardson extrapolation. Options:
// DENSITY_DEVIATION, PRESSURE_DEVIATION
// DENSITY_GRADIENT, PRESSURE_GRADIENT
// DENSITY_SHOCK, PRESSURE_SHOCK
// DENSITY_RICHARDSON, PRESSURE_RICHARDSON
// and combinations thereof...
//
// DENSITY_DEVIATION -- tag regions whose density exceeds a specified
// tolerance, preserving tags from earlier step.
// DENSITY_GRADIENT -- tag around gradients in density.
// DENSITY_SHOCK -- tag around discontinuous regions of density.
// DENSITY_RICHARDSON -- use Richardson extrapolation to tag
// around errors in density.
// PRESSURE_GRADIENT -- tag around gradients in pressure.
// PRESSURE_SHOCK -- tag around discontinuous regions of pressure.
// PRESSURE_RICHARDSON -- use Richardson extrapolation to tag
// around errors in pressure.
//
// Refinement_data {
// refine_criteria -- (string array) contains one or more of the
// tagging options, specified above [REQD]
// DENSITY_DEVIATION {
// density_dev -- (double array) freestream density, i.e. tag cells
// where |density - density_dev| > dev_tol [REQD]
// dev_tol -- (double array) deviation tolerance in density [REQD]
// time_min -- (double array) time on each level at which
// tagging using this criteria is started [0.]
// time_max -- (double array) time on each level at which
// tagging is stopped [DBL_MAX]
// NOTE: For each of the above entries, if a level is NOT specified,
// the value from the next coarser level is used. The time_min
// and time_max options may be used to control whether tagging
// on a level is active (i.e. setting time_max=0 makes it
// }
//
// DENSITY_GRADIENT {
// grad_tol -- (double array) gradient tolerance for each level [REQD]
// time_min -- (double array) time on each level at which
// tagging using this criteria is started [0.]
// time_max -- (double array) time on each level at which
// tagging is stopped [DBL_MAX]
// (see NOTE under DENSITY_DEVIATION above)
// }
//
// DENSITY_SHOCK {
// shock_onset -- (double array) onset tolerance for each level [REQD]
// shock_tol -- (double array) gradient tolerance for each level [REQD]
// time_min -- (double array) time on each level at which
// tagging using this criteria is started [0.]
// time_max -- (double array) time on each level at which
// tagging is stopped [DBL_MAX]
// (see NOTE under DENSITY_DEVIATION above)
// }
//
// DENSITY_RICHARDSON {
// rich_tol -- (double array) Richardson extrapolation tolerance for
// each level [REQD]
// time_min -- (double array) time on each level at which
// tagging using this criteria is started [0.]
// time_max -- (double array) time on each level at which
// tagging is stopped [DBL_MAX]
// (see NOTE under DENSITY_DEVIATION above)
// }
//
// PRESSURE_DEVIATION, PRESSURE_GRADIENT, PRESSURE_SHOCK, PRESSURE_RICHARDSON
// input entries are identical to the DENSITY cases shown above.
//
Refinement_data {
refine_criteria = "PRESSURE_GRADIENT", "PRESSURE_SHOCK"
PRESSURE_GRADIENT {
grad_tol = 10.0, 9.0, 9.5
time_min = 0., 0.5, 0.8
time_max = 10.0, 10.0, 8.0
}
PRESSURE_SHOCK {
shock_tol = 5.0, 8.0, 12.0
shock_onset = 0.85, 0.80
time_min = 0., 0.5, 0.8
time_max = 10.0, 10.0, 8.0
}
}
// Boundary conditions
// Options: FLOW, LINEAR, SYMMETRIC, FIXED
// (see "euler_const.input" for discussion of these options)
// Boundary_data {
// boundary_node_x0y0 -- (above options) lower left corner [FLOW]
// boundary_node_x1y0 -- (above options) lower right corner [FLOW]
// boundary_node_x0y1 -- (above options) upper left corner [FLOW]
// boundary_node_x1y1 -- (above options) upper right corner [FLOW]
// boundary_edge_x0 { -- left edge
// boundary_edge_x1 { -- right edge
// boundary_edge_y0 { -- bottom edge
// boundary_edge_y1 { -- top edge
// boundary_type -- (above options) type at prescribed edge [FLOW]
// density -- (double) density at prescribed edge [0]
// velocity -- (double array) velocity at prescribed edge [0]
// pressure -- (double) pressure at prescribed edge [0]
// }
//
// edge_y1
// node_x0y1 o--------------o node_x1y1
// | |
// | |
// edge_x0 | | edge_x1
// | |
// | |
// node_x0y0 o--------------o node_x1y0
// edge_y0
//
// NOTE: 3D case specifies data at faces as well. See
// /SAMRAI/docs/userdocs/boundaries.ps for further discussion.
Boundary_data {
boundary_node_x0y0 = FLOW
boundary_node_x1y0 = FLOW
boundary_node_x0y1 = FLOW
boundary_node_x1y1 = FLOW
boundary_edge_x0 {
boundary_type = SYMMETRIC
}
boundary_edge_x1 {
boundary_type = SYMMETRIC
}
boundary_edge_y0 {
boundary_type = SYMMETRIC
}
boundary_edge_y1 {
boundary_type = SYMMETRIC
}
}
}
"CartesianGeometry" input section
---------------------------------
CartesianGeometry {
// Specify lower/upper corners of the computational domain and a
// set of non-overlapping boxes defining domain interior. If union
// of boxes is not a parallelpiped, lower/upper corner data corresponds
// to min/max corner indices over all boxes given.
// x_lo -- (double array) lower corner of computational domain [REQD]
// x_up -- (double array) upper corner of computational domain [REQD]
// domain_boxes -- (box array) set of boxes that define interior of
// hysical domain. [REQD]
// periodic_dimension -- (int array) coordinate directions in which
// domain is periodic. Zero indicates not
// periodic, non-zero value indicates periodicity.
// [0]
x_lo = 0.e0 , 0.e0
x_up = 65.e0 , 55.e0
domain_boxes = [(0,9),(14,20)],
[(0,21),(19,35)],
[(0,36),(14,49)],
[(20,0),(49,54)],
[(50,0),(64,35)],
[(50,46),(64,54)],
[(57,36),(64,45)]
periodic_dimension = 0, 1 // periodic in y only
}
"GriddingAlgorithm" input section
---------------------------------
GriddingAlgorithm {
// Information used to create patches in AMR hierarchy.
// max_levels -- (int) max number of mesh levels in hierarchy [REQD]
//
// For most of the following parameters, the number of precribed data
// values need not match the number of levels in the hierarchy
// (determined by max_levels). If more values are given than number
// of levels, extraneous values will be ignored. If less are give, then
// values that correspond to individual levels will apply to those
// levels. Missing values will be taken from those for the finest
// level specified.
//
// ratio_to_coarser {
// level_1 -- (int array) ratio between index spaces on
// level 1 to level 0 [REQD]
// level_2 -- (int array) ratio between index spaces on
// level 2 to level 1 [REQD]
// etc....
// }
// largest_patch_size {
// level_0 -- (int array) largest patch allowed on level 0.
// [REQD]
// level_1 -- (int array) " " " " level 1
// [level 0 entry]
// etc....
// }
// smallest_patch_size {
// level_0 -- (int array) smallest patch allowed on level 0.
// [max ghost width for all variables]
// level_1 -- (int array) " " " " level 1
// [level 0 entry]
// etc....
// }
max_levels = 3
ratio_to_coarser {
level_1 = 4 , 4
level_2 = 4 , 4
level_3 = 4 , 4
}
largest_patch_size {
level_0 = 48 , 48
}
smallest_patch_size {
level_0 = 8 , 8
}
// Tolerances for gridding efficiency and box chopping operations.
// efficiency_tolerance -- (double array) minimum percentage of tagged
// cells allowed in a box [0.8 for each level]
// combine_efficiency -- (double array) threshold specifying when a
// box may be chopped into two smaller boxes.
// If sum of the volumes of smaller boxes is >
// combine efficiency * volume of larger box,
// then larger box will not be chopped.
// [0.8 for each level]
// proper_nesting_buffer -- (int array) number of coarse cells by which
// the next finer level is nested within its
// interior. [1 for each level]
efficiency_tolerance = 0.85e0
combine_efficiency = 0.95e0
proper_nesting_buffer = 1
// Option to read or write gridding information
// write_regrid_boxes -- (bool) Output sequence of refine boxes to file
// read_regrid_boxes -- (bool) Read sequence of refine boxes from file
// regrid_boxes_filename -- (string) file name used for reading or writing
// boxes.
//
// Reading and writing options require some user intervention to assure
// they work properly. Please consult Andy Wissink (awissink@llnl.gov)
// if you are interested in using these options.
write_regrid_boxes = TRUE
regrid_boxes_filename = TRUE
}
"StandardTagAndInitialize" input section
---------------------------------------------
StandardTagAndInitialize {
// Specification of the type of tagging to be performed
// tagging_method -- (string array) one or more entries specifying
// the tagging algorithm used. Options include:
// "GRADIENT_DETECTOR"
// "RICHARDSON_EXTRAPOLATION"
// "REFINE_BOXES"
//
// REFINE_BOXES allows you to prescribe where refinement should occur.
// If this option is used, you must also supply a RefineBoxes database
// entry with the refine boxes on different levels specified. i.e.
// RefineBoxes{
// level_0 -- (BoxArray) Boxes to be refined on level 0
// level_1 -- (BoxArray) Boxes to be refined on level 1
// etc.
//
// Gradient detection option.
tagging_method = "GRADIENT_DETECTOR"
//
// Richarson extrapolation option.
tagging_method = "RICHARDSON_EXTRAPOLATION"
//
// Static refine boxes option.
tagging_method = "REFINE_BOXES"
RefineBoxes {
level0_boxes = [(15,0),(29,14)]
level1_boxes = [(65,10),(114,40)]
}
// Combination.
tagging_method = "RICHARDSON_EXTRAPOLATION", "GRADIENT_DETECTOR"
}
"HyperbolicLevelIntegrator" input section
-----------------------------------------
HyperbolicLevelIntegrator{
// CFL used to determine timestep on each level.
// cfl -- (double) timestep at each level is given by
// dt = cfl*dx/a, where a is max wave speed. [REQD]
// cfl_init -- (double) cfl factor used for initial step [REQD]
cfl = 0.9e0
cfl_init = 0.1e0
// Options for timestep computation.
// lag_dt_computation -- (bool) if true, characteristic data on the
// patch is not recomputed before computing the
// timestep; otherwise it uses data from
// preceding flux calculation [TRUE].
// use_ghosts_to_compute_dt -- (bool) if true, fill ghost data before
// computing timestep on each patch [TRUE].
lag_dt_computation = TRUE
use_ghosts_to_compute_dt = TRUE
}
"TimeRefinementIntegrator" input section
-----------------------------------------
TimeRefinementIntegrator{
// Time stepping specifications.
// start_time -- (double) initial integration time [REQD]
// end_time -- (double) final integration time [REQD]
// grow_dt -- (double) limit factor on maximum growth of timestep
// from one step to the next; useful for ramping up
// timestep size or to prevent rapid growth in timestep
// [REQD]
// max_integrator_steps -- max number timesteps on coarsest level [REQD]
// tag_buffer -- (int array) amount by which cells tagged for error are
// buffered before new fine patches are created to cover
// tagged regions. By default, it is the interval of
// timesteps between successive regids on level.
// [regrid_interval]
start_time = 0.e0
end_time = 100.e0
grow_dt = 1.1e0
max_integrator_steps = 301
tag_buffer = 1,1,1,1 // max of 4 finer levels in hierarchy
}
"LoadBalancer" input section
----------------------------------
LoadBalancer {
// using default uniform load balance configuration (see mesh_LoadBalancerX.h.sed)
}
"DataWriter" input section
----------------------------------
DataWriter {
// Specify which variables to write to visualization file(s).
// Variable_Names -- (string array) list of variable names to
// write to the viz file.
// Variable_Types -- (string array) list of the corresponding
// variable type to write to the viz file.
// Must be either "primitive" or "derived".
// Note that until vizamrai is updated to use
// HDF5, "derived" must always be used.
// Variable_Depths -- (int array) depth, or number of components, of
// each variable.
// Variable_Scales -- (double array) list of scaling factors for each
// variable to be written to the viz file.
Variable_Names = "density", "velocity", "pressure"
Variable_Types = "Derived", "Derived", "Derived"
Variable_Depths = 1, 3, 1
Variable_Scales = 1.0, 1.0, 1.0
}
"TimerManager" input section
----------------------------
TimerManager{
// Specify which timers to invoke
// timer_list -- (string array) list of timers to be turned on - for
// complete list of available timers, see
// SAMRAI/docs/userdocs/timing.pdf [null]
timer_list = "apps::main::*",
"apps::Euler::computeFluxesOnPatch()",
"algs::GriddingAlgorithm::*",
"xfer::*::*"
// Various output options.
// print_total -- (bool) print total time [TRUE]
// print_exclusive -- (bool) output exclusive time - i.e. time spent
// exclusively in each routine. A nice way to
// identify routines with bottlenecks, but is more
// expensive than recording total time [FALSE].
// print_wall -- (bool) output wallclock time [TRUE]
// print_user -- (bool) output user time [FALSE]
// print_sys -- (bool) output system time [FALSE]
// print_summed -- (bool) output sum of times measured across all
// processors during a parallel run [FALSE].
// print_max -- (bool) output maximum time spent on any one
// processor, along with the processor ID [FALSE]
// print_percentage -- (bool) output percentage of total time for each
// printed timer [TRUE]
// print_concurrent -- (bool) in the process of performing exclusive
// time, a nested call tree is created; print the
// tree [FALSE]
// print_timer_overhead -- (bool) each timer start/stop invocation
// incurs about one-half a millionth of a second
// overhead. This option will output the number of
// invocations and the total predicted overhead for
// all the timers [FALSE]
// print_threshold -- (double) do not print any timers that constitute
// < (print_threshold)% of overall run time. Useful
// to prevent gobs of output if many timers are
// turned on. [0.25]
print_exclusive = TRUE
print_max = TRUE
print_timer_overhead = TRUE
print_threshold = 0.0
}
"euler_const.input"
-------------------
// This file specifies some constants used within the Euler code. These
// names and values map to those specified via "#define" statements in
// Euler.C
// 2D X Window plotting options
// Plotting will execute while the program is running for 2D cases.
// Note that since it pipes data to an X-window, it slows the
// performance of the code somewhat, but it sure is PRETTY!
NO_PLOTTING = 0 // no plotting
BOXES_ONLY = 1 // plot bounding boxes only
CONTOURS_ONLY = 2 // plot contours only
SHADE_ONLY = 3 // plot color shade only
CONTOUR_WITH_SHADE = 4 // plot contours and color shade in same window
BOTH_CONTOUR_AND_SHADE = 5 // plot contours and color shade in separate
// window
// Riemann solver used in flux calculation across cell face
APPROX_RIEM_SOLVE = 0 // Colella-Glaz approx Riemann solver
EXACT_RIEM_SOLVE = 1 // Exact Riemann solver
// Flux corner transport options
// CORNER_TRANSPORT_1 is based on an extension of Colella's formulation.
// CORNER_TRANSPORT_2 is a formulation constructed by Trangenstein
// A writeup on these two methods is available from R. Hornung or A. Wissink.
CORNER_TRANSPORT_1 = 1
CORNER_TRANSPORT_2 = 2
//
// Two Different problems can be specified via input
// PIECEWISE_CONSTANT_* : piecewise constant data in *-coordinate direction
// SPHERE : sphere problem
PIECEWISE_CONSTANT_X = 10
PIECEWISE_CONSTANT_Y = 11
PIECEWISE_CONSTANT_Z = 12
SPHERE = 40
// Boundary conditions
// FLOW indicates flow-thru conditions
// LINEAR indicates a linear profile at the boundary
// SYMMETRIC will apply a plane of symmetry
// FIXED designates a solid boundary.
FLOW = 0
LINEAR = 1
SYMMETRIC = 4
FIXED = 5