fix(bssn_rhs)

OutBdLow2Hi_cached 读的是 coarse owned 区域（非 coarse ghost/buffer）
回退旧行为：编译时定义 RP_SYNC_COARSE_AFTER_RESTRICT=1

AMSS_NCKU_source/bssn_class.C

@@ -41,12 +41,6 @@ using namespace std;
 #include "derivatives.h"
 #include "ricci_gamma.h"
 
-// Compile-time switch for per-timestep memory usage collection/printing.
-// Default is OFF to reduce overhead in production runs.
-#ifndef BSSN_ENABLE_MEM_USAGE_LOG
-#define BSSN_ENABLE_MEM_USAGE_LOG 0
-#endif
-
 //================================================================================================
 
 // define bssn_class
@@ -2135,10 +2129,8 @@ void bssn_class::Evolve(int Steps)
   #endif
   */
   
-#if BSSN_ENABLE_MEM_USAGE_LOG
   perf bssn_perf;
   size_t current_min, current_avg, current_max, peak_min, peak_avg, peak_max;
-#endif
 
   for (int lev = 0; lev < GH->levels; lev++)
     GH->Lt[lev] = PhysTime;
@@ -2232,7 +2224,6 @@ void bssn_class::Evolve(int Steps)
 //		fgt(PhysTime-dT_mon,StartTime,dT_mon/2),ErrorMonitor);
 #endif
 
-#if BSSN_ENABLE_MEM_USAGE_LOG
     // Retrieve memory usage information used during computation; master process prints it
     bssn_perf.MemoryUsage(&current_min, &current_avg, &current_max,
                           &peak_min, &peak_avg, &peak_max, nprocs);
@@ -2248,7 +2239,6 @@ void bssn_class::Evolve(int Steps)
              (double)peak_max / (1024.0 * 1024.0));
       cout << endl;
     }
-#endif
     
     // Output puncture positions at each step
     if (myrank == 0)
@@ -5758,6 +5748,12 @@ void bssn_class::SHStep()
 // 0: do not use mixing two levels data for OutBD; 1: do use
 
 #define MIXOUTB 0
+// In the cached Restrict->OutBdLow2Hi path, coarse Sync is usually redundant:
+// OutBdLow2Hi_cached reads coarse owned cells (build_owned_gsl type-4), not coarse ghost/buffer cells.
+// Keep a switch to restore the old behavior if needed for debugging.
+#ifndef RP_SYNC_COARSE_AFTER_RESTRICT
+#define RP_SYNC_COARSE_AFTER_RESTRICT 0
+#endif
 void bssn_class::RestrictProlong(int lev, int YN, bool BB,
                                  MyList<var> *SL, MyList<var> *OL, MyList<var> *corL)
 // we assume
@@ -5821,7 +5817,9 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB,
 //       misc::tillherecheck(GH->Commlev[GH->mylev],GH->start_rank[GH->mylev],a_stream.str());
 #endif
 
+#if (RP_SYNC_COARSE_AFTER_RESTRICT == 1)
       Parallel::Sync_cached(GH->PatL[lev - 1], SynchList_pre, Symmetry, sync_cache_rp_coarse[lev]);
+#endif
 
 #if (PSTR == 1 || PSTR == 2)
 //       a_stream.clear();
@@ -5872,7 +5870,9 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB,
 //       misc::tillherecheck(GH->Commlev[GH->mylev],GH->start_rank[GH->mylev],a_stream.str());
 #endif
 
+#if (RP_SYNC_COARSE_AFTER_RESTRICT == 1)
       Parallel::Sync_cached(GH->PatL[lev - 1], SL, Symmetry, sync_cache_rp_coarse[lev]);
+#endif
 
 #if (PSTR == 1 || PSTR == 2)
 //       a_stream.clear();
@@ -5958,7 +5958,9 @@ void bssn_class::RestrictProlong_aux(int lev, int YN, bool BB,
       Parallel::Restrict_bam(GH->PatL[lev - 1], GH->PatL[lev], SL, SynchList_pre, GH->rsul[lev], Symmetry);
 #endif
 
+#if (RP_SYNC_COARSE_AFTER_RESTRICT == 1)
       Parallel::Sync_cached(GH->PatL[lev - 1], SynchList_pre, Symmetry, sync_cache_rp_coarse[lev]);
+#endif
 
 #if (RPB == 0)
 #if (MIXOUTB == 0)
@@ -5980,7 +5982,9 @@ void bssn_class::RestrictProlong_aux(int lev, int YN, bool BB,
       Parallel::Restrict_bam(GH->PatL[lev - 1], GH->PatL[lev], SL, SL, GH->rsul[lev], Symmetry);
 #endif
 
+#if (RP_SYNC_COARSE_AFTER_RESTRICT == 1)
       Parallel::Sync_cached(GH->PatL[lev - 1], SL, Symmetry, sync_cache_rp_coarse[lev]);
+#endif
 
 #if (RPB == 0)
 #if (MIXOUTB == 0)
@@ -6045,7 +6049,9 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB)
       Parallel::Restrict_bam(GH->PatL[lev - 1], GH->PatL[lev], SynchList_cor, SynchList_pre, GH->rsul[lev], Symmetry);
 #endif
 
+#if (RP_SYNC_COARSE_AFTER_RESTRICT == 1)
       Parallel::Sync_cached(GH->PatL[lev - 1], SynchList_pre, Symmetry, sync_cache_rp_coarse[lev]);
+#endif
 
 #if (RPB == 0)
 #if (MIXOUTB == 0)
@@ -6069,7 +6075,9 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB)
       Parallel::Restrict_bam(GH->PatL[lev - 1], GH->PatL[lev], SynchList_cor, StateList, GH->rsul[lev], Symmetry);
 #endif
 
+#if (RP_SYNC_COARSE_AFTER_RESTRICT == 1)
       Parallel::Sync_cached(GH->PatL[lev - 1], StateList, Symmetry, sync_cache_rp_coarse[lev]);
+#endif
 
 #if (RPB == 0)
 #if (MIXOUTB == 0)
@@ -6146,13 +6154,16 @@ void bssn_class::ProlongRestrict(int lev, int YN, bool BB)
 #else
       Parallel::Restrict_after(GH->PatL[lev - 1], GH->PatL[lev], SynchList_cor, StateList, Symmetry);
 #endif
+#if (RP_SYNC_COARSE_AFTER_RESTRICT == 1)
       Parallel::Sync_cached(GH->PatL[lev - 1], StateList, Symmetry, sync_cache_rp_coarse[lev]);
+#endif
     }
 
     Parallel::Sync_cached(GH->PatL[lev], SynchList_cor, Symmetry, sync_cache_rp_fine[lev]);
   }
 }
 #undef MIXOUTB
+#undef RP_SYNC_COARSE_AFTER_RESTRICT
 
 //================================================================================================
 

AMSS_NCKU_source/bssn_rhs.f90

@@ -62,7 +62,6 @@
   real*8, dimension(ex(1),ex(2),ex(3)),intent(inout) :: Gmx_Res, Gmy_Res, Gmz_Res
 !  gont = 0: success; gont = 1: something wrong
   integer::gont
-  integer :: i,j,k
 
 !~~~~~~> Other variables:
 
@@ -86,13 +85,6 @@
 
   real*8,dimension(3) ::SSS,AAS,ASA,SAA,ASS,SAS,SSA
   real*8            :: dX, dY, dZ, PI
-  real*8            :: divb_loc,det_loc
-  real*8            :: gupxx_loc,gupxy_loc,gupxz_loc,gupyy_loc,gupyz_loc,gupzz_loc
-  real*8            :: Rxx_loc,Rxy_loc,Rxz_loc,Ryy_loc,Ryz_loc,Rzz_loc
-  real*8            :: fxx_loc,fxy_loc,fxz_loc
-  real*8            :: Gamxa_loc,Gamya_loc,Gamza_loc
-  real*8            :: 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
@@ -105,7 +97,7 @@
 #endif
 
 #if (GAUGE == 6 || GAUGE == 7)
-  integer :: BHN
+  integer :: BHN,i,j,k
   real*8, dimension(9) :: Porg
   real*8, dimension(3) :: Mass
   real*8 :: r1,r2,M,A,w1,w2,C1,C2
@@ -153,24 +145,22 @@
   dY = Y(2) - Y(1)
   dZ = Z(2) - Z(1)
 
-  do k=1,ex(3)
+  alpn1 = Lap + ONE
-  do j=1,ex(2)
+  chin1 = chi + ONE
-  do i=1,ex(1)
+  gxx = dxx + ONE
-    alpn1(i,j,k) = Lap(i,j,k) + ONE
+  gyy = dyy + ONE
-    chin1(i,j,k) = chi(i,j,k) + ONE
+  gzz = dzz + ONE
-    gxx(i,j,k) = dxx(i,j,k) + ONE
-    gyy(i,j,k) = dyy(i,j,k) + ONE
-    gzz(i,j,k) = dzz(i,j,k) + ONE
-  enddo
-  enddo
-  enddo
 
   call fderivs(ex,betax,betaxx,betaxy,betaxz,X,Y,Z,ANTI, SYM, SYM,Symmetry,Lev)
   call fderivs(ex,betay,betayx,betayy,betayz,X,Y,Z, SYM,ANTI, SYM,Symmetry,Lev)
   call fderivs(ex,betaz,betazx,betazy,betazz,X,Y,Z, SYM, SYM,ANTI,Symmetry,Lev)
   
