Re-introduce the distinction between vacuum and matter spacetimes.

This is convenient since Kranc cannot easily detect whether the TmunuBase
have storage or not, thus requiring them to have storage at all times.
Differentiating between vacuum and matter spacetime avoids this.

2 files changed, 121 insertions, 273 deletions

diff --git a/ML_BSSN/src/ML_BSSN_RHS.c b/ML_BSSN/src/ML_BSSN_RHS.c
index 4d70d0a..dfa32e4 100644
--- a/ML_BSSN/src/ML_BSSN_RHS.c
+++ b/ML_BSSN/src/ML_BSSN_RHS.c
@@ -135,14 +135,9 @@ void ML_BSSN_RHS_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_INT con
     CCTK_REAL gtu11 = INITVALUE, gtu21 = INITVALUE, gtu22 = INITVALUE, gtu31 = INITVALUE, gtu32 = INITVALUE, gtu33 = INITVALUE;
     CCTK_REAL gu11 = INITVALUE, gu21 = INITVALUE, gu22 = INITVALUE, gu31 = INITVALUE, gu32 = INITVALUE, gu33 = INITVALUE;
     CCTK_REAL R11 = INITVALUE, R12 = INITVALUE, R13 = INITVALUE, R22 = INITVALUE, R23 = INITVALUE, R33 = INITVALUE;
-    CCTK_REAL rho = INITVALUE;
     CCTK_REAL Rphi11 = INITVALUE, Rphi12 = INITVALUE, Rphi13 = INITVALUE, Rphi22 = INITVALUE, Rphi23 = INITVALUE, Rphi33 = INITVALUE;
     CCTK_REAL Rt11 = INITVALUE, Rt12 = INITVALUE, Rt13 = INITVALUE, Rt22 = INITVALUE, Rt23 = INITVALUE, Rt33 = INITVALUE;
-    CCTK_REAL S1 = INITVALUE, S2 = INITVALUE, S3 = INITVALUE;
-    CCTK_REAL T00 = INITVALUE, T01 = INITVALUE, T02 = INITVALUE, T03 = INITVALUE, T11 = INITVALUE, T12 = INITVALUE;
-    CCTK_REAL T13 = INITVALUE, T22 = INITVALUE, T23 = INITVALUE, T33 = INITVALUE;
     CCTK_REAL trAts = INITVALUE;
-    CCTK_REAL trS = INITVALUE;
     CCTK_REAL Xtn1 = INITVALUE, Xtn2 = INITVALUE, Xtn3 = INITVALUE;
     
     /* Declare local copies of grid functions */
@@ -153,16 +148,6 @@ void ML_BSSN_RHS_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_INT con
     CCTK_REAL At22L = INITVALUE, At22rhsL = INITVALUE, At23L = INITVALUE, At23rhsL = INITVALUE, At33L = INITVALUE, At33rhsL = INITVALUE;
     CCTK_REAL B1L = INITVALUE, B1rhsL = INITVALUE, B2L = INITVALUE, B2rhsL = INITVALUE, B3L = INITVALUE, B3rhsL = INITVALUE;
     CCTK_REAL beta1L = INITVALUE, beta1rhsL = INITVALUE, beta2L = INITVALUE, beta2rhsL = INITVALUE, beta3L = INITVALUE, beta3rhsL = INITVALUE;
-    CCTK_REAL eTttL = INITVALUE;
-    CCTK_REAL eTtxL = INITVALUE;
-    CCTK_REAL eTtyL = INITVALUE;
-    CCTK_REAL eTtzL = INITVALUE;
-    CCTK_REAL eTxxL = INITVALUE;
-    CCTK_REAL eTxyL = INITVALUE;
-    CCTK_REAL eTxzL = INITVALUE;
-    CCTK_REAL eTyyL = INITVALUE;
-    CCTK_REAL eTyzL = INITVALUE;
-    CCTK_REAL eTzzL = INITVALUE;
     CCTK_REAL gt11L = INITVALUE, gt11rhsL = INITVALUE, gt12L = INITVALUE, gt12rhsL = INITVALUE, gt13L = INITVALUE, gt13rhsL = INITVALUE;
     CCTK_REAL gt22L = INITVALUE, gt22rhsL = INITVALUE, gt23L = INITVALUE, gt23rhsL = INITVALUE, gt33L = INITVALUE, gt33rhsL = INITVALUE;
     CCTK_REAL phiL = INITVALUE, phirhsL = INITVALUE;
@@ -358,16 +343,6 @@ void ML_BSSN_RHS_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_INT con
     beta1L = beta1[index];
     beta2L = beta2[index];
     beta3L = beta3[index];
-    eTttL = eTtt[index];
-    eTtxL = eTtx[index];
-    eTtyL = eTty[index];
-    eTtzL = eTtz[index];
-    eTxxL = eTxx[index];
-    eTxyL = eTxy[index];
-    eTxzL = eTxz[index];
-    eTyyL = eTyy[index];
-    eTyzL = eTyz[index];
-    eTzzL = eTzz[index];
     gt11L = gt11[index];
     gt12L = gt12[index];
     gt13L = gt13[index];
@@ -1276,38 +1251,6 @@ void ML_BSSN_RHS_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_INT con
     
     R33  =  Rphi33 + Rt33;
     
