1 files changed, 55 insertions, 45 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);