c/*@@ c @file AHFinder.F c @date April 1998 c @author Miguel Alcubierre c @desc c Find an apparent horizon. c @enddesc c@@*/ #include "cctk.h" #include "cctk_parameters.h" #include "cctk_arguments.h" subroutine AHFinder(CCTK_FARGUMENTS) use AHFinder_dat implicit none DECLARE_CCTK_FARGUMENTS DECLARE_CCTK_PARAMETERS logical status,horizon,gauss integer handle,x_index,y_index,z_index,ierror integer ahfso,ahfafter integer i,j,l,m integer NN integer mtype integer getoutpfx integer CCTK_Equals integer CCTK_MyProc integer CCTK_nProcs integer istat CCTK_REAL atime CCTK_REAL zero,half,one CCTK_REAL xmn,ymn,zmn,xmx,ymx,zmx,rmx,r0 CCTK_REAL intexp_h,intexp2_h,intarea_h CCTK_REAL circ_eq,meri_p1,meri_p2 CCTK_REAL aux,expin,expout,small1,small2 CCTK_REAL c0_temp(18) CCTK_REAL maxval(3),minval(3) character*200 almf,logf,areaf,massf,circeqf,merip1f,merip2f save ahfafter,ahfso,atime #include "CactusEinstein/Einstein/src/Einstein.h" ! Description of variables: ! ! status Surface found? ! ! horizon Horizon found? ! ! gauss Ouput Gaussian curvature? ! ! ahfso How often we want to find a horizon? ! ahfafter After how many iterations start looking for horizons? ! ! i,l,m Counters. ! ! NN Total number of terms in expansion. ! ! xmn Minimum value of x over the grid. ! xmx Maximum value of x over the grid. ! ! ymn Minimum value of y over the grid. ! ymx Maximum value of y over the grid. ! ! zmn Minimum value of z over the grid. ! zmx Maximum value of z over the grid. ! ! rmx Radius of largest sphere that fits in the grid. ! ! circ_eq Equatorial circumference. ! meri_p1 Length of meridian at phi=0. ! meri_p2 Length of meridian at phi=pi/2. ! ! mtype Type of surface. ! 1. Integral of H^2 is positive out and negative in. ! 2. Integral of H^2 is negative out and positive in. ! 3. Integral of H^2 is negative out and in. ! 4. Integral of H^2 is positive out and in. ! ************************** ! *** DEFINE NUMBERS *** ! ************************** zero = 0.0D0 half = 0.5D0 one = 1.0D0 dx = cctk_delta_space(1) dy = cctk_delta_space(2) dz = cctk_delta_space(3) ! ******************************** ! *** INITIALIZE {status} *** ! ******************************** status = .false. myproc = CCTK_MyProc(cctkGH) nprocs = CCTK_nProcs(cctkGH) ! ********************************* ! *** INITIALIZE FIRSTCALLS *** ! ********************************* firstleg = .true. firstint = .true. ! *********************************************** ! *** CHECK IF WE WANT TO FIND 3 HORIZONS *** ! *********************************************** if (ahf_find3.ne.0) then find3 = .true. else find3 = .false. end if ! ********************************************* ! *** CHECK IF WE WANT TO FIND A HORIZON *** ! ********************************************* ! See if we want to find a horizon at this time. if (ahf_ncall.eq.0) then atime = ahf_findaftertime ahfafter = ahf_findafter ahfso = ahf_findevery if (atime.le.zero) then if (cctk_iteration.lt.ahfafter) return else if (cctk_time.lt.atime) return ahfafter = cctk_iteration end if end if if (mod(cctk_iteration-ahfafter,ahfso).ne.0) return ! Sorry for this hack, but the call order for BoxInBox ! is screwy in cactus 3.2 #ifdef THORN_BOXINBOX if ((boxinbox.ne.0).and.(ahf_ncall.eq.0)) then return endif #endif ! ******************************************* ! *** ADD TO COUNT OF CALLS TO FINDER *** ! ******************************************* 1 continue ! Say what we are doing. if (find3) then if (mfind.eq.0) then write(*,*) write(*,*) 'AHFinder: Looking for 3 apparent horizons' write(*,*) 'AHFinder: Searching for first horizon' else if (mfind.eq.1) then write(*,*) write(*,*) 'AHFinder: Searching for second horizon' else if (mfind.eq.2) then write(*,*) write(*,*) 'AHFinder: Searching for third horizon' end if else write(*,*) write(*,*) 'AHFinder: Looking for horizon' end if ! Add to count of calls to finder. Notice that this counter ! was initialized to zero, so it must be increased now to 1. ! If we look for 3 horizons, it should be done before looking ! for the first one. if (find3) then if (mfind.eq.0) then ahf_ncall = ahf_ncall + 1 end if else ahf_ncall = ahf_ncall + 1 end if ! ******************************************************* ! *** FIND OUT IF WE WANT TO FIND TRAPPED SURFACE *** ! ******************************************************* if (ahf_trapped_surface.ne.0) then find_trapped_surface = .true. else find_trapped_surface = .false. end if ! ******************************************** ! *** GET THE NAME OF OUTPUT DIRECTORY *** ! ******************************************** call CCTK_FortranString(nfile,outdir,filestr) if (find3) then if (mfind.eq.0) then logf = filestr(1:nfile)//"/ahf_logfile_0" almf = filestr(1:nfile)//"/ahf_coeff_0.alm" areaf = filestr(1:nfile)//"/ahf_area_0.tl" massf = filestr(1:nfile)//"/ahf_mass_0.tl" circeqf = filestr(1:nfile)//"/ahf_circ_eq_0.tl" merip1f = filestr(1:nfile)//"/ahf_meri_p1_0.tl" merip2f = filestr(1:nfile)//"/ahf_meri_p2_0.tl" else if (mfind.eq.1) then logf = filestr(1:nfile)//"/ahf_logfile_1" almf = filestr(1:nfile)//"/ahf_coeff_1.alm" areaf = filestr(1:nfile)//"/ahf_area_1.tl" massf = filestr(1:nfile)//"/ahf_mass_1.tl" circeqf = filestr(1:nfile)//"/ahf_circ_eq_1.tl" merip1f = filestr(1:nfile)//"/ahf_meri_p1_1.tl" merip2f = filestr(1:nfile)//"/ahf_meri_p2_1.tl" else logf = filestr(1:nfile)//"/ahf_logfile_2" almf = filestr(1:nfile)//"/ahf_coeff_2.alm" areaf = filestr(1:nfile)//"/ahf_area_2.tl" massf = filestr(1:nfile)//"/ahf_mass_2.tl" circeqf = filestr(1:nfile)//"/ahf_circ_eq_2.tl" merip1f = filestr(1:nfile)//"/ahf_meri_p1_2.tl" merip2f = filestr(1:nfile)//"/ahf_meri_p2_2.tl" end if else logf = filestr(1:nfile)//"/ahf_logfile" almf = filestr(1:nfile)//"/ahf_coeff.alm" areaf = filestr(1:nfile)//"/ahf_area.tl" massf = filestr(1:nfile)//"/ahf_mass.tl" circeqf = filestr(1:nfile)//"/ahf_circ_eq.tl" merip1f = filestr(1:nfile)//"/ahf_meri_p1.tl" merip2f = filestr(1:nfile)//"/ahf_meri_p2.tl" end if ! ************************************************** ! *** TRANSFORM CONFORMAL TO PHYSICAL METRIC *** ! ************************************************** call CCTK_WARN(4,"Converting metric: conformal -> physical") gxx=gxx*psi**4 gxy=gxy*psi**4 gxz=gxz*psi**4 gyy=gyy*psi**4 gyz=gyz*psi**4 gzz=gzz*psi**4 conformal_state = NOCONFORMAL_METRIC ! *************************************** ! *** FIND OUT THE TYPE OF OUTPUT *** ! *************************************** if (ahf_logfile.ne.0) then logfile = .true. else logfile = .false. end if if (ahf_verbose.ne.0) then verbose = .true. else verbose = .false. end if if (ahf_veryverbose.ne.0) then veryver = .true. else veryver = .false. end if if (ahf_guessverbose.ne.0) then guessverbose = .true. else guessverbose = .false. end if if (logfile.and.(myproc.eq.0)) then if (logfile) then if (myproc.eq.0) then open(1,file=logf,form='formatted',status='replace') write(1,*) write(1,*) 'LOG FILE FOR APPARENT HORIZON FINDER' write(1,*) write(1,*) 'Note: During an evolution, this file will only' write(1,*) 'contain information about the last time that the' write(1,*) 'finder was called.' write(1,*) close(1) end if end if end if ! ******************************************* ! *** FIND {xmn,xmx,ymn,ymx,zmn,zmx} *** ! ******************************************* call CCTK_CoordRange(ierr,cctkGH,xmn,xmx,"x") call CCTK_CoordRange(ierr,cctkGH,ymn,ymx,"y") call CCTK_CoordRange(ierr,cctkGH,zmn,zmx,"z") ! ********************************** ! *** INITIALIZE {xc,yc,zc} *** ! ********************************** if (find3) then if (mfind.eq.0) then xc = ahf_xc_0 yc = ahf_yc_0 zc = ahf_zc_0 else if (mfind.eq.1) then xc = ahf_xc_1 yc = ahf_yc_1 zc = ahf_zc_1 else xc = ahf_xc_2 yc = ahf_yc_2 zc = ahf_zc_2 end if else xc = ahf_xc yc = ahf_yc zc = ahf_zc end if if ((xc.gt.xmx).or.(xc.lt.xmn)) then write(*,*) write(*,*) 'xc is out of the grid, reseting it to zero ...' write(*,*) xc = zero end if if ((yc.gt.ymx).or.(yc.lt.ymn)) then write(*,*) write(*,*) 'yc is out of the grid, reseting it to zero ...' write(*,*) yc = zero end if if ((yc.gt.ymx).or.(yc.lt.ymn)) then write(*,*) write(*,*) 'zc is out of the grid, reseting it to zero ...' write(*,*) zc = zero end if ! ************************** ! *** FIND SYMMETRIES *** ! ************************** ! Find out if we want to move the center. if (ahf_offset.ne.0) then offset = .true. else offset = .false. end if if (ahf_wander.ne.0) then wander = .true. else wander = .false. end if if ((xc.ne.zero).or.(yc.ne.zero).or.(zc.ne.zero).or. . (wander)) offset = .true. ! Find out if we are in non-axisymmetric mode. if (ahf_phi.ne.0) then nonaxi = .true. else nonaxi = .false. end if ! Find out reflection symmetries. if ((CCTK_Equals(ahf_octant,"yes").eq.1).or. . (CCTK_Equals(ahf_octant,"high").eq.1)) then refx = .true. refy = .true. refz = .true. else if (ahf_refx.ne.0) then refx = .true. else refx = .false. end if if (ahf_refy.ne.0) then refy = .true. else refy = .false. end if if (ahf_refz.ne.0) then refz = .true. else refz = .false. end if end if if (ahf_phi.eq.0) then refx = .true. refy = .true. end if ! Force grid symmetries. if (CCTK_Equals(domain,"octant").eq.1) then if (xc.eq.zero) refx = .true. if (yc.eq.zero) refy = .true. if (zc.eq.zero) refz = .true. else if (CCTK_Equals(domain,"quadrant").eq.1) then if (xc.eq.zero) refx = .true. if (yc.eq.zero) refy = .true. end if ! Find {stepx,stepy,stepz}. stepx = 0 stepy = 0 stepz = 0 if (refx) stepx=1 if (refy) stepy=1 if (refz) stepz=1 if (CCTK_Equals(ahf_octant,"high").eq.1) stepx=3 ! Check if we are using Cartoon_2d, in which case ! we override the value of "ahf_phi" since we have ! axisymmetry. cartoon = .false. #ifdef THORN_cartoon_2d if (cartoon_active.ne.0) then cartoon = .true. else cartoon = .false. end if #endif if (cartoon) then nonaxi = .false. refx = .true. refy = .true. end if ! ************************ ! *** FIND {lmax} *** ! ************************ lmax = ahf_lmax if (lmax.ge.20) then write(*,*) write(*,*) 'ahf_lmax must be smaller than 20.' write(*,*) 'STOPPED IN AHFinder.F' write(*,*) STOP end if ! ******************************* ! *** FIND {ntheta,nphi} *** ! ******************************* ntheta = ahf_ntheta nphi = ahf_nphi ! For cartoon, I use only 2 points in the phi direction. ! I could use one, but using 2 minimizes the changes in ! the integration routine. if (cartoon) then nphi = 1 end if ! *************************************************** ! *** ALLOCATE STORAGE FOR COEFFICIENT ARRAYS *** ! *************************************************** if ((ahf_ncall.eq.1).and.(mfind.eq.0)) then allocate(c0(0:lmax),c0_old(0:lmax)) allocate(cc(lmax,lmax),cc_old(lmax,lmax)) allocate(cs(lmax,lmax),cs_old(lmax,lmax)) cs = zero cc = zero cc_old = zero cs_old = zero end if ! ******************************************************** ! *** FIND TOTAL NUMBER OF TERMS IN EXPANSION {NN} *** ! ******************************************************** NN = 0 ! Count {xc,yc,zc}. if (wander) then if (.not.refx) NN = NN + 1 if (.not.refy) NN = NN + 1 if (.not.refz) NN = NN + 1 end if ! Count c0(l). do l=0,lmax,1+stepz NN = NN + 1 end do if (nonaxi) then ! Count cc(l,m). do l=1,lmax do m=1+stepx,l,1+stepx if (stepz*mod(l-m,2).eq.0) NN = NN +1 end do end do ! Count cs(l,m). if (.not.refy) then do l=1,lmax do m=1,l,1+stepx if (stepz*mod(l-m,2).eq.0) NN = NN + 1 end do end do end if end if ! ***************************************** ! *** FIND RADIUS OF INITIAL SPHERE *** ! ***************************************** ! Find largest sphere that fits into the grid, taking ! into account the position of the center. Then multiply ! this number with 0.8 to leave a safety margin at the edge ! of the grid. ! Cartoon. if (cartoon) then xc = zero yc = zero rmx = min(xmx,zmx-zc,zc-zmn) ! Octant. else if (CCTK_Equals(domain,"octant").eq.1) then if ((xc.eq.zero).and.(yc.eq.zero).and.(zc.eq.zero)) then rmx = min(xmx,ymx,zmx) else if ((xc.eq.zero).and.(yc.eq.zero)) then rmx = min(xmx,ymx,zmx-zc,zc) else if ((xc.eq.zero).and.(zc.eq.zero)) then rmx = min(xmx,ymx-yc,zmx,yc) else if ((yc.eq.zero).and.(zc.eq.zero)) then rmx = min(xmx-xc,ymx,zmx,xc) else if (xc.eq.zero) then rmx = min(xmx,ymx-yc,zmx-zc,yc,zc) else if (yc.eq.zero) then rmx = min(xmx-xc,ymx,zmx-zc,xc,zc) else if (zc.eq.zero) then rmx = min(xmx-xc,ymx-yc,zmx,xc,yc) else rmx = min(xmx-xc,ymx-yc,zmx-zc,xc,yc,zc) end if ! Quadrant. else if (CCTK_Equals(domain,"quadrant").eq.1) then if ((xc.eq.zero).and.(yc.eq.zero)) then rmx = min(xmx,ymx,zmx-zc,zc-zmn) else if (xc.eq.zero) then rmx = min(xmx,ymx-yc,zmx-zc,yc,zc-zmn) else if (yc.eq.zero) then rmx = min(xmx-xc,ymx,zmx-zc,xc,zc-zmn) else rmx = min(xmx-xc,ymx-yc,zmx-zc,xc,yc,zc-zmn) end if ! Full. else rmx = min(xmx-xc,ymx-yc,zmx-zc,xc-xmn,yc-ymn,zc-zmn) end if rmx = 0.8*rmx ! Check if a special radius for the initial sphere ! has been forced by the parameter file. if (find3) then if (mfind.eq.0) then r0 = ahf_r0_0 else if (mfind.eq.1) then r0 = ahf_r0_1 else r0 = ahf_r0_2 end if else r0 = ahf_r0 end if if (r0.ne.zero) then if (r0.le.rmx) then rmx = r0 else write (*,*) write (*,*) 'ahf_r0 is too large, rescaling ...' write (*,*) end if end if ! ************************************************* ! *** GET INITIAL GUESS AND LOOK FOR HORIZON *** ! ************************************************* ! initializing flags sloppy = (ahf_sloppyguess.eq.1) if (cartoon) sloppy = .true. flow = (ahf_flow.eq.1) inner = (ahf_inner.eq.1) if (ahf_guessold.ne.0) then guessold = .true. else guessold = .false. end if if (ahf_guess_absmin.ne.0) then guess_absmin = .true. else guess_absmin = .false. end if if (ahf_manual_guess.ne.0) then manual_guess = .true. else manual_guess = .false. end if if (ahf_minarea.ne.0) then minarea = .true. else minarea = .false. end if ! Get initial guess and jump into finder routines if ((.not.flow).and.guessold.and.(ahf_ncall.gt.1)) then ! Minimization algorithm and old guess. call AHFinder_initguess(CCTK_FARGUMENTS,rmx) call AHFinder_min(CCTK_FARGUMENTS,NN,status,logf) else if ((.not.flow).and.manual_guess) then ! Minimization algorithm and manual guess. if (CCTK_Equals(domain,"octant").eq.0) then write(*,*) 'AHFinder: if you want to use ahf_manual_guess' write(*,*) 'without grid=octant, you will have to make' write(*,*) 'it work!' STOP endif c0(0) = ahf_l0_guess c0_temp(2) = ahf_l2_guess c0_temp(4) = ahf_l4_guess c0_temp(6) = ahf_l6_guess c0_temp(8) = ahf_l8_guess c0_temp(10) = ahf_l10_guess c0_temp(12) = ahf_l12_guess c0_temp(14) = ahf_l14_guess c0_temp(16) = ahf_l16_guess c0_temp(18) = ahf_l18_guess do i=2,lmax,2 c0(i) = c0_temp(i) enddo call AHFinder_initguess(CCTK_FARGUMENTS,rmx) call AHFinder_min(CCTK_FARGUMENTS,NN,status,logf) else if (.not.flow) then ! Straight minimization algorithm. call AHFinder_initguess(CCTK_FARGUMENTS,rmx) call AHFinder_min(CCTK_FARGUMENTS,NN,status,logf) else ! Flow algorithm. if ((ahf_ncall.eq.1).and.(mfind.eq.0)) then allocate(hflow0(0:lmax),cflow0(0:lmax),nflow0(0:lmax)) allocate(hflowc(lmax,lmax),cflowc(lmax,lmax),nflowc(lmax,lmax)) allocate(hflows(lmax,lmax),cflows(lmax,lmax),nflows(lmax,lmax)) end if call AHFinder_flow(CCTK_FARGUMENTS,rmx,status,logf) flow = .false. end if ! *************************************** ! *** WRITE MESSAGES AND LOG FILE *** ! *************************************** if (status) then ! Check if the expansion on the surface is small enough ! for a horizon. How small is small enough is of course ! pretty arbitrary, so to reduce the subjectivity I use ! three separate criteria: ! ! 1) The integral of the square of the expansion at the surface ! must be smaller than a given threshold (small1). ! ! 2) The mean value of the expansion on the surface must be ! smaller in absolute value than its standard deviation. ! This means that the mean value of the expansion is ! consistent with zero. ! ! 3) The integral of the square of the expansion at the surface ! must be smaller by a factor of "small2" than the same ! integral slightly outside and slightly inside. horizon = .true. small1 = 1.0D-2 small2 = 0.5D0 ! Surface found. call AHFinder_fun(CCTK_FARGUMENTS) call AHFinder_exp(CCTK_FARGUMENTS) call AHFinder_int(CCTK_FARGUMENTS) if (ahf_gaussout.ne.0) then call AHFinder_gau(CCTK_FARGUMENTS,circ_eq,meri_p1,meri_p2) end if intexp_h = intexp intexp2_h = intexp2 intarea_h = intarea if (intexp2_h.gt.small1) horizon = .false. ! Find standard deviation: sqrt( - ^2) aux = dsqrt(intexp2_h/intarea_h - (intexp_h/intarea_h)**2) if (dabs(intexp_h/intarea_h)-aux.gt.zero) horizon = .false. ! Surface slightly inside aux = c0(0) c0(0) = aux - half*min(dx,dy,dz) call AHFinder_int(CCTK_FARGUMENTS) expin = intexp if (intexp2_h.gt.small2*intexp2) horizon = .false. if (veryver) then write(*,*) write(*,*) 'Surface slightly inside the minimum surface: ' write(*,"(A8,ES14.6)") ' H : ',intexp write(*,"(A8,ES14.6)") ' H^2: ',intexp2 write(*,"(A13,ES14.6)") ' H/area : ',intexp/intarea write(*,"(A13,ES14.6)") ' H^2/area: ',intexp2/intarea write(*,*) 'number of interpolated points: ', . inside_min_count write(*,*) ' number of those that are negative: ', . inside_min_neg_count endif ! Surface slightly outside. c0(0) = aux + half*min(dx,dy,dz) call AHFinder_int(CCTK_FARGUMENTS) expout = intexp if (intexp2_h.gt.small2*intexp2) horizon = .false. if (veryver) then write(*,*) write(*,*) 'Surface slightly outside the minimum surface: ' write(*,"(A8,ES14.6)") ' H : ',intexp write(*,"(A8,ES14.6)") ' H^2: ',intexp2 write(*,"(A13,ES14.6)") ' H/area : ',intexp/intarea write(*,"(A13,ES14.6)") ' H^2/area: ',intexp2/intarea endif c0(0) = aux ! What kind of horizon? if ((expin.le.zero).and.(expout.ge.zero)) then mtype = 1 else if ((expin.ge.zero).and.(expout.le.zero)) then mtype = 2 else if ((expin.le.zero).and.(expout.le.zero)) then mtype = 3 else mtype = 4 end if ! Messages to screen. if (myproc.eq.0) then ! Open logfile. if (logfile) then open(1,file=logf,form='formatted', . status='old',position='append') end if ! Surface information. Notice that here I define ! the surface "mass" as: ! ! M = sqrt( A / (16 pi) ) = 0.141047396 sqrt(A) if (intarea_h.ne.zero) then aux = one/intarea_h else aux = one end if write(*,*) write(*,*) 'AHFinder: Surface found.' write(*,*) write(*,"(A21,ES14.6)") ' Surface area =',intarea_h write(*,"(A21,ES14.6)") ' Surface mass =', . 0.141047396D0*dsqrt(intarea_h) write(*,"(A21,ES14.6)") ' Mean value of H =',aux*intexp_h write(*,"(A21,ES14.6)") ' Mean value of H^2 =',aux*intexp2_h write(*,"(A21,ES14.6)") ' Standard deviation =', . dsqrt(aux*intexp2_h - (aux*intexp_h)**2) write(*,*) write(*,"(A11,ES14.6)") ' circ_eq =',circ_eq write(*,"(A11,ES14.6)") ' meri_p1 =',meri_p1 write(*,"(A11,ES14.6)") ' meri_p2 =',meri_p2 write(*,*) if (.false.) then if (logfile) then write(1,*) write(1,*) 'AHFinder: Surface found.' write(1,*) write(1,"(A21,ES14.6)") ' Surface area =',intarea_h write(1,"(A21,ES14.6)") ' Surface mass =', . 0.141047396D0*dsqrt(intarea_h) write(1,"(A21,ES14.6)") ' Mean value of H =',aux*intexp_h write(1,"(A21,ES14.6)") ' Mean value of H^2 =',aux*intexp2_h write(1,"(A21,ES14.6)") ' Standard deviation =', . dsqrt(aux*intexp2_h - (aux*intexp_h)**2) write(1,*) write(1,"(A11,ES14.6)") ' circ_eq =',circ_eq write(1,"(A11,ES14.6)") ' meri_p1 =',meri_p1 write(1,"(A11,ES14.6)") ' meri_p2 =',meri_p2 write(1,*) end if end if if (mtype.eq.1) then if (horizon) then write(*,*) 'The surface seems to be an outer horizon.' if (logfile) then write(1,*) 'The surface seems to be an outer horizon.' end if else write(*,*) 'The surface is probably an outer horizon, but' write(*,*) 'the mean of the square of the expansion is' write(*,*) 'not small enough. Try a larger lmax, or a' write(*,*) 'smaller grid size.' if (logfile) then write(1,*) 'The surface is probably an outer horizon, but' write(1,*) 'the mean of the square of the expansion is' write(1,*) 'not small enough. Try a larger lmax, or a' write(1,*) 'smaller grid size.' end if end if else if (mtype.eq.2) then if (horizon) then write(*,*) 'The surface seems to be an inner horizon.' if (logfile) then write(1,*) 'The surface seems to be an inner horizon.' end if else write(*,*) 'The surface is probably an inner horizon, but' write(*,*) 'the mean of the square of the expansion is' write(*,*) 'not small enough. Try a larger lmax, or a' write(*,*) 'smaller grid size.' if (logfile) then write(1,*) 'The surface is probably an inner horizon, but' write(1,*) 'the mean of the square of the expansion is' write(1,*) 'not small enough. Try a larger lmax, or a' write(1,*) 'smaller grid size.' end if end if else if (mtype.eq.3)then if (horizon) then write(*,*) 'The surface seems to be a marginal horizon' write(*,*) '(the integral of the expansion is negative' write(*,*) 'on both sides.)' if (logfile) then write(1,*) 'The surface seems to be a marginal horizon' write(1,*) '(the integral of the expansion is negative' write(1,*) 'on both sides.)' end if else write(*,*) 'The surface seems to be a trapped surface.' if (logfile) then write(1,*) 'The surface seems to be a trapped surface.' end if end if else if (mtype.eq.4)then if (horizon) then write(*,*) 'The surface seems to be a marginal horizon' write(*,*) '(the integral of the expansion is positive' write(*,*) 'on both sides.)' if (logfile) then write(1,*) 'The surface seems to be a marginal horizon' write(1,*) '(the integral of the expansion is positive' write(1,*) 'on both sides.)' end if else write(*,*) 'The surface does not seem to be a horizon.' if (logfile) then write(1,*) 'The surface does not seem to be a horizon.' end if end if end if ! Shape. if (verbose) then write(*,*) write(*,*) 'Shape coefficients:' if (offset) then write(*,*) write(*,"(A6,ES14.6)") ' xc =',xc write(*,"(A6,ES14.6)") ' yc =',yc write(*,"(A6,ES14.6)") ' zc =',zc end if write(*,*) write(*,"(A4,I2,A3,ES14.6)") ' c0(',0,') =',c0(0) do l=1+stepz,lmax,1+stepz write(*,"(A4,I2,A3,ES14.6)") ' c0(',l,') =',c0(l) end do end if if (logfile) then write(1,*) write(1,*) 'Shape coefficients:' if (offset) then write(1,*) write(1,"(A6,ES14.6)") ' xc =',xc write(1,"(A6,ES14.6)") ' yc =',yc write(1,"(A6,ES14.6)") ' zc =',zc end if write(1,*) write(1,"(A4,I2,A3,ES14.6)") ' c0(',0,') =',c0(0) do l=1+stepz,lmax,1+stepz write(1,"(A4,I2,A3,ES14.6)") ' c0(',l,') =',c0(l) end do end if if (nonaxi) then if (.not.refy) then if (verbose) then write(*,*) do l=1,lmax do m=1,l if (stepz*mod(l-m,2).eq.0) then write(*,"(A4,I2,A1,I2,A4,ES14.6)") ' cc(', . l,',',m,') =',cc(l,m) end if end do end do write(*,*) do l=1,lmax do m=1,l if (stepz*mod(l-m,2).eq.0) then write(*,"(A4,I2,A1,I2,A4,ES14.6)") ' cs(', . l,',',m,') =',cs(l,m) end if end do end do end if if (logfile) then write(1,*) do l=1,lmax do m=1,l if (stepz*mod(l-m,2).eq.0) then write(1,"(A4,I2,A1,I2,A4,ES14.6)") ' cc(', . l,',',m,') =',cc(l,m) end if end do end do write(1,*) do l=1,lmax do m=1,l if (stepz*mod(l-m,2).eq.0) then write(1,"(A4,I2,A1,I2,A4,ES14.6)") ' cs(', . l,',',m,') =',cs(l,m) end if end do end do end if 10 format(I4,I4,A2,ES14.6,A1,ES14.6) else if (verbose) then write(*,*) do l=1,lmax do m=1,l if (stepz*mod(l-m,2).eq.0) then write(*,"(A4,I2,A1,I2,A4,ES14.6)") ' cc(', . l,',',m,') =',cc(l,m) end if end do end do end if if (logfile) then write(1,*) do l=1,lmax do m=1,l if (stepz*mod(l-m,2).eq.0) then write(1,"(A4,I2,A1,I2,A4,ES14.6)") ' cc(', . l,',',m,') =',cc(l,m) end if end do end do end if 20 format(I4,I4,A2,ES14.6) end if end if ! Close logfile write(*,*) if (logfile) then write(1,*) write(1,*) 'END OF LOG FILE' write(1,*) close(1) end if end if ! No surface found. else if (myproc.eq.0) then ! Open logfile. if (logfile) then open(1,file=logf,form='formatted', . status='old',position='append') end if ! Write messages. write(*,*) write(*,*) 'AHFinder: No horizon found.' write(*,*) if (.false.) then if (logfile) then write(1,*) write(1,*) 'AHFinder: No horizon found.' write(1,*) end if ! Close logfile write(*,*) if (logfile) then write(1,*) write(1,*) 'END OF LOG FILE' write(1,*) close(1) end if endif end if end if ! **************************** ! *** WRITE DATA FILES *** ! **************************** if (.false.) then if (myproc.eq.0) then c Expansion coefficients. if (ahf_ncall.eq.1) then open(1,file=almf,form='formatted', . status='replace') write(1,"(A21)") '# Radial coefficients' write(1,"(A1)") '#' write(1,"(A11,I4)") '# Time step',cctk_iteration write(1,"(A6,ES14.6)") '# Time',cctk_time write(1,"(A6,I4)") '# Call',ahf_ncall else open(1,file=almf,form='formatted', . status='old',position='append') write(1,*) write(1,"(A11,I4)") '# Time step',cctk_iteration write(1,"(A6,ES14.6)") '# Time',cctk_time write(1,"(A6,I4)") '# Call',ahf_ncall end if if (status) then if (mtype.eq.1) then if (horizon) then write(1,"(A30,I4)") '# Surface found: Outer horizon' else write(1,"(A31,I4)") '# Surface found: Outer horizon?' end if else if (mtype.eq.2) then if (horizon) then write(1,"(A30,I4)") '# Surface found: Inner horizon' else write(1,"(A31,I4)") '# Surface found: Inner horizon?' end if else if (mtype.eq.3) then if (horizon) then write(1,"(A34,I4)") '# Surface found: Marginal horizon?' else write(1,"(A32,I4)") '# Surface found: Trapped surface' end if else if (mtype.eq.4) then if (horizon) then write(1,"(A34,I4)") '# Surface found: Marginal horizon?' else write(1,"(A30,I4)") '# Surface found: Not a horizon' end if end if else write(1,"(A18,I4)") '# No surface found' end if if (status) then write(1,"(A7,ES14.6)") '# xc =',xc write(1,"(A7,ES14.6)") '# yc =',yc write(1,"(A7,ES14.6)") '# zc =',zc write(1,"(A1)") '#' write(1,"(A22)") '# a_lm l m' write(1,"(A1)") '#' write(1,"(ES14.6,2I4)") c0(0),0,0 do l=1,lmax do m=l,1,-1 write(1,"(ES14.6,2I4)") cs(l,m),l,-m end do write(1,"(ES14.6,2I4)") c0(l),l,0 do m=1,l write(1,"(ES14.6,2I4)") cc(l,m),l,m end do end do end if close(1) ! Area. if (ahf_ncall.eq.1) then open(1,file=areaf,form='formatted', . status='replace') write(1,"(A14)") '" Horizon area' else open(1,file=areaf,form='formatted', . status='old',position='append') end if if (status) write(1,"(2ES14.6)") cctk_time,intarea_h close(1) ! Mass. if (ahf_ncall.eq.1) then open(1,file=massf,form='formatted', . status='replace') write(1,"(A14)") '" Horizon mass' else open(1,file=massf,form='formatted', . status='old',position='append') end if if (status) write(1,"(2ES14.6)") cctk_time,0.141047396D0*dsqrt(intarea_h) close(1) ! Circumferences. if (ahf_ncall.eq.1) then open(1,file=circeqf,form='formatted', . status='replace') write(1,"(A26)") '" Equatorial circumference' else open(1,file=circeqf,form='formatted', . status='old',position='append') end if if (status) write(1,"(2ES14.6)") cctk_time,circ_eq close(1) if (ahf_ncall.eq.1) then open(1,file=merip1f,form='formatted', . status='replace') write(1,"(A28)") '" Length of meridian, phi=0' else open(1,file=merip1f,form='formatted', . status='old',position='append') end if if (status) write(1,"(2ES14.6)") cctk_time,meri_p1 close(1) if (ahf_ncall.eq.1) then open(1,file=merip2f,form='formatted', . status='replace') write(1,"(A31)") '" Length of meridian, phi=pi/2' else open(1,file=merip2f,form='formatted', . status='old',position='append') end if if (status) write(1,"(2ES14.6)") cctk_time,meri_p2 close(1) end if endif ! ****************************************************** ! *** PREPARING {ahfgrid3,ahf_exp3} FOR OUTPUT *** ! ****************************************************** if (find3) then call AHFinder_find3(CCTK_FARGUMENTS,mtype) end if ! ********************************* ! *** OUTPUT GRID FUNCTIONS *** ! ********************************* if (ahf_2Doutput.ne.0) then if (find3) then if (mfind.eq.2) then call CCTK_OutputVarAsByMethod(ierror,cctkGH, . "ahfinder::ahfgrid3","IOFlexIO_2D","ahfgrid") call CCTK_OutputVarAsByMethod(ierror,cctkGH, . "ahfinder::ahf_exp3","IOFlexIO_2D","ahf_exp") ! call AHFinder_2dio(ahfgrid3,ahf_ncall) ! call AHFinder_2dio(ahf_exp3,ahf_ncall) ! ! call io_write2d(ahfgrid3) ! call io_write2d(ahf_exp3) ! ! call AHFinder_2dio(ahfgrid3,0) ! call AHFinder_2dio(ahf_exp3,0) end if else call CCTK_OutputVarAsByMethod(ierror,cctkGH, . "ahfinder::ahfgrid","IOFlexIO_2D","ahfgrid") call CCTK_OutputVarAsByMethod(ierror,cctkGH, . "ahfinder::ahf_exp","IOFlexIO_2D","ahf_exp") ! call AHFinder_2dio(ahfgrid,ahf_ncall) ! call AHFinder_2dio(ahf_exp,ahf_ncall) ! ! call io_write2d(ahfgrid) ! call io_write2d(ahf_exp) ! ! call AHFinder_2dio(ahfgrid,0) ! call AHFinder_2dio(ahf_exp,0) end if end if if (ahf_3Doutput.ne.0) then ! This is commented out because I still need ! to write the c routine "AHFinder_3dio". ! call AHFinder_2dio(ahfgrid,ahf_ncall) ! call AHFinder_2dio(ahf_exp,ahf_ncall) ! call io_write2d(ahfgrid) ! call io_write2d(ahf_exp) ! call AHFinder_2dio(ahfgrid,0) ! call AHFinder_2dio(ahf_exp,0) end if ! **************** ! *** MASK *** ! **************** if (horizon.and.(ahf_mask.ne.0)) then call AHFinder_mask(CCTK_FARGUMENTS) else if ((mtype.eq.1).and.(ahf_weakmask.ne.0)) then call AHFinder_mask(CCTK_FARGUMENTS) else ahmask = one end if call CCTK_SyncGroup(cctkGH,"ahfinder::ahfmask") ! ************************************************** ! *** TRANSFORM PHYSICAL TO CONFORMAL METRIC *** ! ************************************************** call CCTK_WARN(4,"Converting metric: physical -> conformal") gxx=gxx/psi**4 gxy=gxy/psi**4 gxz=gxz/psi**4 gyy=gyy/psi**4 gyz=gyz/psi**4 gzz=gzz/psi**4 conformal_state = CONFORMAL_METRIC ! ********************************************** ! *** UPDATING COUNTER AND LOOP IF FIND3 *** ! ********************************************** if (find3) then if (mfind.eq.2) then mfind = 0 else mfind = mfind + 1 goto 1 end if end if ! *************** ! *** END *** ! *************** end subroutine AHFinder