-    T00  =  eTttL;
-    
-    T01  =  eTtxL;
-    
-    T02  =  eTtyL;
-    
-    T03  =  eTtzL;
-    
-    T11  =  eTxxL;
-    
-    T12  =  eTxyL;
-    
-    T13  =  eTxzL;
-    
-    T22  =  eTyyL;
-    
-    T23  =  eTyzL;
-    
-    T33  =  eTzzL;
-    
-    rho  =  pow(alphaL,-2)*(T00 - 2*(beta2L*T02 + beta3L*T03) + 
-          2*(beta1L*(-T01 + beta2L*T12 + beta3L*T13) + beta2L*beta3L*T23) + 
-          T11*SQR(beta1L) + T22*SQR(beta2L) + T33*SQR(beta3L));
-    
-    S1  =  (-T01 + beta1L*T11 + beta2L*T12 + beta3L*T13)*INV(alphaL);
-    
-    S2  =  (-T02 + beta1L*T12 + beta2L*T22 + beta3L*T23)*INV(alphaL);
-    
-    S3  =  (-T03 + beta1L*T13 + beta2L*T23 + beta3L*T33)*INV(alphaL);
-    
-    trS  =  gu11*T11 + gu22*T22 + 2*(gu21*T12 + gu31*T13 + gu32*T23) + gu33*T33;
-    
     phirhsL  =  (PDstandardNth1beta1 + PDstandardNth2beta2 + PDstandardNth3beta3 + 
           6*(beta1L*PDstandardNth1phi + beta2L*PDstandardNth2phi + 
              beta3L*PDstandardNth3phi) - alphaL*trKL)/6.;
@@ -1370,9 +1313,7 @@ void ML_BSSN_RHS_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_INT con
                    gtu31*PDstandardNth3trK))) + 
           3*(beta1L*PDstandardNth1Xt1 + gtu22*PDstandardNth22beta1 + 
              beta2L*PDstandardNth2Xt1 + gtu33*PDstandardNth33beta1 + 
-             beta3L*PDstandardNth3Xt1) - 
-          150.7964473723100754462068823974161384415*alphaL*
-           (gtu11*S1 + gtu21*S2 + gtu31*S3) + 
+             beta3L*PDstandardNth3Xt1) + 
           (-3*PDstandardNth1beta1 + 2*
               (PDstandardNth1beta1 + PDstandardNth2beta2 + PDstandardNth3beta3))*
            Xtn1 - 3*(PDstandardNth2beta1*Xtn2 + PDstandardNth3beta1*Xtn3));
@@ -1394,9 +1335,7 @@ void ML_BSSN_RHS_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_INT con
                    gtu32*PDstandardNth3trK))) + 
           3*(gtu11*PDstandardNth11beta2 + beta1L*PDstandardNth1Xt2 + 
              beta2L*PDstandardNth2Xt2 + gtu33*PDstandardNth33beta2 + 
-             beta3L*PDstandardNth3Xt2) - 
-          150.7964473723100754462068823974161384415*alphaL*
-           (gtu21*S1 + gtu22*S2 + gtu32*S3) + 
+             beta3L*PDstandardNth3Xt2) + 
           2*(PDstandardNth1beta1 + PDstandardNth2beta2 + PDstandardNth3beta3)*
            Xtn2 - 3*(PDstandardNth1beta2*Xtn1 + PDstandardNth2beta2*Xtn2 + 
              PDstandardNth3beta2*Xtn3));
@@ -1418,28 +1357,25 @@ void ML_BSSN_RHS_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_INT con
                    gtu33*PDstandardNth3trK))) + 
           3*(gtu11*PDstandardNth11beta3 + beta1L*PDstandardNth1Xt3 + 
              gtu22*PDstandardNth22beta3 + beta2L*PDstandardNth2Xt3 + 
-             beta3L*PDstandardNth3Xt3) - 
-          150.7964473723100754462068823974161384415*alphaL*
-           (gtu31*S1 + gtu32*S2 + gtu33*S3) + 
+             beta3L*PDstandardNth3Xt3) + 
           2*(PDstandardNth1beta1 + PDstandardNth2beta2 + PDstandardNth3beta3)*
            Xtn3 - 3*(PDstandardNth1beta3*Xtn1 + PDstandardNth2beta3*Xtn2 + 
              PDstandardNth3beta3*Xtn3));
     
