New thorn: NoiseX

ADMBaseX/configuration.ccl

@@ -1,3 +1,3 @@
 # Configuration definition for thorn ADMBaseX
 
-REQUIRES Loop
+REQUIRES Loop NoiseX

ADMBaseX/interface.ccl

@@ -3,6 +3,7 @@
 IMPLEMENTS: ADMBaseX
 
 USES INCLUDE HEADER: loop_device.hxx
+USES INCLUDE HEADER: noisex.hxx
 
 
 

ADMBaseX/param.ccl

@@ -47,7 +47,6 @@ CCTK_REAL linear_wave_wavelength "Linear wave wavelength"
   0.0:* :: ""
 } 1.0
 
-CCTK_REAL noise_amplitude "Noise amplitude"
+BOOLEAN add_noise "Add noise to initial data"
 {
-  0.0:* :: ""
-} 0.0
+} no

ADMBaseX/schedule.ccl

@@ -106,9 +106,8 @@ if (CCTK_EQUALS(initial_dtshift, "zero")) {
 
 
 
-if (noise_amplitude != 0) {
+if (add_noise) {
   # TODO: Also add noise during evolution?
-  # TODO: Noise should be added by a separate thorn.
 
   SCHEDULE ADMBaseX_add_noise IN ADMBaseX_PostInitial
   {

ADMBaseX/src/noise.cxx

@@ -4,85 +4,40 @@
 #include <cctk_Arguments.h>
 #include <cctk_Parameters.h>
 
-#include <random>
+#include <noisex.hxx>
 
 namespace ADMBaseX {
 using namespace Loop;
-using namespace std;
 
 extern "C" void ADMBaseX_add_noise(CCTK_ARGUMENTS) {
-  DECLARE_CCTK_ARGUMENTS_ADMBaseX_add_noise;
-  DECLARE_CCTK_PARAMETERS;
-
-  // Hardware random device
-  random_device device;
-  // Create and seed software random number engine from hardware random number
-  default_random_engine engine(device());
-  // Random number distribution
-  uniform_real_distribution<CCTK_REAL> distribution(-noise_amplitude,
-                                                    noise_amplitude);
-  const auto add_noise = [&](CCTK_REAL &restrict var) {
-    var += distribution(engine);
-  };
-
-  const GF3D<CCTK_REAL, 0, 0, 0> gxx_(cctkGH, gxx);
-  const GF3D<CCTK_REAL, 0, 0, 0> gxy_(cctkGH, gxy);
-  const GF3D<CCTK_REAL, 0, 0, 0> gxz_(cctkGH, gxz);
-  const GF3D<CCTK_REAL, 0, 0, 0> gyy_(cctkGH, gyy);
-  const GF3D<CCTK_REAL, 0, 0, 0> gyz_(cctkGH, gyz);
-  const GF3D<CCTK_REAL, 0, 0, 0> gzz_(cctkGH, gzz);
-
-  const GF3D<CCTK_REAL, 0, 0, 0> kxx_(cctkGH, kxx);
-  const GF3D<CCTK_REAL, 0, 0, 0> kxy_(cctkGH, kxy);
-  const GF3D<CCTK_REAL, 0, 0, 0> kxz_(cctkGH, kxz);
-  const GF3D<CCTK_REAL, 0, 0, 0> kyy_(cctkGH, kyy);
-  const GF3D<CCTK_REAL, 0, 0, 0> kyz_(cctkGH, kyz);
-  const GF3D<CCTK_REAL, 0, 0, 0> kzz_(cctkGH, kzz);
-
-  const GF3D<CCTK_REAL, 0, 0, 0> alp_(cctkGH, alp);
-
-  const GF3D<CCTK_REAL, 0, 0, 0> dtalp_(cctkGH, dtalp);
-
-  const GF3D<CCTK_REAL, 0, 0, 0> betax_(cctkGH, betax);
-  const GF3D<CCTK_REAL, 0, 0, 0> betay_(cctkGH, betay);
-  const GF3D<CCTK_REAL, 0, 0, 0> betaz_(cctkGH, betaz);
-
-  const GF3D<CCTK_REAL, 0, 0, 0> dtbetax_(cctkGH, dtbetax);
-  const GF3D<CCTK_REAL, 0, 0, 0> dtbetay_(cctkGH, dtbetay);
-  const GF3D<CCTK_REAL, 0, 0, 0> dtbetaz_(cctkGH, dtbetaz);
-
-  loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) { add_noise(gxx_(p.I)); });
-  loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) { add_noise(gxy_(p.I)); });
-  loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) { add_noise(gxz_(p.I)); });
-  loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) { add_noise(gyy_(p.I)); });
-  loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) { add_noise(gyz_(p.I)); });
-  loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) { add_noise(gzz_(p.I)); });
-
-  loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) { add_noise(kxx_(p.I)); });
-  loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) { add_noise(kxy_(p.I)); });
-  loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) { add_noise(kxz_(p.I)); });
-  loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) { add_noise(kyy_(p.I)); });
-  loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) { add_noise(kyz_(p.I)); });
-  loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) { add_noise(kzz_(p.I)); });
-
-  loop_int<0, 0, 0>(cctkGH, [&](const PointDesc &p) { add_noise(alp_(p.I)); });
-
-  loop_int<0, 0, 0>(cctkGH,
-                    [&](const PointDesc &p) { add_noise(dtalp_(p.I)); });
-
-  loop_int<0, 0, 0>(cctkGH,
-                    [&](const PointDesc &p) { add_noise(betax_(p.I)); });
-  loop_int<0, 0, 0>(cctkGH,
-                    [&](const PointDesc &p) { add_noise(betay_(p.I)); });
-  loop_int<0, 0, 0>(cctkGH,
-                    [&](const PointDesc &p) { add_noise(betaz_(p.I)); });
-
-  loop_int<0, 0, 0>(cctkGH,
-                    [&](const PointDesc &p) { add_noise(dtbetax_(p.I)); });
-  loop_int<0, 0, 0>(cctkGH,
-                    [&](const PointDesc &p) { add_noise(dtbetay_(p.I)); });
-  loop_int<0, 0, 0>(cctkGH,
-                    [&](const PointDesc &p) { add_noise(dtbetaz_(p.I)); });
+  DECLARE_CCTK_ARGUMENTSX_ADMBaseX_add_noise;
+
+  NoiseX::add(cctkGH, gxx);
+  NoiseX::add(cctkGH, gxx);
+  NoiseX::add(cctkGH, gxy);
+  NoiseX::add(cctkGH, gxz);
+  NoiseX::add(cctkGH, gyy);
+  NoiseX::add(cctkGH, gyz);
+  NoiseX::add(cctkGH, gzz);
+
+  NoiseX::add(cctkGH, kxx);
+  NoiseX::add(cctkGH, kxy);
+  NoiseX::add(cctkGH, kxz);
+  NoiseX::add(cctkGH, kyy);
+  NoiseX::add(cctkGH, kyz);
+  NoiseX::add(cctkGH, kzz);
+
+  NoiseX::add(cctkGH, alp);
+
+  NoiseX::add(cctkGH, dtalp);
+
+  NoiseX::add(cctkGH, betax);
+  NoiseX::add(cctkGH, betay);
+  NoiseX::add(cctkGH, betaz);
+
+  NoiseX::add(cctkGH, dtbetax);
+  NoiseX::add(cctkGH, dtbetay);
+  NoiseX::add(cctkGH, dtbetaz);
 }
 
 } // namespace ADMBaseX

