wiki:ModulesOnAstroBear

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.
  1. 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.
  1. 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 six 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."

The six subroutines can be split into two groups

Those that modify an individual AMR patch

Note these may get called multiple times on a given processor or not at all depending on how patches are distributed among the processors. The only assignment statements should be for variables local to each subroutine or the fluid data belonging to the info patch. Do not modify module variables or global variables within these subroutines.

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

And those that get called on by every processor at certain times.

Note here you can adjust module variables or modify the properties of Objects. You can also add global communication here if you desire though this is usually not necessary.

  • ProblemModuleInit(): Module variables are initialized here. This routine is where problem.data are read in and module-level namelists are populated and is called once per processor at the start (or restart) of a simulation. 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.
  • ProblemBeforeGlobalStep(level): This is called by every processor one per step per AMR level. If you need to modify variables other than the grid data - this is the correct place to do so. If you are using threading however - you should only modify variables between root steps (when level == 0).

Here is the bare minimum problem module

MODULE Problem
  IMPLICIT NONE
  SAVE
  PUBLIC ProblemModuleInit, ProblemGridInit, ProblemBeforeStep, &
         ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep
  PRIVATE

CONTAINS

  SUBROUTINE ProblemModuleInit()
  END SUBROUTINE ProblemModuleInit

  SUBROUTINE ProblemGridInit(Info)
    TYPE(InfoDef) :: Info
  END SUBROUTINE ProblemGridInit

  SUBROUTINE ProblemBeforeStep(Info)
    TYPE(InfoDef) :: Info
  END SUBROUTINE ProblemBeforeStep

  SUBROUTINE ProblemAfterStep(Info)
    TYPE(InfoDef) :: Info
  END SUBROUTINE ProblemAfterStep

  SUBROUTINE ProblemSetErrFlag(Info)
    TYPE(InfoDef) :: Info
  END SUBROUTINE ProblemSetErrFlag

  SUBROUTINE ProblemBeforeGlobalStep(n)
     INTEGER :: n
  END SUBROUTINE ProblemBeforeGlobalStep

END MODULE Problem



AstroBEAR Module Basics


Reading in Problem parameters and converting to computational units

Initializing a grid involves setting the fluid values at each cell with spatially averaged values of your initial coniditions. 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 uniform background (density = rhoAmb, Temperature = Tamb) and an overdense spherical clump centered at the origin with radius R, density rhoClump in pressure equilibrium with the ambient. We then want to add a constant wind of density rhoWind, Temperature TWind, and velocity Vwind coming from the left boundary.

First, we need to declare our variables and add them to the ProblemData Namelist. Then in ProblemModuleInit, we read them in from the problem.data file and convert all of the variables to computational units. Here I will assume that input variables are in cgs units (with Temperature in Kelvin). We are also going to declare additional derived variables to store the calculated wind pressure, clump pressure, and ambient pressure.

MODULE Problem
  USE GlobalDeclarations
  USE PhysicsDeclarations
  USE DataDeclarations
  IMPLICIT NONE
  SAVE
  PUBLIC ProblemModuleInit, ProblemGridInit, ProblemBeforeStep, &
         ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep
  PRIVATE 
  REAL(KIND=qPREC) :: rhoAmb, Tamb, rhoClump, RClump, rhoWind, TWind, vWind
  REAL(KIND=qPREC) :: pWind, pClump, pAmb

CONTAINS

  SUBROUTINE ProblemModuleInit()
    NAMELIST/ProblemData/ rhoAmb, Tamb, rhoClump, RClump, rhoWind, TWind, vWind
    OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data', STATUS="OLD")
    READ(PROBLEM_DATA_HANDLE,NML=ProblemData)
    CLOSE(PROBLEM_DATA_HANDLE)
    rhoAmb = rhoAmb / rScale
    rhoClump = rhoClump / rScale
    rhoWind = rhoWind / rScale
    TAmb = TAmb / TempScale
    TWind = TWind / TempScale
    RClump = RClump / lScale
    vWind = vWind / velScale

    pAmb = rhoAmb * TAmb
    pWind = rhoWind * TWind
    pClump = pAmb

  END SUBROUTINE

Using Template Functions

