
SetEnhancedTimes[False];

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

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

prefix = "CT_";
suffix =
  ""
  <> If [derivOrder!=4, "_O" <> ToString[derivOrder], ""]
  ;

CTThorn = prefix <> "Dust" <> 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;

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,
     {DD, EE, SS, 
      rho, eps, prs, V, W, u, T, T0, t, t0, cttS, cttV, rhodp,
      g, gu, detg, detg4, gt, gtu, Gt, Xt, k, a, b, dir}
]; 

Map [AssertSymmetricIncreasing,
     {g[la,lb], gt[la,lb], G[ua,lb,lc], T[la,lb], k[la,lb]}];
Map [AssertSymmetricIncreasing,
     {gu[ua,ub], gtu[ua,ub]}];

(* Use the CartGrid3D variable names *)
x1=x; x2=y; x3=z;

(* Use the ADMBase variable names *)
g11=gxx; g12=gxy; g22=gyy; g13=gxz; g23=gyz; g33=gzz;
a=alp;
b1=betax; b2=betay; b3=betaz;
k11=kxx; k12=kxy; k22=kyy; k13=kxz; k23=kyz; k33=kzz;
bssnphi=phi; bssng11=gt11; bssng12=gt12; bssng13=gt13; bssng22=gt22; bssng23=gt23; bssng33=gt33;
bssntrK=trK; bssnA11=At11; bssnA12=At12; bssnA13=At13; bssnA22=At22; bssnA23=At23; bssnA33=At33;

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

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

evolvedGroups =
  {SetGroupName [CreateGroupFromTensor [DD      ], prefix <> "D"],
   SetGroupName [CreateGroupFromTensor [EE      ], prefix <> "E"],
   SetGroupName [CreateGroupFromTensor [SS[la]  ], prefix <> "S"]};
evaluatedGroups =
  {SetGroupName [CreateGroupFromTensor [rho    ], prefix <> "rho"],
   SetGroupName [CreateGroupFromTensor [eps    ], prefix <> "eps"],
   SetGroupName [CreateGroupFromTensor [prs    ], prefix <> "prs"],
   SetGroupName [CreateGroupFromTensor [u[la]  ], prefix <> "u"],
   SetGroupName [CreateGroupFromTensor [V[ua]  ], prefix <> "V"],
   SetGroupName [CreateGroupFromTensor [W      ], prefix <> "W"],
   SetGroupName [CreateGroupFromTensor [cttS   ], prefix <> "cttS"],
   SetGroupName [CreateGroupFromTensor [cttV[la]], prefix <> "cttV"],
   SetGroupName [CreateGroupFromTensor [rhodp  ], prefix <> "rhodp"]};

declaredGroups = Join [evolvedGroups, evaluatedGroups];
declaredGroupNames = Map [First, declaredGroups];