NoiseX/README

@@ -0,0 +1,11 @@
+Cactus Code Thorn Noise
+Author(s)    : Lucas Timotheo Sanches <lsanches@lsu.edu>
+Maintainer(s): Lucas Timotheo Sanches <lsanches@lsu.edu>
+Licence      : GPL
+--------------------------------------------------------------------------
+
+1. Purpose
+
+This thorn provides a unified way to add pseudo-random noise to CarpetX
+variables. This is useful mainly for testing time evolution codes, i.e.,
+performing robust stability tests.

NoiseX/configuration.ccl

@@ -0,0 +1,7 @@
+# Configuration definition for thorn NoiseX
+
+REQUIRES Loop
+
+PROVIDES NoiseX
+{
+}

NoiseX/interface.ccl

@@ -0,0 +1,7 @@
+# Interface definition for thorn NoiseX
+
+IMPLEMENTS: NoiseX
+
+USES INCLUDE HEADER: loop.hxx
+
+INCLUDE HEADER: noisex.hxx in noisex.hxx

NoiseX/param.ccl

@@ -0,0 +1,29 @@
+# Parameter definitions for thorn NoiseX
+
+KEYWORD seed_type "The type of seed to use"
+{
+  "random" :: "Uses the hardware random device to generate a seed"
+  "fixed" :: "A fixed integer value is used as sseed"
+} "random"
+
+KEYWORD random_engine "The random engine to use"
+{
+  "default" :: "Uses C++'s default random engine"
+  "mersenne twister" :: "32-bit Mersenne Twister by Matsumoto and Nishimura, 1998"
+  "sine wave" :: "Uses the values of a high frequency sine wave to create pseudo-randomness"
+} "default"
+
+INT fixed_seed_value "The value of the fixed seed to use"
+{
+  0:* :: "Positive integer"
+} 100
+
+REAL noise_amplitude "The value of the noise amplitude to use"
+{
+  0.0:* :: "Positive real"
+} 1.0e-10
+
+REAL noise_frequency "The value of the noise frenquecy to use with the sine wave engie"
+{
+  0.0:* :: "Positive real"
+} 5.0

