ML_Kretschmann: add ζ calculation

5 files changed, 42 insertions, 0 deletions

diff --git a/ML_Kretschmann/interface.ccl b/ML_Kretschmann/interface.ccl
index 16d2721..c6e4182 100644
--- a/ML_Kretschmann/interface.ccl
+++ b/ML_Kretschmann/interface.ccl
@@ -32,3 +32,8 @@ CCTK_REAL ML_Kretschmann type=GF timelevels=1 tags='tensortypealias="Scalar" ten
 {
   Kretsch
 } "ML_Kretschmann"
+
+CCTK_REAL ML_zeta type=GF timelevels=1 tags='tensortypealias="Scalar" tensorweight=0'
+{
+  zeta
+} "ML_zeta"
diff --git a/ML_Kretschmann/schedule.ccl b/ML_Kretschmann/schedule.ccl
index 6d8259a..09788ab 100644
--- a/ML_Kretschmann/schedule.ccl
+++ b/ML_Kretschmann/schedule.ccl
@@ -2,6 +2,7 @@
 
 
 STORAGE: ML_Kretschmann[other_timelevels]
+STORAGE: ML_zeta[other_timelevels]
 
 schedule ML_Kretschmann_Startup at STARTUP
 {
@@ -37,6 +38,7 @@ schedule ML_Kretschmann_kretschmann in ML_Kretschmann_kretschmann_group
   READS: ML_BSSN::gt33(Everywhere)
   READS: ML_BSSN::trK(Everywhere)
   WRITES: ML_Kretschmann::Kretsch(Interior)
+  WRITES: ML_zeta::zeta(Interior)
 } "ML_Kretschmann_kretschmann"
 
 schedule ML_Kretschmann_kretschmann_SelectBCs in ML_Kretschmann_kretschmann_bc_group
@@ -44,6 +46,7 @@ schedule ML_Kretschmann_kretschmann_SelectBCs in ML_Kretschmann_kretschmann_bc_g
   LANG: C
   OPTIONS: level
   SYNC: ML_Kretschmann
+  SYNC: ML_zeta
 } "ML_Kretschmann_kretschmann_SelectBCs"
 
 schedule group ApplyBCs as ML_Kretschmann_kretschmann_ApplyBCs in ML_Kretschmann_kretschmann_bc_group after ML_Kretschmann_kretschmann_SelectBCs
diff --git a/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc b/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc
index 83851f4..a43da4c 100644
--- a/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc
+++ b/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc
@@ -38,6 +38,9 @@ extern "C" void ML_Kretschmann_kretschmann_SelectBCs(CCTK_ARGUMENTS)
   ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_Kretschmann::ML_Kretschmann","flat");
   if (ierr < 0)
     CCTK_WARN(1, "Failed to register flat BC for ML_Kretschmann::ML_Kretschmann.");
+  ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, GenericFD_GetBoundaryWidth(cctkGH), -1 /* no table */, "ML_Kretschmann::ML_zeta","flat");
+  if (ierr < 0)
+    CCTK_WARN(1, "Failed to register flat BC for ML_Kretschmann::ML_zeta.");
   return;
 }
 
@@ -2066,6 +2069,31 @@ static void ML_Kretschmann_kretschmann_Body(const cGH* restrict const cctkGH, co
     
     CCTK_REAL_VEC KretschL CCTK_ATTRIBUTE_UNUSED = 
       kadd(term1,kadd(term2,term3));
+
+    {
+        double xx = x[index];
+        double x2 = SQR(xx);
+        double x4 = SQR(x2);
+        double reg_part = SQR(K22) / SQR(g22) - 0.25 * (SQR(Dg221) * gu11 + gu33 * SQR(Dg223) + 2.0 * Dg221 * Dg223 * gu13) / SQR(g22);
+        double lambda = (x2 * g22 * (1.0 - gu11 * g22) - xx * x2 * g22 * (Dg221 * gu11 + gu13 * Dg223)) / (x4 * SQR(g22));
+
+        if (fabs(xx) < 1e-14) {
+            double dphi_x = PDstandardNth1phi;
+            double d2phi_xx = PDstandardNth11phi;
+            double dgt11_x  = PDstandardNth1gt11;
+            double d2gt11_xx = PDstandardNth11gt11;
+            double Dgu131 = -Dg131 * (gu11 * gu33 + SQR(gu13));
+            double Dgu111 = 0.0;
+            double D2g1111 = - 6.0 * SQR(dphi_x) * gt11L / SQR(SQR(phiL))
+                             - 2.0 * d2phi_xx * gt11L / (phiL * SQR(phiL))
+                             - 4.0 * dphi_x * dgt11_x / (phiL * SQR(phiL))
+                             + d2gt11_xx / SQR(phiL);
+            double D2gu1111 = -(SQR(gu11) * D2g1111 + 2.0 * Dg131 * (Dgu111 * gu13 + gu11 * Dgu131) + SQR(gu13) * D2g3311);
+            double d2x_guxx_gyy = D2gu1111 * g22 + gu11 * D2g2211;
+            lambda = -(gu11 * D2g2211 + Dgu131 * Dg223) / g22 - 0.5 * d2x_guxx_gyy / g22;
+        }
+        zeta[index] = reg_part + lambda;
+    }
     
     /* Copy local copies back to grid functions */
     vec_store_partial_prepare(i,vecimin,vecimax);
diff --git a/ML_Kretschmann/src/RegisterSymmetries.cc b/ML_Kretschmann/src/RegisterSymmetries.cc
index d52728a..9e3ce3f 100644
--- a/ML_Kretschmann/src/RegisterSymmetries.cc
+++ b/ML_Kretschmann/src/RegisterSymmetries.cc
@@ -21,4 +21,9 @@ extern "C" void ML_Kretschmann_RegisterSymmetries(CCTK_ARGUMENTS)
   sym[2] = 1;
   SetCartSymVN(cctkGH, sym, "ML_Kretschmann::Kretsch");
   
+  /* Register symmetries of grid functions */
+  sym[0] = 1;
+  sym[1] = 1;
+  sym[2] = 1;
+  SetCartSymVN(cctkGH, sym, "ML_Kretschmann::zeta");
 }
diff --git a/ML_Kretschmann_Helper/src/SetGroupTags.c b/ML_Kretschmann_Helper/src/SetGroupTags.c
index 5bf2bf4..17b939a 100644
--- a/ML_Kretschmann_Helper/src/SetGroupTags.c
+++ b/ML_Kretschmann_Helper/src/SetGroupTags.c
@@ -40,6 +40,7 @@ int ML_Kretschmann_SetGroupTags (void)
   int const checkpoint = 0;
 
   set_group_tags (checkpoint, checkpoint, 0, "ML_Kretschmann::ML_Kretschmann");
+  set_group_tags (checkpoint, checkpoint, 0, "ML_Kretschmann::ML_zeta");
 
   return 0;
 }