1346 lines
50 KiB
Fortran
1346 lines
50 KiB
Fortran
|
|
|
|
#include "macrodef.fh"
|
|
|
|
!-----------------------------------------------------------------------------
|
|
!
|
|
! compute the Newman-Penrose Weyl scalar Psi4
|
|
! for BSSN dynamical variables
|
|
!
|
|
!-----------------------------------------------------------------------------
|
|
|
|
subroutine getnp4(ex, X, Y, Z, &
|
|
chi, trK, &
|
|
dxx,gxy,gxz,dyy,gyz,dzz, &
|
|
Axx,Axy,Axz,Ayy,Ayz,Azz, &
|
|
Gamxxx,Gamxxy,Gamxxz,Gamxyy,Gamxyz,Gamxzz,&
|
|
Gamyxx,Gamyxy,Gamyxz,Gamyyy,Gamyyz,Gamyzz,&
|
|
Gamzxx,Gamzxy,Gamzxz,Gamzyy,Gamzyz,Gamzzz,&
|
|
Rxx,Rxy,Rxz,Ryy,Ryz,Rzz,&
|
|
Rpsi4, Ipsi4, &
|
|
symmetry)
|
|
|
|
implicit none
|
|
|
|
!~~~~~~> Input parameters:
|
|
|
|
integer,intent(in ):: ex(1:3),symmetry
|
|
real*8, intent(in ):: X(1:ex(1)),Y(1:ex(2)),Z(1:ex(3))
|
|
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: dxx,gxy,gxz,dyy,gyz,dzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: Axx,Axy,Axz,Ayy,Ayz,Azz
|
|
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: chi,trK
|
|
! physical second kind of connection
|
|
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Gamxxx, Gamxxy, Gamxxz
|
|
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Gamxyy, Gamxyz, Gamxzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Gamyxx, Gamyxy, Gamyxz
|
|
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Gamyyy, Gamyyz, Gamyzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Gamzxx, Gamzxy, Gamzxz
|
|
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Gamzyy, Gamzyz, Gamzzz
|
|
! physical Ricci tensor
|
|
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Rxx,Rxy,Rxz,Ryy,Ryz,Rzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)), intent(out):: Rpsi4,Ipsi4
|
|
|
|
!~~~~~~> Other variables:
|
|
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: f,fx,fy,fz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: gxx,gyy,gzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: chix,chiy,chiz,chipn1
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: vx,vy,vz,ux,uy,uz,wx,wy,wz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: Exx,Exy,Exz,Eyy,Eyz,Ezz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: Bxx,Bxy,Bxz,Byy,Byz,Bzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: Axxx,Axxy,Axxz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: Axyx,Axyy,Axyz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: Axzx,Axzy,Axzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: Ayyx,Ayyy,Ayyz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: Ayzx,Ayzy,Ayzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: Azzx,Azzy,Azzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: gupxx,gupxy,gupxz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: gupyy,gupyz,gupzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: uuwwxx,uuwwxy,uuwwxz,uuwwyy,uuwwyz,uuwwzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: uwxx,uwxy,uwxz,uwyy,uwyz,uwzz
|
|
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: adm_dxx,adm_gxy,adm_gxz,adm_dyy,adm_gyz,adm_dzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: Kxx,Kxy,Kxz,Kyy,Kyz,Kzz
|
|
|
|
real*8, parameter :: ZEO = 0.d0, ONE = 1.d0, TWO = 2.d0
|
|
real*8, parameter :: F1o3 = 1.d0/3.d0, FOUR = 4.d0
|
|
real*8, parameter :: SYM = 1.D0, ANTI= - 1.D0
|
|
real*8 :: dX, dY, dZ
|
|
integer::i,j,k
|
|
real*8,parameter::TINYRR=1.d-14
|
|
|
|
dX = X(2) - X(1)
|
|
dY = Y(2) - Y(1)
|
|
dZ = Z(2) - Z(1)
|
|
|
|
gxx = dxx + ONE
|
|
gyy = dyy + ONE
|
|
gzz = dzz + ONE
|
|
chipn1= chi + ONE
|
|
|
|
#if (ABV == 1)
|
|
call bssn2adm(ex,chipn1,trK,gxx,gxy,gxz,gyy,gyz,gzz, &
|
|
Axx,Axy,Axz,Ayy,Ayz,Azz, &
|
|
adm_dxx,adm_gxy,adm_gxz,adm_dyy,adm_gyz,adm_dzz, &
|
|
Kxx,Kxy,Kxz,Kyy,Kyz,Kzz)
|
|
adm_dxx = adm_dxx - ONE
|
|
adm_dyy = adm_dyy - ONE
|
|
adm_dzz = adm_dzz - ONE
|
|
call adm_ricci_gamma(ex, X, Y, Z, &
|
|
adm_dxx,adm_gxy,adm_gxz,adm_dyy,adm_gyz,adm_dzz,&
|
|
Gamxxx,Gamxxy,Gamxxz,Gamxyy,Gamxyz,Gamxzz,&
|
|
Gamyxx,Gamyxy,Gamyxz,Gamyyy,Gamyyz,Gamyzz,&
|
|
Gamzxx,Gamzxy,Gamzxz,Gamzyy,Gamzyz,Gamzzz,&
|
|
Rxx,Rxy,Rxz,Ryy,Ryz,Rzz,&
|
|
Symmetry)
|
|
#endif
|
|
|
|
! invert tilted metric
|
|
gupzz = gxx * gyy * gzz + gxy * gyz * gxz + gxz * gxy * gyz - &
|
|
gxz * gyy * gxz - gxy * gxy * gzz - gxx * gyz * gyz
|
|
gupxx = ( gyy * gzz - gyz * gyz ) / gupzz
|
|
gupxy = - ( gxy * gzz - gyz * gxz ) / gupzz
|
|
gupxz = ( gxy * gyz - gyy * gxz ) / gupzz
|
|
gupyy = ( gxx * gzz - gxz * gxz ) / gupzz
|
|
gupyz = - ( gxx * gyz - gxy * gxz ) / gupzz
|
|
gupzz = ( gxx * gyy - gxy * gxy ) / gupzz
|
|
|
|
! initialize U, V, W vetors
|
|
#if (tetradtype == 0)
|
|
do i=1,ex(1)
|
|
do j=1,ex(2)
|
|
do k=1,ex(3)
|
|
if(abs(X(i)) < TINYRR .and. abs(Y(j)) < TINYRR .and. abs(Z(k)) < TINYRR)then
|
|
vx(i,j,k) = TINYRR
|
|
vy(i,j,k) = TINYRR
|
|
vz(i,j,k) = TINYRR
|
|
else
|
|
vx(i,j,k) = X(i)
|
|
vy(i,j,k) = Y(j)
|
|
vz(i,j,k) = Z(k)
|
|
endif
|
|
if(abs(X(i)) < TINYRR .and. abs(Y(j)) < TINYRR)then
|
|
ux(i,j,k) = - TINYRR
|
|
uy(i,j,k) = TINYRR
|
|
uz(i,j,k) = ZEO
|
|
wx(i,j,k) = TINYRR*Z(k)
|
|
wy(i,j,k) = TINYRR*Z(k)
|
|
wz(i,j,k) = -2*TINYRR*TINYRR
|
|
else
|
|
ux(i,j,k) = - Y(j)
|
|
uy(i,j,k) = X(i)
|
|
uz(i,j,k) = ZEO
|
|
wx(i,j,k) = X(i)*Z(k)
|
|
wy(i,j,k) = Y(j)*Z(k)
|
|
wz(i,j,k) = -(X(i)*X(i) + Y(j)*Y(j))
|
|
endif
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
f = 1.d0/chipn1
|
|
|
|
fx = gxx*vx*vx + gyy*vy*vy + gzz*vz*vz &
|
|
+(gxy*vx*vy + gxz*vx*vz + gyz*vy*vz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
vx = vx/fx
|
|
vy = vy/fx
|
|
vz = vz/fx
|
|
|
|
fx = gxx*vx*ux + gxy*vx*uy + gxz*vx*uz + &
|
|
gxy*vy*ux + gyy*vy*uy + gyz*vy*uz + &
|
|
gxz*vz*ux + gyz*vz*uy + gzz*vz*uz
|
|
fx = fx*f
|
|
ux = ux - fx*vx
|
|
uy = uy - fx*vy
|
|
uz = uz - fx*vz
|
|
fx = gxx*ux*ux + gyy*uy*uy + gzz*uz*uz &
|
|
+(gxy*ux*uy + gxz*ux*uz + gyz*uy*uz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
ux = ux/fx
|
|
uy = uy/fx
|
|
uz = uz/fx
|
|
|
|
fx = gxx*vx*wx + gxy*vx*wy + gxz*vx*wz + &
|
|
gxy*vy*wx + gyy*vy*wy + gyz*vy*wz + &
|
|
gxz*vz*wx + gyz*vz*wy + gzz*vz*wz
|
|
fx = fx*f
|
|
wx = wx - fx*vx
|
|
wy = wy - fx*vy
|
|
wz = wz - fx*vz
|
|
fx = gxx*ux*wx + gxy*ux*wy + gxz*ux*wz + &
|
|
gxy*uy*wx + gyy*uy*wy + gyz*uy*wz + &
|
|
gxz*uz*wx + gyz*uz*wy + gzz*uz*wz
|
|
fx = fx*f
|
|
wx = wx - fx*ux
|
|
wy = wy - fx*uy
|
|
wz = wz - fx*uz
|
|
fx = gxx*wx*wx + gyy*wy*wy + gzz*wz*wz &
|
|
+(gxy*wx*wy + gxz*wx*wz + gyz*wy*wz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
wx = wx/fx
|
|
wy = wy/fx
|
|
wz = wz/fx
|
|
#elif (tetradtype == 1)
|
|
do i=1,ex(1)
|
|
do j=1,ex(2)
|
|
do k=1,ex(3)
|
|
if(abs(X(i)) < TINYRR .and. abs(Y(j)) < TINYRR .and. abs(Z(k)) < TINYRR)then
|
|
vx(i,j,k) = TINYRR
|
|
vy(i,j,k) = TINYRR
|
|
vz(i,j,k) = TINYRR
|
|
else
|
|
vx(i,j,k) = X(i)
|
|
vy(i,j,k) = Y(j)
|
|
vz(i,j,k) = Z(k)
|
|
endif
|
|
if(abs(X(i)) < TINYRR .and. abs(Y(j)) < TINYRR)then
|
|
ux(i,j,k) = - TINYRR
|
|
uy(i,j,k) = TINYRR
|
|
uz(i,j,k) = ZEO
|
|
wx(i,j,k) = TINYRR*Z(k)
|
|
wy(i,j,k) = TINYRR*Z(k)
|
|
wz(i,j,k) = -2*TINYRR*TINYRR
|
|
else
|
|
ux(i,j,k) = - Y(j)
|
|
uy(i,j,k) = X(i)
|
|
uz(i,j,k) = ZEO
|
|
wx(i,j,k) = X(i)*Z(k)
|
|
wy(i,j,k) = Y(j)*Z(k)
|
|
wz(i,j,k) = -(X(i)*X(i) + Y(j)*Y(j))
|
|
endif
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
f = 1.d0/chipn1
|
|
|
|
fx = gxx*wx*wx + gyy*wy*wy + gzz*wz*wz &
|
|
+(gxy*wx*wy + gxz*wx*wz + gyz*wy*wz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
wx = wx/fx
|
|
wy = wy/fx
|
|
wz = wz/fx
|
|
|
|
fx = gxx*wx*ux + gxy*wx*uy + gxz*wx*uz + &
|
|
gxy*wy*ux + gyy*wy*uy + gyz*wy*uz + &
|
|
gxz*wz*ux + gyz*wz*uy + gzz*wz*uz
|
|
fx = fx*f
|
|
ux = ux - fx*wx
|
|
uy = uy - fx*wy
|
|
uz = uz - fx*wz
|
|
fx = gxx*ux*ux + gyy*uy*uy + gzz*uz*uz &
|
|
+(gxy*ux*uy + gxz*ux*uz + gyz*uy*uz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
ux = ux/fx
|
|
uy = uy/fx
|
|
uz = uz/fx
|
|
|
|
fx = gxx*vx*wx + gxy*vx*wy + gxz*vx*wz + &
|
|
gxy*vy*wx + gyy*vy*wy + gyz*vy*wz + &
|
|
gxz*vz*wx + gyz*vz*wy + gzz*vz*wz
|
|
fx = fx*f
|
|
vx = vx - fx*wx
|
|
vy = vy - fx*wy
|
|
vz = vz - fx*wz
|
|
fx = gxx*ux*vx + gxy*ux*vy + gxz*ux*vz + &
|
|
gxy*uy*vx + gyy*uy*vy + gyz*uy*vz + &
|
|
gxz*uz*vx + gyz*uz*vy + gzz*uz*vz
|
|
fx = fx*f
|
|
vx = vx - fx*ux
|
|
vy = vy - fx*uy
|
|
vz = vz - fx*uz
|
|
fx = gxx*vx*vx + gyy*vy*vy + gzz*vz*vz &
|
|
+(gxy*vx*vy + gxz*vx*vz + gyz*vy*vz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
vx = vx/fx
|
|
vy = vy/fx
|
|
vz = vz/fx
|
|
#elif (tetradtype == 2)
|
|
do i=1,ex(1)
|
|
do j=1,ex(2)
|
|
do k=1,ex(3)
|
|
if(abs(X(i)) < TINYRR .and. abs(Y(j)) < TINYRR .and. abs(Z(k)) < TINYRR)then
|
|
vx(i,j,k) = TINYRR
|
|
vy(i,j,k) = TINYRR
|
|
vz(i,j,k) = TINYRR
|
|
else
|
|
vx(i,j,k) = X(i)
|
|
vy(i,j,k) = Y(j)
|
|
vz(i,j,k) = Z(k)
|
|
endif
|
|
if(abs(X(i)) < TINYRR .and. abs(Y(j)) < TINYRR)then
|
|
ux(i,j,k) = - TINYRR
|
|
uy(i,j,k) = TINYRR
|
|
uz(i,j,k) = ZEO
|
|
wx(i,j,k) = TINYRR*Z(k)
|
|
wy(i,j,k) = TINYRR*Z(k)
|
|
wz(i,j,k) = -2*TINYRR*TINYRR
|
|
else
|
|
ux(i,j,k) = - Y(j)
|
|
uy(i,j,k) = X(i)
|
|
uz(i,j,k) = ZEO
|
|
wx(i,j,k) = X(i)*Z(k)
|
|
wy(i,j,k) = Y(j)*Z(k)
|
|
wz(i,j,k) = -(X(i)*X(i) + Y(j)*Y(j))
|
|
endif
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
fx = vx
|
|
fy = vy
|
|
fz = vz
|
|
vx = gupxx*fx + gupxy*fy + gupxz*fz
|
|
vy = gupxy*fx + gupyy*fy + gupyz*fz
|
|
vz = gupxz*fx + gupyz*fy + gupzz*fz
|
|
|
|
f = 1.d0/chipn1
|
|
|
|
fx = gxx*vx*vx + gyy*vy*vy + gzz*vz*vz &
|
|
+(gxy*vx*vy + gxz*vx*vz + gyz*vy*vz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
vx = vx/fx
|
|
vy = vy/fx
|
|
vz = vz/fx
|
|
|
|
fx = gxx*vx*ux + gxy*vx*uy + gxz*vx*uz + &
|
|
gxy*vy*ux + gyy*vy*uy + gyz*vy*uz + &
|
|
gxz*vz*ux + gyz*vz*uy + gzz*vz*uz
|
|
fx = fx*f
|
|
ux = ux - fx*vx
|
|
uy = uy - fx*vy
|
|
uz = uz - fx*vz
|
|
fx = gxx*ux*ux + gyy*uy*uy + gzz*uz*uz &
|
|
+(gxy*ux*uy + gxz*ux*uz + gyz*uy*uz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
ux = ux/fx
|
|
uy = uy/fx
|
|
uz = uz/fx
|
|
|
|
fx = gxx*vx*wx + gxy*vx*wy + gxz*vx*wz + &
|
|
gxy*vy*wx + gyy*vy*wy + gyz*vy*wz + &
|
|
gxz*vz*wx + gyz*vz*wy + gzz*vz*wz
|
|
fx = fx*f
|
|
wx = wx - fx*vx
|
|
wy = wy - fx*vy
|
|
wz = wz - fx*vz
|
|
fx = gxx*ux*wx + gxy*ux*wy + gxz*ux*wz + &
|
|
gxy*uy*wx + gyy*uy*wy + gyz*uy*wz + &
|
|
gxz*uz*wx + gyz*uz*wy + gzz*uz*wz
|
|
fx = fx*f
|
|
wx = wx - fx*ux
|
|
wy = wy - fx*uy
|
|
wz = wz - fx*uz
|
|
fx = gxx*wx*wx + gyy*wy*wy + gzz*wz*wz &
|
|
+(gxy*wx*wy + gxz*wx*wz + gyz*wy*wz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
wx = wx/fx
|
|
wy = wy/fx
|
|
wz = wz/fx
|
|
#endif
|
|
|
|
call fderivs(ex,Axx,Axxx,Axxy,Axxz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0)
|
|
call fderivs(ex,Axy,Axyx,Axyy,Axyz,X,Y,Z,ANTI,ANTI,SYM ,Symmetry,0)
|
|
call fderivs(ex,Axz,Axzx,Axzy,Axzz,X,Y,Z,ANTI,SYM ,ANTI,Symmetry,0)
|
|
call fderivs(ex,Ayy,Ayyx,Ayyy,Ayyz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0)
|
|
call fderivs(ex,Ayz,Ayzx,Ayzy,Ayzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,0)
|
|
call fderivs(ex,Azz,Azzx,Azzy,Azzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0)
|
|
|
|
call fderivs(ex,chi,chix,chiy,chiz,X,Y,Z,SYM,SYM,SYM,symmetry,0)
|
|
call fderivs(ex,trK,fx,fy,fz,X,Y,Z,SYM,SYM,SYM,symmetry,0)
|
|
! compute D_k K_ij up to chi^-1
|
|
Axxx = Axxx - (Gamxxx*Axx + Gamyxx*Axy + Gamzxx*Axz)*TWO - chix/chipn1*Axx + F1o3*gxx*fx
|
|
Axxy = Axxy - (Gamxxy*Axx + Gamyxy*Axy + Gamzxy*Axz)*TWO - chiy/chipn1*Axx + F1o3*gxx*fy
|
|
Axxz = Axxz - (Gamxxz*Axx + Gamyxz*Axy + Gamzxz*Axz)*TWO - chiz/chipn1*Axx + F1o3*gxx*fz
|
|
Ayyx = Ayyx - (Gamxxy*Axy + Gamyxy*Ayy + Gamzxy*Ayz)*TWO - chix/chipn1*Ayy + F1o3*gyy*fx
|
|
Ayyy = Ayyy - (Gamxyy*Axy + Gamyyy*Ayy + Gamzyy*Ayz)*TWO - chiy/chipn1*Ayy + F1o3*gyy*fy
|
|
Ayyz = Ayyz - (Gamxyz*Axy + Gamyyz*Ayy + Gamzyz*Ayz)*TWO - chiz/chipn1*Ayy + F1o3*gyy*fz
|
|
Azzx = Azzx - (Gamxxz*Axz + Gamyxz*Ayz + Gamzxz*Azz)*TWO - chix/chipn1*Azz + F1o3*gzz*fx
|
|
Azzy = Azzy - (Gamxyz*Axz + Gamyyz*Ayz + Gamzyz*Azz)*TWO - chiy/chipn1*Azz + F1o3*gzz*fy
|
|
Azzz = Azzz - (Gamxzz*Axz + Gamyzz*Ayz + Gamzzz*Azz)*TWO - chiz/chipn1*Azz + F1o3*gzz*fz
|
|
Axyx = Axyx - (Gamxxy*Axx + Gamyxy*Axy + Gamzxy*Axz + &
|
|
Gamxxx*Axy + Gamyxx*Ayy + Gamzxx*Ayz) - chix/chipn1*Axy + F1o3*gxy*fx
|
|
Axyy = Axyy - (Gamxyy*Axx + Gamyyy*Axy + Gamzyy*Axz + &
|
|
Gamxxy*Axy + Gamyxy*Ayy + Gamzxy*Ayz) - chiy/chipn1*Axy + F1o3*gxy*fy
|
|
Axyz = Axyz - (Gamxyz*Axx + Gamyyz*Axy + Gamzyz*Axz + &
|
|
Gamxxz*Axy + Gamyxz*Ayy + Gamzxz*Ayz) - chiz/chipn1*Axy + F1o3*gxy*fz
|
|
Axzx = Axzx - (Gamxxz*Axx + Gamyxz*Axy + Gamzxz*Axz + &
|
|
Gamxxx*Axz + Gamyxx*Ayz + Gamzxx*Azz) - chix/chipn1*Axz + F1o3*gxz*fx
|
|
Axzy = Axzy - (Gamxyz*Axx + Gamyyz*Axy + Gamzyz*Axz + &
|
|
Gamxxy*Axz + Gamyxy*Ayz + Gamzxy*Azz) - chiy/chipn1*Axz + F1o3*gxz*fy
|
|
Axzz = Axzz - (Gamxzz*Axx + Gamyzz*Axy + Gamzzz*Axz + &
|
|
Gamxxz*Axz + Gamyxz*Ayz + Gamzxz*Azz) - chiz/chipn1*Axz + F1o3*gxz*fz
|
|
Ayzx = Ayzx - (Gamxxz*Axy + Gamyxz*Ayy + Gamzxz*Ayz + &
|
|
Gamxxy*Axz + Gamyxy*Ayz + Gamzxy*Azz) - chix/chipn1*Ayz + F1o3*gyz*fx
|
|
Ayzy = Ayzy - (Gamxyz*Axy + Gamyyz*Ayy + Gamzyz*Ayz + &
|
|
Gamxyy*Axz + Gamyyy*Ayz + Gamzyy*Azz) - chiy/chipn1*Ayz + F1o3*gyz*fy
|
|
Ayzz = Ayzz - (Gamxzz*Axy + Gamyzz*Ayy + Gamzzz*Ayz + &
|
|
Gamxyz*Axz + Gamyyz*Ayz + Gamzyz*Azz) - chiz/chipn1*Ayz + F1o3*gyz*fz
|
|
! symmetrize B_ij = v^k (D_k K_ij -D_j K_ik)
|
|
Bxx =(vy*(Axxy - Axyx) + vz*(Axxz - Axzx))*f
|
|
Byy =(vx*(Ayyx - Axyy) + vz*(Ayyz - Ayzy))*f
|
|
Bzz =(vx*(Azzx - Axzz) + vy*(Azzy - Ayzz))*f
|
|
Bxy =(vx*(Axyx - (Axxy+Axyx)/TWO) + vy*(Axyy-Ayyx)/TWO + vz*(Axyz - (Axzy+Ayzx)/TWO))*f
|
|
Bxz =(vx*(Axzx - (Axxz+Axzx)/TWO) + vy*(Axzy - (Axyz+Ayzx)/TWO) + vz*(Axzz-Azzx)/TWO)*f
|
|
Byz =(vx*(Ayzx - (Axyz+Axzy)/TWO) + vy*(Ayzy - (Ayyz+Ayzy)/TWO) + vz*(Ayzz-Azzy)/TWO)*f
|
|
! E_ij = R_ij - K_ik * K^k_j + K * K_ij
|
|
|
|
! K_ij up to chi^-1
|
|
Axxx = Axx + F1o3*trK*gxx
|
|
Axyx = Axy + F1o3*trK*gxy
|
|
Axzx = Axz + F1o3*trK*gxz
|
|
Ayyx = Ayy + F1o3*trK*gyy
|
|
Ayzx = Ayz + F1o3*trK*gyz
|
|
Azzx = Azz + F1o3*trK*gzz
|
|
! gup and A_ijk cancel a chi^-1
|
|
Exx = gupxx * Axxx * Axxx + gupyy * Axyx * Axyx + gupzz * Axzx * Axzx + &
|
|
TWO * (gupxy * Axxx * Axyx + gupxz * Axxx * Axzx + gupyz * Axyx * Axzx)
|
|
Eyy = gupxx * Axyx * Axyx + gupyy * Ayyx * Ayyx + gupzz * Ayzx * Ayzx + &
|
|
TWO * (gupxy * Axyx * Ayyx + gupxz * Axyx * Ayzx + gupyz * Ayyx * Ayzx)
|
|
Ezz = gupxx * Axzx * Axzx + gupyy * Ayzx * Ayzx + gupzz * Azzx * Azzx + &
|
|
TWO * (gupxy * Axzx * Ayzx + gupxz * Axzx * Azzx + gupyz * Ayzx * Azzx)
|
|
Exy = gupxx * Axxx * Axyx + gupyy * Axyx * Ayyx + gupzz * Axzx * Ayzx + &
|
|
gupxy *(Axxx * Ayyx + Axyx * Axyx) + &
|
|
gupxz *(Axxx * Ayzx + Axzx * Axyx) + &
|
|
gupyz *(Axyx * Ayzx + Axzx * Ayyx)
|
|
Exz = gupxx * Axxx * Axzx + gupyy * Axyx * Ayzx + gupzz * Axzx * Azzx + &
|
|
gupxy *(Axxx * Ayzx + Axyx * Axzx) + &
|
|
gupxz *(Axxx * Azzx + Axzx * Axzx) + &
|
|
gupyz *(Axyx * Azzx + Axzx * Ayzx)
|
|
Eyz = gupxx * Axyx * Axzx + gupyy * Ayyx * Ayzx + gupzz * Ayzx * Azzx + &
|
|
gupxy *(Axyx * Ayzx + Ayyx * Axzx) + &
|
|
gupxz *(Axyx * Azzx + Ayzx * Axzx) + &
|
|
gupyz *(Ayyx * Azzx + Ayzx * Ayzx)
|
|
|
|
Exx = Rxx - (Exx - Axxx*trK)*f - Bxx
|
|
Exy = Rxy - (Exy - Axyx*trK)*f - Bxy
|
|
Exz = Rxz - (Exz - Axzx*trK)*f - Bxz
|
|
Eyy = Ryy - (Eyy - Ayyx*trK)*f - Byy
|
|
Eyz = Ryz - (Eyz - Ayzx*trK)*f - Byz
|
|
Ezz = Rzz - (Ezz - Azzx*trK)*f - Bzz
|
|
!set m = (u - iw)/sqrt(2) following Frans, PRD 75, 124018(2007)
|
|
! compute uuww^ij = u^i * u^j - w^i * w^j
|
|
uuwwxx = ux * ux - wx * wx
|
|
uuwwxy = ux * uy - wx * wy
|
|
uuwwxz = ux * uz - wx * wz
|
|
uuwwyy = uy * uy - wy * wy
|
|
uuwwyz = uy * uz - wy * wz
|
|
uuwwzz = uz * uz - wz * wz
|
|
|
|
! compute uw^ij = u^i * w^j + w^i * u^j
|
|
uwxx = ux * wx + wx * ux
|
|
uwxy = ux * wy + wx * uy
|
|
uwxz = ux * wz + wx * uz
|
|
uwyy = uy * wy + wy * uy
|
|
uwyz = uy * wz + wy * uz
|
|
uwzz = uz * wz + wz * uz
|
|
!the real part of Psi4
|
|
Rpsi4 = Exx * uuwwxx + Eyy * uuwwyy + Ezz * uuwwzz &
|
|
+ (Exy * uuwwxy + Exz * uuwwxz + Eyz * uuwwyz) * TWO
|
|
|
|
!the imaginary part of Psi4
|
|
Ipsi4 = Exx * uwxx + Eyy * uwyy + Ezz * uwzz &
|
|
+ (Exy * uwxy + Exz * uwxz + Eyz * uwyz) * TWO
|
|
|
|
!multiply with -1/2
|
|
Rpsi4 = - Rpsi4/TWO
|
|
Ipsi4 = - Ipsi4/TWO
|
|
|
|
return
|
|
|
|
end subroutine getnp4
|
|
!-----------------------------------------------------------------------------
|
|
!
|
|
! compute the Newman-Penrose Weyl scalar Psi4
|
|
! for BSSN dynamical variables for shell
|
|
!
|
|
!-----------------------------------------------------------------------------
|
|
|
|
subroutine getnp4_ss(ex,crho,sigma,R, X, Y, Z, &
|
|
drhodx, drhody, drhodz, &
|
|
dsigmadx,dsigmady,dsigmadz, &
|
|
dRdx,dRdy,dRdz, &
|
|
drhodxx,drhodxy,drhodxz,drhodyy,drhodyz,drhodzz, &
|
|
dsigmadxx,dsigmadxy,dsigmadxz,dsigmadyy,dsigmadyz,dsigmadzz, &
|
|
dRdxx,dRdxy,dRdxz,dRdyy,dRdyz,dRdzz, &
|
|
chi, trK, &
|
|
dxx,gxy,gxz,dyy,gyz,dzz, &
|
|
Axx,Axy,Axz,Ayy,Ayz,Azz, &
|
|
Gamxxx,Gamxxy,Gamxxz,Gamxyy,Gamxyz,Gamxzz,&
|
|
Gamyxx,Gamyxy,Gamyxz,Gamyyy,Gamyyz,Gamyzz,&
|
|
Gamzxx,Gamzxy,Gamzxz,Gamzyy,Gamzyz,Gamzzz,&
|
|
Rxx,Rxy,Rxz,Ryy,Ryz,Rzz,&
|
|
Rpsi4, Ipsi4, &
|
|
symmetry,sst)
|
|
|
|
implicit none
|
|
|
|
!~~~~~~> Input parameters:
|
|
|
|
integer,intent(in ):: ex(1:3),symmetry,sst
|
|
double precision,intent(in),dimension(ex(1))::crho
|
|
double precision,intent(in),dimension(ex(2))::sigma
|
|
double precision,intent(in),dimension(ex(3))::R
|
|
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: X,Y,Z
|
|
double precision,intent(in),dimension(ex(1),ex(2),ex(3))::drhodx, drhody, drhodz
|
|
double precision,intent(in),dimension(ex(1),ex(2),ex(3))::dsigmadx,dsigmady,dsigmadz
|
|
double precision,intent(in),dimension(ex(1),ex(2),ex(3))::dRdx,dRdy,dRdz
|
|
double precision,intent(in),dimension(ex(1),ex(2),ex(3))::drhodxx,drhodxy,drhodxz,drhodyy,drhodyz,drhodzz
|
|
double precision,intent(in),dimension(ex(1),ex(2),ex(3))::dsigmadxx,dsigmadxy,dsigmadxz,dsigmadyy,dsigmadyz,dsigmadzz
|
|
double precision,intent(in),dimension(ex(1),ex(2),ex(3))::dRdxx,dRdxy,dRdxz,dRdyy,dRdyz,dRdzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: dxx,gxy,gxz,dyy,gyz,dzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: Axx,Axy,Axz,Ayy,Ayz,Azz
|
|
real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: chi,trK
|
|
! physical second kind of connection
|
|
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Gamxxx, Gamxxy, Gamxxz
|
|
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Gamxyy, Gamxyz, Gamxzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Gamyxx, Gamyxy, Gamyxz
|
|
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Gamyyy, Gamyyz, Gamyzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Gamzxx, Gamzxy, Gamzxz
|
|
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Gamzyy, Gamzyz, Gamzzz
|
|
! physical Ricci tensor
|
|
real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Rxx,Rxy,Rxz,Ryy,Ryz,Rzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)), intent(out):: Rpsi4,Ipsi4
|
|
|
|
!~~~~~~> Other variables:
|
|
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: f,fx,fy,fz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: gxx,gyy,gzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: chix,chiy,chiz,chipn1
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: vx,vy,vz,ux,uy,uz,wx,wy,wz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: Exx,Exy,Exz,Eyy,Eyz,Ezz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: Bxx,Bxy,Bxz,Byy,Byz,Bzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: Axxx,Axxy,Axxz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: Axyx,Axyy,Axyz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: Axzx,Axzy,Axzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: Ayyx,Ayyy,Ayyz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: Ayzx,Ayzy,Ayzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: Azzx,Azzy,Azzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: gupxx,gupxy,gupxz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: gupyy,gupyz,gupzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: uuwwxx,uuwwxy,uuwwxz,uuwwyy,uuwwyz,uuwwzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: uwxx,uwxy,uwxz,uwyy,uwyz,uwzz
|
|
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: adm_dxx,adm_gxy,adm_gxz,adm_dyy,adm_gyz,adm_dzz
|
|
real*8, dimension(ex(1),ex(2),ex(3)) :: Kxx,Kxy,Kxz,Kyy,Kyz,Kzz
|
|
|
|
real*8, parameter :: ZEO = 0.d0, ONE = 1.d0, TWO = 2.d0
|
|
real*8, parameter :: F1o3 = 1.d0/3.d0, FOUR = 4.d0
|
|
real*8, parameter :: SYM = 1.D0, ANTI= - 1.D0
|
|
integer::i,j,k
|
|
real*8,parameter::TINYRR=1.d-14
|
|
|
|
gxx = dxx + ONE
|
|
gyy = dyy + ONE
|
|
gzz = dzz + ONE
|
|
chipn1= chi + ONE
|
|
|
|
#if (ABV == 1)
|
|
call bssn2adm(ex,chipn1,trK,gxx,gxy,gxz,gyy,gyz,gzz, &
|
|
Axx,Axy,Axz,Ayy,Ayz,Azz, &
|
|
adm_dxx,adm_gxy,adm_gxz,adm_dyy,adm_gyz,adm_dzz, &
|
|
Kxx,Kxy,Kxz,Kyy,Kyz,Kzz)
|
|
adm_dxx = adm_dxx - ONE
|
|
adm_dyy = adm_dyy - ONE
|
|
adm_dzz = adm_dzz - ONE
|
|
call adm_ricci_gamma_ss(ex,crho,sigma,R,X, Y, Z, &
|
|
drhodx, drhody, drhodz, &
|
|
dsigmadx,dsigmady,dsigmadz, &
|
|
dRdx,dRdy,dRdz, &
|
|
drhodxx,drhodxy,drhodxz,drhodyy,drhodyz,drhodzz, &
|
|
dsigmadxx,dsigmadxy,dsigmadxz,dsigmadyy,dsigmadyz,dsigmadzz, &
|
|
dRdxx,dRdxy,dRdxz,dRdyy,dRdyz,dRdzz, &
|
|
adm_dxx,adm_gxy,adm_gxz,adm_dyy,adm_gyz,adm_dzz,&
|
|
Gamxxx,Gamxxy,Gamxxz,Gamxyy,Gamxyz,Gamxzz,&
|
|
Gamyxx,Gamyxy,Gamyxz,Gamyyy,Gamyyz,Gamyzz,&
|
|
Gamzxx,Gamzxy,Gamzxz,Gamzyy,Gamzyz,Gamzzz,&
|
|
Rxx,Rxy,Rxz,Ryy,Ryz,Rzz,&
|
|
Symmetry,0,sst)
|
|
#endif
|
|
|
|
! invert tilted metric
|
|
gupzz = gxx * gyy * gzz + gxy * gyz * gxz + gxz * gxy * gyz - &
|
|
gxz * gyy * gxz - gxy * gxy * gzz - gxx * gyz * gyz
|
|
gupxx = ( gyy * gzz - gyz * gyz ) / gupzz
|
|
gupxy = - ( gxy * gzz - gyz * gxz ) / gupzz
|
|
gupxz = ( gxy * gyz - gyy * gxz ) / gupzz
|
|
gupyy = ( gxx * gzz - gxz * gxz ) / gupzz
|
|
gupyz = - ( gxx * gyz - gxy * gxz ) / gupzz
|
|
gupzz = ( gxx * gyy - gxy * gxy ) / gupzz
|
|
|
|
! initialize U, V, W vetors
|
|
#if (tetradtype == 0)
|
|
do i=1,ex(1)
|
|
do j=1,ex(2)
|
|
do k=1,ex(3)
|
|
if(abs(X(i,j,k)) < TINYRR .and. abs(Y(i,j,k)) < TINYRR .and. abs(Z(i,j,k)) < TINYRR)then
|
|
vx(i,j,k) = TINYRR
|
|
vy(i,j,k) = TINYRR
|
|
vz(i,j,k) = TINYRR
|
|
else
|
|
vx(i,j,k) = X(i,j,k)
|
|
vy(i,j,k) = Y(i,j,k)
|
|
vz(i,j,k) = Z(i,j,k)
|
|
endif
|
|
if(abs(X(i,j,k)) < TINYRR .and. abs(Y(i,j,k)) < TINYRR)then
|
|
ux(i,j,k) = - TINYRR
|
|
uy(i,j,k) = TINYRR
|
|
uz(i,j,k) = ZEO
|
|
wx(i,j,k) = TINYRR*Z(i,j,k)
|
|
wy(i,j,k) = TINYRR*Z(i,j,k)
|
|
wz(i,j,k) = -2*TINYRR*TINYRR
|
|
else
|
|
ux(i,j,k) = - Y(i,j,k)
|
|
uy(i,j,k) = X(i,j,k)
|
|
uz(i,j,k) = ZEO
|
|
wx(i,j,k) = X(i,j,k)*Z(i,j,k)
|
|
wy(i,j,k) = Y(i,j,k)*Z(i,j,k)
|
|
wz(i,j,k) = -(X(i,j,k)*X(i,j,k) + Y(i,j,k)*Y(i,j,k))
|
|
endif
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
f = 1.d0/chipn1
|
|
|
|
fx = gxx*vx*vx + gyy*vy*vy + gzz*vz*vz &
|
|
+(gxy*vx*vy + gxz*vx*vz + gyz*vy*vz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
vx = vx/fx
|
|
vy = vy/fx
|
|
vz = vz/fx
|
|
|
|
fx = gxx*vx*ux + gxy*vx*uy + gxz*vx*uz + &
|
|
gxy*vy*ux + gyy*vy*uy + gyz*vy*uz + &
|
|
gxz*vz*ux + gyz*vz*uy + gzz*vz*uz
|
|
fx = fx*f
|
|
ux = ux - fx*vx
|
|
uy = uy - fx*vy
|
|
uz = uz - fx*vz
|
|
fx = gxx*ux*ux + gyy*uy*uy + gzz*uz*uz &
|
|
+(gxy*ux*uy + gxz*ux*uz + gyz*uy*uz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
ux = ux/fx
|
|
uy = uy/fx
|
|
uz = uz/fx
|
|
|
|
fx = gxx*vx*wx + gxy*vx*wy + gxz*vx*wz + &
|
|
gxy*vy*wx + gyy*vy*wy + gyz*vy*wz + &
|
|
gxz*vz*wx + gyz*vz*wy + gzz*vz*wz
|
|
fx = fx*f
|
|
wx = wx - fx*vx
|
|
wy = wy - fx*vy
|
|
wz = wz - fx*vz
|
|
fx = gxx*ux*wx + gxy*ux*wy + gxz*ux*wz + &
|
|
gxy*uy*wx + gyy*uy*wy + gyz*uy*wz + &
|
|
gxz*uz*wx + gyz*uz*wy + gzz*uz*wz
|
|
fx = fx*f
|
|
wx = wx - fx*ux
|
|
wy = wy - fx*uy
|
|
wz = wz - fx*uz
|
|
fx = gxx*wx*wx + gyy*wy*wy + gzz*wz*wz &
|
|
+(gxy*wx*wy + gxz*wx*wz + gyz*wy*wz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
wx = wx/fx
|
|
wy = wy/fx
|
|
wz = wz/fx
|
|
#elif (tetradtype == 1)
|
|
do i=1,ex(1)
|
|
do j=1,ex(2)
|
|
do k=1,ex(3)
|
|
if(abs(X(i,j,k)) < TINYRR .and. abs(Y(i,j,k)) < TINYRR .and. abs(Z(i,j,k)) < TINYRR)then
|
|
vx(i,j,k) = TINYRR
|
|
vy(i,j,k) = TINYRR
|
|
vz(i,j,k) = TINYRR
|
|
else
|
|
vx(i,j,k) = X(i,j,k)
|
|
vy(i,j,k) = Y(i,j,k)
|
|
vz(i,j,k) = Z(i,j,k)
|
|
endif
|
|
if(abs(X(i,j,k)) < TINYRR .and. abs(Y(i,j,k)) < TINYRR)then
|
|
ux(i,j,k) = - TINYRR
|
|
uy(i,j,k) = TINYRR
|
|
uz(i,j,k) = ZEO
|
|
wx(i,j,k) = TINYRR*Z(i,j,k)
|
|
wy(i,j,k) = TINYRR*Z(i,j,k)
|
|
wz(i,j,k) = -2*TINYRR*TINYRR
|
|
else
|
|
ux(i,j,k) = - Y(i,j,k)
|
|
uy(i,j,k) = X(i,j,k)
|
|
uz(i,j,k) = ZEO
|
|
wx(i,j,k) = X(i,j,k)*Z(i,j,k)
|
|
wy(i,j,k) = Y(i,j,k)*Z(i,j,k)
|
|
wz(i,j,k) = -(X(i,j,k)*X(i,j,k) + Y(i,j,k)*Y(i,j,k))
|
|
endif
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
f = 1.d0/chipn1
|
|
|
|
fx = gxx*wx*wx + gyy*wy*wy + gzz*wz*wz &
|
|
+(gxy*wx*wy + gxz*wx*wz + gyz*wy*wz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
wx = wx/fx
|
|
wy = wy/fx
|
|
wz = wz/fx
|
|
|
|
fx = gxx*wx*ux + gxy*wx*uy + gxz*wx*uz + &
|
|
gxy*wy*ux + gyy*wy*uy + gyz*wy*uz + &
|
|
gxz*wz*ux + gyz*wz*uy + gzz*wz*uz
|
|
fx = fx*f
|
|
ux = ux - fx*wx
|
|
uy = uy - fx*wy
|
|
uz = uz - fx*wz
|
|
fx = gxx*ux*ux + gyy*uy*uy + gzz*uz*uz &
|
|
+(gxy*ux*uy + gxz*ux*uz + gyz*uy*uz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
ux = ux/fx
|
|
uy = uy/fx
|
|
uz = uz/fx
|
|
|
|
fx = gxx*vx*wx + gxy*vx*wy + gxz*vx*wz + &
|
|
gxy*vy*wx + gyy*vy*wy + gyz*vy*wz + &
|
|
gxz*vz*wx + gyz*vz*wy + gzz*vz*wz
|
|
fx = fx*f
|
|
vx = vx - fx*wx
|
|
vy = vy - fx*wy
|
|
vz = vz - fx*wz
|
|
fx = gxx*ux*vx + gxy*ux*vy + gxz*ux*vz + &
|
|
gxy*uy*vx + gyy*uy*vy + gyz*uy*vz + &
|
|
gxz*uz*vx + gyz*uz*vy + gzz*uz*vz
|
|
fx = fx*f
|
|
vx = vx - fx*ux
|
|
vy = vy - fx*uy
|
|
vz = vz - fx*uz
|
|
fx = gxx*vx*vx + gyy*vy*vy + gzz*vz*vz &
|
|
+(gxy*vx*vy + gxz*vx*vz + gyz*vy*vz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
vx = vx/fx
|
|
vy = vy/fx
|
|
vz = vz/fx
|
|
#elif (tetradtype == 2)
|
|
do i=1,ex(1)
|
|
do j=1,ex(2)
|
|
do k=1,ex(3)
|
|
if(abs(X(i,j,k)) < TINYRR .and. abs(Y(i,j,k)) < TINYRR .and. abs(Z(i,j,k)) < TINYRR)then
|
|
vx(i,j,k) = TINYRR
|
|
vy(i,j,k) = TINYRR
|
|
vz(i,j,k) = TINYRR
|
|
else
|
|
vx(i,j,k) = X(i,j,k)
|
|
vy(i,j,k) = Y(i,j,k)
|
|
vz(i,j,k) = Z(i,j,k)
|
|
endif
|
|
if(abs(X(i,j,k)) < TINYRR .and. abs(Y(i,j,k)) < TINYRR)then
|
|
ux(i,j,k) = - TINYRR
|
|
uy(i,j,k) = TINYRR
|
|
uz(i,j,k) = ZEO
|
|
wx(i,j,k) = TINYRR*Z(i,j,k)
|
|
wy(i,j,k) = TINYRR*Z(i,j,k)
|
|
wz(i,j,k) = -2*TINYRR*TINYRR
|
|
else
|
|
ux(i,j,k) = - Y(i,j,k)
|
|
uy(i,j,k) = X(i,j,k)
|
|
uz(i,j,k) = ZEO
|
|
wx(i,j,k) = X(i,j,k)*Z(i,j,k)
|
|
wy(i,j,k) = Y(i,j,k)*Z(i,j,k)
|
|
wz(i,j,k) = -(X(i,j,k)*X(i,j,k) + Y(i,j,k)*Y(i,j,k))
|
|
endif
|
|
enddo
|
|
enddo
|
|
enddo
|
|
|
|
fx = vx
|
|
fy = vy
|
|
fz = vz
|
|
vx = gupxx*fx + gupxy*fy + gupxz*fz
|
|
vy = gupxy*fx + gupyy*fy + gupyz*fz
|
|
vz = gupxz*fx + gupyz*fy + gupzz*fz
|
|
|
|
f = 1.d0/chipn1
|
|
|
|
fx = gxx*vx*vx + gyy*vy*vy + gzz*vz*vz &
|
|
+(gxy*vx*vy + gxz*vx*vz + gyz*vy*vz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
vx = vx/fx
|
|
vy = vy/fx
|
|
vz = vz/fx
|
|
|
|
fx = gxx*vx*ux + gxy*vx*uy + gxz*vx*uz + &
|
|
gxy*vy*ux + gyy*vy*uy + gyz*vy*uz + &
|
|
gxz*vz*ux + gyz*vz*uy + gzz*vz*uz
|
|
fx = fx*f
|
|
ux = ux - fx*vx
|
|
uy = uy - fx*vy
|
|
uz = uz - fx*vz
|
|
fx = gxx*ux*ux + gyy*uy*uy + gzz*uz*uz &
|
|
+(gxy*ux*uy + gxz*ux*uz + gyz*uy*uz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
ux = ux/fx
|
|
uy = uy/fx
|
|
uz = uz/fx
|
|
|
|
fx = gxx*vx*wx + gxy*vx*wy + gxz*vx*wz + &
|
|
gxy*vy*wx + gyy*vy*wy + gyz*vy*wz + &
|
|
gxz*vz*wx + gyz*vz*wy + gzz*vz*wz
|
|
fx = fx*f
|
|
wx = wx - fx*vx
|
|
wy = wy - fx*vy
|
|
wz = wz - fx*vz
|
|
fx = gxx*ux*wx + gxy*ux*wy + gxz*ux*wz + &
|
|
gxy*uy*wx + gyy*uy*wy + gyz*uy*wz + &
|
|
gxz*uz*wx + gyz*uz*wy + gzz*uz*wz
|
|
fx = fx*f
|
|
wx = wx - fx*ux
|
|
wy = wy - fx*uy
|
|
wz = wz - fx*uz
|
|
fx = gxx*wx*wx + gyy*wy*wy + gzz*wz*wz &
|
|
+(gxy*wx*wy + gxz*wx*wz + gyz*wy*wz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
wx = wx/fx
|
|
wy = wy/fx
|
|
wz = wz/fx
|
|
#endif
|
|
|
|
call fderivs_shc(ex,Axx,Axxx,Axxy,Axxz,crho,sigma,R, SYM, SYM,SYM,Symmetry,0,sst, &
|
|
drhodx, drhody, drhodz, &
|
|
dsigmadx,dsigmady,dsigmadz, &
|
|
dRdx,dRdy,dRdz)
|
|
call fderivs_shc(ex,Axy,Axyx,Axyy,Axyz,crho,sigma,R,ANTI,ANTI,SYM,Symmetry,0,sst, &
|
|
drhodx, drhody, drhodz, &
|
|
dsigmadx,dsigmady,dsigmadz, &
|
|
dRdx,dRdy,dRdz)
|
|
call fderivs_shc(ex,Axz,Axzx,Axzy,Axzz,crho,sigma,R,ANTI,SYM ,ANTI,Symmetry,0,sst, &
|
|
drhodx, drhody, drhodz, &
|
|
dsigmadx,dsigmady,dsigmadz, &
|
|
dRdx,dRdy,dRdz)
|
|
call fderivs_shc(ex,Ayy,Ayyx,Ayyy,Ayyz,crho,sigma,R, SYM, SYM,SYM,Symmetry,0,sst, &
|
|
drhodx, drhody, drhodz, &
|
|
dsigmadx,dsigmady,dsigmadz, &
|
|
dRdx,dRdy,dRdz)
|
|
call fderivs_shc(ex,Ayz,Ayzx,Ayzy,Ayzz,crho,sigma,R,SYM ,ANTI,ANTI,Symmetry,0,sst, &
|
|
drhodx, drhody, drhodz, &
|
|
dsigmadx,dsigmady,dsigmadz, &
|
|
dRdx,dRdy,dRdz)
|
|
call fderivs_shc(ex,Azz,Azzx,Azzy,Azzz,crho,sigma,R, SYM, SYM,SYM,Symmetry,0,sst, &
|
|
drhodx, drhody, drhodz, &
|
|
dsigmadx,dsigmady,dsigmadz, &
|
|
dRdx,dRdy,dRdz)
|
|
|
|
call fderivs_shc(ex,chi,chix,chiy,chiz,crho,sigma,R, SYM, SYM,SYM,Symmetry,0,sst, &
|
|
drhodx, drhody, drhodz, &
|
|
dsigmadx,dsigmady,dsigmadz, &
|
|
dRdx,dRdy,dRdz)
|
|
|
|
call fderivs_shc(ex,trK,fx,fy,fz,crho,sigma,R, SYM, SYM,SYM,Symmetry,0,sst, &
|
|
drhodx, drhody, drhodz, &
|
|
dsigmadx,dsigmady,dsigmadz, &
|
|
dRdx,dRdy,dRdz)
|
|
! compute D_k K_ij up to chi^-1
|
|
Axxx = Axxx - (Gamxxx*Axx + Gamyxx*Axy + Gamzxx*Axz)*TWO - chix/chipn1*Axx + F1o3*gxx*fx
|
|
Axxy = Axxy - (Gamxxy*Axx + Gamyxy*Axy + Gamzxy*Axz)*TWO - chiy/chipn1*Axx + F1o3*gxx*fy
|
|
Axxz = Axxz - (Gamxxz*Axx + Gamyxz*Axy + Gamzxz*Axz)*TWO - chiz/chipn1*Axx + F1o3*gxx*fz
|
|
Ayyx = Ayyx - (Gamxxy*Axy + Gamyxy*Ayy + Gamzxy*Ayz)*TWO - chix/chipn1*Ayy + F1o3*gyy*fx
|
|
Ayyy = Ayyy - (Gamxyy*Axy + Gamyyy*Ayy + Gamzyy*Ayz)*TWO - chiy/chipn1*Ayy + F1o3*gyy*fy
|
|
Ayyz = Ayyz - (Gamxyz*Axy + Gamyyz*Ayy + Gamzyz*Ayz)*TWO - chiz/chipn1*Ayy + F1o3*gyy*fz
|
|
Azzx = Azzx - (Gamxxz*Axz + Gamyxz*Ayz + Gamzxz*Azz)*TWO - chix/chipn1*Azz + F1o3*gzz*fx
|
|
Azzy = Azzy - (Gamxyz*Axz + Gamyyz*Ayz + Gamzyz*Azz)*TWO - chiy/chipn1*Azz + F1o3*gzz*fy
|
|
Azzz = Azzz - (Gamxzz*Axz + Gamyzz*Ayz + Gamzzz*Azz)*TWO - chiz/chipn1*Azz + F1o3*gzz*fz
|
|
Axyx = Axyx - (Gamxxy*Axx + Gamyxy*Axy + Gamzxy*Axz + &
|
|
Gamxxx*Axy + Gamyxx*Ayy + Gamzxx*Ayz) - chix/chipn1*Axy + F1o3*gxy*fx
|
|
Axyy = Axyy - (Gamxyy*Axx + Gamyyy*Axy + Gamzyy*Axz + &
|
|
Gamxxy*Axy + Gamyxy*Ayy + Gamzxy*Ayz) - chiy/chipn1*Axy + F1o3*gxy*fy
|
|
Axyz = Axyz - (Gamxyz*Axx + Gamyyz*Axy + Gamzyz*Axz + &
|
|
Gamxxz*Axy + Gamyxz*Ayy + Gamzxz*Ayz) - chiz/chipn1*Axy + F1o3*gxy*fz
|
|
Axzx = Axzx - (Gamxxz*Axx + Gamyxz*Axy + Gamzxz*Axz + &
|
|
Gamxxx*Axz + Gamyxx*Ayz + Gamzxx*Azz) - chix/chipn1*Axz + F1o3*gxz*fx
|
|
Axzy = Axzy - (Gamxyz*Axx + Gamyyz*Axy + Gamzyz*Axz + &
|
|
Gamxxy*Axz + Gamyxy*Ayz + Gamzxy*Azz) - chiy/chipn1*Axz + F1o3*gxz*fy
|
|
Axzz = Axzz - (Gamxzz*Axx + Gamyzz*Axy + Gamzzz*Axz + &
|
|
Gamxxz*Axz + Gamyxz*Ayz + Gamzxz*Azz) - chiz/chipn1*Axz + F1o3*gxz*fz
|
|
Ayzx = Ayzx - (Gamxxz*Axy + Gamyxz*Ayy + Gamzxz*Ayz + &
|
|
Gamxxy*Axz + Gamyxy*Ayz + Gamzxy*Azz) - chix/chipn1*Ayz + F1o3*gyz*fx
|
|
Ayzy = Ayzy - (Gamxyz*Axy + Gamyyz*Ayy + Gamzyz*Ayz + &
|
|
Gamxyy*Axz + Gamyyy*Ayz + Gamzyy*Azz) - chiy/chipn1*Ayz + F1o3*gyz*fy
|
|
Ayzz = Ayzz - (Gamxzz*Axy + Gamyzz*Ayy + Gamzzz*Ayz + &
|
|
Gamxyz*Axz + Gamyyz*Ayz + Gamzyz*Azz) - chiz/chipn1*Ayz + F1o3*gyz*fz
|
|
! symmetrize B_ij = v^k (D_k K_ij -D_j K_ik)
|
|
Bxx =(vy*(Axxy - Axyx) + vz*(Axxz - Axzx))*f
|
|
Byy =(vx*(Ayyx - Axyy) + vz*(Ayyz - Ayzy))*f
|
|
Bzz =(vx*(Azzx - Axzz) + vy*(Azzy - Ayzz))*f
|
|
Bxy =(vx*(Axyx - (Axxy+Axyx)/TWO) + vy*(Axyy-Ayyx)/TWO + vz*(Axyz - (Axzy+Ayzx)/TWO))*f
|
|
Bxz =(vx*(Axzx - (Axxz+Axzx)/TWO) + vy*(Axzy - (Axyz+Ayzx)/TWO) + vz*(Axzz-Azzx)/TWO)*f
|
|
Byz =(vx*(Ayzx - (Axyz+Axzy)/TWO) + vy*(Ayzy - (Ayyz+Ayzy)/TWO) + vz*(Ayzz-Azzy)/TWO)*f
|
|
! E_ij = R_ij - K_ik * K^k_j + K * K_ij
|
|
|
|
! K_ij up to chi^-1
|
|
Axxx = Axx + F1o3*trK*gxx
|
|
Axyx = Axy + F1o3*trK*gxy
|
|
Axzx = Axz + F1o3*trK*gxz
|
|
Ayyx = Ayy + F1o3*trK*gyy
|
|
Ayzx = Ayz + F1o3*trK*gyz
|
|
Azzx = Azz + F1o3*trK*gzz
|
|
! gup and A_ijk cancel a chi^-1
|
|
Exx = gupxx * Axxx * Axxx + gupyy * Axyx * Axyx + gupzz * Axzx * Axzx + &
|
|
TWO * (gupxy * Axxx * Axyx + gupxz * Axxx * Axzx + gupyz * Axyx * Axzx)
|
|
Eyy = gupxx * Axyx * Axyx + gupyy * Ayyx * Ayyx + gupzz * Ayzx * Ayzx + &
|
|
TWO * (gupxy * Axyx * Ayyx + gupxz * Axyx * Ayzx + gupyz * Ayyx * Ayzx)
|
|
Ezz = gupxx * Axzx * Axzx + gupyy * Ayzx * Ayzx + gupzz * Azzx * Azzx + &
|
|
TWO * (gupxy * Axzx * Ayzx + gupxz * Axzx * Azzx + gupyz * Ayzx * Azzx)
|
|
Exy = gupxx * Axxx * Axyx + gupyy * Axyx * Ayyx + gupzz * Axzx * Ayzx + &
|
|
gupxy *(Axxx * Ayyx + Axyx * Axyx) + &
|
|
gupxz *(Axxx * Ayzx + Axzx * Axyx) + &
|
|
gupyz *(Axyx * Ayzx + Axzx * Ayyx)
|
|
Exz = gupxx * Axxx * Axzx + gupyy * Axyx * Ayzx + gupzz * Axzx * Azzx + &
|
|
gupxy *(Axxx * Ayzx + Axyx * Axzx) + &
|
|
gupxz *(Axxx * Azzx + Axzx * Axzx) + &
|
|
gupyz *(Axyx * Azzx + Axzx * Ayzx)
|
|
Eyz = gupxx * Axyx * Axzx + gupyy * Ayyx * Ayzx + gupzz * Ayzx * Azzx + &
|
|
gupxy *(Axyx * Ayzx + Ayyx * Axzx) + &
|
|
gupxz *(Axyx * Azzx + Ayzx * Axzx) + &
|
|
gupyz *(Ayyx * Azzx + Ayzx * Ayzx)
|
|
|
|
Exx = Rxx - (Exx - Axxx*trK)*f - Bxx
|
|
Exy = Rxy - (Exy - Axyx*trK)*f - Bxy
|
|
Exz = Rxz - (Exz - Axzx*trK)*f - Bxz
|
|
Eyy = Ryy - (Eyy - Ayyx*trK)*f - Byy
|
|
Eyz = Ryz - (Eyz - Ayzx*trK)*f - Byz
|
|
Ezz = Rzz - (Ezz - Azzx*trK)*f - Bzz
|
|
!set m = (u - iw)/sqrt(2) following Frans, PRD 75, 124018(2007)
|
|
! compute uuww^ij = u^i * u^j - w^i * w^j
|
|
uuwwxx = ux * ux - wx * wx
|
|
uuwwxy = ux * uy - wx * wy
|
|
uuwwxz = ux * uz - wx * wz
|
|
uuwwyy = uy * uy - wy * wy
|
|
uuwwyz = uy * uz - wy * wz
|
|
uuwwzz = uz * uz - wz * wz
|
|
|
|
! compute uw^ij = u^i * w^j + w^i * u^j
|
|
uwxx = ux * wx + wx * ux
|
|
uwxy = ux * wy + wx * uy
|
|
uwxz = ux * wz + wx * uz
|
|
uwyy = uy * wy + wy * uy
|
|
uwyz = uy * wz + wy * uz
|
|
uwzz = uz * wz + wz * uz
|
|
!the real part of Psi4
|
|
Rpsi4 = Exx * uuwwxx + Eyy * uuwwyy + Ezz * uuwwzz &
|
|
+ (Exy * uuwwxy + Exz * uuwwxz + Eyz * uuwwyz) * TWO
|
|
|
|
!the imaginary part of Psi4
|
|
Ipsi4 = Exx * uwxx + Eyy * uwyy + Ezz * uwzz &
|
|
+ (Exy * uwxy + Exz * uwxz + Eyz * uwyz) * TWO
|
|
|
|
!multiply with -1/2
|
|
Rpsi4 = - Rpsi4/TWO
|
|
Ipsi4 = - Ipsi4/TWO
|
|
|
|
return
|
|
|
|
end subroutine getnp4_ss
|
|
!-----------------------------------------------------------------------------
|
|
!
|
|
! compute the Newman-Penrose Weyl scalar Psi4
|
|
! for BSSN dynamical variables
|
|
! for single point
|
|
!-----------------------------------------------------------------------------
|
|
|
|
subroutine getnp4_point(X, Y, Z, &
|
|
chi, trK, &
|
|
dxx,gxy,gxz,dyy,gyz,dzz, &
|
|
Axx,Axy,Axz,Ayy,Ayz,Azz, &
|
|
chix,chiy,chiz, &
|
|
trKx,trKy,trKz, &
|
|
Axxx,Axxy,Axxz, &
|
|
Axyx,Axyy,Axyz, &
|
|
Axzx,Axzy,Axzz, &
|
|
Ayyx,Ayyy,Ayyz, &
|
|
Ayzx,Ayzy,Ayzz, &
|
|
Azzx,Azzy,Azzz, &
|
|
Gamxxx,Gamxxy,Gamxxz,Gamxyy,Gamxyz,Gamxzz,&
|
|
Gamyxx,Gamyxy,Gamyxz,Gamyyy,Gamyyz,Gamyzz,&
|
|
Gamzxx,Gamzxy,Gamzxz,Gamzyy,Gamzyz,Gamzzz,&
|
|
Rxx,Rxy,Rxz,Ryy,Ryz,Rzz,&
|
|
Rpsi4, Ipsi4)
|
|
|
|
implicit none
|
|
|
|
!~~~~~~> Input parameters:
|
|
|
|
real*8, intent(in ) :: X,Y,Z
|
|
real*8,intent(in ) :: dxx,gxy,gxz,dyy,gyz,dzz
|
|
real*8,intent(in ) :: Axx,Axy,Axz,Ayy,Ayz,Azz
|
|
real*8,intent(in ) :: chi,trK
|
|
real*8,intent(in ) :: chix,chiy,chiz
|
|
real*8,intent(in ) :: trKx,trKy,trKz
|
|
! covariant derivatives when out
|
|
real*8,intent(inout) :: Axxx,Axxy,Axxz
|
|
real*8,intent(inout) :: Axyx,Axyy,Axyz
|
|
real*8,intent(inout) :: Axzx,Axzy,Axzz
|
|
real*8,intent(inout) :: Ayyx,Ayyy,Ayyz
|
|
real*8,intent(inout) :: Ayzx,Ayzy,Ayzz
|
|
real*8,intent(inout) :: Azzx,Azzy,Azzz
|
|
! physical second kind of connection
|
|
real*8,intent(in) :: Gamxxx, Gamxxy, Gamxxz
|
|
real*8,intent(in) :: Gamxyy, Gamxyz, Gamxzz
|
|
real*8,intent(in) :: Gamyxx, Gamyxy, Gamyxz
|
|
real*8,intent(in) :: Gamyyy, Gamyyz, Gamyzz
|
|
real*8,intent(in) :: Gamzxx, Gamzxy, Gamzxz
|
|
real*8,intent(in) :: Gamzyy, Gamzyz, Gamzzz
|
|
! physical Ricci tensor
|
|
real*8,intent(in) :: Rxx,Rxy,Rxz,Ryy,Ryz,Rzz
|
|
real*8, intent(out):: Rpsi4,Ipsi4
|
|
|
|
!~~~~~~> Other variables:
|
|
|
|
real*8 :: f,fx,fy,fz
|
|
real*8 :: gxx,gyy,gzz,chipn1
|
|
real*8 :: vx,vy,vz,ux,uy,uz,wx,wy,wz
|
|
real*8 :: Exx,Exy,Exz,Eyy,Eyz,Ezz
|
|
real*8 :: Bxx,Bxy,Bxz,Byy,Byz,Bzz
|
|
real*8 :: gupxx,gupxy,gupxz
|
|
real*8 :: gupyy,gupyz,gupzz
|
|
real*8 :: uuwwxx,uuwwxy,uuwwxz,uuwwyy,uuwwyz,uuwwzz
|
|
real*8 :: uwxx,uwxy,uwxz,uwyy,uwyz,uwzz
|
|
|
|
real*8, parameter :: ZEO = 0.d0, ONE = 1.d0, TWO = 2.d0
|
|
real*8, parameter :: F1o3 = 1.d0/3.d0, FOUR = 4.d0
|
|
real*8, parameter :: SYM = 1.D0, ANTI= - 1.D0
|
|
real*8,parameter::TINYRR=1.d-14
|
|
|
|
gxx = dxx + ONE
|
|
gyy = dyy + ONE
|
|
gzz = dzz + ONE
|
|
chipn1= chi + ONE
|
|
|
|
! invert tilted metric
|
|
gupzz = gxx * gyy * gzz + gxy * gyz * gxz + gxz * gxy * gyz - &
|
|
gxz * gyy * gxz - gxy * gxy * gzz - gxx * gyz * gyz
|
|
gupxx = ( gyy * gzz - gyz * gyz ) / gupzz
|
|
gupxy = - ( gxy * gzz - gyz * gxz ) / gupzz
|
|
gupxz = ( gxy * gyz - gyy * gxz ) / gupzz
|
|
gupyy = ( gxx * gzz - gxz * gxz ) / gupzz
|
|
gupyz = - ( gxx * gyz - gxy * gxz ) / gupzz
|
|
gupzz = ( gxx * gyy - gxy * gxy ) / gupzz
|
|
|
|
! initialize U, V, W vetors
|
|
! v:r; u: phi; w: theta
|
|
#if (tetradtype == 0)
|
|
if(abs(X) < TINYRR .and. abs(Y) < TINYRR .and. abs(Z) < TINYRR)then
|
|
vx = TINYRR
|
|
vy = TINYRR
|
|
vz = TINYRR
|
|
else
|
|
vx = X
|
|
vy = Y
|
|
vz = Z
|
|
endif
|
|
if(abs(X) < TINYRR .and. abs(Y) < TINYRR)then
|
|
ux = - TINYRR
|
|
uy = TINYRR
|
|
uz = ZEO
|
|
wx = TINYRR*Z
|
|
wy = TINYRR*Z
|
|
wz = -2*TINYRR*TINYRR
|
|
else
|
|
ux = - Y
|
|
uy = X
|
|
uz = ZEO
|
|
wx = X*Z
|
|
wy = Y*Z
|
|
wz = -(X*X + Y*Y)
|
|
endif
|
|
|
|
f = 1.d0/chipn1
|
|
|
|
fx = gxx*vx*vx + gyy*vy*vy + gzz*vz*vz &
|
|
+(gxy*vx*vy + gxz*vx*vz + gyz*vy*vz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
vx = vx/fx
|
|
vy = vy/fx
|
|
vz = vz/fx
|
|
|
|
fx = gxx*vx*ux + gxy*vx*uy + gxz*vx*uz + &
|
|
gxy*vy*ux + gyy*vy*uy + gyz*vy*uz + &
|
|
gxz*vz*ux + gyz*vz*uy + gzz*vz*uz
|
|
fx = fx*f
|
|
ux = ux - fx*vx
|
|
uy = uy - fx*vy
|
|
uz = uz - fx*vz
|
|
fx = gxx*ux*ux + gyy*uy*uy + gzz*uz*uz &
|
|
+(gxy*ux*uy + gxz*ux*uz + gyz*uy*uz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
ux = ux/fx
|
|
uy = uy/fx
|
|
uz = uz/fx
|
|
|
|
fx = gxx*vx*wx + gxy*vx*wy + gxz*vx*wz + &
|
|
gxy*vy*wx + gyy*vy*wy + gyz*vy*wz + &
|
|
gxz*vz*wx + gyz*vz*wy + gzz*vz*wz
|
|
fx = fx*f
|
|
wx = wx - fx*vx
|
|
wy = wy - fx*vy
|
|
wz = wz - fx*vz
|
|
fx = gxx*ux*wx + gxy*ux*wy + gxz*ux*wz + &
|
|
gxy*uy*wx + gyy*uy*wy + gyz*uy*wz + &
|
|
gxz*uz*wx + gyz*uz*wy + gzz*uz*wz
|
|
fx = fx*f
|
|
wx = wx - fx*ux
|
|
wy = wy - fx*uy
|
|
wz = wz - fx*uz
|
|
fx = gxx*wx*wx + gyy*wy*wy + gzz*wz*wz &
|
|
+(gxy*wx*wy + gxz*wx*wz + gyz*wy*wz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
wx = wx/fx
|
|
wy = wy/fx
|
|
wz = wz/fx
|
|
#elif (tetradtype == 1)
|
|
if(abs(X) < TINYRR .and. abs(Y) < TINYRR .and. abs(Z) < TINYRR)then
|
|
vx = TINYRR
|
|
vy = TINYRR
|
|
vz = TINYRR
|
|
else
|
|
vx = X
|
|
vy = Y
|
|
vz = Z
|
|
endif
|
|
if(abs(X) < TINYRR .and. abs(Y) < TINYRR)then
|
|
ux = - TINYRR
|
|
uy = TINYRR
|
|
uz = ZEO
|
|
wx = TINYRR*Z
|
|
wy = TINYRR*Z
|
|
wz = -2*TINYRR*TINYRR
|
|
else
|
|
ux = - Y
|
|
uy = X
|
|
uz = ZEO
|
|
wx = X*Z
|
|
wy = Y*Z
|
|
wz = -(X*X + Y*Y)
|
|
endif
|
|
|
|
f = 1.d0/chipn1
|
|
|
|
fx = gxx*wx*wx + gyy*wy*wy + gzz*wz*wz &
|
|
+(gxy*wx*wy + gxz*wx*wz + gyz*wy*wz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
wx = wx/fx
|
|
wy = wy/fx
|
|
wz = wz/fx
|
|
|
|
fx = gxx*wx*ux + gxy*wx*uy + gxz*wx*uz + &
|
|
gxy*wy*ux + gyy*wy*uy + gyz*wy*uz + &
|
|
gxz*wz*ux + gyz*wz*uy + gzz*wz*uz
|
|
fx = fx*f
|
|
ux = ux - fx*wx
|
|
uy = uy - fx*wy
|
|
uz = uz - fx*wz
|
|
fx = gxx*ux*ux + gyy*uy*uy + gzz*uz*uz &
|
|
+(gxy*ux*uy + gxz*ux*uz + gyz*uy*uz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
ux = ux/fx
|
|
uy = uy/fx
|
|
uz = uz/fx
|
|
|
|
fx = gxx*vx*wx + gxy*vx*wy + gxz*vx*wz + &
|
|
gxy*vy*wx + gyy*vy*wy + gyz*vy*wz + &
|
|
gxz*vz*wx + gyz*vz*wy + gzz*vz*wz
|
|
fx = fx*f
|
|
vx = vx - fx*wx
|
|
vy = vy - fx*wy
|
|
vz = vz - fx*wz
|
|
fx = gxx*ux*vx + gxy*ux*vy + gxz*ux*vz + &
|
|
gxy*uy*vx + gyy*uy*vy + gyz*uy*vz + &
|
|
gxz*uz*vx + gyz*uz*vy + gzz*uz*vz
|
|
fx = fx*f
|
|
vx = vx - fx*ux
|
|
vy = vy - fx*uy
|
|
vz = vz - fx*uz
|
|
fx = gxx*vx*vx + gyy*vy*vy + gzz*vz*vz &
|
|
+(gxy*vx*vy + gxz*vx*vz + gyz*vy*vz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
vx = vx/fx
|
|
vy = vy/fx
|
|
vz = vz/fx
|
|
#elif (tetradtype == 2)
|
|
if(abs(X) < TINYRR .and. abs(Y) < TINYRR .and. abs(Z) < TINYRR)then
|
|
vx = TINYRR
|
|
vy = TINYRR
|
|
vz = TINYRR
|
|
else
|
|
vx = X
|
|
vy = Y
|
|
vz = Z
|
|
endif
|
|
if(abs(X) < TINYRR .and. abs(Y) < TINYRR)then
|
|
ux = - TINYRR
|
|
uy = TINYRR
|
|
uz = ZEO
|
|
wx = TINYRR*Z
|
|
wy = TINYRR*Z
|
|
wz = -2*TINYRR*TINYRR
|
|
else
|
|
ux = - Y
|
|
uy = X
|
|
uz = ZEO
|
|
wx = X*Z
|
|
wy = Y*Z
|
|
wz = -(X*X + Y*Y)
|
|
endif
|
|
|
|
fx = vx
|
|
fy = vy
|
|
fz = vz
|
|
vx = gupxx*fx + gupxy*fy + gupxz*fz
|
|
vy = gupxy*fx + gupyy*fy + gupyz*fz
|
|
vz = gupxz*fx + gupyz*fy + gupzz*fz
|
|
|
|
f = 1.d0/chipn1
|
|
|
|
fx = gxx*vx*vx + gyy*vy*vy + gzz*vz*vz &
|
|
+(gxy*vx*vy + gxz*vx*vz + gyz*vy*vz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
vx = vx/fx
|
|
vy = vy/fx
|
|
vz = vz/fx
|
|
|
|
fx = gxx*vx*ux + gxy*vx*uy + gxz*vx*uz + &
|
|
gxy*vy*ux + gyy*vy*uy + gyz*vy*uz + &
|
|
gxz*vz*ux + gyz*vz*uy + gzz*vz*uz
|
|
fx = fx*f
|
|
ux = ux - fx*vx
|
|
uy = uy - fx*vy
|
|
uz = uz - fx*vz
|
|
fx = gxx*ux*ux + gyy*uy*uy + gzz*uz*uz &
|
|
+(gxy*ux*uy + gxz*ux*uz + gyz*uy*uz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
ux = ux/fx
|
|
uy = uy/fx
|
|
uz = uz/fx
|
|
|
|
fx = gxx*vx*wx + gxy*vx*wy + gxz*vx*wz + &
|
|
gxy*vy*wx + gyy*vy*wy + gyz*vy*wz + &
|
|
gxz*vz*wx + gyz*vz*wy + gzz*vz*wz
|
|
fx = fx*f
|
|
wx = wx - fx*vx
|
|
wy = wy - fx*vy
|
|
wz = wz - fx*vz
|
|
fx = gxx*ux*wx + gxy*ux*wy + gxz*ux*wz + &
|
|
gxy*uy*wx + gyy*uy*wy + gyz*uy*wz + &
|
|
gxz*uz*wx + gyz*uz*wy + gzz*uz*wz
|
|
fx = fx*f
|
|
wx = wx - fx*ux
|
|
wy = wy - fx*uy
|
|
wz = wz - fx*uz
|
|
fx = gxx*wx*wx + gyy*wy*wy + gzz*wz*wz &
|
|
+(gxy*wx*wy + gxz*wx*wz + gyz*wy*wz)*TWO
|
|
fx = dsqrt(fx*f)
|
|
wx = wx/fx
|
|
wy = wy/fx
|
|
wz = wz/fx
|
|
#endif
|
|
|
|
! compute D_k K_ij up to chi^-1
|
|
Axxx = Axxx - (Gamxxx*Axx + Gamyxx*Axy + Gamzxx*Axz)*TWO - chix/chipn1*Axx + F1o3*gxx*trKx
|
|
Axxy = Axxy - (Gamxxy*Axx + Gamyxy*Axy + Gamzxy*Axz)*TWO - chiy/chipn1*Axx + F1o3*gxx*trKy
|
|
Axxz = Axxz - (Gamxxz*Axx + Gamyxz*Axy + Gamzxz*Axz)*TWO - chiz/chipn1*Axx + F1o3*gxx*trKz
|
|
Ayyx = Ayyx - (Gamxxy*Axy + Gamyxy*Ayy + Gamzxy*Ayz)*TWO - chix/chipn1*Ayy + F1o3*gyy*trKx
|
|
Ayyy = Ayyy - (Gamxyy*Axy + Gamyyy*Ayy + Gamzyy*Ayz)*TWO - chiy/chipn1*Ayy + F1o3*gyy*trKy
|
|
Ayyz = Ayyz - (Gamxyz*Axy + Gamyyz*Ayy + Gamzyz*Ayz)*TWO - chiz/chipn1*Ayy + F1o3*gyy*trKz
|
|
Azzx = Azzx - (Gamxxz*Axz + Gamyxz*Ayz + Gamzxz*Azz)*TWO - chix/chipn1*Azz + F1o3*gzz*trKx
|
|
Azzy = Azzy - (Gamxyz*Axz + Gamyyz*Ayz + Gamzyz*Azz)*TWO - chiy/chipn1*Azz + F1o3*gzz*trKy
|
|
Azzz = Azzz - (Gamxzz*Axz + Gamyzz*Ayz + Gamzzz*Azz)*TWO - chiz/chipn1*Azz + F1o3*gzz*trKz
|
|
Axyx = Axyx - (Gamxxy*Axx + Gamyxy*Axy + Gamzxy*Axz + &
|
|
Gamxxx*Axy + Gamyxx*Ayy + Gamzxx*Ayz) - chix/chipn1*Axy + F1o3*gxy*trKx
|
|
Axyy = Axyy - (Gamxyy*Axx + Gamyyy*Axy + Gamzyy*Axz + &
|
|
Gamxxy*Axy + Gamyxy*Ayy + Gamzxy*Ayz) - chiy/chipn1*Axy + F1o3*gxy*trKy
|
|
Axyz = Axyz - (Gamxyz*Axx + Gamyyz*Axy + Gamzyz*Axz + &
|
|
Gamxxz*Axy + Gamyxz*Ayy + Gamzxz*Ayz) - chiz/chipn1*Axy + F1o3*gxy*trKz
|
|
Axzx = Axzx - (Gamxxz*Axx + Gamyxz*Axy + Gamzxz*Axz + &
|
|
Gamxxx*Axz + Gamyxx*Ayz + Gamzxx*Azz) - chix/chipn1*Axz + F1o3*gxz*trKx
|
|
Axzy = Axzy - (Gamxyz*Axx + Gamyyz*Axy + Gamzyz*Axz + &
|
|
Gamxxy*Axz + Gamyxy*Ayz + Gamzxy*Azz) - chiy/chipn1*Axz + F1o3*gxz*trKy
|
|
Axzz = Axzz - (Gamxzz*Axx + Gamyzz*Axy + Gamzzz*Axz + &
|
|
Gamxxz*Axz + Gamyxz*Ayz + Gamzxz*Azz) - chiz/chipn1*Axz + F1o3*gxz*trKz
|
|
Ayzx = Ayzx - (Gamxxz*Axy + Gamyxz*Ayy + Gamzxz*Ayz + &
|
|
Gamxxy*Axz + Gamyxy*Ayz + Gamzxy*Azz) - chix/chipn1*Ayz + F1o3*gyz*trKx
|
|
Ayzy = Ayzy - (Gamxyz*Axy + Gamyyz*Ayy + Gamzyz*Ayz + &
|
|
Gamxyy*Axz + Gamyyy*Ayz + Gamzyy*Azz) - chiy/chipn1*Ayz + F1o3*gyz*trKy
|
|
Ayzz = Ayzz - (Gamxzz*Axy + Gamyzz*Ayy + Gamzzz*Ayz + &
|
|
Gamxyz*Axz + Gamyyz*Ayz + Gamzyz*Azz) - chiz/chipn1*Ayz + F1o3*gyz*trKz
|
|
! symmetrize B_ij = v^k (D_k K_ij -D_j K_ik)
|
|
Bxx =(vy*(Axxy - Axyx) + vz*(Axxz - Axzx))*f
|
|
Byy =(vx*(Ayyx - Axyy) + vz*(Ayyz - Ayzy))*f
|
|
Bzz =(vx*(Azzx - Axzz) + vy*(Azzy - Ayzz))*f
|
|
Bxy =(vx*(Axyx - (Axxy+Axyx)/TWO) + vy*(Axyy-Ayyx)/TWO + vz*(Axyz - (Axzy+Ayzx)/TWO))*f
|
|
Bxz =(vx*(Axzx - (Axxz+Axzx)/TWO) + vy*(Axzy - (Axyz+Ayzx)/TWO) + vz*(Axzz-Azzx)/TWO)*f
|
|
Byz =(vx*(Ayzx - (Axyz+Axzy)/TWO) + vy*(Ayzy - (Ayyz+Ayzy)/TWO) + vz*(Ayzz-Azzy)/TWO)*f
|
|
! E_ij = R_ij - K_ik * K^k_j + K * K_ij
|
|
|
|
! K_ij up to chi^-1
|
|
Axxx = Axx + F1o3*trK*gxx
|
|
Axyx = Axy + F1o3*trK*gxy
|
|
Axzx = Axz + F1o3*trK*gxz
|
|
Ayyx = Ayy + F1o3*trK*gyy
|
|
Ayzx = Ayz + F1o3*trK*gyz
|
|
Azzx = Azz + F1o3*trK*gzz
|
|
! gup and A_ijk cancel a chi^-1
|
|
Exx = gupxx * Axxx * Axxx + gupyy * Axyx * Axyx + gupzz * Axzx * Axzx + &
|
|
TWO * (gupxy * Axxx * Axyx + gupxz * Axxx * Axzx + gupyz * Axyx * Axzx)
|
|
Eyy = gupxx * Axyx * Axyx + gupyy * Ayyx * Ayyx + gupzz * Ayzx * Ayzx + &
|
|
TWO * (gupxy * Axyx * Ayyx + gupxz * Axyx * Ayzx + gupyz * Ayyx * Ayzx)
|
|
Ezz = gupxx * Axzx * Axzx + gupyy * Ayzx * Ayzx + gupzz * Azzx * Azzx + &
|
|
TWO * (gupxy * Axzx * Ayzx + gupxz * Axzx * Azzx + gupyz * Ayzx * Azzx)
|
|
Exy = gupxx * Axxx * Axyx + gupyy * Axyx * Ayyx + gupzz * Axzx * Ayzx + &
|
|
gupxy *(Axxx * Ayyx + Axyx * Axyx) + &
|
|
gupxz *(Axxx * Ayzx + Axzx * Axyx) + &
|
|
gupyz *(Axyx * Ayzx + Axzx * Ayyx)
|
|
Exz = gupxx * Axxx * Axzx + gupyy * Axyx * Ayzx + gupzz * Axzx * Azzx + &
|
|
gupxy *(Axxx * Ayzx + Axyx * Axzx) + &
|
|
gupxz *(Axxx * Azzx + Axzx * Axzx) + &
|
|
gupyz *(Axyx * Azzx + Axzx * Ayzx)
|
|
Eyz = gupxx * Axyx * Axzx + gupyy * Ayyx * Ayzx + gupzz * Ayzx * Azzx + &
|
|
gupxy *(Axyx * Ayzx + Ayyx * Axzx) + &
|
|
gupxz *(Axyx * Azzx + Ayzx * Axzx) + &
|
|
gupyz *(Ayyx * Azzx + Ayzx * Ayzx)
|
|
|
|
Exx = Rxx - (Exx - Axxx*trK)*f - Bxx
|
|
Exy = Rxy - (Exy - Axyx*trK)*f - Bxy
|
|
Exz = Rxz - (Exz - Axzx*trK)*f - Bxz
|
|
Eyy = Ryy - (Eyy - Ayyx*trK)*f - Byy
|
|
Eyz = Ryz - (Eyz - Ayzx*trK)*f - Byz
|
|
Ezz = Rzz - (Ezz - Azzx*trK)*f - Bzz
|
|
!set m = (u - iw)/sqrt(2) following Frans, PRD 75, 124018(2007)
|
|
! compute uuww^ij = u^i * u^j - w^i * w^j
|
|
uuwwxx = ux * ux - wx * wx
|
|
uuwwxy = ux * uy - wx * wy
|
|
uuwwxz = ux * uz - wx * wz
|
|
uuwwyy = uy * uy - wy * wy
|
|
uuwwyz = uy * uz - wy * wz
|
|
uuwwzz = uz * uz - wz * wz
|
|
|
|
! compute uw^ij = u^i * w^j + w^i * u^j
|
|
uwxx = ux * wx + wx * ux
|
|
uwxy = ux * wy + wx * uy
|
|
uwxz = ux * wz + wx * uz
|
|
uwyy = uy * wy + wy * uy
|
|
uwyz = uy * wz + wy * uz
|
|
uwzz = uz * wz + wz * uz
|
|
!the real part of Psi4
|
|
Rpsi4 = Exx * uuwwxx + Eyy * uuwwyy + Ezz * uuwwzz &
|
|
+ (Exy * uuwwxy + Exz * uuwwxz + Eyz * uuwwyz) * TWO
|
|
|
|
!the imaginary part of Psi4
|
|
Ipsi4 = Exx * uwxx + Eyy * uwyy + Ezz * uwzz &
|
|
+ (Exy * uwxy + Exz * uwxz + Eyz * uwyz) * TWO
|
|
|
|
!multiply with -1/2
|
|
Rpsi4 = - Rpsi4/TWO
|
|
Ipsi4 = - Ipsi4/TWO
|
|
|
|
return
|
|
|
|
end subroutine getnp4_point
|