Move RHS evaluations around to optimize for 32Kb instruction cache machines.

Move RHS evaluations for trK and the lapse to RHS1 to better even out the
instruction count between RHS1 and RHS2 leading to better performance on
32Kb instruction cache machines. Also replace covariant derivatives with
respect to the physical metric with covariant derivatives with respect to
the conformal metric. This avoids having to compute the physical Christoffel
symbols. The performance on 64Kb instruction cache machines is unchanged.
Also get rid of some unused shorthands and generate code for 3 timelevels.

Signed-off-by: Peter Diener <diener@cct.lsu.edu>

5 files changed, 90 insertions, 195 deletions

diff --git a/ML_BSSN/interface.ccl b/ML_BSSN/interface.ccl
index e257aee..0474dc3 100644
--- a/ML_BSSN/interface.ccl
+++ b/ML_BSSN/interface.ccl
@@ -59,7 +59,7 @@ CCTK_REAL ML_mom type=GF timelevels=1 tags='tensortypealias="D" tensorweight=1.0
 } "ML_mom"
 
 public:
-CCTK_REAL ML_curv type=GF timelevels=4 tags='tensortypealias="DD_sym" tensorweight=-0.66666666666666666667'
+CCTK_REAL ML_curv type=GF timelevels=3 tags='tensortypealias="DD_sym" tensorweight=-0.66666666666666666667'
 {
   At11,
   At12,
@@ -70,13 +70,13 @@ CCTK_REAL ML_curv type=GF timelevels=4 tags='tensortypealias="DD_sym" tensorweig
 } "ML_curv"
 
 public:
-CCTK_REAL ML_dtlapse type=GF timelevels=4 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000'
+CCTK_REAL ML_dtlapse type=GF timelevels=3 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000'
 {
   A
 } "ML_dtlapse"
 
 public:
-CCTK_REAL ML_dtshift type=GF timelevels=4 tags='tensortypealias="U" tensorweight=1.0000000000000000000'
+CCTK_REAL ML_dtshift type=GF timelevels=3 tags='tensortypealias="U" tensorweight=1.0000000000000000000'
 {
   B1,
   B2,
@@ -84,7 +84,7 @@ CCTK_REAL ML_dtshift type=GF timelevels=4 tags='tensortypealias="U" tensorweight
 } "ML_dtshift"
 
 public:
-CCTK_REAL ML_Gamma type=GF timelevels=4 tags='tensortypealias="U" tensorweight=0.66666666666666666667'
+CCTK_REAL ML_Gamma type=GF timelevels=3 tags='tensortypealias="U" tensorweight=0.66666666666666666667'
 {
   Xt1,
   Xt2,
@@ -92,19 +92,19 @@ CCTK_REAL ML_Gamma type=GF timelevels=4 tags='tensortypealias="U" tensorweight=0
 } "ML_Gamma"
 
 public:
-CCTK_REAL ML_lapse type=GF timelevels=4 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000'
+CCTK_REAL ML_lapse type=GF timelevels=3 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000'
 {
   alpha
 } "ML_lapse"
 
 public:
-CCTK_REAL ML_log_confac type=GF timelevels=4 tags='tensortypealias="Scalar" tensorweight=0.16666666666666666667'
+CCTK_REAL ML_log_confac type=GF timelevels=3 tags='tensortypealias="Scalar" tensorweight=0.16666666666666666667'
 {
   phi
 } "ML_log_confac"
 
 public:
-CCTK_REAL ML_metric type=GF timelevels=4 tags='tensortypealias="DD_sym" tensorweight=-0.66666666666666666667'
+CCTK_REAL ML_metric type=GF timelevels=3 tags='tensortypealias="DD_sym" tensorweight=-0.66666666666666666667'
 {
   gt11,
   gt12,
@@ -115,7 +115,7 @@ CCTK_REAL ML_metric type=GF timelevels=4 tags='tensortypealias="DD_sym" tensorwe
 } "ML_metric"
 
 public:
-CCTK_REAL ML_shift type=GF timelevels=4 tags='tensortypealias="U" tensorweight=1.0000000000000000000'
+CCTK_REAL ML_shift type=GF timelevels=3 tags='tensortypealias="U" tensorweight=1.0000000000000000000'
 {
   beta1,
   beta2,
@@ -123,13 +123,13 @@ CCTK_REAL ML_shift type=GF timelevels=4 tags='tensortypealias="U" tensorweight=1
 } "ML_shift"
 
 public:
-CCTK_REAL ML_trace_curv type=GF timelevels=4 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000'
+CCTK_REAL ML_trace_curv type=GF timelevels=3 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000'
 {
   trK
 } "ML_trace_curv"
 
 public:
