Merge branch 'master' into CCZ4

Conflicts:
	m/McLachlan_BSSN.m

14 files changed, 214 insertions, 548 deletions

diff --git a/m/Makefile b/m/Makefile
index 48c7197..3af082a 100644
--- a/m/Makefile
+++ b/m/Makefile
@@ -6,60 +6,58 @@ all: McLachlan_ADM.out McLachlan_BSSN.out McLachlan_ADMConstraints.out McLachlan
 	@echo
 
 McLachlan_ADM.out: McLachlan_ADM.m
-	rm -rf ML_ADM*
+	rm -rf ML_ADM
 	./runmath.sh $^
-	for thorn in ML_ADM*; do \
-		./copy-if-changed.sh $$thorn ../$$thorn; \
+	for thorn in ML_ADM; do				\
+		./copy-if-changed.sh $$thorn ../$$thorn;	\
 	done
 
 McLachlan_BSSN.out: McLachlan_BSSN.m
 	rm -rf ML_BSSN*
 	./runmath.sh $^
-	for thorn in ML_BSSN*; do \
-		./copy-if-changed.sh $$thorn ../$$thorn; \
-		./create-helper-thorn.sh $$thorn; \
+	for thorn in ML_BSSN*; do					   \
+		./copy-if-changed.sh $$thorn ../$$thorn &&		   \
+		./create-helper-thorn.sh $$thorn &&			   \
 		./copy-if-changed.sh $${thorn}_Helper ../$${thorn}_Helper; \
 	done
 
 McLachlan_ADMConstraints.out: McLachlan_ADMConstraints.m
 	rm -rf ML_ADMConstraints*
 	./runmath.sh $^
-	perl -pi -e 's/MoL_PostStep/MoL_PseudoEvolution/g' ML_ADMConstraints*/schedule.ccl
-	for thorn in ML_ADMConstraints*; do \
-		./copy-if-changed.sh $$thorn ../$$thorn; \
+	for thorn in ML_ADMConstraints*; do			\
+		./copy-if-changed.sh $$thorn ../$$thorn;	\
 	done
 
 McLachlan_ADMQuantities.out: McLachlan_ADMQuantities.m
 	rm -rf ML_ADMQuantities*
 	./runmath.sh $^
-	perl -pi -e 's/MoL_PostStep/MoL_PseudoEvolution/g' ML_ADMQuantities*/schedule.ccl
-	for thorn in ML_ADMQuantities*; do \
-		./copy-if-changed.sh $$thorn ../$$thorn; \
+	for thorn in ML_ADMQuantities*; do			\
+		./copy-if-changed.sh $$thorn ../$$thorn;	\
 	done
 
 WaveToy.out: WaveToy.m
 	rm -rf ML_WaveToy
 	./runmath.sh $^
-	for thorn in ML_WaveToy; do \
-		./copy-if-changed.sh $$thorn ../$$thorn; \
+	for thorn in ML_WaveToy; do				\
+		./copy-if-changed.sh $$thorn ../$$thorn;	\
 	done
 
 WaveToyFO.out: WaveToyFO.m
 	rm -rf ML_WaveToyFO
 	./runmath.sh $^
-	for thorn in ML_WaveToyFO; do \
-		./copy-if-changed.sh $$thorn ../$$thorn; \
+	for thorn in ML_WaveToyFO; do				\
+		./copy-if-changed.sh $$thorn ../$$thorn;	\
 	done
 
 hydro.out: hydro.m
 	rm -rf ML_hydro
 	./runmath.sh $^
-	for thorn in ML_hydro; do \
-		./copy-if-changed.sh $$thorn ../$$thorn; \
+	for thorn in ML_hydro; do				\
+		./copy-if-changed.sh $$thorn ../$$thorn;	\
 	done
 
 clean:
-	rm -rf ML_ADM* ML_BSSN* ML_ADMConstraints*
+	rm -rf ML_ADM ML_BSSN* ML_ADMConstraints* ML_ADMQuantities*
 	rm -rf ML_WaveToy ML_WaveToyFO
 	rm -rf ML_hydro
 	rm -f McLachlan_ADM.out McLachlan_ADM.err
diff --git a/m/McLachlan_ADM.m b/m/McLachlan_ADM.m
index 663b8b3..92f0590 100644
--- a/m/McLachlan_ADM.m
+++ b/m/McLachlan_ADM.m
@@ -1,8 +1,3 @@
-$Path = Join[$Path, {"../../../repos/Kranc/Tools/CodeGen",
-                     "../../../repos/Kranc/Tools/MathematicaMisc"}];
- 
-Get["KrancThorn`"];
-
 
 SetEnhancedTimes[False];
 SetSourceLanguage["C"];
@@ -11,27 +6,25 @@ SetSourceLanguage["C"];
 (* Options *)
 (******************************************************************************)
 
-(* derivative order: 2, 4, 6, 8, ... *)
-derivOrder = 4;
-
 (* useJacobian: True or False *)
-useJacobian = False;
+useJacobian = True;
 
 (* timelevels: 2 or 3
    (keep this at 3; this is better chosen with a run-time parameter) *)
 evolutionTimelevels = 3;
 
 (* matter: 0 or 1 *)
-addMatter = 0;
+addMatter = 1;
 
 
 
 prefix = "ML_";
 suffix =
-  If [useJacobian, "_MP", ""] <>
-  If [derivOrder!=4, "_O" <> ToString[derivOrder], ""] <>
+  (* If [useJacobian, "_MP", ""] <> *)
+  (* If [derivOrder!=4, "_O" <> ToString[derivOrder], ""] <> *)
   If [evolutionTimelevels!=3, "_TL" <> ToString[evolutionTimelevels], ""] <>
-  If [addMatter!=0, "_M", ""];
+  (* If [addMatter!=0, "_M", ""] <> *)
+  "";
 
 ADM = prefix <> "ADM" <> suffix;
 
@@ -43,10 +36,10 @@ KD = KroneckerDelta;
 
 derivatives =
 {
-  PDstandardNth[i_]     -> StandardCenteredDifferenceOperator[1,derivOrder/2,i],
-  PDstandardNth[i_, i_] -> StandardCenteredDifferenceOperator[2,derivOrder/2,i],
-  PDstandardNth[i_, j_] -> StandardCenteredDifferenceOperator[1,derivOrder/2,i]
-                           StandardCenteredDifferenceOperator[1,derivOrder/2,j]
+  PDstandardNth[i_]     -> StandardCenteredDifferenceOperator[1,fdOrder/2,i],
+  PDstandardNth[i_, i_] -> StandardCenteredDifferenceOperator[2,fdOrder/2,i],
+  PDstandardNth[i_, j_] -> StandardCenteredDifferenceOperator[1,fdOrder/2,i]
+                           StandardCenteredDifferenceOperator[1,fdOrder/2,j]
 };
 
 PD = PDstandardNth;
@@ -278,7 +271,7 @@ constraintsCalc =
     gu[ua,ub] -> 1/detg detgExpr MatrixInverse[g[ua,ub]],
     G[ua,lb,lc] -> 1/2 gu[ua,ud]
                    (PD[g[lb,ld],lc] + PD[g[lc,ld],lb] - PD[g[lb,lc],ld]),
