Writing An AstroBEAR Module

AstroBEAR problem modules are stored in the astrobear/contrib directory, typically under another subdirectory. The U. of R. group, for instance, keeps their work in astrobear/contrib/astro.

Throughout this page we will be using the following convention: the new module you are creating is called NewModule, and the associated file is new_module.f90. We will be also be assuming that the user is working in the $(ASTROBEAR)/contrib/astro directory.

Creating a New Module In AstroBEAR

This creates a module file and makes sure that AstroBEAR can link to it. The actual writing of the module and defining its interfaces is a separate process described below.

Note that these instructions assume that you are in contrib/astro.

1. Create a new Directory NewModule in contrib/astro.

2. Create a file NewModule/new_module.f90. This is the Fortran 90 file where your code will be.

3. Declare a new logical variable lModuleName in globaldeclarations.f90.

4. Add your new flag lModuleName to the ProblemData namelist in globaldeclarations.f90.

5. Add the line USE ModuleName to the list of USE statements in problem.f90.

6. In Makefile, look under BASE_PROBLEM_SOURCES=\ for a group of statements of the form $(BEARCLAW_ASTRO)/ModuleName/ModuleName.f90 \ and create a similar statement for NewModule.

7. Also in Makefile, look under BASE_PROBLEM_OBJECTS=\ for a group of statements of the form $(BEARCLAW_ASTRO)/ModuleName/ModuleName.o \ and create a similar statement for NewModule.

Integrating your Module Into AstroBEAR

Each problem module contains four main subroutines that are required in order for the module to interface with the the AMR engine. These subroutines are listed below, along with the "hooks" needed to link these routines into AstroBEAR.

The basics of writing these modules will be explained below under "Module Basics."

• scaleNewModule(): The scaling subroutine. scaleNewModule is where computational units are obtained from physical units. It is also where AstroBEAR reads in data from problem-specific input files.
  • hook: in i_setprob.f90::SetProb SetProb(), look for a group of statements of the form IF (lModuleName) CALL scaleModuleName and add a similar line for your module.

• qinitNewModule(Info): This is where you initialize the data arrays. Remember that this routine is getting called on a grid-by-grid basis, so attempts to initialize data outside this grid will probably cause a segfault.
  • hook: in i_qinit.f90::qinit(), look for a group of statements of the form IF (lModuleName) CALL PlaceModuleName(Info) and add a similar line for your module.

• b4stepNewModule(Info): The before-step subroutine. This procedure is called before each time-step of the simulation, so that the driving mechanisms specified by the user can be reapplied at each step. For example, a jet simulation would inject new material into the system during b4step(). Note that no integration occurs during this step; this is just where new conditions are introduced (and sometimes renewed).
  • hook: in problem.f90::b4step(), look for a group of statements of the form IF (lModuleName) CALL b4stepModuleName(Info) and add a similar line for your module.

• problemErrFlagNewModule(Info): Flags regions where this module requires additional resolution. If you do not have any special refinement needs, then just leave this routine as a stub.
  • hook: in problem.f90::problemErrFlag(), look for a group of statements of the form IF (lModuleName) CALL problemErrFlagModuleName(Info) and add a similar line for your module.

AstroBEAR Module Basics

Simulation Data

All AstroBEAR modules have at one thing in common: initializing the problem domain. Within our code, the problem domain's data is held in NodeInfo structures, which is why so many module subroutine take a NodeInfo structure as a parameter.

There are two major data arrays in NodeInfo: the q array and the aux array. q holds the cell-centered data and is used by all AstroBEAR simulations. The aux array holds face-centered and edge-centered data, and is only used in MHD problems. If you are running a strictly hydrodynamic problem or a hydrodynamic + elliptic problem, then you will not need aux.

Dimensions

q and aux are both 5D arrays, with the first four dimensions giving the coordinates and the fifth being a variable. The 4D option is an old BearCLAW legacy, but AstroBEAR is not set up to run relativistic problems. Consequently, the 4D index should always be 1:

Info%q(x, y, z, 1, variable)

Currently, AstroBEAR can only run 2D and 3D problems, but a 1D hydro or MHD problem can be simulated by defining a very narrow 2D problem domain and then making sure all the activity is defined in the x-direction (i.e., no py or pz components).

Units and Scaling

Astrophysical problems involve many different physical units and constants with a wide range of scales. To avoid the loss of precision that comes when computers try to work with, say, a 10^-8 variable and a 10²⁴ constant in the same expression, we scale our units into computational units before storing them in the data arrays.

