Add formulation parameter to select between BSSN and CCZ4.

The default selects BSSN for backwards compatibility.

1 files changed, 29 insertions, 19 deletions

diff --git a/m/McLachlan_BSSN.m b/m/McLachlan_BSSN.m
index e6dd50b..16d076a 100644
--- a/m/McLachlan_BSSN.m
+++ b/m/McLachlan_BSSN.m
@@ -164,7 +164,11 @@ T11=eTxx; T12=eTxy; T22=eTyy; T13=eTxz; T23=eTyz; T33=eTzz;
    with the definition of the Cactus parameter conformalMethod *)
 CMBSSNphi = 0;
 CMBSSNW   = 1;
-CMCCZ4    = 2;
+
+(* 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_] =
@@ -427,7 +431,7 @@ convertToADMBaseDtLapseShiftCalc =
     admdtalpha -> - harmonicF alpha^harmonicN
                     (+ LapseACoeff       A
                      + ((1 - LapseACoeff)
-                        (trK - IfThen[conformalMethod==CMCCZ4, 2 Theta, 0])))
+                        (trK - IfThen[formulation==FCCZ4, 2 Theta, 0])))
                   + LapseAdvectionCoeff Upwind[beta[ua], alpha, la],
     admdtbeta[ua] -> IfThen[harmonicShift,
                             - 1/2 gtu[ua,uj] phi alpha
@@ -468,7 +472,7 @@ convertToADMBaseDtLapseShiftBoundaryCalc =
     admdtalpha -> - harmonicF alpha^harmonicN
                     (+ LapseACoeff       A
                      + ((1 - LapseACoeff)
-                        (trK - IfThen[conformalMethod==CMCCZ4, 2 Theta, 0]))),
+                        (trK - IfThen[formulation==FCCZ4, 2 Theta, 0]))),
     admdtbeta[ua] -> IfThen[harmonicShift,
                             0,
                             (* else *)
@@ -505,7 +509,7 @@ convertToADMBaseFakeDtLapseShiftCalc =
     admdtalpha -> - harmonicF alpha^harmonicN
                     (+ LapseACoeff       A
                      + ((1 - LapseACoeff)
-                        (trK - IfThen[conformalMethod==CMCCZ4, 2 Theta, 0]))),
+                        (trK - IfThen[formulation==FCCZ4, 2 Theta, 0]))),
     admdtbeta[ua] -> IfThen[harmonicShift,
                             0,
                             (* else *)
@@ -565,7 +569,7 @@ evolCalc =
 
     (* The Z quantities *)
     (* gr-qc:1106.2254 (2011), eqn. (23) *)
-    Z[ud] -> IfThen[conformalMethod==CMCCZ4,
+    Z[ud] -> IfThen[formulation==FCCZ4,
                     (1/2) gu[ua,ud] (- PD[gt[la,lb],lc] gtu[ub,uc]
                                      + gt[la,lc] Xt[uc]),
                     0],
@@ -596,7 +600,7 @@ evolCalc =
     Atu[ua,ub] -> Atm[ua,lc] gtu[ub,uc],
     
     R[la,lb] -> + Rt[la,lb] + Rphi[la,lb]
-                + IfThen[conformalMethod==CMCCZ4,
+                + 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])
@@ -626,7 +630,7 @@ evolCalc =
     (* 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[conformalMethod==CMCCZ4,
+                                 - IfThen[formulation==FCCZ4,
                                           (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]
@@ -644,7 +648,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[conformalMethod==CMCCZ4,
+                      + IfThen[formulation==FCCZ4,
                                + IfThen[GammaShift,
                                         2 e4phi (- Z[uj] PD[beta[ui],lj]
                                                  + (2/3) Z[ui] PD[beta[uj],lj]),
@@ -659,7 +663,7 @@ evolCalc =
     dot[Xt[ui]]    -> dotXt[ui],
 
     (* gr-qc:1106.2254 (2011), eqn. (18) *)
-    dotTheta       -> IfThen[conformalMethod==CMCCZ4,
+    dotTheta       -> IfThen[formulation==FCCZ4,
                              - PD[alpha,la] Z[ua]
                              + (1/2) alpha (gu[ua,ub] R[la,lb]
                                             - Atm[ua,lb] Atm[ub,la]
@@ -677,7 +681,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[conformalMethod==CMCCZ4,
+                      + IfThen[formulation==FCCZ4,
                                + 2 dotTheta + 2 PD[alpha,la] Z[ua]
                                + dampk1 (1 - dampk2) alpha Theta,
                                0]
@@ -694,7 +698,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[conformalMethod==CMCCZ4,
+                      + alpha (+ ((trK - IfThen[formulation==FCCZ4,
                                                 2 Theta,
                                                 0])
                                   At[la,lb])
@@ -719,11 +723,11 @@ evolCalc =
     dot[alpha] -> - harmonicF alpha^harmonicN
                     (+ LapseACoeff       A
                      + ((1 - LapseACoeff)
-                        (+ trK - IfThen[conformalMethod==CMCCZ4, 2 Theta, 0]
+                        (+ trK - IfThen[formulation==FCCZ4, 2 Theta, 0]
                          + AlphaDriver (alpha - 1)))),
 
     dot[A]     -> + (LapseACoeff
-                     (+ dottrK - IfThen[conformalMethod==CMCCZ4, 2 dotTheta, 0]
+                     (+ dottrK - IfThen[formulation==FCCZ4, 2 dotTheta, 0]
                       - AlphaDriver A)),
     
     eta -> etaExpr,
@@ -772,7 +776,7 @@ advectCalc =
        dot[Xt[ui]] -> dot[Xt[ui]] + Upwind[beta[uj], Xt[ui], lj],
 
         dot[Theta] -> dot[Theta]
-                      + IfThen[conformalMethod==CMCCZ4,
+                      + IfThen[formulation==FCCZ4,
                                Upwind[beta[ua], Theta, la],
                                0],
 
@@ -1036,7 +1040,7 @@ constraintsCalc =
     gu[ua,ub]   -> em4phi gtu[ua,ub],
 
     (* The Z quantities *)
-    Z[ud] -> IfThen[conformalMethod==CMCCZ4,
+    Z[ud] -> IfThen[formulation==FCCZ4,
                     (1/2) gu[ua,ud] (- PD[gt[la,lb],lc] gtu[ub,uc]
                                      + gt[la,lc] Xt[uc]),
                     0],
@@ -1092,7 +1096,7 @@ constraintsCalc =
                                  - (1/3) g[lb,lc] gu[ua,ud] ddetg[ld]),
     
     R[la,lb] -> + Rt[la,lb] + Rphi[la,lb]
-                + IfThen[conformalMethod==CMCCZ4,
+                + 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])
@@ -1262,9 +1266,15 @@ intParameters =
   {
     Name -> conformalMethod,
     Description -> "Treatment of conformal factor",
-    AllowedValues -> {{Value -> "0", Description -> "BSSN phi method"},
-                      {Value -> "1", Description -> "BSSN W method"},
-                      {Value -> "2", Description -> "CCZ4 method"}},
+    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
   },
   {