Contents

Acknowledgments

The FLASH Code Group is supported by the (ASC) FLASH Center at the University of Chicago under U. S. Department of Energy contract B341495. Some of the test calculations described here were performed using the Origin 2000 computer at Argonne National Laboratory and the (ASC) Nirvana computer at Los Alamos National Laboratory.

Part 1
PRELIMINARIES

1  Introduction

FLASH is a modular, adaptive-mesh, parallel simulation code capable of handling general compressible flow problems found in many astrophysical environments. FLASH is designed to allow users to configure initial and boundary conditions, change algorithms, and add new physics modules with minimal effort. It uses the PARAMESH library to manage a block-structured adaptive grid, placing resolution elements only where they are needed most. FLASH uses the Message-Passing Interface (MPI) library to achieve portability and scalability on a variety of different parallel computers.

The Center for Astrophysical Thermonuclear Flashes, or FLASH Center, was founded at the University of Chicago in 1997 under contract to the United States Department of Energy as part of its Accelerated Strategic Computing Initiative ((ASC)). The goal of the Center is to address several problems related to thermonuclear flashes on the surfaces of compact stars (neutron stars and white dwarfs), in particular X-ray bursts, Type Ia supernovae, and novae. To solve these problems requires the participants in the Center to develop new simulation tools capable of handling the extreme resolution and physical requirements imposed by conditions in these explosions, and to do so while making efficient use of the parallel supercomputers developed by the ASC project, the most powerful constructed to date.

1.1  What's new in FLASH 2.0

The FLASH 2.0 code has been greatly expanded from the original FLASH 1.0 (Fryxell et al. 2000) and demonstrates significant progress toward the (ASC) FLASH Center's goal in the form of increased problem-solving capability, increased modularization, and further development of the code. FLASH 2.0 includes the following new features:

• An increased number of physics modules, including MHD, self-gravity, and diffusion;
• New hydro solvers in addition to the existing PPM-based solver;
• Numerous new test problem setups;
• A new setup script written in Python which includes more options, including the ability to produce an HTML listing of all runtime parameters;
• Support for HDF5;
• Use of variable databases to improve encapsulation;
• New fixes for bugs reported in FLASH1.6;
• Improved and greatly expanded documentation.

1.2  About the user's guide

This user's guide is designed to enable individuals unfamiliar with the FLASH code to quickly get acquainted with its structure and to move beyond the simple test problems distributed with FLASH, customizing it to suit their own needs. Section 2 discusses how to get started quickly with FLASH, describing how to configure, build, and run the code with one of the included test problems, then examine the resulting output. Users familiar with the capabilities of FLASH who wish to quickly `get their feet wet' with the code should begin with this section.

Part II begins with an overview of the FLASH code architecture. It also includes a brief overview of the modules which are described individually and in more detail in Part III. Part II also explains how to use the driver-supplied and simulated services.

Part III describes in detail each of the modules included with the code, along with their submodules, runtime parameters, use with included solvers, and the equations and algorithms they use. Important note: We assume that the reader has some familiarity both with the basic physics involved and with numerical methods for solving partial differential equations. This familiarity is absolutely essential in using FLASH (or any other simulation code) to arrive at meaningful solutions to physical problems. The novice reader is directed to an introductory text, examples of which include

       Fletcher, C. A. J.  Computational Techniques for Fluid Dynamics (Springer-Verlag, 1991)

       Laney, C. B. Computational Gasdynamics (Cambridge UP, 1998)

       LeVeque, R. J., Mihalas, D., Dorfi, E. A., and Müller, E., eds. Computational Methods for Astrophysical Fluid Flow (Springer, 1998)

       Roache, P. Fundamentals of Computational Fluid Dynamics (Hermosa, 1998)

       Toro, E. F.  Riemann Solvers and Numerical Methods for Fluid Dynamics, 2nd Edition (Springer, 1999)

The advanced reader who wishes to know more specific information about a given module's algorithm is directed to the literature referenced in the algorithm section of the chapter in question.

Part IV describes the different test problems distributed with FLASH. Part V describes in more detail the use of the configuration and analysis tools distributed with FLASH. Finally, Part VI gives detailed instructions for extending FLASH's capabilities by adding new problem setups, describes how new solvers may be integrated into the code, and Part VII lists publication references.

2  Quick start

This section describes how to quickly get up and running with FLASH, showing how to configure and build it to solve the Sedov explosion problem, how to run it, and how to examine the output using IDL.

2.1  System requirements

You should verify that you have the following:

