diff --git a/Weyl/configuration.ccl b/Weyl/configuration.ccl index 0f28756c..64cc7a6a 100644 --- a/Weyl/configuration.ccl +++ b/Weyl/configuration.ccl @@ -1,3 +1,3 @@ # Configuration definitions for thorn Weyl -REQUIRES Arith Loop +REQUIRES Arith Loop Derivs diff --git a/Weyl/interface.ccl b/Weyl/interface.ccl index 76f01e3a..9b5e575a 100644 --- a/Weyl/interface.ccl +++ b/Weyl/interface.ccl @@ -15,6 +15,10 @@ USES INCLUDE HEADER: ten3.hxx USES INCLUDE HEADER: vec.hxx USES INCLUDE HEADER: vect.hxx +INHERITS: Derivs + +USES INCLUDE HEADER: derivs_spacetimex.hxx + # TODO: Declare these variables without ghost zones? diff --git a/Weyl/src/derivs.hxx b/Weyl/src/derivs.hxx deleted file mode 100644 index eaa58efd..00000000 --- a/Weyl/src/derivs.hxx +++ /dev/null @@ -1,358 +0,0 @@ -#ifndef DERIVS_HXX -#define DERIVS_HXX - -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -namespace Weyl { -using namespace Arith; -using namespace Loop; -using namespace std; - -//////////////////////////////////////////////////////////////////////////////// - -constexpr int deriv_order = 4; - -//////////////////////////////////////////////////////////////////////////////// - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST T pow3(const T &x) { - return pow2(x) * x; -} - -//////////////////////////////////////////////////////////////////////////////// - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv1d(const simdl &mask, const T *restrict const var, const ptrdiff_t di, - const T dx) { - const auto load = [&](const int n) { - return maskz_loadu(mask, &var[n * di]) - maskz_loadu(mask, &var[-n * di]); - }; - if constexpr (deriv_order == 2) - return 1 / T(2) * load(1) / dx; - if constexpr (deriv_order == 4) - return (-1 / T(12) * load(2) + 2 / T(3) * load(1)) / dx; - if constexpr (deriv_order == 6) - return (1 / T(60) * load(3) - 3 / T(20) * load(2) + 3 / T(4) * load(1)) / - dx; -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv2_1d(const simdl &mask, const T *restrict const var, const ptrdiff_t di, - const T dx) { - const auto load = [&](const int n) { - return maskz_loadu(mask, &var[n * di]) + maskz_loadu(mask, &var[-n * di]); - }; - const auto load0 = [&]() { return maskz_loadu(mask, &var[0]); }; - if constexpr (deriv_order == 2) - return (load(1) - 2 * load0()) / pow2(dx); - if constexpr (deriv_order == 4) - return (1 / T(12) * load(2) - 4 / T(3) * load(1) + 5 / T(2) * load0()) / - pow2(dx); - if constexpr (deriv_order == 6) - return (1 / T(90) * load(3) - 3 / T(20) * load(2) + 3 / T(2) * load(1) - - 49 / T(18) * load0()) / - pow2(dx); -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv3_1d(const simdl &mask, const T *restrict const var, const ptrdiff_t di, - const T dx) { - const auto load = [&](const int n) { - return maskz_loadu(mask, &var[n * di]) - maskz_loadu(mask, &var[-n * di]); - }; - if constexpr (deriv_order == 2) - return (1 / T(2) * load(2) - load(1)) / pow3(dx); - if constexpr (deriv_order == 4) - return (-1 / T(8) * load(3) + load(2) - 13 / T(8) * load(1)) / pow3(dx); - if constexpr (deriv_order == 6) - return (7 / T(240) * load(4) - 3 / T(10) * load(3) + - 169 / T(129) * load(2) - 61 / T(30) * load(1)) / - pow3(dx); -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv2_2d(const int vavail, const simdl &mask, const T *restrict const var, - const ptrdiff_t di, const ptrdiff_t dj, const T dx, const T dy) { - constexpr size_t vsize = tuple_size_v >; - if (di == 1) { - assert(vavail > 0); - constexpr int maxnpoints = deriv_order + 1 + vsize - 1; - const int npoints = deriv_order + 1 + min(int(vsize), vavail) - 1; - array, div_ceil(maxnpoints, int(vsize))> arrx; - for (int i = 0; i < maxnpoints; i += vsize) { - if (i < npoints) { - const simdl mask1 = mask_for_loop_tail >(i, npoints); - arrx[div_floor(i, int(vsize))] = - deriv1d(mask1, &var[i - deriv_order / 2], dj, dy); - } - } -#ifdef CCTK_DEBUG - for (int i = npoints; i < align_ceil(maxnpoints, int(vsize)); ++i) - ((T *)&arrx[0])[i] = Arith::nan()(); // unused -#endif - const T *const varx = (T *)&arrx[0] + deriv_order / 2; - return deriv1d(mask, varx, 1, dx); - } else { - assert(dj != 1); - array, deriv_order + 1> arrx; - for (int j = -deriv_order / 2; j <= deriv_order / 2; ++j) - if (j == 0) { -#ifdef CCTK_DEBUG - arrx[deriv_order / 2 + j] = Arith::nan >()(); // unused -#endif - } else { - arrx[deriv_order / 2 + j] = deriv1d(mask, &var[j * dj], di, dx); - } - const T *const varx = (T *)(&arrx[deriv_order / 2]); - return deriv1d(mask, varx, vsize, dy); - } -} - -//////////////////////////////////////////////////////////////////////////////// - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv(const simdl &mask, const GF3D2 &gf_, const vect &I, - const vec &dx) { - static_assert(dir >= 0 && dir < D, ""); - const auto &DI = vect::unit; - const ptrdiff_t di = gf_.delta(DI(dir)); - return deriv1d(mask, &gf_(I), di, dx(dir)); -} - -template -inline ARITH_INLINE - ARITH_DEVICE ARITH_HOST enable_if_t<(dir1 == dir2), simd > - deriv2(const int vavail, const simdl &mask, const GF3D2 &gf_, - const vect &I, const vec &dx) { - static_assert(dir1 >= 0 && dir1 < D, ""); - static_assert(dir2 >= 0 && dir2 < D, ""); - const auto &DI = vect::unit; - const ptrdiff_t di = gf_.delta(DI(dir1)); - return deriv2_1d(mask, &gf_(I), di, dx(dir1)); -} - -template -inline ARITH_INLINE - ARITH_DEVICE ARITH_HOST enable_if_t<(dir1 != dir2), simd > - deriv2(const int vavail, const simdl &mask, const GF3D2 &gf_, - const vect &I, const vec &dx) { - static_assert(dir1 >= 0 && dir1 < D, ""); - static_assert(dir2 >= 0 && dir2 < D, ""); - const auto &DI = vect::unit; - const ptrdiff_t di = gf_.delta(DI(dir1)); - const ptrdiff_t dj = gf_.delta(DI(dir2)); - return deriv2_2d(vavail, mask, &gf_(I), di, dj, dx(dir1), dx(dir2)); -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST - enable_if_t<(dir1 == dir2 && dir1 == dir3), simd > - deriv3(const int vavail, const simdl &mask, const GF3D2 &gf_, - const vect &I, const vec &dx) { - static_assert(dir1 >= 0 && dir1 < D, ""); - static_assert(dir2 >= 0 && dir2 < D, ""); - static_assert(dir3 >= 0 && dir3 < D, ""); - const auto &DI = vect::unit; - const ptrdiff_t di = gf_.delta(DI(dir1)); - return deriv3_1d(mask, &gf_(I), di, dx(dir1)); -} - -//////////////////////////////////////////////////////////////////////////////// - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST vec, dim> -deriv(const simdl &mask, const GF3D2 &gf_, const vect &I, - const vec &dx) { - return {deriv<0>(mask, gf_, I, dx), deriv<1>(mask, gf_, I, dx), - deriv<2>(mask, gf_, I, dx)}; -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST smat, dim> -deriv2(const int vavail, const simdl &mask, const GF3D2 &gf_, - const vect &I, const vec &dx) { - return {deriv2<0, 0>(vavail, mask, gf_, I, dx), - deriv2<0, 1>(vavail, mask, gf_, I, dx), - deriv2<0, 2>(vavail, mask, gf_, I, dx), - deriv2<1, 1>(vavail, mask, gf_, I, dx), - deriv2<1, 2>(vavail, mask, gf_, I, dx), - deriv2<2, 2>(vavail, mask, gf_, I, dx)}; -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST sten3, dim> -deriv3(const int vavail, const simdl &mask, const GF3D2 &gf_, - const vect &I, const vec &dx) { - return {deriv3<0, 0, 0>(vavail, mask, gf_, I, dx), - deriv3<0, 0, 1>(vavail, mask, gf_, I, dx), - deriv3<0, 0, 2>(vavail, mask, gf_, I, dx), - deriv3<0, 1, 1>(vavail, mask, gf_, I, dx), - deriv3<0, 1, 2>(vavail, mask, gf_, I, dx), - deriv3<0, 2, 2>(vavail, mask, gf_, I, dx), - deriv3<1, 1, 1>(vavail, mask, gf_, I, dx), - deriv3<1, 1, 2>(vavail, mask, gf_, I, dx), - deriv3<1, 2, 2>(vavail, mask, gf_, I, dx), - deriv3<2, 2, 2>(vavail, mask, gf_, I, dx)}; -} - -//////////////////////////////////////////////////////////////////////////////// - -template -CCTK_ATTRIBUTE_NOINLINE void -calc_copy(const cGH *restrict const cctkGH, const GF3D2 &gf1, - const GF3D5 &gf0, const GF3D5layout &layout0) { - DECLARE_CCTK_ARGUMENTS; - - typedef simd vreal; - typedef simdl vbool; - constexpr size_t vsize = tuple_size_v; - - const vec dx([&](int a) { return CCTK_DELTA_SPACE(a); }); - - const Loop::GridDescBaseDevice grid(cctkGH); - grid.loop_int_device<0, 0, 0, vsize>( - grid.nghostzones, - [=] ARITH_DEVICE ARITH_HOST(const PointDesc &p) ARITH_INLINE { - const vbool mask = mask_for_loop_tail(p.i, p.imax); - const GF3D5index index0(layout0, p.I); - const auto val = gf1(mask, p.I); - gf0.store(mask, index0, val); - }); -} - -template -CCTK_ATTRIBUTE_NOINLINE void -calc_derivs(const cGH *restrict const cctkGH, const GF3D2 &gf1, - const GF3D5 &gf0, const vec, dim> &dgf0, - const GF3D5layout &layout0) { - DECLARE_CCTK_ARGUMENTS; - - typedef simd vreal; - typedef simdl vbool; - constexpr size_t vsize = tuple_size_v; - - const vec dx([&](int a) { return CCTK_DELTA_SPACE(a); }); - - const Loop::GridDescBaseDevice grid(cctkGH); - grid.loop_int_device<0, 0, 0, vsize>( - grid.nghostzones, - [=] ARITH_DEVICE ARITH_HOST(const PointDesc &p) ARITH_INLINE { - const vbool mask = mask_for_loop_tail(p.i, p.imax); - const GF3D5index index0(layout0, p.I); - const auto val = gf1(mask, p.I); - gf0.store(mask, index0, val); - const auto dval = deriv(mask, gf1, p.I, dx); - dgf0.store(mask, index0, dval); - }); -} - -template -CCTK_ATTRIBUTE_NOINLINE void -calc_derivs2(const cGH *restrict const cctkGH, const GF3D2 &gf1, - const GF3D5 &gf0, const vec, dim> &dgf0, - const smat, dim> &ddgf0, const GF3D5layout &layout0) { - DECLARE_CCTK_ARGUMENTS; - - typedef simd vreal; - typedef simdl vbool; - constexpr size_t vsize = tuple_size_v; - - const vec dx([&](int a) { return CCTK_DELTA_SPACE(a); }); - - const Loop::GridDescBaseDevice grid(cctkGH); - grid.loop_int_device<0, 0, 0, vsize>( - grid.nghostzones, - [=] ARITH_DEVICE ARITH_HOST(const PointDesc &p) ARITH_INLINE { - const vbool mask = mask_for_loop_tail(p.i, p.imax); - const GF3D5index index0(layout0, p.I); - const auto val = gf1(mask, p.I); - gf0.store(mask, index0, val); - const auto dval = deriv(mask, gf1, p.I, dx); - dgf0.store(mask, index0, dval); - const auto ddval = deriv2(p.imax - p.i, mask, gf1, p.I, dx); - ddgf0.store(mask, index0, ddval); - }); -} - -template -void CCTK_ATTRIBUTE_NOINLINE calc_copy(const cGH *restrict const cctkGH, - const vec, dim> &gf0_, - const vec, dim> &gf_, - const GF3D5layout &layout) { - for (int a = 0; a < 3; ++a) - calc_copy(cctkGH, gf0_(a), gf_(a), layout); -} - -template -void CCTK_ATTRIBUTE_NOINLINE calc_derivs( - const cGH *restrict const cctkGH, const vec, dim> &gf0_, - const vec, dim> &gf_, const vec, dim>, dim> &dgf_, - const GF3D5layout &layout) { - for (int a = 0; a < 3; ++a) - calc_derivs(cctkGH, gf0_(a), gf_(a), dgf_(a), layout); -} - -template -void CCTK_ATTRIBUTE_NOINLINE calc_derivs2( - const cGH *restrict const cctkGH, const vec, dim> &gf0_, - const vec, dim> &gf_, const vec, dim>, dim> &dgf_, - const vec, dim>, dim> &ddgf_, const GF3D5layout &layout) { - for (int a = 0; a < 3; ++a) - calc_derivs2(cctkGH, gf0_(a), gf_(a), dgf_(a), ddgf_(a), layout); -} - -template -void CCTK_ATTRIBUTE_NOINLINE calc_copy(const cGH *restrict const cctkGH, - const smat, dim> &gf0_, - const smat, dim> &gf_, - const GF3D5layout &layout) { - for (int a = 0; a < 3; ++a) - for (int b = a; b < 3; ++b) - calc_copy(cctkGH, gf0_(a, b), gf_(a, b), layout); -} - -template -void CCTK_ATTRIBUTE_NOINLINE calc_derivs( - const cGH *restrict const cctkGH, const smat, dim> &gf0_, - const smat, dim> &gf_, const smat, dim>, dim> &dgf_, - const GF3D5layout &layout) { - for (int a = 0; a < 3; ++a) - for (int b = a; b < 3; ++b) - calc_derivs(cctkGH, gf0_(a, b), gf_(a, b), dgf_(a, b), layout); -} - -template -void CCTK_ATTRIBUTE_NOINLINE calc_derivs2( - const cGH *restrict const cctkGH, const smat, dim> &gf0_, - const smat, dim> &gf_, const smat, dim>, dim> &dgf_, - const smat, dim>, dim> &ddgf_, const GF3D5layout &layout) { - for (int a = 0; a < 3; ++a) - for (int b = a; b < 3; ++b) - calc_derivs2(cctkGH, gf0_(a, b), gf_(a, b), dgf_(a, b), ddgf_(a, b), - layout); -} - -} // namespace Weyl - -#endif // #ifndef DERIVS_HXX diff --git a/Weyl/src/weyl.cxx b/Weyl/src/weyl.cxx index deba4f59..cb418e90 100644 --- a/Weyl/src/weyl.cxx +++ b/Weyl/src/weyl.cxx @@ -10,7 +10,7 @@ #endif #endif -#include "derivs.hxx" +#include "derivs_spacetimex.hxx" #include "physics.hxx" #include "weyl_vars.hxx" @@ -27,6 +27,8 @@ #include +using namespace derivs_spacetimex; + namespace Weyl { using namespace Arith; using namespace Loop; diff --git a/Weyl/src/weyl_vars.hxx b/Weyl/src/weyl_vars.hxx index 2c58a4d3..c2e8d8f8 100644 --- a/Weyl/src/weyl_vars.hxx +++ b/Weyl/src/weyl_vars.hxx @@ -1,7 +1,7 @@ #ifndef WEYL_VARS_HXX #define WEYL_VARS_HXX -#include "derivs.hxx" +#include "derivs_spacetimex.hxx" #include "physics.hxx" #include diff --git a/Z4c/interface.ccl b/Z4c/interface.ccl index 350d5408..9619c86c 100644 --- a/Z4c/interface.ccl +++ b/Z4c/interface.ccl @@ -14,6 +14,10 @@ USES INCLUDE HEADER: sum.hxx USES INCLUDE HEADER: vec.hxx USES INCLUDE HEADER: vect.hxx +INHERITS: Derivs + +USES INCLUDE HEADER: derivs_spacetimex.hxx + CCTK_INT FUNCTION GetCallFunctionCount() diff --git a/Z4c/param.ccl b/Z4c/param.ccl index b3d3ba2c..97ba98fc 100644 --- a/Z4c/param.ccl +++ b/Z4c/param.ccl @@ -52,8 +52,3 @@ CCTK_REAL alphaG_floor "Floor for alphaG" STEERABLE=always { (0:* :: "" } 1.0e-10 - -CCTK_REAL epsdiss "Dissipation coefficient " STEERABLE=always -{ - 0.0:* :: "" -} 0.32 diff --git a/Z4c/src/adm2.cxx b/Z4c/src/adm2.cxx index 3fbdc94c..20427b63 100644 --- a/Z4c/src/adm2.cxx +++ b/Z4c/src/adm2.cxx @@ -8,7 +8,7 @@ #endif #endif -#include "derivs.hxx" +#include "derivs_spacetimex.hxx" #include "physics.hxx" #include "z4c_vars.hxx" @@ -27,6 +27,8 @@ #include +using namespace derivs_spacetimex; + namespace Z4c { using namespace Arith; using namespace Loop; diff --git a/Z4c/src/constraints.cxx b/Z4c/src/constraints.cxx index faf32e3a..9425d7df 100644 --- a/Z4c/src/constraints.cxx +++ b/Z4c/src/constraints.cxx @@ -8,7 +8,7 @@ #endif #endif -#include "derivs.hxx" +#include "derivs_spacetimex.hxx" #include "physics.hxx" #include "z4c_vars.hxx" @@ -25,6 +25,8 @@ #include +using namespace derivs_spacetimex; + namespace Z4c { using namespace Arith; using namespace Loop; diff --git a/Z4c/src/derivs.hxx b/Z4c/src/derivs.hxx deleted file mode 100644 index 6fbe6f03..00000000 --- a/Z4c/src/derivs.hxx +++ /dev/null @@ -1,455 +0,0 @@ -#ifndef DERIVS_HXX -#define DERIVS_HXX - -#include -#include -#include -#include -#include - -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -namespace Z4c { -using namespace Arith; -using namespace Loop; -using namespace std; - -//////////////////////////////////////////////////////////////////////////////// - -constexpr int deriv_order = 4; - -//////////////////////////////////////////////////////////////////////////////// - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv1d(const simdl &mask, const T *restrict const var, const ptrdiff_t di, - const T dx) { - if constexpr (deriv_order == 2) - return -1 / T(2) * - (maskz_loadu(mask, &var[-di]) - maskz_loadu(mask, &var[+di])) / dx; - if constexpr (deriv_order == 4) - return (1 / T(12) * - (maskz_loadu(mask, &var[-2 * di]) - - maskz_loadu(mask, &var[+2 * di])) // - - - 2 / T(3) * - (maskz_loadu(mask, &var[-di]) - maskz_loadu(mask, &var[+di]))) / - dx; -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv1d_upwind(const simdl &mask, const T *restrict const var, - const ptrdiff_t di, const simd &vel, const T dx) { - // arXiv:1111.2177 [gr-qc], (71) - if constexpr (deriv_order == 2) { - // if (sign) - // // + [ 0 -1 +1 0 0] - // // + 1/2 [+1 -2 +1 0 0] - // // [+1/2 -2 +3/2 0 0] - // return (1 / T(2) * var[-2 * di] // - // - 2 * var[-di] // - // + 3 / T(2) * var[0]) / - // dx; - // else - // // + [ 0 0 -1 +1 0 ] - // // - 1/2 [ 0 0 +1 -2 +1 ] - // // [ 0 0 -3/2 +2 -1/2] - // return (-3 / T(2) * var[0] // - // + 2 * var[+di] // - // - 1 / T(2) * var[+2 * di]) / - // dx; - - // + 1/2 [+1/2 -2 +3/2 0 0 ] - // + 1/2 [ 0 0 -3/2 +2 -1/2] - // [+1/4 -1 0 +1 -1/4] - const simd symm = - 1 / T(4) * - (maskz_loadu(mask, &var[-2 * di]) - - maskz_loadu(mask, &var[2 * di])) // - - (maskz_loadu(mask, &var[-di]) - maskz_loadu(mask, &var[di])); - // + 1/2 [+1/2 -2 +3/2 0 0 ] - // - 1/2 [ 0 0 -3/2 +2 -1/2] - // [+1/4 -1 +3/2 -1 +1/4] - const simd anti = - 1 / T(4) * - (maskz_loadu(mask, &var[-2 * di]) + - maskz_loadu(mask, &var[2 * di])) // - - (maskz_loadu(mask, &var[-di]) + maskz_loadu(mask, &var[di])) // - + 3 / T(2) * maskz_loadu(mask, &var[0]); - return (vel * symm - fabs(vel) * anti) / dx; - } - if constexpr (deriv_order == 4) { - // A fourth order stencil for a first derivative, shifted by one grid point - - // if (sign) - // return (-1 / T(12) * var[-3 * di] // - // + 1 / T(2) * var[-2 * di] // - // - 3 / T(2) * var[-di] // - // + 5 / T(6) * var[0] // - // + 1 / T(4) * var[+di]) / - // dx; - // else - // return (-1 / T(4) * var[-di] // - // - 5 / T(6) * var[0] // - // + 3 / T(2) * var[+di] // - // - 1 / T(2) * var[+2 * di] // - // + 1 / T(12) * var[+3 * di]) / - // dx; - - const simd symm = - -1 / T(24) * - (maskz_loadu(mask, &var[-3 * di]) - - maskz_loadu(mask, &var[3 * di])) // - + 1 / T(4) * - (maskz_loadu(mask, &var[-2 * di]) - - maskz_loadu(mask, &var[2 * di])) // - - - 7 / T(8) * (maskz_loadu(mask, &var[-di]) - maskz_loadu(mask, &var[di])); - const simd anti = - -1 / T(24) * - (maskz_loadu(mask, &var[-3 * di]) + - maskz_loadu(mask, &var[3 * di])) // - + 1 / T(4) * - (maskz_loadu(mask, &var[-2 * di]) + - maskz_loadu(mask, &var[2 * di])) // - - 5 / T(8) * - (maskz_loadu(mask, &var[-di]) + maskz_loadu(mask, &var[di])) // - + 5 / T(6) * maskz_loadu(mask, &var[0]); - return (vel * symm - fabs(vel) * anti) / dx; - } -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv2_1d(const simdl &mask, const T *restrict const var, const ptrdiff_t di, - const T dx) { - if constexpr (deriv_order == 2) - return ((maskz_loadu(mask, &var[-di]) + maskz_loadu(mask, &var[+di])) // - - 2 * maskz_loadu(mask, &var[0])) / - pow2(dx); - if constexpr (deriv_order == 4) - return (-1 / T(12) * - (maskz_loadu(mask, &var[-2 * di]) + - maskz_loadu(mask, &var[+2 * di])) // - + - 4 / T(3) * - (maskz_loadu(mask, &var[-di]) + maskz_loadu(mask, &var[+di])) // - - 5 / T(2) * maskz_loadu(mask, &var[0])) / - pow2(dx); -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv2_2d(const int vavail, const simdl &mask, const T *restrict const var, - const ptrdiff_t di, const ptrdiff_t dj, const T dx, const T dy) { - constexpr size_t vsize = tuple_size_v >; - if (di == 1) { - assert(vavail > 0); - constexpr int maxnpoints = deriv_order + 1 + vsize - 1; - const int npoints = deriv_order + 1 + min(int(vsize), vavail) - 1; - array, div_ceil(maxnpoints, int(vsize))> arrx; - for (int i = 0; i < maxnpoints; i += vsize) { - if (i < npoints) { - const simdl mask1 = mask_for_loop_tail >(i, npoints); - arrx[div_floor(i, int(vsize))] = - deriv1d(mask1, &var[i - deriv_order / 2], dj, dy); - } - } -#ifdef CCTK_DEBUG - for (int i = npoints; i < align_ceil(maxnpoints, int(vsize)); ++i) - ((T *)&arrx[0])[i] = Arith::nan()(); // unused -#endif - const T *const varx = (T *)&arrx[0] + deriv_order / 2; - return deriv1d(mask, varx, 1, dx); - } else { - assert(dj != 1); - array, deriv_order + 1> arrx; - for (int j = -deriv_order / 2; j <= deriv_order / 2; ++j) - if (j == 0) { -#ifdef CCTK_DEBUG - arrx[deriv_order / 2 + j] = Arith::nan >()(); // unused -#endif - } else { - arrx[deriv_order / 2 + j] = deriv1d(mask, &var[j * dj], di, dx); - } - const T *const varx = (T *)(&arrx[deriv_order / 2]); - return deriv1d(mask, varx, vsize, dy); - } -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv1d_diss(const simdl &mask, const T *restrict const var, - const ptrdiff_t di, const T dx) { - if constexpr (deriv_order == 2) - return ((maskz_loadu(mask, &var[-2 * di]) + - maskz_loadu(mask, &var[+2 * di])) // - - - 4 * (maskz_loadu(mask, &var[-di]) + maskz_loadu(mask, &var[+di])) // - + 6 * maskz_loadu(mask, &var[0])) / - dx; - if constexpr (deriv_order == 4) - return ((maskz_loadu(mask, &var[-3 * di]) + - maskz_loadu(mask, &var[+3 * di])) // - - 6 * (maskz_loadu(mask, &var[-2 * di]) + - maskz_loadu(mask, &var[+2 * di])) // - + 15 * (maskz_loadu(mask, &var[-di]) + - maskz_loadu(mask, &var[+di])) // - - 20 * maskz_loadu(mask, &var[0])) / - dx; -} - -//////////////////////////////////////////////////////////////////////////////// - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv(const simdl &mask, const GF3D2 &gf_, const vect &I, - const vec &dx) { - static_assert(dir >= 0 && dir < D, ""); - const auto &DI = vect::unit; - const ptrdiff_t di = gf_.delta(DI(dir)); - return deriv1d(mask, &gf_(I), di, dx(dir)); -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv_upwind(const simdl &mask, const GF3D2 &gf_, - const vect &I, const vec, D> &vel, - const vec &dx) { - static_assert(dir >= 0 && dir < D, ""); - const auto &DI = vect::unit; - const ptrdiff_t di = gf_.delta(DI(dir)); - return deriv1d_upwind(mask, &gf_(I), di, vel(dir), dx(dir)); -} - -template -inline ARITH_INLINE - ARITH_DEVICE ARITH_HOST enable_if_t<(dir1 == dir2), simd > - deriv2(const int vavail, const simdl &mask, const GF3D2 &gf_, - const vect &I, const vec &dx) { - static_assert(dir1 >= 0 && dir1 < D, ""); - static_assert(dir2 >= 0 && dir2 < D, ""); - const auto &DI = vect::unit; - const ptrdiff_t di = gf_.delta(DI(dir1)); - return deriv2_1d(mask, &gf_(I), di, dx(dir1)); -} - -template -inline ARITH_INLINE - ARITH_DEVICE ARITH_HOST enable_if_t<(dir1 != dir2), simd > - deriv2(const int vavail, const simdl &mask, const GF3D2 &gf_, - const vect &I, const vec &dx) { - static_assert(dir1 >= 0 && dir1 < D, ""); - static_assert(dir2 >= 0 && dir2 < D, ""); - const auto &DI = vect::unit; - const ptrdiff_t di = gf_.delta(DI(dir1)); - const ptrdiff_t dj = gf_.delta(DI(dir2)); - return deriv2_2d(vavail, mask, &gf_(I), di, dj, dx(dir1), dx(dir2)); -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv_diss(const simdl &mask, const GF3D2 &gf_, - const vect &I, const vec &dx) { - static_assert(dir >= 0 && dir < D, ""); - const auto &DI = vect::unit; - const ptrdiff_t di = gf_.delta(DI(dir)); - return deriv1d_diss(mask, &gf_(I), di, dx(dir)); -} - -//////////////////////////////////////////////////////////////////////////////// - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST vec, dim> -deriv(const simdl &mask, const GF3D2 &gf_, const vect &I, - const vec &dx) { - return {deriv<0>(mask, gf_, I, dx), deriv<1>(mask, gf_, I, dx), - deriv<2>(mask, gf_, I, dx)}; -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv_upwind(const simdl &mask, const GF3D2 &gf_, - const vect &I, const vec, dim> &vel, - const vec &dx) { - return deriv_upwind<0>(mask, gf_, I, vel, dx) + - deriv_upwind<1>(mask, gf_, I, vel, dx) + - deriv_upwind<2>(mask, gf_, I, vel, dx); -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST smat, dim> -deriv2(const int vavail, const simdl &mask, const GF3D2 &gf_, - const vect &I, const vec &dx) { - return {deriv2<0, 0>(vavail, mask, gf_, I, dx), - deriv2<0, 1>(vavail, mask, gf_, I, dx), - deriv2<0, 2>(vavail, mask, gf_, I, dx), - deriv2<1, 1>(vavail, mask, gf_, I, dx), - deriv2<1, 2>(vavail, mask, gf_, I, dx), - deriv2<2, 2>(vavail, mask, gf_, I, dx)}; -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -diss(const simdl &mask, const GF3D2 &gf_, const vect &I, - const vec &dx) { - // arXiv:gr-qc/0610128, (63), with r=2 - constexpr int diss_order = deriv_order + 2; - constexpr int sign = diss_order % 4 == 0 ? -1 : +1; - return sign / T(pown(2, deriv_order + 2)) * - (deriv_diss<0>(mask, gf_, I, dx) // - + deriv_diss<1>(mask, gf_, I, dx) // - + deriv_diss<2>(mask, gf_, I, dx)); -} - -//////////////////////////////////////////////////////////////////////////////// - -template -CCTK_ATTRIBUTE_NOINLINE void -calc_derivs(const cGH *restrict const cctkGH, const GF3D2 &gf1, - const GF3D5 &gf0, const vec, dim> &dgf0, - const GF3D5layout &layout0) { - DECLARE_CCTK_ARGUMENTS; - - typedef simd vreal; - typedef simdl vbool; - constexpr size_t vsize = tuple_size_v; - - const vec dx([&](int a) { return CCTK_DELTA_SPACE(a); }); - - const Loop::GridDescBaseDevice grid(cctkGH); - grid.loop_int_device<0, 0, 0, vsize>( - grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE { - const vbool mask = mask_for_loop_tail(p.i, p.imax); - const GF3D5index index0(layout0, p.I); - const auto val = gf1(mask, p.I); - gf0.store(mask, index0, val); - const auto dval = deriv(mask, gf1, p.I, dx); - dgf0.store(mask, index0, dval); - }); -} - -template -CCTK_ATTRIBUTE_NOINLINE void -calc_derivs2(const cGH *restrict const cctkGH, const GF3D2 &gf1, - const GF3D5 &gf0, const vec, dim> &dgf0, - const smat, dim> &ddgf0, const GF3D5layout &layout0) { - DECLARE_CCTK_ARGUMENTS; - - typedef simd vreal; - typedef simdl vbool; - constexpr size_t vsize = tuple_size_v; - - const vec dx([&](int a) { return CCTK_DELTA_SPACE(a); }); - - const Loop::GridDescBaseDevice grid(cctkGH); - grid.loop_int_device<0, 0, 0, vsize>( - grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE { - const vbool mask = mask_for_loop_tail(p.i, p.imax); - const int vavail = p.imax - p.i; - const GF3D5index index0(layout0, p.I); - const auto val = gf1(mask, p.I); - gf0.store(mask, index0, val); - const auto dval = deriv(mask, gf1, p.I, dx); - dgf0.store(mask, index0, dval); - const auto ddval = deriv2(vavail, mask, gf1, p.I, dx); - ddgf0.store(mask, index0, ddval); - }); -} - -template -CCTK_ATTRIBUTE_NOINLINE void -calc_derivs(const cGH *restrict const cctkGH, - const vec, dim> &gf0_, const vec, dim> &gf_, - const vec, dim>, dim> &dgf_, - const GF3D5layout &layout) { - for (int a = 0; a < 3; ++a) - calc_derivs(cctkGH, gf0_(a), gf_(a), dgf_(a), layout); -} - -template -CCTK_ATTRIBUTE_NOINLINE void calc_derivs2( - const cGH *restrict const cctkGH, const vec, dim> &gf0_, - const vec, dim> &gf_, const vec, dim>, dim> &dgf_, - const vec, dim>, dim> &ddgf_, const GF3D5layout &layout) { - for (int a = 0; a < 3; ++a) - calc_derivs2(cctkGH, gf0_(a), gf_(a), dgf_(a), ddgf_(a), layout); -} - -template -CCTK_ATTRIBUTE_NOINLINE void calc_derivs( - const cGH *restrict const cctkGH, const smat, dim> &gf0_, - const smat, dim> &gf_, const smat, dim>, dim> &dgf_, - const GF3D5layout &layout) { - for (int a = 0; a < 3; ++a) - for (int b = a; b < 3; ++b) - calc_derivs(cctkGH, gf0_(a, b), gf_(a, b), dgf_(a, b), layout); -} - -template -CCTK_ATTRIBUTE_NOINLINE void calc_derivs2( - const cGH *restrict const cctkGH, const smat, dim> &gf0_, - const smat, dim> &gf_, const smat, dim>, dim> &dgf_, - const smat, dim>, dim> &ddgf_, const GF3D5layout &layout) { - for (int a = 0; a < 3; ++a) - for (int b = a; b < 3; ++b) - calc_derivs2(cctkGH, gf0_(a, b), gf_(a, b), dgf_(a, b), ddgf_(a, b), - layout); -} - -template -CCTK_ATTRIBUTE_NOINLINE void -apply_upwind_diss(const cGH *restrict const cctkGH, const GF3D2 &gf_, - const vec, dim> &gf_betaG_, - const GF3D2 &gf_rhs_) { - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - const vec dx([&](int a) { return CCTK_DELTA_SPACE(a); }); - - typedef simd vreal; - typedef simdl vbool; - constexpr size_t vsize = tuple_size_v; - - if (epsdiss == 0) { - - const Loop::GridDescBaseDevice grid(cctkGH); - grid.loop_int_device<0, 0, 0, vsize>( - grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE { - const vbool mask = mask_for_loop_tail(p.i, p.imax); - const vec betaG = gf_betaG_(mask, p.I); - const vreal rhs_old = gf_rhs_(mask, p.I); - const vreal rhs_new = - rhs_old + deriv_upwind(mask, gf_, p.I, betaG, dx); - gf_rhs_.store(mask, p.I, rhs_new); - }); - - } else { - - const Loop::GridDescBaseDevice grid(cctkGH); - grid.loop_int_device<0, 0, 0, vsize>( - grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE { - const vbool mask = mask_for_loop_tail(p.i, p.imax); - const vec betaG = gf_betaG_(mask, p.I); - const vreal rhs_old = gf_rhs_(mask, p.I); - const vreal rhs_new = rhs_old + - deriv_upwind(mask, gf_, p.I, betaG, dx) + - epsdiss * diss(mask, gf_, p.I, dx); - gf_rhs_.store(mask, p.I, rhs_new); - }); - } -} - -} // namespace Z4c - -#endif // #ifndef DERIVS_HXX diff --git a/Z4c/src/initial2.cxx b/Z4c/src/initial2.cxx index ec18e82a..1e303625 100644 --- a/Z4c/src/initial2.cxx +++ b/Z4c/src/initial2.cxx @@ -1,4 +1,4 @@ -#include "derivs.hxx" +#include "derivs_spacetimex.hxx" #include "physics.hxx" #include @@ -11,6 +11,8 @@ #include #endif +using namespace derivs_spacetimex; + namespace Z4c { using namespace Arith; using namespace Loop; diff --git a/Z4c/src/rhs.cxx b/Z4c/src/rhs.cxx index c8302786..8dde1ab5 100644 --- a/Z4c/src/rhs.cxx +++ b/Z4c/src/rhs.cxx @@ -8,7 +8,7 @@ #endif #endif -#include "derivs.hxx" +#include "derivs_spacetimex.hxx" #include "physics.hxx" #include "z4c_vars.hxx" @@ -27,6 +27,8 @@ #include +using namespace derivs_spacetimex; + namespace Z4c { using namespace Arith; using namespace Loop; diff --git a/Z4c/src/test.cxx b/Z4c/src/test.cxx index 939647f9..93066ea8 100644 --- a/Z4c/src/test.cxx +++ b/Z4c/src/test.cxx @@ -1,4 +1,4 @@ -#include "derivs.hxx" +#include "derivs_spacetimex.hxx" #include "physics.hxx" #include @@ -13,6 +13,8 @@ #include #include +using namespace derivs_spacetimex; + namespace Z4c { using namespace std; diff --git a/Z4c/src/z4c_vars.hxx b/Z4c/src/z4c_vars.hxx index 0926b0c8..e06ca45a 100644 --- a/Z4c/src/z4c_vars.hxx +++ b/Z4c/src/z4c_vars.hxx @@ -1,7 +1,7 @@ #ifndef Z4C_HXX #define Z4C_HXX -#include "derivs.hxx" +#include "derivs_spacetimex.hxx" #include "physics.hxx" #include @@ -11,6 +11,8 @@ #include #include +using namespace derivs_spacetimex; + namespace Z4c { // See arXiv:1212.2901 [gr-qc] @@ -819,123 +821,6 @@ template struct z4c_vars : z4c_vars_noderivs { // {} - template - ARITH_INLINE ARITH_DEVICE ARITH_HOST - z4c_vars(const T &kappa1, const T &kappa2, const T &f_mu_L, const T &f_mu_S, - const T &eta, - // - const GF3D3ptr &gf_chi_, - // - const GF3D3ptr &gf_gammatxx_, - const GF3D3ptr &gf_gammatxy_, - const GF3D3ptr &gf_gammatxz_, - const GF3D3ptr &gf_gammatyy_, - const GF3D3ptr &gf_gammatyz_, - const GF3D3ptr &gf_gammatzz_, - // - const GF3D3ptr &gf_Kh_, - // - const GF3D3ptr &gf_Atxx_, - const GF3D3ptr &gf_Atxy_, - const GF3D3ptr &gf_Atxz_, - const GF3D3ptr &gf_Atyy_, - const GF3D3ptr &gf_Atyz_, - const GF3D3ptr &gf_Atzz_, - // - const GF3D3ptr &gf_Gamtx_, - const GF3D3ptr &gf_Gamty_, - const GF3D3ptr &gf_Gamtz_, - // - const GF3D3ptr &gf_Theta_, - // - const GF3D3ptr &gf_alphaG_, - // - const GF3D3ptr &gf_betaGx_, - const GF3D3ptr &gf_betaGy_, - const GF3D3ptr &gf_betaGz_, - // - const GF3D3ptr &gf_eTtt_, - // - const GF3D3ptr &gf_eTtx_, - const GF3D3ptr &gf_eTty_, - const GF3D3ptr &gf_eTtz_, - // - const GF3D3ptr &gf_eTxx_, - const GF3D3ptr &gf_eTxy_, - const GF3D3ptr &gf_eTxz_, - const GF3D3ptr &gf_eTyy_, - const GF3D3ptr &gf_eTyz_, - const GF3D3ptr &gf_eTzz_, - // - const vect &I, const vec &dx) - : z4c_vars(kappa1, kappa2, f_mu_L, f_mu_S, eta, - // - gf_chi_(I), deriv(gf_chi_, I, dx), deriv2(gf_chi_, I, dx), - // - smat(gf_gammatxx_, gf_gammatxy_, gf_gammatxz_, - gf_gammatyy_, gf_gammatyz_, gf_gammatzz_, I), - smat, 3>{ - deriv(gf_gammatxx_, I, dx), - deriv(gf_gammatxy_, I, dx), - deriv(gf_gammatxz_, I, dx), - deriv(gf_gammatyy_, I, dx), - deriv(gf_gammatyz_, I, dx), - deriv(gf_gammatzz_, I, dx), - }, - smat, 3>{ - deriv2(gf_gammatxx_, I, dx), - deriv2(gf_gammatxy_, I, dx), - deriv2(gf_gammatxz_, I, dx), - deriv2(gf_gammatyy_, I, dx), - deriv2(gf_gammatyz_, I, dx), - deriv2(gf_gammatzz_, I, dx), - }, - // - gf_Kh_(I), deriv(gf_Kh_, I, dx), - // - smat(gf_Atxx_, gf_Atxy_, gf_Atxz_, gf_Atyy_, gf_Atyz_, - gf_Atzz_, I), - smat, 3>{ - deriv(gf_Atxx_, I, dx), - deriv(gf_Atxy_, I, dx), - deriv(gf_Atxz_, I, dx), - deriv(gf_Atyy_, I, dx), - deriv(gf_Atyz_, I, dx), - deriv(gf_Atzz_, I, dx), - }, - // - vec(gf_Gamtx_, gf_Gamty_, gf_Gamtz_, I), - vec, 3>{ - deriv(gf_Gamtx_, I, dx), - deriv(gf_Gamty_, I, dx), - deriv(gf_Gamtz_, I, dx), - }, - // - gf_Theta_(I), deriv(gf_Theta_, I, dx), - // - gf_alphaG_(I), deriv(gf_alphaG_, I, dx), - deriv2(gf_alphaG_, I, dx), - // - vec(gf_betaGx_, gf_betaGy_, gf_betaGz_, I), - vec, 3>{ - deriv(gf_betaGx_, I, dx), - deriv(gf_betaGy_, I, dx), - deriv(gf_betaGz_, I, dx), - }, - vec, 3>{ - deriv2(gf_betaGx_, I, dx), - deriv2(gf_betaGy_, I, dx), - deriv2(gf_betaGz_, I, dx), - }, - // - gf_eTtt_(I), - // - vec(gf_eTtx_, gf_eTty_, gf_eTtz_, I), - // - smat(gf_eTxx_, gf_eTxy_, gf_eTxz_, gf_eTyy_, gf_eTyz_, - gf_eTzz_, I)) - // - {} }; } // namespace Z4c diff --git a/scripts/spacetimex.th b/scripts/spacetimex.th index 0e44c484..994dce2a 100644 --- a/scripts/spacetimex.th +++ b/scripts/spacetimex.th @@ -695,7 +695,8 @@ Numerical/AEILocalInterp # CarpetX thorns !TARGET = $ARR !TYPE = git -!URL = https://github.com/eschnett/CarpetX.git +!URL = https://github.com/stevenrbrandt/CarpetX.git +!REPO_BRANCH = unify-derivs !REPO_PATH= $2 !CHECKOUT = CarpetX/ADMBaseX