path: root/src/AHFinder_pow.F

/*@@
  @file      AHFinder_pow.F
  @date      April 1998
  @author    Miguel Alcubierre
  @desc 
             Multidimensional minimization routines taken
             from Numerical Recipes.  The routines have some
             modifications to adapt them to the specific problem
             we are trying to solve.
  @enddesc 
@@*/

#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"


      module F1COM

!     Data to be shared by minimization routines.

      implicit none

      logical error

      integer NCOM

      CCTK_REAL XMIN
      CCTK_REAL, allocatable, dimension(:) :: PCOM,XICOM


!     ***************
!     ***   END   ***
!     ***************

      end module F1COM




***********************************************************************
***********************************************************************




      subroutine POWELL(CCTK_ARGUMENTS,P,XI,N,FTOL,ITER,ITMAX,FRET,found)

!     Minimization of a function FUNC on N variables (FUNC is
!     not an argument, it is a fixed function name).  Input consists
!     of an initial starting point P that is a vector of length N;
!     an initial matrix XI of dimensions NxN and whose columns contain
!     the initial set of directions (usually the N unit vectors); and
!     FTOL, the fractional tolerance in the function value such that
!     failure to decrease by more than this ammount on one iteration
!     signals doneness.  On output, P is set to the best point found,
!     XI is the then-current direction set, FRET is the returned
!     function value at P, and ITER is the number of iterations taken.
!
!     Extra output added by myself (Miguel Alcubierre) is a logical
!     variable "found" that is set to true is a minimum was indeed
!     found, and to false if it was not.

      use AHFinder_dat

      implicit none

      DECLARE_CCTK_ARGUMENTS

      logical found

      integer i,j,N
      integer ITER,ITMAX,IBIG

      CCTK_REAL FUNC,DEL,FP,FPTT,FRET,FTOL,T,TOL1,ZEPS

      CCTK_REAL, dimension(1:N) :: P,PT,PTT,XIT
      CCTK_REAL, dimension(1:N,1:N) :: XI

      character(len=200) :: logf

      parameter (ZEPS=1.0D-10)


!     *************************
!     ***   START ROUTINE   ***
!     *************************

      FRET = FUNC(CCTK_ARGUMENTS,P,N)

!     Save initial point.

      do j=1,N
         PT(j) = P(j)
      end do

      ITER = 0

   10 ITER = ITER + 1

      if ((myproc.eq.0).and.veryver) then
         write(*,*)
         write(*,"(A20,I2)") 'POWELL ITERATION  = ',ITER
      end if

      if ((myproc.eq.0).and.logfile) then
         logf = filestr(1:nfile)//"/ahf_logfile"
         open(11,file=logf,form='formatted',status='old',
     .   position='append')
         write(11,*)
         write(11,"(A20,I2)") 'POWELL ITERATION  = ',ITER
         close(11)
      end if

      FP = FRET
      IBIG = 0

!     Initialize DEL.  It will be the biggest function decrease.

      DEL = 0.0D0

!     In each iteration, loop over all directions in the set.

      do i=1,N

!        Copy the direction.

         do j=1,N
            XIT(j) = XI(j,i)
         end do

         FPTT = FRET

!        Minimize along it.

         call LINMIN(CCTK_ARGUMENTS,P,XIT,N,FRET,FTOL)

!        And record it if it is the largest decrease so far.

         if (abs(FPTT-FRET).gt.DEL) then
            DEL  = abs(FPTT-FRET)
            IBIG = I
         end if

   20   continue

      end do

!     Termination criterion.

      TOL1 = 0.5D0*FTOL*(abs(FP)+abs(FRET)) + ZEPS

      if (abs(FP-FRET).le.TOL1) then
         found = .true.
         return
      else if (abs(FRET).lt.ZEPS) then
         found = .true.
         return
      end if

!     Too many iterations?

      if (ITER.ge.ITMAX) return

!     Construct the extrapolated point and the average direction.

      do j=1,N
         PTT(j) = 2.0D0*P(j) - PT(j)
         XIT(j) = P(j) - PT(j)
         PT(j)  = P(j)
      end do

