
SetEnhancedTimes[False];

(******************************************************************************)
(* Options *)
(******************************************************************************)

createCode[derivOrder_, useJacobian_, splitUpwindDerivs_, evolutionTimelevels_, addMatter_] :=
Module[{},

prefix = "CL_";
suffix =
  ""
  <> If [useJacobian, "_MP", ""]
  <> If [derivOrder!=4, "_O" <> ToString[derivOrder], ""]
  <> If [splitUpwindDerivs, "", "_UPW"]
  (* <> If [evolutionTimelevels!=3, "_TL" <> ToString[evolutionTimelevels], ""] *)
  (* <> If [addMatter==1, "_M", ""] *)
  ;

BSSN = prefix <> "BSSN" <> suffix;



(******************************************************************************)
(* 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_] ->
    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,
      J, dJ,
      admg, admK, admalpha, admdtalpha, admbeta, admdtbeta, H, M,
      g, detg, gu, G, R, trR, Km, trK, cdphi,
      phi, dphi, gt, dgt, trdgt, cdgt, At, Xt, Xtn, alpha, dalpha,
      cdalpha,
      beta, dbeta, cdbeta, B,
      Atm, Atu, trA, Ats, trAts,
      dgts, trdgts, trcdgt,
      cXt, cS, cA,
      e4phi, em4phi, ddetg, detgt, gtu, ddetgt, dgtu, ddgtu, Gtl, Gtlu, Gt,
      Rt, Rphi, gK,
      T00, T0, T, rho, S,
      x, y, z, r}];

(* 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],
      gt[la,lb], At[la,lb], Ats[la,lb], Rt[la,lb], Rphi[la,lb], T[la,lb]}];
AssertSymmetricIncreasing [dgt[la,lb,lc], lb, lc];
AssertSymmetricIncreasing [cdgt[la,lb,lc], lb, lc];
AssertSymmetricIncreasing [dJ[ua,lb,lc], lb, lc];
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];
AssertSymmetricIncreasing [dgts[la,lb,lc], lb, lc];
Map [AssertSymmetricDecreasing,
     {gu[ua,ub], gtu[ua,ub], Atu[ua,ub]}];
AssertSymmetricDecreasing [dgtu[ua,ub,lc], ua, ub];
AssertSymmetricDecreasing [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 *)
(******************************************************************************)

pi = N[Pi,40]; 

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]]]]}];



(******************************************************************************)
(* Groups *)
(******************************************************************************)

SetGroupTimelevels[g_,tl_] = Join[g, {Timelevels -> tl}];

evolvedGroups =
  {SetGroupName [CreateGroupFromTensor[phi          ], prefix <> "log_confac" ],
   SetGroupName [CreateGroupFromTensor[dphi[la]     ], prefix <> "dlog_confac"],
   SetGroupName [CreateGroupFromTensor[gt[la,lb]    ], prefix <> "metric"     ],
   SetGroupName [CreateGroupFromTensor[dgt[la,lb,lc]], prefix <> "dmetric"    ],
   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[dalpha[la]   ], prefix <> "dlapse"     ],
   SetGroupName [CreateGroupFromTensor[beta[ua]     ], prefix <> "shift"      ],
   SetGroupName [CreateGroupFromTensor[dbeta[la,ub] ], prefix <> "dshift"     ],
   SetGroupName [CreateGroupFromTensor[B[ua]        ], prefix <> "dtshift"    ]};
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"     ],
   SetGroupName [CreateGroupFromTensor[cdalpha[la]   ], prefix <> "cons_dlapse"     ],
   SetGroupName [CreateGroupFromTensor[cdbeta[la,ub] ], prefix <> "cons_dshift"     ],
   SetGroupName [CreateGroupFromTensor[cdphi[la]     ], prefix <> "cons_dlog_confac"],
   SetGroupName [CreateGroupFromTensor[cdgt[la,lb,lc]], prefix <> "cons_dmetric"    ]
};

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];



(******************************************************************************)
(* Gauge conditions *)
(******************************************************************************)

gaugef[alpha_,phi_] := harmonicF alpha^(harmonicN-2);
gaugeG[alpha_,phi_] := ShiftGammaCoeff alpha^(-2);
gaugeH[alpha_,phi_] := e4phi;



(******************************************************************************)
(* Initial data *)
(******************************************************************************)

