Vectorize

8 files changed, 459 insertions, 336 deletions

diff --git a/ML_BSSN/src/ML_BSSN_Advect.cc b/ML_BSSN/src/ML_BSSN_Advect.cc
index e03d073..c9dd430 100644
--- a/ML_BSSN/src/ML_BSSN_Advect.cc
+++ b/ML_BSSN/src/ML_BSSN_Advect.cc
@@ -1033,57 +1033,67 @@ static void ML_BSSN_Advect_Body(const cGH* restrict const cctkGH, const int dir,
         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);
 
-    PDupwindNthAnti2At11 = 0;
-    PDupwindNthSymm2At11 = 0;
-    PDupwindNthAnti2At12 = origin ? PDupwindNthAnti1At11 - PDupwindNthAnti1At22 : (At11[index] - At22[index]) / xx;
-    PDupwindNthSymm2At12 = origin ? PDupwindNthSymm1At11 - PDupwindNthSymm1At22 : (At11[index] - At22[index]) / xx;
-    PDupwindNthAnti2At13 = 0;
-    PDupwindNthSymm2At13 = 0;
-    PDupwindNthAnti2At22 = 0;
-    PDupwindNthSymm2At22 = 0;
-    PDupwindNthAnti2At23 = origin ? PDupwindNthAnti1At13 : At13[index] / xx;
-    PDupwindNthSymm2At23 = origin ? PDupwindNthSymm1At13 : At13[index] / xx;
-    PDupwindNthAnti2At33 = 0;
-    PDupwindNthSymm2At33 = 0;
+    CCTK_REAL_VEC xinv  = kdiv(kone, xx);
+    CCTK_REAL_VEC x2inv = kdiv(kone, x2);
 
-    PDupwindNthAnti2B1 = 0;
-    PDupwindNthSymm2B1 = 0;
-    PDupwindNthAnti2B2 = origin ? PDupwindNthAnti1B1 : B1[index] / xx;
-    PDupwindNthSymm2B2 = origin ? PDupwindNthSymm1B1 : B1[index] / xx;
-    PDupwindNthAnti2B3 = 0;
-    PDupwindNthSymm2B3 = 0;
+    CCTK_REAL_VEC absx = kfabs(xx);
+    CCTK_REAL_VEC eps  = ToReal(1e-8);
+    CCTK_BOOLEAN_VEC origin = kcmplt(absx, eps);
 
-    PDupwindNthAnti2beta1 = 0;
-    PDupwindNthSymm2beta1 = 0;
-    PDupwindNthAnti2beta2 = origin ? PDupwindNthAnti1beta1 : beta1[index] / xx;
-    PDupwindNthSymm2beta2 = origin ? PDupwindNthSymm1beta1 : beta1[index] / xx;
-    PDupwindNthAnti2beta3 = 0;
-    PDupwindNthSymm2beta3 = 0;
+    PDupwindNthAnti2At11 = kzero;
+    PDupwindNthSymm2At11 = kzero;
+    PDupwindNthAnti2At12 = kifthen(origin, ksub(PDupwindNthAnti1At11, PDupwindNthAnti1At22), kmul(ksub(At11L, At22L), xinv));
+    PDupwindNthSymm2At12 = kifthen(origin, ksub(PDupwindNthSymm1At11, PDupwindNthSymm1At22), kmul(ksub(At11L, At22L), xinv));
+    PDupwindNthAnti2At13 = kzero;
+    PDupwindNthSymm2At13 = kzero;
+    PDupwindNthAnti2At22 = kzero;
+    PDupwindNthSymm2At22 = kzero;
+    PDupwindNthAnti2At23 = kifthen(origin, PDupwindNthAnti1At13, kmul(At13L, xinv));
+    PDupwindNthSymm2At23 = kifthen(origin, PDupwindNthSymm1At13, kmul(At13L, xinv));
+    PDupwindNthAnti2At33 = kzero;
+    PDupwindNthSymm2At33 = kzero;
 