+  div_beta = betaxx + betayy + betazz
+ 
   call fderivs(ex,chi,chix,chiy,chiz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev)
 
+  chi_rhs = F2o3 *chin1*( alpn1 * trK - div_beta ) !rhs for chi
+
   call fderivs(ex,dxx,gxxx,gxxy,gxxz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev)
   call fderivs(ex,gxy,gxyx,gxyy,gxyz,X,Y,Z,ANTI,ANTI,SYM ,Symmetry,Lev)
   call fderivs(ex,gxz,gxzx,gxzy,gxzz,X,Y,Z,ANTI,SYM ,ANTI,Symmetry,Lev)
@@ -178,179 +168,151 @@
   call fderivs(ex,gyz,gyzx,gyzy,gyzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,Lev)
   call fderivs(ex,dzz,gzzx,gzzy,gzzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,Lev)
 
-  do k=1,ex(3)
+  gxx_rhs = - TWO * alpn1 * Axx    -  F2o3 * gxx * div_beta          + &
-  do j=1,ex(2)
+              TWO *(  gxx * betaxx +   gxy * betayx +   gxz * betazx)
-  do i=1,ex(1)
-    divb_loc = betaxx(i,j,k) + betayy(i,j,k) + betazz(i,j,k)
-    div_beta(i,j,k) = divb_loc
 
-    chi_rhs(i,j,k) = F2o3 * chin1(i,j,k) * (alpn1(i,j,k) * trK(i,j,k) - divb_loc)
+  gyy_rhs = - TWO * alpn1 * Ayy    -  F2o3 * gyy * div_beta          + &
+              TWO *(  gxy * betaxy +   gyy * betayy +   gyz * betazy)
 
-    gxx_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Axx(i,j,k) - F2o3 * gxx(i,j,k) * divb_loc + &
+  gzz_rhs = - TWO * alpn1 * Azz    -  F2o3 * gzz * div_beta          + &
-         TWO * ( gxx(i,j,k) * betaxx(i,j,k) + gxy(i,j,k) * betayx(i,j,k) + gxz(i,j,k) * betazx(i,j,k) )
+              TWO *(  gxz * betaxz +   gyz * betayz +   gzz * betazz)
 
-    gyy_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Ayy(i,j,k) - F2o3 * gyy(i,j,k) * divb_loc + &
+  gxy_rhs = - TWO * alpn1 * Axy    +  F1o3 * gxy    * div_beta       + &
-         TWO * ( gxy(i,j,k) * betaxy(i,j,k) + gyy(i,j,k) * betayy(i,j,k) + gyz(i,j,k) * betazy(i,j,k) )
+                      gxx * betaxy                  +   gxz * betazy + &
+                                       gyy * betayx +   gyz * betazx   &
+                                                    -   gxy * betazz
 
-    gzz_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Azz(i,j,k) - F2o3 * gzz(i,j,k) * divb_loc + &
+  gyz_rhs = - TWO * alpn1 * Ayz    +  F1o3 * gyz    * div_beta       + &
-         TWO * ( gxz(i,j,k) * betaxz(i,j,k) + gyz(i,j,k) * betayz(i,j,k) + gzz(i,j,k) * betazz(i,j,k) )
+                      gxy * betaxz +   gyy * betayz                  + &
+                      gxz * betaxy                  +   gzz * betazy   &
+                                                    -   gyz * betaxx
  
-    gxy_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Axy(i,j,k) + F1o3 * gxy(i,j,k) * divb_loc + &
+  gxz_rhs = - TWO * alpn1 * Axz    +  F1o3 * gxz    * div_beta       + &
-         gxx(i,j,k) * betaxy(i,j,k) + gxz(i,j,k) * betazy(i,j,k) + gyy(i,j,k) * betayx(i,j,k) + &
+                      gxx * betaxz +   gxy * betayz                  + &
-         gyz(i,j,k) * betazx(i,j,k) - gxy(i,j,k) * betazz(i,j,k)
+                                       gyz * betayx +   gzz * betazx   &
+                                                    -   gxz * betayy     !rhs for gij
 
-    gyz_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Ayz(i,j,k) + F1o3 * gyz(i,j,k) * divb_loc + &
+! invert tilted metric
-         gxy(i,j,k) * betaxz(i,j,k) + gyy(i,j,k) * betayz(i,j,k) + gxz(i,j,k) * betaxy(i,j,k) + &
+  gupzz =  gxx * gyy * gzz + gxy * gyz * gxz + gxz * gxy * gyz - &
-         gzz(i,j,k) * betazy(i,j,k) - gyz(i,j,k) * betaxx(i,j,k)
+           gxz * gyy * gxz - gxy * gxy * gzz - gxx * gyz * gyz
-
+  gupxx =   ( gyy * gzz - gyz * gyz ) / gupzz
-    gxz_rhs(i,j,k) = - TWO * alpn1(i,j,k) * Axz(i,j,k) + F1o3 * gxz(i,j,k) * divb_loc + &
+  gupxy = - ( gxy * gzz - gyz * gxz ) / gupzz
-         gxx(i,j,k) * betaxz(i,j,k) + gxy(i,j,k) * betayz(i,j,k) + gyz(i,j,k) * betayx(i,j,k) + &
+  gupxz =   ( gxy * gyz - gyy * gxz ) / gupzz
-         gzz(i,j,k) * betazx(i,j,k) - gxz(i,j,k) * betayy(i,j,k)
+  gupyy =   ( gxx * gzz - gxz * gxz ) / gupzz
-
+  gupyz = - ( gxx * gyz - gxy * gxz ) / gupzz
-    det_loc = gxx(i,j,k) * gyy(i,j,k) * gzz(i,j,k) + gxy(i,j,k) * gyz(i,j,k) * gxz(i,j,k) + &
+  gupzz =   ( gxx * gyy - gxy * gxy ) / gupzz
-         gxz(i,j,k) * gxy(i,j,k) * gyz(i,j,k) - gxz(i,j,k) * gyy(i,j,k) * gxz(i,j,k) - &
-         gxy(i,j,k) * gxy(i,j,k) * gzz(i,j,k) - gxx(i,j,k) * gyz(i,j,k) * gyz(i,j,k)
-    gupxx_loc = ( gyy(i,j,k) * gzz(i,j,k) - gyz(i,j,k) * gyz(i,j,k) ) / det_loc
-    gupxy_loc = - ( gxy(i,j,k) * gzz(i,j,k) - gyz(i,j,k) * gxz(i,j,k) ) / det_loc
-    gupxz_loc = ( gxy(i,j,k) * gyz(i,j,k) - gyy(i,j,k) * gxz(i,j,k) ) / det_loc
-    gupyy_loc = ( gxx(i,j,k) * gzz(i,j,k) - gxz(i,j,k) * gxz(i,j,k) ) / det_loc
-    gupyz_loc = - ( gxx(i,j,k) * gyz(i,j,k) - gxy(i,j,k) * gxz(i,j,k) ) / det_loc
-    gupzz_loc = ( gxx(i,j,k) * gyy(i,j,k) - gxy(i,j,k) * gxy(i,j,k) ) / det_loc
-    gupxx(i,j,k) = gupxx_loc
-    gupxy(i,j,k) = gupxy_loc
-    gupxz(i,j,k) = gupxz_loc
-    gupyy(i,j,k) = gupyy_loc
-    gupyz(i,j,k) = gupyz_loc
-    gupzz(i,j,k) = gupzz_loc
 
   if(co == 0)then
