Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

TestDerivs: Add 6th order central finite difference #204

Merged
merged 1 commit into from
Aug 17, 2023
Merged

TestDerivs: Add 6th order central finite difference #204

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
50 changes: 50 additions & 0 deletions Derivs/src/derivs.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,22 @@ calc_derivs(const GF3D5<T> &gf, const vec<GF3D5<T>, dim> &dgf,
});
break;

case 6:
grid.loop_int_device<CI, CJ, CK, vsize>(
grid.nghostzones,
[=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);
// Take account of ghost points
const vbool mask1 =
mask_for_loop_tail<vbool>(p.i, p.imax + deriv_order / 2);
const GF3D5index index(layout, p.I);
const auto val = gf0(mask, p.I);
const auto dval = calc_deriv<6>(gf0, mask1, p.I, dx);
gf.store(mask, index, val);
dgf.store(mask, index, dval);
});
break;

default:
CCTK_VERROR("Unsupported derivative order %d", deriv_order);
}
Expand Down Expand Up @@ -124,6 +140,24 @@ calc_derivs2(const GF3D5<T> &gf, const vec<GF3D5<T>, dim> &dgf,
});
break;

case 6:
grid.loop_int_device<CI, CJ, CK, vsize>(
grid.nghostzones,
[=] CCTK_DEVICE(const PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
const vbool mask = mask_for_loop_tail<vbool>(p.i, p.imax);
// Take account of ghost points
const vbool mask1 =
mask_for_loop_tail<vbool>(p.i, p.imax + deriv_order / 2);
const GF3D5index index(layout, p.I);
const auto val = gf0(mask, p.I);
const auto dval = calc_deriv<6>(gf0, mask1, p.I, dx);
const auto ddval = calc_deriv2<6>(gf0, mask1, p.I, dx);
gf.store(mask, index, val);
dgf.store(mask, index, dval);
ddgf.store(mask, index, ddval);
});
break;

default:
CCTK_VERROR("Unsupported derivative order %d", deriv_order);
}
Expand All @@ -143,6 +177,10 @@ template CCTK_DEVICE CCTK_HOST Arith::vec<Arith::simd<T>, Loop::dim>
calc_deriv<4>(const Loop::GF3D2<const T> &gf, const Arith::simdl<T> &mask,
const Arith::vect<int, Loop::dim> &I,
const Arith::vect<T, Loop::dim> &dx);
template CCTK_DEVICE CCTK_HOST Arith::vec<Arith::simd<T>, Loop::dim>
calc_deriv<6>(const Loop::GF3D2<const T> &gf, const Arith::simdl<T> &mask,
const Arith::vect<int, Loop::dim> &I,
const Arith::vect<T, Loop::dim> &dx);