initialCalc =
{
  Name -> BSSN <> "_Minkowski",
  Schedule -> {"IN ADMBase_InitialData"},
  ConditionalOnKeyword -> {"my_initial_data", "Minkowski"},
  Equations ->
  {
    phi           -> 0,
    gt[la,lb]     -> KD[la,lb],
    trK           -> 0,
    At[la,lb]     -> 0,
    Xt[ua]        -> 0,
    dphi[la]      -> 0,
    dgt[lc,la,lb] -> 0,
    alpha         -> 1,
    dalpha[la]    -> 0,
    beta[ua]      -> 0,
    dbeta[la,ub]  -> 0,
    B[ua]         -> 0
  }
};

(******************************************************************************)
(* Convert from ADMBase *)
(******************************************************************************)

convertFromADMBaseCalc =
{
  Name -> BSSN <> "_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       -> Log[detg]/12,
    em4phi    -> 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]
  }
};

convertFromADMBaseGammaCalc =
{
  Name -> BSSN <> "_convertFromADMBaseGamma",
  Schedule -> {"AT initial AFTER " <> BSSN <> "_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]},
  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]),
    
    dphi[la] -> 12 PD[phi,la],
    
    dgt[lc,la,lb] -> PD[gt[la,lb],lc],
    
    Xt[ua] -> gtu[ub,uc] Gt[ua,lb,lc],
    
    dalpha[la] -> PD[alpha,la],

    dbeta[la,ub] -> PD[beta[ub],la],
    
    B[ua] -> + admdtbeta[ua]
             - Upwind[beta[ub], admbeta[ua], lb]
  }
};



(******************************************************************************)
(* Convert to ADMBase *)
(******************************************************************************)

