Many small changes to fix some bugs in the 3 horizon mode that affected

octant symmetry.  This changes have no effect in the 1 horizon mode, nor
do thye have any effect in 3 horizon mode with a full grid.


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

1 files changed, 104 insertions, 103 deletions

diff --git a/src/AHFinder.F b/src/AHFinder.F
index e6ceb06..0f38823 100644
--- a/src/AHFinder.F
+++ b/src/AHFinder.F
@@ -36,7 +36,7 @@
       integer istat
 
       CCTK_REAL atime
-      CCTK_REAL xmn,ymn,zmn,xmx,ymx,zmx,rmx,r0
+      CCTK_REAL 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
@@ -67,15 +67,6 @@
 !
 !     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.
@@ -129,6 +120,17 @@
       if (mod(int(cctk_iteration)-ahfafter,ahfso).ne.0) return
 
 
+!     *******************************************************
+!     ***   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
+
+
 !     ***************************
 !     ***   GRID PARAMETERS   ***
 !     ***************************
@@ -149,13 +151,6 @@
       call AHFinder_SetSym(CCTK_FARGUMENTS)
 
 
-!     ********************************
-!     ***   INITIALIZE  {status}   ***
-!     ********************************
-
-      status = .false.
-
-
 !     **********************************************
 !     ***   FIND OUT ON WHICH PROCESSOR WE ARE   ***
 !     **********************************************
@@ -163,12 +158,45 @@
       myproc = CCTK_MyProc(cctkGH)
 
 
-!     *********************************
-!     ***   INITIALIZE FIRSTCALLS   ***
-!     *********************************
+!     *******************************************
+!     ***   FIND  {xmn,xmx,ymn,ymx,zmn,zmx}   ***
+!     *******************************************
+    
+      call CCTK_CoordRange(ierr,cctkGH,xmn,xmx,-1,"x","cart3d")
+      call CCTK_CoordRange(ierr,cctkGH,ymn,ymx,-1,"y","cart3d")
+      call CCTK_CoordRange(ierr,cctkGH,zmn,zmx,-1,"z","cart3d")
 
-      firstleg = .true.
-      firstint = .true.
+
+!     **************************************************
+!     ***   TRANSFORM CONFORMAL TO PHYSICAL METRIC   ***
+!     **************************************************
+
+      ahfconformal = .false.
+
+!     The horizon finding routines know nothing about conformal
+!     metrics, so we need to make sure that we are working with
+!     the physical metric.  I change not only the metric, but
+!     also the "conformal state" flag to make sure that any
+!     macros I use have always the physical metric.  However,
+!     I need to remember the original value of the flag to change
+!     it back at the end if necesary.
+
+      if (conformal_state == CONFORMAL_METRIC) then
+
+         call CCTK_WARN(4,"Converting metric: conformal -> physical")
+
+         ahfconformal = .true.
+
+         conformal_state = NOCONFORMAL_METRIC
+
+         gxx = gxx*psi**4
+         gxy = gxy*psi**4
+         gxz = gxz*psi**4
+         gyy = gyy*psi**4
+         gyz = gyz*psi**4
+         gzz = gzz*psi**4
+
+      end if
 
 
 !     ***********************************************
@@ -186,6 +214,14 @@
 !     ***   ADD TO COUNT OF CALLS TO FINDER   ***
 !     *******************************************
 
+!     IMPORTANT:  When looking for 3 horizons, the label 1 below
+!     marks the point where one comes back in order to look for
+!     the next horizon.
+
+!     Because of this, anything that is likely to be different
+!     for the different horizons should NOT be initialized until
+!     after this point!
+
  1    continue
 
 !     Say what we are doing.
@@ -222,20 +258,24 @@
       end if
 
 
-!     *******************************************************
-!     ***   FIND OUT IF WE WANT TO FIND TRAPPED SURFACE   ***
-!     *******************************************************
+!     ********************************
+!     ***   INITIALIZE  {status}   ***
+!     ********************************
+
+      status = .false.
 
-      if (ahf_trapped_surface.ne.0) then
-         find_trapped_surface = .true.
-      else
-         find_trapped_surface = .false.
-      end if
 
+!     *********************************
+!     ***   INITIALIZE FIRSTCALLS   ***
+!     *********************************
 
-!     ********************************************
-!     ***   GET THE NAME OF OUTPUT DIRECTORY   ***
-!     ********************************************
+      firstleg = .true.
+      firstint = .true.
+
+
+!     **********************
+!     ***   FILE NAMES   ***
+!     **********************
 
       call CCTK_FortranString(nfile,outdir,filestr)
 
@@ -298,38 +338,6 @@
       end if
 
 