extraGroups =
  {{"Grid::coordinates", {x, y, z, r}},
   {"ADMBase::metric",  {gxx, gxy, gxz, gyy, gyz, gzz}},
   {"ADMBase::lapse",   {alp}},
   {"ADMBase::shift",   {betax, betay, betaz}},
   {"ADMBase::curv",  {kxx, kxy, kxz, kyy, kyz, kzz}},
   {"ML_BSSN::ML_log_confac",  {phi}},
   {"ML_BSSN::ML_metric",   {gt11, gt12, gt13, gt22, gt23, gt33}},
   {"ML_BSSN::ML_trace_curv",   {trK}},
   {"ML_BSSN::ML_curv",   {At11, At12, At13, At22, At23, At33}},
   {"ML_BSSN::ML_gamma",   {Xt1, Xt2, Xt3}},
   {"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 *)
(******************************************************************************)

initialMinkowskiCalc =
{
  Name -> CTThorn <> "_Minkowski",
  Schedule -> {"IN CCTK_INITIAL after ADMBase_Initial before " <> CTThorn <> "_setCTTrhs"},
  ConditionalOnKeyword -> {"my_initial_data", "Minkowski"},
  Equations -> 
  {
    rho        -> 0,
    rhodp      -> 0,
    eps        -> 0,
    prs        -> w rho,
    u[la]      -> 0
  }
};

initialFLRWPertCalc =
{
  Name -> CTThorn <> "_FLRW_Pert",
  Schedule -> {"IN CCTK_INITIAL after ADMBase_Initial before " <> CTThorn <> "_setCTTrhs"},
  ConditionalOnKeyword -> {"my_initial_data", "pFLRW"},
  Shorthands -> {rhocrit, a0t, a2, H0t, trK2, aa, HH},
  Equations -> 
  {
    HH         -> "pflrw_H0" "pflrw_t0" / t,
    aa         -> "pflrw_a0" ( t / "pflrw_t0" )^(2/3), (*( 1 + (x^2+y^2+z^2) / 4 )^(-2),*)

    rhocrit    -> 3 HH^2 / (8 Pi),
 
    rho        -> "pflrw_omegaM" rhocrit ( 1 + Sum[ "(pflrw_ax["<>ToString[i]<>"])" Sin["pflrw_kx["<>ToString[i]<>"]" x], {i,0,19}]
			        	     + Sum[ "(pflrw_ay["<>ToString[i]<>"])" Sin["pflrw_ky["<>ToString[i]<>"]" y], {i,0,19}]
                                             + Sum[ "(pflrw_az["<>ToString[i]<>"])" Sin["pflrw_kz["<>ToString[i]<>"]" z], {i,0,19}] ),
    rhodp      -> rho aa^3,
    eps        -> 0,
    prs        -> w rho,
    u[la]      -> 0,

    a2         -> aa^2,

    g11        -> a2,
    g12        -> 0,
    g13        -> 0,
    g22        -> a2,
    g23        -> 0,
    g33        -> a2,

    k11        -> - a2 HH,
    k12        -> 0,
    k13        -> 0,
    k22        -> - a2 HH,
    k23        -> 0,
    k33        -> - a2 HH
  }
};

(* Szekeres-Meures-Bruni model from astro-ph:1103.0501 *)
Hub0 = "pflrw_H0";                   (* Initial Hubble *)
OL = Lambda / (3 Hub0^2);   (* Omega_Lambda   *)
OM = 1 - OL;                (* Omega_M        *)
BB = 3/4 Hub0^2 (OL (1 - OL)^2)^(1/3);
kk = 1;
r0b = (3 Hub0^2 OM)/(8 \[Pi]);
BBp[z_] := 1 - Sin[kk z];
BBm[z_] := 0;
Gam[z_] := 0;
Omega[z_] := 0;
Tau[tt_] := Sqrt[(3 Lambda)/4] tt;
fm[tt_] := Cosh[Tau[tt]]/Sinh[Tau[tt]];
FSolve = DSolve[
   D[FF[x], {x, 2}] + 4/3 Coth[x] D[FF[x], x] - 
     2/3 Sinh[x]^-2 FF[x] == 0, FF[x], x];
fp[tt_] := 
 FF[x] /. FSolve[[1, 1]] /. {x -> Tau[tt], C[1] -> 0, C[2] -> 1};
S[tt_] := ((1 - OL)/OL)^(1/3)
   Sinh[3/2 Hub0^2 Sqrt[OL] tt]^(2/3);
dS[tt_] := D[S[ttt], ttt] /. ttt->tt;
A[x_, y_, z_, tt_] := 
  1 + BB BBp[z] ((x + Gam[z])^2 + (y + Omega[z])^2);
F[z_, tt_] := fp[tt] BBp[z] + fm[tt] BBm[z];
Z[x_, y_, z_, tt_] := F[z, tt] + A[x, y, z, tt];
dZ[x_, y_, z_, tt_] := D[Z[x, y, z, ttt], ttt] /.ttt->tt;
dens[x_, y_, z_, tt_] := r0b A[x, y, z, tt]/(
   S[tt]^3 Z[x, y, z, tt]);
tini=ArcSinh[((1 - OL)/OL)^(-1/2)]/(3/2 Hub0 Sqrt[OL]);

cMet[x_, y_, z_, t_] = 
  { {S[t]^2, 0, 0}, 
    {0, S[t]^2, 0}, 
    {0, 0, S[t]^2 Z[x,y,z,t]^2} };
cInvMet[x_, y_, z_, t_] = Simplify[Inverse[cMet[x,y,z,t]],
   Element[x, Reals] && Element[y, Reals] && Element[z, Reals] && 
    Element[t, Reals]];
dcMet[x_, y_, z_, t_] = Det[cMet[x,y,z,t]];
BSSNGamma[n_] := Sum[
    - 1/3 dcMet[x,y,z,t]^(-2/3) cInvMet[x,y,z,t][[n, j]] D[dcMet[x,y,z,t], {x,y,z}[[j]]] -
    dcMet[x,y,z,t]^(1/3) D[cInvMet[x,y,z,t][[n, j]], {x, y, z}[[j]]], {j, 1, 3}];
G1[x_, y_, z_, t_] = BSSNGamma[1];
G2[x_, y_, z_, t_] = BSSNGamma[2];
G3[x_, y_, z_, t_] = BSSNGamma[3];

initialMBCalc =
{
  Name -> CTThorn <> "_MB",
  Schedule -> {"IN CCTK_INITIAL after ADMBase_Initial before " <> CTThorn <> "_setCTTrhs"},
  ConditionalOnKeyword -> {"my_initial_data", "SMB"},
  Shorthands -> {detg},
  Equations -> 
  {
    g11        -> S[t]^2,
    g12        -> 0,
    g13        -> 0,
    g22        -> S[t]^2,
    g23        -> 0,
    g33        -> S[t]^2 Z[x,y,z,t]^2,
    detg       -> Det[{{g11,g12,g13},{g12,g22,g23},{g13,g23,g33}}],

    k11        -> - dS[t] S[t],
    k12        -> 0,
    k13        -> 0,
    k22        -> - dS[t] S[t],
    k23        -> 0,
    k33        -> - Z[x,y,z,t] S[t] ( dZ[x,y,z,t] S[t] + Z[x,y,z,t] dS[t] ),

    a          -> 1,
    b[ua]      -> 0,

    rho        -> dens[x,y,z,t],
    eps        -> 0,
    prs        -> w rho,
    u[la]      -> 0,
    rhodp      -> rho Sqrt[detg]
  }
};

boundaryMBCalc =
{
  Name -> CTThorn <> "_MB_bound",
  Schedule -> {"IN MoL_PostStep"},
  ConditionalOnKeyword -> {"my_initial_data", "SMB"},
  Where -> BoundaryWithGhosts,
  Shorthands -> {detg, det4g, gu[ua,ub], gtu[ua,ub], Gt[ua,lb,lc], f1, f2, u0, h},
  Equations -> 
  {
    g11        -> S[t]^2,
    g12        -> 0,
    g13        -> 0,
    g22        -> S[t]^2,
    g23        -> 0,
    g33        -> S[t]^2 Z[x,y,z,t]^2,
    detg       -> Det[{{g11,g12,g13},{g12,g22,g23},{g13,g23,g33}}],

    k11        -> - dS[t] S[t],
    k12        -> 0,
    k13        -> 0,
    k22        -> - dS[t] S[t],
    k23        -> 0,
    k33        -> - Z[x,y,z,t] S[t] ( dZ[x,y,z,t] S[t] + Z[x,y,z,t] dS[t] ),

    det4g      -> - a detg,
    gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]],
    f1         -> b[ua] u[la] / a^2,
    f2         -> (gu[ua,ub] u[la] u[lb] + 1) / a^2,
    u0         -> (f1 + Sqrt[f1^2 + 4 f2]) / 2,

    h          -> 1 + eps + w,
    W          -> Sqrt[-det4g] u0,
    rho        -> dens[x,y,z,t],
    eps        -> 0,
    prs        -> w rho,
    u[la]      -> 0,
    rhodp      -> rho Sqrt[detg],

    phi        -> detg^(-1/6),
    gt11       -> phi^2 g11,
    gt12       -> phi^2 g12,
    gt13       -> phi^2 g13,
    gt22       -> phi^2 g22,
    gt23       -> phi^2 g23,
    gt33       -> phi^2 g33,

    trK        -> gu[ua,ub] k[la,lb],
    At11       -> phi^2 ( k11 - trK g11 / 3),
    At12       -> phi^2 ( k12 - trK g12 / 3),
    At13       -> phi^2 ( k13 - trK g13 / 3),
    At22       -> phi^2 ( k22 - trK g22 / 3),
    At23       -> phi^2 ( k23 - trK g23 / 3),
    At33       -> phi^2 ( k33 - trK g33 / 3),

    Xt1        -> 0, (* Simplified BCs to speed up compilation -- correct BCs on gamma branch *)
    Xt2        -> 0,
    Xt3        -> 0
  }
};