-      Gmx_Res(i,j,k) = Gamx(i,j,k) - ( &
+! Gam^i_Res = Gam^i + gup^ij_,j
-           gupxx_loc*(gupxx_loc*gxxx(i,j,k)+gupxy_loc*gxyx(i,j,k)+gupxz_loc*gxzx(i,j,k)) + &
+  Gmx_Res = Gamx - (gupxx*(gupxx*gxxx+gupxy*gxyx+gupxz*gxzx)&
-           gupxy_loc*(gupxx_loc*gxyx(i,j,k)+gupxy_loc*gyyx(i,j,k)+gupxz_loc*gyzx(i,j,k)) + &
+                   +gupxy*(gupxx*gxyx+gupxy*gyyx+gupxz*gyzx)&
-           gupxz_loc*(gupxx_loc*gxzx(i,j,k)+gupxy_loc*gyzx(i,j,k)+gupxz_loc*gzzx(i,j,k)) + &
+                   +gupxz*(gupxx*gxzx+gupxy*gyzx+gupxz*gzzx)&
-           gupxx_loc*(gupxy_loc*gxxy(i,j,k)+gupyy_loc*gxyy(i,j,k)+gupyz_loc*gxzy(i,j,k)) + &
+                   +gupxx*(gupxy*gxxy+gupyy*gxyy+gupyz*gxzy)&
-           gupxy_loc*(gupxy_loc*gxyy(i,j,k)+gupyy_loc*gyyy(i,j,k)+gupyz_loc*gyzy(i,j,k)) + &
+                   +gupxy*(gupxy*gxyy+gupyy*gyyy+gupyz*gyzy)&
-           gupxz_loc*(gupxy_loc*gxzy(i,j,k)+gupyy_loc*gyzy(i,j,k)+gupyz_loc*gzzy(i,j,k)) + &
+                   +gupxz*(gupxy*gxzy+gupyy*gyzy+gupyz*gzzy)&
-           gupxx_loc*(gupxz_loc*gxxz(i,j,k)+gupyz_loc*gxyz(i,j,k)+gupzz_loc*gxzz(i,j,k)) + &
+                   +gupxx*(gupxz*gxxz+gupyz*gxyz+gupzz*gxzz)&
-           gupxy_loc*(gupxz_loc*gxyz(i,j,k)+gupyz_loc*gyyz(i,j,k)+gupzz_loc*gyzz(i,j,k)) + &
+                   +gupxy*(gupxz*gxyz+gupyz*gyyz+gupzz*gyzz)&
-           gupxz_loc*(gupxz_loc*gxzz(i,j,k)+gupyz_loc*gyzz(i,j,k)+gupzz_loc*gzzz(i,j,k)))
+                   +gupxz*(gupxz*gxzz+gupyz*gyzz+gupzz*gzzz))
-      Gmy_Res(i,j,k) = Gamy(i,j,k) - ( &
+  Gmy_Res = Gamy - (gupxx*(gupxy*gxxx+gupyy*gxyx+gupyz*gxzx)&
-           gupxx_loc*(gupxy_loc*gxxx(i,j,k)+gupyy_loc*gxyx(i,j,k)+gupyz_loc*gxzx(i,j,k)) + &
+                   +gupxy*(gupxy*gxyx+gupyy*gyyx+gupyz*gyzx)&
-           gupxy_loc*(gupxy_loc*gxyx(i,j,k)+gupyy_loc*gyyx(i,j,k)+gupyz_loc*gyzx(i,j,k)) + &
+                   +gupxz*(gupxy*gxzx+gupyy*gyzx+gupyz*gzzx)&
-           gupxz_loc*(gupxy_loc*gxzx(i,j,k)+gupyy_loc*gyzx(i,j,k)+gupyz_loc*gzzx(i,j,k)) + &
+                   +gupxy*(gupxy*gxxy+gupyy*gxyy+gupyz*gxzy)&
-           gupxy_loc*(gupxy_loc*gxxy(i,j,k)+gupyy_loc*gxyy(i,j,k)+gupyz_loc*gxzy(i,j,k)) + &
+                   +gupyy*(gupxy*gxyy+gupyy*gyyy+gupyz*gyzy)&
-           gupyy_loc*(gupxy_loc*gxyy(i,j,k)+gupyy_loc*gyyy(i,j,k)+gupyz_loc*gyzy(i,j,k)) + &
+                   +gupyz*(gupxy*gxzy+gupyy*gyzy+gupyz*gzzy)&
-           gupyz_loc*(gupxy_loc*gxzy(i,j,k)+gupyy_loc*gyzy(i,j,k)+gupyz_loc*gzzy(i,j,k)) + &
+                   +gupxy*(gupxz*gxxz+gupyz*gxyz+gupzz*gxzz)&
-           gupxy_loc*(gupxz_loc*gxxz(i,j,k)+gupyz_loc*gxyz(i,j,k)+gupzz_loc*gxzz(i,j,k)) + &
+                   +gupyy*(gupxz*gxyz+gupyz*gyyz+gupzz*gyzz)&
-           gupyy_loc*(gupxz_loc*gxyz(i,j,k)+gupyz_loc*gyyz(i,j,k)+gupzz_loc*gyzz(i,j,k)) + &
+                   +gupyz*(gupxz*gxzz+gupyz*gyzz+gupzz*gzzz))
-           gupyz_loc*(gupxz_loc*gxzz(i,j,k)+gupyz_loc*gyzz(i,j,k)+gupzz_loc*gzzz(i,j,k)))
+  Gmz_Res = Gamz - (gupxx*(gupxz*gxxx+gupyz*gxyx+gupzz*gxzx)&
-      Gmz_Res(i,j,k) = Gamz(i,j,k) - ( &
+                   +gupxy*(gupxz*gxyx+gupyz*gyyx+gupzz*gyzx)&
-           gupxx_loc*(gupxz_loc*gxxx(i,j,k)+gupyz_loc*gxyx(i,j,k)+gupzz_loc*gxzx(i,j,k)) + &
+                   +gupxz*(gupxz*gxzx+gupyz*gyzx+gupzz*gzzx)&
-           gupxy_loc*(gupxz_loc*gxyx(i,j,k)+gupyz_loc*gyyx(i,j,k)+gupzz_loc*gyzx(i,j,k)) + &
+                   +gupxy*(gupxz*gxxy+gupyz*gxyy+gupzz*gxzy)&
-           gupxz_loc*(gupxz_loc*gxzx(i,j,k)+gupyz_loc*gyzx(i,j,k)+gupzz_loc*gzzx(i,j,k)) + &
+                   +gupyy*(gupxz*gxyy+gupyz*gyyy+gupzz*gyzy)&
-           gupxy_loc*(gupxz_loc*gxxy(i,j,k)+gupyz_loc*gxyy(i,j,k)+gupzz_loc*gxzy(i,j,k)) + &
+                   +gupyz*(gupxz*gxzy+gupyz*gyzy+gupzz*gzzy)&
-           gupyy_loc*(gupxz_loc*gxyy(i,j,k)+gupyz_loc*gyyy(i,j,k)+gupzz_loc*gyzy(i,j,k)) + &
+                   +gupxz*(gupxz*gxxz+gupyz*gxyz+gupzz*gxzz)&
-           gupyz_loc*(gupxz_loc*gxzy(i,j,k)+gupyz_loc*gyzy(i,j,k)+gupzz_loc*gzzy(i,j,k)) + &
+                   +gupyz*(gupxz*gxyz+gupyz*gyyz+gupzz*gyzz)&
-           gupxz_loc*(gupxz_loc*gxxz(i,j,k)+gupyz_loc*gxyz(i,j,k)+gupzz_loc*gxzz(i,j,k)) + &
+                   +gupzz*(gupxz*gxzz+gupyz*gyzz+gupzz*gzzz))
-           gupyz_loc*(gupxz_loc*gxyz(i,j,k)+gupyz_loc*gyyz(i,j,k)+gupzz_loc*gyzz(i,j,k)) + &
-           gupzz_loc*(gupxz_loc*gxzz(i,j,k)+gupyz_loc*gyzz(i,j,k)+gupzz_loc*gzzz(i,j,k)))
   endif
 
-    Gamxxx(i,j,k)=HALF*( gupxx_loc*gxxx(i,j,k) + gupxy_loc*(TWO*gxyx(i,j,k) - gxxy(i,j,k)) + gupxz_loc*(TWO*gxzx(i,j,k) - gxxz(i,j,k)))
+! second kind of connection
-    Gamyxx(i,j,k)=HALF*( gupxy_loc*gxxx(i,j,k) + gupyy_loc*(TWO*gxyx(i,j,k) - gxxy(i,j,k)) + gupyz_loc*(TWO*gxzx(i,j,k) - gxxz(i,j,k)))
+  Gamxxx =HALF*( gupxx*gxxx + gupxy*(TWO*gxyx - gxxy ) + gupxz*(TWO*gxzx - gxxz ))
-    Gamzxx(i,j,k)=HALF*( gupxz_loc*gxxx(i,j,k) + gupyz_loc*(TWO*gxyx(i,j,k) - gxxy(i,j,k)) + gupzz_loc*(TWO*gxzx(i,j,k) - gxxz(i,j,k)))
+  Gamyxx =HALF*( gupxy*gxxx + gupyy*(TWO*gxyx - gxxy ) + gupyz*(TWO*gxzx - gxxz ))
+  Gamzxx =HALF*( gupxz*gxxx + gupyz*(TWO*gxyx - gxxy ) + gupzz*(TWO*gxzx - gxxz ))
  
-    Gamxyy(i,j,k)=HALF*( gupxx_loc*(TWO*gxyy(i,j,k) - gyyx(i,j,k)) + gupxy_loc*gyyy(i,j,k) + gupxz_loc*(TWO*gyzy(i,j,k) - gyyz(i,j,k)))
+  Gamxyy =HALF*( gupxx*(TWO*gxyy - gyyx ) + gupxy*gyyy + gupxz*(TWO*gyzy - gyyz ))
-    Gamyyy(i,j,k)=HALF*( gupxy_loc*(TWO*gxyy(i,j,k) - gyyx(i,j,k)) + gupyy_loc*gyyy(i,j,k) + gupyz_loc*(TWO*gyzy(i,j,k) - gyyz(i,j,k)))
+  Gamyyy =HALF*( gupxy*(TWO*gxyy - gyyx ) + gupyy*gyyy + gupyz*(TWO*gyzy - gyyz ))
-    Gamzyy(i,j,k)=HALF*( gupxz_loc*(TWO*gxyy(i,j,k) - gyyx(i,j,k)) + gupyz_loc*gyyy(i,j,k) + gupzz_loc*(TWO*gyzy(i,j,k) - gyyz(i,j,k)))
+  Gamzyy =HALF*( gupxz*(TWO*gxyy - gyyx ) + gupyz*gyyy + gupzz*(TWO*gyzy - gyyz ))
 
