TestRKAB: Code refactoring. Added noise initial data

TestRKAB/par/gaussian.par

@@ -25,6 +25,7 @@ $dt = $cfl * $dx
 $itlast = ceil($final_time / $dt)
 $out_every_it = ceil($out_every_time / $dt)
 
+TestRKAB::compute_error = yes
 TestRKAB::initial_condition = "Gaussian"
 TestRKAB::gaussian_width = $cfl / sqrt(2)
 

TestRKAB/par/noise.par

@@ -0,0 +1,79 @@
+ActiveThorns = "
+    CarpetX
+    IOUtil
+    ODESolvers
+    TestRKAB
+"
+
+###### Settings ######
+$final_time = 1.0
+$out_every_time = 0.25
+
+$start = -1.0
+$end = 1.0
+
+$ncells = 160
+$cfl = 0.25
+######################
+
+# We are using a 5 point centered stencil
+$ghost_size = 2
+
+$dx = ($end - $start)/$ncells
+$dt = $cfl * $dx
+
+$itlast = ceil($final_time / $dt)
+$out_every_it = ceil($out_every_time / $dt)
+
+TestRKAB::compute_error = no
+TestRKAB::initial_condition = "noise"
+TestRKAB::noise_seed = 100.0
+TestRKAB::noise_boundary = 1.0e-7
+
+CarpetX::poison_undefined_values = no
+CarpetX::verbose = no
+
+Cactus::cctk_show_schedule = yes
+Cactus::presync_mode = "mixed-error"
+
+CarpetX::xmin = $start
+CarpetX::ymin = $start
+CarpetX::zmin = $start
+
+CarpetX::xmax = $end
+CarpetX::ymax = $end
+CarpetX::zmax = $end
+
+CarpetX::ncells_x = $ncells
+CarpetX::ncells_y = $ncells
+CarpetX::ncells_z = $ncells
+
+CarpetX::periodic = yes
+CarpetX::periodic_x = yes
+CarpetX::periodic_y = yes
+CarpetX::periodic_z = yes
+
+CarpetX::ghost_size = $ghost_size
+
+CarpetX::dtfac = $cfl
+
+CarpetX::blocking_factor_x = 1
+CarpetX::blocking_factor_y = 1
+CarpetX::blocking_factor_z = 1
+
+Cactus::terminate = "iteration"
+Cactus::cctk_itlast = $itlast
+
+ODESolvers::method = "RKAB"
+ODESolvers::verbose = yes
+
+IO::out_dir = $parfile
+IO::out_every = $out_every_it
+
+CarpetX::out_norm_vars = ""
+CarpetX::out_tsv_vars = ""
+
+CarpetX::out_silo_vars = "
+    TestRKAB::state
+    TestRKAB::energy
+"

TestRKAB/par/standing.par

@@ -28,6 +28,7 @@ $dt = $cfl * $dx
 $itlast = ceil($final_time / $dt)
 $out_every_it = ceil($out_every_time / $dt)
 
+TestRKAB::compute_error = yes
 TestRKAB::initial_condition = "standing wave"
 TestRKAB::amplitude = 1.0
 TestRKAB::standing_wave_kx = $wave_k

TestRKAB/param.ccl

@@ -4,8 +4,13 @@ KEYWORD initial_condition "Initial condition"
 {
   "standing wave" :: "Standing wave"
   "Gaussian" :: "Gaussian"
+  "noise" :: "noise"
 } "standing wave"
 
