Generate a separate thorn for CCZ4.

This still uses the same Kranc script as the standard BSSN code, but having a separate thorn means that the BSSN is not negatively affected performance-wise.

2 files changed, 136 insertions, 150 deletions

diff --git a/m/Makefile b/m/Makefile
index 3af082a..cf8260b 100644
--- a/m/Makefile
+++ b/m/Makefile
@@ -13,13 +13,14 @@ McLachlan_ADM.out: McLachlan_ADM.m
 	done
 
 McLachlan_BSSN.out: McLachlan_BSSN.m
-	rm -rf ML_BSSN*
+	rm -rf ML_BSSN* ML_CCZ4*
 	./runmath.sh $^
 	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
+	./copy-if-changed.sh ML_CCZ4 ../ML_CCZ4
 
 McLachlan_ADMConstraints.out: McLachlan_ADMConstraints.m
 	rm -rf ML_ADMConstraints*
diff --git a/m/McLachlan_BSSN.m b/m/McLachlan_BSSN.m
index 16d076a..0ce97b7 100644
--- a/m/McLachlan_BSSN.m
+++ b/m/McLachlan_BSSN.m
@@ -6,8 +6,8 @@ SetSourceLanguage["C"];
 (* Options *)
 (******************************************************************************)
 
-createCode[derivOrder_, useJacobian_, splitUpwindDerivs_, evolutionTimelevels_, addMatter_] :=
-Module[{},
+createCode[derivOrder_, useJacobian_, splitUpwindDerivs_, evolutionTimelevels_, addMatter_, formulation_] :=
+Module[{prefix, suffix, thorn},
 
 prefix = "ML_";
 suffix =
@@ -19,9 +19,10 @@ suffix =
   (* <> If [addMatter==1, "_M", ""] *)
   ;
 
-BSSN = prefix <> "BSSN" <> suffix;
-
+thorn = prefix <> formulation <> suffix;
 
+SetAttributes[IfCCZ4, HoldAll];
+IfCCZ4[expr_, else_:Sequence[]] := If[formulation === "CCZ4", expr, Unevaluated[else]];
 
 (******************************************************************************)
 (* Derivatives *)
@@ -165,11 +166,6 @@ T11=eTxx; T12=eTxy; T22=eTyy; T13=eTxz; T23=eTyz; T33=eTzz;
 CMBSSNphi = 0;
 CMBSSNW   = 1;
 
-(* enum constants for formulation; these must be consistent
-   with the definition of the Cactus parameter formulation *)
-FBSSN    = 0;
-FCCZ4    = 1;
-
 detgExpr  = Det [MatrixOfComponents [g [la,lb]]];
 ddetgExpr[la_] =
   Sum [D[Det[MatrixOfComponents[g[la, lb]]], X] PD[X, la],
@@ -193,13 +189,13 @@ evolvedGroups =
   {SetGroupName [CreateGroupFromTensor [phi      ], prefix <> "log_confac"],
    SetGroupName [CreateGroupFromTensor [gt[la,lb]], prefix <> "metric"    ],
    SetGroupName [CreateGroupFromTensor [Xt[ua]   ], prefix <> "Gamma"     ],
-   SetGroupName [CreateGroupFromTensor [Theta    ], prefix <> "Theta"     ],
    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"   ]};
+   SetGroupName [CreateGroupFromTensor [B[ua]    ], prefix <> "dtshift"   ],
+   IfCCZ4[SetGroupName[CreateGroupFromTensor[Theta], prefix <> "Theta"]]};
 evaluatedGroups =
   {SetGroupName [CreateGroupFromTensor [H      ], prefix <> "Ham"],
    SetGroupName [CreateGroupFromTensor [M[la]  ], prefix <> "mom"],
@@ -235,7 +231,7 @@ groups = Join [declaredGroups, extraGroups];
 
 initialCalc =
 {
-  Name -> BSSN <> "_Minkowski",
+  Name -> thorn <> "_Minkowski",
   Schedule -> {"IN ADMBase_InitialData"},
   ConditionalOnKeyword -> {"my_initial_data", "Minkowski"},
   Equations -> 
@@ -245,11 +241,11 @@ initialCalc =
     trK       -> 0,
     At[la,lb] -> 0,
     Xt[ua]    -> 0,
-    Theta     -> 0,
     alpha     -> 1,
     A         -> 0,
     beta[ua]  -> 0,
-    B[ua]     -> 0
+    B[ua]     -> 0,
+    IfCCZ4[Theta -> 0]
   }
 };
 
@@ -286,7 +282,7 @@ Module[
 
 convertFromADMBaseCalc =
 {
-  Name -> BSSN <> "_convertFromADMBase",
+  Name -> thorn <> "_convertFromADMBase",
   Schedule -> {"AT initial AFTER ADMBase_PostInitial"},
   ConditionalOnKeyword -> {"my_initial_data", "ADMBase"},
   Shorthands -> {g[la,lb], detg, gu[ua,ub], em4phi},
@@ -303,18 +299,18 @@ convertFromADMBaseCalc =
     trK       -> gu[ua,ub] admK[la,lb],
     At[la,lb] -> em4phi (admK[la,lb] - (1/3) g[la,lb] trK),
     
-    Theta     -> 0,
-    
     alpha     -> admalpha,
     
-    beta[ua]  -> admbeta[ua]
+    beta[ua]  -> admbeta[ua],
+
+    IfCCZ4[Theta -> 0]
   }
 };
 
 convertFromADMBaseGammaCalc =
 {
-  Name -> BSSN <> "_convertFromADMBaseGamma",
-  Schedule -> {"AT initial AFTER " <> BSSN <> "_convertFromADMBase"},
+  Name -> thorn <> "_convertFromADMBaseGamma",
+  Schedule -> {"AT initial AFTER " <> thorn <> "_convertFromADMBase"},
   ConditionalOnKeyword -> {"my_initial_data", "ADMBase"},
   (*
   Where -> InteriorNoSync,
@@ -370,8 +366,8 @@ convertFromADMBaseGammaCalc =
    initialised. *)
 initGammaCalc =
 {
-  Name -> BSSN <> "_InitGamma",
-  Schedule -> {"AT initial BEFORE " <> BSSN <> "_convertFromADMBaseGamma"},
+  Name -> thorn <> "_InitGamma",
+  Schedule -> {"AT initial BEFORE " <> thorn <> "_convertFromADMBaseGamma"},
   ConditionalOnKeyword -> {"my_initial_data", "ADMBase"},
   Where -> Everywhere,
   Equations -> 
@@ -390,13 +386,13 @@ initGammaCalc =
 
 convertToADMBaseCalc =
 {
-  Name -> BSSN <> "_convertToADMBase",
-  Schedule -> {"IN " <> BSSN <> "_convertToADMBaseGroup"},
+  Name -> thorn <> "_convertToADMBase",
+  Schedule -> {"IN " <> thorn <> "_convertToADMBaseGroup"},
   Where -> Everywhere,
   Shorthands -> {e4phi},
   Equations -> 
   {
-    e4phi       -> IfThen[conformalMethod==CMBSSNW, 1/phi^2, Exp[4 phi]],
+    e4phi       -> IfThen[conformalMethod, 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,
@@ -406,8 +402,8 @@ convertToADMBaseCalc =
 
 convertToADMBaseDtLapseShiftCalc =
 {
-  Name -> BSSN <> "_convertToADMBaseDtLapseShift",
-  Schedule -> {"IN " <> BSSN <> "_convertToADMBaseGroup"},
+  Name -> thorn <> "_convertToADMBaseDtLapseShift",
+  Schedule -> {"IN " <> thorn <> "_convertToADMBaseGroup"},
   ConditionalOnKeyword -> {"dt_lapse_shift_method", "correct"},
   Where -> Interior,
   Shorthands -> {dir[ua], detgt, gtu[ua,ub], eta, theta},
@@ -431,7 +427,7 @@ convertToADMBaseDtLapseShiftCalc =
     admdtalpha -> - harmonicF alpha^harmonicN
                     (+ LapseACoeff       A
                      + ((1 - LapseACoeff)
-                        (trK - IfThen[formulation==FCCZ4, 2 Theta, 0])))
+                        (trK - IfCCZ4[2 Theta, 0])))
                   + LapseAdvectionCoeff Upwind[beta[ua], alpha, la],
     admdtbeta[ua] -> IfThen[harmonicShift,
                             - 1/2 gtu[ua,uj] phi alpha
@@ -450,8 +446,8 @@ convertToADMBaseDtLapseShiftCalc =
 
 convertToADMBaseDtLapseShiftBoundaryCalc =
 {
-  Name -> BSSN <> "_convertToADMBaseDtLapseShiftBoundary",
-  Schedule -> {"IN " <> BSSN <> "_convertToADMBaseGroup"},
+  Name -> thorn <> "_convertToADMBaseDtLapseShiftBoundary",
+  Schedule -> {"IN " <> thorn <> "_convertToADMBaseGroup"},
   ConditionalOnKeyword -> {"dt_lapse_shift_method", "correct"},
   Where -> BoundaryWithGhosts,
   Shorthands -> {detgt, gtu[ua,ub], eta, theta},
@@ -472,7 +468,7 @@ convertToADMBaseDtLapseShiftBoundaryCalc =
     admdtalpha -> - harmonicF alpha^harmonicN
                     (+ LapseACoeff       A
                      + ((1 - LapseACoeff)
-                        (trK - IfThen[formulation==FCCZ4, 2 Theta, 0]))),
+                        (trK - IfCCZ4[2 Theta, 0]))),
     admdtbeta[ua] -> IfThen[harmonicShift,
                             0,
                             (* else *)
@@ -485,8 +481,8 @@ convertToADMBaseDtLapseShiftBoundaryCalc =
 
 convertToADMBaseFakeDtLapseShiftCalc =
 {
-  Name -> BSSN <> "_convertToADMBaseFakeDtLapseShift",
-  Schedule -> {"IN " <> BSSN <> "_convertToADMBaseGroup"},
+  Name -> thorn <> "_convertToADMBaseFakeDtLapseShift",
+  Schedule -> {"IN " <> thorn <> "_convertToADMBaseGroup"},
   ConditionalOnKeyword -> {"dt_lapse_shift_method", "noLapseShiftAdvection"},
   Where -> Everywhere,
   Shorthands -> {detgt, gtu[ua,ub], eta, theta},
@@ -509,7 +505,7 @@ convertToADMBaseFakeDtLapseShiftCalc =
     admdtalpha -> - harmonicF alpha^harmonicN
                     (+ LapseACoeff       A
                      + ((1 - LapseACoeff)
-                        (trK - IfThen[formulation==FCCZ4, 2 Theta, 0]))),
+                        (trK - IfCCZ4[2 Theta, 0]))),
     admdtbeta[ua] -> IfThen[harmonicShift,
                             0,
                             (* else *)
@@ -526,8 +522,8 @@ convertToADMBaseFakeDtLapseShiftCalc =
 
 evolCalc =
 {
-  Name -> BSSN <> "_RHS",
-  Schedule -> {"IN " <> BSSN <> "_evolCalcGroup"},
+  Name -> thorn <> "_RHS",
+  Schedule -> {"IN " <> thorn <> "_evolCalcGroup"},
   (*
   Where -> Interior,
   *)
@@ -536,14 +532,14 @@ evolCalc =
      radiative boundary conditions.  *)
   Where -> InteriorNoSync,
   Shorthands -> {dir[ua],
-                 detgt, gtu[ua,ub], Z[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], dotTheta,
-                 epsdiss[ua]},
+                 rho, S[la], trS, fac1, fac2, dottrK, dotXt[ua],
+                 epsdiss[ua], IfCCZ4[Z[ua]], IfCCZ4[dotTheta]},
   Equations -> 
   {
     dir[ua] -> Sign[beta[ua]],
@@ -569,10 +565,9 @@ evolCalc =
 
     (* The Z quantities *)
     (* gr-qc:1106.2254 (2011), eqn. (23) *)
-    Z[ud] -> IfThen[formulation==FCCZ4,
-                    (1/2) gu[ua,ud] (- PD[gt[la,lb],lc] gtu[ub,uc]
-                                     + gt[la,lc] Xt[uc]),
-                    0],
+    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 *)
@@ -599,13 +594,12 @@ evolCalc =
     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]
-                + IfThen[formulation==FCCZ4,
-                         + (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],
-                         0],
+    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 *)
     
@@ -629,10 +623,7 @@ evolCalc =
     (* 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]
-                                 - IfThen[formulation==FCCZ4,
-                                          (1/3) At[lc,ld] gtu[uc,ud] gt[la,lb],
-                                          0])
+    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) *)
@@ -648,7 +639,7 @@ evolCalc =
                       + (1/3) gtu[ui,uj] PD[beta[ul],lj,ll]
                       - Xtn[uj] PD[beta[ui],lj] 
                       + (2/3) Xtn[ui] PD[beta[uj],lj]
-                      + IfThen[formulation==FCCZ4,
+                      + IfCCZ4[
                                + IfThen[GammaShift,
                                         2 e4phi (- Z[uj] PD[beta[ui],lj]
                                                  + (2/3) Z[ui] PD[beta[uj],lj]),
@@ -657,21 +648,21 @@ evolCalc =
                                + 2 gtu[ui,uj] (+ alpha PD[Theta,lj]
                                                - Theta PD[alpha,lj])
                                - 2 alpha e4phi dampk1 Z[ui],
-                               0]
+                        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) *)
-    dotTheta       -> IfThen[formulation==FCCZ4,
-                             - PD[alpha,la] Z[ua]
-                             + (1/2) alpha (gu[ua,ub] R[la,lb]
-                                            - Atm[ua,lb] Atm[ub,la]
-                                            + (2/3) trK^2 - 2 trK Theta)
-                             - dampk1 (2 + dampk2) alpha Theta,
-                             0],
+    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)
+    ],
 