Usually, the physical scales are defined in the physics.data file—you simply enter the scales for density, temperature, velocity, etc in that file, and AstroBEAR will read them in. More complicated scaling would be defined in the scaleNewModule() routine (see below).

You have two options for making sure that you only put scaled quantities in the data arrays: you can scale your input values before you enter them into your input file (and then assume that you are reading in scaled quantities), or you can use physical quantities in your input files and then scale them within your problem module:

scaled_qty = physical_qty / physical_scale

Either way, a good sanity check is to print out the physical quantities your program uses after the problem is set up. This verifies that the values you think are going in are the values that are actually getting used.

Initializing a Grid

Initializing a grid involves taking a spatially-constructed problem setup and discretizing it so that it fits nicely in an array. This process is easiest to explain by dissecting an example, such as the one below, where we are trying to initialize the grid with a density distribution given by rho(x,y,z):

    q=>Info%q

    rmbc=Info%r *Info%mbc

    mx=Info%mX(1); dx=Info%dX(1); xlower=Info%Xlower(1)
    my=Info%mX(2); dy=Info%dX(2); ylower=Info%Xlower(2)
    mz=Info%mX(3); dz=Info%dX(3); zlower=Info%Xlower(3)

    SELECT CASE(nDim)
    CASE(2)
       zrmbc=0;mz=1;zl=0;dz=0
    CASE(3)
       zrmbc=rmbc
    END SELECT

    ! Initialize environment
    DO k=1-zrmbc, mz+zrmbc
        DO j=1-rmbc, my+rmbc
            DO i=1-rmbc, mx+rmbc

                x=(xlower + (REAL(i,xPrec)-half) * dx)
                y=(ylower + (REAL(j,xPrec)-half) * dy)
                z=(zlower + (REAL(k,xPrec)-half) * dz)

                q(i,j,k) = rho(x,y,z)
            END DO
        END DO
    END DO

Each grid (or NodeInfo structure) comes with several arrays that describe its geometry:

• mX(4): The grid's overall size. Each element mx(n) represents the number of cells along the nth dimension.
• Xlower(4): The spatial coordinates of the grid's lower bound in each dimension. These values are always given in computational units.
• dX(4): The size of a spatial step in a given dimension. dX(n) can also be thought of as the size of a cell along dimension n.

To make this example more readable (and therefore easier to debug), we assigned the various array values more Cartesian-sounding names (ie., Info%dX(2) = dy).

Similarly, when calculating the spatial equivalents of an index, we went with x, y and z. The half-step in the position calculations represents the fact that the spatial position is at the center of the cell, not the left side.

The check on the number of dimensions makes sure that 2D problem does not accidentally get initialized with ghost cells (see below) or a spatial step. This allows us to use the same initialization code for 2D or 3D problems.

You may have noticed the nested loops don't go from 1 to mX(n). This is because the data array size is not the same as the grid size. The size of each grid along dimension n is mX(n). The data arrays, however, have a number of "ghost cells" associated with them. The number of ghost cells on the end of a grid is given by the quantity

rmbc = Info%r * Info%mbc

These ghost cells only appear along the dimensions of the problem, though, so the data arrays for a 2D problem will not have ghost cells along the third dimension. So the real extents of the data arrays in a 3D problem are:

Info%q(1 - rmbc : mX(1) + rmbc, &
       1 - rmbc : mX(2) + rmbc, &
       1 - rmbc : mX(2) + rmbc, &
       1, &
       NrVars)

Info%aux(1 - rmbc : mX(1) + rmbc + 1, &
         1 - rmbc : mX(2) + rmbc + 1, &
         1 - rmbc : mX(2) + rmbc + 1, &
         1, &
         NrVars)

Flagging Cells for Refinement

Some modules may need specific regions refined, regardless of whether or not there is any obvious error there. AstroBEAR flags cells for refinement using the array

Info%ErrorFlags(1 - rmbc : mX(1) + rmbc, &
                1 - rmbc : mX(2) + rmbc, &
                1 - rmbc : mX(2) + rmbc, &
                1, &
                NrVars)

To clear the cell at (i,j,k), simply set Info%ErrorFlags(i,j,k) to 0. An error flag of 0 means that the cell does not need to be refined (although it might happen anyway). To mark a cell for refinement, set Info%ErrorFlags(i,j,k) to 1. This cell will be refined the next time NewSubGrids() is called.