(* LTB model from gr-qc:9812046, section 5.4 *)
r[x_, y_, z_] := Sqrt[x^2 + y^2 + z^2];
th[x_, y_, z_] := ArcTan[Sqrt[x^2 + y^2]/Sqrt[z^2]];
ph[x_, y_, z_] := ArcTan[Sqrt[y^2]/Sqrt[x^2]];
L[x_, y_, z_] := {
   {D[r[x, y, z], x], D[th[x, y, z], x], D[ph[x, y, z], x]},
   {D[r[x, y, z], y], D[th[x, y, z], y], D[ph[x, y, z], y]},
   {D[r[x, y, z], z], D[th[x, y, z], z], D[ph[x, y, z], z]}
   };
g[X_, Y_] := {
   {X^2, 0, 0},
   {0, Y^2, 0},
   {0, 0, Y^2 Sin[th[x, y, z]]^2}
   };
cartMet[X_, Y_, x_, y_, z_] = 
  L[x, y, z].g[X, Y].Transpose[L[x, y, z]] // Simplify;
M[r_] := 4/3 1/6 r^3;
a[r_] := 5 - 1/(1 + r^2/"ltb_rScal");
Y[r_, t_] := (9/2 M[r])^(1/3) (a[r] - t)^(2/3);
X[r_, t_] = 
  D[Y[r, t], r];
