TestRKAB: fixed noise initial data to not repeat itself while using m…

…ultiple threads
EinsteinToolkit · Oct 30, 2024 · 182cec9 · 182cec9
1 parent e05f124
commit 182cec9
Show file tree

Hide file tree

Showing 2 changed files with 20 additions and 14 deletions.
diff --git a/TestRKAB/param.ccl b/TestRKAB/param.ccl
@@ -36,10 +36,10 @@ CCTK_REAL gaussian_width "width of Gaussian"
   (0:* :: ""
 } 1.0
 
-CCTK_REAL noise_seed "Seed of the RNG for noise initial data"
+CCTK_INT noise_seed "Seed of the RNG for noise initial data"
 {
   *:* :: ""
-} 100.0
+} 100
 
 CCTK_REAL noise_boundary "Nois will be generated in the [-noise_boundary, noise_boundary] range"
 {

diff --git a/TestRKAB/src/initial.cxx b/TestRKAB/src/initial.cxx
@@ -48,19 +48,25 @@ extern "C" void TestRKAB_Initial(CCTK_ARGUMENTS) {
         });
 
   } 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};
+    // Make sure our RNG structures are initialized only once, regardless of the
+    // number of threads used.
+    static std::mt19937_64 noise_engine{noise_seed};
+    static std::uniform_real_distribution<CCTK_REAL> noise_distrib{
+        -noise_boundary, noise_boundary};
 
-    grid.loop_int<0, 0, 0>(grid.nghostzones,
-                           [&] CCTK_HOST(const Loop::PointDesc &p)
-                               CCTK_ATTRIBUTE_ALWAYS_INLINE {
-                                 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);
-                               });
+    grid.loop_int<0, 0, 0>(
+        grid.nghostzones,
+        [&] CCTK_HOST(const Loop::PointDesc &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+#pragma omp critical
+          {
+            const auto value{noise_distrib(noise_engine)};
+            phi(p.I) = value;
+            Pi(p.I) = value;
+            Dx(p.I) = value;
+            Dy(p.I) = value;
+            Dz(p.I) = value;
+          }
+        });
   } else {
     CCTK_VERROR("Unknown initial condition \"%s\"", initial_condition);
   }