-    PDupwindNthAnti2gt11 = 0;
-    PDupwindNthSymm2gt11 = 0;
-    PDupwindNthAnti2gt12 = origin ? PDupwindNthAnti1gt11 - PDupwindNthAnti1gt22 : (gt11[index] - gt22[index]) / xx;
-    PDupwindNthSymm2gt12 = origin ? PDupwindNthSymm1gt11 - PDupwindNthSymm1gt22 : (gt11[index] - gt22[index]) / xx;
-    PDupwindNthAnti2gt13 = 0;
-    PDupwindNthSymm2gt13 = 0;
-    PDupwindNthAnti2gt22 = 0;
-    PDupwindNthSymm2gt22 = 0;
-    PDupwindNthAnti2gt23 = origin ? PDupwindNthAnti1gt13 : gt13[index] / xx;
-    PDupwindNthSymm2gt23 = origin ? PDupwindNthSymm1gt13 : gt13[index] / xx;
-    PDupwindNthAnti2gt33 = 0;
-    PDupwindNthSymm2gt33 = 0;
+    PDupwindNthAnti2B1 = kzero;
+    PDupwindNthSymm2B1 = kzero;
+    PDupwindNthAnti2B2 = kifthen(origin, PDupwindNthAnti1B1, kmul(B1L, xinv));
+    PDupwindNthSymm2B2 = kifthen(origin, PDupwindNthSymm1B1, kmul(B1L, xinv));
+    PDupwindNthAnti2B3 = kzero;
+    PDupwindNthSymm2B3 = kzero;
+
+    PDupwindNthAnti2beta1 = kzero;
+    PDupwindNthSymm2beta1 = kzero;
+    PDupwindNthAnti2beta2 = kifthen(origin, PDupwindNthAnti1beta1, kmul(beta1L, xinv));
+    PDupwindNthSymm2beta2 = kifthen(origin, PDupwindNthSymm1beta1, kmul(beta1L, xinv));
+    PDupwindNthAnti2beta3 = kzero;
+    PDupwindNthSymm2beta3 = kzero;
+
+    PDupwindNthAnti2gt11 = kzero;
+    PDupwindNthSymm2gt11 = kzero;
+    PDupwindNthAnti2gt12 = kifthen(origin, ksub(PDupwindNthAnti1gt11, PDupwindNthAnti1gt22), kmul(ksub(gt11L, gt22L), xinv));
+    PDupwindNthSymm2gt12 = kifthen(origin, ksub(PDupwindNthSymm1gt11, PDupwindNthSymm1gt22), kmul(ksub(gt11L, gt22L), xinv));
+    PDupwindNthAnti2gt13 = kzero;
+    PDupwindNthSymm2gt13 = kzero;
+    PDupwindNthAnti2gt22 = kzero;
+    PDupwindNthSymm2gt22 = kzero;
+    PDupwindNthAnti2gt23 = kifthen(origin, PDupwindNthAnti1gt13, kmul(gt13L, xinv));
+    PDupwindNthSymm2gt23 = kifthen(origin, PDupwindNthSymm1gt13, kmul(gt13L, xinv));
+    PDupwindNthAnti2gt33 = kzero;
+    PDupwindNthSymm2gt33 = kzero;
+
+    PDupwindNthAnti2phi = kzero;
+    PDupwindNthSymm2phi = kzero;
+    PDupwindNthAnti2trK = kzero;
+    PDupwindNthSymm2trK = kzero;
+    PDupwindNthAnti2Xt2 = kifthen(origin, PDupwindNthAnti1Xt1, kmul(Xt1L, xinv));
+    PDupwindNthSymm2Xt2 = kifthen(origin, PDupwindNthSymm1Xt1, kmul(Xt1L, xinv));
 
-    PDupwindNthAnti2phi = 0;
-    PDupwindNthSymm2phi = 0;
-    PDupwindNthAnti2trK = 0;
-    PDupwindNthSymm2trK = 0;
-    PDupwindNthAnti2Xt2 = origin ? PDupwindNthAnti1Xt1 : Xt1[index] / xx;
-    PDupwindNthSymm2Xt2 = origin ? PDupwindNthSymm1Xt1 : Xt1[index] / xx;
 
-    
     /* Calculate temporaries and grid functions */
     ptrdiff_t dir1 CCTK_ATTRIBUTE_UNUSED = kisgn(beta1L);
     
diff --git a/ML_BSSN/src/ML_BSSN_RHS.cc b/ML_BSSN/src/ML_BSSN_RHS.cc
index 098f585..a3bb92e 100644
--- a/ML_BSSN/src/ML_BSSN_RHS.cc
+++ b/ML_BSSN/src/ML_BSSN_RHS.cc
@@ -946,60 +946,79 @@ static void ML_BSSN_RHS_Body(const cGH* restrict const cctkGH, const int dir, co
         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);
 
