!!****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 !! !! t_init Initial time (use to generate analytical !! solution if >0) !! !!*** 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_rho, sum_p, sum_vx, sum_vy, sum_vz real :: vel, pres, dens integer, parameter :: q = MAXCELLS real, DIMENSION(q) :: x, y, z, vx, vy, vz, p, rho real, DIMENSION(q) :: e, ek, gam real, save, DIMENSION(ionmax) :: xn ! ! A calculated 1-d initial model ! integer, parameter :: n_profile = 1000 integer, save :: nsubzones real , save :: insubzones, insubzm1 real , save :: inszd real, dimension(n_profile), save :: r_prof, rho_prof, p_prof real, dimension(n_profile), save :: v_prof real, save :: dr_prof character(len=80), save :: prof_file real :: xmin, xmax, ymin, ymax, zmin, zmax, diagonal real :: xx, dxx, yy, dyy, zz, dzz, frac ! ! Parameters needed by the initial model ! real, save :: pi, p_ambient, rho_ambient, exp_energy, r_init real, save :: t_init, 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 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 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') 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,'prof_file',prof_file) diagonal = (xmax-xmin)**2 diagonal = diagonal + k2d*(ymax-ymin)**2 diagonal = diagonal + k3d*(zmax-zmin)**2 diagonal = sqrt(diagonal) call get_parm_from_context(global_parm_context, & 'nsubzones',nsubzones) ! For now set t_init=0.0 t_init=0.0 if (nsubzones .le. 1) nsubzones = 2 insubzones = 1./real(nsubzones) insubzm1 = 1./real(nsubzones-1) inszd = insubzones**ndim ! ! Calculate the initial volume and interior pressure. ! if (ndim .eq. 1) then vctr = 2. * r_init elseif (ndim .eq. 2) then vctr = pi * r_init**2 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 (*,*) if (t_init .gt. 0.) then call stamp_logfile & ('init_block: t_init > 0, sedov solution currently broken', "warning") endif endif ! ! Construct the radial samples needed for the initialization. ! dr_prof = diagonal / (n_profile-1) do i = 1, n_profile r_prof(i) = (i-1) * dr_prof enddo ! ! If t>0, use the analytic Sedov solution to initialize the ! code. Otherwise, just use a top-hat. ! if (t_init .gt. 0.) then call set_analytic_sedov (n_profile, r_prof, rho_prof, p_prof, & v_prof, t_init, gamma, exp_energy, & p_ambient, rho_ambient) else do i = 1, n_profile rho_prof(i) = rho_ambient p_prof(i) = p_ambient v_prof(i) = 0. if (r_prof(i) .le. r_init) p_prof(i) = p_exp enddo endif open(unit=1,file=prof_file,status='unknown') write(1,*) dr_prof*1000. - 1.41562919156465981 write(1,*) p_exp - 681.16756574500903 do i = 1, n_profile write(1,*) i, dr_prof*i, rho_prof(i),v_prof(i),p_prof(i) enddo close(1) endif !------------------------------------------------------------------------------- ! ! End of initialization. We now have a one-d profile and need only ! apply it. ! !------------------------------------------------------------------------------- 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) ! ! For each cell ! do k = 1, 2*nguard*k3d+nzb if (k .eq. 1) then dzz = z(2) - z(1) else dzz = z(k) - z(k-1) endif zz = z(k) do j = 1, 2*nguard*k2d+nyb if (j .eq. 1) then dyy = y(2) - y(1) else dyy = y(j) - y(j-1) endif yy = y(j) do i = 1, 2*nguard+nxb xx = x(i) if (i .eq. 1) then dxx = x(2) - x(1) else dxx = x(i) - x(i-1) endif sum_rho = 0. sum_p = 0. sum_vx = 0. sum_vy = 0. sum_vz = 0. ! ! Break the cell into nsubzones^ndim sub-zones, and look up the ! appropriate quantities along the 1d profile for that subzone. ! ! Have the final values for the zone be equal to the average of ! the subzone values. ! do kk = 0, (nsubzones-1)*k3d zz = z(k) + (kk*insubzm1-.5)*dzz zdist = (zz - zctr) * k3d do jj = 0, (nsubzones-1)*k2d yy = y(j) + (jj*insubzm1-.5)*dyy ydist = (yy - yctr) * k2d do ii = 0, (nsubzones-1) xx = x(i) + (ii*insubzm1-.5)*dxx xdist = xx - xctr dist = sqrt( xdist**2 + ydist**2 + zdist**2 ) distinv = 1. / max( dist, 1.E-10 ) call find (r_prof, n_profile, dist, jlo) ! ! a point at `dist' is frac-way between jlo and jhi. We do a ! linear interpolation of the quantities at jlo and jhi and sum those. ! if (jlo .eq. 0) then jlo = 1 jhi = 1 frac = 0. else if (jlo .eq. n_profile) then jlo = n_profile jhi = n_profile frac = 0. else jhi = jlo + 1 frac = (dist - r_prof(jlo)) / & (r_prof(jhi)-r_prof(jlo)) endif ! ! Now total these quantities. Note that v is a radial velocity; ! we multiply by the tangents of the appropriate angles to get ! the projections in the x, y and z directions. ! sum_p = sum_p + & p_prof(jlo) + frac*(p_prof(jhi) - p_prof(jlo)) sum_rho = sum_rho + & rho_prof(jlo) + frac*(rho_prof(jhi)- rho_prof(jlo)) vel = v_prof(jlo) + frac*(v_prof(jhi) - v_prof(jlo)) sum_vx = sum_vx + vel*xdist*distinv sum_vy = sum_vy + vel*ydist*distinv sum_vz = sum_vz + vel*zdist*distinv enddo enddo enddo p (i) = max(sum_p * inszd, smallp) rho (i) = max(sum_rho * inszd, smlrho) vx (i) = sum_vx * inszd vy (i) = sum_vy * inszd vz (i) = sum_vz * inszd 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 ! Subroutine: set_analytic_sedov() ! Description: Given a set of arrays to store radius, density, pressure, and ! velocity, together with a time and a ratio of specific heats, ! generate the analytical solution to the Sedov problem. ! Currently this routine is limited to gamma=1.4 (I've hardwired ! the dimensionless constant beta) and will complain if it gets ! anything else. subroutine set_analytic_sedov (N, r, rho, p, v, t, gamma, E, & p_ambient, rho_ambient) !============================================================================== integer N real r(N), rho(N), p(N), v(N), t, gamma, E, p_ambient, rho_ambient real beta, R0, dr, xi, nu1, nu2, nu3, nu4, nu5, VV, G, Z, & kappa, zeta, epsilon, c2sqr, k, gamp1, gam7, gamm1 integer i !============================================================================== if (gamma .ne. 1.4) then write (*,*) 'Warning! init_block() found gamma<>1.4 and t>0.' write (*,*) ' Analytical initial conditions will be' write (*,*) ' wrong. Assuming beta=1.033...' endif ! Compute dimensionless scaling constant and explosion radius. beta = 1.033 R0 = beta * (E*t*t/rho_ambient)**0.2 ! Compute exponents for self-similar solution. nu1 = - (13.*gamma*gamma - 7.*gamma + 12.) / & ((3.*gamma - 1.) * (2.*gamma+1.)) nu2 = 5. * (gamma - 1.) / (2.*gamma + 1.) nu3 = 3. / (2.*gamma + 1.) nu4 = - nu1 / (2. - gamma) nu5 = - 2. / (2. - gamma) ! Other useful combinations of gamma. gamp1 = gamma + 1.E0 gamm1 = gamma - 1.E0 gam7 = 7.E0 - gamma k = gamp1 / gamm1 !============================================================================== ! Generate the solution. do i = 1, N xi = r(i) / R0 ! Fraction of explosion radius. if (xi .le. 1.) then ! Inside the explosion. ! Compute V(xi) using bisection. ! See Landau & Lifshitz for notation. call compute_sedov_v (xi, gamma, nu1, nu2, VV) G = k * (k*(gamma*VV-1.))**nu3 * & (gamp1/gam7*(5.-(3.*gamma-1.)*VV))**nu4 * & (k*(1.-VV))**nu5 rho(i) = rho_ambient * G v(i) = 2.*r(i)*VV / (5.*t) if (xi .le. 1.E-6) then ! Use asymptotic r->0 solution. kappa = ( (0.5*gamp1/gamma)**2. * & (gamp1/gam7* & (5.-(3.*gamma-1.)/gamma))**nu1 )**(1./nu2) * & gamm1/gamp1 / gamma epsilon = k**(nu5+1.) * (k*gamma*kappa)**nu3 * & (gamp1/gam7*(3.*gamma-1.))**nu4 * & ((2.*gamma+1)/gamma/(3.*gamma-1.)) * & (gamm1/gamma) zeta = gamm1*gamm1/(2.*gamma*gamma*kappa) p(i) = rho_ambient/gamma * 0.16*(R0/t)**2 * epsilon*zeta else Z = gamma * gamm1 * (1.-VV) * VV**2 / (2.*(gamma*VV-1.)) c2sqr = 0.16*(r(i)/t)**2 * Z p(i) = rho(i) * c2sqr / gamma endif else ! Outside the explosion. rho(i) = rho_ambient p(i) = p_ambient v(i) = 0. endif enddo return end subroutine set_analytic_sedov ! Subroutine: compute_sedov_v() ! Description: Compute the dimensionless velocity function V(xi) in the ! Sedov problem using bisection. See Landau & Lifshitz for ! notation. subroutine compute_sedov_v (xi, gamma, nu1, nu2, V) !============================================================================== real xi, gamma, nu1, nu2, V real VL, VR, xiL, xiR, Vmid, ximid, sedov_vfunc, tolerance, logxi integer n_iter, n_iter_max parameter (n_iter_max = 500, tolerance = 1.E-6) !============================================================================== if (xi .le. 1.E-6) then ! Use asymptotic xi->0 solution. V = 1./gamma else ! Do bisection root search. logxi = alog(xi) VL = 1./gamma VR = 2./gamma xiL = sedov_vfunc(VL, gamma, nu1, nu2) xiR = sedov_vfunc(VR, gamma, nu1, nu2) n_iter = 1 10 Vmid = 0.5 * (VL + VR) ximid = sedov_vfunc(Vmid, gamma, nu1, nu2) if ((abs(ximid - logxi) .le. tolerance) .or. & (n_iter .gt. n_iter_max)) goto 20 n_iter = n_iter + 1 if (ximid .gt. logxi) then VR = Vmid else VL = Vmid endif goto 10 20 if (n_iter .gt. n_iter_max) & write (*,*) 'compute_sedov_v: did not reach ', & 'max precision for xi = ', xi V = Vmid endif return end subroutine compute_sedov_v !****************************************************************************** ! Function: sedov_vfunc() ! Description: Function to use in bisection root search (compute_sedov_v()). real function sedov_vfunc (V, gamma, nu1, nu2) real V, gamma, nu1, nu2 real gamp1, gamm1, gam7, k, xi gamp1 = gamma + 1. gamm1 = gamma - 1. gam7 = 7. - gamma k = gamp1 / gamm1 xi = nu1*alog(5.-(3.*gamma-1.)*V) + & nu2*alog(gamma*V-1.) - & nu1*alog(gam7/gamp1) - nu2*alog(gamm1/gamp1) - & 2.*alog(0.5*gamp1) sedov_vfunc = 0.2 * xi return end function sedov_vfunc !****************************************************************************** ! Routine: find() ! Description: Given a monotonically increasing table x(N) and a test value ! x0, return the index i of the largest table value less than ! or equal to x0 (or 0 if x0 < x(1)). Use binary search. subroutine find (x, N, x0, i) integer N, i, il, ir, im real x(N), x0 if (x0 .lt. x(1)) then i = 0 elseif (x0 .gt. x(N)) then i = N else il = 1 ir = N 10 if (ir .eq. il+1) goto 20 im = (il + ir) / 2 if (x(im) .gt. x0) then ir = im else il = im endif goto 10 20 i = il endif return end subroutine find