!******************************************************************************* ! Subroutine: init_block() ! Description: Initializes fluid data (density, pressure, velocity, etc.) for ! a specified block. This version sets up the dust cloud ! collapse problem. ! References: Colgate, S. A., & White, R. H. 1966, ApJ, 143, 626 ! Monchmeyer, R., & Muller, E. 1989, A&A, 217, 351 ! Parameters: block_no The number of the block to initialize ! Runtime: R_init Initial radius of cloud ! rho_0 Initial density of cloud ! T_ambient Initial ambient temperature (everywhere) ! x/y/zctr Coordinates of the center of the cloud subroutine init_block (block_no) !=============================================================================== use logfile use physical_constants use runtime_parameters use dBase, ONLY: k2d, k3d, & & ionmax, nguard, & & dBasePropertyInteger, & & dBasePropertyReal, dBaseKeyNumber, & & dBaseBlockSize, dBaseBlockCoord, & & dBaseGetDataPtrAllBlocks, & & dBaseSpecies implicit none real, pointer, dimension(:,:,:,:,:) :: data integer, PARAMETER :: MAXDIMS = 3 logical, save :: first_call = .true. integer, save :: nxb, nyb, nzb integer, save :: ivelx, ively, ivelz, itemp integer, save :: ipres, iener, idens, igame, igamc integer, save :: inuc_begin integer, save :: ndim, mype, masterpe real, save :: xmin, xmax, ymin, ymax, zmin, zmax real, save :: smallp, r_init, gamma, t_ambient real, save :: rho_0, smlrho, xctr, yctr, zctr real :: size(MAXDIMS), coord(MAXDIMS) integer block_no, i, j, k, n, imax, jmax, kmax, jlo integer Nint, ii, jj, kk real delx, xx, dely, yy, delz, zz, velocity, distinv real vfrac, xdist, ydist, zdist, dist, vctr real Nintinv, sum_rho, sum_p, sum_vx, sum_vy, sum_vz, & & Nintinv1 real xxmin, xxmax, yymin, yymax, zzmin, zzmax real pi, gascon, Newton, m_p, boltz real sndspd, presfac, ek integer N_ext integer N_prof parameter (N_prof = 10000) real r_prof(N_prof), rho_prof(N_prof), p_prof(N_prof) real v_prof(N_prof), dr_prof save r_prof, rho_prof, p_prof, v_prof save sndspd, presfac real :: l, dldr, dmdr, r, m, hdr, const, mconst, deltar, lpred, mpred, B logical :: iso_set = .true. integer :: nsubstep !=============================================================================== ! Initialize scalar quantities we will need. data => dBaseGetDataPtrAllBlocks() ! get soln vector call get_constant_from_db ("pi", pi) call get_constant_from_db ("ideal gas constant", gascon) call get_constant_from_db ("Newton", Newton) call get_constant_from_db ("proton mass", m_p) call get_constant_from_db ("Boltzmann", boltz) if (first_call) then call get_parm_from_context("xmin", xmin) call get_parm_from_context("xmax", xmax) call get_parm_from_context("ymin", ymin) call get_parm_from_context("ymax", ymax) call get_parm_from_context("zmax", zmax) call get_parm_from_context("zmin", zmin) call get_parm_from_context("smallp", smallp) call get_parm_from_context("r_init", r_init) call get_parm_from_context("gamma", gamma) call get_parm_from_context("t_ambient", t_ambient) call get_parm_from_context("rho_0", rho_0) call get_parm_from_context("smlrho", smlrho) call get_parm_from_context("xctr", xctr) call get_parm_from_context("yctr", yctr) call get_parm_from_context("zctr", zctr) nxb = dBasePropertyInteger("xBlockSize") nyb = dBasePropertyInteger("yBlockSize") nzb = dBasePropertyInteger("zBlockSize") myPE = dBasePropertyInteger("MyProcessor") masterPE = dBasePropertyInteger("MasterProcessor") ndim = dBasePropertyInteger("Dimensionality") ivelx = dBaseKeyNumber("velx") ively = dBaseKeyNumber("vely") ivelz = dBaseKeyNumber("velz") idens = dBaseKeyNumber("dens") ipres = dBaseKeyNumber("pres") itemp = dBaseKeyNumber("temp") igame = dBaseKeyNumber("game") igamc = dBaseKeyNumber("gamc") iener = dBaseKeyNumber("ener") inuc_begin= dBaseSpecies(1) ! first call set to false later ... end if imax = 2*nguard + nxb ! Maximum values of the cell index ranges jmax = 2*nguard*k2d + nyb kmax = 2*nguard*k3d + nzb ! Coordinates of the edges of the block size = dBaseBlockSize(block_no) coord = dBaseBlockCoord(block_no) xxmax = coord(1) + 0.5*size(1) xxmin = coord(1) - 0.5*size(1) yymax = coord(2) + 0.5*size(2) yymin = coord(2) - 0.5*size(2) zzmax = coord(3) + 0.5*size(3) zzmin = coord(3) - 0.5*size(3) ! Cell size delx = size(1) / real(nxb) dely = size(2) / real(nyb) delz = size(3) / real(nzb) !------------------------------------------------------------------------------- ! Write a message to stdout describing the problem setup. if (first_call .and. MyPE .eq. MasterPE) then write (*,*) call write_logfile & & ("flash: initializing for gas collapse problem.") write (*,*) & & 'flash: initializing for gas collapse problem.' write (*,*) 1 format (1X, 1P, 4(A13, E12.6, :, 1X)) 2 format (1X, 1P, A13, I12) endif ! Construct the radial samples needed for the initialization. if (first_call) then first_call = .false. ! rho2 = fltarr(n+1) ! rad2 = fltarr(n+1) ! mr2 = fltarr(n+1) deltar = sqrt(xctr*xctr*2)/N_prof const = 4*pi*Newton*(rho_0)**2/3.3e13/(gamma*(1)**2) mconst = 4*pi*rho_0*(xctr*.7)**3 nsubstep = 4 l = 0.0 dldr = 0.0 dmdr = 0 r = 0.0 m = 0.0 hdr = deltar/nsubstep*.5 ! rho_prof(0) = rho_0*exp(l) ! r_prof(0) = r !print*,'calculating table' ! do i=1,N_prof ! do j=1,nsubstep ! r = r + hdr ! mpred = m + hdr*dmdr ! lpred = l + hdr*dldr ! dmdr = r**2*exp(lpred) ! dldr = -const*mpred*exp(lpred*(1-gamma))/(r**2) ! r = r + hdr ! m = m + 2*hdr*dmdr ! l = l + 2*hdr*dldr ! dldr = -const*m*exp(l*(1-gamma))/(r**2) ! dmdr = r**2*exp(l) ! enddo ! rho_prof(i) = rho_0*exp(l) ! r_prof(i) = r ! enddo r= 1 r_prof(1) = 0 rho_prof(1) = 19.5 m = 4.0/3.0*pi*(1)*19.5 hdr=hdr*2 dldr = 0 l=rho_0 do i=1,N_prof do j=1,nsubstep if (r < xctr*.4) then dldr = -1/r**2!-m*newton*(l**(2-gamma))/(gamma*(3.3e13/(rho_0**gamma))*r**2)! *(l**(gamma-2))) l = 1e10/r else !isothermal solution if (iso_set) then B = -l*newton*m/(dldr*r**2) l = l*1e-4 ! l=1e-1 print*,'switched to isothermal, B = ',B iso_set = .false. endif dldr =-l/(B)*newton*m/r**2 l = hdr*dldr+l ! print*,'dldr = ',dldr endif m = m + 4./3.*pi*((r)**3-(r-hdr)**3)*l r = hdr+r !!print*,i,' ',j,' dldr ',dldr,' m ',m,' hdr ',hdr,' l ',l enddo !!print*,l,' ',r !!call abort !! print*,'dPdr ', gamma*(3.3e13/(rho_0**gamma))*l**(gamma-1)*dldr !! print*,'other side ',l*newton*m/r**2 rho_prof(i) = l r_prof(i) = r enddo endif !! !!------------------------------------------------------------------------------- ! 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 whether it is inside the ! initial radius or outside and initialize the hydro variables ! appropriately. Nint = 7 Nintinv = 1./float(Nint) Nintinv1= 1./(float(Nint)-1.) data(idens,:,:,:,block_no) = smlrho do k = 1, kmax do j = 1, jmax do i = 1, imax sum_rho = 1e-10 sum_p = 0 sum_vx = 0. sum_vy = 0. sum_vz = 0. do kk = 0, (Nint-1)*k3d zz = zzmin + delz*(-.5+real(k-nguard-1)+kk*Nintinv1) zdist = (zz - zctr) * k3d do jj = 0, (Nint-1)*k2d yy = yymin + dely*(-.5+real(j-nguard-1)+jj*Nintinv1) ydist = (yy - yctr) * k2d do ii = 0, Nint-1 xx = xxmin + delx*(-.5+real(i-nguard-1)+ii*Nintinv1) xdist = xx - xctr dist = sqrt( xdist**2 + ydist**2 + zdist**2 ) distinv = 1. / max( dist, 1.E-10 ) ! if (dist < (xctr*.8)) then ! sum_rho = 9947*distinv call find(r_prof,N_prof,dist,jlo) if(jlo>0) then sum_rho = sum_rho+rho_prof(jlo)*Nintinv*Nintinv*Nintinv else sum_rho = sum_rho+rho_prof(1)*4*Nintinv*Nintinv*Nintinv endif ! data(itemp,i,j,k,block_no) = 1.5e7 ! else ! data(itemp,i,j,k,block_no) = 1e-30 ! sum_rho = 1e-5 ! endif ! if(distinv == 1/1e-10) then ! sum_rho = 8*Nintinv*Nintinv*Nintinv ! endif enddo enddo enddo if (dist < (xctr*.4)) then ! data(itemp,i,j,k,block_no) = 1.7e9 ! data(ipres,i,j,k,block_no) = (data(idens,i,j,k,block_no)*gascon)*data(itemp,i,j,k,block_no) data(ipres,i,j,k,block_no) = sum_rho*(.01884)*(gamma-1) data(itemp,i,j,k,block_no) = data(ipres,i,j,k,block_no) / & & (data(idens,i,j,k,block_no)*gascon) else !isotherma9l9 call find(r_prof,N_prof,xctr*.3999,jlo) data(ipres,i,j,k,block_no) = rho_prof(jlo)*(.01884)*(gamma-1) !1e-10!B* sum_rho data(itemp,i,j,k,block_no) = data(ipres,i,j,k,block_no) / & & (data(idens,i,j,k,block_no)*gascon) endif data(itemp,i,j,k,block_no) = 1e-8 data(idens,i,j,k,block_no) = max(sum_rho, & & smlrho) ! data(ipres,i,j,k,block_no) = max(sum_p, & ! & smallp) data(ipres,i,j,k,block_no) = (data(idens,i,j,k,block_no)*gascon)*data(itemp,i,j,k,block_no) !************************** for constant pressure test ! data(ipres,i,j,k,block_no) = smallp * 100. ! data(itemp,i,j,k,block_no) = data(ipres,i,j,k,block_no) / & ! & (data(idens,i,j,k,block_no)*gascon) !************************** data(ivelx,i,j,k,block_no) = sum_vx data(ively,i,j,k,block_no) = sum_vy data(ivelz,i,j,k,block_no) = sum_vz !************************** for constant velocity gradient/uniform density test ! data(idens,i,j,k,block_no) = rho_0 ! data(ipres,i,j,k,block_no) = rho_0*presfac ! data(ivelx,i,j,k,block_no) = 3.4E9 * ! & (1.-(xxmin+delx*(i-nguard-0.5))/xmax) ! data(ively,i,j,k,block_no) = 0. ! data(ivelz,i,j,k,block_no) = 0. !************************** enddo enddo enddo !------------------------------------------------------------------------------- ! Initialize the nuclear abundances. These are not of interest ! in this problem, so we set them to 1. everywhere. do n = 1, ionmax do k = 1, kmax do j = 1, jmax do i = 1, imax data(inuc_begin+n-1,i,j,k,block_no) = 1. enddo enddo enddo enddo ! Compute the gas energy and set the gamma-values needed for ! the equation of state. do k = 1, kmax do j = 1, jmax do i = 1, imax data(igame,i,j,k,block_no) = gamma data(igamc,i,j,k,block_no) = gamma ek = 0.5 * (data(ivelx,i,j,k,block_no)**2 + & & data(ively,i,j,k,block_no)**2 + & & data(ivelz,i,j,k,block_no)**2) data(iener,i,j,k,block_no) = data(ipres,i,j,k,block_no) / & & (data(igame,i,j,k,block_no)-1.) data(iener,i,j,k,block_no) = data(iener,i,j,k,block_no) / & & data(idens,i,j,k,block_no) data(iener,i,j,k,block_no) = data(iener,i,j,k,block_no) + ek ! data(iener,i,j,k,block_no) = max(data(iener,i,j,k,block_no),smallp) enddo enddo enddo !=============================================================================== end !******************************************************************************* ! 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