!!****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(:,:,:,:) :: facexData,faceyData #if NDIM > 0 real, pointer, dimension(:,:,:,:) :: facezData #endif #ifndef CHAM_SPEC integer :: CHAM_SPEC = 1, PLAS_SPEC = 2, GOLD_SPEC = 3, PART_SPEC = 4, TARG_SPEC = 5 #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)) #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 = blkLimits(LOW,KAXIS),blkLimits(HIGH,KAXIS) do j = blkLimits(LOW,JAXIS),blkLimits(HIGH,JAXIS) do i = blkLimits(LOW,IAXIS),blkLimits(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_goldHeight + sim_plasticHeight + sim_vacuumHeight .and. & xcent(i) >= sim_vacuumHeight + sim_plasticHeight + sim_goldHeight .and. & (sqrt((xcent(i) - sim_particleXcenter)**2 + (ycent(j) - sim_particleYcenter)**2) > sim_particleRadius)) then species = TARG_SPEC elseif ( xcent(i) <= sim_goldHeight + sim_plasticHeight + sim_vacuumHeight .and. & xcent(i) >= sim_vacuumHeight + sim_plasticHeight ) then species = GOLD_SPEC elseif ( xcent(i) <= sim_plasticHeight + sim_vacuumHeight .and. & xcent(i) >= sim_vacuumHeight ) then species = PLAS_SPEC end if elseif(NDIM == 2 .or. NDIM == 3) then if ( xcent(i) <= sim_targetRadius .and. & ycent(j) <= sim_targetHeight + sim_plasticHeight + sim_goldHeight + sim_vacuumHeight .and. & ycent(j) >= sim_vacuumHeight + sim_plasticHeight + sim_goldHeight .and. & (sqrt((xcent(i) - sim_particleXcenter)**2 + (ycent(j) - sim_particleYcenter)**2) > sim_particleRadius)) then species = TARG_SPEC elseif ( xcent(i) <= sim_targetRadius .and. & ycent(j) <= sim_vacuumHeight + sim_plasticHeight + sim_goldHeight .and. & ycent(j) >= sim_vacuumHeight + sim_plasticHeight ) then species = GOLD_SPEC elseif (xcent(i) <= sim_targetRadius .and. & ycent(j) <= sim_plasticHeight + sim_vacuumHeight .and. & ycent(j) >= sim_vacuumHeight ) then species = PLAS_SPEC elseif ((sqrt((xcent(i) - sim_particleXcenter)**2 + (ycent(j) - sim_particleYcenter)**2))<= sim_particleRadius) then species = PART_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 == PART_SPEC) then rho = sim_rhoPart tele = sim_telePart tion = sim_tionPart trad = sim_tradPart elseif (species == GOLD_SPEC) then rho = sim_rhoGold tele = sim_teleGold tion = sim_tionGold trad = sim_tradGold elseif (species == PLAS_SPEC) then rho = sim_rhoPlas tele = sim_telePlas tion = sim_tionPlas trad = sim_tradPlas elseif (species == CHAM_SPEC) then rho = sim_rhoCham tele = sim_teleCham tion = sim_tionCham trad = sim_tradCham else rho = sim_rhoTarg tele = sim_teleTarg tion = sim_tionTarg trad = sim_tradTarg 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) #endif #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)= 0.0 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 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