Vectorize 2

5 files changed, 377 insertions, 342 deletions

diff --git a/ML_BSSN/src/ML_BSSN_RHS.cc b/ML_BSSN/src/ML_BSSN_RHS.cc
index a3bb92e..1152240 100644
--- a/ML_BSSN/src/ML_BSSN_RHS.cc
+++ b/ML_BSSN/src/ML_BSSN_RHS.cc
@@ -2518,55 +2518,47 @@ static void ML_BSSN_RHS_Body(const cGH* restrict const cctkGH, const int dir, co
     }
 #else
     {
-        double gu[3][3] = {{ gu11, gu12, gu13 },
-                           { gu12, gu22, gu23 },
-                           { gu13, gu23, gu33 }};
-        double G[3][3][3] = {{{ G111, G112, G113 },
-                              { G112, G122, G123 },
-                              { G113, G123, G133 }},
-                             {{ G211, G212, G213 },
-                              { G212, G222, G223 },
-                              { G213, G223, G233 }},
-                             {{ G311, G312, G313 },
-                              { G312, G322, G323 },
-                              { G313, G323, G333 }}};
+    //    double gu[3][3] = {{ gu11, gu12, gu13 },
+    //                       { gu12, gu22, gu23 },
+    //                       { gu13, gu23, gu33 }};
+    //    double G[3][3][3] = {{{ G111, G112, G113 },
+    //                          { G112, G122, G123 },
+    //                          { G113, G123, G133 }},
+    //                         {{ G211, G212, G213 },
+    //                          { G212, G222, G223 },
+    //                          { G213, G223, G233 }},
+    //                         {{ G311, G312, G313 },
+    //                          { G312, G322, G323 },
+    //                          { G313, G323, G333 }}};
 
-        double X1 = 0.0, X3 = 0.0;
+    //    double X1 = 0.0, X3 = 0.0;
 
-        double dux_alpha, duz_alpha;
+    //    double dux_alpha, duz_alpha;
 
-        for (int fn = 0; fn < 3; fn++)
-            for (int bz = 0; bz < 3; bz++) {
-                X1 += gu[fn][bz] * G[0][fn][bz];
-                X3 += gu[fn][bz] * G[2][fn][bz];
-            }
+    //    for (int fn = 0; fn < 3; fn++)
+    //        for (int bz = 0; bz < 3; bz++) {
+    //            X1 += gu[fn][bz] * G[0][fn][bz];
+    //            X3 += gu[fn][bz] * G[2][fn][bz];
+    //        }
 
-        dux_alpha = gu[0][0] * PDstandardNth1alpha + gu[0][2] * PDstandardNth3alpha;
-        duz_alpha = gu[2][0] * PDstandardNth1alpha + gu[2][2] * PDstandardNth3alpha;
+    //    dux_alpha = gu[0][0] * PDstandardNth1alpha + gu[0][2] * PDstandardNth3alpha;
+    //    duz_alpha = gu[2][0] * PDstandardNth1alpha + gu[2][2] * PDstandardNth3alpha;
 
-        beta1rhsL = alphaL * alphaL * X1 - alphaL * dux_alpha - alphaL * BetaDriver * beta1L;
-        beta2rhsL = 0.0;
-        beta3rhsL = alphaL * alphaL * X3 - alphaL * duz_alpha - alphaL * BetaDriver * beta3L;
+    //    beta1rhsL = alphaL * alphaL * X1 - alphaL * dux_alpha - alphaL * BetaDriver * beta1L;
+    //    beta2rhsL = 0.0;
+    //    beta3rhsL = alphaL * alphaL * X3 - alphaL * duz_alpha - alphaL * BetaDriver * beta3L;
     }
 #endif
     else
     {
-        double xx = x[index];
-        double zz = z[index];
-        double t  = cctkGH->cctk_time;
-        //double rr = sqrt(xx * xx + zz * zz);
-        //double phi = atan2(zz, xx);
+      beta1rhsL = theta * ToReal(ShiftGammaCoeff) *
+          (ToReal(ShiftBCoeff) * B1L + (kone - ToReal(ShiftBCoeff)) * (Xt1L - eta * ToReal(BetaDriver) * beta1L));
 
-        //double gauge_r = -20.0 * (1.0 + 0.0 * cos(2 * phi) * rr * rr) * rr / (1.0 + rr * rr * rr * rr);
+      beta2rhsL = theta * ToReal(ShiftGammaCoeff) *
+          (ToReal(ShiftBCoeff) * B2L + (kone - ToReal(ShiftBCoeff)) * (Xt2L - eta * ToReal(BetaDriver) * beta2L));
 
-      beta1rhsL = theta * ShiftGammaCoeff *
-          (ShiftBCoeff * B1L + (1.0 - ShiftBCoeff) * (Xt1L - eta * BetaDriver * beta1L));
-
-      beta2rhsL = theta * ShiftGammaCoeff *
-          (ShiftBCoeff * B2L + (1.0 - ShiftBCoeff) * (Xt2L - eta * BetaDriver * beta2L));
-
-      beta3rhsL = theta * ShiftGammaCoeff *
-          (ShiftBCoeff * B3L + (1.0 - ShiftBCoeff) * (Xt3L - eta * BetaDriver * beta3L));
+      beta3rhsL = theta * ToReal(ShiftGammaCoeff) *
+          (ToReal(ShiftBCoeff) * B3L + (kone - ToReal(ShiftBCoeff)) * (Xt3L - eta * ToReal(BetaDriver) * beta3L));
     }
     
     CCTK_REAL_VEC B1rhsL CCTK_ATTRIBUTE_UNUSED = 