NoiseX/schedule.ccl

@@ -0,0 +1,6 @@
+# Schedule definitions for thorn NoiseX
+
+SCHEDULE NoiseX_setup_globals AT startup BEFORE Driver_Startup
+{
+  LANG: C
+} "Create NoiseX global data"

NoiseX/src/make.code.defn

@@ -0,0 +1,7 @@
+# Main make.code.defn file for thorn NoiseX
+
+# Source files in this directory
+SRCS = noisex.cxx
+
+# Subdirectories containing source files
+SUBDIRS =

NoiseX/src/noisex.cxx

@@ -0,0 +1,78 @@
+#include <loop.hxx>
+
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+#include "noisex.hxx"
+
+#include <cmath>
+#include <random>
+
+namespace NoiseX {
+
+struct NoiseData {
+  std::default_random_engine default_engine{};
+  std::mt19937 mt_engine{};
+  std::uniform_real_distribution<CCTK_REAL> distrib{};
+};
+
+static NoiseData noise_data{};
+
+extern "C" int NoiseX_setup_globals() {
+  DECLARE_CCTK_PARAMETERS;
+
+  if (CCTK_EQUALS(seed_type, "random")) {
+    std::random_device device{};
+    const auto random_seed_value{device()};
+    noise_data.default_engine = std::default_random_engine(random_seed_value);
+    noise_data.mt_engine = std::mt19937(random_seed_value);
+  } else if (CCTK_EQUALS(seed_type, "fixed")) {
+    noise_data.default_engine = std::default_random_engine(fixed_seed_value);
+    noise_data.mt_engine = std::mt19937(fixed_seed_value);
+  } else {
+    CCTK_VERROR("Unrecognized seed type \"%s\"", seed_type);
+  }
+
+  noise_data.distrib = std::uniform_real_distribution<CCTK_REAL>(
+      -noise_amplitude, noise_amplitude);
+
+  return 0;
+}
+
+void add(const cGH *restrict const cctkGH, Loop::GF3D2<CCTK_REAL> var) {
+  DECLARE_CCTK_ARGUMENTS;
+  DECLARE_CCTK_PARAMETERS;
+
+  const Loop::GridDescBase grid(cctkGH);
+
+  if (CCTK_EQUALS(random_engine, "default")) {
+    grid.loop_int<0, 0, 0>(
+        grid.nghostzones,
+        [&] CCTK_HOST(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+          var(p.I) += noise_data.distrib(noise_data.default_engine);
+        });
+
+  } else if (CCTK_EQUALS(random_engine, "mersenne twister")) {
+    grid.loop_int<0, 0, 0>(
+        grid.nghostzones,
+        [&] CCTK_HOST(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+          var(p.I) += noise_data.distrib(noise_data.mt_engine);
+        });
+
+  } else if (CCTK_EQUALS(random_engine, "sine wave")) {
+    grid.loop_int<0, 0, 0>(
+        grid.nghostzones,
+        [=] CCTK_HOST(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+          using std::sin;
+          var(p.I) += noise_amplitude * sin(noise_frequency * p.x / p.dx) *
+                      sin(noise_frequency * p.y / p.dy) *
+                      sin(noise_frequency * p.z / p.dz);
+        });
+
+  } else {
+    CCTK_VERROR("Unrecognized random engine \"%s\"", random_engine);
+  }
+}
+
+} // namespace NoiseX

NoiseX/src/noisex.hxx

@@ -0,0 +1,9 @@
+#include <cctk.h>
+
+#include <loop.hxx>
+
+namespace NoiseX {
+
+void add(const cGH *restrict const cctkGH, Loop::GF3D2<CCTK_REAL> var);
+
+} // namespace NoiseX

scripts/carpetx.th

@@ -709,6 +709,7 @@ CarpetX/FluxWaveToyX
 CarpetX/HydroBaseX
 CarpetX/Loop
 CarpetX/MovingBoxToy
+CarpetX/NoiseX
 CarpetX/ODESolvers
 CarpetX/PDESolvers
 CarpetX/PoissonX