diff options
Diffstat (limited to 'ML_BSSN/src/ML_BSSN_RHS.cc')
| -rw-r--r-- | ML_BSSN/src/ML_BSSN_RHS.cc | 119 |
1 files changed, 69 insertions, 50 deletions
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); |