-    R[la,lb] -> G[uc,ld,la] G[lc,lb,ud] - G[uc,la,lb] G[lc,ld,ud]
+    R[la,lb] -> gu[us,ur](G[um,la,lr] G[uk,ls,lb] g[lk,lm] - G[um,la,lb] G[uk,ls,lr] g[lk,lm])
                 + 1/2 gu[uc,ud] (- PD[g[lc,ld],la,lb] + PD[g[lc,la],ld,lb]
                                  - PD[g[la,lb],lc,ld] + PD[g[ld,lb],lc,la]),
     trR -> R[la,lb] gu[ua,ub],
@@ -349,6 +342,15 @@ keywordParameters =
   }
 };
 
+intParameters =
+{
+  {
+    Name -> fdOrder,
+    Default -> 4,
+    AllowedValues -> {2,4,6,8}
+  }
+};
+
 (******************************************************************************)
 (* Construct the thorns *)
 (******************************************************************************)
@@ -371,6 +373,8 @@ CreateKrancThornTT [groups, ".", ADM,
   EvolutionTimelevels -> evolutionTimelevels,
   UseJacobian -> useJacobian,
   UseLoopControl -> True,
+  UseVectors -> True,
   InheritedImplementations -> inheritedImplementations,
-  KeywordParameters -> keywordParameters
+  KeywordParameters -> keywordParameters,
+  IntParameters -> intParameters
 ];
diff --git a/m/McLachlan_ADMConstraints.m b/m/McLachlan_ADMConstraints.m
index 1abb5ae..19346e7 100644
--- a/m/McLachlan_ADMConstraints.m
+++ b/m/McLachlan_ADMConstraints.m
@@ -1,7 +1,3 @@
-$Path = Join[$Path, {"../../../repos/Kranc/Tools/CodeGen",
-                     "../../../repos/Kranc/Tools/MathematicaMisc"}];
-
-Get["KrancThorn`"];
 
 SetEnhancedTimes[False];
 SetSourceLanguage["C"];
@@ -138,8 +134,9 @@ ADMConstraintsCalc =
   Name -> ADMConstraints,
   Schedule -> Automatic,
   Where -> Interior,