-    Gamxzz(i,j,k)=HALF*( gupxx_loc*(TWO*gxzz(i,j,k) - gzzx(i,j,k)) + gupxy_loc*(TWO*gyzz(i,j,k) - gzzy(i,j,k)) + gupxz_loc*gzzz(i,j,k))
+  Gamxzz =HALF*( gupxx*(TWO*gxzz - gzzx ) + gupxy*(TWO*gyzz - gzzy ) + gupxz*gzzz)
-    Gamyzz(i,j,k)=HALF*( gupxy_loc*(TWO*gxzz(i,j,k) - gzzx(i,j,k)) + gupyy_loc*(TWO*gyzz(i,j,k) - gzzy(i,j,k)) + gupyz_loc*gzzz(i,j,k))
+  Gamyzz =HALF*( gupxy*(TWO*gxzz - gzzx ) + gupyy*(TWO*gyzz - gzzy ) + gupyz*gzzz)
-    Gamzzz(i,j,k)=HALF*( gupxz_loc*(TWO*gxzz(i,j,k) - gzzx(i,j,k)) + gupyz_loc*(TWO*gyzz(i,j,k) - gzzy(i,j,k)) + gupzz_loc*gzzz(i,j,k))
+  Gamzzz =HALF*( gupxz*(TWO*gxzz - gzzx ) + gupyz*(TWO*gyzz - gzzy ) + gupzz*gzzz)
 
-    Gamxxy(i,j,k)=HALF*( gupxx_loc*gxxy(i,j,k) + gupxy_loc*gyyx(i,j,k) + gupxz_loc*(gxzy(i,j,k) + gyzx(i,j,k) - gxyz(i,j,k)) )
+  Gamxxy =HALF*( gupxx*gxxy + gupxy*gyyx + gupxz*( gxzy + gyzx - gxyz ) )
-    Gamyxy(i,j,k)=HALF*( gupxy_loc*gxxy(i,j,k) + gupyy_loc*gyyx(i,j,k) + gupyz_loc*(gxzy(i,j,k) + gyzx(i,j,k) - gxyz(i,j,k)) )
+  Gamyxy =HALF*( gupxy*gxxy + gupyy*gyyx + gupyz*( gxzy + gyzx - gxyz ) )
-    Gamzxy(i,j,k)=HALF*( gupxz_loc*gxxy(i,j,k) + gupyz_loc*gyyx(i,j,k) + gupzz_loc*(gxzy(i,j,k) + gyzx(i,j,k) - gxyz(i,j,k)) )
+  Gamzxy =HALF*( gupxz*gxxy + gupyz*gyyx + gupzz*( gxzy + gyzx - gxyz ) )
 
-    Gamxxz(i,j,k)=HALF*( gupxx_loc*gxxz(i,j,k) + gupxy_loc*(gxyz(i,j,k) + gyzx(i,j,k) - gxzy(i,j,k)) + gupxz_loc*gzzx(i,j,k) )
+  Gamxxz =HALF*( gupxx*gxxz + gupxy*( gxyz + gyzx - gxzy ) + gupxz*gzzx )
-    Gamyxz(i,j,k)=HALF*( gupxy_loc*gxxz(i,j,k) + gupyy_loc*(gxyz(i,j,k) + gyzx(i,j,k) - gxzy(i,j,k)) + gupyz_loc*gzzx(i,j,k) )
+  Gamyxz =HALF*( gupxy*gxxz + gupyy*( gxyz + gyzx - gxzy ) + gupyz*gzzx )
-    Gamzxz(i,j,k)=HALF*( gupxz_loc*gxxz(i,j,k) + gupyz_loc*(gxyz(i,j,k) + gyzx(i,j,k) - gxzy(i,j,k)) + gupzz_loc*gzzx(i,j,k) )
+  Gamzxz =HALF*( gupxz*gxxz + gupyz*( gxyz + gyzx - gxzy ) + gupzz*gzzx )
 
-    Gamxyz(i,j,k)=HALF*( gupxx_loc*(gxyz(i,j,k) + gxzy(i,j,k) - gyzx(i,j,k)) + gupxy_loc*gyyz(i,j,k) + gupxz_loc*gzzy(i,j,k) )
+  Gamxyz =HALF*( gupxx*( gxyz + gxzy - gyzx ) + gupxy*gyyz + gupxz*gzzy )
-    Gamyyz(i,j,k)=HALF*( gupxy_loc*(gxyz(i,j,k) + gxzy(i,j,k) - gyzx(i,j,k)) + gupyy_loc*gyyz(i,j,k) + gupyz_loc*gzzy(i,j,k) )
+  Gamyyz =HALF*( gupxy*( gxyz + gxzy - gyzx ) + gupyy*gyyz + gupyz*gzzy )
-    Gamzyz(i,j,k)=HALF*( gupxz_loc*(gxyz(i,j,k) + gxzy(i,j,k) - gyzx(i,j,k)) + gupyz_loc*gyyz(i,j,k) + gupzz_loc*gzzy(i,j,k) )
+  Gamzyz =HALF*( gupxz*( gxyz + gxzy - gyzx ) + gupyz*gyyz + gupzz*gzzy )
-  enddo
-  enddo
-  enddo
 ! Raise indices of \tilde A_{ij} and store in R_ij
 
+  Rxx =    gupxx * gupxx * Axx + gupxy * gupxy * Ayy + gupxz * gupxz * Azz + &
+      TWO*(gupxx * gupxy * Axy + gupxx * gupxz * Axz + gupxy * gupxz * Ayz)
+
+  Ryy =    gupxy * gupxy * Axx + gupyy * gupyy * Ayy + gupyz * gupyz * Azz + &
+      TWO*(gupxy * gupyy * Axy + gupxy * gupyz * Axz + gupyy * gupyz * Ayz)
+
+  Rzz =    gupxz * gupxz * Axx + gupyz * gupyz * Ayy + gupzz * gupzz * Azz + &
+      TWO*(gupxz * gupyz * Axy + gupxz * gupzz * Axz + gupyz * gupzz * Ayz)
+
+  Rxy =    gupxx * gupxy * Axx + gupxy * gupyy * Ayy + gupxz * gupyz * Azz + &
+          (gupxx * gupyy       + gupxy * gupxy)* Axy                       + &
+          (gupxx * gupyz       + gupxz * gupxy)* Axz                       + &
+          (gupxy * gupyz       + gupxz * gupyy)* Ayz
+
+  Rxz =    gupxx * gupxz * Axx + gupxy * gupyz * Ayy + gupxz * gupzz * Azz + &
+          (gupxx * gupyz       + gupxy * gupxz)* Axy                       + &
+          (gupxx * gupzz       + gupxz * gupxz)* Axz                       + &
+          (gupxy * gupzz       + gupxz * gupyz)* Ayz
+
+  Ryz =    gupxy * gupxz * Axx + gupyy * gupyz * Ayy + gupyz * gupzz * Azz + &
+          (gupxy * gupyz       + gupyy * gupxz)* Axy                       + &
+          (gupxy * gupzz       + gupyz * gupxz)* Axz                       + &
+          (gupyy * gupzz       + gupyz * gupyz)* Ayz
+
 ! Right hand side for Gam^i without shift terms...
   call fderivs(ex,Lap,Lapx,Lapy,Lapz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev)
   call fderivs(ex,trK,Kx,Ky,Kz,X,Y,Z,SYM,SYM,SYM,symmetry,Lev)
-  do k=1,ex(3)
-  do j=1,ex(2)
-  do i=1,ex(1)
-    gupxx_loc = gupxx(i,j,k)
-    gupxy_loc = gupxy(i,j,k)
-    gupxz_loc = gupxz(i,j,k)
-    gupyy_loc = gupyy(i,j,k)
-    gupyz_loc = gupyz(i,j,k)
-    gupzz_loc = gupzz(i,j,k)
 
