Writing Modules in AstroBEAR 2.0

AstroBEAR problem modules are stored in the modules directory. Each problem module gets its own directory; when compiling the code, the user creates a symbolic link Problem to the appropriate directory. For more information about setting up a problem directory, see Setting up and Compiling A Problem for more details.

The contents of the problem directory are up to the author, but at least one file must be present: problem.f90. This is the source file that contains the routines referenced by module_control.f90. The tutorial below is basically a discussion of how to fill this module out.

When a problem directory is checked into AstroBEAR, it is usually checked in with a number of .data files. From the standpoint of this tutorial, the most important is problem.data. This is where all of the user-defined variables are stored; these are typically read in by the ProblemModuleInit() routine in problem.f90.

Throughout this page we will be using the following convention: the new module file you are creating is called problem.f90. We will be also be assuming that the user is working in the $(ASTROBEAR)/modules 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 modules.

1. Create a new Directory NewModule in modules.

2. Create a symbolic link to this file with the command ln -s NewModule Problem. Note, that if the Problem symbolic link already exists, you may have to first remove it and then make it anew.

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

Integrating your Module Into AstroBEAR

Each problem module contains five main subroutines that are required in order for the module to interface with the the AMR engine. These subroutines are listed below.

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

• ProblemModuleInit(): Module variables are initialized here. This routine is where problem.data are read in and module-level namelists are populated. ProblemModuleInit() is also a popular place to put sanity checks that verify the correctness of all inputs. ProblemModuleInit() is also a common place to initialize source terms.

• ProblemGridInit(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. For an example routine see ProblemGridInit

• ProblemBeforeStep(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 ProblemBeforeStep(). Note that no integration occurs during this step; this is just where new conditions are introduced (and sometimes renewed). If you have no special pre-step needs, then leave this routine as a stub.

• ProblemAfterStep(Info): The after-step subroutine. This procedure is called following each time step of a system. This is perhaps the least commonly used of the control subroutines, but divergence cleaning, special output files, and other post-processing operations could be performed here. If you have no special post-step instructions, then just leave this routine as a stub.

• ProblemSetErrFlag(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.

AstroBEAR Module Basics

The best way to understand how a module works is by example. The following example is a module for the Rayleigh-Taylor instability.

Place sample module here

Simulation Data

All AstroBEAR modules have at least one thing in common: initializing the problem domain. Within our code, the problem domain's data is held in InfoDef structures, which is why so many module subroutines take an InfoDef structure as a parameter. To make use of the InfoDef structure this statement is required at the beginning of problem.f90:

USE DataDeclarations

There are two major data arrays in InfoDef: the q array and the aux array. q holds the cell-centered data and is used by all AstroBEAR simulations. The q array takes the form q(x,y,z,variable) where variable is itself a 1D array that holds the physical quantities such as density, momentum, energy, etc. The order of the quantities in the variable array is (rho, px, py, pz, E). The aux array holds face-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

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).

The core region of the Info%q array (which does not include ghost zones) is a 1:mx by 1:my by 1:mz box. Mx, my, and mz denote the number of cells in the x, y, and z directions, respectively. In two dimensions, mz = 1, reducing the box to a rectangle. Info%q is cell-centered, so the values are assumed to be taken from the midpoint of the cell. Thus, the cell-to-space conversion is:

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

The Info%aux array is a little different. The aux array holds magnetic flux values, which are face-centered. This means that every cell-centered value in Info%q is bracked in each dimension by two Info%aux values. To accommodate the extra values, Info%aux is a 1:mx+1 by 1:my+1 by 1:mz+1 box, but the aux dimensions are actually different for each variable:

Bx = Info%aux(1:mx+1, 1:my, 1:mz, 1)
By = Info%aux(1:mx, 1:my+1, 1:mz, 2)
Bz = Info%aux(1:mx, 1:my, 1:mz+1, 3)

The additional cells (the ones in the "upper-front right" corner of the aux array) are not used.

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 ProblemModuleInit() routine (see above).

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):

Note that during the ProblemGridInit routine, ghost zones do not need to be initialized (rmbc = 0) - however - during beforestep calculations they should be

    q=>Info%q

    rmbc=levels(Info%level)%gmbc(levels(Info%level)%step)

    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 InfoDef structure) comes with several arrays that describe its geometry:

• mX(3): The grid's overall size. Each element mx(n) represents the number of cells along the nth dimension.
• Xlower(3): The spatial coordinates of the grid's lower bound in each dimension. These values are always given in computational units.
• dX(3): 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.

The levels() array is a global array structures that contain information specific to each level. You can see in this example that we reference the appropriate level by using the Info%level attribute.

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

levels(Info%level)%CoarsenRatio * levels(Info%level)%gmbc

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(3) + rmbc, &
       NrVars)