+BOOLEAN compute_error "Turn on error computation for Gaussian and standing wave initial data"
+{
+} no
+
 CCTK_REAL amplitude "Initial amplitude"
 {
   *:* :: ""
@@ -29,4 +34,14 @@ CCTK_REAL standing_wave_kz "ky for standing wave"
 CCTK_REAL gaussian_width "width of Gaussian"
 {
   (0:* :: ""
-} 1.0
+} 1.0
+
+CCTK_REAL noise_seed "Seed of the RNG for noise initial data"
+{
+  *:* :: ""
+} 100.0
+
+CCTK_REAL noise_boundary "Nois will be generated in the [-noise_boundary, noise_boundary] range"
+{
+  *:* :: ""
+} 1.0e-7

TestRKAB/schedule.ccl

@@ -37,13 +37,16 @@ SCHEDULE TestRKAB_Sync IN ODESolvers_PostStep
 
 
 # Analysis
-SCHEDULE TestRKAB_Error IN ODESolvers_PostStep AFTER TestRKAB_Sync
+if (CCTK_EQUALS(initial_condition, "Gaussian") || CCTK_EQUALS(initial_condition, "standing wave"))
 {
-  LANG: C
-  READS: state(interior)
-  WRITES: error(interior)
-  SYNC: error
-} "Calculate error in scalar wave state"
+  SCHEDULE TestRKAB_Error IN ODESolvers_PostStep AFTER TestRKAB_Sync
+  {
+    LANG: C
+    READS: state(interior)
+    WRITES: error(interior)
+    SYNC: error
+  } "Calculate error in scalar wave state"
+}
 
 SCHEDULE TestRKAB_Energy IN ODESolvers_PostStep AFTER TestRKAB_Sync
 {

TestRKAB/src/energy.cxx

@@ -0,0 +1,27 @@
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+#include <loop_device.hxx>
+
+namespace TestRKAB {
+using namespace Arith;
+
+extern "C" void TestRKAB_Energy(CCTK_ARGUMENTS) {
+  DECLARE_CCTK_ARGUMENTSX_TestRKAB_Energy;
+  DECLARE_CCTK_PARAMETERS;
+
+  grid.loop_int_device<0, 0, 0>(
+      grid.nghostzones,
+      [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+        const auto local_Pi{Pi(p.I)};
+        const auto local_Dx{Dx(p.I)};
+        const auto local_Dy{Dy(p.I)};
+        const auto local_Dz{Dz(p.I)};
+
+        energy_density(p.I) = 0.5 * (local_Pi * local_Pi + local_Dx * local_Dx +
+                                     local_Dy * local_Dy + local_Dz * local_Dz);
+      });
+}
+
+} // namespace TestRKAB

TestRKAB/src/error.cxx

@@ -0,0 +1,85 @@
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+#include <loop_device.hxx>
+
+#include "standing_wave.hxx"
+#include "gaussian.hxx"
+
+namespace TestRKAB {
+using namespace Arith;
+
+extern "C" void TestRKAB_Error(CCTK_ARGUMENTS) {
+  DECLARE_CCTK_ARGUMENTSX_TestRKAB_Error;
+  DECLARE_CCTK_PARAMETERS;
+
+  if (CCTK_EQUALS(initial_condition, "standing wave")) {
+    grid.loop_int_device<0, 0, 0>(
+        grid.nghostzones,
+        [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+          using std::fabs;
+
+          const auto t{cctk_time};
+
+          const auto A{amplitude};
+          const auto kx{standing_wave_kx};
+          const auto ky{standing_wave_ky};
+          const auto kz{standing_wave_kz};
+
+          const auto expected_phi{sw::phi(A, kx, ky, kz, t, p.x, p.y, p.z)};
+          const auto expected_Pi{sw::Pi(A, kx, ky, kz, t, p.x, p.y, p.z)};
+          const auto expected_Dx{sw::Dx(A, kx, ky, kz, t, p.x, p.y, p.z)};
+          const auto expected_Dy{sw::Dy(A, kx, ky, kz, t, p.x, p.y, p.z)};
+          const auto expected_Dz{sw::Dz(A, kx, ky, kz, t, p.x, p.y, p.z)};
+
+          const auto actual_phi{phi(p.I)};
+          const auto actual_Pi{Pi(p.I)};
+          const auto actual_Dx{Dx(p.I)};
+          const auto actual_Dy{Dy(p.I)};
+          const auto actual_Dz{Dz(p.I)};
+
+          phi_err(p.I) = fabs(expected_phi - actual_phi);
+          Pi_err(p.I) = fabs(expected_Pi - actual_Pi);
+          Dx_err(p.I) = fabs(expected_Dx - actual_Dx);
+          Dy_err(p.I) = fabs(expected_Dy - actual_Dy);
+          Dz_err(p.I) = fabs(expected_Dz - actual_Dz);
+        });
+
+  } else if (CCTK_EQUALS(initial_condition, "Gaussian")) {
+    grid.loop_int_device<0, 0, 0>(
+        grid.nghostzones,
+        [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+          using std::fabs;
+
+          const auto t{cctk_time};
+
+          const auto expected_phi{
+              gauss::phi(amplitude, gaussian_width, t, p.x, p.y, p.z)};
+          const auto expected_Pi{
+              gauss::Pi(amplitude, gaussian_width, t, p.x, p.y, p.z)};
+          const auto expected_Dx{
+              gauss::Dx(amplitude, gaussian_width, t, p.x, p.y, p.z)};
+          const auto expected_Dy{
+              gauss::Dy(amplitude, gaussian_width, t, p.x, p.y, p.z)};
+          const auto expected_Dz{
+              gauss::Dz(amplitude, gaussian_width, t, p.x, p.y, p.z)};
+
+          const auto actual_phi{phi(p.I)};
+          const auto actual_Pi{Pi(p.I)};
+          const auto actual_Dx{Dx(p.I)};
+          const auto actual_Dy{Dy(p.I)};
+          const auto actual_Dz{Dz(p.I)};
+
+          phi_err(p.I) = fabs(expected_phi - actual_phi);
+          Pi_err(p.I) = fabs(expected_Pi - actual_Pi);
+          Dx_err(p.I) = fabs(expected_Dx - actual_Dx);
+          Dy_err(p.I) = fabs(expected_Dy - actual_Dy);
+          Dz_err(p.I) = fabs(expected_Dz - actual_Dz);
+        });
+  } else {
+    CCTK_VERROR("Unknown initial condition \"%s\"", initial_condition);
+  }
+}
+
+} // namespace TestRKAB

TestRKAB/src/initial.cxx

@@ -0,0 +1,71 @@
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+#include <loop_device.hxx>
+
+#include "standing_wave.hxx"
+#include "gaussian.hxx"
+
+#include <random>
+
+namespace TestRKAB {
+using namespace Arith;
+
+extern "C" void TestRKAB_Initial(CCTK_ARGUMENTS) {
+  DECLARE_CCTK_ARGUMENTSX_TestRKAB_Initial;
+  DECLARE_CCTK_PARAMETERS;
+
+  if (CCTK_EQUALS(initial_condition, "standing wave")) {
+    grid.loop_int_device<0, 0, 0>(
+        grid.nghostzones,
+        [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+          const auto t{cctk_time};
+
+          const auto A{amplitude};
+          const auto kx{standing_wave_kx};
+          const auto ky{standing_wave_ky};
+          const auto kz{standing_wave_kz};
+
+          phi(p.I) = sw::phi(A, kx, ky, kz, t, p.x, p.y, p.z);
+          Pi(p.I) = sw::Pi(A, kx, ky, kz, t, p.x, p.y, p.z);
+          Dx(p.I) = sw::Dx(A, kx, ky, kz, t, p.x, p.y, p.z);
+          Dy(p.I) = sw::Dy(A, kx, ky, kz, t, p.x, p.y, p.z);
+          Dz(p.I) = sw::Dz(A, kx, ky, kz, t, p.x, p.y, p.z);
+        });
+
+  } else if (CCTK_EQUALS(initial_condition, "Gaussian")) {
+    grid.loop_int_device<0, 0, 0>(
+        grid.nghostzones,
+        [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+          const auto t{cctk_time};
+
+          phi(p.I) = gauss::phi(amplitude, gaussian_width, t, p.x, p.y, p.z);
+          Pi(p.I) = gauss::Pi(amplitude, gaussian_width, t, p.x, p.y, p.z);
+          Dx(p.I) = gauss::Dx(amplitude, gaussian_width, t, p.x, p.y, p.z);
+          Dy(p.I) = gauss::Dy(amplitude, gaussian_width, t, p.x, p.y, p.z);
+          Dz(p.I) = gauss::Dz(amplitude, gaussian_width, t, p.x, p.y, p.z);
+        });
+
+  } else if (CCTK_EQUALS(initial_condition, "noise")) {
+    std::uniform_real_distribution<CCTK_REAL> noise_distrib{-noise_boundary,
+                                                            noise_boundary};
+    std::mt19937 noise_engine{noise_seed};
+
+    grid.loop_int_device<0, 0, 0>(grid.nghostzones,
+                                  [&] CCTK_DEVICE(const Loop::PointDesc &p)
+                                      CCTK_ATTRIBUTE_ALWAYS_INLINE {
+                                        const auto t{cctk_time};
+
+                                        phi(p.I) = noise_distrib(noise_engine);
+                                        Pi(p.I) = noise_distrib(noise_engine);
+                                        Dx(p.I) = noise_distrib(noise_engine);
+                                        Dy(p.I) = noise_distrib(noise_engine);
+                                        Dz(p.I) = noise_distrib(noise_engine);
+                                      });
+  } else {
+    CCTK_VERROR("Unknown initial condition \"%s\"", initial_condition);
+  }
+}
+
+} // namespace TestRKAB

TestRKAB/src/make.code.defn

@@ -1,7 +1,7 @@
 # Main make.code.defn file for thorn TestSubcyclingMC
 
 # Source files in this directory
-SRCS = testrkab.cxx
+SRCS = energy.cxx error.cxx initial.cxx rhs.cxx
 
 # Subdirectories containing source files
 SUBDIRS =

TestRKAB/src/rhs.cxx

@@ -0,0 +1,44 @@
+
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+#include <loop_device.hxx>
+
+namespace TestRKAB {
+using namespace Arith;
+
+// Finite difference helpers
+enum class fd_dir : std::size_t { x = 0, y = 1, z = 2 };
+
+template <fd_dir direction, typename T>
+static inline auto fd_c_1_4(const Loop::PointDesc &p,
+                            const Loop::GF3D2<T> &gf) noexcept -> T {
+  constexpr auto d{static_cast<size_t>(direction)};
+  const auto num{gf(p.I - 2 * p.DI[d]) - 8.0 * gf(p.I - 1 * p.DI[d]) +
+                 8.0 * gf(p.I + 1 * p.DI[d]) - 1.0 * gf(p.I + 2 * p.DI[d])};
+  const auto den{1.0 / (12.0 * p.DX[d])};
+  return num * den;
+}
+
+extern "C" void TestRKAB_RHS(CCTK_ARGUMENTS) {
+  DECLARE_CCTK_ARGUMENTSX_TestRKAB_RHS;
+  DECLARE_CCTK_PARAMETERS;
+
+  grid.loop_int_device<0, 0, 0>(
+      grid.nghostzones,
+      [=] CCTK_DEVICE(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+        phi_rhs(p.I) = Pi(p.I);
+        Pi_rhs(p.I) = fd_c_1_4<fd_dir::x>(p, Dx) + fd_c_1_4<fd_dir::y>(p, Dy) +
+                      fd_c_1_4<fd_dir::z>(p, Dz);
+        Dx_rhs(p.I) = fd_c_1_4<fd_dir::x>(p, Pi);
+        Dy_rhs(p.I) = fd_c_1_4<fd_dir::y>(p, Pi);
+        Dz_rhs(p.I) = fd_c_1_4<fd_dir::z>(p, Pi);
+      });
+}
+
+extern "C" void TestRKAB_Sync(CCTK_ARGUMENTS) {
+  // do nothing
+}
+
+} // namespace TestRKAB