-    Rxx_loc = gupxx_loc * gupxx_loc * Axx(i,j,k) + gupxy_loc * gupxy_loc * Ayy(i,j,k) + gupxz_loc * gupxz_loc * Azz(i,j,k) + &
+   Gamx_rhs = - TWO * (   Lapx * Rxx +   Lapy * Rxy +   Lapz * Rxz ) + &
-         TWO * (gupxx_loc * gupxy_loc * Axy(i,j,k) + gupxx_loc * gupxz_loc * Axz(i,j,k) + gupxy_loc * gupxz_loc * Ayz(i,j,k))
+        TWO * alpn1 * (                                                &
-    Ryy_loc = gupxy_loc * gupxy_loc * Axx(i,j,k) + gupyy_loc * gupyy_loc * Ayy(i,j,k) + gupyz_loc * gupyz_loc * Azz(i,j,k) + &
+        -F3o2/chin1 * (   chix * Rxx +   chiy * Rxy +   chiz * Rxz ) - &
-         TWO * (gupxy_loc * gupyy_loc * Axy(i,j,k) + gupxy_loc * gupyz_loc * Axz(i,j,k) + gupyy_loc * gupyz_loc * Ayz(i,j,k))
+              gupxx * (   F2o3 * Kx  +  EIGHT * PI * Sx            ) - &
-    Rzz_loc = gupxz_loc * gupxz_loc * Axx(i,j,k) + gupyz_loc * gupyz_loc * Ayy(i,j,k) + gupzz_loc * gupzz_loc * Azz(i,j,k) + &
+              gupxy * (   F2o3 * Ky  +  EIGHT * PI * Sy            ) - &
-         TWO * (gupxz_loc * gupyz_loc * Axy(i,j,k) + gupxz_loc * gupzz_loc * Axz(i,j,k) + gupyz_loc * gupzz_loc * Ayz(i,j,k))
+              gupxz * (   F2o3 * Kz  +  EIGHT * PI * Sz            ) + &
-    Rxy_loc = gupxx_loc * gupxy_loc * Axx(i,j,k) + gupxy_loc * gupyy_loc * Ayy(i,j,k) + gupxz_loc * gupyz_loc * Azz(i,j,k) + &
+                        Gamxxx * Rxx + Gamxyy * Ryy + Gamxzz * Rzz   + &
-         (gupxx_loc * gupyy_loc + gupxy_loc * gupxy_loc) * Axy(i,j,k) + &
+                TWO * ( Gamxxy * Rxy + Gamxxz * Rxz + Gamxyz * Ryz ) )
-         (gupxx_loc * gupyz_loc + gupxz_loc * gupxy_loc) * Axz(i,j,k) + &
-         (gupxy_loc * gupyz_loc + gupxz_loc * gupyy_loc) * Ayz(i,j,k)
-    Rxz_loc = gupxx_loc * gupxz_loc * Axx(i,j,k) + gupxy_loc * gupyz_loc * Ayy(i,j,k) + gupxz_loc * gupzz_loc * Azz(i,j,k) + &
-         (gupxx_loc * gupyz_loc + gupxy_loc * gupxz_loc) * Axy(i,j,k) + &
-         (gupxx_loc * gupzz_loc + gupxz_loc * gupxz_loc) * Axz(i,j,k) + &
-         (gupxy_loc * gupzz_loc + gupxz_loc * gupyz_loc) * Ayz(i,j,k)
-    Ryz_loc = gupxy_loc * gupxz_loc * Axx(i,j,k) + gupyy_loc * gupyz_loc * Ayy(i,j,k) + gupyz_loc * gupzz_loc * Azz(i,j,k) + &
-         (gupxy_loc * gupyz_loc + gupyy_loc * gupxz_loc) * Axy(i,j,k) + &
-         (gupxy_loc * gupzz_loc + gupyz_loc * gupxz_loc) * Axz(i,j,k) + &
-         (gupyy_loc * gupzz_loc + gupyz_loc * gupyz_loc) * Ayz(i,j,k)
-    Rxx(i,j,k) = Rxx_loc
-    Ryy(i,j,k) = Ryy_loc
-    Rzz(i,j,k) = Rzz_loc
-    Rxy(i,j,k) = Rxy_loc
-    Rxz(i,j,k) = Rxz_loc
-    Ryz(i,j,k) = Ryz_loc
 
-    Gamx_rhs(i,j,k) = - TWO * (Lapx(i,j,k) * Rxx_loc + Lapy(i,j,k) * Rxy_loc + Lapz(i,j,k) * Rxz_loc) + &
+   Gamy_rhs = - TWO * (   Lapx * Rxy +   Lapy * Ryy +   Lapz * Ryz ) + &
-         TWO * alpn1(i,j,k) * ( &
+        TWO * alpn1 * (                                                &
-         -F3o2/chin1(i,j,k) * (chix(i,j,k) * Rxx_loc + chiy(i,j,k) * Rxy_loc + chiz(i,j,k) * Rxz_loc) - &
+        -F3o2/chin1 * (   chix * Rxy +  chiy * Ryy +    chiz * Ryz ) - &
-         gupxx_loc * (F2o3 * Kx(i,j,k) + EIGHT * PI * Sx(i,j,k)) - &
+              gupxy * (   F2o3 * Kx  +  EIGHT * PI * Sx            ) - &
-         gupxy_loc * (F2o3 * Ky(i,j,k) + EIGHT * PI * Sy(i,j,k)) - &
+              gupyy * (   F2o3 * Ky  +  EIGHT * PI * Sy            ) - &
-         gupxz_loc * (F2o3 * Kz(i,j,k) + EIGHT * PI * Sz(i,j,k)) + &
+              gupyz * (   F2o3 * Kz  +  EIGHT * PI * Sz            ) + &
-         Gamxxx(i,j,k) * Rxx_loc + Gamxyy(i,j,k) * Ryy_loc + Gamxzz(i,j,k) * Rzz_loc + &
+                        Gamyxx * Rxx + Gamyyy * Ryy + Gamyzz * Rzz   + &
-         TWO * (Gamxxy(i,j,k) * Rxy_loc + Gamxxz(i,j,k) * Rxz_loc + Gamxyz(i,j,k) * Ryz_loc))
+                TWO * ( Gamyxy * Rxy + Gamyxz * Rxz + Gamyyz * Ryz ) )
 
-    Gamy_rhs(i,j,k) = - TWO * (Lapx(i,j,k) * Rxy_loc + Lapy(i,j,k) * Ryy_loc + Lapz(i,j,k) * Ryz_loc) + &
+   Gamz_rhs = - TWO * (   Lapx * Rxz +   Lapy * Ryz +   Lapz * Rzz ) + &
-         TWO * alpn1(i,j,k) * ( &
+        TWO * alpn1 * (                                                &
-         -F3o2/chin1(i,j,k) * (chix(i,j,k) * Rxy_loc + chiy(i,j,k) * Ryy_loc + chiz(i,j,k) * Ryz_loc) - &
+        -F3o2/chin1 * (   chix * Rxz +  chiy * Ryz +    chiz * Rzz ) - &
-         gupxy_loc * (F2o3 * Kx(i,j,k) + EIGHT * PI * Sx(i,j,k)) - &
+              gupxz * (   F2o3 * Kx  +  EIGHT * PI * Sx            ) - &
-         gupyy_loc * (F2o3 * Ky(i,j,k) + EIGHT * PI * Sy(i,j,k)) - &
+              gupyz * (   F2o3 * Ky  +  EIGHT * PI * Sy            ) - &
-         gupyz_loc * (F2o3 * Kz(i,j,k) + EIGHT * PI * Sz(i,j,k)) + &
+              gupzz * (   F2o3 * Kz  +  EIGHT * PI * Sz            ) + &
-         Gamyxx(i,j,k) * Rxx_loc + Gamyyy(i,j,k) * Ryy_loc + Gamyzz(i,j,k) * Rzz_loc + &
+                        Gamzxx * Rxx + Gamzyy * Ryy + Gamzzz * Rzz   + &
-         TWO * (Gamyxy(i,j,k) * Rxy_loc + Gamyxz(i,j,k) * Rxz_loc + Gamyyz(i,j,k) * Ryz_loc))
+                TWO * ( Gamzxy * Rxy + Gamzxz * Rxz + Gamzyz * Ryz ) )
-
-    Gamz_rhs(i,j,k) = - TWO * (Lapx(i,j,k) * Rxz_loc + Lapy(i,j,k) * Ryz_loc + Lapz(i,j,k) * Rzz_loc) + &
-         TWO * alpn1(i,j,k) * ( &
-         -F3o2/chin1(i,j,k) * (chix(i,j,k) * Rxz_loc + chiy(i,j,k) * Ryz_loc + chiz(i,j,k) * Rzz_loc) - &
-         gupxz_loc * (F2o3 * Kx(i,j,k) + EIGHT * PI * Sx(i,j,k)) - &
-         gupyz_loc * (F2o3 * Ky(i,j,k) + EIGHT * PI * Sy(i,j,k)) - &
-         gupzz_loc * (F2o3 * Kz(i,j,k) + EIGHT * PI * Sz(i,j,k)) + &
-         Gamzxx(i,j,k) * Rxx_loc + Gamzyy(i,j,k) * Ryy_loc + Gamzzz(i,j,k) * Rzz_loc + &
-         TWO * (Gamzxy(i,j,k) * Rxy_loc + Gamzxz(i,j,k) * Rxz_loc + Gamzyz(i,j,k) * Ryz_loc))
-  enddo
-  enddo
-  enddo
 
   call fdderivs(ex,betax,gxxx,gxyx,gxzx,gyyx,gyzx,gzzx,&
                 X,Y,Z,ANTI,SYM, SYM ,Symmetry,Lev)
@@ -359,54 +321,38 @@
   call fdderivs(ex,betaz,gxxz,gxyz,gxzz,gyyz,gyzz,gzzz,&
                 X,Y,Z,SYM ,SYM, ANTI,Symmetry,Lev)
 
