From 779de67ad5a8ffc90fd5b908941a492461473ebc Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Wed, 12 Aug 2020 17:29:21 +0200 Subject: Vectorize kretschmann --- ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc | 48 +++++++++++++----------- 1 file changed, 27 insertions(+), 21 deletions(-) diff --git a/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc b/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc index 6a8f0e6..89461ea 100644 --- a/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc +++ b/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc @@ -2084,34 +2084,40 @@ static void ML_Kretschmann_kretschmann_Body(const cGH* restrict const cctkGH, co CCTK_REAL_VEC KretschL CCTK_ATTRIBUTE_UNUSED = kadd(term1,kadd(term2,term3)); + CCTK_REAL_VEC zetaL; { - double xx = x[index]; - double x2 = SQR(xx); - double x4 = SQR(x2); - double reg_part = SQR(K22) / SQR(g22) - 0.25 * (SQR(Dg221) * gu11 + gu33 * SQR(Dg223) + 2.0 * Dg221 * Dg223 * gu13) / SQR(g22); - double lambda = (x2 * g22 * (1.0 - gu11 * g22) - xx * x2 * g22 * (Dg221 * gu11 + gu13 * Dg223)) / (x4 * SQR(g22)); + CCTK_REAL_VEC kquart = ToReal(0.25); + CCTK_REAL_VEC kfour = ToReal(4.0); + CCTK_REAL_VEC ksix = ToReal(6.0); - if (fabs(xx) < 1e-14) { - double dphi_x = PDstandardNth1phi; - double d2phi_xx = PDstandardNth11phi; - double dgt11_x = PDstandardNth1gt11; - double d2gt11_xx = PDstandardNth11gt11; - double Dgu131 = -Dg131 * (gu11 * gu33 + SQR(gu13)); - double Dgu111 = 0.0; - double D2g1111 = - 6.0 * SQR(dphi_x) * gt11L / SQR(SQR(phiL)) - - 2.0 * d2phi_xx * gt11L / (phiL * SQR(phiL)) - - 4.0 * dphi_x * dgt11_x / (phiL * SQR(phiL)) - + d2gt11_xx / SQR(phiL); - double D2gu1111 = -(SQR(gu11) * D2g1111 + 2.0 * Dg131 * (Dgu111 * gu13 + gu11 * Dgu131) + SQR(gu13) * D2g3311); - double d2x_guxx_gyy = D2gu1111 * g22 + gu11 * D2g2211; - lambda = -(gu11 * D2g2211 + Dgu131 * Dg223) / g22 - 0.5 * d2x_guxx_gyy / g22; - } - zeta[index] = reg_part + lambda; + CCTK_REAL_VEC x2 = SQR(xx); + CCTK_REAL_VEC x4 = SQR(x2); + + CCTK_REAL_VEC reg_part = SQR(K22) / SQR(g22) - kquart * (SQR(Dg221) * gu11 + gu33 * SQR(Dg223) + ktwo * Dg221 * Dg223 * gu13) / SQR(g22); + CCTK_REAL_VEC lambda_reg = (x2 * g22 * (kone - gu11 * g22) - xx * x2 * g22 * (Dg221 * gu11 + gu13 * Dg223)) / (x4 * SQR(g22)); + + CCTK_REAL_VEC dphi_x = PDstandardNth1phi; + CCTK_REAL_VEC d2phi_xx = PDstandardNth11phi; + CCTK_REAL_VEC dgt11_x = PDstandardNth1gt11; + CCTK_REAL_VEC d2gt11_xx = PDstandardNth11gt11; + CCTK_REAL_VEC Dgu131 = -Dg131 * (gu11 * gu33 + SQR(gu13)); + CCTK_REAL_VEC Dgu111 = kzero; + CCTK_REAL_VEC D2g1111 = - ksix * SQR(dphi_x) * gt11L / SQR(SQR(phiL)) + - ktwo * d2phi_xx * gt11L / (phiL * SQR(phiL)) + - kfour * dphi_x * dgt11_x / (phiL * SQR(phiL)) + + d2gt11_xx / SQR(phiL); + CCTK_REAL_VEC D2gu1111 = -(SQR(gu11) * D2g1111 + ktwo * Dg131 * (Dgu111 * gu13 + gu11 * Dgu131) + SQR(gu13) * D2g3311); + CCTK_REAL_VEC d2x_guxx_gyy = D2gu1111 * g22 + gu11 * D2g2211; + CCTK_REAL_VEC lambda_sing = -(gu11 * D2g2211 + Dgu131 * Dg223) / g22 - khalf * d2x_guxx_gyy / g22; + CCTK_REAL_VEC lambda = kifthen(origin, lambda_sing, lambda_reg); + + zetaL = reg_part + lambda; } /* Copy local copies back to grid functions */ vec_store_partial_prepare(i,vecimin,vecimax); vec_store_nta_partial(Kretsch[index],KretschL); + vec_store_nta_partial(zeta[index], zetaL); } CCTK_ENDLOOP3STR(ML_Kretschmann_kretschmann); } -- cgit v1.2.3