/*@@ @file AHFinder_flow.F @date November 1998 @author Miguel Alcubierre @desc Master routine to control the iterations for the flow algorithm. The flow algorithm used here is that of Carsten Gunlach: Phys. Rev. D 57, 863 (1998). @enddesc @@*/ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" subroutine AHFinder_flow(CCTK_ARGUMENTS,rmx,status,logf) use AHFinder_dat implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS logical status,new,update integer i,j,l,m integer itmax CCTK_REAL rmx,ahftol CCTK_REAL A,B CCTK_REAL ahfsum,ahfdiff,ahfdiff_old CCTK_REAL reldiff CCTK_REAL maxchange,minchange CCTK_REAL zero,one CCTK_REAL dd,aux character(len=200) :: logf ! Description of variables: ! ! new Flag to indicate that stepsize has changed. ! update Flag to indicate if coefficients will be updated. ! ! i,j,l,m Counters. ! ! itmax Maximum number of iterations. ! ! ahftol Tolerance used to decide if we have converged. ! ! 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). ! ! ahfsum Square root of the sum of the squares of the ! expansion coefficients. ! ! ahfdiff Square root of the sum of the squares of the ! change in the coefficients. ! ! rediff Relative difference between having done one large ! step and two small ones. ! ! maxchange A value of "reldiff" larger than this forces a ! reduction of the stepsize. ! ! minchange A value of "reldiff" smaller than this forces an ! increase of the stepsize. ! ! dd Maximum of (dx,dy,dz). ! *************************** ! *** DEFINE NUMBERS *** ! *************************** zero = 0.0D0 one = 1.0D0 dd = max(dx,dy,dz) ! *************************** ! *** FIND PARAMETERS *** ! *************************** A = ahf_flowa B = ahf_flowb hw = ahf_flowh cw = ahf_flowc nw = ahf_flown maxchange = ahf_maxchange minchange = ahf_minchange ! ******************************* ! *** MAIN ITERATION LOOP *** ! ******************************* ! Find itmax and ahftol. 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_ARGUMENTS) call AHFinder_exp(CCTK_ARGUMENTS) call AHFinder_int(CCTK_ARGUMENTS) ! ***************************** ! *** CHECK ERROR FLAGS *** ! ***************************** update = .true. ! 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 three reasons: ! a) The radius is negative somewhere. This can mean ! one of two things: ! ! * There is no horizon, and the flow has reached ! the origin. ! ! * The origin for our expansion is outside the horizon. ! This can happen if the black-hole is not centered, ! or if we have more than one horizon. if (interror1.ne.0) then write(*,*) write(*,*) 'Negative radius.' return ! b) We are out of the computational domain, pity. else if (interror2.ne.0) then write(*,*) write(*,*) 'Out of bounds, giving up.' return ! c) We are inside the mask. else if (interror3.ne.0) then write(*,*) write(*,*) 'Inside mask, giving up.' return ! We should never get here! else write(*,*) write(*,*) 'PANIC! 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 return the coefficients ! to their original values and, and do not update them. ! ! Since I will not update the coefficients, I need to ! assign values to "afhsum" and "ahfdiff". In order ! to force a reduction of the step size, I assign ! a very large number to "ahfdiff", which will make ! "reldiff" close to 1.0, far too large. I also ! assign a value of 1.0 to "ahfsum", so it is much ! smaller than "ahfdiff". This is important in order ! to prevent the algorithm from thinking we are done ! (see below for termination criteria). else if (interror) then interror = .false. c0 = c0_old cc = cc_old cs = cs_old ahfsum = one ahfdiff = 1.0D10 update = .false. end if ! **************************************** ! *** UPDATE SPECTRAL COEFFICIENTS *** ! **************************************** if (update) then ! 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 = sqrt(ahfsum) ahfdiff = sqrt(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 if end do ! ************************************ ! *** FIND RELATIVE DIFFERENCE *** ! ************************************ ! Here we find the relative difference between having done two ! small steps and one large step. reldiff = abs(ahfdiff - ahfdiff_old) . / (abs(ahfdiff) + abs(ahfdiff_old)) ! Open logfile. logf = filestr(1:nfile)//"/ahf_logfile" if (logfile.and.(myproc.eq.0)) then open(11,file=logf,form='formatted',status='old', . position='append') end if ! *************************** ! *** ADJUST STEPSIZE *** ! *************************** new = .false. ! If the relative difference is too large, we reduce the ! stepsize for the next iteration. if ((reldiff.gt.maxchange).or. . (ahfdiff.gt.maxchange*ahfsum)) then new = .true. A = 0.5D0*A if (veryver) then write(*,*) write(*,"(A15,ES14.6)") ' New stepsize =',A end if if (logfile) then write(11,*) write(11,"(A15,ES14.6)") ' New stepsize =',A end if ! If the difference 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 if (veryver) then write(*,*) write(*,"(A15,ES14.6)") ' New stepsize =',A end if if (logfile) then write(11,*) write(11,"(A15,ES14.6)") ' New stepsize =',A end if end if ! ****************************************** ! *** LOGFILE AND MESSAGES TO SCREEN *** ! ****************************************** ! Write messages. if (myproc.eq.0) then if (veryver) then write(*,*) write(*,*) write(*,"(A16,I3)") ' 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 (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(11,*) write(11,"(A16,I3)") ' FLOW ITERATION ',i if (intarea.ne.zero) then aux = one/intarea else aux = one end if write(11,*) write(11,"(A21,ES14.6)") ' Surface area =',intarea write(11,"(A21,ES14.6)") ' Mean value of H =',aux*intexp write(11,"(A21,ES14.6)") ' Mean value of H^2 =',aux*intexp2 write(11,"(A21,ES14.6)") ' ahfdiff/ahfsum =',ahfdiff/ahfsum if (offset) then write(11,*) write(11,"(A6,ES14.6)") ' xc =',xc write(11,"(A6,ES14.6)") ' yc =',yc write(11,"(A6,ES14.6)") ' zc =',zc end if write(11,*) write(11,"(A20)") ' Shape coefficients:' write(11,*) write(11,"(A4,I2,A3,ES14.6)") ' c0(',0,') =',c0(0) do l=1+stepz,lmax,1+stepz write(11,"(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(11,*) do l=1,lmax do m=1,l if (stepz*mod(l-m,2).eq.0) then write(11,"(A4,I2,A1,I2,A4,ES14.6)") ' cc(', . l,',',m,') =',cc(l,m) end if end do end do write(11,*) do l=1,lmax do m=1,l if (stepz*mod(l-m,2).eq.0) then write(11,"(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(11,*) do l=1,lmax do m=1,l if (stepz*mod(l-m,2).eq.0) then write(11,"(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(11) 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(*,*) 'Surface radius 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(*,*) 'Too many iterations' end if ! *************** ! *** END *** ! *************** end subroutine AHFinder_flow