!     Function value at extrapolated point.

      FPTT = FUNC(CCTK_ARGUMENTS,PTT,N)

!     First reason not to use new direction.

      if (FPTT.ge.FP) goto 10

!     Second reason not to use new direction.

      T = 2.0D0*(FP-2.0D0*FRET+FPTT)*(FP-FRET-DEL)**2
     .  - DEL*(FP-FPTT)**2

      if (T.gt.0.0D0) goto 10

!     Move to the minimum of the new direction.

      call LINMIN(CCTK_ARGUMENTS,P,XIT,N,FRET,FTOL)

!     Save new direction.

      do j=1,N
         XI(j,IBIG) = XIT(j)
      end do

!     Back for another iteration.

      goto 10


!     ***************
!     ***   END   ***
!     ***************

      end subroutine POWELL




***********************************************************************
***********************************************************************



      subroutine LINMIN(CCTK_ARGUMENTS,P,XI,N,FRET,TOL)

!     Given an N dimensional point P and an N dimensional direction XI,
!     moves and resets P to where the function FUNC takes on a
!     minimum along the direction XI from P, and replaces XI by the
!     actual vector displacement that P was moved.  Also returns as
!     FRET the value of FUNC at the returned location P.  This is
!     actually all acomplished by calling the routines MNBRAK and
!     PARABOLA.

      use F1COM

      implicit none

      DECLARE_CCTK_ARGUMENTS

      integer i,j,N
      integer l,m

      CCTK_REAL XA,XB,XX,FA,FB,FX
      CCTK_REAL FRET,PARABOLA,TOL
      CCTK_REAL zero
      CCTK_REAL dx,dy,dz

      CCTK_REAL, dimension(1:N) :: P,XI

      CCTK_REAL F1DIM

      EXTERNAL  F1DIM


!     **********************************************
!     ***   ALLOCATE STORAGE FOR SHARED ARRAYS   ***
!     **********************************************

      zero = 0.0d0

      allocate(PCOM(1:N),XICOM(1:N))


!     *************************
!     ***   START ROUTINE   ***
!     *************************

      dx = cctk_delta_space(1)
      dy = cctk_delta_space(2)
      dz = cctk_delta_space(3)

      NCOM = N

      do j=1,N
         PCOM(j)  = P(j)
         XICOM(j) = XI(j)
      end do

!     Initial guess for brackets.

      XA = 0.0D0
      XX = 0.1D0*min(dx,dy,dz)

!     Bracket the minimum.

      call MNBRAK(CCTK_ARGUMENTS,XA,XX,XB,FA,FX,FB,F1DIM)

!     Check if function is constant.

      if ((FA.ne.FX).or.(FB.ne.FX)) then

!        Call 1D minimization.

         FRET = PARABOLA(CCTK_ARGUMENTS,XA,XX,XB,FA,FX,FB,F1DIM,TOL)

!        Construct the vector results to return.

         XI = XMIN*XI
         P  = P + XI

      end if

      
!     ******************************
!     ***   DEALLOCATE STORAGE   ***
!     ******************************

      deallocate(PCOM,XICOM)


!     ***************
!     ***   END   ***
!     ***************

      end subroutine LINMIN




***********************************************************************
***********************************************************************




      CCTK_REAL function F1DIM(CCTK_ARGUMENTS,XX)

!     Must accompany LINMIN

      use F1COM

      implicit none

      DECLARE_CCTK_ARGUMENTS

      integer j

      CCTK_REAL XX,FUNC

      CCTK_REAL, dimension(1:NCOM) :: XT


!     *************************
!     ***   START ROUTINE   ***
!     *************************

      do j=1,NCOM
         XT(j) = PCOM(j) + XX*XICOM(j)
      end do

      F1DIM = FUNC(CCTK_ARGUMENTS,XT,NCOM)


!     ***************
!     ***   END   ***
!     ***************

      end function F1DIM




***********************************************************************
***********************************************************************



      subroutine MNBRAK(CCTK_ARGUMENTS,XA,XB,XC,FA,FB,FC,FUNC)

