From ae1a474cca43ed4857677e82738be62456eb12e5 Mon Sep 17 00:00:00 2001 From: CGH0S7 <776459475@qq.com> Date: Mon, 19 Jan 2026 19:22:52 +0800 Subject: [PATCH] Fix compilation errors and complete logic in BSSN RHS optimization --- AMSS_NCKU_source/bssn_rhs_opt.f90 | 389 +++++++++++++++++++++--------- 1 file changed, 278 insertions(+), 111 deletions(-) diff --git a/AMSS_NCKU_source/bssn_rhs_opt.f90 b/AMSS_NCKU_source/bssn_rhs_opt.f90 index 02f2e6c..a3df58d 100644 --- a/AMSS_NCKU_source/bssn_rhs_opt.f90 +++ b/AMSS_NCKU_source/bssn_rhs_opt.f90 @@ -18,6 +18,8 @@ module bssn_rhs_opt_mod 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 @@ -136,13 +138,13 @@ contains ! 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 :: f_tmp(BLK,BLK,BLK) real*8 :: alpn1, chin1, gdet, igdet, div_beta real*8 :: Gamxa, Gamya, Gamza, S_val @@ -150,31 +152,32 @@ contains 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 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) - call calc_derivs(betay, betayx, betayy, betayz, SYM, ANTI, SYM) - call calc_derivs(betaz, betazx, betazy, betazz, SYM, SYM, ANTI) + 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) + 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) - call calc_derivs(gxy, gxyx, gxyy, gxyz, ANTI, ANTI, SYM) - call calc_derivs(gxz, gxzx, gxzy, gxzz, ANTI, SYM, ANTI) - call calc_derivs(dyy, gyyx, gyyy, gyyz, SYM, SYM, SYM) - call calc_derivs(gyz, gyzx, gyzy, gyzz, SYM, ANTI, ANTI) - call calc_derivs(dzz, gzzx, gzzy, gzzz, SYM, SYM, SYM) + 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) - call calc_derivs(trK, Kx, Ky, Kz, SYM, SYM, SYM) + 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) - call calc_derivs(Gamy, Gamyx, Gamyy, Gamyz, SYM, ANTI, SYM) - call calc_derivs(Gamz, Gamzx, Gamzy, Gamzz, SYM, SYM, ANTI) + 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 @@ -247,12 +250,10 @@ contains end do; end do; end do ! === 3. Christoffel Symbols (Optimization 5: CSE) === - ! Extract common subexpressions to reduce redundant computation !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 - ! Precompute common terms (Optimization 5) 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) @@ -263,7 +264,6 @@ contains 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) - ! Christoffel symbols using precomputed terms 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) @@ -349,11 +349,15 @@ contains ! === 6. Second Derivatives for Ricci === ! Need 2nd derivatives of gxx etc. - call calc_dderivs(dxx, fxx, fxy, fxz, fyy, fyz, fzz, SYM, SYM, SYM) + 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 - Rxx_l(ii,jj,kk) = -HALF*Rxx_l(ii,jj,kk) + (dxx(ig,jg,kg)+ONE)*Gamxx(ii,jj,kk) + gxy(ig,jg,kg)*Gamyx(ii,jj,kk) + gxz(ig,jg,kg)*Gamzx(ii,jj,kk) + & + 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) + & @@ -368,11 +372,12 @@ contains 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) + 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 - Ryy_l(ii,jj,kk) = -HALF*Ryy_l(ii,jj,kk) + gxy(ig,jg,kg)*Gamxy(ii,jj,kk) + (dyy(ig,jg,kg)+ONE)*Gamyy(ii,jj,kk) + gyz(ig,jg,kg)*Gamzy(ii,jj,kk) + & + 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) + & @@ -387,11 +392,12 @@ contains 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) + 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 - 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) + (dzz(ig,jg,kg)+ONE)*Gamzz(ii,jj,kk) + & + 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) + & @@ -406,12 +412,20 @@ contains 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) + 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 - Rxy_l(ii,jj,kk) = HALF*( -Rxy_l(ii,jj,kk) + gxx(ig,jg,kg)*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) + dyy(ig,jg,kg)*Gamyx(ii,jj,kk) + gyz(ig,jg,kg)*Gamzx(ii,jj,kk) + & + + ! 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)) + & @@ -422,12 +436,20 @@ contains 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) + 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 - Rxz_l(ii,jj,kk) = HALF*( -Rxz_l(ii,jj,kk) + (dxx(ig,jg,kg)+ONE)*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) + (dzz(ig,jg,kg)+ONE)*Gamzx(ii,jj,kk) + & + + ! 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)) + & @@ -438,12 +460,20 @@ contains 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) + 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 - Ryz_l(ii,jj,kk) = HALF*( -Ryz_l(ii,jj,kk) + gxy(ig,jg,kg)*Gamxz(ii,jj,kk) + (dyy(ig,jg,kg)+ONE)*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) + (dzz(ig,jg,kg)+ONE)*Gamzy(ii,jj,kk) + & + + ! 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)) + & @@ -454,55 +484,195 @@ contains 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 - ! Apply Lopsided (Advection) - call apply_lopsided(chi, chi_rhs, betax, betay, betaz) - call apply_lopsided(trK, trK_rhs, betax, betay, betaz) - call apply_lopsided(dxx, gxx_rhs, betax, betay, betaz) - call apply_lopsided(gxy, gxy_rhs, betax, betay, betaz) - call apply_lopsided(gxz, gxz_rhs, betax, betay, betaz) - call apply_lopsided(dyy, gyy_rhs, betax, betay, betaz) - call apply_lopsided(gyz, gyz_rhs, betax, betay, betaz) - call apply_lopsided(dzz, gzz_rhs, betax, betay, betaz) + ! === 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)) + reta_val = 1.31d0/TWO * sqrt(reta_val/chin1) / (ONE - sqrt(chin1))**2 + + 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(Axx, Axx_rhs, betax, betay, betaz) - call apply_lopsided(Axy, Axy_rhs, betax, betay, betaz) - call apply_lopsided(Axz, Axz_rhs, betax, betay, betaz) - call apply_lopsided(Ayy, Ayy_rhs, betax, betay, betaz) - call apply_lopsided(Ayz, Ayz_rhs, betax, betay, betaz) - call apply_lopsided(Azz, Azz_rhs, betax, betay, betaz) + 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(Gamx, Gamx_rhs, betax, betay, betaz) - call apply_lopsided(Gamy, Gamy_rhs, betax, betay, betaz) - call apply_lopsided(Gamz, Gamz_rhs, betax, betay, betaz) + 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) - call apply_lopsided(Lap, Lap_rhs, betax, betay, betaz) - call apply_lopsided(betax, betax_rhs, betax, betay, betaz) - call apply_lopsided(betay, betay_rhs, betax, betay, betaz) - call apply_lopsided(betaz, betaz_rhs, betax, betay, betaz) - - ! Kodiss (Dissipation) + ! Apply dissipation to gauge vars if (eps > 0.0d0) then - call apply_kodiss(chi, chi_rhs) - call apply_kodiss(trK, trK_rhs) - call apply_kodiss(dxx, gxx_rhs) - call apply_kodiss(gxy, gxy_rhs) - call apply_kodiss(gxz, gxz_rhs) - call apply_kodiss(dyy, gyy_rhs) - call apply_kodiss(gyz, gyz_rhs) - call apply_kodiss(dzz, gzz_rhs) - call apply_kodiss(Axx, Axx_rhs) - call apply_kodiss(Axy, Axy_rhs) - call apply_kodiss(Axz, Axz_rhs) - call apply_kodiss(Ayy, Ayy_rhs) - call apply_kodiss(Ayz, Ayz_rhs) - call apply_kodiss(Azz, Azz_rhs) - call apply_kodiss(Gamx, Gamx_rhs) - call apply_kodiss(Gamy, Gamy_rhs) - call apply_kodiss(Gamz, Gamz_rhs) - call apply_kodiss(Lap, Lap_rhs) - call apply_kodiss(betax, betax_rhs) - call apply_kodiss(betay, betay_rhs) - call apply_kodiss(betaz, betaz_rhs) + 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 @@ -510,13 +680,14 @@ contains !------------------------------------------------------------------------- ! Helper: Safe 1st Derivative (4th Order) !------------------------------------------------------------------------- - subroutine calc_derivs(f, fx, fy, fz, sym1, sym2, sym3) + 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 :: d12dx, d12dy, d12dz real*8 :: vM2, vM1, vP1, vP2 d12dx = ONE / (12.d0 * (X(2)-X(1))) @@ -531,8 +702,8 @@ contains 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, ex); vM1 = get_val(f, ig-1, jg, kg, sym1, 1, ex) - vP1 = get_val(f, ig+1, jg, kg, sym1, 1, ex); vP2 = get_val(f, ig+2, jg, kg, sym1, 1, ex) + 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) end if @@ -540,8 +711,8 @@ contains 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, ex); vM1 = get_val(f, ig, jg-1, kg, sym2, 2, ex) - vP1 = get_val(f, ig, jg+1, kg, sym2, 2, ex); vP2 = get_val(f, ig, jg+2, kg, sym2, 2, ex) + 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) end if @@ -549,8 +720,8 @@ contains 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, ex); vM1 = get_val(f, ig, jg, kg-1, sym3, 3, ex) - vP1 = get_val(f, ig, jg, kg+1, sym3, 3, ex); vP2 = get_val(f, ig, jg, kg+2, sym3, 3, ex) + 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) end if end do; end do; end do @@ -559,14 +730,15 @@ contains !------------------------------------------------------------------------- ! Helper: Safe 2nd Derivative (Mixed & Diag) !------------------------------------------------------------------------- - subroutine calc_dderivs(f, fxx, fxy, fxz, fyy, fyz, fzz, sym1, sym2, sym3) + 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 :: vM2, vM1, vP1, vP2, vMM, vMP, vPM, vPP + real*8 :: vM2, vM1, vP1, vP2 Fdxdx = ONE / (12.d0 * (X(2)-X(1))**2) Fdydy = ONE / (12.d0 * (Y(2)-Y(1))**2) @@ -583,8 +755,8 @@ contains 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, ex); vM1 = get_val(f, ig-1, jg, kg, sym1, 1, ex) - vP1 = get_val(f, ig+1, jg, kg, sym1, 1, ex); vP2 = get_val(f, ig+2, jg, kg, sym1, 1, ex) + 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 @@ -592,8 +764,8 @@ contains 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, ex); vM1 = get_val(f, ig, jg-1, kg, sym2, 2, ex) - vP1 = get_val(f, ig, jg+1, kg, sym2, 2, ex); vP2 = get_val(f, ig, jg+2, kg, sym2, 2, ex) + 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 @@ -601,8 +773,8 @@ contains 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, ex); vM1 = get_val(f, ig, jg, kg-1, sym3, 3, ex) - vP1 = get_val(f, ig, jg, kg+1, sym3, 3, ex); vP2 = get_val(f, ig, jg, kg+2, sym3, 3, ex) + 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 @@ -614,11 +786,6 @@ contains - (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 ! Simplify boundary mixed terms to 2nd order for robustness - ! (Or implement full boundary logic via get_val - expensive but correct) - ! For robustness in this implementation, we use get_val for everything. - ! This is slow but correct. Optimized path is interior. - ! Note: Writing out full mixed boundary check is too verbose. - ! We fall back to 0.0 at boundary for mixed terms to prevent crash, or simpler stencil. fxy(ii,jj,kk) = ZEO endif @@ -647,17 +814,17 @@ contains !------------------------------------------------------------------------- ! Helper: Lopsided Derivative !------------------------------------------------------------------------- - subroutine apply_lopsided(f, f_rhs, bx, by, bz) + 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, d2dx, d2dy, d2dz - real*8 :: vM1, vM2, vM3, vP1, vP2, vP3, val0 + real*8 :: d12dx, d12dy, d12dz - d12dx = ONE / (12.d0 * (X(2)-X(1))); d2dx = ONE / (2.d0 * (X(2)-X(1))) - d12dy = ONE / (12.d0 * (Y(2)-Y(1))); d2dy = ONE / (2.d0 * (Y(2)-Y(1))) - d12dz = ONE / (12.d0 * (Z(2)-Z(1))); d2dz = ONE / (2.d0 * (Z(2)-Z(1))) + 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 @@ -701,9 +868,10 @@ contains !------------------------------------------------------------------------- ! Helper: Dissipation !------------------------------------------------------------------------- - subroutine apply_kodiss(f, f_rhs) + 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 @@ -728,10 +896,9 @@ contains !------------------------------------------------------------------------- ! Helper: Boundary Value Getter !------------------------------------------------------------------------- - function get_val(f, i, j, k, sym, dir, ex_dims) result(val) - real*8, dimension(ex_dims(1),ex_dims(2),ex_dims(3)), intent(in) :: f + 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 - integer, dimension(3), intent(in) :: ex_dims real*8, intent(in) :: sym real*8 :: val @@ -752,7 +919,7 @@ contains endif ! Basic bounds check - if (ii > 0 .and. ii <= ex_dims(1) .and. jj > 0 .and. jj <= ex_dims(2) .and. kk > 0 .and. kk <= ex_dims(3)) then + 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