c/*@@ c @file AHFinder_exp.F c @date April 1998 c @author Miguel Alcubierre c @desc c Here I calculate the (normalized) expansion of the c horizon function. The expansion is defined by the c following expression (a=lapse, b=shift): c c __2 2 / ab __a__b \ c exp = \/ f / u + d f d f / u | K - \/ \/ f / u | - trK c a b \ / c c where: c c / mn \ 1/2 c u = | g d f d f | c \ m n / c c @enddesc c@@*/ c Note that including cactus.h will also include AHFinder.h !#include "cactus.h" #include "cctk.h" #include "cctk_parameters.h" #include "cctk_arguments.h" subroutine AHFinder_exp(CCTK_FARGUMENTS) use AHFinder_dat implicit none DECLARE_CCTK_FARGUMENTS DECLARE_CCTK_PARAMETERS integer i,j,k INTEGER CCTK_Equals CCTK_REAL det,idet CCTK_REAL gdxx,gdyy,gdzz,gdxy,gdxz,gdyz CCTK_REAL guxx,guyy,guzz,guxy,guxz,guyz CCTK_REAL kdxx,kdyy,kdzz,kdxy,kdxz,kdyz CCTK_REAL ddxxx,ddxyy,ddxzz,ddxxy,ddxxz,ddxyz CCTK_REAL ddyxx,ddyyy,ddyzz,ddyxy,ddyxz,ddyyz CCTK_REAL ddzxx,ddzyy,ddzzz,ddzxy,ddzxz,ddzyz CCTK_REAL d2fxx,d2fyy,d2fzz,d2fxy,d2fxz,d2fyz 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 zero,half,one,two INTEGER istat ! Description of variables: ! ! i,j,k, counters ! ! gdij g ! ij ! ij ! guij g ! ! ! det det(g) ! ! idet 1/det ! ! ! kdij K ! ij ! ! ddijk 1/2 d g ! i jk ! ! dfi d f ! i ! ij ! dfui g d f ! j ! ! d2fij d d f ! i j ! ! c2fij f ! ;ij ! ! aux,T* auxilliary variables ! *************************************** ! *** DEFINE {zero,half,one,two} *** ! *************************************** zero = 0.0d0 half = 0.5d0 one = 1.0d0 two = 2.0d0 ! ************************** ! *** FIND EXPANSION *** ! ************************** idx = half/dx 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 do k=2,nz-1 do j=2,ny-1 do i=2,nx-1 ! Find spatial metric. gdxx = gxx(i,j,k) gdyy = gyy(i,j,k) gdzz = gzz(i,j,k) gdxy = gxy(i,j,k) gdxz = gxz(i,j,k) gdyz = gyz(i,j,k) ! Find extrinsic curvature. kdxx = kxx(i,j,k) kdyy = kyy(i,j,k) kdzz = kzz(i,j,k) kdxy = kxy(i,j,k) kdxz = kxz(i,j,k) kdyz = kyz(i,j,k) ! Find Ds using finite differences. ddxxx = idx*(gxx(i+1,j,k) - gxx(i-1,j,k)) ddxyy = idx*(gyy(i+1,j,k) - gyy(i-1,j,k)) ddxzz = idx*(gzz(i+1,j,k) - gzz(i-1,j,k)) ddxxy = idx*(gxy(i+1,j,k) - gxy(i-1,j,k)) ddxxz = idx*(gxz(i+1,j,k) - gxz(i-1,j,k)) ddxyz = idx*(gyz(i+1,j,k) - gyz(i-1,j,k)) ddyxx = idy*(gxx(i,j+1,k) - gxx(i,j-1,k)) ddyyy = idy*(gyy(i,j+1,k) - gyy(i,j-1,k)) ddyzz = idy*(gzz(i,j+1,k) - gzz(i,j-1,k)) ddyxy = idy*(gxy(i,j+1,k) - gxy(i,j-1,k)) ddyxz = idy*(gxz(i,j+1,k) - gxz(i,j-1,k)) ddyyz = idy*(gyz(i,j+1,k) - gyz(i,j-1,k)) ddzxx = idz*(gxx(i,j,k+1) - gxx(i,j,k-1)) ddzyy = idz*(gyy(i,j,k+1) - gyy(i,j,k-1)) ddzzz = idz*(gzz(i,j,k+1) - gzz(i,j,k-1)) ddzxy = idz*(gxy(i,j,k+1) - gxy(i,j,k-1)) ddzxz = idz*(gxz(i,j,k+1) - gxz(i,j,k-1)) ddzyz = idz*(gyz(i,j,k+1) - gyz(i,j,k-1)) ! Find determinant of spatial metric. det = gdxx*gdyy*gdzz + two*gdxy*gdxz*gdyz . - gdxx*gdyz**2 - gdyy*gdxz**2 - gdzz*gdxy**2 if (det.le.zero) then ahf_exp(i,j,k) = zero goto 10 end if idet = one/det ! Find inverse spatial metric. guxx = idet*(gdyy*gdzz - gdyz**2) guyy = idet*(gdxx*gdzz - gdxz**2) guzz = idet*(gdxx*gdyy - gdxy**2) guxy = idet*(gdxz*gdyz - gdzz*gdxy) guxz = idet*(gdxy*gdyz - gdyy*gdxz) guyz = idet*(gdxy*gdxz - gdxx*gdyz) ! Find spatial derivatives of f. T0 = two*ahfgrid(i,j,k) T1 = ahfgrid(i+1,j,k) T2 = ahfgrid(i-1,j,k) dfx = (T1 - T2)*idx d2fxx = (T1 - T0 + T2)*idx2 T1 = ahfgrid(i,j+1,k) T2 = ahfgrid(i,j-1,k) dfy = (T1 - T2)*idy d2fyy = (T1 - T0 + T2)*idy2 T1 = ahfgrid(i,j,k+1) T2 = ahfgrid(i,j,k-1) 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. ahfgradx(i,j,k) = dfx ahfgrady(i,j,k) = dfy ahfgradz(i,j,k) = dfz aux = guxx*dfx**2 + guyy*dfy**2 + guzz*dfz**2 . + two*(guxy*dfx*dfy + guxz*dfx*dfz + guyz*dfy*dfz) ahfgradn(i,j,k) = dsqrt(aux) if (ahfgradn(i,j,k).eq.zero) ahfgradn(i,j,k) = one ! 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 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 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. 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. c2fxx = d2fxx - half*(dfux*ddxxx . + dfuy*(two*ddxxy - ddyxx) . + dfuz*(two*ddxxz - ddzxx)) c2fyy = d2fyy - half*(dfuy*ddyyy . + dfux*(two*ddyxy - ddxyy) . + dfuz*(two*ddyyz - ddzyy)) c2fzz = d2fzz - half*(dfuz*ddzzz . + dfux*(two*ddzxz - ddxzz) . + dfuy*(two*ddzyz - ddyzz)) c2fxy = d2fxy - half*(dfux*ddyxx + dfuy*ddxyy . + dfuz*(ddxyz + ddyxz - ddzxy)) c2fxz = d2fxz - half*(dfux*ddzxx + dfuz*ddxzz . + dfuy*(ddxyz + ddzxy - ddyxz)) c2fyz = d2fyz - half*(dfuy*ddzyy + dfuz*ddyzz . + dfux*(ddyxz + ddzxy - ddxyz)) ! / m \ 1/2 ! Find: u = | d f d f | ! \ m / T0 = dfx*dfux + dfy*dfuy + dfz*dfuz if (T0.gt.zero) then T0 = one/dsqrt(T0) else ahf_exp(i,j,k) = zero goto 10 end if ! __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 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 T3 = c2fxx*dfux**2 + c2fyy*dfuy**2 + c2fzz*dfuz**2 . + two*(dfux*(c2fxy*dfuy + c2fxz*dfuz) . + c2fyz*dfuy*dfuz) T3 = T3*T0**3 ! Find: trK T4 = guxx*kdxx + guyy*kdyy + guzz*kdzz . + two*(guxy*kdxy + guxz*kdxz + guyz*kdyz) ! Find the expansion. ahf_exp(i,j,k) = T1 + T2 - T3 - T4 ! Check for NaNs. if (ahf_exp(i,j,k).ge.zero) then else if (ahf_exp(i,j,k).lt.zero) then else if (myproc.eq.0) then write(*,*) 'NaN in ahf_exp at point',i,j,k end if STOP end if 10 continue end do end do end do ! Boundaries. do k=1,nz do j=1,ny ahf_exp(1 ,j,k) = ahf_exp(2 ,j,k) ahf_exp(nx,j,k) = ahf_exp(nx-1,j,k) ahfgradx(nx,j,k) = ahfgradx(nx-1,j,k) ahfgrady(nx,j,k) = ahfgrady(nx-1,j,k) ahfgradz(nx,j,k) = ahfgradz(nx-1,j,k) ahfgradn(nx,j,k) = ahfgradn(nx-1,j,k) ahfgrady(1 ,j,k) = ahfgrady(2,j,k) ahfgradz(1 ,j,k) = ahfgradz(2,j,k) ahfgradn(1 ,j,k) = ahfgradn(2,j,k) ! if ((contains("grid","octant").ne.0).or. ! . (contains("grid","quadrant").ne.0)) then if ((CCTK_Equals(domain,"octant").eq.1).or. . (CCTK_Equals(domain,"quadrant").eq.1)) then ahfgradx(1,j,k) = -ahfgradx(2,j,k) else ahfgradx(1,j,k) = +ahfgradx(2,j,k) end if end do end do do k=1,nz do i=1,nx ahf_exp(i,1 ,k) = ahf_exp(i,2 ,k) ahf_exp(i,ny,k) = ahf_exp(i,ny-1,k) ahfgradx(i,ny,k) = ahfgradx(i,ny-1,k) ahfgrady(i,ny,k) = ahfgrady(i,ny-1,k) ahfgradz(i,ny,k) = ahfgradz(i,ny-1,k) ahfgradn(i,ny,k) = ahfgradn(i,ny-1,k) ahfgradx(i,1 ,k) = ahfgradx(i,2,k) ahfgradz(i,1 ,k) = ahfgradz(i,2,k) ahfgradn(i,1 ,k) = ahfgradn(i,2,k) ! if ((contains("grid","octant").ne.0).or. ! . (contains("grid","quadrant").ne.0)) then if ((CCTK_Equals(domain,"octant").eq.1).or. . (CCTK_Equals(domain,"quadrant").eq.1)) then ahfgrady(i,1,k) = -ahfgrady(i,2,k) else ahfgrady(i,1,k) = +ahfgrady(i,2,k) end if end do end do do j=1,ny do i=1,nx ahf_exp(i,j,1 ) = ahf_exp(i,j,2 ) ahf_exp(i,j,nz) = ahf_exp(i,j,nz-1) ahfgradx(i,j,nz) = ahfgradx(i,j,nz-1) ahfgrady(i,j,nz) = ahfgrady(i,j,nz-1) ahfgradz(i,j,nz) = ahfgradz(i,j,nz-1) ahfgradn(i,j,nz) = ahfgradn(i,j,nz-1) ahfgradx(i,j,1 ) = ahfgradx(i,j,2) ahfgrady(i,j,1 ) = ahfgrady(i,j,2) ahfgradn(i,j,1 ) = ahfgradn(i,j,2) ! if (contains("grid","octant").ne.0) then if (CCTK_Equals(domain,"octant").eq.1) then ahfgradz(i,j,1) = -ahfgradz(i,j,2) else ahfgradz(i,j,1) = +ahfgradz(i,j,2) end if end do end do ! Synchronize. call CCTK_SyncGroup(cctkGH,"ahfinder::ahfgradient") call CCTK_SyncGroup(cctkGH,"ahfinder::ahfindergrid") ! call synconefunc(ahf_exp) ! call synconefunc(ahfgradx) ! call synconefunc(ahfgrady) ! call synconefunc(ahfgradz) ! call synconefunc(ahfgradn) ! *************** ! *** END *** ! *************** return end