!!****f* setups/sedov/init_block !! !! NAME !! !! init_block !! !! !! SYNOPSIS !! !! init_block(block_no) !! !! init_block(integer) !! !! !! DESCRIPTION !! !! Initializes fluid data (density, pressure, velocity, etc.) for !! a specified block. This version sets up the Sedov spherical !! explosion problem. !! !! References: Sedov, L. I., 1959, Similarity and Dimensional Methods !! in Mechanics (New York: Academic) !! !! Landau, L. D., & Lifshitz, E. M., 1987, Fluid Mechanics, !! 2d ed. (Oxford: Pergamon) !! !! ARGUMENTS !! !! block_no The number of the block to initialize !! !! !! PARAMETERS !! !! p_ambient Initial ambient pressure !! !! rho_ambient Initial ambient density !! !! exp_energy Explosion energy (initially all internal) !! !! r_init Radius of initial high-pressure region !! !! x/y/zctr Coordinates of the center of the explosion !! !! !!*** subroutine init_block (block_no) use logfile, ONLY: stamp_logfile use physical_constants, ONLY: get_constant_from_db use runtime_parameters use multifluid_database, ONLY: find_fluid_index use dBase, ONLY: nxb, nyb, nzb, nguard, ionmax, & MAXCELLS, k2d, k3d, ndim, & dBasePropertyInteger, & dBaseKeyNumber, dBaseSpecies, & dBaseGetCoords, dBasePutData, & CARTESIAN, CYLINDRICAL implicit none integer :: block_no integer :: i, j, k, n, jlo, jhi integer :: Nint, ii, jj, kk real :: delx, dely, delz, distinv real :: xdist, ydist, zdist real :: sum_p, vol, vol_amb, vol_pert, vol_sub integer, parameter :: q = MAXCELLS real, DIMENSION(q) :: x, xl, xr, y, yl, yr, z, zl, zr, vx, vy, vz, p, rho real, DIMENSION(q) :: e, ek, gam real, save, DIMENSION(ionmax) :: xn real :: xmin, xmax, ymin, ymax, zmin, zmax real :: xx, dxx, yy, dyy, zz, dzz, frac, xxl, xxr ! ! Parameters needed by the initial model ! real, save :: pi, p_ambient, rho_ambient, exp_energy, r_init real, save :: gamma, xctr, yctr, zctr, p_exp, vctr real, save :: smallx, smlrho, smallp integer, save :: iXvector integer, save :: iXcoord, iYcoord, iZcoord integer, save :: iPoint integer, save :: izn, iznl, iznr real :: dist integer, save :: idens, ipres, iener, igame, igamc integer, save :: ivelx, ively, ivelz, inuc_begin, ifuel integer, save :: MyPE, MasterPE integer :: meshGeom logical :: validGeom integer, save :: nsubzones logical, save :: first_call = .true. !------------------------------------------------------------------------------ ! ! initialization ! !------------------------------------------------------------------------------ if (first_call) then first_call = .false. MyPE = dBasePropertyInteger('MyProcessor') MasterPE = dBasePropertyInteger('MasterProcessor') meshGeom = dBasePropertyInteger("MeshGeometry") validGeom = ( meshGeom == CARTESIAN .OR. & (meshGeom == CYLINDRICAL .AND. ndim == 2)) if (.NOT. validGeom) then print *, "ERROR: sedov only works for Cartesian or 2-d cylindrical geometries" call abort_flash("ERROR: sedov only works for Cartesian or 2-d cylindrical geometries") endif ! ! Grab parameters from databases ! idens = dBaseKeyNumber('dens') ivelx = dBaseKeyNumber('velx') ively = dBaseKeyNumber('vely') ivelz = dBaseKeyNumber('velz') iener = dBaseKeyNumber('ener') ipres = dBaseKeyNumber('pres') igame = dBaseKeyNumber('game') igamc = dBaseKeyNumber('gamc') if (idens < 0) call abort_flash("init_block: no key dens") if (ivelx < 0) call abort_flash("init_block: no key velx") if (ively < 0) call abort_flash("init_block: no key vely") if (ivelz < 0) call abort_flash("init_block: no key velz") if (iener < 0) call abort_flash("init_block: no key ener") if (ipres < 0) call abort_flash("init_block: no key pres") if (igame < 0) call abort_flash("init_block: no key game") if (igamc < 0) call abort_flash("init_block: no key gamc") inuc_begin = dBaseSpecies(1) ifuel = 1 iXvector = dBaseKeyNumber('xVector') iPoint = dBaseKeyNumber('Point') iXcoord = dBaseKeyNumber('xCoord') iYcoord = dBaseKeyNumber('yCoord') iZcoord = dBaseKeyNumber('zCoord') izn = dBaseKeyNumber('zn') iznl = dBaseKeyNumber('znl') iznr = dBaseKeyNumber('znr') call get_constant_from_db ("pi", pi) call get_parm_from_context(global_parm_context, 'p_ambient', p_ambient) call get_parm_from_context(global_parm_context, & 'rho_ambient', rho_ambient) call get_parm_from_context(global_parm_context, 'exp_energy', exp_energy) call get_parm_from_context(global_parm_context, 'r_init', r_init) call get_parm_from_context(global_parm_context, 'gamma', gamma) call get_parm_from_context(global_parm_context, 'smallx', smallx) call get_parm_from_context(global_parm_context, 'smlrho', smlrho) call get_parm_from_context(global_parm_context, 'smallp', smallp) call get_parm_from_context(global_parm_context,'xctr',xctr) call get_parm_from_context(global_parm_context,'yctr',yctr) call get_parm_from_context(global_parm_context,'zctr',zctr) xn(:) = smallx xn(ifuel) = 1.0 - (ionmax-1)*smallx call get_parm_from_context(global_parm_context,'xmin',xmin) call get_parm_from_context(global_parm_context,'ymin',ymin) call get_parm_from_context(global_parm_context,'zmin',zmin) 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, & 'nsubzones',nsubzones) ! ! Calculate the initial volume and interior pressure. ! if (ndim .eq. 1) then vctr = 2. * r_init elseif (ndim .eq. 2) then if (meshGeom == CARTESIAN) then vctr = pi * r_init**2 else if (meshGeom == CYLINDRICAL) then vctr = (4./3.)*pi*r_init**3 endif else vctr = 4./3.*pi*r_init**3 endif p_exp = (gamma-1.) * exp_energy / vctr ! ! Write a message to stdout describing the problem setup. ! if (MyPE .eq. MasterPE) then write (*,*) call stamp_logfile("initializing for sedov problem", "run_init") write (*,*) 'flash: initializing for sedov problem.' write (*,*) write (*,*) 'p_ambient = ', p_ambient write (*,*) 'rho_ambient= ', rho_ambient write (*,*) 'gamma = ', gamma write (*,*) 'exp_energy = ', exp_energy write (*,*) 'r_init = ', r_init write (*,*) 'p_exp = ', p_exp write (*,*) 'xctr = ', xctr write (*,*) 'yctr = ', yctr write (*,*) 'zctr = ', zctr write (*,*) 'ndim = ', ndim write (*,*) 'nsubzones = ', nsubzones write (*,*) endif endif x(:) = 0.0 xl(:) = 0.0 xr(:) = 0.0 y(:) = 0.0 yl(:) = 0.0 yr(:) = 0.0 z(:) = 0.0 zl(:) = 0.0 zr(:) = 0.0 if (ndim == 3) then call dBaseGetCoords(izn, iZcoord, block_no, z) call dBaseGetCoords(iznl, iZcoord, block_no, zl) call dBaseGetCoords(iznr, iZcoord, block_no, zr) endif if (ndim >= 2) then call dBaseGetCoords(izn, iYcoord, block_no, y) call dBaseGetCoords(iznl, iYcoord, block_no, yl) call dBaseGetCoords(iznr, iYcoord, block_no, yr) endif call dBaseGetCoords(izn, iXcoord, block_no, x) call dBaseGetCoords(iznl, iXcoord, block_no, xl) call dBaseGetCoords(iznr, iXcoord, block_no, xr) ! ! For each cell ! do k = nguard*k3d+1, nguard*k3d+nzb dzz = zr(k) - zl(k) do j = nguard*k2d+1, nguard*k2d+nyb dyy = yr(j) - yl(j) do i = nguard+1, nguard+nxb dxx = xr(i) - xl(i) sum_p = 0.0 vol = 0.0 vol_amb = 0.0 vol_pert = 0.0 ! Break the cell into nsubzones^ndim sub-zones to more accurately ! initialize the perturbation do kk = 0, (nsubzones-1)*k3d zz = zl(k) + (dzz/nsubzones)*(kk + 0.5) zdist = (zz - zctr) * k3d do jj = 0, (nsubzones-1)*k2d yy = yl(j) + (dyy/nsubzones)*(jj + 0.5) ydist = (yy - yctr) * k2d do ii = 0, (nsubzones-1) xxl = xl(i) + (dxx/nsubzones)*ii xxr = xl(i) + (dxx/nsubzones)*(ii + 1.0) xx = 0.5*(xxl + xxr) xdist = xx - xctr dist = sqrt( xdist**2 + ydist**2 + zdist**2 ) if (meshGeom == CARTESIAN) then vol_sub = 1.0 else if (meshGeom == CYLINDRICAL) then vol_sub = pi*(xxr**2 - xxl**2) endif if (dist <= r_init) then vol_pert = vol_pert + vol_sub else vol_amb = vol_amb + vol_sub endif vol = vol + vol_sub enddo enddo enddo p (i) = max(((vol_pert*p_exp) + & (vol_amb*p_ambient))/vol, smallp) rho (i) = rho_ambient vx (i) = 0.0 vy (i) = 0.0 vz (i) = 0.0 ek (i) = 0.5*(vx(i)*vx(i) + vy(i)*vy(i) + vz(i)*vz(i)) ! ! assume gamma-law equation of state ! gam (i) = gamma e (i) = p(i)/(gam(i)-1.) e (i) = e(i)/rho(i) + ek(i) e (i) = max (e(i), smallp) do n=1,ionmax call dBasePutData(inuc_begin-1+n,ipoint, & i, j, k, block_no, xn(n), 1) enddo enddo call dBasePutData(idens, iXvector, j, k, block_no, rho, 1) call dBasePutData(ipres, iXvector, j, k, block_no, p, 1 ) call dBasePutData(ivelx, iXvector, j, k, block_no, vx, 1 ) call dBasePutData(ively, iXvector, j, k, block_no, vy, 1 ) call dBasePutData(ivelz, iXvector, j, k, block_no, vz, 1 ) call dBasePutData(igame, iXvector, j, k, block_no, gam, 1) call dBasePutData(igamc, iXvector, j, k, block_no, gam, 1) call dBasePutData(iener, iXvector, j, k, block_no, e, 1) enddo enddo return end subroutine init_block