-    PDstandardNth2trK    = 0;
-    PDstandardNth2alpha  = 0;
-    PDstandardNth22alpha = origin ? PDstandardNth11alpha : PDstandardNth1alpha / xx;
-    PDstandardNth12alpha = 0;
-    PDstandardNth23alpha = 0;
+    CCTK_REAL_VEC xinv  = kdiv(kone, xx);
+    CCTK_REAL_VEC x2inv = kdiv(kone, x2);
 
-    PDstandardNth2beta1  = 0;
-    PDstandardNth22beta1 = origin ? PDstandardNth11beta1 / 2 : PDstandardNth1beta1 / xx - beta1[index] / (xx * xx);
-    PDstandardNth12beta1 = 0;
-    PDstandardNth23beta1 = 0;
-    PDstandardNth2beta2  = origin ? PDstandardNth1beta1 : beta1[index] / xx;
-    PDstandardNth22beta2 = 0;
-    PDstandardNth12beta2 = origin ? PDstandardNth11beta1 / 2 : PDstandardNth1beta1 / xx - beta1[index] / (xx * xx);
-    PDstandardNth23beta2 = origin ? PDstandardNth13beta1 : PDstandardNth3beta1 / xx;
-    PDstandardNth2beta3  = 0;
-    PDstandardNth22beta3 = origin ? PDstandardNth11beta3 : PDstandardNth1beta3 / xx;
-    PDstandardNth12beta3 = 0;
-    PDstandardNth23beta3 = 0;
+    CCTK_REAL_VEC absx = kfabs(xx);
+    CCTK_REAL_VEC eps  = ToReal(1e-8);
+    CCTK_BOOLEAN_VEC origin = kcmplt(absx, eps);
 
-    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;
+    PDstandardNth2trK    = kzero;
+    PDstandardNth2alpha  = kzero;
+    PDstandardNth22alpha = kifthen(origin, PDstandardNth11alpha, kmul(PDstandardNth1alpha, xinv));
+    PDstandardNth12alpha = kzero;
+    PDstandardNth23alpha = kzero;
 
-    PDstandardNth2phi  = 0;
-    PDstandardNth22phi = origin ? PDstandardNth11phi : PDstandardNth1phi / xx;
-    PDstandardNth12phi = 0;
-    PDstandardNth23phi = 0;
-    PDstandardNth2Xt1  = 0;
-    PDstandardNth2Xt2  = origin ? PDstandardNth1Xt1 : Xt1[index] / xx;
-    PDstandardNth2Xt3  = 0;
+    PDstandardNth2beta1  = kzero;
+    PDstandardNth22beta1 = kifthen(origin, kmul(PDstandardNth11beta1, khalf),
+                                           ksub(kmul(PDstandardNth1beta1, xinv), kmul(beta1L, x2inv)));
+    PDstandardNth12beta1 = kzero;
+    PDstandardNth23beta1 = kzero;
+    PDstandardNth2beta2  = kifthen(origin, PDstandardNth1beta1, kmul(beta1L, xinv));
+    PDstandardNth22beta2 = kzero;
+    PDstandardNth12beta2 = kifthen(origin, kmul(PDstandardNth11beta1, khalf),
+                                           ksub(kmul(PDstandardNth1beta1, xinv), kmul(beta1L, x2inv)));
+    PDstandardNth23beta2 = kifthen(origin, PDstandardNth13beta1, kmul(PDstandardNth3beta1, xinv));
+    PDstandardNth2beta3  = kzero;
+    PDstandardNth22beta3 = kifthen(origin, PDstandardNth11beta3, kmul(PDstandardNth1beta3, xinv));
+    PDstandardNth12beta3 = kzero;
+    PDstandardNth23beta3 = kzero;
+
+    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;
+
+    PDstandardNth2phi  = kzero;
+    PDstandardNth22phi = kifthen(origin, PDstandardNth11phi, kmul(PDstandardNth1phi, xinv));
+    PDstandardNth12phi = kzero;
+    PDstandardNth23phi = kzero;
+    PDstandardNth2Xt1  = kzero;
+    PDstandardNth2Xt2  = kifthen(origin, PDstandardNth1Xt1, kmul(Xt1L, xinv));
+    PDstandardNth2Xt3  = kzero;
     /* Calculate temporaries and grid functions */
     ptrdiff_t dir1 CCTK_ATTRIBUTE_UNUSED = kisgn(beta1L);
     
