1 files changed, 61 insertions, 43 deletions
diff --git a/ML_BSSN/src/ML_BSSN_constraints.cc b/ML_BSSN/src/ML_BSSN_constraints.cc
index 62f8953..50c2e48 100644
--- a/ML_BSSN/src/ML_BSSN_constraints.cc
+++ b/ML_BSSN/src/ML_BSSN_constraints.cc
@@ -240,8 +240,6 @@ static void ML_BSSN_constraints_Body(const cGH* restrict const cctkGH, const int
   {
     const ptrdiff_t index CCTK_ATTRIBUTE_UNUSED = di*i + dj*j + dk*k;
     // vec_iter_counter+=CCTK_REAL_VEC_SIZE;
-    if (fabs(y[index]) > 1e-8)
-        continue;
     
     /* Assign local copies of grid functions */
     
@@ -822,51 +820,71 @@ static void ML_BSSN_constraints_Body(const cGH* restrict const cctkGH, const int
       default:
         CCTK_BUILTIN_UNREACHABLE();
     }
-    
-    double xx = x[index];
-    int origin = fabs(xx) < 1e-8;
+    CCTK_REAL_VEC kzero = ToReal(0.0);
+    CCTK_REAL_VEC kone = ToReal(1.0);
+    CCTK_REAL_VEC ktwo = ToReal(2.0);
+    CCTK_REAL_VEC xx = vec_load(x[index]);
+    CCTK_REAL_VEC x2 = kmul(xx, xx);
+
+    CCTK_REAL_VEC xinv  = kdiv(kone, xx);
+    CCTK_REAL_VEC x2inv = kdiv(kone, x2);
 
-    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;
+    CCTK_REAL_VEC absx = kfabs(xx);
+    CCTK_REAL_VEC eps  = ToReal(1e-8);
+    CCTK_BOOLEAN_VEC origin = kcmplt(absx, eps);
 
-    PDstandardNth2phi = 0;
-    PDstandardNth22phi = origin ? PDstandardNth11phi : PDstandardNth1phi / xx;
-    PDstandardNth12phi = 0;
-    PDstandardNth23phi = 0;
-    PDstandardNth2trK = 0;
-    PDstandardNth2Xt1 = 0;
-    PDstandardNth2Xt2 = origin ? PDstandardNth1Xt1 : Xt1[index] / xx;
-    PDstandardNth2Xt3 = 0;
+    PDstandardNth2gt11  = kzero;
+    PDstandardNth22gt11 = kifthen(origin, PDstandardNth11gt22,
+                                          ksub(kmul(PDstandardNth1gt11, xinv), kmul(ktwo, kmul(ksub(gt11L, gt22L), x2inv))));
+    PDstandardNth12gt11 = kzero;
+    PDstandardNth23gt11 = kzero;
+    PDstandardNth2gt12  = kifthen(origin, ksub(PDstandardNth1gt11, PDstandardNth1gt22),
+                                          kmul(ksub(gt11L, gt22L), xinv));
+    PDstandardNth22gt12 = kzero;
+    PDstandardNth12gt12 = kifthen(origin, kmul(ksub(PDstandardNth11gt11, PDstandardNth11gt22), khalf),
+                                          ksub(kmul(ksub(PDstandardNth1gt11, PDstandardNth1gt22), xinv), kmul(ksub(gt11L, gt22L), x2inv)));
+    PDstandardNth23gt12 = kifthen(origin, ksub(PDstandardNth13gt11, PDstandardNth13gt22),
+                                          kmul(ksub(PDstandardNth3gt11, PDstandardNth3gt22), xinv));
+    PDstandardNth2gt13  = kzero;
+    PDstandardNth22gt13 = kifthen(origin, kmul(PDstandardNth11gt13, khalf),
+                                          ksub(kmul(PDstandardNth1gt13, xinv), kmul(gt13L, x2inv)));
+    PDstandardNth12gt13 = kzero;
+    PDstandardNth23gt13 = kzero;
+    PDstandardNth2gt22  = kzero;
+    PDstandardNth22gt22 = kifthen(origin, PDstandardNth11gt11,
+                                          kadd(kmul(PDstandardNth1gt22, xinv), kmul(ktwo, kmul(ksub(gt11L, gt22L), x2inv))));
+    PDstandardNth12gt22 = kzero;
+    PDstandardNth23gt22 = kzero;
+    PDstandardNth2gt23  = kifthen(origin, PDstandardNth1gt13, kmul(gt13L, xinv));
+    PDstandardNth22gt23 = kzero;
+    PDstandardNth12gt23 = kifthen(origin, kmul(PDstandardNth11gt13, khalf),
+                                          ksub(kmul(PDstandardNth1gt13, xinv), kmul(gt13L, x2inv)));
+    PDstandardNth23gt23 = kifthen(origin, PDstandardNth13gt13, kmul(PDstandardNth3gt13, xinv));
+    PDstandardNth2gt33  = kzero;
+    PDstandardNth22gt33 = kifthen(origin, PDstandardNth11gt33, kmul(PDstandardNth1gt33, xinv));
+    PDstandardNth12gt33 = kzero;
+    PDstandardNth23gt33 = kzero;
 
-    PDstandardNth2At11 = 0; //origin ? -2 * PDstandardNth1At12 : -2 * At12[index] / xx;
-    PDstandardNth2At12 = origin ? PDstandardNth1At11 - PDstandardNth1At22 : (At11[index] - At22[index]) / xx;
-    PDstandardNth2At13 = 0; //origin ? -PDstandardNth1At23 : - At23[index] / xx;
-    PDstandardNth2At22 = 0; //origin ? 2 * PDstandardNth1At12 : 2 * At12[index] / xx;
-    PDstandardNth2At23 = origin ? PDstandardNth1At13 : At13[index] / xx;
-    PDstandardNth2At33 = 0;
+    PDstandardNth2phi  = kzero;
+    PDstandardNth22phi = kifthen(origin, PDstandardNth11phi, kmul(PDstandardNth1phi, xinv));
+    PDstandardNth12phi = kzero;
+    PDstandardNth23phi = kzero;
+    PDstandardNth2Xt1  = kzero;
+    PDstandardNth2Xt2  = kifthen(origin, PDstandardNth1Xt1, kmul(Xt1L, xinv));
+    PDstandardNth2Xt3  = kzero;
+    PDstandardNth2trK  = kzero;
 
+    //PDstandardNth2At11 = origin ? -2 * PDstandardNth1At12 : -2 * At12[index] / xx;
+    //PDstandardNth2At13 = origin ? -PDstandardNth1At23 : - At23[index] / xx;
+    //PDstandardNth2At22 = origin ? 2 * PDstandardNth1At12 : 2 * At12[index] / xx;
+    PDstandardNth2At11 = kzero;
+    PDstandardNth2At12 = kifthen(origin, ksub(PDstandardNth1At11, PDstandardNth1At22),
+                                         kmul(ksub(At11L, At22L), xinv));
+    PDstandardNth2At13 = kzero;
+    PDstandardNth2At22 = kzero;
+    PDstandardNth2At23 = kifthen(origin, PDstandardNth1At13, kmul(At13L, xinv));
+    PDstandardNth2At33 = kzero;
+    
     /* Calculate temporaries and grid functions */
     CCTK_REAL_VEC JacPDstandardNth11gt11 CCTK_ATTRIBUTE_UNUSED;
     CCTK_REAL_VEC JacPDstandardNth11gt12 CCTK_ATTRIBUTE_UNUSED;