64-BitBrainstorm_2026/AMSS-NCKU
8
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
ed89bc029b7945d2c5790027f679ece987506053
AMSS-NCKU/AMSS_NCKU_source/bssn_rhs_opt.f90
CGH0S7 ed89bc029b Fix potential division by zero in reta_val calculation and enable NaN checks 
Added a safety check for the denominator in the reta_val calculation to
prevent division by zero when chi approaches zero (e.g., at far-field
boundaries). Also enabled DEBUG_NAN_CHECK macro to catch invalid inputs early.
Initialized output arrays to zero to prevent uninitialized memory access.
2026-01-19 20:29:48 +08:00

1126 lines
78 KiB
Fortran
Raw Blame History

#include "macrodef.fh"
module bssn_rhs_opt_mod
implicit none
! Block Size Configuration
! 8^3 * 8 bytes * ~50 arrays ~= 200KB, fitting comfortably in L2 Cache
integer, parameter :: BLK = 8
! Physics Constants
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 :: PI = 3.14159265358979323846d0
real*8, parameter :: F1o3 = 1.D0/3.D0, F2o3 = 2.D0/3.D0, F3o2=1.5d0
real*8, parameter :: F16=1.6d1, F8=8.d0, F1o12=1.0d0/12.0d0, F1o144=1.0d0/144.0d0
real*8, parameter :: F1o4=0.25d0
real*8, parameter :: F30=3.d1, F18=1.8d1, F10=1.d1, F6=6.d0, F3=3.d0
real*8, parameter :: F64=6.4d1, FIT=1.5d1, TWT=2.d1, SIX=6.d0
real*8, parameter :: FF = 0.75d0, eta=2.d0
contains
! Main Entry Point for Optimized Computation
subroutine compute_rhs_bssn_opt(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, gont)
implicit none
! Arguments
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,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
real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Gamxxx,Gamxxy,Gamxxz,Gamxyy,Gamxyz,Gamxzz
real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Gamyxx,Gamyxy,Gamyxz,Gamyyy,Gamyyz,Gamyzz
real*8, dimension(ex(1),ex(2),ex(3)),intent(out) :: Gamzxx,Gamzxy,Gamzxz,Gamzyy,Gamzyz,Gamzzz
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
integer, intent(out) :: gont
! Local variables
integer :: kb, jb, ib, ie, je, ke
real*8 :: dX_val, dY_val, dZ_val
real*8 :: idx, idy, idz
real*8 :: dX_check
real*8, parameter :: SYM = 1.D0, ANTI= -1.D0
dX_val = X(2) - X(1)
dY_val = Y(2) - Y(1)
dZ_val = Z(2) - Z(1)
idx = ONE/dX_val
idy = ONE/dY_val
idz = ONE/dZ_val
! Initialize outputs to zero
chi_rhs = ZEO; trK_rhs = ZEO
gxx_rhs = ZEO; gxy_rhs = ZEO; gxz_rhs = ZEO
gyy_rhs = ZEO; gyz_rhs = ZEO; gzz_rhs = ZEO
Axx_rhs = ZEO; Axy_rhs = ZEO; Axz_rhs = ZEO
Ayy_rhs = ZEO; Ayz_rhs = ZEO; Azz_rhs = ZEO
Gamx_rhs = ZEO; Gamy_rhs = ZEO; Gamz_rhs = ZEO
Lap_rhs = ZEO; betax_rhs = ZEO; betay_rhs = ZEO; betaz_rhs = ZEO
dtSfx_rhs = ZEO; dtSfy_rhs = ZEO; dtSfz_rhs = ZEO
Gamxxx = ZEO; Gamxxy = ZEO; Gamxxz = ZEO
Gamxyy = ZEO; Gamxyz = ZEO; Gamxzz = ZEO
Gamyxx = ZEO; Gamyxy = ZEO; Gamyxz = ZEO
Gamyyy = ZEO; Gamyyz = ZEO; Gamyzz = ZEO
Gamzxx = ZEO; Gamzxy = ZEO; Gamzxz = ZEO
Gamzyy = ZEO; Gamzyz = ZEO; Gamzzz = ZEO
Rxx = ZEO; Rxy = ZEO; Rxz = ZEO
Ryy = ZEO; Ryz = ZEO; Rzz = ZEO
! Optimization 1: NaN Check
#define DEBUG_NAN_CHECK
#ifdef DEBUG_NAN_CHECK
dX_check = sum(chi)+sum(trK)+sum(dxx)+sum(gxy)+sum(gxz)+sum(dyy)+sum(gyz)+sum(dzz)
if(dX_check.ne.dX_check) then
write(*,*) "bssn_rhs_opt: NaN detected in input variables!"
gont = 1
return
endif
#endif
gont = 0
! Optimization 2: Blocked Execution
do kb = 1, ex(3), BLK
ke = min(kb+BLK-1, ex(3))
do jb = 1, ex(2), BLK
je = min(jb+BLK-1, ex(2))
do ib = 1, ex(1), BLK
ie = min(ib+BLK-1, ex(1))
call compute_block_kernel(ib, ie, jb, je, kb, ke)
end do
end do
end do
contains
subroutine compute_block_kernel(is, ie, js, je, ks, ke)
integer, intent(in) :: is, ie, js, je, ks, ke
integer :: nx, ny, nz
integer :: ii, jj, kk, ig, jg, kg
! Local Arrays (Stack Allocated)
! Derivatives
real*8 :: gxxx(BLK,BLK,BLK), gxxy(BLK,BLK,BLK), gxxz(BLK,BLK,BLK)
real*8 :: gxyx(BLK,BLK,BLK), gxyy(BLK,BLK,BLK), gxyz(BLK,BLK,BLK)
real*8 :: gxzx(BLK,BLK,BLK), gxzy(BLK,BLK,BLK), gxzz(BLK,BLK,BLK)
real*8 :: gyyx(BLK,BLK,BLK), gyyy(BLK,BLK,BLK), gyyz(BLK,BLK,BLK)
real*8 :: gyzx(BLK,BLK,BLK), gyzy(BLK,BLK,BLK), gyzz(BLK,BLK,BLK)
real*8 :: gzzx(BLK,BLK,BLK), gzzy(BLK,BLK,BLK), gzzz(BLK,BLK,BLK)
real*8 :: chix(BLK,BLK,BLK), chiy(BLK,BLK,BLK), chiz(BLK,BLK,BLK)
real*8 :: betaxx(BLK,BLK,BLK), betaxy(BLK,BLK,BLK), betaxz(BLK,BLK,BLK)
real*8 :: betayx(BLK,BLK,BLK), betayy(BLK,BLK,BLK), betayz(BLK,BLK,BLK)
real*8 :: betazx(BLK,BLK,BLK), betazy(BLK,BLK,BLK), betazz(BLK,BLK,BLK)
real*8 :: Lapx(BLK,BLK,BLK), Lapy(BLK,BLK,BLK), Lapz(BLK,BLK,BLK)
real*8 :: Kx(BLK,BLK,BLK), Ky(BLK,BLK,BLK), Kz(BLK,BLK,BLK)
real*8 :: Gamxx(BLK,BLK,BLK), Gamxy(BLK,BLK,BLK), Gamxz(BLK,BLK,BLK)
real*8 :: Gamyx(BLK,BLK,BLK), Gamyy(BLK,BLK,BLK), Gamyz(BLK,BLK,BLK)
real*8 :: Gamzx(BLK,BLK,BLK), Gamzy(BLK,BLK,BLK), Gamzz(BLK,BLK,BLK)
! 2nd Derivatives
real*8 :: fxx(BLK,BLK,BLK), fxy(BLK,BLK,BLK), fxz(BLK,BLK,BLK)
real*8 :: fyy(BLK,BLK,BLK), fyz(BLK,BLK,BLK), fzz(BLK,BLK,BLK)
real*8 :: f(BLK,BLK,BLK)
! Physics Intermediates
real*8 :: gupxx(BLK,BLK,BLK), gupxy(BLK,BLK,BLK), gupxz(BLK,BLK,BLK)
real*8 :: gupyy(BLK,BLK,BLK), gupyz(BLK,BLK,BLK), gupzz(BLK,BLK,BLK)
real*8 :: Rxx_l(BLK,BLK,BLK), Rxy_l(BLK,BLK,BLK), Rxz_l(BLK,BLK,BLK)
real*8 :: Ryy_l(BLK,BLK,BLK), Ryz_l(BLK,BLK,BLK), Rzz_l(BLK,BLK,BLK)
real*8 :: alpn1, chin1, gdet, igdet, div_beta
real*8 :: Gamxa, Gamya, Gamza, S_val
real*8 :: val_gxx, val_gyy, val_gzz
real*8 :: gxy2, gxz2, gyz2 ! For optimized matrix inversion
real*8 :: term_xx1, term_xx2, term_yy1, term_yy2, term_zz1, term_zz2 ! For Christoffel CSE
real*8 :: term_xy, term_xz, term_yz
real*8 :: reta_val, denom
nx = ie - is + 1
ny = je - js + 1
nz = ke - ks + 1
! === 1. Compute 1st Derivatives ===
call calc_derivs(betax, betaxx, betaxy, betaxz, ANTI, SYM, SYM, is, ie, js, je, ks, ke)
call calc_derivs(betay, betayx, betayy, betayz, SYM, ANTI, SYM, is, ie, js, je, ks, ke)
call calc_derivs(betaz, betazx, betazy, betazz, SYM, SYM, ANTI, is, ie, js, je, ks, ke)
call calc_derivs(chi, chix, chiy, chiz, SYM, SYM, SYM, is, ie, js, je, ks, ke)
call calc_derivs(dxx, gxxx, gxxy, gxxz, SYM, SYM, SYM, is, ie, js, je, ks, ke)
call calc_derivs(gxy, gxyx, gxyy, gxyz, ANTI, ANTI, SYM, is, ie, js, je, ks, ke)
call calc_derivs(gxz, gxzx, gxzy, gxzz, ANTI, SYM, ANTI, is, ie, js, je, ks, ke)
call calc_derivs(dyy, gyyx, gyyy, gyyz, SYM, SYM, SYM, is, ie, js, je, ks, ke)
call calc_derivs(gyz, gyzx, gyzy, gyzz, SYM, ANTI, ANTI, is, ie, js, je, ks, ke)
call calc_derivs(dzz, gzzx, gzzy, gzzz, SYM, SYM, SYM, is, ie, js, je, ks, ke)
call calc_derivs(Lap, Lapx, Lapy, Lapz, SYM, SYM, SYM, is, ie, js, je, ks, ke)
call calc_derivs(trK, Kx, Ky, Kz, SYM, SYM, SYM, is, ie, js, je, ks, ke)
call calc_derivs(Gamx, Gamxx, Gamxy, Gamxz, ANTI, SYM, SYM, is, ie, js, je, ks, ke)
call calc_derivs(Gamy, Gamyx, Gamyy, Gamyz, SYM, ANTI, SYM, is, ie, js, je, ks, ke)
call calc_derivs(Gamz, Gamzx, Gamzy, Gamzz, SYM, SYM, ANTI, is, ie, js, je, ks, ke)
! === 2. Compute Evolution RHS (Part 1) ===
!DIR$ SIMD
do kk=1,nz; do jj=1,ny; do ii=1,nx
ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
div_beta = betaxx(ii,jj,kk) + betayy(ii,jj,kk) + betazz(ii,jj,kk)
alpn1 = Lap(ig,jg,kg) + ONE
chin1 = chi(ig,jg,kg) + ONE
! chi_rhs
chi_rhs(ig,jg,kg) = F2o3 * chin1 * (alpn1 * trK(ig,jg,kg) - div_beta)
! Metric RHS
val_gxx = dxx(ig,jg,kg) + ONE
val_gyy = dyy(ig,jg,kg) + ONE
val_gzz = dzz(ig,jg,kg) + ONE
gxx_rhs(ig,jg,kg) = -TWO*alpn1*Axx(ig,jg,kg) - F2o3*val_gxx*div_beta + &
TWO*(val_gxx*betaxx(ii,jj,kk) + gxy(ig,jg,kg)*betayx(ii,jj,kk) + gxz(ig,jg,kg)*betazx(ii,jj,kk))
gyy_rhs(ig,jg,kg) = -TWO*alpn1*Ayy(ig,jg,kg) - F2o3*val_gyy*div_beta + &
TWO*(gxy(ig,jg,kg)*betaxy(ii,jj,kk) + val_gyy*betayy(ii,jj,kk) + gyz(ig,jg,kg)*betazy(ii,jj,kk))
gzz_rhs(ig,jg,kg) = -TWO*alpn1*Azz(ig,jg,kg) - F2o3*val_gzz*div_beta + &
TWO*(gxz(ig,jg,kg)*betaxz(ii,jj,kk) + gyz(ig,jg,kg)*betayz(ii,jj,kk) + val_gzz*betazz(ii,jj,kk))
gxy_rhs(ig,jg,kg) = -TWO*alpn1*Axy(ig,jg,kg) + F1o3*gxy(ig,jg,kg)*div_beta + &
val_gxx*betaxy(ii,jj,kk) + gxz(ig,jg,kg)*betazy(ii,jj,kk) + &
val_gyy*betayx(ii,jj,kk) + gyz(ig,jg,kg)*betazx(ii,jj,kk) - gxy(ig,jg,kg)*betazz(ii,jj,kk)
gyz_rhs(ig,jg,kg) = -TWO*alpn1*Ayz(ig,jg,kg) + F1o3*gyz(ig,jg,kg)*div_beta + &
gxy(ig,jg,kg)*betaxz(ii,jj,kk) + val_gyy*betayz(ii,jj,kk) + &
gxz(ig,jg,kg)*betaxy(ii,jj,kk) + val_gzz*betazy(ii,jj,kk) - gyz(ig,jg,kg)*betaxx(ii,jj,kk)
gxz_rhs(ig,jg,kg) = -TWO*alpn1*Axz(ig,jg,kg) + F1o3*gxz(ig,jg,kg)*div_beta + &
val_gxx*betaxz(ii,jj,kk) + gxy(ig,jg,kg)*betayz(ii,jj,kk) + &
gyz(ig,jg,kg)*betayx(ii,jj,kk) + val_gzz*betazx(ii,jj,kk) - gxz(ig,jg,kg)*betayy(ii,jj,kk)
! Inverse Metric (Optimization 4: Optimized 3x3 matrix inversion)
! Precompute squared terms to avoid redundant multiplications
gxy2 = gxy(ig,jg,kg)**2
gxz2 = gxz(ig,jg,kg)**2
gyz2 = gyz(ig,jg,kg)**2
! Compute determinant using precomputed squares
gdet = val_gxx*val_gyy*val_gzz + TWO*gxy(ig,jg,kg)*gxz(ig,jg,kg)*gyz(ig,jg,kg) - &
val_gxx*gyz2 - val_gyy*gxz2 - val_gzz*gxy2
igdet = ONE / gdet
! Compute cofactors and multiply by inverse determinant
gupxx(ii,jj,kk) = (val_gyy*val_gzz - gyz2) * igdet
gupyy(ii,jj,kk) = (val_gxx*val_gzz - gxz2) * igdet
gupzz(ii,jj,kk) = (val_gxx*val_gyy - gxy2) * igdet
gupxy(ii,jj,kk) = (gxz(ig,jg,kg)*gyz(ig,jg,kg) - gxy(ig,jg,kg)*val_gzz) * igdet
gupxz(ii,jj,kk) = (gxy(ig,jg,kg)*gyz(ig,jg,kg) - gxz(ig,jg,kg)*val_gyy) * igdet
gupyz(ii,jj,kk) = (gxy(ig,jg,kg)*gxz(ig,jg,kg) - gyz(ig,jg,kg)*val_gxx) * igdet
! Compute Gam^i Residual (Constraint)
if(co == 0) then
Gmx_Res(ig,jg,kg) = Gamx(ig,jg,kg) - ( &
gupxx(ii,jj,kk)*(gupxx(ii,jj,kk)*gxxx(ii,jj,kk)+gupxy(ii,jj,kk)*gxyx(ii,jj,kk)+gupxz(ii,jj,kk)*gxzx(ii,jj,kk)) + &
gupxy(ii,jj,kk)*(gupxx(ii,jj,kk)*gxyx(ii,jj,kk)+gupxy(ii,jj,kk)*gyyx(ii,jj,kk)+gupxz(ii,jj,kk)*gyzx(ii,jj,kk)) + &
gupxz(ii,jj,kk)*(gupxx(ii,jj,kk)*gxzx(ii,jj,kk)+gupxy(ii,jj,kk)*gyzx(ii,jj,kk)+gupxz(ii,jj,kk)*gzzx(ii,jj,kk)) + &
gupxx(ii,jj,kk)*(gupxy(ii,jj,kk)*gxxy(ii,jj,kk)+gupyy(ii,jj,kk)*gxyy(ii,jj,kk)+gupyz(ii,jj,kk)*gxzy(ii,jj,kk)) + &
gupxy(ii,jj,kk)*(gupxy(ii,jj,kk)*gxyy(ii,jj,kk)+gupyy(ii,jj,kk)*gyyy(ii,jj,kk)+gupyz(ii,jj,kk)*gyzy(ii,jj,kk)) + &
gupxz(ii,jj,kk)*(gupxy(ii,jj,kk)*gxzy(ii,jj,kk)+gupyy(ii,jj,kk)*gyzy(ii,jj,kk)+gupyz(ii,jj,kk)*gzzy(ii,jj,kk)) + &
gupxx(ii,jj,kk)*(gupxz(ii,jj,kk)*gxxz(ii,jj,kk)+gupyz(ii,jj,kk)*gxyz(ii,jj,kk)+gupzz(ii,jj,kk)*gxzz(ii,jj,kk)) + &
gupxy(ii,jj,kk)*(gupxz(ii,jj,kk)*gxyz(ii,jj,kk)+gupyz(ii,jj,kk)*gyyz(ii,jj,kk)+gupzz(ii,jj,kk)*gyzz(ii,jj,kk)) + &
gupxz(ii,jj,kk)*(gupxz(ii,jj,kk)*gxzz(ii,jj,kk)+gupyz(ii,jj,kk)*gyzz(ii,jj,kk)+gupzz(ii,jj,kk)*gzzz(ii,jj,kk)) )
end if
end do; end do; end do
! === 3. Christoffel Symbols (Optimization 5: CSE) ===
!DIR$ SIMD
do kk=1,nz; do jj=1,ny; do ii=1,nx
ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
term_xx1 = TWO*gxyx(ii,jj,kk) - gxxy(ii,jj,kk)
term_xx2 = TWO*gxzx(ii,jj,kk) - gxxz(ii,jj,kk)
term_yy1 = TWO*gxyy(ii,jj,kk) - gyyx(ii,jj,kk)
term_yy2 = TWO*gyzy(ii,jj,kk) - gyyz(ii,jj,kk)
term_zz1 = TWO*gxzz(ii,jj,kk) - gzzx(ii,jj,kk)
term_zz2 = TWO*gyzz(ii,jj,kk) - gzzy(ii,jj,kk)
term_xy = gxzy(ii,jj,kk) + gyzx(ii,jj,kk) - gxyz(ii,jj,kk)
term_xz = gxyz(ii,jj,kk) + gyzx(ii,jj,kk) - gxzy(ii,jj,kk)
term_yz = gxyz(ii,jj,kk) + gxzy(ii,jj,kk) - gyzx(ii,jj,kk)
Gamxxx(ig,jg,kg) = HALF*(gupxx(ii,jj,kk)*gxxx(ii,jj,kk) + gupxy(ii,jj,kk)*term_xx1 + gupxz(ii,jj,kk)*term_xx2)
Gamyxx(ig,jg,kg) = HALF*(gupxy(ii,jj,kk)*gxxx(ii,jj,kk) + gupyy(ii,jj,kk)*term_xx1 + gupyz(ii,jj,kk)*term_xx2)
Gamzxx(ig,jg,kg) = HALF*(gupxz(ii,jj,kk)*gxxx(ii,jj,kk) + gupyz(ii,jj,kk)*term_xx1 + gupzz(ii,jj,kk)*term_xx2)
Gamxyy(ig,jg,kg) = HALF*(gupxx(ii,jj,kk)*term_yy1 + gupxy(ii,jj,kk)*gyyy(ii,jj,kk) + gupxz(ii,jj,kk)*term_yy2)
Gamyyy(ig,jg,kg) = HALF*(gupxy(ii,jj,kk)*term_yy1 + gupyy(ii,jj,kk)*gyyy(ii,jj,kk) + gupyz(ii,jj,kk)*term_yy2)
Gamzyy(ig,jg,kg) = HALF*(gupxz(ii,jj,kk)*term_yy1 + gupyz(ii,jj,kk)*gyyy(ii,jj,kk) + gupzz(ii,jj,kk)*term_yy2)
Gamxzz(ig,jg,kg) = HALF*(gupxx(ii,jj,kk)*term_zz1 + gupxy(ii,jj,kk)*term_zz2 + gupxz(ii,jj,kk)*gzzz(ii,jj,kk))
Gamyzz(ig,jg,kg) = HALF*(gupxy(ii,jj,kk)*term_zz1 + gupyy(ii,jj,kk)*term_zz2 + gupyz(ii,jj,kk)*gzzz(ii,jj,kk))
Gamzzz(ig,jg,kg) = HALF*(gupxz(ii,jj,kk)*term_zz1 + gupyz(ii,jj,kk)*term_zz2 + gupzz(ii,jj,kk)*gzzz(ii,jj,kk))
Gamxxy(ig,jg,kg) = HALF*(gupxx(ii,jj,kk)*gxxy(ii,jj,kk) + gupxy(ii,jj,kk)*gyyx(ii,jj,kk) + gupxz(ii,jj,kk)*term_xy)
Gamyxy(ig,jg,kg) = HALF*(gupxy(ii,jj,kk)*gxxy(ii,jj,kk) + gupyy(ii,jj,kk)*gyyx(ii,jj,kk) + gupyz(ii,jj,kk)*term_xy)
Gamzxy(ig,jg,kg) = HALF*(gupxz(ii,jj,kk)*gxxy(ii,jj,kk) + gupyz(ii,jj,kk)*gyyx(ii,jj,kk) + gupzz(ii,jj,kk)*term_xy)
Gamxxz(ig,jg,kg) = HALF*(gupxx(ii,jj,kk)*gxxz(ii,jj,kk) + gupxy(ii,jj,kk)*term_xz + gupxz(ii,jj,kk)*gzzx(ii,jj,kk))
Gamyxz(ig,jg,kg) = HALF*(gupxy(ii,jj,kk)*gxxz(ii,jj,kk) + gupyy(ii,jj,kk)*term_xz + gupyz(ii,jj,kk)*gzzx(ii,jj,kk))
Gamzxz(ig,jg,kg) = HALF*(gupxz(ii,jj,kk)*gxxz(ii,jj,kk) + gupyz(ii,jj,kk)*term_xz + gupzz(ii,jj,kk)*gzzx(ii,jj,kk))
Gamxyz(ig,jg,kg) = HALF*(gupxx(ii,jj,kk)*term_yz + gupxy(ii,jj,kk)*gyyz(ii,jj,kk) + gupxz(ii,jj,kk)*gzzy(ii,jj,kk))
Gamyyz(ig,jg,kg) = HALF*(gupxy(ii,jj,kk)*term_yz + gupyy(ii,jj,kk)*gyyz(ii,jj,kk) + gupyz(ii,jj,kk)*gzzy(ii,jj,kk))
Gamzyz(ig,jg,kg) = HALF*(gupxz(ii,jj,kk)*term_yz + gupyz(ii,jj,kk)*gyyz(ii,jj,kk) + gupzz(ii,jj,kk)*gzzy(ii,jj,kk))
end do; end do; end do
! === 4. Ricci Tensor (Part 1 - Terms from A_ij) ===
!DIR$ SIMD
do kk=1,nz; do jj=1,ny; do ii=1,nx
ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
Rxx_l(ii,jj,kk) = gupxx(ii,jj,kk)**2*Axx(ig,jg,kg) + gupxy(ii,jj,kk)**2*Ayy(ig,jg,kg) + gupxz(ii,jj,kk)**2*Azz(ig,jg,kg) + &
TWO*(gupxx(ii,jj,kk)*gupxy(ii,jj,kk)*Axy(ig,jg,kg) + gupxx(ii,jj,kk)*gupxz(ii,jj,kk)*Axz(ig,jg,kg) + gupxy(ii,jj,kk)*gupxz(ii,jj,kk)*Ayz(ig,jg,kg))
Ryy_l(ii,jj,kk) = gupxy(ii,jj,kk)**2*Axx(ig,jg,kg) + gupyy(ii,jj,kk)**2*Ayy(ig,jg,kg) + gupyz(ii,jj,kk)**2*Azz(ig,jg,kg) + &
TWO*(gupxy(ii,jj,kk)*gupyy(ii,jj,kk)*Axy(ig,jg,kg) + gupxy(ii,jj,kk)*gupyz(ii,jj,kk)*Axz(ig,jg,kg) + gupyy(ii,jj,kk)*gupyz(ii,jj,kk)*Ayz(ig,jg,kg))
Rzz_l(ii,jj,kk) = gupxz(ii,jj,kk)**2*Axx(ig,jg,kg) + gupyz(ii,jj,kk)**2*Ayy(ig,jg,kg) + gupzz(ii,jj,kk)**2*Azz(ig,jg,kg) + &
TWO*(gupxz(ii,jj,kk)*gupyz(ii,jj,kk)*Axy(ig,jg,kg) + gupxz(ii,jj,kk)*gupzz(ii,jj,kk)*Axz(ig,jg,kg) + gupyz(ii,jj,kk)*gupzz(ii,jj,kk)*Ayz(ig,jg,kg))
Rxy_l(ii,jj,kk) = gupxx(ii,jj,kk)*gupxy(ii,jj,kk)*Axx(ig,jg,kg) + gupxy(ii,jj,kk)*gupyy(ii,jj,kk)*Ayy(ig,jg,kg) + gupxz(ii,jj,kk)*gupyz(ii,jj,kk)*Azz(ig,jg,kg) + &
(gupxx(ii,jj,kk)*gupyy(ii,jj,kk) + gupxy(ii,jj,kk)**2)*Axy(ig,jg,kg) + &
(gupxx(ii,jj,kk)*gupyz(ii,jj,kk) + gupxz(ii,jj,kk)*gupxy(ii,jj,kk))*Axz(ig,jg,kg) + &
(gupxy(ii,jj,kk)*gupyz(ii,jj,kk) + gupxz(ii,jj,kk)*gupyy(ii,jj,kk))*Ayz(ig,jg,kg)
Rxz_l(ii,jj,kk) = gupxx(ii,jj,kk)*gupxz(ii,jj,kk)*Axx(ig,jg,kg) + gupxy(ii,jj,kk)*gupyz(ii,jj,kk)*Ayy(ig,jg,kg) + gupxz(ii,jj,kk)*gupzz(ii,jj,kk)*Azz(ig,jg,kg) + &
(gupxx(ii,jj,kk)*gupyz(ii,jj,kk) + gupxy(ii,jj,kk)*gupxz(ii,jj,kk))*Axy(ig,jg,kg) + &
(gupxx(ii,jj,kk)*gupzz(ii,jj,kk) + gupxz(ii,jj,kk)**2)*Axz(ig,jg,kg) + &
(gupxy(ii,jj,kk)*gupzz(ii,jj,kk) + gupxz(ii,jj,kk)*gupyz(ii,jj,kk))*Ayz(ig,jg,kg)
Ryz_l(ii,jj,kk) = gupxy(ii,jj,kk)*gupxz(ii,jj,kk)*Axx(ig,jg,kg) + gupyy(ii,jj,kk)*gupyz(ii,jj,kk)*Ayy(ig,jg,kg) + gupyz(ii,jj,kk)*gupzz(ii,jj,kk)*Azz(ig,jg,kg) + &
(gupxy(ii,jj,kk)*gupyz(ii,jj,kk) + gupyy(ii,jj,kk)*gupxz(ii,jj,kk))*Axy(ig,jg,kg) + &
(gupxy(ii,jj,kk)*gupzz(ii,jj,kk) + gupyz(ii,jj,kk)*gupxz(ii,jj,kk))*Axz(ig,jg,kg) + &
(gupyy(ii,jj,kk)*gupzz(ii,jj,kk) + gupyz(ii,jj,kk)**2)*Ayz(ig,jg,kg)
end do; end do; end do
! === 5. Gamma^i RHS ===
!DIR$ SIMD
do kk=1,nz; do jj=1,ny; do ii=1,nx
ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
alpn1 = Lap(ig,jg,kg) + ONE
chin1 = chi(ig,jg,kg) + ONE
Gamx_rhs(ig,jg,kg) = -TWO*(Lapx(ii,jj,kk)*Rxx_l(ii,jj,kk) + Lapy(ii,jj,kk)*Rxy_l(ii,jj,kk) + Lapz(ii,jj,kk)*Rxz_l(ii,jj,kk)) + &
TWO*alpn1*( -F3o2/chin1*(chix(ii,jj,kk)*Rxx_l(ii,jj,kk) + chiy(ii,jj,kk)*Rxy_l(ii,jj,kk) + chiz(ii,jj,kk)*Rxz_l(ii,jj,kk)) - &
gupxx(ii,jj,kk)*(F2o3*Kx(ii,jj,kk) + EIGHT*PI*Sx(ig,jg,kg)) - &
gupxy(ii,jj,kk)*(F2o3*Ky(ii,jj,kk) + EIGHT*PI*Sy(ig,jg,kg)) - &
gupxz(ii,jj,kk)*(F2o3*Kz(ii,jj,kk) + EIGHT*PI*Sz(ig,jg,kg)) + &
Gamxxx(ig,jg,kg)*Rxx_l(ii,jj,kk) + Gamxyy(ig,jg,kg)*Ryy_l(ii,jj,kk) + Gamxzz(ig,jg,kg)*Rzz_l(ii,jj,kk) + &
TWO*(Gamxxy(ig,jg,kg)*Rxy_l(ii,jj,kk) + Gamxxz(ig,jg,kg)*Rxz_l(ii,jj,kk) + Gamxyz(ig,jg,kg)*Ryz_l(ii,jj,kk)) )
Gamy_rhs(ig,jg,kg) = -TWO*(Lapx(ii,jj,kk)*Rxy_l(ii,jj,kk) + Lapy(ii,jj,kk)*Ryy_l(ii,jj,kk) + Lapz(ii,jj,kk)*Ryz_l(ii,jj,kk)) + &
TWO*alpn1*( -F3o2/chin1*(chix(ii,jj,kk)*Rxy_l(ii,jj,kk) + chiy(ii,jj,kk)*Ryy_l(ii,jj,kk) + chiz(ii,jj,kk)*Ryz_l(ii,jj,kk)) - &
gupxy(ii,jj,kk)*(F2o3*Kx(ii,jj,kk) + EIGHT*PI*Sx(ig,jg,kg)) - &
gupyy(ii,jj,kk)*(F2o3*Ky(ii,jj,kk) + EIGHT*PI*Sy(ig,jg,kg)) - &
gupyz(ii,jj,kk)*(F2o3*Kz(ii,jj,kk) + EIGHT*PI*Sz(ig,jg,kg)) + &
Gamyxx(ig,jg,kg)*Rxx_l(ii,jj,kk) + Gamyyy(ig,jg,kg)*Ryy_l(ii,jj,kk) + Gamyzz(ig,jg,kg)*Rzz_l(ii,jj,kk) + &
TWO*(Gamyxy(ig,jg,kg)*Rxy_l(ii,jj,kk) + Gamyxz(ig,jg,kg)*Rxz_l(ii,jj,kk) + Gamyyz(ig,jg,kg)*Ryz_l(ii,jj,kk)) )
Gamz_rhs(ig,jg,kg) = -TWO*(Lapx(ii,jj,kk)*Rxz_l(ii,jj,kk) + Lapy(ii,jj,kk)*Ryz_l(ii,jj,kk) + Lapz(ii,jj,kk)*Rzz_l(ii,jj,kk)) + &
TWO*alpn1*( -F3o2/chin1*(chix(ii,jj,kk)*Rxz_l(ii,jj,kk) + chiy(ii,jj,kk)*Ryz_l(ii,jj,kk) + chiz(ii,jj,kk)*Rzz_l(ii,jj,kk)) - &
gupxz(ii,jj,kk)*(F2o3*Kx(ii,jj,kk) + EIGHT*PI*Sx(ig,jg,kg)) - &
gupyz(ii,jj,kk)*(F2o3*Ky(ii,jj,kk) + EIGHT*PI*Sy(ig,jg,kg)) - &
gupzz(ii,jj,kk)*(F2o3*Kz(ii,jj,kk) + EIGHT*PI*Sz(ig,jg,kg)) + &
Gamzxx(ig,jg,kg)*Rxx_l(ii,jj,kk) + Gamzyy(ig,jg,kg)*Ryy_l(ii,jj,kk) + Gamzzz(ig,jg,kg)*Rzz_l(ii,jj,kk) + &
TWO*(Gamzxy(ig,jg,kg)*Rxy_l(ii,jj,kk) + Gamzxz(ig,jg,kg)*Rxz_l(ii,jj,kk) + Gamzyz(ig,jg,kg)*Ryz_l(ii,jj,kk)) )
end do; end do; end do
! === 6. Second Derivatives for Ricci ===
! Need 2nd derivatives of gxx etc.
call calc_dderivs(dxx, fxx, fxy, fxz, fyy, fyz, fzz, SYM, SYM, SYM, is, ie, js, je, ks, ke)
!DIR$ SIMD
do kk=1,nz; do jj=1,ny; do ii=1,nx
ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
val_gxx = dxx(ig,jg,kg)+ONE
val_gyy = dyy(ig,jg,kg)+ONE
val_gzz = dzz(ig,jg,kg)+ONE
Rxx_l(ii,jj,kk) = -HALF*Rxx_l(ii,jj,kk) + val_gxx*Gamxx(ii,jj,kk) + gxy(ig,jg,kg)*Gamyx(ii,jj,kk) + gxz(ig,jg,kg)*Gamzx(ii,jj,kk) + &
(gupxx(ii,jj,kk)*Gamxxx(ig,jg,kg)+gupyy(ii,jj,kk)*Gamxyy(ig,jg,kg)+gupzz(ii,jj,kk)*Gamxzz(ig,jg,kg) + &
TWO*(gupxy(ii,jj,kk)*Gamxxy(ig,jg,kg)+gupxz(ii,jj,kk)*Gamxxz(ig,jg,kg)+gupyz(ii,jj,kk)*Gamxyz(ig,jg,kg)))*gxxx(ii,jj,kk) + &
(gupxx(ii,jj,kk)*Gamyxx(ig,jg,kg)+gupyy(ii,jj,kk)*Gamyyy(ig,jg,kg)+gupzz(ii,jj,kk)*Gamyzz(ig,jg,kg) + &
TWO*(gupxy(ii,jj,kk)*Gamyxy(ig,jg,kg)+gupxz(ii,jj,kk)*Gamyxz(ig,jg,kg)+gupyz(ii,jj,kk)*Gamyyz(ig,jg,kg)))*gxyx(ii,jj,kk) + &
(gupxx(ii,jj,kk)*Gamzxx(ig,jg,kg)+gupyy(ii,jj,kk)*Gamzyy(ig,jg,kg)+gupzz(ii,jj,kk)*Gamzzz(ig,jg,kg) + &
TWO*(gupxy(ii,jj,kk)*Gamzxy(ig,jg,kg)+gupxz(ii,jj,kk)*Gamzxz(ig,jg,kg)+gupyz(ii,jj,kk)*Gamzyz(ig,jg,kg)))*gxzx(ii,jj,kk) + &
gupxx(ii,jj,kk)*(fxx(ii,jj,kk) + TWO*(Gamxxx(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gxyx(ii,jj,kk) + Gamzxx(ig,jg,kg)*gxzx(ii,jj,kk)) + Gamxxx(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gxxy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gxxz(ii,jj,kk)) + &
gupxy(ii,jj,kk)*(fxy(ii,jj,kk) + TWO*(Gamxxx(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gyyx(ii,jj,kk) + Gamzxx(ig,jg,kg)*gyzx(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxyx(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxzx(ii,jj,kk)) + Gamxxy(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxxy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxxz(ii,jj,kk) + Gamxxx(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gxyz(ii,jj,kk)) + &
gupxz(ii,jj,kk)*(fxz(ii,jj,kk) + TWO*(Gamxxx(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gyzx(ii,jj,kk) + Gamzxx(ig,jg,kg)*gzzx(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxzx(ii,jj,kk)) + Gamxxz(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxxy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxxz(ii,jj,kk) + Gamxxx(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gxzz(ii,jj,kk)) + &
gupyy(ii,jj,kk)*(fyy(ii,jj,kk) + TWO*(Gamxxy(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyyx(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyzx(ii,jj,kk)) + Gamxxy(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxyz(ii,jj,kk)) + &
gupyz(ii,jj,kk)*(fyz(ii,jj,kk) + TWO*(Gamxxy(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyzx(ii,jj,kk) + Gamzxy(ig,jg,kg)*gzzx(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyyx(ii,jj,kk) + Gamzxz(ig,jg,kg)*gyzx(ii,jj,kk)) + Gamxxz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxzz(ii,jj,kk)) + &
gupzz(ii,jj,kk)*(fzz(ii,jj,kk) + TWO*(Gamxxz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamzxz(ig,jg,kg)*gzzx(ii,jj,kk)) + Gamxxz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxzz(ii,jj,kk))
end do; end do; end do
call calc_dderivs(dyy, fxx, fxy, fxz, fyy, fyz, fzz, SYM, SYM, SYM, is, ie, js, je, ks, ke)
!DIR$ SIMD
do kk=1,nz; do jj=1,ny; do ii=1,nx
ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
val_gyy = dyy(ig,jg,kg)+ONE
Ryy_l(ii,jj,kk) = -HALF*Ryy_l(ii,jj,kk) + gxy(ig,jg,kg)*Gamxy(ii,jj,kk) + val_gyy*Gamyy(ii,jj,kk) + gyz(ig,jg,kg)*Gamzy(ii,jj,kk) + &
(gupxx(ii,jj,kk)*Gamxxx(ig,jg,kg)+gupyy(ii,jj,kk)*Gamxyy(ig,jg,kg)+gupzz(ii,jj,kk)*Gamxzz(ig,jg,kg) + &
TWO*(gupxy(ii,jj,kk)*Gamxxy(ig,jg,kg)+gupxz(ii,jj,kk)*Gamxxz(ig,jg,kg)+gupyz(ii,jj,kk)*Gamxyz(ig,jg,kg)))*gxyy(ii,jj,kk) + &
(gupxx(ii,jj,kk)*Gamyxx(ig,jg,kg)+gupyy(ii,jj,kk)*Gamyyy(ig,jg,kg)+gupzz(ii,jj,kk)*Gamyzz(ig,jg,kg) + &
TWO*(gupxy(ii,jj,kk)*Gamyxy(ig,jg,kg)+gupxz(ii,jj,kk)*Gamyxz(ig,jg,kg)+gupyz(ii,jj,kk)*Gamyyz(ig,jg,kg)))*gyyy(ii,jj,kk) + &
(gupxx(ii,jj,kk)*Gamzxx(ig,jg,kg)+gupyy(ii,jj,kk)*Gamzyy(ig,jg,kg)+gupzz(ii,jj,kk)*Gamzzz(ig,jg,kg) + &
TWO*(gupxy(ii,jj,kk)*Gamzxy(ig,jg,kg)+gupxz(ii,jj,kk)*Gamzxz(ig,jg,kg)+gupyz(ii,jj,kk)*Gamzyz(ig,jg,kg)))*gyzy(ii,jj,kk) + &
gupxx(ii,jj,kk)*(fxx(ii,jj,kk) + TWO*(Gamxxy(ig,jg,kg)*gxxy(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxzy(ii,jj,kk)) + Gamxxy(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxyz(ii,jj,kk)) + &
gupxy(ii,jj,kk)*(fxy(ii,jj,kk) + TWO*(Gamxxy(ig,jg,kg)*gxyy(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyzy(ii,jj,kk) + Gamxyy(ig,jg,kg)*gxxy(ii,jj,kk) + Gamyyy(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzyy(ig,jg,kg)*gxzy(ii,jj,kk)) + Gamxyy(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyyy(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzyy(ig,jg,kg)*gxyz(ii,jj,kk) + Gamxxy(ig,jg,kg)*gyyx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyyz(ii,jj,kk)) + &
gupxz(ii,jj,kk)*(fxz(ii,jj,kk) + TWO*(Gamxxy(ig,jg,kg)*gxzy(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gzzy(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxxy(ii,jj,kk) + Gamyyz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gxzy(ii,jj,kk)) + Gamxyz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamxxy(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyzz(ii,jj,kk)) + &
gupyy(ii,jj,kk)*(fyy(ii,jj,kk) + TWO*(Gamxyy(ig,jg,kg)*gxyy(ii,jj,kk) + Gamyyy(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzyy(ig,jg,kg)*gyzy(ii,jj,kk)) + Gamxyy(ig,jg,kg)*gyyx(ii,jj,kk) + Gamyyy(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzyy(ig,jg,kg)*gyyz(ii,jj,kk)) + &
gupyz(ii,jj,kk)*(fyz(ii,jj,kk) + TWO*(Gamxyy(ig,jg,kg)*gxzy(ii,jj,kk) + Gamyyy(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzyy(ig,jg,kg)*gzzy(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gyzy(ii,jj,kk)) + Gamxyz(ig,jg,kg)*gyyx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gyyz(ii,jj,kk) + Gamxyy(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyyy(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzyy(ig,jg,kg)*gyzz(ii,jj,kk)) + &
gupzz(ii,jj,kk)*(fzz(ii,jj,kk) + TWO*(Gamxyz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gzzy(ii,jj,kk)) + Gamxyz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gyzz(ii,jj,kk))
end do; end do; end do
call calc_dderivs(dzz, fxx, fxy, fxz, fyy, fyz, fzz, SYM, SYM, SYM, is, ie, js, je, ks, ke)
!DIR$ SIMD
do kk=1,nz; do jj=1,ny; do ii=1,nx
ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
val_gzz = dzz(ig,jg,kg)+ONE
Rzz_l(ii,jj,kk) = -HALF*Rzz_l(ii,jj,kk) + gxz(ig,jg,kg)*Gamxz(ii,jj,kk) + gyz(ig,jg,kg)*Gamyz(ii,jj,kk) + val_gzz*Gamzz(ii,jj,kk) + &
(gupxx(ii,jj,kk)*Gamxxx(ig,jg,kg)+gupyy(ii,jj,kk)*Gamxyy(ig,jg,kg)+gupzz(ii,jj,kk)*Gamxzz(ig,jg,kg) + &
TWO*(gupxy(ii,jj,kk)*Gamxxy(ig,jg,kg)+gupxz(ii,jj,kk)*Gamxxz(ig,jg,kg)+gupyz(ii,jj,kk)*Gamxyz(ig,jg,kg)))*gxzz(ii,jj,kk) + &
(gupxx(ii,jj,kk)*Gamyxx(ig,jg,kg)+gupyy(ii,jj,kk)*Gamyyy(ig,jg,kg)+gupzz(ii,jj,kk)*Gamyzz(ig,jg,kg) + &
TWO*(gupxy(ii,jj,kk)*Gamyxy(ig,jg,kg)+gupxz(ii,jj,kk)*Gamyxz(ig,jg,kg)+gupyz(ii,jj,kk)*Gamyyz(ig,jg,kg)))*gyzz(ii,jj,kk) + &
(gupxx(ii,jj,kk)*Gamzxx(ig,jg,kg)+gupyy(ii,jj,kk)*Gamzyy(ig,jg,kg)+gupzz(ii,jj,kk)*Gamzzz(ig,jg,kg) + &
TWO*(gupxy(ii,jj,kk)*Gamzxy(ig,jg,kg)+gupxz(ii,jj,kk)*Gamzxz(ig,jg,kg)+gupyz(ii,jj,kk)*Gamzyz(ig,jg,kg)))*gzzz(ii,jj,kk) + &
gupxx(ii,jj,kk)*(fxx(ii,jj,kk) + TWO*(Gamxxz(ig,jg,kg)*gxxz(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxzz(ii,jj,kk)) + Gamxxz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxzz(ii,jj,kk)) + &
gupxy(ii,jj,kk)*(fxy(ii,jj,kk) + TWO*(Gamxxz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyyz(ii,jj,kk) + Gamzxz(ig,jg,kg)*gyzz(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxxz(ii,jj,kk) + Gamyyz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamzyz(ig,jg,kg)*gxzz(ii,jj,kk)) + Gamxyz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gyzz(ii,jj,kk)) + &
gupxz(ii,jj,kk)*(fxz(ii,jj,kk) + TWO*(Gamxxz(ig,jg,kg)*gxzz(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyzz(ii,jj,kk) + Gamzxz(ig,jg,kg)*gzzz(ii,jj,kk) + Gamxzz(ig,jg,kg)*gxxz(ii,jj,kk) + Gamyzz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamzzz(ig,jg,kg)*gxzz(ii,jj,kk)) + Gamxzz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyzz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzzz(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gzzx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gzzy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gzzz(ii,jj,kk)) + &
gupyy(ii,jj,kk)*(fyy(ii,jj,kk) + TWO*(Gamxyz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyyz(ii,jj,kk) + Gamzyz(ig,jg,kg)*gyzz(ii,jj,kk)) + Gamxyz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gyzz(ii,jj,kk)) + &
gupyz(ii,jj,kk)*(fyz(ii,jj,kk) + TWO*(Gamxyz(ig,jg,kg)*gxzz(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyzz(ii,jj,kk) + Gamzyz(ig,jg,kg)*gzzz(ii,jj,kk) + Gamxzz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamyzz(ig,jg,kg)*gyyz(ii,jj,kk) + Gamzzz(ig,jg,kg)*gyzz(ii,jj,kk)) + Gamxzz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyzz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzzz(ig,jg,kg)*gyzz(ii,jj,kk) + Gamxyz(ig,jg,kg)*gzzx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gzzy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gzzz(ii,jj,kk)) + &
gupzz(ii,jj,kk)*(fzz(ii,jj,kk) + TWO*(Gamxzz(ig,jg,kg)*gxzz(ii,jj,kk) + Gamyzz(ig,jg,kg)*gyzz(ii,jj,kk) + Gamzzz(ig,jg,kg)*gzzz(ii,jj,kk)) + Gamxzz(ig,jg,kg)*gzzx(ii,jj,kk) + Gamyzz(ig,jg,kg)*gzzy(ii,jj,kk) + Gamzzz(ig,jg,kg)*gzzz(ii,jj,kk))
end do; end do; end do
call calc_dderivs(gxy, fxx, fxy, fxz, fyy, fyz, fzz, ANTI, ANTI, SYM, is, ie, js, je, ks, ke)
!DIR$ SIMD
do kk=1,nz; do jj=1,ny; do ii=1,nx
ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
! Compute Gam^i at local point
Gamxa = Gamx(ig,jg,kg)
Gamya = Gamy(ig,jg,kg)
Gamza = Gamz(ig,jg,kg)
val_gxx = dxx(ig,jg,kg)+ONE
val_gyy = dyy(ig,jg,kg)+ONE
Rxy_l(ii,jj,kk) = HALF*( -Rxy_l(ii,jj,kk) + val_gxx*Gamxy(ii,jj,kk) + gxy(ig,jg,kg)*Gamyy(ii,jj,kk) + gxz(ig,jg,kg)*Gamzy(ii,jj,kk) + &
gxy(ig,jg,kg)*Gamxx(ii,jj,kk) + val_gyy*Gamyx(ii,jj,kk) + gyz(ig,jg,kg)*Gamzx(ii,jj,kk) + &
Gamxa*gxyx(ii,jj,kk) + Gamya*gyyx(ii,jj,kk) + Gamza*gyzx(ii,jj,kk) + &
Gamxa*gxxy(ii,jj,kk) + Gamya*gxyy(ii,jj,kk) + Gamza*gxzy(ii,jj,kk) ) + &
gupxx(ii,jj,kk)*(Gamxxx(ig,jg,kg)*gxxy(ii,jj,kk) + Gamyxx(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gxzy(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxyx(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxzx(ii,jj,kk) + Gamxxx(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gxyz(ii,jj,kk)) + &
gupxy(ii,jj,kk)*(fxy(ii,jj,kk) + Gamxxx(ig,jg,kg)*gxyy(ii,jj,kk) + Gamyxx(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gyzy(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyyx(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyzx(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxyz(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxxy(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxzy(ii,jj,kk) + Gamxyy(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyyy(ig,jg,kg)*gxyx(ii,jj,kk) + Gamzyy(ig,jg,kg)*gxzx(ii,jj,kk) + Gamxxx(ig,jg,kg)*gyyx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gyyz(ii,jj,kk)) + &
gupxz(ii,jj,kk)*(fxz(ii,jj,kk) + Gamxxx(ig,jg,kg)*gxzy(ii,jj,kk) + Gamyxx(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gzzy(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyzx(ii,jj,kk) + Gamzxy(ig,jg,kg)*gzzx(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxxy(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamzyz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamxxx(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gyzz(ii,jj,kk)) + &
gupyy(ii,jj,kk)*(Gamxxy(ig,jg,kg)*gxyy(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyzy(ii,jj,kk) + Gamxyy(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyyy(ig,jg,kg)*gyyx(ii,jj,kk) + Gamzyy(ig,jg,kg)*gyzx(ii,jj,kk) + Gamxxy(ig,jg,kg)*gyyx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyyz(ii,jj,kk)) + &
gupyz(ii,jj,kk)*(fyz(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxzy(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gzzy(ii,jj,kk) + Gamxyy(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyyy(ig,jg,kg)*gyzx(ii,jj,kk) + Gamzyy(ig,jg,kg)*gzzx(ii,jj,kk) + Gamxxz(ig,jg,kg)*gyyx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gyyz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyyx(ii,jj,kk) + Gamzyz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamxxy(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyzz(ii,jj,kk)) + &
gupzz(ii,jj,kk)*(Gamxxz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gzzy(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamzyz(ig,jg,kg)*gzzx(ii,jj,kk) + Gamxxz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gyzz(ii,jj,kk))
end do; end do; end do
call calc_dderivs(gxz, fxx, fxy, fxz, fyy, fyz, fzz, ANTI, SYM, ANTI, is, ie, js, je, ks, ke)
!DIR$ SIMD
do kk=1,nz; do jj=1,ny; do ii=1,nx
ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
! Compute Gam^i at local point
Gamxa = Gamx(ig,jg,kg)
Gamya = Gamy(ig,jg,kg)
Gamza = Gamz(ig,jg,kg)
val_gxx = dxx(ig,jg,kg)+ONE
val_gzz = dzz(ig,jg,kg)+ONE
Rxz_l(ii,jj,kk) = HALF*( -Rxz_l(ii,jj,kk) + val_gxx*Gamxz(ii,jj,kk) + gxy(ig,jg,kg)*Gamyz(ii,jj,kk) + gxz(ig,jg,kg)*Gamzz(ii,jj,kk) + &
gxz(ig,jg,kg)*Gamxx(ii,jj,kk) + gyz(ig,jg,kg)*Gamyx(ii,jj,kk) + val_gzz*Gamzx(ii,jj,kk) + &
Gamxa*gxzx(ii,jj,kk) + Gamya*gyzx(ii,jj,kk) + Gamza*gzzx(ii,jj,kk) + &
Gamxa*gxxz(ii,jj,kk) + Gamya*gxyz(ii,jj,kk) + Gamza*gxzz(ii,jj,kk) ) + &
gupxx(ii,jj,kk)*(fxx(ii,jj,kk) + Gamxxx(ig,jg,kg)*gxxz(ii,jj,kk) + Gamyxx(ig,jg,kg)*gxyz(ii,jj,kk) + Gamzxx(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamxxx(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gxzz(ii,jj,kk)) + &
gupxy(ii,jj,kk)*(fxy(ii,jj,kk) + Gamxxx(ig,jg,kg)*gxyz(ii,jj,kk) + Gamyxx(ig,jg,kg)*gyyz(ii,jj,kk) + Gamzxx(ig,jg,kg)*gyzz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyyx(ii,jj,kk) + Gamzxz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxxz(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxyz(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamzyz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamxxx(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gyzz(ii,jj,kk)) + &
gupxz(ii,jj,kk)*(fxz(ii,jj,kk) + Gamxxx(ig,jg,kg)*gxzz(ii,jj,kk) + Gamyxx(ig,jg,kg)*gyzz(ii,jj,kk) + Gamzxx(ig,jg,kg)*gzzz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamzxz(ig,jg,kg)*gzzx(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxxz(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxzz(ig,jg,kg)*gxxx(ii,jj,kk) + Gamyzz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamzzz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamxxx(ig,jg,kg)*gzzx(ii,jj,kk) + Gamyxx(ig,jg,kg)*gzzy(ii,jj,kk) + Gamzxx(ig,jg,kg)*gzzz(ii,jj,kk)) + &
gupyy(ii,jj,kk)*(Gamxxy(ig,jg,kg)*gxyz(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyyz(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyzz(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyyx(ii,jj,kk) + Gamzyz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamxxy(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyzz(ii,jj,kk)) + &
gupyz(ii,jj,kk)*(fyz(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxzz(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyzz(ii,jj,kk) + Gamzxy(ig,jg,kg)*gzzz(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamzyz(ig,jg,kg)*gzzx(ii,jj,kk) + Gamxxz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gyzz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyyz(ii,jj,kk) + Gamzxz(ig,jg,kg)*gyzz(ii,jj,kk) + Gamxzz(ig,jg,kg)*gxyx(ii,jj,kk) + Gamyzz(ig,jg,kg)*gyyx(ii,jj,kk) + Gamzzz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamxxy(ig,jg,kg)*gzzx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gzzy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gzzz(ii,jj,kk)) + &
gupzz(ii,jj,kk)*(Gamxxz(ig,jg,kg)*gxzz(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyzz(ii,jj,kk) + Gamzxz(ig,jg,kg)*gzzz(ii,jj,kk) + Gamxzz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyzz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamzzz(ig,jg,kg)*gzzx(ii,jj,kk) + Gamxxz(ig,jg,kg)*gzzx(ii,jj,kk) + Gamyxz(ig,jg,kg)*gzzy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gzzz(ii,jj,kk))
end do; end do; end do
call calc_dderivs(gyz, fxx, fxy, fxz, fyy, fyz, fzz, SYM, ANTI, ANTI, is, ie, js, je, ks, ke)
!DIR$ SIMD
do kk=1,nz; do jj=1,ny; do ii=1,nx
ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
! Compute Gam^i at local point
Gamxa = Gamx(ig,jg,kg)
Gamya = Gamy(ig,jg,kg)
Gamza = Gamz(ig,jg,kg)
val_gyy = dyy(ig,jg,kg)+ONE
val_gzz = dzz(ig,jg,kg)+ONE
Ryz_l(ii,jj,kk) = HALF*( -Ryz_l(ii,jj,kk) + gxy(ig,jg,kg)*Gamxz(ii,jj,kk) + val_gyy*Gamyz(ii,jj,kk) + gyz(ig,jg,kg)*Gamzz(ii,jj,kk) + &
gxz(ig,jg,kg)*Gamxy(ii,jj,kk) + gyz(ig,jg,kg)*Gamyy(ii,jj,kk) + val_gzz*Gamzy(ii,jj,kk) + &
Gamxa*gxzy(ii,jj,kk) + Gamya*gyzy(ii,jj,kk) + Gamza*gzzy(ii,jj,kk) + &
Gamxa*gxyz(ii,jj,kk) + Gamya*gyyz(ii,jj,kk) + Gamza*gyzz(ii,jj,kk) ) + &
gupxx(ii,jj,kk)*(Gamxxy(ig,jg,kg)*gxxz(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxyz(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxxy(ii,jj,kk) + Gamyxz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gxzz(ii,jj,kk)) + &
gupxy(ii,jj,kk)*(fxy(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxyz(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyyz(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyzz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamxyy(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyyy(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzyy(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxyy(ig,jg,kg)*gxxz(ii,jj,kk) + Gamyyy(ig,jg,kg)*gxyz(ii,jj,kk) + Gamzyy(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxxy(ii,jj,kk) + Gamyyz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamxxy(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gyzz(ii,jj,kk)) + &
gupxz(ii,jj,kk)*(fxz(ii,jj,kk) + Gamxxy(ig,jg,kg)*gxzz(ii,jj,kk) + Gamyxy(ig,jg,kg)*gyzz(ii,jj,kk) + Gamzxy(ig,jg,kg)*gzzz(ii,jj,kk) + Gamxxz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamyxz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzxz(ig,jg,kg)*gzzy(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxzx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxxz(ii,jj,kk) + Gamyyz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamzyz(ig,jg,kg)*gxzz(ii,jj,kk) + Gamxzz(ig,jg,kg)*gxxy(ii,jj,kk) + Gamyzz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamzzz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamxxy(ig,jg,kg)*gzzx(ii,jj,kk) + Gamyxy(ig,jg,kg)*gzzy(ii,jj,kk) + Gamzxy(ig,jg,kg)*gzzz(ii,jj,kk)) + &
gupyy(ii,jj,kk)*(fyy(ii,jj,kk) + Gamxyy(ig,jg,kg)*gxyz(ii,jj,kk) + Gamyyy(ig,jg,kg)*gyyz(ii,jj,kk) + Gamzyy(ig,jg,kg)*gyzz(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamxyy(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyyy(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzyy(ig,jg,kg)*gyzz(ii,jj,kk)) + &
gupyz(ii,jj,kk)*(fyz(ii,jj,kk) + Gamxyy(ig,jg,kg)*gxzz(ii,jj,kk) + Gamyyy(ig,jg,kg)*gyzz(ii,jj,kk) + Gamzyy(ig,jg,kg)*gzzz(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gzzy(ii,jj,kk) + Gamxyz(ig,jg,kg)*gyzx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gyzz(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxyz(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyyz(ii,jj,kk) + Gamzyz(ig,jg,kg)*gyzz(ii,jj,kk) + Gamxzz(ig,jg,kg)*gxyy(ii,jj,kk) + Gamyzz(ig,jg,kg)*gyyy(ii,jj,kk) + Gamzzz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamxyy(ig,jg,kg)*gzzx(ii,jj,kk) + Gamyyy(ig,jg,kg)*gzzy(ii,jj,kk) + Gamzyy(ig,jg,kg)*gzzz(ii,jj,kk)) + &
gupzz(ii,jj,kk)*(fzz(ii,jj,kk) + Gamxyz(ig,jg,kg)*gxzz(ii,jj,kk) + Gamyyz(ig,jg,kg)*gyzz(ii,jj,kk) + Gamzyz(ig,jg,kg)*gzzz(ii,jj,kk) + Gamxzz(ig,jg,kg)*gxzy(ii,jj,kk) + Gamyzz(ig,jg,kg)*gyzy(ii,jj,kk) + Gamzzz(ig,jg,kg)*gzzy(ii,jj,kk) + Gamxyz(ig,jg,kg)*gzzx(ii,jj,kk) + Gamyyz(ig,jg,kg)*gzzy(ii,jj,kk) + Gamzyz(ig,jg,kg)*gzzz(ii,jj,kk))
end do; end do; end do
! === 7. Chi Terms for Ricci ===
call calc_dderivs(chi, fxx, fxy, fxz, fyy, fyz, fzz, SYM, SYM, SYM, is, ie, js, je, ks, ke)
!DIR$ SIMD
do kk=1,nz; do jj=1,ny; do ii=1,nx
ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
chin1 = chi(ig,jg,kg) + ONE
! Compute trace part f (D^2 chi terms)
f(ii,jj,kk) = gupxx(ii,jj,kk)*(fxx(ii,jj,kk) - F3o2/chin1*chix(ii,jj,kk)**2) + &
gupyy(ii,jj,kk)*(fyy(ii,jj,kk) - F3o2/chin1*chiy(ii,jj,kk)**2) + &
gupzz(ii,jj,kk)*(fzz(ii,jj,kk) - F3o2/chin1*chiz(ii,jj,kk)**2) + &
TWO*(gupxy(ii,jj,kk)*(fxy(ii,jj,kk) - F3o2/chin1*chix(ii,jj,kk)*chiy(ii,jj,kk)) + &
gupxz(ii,jj,kk)*(fxz(ii,jj,kk) - F3o2/chin1*chix(ii,jj,kk)*chiz(ii,jj,kk)) + &
gupyz(ii,jj,kk)*(fyz(ii,jj,kk) - F3o2/chin1*chiy(ii,jj,kk)*chiz(ii,jj,kk)))
! Add chi terms to Ricci
val_gxx = dxx(ig,jg,kg)+ONE; val_gyy = dyy(ig,jg,kg)+ONE; val_gzz = dzz(ig,jg,kg)+ONE
Rxx_l(ii,jj,kk) = Rxx_l(ii,jj,kk) + (fxx(ii,jj,kk) - chix(ii,jj,kk)**2/chin1/TWO + val_gxx*f(ii,jj,kk))/chin1/TWO
Ryy_l(ii,jj,kk) = Ryy_l(ii,jj,kk) + (fyy(ii,jj,kk) - chiy(ii,jj,kk)**2/chin1/TWO + val_gyy*f(ii,jj,kk))/chin1/TWO
Rzz_l(ii,jj,kk) = Rzz_l(ii,jj,kk) + (fzz(ii,jj,kk) - chiz(ii,jj,kk)**2/chin1/TWO + val_gzz*f(ii,jj,kk))/chin1/TWO
Rxy_l(ii,jj,kk) = Rxy_l(ii,jj,kk) + (fxy(ii,jj,kk) - chix(ii,jj,kk)*chiy(ii,jj,kk)/chin1/TWO + gxy(ig,jg,kg)*f(ii,jj,kk))/chin1/TWO
Rxz_l(ii,jj,kk) = Rxz_l(ii,jj,kk) + (fxz(ii,jj,kk) - chix(ii,jj,kk)*chiz(ii,jj,kk)/chin1/TWO + gxz(ig,jg,kg)*f(ii,jj,kk))/chin1/TWO
Ryz_l(ii,jj,kk) = Ryz_l(ii,jj,kk) + (fyz(ii,jj,kk) - chiy(ii,jj,kk)*chiz(ii,jj,kk)/chin1/TWO + gyz(ig,jg,kg)*f(ii,jj,kk))/chin1/TWO
! Also need f for A_ij RHS later? No, f is recalculated from Lap later.
end do; end do; end do
! === 8. Lap Terms ===
call calc_dderivs(Lap, fxx, fxy, fxz, fyy, fyz, fzz, SYM, SYM, SYM, is, ie, js, je, ks, ke)
!DIR$ SIMD
do kk=1,nz; do jj=1,ny; do ii=1,nx
ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
! Subtract connection terms from D_i D_j Lap (fxx is just d_i d_j Lap initially)
fxx(ii,jj,kk) = fxx(ii,jj,kk) - Gamxxx(ig,jg,kg)*Lapx(ii,jj,kk) - Gamyxx(ig,jg,kg)*Lapy(ii,jj,kk) - Gamzxx(ig,jg,kg)*Lapz(ii,jj,kk)
fyy(ii,jj,kk) = fyy(ii,jj,kk) - Gamxyy(ig,jg,kg)*Lapx(ii,jj,kk) - Gamyyy(ig,jg,kg)*Lapy(ii,jj,kk) - Gamzyy(ig,jg,kg)*Lapz(ii,jj,kk)
fzz(ii,jj,kk) = fzz(ii,jj,kk) - Gamxzz(ig,jg,kg)*Lapx(ii,jj,kk) - Gamyzz(ig,jg,kg)*Lapy(ii,jj,kk) - Gamzzz(ig,jg,kg)*Lapz(ii,jj,kk)
fxy(ii,jj,kk) = fxy(ii,jj,kk) - Gamxxy(ig,jg,kg)*Lapx(ii,jj,kk) - Gamyxy(ig,jg,kg)*Lapy(ii,jj,kk) - Gamzxy(ig,jg,kg)*Lapz(ii,jj,kk)
fxz(ii,jj,kk) = fxz(ii,jj,kk) - Gamxxz(ig,jg,kg)*Lapx(ii,jj,kk) - Gamyxz(ig,jg,kg)*Lapy(ii,jj,kk) - Gamzxz(ig,jg,kg)*Lapz(ii,jj,kk)
fyz(ii,jj,kk) = fyz(ii,jj,kk) - Gamxyz(ig,jg,kg)*Lapx(ii,jj,kk) - Gamyyz(ig,jg,kg)*Lapy(ii,jj,kk) - Gamzyz(ig,jg,kg)*Lapz(ii,jj,kk)
! trK_rhs partial update: g^ij D_i D_j Lap
trK_rhs(ig,jg,kg) = chin1 * (gupxx(ii,jj,kk)*fxx(ii,jj,kk) + gupyy(ii,jj,kk)*fyy(ii,jj,kk) + gupzz(ii,jj,kk)*fzz(ii,jj,kk) + &
TWO*(gupxy(ii,jj,kk)*fxy(ii,jj,kk) + gupxz(ii,jj,kk)*fxz(ii,jj,kk) + gupyz(ii,jj,kk)*fyz(ii,jj,kk)))
! Prepare fxx etc for A_ij RHS (Add R_ij and S_ij parts)
! Need Lap value
alpn1 = Lap(ig,jg,kg) + ONE
! Store fxx_new = alpn1 * (Rxx - 8PI Sxx) - D_i D_j Lap
fxx(ii,jj,kk) = alpn1 * (Rxx_l(ii,jj,kk) - EIGHT*PI*Sxx(ig,jg,kg)) - fxx(ii,jj,kk)
fyy(ii,jj,kk) = alpn1 * (Ryy_l(ii,jj,kk) - EIGHT*PI*Syy(ig,jg,kg)) - fyy(ii,jj,kk)
fzz(ii,jj,kk) = alpn1 * (Rzz_l(ii,jj,kk) - EIGHT*PI*Szz(ig,jg,kg)) - fzz(ii,jj,kk)
fxy(ii,jj,kk) = alpn1 * (Rxy_l(ii,jj,kk) - EIGHT*PI*Sxy(ig,jg,kg)) - fxy(ii,jj,kk)
fxz(ii,jj,kk) = alpn1 * (Rxz_l(ii,jj,kk) - EIGHT*PI*Sxz(ig,jg,kg)) - fxz(ii,jj,kk)
fyz(ii,jj,kk) = alpn1 * (Ryz_l(ii,jj,kk) - EIGHT*PI*Syz(ig,jg,kg)) - fyz(ii,jj,kk)
end do; end do; end do
! === 9. A_ij RHS and Final trK RHS ===
!DIR$ SIMD
do kk=1,nz; do jj=1,ny; do ii=1,nx
ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
val_gxx = dxx(ig,jg,kg)+ONE; val_gyy = dyy(ig,jg,kg)+ONE; val_gzz = dzz(ig,jg,kg)+ONE
alpn1 = Lap(ig,jg,kg) + ONE
chin1 = chi(ig,jg,kg) + ONE
div_beta = betaxx(ii,jj,kk) + betayy(ii,jj,kk) + betazz(ii,jj,kk)
! Trace-free part f
f(ii,jj,kk) = F1o3 * (gupxx(ii,jj,kk)*fxx(ii,jj,kk) + gupyy(ii,jj,kk)*fyy(ii,jj,kk) + gupzz(ii,jj,kk)*fzz(ii,jj,kk) + &
TWO*(gupxy(ii,jj,kk)*fxy(ii,jj,kk) + gupxz(ii,jj,kk)*fxz(ii,jj,kk) + gupyz(ii,jj,kk)*fyz(ii,jj,kk)))
! Initial A_ij RHS (TF(R_ij part))
Axx_rhs(ig,jg,kg) = fxx(ii,jj,kk) - val_gxx * f(ii,jj,kk)
Ayy_rhs(ig,jg,kg) = fyy(ii,jj,kk) - val_gyy * f(ii,jj,kk)
Azz_rhs(ig,jg,kg) = fzz(ii,jj,kk) - val_gzz * f(ii,jj,kk)
Axy_rhs(ig,jg,kg) = fxy(ii,jj,kk) - gxy(ig,jg,kg) * f(ii,jj,kk)
Axz_rhs(ig,jg,kg) = fxz(ii,jj,kk) - gxz(ig,jg,kg) * f(ii,jj,kk)
Ayz_rhs(ig,jg,kg) = fyz(ii,jj,kk) - gyz(ig,jg,kg) * f(ii,jj,kk)
! Add A_il A^l_j terms (Calculate fxx as A_ik A^k_j)
! fxx = A_xx^2 ...
! Re-calculate fxx... reusing local arrays
! fxx = g^lm A_xl A_mj
fxx(ii,jj,kk) = gupxx(ii,jj,kk)*Axx(ig,jg,kg)**2 + gupyy(ii,jj,kk)*Axy(ig,jg,kg)**2 + gupzz(ii,jj,kk)*Axz(ig,jg,kg)**2 + &
TWO*(gupxy(ii,jj,kk)*Axx(ig,jg,kg)*Axy(ig,jg,kg) + gupxz(ii,jj,kk)*Axx(ig,jg,kg)*Axz(ig,jg,kg) + gupyz(ii,jj,kk)*Axy(ig,jg,kg)*Axz(ig,jg,kg))
fyy(ii,jj,kk) = gupxx(ii,jj,kk)*Axy(ig,jg,kg)**2 + gupyy(ii,jj,kk)*Ayy(ig,jg,kg)**2 + gupzz(ii,jj,kk)*Ayz(ig,jg,kg)**2 + &
TWO*(gupxy(ii,jj,kk)*Axy(ig,jg,kg)*Ayy(ig,jg,kg) + gupxz(ii,jj,kk)*Axy(ig,jg,kg)*Ayz(ig,jg,kg) + gupyz(ii,jj,kk)*Ayy(ig,jg,kg)*Ayz(ig,jg,kg))
fzz(ii,jj,kk) = gupxx(ii,jj,kk)*Axz(ig,jg,kg)**2 + gupyy(ii,jj,kk)*Ayz(ig,jg,kg)**2 + gupzz(ii,jj,kk)*Azz(ig,jg,kg)**2 + &
TWO*(gupxy(ii,jj,kk)*Axz(ig,jg,kg)*Ayz(ig,jg,kg) + gupxz(ii,jj,kk)*Axz(ig,jg,kg)*Azz(ig,jg,kg) + gupyz(ii,jj,kk)*Ayz(ig,jg,kg)*Azz(ig,jg,kg))
fxy(ii,jj,kk) = gupxx(ii,jj,kk)*Axx(ig,jg,kg)*Axy(ig,jg,kg) + gupyy(ii,jj,kk)*Axy(ig,jg,kg)*Ayy(ig,jg,kg) + gupzz(ii,jj,kk)*Axz(ig,jg,kg)*Ayz(ig,jg,kg) + &
gupxy(ii,jj,kk)*(Axx(ig,jg,kg)*Ayy(ig,jg,kg) + Axy(ig,jg,kg)**2) + &
gupxz(ii,jj,kk)*(Axx(ig,jg,kg)*Ayz(ig,jg,kg) + Axz(ig,jg,kg)*Axy(ig,jg,kg)) + &
gupyz(ii,jj,kk)*(Axy(ig,jg,kg)*Ayz(ig,jg,kg) + Axz(ig,jg,kg)*Ayy(ig,jg,kg))
fxz(ii,jj,kk) = gupxx(ii,jj,kk)*Axx(ig,jg,kg)*Axz(ig,jg,kg) + gupyy(ii,jj,kk)*Axy(ig,jg,kg)*Ayz(ig,jg,kg) + gupzz(ii,jj,kk)*Axz(ig,jg,kg)*Azz(ig,jg,kg) + &
gupxy(ii,jj,kk)*(Axx(ig,jg,kg)*Ayz(ig,jg,kg) + Axy(ig,jg,kg)*Axz(ig,jg,kg)) + &
gupxz(ii,jj,kk)*(Axx(ig,jg,kg)*Azz(ig,jg,kg) + Axz(ig,jg,kg)**2) + &
gupyz(ii,jj,kk)*(Axy(ig,jg,kg)*Azz(ig,jg,kg) + Axz(ig,jg,kg)*Ayz(ig,jg,kg))
fyz(ii,jj,kk) = gupxx(ii,jj,kk)*Axy(ig,jg,kg)*Axz(ig,jg,kg) + gupyy(ii,jj,kk)*Ayy(ig,jg,kg)*Ayz(ig,jg,kg) + gupzz(ii,jj,kk)*Ayz(ig,jg,kg)*Azz(ig,jg,kg) + &
gupxy(ii,jj,kk)*(Axy(ig,jg,kg)*Ayz(ig,jg,kg) + Ayy(ig,jg,kg)*Axz(ig,jg,kg)) + &
gupxz(ii,jj,kk)*(Axy(ig,jg,kg)*Azz(ig,jg,kg) + Ayz(ig,jg,kg)*Axz(ig,jg,kg)) + &
gupyz(ii,jj,kk)*(Ayy(ig,jg,kg)*Azz(ig,jg,kg) + Ayz(ig,jg,kg)**2)
! Update A_ij RHS with trK and A^2 terms and Lie derivative terms (shift part)
Axx_rhs(ig,jg,kg) = chin1*Axx_rhs(ig,jg,kg) + alpn1*(trK(ig,jg,kg)*Axx(ig,jg,kg) - TWO*fxx(ii,jj,kk)) + &
TWO*(Axx(ig,jg,kg)*betaxx(ii,jj,kk) + Axy(ig,jg,kg)*betayx(ii,jj,kk) + Axz(ig,jg,kg)*betazx(ii,jj,kk)) - F2o3*Axx(ig,jg,kg)*div_beta
Ayy_rhs(ig,jg,kg) = chin1*Ayy_rhs(ig,jg,kg) + alpn1*(trK(ig,jg,kg)*Ayy(ig,jg,kg) - TWO*fyy(ii,jj,kk)) + &
TWO*(Axy(ig,jg,kg)*betaxy(ii,jj,kk) + Ayy(ig,jg,kg)*betayy(ii,jj,kk) + Ayz(ig,jg,kg)*betazy(ii,jj,kk)) - F2o3*Ayy(ig,jg,kg)*div_beta
Azz_rhs(ig,jg,kg) = chin1*Azz_rhs(ig,jg,kg) + alpn1*(trK(ig,jg,kg)*Azz(ig,jg,kg) - TWO*fzz(ii,jj,kk)) + &
TWO*(Axz(ig,jg,kg)*betaxz(ii,jj,kk) + Ayz(ig,jg,kg)*betayz(ii,jj,kk) + Azz(ig,jg,kg)*betazz(ii,jj,kk)) - F2o3*Azz(ig,jg,kg)*div_beta
Axy_rhs(ig,jg,kg) = chin1*Axy_rhs(ig,jg,kg) + alpn1*(trK(ig,jg,kg)*Axy(ig,jg,kg) - TWO*fxy(ii,jj,kk)) + &
Axx(ig,jg,kg)*betaxy(ii,jj,kk) + Axz(ig,jg,kg)*betazy(ii,jj,kk) + &
Ayy(ig,jg,kg)*betayx(ii,jj,kk) + Ayz(ig,jg,kg)*betazx(ii,jj,kk) + &
F1o3*Axy(ig,jg,kg)*div_beta - Axy(ig,jg,kg)*betazz(ii,jj,kk)
Axz_rhs(ig,jg,kg) = chin1*Axz_rhs(ig,jg,kg) + alpn1*(trK(ig,jg,kg)*Axz(ig,jg,kg) - TWO*fxz(ii,jj,kk)) + &
Axx(ig,jg,kg)*betaxz(ii,jj,kk) + Axy(ig,jg,kg)*betayz(ii,jj,kk) + &
Ayz(ig,jg,kg)*betayx(ii,jj,kk) + Azz(ig,jg,kg)*betazx(ii,jj,kk) + &
F1o3*Axz(ig,jg,kg)*div_beta - Axz(ig,jg,kg)*betayy(ii,jj,kk)
Ayz_rhs(ig,jg,kg) = chin1*Ayz_rhs(ig,jg,kg) + alpn1*(trK(ig,jg,kg)*Ayz(ig,jg,kg) - TWO*fyz(ii,jj,kk)) + &
Axy(ig,jg,kg)*betaxz(ii,jj,kk) + Ayy(ig,jg,kg)*betayz(ii,jj,kk) + &
Axz(ig,jg,kg)*betaxy(ii,jj,kk) + Azz(ig,jg,kg)*betazy(ii,jj,kk) + &
F1o3*Ayz(ig,jg,kg)*div_beta - Ayz(ig,jg,kg)*betaxx(ii,jj,kk)
! Compute trace of S_ij (f is chin1 here)
S_val = chin1 * (gupxx(ii,jj,kk)*Sxx(ig,jg,kg) + gupyy(ii,jj,kk)*Syy(ig,jg,kg) + gupzz(ii,jj,kk)*Szz(ig,jg,kg) + &
TWO*(gupxy(ii,jj,kk)*Sxy(ig,jg,kg) + gupxz(ii,jj,kk)*Sxz(ig,jg,kg) + gupyz(ii,jj,kk)*Syz(ig,jg,kg)))
! Final trK RHS (Reuse fxx as A^2 trace)
! fxx (trace of A^2) = g^ij g^kl A_ik A_jl = f (computed earlier)
! Wait, I computed fxx...fzz as A^2 components.
! Trace is g^ij f_ij.
f(ii,jj,kk) = gupxx(ii,jj,kk)*fxx(ii,jj,kk) + gupyy(ii,jj,kk)*fyy(ii,jj,kk) + gupzz(ii,jj,kk)*fzz(ii,jj,kk) + &
TWO*(gupxy(ii,jj,kk)*fxy(ii,jj,kk) + gupxz(ii,jj,kk)*fxz(ii,jj,kk) + gupyz(ii,jj,kk)*fyz(ii,jj,kk))
trK_rhs(ig,jg,kg) = -trK_rhs(ig,jg,kg) + alpn1*(F1o3*trK(ig,jg,kg)**2 + f(ii,jj,kk) + FOUR*PI*(rho(ig,jg,kg)+S_val))
end do; end do; end do
! === 10. Gauge Conditions and Output Copy ===
!DIR$ SIMD
do kk=1,nz; do jj=1,ny; do ii=1,nx
ig=is+ii-1; jg=js+jj-1; kg=ks+kk-1
alpn1 = Lap(ig,jg,kg) + ONE
chin1 = chi(ig,jg,kg) + ONE
Lap_rhs(ig,jg,kg) = -TWO*alpn1*trK(ig,jg,kg)
! Betax RHS (GAUGE == 2)
betax_rhs(ig,jg,kg) = FF*dtSfx(ig,jg,kg)
betay_rhs(ig,jg,kg) = FF*dtSfy(ig,jg,kg)
betaz_rhs(ig,jg,kg) = FF*dtSfz(ig,jg,kg)
! dtSfx RHS (GAUGE == 2)
! reta calculation using chix... computed in Step 1
reta_val = gupxx(ii,jj,kk)*chix(ii,jj,kk)**2 + gupyy(ii,jj,kk)*chiy(ii,jj,kk)**2 + gupzz(ii,jj,kk)*chiz(ii,jj,kk)**2 + &
TWO*(gupxy(ii,jj,kk)*chix(ii,jj,kk)*chiy(ii,jj,kk) + gupxz(ii,jj,kk)*chix(ii,jj,kk)*chiz(ii,jj,kk) + gupyz(ii,jj,kk)*chiy(ii,jj,kk)*chiz(ii,jj,kk))
if (chin1 > 1.0d-10) then
denom = ONE - sqrt(chin1)
if (abs(denom) < 1.0d-15) denom = 1.0d-15
reta_val = 1.31d0/TWO * sqrt(max(reta_val, ZEO)/chin1) / denom**2
else
reta_val = ZEO
endif
dtSfx_rhs(ig,jg,kg) = Gamx_rhs(ig,jg,kg) - reta_val*dtSfx(ig,jg,kg)
dtSfy_rhs(ig,jg,kg) = Gamy_rhs(ig,jg,kg) - reta_val*dtSfy(ig,jg,kg)
dtSfz_rhs(ig,jg,kg) = Gamz_rhs(ig,jg,kg) - reta_val*dtSfz(ig,jg,kg)
! Copy Rxx_l back to global Rxx (required by interface)
Rxx(ig,jg,kg) = Rxx_l(ii,jj,kk)
Ryy(ig,jg,kg) = Ryy_l(ii,jj,kk)
Rzz(ig,jg,kg) = Rzz_l(ii,jj,kk)
Rxy(ig,jg,kg) = Rxy_l(ii,jj,kk)
Rxz(ig,jg,kg) = Rxz_l(ii,jj,kk)
Ryz(ig,jg,kg) = Ryz_l(ii,jj,kk)
end do; end do; end do
call apply_lopsided(dtSfx, dtSfx_rhs, betax, betay, betaz, is, ie, js, je, ks, ke)
call apply_lopsided(dtSfy, dtSfy_rhs, betax, betay, betaz, is, ie, js, je, ks, ke)
call apply_lopsided(dtSfz, dtSfz_rhs, betax, betay, betaz, is, ie, js, je, ks, ke)
call apply_lopsided(betax, betax_rhs, betax, betay, betaz, is, ie, js, je, ks, ke)
call apply_lopsided(betay, betay_rhs, betax, betay, betaz, is, ie, js, je, ks, ke)
call apply_lopsided(betaz, betaz_rhs, betax, betay, betaz, is, ie, js, je, ks, ke)
! Apply dissipation to gauge vars
if (eps > 0.0d0) then
call apply_kodiss(dtSfx, dtSfx_rhs, is, ie, js, je, ks, ke)
call apply_kodiss(dtSfy, dtSfy_rhs, is, ie, js, je, ks, ke)
call apply_kodiss(dtSfz, dtSfz_rhs, is, ie, js, je, ks, ke)
call apply_kodiss(betax, betax_rhs, is, ie, js, je, ks, ke)
call apply_kodiss(betay, betay_rhs, is, ie, js, je, ks, ke)
call apply_kodiss(betaz, betaz_rhs, is, ie, js, je, ks, ke)
end if
end subroutine compute_block_kernel
!-------------------------------------------------------------------------
! Helper: Safe 1st Derivative (4th Order with 2nd Order Fallback)
!-------------------------------------------------------------------------
subroutine calc_derivs(f, fx, fy, fz, sym1, sym2, sym3, is, ie, js, je, ks, ke)
real*8, dimension(ex(1),ex(2),ex(3)), intent(in) :: f
real*8, dimension(BLK,BLK,BLK), intent(out) :: fx, fy, fz
real*8, intent(in) :: sym1, sym2, sym3
integer, intent(in) :: is, ie, js, je, ks, ke
integer :: ii, jj, kk, ig, jg, kg
real*8 :: d12dx, d12dy, d12dz, d2dx, d2dy, d2dz
real*8 :: vM2, vM1, vP1, vP2
! Boundary Handling
integer, parameter :: NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2
integer :: imin, jmin, kmin, imax, jmax, kmax
real*8 :: dX_l, dY_l, dZ_l
dX_l = X(2)-X(1)
dY_l = Y(2)-Y(1)
dZ_l = Z(2)-Z(1)
d12dx = ONE / (12.d0 * dX_l)
d12dy = ONE / (12.d0 * dY_l)
d12dz = ONE / (12.d0 * dZ_l)
d2dx = ONE / (TWO * dX_l)
d2dy = ONE / (TWO * dY_l)
d2dz = ONE / (TWO * dZ_l)
imax = ex(1); jmax = ex(2); kmax = ex(3)
imin = 1; jmin = 1; kmin = 1
if(Symmetry > NO_SYMM .and. dabs(Z(1)) < dZ_l) kmin = -1
if(Symmetry > EQ_SYMM .and. dabs(X(1)) < dX_l) imin = -1
if(Symmetry > EQ_SYMM .and. dabs(Y(1)) < dY_l) jmin = -1
! Initialize to zero (for points where derivative cannot be computed)
fx = ZEO; fy = ZEO; fz = ZEO
!DIR$ SIMD
do kk=1,ke-ks+1; do jj=1,je-js+1; do ii=1,ie-is+1
ig = is + ii - 1; jg = js + jj - 1; kg = ks + kk - 1
! X-Derivative
if (ig+2 <= imax .and. ig-2 >= imin) then
! 4th Order
if (ig > 2 .and. ig < ex(1)-1) then
fx(ii,jj,kk) = d12dx * (f(ig-2,jg,kg) - EIGHT*f(ig-1,jg,kg) + EIGHT*f(ig+1,jg,kg) - f(ig+2,jg,kg))
else
vM2 = get_val(f, ig-2, jg, kg, sym1, 1); vM1 = get_val(f, ig-1, jg, kg, sym1, 1)
vP1 = get_val(f, ig+1, jg, kg, sym1, 1); vP2 = get_val(f, ig+2, jg, kg, sym1, 1)
fx(ii,jj,kk) = d12dx * (vM2 - EIGHT*vM1 + EIGHT*vP1 - vP2)
endif
elseif (ig+1 <= imax .and. ig-1 >= imin) then
! 2nd Order Fallback
if (ig > 1 .and. ig < ex(1)) then
fx(ii,jj,kk) = d2dx * (-f(ig-1,jg,kg) + f(ig+1,jg,kg))
else
vM1 = get_val(f, ig-1, jg, kg, sym1, 1)
vP1 = get_val(f, ig+1, jg, kg, sym1, 1)
fx(ii,jj,kk) = d2dx * (-vM1 + vP1)
endif
endif
! Y-Derivative
if (jg+2 <= jmax .and. jg-2 >= jmin) then
if (jg > 2 .and. jg < ex(2)-1) then
fy(ii,jj,kk) = d12dy * (f(ig,jg-2,kg) - EIGHT*f(ig,jg-1,kg) + EIGHT*f(ig,jg+1,kg) - f(ig,jg+2,kg))
else
vM2 = get_val(f, ig, jg-2, kg, sym2, 2); vM1 = get_val(f, ig, jg-1, kg, sym2, 2)
vP1 = get_val(f, ig, jg+1, kg, sym2, 2); vP2 = get_val(f, ig, jg+2, kg, sym2, 2)
fy(ii,jj,kk) = d12dy * (vM2 - EIGHT*vM1 + EIGHT*vP1 - vP2)
endif
elseif (jg+1 <= jmax .and. jg-1 >= jmin) then
if (jg > 1 .and. jg < ex(2)) then
fy(ii,jj,kk) = d2dy * (-f(ig,jg-1,kg) + f(ig,jg+1,kg))
else
vM1 = get_val(f, ig, jg-1, kg, sym2, 2)
vP1 = get_val(f, ig, jg+1, kg, sym2, 2)
fy(ii,jj,kk) = d2dy * (-vM1 + vP1)
endif
endif
! Z-Derivative
if (kg+2 <= kmax .and. kg-2 >= kmin) then
if (kg > 2 .and. kg < ex(3)-1) then
fz(ii,jj,kk) = d12dz * (f(ig,jg,kg-2) - EIGHT*f(ig,jg,kg-1) + EIGHT*f(ig,jg,kg+1) - f(ig,jg,kg+2))
else
vM2 = get_val(f, ig, jg, kg-2, sym3, 3); vM1 = get_val(f, ig, jg, kg-1, sym3, 3)
vP1 = get_val(f, ig, jg, kg+1, sym3, 3); vP2 = get_val(f, ig, jg, kg+2, sym3, 3)
fz(ii,jj,kk) = d12dz * (vM2 - EIGHT*vM1 + EIGHT*vP1 - vP2)
endif
elseif (kg+1 <= kmax .and. kg-1 >= kmin) then
if (kg > 1 .and. kg < ex(3)) then
fz(ii,jj,kk) = d2dz * (-f(ig,jg,kg-1) + f(ig,jg,kg+1))
else
vM1 = get_val(f, ig, jg, kg-1, sym3, 3)
vP1 = get_val(f, ig, jg, kg+1, sym3, 3)
fz(ii,jj,kk) = d2dz * (-vM1 + vP1)
endif
endif
end do; end do; end do
end subroutine calc_derivs
!-------------------------------------------------------------------------
! Helper: Safe 2nd Derivative (Mixed & Diag)
!-------------------------------------------------------------------------
subroutine calc_dderivs(f, fxx, fxy, fxz, fyy, fyz, fzz, sym1, sym2, sym3, is, ie, js, je, ks, ke)
real*8, dimension(ex(1),ex(2),ex(3)), intent(in) :: f
real*8, dimension(BLK,BLK,BLK), intent(out) :: fxx, fxy, fxz, fyy, fyz, fzz
real*8, intent(in) :: sym1, sym2, sym3
integer, intent(in) :: is, ie, js, je, ks, ke
integer :: ii, jj, kk, ig, jg, kg
real*8 :: Fdxdx, Fdydy, Fdzdz, Fdxdy, Fdxdz, Fdydz
real*8 :: Sdxdx, Sdydy, Sdzdz, Sdxdy, Sdxdz, Sdydz
real*8 :: vM2, vM1, vP1, vP2
real*8 :: vM1M1, vP1M1, vM1P1, vP1P1
! Boundary Handling
integer, parameter :: NO_SYMM = 0, EQ_SYMM = 1, OCTANT = 2
integer :: imin, jmin, kmin, imax, jmax, kmax
real*8 :: dX_l, dY_l, dZ_l
dX_l = X(2)-X(1)
dY_l = Y(2)-Y(1)
dZ_l = Z(2)-Z(1)
Fdxdx = ONE / (12.d0 * dX_l**2)
Fdydy = ONE / (12.d0 * dY_l**2)
Fdzdz = ONE / (12.d0 * dZ_l**2)
Fdxdy = ONE / (144.d0 * dX_l*dY_l)
Fdxdz = ONE / (144.d0 * dX_l*dZ_l)
Fdydz = ONE / (144.d0 * dY_l*dZ_l)
Sdxdx = ONE / (dX_l**2)
Sdydy = ONE / (dY_l**2)
Sdzdz = ONE / (dZ_l**2)
Sdxdy = F1o4 / (dX_l*dY_l)
Sdxdz = F1o4 / (dX_l*dZ_l)
Sdydz = F1o4 / (dY_l*dZ_l)
imax = ex(1); jmax = ex(2); kmax = ex(3)
imin = 1; jmin = 1; kmin = 1
if(Symmetry > NO_SYMM .and. dabs(Z(1)) < dZ_l) kmin = -1
if(Symmetry > EQ_SYMM .and. dabs(X(1)) < dX_l) imin = -1
if(Symmetry > EQ_SYMM .and. dabs(Y(1)) < dY_l) jmin = -1
! Initialize to zero
fxx=ZEO; fyy=ZEO; fzz=ZEO
fxy=ZEO; fxz=ZEO; fyz=ZEO
!DIR$ SIMD
do kk=1,ke-ks+1; do jj=1,je-js+1; do ii=1,ie-is+1
ig = is + ii - 1; jg = js + jj - 1; kg = ks + kk - 1
! FXX
if (ig+2 <= imax .and. ig-2 >= imin) then
if (ig > 2 .and. ig < ex(1)-1) then
fxx(ii,jj,kk) = Fdxdx * (-f(ig-2,jg,kg) + F16*f(ig-1,jg,kg) - F30*f(ig,jg,kg) + F16*f(ig+1,jg,kg) - f(ig+2,jg,kg))
else
vM2 = get_val(f, ig-2, jg, kg, sym1, 1); vM1 = get_val(f, ig-1, jg, kg, sym1, 1)
vP1 = get_val(f, ig+1, jg, kg, sym1, 1); vP2 = get_val(f, ig+2, jg, kg, sym1, 1)
fxx(ii,jj,kk) = Fdxdx * (-vM2 + F16*vM1 - F30*f(ig,jg,kg) + F16*vP1 - vP2)
endif
elseif (ig+1 <= imax .and. ig-1 >= imin) then
if (ig > 1 .and. ig < ex(1)) then
fxx(ii,jj,kk) = Sdxdx * (f(ig-1,jg,kg) - TWO*f(ig,jg,kg) + f(ig+1,jg,kg))
else
vM1 = get_val(f, ig-1, jg, kg, sym1, 1); vP1 = get_val(f, ig+1, jg, kg, sym1, 1)
fxx(ii,jj,kk) = Sdxdx * (vM1 - TWO*f(ig,jg,kg) + vP1)
endif
endif
! FYY
if (jg+2 <= jmax .and. jg-2 >= jmin) then
if (jg > 2 .and. jg < ex(2)-1) then
fyy(ii,jj,kk) = Fdydy * (-f(ig,jg-2,kg) + F16*f(ig,jg-1,kg) - F30*f(ig,jg,kg) + F16*f(ig,jg+1,kg) - f(ig,jg+2,kg))
else
vM2 = get_val(f, ig, jg-2, kg, sym2, 2); vM1 = get_val(f, ig, jg-1, kg, sym2, 2)
vP1 = get_val(f, ig, jg+1, kg, sym2, 2); vP2 = get_val(f, ig, jg+2, kg, sym2, 2)
fyy(ii,jj,kk) = Fdydy * (-vM2 + F16*vM1 - F30*f(ig,jg,kg) + F16*vP1 - vP2)
endif
elseif (jg+1 <= jmax .and. jg-1 >= jmin) then
if (jg > 1 .and. jg < ex(2)) then
fyy(ii,jj,kk) = Sdydy * (f(ig,jg-1,kg) - TWO*f(ig,jg,kg) + f(ig,jg+1,kg))
else
vM1 = get_val(f, ig, jg-1, kg, sym2, 2); vP1 = get_val(f, ig, jg+1, kg, sym2, 2)
fyy(ii,jj,kk) = Sdydy * (vM1 - TWO*f(ig,jg,kg) + vP1)
endif
endif
! FZZ
if (kg+2 <= kmax .and. kg-2 >= kmin) then
if (kg > 2 .and. kg < ex(3)-1) then
fzz(ii,jj,kk) = Fdzdz * (-f(ig,jg,kg-2) + F16*f(ig,jg,kg-1) - F30*f(ig,jg,kg) + F16*f(ig,jg,kg+1) - f(ig,jg,kg+2))
else
vM2 = get_val(f, ig, jg, kg-2, sym3, 3); vM1 = get_val(f, ig, jg, kg-1, sym3, 3)
vP1 = get_val(f, ig, jg, kg+1, sym3, 3); vP2 = get_val(f, ig, jg, kg+2, sym3, 3)
fzz(ii,jj,kk) = Fdzdz * (-vM2 + F16*vM1 - F30*f(ig,jg,kg) + F16*vP1 - vP2)
endif
elseif (kg+1 <= kmax .and. kg-1 >= kmin) then
if (kg > 1 .and. kg < ex(3)) then
fzz(ii,jj,kk) = Sdzdz * (f(ig,jg,kg-1) - TWO*f(ig,jg,kg) + f(ig,jg,kg+1))
else
vM1 = get_val(f, ig, jg, kg-1, sym3, 3); vP1 = get_val(f, ig, jg, kg+1, sym3, 3)
fzz(ii,jj,kk) = Sdzdz * (vM1 - TWO*f(ig,jg,kg) + vP1)
endif
endif
! FXY
if (ig+2 <= imax .and. ig-2 >= imin .and. jg+2 <= jmax .and. jg-2 >= jmin) then
! 4th order mixed (simplification: skip detailed get_val for mixed 4th order bounds to save complexity,
! most points are inner. If boundary, fallback to 2nd order)
if (ig > 2 .and. ig < ex(1)-1 .and. jg > 2 .and. jg < ex(2)-1) then
fxy(ii,jj,kk) = Fdxdy * ( (f(ig-2,jg-2,kg)-EIGHT*f(ig-1,jg-2,kg)+EIGHT*f(ig+1,jg-2,kg)-f(ig+2,jg-2,kg)) &
-EIGHT*(f(ig-2,jg-1,kg)-EIGHT*f(ig-1,jg-1,kg)+EIGHT*f(ig+1,jg-1,kg)-f(ig+2,jg-1,kg)) &
+EIGHT*(f(ig-2,jg+1,kg)-EIGHT*f(ig-1,jg+1,kg)+EIGHT*f(ig+1,jg+1,kg)-f(ig+2,jg+1,kg)) &
- (f(ig-2,jg+2,kg)-EIGHT*f(ig-1,jg+2,kg)+EIGHT*f(ig+1,jg+2,kg)-f(ig+2,jg+2,kg)) )
else
! Fallback to 2nd order if near boundary (simpler than full 4th order get_val expansion)
! Check 2nd order bounds
if (ig+1 <= imax .and. ig-1 >= imin .and. jg+1 <= jmax .and. jg-1 >= jmin) then
! 2nd order mixed
vM1M1 = get_val(f, ig-1, jg-1, kg, sym1, 1) ! Sym? Mixed is tricky. Use simplified get_val or assume symmetry handled roughly
! Actually get_val logic needs dir. Mixed deriv uses 'corners'.
! Legacy uses fh.
! Let's implement safe 2nd order using get_val calls (corners might need double get_val?)
! get_val handles boundaries one by one.
! Corner f(i-1, j-1) -> if i<1, reflect i. Then if j<1, reflect j.
! get_val only handles ONE direction.
! However, corners are rare (only at edges/corners of grid).
! For safety, we can return 0 or do best effort.
! Let's do best effort: Assume inner or simple bounds.
! If strictly inner 2nd order:
if (ig > 1 .and. ig < ex(1) .and. jg > 1 .and. jg < ex(2)) then
fxy(ii,jj,kk) = Sdxdy * (f(ig-1,jg-1,kg) - f(ig+1,jg-1,kg) - f(ig-1,jg+1,kg) + f(ig+1,jg+1,kg))
else
fxy(ii,jj,kk) = ZEO ! Skip corners for safety/simplicity to avoid NaN
endif
endif
endif
elseif (ig+1 <= imax .and. ig-1 >= imin .and. jg+1 <= jmax .and. jg-1 >= jmin) then
if (ig > 1 .and. ig < ex(1) .and. jg > 1 .and. jg < ex(2)) then
fxy(ii,jj,kk) = Sdxdy * (f(ig-1,jg-1,kg) - f(ig+1,jg-1,kg) - f(ig-1,jg+1,kg) + f(ig+1,jg+1,kg))
else
fxy(ii,jj,kk) = ZEO
endif
endif
! FXZ
if (ig+2 <= imax .and. ig-2 >= imin .and. kg+2 <= kmax .and. kg-2 >= kmin) then
if (ig > 2 .and. ig < ex(1)-1 .and. kg > 2 .and. kg < ex(3)-1) then
fxz(ii,jj,kk) = Fdxdz * ( (f(ig-2,jg,kg-2)-EIGHT*f(ig-1,jg,kg-2)+EIGHT*f(ig+1,jg,kg-2)-f(ig+2,jg,kg-2)) &
-EIGHT*(f(ig-2,jg,kg-1)-EIGHT*f(ig-1,jg,kg-1)+EIGHT*f(ig+1,jg,kg-1)-f(ig+2,jg,kg-1)) &
+EIGHT*(f(ig-2,jg,kg+1)-EIGHT*f(ig-1,jg,kg+1)+EIGHT*f(ig+1,jg,kg+1)-f(ig+2,jg,kg+1)) &
- (f(ig-2,jg,kg+2)-EIGHT*f(ig-1,jg,kg+2)+EIGHT*f(ig+1,jg,kg+2)-f(ig+2,jg,kg+2)) )
else
if (ig+1 <= imax .and. ig-1 >= imin .and. kg+1 <= kmax .and. kg-1 >= kmin) then
if (ig > 1 .and. ig < ex(1) .and. kg > 1 .and. kg < ex(3)) then
fxz(ii,jj,kk) = Sdxdz * (f(ig-1,jg,kg-1) - f(ig+1,jg,kg-1) - f(ig-1,jg,kg+1) + f(ig+1,jg,kg+1))
else
fxz(ii,jj,kk) = ZEO
endif
endif
endif
elseif (ig+1 <= imax .and. ig-1 >= imin .and. kg+1 <= kmax .and. kg-1 >= kmin) then
if (ig > 1 .and. ig < ex(1) .and. kg > 1 .and. kg < ex(3)) then
fxz(ii,jj,kk) = Sdxdz * (f(ig-1,jg,kg-1) - f(ig+1,jg,kg-1) - f(ig-1,jg,kg+1) + f(ig+1,jg,kg+1))
else
fxz(ii,jj,kk) = ZEO
endif
endif
! FYZ
if (jg+2 <= jmax .and. jg-2 >= jmin .and. kg+2 <= kmax .and. kg-2 >= kmin) then
if (jg > 2 .and. jg < ex(2)-1 .and. kg > 2 .and. kg < ex(3)-1) then
fyz(ii,jj,kk) = Fdydz * ( (f(ig,jg-2,kg-2)-EIGHT*f(ig,jg-1,kg-2)+EIGHT*f(ig,jg+1,kg-2)-f(ig,jg+2,kg-2)) &
-EIGHT*(f(ig,jg-2,kg-1)-EIGHT*f(ig,jg-1,kg-1)+EIGHT*f(ig,jg+1,kg-1)-f(ig,jg+2,kg-1)) &
+EIGHT*(f(ig,jg-2,kg+1)-EIGHT*f(ig,jg-1,kg+1)+EIGHT*f(ig,jg+1,kg+1)-f(ig,jg+2,kg+1)) &
- (f(ig,jg-2,kg+2)-EIGHT*f(ig,jg-1,kg+2)+EIGHT*f(ig,jg+1,kg+2)-f(ig,jg+2,kg+2)) )
else
if (jg+1 <= jmax .and. jg-1 >= jmin .and. kg+1 <= kmax .and. kg-1 >= kmin) then
if (jg > 1 .and. jg < ex(2) .and. kg > 1 .and. kg < ex(3)) then
fyz(ii,jj,kk) = Sdydz * (f(ig,jg-1,kg-1) - f(ig,jg+1,kg-1) - f(ig,jg-1,kg+1) + f(ig,jg+1,kg+1))
else
fyz(ii,jj,kk) = ZEO
endif
endif
endif
elseif (jg+1 <= jmax .and. jg-1 >= jmin .and. kg+1 <= kmax .and. kg-1 >= kmin) then
if (jg > 1 .and. jg < ex(2) .and. kg > 1 .and. kg < ex(3)) then
fyz(ii,jj,kk) = Sdydz * (f(ig,jg-1,kg-1) - f(ig,jg+1,kg-1) - f(ig,jg-1,kg+1) + f(ig,jg+1,kg+1))
else
fyz(ii,jj,kk) = ZEO
endif
endif
end do; end do; end do
end subroutine calc_dderivs
!-------------------------------------------------------------------------
! Helper: Lopsided Derivative
!-------------------------------------------------------------------------
subroutine apply_lopsided(f, f_rhs, bx, by, bz, is, ie, js, je, ks, ke)
real*8, dimension(ex(1),ex(2),ex(3)), intent(in) :: f, bx, by, bz
real*8, dimension(ex(1),ex(2),ex(3)), intent(inout) :: f_rhs
integer, intent(in) :: is, ie, js, je, ks, ke
integer :: ii, jj, kk, ig, jg, kg
real*8 :: d12dx, d12dy, d12dz
d12dx = ONE / (12.d0 * (X(2)-X(1)))
d12dy = ONE / (12.d0 * (Y(2)-Y(1)))
d12dz = ONE / (12.d0 * (Z(2)-Z(1)))
!DIR$ SIMD
do kk=1,ke-ks+1; do jj=1,je-js+1; do ii=1,ie-is+1
ig = is + ii - 1; jg = js + jj - 1; kg = ks + kk - 1
! X-Direction
if (bx(ii,jj,kk) > ZEO) then
if (ig > 3 .and. ig < ex(1)) then
f_rhs(ig,jg,kg) = f_rhs(ig,jg,kg) + bx(ii,jj,kk)*d12dx*(-F3*f(ig-1,jg,kg)-F10*f(ig,jg,kg)+F18*f(ig+1,jg,kg)-F6*f(ig+2,jg,kg)+f(ig+3,jg,kg))
endif
elseif (bx(ii,jj,kk) < ZEO) then
if (ig > 0 .and. ig < ex(1)-3) then
f_rhs(ig,jg,kg) = f_rhs(ig,jg,kg) - bx(ii,jj,kk)*d12dx*(-F3*f(ig+1,jg,kg)-F10*f(ig,jg,kg)+F18*f(ig-1,jg,kg)-F6*f(ig-2,jg,kg)+f(ig-3,jg,kg))
endif
endif
! Y-Direction
if (by(ii,jj,kk) > ZEO) then
if (jg > 3 .and. jg < ex(2)) then
f_rhs(ig,jg,kg) = f_rhs(ig,jg,kg) + by(ii,jj,kk)*d12dy*(-F3*f(ig,jg-1,kg)-F10*f(ig,jg,kg)+F18*f(ig,jg+1,kg)-F6*f(ig,jg+2,kg)+f(ig,jg+3,kg))
endif
elseif (by(ii,jj,kk) < ZEO) then
if (jg > 0 .and. jg < ex(2)-3) then
f_rhs(ig,jg,kg) = f_rhs(ig,jg,kg) - by(ii,jj,kk)*d12dy*(-F3*f(ig,jg+1,kg)-F10*f(ig,jg,kg)+F18*f(ig,jg-1,kg)-F6*f(ig,jg-2,kg)+f(ig,jg-3,kg))
endif
endif
! Z-Direction
if (bz(ii,jj,kk) > ZEO) then
if (kg > 3 .and. kg < ex(3)) then
f_rhs(ig,jg,kg) = f_rhs(ig,jg,kg) + bz(ii,jj,kk)*d12dz*(-F3*f(ig,jg,kg-1)-F10*f(ig,jg,kg)+F18*f(ig,jg,kg+1)-F6*f(ig,jg,kg+2)+f(ig,jg,kg+3))
endif
elseif (bz(ii,jj,kk) < ZEO) then
if (kg > 0 .and. kg < ex(3)-3) then
f_rhs(ig,jg,kg) = f_rhs(ig,jg,kg) - bz(ii,jj,kk)*d12dz*(-F3*f(ig,jg,kg+1)-F10*f(ig,jg,kg)+F18*f(ig,jg,kg-1)-F6*f(ig,jg,kg-2)+f(ig,jg,kg-3))
endif
endif
end do; end do; end do
end subroutine apply_lopsided
!-------------------------------------------------------------------------
! Helper: Dissipation
!-------------------------------------------------------------------------
subroutine apply_kodiss(f, f_rhs, is, ie, js, je, ks, ke)
real*8, dimension(ex(1),ex(2),ex(3)), intent(in) :: f
real*8, dimension(ex(1),ex(2),ex(3)), intent(inout) :: f_rhs
integer, intent(in) :: is, ie, js, je, ks, ke
integer :: ii, jj, kk, ig, jg, kg
real*8 :: dX, dY, dZ
dX=X(2)-X(1); dY=Y(2)-Y(1); dZ=Z(2)-Z(1)
!DIR$ SIMD
do kk=1,ke-ks+1; do jj=1,je-js+1; do ii=1,ie-is+1
ig = is + ii - 1; jg = js + jj - 1; kg = ks + kk - 1
if (ig > 3 .and. ig < ex(1)-3 .and. jg > 3 .and. jg < ex(2)-3 .and. kg > 3 .and. kg < ex(3)-3) then
! X
f_rhs(ig,jg,kg) = f_rhs(ig,jg,kg) + eps/dX/F64 * ((f(ig-3,jg,kg)+f(ig+3,jg,kg)) - SIX*(f(ig-2,jg,kg)+f(ig+2,jg,kg)) + FIT*(f(ig-1,jg,kg)+f(ig+1,jg,kg)) - TWT*f(ig,jg,kg))
! Y
f_rhs(ig,jg,kg) = f_rhs(ig,jg,kg) + eps/dY/F64 * ((f(ig,jg-3,kg)+f(ig,jg+3,kg)) - SIX*(f(ig,jg-2,kg)+f(ig,jg+2,kg)) + FIT*(f(ig,jg-1,kg)+f(ig,jg+1,kg)) - TWT*f(ig,jg,kg))
! Z
f_rhs(ig,jg,kg) = f_rhs(ig,jg,kg) + eps/dZ/F64 * ((f(ig,jg,kg-3)+f(ig,jg,kg+3)) - SIX*(f(ig,jg,kg-2)+f(ig,jg,kg+2)) + FIT*(f(ig,jg,kg-1)+f(ig,jg,kg+1)) - TWT*f(ig,jg,kg))
endif
end do; end do; end do
end subroutine apply_kodiss
!-------------------------------------------------------------------------
! Helper: Boundary Value Getter
!-------------------------------------------------------------------------
function get_val(f, i, j, k, sym, dir) result(val)
real*8, dimension(ex(1),ex(2),ex(3)), intent(in) :: f
integer, intent(in) :: i, j, k, dir
real*8, intent(in) :: sym
real*8 :: val
integer :: ii, jj, kk
ii = i; jj = j; kk = k
val = 0.0d0
! Handle Lower Boundary
if (dir == 1 .and. i < 1) then
ii = 1 + (1 - i); val = sym
elseif (dir == 2 .and. j < 1) then
jj = 1 + (1 - j); val = sym
elseif (dir == 3 .and. k < 1) then
kk = 1 + (1 - k); val = sym
else
val = 1.0d0 ! No symmetry flip needed
endif
! Basic bounds check
if (ii > 0 .and. ii <= ex(1) .and. jj > 0 .and. jj <= ex(2) .and. kk > 0 .and. kk <= ex(3)) then
val = val * f(ii, jj, kk)
else
val = 0.0d0 ! Fallback
endif
end function get_val
end subroutine compute_rhs_bssn_opt
end module bssn_rhs_opt_mod
Reference in New Issue View Git Blame Copy Permalink