3.5 Creating a Simulation

The routine Simulation_initBlock is called by the Grid unit to apply initial conditions to the physical domain. If the AMR grid PARAMESH is being used, the formation of the physical domain starts at the lowest level of refinement. Initial conditions are applied to each block at this level by calling Simulation_initBlock. The Grid unit then checks the refinement criteria in the blocks it has created and refines the blocks if the criteria are met. It then calls Simulation_initBlock to initialize the newly created blocks. This process repeats until the grid reaches the required refinement level in the areas marked for refinement. The Uniform Grid has only one level, with same resolution everywhere. Therefore, only one block per processor is created and Simulation_initBlock is called to initialize this single block. It is important to note that a problem's Simulation_initBlock routine is the same regardless of whether PARAMESH or Uniform Grid is being used. The Grid unit handles these differences, not the Simulation unit.

A use statement for the Simulation_data
One of more use statement to access other unit interfaces being used, for example use Grid_interface, ONLY: Grid_putPointData
Variable typing implicit none statement
Necessary #include header files
Declaration of arguments and local variables.
Generation of initial conditions either from a file, or directly calculated in the routine
Calls to the various Grid_putData routines to store the values of solution variables.

We continue to look at the Sod setup and describe its Simulation_initBlock in detail. The first part of the routine contains all the declarations as shown below. The first statement in routine is the use statement, which provides access to the runtime parameters and other unit scope variables initialized in the Simulation_init routine. The include files bring in the needed constants, and then the arguments are defined. The declaration of the local variables is next, with allocatable arrays for each block.



subroutine Simulation_initBlock(blockID)
  ! get the needed unit scope data
  use Simulation_data, ONLY: sim_posn, sim_xCos, sim_yCos, sim_zCos,&
                             sim_rhoLeft, sim_pLeft, sim_uLeft,     &
                             sim_rhoRight, sim_pRight, sim_uRight,  &
                             sim_smallX, sim_gamma, sim_smallP
  use Grid_interfaces, ONLY : Grid_getBlkIndexLimits, Grid_getCellCoords, &
                              Grid_putPointData
  implicit none


! get all the constants
#include "constants.h"
#include "Flash.h"


! define arguments and indicate whether they are input or output
  integer, intent(in) :: blockID


! declare all local variables.
  integer :: i, j, k, n
  integer :: iMax, jMax, kMax
  real :: xx, yy,  zz, xxL, xxR
  real :: lPosn0, lPosn


! arrays to hold coordinate information for the block
  real, allocatable, dimension(:) :: xCenter, xLeft, xRight, yCoord, zCoord


! array to get integer indices defining the beginning and the end
  ! of a block.
  integer, dimension(2,MDIM) :: blkLimits, blkLimitsGC


! the number of grid points along each dimension
  integer :: sizeX, sizeY, sizeZ


integer, dimension(MDIM) :: axis
  integer :: dataSize
  logical :: gcell = .true.


! these variables store the calculated initial values of physical
  ! variables a grid point at a time.
  real :: rhoZone, velxZone, velyZone, velzZone, presZone, &
       enerZone, ekinZone

Note that FLASH promotes all floating point variables to double precision at compile time for maximum portability. We therefore declare all floating point variables with real in the source code. In the next part of the code we allocate the arrays that will hold the coordinates.

FLASH3 Transition: FLASH4-alpha supports blocks that are not sized at compile time to generalize the Uniform Grid, and to be able to support different AMR packages in future. For this reason, the arrays are not sized with the static NXB etc. as was the case in FLASH2. Instead they are allocated on a block by block basis in Simulation_initBlock. Performance is compromised by the use of allocatable arrays, however, since this part of the code is executed only at the beginning of the simulation, it has negligible impact on the overall execution time in production runs.

! get the integer endpoints of the block in all dimensions
  ! the array blkLimits returns the interior end points
  ! whereas array blkLimitsGC returns endpoints including guardcells
  call Grid_getBlkIndexLimits(blockId,blkLimits,blkLimitsGC)


! get the size along each dimension for allocation and then allocate
  sizeX = blkLimitsGC(HIGH,IAXIS)
  sizeY = blkLimitsGC(HIGH,JAXIS)
  sizeZ = blkLimitsGC(HIGH,KAXIS)
  allocate(xLeft(sizeX))
  allocate(xRight(sizeX))
  allocate(xCenter(sizeX))
  allocate(yCoord(sizeY))
  allocate(zCoord(sizeZ))

The next part of the routine involves setting up the initial conditions. This section could be code for interpolating a given set of initial conditions, constructing some analytic model, or reading in a table of initial values. In the present example, we begin by getting the coordinates for the cells in the current block. This is done by a set of calls to Grid_getCellCoords . Next we create loops that compute appropriate values for each grid point, since we are constructing initial conditions from a model. Note that we use the blkLimits array from Grid_getBlkIndexLimits in looping over the spatial indices to initialize only the interior cells in the block. To initialize the entire block, including the guardcells, the blkLimitsGC array should be used.