-CCTK_REAL ML_curvrhs type=GF timelevels=4 tags='tensortypealias="DD_sym" tensorweight=-0.66666666666666666667'
+CCTK_REAL ML_curvrhs type=GF timelevels=3 tags='tensortypealias="DD_sym" tensorweight=-0.66666666666666666667'
 {
   At11rhs,
   At12rhs,
@@ -140,13 +140,13 @@ CCTK_REAL ML_curvrhs type=GF timelevels=4 tags='tensortypealias="DD_sym" tensorw
 } "ML_curvrhs"
 
 public:
-CCTK_REAL ML_dtlapserhs type=GF timelevels=4 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000'
+CCTK_REAL ML_dtlapserhs type=GF timelevels=3 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000'
 {
   Arhs
 } "ML_dtlapserhs"
 
 public:
-CCTK_REAL ML_dtshiftrhs type=GF timelevels=4 tags='tensortypealias="U" tensorweight=1.0000000000000000000'
+CCTK_REAL ML_dtshiftrhs type=GF timelevels=3 tags='tensortypealias="U" tensorweight=1.0000000000000000000'
 {
   B1rhs,
   B2rhs,
@@ -154,7 +154,7 @@ CCTK_REAL ML_dtshiftrhs type=GF timelevels=4 tags='tensortypealias="U" tensorwei
 } "ML_dtshiftrhs"
 
 public:
-CCTK_REAL ML_Gammarhs type=GF timelevels=4 tags='tensortypealias="U" tensorweight=0.66666666666666666667'
+CCTK_REAL ML_Gammarhs type=GF timelevels=3 tags='tensortypealias="U" tensorweight=0.66666666666666666667'
 {
   Xt1rhs,
   Xt2rhs,
@@ -162,19 +162,19 @@ CCTK_REAL ML_Gammarhs type=GF timelevels=4 tags='tensortypealias="U" tensorweigh
 } "ML_Gammarhs"
 
 public:
-CCTK_REAL ML_lapserhs type=GF timelevels=4 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000'
+CCTK_REAL ML_lapserhs type=GF timelevels=3 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000'
 {
   alpharhs
 } "ML_lapserhs"
 
 public:
-CCTK_REAL ML_log_confacrhs type=GF timelevels=4 tags='tensortypealias="Scalar" tensorweight=0.16666666666666666667'
+CCTK_REAL ML_log_confacrhs type=GF timelevels=3 tags='tensortypealias="Scalar" tensorweight=0.16666666666666666667'
 {
   phirhs
 } "ML_log_confacrhs"
 
 public:
-CCTK_REAL ML_metricrhs type=GF timelevels=4 tags='tensortypealias="DD_sym" tensorweight=-0.66666666666666666667'
+CCTK_REAL ML_metricrhs type=GF timelevels=3 tags='tensortypealias="DD_sym" tensorweight=-0.66666666666666666667'
 {
   gt11rhs,
   gt12rhs,
@@ -185,7 +185,7 @@ CCTK_REAL ML_metricrhs type=GF timelevels=4 tags='tensortypealias="DD_sym" tenso
 } "ML_metricrhs"
 
 public:
-CCTK_REAL ML_shiftrhs type=GF timelevels=4 tags='tensortypealias="U" tensorweight=1.0000000000000000000'
+CCTK_REAL ML_shiftrhs type=GF timelevels=3 tags='tensortypealias="U" tensorweight=1.0000000000000000000'
 {
   beta1rhs,
   beta2rhs,
@@ -193,7 +193,7 @@ CCTK_REAL ML_shiftrhs type=GF timelevels=4 tags='tensortypealias="U" tensorweigh
 } "ML_shiftrhs"
 
 public:
-CCTK_REAL ML_trace_curvrhs type=GF timelevels=4 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000'
+CCTK_REAL ML_trace_curvrhs type=GF timelevels=3 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000'
 {
   trKrhs
 } "ML_trace_curvrhs"
diff --git a/ML_BSSN/param.ccl b/ML_BSSN/param.ccl
index 8fb0583..716c648 100644
--- a/ML_BSSN/param.ccl
+++ b/ML_BSSN/param.ccl
@@ -179,13 +179,13 @@ CCTK_INT ML_BSSN_MaxNumEvolvedVars "Number of evolved variables used by this tho
 restricted:
 CCTK_INT timelevels "Number of active timelevels"
 {
-  0:4 :: ""
+  0:3 :: ""
 } 3
 
 restricted:
 CCTK_INT rhs_timelevels "Number of active RHS timelevels"
 {
-  0:4 :: ""
+  0:3 :: ""
 } 1
 
 restricted:
