From 9882f6095edcfe85002713856243046620f1df17 Mon Sep 17 00:00:00 2001 From: miguel Date: Tue, 8 Oct 2002 20:43:33 +0000 Subject: minor cleaning up, plus a stupid fix to a problem with the absoft compiler. It turns out I needed to change the order of two definitions to stop the damm compiler from inventing a non-existent NaN. &%$#*! git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinder/trunk@316 89daf98e-ef62-4674-b946-b8ff9de2216c --- src/AHFinder_dat.F | 10 +- src/AHFinder_exp.F | 311 ++++++++++++++++++++++++++++++----------------------- src/AHFinder_gau.F | 23 +++- 3 files changed, 202 insertions(+), 142 deletions(-) (limited to 'src') diff --git a/src/AHFinder_dat.F b/src/AHFinder_dat.F index 4e31ed8..a7fc589 100644 --- a/src/AHFinder_dat.F +++ b/src/AHFinder_dat.F @@ -49,6 +49,8 @@ CCTK_REAL hw,cw,nw CCTK_REAL dx,dy,dz + CCTK_REAL avgx,avgy,avgz + CCTK_REAL, allocatable, dimension (:) :: c0,c0_0,c0_1,c0_2,c0_old CCTK_REAL, allocatable, dimension (:,:) :: cc,cc_0,cc_1,cc_2,cc_old CCTK_REAL, allocatable, dimension (:,:) :: cs,cs_0,cs_1,cs_2,cs_old @@ -59,8 +61,6 @@ character*200 filestr - CCTK_REAL avgx, avgy - data ahf_ncall / 0 / data status_old / .false. / data status_old_0 / .false. / @@ -163,9 +163,11 @@ ! ! inside_min_neg_count Number of elements in integral that have ! negative expansion. +! +! avgx X-position of centroid of apparent horizon +! avgy Y-position of centroid of apparent horizon +! avgz Z-position of centroid of apparent horizon -! avgx X-position of centroid of apparent horizon -! avgy Y-position of centroid of apparent horizon ! *************** ! *** END *** diff --git a/src/AHFinder_exp.F b/src/AHFinder_exp.F index e990501..388b30f 100644 --- a/src/AHFinder_exp.F +++ b/src/AHFinder_exp.F @@ -48,7 +48,8 @@ CCTK_REAL c2fxx,c2fyy,c2fzz,c2fxy,c2fxz,c2fyz CCTK_REAL dfx,dfy,dfz,dfux,dfuy,dfuz CCTK_REAL aux,T0,T1,T2,T3,T4 - CCTK_REAL idx,idy,idz,idx2,idy2,idz2,idxy,idxz,idyz + CCTK_REAL idx,idy,idz + CCTK_REAL idx2,idy2,idz2,idxy,idxz,idyz CCTK_REAL zero,half,one,two ! Description of variables: @@ -105,14 +106,14 @@ idy = half/dy idz = half/dz - idxy = idx*idy - idxz = idx*idz - idyz = idy*idz - idx2 = one/dx**2 idy2 = one/dy**2 idz2 = one/dz**2 + idxy = idx*idy + idxz = idx*idz + idyz = idy*idz + do k=2,nz-1 do j=2,ny-1 do i=2,nx-1 @@ -163,11 +164,13 @@ det = gdxx*gdyy*gdzz + two*gdxy*gdxz*gdyz . - gdxx*gdyz**2 - gdyy*gdxz**2 - gdzz*gdxy**2 +! If determinant is not zero proceed. + if (det.gt.zero) then idet = one/det -! Find inverse spatial metric. +! Find inverse spatial metric. guxx = idet*(gdyy*gdzz - gdyz**2) guyy = idet*(gdxx*gdzz - gdxz**2) @@ -177,7 +180,7 @@ guxz = idet*(gdxy*gdyz - gdyy*gdxz) guyz = idet*(gdxy*gdxz - gdxx*gdyz) -! Find spatial derivatives of f. +! Find spatial derivatives of f. T0 = two*ahfgrid(i,j,k) @@ -199,13 +202,13 @@ dfz = (T1 - T2)*idz d2fzz = (T1 - T0 + T2)*idz2 -! Save gradient of horizon function and its norm -! (they will be needed later in the flow algorithm -! and in the gaussian curvature). Here I also try -! to avoid possible division by zero later by not -! allowing ahfgradn(i,j,k) = 0. This should only -! ever happen far from the horizon, so resetting this -! to 1 should have no important effects. +! Save gradient of horizon function and its norm +! (they will be needed later in the flow algorithm +! and in the gaussian curvature). Here I also try +! to avoid possible division by zero later by not +! allowing ahfgradn(i,j,k) = 0. This should only +! ever happen far from the horizon, so resetting this +! to 1 should have no important effects. ahfgradx(i,j,k) = dfx ahfgrady(i,j,k) = dfy @@ -217,22 +220,22 @@ ahfgradn(i,j,k) = sqrt(aux) if (ahfgradn(i,j,k).eq.zero) ahfgradn(i,j,k) = one -! Find crossed derivatives. +! Find crossed derivatives. d2fxy = (ahfgrid(i+1,j+1,k) + ahfgrid(i-1,j-1,k) - . - ahfgrid(i+1,j-1,k) - ahfgrid(i-1,j+1,k))*idxy + . - ahfgrid(i+1,j-1,k) - ahfgrid(i-1,j+1,k))*idxy d2fxz = (ahfgrid(i+1,j,k+1) + ahfgrid(i-1,j,k-1) - . - ahfgrid(i+1,j,k-1) - ahfgrid(i-1,j,k+1))*idxz + . - ahfgrid(i+1,j,k-1) - ahfgrid(i-1,j,k+1))*idxz d2fyz = (ahfgrid(i,j+1,k+1) + ahfgrid(i,j-1,k-1) - . - ahfgrid(i,j+1,k-1) - ahfgrid(i,j-1,k+1))*idyz - -! Raise indices in first derivatives. + . - ahfgrid(i,j+1,k-1) - ahfgrid(i,j-1,k+1))*idyz + +! Raise indices in first derivatives. dfux = guxx*dfx + guxy*dfy + guxz*dfz dfuy = guxy*dfx + guyy*dfy + guyz*dfz dfuz = guxz*dfx + guyz*dfy + guzz*dfz -! Find second covariant derivatives of f. +! Find second covariant derivatives of f. c2fxx = d2fxx - half*(dfux*ddxxx . + dfuy*(two*ddxxy - ddyxx) @@ -255,35 +258,38 @@ c2fyz = d2fyz - half*(dfuy*ddzyy + dfuz*ddyzz . + dfux*(ddyxz + ddzxy - ddxyz)) -! / m \ 1/2 -! Find: u = | d f d f | -! \ m / +! / m \ 1/2 +! Find: u = | d f d f | +! \ m / T0 = dfx*dfux + dfy*dfuy + dfz*dfuz +! if T0 is positive proceed. + if (T0.gt.zero) then + T0 = one/sqrt(T0) -! __2 -! Find: \/ f / u +! __2 +! Find: \/ f / u T1 = guxx*c2fxx + guyy*c2fyy + guzz*c2fzz . + two*(guxy*c2fxy + guxz*c2fxz + guyz*c2fyz) T1 = T1*T0 -! a b 2 -! Find: K d f d f / u -! ab +! a b 2 +! Find: K d f d f / u +! ab T2 = kdxx*dfux**2 + kdyy*dfuy**2 + kdzz*dfuz**2 . + two*(dfux*(kdxy*dfuy + kdxz*dfuz) + kdyz*dfuy*dfuz) T2 = T2*T0**2 -! __a__b 3 -! Find: \/ \/ f d f d f / u -! a b +! __a__b 3 +! Find: \/ \/ f d f d f / u +! a b T3 = c2fxx*dfux**2 + c2fyy*dfuy**2 + c2fzz*dfuz**2 . + two*(dfux*(c2fxy*dfuy + c2fxz*dfuz) @@ -291,132 +297,169 @@ T3 = T3*T0**3 -! Find: trK +! Find: trK T4 = guxx*kdxx + guyy*kdyy + guzz*kdzz . + two*(guxy*kdxy + guxz*kdxz + guyz*kdyz) -! Find the expansion. +! Find the expansion. ahf_exp(i,j,k) = T1 + T2 - T3 - T4 +! T0 was not positive. + else ahf_exp(i,j,k) = zero - endif + end if + +! Determinant was zero. + else ahf_exp(i,j,k) = zero - endif + end if + end do end do end do -! Zero out the edges before doing the boundaries - ahf_exp( 1, 1, :) = 0 - ahf_exp( 1,ny, :) = 0 - ahf_exp( 1, :, 1) = 0 - ahf_exp( 1, :,nz) = 0 - ahf_exp(nx, 1, :) = 0 - ahf_exp(nx,ny, :) = 0 - ahf_exp(nx, :, 1) = 0 - ahf_exp(nx, :,nz) = 0 - ahf_exp( :, 1, 1) = 0 - ahf_exp( :,ny, 1) = 0 - ahf_exp( :, 1,nz) = 0 - ahf_exp( :,ny,nz) = 0 - - ahfgradx( 1, 1, :) = 0 - ahfgradx( 1,ny, :) = 0 - ahfgradx( 1, :, 1) = 0 - ahfgradx( 1, :,nz) = 0 - ahfgradx(nx, 1, :) = 0 - ahfgradx(nx,ny, :) = 0 - ahfgradx(nx, :, 1) = 0 - ahfgradx(nx, :,nz) = 0 - ahfgradx( :, 1, 1) = 0 - ahfgradx( :,ny, 1) = 0 - ahfgradx( :, 1,nz) = 0 - ahfgradx( :,ny,nz) = 0 - - ahfgrady( 1, 1, :) = 0 - ahfgrady( 1,ny, :) = 0 - ahfgrady( 1, :, 1) = 0 - ahfgrady( 1, :,nz) = 0 - ahfgrady(nx, 1, :) = 0 - ahfgrady(nx,ny, :) = 0 - ahfgrady(nx, :, 1) = 0 - ahfgrady(nx, :,nz) = 0 - ahfgrady( :, 1, 1) = 0 - ahfgrady( :,ny, 1) = 0 - ahfgrady( :, 1,nz) = 0 - ahfgrady( :,ny,nz) = 0 - - ahfgradz( 1, 1, :) = 0 - ahfgradz( 1,ny, :) = 0 - ahfgradz( 1, :, 1) = 0 - ahfgradz( 1, :,nz) = 0 - ahfgradz(nx, 1, :) = 0 - ahfgradz(nx,ny, :) = 0 - ahfgradz(nx, :, 1) = 0 - ahfgradz(nx, :,nz) = 0 - ahfgradz( :, 1, 1) = 0 - ahfgradz( :,ny, 1) = 0 - ahfgradz( :, 1,nz) = 0 - ahfgradz( :,ny,nz) = 0 - - ahfgradn( 1, 1, :) = 0 - ahfgradn( 1,ny, :) = 0 - ahfgradn( 1, :, 1) = 0 - ahfgradn( 1, :,nz) = 0 - ahfgradn(nx, 1, :) = 0 - ahfgradn(nx,ny, :) = 0 - ahfgradn(nx, :, 1) = 0 - ahfgradn(nx, :,nz) = 0 - ahfgradn( :, 1, 1) = 0 - ahfgradn( :,ny, 1) = 0 - ahfgradn( :, 1,nz) = 0 - ahfgradn( :,ny,nz) = 0 +! Zero out the edges before doing the boundaries. + + ahf_exp( 1, 1, :) = 0.0d0 + ahf_exp( 1,ny, :) = 0.0d0 + ahf_exp( 1, :, 1) = 0.0d0 + ahf_exp( 1, :,nz) = 0.0d0 + ahf_exp(nx, 1, :) = 0.0d0 + ahf_exp(nx,ny, :) = 0.0d0 + ahf_exp(nx, :, 1) = 0.0d0 + ahf_exp(nx, :,nz) = 0.0d0 + ahf_exp( :, 1, 1) = 0.0d0 + ahf_exp( :,ny, 1) = 0.0d0 + ahf_exp( :, 1,nz) = 0.0d0 + ahf_exp( :,ny,nz) = 0.0d0 + + ahfgradx( 1, 1, :) = 0.0d0 + ahfgradx( 1,ny, :) = 0.0d0 + ahfgradx( 1, :, 1) = 0.0d0 + ahfgradx( 1, :,nz) = 0.0d0 + ahfgradx(nx, 1, :) = 0.0d0 + ahfgradx(nx,ny, :) = 0.0d0 + ahfgradx(nx, :, 1) = 0.0d0 + ahfgradx(nx, :,nz) = 0.0d0 + ahfgradx( :, 1, 1) = 0.0d0 + ahfgradx( :,ny, 1) = 0.0d0 + ahfgradx( :, 1,nz) = 0.0d0 + ahfgradx( :,ny,nz) = 0.0d0 + + ahfgrady( 1, 1, :) = 0.0d0 + ahfgrady( 1,ny, :) = 0.0d0 + ahfgrady( 1, :, 1) = 0.0d0 + ahfgrady( 1, :,nz) = 0.0d0 + ahfgrady(nx, 1, :) = 0.0d0 + ahfgrady(nx,ny, :) = 0.0d0 + ahfgrady(nx, :, 1) = 0.0d0 + ahfgrady(nx, :,nz) = 0.0d0 + ahfgrady( :, 1, 1) = 0.0d0 + ahfgrady( :,ny, 1) = 0.0d0 + ahfgrady( :, 1,nz) = 0.0d0 + ahfgrady( :,ny,nz) = 0.0d0 + + ahfgradz( 1, 1, :) = 0.0d0 + ahfgradz( 1,ny, :) = 0.0d0 + ahfgradz( 1, :, 1) = 0.0d0 + ahfgradz( 1, :,nz) = 0.0d0 + ahfgradz(nx, 1, :) = 0.0d0 + ahfgradz(nx,ny, :) = 0.0d0 + ahfgradz(nx, :, 1) = 0.0d0 + ahfgradz(nx, :,nz) = 0.0d0 + ahfgradz( :, 1, 1) = 0.0d0 + ahfgradz( :,ny, 1) = 0.0d0 + ahfgradz( :, 1,nz) = 0.0d0 + ahfgradz( :,ny,nz) = 0.0d0 + + ahfgradn( 1, 1, :) = 0.0d0 + ahfgradn( 1,ny, :) = 0.0d0 + ahfgradn( 1, :, 1) = 0.0d0 + ahfgradn( 1, :,nz) = 0.0d0 + ahfgradn(nx, 1, :) = 0.0d0 + ahfgradn(nx,ny, :) = 0.0d0 + ahfgradn(nx, :, 1) = 0.0d0 + ahfgradn(nx, :,nz) = 0.0d0 + ahfgradn( :, 1, 1) = 0.0d0 + ahfgradn( :,ny, 1) = 0.0d0 + ahfgradn( :, 1,nz) = 0.0d0 + ahfgradn( :,ny,nz) = 0.0d0 ! Boundaries on x direction. - ahf_exp(1,2:ny-1,2:nz-1) = 2.0D0*ahf_exp(2,2:ny-1,2:nz-1) - ahf_exp(3,2:ny-1,2:nz-1) - ahfgradx(1,2:ny-1,2:nz-1) = 2.0D0*ahfgradx(2,2:ny-1,2:nz-1) - ahfgradx(3,2:ny-1,2:nz-1) - ahfgrady(1,2:ny-1,2:nz-1) = 2.0D0*ahfgrady(2,2:ny-1,2:nz-1) - ahfgrady(3,2:ny-1,2:nz-1) - ahfgradz(1,2:ny-1,2:nz-1) = 2.0D0*ahfgradz(2,2:ny-1,2:nz-1) - ahfgradz(3,2:ny-1,2:nz-1) - ahfgradn(1,2:ny-1,2:nz-1) = 2.0D0*ahfgradn(2,2:ny-1,2:nz-1) - ahfgradn(3,2:ny-1,2:nz-1) - - ahf_exp(nx,2:ny-1,2:nz-1) = 2.0D0*ahf_exp(nx-1,2:ny-1,2:nz-1) - ahf_exp(nx-2,2:ny-1,2:nz-1) - ahfgradx(nx,2:ny-1,2:nz-1) = 2.0D0*ahfgradx(nx-1,2:ny-1,2:nz-1) - ahfgradx(nx-2,2:ny-1,2:nz-1) - ahfgrady(nx,2:ny-1,2:nz-1) = 2.0D0*ahfgrady(nx-1,2:ny-1,2:nz-1) - ahfgrady(nx-2,2:ny-1,2:nz-1) - ahfgradz(nx,2:ny-1,2:nz-1) = 2.0D0*ahfgradz(nx-1,2:ny-1,2:nz-1) - ahfgradz(nx-2,2:ny-1,2:nz-1) - ahfgradn(nx,2:ny-1,2:nz-1) = 2.0D0*ahfgradn(nx-1,2:ny-1,2:nz-1) - ahfgradn(nx-2,2:ny-1,2:nz-1) + ahf_exp(1,2:ny-1,2:nz-1) = 2.0D0*ahf_exp(2,2:ny-1,2:nz-1) + . - ahf_exp(3,2:ny-1,2:nz-1) + ahfgradx(1,2:ny-1,2:nz-1) = 2.0D0*ahfgradx(2,2:ny-1,2:nz-1) + . - ahfgradx(3,2:ny-1,2:nz-1) + ahfgrady(1,2:ny-1,2:nz-1) = 2.0D0*ahfgrady(2,2:ny-1,2:nz-1) + . - ahfgrady(3,2:ny-1,2:nz-1) + ahfgradz(1,2:ny-1,2:nz-1) = 2.0D0*ahfgradz(2,2:ny-1,2:nz-1) + . - ahfgradz(3,2:ny-1,2:nz-1) + ahfgradn(1,2:ny-1,2:nz-1) = 2.0D0*ahfgradn(2,2:ny-1,2:nz-1) + . - ahfgradn(3,2:ny-1,2:nz-1) + + ahf_exp(nx,2:ny-1,2:nz-1) = 2.0D0*ahf_exp(nx-1,2:ny-1,2:nz-1) + . - ahf_exp(nx-2,2:ny-1,2:nz-1) + ahfgradx(nx,2:ny-1,2:nz-1) = 2.0D0*ahfgradx(nx-1,2:ny-1,2:nz-1) + . - ahfgradx(nx-2,2:ny-1,2:nz-1) + ahfgrady(nx,2:ny-1,2:nz-1) = 2.0D0*ahfgrady(nx-1,2:ny-1,2:nz-1) + . - ahfgrady(nx-2,2:ny-1,2:nz-1) + ahfgradz(nx,2:ny-1,2:nz-1) = 2.0D0*ahfgradz(nx-1,2:ny-1,2:nz-1) + . - ahfgradz(nx-2,2:ny-1,2:nz-1) + ahfgradn(nx,2:ny-1,2:nz-1) = 2.0D0*ahfgradn(nx-1,2:ny-1,2:nz-1) + . - ahfgradn(nx-2,2:ny-1,2:nz-1) ! Boundaries on y direction. - ahf_exp(2:nx-1,1,2:nz-1) = 2.0D0*ahf_exp(2:nx-1,2,2:nz-1) - ahf_exp(2:nx-1,3,2:nz-1) - ahfgradx(2:nx-1,1,2:nz-1) = 2.0D0*ahfgradx(2:nx-1,2,2:nz-1) - ahfgradx(2:nx-1,3,2:nz-1) - ahfgrady(2:nx-1,1,2:nz-1) = 2.0D0*ahfgrady(2:nx-1,2,2:nz-1) - ahfgrady(2:nx-1,3,2:nz-1) - ahfgradz(2:nx-1,1,2:nz-1) = 2.0D0*ahfgradz(2:nx-1,2,2:nz-1) - ahfgradz(2:nx-1,3,2:nz-1) - ahfgradn(2:nx-1,1,2:nz-1) = 2.0D0*ahfgradn(2:nx-1,2,2:nz-1) - ahfgradn(2:nx-1,3,2:nz-1) - - ahf_exp(2:nx-1,ny,2:nz-1) = 2.0D0*ahf_exp(2:nx-1,ny-1,2:nz-1) - ahf_exp(2:nx-1,ny-2,2:nz-1) - ahfgradx(2:nx-1,ny,2:nz-1) = 2.0D0*ahfgradx(2:nx-1,ny-1,2:nz-1) - ahfgradx(2:nx-1,ny-2,2:nz-1) - ahfgrady(2:nx-1,ny,2:nz-1) = 2.0D0*ahfgrady(2:nx-1,ny-1,2:nz-1) - ahfgrady(2:nx-1,ny-2,2:nz-1) - ahfgradz(2:nx-1,ny,2:nz-1) = 2.0D0*ahfgradz(2:nx-1,ny-1,2:nz-1) - ahfgradz(2:nx-1,ny-2,2:nz-1) - ahfgradn(2:nx-1,ny,2:nz-1) = 2.0D0*ahfgradn(2:nx-1,ny-1,2:nz-1) - ahfgradn(2:nx-1,ny-2,2:nz-1) + ahf_exp(2:nx-1,1,2:nz-1) = 2.0D0*ahf_exp(2:nx-1,2,2:nz-1) + . - ahf_exp(2:nx-1,3,2:nz-1) + ahfgradx(2:nx-1,1,2:nz-1) = 2.0D0*ahfgradx(2:nx-1,2,2:nz-1) + . - ahfgradx(2:nx-1,3,2:nz-1) + ahfgrady(2:nx-1,1,2:nz-1) = 2.0D0*ahfgrady(2:nx-1,2,2:nz-1) + . - ahfgrady(2:nx-1,3,2:nz-1) + ahfgradz(2:nx-1,1,2:nz-1) = 2.0D0*ahfgradz(2:nx-1,2,2:nz-1) + . - ahfgradz(2:nx-1,3,2:nz-1) + ahfgradn(2:nx-1,1,2:nz-1) = 2.0D0*ahfgradn(2:nx-1,2,2:nz-1) + . - ahfgradn(2:nx-1,3,2:nz-1) + + ahf_exp(2:nx-1,ny,2:nz-1) = 2.0D0*ahf_exp(2:nx-1,ny-1,2:nz-1) + . - ahf_exp(2:nx-1,ny-2,2:nz-1) + ahfgradx(2:nx-1,ny,2:nz-1) = 2.0D0*ahfgradx(2:nx-1,ny-1,2:nz-1) + . - ahfgradx(2:nx-1,ny-2,2:nz-1) + ahfgrady(2:nx-1,ny,2:nz-1) = 2.0D0*ahfgrady(2:nx-1,ny-1,2:nz-1) + . - ahfgrady(2:nx-1,ny-2,2:nz-1) + ahfgradz(2:nx-1,ny,2:nz-1) = 2.0D0*ahfgradz(2:nx-1,ny-1,2:nz-1) + . - ahfgradz(2:nx-1,ny-2,2:nz-1) + ahfgradn(2:nx-1,ny,2:nz-1) = 2.0D0*ahfgradn(2:nx-1,ny-1,2:nz-1) + . - ahfgradn(2:nx-1,ny-2,2:nz-1) ! Boundaries on z direction. - ahf_exp(2:nx-1,2:ny-1,1) = 2.0D0*ahf_exp(2:nx-1,2:ny-1,2) - ahf_exp(2:nx-1,2:ny-1,3) - ahfgradx(2:nx-1,2:ny-1,1) = 2.0D0*ahfgradx(2:nx-1,2:ny-1,2) - ahfgradx(2:nx-1,2:ny-1,3) - ahfgrady(2:nx-1,2:ny-1,1) = 2.0D0*ahfgrady(2:nx-1,2:ny-1,2) - ahfgrady(2:nx-1,2:ny-1,3) - ahfgradz(2:nx-1,2:ny-1,1) = 2.0D0*ahfgradz(2:nx-1,2:ny-1,2) - ahfgradz(2:nx-1,2:ny-1,3) - ahfgradn(2:nx-1,2:ny-1,1) = 2.0D0*ahfgradn(2:nx-1,2:ny-1,2) - ahfgradn(2:nx-1,2:ny-1,3) - - ahf_exp(2:nx-1,2:ny-1,nz) = 2.0D0*ahf_exp(2:nx-1,2:ny-1,nz-1) - ahf_exp(2:nx-1,2:ny-1,nz-2) - ahfgradx(2:nx-1,2:ny-1,nz) = 2.0D0*ahfgradx(2:nx-1,2:ny-1,nz-1) - ahfgradx(2:nx-1,2:ny-1,nz-2) - ahfgrady(2:nx-1,2:ny-1,nz) = 2.0D0*ahfgrady(2:nx-1,2:ny-1,nz-1) - ahfgrady(2:nx-1,2:ny-1,nz-2) - ahfgradz(2:nx-1,2:ny-1,nz) = 2.0D0*ahfgradz(2:nx-1,2:ny-1,nz-1) - ahfgradz(2:nx-1,2:ny-1,nz-2) - ahfgradn(2:nx-1,2:ny-1,nz) = 2.0D0*ahfgradn(2:nx-1,2:ny-1,nz-1) - ahfgradn(2:nx-1,2:ny-1,nz-2) + ahf_exp(2:nx-1,2:ny-1,1) = 2.0D0*ahf_exp(2:nx-1,2:ny-1,2) + . - ahf_exp(2:nx-1,2:ny-1,3) + ahfgradx(2:nx-1,2:ny-1,1) = 2.0D0*ahfgradx(2:nx-1,2:ny-1,2) + . - ahfgradx(2:nx-1,2:ny-1,3) + ahfgrady(2:nx-1,2:ny-1,1) = 2.0D0*ahfgrady(2:nx-1,2:ny-1,2) + . - ahfgrady(2:nx-1,2:ny-1,3) + ahfgradz(2:nx-1,2:ny-1,1) = 2.0D0*ahfgradz(2:nx-1,2:ny-1,2) + . - ahfgradz(2:nx-1,2:ny-1,3) + ahfgradn(2:nx-1,2:ny-1,1) = 2.0D0*ahfgradn(2:nx-1,2:ny-1,2) + . - ahfgradn(2:nx-1,2:ny-1,3) + + ahf_exp(2:nx-1,2:ny-1,nz) = 2.0D0*ahf_exp(2:nx-1,2:ny-1,nz-1) + . - ahf_exp(2:nx-1,2:ny-1,nz-2) + ahfgradx(2:nx-1,2:ny-1,nz) = 2.0D0*ahfgradx(2:nx-1,2:ny-1,nz-1) + . - ahfgradx(2:nx-1,2:ny-1,nz-2) + ahfgrady(2:nx-1,2:ny-1,nz) = 2.0D0*ahfgrady(2:nx-1,2:ny-1,nz-1) + . - ahfgrady(2:nx-1,2:ny-1,nz-2) + ahfgradz(2:nx-1,2:ny-1,nz) = 2.0D0*ahfgradz(2:nx-1,2:ny-1,nz-1) + . - ahfgradz(2:nx-1,2:ny-1,nz-2) + ahfgradn(2:nx-1,2:ny-1,nz) = 2.0D0*ahfgradn(2:nx-1,2:ny-1,nz-1) + . - ahfgradn(2:nx-1,2:ny-1,nz-2) ! Synchronize. diff --git a/src/AHFinder_gau.F b/src/AHFinder_gau.F index 05160c3..28ab392 100644 --- a/src/AHFinder_gau.F +++ b/src/AHFinder_gau.F @@ -505,6 +505,7 @@ avgx = 0.0D0 avgy = 0.0D0 + avgz = 0.0D0 do j=1,nphi+1 do i=1,ntheta+1 @@ -573,9 +574,6 @@ xp = xa(i,j) yp = ya(i,j) zp = za(i,j) - - avgx = avgx + xp - avgy = avgy + yp if ((xp.gt.xmx).or.(xp.lt.xmn).or. . (yp.gt.ymx).or.(yp.lt.ymn).or. @@ -583,22 +581,39 @@ error2 = 1 end if +! Add to sums for average position. + + avgx = avgx + xp + avgy = avgy + yp + avgz = avgz + zp + end do end do +! Find position of horizon centroid. + avgx = avgx/dble(npoints) avgy = avgy/dble(npoints) + avgz = avgz/dble(npoints) if ((.not.find3).or. . (find3.and.(mfind-1.eq.horizon_centroid))) then + ahf_centroid_x = avgx ahf_centroid_y = avgy +! Take symmetries into account. + + if (refx) ahf_centroid_x = xc + if (refy) ahf_centroid_y = yc + if (verbose) then + write(*,*) write(6,*) 'AHFinder_gau: ahf_centroid_x = ', ahf_centroid_x write(6,*) 'AHFinder_gau: ahf_centroid_y = ', ahf_centroid_y end if - endif + + end if ! Other processors. -- cgit v1.2.3