Compute a function Kzeta = (K - 12 zeta)/rho^2cartoon

4 files changed, 120 insertions, 1 deletions

diff --git a/ML_Kretschmann/interface.ccl b/ML_Kretschmann/interface.ccl
index c6e4182..eed9f3c 100644
--- a/ML_Kretschmann/interface.ccl
+++ b/ML_Kretschmann/interface.ccl
@@ -35,5 +35,6 @@ CCTK_REAL ML_Kretschmann type=GF timelevels=1 tags='tensortypealias="Scalar" ten
 
 CCTK_REAL ML_zeta type=GF timelevels=1 tags='tensortypealias="Scalar" tensorweight=0'
 {
-  zeta
+  zeta,
+  Kzeta
 } "ML_zeta"
diff --git a/ML_Kretschmann/schedule.ccl b/ML_Kretschmann/schedule.ccl
index 09788ab..416b188 100644
--- a/ML_Kretschmann/schedule.ccl
+++ b/ML_Kretschmann/schedule.ccl
@@ -39,6 +39,7 @@ schedule ML_Kretschmann_kretschmann in ML_Kretschmann_kretschmann_group
   READS: ML_BSSN::trK(Everywhere)
   WRITES: ML_Kretschmann::Kretsch(Interior)
   WRITES: ML_zeta::zeta(Interior)
+  WRITES: ML_zeta::Kzeta(Interior)
 } "ML_Kretschmann_kretschmann"
 
 schedule ML_Kretschmann_kretschmann_SelectBCs in ML_Kretschmann_kretschmann_bc_group
diff --git a/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc b/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc
index 89461ea..deb01dd 100644
--- a/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc
+++ b/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc
@@ -2122,6 +2122,116 @@ static void ML_Kretschmann_kretschmann_Body(const cGH* restrict const cctkGH, co
   CCTK_ENDLOOP3STR(ML_Kretschmann_kretschmann);
 }
 
