SetEnhancedTimes[False];
SetSourceLanguage["C"];

prefix = "ML_";
thorn  = prefix <> "Kretschmann";

Map [DefineTensor, {
     g, gu, G, R, Km, trK, phi, gt, At, e4phi, Dphi, D2phi,
     Gl, Riem, Riemu, K, Ku, CDK, CDKu, Kretsch, Dg, D2g
     }
];

declaredGroups = {
    SetGroupName[CreateGroupFromTensor[Kretsch], prefix <> "Kretschmann"]
};

declaredGroupNames = Map[First, declaredGroups];

extraGroups = {
    { "ML_BSSN::ML_metric",     { gt11, gt12, gt13, gt22, gt23, gt33 } },
    { "ML_BSSN::ML_curv",       { At11, At12, At13, At22, At23, At33 } },
    { "ML_BSSN::ML_trace_curv", { trK } },
    { "ML_BSSN::ML_log_confac", { phi } }
};

groups = Join[declaredGroups, extraGroups];

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

PD = PDstandardNth;

Map[AssertSymmetricIncreasing,
     { g[la, lb], K[la, lb], R[la, lb],
       gt[la, lb], At[la, lb], D2phi[la, lb]}];
Map [AssertSymmetricIncreasing,
     { gu[ua, ub], gtu[ua, ub], Ku[ua, ub]}];
AssertSymmetricIncreasing [G[ua,lb,lc], lb, lc];
AssertSymmetricIncreasing [Gl[la,lb,lc], lb, lc];
AssertSymmetricIncreasing [Dg[la, lb,lc], la, lb];
AssertSymmetricIncreasing [CDK[la, lb, lc], la, lb];
AssertSymmetricIncreasing [CDKu[ua, ub, uc], ua, ub];
AssertSymmetricIncreasing [D2g[la, lb, lc, ld], la, lb];
AssertSymmetricIncreasing [D2g[la, lb, lc, ld], lc, ld];

kretschmannCalc =
{
    Name -> thorn <> "_kretschmann",
    Schedule -> Automatic,
    Where -> Interior,
    After -> "MoL_PostStep",
    Shorthands -> {term1, term2, term3,
                   e4phi, Dphi[la], D2phi[la, lb],
                   Dg[la, lb, lc], D2g[la, lb, lc, ld],
                   g[la, lb], gu[ua, ub],
                   Gl[la, lb, lc], G[ua, lb, lc],
                   Riem[la, lb, lc, ld], Riemu[ua, ub, uc, ud], R[la, lb],
                   K[la, lb], Km[ua, lb], Ku[ua, ub],
                   CDK[la, lb, lc], CDKu[ua, ub, uc]},
    Equations -> {
        e4phi          -> IfThen[conformalMethod == 1, 1/phi^2, Exp[4 phi]],
        Dphi[la]       -> IfThen[conformalMethod == 1,
                                 - PD[phi, la] / (2 phi ),
                                 PD[phi, la]],
        D2phi[la, lb]  -> IfThen[conformalMethod == 1,
                                 - 1/ (2 phi) (PD[phi, la, lb] -
                                               PD[phi, la] PD[phi, lb] / phi),
                                 PD[phi, la, lb]],
        Dg[la, lb, lc] -> e4phi * (4 Dphi[lc] gt[la, lb] + PD[gt[la, lb], lc]),
        D2g[la, lb, lc, ld] -> e4phi * (4 gt[la, lb] (4 Dphi[lc] Dphi[ld] +
                                                      D2phi[lc, ld]) +
                                        4 Dphi[lc] PD[gt[la, lb], ld] +
                                        4 Dphi[ld] PD[gt[la, lb], lc] +
                                        PD[gt[la, lb], lc, ld]),

        g[la, lb]     -> e4phi gt[la, lb],
        gu[ua, ub]    -> MatrixInverse[g[ua, ub]],

        Gl[la, lb, lc]  -> 1/2 (Dg[la, lb, lc] + Dg[la, lc, lb] - Dg[lb, lc, la]),
        G[ua, lb, lc]   -> gu[ua, ud] Gl[ld, lb, lc],

        Riem[la, lb, lc, ld]  -> 1/2 (D2g[la, ld, lb, lc] + D2g[lb, lc, la, ld] -
                                      D2g[lb, ld, la, lc] - D2g[la, lc, lb, ld]) +
                                 Gl[le, la, ld] G[ue, lb, lc] -
                                 Gl[le, la, lc] G[ue, lb, ld],
        (*
         * kranc properly zeroes the zero-from-antisymmetry terms from Riemann above,
         * but fails to do so when raising indices;
         * therefore indicate the (anti)symmetries explicitly.
         * results in much faster code
         *)
        Riemu[ua, ub, uc, ud] -> gu[ua, ue] gu[ub, uf] gu[uc, ug] gu[ud, uh]
                                 1/4 (Riem[le, lf, lg, lh] - Riem[lf, le, lg, lh] -
                                      Riem[le, lf, lh, lg] + Riem[lf, le, lh, lg]),
        R[la, lb]             -> gu[uc, ud] Riem[lc, la, ld, lb],

        K[la, lb]  -> e4phi At[la, lb] + 1/3 trK g[la, lb],
        Km[ua, lb] -> gu[ua, uc] K[lc, lb],
        Ku[ua, ub] -> gu[ub, uc] Km[ua, lc],

        CDK[la, lb, lc]  -> e4phi * (4 Dphi[lc] At[la, lb] + PD[At[la, lb], lc]) +
                            1/3 (PD[trK, lc] g[la, lb] + trK Dg[la, lb, lc]) -
                            G[ud, la, lc] K[ld, lb] - G[ud, lb, lc] K[ld, la],
        CDKu[ua, ub, uc] -> gu[ua, ud] gu[ub, ue] gu[uc, uf] CDK[ld, le, lf],

        term1  -> (Riem[la, lb, lc, ld] + K[la, lc] K[lb, ld] - K[la, ld] K[lb, lc])
                  (Riemu[ua, ub, uc, ud] + Ku[ua, uc] Ku[ub, ud] - Ku[ua, ud] Ku[ub, uc]),
        term2  -> 8 (CDK[la, lb, lc] CDKu[ua, uc, ub] - CDK[la, lb, lc] CDKu[ua, ub, uc]),
        term3  -> 4 gu[ua, uc] gu[ub, ud]
                  (R[la, lb] + trK K[la, lb] - K[la, le] Km[ue, lb])
                  (R[lc, ld] + trK K[lc, ld] - K[lc, le] Km[ue, ld]),

        Kretsch -> term1 + term2 + term3
    }
};

intParameters = {
    (* FIXME: pull this parameter value from ML_BSSN ? *)
    {
      Name -> conformalMethod,
      Description -> "Treatment of conformal factor",
      AllowedValues -> {{Value -> "0", Description -> "phi method"},
                        {Value -> "1", Description -> "W method"}},
      Default -> 0
    },
    {
      Name          -> fdOrder,
      Default       -> 4,
      AllowedValues -> {2, 4, 6, 8}
    }
};

CreateKrancThornTT[ groups, ".", thorn,
    Calculations               -> { kretschmannCalc },
    PartialDerivatives         -> derivatives,
    EvolutionTimelevels        -> 3,
    DeclaredGroups -> declaredGroupNames,
    DefaultEvolutionTimelevels -> 3,
    UseJacobian                -> True,
    UseLoopControl             -> True,
    UseVectors                 -> True,
    IntParameters              -> intParameters,
    InheritedImplementations   -> { "ML_BSSN" }
];