convertToADMBaseCalc =
{
  Name -> BSSN <> "_convertToADMBase",
  Schedule -> {"IN " <> BSSN <> "_convertToADMBaseGroup"},
  Where -> Everywhere,
  Shorthands -> {e4phi},
  Equations ->
  {
    e4phi       -> 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 -> BSSN <> "_convertToADMBaseDtLapseShift",
  Schedule -> {"IN " <> BSSN <> "_convertToADMBaseGroup"},
  Where -> Interior,
  Shorthands -> {dir[ua]},
  Equations ->
  {
    dir[ua] -> Sign[beta[ua]],
    
    (* see RHS *)
    admdtalpha -> - alpha^2 gaugef[alpha,phi] trK
                  + Upwind[beta[ua], alpha, la],
    admdtbeta[ua] -> + alpha^2 gaugeG[alpha,phi] B[ua]
                     + Upwind[beta[ub], beta[ua], lb]
  }
};

convertToADMBaseDtLapseShiftBoundaryCalc =
{
  Name -> BSSN <> "_convertToADMBaseDtLapseShiftBoundary",
  Schedule -> {"IN " <> BSSN <> "_convertToADMBaseGroup"},
  Where -> BoundaryWithGhosts,
  Equations ->
  {
    (* see RHS, but omit derivatives near the boundary *)
    admdtalpha -> - alpha^2 gaugef[alpha,phi] trK,
    admdtbeta[ua] -> + alpha^2 gaugeG[alpha,phi] B[ua]
  }
};

(******************************************************************************)
(* Evolution equations *)
(******************************************************************************)

evolCalc =
{
  Name -> BSSN <> "_RHS",
  Schedule -> {"IN " <> BSSN <> "_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, e4phi, em4phi,
                 gtu[ua,ub], Gtl[la,lb,lc], Gtlu[la,lb,uc], Gt[ua,lb,lc],
                 Xtn[ua], Rt[la,lb], Rphi[la,lb], R[la,lb],
                 Atm[ua,lb], Atu[ua,ub], Ats[la,lb], trAts,
                 dgts[lc,la,lb], trdgts[lc], trcdgt[lc],
                 g[la,lb], gu[ua,ub],
                 rho, S[la], trS},
  Equations ->
  {
    dir[ua] -> Sign[beta[ua]],
    
    detgt  -> 1 (* detgtExpr *),
    e4phi  -> Exp[4 phi],
    em4phi -> 1 / e4phi,
    
    (* This leads to simpler code... *)
    gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
    
    Gtl[la,lb,lc] -> 1/2 (dgt[lc,la,lb] + dgt[lb,la,lc] - dgt[la,lb,lc]),

    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],
    
    Rt[la,lb] -> - 1/2 gtu[uc,ud] PD[dgt[ld,la,lb],lc]
                 + 1/2 gt[lc,la] PD[Xt[uc],lb]
                 + 1/2 gt[lc,lb] PD[Xt[uc],la]
                 + 1/2 Xtn[ud] Gtl[la,lb,ld]
                 + 1/2 Xtn[ud] Gtl[lb,la,ld]
                 + (+ Gt[uc,la,ld] Gtlu[lb,lc,ud]
                    + Gt[uc,lb,ld] Gtlu[la,lc,ud]
                    + Gt[uc,la,ld] Gtlu[lc,lb,ud]),
    
    Rphi[la,lb] -> - 1/6 CDt[dphi[la],lb]
                   - 1/6 gt[la,lb] gtu[uc,ud] CDt[dphi[lc],ld]
                   + 1/36 dphi[la] dphi[lb]
                   - 1/36 gt[la,lb] gtu[uc,ud] dphi[lc] dphi[ld],
    
    R[la,lb] -> Rt[la,lb] + Rphi[la,lb],

    Atm[ua,lb] -> gtu[ua,uc] At[lc,lb],
    Atu[ua,ub] -> Atm[ua,lc] gtu[ub,uc],
    
    g[la,lb]    -> e4phi gt[la,lb],
    gu[ua,ub]   -> em4phi gtu[ua,ub],

    (* 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 (gu[ui,uj] T[li,lj]),

    (* Constraints *)
    cdalpha[la]    -> dalpha[la] - PD[alpha,la],
    cdbeta[la,ub]  -> dbeta[la,ub] - PD[beta[ub],la],
    cdphi[la]      -> dphi[la] - 12 PD[phi,la],
    cdgt[lc,la,lb] -> dgt[lc,la,lb] - PD[gt[la,lb],lc],
    
    (* Evolution equations *)
    dot[phi] -> - 1/6 alpha trK
                + 1/6 dbeta[lc,uc]
                + Upwind[beta[ua], phi, la],
    
    dot[trK] -> - em4phi gtu[ua,ub] (+ CDt[dalpha[lb],la]
                                     + 1/6 dphi[la] dalpha[lb])
                + alpha (Atu[ua,ub] At[la,lb] + 1/3 trK^2)
                + Upwind[beta[ua], trK, la],
    (* Matter term *)
                + addMatter (4 pi alpha (rho + trS)),
    
    dot[dphi[la]] -> - 2 alpha PD[trK,la]
                     + 2 PD[dbeta[lb,ub],la]
                     - 2 dalpha[la] trK
                     + dbeta[la,ub] dphi[lb]
                     - DphiDriver cdphi[la]
                     + Upwind[beta[ub], dphi[la], lb],
    
    dot[gt[la,lb]] -> - 2 alpha At[la,lb]
                      + gt[la,lc] dbeta[lb,uc]
                      + gt[lb,lc] dbeta[la,uc]
                      - 2/3 gt[la,lb] dbeta[lc,uc]
                      + Upwind[beta[uc], gt[la,lb], lc],
    
    Ats[la,lb] -> + alpha R[la,lb]
                  - CDt[dalpha[lb],la]
                  + 1/6 (dphi[la] dalpha[lb] + dphi[lb] dalpha[la]),
    
    trAts -> gu[ua,ub] Ats[la,lb],
    
    dot[At[la,lb]] -> + em4phi (Ats[la,lb] - 1/3 g[la,lb] trAts)
                      + alpha (trK At[la,lb] - 2 At[la,lc] Atm[uc,lb])
                      + At[la,lc] dbeta[lb,uc] + At[lb,lc] dbeta[la,uc]
                      - 2/3 At[la,lb] dbeta[lc,uc]
                      + Upwind[beta[uc], At[la,lb], lc],
    (* Matter term *)
                      + addMatter (- em4phi alpha 8 pi
                                     (T[la,lb] - (1/3) g[la,lb] trS)),
    
    dgts[lc,la,lb] -> + PD[dbeta[la,ud],lc] gt[lb,ld]
                      + PD[dbeta[lb,ud],lc] gt[la,ld],
    
    trdgts[lc] -> gtu[ua,ub] dgts[lc,la,lb],
    
    trcdgt[lc] -> gtu[ua,ub] cdgt[lc,la,lb],
    
    dot[dgt[lc,la,lb]] -> - 2 alpha PD[At[la,lb],lc]
                          + (dgts[lc,la,lb] - 1/3 gt[la,lb] trdgts[lc])
                          - 2 dalpha[lc] At[la,lb]
                          + dbeta[lc,ud] dgt[ld,la,lb]
                          + dgt[lc,ld,la] dbeta[lb,ud]
                          + dgt[lc,ld,lb] dbeta[la,ud]
                          - 2/3 dgt[lc,la,lb] dbeta[ld,ud]
                          - DgtDriver (+ cdgt[lc,la,lb]
                                       - 1/3 gt[la,lb] trcdgt[lc])
                          + Upwind[beta[ud], dgt[lc,la,lb], ld],
    
    dot[Xt[ua]] -> + gtu[uc,ud] PD[dbeta[ld,ua],lc]
                   + 1/3 gtu[ua,ub] PD[dbeta[lc,uc],lb]
                   + sigma gtu[ua,ub] (+ PD[dbeta[lc,uc],lb]
                                       - PD[dbeta[lb,uc],lc])
                   - 4/3 alpha gtu[ua,ub] PD[trK,lb]
                   - Xtn[ub] dbeta[lb,ua]
                   + 2/3 Xtn[ua] dbeta[lb,ub]
                   - 2 Atu[ua,ub] dalpha[lb]
                   + 2 alpha (+ Gt[ua,lc,ld] Atu[uc,ud]
                              + 1/2 Atu[ua,ub] dphi[lb])
                   + Upwind[beta[ub], Xt[ua], lb],
    (* Matter tetm *)
                   + addMatter (- 16 pi alpha gtu[ua,ub] S[lb]),
    
    dot[alpha] -> - alpha^2 gaugef[alpha,phi] trK
                  + Upwind[beta[ua], alpha, la],
    
    dot[dalpha[la]] -> - alpha^2 gaugef[alpha,phi] PD[trK,la]
                       - 2 gaugef[alpha,phi] alpha dalpha[la] trK
                       - alpha^2 D[gaugef[alpha,phi],alpha] dalpha[la] trK
                       - 1/12 alpha^2 D[gaugef[alpha,phi],phi] dphi[la] trK
                       + dbeta[la,ub] dalpha[lb]
                       - DAlphaDriver cdalpha[la]
                       + Updinw[beta[ub], dalpha[la], lb],
    
    dot[beta[ua]] -> + alpha^2 gaugeG[alpha,phi] B[ua]
                     + Upwind[beta[ub], beta[ua], lb],
    
    dot[B[ua]] -> + em4phi gaugeH[alpha,phi] dot[Xt[ua]]
                  - BetaDriver B[ua]
                  + Upwind[beta[ub], B[ua], lb]
                  - Upwind[beta[ub], Xt[ua], lb],
    
    dot[dbeta[la,ub]] -> + alpha^2 gaugeG[alpha,phi] PD[B[ub],la]
                         + 2 gaugeG[alpha,phi] alpha dalpha[la] B[ub]
                         + alpha^2 D[gaugeG[alpha,phi],alpha] dalpha[la] B[ub]
                         + 1/12 alpha^2 D[gaugeG[alpha,phi],phi] dphi[la] B[ub]
                         + dbeta[la,uc] dbeta[lc,ub]
                         - DBetaDriver cdbeta[la,ub]
                         + Upwind[beta[uc], dbeta[la,ub], lc]
  }
};

evol1Calc =
{
  Name -> BSSN <> "_RHS1",
  Schedule -> {"IN " <> BSSN <> "_evolCalcGroup BEFORE " <> BSSN <> "_RHS3"},
  (*
  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, e4phi, em4phi,
                 gtu[ua,ub], Gtl[la,lb,lc], Gt[ua,lb,lc],
                 Xtn[ua], Rt[la,lb], Rphi[la,lb], R[la,lb],
                 Atm[ua,lb], Atu[ua,ub], Ats[la,lb],
                 dgts[lc,la,lb], trdgts[lc], trcdgt[lc],
                 S[la]},
  Equations ->
  {
    dir[ua] -> Sign[beta[ua]],
    
    detgt  -> 1 (* detgtExpr *),
    e4phi  -> Exp[4 phi],
    em4phi -> 1 / e4phi,
    
    (* This leads to simpler code... *)
    gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
    
(*    Gtl[la,lb,lc] -> 1/2 (dgt[lc,la,lb] + dgt[lb,la,lc] - dgt[la,lb,lc]),

    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], *)

    Gt[ua,lb,lc] -> 1/2 gtu[ua,ud] 
                       (dgt[lc,ld,lb] + dgt[lb,ld,lc] - dgt[ld,lb,lc]),

    Xtn[ua] -> gtu[ub,uc] Gt[ua,lb,lc], 
    
    Atm[ua,lb] -> gtu[ua,uc] At[lc,lb],
    Atu[ua,ub] -> Atm[ua,lc] gtu[ub,uc],
    
    (* Matter terms *)

    (* 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])),

    (* Constraints *)
    cdgt[lc,la,lb] -> dgt[lc,la,lb] - PD[gt[la,lb],lc],
    
    (* Evolution equations *)
    dgts[lc,la,lb] -> + PD[dbeta[la,ud],lc] gt[lb,ld]
                      + PD[dbeta[lb,ud],lc] gt[la,ld],
    
    trdgts[lc] -> gtu[ua,ub] dgts[lc,la,lb],
    
    trcdgt[lc] -> gtu[ua,ub] cdgt[lc,la,lb],
    
    dot[dgt[lc,la,lb]] -> - 2 alpha PD[At[la,lb],lc]
                          + (dgts[lc,la,lb] - 1/3 gt[la,lb] trdgts[lc])
                          - 2 dalpha[lc] At[la,lb]
                          + dbeta[lc,ud] dgt[ld,la,lb]
                          + dgt[lc,ld,la] dbeta[lb,ud]
                          + dgt[lc,ld,lb] dbeta[la,ud]
                          - 2/3 dgt[lc,la,lb] dbeta[ld,ud]
                          - DgtDriver (+ cdgt[lc,la,lb]
                                       - 1/3 gt[la,lb] trcdgt[lc])
                          + Upwind[beta[ud], dgt[lc,la,lb], ld],
    
    dot[Xt[ua]] -> + gtu[uc,ud] PD[dbeta[ld,ua],lc]
                   + 1/3 gtu[ua,ub] PD[dbeta[lc,uc],lb]
                   + sigma gtu[ua,ub] (+ PD[dbeta[lc,uc],lb]
                                       - PD[dbeta[lb,uc],lc])
                   - 4/3 alpha gtu[ua,ub] PD[trK,lb]
                   - Xtn[ub] dbeta[lb,ua]
                   + 2/3 Xtn[ua] dbeta[lb,ub]
                   - 2 Atu[ua,ub] dalpha[lb]
                   + 2 alpha (+ Gt[ua,lc,ld] Atu[uc,ud]
                              + 1/2 Atu[ua,ub] dphi[lb])
                   + Upwind[beta[ub], Xt[ua], lb]
    (* Matter term *)
                   + addMatter (- 16 pi alpha gtu[ua,ub] S[lb])
  }
};

