Merge branch 'master' of /Users/eschnett/Cbeta/arrangements/McLachlan

Conflicts:
	ML_BSSN/src/ML_BSSN_RHS1.c
	ML_BSSN_MP/src/ML_BSSN_MP_RHS1.c
	ML_BSSN_O2/src/ML_BSSN_O2_RHS1.c
	ML_BSSN_O8/src/ML_BSSN_O8_RHS1.c
	m/McLachlan_BSSN.m

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 188e0f1..c0e80de 100644
--- a/ML_BSSN/param.ccl
+++ b/ML_BSSN/param.ccl
@@ -185,13 +185,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 932cdc1..1378a36 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 + 
@@ -489,6 +530,14 @@ void ML_BSSN_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int con
     CCTK_REAL theta = IfThen(rL > SpatialShiftGammaCoeffRadius,exp(1 - 
       rL*INV(SpatialShiftGammaCoeffRadius)),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*theta;
@@ -518,6 +567,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;
@@ -531,6 +582,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);
 }