+static void ML_Kretschmann_kzeta(const cGH* restrict const cctkGH, const int dir, const int face, const CCTK_REAL normal[3], const CCTK_REAL tangentA[3], const CCTK_REAL tangentB[3], const int imin[3], const int imax[3], const int n_subblock_gfs, CCTK_REAL* restrict const subblock_gfs[])
+{
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+  
+  
+  /* Include user-supplied include files */
+  
+  /* Initialise finite differencing variables */
+  const ptrdiff_t di CCTK_ATTRIBUTE_UNUSED = 1;
+  const ptrdiff_t dj CCTK_ATTRIBUTE_UNUSED = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  const ptrdiff_t dk CCTK_ATTRIBUTE_UNUSED = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0);
+  
+
+  /* Loop over the grid points */
+  const int imin0=imin[0];
+  const int imin1=imin[1];
+  const int imin2=imin[2];
+  const int imax0=imax[0];
+  const int imax1=imax[1];
+  const int imax2=imax[2];
+  #pragma omp parallel // reduction(+: vec_iter_counter, vec_op_counter, vec_mem_counter)
+  CCTK_LOOP3STR(ML_Kretschmann_kzeta,
+    i,j,k, imin0,imin1,imin2, imax0,imax1,imax2,
+    cctk_ash[0],cctk_ash[1],cctk_ash[2],
+    vecimin,vecimax, CCTK_REAL_VEC_SIZE)
+  {
+    const ptrdiff_t index CCTK_ATTRIBUTE_UNUSED = di*i + dj*j + dk*k;
+    // vec_iter_counter+=CCTK_REAL_VEC_SIZE;
+    
+    /* Assign local copies of grid functions */
+    
+    CCTK_REAL_VEC gt22L CCTK_ATTRIBUTE_UNUSED = vec_load(gt22[index]);
+    CCTK_REAL_VEC gt22p1 CCTK_ATTRIBUTE_UNUSED = vec_load(gt22[index + 1]);
+    CCTK_REAL_VEC gt22p2 CCTK_ATTRIBUTE_UNUSED = vec_load(gt22[index + 2]);
+    CCTK_REAL_VEC phiL CCTK_ATTRIBUTE_UNUSED = vec_load(phi[index]);
+    CCTK_REAL_VEC phip1 CCTK_ATTRIBUTE_UNUSED = vec_load(phi[index + 1]);
+    CCTK_REAL_VEC phip2 CCTK_ATTRIBUTE_UNUSED = vec_load(phi[index + 2]);
+
+    CCTK_REAL_VEC KretschL CCTK_ATTRIBUTE_UNUSED = vec_load(Kretsch[index]);
+    CCTK_REAL_VEC zetaL CCTK_ATTRIBUTE_UNUSED = vec_load(zeta[index]);
+
+    CCTK_REAL_VEC Kretschp1 CCTK_ATTRIBUTE_UNUSED = vec_load(Kretsch[index + 1]);
+    CCTK_REAL_VEC Kretschp2 CCTK_ATTRIBUTE_UNUSED = vec_load(Kretsch[index + 2]);
+    CCTK_REAL_VEC zetap1 CCTK_ATTRIBUTE_UNUSED = vec_load(zeta[index + 1]);
+    CCTK_REAL_VEC zetap2 CCTK_ATTRIBUTE_UNUSED = vec_load(zeta[index + 2]);
+    
+    CCTK_REAL_VEC ktwo = ToReal(2.0);
+
+    CCTK_REAL_VEC xx  = vec_load(x[index]);
+    CCTK_REAL_VEC xp1 = vec_load(x[index + 1]);
+    CCTK_REAL_VEC xp2 = vec_load(x[index + 2]);
+
+    CCTK_REAL_VEC absx = kfabs(xx);
+    CCTK_REAL_VEC eps  = ToReal(1e-8);
+    CCTK_BOOLEAN_VEC origin = kcmplt(absx, eps);
+
+    
+    CCTK_REAL_VEC e4phi CCTK_ATTRIBUTE_UNUSED;
+    CCTK_REAL_VEC e4phip1 CCTK_ATTRIBUTE_UNUSED;
+    CCTK_REAL_VEC e4phip2 CCTK_ATTRIBUTE_UNUSED;
+    
+    if (conformalMethod == 1)
+    {
+      e4phi   = kdiv(ToReal(1),SQR(phiL));
+      e4phip1 = kdiv(ToReal(1),SQR(phip1));
+      e4phip2 = kdiv(ToReal(1),SQR(phip2));
+    }
+    else
+    {
+      e4phi   = kexp(kmul(phiL,ToReal(4)));
+      e4phip1 = kexp(kmul(phip1,ToReal(4)));
+      e4phip2 = kexp(kmul(phip2,ToReal(4)));
+    }
+
+    CCTK_REAL_VEC g22   = kmul(gt22L, e4phi);
+    CCTK_REAL_VEC g22p1 = kmul(gt22p1,e4phip1);
+    CCTK_REAL_VEC g22p2 = kmul(gt22p2,e4phip2);
+
+    CCTK_REAL_VEC KzetaL;
+    {
+        CCTK_REAL_VEC fact_interp_p1 = ToReal(0.5625);
+        CCTK_REAL_VEC fact_interp_p2 = ToReal(-0.0625);
+
+        CCTK_REAL_VEC ksix   = ToReal(6.0);
+        CCTK_REAL_VEC ktwelve = ToReal(12.0);
+
+        CCTK_REAL_VEC x2   = SQR(xx);
+        CCTK_REAL_VEC x2p1 = SQR(xp1);
+        CCTK_REAL_VEC x2p2 = SQR(xp2);
+
+        CCTK_REAL_VEC rho2   = g22   * x2;
+        CCTK_REAL_VEC rho2p1 = g22p1 * x2p1;
+        CCTK_REAL_VEC rho2p2 = g22p2 * x2p2;
+
+        CCTK_REAL_VEC kzeta_reg  = (KretschL - ktwelve * SQR(zetaL)) / rho2;
+        CCTK_REAL_VEC kzetap1    = (Kretschp1 - ktwelve * SQR(zetap1)) / rho2p1;
+        CCTK_REAL_VEC kzetap2    = (Kretschp2 - ktwelve * SQR(zetap2)) / rho2p2;
+        CCTK_REAL_VEC kzeta_sing = ktwo * (fact_interp_p1 * kzetap1 + fact_interp_p2 * kzetap2);
+
+        KzetaL = kifthen(origin, kzeta_sing, kzeta_reg);
+    }
+    
+    /* Copy local copies back to grid functions */
+    vec_store_partial_prepare(i,vecimin,vecimax);
+    vec_store_nta_partial(Kzeta[index], KzetaL);
+  }
+  CCTK_ENDLOOP3STR(ML_Kretschmann_kzeta);
+}
+
 extern "C" void ML_Kretschmann_kretschmann(CCTK_ARGUMENTS)
 {
   DECLARE_CCTK_ARGUMENTS;
@@ -2139,6 +2249,7 @@ extern "C" void ML_Kretschmann_kretschmann(CCTK_ARGUMENTS)
   }
   
   GenericFD_LoopOverInterior(cctkGH, ML_Kretschmann_kretschmann_Body);
+  GenericFD_LoopOverInterior(cctkGH, ML_Kretschmann_kzeta);
   
   if (verbose > 1)
   {
diff --git a/ML_Kretschmann/src/RegisterSymmetries.cc b/ML_Kretschmann/src/RegisterSymmetries.cc
index 9e3ce3f..ad4d1ce 100644
--- a/ML_Kretschmann/src/RegisterSymmetries.cc
+++ b/ML_Kretschmann/src/RegisterSymmetries.cc
@@ -26,4 +26,10 @@ extern "C" void ML_Kretschmann_RegisterSymmetries(CCTK_ARGUMENTS)
   sym[1] = 1;
   sym[2] = 1;
   SetCartSymVN(cctkGH, sym, "ML_Kretschmann::zeta");
+
+  /* Register symmetries of grid functions */
+  sym[0] = 1;
+  sym[1] = 1;
+  sym[2] = 1;
+  SetCartSymVN(cctkGH, sym, "ML_Kretschmann::Kzeta");
 }