The preferred way to apply our initial/boundary functions to grid patches is to use template subroutines. This involves writing IC and BC subroutines that adhere to a certain form and then passing them to the ApplyTemplate subroutine. The ApplyTemplate loops over the grid cells and converts the q-array to primitive form (density, velocity, pressure) before calling the template function (in this case InitialCondition or BoundaryCondition.

 SUBROUTINE ProblemGridInit(Info)
    TYPE(InfoDef) :: Info
    CALL ApplyFunction(InitialCondition, Info)
  END SUBROUTINE ProblemGridInit

  SUBROUTINE ProblemBeforeStep(Info)
    TYPE(InfoDef) :: Info
    CALL ApplyFunction(BoundaryCondition, Info)
  END SUBROUTINE ProblemBeforeStep


  SUBROUTINE InitialCondition(pos, t, dt, q, Ghost)
    REAL(KIND=qPREC), DIMENSION(:) :: pos
    REAL(KIND=qPREC), DIMENSION(:) :: q
    REAL(KIND=qPREC) :: t, dt
    LOGICAL :: Ghost
    IF (sqrt(sum(pos**2)) < RClump) THEN
       q(1) = rhoClump
       q(imom(1:nDim))=0
       q(iE) = pClump
    ELSE
       q(1) = rhoAmb
       q(imom(1:nDim)) = 0
       q(iE) = pAmb
    END IF
  END SUBROUTINE InitialCondition


  SUBROUTINE BoundaryCondition(pos, t, dt, q, Ghost)
    REAL(KIND=qPREC), DIMENSION(:) :: pos
    REAL(KIND=qPREC), DIMENSION(:) :: q
    REAL(KIND=qPREC) :: t, dt
    LOGICAL :: Ghost
    IF (pos(1) < GxBounds(1,1)) THEN
       q(1) = rhoWind
       q(imom(1:nDim))=0
       q(ivx) = vWind
       q(iE) = pWind
    END IF
  END SUBROUTINE BoundaryCondition

 

Manually Looping Over Zones

  SUBROUTINE ProblemGridInit(Info)
    TYPE(InfoDef) :: Info
    INTEGER :: i,j,k
    REAL(KIND=xPREC) :: pos(3)

    ! Loop over cells in patch
    DO i=1, Info%mX(1)
      DO j=1, Info%mX(2)
        DO k=1, Info%mX(3)
          pos=CellPos(Info, i, j, k)
          IF (sqrt(sum(pos**2)) < RClump) THEN
            Info%q(i,j,k,irho) = rhoClump
            Info%q(i,j,k,imom(1:nDim)) = 0
            Info%q(i,j,k,iE) = gamma7*pClump
          ELSE
            Info%q(i,j,k,irho) = rhoAmb
            Info%q(i,j,k,imom(1:nDim)) = 0
            Info%q(i,j,k,iE) = gamma7*pAmb
          END IF                  
        END DO
      END DO
    END DO
  END SUBROUTINE

  SUBROUTINE ProblemBeforeStep(Info)
    TYPE(InfoDef) :: Info
    INTEGER :: i,j,k, mbc(3)
    REAL(KIND=xPREC) :: pos(3)
    
    !determine number of ghost zones for each dimension
    mbc=levels(Info%level)%gmbc(levels(Info%level)%step)*merge((/1,1,1/),(/0,0,0/),nDim>=(/1,2,3/))

    ! Initialize wind in leftmost boundary
    DO i=1-mbc(1), Info%mX(1)+mbc(1)
      DO j=1-mbc(2), Info%mX(2)+mbc(2)
        DO k=1-mbc(3), Info%mX(3)+mbc(3)
          pos=CellPos(Info, i, j, k)
          IF (pos(1) < GxBounds(1,1)) THEN
            Info%q(i,j,k,irho) = rhoWind
            Info%q(i,j,k,ivx)=vWind*rhoWind
            IF (iE /= 0) Info%q(i,j,k, iE)=gamma7*pWind+.5d0*rhoWind*vWind**2
            IF (ivy /= 0) Info%q(i,j,k, ivy)=0d0
            IF (ivz /= 0) Info%q(i,j,k, ivz)=0d0
          END IF                  
        END DO
      END DO
    END DO
  END SUBROUTINE

Using Prebuilt Objects

There are a lot of other ways we could modify this simple example to have the clump be located anywhere, to have a density profile that is smoothed at the edge, to be a different temperature, or move with a particular velocity, etc… Fortunately, clumps are a commonly used object (as are uniform backgrounds and winds) and there are modules designed to assist users in easily creating clumps, uniform backgrounds, and winds. For example the above module could be rewritten using an Ambient Object, a Clump Object, and a Wind Object.

MODULE Problem
  USE GlobalDeclarations
  USE DataDeclarations
  USE Clumps
  USE Ambients
  USE Winds
  IMPLICIT NONE
  SAVE
  PUBLIC ProblemModuleInit, ProblemGridInit, ProblemBeforeStep, &
         ProblemAfterStep, ProblemSetErrFlag, ProblemBeforeGlobalStep
  PRIVATE 
  REAL(KIND=qPREC) :: rhoAmb, Tamb, rhoClump, RClump, rhoWind, TWind, vWind
  REAL(KIND=qPREC) :: pWind, pClump, pAmb

CONTAINS

  SUBROUTINE ProblemModuleInit()
    TYPE(AmbientDef), POINTER :: Ambient
    TYPE(ClumpDef), POINTER :: Clump
    TYPE(WindDef), POINTER :: Wind
    NAMELIST/ProblemData/ rho, radius
    OPEN(UNIT=PROBLEM_DATA_HANDLE, FILE='problem.data', STATUS="OLD")
    READ(PROBLEM_DATA_HANDLE,NML=ProblemData)
    CLOSE(PROBLEM_DATA_HANDLE)

    rhoAmb = rhoAmb / rScale
    rhoClump = rhoClump / rScale
    rhoWind = rhoWind / rScale
    TAmb = TAmb / TempScale
    TWind = TWind / TempScale
    RClump = RClump / lScale
    vWind = vWind / velScale

    pAmb = rhoAmb * TAmb
    pWind = rhoWind * TWind
    pClump = pAmb

    CALL CreateAmbient(Ambient)
    Ambient%density = rhoAmb
    Ambient%pressure = pAmb
    CALL UpdateAmbient(Ambient)

    CALL CreateClump(Clump)
    Clump%density = rhoClump
    Clump%radius = RClump
    Clump%pressure = pClump
    CALL UpdateClump(Clump)

    CALL CreateWind(Wind)
    Wind%density = rhoWind
    Wind%velocity = vWind
    Wind%temperature = TWind
    Wind%direction = 1
    Wind%edge = 1
    CALL UpdateWind(Wind)

  END SUBROUTINE
  
  SUBROUTINE ProblemGridInit(Info)
    TYPE(InfoDef) :: Info
  END SUBROUTINE
 
  SUBROUTINE ProblemBeforeStep(Info)
    TYPE(InfoDef) :: Info
  END SUBROUTINE ProblemBeforeStep

  SUBROUTINE ProblemAfterStep(Info)
    TYPE(InfoDef) :: Info
  END SUBROUTINE ProblemAfterStep

  SUBROUTINE ProblemSetErrFlag(Info)
    TYPE(InfoDef) :: Info
  END SUBROUTINE ProblemSetErrFlag

  SUBROUTINE ProblemBeforeGlobalStep(n)
     INTEGER :: n
  END SUBROUTINE ProblemBeforeGlobalStep

END MODULE

So what is going on here?

And we're done. We don't have to worry about dx or mx or cell positions, or the number of ghost zones etc… All of that detailed work is done for us. It is important to note that the order that objects are created is the same order they are placed on the grid. So had we created the clump object first - the clump data would have been overwritten by the ambient module and there would be no clump.

If we want to get more complicated - we can modify other clump/wind/ambient attributes. All of the available (and default) options for the various objects should be documented on the AstroBearObjects page.

To try the problem module we just built, we can copy the data files from $(ASTROBEAR)/modules/Template directory to the current problem directory $(ASTROBEAR)/modules/Problem and add these three lines to problem.data

&ProblemData
rhoAmb = 1.0
TAmb = 1.0
rhoClump = 10.0
rClump = 1.0
rhoWind = 1.0
TWind = 1.0
vWind = 10.0
/

If we follow the procedure Setting up and Compiling A Problem to try to compile and run our new problem module. If you set the dimensions and final time in global.data as

nDim     = 2                ! number of dimensions for this problem (1-3)
GmX      = 30,30,1          ! Base grid resolution [x,y,z]
MaxLevel = 0                        ! Maximum level for this simulation (0 is fixed grid)

and

final_time        = 15d-1               ! The final time in computational units.
final_frame       = 20              ! The final frame [10]

and follow First Run to run it. Following Chapter 2 to use VisIT to analyze the results, we can get a movie like this ClumpMovie


Flagging Cells for Refinement

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

  Info%ErrFlag(1:mx,1:my,1:mz)

To clear the cell at (i,j,k), simply set Info%ErrFlag(i,j,k) to 0. An error flag of 0 means that the cell does not need to be refined. This is in general, not a good idea since a previous routine might have already flagged that cell for refinement with good reason. To mark a cell for refinement, set Info%ErrFlag(i,j,k) to 1. The place to do this is in the ProblemSetErrFlag() routine; most conventional physical criteria for refinement are already handled by AstroBEAR itself. For more information see ControllingRefinement for more information as well as ProblemSetErrFlag for an example of how to use this subroutine.


A few notes on magnetic aux fields

Aux fields are particularly difficult to work with - and if initialized improperly (either form a non-divergenceless physical model or from 2nd order errors due to estimating face averages with face-centered values) will produce a probably small but not insignificant divergence that will stick around for the course of the simulation. The easiest way to avoid divergence in your B-fields is to first calculate the vector potential and then to take the curl discretely. In 2D, this means calculating the value of the vector potential at cell corners (which only has a z-component)- and then differencing them in y to get Bx and differencing them in x to get -By. In 3D, the vector potential should be averaged along each edge of the component parallel to that edge. For example Ax should be averaged along each edge that is parallel to the x-axis. The 2nd order errors due to estimating the average value along an edge by the midpoint will not produce divergence in the resulting B-field.

Let's suppose in we want to initialize the grid with a B-field By=sin(x). Well the vector potential would just be Az=cos(x) and our ProblemGridInit routine would look like:

  SUBROUTINE ProblemGridInit(Info)
    TYPE(InfoDef) :: Info
    INTEGER :: i,j,k
    REAL(KIND=qPREC) :: pos(3), dx
    dx = levels(Info%level)%dx
    IF (MaintainAuxArrays) THEN
      Info%aux(1:Info%mX(1)+1,1:Info%mX(2), 1:Info%mX(3), 1)=0d0
      IF (nDim == 3) Info%aux(1:Info%mX(1),1:Info%mX(2), 1:Info%mX(3)+1, 3)=0d0
      DO i=1, Info%mX(1)
        DO j=1, Info%mX(2)+1
          DO k=1, Info%mX(3)
            pos=CellPos(Info, i, j, k)          
            Info%aux(i,j,k,2)=(cos(pos(1)+half*dx)-cos(pos(1)-half*dx))/dx
          END DO
        END DO
      END DO
      DO i=1, Info%mX(1)
        DO j=1, Info%mX(2)
          DO k=1, Info%mX(3)
            Info%q(i,j,k,iBx)=0d0
            Info%q(i,j,k,iBy)=half*(Info%aux(i,j,k,2)+Info%aux(i,j+1,k,2))
            Info%q(i,j,k,iBz)=0d0
            Info%q(i,j,k,iE)=gamma7+half*Info%q(i,j,k,iBy)**2
          END DO
        END DO
      END DO
    ELSE
      DO i=1, Info%mX(1)
        DO j=1, Info%mX(2)
          DO k=1, Info%mX(3)
            pos=CellPos(Info, i, j, k)
            Info%q(i,j,k,iBx)=0d0
            Info%q(i,j,k,iBy)=sin((pos(1))
            Info%q(i,j,k,iBz)=0d0
          END DO
        END DO
      END DO
    END IF
  END SUBROUTINE

We can also use the ApplyEMF or ApplyBField routines just as we would use ApplyTemplate but the subroutine should return either the Vector Potential (EMF) or the B field respectively. You should only use ApplyBField if your B-Field is constant - otherwise use the vector potential with ApplyEMF.

  SUBROUTINE ProblemGridInit(Info)
    CALL ApplyEMF(myEMF, Info)
  END SUBROUTINE

  SUBROUTINE myEMF(pos, t, dt, q)
    REAL(KIND=qPREC), DIMENSION(:) :: pos
    REAL(KIND=qPREC), DIMENSION(:) :: q
    REAL(KIND=qPREC) :: t, dt
    q=(/0d0,0d0,cos(pos(1))/)
  END SUBROUTINE myEMF

Summary

Writing a problem module can be extremely complicated, but for first time users writing relatively simple modules, here are a few helpful tips:


Last modified 8 years ago Last modified on 04/24/17 16:24:10

Attachments (3)

Note: See TracWiki for help on using the wiki.