diff --git a/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.cc b/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.cc
index 693b384..eff3b7f 100644
--- a/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.cc
+++ b/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.cc
@@ -473,15 +473,26 @@ static void ML_BSSN_convertFromADMBaseGamma_Body(const cGH* restrict const cctkG
         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);
 
-    PDstandardNth2gt11  = 0;
-    PDstandardNth2gt12  = origin ? PDstandardNth1gt11 - PDstandardNth1gt22 : (gt11[index] - gt22[index]) / xx;
-    PDstandardNth2gt13  = 0;
-    PDstandardNth2gt22  = 0;
-    PDstandardNth2gt23  = origin ? PDstandardNth1gt13 : gt13[index] / xx;
-    PDstandardNth2gt33  = 0;
+    CCTK_REAL_VEC xinv  = kdiv(kone, xx);
+    CCTK_REAL_VEC x2inv = kdiv(kone, x2);
+
+    CCTK_REAL_VEC absx = kfabs(xx);
+    CCTK_REAL_VEC eps  = ToReal(1e-8);
+    CCTK_BOOLEAN_VEC origin = kcmplt(absx, eps);
+
+    PDstandardNth2gt11  = kzero;
+    PDstandardNth2gt12  = kifthen(origin, ksub(PDstandardNth1gt11, PDstandardNth1gt22),
+                                          kmul(ksub(gt11L, gt22L), xinv));
+    PDstandardNth2gt13  = kzero;
+    PDstandardNth2gt22  = kzero;
+    PDstandardNth2gt23  = kifthen(origin, PDstandardNth1gt13, kmul(gt13L, xinv));
+    PDstandardNth2gt33  = kzero;
     
     /* Calculate temporaries and grid functions */
     ptrdiff_t dir1 CCTK_ATTRIBUTE_UNUSED = kisgn(beta1L);
diff --git a/ML_BSSN/src/eval_timediff.cc b/ML_BSSN/src/eval_timediff.cc
index 46f7599..ec35923 100644
--- a/ML_BSSN/src/eval_timediff.cc
+++ b/ML_BSSN/src/eval_timediff.cc
@@ -867,60 +867,79 @@ static void eval_timediff_body(const cGH* restrict const cctkGH, const int dir,
         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);
 
-    PDstandardNth2trK    = 0;
-    PDstandardNth2alpha  = 0;
-    PDstandardNth22alpha = origin ? PDstandardNth11alpha : PDstandardNth1alpha / xx;
-    PDstandardNth12alpha = 0;
-    PDstandardNth23alpha = 0;
+    CCTK_REAL_VEC xinv  = kdiv(kone, xx);
+    CCTK_REAL_VEC x2inv = kdiv(kone, x2);
 
-    PDstandardNth2beta1  = 0;
-    PDstandardNth22beta1 = origin ? PDstandardNth11beta1 / 2 : PDstandardNth1beta1 / xx - beta1[index] / (xx * xx);
-    PDstandardNth12beta1 = 0;
-    PDstandardNth23beta1 = 0;
-    PDstandardNth2beta2  = origin ? PDstandardNth1beta1 : beta1[index] / xx;
-    PDstandardNth22beta2 = 0;
-    PDstandardNth12beta2 = origin ? PDstandardNth11beta1 / 2 : PDstandardNth1beta1 / xx - beta1[index] / (xx * xx);
-    PDstandardNth23beta2 = origin ? PDstandardNth13beta1 : PDstandardNth3beta1 / xx;
-    PDstandardNth2beta3  = 0;
-    PDstandardNth22beta3 = origin ? PDstandardNth11beta3 : PDstandardNth1beta3 / xx;
-    PDstandardNth12beta3 = 0;
-    PDstandardNth23beta3 = 0;
+    CCTK_REAL_VEC absx = kfabs(xx);
+    CCTK_REAL_VEC eps  = ToReal(1e-8);
+    CCTK_BOOLEAN_VEC origin = kcmplt(absx, eps);
 
-    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;
+    PDstandardNth2trK    = kzero;
+    PDstandardNth2alpha  = kzero;
+    PDstandardNth22alpha = kifthen(origin, PDstandardNth11alpha, kmul(PDstandardNth1alpha, xinv));
+    PDstandardNth12alpha = kzero;
+    PDstandardNth23alpha = kzero;
 