\[Rho][r_, t_] = 
  D[M[r], r]/(4 \[Pi] D[Y[r, t], r] Y[r, t]^2) // Simplify;
cMet[x_, y_, z_, t_] = 
  Simplify[cartMet[X[r, t], Y[r, t], x, y, z] /. r -> r[x, y, z], 
   Element[x, Reals] && Element[y, Reals] && Element[z, Reals] && 
    Element[t, Reals]];
cInvMet[x_, y_, z_, t_] = Simplify[Inverse[cMet[x,y,z,t]],
   Element[x, Reals] && Element[y, Reals] && Element[z, Reals] && 
    Element[t, Reals]];
dcMet[x_, y_, z_, t_] = Det[cMet[x,y,z,t]];
BSSNGamma[n_] := Sum[
    - 1/3 dcMet[x,y,z,t]^(-2/3) cInvMet[x,y,z,t][[n, j]] D[dcMet[x,y,z,t], {x,y,z}[[j]]] -
    dcMet[x,y,z,t]^(1/3) D[cInvMet[x,y,z,t][[n, j]], {x, y, z}[[j]]], {j, 1, 3}];
m11[x_, y_, z_, t_] = cMet[x, y, z, t][[1, 1]];
m12[x_, y_, z_, t_] = cMet[x, y, z, t][[1, 2]];
m13[x_, y_, z_, t_] = cMet[x, y, z, t][[1, 3]];
m22[x_, y_, z_, t_] = cMet[x, y, z, t][[2, 2]];
m23[x_, y_, z_, t_] = cMet[x, y, z, t][[2, 3]];
m33[x_, y_, z_, t_] = cMet[x, y, z, t][[3, 3]];
G1[x_, y_, z_, t_] = BSSNGamma[1];
G2[x_, y_, z_, t_] = BSSNGamma[2];
G3[x_, y_, z_, t_] = BSSNGamma[3];
c11[x_, y_, z_, 
   t_] = -1/2 Simplify[D[cMet[x, y, z, t][[1, 1]], t], t >= 1];
c12[x_, y_, z_, 
   t_] = -1/2 Simplify[D[cMet[x, y, z, t][[1, 2]], t], t >= 1];
c13[x_, y_, z_, 
   t_] = -1/2 Simplify[D[cMet[x, y, z, t][[1, 3]], t], t >= 1];
c22[x_, y_, z_, 
   t_] = -1/2 Simplify[D[cMet[x, y, z, t][[2, 2]], t], t >= 1];
c23[x_, y_, z_, 
   t_] = -1/2 Simplify[D[cMet[x, y, z, t][[2, 3]], t], t >= 1];
c33[x_, y_, z_, 
   t_] = -1/2 Simplify[D[cMet[x, y, z, t][[3, 3]], t], t >= 1];
dens[x_, y_, z_, t_] = \[Rho][r[x, y, z], t];
ord = 2;
met11[x_, y_, z_, t_] := 
  IfThen[x^2 + y^2 < 10^-6, 
   Normal[Series[
      m11[r1, r2, z, t], {r1, 0, ord}, {r2, 0, ord}]] /. {r1 -> x, 
     r2 -> y}, m11[x, y, z, t]];
met12[x_, y_, z_, t_] := 
  IfThen[x^2 + y^2 < 10^-6, 
   Normal[Series[
      m12[r1, r2, z, t], {r1, 0, ord}, {r2, 0, ord}]] /. {r1 -> x, 
     r2 -> y}, m12[x, y, z, t]];
met13[x_, y_, z_, t_] := 
  IfThen[x^2 + y^2 < 10^-6, 
   Normal[Series[
      m13[r1, r2, z, t], {r1, 0, ord}, {r2, 0, ord}]] /. {r1 -> x, 
     r2 -> y}, m13[x, y, z, t]];
met22[x_, y_, z_, t_] := 
  IfThen[x^2 + y^2 < 10^-6, 
   Normal[Series[
      m22[r1, r2, z, t], {r1, 0, ord}, {r2, 0, ord}]] /. {r1 -> x, 
     r2 -> y}, m22[x, y, z, t]];
met23[x_, y_, z_, t_] := 
  IfThen[x^2 + y^2 < 10^-6, 
   Normal[Series[
      m23[r1, r2, z, t], {r1, 0, ord}, {r2, 0, ord}]] /. {r1 -> x, 
     r2 -> y}, m23[x, y, z, t]];
met33[x_, y_, z_, t_] := 
  IfThen[x^2 + y^2 < 10^-6, 
   Normal[Series[
      m33[r1, r2, z, t], {r1, 0, ord}, {r2, 0, ord}]] /. {r1 -> x, 
     r2 -> y}, m33[x, y, z, t]];