!     Given a function FUNC, and given distinct initial points
!     XA and XB, this routine searches in the downhill direction
!     (defined by the function as evaluated at the initial points)
!     and returns new points XA, XB, XC which bracket a minimum
!     of the function.  Also returned are the function values
!     at the three points, FA, FB and FC.

      implicit none

      DECLARE_CCTK_ARGUMENTS

      CCTK_REAL XA,XB,XC
      CCTK_REAL FA,FB,FC,FU,FUNC
      CCTK_REAL DUM,Q,S,U,ULIM
      CCTK_REAL GOLD,GLIMIT,TINY

      parameter (GOLD=1.618034D0,GLIMIT=4.0D0,TINY=1.0D-20)


!     *************************
!     ***   START ROUTINE   ***
!     *************************

      FA = FUNC(CCTK_ARGUMENTS,XA)
      FB = FUNC(CCTK_ARGUMENTS,XB)

!     Switch roles of A and B so that we can go downhill
!     in the direction from A to B.

      if (FB.gt.FA) then

         DUM = XA
         XA  = XB
         XB  = DUM

         DUM = FB
         FB  = FA
         FA  = DUM

      end if

!     First guess for C.

      XC = XB + GOLD*(XB-XA)

      FC = FUNC(CCTK_ARGUMENTS,XC)

!     Check if the function is constant.

      if ((FA.eq.FB).and.(FB.eq.FC)) then
         return
      end if

!     Keep returning here until we bracket.

      do while(FB.gt.FC)

!        Compute U by parabollic extrapolation from A, B, C.
!        TINY is used to prevent any possible division by zero.

         S = (XB-XA)*(FB-FC)
         Q = (XB-XC)*(FB-FA)
         U = XB - ((XB-XC)*Q - (XB-XA)*S)
     .     / (2.0D0*SIGN(MAX(ABS(Q-S),TINY),Q-S))

!        We will not go further than this.

         ULIM = XB + GLIMIT*(XC-XB)

!        Now to test various possibilities.

         if ((XB-U)*(U-XC).gt.0.0D0) then

!           Parabolic U is between B and C, try it.

            FU = FUNC(CCTK_ARGUMENTS,U)

            if (FU.lt.FC) then

!              Got a minimum between B and C.

               XA = XB
               FA = FB
               XB = U
               FB = FU
               return

            else if (FU.gt.FB) then

!              Got a minimum between A and U.

               XC = U
               FC = FU
               return

            end if

!           Parabolic fit was no use.  Use default magnification.

            U = XC + GOLD*(XC-XB)

            FU = FUNC(CCTK_ARGUMENTS,U)

         else if ((XC-U)*(U-ULIM).gt.0.0D0) then

!           Parabolic fit is between C and its allowed limit.

            FU = FUNC(CCTK_ARGUMENTS,U)

            if (FU.lt.FC) then

               XB = XC
               XC = U
               U  = XC + GOLD*(XC-XB)
               FB = FC
               FC = FU

               FU = FUNC(CCTK_ARGUMENTS,U)

            end if

         else if ((U-ULIM)*(ULIM-XC).ge.0.0D0) then

!           Limit parabolic U to maximum allowed value.

            U = ULIM

            FU = FUNC(CCTK_ARGUMENTS,U)

         else

!           Reject parabolic U, use default magnification.

            U = XC + GOLD*(XC-XB)
            FU = FUNC(CCTK_ARGUMENTS,U)

         end if

!        Eliminate oldest point and continue.

         XA = XB
         XB = XC
         XC = U
         FA = FB
         FB = FC
         FC = FU

      end do


!     ***************
!     ***   END   ***
!     ***************

      end subroutine MNBRAK




***********************************************************************
***********************************************************************




      CCTK_REAL function PARABOLA(CCTK_ARGUMENTS,XA,XB,XC,FXA,FXB,FXC,
     .  FUNC,TOL)

