diff options
| -rw-r--r-- | interface.ccl | 38 | ||||
| -rw-r--r-- | param.ccl | 5 | ||||
| -rw-r--r-- | src/AHFinder_gau.F | 91 |
3 files changed, 92 insertions, 42 deletions
diff --git a/interface.ccl b/interface.ccl index 3dbac6a..6eb1977 100644 --- a/interface.ccl +++ b/interface.ccl @@ -88,24 +88,46 @@ ahfgrid3, ahf_exp3 } "Grid functions to use in find3 algorithm" +real drift_correct_vars type=SCALAR +{ + ahf_centroid_x, + ahf_centroid_y, + ahf_omega_corr, + ahf_omega_cum +} "Variables for drift correction" + + +####################################### +### SCALARS AND ARRAYS FOR OUTPUT ### +####################################### + REAL out_scalars TYPE=SCALAR { - out_mass, out_radius, out_meridian1, out_meridian2, - out_area, out_perimeter, out_asymx, out_asymy, out_asymz, - out_centerx, out_centery, out_centerz -} "AHFinder scalar variables to output" + out_mass, + out_radius, + out_area, + out_perimeter, + out_meridian1, + out_meridian2, + out_asymx, + out_asymy, + out_asymz, + out_centerx, + out_centery, + out_centerz +} "Output of scalar variables" REAL out_1d_legen TYPE=ARRAY DISTRIB=CONSTANT DIM=1 SIZE=ahf_lmax { out_c0 -} "Coefficients for LEGEN(l,0,)" +} "Output of c0 coefficients" REAL out_2d_legen TYPE=ARRAY DISTRIB=CONSTANT DIM=2 SIZE=ahf_lmax,ahf_lmax { - out_cc, out_cs -} "Coefficients for LEGEN(l,m,)" + out_cc,out_cs +} "Output of cc and cs coefficients" #REAL out_2d_gauss TYPE=ARRAY DISTRIB=CONSTANT DIM=2 SIZE=ahf_ntheta+1,ahf_nphi+1 #{ # out_gauss -#} "Gaussian Curvature" +#} "Output of gaussian curvature" @@ -84,7 +84,7 @@ BOOLEAN ahf_wander "Allow the center to wander?" { } "no" -INT ahf_lmax "Maximum number of terms in theta expansion" # STEERABLE = ALWAYS +INT ahf_lmax "Maximum number of terms in theta expansion" { 0:20 :: "Range from 0 to 20" } 8 @@ -511,5 +511,8 @@ EXTENDS KEYWORD shift "Which shift condition to use" "ahfinder" :: "Set the shift after finding horizon" } "" +USES REAL rotation_omega + + diff --git a/src/AHFinder_gau.F b/src/AHFinder_gau.F index 37be435..a345717 100644 --- a/src/AHFinder_gau.F +++ b/src/AHFinder_gau.F @@ -574,62 +574,87 @@ end do end do - avgx = avgx/(dble(nphi+1)*(ntheta+1)) - avgy = avgy/(dble(nphi+1)*(ntheta+1)) + avgx = avgx/dble(npoints) + avgy = avgy/dble(npoints) + + ahf_centroid_x = avgx + ahf_centroid_y = avgy + +! Calculate correction in omega to stop "drift" of apparent horizon. -! Alter shift to correct for rotational "drift" of apparent horizon. - if (drift_correct_on.eq.1) then - write(0,*) '**** AHFinder_gau: avgx = ',avgx - write(0,*) '**** AHFinder_gau: avgy = ',avgy + write(0,*) 'AHFinder_gau: ahf_centroid_x = ',ahf_centroid_x + write(0,*) 'AHFinder_gau: ahf_centroid_y = ',ahf_centroid_y xi_new = atan2(avgx,avgy) - if (ahf_ncall .gt. 1) then + if (ahf_ncall.gt.1) then + xi_drift = xi_new - xi_old - else - xi_drift = xi_new - endif - omega_corr = - xi_drift / (ahf_findevery * cctk_delta_time) + 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 - if (ahf_ncall .gt. 1) then - omega_cum = omega_cum + omega_corr else - omega_cum = omega_corr - endif - 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: omega_cum = ',omega_cum + omega_cum = rotation_omega - 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 + end if + + write(0,*) '**** AHFinder_gau: omega_cum = ',omega_cum xi_old = xi_new end if +! Other processors. + else xa = half*(xmx+xmn) ya = half*(ymx+ymn) za = half*(zmx+zmn) +! We set omega_corr to zero of 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 = 0.0D0 + + end if + +! Actually apply drift correction to shift. + + if (drift_correct_on.eq.1) then + + call CCTK_ReduceLocalScalar(omega_corr,cctkGH,-1,sum_handle, + . error1,rerror,CCTK_VARIABLE_INT) + + 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 + end if ! Reduce the errors across processors (all processors must |