cur11[x_, y_, z_, t_] := 
  IfThen[x^2 + y^2 < 10^-6, 
   Normal[Series[
      c11[r1, r2, z, t], {r1, 0, ord}, {r2, 0, ord}]] /. {r1 -> x, 
     r2 -> y}, c11[x, y, z, t]];
cur12[x_, y_, z_, t_] := 
  IfThen[x^2 + y^2 < 10^-6, 
   Normal[Series[
      c12[r1, r2, z, t], {r1, 0, ord}, {r2, 0, ord}]] /. {r1 -> x, 
     r2 -> y}, c12[x, y, z, t]];
cur13[x_, y_, z_, t_] := 
  IfThen[x^2 + y^2 < 10^-6, 
   Normal[Series[
      c13[r1, r2, z, t], {r1, 0, ord}, {r2, 0, ord}]] /. {r1 -> x, 
     r2 -> y}, c13[x, y, z, t]];
cur22[x_, y_, z_, t_] := 
  IfThen[x^2 + y^2 < 10^-6, 
   Normal[Series[
      c22[r1, r2, z, t], {r1, 0, ord}, {r2, 0, ord}]] /. {r1 -> x, 
     r2 -> y}, c22[x, y, z, t]];
cur23[x_, y_, z_, t_] := 
  IfThen[x^2 + y^2 < 10^-6, 
   Normal[Series[
      c23[r1, r2, z, t], {r1, 0, ord}, {r2, 0, ord}]] /. {r1 -> x, 
     r2 -> y}, c23[x, y, z, t]];
cur33[x_, y_, z_, t_] := 
  IfThen[x^2 + y^2 < 10^-6, 
   Normal[Series[
      c33[r1, r2, z, t], {r1, 0, ord}, {r2, 0, ord}]] /. {r1 -> x, 
     r2 -> y}, c33[x, y, z, t]];
dns[x_, y_, z_, t_] := 
  IfThen[x^2 + y^2 < 10^-6, 
   Normal[Series[
      dens[r1, r2, z, t], {r1, 0, ord}, {r2, 0, ord}]] /. {r1 -> x, 
     r2 -> y}, dens[x, y, z, t]];

initialLTBCalc =
{
  Name -> CTThorn <> "_LTB",
  Schedule -> {"IN CCTK_INITIAL after ADMBase_Initial before " <> CTThorn <> "_setCTTrhs"},
  ConditionalOnKeyword -> {"my_initial_data", "LTB"},
  Shorthands -> {detg},
  Equations -> 
  {
    g11        -> met11[x,y,z,1],
    g12        -> met12[x,y,z,1],
    g13        -> met13[x,y,z,1],
    g22        -> met22[x,y,z,1],
    g23        -> met23[x,y,z,1],
    g33        -> met33[x,y,z,1],
    detg       -> Det[{{g11,g12,g13},{g12,g22,g23},{g13,g23,g33}}],

    k11        -> cur11[x,y,z,1],
    k12        -> cur12[x,y,z,1],
    k13        -> cur13[x,y,z,1],
    k22        -> cur22[x,y,z,1],
    k23        -> cur23[x,y,z,1],
    k33        -> cur33[x,y,z,1],

    rho        -> dns[x,y,z,1],
    eps        -> 0,
    prs        -> w rho,
    u[la]      -> 0,
    rhodp      -> dns[x,y,z,1] Sqrt[detg]
  }
};

boundaryLTBCalc =
{
  Name -> CTThorn <> "_LTB_bound",
  Schedule -> {"IN MoL_PostStep"},
  ConditionalOnKeyword -> {"my_initial_data", "LTB"},
  Where -> BoundaryWithGhosts,
  NoSimplify -> True,
  Shorthands -> {detg, det4g, gu[ua,ub], gtu[ua,ub], Gt[ua,lb,lc], f1, f2, u0, h},
  Equations -> 
  {
    g11        -> met11[x,y,z,t],
    g12        -> met12[x,y,z,t],
    g13        -> met13[x,y,z,t],
    g22        -> met22[x,y,z,t],
    g23        -> met23[x,y,z,t],
    g33        -> met33[x,y,z,t],

    k11        -> cur11[x,y,z,t],
    k12        -> cur12[x,y,z,t],
    k13        -> cur13[x,y,z,t],
    k22        -> cur22[x,y,z,t],
    k23        -> cur23[x,y,z,t],
    k33        -> cur33[x,y,z,t],

    detg       -> Det[{{g11,g12,g13},{g12,g22,g23},{g13,g23,g33}}],
    det4g      -> - a detg,
    gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]],
    f1         -> b[ua] u[la] / a^2,
    f2         -> (gu[ua,ub] u[la] u[lb] + 1) / a^2,
    u0         -> (f1 + Sqrt[f1^2 + 4 f2]) / 2,

    h          -> 1 + eps + w,
    W          -> Sqrt[-det4g] u0,
    rho        -> dens[x,y,z,t],
    eps        -> 0,
    prs        -> w rho,
    u[la]      -> 0,
    rhodp      -> rho Sqrt[detg],

    phi        -> detg^(-1/6),
    gt11       -> phi^2 g11,
    gt12       -> phi^2 g12,
    gt13       -> phi^2 g13,
    gt22       -> phi^2 g22,
    gt23       -> phi^2 g23,
    gt33       -> phi^2 g33,

    trK        -> gu[ua,ub] k[la,lb],
    At11       -> phi^2 ( k11 - trK g11 / 3),
    At12       -> phi^2 ( k12 - trK g12 / 3),
    At13       -> phi^2 ( k13 - trK g13 / 3),
    At22       -> phi^2 ( k22 - trK g22 / 3),
    At23       -> phi^2 ( k23 - trK g23 / 3),
    At33       -> phi^2 ( k33 - trK g33 / 3),

    Xt1        -> 0, (* Simplified BCs to speed up compilation -- correct BCs on gamma branch *)
    Xt2        -> 0,
    Xt3        -> 0
  }
};