+  After -> "MoL_PostStep",
   Shorthands -> {detg, gu[ua,ub], G[ua,lb,lc],
-                 R[la,lb], trR, Km[la,lb], trK,
+                 R[la,lb], trR, Km[ua,lb], trK,
                  rho, S[la]},
   Equations -> 
   {
@@ -147,7 +144,7 @@ ADMConstraintsCalc =
     gu[ua,ub] -> 1/detg detgExpr MatrixInverse[g[ua,ub]],
     G[ua,lb,lc] -> 1/2 gu[ua,ud]
                    (PD[g[lb,ld],lc] + PD[g[lc,ld],lb] - PD[g[lb,lc],ld]),
-    R[la,lb] -> G[uc,ld,la] G[lc,lb,ud] - G[uc,la,lb] G[lc,ld,ud]
+    R[la,lb] -> gu[us,ur](G[um,la,lr] G[uk,ls,lb] g[lk,lm] - G[um,la,lb] G[uk,ls,lr] g[lk,lm])
                 + 1/2 gu[uc,ud] (- PD[g[lc,ld],la,lb] + PD[g[lc,la],ld,lb]
                                  - PD[g[la,lb],lc,ld] + PD[g[ld,lb],lc,la]),
     trR -> R[la,lb] gu[ua,ub],
diff --git a/m/McLachlan_ADMQuantities.m b/m/McLachlan_ADMQuantities.m
index dfa5d22..003a9f8 100644
--- a/m/McLachlan_ADMQuantities.m
+++ b/m/McLachlan_ADMQuantities.m
@@ -1,7 +1,3 @@
-$Path = Join[$Path, {"../../../repos/Kranc/Tools/CodeGen",
-                     "../../../repos/Kranc/Tools/MathematicaMisc"}];
-
-Get["KrancThorn`"];
 
 SetEnhancedTimes[False];
 SetSourceLanguage["C"];
@@ -212,6 +208,7 @@ ADMQuantitiesCalc =
   Name -> ADMQuantities,
   Schedule -> Automatic,
   Where -> Interior,
+  After -> "MoL_PostStep",
   Shorthands -> {detgt, gtu[ua,ub], dgtu[ua,ub,lc],
                  Gtl[la,lb,lc], Gtlu[la,lb,uc], Gt[ua,lb,lc],
                  Xtn[ua], Rt[la,lb], trRt,
diff --git a/m/McLachlan_BSSN.m b/m/McLachlan_BSSN.m
index fb35535..30cda03 100644
--- a/m/McLachlan_BSSN.m
+++ b/m/McLachlan_BSSN.m
@@ -1,7 +1,3 @@
-$Path = Join[$Path, {"../../../repos/Kranc/Tools/CodeGen",
-                     "../../../repos/Kranc/Tools/MathematicaMisc"}];
-
-Get["KrancThorn`"];
 
 SetEnhancedTimes[False];
 SetSourceLanguage["C"];
@@ -25,6 +21,8 @@ suffix =
 
 BSSN = prefix <> "BSSN" <> suffix;
 
+
+
 (******************************************************************************)
 (* Derivatives *)
 (******************************************************************************)
@@ -33,157 +31,54 @@ KD = KroneckerDelta;
 
 derivatives =
 {
-  PDstandardNth[i_]    -> StandardCenteredDifferenceOperator[1,derivOrder/2,i],
-  PDstandardNth[i_,i_] -> StandardCenteredDifferenceOperator[2,derivOrder/2,i],
-  PDstandardNth[i_,j_] -> StandardCenteredDifferenceOperator[1,derivOrder/2,i] *
-                          StandardCenteredDifferenceOperator[1,derivOrder/2,j],
+  PDstandardNth[i_]    -> StandardCenteredDifferenceOperator[1,fdOrder/2,i],
+  PDstandardNth[i_,i_] -> StandardCenteredDifferenceOperator[2,fdOrder/2,i],
+  PDstandardNth[i_,j_] -> StandardCenteredDifferenceOperator[1,fdOrder/2,i] *
+                          StandardCenteredDifferenceOperator[1,fdOrder/2,j],
   PDdissipationNth[i_] ->
-    spacing[i]^(derivOrder+1) / 2^(derivOrder+2) *
-    StandardCenteredDifferenceOperator[derivOrder+2,derivOrder/2+1,i],
+    spacing[i]^(fdOrder+1) / 2^(fdOrder+2) *
+    StandardCenteredDifferenceOperator[fdOrder+2,fdOrder/2+1,i],
   
 (* PD: These come from my mathematica notebook
    "Upwind-Kranc-Convert.nb" that converts upwinding finite
    differencing operators generated by
    StandardUpwindDifferenceOperator into this form *)
 
-  Switch[derivOrder,
-    2,
-    PDupwindNth[1] -> (dir[1]*(-3 + 4*shift[1]^dir[1] - 
-                         shift[1]^(2*dir[1])))/(2*spacing[1]),
-    4,
-    PDupwindNth[1] -> (dir[1]*(-10 - 3/shift[1]^dir[1] + 18*shift[1]^dir[1] -
-                         6*shift[1]^(2*dir[1]) + 
-                         shift[1]^(3*dir[1])))/(12*spacing[1]),
-    6,
-    PDupwindNth[1] -> (dir[1]*(-35 + 2/shift[1]^(2*dir[1]) - 
-                         24/shift[1]^dir[1] + 80*shift[1]^dir[1] -
-                         30*shift[1]^(2*dir[1]) + 8*shift[1]^(3*dir[1]) -
-                         shift[1]^(4*dir[1])))/(60*spacing[1]),
-    8,
-    PDupwindNth[1] -> (dir[1]*(-378 - 5/shift[1]^(3*dir[1]) +
-                         60/shift[1]^(2*dir[1]) - 420/shift[1]^dir[1] +
-                         1050*shift[1]^dir[1] - 420*shift[1]^(2*dir[1]) + 
-                         140*shift[1]^(3*dir[1]) - 30*shift[1]^(4*dir[1]) +
-                         3*shift[1]^(5*dir[1])))/(840*spacing[1])],
-  Switch[derivOrder,
-    2,
-    PDupwindNth[2] -> (dir[2]*(-3 + 4*shift[2]^dir[2] - 
-                         shift[2]^(2*dir[2])))/(2*spacing[2]),
-    4,
-    PDupwindNth[2] -> (dir[2]*(-10 - 3/shift[2]^dir[2] + 18*shift[2]^dir[2] -
-                         6*shift[2]^(2*dir[2]) + 
-                         shift[2]^(3*dir[2])))/(12*spacing[2]),
-    6,
-    PDupwindNth[2] -> (dir[2]*(-35 + 2/shift[2]^(2*dir[2]) - 
-                         24/shift[2]^dir[2] + 80*shift[2]^dir[2] -
-                         30*shift[2]^(2*dir[2]) + 8*shift[2]^(3*dir[2]) -
-                         shift[2]^(4*dir[2])))/(60*spacing[2]),
-    8,
-    PDupwindNth[2] -> (dir[2]*(-378 - 5/shift[2]^(3*dir[2]) +
-                         60/shift[2]^(2*dir[2]) - 420/shift[2]^dir[2] +
-                         1050*shift[2]^dir[2] - 420*shift[2]^(2*dir[2]) + 
-                         140*shift[2]^(3*dir[2]) - 30*shift[2]^(4*dir[2]) +
-                         3*shift[2]^(5*dir[2])))/(840*spacing[2])],
-  Switch[derivOrder,
-    2,
-    PDupwindNth[3] -> (dir[3]*(-3 + 4*shift[3]^dir[3] - 
-                         shift[3]^(2*dir[3])))/(2*spacing[3]),
-    4,
-    PDupwindNth[3] -> (dir[3]*(-10 - 3/shift[3]^dir[3] + 18*shift[3]^dir[3] -
-                         6*shift[3]^(2*dir[3]) + 
-                         shift[3]^(3*dir[3])))/(12*spacing[3]),
-    6,
-    PDupwindNth[3] -> (dir[3]*(-35 + 2/shift[3]^(2*dir[3]) - 
-                         24/shift[3]^dir[3] + 80*shift[3]^dir[3] -
-                         30*shift[3]^(2*dir[3]) + 8*shift[3]^(3*dir[3]) -
-                         shift[3]^(4*dir[3])))/(60*spacing[3]),
-    8,
-    PDupwindNth[3] -> (dir[3]*(-378 - 5/shift[3]^(3*dir[3]) +
-                         60/shift[3]^(2*dir[3]) - 420/shift[3]^dir[3] +
-                         1050*shift[3]^dir[3] - 420*shift[3]^(2*dir[3]) + 
-                         140*shift[3]^(3*dir[3]) - 30*shift[3]^(4*dir[3]) +
-                         3*shift[3]^(5*dir[3])))/(840*spacing[3])],
-
-  (* TODO: make these higher order stencils *)
-  PDonesided[1] -> dir[1] (-1 + shift[1]^dir[1]) / spacing[1],
-  PDonesided[2] -> dir[2] (-1 + shift[2]^dir[2]) / spacing[2],
-  PDonesided[3] -> dir[3] (-1 + shift[3]^dir[3]) / spacing[3]
+  Sequence@@Flatten[Table[
+   {PDupwindNth[i] -> Switch[fdOrder,
+      2, (dir[i]*(-3 + 4*shift[i]^dir[i] - shift[i]^(2*dir[i])))/(2*spacing[i]),
+      4, (dir[i]*(-10 - 3/shift[i]^dir[i] + 18*shift[i]^dir[i] -
+          6*shift[i]^(2*dir[i]) + shift[i]^(3*dir[i])))/(12*spacing[i]),
+      6, (dir[i]*(-35 + 2/shift[i]^(2*dir[i]) - 24/shift[i]^dir[i] + 80*shift[i]^dir[i] -
+          30*shift[i]^(2*dir[i]) + 8*shift[i]^(3*dir[i]) - shift[i]^(4*dir[i])))/(60*spacing[i]),
+      8, (dir[i]*(-378 - 5/shift[i]^(3*dir[i]) + 60/shift[i]^(2*dir[i]) - 420/shift[i]^dir[i] +
+          1050*shift[i]^dir[i] - 420*shift[i]^(2*dir[i]) + 140*shift[i]^(3*dir[i]) - 30*shift[i]^(4*dir[i]) +
+          3*shift[i]^(5*dir[i])))/(840*spacing[i])],
+
+    PDupwindNthAnti[i] -> Switch[fdOrder,
+      2, (+1 shift[i]^(-2) -4 shift[i]^(-1) +0 shift[i]^( 0) +4 shift[i]^(+1) -1 shift[i]^(+2)) / (4 spacing[i]),
+      4, (-1 shift[i]^(-3) +6 shift[i]^(-2) -21 shift[i]^(-1 )+0 shift[i]^( 0) +21 shift[i]^(+1)
+          -6 shift[i]^(+2) +1 shift[i]^(+3)) / (24 spacing[i]),
+      6, (+1 shift[i]^(-4) -8 shift[i]^(-3) +32 shift[i]^(-2) -104 shift[i]^(-1) +0 shift[i]^( 0)
+          +104 shift[i]^(+1) -32 shift[i]^(+2) +8 shift[i]^(+3) -1 shift[i]^(+4)) / (120 spacing[i]),
+      8, (-3 shift[i]^(-5) +30 shift[i]^(-4) -145 shift[i]^(-3) +480 shift[i]^(-2) -1470 shift[i]^(-1)
+          +0 shift[i]^( 0) +1470 shift[i]^(+1) -480 shift[i]^(+2) +145 shift[i]^(+3) -30 shift[i]^(+4)
+          +3 shift[i]^(+5)) / (1680 spacing[i])],
+
+    PDupwindNthSymm[i] -> Switch[fdOrder,
+     2, (-1 shift[i]^(-2) +4 shift[i]^(-1) -6 shift[i]^( 0) +4 shift[i]^(+1) -1 shift[i]^(+2)) / (4 spacing[i]),
+     4, (+1 shift[i]^(-3) -6 shift[i]^(-2) +15 shift[i]^(-1) -20 shift[i]^( 0) +15 shift[i]^(+1)
+         -6 shift[i]^(+2) +1 shift[i]^(+3)) / (24 spacing[i]),
+     6, (-1 shift[i]^(-4) +8 shift[i]^(-3) - 28 shift[i]^(-2)+56 shift[i]^(-1)-70 shift[i]^( 0)
+         +56 shift[i]^(+1) -28 shift[i]^(+2) +8 shift[i]^(+3) -1 shift[i]^(+4)) / (120 spacing[i]),
+     8, (+3 shift[i]^(-5) -30 shift[i]^(-4) +135 shift[i]^(-3) -360 shift[i]^(-2) +630 shift[i]^(-1)
+         -756 shift[i]^( 0) +630 shift[i]^(+1) -360 shift[i]^(+2) +135 shift[i]^(+3) -30 shift[i]^(+4)
+         +3 shift[i]^(+5)) / (1680 spacing[i])],
+
+    (* TODO: make these higher order stencils *)
+    PDonesided[i] -> dir[i] (-1 + shift[i]^dir[i]) / spacing[i]} /. i->j, {j,1,3}],1]
 };
 
-derivatives = Join[derivatives,
-  Flatten[Table[
-  Switch[derivOrder,
-    2,
-    {PDupwindNthAnti[i] -> (  +1 shift[i]^(-2)
-                              -4 shift[i]^(-1)
-                              +0 shift[i]^( 0)
-                              +4 shift[i]^(+1)
-                              -1 shift[i]^(+2)) / (4 spacing[i]),
-     PDupwindNthSymm[i] -> (  -1 shift[i]^(-2)
-                              +4 shift[i]^(-1)
-                              -6 shift[i]^( 0)
-                              +4 shift[i]^(+1)
-                              -1 shift[i]^(+2)) / (4 spacing[i])},
-    4,
-    {PDupwindNthAnti[i] -> (  -1 shift[i]^(-3)
-                              +6 shift[i]^(-2)
-                             -21 shift[i]^(-1)
-                              +0 shift[i]^( 0)
-                             +21 shift[i]^(+1)
-                              -6 shift[i]^(+2)
-                              +1 shift[i]^(+3)) / (24 spacing[i]),
-     PDupwindNthSymm[i] -> (  +1 shift[i]^(-3)
-                              -6 shift[i]^(-2)
-                             +15 shift[i]^(-1)
-                             -20 shift[i]^( 0)
-                             +15 shift[i]^(+1)
-                              -6 shift[i]^(+2)
-                              +1 shift[i]^(+3)) / (24 spacing[i])},
-    6,
-    {PDupwindNthAnti[i] -> (  +1 shift[i]^(-4)
-                              -8 shift[i]^(-3)
-                             +32 shift[i]^(-2)
-                            -104 shift[i]^(-1)
-                              +0 shift[i]^( 0)
-                            +104 shift[i]^(+1)
-                             -32 shift[i]^(+2)
-                              +8 shift[i]^(+3)
-                              -1 shift[i]^(+4)) / (120 spacing[i]),
-     PDupwindNthSymm[i] -> (  -1 shift[i]^(-4)
-                              +8 shift[i]^(-3)
-                             -28 shift[i]^(-2)
-                             +56 shift[i]^(-1)
-                             -70 shift[i]^( 0)
-                             +56 shift[i]^(+1)
-                             -28 shift[i]^(+2)
-                              +8 shift[i]^(+3)
-                              -1 shift[i]^(+4)) / (120 spacing[i])},
-    8,
-    {PDupwindNthAnti[i] -> (   -3 shift[i]^(-5)
-                              +30 shift[i]^(-4)
-                             -145 shift[i]^(-3)
-                             +480 shift[i]^(-2)
-                            -1470 shift[i]^(-1)
-                               +0 shift[i]^( 0)
-                            +1470 shift[i]^(+1)
-                             -480 shift[i]^(+2)
-                             +145 shift[i]^(+3)
-                              -30 shift[i]^(+4)
-                               +3 shift[i]^(+5)) / (1680 spacing[i]),
-     PDupwindNthSymm[i] -> (   +3 shift[i]^(-5)
-                              -30 shift[i]^(-4)
-                             +135 shift[i]^(-3)
-                             -360 shift[i]^(-2)
-                             +630 shift[i]^(-1)
-                             -756 shift[i]^( 0)
-                             +630 shift[i]^(+1)
-                             -360 shift[i]^(+2)
-                             +135 shift[i]^(+3)
-                              -30 shift[i]^(+4)
-                               +3 shift[i]^(+5)) / (1680 spacing[i])}],
-  {i,1,3}]]
-];
-
 PD     = PDstandardNth;
 PDu    = PDupwindNth;
 PDua   = PDupwindNthAnti;
@@ -323,9 +218,9 @@ extraGroups =
    {"TmunuBase::stress_energy_tensor", {eTxx, eTxy, eTxz, eTyy, eTyz, eTzz}}
 };
 
+groups = Join [declaredGroups, extraGroups];
 
 
-groups = Join [declaredGroups, extraGroups];
 
 (******************************************************************************)
 (* Initial data *)
@@ -473,11 +368,14 @@ initGammaCalc =
   Equations -> 
   {
     Xt[ua] -> 0,
+    A      -> 0,
+    B[ua]  -> 0,
     Tet    -> 0
   }
 };
 
 
