Regenerate code

5 files changed, 380 insertions, 38 deletions

diff --git a/ML_BSSN/param.ccl b/ML_BSSN/param.ccl
index e73534d..acf1f14 100644
--- a/ML_BSSN/param.ccl
+++ b/ML_BSSN/param.ccl
@@ -149,6 +149,12 @@ CCTK_INT fdOrder "fdOrder"
   *:* :: ""
 } 4
 
+restricted:
+CCTK_INT harmonicShift "Whether to use the harmonic shift"
+{
+  *:* :: ""
+} 0
+
 private:
 KEYWORD my_initial_data "my_initial_data"
 {
diff --git a/ML_BSSN/src/ML_BSSN_RHS1.cc b/ML_BSSN/src/ML_BSSN_RHS1.cc
index 9e64eab..5b5f80a 100644
--- a/ML_BSSN/src/ML_BSSN_RHS1.cc
+++ b/ML_BSSN/src/ML_BSSN_RHS1.cc
@@ -1338,17 +1338,35 @@ static void ML_BSSN_RHS1_Body(cGH const * restrict const cctkGH, int const dir,
     CCTK_REAL_VEC theta = 
       kfmin(ToReal(1),kexp(knmsub(rL,INV(ToReal(SpatialShiftGammaCoeffRadius)),ToReal(1))));
     
-    CCTK_REAL_VEC beta1rhsL = 
-      kmul(theta,kmul(kadd(Xt1L,kmadd(beta1L,kmul(eta,ToReal(BetaDriver*(-1 + 
-      ShiftBCoeff))),kmul(ksub(B1L,Xt1L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)));
+    CCTK_REAL_VEC beta1rhsL;
+    CCTK_REAL_VEC beta2rhsL;
+    CCTK_REAL_VEC beta3rhsL;
     
-    CCTK_REAL_VEC beta2rhsL = 
-      kmul(theta,kmul(kadd(Xt2L,kmadd(beta2L,kmul(eta,ToReal(BetaDriver*(-1 + 
-      ShiftBCoeff))),kmul(ksub(B2L,Xt2L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)));
-    
-    CCTK_REAL_VEC beta3rhsL = 
-      kmul(theta,kmul(kadd(Xt3L,kmadd(beta3L,kmul(eta,ToReal(BetaDriver*(-1 + 
-      ShiftBCoeff))),kmul(ksub(B3L,Xt3L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)));
+    if (harmonicShift)
+    {
+      beta1rhsL = 
+        kmul(alphaL,kmul(phiL,kmul(ToReal(0.5),kmadd(kmadd(gtu11,JacPDstandardNth1alpha,kmadd(gtu12,JacPDstandardNth2alpha,kmul(gtu13,JacPDstandardNth3alpha))),kmul(phiL,ToReal(-2)),kmul(alphaL,kmadd(phiL,kmadd(JacPDstandardNth1gt11,SQR(gtu11),kmul(JacPDstandardNth1gt22,kmul(SQR(gtu12),ToReal(2)))),kmadd(gtu13,kmadd(JacPDstandardNth3phi,ToReal(2),kmul(phiL,kmadd(gtu33,JacPDstandardNth3gt33,kmsub(kmadd(gtu13,JacPDstandardNth1gt33,kmadd(gtu22,JacPDstandardNth2gt23,kmul(gtu23,JacPDstandardNth2gt33))),ToReal(2),kmul(gtu22,JacPDstandardNth3gt22))))),kmadd(gtu11,kmadd(JacPDstandardNth1phi,ToReal(2),kmul(phiL,kmadd(gtu12,JacPDstandardNth2gt11,kmadd(gtu13,JacPDstandardNth3gt11,kmadd(gtu23,kmul(JacPDstandardNth1gt23,ToReal(-2)),knmsub(gtu22,JacPDstandardNth1gt22,kmadd(kmadd(gtu12,JacPDstandardNth1gt12,kmadd(gtu13,JacPDstandardNth1gt13,kmul(gtu22,JacPDstandardNth2gt12))),ToReal(2),kmadd(gtu23,kmul(JacPDstandardNth2gt13,ToReal(2)),kmadd(gtu23,kmul(JacPDstandardNth3gt12,ToReal(2)),kmul(gtu33,kmsub(JacPDstandardNth3gt13,ToReal(2),JacPDstandardNth1gt33))))))))))),kmul(gtu12,kmadd(JacPDstandardNth2phi,ToReal(2),kmul(phiL,kmadd(gtu22,JacPDstandardNth2gt22,kmadd(gtu23,kmul(JacPDstandardNth3gt22,ToReal(2)),kmadd(gtu33,kmsub(JacPDstandardNth3gt23,ToReal(2),JacPDstandardNth2gt33),kmul(gtu13,kmul(JacPDstandardNth1gt23,ToReal(4)))))))))))))))));
+      
+      beta2rhsL = 
+        kmul(alphaL,kmul(phiL,kmul(ToReal(0.5),kmadd(kmadd(gtu12,JacPDstandardNth1alpha,kmadd(gtu22,JacPDstandardNth2alpha,kmul(gtu23,JacPDstandardNth3alpha))),kmul(phiL,ToReal(-2)),kmul(alphaL,kmadd(phiL,kmadd(JacPDstandardNth2gt22,SQR(gtu22),kmul(JacPDstandardNth2gt11,kmul(SQR(gtu12),ToReal(2)))),kmadd(gtu23,kmadd(JacPDstandardNth3phi,ToReal(2),kmul(phiL,kmadd(gtu33,JacPDstandardNth3gt33,kmsub(kmadd(gtu11,JacPDstandardNth1gt13,kmadd(gtu13,JacPDstandardNth1gt33,kmul(gtu23,JacPDstandardNth2gt33))),ToReal(2),kmul(gtu11,JacPDstandardNth3gt11))))),kmadd(gtu22,kmadd(JacPDstandardNth2phi,ToReal(2),kmul(phiL,kmadd(gtu23,JacPDstandardNth3gt22,kmadd(kmadd(gtu23,JacPDstandardNth2gt23,kmul(gtu13,kadd(JacPDstandardNth1gt23,ksub(JacPDstandardNth3gt12,JacPDstandardNth2gt13)))),ToReal(2),kmadd(gtu11,kmsub(JacPDstandardNth1gt12,ToReal(2),JacPDstandardNth2gt11),kmul(gtu33,kmsub(JacPDstandardNth3gt23,ToReal(2),JacPDstandardNth2gt33))))))),kmul(gtu12,kmadd(JacPDstandardNth1phi,ToReal(2),kmul(phiL,kmadd(gtu11,JacPDstandardNth1gt11,kmadd(gtu13,kmul(JacPDstandardNth3gt11,ToReal(2)),kmadd(gtu22,kmadd(JacPDstandardNth2gt12,ToReal(2),JacPDstandardNth1gt22),kmadd(gtu33,kmsub(JacPDstandardNth3gt13,ToReal(2),JacPDstandardNth1gt33),kmul(gtu23,kmul(JacPDstandardNth2gt13,ToReal(4))))))))))))))))));
+      
+      beta3rhsL = 
+        kmul(alphaL,kmul(phiL,kmul(ToReal(0.5),kmadd(kmadd(gtu13,JacPDstandardNth1alpha,kmadd(gtu23,JacPDstandardNth2alpha,kmul(gtu33,JacPDstandardNth3alpha))),kmul(phiL,ToReal(-2)),kmul(alphaL,kmadd(phiL,kmul(kmadd(JacPDstandardNth3gt11,SQR(gtu13),kmul(JacPDstandardNth3gt22,SQR(gtu23))),ToReal(2)),kmadd(gtu23,kmadd(JacPDstandardNth2phi,ToReal(2),kmul(phiL,kmadd(gtu22,JacPDstandardNth2gt22,kmadd(gtu33,JacPDstandardNth2gt33,kmsub(kmadd(gtu11,JacPDstandardNth1gt12,kmadd(gtu12,JacPDstandardNth1gt22,kmul(gtu33,JacPDstandardNth3gt23))),ToReal(2),kmul(gtu11,JacPDstandardNth2gt11)))))),kmadd(gtu33,kmadd(JacPDstandardNth3phi,ToReal(2),kmul(phiL,kmadd(gtu33,JacPDstandardNth3gt33,knmsub(gtu22,JacPDstandardNth3gt22,kmadd(kmadd(gtu22,JacPDstandardNth2gt23,kmul(gtu12,kadd(JacPDstandardNth1gt23,ksub(JacPDstandardNth2gt13,JacPDstandardNth3gt12)))),ToReal(2),kmul(gtu11,kmsub(JacPDstandardNth1gt13,ToReal(2),JacPDstandardNth3gt11))))))),kmul(gtu13,kmadd(JacPDstandardNth1phi,ToReal(2),kmul(phiL,kmadd(gtu11,JacPDstandardNth1gt11,kmadd(gtu12,kmul(JacPDstandardNth2gt11,ToReal(2)),kmadd(gtu22,kmsub(JacPDstandardNth2gt12,ToReal(2),JacPDstandardNth1gt22),kmadd(gtu33,kmadd(JacPDstandardNth3gt13,ToReal(2),JacPDstandardNth1gt33),kmul(gtu23,kmul(JacPDstandardNth3gt12,ToReal(4))))))))))))))))));
+    }
+    else
+    {
+      beta1rhsL = 
+        kmul(theta,kmul(kadd(Xt1L,kmadd(beta1L,kmul(eta,ToReal(BetaDriver*(-1 + 
+        ShiftBCoeff))),kmul(ksub(B1L,Xt1L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)));
+      
+      beta2rhsL = 
+        kmul(theta,kmul(kadd(Xt2L,kmadd(beta2L,kmul(eta,ToReal(BetaDriver*(-1 + 
+        ShiftBCoeff))),kmul(ksub(B2L,Xt2L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)));
+      
+      beta3rhsL = 
+        kmul(theta,kmul(kadd(Xt3L,kmadd(beta3L,kmul(eta,ToReal(BetaDriver*(-1 + 
+        ShiftBCoeff))),kmul(ksub(B3L,Xt3L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)));
+    }
     
     CCTK_REAL_VEC B1rhsL = 
       kmul(knmsub(B1L,kmul(eta,ToReal(BetaDriver)),dotXt1),ToReal(ShiftBCoeff));
diff --git a/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShift.cc b/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShift.cc
index 6a76e6c..2abc29e 100644
--- a/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShift.cc
+++ b/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShift.cc
@@ -223,6 +223,13 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const
     CCTK_REAL_VEC beta1L = vec_load(beta1[index]);
     CCTK_REAL_VEC beta2L = vec_load(beta2[index]);
     CCTK_REAL_VEC beta3L = vec_load(beta3[index]);
+    CCTK_REAL_VEC gt11L = vec_load(gt11[index]);
+    CCTK_REAL_VEC gt12L = vec_load(gt12[index]);
+    CCTK_REAL_VEC gt13L = vec_load(gt13[index]);
+    CCTK_REAL_VEC gt22L = vec_load(gt22[index]);
+    CCTK_REAL_VEC gt23L = vec_load(gt23[index]);
+    CCTK_REAL_VEC gt33L = vec_load(gt33[index]);
+    CCTK_REAL_VEC phiL = vec_load(phi[index]);
     CCTK_REAL_VEC rL = vec_load(r[index]);
     CCTK_REAL_VEC trKL = vec_load(trK[index]);
     CCTK_REAL_VEC Xt1L = vec_load(Xt1[index]);
@@ -248,6 +255,9 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const
     /* Include user supplied include files */
     
     /* Precompute derivatives */
+    CCTK_REAL_VEC PDstandardNth1alpha;
+    CCTK_REAL_VEC PDstandardNth2alpha;
+    CCTK_REAL_VEC PDstandardNth3alpha;
     CCTK_REAL_VEC PDupwindNthAnti1alpha;
     CCTK_REAL_VEC PDupwindNthSymm1alpha;
     CCTK_REAL_VEC PDupwindNthAnti2alpha;
@@ -272,10 +282,34 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const
     CCTK_REAL_VEC PDupwindNthSymm2beta3;
     CCTK_REAL_VEC PDupwindNthAnti3beta3;
     CCTK_REAL_VEC PDupwindNthSymm3beta3;
+    CCTK_REAL_VEC PDstandardNth1gt11;
+    CCTK_REAL_VEC PDstandardNth2gt11;
+    CCTK_REAL_VEC PDstandardNth3gt11;
+    CCTK_REAL_VEC PDstandardNth1gt12;
+    CCTK_REAL_VEC PDstandardNth2gt12;
+    CCTK_REAL_VEC PDstandardNth3gt12;
+    CCTK_REAL_VEC PDstandardNth1gt13;
+    CCTK_REAL_VEC PDstandardNth2gt13;
+    CCTK_REAL_VEC PDstandardNth3gt13;
+    CCTK_REAL_VEC PDstandardNth1gt22;
+    CCTK_REAL_VEC PDstandardNth2gt22;
+    CCTK_REAL_VEC PDstandardNth3gt22;
+    CCTK_REAL_VEC PDstandardNth1gt23;
+    CCTK_REAL_VEC PDstandardNth2gt23;
+    CCTK_REAL_VEC PDstandardNth3gt23;
+    CCTK_REAL_VEC PDstandardNth1gt33;
+    CCTK_REAL_VEC PDstandardNth2gt33;
+    CCTK_REAL_VEC PDstandardNth3gt33;
+    CCTK_REAL_VEC PDstandardNth1phi;
+    CCTK_REAL_VEC PDstandardNth2phi;
+    CCTK_REAL_VEC PDstandardNth3phi;
     
     switch(fdOrder)
     {
       case 2:
+        PDstandardNth1alpha = PDstandardNthfdOrder21(&alpha[index]);
+        PDstandardNth2alpha = PDstandardNthfdOrder22(&alpha[index]);
+        PDstandardNth3alpha = PDstandardNthfdOrder23(&alpha[index]);
         PDupwindNthAnti1alpha = PDupwindNthAntifdOrder21(&alpha[index]);
         PDupwindNthSymm1alpha = PDupwindNthSymmfdOrder21(&alpha[index]);
         PDupwindNthAnti2alpha = PDupwindNthAntifdOrder22(&alpha[index]);
@@ -300,9 +334,33 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const
         PDupwindNthSymm2beta3 = PDupwindNthSymmfdOrder22(&beta3[index]);
         PDupwindNthAnti3beta3 = PDupwindNthAntifdOrder23(&beta3[index]);
         PDupwindNthSymm3beta3 = PDupwindNthSymmfdOrder23(&beta3[index]);
+        PDstandardNth1gt11 = PDstandardNthfdOrder21(&gt11[index]);
+        PDstandardNth2gt11 = PDstandardNthfdOrder22(&gt11[index]);
+        PDstandardNth3gt11 = PDstandardNthfdOrder23(&gt11[index]);
+        PDstandardNth1gt12 = PDstandardNthfdOrder21(&gt12[index]);
+        PDstandardNth2gt12 = PDstandardNthfdOrder22(&gt12[index]);
+        PDstandardNth3gt12 = PDstandardNthfdOrder23(&gt12[index]);
+        PDstandardNth1gt13 = PDstandardNthfdOrder21(&gt13[index]);
+        PDstandardNth2gt13 = PDstandardNthfdOrder22(&gt13[index]);
+        PDstandardNth3gt13 = PDstandardNthfdOrder23(&gt13[index]);
+        PDstandardNth1gt22 = PDstandardNthfdOrder21(&gt22[index]);
+        PDstandardNth2gt22 = PDstandardNthfdOrder22(&gt22[index]);
+        PDstandardNth3gt22 = PDstandardNthfdOrder23(&gt22[index]);
+        PDstandardNth1gt23 = PDstandardNthfdOrder21(&gt23[index]);
+        PDstandardNth2gt23 = PDstandardNthfdOrder22(&gt23[index]);
+        PDstandardNth3gt23 = PDstandardNthfdOrder23(&gt23[index]);
+        PDstandardNth1gt33 = PDstandardNthfdOrder21(&gt33[index]);
+        PDstandardNth2gt33 = PDstandardNthfdOrder22(&gt33[index]);
+        PDstandardNth3gt33 = PDstandardNthfdOrder23(&gt33[index]);
+        PDstandardNth1phi = PDstandardNthfdOrder21(&phi[index]);
+        PDstandardNth2phi = PDstandardNthfdOrder22(&phi[index]);
+        PDstandardNth3phi = PDstandardNthfdOrder23(&phi[index]);
         break;
       
       case 4:
+        PDstandardNth1alpha = PDstandardNthfdOrder41(&alpha[index]);
+        PDstandardNth2alpha = PDstandardNthfdOrder42(&alpha[index]);
+        PDstandardNth3alpha = PDstandardNthfdOrder43(&alpha[index]);
         PDupwindNthAnti1alpha = PDupwindNthAntifdOrder41(&alpha[index]);
         PDupwindNthSymm1alpha = PDupwindNthSymmfdOrder41(&alpha[index]);
         PDupwindNthAnti2alpha = PDupwindNthAntifdOrder42(&alpha[index]);
@@ -327,9 +385,33 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const
         PDupwindNthSymm2beta3 = PDupwindNthSymmfdOrder42(&beta3[index]);
         PDupwindNthAnti3beta3 = PDupwindNthAntifdOrder43(&beta3[index]);
         PDupwindNthSymm3beta3 = PDupwindNthSymmfdOrder43(&beta3[index]);
+        PDstandardNth1gt11 = PDstandardNthfdOrder41(&gt11[index]);
+        PDstandardNth2gt11 = PDstandardNthfdOrder42(&gt11[index]);
+        PDstandardNth3gt11 = PDstandardNthfdOrder43(&gt11[index]);
+        PDstandardNth1gt12 = PDstandardNthfdOrder41(&gt12[index]);
+        PDstandardNth2gt12 = PDstandardNthfdOrder42(&gt12[index]);
+        PDstandardNth3gt12 = PDstandardNthfdOrder43(&gt12[index]);
+        PDstandardNth1gt13 = PDstandardNthfdOrder41(&gt13[index]);
+        PDstandardNth2gt13 = PDstandardNthfdOrder42(&gt13[index]);
+        PDstandardNth3gt13 = PDstandardNthfdOrder43(&gt13[index]);
+        PDstandardNth1gt22 = PDstandardNthfdOrder41(&gt22[index]);
+        PDstandardNth2gt22 = PDstandardNthfdOrder42(&gt22[index]);
+        PDstandardNth3gt22 = PDstandardNthfdOrder43(&gt22[index]);
+        PDstandardNth1gt23 = PDstandardNthfdOrder41(&gt23[index]);
+        PDstandardNth2gt23 = PDstandardNthfdOrder42(&gt23[index]);
+        PDstandardNth3gt23 = PDstandardNthfdOrder43(&gt23[index]);
+        PDstandardNth1gt33 = PDstandardNthfdOrder41(&gt33[index]);
+        PDstandardNth2gt33 = PDstandardNthfdOrder42(&gt33[index]);
+        PDstandardNth3gt33 = PDstandardNthfdOrder43(&gt33[index]);
+        PDstandardNth1phi = PDstandardNthfdOrder41(&phi[index]);
+        PDstandardNth2phi = PDstandardNthfdOrder42(&phi[index]);
+        PDstandardNth3phi = PDstandardNthfdOrder43(&phi[index]);
         break;
       
       case 6:
+        PDstandardNth1alpha = PDstandardNthfdOrder61(&alpha[index]);
+        PDstandardNth2alpha = PDstandardNthfdOrder62(&alpha[index]);
+        PDstandardNth3alpha = PDstandardNthfdOrder63(&alpha[index]);
         PDupwindNthAnti1alpha = PDupwindNthAntifdOrder61(&alpha[index]);
         PDupwindNthSymm1alpha = PDupwindNthSymmfdOrder61(&alpha[index]);
         PDupwindNthAnti2alpha = PDupwindNthAntifdOrder62(&alpha[index]);
@@ -354,9 +436,33 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const
         PDupwindNthSymm2beta3 = PDupwindNthSymmfdOrder62(&beta3[index]);
         PDupwindNthAnti3beta3 = PDupwindNthAntifdOrder63(&beta3[index]);
         PDupwindNthSymm3beta3 = PDupwindNthSymmfdOrder63(&beta3[index]);
+        PDstandardNth1gt11 = PDstandardNthfdOrder61(&gt11[index]);
+        PDstandardNth2gt11 = PDstandardNthfdOrder62(&gt11[index]);
+        PDstandardNth3gt11 = PDstandardNthfdOrder63(&gt11[index]);
+        PDstandardNth1gt12 = PDstandardNthfdOrder61(&gt12[index]);
+        PDstandardNth2gt12 = PDstandardNthfdOrder62(&gt12[index]);
+        PDstandardNth3gt12 = PDstandardNthfdOrder63(&gt12[index]);
+        PDstandardNth1gt13 = PDstandardNthfdOrder61(&gt13[index]);
+        PDstandardNth2gt13 = PDstandardNthfdOrder62(&gt13[index]);
+        PDstandardNth3gt13 = PDstandardNthfdOrder63(&gt13[index]);
+        PDstandardNth1gt22 = PDstandardNthfdOrder61(&gt22[index]);
+        PDstandardNth2gt22 = PDstandardNthfdOrder62(&gt22[index]);
+        PDstandardNth3gt22 = PDstandardNthfdOrder63(&gt22[index]);
+        PDstandardNth1gt23 = PDstandardNthfdOrder61(&gt23[index]);
+        PDstandardNth2gt23 = PDstandardNthfdOrder62(&gt23[index]);
+        PDstandardNth3gt23 = PDstandardNthfdOrder63(&gt23[index]);
+        PDstandardNth1gt33 = PDstandardNthfdOrder61(&gt33[index]);
+        PDstandardNth2gt33 = PDstandardNthfdOrder62(&gt33[index]);
+        PDstandardNth3gt33 = PDstandardNthfdOrder63(&gt33[index]);
+        PDstandardNth1phi = PDstandardNthfdOrder61(&phi[index]);
+        PDstandardNth2phi = PDstandardNthfdOrder62(&phi[index]);
+        PDstandardNth3phi = PDstandardNthfdOrder63(&phi[index]);
         break;
       
       case 8:
+        PDstandardNth1alpha = PDstandardNthfdOrder81(&alpha[index]);
+        PDstandardNth2alpha = PDstandardNthfdOrder82(&alpha[index]);
+        PDstandardNth3alpha = PDstandardNthfdOrder83(&alpha[index]);
         PDupwindNthAnti1alpha = PDupwindNthAntifdOrder81(&alpha[index]);
         PDupwindNthSymm1alpha = PDupwindNthSymmfdOrder81(&alpha[index]);
         PDupwindNthAnti2alpha = PDupwindNthAntifdOrder82(&alpha[index]);
@@ -381,6 +487,27 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const
         PDupwindNthSymm2beta3 = PDupwindNthSymmfdOrder82(&beta3[index]);
         PDupwindNthAnti3beta3 = PDupwindNthAntifdOrder83(&beta3[index]);
         PDupwindNthSymm3beta3 = PDupwindNthSymmfdOrder83(&beta3[index]);
+        PDstandardNth1gt11 = PDstandardNthfdOrder81(&gt11[index]);
+        PDstandardNth2gt11 = PDstandardNthfdOrder82(&gt11[index]);
+        PDstandardNth3gt11 = PDstandardNthfdOrder83(&gt11[index]);
+        PDstandardNth1gt12 = PDstandardNthfdOrder81(&gt12[index]);
+        PDstandardNth2gt12 = PDstandardNthfdOrder82(&gt12[index]);
+        PDstandardNth3gt12 = PDstandardNthfdOrder83(&gt12[index]);
+        PDstandardNth1gt13 = PDstandardNthfdOrder81(&gt13[index]);
+        PDstandardNth2gt13 = PDstandardNthfdOrder82(&gt13[index]);
+        PDstandardNth3gt13 = PDstandardNthfdOrder83(&gt13[index]);
+        PDstandardNth1gt22 = PDstandardNthfdOrder81(&gt22[index]);
+        PDstandardNth2gt22 = PDstandardNthfdOrder82(&gt22[index]);
+        PDstandardNth3gt22 = PDstandardNthfdOrder83(&gt22[index]);
+        PDstandardNth1gt23 = PDstandardNthfdOrder81(&gt23[index]);
+        PDstandardNth2gt23 = PDstandardNthfdOrder82(&gt23[index]);
+        PDstandardNth3gt23 = PDstandardNthfdOrder83(&gt23[index]);
+        PDstandardNth1gt33 = PDstandardNthfdOrder81(&gt33[index]);
+        PDstandardNth2gt33 = PDstandardNthfdOrder82(&gt33[index]);
+        PDstandardNth3gt33 = PDstandardNthfdOrder83(&gt33[index]);
+        PDstandardNth1phi = PDstandardNthfdOrder81(&phi[index]);
+        PDstandardNth2phi = PDstandardNthfdOrder82(&phi[index]);
+        PDstandardNth3phi = PDstandardNthfdOrder83(&phi[index]);
         break;
     }
     
@@ -391,6 +518,30 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const
     
     ptrdiff_t dir3 = Sign(beta3L);
     
+    CCTK_REAL_VEC JacPDstandardNth1alpha;
+    CCTK_REAL_VEC JacPDstandardNth1gt11;
+    CCTK_REAL_VEC JacPDstandardNth1gt12;
+    CCTK_REAL_VEC JacPDstandardNth1gt13;
+    CCTK_REAL_VEC JacPDstandardNth1gt22;
+    CCTK_REAL_VEC JacPDstandardNth1gt23;
+    CCTK_REAL_VEC JacPDstandardNth1gt33;
+    CCTK_REAL_VEC JacPDstandardNth1phi;
+    CCTK_REAL_VEC JacPDstandardNth2alpha;
+    CCTK_REAL_VEC JacPDstandardNth2gt11;
+    CCTK_REAL_VEC JacPDstandardNth2gt12;
+    CCTK_REAL_VEC JacPDstandardNth2gt13;
+    CCTK_REAL_VEC JacPDstandardNth2gt22;
+    CCTK_REAL_VEC JacPDstandardNth2gt23;
+    CCTK_REAL_VEC JacPDstandardNth2gt33;
+    CCTK_REAL_VEC JacPDstandardNth2phi;
+    CCTK_REAL_VEC JacPDstandardNth3alpha;
+    CCTK_REAL_VEC JacPDstandardNth3gt11;
+    CCTK_REAL_VEC JacPDstandardNth3gt12;
+    CCTK_REAL_VEC JacPDstandardNth3gt13;
+    CCTK_REAL_VEC JacPDstandardNth3gt22;
+    CCTK_REAL_VEC JacPDstandardNth3gt23;
+    CCTK_REAL_VEC JacPDstandardNth3gt33;
+    CCTK_REAL_VEC JacPDstandardNth3phi;
     CCTK_REAL_VEC JacPDupwindNthAnti1alpha;
     CCTK_REAL_VEC JacPDupwindNthAnti1beta1;
     CCTK_REAL_VEC JacPDupwindNthAnti1beta2;
@@ -418,6 +569,78 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const
     
     if (use_jacobian)
     {
+      JacPDstandardNth1alpha = 
+        kmadd(J11L,PDstandardNth1alpha,kmadd(J21L,PDstandardNth2alpha,kmul(J31L,PDstandardNth3alpha)));
+      
+      JacPDstandardNth1gt11 = 
+        kmadd(J11L,PDstandardNth1gt11,kmadd(J21L,PDstandardNth2gt11,kmul(J31L,PDstandardNth3gt11)));
+      
+      JacPDstandardNth1gt12 = 
+        kmadd(J11L,PDstandardNth1gt12,kmadd(J21L,PDstandardNth2gt12,kmul(J31L,PDstandardNth3gt12)));
+      
+      JacPDstandardNth1gt13 = 
+        kmadd(J11L,PDstandardNth1gt13,kmadd(J21L,PDstandardNth2gt13,kmul(J31L,PDstandardNth3gt13)));
+      
+      JacPDstandardNth1gt22 = 
+        kmadd(J11L,PDstandardNth1gt22,kmadd(J21L,PDstandardNth2gt22,kmul(J31L,PDstandardNth3gt22)));
+      
+      JacPDstandardNth1gt23 = 
+        kmadd(J11L,PDstandardNth1gt23,kmadd(J21L,PDstandardNth2gt23,kmul(J31L,PDstandardNth3gt23)));
+      
+      JacPDstandardNth1gt33 = 
+        kmadd(J11L,PDstandardNth1gt33,kmadd(J21L,PDstandardNth2gt33,kmul(J31L,PDstandardNth3gt33)));
+      
+      JacPDstandardNth1phi = 
+        kmadd(J11L,PDstandardNth1phi,kmadd(J21L,PDstandardNth2phi,kmul(J31L,PDstandardNth3phi)));
+      
+      JacPDstandardNth2alpha = 
+        kmadd(J12L,PDstandardNth1alpha,kmadd(J22L,PDstandardNth2alpha,kmul(J32L,PDstandardNth3alpha)));
+      
+      JacPDstandardNth2gt11 = 
+        kmadd(J12L,PDstandardNth1gt11,kmadd(J22L,PDstandardNth2gt11,kmul(J32L,PDstandardNth3gt11)));
+      
+      JacPDstandardNth2gt12 = 
+        kmadd(J12L,PDstandardNth1gt12,kmadd(J22L,PDstandardNth2gt12,kmul(J32L,PDstandardNth3gt12)));
+      
+      JacPDstandardNth2gt13 = 
+        kmadd(J12L,PDstandardNth1gt13,kmadd(J22L,PDstandardNth2gt13,kmul(J32L,PDstandardNth3gt13)));
+      
+      JacPDstandardNth2gt22 = 
+        kmadd(J12L,PDstandardNth1gt22,kmadd(J22L,PDstandardNth2gt22,kmul(J32L,PDstandardNth3gt22)));
+      
+      JacPDstandardNth2gt23 = 
+        kmadd(J12L,PDstandardNth1gt23,kmadd(J22L,PDstandardNth2gt23,kmul(J32L,PDstandardNth3gt23)));
+      
+      JacPDstandardNth2gt33 = 
+        kmadd(J12L,PDstandardNth1gt33,kmadd(J22L,PDstandardNth2gt33,kmul(J32L,PDstandardNth3gt33)));
+      
+      JacPDstandardNth2phi = 
+        kmadd(J12L,PDstandardNth1phi,kmadd(J22L,PDstandardNth2phi,kmul(J32L,PDstandardNth3phi)));
+      
+      JacPDstandardNth3alpha = 
+        kmadd(J13L,PDstandardNth1alpha,kmadd(J23L,PDstandardNth2alpha,kmul(J33L,PDstandardNth3alpha)));
+      
+      JacPDstandardNth3gt11 = 
+        kmadd(J13L,PDstandardNth1gt11,kmadd(J23L,PDstandardNth2gt11,kmul(J33L,PDstandardNth3gt11)));
+      
+      JacPDstandardNth3gt12 = 
+        kmadd(J13L,PDstandardNth1gt12,kmadd(J23L,PDstandardNth2gt12,kmul(J33L,PDstandardNth3gt12)));
+      
+      JacPDstandardNth3gt13 = 
+        kmadd(J13L,PDstandardNth1gt13,kmadd(J23L,PDstandardNth2gt13,kmul(J33L,PDstandardNth3gt13)));
+      
+      JacPDstandardNth3gt22 = 
+        kmadd(J13L,PDstandardNth1gt22,kmadd(J23L,PDstandardNth2gt22,kmul(J33L,PDstandardNth3gt22)));
+      
+      JacPDstandardNth3gt23 = 
+        kmadd(J13L,PDstandardNth1gt23,kmadd(J23L,PDstandardNth2gt23,kmul(J33L,PDstandardNth3gt23)));
+      
+      JacPDstandardNth3gt33 = 
+        kmadd(J13L,PDstandardNth1gt33,kmadd(J23L,PDstandardNth2gt33,kmul(J33L,PDstandardNth3gt33)));
+      
+      JacPDstandardNth3phi = 
+        kmadd(J13L,PDstandardNth1phi,kmadd(J23L,PDstandardNth2phi,kmul(J33L,PDstandardNth3phi)));
+      
       JacPDupwindNthAnti1alpha = 
         kmadd(J11L,PDupwindNthAnti1alpha,kmadd(J21L,PDupwindNthAnti2alpha,kmul(J31L,PDupwindNthAnti3alpha)));
       
@@ -492,6 +715,54 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const
     }
     else
     {
+      JacPDstandardNth1alpha = PDstandardNth1alpha;
+      
+      JacPDstandardNth1gt11 = PDstandardNth1gt11;
+      
+      JacPDstandardNth1gt12 = PDstandardNth1gt12;
+      
+      JacPDstandardNth1gt13 = PDstandardNth1gt13;
+      
+      JacPDstandardNth1gt22 = PDstandardNth1gt22;
+      
+      JacPDstandardNth1gt23 = PDstandardNth1gt23;
+      
+      JacPDstandardNth1gt33 = PDstandardNth1gt33;
+      
+      JacPDstandardNth1phi = PDstandardNth1phi;
+      
+      JacPDstandardNth2alpha = PDstandardNth2alpha;
+      
+      JacPDstandardNth2gt11 = PDstandardNth2gt11;
+      
+      JacPDstandardNth2gt12 = PDstandardNth2gt12;
+      
+      JacPDstandardNth2gt13 = PDstandardNth2gt13;
+      
+      JacPDstandardNth2gt22 = PDstandardNth2gt22;
+      
+      JacPDstandardNth2gt23 = PDstandardNth2gt23;
+      
+      JacPDstandardNth2gt33 = PDstandardNth2gt33;
+      
+      JacPDstandardNth2phi = PDstandardNth2phi;
+      
+      JacPDstandardNth3alpha = PDstandardNth3alpha;
+      
+      JacPDstandardNth3gt11 = PDstandardNth3gt11;
+      
+      JacPDstandardNth3gt12 = PDstandardNth3gt12;
+      
+      JacPDstandardNth3gt13 = PDstandardNth3gt13;
+      
+      JacPDstandardNth3gt22 = PDstandardNth3gt22;
+      
+      JacPDstandardNth3gt23 = PDstandardNth3gt23;
+      
+      JacPDstandardNth3gt33 = PDstandardNth3gt33;
+      
+      JacPDstandardNth3phi = PDstandardNth3phi;
+      
       JacPDupwindNthAnti1alpha = PDupwindNthAnti1alpha;
       
       JacPDupwindNthAnti1beta1 = PDupwindNthAnti1beta1;
@@ -541,6 +812,23 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const
       JacPDupwindNthSymm3beta3 = PDupwindNthSymm3beta3;
     }
     
+    CCTK_REAL_VEC detgt = ToReal(1);
+    
+    CCTK_REAL_VEC gtu11 = kmul(INV(detgt),kmsub(gt22L,gt33L,SQR(gt23L)));
+    
+    CCTK_REAL_VEC gtu12 = 
+      kmul(INV(detgt),kmsub(gt13L,gt23L,kmul(gt12L,gt33L)));
+    
+    CCTK_REAL_VEC gtu13 = 
+      kmul(INV(detgt),kmsub(gt12L,gt23L,kmul(gt13L,gt22L)));
+    
+    CCTK_REAL_VEC gtu22 = kmul(INV(detgt),kmsub(gt11L,gt33L,SQR(gt13L)));
+    
+    CCTK_REAL_VEC gtu23 = 
+      kmul(INV(detgt),kmsub(gt12L,gt13L,kmul(gt11L,gt23L)));
+    
+    CCTK_REAL_VEC gtu33 = kmul(INV(detgt),kmsub(gt11L,gt22L,SQR(gt12L)));
+    
     CCTK_REAL_VEC eta = 
       kfmin(ToReal(1),kmul(INV(rL),ToReal(SpatialBetaDriverRadius)));
     
@@ -551,19 +839,19 @@ static void ML_BSSN_convertToADMBaseDtLapseShift_Body(cGH const * restrict const
       kmsub(kmadd(beta1L,JacPDupwindNthAnti1alpha,kmadd(beta2L,JacPDupwindNthAnti2alpha,kmadd(beta3L,JacPDupwindNthAnti3alpha,kmadd(JacPDupwindNthSymm1alpha,kfabs(beta1L),kmadd(JacPDupwindNthSymm2alpha,kfabs(beta2L),kmul(JacPDupwindNthSymm3alpha,kfabs(beta3L))))))),ToReal(LapseAdvectionCoeff),kmul(kpow(alphaL,harmonicN),kmul(ToReal(harmonicF),kmadd(ksub(AL,trKL),ToReal(LapseACoeff),trKL))));
     
     CCTK_REAL_VEC dtbetaxL = 
-      kmadd(kmadd(beta1L,JacPDupwindNthAnti1beta1,kmadd(beta2L,JacPDupwindNthAnti2beta1,kmadd(beta3L,JacPDupwindNthAnti3beta1,kmadd(JacPDupwindNthSymm1beta1,kfabs(beta1L),kmadd(JacPDupwindNthSymm2beta1,kfabs(beta2L),kmul(JacPDupwindNthSymm3beta1,kfabs(beta3L))))))),ToReal(ShiftAdvectionCoeff),kmul(theta,kmul(kadd(Xt1L,kmadd(beta1L,kmul(eta,ToReal(BetaDriver*(-1 
+      kmadd(kmadd(beta1L,JacPDupwindNthAnti1beta1,kmadd(beta2L,JacPDupwindNthAnti2beta1,kmadd(beta3L,JacPDupwindNthAnti3beta1,kmadd(JacPDupwindNthSymm1beta1,kfabs(beta1L),kmadd(JacPDupwindNthSymm2beta1,kfabs(beta2L),kmul(JacPDupwindNthSymm3beta1,kfabs(beta3L))))))),ToReal(ShiftAdvectionCoeff),IfThen(harmonicShift,kmul(alphaL,kmul(phiL,kmul(ToReal(0.5),kmadd(kmadd(gtu11,JacPDstandardNth1alpha,kmadd(gtu12,JacPDstandardNth2alpha,kmul(gtu13,JacPDstandardNth3alpha))),kmul(phiL,ToReal(-2)),kmul(alphaL,kmadd(phiL,kmadd(JacPDstandardNth1gt11,SQR(gtu11),kmul(JacPDstandardNth1gt22,kmul(SQR(gtu12),ToReal(2)))),kmadd(gtu13,kmadd(JacPDstandardNth3phi,ToReal(2),kmul(phiL,kmadd(gtu33,JacPDstandardNth3gt33,kmsub(kmadd(gtu13,JacPDstandardNth1gt33,kmadd(gtu22,JacPDstandardNth2gt23,kmul(gtu23,JacPDstandardNth2gt33))),ToReal(2),kmul(gtu22,JacPDstandardNth3gt22))))),kmadd(gtu11,kmadd(JacPDstandardNth1phi,ToReal(2),kmul(phiL,kmadd(gtu12,JacPDstandardNth2gt11,kmadd(gtu13,JacPDstandardNth3gt11,kmadd(gtu23,kmul(JacPDstandardNth1gt23,ToReal(-2)),knmsub(gtu22,JacPDstandardNth1gt22,kmadd(kmadd(gtu12,JacPDstandardNth1gt12,kmadd(gtu13,JacPDstandardNth1gt13,kmul(gtu22,JacPDstandardNth2gt12))),ToReal(2),kmadd(gtu23,kmul(JacPDstandardNth2gt13,ToReal(2)),kmadd(gtu23,kmul(JacPDstandardNth3gt12,ToReal(2)),kmul(gtu33,kmsub(JacPDstandardNth3gt13,ToReal(2),JacPDstandardNth1gt33))))))))))),kmul(gtu12,kmadd(JacPDstandardNth2phi,ToReal(2),kmul(phiL,kmadd(gtu22,JacPDstandardNth2gt22,kmadd(gtu23,kmul(JacPDstandardNth3gt22,ToReal(2)),kmadd(gtu33,kmsub(JacPDstandardNth3gt23,ToReal(2),JacPDstandardNth2gt33),kmul(gtu13,kmul(JacPDstandardNth1gt23,ToReal(4))))))))))))))))),kmul(theta,kmul(kadd(Xt1L,kmadd(beta1L,kmul(eta,ToReal(BetaDriver*(-1 
       + 
-      ShiftBCoeff))),kmul(ksub(B1L,Xt1L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))));
+      ShiftBCoeff))),kmul(ksub(B1L,Xt1L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)))));
     
     CCTK_REAL_VEC dtbetayL = 
-      kmadd(kmadd(beta1L,JacPDupwindNthAnti1beta2,kmadd(beta2L,JacPDupwindNthAnti2beta2,kmadd(beta3L,JacPDupwindNthAnti3beta2,kmadd(JacPDupwindNthSymm1beta2,kfabs(beta1L),kmadd(JacPDupwindNthSymm2beta2,kfabs(beta2L),kmul(JacPDupwindNthSymm3beta2,kfabs(beta3L))))))),ToReal(ShiftAdvectionCoeff),kmul(theta,kmul(kadd(Xt2L,kmadd(beta2L,kmul(eta,ToReal(BetaDriver*(-1 
+      kmadd(kmadd(beta1L,JacPDupwindNthAnti1beta2,kmadd(beta2L,JacPDupwindNthAnti2beta2,kmadd(beta3L,JacPDupwindNthAnti3beta2,kmadd(JacPDupwindNthSymm1beta2,kfabs(beta1L),kmadd(JacPDupwindNthSymm2beta2,kfabs(beta2L),kmul(JacPDupwindNthSymm3beta2,kfabs(beta3L))))))),ToReal(ShiftAdvectionCoeff),IfThen(harmonicShift,kmul(alphaL,kmul(phiL,kmul(ToReal(0.5),kmadd(kmadd(gtu12,JacPDstandardNth1alpha,kmadd(gtu22,JacPDstandardNth2alpha,kmul(gtu23,JacPDstandardNth3alpha))),kmul(phiL,ToReal(-2)),kmul(alphaL,kmadd(phiL,kmadd(JacPDstandardNth2gt22,SQR(gtu22),kmul(JacPDstandardNth2gt11,kmul(SQR(gtu12),ToReal(2)))),kmadd(gtu23,kmadd(JacPDstandardNth3phi,ToReal(2),kmul(phiL,kmadd(gtu33,JacPDstandardNth3gt33,kmsub(kmadd(gtu11,JacPDstandardNth1gt13,kmadd(gtu13,JacPDstandardNth1gt33,kmul(gtu23,JacPDstandardNth2gt33))),ToReal(2),kmul(gtu11,JacPDstandardNth3gt11))))),kmadd(gtu22,kmadd(JacPDstandardNth2phi,ToReal(2),kmul(phiL,kmadd(gtu23,JacPDstandardNth3gt22,kmadd(kmadd(gtu23,JacPDstandardNth2gt23,kmul(gtu13,kadd(JacPDstandardNth1gt23,ksub(JacPDstandardNth3gt12,JacPDstandardNth2gt13)))),ToReal(2),kmadd(gtu11,kmsub(JacPDstandardNth1gt12,ToReal(2),JacPDstandardNth2gt11),kmul(gtu33,kmsub(JacPDstandardNth3gt23,ToReal(2),JacPDstandardNth2gt33))))))),kmul(gtu12,kmadd(JacPDstandardNth1phi,ToReal(2),kmul(phiL,kmadd(gtu11,JacPDstandardNth1gt11,kmadd(gtu13,kmul(JacPDstandardNth3gt11,ToReal(2)),kmadd(gtu22,kmadd(JacPDstandardNth2gt12,ToReal(2),JacPDstandardNth1gt22),kmadd(gtu33,kmsub(JacPDstandardNth3gt13,ToReal(2),JacPDstandardNth1gt33),kmul(gtu23,kmul(JacPDstandardNth2gt13,ToReal(4)))))))))))))))))),kmul(theta,kmul(kadd(Xt2L,kmadd(beta2L,kmul(eta,ToReal(BetaDriver*(-1 
       + 
-      ShiftBCoeff))),kmul(ksub(B2L,Xt2L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))));
+      ShiftBCoeff))),kmul(ksub(B2L,Xt2L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)))));
     
     CCTK_REAL_VEC dtbetazL = 
-      kmadd(kmadd(beta1L,JacPDupwindNthAnti1beta3,kmadd(beta2L,JacPDupwindNthAnti2beta3,kmadd(beta3L,JacPDupwindNthAnti3beta3,kmadd(JacPDupwindNthSymm1beta3,kfabs(beta1L),kmadd(JacPDupwindNthSymm2beta3,kfabs(beta2L),kmul(JacPDupwindNthSymm3beta3,kfabs(beta3L))))))),ToReal(ShiftAdvectionCoeff),kmul(theta,kmul(kadd(Xt3L,kmadd(beta3L,kmul(eta,ToReal(BetaDriver*(-1 
+      kmadd(kmadd(beta1L,JacPDupwindNthAnti1beta3,kmadd(beta2L,JacPDupwindNthAnti2beta3,kmadd(beta3L,JacPDupwindNthAnti3beta3,kmadd(JacPDupwindNthSymm1beta3,kfabs(beta1L),kmadd(JacPDupwindNthSymm2beta3,kfabs(beta2L),kmul(JacPDupwindNthSymm3beta3,kfabs(beta3L))))))),ToReal(ShiftAdvectionCoeff),IfThen(harmonicShift,kmul(alphaL,kmul(phiL,kmul(ToReal(0.5),kmadd(kmadd(gtu13,JacPDstandardNth1alpha,kmadd(gtu23,JacPDstandardNth2alpha,kmul(gtu33,JacPDstandardNth3alpha))),kmul(phiL,ToReal(-2)),kmul(alphaL,kmadd(phiL,kmul(kmadd(JacPDstandardNth3gt11,SQR(gtu13),kmul(JacPDstandardNth3gt22,SQR(gtu23))),ToReal(2)),kmadd(gtu23,kmadd(JacPDstandardNth2phi,ToReal(2),kmul(phiL,kmadd(gtu22,JacPDstandardNth2gt22,kmadd(gtu33,JacPDstandardNth2gt33,kmsub(kmadd(gtu11,JacPDstandardNth1gt12,kmadd(gtu12,JacPDstandardNth1gt22,kmul(gtu33,JacPDstandardNth3gt23))),ToReal(2),kmul(gtu11,JacPDstandardNth2gt11)))))),kmadd(gtu33,kmadd(JacPDstandardNth3phi,ToReal(2),kmul(phiL,kmadd(gtu33,JacPDstandardNth3gt33,knmsub(gtu22,JacPDstandardNth3gt22,kmadd(kmadd(gtu22,JacPDstandardNth2gt23,kmul(gtu12,kadd(JacPDstandardNth1gt23,ksub(JacPDstandardNth2gt13,JacPDstandardNth3gt12)))),ToReal(2),kmul(gtu11,kmsub(JacPDstandardNth1gt13,ToReal(2),JacPDstandardNth3gt11))))))),kmul(gtu13,kmadd(JacPDstandardNth1phi,ToReal(2),kmul(phiL,kmadd(gtu11,JacPDstandardNth1gt11,kmadd(gtu12,kmul(JacPDstandardNth2gt11,ToReal(2)),kmadd(gtu22,kmsub(JacPDstandardNth2gt12,ToReal(2),JacPDstandardNth1gt22),kmadd(gtu33,kmadd(JacPDstandardNth3gt13,ToReal(2),JacPDstandardNth1gt33),kmul(gtu23,kmul(JacPDstandardNth3gt12,ToReal(4)))))))))))))))))),kmul(theta,kmul(kadd(Xt3L,kmadd(beta3L,kmul(eta,ToReal(BetaDriver*(-1 
       + 
-      ShiftBCoeff))),kmul(ksub(B3L,Xt3L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff))));
+      ShiftBCoeff))),kmul(ksub(B3L,Xt3L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)))));
     
     /* If necessary, store only partial vectors after the first iteration */
     
@@ -625,8 +913,8 @@ extern "C" void ML_BSSN_convertToADMBaseDtLapseShift(CCTK_ARGUMENTS)
     return;
   }
   
-  const char *groups[] = {"ADMBase::dtlapse","ADMBase::dtshift","grid::coordinates","Grid::coordinates","ML_BSSN::ML_dtlapse","ML_BSSN::ML_dtshift","ML_BSSN::ML_Gamma","ML_BSSN::ML_lapse","ML_BSSN::ML_shift","ML_BSSN::ML_trace_curv"};
-  GenericFD_AssertGroupStorage(cctkGH, "ML_BSSN_convertToADMBaseDtLapseShift", 10, groups);
+  const char *groups[] = {"ADMBase::dtlapse","ADMBase::dtshift","grid::coordinates","Grid::coordinates","ML_BSSN::ML_dtlapse","ML_BSSN::ML_dtshift","ML_BSSN::ML_Gamma","ML_BSSN::ML_lapse","ML_BSSN::ML_log_confac","ML_BSSN::ML_metric","ML_BSSN::ML_shift","ML_BSSN::ML_trace_curv"};
+  GenericFD_AssertGroupStorage(cctkGH, "ML_BSSN_convertToADMBaseDtLapseShift", 12, groups);
   
   switch(fdOrder)
   {
diff --git a/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShiftBoundary.cc b/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShiftBoundary.cc
index db5de18..a957b55 100644
--- a/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShiftBoundary.cc
+++ b/ML_BSSN/src/ML_BSSN_convertToADMBaseDtLapseShiftBoundary.cc
@@ -260,17 +260,32 @@ static void ML_BSSN_convertToADMBaseDtLapseShiftBoundary_Body(cGH const * restri
     CCTK_REAL_VEC dtalpL = 
       kneg(kmul(kpow(alphaL,harmonicN),kmul(ToReal(harmonicF),kmadd(ksub(AL,trKL),ToReal(LapseACoeff),trKL))));
     
-    CCTK_REAL_VEC dtbetaxL = 
-      kmul(theta,kmul(kadd(Xt1L,kmadd(beta1L,kmul(eta,ToReal(BetaDriver*(-1 + 
-      ShiftBCoeff))),kmul(ksub(B1L,Xt1L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)));
+    CCTK_REAL_VEC dtbetaxL;
+    CCTK_REAL_VEC dtbetayL;
+    CCTK_REAL_VEC dtbetazL;
     
-    CCTK_REAL_VEC dtbetayL = 
-      kmul(theta,kmul(kadd(Xt2L,kmadd(beta2L,kmul(eta,ToReal(BetaDriver*(-1 + 
-      ShiftBCoeff))),kmul(ksub(B2L,Xt2L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)));
-    
-    CCTK_REAL_VEC dtbetazL = 
-      kmul(theta,kmul(kadd(Xt3L,kmadd(beta3L,kmul(eta,ToReal(BetaDriver*(-1 + 
-      ShiftBCoeff))),kmul(ksub(B3L,Xt3L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)));
+    if (harmonicShift)
+    {
+      dtbetaxL = ToReal(0);
+      
+      dtbetayL = ToReal(0);
+      
+      dtbetazL = ToReal(0);
+    }
+    else
+    {
+      dtbetaxL = 
+        kmul(theta,kmul(kadd(Xt1L,kmadd(beta1L,kmul(eta,ToReal(BetaDriver*(-1 + 
+        ShiftBCoeff))),kmul(ksub(B1L,Xt1L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)));
+      
+      dtbetayL = 
+        kmul(theta,kmul(kadd(Xt2L,kmadd(beta2L,kmul(eta,ToReal(BetaDriver*(-1 + 
+        ShiftBCoeff))),kmul(ksub(B2L,Xt2L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)));
+      
+      dtbetazL = 
+        kmul(theta,kmul(kadd(Xt3L,kmadd(beta3L,kmul(eta,ToReal(BetaDriver*(-1 + 
+        ShiftBCoeff))),kmul(ksub(B3L,Xt3L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)));
+    }
     
     /* If necessary, store only partial vectors after the first iteration */
     
diff --git a/ML_BSSN/src/ML_BSSN_convertToADMBaseFakeDtLapseShift.cc b/ML_BSSN/src/ML_BSSN_convertToADMBaseFakeDtLapseShift.cc
index a93eb01..0d8065d 100644
--- a/ML_BSSN/src/ML_BSSN_convertToADMBaseFakeDtLapseShift.cc
+++ b/ML_BSSN/src/ML_BSSN_convertToADMBaseFakeDtLapseShift.cc
@@ -245,17 +245,32 @@ static void ML_BSSN_convertToADMBaseFakeDtLapseShift_Body(cGH const * restrict c
     CCTK_REAL_VEC dtalpL = 
       kneg(kmul(kpow(alphaL,harmonicN),kmul(ToReal(harmonicF),kmadd(ksub(AL,trKL),ToReal(LapseACoeff),trKL))));
     
-    CCTK_REAL_VEC dtbetaxL = 
-      kmul(theta,kmul(kadd(Xt1L,kmadd(beta1L,kmul(eta,ToReal(BetaDriver*(-1 + 
-      ShiftBCoeff))),kmul(ksub(B1L,Xt1L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)));
+    CCTK_REAL_VEC dtbetaxL;
+    CCTK_REAL_VEC dtbetayL;
+    CCTK_REAL_VEC dtbetazL;
     
-    CCTK_REAL_VEC dtbetayL = 
-      kmul(theta,kmul(kadd(Xt2L,kmadd(beta2L,kmul(eta,ToReal(BetaDriver*(-1 + 
-      ShiftBCoeff))),kmul(ksub(B2L,Xt2L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)));
-    
-    CCTK_REAL_VEC dtbetazL = 
-      kmul(theta,kmul(kadd(Xt3L,kmadd(beta3L,kmul(eta,ToReal(BetaDriver*(-1 + 
-      ShiftBCoeff))),kmul(ksub(B3L,Xt3L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)));
+    if (harmonicShift)
+    {
+      dtbetaxL = ToReal(0);
+      
+      dtbetayL = ToReal(0);
+      
+      dtbetazL = ToReal(0);
+    }
+    else
+    {
+      dtbetaxL = 
+        kmul(theta,kmul(kadd(Xt1L,kmadd(beta1L,kmul(eta,ToReal(BetaDriver*(-1 + 
+        ShiftBCoeff))),kmul(ksub(B1L,Xt1L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)));
+      
+      dtbetayL = 
+        kmul(theta,kmul(kadd(Xt2L,kmadd(beta2L,kmul(eta,ToReal(BetaDriver*(-1 + 
+        ShiftBCoeff))),kmul(ksub(B2L,Xt2L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)));
+      
+      dtbetazL = 
+        kmul(theta,kmul(kadd(Xt3L,kmadd(beta3L,kmul(eta,ToReal(BetaDriver*(-1 + 
+        ShiftBCoeff))),kmul(ksub(B3L,Xt3L),ToReal(ShiftBCoeff)))),ToReal(ShiftGammaCoeff)));
+    }
     
     /* If necessary, store only partial vectors after the first iteration */