+  fxx = gxxx + gxyy + gxzz
+  fxy = gxyx + gyyy + gyzz
+  fxz = gxzx + gyzy + gzzz
+
+  Gamxa =       gupxx * Gamxxx + gupyy * Gamxyy + gupzz * Gamxzz + &
+          TWO*( gupxy * Gamxxy + gupxz * Gamxxz + gupyz * Gamxyz )
+  Gamya =       gupxx * Gamyxx + gupyy * Gamyyy + gupzz * Gamyzz + &
+          TWO*( gupxy * Gamyxy + gupxz * Gamyxz + gupyz * Gamyyz )
+  Gamza =       gupxx * Gamzxx + gupyy * Gamzyy + gupzz * Gamzzz + &
+          TWO*( gupxy * Gamzxy + gupxz * Gamzxz + gupyz * Gamzyz )
+
   call fderivs(ex,Gamx,Gamxx,Gamxy,Gamxz,X,Y,Z,ANTI,SYM ,SYM ,Symmetry,Lev)
   call fderivs(ex,Gamy,Gamyx,Gamyy,Gamyz,X,Y,Z,SYM ,ANTI,SYM ,Symmetry,Lev)
   call fderivs(ex,Gamz,Gamzx,Gamzy,Gamzz,X,Y,Z,SYM ,SYM ,ANTI,Symmetry,Lev)
-  do k=1,ex(3)
-  do j=1,ex(2)
-  do i=1,ex(1)
-    divb_loc = div_beta(i,j,k)
-    fxx_loc = gxxx(i,j,k) + gxyy(i,j,k) + gxzz(i,j,k)
-    fxy_loc = gxyx(i,j,k) + gyyy(i,j,k) + gyzz(i,j,k)
-    fxz_loc = gxzx(i,j,k) + gyzy(i,j,k) + gzzz(i,j,k)
 
-    gupxx_loc = gupxx(i,j,k)
+  Gamx_rhs =               Gamx_rhs +  F2o3 *  Gamxa * div_beta        - &
-    gupxy_loc = gupxy(i,j,k)
+                     Gamxa * betaxx - Gamya * betaxy - Gamza * betaxz  + &
-    gupxz_loc = gupxz(i,j,k)
+             F1o3 * (gupxx * fxx    + gupxy * fxy    + gupxz * fxz    ) + &
-    gupyy_loc = gupyy(i,j,k)
+                     gupxx * gxxx   + gupyy * gyyx   + gupzz * gzzx    + &
-    gupyz_loc = gupyz(i,j,k)
+              TWO * (gupxy * gxyx   + gupxz * gxzx   + gupyz * gyzx  )
-    gupzz_loc = gupzz(i,j,k)
 
-    Gamxa_loc = gupxx_loc * Gamxxx(i,j,k) + gupyy_loc * Gamxyy(i,j,k) + gupzz_loc * Gamxzz(i,j,k) + &
+  Gamy_rhs =               Gamy_rhs +  F2o3 *  Gamya * div_beta        - &
-         TWO * (gupxy_loc * Gamxxy(i,j,k) + gupxz_loc * Gamxxz(i,j,k) + gupyz_loc * Gamxyz(i,j,k))
+                     Gamxa * betayx - Gamya * betayy - Gamza * betayz  + &
-    Gamya_loc = gupxx_loc * Gamyxx(i,j,k) + gupyy_loc * Gamyyy(i,j,k) + gupzz_loc * Gamyzz(i,j,k) + &
+             F1o3 * (gupxy * fxx    + gupyy * fxy    + gupyz * fxz    ) + &
-         TWO * (gupxy_loc * Gamyxy(i,j,k) + gupxz_loc * Gamyxz(i,j,k) + gupyz_loc * Gamyyz(i,j,k))
+                     gupxx * gxxy   + gupyy * gyyy   + gupzz * gzzy    + &
-    Gamza_loc = gupxx_loc * Gamzxx(i,j,k) + gupyy_loc * Gamzyy(i,j,k) + gupzz_loc * Gamzzz(i,j,k) + &
+              TWO * (gupxy * gxyy   + gupxz * gxzy   + gupyz * gyzy  )
-         TWO * (gupxy_loc * Gamzxy(i,j,k) + gupxz_loc * Gamzxz(i,j,k) + gupyz_loc * Gamzyz(i,j,k))
-    Gamxa(i,j,k) = Gamxa_loc
-    Gamya(i,j,k) = Gamya_loc
-    Gamza(i,j,k) = Gamza_loc
 
-    Gamx_rhs(i,j,k) = Gamx_rhs(i,j,k) + F2o3 * Gamxa_loc * divb_loc - &
+  Gamz_rhs =               Gamz_rhs +  F2o3 *  Gamza * div_beta        - &
-         Gamxa_loc * betaxx(i,j,k) - Gamya_loc * betaxy(i,j,k) - Gamza_loc * betaxz(i,j,k) + &
+                     Gamxa * betazx - Gamya * betazy - Gamza * betazz  + &
-         F1o3 * (gupxx_loc * fxx_loc + gupxy_loc * fxy_loc + gupxz_loc * fxz_loc) + &
+             F1o3 * (gupxz * fxx    + gupyz * fxy    + gupzz * fxz    ) + &
-         gupxx_loc * gxxx(i,j,k) + gupyy_loc * gyyx(i,j,k) + gupzz_loc * gzzx(i,j,k) + &
+                     gupxx * gxxz   + gupyy * gyyz   + gupzz * gzzz    + &
-         TWO * (gupxy_loc * gxyx(i,j,k) + gupxz_loc * gxzx(i,j,k) + gupyz_loc * gyzx(i,j,k))
+              TWO * (gupxy * gxyz   + gupxz * gxzz   + gupyz * gyzz  )    !rhs for Gam^i
-
-    Gamy_rhs(i,j,k) = Gamy_rhs(i,j,k) + F2o3 * Gamya_loc * divb_loc - &
-         Gamxa_loc * betayx(i,j,k) - Gamya_loc * betayy(i,j,k) - Gamza_loc * betayz(i,j,k) + &
-         F1o3 * (gupxy_loc * fxx_loc + gupyy_loc * fxy_loc + gupyz_loc * fxz_loc) + &
-         gupxx_loc * gxxy(i,j,k) + gupyy_loc * gyyy(i,j,k) + gupzz_loc * gzzy(i,j,k) + &
-         TWO * (gupxy_loc * gxyy(i,j,k) + gupxz_loc * gxzy(i,j,k) + gupyz_loc * gyzy(i,j,k))
-
-    Gamz_rhs(i,j,k) = Gamz_rhs(i,j,k) + F2o3 * Gamza_loc * divb_loc - &
-         Gamxa_loc * betazx(i,j,k) - Gamya_loc * betazy(i,j,k) - Gamza_loc * betazz(i,j,k) + &
-         F1o3 * (gupxz_loc * fxx_loc + gupyz_loc * fxy_loc + gupzz_loc * fxz_loc) + &
-         gupxx_loc * gxxz(i,j,k) + gupyy_loc * gyyz(i,j,k) + gupzz_loc * gzzz(i,j,k) + &
-         TWO * (gupxy_loc * gxyz(i,j,k) + gupxz_loc * gxzz(i,j,k) + gupyz_loc * gyzz(i,j,k))
-  enddo
-  enddo
-  enddo
 
 !first kind of connection stored in gij,k
   gxxx = gxx * Gamxxx + gxy * Gamyxx + gxz * Gamzxx
@@ -658,187 +604,189 @@
 !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)
+  fxx = fxx - Gamxxx * chix - Gamyxx * chiy - Gamzxx * chiz
-  do j=1,ex(2)
+  fxy = fxy - Gamxxy * chix - Gamyxy * chiy - Gamzxy * chiz
-  do i=1,ex(1)
+  fxz = fxz - Gamxxz * chix - Gamyxz * chiy - Gamzxz * chiz
-    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)
+  fyy = fyy - Gamxyy * chix - Gamyyy * chiy - Gamzyy * chiz
-    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)
+  fyz = fyz - Gamxyz * chix - Gamyyz * chiy - Gamzyz * chiz
-    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)
+  fzz = fzz - Gamxzz * chix - Gamyzz * chiy - Gamzzz * chiz
-    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)
+! Store D^l D_l chi - 3/(2*chi) D^l chi D_l chi in f
-    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 =        gupxx * ( fxx - F3o2/chin1 * chix * chix ) + &
-    f_loc = gupxx(i,j,k) * (fxx(i,j,k) - F3o2/chin_loc * chix(i,j,k) * chix(i,j,k)) + &
+             gupyy * ( fyy - F3o2/chin1 * chiy * chiy ) + &
-            gupyy(i,j,k) * (fyy(i,j,k) - F3o2/chin_loc * chiy(i,j,k) * chiy(i,j,k)) + &
+             gupzz * ( fzz - F3o2/chin1 * chiz * chiz ) + &
-            gupzz(i,j,k) * (fzz(i,j,k) - F3o2/chin_loc * chiz(i,j,k) * chiz(i,j,k)) + &
+       TWO * gupxy * ( fxy - F3o2/chin1 * chix * chiy ) + &
-            TWO * gupxy(i,j,k) * (fxy(i,j,k) - F3o2/chin_loc * chix(i,j,k) * chiy(i,j,k)) + &
+       TWO * gupxz * ( fxz - F3o2/chin1 * chix * chiz ) + &
-            TWO * gupxz(i,j,k) * (fxz(i,j,k) - F3o2/chin_loc * chix(i,j,k) * chiz(i,j,k)) + &
+       TWO * gupyz * ( fyz - F3o2/chin1 * chiy * chiz ) 
-            TWO * gupyz(i,j,k) * (fyz(i,j,k) - F3o2/chin_loc * chiy(i,j,k) * chiz(i,j,k))
+! Add chi part to Ricci tensor:
-    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
+  Rxx = Rxx + (fxx - chix*chix/chin1/TWO + gxx * f)/chin1/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
+  Ryy = Ryy + (fyy - chiy*chiy/chin1/TWO + gyy * f)/chin1/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
+  Rzz = Rzz + (fzz - chiz*chiz/chin1/TWO + gzz * f)/chin1/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
+  Rxy = Rxy + (fxy - chix*chiy/chin1/TWO + gxy * f)/chin1/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
+  Rxz = Rxz + (fxz - chix*chiz/chin1/TWO + gxz * f)/chin1/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
+  Ryz = Ryz + (fyz - chiy*chiz/chin1/TWO + gyz * f)/chin1/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)
+  gxxx = (gupxx * chix + gupxy * chiy + gupxz * chiz)/chin1
-  do j=1,ex(2)
+  gxxy = (gupxy * chix + gupyy * chiy + gupyz * chiz)/chin1
-  do i=1,ex(1)
+  gxxz = (gupxz * chix + gupyz * chiy + gupzz * chiz)/chin1
-    chin_loc = chin1(i,j,k)
+! now get physical second kind of connection
-    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
+  Gamxxx = Gamxxx - ( (chix + chix)/chin1 - gxx * gxxx )*HALF
-    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
+  Gamyxx = Gamyxx - (                     - gxx * gxxy )*HALF
-    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
+  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
 