-    PDstandardNth2phi  = 0;
-    PDstandardNth22phi = origin ? PDstandardNth11phi : PDstandardNth1phi / xx;
-    PDstandardNth12phi = 0;
-    PDstandardNth23phi = 0;
-    PDstandardNth2Xt1  = 0;
-    PDstandardNth2Xt2  = origin ? PDstandardNth1Xt1 : Xt1[index] / xx;
-    PDstandardNth2Xt3  = 0;
+    PDstandardNth2beta1  = kzero;
+    PDstandardNth22beta1 = kifthen(origin, kmul(PDstandardNth11beta1, khalf),
+                                           ksub(kmul(PDstandardNth1beta1, xinv), kmul(beta1L, x2inv)));
+    PDstandardNth12beta1 = kzero;
+    PDstandardNth23beta1 = kzero;
+    PDstandardNth2beta2  = kifthen(origin, PDstandardNth1beta1, kmul(beta1L, xinv));
+    PDstandardNth22beta2 = kzero;
+    PDstandardNth12beta2 = kifthen(origin, kmul(PDstandardNth11beta1, khalf),
+                                           ksub(kmul(PDstandardNth1beta1, xinv), kmul(beta1L, x2inv)));
+    PDstandardNth23beta2 = kifthen(origin, PDstandardNth13beta1, kmul(PDstandardNth3beta1, xinv));
+    PDstandardNth2beta3  = kzero;
+    PDstandardNth22beta3 = kifthen(origin, PDstandardNth11beta3, kmul(PDstandardNth1beta3, xinv));
+    PDstandardNth12beta3 = kzero;
+    PDstandardNth23beta3 = kzero;
+
+    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;
+
+    PDstandardNth2phi  = kzero;
+    PDstandardNth22phi = kifthen(origin, PDstandardNth11phi, kmul(PDstandardNth1phi, xinv));
+    PDstandardNth12phi = kzero;
+    PDstandardNth23phi = kzero;
+    PDstandardNth2Xt1  = kzero;
+    PDstandardNth2Xt2  = kifthen(origin, PDstandardNth1Xt1, kmul(Xt1L, xinv));
+    PDstandardNth2Xt3  = kzero;
     /* Calculate temporaries and grid functions */
     
     CCTK_REAL_VEC JacPDstandardNth11alpha CCTK_ATTRIBUTE_UNUSED;
diff --git a/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc b/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc
index a43da4c..6a8f0e6 100644
--- a/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc
+++ b/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc
@@ -690,50 +690,64 @@ static void ML_Kretschmann_kretschmann_Body(const cGH* restrict const cctkGH, co
         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);
 
-    PDstandardNth2trK    = 0;
+    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;
+    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 = 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;
+    PDstandardNth2phi  = kzero;
+    PDstandardNth22phi = kifthen(origin, PDstandardNth11phi, kmul(PDstandardNth1phi, xinv));
+    PDstandardNth12phi = kzero;
+    PDstandardNth23phi = kzero;
+    PDstandardNth2trK    = kzero;
+
+    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 JacPDstandardNth11gt22 CCTK_ATTRIBUTE_UNUSED;
diff --git a/cartoon.c b/cartoon.c
index 13e1b71..9f3c03b 100644
--- a/cartoon.c
+++ b/cartoon.c
@@ -1,64 +1,84 @@
-    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);
 
-    PDstandardNth2trK    = 0;
-    PDstandardNth2alp  = 0;
-    PDstandardNth22alp = origin ? PDstandardNth11alp : PDstandardNth1alp / xx;
-    PDstandardNth12alp = 0;
-    PDstandardNth23alp = 0;
+    CCTK_REAL_VEC xinv  = kdiv(kone, xx);
+    CCTK_REAL_VEC x2inv = kdiv(kone, x2);
 
