support a MUSIC2 compatible RNG

cosmo-sims · Oct 25, 2024 · 847a8f0 · 847a8f0
1 parent eb5667c
commit 847a8f0
Show file tree

Hide file tree

Showing 4 changed files with 37 additions and 15 deletions.
diff --git a/example.conf b/example.conf
@@ -107,7 +107,8 @@ seed           = 12345
 
 ##> The MUSIC1 multi-scale random number generator is provided for convenience
 ## warning: MUSIC1 generator is not MPI parallel (yet) (memory is needed for full field on each task)
-# generator      = MUSIC1  
+# generator      = MUSIC1
+# music2_rng     = false
 # seed[7]        = 12345
 # seed[8]        = 23456
 # seed[9]        = 34567

diff --git a/src/plugins/random_music.cc b/src/plugins/random_music.cc
@@ -216,6 +216,11 @@ void RNG_music::parse_random_parameters(void)
 void RNG_music::compute_random_numbers(void)
 {
   bool rndsign = pcf_->get_value_safe<bool>("random", "grafic_sign", false);
+  bool music2_rng = pcf_->get_value_safe<bool>("random", "music2_rng", false);
+  if (music2_rng)
+  {
+    music::ilog.Print("MUSIC1 RNG plugin: using MUSIC2 compatible seeds");
+  }
 
   //--- FILL ALL WHITE NOISE ARRAYS FOR WHICH WE NEED THE FULL FIELD ---//
 
@@ -235,9 +240,9 @@ void RNG_music::compute_random_numbers(void)
       //#warning add possibility to read noise from file also here!
 
       if (rngfnames_[i].size() > 0)
-        music::ilog.Print("Warning: Cannot use filenames for higher levels currently! Ignoring!");
+        music::wlog.Print("Warning: Cannot use filenames for higher levels currently! Ignoring!");
 
-      randc_[i] = new rng(*randc_[i - 1], ran_cube_size_, rngseeds_[i], true);
+      randc_[i] = new rng(*randc_[i - 1], ran_cube_size_, rngseeds_[i], music2_rng, true);
       delete randc_[i - 1];
       randc_[i - 1] = NULL;
     }

diff --git a/src/plugins/random_music_wnoise_generator.cc b/src/plugins/random_music_wnoise_generator.cc
@@ -22,6 +22,21 @@
 #include <random_plugin.hh>
 #include "random_music_wnoise_generator.hh"
 
+inline double Meyer_scaling_function( double k, double kmax )
+{
+  constexpr double twopithirds{2.0*M_PI/3.0};
+  constexpr double fourpithirds{4.0*M_PI/3.0};
+  auto nu = []( double x ){ return x<0.0?0.0:(x<1.0?x:1.0); };
+
+  k = std::abs(k)/kmax * 2 * M_PI;
+
+  if( k < twopithirds ) return 1.0;
+  else if( k< fourpithirds ){
+    return std::cos( 0.5*M_PI * nu(3*k/(2*M_PI)-1.0) );
+  }
+  return 0.0;
+}
+
 template <typename T>
 music_wnoise_generator<T>::music_wnoise_generator(unsigned res, unsigned cubesize, long baseseed, int *x0, int *lx)
     : res_(res), cubesize_(cubesize), ncubes_(1), baseseed_(baseseed)
@@ -43,7 +58,7 @@ music_wnoise_generator<T>::music_wnoise_generator(unsigned res, unsigned cubesiz
 
   initialize();
   mean = fill_all();
-  
+
   if (zeromean)
   {
     mean = 0.0;
@@ -389,7 +404,7 @@ music_wnoise_generator<T>::music_wnoise_generator(/*const*/ music_wnoise_generat
 
 template <typename T>
 music_wnoise_generator<T>::music_wnoise_generator(music_wnoise_generator<T> &rc, unsigned cubesize, long baseseed,
-                                                  bool kspace, bool isolated, int *x0_, int *lx_, bool zeromean)
+                                                  bool music2_rng, bool kspace, bool isolated, int *x0_, int *lx_, bool zeromean)
     : res_(2 * rc.res_), cubesize_(cubesize), ncubes_(1), baseseed_(baseseed)
 {
   initialize();
@@ -496,15 +511,16 @@ music_wnoise_generator<T>::music_wnoise_generator(music_wnoise_generator<T> &rc,
 
           val *= val_phas * sqrt8;
 
-          if (i != (int)nxc / 2 && j != (int)nyc / 2 && k != (int)nzc / 2)
-          {
-            cfine[qf][0] = val.real();
-            cfine[qf][1] = val.imag();
-          }
-          else
-          {
-            //RE(cfine[qf]) = val.real();
-            //IM(cfine[qf]) = 0.0;
+          if(i != (int)nxc / 2 && j != (int)nyc / 2 && k != (int)nzc / 2){
+            double blend_coarse_x = Meyer_scaling_function(kx, nxc / 2);
+            double blend_coarse_y = Meyer_scaling_function(ky, nyc / 2);
+            double blend_coarse_z = Meyer_scaling_function(kz, nzc / 2);
+
+            double blend_coarse = blend_coarse_x*blend_coarse_y*blend_coarse_z;
+            double blend_fine = std::sqrt(1.0-blend_coarse*blend_coarse);
+
+            cfine[qf][0] = blend_fine * cfine[qf][0] + blend_coarse * val.real();
+            cfine[qf][1] = blend_fine * cfine[qf][1] + blend_coarse * val.imag();
           }
         }
 

diff --git a/src/plugins/random_music_wnoise_generator.hh b/src/plugins/random_music_wnoise_generator.hh
@@ -162,7 +162,7 @@ public:
   music_wnoise_generator(unsigned res, unsigned cubesize, long baseseed, int *x0, int *lx);
 
   //! constructor for constrained fine field
-  music_wnoise_generator(music_wnoise_generator<T> &rc, unsigned cubesize, long baseseed,
+  music_wnoise_generator(music_wnoise_generator<T> &rc, unsigned cubesize, long baseseed, bool music2_rng,
                          bool kspace = false, bool isolated = false, int *x0_ = NULL, int *lx_ = NULL, bool zeromean = true);
 
   //! constructor