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