-    PDstandardNth2beta1  = 0;
-    PDstandardNth22beta1 = origin ? PDstandardNth11beta1 / 2 : PDstandardNth1beta1 / xx - beta1[index] / (xx * xx);
-    PDstandardNth12beta1 = 0;
-    PDstandardNth23beta1 = 0;
-    PDstandardNth2beta2  = origin ? PDstandardNth1beta1 : beta1[index] / xx;
-    PDstandardNth22beta2 = 0;
-    PDstandardNth12beta2 = origin ? PDstandardNth11beta1 / 2 : PDstandardNth1beta1 / xx - beta1[index] / (xx * xx);
-    PDstandardNth23beta2 = origin ? PDstandardNth13beta1 : PDstandardNth3beta1 / xx;
-    PDstandardNth2beta3  = 0;
-    PDstandardNth22beta3 = origin ? PDstandardNth11beta3 : PDstandardNth1beta3 / xx;
-    PDstandardNth12beta3 = 0;
-    PDstandardNth23beta3 = 0;
+    CCTK_REAL_VEC absx = kfabs(xx);
+    CCTK_REAL_VEC eps  = ToReal(1e-8);
+    CCTK_BOOLEAN_VEC origin = kcmplt(absx, eps);
 
-    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;
+    PDstandardNth2trK    = kzero;
+    PDstandardNth2alpha  = kzero;
+    PDstandardNth22alpha = kifthen(origin, PDstandardNth11alpha, kmul(PDstandardNth1alpha, xinv));
+    PDstandardNth12alpha = kzero;
+    PDstandardNth23alpha = kzero;
 
-    PDstandardNth2phi  = 0;
-    PDstandardNth22phi = origin ? PDstandardNth11phi : PDstandardNth1phi / xx;
-    PDstandardNth12phi = 0;
-    PDstandardNth23phi = 0;
-    PDstandardNth2trK  = 0;
-    PDstandardNth2Xt1  = 0;
-    PDstandardNth2Xt2  = origin ? PDstandardNth1Xt1 : Xt1[index] / xx;
-    PDstandardNth2Xt3  = 0;
+    PDstandardNth2beta1  = kzero;
+    PDstandardNth22beta1 = kifthen(origin, kmul(PDstandardNth11beta1, khalf),
+                                           ksub(kmul(PDstandardNth1beta1, xinv), kmul(beta1L, x2inv)));
+    PDstandardNth12beta1 = kzero;
+    PDstandardNth23beta1 = kzero;
+    PDstandardNth2beta2  = kifthen(origin, PDstandardNth1beta1, kmul(beta1L, xinv));
+    PDstandardNth22beta2 = kzero;
+    PDstandardNth12beta2 = kifthen(origin, kmul(PDstandardNth11beta1, khalf),
+                                           ksub(kmul(PDstandardNth1beta1, xinv), kmul(beta1L, x2inv)));
+    PDstandardNth23beta2 = kifthen(origin, PDstandardNth13beta1, kmul(PDstandardNth3beta1, xinv));
+    PDstandardNth2beta3  = kzero;
+    PDstandardNth22beta3 = kifthen(origin, PDstandardNth11beta3, kmul(PDstandardNth1beta3, xinv));
+    PDstandardNth12beta3 = kzero;
+    PDstandardNth23beta3 = kzero;
 
-    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;
+    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;
+
+    PDstandardNth2phi  = kzero;
+    PDstandardNth22phi = kifthen(origin, PDstandardNth11phi, kmul(PDstandardNth1phi, xinv));
+    PDstandardNth12phi = kzero;
+    PDstandardNth23phi = kzero;
+    PDstandardNth2Xt1  = kzero;
+    PDstandardNth2Xt2  = kifthen(origin, PDstandardNth1Xt1, kmul(Xt1L, xinv));
+    PDstandardNth2Xt3  = 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;
diff --git a/cartoon_advect.c b/cartoon_advect.c
index 91da4da..74ff90a 100644
--- a/cartoon_advect.c
+++ b/cartoon_advect.c
@@ -1,50 +1,59 @@
-    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);
 
