TestRKAB: initialize pre using analytic expression

EinsteinToolkit · Sep 20, 2024 · cbcea62 · cbcea62
1 parent 2586ba9
commit cbcea62
Show file tree

Hide file tree

Showing 2 changed files with 27 additions and 2 deletions.
diff --git a/TestRKAB/schedule.ccl b/TestRKAB/schedule.ccl
@@ -6,8 +6,8 @@ STORAGE: state pre rhs
 SCHEDULE TestRKAB_Initial AT initial
 {
   LANG: C
-  WRITES: state(interior)
-  SYNC: state
+  WRITES: state(interior) pre(interior)
+  SYNC: state pre
 } "Initialize scalar wave state"
 
 SCHEDULE TestRKAB_Sync AT postregrid

diff --git a/TestRKAB/src/testrkab.cxx b/TestRKAB/src/testrkab.cxx
@@ -56,10 +56,33 @@ constexpr void gaussian(const T A, const T W, const T t, const T x, const T y,
   }
 }
 
+template <typename T>
+constexpr void dtgaussian(const T A, const T W, const T t, const T x, const T y,
+                          const T z, T &dtu, T &dtrho) {
+  using std::exp, std::pow, std::sqrt;
+
+  const T r = sqrt(pow(x, 2) + pow(y, 2) + pow(z, 2));
+  const auto f = [&](const T v) {
+    return A * exp(-pow(v, 2) / (2 * pow(W, 2)));
+  };
+
+  if (r < sqrt(std::numeric_limits<T>::epsilon())) {
+    dtu = -2 / pow(W, 4) * f(t) * (pow(t, 2) - pow(W, 2));
+    dtrho = (2 * t * (pow(t, 2) - 3 * pow(W, 2)) * f(t)) / pow(W, 6);
+  } else {
+    dtu = -(f(t - r) * (t - r) - f(t + r) * (t + r)) / (pow(W, 2) * r);
+    dtrho = -(f(t - r) * (1 - pow(t - r, 2) / pow(W, 2)) -
+              f(t + r) * (1 - pow(t + r, 2) / pow(W, 2))) /
+            (pow(W, 2) * r);
+  }
+}
+
 extern "C" void TestRKAB_Initial(CCTK_ARGUMENTS) {
   DECLARE_CCTK_ARGUMENTSX_TestRKAB_Initial;
   DECLARE_CCTK_PARAMETERS;
 
+  const CCTK_REAL dt = cctk_delta_time;
+
   if (CCTK_EQUALS(initial_condition, "standing wave")) {
     grid.loop_int_device<0, 0, 0>(
         grid.nghostzones,
@@ -75,6 +98,8 @@ extern "C" void TestRKAB_Initial(CCTK_ARGUMENTS) {
         [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
           gaussian(amplitude, gaussian_width, cctk_time, p.x, p.y, p.z, u(p.I),
                    rho(p.I));
+          dtgaussian(amplitude, gaussian_width, cctk_time - dt, p.x, p.y, p.z,
+                     u_pre(p.I), rho_pre(p.I));
         });
 
   } else {