evol2Calc =
{
  Name -> BSSN <> "_RHS2",
  Schedule -> {"IN " <> BSSN <> "_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, e4phi, em4phi,
                 gtu[ua,ub], Gtl[la,lb,lc], Gtlu[la,lb,uc], Gt[ua,lb,lc],
                 Xtn[ua], Rt[la,lb], Rphi[la,lb], R[la,lb],
                 Atm[ua,lb], Atu[ua,ub], Ats[la,lb], trAts,
                 dgts[lc,la,lb], trdgts[lc], trcdgt[lc],
                 g[la,lb], gu[ua,ub],
                 rho, trS},
  Equations ->
  {
    dir[ua] -> Sign[beta[ua]],
    
    detgt  -> 1 (* detgtExpr *),
    e4phi  -> Exp[4 phi],
    em4phi -> 1 / e4phi,
    
    (* This leads to simpler code... *)
    gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],
    
    Gtl[la,lb,lc] -> 1/2 (dgt[lc,la,lb] + dgt[lb,la,lc] - dgt[la,lb,lc]),

    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],
    
    Rt[la,lb] -> - 1/2 gtu[uc,ud] PD[dgt[ld,la,lb],lc]
                 + 1/2 gt[lc,la] PD[Xt[uc],lb]
                 + 1/2 gt[lc,lb] PD[Xt[uc],la]
                 + 1/2 Xtn[ud] Gtl[la,lb,ld]
                 + 1/2 Xtn[ud] Gtl[lb,la,ld]
                 + (+ Gt[uc,la,ld] Gtlu[lb,lc,ud]
                    + Gt[uc,lb,ld] Gtlu[la,lc,ud]
                    + Gt[uc,la,ld] Gtlu[lc,lb,ud]),
    
    Rphi[la,lb] -> - 1/6 CDt[dphi[la],lb]
                   - 1/6 gt[la,lb] gtu[uc,ud] CDt[dphi[lc],ld]
                   + 1/36 dphi[la] dphi[lb]
                   - 1/36 gt[la,lb] gtu[uc,ud] dphi[lc] dphi[ld],
    
    R[la,lb] -> Rt[la,lb] + Rphi[la,lb],

    Atm[ua,lb] -> gtu[ua,uc] At[lc,lb],
    Atu[ua,ub] -> Atm[ua,lc] gtu[ub,uc],
    
    g[la,lb]    -> e4phi gt[la,lb],
    gu[ua,ub]   -> em4phi gtu[ua,ub],

    (* 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])),

    (* trS = gamma^ij T_ij  *)
    trS -> addMatter (gu[ui,uj] T[li,lj]),

    (* Evolution equations *)
    
    dot[trK] -> - em4phi gtu[ua,ub] (+ CDt[dalpha[lb],la]
                                     + 1/6 dphi[la] dalpha[lb])
                + alpha (Atu[ua,ub] At[la,lb] + 1/3 trK^2)
                + Upwind[beta[ua], trK, la]
    (* Matter term *)
                + addMatter (4 pi alpha (rho + trS)),
    
    Ats[la,lb] -> + alpha R[la,lb]
                  - CDt[dalpha[lb],la]
                  + 1/6 (dphi[la] dalpha[lb] + dphi[lb] dalpha[la]),
    
    trAts -> gu[ua,ub] Ats[la,lb],
    
    dot[At[la,lb]] -> + em4phi (Ats[la,lb] - 1/3 g[la,lb] trAts)
                      + alpha (trK At[la,lb] - 2 At[la,lc] Atm[uc,lb])
                      + At[la,lc] dbeta[lb,uc] + At[lb,lc] dbeta[la,uc]
                      - 2/3 At[la,lb] dbeta[lc,uc]
                      + Upwind[beta[uc], At[la,lb], lc]
    (* Matter term *)
                      + addMatter (- em4phi alpha 8 pi
                                     (T[la,lb] - (1/3) g[la,lb] trS))
  }
};

