diff --git a/ADMBaseX/configuration.ccl b/ADMBaseX/configuration.ccl index 198e70337..caaa54e65 100644 --- a/ADMBaseX/configuration.ccl +++ b/ADMBaseX/configuration.ccl @@ -1,3 +1,3 @@ # Configuration definition for thorn ADMBaseX -REQUIRES Loop +REQUIRES Loop NoiseX diff --git a/ADMBaseX/interface.ccl b/ADMBaseX/interface.ccl index adaf3b270..ae7a6438c 100644 --- a/ADMBaseX/interface.ccl +++ b/ADMBaseX/interface.ccl @@ -3,6 +3,7 @@ IMPLEMENTS: ADMBaseX USES INCLUDE HEADER: loop_device.hxx +USES INCLUDE HEADER: noisex.hxx diff --git a/ADMBaseX/param.ccl b/ADMBaseX/param.ccl index 725b64ba3..0d7a267f6 100644 --- a/ADMBaseX/param.ccl +++ b/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 diff --git a/ADMBaseX/schedule.ccl b/ADMBaseX/schedule.ccl index 9d0db9191..a00927403 100644 --- a/ADMBaseX/schedule.ccl +++ b/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 { diff --git a/ADMBaseX/src/noise.cxx b/ADMBaseX/src/noise.cxx index e805b13d7..6f60c4e5e 100644 --- a/ADMBaseX/src/noise.cxx +++ b/ADMBaseX/src/noise.cxx @@ -4,85 +4,40 @@ #include #include -#include +#include 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 distribution(-noise_amplitude, - noise_amplitude); - const auto add_noise = [&](CCTK_REAL &restrict var) { - var += distribution(engine); - }; - - const GF3D gxx_(cctkGH, gxx); - const GF3D gxy_(cctkGH, gxy); - const GF3D gxz_(cctkGH, gxz); - const GF3D gyy_(cctkGH, gyy); - const GF3D gyz_(cctkGH, gyz); - const GF3D gzz_(cctkGH, gzz); - - const GF3D kxx_(cctkGH, kxx); - const GF3D kxy_(cctkGH, kxy); - const GF3D kxz_(cctkGH, kxz); - const GF3D kyy_(cctkGH, kyy); - const GF3D kyz_(cctkGH, kyz); - const GF3D kzz_(cctkGH, kzz); - - const GF3D alp_(cctkGH, alp); - - const GF3D dtalp_(cctkGH, dtalp); - - const GF3D betax_(cctkGH, betax); - const GF3D betay_(cctkGH, betay); - const GF3D betaz_(cctkGH, betaz); - - const GF3D dtbetax_(cctkGH, dtbetax); - const GF3D dtbetay_(cctkGH, dtbetay); - const GF3D 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 diff --git a/NoiseX/README b/NoiseX/README new file mode 100644 index 000000000..721e66bfc --- /dev/null +++ b/NoiseX/README @@ -0,0 +1,11 @@ +Cactus Code Thorn Noise +Author(s) : Lucas Timotheo Sanches +Maintainer(s): Lucas Timotheo Sanches +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. diff --git a/NoiseX/configuration.ccl b/NoiseX/configuration.ccl new file mode 100644 index 000000000..69368e8bd --- /dev/null +++ b/NoiseX/configuration.ccl @@ -0,0 +1,7 @@ +# Configuration definition for thorn NoiseX + +REQUIRES Loop + +PROVIDES NoiseX +{ +} \ No newline at end of file diff --git a/NoiseX/interface.ccl b/NoiseX/interface.ccl new file mode 100644 index 000000000..8107f55b4 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/NoiseX/param.ccl b/NoiseX/param.ccl new file mode 100644 index 000000000..055f664f8 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/NoiseX/schedule.ccl b/NoiseX/schedule.ccl new file mode 100644 index 000000000..4cd9abb97 --- /dev/null +++ b/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" \ No newline at end of file diff --git a/NoiseX/src/make.code.defn b/NoiseX/src/make.code.defn new file mode 100644 index 000000000..fd845890d --- /dev/null +++ b/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 = \ No newline at end of file diff --git a/NoiseX/src/noisex.cxx b/NoiseX/src/noisex.cxx new file mode 100644 index 000000000..ae68aba97 --- /dev/null +++ b/NoiseX/src/noisex.cxx @@ -0,0 +1,78 @@ +#include + +#include +#include +#include + +#include "noisex.hxx" + +#include +#include + +namespace NoiseX { + +struct NoiseData { + std::default_random_engine default_engine{}; + std::mt19937 mt_engine{}; + std::uniform_real_distribution 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( + -noise_amplitude, noise_amplitude); + + return 0; +} + +void add(const cGH *restrict const cctkGH, Loop::GF3D2 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 \ No newline at end of file diff --git a/NoiseX/src/noisex.hxx b/NoiseX/src/noisex.hxx new file mode 100644 index 000000000..36235cb13 --- /dev/null +++ b/NoiseX/src/noisex.hxx @@ -0,0 +1,9 @@ +#include + +#include + +namespace NoiseX { + +void add(const cGH *restrict const cctkGH, Loop::GF3D2 var); + +} // namespace NoiseX \ No newline at end of file diff --git a/scripts/carpetx.th b/scripts/carpetx.th index 7b36eb5f4..9cedc4317 100644 --- a/scripts/carpetx.th +++ b/scripts/carpetx.th @@ -709,6 +709,7 @@ CarpetX/FluxWaveToyX CarpetX/HydroBaseX CarpetX/Loop CarpetX/MovingBoxToy +CarpetX/NoiseX CarpetX/ODESolvers CarpetX/PDESolvers CarpetX/PoissonX