Info%aux(1 - rmbc : mX(1) + rmbc + 1, &
         1 - rmbc : mX(2) + rmbc + 1, &
         1 - rmbc : mX(3) + rmbc + 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, &
                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. The best place to do this is in the ProblemSetErrFlags() routine; most conventional physical criteria for refinement are already handled by AstroBEAR itself.

Additional Physics

AstroBEAR supports hydrodynamics and AMR by default. Other physical processes such as magnetic fields and source terms require extra overhead, so they must be enabled by the user.

MHD

Enabling MHD requires two changes to the data files:

• In global.data, set the MaintainAuxArrays flag to T.
• In physics.data, set lMHD to T.

You will also need to initialize the aux arrays in your module file. In fact, magnetic fields are often initialized by setting the aux values first, and then calculating the cell-centered magnetic fields in q by averaging the aux values on either side of a cell:

Info%q(i,j,k,iBx) = half * (Info%aux(i,j,k,1) + Info%aux(i+1,j,k,1))
Info%q(i,j,k,iBy) = half * (Info%aux(i,j,k,1) + Info%aux(i,j+1,k,2))
Info%q(i,j,k,iBz) = half * (Info%aux(i,j,k,1) + Info%aux(i,j,k+1,3))

Cooling

Two things are required to turn on cooling: the lCooling flag to indicate that cooling is active in this simulation, and iCooling to specify the type of cooling to use. These values are usually included in the problem.data file. The user must also create a cooling object in ProblemModuleInit() to manage the cooling settings. An example of cooling object creation can be seen below:

    IF(iCooling>0) THEN
       IF (.NOT. lRestart) THEN
           ! see sources/cooling.f90::CreateCoolingObject for
           ! default values of a cooling source term
           CALL CreateCoolingObject(coolingobj)
       ELSE
           coolingobj => firstcoolingobj
       END IF
    END IF

    coolingobj%iCooling=iCooling
    SELECT CASE(iCooling) ! cases defined in sources/cooling.f90
    CASE(NoCool)
    CASE(AnalyticCool)
       coolingobj%alpha=alpha
       coolingobj%beta=beta
    CASE(DMCool)
    CASE(IICool)
    CASE DEFAULT
    END SELECT

    coolingobj%floortemp=1d0
    coolingobj%mintemp=0.001

The .NOT. lRestart conditional prevents AstroBEAR from creating a new cooling object on restarts; this is because the cooling objects will be read in from the restart files.

Self-Gravity

AstroBEAR uses the ​hypre library to solve the self-gravity equations. To use self-gravity:

1. Look for the HYPREFLAG variable in Makefile.inc and make sure that it is set to 1.
2. Set the lSelfGravity flag in your physics.data file and set it to T.

Hypre will automatically initialize the potential field using the density. The only caveat is that the initial density cannot be uniform. When the density is uniform, hypre produces a ​singular matrix that it can't solve. Fortunately, a small density perturbation takes care of this problem without substantially affecting the dynamics of the domain. AstroBEAR comes with a Perturbation object type that can be used for this.

Sink Particles

The ability to form sink particles in AstroBEAR is tied to self-gravity. To simply enable sink particles:

1. Look for the HYPREFLAG variable in Makefile.inc and make sure that it is set to 1.
2. Set the lSelfGravity flag in your physics.data file and set it to T.

If you just want your simulation to have the option of forming sink particles, no further action is required. If you want to start your simulation off with sink particles, then you will have to create one in problem.f90::ProblemModuleInit():

      NAMELIST /ProblemData/ nParticles
      NAMELIST /ParticleData/ mass,xloc,vel
      OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data', STATUS="OLD")
      READ(PROBLEM_DATA_HANDLE,NML=ProblemData)

      IF (.NOT. lRestart) THEN

         DO i=1,nParticles

            READ(PROBLEM_DATA_HANDLE,NML=ParticleData)
            NULLIFY(Particle)
            CALL CreateParticle(Particle)
            Particle%mass=mass
            Particle%xloc=xloc
            Particle%vel=vel         
            CALL AddSinkParticle(Particle)
 
         END DO

         CLOSE(PROBLEM_DATA_HANDLE)      
         OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='restart.data', STATUS="UNKNOWN")
         WRITE(PROBLEM_DATA_HANDLE,NML=RestartData)
         CLOSE(PROBLEM_DATA_HANDLE)

      END IF

Depending on the features of your simulation, more objects might have to be declared in conjunction with the sink particle. The .NOT. lRestart conditional is important, as it prevents AstroBEAR from adding the same particle again on a restart.

Domain Objects

While they do not implement actual physics, AstroBEAR domain objects are useful enough to mention here. AstroBEAR 2.0 comes with a Shape object that can be used to define regions of space within a domain. These shape objects can be used by other objects and processes to impose special conditions within regions of the domain, even if they span multiple grids or processors. For more information, visit the AstroBEAR objects tutorial.