!!****if* source/Simulation/SimulationMain/LaserSlab/Simulation_initBlock !! !! NAME !! !! Simulation_initBlock !! !! !! SYNOPSIS !! !! call Simulation_initBlock(integer(IN) :: blockID) !! !! !! !! !! DESCRIPTION !! !! Initializes fluid data (density, pressure, velocity, etc.) for !! a specified block. !! !! ARGUMENTS !! !! blockID - the number of the block to initialize !! !! !! !!*** subroutine Simulation_initBlock(blockId) use Simulation_data use Grid_interface, ONLY : Grid_getBlkIndexLimits, & Grid_getCellCoords, Grid_putPointData,& Grid_getBlkPtr, & Grid_releaseBlkPtr use Driver_interface, ONLY: Driver_abortFlash use RadTrans_interface, ONLY: RadTrans_mgdEFromT implicit none #include "constants.h" #include "Flash.h" ! compute the maximum length of a vector in each coordinate direction ! (including guardcells) integer, intent(in) :: blockId integer :: i, j, k, n integer :: blkLimits(2, MDIM) integer :: blkLimitsGC(2, MDIM) integer :: axis(MDIM) real, allocatable :: xcent(:), ycent(:), zcent(:) real :: tradActual real :: rho, tele, trad, tion, zbar, abar integer :: species real, pointer, dimension(:,:,:,:) :: solnData, facexData,faceyData #if NDIM > 0 real, pointer, dimension(:,:,:,:) :: facezData #endif #ifndef CHAM_SPEC integer :: CHAM_SPEC = 1, TARG_SPEC = 2 #endif ! get the coordinate information for the current block from the database call Grid_getBlkIndexLimits(blockId,blkLimits,blkLimitsGC) ! get the coordinate information for the current block from the database call Grid_getBlkIndexLimits(blockId,blkLimits,blkLimitsGC) allocate(xcent(blkLimitsGC(HIGH, IAXIS))) call Grid_getCellCoords(IAXIS, blockId, CENTER, .true., & xcent, blkLimitsGC(HIGH, IAXIS)) allocate(ycent(blkLimitsGC(HIGH, JAXIS))) call Grid_getCellCoords(JAXIS, blockId, CENTER, .true., & ycent, blkLimitsGC(HIGH, JAXIS)) allocate(zcent(blkLimitsGC(HIGH, KAXIS))) call Grid_getCellCoords(KAXIS, blockId, CENTER, .true., & zcent, blkLimitsGC(HIGH, KAXIS)) call Grid_getBlkPtr(blockID,solnData,CENTER) #if NFACE_VARS > 0 if (sim_killdivb) then call Grid_getBlkPtr(blockID,facexData,FACEX) call Grid_getBlkPtr(blockID,faceyData,FACEY) if (NDIM>2) call Grid_getBlkPtr(blockID,facezData,FACEZ) endif #endif !------------------------------------------------------------------------------ ! Loop over cells and set the initial state do k = blkLimitsGC(LOW,KAXIS),blkLimitsGC(HIGH,KAXIS) do j = blkLimitsGC(LOW,JAXIS),blkLimitsGC(HIGH,JAXIS) do i = blkLimitsGC(LOW,IAXIS),blkLimitsGC(HIGH,IAXIS) axis(IAXIS) = i axis(JAXIS) = j axis(KAXIS) = k species = CHAM_SPEC if (sim_initGeom == "slab") then if(NDIM == 1) then if ( xcent(i) <= sim_targetHeight + sim_vacuumHeight .and. & xcent(i) >= sim_vacuumHeight ) then species = TARG_SPEC end if elseif(NDIM == 2 .or. NDIM == 3) then if (xcent(i) <= sim_targetRadius .and. & ycent(j) <= sim_targetHeight + sim_vacuumHeight .and. & ycent(j) >= sim_vacuumHeight ) then species = TARG_SPEC end if if ( xcent(i) <= 0.02 .and. & ycent(j) <= 0.024 + sim_vacuumHeight .and. & ycent(j) >= sim_vacuumHeight ) then species = CHAM_SPEC end if end if else if (sqrt(xcent(i)**2+ycent(j)**2+zcent(k)**2)<= sim_targetRadius) then species = TARG_SPEC end if end if if(species == TARG_SPEC) then rho = sim_rhoTarg tele = sim_teleTarg tion = sim_tionTarg trad = sim_tradTarg else rho = sim_rhoCham tele = sim_teleCham tion = sim_tionCham trad = sim_tradCham end if call Grid_putPointData(blockId, CENTER, DENS_VAR, EXTERIOR, axis, rho) call Grid_putPointData(blockId, CENTER, TEMP_VAR, EXTERIOR, axis, tele) #ifdef FLASH_3T call Grid_putPointData(blockId, CENTER, TION_VAR, EXTERIOR, axis, tion) call Grid_putPointData(blockId, CENTER, TELE_VAR, EXTERIOR, axis, tele) ! Set up radiation energy density: call RadTrans_mgdEFromT(blockId, axis, trad, tradActual) call Grid_putPointData(blockId, CENTER, TRAD_VAR, EXTERIOR, axis, tradActual) #endif if (NSPECIES > 0) then ! Fill mass fractions in solution array if we have any SPECIES defined. ! We put nearly all the mass into either the Xe material if XE_SPEC is defined, ! or else into the first species. do n = SPECIES_BEGIN,SPECIES_END if (n==species) then call Grid_putPointData(blockID, CENTER, n, EXTERIOR, axis, 1.0e0-(NSPECIES-1)*sim_smallX) else call Grid_putPointData(blockID, CENTER, n, EXTERIOR, axis, sim_smallX) end if enddo end if #ifdef BDRY_VAR call Grid_putPointData(blockId, CENTER, BDRY_VAR, EXTERIOR, axis, -1.0) if(NDIM == 2) then if(xcent(i) >= 0.015 .and. xcent(i) <= 0.045 .and. & ycent(j) >= 0.2 + sim_targetHeight + sim_vacuumHeight .and. & ycent(j) <= 0.2 + sim_targetHeight + sim_vacuumHeight + sim_grid) then call Grid_putPointData(blockId, CENTER, BDRY_VAR, EXTERIOR, axis, 1.0) else if(xcent(i) >= 0.075 .and. xcent(i) <= 0.105 .and. & ycent(j) >= 0.2 + sim_targetHeight + sim_vacuumHeight .and. & ycent(j) <= 0.2 + sim_targetHeight + sim_vacuumHeight + sim_grid) then call Grid_putPointData(blockId, CENTER, BDRY_VAR, EXTERIOR, axis, 1.0) else if(xcent(i) >= 0.135 .and. xcent(i) <= 0.165 .and. & ycent(j) >= 0.2 + sim_targetHeight + sim_vacuumHeight .and. & ycent(j) <= 0.2 + sim_targetHeight + sim_vacuumHeight + sim_grid) then call Grid_putPointData(blockId, CENTER, BDRY_VAR, EXTERIOR, axis, 1.0) else if(xcent(i) >= 0.195 .and. xcent(i) <= 0.225 .and. & ycent(j) >= 0.2 + sim_targetHeight + sim_vacuumHeight .and. & ycent(j) <= 0.2 + sim_targetHeight + sim_vacuumHeight + sim_grid) then call Grid_putPointData(blockId, CENTER, BDRY_VAR, EXTERIOR, axis, 1.0) else if(xcent(i) >= 0.255 .and. xcent(i) <= 0.285 .and. & ycent(j) >= 0.2 + sim_targetHeight + sim_vacuumHeight .and. & ycent(j) <= 0.2 + sim_targetHeight + sim_vacuumHeight + sim_grid) then call Grid_putPointData(blockId, CENTER, BDRY_VAR, EXTERIOR, axis, 1.0) end if end if #endif solnData(MAGX_VAR,i,j,k) = 0.0 solnData(MAGY_VAR,i,j,k) = sim_magy !!solnData(MAGZ_VAR,i,j,k) = 0.0 #if NFACE_VARS > 0 !! In this case we initialized Az using the cell-cornered coordinates. if (sim_killdivb) then if (NDIM >= 2) then facexData(MAG_FACE_VAR,i,j,k)= 0.0 faceyData(MAG_FACE_VAR,i,j,k)= sim_magy if (NDIM>2) facezData(MAG_FACE_VAR,i,j,k)= 0.0 endif endif #endif enddo enddo enddo #if NFACE_VARS > 0 if (sim_killdivb) then do k=blkLimitsGC(LOW,KAXIS),blkLimitsGC(HIGH,KAXIS) i=blkLimitsGC(HIGH,IAXIS)+1 do j=blkLimitsGC(LOW,JAXIS),blkLimitsGC(HIGH,JAXIS) facexData(MAG_FACE_VAR,i,j,k)= 0.0 enddo j=blkLimitsGC(HIGH,JAXIS)+1 do i=blkLimitsGC(LOW,IAXIS),blkLimitsGC(HIGH,IAXIS) faceyData(MAG_FACE_VAR,i,j,k)= sim_magy enddo enddo endif #endif call Grid_releaseBlkPtr(blockID,solnData,CENTER) #if NFACE_VARS > 0 if (sim_killdivb) then call Grid_releaseBlkPtr(blockID,facexData,FACEX) call Grid_releaseBlkPtr(blockID,faceyData,FACEY) if (NDIM>2) call Grid_releaseBlkPtr(blockID,facezData,FACEZ) endif #endif deallocate(xcent) deallocate(ycent) deallocate(zcent) return end subroutine Simulation_initBlock