From 7bb9042b1859fdd4c828681cd9226921fa4b6a78 Mon Sep 17 00:00:00 2001 From: CGH0S7 <776459475@qq.com> Date: Tue, 3 Mar 2026 15:41:26 +0800 Subject: [PATCH] bssn_rhs(fortran): migrate C kernel loop-fusion optimizations --- AMSS_NCKU_source/bssn_rhs.f90 | 484 +++++++++++++++++++--------------- 1 file changed, 268 insertions(+), 216 deletions(-) diff --git a/AMSS_NCKU_source/bssn_rhs.f90 b/AMSS_NCKU_source/bssn_rhs.f90 index c8c044e..a18be11 100644 --- a/AMSS_NCKU_source/bssn_rhs.f90 +++ b/AMSS_NCKU_source/bssn_rhs.f90 @@ -59,9 +59,10 @@ 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 + 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: @@ -83,11 +84,16 @@ 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, 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 + 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 @@ -96,11 +102,11 @@ real*8, dimension(ex(1),ex(2),ex(3)) :: reta #endif -#if (GAUGE == 6 || GAUGE == 7) - integer :: BHN,i,j,k - real*8, dimension(9) :: Porg - real*8, dimension(3) :: Mass - real*8 :: r1,r2,M,A,w1,w2,C1,C2 +#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) @@ -145,174 +151,204 @@ dY = Y(2) - Y(1) dZ = Z(2) - Z(1) - alpn1 = Lap + ONE - chin1 = chi + ONE - gxx = dxx + ONE - gyy = dyy + ONE - gzz = dzz + ONE + 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) - div_beta = betaxx + betayy + betazz - - call fderivs(ex,chi,chix,chiy,chiz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev) + call fderivs(ex,chi,chix,chiy,chiz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev) - chi_rhs = F2o3 *chin1*( alpn1 * trK - div_beta ) !rhs for chi - - 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) - - gxx_rhs = - TWO * alpn1 * Axx - F2o3 * gxx * div_beta + & - TWO *( gxx * betaxx + gxy * betayx + gxz * betazx) - - gyy_rhs = - TWO * alpn1 * Ayy - F2o3 * gyy * div_beta + & - TWO *( gxy * betaxy + gyy * betayy + gyz * betazy) - - gzz_rhs = - TWO * alpn1 * Azz - F2o3 * gzz * div_beta + & - TWO *( gxz * betaxz + gyz * betayz + gzz * betazz) - - gxy_rhs = - TWO * alpn1 * Axy + F1o3 * gxy * div_beta + & - gxx * betaxy + gxz * betazy + & - gyy * betayx + gyz * betazx & - - gxy * betazz - - gyz_rhs = - TWO * alpn1 * Ayz + F1o3 * gyz * div_beta + & - gxy * betaxz + gyy * betayz + & - gxz * betaxy + gzz * betazy & - - gyz * betaxx - - gxz_rhs = - TWO * alpn1 * Axz + F1o3 * gxz * div_beta + & - gxx * betaxz + gxy * betayz + & - gyz * betayx + gzz * betazx & - - gxz * betayy !rhs for gij - -! 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 - - if(co == 0)then -! Gam^i_Res = Gam^i + gup^ij_,j - Gmx_Res = Gamx - (gupxx*(gupxx*gxxx+gupxy*gxyx+gupxz*gxzx)& - +gupxy*(gupxx*gxyx+gupxy*gyyx+gupxz*gyzx)& - +gupxz*(gupxx*gxzx+gupxy*gyzx+gupxz*gzzx)& - +gupxx*(gupxy*gxxy+gupyy*gxyy+gupyz*gxzy)& - +gupxy*(gupxy*gxyy+gupyy*gyyy+gupyz*gyzy)& - +gupxz*(gupxy*gxzy+gupyy*gyzy+gupyz*gzzy)& - +gupxx*(gupxz*gxxz+gupyz*gxyz+gupzz*gxzz)& - +gupxy*(gupxz*gxyz+gupyz*gyyz+gupzz*gyzz)& - +gupxz*(gupxz*gxzz+gupyz*gyzz+gupzz*gzzz)) - Gmy_Res = Gamy - (gupxx*(gupxy*gxxx+gupyy*gxyx+gupyz*gxzx)& - +gupxy*(gupxy*gxyx+gupyy*gyyx+gupyz*gyzx)& - +gupxz*(gupxy*gxzx+gupyy*gyzx+gupyz*gzzx)& - +gupxy*(gupxy*gxxy+gupyy*gxyy+gupyz*gxzy)& - +gupyy*(gupxy*gxyy+gupyy*gyyy+gupyz*gyzy)& - +gupyz*(gupxy*gxzy+gupyy*gyzy+gupyz*gzzy)& - +gupxy*(gupxz*gxxz+gupyz*gxyz+gupzz*gxzz)& - +gupyy*(gupxz*gxyz+gupyz*gyyz+gupzz*gyzz)& - +gupyz*(gupxz*gxzz+gupyz*gyzz+gupzz*gzzz)) - Gmz_Res = Gamz - (gupxx*(gupxz*gxxx+gupyz*gxyx+gupzz*gxzx)& - +gupxy*(gupxz*gxyx+gupyz*gyyx+gupzz*gyzx)& - +gupxz*(gupxz*gxzx+gupyz*gyzx+gupzz*gzzx)& - +gupxy*(gupxz*gxxy+gupyz*gxyy+gupzz*gxzy)& - +gupyy*(gupxz*gxyy+gupyz*gyyy+gupzz*gyzy)& - +gupyz*(gupxz*gxzy+gupyz*gyzy+gupzz*gzzy)& - +gupxz*(gupxz*gxxz+gupyz*gxyz+gupzz*gxzz)& - +gupyz*(gupxz*gxyz+gupyz*gyyz+gupzz*gyzz)& - +gupzz*(gupxz*gxzz+gupyz*gyzz+gupzz*gzzz)) - endif - -! second kind of connection - Gamxxx =HALF*( gupxx*gxxx + gupxy*(TWO*gxyx - gxxy ) + gupxz*(TWO*gxzx - gxxz )) - Gamyxx =HALF*( gupxy*gxxx + gupyy*(TWO*gxyx - gxxy ) + gupyz*(TWO*gxzx - gxxz )) - Gamzxx =HALF*( gupxz*gxxx + gupyz*(TWO*gxyx - gxxy ) + gupzz*(TWO*gxzx - gxxz )) - - Gamxyy =HALF*( gupxx*(TWO*gxyy - gyyx ) + gupxy*gyyy + gupxz*(TWO*gyzy - gyyz )) - Gamyyy =HALF*( gupxy*(TWO*gxyy - gyyx ) + gupyy*gyyy + gupyz*(TWO*gyzy - gyyz )) - Gamzyy =HALF*( gupxz*(TWO*gxyy - gyyx ) + gupyz*gyyy + gupzz*(TWO*gyzy - gyyz )) - - Gamxzz =HALF*( gupxx*(TWO*gxzz - gzzx ) + gupxy*(TWO*gyzz - gzzy ) + gupxz*gzzz) - Gamyzz =HALF*( gupxy*(TWO*gxzz - gzzx ) + gupyy*(TWO*gyzz - gzzy ) + gupyz*gzzz) - Gamzzz =HALF*( gupxz*(TWO*gxzz - gzzx ) + gupyz*(TWO*gyzz - gzzy ) + gupzz*gzzz) - - Gamxxy =HALF*( gupxx*gxxy + gupxy*gyyx + gupxz*( gxzy + gyzx - gxyz ) ) - Gamyxy =HALF*( gupxy*gxxy + gupyy*gyyx + gupyz*( gxzy + gyzx - gxyz ) ) - Gamzxy =HALF*( gupxz*gxxy + gupyz*gyyx + gupzz*( gxzy + gyzx - gxyz ) ) - - Gamxxz =HALF*( gupxx*gxxz + gupxy*( gxyz + gyzx - gxzy ) + gupxz*gzzx ) - Gamyxz =HALF*( gupxy*gxxz + gupyy*( gxyz + gyzx - gxzy ) + gupyz*gzzx ) - Gamzxz =HALF*( gupxz*gxxz + gupyz*( gxyz + gyzx - gxzy ) + gupzz*gzzx ) - - Gamxyz =HALF*( gupxx*( gxyz + gxzy - gyzx ) + gupxy*gyyz + gupxz*gzzy ) - Gamyyz =HALF*( gupxy*( gxyz + gxzy - gyzx ) + gupyy*gyyz + gupyz*gzzy ) - Gamzyz =HALF*( gupxz*( gxyz + gxzy - gyzx ) + gupyz*gyyz + gupzz*gzzy ) -! Raise indices of \tilde A_{ij} and store in R_ij - - Rxx = gupxx * gupxx * Axx + gupxy * gupxy * Ayy + gupxz * gupxz * Azz + & - TWO*(gupxx * gupxy * Axy + gupxx * gupxz * Axz + gupxy * gupxz * Ayz) - - Ryy = gupxy * gupxy * Axx + gupyy * gupyy * Ayy + gupyz * gupyz * Azz + & - TWO*(gupxy * gupyy * Axy + gupxy * gupyz * Axz + gupyy * gupyz * Ayz) - - Rzz = gupxz * gupxz * Axx + gupyz * gupyz * Ayy + gupzz * gupzz * Azz + & - TWO*(gupxz * gupyz * Axy + gupxz * gupzz * Axz + gupyz * gupzz * Ayz) - - Rxy = gupxx * gupxy * Axx + gupxy * gupyy * Ayy + gupxz * gupyz * Azz + & - (gupxx * gupyy + gupxy * gupxy)* Axy + & - (gupxx * gupyz + gupxz * gupxy)* Axz + & - (gupxy * gupyz + gupxz * gupyy)* Ayz - - Rxz = gupxx * gupxz * Axx + gupxy * gupyz * Ayy + gupxz * gupzz * Azz + & - (gupxx * gupyz + gupxy * gupxz)* Axy + & - (gupxx * gupzz + gupxz * gupxz)* Axz + & - (gupxy * gupzz + gupxz * gupyz)* Ayz - - Ryz = gupxy * gupxz * Axx + gupyy * gupyz * Ayy + gupyz * gupzz * Azz + & - (gupxy * gupyz + gupyy * gupxz)* Axy + & - (gupxy * gupzz + gupyz * gupxz)* Axz + & - (gupyy * gupzz + gupyz * gupyz)* Ayz - -! 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) - - Gamx_rhs = - TWO * ( Lapx * Rxx + Lapy * Rxy + Lapz * Rxz ) + & - TWO * alpn1 * ( & - -F3o2/chin1 * ( chix * Rxx + chiy * Rxy + chiz * Rxz ) - & - gupxx * ( F2o3 * Kx + EIGHT * PI * Sx ) - & - gupxy * ( F2o3 * Ky + EIGHT * PI * Sy ) - & - gupxz * ( F2o3 * Kz + EIGHT * PI * Sz ) + & - Gamxxx * Rxx + Gamxyy * Ryy + Gamxzz * Rzz + & - TWO * ( Gamxxy * Rxy + Gamxxz * Rxz + Gamxyz * Ryz ) ) - - Gamy_rhs = - TWO * ( Lapx * Rxy + Lapy * Ryy + Lapz * Ryz ) + & - TWO * alpn1 * ( & - -F3o2/chin1 * ( chix * Rxy + chiy * Ryy + chiz * Ryz ) - & - gupxy * ( F2o3 * Kx + EIGHT * PI * Sx ) - & - gupyy * ( F2o3 * Ky + EIGHT * PI * Sy ) - & - gupyz * ( F2o3 * Kz + EIGHT * PI * Sz ) + & - Gamyxx * Rxx + Gamyyy * Ryy + Gamyzz * Rzz + & - TWO * ( Gamyxy * Rxy + Gamyxz * Rxz + Gamyyz * Ryz ) ) - - Gamz_rhs = - TWO * ( Lapx * Rxz + Lapy * Ryz + Lapz * Rzz ) + & - TWO * alpn1 * ( & - -F3o2/chin1 * ( chix * Rxz + chiy * Ryz + chiz * Rzz ) - & - gupxz * ( F2o3 * Kx + EIGHT * PI * Sx ) - & - gupyz * ( F2o3 * Ky + EIGHT * PI * Sy ) - & - gupzz * ( F2o3 * Kz + EIGHT * PI * Sz ) + & - Gamzxx * Rxx + Gamzyy * Ryy + Gamzzz * Rzz + & - TWO * ( Gamzxy * Rxy + Gamzxz * Rxz + Gamzyz * Ryz ) ) + 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) @@ -321,38 +357,54 @@ call fdderivs(ex,betaz,gxxz,gxyz,gxzz,gyyz,gyzz,gzzz,& X,Y,Z,SYM ,SYM, ANTI,Symmetry,Lev) - fxx = gxxx + gxyy + gxzz - fxy = gxyx + gyyy + gyzz - fxz = gxzx + gyzy + gzzz - - Gamxa = gupxx * Gamxxx + gupyy * Gamxyy + gupzz * Gamxzz + & - TWO*( gupxy * Gamxxy + gupxz * Gamxxz + gupyz * Gamxyz ) - Gamya = gupxx * Gamyxx + gupyy * Gamyyy + gupzz * Gamyzz + & - TWO*( gupxy * Gamyxy + gupxz * Gamyxz + gupyz * Gamyyz ) - Gamza = gupxx * Gamzxx + gupyy * Gamzyy + gupzz * Gamzzz + & - TWO*( gupxy * Gamzxy + gupxz * Gamzxz + gupyz * Gamzyz ) - - 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) - - Gamx_rhs = Gamx_rhs + F2o3 * Gamxa * div_beta - & - Gamxa * betaxx - Gamya * betaxy - Gamza * betaxz + & - F1o3 * (gupxx * fxx + gupxy * fxy + gupxz * fxz ) + & - gupxx * gxxx + gupyy * gyyx + gupzz * gzzx + & - TWO * (gupxy * gxyx + gupxz * gxzx + gupyz * gyzx ) - - Gamy_rhs = Gamy_rhs + F2o3 * Gamya * div_beta - & - Gamxa * betayx - Gamya * betayy - Gamza * betayz + & - F1o3 * (gupxy * fxx + gupyy * fxy + gupyz * fxz ) + & - gupxx * gxxy + gupyy * gyyy + gupzz * gzzy + & - TWO * (gupxy * gxyy + gupxz * gxzy + gupyz * gyzy ) - - Gamz_rhs = Gamz_rhs + F2o3 * Gamza * div_beta - & - Gamxa * betazx - Gamya * betazy - Gamza * betazz + & - F1o3 * (gupxz * fxx + gupyz * fxy + gupzz * fxz ) + & - gupxx * gxxz + gupyy * gyyz + gupzz * gzzz + & - TWO * (gupxy * gxyz + gupxz * gxzz + gupyz * gyzz ) !rhs for Gam^i + 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 @@ -948,15 +1000,15 @@ !!!!!!!!!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): kodis stencil coefficients sum to zero, -! so the constant offset has no effect on dissipation. - - call lopsided_kodis(ex,X,Y,Z,gxx,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,gyy,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,gzz,gzz_rhs,betax,betay,betaz,Symmetry,SSS,eps) +! 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)