Moved drift correction code to AEIDevelopment/DriftCorrect.

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

1 files changed, 6 insertions, 217 deletions

diff --git a/src/AHFinder_gau.F b/src/AHFinder_gau.F
index 2978e71..4baad67 100644
--- a/src/AHFinder_gau.F
+++ b/src/AHFinder_gau.F
@@ -47,7 +47,6 @@
       CCTK_REAL trxi,xi2
       CCTK_REAL zero,half,one,two,three,four,pi
       CCTK_REAL aux,sina,cosa
-      CCTK_REAL omega_corr_all, radvel_corr_all
 
       CCTK_REAL, dimension(3,3)     :: ug,xi
       CCTK_REAL, dimension(2,2)     :: ga,ua
@@ -200,48 +199,6 @@
          npoints = 1
       end if
 
-!     *****************************************************************
-!     ***   Get parameters for drift-correction with attenutation   ***
-!     ***   function for puncture data if necessary.                ***
-!     *****************************************************************
-
-      if ((drift_correct_on.eq.1).or.(rad_drift_correct_on.eq.1)) then
-        use_att = 0
-        call CCTK_ParameterValString(nchars,
-     .           "rot_shift_att","Einstein",tmpstr)
-        if ( str_comp(tmpstr,"yes") .or. str_comp(tmpstr,"1") ) then
-          call CCTK_ParameterValString(nchars,
-     .             "initial_data","Einstein",tmpstr)
-          if ( str_comp(tmpstr,"BrBr") ) then
-            use_rot_att = .true.
-            call CCTK_ParameterValString(nchars,"bhm1","BAM_Elliptic",tmpstr)
-            read (tmpstr,*) m1p
-            call CCTK_ParameterValString(nchars,"bhx1","BAM_Elliptic",tmpstr)
-            read (tmpstr,*) x1p
-            call CCTK_ParameterValString(nchars,"bhy1","BAM_Elliptic",tmpstr)
-            read (tmpstr,*) y1p
-            call CCTK_ParameterValString(nchars,"bhz1","BAM_Elliptic",tmpstr)
-            read (tmpstr,*) z1p
-            call CCTK_ParameterValString(nchars,"bhm2","BAM_Elliptic",tmpstr)
-            read (tmpstr,*) m2p
-            call CCTK_ParameterValString(nchars,"bhx2","BAM_Elliptic",tmpstr)
-            read (tmpstr,*) x2p
-            call CCTK_ParameterValString(nchars,"bhy2","BAM_Elliptic",tmpstr)
-            read (tmpstr,*) y2p
-            call CCTK_ParameterValString(nchars,"bhz2","BAM_Elliptic",tmpstr)
-            read (tmpstr,*) z2p
-            call CCTK_ParameterValString(nchars,"rot_shift_att_pow",
-     .                                               "Einstein",tmpstr)
-            read (tmpstr,*) apower
-            call CCTK_ParameterValString(nchars,"rot_shift_att_sigma",
-     .                                               "Einstein",tmpstr)
-            read (tmpstr,*) sigma
-            if ( m1p .gt. zero ) use_att = use_att + 1
-            if ( m2p .gt. zero ) use_att = use_att + 1
-          end if
-        end if
-      end if
-
 !     ******************************************************
 !     ***   Get the reduction handle for sum operation   ***
 !     ******************************************************
@@ -633,89 +590,13 @@
          avgx = avgx/dble(npoints)
          avgy = avgy/dble(npoints)
 