-    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
+  fxx = fxx - Gamxxx*Lapx - Gamyxx*Lapy - Gamzxx*Lapz
-    Gamyxx(i,j,k) = Gamyxx(i,j,k) - (                                   - gxx(i,j,k) * gxxy(i,j,k) )*HALF
+  fyy = fyy - Gamxyy*Lapx - Gamyyy*Lapy - Gamzyy*Lapz
-    Gamzxx(i,j,k) = Gamzxx(i,j,k) - (                                   - gxx(i,j,k) * gxxz(i,j,k) )*HALF
+  fzz = fzz - Gamxzz*Lapx - Gamyzz*Lapy - Gamzzz*Lapz
-    Gamxyy(i,j,k) = Gamxyy(i,j,k) - (                                   - gyy(i,j,k) * gxxx(i,j,k) )*HALF
+  fxy = fxy - Gamxxy*Lapx - Gamyxy*Lapy - Gamzxy*Lapz
-    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
+  fxz = fxz - Gamxxz*Lapx - Gamyxz*Lapy - Gamzxz*Lapz
-    Gamzyy(i,j,k) = Gamzyy(i,j,k) - (                                   - gyy(i,j,k) * gxxz(i,j,k) )*HALF
+  fyz = fyz - Gamxyz*Lapx - Gamyyz*Lapy - Gamzyz*Lapz
-    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)
+! store D^i D_i Lap in trK_rhs upto chi
-    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)
+  trK_rhs =    gupxx * fxx + gupyy * fyy + gupzz * fzz + &
-    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)
+        TWO* ( gupxy * fxy + gupxz * fxz + gupyz * fyz )
-    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)
+#if 1        
-    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)
+!! follow bam code
-    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)
+  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)
   
-    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) + &
+  fxx = alpn1 * (Rxx - EIGHT * PI * Sxx) - fxx
-                     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))
+  fxy = alpn1 * (Rxy - EIGHT * PI * Sxy) - fxy
-  enddo
+  fxz = alpn1 * (Rxz - EIGHT * PI * Sxz) - fxz
-  enddo
+  fyy = alpn1 * (Ryy - EIGHT * PI * Syy) - fyy
-  enddo
+  fyz = alpn1 * (Ryz - EIGHT * PI * Syz) - fyz
-  do k=1,ex(3)
+  fzz = alpn1 * (Rzz - EIGHT * PI * Szz) - fzz
-  do j=1,ex(2)
+#else        
-  do i=1,ex(1)
+! Add lapse and S_ij parts to Ricci tensor:
-    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) + &
+  fxx = alpn1 * (Rxx - EIGHT * PI * Sxx) - fxx
-           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)) )
+  fxy = alpn1 * (Rxy - EIGHT * PI * Sxy) - fxy
-    S(i,j,k) = S_loc
+  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
 
-    f_loc = F2o3 * trK(i,j,k) * trK(i,j,k) - ( &
+! Compute trace-free part (note: chi^-1 and chi cancel!):
-            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) + &
+  f = F1o3 *(  gupxx * fxx + gupyy * fyy + gupzz * fzz + &
-            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)) + &
+        TWO* ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) )
-            alpn1(i,j,k)/chin_loc * f_loc )
+#endif
-    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)
+  Axx_rhs = fxx - gxx * f
-    l_fxy = alpn1(i,j,k) * (Rxy(i,j,k) - EIGHT * PI * Sxy(i,j,k)) - fxy(i,j,k)
+  Ayy_rhs = fyy - gyy * f
-    l_fxz = alpn1(i,j,k) * (Rxz(i,j,k) - EIGHT * PI * Sxz(i,j,k)) - fxz(i,j,k)
+  Azz_rhs = fzz - gzz * f
-    l_fyy = alpn1(i,j,k) * (Ryy(i,j,k) - EIGHT * PI * Syy(i,j,k)) - fyy(i,j,k)
+  Axy_rhs = fxy - gxy * f
-    l_fyz = alpn1(i,j,k) * (Ryz(i,j,k) - EIGHT * PI * Syz(i,j,k)) - fyz(i,j,k)
+  Axz_rhs = fxz - gxz * f
-    l_fzz = alpn1(i,j,k) * (Rzz(i,j,k) - EIGHT * PI * Szz(i,j,k)) - fzz(i,j,k)
+  Ayz_rhs = fyz - gyz * f
 
-    Axx_rhs(i,j,k) = l_fxx - gxx(i,j,k) * f_loc
+! Now: store A_il A^l_j into fij:
-    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) + &
+  fxx =       gupxx * Axx * Axx + gupyy * Axy * Axy + gupzz * Axz * Axz + &
-                 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) + &
+       TWO * (gupxy * Axx * Axy + gupxz * Axx * Axz + gupyz * Axy * Axz)
-                 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 =       gupxx * Axy * Axy + gupyy * Ayy * Ayy + gupzz * Ayz * Ayz + &
-    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) + &
+       TWO * (gupxy * Axy * Ayy + gupxz * Axy * Ayz + gupyz * Ayy * Ayz)
-                 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) + &
+  fzz =       gupxx * Axz * Axz + gupyy * Ayz * Ayz + gupzz * Azz * Azz + &
-                 gupxz(i,j,k) * Axy(i,j,k) * Ayz(i,j,k) + gupyz(i,j,k) * Ayy(i,j,k) * Ayz(i,j,k))
+       TWO * (gupxy * Axz * Ayz + gupxz * Axz * Azz + gupyz * Ayz * Azz)
-    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) + &
+  fxy =       gupxx * Axx * Axy + gupyy * Axy * Ayy + gupzz * Axz * Ayz + &
-                 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) + &
+              gupxy *(Axx * Ayy + Axy * Axy)                            + &
-                 gupxz(i,j,k) * Axz(i,j,k) * Azz(i,j,k) + gupyz(i,j,k) * Ayz(i,j,k) * Azz(i,j,k))
+              gupxz *(Axx * Ayz + Axz * Axy)                            + &
-    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) + &
+              gupyz *(Axy * Ayz + Axz * Ayy)
-                 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)) + &
+  fxz =       gupxx * Axx * Axz + gupyy * Axy * Ayz + gupzz * Axz * Azz + &
-                 gupxz(i,j,k) * (Axx(i,j,k) * Ayz(i,j,k) + Axz(i,j,k) * Axy(i,j,k)) + &
+              gupxy *(Axx * Ayz + Axy * Axz)                            + &
-                 gupyz(i,j,k) * (Axy(i,j,k) * Ayz(i,j,k) + Axz(i,j,k) * Ayy(i,j,k))
+              gupxz *(Axx * Azz + Axz * Axz)                            + &
-    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) + &
+              gupyz *(Axy * Azz + Axz * Ayz)
-                 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)) + &
+  fyz =       gupxx * Axy * Axz + gupyy * Ayy * Ayz + gupzz * Ayz * Azz + &
-                 gupxz(i,j,k) * (Axx(i,j,k) * Azz(i,j,k) + Axz(i,j,k) * Axz(i,j,k)) + &
+              gupxy *(Axy * Ayz + Ayy * Axz)                            + &
-                 gupyz(i,j,k) * (Axy(i,j,k) * Azz(i,j,k) + Axz(i,j,k) * Ayz(i,j,k))
+              gupxz *(Axy * Azz + Ayz * Axz)                            + &
-    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) + &
+              gupyz *(Ayy * Azz + Ayz * Ayz)
-                 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)
+  f = chin1
+! store D^i D_i Lap in trK_rhs
+  trK_rhs = f*trK_rhs
           
