From f1fe9fd44365f1c29315793ac917acbd5119765e Mon Sep 17 00:00:00 2001 From: CGH0S7 <776459475@qq.com> Date: Tue, 3 Mar 2026 15:57:10 +0800 Subject: [PATCH] =?UTF-8?q?=E8=BF=81=E7=A7=BBC=E7=AE=97=E5=AD=90=E7=9A=84?= =?UTF-8?q?=E5=BE=AA=E7=8E=AF=E8=9E=8D=E5=90=88=E5=92=8C=E4=B8=B4=E6=97=B6?= =?UTF-8?q?=E9=87=8F=E6=B6=88=E9=99=A4?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- AMSS_NCKU_source/bssn_rhs.f90 | 372 +++++++++++++++++----------------- 1 file changed, 186 insertions(+), 186 deletions(-) diff --git a/AMSS_NCKU_source/bssn_rhs.f90 b/AMSS_NCKU_source/bssn_rhs.f90 index a18be11..2f386f6 100644 --- a/AMSS_NCKU_source/bssn_rhs.f90 +++ b/AMSS_NCKU_source/bssn_rhs.f90 @@ -91,6 +91,8 @@ real*8 :: Rxx_loc,Rxy_loc,Rxz_loc,Ryy_loc,Ryz_loc,Rzz_loc real*8 :: fxx_loc,fxy_loc,fxz_loc real*8 :: Gamxa_loc,Gamya_loc,Gamza_loc + real*8 :: f_loc,chin_loc + real*8 :: l_fxx,l_fxy,l_fxz,l_fyy,l_fyz,l_fzz,S_loc real*8, parameter :: ZEO = 0.d0,ONE = 1.D0, TWO = 2.D0, FOUR = 4.D0 real*8, parameter :: EIGHT = 8.D0, HALF = 0.5D0, THR = 3.d0 real*8, parameter :: SYM = 1.D0, ANTI= - 1.D0 @@ -653,192 +655,190 @@ Gamxyz * gxzz + Gamyyz * gyzz + Gamzyz * gzzz + & Gamxzz * gxzy + Gamyzz * gyzy + Gamzzz * gzzy + & Gamxyz * gzzx + Gamyyz * gzzy + Gamzyz * gzzz ) -!covariant second derivative of chi respect to tilted metric - call fdderivs(ex,chi,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) - - fxx = fxx - Gamxxx * chix - Gamyxx * chiy - Gamzxx * chiz - fxy = fxy - Gamxxy * chix - Gamyxy * chiy - Gamzxy * chiz - fxz = fxz - Gamxxz * chix - Gamyxz * chiy - Gamzxz * chiz - fyy = fyy - Gamxyy * chix - Gamyyy * chiy - Gamzyy * chiz - fyz = fyz - Gamxyz * chix - Gamyyz * chiy - Gamzyz * chiz - fzz = fzz - Gamxzz * chix - Gamyzz * chiy - Gamzzz * chiz -! Store D^l D_l chi - 3/(2*chi) D^l chi D_l chi in f - - f = gupxx * ( fxx - F3o2/chin1 * chix * chix ) + & - gupyy * ( fyy - F3o2/chin1 * chiy * chiy ) + & - gupzz * ( fzz - F3o2/chin1 * chiz * chiz ) + & - TWO * gupxy * ( fxy - F3o2/chin1 * chix * chiy ) + & - TWO * gupxz * ( fxz - F3o2/chin1 * chix * chiz ) + & - TWO * gupyz * ( fyz - F3o2/chin1 * chiy * chiz ) -! Add chi part to Ricci tensor: - - Rxx = Rxx + (fxx - chix*chix/chin1/TWO + gxx * f)/chin1/TWO - Ryy = Ryy + (fyy - chiy*chiy/chin1/TWO + gyy * f)/chin1/TWO - Rzz = Rzz + (fzz - chiz*chiz/chin1/TWO + gzz * f)/chin1/TWO - Rxy = Rxy + (fxy - chix*chiy/chin1/TWO + gxy * f)/chin1/TWO - Rxz = Rxz + (fxz - chix*chiz/chin1/TWO + gxz * f)/chin1/TWO - Ryz = Ryz + (fyz - chiy*chiz/chin1/TWO + gyz * f)/chin1/TWO - -! covariant second derivatives of the lapse respect to physical metric - call fdderivs(ex,Lap,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z, & - SYM,SYM,SYM,symmetry,Lev) - - gxxx = (gupxx * chix + gupxy * chiy + gupxz * chiz)/chin1 - gxxy = (gupxy * chix + gupyy * chiy + gupyz * chiz)/chin1 - gxxz = (gupxz * chix + gupyz * chiy + gupzz * chiz)/chin1 -! now get physical second kind of connection - Gamxxx = Gamxxx - ( (chix + chix)/chin1 - gxx * gxxx )*HALF - Gamyxx = Gamyxx - ( - gxx * gxxy )*HALF - Gamzxx = Gamzxx - ( - gxx * gxxz )*HALF - Gamxyy = Gamxyy - ( - gyy * gxxx )*HALF - Gamyyy = Gamyyy - ( (chiy + chiy)/chin1 - gyy * gxxy )*HALF - Gamzyy = Gamzyy - ( - gyy * gxxz )*HALF - Gamxzz = Gamxzz - ( - gzz * gxxx )*HALF - Gamyzz = Gamyzz - ( - gzz * gxxy )*HALF - Gamzzz = Gamzzz - ( (chiz + chiz)/chin1 - gzz * gxxz )*HALF - Gamxxy = Gamxxy - ( chiy /chin1 - gxy * gxxx )*HALF - Gamyxy = Gamyxy - ( chix /chin1 - gxy * gxxy )*HALF - Gamzxy = Gamzxy - ( - gxy * gxxz )*HALF - Gamxxz = Gamxxz - ( chiz /chin1 - gxz * gxxx )*HALF - Gamyxz = Gamyxz - ( - gxz * gxxy )*HALF - Gamzxz = Gamzxz - ( chix /chin1 - gxz * gxxz )*HALF - Gamxyz = Gamxyz - ( - gyz * gxxx )*HALF - Gamyyz = Gamyyz - ( chiz /chin1 - gyz * gxxy )*HALF - Gamzyz = Gamzyz - ( chiy /chin1 - gyz * gxxz )*HALF - - fxx = fxx - Gamxxx*Lapx - Gamyxx*Lapy - Gamzxx*Lapz - fyy = fyy - Gamxyy*Lapx - Gamyyy*Lapy - Gamzyy*Lapz - fzz = fzz - Gamxzz*Lapx - Gamyzz*Lapy - Gamzzz*Lapz - fxy = fxy - Gamxxy*Lapx - Gamyxy*Lapy - Gamzxy*Lapz - fxz = fxz - Gamxxz*Lapx - Gamyxz*Lapy - Gamzxz*Lapz - fyz = fyz - Gamxyz*Lapx - Gamyyz*Lapy - Gamzyz*Lapz - -! store D^i D_i Lap in trK_rhs upto chi - trK_rhs = gupxx * fxx + gupyy * fyy + gupzz * fzz + & - TWO* ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) -#if 1 -!! follow bam code - S = chin1 * ( gupxx * Sxx + gupyy * Syy + gupzz * Szz + & - TWO * ( gupxy * Sxy + gupxz * Sxz + gupyz * Syz ) ) - f = F2o3 * trK * trK -(& - gupxx * ( & - gupxx * Axx * Axx + gupyy * Axy * Axy + gupzz * Axz * Axz + & - TWO * (gupxy * Axx * Axy + gupxz * Axx * Axz + gupyz * Axy * Axz) ) + & - gupyy * ( & - gupxx * Axy * Axy + gupyy * Ayy * Ayy + gupzz * Ayz * Ayz + & - TWO * (gupxy * Axy * Ayy + gupxz * Axy * Ayz + gupyz * Ayy * Ayz) ) + & - gupzz * ( & - gupxx * Axz * Axz + gupyy * Ayz * Ayz + gupzz * Azz * Azz + & - TWO * (gupxy * Axz * Ayz + gupxz * Axz * Azz + gupyz * Ayz * Azz) ) + & - TWO * ( & - gupxy * ( & - gupxx * Axx * Axy + gupyy * Axy * Ayy + gupzz * Axz * Ayz + & - gupxy * (Axx * Ayy + Axy * Axy) + & - gupxz * (Axx * Ayz + Axz * Axy) + & - gupyz * (Axy * Ayz + Axz * Ayy) ) + & - gupxz * ( & - gupxx * Axx * Axz + gupyy * Axy * Ayz + gupzz * Axz * Azz + & - gupxy * (Axx * Ayz + Axy * Axz) + & - gupxz * (Axx * Azz + Axz * Axz) + & - gupyz * (Axy * Azz + Axz * Ayz) ) + & - gupyz * ( & - gupxx * Axy * Axz + gupyy * Ayy * Ayz + gupzz * Ayz * Azz + & - gupxy * (Axy * Ayz + Ayy * Axz) + & - gupxz * (Axy * Azz + Ayz * Axz) + & - gupyz * (Ayy * Azz + Ayz * Ayz) ) )) -1.6d1*PI*rho + EIGHT * PI * S - f = - F1o3 *( gupxx * fxx + gupyy * fyy + gupzz * fzz + & - TWO* ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) + alpn1/chin1*f) - - fxx = alpn1 * (Rxx - EIGHT * PI * Sxx) - fxx - fxy = alpn1 * (Rxy - EIGHT * PI * Sxy) - fxy - fxz = alpn1 * (Rxz - EIGHT * PI * Sxz) - fxz - fyy = alpn1 * (Ryy - EIGHT * PI * Syy) - fyy - fyz = alpn1 * (Ryz - EIGHT * PI * Syz) - fyz - fzz = alpn1 * (Rzz - EIGHT * PI * Szz) - fzz -#else -! Add lapse and S_ij parts to Ricci tensor: - - fxx = alpn1 * (Rxx - EIGHT * PI * Sxx) - fxx - fxy = alpn1 * (Rxy - EIGHT * PI * Sxy) - fxy - fxz = alpn1 * (Rxz - EIGHT * PI * Sxz) - fxz - fyy = alpn1 * (Ryy - EIGHT * PI * Syy) - fyy - fyz = alpn1 * (Ryz - EIGHT * PI * Syz) - fyz - fzz = alpn1 * (Rzz - EIGHT * PI * Szz) - fzz - -! Compute trace-free part (note: chi^-1 and chi cancel!): - - f = F1o3 *( gupxx * fxx + gupyy * fyy + gupzz * fzz + & - TWO* ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) ) -#endif - - Axx_rhs = fxx - gxx * f - Ayy_rhs = fyy - gyy * f - Azz_rhs = fzz - gzz * f - Axy_rhs = fxy - gxy * f - Axz_rhs = fxz - gxz * f - Ayz_rhs = fyz - gyz * f - -! Now: store A_il A^l_j into fij: - - fxx = gupxx * Axx * Axx + gupyy * Axy * Axy + gupzz * Axz * Axz + & - TWO * (gupxy * Axx * Axy + gupxz * Axx * Axz + gupyz * Axy * Axz) - fyy = gupxx * Axy * Axy + gupyy * Ayy * Ayy + gupzz * Ayz * Ayz + & - TWO * (gupxy * Axy * Ayy + gupxz * Axy * Ayz + gupyz * Ayy * Ayz) - fzz = gupxx * Axz * Axz + gupyy * Ayz * Ayz + gupzz * Azz * Azz + & - TWO * (gupxy * Axz * Ayz + gupxz * Axz * Azz + gupyz * Ayz * Azz) - fxy = gupxx * Axx * Axy + gupyy * Axy * Ayy + gupzz * Axz * Ayz + & - gupxy *(Axx * Ayy + Axy * Axy) + & - gupxz *(Axx * Ayz + Axz * Axy) + & - gupyz *(Axy * Ayz + Axz * Ayy) - fxz = gupxx * Axx * Axz + gupyy * Axy * Ayz + gupzz * Axz * Azz + & - gupxy *(Axx * Ayz + Axy * Axz) + & - gupxz *(Axx * Azz + Axz * Axz) + & - gupyz *(Axy * Azz + Axz * Ayz) - fyz = gupxx * Axy * Axz + gupyy * Ayy * Ayz + gupzz * Ayz * Azz + & - gupxy *(Axy * Ayz + Ayy * Axz) + & - gupxz *(Axy * Azz + Ayz * Axz) + & - gupyz *(Ayy * Azz + Ayz * Ayz) - - f = chin1 -! store D^i D_i Lap in trK_rhs - trK_rhs = f*trK_rhs - - Axx_rhs = f * Axx_rhs+ alpn1 * (trK * Axx - TWO * fxx) + & - TWO * ( Axx * betaxx + Axy * betayx + Axz * betazx )- & - F2o3 * Axx * div_beta - - Ayy_rhs = f * Ayy_rhs+ alpn1 * (trK * Ayy - TWO * fyy) + & - TWO * ( Axy * betaxy + Ayy * betayy + Ayz * betazy )- & - F2o3 * Ayy * div_beta - - Azz_rhs = f * Azz_rhs+ alpn1 * (trK * Azz - TWO * fzz) + & - TWO * ( Axz * betaxz + Ayz * betayz + Azz * betazz )- & - F2o3 * Azz * div_beta - - Axy_rhs = f * Axy_rhs+ alpn1 *( trK * Axy - TWO * fxy )+ & - Axx * betaxy + Axz * betazy + & - Ayy * betayx + Ayz * betazx + & - F1o3 * Axy * div_beta - Axy * betazz - - Ayz_rhs = f * Ayz_rhs+ alpn1 *( trK * Ayz - TWO * fyz )+ & - Axy * betaxz + Ayy * betayz + & - Axz * betaxy + Azz * betazy + & - F1o3 * Ayz * div_beta - Ayz * betaxx - - Axz_rhs = f * Axz_rhs+ alpn1 *( trK * Axz - TWO * fxz )+ & - Axx * betaxz + Axy * betayz + & - Ayz * betayx + Azz * betazx + & - F1o3 * Axz * div_beta - Axz * betayy !rhs for Aij - -! Compute trace of S_ij - - S = f * ( gupxx * Sxx + gupyy * Syy + gupzz * Szz + & - TWO * ( gupxy * Sxy + gupxz * Sxz + gupyz * Syz ) ) - - trK_rhs = - trK_rhs + alpn1 *( F1o3 * trK * trK + & - gupxx * fxx + gupyy * fyy + gupzz * fzz + & - TWO * ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) + & - FOUR * PI * ( rho + S )) !rhs for trK +!covariant second derivative of chi respect to tilted metric + call fdderivs(ex,chi,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev) + + do k=1,ex(3) + do j=1,ex(2) + do i=1,ex(1) + fxx(i,j,k) = fxx(i,j,k) - Gamxxx(i,j,k) * chix(i,j,k) - Gamyxx(i,j,k) * chiy(i,j,k) - Gamzxx(i,j,k) * chiz(i,j,k) + fxy(i,j,k) = fxy(i,j,k) - Gamxxy(i,j,k) * chix(i,j,k) - Gamyxy(i,j,k) * chiy(i,j,k) - Gamzxy(i,j,k) * chiz(i,j,k) + fxz(i,j,k) = fxz(i,j,k) - Gamxxz(i,j,k) * chix(i,j,k) - Gamyxz(i,j,k) * chiy(i,j,k) - Gamzxz(i,j,k) * chiz(i,j,k) + fyy(i,j,k) = fyy(i,j,k) - Gamxyy(i,j,k) * chix(i,j,k) - Gamyyy(i,j,k) * chiy(i,j,k) - Gamzyy(i,j,k) * chiz(i,j,k) + fyz(i,j,k) = fyz(i,j,k) - Gamxyz(i,j,k) * chix(i,j,k) - Gamyyz(i,j,k) * chiy(i,j,k) - Gamzyz(i,j,k) * chiz(i,j,k) + fzz(i,j,k) = fzz(i,j,k) - Gamxzz(i,j,k) * chix(i,j,k) - Gamyzz(i,j,k) * chiy(i,j,k) - Gamzzz(i,j,k) * chiz(i,j,k) + + chin_loc = chin1(i,j,k) + f_loc = gupxx(i,j,k) * (fxx(i,j,k) - F3o2/chin_loc * chix(i,j,k) * chix(i,j,k)) + & + gupyy(i,j,k) * (fyy(i,j,k) - F3o2/chin_loc * chiy(i,j,k) * chiy(i,j,k)) + & + gupzz(i,j,k) * (fzz(i,j,k) - F3o2/chin_loc * chiz(i,j,k) * chiz(i,j,k)) + & + TWO * gupxy(i,j,k) * (fxy(i,j,k) - F3o2/chin_loc * chix(i,j,k) * chiy(i,j,k)) + & + TWO * gupxz(i,j,k) * (fxz(i,j,k) - F3o2/chin_loc * chix(i,j,k) * chiz(i,j,k)) + & + TWO * gupyz(i,j,k) * (fyz(i,j,k) - F3o2/chin_loc * chiy(i,j,k) * chiz(i,j,k)) + f(i,j,k) = f_loc + + Rxx(i,j,k) = Rxx(i,j,k) + (fxx(i,j,k) - chix(i,j,k)*chix(i,j,k)/chin_loc/TWO + gxx(i,j,k) * f_loc)/chin_loc/TWO + Ryy(i,j,k) = Ryy(i,j,k) + (fyy(i,j,k) - chiy(i,j,k)*chiy(i,j,k)/chin_loc/TWO + gyy(i,j,k) * f_loc)/chin_loc/TWO + Rzz(i,j,k) = Rzz(i,j,k) + (fzz(i,j,k) - chiz(i,j,k)*chiz(i,j,k)/chin_loc/TWO + gzz(i,j,k) * f_loc)/chin_loc/TWO + Rxy(i,j,k) = Rxy(i,j,k) + (fxy(i,j,k) - chix(i,j,k)*chiy(i,j,k)/chin_loc/TWO + gxy(i,j,k) * f_loc)/chin_loc/TWO + Rxz(i,j,k) = Rxz(i,j,k) + (fxz(i,j,k) - chix(i,j,k)*chiz(i,j,k)/chin_loc/TWO + gxz(i,j,k) * f_loc)/chin_loc/TWO + Ryz(i,j,k) = Ryz(i,j,k) + (fyz(i,j,k) - chiy(i,j,k)*chiz(i,j,k)/chin_loc/TWO + gyz(i,j,k) * f_loc)/chin_loc/TWO + enddo + enddo + enddo + +! covariant second derivatives of the lapse respect to physical metric + call fdderivs(ex,Lap,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z, & + SYM,SYM,SYM,symmetry,Lev) + + do k=1,ex(3) + do j=1,ex(2) + do i=1,ex(1) + chin_loc = chin1(i,j,k) + gxxx(i,j,k) = (gupxx(i,j,k) * chix(i,j,k) + gupxy(i,j,k) * chiy(i,j,k) + gupxz(i,j,k) * chiz(i,j,k)) / chin_loc + gxxy(i,j,k) = (gupxy(i,j,k) * chix(i,j,k) + gupyy(i,j,k) * chiy(i,j,k) + gupyz(i,j,k) * chiz(i,j,k)) / chin_loc + gxxz(i,j,k) = (gupxz(i,j,k) * chix(i,j,k) + gupyz(i,j,k) * chiy(i,j,k) + gupzz(i,j,k) * chiz(i,j,k)) / chin_loc + + Gamxxx(i,j,k) = Gamxxx(i,j,k) - ( (chix(i,j,k) + chix(i,j,k))/chin_loc - gxx(i,j,k) * gxxx(i,j,k) )*HALF + Gamyxx(i,j,k) = Gamyxx(i,j,k) - ( - gxx(i,j,k) * gxxy(i,j,k) )*HALF + Gamzxx(i,j,k) = Gamzxx(i,j,k) - ( - gxx(i,j,k) * gxxz(i,j,k) )*HALF + Gamxyy(i,j,k) = Gamxyy(i,j,k) - ( - gyy(i,j,k) * gxxx(i,j,k) )*HALF + Gamyyy(i,j,k) = Gamyyy(i,j,k) - ( (chiy(i,j,k) + chiy(i,j,k))/chin_loc - gyy(i,j,k) * gxxy(i,j,k) )*HALF + Gamzyy(i,j,k) = Gamzyy(i,j,k) - ( - gyy(i,j,k) * gxxz(i,j,k) )*HALF + Gamxzz(i,j,k) = Gamxzz(i,j,k) - ( - gzz(i,j,k) * gxxx(i,j,k) )*HALF + Gamyzz(i,j,k) = Gamyzz(i,j,k) - ( - gzz(i,j,k) * gxxy(i,j,k) )*HALF + Gamzzz(i,j,k) = Gamzzz(i,j,k) - ( (chiz(i,j,k) + chiz(i,j,k))/chin_loc - gzz(i,j,k) * gxxz(i,j,k) )*HALF + Gamxxy(i,j,k) = Gamxxy(i,j,k) - ( chiy(i,j,k) /chin_loc - gxy(i,j,k) * gxxx(i,j,k) )*HALF + Gamyxy(i,j,k) = Gamyxy(i,j,k) - ( chix(i,j,k) /chin_loc - gxy(i,j,k) * gxxy(i,j,k) )*HALF + Gamzxy(i,j,k) = Gamzxy(i,j,k) - ( - gxy(i,j,k) * gxxz(i,j,k) )*HALF + Gamxxz(i,j,k) = Gamxxz(i,j,k) - ( chiz(i,j,k) /chin_loc - gxz(i,j,k) * gxxx(i,j,k) )*HALF + Gamyxz(i,j,k) = Gamyxz(i,j,k) - ( - gxz(i,j,k) * gxxy(i,j,k) )*HALF + Gamzxz(i,j,k) = Gamzxz(i,j,k) - ( chix(i,j,k) /chin_loc - gxz(i,j,k) * gxxz(i,j,k) )*HALF + Gamxyz(i,j,k) = Gamxyz(i,j,k) - ( - gyz(i,j,k) * gxxx(i,j,k) )*HALF + Gamyyz(i,j,k) = Gamyyz(i,j,k) - ( chiz(i,j,k) /chin_loc - gyz(i,j,k) * gxxy(i,j,k) )*HALF + Gamzyz(i,j,k) = Gamzyz(i,j,k) - ( chiy(i,j,k) /chin_loc - gyz(i,j,k) * gxxz(i,j,k) )*HALF + + fxx(i,j,k) = fxx(i,j,k) - Gamxxx(i,j,k)*Lapx(i,j,k) - Gamyxx(i,j,k)*Lapy(i,j,k) - Gamzxx(i,j,k)*Lapz(i,j,k) + fyy(i,j,k) = fyy(i,j,k) - Gamxyy(i,j,k)*Lapx(i,j,k) - Gamyyy(i,j,k)*Lapy(i,j,k) - Gamzyy(i,j,k)*Lapz(i,j,k) + fzz(i,j,k) = fzz(i,j,k) - Gamxzz(i,j,k)*Lapx(i,j,k) - Gamyzz(i,j,k)*Lapy(i,j,k) - Gamzzz(i,j,k)*Lapz(i,j,k) + fxy(i,j,k) = fxy(i,j,k) - Gamxxy(i,j,k)*Lapx(i,j,k) - Gamyxy(i,j,k)*Lapy(i,j,k) - Gamzxy(i,j,k)*Lapz(i,j,k) + fxz(i,j,k) = fxz(i,j,k) - Gamxxz(i,j,k)*Lapx(i,j,k) - Gamyxz(i,j,k)*Lapy(i,j,k) - Gamzxz(i,j,k)*Lapz(i,j,k) + fyz(i,j,k) = fyz(i,j,k) - Gamxyz(i,j,k)*Lapx(i,j,k) - Gamyyz(i,j,k)*Lapy(i,j,k) - Gamzyz(i,j,k)*Lapz(i,j,k) + + trK_rhs(i,j,k) = gupxx(i,j,k) * fxx(i,j,k) + gupyy(i,j,k) * fyy(i,j,k) + gupzz(i,j,k) * fzz(i,j,k) + & + TWO * (gupxy(i,j,k) * fxy(i,j,k) + gupxz(i,j,k) * fxz(i,j,k) + gupyz(i,j,k) * fyz(i,j,k)) + enddo + enddo + enddo + do k=1,ex(3) + do j=1,ex(2) + do i=1,ex(1) + divb_loc = div_beta(i,j,k) + chin_loc = chin1(i,j,k) + + S_loc = chin_loc * ( gupxx(i,j,k) * Sxx(i,j,k) + gupyy(i,j,k) * Syy(i,j,k) + gupzz(i,j,k) * Szz(i,j,k) + & + TWO * (gupxy(i,j,k) * Sxy(i,j,k) + gupxz(i,j,k) * Sxz(i,j,k) + gupyz(i,j,k) * Syz(i,j,k)) ) + S(i,j,k) = S_loc + + f_loc = F2o3 * trK(i,j,k) * trK(i,j,k) - ( & + gupxx(i,j,k) * ( gupxx(i,j,k) * Axx(i,j,k) * Axx(i,j,k) + gupyy(i,j,k) * Axy(i,j,k) * Axy(i,j,k) + & + gupzz(i,j,k) * Axz(i,j,k) * Axz(i,j,k) + & + TWO * (gupxy(i,j,k) * Axx(i,j,k) * Axy(i,j,k) + gupxz(i,j,k) * Axx(i,j,k) * Axz(i,j,k) + & + gupyz(i,j,k) * Axy(i,j,k) * Axz(i,j,k)) ) + & + gupyy(i,j,k) * ( gupxx(i,j,k) * Axy(i,j,k) * Axy(i,j,k) + gupyy(i,j,k) * Ayy(i,j,k) * Ayy(i,j,k) + & + gupzz(i,j,k) * Ayz(i,j,k) * Ayz(i,j,k) + & + TWO * (gupxy(i,j,k) * Axy(i,j,k) * Ayy(i,j,k) + gupxz(i,j,k) * Axy(i,j,k) * Ayz(i,j,k) + & + gupyz(i,j,k) * Ayy(i,j,k) * Ayz(i,j,k)) ) + & + gupzz(i,j,k) * ( gupxx(i,j,k) * Axz(i,j,k) * Axz(i,j,k) + gupyy(i,j,k) * Ayz(i,j,k) * Ayz(i,j,k) + & + gupzz(i,j,k) * Azz(i,j,k) * Azz(i,j,k) + & + TWO * (gupxy(i,j,k) * Axz(i,j,k) * Ayz(i,j,k) + gupxz(i,j,k) * Axz(i,j,k) * Azz(i,j,k) + & + gupyz(i,j,k) * Ayz(i,j,k) * Azz(i,j,k)) ) + & + TWO * ( gupxy(i,j,k) * ( gupxx(i,j,k) * Axx(i,j,k) * Axy(i,j,k) + gupyy(i,j,k) * Axy(i,j,k) * Ayy(i,j,k) + & + gupzz(i,j,k) * Axz(i,j,k) * Ayz(i,j,k) + & + gupxy(i,j,k) * (Axx(i,j,k) * Ayy(i,j,k) + Axy(i,j,k) * Axy(i,j,k)) + & + gupxz(i,j,k) * (Axx(i,j,k) * Ayz(i,j,k) + Axz(i,j,k) * Axy(i,j,k)) + & + gupyz(i,j,k) * (Axy(i,j,k) * Ayz(i,j,k) + Axz(i,j,k) * Ayy(i,j,k)) ) + & + gupxz(i,j,k) * ( gupxx(i,j,k) * Axx(i,j,k) * Axz(i,j,k) + gupyy(i,j,k) * Axy(i,j,k) * Ayz(i,j,k) + & + gupzz(i,j,k) * Axz(i,j,k) * Azz(i,j,k) + & + gupxy(i,j,k) * (Axx(i,j,k) * Ayz(i,j,k) + Axy(i,j,k) * Axz(i,j,k)) + & + gupxz(i,j,k) * (Axx(i,j,k) * Azz(i,j,k) + Axz(i,j,k) * Axz(i,j,k)) + & + gupyz(i,j,k) * (Axy(i,j,k) * Azz(i,j,k) + Axz(i,j,k) * Ayz(i,j,k)) ) + & + gupyz(i,j,k) * ( gupxx(i,j,k) * Axy(i,j,k) * Axz(i,j,k) + gupyy(i,j,k) * Ayy(i,j,k) * Ayz(i,j,k) + & + gupzz(i,j,k) * Ayz(i,j,k) * Azz(i,j,k) + & + gupxy(i,j,k) * (Axy(i,j,k) * Ayz(i,j,k) + Ayy(i,j,k) * Axz(i,j,k)) + & + gupxz(i,j,k) * (Axy(i,j,k) * Azz(i,j,k) + Ayz(i,j,k) * Axz(i,j,k)) + & + gupyz(i,j,k) * (Ayy(i,j,k) * Azz(i,j,k) + Ayz(i,j,k) * Ayz(i,j,k)) ) ) ) - & + F16 * PI * rho(i,j,k) + EIGHT * PI * S_loc + + f_loc = -F1o3 * ( gupxx(i,j,k) * fxx(i,j,k) + gupyy(i,j,k) * fyy(i,j,k) + gupzz(i,j,k) * fzz(i,j,k) + & + TWO * (gupxy(i,j,k) * fxy(i,j,k) + gupxz(i,j,k) * fxz(i,j,k) + gupyz(i,j,k) * fyz(i,j,k)) + & + alpn1(i,j,k)/chin_loc * f_loc ) + f(i,j,k) = f_loc + + l_fxx = alpn1(i,j,k) * (Rxx(i,j,k) - EIGHT * PI * Sxx(i,j,k)) - fxx(i,j,k) + l_fxy = alpn1(i,j,k) * (Rxy(i,j,k) - EIGHT * PI * Sxy(i,j,k)) - fxy(i,j,k) + l_fxz = alpn1(i,j,k) * (Rxz(i,j,k) - EIGHT * PI * Sxz(i,j,k)) - fxz(i,j,k) + l_fyy = alpn1(i,j,k) * (Ryy(i,j,k) - EIGHT * PI * Syy(i,j,k)) - fyy(i,j,k) + l_fyz = alpn1(i,j,k) * (Ryz(i,j,k) - EIGHT * PI * Syz(i,j,k)) - fyz(i,j,k) + l_fzz = alpn1(i,j,k) * (Rzz(i,j,k) - EIGHT * PI * Szz(i,j,k)) - fzz(i,j,k) + + Axx_rhs(i,j,k) = l_fxx - gxx(i,j,k) * f_loc + Ayy_rhs(i,j,k) = l_fyy - gyy(i,j,k) * f_loc + Azz_rhs(i,j,k) = l_fzz - gzz(i,j,k) * f_loc + Axy_rhs(i,j,k) = l_fxy - gxy(i,j,k) * f_loc + Axz_rhs(i,j,k) = l_fxz - gxz(i,j,k) * f_loc + Ayz_rhs(i,j,k) = l_fyz - gyz(i,j,k) * f_loc + + fxx(i,j,k) = gupxx(i,j,k) * Axx(i,j,k) * Axx(i,j,k) + gupyy(i,j,k) * Axy(i,j,k) * Axy(i,j,k) + & + gupzz(i,j,k) * Axz(i,j,k) * Axz(i,j,k) + TWO * (gupxy(i,j,k) * Axx(i,j,k) * Axy(i,j,k) + & + gupxz(i,j,k) * Axx(i,j,k) * Axz(i,j,k) + gupyz(i,j,k) * Axy(i,j,k) * Axz(i,j,k)) + fyy(i,j,k) = gupxx(i,j,k) * Axy(i,j,k) * Axy(i,j,k) + gupyy(i,j,k) * Ayy(i,j,k) * Ayy(i,j,k) + & + gupzz(i,j,k) * Ayz(i,j,k) * Ayz(i,j,k) + TWO * (gupxy(i,j,k) * Axy(i,j,k) * Ayy(i,j,k) + & + gupxz(i,j,k) * Axy(i,j,k) * Ayz(i,j,k) + gupyz(i,j,k) * Ayy(i,j,k) * Ayz(i,j,k)) + fzz(i,j,k) = gupxx(i,j,k) * Axz(i,j,k) * Axz(i,j,k) + gupyy(i,j,k) * Ayz(i,j,k) * Ayz(i,j,k) + & + gupzz(i,j,k) * Azz(i,j,k) * Azz(i,j,k) + TWO * (gupxy(i,j,k) * Axz(i,j,k) * Ayz(i,j,k) + & + gupxz(i,j,k) * Axz(i,j,k) * Azz(i,j,k) + gupyz(i,j,k) * Ayz(i,j,k) * Azz(i,j,k)) + fxy(i,j,k) = gupxx(i,j,k) * Axx(i,j,k) * Axy(i,j,k) + gupyy(i,j,k) * Axy(i,j,k) * Ayy(i,j,k) + & + gupzz(i,j,k) * Axz(i,j,k) * Ayz(i,j,k) + gupxy(i,j,k) * (Axx(i,j,k) * Ayy(i,j,k) + Axy(i,j,k) * Axy(i,j,k)) + & + gupxz(i,j,k) * (Axx(i,j,k) * Ayz(i,j,k) + Axz(i,j,k) * Axy(i,j,k)) + & + gupyz(i,j,k) * (Axy(i,j,k) * Ayz(i,j,k) + Axz(i,j,k) * Ayy(i,j,k)) + fxz(i,j,k) = gupxx(i,j,k) * Axx(i,j,k) * Axz(i,j,k) + gupyy(i,j,k) * Axy(i,j,k) * Ayz(i,j,k) + & + gupzz(i,j,k) * Axz(i,j,k) * Azz(i,j,k) + gupxy(i,j,k) * (Axx(i,j,k) * Ayz(i,j,k) + Axy(i,j,k) * Axz(i,j,k)) + & + gupxz(i,j,k) * (Axx(i,j,k) * Azz(i,j,k) + Axz(i,j,k) * Axz(i,j,k)) + & + gupyz(i,j,k) * (Axy(i,j,k) * Azz(i,j,k) + Axz(i,j,k) * Ayz(i,j,k)) + fyz(i,j,k) = gupxx(i,j,k) * Axy(i,j,k) * Axz(i,j,k) + gupyy(i,j,k) * Ayy(i,j,k) * Ayz(i,j,k) + & + gupzz(i,j,k) * Ayz(i,j,k) * Azz(i,j,k) + gupxy(i,j,k) * (Axy(i,j,k) * Ayz(i,j,k) + Ayy(i,j,k) * Axz(i,j,k)) + & + gupxz(i,j,k) * (Axy(i,j,k) * Azz(i,j,k) + Ayz(i,j,k) * Axz(i,j,k)) + & + gupyz(i,j,k) * (Ayy(i,j,k) * Azz(i,j,k) + Ayz(i,j,k) * Ayz(i,j,k)) + + trK_rhs(i,j,k) = chin_loc * trK_rhs(i,j,k) + + Axx_rhs(i,j,k) = chin_loc * Axx_rhs(i,j,k) + alpn1(i,j,k) * (trK(i,j,k) * Axx(i,j,k) - TWO * fxx(i,j,k)) + & + TWO * (Axx(i,j,k) * betaxx(i,j,k) + Axy(i,j,k) * betayx(i,j,k) + Axz(i,j,k) * betazx(i,j,k)) - & + F2o3 * Axx(i,j,k) * divb_loc + Ayy_rhs(i,j,k) = chin_loc * Ayy_rhs(i,j,k) + alpn1(i,j,k) * (trK(i,j,k) * Ayy(i,j,k) - TWO * fyy(i,j,k)) + & + TWO * (Axy(i,j,k) * betaxy(i,j,k) + Ayy(i,j,k) * betayy(i,j,k) + Ayz(i,j,k) * betazy(i,j,k)) - & + F2o3 * Ayy(i,j,k) * divb_loc + Azz_rhs(i,j,k) = chin_loc * Azz_rhs(i,j,k) + alpn1(i,j,k) * (trK(i,j,k) * Azz(i,j,k) - TWO * fzz(i,j,k)) + & + TWO * (Axz(i,j,k) * betaxz(i,j,k) + Ayz(i,j,k) * betayz(i,j,k) + Azz(i,j,k) * betazz(i,j,k)) - & + F2o3 * Azz(i,j,k) * divb_loc + Axy_rhs(i,j,k) = chin_loc * Axy_rhs(i,j,k) + alpn1(i,j,k) * (trK(i,j,k) * Axy(i,j,k) - TWO * fxy(i,j,k)) + & + Axx(i,j,k) * betaxy(i,j,k) + Axz(i,j,k) * betazy(i,j,k) + Ayy(i,j,k) * betayx(i,j,k) + & + Ayz(i,j,k) * betazx(i,j,k) + F1o3 * Axy(i,j,k) * divb_loc - Axy(i,j,k) * betazz(i,j,k) + Ayz_rhs(i,j,k) = chin_loc * Ayz_rhs(i,j,k) + alpn1(i,j,k) * (trK(i,j,k) * Ayz(i,j,k) - TWO * fyz(i,j,k)) + & + Axy(i,j,k) * betaxz(i,j,k) + Ayy(i,j,k) * betayz(i,j,k) + Axz(i,j,k) * betaxy(i,j,k) + & + Azz(i,j,k) * betazy(i,j,k) + F1o3 * Ayz(i,j,k) * divb_loc - Ayz(i,j,k) * betaxx(i,j,k) + Axz_rhs(i,j,k) = chin_loc * Axz_rhs(i,j,k) + alpn1(i,j,k) * (trK(i,j,k) * Axz(i,j,k) - TWO * fxz(i,j,k)) + & + Axx(i,j,k) * betaxz(i,j,k) + Axy(i,j,k) * betayz(i,j,k) + Ayz(i,j,k) * betayx(i,j,k) + & + Azz(i,j,k) * betazx(i,j,k) + F1o3 * Axz(i,j,k) * divb_loc - Axz(i,j,k) * betayy(i,j,k) + + trK_rhs(i,j,k) = - trK_rhs(i,j,k) + alpn1(i,j,k) * ( F1o3 * trK(i,j,k) * trK(i,j,k) + & + gupxx(i,j,k) * fxx(i,j,k) + gupyy(i,j,k) * fyy(i,j,k) + gupzz(i,j,k) * fzz(i,j,k) + & + TWO * (gupxy(i,j,k) * fxy(i,j,k) + gupxz(i,j,k) * fxz(i,j,k) + gupyz(i,j,k) * fyz(i,j,k)) + & + FOUR * PI * (rho(i,j,k) + S_loc) ) + enddo + enddo + enddo !!!! gauge variable part