+
 (******************************************************************************)
 (* Convert to ADMBase *)
 (******************************************************************************)
@@ -504,11 +402,15 @@ convertToADMBaseDtLapseShiftCalc =
   Schedule -> {"IN " <> BSSN <> "_convertToADMBaseGroup"},
   ConditionalOnKeyword -> {"dt_lapse_shift_method", "correct"},
   Where -> Interior,
-  Shorthands -> {dir[ua], eta, theta},
+  Shorthands -> {dir[ua], detgt, gtu[ua,ub], eta, theta},
   Equations -> 
   {
     dir[ua] -> Sign[beta[ua]],
     
+    detgt -> 1 (* detgtExpr *),
+    (* This leads to simpler code... *)
+    gtu[ua,ub]   -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
+    
     eta -> etaExpr,
     theta -> thetaExpr,
     
@@ -523,9 +425,17 @@ convertToADMBaseDtLapseShiftCalc =
                     (+ LapseACoeff       A
                      + (1 - LapseACoeff) (trK - 2 Tet))
                   + LapseAdvectionCoeff Upwind[beta[ua], alpha, la],
-    admdtbeta[ua] -> + theta ShiftGammaCoeff
-                       (+ ShiftBCoeff        B[ua]
-                        + (1 - ShiftBCoeff) (Xt[ua] - eta BetaDriver beta[ua]))
+    admdtbeta[ua] -> IfThen[harmonicShift,
+                            - 1/2 gtu[ua,uj] phi alpha
+                              (- 2 alpha PD[phi,lj]
+                               + 2 phi PD[alpha,lj]
+                               + gtu[uk,ul] phi alpha
+                                 (PD[gt[lk,ll],lj] - 2 PD[gt[lj,lk],ll])),
+                            (* else *)
+                            + theta ShiftGammaCoeff
+                              (+ ShiftBCoeff B[ua]
+                               + (1 - ShiftBCoeff)
+                                 (Xt[ua] - eta BetaDriver beta[ua]))]
                      + ShiftAdvectionCoeff Upwind[beta[ub], beta[ua], lb]
   }
 };
@@ -536,9 +446,13 @@ convertToADMBaseDtLapseShiftBoundaryCalc =
   Schedule -> {"IN " <> BSSN <> "_convertToADMBaseGroup"},
   ConditionalOnKeyword -> {"dt_lapse_shift_method", "correct"},
   Where -> BoundaryWithGhosts,
