#include "macrodef.fh" function compute_rhs_bssn(ex, T,X, Y, Z, & chi , trK , & dxx , gxy , gxz , dyy , gyz , dzz, & Axx , Axy , Axz , Ayy , Ayz , Azz, & Gamx , Gamy , Gamz , & Lap , betax , betay , betaz , & dtSfx , dtSfy , dtSfz , & chi_rhs, trK_rhs, & gxx_rhs, gxy_rhs, gxz_rhs, gyy_rhs, gyz_rhs, gzz_rhs, & Axx_rhs, Axy_rhs, Axz_rhs, Ayy_rhs, Ayz_rhs, Azz_rhs, & Gamx_rhs, Gamy_rhs, Gamz_rhs, & Lap_rhs, betax_rhs, betay_rhs, betaz_rhs, & dtSfx_rhs, dtSfy_rhs, dtSfz_rhs, & rho,Sx,Sy,Sz,Sxx,Sxy,Sxz,Syy,Syz,Szz, & Gamxxx,Gamxxy,Gamxxz,Gamxyy,Gamxyz,Gamxzz, & Gamyxx,Gamyxy,Gamyxz,Gamyyy,Gamyyz,Gamyzz, & Gamzxx,Gamzxy,Gamzxz,Gamzyy,Gamzyz,Gamzzz, & Rxx,Rxy,Rxz,Ryy,Ryz,Rzz, & ham_Res, movx_Res, movy_Res, movz_Res, & Gmx_Res, Gmy_Res, Gmz_Res, & Symmetry,Lev,eps,co) result(gont) ! calculate constraint violation when co=0 implicit none !~~~~~~> Input parameters: integer,intent(in ):: ex(1:3), Symmetry,Lev,co real*8, intent(in ):: T 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(inout) :: chi,dxx,dyy,dzz real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: trK real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: gxy,gxz,gyz 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 ) :: Gamx,Gamy,Gamz real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Lap, betax, betay, betaz real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: dtSfx, dtSfy, dtSfz real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: chi_rhs,trK_rhs real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: gxx_rhs,gxy_rhs,gxz_rhs real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: gyy_rhs,gyz_rhs,gzz_rhs real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Axx_rhs,Axy_rhs,Axz_rhs real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Ayy_rhs,Ayz_rhs,Azz_rhs real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Gamx_rhs,Gamy_rhs,Gamz_rhs real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Lap_rhs, betax_rhs, betay_rhs, betaz_rhs real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: dtSfx_rhs,dtSfy_rhs,dtSfz_rhs real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: rho,Sx,Sy,Sz real*8, dimension(ex(1),ex(2),ex(3)),intent(in ) :: Sxx,Sxy,Sxz,Syy,Syz,Szz ! when out, physical second kind of connection real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Gamxxx, Gamxxy, Gamxxz real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Gamxyy, Gamxyz, Gamxzz real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Gamyxx, Gamyxy, Gamyxz real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Gamyyy, Gamyyz, Gamyzz real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Gamzxx, Gamzxy, Gamzxz real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Gamzyy, Gamzyz, Gamzzz ! when out, physical Ricci tensor real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Rxx,Rxy,Rxz,Ryy,Ryz,Rzz real*8,intent(in) :: eps real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: ham_Res, movx_Res, movy_Res, movz_Res real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Gmx_Res, Gmy_Res, Gmz_Res ! gont = 0: success; gont = 1: something wrong integer::gont integer :: i,j,k !~~~~~~> Other variables: real*8, dimension(ex(1),ex(2),ex(3)) :: gxx,gyy,gzz real*8, dimension(ex(1),ex(2),ex(3)) :: chix,chiy,chiz real*8, dimension(ex(1),ex(2),ex(3)) :: gxxx,gxyx,gxzx,gyyx,gyzx,gzzx real*8, dimension(ex(1),ex(2),ex(3)) :: gxxy,gxyy,gxzy,gyyy,gyzy,gzzy real*8, dimension(ex(1),ex(2),ex(3)) :: gxxz,gxyz,gxzz,gyyz,gyzz,gzzz real*8, dimension(ex(1),ex(2),ex(3)) :: Lapx,Lapy,Lapz real*8, dimension(ex(1),ex(2),ex(3)) :: betaxx,betaxy,betaxz real*8, dimension(ex(1),ex(2),ex(3)) :: betayx,betayy,betayz real*8, dimension(ex(1),ex(2),ex(3)) :: betazx,betazy,betazz real*8, dimension(ex(1),ex(2),ex(3)) :: Gamxx,Gamxy,Gamxz real*8, dimension(ex(1),ex(2),ex(3)) :: Gamyx,Gamyy,Gamyz real*8, dimension(ex(1),ex(2),ex(3)) :: Gamzx,Gamzy,Gamzz real*8, dimension(ex(1),ex(2),ex(3)) :: Kx,Ky,Kz,div_beta,S real*8, dimension(ex(1),ex(2),ex(3)) :: f,fxx,fxy,fxz,fyy,fyz,fzz real*8, dimension(ex(1),ex(2),ex(3)) :: Gamxa,Gamya,Gamza,alpn1,chin1 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(3) ::SSS,AAS,ASA,SAA,ASS,SAS,SSA real*8 :: dX, dY, dZ, PI real*8 :: divb_loc,det_loc real*8 :: gupxx_loc,gupxy_loc,gupxz_loc,gupyy_loc,gupyz_loc,gupzz_loc real*8 :: Rxx_loc,Rxy_loc,Rxz_loc,Ryy_loc,Ryz_loc,Rzz_loc real*8 :: fxx_loc,fxy_loc,fxz_loc real*8 :: Gamxa_loc,Gamya_loc,Gamza_loc real*8, parameter :: ZEO = 0.d0,ONE = 1.D0, TWO = 2.D0, FOUR = 4.D0 real*8, parameter :: EIGHT = 8.D0, HALF = 0.5D0, THR = 3.d0 real*8, parameter :: SYM = 1.D0, ANTI= - 1.D0 double precision,parameter::FF = 0.75d0,eta=2.d0 real*8, parameter :: F1o3 = 1.D0/3.D0, F2o3 = 2.D0/3.D0,F3o2=1.5d0, F1o6 = 1.D0/6.D0 real*8, parameter :: F16=1.6d1,F8=8.d0 #if (GAUGE == 2 || GAUGE == 3 || GAUGE == 4 || GAUGE == 5) real*8, dimension(ex(1),ex(2),ex(3)) :: reta #endif #if (GAUGE == 6 || GAUGE == 7) integer :: BHN real*8, dimension(9) :: Porg real*8, dimension(3) :: Mass real*8 :: r1,r2,M,A,w1,w2,C1,C2 real*8, dimension(ex(1),ex(2),ex(3)) :: reta call getpbh(BHN,Porg,Mass) #endif !!! sanity check (disabled in production builds for performance) #ifdef DEBUG dX = sum(chi)+sum(trK)+sum(dxx)+sum(gxy)+sum(gxz)+sum(dyy)+sum(gyz)+sum(dzz) & +sum(Axx)+sum(Axy)+sum(Axz)+sum(Ayy)+sum(Ayz)+sum(Azz) & +sum(Gamx)+sum(Gamy)+sum(Gamz) & +sum(Lap)+sum(betax)+sum(betay)+sum(betaz) if(dX.ne.dX) then if(sum(chi).ne.sum(chi))write(*,*)"bssn.f90: find NaN in chi" if(sum(trK).ne.sum(trK))write(*,*)"bssn.f90: find NaN in trk" if(sum(dxx).ne.sum(dxx))write(*,*)"bssn.f90: find NaN in dxx" if(sum(gxy).ne.sum(gxy))write(*,*)"bssn.f90: find NaN in gxy" if(sum(gxz).ne.sum(gxz))write(*,*)"bssn.f90: find NaN in gxz" if(sum(dyy).ne.sum(dyy))write(*,*)"bssn.f90: find NaN in dyy" if(sum(gyz).ne.sum(gyz))write(*,*)"bssn.f90: find NaN in gyz" if(sum(dzz).ne.sum(dzz))write(*,*)"bssn.f90: find NaN in dzz" if(sum(Axx).ne.sum(Axx))write(*,*)"bssn.f90: find NaN in Axx" if(sum(Axy).ne.sum(Axy))write(*,*)"bssn.f90: find NaN in Axy" if(sum(Axz).ne.sum(Axz))write(*,*)"bssn.f90: find NaN in Axz" if(sum(Ayy).ne.sum(Ayy))write(*,*)"bssn.f90: find NaN in Ayy" if(sum(Ayz).ne.sum(Ayz))write(*,*)"bssn.f90: find NaN in Ayz" if(sum(Azz).ne.sum(Azz))write(*,*)"bssn.f90: find NaN in Azz" if(sum(Gamx).ne.sum(Gamx))write(*,*)"bssn.f90: find NaN in Gamx" if(sum(Gamy).ne.sum(Gamy))write(*,*)"bssn.f90: find NaN in Gamy" if(sum(Gamz).ne.sum(Gamz))write(*,*)"bssn.f90: find NaN in Gamz" if(sum(Lap).ne.sum(Lap))write(*,*)"bssn.f90: find NaN in Lap" if(sum(betax).ne.sum(betax))write(*,*)"bssn.f90: find NaN in betax" if(sum(betay).ne.sum(betay))write(*,*)"bssn.f90: find NaN in betay" if(sum(betaz).ne.sum(betaz))write(*,*)"bssn.f90: find NaN in betaz" gont = 1 return endif #endif PI = dacos(-ONE) dX = X(2) - X(1) dY = Y(2) - Y(1) dZ = Z(2) - Z(1) do k=1,ex(3) do j=1,ex(2) do i=1,ex(1) alpn1(i,j,k) = Lap(i,j,k) + ONE chin1(i,j,k) = chi(i,j,k) + ONE gxx(i,j,k) = dxx(i,j,k) + ONE gyy(i,j,k) = dyy(i,j,k) + ONE gzz(i,j,k) = dzz(i,j,k) + ONE enddo enddo enddo call fderivs(ex,betax,betaxx,betaxy,betaxz,X,Y,Z,ANTI, SYM, SYM,Symmetry,Lev) call fderivs(ex,betay,betayx,betayy,betayz,X,Y,Z, SYM,ANTI, SYM,Symmetry,Lev) call fderivs(ex,betaz,betazx,betazy,betazz,X,Y,Z, SYM, SYM,ANTI,Symmetry,Lev) call fderivs(ex,chi,chix,chiy,chiz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev) call fderivs(ex,dxx,gxxx,gxxy,gxxz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev) call fderivs(ex,gxy,gxyx,gxyy,gxyz,X,Y,Z,ANTI,ANTI,SYM ,Symmetry,Lev) call fderivs(ex,gxz,gxzx,gxzy,gxzz,X,Y,Z,ANTI,SYM ,ANTI,Symmetry,Lev) call fderivs(ex,dyy,gyyx,gyyy,gyyz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev) call fderivs(ex,gyz,gyzx,gyzy,gyzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,Lev) call fderivs(ex,dzz,gzzx,gzzy,gzzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev) do k=1,ex(3) do j=1,ex(2) do i=1,ex(1) divb_loc = betaxx(i,j,k) + betayy(i,j,k) + betazz(i,j,k) div_beta(i,j,k) = divb_loc chi_rhs(i,j,k) = F2o3 * chin1(i,j,k) * (alpn1(i,j,k) * trK(i,j,k) - divb_loc) gxx_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Axx(i,j,k) - F2o3 * gxx(i,j,k) * divb_loc + & TWO * ( gxx(i,j,k) * betaxx(i,j,k) + gxy(i,j,k) * betayx(i,j,k) + gxz(i,j,k) * betazx(i,j,k) ) gyy_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Ayy(i,j,k) - F2o3 * gyy(i,j,k) * divb_loc + & TWO * ( gxy(i,j,k) * betaxy(i,j,k) + gyy(i,j,k) * betayy(i,j,k) + gyz(i,j,k) * betazy(i,j,k) ) gzz_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Azz(i,j,k) - F2o3 * gzz(i,j,k) * divb_loc + & TWO * ( gxz(i,j,k) * betaxz(i,j,k) + gyz(i,j,k) * betayz(i,j,k) + gzz(i,j,k) * betazz(i,j,k) ) gxy_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Axy(i,j,k) + F1o3 * gxy(i,j,k) * divb_loc + & gxx(i,j,k) * betaxy(i,j,k) + gxz(i,j,k) * betazy(i,j,k) + gyy(i,j,k) * betayx(i,j,k) + & gyz(i,j,k) * betazx(i,j,k) - gxy(i,j,k) * betazz(i,j,k) gyz_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Ayz(i,j,k) + F1o3 * gyz(i,j,k) * divb_loc + & gxy(i,j,k) * betaxz(i,j,k) + gyy(i,j,k) * betayz(i,j,k) + gxz(i,j,k) * betaxy(i,j,k) + & gzz(i,j,k) * betazy(i,j,k) - gyz(i,j,k) * betaxx(i,j,k) gxz_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Axz(i,j,k) + F1o3 * gxz(i,j,k) * divb_loc + & gxx(i,j,k) * betaxz(i,j,k) + gxy(i,j,k) * betayz(i,j,k) + gyz(i,j,k) * betayx(i,j,k) + & gzz(i,j,k) * betazx(i,j,k) - gxz(i,j,k) * betayy(i,j,k) det_loc = gxx(i,j,k) * gyy(i,j,k) * gzz(i,j,k) + gxy(i,j,k) * gyz(i,j,k) * gxz(i,j,k) + & gxz(i,j,k) * gxy(i,j,k) * gyz(i,j,k) - gxz(i,j,k) * gyy(i,j,k) * gxz(i,j,k) - & gxy(i,j,k) * gxy(i,j,k) * gzz(i,j,k) - gxx(i,j,k) * gyz(i,j,k) * gyz(i,j,k) gupxx_loc = ( gyy(i,j,k) * gzz(i,j,k) - gyz(i,j,k) * gyz(i,j,k) ) / det_loc gupxy_loc = - ( gxy(i,j,k) * gzz(i,j,k) - gyz(i,j,k) * gxz(i,j,k) ) / det_loc gupxz_loc = ( gxy(i,j,k) * gyz(i,j,k) - gyy(i,j,k) * gxz(i,j,k) ) / det_loc gupyy_loc = ( gxx(i,j,k) * gzz(i,j,k) - gxz(i,j,k) * gxz(i,j,k) ) / det_loc gupyz_loc = - ( gxx(i,j,k) * gyz(i,j,k) - gxy(i,j,k) * gxz(i,j,k) ) / det_loc gupzz_loc = ( gxx(i,j,k) * gyy(i,j,k) - gxy(i,j,k) * gxy(i,j,k) ) / det_loc gupxx(i,j,k) = gupxx_loc gupxy(i,j,k) = gupxy_loc gupxz(i,j,k) = gupxz_loc gupyy(i,j,k) = gupyy_loc gupyz(i,j,k) = gupyz_loc gupzz(i,j,k) = gupzz_loc if(co == 0)then Gmx_Res(i,j,k) = Gamx(i,j,k) - ( & gupxx_loc*(gupxx_loc*gxxx(i,j,k)+gupxy_loc*gxyx(i,j,k)+gupxz_loc*gxzx(i,j,k)) + & gupxy_loc*(gupxx_loc*gxyx(i,j,k)+gupxy_loc*gyyx(i,j,k)+gupxz_loc*gyzx(i,j,k)) + & gupxz_loc*(gupxx_loc*gxzx(i,j,k)+gupxy_loc*gyzx(i,j,k)+gupxz_loc*gzzx(i,j,k)) + & gupxx_loc*(gupxy_loc*gxxy(i,j,k)+gupyy_loc*gxyy(i,j,k)+gupyz_loc*gxzy(i,j,k)) + & gupxy_loc*(gupxy_loc*gxyy(i,j,k)+gupyy_loc*gyyy(i,j,k)+gupyz_loc*gyzy(i,j,k)) + & gupxz_loc*(gupxy_loc*gxzy(i,j,k)+gupyy_loc*gyzy(i,j,k)+gupyz_loc*gzzy(i,j,k)) + & gupxx_loc*(gupxz_loc*gxxz(i,j,k)+gupyz_loc*gxyz(i,j,k)+gupzz_loc*gxzz(i,j,k)) + & gupxy_loc*(gupxz_loc*gxyz(i,j,k)+gupyz_loc*gyyz(i,j,k)+gupzz_loc*gyzz(i,j,k)) + & gupxz_loc*(gupxz_loc*gxzz(i,j,k)+gupyz_loc*gyzz(i,j,k)+gupzz_loc*gzzz(i,j,k))) Gmy_Res(i,j,k) = Gamy(i,j,k) - ( & gupxx_loc*(gupxy_loc*gxxx(i,j,k)+gupyy_loc*gxyx(i,j,k)+gupyz_loc*gxzx(i,j,k)) + & gupxy_loc*(gupxy_loc*gxyx(i,j,k)+gupyy_loc*gyyx(i,j,k)+gupyz_loc*gyzx(i,j,k)) + & gupxz_loc*(gupxy_loc*gxzx(i,j,k)+gupyy_loc*gyzx(i,j,k)+gupyz_loc*gzzx(i,j,k)) + & gupxy_loc*(gupxy_loc*gxxy(i,j,k)+gupyy_loc*gxyy(i,j,k)+gupyz_loc*gxzy(i,j,k)) + & gupyy_loc*(gupxy_loc*gxyy(i,j,k)+gupyy_loc*gyyy(i,j,k)+gupyz_loc*gyzy(i,j,k)) + & gupyz_loc*(gupxy_loc*gxzy(i,j,k)+gupyy_loc*gyzy(i,j,k)+gupyz_loc*gzzy(i,j,k)) + & gupxy_loc*(gupxz_loc*gxxz(i,j,k)+gupyz_loc*gxyz(i,j,k)+gupzz_loc*gxzz(i,j,k)) + & gupyy_loc*(gupxz_loc*gxyz(i,j,k)+gupyz_loc*gyyz(i,j,k)+gupzz_loc*gyzz(i,j,k)) + & gupyz_loc*(gupxz_loc*gxzz(i,j,k)+gupyz_loc*gyzz(i,j,k)+gupzz_loc*gzzz(i,j,k))) Gmz_Res(i,j,k) = Gamz(i,j,k) - ( & gupxx_loc*(gupxz_loc*gxxx(i,j,k)+gupyz_loc*gxyx(i,j,k)+gupzz_loc*gxzx(i,j,k)) + & gupxy_loc*(gupxz_loc*gxyx(i,j,k)+gupyz_loc*gyyx(i,j,k)+gupzz_loc*gyzx(i,j,k)) + & gupxz_loc*(gupxz_loc*gxzx(i,j,k)+gupyz_loc*gyzx(i,j,k)+gupzz_loc*gzzx(i,j,k)) + & gupxy_loc*(gupxz_loc*gxxy(i,j,k)+gupyz_loc*gxyy(i,j,k)+gupzz_loc*gxzy(i,j,k)) + & gupyy_loc*(gupxz_loc*gxyy(i,j,k)+gupyz_loc*gyyy(i,j,k)+gupzz_loc*gyzy(i,j,k)) + & gupyz_loc*(gupxz_loc*gxzy(i,j,k)+gupyz_loc*gyzy(i,j,k)+gupzz_loc*gzzy(i,j,k)) + & gupxz_loc*(gupxz_loc*gxxz(i,j,k)+gupyz_loc*gxyz(i,j,k)+gupzz_loc*gxzz(i,j,k)) + & gupyz_loc*(gupxz_loc*gxyz(i,j,k)+gupyz_loc*gyyz(i,j,k)+gupzz_loc*gyzz(i,j,k)) + & gupzz_loc*(gupxz_loc*gxzz(i,j,k)+gupyz_loc*gyzz(i,j,k)+gupzz_loc*gzzz(i,j,k))) endif Gamxxx(i,j,k)=HALF*( gupxx_loc*gxxx(i,j,k) + gupxy_loc*(TWO*gxyx(i,j,k) - gxxy(i,j,k)) + gupxz_loc*(TWO*gxzx(i,j,k) - gxxz(i,j,k))) Gamyxx(i,j,k)=HALF*( gupxy_loc*gxxx(i,j,k) + gupyy_loc*(TWO*gxyx(i,j,k) - gxxy(i,j,k)) + gupyz_loc*(TWO*gxzx(i,j,k) - gxxz(i,j,k))) Gamzxx(i,j,k)=HALF*( gupxz_loc*gxxx(i,j,k) + gupyz_loc*(TWO*gxyx(i,j,k) - gxxy(i,j,k)) + gupzz_loc*(TWO*gxzx(i,j,k) - gxxz(i,j,k))) Gamxyy(i,j,k)=HALF*( gupxx_loc*(TWO*gxyy(i,j,k) - gyyx(i,j,k)) + gupxy_loc*gyyy(i,j,k) + gupxz_loc*(TWO*gyzy(i,j,k) - gyyz(i,j,k))) Gamyyy(i,j,k)=HALF*( gupxy_loc*(TWO*gxyy(i,j,k) - gyyx(i,j,k)) + gupyy_loc*gyyy(i,j,k) + gupyz_loc*(TWO*gyzy(i,j,k) - gyyz(i,j,k))) Gamzyy(i,j,k)=HALF*( gupxz_loc*(TWO*gxyy(i,j,k) - gyyx(i,j,k)) + gupyz_loc*gyyy(i,j,k) + gupzz_loc*(TWO*gyzy(i,j,k) - gyyz(i,j,k))) Gamxzz(i,j,k)=HALF*( gupxx_loc*(TWO*gxzz(i,j,k) - gzzx(i,j,k)) + gupxy_loc*(TWO*gyzz(i,j,k) - gzzy(i,j,k)) + gupxz_loc*gzzz(i,j,k)) Gamyzz(i,j,k)=HALF*( gupxy_loc*(TWO*gxzz(i,j,k) - gzzx(i,j,k)) + gupyy_loc*(TWO*gyzz(i,j,k) - gzzy(i,j,k)) + gupyz_loc*gzzz(i,j,k)) Gamzzz(i,j,k)=HALF*( gupxz_loc*(TWO*gxzz(i,j,k) - gzzx(i,j,k)) + gupyz_loc*(TWO*gyzz(i,j,k) - gzzy(i,j,k)) + gupzz_loc*gzzz(i,j,k)) Gamxxy(i,j,k)=HALF*( gupxx_loc*gxxy(i,j,k) + gupxy_loc*gyyx(i,j,k) + gupxz_loc*(gxzy(i,j,k) + gyzx(i,j,k) - gxyz(i,j,k)) ) Gamyxy(i,j,k)=HALF*( gupxy_loc*gxxy(i,j,k) + gupyy_loc*gyyx(i,j,k) + gupyz_loc*(gxzy(i,j,k) + gyzx(i,j,k) - gxyz(i,j,k)) ) Gamzxy(i,j,k)=HALF*( gupxz_loc*gxxy(i,j,k) + gupyz_loc*gyyx(i,j,k) + gupzz_loc*(gxzy(i,j,k) + gyzx(i,j,k) - gxyz(i,j,k)) ) Gamxxz(i,j,k)=HALF*( gupxx_loc*gxxz(i,j,k) + gupxy_loc*(gxyz(i,j,k) + gyzx(i,j,k) - gxzy(i,j,k)) + gupxz_loc*gzzx(i,j,k) ) Gamyxz(i,j,k)=HALF*( gupxy_loc*gxxz(i,j,k) + gupyy_loc*(gxyz(i,j,k) + gyzx(i,j,k) - gxzy(i,j,k)) + gupyz_loc*gzzx(i,j,k) ) Gamzxz(i,j,k)=HALF*( gupxz_loc*gxxz(i,j,k) + gupyz_loc*(gxyz(i,j,k) + gyzx(i,j,k) - gxzy(i,j,k)) + gupzz_loc*gzzx(i,j,k) ) Gamxyz(i,j,k)=HALF*( gupxx_loc*(gxyz(i,j,k) + gxzy(i,j,k) - gyzx(i,j,k)) + gupxy_loc*gyyz(i,j,k) + gupxz_loc*gzzy(i,j,k) ) Gamyyz(i,j,k)=HALF*( gupxy_loc*(gxyz(i,j,k) + gxzy(i,j,k) - gyzx(i,j,k)) + gupyy_loc*gyyz(i,j,k) + gupyz_loc*gzzy(i,j,k) ) Gamzyz(i,j,k)=HALF*( gupxz_loc*(gxyz(i,j,k) + gxzy(i,j,k) - gyzx(i,j,k)) + gupyz_loc*gyyz(i,j,k) + gupzz_loc*gzzy(i,j,k) ) enddo enddo enddo ! Raise indices of \tilde A_{ij} and store in R_ij ! Right hand side for Gam^i without shift terms... call fderivs(ex,Lap,Lapx,Lapy,Lapz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) call fderivs(ex,trK,Kx,Ky,Kz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev) do k=1,ex(3) do j=1,ex(2) do i=1,ex(1) gupxx_loc = gupxx(i,j,k) gupxy_loc = gupxy(i,j,k) gupxz_loc = gupxz(i,j,k) gupyy_loc = gupyy(i,j,k) gupyz_loc = gupyz(i,j,k) gupzz_loc = gupzz(i,j,k) Rxx_loc = gupxx_loc * gupxx_loc * Axx(i,j,k) + gupxy_loc * gupxy_loc * Ayy(i,j,k) + gupxz_loc * gupxz_loc * Azz(i,j,k) + & TWO * (gupxx_loc * gupxy_loc * Axy(i,j,k) + gupxx_loc * gupxz_loc * Axz(i,j,k) + gupxy_loc * gupxz_loc * Ayz(i,j,k)) Ryy_loc = gupxy_loc * gupxy_loc * Axx(i,j,k) + gupyy_loc * gupyy_loc * Ayy(i,j,k) + gupyz_loc * gupyz_loc * Azz(i,j,k) + & TWO * (gupxy_loc * gupyy_loc * Axy(i,j,k) + gupxy_loc * gupyz_loc * Axz(i,j,k) + gupyy_loc * gupyz_loc * Ayz(i,j,k)) Rzz_loc = gupxz_loc * gupxz_loc * Axx(i,j,k) + gupyz_loc * gupyz_loc * Ayy(i,j,k) + gupzz_loc * gupzz_loc * Azz(i,j,k) + & TWO * (gupxz_loc * gupyz_loc * Axy(i,j,k) + gupxz_loc * gupzz_loc * Axz(i,j,k) + gupyz_loc * gupzz_loc * Ayz(i,j,k)) Rxy_loc = gupxx_loc * gupxy_loc * Axx(i,j,k) + gupxy_loc * gupyy_loc * Ayy(i,j,k) + gupxz_loc * gupyz_loc * Azz(i,j,k) + & (gupxx_loc * gupyy_loc + gupxy_loc * gupxy_loc) * Axy(i,j,k) + & (gupxx_loc * gupyz_loc + gupxz_loc * gupxy_loc) * Axz(i,j,k) + & (gupxy_loc * gupyz_loc + gupxz_loc * gupyy_loc) * Ayz(i,j,k) Rxz_loc = gupxx_loc * gupxz_loc * Axx(i,j,k) + gupxy_loc * gupyz_loc * Ayy(i,j,k) + gupxz_loc * gupzz_loc * Azz(i,j,k) + & (gupxx_loc * gupyz_loc + gupxy_loc * gupxz_loc) * Axy(i,j,k) + & (gupxx_loc * gupzz_loc + gupxz_loc * gupxz_loc) * Axz(i,j,k) + & (gupxy_loc * gupzz_loc + gupxz_loc * gupyz_loc) * Ayz(i,j,k) Ryz_loc = gupxy_loc * gupxz_loc * Axx(i,j,k) + gupyy_loc * gupyz_loc * Ayy(i,j,k) + gupyz_loc * gupzz_loc * Azz(i,j,k) + & (gupxy_loc * gupyz_loc + gupyy_loc * gupxz_loc) * Axy(i,j,k) + & (gupxy_loc * gupzz_loc + gupyz_loc * gupxz_loc) * Axz(i,j,k) + & (gupyy_loc * gupzz_loc + gupyz_loc * gupyz_loc) * Ayz(i,j,k) Rxx(i,j,k) = Rxx_loc Ryy(i,j,k) = Ryy_loc Rzz(i,j,k) = Rzz_loc Rxy(i,j,k) = Rxy_loc Rxz(i,j,k) = Rxz_loc Ryz(i,j,k) = Ryz_loc Gamx_rhs(i,j,k) = - TWO * (Lapx(i,j,k) * Rxx_loc + Lapy(i,j,k) * Rxy_loc + Lapz(i,j,k) * Rxz_loc) + & TWO * alpn1(i,j,k) * ( & -F3o2/chin1(i,j,k) * (chix(i,j,k) * Rxx_loc + chiy(i,j,k) * Rxy_loc + chiz(i,j,k) * Rxz_loc) - & gupxx_loc * (F2o3 * Kx(i,j,k) + EIGHT * PI * Sx(i,j,k)) - & gupxy_loc * (F2o3 * Ky(i,j,k) + EIGHT * PI * Sy(i,j,k)) - & gupxz_loc * (F2o3 * Kz(i,j,k) + EIGHT * PI * Sz(i,j,k)) + & Gamxxx(i,j,k) * Rxx_loc + Gamxyy(i,j,k) * Ryy_loc + Gamxzz(i,j,k) * Rzz_loc + & TWO * (Gamxxy(i,j,k) * Rxy_loc + Gamxxz(i,j,k) * Rxz_loc + Gamxyz(i,j,k) * Ryz_loc)) Gamy_rhs(i,j,k) = - TWO * (Lapx(i,j,k) * Rxy_loc + Lapy(i,j,k) * Ryy_loc + Lapz(i,j,k) * Ryz_loc) + & TWO * alpn1(i,j,k) * ( & -F3o2/chin1(i,j,k) * (chix(i,j,k) * Rxy_loc + chiy(i,j,k) * Ryy_loc + chiz(i,j,k) * Ryz_loc) - & gupxy_loc * (F2o3 * Kx(i,j,k) + EIGHT * PI * Sx(i,j,k)) - & gupyy_loc * (F2o3 * Ky(i,j,k) + EIGHT * PI * Sy(i,j,k)) - & gupyz_loc * (F2o3 * Kz(i,j,k) + EIGHT * PI * Sz(i,j,k)) + & Gamyxx(i,j,k) * Rxx_loc + Gamyyy(i,j,k) * Ryy_loc + Gamyzz(i,j,k) * Rzz_loc + & TWO * (Gamyxy(i,j,k) * Rxy_loc + Gamyxz(i,j,k) * Rxz_loc + Gamyyz(i,j,k) * Ryz_loc)) Gamz_rhs(i,j,k) = - TWO * (Lapx(i,j,k) * Rxz_loc + Lapy(i,j,k) * Ryz_loc + Lapz(i,j,k) * Rzz_loc) + & TWO * alpn1(i,j,k) * ( & -F3o2/chin1(i,j,k) * (chix(i,j,k) * Rxz_loc + chiy(i,j,k) * Ryz_loc + chiz(i,j,k) * Rzz_loc) - & gupxz_loc * (F2o3 * Kx(i,j,k) + EIGHT * PI * Sx(i,j,k)) - & gupyz_loc * (F2o3 * Ky(i,j,k) + EIGHT * PI * Sy(i,j,k)) - & gupzz_loc * (F2o3 * Kz(i,j,k) + EIGHT * PI * Sz(i,j,k)) + & Gamzxx(i,j,k) * Rxx_loc + Gamzyy(i,j,k) * Ryy_loc + Gamzzz(i,j,k) * Rzz_loc + & TWO * (Gamzxy(i,j,k) * Rxy_loc + Gamzxz(i,j,k) * Rxz_loc + Gamzyz(i,j,k) * Ryz_loc)) enddo enddo enddo call fdderivs(ex,betax,gxxx,gxyx,gxzx,gyyx,gyzx,gzzx,& X,Y,Z,ANTI,SYM, SYM ,Symmetry,Lev) call fdderivs(ex,betay,gxxy,gxyy,gxzy,gyyy,gyzy,gzzy,& X,Y,Z,SYM ,ANTI,SYM ,Symmetry,Lev) call fdderivs(ex,betaz,gxxz,gxyz,gxzz,gyyz,gyzz,gzzz,& X,Y,Z,SYM ,SYM, ANTI,Symmetry,Lev) call fderivs(ex,Gamx,Gamxx,Gamxy,Gamxz,X,Y,Z,ANTI,SYM ,SYM ,Symmetry,Lev) call fderivs(ex,Gamy,Gamyx,Gamyy,Gamyz,X,Y,Z,SYM ,ANTI,SYM ,Symmetry,Lev) call fderivs(ex,Gamz,Gamzx,Gamzy,Gamzz,X,Y,Z,SYM ,SYM ,ANTI,Symmetry,Lev) do k=1,ex(3) do j=1,ex(2) do i=1,ex(1) divb_loc = div_beta(i,j,k) fxx_loc = gxxx(i,j,k) + gxyy(i,j,k) + gxzz(i,j,k) fxy_loc = gxyx(i,j,k) + gyyy(i,j,k) + gyzz(i,j,k) fxz_loc = gxzx(i,j,k) + gyzy(i,j,k) + gzzz(i,j,k) gupxx_loc = gupxx(i,j,k) gupxy_loc = gupxy(i,j,k) gupxz_loc = gupxz(i,j,k) gupyy_loc = gupyy(i,j,k) gupyz_loc = gupyz(i,j,k) gupzz_loc = gupzz(i,j,k) Gamxa_loc = gupxx_loc * Gamxxx(i,j,k) + gupyy_loc * Gamxyy(i,j,k) + gupzz_loc * Gamxzz(i,j,k) + & TWO * (gupxy_loc * Gamxxy(i,j,k) + gupxz_loc * Gamxxz(i,j,k) + gupyz_loc * Gamxyz(i,j,k)) Gamya_loc = gupxx_loc * Gamyxx(i,j,k) + gupyy_loc * Gamyyy(i,j,k) + gupzz_loc * Gamyzz(i,j,k) + & TWO * (gupxy_loc * Gamyxy(i,j,k) + gupxz_loc * Gamyxz(i,j,k) + gupyz_loc * Gamyyz(i,j,k)) Gamza_loc = gupxx_loc * Gamzxx(i,j,k) + gupyy_loc * Gamzyy(i,j,k) + gupzz_loc * Gamzzz(i,j,k) + & TWO * (gupxy_loc * Gamzxy(i,j,k) + gupxz_loc * Gamzxz(i,j,k) + gupyz_loc * Gamzyz(i,j,k)) Gamxa(i,j,k) = Gamxa_loc Gamya(i,j,k) = Gamya_loc Gamza(i,j,k) = Gamza_loc Gamx_rhs(i,j,k) = Gamx_rhs(i,j,k) + F2o3 * Gamxa_loc * divb_loc - & Gamxa_loc * betaxx(i,j,k) - Gamya_loc * betaxy(i,j,k) - Gamza_loc * betaxz(i,j,k) + & F1o3 * (gupxx_loc * fxx_loc + gupxy_loc * fxy_loc + gupxz_loc * fxz_loc) + & gupxx_loc * gxxx(i,j,k) + gupyy_loc * gyyx(i,j,k) + gupzz_loc * gzzx(i,j,k) + & TWO * (gupxy_loc * gxyx(i,j,k) + gupxz_loc * gxzx(i,j,k) + gupyz_loc * gyzx(i,j,k)) Gamy_rhs(i,j,k) = Gamy_rhs(i,j,k) + F2o3 * Gamya_loc * divb_loc - & Gamxa_loc * betayx(i,j,k) - Gamya_loc * betayy(i,j,k) - Gamza_loc * betayz(i,j,k) + & F1o3 * (gupxy_loc * fxx_loc + gupyy_loc * fxy_loc + gupyz_loc * fxz_loc) + & gupxx_loc * gxxy(i,j,k) + gupyy_loc * gyyy(i,j,k) + gupzz_loc * gzzy(i,j,k) + & TWO * (gupxy_loc * gxyy(i,j,k) + gupxz_loc * gxzy(i,j,k) + gupyz_loc * gyzy(i,j,k)) Gamz_rhs(i,j,k) = Gamz_rhs(i,j,k) + F2o3 * Gamza_loc * divb_loc - & Gamxa_loc * betazx(i,j,k) - Gamya_loc * betazy(i,j,k) - Gamza_loc * betazz(i,j,k) + & F1o3 * (gupxz_loc * fxx_loc + gupyz_loc * fxy_loc + gupzz_loc * fxz_loc) + & gupxx_loc * gxxz(i,j,k) + gupyy_loc * gyyz(i,j,k) + gupzz_loc * gzzz(i,j,k) + & TWO * (gupxy_loc * gxyz(i,j,k) + gupxz_loc * gxzz(i,j,k) + gupyz_loc * gyzz(i,j,k)) enddo enddo enddo !first kind of connection stored in gij,k gxxx = gxx * Gamxxx + gxy * Gamyxx + gxz * Gamzxx gxyx = gxx * Gamxxy + gxy * Gamyxy + gxz * Gamzxy gxzx = gxx * Gamxxz + gxy * Gamyxz + gxz * Gamzxz gyyx = gxx * Gamxyy + gxy * Gamyyy + gxz * Gamzyy gyzx = gxx * Gamxyz + gxy * Gamyyz + gxz * Gamzyz gzzx = gxx * Gamxzz + gxy * Gamyzz + gxz * Gamzzz gxxy = gxy * Gamxxx + gyy * Gamyxx + gyz * Gamzxx gxyy = gxy * Gamxxy + gyy * Gamyxy + gyz * Gamzxy gxzy = gxy * Gamxxz + gyy * Gamyxz + gyz * Gamzxz gyyy = gxy * Gamxyy + gyy * Gamyyy + gyz * Gamzyy gyzy = gxy * Gamxyz + gyy * Gamyyz + gyz * Gamzyz gzzy = gxy * Gamxzz + gyy * Gamyzz + gyz * Gamzzz gxxz = gxz * Gamxxx + gyz * Gamyxx + gzz * Gamzxx gxyz = gxz * Gamxxy + gyz * Gamyxy + gzz * Gamzxy gxzz = gxz * Gamxxz + gyz * Gamyxz + gzz * Gamzxz gyyz = gxz * Gamxyy + gyz * Gamyyy + gzz * Gamzyy gyzz = gxz * Gamxyz + gyz * Gamyyz + gzz * Gamzyz gzzz = gxz * Gamxzz + gyz * Gamyzz + gzz * Gamzzz !compute Ricci tensor for tilted metric call fdderivs(ex,dxx,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,SYM ,SYM ,symmetry,Lev) Rxx = gupxx * fxx + gupyy * fyy + gupzz * fzz + & ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) * TWO call fdderivs(ex,dyy,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,SYM ,SYM ,symmetry,Lev) Ryy = gupxx * fxx + gupyy * fyy + gupzz * fzz + & ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) * TWO call fdderivs(ex,dzz,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,SYM ,SYM ,symmetry,Lev) Rzz = gupxx * fxx + gupyy * fyy + gupzz * fzz + & ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) * TWO call fdderivs(ex,gxy,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,ANTI, ANTI,SYM ,symmetry,Lev) Rxy = gupxx * fxx + gupyy * fyy + gupzz * fzz + & ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) * TWO call fdderivs(ex,gxz,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,ANTI ,SYM ,ANTI,symmetry,Lev) Rxz = gupxx * fxx + gupyy * fyy + gupzz * fzz + & ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) * TWO call fdderivs(ex,gyz,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM ,ANTI ,ANTI,symmetry,Lev) Ryz = gupxx * fxx + gupyy * fyy + gupzz * fzz + & ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) * TWO Rxx = - HALF * Rxx + & gxx * Gamxx+ gxy * Gamyx + gxz * Gamzx + & Gamxa * gxxx + Gamya * gxyx + Gamza * gxzx + & gupxx *( & TWO*(Gamxxx * gxxx + Gamyxx * gxyx + Gamzxx * gxzx) + & Gamxxx * gxxx + Gamyxx * gxxy + Gamzxx * gxxz )+ & gupxy *( & TWO*(Gamxxx * gxyx + Gamyxx * gyyx + Gamzxx * gyzx + & Gamxxy * gxxx + Gamyxy * gxyx + Gamzxy * gxzx) + & Gamxxy * gxxx + Gamyxy * gxxy + Gamzxy * gxxz + & Gamxxx * gxyx + Gamyxx * gxyy + Gamzxx * gxyz )+ & gupxz *( & TWO*(Gamxxx * gxzx + Gamyxx * gyzx + Gamzxx * gzzx + & Gamxxz * gxxx + Gamyxz * gxyx + Gamzxz * gxzx) + & Gamxxz * gxxx + Gamyxz * gxxy + Gamzxz * gxxz + & Gamxxx * gxzx + Gamyxx * gxzy + Gamzxx * gxzz )+ & gupyy *( & TWO*(Gamxxy * gxyx + Gamyxy * gyyx + Gamzxy * gyzx) + & Gamxxy * gxyx + Gamyxy * gxyy + Gamzxy * gxyz )+ & gupyz *( & TWO*(Gamxxy * gxzx + Gamyxy * gyzx + Gamzxy * gzzx + & Gamxxz * gxyx + Gamyxz * gyyx + Gamzxz * gyzx) + & Gamxxz * gxyx + Gamyxz * gxyy + Gamzxz * gxyz + & Gamxxy * gxzx + Gamyxy * gxzy + Gamzxy * gxzz )+ & gupzz *( & TWO*(Gamxxz * gxzx + Gamyxz * gyzx + Gamzxz * gzzx) + & Gamxxz * gxzx + Gamyxz * gxzy + Gamzxz * gxzz ) Ryy = - HALF * Ryy + & gxy * Gamxy+ gyy * Gamyy + gyz * Gamzy + & Gamxa * gxyy + Gamya * gyyy + Gamza * gyzy + & gupxx *( & TWO*(Gamxxy * gxxy + Gamyxy * gxyy + Gamzxy * gxzy) + & Gamxxy * gxyx + Gamyxy * gxyy + Gamzxy * gxyz )+ & gupxy *( & TWO*(Gamxxy * gxyy + Gamyxy * gyyy + Gamzxy * gyzy + & Gamxyy * gxxy + Gamyyy * gxyy + Gamzyy * gxzy) + & Gamxyy * gxyx + Gamyyy * gxyy + Gamzyy * gxyz + & Gamxxy * gyyx + Gamyxy * gyyy + Gamzxy * gyyz )+ & gupxz *( & TWO*(Gamxxy * gxzy + Gamyxy * gyzy + Gamzxy * gzzy + & Gamxyz * gxxy + Gamyyz * gxyy + Gamzyz * gxzy) + & Gamxyz * gxyx + Gamyyz * gxyy + Gamzyz * gxyz + & Gamxxy * gyzx + Gamyxy * gyzy + Gamzxy * gyzz )+ & gupyy *( & TWO*(Gamxyy * gxyy + Gamyyy * gyyy + Gamzyy * gyzy) + & Gamxyy * gyyx + Gamyyy * gyyy + Gamzyy * gyyz )+ & gupyz *( & TWO*(Gamxyy * gxzy + Gamyyy * gyzy + Gamzyy * gzzy + & Gamxyz * gxyy + Gamyyz * gyyy + Gamzyz * gyzy) + & Gamxyz * gyyx + Gamyyz * gyyy + Gamzyz * gyyz + & Gamxyy * gyzx + Gamyyy * gyzy + Gamzyy * gyzz )+ & gupzz *( & TWO*(Gamxyz * gxzy + Gamyyz * gyzy + Gamzyz * gzzy) + & Gamxyz * gyzx + Gamyyz * gyzy + Gamzyz * gyzz ) Rzz = - HALF * Rzz + & gxz * Gamxz+ gyz * Gamyz + gzz * Gamzz + & Gamxa * gxzz + Gamya * gyzz + Gamza * gzzz + & gupxx *( & TWO*(Gamxxz * gxxz + Gamyxz * gxyz + Gamzxz * gxzz) + & Gamxxz * gxzx + Gamyxz * gxzy + Gamzxz * gxzz )+ & gupxy *( & TWO*(Gamxxz * gxyz + Gamyxz * gyyz + Gamzxz * gyzz + & Gamxyz * gxxz + Gamyyz * gxyz + Gamzyz * gxzz) + & Gamxyz * gxzx + Gamyyz * gxzy + Gamzyz * gxzz + & Gamxxz * gyzx + Gamyxz * gyzy + Gamzxz * gyzz )+ & gupxz *( & TWO*(Gamxxz * gxzz + Gamyxz * gyzz + Gamzxz * gzzz + & Gamxzz * gxxz + Gamyzz * gxyz + Gamzzz * gxzz) + & Gamxzz * gxzx + Gamyzz * gxzy + Gamzzz * gxzz + & Gamxxz * gzzx + Gamyxz * gzzy + Gamzxz * gzzz )+ & gupyy *( & TWO*(Gamxyz * gxyz + Gamyyz * gyyz + Gamzyz * gyzz) + & Gamxyz * gyzx + Gamyyz * gyzy + Gamzyz * gyzz )+ & gupyz *( & TWO*(Gamxyz * gxzz + Gamyyz * gyzz + Gamzyz * gzzz + & Gamxzz * gxyz + Gamyzz * gyyz + Gamzzz * gyzz) + & Gamxzz * gyzx + Gamyzz * gyzy + Gamzzz * gyzz + & Gamxyz * gzzx + Gamyyz * gzzy + Gamzyz * gzzz )+ & gupzz *( & TWO*(Gamxzz * gxzz + Gamyzz * gyzz + Gamzzz * gzzz) + & Gamxzz * gzzx + Gamyzz * gzzy + Gamzzz * gzzz ) Rxy = HALF*( - Rxy + & gxx * Gamxy + gxy * Gamyy + gxz * Gamzy + & gxy * Gamxx + gyy * Gamyx + gyz * Gamzx + & Gamxa * gxyx + Gamya * gyyx + Gamza * gyzx + & Gamxa * gxxy + Gamya * gxyy + Gamza * gxzy )+ & gupxx *( & Gamxxx * gxxy + Gamyxx * gxyy + Gamzxx * gxzy + & Gamxxy * gxxx + Gamyxy * gxyx + Gamzxy * gxzx + & Gamxxx * gxyx + Gamyxx * gxyy + Gamzxx * gxyz )+ & gupxy *( & Gamxxx * gxyy + Gamyxx * gyyy + Gamzxx * gyzy + & Gamxxy * gxyx + Gamyxy * gyyx + Gamzxy * gyzx + & Gamxxy * gxyx + Gamyxy * gxyy + Gamzxy * gxyz + & Gamxxy * gxxy + Gamyxy * gxyy + Gamzxy * gxzy + & Gamxyy * gxxx + Gamyyy * gxyx + Gamzyy * gxzx + & Gamxxx * gyyx + Gamyxx * gyyy + Gamzxx * gyyz )+ & gupxz *( & Gamxxx * gxzy + Gamyxx * gyzy + Gamzxx * gzzy + & Gamxxy * gxzx + Gamyxy * gyzx + Gamzxy * gzzx + & Gamxxz * gxyx + Gamyxz * gxyy + Gamzxz * gxyz + & Gamxxz * gxxy + Gamyxz * gxyy + Gamzxz * gxzy + & Gamxyz * gxxx + Gamyyz * gxyx + Gamzyz * gxzx + & Gamxxx * gyzx + Gamyxx * gyzy + Gamzxx * gyzz )+ & gupyy *( & Gamxxy * gxyy + Gamyxy * gyyy + Gamzxy * gyzy + & Gamxyy * gxyx + Gamyyy * gyyx + Gamzyy * gyzx + & Gamxxy * gyyx + Gamyxy * gyyy + Gamzxy * gyyz )+ & gupyz *( & Gamxxy * gxzy + Gamyxy * gyzy + Gamzxy * gzzy + & Gamxyy * gxzx + Gamyyy * gyzx + Gamzyy * gzzx + & Gamxxz * gyyx + Gamyxz * gyyy + Gamzxz * gyyz + & Gamxxz * gxyy + Gamyxz * gyyy + Gamzxz * gyzy + & Gamxyz * gxyx + Gamyyz * gyyx + Gamzyz * gyzx + & Gamxxy * gyzx + Gamyxy * gyzy + Gamzxy * gyzz )+ & gupzz *( & Gamxxz * gxzy + Gamyxz * gyzy + Gamzxz * gzzy + & Gamxyz * gxzx + Gamyyz * gyzx + Gamzyz * gzzx + & Gamxxz * gyzx + Gamyxz * gyzy + Gamzxz * gyzz ) Rxz = HALF*( - Rxz + & gxx * Gamxz + gxy * Gamyz + gxz * Gamzz + & gxz * Gamxx + gyz * Gamyx + gzz * Gamzx + & Gamxa * gxzx + Gamya * gyzx + Gamza * gzzx + & Gamxa * gxxz + Gamya * gxyz + Gamza * gxzz )+ & gupxx *( & Gamxxx * gxxz + Gamyxx * gxyz + Gamzxx * gxzz + & Gamxxz * gxxx + Gamyxz * gxyx + Gamzxz * gxzx + & Gamxxx * gxzx + Gamyxx * gxzy + Gamzxx * gxzz )+ & gupxy *( & Gamxxx * gxyz + Gamyxx * gyyz + Gamzxx * gyzz + & Gamxxz * gxyx + Gamyxz * gyyx + Gamzxz * gyzx + & Gamxxy * gxzx + Gamyxy * gxzy + Gamzxy * gxzz + & Gamxxy * gxxz + Gamyxy * gxyz + Gamzxy * gxzz + & Gamxyz * gxxx + Gamyyz * gxyx + Gamzyz * gxzx + & Gamxxx * gyzx + Gamyxx * gyzy + Gamzxx * gyzz )+ & gupxz *( & Gamxxx * gxzz + Gamyxx * gyzz + Gamzxx * gzzz + & Gamxxz * gxzx + Gamyxz * gyzx + Gamzxz * gzzx + & Gamxxz * gxzx + Gamyxz * gxzy + Gamzxz * gxzz + & Gamxxz * gxxz + Gamyxz * gxyz + Gamzxz * gxzz + & Gamxzz * gxxx + Gamyzz * gxyx + Gamzzz * gxzx + & Gamxxx * gzzx + Gamyxx * gzzy + Gamzxx * gzzz )+ & gupyy *( & Gamxxy * gxyz + Gamyxy * gyyz + Gamzxy * gyzz + & Gamxyz * gxyx + Gamyyz * gyyx + Gamzyz * gyzx + & Gamxxy * gyzx + Gamyxy * gyzy + Gamzxy * gyzz )+ & gupyz *( & Gamxxy * gxzz + Gamyxy * gyzz + Gamzxy * gzzz + & Gamxyz * gxzx + Gamyyz * gyzx + Gamzyz * gzzx + & Gamxxz * gyzx + Gamyxz * gyzy + Gamzxz * gyzz + & Gamxxz * gxyz + Gamyxz * gyyz + Gamzxz * gyzz + & Gamxzz * gxyx + Gamyzz * gyyx + Gamzzz * gyzx + & Gamxxy * gzzx + Gamyxy * gzzy + Gamzxy * gzzz )+ & gupzz *( & Gamxxz * gxzz + Gamyxz * gyzz + Gamzxz * gzzz + & Gamxzz * gxzx + Gamyzz * gyzx + Gamzzz * gzzx + & Gamxxz * gzzx + Gamyxz * gzzy + Gamzxz * gzzz ) Ryz = HALF*( - Ryz + & gxy * Gamxz + gyy * Gamyz + gyz * Gamzz + & gxz * Gamxy + gyz * Gamyy + gzz * Gamzy + & Gamxa * gxzy + Gamya * gyzy + Gamza * gzzy + & Gamxa * gxyz + Gamya * gyyz + Gamza * gyzz )+ & gupxx *( & Gamxxy * gxxz + Gamyxy * gxyz + Gamzxy * gxzz + & Gamxxz * gxxy + Gamyxz * gxyy + Gamzxz * gxzy + & Gamxxy * gxzx + Gamyxy * gxzy + Gamzxy * gxzz )+ & gupxy *( & Gamxxy * gxyz + Gamyxy * gyyz + Gamzxy * gyzz + & Gamxxz * gxyy + Gamyxz * gyyy + Gamzxz * gyzy + & Gamxyy * gxzx + Gamyyy * gxzy + Gamzyy * gxzz + & Gamxyy * gxxz + Gamyyy * gxyz + Gamzyy * gxzz + & Gamxyz * gxxy + Gamyyz * gxyy + Gamzyz * gxzy + & Gamxxy * gyzx + Gamyxy * gyzy + Gamzxy * gyzz )+ & gupxz *( & Gamxxy * gxzz + Gamyxy * gyzz + Gamzxy * gzzz + & Gamxxz * gxzy + Gamyxz * gyzy + Gamzxz * gzzy + & Gamxyz * gxzx + Gamyyz * gxzy + Gamzyz * gxzz + & Gamxyz * gxxz + Gamyyz * gxyz + Gamzyz * gxzz + & Gamxzz * gxxy + Gamyzz * gxyy + Gamzzz * gxzy + & Gamxxy * gzzx + Gamyxy * gzzy + Gamzxy * gzzz )+ & gupyy *( & Gamxyy * gxyz + Gamyyy * gyyz + Gamzyy * gyzz + & Gamxyz * gxyy + Gamyyz * gyyy + Gamzyz * gyzy + & Gamxyy * gyzx + Gamyyy * gyzy + Gamzyy * gyzz )+ & gupyz *( & Gamxyy * gxzz + Gamyyy * gyzz + Gamzyy * gzzz + & Gamxyz * gxzy + Gamyyz * gyzy + Gamzyz * gzzy + & Gamxyz * gyzx + Gamyyz * gyzy + Gamzyz * gyzz + & Gamxyz * gxyz + Gamyyz * gyyz + Gamzyz * gyzz + & Gamxzz * gxyy + Gamyzz * gyyy + Gamzzz * gyzy + & Gamxyy * gzzx + Gamyyy * gzzy + Gamzyy * gzzz )+ & gupzz *( & Gamxyz * gxzz + Gamyyz * gyzz + Gamzyz * gzzz + & Gamxzz * gxzy + Gamyzz * gyzy + Gamzzz * gzzy + & Gamxyz * gzzx + Gamyyz * gzzy + Gamzyz * gzzz ) !covariant second derivative of chi respect to tilted metric call fdderivs(ex,chi,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) fxx = fxx - Gamxxx * chix - Gamyxx * chiy - Gamzxx * chiz fxy = fxy - Gamxxy * chix - Gamyxy * chiy - Gamzxy * chiz fxz = fxz - Gamxxz * chix - Gamyxz * chiy - Gamzxz * chiz fyy = fyy - Gamxyy * chix - Gamyyy * chiy - Gamzyy * chiz fyz = fyz - Gamxyz * chix - Gamyyz * chiy - Gamzyz * chiz fzz = fzz - Gamxzz * chix - Gamyzz * chiy - Gamzzz * chiz ! Store D^l D_l chi - 3/(2*chi) D^l chi D_l chi in f f = gupxx * ( fxx - F3o2/chin1 * chix * chix ) + & gupyy * ( fyy - F3o2/chin1 * chiy * chiy ) + & gupzz * ( fzz - F3o2/chin1 * chiz * chiz ) + & TWO * gupxy * ( fxy - F3o2/chin1 * chix * chiy ) + & TWO * gupxz * ( fxz - F3o2/chin1 * chix * chiz ) + & TWO * gupyz * ( fyz - F3o2/chin1 * chiy * chiz ) ! Add chi part to Ricci tensor: Rxx = Rxx + (fxx - chix*chix/chin1/TWO + gxx * f)/chin1/TWO Ryy = Ryy + (fyy - chiy*chiy/chin1/TWO + gyy * f)/chin1/TWO Rzz = Rzz + (fzz - chiz*chiz/chin1/TWO + gzz * f)/chin1/TWO Rxy = Rxy + (fxy - chix*chiy/chin1/TWO + gxy * f)/chin1/TWO Rxz = Rxz + (fxz - chix*chiz/chin1/TWO + gxz * f)/chin1/TWO Ryz = Ryz + (fyz - chiy*chiz/chin1/TWO + gyz * f)/chin1/TWO ! covariant second derivatives of the lapse respect to physical metric call fdderivs(ex,Lap,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z, & SYM,SYM,SYM,symmetry,Lev) gxxx = (gupxx * chix + gupxy * chiy + gupxz * chiz)/chin1 gxxy = (gupxy * chix + gupyy * chiy + gupyz * chiz)/chin1 gxxz = (gupxz * chix + gupyz * chiy + gupzz * chiz)/chin1 ! now get physical second kind of connection Gamxxx = Gamxxx - ( (chix + chix)/chin1 - gxx * gxxx )*HALF Gamyxx = Gamyxx - ( - gxx * gxxy )*HALF Gamzxx = Gamzxx - ( - gxx * gxxz )*HALF Gamxyy = Gamxyy - ( - gyy * gxxx )*HALF Gamyyy = Gamyyy - ( (chiy + chiy)/chin1 - gyy * gxxy )*HALF Gamzyy = Gamzyy - ( - gyy * gxxz )*HALF Gamxzz = Gamxzz - ( - gzz * gxxx )*HALF Gamyzz = Gamyzz - ( - gzz * gxxy )*HALF Gamzzz = Gamzzz - ( (chiz + chiz)/chin1 - gzz * gxxz )*HALF Gamxxy = Gamxxy - ( chiy /chin1 - gxy * gxxx )*HALF Gamyxy = Gamyxy - ( chix /chin1 - gxy * gxxy )*HALF Gamzxy = Gamzxy - ( - gxy * gxxz )*HALF Gamxxz = Gamxxz - ( chiz /chin1 - gxz * gxxx )*HALF Gamyxz = Gamyxz - ( - gxz * gxxy )*HALF Gamzxz = Gamzxz - ( chix /chin1 - gxz * gxxz )*HALF Gamxyz = Gamxyz - ( - gyz * gxxx )*HALF Gamyyz = Gamyyz - ( chiz /chin1 - gyz * gxxy )*HALF Gamzyz = Gamzyz - ( chiy /chin1 - gyz * gxxz )*HALF fxx = fxx - Gamxxx*Lapx - Gamyxx*Lapy - Gamzxx*Lapz fyy = fyy - Gamxyy*Lapx - Gamyyy*Lapy - Gamzyy*Lapz fzz = fzz - Gamxzz*Lapx - Gamyzz*Lapy - Gamzzz*Lapz fxy = fxy - Gamxxy*Lapx - Gamyxy*Lapy - Gamzxy*Lapz fxz = fxz - Gamxxz*Lapx - Gamyxz*Lapy - Gamzxz*Lapz fyz = fyz - Gamxyz*Lapx - Gamyyz*Lapy - Gamzyz*Lapz ! store D^i D_i Lap in trK_rhs upto chi trK_rhs = gupxx * fxx + gupyy * fyy + gupzz * fzz + & TWO* ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) #if 1 !! follow bam code S = chin1 * ( gupxx * Sxx + gupyy * Syy + gupzz * Szz + & TWO * ( gupxy * Sxy + gupxz * Sxz + gupyz * Syz ) ) f = F2o3 * trK * trK -(& gupxx * ( & gupxx * Axx * Axx + gupyy * Axy * Axy + gupzz * Axz * Axz + & TWO * (gupxy * Axx * Axy + gupxz * Axx * Axz + gupyz * Axy * Axz) ) + & gupyy * ( & gupxx * Axy * Axy + gupyy * Ayy * Ayy + gupzz * Ayz * Ayz + & TWO * (gupxy * Axy * Ayy + gupxz * Axy * Ayz + gupyz * Ayy * Ayz) ) + & gupzz * ( & gupxx * Axz * Axz + gupyy * Ayz * Ayz + gupzz * Azz * Azz + & TWO * (gupxy * Axz * Ayz + gupxz * Axz * Azz + gupyz * Ayz * Azz) ) + & TWO * ( & gupxy * ( & gupxx * Axx * Axy + gupyy * Axy * Ayy + gupzz * Axz * Ayz + & gupxy * (Axx * Ayy + Axy * Axy) + & gupxz * (Axx * Ayz + Axz * Axy) + & gupyz * (Axy * Ayz + Axz * Ayy) ) + & gupxz * ( & gupxx * Axx * Axz + gupyy * Axy * Ayz + gupzz * Axz * Azz + & gupxy * (Axx * Ayz + Axy * Axz) + & gupxz * (Axx * Azz + Axz * Axz) + & gupyz * (Axy * Azz + Axz * Ayz) ) + & gupyz * ( & gupxx * Axy * Axz + gupyy * Ayy * Ayz + gupzz * Ayz * Azz + & gupxy * (Axy * Ayz + Ayy * Axz) + & gupxz * (Axy * Azz + Ayz * Axz) + & gupyz * (Ayy * Azz + Ayz * Ayz) ) )) -1.6d1*PI*rho + EIGHT * PI * S f = - F1o3 *( gupxx * fxx + gupyy * fyy + gupzz * fzz + & TWO* ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) + alpn1/chin1*f) fxx = alpn1 * (Rxx - EIGHT * PI * Sxx) - fxx fxy = alpn1 * (Rxy - EIGHT * PI * Sxy) - fxy fxz = alpn1 * (Rxz - EIGHT * PI * Sxz) - fxz fyy = alpn1 * (Ryy - EIGHT * PI * Syy) - fyy fyz = alpn1 * (Ryz - EIGHT * PI * Syz) - fyz fzz = alpn1 * (Rzz - EIGHT * PI * Szz) - fzz #else ! Add lapse and S_ij parts to Ricci tensor: fxx = alpn1 * (Rxx - EIGHT * PI * Sxx) - fxx fxy = alpn1 * (Rxy - EIGHT * PI * Sxy) - fxy fxz = alpn1 * (Rxz - EIGHT * PI * Sxz) - fxz fyy = alpn1 * (Ryy - EIGHT * PI * Syy) - fyy fyz = alpn1 * (Ryz - EIGHT * PI * Syz) - fyz fzz = alpn1 * (Rzz - EIGHT * PI * Szz) - fzz ! Compute trace-free part (note: chi^-1 and chi cancel!): f = F1o3 *( gupxx * fxx + gupyy * fyy + gupzz * fzz + & TWO* ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) ) #endif Axx_rhs = fxx - gxx * f Ayy_rhs = fyy - gyy * f Azz_rhs = fzz - gzz * f Axy_rhs = fxy - gxy * f Axz_rhs = fxz - gxz * f Ayz_rhs = fyz - gyz * f ! Now: store A_il A^l_j into fij: fxx = gupxx * Axx * Axx + gupyy * Axy * Axy + gupzz * Axz * Axz + & TWO * (gupxy * Axx * Axy + gupxz * Axx * Axz + gupyz * Axy * Axz) fyy = gupxx * Axy * Axy + gupyy * Ayy * Ayy + gupzz * Ayz * Ayz + & TWO * (gupxy * Axy * Ayy + gupxz * Axy * Ayz + gupyz * Ayy * Ayz) fzz = gupxx * Axz * Axz + gupyy * Ayz * Ayz + gupzz * Azz * Azz + & TWO * (gupxy * Axz * Ayz + gupxz * Axz * Azz + gupyz * Ayz * Azz) fxy = gupxx * Axx * Axy + gupyy * Axy * Ayy + gupzz * Axz * Ayz + & gupxy *(Axx * Ayy + Axy * Axy) + & gupxz *(Axx * Ayz + Axz * Axy) + & gupyz *(Axy * Ayz + Axz * Ayy) fxz = gupxx * Axx * Axz + gupyy * Axy * Ayz + gupzz * Axz * Azz + & gupxy *(Axx * Ayz + Axy * Axz) + & gupxz *(Axx * Azz + Axz * Axz) + & gupyz *(Axy * Azz + Axz * Ayz) fyz = gupxx * Axy * Axz + gupyy * Ayy * Ayz + gupzz * Ayz * Azz + & gupxy *(Axy * Ayz + Ayy * Axz) + & gupxz *(Axy * Azz + Ayz * Axz) + & gupyz *(Ayy * Azz + Ayz * Ayz) f = chin1 ! store D^i D_i Lap in trK_rhs trK_rhs = f*trK_rhs Axx_rhs = f * Axx_rhs+ alpn1 * (trK * Axx - TWO * fxx) + & TWO * ( Axx * betaxx + Axy * betayx + Axz * betazx )- & F2o3 * Axx * div_beta Ayy_rhs = f * Ayy_rhs+ alpn1 * (trK * Ayy - TWO * fyy) + & TWO * ( Axy * betaxy + Ayy * betayy + Ayz * betazy )- & F2o3 * Ayy * div_beta Azz_rhs = f * Azz_rhs+ alpn1 * (trK * Azz - TWO * fzz) + & TWO * ( Axz * betaxz + Ayz * betayz + Azz * betazz )- & F2o3 * Azz * div_beta Axy_rhs = f * Axy_rhs+ alpn1 *( trK * Axy - TWO * fxy )+ & Axx * betaxy + Axz * betazy + & Ayy * betayx + Ayz * betazx + & F1o3 * Axy * div_beta - Axy * betazz Ayz_rhs = f * Ayz_rhs+ alpn1 *( trK * Ayz - TWO * fyz )+ & Axy * betaxz + Ayy * betayz + & Axz * betaxy + Azz * betazy + & F1o3 * Ayz * div_beta - Ayz * betaxx Axz_rhs = f * Axz_rhs+ alpn1 *( trK * Axz - TWO * fxz )+ & Axx * betaxz + Axy * betayz + & Ayz * betayx + Azz * betazx + & F1o3 * Axz * div_beta - Axz * betayy !rhs for Aij ! Compute trace of S_ij S = f * ( gupxx * Sxx + gupyy * Syy + gupzz * Szz + & TWO * ( gupxy * Sxy + gupxz * Sxz + gupyz * Syz ) ) trK_rhs = - trK_rhs + alpn1 *( F1o3 * trK * trK + & gupxx * fxx + gupyy * fyy + gupzz * fzz + & TWO * ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) + & FOUR * PI * ( rho + S )) !rhs for trK !!!! gauge variable part Lap_rhs = -TWO*alpn1*trK #if (GAUGE == 0) betax_rhs = FF*dtSfx betay_rhs = FF*dtSfy betaz_rhs = FF*dtSfz dtSfx_rhs = Gamx_rhs - eta*dtSfx dtSfy_rhs = Gamy_rhs - eta*dtSfy dtSfz_rhs = Gamz_rhs - eta*dtSfz #elif (GAUGE == 1) betax_rhs = Gamx - eta*betax betay_rhs = Gamy - eta*betay betaz_rhs = Gamz - eta*betaz dtSfx_rhs = ZEO dtSfy_rhs = ZEO dtSfz_rhs = ZEO #elif (GAUGE == 2) betax_rhs = FF*dtSfx betay_rhs = FF*dtSfy betaz_rhs = FF*dtSfz call fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) reta = gupxx * dtSfx_rhs * dtSfx_rhs + gupyy * dtSfy_rhs * dtSfy_rhs + gupzz * dtSfz_rhs * dtSfz_rhs + & TWO * (gupxy * dtSfx_rhs * dtSfy_rhs + gupxz * dtSfx_rhs * dtSfz_rhs + gupyz * dtSfy_rhs * dtSfz_rhs) reta = 1.31d0/2*dsqrt(reta/chin1)/(1-dsqrt(chin1))**2 dtSfx_rhs = Gamx_rhs - reta*dtSfx dtSfy_rhs = Gamy_rhs - reta*dtSfy dtSfz_rhs = Gamz_rhs - reta*dtSfz #elif (GAUGE == 3) betax_rhs = FF*dtSfx betay_rhs = FF*dtSfy betaz_rhs = FF*dtSfz call fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) reta = gupxx * dtSfx_rhs * dtSfx_rhs + gupyy * dtSfy_rhs * dtSfy_rhs + gupzz * dtSfz_rhs * dtSfz_rhs + & TWO * (gupxy * dtSfx_rhs * dtSfy_rhs + gupxz * dtSfx_rhs * dtSfz_rhs + gupyz * dtSfy_rhs * dtSfz_rhs) reta = 1.31d0/2*dsqrt(reta/chin1)/(1-chin1)**2 dtSfx_rhs = Gamx_rhs - reta*dtSfx dtSfy_rhs = Gamy_rhs - reta*dtSfy dtSfz_rhs = Gamz_rhs - reta*dtSfz #elif (GAUGE == 4) call fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) reta = gupxx * dtSfx_rhs * dtSfx_rhs + gupyy * dtSfy_rhs * dtSfy_rhs + gupzz * dtSfz_rhs * dtSfz_rhs + & TWO * (gupxy * dtSfx_rhs * dtSfy_rhs + gupxz * dtSfx_rhs * dtSfz_rhs + gupyz * dtSfy_rhs * dtSfz_rhs) reta = 1.31d0/2*dsqrt(reta/chin1)/(1-dsqrt(chin1))**2 betax_rhs = FF*Gamx - reta*betax betay_rhs = FF*Gamy - reta*betay betaz_rhs = FF*Gamz - reta*betaz dtSfx_rhs = ZEO dtSfy_rhs = ZEO dtSfz_rhs = ZEO #elif (GAUGE == 5) call fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) reta = gupxx * dtSfx_rhs * dtSfx_rhs + gupyy * dtSfy_rhs * dtSfy_rhs + gupzz * dtSfz_rhs * dtSfz_rhs + & TWO * (gupxy * dtSfx_rhs * dtSfy_rhs + gupxz * dtSfx_rhs * dtSfz_rhs + gupyz * dtSfy_rhs * dtSfz_rhs) reta = 1.31d0/2*dsqrt(reta/chin1)/(1-chin1)**2 betax_rhs = FF*Gamx - reta*betax betay_rhs = FF*Gamy - reta*betay betaz_rhs = FF*Gamz - reta*betaz dtSfx_rhs = ZEO dtSfy_rhs = ZEO dtSfz_rhs = ZEO #elif (GAUGE == 6) if(BHN==2)then M = Mass(1)+Mass(2) A = 2.d0/M w1 = 1.2d1 w2 = w1 C1 = 1.d0/Mass(1) - A C2 = 1.d0/Mass(2) - A do k=1,ex(3) do j=1,ex(2) do i=1,ex(1) r1 = ((Porg(1)-X(i))**2+(Porg(2)-Y(j))**2+(Porg(3)-Z(k))**2)/ & ((Porg(1)-Porg(4))**2+(Porg(2)-Porg(5))**2+(Porg(3)-Porg(6))**2) r2 = ((Porg(4)-X(i))**2+(Porg(5)-Y(j))**2+(Porg(6)-Z(k))**2)/ & ((Porg(1)-Porg(4))**2+(Porg(2)-Porg(5))**2+(Porg(3)-Porg(6))**2) reta(i,j,k) = A + C1/(ONE+w1*r1) + C2/(ONE+w2*r2) enddo enddo enddo else write(*,*) "not support BH_num in Jason's form 1",BHN endif betax_rhs = FF*dtSfx betay_rhs = FF*dtSfy betaz_rhs = FF*dtSfz dtSfx_rhs = Gamx_rhs - reta*dtSfx dtSfy_rhs = Gamy_rhs - reta*dtSfy dtSfz_rhs = Gamz_rhs - reta*dtSfz #elif (GAUGE == 7) if(BHN==2)then M = Mass(1)+Mass(2) A = 2.d0/M w1 = 1.2d1 w2 = w1 C1 = 1.d0/Mass(1) - A C2 = 1.d0/Mass(2) - A do k=1,ex(3) do j=1,ex(2) do i=1,ex(1) r1 = ((Porg(1)-X(i))**2+(Porg(2)-Y(j))**2+(Porg(3)-Z(k))**2)/ & ((Porg(1)-Porg(4))**2+(Porg(2)-Porg(5))**2+(Porg(3)-Porg(6))**2) r2 = ((Porg(4)-X(i))**2+(Porg(5)-Y(j))**2+(Porg(6)-Z(k))**2)/ & ((Porg(1)-Porg(4))**2+(Porg(2)-Porg(5))**2+(Porg(3)-Porg(6))**2) reta(i,j,k) = A + C1*dexp(-w1*r1) + C2*dexp(-w2*r2) enddo enddo enddo else write(*,*) "not support BH_num in Jason's form 2",BHN endif betax_rhs = FF*dtSfx betay_rhs = FF*dtSfy betaz_rhs = FF*dtSfz dtSfx_rhs = Gamx_rhs - reta*dtSfx dtSfy_rhs = Gamy_rhs - reta*dtSfy dtSfz_rhs = Gamz_rhs - reta*dtSfz #endif SSS(1)=SYM SSS(2)=SYM SSS(3)=SYM AAS(1)=ANTI AAS(2)=ANTI AAS(3)=SYM ASA(1)=ANTI ASA(2)=SYM ASA(3)=ANTI SAA(1)=SYM SAA(2)=ANTI SAA(3)=ANTI ASS(1)=ANTI ASS(2)=SYM ASS(3)=SYM SAS(1)=SYM SAS(2)=ANTI SAS(3)=SYM SSA(1)=SYM SSA(2)=SYM SSA(3)=ANTI !!!!!!!!!advection term + Kreiss-Oliger dissipation (merged for cache efficiency) ! lopsided_kodis shares the symmetry_bd buffer between advection and ! dissipation, eliminating redundant full-grid copies. For metric variables ! gxx/gyy/gzz (=dxx/dyy/dzz+1): stencil coefficients sum to zero, ! so the constant offset has no effect on dissipation. call lopsided_kodis(ex,X,Y,Z,dxx,gxx_rhs,betax,betay,betaz,Symmetry,SSS,eps) call lopsided_kodis(ex,X,Y,Z,gxy,gxy_rhs,betax,betay,betaz,Symmetry,AAS,eps) call lopsided_kodis(ex,X,Y,Z,gxz,gxz_rhs,betax,betay,betaz,Symmetry,ASA,eps) call lopsided_kodis(ex,X,Y,Z,dyy,gyy_rhs,betax,betay,betaz,Symmetry,SSS,eps) call lopsided_kodis(ex,X,Y,Z,gyz,gyz_rhs,betax,betay,betaz,Symmetry,SAA,eps) call lopsided_kodis(ex,X,Y,Z,dzz,gzz_rhs,betax,betay,betaz,Symmetry,SSS,eps) call lopsided_kodis(ex,X,Y,Z,Axx,Axx_rhs,betax,betay,betaz,Symmetry,SSS,eps) call lopsided_kodis(ex,X,Y,Z,Axy,Axy_rhs,betax,betay,betaz,Symmetry,AAS,eps) call lopsided_kodis(ex,X,Y,Z,Axz,Axz_rhs,betax,betay,betaz,Symmetry,ASA,eps) call lopsided_kodis(ex,X,Y,Z,Ayy,Ayy_rhs,betax,betay,betaz,Symmetry,SSS,eps) call lopsided_kodis(ex,X,Y,Z,Ayz,Ayz_rhs,betax,betay,betaz,Symmetry,SAA,eps) call lopsided_kodis(ex,X,Y,Z,Azz,Azz_rhs,betax,betay,betaz,Symmetry,SSS,eps) call lopsided_kodis(ex,X,Y,Z,chi,chi_rhs,betax,betay,betaz,Symmetry,SSS,eps) call lopsided_kodis(ex,X,Y,Z,trK,trK_rhs,betax,betay,betaz,Symmetry,SSS,eps) call lopsided_kodis(ex,X,Y,Z,Gamx,Gamx_rhs,betax,betay,betaz,Symmetry,ASS,eps) call lopsided_kodis(ex,X,Y,Z,Gamy,Gamy_rhs,betax,betay,betaz,Symmetry,SAS,eps) call lopsided_kodis(ex,X,Y,Z,Gamz,Gamz_rhs,betax,betay,betaz,Symmetry,SSA,eps) #if 1 !! bam does not apply dissipation on gauge variables call lopsided_kodis(ex,X,Y,Z,Lap,Lap_rhs,betax,betay,betaz,Symmetry,SSS,eps) #if (GAUGE == 0 || GAUGE == 1 || GAUGE == 2 || GAUGE == 3 || GAUGE == 4 || GAUGE == 5 || GAUGE == 6 || GAUGE == 7) call lopsided_kodis(ex,X,Y,Z,betax,betax_rhs,betax,betay,betaz,Symmetry,ASS,eps) call lopsided_kodis(ex,X,Y,Z,betay,betay_rhs,betax,betay,betaz,Symmetry,SAS,eps) call lopsided_kodis(ex,X,Y,Z,betaz,betaz_rhs,betax,betay,betaz,Symmetry,SSA,eps) #endif #if (GAUGE == 0 || GAUGE == 2 || GAUGE == 3 || GAUGE == 6 || GAUGE == 7) call lopsided_kodis(ex,X,Y,Z,dtSfx,dtSfx_rhs,betax,betay,betaz,Symmetry,ASS,eps) call lopsided_kodis(ex,X,Y,Z,dtSfy,dtSfy_rhs,betax,betay,betaz,Symmetry,SAS,eps) call lopsided_kodis(ex,X,Y,Z,dtSfz,dtSfz_rhs,betax,betay,betaz,Symmetry,SSA,eps) #endif #else ! No dissipation on gauge variables (advection only) call lopsided(ex,X,Y,Z,Lap,Lap_rhs,betax,betay,betaz,Symmetry,SSS) #if (GAUGE == 0 || GAUGE == 1 || GAUGE == 2 || GAUGE == 3 || GAUGE == 4 || GAUGE == 5 || GAUGE == 6 || GAUGE == 7) call lopsided(ex,X,Y,Z,betax,betax_rhs,betax,betay,betaz,Symmetry,ASS) call lopsided(ex,X,Y,Z,betay,betay_rhs,betax,betay,betaz,Symmetry,SAS) call lopsided(ex,X,Y,Z,betaz,betaz_rhs,betax,betay,betaz,Symmetry,SSA) #endif #if (GAUGE == 0 || GAUGE == 2 || GAUGE == 3 || GAUGE == 6 || GAUGE == 7) call lopsided(ex,X,Y,Z,dtSfx,dtSfx_rhs,betax,betay,betaz,Symmetry,ASS) call lopsided(ex,X,Y,Z,dtSfy,dtSfy_rhs,betax,betay,betaz,Symmetry,SAS) call lopsided(ex,X,Y,Z,dtSfz,dtSfz_rhs,betax,betay,betaz,Symmetry,SSA) #endif #endif if(co == 0)then ! ham_Res = trR + 2/3 * K^2 - A_ij * A^ij - 16 * PI * rho ! here trR is respect to physical metric ham_Res = gupxx * Rxx + gupyy * Ryy + gupzz * Rzz + & TWO* ( gupxy * Rxy + gupxz * Rxz + gupyz * Ryz ) ham_Res = chin1*ham_Res + F2o3 * trK * trK -(& gupxx * ( & gupxx * Axx * Axx + gupyy * Axy * Axy + gupzz * Axz * Axz + & TWO * (gupxy * Axx * Axy + gupxz * Axx * Axz + gupyz * Axy * Axz) ) + & gupyy * ( & gupxx * Axy * Axy + gupyy * Ayy * Ayy + gupzz * Ayz * Ayz + & TWO * (gupxy * Axy * Ayy + gupxz * Axy * Ayz + gupyz * Ayy * Ayz) ) + & gupzz * ( & gupxx * Axz * Axz + gupyy * Ayz * Ayz + gupzz * Azz * Azz + & TWO * (gupxy * Axz * Ayz + gupxz * Axz * Azz + gupyz * Ayz * Azz) ) + & TWO * ( & gupxy * ( & gupxx * Axx * Axy + gupyy * Axy * Ayy + gupzz * Axz * Ayz + & gupxy * (Axx * Ayy + Axy * Axy) + & gupxz * (Axx * Ayz + Axz * Axy) + & gupyz * (Axy * Ayz + Axz * Ayy) ) + & gupxz * ( & gupxx * Axx * Axz + gupyy * Axy * Ayz + gupzz * Axz * Azz + & gupxy * (Axx * Ayz + Axy * Axz) + & gupxz * (Axx * Azz + Axz * Axz) + & gupyz * (Axy * Azz + Axz * Ayz) ) + & gupyz * ( & gupxx * Axy * Axz + gupyy * Ayy * Ayz + gupzz * Ayz * Azz + & gupxy * (Axy * Ayz + Ayy * Axz) + & gupxz * (Axy * Azz + Ayz * Axz) + & gupyz * (Ayy * Azz + Ayz * Ayz) ) ))- F16 * PI * rho ! mov_Res_j = gupkj*(-1/chi d_k chi*A_ij + D_k A_ij) - 2/3 d_j trK - 8 PI s_j where D respect to physical metric ! store D_i A_jk - 1/chi d_i chi*A_jk in gjk_i call fderivs(ex,Axx,gxxx,gxxy,gxxz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0) call fderivs(ex,Axy,gxyx,gxyy,gxyz,X,Y,Z,ANTI,ANTI,SYM ,Symmetry,0) call fderivs(ex,Axz,gxzx,gxzy,gxzz,X,Y,Z,ANTI,SYM ,ANTI,Symmetry,0) call fderivs(ex,Ayy,gyyx,gyyy,gyyz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0) call fderivs(ex,Ayz,gyzx,gyzy,gyzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,0) call fderivs(ex,Azz,gzzx,gzzy,gzzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0) gxxx = gxxx - ( Gamxxx * Axx + Gamyxx * Axy + Gamzxx * Axz & + Gamxxx * Axx + Gamyxx * Axy + Gamzxx * Axz) - chix*Axx/chin1 gxyx = gxyx - ( Gamxxy * Axx + Gamyxy * Axy + Gamzxy * Axz & + Gamxxx * Axy + Gamyxx * Ayy + Gamzxx * Ayz) - chix*Axy/chin1 gxzx = gxzx - ( Gamxxz * Axx + Gamyxz * Axy + Gamzxz * Axz & + Gamxxx * Axz + Gamyxx * Ayz + Gamzxx * Azz) - chix*Axz/chin1 gyyx = gyyx - ( Gamxxy * Axy + Gamyxy * Ayy + Gamzxy * Ayz & + Gamxxy * Axy + Gamyxy * Ayy + Gamzxy * Ayz) - chix*Ayy/chin1 gyzx = gyzx - ( Gamxxz * Axy + Gamyxz * Ayy + Gamzxz * Ayz & + Gamxxy * Axz + Gamyxy * Ayz + Gamzxy * Azz) - chix*Ayz/chin1 gzzx = gzzx - ( Gamxxz * Axz + Gamyxz * Ayz + Gamzxz * Azz & + Gamxxz * Axz + Gamyxz * Ayz + Gamzxz * Azz) - chix*Azz/chin1 gxxy = gxxy - ( Gamxxy * Axx + Gamyxy * Axy + Gamzxy * Axz & + Gamxxy * Axx + Gamyxy * Axy + Gamzxy * Axz) - chiy*Axx/chin1 gxyy = gxyy - ( Gamxyy * Axx + Gamyyy * Axy + Gamzyy * Axz & + Gamxxy * Axy + Gamyxy * Ayy + Gamzxy * Ayz) - chiy*Axy/chin1 gxzy = gxzy - ( Gamxyz * Axx + Gamyyz * Axy + Gamzyz * Axz & + Gamxxy * Axz + Gamyxy * Ayz + Gamzxy * Azz) - chiy*Axz/chin1 gyyy = gyyy - ( Gamxyy * Axy + Gamyyy * Ayy + Gamzyy * Ayz & + Gamxyy * Axy + Gamyyy * Ayy + Gamzyy * Ayz) - chiy*Ayy/chin1 gyzy = gyzy - ( Gamxyz * Axy + Gamyyz * Ayy + Gamzyz * Ayz & + Gamxyy * Axz + Gamyyy * Ayz + Gamzyy * Azz) - chiy*Ayz/chin1 gzzy = gzzy - ( Gamxyz * Axz + Gamyyz * Ayz + Gamzyz * Azz & + Gamxyz * Axz + Gamyyz * Ayz + Gamzyz * Azz) - chiy*Azz/chin1 gxxz = gxxz - ( Gamxxz * Axx + Gamyxz * Axy + Gamzxz * Axz & + Gamxxz * Axx + Gamyxz * Axy + Gamzxz * Axz) - chiz*Axx/chin1 gxyz = gxyz - ( Gamxyz * Axx + Gamyyz * Axy + Gamzyz * Axz & + Gamxxz * Axy + Gamyxz * Ayy + Gamzxz * Ayz) - chiz*Axy/chin1 gxzz = gxzz - ( Gamxzz * Axx + Gamyzz * Axy + Gamzzz * Axz & + Gamxxz * Axz + Gamyxz * Ayz + Gamzxz * Azz) - chiz*Axz/chin1 gyyz = gyyz - ( Gamxyz * Axy + Gamyyz * Ayy + Gamzyz * Ayz & + Gamxyz * Axy + Gamyyz * Ayy + Gamzyz * Ayz) - chiz*Ayy/chin1 gyzz = gyzz - ( Gamxzz * Axy + Gamyzz * Ayy + Gamzzz * Ayz & + Gamxyz * Axz + Gamyyz * Ayz + Gamzyz * Azz) - chiz*Ayz/chin1 gzzz = gzzz - ( Gamxzz * Axz + Gamyzz * Ayz + Gamzzz * Azz & + Gamxzz * Axz + Gamyzz * Ayz + Gamzzz * Azz) - chiz*Azz/chin1 movx_Res = gupxx*gxxx + gupyy*gxyy + gupzz*gxzz & +gupxy*gxyx + gupxz*gxzx + gupyz*gxzy & +gupxy*gxxy + gupxz*gxxz + gupyz*gxyz movy_Res = gupxx*gxyx + gupyy*gyyy + gupzz*gyzz & +gupxy*gyyx + gupxz*gyzx + gupyz*gyzy & +gupxy*gxyy + gupxz*gxyz + gupyz*gyyz movz_Res = gupxx*gxzx + gupyy*gyzy + gupzz*gzzz & +gupxy*gyzx + gupxz*gzzx + gupyz*gzzy & +gupxy*gxzy + gupxz*gxzz + gupyz*gyzz movx_Res = movx_Res - F2o3*Kx - F8*PI*sx movy_Res = movy_Res - F2o3*Ky - F8*PI*sy movz_Res = movz_Res - F2o3*Kz - F8*PI*sz endif #if (ABV == 1) call ricci_gamma(ex, X, Y, Z, & chi, & dxx , gxy , gxz , dyy , gyz , dzz,& Gamx , Gamy , Gamz , & 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) call constraint_bssn(ex, X, Y, Z,& chi,trK, & dxx,gxy,gxz,dyy,gyz,dzz, & Axx,Axy,Axz,Ayy,Ayz,Azz, & Gamx,Gamy,Gamz,& Lap,betax,betay,betaz,rho,Sx,Sy,Sz,& Gamxxx, Gamxxy, Gamxxz,Gamxyy, Gamxyz, Gamxzz, & Gamyxx, Gamyxy, Gamyxz,Gamyyy, Gamyyz, Gamyzz, & Gamzxx, Gamzxy, Gamzxz,Gamzyy, Gamzyz, Gamzzz, & Rxx,Rxy,Rxz,Ryy,Ryz,Rzz, & ham_Res,movx_Res,movy_Res,movz_Res,Gmx_Res,Gmy_Res,Gmz_Res, & Symmetry) #endif #if 0 #define i 2 if(Lev == 1)then write(*,*) X(i),Y(i),Z(i) write(*,*) Axx(i,i,i),Axy(i,i,i),Axz(i,i,i),Ayy(i,i,i),Ayz(i,i,i),Azz(i,i,i) write(*,*) 1+Lap(i,i,i),dtSfx(i,i,i),dtSfy(i,i,i),dtSfz(i,i,i) write(*,*) betax(i,i,i),betay(i,i,i),betaz(i,i,i) write(*,*) 1+chi(i,i,i),Gamx(i,i,i),Gamy(i,i,i),Gamz(i,i,i) write(*,*) gxx(i,i,i),gxy(i,i,i),gxz(i,i,i),gyy(i,i,i),gyz(i,i,i),gzz(i,i,i) write(*,*) trK(i,i,i) write(*,*) "=====" write(*,*) Axx_rhs(i,i,i),Axy_rhs(i,i,i),Axz_rhs(i,i,i),Ayy_rhs(i,i,i),Ayz_rhs(i,i,i),Azz_rhs(i,i,i) write(*,*) Lap_rhs(i,i,i),dtSfx_rhs(i,i,i),dtSfy_rhs(i,i,i),dtSfz_rhs(i,i,i) write(*,*) betax_rhs(i,i,i),betay_rhs(i,i,i),betaz_rhs(i,i,i) write(*,*) chi_rhs(i,i,i),Gamx_rhs(i,i,i),Gamy_rhs(i,i,i),Gamz_rhs(i,i,i) write(*,*) gxx_rhs(i,i,i),gxy_rhs(i,i,i),gxz_rhs(i,i,i),gyy_rhs(i,i,i),gyz_rhs(i,i,i),gzz_rhs(i,i,i) write(*,*) trK_rhs(i,i,i) endif #undef i !!stop #endif gont = 0 return end function compute_rhs_bssn