-    PDupwindNthAnti2At11 = origin ? -2 * PDupwindNthAnti1At12 : -2 * At12[index] / xx;
-    PDupwindNthSymm2At11 = origin ? -2 * PDupwindNthSymm1At12 : -2 * At12[index] / xx;
-    PDupwindNthAnti2At12 = origin ? PDupwindNthAnti1At11 - PDupwindNthAntiAt22 : (At11[index] - At22[index]) / xx;
-    PDupwindNthSymm2At12 = origin ? PDupwindNthSymm1At11 - PDupwindNthSymmAt22 : (At11[index] - At22[index]) / xx;
-    PDupwindNthAnti2At13 = origin ? -PDupwindNthAnti1At23 : - At23[index] / xx;
-    PDupwindNthSymm2At13 = origin ? -PDupwindNthSymm1At23 : - At23[index] / xx;
-    PDupwindNthAnti2At22 = origin ? 2 * PDupwindNthAnti1At12 : 2 * At12[index] / xx;
-    PDupwindNthSymm2At22 = origin ? 2 * PDupwindNthSymm1At12 : 2 * At12[index] / xx;
-    PDupwindNthAnti2At23 = origin ? PDupwindNth1AntiAt13 : At13[index] / xx;
-    PDupwindNthSymm2At23 = origin ? PDupwindNth1SymmAt13 : At13[index] / xx;
-    PDupwindNthAnti2At33 = 0;
-    PDupwindNthSymm2At33 = 0;
+    CCTK_REAL_VEC xinv  = kdiv(kone, xx);
+    CCTK_REAL_VEC x2inv = kdiv(kone, x2);
 
-    PDupwindNthAnti2B1 = 0;
-    PDupwindNthSymm2B1 = 0;
-    PDupwindNthAnti2B2 = origin ? PDstandardNth1B1 : B1[index] / xx;
-    PDupwindNthSymm2B2 = origin ? PDstandardNth1B1 : B1[index] / xx;
-    PDupwindNthAnti2B3 = 0;
-    PDupwindNthSymm2B3 = 0;
+    CCTK_REAL_VEC absx = kfabs(xx);
+    CCTK_REAL_VEC eps  = ToReal(1e-8);
+    CCTK_BOOLEAN_VEC origin = kcmplt(absx, eps);
 
-    PDupwindNthAnti2beta1 = 0;
-    PDupwindNthSymm2beta1 = 0;
-    PDupwindNthAnti2beta2 = PDupwindNthAntifdOrder82(&beta2[index]);
-    PDupwindNthSymm2beta2 = PDupwindNthSymmfdOrder82(&beta2[index]);
-    PDupwindNthAnti2beta3 = 0;
-    PDupwindNthSymm2beta3 = 0;
+    PDupwindNthAnti2At11 = kzero;
+    PDupwindNthSymm2At11 = kzero;
+    PDupwindNthAnti2At12 = kifthen(origin, ksub(PDupwindNthAnti1At11, PDupwindNthAntiAt22), kmul(ksub(At11L, At22L), xinv));
+    PDupwindNthSymm2At12 = kifthen(origin, ksub(PDupwindNthSymm1At11, PDupwindNthSymmAt22), kmul(ksub(At11L, At22L), xinv));
+    PDupwindNthAnti2At13 = kzero;
+    PDupwindNthSymm2At13 = kzero;
+    PDupwindNthAnti2At22 = kzero;
+    PDupwindNthSymm2At22 = kzero;
+    PDupwindNthAnti2At23 = kifthen(origin, PDupwindNth1AntiAt13, kmul(At13L, xinv));
+    PDupwindNthSymm2At23 = kifthen(origin, PDupwindNth1SymmAt13, kmul(At13L, xinv));
+    PDupwindNthAnti2At33 = kzero;
+    PDupwindNthSymm2At33 = kzero;
 
-    PDupwindNthAnti2gt11 = 0;
-    PDupwindNthSymm2gt11 = 0;
-    PDupwindNthAnti2gt12 = origin ? PDstandardNthAnti1gt11 - PDstandardNthAnti1gt22 : (gt11[index] - gt22[index]) / xx;
-    PDupwindNthSymm2gt12 = origin ? PDstandardNthSymm1gt11 - PDstandardNthSymm1gt22 : (gt11[index] - gt22[index]) / xx;
-    PDupwindNthAnti2gt13 = 0;
-    PDupwindNthSymm2gt13 = 0;
-    PDupwindNthAnti2gt22 = 0;
-    PDupwindNthSymm2gt22 = 0;
-    PDupwindNthAnti2gt23 = origin ? PDstandardNthAnti1gt13 : gt13[index] / xx;
-    PDupwindNthSymm2gt23 = origin ? PDstandardNthSymm1gt13 : gt13[index] / xx;
-    PDupwindNthAnti2gt33 = 0;
-    PDupwindNthSymm2gt33 = 0;
+    PDupwindNthAnti2B1 = kzero;
+    PDupwindNthSymm2B1 = kzero;
+    PDupwindNthAnti2B2 = kifthen(origin, PDstandardNthAnti1B1, kmul(B1L, xinv));
+    PDupwindNthSymm2B2 = kifthen(origin, PDstandardNthSymm1B1, kmul(B1L, xinv));
+    PDupwindNthAnti2B3 = kzero;
+    PDupwindNthSymm2B3 = kzero;
 