(******************************************************************************)
(* Convert from primitives *)
(******************************************************************************)

convertFromPrimitivesCalc =
{
  Name -> CTThorn <> "_convertFromPrimitives",
  Schedule -> {"IN CCTK_INITIAL after ADMBase_Initial after " <> CTThorn <> "_FLRW_Pert after " <> CTThorn <> "_Minkowski after " <> CTThorn <> "_MB_bound after " <> CTThorn <> "_LTB_bound after CT_MultiLevel"},
  Shorthands -> {detg, det4g, gu[ua,ub], f1, f2, u0, h},
  Equations -> 
  {
    detg      -> detgExpr,
    det4g     -> - a detg,
    gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]],
    f1        -> b[ua] u[la] / a^2,
    f2        -> (gu[ua,ub] u[la] u[lb] + 1) / a^2,
    u0        -> (f1 + Sqrt[f1^2 + 4 f2]) / 2,
    h         -> 1 + eps + w,

    W  -> Sqrt[-det4g] u0,
    DD -> W rho,
    EE -> W rho eps,
    SS[la] -> W rho h u[la],
    V[ua]  -> gu[ua, ub] u[lb] u0
  }
};


(******************************************************************************)
(* Convert to primitives *)
(******************************************************************************)

usolve = Solve[
       W == Sqrt[-det4g] (a^-2 (b1 u1 + b2 u2 + b3 u3) + 
            Sqrt[a^-4 (b1 u1 + b2 u2 + b3 u3)^2 + 
               4 a^-2 (1 + gu11 u1 u1 + 2 gu12 u1 u2 + 2 gu13 u3 u3 + 
               gu22 u2 u2 + 2 gu23 u2 u3 + gu33 u3 u3)]) / 2 &&
    
       DD == rho W &&
       EE == rho eps W &&
       SS1 == rho h u1 W &&
       SS2 == rho h u2 W &&
       SS3 == rho h u3 W,
     {eps, rho, u1, u2, u3, W}
];

convertToPrimitivesCalc =
{
  Name -> CTThorn <> "_convertToPrimitives",
  Schedule -> {"in MoL_PostStep before SetTmunu after ADMBase_SetADMVars"},
  Where -> Everywhere,
  Shorthands -> {detg, det4g, gu[ua,ub], f1, f2, u0, h},
  Equations -> 
  {
    detg      -> detgExpr,
    det4g     -> - a detg,
    gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]],

    usolve[[1,1]], (* eps *)

    h         -> 1 + eps + w,

    usolve[[1,2]], (* rho *)
    usolve[[1,3]], (* u1 *)
    usolve[[1,4]], (* u2 *)
    usolve[[1,5]], (* u3 *)
    usolve[[1,6]], (* W *)

    prs -> w rho,
    rhodp -> rho Sqrt[detg],

    f1        -> b[ua] u[la] / a^2,
    f2        -> (gu[ua,ub] u[la] u[lb] + 1) / a^2,
    u0        -> (f1 + Sqrt[f1^2 + 4 f2]) / 2,
    V[ua]     -> gu[ua, ub] u[lb] u0
  }
};

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

