Adding parameter "drift_correct_horizon" to apply the drift only to the

chosen horizon.


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

2 files changed, 49 insertions, 43 deletions

diff --git a/param.ccl b/param.ccl
index d0bae38..7c6555c 100644
--- a/param.ccl
+++ b/param.ccl
@@ -64,10 +64,17 @@ REAL trapped_surface_delta "find (expansion = delta) surface" STEERABLE = ALWAYS
 ################################################
 ### Parameters for drift correction/tracking ###
 ################################################
+
 BOOLEAN drift_correct_on "Apply drift correction (horizon-tracking shift rotation)?"
 {
 } "no"
 
+INT drift_correct_horizon "Which horizon to use for drift correction"
+{
+0:2  :: "There are only three horizons"
+} 0
+
+
 ########################################
 ### Parameters for surface expansion ###
 ########################################
diff --git a/src/AHFinder_gau.F b/src/AHFinder_gau.F
index 3c1f75c..c9cbf4b 100644
--- a/src/AHFinder_gau.F
+++ b/src/AHFinder_gau.F
@@ -592,41 +592,44 @@
 !        Calculate correction in omega to stop "drift" of apparent horizon.
 
          if (drift_correct_on.eq.1) then
+            if ((.not.find3).or.
+     .          (find3.and.(mfind.eq.drift_correct_horizon))) then
 
-           xi_new = atan2(avgx,avgy)
+               xi_new = atan2(avgx,avgy)
 
-           write(0,*) 'AHFinder_gau: ahf_centroid_x = ',ahf_centroid_x
-           write(0,*) 'AHFinder_gau: ahf_centroid_y = ',ahf_centroid_y
-           write(0,*) 'AHFinder_gau: xi_new = ',xi_new
+               write(0,*) 'AHFinder_gau: ahf_centroid_x = ',ahf_centroid_x
+               write(0,*) 'AHFinder_gau: ahf_centroid_y = ',ahf_centroid_y
+               write(0,*) 'AHFinder_gau: xi_new = ',xi_new
 
-           if (ahf_ncall.gt.1) then
+               if (ahf_ncall.gt.1) then
 
-              xi_drift = xi_new - xi_old
+                  xi_drift = xi_new - xi_old
 
-              write(0,*) 'AHFinder_gau: xi_drift = ',xi_drift
+                  write(0,*) 'AHFinder_gau: xi_drift = ',xi_drift
 
-              omega_corr = - 1.1D0 * xi_drift
-     .           /(dble(ahf_findevery)*cctk_delta_time)
-              omega_cum  = omega_cum + omega_corr
+                  omega_corr = - 1.1D0 * xi_drift
+     .               /(dble(ahf_findevery)*cctk_delta_time)
+                  omega_cum  = omega_cum + omega_corr
 
-              write(0,*) 'AHFinder_gau: ncall = ',ahf_ncall
-              write(0,*) 'AHFinder_gau: xi_drift = ',xi_drift
-              write(0,*) 'AHFinder_gau: ahf_findevery = ',
-     .         ahf_findevery
-              write(0,*) 'AHFinder_gau: cctk_delta_time = ',
-     .         cctk_delta_time
-              write(0,*) 'AHFinder_gau: omega_corr = ',omega_corr
+                  write(0,*) 'AHFinder_gau: ncall = ',ahf_ncall
+                  write(0,*) 'AHFinder_gau: xi_drift = ',xi_drift
+                  write(0,*) 'AHFinder_gau: ahf_findevery = ',
+     .             ahf_findevery
+                  write(0,*) 'AHFinder_gau: cctk_delta_time = ',
+     .             cctk_delta_time
+                  write(0,*) 'AHFinder_gau: omega_corr = ',omega_corr
 
-           else
+               else
 
-               omega_cum = rotation_omega
+                   omega_cum = rotation_omega
 
-           end if
+               end if
 
-           write(0,*) '**** AHFinder_gau: omega_cum = ',omega_cum
+               write(0,*) '**** AHFinder_gau: omega_cum = ',omega_cum
 
-           xi_old = xi_new
+               xi_old = xi_new
 
+            end if
          end if
 
 !     Other processors.
@@ -637,7 +640,7 @@
          ya = half*(ymx+ymn)
          za = half*(zmx+zmn)
 
-!        We set omega_corr to zero of myproc.ne.0, and then sum omega_corr
+!        We set omega_corr to zero for (myproc.ne.0), and then sum omega_corr
 !        over all processors.  This is a cheap way of communicating
 !        omega_corr from processor 0 to all other processors.
 !        (Recommended by Thomas Radke)
@@ -649,28 +652,24 @@
 !     Actually apply drift correction to shift.
 
       if (drift_correct_on.eq.1) then
+         if ((.not.find3).or.
+     .       (find3.and.(mfind.eq.drift_correct_horizon))) then
 
-          call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,sum_handle,
-     .       omega_corr,omega_corr_all,CCTK_VARIABLE_REAL)
-          if (ierror.ne.0) then
-             call CCTK_WARN(1,"Reduction failed!")
-          end if
-          omega_corr = omega_corr_all
-
-          ahf_omega_corr = omega_corr
-          ahf_omega_cum  = omega_cum
-
-          do k=1,nz
-            do j=1,ny
-              do i=1,nx
-                betax(i,j,k) = betax(i,j,k) 
-     .             - y(i,j,k)*omega_corr/psi(i,j,k)**2
-                betay(i,j,k) = betay(i,j,k) 
-     .             + x(i,j,k)*omega_corr/psi(i,j,k)**2
-             end do
-            end do
-          end do
+            call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,sum_handle,
+     .         omega_corr,omega_corr_all,CCTK_VARIABLE_REAL)
+            if (ierror.ne.0) then
+               call CCTK_WARN(1,"Reduction failed!")
+            end if
 
+            omega_corr = omega_corr_all
+
+            ahf_omega_corr = omega_corr
+            ahf_omega_cum  = omega_cum
+
+            betax = betax - y*omega_corr/psi**2
+            betay = betay + x*omega_corr/psi**2
+
+         end if
       end if
 
 !     Reduce the errors across processors (all processors must