!     Given a function F, and given a bracketing triplet of abscissas
!     XA, XB, XC (such that XB is between XA and XC, and F(XB) is less
!     than both F(XA) and F(XC)), this routine isolates the minimum
!     to a fractional precision of about TOL using inverse parabolic
!     interpolation. The abscissa of the minimum is returned as XMIN,
!     and the minimum function value is returned as PARABOLA, the
!     returned function value.

      use F1COM

      implicit none

      DECLARE_CCTK_ARGUMENTS

      integer ITER,ITMAX

      CCTK_REAL XA,XB,XC,FXA,FXB,FXC,FUNC,TOL
      CCTK_REAL A,B,D,P,Q,S,U,XX,XM,TOL1
      CCTK_REAL FA,FB,FU,FX,FP,ZEPS

      parameter (ITMAX=100,ZEPS=1.0D-10)


!     *************************
!     ***   START ROUTINE   ***
!     *************************

      error = .false.

!     A and B must be in ascending order, though the input
!     abscissas need not be.

      if (XC.lt.XA) then
         A = XC
         B = XA
         FA = FXC
         FB = FXA
      else
         A = XA
         B = XC
         FA = FXA
         FB = FXC
      end if

      XX = XB
      FX = FXB
      FP = FA

!     Main program loop.

      do ITER=1,ITMAX

         XM = 0.5D0*(A + B)

!        Test for done here.

         TOL1 = 0.5D0*TOL*(abs(FP)+abs(FX)) + ZEPS
         if (abs(FP-FX).le.TOL1) goto 3

!        Construct a trial parabolic fit.

         S = (XX-A)*(FX-FB)
         Q = (XX-B)*(FX-FA)
         P = (XX-B)*Q - (XX-A)*S
         Q = 2.0D0*(Q-S)

!        Check if Q is zero.  This can only happen if the
!        three points are colinear (in which case they where
!        not bracketing a minimum in the first place), or
!        if two of them are in fact the same point.

         if (Q.eq.0.0D0) then
            write(*,*) 'AHFinder: problem in PARABOLA'
            error = .true.
            return            
         end if

!        Find new point.

         D = - P/Q
         U = XX + D

!        Check if X was already the minimum.

         if (abs(U-XX).lt.ZEPS) goto 3

!        This is the one function evaluation per iteration.

         FU = FUNC(CCTK_ARGUMENTS,U)

!        And now we have to decide what to do with our function
!        evaluation.

         if (FU.le.FX) then

            if (U.ge.XX) then
               A = XX
               FA = FX
            else
               B = XX
               FB = FX
            end if

            XX = U

            FP = FX
            FX = FU

         else

            if (U.lt.XX) then
               A = U
               FA = FU
            else
               B = U
               FB = FU
            end if

         end if

      end do

      write(*,*) 'PARABOLA exceeded maximum number of iterations.'

    3 XMIN = XX
      PARABOLA = FX


!     ***************
!     ***   END   ***
!     ***************

      end function PARABOLA




***********************************************************************
***********************************************************************



      CCTK_REAL function FUNC(CCTK_ARGUMENTS,P,N)

!     This is a shell to find the surface integral of the expansion.

      use AHFinder_dat

      implicit none

      DECLARE_CCTK_ARGUMENTS
      DECLARE_CCTK_PARAMETERS

      logical flag

      integer i,l,m,N
      integer iter

      CCTK_REAL zero,half,one
      CCTK_REAL aux

      CCTK_REAL, dimension(1:N) :: P

      character(len=200) :: logf

      save iter,logf


!     **************************
!     ***   DEFINE NUMBERS   ***
!     **************************

      zero = 0.0D0
      half = 0.5D0
      one  = 1.0D0


!     ************************
!     ***   OPEN LOGFILE   ***
!     ************************

      logf = filestr(1:nfile)//"/ahf_logfile"

      if (logfile.and.(myproc.eq.0)) then
         open(12,file=logf,form='formatted',status='old',
     .   position='append')
      end if


!     ****************************
!     ***   FIND  {c0,cc,cs}   ***
!     ****************************

      i = 0

!     Find {xc,yc,zc}.

      if (wander) then
         if (.not.refx) then
            i = i+1
            xc = P(1)
         end if
         if (.not.refy) then
            i = i+1
            yc = P(2)
         end if
         if (.not.refz) then
            i = i+1
            zc = P(3)
         end if
      end if

