ML_Kretschmann: update for cartoon

1 files changed, 45 insertions, 37 deletions

diff --git a/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc b/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc
index 94b138c..83851f4 100644
--- a/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc
+++ b/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc
@@ -686,6 +686,51 @@ static void ML_Kretschmann_kretschmann_Body(const cGH* restrict const cctkGH, co
       default:
         CCTK_BUILTIN_UNREACHABLE();
     }
+
+    double xx = x[index];
+    int origin = fabs(xx) < 1e-8;
+
+    PDstandardNth2trK    = 0;
+
+    PDstandardNth2gt11  = 0;
+    PDstandardNth22gt11 = origin ? PDstandardNth11gt22 : PDstandardNth1gt11 / xx - 2 * (gt11[index] - gt22[index]) / (xx * xx);
+    PDstandardNth12gt11 = 0;
+    PDstandardNth23gt11 = 0;
+    PDstandardNth2gt12  = origin ? PDstandardNth1gt11 - PDstandardNth1gt22 : (gt11[index] - gt22[index]) / xx;
+    PDstandardNth22gt12 = 0;
+    PDstandardNth12gt12 = origin ? (PDstandardNth11gt11 - PDstandardNth11gt22) / 2 : (PDstandardNth1gt11 - PDstandardNth1gt22) / xx - (gt11[index] - gt22[index]) / (xx * xx);
+    PDstandardNth23gt12 = origin ? PDstandardNth13gt11 - PDstandardNth13gt22 : (PDstandardNth3gt11  - PDstandardNth3gt22) / xx;
+    PDstandardNth2gt13  = 0;
+    PDstandardNth22gt13 = origin ? PDstandardNth11gt13 / 2 : PDstandardNth1gt13 / xx - gt13[index] / (xx * xx);
+    PDstandardNth12gt13 = 0;
+    PDstandardNth23gt13 = 0;
+    PDstandardNth2gt22  = 0;
+    PDstandardNth22gt22 = origin ? PDstandardNth11gt11 : PDstandardNth1gt22 / xx + 2 * (gt11[index] - gt22[index]) / (xx * xx);
+    PDstandardNth12gt22 = 0;
+    PDstandardNth23gt22 = 0;
+    PDstandardNth2gt23  = origin ? PDstandardNth1gt13 : gt13[index] / xx;
+    PDstandardNth22gt23 = 0;
+    PDstandardNth12gt23 = origin ? PDstandardNth11gt13 / 2 : PDstandardNth1gt13 / xx - gt13[index] / (xx * xx);
+    PDstandardNth23gt23 = origin ? PDstandardNth13gt13 : PDstandardNth3gt13 / xx;
+    PDstandardNth2gt33  = 0;
+    PDstandardNth22gt33 = origin ? PDstandardNth11gt33 : PDstandardNth1gt33 / xx;
+    PDstandardNth12gt33 = 0;
+    PDstandardNth23gt33 = 0;
+
+    PDstandardNth2phi  = 0;
+    PDstandardNth22phi = origin ? PDstandardNth11phi : PDstandardNth1phi / xx;
+    PDstandardNth12phi = 0;
+    PDstandardNth23phi = 0;
+    PDstandardNth2trK  = 0;
+
+    PDstandardNth2At11 = origin ? -2 * PDstandardNth1At12 : -2 * At12[index] / xx;
+    PDstandardNth2At12 = origin ? PDstandardNth1At11 - PDstandardNth1At22 : (At11[index] - At22[index]) / xx;
+    PDstandardNth2At13 = origin ? -PDstandardNth1At23 : - At23[index] / xx;
+    PDstandardNth2At22 = origin ? 2 * PDstandardNth1At12 : 2 * At12[index] / xx;
+    PDstandardNth2At23 = origin ? PDstandardNth1At13 : At13[index] / xx;
+    PDstandardNth2At33 = 0;
+    if (fabs(y[index]) > 1e-8)
+        continue;
     
     /* Calculate temporaries and grid functions */
     CCTK_REAL_VEC JacPDstandardNth11gt22 CCTK_ATTRIBUTE_UNUSED;
@@ -2045,43 +2090,6 @@ extern "C" void ML_Kretschmann_kretschmann(CCTK_ARGUMENTS)
     return;
   }
   
-  const char* const groups[] = {
-    "ML_BSSN::ML_curv",
-    "ML_BSSN::ML_log_confac",
-    "ML_BSSN::ML_metric",
-    "ML_BSSN::ML_trace_curv",
-    "ML_Kretschmann::ML_Kretschmann"};
-  GenericFD_AssertGroupStorage(cctkGH, "ML_Kretschmann_kretschmann", 5, groups);
-  
-  switch (fdOrder)
-  {
-    case 2:
-    {
-      GenericFD_EnsureStencilFits(cctkGH, "ML_Kretschmann_kretschmann", 1, 1, 1);
-      break;
-    }
-    
-    case 4:
-    {
-      GenericFD_EnsureStencilFits(cctkGH, "ML_Kretschmann_kretschmann", 2, 2, 2);
-      break;
-    }
-    
-    case 6:
-    {
-      GenericFD_EnsureStencilFits(cctkGH, "ML_Kretschmann_kretschmann", 3, 3, 3);
-      break;
-    }
-    
-    case 8:
-    {
-      GenericFD_EnsureStencilFits(cctkGH, "ML_Kretschmann_kretschmann", 4, 4, 4);
-      break;
-    }
-    default:
-      CCTK_BUILTIN_UNREACHABLE();
-  }
-  
   GenericFD_LoopOverInterior(cctkGH, ML_Kretschmann_kretschmann_Body);
   
   if (verbose > 1)