-  Shorthands -> {eta, theta},
+  Shorthands -> {detgt, gtu[ua,ub], eta, theta},
   Equations ->
   {
+    detgt -> 1 (* detgtExpr *),
+    (* This leads to simpler code... *)
+    gtu[ua,ub]   -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
+    
     eta -> etaExpr,
     theta -> thetaExpr,
     
@@ -551,9 +465,13 @@ convertToADMBaseDtLapseShiftBoundaryCalc =
     admdtalpha -> - harmonicF alpha^harmonicN
                     (+ LapseACoeff       A
                      + (1 - LapseACoeff) (trK - 2 Tet)),
-    admdtbeta[ua] -> + theta ShiftGammaCoeff
-                       (+ ShiftBCoeff        B[ua]
-                        + (1 - ShiftBCoeff) (Xt[ua] - eta BetaDriver beta[ua]))
+    admdtbeta[ua] -> IfThen[harmonicShift,
+                            0,
+                            (* else *)
+                            + theta ShiftGammaCoeff
+                              (+ ShiftBCoeff B[ua]
+                               + (1 - ShiftBCoeff)
+                                 (Xt[ua] - eta BetaDriver beta[ua]))]
   }
 };
 
@@ -563,9 +481,13 @@ convertToADMBaseFakeDtLapseShiftCalc =
   Schedule -> {"IN " <> BSSN <> "_convertToADMBaseGroup"},
   ConditionalOnKeyword -> {"dt_lapse_shift_method", "noLapseShiftAdvection"},
   Where -> Everywhere,
-  Shorthands -> {eta, theta},
+  Shorthands -> {detgt, gtu[ua,ub], eta, theta},
   Equations ->
   {
+    detgt -> 1 (* detgtExpr *),
+    (* This leads to simpler code... *)
+    gtu[ua,ub]   -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
+    
     eta -> etaExpr,
     theta -> thetaExpr,
     
@@ -580,9 +502,13 @@ convertToADMBaseFakeDtLapseShiftCalc =
     admdtalpha -> - harmonicF alpha^harmonicN
                     (+ LapseACoeff       A
                      + (1 - LapseACoeff) (trK - 2 Tet)),
-    admdtbeta[ua] -> + theta ShiftGammaCoeff
-                       (+ ShiftBCoeff        B[ua]
-                        + (1 - ShiftBCoeff) (Xt[ua] - eta BetaDriver beta[ua]))
+    admdtbeta[ua] -> IfThen[harmonicShift,
+                            0,
+                            (* else *)
+                            + theta ShiftGammaCoeff
+                              (+ ShiftBCoeff B[ua]
+                               + (1 - ShiftBCoeff)
+                                 (Xt[ua] - eta BetaDriver beta[ua]))]
   }
 };
 
@@ -784,9 +710,17 @@ evolCalc =
     
     (* dot[beta[ua]] -> eta Xt[ua], *)
     (* dot[beta[ua]] -> ShiftGammaCoeff alpha^ShiftAlphaPower B[ua], *)
-    dot[beta[ua]] -> + theta ShiftGammaCoeff
-                       (+ ShiftBCoeff        B[ua]
-                        + (1 - ShiftBCoeff) (Xt[ua] - eta BetaDriver beta[ua])),
+    dot[beta[ua]] -> IfThen[harmonicShift,
+                            - 1/2 gtu[ua,uj] phi alpha
+                              (- 2 alpha PD[phi,lj]
+                               + 2 phi PD[alpha,lj]
+                               + gtu[uk,ul] phi alpha
+                                 (PD[gt[lk,ll],lj] - 2 PD[gt[lj,lk],ll])),
+                            (* else *)
+                            + theta ShiftGammaCoeff
+                              (+ ShiftBCoeff B[ua]
+                               + (1 - ShiftBCoeff)
+                                 (Xt[ua] - eta BetaDriver beta[ua]))],
 
     dot[B[ua]]    -> + ShiftBCoeff (dotXt[ua] - eta BetaDriver B[ua])
 
@@ -859,7 +793,6 @@ evolCalc2 = PartialCalculation[evolCalc, "2",
     dot[At[la,lb]]
   }];
 
-
 dissCalc =
 {
   Name -> BSSN <> "_Dissipation",
@@ -877,6 +810,24 @@ dissCalc =
   }
 };
 
