From 248a4737bc51b9739881564eea786cf350ae18e1 Mon Sep 17 00:00:00 2001 From: Anton Khirnov Date: Fri, 18 Sep 2020 12:57:45 +0200 Subject: Compute a function Kzeta = (K - 12 zeta)/rho^2 --- ML_Kretschmann/interface.ccl | 3 +- ML_Kretschmann/schedule.ccl | 1 + ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc | 111 +++++++++++++++++++++++ ML_Kretschmann/src/RegisterSymmetries.cc | 6 ++ 4 files changed, 120 insertions(+), 1 deletion(-) diff --git a/ML_Kretschmann/interface.ccl b/ML_Kretschmann/interface.ccl index c6e4182..eed9f3c 100644 --- a/ML_Kretschmann/interface.ccl +++ b/ML_Kretschmann/interface.ccl @@ -35,5 +35,6 @@ CCTK_REAL ML_Kretschmann type=GF timelevels=1 tags='tensortypealias="Scalar" ten CCTK_REAL ML_zeta type=GF timelevels=1 tags='tensortypealias="Scalar" tensorweight=0' { - zeta + zeta, + Kzeta } "ML_zeta" diff --git a/ML_Kretschmann/schedule.ccl b/ML_Kretschmann/schedule.ccl index 09788ab..416b188 100644 --- a/ML_Kretschmann/schedule.ccl +++ b/ML_Kretschmann/schedule.ccl @@ -39,6 +39,7 @@ schedule ML_Kretschmann_kretschmann in ML_Kretschmann_kretschmann_group READS: ML_BSSN::trK(Everywhere) WRITES: ML_Kretschmann::Kretsch(Interior) WRITES: ML_zeta::zeta(Interior) + WRITES: ML_zeta::Kzeta(Interior) } "ML_Kretschmann_kretschmann" schedule ML_Kretschmann_kretschmann_SelectBCs in ML_Kretschmann_kretschmann_bc_group diff --git a/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc b/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc index 89461ea..deb01dd 100644 --- a/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc +++ b/ML_Kretschmann/src/ML_Kretschmann_kretschmann.cc @@ -2122,6 +2122,116 @@ static void ML_Kretschmann_kretschmann_Body(const cGH* restrict const cctkGH, co CCTK_ENDLOOP3STR(ML_Kretschmann_kretschmann); } +static void ML_Kretschmann_kzeta(const cGH* restrict const cctkGH, const int dir, const int face, const CCTK_REAL normal[3], const CCTK_REAL tangentA[3], const CCTK_REAL tangentB[3], const int imin[3], const int imax[3], const int n_subblock_gfs, CCTK_REAL* restrict const subblock_gfs[]) +{ + DECLARE_CCTK_ARGUMENTS; + DECLARE_CCTK_PARAMETERS; + + + /* Include user-supplied include files */ + + /* Initialise finite differencing variables */ + const ptrdiff_t di CCTK_ATTRIBUTE_UNUSED = 1; + const ptrdiff_t dj CCTK_ATTRIBUTE_UNUSED = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0); + const ptrdiff_t dk CCTK_ATTRIBUTE_UNUSED = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0); + + + /* Loop over the grid points */ + const int imin0=imin[0]; + const int imin1=imin[1]; + const int imin2=imin[2]; + const int imax0=imax[0]; + const int imax1=imax[1]; + const int imax2=imax[2]; + #pragma omp parallel // reduction(+: vec_iter_counter, vec_op_counter, vec_mem_counter) + CCTK_LOOP3STR(ML_Kretschmann_kzeta, + i,j,k, imin0,imin1,imin2, imax0,imax1,imax2, + cctk_ash[0],cctk_ash[1],cctk_ash[2], + vecimin,vecimax, CCTK_REAL_VEC_SIZE) + { + const ptrdiff_t index CCTK_ATTRIBUTE_UNUSED = di*i + dj*j + dk*k; + // vec_iter_counter+=CCTK_REAL_VEC_SIZE; + + /* Assign local copies of grid functions */ + + CCTK_REAL_VEC gt22L CCTK_ATTRIBUTE_UNUSED = vec_load(gt22[index]); + CCTK_REAL_VEC gt22p1 CCTK_ATTRIBUTE_UNUSED = vec_load(gt22[index + 1]); + CCTK_REAL_VEC gt22p2 CCTK_ATTRIBUTE_UNUSED = vec_load(gt22[index + 2]); + CCTK_REAL_VEC phiL CCTK_ATTRIBUTE_UNUSED = vec_load(phi[index]); + CCTK_REAL_VEC phip1 CCTK_ATTRIBUTE_UNUSED = vec_load(phi[index + 1]); + CCTK_REAL_VEC phip2 CCTK_ATTRIBUTE_UNUSED = vec_load(phi[index + 2]); + + CCTK_REAL_VEC KretschL CCTK_ATTRIBUTE_UNUSED = vec_load(Kretsch[index]); + CCTK_REAL_VEC zetaL CCTK_ATTRIBUTE_UNUSED = vec_load(zeta[index]); + + CCTK_REAL_VEC Kretschp1 CCTK_ATTRIBUTE_UNUSED = vec_load(Kretsch[index + 1]); + CCTK_REAL_VEC Kretschp2 CCTK_ATTRIBUTE_UNUSED = vec_load(Kretsch[index + 2]); + CCTK_REAL_VEC zetap1 CCTK_ATTRIBUTE_UNUSED = vec_load(zeta[index + 1]); + CCTK_REAL_VEC zetap2 CCTK_ATTRIBUTE_UNUSED = vec_load(zeta[index + 2]); + + CCTK_REAL_VEC ktwo = ToReal(2.0); + + CCTK_REAL_VEC xx = vec_load(x[index]); + CCTK_REAL_VEC xp1 = vec_load(x[index + 1]); + CCTK_REAL_VEC xp2 = vec_load(x[index + 2]); + + CCTK_REAL_VEC absx = kfabs(xx); + CCTK_REAL_VEC eps = ToReal(1e-8); + CCTK_BOOLEAN_VEC origin = kcmplt(absx, eps); + + + CCTK_REAL_VEC e4phi CCTK_ATTRIBUTE_UNUSED; + CCTK_REAL_VEC e4phip1 CCTK_ATTRIBUTE_UNUSED; + CCTK_REAL_VEC e4phip2 CCTK_ATTRIBUTE_UNUSED; + + if (conformalMethod == 1) + { + e4phi = kdiv(ToReal(1),SQR(phiL)); + e4phip1 = kdiv(ToReal(1),SQR(phip1)); + e4phip2 = kdiv(ToReal(1),SQR(phip2)); + } + else + { + e4phi = kexp(kmul(phiL,ToReal(4))); + e4phip1 = kexp(kmul(phip1,ToReal(4))); + e4phip2 = kexp(kmul(phip2,ToReal(4))); + } + + CCTK_REAL_VEC g22 = kmul(gt22L, e4phi); + CCTK_REAL_VEC g22p1 = kmul(gt22p1,e4phip1); + CCTK_REAL_VEC g22p2 = kmul(gt22p2,e4phip2); + + CCTK_REAL_VEC KzetaL; + { + CCTK_REAL_VEC fact_interp_p1 = ToReal(0.5625); + CCTK_REAL_VEC fact_interp_p2 = ToReal(-0.0625); + + CCTK_REAL_VEC ksix = ToReal(6.0); + CCTK_REAL_VEC ktwelve = ToReal(12.0); + + CCTK_REAL_VEC x2 = SQR(xx); + CCTK_REAL_VEC x2p1 = SQR(xp1); + CCTK_REAL_VEC x2p2 = SQR(xp2); + + CCTK_REAL_VEC rho2 = g22 * x2; + CCTK_REAL_VEC rho2p1 = g22p1 * x2p1; + CCTK_REAL_VEC rho2p2 = g22p2 * x2p2; + + CCTK_REAL_VEC kzeta_reg = (KretschL - ktwelve * SQR(zetaL)) / rho2; + CCTK_REAL_VEC kzetap1 = (Kretschp1 - ktwelve * SQR(zetap1)) / rho2p1; + CCTK_REAL_VEC kzetap2 = (Kretschp2 - ktwelve * SQR(zetap2)) / rho2p2; + CCTK_REAL_VEC kzeta_sing = ktwo * (fact_interp_p1 * kzetap1 + fact_interp_p2 * kzetap2); + + KzetaL = kifthen(origin, kzeta_sing, kzeta_reg); + } + + /* Copy local copies back to grid functions */ + vec_store_partial_prepare(i,vecimin,vecimax); + vec_store_nta_partial(Kzeta[index], KzetaL); + } + CCTK_ENDLOOP3STR(ML_Kretschmann_kzeta); +} + extern "C" void ML_Kretschmann_kretschmann(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS; @@ -2139,6 +2249,7 @@ extern "C" void ML_Kretschmann_kretschmann(CCTK_ARGUMENTS) } GenericFD_LoopOverInterior(cctkGH, ML_Kretschmann_kretschmann_Body); + GenericFD_LoopOverInterior(cctkGH, ML_Kretschmann_kzeta); if (verbose > 1) { diff --git a/ML_Kretschmann/src/RegisterSymmetries.cc b/ML_Kretschmann/src/RegisterSymmetries.cc index 9e3ce3f..ad4d1ce 100644 --- a/ML_Kretschmann/src/RegisterSymmetries.cc +++ b/ML_Kretschmann/src/RegisterSymmetries.cc @@ -26,4 +26,10 @@ extern "C" void ML_Kretschmann_RegisterSymmetries(CCTK_ARGUMENTS) sym[1] = 1; sym[2] = 1; SetCartSymVN(cctkGH, sym, "ML_Kretschmann::zeta"); + + /* Register symmetries of grid functions */ + sym[0] = 1; + sym[1] = 1; + sym[2] = 1; + SetCartSymVN(cctkGH, sym, "ML_Kretschmann::Kzeta"); } -- cgit v1.2.3