From 5209fc8f9afe5afc9a1315839bc7501aa103e34d Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Wed, 12 Aug 2020 13:17:20 +0200 Subject: Vectorize --- ML_BSSN/src/ML_BSSN_Advect.cc | 100 ++++++++-------- ML_BSSN/src/ML_BSSN_RHS.cc | 119 +++++++++++-------- ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.cc | 27 +++-- ML_BSSN/src/eval_timediff.cc | 119 +++++++++++-------- ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc | 94 ++++++++------- cartoon.c | 138 +++++++++++++---------- cartoon_advect.c | 97 ++++++++-------- cartoon_constr.c | 101 ++++++++++------- 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; -- cgit v1.2.3