-    PDupwindNthAnti2phi = 0;
-    PDupwindNthSymm2phi = 0;
-    PDupwindNthAnti2trK = 0;
-    PDupwindNthSymm2trK = 0;
-    PDupwindNthAnti2Xt2 = origin ? PDstandardNthAnti1Xt1 : Xt1[index] / xx;
-    PDupwindNthSymm2Xt2 = origin ? PDstandardNthSymm1Xt1 : Xt1[index] / xx;
+    PDupwindNthAnti2beta1 = kzero;
+    PDupwindNthSymm2beta1 = kzero;
+    PDupwindNthAnti2beta2 = kifthen(origin, PDstandardNthAnti1beta1, kmul(beta1L, xinv));
+    PDupwindNthSymm2beta2 = kifthen(origin, PDstandardNthSymm1beta1, kmul(beta1L, xinv));
+    PDupwindNthAnti2beta3 = kzero;
+    PDupwindNthSymm2beta3 = kzero;
 
+    PDupwindNthAnti2gt11 = kzero;
+    PDupwindNthSymm2gt11 = kzero;
+    PDupwindNthAnti2gt12 = kifthen(origin, ksub(PDstandardNthAnti1gt11, PDstandardNthAnti1gt22), kmul(ksub(gt11L, gt22L), xinv));
+    PDupwindNthSymm2gt12 = kifthen(origin, ksub(PDstandardNthSymm1gt11, PDstandardNthSymm1gt22), kmul(ksub(gt11L, gt22L), xinv));
+    PDupwindNthAnti2gt13 = kzero;
+    PDupwindNthSymm2gt13 = kzero;
+    PDupwindNthAnti2gt22 = kzero;
+    PDupwindNthSymm2gt22 = kzero;
+    PDupwindNthAnti2gt23 = kifthen(origin, PDstandardNthAnti1gt13, kmul(gt13L, xinv));
+    PDupwindNthSymm2gt23 = kifthen(origin, PDstandardNthSymm1gt13, kmul(gt13L, xinv));
+    PDupwindNthAnti2gt33 = kzero;
+    PDupwindNthSymm2gt33 = kzero;
+
+    PDupwindNthAnti2phi = kzero;
+    PDupwindNthSymm2phi = kzero;
+    PDupwindNthAnti2trK = kzero;
+    PDupwindNthSymm2trK = kzero;
+    PDupwindNthAnti2Xt2 = kifthen(origin, PDstandardNthAnti1Xt1, kmul(Xt1L, xinv));
+    PDupwindNthSymm2Xt2 = kifthen(origin, PDstandardNthSymm1Xt1, kmul(Xt1L, xinv));
diff --git a/cartoon_constr.c b/cartoon_constr.c
index badc89c..91b3445 100644
--- a/cartoon_constr.c
+++ b/cartoon_constr.c
@@ -1,43 +1,64 @@
-    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);
 
-    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 xinv  = kdiv(kone, xx);
+    CCTK_REAL_VEC x2inv = kdiv(kone, x2);
 
-    PDstandardNth2phi = 0;
-    PDstandardNth22phi = origin ? PDstandardNth11phi : PDstandardNth1phi / xx;
-    PDstandardNth12phi = 0;
-    PDstandardNth23phi = 0;
-    PDstandardNth2trK = 0;
-    PDstandardNth2Xt1 = 0;
-    PDstandardNth2Xt2 = origin ? PDstandardNth1Xt1 : Xt1[index] / xx;
-    PDstandardNth2Xt3 = 0;
+    CCTK_REAL_VEC absx = kfabs(xx);
+    CCTK_REAL_VEC eps  = ToReal(1e-8);
+    CCTK_BOOLEAN_VEC origin = kcmplt(absx, eps);
 
-    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;
+    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;
+
+    PDstandardNth2phi  = kzero;
+    PDstandardNth22phi = kifthen(origin, PDstandardNth11phi, kmul(PDstandardNth1phi, vinx));
+    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;