!     Find c0(l).

      do l=0,lmax,1+stepz
         i = i+1
         c0(l) = P(i)
      end do

!     Find {cc,cs}.

      if (nonaxi) then

!        Find cc(l,m).

         do l=1,lmax
            do m=1+stepx,l,1+stepx
               if (stepz*mod(l-m,2).eq.0) then
                  i = i+1
                  cc(l,m) = P(i)
               end if
            end do
         end do

!        Find 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
                     i = i+1
                     cs(l,m) = P(i)
                  end if
               end do
            end do
         end if

      end if


!     *********************************************
!     ***   CALCULATE EXPANSION IF WE NEED TO   ***
!     *********************************************

!     Firstcall?

      if (firstfun) then

         iter = 1

!        Find horizon function and expansion.

         call AHFinder_fun(CCTK_ARGUMENTS)
         call AHFinder_exp(CCTK_ARGUMENTS)

!        Save old values of coefficients.

         xc_old = xc
         yc_old = yc
         zc_old = zc

         c0_old = c0
         cc_old = cc
         cs_old = cs

         firstfun = .false.

!     Not first call.

      else

         iter = iter + 1

!        Check if coefficients changed from last call.
!        I check this because the minimization algorithm
!        sometimes reevaluates the function at the same
!        point twice.  Also, if only c0(0) changed, the
!        expansion will not be affected.

         flag = .false.

         if (wander) then
            if (xc.ne.xc_old) flag = .true.
            if (yc.ne.yc_old) flag = .true.
            if (zc.ne.zc_old) flag = .true.
         end if

         do l=1+stepz,lmax,1+stepz
            if (c0(l).ne.c0_old(l)) flag = .true.
         end do

         if (nonaxi) then
            do l=1,lmax
               do m=1+stepx,l,1+stepx
                  if (stepz*mod(l-m,2).eq.0) then
                     if (cc(l,m).ne.cc_old(l,m)) flag = .true.
                  end if
               end do
               if (.not.refy) then
                  do m=1,l,1+stepx
                     if (stepz*mod(l-m,2).eq.0) then
                        if (cs(l,m).ne.cs_old(l,m)) flag = .true.
                     end if
                  end do
               end if
            end do
         end if

!        If necessary, find horizon function and expansion.

         if (flag) then
            call AHFinder_fun(CCTK_ARGUMENTS)
            call AHFinder_exp(CCTK_ARGUMENTS)
         end if

!        Save old values of coefficients.

         c0_old = c0
         cc_old = cc
         cs_old = cs

      end if

      if (myproc.eq.0) then

         if (veryver) then
            write(*,*)
            write(*,"(A13,I5)") 'FUNC CALL  =',iter
         end if

         if (logfile) then
            write(12,*)
            write(12,"(A13,I5)") 'FUNC CALL  =',iter
         end if

      end if


!     **********************************
!     ***   FIND SURFACE INTEGRALS   ***
!     **********************************

      call AHFinder_int(CCTK_ARGUMENTS)

!     Messages to screen and logfile.

      if (myproc.eq.0) then

!        Find inverse of area.  Be carefull not to divide by zero!
!        Also, a value of 1.0D10 indicates an error in the integration
!        routine, and all integrals end up with the same value.
!        So dont divide by it or I will be confused later.

         if ((intarea.ne.zero).and.(intarea.ne.1.0D10)) then
            aux = one/intarea
         else
            aux = one
         end if

!        Screen.

         if (veryver) then
            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(*,*) 'Number of interpolated points:     ',
     .         inside_min_count
            write(*,*) 'Number of those that are negative: ',
     .         inside_min_neg_count
         end if

!        Logfile.

         if (logfile) then
            write(12,*)
            write(12,"(A21,ES14.6)") ' Surface area       =',intarea
            write(12,"(A21,ES14.6)") ' Mean value of H    =',aux*intexp
            write(12,"(A21,ES14.6)") ' Mean value of H^2  =',aux*intexp2
            write(12,*) 'Number of interpolated points:     ',
     .         inside_min_count
            write(12,*) 'Number of those that are negative: ',
     .         inside_min_neg_count
         end if

      end if