diff --git a/ML_BSSN/src/ML_BSSN_constraints.cc b/ML_BSSN/src/ML_BSSN_constraints.cc
index 62f8953..50c2e48 100644
--- a/ML_BSSN/src/ML_BSSN_constraints.cc
+++ b/ML_BSSN/src/ML_BSSN_constraints.cc
@@ -240,8 +240,6 @@ static void ML_BSSN_constraints_Body(const cGH* restrict const cctkGH, const int
   {
     const ptrdiff_t index CCTK_ATTRIBUTE_UNUSED = di*i + dj*j + dk*k;
     // vec_iter_counter+=CCTK_REAL_VEC_SIZE;
-    if (fabs(y[index]) > 1e-8)
-        continue;
     
     /* Assign local copies of grid functions */
     
@@ -822,51 +820,71 @@ static void ML_BSSN_constraints_Body(const cGH* restrict const cctkGH, const int
       default:
         CCTK_BUILTIN_UNREACHABLE();
     }
-    
-    double xx = x[index];
-    int origin = fabs(xx) < 1e-8;
+    CCTK_REAL_VEC kzero = ToReal(0.0);
+    CCTK_REAL_VEC kone = ToReal(1.0);
+    CCTK_REAL_VEC ktwo = ToReal(2.0);
+    CCTK_REAL_VEC xx = vec_load(x[index]);
+    CCTK_REAL_VEC x2 = kmul(xx, xx);
+
+    CCTK_REAL_VEC xinv  = kdiv(kone, xx);
+    CCTK_REAL_VEC x2inv = kdiv(kone, x2);
 
-    PDstandardNth2gt11  = 0;
-    PDstandardNth22gt11 = origin ? PDstandardNth11gt22 : PDstandardNth1gt11 / xx - 2 * (gt11[index] - gt22[index]) / (xx * xx);
-    PDstandardNth12gt11 = 0;
-    PDstandardNth23gt11 = 0;
-    PDstandardNth2gt12  = origin ? PDstandardNth1gt11 - PDstandardNth1gt22 : (gt11[index] - gt22[index]) / xx;
-    PDstandardNth22gt12 = 0;
-    PDstandardNth12gt12 = origin ? (PDstandardNth11gt11 - PDstandardNth11gt22) / 2 : (PDstandardNth1gt11 - PDstandardNth1gt22) / xx - (gt11[index] - gt22[index]) / (xx * xx);
-    PDstandardNth23gt12 = origin ? PDstandardNth13gt11 - PDstandardNth13gt22 : (PDstandardNth3gt11  - PDstandardNth3gt22) / xx;
-    PDstandardNth2gt13  = 0;
-    PDstandardNth22gt13 = origin ? PDstandardNth11gt13 / 2 : PDstandardNth1gt13 / xx - gt13[index] / (xx * xx);
-    PDstandardNth12gt13 = 0;
-    PDstandardNth23gt13 = 0;
-    PDstandardNth2gt22  = 0;
-    PDstandardNth22gt22 = origin ? PDstandardNth11gt11 : PDstandardNth1gt22 / xx + 2 * (gt11[index] - gt22[index]) / (xx * xx);
-    PDstandardNth12gt22 = 0;
-    PDstandardNth23gt22 = 0;
-    PDstandardNth2gt23  = origin ? PDstandardNth1gt13 : gt13[index] / xx;
-    PDstandardNth22gt23 = 0;
-    PDstandardNth12gt23 = origin ? PDstandardNth11gt13 / 2 : PDstandardNth1gt13 / xx - gt13[index] / (xx * xx);
-    PDstandardNth23gt23 = origin ? PDstandardNth13gt13 : PDstandardNth3gt13 / xx;
-    PDstandardNth2gt33  = 0;
-    PDstandardNth22gt33 = origin ? PDstandardNth11gt33 : PDstandardNth1gt33 / xx;
-    PDstandardNth12gt33 = 0;
-    PDstandardNth23gt33 = 0;
+    CCTK_REAL_VEC absx = kfabs(xx);
+    CCTK_REAL_VEC eps  = ToReal(1e-8);
+    CCTK_BOOLEAN_VEC origin = kcmplt(absx, eps);
 
-    PDstandardNth2phi = 0;
-    PDstandardNth22phi = origin ? PDstandardNth11phi : PDstandardNth1phi / xx;
-    PDstandardNth12phi = 0;
-    PDstandardNth23phi = 0;
-    PDstandardNth2trK = 0;
-    PDstandardNth2Xt1 = 0;
-    PDstandardNth2Xt2 = origin ? PDstandardNth1Xt1 : Xt1[index] / xx;
-    PDstandardNth2Xt3 = 0;
+    PDstandardNth2gt11  = kzero;
+    PDstandardNth22gt11 = kifthen(origin, PDstandardNth11gt22,
+                                          ksub(kmul(PDstandardNth1gt11, xinv), kmul(ktwo, kmul(ksub(gt11L, gt22L), x2inv))));
+    PDstandardNth12gt11 = kzero;
+    PDstandardNth23gt11 = kzero;
+    PDstandardNth2gt12  = kifthen(origin, ksub(PDstandardNth1gt11, PDstandardNth1gt22),
+                                          kmul(ksub(gt11L, gt22L), xinv));
+    PDstandardNth22gt12 = kzero;
+    PDstandardNth12gt12 = kifthen(origin, kmul(ksub(PDstandardNth11gt11, PDstandardNth11gt22), khalf),
+                                          ksub(kmul(ksub(PDstandardNth1gt11, PDstandardNth1gt22), xinv), kmul(ksub(gt11L, gt22L), x2inv)));
+    PDstandardNth23gt12 = kifthen(origin, ksub(PDstandardNth13gt11, PDstandardNth13gt22),
+                                          kmul(ksub(PDstandardNth3gt11, PDstandardNth3gt22), xinv));
+    PDstandardNth2gt13  = kzero;
+    PDstandardNth22gt13 = kifthen(origin, kmul(PDstandardNth11gt13, khalf),
+                                          ksub(kmul(PDstandardNth1gt13, xinv), kmul(gt13L, x2inv)));
+    PDstandardNth12gt13 = kzero;
+    PDstandardNth23gt13 = kzero;
+    PDstandardNth2gt22  = kzero;
+    PDstandardNth22gt22 = kifthen(origin, PDstandardNth11gt11,
+                                          kadd(kmul(PDstandardNth1gt22, xinv), kmul(ktwo, kmul(ksub(gt11L, gt22L), x2inv))));
+    PDstandardNth12gt22 = kzero;
+    PDstandardNth23gt22 = kzero;
+    PDstandardNth2gt23  = kifthen(origin, PDstandardNth1gt13, kmul(gt13L, xinv));
+    PDstandardNth22gt23 = kzero;
+    PDstandardNth12gt23 = kifthen(origin, kmul(PDstandardNth11gt13, khalf),
+                                          ksub(kmul(PDstandardNth1gt13, xinv), kmul(gt13L, x2inv)));
+    PDstandardNth23gt23 = kifthen(origin, PDstandardNth13gt13, kmul(PDstandardNth3gt13, xinv));
+    PDstandardNth2gt33  = kzero;
+    PDstandardNth22gt33 = kifthen(origin, PDstandardNth11gt33, kmul(PDstandardNth1gt33, xinv));
+    PDstandardNth12gt33 = kzero;
+    PDstandardNth23gt33 = kzero;
 
-    PDstandardNth2At11 = 0; //origin ? -2 * PDstandardNth1At12 : -2 * At12[index] / xx;
-    PDstandardNth2At12 = origin ? PDstandardNth1At11 - PDstandardNth1At22 : (At11[index] - At22[index]) / xx;
-    PDstandardNth2At13 = 0; //origin ? -PDstandardNth1At23 : - At23[index] / xx;
-    PDstandardNth2At22 = 0; //origin ? 2 * PDstandardNth1At12 : 2 * At12[index] / xx;
-    PDstandardNth2At23 = origin ? PDstandardNth1At13 : At13[index] / xx;
-    PDstandardNth2At33 = 0;
+    PDstandardNth2phi  = kzero;
+    PDstandardNth22phi = kifthen(origin, PDstandardNth11phi, kmul(PDstandardNth1phi, xinv));
+    PDstandardNth12phi = kzero;
+    PDstandardNth23phi = kzero;
+    PDstandardNth2Xt1  = kzero;
+    PDstandardNth2Xt2  = kifthen(origin, PDstandardNth1Xt1, kmul(Xt1L, xinv));
+    PDstandardNth2Xt3  = kzero;
+    PDstandardNth2trK  = kzero;
 
+    //PDstandardNth2At11 = origin ? -2 * PDstandardNth1At12 : -2 * At12[index] / xx;
+    //PDstandardNth2At13 = origin ? -PDstandardNth1At23 : - At23[index] / xx;
+    //PDstandardNth2At22 = origin ? 2 * PDstandardNth1At12 : 2 * At12[index] / xx;
+    PDstandardNth2At11 = kzero;
+    PDstandardNth2At12 = kifthen(origin, ksub(PDstandardNth1At11, PDstandardNth1At22),
+                                         kmul(ksub(At11L, At22L), xinv));
+    PDstandardNth2At13 = kzero;
+    PDstandardNth2At22 = kzero;
+    PDstandardNth2At23 = kifthen(origin, PDstandardNth1At13, kmul(At13L, xinv));
+    PDstandardNth2At33 = kzero;
+    
     /* Calculate temporaries and grid functions */
     CCTK_REAL_VEC JacPDstandardNth11gt11 CCTK_ATTRIBUTE_UNUSED;
     CCTK_REAL_VEC JacPDstandardNth11gt12 CCTK_ATTRIBUTE_UNUSED;
diff --git a/ML_BSSN/src/ML_BSSN_lapse_evol.cc b/ML_BSSN/src/ML_BSSN_lapse_evol.cc
index e688b5a..0388510 100644
--- a/ML_BSSN/src/ML_BSSN_lapse_evol.cc
+++ b/ML_BSSN/src/ML_BSSN_lapse_evol.cc
@@ -211,22 +211,47 @@ static void ML_BSSN_lapse_evol_Body(const cGH* restrict const cctkGH, const int
             break;
         }
 
-        PDupwindNthAnti2alpha = 0.0;
-        PDupwindNthSymm2alpha = 0.0;
 
+        PDupwindNthAnti2alpha = ToReal(0.0);
+        PDupwindNthSymm2alpha = ToReal(0.0);
+
+        CCTK_REAL_VEC alpharhsL;
+
+#if 1
         if (lapse_gh) {
-            const double e4phi = conformalMethod == 1 ? 1.0 / (phiL * phiL) : exp(4.0 * phiL);
-            const double det_gamma = pow(e4phi, 3.0);
+            CCTK_REAL_VEC e4phi = kifthen(conformalMethod == 1, ToReal(1.0) / (phiL * phiL), kexp(ToReal(4.0) * phiL));
+            CCTK_REAL_VEC det_gamma = kpow(e4phi, 3.0);
 
-            alpharhs[index] = -(alphaL * alphaL) * trKL + gh_eta_l_bar * pow(alphaL, gh_q) * (alphaL * alphaL) * log(pow(det_gamma, 0.5 * gh_p) / alphaL)
-                + LapseAdvectionCoeff * (beta1L * PDupwindNthAnti1alpha + beta3L * PDupwindNthAnti3alpha +
-                                         fabs(beta1L) * PDupwindNthSymm1alpha + fabs(beta3L) * PDupwindNthSymm3alpha);
+            alpharhsL = -(alphaL * alphaL) * trKL + ToReal(gh_eta_l_bar) * kpow(alphaL, gh_q) * (alphaL * alphaL) * klog(kpow(det_gamma, 0.5 * gh_p) / alphaL)
+                + ToReal(LapseAdvectionCoeff) * (beta1L * PDupwindNthAnti1alpha + beta3L * PDupwindNthAnti3alpha +
+                                         kfabs(beta1L) * PDupwindNthSymm1alpha + kfabs(beta3L) * PDupwindNthSymm3alpha);
         } else {
-            alpharhs[index] = -harmonicF * pow(alphaL, harmonicN) * trKL + WFactor * WL
+            alpharhsL = -ToReal(harmonicF) * kpow(alphaL, harmonicN) * trKL + ToReal(WFactor) * WL
+                + ToReal(LapseAdvectionCoeff) * (beta1L * PDupwindNthAnti1alpha + beta3L * PDupwindNthAnti3alpha +
+                                         kfabs(beta1L) * PDupwindNthSymm1alpha + kfabs(beta3L) * PDupwindNthSymm3alpha);
+        }
+#else
+        {
+            const double t = cctkGH->cctk_time;
+            const double fact_harmonic = alphaL * alphaL;
+            const double fact_1plog = 2.0 * alphaL;
+            double fact;
+
+            if (t < 5.0)
+                fact = fact_harmonic;
+            else if (t > 5.5)
+                fact = fact_1plog;
+            else
+                fact = fact_harmonic * exp(-36.0 * pow((t - 5.0) / 0.5, 4)) + fact_1plog * (1.0 -  exp(-36.0 * pow((t - 5.0) / 0.5, 4)));
+
+            alpharhs[index] = -fact * trKL
                 + LapseAdvectionCoeff * (beta1L * PDupwindNthAnti1alpha + beta3L * PDupwindNthAnti3alpha +
                                          fabs(beta1L) * PDupwindNthSymm1alpha + fabs(beta3L) * PDupwindNthSymm3alpha);
         }
+#endif
 
+        vec_store_partial_prepare(i,vecimin,vecimax);
+        vec_store_nta_partial(alpharhs[index],alpharhsL);
     }
     CCTK_ENDLOOP3STR(ML_BSSN_lapse_evol);
 }