-    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)) + &
+  Axx_rhs =           f * Axx_rhs+ alpn1 * (trK * Axx - TWO * fxx)  + &
-                     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)) - &
+           TWO * (  Axx * betaxx +   Axy * betayx +   Axz * betazx )- &
-                     F2o3 * Axx(i,j,k) * divb_loc
+             F2o3 * Axx * div_beta
-    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) + &
+  Ayy_rhs =           f * Ayy_rhs+ alpn1 * (trK * Ayy - TWO * fyy)  + &
-                    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 * (  Axy * betaxy +   Ayy * betayy +   Ayz * betazy )- &
-                    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)) + &
+             F2o3 * Ayy * div_beta
-                    FOUR * PI * (rho(i,j,k) + S_loc) )
+
-  enddo
+  Azz_rhs =           f * Azz_rhs+ alpn1 * (trK * Azz - TWO * fzz)  + &
-  enddo
+           TWO * (  Axz * betaxz +   Ayz * betayz +   Azz * betazz )- &
-  enddo
+             F2o3 * Azz * div_beta
+
+  Axy_rhs =           f * Axy_rhs+ alpn1 *( trK * Axy  - TWO * fxy )+ &
+                    Axx * betaxy                  +   Axz * betazy  + &
+                                     Ayy * betayx +   Ayz * betazx  + &
+             F1o3 * Axy * div_beta                -   Axy * betazz
+
+  Ayz_rhs =           f * Ayz_rhs+ alpn1 *( trK * Ayz  - TWO * fyz )+ &
+                    Axy * betaxz +   Ayy * betayz                   + &
+                    Axz * betaxy                  +   Azz * betazy  + &
+             F1o3 * Ayz * div_beta                -   Ayz * betaxx
+ 
+  Axz_rhs =           f * Axz_rhs+ alpn1 *( trK * Axz  - TWO * fxz )+ &
+                    Axx * betaxz +   Axy * betayz                   + &
+                                     Ayz * betayx +   Azz * betazx  + &
+             F1o3 * Axz * div_beta                -   Axz * betayy      !rhs for Aij
+
+! Compute trace of S_ij
+
+  S =  f * ( gupxx * Sxx + gupyy * Syy + gupzz * Szz + &
+     TWO * ( gupxy * Sxy + gupxz * Sxz + gupyz * Syz ) )
+
+  trK_rhs = - trK_rhs + alpn1 *( F1o3 * trK * trK         + &
+                gupxx * fxx + gupyy * fyy + gupzz * fzz   + &
+        TWO * ( gupxy * fxy + gupxz * fxz + gupyz * fyz ) + &
+       FOUR * PI * ( rho + S ))                                !rhs for trK
   
 !!!! gauge variable part
 
@@ -1000,15 +948,15 @@
 !!!!!!!!!advection term + Kreiss-Oliger dissipation (merged for cache efficiency)
 ! lopsided_kodis shares the symmetry_bd buffer between advection and
 ! dissipation, eliminating redundant full-grid copies. For metric variables
-! gxx/gyy/gzz (=dxx/dyy/dzz+1): stencil coefficients sum to zero,
+! gxx/gyy/gzz (=dxx/dyy/dzz+1): kodis stencil coefficients sum to zero,
 ! so the constant offset has no effect on dissipation.
 
-  call lopsided_kodis(ex,X,Y,Z,dxx,gxx_rhs,betax,betay,betaz,Symmetry,SSS,eps)
+  call lopsided_kodis(ex,X,Y,Z,gxx,gxx_rhs,betax,betay,betaz,Symmetry,SSS,eps)
   call lopsided_kodis(ex,X,Y,Z,gxy,gxy_rhs,betax,betay,betaz,Symmetry,AAS,eps)
   call lopsided_kodis(ex,X,Y,Z,gxz,gxz_rhs,betax,betay,betaz,Symmetry,ASA,eps)
-  call lopsided_kodis(ex,X,Y,Z,dyy,gyy_rhs,betax,betay,betaz,Symmetry,SSS,eps)
+  call lopsided_kodis(ex,X,Y,Z,gyy,gyy_rhs,betax,betay,betaz,Symmetry,SSS,eps)
   call lopsided_kodis(ex,X,Y,Z,gyz,gyz_rhs,betax,betay,betaz,Symmetry,SAA,eps)
-  call lopsided_kodis(ex,X,Y,Z,dzz,gzz_rhs,betax,betay,betaz,Symmetry,SSS,eps)
+  call lopsided_kodis(ex,X,Y,Z,gzz,gzz_rhs,betax,betay,betaz,Symmetry,SSS,eps)
 
   call lopsided_kodis(ex,X,Y,Z,Axx,Axx_rhs,betax,betay,betaz,Symmetry,SSS,eps)
   call lopsided_kodis(ex,X,Y,Z,Axy,Axy_rhs,betax,betay,betaz,Symmetry,AAS,eps)

AMSS_NCKU_source/bssn_rhs_c.C

@@ -716,7 +716,6 @@ int f_compute_rhs_bssn(int *ex, double &T,
 
         // 24ms //
         fdderivs(ex,Lap,fxx,fxy,fxz,fyy,fyz,fzz,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev);
-        fderivs(ex,chi,dtSfx_rhs,dtSfy_rhs,dtSfz_rhs,X,Y,Z,SYM,SYM,SYM,Symmetry,Lev);
 
         // 6ms //
         for (int i=0;i<all;i+=1) {
@@ -1014,12 +1013,12 @@ int f_compute_rhs_bssn(int *ex, double &T,
             betaz_rhs[i] = FF * dtSfz[i];
 
             reta[i] =
-                gupxx[i] * dtSfx_rhs[i] * dtSfx_rhs[i]
+                gupxx[i] * chix[i] * chix[i]
-                + gupyy[i] * dtSfy_rhs[i] * dtSfy_rhs[i]
+                + gupyy[i] * chiy[i] * chiy[i]
-                + gupzz[i] * dtSfz_rhs[i] * dtSfz_rhs[i]
+                + gupzz[i] * chiz[i] * chiz[i]
-                + TWO * ( gupxy[i] * dtSfx_rhs[i] * dtSfy_rhs[i]
+                + TWO * ( gupxy[i] * chix[i] * chiy[i]
-                        + gupxz[i] * dtSfx_rhs[i] * dtSfz_rhs[i]
+                        + gupxz[i] * chix[i] * chiz[i]
-                        + gupyz[i] * dtSfy_rhs[i] * dtSfz_rhs[i] );
+                        + gupyz[i] * chiy[i] * chiz[i] );
 
             #if (GAUGE == 2)
             reta[i] = 1.31 / 2.0 * sqrt( reta[i] / chin1[i] ) / pow( (ONE - sqrt(chin1[i])), 2.0 );
@@ -1032,12 +1031,12 @@ int f_compute_rhs_bssn(int *ex, double &T,
             dtSfz_rhs[i] = Gamz_rhs[i] - reta[i] * dtSfz[i];
             #elif (GAUGE == 4 || GAUGE == 5)
             reta[i] =
-                gupxx[i] * dtSfx_rhs[i] * dtSfx_rhs[i]
+                gupxx[i] * chix[i] * chix[i]
-                + gupyy[i] * dtSfy_rhs[i] * dtSfy_rhs[i]
+                + gupyy[i] * chiy[i] * chiy[i]
-                + gupzz[i] * dtSfz_rhs[i] * dtSfz_rhs[i]
+                + gupzz[i] * chiz[i] * chiz[i]
-                + TWO * ( gupxy[i] * dtSfx_rhs[i] * dtSfy_rhs[i]
+                + TWO * ( gupxy[i] * chix[i] * chiy[i]
-                        + gupxz[i] * dtSfx_rhs[i] * dtSfz_rhs[i]
+                        + gupxz[i] * chix[i] * chiz[i]
-                        + gupyz[i] * dtSfy_rhs[i] * dtSfz_rhs[i] );
+                        + gupyz[i] * chiy[i] * chiz[i] );
 
             #if (GAUGE == 4)
             reta[i] = 1.31 / 2.0 * sqrt( reta[i] / chin1[i] ) / pow( (ONE - sqrt(chin1[i])), 2.0 );

AMSS_NCKU_source/makefile

@@ -64,8 +64,8 @@ lopsided_c.o: lopsided_c.C
 lopsided_kodis_c.o: lopsided_kodis_c.C
 	${CXX} $(CXXAPPFLAGS) -c $< $(filein) -o $@
 
-#interp_lb_profile.o: interp_lb_profile.C interp_lb_profile.h
+interp_lb_profile.o: interp_lb_profile.C interp_lb_profile.h
-#	${CXX} $(CXXAPPFLAGS) -c $< $(filein) -o $@
+	${CXX} $(CXXAPPFLAGS) -c $< $(filein) -o $@
 
 ## TwoPunctureABE uses fixed optimal flags with its own PGO profile, independent of CXXAPPFLAGS
 TP_PROFDATA = /home/$(shell whoami)/AMSS-NCKU/pgo_profile/TwoPunctureABE.profdata