evol3Calc =
{
  Name -> BSSN <> "_RHS3",
  Schedule -> {"IN " <> BSSN <> "_evolCalcGroup AFTER " <> BSSN <> "_RHS1"},
  (*
  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, e4phi, em4phi, gtu[ua,ub]},
  Equations ->
  {
    dir[ua] -> Sign[beta[ua]],
    
    detgt  -> 1 (* detgtExpr *),
    e4phi  -> Exp[4 phi],
    em4phi -> 1 / e4phi,

    (* This leads to simpler code... *)
    gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],

    (* Constraints *)
    cdalpha[la]    -> dalpha[la] - PD[alpha,la],
    cdbeta[la,ub]  -> dbeta[la,ub] - PD[beta[ub],la],
    cdphi[la]      -> dphi[la] - 12 PD[phi,la],

    (* Evolution equations *)

    dot[phi] -> - 1/6 alpha trK
                + 1/6 dbeta[lc,uc]
                + Upwind[beta[ua], phi, la],
    
    dot[dphi[la]] -> - 2 alpha PD[trK,la]
                     + 2 PD[dbeta[lb,ub],la]
                     - 2 dalpha[la] trK
                     + dbeta[la,ub] dphi[lb]
                     - DphiDriver cdphi[la]
                     + Upwind[beta[ub], dphi[la], lb],
    
    dot[gt[la,lb]] -> - 2 alpha At[la,lb]
                      + gt[la,lc] dbeta[lb,uc]
                      + gt[lb,lc] dbeta[la,uc]
                      - 2/3 gt[la,lb] dbeta[lc,uc]
                      + Upwind[beta[uc], gt[la,lb], lc],
    
    dot[alpha] -> - alpha^2 gaugef[alpha,phi] trK
                  + Upwind[beta[ua], alpha, la],
    
    dot[dalpha[la]] -> - alpha^2 gaugef[alpha,phi] PD[trK,la]
                       - 2 gaugef[alpha,phi] alpha dalpha[la] trK
                       - alpha^2 D[gaugef[alpha,phi],alpha] dalpha[la] trK
                       - 1/12 alpha^2 D[gaugef[alpha,phi],phi] dphi[la] trK
                       + dbeta[la,ub] dalpha[lb]
                       - DAlphaDriver cdalpha[la]
                       + Upwind[beta[ub], dalpha[la], lb],

    dot[beta[ua]] -> IfThen[harmonicShift,
                            - 1/2 gtu[ua,uj] em4phi alpha
                              (- 2 alpha (-2) PD[phi,lj]
                               + 2 PD[alpha,lj]
                               + alpha ((* Ddetgt[lj] *) - 2 gtu[uk,ul] PD[gt[lj,lk],ll])) + Upwind[beta[ub], beta[ua], lb],
                            (* else *)

                     + alpha^2 gaugeG[alpha,phi] B[ua]
                     + Upwind[beta[ub], beta[ua], lb]],
    
    dot[B[ua]] -> + em4phi gaugeH[alpha,phi] dot[Xt[ua]]
                  - BetaDriver B[ua]
                  + Upwind[beta[ub], B[ua], lb]
                  - Upwind[beta[ub], Xt[ua], lb],
    
    dot[dbeta[la,ub]] -> IfThen[harmonicShift,
                         Get["dot-dbeta.m"],
                         (* else *)
                         + alpha^2 gaugeG[alpha,phi] PD[B[ub],la]
                         + 2 gaugeG[alpha,phi] alpha dalpha[la] B[ub]
                         + alpha^2 D[gaugeG[alpha,phi],alpha] dalpha[la] B[ub]
                         + 1/12 alpha^2 D[gaugeG[alpha,phi],phi] dphi[la] B[ub]
                         + dbeta[la,uc] dbeta[lc,ub]
                         - DBetaDriver cdbeta[la,ub]
                         + Upwind[beta[uc], dbeta[la,ub], lc]]
  }
};

