c/*@@ c @file AHFinder_flow.F c @date November 1998 c @author Miguel Alcubierre c @desc c Master routine to control the iterations for c the flow algorithm. c c The flow algorithm used here is that of Carsten Gunlach: c Phys. Rev. D 57, 863 (1998). c c @enddesc c@@*/ c Note that including cactus.h will also include AHFinder.h !#include "cactus.h" #include "cctk.h" #include "cctk_parameters.h" #include "cctk_arguments.h" subroutine AHFinder_flow(CCTK_FARGUMENTS,rmx,status,logf) use AHFinder_dat implicit none DECLARE_CCTK_FARGUMENTS DECLARE_CCTK_PARAMETERS logical status,new CCTK_INT ITMAX integer i,j,l,m CCTK_REAL rmx,ahftol CCTK_REAL A,B CCTK_REAL ahfsum,ahfdiff,ahfdiff_old CCTK_REAL maxchange,minchange CCTK_REAL reldiff CCTK_REAL spher CCTK_REAL zero,one CCTK_REAL dd,aux character*200 logf ! Description of variables: ! ! i,l,m Counters. ! ! ITMAX Maximum number of iterations. ! ! A,B Flow parameters. ! ! hw Weigth of `H' flow (see Carsten's paper). ! cw Weigth of `C' flow (see Carsten's paper). ! nw Weigth of `N' flow (see Carsten's paper). ! ! ahftol Tolerance used to decide if we have converged. ! *************************** ! *** DEFINE NUMBERS *** ! *************************** zero = 0.0D0 one = 1.0D0 dd = min(dx,dy,dz) ! *************************** ! *** FIND PARAMETERS *** ! *************************** ! A = getr8("ahf_flowa") ! B = getr8("ahf_flowb") A = ahf_flowa B = ahf_flowb ! hw = getr8("ahf_flowh") ! cw = getr8("ahf_flowc") ! nw = getr8("ahf_flown") hw = ahf_flowh cw = ahf_flowc nw = ahf_flown ! maxchange = getr8("ahf_maxchange") ! minchange = getr8("ahf_minchange") maxchange = ahf_maxchange minchange = ahf_minchange ! ************************* ! *** INITIAL GUESS *** ! ************************* 5 continue ! Initialize arrays. c0 = zero cc = zero cs = zero hflow0 = zero hflowc = zero hflows = zero cflow0 = zero cflowc = zero cflows = zero nflow0 = zero nflowc = zero nflows = zero ! As initial guess I take the largest sphere that ! can fit in the grid. if ((myproc.eq.0).and.verbose) then write(*,*) write(*,"(A35,ES11.3)") ' Initial guess sphere with radius: ',rmx write(*,"(A13,3ES11.3)") ' centered on:',real(xc),real(yc),real(zc) end if c0(0) = rmx ! ******************************* ! *** MAIN ITERATION LOOP *** ! ******************************* ! Find ITMAX and ahftol. ! ITMAX = geti("ahf_flowiter") ! ahftol = getr8("ahf_flowtol") ITMAX = ahf_flowiter ahftol = ahf_flowtol ! Start the iterations. do i=1,ITMAX ! Save the old values of the coefficients. c0_old = c0 cc_old = cc cs_old = cs ! In order to adapt the stepsize, I do two steps ! with the current stepsize, and compare them with ! one step with double the stepsize. do j=1,2 ahfsum = zero ahfdiff = zero ! ************************************ ! *** FIND SPECTRAL COMPONENTS *** ! ************************************ call AHFinder_fun(CCTK_FARGUMENTS) call AHFinder_exp(CCTK_FARGUMENTS) call AHFinder_int(CCTK_FARGUMENTS) ! ***************************** ! *** CHECK ERROR FLAGS *** ! ***************************** ! If interror is true there was an error in the integral. ! There are two possibilities here: ! 1) This is the first of the two small steps (j=1). ! We then went out of bounds in the previous iteration, ! so there is nothing we can do to fix it. if (interror.and.(j.eq.1)) then ! We can be out of bounds for two reasons: ! a) We are out of the computational domain, pitty. if (interror2.ne.0) then write(*,*) write(*,*) 'AHFinder: Out of bounds, giving up.' return ! b) The radius is negative somewhere. This is a really ! interesting case, since it probably means that the ! origin of our expansion is outside the horizon. ! This can happen is the black-hole is not centered, ! or if we have more than one horizon. else if (interror1.ne.0) then write(*,*) write(*,*) 'AHFinder: Negative radius.' return ! We should never get here! else write(*,*) write(*,*) 'Something is very wrong!' return end if ! 2) This is the second of the two small steps (j=2). ! We can still try to reduce the stepsize and start ! the iteration again. We then jump to the place ! where the stepsize is adapted, and force it to ! be reduced. else if (interror) then interror = .false. reldiff = one ahfsum = one ahfdiff = one intexp2 = one goto 100 end if ! **************************************** ! *** UPDATE SPECTRAL COEFFICIENTS *** ! **************************************** ! Update c0. do l=0,lmax,1+stepz aux = c0(l) c0(l) = aux - A/(one + B*dble(l)*(dble(l) + one)) . *(hw*hflow0(l) + cw*cflow0(l) + nw*nflow0(l)) ahfsum = ahfsum + c0_old(l)**2 ahfdiff = ahfdiff + (c0(l) - c0_old(l))**2 end do ! Update {cc,cs}. if (nonaxi) then ! Update cc(l,m). do l=1,lmax do m=1+stepx,l,1+stepx if (stepz*mod(l-m,2).eq.0) then aux = cc(l,m) cc(l,m) = aux - A/(one . + B*dble(l)*(dble(l) + one)) . *(hw*hflowc(l,m) + cw*cflowc(l,m) . + nw*nflowc(l,m)) ahfsum = ahfsum + cc_old(l,m)**2 ahfdiff = ahfdiff + (cc(l,m)-cc_old(l,m))**2 end if end do end do ! Update 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) then cs(l,m) = cs(l,m) - A/(one . + B*dble(l)*(dble(l) + one)) . *(hw*hflows(l,m) + cw*cflows(l,m) . + nw*nflows(l,m)) ahfsum = ahfsum + cs_old(l,m)**2 ahfdiff = ahfdiff + (cs(l,m)-cs_old(l,m))**2 end if end do end do end if end if ahfsum = dsqrt(ahfsum) ahfdiff = dsqrt(ahfdiff) ! If this is the first of the two small steps, save ! twice the difference. Why twice? Because this would ! have been the resulting difference if we had done only ! one step that was twice as big. if (j.eq.1) then ahfdiff_old = 2.0D0*ahfdiff end if end do ! *************************** ! *** ADJUST STEPSIZE *** ! *************************** ! Here we find the relative difference between having done two ! small steps and one large step. reldiff = dabs(ahfdiff - ahfdiff_old) . / dabs(ahfdiff + ahfdiff_old) 100 new = .false. ! If the relative difference is too large, we reduce the ! stepsize and ignore this iteration. if ((reldiff.gt.maxchange).or. . (ahfdiff.gt.maxchange*ahfsum)) then new = .true. A = 0.5D0*A c0 = c0_old cc = cc_old cs = cs_old ! If the diference is too small, we can safely increase ! the stepsize for the next iteration. else if (reldiff.lt.minchange) then new = .true. A = 2.0D0*A end if ! ****************************************** ! *** LOGFILE AND MESSAGES TO SCREEN *** ! ****************************************** ! Open logfile. logf = filestr(1:nfile)//"ahf_logfile" if (logfile.and.(myproc.eq.0)) then open(1,file=logf,form='formatted',status='old', . position='append') end if ! Write messages. if (myproc.eq.0) then if (veryver) then write(*,*) write(*,*) write(*,"(A28,I3)") ' AHFinder: FLOW ITERATION ',i if (intarea.ne.zero) then aux = one/intarea else aux = one end if write(*,*) write(*,"(A21,ES14.6)") ' Surface area =',intarea write(*,"(A21,ES14.6)") ' Mean value of H =',aux*intexp write(*,"(A21,ES14.6)") ' Mean value of H^2 =',aux*intexp2 write(*,"(A21,ES14.6)") ' ahfdiff/ahfsum =',ahfdiff/ahfsum if (new) then write(*,*) write(*,"(A15,ES14.6)") ' New stepsize =',A end if 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(*,"(A20)") ' Shape coefficients:' 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,"(A28,I3)") ' AHFinder: FLOW ITERATION ',i if (intarea.ne.zero) then aux = one/intarea else aux = one end if write(1,*) write(1,"(A21,ES14.6)") ' Surface area =',intarea write(1,"(A21,ES14.6)") ' Mean value of H =',aux*intexp write(1,"(A21,ES14.6)") ' Mean value of H^2 =',aux*intexp2 write(1,"(A21,ES14.6)") ' ahfdiff/ahfsum =',ahfdiff/ahfsum if (new) then write(1,*) write(1,"(A15,ES14.6)") ' New stepsize =',A end if 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,"(A20)") ' Shape coefficients:' 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 (veryver) 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 (veryver) 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 end if ! Close logfile. if (logfile.and.(myproc.eq.0)) then close(1) end if ! *********************************** ! *** CHECK IF WE HAVE FAILED *** ! *********************************** ! The finder needs a way to give up if there is no ! horizon. At the moment, the way in which I do this ! is that I just give up if the monopole term is of ! the order of 2 grid points in size. Maybe this ! could be improved on? if (c0(0).lt.2.0D0*dd) then write(*,*) write(*,"(A43)") ' AHFinder: c0(0) is too small, giving up.' return end if ! ******************************************** ! *** CHECK IF WE HAVE FOUND A HORIZON *** ! ******************************************** ! For this I use a straigthforward test to check if the ! change in the surface is below a given tolerance. This ! usually works well, but if the tolerance is too small then ! the convergence becomes very slow. Also, since I use an ! adaptive stepsize, I run the risk of having a small change ! simply because the stepsize itself is very small, and not ! because we found a horizon. This can happen is the tolerance ! is not small enough. So we have a balancing act here. if (ahfdiff.lt.ahftol*ahfsum) then status = .true. return end if end do ! ****************************** ! *** TOO MANY IERATIONS *** ! ****************************** ! If we ever get here, we did too many iterations and ! failed to converge. status = .true. if (veryver) then write(*,*) write(*,"(A32,I2)") ' AHFinder: Too many iterations' end if ! *************** ! *** END *** ! *************** end subroutine AHFinder_flow real*8 function spher(theta,phi) ! Function to reconstruct a function pointwise ! from the expansion coefficients {c0,cc,cs}. use AHFinder_dat implicit none integer l,m REAL LEGEN REAL theta,phi REAL F,cost,cosa,sina REAL aux ! ****************************** ! *** Axisymmetric terms *** ! ****************************** F = 0 cost = cos(theta) do l=0,lmax F = F + c0(l)*legen(l,0,cost) end do ! ********************************** ! *** Non-axisymmetric terms *** ! ********************************** do m=1,lmax aux = m*phi sina = sin(aux) cosa = cos(aux) do l=m,lmax F = F + LEGEN(l,m,cost)*(cc(l,m)*cosa + cs(l,m)*sina) end do end do spher = F ! *************** ! *** END *** ! *************** end function spher