Making sure the flow algorithm knows about different initial guess options.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinder/trunk@142 89daf98e-ef62-4674-b946-b8ff9de2216c

1 files changed, 65 insertions, 24 deletions

diff --git a/src/AHFinder.F b/src/AHFinder.F
index c270b76..178ab05 100644
--- a/src/AHFinder.F
+++ b/src/AHFinder.F
@@ -40,7 +40,7 @@
       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 c0_temp(0:18)
       CCTK_REAL maxval(3),minval(3)
       CCTK_REAL rhor
 
@@ -714,9 +714,9 @@
       end if
 
 
-!     *************************************************
-!     ***   GET INITIAL GUESS AND LOOK FOR HORIZON  ***
-!     *************************************************
+!     *****************************
+!     ***   GET INITIAL GUESS   ***
+!     *****************************
 
 !     Initializing flags
 
@@ -752,26 +752,27 @@
          minarea = .false.
       end if
 
-!     Get initial guess and jump into finder routines
+!     Get initial guess.
 
-      if ((.not.flow).and.guessold.and.(ahf_ncall.gt.1)) then
+      if (guessold.and.status_old) then
 
-!        Minimization algorithm and old horizon as initial guess.
+!        Old horizon as initial guess.
 
-         call AHFinder_min(CCTK_FARGUMENTS,NN,status,logf)
+         if ((myproc.eq.0).and.verbose) then
+            write(*,*)
+            write(*,*)  'AHFinder: Using old horizon as initial guess'
+         end if
 
-      else if ((.not.flow).and.manual_guess) then
+      else if (manual_guess) then
 
-!        Minimization algorithm and manual guess.
+!        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
+         if ((myproc.eq.0).and.verbose) then
+            write(*,*)
+            write(*,*)  'AHFinder: Manual guess'
+         end if
 
-         c0(0) = ahf_l0_guess
+         c0_temp(0) = ahf_l0_guess
          c0_temp(2) = ahf_l2_guess
          c0_temp(4) = ahf_l4_guess
          c0_temp(6) = ahf_l6_guess
@@ -782,17 +783,55 @@
          c0_temp(16) = ahf_l16_guess
          c0_temp(18) = ahf_l18_guess
 
-         do i=2,lmax,2
-            c0(i) = c0_temp(i)
-         enddo
+         if (lmax.le.18) then
+            do i=0,lmax
+               c0(i) = c0_temp(i)
+            end do
+         else
+            do i=0,18
+               c0(i) = c0_temp(i)
+            end do
+         end if
 
-         call AHFinder_min(CCTK_FARGUMENTS,NN,status,logf)
+      else
+
+!        Standard initial guess.
 
-      else if (.not.flow) then
+         if (.not.flow) then
+
+!           For minimization, call initial guess routine.
+
+            call AHFinder_initguess(CCTK_FARGUMENTS,rmx)
+
+         else
 
-!        Straight minimization algorithm.
+!           For flow, take as initial guess a large sphere.
+
+            c0 = zero
+            cc = zero
+            cs = zero
+
+            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
+
+         end if
+
+      end if
+
+
+!     ***************************
+!     ***   LOOK FOR HORIZON  ***
+!     ***************************
+
+      if (.not.flow) then
+
+!        Minimization algorithm.
 
-         call AHFinder_initguess(CCTK_FARGUMENTS,rmx)
          call AHFinder_min(CCTK_FARGUMENTS,NN,status,logf)
 
       else
@@ -815,6 +854,8 @@
 !     ***   WRITE MESSAGES AND LOG FILE   ***
 !     ***************************************
 
+      status_old = status
+
       if (status) then
 
 !        Check if the expansion on the surface is small enough