-    trKrhsL  =  (G111*gu11 + G122*gu22 + 2.*(G112*gu21 + G113*gu31 + G123*gu32) + 
+    trKrhsL  =  -(gu11*PDstandardNth11alpha) - 2*gu21*PDstandardNth12alpha - 
+        2*gu31*PDstandardNth13alpha + 
+        (G111*gu11 + 2*G112*gu21 + G122*gu22 + 2*G113*gu31 + 2*G123*gu32 + 
            G133*gu33)*PDstandardNth1alpha + beta1L*PDstandardNth1trK - 
-        2.*(gu21*PDstandardNth12alpha + gu31*PDstandardNth13alpha + 
-           gu32*PDstandardNth23alpha) + 
-        (G211*gu11 + G222*gu22 + 2.*(G212*gu21 + G213*gu31 + G223*gu32) + 
+        gu22*PDstandardNth22alpha - 2*gu32*PDstandardNth23alpha + 
+        (G211*gu11 + 2*G212*gu21 + G222*gu22 + 2*G213*gu31 + 2*G223*gu32 + 
            G233*gu33)*PDstandardNth2alpha + beta2L*PDstandardNth2trK - 
-        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) + 
+        gu33*PDstandardNth33alpha + G311*gu11*PDstandardNth3alpha + 
+        G322*gu22*PDstandardNth3alpha + 2*G313*gu31*PDstandardNth3alpha + 
+        2*G323*gu32*PDstandardNth3alpha + G333*gu33*PDstandardNth3alpha + 
+        2*(alphaL*Atm12*Atm21 + alphaL*Atm13*Atm31 + alphaL*Atm23*Atm32 + 
            G312*gu21*PDstandardNth3alpha) + beta3L*PDstandardNth3trK + 
-        alphaL*(12.56637061435917295385057353311801153679*(rho + trS) + 
-           SQR(Atm11) + SQR(Atm22) + SQR(Atm33) + 
-           0.3333333333333333333333333333333333333333*SQR(trKL));
+        alphaL*SQR(Atm11) + alphaL*SQR(Atm22) + alphaL*SQR(Atm33) + 
+        alphaL*kthird*SQR(trKL);
     
     Ats11  =  -PDstandardNth11alpha + G111*PDstandardNth1alpha + 
         G211*PDstandardNth2alpha + G311*PDstandardNth3alpha + alphaL*R11;
@@ -1462,76 +1398,51 @@ void ML_BSSN_RHS_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK_INT con
     trAts  =  Ats11*gu11 + Ats22*gu22 + 2*(Ats12*gu21 + Ats13*gu31 + Ats23*gu32) + 
         Ats33*gu33;
     
-    At11rhsL  =  -2.*alphaL*(At11L*Atm11 + At12L*Atm21 + At13L*Atm31) + 
-        beta1L*PDstandardNth1At11 + 2.*
-         (At12L*PDstandardNth1beta2 + At13L*PDstandardNth1beta3) + 
-        beta2L*PDstandardNth2At11 + beta3L*PDstandardNth3At11 + 
-        At11L*(1.333333333333333333333333333333333333333*PDstandardNth1beta1 - 
-           0.6666666666666666666666666666666666666667*
-            (PDstandardNth2beta2 + PDstandardNth3beta3) + alphaL*trKL) + 
-        em4phi*(Ats11 - 0.3333333333333333333333333333333333333333*g11*trAts + 
-           alphaL*(-25.13274122871834590770114706623602307358*T11 + 
-              8.377580409572781969233715688745341024526*g11*trS));
-    
-    At12rhsL  =  -2.*alphaL*(At11L*Atm12 + At12L*Atm22 + At13L*Atm32) + 
+    At11rhsL  =  beta1L*PDstandardNth1At11 + 
+        2*(At11L*PDstandardNth1beta1 + At12L*PDstandardNth1beta2 + 
+           At13L*PDstandardNth1beta3) + beta2L*PDstandardNth2At11 + 
+        beta3L*PDstandardNth3At11 - At11L*ktwothird*
+         (PDstandardNth1beta1 + PDstandardNth2beta2 + PDstandardNth3beta3) + 
+        em4phi*(Ats11 - g11*kthird*trAts) + 
+        alphaL*(-2*(At11L*Atm11 + At12L*Atm21 + At13L*Atm31) + At11L*trKL);
+    
+    At12rhsL  =  -2*alphaL*(At11L*Atm12 + At12L*Atm22 + At13L*Atm32) + 
         beta1L*PDstandardNth1At12 + At22L*PDstandardNth1beta2 + 
         At23L*PDstandardNth1beta3 + beta2L*PDstandardNth2At12 + 
         At11L*PDstandardNth2beta1 + At13L*PDstandardNth2beta3 + 
-        beta3L*PDstandardNth3At12 + At12L*
-         (0.3333333333333333333333333333333333333333*
-            (PDstandardNth1beta1 + PDstandardNth2beta2) - 
-           0.6666666666666666666666666666666666666667*PDstandardNth3beta3 + 
-           alphaL*trKL) + em4phi*(Ats12 - 
-           0.3333333333333333333333333333333333333333*g12*trAts + 
-           alphaL*(-25.13274122871834590770114706623602307358*T12 + 
-              8.377580409572781969233715688745341024526*g12*trS));
-    
-    At13rhsL  =  -2.*alphaL*(At11L*Atm13 + At12L*Atm23 + At13L*Atm33) + 
+        beta3L*PDstandardNth3At12 + em4phi*(Ats12 - g12*kthird*trAts) + 
+        At12L*(kthird*(PDstandardNth1beta1 + PDstandardNth2beta2) - 
+           ktwothird*PDstandardNth3beta3 + alphaL*trKL);
+    
+    At13rhsL  =  -2*alphaL*(At11L*Atm13 + At12L*Atm23 + At13L*Atm33) + 
         beta1L*PDstandardNth1At13 + At23L*PDstandardNth1beta2 + 
         At33L*PDstandardNth1beta3 + beta2L*PDstandardNth2At13 + 
         beta3L*PDstandardNth3At13 + At11L*PDstandardNth3beta1 + 
-        At12L*PDstandardNth3beta2 + At13L*
-         (-0.6666666666666666666666666666666666666667*PDstandardNth2beta2 + 
-           0.3333333333333333333333333333333333333333*
-            (PDstandardNth1beta1 + PDstandardNth3beta3) + alphaL*trKL) + 
-        em4phi*(Ats13 - 0.3333333333333333333333333333333333333333*g13*trAts + 
-           alphaL*(-25.13274122871834590770114706623602307358*T13 + 
-              8.377580409572781969233715688745341024526*g13*trS));
-    
-    At22rhsL  =  -2.*alphaL*(At12L*Atm12 + At22L*Atm22 + At23L*Atm32) + 
-        beta1L*PDstandardNth1At22 + beta2L*PDstandardNth2At22 + 
-        2.*(At12L*PDstandardNth2beta1 + At23L*PDstandardNth2beta3) + 
-        beta3L*PDstandardNth3At22 + At22L*
-         (1.333333333333333333333333333333333333333*PDstandardNth2beta2 - 
-           0.6666666666666666666666666666666666666667*
-            (PDstandardNth1beta1 + PDstandardNth3beta3) + alphaL*trKL) + 
-        em4phi*(Ats22 - 0.3333333333333333333333333333333333333333*g22*trAts + 
-           alphaL*(-25.13274122871834590770114706623602307358*T22 + 
-              8.377580409572781969233715688745341024526*g22*trS));
-    
-    At23rhsL  =  -2.*alphaL*(At12L*Atm13 + At22L*Atm23 + At23L*Atm33) + 
+        At12L*PDstandardNth3beta2 + em4phi*(Ats13 - g13*kthird*trAts) + 
+        At13L*(-(ktwothird*PDstandardNth2beta2) + 
+           kthird*(PDstandardNth1beta1 + PDstandardNth3beta3) + alphaL*trKL);
+    
+    At22rhsL  =  beta1L*PDstandardNth1At22 + beta2L*PDstandardNth2At22 + 
+        2*(At12L*PDstandardNth2beta1 + At22L*PDstandardNth2beta2 + 
+           At23L*PDstandardNth2beta3) + beta3L*PDstandardNth3At22 - 
+        At22L*ktwothird*(PDstandardNth1beta1 + PDstandardNth2beta2 + 
+           PDstandardNth3beta3) + em4phi*(Ats22 - g22*kthird*trAts) + 
+        alphaL*(-2*(At12L*Atm12 + At22L*Atm22 + At23L*Atm32) + At22L*trKL);
+    
+    At23rhsL  =  -2*alphaL*(At12L*Atm13 + At22L*Atm23 + At23L*Atm33) + 
         beta1L*PDstandardNth1At23 + beta2L*PDstandardNth2At23 + 
         At13L*PDstandardNth2beta1 + At33L*PDstandardNth2beta3 + 
         beta3L*PDstandardNth3At23 + At12L*PDstandardNth3beta1 + 
-        At22L*PDstandardNth3beta2 + At23L*
-         (-0.6666666666666666666666666666666666666667*PDstandardNth1beta1 + 
-           0.3333333333333333333333333333333333333333*
-            (PDstandardNth2beta2 + PDstandardNth3beta3) + alphaL*trKL) + 
-        em4phi*(Ats23 - 0.3333333333333333333333333333333333333333*g23*trAts + 
-           alphaL*(-25.13274122871834590770114706623602307358*T23 + 
-              8.377580409572781969233715688745341024526*g23*trS));
-    
-    At33rhsL  =  -2.*alphaL*(At13L*Atm13 + At23L*Atm23 + At33L*Atm33) + 
-        beta1L*PDstandardNth1At33 + beta2L*PDstandardNth2At33 + 
-        beta3L*PDstandardNth3At33 + 2.*
-         (At13L*PDstandardNth3beta1 + At23L*PDstandardNth3beta2) + 
-        At33L*(-0.6666666666666666666666666666666666666667*
-            (PDstandardNth1beta1 + PDstandardNth2beta2) + 
-           1.333333333333333333333333333333333333333*PDstandardNth3beta3 + 
-           alphaL*trKL) + em4phi*(Ats33 - 
-           0.3333333333333333333333333333333333333333*g33*trAts + 
-           alphaL*(-25.13274122871834590770114706623602307358*T33 + 
-              8.377580409572781969233715688745341024526*g33*trS));
+        At22L*PDstandardNth3beta2 + em4phi*(Ats23 - g23*kthird*trAts) + 
+        At23L*(-(ktwothird*PDstandardNth1beta1) + 
+           kthird*(PDstandardNth2beta2 + PDstandardNth3beta3) + alphaL*trKL);
+    
+    At33rhsL  =  beta1L*PDstandardNth1At33 + beta2L*PDstandardNth2At33 + 
+        beta3L*PDstandardNth3At33 - At33L*ktwothird*
+         (PDstandardNth1beta1 + PDstandardNth2beta2 + PDstandardNth3beta3) + 
+        2*(At13L*PDstandardNth3beta1 + At23L*PDstandardNth3beta2 + 
+           At33L*PDstandardNth3beta3) + em4phi*(Ats33 - g33*kthird*trAts) + 
+        alphaL*(-2*(At13L*Atm13 + At23L*Atm23 + At33L*Atm33) + At33L*trKL);
     
     alpharhsL  =  LapseAdvectionCoeff*
          (beta1L*PDstandardNth1alpha + beta2L*PDstandardNth2alpha + 
diff --git a/ML_BSSN/src/ML_BSSN_constraints.c b/ML_BSSN/src/ML_BSSN_constraints.c
index 5aea55e..4d86763 100644
--- a/ML_BSSN/src/ML_BSSN_constraints.c
+++ b/ML_BSSN/src/ML_BSSN_constraints.c
@@ -121,31 +121,15 @@ void ML_BSSN_constraints_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK
     CCTK_REAL gtu11 = INITVALUE, gtu21 = INITVALUE, gtu22 = INITVALUE, gtu31 = INITVALUE, gtu32 = INITVALUE, gtu33 = INITVALUE;
     CCTK_REAL gu11 = INITVALUE, gu21 = INITVALUE, gu22 = INITVALUE, gu31 = INITVALUE, gu32 = INITVALUE, gu33 = INITVALUE;
     CCTK_REAL R11 = INITVALUE, R12 = INITVALUE, R13 = INITVALUE, R22 = INITVALUE, R23 = INITVALUE, R33 = INITVALUE;
-    CCTK_REAL rho = INITVALUE;
     CCTK_REAL Rphi11 = INITVALUE, Rphi12 = INITVALUE, Rphi13 = INITVALUE, Rphi22 = INITVALUE, Rphi23 = INITVALUE, Rphi33 = INITVALUE;
     CCTK_REAL Rt11 = INITVALUE, Rt12 = INITVALUE, Rt13 = INITVALUE, Rt22 = INITVALUE, Rt23 = INITVALUE, Rt33 = INITVALUE;
-    CCTK_REAL S1 = INITVALUE, S2 = INITVALUE, S3 = INITVALUE;
-    CCTK_REAL T00 = INITVALUE, T01 = INITVALUE, T02 = INITVALUE, T03 = INITVALUE, T11 = INITVALUE, T12 = INITVALUE;
-    CCTK_REAL T13 = INITVALUE, T22 = INITVALUE, T23 = INITVALUE, T33 = INITVALUE;
     CCTK_REAL trR = INITVALUE;
     
     /* Declare local copies of grid functions */
-    CCTK_REAL alphaL = INITVALUE;
     CCTK_REAL At11L = INITVALUE, At12L = INITVALUE, At13L = INITVALUE, At22L = INITVALUE, At23L = INITVALUE, At33L = INITVALUE;
-    CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE;
     CCTK_REAL cAL = INITVALUE;
     CCTK_REAL cSL = INITVALUE;
     CCTK_REAL cXt1L = INITVALUE, cXt2L = INITVALUE, cXt3L = INITVALUE;
-    CCTK_REAL eTttL = INITVALUE;
-    CCTK_REAL eTtxL = INITVALUE;
-    CCTK_REAL eTtyL = INITVALUE;
-    CCTK_REAL eTtzL = INITVALUE;
-    CCTK_REAL eTxxL = INITVALUE;
-    CCTK_REAL eTxyL = INITVALUE;
-    CCTK_REAL eTxzL = INITVALUE;
-    CCTK_REAL eTyyL = INITVALUE;
-    CCTK_REAL eTyzL = INITVALUE;
-    CCTK_REAL eTzzL = INITVALUE;
     CCTK_REAL gt11L = INITVALUE, gt12L = INITVALUE, gt13L = INITVALUE, gt22L = INITVALUE, gt23L = INITVALUE, gt33L = INITVALUE;
     CCTK_REAL HL = INITVALUE;
     CCTK_REAL M1L = INITVALUE, M2L = INITVALUE, M3L = INITVALUE;
@@ -271,26 +255,12 @@ void ML_BSSN_constraints_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK
     CCTK_REAL PDstandardNth3Xt3 = INITVALUE;
     
     /* Assign local copies of grid functions */
-    alphaL = alpha[index];
     At11L = At11[index];
     At12L = At12[index];
     At13L = At13[index];
     At22L = At22[index];
     At23L = At23[index];
     At33L = At33[index];
-    beta1L = beta1[index];
-    beta2L = beta2[index];
-    beta3L = beta3[index];
-    eTttL = eTtt[index];
-    eTtxL = eTtx[index];
-    eTtyL = eTty[index];
-    eTtzL = eTtz[index];
-    eTxxL = eTxx[index];
-    eTxyL = eTxy[index];
-    eTxzL = eTxz[index];
-    eTyyL = eTyy[index];
-    eTyzL = eTyz[index];
-    eTzzL = eTzz[index];
     gt11L = gt11[index];
     gt12L = gt12[index];
     gt13L = gt13[index];
@@ -1062,110 +1032,77 @@ void ML_BSSN_constraints_Body(cGH const * const cctkGH, CCTK_INT const dir, CCTK
     
     Atm33  =  At13L*gtu31 + At23L*gtu32 + At33L*gtu33;
     
-    T00  =  eTttL;
-    
-    T01  =  eTtxL;
-    
-    T02  =  eTtyL;
-    
-    T03  =  eTtzL;
-    
-    T11  =  eTxxL;
-    
-    T12  =  eTxyL;
-    
-    T13  =  eTxzL;
-    
-    T22  =  eTyyL;
-    
-    T23  =  eTyzL;
-    
-    T33  =  eTzzL;
-    
-    rho  =  pow(alphaL,-2)*(T00 - 2*(beta2L*T02 + beta3L*T03) + 
-          2*(beta1L*(-T01 + beta2L*T12 + beta3L*T13) + beta2L*beta3L*T23) + 
-          T11*SQR(beta1L) + T22*SQR(beta2L) + T33*SQR(beta3L));
-    
-    S1  =  (-T01 + beta1L*T11 + beta2L*T12 + beta3L*T13)*INV(alphaL);
-    
-    S2  =  (-T02 + beta1L*T12 + beta2L*T22 + beta3L*T23)*INV(alphaL);
-    
-    S3  =  (-T03 + beta1L*T13 + beta2L*T23 + beta3L*T33)*INV(alphaL);
-    
-    HL  =  -2.*(Atm12*Atm21 + Atm13*Atm31 + Atm23*Atm32) - 
-        50.26548245743669181540229413247204614715*rho + trR - 
-        1.*(SQR(Atm11) + SQR(Atm22) + SQR(Atm33)) + 
-        0.6666666666666666666666666666666666666667*SQR(trKL);
-    
-    M1L  =  -2.*((At11L*Gt111 + At12L*Gt211 + At13L*Gt311)*gtu11 + 
-           At11L*Gt123*gtu32) - 1.*((At11L*Gt122 + At22L*Gt212 + 
-              At12L*(Gt112 + Gt222) + At23L*Gt312 + At13L*Gt322)*gtu22 + 
-           (At12L*Gt113 + At23L*Gt212)*gtu32 + 
-           (At11L*Gt133 + At23L*Gt213 + At12L*Gt233 + At33L*Gt313 + At13L*Gt333)*
-            gtu33 + At13L*(Gt111*gtu31 + Gt112*gtu32 + Gt113*gtu33)) + 
-        gtu11*(PDstandardNth1At11 + 6.*At11L*PDstandardNth1phi) - 
-        0.6666666666666666666666666666666666666667*PDstandardNth1trK + 
-        gtu22*(PDstandardNth2At12 + 6.*At12L*PDstandardNth2phi) + 
-        gtu21*(-1.*(At12L*Gt111 + At22L*Gt211 + At23L*Gt311) - 
-           3.*(At11L*Gt112 + At12L*Gt212 + At13L*Gt312) + PDstandardNth1At12 + 
-           PDstandardNth2At11 + 6.*(At12L*PDstandardNth1phi + 
-              At11L*PDstandardNth2phi)) + 
-        gtu33*(PDstandardNth3At13 + 6.*At13L*PDstandardNth3phi) + 
-        gtu31*(-1.*(At23L*Gt211 + At33L*Gt311) - 
-           3.*(At11L*Gt113 + At12L*Gt213 + At13L*Gt313) + PDstandardNth1At13 + 
-           PDstandardNth3At11 + 6.*(At13L*PDstandardNth1phi + 
-              At11L*PDstandardNth3phi)) + 
-        gtu32*(-1.*(At22L*Gt213 + At33L*Gt312 + At23L*Gt313) - 
-           2.*(At12L*Gt223 + At13L*Gt323) + PDstandardNth2At13 + 
-           PDstandardNth3At12 + 6.*(At13L*PDstandardNth2phi + 
-              At12L*PDstandardNth3phi)) - 
-        25.13274122871834590770114706623602307358*S1;
-    
-    M2L  =  -2.*((At12L*Gt122 + At22L*Gt222 + At23L*Gt322)*gtu22 + 
-           At12L*Gt113*gtu31) - 3.*((At12L*Gt112 + At22L*Gt212)*gtu21 + 
-           At12L*Gt123*gtu32) - 1.*((At11L*Gt112 + At22L*Gt211 + 
-              At12L*(Gt111 + Gt212) + At23L*Gt311 + At13L*Gt312)*gtu11 + 
-           (At13L*Gt112 + At23L*Gt212)*gtu31 + At11L*(Gt122*gtu21 + Gt123*gtu31) + 
-           (At12L*Gt133 + At22L*Gt233 + At33L*Gt323 + At23L*(Gt223 + Gt333))*
-            gtu33 + At13L*(Gt122*gtu32 + Gt123*gtu33)) + 
-        gtu11*(PDstandardNth1At12 + 6.*At12L*PDstandardNth1phi) + 
-        gtu22*(PDstandardNth2At22 + 6.*At22L*PDstandardNth2phi) + 
-        gtu21*(-3.*At23L*Gt312 - 1.*(At12L*Gt222 + At13L*Gt322) + 
-           PDstandardNth1At22 + PDstandardNth2At12 + 
-           6.*(At22L*PDstandardNth1phi + At12L*PDstandardNth2phi)) - 
-        0.6666666666666666666666666666666666666667*PDstandardNth2trK + 
-        gtu33*(PDstandardNth3At23 + 6.*At23L*PDstandardNth3phi) + 
-        gtu31*(-2.*(At22L*Gt213 + At23L*Gt313) - 
-           1.*(At12L*Gt223 + At33L*Gt312 + At13L*Gt323) + PDstandardNth1At23 + 
-           PDstandardNth3At12 + 6.*(At23L*PDstandardNth1phi + 
-              At12L*PDstandardNth3phi)) + 
-        gtu32*(-1.*(At23L*Gt222 + At33L*Gt322) - 3.*(At22L*Gt223 + At23L*Gt323) + 
-           PDstandardNth2At23 + PDstandardNth3At22 + 
-           6.*(At23L*PDstandardNth2phi + At22L*PDstandardNth3phi)) - 
-        25.13274122871834590770114706623602307358*S2;
-    
-    M3L  =  -3.*((At13L*Gt113 + At23L*Gt213)*gtu31 + At13L*Gt123*gtu32) - 
-        1.*((At11L*Gt113 + At23L*Gt211 + At12L*Gt213 + At33L*Gt311 + 
-              At13L*(Gt111 + Gt313))*gtu11 + 
-           (At13L*Gt122 + At22L*Gt223 + At33L*Gt322 + At23L*(Gt222 + Gt323))*
-            gtu22 + At11L*Gt133*gtu31 + 
-           At12L*(Gt113*gtu21 + Gt123*gtu22 + Gt133*gtu32)) + 
-        gtu11*(PDstandardNth1At13 + 6.*At13L*PDstandardNth1phi) + 
-        gtu22*(PDstandardNth2At23 + 6.*At23L*PDstandardNth2phi) + 
-        gtu21*(-2.*(At13L*Gt112 + At23L*Gt212 + At33L*Gt312) - 
-           1.*(At11L*Gt123 + At22L*Gt213 + At12L*Gt223 + At23L*Gt313 + 
-              At13L*Gt323) + PDstandardNth1At23 + PDstandardNth2At13 + 
-           6.*(At23L*PDstandardNth1phi + At13L*PDstandardNth2phi)) + 
-        gtu31*(-1.*(At12L*Gt233 + At13L*Gt333) + PDstandardNth1At33 + 
-           At33L*(-3.*Gt313 + 6.*PDstandardNth1phi) + PDstandardNth3At13 + 
-           6.*At13L*PDstandardNth3phi) + 
-        gtu33*(-2.*(At13L*Gt133 + At23L*Gt233 + At33L*Gt333) + PDstandardNth3At33 + 
-           6.*At33L*PDstandardNth3phi) + 
-        gtu32*(-3.*(At23L*Gt223 + At33L*Gt323) - 1.*(At22L*Gt233 + At23L*Gt333) + 
-           PDstandardNth2At33 + PDstandardNth3At23 + 
-           6.*(At33L*PDstandardNth2phi + At23L*PDstandardNth3phi)) - 
-        0.6666666666666666666666666666666666666667*PDstandardNth3trK - 
-        25.13274122871834590770114706623602307358*S3;
+    HL  =  -2*(Atm12*Atm21 + Atm13*Atm31 + Atm23*Atm32) + trR - SQR(Atm11) - 
+        SQR(Atm22) - SQR(Atm33) + ktwothird*SQR(trKL);
+    
+    M1L  =  -(At12L*Gt112*gtu22) - At11L*Gt122*gtu22 - At22L*Gt212*gtu22 - 
+        At12L*Gt222*gtu22 - At23L*Gt312*gtu22 - At13L*Gt322*gtu22 - 
+        At13L*Gt111*gtu31 - 3*At11L*Gt113*gtu31 - At23L*Gt211*gtu31 - 
+        3*At12L*Gt213*gtu31 - At33L*Gt311*gtu31 - 3*At13L*Gt313*gtu31 - 
+        At13L*Gt112*gtu32 - At12L*Gt113*gtu32 - At23L*Gt212*gtu32 - 
+        At22L*Gt213*gtu32 - 2*At12L*Gt223*gtu32 - At33L*Gt312*gtu32 - 
+        At23L*Gt313*gtu32 - 2*At13L*Gt323*gtu32 - 
+        2*(At11L*Gt111*gtu11 + At12L*Gt211*gtu11 + At13L*Gt311*gtu11 + 
+           At11L*Gt123*gtu32) - At13L*Gt113*gtu33 - At11L*Gt133*gtu33 - 
+        At23L*Gt213*gtu33 - At12L*Gt233*gtu33 - At33L*Gt313*gtu33 - 
+        At13L*Gt333*gtu33 + gtu11*PDstandardNth1At11 + 
+        gtu21*(-(At12L*Gt111) - 3*At11L*Gt112 - At22L*Gt211 - 3*At12L*Gt212 - 
+           At23L*Gt311 - 3*At13L*Gt312 + PDstandardNth1At12) + 
+        gtu31*PDstandardNth1At13 + 6*At11L*gtu11*PDstandardNth1phi + 
+        6*At12L*gtu21*PDstandardNth1phi + 6*At13L*gtu31*PDstandardNth1phi - 
+        ktwothird*PDstandardNth1trK + gtu21*PDstandardNth2At11 + 
+        gtu22*PDstandardNth2At12 + gtu32*PDstandardNth2At13 + 
+        6*At11L*gtu21*PDstandardNth2phi + 6*At12L*gtu22*PDstandardNth2phi + 
+        6*At13L*gtu32*PDstandardNth2phi + gtu31*PDstandardNth3At11 + 
+        gtu32*PDstandardNth3At12 + gtu33*PDstandardNth3At13 + 
+        6*At11L*gtu31*PDstandardNth3phi + 6*At12L*gtu32*PDstandardNth3phi + 
+        6*At13L*gtu33*PDstandardNth3phi;
+    
+    M2L  =  -(At12L*Gt111*gtu11) - At11L*Gt112*gtu11 - At22L*Gt211*gtu11 - 
+        At12L*Gt212*gtu11 - At23L*Gt311*gtu11 - At13L*Gt312*gtu11 - 
+        3*At12L*Gt112*gtu21 - At11L*Gt122*gtu21 - 3*At22L*Gt212*gtu21 - 
+        At12L*Gt222*gtu21 - 3*At23L*Gt312*gtu21 - At13L*Gt322*gtu21 - 
+        2*At12L*Gt122*gtu22 - 2*At22L*Gt222*gtu22 - 2*At23L*Gt322*gtu22 - 
+        At13L*Gt112*gtu31 - 2*At12L*Gt113*gtu31 - At11L*Gt123*gtu31 - 
+        At23L*Gt212*gtu31 - 2*At22L*Gt213*gtu31 - At12L*Gt223*gtu31 - 
+        At33L*Gt312*gtu31 - 2*At23L*Gt313*gtu31 - At13L*Gt323*gtu31 - 
+        At13L*Gt122*gtu32 - 3*At12L*Gt123*gtu32 - At23L*Gt222*gtu32 - 
+        3*At22L*Gt223*gtu32 - At33L*Gt322*gtu32 - 3*At23L*Gt323*gtu32 - 
+        At13L*Gt123*gtu33 - At12L*Gt133*gtu33 - At23L*Gt223*gtu33 - 
+        At22L*Gt233*gtu33 - At33L*Gt323*gtu33 - At23L*Gt333*gtu33 + 
+        gtu11*PDstandardNth1At12 + gtu21*PDstandardNth1At22 + 
+        gtu31*PDstandardNth1At23 + 6*At12L*gtu11*PDstandardNth1phi + 
+        6*At22L*gtu21*PDstandardNth1phi + 6*At23L*gtu31*PDstandardNth1phi + 
+        gtu21*PDstandardNth2At12 + gtu22*PDstandardNth2At22 + 
+        gtu32*PDstandardNth2At23 + 6*At12L*gtu21*PDstandardNth2phi + 
+        6*At22L*gtu22*PDstandardNth2phi + 6*At23L*gtu32*PDstandardNth2phi - 
+        ktwothird*PDstandardNth2trK + gtu31*PDstandardNth3At12 + 
+        gtu32*PDstandardNth3At22 + gtu33*PDstandardNth3At23 + 
+        6*At12L*gtu31*PDstandardNth3phi + 6*At22L*gtu32*PDstandardNth3phi + 
+        6*At23L*gtu33*PDstandardNth3phi;
+    
+    M3L  =  -(At13L*Gt111*gtu11) - At11L*Gt113*gtu11 - At23L*Gt211*gtu11 - 
+        At12L*Gt213*gtu11 - At33L*Gt311*gtu11 - At13L*Gt313*gtu11 - 
+        2*At13L*Gt112*gtu21 - At12L*Gt113*gtu21 - At11L*Gt123*gtu21 - 
+        2*At23L*Gt212*gtu21 - At22L*Gt213*gtu21 - At12L*Gt223*gtu21 - 
+        2*At33L*Gt312*gtu21 - At23L*Gt313*gtu21 - At13L*Gt323*gtu21 - 
+        At13L*Gt122*gtu22 - At12L*Gt123*gtu22 - At23L*Gt222*gtu22 - 
+        At22L*Gt223*gtu22 - At33L*Gt322*gtu22 - At23L*Gt323*gtu22 - 
+        3*At13L*Gt113*gtu31 - At11L*Gt133*gtu31 - 3*At23L*Gt213*gtu31 - 
+        At12L*Gt233*gtu31 - 3*At33L*Gt313*gtu31 - At13L*Gt333*gtu31 - 
+        3*At13L*Gt123*gtu32 - At12L*Gt133*gtu32 - 3*At23L*Gt223*gtu32 - 
+        At22L*Gt233*gtu32 - 3*At33L*Gt323*gtu32 - At23L*Gt333*gtu32 - 
+        2*At13L*Gt133*gtu33 - 2*At23L*Gt233*gtu33 - 2*At33L*Gt333*gtu33 + 
+        gtu11*PDstandardNth1At13 + gtu21*PDstandardNth1At23 + 
+        gtu31*PDstandardNth1At33 + 6*At13L*gtu11*PDstandardNth1phi + 
+        6*At23L*gtu21*PDstandardNth1phi + 6*At33L*gtu31*PDstandardNth1phi + 
+        gtu21*PDstandardNth2At13 + gtu22*PDstandardNth2At23 + 
+        gtu32*PDstandardNth2At33 + 6*At13L*gtu21*PDstandardNth2phi + 
+        6*At23L*gtu22*PDstandardNth2phi + 6*At33L*gtu32*PDstandardNth2phi + 
+        gtu31*PDstandardNth3At13 + gtu32*PDstandardNth3At23 + 
+        gtu33*PDstandardNth3At33 + 6*At13L*gtu31*PDstandardNth3phi + 
+        6*At23L*gtu32*PDstandardNth3phi + 6*At33L*gtu33*PDstandardNth3phi - 
+        ktwothird*PDstandardNth3trK;
     
     cSL  =  Log(detgt);