!******************************************************************************* ! ! Subroutine: init_block() ! ! 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) ! ! Parameters: block_no The number of the block to initialize ! ! Runtime: 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: write_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 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 logical, save :: first_call = .true. 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 = 10000 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 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 ! *** durham *** real :: newtonG ! Newton const. G (FROM DB) real :: r200 ! virial radius (DERIVED) real :: conC ! concentration parameter (FROM GP) real :: eta ! gas profile exp. (FROM GP) real :: halo_m ! Halo mass (FROM GP) real :: km, kpc, mpc, sol, prot_mass, k_B, mu real, save :: rho_crit ! critical density (DERIVED) real, save :: alpha ! NFW charact. pot.(DERIVED) real, save :: spec_e ! specific energy (DERIVED) real, save :: r_s ! shape radius (DERIVED) real, save :: rho_s ! shape density (DERIVED) real, save :: delta_c ! shape parameter (DERIVED) ! *** durham *** !---------------------------------------------------------------------------- ! ! initialization ! !---------------------------------------------------------------------------- #ifdef DEBUG write(*,*) ' ** init_block.F90 <----- IN' #endif firstcall: if (first_call) then first_call = .false. MyPE = dBasePropertyInteger('MyProcessor') MasterPE = dBasePropertyInteger('MasterProcessor') ! ! 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') 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, '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, 't_init', t_init) 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) ! *** durham *** call get_constant_from_db ("Newton", newtonG) call get_constant_from_db ("proton mass", prot_mass) call get_constant_from_db ("Boltzmann", k_B) km = 1.E5 kpc = 3.0856775807E21 mpc = 3.0856775807E24 sol = 1.9889225E33 mu = 0.59 call get_parm_from_context(global_parm_context, 'halo_m', halo_m) call get_parm_from_context(global_parm_context, 'conC', conC) call get_parm_from_context(global_parm_context, 'eta', eta) ! SET UP INITIAL VALUES: rho_crit = 3.0E+04 / (8.0 * pi * newtonG) * (km/mpc)**2 if (MyPE .eq. MasterPE) write(*,*) "Critical density :", rho_crit if (MyPE .eq. MasterPE) write(*,*) " ", & rho_crit * mpc**3 / sol, "Msun/Mpc^3" r200 = 1.63E-02 * halo_m**(1.0/3.0) * kpc if (MyPE .eq. MasterPE) write(*,*) "Virial radius :", r200 if (MyPE .eq. MasterPE) write(*,*) " ", & r200 / mpc, "Mpc" r_s = r200/6.0 if (MyPE .eq. MasterPE) write(*,*) "Shape radius :", r_s if (MyPE .eq. MasterPE) write(*,*) " ", & r_s / mpc, "Mpc" delta_c = 200.0/3.0 * conC**3 / ( log(1.0 + conC) - conC/(1.0 + conC) ) if (MyPE .eq. MasterPE) write(*,*) "Delta_C :", delta_c ! *** ++++++ *** alpha = 4.0 * pi * newtonG * delta_c * rho_crit * r_s**2 if (MyPE .eq. MasterPE) write(*,*) "Caracteristic pot. :", alpha spec_e = k_B / (mu * prot_mass) if (MyPE .eq. MasterPE) write(*,*) "Temperature :", alpha / & (eta * spec_e) spec_e = alpha / eta ! END SET UP ! *** durham *** 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) 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 write_logfile("flash: Initializing for cluster stability problem.") write (*,*) 'flash: initializing for cluster stability 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 write_logfile & ! & ("warning: t_init>0: sedov solution currently broken!") ! 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_init .gt. 0.) then ! *** durham *** write(*,*) 'Initial time is not 0.0!' stop ! *** durham *** else do i = 1, n_profile ! *** durham *** if (r_prof(i) == 0.0) then rho_prof(i) = 1.0E-03 * rho_crit * exp(eta) else rho_prof(i) = 1.0E-03 * rho_crit * & ( 1.0 + r_prof(i)/r_s )**( eta * r_s/r_prof(i) ) end if p_prof(i) = rho_prof(i) * spec_e v_prof(i) = 0.0 ! *** durham *** enddo endif open(unit=1,file='profile.dat',status='unknown') write(1,*) dr_prof*1000. - 1.41562919156465981 write(1,*) p_exp - 681.16756574500903 do i = 1, n_profile write(1,'(i6,5es12.4)') & i, dr_prof*i, rho_prof(i), v_prof(i), p_prof(i), & rho_crit * delta_c / & ( (r_prof(i) + r_s/20.0)/r_s * (1.0 + r_prof(i)/r_s)**2 ) enddo close(1) endif firstcall !------------------------------------------------------------------------------- ! ! 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, smallx ) ! changed 10^-10 to smallx 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 #ifdef DEBUG write(*,*) ' ** init_block.F90 <----- OUT' #endif end subroutine init_block !=============================================================================== ! 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 end subroutine find