-    dot[Theta] -> dotTheta,
+    IfCCZ4[
+      dot[Theta] -> dotTheta
+    ],
 
     (* PRD 62, 044034 (2000), eqn. (11) *)
     (* gr-qc:1106.2254 (2011), eqn. (17) *)
@@ -681,7 +672,7 @@ evolCalc =
                                 + 2 cdphi[la] PD[alpha,lb] )
                                 - Xtn[ua] PD[alpha,la] )
                       + alpha (Atm[ua,lb] Atm[ub,la] + (1/3) trK^2)
-                      + IfThen[formulation==FCCZ4,
+                      + IfCCZ4[
                                + 2 dotTheta + 2 PD[alpha,la] Z[ua]
                                + dampk1 (1 - dampk2) alpha Theta,
                                0]
@@ -698,9 +689,7 @@ evolCalc =
                       + 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 - IfThen[formulation==FCCZ4,
-                                                2 Theta,
-                                                0])
+                      + 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]
@@ -723,11 +712,11 @@ evolCalc =
     dot[alpha] -> - harmonicF alpha^harmonicN
                     (+ LapseACoeff       A
                      + ((1 - LapseACoeff)
-                        (+ trK - IfThen[formulation==FCCZ4, 2 Theta, 0]
+                        (+ trK - IfCCZ4[2 Theta, 0]
                          + AlphaDriver (alpha - 1)))),
 
     dot[A]     -> + (LapseACoeff
-                     (+ dottrK - IfThen[formulation==FCCZ4, 2 dotTheta, 0]
+                     (+ dottrK - IfCCZ4[2 dotTheta, 0]
                       - AlphaDriver A)),
     
     eta -> etaExpr,
@@ -754,9 +743,9 @@ evolCalc =
 
 advectCalc =
 {
-  Name -> BSSN <> "_Advect",
-  Schedule -> {"IN " <> BSSN <> "_evolCalcGroup " <>
-               "AFTER (" <> BSSN <> "_RHS1 " <> BSSN <> "_RHS2)"},
+  Name -> thorn <> "_Advect",
+  Schedule -> {"IN " <> thorn <> "_evolCalcGroup " <>
+               "AFTER (" <> thorn <> "_RHS1 " <> thorn <> "_RHS2)"},
   (*
   Where -> Interior,
   *)
@@ -775,10 +764,9 @@ advectCalc =
 
        dot[Xt[ui]] -> dot[Xt[ui]] + Upwind[beta[uj], Xt[ui], lj],
 
-        dot[Theta] -> dot[Theta]
-                      + IfThen[formulation==FCCZ4,
-                               Upwind[beta[ua], Theta, la],
-                               0],
+    IfCCZ4[
+        dot[Theta] -> dot[Theta] + Upwind[beta[ua], Theta, la]
+    ],
 
           dot[trK] -> dot[trK] + Upwind[beta[ua], trK, la],
 
@@ -811,12 +799,12 @@ evolCalc1 = PartialCalculation[evolCalc, "1",
     dot[phi],
     dot[gt[la,lb]],
     dot[Xt[ui]],
-    dot[Theta],
     dot[trK],
     dot[alpha],
     dot[A],
     dot[beta[ua]],
-    dot[B[ua]]
+    dot[B[ua]],
+    IfCCZ4[dot[Theta]]
   }];
 
 evolCalc2 = PartialCalculation[evolCalc, "2",
@@ -829,9 +817,9 @@ evolCalc2 = PartialCalculation[evolCalc, "2",
 
 dissCalc =
 {
-  Name -> BSSN <> "_Dissipation",
-  Schedule -> {"IN " <> BSSN <> "_evolCalcGroup " <>
-               "AFTER (" <> BSSN <> "_RHS1 " <> BSSN <> "_RHS2)"},
+  Name -> thorn <> "_Dissipation",
+  Schedule -> {"IN " <> thorn <> "_evolCalcGroup " <>
+               "AFTER (" <> thorn <> "_RHS1 " <> thorn <> "_RHS2)"},
   ConditionalOnKeyword -> {"apply_dissipation", "always"},
   Where -> InteriorNoSync,
   Shorthands -> {epsdiss[ua]},
@@ -840,7 +828,7 @@ dissCalc =
     epsdiss[ua] -> EpsDiss,
     Sequence@@Table[
       dot[var]       -> dot[var] + epsdiss[ux] PDdiss[var,lx],
-      {var, {phi, gt[la,lb], Xt[ui], Theta, trK, At[la,lb],
+      {var, {phi, gt[la,lb], Xt[ui], IfCCZ4[Theta], trK, At[la,lb],
              alpha, A, beta[ua], B[ua]}}]
   }
 };
@@ -848,9 +836,9 @@ dissCalc =
 dissCalcs =
 Table[
 {
-  Name -> BSSN <> "_Dissipation_" <> ToString[var /. {Tensor[n_,__] -> n}],
-  Schedule -> {"IN " <> BSSN <> "_evolCalcGroup " <>
-               "AFTER (" <> BSSN <> "_RHS1 " <> BSSN <> "_RHS2)"},
+  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]},
@@ -860,13 +848,13 @@ Table[
     dot[var]    -> dot[var] + epsdiss[ux] PDdiss[var,lx]
   }
 },
-  {var, {phi, gt[la,lb], Xt[ui], Theta, trK, At[la,lb],
+  {var, {phi, gt[la,lb], Xt[ui], IfCCZ4[Theta], trK, At[la,lb],
          alpha, A, beta[ua], B[ua]}}
 ];
 
 RHSStaticBoundaryCalc =
 {
-  Name -> BSSN <> "_RHSStaticBoundary",
+  Name -> thorn <> "_RHSStaticBoundary",
   Schedule -> {"IN MoL_CalcRHS"},
   ConditionalOnKeyword -> {"my_rhs_boundary_condition", "static"},
   Where -> Boundary,
@@ -877,11 +865,11 @@ RHSStaticBoundaryCalc =
     dot[trK]       -> 0,
     dot[At[la,lb]] -> 0,
     dot[Xt[ua]]    -> 0,
-    dot[Theta]     -> 0,
     dot[alpha]     -> 0,
     dot[A]         -> 0,
     dot[beta[ua]]  -> 0,
-    dot[B[ua]]     -> 0
+    dot[B[ua]]     -> 0,
+    IfCCZ4[dot[Theta] -> 0]
   }
 };
 
@@ -890,8 +878,8 @@ RHSStaticBoundaryCalc =
    them to be zero *)
 initRHSCalc =
 {
-  Name -> BSSN <> "_InitRHS",
-  Schedule -> {"AT analysis BEFORE " <> BSSN <> "_evolCalcGroup"},
+  Name -> thorn <> "_InitRHS",
+  Schedule -> {"AT analysis BEFORE " <> thorn <> "_evolCalcGroup"},
   Where -> Everywhere,
   Equations -> 
   {
@@ -900,17 +888,17 @@ initRHSCalc =
     dot[trK]       -> 0,
     dot[At[la,lb]] -> 0,
     dot[Xt[ua]]    -> 0,
-    dot[Theta]     -> 0,
     dot[alpha]     -> 0,
     dot[A]         -> 0,
     dot[beta[ua]]  -> 0,
-    dot[B[ua]]     -> 0
+    dot[B[ua]]     -> 0,
+    IfCCZ4[dot[Theta] -> 0]
   }
 };
 
 RHSRadiativeBoundaryCalc =
 {
-  Name -> BSSN <> "_RHSRadiativeBoundary",
+  Name -> thorn <> "_RHSRadiativeBoundary",
   Schedule -> {"IN MoL_CalcRHS"},
   ConditionalOnKeyword -> {"my_rhs_boundary_condition", "radiative"},
   Where -> Boundary,
@@ -940,17 +928,19 @@ RHSRadiativeBoundaryCalc =
     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[Theta]     -> - vg su[uc] PDo[Theta    ,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]
+    dot[B[ua]]     -> -    su[uc] PDo[B[ua]    ,lc],
+    IfCCZ4[
+      dot[Theta]   -> - vg su[uc] PDo[Theta    ,lc]
+    ]
   }
 };
 
 enforceCalc =
 {
-  Name -> BSSN <> "_enforce",
+  Name -> thorn <> "_enforce",
   Schedule -> {"IN MoL_PostStepModify"},
   Shorthands -> {detgt, gtu[ua,ub], trAt},
   Equations -> 
@@ -981,7 +971,7 @@ enforceCalc =
 
 boundaryCalc =
 {
-  Name -> BSSN <> "_boundary",
+  Name -> thorn <> "_boundary",
   Schedule -> {"IN MoL_PostStep"},
   ConditionalOnKeyword -> {"my_boundary_condition", "Minkowski"},
   Where -> BoundaryWithGhosts,
@@ -992,11 +982,11 @@ boundaryCalc =
     trK       -> 0,
     At[la,lb] -> 0,
     Xt[ua]    -> 0,
-    Theta     -> 0,
     alpha     -> 1,
     A         -> 0,
     beta[ua]  -> 0,
-    B[ua]     -> 0
+    B[ua]     -> 0,
+    IfCCZ4[Theta -> 0]
   }
 };
 
@@ -1006,7 +996,7 @@ boundaryCalc =
 
 constraintsCalc =
 {
-  Name -> BSSN <> "_constraints",
+  Name -> thorn <> "_constraints",
   Schedule -> Automatic,
   After -> "MoL_PostStep",
   Where -> Interior,
@@ -1040,10 +1030,9 @@ constraintsCalc =
     gu[ua,ub]   -> em4phi gtu[ua,ub],
 
     (* The Z quantities *)
-    Z[ud] -> IfThen[formulation==FCCZ4,
-                    (1/2) gu[ua,ud] (- PD[gt[la,lb],lc] gtu[ub,uc]
-                                     + gt[la,lc] Xt[uc]),
-                    0],
+    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]
@@ -1095,13 +1084,14 @@ constraintsCalc =
                    + 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]
-                + IfThen[formulation==FCCZ4,
-                         + (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],
-                         0],
+    R[la,lb] -> + Rt[la,lb] + Rphi[la,lb],
+
+    IfCCZ4[
+      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, *)
@@ -1177,23 +1167,23 @@ extendedKeywordParameters =
 {
   {
     Name -> "ADMBase::evolution_method",
-    AllowedValues -> {BSSN}
+    AllowedValues -> {thorn}
   },
   {
     Name -> "ADMBase::lapse_evolution_method",
-    AllowedValues -> {BSSN}
+    AllowedValues -> {thorn}
   },
   {
     Name -> "ADMBase::shift_evolution_method",
-    AllowedValues -> {BSSN}
+    AllowedValues -> {thorn}
   },
   {
     Name -> "ADMBase::dtlapse_evolution_method",
-    AllowedValues -> {BSSN}
+    AllowedValues -> {thorn}
   },
   {
     Name -> "ADMBase::dtshift_evolution_method",
-    AllowedValues -> {BSSN}
+    AllowedValues -> {thorn}
   }
 };
 
@@ -1263,25 +1253,11 @@ keywordParameters =
 
 intParameters =
 {
-  {
-    Name -> conformalMethod,
-    Description -> "Treatment of conformal factor",
-    AllowedValues -> {{Value -> "0", Description -> "phi method"},
-                      {Value -> "1", Description -> "W method"}},
-    Default -> 0
-  },
-  {
-    Name -> formulation,
-    Description -> "Formulation of the equations",
-    AllowedValues -> {{Value -> "0", Description -> "BSSN"},
-                      {Value -> "1", Description -> "CCZ4"}},
-    Default -> 0
-  },
-  {
+  IfCCZ4[{
     Name -> GammaShift,
     Description -> "Covariant shift term in Gamma",
     Default -> 0
-  },
+  }],
   {
     Name -> harmonicN,
     Description -> "d/dt alpha = - f alpha^n K  (harmonic=2, 1+log=1)",
@@ -1292,31 +1268,38 @@ intParameters =
     Default -> 0
   },
   {
-    Name -> harmonicShift,
-    Description -> "Whether to use the harmonic shift",
-    AllowedValues -> {{Value -> "0", Description -> "Gamma driver shift"},
-                      {Value -> "1", Description -> "Harmonic shift"}},
+    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 -> 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",
@@ -1410,7 +1393,7 @@ Join[
   {} (*dissCalcs*)
 ];
 
-CreateKrancThornTT [groups, ".", BSSN,
+CreateKrancThornTT [groups, ".", thorn,
   Calculations -> calculations,
   DeclaredGroups -> declaredGroupNames,
   PartialDerivatives -> derivatives,
@@ -1443,9 +1426,11 @@ CreateKrancThornTT [groups, ".", BSSN,
 (* matter: 0 or 1
    (matter seems cheap; it should be always enabled) *)
 
-createCode[2, False, True , 3, 1];
-createCode[4, False, True , 3, 1];
-createCode[4, False, False, 3, 1];
-createCode[4, True , True , 3, 1];
-createCode[8, False, True , 3, 1];
-createCode[8, True , True , 3, 1];
+createCode[2, False, True , 3, 1, "BSSN"];
+createCode[4, False, True , 3, 1, "BSSN"];
+createCode[4, False, False, 3, 1, "BSSN"];
+createCode[4, True , True , 3, 1, "BSSN"];
+createCode[8, False, True , 3, 1, "BSSN"];
+createCode[8, True , True , 3, 1, "BSSN"];
+
+createCode[4, False, True , 3, 1, "CCZ4"];