diff options
Diffstat (limited to 'src/ADMmass_integrand3D.F')
| -rw-r--r-- | src/ADMmass_integrand3D.F | 170 |
1 files changed, 170 insertions, 0 deletions
diff --git a/src/ADMmass_integrand3D.F b/src/ADMmass_integrand3D.F new file mode 100644 index 0000000..06b1198 --- /dev/null +++ b/src/ADMmass_integrand3D.F @@ -0,0 +1,170 @@ +#include "cctk.h" + +c ======================================================================== + + SUBROUTINE ADMmass_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,gxy,gxz, + & gyy,gyz,gzz,ADMmass_int,Psi,Psi_power) + +c ------------------------------------------------------------------------ +c +c Estimates the ADM mass at a given radius using Equation (7) from +c O Murchadha and York, Phys Rev D, 10, 1974 page 2345 +c +c ------------------------------------------------------------------------ + + IMPLICIT NONE + +c Input variables + INTEGER,INTENT(IN) :: + & Psi_power + CCTK_REAL,INTENT(IN) :: + & Dx,Dy,Dz,origin(3) + CCTK_REAL,DIMENSION(:),INTENT(IN) :: + & x,y,z + CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) :: + & gxx,gxy,gxz,gyy,gyz,gzz,Psi + +c Output variables + CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) :: + & ADMmass_int + +c Local variables, here only + INTEGER :: + & i,j,k,ip + CCTK_REAL,PARAMETER :: + & half = 0.5D0 + CCTK_REAL :: + & rad,ux,uy,uz,det,dxx,dxy,dxz,dyy,dyz,dzz,uxx,uxy,uxz,uyy, + & uyz,uzz,term1,term2,term3,dxx_y,dxx_z,dxy_x,dxy_y,dxy_z, + & dyy_x,dxz_x,dxz_y,dxz_z,dyz_x,dzz_x,dyy_z,dyz_y,dyz_z,dzz_y, + & Pi,idet,p,pip,pim,pjp,pjm,pkp,pkm + +c ------------------------------------------------------------------------ + + Pi = ACOS(-1D0) + +c Because other codes evolve Psi**4 + SELECT CASE (Psi_power) + + CASE (1) + ip = 4 + + CASE (4) + ip = 1 + + CASE DEFAULT + WRITE(*,*) "This value of Psi_power is not supported" + + END SELECT + + + DO k = 2, SIZE(z)-1 + DO j = 2, SIZE(y)-1 + DO i = 2, SIZE(x)-1 + + rad = SQRT((x(i)-origin(1))**2 + & +(y(j)-origin(2))**2 + & +(z(k)-origin(3))**2) + + IF (rad.NE.0) THEN + + ux = (x(i)-origin(1))/rad + uy = (y(j)-origin(2))/rad + uz = (z(k)-origin(3))/rad + +c Abbreviations for metric coefficients +c ------------------------------------- + p = psi(i,j,k)**ip + + dxx = p*gxx(i,j,k); dxy = p*gxy(i,j,k) + dxz = p*gxz(i,j,k); dyy = p*gyy(i,j,k) + dyz = p*gyz(i,j,k); dzz = p*gzz(i,j,k) + +c Determinant of 3-metric +c ----------------------- + det = (dxx*dyy*dzz + 2.0D0*dxy*dxz*dyz + & - (dxx*dyz**2 + dyy*dxz**2 + dzz*dxy**2)) + + idet = 1.0/det + +c Inverse 3-metric +c ---------------- + uxx = idet*(dyy*dzz - dyz**2) + uyy = idet*(dxx*dzz - dxz**2) + uzz = idet*(dxx*dyy - dxy**2) + uxy = idet*(dxz*dyz - dzz*dxy) + uxz = idet*(dxy*dyz - dyy*dxz) + uyz = idet*(dxy*dxz - dxx*dyz) + +c Derivatives of the 3-metric +c --------------------------- + pip = psi(i+1,j,k)**ip + pim = psi(i-1,j,k)**ip + pjp = psi(i,j+1,k)**ip + pjm = psi(i,j-1,k)**ip + pkp = psi(i,j,k+1)**ip + pkm = psi(i,j,k-1)**ip + + dxx_y = half/Dy*(pjp*gxx(i,j+1,k)-pjm*gxx(i,j-1,k)) + dxx_z = half/Dz*(pkp*gxx(i,j,k+1)-pkm*gxx(i,j,k-1)) + dxy_x = half/Dx*(pip*gxy(i+1,j,k)-pim*gxy(i-1,j,k)) + dxy_y = half/Dy*(pjp*gxy(i,j+1,k)-pjm*gxy(i,j-1,k)) + dxy_z = half/Dz*(pkp*gxy(i,j,k+1)-pkm*gxy(i,j,k-1)) + dyy_x = half/Dx*(pip*gyy(i+1,j,k)-pim*gyy(i-1,j,k)) + dyy_z = half/Dz*(pkp*gyy(i,j,k+1)-pkm*gyy(i,j,k-1)) + dxz_x = half/Dx*(pip*gxz(i+1,j,k)-pim*gxz(i-1,j,k)) + dxz_y = half/Dy*(pjp*gxz(i,j+1,k)-pjm*gxz(i,j-1,k)) + dxz_z = half/Dz*(pkp*gxz(i,j,k+1)-pkm*gxz(i,j,k-1)) + dyz_x = half/Dx*(pip*gyz(i+1,j,k)-pim*gyz(i-1,j,k)) + dyz_y = half/Dy*(pjp*gyz(i,j+1,k)-pjm*gyz(i,j-1,k)) + dyz_z = half/Dz*(pkp*gyz(i,j,k+1)-pkm*gyz(i,j,k-1)) + dzz_x = half/Dx*(pip*gzz(i+1,j,k)-pim*gzz(i-1,j,k)) + dzz_y = half/Dy*(pjp*gzz(i,j+1,k)-pjm*gzz(i,j-1,k)) + + term1 = uxy*(dxx_y-dxy_x)+uxz*(dxx_z-dxz_x) + & +uyy*(dxy_y-dyy_x)+uyz*(dxy_z-dyz_x) + & +uyz*(dxz_y-dyz_x)+uzz*(dxz_z-dzz_x) + + term2 = uyz*(dyy_z-dyz_y)+uxy*(dyy_x-dxy_y) + & +uzz*(dyz_z-dzz_y)+uxz*(dyz_x-dxz_y) + & +uxz*(dxy_z-dxz_y)+uxx*(dxy_x-dxx_y) + + term3 = uxz*(dzz_x-dxz_z)+uyz*(dzz_y-dyz_z) + & +uxx*(dxz_x-dxx_z)+uxy*(dxz_y-dxy_z) + & +uxy*(dyz_x-dxy_z)+uyy*(dyz_y-dyy_z) + + ADMmass_int(i,j,k) = 1.0D0/16.0D0/Pi*(ux*term1+ + & uy*term2+uz*term3)*SQRT(det)*rad**2 + + ELSE + + ADMmass_int(i,j,k) = 0.0D0 + + ENDIF + + ENDDO + ENDDO + ENDDO + + +c This is needed when the grid is an octant, but it does not hurt +c if it is not + + DO k = 2, size(z)-1 + DO j = 2, size(y)-1 + ADMmass_int(1,j,k) = ADMmass_int(2,j,k) + ENDDO + ENDDO + DO k = 2, size(z)-1 + DO i = 1, size(x)-1 + ADMmass_int(i,1,k) = ADMmass_int(i,2,k) + ENDDO + ENDDO + DO j = 1, size(y)-1 + DO i = 1, size(x)-1 + ADMmass_int(i,j,1) = ADMmass_int(i,j,2) + ENDDO + ENDDO + + + END SUBROUTINE ADMmass_integrand3D |