Vectorize kretschmann

1 files changed, 27 insertions, 21 deletions

diff --git a/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc b/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc
index 6a8f0e6..89461ea 100644
--- a/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc
+++ b/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc
@@ -2084,34 +2084,40 @@ static void ML_Kretschmann_kretschmann_Body(const cGH* restrict const cctkGH, co
     CCTK_REAL_VEC KretschL CCTK_ATTRIBUTE_UNUSED = 
       kadd(term1,kadd(term2,term3));
 
+    CCTK_REAL_VEC zetaL;
     {
-        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));
+        CCTK_REAL_VEC kquart = ToReal(0.25);
+        CCTK_REAL_VEC kfour  = ToReal(4.0);
+        CCTK_REAL_VEC ksix   = ToReal(6.0);
 
-        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;
+        CCTK_REAL_VEC x2 = SQR(xx);
+        CCTK_REAL_VEC x4 = SQR(x2);
+
+        CCTK_REAL_VEC reg_part = SQR(K22) / SQR(g22) - kquart * (SQR(Dg221) * gu11 + gu33 * SQR(Dg223) + ktwo * Dg221 * Dg223 * gu13) / SQR(g22);
+        CCTK_REAL_VEC lambda_reg = (x2 * g22 * (kone - gu11 * g22) - xx * x2 * g22 * (Dg221 * gu11 + gu13 * Dg223)) / (x4 * SQR(g22));
+
+	CCTK_REAL_VEC dphi_x = PDstandardNth1phi;
+	CCTK_REAL_VEC d2phi_xx = PDstandardNth11phi;
+	CCTK_REAL_VEC dgt11_x  = PDstandardNth1gt11;
+	CCTK_REAL_VEC d2gt11_xx = PDstandardNth11gt11;
+	CCTK_REAL_VEC Dgu131 = -Dg131 * (gu11 * gu33 + SQR(gu13));
+	CCTK_REAL_VEC Dgu111 = kzero;
+	CCTK_REAL_VEC D2g1111 = - ksix * SQR(dphi_x) * gt11L / SQR(SQR(phiL))
+			 - ktwo * d2phi_xx * gt11L / (phiL * SQR(phiL))
+			 - kfour * dphi_x * dgt11_x / (phiL * SQR(phiL))
+			 + d2gt11_xx / SQR(phiL);
+	CCTK_REAL_VEC D2gu1111 = -(SQR(gu11) * D2g1111 + ktwo * Dg131 * (Dgu111 * gu13 + gu11 * Dgu131) + SQR(gu13) * D2g3311);
+	CCTK_REAL_VEC d2x_guxx_gyy = D2gu1111 * g22 + gu11 * D2g2211;
+        CCTK_REAL_VEC lambda_sing = -(gu11 * D2g2211 + Dgu131 * Dg223) / g22 - khalf * d2x_guxx_gyy / g22;
+        CCTK_REAL_VEC lambda = kifthen(origin, lambda_sing, lambda_reg);
+
+        zetaL = reg_part + lambda;
     }
     
     /* Copy local copies back to grid functions */
     vec_store_partial_prepare(i,vecimin,vecimax);
     vec_store_nta_partial(Kretsch[index],KretschL);
+    vec_store_nta_partial(zeta[index], zetaL);
   }
   CCTK_ENDLOOP3STR(ML_Kretschmann_kretschmann);
 }