evolCalc =
{
  Name -> CTThorn <> "_RHS",
  Schedule -> {"IN MoL_CalcRHS"},
  Where -> InteriorNoSync,
  Shorthands -> {dir[ua], detg, det4g, gu[ua,ub], f1, f2, u0, dtW, SS0, h},
  Equations -> 
  {
    dir[ua] -> Sign[b[ua]],

    detg      -> detgExpr,
    det4g     -> -a detg,
    gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]],
    f1        -> b[ua] u[la] / a^2,
    f2        -> (gu[ua,ub] u[la] u[lb] + 1) / a^2,
    u0        -> (f1 + Sqrt[f1^2 + 4 f2]) / 2,
    dtW       -> 0, (*TODO: this has to be constructed from a number of other rhs's.
                            Irrelevant as long as prs=0*) 
    h         -> 1 + eps + prs / rho,
    SS0       -> W rho h u0,

    dot[DD] -> - PD[DD V[ua], la],
    dot[EE] -> - PD[EE V[ua], la] - prs dtW - prs PD[W V[ua], la],
    dot[SS[la]] -> - PD[SS[la] V[ub], lb]  
                   + SS0 PD[-a^2 + g[la,lb] b[ua] b[ub],la] / 2
                   + gu[uc,ue] SS[le] PD[g[lb,lc] b[ub],la]
                   + gu[ub,ud] gu[uc,ue] SS[ld] SS[le] PD[g[lb,lc],la] / (2 SS0)
                   - Sqrt[-det4g] PD[prs, la]
  }
};

RHSStaticBoundaryCalc =
{
  Name -> CTThorn <> "_RHSStaticBoundary",
  Schedule -> {"IN MoL_CalcRHS"},
  ConditionalOnKeyword -> {"my_rhs_boundary_condition", "static"},
  Where -> Boundary,
  Equations -> 
  {
    dot[DD]       -> 0,
    dot[EE]       -> 0,
    dot[SS[la]]   -> 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 -> CTThorn <> "_InitRHS",
  Schedule -> {"AT analysis BEFORE " <> CTThorn <> "_RHS"},
  Where -> Everywhere,
  Equations -> 
  {
    dot[DD]       -> 0,
    dot[EE]       -> 0,
    dot[SS[la]]   -> 0
  }
};

(******************************************************************************)
(* Populate the energy-momentum tensor *)
(******************************************************************************)

addToTmunuCalc =
{
  Name -> CTThorn <> "_addToTmunu",
  Schedule -> {"IN AddToTmunu"},
  ConditionalOnKeyword -> {"coupling", "yes"},
  Where -> Everywhere,
  Shorthands -> {detg, det4g, gu[ua,ub], f1, f2, u0u, u0l, bsq, rhocrit, rhof},
  Equations -> 
  {
    detg      -> detgExpr,
    det4g     -> - a detg,
    gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]],
    f1        -> b[ua] u[la] / a^2,
    f2        -> (gu[ua,ub] u[la] u[lb] + 1) / a^2,
    u0u       -> (f1 + Sqrt[f1^2 + 4 f2]) / 2,
    bsq       -> g[la,lb] b[ua] b[ub],
    u0l       -> (-a^2 + bsq) u0u + b[ua] u[la],

    T00       -> (rho + prs) u0l^2 + (prs - Lambda/(8 Pi)) (- a^2 + bsq),
    T0[la]    -> (rho + prs) u0l u[la] + (prs - Lambda/(8 Pi)) g[la,lb] b[ub],
    T[la,lb]  -> (rho + prs) u[la] u[lb] + (prs - Lambda/(8 Pi)) g[la,lb]
  }
};

setCTTrhs =
{
  Name -> CTThorn <> "_setCTTrhs",
  Schedule -> {"AT INITIAL before CT_MultiLevel"},
  Where -> Interior,
  Shorthands -> {detg, det4g, gu[ua,ub], f1, f2, u0u, u0l, bsq, t00, t0[la], t[la,lb]},
  Equations -> 
  {
    detg      -> detgExpr,
    det4g     -> - a detg,
    gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]],
    f1        -> b[ua] u[la] / a^2,
    f2        -> (gu[ua,ub] u[la] u[lb] + 1) / a^2,
    u0u       -> (f1 + Sqrt[f1^2 + 4 f2]) / 2,
    bsq       -> g[la,lb] b[ua] b[ub],
    u0l       -> (-a^2 + bsq) u0u + b[ua] u[la],

    t00       -> (rho + prs) u0l^2 + (prs - Lambda/(8 Pi)) (- a^2 + bsq),
    t0[la]    -> (rho + prs) u0l u[la] + (prs - Lambda/(8 Pi)) g[la,lb] b[ub],
    t[la,lb]  -> (rho + prs) u[la] u[lb] + (prs - Lambda/(8 Pi)) g[la,lb],

    cttS      -> 2 pi (t00 - 2 b[ua] t0[la] + b[ua] b[ub] t[la,lb]) / a^2,
    cttV[lb]  -> 8 pi (t0[lb] - b[ua] t[la,lb]) / a
  }
};

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

