!****************************************************************************** ! ! Subroutine: init_block() ! ! Description: Initializes fluid data (density, pressure, velocity, etc.) for ! a specified block. This version sets up the RM 1gamma shock-tube ! problem. ! Shuang Zhang @ Vizlab, 10/02/02 ! Parameters: block_no The number of the block to initialize subroutine init_block(block_no) !============================================================================== use logfile, ONLY: write_logfile use dBase, ONLY: k2d, k3d, nxb, nyb, nzb, nguard, & ndim, ionmax, & dBaseKeyNumber, dBaseSpecies, & dBaseGetCoords, dBasePutData, & dBasePropertyInteger use runtime_parameters, ONLY: get_parm_from_context, & & GLOBAL_PARM_CONTEXT implicit none ! compute the maximum length of a vector in each coordinate direction ! (including guardcells) integer, parameter :: q = max(nxb+2*nguard, & & nyb+2*nguard*k2d, & & nzb+2*nguard*k3d) integer :: block_no integer :: i, j, k, n integer :: imax, jmax, kmax integer, save :: MyPE, MasterPE real :: xx, yy, zz, xxL, xxR real, save :: xcos, ycos, zcos real :: lposn0, lposn real, dimension(q) :: x, y, z, xl, xr integer, save :: idens, itemp, ipres, iener, igame, igamc integer, save :: ivelx, ively, ivelz, inuc_begin real, save :: pamb, vamb, damb real :: pshock, vshock, dshock,rr,dnew real :: ishock, icon1, ngs, sigma, gas1, gas2, pratio real, save :: mach,radius,trans,dratio,shock,& & xmax, ymax,zmax,x_cen, y_cen, z_cen real, external :: erfunc real :: rho_zone, velx_zone, vely_zone, velz_zone, pres_zone, & & ener_zone, ekin_zone integer, save :: iXvector, iYvector, iZvector integer, save :: iXcoord, iYcoord, iZcoord integer, save :: iPt integer, save :: izn, iznl, iznr real, save :: smallp,smallx real, save :: gamma logical, save :: firstCall = .TRUE. call write_logfile & & ("flash: initializing for 1gamma rm problem.") !============================================================================== if (firstCall) then MyPE = dBasePropertyInteger('MyProcessor') MasterPE = dBasePropertyInteger('MasterProcessor') ! grab pointers into the database, so we don't need expensive string compares idens = dBaseKeyNumber('dens') ivelx = dBaseKeyNumber('velx') ively = dBaseKeyNumber('vely') ivelz = dBaseKeyNumber('velz') iener = dBaseKeyNumber('ener') ipres = dBaseKeyNumber('pres') itemp = dBaseKeyNumber('temp') igame = dBaseKeyNumber('game') igamc = dBaseKeyNumber('gamc') inuc_begin = dBaseSpecies(1) iXcoord = dBaseKeyNumber('xCoord') iYcoord = dBaseKeyNumber('yCoord') iZcoord = dBaseKeyNumber('zCoord') iPt = dBaseKeyNumber('Point') izn = dBaseKeyNumber('zn') iznl = dBaseKeyNumber('znl') iznr = dBaseKeyNumber('znr') ! get any runtime parameters we need for this setup call get_parm_from_context(GLOBAL_PARM_CONTEXT, & & 'smallp', smallp) call get_parm_from_context(GLOBAL_PARM_CONTEXT, & & 'smallx', smallx) call get_parm_from_context(GLOBAL_PARM_CONTEXT, & & 'gamma', gamma) call get_parm_from_context(GLOBAL_PARM_CONTEXT,'mach', mach ) call get_parm_from_context(GLOBAL_PARM_CONTEXT,'shock', shock ) call get_parm_from_context(GLOBAL_PARM_CONTEXT,'damb', damb ) call get_parm_from_context(GLOBAL_PARM_CONTEXT,'vamb', vamb ) call get_parm_from_context(GLOBAL_PARM_CONTEXT,'radius', radius ) call get_parm_from_context(GLOBAL_PARM_CONTEXT,'trans', trans ) call get_parm_from_context(GLOBAL_PARM_CONTEXT,'pamb', pamb ) call get_parm_from_context(GLOBAL_PARM_CONTEXT,'xmax', xmax ) call get_parm_from_context(GLOBAL_PARM_CONTEXT,'ymax', ymax ) call get_parm_from_context(GLOBAL_PARM_CONTEXT,'zmax', zmax ) call get_parm_from_context(GLOBAL_PARM_CONTEXT,'x_cen', x_cen ) call get_parm_from_context(GLOBAL_PARM_CONTEXT,'y_cen', y_cen ) call get_parm_from_context(GLOBAL_PARM_CONTEXT,'z_cen', z_cen ) call get_parm_from_context(GLOBAL_PARM_CONTEXT,'dratio', dratio ) call get_parm_from_context(GLOBAL_PARM_CONTEXT,'gamma', gamma ) ! dump some output to stdout listing the paramters if (MyPE == MasterPE) then call write_logfile & & ("flash: initializing for 1gamma rm problem.") endif firstCall = .FALSE. endif ngs = trans*ymax ishock = shock*xmax icon1 = ishock+2*ngs+0.01*ymax dnew = dratio*damb gas1 = (gamma+1.)/(gamma-1.) gas2 = 2.0*gamma/(gamma+1.) pratio = 1.+gas2*(mach*mach-1.) pshock = pamb * pratio dshock = damb * (1.+gas1*pratio)/(gas1+pratio) vshock = sqrt(pamb/damb/gamma)*(pratio-1.) & & *sqrt(gas2/(pratio+1./gas1)) sigma = ngs/1. imax = 2*nguard + nxb ! Maximum values of the cell index ranges jmax = 2*nguard*k2d + nyb kmax = 2*nguard*k3d + nzb ! get the coordinate information for the current block from the database x(:) = 0.0 y(:) = 0.0 z(:) = 0.0 if (ndim == 3) call dBaseGetCoords(izn, iZcoord, block_no, z) if (ndim >= 2) call dBaseGetCoords(izn, iYcoord, block_no, y) call dBaseGetCoords(izn, iXcoord, block_no, x) !------------------------------------------------------------------------------ ! Loop over cells in the block. For each, compute the physical position of ! its left and right edge and its center as well as its physical width. ! Then decide which side of the initial discontinuity it is on and initialize ! the hydro variables appropriately. ! mz: I don't know why we initialize the guardcells here. do k = 1, kmax do j = 1, jmax do i = 1, imax if (x(i) .lt. ishock) then ! in shock region rho_zone = dshock pres_zone = pshock velx_zone = vshock-vamb*vshock vely_zone = 0. velz_zone = 0. else rho_zone = damb vely_zone = 0. velz_zone = 0. rr = (x(i)-(icon1+radius*ymax))**2+ & ! Didn't use x_cen, depend on shock position & (y(j)-y_cen)**2 &! Cylinder & +(z(k)-z_cen)**2 ! Sphere if (rr.le.(radius*ymax+2*ngs)**2.and. & & rr.gt.(radius*ymax+ngs)**2) then rho_zone=((dnew+damb)/2.)* & & (1.-(dnew-damb)/(dnew+damb)* & & erfunc(sqrt(ABS(rr-(radius*ymax+ngs)**2))/sigma,0.001)) else if(rr.le.(radius*ymax+ngs)**2.and. & & rr.gt.(radius*ymax)**2) then rho_zone = ((dnew+damb)/2.)* & & (1.+(dnew-damb)/(dnew+damb)* & & erfunc(sqrt(ABS(rr-(radius*ymax+ngs)**2))/sigma,0.001)) else if (rr.le.(radius*ymax)**2) then rho_zone = dnew endif pres_zone = pamb velx_zone = -vamb*vshock endif call dBasePutData(inuc_begin,iPt, & & i,j,k, block_no, 1.0e0-(n-1)*smallx) do n = 2, ionmax call dBasePutData(inuc_begin-1+n,iPt, & & i,j,k, block_no, smallx) enddo ! Compute the gas energy and set the gamma-values needed for the equation of ! state. ekin_zone = 0.5 * (velx_zone**2 + & & vely_zone**2 + & & velz_zone**2) ener_zone = pres_zone / (gamma-1.) ener_zone = ener_zone / rho_zone ener_zone = ener_zone + ekin_zone ener_zone = max(ener_zone, smallp) ! store the variables in the current zone via the database put methods call dBasePutData(idens, iPt, i,j,k, block_no, rho_zone) call dBasePutData(ipres, iPt, i,j,k, block_no, pres_zone) call dBasePutData(iener, iPt, i,j,k, block_no, ener_zone) call dBasePutData(igamc, iPt, i,j,k, block_no, gamma) call dBasePutData(igame, iPt, i,j,k, block_no, gamma) call dBasePutData(ivelx, iPt, i,j,k, block_no, velx_zone) call dBasePutData(ively, iPt, i,j,k, block_no, vely_zone) call dBasePutData(ivelz, iPt, i,j,k, block_no, velz_zone) enddo enddo enddo !============================================================================== return end ! Error function definition real function erfunc(xx,dxx) real xx,dxx real xnn integer nn,ie xnn = xx/dxx nn = int(xnn) erfunc=1.0+exp(-xx**2) do ie = 1, nn-1 erfunc = erfunc + 2.*exp(-(dxx*ie)**2) enddo erfunc=(2./sqrt(3.1415926))*erfunc*dxx/2.0 return end