RHSStaticBoundaryCalc =
{
  Name -> BSSN <> "_RHSStaticBoundary",
  Schedule -> {"IN MoL_CalcRHS"},
  ConditionalOnKeyword -> {"my_rhs_boundary_condition", "static"},
  Where -> Boundary,
  Equations ->
  {
    phi           -> 0,
    gt[la,lb]     -> KD[la,lb],
    trK           -> 0,
    At[la,lb]     -> 0,
    Xt[ua]        -> 0,
    dphi[la]      -> 0,
    dgt[lc,la,lb] -> 0,
    alpha         -> 1,
    dalpha[la]    -> 0,
    beta[ua]      -> 0,
    dbeta[la,ub]  -> 0,
    B[ua]         -> 0
  }
};

enforceCalc =
{
  Name -> BSSN <> "_enforce",
  Schedule -> {"IN MoL_PostStepModify"},
  Shorthands -> {detgt, gtu[ua,ub], trAt, trdgt[la]},
  Equations ->
  {
    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,
    
    trdgt[lc] -> gtu[ua,ub] dgt[lc,la,lb],

    dgt[lc,la,lb] -> dgt[lc,la,lb] - (1/3) gt[la,lb] trdgt[lc],

    alpha -> Max[alpha, MinimumLapse]
  }
};

