Tests: Add McLachlan.mt

This exercises a large amount of code, but takes about 2 minutes to run, even with NoSimplify -> True.

2 files changed, 1461 insertions, 1 deletions

diff --git a/Tests/McLachlan.mt b/Tests/McLachlan.mt
new file mode 100644
index 0000000..5eeea8d
--- /dev/null
+++ b/Tests/McLachlan.mt
@@ -0,0 +1,1459 @@
+
+(* Mathematica Test File *)
+
+SetEnhancedTimes[False];
+SetSourceLanguage["C"];
+
+(******************************************************************************)
+(* Options *)
+(******************************************************************************)
+
+createCode[derivOrder_, useJacobian_, splitUpwindDerivs_, evolutionTimelevels_, addMatter_, formulation_] :=
+Module[{prefix, suffix, thorn},
+
+prefix = "ML_";
+suffix =
+  ""
+  <> If [useJacobian, "_MP", ""]
+  <> If [derivOrder!=4, "_O" <> ToString[derivOrder], ""]
+  <> If [splitUpwindDerivs, "", "_UPW"]
+  (* <> If [evolutionTimelevels!=3, "_TL" <> ToString[evolutionTimelevels], ""] *)
+  (* <> If [addMatter==1, "_M", ""] *)
+  ;
+
+thorn = prefix <> formulation <> suffix;
+
+SetAttributes[IfCCZ4, HoldAll];
+IfCCZ4[expr_, else_:Sequence[]] := If[formulation === "CCZ4", expr, Unevaluated[else]];
+
+(******************************************************************************)
+(* Derivatives *)
+(******************************************************************************)
+
+KD = KroneckerDelta;
+
+derivatives =
+{
+  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_] ->
+    (-1)^(fdOrder/2) *
+    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 *)
+
+  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]
+};
+
+PD     = PDstandardNth;
+PDu    = PDupwindNth;
+PDua   = PDupwindNthAnti;
+PDus   = PDupwindNthSymm;
+(* PDo    = PDonesided; *)
+PDdiss = PDdissipationNth;
+
+If [splitUpwindDerivs,
+    Upwind[dir_, var_, idx_] := dir PDua[var,idx] + Abs[dir] PDus[var,idx],
+    Upwind[dir_, var_, idx_] := dir PDu[var,idx]];
+
+
+
+(******************************************************************************)
+(* Tensors *)
+(******************************************************************************)
+
+(* Register the tensor quantities with the TensorTools package *)
+Map [DefineTensor,
+     {normal, tangentA, tangentB, dir,
+      nn, nu, nlen, nlen2, su, vg,
+      xx, rr, th, ph,
+      admg, admK, admalpha, admdtalpha, admbeta, admdtbeta, H, M,
+      g, detg, gu, G, R, trR, Km, trK, cdphi, cdphi2,
+      phi, gt, At, Xt, Xtn, Theta, Z,
+      alpha, A, beta, B, Atm, Atu, trA, Ats, trAts,
+      dottrK, dotXt,
+      cXt, cS, cA,
+      e4phi, em4phi, ddetg, detgt, gtu, ddetgt, dgtu, ddgtu, Gtl, Gtlu, Gt, Ddetgt,
+      Rt, Rphi, gK,
+      T00, T0, T, rho, S,
+      x, y, z, r,
+      epsdiss}];
+
+(* NOTE: It seems as if Lie[.,.] did not take these tensor weights
+   into account.  Presumably, CD[.,.] and CDt[.,.] don't do this either.  *)
+SetTensorAttribute[phi, TensorWeight, +1/6];
+SetTensorAttribute[gt,  TensorWeight, -2/3];
+SetTensorAttribute[Xt,  TensorWeight, +2/3];
+SetTensorAttribute[At,  TensorWeight, -2/3];
+SetTensorAttribute[cXt, TensorWeight, +2/3];
+SetTensorAttribute[cS,  TensorWeight, +2  ];
+
+Map [AssertSymmetricIncreasing,
+     {admg[la,lb], admK[la,lb], g[la,lb], K[la,lb], R[la,lb], cdphi2[la,lb],
+      gt[la,lb], At[la,lb], Ats[la,lb], Rt[la,lb], Rphi[la,lb], T[la,lb]}];
+AssertSymmetricIncreasing [G[ua,lb,lc], lb, lc];
+AssertSymmetricIncreasing [Gtl[la,lb,lc], lb, lc];
+AssertSymmetricIncreasing [Gt[ua,lb,lc], lb, lc];
+AssertSymmetricIncreasing [gK[la,lb,lc], la, lb];
+Map [AssertSymmetricIncreasing,
+     {gu[ua,ub], gtu[ua,ub], Atu[ua,ub]}];
+AssertSymmetricIncreasing [dgtu[ua,ub,lc], ua, ub];
+AssertSymmetricIncreasing [ddgtu[ua,ub,lc,ld], ua, ub];
+AssertSymmetricIncreasing [ddgtu[ua,ub,lc,ld], lc, ld];
+
+DefineConnection [CD, PD, G];
+DefineConnection [CDt, PD, Gt];
+
+(* Use the CartGrid3D variable names *)
+x1=x; x2=y; x3=z;
+
+(* Use the ADMBase variable names *)
+admg11=gxx; admg12=gxy; admg22=gyy; admg13=gxz; admg23=gyz; admg33=gzz;
+admK11=kxx; admK12=kxy; admK22=kyy; admK13=kxz; admK23=kyz; admK33=kzz;
+admalpha=alp;
+admdtalpha=dtalp;
+admbeta1=betax; admbeta2=betay; admbeta3=betaz;
+admdtbeta1=dtbetax; admdtbeta2=dtbetay; admdtbeta3=dtbetaz;
+
+(* Use the TmunuBase variable names *)
+T00=eTtt;
+T01=eTtx; T02=eTty; T03=eTtz;
+T11=eTxx; T12=eTxy; T22=eTyy; T13=eTxz; T23=eTyz; T33=eTzz;
+
+
+
+(******************************************************************************)
+(* Expressions *)
+(******************************************************************************)
+
+(* enum constants for conformalMethod; these must be consistent
+   with the definition of the Cactus parameter conformalMethod *)
+CMphi = 0;
+CMW   = 1;
+
+detgExpr  = Det [MatrixOfComponents [g [la,lb]]];
+ddetgExpr[la_] =
+  Sum [D[Det[MatrixOfComponents[g[la, lb]]], X] PD[X, la],
+       {X, Union[Flatten[MatrixOfComponents[g[la, lb]]]]}];
+
+detgtExpr = Det [MatrixOfComponents [gt[la,lb]]];
+ddetgtExpr[la_] =
+  Sum [D[Det[MatrixOfComponents[gt[la, lb]]], X] PD[X, la],
+       {X, Union[Flatten[MatrixOfComponents[gt[la, lb]]]]}];
+
+etaExpr = SpatialBetaDriverRadius / Max [r, SpatialBetaDriverRadius];
+thetaExpr = Min [Exp [1 - r / SpatialShiftGammaCoeffRadius], 1];
+
+
+
+(******************************************************************************)
+(* Groups *)
+(******************************************************************************)
+
+evolvedGroups =
+  {SetGroupName [CreateGroupFromTensor [phi      ], prefix <> "log_confac"],
+   SetGroupName [CreateGroupFromTensor [gt[la,lb]], prefix <> "metric"    ],
+   SetGroupName [CreateGroupFromTensor [Xt[ua]   ], prefix <> "Gamma"     ],
+   SetGroupName [CreateGroupFromTensor [trK      ], prefix <> "trace_curv"],
+   SetGroupName [CreateGroupFromTensor [At[la,lb]], prefix <> "curv"      ],
+   SetGroupName [CreateGroupFromTensor [alpha    ], prefix <> "lapse"     ],
+   SetGroupName [CreateGroupFromTensor [A        ], prefix <> "dtlapse"   ],
+   SetGroupName [CreateGroupFromTensor [beta[ua] ], prefix <> "shift"     ],
+   SetGroupName [CreateGroupFromTensor [B[ua]    ], prefix <> "dtshift"   ],
+   IfCCZ4[SetGroupName[CreateGroupFromTensor[Theta], prefix <> "Theta"]]};
+evaluatedGroups =
+  {SetGroupName [CreateGroupFromTensor [H      ], prefix <> "Ham"],
+   SetGroupName [CreateGroupFromTensor [M[la]  ], prefix <> "mom"],
+   SetGroupName [CreateGroupFromTensor [cS     ], prefix <> "cons_detg"],
+   SetGroupName [CreateGroupFromTensor [cXt[ua]], prefix <> "cons_Gamma"],
+   SetGroupName [CreateGroupFromTensor [cA     ], prefix <> "cons_traceA"]};
+
+declaredGroups = Join [evolvedGroups, evaluatedGroups];
+declaredGroupNames = Map [First, declaredGroups];
+
+
+
+extraGroups =
+  {{"Grid::coordinates", {x, y, z, r}},
+   {"ADMBase::metric",  {gxx, gxy, gxz, gyy, gyz, gzz}},
+   {"ADMBase::curv",    {kxx, kxy, kxz, kyy, kyz, kzz}},
+   {"ADMBase::lapse",   {alp}},
+   {"ADMBase::dtlapse", {dtalp}},
+   {"ADMBase::shift",   {betax, betay, betaz}},
+   {"ADMBase::dtshift", {dtbetax, dtbetay, dtbetaz}},
+   {"TmunuBase::stress_energy_scalar", {eTtt}},
+   {"TmunuBase::stress_energy_vector", {eTtx, eTty, eTtz}},
+   {"TmunuBase::stress_energy_tensor", {eTxx, eTxy, eTxz, eTyy, eTyz, eTzz}}
+};
+
+groups = Join [declaredGroups, extraGroups];
+
+
+
+(******************************************************************************)
+(* Initial data *)
+(******************************************************************************)
+
+initialCalc =
+{
+  Name -> thorn <> "_Minkowski",
+  Schedule -> {"IN ADMBase_InitialData"},
+  ConditionalOnKeyword -> {"my_initial_data", "Minkowski"},
+  Equations -> 
+  {
+    phi       -> IfThen[conformalMethod==CMW, 1, 0],
+    gt[la,lb] -> KD[la,lb],
+    trK       -> 0,
+    At[la,lb] -> 0,
+    Xt[ua]    -> 0,
+    alpha     -> 1,
+    A         -> 0,
+    beta[ua]  -> 0,
+    B[ua]     -> 0,
+    IfCCZ4[Theta -> 0]
+  }
+};
+
+
+
+(******************************************************************************)
+(* Split a calculation *)
+(******************************************************************************)
+
+PartialCalculation[calc_, suffix_, updates_, evolVars_] :=
+Module[
+  {name, calc1, replaces, calc2, vars, patterns, eqs, calc3},
+  (* Add suffix to name *)
+  name     = lookup[calc, Name] <> suffix;
+  calc1    = mapReplace[calc, Name, name];
+  (* Replace some entries in the calculation *)
+  (* replaces = Map[Function[rule, mapReplace[#, rule[[1]], rule[[2]]]&], updates]; *)
+  replaces = updates //. (lhs_ -> rhs_) -> (mapReplace[#, lhs, rhs]&);
+  calc2 = Apply[Composition, replaces][calc1];
+  (* Remove unnecessary equations *)
+  vars     = Join[evolVars, lookup[calc2, Shorthands]];
+  patterns = Replace[vars, {    Tensor[n_,__]  ->     Tensor[n,__] ,
+                            dot[Tensor[n_,__]] -> dot[Tensor[n,__]]}, 1];
+  eqs      = FilterRules[lookup[calc, Equations], patterns];
+  calc3    = mapReplace[calc2, Equations, eqs];
+  calc3
+];
+
+
+
+(******************************************************************************)
+(* Convert from ADMBase *)
+(******************************************************************************)
+
+convertFromADMBaseCalc =
+{
+  Name -> thorn <> "_convertFromADMBase",
+  Schedule -> {"AT initial AFTER ADMBase_PostInitial"},
+  ConditionalOnKeyword -> {"my_initial_data", "ADMBase"},
+  Shorthands -> {g[la,lb], detg, gu[ua,ub], em4phi},
+  Equations -> 
+  {
+    g[la,lb]  -> admg[la,lb],
+    detg      -> detgExpr,
+    gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]],
+    
+    phi       -> IfThen[conformalMethod==CMW, detg^(-1/6), Log[detg]/12],
+    em4phi    -> IfThen[conformalMethod==CMW, phi^2, Exp[-4 phi]],
+    gt[la,lb] -> em4phi g[la,lb],
+    
+    trK       -> gu[ua,ub] admK[la,lb],
+    At[la,lb] -> em4phi (admK[la,lb] - (1/3) g[la,lb] trK),
+    
+    alpha     -> admalpha,
+    
+    beta[ua]  -> admbeta[ua],
+
+    IfCCZ4[Theta -> 0]
+  }
+};
+
+convertFromADMBaseGammaCalc =
+{
+  Name -> thorn <> "_convertFromADMBaseGamma",
+  Schedule -> {"AT initial AFTER " <> thorn <> "_convertFromADMBase"},
+  ConditionalOnKeyword -> {"my_initial_data", "ADMBase"},
+  (*
+  Where -> InteriorNoSync,
+  *)
+  (* Do not synchronise right after this routine; instead, synchronise
+     after extrapolating *)
+  Where -> Interior,
+  (* Synchronise after this routine, so that the refinement boundaries
+     are set correctly before extrapolating.  (We will need to
+     synchronise again after extrapolating because extrapolation does
+     not fill ghost zones, but this is irrelevant here.)  *)
+  Shorthands -> {dir[ua],
+                 detgt, gtu[ua,ub], Gt[ua,lb,lc], theta},
+  Equations -> 
+  {
+    dir[ua] -> Sign[beta[ua]],
+    
+    detgt        -> 1 (* detgtExpr *),
+    gtu[ua,ub]   -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
+    Gt[ua,lb,lc] -> 1/2 gtu[ua,ud]
+                    (PD[gt[lb,ld],lc] + PD[gt[lc,ld],lb] - PD[gt[lb,lc],ld]),
+    Xt[ua] -> gtu[ub,uc] Gt[ua,lb,lc],
+    
+(*
+    A -> - admdtalpha / (harmonicF alpha^harmonicN) (LapseAdvectionCoeff - 1),
+*)
+    (* If LapseACoeff=0, then A is not evolved, in the sense that it
+       does not influence the time evolution of other variables.  *)
+    A -> IfThen [LapseACoeff != 0,
+                 1 / (- harmonicF alpha^harmonicN)
+                 (+ admdtalpha
+                  - LapseAdvectionCoeff Upwind[beta[ua], alpha, la]),
+                 0],
+    
+    theta -> thetaExpr,
+    
+    (* If ShiftBCoeff=0 or theta ShiftGammaCoeff=0, then B^i is not
+       evolved, in the sense that it does not influence the time
+       evolution of other variables.  *)
+    B[ua] -> IfThen [ShiftGammaCoeff ShiftBCoeff != 0,
+                     1 / (theta ShiftGammaCoeff)
+                     (+ admdtbeta[ua]
+                      - ShiftAdvectionCoeff Upwind[beta[ub], beta[ua], lb]),
+                     0]
+  }
+};
+
+(* Initialise the Gamma variables to 0.  This is necessary with
+   multipatch because convertFromADMBaseGamma does not perform the
+   conversion in the boundary points, and the order in which symmetry
+   (interpatch) and outer boundary conditions is applied means that
+   points which are both interpatch and symmetry points are never
+   initialised. *)
+initGammaCalc =
+{
+  Name -> thorn <> "_InitGamma",
+  Schedule -> {"AT initial BEFORE " <> thorn <> "_convertFromADMBaseGamma"},
+  ConditionalOnKeyword -> {"my_initial_data", "ADMBase"},
+  Where -> Everywhere,
+  Equations -> 
+  {
+    Xt[ua] -> 0,
+    A      -> 0,
+    B[ua]  -> 0
+  }
+};
+
+
+
+(******************************************************************************)
+(* Convert to ADMBase *)
+(******************************************************************************)
+
+convertToADMBaseCalc =
+{
+  Name -> thorn <> "_convertToADMBase",
+  Schedule -> {"IN " <> thorn <> "_convertToADMBaseGroup"},
+  Where -> Everywhere,
+  Shorthands -> {e4phi},
+  Equations -> 
+  {
+    e4phi       -> IfThen[conformalMethod==CMW, 1/phi^2, Exp[4 phi]],
+    admg[la,lb] -> e4phi gt[la,lb],
+    admK[la,lb] -> e4phi At[la,lb] + (1/3) admg[la,lb] trK,
+    admalpha    -> alpha,
+    admbeta[ua] -> beta[ua]
+  }
+};
+
+convertToADMBaseDtLapseShiftCalc =
+{
+  Name -> thorn <> "_convertToADMBaseDtLapseShift",
+  Schedule -> {"IN " <> thorn <> "_convertToADMBaseGroup"},
+  ConditionalOnKeyword -> {"dt_lapse_shift_method", "correct"},
+  Where -> Interior,
+  Shorthands -> {dir[ua], detgt, gtu[ua,ub], eta, theta, em4phi, Ddetgt[la]},
+  Equations -> 
+  {
+    dir[ua] -> Sign[beta[ua]],
+    
+    detgt -> 1 (* detgtExpr *),
+    (* This leads to simpler code... *)
+    gtu[ua,ub]   -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
+    em4phi       -> IfThen[conformalMethod==CMW, phi^2, Exp[-4 phi]],
+    
+    eta -> etaExpr,
+    theta -> thetaExpr,
+
+    (* Ddetgt should be zero analytically, but we're not assuming it here. Change commenting to assume it.*)
+    Ddetgt[la] -> gtu[uk,ul] PD[gt[lk,ll],la],
+    (*Ddetgt[la] -> 0,*)
+ 
+    (* see RHS *)
+(*
+    admdtalpha -> - harmonicF alpha^harmonicN
+                    ((1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK)
+                  + LapseAdvectionCoeff beta[ua] PDu[alpha,la],
+*)
+    admdtalpha -> - harmonicF alpha^harmonicN
+                    (+ LapseACoeff       A
+                     + ((1 - LapseACoeff)
+                        (trK - IfCCZ4[2 Theta, 0])))
+                  + LapseAdvectionCoeff Upwind[beta[ua], alpha, la],
+    admdtbeta[ua] -> IfThen[harmonicShift,
+                            - 1/2 gtu[ua,uj] em4phi alpha
+                              (- 2 alpha IfThen[conformalMethod==CMW,1/phi,-2] PD[phi,lj]
+                               + 2 PD[alpha,lj]
+                               + alpha (Ddetgt[lj] - 2 gtu[uk,ul] 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]
+  }
+};
+
+convertToADMBaseDtLapseShiftBoundaryCalc =
+{
+  Name -> thorn <> "_convertToADMBaseDtLapseShiftBoundary",
+  Schedule -> {"IN " <> thorn <> "_convertToADMBaseGroup"},
+  ConditionalOnKeyword -> {"dt_lapse_shift_method", "correct"},
+  Where -> BoundaryWithGhosts,
+  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,
+    
+    (* see RHS, but omit derivatives near the boundary *)
+(*
+    admdtalpha -> - harmonicF alpha^harmonicN
+                    ((1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK),
+*)
+    admdtalpha -> - harmonicF alpha^harmonicN
+                    (+ LapseACoeff       A
+                     + ((1 - LapseACoeff)
+                        (trK - IfCCZ4[2 Theta, 0]))),
+    admdtbeta[ua] -> IfThen[harmonicShift,
+                            0,
+                            (* else *)
+                            + theta ShiftGammaCoeff
+                              (+ ShiftBCoeff B[ua]
+                               + (1 - ShiftBCoeff)
+                                 (Xt[ua] - eta BetaDriver beta[ua]))]
+  }
+};
+
+convertToADMBaseFakeDtLapseShiftCalc =
+{
+  Name -> thorn <> "_convertToADMBaseFakeDtLapseShift",
+  Schedule -> {"IN " <> thorn <> "_convertToADMBaseGroup"},
+  ConditionalOnKeyword -> {"dt_lapse_shift_method", "noLapseShiftAdvection"},
+  Where -> Everywhere,
+  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,
+    
+    (* see RHS, but omit derivatives everywhere (which is wrong, but
+       faster, since it does not require synchronisation or boundary
+       conditions) *)
+(*
+    admdtalpha -> - harmonicF alpha^harmonicN
+                    ((1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK),
+*)
+    admdtalpha -> - harmonicF alpha^harmonicN
+                    (+ LapseACoeff       A
+                     + ((1 - LapseACoeff)
+                        (trK - IfCCZ4[2 Theta, 0]))),
+    admdtbeta[ua] -> IfThen[harmonicShift,
+                            0,
+                            (* else *)
+                            + theta ShiftGammaCoeff
+                              (+ ShiftBCoeff B[ua]
+                               + (1 - ShiftBCoeff)
+                                 (Xt[ua] - eta BetaDriver beta[ua]))]
+  }
+};
+
+(******************************************************************************)
+(* Evolution equations *)
+(******************************************************************************)
+
+evolCalc =
+{
+  Name -> thorn <> "_RHS",
+  Schedule -> {"IN " <> thorn <> "_evolCalcGroup"},
+  (*
+  Where -> Interior,
+  *)
+  (* Synchronise the RHS grid functions after this routine, so that
+     the refinement boundaries are set correctly before applying the
+     radiative boundary conditions.  *)
+  Where -> InteriorNoSync,
+  Shorthands -> {dir[ua],
+                 detgt, gtu[ua,ub],
+                 Gt[ua,lb,lc], Gtl[la,lb,lc], Gtlu[la,lb,uc], Xtn[ua],
+                 Rt[la,lb], Rphi[la,lb], R[la,lb],
+                 Atm[ua,lb], Atu[ua,ub],
+                 e4phi, em4phi, cdphi[la], cdphi2[la,lb], g[la,lb], detg,
+                 gu[ua,ub], Ats[la,lb], trAts, eta, theta,
+                 rho, S[la], trS, fac1, fac2, dottrK, dotXt[ua],
+                 epsdiss[ua], IfCCZ4[Z[ua]], IfCCZ4[dotTheta], Ddetgt[la]},
+  Equations -> 
+  {
+    dir[ua] -> Sign[beta[ua]],
+    
+    detgt -> 1 (* detgtExpr *),
+    
+    (* This leads to simpler code... *)
+    gtu[ua,ub]   -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
+    Gtl[la,lb,lc]  -> 1/2
+                      (PD[gt[lb,la],lc] + PD[gt[lc,la],lb] - PD[gt[lb,lc],la]),
+    Gtlu[la,lb,uc] -> gtu[uc,ud] Gtl[la,lb,ld],
+    Gt[ua,lb,lc]   -> gtu[ua,ud] Gtl[ld,lb,lc],
+ 
+    (* The conformal connection functions calculated from the conformal metric,
+       used instead of Xt where no derivatives of Xt are taken *)
+    Xtn[ui] -> gtu[uj,uk] Gt[ui,lj,lk],
+
+    e4phi       -> IfThen[conformalMethod==CMW, 1/phi^2, Exp[4 phi]],
+    em4phi      -> 1 / e4phi,
+    g[la,lb]    -> e4phi gt[la,lb],
+    detg        -> detgExpr,
+    gu[ua,ub]   -> em4phi gtu[ua,ub],
+
+    (* The Z quantities *)
+    (* gr-qc:1106.2254 (2011), eqn. (23) *)
+    IfCCZ4[
+      Z[ud] -> (1/2) gu[ua,ud] (- PD[gt[la,lb],lc] gtu[ub,uc] + gt[la,lc] Xt[uc])
+    ],
+
+    (* PRD 62, 044034 (2000), eqn. (18) *)
+    (* Adding Z term by changing Xtn to Xt *)
+    Rt[li,lj] -> - (1/2) gtu[ul,um] PD[gt[li,lj],ll,lm]
+                 + (1/2) gt[lk,li] PD[Xt[uk],lj]
+                 + (1/2) gt[lk,lj] PD[Xt[uk],li]
+                 + (1/2) Xtn[uk] Gtl[li,lj,lk]
+                 + (1/2) Xtn[uk] Gtl[lj,li,lk]
+                 + (+ Gt[uk,li,ll] Gtlu[lj,lk,ul]
+                    + Gt[uk,lj,ll] Gtlu[li,lk,ul]
+                    + Gt[uk,li,ll] Gtlu[lk,lj,ul]),
+
+    fac1 -> IfThen[conformalMethod==CMW, -1/(2 phi), 1],
+    cdphi[la] -> fac1 CDt[phi,la],
+    fac2 -> IfThen[conformalMethod==CMW, 1/(2 phi^2), 0],
+    cdphi2[la,lb] -> fac1 CDt[phi,la,lb] + fac2 CDt[phi,la] CDt[phi,lb],
+
+    (* PRD 62, 044034 (2000), eqn. (15) *)
+    Rphi[li,lj] -> - 2 cdphi2[lj,li]
+                   - 2 gt[li,lj] gtu[ul,un] cdphi2[ll,ln]
+                   + 4 cdphi[li] cdphi[lj]
+                   - 4 gt[li,lj] gtu[ul,un] cdphi[ln] cdphi[ll],
+    
+    Atm[ua,lb] -> gtu[ua,uc] At[lc,lb],
+    Atu[ua,ub] -> Atm[ua,lc] gtu[ub,uc],
+    
+    R[la,lb] -> Rt[la,lb] + Rphi[la,lb],
+    IfCCZ4[
+      R[la,lb] -> R[la,lb] + (2/phi) (+ g[la,lc] Z[uc] PD[phi,lb]
+        + g[lb,lc] Z[uc] PD[phi,la] - g[la,lb] Z[uc] PD[phi,lc])
+        + e4phi Z[uc] PD[gt[la,lb],lc]
+    ],
+    
+    (* Matter terms *)
+    
+    (* rho = n^a n^b T_ab *)
+    rho -> addMatter
+           (1/alpha^2 (T00 - 2 beta[ui] T0[li] + beta[ui] beta[uj] T[li,lj])),
+    
+    (* S_i = -p^a_i n^b T_ab, where p^a_i = delta^a_i + n^a n_i *)
+    S[li] -> addMatter (-1/alpha (T0[li] - beta[uj] T[li,lj])),
+    
+    (* trS = gamma^ij T_ij  *)
+    trS -> addMatter (em4phi gtu[ui,uj] T[li,lj]),
+    
+    (* RHS terms *)
+    
+    (* PRD 62, 044034 (2000), eqn. (10) *)
+    (* PRD 67 084023 (2003), eqn. (16) and (23) *)
+    dot[phi]       -> IfThen[conformalMethod==CMW, 1/3 phi, -1/6]
+                      (alpha trK - PD[beta[ua],la]),
+    
+    (* PRD 62, 044034 (2000), eqn. (9) *)
+    (* gr-qc:1106.2254 (2011), eqn. (14) *)
+    (* removing trA from Aij ensures that detg = 1 *)
+    dot[gt[la,lb]] -> - 2 alpha (At[la,lb] - IfCCZ4[(1/3) At[lc,ld] gtu[uc,ud] gt[la,lb], 0])
+                      + gt[la,lc] PD[beta[uc],lb] + gt[lb,lc] PD[beta[uc],la]
+                      - (2/3) gt[la,lb] PD[beta[uc],lc],
+    (* PRD 62, 044034 (2000), eqn. (20) *)
+    (* PRD 67 084023 (2003), eqn (26) *)
+    (* gr-qc:1106.2254 (2011), eqn. (19) *)
+    (* Adding Z terms by changing Xtn to Xt,
+       also adding extra Z and Theta terms *)
+    dotXt[ui]      -> - 2 Atu[ui,uj] PD[alpha,lj]
+                      + 2 alpha (+ Gt[ui,lj,lk] Atu[uk,uj]
+                                 - (2/3) gtu[ui,uj] PD[trK,lj]
+                                 + 6 Atu[ui,uj] cdphi[lj])
+                      + gtu[uj,ul] PD[beta[ui],lj,ll]
+                      + (1/3) gtu[ui,uj] PD[beta[ul],lj,ll]
+                      - Xtn[uj] PD[beta[ui],lj] 
+                      + (2/3) Xtn[ui] PD[beta[uj],lj]
+                      + IfCCZ4[
+                               + GammaShift 2 e4phi (- Z[uj] PD[beta[ui],lj]
+                                                 + (2/3) Z[ui] PD[beta[uj],lj])
+                               - (4/3) alpha e4phi Z[ui] trK
+                               + 2 gtu[ui,uj] (+ alpha PD[Theta,lj]
+                                               - Theta PD[alpha,lj])
+                               - 2 alpha e4phi dampk1 Z[ui],
+                        0]
+    (* Equation (4.28) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *)
+                      + addMatter (- 16 Pi alpha gtu[ui,uj] S[lj]),
+    dot[Xt[ui]]    -> dotXt[ui],
+
+    (* gr-qc:1106.2254 (2011), eqn. (18) *)
+    IfCCZ4[
+      dotTheta ->
+        - PD[alpha,la] Z[ua] - dampk1 (2 + dampk2) alpha Theta
+        + (1/2) alpha (gu[ua,ub] R[la,lb] - Atm[ua,lb] Atm[ub,la] + (2/3) trK^2 - 2 trK Theta) 
+        + addMatter (- 8 Pi alpha rho)
+    ],
+
+    IfCCZ4[
+      dot[Theta] -> dotTheta
+    ],
+
+    (* PRD 62, 044034 (2000), eqn. (11) *)
+    (* gr-qc:1106.2254 (2011), eqn. (17) *)
+    (* Adding the RHS of Theta to K, because K_Z4 = K_BSSN + 2 Theta *)
+    (* Also adding the Z term, as it has to cancel with the one in Theta *)
+    dottrK         -> - em4phi ( gtu[ua,ub] ( PD[alpha,la,lb]
+                                + 2 cdphi[la] PD[alpha,lb] )
+                                - Xtn[ua] PD[alpha,la] )
+                      + alpha (Atm[ua,lb] Atm[ub,la] + (1/3) trK^2)
+                      + IfCCZ4[
+                               + 2 dotTheta + 2 PD[alpha,la] Z[ua]
+                               + dampk1 (1 - dampk2) alpha Theta,
+                               0]
+    (* Equation (4.21) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *)
+                      + addMatter (4 Pi alpha (rho + trS)),
+    dot[trK]       -> dottrK,
+
+    (* PRD 62, 044034 (2000), eqn. (12) *)
+    (* TODO: Should we use the Hamiltonian constraint to make Rij tracefree? *)
+    (* gr-qc:1106.2254 (2011), eqn. (15) *)
+    (* Adding Z terms in the Ricci and Theta terms *)
+    Ats[la,lb]     -> - CDt[alpha,la,lb] +
+                      + 2 (PD[alpha,la] cdphi[lb] + PD[alpha,lb] cdphi[la] )
+                      + alpha R[la,lb],
+    trAts          -> gu[ua,ub] Ats[la,lb],
+    dot[At[la,lb]] -> + em4phi (+ Ats[la,lb] - (1/3) g[la,lb] trAts )
+                      + alpha (+ ((trK - IfCCZ4[2 Theta, 0])
+                                  At[la,lb])
+                               - 2 At[la,lc] Atm[uc,lb])
+                      + At[la,lc] PD[beta[uc],lb] + At[lb,lc] PD[beta[uc],la]
+                      - (2/3) At[la,lb] PD[beta[uc],lc]
+    (* Equation (4.23) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *)
+                      + addMatter (- em4phi alpha 8 Pi
+                                     (T[la,lb] - (1/3) g[la,lb] trS)),
+    
+    (* dot[alpha] -> - harmonicF alpha^harmonicN trK, *)
+    (* dot[alpha] -> - harmonicF alpha^harmonicN A + Lie[alpha, beta], *)
+(*
+    dot[alpha] -> - harmonicF alpha^harmonicN (
+                    (1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK)
+                  + LapseAdvectionCoeff beta[ua] PDu[alpha,la],
+
+
+    dot[A]     -> (1 - LapseAdvectionCoeff) (dottrK - AlphaDriver A),
+
+*)
+    dot[alpha] -> - harmonicF alpha^harmonicN
+                    (+ LapseACoeff       A
+                     + ((1 - LapseACoeff)
+                        (+ trK - IfCCZ4[2 Theta, 0]
+                         + AlphaDriver (alpha - 1)))),
+
+    dot[A]     -> + (LapseACoeff
+                     (+ dottrK - IfCCZ4[2 dotTheta, 0]
+                      - AlphaDriver A)),
+    
+    eta -> etaExpr,
+    theta -> thetaExpr,
+
+    (* Ddetgt should be zero analytically, but we're not assuming it here. Change commenting to assume it.*)
+    Ddetgt[la] -> gtu[uk,ul] PD[gt[lk,ll],la],
+    (*Ddetgt[la] -> 0,*)
+
+    (* dot[beta[ua]] -> eta Xt[ua], *)
+    (* dot[beta[ua]] -> ShiftGammaCoeff alpha^ShiftAlphaPower B[ua], *)
+    dot[beta[ua]] -> IfThen[harmonicShift,
+                            - 1/2 gtu[ua,uj] em4phi alpha
+                              (- 2 alpha IfThen[conformalMethod==CMW,1/phi,-2] PD[phi,lj]
+                               + 2 PD[alpha,lj]
+                               + alpha (Ddetgt[lj] - 2 gtu[uk,ul] 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])
+      (* Note that this dotXt[ua] is not yet \partial_t \tilde \Gamma^i, because the
+         advection term has not yet been added.  It is actually
+         \partial_t \tilde \Gamma^i - \beta^j \partial_j \tilde \Gamma^i *)
+  }
+};
+
+advectCalc =
+{
+  Name -> thorn <> "_Advect",
+  Schedule -> {"IN " <> thorn <> "_evolCalcGroup " <>
+               "AFTER (" <> thorn <> "_RHS1 " <> thorn <> "_RHS2)"},
+  (*
+  Where -> Interior,
+  *)
+  (* Synchronise the RHS grid functions after this routine, so that
+     the refinement boundaries are set correctly before applying the
+     radiative boundary conditions.  *)
+  Where -> InteriorNoSync,
+  Shorthands -> {dir[ua]},
+  Equations ->
+  {
+           dir[ua] -> Sign[beta[ua]],
+
+          dot[phi] -> dot[phi] + Upwind[beta[ua], phi, la],
+
+    dot[gt[la,lb]] -> dot[gt[la,lb]] + Upwind[beta[uc], gt[la,lb], lc],
+
+       dot[Xt[ui]] -> dot[Xt[ui]] + Upwind[beta[uj], Xt[ui], lj],
+
+    IfCCZ4[
+        dot[Theta] -> dot[Theta] + Upwind[beta[ua], Theta, la]
+    ],
+
+          dot[trK] -> dot[trK] + Upwind[beta[ua], trK, la],
+
+    dot[At[la,lb]] -> dot[At[la,lb]] + Upwind[beta[uc], At[la,lb], lc],
+
+        dot[alpha] -> dot[alpha]
+                      + LapseAdvectionCoeff Upwind[beta[ua], alpha, la],
+
+            dot[A] -> dot[A]
+                      + LapseACoeff (
+                        + LapseAdvectionCoeff       Upwind[beta[ua], A, la]
+                        + (1 - LapseAdvectionCoeff) Upwind[beta[ua], trK, la]),
+
+     dot[beta[ua]] -> dot[beta[ua]]
+                      + ShiftAdvectionCoeff Upwind[beta[ub], beta[ua], lb],
+
+        dot[B[ua]] -> dot[B[ua]]
+                      + ShiftBCoeff (
+                        + ShiftAdvectionCoeff Upwind[beta[ub], B[ua], lb]
+                        + ((1 - ShiftAdvectionCoeff)
+                           Upwind[beta[ub], Xt[ua], lb]))
+                      (* Note that the advection term \beta^j \partial_j \tilde \Gamma^i is not
+                         subtracted here when ShiftAdvectionCoefficient == 1 because it was
+                         implicitly subtracted before (see comment in previous calculation of
+                         dot[B[ua]]. *)
+  }
+};
+
+evolCalc1 = PartialCalculation[evolCalc, "1",
+  {
+    ConditionalOnKeyword -> {"RHS_calculation", "split"}
+  },
+  {
+    dot[phi],
+    dot[gt[la,lb]],
+    dot[Xt[ui]],
+    dot[trK],
+    dot[alpha],
+    dot[A],
+    dot[beta[ua]],
+    dot[B[ua]],
+    IfCCZ4[dot[Theta]]
+  }];
+
+evolCalc2 = PartialCalculation[evolCalc, "2",
+  {
+    ConditionalOnKeyword -> {"RHS_calculation", "split"}
+  },
+  {
+    dot[At[la,lb]]
+  }];
+
+dissCalc =
+{
+  Name -> thorn <> "_Dissipation",
+  Schedule -> {"IN " <> thorn <> "_evolCalcGroup " <>
+               "AFTER (" <> thorn <> "_RHS1 " <> thorn <> "_RHS2)"},
+  ConditionalOnKeyword -> {"apply_dissipation", "always"},
+  Where -> InteriorNoSync,
+  Shorthands -> {epsdiss[ua]},
+  Equations ->
+  {
+    epsdiss[ua] -> EpsDiss,
+    Sequence@@Table[
+      dot[var]       -> dot[var] + epsdiss[ux] PDdiss[var,lx],
+      {var, {phi, gt[la,lb], Xt[ui], IfCCZ4[Theta], trK, At[la,lb],
+             alpha, A, beta[ua], B[ua]}}]
+  }
+};
+
+dissCalcs =
+Table[
+{
+  Name -> thorn <> "_Dissipation_" <> ToString[var /. {Tensor[n_,__] -> n}],
+  Schedule -> {"IN " <> thorn <> "_evolCalcGroup " <>
+               "AFTER (" <> thorn <> "_RHS1 " <> thorn <> "_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], IfCCZ4[Theta], trK, At[la,lb],
+         alpha, A, beta[ua], B[ua]}}
+];
+
+RHSStaticBoundaryCalc =
+{
+  Name -> thorn <> "_RHSStaticBoundary",
+  Schedule -> {"IN MoL_CalcRHS"},
+  ConditionalOnKeyword -> {"my_rhs_boundary_condition", "static"},
+  Where -> Boundary,
+  Equations -> 
+  {
+    dot[phi]       -> 0,
+    dot[gt[la,lb]] -> 0,
+    dot[trK]       -> 0,
+    dot[At[la,lb]] -> 0,
+    dot[Xt[ua]]    -> 0,
+    dot[alpha]     -> 0,
+    dot[A]         -> 0,
+    dot[beta[ua]]  -> 0,
+    dot[B[ua]]     -> 0,
+    IfCCZ4[dot[Theta] -> 0]
+  }
+};
+
+(* Initialise the RHS variables in analysis in case they are going to
+   be output - the noninterior points cannot be filled, so we define
+   them to be zero *)
+initRHSCalc =
+{
+  Name -> thorn <> "_InitRHS",
+  Schedule -> {"AT analysis BEFORE " <> thorn <> "_evolCalcGroup"},
+  Where -> Everywhere,
+  Equations -> 
+  {
+    dot[phi]       -> 0,
+    dot[gt[la,lb]] -> 0,
+    dot[trK]       -> 0,
+    dot[At[la,lb]] -> 0,
+    dot[Xt[ua]]    -> 0,
+    dot[alpha]     -> 0,
+    dot[A]         -> 0,
+    dot[beta[ua]]  -> 0,
+    dot[B[ua]]     -> 0,
+    IfCCZ4[dot[Theta] -> 0]
+  }
+};
+
+RHSRadiativeBoundaryCalc =
+{
+  Name -> thorn <> "_RHSRadiativeBoundary",
+  Schedule -> {"IN MoL_CalcRHS"},
+  ConditionalOnKeyword -> {"my_rhs_boundary_condition", "radiative"},
+  Where -> Boundary,
+  Shorthands -> {dir[ua],
+                 detgt, gtu[ua,ub], em4phi, gu[ua,ub],
+                 nn[la], nu[ua], nlen, nlen2, su[ua],
+                 vg},
+  Equations -> 
+  {
+    dir[ua] -> Sign[normal[ua]],
+    
+    detgt      -> 1 (* detgtExpr *),
+    gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
+    em4phi     -> IfThen[conformalMethod==CMW, phi^2, Exp[-4 phi]],
+    gu[ua,ub]  -> em4phi gtu[ua,ub],
+    
+    nn[la] -> Euc[la,lb] normal[ub],
+    nu[ua] -> gu[ua,ub] nn[lb],
+    nlen2  -> nu[ua] nn[la],
+    nlen   -> Sqrt[nlen2],
+    su[ua] -> nu[ua] / nlen,
+    
+    vg -> Sqrt[harmonicF],
+    
+    dot[phi]       -> - vg su[uc] PDo[phi      ,lc],
+    dot[gt[la,lb]] -> -    su[uc] PDo[gt[la,lb],lc],
+    dot[trK]       -> - vg su[uc] PDo[trK      ,lc],
+    dot[At[la,lb]] -> -    su[uc] PDo[At[la,lb],lc],
+    dot[Xt[ua]]    -> -    su[uc] PDo[Xt[ua]   ,lc],
+    dot[alpha]     -> - vg su[uc] PDo[alpha    ,lc],
+    dot[A]         -> - vg su[uc] PDo[A        ,lc],
+    dot[beta[ua]]  -> -    su[uc] PDo[beta[ua] ,lc],
+    dot[B[ua]]     -> -    su[uc] PDo[B[ua]    ,lc],
+    IfCCZ4[
+      dot[Theta]   -> - vg su[uc] PDo[Theta    ,lc]
+    ]
+  }
+};
+
+enforceCalc =
+{
+  Name -> thorn <> "_enforce",
+  Schedule -> {"IN MoL_PostStepModify"},
+  Shorthands -> {detgt, gtu[ua,ub], trAt},
+  Equations -> 
+  {
+    (* 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,
+       which is not always the case (since we don't enforce it). On
+       the other hand, this may not be so important... *)
+    detgt -> 1 (* detgtExpr *),
+    gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
+    
+    trAt -> gtu[ua,ub] At[la,lb],
+    
+    At[la,lb] -> At[la,lb] - (1/3) gt[la,lb] trAt,
+    
+    alpha -> Max[alpha, MinimumLapse]
+  }
+};
+
+(******************************************************************************)
+(* Boundary conditions *)
+(******************************************************************************)
+
+boundaryCalc =
+{
+  Name -> thorn <> "_boundary",
+  Schedule -> {"IN MoL_PostStep"},
+  ConditionalOnKeyword -> {"my_boundary_condition", "Minkowski"},
+  Where -> BoundaryWithGhosts,
+  Equations -> 
+  {
+    phi       -> IfThen[conformalMethod==CMW, 1, 0],
+    gt[la,lb] -> KD[la,lb],
+    trK       -> 0,
+    At[la,lb] -> 0,
+    Xt[ua]    -> 0,
+    alpha     -> 1,
+    A         -> 0,
+    beta[ua]  -> 0,
+    B[ua]     -> 0,
+    IfCCZ4[Theta -> 0]
+  }
+};
+
+(******************************************************************************)
+(* Constraint equations *)
+(******************************************************************************)
+
+constraintsCalc =
+{
+  Name -> thorn <> "_constraints",
+  Schedule -> Automatic,
+  After -> "MoL_PostStep",
+  Where -> Interior,
+  Shorthands -> {detgt, ddetgt[la], gtu[ua,ub], Z[ua],
+                 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[ua,lb],
+                 gK[la,lb,lc], cdphi[la], cdphi2[la,lb],
+                 rho, S[la], fac1, fac2},
+  Equations -> 
+  {
+    detgt        -> 1 (* detgtExpr *),
+    ddetgt[la]   -> 0 (* ddetgtExpr[la] *),
+    
+    (* This leads to simpler code... *)
+    gtu[ua,ub]   -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
+    Gtl[la,lb,lc]  -> 1/2
+                      (PD[gt[lb,la],lc] + PD[gt[lc,la],lb] - PD[gt[lb,lc],la]),
+    Gtlu[la,lb,uc] -> gtu[uc,ud] Gtl[la,lb,ld],
+    Gt[ua,lb,lc]   -> gtu[ua,ud] Gtl[ld,lb,lc],
+    
+    (* The conformal connection functions calculated from the conformal metric,
+       used instead of Xt where no derivatives of Xt are taken *)
+    Xtn[ui] -> gtu[uj,uk] Gt[ui,lj,lk],
+
+    e4phi       -> IfThen[conformalMethod==CMW, 1/phi^2, Exp[4 phi]],
+    em4phi      -> 1 / e4phi,
+    g[la,lb]    -> e4phi gt[la,lb],
+    detg        -> e4phi^3,
+    gu[ua,ub]   -> em4phi gtu[ua,ub],
+
+    (* The Z quantities *)
+    IfCCZ4[
+      Z[ud] -> (1/2) gu[ua,ud] (- PD[gt[la,lb],lc] gtu[ub,uc] + gt[la,lc] Xt[uc])
+    ],
+    
+    (* PRD 62, 044034 (2000), eqn. (18) *)
+    Rt[li,lj] -> - (1/2) gtu[ul,um] PD[gt[li,lj],ll,lm]
+                 + (1/2) gt[lk,li] PD[Xt[uk],lj]
+                 + (1/2) gt[lk,lj] PD[Xt[uk],li]
+                 + (1/2) Xtn[uk] Gtl[li,lj,lk]
+                 + (1/2) Xtn[uk] Gtl[lj,li,lk]
+                 + (+ Gt[uk,li,ll] Gtlu[lj,lk,ul]
+                    + Gt[uk,lj,ll] Gtlu[li,lk,ul]
+                    + Gt[uk,li,ll] Gtlu[lk,lj,ul]),
+
+    (* From the long turducken paper.
+       This expression seems to give the same result as the one from 044034.  *)
+    (* TODO: symmetrise correctly: (ij) = (1/2) [i+j] *)
+(*
+    Rt[li,lj] -> - (1/2) gtu[uk,ul] PD[gt[li,lj],lk,ll]
+                 + gt[lk,li] PD[Xt[uk],lj] + gt[lk,lj] PD[Xt[uk],li]
+                 + gt[li,ln] Gt[un,lj,lk] gtu[um,ua] gtu[uk,ub] PD[gt[la,lb],lm]
+                 + gt[lj,ln] Gt[un,li,lk] gtu[um,ua] gtu[uk,ub] PD[gt[la,lb],lm]
+                 + gtu[ul,us] (+ 2 Gt[uk,ll,li] gt[lj,ln] Gt[un,lk,ls]
+                               + 2 Gt[uk,ll,lj] gt[li,ln] Gt[un,lk,ls]
+                               + Gt[uk,li,ls] gt[lk,ln] Gt[un,ll,lj]),
+*)
+
+    (* Below would be a straightforward calculation,
+       without taking any Gamma^i into account.
+       This expression gives a different answer!  *)
+(*
+    Rt[la,lb] -> + Gt[u1,l2,la] Gt[l1,lb,u2] - Gt[u1,la,lb] Gt[l1,l2,u2]
+                 + 1/2 gtu[u1,u2] (- PD[gt[l1,l2],la,lb] + PD[gt[l1,la],l2,lb]
+                                   - PD[gt[la,lb],l1,l2] + PD[gt[l2,lb],l1,la]),
+*)
+
+    fac1 -> IfThen[conformalMethod==CMW, -1/(2 phi), 1],
+    cdphi[la] -> fac1 CDt[phi,la],
+    fac2 -> IfThen[conformalMethod==CMW, 1/(2 phi^2), 0],
+    cdphi2[la,lb] -> fac1 CDt[phi,la,lb] + fac2 CDt[phi,la] CDt[phi,lb],
+
+    (* PRD 62, 044034 (2000), eqn. (15) *)
+    Rphi[li,lj] -> - 2 cdphi2[lj,li]
+                   - 2 gt[li,lj] gtu[ul,un] cdphi2[ll,ln]
+                   + 4 cdphi[li] cdphi[lj]
+                   - 4 gt[li,lj] gtu[ul,un] cdphi[ln] cdphi[ll],
+    
+    (* ddetg[la] -> PD[e4phi detg,la], *)
+    ddetg[la]   -> e4phi ddetgt[la] + 4 detgt e4phi PD[phi,la],
+    (* TODO: check this equation, maybe simplify it by omitting ddetg *)
+    G[ua,lb,lc] -> Gt[ua,lb,lc]
+                   + 1/(2 detg) (+ KD[ua,lb] ddetg[lc] + KD[ua,lc] ddetg[lb]
+                                 - (1/3) g[lb,lc] gu[ua,ud] ddetg[ld]),
+    
+    R[la,lb] -> + Rt[la,lb] + Rphi[la,lb],
+
+    IfCCZ4[
+      R[la,lb] -> R[la, lb] + (2/phi) (+ g[la,lc] Z[uc] PD[phi,lb]
+        + g[lb,lc] Z[uc] PD[phi,la] - g[la,lb] Z[uc] PD[phi,lc])
+        + e4phi Z[uc] PD[gt[la,lb],lc]
+    ],
+
+    trR      -> gu[ua,ub] R[la,lb],
+    
+    (* K[la,lb] -> e4phi At[la,lb] + (1/3) g[la,lb] trK, *)
+    (* Km[ua,lb] -> gu[ua,uc] K[lc,lb], *)
+    Atm[ua,lb] -> gtu[ua,uc] At[lc,lb],
+    
+    (* Matter terms *)
+    
+    (* rho = n^a n^b T_ab *)
+    rho -> 1/alpha^2 (T00 - 2 beta[ui] T0[li] + beta[ui] beta[uj] T[li,lj]),
+    
+    (* S_i = -p^a_i n^b T_ab, where p^a_i = delta^a_i + n^a n_i *)
+    S[li] -> -1/alpha (T0[li] - beta[uj] T[li,lj]),
+    
+    (* Constraints *)
+    
+    (* H -> trR - Km[ua,lb] Km[ub,la] + trK^2, *)
+    (* PRD 67, 084023 (2003), eqn. (19) *)
+    H -> trR - Atm[ua,lb] Atm[ub,la] + (2/3) trK^2 - addMatter 16 Pi rho,
+    
+    (* gK[la,lb,lc] -> CD[K[la,lb],lc], *)
+(*    gK[la,lb,lc] -> + 4 e4phi PD[phi,lc] At[la,lb] + e4phi CD[At[la,lb],lc]
+                    + (1/3) g[la,lb] PD[trK,lc],
+
+    M[la] -> gu[ub,uc] (gK[lc,la,lb] - gK[lc,lb,la]), *)
+
+    M[li] -> + gtu[uj,uk] (CDt[At[li,lj],lk] + 6 At[li,lj] cdphi[lk])
+             - (2/3) PD[trK,li]
+             - addMatter 8 Pi S[li],
+    (* TODO: use PRD 67, 084023 (2003), eqn. (20) *)
+    
+    (* det gamma-tilde *)
+    cS -> Log[detgt],
+    
+    (* Gamma constraint *)
+    cXt[ua] -> gtu[ub,uc] Gt[ua,lb,lc] - Xt[ua],
+    
+    (* trace A-tilde *)
+    cA -> gtu[ua,ub] At[la,lb]
+  }
+};
+
+constraintsCalc1 = PartialCalculation[constraintsCalc, "1",
+  {},
+  {
+    H
+  }];
+
+constraintsCalc2 = PartialCalculation[constraintsCalc, "2",
+  {},
+  {
+    M[li],
+    cS,
+    cXt[ua],
+    cA
+  }];
+
+(******************************************************************************)
+(* Implementations *)
+(******************************************************************************)
+
+inheritedImplementations =
+  Join[{"ADMBase"},
+       If [addMatter!=0, {"TmunuBase"}, {}]];
+
+(******************************************************************************)
+(* Parameters *)
+(******************************************************************************)
+
+inheritedKeywordParameters = {};
+
+extendedKeywordParameters =
+{
+  {
+    Name -> "ADMBase::evolution_method",
+    AllowedValues -> {thorn}
+  },
+  {
+    Name -> "ADMBase::lapse_evolution_method",
+    AllowedValues -> {thorn}
+  },
+  {
+    Name -> "ADMBase::shift_evolution_method",
+    AllowedValues -> {thorn}
+  },
+  {
+    Name -> "ADMBase::dtlapse_evolution_method",
+    AllowedValues -> {thorn}
+  },
+  {
+    Name -> "ADMBase::dtshift_evolution_method",
+    AllowedValues -> {thorn}
+  }
+};
+
+keywordParameters =
+{
+  {
+    Name -> "my_initial_data",
+    (* Visibility -> "restricted", *)
+    (* Description -> "ddd", *)
+    AllowedValues -> {"ADMBase", "Minkowski"},
+    Default -> "ADMBase"
+  },
+  {
+    Name -> "my_initial_boundary_condition",
+    Visibility -> "restricted",
+    (* Description -> "ddd", *)
+    AllowedValues -> {"none"},
+    Default -> "none"
+  },
+  {
+    Name -> "my_rhs_boundary_condition",
+    Visibility -> "restricted",
+    (* Description -> "ddd", *)
+    AllowedValues -> {"none", "static", "radiative"},
+    Default -> "none"
+  },
+  {
+    Name -> "my_boundary_condition",
+    (* Visibility -> "restricted", *)
+    (* Description -> "ddd", *)
+    AllowedValues -> {"none", "Minkowski"},
+    Default -> "none"
+  },
+  {
+    Name -> "calculate_ADMBase_variables_at",
+    Visibility -> "restricted",
+    (* Description -> "ddd", *)
+    AllowedValues -> {"MoL_PostStep", "CCTK_EVOL", "CCTK_ANALYSIS"},
+    Default -> "MoL_PostStep"
+  },
+  {
+    Name -> "UseSpatialBetaDriver",
+    Visibility -> "restricted",
+    (* Description -> "ddd", *)
+    AllowedValues -> {"no", "yes"},
+    Default -> "no"
+  },
+  {
+    Name -> "dt_lapse_shift_method",
+    Description -> "Treatment of ADMBase dtlapse and dtshift",
+    AllowedValues -> {"correct",
+                      "noLapseShiftAdvection" (* omit lapse and shift advection terms (faster) *)
+                     },
+    Default -> "correct"
+  },
+  {
+    Name -> "apply_dissipation",
+    Description -> "Whether to apply dissipation to the RHSs",
+    AllowedValues -> {"always",
+                      "never" (* yes and no keyword values confuse Cactus, and Kranc
+                                 doesn't support boolean parameters *)
+                     },
+    Default -> "never"
+  }
+
+};
+
+intParameters =
+{
+  {
+    Name -> harmonicN,
+    Description -> "d/dt alpha = - f alpha^n K  (harmonic=2, 1+log=1)",
+    Default -> 2
+  },
+  {
+    Name -> ShiftAlphaPower,
+    Default -> 0
+  },
+  {
+    Name -> conformalMethod,
+    Description -> "Treatment of conformal factor",
+    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
+  }
+};
+
+realParameters =
+{
+  IfCCZ4[{
+    Name -> GammaShift,
+    Description -> "Covariant shift term in Gamma",
+    Default -> 0.5
+  }],   
+  IfCCZ4[{
+    Name -> dampk1,
+    Description -> "CCZ4 damping term 1 for Theta and Z",
+    Default -> 0
+  }],
+  IfCCZ4[{
+    Name -> dampk2,
+    Description -> "CCZ4 damping term 2 for Theta and Z",
+    Default -> 0
+  }],
+  {
+    Name -> LapseACoeff,
+    Description -> "Whether to evolve A in time",
+    Default -> 0
+  },
+  {
+    Name -> harmonicF,
+    Description -> "d/dt alpha = - f alpha^n K   (harmonic=1, 1+log=2)",
+    Default -> 1
+  },
+  {
+    Name -> AlphaDriver,
+    Default -> 0
+  },
+  {
+    Name -> ShiftBCoeff,
+    Description -> "Whether to evolve B^i in time",
+    Default -> 1
+  },
+  {
+    Name -> ShiftGammaCoeff,
+    Default -> 0
+  },
+  {
+    Name -> BetaDriver,
+    Default -> 0
+  },
+  {
+    Name -> LapseAdvectionCoeff,
+    Description -> "Factor in front of the lapse advection terms in 1+log",
+    Default -> 1
+  },
+  {
+    Name -> ShiftAdvectionCoeff,
+    Description -> "Factor in front of the shift advection terms in gamma driver",
+    Default -> 1
+  },
+  {
+    Name -> MinimumLapse,
+    Description -> "Minimum value of the lapse function",
+    Default -> -1
+  },
+  {
+    Name -> SpatialBetaDriverRadius,
+    Description -> "Radius at which the BetaDriver starts to be reduced",
+    AllowedValues -> {{Value -> "(0:*", Description -> "Positive"}},
+    Default -> 10^12
+  },
+  {
+    Name -> SpatialShiftGammaCoeffRadius,
+    Description -> "Radius at which the ShiftGammaCoefficient starts to be reduced",
+    AllowedValues -> {{Value -> "(0:*", Description -> "Positive"}},
+    Default -> 10^12
+  },
+  {
+    Name -> EpsDiss,
+    Description -> "Dissipation strength",
+    AllowedValues -> {{Value -> "(0:*", Description -> "Positive"}},
+    Default -> 0
+  }
+};
+
+(******************************************************************************)
+(* Construct the thorns *)
+(******************************************************************************)
+
+calculations =
+Join[
+{
+  initialCalc,
+  convertFromADMBaseCalc,
+  initGammaCalc,
+  convertFromADMBaseGammaCalc,
+  (* evolCalc, *)
+  evolCalc1, evolCalc2,
+  dissCalc,
+  advectCalc,
+  initRHSCalc,
+  (* evol1Calc, evol2Calc, *)
+  RHSStaticBoundaryCalc,
+  (* RHSRadiativeBoundaryCalc, *)
+  enforceCalc,
+  boundaryCalc,
+  convertToADMBaseCalc,
+  convertToADMBaseDtLapseShiftCalc,
+  convertToADMBaseDtLapseShiftBoundaryCalc,
+  convertToADMBaseFakeDtLapseShiftCalc,
+  (* constraintsCalc, *)
+  constraintsCalc1, constraintsCalc2
+},
+  {} (*dissCalcs*)
+];
+
+calculations = Map[Append[#, NoSimplify -> True] &, calculations];
+
+CreateKrancThornTT [groups, "TestThorns", thorn,
+  Calculations -> calculations,
+  DeclaredGroups -> declaredGroupNames,
+  PartialDerivatives -> derivatives,
+  EvolutionTimelevels -> evolutionTimelevels,
+  DefaultEvolutionTimelevels -> 3,
+  UseJacobian -> True,
+  UseLoopControl -> True,
+  UseVectors -> True,
+  InheritedImplementations -> inheritedImplementations,
+  InheritedKeywordParameters -> inheritedKeywordParameters,
+  ExtendedKeywordParameters -> extendedKeywordParameters,
+  KeywordParameters -> keywordParameters,
+  IntParameters -> intParameters,
+  RealParameters -> realParameters
+];
+
+];
+
+(****************************************************************)
+(* McLachlan *)
+(****************************************************************)
+
+Test[
+
+(******************************************************************************)
+(* Options *)
+(******************************************************************************)
+
+(* These are the arguments to createCode:
+   - 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
+     (matter seems cheap; it should be always enabled)
+   - thorn base name
+*)
+
+createCode[4, False, True , 3, 1, "BSSN"];
+(* createCode[4, False, True , 3, 1, "CCZ4"]; *)
+	,
+	Null
+	,
+	TestID->"McLachlan"
+]
diff --git a/Tests/RunTests.m b/Tests/RunTests.m
index 46fee00..5d749ac 100755
--- a/Tests/RunTests.m
+++ b/Tests/RunTests.m
@@ -21,7 +21,8 @@ Needs["KrancThorn`"];
 SetDebugLevel[DebugQuiet];
 
 tests = {
-		"Kranc"
+  "Kranc",
+  "McLachlan"
 };
 
 (Print["\n"]; TestRun[#<>".mt", Loggers -> {VerbosePrintLogger[]}, TestRunTitle -> #]) & /@ tests;