/*@@ @file exa_IJ.F @date October 1999 @author Miguel Alcubierre @desc Calculate the curvature invariants I and J at the positions of the geodesics in exact spacerime. @enddesc @@*/ #include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_Functions.h" /*@@ @routine exa_IJ @date December 1997 @author Miguel Alcubierre @desc @enddesc @calls exactdata @calledby @history @endhistory @@*/ subroutine exa_IJ(CCTK_ARGUMENTS) use arrays implicit none DECLARE_CCTK_PARAMETERS DECLARE_CCTK_ARGUMENTS integer i,j,k,l,m,o,p,q CCTK_REAL dx,dy,dz CCTK_REAL xs,ys,zs,ts CCTK_REAL xp,xm,yp,ym,zp,zm CCTK_REAL DET,KT CCTK_REAL zero,sixth,half,one,two,three CCTK_REAL teps,ieps CCTK_REAL aux CCTK_REAL gxx_junk,gyy_junk,gzz_junk CCTK_REAL gxy_junk,gxz_junk,gyz_junk CCTK_REAL alp_junk,ax_junk,ay_junk,az_junk CCTK_REAL betax_junk,betay_junk,betaz_junk CCTK_REAL bxx_junk,bxy_junk,bxz_junk CCTK_REAL byx_junk,byy_junk,byz_junk CCTK_REAL bzx_junk,bzy_junk,bzz_junk CCTK_REAL, dimension (3,3) :: GL,GU CCTK_REAL, dimension (3,3) :: KL,KP,KM,K2 CCTK_REAL, dimension (3,3) :: RL,EL,BL CCTK_REAL, dimension (3,3,3) :: e,DK,CDK CCTK_REAL, dimension (3,3,3) :: DG,DP,DM,GAMMA CCTK_REAL, dimension (3,3,3,3) :: DD,DGAM ! Description of variables: ! ! i,j,k,l Counters. ! m,o,p,q Counters. ! ! ts Coordinate time. ! ! xs,ys,zs Position. ! ! ! e(i,j,k) Levi-Civita totally antisymmetric tensor. ! ! ! GL(i,j) Spatial metric. ! ! ! GU(i,j) Inverse spatial metric. ! ! ! KL(i,j) Extrinsic curvature. ! ! ! DET det(g) ! ! ! DG(i,j,k) d g / 2 ! i kj ! _ i ! GAMMA(i,j,k) Chistoffel symbols | ! jk ! ! KT Trace of the extrinsic curvature. ! ! k ! K2(i,j) K K ! ik j ! ! DK(i,j,k) d K ! i jk ! _ m _m ! CDK(i,j,k) d K - | K - | K ! i jk ij mk ik mj ! ! DD(i,j,k,l) d DG ! i jkl ! _ j ! DGAM(i,j,k,l) d | ! i kl ! ! RL(i,j) Ricci tensor. ! ! ! EL(i,j) Electric part of the Riemann tensor. ! ! ! BL(i,j) Magnetic part of the Riemann tensor. ! ! ! *P Variables calculated at coordinate + teps. ! ! ! *M Variables calculated at coordinate - teps. ! ! ! *_junk Variables imported from exact solution ! that are not needed later for anything. ! ! teps Small number for calculating derivatives ! using finite differences. ! ! ieps 1 / (2 teps) ! ************************** ! *** DEFINE NUMBERS *** ! ************************** zero = 0.0d0 sixth = 1.0D0/6.0D0 half = 0.5d0 one = 1.0d0 two = 2.0d0 three = 3.0d0 ! Grid parameters. dx = cctk_delta_space(1) dy = cctk_delta_space(2) dz = cctk_delta_space(3) ! **************************************** ! *** CHOOSE teps AND FIND ieps *** ! **************************************** ! Do not be tempted to use a very small "teps" here just ! because you have an exact solution. In fact, if the ! exact solution involves square roots (and many do) ! then it will only be accurate to single precision. ! Moreover, if it does not change very fast then a small ! value of "teps" (which might be as large as 0.0001) will ! rapidly get you to the level of the roundoff error. ! SO BE WARNED! teps = dx ieps = half/teps ! ************************************* ! *** LOOP OVER THE GRID POINTS *** ! ************************************* do i=1,lgeos ! ************************************************** ! *** FIND EXACT VALUES OF GEOMETRIC OBJECTS *** ! ************************************************** ts = geote(i) xs = geoxe(i) ys = geoye(i) zs = geoze(i) #ifdef THORN_EXACT call exactdata(xs,ys,zs,ts, & GL(1,1),GL(2,2),GL(3,3),GL(1,2),GL(2,3),GL(1,3), & KL(1,1),KL(2,2),KL(3,3),KL(1,2),KL(2,3),KL(1,3), & DG(1,1,1),DG(1,2,2),DG(1,3,3),DG(1,1,2),DG(1,2,3),DG(1,1,3), & DG(2,1,1),DG(2,2,2),DG(2,3,3),DG(2,1,2),DG(2,2,3),DG(2,1,3), & DG(3,1,1),DG(3,2,2),DG(3,3,3),DG(3,1,2),DG(3,2,3),DG(3,1,3), & alp_junk,ax_junk,ay_junk,az_junk, & betax_junk,betay_junk,betaz_junk, & bxx_junk,bxy_junk,bxz_junk, & byx_junk,byy_junk,byz_junk, & bzx_junk,bzy_junk,bzz_junk) #else GL = 0 KL = 0 DG = 0 KM = 0 KP = 0 DM = 0 DP = 0 #endif GL(2,1) = GL(1,2) GL(3,1) = GL(1,3) GL(3,2) = GL(2,3) KL(2,1) = KL(1,2) KL(3,1) = KL(1,3) KL(3,2) = KL(2,3) do l=1,3 DG(l,2,1) = DG(l,1,2) DG(l,3,1) = DG(l,1,3) DG(l,3,2) = DG(l,2,3) end do ! ******************************* ! *** FIND INVERSE METRIC *** ! ******************************* ! Find the determinant of the 3-metric. DET = GL(1,1)*GL(2,2)*GL(3,3) & + two*GL(1,2)*GL(1,3)*GL(2,3) & - (GL(1,1)*GL(2,3)**2 + GL(2,2)*GL(1,3)**2 & + GL(3,3)*GL(1,2)**2) aux = one/DET ! Find inverse metric components. GU(1,1) = aux*(GL(2,2)*GL(3,3) - GL(2,3)**2) GU(2,2) = aux*(GL(1,1)*GL(3,3) - GL(1,3)**2) GU(3,3) = aux*(GL(1,1)*GL(2,2) - GL(1,2)**2) GU(1,2) = aux*(GL(1,3)*GL(2,3) - GL(3,3)*GL(1,2)) GU(1,3) = aux*(GL(1,2)*GL(2,3) - GL(2,2)*GL(1,3)) GU(2,3) = aux*(GL(1,2)*GL(1,3) - GL(1,1)*GL(2,3)) GU(2,1) = GU(1,2) GU(3,1) = GU(1,3) GU(3,2) = GU(2,3) ! ************************** ! *** FIND e(i,j,k) *** ! ************************** ! The Levi-Civita tensor is defined as: ! ! / +1 if {i,j,k} is an even permutation of {1,2,3} ! e = sqrt(g) | ! ijk \ -1 if {i,j,k} is an odd permutation of {1,2,3} e = zero aux = dsqrt(DET) e(1,2,3) = + aux e(1,3,2) = - aux e(2,1,3) = - aux e(2,3,1) = + aux e(3,1,2) = + aux e(3,2,1) = - aux ! ************************* ! *** FIND {KT,K2} *** ! ************************* KT = zero do m=1,3 do l=1,3 KT = KT + GU(l,m)*KL(l,m) end do end do do m=1,3 do l=1,3 aux = zero do p=1,3 do q=1,3 aux = aux + GU(p,q)*KL(l,p)*KL(m,q) end do end do K2(l,m) = aux end do end do ! ************************************ ! *** FIND CHRISTOFFEL SYMBOLS *** ! ************************************ ! Remember that the D are 1/2 of the metric derivatives! do l=1,3 do k=1,3 do j=1,3 aux = zero do m=1,3 aux = aux + GU(j,m) & *(DG(k,l,m) + DG(l,k,m) - DG(m,k,l)) end do GAMMA(j,k,l) = aux end do end do end do ! ***************************************** ! *** FIND DERIVATIVES OF K AND D *** ! ***************************************** xp = xs + teps xm = xs - teps yp = ys + teps ym = ys - teps zp = zs + teps zm = zs - teps ! Find x derivatives. #ifdef THORN_EXACT call exactdata(xp,ys,zs,ts, & gxx_junk,gyy_junk,gzz_junk,gxy_junk,gyz_junk,gxz_junk, & KP(1,1),KP(2,2),KP(3,3),KP(1,2),KP(2,3),KP(1,3), & DP(1,1,1),DP(1,2,2),DP(1,3,3),DP(1,1,2),DP(1,2,3),DP(1,1,3), & DP(2,1,1),DP(2,2,2),DP(2,3,3),DP(2,1,2),DP(2,2,3),DP(2,1,3), & DP(3,1,1),DP(3,2,2),DP(3,3,3),DP(3,1,2),DP(3,2,3),DP(3,1,3), & alp_junk,ax_junk,ay_junk,az_junk, & betax_junk,betay_junk,betaz_junk, & bxx_junk,bxy_junk,bxz_junk, & byx_junk,byy_junk,byz_junk, & bzx_junk,bzy_junk,bzz_junk) KP(2,1) = KP(1,2) KP(3,1) = KP(1,3) KP(3,2) = KP(2,3) do l=1,3 DP(l,2,1) = DP(l,1,2) DP(l,3,1) = DP(l,1,3) DP(l,3,2) = DP(l,2,3) end do call exactdata(xm,ys,zs,ts, & gxx_junk,gyy_junk,gzz_junk,gxy_junk,gyz_junk,gxz_junk, & KM(1,1),KM(2,2),KM(3,3),KM(1,2),KM(2,3),KM(1,3), & DM(1,1,1),DM(1,2,2),DM(1,3,3),DM(1,1,2),DM(1,2,3),DM(1,1,3), & DM(2,1,1),DM(2,2,2),DM(2,3,3),DM(2,1,2),DM(2,2,3),DM(2,1,3), & DM(3,1,1),DM(3,2,2),DM(3,3,3),DM(3,1,2),DM(3,2,3),DM(3,1,3), & alp_junk,ax_junk,ay_junk,az_junk, & betax_junk,betay_junk,betaz_junk, & bxx_junk,bxy_junk,bxz_junk, & byx_junk,byy_junk,byz_junk, & bzx_junk,bzy_junk,bzz_junk) KM(2,1) = KM(1,2) KM(3,1) = KM(1,3) KM(3,2) = KM(2,3) do l=1,3 DM(l,2,1) = DM(l,1,2) DM(l,3,1) = DM(l,1,3) DM(l,3,2) = DM(l,2,3) end do #endif do m=1,3 do l=1,3 DK(1,l,m) = ieps*(KP(l,m) - KM(l,m)) end do end do do p=1,3 do m=1,3 do l=1,3 DD(1,l,m,p) = ieps*(DP(l,m,p) - DM(l,m,p)) end do end do end do ! Find y derivatives. #ifdef THORN_EXACT call exactdata(xs,yp,zs,ts, & gxx_junk,gyy_junk,gzz_junk,gxy_junk,gyz_junk,gxz_junk, & KP(1,1),KP(2,2),KP(3,3),KP(1,2),KP(2,3),KP(1,3), & DP(1,1,1),DP(1,2,2),DP(1,3,3),DP(1,1,2),DP(1,2,3),DP(1,1,3), & DP(2,1,1),DP(2,2,2),DP(2,3,3),DP(2,1,2),DP(2,2,3),DP(2,1,3), & DP(3,1,1),DP(3,2,2),DP(3,3,3),DP(3,1,2),DP(3,2,3),DP(3,1,3), & alp_junk,ax_junk,ay_junk,az_junk, & betax_junk,betay_junk,betaz_junk, & bxx_junk,bxy_junk,bxz_junk, & byx_junk,byy_junk,byz_junk, & bzx_junk,bzy_junk,bzz_junk) KP(2,1) = KP(1,2) KP(3,1) = KP(1,3) KP(3,2) = KP(2,3) do l=1,3 DP(l,2,1) = DP(l,1,2) DP(l,3,1) = DP(l,1,3) DP(l,3,2) = DP(l,2,3) end do call exactdata(xs,ym,zs,ts, & gxx_junk,gyy_junk,gzz_junk,gxy_junk,gyz_junk,gxz_junk, & KM(1,1),KM(2,2),KM(3,3),KM(1,2),KM(2,3),KM(1,3), & DM(1,1,1),DM(1,2,2),DM(1,3,3),DM(1,1,2),DM(1,2,3),DM(1,1,3), & DM(2,1,1),DM(2,2,2),DM(2,3,3),DM(2,1,2),DM(2,2,3),DM(2,1,3), & DM(3,1,1),DM(3,2,2),DM(3,3,3),DM(3,1,2),DM(3,2,3),DM(3,1,3), & alp_junk,ax_junk,ay_junk,az_junk, & betax_junk,betay_junk,betaz_junk, & bxx_junk,bxy_junk,bxz_junk, & byx_junk,byy_junk,byz_junk, & bzx_junk,bzy_junk,bzz_junk) KM(2,1) = KM(1,2) KM(3,1) = KM(1,3) KM(3,2) = KM(2,3) do l=1,3 DM(l,2,1) = DM(l,1,2) DM(l,3,1) = DM(l,1,3) DM(l,3,2) = DM(l,2,3) end do #endif do m=1,3 do l=1,3 DK(2,l,m) = ieps*(KP(l,m) - KM(l,m)) end do end do do p=1,3 do m=1,3 do l=1,3 DD(2,l,m,p) = ieps*(DP(l,m,p) - DM(l,m,p)) end do end do end do ! Find z derivatives. #ifdef THORN_EXACT call exactdata(xs,ys,zp,ts, & gxx_junk,gyy_junk,gzz_junk,gxy_junk,gyz_junk,gxz_junk, & KP(1,1),KP(2,2),KP(3,3),KP(1,2),KP(2,3),KP(1,3), & DP(1,1,1),DP(1,2,2),DP(1,3,3),DP(1,1,2),DP(1,2,3),DP(1,1,3), & DP(2,1,1),DP(2,2,2),DP(2,3,3),DP(2,1,2),DP(2,2,3),DP(2,1,3), & DP(3,1,1),DP(3,2,2),DP(3,3,3),DP(3,1,2),DP(3,2,3),DP(3,1,3), & alp_junk,ax_junk,ay_junk,az_junk, & betax_junk,betay_junk,betaz_junk, & bxx_junk,bxy_junk,bxz_junk, & byx_junk,byy_junk,byz_junk, & bzx_junk,bzy_junk,bzz_junk) KP(2,1) = KP(1,2) KP(3,1) = KP(1,3) KP(3,2) = KP(2,3) do l=1,3 DP(l,2,1) = DP(l,1,2) DP(l,3,1) = DP(l,1,3) DP(l,3,2) = DP(l,2,3) end do call exactdata(xs,ys,zm,ts, & gxx_junk,gyy_junk,gzz_junk,gxy_junk,gyz_junk,gxz_junk, & KM(1,1),KM(2,2),KM(3,3),KM(1,2),KM(2,3),KM(1,3), & DM(1,1,1),DM(1,2,2),DM(1,3,3),DM(1,1,2),DM(1,2,3),DM(1,1,3), & DM(2,1,1),DM(2,2,2),DM(2,3,3),DM(2,1,2),DM(2,2,3),DM(2,1,3), & DM(3,1,1),DM(3,2,2),DM(3,3,3),DM(3,1,2),DM(3,2,3),DM(3,1,3), & alp_junk,ax_junk,ay_junk,az_junk, & betax_junk,betay_junk,betaz_junk, & bxx_junk,bxy_junk,bxz_junk, & byx_junk,byy_junk,byz_junk, & bzx_junk,bzy_junk,bzz_junk) KM(2,1) = KM(1,2) KM(3,1) = KM(1,3) KM(3,2) = KM(2,3) do l=1,3 DM(l,2,1) = DM(l,1,2) DM(l,3,1) = DM(l,1,3) DM(l,3,2) = DM(l,2,3) end do #endif do m=1,3 do l=1,3 DK(3,l,m) = ieps*(KP(l,m) - KM(l,m)) end do end do do p=1,3 do m=1,3 do l=1,3 DD(3,l,m,p) = ieps*(DP(l,m,p) - DM(l,m,p)) end do end do end do ! ******************************************** ! *** FIND COVARIANT DERIVATIVES OF K *** ! ******************************************** do m=1,3 do k=1,3 do j=1,3 aux = zero do p=1,3 aux = aux & - GAMMA(p,k,j)*KL(p,m) & - GAMMA(p,m,j)*KL(k,p) end do CDK(j,k,m) = DK(j,k,m) + aux end do end do end do ! *************************************************** ! *** FIND DERIVATIVES OF CHRISTOFFEL SYMBOLS *** ! *************************************************** ! Remember that the Dijk are 1/2 of the metric derivatives! do m=1,3 do l=1,3 do k=1,3 do j=1,3 aux = zero do p=1,3 aux = aux + GU(k,p) & *(DD(j,l,m,p) + DD(j,m,l,p) - DD(j,p,l,m)) end do do o=1,3 do p=1,3 do q=1,3 aux = aux & - two*GU(k,p)*GU(o,q)*DG(j,p,q) & *(DG(l,m,o) + DG(m,l,o) - DG(o,l,m)) end do end do end do DGAM(j,k,l,m) = aux end do end do end do end do ! ********************************* ! *** FIND THE RICCI TENSOR *** ! ********************************* do m=1,3 do l=1,3 aux = zero do p=1,3 aux = aux + DGAM(p,p,l,m) - DGAM(m,p,p,l) end do do p=1,3 do q=1,3 aux = aux & + GAMMA(p,p,q)*GAMMA(q,l,m) & - GAMMA(p,q,l)*GAMMA(q,p,m) end do end do RL(l,m) = aux end do end do ! ************************************************** ! *** FIND THE ELECTRIC AND MAGNETIC TENSORS *** ! ************************************************** ! The electric part of the Riemman tensor is defined as: ! ! k ! E = - R - trK K + K K ! ij ij ij ik j do m=1,3 do l=1,3 EL(l,m) = - RL(l,m) - KL(l,m)*KT + K2(l,m) end do end do ! The magnetic part of the Riemman tensor is defined as: ! ! pq ! B = e K ! ij i pq;j do m=1,3 do l=1,3 aux = zero do j=1,3 do k=1,3 do p=1,3 do q=1,3 aux = aux + GU(j,p)*GU(k,q) & *e(l,p,q)*CDK(j,k,m) end do end do end do end do BL(l,m) = aux end do end do ! *********************************** ! *** FIND INVARIANTS I and J *** ! *********************************** ! Notice that the normalization of the invariants ! I and J is completly arbitrary. The particular ! factors used here might seem a bit odd, but they ! were chosen to be compatible with those used in ! thorn_PsiKadelia. ! Real part of I: ! ! ij ij ! Re(I) = 1/2 ( E E - B B ) ! ij ij aux = zero do l=1,3 do m=1,3 do p=1,3 do q=1,3 aux = aux + GU(l,p)*GU(m,q) & *(EL(l,m)*EL(p,q) - BL(l,m)*BL(p,q)) end do end do end do end do geoexaI_r(i) = half*aux ! Imaginary part of I: ! ! ij ! Im(I) = E B ! ij aux = zero do l=1,3 do m=1,3 do p=1,3 do q=1,3 aux = aux + GU(l,p)*GU(m,q)*EL(l,m)*BL(p,q) end do end do end do end do geoexaI_i(i) = aux ! Real part of J: ! ! i jk i jk ! Re(J) = + 1/6 E ( E E - 3 B B ) ! ij k k aux = zero do j=1,3 do k=1,3 do l=1,3 do m=1,3 do p=1,3 do q=1,3 aux = aux + GU(j,m)*GU(k,p)*GU(l,q) & * EL(j,k)*(EL(m,l)*EL(p,q) & - three*BL(m,l)*BL(p,q)) end do end do end do end do end do end do geoexaJ_r(i) = sixth*aux ! Imaginary part of J: ! ! jk i i ! Im(J) = - 1/6 B ( B B - 3 E E ) ! ij k ij k aux = zero do j=1,3 do k=1,3 do l=1,3 do m=1,3 do p=1,3 do q=1,3 aux = aux - GU(j,m)*GU(k,p)*GU(l,q) & * BL(p,q)*(BL(j,k)*BL(m,l) & - three*EL(j,k)*EL(m,l)) end do end do end do end do end do end do geoexaJ_i(i) = sixth*aux end do ! *************** ! *** END *** ! *************** return end