(******************************************************************************)
(* Boundary conditions *)
(******************************************************************************)

boundaryCalc =
{
  Name -> BSSN <> "_boundary",
  Schedule -> {"IN MoL_PostStep"},
  ConditionalOnKeyword -> {"my_boundary_condition", "Minkowski"},
  Where -> BoundaryWithGhosts,
  Equations ->
  {
    phi           -> 0,
    gt[la,lb]     -> KD[la,lb],
    trK           -> 0,
    At[la,lb]     -> 0,
    Xt[ua]        -> 0,
    dphi[la]      -> 0,
    dgt[lc,la,lb] -> 0,
    alpha         -> 1,
    dalpha[la]    -> 0,
    beta[ua]      -> 0,
    dbeta[la,ub]  -> 0,
    B[ua]         -> 0
  }
};

(******************************************************************************)
(* Constraint equations *)
(******************************************************************************)

constraintsCalc =
{
  Name -> BSSN <> "_constraints",
  Schedule -> Automatic,
  After -> "MoL_PostStep",
  Where -> Interior,
  Shorthands -> {detgt, gtu[ua,ub], Gtl[la,lb,lc], Gtlu[la,lb,uc],
                 Gt[ua,lb,lc], 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],
                 rho, S[la]},
  Equations ->
  {
    detgt        -> 1 (* detgtExpr *),
    gtu[ua,ub]   -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]],

    Gtl[la,lb,lc] -> 1/2 (dgt[lc,la,lb] + dgt[lb,la,lc] - dgt[la,lb,lc]),

    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],

    Rt[la,lb] -> - 1/2 gtu[uc,ud] PD[dgt[ld,la,lb],lc]
                 + 1/2 gt[lc,la] PD[Xt[uc],lb]
                 + 1/2 gt[lc,lb] PD[Xt[uc],la]
                 + 1/2 Xtn[ud] Gtl[la,lb,ld]
                 + 1/2 Xtn[ud] Gtl[lb,la,ld]
                 + (+ Gt[uc,la,ld] Gtlu[lb,lc,ud]
                    + Gt[uc,lb,ld] Gtlu[la,lc,ud]
                    + Gt[uc,la,ld] Gtlu[lc,lb,ud]),

    Rphi[la,lb] -> - 1/6 CDt[dphi[la],lb]
                   - 1/6 gt[la,lb] gtu[uc,ud] CDt[dphi[lc],ld]
                   + 1/36 dphi[la] dphi[lb]
                   - 1/36 gt[la,lb] gtu[uc,ud] dphi[lc] dphi[ld],

    
    e4phi       -> Exp[4 phi],
    em4phi      -> 1 / e4phi,
    
    R[la,lb] -> Rt[la,lb] + Rphi[la,lb],
    trR      -> em4phi gtu[ua,ub] R[la,lb],
    
    Atm[ua,lb] -> gtu[ua,uc] At[lc,lb],
    
    (* Matter terms (Commented out for now since we don't have matter terms
       in the evolution equations either *)
    
    (* 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 - Atm[ua,lb] Atm[ub,la] + (2/3) trK^2 - addMatter 16 pi rho,
    
    M[la] -> gtu[ub,uc] ( CDt[At[lc,la],lb] + (1/2) At[lb,la] dphi[lc] )
             - (2/3) PD[trK,la] - addMatter 8 pi S[la],
    
    (* 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]
  }
};

