fix(bssn_rhs)

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

AMSS_NCKU_source/bssn_class.C

@@ -5745,11 +5745,17 @@ void bssn_class::SHStep()
 
 //================================================================================================
 
-// 0: do not use mixing two levels data for OutBD; 1: do use
-
-#define MIXOUTB 0
-void bssn_class::RestrictProlong(int lev, int YN, bool BB,
-                                 MyList<var> *SL, MyList<var> *OL, MyList<var> *corL)
+// 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
 // StateList      1   -----------
 //
@@ -5811,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
 
-      Parallel::Sync_cached(GH->PatL[lev - 1], SynchList_pre, Symmetry, sync_cache_rp_coarse[lev]);
+#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();
@@ -5862,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
 
-      Parallel::Sync_cached(GH->PatL[lev - 1], SL, Symmetry, sync_cache_rp_coarse[lev]);
+#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();
@@ -5948,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
 
-      Parallel::Sync_cached(GH->PatL[lev - 1], SynchList_pre, Symmetry, sync_cache_rp_coarse[lev]);
+#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)
@@ -5970,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
 
-      Parallel::Sync_cached(GH->PatL[lev - 1], SL, Symmetry, sync_cache_rp_coarse[lev]);
+#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)
@@ -6035,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
 
-      Parallel::Sync_cached(GH->PatL[lev - 1], SynchList_pre, Symmetry, sync_cache_rp_coarse[lev]);
+#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)
@@ -6059,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
 
-      Parallel::Sync_cached(GH->PatL[lev - 1], StateList, Symmetry, sync_cache_rp_coarse[lev]);
+#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)
@@ -6083,7 +6101,7 @@ void bssn_class::RestrictProlong(int lev, int YN, bool BB)
 
 //================================================================================================
 
