!---------------------------------------------------------------------- ! This routine updates the following grid functions ! pp pi !---------------------------------------------------------------------- subroutine update0(pp_np1,pp_n,pi_np1,pi_n,a_np1,a_n,l_np1,l_n, & g1_Nr,r,dr,dt,epsdis) implicit none include 'globals.inc' integer g1_Nr real*8 pp_np1(1:g1_Nr) real*8 pp_n(1:g1_Nr) real*8 pi_np1(1:g1_Nr) real*8 pi_n(1:g1_Nr) real*8 a_np1(1:g1_Nr) real*8 a_n(1:g1_Nr) real*8 l_np1(1:g1_Nr) real*8 l_n(1:g1_Nr) real*8 r(*) real*8 dr real*8 dt real*8 epsdis integer i,j,k i=1 pp_np1(i)=pp_np1(i)-pp_np1(i) i=2 pp_np1(i)=pp_np1(i)-((pp_np1(i)-pp_n(i))/dt-((l_np1(i+1)*pi & _np1(i+1)/a_np1(i+1)-l_np1(i-1)*pi_np1(i-1)/a_np1(i-1))/(2* & dr)+(l_n(i+1)*pi_n(i+1)/a_n(i+1)-l_n(i-1)*pi_n(i-1)/a_n(i-1) & )/(2*dr))/2)/(dt/(dt*dt)) do i=3, g1_Nr-2, 1 pp_np1(i)=pp_np1(i)-((pp_np1(i)-pp_n(i))/dt-(((-epsdis)/(16*dt & )*(6*pp_np1(i)+pp_np1(i+2)+pp_np1(i-2)-4*(pp_np1(i+1)+pp_np1 & (i-1)))+(-epsdis)/(16*dt)*(6*pp_n(i)+pp_n(i+2)+pp_n(i-2)-4*( & pp_n(i+1)+pp_n(i-1))))/2+((l_np1(i+1)*pi_np1(i+1)/a_np1(i+1) & -l_np1(i-1)*pi_np1(i-1)/a_np1(i-1))/(2*dr)+(l_n(i+1)*pi_n( & i+1)/a_n(i+1)-l_n(i-1)*pi_n(i-1)/a_n(i-1))/(2*dr))/2))/(dt/( & dt*dt)-2*(-epsdis)/(16*dt)*(6D0)/(4D0)) end do i=g1_Nr-1 pp_np1(i)=pp_np1(i)-((pp_np1(i)-pp_n(i))/dt-((l_np1(i+1)*pi & _np1(i+1)/a_np1(i+1)-l_np1(i-1)*pi_np1(i-1)/a_np1(i-1))/(2* & dr)+(l_n(i+1)*pi_n(i+1)/a_n(i+1)-l_n(i-1)*pi_n(i-1)/a_n(i-1) & )/(2*dr))/2)/(dt/(dt*dt)) i=g1_Nr pp_np1(i)=pp_np1(i)-((pp_np1(i)-pp_n(i))/dt+((3*pp_np1(i)-4*pp & _np1(i-1)+pp_np1(i-2))/(2*dr)+pp_np1(i)/r(i)+(3*pp_n(i)-4*pp & _n(i-1)+pp_n(i-2))/(2*dr)+pp_n(i)/r(i))/2)/(dt/(dt*dt)+2*(2* & dr*(3D0)/(2*dr*2*dr)+r(i)/(r(i)*r(i)))/(4D0)) i=1 pi_np1(i)=pi_np1(i)-(+pi_np1(i)-4*pi_np1(i+1)/3+pi_np1(i+2)/3+ & +pi_n(i)-4*pi_n(i+1)/3+pi_n(i+2)/3)/2/(0.5D0) i=2 pi_np1(i)=pi_np1(i)-((pi_np1(i)-pi_n(i))/dt-(3*(r(i+1)**2*l & _np1(i+1)*pp_np1(i+1)/a_np1(i+1)-r(i-1)**2*l_np1(i-1)*pp_np1 & (i-1)/a_np1(i-1))/(r(i+1)**3-r(i-1)**3)+3*(r(i+1)**2*l_n(i+1 & )*pp_n(i+1)/a_n(i+1)-r(i-1)**2*l_n(i-1)*pp_n(i-1)/a_n(i-1))/ & (r(i+1)**3-r(i-1)**3))/2)/(dt/(dt*dt)) do i=3, g1_Nr-2, 1 pi_np1(i)=pi_np1(i)-((pi_np1(i)-pi_n(i))/dt-(((-epsdis)/(16*dt & )*(6*pi_np1(i)+pi_np1(i+2)+pi_np1(i-2)-4*(pi_np1(i+1)+pi_np1 & (i-1)))+(-epsdis)/(16*dt)*(6*pi_n(i)+pi_n(i+2)+pi_n(i-2)-4*( & pi_n(i+1)+pi_n(i-1))))/2+(3*(r(i+1)**2*l_np1(i+1)*pp_np1(i+1 & )/a_np1(i+1)-r(i-1)**2*l_np1(i-1)*pp_np1(i-1)/a_np1(i-1))/( & r(i+1)**3-r(i-1)**3)+3*(r(i+1)**2*l_n(i+1)*pp_n(i+1)/a_n(i+1 & )-r(i-1)**2*l_n(i-1)*pp_n(i-1)/a_n(i-1))/(r(i+1)**3-r(i-1)** & 3))/2))/(dt/(dt*dt)-2*(-epsdis)/(16*dt)*(6D0)/(4D0)) end do i=g1_Nr-1 pi_np1(i)=pi_np1(i)-((pi_np1(i)-pi_n(i))/dt-(3*(r(i+1)**2*l & _np1(i+1)*pp_np1(i+1)/a_np1(i+1)-r(i-1)**2*l_np1(i-1)*pp_np1 & (i-1)/a_np1(i-1))/(r(i+1)**3-r(i-1)**3)+3*(r(i+1)**2*l_n(i+1 & )*pp_n(i+1)/a_n(i+1)-r(i-1)**2*l_n(i-1)*pp_n(i-1)/a_n(i-1))/ & (r(i+1)**3-r(i-1)**3))/2)/(dt/(dt*dt)) i=g1_Nr pi_np1(i)=pi_np1(i)-((pi_np1(i)-pi_n(i))/dt+((3*pi_np1(i)-4*pi & _np1(i-1)+pi_np1(i-2))/(2*dr)+pi_np1(i)/r(i)+(3*pi_n(i)-4*pi & _n(i-1)+pi_n(i-2))/(2*dr)+pi_n(i)/r(i))/2)/(dt/(dt*dt)+2*(2* & dr*(3D0)/(2*dr*2*dr)+r(i)/(r(i)*r(i)))/(4D0)) return end !---------------------------------------------------------------------- ! This routine updates the following grid functions ! a l tmr !---------------------------------------------------------------------- subroutine calc_al(pp,ppnm1,pi,pinm1,a,anm1,l,lnm1,tmr,f0,f2,g0, & g1_Nr,r,t,flatspace,hcdlnatol,hcmxiter,hcrc,tmrbh) implicit none include 'globals.inc' integer g1_Nr real*8 pp(1:g1_Nr) real*8 ppnm1(1:g1_Nr) real*8 pi(1:g1_Nr) real*8 pinm1(1:g1_Nr) real*8 a(1:g1_Nr) real*8 anm1(1:g1_Nr) real*8 l(1:g1_Nr) real*8 lnm1(1:g1_Nr) real*8 tmr(1:g1_Nr) real*8 f0(1:g1_Nr) real*8 f2(1:g1_Nr) real*8 g0(1:g1_Nr) real*8 r(*) real*8 t integer flatspace real*8 hcdlnatol integer hcmxiter integer hcrc real*8 tmrbh c----------------------------------------------------------------------- c Included code for Hamiltonian constraint/slicing c equation solvers and computation of 2m/r. c c The RNPL generated f77 code (in updates.f) which 'include's c update segments such as the following also 'include's c the file 'global.inc' that defines some global variables c maintained by a canonical RNPL f77 code. c Be careful that variables you use in hand-coded c update segments do not collide with those defined in c 'global.inc' (examples: 'iter', 'done', 'nr0') c----------------------------------------------------------------------- integer getu real*8 cpi parameter ( cpi = 3.141 5926 5358 9793 d0 ) real*8 pph, pih, ah, wh, rh real*8 lnlj, lnljp1 real*8 lnaj, lnajp1, dlnaj, dlnajp1, & drm1, resid, jacel integer hciter integer nr, j, ubh logical hcdone logical vstrace parameter ( vstrace = .false. ) logical hctrace parameter ( hctrace = .false. ) logical ltrace parameter ( ltrace = .false. ) logical firstflatspace data firstflatspace / .true. / save c----------------------------------------------------------------------- c Use a more convenient variable for the number of spatial c grid points. c----------------------------------------------------------------------- nr = g1_nr if( vstrace ) then call vsxynt('pp',t,r,pp,nr) call vsxynt('pi',t,r,pi,nr) call vsxynt('ppnm1',t,r,ppnm1,nr) call vsxynt('pinm1',t,r,pinm1,nr) end if if( ltrace ) then write(0,*) 'calc_al: In routine with g1_nr = ', g1_nr write(0,*) 'calc_al: t = ', t write(0,*) 'calc_al: hcdlnatol = ', hcdlnatol write(0,*) 'calc_al: hcmxiter = ', hcmxiter write(0,*) 'calc_al: flatspace = ', flatspace write(0,*) end if c----------------------------------------------------------------------- c If flatspace option enabled, set vars. to flat spacetime c values and return. c----------------------------------------------------------------------- if( flatspace .ne. 0 ) then if( firstflatspace ) then write(0,*) 'calc_al: *** Flatspace option enabled' firstflatspace = .false. end if do j = 1 , nr a(j) = 1.0d0 l(j) = 1.0d0 tmr(j) = 0.0d0 end do return end if c----------------------------------------------------------------------- c Compute coeficient functions for Hamiltonian constraint c----------------------------------------------------------------------- do j = 1 , nr - 1 pph = 0.5d0 * (pp(j+1) + pp(j)) pih = 0.5d0 * (pi(j+1) + pi(j)) rh = 0.5d0 * (r(j+1) + r(j)) f0(j) = -0.5d0 / rh - 2.0d0 * cpi * rh * (pph**2 + pih**2) f2(j) = 0.5d0 / rh end do if( vstrace ) then call vsxynt('f0',t,r,f0,nr-1) call vsxynt('f2',t,r,f0,nr-1) end if c----------------------------------------------------------------------- c Solve Hamiltonian constraint for A := ln(a) outwards from r=0, c with boundary condition A(0) = 0, using pointwise c Newton iteration. Simultaneously compute a = exp(A) and c 2m/r = 1 - a^{-2}. c----------------------------------------------------------------------- a(1) = 1.0d0 tmr(1) = 0.0d0 lnaj = 0.0d0 do j = 1 , nr - 1 if( j .eq. 1 ) then lnajp1 = lnaj else lnajp1 = lnaj + dlnaj end if hcdone = .false. hciter = 0 drm1 = 1.0d0 / (r(j+1) - r(j)) c----------------------------------------------------------------------- c Note that the 'do while()' construct, although not f77 c *is* f90, and supported by our f77 compilers. I consider c is a perfectly safe extension to f77. c----------------------------------------------------------------------- do while( .not. hcdone ) resid = drm1 * (lnajp1 - lnaj) + & f2(j) * exp(lnajp1 + lnaj) + & f0(j) jacel = drm1 + f2(j) * exp(lnajp1 + lnaj) dlnajp1 = resid / jacel hciter = hciter + 1 lnajp1 = lnajp1 - dlnajp1 if( hctrace ) then write(0,1000) j, hciter, resid, dlnajp1 1000 format(' j-1',i5,': Step',i3,' res: ', & 1pe10.2,' d(ln a): ',1pE10.2) end if if( abs(dlnajp1) .le. hcdlnatol ) then hcdone = .true. else if( hciter .gt. hcmxiter ) then c----------------------------------------------------------------------- c Exit program if Newton iteration didn't converge. c----------------------------------------------------------------------- write(0,*) 'calc_al: Newton iteration did not ', & 'converge at j = ', j, ' r = ', r(j) stop end if end do dlnaj = lnajp1 - lnaj a(j+1) = exp(lnajp1) tmr(j+1) = 1.0d0 - 1.0d0 / a(j+1)**2 c----------------------------------------------------------------------- c Check for black hole formation. c----------------------------------------------------------------------- if( tmr(j+1) .ge. tmrbh ) then write(0,2000) t, r(j+1), tmr(j+1) 2000 format(' calc_al: Black hole detected at t = ',f7.2, & ' r = ',f8.3,' 2m/r = ',f7.3) c----------------------------------------------------------------------- c When collapse is detected, create a file called BH then c exit. Provides simple mechanism to communicate black hole c formation to the "outside world" (e.g. to a shell c script. c----------------------------------------------------------------------- ubh = getu() open(ubh,file='BH',form='formatted') write(ubh,*) 1 close(ubh) stop end if lnaj = lnajp1 if( hctrace ) then write(0,*) end if end do c----------------------------------------------------------------------- c Compute coeficient functions for polar slicing equation. c----------------------------------------------------------------------- do j = 1 , nr - 1 pph = 0.5d0 * (pp(j+1) + pp(j)) pih = 0.5d0 * (pi(j+1) + pi(j)) ah = 0.5d0 * (a(j+1) + a(j)) rh = 0.5d0 * (r(j+1) + r(j)) g0(j) = -0.5d0 * (ah**2 - 1.0d0) / rh - & 2.0d0 * cpi * rh * (pph**2 + pih**2) end do c----------------------------------------------------------------------- c Solve slicing equation for lapse inwards from r=rmax, with c boundary condition l(rmax) = 1 / a(rmax). c----------------------------------------------------------------------- l(nr) = 1.0d0 / a(nr) lnljp1 = log(l(nr)) do j = nr - 1 , 1 , -1 lnlj = lnljp1 + (r(j+1) - r(j)) * g0(j) l(j) = exp(lnlj) lnljp1 = lnlj end do return end