Derivs: add 8th order finite difference derivative operator

lwJi · Sep 6, 2024 · 9bbcf5a · 9bbcf5a
1 parent 293dcec
commit 9bbcf5a
Show file tree

Hide file tree

Showing 3 changed files with 64 additions and 0 deletions.
diff --git a/Derivs/src/derivs.cxx b/Derivs/src/derivs.cxx
@@ -104,6 +104,22 @@ calc_derivs(const GF3D5<T> &gf, const vec<GF3D5<T>, dim> &dgf,
         });
     break;
 
+  case 8:
+    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<8>(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);
   }
@@ -199,6 +215,24 @@ calc_derivs2(const GF3D5<T> &gf, const vec<GF3D5<T>, dim> &dgf,
         });
     break;
 
+  case 8:
+    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<8>(gf0, mask1, p.I, dx);
+          const auto ddval = calc_deriv2<8>(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);
   }

diff --git a/Derivs/src/derivs.hxx b/Derivs/src/derivs.hxx
@@ -142,6 +142,19 @@ inline CCTK_ATTRIBUTE_ALWAYS_INLINE
          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 == 8, R>
+    deriv1d(const TS var, const T dx) {
+  const T c1 = 4 / (5 * dx);
+  const T c2 = -1 / (5 * dx);
+  const T c3 = 4 / (105 * dx);
+  const T c4 = -1 / (280 * dx);
+  return c4 * (var(4) - var(-4)) + c3 * (var(3) - var(-3)) +
+         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
@@ -218,6 +231,22 @@ inline CCTK_ATTRIBUTE_ALWAYS_INLINE
          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 == 8, R>
+    deriv2_1d(const TS var, const T dx) {
+  constexpr T c1 = 8 / T(5);
+  constexpr T c2 = -1 / T(5);
+  constexpr T c3 = 8 / T(315);
+  constexpr T c4 = -1 / T(560);
+  return (c4 * ((var(+4) - var(+0)) - (var(-0) - var(-4))) +
+          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);
+}
+
 template <int deriv_order, bool vectorize_di, typename T, typename TS,
           typename R = std::result_of_t<TS(int, int)>>
 inline CCTK_ATTRIBUTE_ALWAYS_INLINE CCTK_DEVICE CCTK_HOST R

diff --git a/TestDerivs/param.ccl b/TestDerivs/param.ccl
@@ -10,6 +10,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"
+  8 :: "Eighth order finite difference"
 } 4
 
 CCTK_REAL kxx "par for polynomial"