Added support for attenuated co-rotating shift in drift correction.

Fixed problem with ahf_omega_cum being reset to initial rotation parameter
after restart from a checkpoint.


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

1 files changed, 149 insertions, 30 deletions

diff --git a/src/AHFinder_gau.F b/src/AHFinder_gau.F
index 021d612..9a5b58e 100644
--- a/src/AHFinder_gau.F
+++ b/src/AHFinder_gau.F
@@ -25,6 +25,7 @@
 
       DECLARE_CCTK_ARGUMENTS
       DECLARE_CCTK_PARAMETERS
+      DECLARE_CCTK_FUNCTIONS
 
       logical firstgau
       logical firstcal(4)
@@ -60,7 +61,17 @@
       CCTK_REAL, allocatable, dimension(:,:) :: txx,tyy,tzz,txy,txz,tyz
       CCTK_REAL, allocatable, dimension(:,:) :: g11,g22,g12
 
-      character*200 gaussf
+      CCTK_REAL :: m1p, m2p, x1p, x2p,  y1p, y2p, z1p, z2p
+      CCTK_REAL :: r12, r22, sigma
+      CCTK_REAL, dimension(nx,ny,nz) :: fac
+      CCTK_INT :: use_att, apower, nchars
+      logical :: use_rot_att
+      logical :: str_comp
+      external :: str_comp
+
+      character(len=32) :: tmpstr
+
+      character(len=200) :: gaussf
 
 !     Declarations for macros.
 
@@ -189,6 +200,47 @@
          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.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   ***
@@ -581,28 +633,31 @@
          avgx = avgx/dble(npoints)
          avgy = avgy/dble(npoints)
 
-         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
-
 !        Calculate correction in omega to stop "drift" of apparent horizon.
+
          if ((.not.find3).or.
-     .          (find3.and.(mfind-1.eq.drift_correct_horizon))) then
+     .       (find3.and.(mfind-1.eq.drift_correct_horizon))) then
 
             if (drift_correct_on.eq.1) then
 
+               ahf_centroid_x = avgx
+               ahf_centroid_y = 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
 
-               if (ahf_ncall.gt.1) then
+               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 = - 1.1 * xi_drift
+                  omega_corr = - drift_correct_factor * xi_drift
      .               /(dble(ahf_findevery)*cctk_delta_time)
                   omega_cum  = omega_cum + omega_corr
 
@@ -614,11 +669,24 @@
      .             cctk_delta_time
                   write(0,*) 'AHFinder_gau: omega_corr = ',omega_corr
 
+                  ahf_omega_corr = omega_corr
+                  ahf_omega_cum  = omega_cum
+
                else
 
-                   omega_cum = rotation_omega
+                   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
@@ -628,23 +696,26 @@
             end if
 
 !        Do drift correction in the radial direction
-           if (rad_drift_correct.eq.1) then
-                r_new = sqrt(avgx**2 + avgy**2)
 
-                if (ahf_ncall.gt.1) then
+            if (rad_drift_correct.eq.1) then
+
+               r_new = sqrt(avgx**2 + avgy**2)
 
-                   radvel_corr = - rad_drift_correct_factor * 
-     .                (r_new - r_old) /(dble(ahf_findevery)*cctk_delta_time)
+               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
+                  write(0,*) 'AHFinder_gau: radvel_corr = ',radvel_corr
+               else
+                  radvel_corr = 0
+               endif
 
-                r_old = r_new
+               r_old = r_new
 
             endif
-        endif
+         endif
+
 
 !     Other processors.
 
@@ -678,19 +749,35 @@
 
             omega_corr = omega_corr_all
 
-            ahf_omega_corr = omega_corr
-            ahf_omega_cum  = omega_cum
+!            ahf_omega_corr = omega_corr
+!            ahf_omega_cum  = omega_cum
 
             if (ahf_ncall.gt.1) then
 
-               betax = betax - y*omega_corr/psi**rotation_psipower
-               betay = betay + x*omega_corr/psi**rotation_psipower
+               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.eq.1) then
+         if (rad_drift_correct.eq.1) then
 
             call CCTK_ReduceLocalScalar(ierror,cctkGH,-1,sum_handle,
      .         radvel_corr,radvel_corr_all,CCTK_VARIABLE_REAL)
@@ -704,12 +791,27 @@
 
             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) /
-     .                     psi**rotation_psipower
-                  betay = betay + radvel_corr * y/sqrt(x**2+y**2) / 
-     .                     psi**rotation_psipower
+                  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
@@ -1329,3 +1431,20 @@
       end subroutine AHFinder_gau
 
 
+!     ******************************************
+!     *** String comparison utility function ***
+!     ******************************************
+
+      function str_comp(a1,a2)
+
+      logical :: str_comp
+      character(len=*),intent(in) :: a1, a2
+
+      if ( ( verify(trim(a1),trim(a2)).eq.0 )  .and.
+     .                     ( len_trim(a1) .eq. len_trim(a2) ) ) then
+        str_comp = .true.
+      else
+        str_comp = .false.
+      end if
+
+      end function str_comp