####################### Setting Up a Simulation ####################### pySimulation ============ .. toctree:: pySimulation A simulation in **pyFlash** is specified by providing a ``Simulation`` class definition in `pySimulation.py`` that inherits from the :py:class:`~flash.Simulation_base.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. .. code-block:: python :caption: 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. .. code-block:: python def init(self): xmin = parm.xmin.getVal() xmax = parm.xmax.getVal() self.posn = 0.5*(xmax-xmin) This implementation makes use of the :py:meth:`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 :py:meth:`~flash.Simulation_base.pySimulation.initBlock` method initializes the variables on a single grid block. .. code-block:: python 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 :py:data:`~flash.pyFlash4.Grid.Block.hiGC` attribute to get the index ranges of the block dimensions and the :py:meth:`~flash.pyFlash4.Grid.Block.__setitem__` to directly set the values in the block pointer. Coordinate information is directly taken from the :py:data:`~flash.pyFlash4.Grid.Block.xCoord` array. As an alternative :py:class:`~flash.pyFlash4.Grid.Block` objects can be cast to `numpy `_ arrays and the grid block pointer set through the array by passing the ``copy=False`` argument. .. code-block:: python :caption: 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. .. code-block:: python :caption: 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. .. code-block:: python :caption: 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 :py:class:`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.