!     ******************************
!     ***   WRITE COEFFICIENTS   ***
!     ******************************

      if (myproc.eq.0) then

!        Headers to screen.

         if (veryver) then

            write(*,*)
            write(*,"(A20)") ' Shape coefficients:'

            if (offset) then
               write(*,*)
               write(*,"(A8,ES14.6)") '   xc  =',xc
               write(*,"(A8,ES14.6)") '   yc  =',yc
               write(*,"(A8,ES14.6)") '   zc  =',zc
            end if

            write(*,*)
            write(*,*) '  l            c0_l'
            write(*,*)
            write(*,"(I4,A6,ES14.6)") 0,'      ',c0(0)

         end if

!        Headers to logfile.

         if (logfile) then

            write(12,*)
            write(12,"(A20)") ' Shape coefficients:'

            if (offset) then
               write(12,*)
               write(12,"(A8,ES14.6)") '   xc  =',xc
               write(12,"(A8,ES14.6)") '   yc  =',yc
               write(12,"(A8,ES14.6)") '   zc  =',zc
            end if

            write(12,*)
            write(12,*) '  l            c0_l'
            write(12,*)
            write(12,"(I4,A6,ES14.6)") 0,'      ',c0(0)

         end if

!        c0 coefficients.

         do l=1+stepz,lmax,1+stepz
            if (veryver) then
               write(*,"(I4,A6,ES14.6)") l,'      ',c0(l)
            end if
            if (logfile) then
               write(12,"(I4,A6,ES14.6)") l,'      ',c0(l)
            end if
         end do

!        cc and cs coefficients (only for nonaxisymmetric case).

         if (nonaxi) then

!           Reflection symmetry on y:  some coefficients are not there.

            if (.not.refy) then

               if (veryver) then
                  write(*,*)
                  write(*,*) '  l   m       cc_lm          cs_lm'
                  write(*,*)
               end if

               if (logfile) then
                  write(12,*)
                  write(12,*) '  l   m       cc_lm          cs_lm'
                  write(12,*)
               end if

               do l=1,lmax
                  do m=1,l
                     if (stepz*mod(l-m,2).eq.0) then
                        if (veryver) then
                           write(*,10) l,m,'  ',cc(l,m),' ',cs(l,m)
                        end if
                        if (logfile) then
                           write(12,10) l,m,'  ',cc(l,m),' ',cs(l,m)
                        end if
   10                   format(I4,I4,A2,ES14.6,A1,ES14.6)
                     end if
                  end do
               end do

!           No reflection symmetry on y.

            else

               if (veryver) then
                  write(*,*)
                  write(*,*) '  l   m       cc_lm'
                  write(*,*)
               end if

               if (logfile) then
                  write(12,*)
                  write(12,*) '  l   m       cc_lm'
                  write(12,*)
               end if

               do l=1,lmax
                  do m=1+stepx,l,1+stepx
                     if (stepz*mod(l-m,2).eq.0) then
                        if (veryver) then
                           write(*,20) l,m,'  ',cc(l,m)
                        end if
                        if (logfile) then
                           write(12,20) l,m,'  ',cc(l,m)
                        end if
   20                   format(I4,I4,A2,ES14.6)
                     end if
                  end do
               end do

            end if
         end if
      end if


!     *********************
!     ***   FIND FUNC   ***
!     *********************

!     The value of FUNC will be equal to the surface integral
!     of the square of the expansion or to the surface area,
!     depending on what we want to minimize.

      if (minarea) then
         FUNC = intarea
      else
         if(find_trapped_surface) then
            FUNC = intexpdel2
         else
            FUNC = intexp2
         endif
      end if


!     *************************
!     ***   CLOSE LOGFILE   ***
!     *************************

      if (logfile.and.(myproc.eq.0)) then
         close(12)
      end if


!     ***************
!     ***   END   ***
!     ***************

      end function FUNC