From 91985eefce633a52daf8e9a56db73eb40b0148bd Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Sun, 19 May 2019 17:19:41 +0200 Subject: [PATCH] added N-GenIC compatible RNG --- example.conf | 4 + include/random_plugin.hh | 3 + src/main.cc | 5 +- src/plugins/random_music.cc | 2 + src/plugins/random_ngenic.cc | 162 +++++++++++++++++++++++++++++++++++ 5 files changed, 174 insertions(+), 2 deletions(-) create mode 100644 src/plugins/random_ngenic.cc diff --git a/example.conf b/example.conf index 078e066..0ceec1b 100644 --- a/example.conf +++ b/example.conf @@ -10,6 +10,10 @@ fbase_analysis = output format = gadget2 filename = ics_gadget.dat +[random] +generator = NGENIC +seed = 9001 + [cosmology] transfer = CLASS #eisenstein Omega_m = 0.302 diff --git a/include/random_plugin.hh b/include/random_plugin.hh index 4bc8f6d..b2a0212 100644 --- a/include/random_plugin.hh +++ b/include/random_plugin.hh @@ -2,6 +2,8 @@ #include #include +#include +#include #define DEF_RAN_CUBE_SIZE 32 @@ -16,6 +18,7 @@ class RNG_plugin } virtual ~RNG_plugin() {} virtual bool isMultiscale() const = 0; + virtual void Fill_Grid( Grid_FFT& g ) const = 0; //virtual void FillGrid(int level, DensityGrid &R) = 0; }; diff --git a/src/main.cc b/src/main.cc index 551f6c0..afda857 100644 --- a/src/main.cc +++ b/src/main.cc @@ -71,7 +71,7 @@ int main( int argc, char** argv ) { // print_region_generator_plugins(); print_TransferFunction_plugins(); - // print_RNG_plugins(); + print_RNG_plugins(); print_output_plugins(); csoca::elog << "In order to run, you need to specify a parameter file!" << std::endl; @@ -147,7 +147,8 @@ int main( int argc, char** argv ) OrszagConvolver Conv({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); - phi.FillRandomReal(6519); + //phi.FillRandomReal(6519); + the_random_number_generator->Fill_Grid( phi ); //====================================================================== //... compute 1LPT displacement potential .... diff --git a/src/plugins/random_music.cc b/src/plugins/random_music.cc index b0c3fd7..a13726f 100644 --- a/src/plugins/random_music.cc +++ b/src/plugins/random_music.cc @@ -40,6 +40,8 @@ public: bool isMultiscale() const { return true; } + void Fill_Grid( Grid_FFT& g ) const { } + void initialize_for_grid_structure()//const refinement_hierarchy &refh) { //prefh_ = &refh; diff --git a/src/plugins/random_ngenic.cc b/src/plugins/random_ngenic.cc new file mode 100644 index 0000000..1400606 --- /dev/null +++ b/src/plugins/random_ngenic.cc @@ -0,0 +1,162 @@ + +#include +#include +#include + +#include +#include + +#include +#include + +class RNG_ngenic : public RNG_plugin +{ +private: + gsl_rng *pRandomGenerator_; + long RandomSeed_; + size_t nres_; + std::vector SeedTable_; + +public: + explicit RNG_ngenic(ConfigFile &cf) : RNG_plugin(cf) + { + + RandomSeed_ = cf.GetValue("random", "seed"); + nres_ = cf.GetValue("setup", "GridRes"); + pRandomGenerator_ = gsl_rng_alloc(gsl_rng_ranlxd1); + + SeedTable_.assign(nres_ * nres_, 0u); + + for (size_t i = 0; i < nres_ / 2; ++i) + { + for (size_t j = 0; j < i; j++) + SeedTable_[i * nres_ + j] = 0x7fffffff * gsl_rng_uniform(pRandomGenerator_); + + for (size_t j = 0; j < i + 1; j++) + SeedTable_[j * nres_ + i] = 0x7fffffff * gsl_rng_uniform(pRandomGenerator_); + + for (size_t j = 0; j < i; j++) + SeedTable_[(nres_ - 1 - i) * nres_ + j] = 0x7fffffff * gsl_rng_uniform(pRandomGenerator_); + + for (size_t j = 0; j < i + 1; j++) + SeedTable_[(nres_ - 1 - j) * nres_ + i] = 0x7fffffff * gsl_rng_uniform(pRandomGenerator_); + + for (size_t j = 0; j < i; j++) + SeedTable_[i * nres_ + (nres_ - 1 - j)] = 0x7fffffff * gsl_rng_uniform(pRandomGenerator_); + + for (size_t j = 0; j < i + 1; j++) + SeedTable_[j * nres_ + (nres_ - 1 - i)] = 0x7fffffff * gsl_rng_uniform(pRandomGenerator_); + + for (size_t j = 0; j < i; j++) + SeedTable_[(nres_ - 1 - i) * nres_ + (nres_ - 1 - j)] = 0x7fffffff * gsl_rng_uniform(pRandomGenerator_); + + for (size_t j = 0; j < i + 1; j++) + SeedTable_[(nres_ - 1 - j) * nres_ + (nres_ - 1 - i)] = 0x7fffffff * gsl_rng_uniform(pRandomGenerator_); + } + } + + virtual ~RNG_ngenic() + { + gsl_rng_free(pRandomGenerator_); + } + + bool isMultiscale() const { return false; } + + void Fill_Grid(Grid_FFT &g) const + { + double fnorm = std::pow((double)nres_, -1.5); + g.FourierTransformForward(false); + +#ifdef USE_MPI + // transform is transposed! + for (size_t j = 0; j < g.size(0); ++j) + { + for (size_t i = 0; i < g.size(1); ++i) + { +#else + for (size_t i = 0; i < g.size(0); ++i) + { + for (size_t j = 0; j < g.size(1); ++j) + { +#endif + ptrdiff_t ii = (i>0)? g.size(1) - i : 0; + gsl_rng_set( pRandomGenerator_, SeedTable_[i * nres_ + j]); + for (size_t k = 0; k < g.size(2); ++k) + { + double phase = gsl_rng_uniform(pRandomGenerator_) * 2 * M_PI; + double ampl; + do + { + ampl = gsl_rng_uniform(pRandomGenerator_); + } while (ampl == 0); + if (i == nres_ / 2 || j == nres_ / 2 || k == nres_ / 2) + continue; + if (i == 0 && j == 0 && k == 0) + continue; + + real_t rp = -std::sqrt(-std::log(ampl)) * std::cos(phase);// * fnorm; + real_t ip = -std::sqrt(-std::log(ampl)) * std::sin(phase);// * fnorm; + ccomplex_t zrand(rp,ip); + + if (k > 0) + { + g.kelem(i,j,k) = zrand; + // RE(knoise[(i * res + j) * (res / 2 + 1) + k]) = rp; + // IM(knoise[(i * res + j) * (res / 2 + 1) + k]) = ip; + } + else /* k=0 plane needs special treatment */ + { + if (i == 0) + { + if (j >= nres_ / 2) + { + continue; + } + else + { + int jj = (int)nres_ - (int)j; /* note: j!=0 surely holds at this point */ + g.kelem(i,j,k) = zrand; + g.kelem(i,jj,k) = std::conj(zrand); + + // RE(knoise[(i * res + j) * (res / 2 + 1) + k]) = rp; + // IM(knoise[(i * res + j) * (res / 2 + 1) + k]) = ip; + + // RE(knoise[(i * res + jj) * (res / 2 + 1) + k]) = rp; + // IM(knoise[(i * res + jj) * (res / 2 + 1) + k]) = -ip; + } + } + else + { + if (i >= nres_ / 2) + { + continue; + } + else + { + ptrdiff_t ii = (i>0)? nres_ - i : 0; + ptrdiff_t jj = (j>0)? nres_ - j : 0; + + g.kelem(i,j,k) = zrand; + + // RE(knoise[(i * res + j) * (res / 2 + 1) + k]) = rp; + // IM(knoise[(i * res + j) * (res / 2 + 1) + k]) = ip; + + if (ii >= 0 && ii < (int)nres_) + { + // RE(knoise[(ii * res + jj) * (res / 2 + 1) + k]) = rp; + // IM(knoise[(ii * res + jj) * (res / 2 + 1) + k]) = -ip; + g.kelem(ii,jj,k) = std::conj(zrand); + } + } + } + } + } + } + } + } +}; + +namespace +{ +RNG_plugin_creator_concrete creator("NGENIC"); +}