diff options
author | shawley <shawley@89daf98e-ef62-4674-b946-b8ff9de2216c> | 2001-12-18 18:16:06 +0000 |
---|---|---|
committer | shawley <shawley@89daf98e-ef62-4674-b946-b8ff9de2216c> | 2001-12-18 18:16:06 +0000 |
commit | 4c0d48910d68993a129b73d40a862f263cae34a1 (patch) | |
tree | 47f6b14159cd2762156262fbad766c71d58d80b1 /src/AHFinder_gau.F | |
parent | 4ae20305102bd9ad0cb706279a7360e127cc2142 (diff) |
added radial shift drift-correction, changed 1.1 to cactus parameter
git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinder/trunk@266 89daf98e-ef62-4674-b946-b8ff9de2216c
Diffstat (limited to 'src/AHFinder_gau.F')
-rw-r--r-- | src/AHFinder_gau.F | 70 |
1 files changed, 58 insertions, 12 deletions
diff --git a/src/AHFinder_gau.F b/src/AHFinder_gau.F index 803c7da..063dd22 100644 --- a/src/AHFinder_gau.F +++ b/src/AHFinder_gau.F @@ -46,7 +46,7 @@ CCTK_REAL trxi,xi2 CCTK_REAL zero,half,one,two,three,four,pi CCTK_REAL aux,sina,cosa - CCTK_REAL omega_corr_all + CCTK_REAL omega_corr_all, radvel_corr_all CCTK_REAL, dimension(3,3) :: ug,xi CCTK_REAL, dimension(2,2) :: ga,ua @@ -583,18 +583,17 @@ 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 (drift_correct_on.eq.1) then - - if ((.not.find3).or. + if ((.not.find3).or. . (find3.and.(mfind-1.eq.drift_correct_horizon))) then + if (drift_correct_on.eq.1) then + 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 @@ -603,7 +602,7 @@ write(0,*) 'AHFinder_gau: xi_drift = ',xi_drift - omega_corr = - 1.1D0 * xi_drift + omega_corr = - drift_correct_factor * xi_drift . /(dble(ahf_findevery)*cctk_delta_time) omega_cum = omega_cum + omega_corr @@ -622,12 +621,31 @@ end if - write(0,*) '**** AHFinder_gau: omega_cum = ',omega_cum + write(0,*) 'AHFinder_gau: omega_cum = ',omega_cum xi_old = xi_new end if - 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 + + 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 ! Other processors. @@ -643,14 +661,15 @@ ! (Recommended by Thomas Radke) omega_corr = zero + radvel_corr = zero end if ! Actually apply drift correction to shift. - if (drift_correct_on.eq.1) then - if ((.not.find3).or. + 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) @@ -671,8 +690,35 @@ end if end if + + if (rad_drift_correct.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 + + 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 + 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). |