• A copy of the FLASH source code distribution. This is most likely available either as a Unix tar file or as a local Concurrent Versions System (CVS) source tree. To request a copy of the distribution, click on the ``Code Request'' link at the FLASH Code Group web site, http://flash.uchicago.edu/flashcode/. You will be asked to fill out a short form before receiving download instructions. Please remember the user name and password you use to download the code: you will need these to get bug fixes and updates to FLASH.
• A Fortran 90 compiler and a C compiler. Most of FLASH is written in Fortran 90. Information available at the Fortran Market web site (http://www.fortran.com/) can help you select a Fortran 90 compiler for your system.
• An installed copy of the Message-Passing Interface (MPI) library. A freely available implementation of MPI called MPICH has been created at Argonne National Laboratory and can be accessed on the World Wide Web at
  http://www-unix.mcs.anl.gov/mpi/mpich/.
• If you plan to use the Hierarchical Data Format (HDF) for output files (the default), you will need an installed copy of the freely available HDF library. Currently FLASH supports HDF versions 4 and 5 (the two formats are not compatible). HDF is available from the HDF Project of the National Computational Science Alliance at http://hdf.ncsa.uiuc.edu/. The contents of HDF output files produced by the FLASH io/hdf* modules are described in detail in section 5.
• To use the output analysis tools described in this section, a copy of the IDL language from Research Systems, Inc. (http://www.rsinc.com/). IDL is a commercial product. It is not required for the analysis of FLASH output, but without it the user will need to write his or her own analysis tools. (FLASH output formats are described in Section 5.) The FLASH code group is currently working with members of the Mathematics and Computer Science Division at Argonne National Laboratory to develop more sophisticated visualization tools to distribute either as part of or alongside FLASH.
• The GNU make utility, gmake. This utility is freely available and has been ported to a wide variety of different systems. (For more information, see the entry for make in the development software listing at http://www.gnu.org/.) On some systems make is an alias for gmake. GNU make is required because FLASH 1.6 and higher uses macro concatenation when constructing Makefiles.
• To run the FLASH 2.0 Python setup script, a copy of the Python language, version 1.5.2 or later. Several different versions of Python are freely available at
  http://www.python.org.

FLASH has been tested on the following Unix-based platforms. In addition, it may work with others not listed (see Section ).

• SGI single- and multiprocessor systems running IRIX
• Intel- and Alpha-based single and multiprocessor systems running Linux, including clusters
• Cray/SGI T3E running UNICOS
• The (ASC) Nirvana machine, built by SGI
• IBM SP2 systems including the (ASC) Blue Pacific, Frost, Quad (at Argonne), and Blue Horizon machines
• The (ASC) Red machine, built by Intel

2.2  Unpacking and configuring FLASH for quick start

Next, configure the FLASH source tree for the Sedov explosion problem. Type

./setup sedov -auto

From the FLASH root directory (i.e. the directory from which you ran setup), execute gmake. This will compile the FLASH code. If you should have problems and need to recompile, `gmake clean' will remove all object files from the object/ directory, leaving the source configuration intact; `gmake realclean' will remove all files and links from object/. After `gmake realclean,' a new invocation of setup is required before the code can be built.

Assuming compilation and linking were successful, you should now find an executable named flashX in the object/ directory, where X is the major version number (e.g., 2 for X.Y = 2.0). You may wish to check that this is the case.

FLASH expects to find a flat text file named flash.par in the directory from which it is run. This file sets the values of various runtime parameters which determine the behavior of FLASH. If it is not present, FLASH will abort; flash.par must be created in order for the program to run. Here we will create a simple flash.par which sets a few parameters and allows the rest to take on default values. With your text editor, create flash.par in the main FLASH directory with the contents of Figure 1.

# runtime parameters lrefine_max = 4 basenm = ßedov_4_" restart = .false. trstrt = 0.005 nend = 1000 tmax = 0.02 gamma = 1.4 xl_boundary_type = öutflow" xr_boundary_type = öutflow" yl_boundary_type = öutflow" yr_boundary_type = öutflow" plot_var_1 = "dens" plot_var_2 = "temp" plot_var_3 = "pres"

This instructs FLASH to use up to four levels of adaptive mesh refinement (AMR) and to name the output files appropriately. We will not be starting from a checkpoint file (this is the default, but here it is explicitly set for purposes of illustration). Output files are to be written every 0.005 time units and will be created until t=0.02 or 1000 timesteps have been taken, whichever comes first. The ratio of specific heats for the gas (g) is taken to be 1.4, and all four boundaries of the two-dimensional grid have outflow (zero-gradient or Neumann) boundary conditions. Note the format of the file: each line is a comment (denoted by a hash mark, #), blank, or of the form variable = value. String values are enclosed in double quotes ("). Boolean values are indicated in the Fortran style, .true. or .false. Be sure to insert a carriage return after the last line of text. A full list of the parameters available for your current setup is contained in paramFile.html, which also includes brief comments for each parameter.

2.3  Running FLASH

We are now ready to run FLASH. To run FLASH on N processors, type

mpirun -np N object/flashX

You should see a number of lines of output indicating that FLASH is initializing the Sedov problem, listing the initial parameters, and giving the timestep chosen at each step. After the run is finished, in the current directory you should find several files:

• flash.log echoes the runtime parameter settings and indicates the run time, the build time, and the build machine. During the run, a line is written for each timestep, along with any warning messages. If the run terminates normally, a performance summary is written to this file. Messages indicating when the code refined and what output resulted are also contained in flash.log.
• flash.dat contains a number of integral quantities as functions of time: total mass, total energy, total momentum, etc. This file can be used directly by plotting programs such as gnuplot; note that the first line begins with a hash (#) and thus is ignored by gnuplot.
• sedov_4_hdf_chk_000* are the different checkpoint files. These are complete dumps of the entire simulation at intervals of trstrt, suitable for use in restarting the simulation. They are also the primary output products of FLASH.
• sedov_4_hdf_plt_cnt_000* are plot files. These are files containing only density, temperature, and pressure (in single precision, for some I/O modules). They are designed to be written more frequently than checkpoint files for the purpose of making simulation movies (or for analyses that do not require all of the checkpointed quantities).
• amr_log includes various messages from the PARAMESH package.

setenv XFLASH_DIR "$PWD/tools/idl"
setenv IDL_PATH "${XFLASH_DIR}:$IDL_PATH"

setenv IDL_PATH "$XFLASH_DIR":idl-root-path

When the control panel settings are to your satisfaction, click the `Plot' button to generate the plot. For Postscript and GIF output, a file is created in the current directory. The result should look something like Figure 2, although this figure was generated from a run with eight levels of refinement rather than the four used in the quick start example run. With fewer levels of refinement the Cartesian grid causes the explosion to appear somewhat diamond-shaped.

FLASH is intended to be customized by the user to work with interesting initial and boundary conditions. In the following sections we will cover in more detail the algorithms and structure of FLASH and the sample problems and tools distributed with it.

Part 2
STRUCTURE

3  Overview of FLASH architecture

The FLASH source code is a collection of components called FLASH modules. FLASH modules can be combined in a variety of ways to form specific FLASH applications. Of course, not all available FLASH modules are necessarily used when solving any one particular problem. Thus, it is important to distinguish between the entire FLASH source code and a given FLASH application.

3.1  Structure of a FLASH module

Most generally, a FLASH module represents some well-defined, top-level unit of functionality useful for a given class of problems. Its structure conforms to a small set of rules that facilitate its interactions with other modules in the creation of an application. Primary among these are the rules governing the retrieval and modification of data on the solution grid. A module must also announce a general set of requirements to the framework as well as publish a public interface of its services. Here we focus on the internal structure of a FLASH module appropriate for users wishing to extend the current FLASH functionality.

First, it is important to recall how a selected group of FLASH modules is combined to form a particular application. This process is carried out entirely by the FLASH setup tool, which uses configuration information provided by the modules and problem setup to properly parse the source tree and isolate the source files needed to carry set-up a specific problem. For performance reasons, setup ties modules together statically before the application is compiled.

Each FLASH module is divided into three principal components:

a) Configuration layer
b) ``Wrapper'' or ``interface'' layer
c) Algorithm

Additionally, a module may contain sub-modules which inherit from and override the functionality in the parent module. Each of these components is discussed in detail in the following sections.

3.1.1  Configuration layer

Information about module dependencies, default sub-modules, runtime parameter definitions, library requirements, and so on is stored in plain text files named Config in the different module directories. These are parsed by setup when configuring the source tree and are used to create the code needed to register module variables, implement the runtime parameters, choose sub-modules when only a generic module has been specified, prevent mutually exclusive modules from being included together, and to flag problems when dependencies are not resolved by some included module. In the future they may contain additional information about module interrelationships.

3.1.1.1  Configuration file syntax  

The syntax of the configuration files is as follows. Arbitrarily many spaces and/or tabs may be used, but all keywords must be in uppercase. Lines not matching an admissible pattern are ignored. (Someday soon they will generate a syntax error.)

• # comment
  A comment. Can appear at the end of a line.
• DEFAULT sub-module
  Specifies which sub-module of the current module is to be used as a default if a specific sub-module has not been indicated in the modules file (Section 6.1.3). For example, the Config file for the materials/eos module specifies gamma as the default sub-module. If no sub-module is explicitly included (ie. INCLUDE materials/eos is placed in modules), then this command instructs setup to assume that the gamma submodule was meant (as though INCLUDE materials/eos/gamma had been placed in modules).
• EXCLUSIVE sub-module...
  Specify a list of sub-modules that cannot be included together. If no EXCLUSIVE instruction is given, it is perfectly legal to simultaneously include more than one sub-module in the code.
• REQUIRES module[/sub-module[/sub-module...]]
  Specify a module requirement. Module requirements can be general, not specifying sub-modules, so that module dependencies can be independent of particular algorithms. For example, the statement REQUIRES materials/eos in a module's Config file indicates to setup that the materials/eos module is needed by this module. No particular equation of state is specified, but some EOS is needed, and as long as one is included by modules, the dependency will be satisfied. More specific dependencies can be indicated by specifying sub-modules. For example, materials/eos is a module with several possible sub-modules, each corresponding to a different equation of state. To specify a requirement for, eg., the Helmholtz EOS, use REQUIRES materials/eos/helmholtz. Giving a complete set of module requirements is helpful to the end user, because setup uses them to generate the modules file when invoked with the -auto option.
• PARAMETER name type default
  Specify a runtime parameter. Parameter names are unique up to 20 characters and may not contain spaces. Admissible types include REAL (or DOUBLE), INTEGER, STRING, and BOOLEAN (or LOGICAL). Default values for REAL and INTEGER parameters must be valid numbers, or compilation will fail. Default STRING values can be unquoted (in which case only the first word is recognized) or enclosed in double quotes ("). Default BOOLEAN values must be TRUE or FALSE to avoid compilation errors. Once defined, runtime parameters are available to the entire code.
• VARIABLE name [attribute_list]
  Register variable with the framework with name name and attributes defined by attribute_list. These variables can later be accessed programatically using the database accessor methods (see 3.1.2). Valid attributes are as follows:

  • ADVECT/NOADVECT
    A variable, Q, with the ADVECT property obeys an advection equation,
    

    ¶Q ¶t + Ñ ·(Q v) = 0

    (1)

  • RENORM/NORENORM
    Variables, {Q_i}, marked with the RENORM property obey the constraint:
    

    å i  Qi = 0

    (2)

  • CONSERVE/NOCONSERVE
    Variables marked with the CONSERVE property obey conservation laws ( e.g. momentum vs. velocity).
• LIBRARY name
  Specify a library requirement. Different flash modules require different external libraries and they must inform setup so it can link them into the executable. Valid library names are HDF4 and HDF5. Support for more libraries can be added by modifying the site-specific Makefile.h files to include appropriate makefile macros.
• FLUX name
  Register flux variable name with the framework.

#! SOME_PARAMETER The purpose of this parameter is whatever

#! SOME_PARAMETER The purpose of this parameter is whatever
#! &              This is a second line

Parameter comment lines are special because they are used by setup to build a formatted list of commented runtime parameters for a particular problem setup. This information is generated in the file paramFile.html in the $FLASH_HOME directory. A file paramFile.txt is also generated.

3.1.2  Interface layer and database module

After the module Config and Makefile are written, the source code files that carry out the specific work of the modules must be added. These source files can be separated into two broad categories: what we term ``wrapper functions'' and ``algorithms.'' In this section we discuss how to construct wrapper functions.

When constructing a FLASH module the designer must define a public interface of procedures that the module exposes to clients (ie other modules in an application). This is true regardless of the specific development language or syntactic features chosen to organize the procedures. These public functions are then defined in one or more source code files that form what we refer to as the interface layer.

Currently, there is no language-level formality in FLASH for enforcing the distinction between the public interface and private module functions that the interface harnesses. The developer is certainly encouraged to implement this within the chosen development language - static functions in C; private class functions in C++; private subroutines in a Fortran module, etc. However, nothing in FLASH will force this distinction and carry out the associated name-hiding within an application.

The most important aspect of the interface/algorithm distinction is related to the rules for data access. Wrapper functions communicate directly with the FLASH database module to access grid data (see below). However, algorithms are not permitted to query the database module directly. Instead, they must receive all data via a formal function argument list. Thus, when a module A wishes to request the services of module B, A calls one of B's public wrapper functions. Rather than being required to pass all necessary data to B through a procedure argument list, B may ``pull'' data it needs access to from the Database, marshal it as necessary, call the module algorithm(s), receive the updated data, and update the database.

The following subsections describe the methods provided by the database module in more detail.

3.1.2.1  dBaseGetData/dBasePutData   

Usage 

  call dBaseGetData([variable, [direction, [q1, [q2, [q3,]]]]] block, storage)
  call dBasePutData([variable, [direction, [q1, [q2, [q3,]]]]] block, storage)

Description 

Data exchange with grid variables. All variables registered with the framework via the VARIABLE keywords within a module Config file can be read/written with this pair of functions. The data is assumed to be composed of one or more structured blocks each with an integer block id = [1,num_blocks], where num_blocks is the total number of blocks on a given processor.

Example 

Given a Config file with the following variable registration specification:

  VARIABLE dens

  real, dimension(nxb,nyb,nzb) :: density
  do this_block = 1, total_blocks
    call dBaseGetData("dens", this_block, density)
    call foo(density)
  end do

Arguments 

Strings 

Integers 

Integers for variables and directions are not publicly available, but can be accessed through dBaseKeyNumber().

Note 

The ``xyzCube'' keyword is obsolete. When direction keyword is not specified, whole variable will be exchanged. When, in addition, variable keyword is omitted all variables for the block will be exchanged.

3.1.2.2  dBaseGetCoords/dBasePutCoords   

Usage 

  call dBaseGetCoords(variable, direction, [q,] block, storage)
  call dBasePutCoords(variable, direction, [q,] block, storage)

Description 

Access global coordinate information for a given block, including ghost points.

Example 

  real, DIMENSION(block_size) :: xCoords, data
  do i = 1, lnblocks
    call dBaseGetCoords("zn", "xCoord", i, xCoords)
    call dBaseGetData("dens", "xVector", 0, 0, i, data)
    call foo(data,xCoords)
  enddo

Arguments 

Strings 

Integers 

Integers for variables and directions are not publicly available, but can be accessed through dBaseKeyNumber().

3.1.2.3  dBaseGetBoundaryFluxes/dBasePutBoundaryFluxes   

Usage 

  call dBaseGetBoundaryFluxes(position, direction, block, storage)
  call dBasePutBoundaryFluxes(position, direction, block, storage)

Description 

Access boundary fluxes for all flux variables on a specified block associated with a given time-level. Currently, only the current and previous time-step are supported.

Example 

  real, dimension(nfluxes,ny,nz): xl_bound_fluxes, xr_bound_fluxes
  do this_block = 1, num_blocks
    call dBaseGetBoundaryFluxes(0,0,"xCoord", this_block, xl_bound_fluxes)
    call dBaseGetBoundaryFluxes(0,1,"xCoord", this_block, xr_bound_fluxes)
    call foo(xl_bound_fluxes, xr_bound_fluxes)
    call dBasePutBoundaryFluxes(0,0,"xCoord",this_block,xl_bound_fluxes)
    call dBasePutBoundaryFluxes(0,1,"xCoord",this_block,xl_bound_fluxes)
  end do

Arguments 

Strings 

Integers 

3.1.2.4  dBaseKeyNumber   

Usage 

  result = dBaseKeyNumber(keyname)

Description 

For faster performance, dBase{Get,Put}{Data,Coords} can be called with integer arguments instead of strings. Each of the string arguments accepted by the Get/Put methods can be replaced by a corresponding integer. However, these integers are not publicly available. To obtain them, one must call dBaseKeyNumber().

Example 

  integer :: idens, ixCoord
  idens   = dBaseKeyNumber("dens")
  ixCoord = dBaseKeyNumber("xCoord")
  call dBasePutData(idens, ixCoord, block_no, data)

Arguments and return type 

Strings 

3.1.2.5  dBaseSpecies    

Usage 

  result = dBaseSpecies(index)

Description 

Maps species number (from one to maximum number of species) to variable number (actual index in ``unk'' array).

Arguments and return type 

At present, all of the species are stored with adjacent indices in the solution array, thus one can find the index of the first isotope with a call to dBaseSpecies(1), and increment this value by 1 to get the next species.

3.1.2.6  dBaseVarName    

Usage 

  result = dBaseVarName(keynumber)

Description 

Given a key number, return the associated variable name.

Example 

  integer       :: idens
  char(len = 4) :: name
  idens = dBaseKeyNumber("dens")
  name  = dBaseVarName(idens)     ! name now = "dens"

Arguments and return type 

3.1.2.7  dBasePropertyInteger   

Usage 

  result = dBasePropertyInteger(property)

Description 

Accessor methods for integer-valued scalar variables.

Arguments and return type 

Strings 

3.1.2.8  dBasePropertyReal   

Usage 

  result = dBasePropertyReal(property)

Description 

Accessor methods for real-valued scalar variables.

Arguments and return type 

Strings 

3.1.2.9  dBaseSetProperty   

Usage 

  call dBaseSetProperty(property, value)

Description

Mutator methods for writeable scalar variables. See dBasePropertyInteger/Real for documentation on property names.

Arguments

Strings

3.1.2.10  dBaseVarIndex   

Usage

  result = dBaseVarIndex(property)

Description

Returns pointer to array of integer keys of variables with specified property.

Arguments and return type

Strings

3.1.2.11  Various pointer-returning functions   

Each of these functions allows FLASH developers to hook directly into an internal data structure in the database. In general these functions will offer better performance then their corresponding dBaseGet/Put counterparts, and will require less memory overhead. However, the interfaces are more complicated and the functions are less flexible, and less safe, so it is suggested that developers strongly consider using dBaseGet/PutData when performance differences are small.

Each function returns a Fortran 90 pointer to the solution vector on the specified block. If no block is specified, a pointer is returned to all blocks on the calling processor. Currently the array index layout is assumed to be (var,nx,ny,nz,block) in row-major ordering. The scratch (unksm) array stores variables with no guardcells; this name should probably be changed in the future.

Argument and return type

3.1.2.12  dBaseGetDataPtrAllBlocks()   

Return an F90 pointer to the left-hand-side solution vector for all blocks on a given processor, arranged in row-major order as: (var,nx,ny,nz,block). dBaseKeyNumber must still be called to access the elements of the array.

3.1.2.13  dBaseGetDataPtrSingleBlock(block_no)   

Return an F90 pointer to the left-hand-side solution vector on a specified block, arranged as (var,nx,ny,nz).

3.1.2.14  dBaseGetPtrToXCoords()   

Return an F90 pointer to an array containing information on the x-coordinates of the AMR blocks. The array returned is arranged as (block_position, i, block_number), where block_position values denote center, left, or right coordinates, and are obtained by calling dBaseKeyNumber with ``zn'', ``znl'', ``znr'' and using the corresponding index to access the approriate row in the array. For example:

  real, pointer, dimension(:,:,:) :: xCoords
  real                            :: x, xl, xr
  integer                         :: izn, iznl, iznr
  izn  = dBaseKeyNumber("izn")
  iznl = dBaseKeyNumber("iznl")
  iznr = dBaseKeyNumber("iznr")
  xCoord => dBaseGetPtrToXCoords()
  do this_block = 1, num_blocks
    do i = 1, blocksize
      x  = xCoord(izn, i, this_block)     ! get first center coord
      xl = xCoord(iznl, i, this_block)    ! get first left coord
      xr = xCoord(iznr, i, this_block)    ! get first right coord
      ...
    enddo
  enddo

3.1.2.15  dBaseGetPtrToYCoords()   

See dBaseGetPtrToXCoords().

3.1.2.16  dBaseGetPtrToZCoords()   

See dBaseGetPtrToXCoords().

3.1.2.17  dBaseGetScratchPtrAllBlocks()   

Return an F90 pointer to a scratch array of size (2, nxb, nyb, nzb, maxblocks).

3.1.2.18  dBaseGetScratchPtrSingleBlock(block_no)   

Return an F90 pointer to a scratch array of size (2, nxb, nyb, nzb).

3.1.2.19  AMR tree interface functions    These functions enable FLASH developers to directly access the data structures used by PARAMESH to describe the adaptive mesh. In general they should not be needed by developers of physics modules. Also, they may not be available in future versions of FLASH.

Arguments and return types

dBaseNeighborBlockList (block)

Given a block ID, return an array of block ID's which are the neighbors of the specified block. The returned array is of size max_faces = 6, but not all of the six elements will have meaningful values if the problem is run in fewer than three dimensions. Assuming the function is called as NEIGH = dBaseNeighborBlockList(), the ordering is as follows. The neighbor on the lower x face of block L is at NEIGH(1,L), the neighbor on the upper x face at NEIGH(2,L), the lower y face at NEIGH(3,L), the upper y face at NEIGH(4,L), the lower z face at NEIGH(5,L), and the upper z face at NEIGH(6,L). If any of these values are set to -1 or lower, there is no neighbor to this block at its refinement level. However there may be a neighbor to this block's parent. If the value is -20 or lower then this face represents an external boundary, and the user is required to apply some boundary condition on this face. The exact value below -20 can be used to distinguish between the different boundary conditions which the user may wish to implement.

dbaseNeighborBlockProcList (block)

Given a block ID, return an array of size max_faces = 6 elements containing processor ID's identifying the processor that a given neighbor resides on. Ordering is identical to dBaseNeighborBlockList().

dBaseChildBlockList (block)

Given a block ID, return an array of size max_child = 2 * max_dim elements containing the block ID's of the child blocks of the specified block. The children of a parent are numbered according to the Fortran array ordering convention, ie. child 1 is at the lower x, y, and z corner of the parent, child 2 at the higher x coordinate but lower y and z, child 3 at lower x, higher y and lower z, child 4 at higher x and y and lower z, and so on.

dBaseChildBlockProcList (block)

Given a block ID, return an array of size max_child elements containing processor ID's of the children of the pecified block. Ordering is identical to dBaseChildBlockList().

dBaseParentBlockList (block)

Given a block ID, return the ID of the block's parent block.

dBaseParentBlockProcList (block)

Given a block ID, return the processor ID upon which the block's parent resides.

dBaseRefinementLevel (block)

Given a block ID, return that block's integer level of refinement.

dBaseNodeType (block)

Given a block ID, return the block's node type. If 1 then the node is a leaf node, if 2 then the node is a parent but with at least 1 leaf child, otherwise it is set to 3 and it does not have any up-to-date data.

dbaseNeighborType (block)

Given a block ID, return an array of size (mfaces, maxblocks_tr), containing the type ID's of the neighbors of the specified block. mfaces = mdim * 2, where mdim is the maximum possible dimensionality (3).

dBaseBlockCoord (block)

Given a block ID, return an array of size ndim containing the x,y,z coordinates of the center of the block.

dBaseBlockSize (block)

Given a block ID, return an array of size ndim containing the block size in the x, y, and z directions.

dBaseRefine (block)

Given a block ID, return .true. if that block is set for refinement in the next call to amr_refine_derefine(), and .false. otherwise.

dBaseDeRefine (block)

Given a block ID, return .true. if that block is set for refinement in the next call to amr_refine_derefine(), and .false. otherwise.

3.1.3  Algorithms

Within each module are one or more procedures which perform the bulk of the computational work for the module. A principal strategy behind the FLASH architecture is to decouple these procedures as much as possible from the details of the framework in which they are embedded. This is accomplished by requiring that all module algorithms communicate data only through function argument lists. That is, algorithms may not query the database directly nor may they depend on the existence of externally defined or global variables. This design ensures that algorithms can be tested, developed, and interchanged in complete isolation from the larger, more complicated framework.

Thus, each algorithm in a module should have a well defined argument list. It is up to the algorithm developer to make this as general or restrictive as he/she sees fit. However, it is important to keep in mind that, the more rigid the argument list, the less chance that another algorithm can share its interface. The consequence is that the developer would have to add an entirely new wrapper function for just slightly different functionality.

3.2  The FLASH source tree

An abstract representation of the FLASH architecture appears in Figure 3. Each box in this figure represents a component (FLASH module), which publishes a small set of public methods to its clients. These public methods are expressed through virtual function definitions (stubs under Fortran 90), which are implemented by real functions supplied by sub-modules. Typically each component represents a different class of solver. For time-dependent problems, the driver uses time-splitting techniques to compose the different solvers. The solvers are divided into different classes on the basis of their ability to be composed in this fashion and upon natural differences in solution method (e.g., hyperbolic solvers for hydrodynamics, elliptic solvers for radiation and gravity, ODE solvers for source etc.).

The adaptive mesh refinement module is treated in the same way as the solvers. The means by which the driver shares data with the solver objects is the primary way in which the architecture affects the overall performance of the code. Choices here, in order of decreasing performance and increasing flexibility, include global memory, argument-passing, and messaging. FLASH 2.0 has eliminated global variable access in favor of a well-defined set of accessor and mutator methods managed by the centralized database module. When done with an eye toward optimization, the effects on performance are tolerable, and the benefits for maintainability and extensibility are significant. This is discussed in greater detail below.

3.2.1  Code infrastructure

The structure of the FLASH source tree reflects the module structure of the code, as shown in Figure 4 for an older version of FLASH. The general plan is that source code is organized into one set of directories, while the code is built in a separate directory using links to the appropriate source files. The links are created by a source configuration script called setup, which makes the links using options selected by the user and then creates an appropriate makefile in the build directory. The user then builds the executable with a single invocation of gmake.

Source code for each of the different code modules is stored in subdirectories under source/. The code modules implement different physics, such as hydrodynamics, nuclear burning, and gravity, or different major program components, such as the main driver code and the input/output code. Each module directory contains source code files, makefile fragments indicating how to build the module, and a configuration file (see Section 6.1.2).

Each module subdirectory may also have additional sub-module directories underneath it. These contain code, makefiles, and configuration files specific to different variants of the module. For example, the hydro/ module directory can contain files which are generic to hydrodynamical solvers, while its explicit/ subdirectory contains files specific to explicit hydro schemes and its implicit/ subdirectory contains files specific to implicit solvers. Configuration files for other modules which need hydrodynamics can specify hydro as a requirement without mentioning a specific solver; the user can then choose one solver or the other when building the code (via the modules file; Section 6.1.3). When setup configures the source tree it treats each sub-module as inheriting all of the source code, configuration files, and makefiles in its parent module's directory, so generic code does not have to be duplicated. Sub-modules can themselves have sub-modules, so for example one might have hydro/explicit/split/ppm and hydro/implicit/ppm. Source files at a given level of the directory hierarchy override files with the same name at higher levels, whereas makefiles and configuration files are cumulative. This permits modules to supply stub routines that are treated as `virtual functions' to be overridden by specific sub-modules, and it permits sub-module directories to be self-contained.

When a module is not explicitly included by Modules, only one thing is done differently by setup: sub-modules are not included, except for a null sub-module, if it is present. Most top-level modules should contain only stub files to be overridden by sub-modules, so this behavior allows the module to be `not included' without extra machinery (such as the special stub files and makefiles required by earlier versions of FLASH). In those cases in which the module's files are not appropriate for the `not included' case, the null sub-module allows one to override them with appropriate versions.

New solvers and new physics can be added. At the current stage of development of FLASH it is probably best to consult the authors of FLASH (see Section 7) for assistance in this. Some general guidelines for adding solvers to FLASH 2.0 may be found in Section 17.

The setups/ directory has a structure similar to that of source/. In this case, however, each of the "modules" represents a different initial model or problem, and the problems are mutually exclusive; only one is included in the code at a time. Also, the problem directories have no equivalent of sub-modules. A makefile fragment specific to a problem need not be present, but if it is, it is called Makefile. Section 16 describes how to add new problem directories.

The setup script creates a directory called object/ in which the executable is built. In this directory setup creates links to all of the source code files found in the specified module and sub-module directories as well as the specified problem directory. (A source code file has the extension .c, .C, .f, .f90, .F90, .F, .fh, or .h.) Because the problem setup directory and the machine-dependent directory are scanned last, links to files in these directories override the ``defaults'' taken from the source/ tree. Hence special variants of routines needed for a particular problem can be used in place of the standard versions by simply giving the files containing them the same names.

Using information from the configuration files in the specified module and problem directories, setup creates a file named init_global_parms.F90 to parse the runtime parameter file and initialize the runtime parameter database. It also creates a file named rt_parms.txt, which concatenates all of the PARAMETER statements found in the appropriate configuration files and so can be used as a ``master list'' of all of the runtime parameters available to the executable.

setup also creates makefiles in object/ for each of the included modules. Each copy is named Makefile.module, where module is driver, hydro, gravity, and so forth. Each of these files is constructed by concatenating the makefiles found in each included module path. So, for example, including hydro/explicit/split/ppm causes Makefile.hydro to be generated from files named Makefile in hydro/, hydro/explicit/, hydro/explicit/split/, and hydro/explicit/split/ppm/. If the module is not explicitly included, then only hydro/Makefile is used, under the assumption that the subroutines at this level are to be used when the module is not included. The setup script creates a master makefile (Makefile) in object/ which includes all of the different modules' makefile fragments together with the site- or operating system-dependent Makefile.h.

The master Makefile created by setup creates a temporary subroutine, buildstamp.F90, which echoes the date, time, and location of the build to the FLASH log file when FLASH is run. To ensure that this subroutine is regenerated each time the executable is linked, the Makefile deletes buildstamp.F90 immediately after compiling it.

The setup script can be run with the -portable option to create a directory with real files which can be collected together with tar and moved elsewhere for building. In this case the build directory is assigned the name object_problem/. Further information on the options available with setup may be found in Section 13. .

Additional directories included with FLASH are tools/, which contains tools for working with FLASH and its output (Sec:FLASH output comparison utility), and docs/, which contains documentation for FLASH (including this user's guide) and the PARAMESH library.

3.3  Modules included with FLASH: a brief overview

The current FLASH distribution comes with a set of core components that form the backbone of many common problems, namely: database, driver, hydro, io, mesh, particles, source_terms, gravity, and materials. A detailed discussion of the role of each of these modules is presented in Part IV. Here we give a brief overview of each.

Section 4 describes in detail the various driver modules which may be implemented with FLASH. In addition to the default driver, which controls the initialization, evolution, and output of a FLASH simulation, four new driver modules have been written to implement different explicit time advancement algorithms. Three are written in the delta formulation: euler1, rk3, and strang_delta. The fourth, strang_state, is written in the state-vector formulation. A subsection concerning simulation services, runtime parameters, and logfiles is also included in this section.

Section 5 describe the FLASH I/O modules, which control how FLASH data structures are stored on different platforms and in different formats. Discussed in this section are the main I/O module hdf4, which uses the Hierarchical Data Format (HDF) for storing simulation data, two major HDF5 modules (serial and parallel versions), and Fortran 77 (f77) modules.

Section 6 describes the mesh module, together with the PARAMESH package of subroutines for the parallelization and adaptive mesh refinement (AMR) portion of FLASH.

Section 7 describes the two hydrodynamic modules included in FLASH 2.x. The first is based on the PROMETHEUS code (Fryxell, Müller, and Arnett 1989); the second is based on Kurganov numerical methods.

Section 8 describes the magnetohydrodynamics module included with the FLASH code, which solves the equations of ideal MHD.

Section 9 discusses the material properties module, which handles the tracking of multiple fluids in FLASH simulations. It includes the equation of state module, which implements the EOS for the hydrodynamical and nuclear burning solvers; the composition submodule, which sets up the different compositions needed by FLASH; and the stellar conductivity module, which may be used for computing the opacity of stellar material. Section describes source terms, including the nuclear burning module, which calculates the nuclear burning rate of a hydrodynamical simulation, and the stirring module, which adds a divergence-free, time-correlated `stirring' velocity at selected modes in a given hydrodynamical simulation.

Section 11 describes the gravity module, which computes gravitational potential or gravitational acceleration source terms for the code. It includes several sub-modules: the constant submodule, the plane parallel sub-module, the ptmass submodule, and the Poisson submodules.

Part 3
MODULES

4  Driver modules

The driver module controls the initialization, evolution, and output of a FLASH simulation. Initialization can be from scratch or from a stored checkpoint file. Evolution can use any of several different operator-splitting techniques to combine different physics operators and integrate them in time, or call a single operator if the problem of interest is not time-dependent. Output involves the production of checkpoint files, plot files, analysis data, and log file time stamps. In addition to these functions, the driver supplies important simulation services to the rest of the FLASH framework, including Fortran modules to handle runtime parameters, physical constants, memory usage reports, and log file management (these are discussed further in Section ).

The initialization and termination routines and simulation services modules are common to both time-dependent and time-independent drivers and thus are included at the highest level of driver. The file flash.F90 contains the main FLASH program (equivalent to main() in C) and calls these routines as needed. The default flash.F90 is empty and is intended to be overridden by submodules of driver. At this time only time-dependent drivers are supplied with FLASH; these are submodules of the driver/time_dep module. The time_dep version of flash.F90 calls the FLASH initialization routine, loops over timesteps, and then calls the FLASH termination routine. During the time loop it computes new timesteps, calls an evolution routine (evolve()), and calls output routines as necessary.

The details of each available time integration method are completely determined by the version of evolve() supplied by that method. The default time update method is to call each physics module's update routine for two equal timesteps - thus, hydro, source terms, gravity, hydro, source terms, gravity. The hydrodynamics update routines take a ``sweep order'' argument in case they are directionally split; in this case, the first call uses the ordering x-y-z, and the second call uses z-y-x. Each of the update routines is assumed to directly modify the solution variables. At the end of each pair of timesteps, the condition for updating the mesh refinement pattern is tested, and a refinement update is carried out if required.

The alternative ``delta formulation'' drivers (driver/time_dep/delta_form) modify a set of variables containing the change in the solution during the timestep. The change is only applied to the solution variables after all operators have been invoked. This technique permits more general time integration methods, such as Runge-Kutta methods, to be employed, and it provides a more flexible method for composing operators. However, only a few physics modules can make use of it as yet. More details on the delta formulation drivers appear in Section 4.1.

The driver module supplies certain runtime parameters regardless of which type of driver is chosen. These are described in Table 1.

4.1  New driver modules euler1, rk3, and strang_delta

New driver modules have been written to implement different explicit time advancement algorithms. This usage of the driver modules is slightly different than that of the default driver module, which does not directly implement a time advancement algorithm; the default driver and hydro modules each implement parts of the Strang splitting time advancement. First, a listing of the time advancement tasks common to all of the new drivers will be given. In following subsections, the details of each time advancement method will be described.

The three driver modules written in the delta formulation are euler1, rk3, and strang_delta. They make appropriate calls to the physics modules, and update the solution by calling functions provided by the formulation module. The strang_state driver is written in the state-vector formulation; it also calls the physics modules, but does not update the solution. To use the new modules, first choose the driver by including one of the following lines into the Modules file.

/driver/time_dep/delta_form/euler1
/driver/time_dep/delta_form/rk3
/driver/time_dep/delta_form/strang_delta
/driver/time_dep/delta_form/strang_state

/formulation/delta_form

/formulation

For delta formulation time advancements, the new driver modules use the formulation module to hold and access the global DU array, and to update the solution. In the state-vector formulation, the formulation module is not directly used by the driver; instead, the physics modules call the update subroutines the formulation module provides.

The new driver modules discretize the left-hand side of

¶V ¶t = (spatial difference terms)+ (source terms).

(3)

The time advancement algorithm is contained in a subroutine named evolve. Each call to evolve advances the solution through one time step, Dt. Each time advancement algorithm begins with a vector of primary variables, Vⁿ at time tⁿ, and an associated set of auxiliary variables, Wⁿ. Primarily through calls to other physics modules, evolve applies a set of operations to the variables to produce an updated vector, Vⁿ⁺¹ at time t_n+1 = t_n + Dt. Depending on the formulation, the time advancement may or may not update the auxiliary variables - in the state-vector formulation, the other physics modules update them.

The distinction between V and W is that time-dependent differential equations are solved to determine the primary variables. The auxiliary variables are obtained from the primary variables through algebraic relations. Often the primary variables are the conserved variables, U, and in the rest of this section U will replace V. However, the time advancement algorithms implemented do not require this correspondence.

The time advancement algorithms are written generally, in that each differential equation is treated in the same way. The distinction between the equations (for example, between the x-momentum equation and the total energy equation) is expressed in the other physics modules. The time advancement algorithm does not need to know the identity of the variables on which it operates, except possibly to update the auxiliary variables from the primary variables; but this update is handled by a call to a subroutine provided by the formulation module.

4.1.1  The euler1 module

The euler1 module implements the first-order, Euler explicit scheme:

Un+1 = Un + Dt L(Un)

(4)

where L(U) represents all of the physics modules. The euler explicit method is implemented in the delta formulation. No runtime parameters are defined for this module.

At the beginning of a time step, DU is set to zero. Each of the physics modules is called, with Uⁿ as the initial state, and adds its contribution to DU. After all the physics modules have been called, the global DU array holds L(Uⁿ). Equation (4) yields Uⁿ⁺¹. Finally, the auxiliary variables are updated from the conserved variables with a call to the global update subroutine provided by the formulation module.

Note that because all the physics modules start with the same initial state, the order in which the physics modules are called does not affect the results (except possibly through floating point roundoff differences when contributing to DU.)

The set of steps, consisting of calls to physics modules, updating the conserved variables, and updating the auxiliary variables, is often called a stage. The majority of the computational cost of a stage is the in the calls to the other physics modules; this component corresponds to a ``function evaluation'' for ordinary differential equation solvers. In the Euler explicit algorithm, there is one stage per time step.

4.1.2  The rk3 module

Runge-Kutta schemes are a class ordinary differential equation solvers which are appreciated for their higher order of accuracy, ease of implementation, and relatively low storage requirements. There are many third-order Runge-Kutta methods; all require at least three stages. Most require at least three storage locations per primary variable, but the one implemented in the delta formulation in FLASH, derived by Williamson (J. Comp. Phys. 35:48, 1980) requires only two.

The two storage registers will be referred to as U and DU. The global solution vector, U, holds Uⁿ at the beginning of the time step, then intermediate solutions U^(·) at the end of each stage. The manipulation of the global DU array is more complicated. DU accumulates contributions from the physics modules during a stage, but it also holds results from previous stages; it is important to distinguish between the results of the physics modules, L(U), and the quantity held in the global DU array. First the algorithm will be shown; then the usage of the global DU array will be discussed. For the equation

¶U ¶t = L(U),

(5)

Williamson's algorithm is, starting with Uⁿ,

U(1) = Un +  1 3 Dt é ë L(Un) ù û (6) U(2) = U(1) +  15 16 Dt é ë L(U(1)) -  5 9 L(Un) ù û (7) Un+1 = U(3) = U(2) +  8 15 Dt é ë L(U(2)) -  153 128 æ è L(U(1)) -  5 9 L(Un) ö ø ù û . (8)

U^(m) is the result of the m-th stage, and the auxiliary variables are updated each time a new U^(m) is computed.

The algorithm is implemented using the following steps to attain the low storage. At the beginning of the time step, DU is set to zero. During the first stage each physics module contributes to DU, so after all have contributed, DU holds the bracketed term in eq. (6), L(Uⁿ). U⁽¹⁾ is then computed using eq. (6). Stage 1 is completed by multiplying DU by -5/9, which is required for the following stages. The process is repeated for stage 2: after the physics modules have contributed, DU holds the bracketed term in eq. (7); U⁽²⁾ is computed by eq. (7) and stored in U; then DU is multiplied by -153/128. Stage 3 is similar, but ends after U⁽³⁾ = Uⁿ⁺¹ is computed and stored. It is critical that the only changes made to the DU array are those just listed; no physics module should change the value of DU, except to add its contribution, and since DU holds information from previous stages, it should not be reset to zero except at the beginning of the time step.

No runtime parameters are defined for this module.

4.1.3  strang_state and strang_delta modules

The second-order accurate splitting method (Strang 1968) is attractive because of its low memory requirements. The algorithm is based on the operator splitting approach, in which a set of subproblems is solved rather a single complicated problem. Each subproblem typically accounts for one term in a system of partial differential equations, representing a particular type of physics and for which an appropriate (specialized) numerical method is available. The basic operator splitting method is first-order accurate, but the Strang splitting scheme is second-order accurate over two time steps. In the first time step, the subproblems are solved in a given sequence. Second-order accuracy is obtained by reversing the sequence in the second time step.

A key feature of the operator splitting approach is that the output of one subproblem is the input to the next subproblem. This allows an implementation that, globally, stores only th