constraintsBoundaryCalc =
{
  Name -> BSSN <> "_constraints_boundary",
  Schedule -> Automatic,
  After -> "MoL_PostStep",
  Where -> BoundaryWithGhosts,
  Equations ->
  {
    H              -> 0,
    M[la]          -> 0,
    cS             -> 0,
    cXt[ua]        -> 0,
    cA             -> 0,
    cdalpha[la]    -> 0,
    cdbeta[la,ub]  -> 0,
    cdphi[la]      -> 0,
    cdgt[lc,la,lb] -> 0
  }
};

(******************************************************************************)
(* Implementations *)
(******************************************************************************)

inheritedImplementations =
  Join[{"ADMBase"},
       If [addMatter!=0, {"TmunuBase"}, {}]];

(******************************************************************************)
(* Parameters *)
(******************************************************************************)

inheritedKeywordParameters = {};

extendedKeywordParameters =
{
  {
    Name -> "ADMBase::evolution_method",
    AllowedValues -> {BSSN}
  },
  {
    Name -> "ADMBase::lapse_evolution_method",
    AllowedValues -> {BSSN}
  },
  {
    Name -> "ADMBase::shift_evolution_method",
    AllowedValues -> {BSSN}
  },
  {
    Name -> "ADMBase::dtlapse_evolution_method",
    AllowedValues -> {BSSN}
  },
  {
    Name -> "ADMBase::dtshift_evolution_method",
    AllowedValues -> {BSSN}
  }
};

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"},
    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"
  }
};

intParameters =
{
  {
    Name -> harmonicN,
    Description -> "d/dt alpha = - f alpha^n K  (harmonic=2, 1+log=1)",
    Default -> 2
  },
  {
    Name -> ShiftAlphaPower,
    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 =
{
  {
    Name -> harmonicF,
    Description -> "d/dt alpha = - f alpha^n K   (harmonic=1, 1+log=2)",
    Default -> 1
  },
  {
    Name -> DAlphaDriver,
    Default -> 0
  },
  {
    Name -> ShiftGammaCoeff,
    Default -> 0
  },
  {
    Name -> BetaDriver,
    Default -> 0
  },
  {
    Name -> DBetaDriver,
    Default -> 0
  },
  {
    Name -> DphiDriver,
    Default -> 0
  },
  {
    Name -> DgtDriver,
    Default -> 0
  },
  {
    Name -> sigma,
    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
  }
};

(******************************************************************************)
(* Construct the thorns *)
(******************************************************************************)

calculations =
{
  initialCalc,
  convertFromADMBaseCalc,
  convertFromADMBaseGammaCalc,
(*  evolCalc, *)
  evol1Calc, evol2Calc, evol3Calc,
  RHSStaticBoundaryCalc,
  enforceCalc,
  boundaryCalc,
  convertToADMBaseCalc,
  convertToADMBaseDtLapseShiftCalc,
  convertToADMBaseDtLapseShiftBoundaryCalc,
  constraintsCalc,
  constraintsBoundaryCalc
};

CreateKrancThornTT [groups, ".", BSSN,
  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
];

];



(******************************************************************************)
(* Options *)
(******************************************************************************)

(* 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) *)

createCode[4, False, True, 3, 1];