diff --git a/ML_BSSN/src/eval_timediff.cc b/ML_BSSN/src/eval_timediff.cc
index ec35923..ae4f827 100644
--- a/ML_BSSN/src/eval_timediff.cc
+++ b/ML_BSSN/src/eval_timediff.cc
@@ -247,7 +247,7 @@ static void eval_timediff_body(const cGH* restrict const cctkGH, const int dir,
     
     CCTK_REAL_VEC eTttL, eTtxL, eTtyL, eTtzL, eTxxL, eTxyL, eTxzL, eTyyL, eTyzL, eTzzL CCTK_ATTRIBUTE_UNUSED ;
 
-    double beta_term_val;
+    //double beta_term_val;
     
     
     CCTK_REAL_VEC dJ111L, dJ112L, dJ113L, dJ122L, dJ123L, dJ133L, dJ211L, dJ212L, dJ213L, dJ222L, dJ223L, dJ233L, dJ311L, dJ312L, dJ313L, dJ322L, dJ323L, dJ333L, J11L, J12L, J13L, J21L, J22L, J23L, J31L, J32L, J33L CCTK_ATTRIBUTE_UNUSED ;
@@ -2107,65 +2107,65 @@ static void eval_timediff_body(const cGH* restrict const cctkGH, const int dir,
     CCTK_REAL_VEC G333 CCTK_ATTRIBUTE_UNUSED = 
       kadd(Gt333,kdiv(kmul(ToReal(0.5),kmadd(g33,kmul(kmadd(ddetg1,gu13,kmadd(ddetg2,gu23,kmul(ddetg3,gu33))),ToReal(-0.333333333333333333333333333333)),kmul(ddetg3,ToReal(2)))),detg));
 
-    {
-        double Gt[3][3][3] = {
-            {
-                { Gt111, Gt112, Gt113 },
-                { Gt112, Gt122, Gt123 },
-                { Gt113, Gt123, Gt133 },
-            },
-            {
-                { Gt211, Gt212, Gt213 },
-                { Gt212, Gt222, Gt223 },
-                { Gt213, Gt223, Gt233 },
-            },
-            {
-                { Gt311, Gt312, Gt313 },
-                { Gt312, Gt322, Gt323 },
-                { Gt313, Gt323, Gt333 },
-            }
-        };
-        double gt[3][3] = {
-            { gt11L, gt12L, gt13L },
-            { gt12L, gt22L, gt23L },
-            { gt13L, gt23L, gt33L },
-        };
-        double gtu[3][3] = {
-            { gtu11, gtu12, gtu13 },
-            { gtu12, gtu22, gtu23 },
-            { gtu13, gtu23, gtu33 },
-        };
+    //{
+    //    double Gt[3][3][3] = {
+    //        {
+    //            { Gt111, Gt112, Gt113 },
+    //            { Gt112, Gt122, Gt123 },
+    //            { Gt113, Gt123, Gt133 },
+    //        },
+    //        {
+    //            { Gt211, Gt212, Gt213 },
+    //            { Gt212, Gt222, Gt223 },
+    //            { Gt213, Gt223, Gt233 },
+    //        },
+    //        {
+    //            { Gt311, Gt312, Gt313 },
+    //            { Gt312, Gt322, Gt323 },
+    //            { Gt313, Gt323, Gt333 },
+    //        }
+    //    };
+    //    double gt[3][3] = {
+    //        { gt11L, gt12L, gt13L },
+    //        { gt12L, gt22L, gt23L },
+    //        { gt13L, gt23L, gt33L },
+    //    };
+    //    double gtu[3][3] = {
+    //        { gtu11, gtu12, gtu13 },
+    //        { gtu12, gtu22, gtu23 },
+    //        { gtu13, gtu23, gtu33 },
+    //    };
 
-        double dphi[3] = { PDstandardNth1phi, 0.0, PDstandardNth3phi };
+    //    double dphi[3] = { PDstandardNth1phi, 0.0, PDstandardNth3phi };
 
-        double G[3][3][3];
+    //    double G[3][3][3];
 
-        for (int idx0 = 0; idx0 < 3; idx0++)
-            for (int idx1 = 0; idx1 < 3; idx1++)
-                for (int idx2 = 0; idx2 < 3; idx2++) {
-                    G[idx0][idx1][idx2] = Gt[idx0][idx1][idx2] - (dphi[idx1] * (idx0 == idx2) + dphi[idx2] * (idx0 == idx1) - gt[idx1][idx2] * (gtu[idx0][0] * dphi[0] + gtu[idx0][2] * dphi[2])) / phiL;
-                }
-        G111 = G[0][0][0];
-        G122 = G[0][1][1];
-        G133 = G[0][2][2];
-        G112 = G[0][0][1];
-        G113 = G[0][0][2];
-        G123 = G[0][1][2];
+    //    for (int idx0 = 0; idx0 < 3; idx0++)
+    //        for (int idx1 = 0; idx1 < 3; idx1++)
+    //            for (int idx2 = 0; idx2 < 3; idx2++) {
+    //                G[idx0][idx1][idx2] = Gt[idx0][idx1][idx2] - (dphi[idx1] * (idx0 == idx2) + dphi[idx2] * (idx0 == idx1) - gt[idx1][idx2] * (gtu[idx0][0] * dphi[0] + gtu[idx0][2] * dphi[2])) / phiL;
+    //            }
+    //    G111 = G[0][0][0];
+    //    G122 = G[0][1][1];
+    //    G133 = G[0][2][2];
+    //    G112 = G[0][0][1];
+    //    G113 = G[0][0][2];
+    //    G123 = G[0][1][2];
 
-        G211 = G[1][0][0];
-        G222 = G[1][1][1];
-        G233 = G[1][2][2];
-        G212 = G[1][0][1];
-        G213 = G[1][0][2];
-        G223 = G[1][1][2];
+    //    G211 = G[1][0][0];
+    //    G222 = G[1][1][1];
+    //    G233 = G[1][2][2];
+    //    G212 = G[1][0][1];
+    //    G213 = G[1][0][2];
+    //    G223 = G[1][1][2];
 
-        G311 = G[2][0][0];
-        G322 = G[2][1][1];
-        G333 = G[2][2][2];
-        G312 = G[2][0][1];
-        G313 = G[2][0][2];
-        G323 = G[2][1][2];
-    }
+    //    G311 = G[2][0][0];
+    //    G322 = G[2][1][1];
+    //    G333 = G[2][2][2];
+    //    G312 = G[2][0][1];
+    //    G313 = G[2][0][2];
+    //    G323 = G[2][1][2];
+    //}
     
     CCTK_REAL_VEC K11 CCTK_ATTRIBUTE_UNUSED = 
       kmadd(At11L,e4phi,kmul(trKL,kmul(g11,ToReal(0.333333333333333333333333333333))));
@@ -2366,237 +2366,237 @@ static void eval_timediff_body(const cGH* restrict const cctkGH, const int dir,
     CCTK_REAL_VEC Xtdot3L CCTK_ATTRIBUTE_UNUSED = dotXt3;
     
     CCTK_REAL_VEC Kdot11L CCTK_ATTRIBUTE_UNUSED = 
-      G111 * JacPDstandardNth1alpha + G211 * JacPDstandardNth2alpha + G311 * JacPDstandardNth3alpha + alphaL * (trKL * K11 + R11 - 2.0 * (K11 * Km11 + K12 * Km12  + K13 * Km13)) - JacPDstandardNth11alpha;
+      G111 * JacPDstandardNth1alpha + G211 * JacPDstandardNth2alpha + G311 * JacPDstandardNth3alpha + alphaL * (trKL * K11 + R11 - ktwo * (K11 * Km11 + K12 * Km12  + K13 * Km13)) - JacPDstandardNth11alpha;
     
     CCTK_REAL_VEC Kdot12L CCTK_ATTRIBUTE_UNUSED = 
-      G112 * JacPDstandardNth1alpha + G212 * JacPDstandardNth2alpha + G312 * JacPDstandardNth3alpha + alphaL * (trKL * K12 + R12 - 2.0 * (K11 * Km21 + K12 * Km22  + K13 * Km23)) - JacPDstandardNth12alpha;
+      G112 * JacPDstandardNth1alpha + G212 * JacPDstandardNth2alpha + G312 * JacPDstandardNth3alpha + alphaL * (trKL * K12 + R12 - ktwo * (K11 * Km21 + K12 * Km22  + K13 * Km23)) - JacPDstandardNth12alpha;
     
     CCTK_REAL_VEC Kdot13L CCTK_ATTRIBUTE_UNUSED = 
-      G113 * JacPDstandardNth1alpha + G213 * JacPDstandardNth2alpha + G313 * JacPDstandardNth3alpha + alphaL * (trKL * K13 + R13 - 2.0 * (K11 * Km31 + K12 * Km32  + K13 * Km33)) - JacPDstandardNth13alpha;
+      G113 * JacPDstandardNth1alpha + G213 * JacPDstandardNth2alpha + G313 * JacPDstandardNth3alpha + alphaL * (trKL * K13 + R13 - ktwo * (K11 * Km31 + K12 * Km32  + K13 * Km33)) - JacPDstandardNth13alpha;
     
     CCTK_REAL_VEC Kdot22L CCTK_ATTRIBUTE_UNUSED = 
-      G122 * JacPDstandardNth1alpha + G222 * JacPDstandardNth2alpha + G322 * JacPDstandardNth3alpha + alphaL * (trKL * K22 + R22 - 2.0 * (K12 * Km21 + K22 * Km22  + K23 * Km23)) - JacPDstandardNth22alpha;
+      G122 * JacPDstandardNth1alpha + G222 * JacPDstandardNth2alpha + G322 * JacPDstandardNth3alpha + alphaL * (trKL * K22 + R22 - ktwo * (K12 * Km21 + K22 * Km22  + K23 * Km23)) - JacPDstandardNth22alpha;
     
     CCTK_REAL_VEC Kdot23L CCTK_ATTRIBUTE_UNUSED = 
-      G123 * JacPDstandardNth1alpha + G223 * JacPDstandardNth2alpha + G323 * JacPDstandardNth3alpha + alphaL * (trKL * K23 + R23 - 2.0 * (K12 * Km31 + K22 * Km32  + K23 * Km33)) - JacPDstandardNth23alpha;
+      G123 * JacPDstandardNth1alpha + G223 * JacPDstandardNth2alpha + G323 * JacPDstandardNth3alpha + alphaL * (trKL * K23 + R23 - ktwo * (K12 * Km31 + K22 * Km32  + K23 * Km33)) - JacPDstandardNth23alpha;
     
     CCTK_REAL_VEC Kdot33L CCTK_ATTRIBUTE_UNUSED = 
-      G133 * JacPDstandardNth1alpha + G233 * JacPDstandardNth2alpha + G333 * JacPDstandardNth3alpha + alphaL * (trKL * K33 + R33 - 2.0 * (K13 * Km31 + K23 * Km32  + K33 * Km33)) - JacPDstandardNth33alpha;
+      G133 * JacPDstandardNth1alpha + G233 * JacPDstandardNth2alpha + G333 * JacPDstandardNth3alpha + alphaL * (trKL * K33 + R33 - ktwo * (K13 * Km31 + K23 * Km32  + K33 * Km33)) - JacPDstandardNth33alpha;
 
     /* compute the beta term */
-    {
-        const double gt[3][3] = {
-            { gt11L, gt12L, gt13L },
-            { gt12L, gt22L, gt23L },
-            { gt13L, gt23L, gt33L },
-        };
+    //{
+    //    const double gt[3][3] = {
+    //        { gt11L, gt12L, gt13L },
+    //        { gt12L, gt22L, gt23L },
+    //        { gt13L, gt23L, gt33L },
+    //    };
 
-        const double gtu[3][3] = {
-            { gtu11, gtu12, gtu13 },
-            { gtu12, gtu22, gtu23 },
-            { gtu13, gtu23, gtu33 },
-        };
+    //    const double gtu[3][3] = {
+    //        { gtu11, gtu12, gtu13 },
+    //        { gtu12, gtu22, gtu23 },
+    //        { gtu13, gtu23, gtu33 },
+    //    };
 
-        const double dgt[3][3][3] = {
-            {
-                { PDstandardNth1gt11, PDstandardNth2gt11, PDstandardNth3gt11 },
-                { PDstandardNth1gt12, PDstandardNth2gt12, PDstandardNth3gt12 },
-                { PDstandardNth1gt13, PDstandardNth2gt13, PDstandardNth3gt13 },
-            },
-            {
-                { PDstandardNth1gt12, PDstandardNth2gt12, PDstandardNth3gt12 },
-                { PDstandardNth1gt22, PDstandardNth2gt22, PDstandardNth3gt22 },
-                { PDstandardNth1gt23, PDstandardNth2gt23, PDstandardNth3gt23 },
-            },
-            {
-                { PDstandardNth1gt13, PDstandardNth2gt13, PDstandardNth3gt13 },
-                { PDstandardNth1gt23, PDstandardNth2gt23, PDstandardNth3gt23 },
-                { PDstandardNth1gt33, PDstandardNth2gt33, PDstandardNth3gt33 },
-            },
-        };
+    //    const double dgt[3][3][3] = {
+    //        {
+    //            { PDstandardNth1gt11, PDstandardNth2gt11, PDstandardNth3gt11 },
+    //            { PDstandardNth1gt12, PDstandardNth2gt12, PDstandardNth3gt12 },
+    //            { PDstandardNth1gt13, PDstandardNth2gt13, PDstandardNth3gt13 },
+    //        },
+    //        {
+    //            { PDstandardNth1gt12, PDstandardNth2gt12, PDstandardNth3gt12 },
+    //            { PDstandardNth1gt22, PDstandardNth2gt22, PDstandardNth3gt22 },
+    //            { PDstandardNth1gt23, PDstandardNth2gt23, PDstandardNth3gt23 },
+    //        },
+    //        {
+    //            { PDstandardNth1gt13, PDstandardNth2gt13, PDstandardNth3gt13 },
+    //            { PDstandardNth1gt23, PDstandardNth2gt23, PDstandardNth3gt23 },
+    //            { PDstandardNth1gt33, PDstandardNth2gt33, PDstandardNth3gt33 },
+    //        },
+    //    };
 
-        const double d2gt[3][3][3][3] = {
-            {
-                {
-                    { PDstandardNth11gt11, PDstandardNth12gt11, PDstandardNth13gt11 },
-                    { PDstandardNth12gt11, PDstandardNth22gt11, PDstandardNth23gt11 },
-                    { PDstandardNth13gt11, PDstandardNth23gt11, PDstandardNth33gt11 },
-                },
-                {
-                    { PDstandardNth11gt12, PDstandardNth12gt12, PDstandardNth13gt12 },
-                    { PDstandardNth12gt12, PDstandardNth22gt12, PDstandardNth23gt12 },
-                    { PDstandardNth13gt12, PDstandardNth23gt12, PDstandardNth33gt12 },
-                },
-                {
-                    { PDstandardNth11gt13, PDstandardNth12gt13, PDstandardNth13gt13 },
-                    { PDstandardNth12gt13, PDstandardNth22gt13, PDstandardNth23gt13 },
-                    { PDstandardNth13gt13, PDstandardNth23gt13, PDstandardNth33gt13 },
-                },
-            },
-            {
-                {
-                    { PDstandardNth11gt12, PDstandardNth12gt12, PDstandardNth13gt12 },
-                    { PDstandardNth12gt12, PDstandardNth22gt12, PDstandardNth23gt12 },
-                    { PDstandardNth13gt12, PDstandardNth23gt12, PDstandardNth33gt12 },
-                },
-                {
-                    { PDstandardNth11gt22, PDstandardNth12gt22, PDstandardNth13gt22 },
-                    { PDstandardNth12gt22, PDstandardNth22gt22, PDstandardNth23gt22 },
-                    { PDstandardNth13gt22, PDstandardNth23gt22, PDstandardNth33gt22 },
-                },
-                {
-                    { PDstandardNth11gt23, PDstandardNth12gt23, PDstandardNth13gt23 },
-                    { PDstandardNth12gt23, PDstandardNth22gt23, PDstandardNth23gt23 },
-                    { PDstandardNth13gt23, PDstandardNth23gt23, PDstandardNth33gt23 },
-                },
-            },
-            {
-                {
-                    { PDstandardNth11gt13, PDstandardNth12gt13, PDstandardNth13gt13 },
-                    { PDstandardNth12gt13, PDstandardNth22gt13, PDstandardNth23gt13 },
-                    { PDstandardNth13gt13, PDstandardNth23gt13, PDstandardNth33gt13 },
-                },
-                {
-                    { PDstandardNth11gt23, PDstandardNth12gt23, PDstandardNth13gt23 },
-                    { PDstandardNth12gt23, PDstandardNth22gt23, PDstandardNth23gt23 },
-                    { PDstandardNth13gt23, PDstandardNth23gt23, PDstandardNth33gt23 },
-                },
-                {
-                    { PDstandardNth11gt33, PDstandardNth12gt33, PDstandardNth13gt33 },
-                    { PDstandardNth12gt33, PDstandardNth22gt33, PDstandardNth23gt33 },
-                    { PDstandardNth13gt33, PDstandardNth23gt33, PDstandardNth33gt33 },
-                },
-            },
-        };
+    //    const double d2gt[3][3][3][3] = {
+    //        {
+    //            {
+    //                { PDstandardNth11gt11, PDstandardNth12gt11, PDstandardNth13gt11 },
+    //                { PDstandardNth12gt11, PDstandardNth22gt11, PDstandardNth23gt11 },
+    //                { PDstandardNth13gt11, PDstandardNth23gt11, PDstandardNth33gt11 },
+    //            },
+    //            {
+    //                { PDstandardNth11gt12, PDstandardNth12gt12, PDstandardNth13gt12 },
+    //                { PDstandardNth12gt12, PDstandardNth22gt12, PDstandardNth23gt12 },
+    //                { PDstandardNth13gt12, PDstandardNth23gt12, PDstandardNth33gt12 },
+    //            },
+    //            {
+    //                { PDstandardNth11gt13, PDstandardNth12gt13, PDstandardNth13gt13 },
+    //                { PDstandardNth12gt13, PDstandardNth22gt13, PDstandardNth23gt13 },
+    //                { PDstandardNth13gt13, PDstandardNth23gt13, PDstandardNth33gt13 },
+    //            },
+    //        },
+    //        {
+    //            {
+    //                { PDstandardNth11gt12, PDstandardNth12gt12, PDstandardNth13gt12 },
+    //                { PDstandardNth12gt12, PDstandardNth22gt12, PDstandardNth23gt12 },
+    //                { PDstandardNth13gt12, PDstandardNth23gt12, PDstandardNth33gt12 },
+    //            },
+    //            {
+    //                { PDstandardNth11gt22, PDstandardNth12gt22, PDstandardNth13gt22 },
+    //                { PDstandardNth12gt22, PDstandardNth22gt22, PDstandardNth23gt22 },
+    //                { PDstandardNth13gt22, PDstandardNth23gt22, PDstandardNth33gt22 },
+    //            },
+    //            {
+    //                { PDstandardNth11gt23, PDstandardNth12gt23, PDstandardNth13gt23 },
+    //                { PDstandardNth12gt23, PDstandardNth22gt23, PDstandardNth23gt23 },
+    //                { PDstandardNth13gt23, PDstandardNth23gt23, PDstandardNth33gt23 },
+    //            },
+    //        },
+    //        {
+    //            {
+    //                { PDstandardNth11gt13, PDstandardNth12gt13, PDstandardNth13gt13 },
+    //                { PDstandardNth12gt13, PDstandardNth22gt13, PDstandardNth23gt13 },
+    //                { PDstandardNth13gt13, PDstandardNth23gt13, PDstandardNth33gt13 },
+    //            },
+    //            {
+    //                { PDstandardNth11gt23, PDstandardNth12gt23, PDstandardNth13gt23 },
+    //                { PDstandardNth12gt23, PDstandardNth22gt23, PDstandardNth23gt23 },
+    //                { PDstandardNth13gt23, PDstandardNth23gt23, PDstandardNth33gt23 },
+    //            },
+    //            {
+    //                { PDstandardNth11gt33, PDstandardNth12gt33, PDstandardNth13gt33 },
+    //                { PDstandardNth12gt33, PDstandardNth22gt33, PDstandardNth23gt33 },
+    //                { PDstandardNth13gt33, PDstandardNth23gt33, PDstandardNth33gt33 },
+    //            },
+    //        },
+    //    };
 
-        const double Ric[3][3] = {
-            { R11, R12, R13 },
-            { R12, R22, R23 },
-            { R13, R23, R33 },
-        };
+    //    const double Ric[3][3] = {
+    //        { R11, R12, R13 },
+    //        { R12, R22, R23 },
+    //        { R13, R23, R33 },
+    //    };
 
-        const double betaval[3] = {beta1L, 0.0, beta3L};
-        const double dbeta[3][3] = {{ PDstandardNth1beta1, PDstandardNth2beta1, PDstandardNth3beta1 },
-                                    { PDstandardNth1beta2, PDstandardNth2beta2, PDstandardNth3beta2 },
-                                    { PDstandardNth1beta3, PDstandardNth2beta3, PDstandardNth3beta3 }};
+    //    const double betaval[3] = {beta1L, 0.0, beta3L};
+    //    const double dbeta[3][3] = {{ PDstandardNth1beta1, PDstandardNth2beta1, PDstandardNth3beta1 },
+    //                                { PDstandardNth1beta2, PDstandardNth2beta2, PDstandardNth3beta2 },
+    //                                { PDstandardNth1beta3, PDstandardNth2beta3, PDstandardNth3beta3 }};
 
-        const double d2beta[3][3][3] = {
-            {
-                { PDstandardNth11beta1, PDstandardNth12beta1, PDstandardNth13beta1 },
-                { PDstandardNth12beta1, PDstandardNth22beta1, PDstandardNth23beta1 },
-                { PDstandardNth13beta1, PDstandardNth23beta1, PDstandardNth33beta1 },
-            },
-            {
-                { PDstandardNth11beta2, PDstandardNth12beta2, PDstandardNth13beta2 },
-                { PDstandardNth12beta2, PDstandardNth22beta2, PDstandardNth23beta2 },
-                { PDstandardNth13beta2, PDstandardNth23beta2, PDstandardNth33beta2 },
-            },
-            {
-                { PDstandardNth11beta3, PDstandardNth12beta3, PDstandardNth13beta3 },
-                { PDstandardNth12beta3, PDstandardNth22beta3, PDstandardNth23beta3 },
-                { PDstandardNth13beta3, PDstandardNth23beta3, PDstandardNth33beta3 },
-            },
-        };
+    //    const double d2beta[3][3][3] = {
+    //        {
+    //            { PDstandardNth11beta1, PDstandardNth12beta1, PDstandardNth13beta1 },
+    //            { PDstandardNth12beta1, PDstandardNth22beta1, PDstandardNth23beta1 },
+    //            { PDstandardNth13beta1, PDstandardNth23beta1, PDstandardNth33beta1 },
+    //        },
+    //        {
+    //            { PDstandardNth11beta2, PDstandardNth12beta2, PDstandardNth13beta2 },
+    //            { PDstandardNth12beta2, PDstandardNth22beta2, PDstandardNth23beta2 },
+    //            { PDstandardNth13beta2, PDstandardNth23beta2, PDstandardNth33beta2 },
+    //        },
+    //        {
+    //            { PDstandardNth11beta3, PDstandardNth12beta3, PDstandardNth13beta3 },
+    //            { PDstandardNth12beta3, PDstandardNth22beta3, PDstandardNth23beta3 },
+    //            { PDstandardNth13beta3, PDstandardNth23beta3, PDstandardNth33beta3 },
+    //        },
+    //    };
 
-        const double G[3][3][3] = {
-            {
-                { G111, G112, G113 },
-                { G112, G122, G123 },
-                { G113, G123, G133 },
-            },
-            {
-                { G211, G212, G213 },
-                { G212, G222, G223 },
-                { G213, G223, G233 },
-            },
-            {
-                { G311, G312, G313 },
-                { G312, G322, G323 },
-                { G313, G323, G333 },
-            }
-        };
+    //    const double G[3][3][3] = {
+    //        {
+    //            { G111, G112, G113 },
+    //            { G112, G122, G123 },
+    //            { G113, G123, G133 },
+    //        },
+    //        {
+    //            { G211, G212, G213 },
+    //            { G212, G222, G223 },
+    //            { G213, G223, G233 },
+    //        },
+    //        {
+    //            { G311, G312, G313 },
+    //            { G312, G322, G323 },
+    //            { G313, G323, G333 },
+    //        }
+    //    };
 
-        const double dalpha[3] = { PDstandardNth1alpha, 0.0, PDstandardNth3alpha };
+    //    const double dalpha[3] = { PDstandardNth1alpha, 0.0, PDstandardNth3alpha };
 
-        const double dphi[3] = { PDstandardNth1phi, 0.0, PDstandardNth3phi };
-        const double d2phi[3][3] = {{ PDstandardNth11phi, PDstandardNth12phi, PDstandardNth13phi },
-                                    { PDstandardNth12phi, PDstandardNth22phi, PDstandardNth23phi },
-                                    { PDstandardNth13phi, PDstandardNth23phi, PDstandardNth33phi }};
+    //    const double dphi[3] = { PDstandardNth1phi, 0.0, PDstandardNth3phi };
+    //    const double d2phi[3][3] = {{ PDstandardNth11phi, PDstandardNth12phi, PDstandardNth13phi },
+    //                                { PDstandardNth12phi, PDstandardNth22phi, PDstandardNth23phi },
+    //                                { PDstandardNth13phi, PDstandardNth23phi, PDstandardNth33phi }};
 
-        double gu[3][3], Ricm[3][3], dg[3][3][3], dgu[3][3][3], d2g[3][3][3][3], dG[3][3][3][3];
-        double d2beta_term[3] = { 0.0 };
+    //    double gu[3][3], Ricm[3][3], dg[3][3][3], dgu[3][3][3], d2g[3][3][3][3], dG[3][3][3][3];
+    //    double d2beta_term[3] = { 0.0 };
 
-        double ric_term;
+    //    double ric_term;
 
-        for (int i = 0; i < 3; i++)
-            for (int j = 0; j < 3; j++)
-                gu[i][j] = phiL * phiL * gtu[i][j];
+    //    for (int i = 0; i < 3; i++)
+    //        for (int j = 0; j < 3; j++)
+    //            gu[i][j] = phiL * phiL * gtu[i][j];
 
 
-        for (int i = 0; i < 3; i++)
-            for (int j = 0; j < 3; j++) {
-                double val = 0.0;
-                for (int k = 0; k < 3; k++)
-                    val += gu[i][k] * Ric[k][j];
-                Ricm[i][j] = val;
-            }
+    //    for (int i = 0; i < 3; i++)
+    //        for (int j = 0; j < 3; j++) {
+    //            double val = 0.0;
+    //            for (int k = 0; k < 3; k++)
+    //                val += gu[i][k] * Ric[k][j];
+    //            Ricm[i][j] = val;
+    //        }
 
-        for (int i = 0; i < 3; i++)
-            for (int j = 0; j < 3; j++)
-                for (int k = 0; k < 3; k++) {
-                    dg[j][k][i] = -2.0 * dphi[i] * gt[j][k] / (phiL * SQR(phiL)) + dgt[j][k][i] / SQR(phiL);
-                }
+    //    for (int i = 0; i < 3; i++)
+    //        for (int j = 0; j < 3; j++)
+    //            for (int k = 0; k < 3; k++) {
+    //                dg[j][k][i] = -2.0 * dphi[i] * gt[j][k] / (phiL * SQR(phiL)) + dgt[j][k][i] / SQR(phiL);
+    //            }
 
-        for (int i = 0; i < 3; i++)
-            for (int k = 0; k < 3; k++)
-                for (int m = 0; m < 3; m++) {
-                    double val = 0.0;
-                    for (int j = 0; j < 3; j++)
-                        for (int l = 0; l < 3; l++)
-                            val += - gu[k][j] * gu[l][m] * dg[j][l][i];
-                    dgu[k][m][i] = val;
-                }
+    //    for (int i = 0; i < 3; i++)
+    //        for (int k = 0; k < 3; k++)
+    //            for (int m = 0; m < 3; m++) {
+    //                double val = 0.0;
+    //                for (int j = 0; j < 3; j++)
+    //                    for (int l = 0; l < 3; l++)
+    //                        val += - gu[k][j] * gu[l][m] * dg[j][l][i];
+    //                dgu[k][m][i] = val;
+    //            }
 
-        for (int i = 0; i < 3; i++)
-            for (int j = 0; j < 3; j++)
-                for (int k = 0; k < 3; k++)
-                    for (int l = 0; l < 3; l++) {
-                        d2g[k][l][i][j] = 6.0 * gt[k][l] * dphi[i] * dphi[j] / SQR(SQR(phiL)) - 2.0 * gt[k][l] * d2phi[i][j] / (phiL * SQR(phiL))
-                                          - 2.0 * dgt[k][l][i] * dphi[j] / (phiL * SQR(phiL)) - 2.0 * dgt[k][l][j] * dphi[i] / (phiL * SQR(phiL)) + d2gt[k][l][i][j] / SQR(phiL);
-                    }
+    //    for (int i = 0; i < 3; i++)
+    //        for (int j = 0; j < 3; j++)
+    //            for (int k = 0; k < 3; k++)
+    //                for (int l = 0; l < 3; l++) {
+    //                    d2g[k][l][i][j] = 6.0 * gt[k][l] * dphi[i] * dphi[j] / SQR(SQR(phiL)) - 2.0 * gt[k][l] * d2phi[i][j] / (phiL * SQR(phiL))
+    //                                      - 2.0 * dgt[k][l][i] * dphi[j] / (phiL * SQR(phiL)) - 2.0 * dgt[k][l][j] * dphi[i] / (phiL * SQR(phiL)) + d2gt[k][l][i][j] / SQR(phiL);
+    //                }
 
-        for (int i = 0; i < 3; i++)
-            for (int j = 0; j < 3; j++)
-                for (int k = 0; k < 3; k++)
-                    for (int l = 0; l < 3; l++) {
-                        double val = 0.0;
-                        for (int m = 0; m < 3; m++) {
-                            val += dgu[k][m][i] * (dg[l][m][j] + dg[j][m][l] - dg[j][l][m]) + gu[k][m] * (d2g[l][m][i][j] + d2g[j][m][l][i] - d2g[j][l][m][i]);
-                        }
-                        dG[k][j][l][i] = 0.5 * val;
-                    }
+    //    for (int i = 0; i < 3; i++)
+    //        for (int j = 0; j < 3; j++)
+    //            for (int k = 0; k < 3; k++)
+    //                for (int l = 0; l < 3; l++) {
+    //                    double val = 0.0;
+    //                    for (int m = 0; m < 3; m++) {
+    //                        val += dgu[k][m][i] * (dg[l][m][j] + dg[j][m][l] - dg[j][l][m]) + gu[k][m] * (d2g[l][m][i][j] + d2g[j][m][l][i] - d2g[j][l][m][i]);
+    //                    }
+    //                    dG[k][j][l][i] = 0.5 * val;
+    //                }
 
-        for (int i = 0; i < 3; i++)
-            for (int j = 0; j < 3; j++)
-                for (int k = 0; k < 3; k++) {
-                    d2beta_term[k] += gu[i][j] * d2beta[k][i][j];
+    //    for (int i = 0; i < 3; i++)
+    //        for (int j = 0; j < 3; j++)
+    //            for (int k = 0; k < 3; k++) {
+    //                d2beta_term[k] += gu[i][j] * d2beta[k][i][j];
 
-                    for (int l = 0; l < 3; l++) {
-                        d2beta_term[k] += gu[i][j] * (G[k][j][l] * dbeta[l][i] + dG[k][j][l][i] * betaval[l] + G[k][i][l] * dbeta[l][j] - G[l][i][j] * dbeta[k][l]);
+    //                for (int l = 0; l < 3; l++) {
+    //                    d2beta_term[k] += gu[i][j] * (G[k][j][l] * dbeta[l][i] + dG[k][j][l][i] * betaval[l] + G[k][i][l] * dbeta[l][j] - G[l][i][j] * dbeta[k][l]);
 
-                        for (int m = 0; m < 3; m++) {
-                            d2beta_term[k] += gu[i][j] * (G[k][i][l] * G[l][j][m] * betaval[m] - G[l][i][j] * G[k][l][m] * betaval[m]);
-                        }
-                    }
-                }
+    //                    for (int m = 0; m < 3; m++) {
+    //                        d2beta_term[k] += gu[i][j] * (G[k][i][l] * G[l][j][m] * betaval[m] - G[l][i][j] * G[k][l][m] * betaval[m]);
+    //                    }
+    //                }
+    //            }
 
-        ric_term = 0.0;
-        for (int i = 0; i < 3; i++)
-            for (int j = 0; j < 3; j++)
-                ric_term += betaval[j] * Ricm[i][j] * dalpha[i];
+    //    ric_term = 0.0;
+    //    for (int i = 0; i < 3; i++)
+    //        for (int j = 0; j < 3; j++)
+    //            ric_term += betaval[j] * Ricm[i][j] * dalpha[i];
 
-        beta_term_val = -(d2beta_term[0] * dalpha[0] + d2beta_term[2] * dalpha[2]) - ric_term;
-    }
+    //    beta_term_val = -(d2beta_term[0] * dalpha[0] + d2beta_term[2] * dalpha[2]) - ric_term;
+    //}
     
     /* Copy local copies back to grid functions */
     vec_store_partial_prepare(i,vecimin,vecimax);
@@ -2606,7 +2606,7 @@ static void eval_timediff_body(const cGH* restrict const cctkGH, const int dir,
     vec_store_nta_partial(Kdot22[index],Kdot22L);
     vec_store_nta_partial(Kdot23[index],Kdot23L);
     vec_store_nta_partial(Kdot33[index],Kdot33L);
-    vec_store_nta_partial(beta_term[index], beta_term_val);
+    //vec_store_nta_partial(beta_term[index], beta_term_val);
     vec_store_nta_partial(Xtdot1[index],Xtdot1L);
     vec_store_nta_partial(Xtdot2[index],Xtdot2L);
     vec_store_nta_partial(Xtdot3[index],Xtdot3L);
diff --git a/cartoon_constr.c b/cartoon_constr.c
index 91b3445..6f87960 100644
--- a/cartoon_constr.c
+++ b/cartoon_constr.c
@@ -44,7 +44,7 @@
     PDstandardNth23gt33 = kzero;
 
     PDstandardNth2phi  = kzero;
-    PDstandardNth22phi = kifthen(origin, PDstandardNth11phi, kmul(PDstandardNth1phi, vinx));
+    PDstandardNth22phi = kifthen(origin, PDstandardNth11phi, kmul(PDstandardNth1phi, xinv));
     PDstandardNth12phi = kzero;
     PDstandardNth23phi = kzero;
     PDstandardNth2Xt1  = kzero;