-void bssn_class::ProlongRestrict(int lev, int YN, bool BB)
+void bssn_class::ProlongRestrict(int lev, int YN, bool BB)
 {
   if (lev > 0)
   {
@@ -6136,15 +6154,18 @@ 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
-      Parallel::Sync_cached(GH->PatL[lev - 1], StateList, Symmetry, sync_cache_rp_coarse[lev]);
+#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 MIXOUTB
+#undef RP_SYNC_COARSE_AFTER_RESTRICT
+
+//================================================================================================
 
 
 

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]
-                + gupyy[i] * dtSfy_rhs[i] * dtSfy_rhs[i]
-                + gupzz[i] * dtSfz_rhs[i] * dtSfz_rhs[i]
-                + TWO * ( gupxy[i] * dtSfx_rhs[i] * dtSfy_rhs[i]
-                        + gupxz[i] * dtSfx_rhs[i] * dtSfz_rhs[i]
-                        + gupyz[i] * dtSfy_rhs[i] * dtSfz_rhs[i] );
+                gupxx[i] * chix[i] * chix[i]
+                + gupyy[i] * chiy[i] * chiy[i]
+                + gupzz[i] * chiz[i] * chiz[i]
+                + TWO * ( gupxy[i] * chix[i] * chiy[i]
+                        + gupxz[i] * chix[i] * chiz[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]
-                + gupyy[i] * dtSfy_rhs[i] * dtSfy_rhs[i]
-                + gupzz[i] * dtSfz_rhs[i] * dtSfz_rhs[i]
-                + TWO * ( gupxy[i] * dtSfx_rhs[i] * dtSfy_rhs[i]
-                        + gupxz[i] * dtSfx_rhs[i] * dtSfz_rhs[i]
-                        + gupyz[i] * dtSfy_rhs[i] * dtSfz_rhs[i] );
+                gupxx[i] * chix[i] * chix[i]
+                + gupyy[i] * chiy[i] * chiy[i]
+                + gupzz[i] * chiz[i] * chiz[i]
+                + TWO * ( gupxy[i] * chix[i] * chiy[i]
+                        + gupxz[i] * chix[i] * chiz[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 );
@@ -1139,59 +1138,59 @@ int f_compute_rhs_bssn(int *ex, double &T,
             fderivs(ex,Ayy,gyyx,gyyy,gyyz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0);
             fderivs(ex,Ayz,gyzx,gyzy,gyzz,X,Y,Z,SYM ,ANTI,ANTI,Symmetry,0);
             fderivs(ex,Azz,gzzx,gzzy,gzzz,X,Y,Z,SYM ,SYM ,SYM ,Symmetry,0);
-        }
-        // 7ms //
-        for (int i=0;i<all;i+=1) {
-            gxxx[i] = gxxx[i] - (  Gamxxx[i] * Axx[i] + Gamyxx[i] * Axy[i] + Gamzxx[i] * Axz[i]
-                                  + Gamxxx[i] * Axx[i] + Gamyxx[i] * Axy[i] + Gamzxx[i] * Axz[i]) - chix[i]*Axx[i]/chin1[i];
-            gxyx[i] = gxyx[i] - (  Gamxxy[i] * Axx[i] + Gamyxy[i] * Axy[i] + Gamzxy[i] * Axz[i]
-                                  + Gamxxx[i] * Axy[i] + Gamyxx[i] * Ayy[i] + Gamzxx[i] * Ayz[i]) - chix[i]*Axy[i]/chin1[i];
-            gxzx[i] = gxzx[i] - (  Gamxxz[i] * Axx[i] + Gamyxz[i] * Axy[i] + Gamzxz[i] * Axz[i]
-                                  + Gamxxx[i] * Axz[i] + Gamyxx[i] * Ayz[i] + Gamzxx[i] * Azz[i]) - chix[i]*Axz[i]/chin1[i];
-            gyyx[i] = gyyx[i] - (  Gamxxy[i] * Axy[i] + Gamyxy[i] * Ayy[i] + Gamzxy[i] * Ayz[i]
-                                  + Gamxxy[i] * Axy[i] + Gamyxy[i] * Ayy[i] + Gamzxy[i] * Ayz[i]) - chix[i]*Ayy[i]/chin1[i];
-            gyzx[i] = gyzx[i] - (  Gamxxz[i] * Axy[i] + Gamyxz[i] * Ayy[i] + Gamzxz[i] * Ayz[i]
-                                  + Gamxxy[i] * Axz[i] + Gamyxy[i] * Ayz[i] + Gamzxy[i] * Azz[i]) - chix[i]*Ayz[i]/chin1[i];
-            gzzx[i] = gzzx[i] - (  Gamxxz[i] * Axz[i] + Gamyxz[i] * Ayz[i] + Gamzxz[i] * Azz[i]
-                                  + Gamxxz[i] * Axz[i] + Gamyxz[i] * Ayz[i] + Gamzxz[i] * Azz[i]) - chix[i]*Azz[i]/chin1[i];
-            gxxy[i] = gxxy[i] - (  Gamxxy[i] * Axx[i] + Gamyxy[i] * Axy[i] + Gamzxy[i] * Axz[i]
-                                  + Gamxxy[i] * Axx[i] + Gamyxy[i] * Axy[i] + Gamzxy[i] * Axz[i]) - chiy[i]*Axx[i]/chin1[i];
-            gxyy[i] = gxyy[i] - (  Gamxyy[i] * Axx[i] + Gamyyy[i] * Axy[i] + Gamzyy[i] * Axz[i]
-                                  + Gamxxy[i] * Axy[i] + Gamyxy[i] * Ayy[i] + Gamzxy[i] * Ayz[i]) - chiy[i]*Axy[i]/chin1[i];
-            gxzy[i] = gxzy[i] - (  Gamxyz[i] * Axx[i] + Gamyyz[i] * Axy[i] + Gamzyz[i] * Axz[i]
-                                  + Gamxxy[i] * Axz[i] + Gamyxy[i] * Ayz[i] + Gamzxy[i] * Azz[i]) - chiy[i]*Axz[i]/chin1[i];
-            gyyy[i] = gyyy[i] - (  Gamxyy[i] * Axy[i] + Gamyyy[i] * Ayy[i] + Gamzyy[i] * Ayz[i]
-                                  + Gamxyy[i] * Axy[i] + Gamyyy[i] * Ayy[i] + Gamzyy[i] * Ayz[i]) - chiy[i]*Ayy[i]/chin1[i];
-            gyzy[i] = gyzy[i] - (  Gamxyz[i] * Axy[i] + Gamyyz[i] * Ayy[i] + Gamzyz[i] * Ayz[i]
-                                  + Gamxyy[i] * Axz[i] + Gamyyy[i] * Ayz[i] + Gamzyy[i] * Azz[i]) - chiy[i]*Ayz[i]/chin1[i];
-            gzzy[i] = gzzy[i] - (  Gamxyz[i] * Axz[i] + Gamyyz[i] * Ayz[i] + Gamzyz[i] * Azz[i]
-                                  + Gamxyz[i] * Axz[i] + Gamyyz[i] * Ayz[i] + Gamzyz[i] * Azz[i]) - chiy[i]*Azz[i]/chin1[i];
-            gxxz[i] = gxxz[i] - (  Gamxxz[i] * Axx[i] + Gamyxz[i] * Axy[i] + Gamzxz[i] * Axz[i]
-                                  + Gamxxz[i] * Axx[i] + Gamyxz[i] * Axy[i] + Gamzxz[i] * Axz[i]) - chiz[i]*Axx[i]/chin1[i];
-            gxyz[i] = gxyz[i] - (  Gamxyz[i] * Axx[i] + Gamyyz[i] * Axy[i] + Gamzyz[i] * Axz[i]
-                                  + Gamxxz[i] * Axy[i] + Gamyxz[i] * Ayy[i] + Gamzxz[i] * Ayz[i]) - chiz[i]*Axy[i]/chin1[i];
-            gxzz[i] = gxzz[i] - (  Gamxzz[i] * Axx[i] + Gamyzz[i] * Axy[i] + Gamzzz[i] * Axz[i]
-                                  + Gamxxz[i] * Axz[i] + Gamyxz[i] * Ayz[i] + Gamzxz[i] * Azz[i]) - chiz[i]*Axz[i]/chin1[i];
-            gyyz[i] = gyyz[i] - (  Gamxyz[i] * Axy[i] + Gamyyz[i] * Ayy[i] + Gamzyz[i] * Ayz[i]
-                                  + Gamxyz[i] * Axy[i] + Gamyyz[i] * Ayy[i] + Gamzyz[i] * Ayz[i]) - chiz[i]*Ayy[i]/chin1[i];
-            gyzz[i] = gyzz[i] - (  Gamxzz[i] * Axy[i] + Gamyzz[i] * Ayy[i] + Gamzzz[i] * Ayz[i]
-                                  + Gamxyz[i] * Axz[i] + Gamyyz[i] * Ayz[i] + Gamzyz[i] * Azz[i]) - chiz[i]*Ayz[i]/chin1[i];
-            gzzz[i] = gzzz[i] - (  Gamxzz[i] * Axz[i] + Gamyzz[i] * Ayz[i] + Gamzzz[i] * Azz[i]
-                                  + Gamxzz[i] * Axz[i] + Gamyzz[i] * Ayz[i] + Gamzzz[i] * Azz[i]) - chiz[i]*Azz[i]/chin1[i];
+            // 7ms //
+            for (int i=0;i<all;i+=1) {
+                gxxx[i] = gxxx[i] - (  Gamxxx[i] * Axx[i] + Gamyxx[i] * Axy[i] + Gamzxx[i] * Axz[i]
+                                      + Gamxxx[i] * Axx[i] + Gamyxx[i] * Axy[i] + Gamzxx[i] * Axz[i]) - chix[i]*Axx[i]/chin1[i];
+                gxyx[i] = gxyx[i] - (  Gamxxy[i] * Axx[i] + Gamyxy[i] * Axy[i] + Gamzxy[i] * Axz[i]
+                                      + Gamxxx[i] * Axy[i] + Gamyxx[i] * Ayy[i] + Gamzxx[i] * Ayz[i]) - chix[i]*Axy[i]/chin1[i];
+                gxzx[i] = gxzx[i] - (  Gamxxz[i] * Axx[i] + Gamyxz[i] * Axy[i] + Gamzxz[i] * Axz[i]
+                                      + Gamxxx[i] * Axz[i] + Gamyxx[i] * Ayz[i] + Gamzxx[i] * Azz[i]) - chix[i]*Axz[i]/chin1[i];
+                gyyx[i] = gyyx[i] - (  Gamxxy[i] * Axy[i] + Gamyxy[i] * Ayy[i] + Gamzxy[i] * Ayz[i]
+                                      + Gamxxy[i] * Axy[i] + Gamyxy[i] * Ayy[i] + Gamzxy[i] * Ayz[i]) - chix[i]*Ayy[i]/chin1[i];
+                gyzx[i] = gyzx[i] - (  Gamxxz[i] * Axy[i] + Gamyxz[i] * Ayy[i] + Gamzxz[i] * Ayz[i]
+                                      + Gamxxy[i] * Axz[i] + Gamyxy[i] * Ayz[i] + Gamzxy[i] * Azz[i]) - chix[i]*Ayz[i]/chin1[i];
+                gzzx[i] = gzzx[i] - (  Gamxxz[i] * Axz[i] + Gamyxz[i] * Ayz[i] + Gamzxz[i] * Azz[i]
+                                      + Gamxxz[i] * Axz[i] + Gamyxz[i] * Ayz[i] + Gamzxz[i] * Azz[i]) - chix[i]*Azz[i]/chin1[i];
+                gxxy[i] = gxxy[i] - (  Gamxxy[i] * Axx[i] + Gamyxy[i] * Axy[i] + Gamzxy[i] * Axz[i]
+                                      + Gamxxy[i] * Axx[i] + Gamyxy[i] * Axy[i] + Gamzxy[i] * Axz[i]) - chiy[i]*Axx[i]/chin1[i];
+                gxyy[i] = gxyy[i] - (  Gamxyy[i] * Axx[i] + Gamyyy[i] * Axy[i] + Gamzyy[i] * Axz[i]
+                                      + Gamxxy[i] * Axy[i] + Gamyxy[i] * Ayy[i] + Gamzxy[i] * Ayz[i]) - chiy[i]*Axy[i]/chin1[i];
+                gxzy[i] = gxzy[i] - (  Gamxyz[i] * Axx[i] + Gamyyz[i] * Axy[i] + Gamzyz[i] * Axz[i]
+                                      + Gamxxy[i] * Axz[i] + Gamyxy[i] * Ayz[i] + Gamzxy[i] * Azz[i]) - chiy[i]*Axz[i]/chin1[i];
+                gyyy[i] = gyyy[i] - (  Gamxyy[i] * Axy[i] + Gamyyy[i] * Ayy[i] + Gamzyy[i] * Ayz[i]
+                                      + Gamxyy[i] * Axy[i] + Gamyyy[i] * Ayy[i] + Gamzyy[i] * Ayz[i]) - chiy[i]*Ayy[i]/chin1[i];
+                gyzy[i] = gyzy[i] - (  Gamxyz[i] * Axy[i] + Gamyyz[i] * Ayy[i] + Gamzyz[i] * Ayz[i]
+                                      + Gamxyy[i] * Axz[i] + Gamyyy[i] * Ayz[i] + Gamzyy[i] * Azz[i]) - chiy[i]*Ayz[i]/chin1[i];
+                gzzy[i] = gzzy[i] - (  Gamxyz[i] * Axz[i] + Gamyyz[i] * Ayz[i] + Gamzyz[i] * Azz[i]
+                                      + Gamxyz[i] * Axz[i] + Gamyyz[i] * Ayz[i] + Gamzyz[i] * Azz[i]) - chiy[i]*Azz[i]/chin1[i];
+                gxxz[i] = gxxz[i] - (  Gamxxz[i] * Axx[i] + Gamyxz[i] * Axy[i] + Gamzxz[i] * Axz[i]
+                                      + Gamxxz[i] * Axx[i] + Gamyxz[i] * Axy[i] + Gamzxz[i] * Axz[i]) - chiz[i]*Axx[i]/chin1[i];
+                gxyz[i] = gxyz[i] - (  Gamxyz[i] * Axx[i] + Gamyyz[i] * Axy[i] + Gamzyz[i] * Axz[i]
+                                      + Gamxxz[i] * Axy[i] + Gamyxz[i] * Ayy[i] + Gamzxz[i] * Ayz[i]) - chiz[i]*Axy[i]/chin1[i];
+                gxzz[i] = gxzz[i] - (  Gamxzz[i] * Axx[i] + Gamyzz[i] * Axy[i] + Gamzzz[i] * Axz[i]
+                                      + Gamxxz[i] * Axz[i] + Gamyxz[i] * Ayz[i] + Gamzxz[i] * Azz[i]) - chiz[i]*Axz[i]/chin1[i];
+                gyyz[i] = gyyz[i] - (  Gamxyz[i] * Axy[i] + Gamyyz[i] * Ayy[i] + Gamzyz[i] * Ayz[i]
+                                      + Gamxyz[i] * Axy[i] + Gamyyz[i] * Ayy[i] + Gamzyz[i] * Ayz[i]) - chiz[i]*Ayy[i]/chin1[i];
+                gyzz[i] = gyzz[i] - (  Gamxzz[i] * Axy[i] + Gamyzz[i] * Ayy[i] + Gamzzz[i] * Ayz[i]
+                                      + Gamxyz[i] * Axz[i] + Gamyyz[i] * Ayz[i] + Gamzyz[i] * Azz[i]) - chiz[i]*Ayz[i]/chin1[i];
+                gzzz[i] = gzzz[i] - (  Gamxzz[i] * Axz[i] + Gamyzz[i] * Ayz[i] + Gamzzz[i] * Azz[i]
+                                      + Gamxzz[i] * Axz[i] + Gamyzz[i] * Ayz[i] + Gamzzz[i] * Azz[i]) - chiz[i]*Azz[i]/chin1[i];
 
-            movx_Res[i] = gupxx[i]*gxxx[i] + gupyy[i]*gxyy[i] + gupzz[i]*gxzz[i]
-                        + gupxy[i]*gxyx[i] + gupxz[i]*gxzx[i] + gupyz[i]*gxzy[i]
-                        + gupxy[i]*gxxy[i] + gupxz[i]*gxxz[i] + gupyz[i]*gxyz[i];
-            movy_Res[i] = gupxx[i]*gxyx[i] + gupyy[i]*gyyy[i] + gupzz[i]*gyzz[i]
-                        + gupxy[i]*gyyx[i] + gupxz[i]*gyzx[i] + gupyz[i]*gyzy[i]
-                        + gupxy[i]*gxyy[i] + gupxz[i]*gxyz[i] + gupyz[i]*gyyz[i];
-            movz_Res[i] = gupxx[i]*gxzx[i] + gupyy[i]*gyzy[i] + gupzz[i]*gzzz[i]
-                        + gupxy[i]*gyzx[i] + gupxz[i]*gzzx[i] + gupyz[i]*gzzy[i]
-                        + gupxy[i]*gxzy[i] + gupxz[i]*gxzz[i] + gupyz[i]*gyzz[i];
+                movx_Res[i] = gupxx[i]*gxxx[i] + gupyy[i]*gxyy[i] + gupzz[i]*gxzz[i]
+                            + gupxy[i]*gxyx[i] + gupxz[i]*gxzx[i] + gupyz[i]*gxzy[i]
+                            + gupxy[i]*gxxy[i] + gupxz[i]*gxxz[i] + gupyz[i]*gxyz[i];
+                movy_Res[i] = gupxx[i]*gxyx[i] + gupyy[i]*gyyy[i] + gupzz[i]*gyzz[i]
+                            + gupxy[i]*gyyx[i] + gupxz[i]*gyzx[i] + gupyz[i]*gyzy[i]
+                            + gupxy[i]*gxyy[i] + gupxz[i]*gxyz[i] + gupyz[i]*gyyz[i];
+                movz_Res[i] = gupxx[i]*gxzx[i] + gupyy[i]*gyzy[i] + gupzz[i]*gzzz[i]
+                            + gupxy[i]*gyzx[i] + gupxz[i]*gzzx[i] + gupyz[i]*gzzy[i]
+                            + gupxy[i]*gxzy[i] + gupxz[i]*gxzz[i] + gupyz[i]*gyzz[i];
 
-            movx_Res[i] = movx_Res[i] - F2o3*Kx[i] - F8*PI*Sx[i];
-            movy_Res[i] = movy_Res[i] - F2o3*Ky[i] - F8*PI*Sy[i];
-            movz_Res[i] = movz_Res[i] - F2o3*Kz[i] - F8*PI*Sz[i];
+                movx_Res[i] = movx_Res[i] - F2o3*Kx[i] - F8*PI*Sx[i];
+                movy_Res[i] = movy_Res[i] - F2o3*Ky[i] - F8*PI*Sy[i];
+                movz_Res[i] = movz_Res[i] - F2o3*Kz[i] - F8*PI*Sz[i];
+            }
         }