template CCTK_DEVICE CCTK_HOST Arith::vec<T, Loop::dim>
calc_deriv<2>(const Loop::GF3D2<const T> &gf,
Expand All @@ -152,6 +190,10 @@ template CCTK_DEVICE CCTK_HOST Arith::vec<T, Loop::dim>
calc_deriv<4>(const Loop::GF3D2<const T> &gf,
const Arith::vect<int, Loop::dim> &I,
const Arith::vect<T, Loop::dim> &dx);
template CCTK_DEVICE CCTK_HOST Arith::vec<T, Loop::dim>
calc_deriv<6>(const Loop::GF3D2<const T> &gf,
const Arith::vect<int, Loop::dim> &I,
const Arith::vect<T, Loop::dim> &dx);

template CCTK_DEVICE CCTK_HOST Arith::smat<Arith::simd<T>, Loop::dim>
calc_deriv2<2>(const Loop::GF3D2<const T> &gf, const Arith::simdl<T> &mask,
Expand All @@ -161,6 +203,10 @@ template CCTK_DEVICE CCTK_HOST Arith::smat<Arith::simd<T>, Loop::dim>
calc_deriv2<4>(const Loop::GF3D2<const T> &gf, const Arith::simdl<T> &mask,
const Arith::vect<int, Loop::dim> &I,
const Arith::vect<T, Loop::dim> &dx);
template CCTK_DEVICE CCTK_HOST Arith::smat<Arith::simd<T>, Loop::dim>
calc_deriv2<6>(const Loop::GF3D2<const T> &gf, const Arith::simdl<T> &mask,
const Arith::vect<int, Loop::dim> &I,
const Arith::vect<T, Loop::dim> &dx);

template CCTK_DEVICE CCTK_HOST Arith::smat<T, Loop::dim>
calc_deriv2<2>(const Loop::GF3D2<const T> &gf,
Expand All @@ -170,6 +216,10 @@ template CCTK_DEVICE CCTK_HOST Arith::smat<T, Loop::dim>
calc_deriv2<4>(const Loop::GF3D2<const T> &gf,
const Arith::vect<int, Loop::dim> &I,
const Arith::vect<T, Loop::dim> &dx);
template CCTK_DEVICE CCTK_HOST Arith::smat<T, Loop::dim>
calc_deriv2<6>(const Loop::GF3D2<const T> &gf,
const Arith::vect<int, Loop::dim> &I,
const Arith::vect<T, Loop::dim> &dx);

template void calc_copy<0, 0, 0>(const GF3D5<T> &gf, const GF3D5layout layout,
const GridDescBaseDevice &grid,
Expand Down
46 changes: 40 additions & 6 deletions Derivs/src/derivs.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -130,6 +130,18 @@ inline CCTK_ATTRIBUTE_ALWAYS_INLINE
return c2 * (var(2) - var(-2)) + c1 * (var(1) - var(-1));
}

template <int deriv_order, typename T, typename TS,
typename R = std::result_of_t<TS(int)> >
inline CCTK_ATTRIBUTE_ALWAYS_INLINE
CCTK_DEVICE CCTK_HOST std::enable_if_t<deriv_order == 6, R>
deriv1d(const TS var, const T dx) {
const T c1 = 3 / (4 * dx);
const T c2 = -3 / (20 * dx);
const T c3 = 1 / (60 * dx);
return c3 * (var(3) - var(-3)) + c2 * (var(2) - var(-2)) +
c1 * (var(1) - var(-1));
}

template <int deriv_order, typename T>
inline CCTK_ATTRIBUTE_ALWAYS_INLINE
CCTK_DEVICE CCTK_HOST std::enable_if_t<deriv_order == 2, simd<T> >
Expand Down Expand Up @@ -234,8 +246,8 @@ inline CCTK_ATTRIBUTE_ALWAYS_INLINE
// constexpr T c0 = -2 / pow2(dx);
// constexpr T c1 = 1 / pow2(dx);
// return c1 * (var(-1) + var(1)) + c0 * var(0);
const T c0 = 1 / pow2(dx);
return c0 * ((var(1) - var(0)) - (var(0) - var(-1)));
const T c1 = 1 / pow2(dx);
return c1 * ((var(1) - var(0)) - (var(0) - var(-1)));
}

template <int deriv_order, typename T, typename TS,
Expand All @@ -248,10 +260,32 @@ inline CCTK_ATTRIBUTE_ALWAYS_INLINE
// constexpr T c2 = -1 / T(12);
// return (c2 * (var(-2) + var(2)) + c1 * (var(-1) + var(1)) + c0 * var(0)) /
// pow2(dx);
constexpr T c0 = 4 / T(3);
constexpr T c1 = -1 / T(12);
return (c1 * ((var(+2) - var(+0)) - (var(-0) - var(-2))) +
c0 * ((var(+1) - var(+0)) - (var(-0) - var(-1)))) /
constexpr T c1 = 4 / T(3);
constexpr T c2 = -1 / T(12);
return (c2 * ((var(+2) - var(+0)) - (var(-0) - var(-2))) +
c1 * ((var(+1) - var(+0)) - (var(-0) - var(-1)))) /
pow2(dx);
}

template <int deriv_order, typename T, typename TS,
typename R = std::result_of_t<TS(int)> >
inline CCTK_ATTRIBUTE_ALWAYS_INLINE
CCTK_DEVICE CCTK_HOST std::enable_if_t<deriv_order == 6, R>
deriv2_1d(const TS var, const T dx) {
// constexpr T c0 = -49 / T(18);
// constexpr T c1 = 3 / T(2);
// constexpr T c2 = -3 / T(20);
// constexpr T c3 = 1 / T(90);
// return (c3 * (var(-3) + var(3)) + c2 * (var(-2) + var(2)) +
// c1 * (var(-1) + var(1)) + c0 * var(0)) /
// pow2(dx);
// constexpr T c0 = -49 / T(18);
constexpr T c1 = 3 / T(2);
constexpr T c2 = -3 / T(20);
constexpr T c3 = 1 / T(90);
return (c3 * ((var(+3) - var(+0)) - (var(-0) - var(-3))) +
c2 * ((var(+2) - var(+0)) - (var(-0) - var(-2))) +
c1 * ((var(+1) - var(+0)) - (var(-0) - var(-1)))) /
pow2(dx);
}

Expand Down
1 change: 1 addition & 0 deletions TestDerivs/param.ccl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ CCTK_INT deriv_order "Order of spatial finite differencing" STEERABLE=never
{
2 :: "Second order finite difference"
4 :: "Fourth order finite difference"
6 :: "Sixth order finite difference"
} 4

CCTK_REAL kxx "par for polynomial"
Expand Down