boundaryCalc =
{
  Name -> CTThorn <> "_boundary",
  Schedule -> {"IN MoL_PostStep"},
  ConditionalOnKeyword -> {"my_boundary_condition", "Minkowski"},
  Where -> BoundaryWithGhosts,
  Equations -> 
  {
    DD         -> 0,
    EE         -> 0,
    SS[la]     -> 0
  }
};


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

inheritedImplementations =
  {"ADMBase", "TmunuBase", "ML_BSSN"};

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


extendedKeywordParameters =
{
  {
    Name -> "ADMBase::initial_data",
    AllowedValues -> {CTThorn}
  },
  {
    Name -> "ADMBase::initial_lapse",
    AllowedValues -> {CTThorn}
  },
  {
    Name -> "ADMBase::initial_shift",
    AllowedValues -> {CTThorn} 
  }
};

keywordParameters =
{
  {
    Name -> "my_initial_data",
    AllowedValues -> {"Minkowski", "pFLRW", "SMB", "LTB"},
    Default -> "Minkowski"
  },
  {
    Name -> "my_rhs_boundary_condition",
    Visibility -> "restricted",
    AllowedValues -> {"none", "static", "radiative"},
    Default -> "none"
  },
  {
    Name -> "my_boundary_condition",
    AllowedValues -> {"none", "Minkowski"},
    Default -> "none"
  },
  {
    Name -> "coupling",
    AllowedValues -> {"yes", "no"},
    Default -> "no"
  }
};

intParameters =
{
  {
    Name -> nmodes,
    Description -> "Number of Fourier modes in each direction", 
    Default -> 0
  },
  {
    Name -> fdOrder,
    Default -> derivOrder,
    AllowedValues -> {2,4,6,8}
  }
};

realParameters =
{
  {
    Name -> w,
    Description -> "Ratio of pressure to density",
    AllowedValues -> {{Value -> "0:0", Description -> "Only dust is supported at this time"}},
    Default -> 0
  },
  {
    Name -> Lambda,
    Description -> "Cosmological constant",
    Default -> 0
  },
  {
    Name -> "pflrw_t0",
    Description -> "Initial time in RW models",
    Default -> 1
  },
  {
    Name -> "pflrw_a0",
    Description -> "Initial value of the scale factor in RW models",
    Default -> 1
  },
  {
    Name -> "pflrw_H0",
    Description -> "Initial value of the Hubble rate in RW models",
    Default -> 1
  },
  {
    Name -> "pflrw_omegaM",
    Description -> "Initial value of the matter density parameter",
    Default -> 1(*,
    AllowedValues -> {1}*)
  },
  {
    Name -> "pflrw_kx[20]",
    Description -> "Wave numbers of the x-modes",
    Default -> 0
  },
  {
    Name -> "pflrw_ky[20]",
    Description -> "Wave numbers of the y-modes",
    Default -> 0
  },
  {
    Name -> "pflrw_kz[20]",
    Description -> "Wave numbers of the z-modes",
    Default -> 0
  },
  {
    Name -> "pflrw_ax[20]",
    Description -> "Amplitudes of the x-modes",
    Default -> 0
  },
  {
    Name -> "pflrw_ay[20]",
    Description -> "Amplitudes of the y-modes",
    Default -> 0
  },
  {
    Name -> "pflrw_az[20]",
    Description -> "Amplitudes of the z-modes",
    Default -> 0
  },
  {
    Name -> "ltb_rScal",
    Description -> "Rescaling of radial coordinate in spherically-symmetric models",
    Default -> 1
  }
};

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

calculations =
{
  initialMinkowskiCalc,
  initialFLRWPertCalc,
  initialMBCalc,
  boundaryMBCalc,
  initialLTBCalc,
  boundaryLTBCalc,
  convertFromPrimitivesCalc,
  evolCalc,
  initRHSCalc,
  RHSStaticBoundaryCalc,
  boundaryCalc,
  convertToPrimitivesCalc,
  addToTmunuCalc,
  setCTTrhs
};

CreateKrancThornTT [groups, ".", CTThorn,
  Calculations -> calculations,
  DeclaredGroups -> declaredGroupNames,
  PartialDerivatives -> derivatives,
  EvolutionTimelevels -> evolutionTimelevels,
  DefaultEvolutionTimelevels -> 3,
  UseJacobian -> useJacobian,
  UseLoopControl -> True,
  UseVectors -> False,
  InheritedImplementations -> inheritedImplementations,
  ExtendedKeywordParameters -> extendedKeywordParameters,
  KeywordParameters -> keywordParameters,
  IntParameters -> intParameters,
  RealParameters -> realParameters
];

];

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

(* derivative order: 2, 4, 6, 8, ... *)
(* useJacobian: False or True *)
(* timelevels: 2 or 3
   (keep this at 3; this is better chosen with a run-time parameter) *)

createCode[4, True, 3];