xCoord(:) = 0.0
  yCoord(:) = 0.0
  zCoord(:) = 0.0
  call Grid_getCellCoords(IAXIS, blockID, LEFT_EDGE,  gcell, xLeft,   sizeX)
  call Grid_getCellCoords(IAXIS, blockID, CENTER,     gcell, xCenter, sizeX)
  call Grid_getCellCoords(IAXIS, blockID, RIGHT_EDGE, gcell, xRight,  sizeX)
  call Grid_getCellCoords(JAXIS, blockID, CENTER,     gcell, yCoord,  sizeY)
  call Grid_getCellCoords(KAXIS, blockID, CENTER,     gcell, zCoord,  sizeZ)


!---------------------------------------
! loop over all of the zones in the current block and set the variables.
!---------------------------------------
  do k = blkLimits(LOW,KAXIS),blkLimits(HIGH,KAXIS)
     zz = zCoord(k) ! coordinates of the cell center in the z-direction
     lPosn0 = sim_posn - zz*sim_zCos/sim_xCos ! Where along the x-axis
                                              ! the shock intersects
                                              ! the xz-plane at the current z.


do j = blkLimits(LOW,JAXIS),blkLimits(HIGH,JAXIS)
        yy = yCoord(j) ! center coordinates in the y-direction
        lPosn = lPosn0 - yy*sim_yCos/sim_xCos ! The position of the
                                              ! shock in the current yz-row.
        dataSize = 1  ! for Grid put data function, we are
                      ! initializing a single cell at a time and
                      ! sending it to Grid


do i = blkLimits(LOW,IAXIS),blkLimits(HIGH,IAXIS)


xx  = xCenter(i) ! center coordinate along x
           xxL = xLeft(i)   ! left coordinate along y
           xxR = xRight(i)  ! right coordinate along z

For the present problem, we create a discontinuity along the shock plane. We do this by initializing the grid points to the left of the shock plane with one value, and the grid points to the right of the shock plane with another value. Recall that the runtime parameters which provide these values are available to us through the Simulation_data module. At this point we can initialize all independent physical variables at each grid point. The following code shows the contents of the loops. Don't forget to store the calculated values in the Grid data structure!

if (xxR <= lPosn) then
     rhoZone = sim_rhoLeft
     presZone = sim_pLeft


velxZone = sim_uLeft * sim_xCos
     velyZone = sim_uLeft * sim_yCos
     velzZone = sim_uLeft * sim_zCos


! initialize cells which straddle the shock.  Treat them as though 1/2 of
     ! the cell lay to the left and 1/2 lay to the right.
  elseif ((xxL < lPosn) .and. (xxR > lPosn)) then


rhoZone = 0.5 * (sim_rhoLeft+sim_rhoRight)
     presZone = 0.5 * (sim_pLeft+sim_pRight)


velxZone = 0.5 *(sim_uLeft+sim_uRight) * sim_xCos
     velyZone = 0.5 *(sim_uLeft+sim_uRight) * sim_yCos
     velzZone = 0.5 *(sim_uLeft+sim_uRight) * sim_zCos


! initialize cells to the right of the initial shock.
  else


rhoZone = sim_rhoRight
     presZone = sim_pRight


velxZone = sim_uRight * sim_xCos
     velyZone = sim_uRight * sim_yCos
     velzZone = sim_uRight * sim_zCos


endif
   axis(IAXIS) = i   ! Get the position of the cell in the block
   axis(JAXIS) = j
   axis(KAXIS) = k


! Compute the gas energy and set the gamma-values
   ! needed for the equation of  state.


ekinZone = 0.5 * (velxZone**2 + velyZone**2 + velzZone**2)


enerZone = presZone / (sim_gamma-1.)
   enerZone = enerZone / rhoZone
   enerZone = enerZone + ekinZone
   enerZone = max(enerZone, sim_smallP)


! store the variables in the current zone via the Grid_putPointData method


call Grid_putPointData(blockId, CENTER, DENS_VAR, EXTERIOR, axis, rhoZone)
   call Grid_putPointData(blockId, CENTER, PRES_VAR, EXTERIOR, axis, presZone)
   call Grid_putPointData(blockId, CENTER, VELX_VAR, EXTERIOR, axis, velxZone)
   call Grid_putPointData(blockId, CENTER, VELY_VAR, EXTERIOR, axis, velyZone)
   call Grid_putPointData(blockId, CENTER, VELZ_VAR, EXTERIOR, axis, velzZone)
   call Grid_putPointData(blockId, CENTER, ENER_VAR, EXTERIOR, axis, enerZone)
   call Grid_putPointData(blockId, CENTER, GAME_VAR, EXTERIOR, axis, sim_gamma)
   call Grid_putPointData(blockId, CENTER, GAMC_VAR, EXTERIOR, axis, sim_gamma)

When Simulation_initBlock returns, the Grid data structures for physical variables have the values of the initial model for the current block. As mentioned before, Simulation_initBlock is called for every block that is created as the code refines the initial model.

3.5 Creating a Simulation_initBlock.F90

3.5 Creating a `Simulation_initBlock.F90`