Add a thorn for computing the Kretschmann scalar.

7 files changed, 205 insertions, 0 deletions

diff --git a/ML_Kretschmann_Helper/configuration.ccl b/ML_Kretschmann_Helper/configuration.ccl
new file mode 100644
index 0000000..edaa5b7
--- /dev/null
+++ b/ML_Kretschmann_Helper/configuration.ccl
@@ -0,0 +1 @@
+REQUIRES THORNS: CartGrid3D CoordGauge
diff --git a/ML_Kretschmann_Helper/interface.ccl b/ML_Kretschmann_Helper/interface.ccl
new file mode 100644
index 0000000..ded78e0
--- /dev/null
+++ b/ML_Kretschmann_Helper/interface.ccl
@@ -0,0 +1,3 @@
+IMPLEMENTS: ML_Kretschmann_Helper
+
+INHERITS: ADMBase CoordGauge ML_Kretschmann
diff --git a/ML_Kretschmann_Helper/param.ccl b/ML_Kretschmann_Helper/param.ccl
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/ML_Kretschmann_Helper/param.ccl
diff --git a/ML_Kretschmann_Helper/schedule.ccl b/ML_Kretschmann_Helper/schedule.ccl
new file mode 100644
index 0000000..09e44d3
--- /dev/null
+++ b/ML_Kretschmann_Helper/schedule.ccl
@@ -0,0 +1,5 @@
+SCHEDULE ML_Kretschmann_SetGroupTags AT startup BEFORE Driver_Startup
+{
+    LANG: C
+    OPTIONS: meta
+} "Set checkpointing and prolongation group tags"
diff --git a/ML_Kretschmann_Helper/src/SetGroupTags.c b/ML_Kretschmann_Helper/src/SetGroupTags.c
new file mode 100644
index 0000000..5bf2bf4
--- /dev/null
+++ b/ML_Kretschmann_Helper/src/SetGroupTags.c
@@ -0,0 +1,46 @@
+#include <cctk.h>
+#include <cctk_Parameters.h>
+#include <util_Table.h>
+
+#include <assert.h>
+
+static void set_group_tags (int const checkpoint, int const persistent,
+                            int const prolongate, char const * restrict const gn)
+{
+  assert (gn);
+
+  int const gi = CCTK_GroupIndex(gn);
+  assert (gi >= 0);
+
+  int const table = CCTK_GroupTagsTableI(gi);
+  assert(table >= 0);
+
+  if (!checkpoint) {
+    int const ierr = Util_TableSetString (table, "no", "Checkpoint");
+    assert (!ierr);
+  }
+
+  if (!persistent) {
+    int const ierr = Util_TableSetString (table, "no", "Persistent");
+    assert (!ierr);
+  }
+
+  if (!prolongate) {
+    int const iret = Util_TableDeleteKey (table, "ProlongationParameter");
+    assert (iret == 0 || iret == UTIL_ERROR_TABLE_NO_SUCH_KEY);
+    int const ierr = Util_TableSetString (table, "none", "Prolongation");
+    assert (!ierr);
+  }
+}
+
+int ML_Kretschmann_SetGroupTags (void)
+{
+  DECLARE_CCTK_PARAMETERS;
+
+  int const checkpoint = 0;
+
+  set_group_tags (checkpoint, checkpoint, 0, "ML_Kretschmann::ML_Kretschmann");
+
+  return 0;
+}
+
diff --git a/ML_Kretschmann_Helper/src/make.code.defn b/ML_Kretschmann_Helper/src/make.code.defn
new file mode 100644
index 0000000..df06acc
--- /dev/null
+++ b/ML_Kretschmann_Helper/src/make.code.defn
@@ -0,0 +1,2 @@
+# -*-Makefile-*-
+SRCS = SetGroupTags.c
diff --git a/m/McLachlan_Kretschmann.m b/m/McLachlan_Kretschmann.m
new file mode 100644
index 0000000..a8ab872
--- /dev/null
+++ b/m/McLachlan_Kretschmann.m
@@ -0,0 +1,148 @@
+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" }
+];