-!     **************************************************
-!     ***   TRANSFORM CONFORMAL TO PHYSICAL METRIC   ***
-!     **************************************************
-
-      ahfconformal = .false.
-
-!     The horizon finding routines know nothing about conformal
-!     metrics, so we need to make sure that we are working with
-!     the physical metric.  I change not only the metric, but
-!     also the "conformal state" flag to make sure that any
-!     macros I use have always the physical metric.  However,
-!     I need to remember the original value of the flag to change
-!     it back at the end if necesary.
-
-      if (conformal_state == CONFORMAL_METRIC) then
-
-         call CCTK_WARN(4,"Converting metric: conformal -> physical")
-
-         ahfconformal = .true.
-
-         conformal_state = NOCONFORMAL_METRIC
-
-         gxx = gxx*psi**4
-         gxy = gxy*psi**4
-         gxz = gxz*psi**4
-         gyy = gyy*psi**4
-         gyz = gyz*psi**4
-         gzz = gzz*psi**4
-
-      end if
-
-
 !     ***************************************
 !     ***   FIND OUT THE TYPE OF OUTPUT   ***
 !     ***************************************
@@ -375,15 +383,6 @@
       end if
 
 
-!     *******************************************
-!     ***   FIND  {xmn,xmx,ymn,ymx,zmn,zmx}   ***
-!     *******************************************
-    
-      call CCTK_CoordRange(ierr,cctkGH,xmn,xmx,-1,"x","cart3d")
-      call CCTK_CoordRange(ierr,cctkGH,ymn,ymx,-1,"y","cart3d")
-      call CCTK_CoordRange(ierr,cctkGH,zmn,zmx,-1,"z","cart3d")
-
-
 !     **********************************
 !     ***   INITIALIZE  {xc,yc,zc}   ***
 !     **********************************
@@ -410,23 +409,20 @@
 
       if ((xc.gt.xmx).or.(xc.lt.xmn)) then
          write(*,*)
-         write(*,*) 'xc is out of the grid, reseting it to zero ...'
-         write(*,*)
-         xc = zero
+         write(*,*) 'xc is out of the grid, giving up ...'
+         goto 5
       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
+         write(*,*) 'yc is out of the grid, giving up ...'
+         goto 5
       end if
 
-      if ((yc.gt.ymx).or.(yc.lt.ymn)) then
+      if ((zc.gt.zmx).or.(zc.lt.zmn)) then
          write(*,*)
-         write(*,*) 'zc is out of the grid, reseting it to zero ...'
-         write(*,*)
-         zc = zero
+         write(*,*) 'zc is out of the grid, giving up ...'
+         goto 5
       end if
 
 
@@ -449,7 +445,7 @@
       end if
 
       if ((xc.ne.zero).or.(yc.ne.zero).or.(zc.ne.zero).or.
-     .     (wander)) offset = .true.
+     .    (wander)) offset = .true.
 
 !     Find out if we are in non-axisymmetric mode.
 
@@ -528,8 +524,6 @@
 #ifdef THORN_cartoon_2d
       if (cartoon_active.ne.0) then
          cartoon = .true.
-      else
-         cartoon = .false.
       end if
 #endif
 
@@ -581,9 +575,11 @@
          allocate(cc(lmax,lmax),cc_old(lmax,lmax))
          allocate(cs(lmax,lmax),cs_old(lmax,lmax))
 
-         cs = zero
+         c0 = zero
          cc = zero
+         cs = zero
 
+         c0_old = zero
          cc_old = zero
          cs_old = zero
 
@@ -747,7 +743,6 @@
          else
             write (*,*)
             write (*,*) 'ahf_r0 is too large, rescaling ...'
-            write (*,*)
          end if
       end if
 
@@ -849,6 +844,10 @@
             cc = zero
             cs = zero
 
+            c0_old = zero
+            cc_old = zero
+            cs_old = zero
+
             if ((myproc.eq.0).and.verbose) then
                write(*,*)
                write(*,"(A35,ES11.3)")  ' Initial guess sphere with radius: ',rmx
@@ -892,6 +891,8 @@
 !     ***   WRITE MESSAGES AND LOG FILE   ***
 !     ***************************************
 
+ 5    continue
+
       status_old = status
 
       if (status) then
@@ -1599,6 +1600,20 @@
       end if
 
 
+!     **********************************************
+!     ***   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
+
+
 !     **************************************************
 !     ***   TRANSFORM PHYSICAL TO CONFORMAL METRIC   ***
 !     **************************************************
@@ -1622,20 +1637,6 @@
       end if
 
 
-!     **********************************************
-!     ***   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   ***
 !     ***************