diff --git a/ML_BSSN/schedule.ccl b/ML_BSSN/schedule.ccl
index 408a345..a7da3b0 100644
--- a/ML_BSSN/schedule.ccl
+++ b/ML_BSSN/schedule.ccl
@@ -23,10 +23,6 @@ if (timelevels == 3)
 {
   STORAGE: ML_curv[3]
 }
-if (timelevels == 4)
-{
-  STORAGE: ML_curv[4]
-}
 
 if (timelevels == 1)
 {
@@ -40,10 +36,6 @@ if (timelevels == 3)
 {
   STORAGE: ML_dtlapse[3]
 }
-if (timelevels == 4)
-{
-  STORAGE: ML_dtlapse[4]
-}
 
 if (timelevels == 1)
 {
@@ -57,10 +49,6 @@ if (timelevels == 3)
 {
   STORAGE: ML_dtshift[3]
 }
-if (timelevels == 4)
-{
-  STORAGE: ML_dtshift[4]
-}
 
 if (timelevels == 1)
 {
@@ -74,10 +62,6 @@ if (timelevels == 3)
 {
   STORAGE: ML_Gamma[3]
 }
-if (timelevels == 4)
-{
-  STORAGE: ML_Gamma[4]
-}
 
 if (timelevels == 1)
 {
@@ -91,10 +75,6 @@ if (timelevels == 3)
 {
   STORAGE: ML_lapse[3]
 }
-if (timelevels == 4)
-{
-  STORAGE: ML_lapse[4]
-}
 
 if (timelevels == 1)
 {
@@ -108,10 +88,6 @@ if (timelevels == 3)
 {
   STORAGE: ML_log_confac[3]
 }
-if (timelevels == 4)
-{
-  STORAGE: ML_log_confac[4]
-}
 
 if (timelevels == 1)
 {
@@ -125,10 +101,6 @@ if (timelevels == 3)
 {
   STORAGE: ML_metric[3]
 }
-if (timelevels == 4)
-{
-  STORAGE: ML_metric[4]
-}
 
 if (timelevels == 1)
 {
@@ -142,10 +114,6 @@ if (timelevels == 3)
 {
   STORAGE: ML_shift[3]
 }
-if (timelevels == 4)
-{
-  STORAGE: ML_shift[4]
-}
 
 if (timelevels == 1)
 {
@@ -159,10 +127,6 @@ if (timelevels == 3)
 {
   STORAGE: ML_trace_curv[3]
 }
-if (timelevels == 4)
-{
-  STORAGE: ML_trace_curv[4]
-}
 
 if (rhs_timelevels == 1)
 {
@@ -176,10 +140,6 @@ if (rhs_timelevels == 3)
 {
   STORAGE: ML_curvrhs[3]
 }
-if (rhs_timelevels == 4)
-{
-  STORAGE: ML_curvrhs[4]
-}
 
 if (rhs_timelevels == 1)
 {
@@ -193,10 +153,6 @@ if (rhs_timelevels == 3)
 {
   STORAGE: ML_dtlapserhs[3]
 }
-if (rhs_timelevels == 4)
-{
-  STORAGE: ML_dtlapserhs[4]
-}
 
 if (rhs_timelevels == 1)
 {
@@ -210,10 +166,6 @@ if (rhs_timelevels == 3)
 {
   STORAGE: ML_dtshiftrhs[3]
 }
-if (rhs_timelevels == 4)
-{
-  STORAGE: ML_dtshiftrhs[4]
-}
 
 if (rhs_timelevels == 1)
 {
@@ -227,10 +179,6 @@ if (rhs_timelevels == 3)
 {
   STORAGE: ML_Gammarhs[3]
 }
-if (rhs_timelevels == 4)
-{
-  STORAGE: ML_Gammarhs[4]
-}
 
 if (rhs_timelevels == 1)
 {
@@ -244,10 +192,6 @@ if (rhs_timelevels == 3)
 {
   STORAGE: ML_lapserhs[3]
 }
-if (rhs_timelevels == 4)
-{
-  STORAGE: ML_lapserhs[4]
-}
 
 if (rhs_timelevels == 1)
 {
@@ -261,10 +205,6 @@ if (rhs_timelevels == 3)
 {
   STORAGE: ML_log_confacrhs[3]
 }
-if (rhs_timelevels == 4)
-{
-  STORAGE: ML_log_confacrhs[4]
-}
 
 if (rhs_timelevels == 1)
 {
@@ -278,10 +218,6 @@ if (rhs_timelevels == 3)
 {
   STORAGE: ML_metricrhs[3]
 }
-if (rhs_timelevels == 4)
-{
-  STORAGE: ML_metricrhs[4]
-}
 
 if (rhs_timelevels == 1)
 {
@@ -295,10 +231,6 @@ if (rhs_timelevels == 3)
 {
   STORAGE: ML_shiftrhs[3]
 }
-if (rhs_timelevels == 4)
-{
-  STORAGE: ML_shiftrhs[4]
-}
 
 if (rhs_timelevels == 1)
 {
@@ -312,10 +244,6 @@ if (rhs_timelevels == 3)
 {
   STORAGE: ML_trace_curvrhs[3]
 }
-if (rhs_timelevels == 4)
-{
-  STORAGE: ML_trace_curvrhs[4]
-}
 
 schedule ML_BSSN_Startup at STARTUP
 {
diff --git a/ML_BSSN/src/ML_BSSN_RHS1.c b/ML_BSSN/src/ML_BSSN_RHS1.c
index 98f904a..c259447 100644
--- a/ML_BSSN/src/ML_BSSN_RHS1.c
+++ b/ML_BSSN/src/ML_BSSN_RHS1.c
@@ -85,6 +85,12 @@ void ML_BSSN_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int con
     // CCTK_REAL PDstandardNth1alpha = INITVALUE;
     // CCTK_REAL PDstandardNth2alpha = INITVALUE;
     // CCTK_REAL PDstandardNth3alpha = INITVALUE;
+    // CCTK_REAL PDstandardNth11alpha = INITVALUE;
+    // CCTK_REAL PDstandardNth22alpha = INITVALUE;
+    // CCTK_REAL PDstandardNth33alpha = INITVALUE;
+    // CCTK_REAL PDstandardNth12alpha = INITVALUE;
+    // CCTK_REAL PDstandardNth13alpha = INITVALUE;
+    // CCTK_REAL PDstandardNth23alpha = INITVALUE;
     // CCTK_REAL PDstandardNth1beta1 = INITVALUE;
     // CCTK_REAL PDstandardNth2beta1 = INITVALUE;
     // CCTK_REAL PDstandardNth3beta1 = INITVALUE;
@@ -138,6 +144,7 @@ void ML_BSSN_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int con
     // CCTK_REAL PDstandardNth3trK = INITVALUE;
     
     /* Assign local copies of grid functions */
+    CCTK_REAL  AL = A[index];
     CCTK_REAL  alphaL = alpha[index];
     CCTK_REAL  At11L = At11[index];
     CCTK_REAL  At12L = At12[index];
@@ -151,6 +158,7 @@ void ML_BSSN_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int con
     CCTK_REAL  beta1L = beta1[index];
     CCTK_REAL  beta2L = beta2[index];
     CCTK_REAL  beta3L = beta3[index];
+    CCTK_REAL  eTttL = (*stress_energy_state) ? (eTtt[index]) : 0.0;
     CCTK_REAL  eTtxL = (*stress_energy_state) ? (eTtx[index]) : 0.0;
     CCTK_REAL  eTtyL = (*stress_energy_state) ? (eTty[index]) : 0.0;
     CCTK_REAL  eTtzL = (*stress_energy_state) ? (eTtz[index]) : 0.0;
@@ -169,6 +177,7 @@ void ML_BSSN_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int con
     CCTK_REAL  phiL = phi[index];
     CCTK_REAL  rL = r[index];
     CCTK_REAL  trKL = trK[index];
+    CCTK_REAL  trKrhsL = trKrhs[index];
     CCTK_REAL  Xt1L = Xt1[index];
     CCTK_REAL  Xt1rhsL = Xt1rhs[index];
     CCTK_REAL  Xt2L = Xt2[index];
@@ -182,6 +191,12 @@ void ML_BSSN_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int con
     CCTK_REAL const PDstandardNth1alpha = PDstandardNth1(alpha, i, j, k);
     CCTK_REAL const PDstandardNth2alpha = PDstandardNth2(alpha, i, j, k);
     CCTK_REAL const PDstandardNth3alpha = PDstandardNth3(alpha, i, j, k);
+    CCTK_REAL const PDstandardNth11alpha = PDstandardNth11(alpha, i, j, k);
+    CCTK_REAL const PDstandardNth22alpha = PDstandardNth22(alpha, i, j, k);
+    CCTK_REAL const PDstandardNth33alpha = PDstandardNth33(alpha, i, j, k);
+    CCTK_REAL const PDstandardNth12alpha = PDstandardNth12(alpha, i, j, k);
+    CCTK_REAL const PDstandardNth13alpha = PDstandardNth13(alpha, i, j, k);
+    CCTK_REAL const PDstandardNth23alpha = PDstandardNth23(alpha, i, j, k);
     CCTK_REAL const PDstandardNth1beta1 = PDstandardNth1(beta1, i, j, k);
     CCTK_REAL const PDstandardNth2beta1 = PDstandardNth2(beta1, i, j, k);
     CCTK_REAL const PDstandardNth3beta1 = PDstandardNth3(beta1, i, j, k);
@@ -374,6 +389,18 @@ void ML_BSSN_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int con
     
     CCTK_REAL Atu33 = Atm31*gtu31 + Atm32*gtu32 + Atm33*gtu33;
     
+    CCTK_REAL e4phi = IfThen(conformalMethod,pow(phiL,-2),exp(4*phiL));
+    
+    CCTK_REAL em4phi = INV(e4phi);
+    
+    CCTK_REAL rho = pow(alphaL,-2)*(eTttL - 2*(beta2L*eTtyL + 
+      beta3L*eTtzL) + 2*(beta1L*(-eTtxL + beta2L*eTxyL + beta3L*eTxzL) + 
+      beta2L*beta3L*eTyzL) + eTxxL*SQR(beta1L) + eTyyL*SQR(beta2L) + 
+      eTzzL*SQR(beta3L));
+    
+    CCTK_REAL trS = em4phi*(eTxxL*gtu11 + eTyyL*gtu22 + 2*(eTxyL*gtu21 + 
+      eTxzL*gtu31 + eTyzL*gtu32) + eTzzL*gtu33);
+    
     CCTK_REAL S1 = (-eTtxL + beta1L*eTxxL + beta2L*eTxyL + 
       beta3L*eTxzL)*INV(alphaL);
     
@@ -383,6 +410,20 @@ void ML_BSSN_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int con
     CCTK_REAL S3 = (-eTtzL + beta1L*eTxzL + beta2L*eTyzL + 
       beta3L*eTzzL)*INV(alphaL);
     
+    trKrhsL = PDupwindNth1(trK, i, j, k)*beta1L + PDupwindNth2(trK, i, 
+      j, k)*beta2L + PDupwindNth3(trK, i, j, k)*beta3L - 
+      em4phi*(gtu11*PDstandardNth11alpha + gtu22*PDstandardNth22alpha + 
+      gtu33*(PDstandardNth33alpha + 2*cdphi3*PDstandardNth3alpha) + 
+      2*(gtu21*PDstandardNth12alpha + gtu31*(PDstandardNth13alpha + 
+      cdphi1*PDstandardNth3alpha) + gtu32*(PDstandardNth23alpha + 
+      cdphi2*PDstandardNth3alpha)) + PDstandardNth1alpha*(2*(cdphi1*gtu11 + 
+      cdphi2*gtu21 + cdphi3*gtu31) - Xtn1) + 
+      PDstandardNth2alpha*(2*(cdphi1*gtu21 + cdphi2*gtu22 + cdphi3*gtu32) - 
+      Xtn2) - PDstandardNth3alpha*Xtn3) + alphaL*(2*(Atm12*Atm21 + 
+      Atm13*Atm31 + Atm23*Atm32) + 
+      12.56637061435917295385057353311801153679*(rho + trS) + SQR(Atm11) + 
+      SQR(Atm22) + SQR(Atm33) + kthird*SQR(trKL));
+    
     CCTK_REAL phirhsL = PDupwindNth1(phi, i, j, k)*beta1L + 
       PDupwindNth2(phi, i, j, k)*beta2L + PDupwindNth3(phi, i, j, 
       k)*beta3L + (PDstandardNth1beta1 + PDstandardNth2beta2 + 
@@ -486,6 +527,14 @@ void ML_BSSN_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int con
     CCTK_REAL eta = BetaDriver*IfThen(rL > 
       SpatialBetaDriverRadius,SpatialBetaDriverRadius*INV(rL),1);
     
+    CCTK_REAL alpharhsL = (PDupwindNth1(alpha, i, j, k)*beta1L + 
+      PDupwindNth2(alpha, i, j, k)*beta2L + PDupwindNth3(alpha, i, j, 
+      k)*beta3L)*LapseAdvectionCoeff + harmonicF*(AL*(-1 + 
+      LapseAdvectionCoeff) - LapseAdvectionCoeff*trKL)*pow(alphaL,harmonicN);
+    
+    CCTK_REAL ArhsL = (-1 + LapseAdvectionCoeff)*(AL*AlphaDriver - 
+      trKrhsL);
+    
     CCTK_REAL beta1rhsL = (PDupwindNth1(beta1, i, j, k)*beta1L + 
       PDupwindNth2(beta1, i, j, k)*beta2L + PDupwindNth3(beta1, i, j, 
       k)*beta3L)*ShiftAdvectionCoeff + B1L*ShiftGammaCoeff;
@@ -515,6 +564,8 @@ void ML_BSSN_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int con
     
     
     /* Copy local copies back to grid functions */
+    alpharhs[index] = alpharhsL;
+    Arhs[index] = ArhsL;
     B1rhs[index] = B1rhsL;
     B2rhs[index] = B2rhsL;
     B3rhs[index] = B3rhsL;
@@ -528,6 +579,7 @@ void ML_BSSN_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int con
     gt23rhs[index] = gt23rhsL;
     gt33rhs[index] = gt33rhsL;
     phirhs[index] = phirhsL;
+    trKrhs[index] = trKrhsL;
     Xt1rhs[index] = Xt1rhsL;
     Xt2rhs[index] = Xt2rhsL;
     Xt3rhs[index] = Xt3rhsL;
diff --git a/ML_BSSN/src/ML_BSSN_RHS2.c b/ML_BSSN/src/ML_BSSN_RHS2.c
index c73c3b9..7972cc7 100644
--- a/ML_BSSN/src/ML_BSSN_RHS2.c
+++ b/ML_BSSN/src/ML_BSSN_RHS2.c
@@ -174,7 +174,6 @@ void ML_BSSN_RHS2_Body(cGH const * restrict const cctkGH, int const dir, int con
     // CCTK_REAL PDstandardNth3Xt3 = INITVALUE;
     
     /* Assign local copies of grid functions */
-    CCTK_REAL  AL = A[index];
     CCTK_REAL  alphaL = alpha[index];
     CCTK_REAL  At11L = At11[index];
     CCTK_REAL  At12L = At12[index];
@@ -185,10 +184,6 @@ void ML_BSSN_RHS2_Body(cGH const * restrict const cctkGH, int const dir, int con
     CCTK_REAL  beta1L = beta1[index];
     CCTK_REAL  beta2L = beta2[index];
     CCTK_REAL  beta3L = beta3[index];
-    CCTK_REAL  eTttL = (*stress_energy_state) ? (eTtt[index]) : 0.0;
-    CCTK_REAL  eTtxL = (*stress_energy_state) ? (eTtx[index]) : 0.0;
-    CCTK_REAL  eTtyL = (*stress_energy_state) ? (eTty[index]) : 0.0;
-    CCTK_REAL  eTtzL = (*stress_energy_state) ? (eTtz[index]) : 0.0;
     CCTK_REAL  eTxxL = (*stress_energy_state) ? (eTxx[index]) : 0.0;
     CCTK_REAL  eTxyL = (*stress_energy_state) ? (eTxy[index]) : 0.0;
     CCTK_REAL  eTxzL = (*stress_energy_state) ? (eTxz[index]) : 0.0;
@@ -203,7 +198,6 @@ void ML_BSSN_RHS2_Body(cGH const * restrict const cctkGH, int const dir, int con
     CCTK_REAL  gt33L = gt33[index];
     CCTK_REAL  phiL = phi[index];
     CCTK_REAL  trKL = trK[index];
-    CCTK_REAL  trKrhsL = trKrhs[index];
     CCTK_REAL  Xt1L = Xt1[index];
     CCTK_REAL  Xt2L = Xt2[index];
     CCTK_REAL  Xt3L = Xt3[index];
@@ -652,60 +646,6 @@ void ML_BSSN_RHS2_Body(cGH const * restrict const cctkGH, int const dir, int con
     
     CCTK_REAL gu33 = em4phi*gtu33;
     
-    CCTK_REAL G111 = Gt111 + cdphi1*(4 - 2*gt11L*gtu11) - 
-      2*gt11L*(cdphi2*gtu21 + cdphi3*gtu31);
-    
-    CCTK_REAL G211 = Gt211 - 2*gt11L*(cdphi1*gtu21 + cdphi2*gtu22 + 
-      cdphi3*gtu32);
-    
-    CCTK_REAL G311 = Gt311 - 2*gt11L*(cdphi1*gtu31 + cdphi2*gtu32 + 
-      cdphi3*gtu33);
-    
-    CCTK_REAL G112 = Gt112 + cdphi2*(2 - 2*gt12L*gtu21) - 
-      2*gt12L*(cdphi1*gtu11 + cdphi3*gtu31);
-    
-    CCTK_REAL G212 = Gt212 + cdphi1*(2 - 2*gt12L*gtu21) - 
-      2*gt12L*(cdphi2*gtu22 + cdphi3*gtu32);
-    
-    CCTK_REAL G312 = Gt312 - 2*gt12L*(cdphi1*gtu31 + cdphi2*gtu32 + 
-      cdphi3*gtu33);
-    
-    CCTK_REAL G113 = Gt113 - 2*gt13L*(cdphi1*gtu11 + cdphi2*gtu21) + 
-      cdphi3*(2 - 2*gt13L*gtu31);
-    
-    CCTK_REAL G213 = Gt213 - 2*gt13L*(cdphi1*gtu21 + cdphi2*gtu22 + 
-      cdphi3*gtu32);
-    
-    CCTK_REAL G313 = Gt313 + cdphi1*(2 - 2*gt13L*gtu31) - 
-      2*gt13L*(cdphi2*gtu32 + cdphi3*gtu33);
-    
-    CCTK_REAL G122 = Gt122 - 2*gt22L*(cdphi1*gtu11 + cdphi2*gtu21 + 
-      cdphi3*gtu31);
-    
-    CCTK_REAL G222 = Gt222 + cdphi2*(4 - 2*gt22L*gtu22) - 
-      2*gt22L*(cdphi1*gtu21 + cdphi3*gtu32);
-    
-    CCTK_REAL G322 = Gt322 - 2*gt22L*(cdphi1*gtu31 + cdphi2*gtu32 + 
-      cdphi3*gtu33);
-    
-    CCTK_REAL G123 = Gt123 - 2*gt23L*(cdphi1*gtu11 + cdphi2*gtu21 + 
-      cdphi3*gtu31);
-    
-    CCTK_REAL G223 = Gt223 - 2*gt23L*(cdphi1*gtu21 + cdphi2*gtu22) + 
-      cdphi3*(2 - 2*gt23L*gtu32);
-    
-    CCTK_REAL G323 = Gt323 + cdphi2*(2 - 2*gt23L*gtu32) - 
-      2*gt23L*(cdphi1*gtu31 + cdphi3*gtu33);
-    
-    CCTK_REAL G133 = Gt133 - 2*gt33L*(cdphi1*gtu11 + cdphi2*gtu21 + 
-      cdphi3*gtu31);
-    
-    CCTK_REAL G233 = Gt233 - 2*gt33L*(cdphi1*gtu21 + cdphi2*gtu22 + 
-      cdphi3*gtu32);
-    
-    CCTK_REAL G333 = Gt333 - 2*gt33L*(cdphi1*gtu31 + cdphi2*gtu32) + 
-      cdphi3*(4 - 2*gt33L*gtu33);
-    
     CCTK_REAL R11 = Rphi11 + Rt11;
     
     CCTK_REAL R12 = Rphi12 + Rt12;
@@ -718,46 +658,32 @@ void ML_BSSN_RHS2_Body(cGH const * restrict const cctkGH, int const dir, int con
     
     CCTK_REAL R33 = Rphi33 + Rt33;
     
-    CCTK_REAL rho = pow(alphaL,-2)*(eTttL - 2*(beta2L*eTtyL + 
-      beta3L*eTtzL) + 2*(beta1L*(-eTtxL + beta2L*eTxyL + beta3L*eTxzL) + 
-      beta2L*beta3L*eTyzL) + eTxxL*SQR(beta1L) + eTyyL*SQR(beta2L) + 
-      eTzzL*SQR(beta3L));
-    
     CCTK_REAL trS = eTxxL*gu11 + eTyyL*gu22 + 2*(eTxyL*gu21 + eTxzL*gu31 + 
       eTyzL*gu32) + eTzzL*gu33;
     
-    trKrhsL = PDupwindNth1(trK, i, j, k)*beta1L + PDupwindNth2(trK, i, 
-      j, k)*beta2L + PDupwindNth3(trK, i, j, k)*beta3L + (G111*gu11 + 
-      G122*gu22 + 2.*(G112*gu21 + G113*gu31 + G123*gu32) + 
-      G133*gu33)*PDstandardNth1alpha - 2.*(gu21*PDstandardNth12alpha + 
-      gu31*PDstandardNth13alpha + gu32*PDstandardNth23alpha) + (G211*gu11 + 
-      G222*gu22 + 2.*(G212*gu21 + G213*gu31 + G223*gu32) + 
-      G233*gu33)*PDstandardNth2alpha - 1.*(gu11*PDstandardNth11alpha + 
-      gu22*PDstandardNth22alpha + gu33*PDstandardNth33alpha) + (G311*gu11 + 
-      G322*gu22 + 2.*(G313*gu31 + G323*gu32) + G333*gu33)*PDstandardNth3alpha 
-      + 2.*(alphaL*(Atm12*Atm21 + Atm13*Atm31 + Atm23*Atm32) + 
-      G312*gu21*PDstandardNth3alpha) + 
-      alphaL*(12.56637061435917295385057353311801153679*(rho + trS) + 
-      SQR(Atm11) + SQR(Atm22) + SQR(Atm33) + 
-      0.3333333333333333333333333333333333333333*SQR(trKL));
-    
-    CCTK_REAL Ats11 = -PDstandardNth11alpha + G111*PDstandardNth1alpha + 
-      G211*PDstandardNth2alpha + G311*PDstandardNth3alpha + alphaL*R11;
+    CCTK_REAL Ats11 = -PDstandardNth11alpha + (4*cdphi1 + 
+      Gt111)*PDstandardNth1alpha + Gt211*PDstandardNth2alpha + 
+      Gt311*PDstandardNth3alpha + alphaL*R11;
     
-    CCTK_REAL Ats12 = -PDstandardNth12alpha + G112*PDstandardNth1alpha + 
-      G212*PDstandardNth2alpha + G312*PDstandardNth3alpha + alphaL*R12;
+    CCTK_REAL Ats12 = -PDstandardNth12alpha + (2*cdphi2 + 
+      Gt112)*PDstandardNth1alpha + (2*cdphi1 + Gt212)*PDstandardNth2alpha + 
+      Gt312*PDstandardNth3alpha + alphaL*R12;
     
-    CCTK_REAL Ats13 = -PDstandardNth13alpha + G113*PDstandardNth1alpha + 
-      G213*PDstandardNth2alpha + G313*PDstandardNth3alpha + alphaL*R13;
+    CCTK_REAL Ats13 = -PDstandardNth13alpha + (2*cdphi3 + 
+      Gt113)*PDstandardNth1alpha + Gt213*PDstandardNth2alpha + (2*cdphi1 + 
+      Gt313)*PDstandardNth3alpha + alphaL*R13;
     
-    CCTK_REAL Ats22 = G122*PDstandardNth1alpha - PDstandardNth22alpha + 
-      G222*PDstandardNth2alpha + G322*PDstandardNth3alpha + alphaL*R22;
+    CCTK_REAL Ats22 = Gt122*PDstandardNth1alpha - PDstandardNth22alpha + 
+      (4*cdphi2 + Gt222)*PDstandardNth2alpha + Gt322*PDstandardNth3alpha + 
+      alphaL*R22;
     
-    CCTK_REAL Ats23 = G123*PDstandardNth1alpha - PDstandardNth23alpha + 
-      G223*PDstandardNth2alpha + G323*PDstandardNth3alpha + alphaL*R23;
+    CCTK_REAL Ats23 = Gt123*PDstandardNth1alpha - PDstandardNth23alpha + 
+      (2*cdphi3 + Gt223)*PDstandardNth2alpha + (2*cdphi2 + 
+      Gt323)*PDstandardNth3alpha + alphaL*R23;
     
-    CCTK_REAL Ats33 = G133*PDstandardNth1alpha + G233*PDstandardNth2alpha 
-      - PDstandardNth33alpha + G333*PDstandardNth3alpha + alphaL*R33;
+    CCTK_REAL Ats33 = Gt133*PDstandardNth1alpha + 
+      Gt233*PDstandardNth2alpha - PDstandardNth33alpha + (4*cdphi3 + 
+      Gt333)*PDstandardNth3alpha + alphaL*R33;
     
     CCTK_REAL trAts = Ats11*gu11 + Ats22*gu22 + 2*(Ats12*gu21 + Ats13*gu31 
       + Ats23*gu32) + Ats33*gu33;
@@ -834,25 +760,14 @@ void ML_BSSN_RHS2_Body(cGH const * restrict const cctkGH, int const dir, int con
       alphaL*(-25.13274122871834590770114706623602307358*eTzzL + 
       8.377580409572781969233715688745341024526*g33*trS));
     
-    CCTK_REAL alpharhsL = (PDupwindNth1(alpha, i, j, k)*beta1L + 
-      PDupwindNth2(alpha, i, j, k)*beta2L + PDupwindNth3(alpha, i, j, 
-      k)*beta3L)*LapseAdvectionCoeff + harmonicF*(AL*(-1 + 
-      LapseAdvectionCoeff) - LapseAdvectionCoeff*trKL)*pow(alphaL,harmonicN);
-    
-    CCTK_REAL ArhsL = (-1 + LapseAdvectionCoeff)*(AL*AlphaDriver - 
-      trKrhsL);
-    
     
     /* Copy local copies back to grid functions */
-    alpharhs[index] = alpharhsL;
-    Arhs[index] = ArhsL;
     At11rhs[index] = At11rhsL;
     At12rhs[index] = At12rhsL;
     At13rhs[index] = At13rhsL;
     At22rhs[index] = At22rhsL;
     At23rhs[index] = At23rhsL;
     At33rhs[index] = At33rhsL;
-    trKrhs[index] = trKrhsL;
   }
   LC_ENDLOOP3 (ML_BSSN_RHS2);
 }