Setting Up a Simulation

pySimulation

A simulation in pyFlash is specified by providing a Simulation class definition in pySimulation.py` that inherits from the pySimulation and overriding the stub implementations of the methods you wish to use. Your simulation files should be put in the Simulation/SimulationMain/python directory to get the correct bindings to the Simulation_*.F90 FLASH routines.

Simulation/SimulationMain/python/Sod/pySimulation.py

import flash.pyFlash4.Grid as gr
import flash.parm as parm
from   flash.FlashH import *
from   flash.Simulation_base import pySimulation

class Simulation(pySimulation):
    posn = 0.5

    def init(self):
        xmin = parm.xmin.getVal()
        xmax = parm.xmax.getVal()
        self.posn = 0.5*(xmax-xmin)

    def initBlock(self, blockID):
        myBlock = gr.Block(blockID)
        for k in range(myBlock.hiGC[2]):
            for j in range(myBlock.hiGC[1]):
                for i in range(myBlock.hiGC[0]):
                    x = myBlock.xCoord[i]
                    if x < self.posn:
                        myBlock[DENS_VAR,i,j,k] = 1.
                        myBlock[PRES_VAR,i,j,k] = 1.
                    else:
                        myBlock[DENS_VAR,i,j,k] = 0.1
                        myBlock[PRES_VAR,i,j,k] = 0.125

Simulation.init()

This method drops in for the Simulation_init.F90 subroutine and is separate from a python class constructor __init__, such that it won’t be called to initialize the class member attributes until after the rest of FLASH may be initialized.

def init(self):
    xmin = parm.xmin.getVal()
    xmax = parm.xmax.getVal()
    self.posn = 0.5*(xmax-xmin)

This implementation makes use of the flash.pyFlash4.RP.rpReal() API to get the values of runtimes that can be used in the rest of the Simulation.

Simulation.initBlock(int: blockID)

The initBlock() method initializes the variables on a single grid block.

def initBlock(self, blockID):
    myBlock = gr.Block(blockID)
    for k in range(myBlock.hiGC[2]):
        for j in range(myBlock.hiGC[1]):
            for i in range(myBlock.hiGC[0]):
                x = myBlock.xCoord[i]
                if x < self.posn:
                    myBlock[DENS_VAR,i,j,k] = 1.
                    myBlock[PRES_VAR,i,j,k] = 1.
                else:
                    myBlock[DENS_VAR,i,j,k] = 0.1
                    myBlock[PRES_VAR,i,j,k] = 0.125

Here we make use of the hiGC attribute to get the index ranges of the block dimensions and the __setitem__() to directly set the values in the block pointer. Coordinate information is directly taken from the xCoord array.

As an alternative Block objects can be cast to numpy arrays and the grid block pointer set through the array by passing the copy=False argument.

Simulation/SimulationMain/python/magnetoHD/rotor/pySimulation.py

 # use np.ndindex to get the iterator
 faceX = gr.Block(blockID, FACEX)
 solnX = np.array(faceX, copy=False)
 for (k,j,i) in np.ndindex(solnX.shape[:3]):
     solnX[MAG_FACE_VAR,i,j,k] = 5./np.sqrt(4.*np.pi)

Note

The *_VAR definitions in flash.FlashH.py are 0-indexed to match the indexing in python as opposed to being 1-indexed like in the setup generated Flash.h header for the fortran source code.

Simulation.adjustEvolution()

Users can interact with the grid data during each time step similarly to the Simulation_adjustEvolution.F90 method.

Simulation/SimulationMain/python/magnetoHD/rotor/pySimulation.py

 def adjustEvolution(self):
     blockList = gr.BlockList(LEAF)
     for lb in range(blockList.nBlocks):
         myBlock = blockList(lb, CENTER)

         x = np.array(myBlock.xCoord)
         y = np.array(myBlock.yCoord)
         z = np.array(myBlock.zCoord)
         xx, yy, zz = np.meshgrid(x,y,z,indexing='ij')

         solnVec = np.array(myBlock, copy=False)
         testVar = solnVec[TEST_VAR,:,:,:]
         magpVar = solnVec[MAGP_VAR,:,:,:]

         testVar[xx<0.5] = magpVar[xx<0.5]

Boundary Conditions

User defined hydrodynamic boundary conditions can be implemented through the bc[IJK][LO|HI] methods.

Simulation/SimulationMain/python/DoubleMachReflection/pySimulation.py

import flash.pyFlash4.Grid as grid
import flash.parm as p
from   flash.FlashH import *

 def init(self):
     gamma = p.gamma.getVal()
     UL = np.zeros(NUNK_VARS)
     UL[DENS_VAR] = self.dL
     UL[PRES_VAR] = self.pL
     UL[VELX_VAR] = self.uL
     UL[VELY_VAR] = self.vL
     UL[GAMC_VAR] = gamma
     UL[GAME_VAR] = gamma
     UL[EINT_VAR] = self.pL/(self.dL*(gamma-1.))
     UL[ENER_VAR] = UL[EINT_VAR] + 0.5*(self.uL**2 + self.vL**2)
     self.UL = UL

     UR = np.zeros(NUNK_VARS)
     UR[DENS_VAR] = self.dR
     UR[PRES_VAR] = self.pR
     UR[VELX_VAR] = self.uR
     UR[VELY_VAR] = self.vR
     UR[GAMC_VAR] = gamma
     UR[GAME_VAR] = gamma
     UR[EINT_VAR] = self.pR/(self.dR*(gamma-1.))
     UR[ENER_VAR] = UR[EINT_VAR] + 0.5*(self.uR**2 + self.vR**2)
     self.UR = UR

     self.ymaxRot = p.ymax.getVal()/np.tan(self.xAngle)

 def bcJLO(self,region):
     # left of the initial shock position we maintain the inflow
     # reflecting everywhere else
     solnVec = np.array(region,copy=False)
     xx = np.array(region.coord2)
     for var in range(solnVec.shape[3]):
         ref = 1.
         if (var == VELY_VAR):
             ref = -1.
         for (i,j,k) in np.ndindex(solnVec.shape[:3]):
             if xx[k] < self.posn:
                 solnVec[i,j,k,var] = self.UL[var]
             else:
                 solnVec[i,j,k,var] = ref*solnVec[2*region.nGuard-1-i,j,k,var]

These routines are passed an object of type flash.pyFlash4.Grid.bcRegion. The underlying region pointers are always layed out in memory so that the outer most index corresponds to the variable and the innermoste corresponds to the direction normal the the boundary being filled, (i1,i2,i3,var).

Warning

These routines will be passed both CENTER regions as well as FACE[XYZ] regions when NFACE_VARS>0 in the simulation. Users should check the value of flash.pyFlash4.Grid.bcRegion.dataStruct to know wheter data is on faces or not, and which one. This is particularly important in MHD simulations.