+dissCalcs =
+Table[
+{
+  Name -> BSSN <> "_Dissipation_" <> ToString[var /. {Tensor[n_,__] -> n}],
+  Schedule -> {"IN " <> BSSN <> "_evolCalcGroup " <>
+               "AFTER (" <> BSSN <> "_RHS1 " <> BSSN <> "_RHS2)"},
+  ConditionalOnKeyword -> {"apply_dissipation", "always"},
+  Where -> InteriorNoSync,
+  Shorthands -> {epsdiss[ua]},
+  Equations ->
+  {
+    epsdiss[ua] -> EpsDiss,
+    dot[var]    -> dot[var] + epsdiss[ux] PDdiss[var,lx]
+  }
+},
+  {var, {phi, gt[la,lb], Xt[ui], trK, At[la,lb], alpha, A, beta[ua], B[ua]}}
+];
+
 RHSStaticBoundaryCalc =
 {
   Name -> BSSN <> "_RHSStaticBoundary",
@@ -940,7 +891,7 @@ RHSRadiativeBoundaryCalc =
     em4phi     -> IfThen [conformalMethod, phi^2, Exp[-4 phi]],
     gu[ua,ub]  -> em4phi gtu[ua,ub],
     
-    nn[la] -> normal[la],
+    nn[la] -> Euc[la,lb] normal[ub],
     nu[ua] -> gu[ua,ub] nn[lb],
     nlen2  -> nu[ua] nn[la],
     nlen   -> Sqrt[nlen2],
@@ -964,11 +915,14 @@ RHSRadiativeBoundaryCalc =
 enforceCalc =
 {
   Name -> BSSN <> "_enforce",
-  Schedule -> {"IN MoL_PostStep BEFORE " <> BSSN <> "_SelectBoundConds"},
+  Schedule -> {"IN MoL_PostStepModify"},
   Shorthands -> {detgt, gtu[ua,ub], trAt},
   Equations -> 
   {
-    (* Enforcing the constraints needs to be a projection, because it
+    (* The following comment is still interesting, but is not correct
+       any more since it is now scheduled in MoL_PostStepModify instead:
+
+       Enforcing the constraints needs to be a projection, because it
        is applied in MoL_PostStep and may thus be applied multiple
        times, not only during time evolution. Therefore detgt has to
        be calculated correctly, without assuming that det gt_ij = 1,
@@ -1018,12 +972,13 @@ constraintsCalc =
 {
   Name -> BSSN <> "_constraints",
   Schedule -> Automatic,
+  After -> "MoL_PostStep",
   Where -> Interior,
   Shorthands -> {detgt, ddetgt[la], gtu[ua,ub],
                  Gt[ua,lb,lc], Gtl[la,lb,lc], Gtlu[la,lb,uc], Xtn[ua],
                  e4phi, em4phi,
                  g[la,lb], detg, gu[ua,ub], ddetg[la], G[ua,lb,lc],
-                 Rt[la,lb], Rphi[la,lb], R[la,lb], trR, Atm[la,lb],
+                 Rt[la,lb], Rphi[la,lb], R[la,lb], trR, Atm[ua,lb],
                  gK[la,lb,lc], cdphi[la], cdphi2[la,lb],
                  rho, S[la], fac1, fac2, Zet[ua]},
   Equations -> 
@@ -1285,6 +1240,18 @@ intParameters =
     AllowedValues -> {{Value -> "0", Description -> "phi method"},
                       {Value -> "1", Description -> "W method"}},
     Default -> 0
+  },
+  {
+    Name -> fdOrder,
+    Default -> derivOrder,
+    AllowedValues -> {2,4,6,8}
+  },
+  {
+    Name -> harmonicShift,
+    Description -> "Whether to use the harmonic shift",
+    AllowedValues -> {{Value -> "0", Description -> "Gamma driver shift"},
+                      {Value -> "1", Description -> "Harmonic shift"}},
+    Default -> 0
   }
 };
 
@@ -1319,7 +1286,7 @@ realParameters =
   },
   {
     Name -> LapseAdvectionCoeff,
-    Description -> "Factor in front of the shift advection terms in 1+log",
+    Description -> "Factor in front of the lapse advection terms in 1+log",
     Default -> 1
   },
   {
@@ -1377,6 +1344,7 @@ realParameters =
 (******************************************************************************)
 
 calculations =
+Join[
 {
   initialCalc,
   convertFromADMBaseCalc,
@@ -1398,7 +1366,9 @@ calculations =
   convertToADMBaseFakeDtLapseShiftCalc,
   (* constraintsCalc, *)
   constraintsCalc1, constraintsCalc2
-};
+},
+  {} (*dissCalcs*)
+];
 
 CreateKrancThornTT [groups, ".", BSSN,
   Calculations -> calculations,
@@ -1406,7 +1376,7 @@ CreateKrancThornTT [groups, ".", BSSN,
   PartialDerivatives -> derivatives,
   EvolutionTimelevels -> evolutionTimelevels,
   DefaultEvolutionTimelevels -> 3,
-  UseJacobian -> useJacobian,
+  UseJacobian -> True,
   UseLoopControl -> True,
   UseVectors -> True,
   InheritedImplementations -> inheritedImplementations,
@@ -1427,6 +1397,7 @@ CreateKrancThornTT [groups, ".", BSSN,
 
 (* derivative order: 2, 4, 6, 8, ... *)
 (* useJacobian: False or True *)
+(* split upwind derivatives: False or True *)
 (* timelevels: 2 or 3
    (keep this at 3; this is better chosen with a run-time parameter) *)
 (* matter: 0 or 1
diff --git a/m/McLachlan_BSSN_Peter.m b/m/McLachlan_BSSN_Peter.m
index 3631867..c5b5b7d 100644
--- a/m/McLachlan_BSSN_Peter.m
+++ b/m/McLachlan_BSSN_Peter.m
@@ -1,7 +1,3 @@
-$Path = Join[$Path, {"../../../repos/Kranc/Tools/CodeGen",

-                     "../../../repos/Kranc/Tools/MathematicaMisc"}];

-

-Get["KrancThorn`"];

 

 SetEnhancedTimes[False];

 SetSourceLanguage["C"];

diff --git a/m/McLachlan_WeylScalars.m b/m/McLachlan_WeylScalars.m
index 86b4a36..5e8d805 100644
--- a/m/McLachlan_WeylScalars.m
+++ b/m/McLachlan_WeylScalars.m
@@ -1,7 +1,3 @@
-$Path = Join[$Path, {"../../../repos/Kranc/Tools/CodeGen",
-                     "../../../repos/Kranc/Tools/MathematicaMisc"}];
-
-Get["KrancThorn`"];
 
 SetEnhancedTimes[False];
 SetSourceLanguage["C"];
diff --git a/m/McLachlantmp.m b/m/McLachlantmp.m
index 1d652c9..54581e2 100644
--- a/m/McLachlantmp.m
+++ b/m/McLachlantmp.m
@@ -1,7 +1,3 @@
-$Path = Join[$Path, {"../../../repos/Kranc/Tools/CodeGen",
-                     "../../../repos/Kranc/Tools/MathematicaMisc"}];
-
-Get["KrancThorn`"];
 
 SetEnhancedTimes[False];
 SetSourceLanguage["C"];
diff --git a/m/WaveToy.m b/m/WaveToy.m
index dfc0862..de2c73e 100644
--- a/m/WaveToy.m
+++ b/m/WaveToy.m
@@ -1,7 +1,3 @@
-$Path = Join[$Path, {"../../../repos/Kranc/Tools/CodeGen",
-                     "../../../repos/Kranc/Tools/MathematicaMisc"}];
-
-Get["KrancThorn`"];
 
 SetEnhancedTimes[False];
 SetSourceLanguage["C"];
diff --git a/m/WaveToyFO.m b/m/WaveToyFO.m
index 4864131..f80b0b2 100644
--- a/m/WaveToyFO.m
+++ b/m/WaveToyFO.m
@@ -1,7 +1,3 @@
-$Path = Join[$Path, {"../../../repos/Kranc/Tools/CodeGen",
-                     "../../../repos/Kranc/Tools/MathematicaMisc"}];
-
-Get["KrancThorn`"];
 
 SetEnhancedTimes[False];
 SetSourceLanguage["C"];
diff --git a/m/WaveToyMP.m b/m/WaveToyMP.m
deleted file mode 100644
index 546a618..0000000
--- a/m/WaveToyMP.m
+++ /dev/null
@@ -1,271 +0,0 @@
-$Path = Join[$Path, {"../../../repos/Kranc/Tools/CodeGen",
-                     "../../../repos/Kranc/Tools/MathematicaMisc"}];
-
-Get["KrancThorn`"];
-
-SetEnhancedTimes[False];
-SetSourceLanguage["C"];
-
-(******************************************************************************)
-(* Derivatives *)
-(******************************************************************************)
-
-derivOrder = 4;
-
-derivatives =
-{
-  PDstandardNth[i_]     -> StandardCenteredDifferenceOperator[1,derivOrder/2,i],
-  PDstandardNth[i_, i_] -> StandardCenteredDifferenceOperator[2,derivOrder/2,i],
-  PDstandardNth[i_, j_] -> StandardCenteredDifferenceOperator[1,derivOrder/2,i]
-                           StandardCenteredDifferenceOperator[1,derivOrder/2,j],
-  PDstandardNth[i_, i_, i_] ->
-    StandardCenteredDifferenceOperator[3,derivOrder/2,i],
-  PDstandardNth[i_, i_, j_] ->
-    StandardCenteredDifferenceOperator[2,derivOrder/2,i]
-    StandardCenteredDifferenceOperator[1,derivOrder/2,j],
-  PDstandardNth[i_, j_, i_] ->
-    StandardCenteredDifferenceOperator[2,derivOrder/2,i]
-    StandardCenteredDifferenceOperator[1,derivOrder/2,j],
-  PDstandardNth[j_, i_, i_] ->
-    StandardCenteredDifferenceOperator[2,derivOrder/2,i]
-    StandardCenteredDifferenceOperator[1,derivOrder/2,j],
-  PDstandardNth[i_, j_, k_] ->
-    StandardCenteredDifferenceOperator[1,derivOrder/2,i]
-    StandardCenteredDifferenceOperator[1,derivOrder/2,j]
-    StandardCenteredDifferenceOperator[1,derivOrder/2,k]
-};
-
-(* local derivatives *)
-PDloc = PDstandardNth;
-
-(* global derivatives *)
-PDglob[var_,lx_] := Jinv[u1,lx] PDloc[var,l1];
-PDglob[var_,lx_,ly_] :=
-  dJinv[u1,lx,ly] PDloc[var,l1] + Jinv[u1,lx] Jinv[u2,ly] PDloc[var,l1,l2];
-
-UseGlobalDerivs = False;
-PD := If [UseGlobalDerivs, PDglob, PDloc];
-
-(* timelevels *)
-evolutionTimelevels = 2;
-
-KD = KroneckerDelta;
-
-(******************************************************************************)
-(* Tensors *)
-(******************************************************************************)
-
-(* Register the tensor quantities with the TensorTools package *)
-Map [DefineTensor, {u, rho, v, w, J, Jinv, dJ, dJinv}];
-
-AssertSymmetricIncreasing [dJ[uA,lb,lc], lb,lc];
-AssertSymmetricIncreasing [dJinv[ua,lx,ly], lx,ly];
-
-(******************************************************************************)
-(* Groups *)
-(******************************************************************************)
-
-evolvedGroups =
-  {SetGroupName [CreateGroupFromTensor [u  ], "WT_u"  ],
-   SetGroupName [CreateGroupFromTensor [rho], "WT_rho"]};
-evaluatedGroups = {};
-
-declaredGroups = Join [evolvedGroups, evaluatedGroups];
-declaredGroupNames = Map [First, declaredGroups];
-
-evolvedGroupsFO =
-  {SetGroupName [CreateGroupFromTensor [u    ], "WT_u"  ],
-   SetGroupName [CreateGroupFromTensor [v[la]], "WT_v"  ],
-   SetGroupName [CreateGroupFromTensor [rho  ], "WT_rho"]};
-evaluatedGroupsFO =
-  {SetGroupName [CreateGroupFromTensor [w[ua]], "WT_w"]};
-
-declaredGroupsFO = Join [evolvedGroupsFO, evaluatedGroupsFO];
-declaredGroupNamesFO = Map [First, declaredGroupsFO];
-
-
-
-extraGroups =
-  {{"MultiPatch::transformation",
-    {dxda,dxdb,dxdc, dyda,dydb,dydc, dzda,dzdb,dzdc}},
-   {"MultiPatch::transformation_inv",
-    {dadx,dady,dadz, dbdx,dbdy,dbdz, dcdx,dcdy,dcdz}},
-   {"MultiPatch::transformation_derivs",
-    {ddxdada,ddxdadb,ddxdadc,ddxdbdb,ddxdbdc,ddxdcdc,
-     ddydada,ddydadb,ddydadc,ddydbdb,ddydbdc,ddydcdc,
-     ddzdada,ddzdadb,ddzdadc,ddzdbdb,ddzdbdc,ddzdcdc}},
-   {"MultiPatch::transformation_inv_derivs",
-    {ddadxdx,ddadxdy,ddadxdz,ddadydy,ddadzdz,ddadydz,
-     ddbdxdx,ddbdxdy,ddbdxdz,ddbdydy,ddbdzdz,ddbdydz,
-     ddcdxdx,ddcdxdy,ddcdxdz,ddcdydy,ddcdzdz,ddcdydz}}};
-
-
-
-groups = Join [declaredGroups, extraGroups];
-groupsFO = Join [declaredGroupsFO, extraGroups];
-
-(******************************************************************************)
-(* Initial data *)
-(******************************************************************************)
-
-initialCalc =
-{
-  Name -> "WT_Gaussian",
-  Schedule -> {"AT initial"},
-  (* Where -> Boundary, *)
-  (* Where -> Interior, *)
-  Equations -> 
-  {
-    u -> 0,
-    rho -> 0
-  }
-};
-
-initialCalcFO =
-{
-  Name -> "WTFO_Gaussian",
-  Schedule -> {"AT initial"},
-  (* Where -> Boundary, *)
-  (* Where -> Interior, *)
-  Equations -> 
-  {
-    u -> 0,
-    v[la] -> 0,
-    rho -> 0
-  }
-};
-
-(******************************************************************************)
-(* Evolution equations *)
-(******************************************************************************)
-
-evolCalc =
-{
-  Name -> "WT_RHS",
-  Schedule -> {"IN MoL_CalcRHS", "AT analysis"},
-  Where -> Interior,
-  Shorthands -> {Jinv[ua,lx], dJinv[ua,lx,ly]},
-  Equations -> 
-  {
-    Jinv11 -> dadx,
-    Jinv12 -> dady,
-    Jinv13 -> dadz,
-    Jinv21 -> dbdx,
-    Jinv22 -> dbdy,
-    Jinv23 -> dbdz,
-    Jinv31 -> dcdx,
-    Jinv32 -> dcdy,
-    Jinv33 -> dcdz,
-    dJinv111 -> ddadxdx,
-    dJinv112 -> ddadxdy,
-    dJinv113 -> ddadxdz,
-    dJinv122 -> ddadydy,
-    dJinv123 -> ddadydz,
-    dJinv133 -> ddadzdz,
-    dJinv211 -> ddadxdx,
-    dJinv212 -> ddadxdy,
-    dJinv213 -> ddadxdz,
-    dJinv222 -> ddadydy,
-    dJinv223 -> ddadydz,
-    dJinv233 -> ddadzdz,
-    dJinv311 -> ddadxdx,
-    dJinv312 -> ddadxdy,
-    dJinv313 -> ddadxdz,
-    dJinv322 -> ddadydy,
-    dJinv323 -> ddadydz,
-    dJinv333 -> ddadzdz,
-    dot[u] -> rho,
-    dot[rho] -> KD[ua,ub] PD[u,la,lb]
-  }
-};
-
-evolCalcFO =
-{
-  Name -> "WTFO_RHS",
-  Schedule -> {"IN MoL_CalcRHS", "AT analysis"},
-  Where -> Interior,
-  Shorthands -> {Jinv[ua,lx]},
-  Equations -> 
-  {
-    Jinv11 -> dadx,
-    Jinv12 -> dady,
-    Jinv13 -> dadz,
-    Jinv21 -> dbdx,
-    Jinv22 -> dbdy,
-    Jinv23 -> dbdz,
-    Jinv31 -> dcdx,
-    Jinv32 -> dcdy,
-    Jinv33 -> dcdz,
-    dot[u] -> rho,
-    dot[rho] -> KD[ua,ub] PD[v[la],lb],
-    dot[v[la]] -> PD[rho,la]
-  }
-};
-
-(******************************************************************************)
-(* Constraint equations *)
-(******************************************************************************)
-
-constraintsCalcFO =
-{
-  Name -> "WTFO_constraints",
-  Schedule -> {"AT analysis"},
-  Where -> Interior,
-  Shorthands -> {Jinv[ua,lx]},
-  Equations -> 
-  {
-    Jinv11 -> dadx,
-    Jinv12 -> dady,
-    Jinv13 -> dadz,
-    Jinv21 -> dbdx,
-    Jinv22 -> dbdy,
-    Jinv23 -> dbdz,
-    Jinv31 -> dcdx,
-    Jinv32 -> dcdy,
-    Jinv33 -> dcdz,
-    w[ua] -> Eps[ua,ub,uc] PD[v[lb],lc]
-  }
-};
-
-(******************************************************************************)
-(* Implementations *)
-(******************************************************************************)
-
-inheritedImplementations = {"MultiPatch"};
-
-(******************************************************************************)
-(* Construct the thorns *)
-(******************************************************************************)
-
-calculations = 
-{
-  initialCalc,
-  evolCalc
-};
-
-CreateKrancThornTT [groups, ".", "ML_WaveToy",
-  Calculations -> calculations,
-  DeclaredGroups -> declaredGroupNames,
-  PartialDerivatives -> derivatives,
-  UseLoopControl -> True,
-  EvolutionTimelevels -> evolutionTimelevels,
-  InheritedImplementations -> inheritedImplementations
-];
-
-
-
-calculationsFO = 
-{
-  initialCalcFO,
-  evolCalcFO,
-  constraintsCalcFO
-};
-
-CreateKrancThornTT [groupsFO, ".", "ML_FOWaveToy",
-  Calculations -> calculationsFO,
-  DeclaredGroups -> declaredGroupNamesFO,
-  PartialDerivatives -> derivatives,
-  UseLoopControl -> True,
-  EvolutionTimelevels -> evolutionTimelevels,
-  InheritedImplementations -> inheritedImplementations
-];
diff --git a/m/hydro.m b/m/hydro.m
index 85b9e45..a69194f 100644
--- a/m/hydro.m
+++ b/m/hydro.m
@@ -113,8 +113,8 @@ prim2conCalc =
   {
     vol     -> h^3,
     mass    -> vol rho,
-    mom[la] -> mass vel[la],
-    ene     -> (1/2) mass vel[ua] vel[la] + mass eps
+    mom[la] -> mass Euc[la,lb] vel[ub],
+    ene     -> (1/2) mass Euc[la,lb] vel[ub] vel[ua] + mass eps
   }
 };
 
@@ -129,8 +129,8 @@ con2primCalc =
   Equations -> 
   {
     rho     -> mass / vol,
-    vel[ua] -> mom[ua] / mass,
-    eps     -> ene / mass - (1/2) vel[ua] vel[la],
+    vel[ua] -> Euc[ua,ub] mom[lb] / mass,
+    eps     -> ene / mass - (1/2) Euc[la,lc] vel[ua] vel[uc],
     
     press   -> Gamma rho eps
   }
diff --git a/m/prototype/ML_BSSN_Helper/src/SetGroupTags.c b/m/prototype/ML_BSSN_Helper/src/SetGroupTags.c
index 4b24e56..7ae24bb 100644
--- a/m/prototype/ML_BSSN_Helper/src/SetGroupTags.c
+++ b/m/prototype/ML_BSSN_Helper/src/SetGroupTags.c
@@ -15,29 +15,30 @@ ML_BSSN_SetGroupTags (void)
 {
   DECLARE_CCTK_PARAMETERS;
   
-  set_group_tags (0, 0, 1, "ADMBase::metric");
-  set_group_tags (0, 0, 1, "ADMBase::curv");
-  set_group_tags (0, 0, 1, "ADMBase::lapse");
-  set_group_tags (0, 0, 1, "ADMBase::shift");
-  set_group_tags (0, 0, 1, "ADMBase::dtlapse");
-  set_group_tags (0, 0, 1, "ADMBase::dtshift");
+  int const checkpoint = timelevels > 1;
+  set_group_tags (checkpoint, checkpoint, 1, "ADMBase::metric");
+  set_group_tags (checkpoint, checkpoint, 1, "ADMBase::curv");
+  set_group_tags (checkpoint, checkpoint, 1, "ADMBase::lapse");
+  set_group_tags (checkpoint, checkpoint, 1, "ADMBase::shift");
+  set_group_tags (checkpoint, checkpoint, 1, "ADMBase::dtlapse");
+  set_group_tags (checkpoint, checkpoint, 1, "ADMBase::dtshift");
   
-  set_group_tags (0, 0, 0, "ML_BSSN::ML_cons_detg");
-  set_group_tags (0, 0, 0, "ML_BSSN::ML_cons_Gamma");
-  set_group_tags (0, 0, 0, "ML_BSSN::ML_cons_traceA");
-  set_group_tags (0, 0, 0, "ML_BSSN::ML_Ham");
-  set_group_tags (0, 0, 0, "ML_BSSN::ML_mom");
+  set_group_tags (checkpoint, checkpoint, 0, "ML_BSSN::ML_cons_detg");
+  set_group_tags (checkpoint, checkpoint, 0, "ML_BSSN::ML_cons_Gamma");
+  set_group_tags (checkpoint, checkpoint, 0, "ML_BSSN::ML_cons_traceA");
+  set_group_tags (checkpoint, checkpoint, 0, "ML_BSSN::ML_Ham");
+  set_group_tags (checkpoint, checkpoint, 0, "ML_BSSN::ML_mom");
   
-  int const checkpoint = rhs_timelevels > 1;
-  set_group_tags (checkpoint, checkpoint, 0, "ML_BSSN::ML_log_confacrhs");
-  set_group_tags (checkpoint, checkpoint, 0, "ML_BSSN::ML_metricrhs");
-  set_group_tags (checkpoint, checkpoint, 0, "ML_BSSN::ML_Gammarhs");
-  set_group_tags (checkpoint, checkpoint, 0, "ML_BSSN::ML_trace_curvrhs");
-  set_group_tags (checkpoint, checkpoint, 0, "ML_BSSN::ML_curvrhs");
-  set_group_tags (checkpoint, checkpoint, 0, "ML_BSSN::ML_lapserhs");
-  set_group_tags (checkpoint, checkpoint, 0, "ML_BSSN::ML_dtlapserhs");
-  set_group_tags (checkpoint, checkpoint, 0, "ML_BSSN::ML_shiftrhs");
-  set_group_tags (checkpoint, checkpoint, 0, "ML_BSSN::ML_dtshiftrhs");
+  int const rhs_checkpoint = rhs_timelevels > 1;
+  set_group_tags (rhs_checkpoint, rhs_checkpoint, 0, "ML_BSSN::ML_log_confacrhs");
+  set_group_tags (rhs_checkpoint, rhs_checkpoint, 0, "ML_BSSN::ML_metricrhs");
+  set_group_tags (rhs_checkpoint, rhs_checkpoint, 0, "ML_BSSN::ML_Gammarhs");
+  set_group_tags (rhs_checkpoint, rhs_checkpoint, 0, "ML_BSSN::ML_trace_curvrhs");
+  set_group_tags (rhs_checkpoint, rhs_checkpoint, 0, "ML_BSSN::ML_curvrhs");
+  set_group_tags (rhs_checkpoint, rhs_checkpoint, 0, "ML_BSSN::ML_lapserhs");
+  set_group_tags (rhs_checkpoint, rhs_checkpoint, 0, "ML_BSSN::ML_dtlapserhs");
+  set_group_tags (rhs_checkpoint, rhs_checkpoint, 0, "ML_BSSN::ML_shiftrhs");
+  set_group_tags (rhs_checkpoint, rhs_checkpoint, 0, "ML_BSSN::ML_dtshiftrhs");
   
   return 0;
 }
diff --git a/m/runmath.sh b/m/runmath.sh
index 306b2e1..18196d1 100755
--- a/m/runmath.sh
+++ b/m/runmath.sh
@@ -3,8 +3,6 @@
 # Abort on errors
 set -e
 
-MATHEMATICA="math"
-
 script=$1
 
 if test -z "$script"; then
@@ -18,17 +16,8 @@ output=$(basename $script .m).out
 
 rm -f $output
 
-# Run Mathematica to regenerate the code
-< $script "$MATHEMATICA" | tee $error
-
-if grep 'KrancError' $error; then
-    echo
-    echo "There was an error when running Kranc on $script."
-    echo "The file $error contains details."
-    echo
-    echo "*** The Cactus thorns have NOT been updated. ***"
-    echo
-    exit 1
-fi
+# Run Kranc to regenerate the code
+../../../repos/Kranc/Bin/kranc $script | tee $error
+[ $PIPESTATUS -eq 0 ] || exit $PIPESTATUS
 
 mv $error $output