mirror/monofonIC
1
0
Fork 0
mirror of https://github.com/cosmo-sims/monofonIC.git synced 2024-09-19 17:03:45 +02:00
Code Releases Activity

added N-GenIC compatible RNG

Browse source
This commit is contained in:
Oliver Hahn 2019-05-19 17:19:41 +02:00
parent e683fdecd0
commit 91985eefce
5 changed files with 174 additions and 2 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
example.conf
View file

@ -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

3
include/random_plugin.hh
View file

@ -2,6 +2,8 @@
#include <map>
#include <config_file.hh>
#include <general.hh>
#include <grid_fft.hh>
#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<real_t>& g ) const = 0;
//virtual void FillGrid(int level, DensityGrid<real_t> &R) = 0;
};

5
src/main.cc
View file

@ -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<real_t> Conv({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
phi.FillRandomReal(6519);
//phi.FillRandomReal(6519);
the_random_number_generator->Fill_Grid( phi );
//======================================================================
//... compute 1LPT displacement potential ....

2
src/plugins/random_music.cc
View file

@ -40,6 +40,8 @@ public:
bool isMultiscale() const { return true; }
void Fill_Grid( Grid_FFT<real_t>& g ) const { }
void initialize_for_grid_structure()//const refinement_hierarchy &refh)
{
//prefh_ = &refh;

162
src/plugins/random_ngenic.cc Normal file
View file

@ -0,0 +1,162 @@
#include <general.hh>
#include <random_plugin.hh>
#include <config_file.hh>
#include <vector>
#include <cmath>
#include <grid_fft.hh>
#include <gsl/gsl_rng.h>
class RNG_ngenic : public RNG_plugin
{
private:
gsl_rng *pRandomGenerator_;
long RandomSeed_;
size_t nres_;
std::vector<unsigned int> SeedTable_;
public:
explicit RNG_ngenic(ConfigFile &cf) : RNG_plugin(cf)
{
RandomSeed_ = cf.GetValue<long>("random", "seed");
nres_ = cf.GetValue<size_t>("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<real_t> &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<RNG_ngenic> creator("NGENIC");
}