-!        Calculate correction in omega to stop "drift" of apparent horizon.
-
-         if ((.not.find3).or.
-     .       (find3.and.(mfind-1.eq.drift_correct_horizon))) then
-
-             ahf_centroid_x = avgx
-             ahf_centroid_y = avgy
-             write(0,*) 'AHFinder_gau: ahf_centroid_x = ',ahf_centroid_x
-             write(0,*) 'AHFinder_gau: ahf_centroid_y = ',ahf_centroid_y
-
-            if (drift_correct_on.eq.1) then
-
-
-               xi_new = atan2(avgx,avgy)
-
-               write(0,*) 'AHFinder_gau: xi_new = ',xi_new
-
-               write(0,*) 'AHFinder_gau: ahf_omega_cum = ',ahf_omega_cum
-
-               if (not_dc_first) then
-
-                  xi_drift = xi_new - xi_old
-
-                  write(0,*) 'AHFinder_gau: xi_drift = ',xi_drift
-
-                  omega_corr = - drift_correct_factor * 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
-
-                  ahf_omega_corr = omega_corr
-                  ahf_omega_cum  = omega_cum
-
-               else
-
-                   if ( drift_first .eq. 1 ) then
-                      ahf_omega_cum = rotation_omega
-                      omega_cum = rotation_omega
-                      drift_first = 0
-                   else
-                      omega_cum = ahf_omega_cum
-                   end if
-                   omega_corr = zero
-
-                   not_dc_first = .true.
-
-                   write(0,*) 'AHFinder_gau: ahf_omega_cum = ',ahf_omega_cum
-
-               end if
-
-               write(0,*) 'AHFinder_gau: omega_cum = ',omega_cum
-
-               xi_old = xi_new
-
-            end if
-
-!        Do drift correction in the radial direction
-           if (rad_drift_correct_on.eq.1) then
-
-               r_new = sqrt(avgx**2 + avgy**2)
-
-               if (ahf_ncall.gt.1) then
-
-                  radvel_corr = - rad_drift_correct_factor * 
-     .               (r_new - r_old) /(dble(ahf_findevery)*cctk_delta_time)
-
-                  write(0,*) 'AHFinder_gau: radvel_corr = ',radvel_corr
-               else
-                  radvel_corr = 0
-               endif
-
-               r_old = r_new
-
-            endif
-         endif
+         ahf_centroid_x = avgx
+         ahf_centroid_y = avgy
 
+         if (verbose) then
+           write(1,*) 'AHFinder_gau: ahf_centroid_x = ', ahf_centroid_x
+           write(1,*) 'AHFinder_gau: ahf_centroid_y = ', ahf_centroid_y
+         end if
 
 !     Other processors.
 
@@ -725,101 +606,9 @@
          ya = half*(ymx+ymn)
          za = half*(zmx+zmn)
 
-!        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)
-
-         omega_corr = zero
-         radvel_corr = zero
-
-      end if
-
-!     Actually apply drift correction to shift.
-
-      if ((.not.find3).or.
-     .       (find3.and.(mfind-1.eq.drift_correct_horizon))) then
-        if (drift_correct_on.eq.1) 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
-
-            if (ahf_ncall.gt.1) then
-
-               if ( use_att .eq. 0 ) then
-                  fac = one / psi**rotation_psipower
-               else
-                  fac = one
-                  if ( m1p .gt. zero ) then
-                     fac = fac - exp( - (
-     .                       ((x-x1p)**2+(y-y1p)**2+(z-z1p)**2)/sigma**2
-     .                                                            )**apower )
-                  end if
-                  if ( m2p .gt. zero ) then
-                     fac = fac - exp( - (
-     .                       ((x-x2p)**2+(y-y2p)**2+(z-z2p)**2)/sigma**2
-     .                                                            )**apower )
-                  end if
-               end if
-
-               betax = betax - y*omega_corr*fac
-               betay = betay + x*omega_corr*fac
-
-            end if
-
-         end if
-
-        if (rad_drift_correct_on.eq.1) then
-
-            call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,sum_handle,
-     .         radvel_corr,radvel_corr_all,CCTK_VARIABLE_REAL)
-            if (ierror.ne.0) then
-               call CCTK_WARN(1,"Reduction failed!")
-            end if
-
-            radvel_corr = radvel_corr_all
-
-            ahf_radvel_corr = radvel_corr
-
-            if (ahf_ncall.gt.1) then
-               
-               if ( use_att .eq. 0 ) then
-                  fac = one / psi**rotation_psipower
-               else
-                  fac = one
-                  if ( m1p .gt. zero ) then
-                     fac = fac - exp( - (
-     .                       ((x-x1p)**2+(y-y1p)**2+(z-z1p)**2)/sigma**2
-     .                                                            )**apower )
-                  end if
-                  if ( m2p .gt. zero ) then
-                     fac = fac - exp( - (
-     .                       ((x-x2p)**2+(y-y2p)**2+(z-z2p)**2)/sigma**2
-     .                                                            )**apower )
-                  end if
-               end if
-
-               where (x**2+y**2 .gt. 0.001)
-                  betax = betax + radvel_corr * x/sqrt(x**2+y**2) * fac
-                  betay = betay + radvel_corr * y/sqrt(x**2+y**2) * fac
-               end where
-
-            end if
-
-         end if
       end if
 
 
-
-
 !     Reduce the errors across processors (all processors must
 !     know about this since all will participate on the interpolation
 !     below).