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

WIP panphasia integration. compiles and runs, but output not tested

Browse source
This commit is contained in:
Oliver Hahn 2023-02-12 14:42:51 -08:00
parent 040324f346
commit bd12e40e22
11 changed files with 2288 additions and 2644 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

3
src/convolution_kernel.hh
View file

@ -4,7 +4,7 @@
a code to generate multi-scale initial conditions
for cosmological simulations
Copyright (C) 2010-19 Oliver Hahn
Copyright (C) 2010-23 Oliver Hahn
*/
@ -16,6 +16,7 @@
#include "config_file.hh"
#include "densities.hh"
#include "density_grid.hh"
#include "transfer_function.hh"
#define ACC_RF(i, j, k) (((((size_t)(i) + nx) % nx) * ny + (((size_t)(j) + ny) % ny)) * 2 * (nz / 2 + 1) + (((size_t)(k) + nz) % nz))

5
src/densities.cc
View file

@ -11,6 +11,7 @@
#include <cstring>
#include "densities.hh"
#include "random.hh"
#include "convolution_kernel.hh"
//TODO: this should be a larger number by default, just to maintain consistency with old default
@ -335,7 +336,7 @@ void fft_interpolate(m1 &V, m2 &v, bool from_basegrid = false)
/*******************************************************************************************/
void GenerateDensityUnigrid(config_file &cf, transfer_function *ptf, tf_type type,
refinement_hierarchy &refh, rand_gen &rand, grid_hierarchy &delta, bool smooth, bool shift)
refinement_hierarchy &refh, noise_generator &rand, grid_hierarchy &delta, bool smooth, bool shift)
{
unsigned levelmin, levelmax, levelminPoisson;
@ -416,7 +417,7 @@ void GenerateDensityUnigrid(config_file &cf, transfer_function *ptf, tf_type typ
/*******************************************************************************************/
void GenerateDensityHierarchy(config_file &cf, transfer_function *ptf, tf_type type,
refinement_hierarchy &refh, rand_gen &rand,
refinement_hierarchy &refh, noise_generator &rand,
grid_hierarchy &delta, bool smooth, bool shift)
{
unsigned levelmin, levelmax, levelminPoisson;

284
src/densities.hh
View file

@ -15,298 +15,20 @@
#include "general.hh"
#include "config_file.hh"
// #include "density_grid.hh"
#include "random.hh"
#include "cosmology.hh"
#include "transfer_function.hh"
#include "general.hh"
void GenerateDensityHierarchy(config_file &cf, transfer_function *ptf, tf_type type,
refinement_hierarchy &refh, rand_gen &rand, grid_hierarchy &delta, bool smooth, bool shift);
refinement_hierarchy &refh, noise_generator &rand, grid_hierarchy &delta, bool smooth, bool shift);
void GenerateDensityUnigrid(config_file &cf, transfer_function *ptf, tf_type type,
refinement_hierarchy &refh, rand_gen &rand, grid_hierarchy &delta, bool smooth, bool shift);
refinement_hierarchy &refh, noise_generator &rand, grid_hierarchy &delta, bool smooth, bool shift);
void normalize_density(grid_hierarchy &delta);
/*!
* @class DensityGrid
* @brief provides infrastructure for computing the initial density field
*
* This class provides access and data management member functions that
* are used when computing the initial density field by convolution with
* transfer functions.
*/
template <typename real_t>
class DensityGrid
{
public:
size_t nx_; //!< number of grid cells in x-direction
size_t ny_; //!< number of grid cells in y-direction
size_t nz_; //!< number of grid cells in z-direction
size_t nzp_; //!< number of cells in memory (z-dir), used for Nyquist padding
size_t nv_[3];
int ox_; //!< offset of grid in x-direction
int oy_; //!< offset of grid in y-direction
int oz_; //!< offset of grid in z-direction
size_t ov_[3];
//! the actual data container in the form of a 1D array
std::vector<real_t> data_;
//! constructor
/*! constructs an instance given the dimensions of the density field
* @param nx the number of cells in x
* @param ny the number of cells in y
* @param nz the number of cells in z
*/
DensityGrid(unsigned nx, unsigned ny, unsigned nz)
: nx_(nx), ny_(ny), nz_(nz), nzp_(2 * (nz_ / 2 + 1)), ox_(0), oy_(0), oz_(0)
{
data_.assign((size_t)nx_ * (size_t)ny_ * (size_t)nzp_, 0.0);
nv_[0] = nx_;
nv_[1] = ny_;
nv_[2] = nz_;
ov_[0] = ox_;
ov_[1] = oy_;
ov_[2] = oz_;
}
DensityGrid(unsigned nx, unsigned ny, unsigned nz, int ox, int oy, int oz)
: nx_(nx), ny_(ny), nz_(nz), nzp_(2 * (nz_ / 2 + 1)), ox_(ox), oy_(oy), oz_(oz)
{
data_.assign((size_t)nx_ * (size_t)ny_ * (size_t)nzp_, 0.0);
nv_[0] = nx_;
nv_[1] = ny_;
nv_[2] = nz_;
ov_[0] = ox_;
ov_[1] = oy_;
ov_[2] = oz_;
}
//! copy constructor
explicit DensityGrid(const DensityGrid<real_t> &g)
: nx_(g.nx_), ny_(g.ny_), nz_(g.nz_), nzp_(g.nzp_),
ox_(g.ox_), oy_(g.oy_), oz_(g.oz_)
{
data_ = g.data_;
nv_[0] = nx_;
nv_[1] = ny_;
nv_[2] = nz_;
ov_[0] = ox_;
ov_[1] = oy_;
ov_[2] = oz_;
}
//! destructor
~DensityGrid()
{
}
//! clears the density object
/*! sets all dimensions to zero and frees the memory
*/
void clear(void)
{
nx_ = ny_ = nz_ = nzp_ = 0;
ox_ = oy_ = oz_ = 0;
nv_[0] = nv_[1] = nv_[2] = 0;
ov_[0] = ov_[1] = ov_[2] = 0;
data_.clear();
std::vector<real_t>().swap(data_);
}
//! query the 3D array sizes of the density object
/*! returns the size of the 3D density object along a specified dimension
* @param i the dimension for which size is to be returned
* @returns array size along dimension i
*/
size_t size(int i)
{
return nv_[i];
}
int offset(int i)
{
return ov_[i];
}
//! zeroes the density object
/*! sets all values to 0.0
*/
void zero(void)
{
data_.assign(data_.size(), 0.0);
}
//! assigns the contents of another DensityGrid to this
DensityGrid &operator=(const DensityGrid<real_t> &g)
{
nx_ = g.nx_;
ny_ = g.ny_;
nz_ = g.nz_;
nzp_ = g.nzp_;
ox_ = g.ox_;
oy_ = g.oy_;
oz_ = g.oz_;
data_ = g.data_;
return *this;
}
//! 3D index based data access operator
inline real_t &operator()(size_t i, size_t j, size_t k)
{
return data_[((size_t)i * ny_ + (size_t)j) * nzp_ + (size_t)k];
}
//! 3D index based const data access operator
inline const real_t &operator()(size_t i, size_t j, size_t k) const
{
return data_[((size_t)i * ny_ + (size_t)j) * nzp_ + (size_t)k];
}
//! recover the pointer to the 1D data array
inline real_t *get_data_ptr(void)
{
return &data_[0];
}
//! fills the density field with random number values
/*! given a pointer to a random_numbers object, fills the field with random values
* @param prc pointer to a random_numbers object
* @param variance the variance of the random numbers (the values returned by prc are multiplied by this)
* @param i0 x-offset (shift) in cells of the density field with respect to the random number field
* @param j0 y-offset (shift) in cells of the density field with respect to the random number field
* @param k0 z-offset (shift) in cells of the density field with respect to the random number field
* @param setzero boolean, if true, the global mean will be subtracted
*/
void fill_rand(/*const*/ random_numbers<real_t> *prc, real_t variance, int i0, int j0, int k0, bool setzero = false)
{
long double sum = 0.0;
#pragma omp parallel for reduction(+ \
: sum)
for (int i = 0; i < nx_; ++i)
for (int j = 0; j < ny_; ++j)
for (int k = 0; k < nz_; ++k)
{
(*this)(i, j, k) = (*prc)(i0 + i, j0 + j, k0 + k) * variance;
sum += (*this)(i, j, k);
}
sum /= nx_ * ny_ * nz_;
if (setzero)
{
#pragma omp parallel for
for (int i = 0; i < nx_; ++i)
for (int j = 0; j < ny_; ++j)
for (int k = 0; k < nz_; ++k)
(*this)(i, j, k) -= sum;
}
}
//! copies the data from another field with 3D index-based access operator
template <class array3>
void copy(array3 &v)
{
#pragma omp parallel for
for (int ix = 0; ix < (int)nx_; ++ix)
for (int iy = 0; iy < (int)ny_; ++iy)
for (int iz = 0; iz < (int)nz_; ++iz)
v(ix, iy, iz) = (*this)(ix, iy, iz);
}
//! adds the data from another field with 3D index-based access operator
template <class array3>
void copy_add(array3 &v)
{
#pragma omp parallel for
for (int ix = 0; ix < (int)nx_; ++ix)
for (int iy = 0; iy < (int)ny_; ++iy)
for (int iz = 0; iz < (int)nz_; ++iz)
v(ix, iy, iz) += (*this)(ix, iy, iz);
}
};
template <typename real_t>
class PaddedDensitySubGrid : public DensityGrid<real_t>
{
public:
using DensityGrid<real_t>::nx_;
using DensityGrid<real_t>::ny_;
using DensityGrid<real_t>::nz_;
using DensityGrid<real_t>::ox_;
using DensityGrid<real_t>::oy_;
using DensityGrid<real_t>::oz_;
using DensityGrid<real_t>::data_;
std::array<size_t,3> pad_;
using DensityGrid<real_t>::fill_rand;
using DensityGrid<real_t>::get_data_ptr;
public:
PaddedDensitySubGrid(int ox, int oy, int oz, unsigned nx, unsigned ny, unsigned nz )
: DensityGrid<real_t>(nx*2, ny*2, nz*2, ox, oy, oz),
pad_{{ nx / 2, ny / 2, nz / 2 }}
{ }
PaddedDensitySubGrid(int ox, int oy, int oz, unsigned nx, unsigned ny, unsigned nz, unsigned padx, unsigned pady, unsigned padz )
: DensityGrid<real_t>(nx + 2 * padx, ny + 2 * pady, nz + 2 * padz, ox, oy, oz),
pad_{{ padx, pady, padz }}
{ }
PaddedDensitySubGrid(const PaddedDensitySubGrid<real_t> &o)
: DensityGrid<real_t>(o)
{ }
size_t margin(int i) const
{
return pad_[i];
}
template <class array3>
void copy_unpad(array3 &v)
{
#pragma omp parallel for
for (size_t ix = pad_[0]; ix < nx_-pad_[0]; ++ix){
const size_t ixu = ix - pad_[0];
for (size_t iy = pad_[1], iyu = 0; iy < ny_-pad_[1]; ++iy, ++iyu)
for (size_t iz = pad_[2], izu = 0; iz < nz_-pad_[2]; ++iz, ++izu)
v(ixu, iyu, izu) = (*this)(ix, iy, iz);
}
}
template <class array3>
void copy_add_unpad(array3 &v)
{
#pragma omp parallel for
for (size_t ix = pad_[0]; ix < nx_-pad_[0]; ++ix){
const size_t ixu = ix - pad_[0];
for (size_t iy = pad_[1], iyu = 0; iy < ny_-pad_[1]; ++iy, ++iyu)
for (size_t iz = pad_[2], izu = 0; iz < nz_-pad_[2]; ++iz, ++izu)
v(ixu, iyu, izu) += (*this)(ix, iy, iz);
}
}
template <class array3>
void copy_subtract_unpad(array3 &v)
{
#pragma omp parallel for
for (size_t ix = pad_[0]; ix < nx_-pad_[0]; ++ix){
const size_t ixu = ix - pad_[0];
for (size_t iy = pad_[1], iyu = 0; iy < ny_-pad_[1]; ++iy, ++iyu)
for (size_t iz = pad_[2], izu = 0; iz < nz_-pad_[2]; ++iz, ++izu)
v(ixu, iyu, izu) -= (*this)(ix, iy, iz);
}
}
};
void coarsen_density(const refinement_hierarchy &rh, GridHierarchy<real_t> &u, bool kspace);
#endif

297
src/density_grid.hh Normal file
View file

@ -0,0 +1,297 @@
#pragma once
#include <vector>
#include <array>
/*!
* @class DensityGrid
* @brief provides infrastructure for computing the initial density field
*
* This class provides access and data management member functions that
* are used when computing the initial density field by convolution with
* transfer functions.
*/
template <typename real_t>
class DensityGrid
{
public:
size_t nx_; //!< number of grid cells in x-direction
size_t ny_; //!< number of grid cells in y-direction
size_t nz_; //!< number of grid cells in z-direction
size_t nzp_; //!< number of cells in memory (z-dir), used for Nyquist padding
std::array<size_t,3> nv_;
int ox_; //!< offset of grid in x-direction
int oy_; //!< offset of grid in y-direction
int oz_; //!< offset of grid in z-direction
std::array<int,3> ov_;
//! the actual data container in the form of a 1D array
std::vector<real_t> data_;
//! constructor
/*! constructs an instance given the dimensions of the density field
* @param nx the number of cells in x
* @param ny the number of cells in y
* @param nz the number of cells in z
*/
DensityGrid(unsigned nx, unsigned ny, unsigned nz)
: nx_(nx), ny_(ny), nz_(nz), nzp_(2 * (nz_ / 2 + 1)), ox_(0), oy_(0), oz_(0)
{
data_.assign((size_t)nx_ * (size_t)ny_ * (size_t)nzp_, 0.0);
nv_ = {nx_,ny_,nz_};
ov_ = {ox_,oy_,oz_};
// nv_[0] = nx_;
// nv_[1] = ny_;
// nv_[2] = nz_;
// ov_[0] = ox_;
// ov_[1] = oy_;
// ov_[2] = oz_;
}
DensityGrid(unsigned nx, unsigned ny, unsigned nz, int ox, int oy, int oz)
: nx_(nx), ny_(ny), nz_(nz), nzp_(2 * (nz_ / 2 + 1)), ox_(ox), oy_(oy), oz_(oz)
{
data_.assign((size_t)nx_ * (size_t)ny_ * (size_t)nzp_, 0.0);
nv_ = {nx_,ny_,nz_};
ov_ = {ox_,oy_,oz_};
// nv_[0] = nx_;
// nv_[1] = ny_;
// nv_[2] = nz_;
// ov_[0] = ox_;
// ov_[1] = oy_;
// ov_[2] = oz_;
}
//! copy constructor
explicit DensityGrid(const DensityGrid<real_t> &g)
: nx_(g.nx_), ny_(g.ny_), nz_(g.nz_), nzp_(g.nzp_),
ox_(g.ox_), oy_(g.oy_), oz_(g.oz_)
{
data_ = g.data_;
nv_ = {nx_,ny_,nz_};
ov_ = {ox_,oy_,oz_};
// nv_[0] = nx_;
// nv_[1] = ny_;
// nv_[2] = nz_;
// ov_[0] = ox_;
// ov_[1] = oy_;
// ov_[2] = oz_;
}
//! destructor
~DensityGrid()
{
}
//! clears the density object
/*! sets all dimensions to zero and frees the memory
*/
void clear(void)
{
nx_ = ny_ = nz_ = nzp_ = 0;
ox_ = oy_ = oz_ = 0;
nv_[0] = nv_[1] = nv_[2] = 0;
ov_[0] = ov_[1] = ov_[2] = 0;
data_.clear();
std::vector<real_t>().swap(data_);
}
//! query the 3D array sizes of the density object
/*! returns the size of the 3D density object along a specified dimension
* @param i the dimension for which size is to be returned
* @returns array size along dimension i
*/
size_t size(int i)
{
return nv_[i];
}
int offset(int i)
{
return ov_[i];
}
//! zeroes the density object
/*! sets all values to 0.0
*/
void zero(void)
{
data_.assign(data_.size(), 0.0);
}
//! assigns the contents of another DensityGrid to this
DensityGrid &operator=(const DensityGrid<real_t> &g)
{
nx_ = g.nx_;
ny_ = g.ny_;
nz_ = g.nz_;
nzp_ = g.nzp_;
ox_ = g.ox_;
oy_ = g.oy_;
oz_ = g.oz_;
data_ = g.data_;
return *this;
}
//! 3D index based data access operator
inline real_t &operator()(size_t i, size_t j, size_t k)
{
return data_[((size_t)i * ny_ + (size_t)j) * nzp_ + (size_t)k];
}
//! 3D index based const data access operator
inline const real_t &operator()(size_t i, size_t j, size_t k) const
{
return data_[((size_t)i * ny_ + (size_t)j) * nzp_ + (size_t)k];
}
//! recover the pointer to the 1D data array
inline real_t *get_data_ptr(void)
{
return &data_[0];
}
#if 0
//! fills the density field with random number values
/*! given a pointer to a random_numbers object, fills the field with random values
* @param prc pointer to a random_numbers object
* @param variance the variance of the random numbers (the values returned by prc are multiplied by this)
* @param i0 x-offset (shift) in cells of the density field with respect to the random number field
* @param j0 y-offset (shift) in cells of the density field with respect to the random number field
* @param k0 z-offset (shift) in cells of the density field with respect to the random number field
* @param setzero boolean, if true, the global mean will be subtracted
*/
void fill_rand(/*const*/ random_numbers<real_t> *prc, real_t variance, int i0, int j0, int k0, bool setzero = false)
{
long double sum = 0.0;
#pragma omp parallel for reduction(+ \
: sum)
for (int i = 0; i < nx_; ++i)
for (int j = 0; j < ny_; ++j)
for (int k = 0; k < nz_; ++k)
{
(*this)(i, j, k) = (*prc)(i0 + i, j0 + j, k0 + k) * variance;
sum += (*this)(i, j, k);
}
sum /= nx_ * ny_ * nz_;
if (setzero)
{
#pragma omp parallel for
for (int i = 0; i < nx_; ++i)
for (int j = 0; j < ny_; ++j)
for (int k = 0; k < nz_; ++k)
(*this)(i, j, k) -= sum;
}
}
#endif
//! copies the data from another field with 3D index-based access operator
template <class array3>
void copy(array3 &v)
{
#pragma omp parallel for
for (int ix = 0; ix < (int)nx_; ++ix)
for (int iy = 0; iy < (int)ny_; ++iy)
for (int iz = 0; iz < (int)nz_; ++iz)
v(ix, iy, iz) = (*this)(ix, iy, iz);
}
//! adds the data from another field with 3D index-based access operator
template <class array3>
void copy_add(array3 &v)
{
#pragma omp parallel for
for (int ix = 0; ix < (int)nx_; ++ix)
for (int iy = 0; iy < (int)ny_; ++iy)
for (int iz = 0; iz < (int)nz_; ++iz)
v(ix, iy, iz) += (*this)(ix, iy, iz);
}
};
template <typename real_t>
class PaddedDensitySubGrid : public DensityGrid<real_t>
{
public:
using DensityGrid<real_t>::nx_;
using DensityGrid<real_t>::ny_;
using DensityGrid<real_t>::nz_;
using DensityGrid<real_t>::ox_;
using DensityGrid<real_t>::oy_;
using DensityGrid<real_t>::oz_;
using DensityGrid<real_t>::data_;
std::array<size_t, 3> pad_;
// using DensityGrid<real_t>::fill_rand;
using DensityGrid<real_t>::get_data_ptr;
public:
PaddedDensitySubGrid(int ox, int oy, int oz, unsigned nx, unsigned ny, unsigned nz)
: DensityGrid<real_t>(nx * 2, ny * 2, nz * 2, ox, oy, oz),
pad_{{nx / 2, ny / 2, nz / 2}}
{
}
PaddedDensitySubGrid(int ox, int oy, int oz, unsigned nx, unsigned ny, unsigned nz, unsigned padx, unsigned pady, unsigned padz)
: DensityGrid<real_t>(nx + 2 * padx, ny + 2 * pady, nz + 2 * padz, ox, oy, oz),
pad_{{padx, pady, padz}}
{
}
PaddedDensitySubGrid(const PaddedDensitySubGrid<real_t> &o)
: DensityGrid<real_t>(o)
{
}
size_t margin(int i) const
{
return pad_[i];
}
template <class array3>
void copy_unpad(array3 &v)
{
#pragma omp parallel for
for (size_t ix = pad_[0]; ix < nx_ - pad_[0]; ++ix)
{
const size_t ixu = ix - pad_[0];
for (size_t iy = pad_[1], iyu = 0; iy < ny_ - pad_[1]; ++iy, ++iyu)
for (size_t iz = pad_[2], izu = 0; iz < nz_ - pad_[2]; ++iz, ++izu)
v(ixu, iyu, izu) = (*this)(ix, iy, iz);
}
}
template <class array3>
void copy_add_unpad(array3 &v)
{
#pragma omp parallel for
for (size_t ix = pad_[0]; ix < nx_ - pad_[0]; ++ix)
{
const size_t ixu = ix - pad_[0];
for (size_t iy = pad_[1], iyu = 0; iy < ny_ - pad_[1]; ++iy, ++iyu)
for (size_t iz = pad_[2], izu = 0; iz < nz_ - pad_[2]; ++iz, ++izu)
v(ixu, iyu, izu) += (*this)(ix, iy, iz);
}
}
template <class array3>
void copy_subtract_unpad(array3 &v)
{
#pragma omp parallel for
for (size_t ix = pad_[0]; ix < nx_ - pad_[0]; ++ix)
{
const size_t ixu = ix - pad_[0];
for (size_t iy = pad_[1], iyu = 0; iy < ny_ - pad_[1]; ++iy, ++iyu)
for (size_t iz = pad_[2], izu = 0; iz < nz_ - pad_[2]; ++iz, ++izu)
v(ixu, iyu, izu) -= (*this)(ix, iy, iz);
}
}
};

35
src/main.cc
View file

@ -43,7 +43,7 @@ extern "C"
#include "transfer_function.hh"
#define THE_CODE_NAME "music!"
#define THE_CODE_VERSION "1.53"
#define THE_CODE_VERSION "2.0a"
namespace music
{
@ -249,7 +249,7 @@ double compute_finest_sigma(grid_hierarchy &u)
return sqrt(sum2 - sum * sum);
}
double compute_finest_max(grid_hierarchy &u)
double compute_finest_absmax(grid_hierarchy &u)
{
double valmax = 0.0;
#pragma omp parallel for reduction(max:valmax)
@ -257,8 +257,8 @@ double compute_finest_max(grid_hierarchy &u)
for (int iy = 0; iy < (int)(*u.get_grid(u.levelmax())).size(1); ++iy)
for (int iz = 0; iz < (int)(*u.get_grid(u.levelmax())).size(2); ++iz)
{
if (fabs((*u.get_grid(u.levelmax()))(ix, iy, iz)) > fabs(valmax))
valmax = (*u.get_grid(u.levelmax()))(ix, iy, iz);
if (std::fabs((*u.get_grid(u.levelmax()))(ix, iy, iz)) > valmax)
valmax = std::fabs((*u.get_grid(u.levelmax()))(ix, iy, iz));
}
return valmax;
@ -299,7 +299,6 @@ void add_constant_value( grid_hierarchy &u, const double val )
/*****************************************************************************************************/
region_generator_plugin *the_region_generator;
RNG_plugin *the_random_number_generator;
int main(int argc, const char *argv[])
{
@ -461,8 +460,7 @@ int main(int argc, const char *argv[])
}
the_region_generator = select_region_generator_plugin(cf);
the_random_number_generator = select_RNG_plugin(cf);
//------------------------------------------------------------------------------
//... determine run parameters
//------------------------------------------------------------------------------
@ -472,6 +470,12 @@ int main(int argc, const char *argv[])
<< " distinct amplitudes for baryon and DM fields!\n"
<< " Perturbation amplitudes will be identical!" << std::endl;
//------------------------------------------------------------------------------
//... start up the random number generator plugin
//... see if we need to set some grid building constraints
noise_generator rand( cf, the_transfer_function_plugin );
//------------------------------------------------------------------------------
//... determine the refinement hierarchy
//------------------------------------------------------------------------------
@ -505,7 +509,8 @@ int main(int argc, const char *argv[])
std::cout << " GENERATING WHITE NOISE\n";
std::cout << "-------------------------------------------------------------\n";
LOGUSER("Computing white noise...");
rand_gen rand(cf, rh_TF, the_transfer_function_plugin);
// rand_gen rand(cf, rh_TF, the_transfer_function_plugin);
rand.initialize_for_grid_structure( rh_TF );
//------------------------------------------------------------------------------
//... initialize the Poisson solver
@ -608,7 +613,7 @@ int main(int argc, const char *argv[])
else
//... displacement
the_poisson_solver->gradient(icoord, u, data_forIO);
double dispmax = compute_finest_max(data_forIO);
double dispmax = compute_finest_absmax(data_forIO);
LOGINFO("max. %c-displacement of HR particles is %f [mean dx]", 'x' + icoord, dispmax * (double)(1ll << data_forIO.levelmax()));
coarsen_density(rh_Poisson, data_forIO, false);
@ -745,7 +750,7 @@ int main(int argc, const char *argv[])
LOGINFO("mean of %c-velocity of high-res particles is %f", 'x' + icoord, meanv);
LOGUSER("mean of %c-velocity of high-res particles is %f", 'x' + icoord, meanv);
double maxv = compute_finest_max(data_forIO);
double maxv = compute_finest_absmax(data_forIO);
LOGINFO("max of abs of %c-velocity of high-res particles is %f", 'x' + icoord, maxv);
coarsen_density(rh_Poisson, data_forIO, false);
@ -815,7 +820,7 @@ int main(int argc, const char *argv[])
LOGINFO("mean of %c-velocity of high-res particles is %f", 'x' + icoord, meanv);
LOGUSER("mean of %c-velocity of high-res particles is %f", 'x' + icoord, meanv);
double maxv = compute_finest_max(data_forIO);
double maxv = compute_finest_absmax(data_forIO);
LOGINFO("max of abs of %c-velocity of high-res particles is %f", 'x' + icoord, maxv);
coarsen_density(rh_Poisson, data_forIO, false);
@ -874,7 +879,7 @@ int main(int argc, const char *argv[])
LOGINFO("mean of %c-velocity of high-res baryons is %f", 'x' + icoord, meanv);
LOGUSER("mean of %c-velocity of high-res baryons is %f", 'x' + icoord, meanv);
double maxv = compute_finest_max(data_forIO);
double maxv = compute_finest_absmax(data_forIO);
LOGINFO("max of abs of %c-velocity of high-res baryons is %f", 'x' + icoord, maxv);
coarsen_density(rh_Poisson, data_forIO, false);
@ -993,7 +998,7 @@ int main(int argc, const char *argv[])
LOGINFO("mean of %c-velocity of high-res particles is %f", 'x' + icoord, meanv);
LOGUSER("mean of %c-velocity of high-res particles is %f", 'x' + icoord, meanv);
double maxv = compute_finest_max(data_forIO);
double maxv = compute_finest_absmax(data_forIO);
LOGINFO("max of abs of %c-velocity of high-res particles is %f", 'x' + icoord, maxv);
std::cerr << " - velocity component " << icoord << " : sigma = " << sigv << std::endl;
@ -1091,7 +1096,7 @@ int main(int argc, const char *argv[])
LOGINFO("mean of %c-velocity of high-res baryons is %f", 'x' + icoord, meanv);
LOGUSER("mean of %c-velocity of high-res baryons is %f", 'x' + icoord, meanv);
double maxv = compute_finest_max(data_forIO);
double maxv = compute_finest_absmax(data_forIO);
LOGINFO("max of abs of %c-velocity of high-res baryons is %f", 'x' + icoord, maxv);
std::cerr << " - velocity component " << icoord << " : sigma = " << sigv << std::endl;
@ -1199,7 +1204,7 @@ int main(int argc, const char *argv[])
else
the_poisson_solver->gradient(icoord, u1, data_forIO);
double dispmax = compute_finest_max(data_forIO);
double dispmax = compute_finest_absmax(data_forIO);
LOGINFO("max. %c-displacement of HR particles is %f [mean dx]", 'x' + icoord, dispmax * (double)(1ll << data_forIO.levelmax()));
coarsen_density(rh_Poisson, data_forIO, false);

613
src/plugins/random_music.cc
View file

@ -1,20 +1,607 @@
#include "random.hh"
#include "random_music_wnoise_generator.hh"
typedef music_wnoise_generator<real_t> rng;
class RNG_music : public RNG_plugin
{
protected:
std::vector<long> rngseeds_;
std::vector<std::string> rngfnames_;
unsigned ran_cube_size_;
int levelmin_, levelmax_, levelmin_seed_;
bool disk_cached_;
bool restart_;
bool initialized_;
std::vector<std::vector<real_t> *> mem_cache_;
//! checks if the specified string is numeric
bool is_number(const std::string &s);
//! parses the random number parameters in the conf file
void parse_random_parameters(void);
//! computes the white noise fields and keeps them either in memory or on disk
void compute_random_numbers(void);
//! adjusts averages
void correct_avg(int icoarse, int ifine);
//! store the white noise fields in memory or on disk
void store_rnd(int ilevel, rng *prng);
class RNG_music : public RNG_plugin{
public:
explicit RNG_music( config_file& cf )
: RNG_plugin( cf )
{ }
~RNG_music() { }
bool is_multiscale() const
{ return true; }
explicit RNG_music(config_file &cf) : RNG_plugin(cf), initialized_(false) {}
~RNG_music() {}
bool is_multiscale() const { return true; }
void initialize_for_grid_structure(const refinement_hierarchy &refh)
{
prefh_ = &refh;
levelmin_ = prefh_->levelmin();
levelmax_ = prefh_->levelmax();
ran_cube_size_ = pcf_->getValueSafe<unsigned>("random", "cubesize", DEF_RAN_CUBE_SIZE);
disk_cached_ = pcf_->getValueSafe<bool>("random", "disk_cached", true);
restart_ = pcf_->getValueSafe<bool>("random", "restart", false);
mem_cache_.assign(levelmax_ - levelmin_ + 1, (std::vector<real_t> *)NULL);
if (restart_ && !disk_cached_)
{
LOGERR("Cannot restart from mem cached random numbers.");
throw std::runtime_error("Cannot restart from mem cached random numbers.");
}
//... determine seed/white noise file data to be applied
parse_random_parameters();
if (!restart_)
{
//... compute the actual random numbers
compute_random_numbers();
}
initialized_ = true;
}
void fill_grid(int level, DensityGrid<real_t> &R);
};
bool RNG_music::is_number(const std::string &s)
{
for (size_t i = 0; i < s.length(); i++)
if (!std::isdigit(s[i]) && s[i] != '-')
return false;
namespace{
RNG_plugin_creator_concrete< RNG_music > creator("MUSIC");
return true;
}
void RNG_music::parse_random_parameters(void)
{
//... parse random number options
for (int i = 0; i <= 100; ++i)
{
char seedstr[128];
std::string tempstr;
bool noseed = false;
sprintf(seedstr, "seed[%d]", i);
if (pcf_->containsKey("random", seedstr))
tempstr = pcf_->getValue<std::string>("random", seedstr);
else
{
// "-2" means that no seed entry was found for that level
tempstr = std::string("-2");
noseed = true;
}
if (is_number(tempstr))
{
long ltemp;
pcf_->convert(tempstr, ltemp);
rngfnames_.push_back("");
if (noseed) // ltemp < 0 )
//... generate some dummy seed which only depends on the level, negative so we know it's not
//... an actual seed and thus should not be used as a constraint for coarse levels
// rngseeds_.push_back( -abs((unsigned)(ltemp-i)%123+(unsigned)(ltemp+827342523521*i)%123456789) );
rngseeds_.push_back(-abs((long)(ltemp - i) % 123 + (long)(ltemp + 7342523521 * i) % 123456789));
else
{
if (ltemp <= 0)
{
LOGERR("Specified seed [random]/%s needs to be a number >0!", seedstr);
throw std::runtime_error("Seed values need to be >0");
}
rngseeds_.push_back(ltemp);
}
}
else
{
rngfnames_.push_back(tempstr);
rngseeds_.push_back(-1);
LOGINFO("Random numbers for level %3d will be read from file.", i);
}
}
//.. determine for which levels random seeds/random number files are given
levelmin_seed_ = -1;
for (unsigned ilevel = 0; ilevel < rngseeds_.size(); ++ilevel)
{
if (levelmin_seed_ < 0 && (rngfnames_[ilevel].size() > 0 || rngseeds_[ilevel] >= 0))
levelmin_seed_ = ilevel;
}
}
void RNG_music::compute_random_numbers(void)
{
bool rndsign = pcf_->getValueSafe<bool>("random", "grafic_sign", false);
std::vector<rng *> randc(std::max(levelmax_, levelmin_seed_) + 1, (rng *)NULL);
//--- FILL ALL WHITE NOISE ARRAYS FOR WHICH WE NEED THE FULL FIELD ---//
//... seeds are given for a level coarser than levelmin
if (levelmin_seed_ < levelmin_)
{
if (rngfnames_[levelmin_seed_].size() > 0)
randc[levelmin_seed_] = new rng(1 << levelmin_seed_, rngfnames_[levelmin_seed_], rndsign);
else
randc[levelmin_seed_] = new rng(1 << levelmin_seed_, ran_cube_size_, rngseeds_[levelmin_seed_]);
for (int i = levelmin_seed_ + 1; i <= levelmin_; ++i)
{
// #warning add possibility to read noise from file also here!
if (rngfnames_[i].size() > 0)
LOGINFO("Warning: Cannot use filenames for higher levels currently! Ignoring!");
randc[i] = new rng(*randc[i - 1], ran_cube_size_, rngseeds_[i]);
delete randc[i - 1];
randc[i - 1] = NULL;
}
}
//... seeds are given for a level finer than levelmin, obtain by averaging
if (levelmin_seed_ > levelmin_)
{
if (rngfnames_[levelmin_seed_].size() > 0)
randc[levelmin_seed_] = new rng(1 << levelmin_seed_, rngfnames_[levelmin_seed_], rndsign);
else
randc[levelmin_seed_] =
new rng(1 << levelmin_seed_, ran_cube_size_, rngseeds_[levelmin_seed_]); //, x0, lx );
for (int ilevel = levelmin_seed_ - 1; ilevel >= (int)levelmin_; --ilevel)
{
if (rngseeds_[ilevel - levelmin_] > 0)
LOGINFO("Warning: random seed for level %d will be ignored.\n"
" consistency requires that it is obtained by restriction from level %d",
ilevel, levelmin_seed_);
randc[ilevel] = new rng(*randc[ilevel + 1]);
if (ilevel + 1 > levelmax_)
{
delete randc[ilevel + 1];
randc[ilevel + 1] = NULL;
}
}
}
//--- GENERATE AND STORE ALL LEVELS, INCLUDING REFINEMENTS ---//
//... levelmin
if (randc[levelmin_] == NULL)
{
if (rngfnames_[levelmin_].size() > 0)
randc[levelmin_] = new rng(1 << levelmin_, rngfnames_[levelmin_], rndsign);
else
randc[levelmin_] = new rng(1 << levelmin_, ran_cube_size_, rngseeds_[levelmin_]);
}
store_rnd(levelmin_, randc[levelmin_]);
//... refinement levels
for (int ilevel = levelmin_ + 1; ilevel <= levelmax_; ++ilevel)
{
int lx[3], x0[3];
int shift[3], levelmin_poisson;
shift[0] = pcf_->getValue<int>("setup", "shift_x");
shift[1] = pcf_->getValue<int>("setup", "shift_y");
shift[2] = pcf_->getValue<int>("setup", "shift_z");
levelmin_poisson = pcf_->getValue<unsigned>("setup", "levelmin");
int lfac = 1 << (ilevel - levelmin_poisson);
std::array<int,3> margin;
if( prefh_->get_margin()>0 ){
margin[0] = prefh_->get_margin();
margin[1] = prefh_->get_margin();
margin[2] = prefh_->get_margin();
}else{
margin[0] = prefh_->size(ilevel, 0)/2;
margin[1] = prefh_->size(ilevel, 1)/2;
margin[2] = prefh_->size(ilevel, 2)/2;
}
lx[0] = prefh_->size(ilevel, 0) + 2 * margin[0];
lx[1] = prefh_->size(ilevel, 1) + 2 * margin[1];
lx[2] = prefh_->size(ilevel, 2) + 2 * margin[2];
x0[0] = prefh_->offset_abs(ilevel, 0) - lfac * shift[0] - margin[0];
x0[1] = prefh_->offset_abs(ilevel, 1) - lfac * shift[1] - margin[1];
x0[2] = prefh_->offset_abs(ilevel, 2) - lfac * shift[2] - margin[2];
if (randc[ilevel] == NULL)
randc[ilevel] =
new rng(*randc[ilevel - 1], ran_cube_size_, rngseeds_[ilevel], x0, lx);
delete randc[ilevel - 1];
randc[ilevel - 1] = NULL;
//... store numbers
store_rnd(ilevel, randc[ilevel]);
}
delete randc[levelmax_];
randc[levelmax_] = NULL;
}
void RNG_music::store_rnd(int ilevel, rng *prng)
{
int shift[3], levelmin_poisson;
shift[0] = pcf_->getValue<int>("setup", "shift_x");
shift[1] = pcf_->getValue<int>("setup", "shift_y");
shift[2] = pcf_->getValue<int>("setup", "shift_z");
levelmin_poisson = pcf_->getValue<unsigned>("setup", "levelmin");
int lfac = 1 << (ilevel - levelmin_poisson);
bool grafic_out = false;
if (grafic_out)
{
std::vector<float> data;
if (ilevel == levelmin_)
{
int N = 1 << levelmin_;
int i0, j0, k0;
i0 = -lfac * shift[0];
j0 = -lfac * shift[1];
k0 = -lfac * shift[2];
char fname[128];
sprintf(fname, "grafic_wnoise_%04d.bin", ilevel);
LOGUSER("Storing white noise field for grafic in file \'%s\'...", fname);
std::ofstream ofs(fname, std::ios::binary | std::ios::trunc);
data.assign(N * N, 0.0);
int blksize = 4 * sizeof(int);
int iseed = 0;
ofs.write(reinterpret_cast<char *>(&blksize), sizeof(int));
ofs.write(reinterpret_cast<char *>(&N), sizeof(int));
ofs.write(reinterpret_cast<char *>(&N), sizeof(int));
ofs.write(reinterpret_cast<char *>(&N), sizeof(int));
ofs.write(reinterpret_cast<char *>(&iseed), sizeof(int));
ofs.write(reinterpret_cast<char *>(&blksize), sizeof(int));
for (int k = 0; k < N; ++k)
{
#pragma omp parallel for
for (int j = 0; j < N; ++j)
for (int i = 0; i < N; ++i)
data[j * N + i] = -(*prng)(i + i0, j + j0, k + k0);
blksize = N * N * sizeof(float);
ofs.write(reinterpret_cast<char *>(&blksize), sizeof(int));
ofs.write(reinterpret_cast<char *>(&data[0]), N * N * sizeof(float));
ofs.write(reinterpret_cast<char *>(&blksize), sizeof(int));
}
ofs.close();
}
else
{
int nx, ny, nz;
int i0, j0, k0;
nx = prefh_->size(ilevel, 0);
ny = prefh_->size(ilevel, 1);
nz = prefh_->size(ilevel, 2);
i0 = prefh_->offset_abs(ilevel, 0) - lfac * shift[0];
j0 = prefh_->offset_abs(ilevel, 1) - lfac * shift[1];
k0 = prefh_->offset_abs(ilevel, 2) - lfac * shift[2];
char fname[128];
sprintf(fname, "grafic_wnoise_%04d.bin", ilevel);
LOGUSER("Storing white noise field for grafic in file \'%s\'...", fname);
LOGDEBUG("(%d,%d,%d) -- (%d,%d,%d) -- lfac = %d", nx, ny, nz, i0, j0, k0, lfac);
std::ofstream ofs(fname, std::ios::binary | std::ios::trunc);
data.assign(nx * ny, 0.0);
int blksize = 4 * sizeof(int);
int iseed = 0;
ofs.write(reinterpret_cast<char *>(&blksize), sizeof(int));
ofs.write(reinterpret_cast<char *>(&nz), sizeof(unsigned));
ofs.write(reinterpret_cast<char *>(&ny), sizeof(unsigned));
ofs.write(reinterpret_cast<char *>(&nx), sizeof(unsigned));
ofs.write(reinterpret_cast<char *>(&iseed), sizeof(int));
ofs.write(reinterpret_cast<char *>(&blksize), sizeof(int));
for (int k = 0; k < nz; ++k)
{
#pragma omp parallel for
for (int j = 0; j < ny; ++j)
for (int i = 0; i < nx; ++i)
data[j * nx + i] = -(*prng)(i + i0, j + j0, k + k0);
blksize = nx * ny * sizeof(float);
ofs.write(reinterpret_cast<char *>(&blksize), sizeof(int));
ofs.write(reinterpret_cast<char *>(&data[0]), nx * ny * sizeof(float));
ofs.write(reinterpret_cast<char *>(&blksize), sizeof(int));
}
ofs.close();
}
}
if (disk_cached_)
{
std::vector<real_t> data;
if (ilevel == levelmin_)
{
int N = 1 << levelmin_;
int i0, j0, k0;
i0 = -lfac * shift[0];
j0 = -lfac * shift[1];
k0 = -lfac * shift[2];
char fname[128];
sprintf(fname, "wnoise_%04d.bin", ilevel);
LOGUSER("Storing white noise field in file \'%s\'...", fname);
std::ofstream ofs(fname, std::ios::binary | std::ios::trunc);
ofs.write(reinterpret_cast<char *>(&N), sizeof(unsigned));
ofs.write(reinterpret_cast<char *>(&N), sizeof(unsigned));
ofs.write(reinterpret_cast<char *>(&N), sizeof(unsigned));
data.assign(N * N, 0.0);
for (int i = 0; i < N; ++i)
{
#pragma omp parallel for
for (int j = 0; j < N; ++j)
for (int k = 0; k < N; ++k)
data[j * N + k] = (*prng)(i + i0, j + j0, k + k0);
ofs.write(reinterpret_cast<char *>(&data[0]), N * N * sizeof(real_t));
}
ofs.close();
}
else
{
int nx, ny, nz;
int i0, j0, k0;
std::array<int,3> margin;
if( prefh_->get_margin()>0 ){
margin[0] = prefh_->get_margin();
margin[1] = prefh_->get_margin();
margin[2] = prefh_->get_margin();
}else{
margin[0] = prefh_->size(ilevel, 0)/2;
margin[1] = prefh_->size(ilevel, 1)/2;
margin[2] = prefh_->size(ilevel, 2)/2;
}
nx = prefh_->size(ilevel, 0) + 2 * margin[0];
ny = prefh_->size(ilevel, 1) + 2 * margin[1];
nz = prefh_->size(ilevel, 2) + 2 * margin[2];
i0 = prefh_->offset_abs(ilevel, 0) - lfac * shift[0] - margin[0];
j0 = prefh_->offset_abs(ilevel, 1) - lfac * shift[1] - margin[1];
k0 = prefh_->offset_abs(ilevel, 2) - lfac * shift[2] - margin[2];
char fname[128];
sprintf(fname, "wnoise_%04d.bin", ilevel);
LOGUSER("Storing white noise field in file \'%s\'...", fname);
std::ofstream ofs(fname, std::ios::binary | std::ios::trunc);
ofs.write(reinterpret_cast<char *>(&nx), sizeof(unsigned));
ofs.write(reinterpret_cast<char *>(&ny), sizeof(unsigned));
ofs.write(reinterpret_cast<char *>(&nz), sizeof(unsigned));
data.assign(ny * nz, 0.0);
for (int i = 0; i < nx; ++i)
{
#pragma omp parallel for
for (int j = 0; j < ny; ++j)
for (int k = 0; k < nz; ++k)
data[j * nz + k] = (*prng)(i + i0, j + j0, k + k0);
ofs.write(reinterpret_cast<char *>(&data[0]), ny * nz * sizeof(real_t));
}
ofs.close();
}
}
else
{
int nx, ny, nz;
int i0, j0, k0;
if (ilevel == levelmin_)
{
i0 = -lfac * shift[0];
j0 = -lfac * shift[1];
k0 = -lfac * shift[2];
nx = ny = nz = 1 << levelmin_;
}
else
{
std::array<int,3> margin;
if( prefh_->get_margin()>0 ){
margin[0] = prefh_->get_margin();
margin[1] = prefh_->get_margin();
margin[2] = prefh_->get_margin();
}else{
margin[0] = prefh_->size(ilevel, 0)/2;
margin[1] = prefh_->size(ilevel, 1)/2;
margin[2] = prefh_->size(ilevel, 2)/2;
}
nx = prefh_->size(ilevel, 0) + 2 * margin[0];
ny = prefh_->size(ilevel, 1) + 2 * margin[1];
nz = prefh_->size(ilevel, 2) + 2 * margin[2];
i0 = prefh_->offset_abs(ilevel, 0) - lfac * shift[0] - margin[0];
j0 = prefh_->offset_abs(ilevel, 1) - lfac * shift[1] - margin[1];
k0 = prefh_->offset_abs(ilevel, 2) - lfac * shift[2] - margin[2];
}
mem_cache_[ilevel - levelmin_] = new std::vector<real_t>(nx * ny * nz, 0.0);
LOGUSER("Copying white noise to mem cache...");
#pragma omp parallel for
for (int i = 0; i < nx; ++i)
for (int j = 0; j < ny; ++j)
for (int k = 0; k < nz; ++k)
(*mem_cache_[ilevel - levelmin_])[((size_t)i * ny + (size_t)j) * nz + (size_t)k] =
(*prng)(i + i0, j + j0, k + k0);
}
}
void RNG_music::fill_grid(int ilevel, DensityGrid<real_t> &A)
{
if (!initialized_)
{
LOGERR("Call to RNG_music::fill_grid before call to RNG_music::initialize_for_grid_structure");
throw std::runtime_error("invalid call order for random number generator");
}
if (restart_)
LOGINFO("Attempting to restart using random numbers for level %d\n from file \'wnoise_%04d.bin\'.", ilevel,
ilevel);
if (disk_cached_)
{
char fname[128];
sprintf(fname, "wnoise_%04d.bin", ilevel);
LOGUSER("Loading white noise from file \'%s\'...", fname);
std::ifstream ifs(fname, std::ios::binary);
if (!ifs.good())
{
LOGERR("White noise file \'%s\'was not found.", fname);
throw std::runtime_error("A white noise file was not found. This is an internal inconsistency and bad.");
}
int nx, ny, nz;
ifs.read(reinterpret_cast<char *>(&nx), sizeof(int));
ifs.read(reinterpret_cast<char *>(&ny), sizeof(int));
ifs.read(reinterpret_cast<char *>(&nz), sizeof(int));
if (nx != (int)A.size(0) || ny != (int)A.size(1) || nz != (int)A.size(2))
{
std::array<int,3> margin;
if( prefh_->get_margin()>0 ){
margin[0] = prefh_->get_margin();
margin[1] = prefh_->get_margin();
margin[2] = prefh_->get_margin();
}else{
margin[0] = prefh_->size(ilevel, 0)/2;
margin[1] = prefh_->size(ilevel, 1)/2;
margin[2] = prefh_->size(ilevel, 2)/2;
}
if (nx == (int)A.size(0) + 2 * margin[0] && ny == (int)A.size(1) + 2 * margin[1] && nz == (int)A.size(2) + 2 * margin[2])
{
int ox = margin[0], oy = margin[1], oz = margin[2];
std::vector<real_t> slice(ny * nz, 0.0);
for (int i = 0; i < nx; ++i)
{
ifs.read(reinterpret_cast<char *>(&slice[0]), ny * nz * sizeof(real_t));
if (i < ox)
continue;
if (i >= 3 * ox)
break;
#pragma omp parallel for
for (int j = oy; j < 3 * oy; ++j)
for (int k = oz; k < 3 * oz; ++k)
A(i - ox, j - oy, k - oz) = slice[j * nz + k];
}
ifs.close();
}
else
{
LOGERR("White noise file is not aligned with array. File: [%d,%d,%d]. Mem: [%d,%d,%d].", nx, ny, nz, A.size(0),
A.size(1), A.size(2));
throw std::runtime_error(
"White noise file is not aligned with array. This is an internal inconsistency and bad.");
}
}
else
{
for (int i = 0; i < nx; ++i)
{
std::vector<real_t> slice(ny * nz, 0.0);
ifs.read(reinterpret_cast<char *>(&slice[0]), ny * nz * sizeof(real_t));
#pragma omp parallel for
for (int j = 0; j < ny; ++j)
for (int k = 0; k < nz; ++k)
A(i, j, k) = slice[j * nz + k];
}
ifs.close();
}
}
else
{
LOGUSER("Copying white noise from memory cache...");
if (mem_cache_[ilevel - levelmin_] == NULL)
LOGERR("Tried to access mem-cached random numbers for level %d. But these are not available!\n", ilevel);
int nx(A.size(0)), ny(A.size(1)), nz(A.size(2));
if ((size_t)nx * (size_t)ny * (size_t)nz != mem_cache_[ilevel - levelmin_]->size())
{
LOGERR("White noise file is not aligned with array. File: [%d,%d,%d]. Mem: [%d,%d,%d].", nx, ny, nz, A.size(0),
A.size(1), A.size(2));
throw std::runtime_error("White noise file is not aligned with array. This is an internal inconsistency and bad");
}
#pragma omp parallel for
for (int i = 0; i < nx; ++i)
for (int j = 0; j < ny; ++j)
for (int k = 0; k < nz; ++k)
A(i, j, k) = (*mem_cache_[ilevel - levelmin_])[((size_t)i * ny + (size_t)j) * nz + (size_t)k];
std::vector<real_t>().swap(*mem_cache_[ilevel - levelmin_]);
delete mem_cache_[ilevel - levelmin_];
mem_cache_[ilevel - levelmin_] = NULL;
}
};
namespace
{
RNG_plugin_creator_concrete<RNG_music> creator("MUSIC");
}

895
src/plugins/random_music_wnoise_generator.cc Normal file
View file

@ -0,0 +1,895 @@
#include <complex>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include "random.hh"
#include "random_music_wnoise_generator.hh"
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)
{
LOGINFO("Generating random numbers (1) with seed %ld", baseseed);
initialize();
fill_subvolume(x0, lx);
}
template <typename T>
music_wnoise_generator<T>::music_wnoise_generator(unsigned res, unsigned cubesize, long baseseed, bool zeromean)
: res_(res), cubesize_(cubesize), ncubes_(1), baseseed_(baseseed)
{
LOGINFO("Generating random numbers (2) with seed %ld", baseseed);
double mean = 0.0;
size_t res_l = res;
bool musicnoise = true;
if (!musicnoise)
cubesize_ = res_;
if (!musicnoise)
LOGERR("This currently breaks compatibility. Need to disable by hand! Make sure to not check into repo");
initialize();
if (musicnoise)
mean = fill_all();
else
{
rnums_.push_back(new Meshvar<T>(res, 0, 0, 0));
cubemap_[0] = 0; // create dummy map index
register_cube(0, 0, 0);
// rapid_proto_ngenic_rng( res_, baseseed_, *this );
}
/*
if( musicnoise )
mean = fill_all();
else
{
rnums_.push_back( new Meshvar<T>( res, 0, 0, 0 ) );
cubemap_[0] = 0; // create dummy map index
register_cube(0,0,0);
rapid_proto_ngenic_rng( res_, baseseed_, *this );
}
*/
if (zeromean)
{
mean = 0.0;
#pragma omp parallel for reduction(+ \
: mean)
for (int i = 0; i < (int)res_; ++i)
for (unsigned j = 0; j < res_; ++j)
for (unsigned k = 0; k < res_; ++k)
mean += (*this)(i, j, k);
mean *= 1.0 / (double)(res_l * res_l * res_l);
#pragma omp parallel for
for (int i = 0; i < (int)res_; ++i)
for (unsigned j = 0; j < res_; ++j)
for (unsigned k = 0; k < res_; ++k)
(*this)(i, j, k) = (*this)(i, j, k) - mean;
}
}
template <typename T>
music_wnoise_generator<T>::music_wnoise_generator(unsigned res, std::string randfname, bool randsign)
: res_(res), cubesize_(res), ncubes_(1)
{
rnums_.push_back(new Meshvar<T>(res, 0, 0, 0));
cubemap_[0] = 0; // create dummy map index
std::ifstream ifs(randfname.c_str(), std::ios::binary);
if (!ifs)
{
LOGERR("Could not open random number file \'%s\'!", randfname.c_str());
throw std::runtime_error(std::string("Could not open random number file \'") + randfname + std::string("\'!"));
}
unsigned vartype;
unsigned nx, ny, nz, blksz32;
size_t blksz64;
int iseed;
// long seed;
float sign4 = -1.0f;
double sign8 = -1.0;
int addrtype = 32;
if (randsign) // use grafic2 sign convention
{
sign4 = 1.0f;
sign8 = 1.0;
}
//... read header and check if 32bit or 64bit block size .../
ifs.read(reinterpret_cast<char *>(&blksz32), sizeof(int));
ifs.read(reinterpret_cast<char *>(&nx), sizeof(unsigned));
if (blksz32 != 4 * sizeof(int) || nx != res_)
{
addrtype = 64;
ifs.seekg(0);
ifs.read(reinterpret_cast<char *>(&blksz64), sizeof(size_t));
ifs.read(reinterpret_cast<char *>(&nx), sizeof(unsigned));
if (blksz64 != 4 * sizeof(int) || nx != res_)
addrtype = -1;
}
ifs.seekg(0);
if (addrtype < 0)
throw std::runtime_error("corrupt random number file");
if (addrtype == 32)
ifs.read(reinterpret_cast<char *>(&blksz32), sizeof(int));
else
ifs.read(reinterpret_cast<char *>(&blksz64), sizeof(size_t));
ifs.read(reinterpret_cast<char *>(&nx), sizeof(unsigned));
ifs.read(reinterpret_cast<char *>(&ny), sizeof(unsigned));
ifs.read(reinterpret_cast<char *>(&nz), sizeof(unsigned));
ifs.read(reinterpret_cast<char *>(&iseed), sizeof(int));
// seed = (long)iseed;
if (nx != res_ || ny != res_ || nz != res_)
{
char errmsg[128];
sprintf(errmsg, "White noise file dimensions do not match level dimensions: %ux%ux%u vs. %u**3", nx, ny, nz, res_);
throw std::runtime_error(errmsg);
}
if (addrtype == 32)
ifs.read(reinterpret_cast<char *>(&blksz32), sizeof(int));
else
ifs.read(reinterpret_cast<char *>(&blksz64), sizeof(size_t));
//... read data ...//
// check whether random numbers are single or double precision numbers
if (addrtype == 32)
{
ifs.read(reinterpret_cast<char *>(&blksz32), sizeof(int));
if (blksz32 == nx * ny * sizeof(float))
vartype = 4;
else if (blksz32 == nx * ny * sizeof(double))
vartype = 8;
else
throw std::runtime_error("corrupt random number file");
}
else
{
ifs.read(reinterpret_cast<char *>(&blksz64), sizeof(size_t));
if (blksz64 == nx * ny * sizeof(float))
vartype = 4;
else if (blksz64 == nx * ny * sizeof(double))
vartype = 8;
else
throw std::runtime_error("corrupt random number file");
}
// rewind to beginning of block
if (addrtype == 32)
ifs.seekg(-sizeof(int), std::ios::cur);
else
ifs.seekg(-sizeof(size_t), std::ios::cur);
std::vector<float> in_float;
std::vector<double> in_double;
LOGINFO("Random number file \'%s\'\n contains %ld numbers. Reading...", randfname.c_str(), nx * ny * nz);
long double sum = 0.0, sum2 = 0.0;
size_t count = 0;
// perform actual reading
if (vartype == 4)
{
for (int ii = 0; ii < (int)nz; ++ii)
{
if (addrtype == 32)
{
ifs.read(reinterpret_cast<char *>(&blksz32), sizeof(int));
if (blksz32 != nx * ny * sizeof(float))
throw std::runtime_error("corrupt random number file");
}
else
{
ifs.read(reinterpret_cast<char *>(&blksz64), sizeof(size_t));
if (blksz64 != nx * ny * sizeof(float))
throw std::runtime_error("corrupt random number file");
}
in_float.assign(nx * ny, 0.0f);
ifs.read((char *)&in_float[0], nx * ny * sizeof(float));
for (int jj = 0, q = 0; jj < (int)ny; ++jj)
for (int kk = 0; kk < (int)nx; ++kk)
{
sum += in_float[q];
sum2 += in_float[q] * in_float[q];
++count;
(*rnums_[0])(kk, jj, ii) = sign4 * in_float[q++];
}
if (addrtype == 32)
{
ifs.read(reinterpret_cast<char *>(&blksz32), sizeof(int));
if (blksz32 != nx * ny * sizeof(float))
throw std::runtime_error("corrupt random number file");
}
else
{
ifs.read(reinterpret_cast<char *>(&blksz64), sizeof(size_t));
if (blksz64 != nx * ny * sizeof(float))
throw std::runtime_error("corrupt random number file");
}
}
}
else if (vartype == 8)
{
for (int ii = 0; ii < (int)nz; ++ii)
{
if (addrtype == 32)
{
ifs.read(reinterpret_cast<char *>(&blksz32), sizeof(int));
if (blksz32 != nx * ny * sizeof(double))
throw std::runtime_error("corrupt random number file");
}
else
{
ifs.read(reinterpret_cast<char *>(&blksz64), sizeof(size_t));
if (blksz64 != nx * ny * sizeof(double))
throw std::runtime_error("corrupt random number file");
}
in_double.assign(nx * ny, 0.0f);
ifs.read((char *)&in_double[0], nx * ny * sizeof(double));
for (int jj = 0, q = 0; jj < (int)ny; ++jj)
for (int kk = 0; kk < (int)nx; ++kk)
{
sum += in_double[q];
sum2 += in_double[q] * in_double[q];
++count;
(*rnums_[0])(kk, jj, ii) = sign8 * in_double[q++];
}
if (addrtype == 32)
{
ifs.read(reinterpret_cast<char *>(&blksz32), sizeof(int));
if (blksz32 != nx * ny * sizeof(double))
throw std::runtime_error("corrupt random number file");
}
else
{
ifs.read(reinterpret_cast<char *>(&blksz64), sizeof(size_t));
if (blksz64 != nx * ny * sizeof(double))
throw std::runtime_error("corrupt random number file");
}
}
}
double mean, var;
mean = sum / count;
var = sum2 / count - mean * mean;
LOGINFO("Random numbers in file have \n mean = %f and var = %f", mean, var);
}
//... copy construct by averaging down
template <typename T>
music_wnoise_generator<T>::music_wnoise_generator(/*const*/ music_wnoise_generator<T> &rc)
{
// if( res > rc.m_res || (res/rc.m_res)%2 != 0 )
// throw std::runtime_error("Invalid restriction in random number container copy constructor.");
long double sum = 0.0, sum2 = 0.0;
size_t count = 0;
LOGINFO("Generating a coarse white noise field by k-space degrading");
//... initialize properties of container
res_ = rc.res_ / 2;
cubesize_ = res_;
ncubes_ = 1;
baseseed_ = -2;
if (sizeof(fftw_real) != sizeof(T))
{
LOGERR("type mismatch with fftw_real in k-space averaging");
throw std::runtime_error("type mismatch with fftw_real in k-space averaging");
}
fftw_real
*rfine = new fftw_real[(size_t)rc.res_ * (size_t)rc.res_ * 2 * ((size_t)rc.res_ / 2 + 1)],
*rcoarse = new fftw_real[(size_t)res_ * (size_t)res_ * 2 * ((size_t)res_ / 2 + 1)];
fftw_complex
*ccoarse = reinterpret_cast<fftw_complex *>(rcoarse),
*cfine = reinterpret_cast<fftw_complex *>(rfine);
int nx(rc.res_), ny(rc.res_), nz(rc.res_), nxc(res_), nyc(res_), nzc(res_);
#ifdef FFTW3
#ifdef SINGLE_PRECISION
fftwf_plan
pf = fftwf_plan_dft_r2c_3d(nx, ny, nz, rfine, cfine, FFTW_ESTIMATE),
ipc = fftwf_plan_dft_c2r_3d(nxc, nyc, nzc, ccoarse, rcoarse, FFTW_ESTIMATE);
#else
fftw_plan
pf = fftw_plan_dft_r2c_3d(nx, ny, nz, rfine, cfine, FFTW_ESTIMATE),
ipc = fftw_plan_dft_c2r_3d(nxc, nyc, nzc, ccoarse, rcoarse, FFTW_ESTIMATE);
#endif
#else
rfftwnd_plan
pf = rfftw3d_create_plan(nx, ny, nz, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE | FFTW_IN_PLACE),
ipc = rfftw3d_create_plan(nxc, nyc, nzc, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE | FFTW_IN_PLACE);
#endif
#pragma omp parallel for
for (int i = 0; i < nx; i++)
for (int j = 0; j < ny; j++)
for (int k = 0; k < nz; k++)
{
size_t q = ((size_t)i * ny + (size_t)j) * (nz + 2) + (size_t)k;
rfine[q] = rc(i, j, k);
}
#ifdef FFTW3
#ifdef SINGLE_PRECISION
fftwf_execute(pf);
#else
fftw_execute(pf);
#endif
#else
#ifndef SINGLETHREAD_FFTW
rfftwnd_threads_one_real_to_complex(omp_get_max_threads(), pf, rfine, NULL);
#else
rfftwnd_one_real_to_complex(pf, rfine, NULL);
#endif
#endif
double fftnorm = 1.0 / ((double)nxc * (double)nyc * (double)nzc);
#pragma omp parallel for
for (int i = 0; i < nxc; i++)
for (int j = 0; j < nyc; j++)
for (int k = 0; k < nzc / 2 + 1; k++)
{
int ii(i), jj(j), kk(k);
if (i > nxc / 2)
ii += nx / 2;
if (j > nyc / 2)
jj += ny / 2;
size_t qc, qf;
double kx = (i <= (int)nxc / 2) ? (double)i : (double)(i - (int)nxc);
double ky = (j <= (int)nyc / 2) ? (double)j : (double)(j - (int)nyc);
double kz = (k <= (int)nzc / 2) ? (double)k : (double)(k - (int)nzc);
qc = ((size_t)i * nyc + (size_t)j) * (nzc / 2 + 1) + (size_t)k;
qf = ((size_t)ii * ny + (size_t)jj) * (nz / 2 + 1) + (size_t)kk;
std::complex<double> val_fine(RE(cfine[qf]), IM(cfine[qf]));
double phase = (kx / nxc + ky / nyc + kz / nzc) * 0.5 * M_PI;
std::complex<double> val_phas(cos(phase), sin(phase));
val_fine *= val_phas * fftnorm / sqrt(8.0);
RE(ccoarse[qc]) = val_fine.real();
IM(ccoarse[qc]) = val_fine.imag();
}
delete[] rfine;
#ifdef FFTW3
#ifdef SINGLE_PRECISION
fftwf_execute(ipc);
#else
fftw_execute(ipc);
#endif
#else
#ifndef SINGLETHREAD_FFTW
rfftwnd_threads_one_complex_to_real(omp_get_max_threads(), ipc, ccoarse, NULL);
#else
rfftwnd_one_complex_to_real(ipc, ccoarse, NULL);
#endif
#endif
rnums_.push_back(new Meshvar<T>(res_, 0, 0, 0));
cubemap_[0] = 0; // map all to single array
#pragma omp parallel for reduction(+ \
: sum, sum2, count)
for (int i = 0; i < nxc; i++)
for (int j = 0; j < nyc; j++)
for (int k = 0; k < nzc; k++)
{
size_t q = ((size_t)i * nyc + (size_t)j) * (nzc + 2) + (size_t)k;
(*rnums_[0])(i, j, k) = rcoarse[q];
sum += (*rnums_[0])(i, j, k);
sum2 += (*rnums_[0])(i, j, k) * (*rnums_[0])(i, j, k);
++count;
}
delete[] rcoarse;
#ifdef FFTW3
#ifdef SINGLE_PRECISION
fftwf_destroy_plan(pf);
fftwf_destroy_plan(ipc);
#else
fftw_destroy_plan(pf);
fftw_destroy_plan(ipc);
#endif
#else
rfftwnd_destroy_plan(pf);
rfftwnd_destroy_plan(ipc);
#endif
double rmean, rvar;
rmean = sum / count;
rvar = sum2 / count - rmean * rmean;
LOGINFO("Restricted random numbers have\n mean = %f, var = %f", rmean, rvar);
}
template <typename T>
music_wnoise_generator<T>::music_wnoise_generator(music_wnoise_generator<T> &rc, unsigned cubesize, long baseseed, int *x0_, int *lx_, bool zeromean)
: res_(2 * rc.res_), cubesize_(cubesize), ncubes_(1), baseseed_(baseseed)
{
initialize();
int x0[3], lx[3];
if (x0_ == NULL || lx_ == NULL)
{
for (int i = 0; i < 3; ++i)
{
x0[i] = 0;
lx[i] = res_;
}
fill_all();
}
else
{
for (int i = 0; i < 3; ++i)
{
x0[i] = x0_[i];
lx[i] = lx_[i];
}
fill_subvolume(x0, lx);
}
LOGINFO("Generating a constrained random number set with seed %ld\n using coarse mode replacement...", baseseed);
assert(lx[0] % 2 == 0 && lx[1] % 2 == 0 && lx[2] % 2 == 0);
size_t nx = lx[0], ny = lx[1], nz = lx[2],
nxc = lx[0] / 2, nyc = lx[1] / 2, nzc = lx[2] / 2;
fftw_real *rfine = new fftw_real[nx * ny * (nz + 2l)];
fftw_complex *cfine = reinterpret_cast<fftw_complex *>(rfine);
#ifdef FFTW3
#ifdef SINGLE_PRECISION
fftwf_plan
pf = fftwf_plan_dft_r2c_3d(nx, ny, nz, rfine, cfine, FFTW_ESTIMATE),
ipf = fftwf_plan_dft_c2r_3d(nx, ny, nz, cfine, rfine, FFTW_ESTIMATE);
#else
fftw_plan
pf = fftw_plan_dft_r2c_3d(nx, ny, nz, rfine, cfine, FFTW_ESTIMATE),
ipf = fftw_plan_dft_c2r_3d(nx, ny, nz, cfine, rfine, FFTW_ESTIMATE);
#endif
#else
rfftwnd_plan
pf = rfftw3d_create_plan(nx, ny, nz, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE | FFTW_IN_PLACE),
ipf = rfftw3d_create_plan(nx, ny, nz, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE | FFTW_IN_PLACE);
#endif
#pragma omp parallel for
for (int i = 0; i < (int)nx; i++)
for (int j = 0; j < (int)ny; j++)
for (int k = 0; k < (int)nz; k++)
{
size_t q = ((size_t)i * (size_t)ny + (size_t)j) * (size_t)(nz + 2) + (size_t)k;
rfine[q] = (*this)(x0[0] + i, x0[1] + j, x0[2] + k);
}
// this->free_all_mem(); // temporarily free memory, allocate again later
fftw_real *rcoarse = new fftw_real[nxc * nyc * (nzc + 2)];
fftw_complex *ccoarse = reinterpret_cast<fftw_complex *>(rcoarse);
#ifdef FFTW3
#ifdef SINGLE_PRECISION
fftwf_plan pc = fftwf_plan_dft_r2c_3d(nxc, nyc, nzc, rcoarse, ccoarse, FFTW_ESTIMATE);
#else
fftw_plan pc = fftw_plan_dft_r2c_3d(nxc, nyc, nzc, rcoarse, ccoarse, FFTW_ESTIMATE);
#endif
#else
rfftwnd_plan pc = rfftw3d_create_plan(nxc, nyc, nzc, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE | FFTW_IN_PLACE);
#endif
#pragma omp parallel for
for (int i = 0; i < (int)nxc; i++)
for (int j = 0; j < (int)nyc; j++)
for (int k = 0; k < (int)nzc; k++)
{
size_t q = ((size_t)i * (size_t)nyc + (size_t)j) * (size_t)(nzc + 2) + (size_t)k;
rcoarse[q] = rc(x0[0] / 2 + i, x0[1] / 2 + j, x0[2] / 2 + k);
}
#ifdef FFTW3
#ifdef SINGLE_PRECISION
fftwf_execute(pc);
fftwf_execute(pf);
#else
fftw_execute(pc);
fftw_execute(pf);
#endif
#else
#ifndef SINGLETHREAD_FFTW
rfftwnd_threads_one_real_to_complex(omp_get_max_threads(), pc, rcoarse, NULL);
rfftwnd_threads_one_real_to_complex(omp_get_max_threads(), pf, rfine, NULL);
#else
rfftwnd_one_real_to_complex(pc, rcoarse, NULL);
rfftwnd_one_real_to_complex(pf, rfine, NULL);
#endif
#endif
double fftnorm = 1.0 / ((double)nx * (double)ny * (double)nz);
double sqrt8 = sqrt(8.0);
double phasefac = -0.5;
// embedding of coarse white noise by fourier interpolation
#pragma omp parallel for
for (int i = 0; i < (int)nxc; i++)
for (int j = 0; j < (int)nyc; j++)
for (int k = 0; k < (int)nzc / 2 + 1; k++)
{
int ii(i), jj(j), kk(k);
// if( i==(int)nxc/2 ) continue;
// if( j==(int)nyc/2 ) continue;
if (i > (int)nxc / 2)
ii += (int)nx / 2;
if (j > (int)nyc / 2)
jj += (int)ny / 2;
size_t qc, qf;
double kx = (i <= (int)nxc / 2) ? (double)i : (double)(i - (int)nxc);
double ky = (j <= (int)nyc / 2) ? (double)j : (double)(j - (int)nyc);
double kz = (k <= (int)nzc / 2) ? (double)k : (double)(k - (int)nzc);
qc = ((size_t)i * nyc + (size_t)j) * (nzc / 2 + 1) + (size_t)k;
qf = ((size_t)ii * ny + (size_t)jj) * (nz / 2 + 1) + (size_t)kk;
std::complex<double> val(RE(ccoarse[qc]), IM(ccoarse[qc]));
double phase = (kx / nxc + ky / nyc + kz / nzc) * phasefac * M_PI;
std::complex<double> val_phas(cos(phase), sin(phase));
val *= val_phas * sqrt8;
if (i != (int)nxc / 2 && j != (int)nyc / 2 && k != (int)nzc / 2)
{
RE(cfine[qf]) = val.real();
IM(cfine[qf]) = val.imag();
}
else
{
// RE(cfine[qf]) = val.real();
// IM(cfine[qf]) = 0.0;
}
}
delete[] rcoarse;
#pragma omp parallel for
for (int i = 0; i < (int)nx; i++)
for (int j = 0; j < (int)ny; j++)
for (int k = 0; k < (int)nz / 2 + 1; k++)
{
size_t q = ((size_t)i * ny + (size_t)j) * (nz / 2 + 1) + (size_t)k;
RE(cfine[q]) *= fftnorm;
IM(cfine[q]) *= fftnorm;
}
#ifdef FFTW3
#ifdef SINGLE_PRECISION
fftwf_execute(ipf);
#else
fftw_execute(ipf);
#endif
#else
#ifndef SINGLETHREAD_FFTW
rfftwnd_threads_one_complex_to_real(omp_get_max_threads(), ipf, cfine, NULL);
#else
rfftwnd_one_complex_to_real(ipf, cfine, NULL);
#endif
#endif
#pragma omp parallel for
for (int i = 0; i < (int)nx; i++)
for (int j = 0; j < (int)ny; j++)
for (int k = 0; k < (int)nz; k++)
{
size_t q = ((size_t)i * ny + (size_t)j) * (nz + 2) + (size_t)k;
(*this)(x0[0] + i, x0[1] + j, x0[2] + k, false) = rfine[q];
}
delete[] rfine;
#ifdef FFTW3
#ifdef SINGLE_PRECISION
fftwf_destroy_plan(pf);
fftwf_destroy_plan(pc);
fftwf_destroy_plan(ipf);
#else
fftw_destroy_plan(pf);
fftw_destroy_plan(pc);
fftw_destroy_plan(ipf);
#endif
#else
fftwnd_destroy_plan(pf);
fftwnd_destroy_plan(pc);
fftwnd_destroy_plan(ipf);
#endif
}
template <typename T>
void music_wnoise_generator<T>::register_cube(int i, int j, int k)
{
i = (i + ncubes_) % ncubes_;
j = (j + ncubes_) % ncubes_;
k = (k + ncubes_) % ncubes_;
size_t icube = ((size_t)i * ncubes_ + (size_t)j) * ncubes_ + (size_t)k;
cubemap_iterator it = cubemap_.find(icube);
if (it == cubemap_.end())
{
rnums_.push_back(NULL);
cubemap_[icube] = rnums_.size() - 1;
#ifdef DEBUG
LOGDEBUG("registering new cube %d,%d,%d . ID = %ld, memloc = %ld", i, j, k, icube, cubemap_[icube]);
#endif
}
}
template <typename T>
double music_wnoise_generator<T>::fill_cube(int i, int j, int k)
{
gsl_rng *RNG = gsl_rng_alloc(gsl_rng_mt19937);
i = (i + ncubes_) % ncubes_;
j = (j + ncubes_) % ncubes_;
k = (k + ncubes_) % ncubes_;
size_t icube = ((size_t)i * ncubes_ + (size_t)j) * ncubes_ + (size_t)k;
long cubeseed = baseseed_ + icube; //... each cube gets its unique seed
gsl_rng_set(RNG, cubeseed);
cubemap_iterator it = cubemap_.find(icube);
if (it == cubemap_.end())
{
LOGERR("Attempt to access non-registered random number cube!");
throw std::runtime_error("Attempt to access non-registered random number cube!");
}
size_t cubeidx = it->second;
if (rnums_[cubeidx] == NULL)
rnums_[cubeidx] = new Meshvar<T>(cubesize_, 0, 0, 0);
double mean = 0.0;
for (int ii = 0; ii < (int)cubesize_; ++ii)
for (int jj = 0; jj < (int)cubesize_; ++jj)
for (int kk = 0; kk < (int)cubesize_; ++kk)
{
(*rnums_[cubeidx])(ii, jj, kk) = gsl_ran_ugaussian_ratio_method(RNG);
mean += (*rnums_[cubeidx])(ii, jj, kk);
}
gsl_rng_free(RNG);
return mean / (cubesize_ * cubesize_ * cubesize_);
}
template <typename T>
void music_wnoise_generator<T>::subtract_from_cube(int i, int j, int k, double val)
{
i = (i + ncubes_) % ncubes_;
j = (j + ncubes_) % ncubes_;
k = (k + ncubes_) % ncubes_;
size_t icube = ((size_t)i * ncubes_ + (size_t)j) * ncubes_ + (size_t)k;
cubemap_iterator it = cubemap_.find(icube);
if (it == cubemap_.end())
{
LOGERR("Attempt to access unallocated RND cube %d,%d,%d in music_wnoise_generator::subtract_from_cube", i, j, k);
throw std::runtime_error("Attempt to access unallocated RND cube in music_wnoise_generator::subtract_from_cube");
}
size_t cubeidx = it->second;
for (int ii = 0; ii < (int)cubesize_; ++ii)
for (int jj = 0; jj < (int)cubesize_; ++jj)
for (int kk = 0; kk < (int)cubesize_; ++kk)
(*rnums_[cubeidx])(ii, jj, kk) -= val;
}
template <typename T>
void music_wnoise_generator<T>::free_cube(int i, int j, int k)
{
i = (i + ncubes_) % ncubes_;
j = (j + ncubes_) % ncubes_;
k = (k + ncubes_) % ncubes_;
size_t icube = ((size_t)i * (size_t)ncubes_ + (size_t)j) * (size_t)ncubes_ + (size_t)k;
cubemap_iterator it = cubemap_.find(icube);
if (it == cubemap_.end())
{
LOGERR("Attempt to access unallocated RND cube %d,%d,%d in music_wnoise_generator::free_cube", i, j, k);
throw std::runtime_error("Attempt to access unallocated RND cube in music_wnoise_generator::free_cube");
}
size_t cubeidx = it->second;
if (rnums_[cubeidx] != NULL)
{
delete rnums_[cubeidx];
rnums_[cubeidx] = NULL;
}
}
template <typename T>
void music_wnoise_generator<T>::initialize(void)
{
ncubes_ = std::max((int)((double)res_ / cubesize_), 1);
if (res_ < cubesize_)
{
ncubes_ = 1;
cubesize_ = res_;
}
LOGINFO("Generating random numbers w/ sample cube size of %d", cubesize_);
}
template <typename T>
double music_wnoise_generator<T>::fill_subvolume(int *i0, int *n)
{
int i0cube[3], ncube[3];
i0cube[0] = (int)((double)(res_ + i0[0]) / cubesize_);
i0cube[1] = (int)((double)(res_ + i0[1]) / cubesize_);
i0cube[2] = (int)((double)(res_ + i0[2]) / cubesize_);
ncube[0] = (int)(n[0] / cubesize_) + 2;
ncube[1] = (int)(n[1] / cubesize_) + 2;
ncube[2] = (int)(n[2] / cubesize_) + 2;
#ifdef DEBUG
LOGDEBUG("random numbers needed for region %d,%d,%d ..+ %d,%d,%d", i0[0], i0[1], i0[2], n[0], n[1], n[2]);
LOGDEBUG("filling cubes %d,%d,%d ..+ %d,%d,%d", i0cube[0], i0cube[1], i0cube[2], ncube[0], ncube[1], ncube[2]);
#endif
double mean = 0.0;
for (int i = i0cube[0]; i < i0cube[0] + ncube[0]; ++i)
for (int j = i0cube[1]; j < i0cube[1] + ncube[1]; ++j)
for (int k = i0cube[2]; k < i0cube[2] + ncube[2]; ++k)
{
int ii(i), jj(j), kk(k);
ii = (ii + ncubes_) % ncubes_;
jj = (jj + ncubes_) % ncubes_;
kk = (kk + ncubes_) % ncubes_;
register_cube(ii, jj, kk);
}
#pragma omp parallel for reduction(+ : mean)
for (int i = i0cube[0]; i < i0cube[0] + ncube[0]; ++i)
for (int j = i0cube[1]; j < i0cube[1] + ncube[1]; ++j)
for (int k = i0cube[2]; k < i0cube[2] + ncube[2]; ++k)
{
int ii(i), jj(j), kk(k);
ii = (ii + ncubes_) % ncubes_;
jj = (jj + ncubes_) % ncubes_;
kk = (kk + ncubes_) % ncubes_;
mean += fill_cube(ii, jj, kk);
}
return mean / (ncube[0] * ncube[1] * ncube[2]);
}
template <typename T>
double music_wnoise_generator<T>::fill_all(void)
{
double sum = 0.0;
for (int i = 0; i < (int)ncubes_; ++i)
for (int j = 0; j < (int)ncubes_; ++j)
for (int k = 0; k < (int)ncubes_; ++k)
{
int ii(i), jj(j), kk(k);
ii = (ii + ncubes_) % ncubes_;
jj = (jj + ncubes_) % ncubes_;
kk = (kk + ncubes_) % ncubes_;
register_cube(ii, jj, kk);
}
#pragma omp parallel for reduction(+ \
: sum)
for (int i = 0; i < (int)ncubes_; ++i)
for (int j = 0; j < (int)ncubes_; ++j)
for (int k = 0; k < (int)ncubes_; ++k)
{
int ii(i), jj(j), kk(k);
ii = (ii + ncubes_) % ncubes_;
jj = (jj + ncubes_) % ncubes_;
kk = (kk + ncubes_) % ncubes_;
sum += fill_cube(ii, jj, kk);
}
//... subtract mean
#pragma omp parallel for reduction(+ \
: sum)
for (int i = 0; i < (int)ncubes_; ++i)
for (int j = 0; j < (int)ncubes_; ++j)
for (int k = 0; k < (int)ncubes_; ++k)
{
int ii(i), jj(j), kk(k);
ii = (ii + ncubes_) % ncubes_;
jj = (jj + ncubes_) % ncubes_;
kk = (kk + ncubes_) % ncubes_;
subtract_from_cube(ii, jj, kk, sum / (ncubes_ * ncubes_ * ncubes_));
}
return sum / (ncubes_ * ncubes_ * ncubes_);
}
template <typename T>
void music_wnoise_generator<T>::print_allocated(void)
{
unsigned ncount = 0, ntot = rnums_.size();
for (size_t i = 0; i < rnums_.size(); ++i)
if (rnums_[i] != NULL)
ncount++;
LOGINFO(" -> %d of %d random number cubes currently allocated", ncount, ntot);
}
template class music_wnoise_generator<float>;
template class music_wnoise_generator<double>;

205
src/plugins/random_music_wnoise_generator.hh Normal file
View file

@ -0,0 +1,205 @@
#ifndef __RANDOM_MUSIC_WNOISE_GENERATOR_HH
#define __RANDOM_MUSIC_WNOISE_GENERATOR_HH
#include <vector>
#include <map>
#include "general.hh"
#include "mesh.hh"
#define DEF_RAN_CUBE_SIZE 32
/*!
* @brief encapsulates all things random number generator related
*/
template< typename T >
class music_wnoise_generator
{
public:
unsigned
res_, //!< resolution of the full mesh
cubesize_, //!< size of one independent random number cube
ncubes_; //!< number of random number cubes to cover the full mesh
long baseseed_; //!< base seed from which cube seeds are computed
protected:
//! vector of 3D meshes (the random number cubes) with random numbers
std::vector< Meshvar<T>* > rnums_;
//! map of 3D indices to cube index
std::map<size_t,size_t> cubemap_;
typedef std::map<size_t,size_t>::iterator cubemap_iterator;
protected:
//! register a cube with the hash map
void register_cube( int i, int j, int k);
//! fills a subcube with random numbers
double fill_cube( int i, int j, int k);
//! subtract a constant from an entire cube
void subtract_from_cube( int i, int j, int k, double val );
//! copy random numbers from a cube to a full grid array
template< class C >
void copy_cube( int i, int j, int k, C& dat )
{
int offi, offj, offk;
offi = i*cubesize_;
offj = j*cubesize_;
offk = k*cubesize_;
i = (i+ncubes_)%ncubes_;
j = (j+ncubes_)%ncubes_;
k = (k+ncubes_)%ncubes_;
size_t icube = (i*ncubes_+j)*ncubes_+k;
cubemap_iterator it = cubemap_.find( icube );
if( it == cubemap_.end() )
{
LOGERR("attempting to copy data from non-existing RND cube %d,%d,%d",i,j,k);
throw std::runtime_error("attempting to copy data from non-existing RND cube");
}
size_t cubeidx = it->second;
for( int ii=0; ii<(int)cubesize_; ++ii )
for( int jj=0; jj<(int)cubesize_; ++jj )
for( int kk=0; kk<(int)cubesize_; ++kk )
dat(offi+ii,offj+jj,offk+kk) = (*rnums_[cubeidx])(ii,jj,kk);
}
//! free the memory associated with a subcube
void free_cube( int i, int j, int k );
//! initialize member variables and allocate memory
void initialize( void );
//! fill a cubic subvolume of the full grid with random numbers
double fill_subvolume( int *i0, int *n );
//! fill an entire grid with random numbers
double fill_all( void );
//! fill an external array instead of the internal field
template< class C >
double fill_all( C& dat )
{
double sum = 0.0;
for( int i=0; i<(int)ncubes_; ++i )
for( int j=0; j<(int)ncubes_; ++j )
for( int k=0; k<(int)ncubes_; ++k )
{
int ii(i),jj(j),kk(k);
register_cube(ii,jj,kk);
}
#pragma omp parallel for reduction(+:sum)
for( int i=0; i<(int)ncubes_; ++i )
for( int j=0; j<(int)ncubes_; ++j )
for( int k=0; k<(int)ncubes_; ++k )
{
int ii(i),jj(j),kk(k);
ii = (ii+ncubes_)%ncubes_;
jj = (jj+ncubes_)%ncubes_;
kk = (kk+ncubes_)%ncubes_;
sum+=fill_cube(ii, jj, kk);
copy_cube(ii,jj,kk,dat);
free_cube(ii, jj, kk);
}
return sum/(ncubes_*ncubes_*ncubes_);
}
//! write the number of allocated random number cubes to stdout
void print_allocated( void );
public:
//! constructor
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, int *x0_=NULL, int *lx_=NULL, bool zeromean=true );
//! constructor
music_wnoise_generator( unsigned res, unsigned cubesize, long baseseed, bool zeromean=true );
//! constructor to read white noise from file
music_wnoise_generator( unsigned res, std::string randfname, bool rndsign );
//! copy constructor for averaged field (not copying) hence explicit!
explicit music_wnoise_generator( /*const*/ music_wnoise_generator <T>& rc );
//! destructor
~music_wnoise_generator()
{
for( unsigned i=0; i<rnums_.size(); ++i )
if( rnums_[i] != NULL )
delete rnums_[i];
rnums_.clear();
}
//! access a random number, this allocates a cube and fills it with consistent random numbers
inline T& operator()( int i, int j, int k, bool fillrand=true )
{
int ic, jc, kc, is, js, ks;
if( ncubes_ == 0 )
throw std::runtime_error("music_wnoise_generator: internal error, not properly initialized");
//... determine cube
ic = (int)((double)i/cubesize_ + ncubes_) % ncubes_;
jc = (int)((double)j/cubesize_ + ncubes_) % ncubes_;
kc = (int)((double)k/cubesize_ + ncubes_) % ncubes_;
size_t icube = ((size_t)ic*ncubes_+(size_t)jc)*ncubes_+(size_t)kc;
cubemap_iterator it = cubemap_.find( icube );
if( it == cubemap_.end() )
{
LOGERR("Attempting to copy data from non-existing RND cube %d,%d,%d @ %d,%d,%d",ic,jc,kc,i,j,k);
throw std::runtime_error("attempting to copy data from non-existing RND cube");
}
size_t cubeidx = it->second;
if( rnums_[ cubeidx ] == NULL )
{
LOGERR("Attempting to access data from non-allocated RND cube %d,%d,%d",ic,jc,kc);
throw std::runtime_error("attempting to access data from non-allocated RND cube");
}
//... determine cell in cube
is = (i - ic * cubesize_ + cubesize_) % cubesize_;
js = (j - jc * cubesize_ + cubesize_) % cubesize_;
ks = (k - kc * cubesize_ + cubesize_) % cubesize_;
return (*rnums_[ cubeidx ])(is,js,ks);
}
//! free all cubes
void free_all_mem( void )
{
for( unsigned i=0; i<rnums_.size(); ++i )
if( rnums_[i] != NULL )
{
delete rnums_[i];
rnums_[i] = NULL;
}
}
};
#endif //__MUSIC_WNOISE_GENERATOR_HH

387
src/plugins/random_panphasia.cc
View file

@ -9,18 +9,20 @@
const int maxdim = 60, maxlev = 50, maxpow = 3 * maxdim;
typedef int rand_offset_[5];
typedef struct {
typedef struct
{
int state[133]; // Nstore = Nstate (=5) + Nbatch (=128)
int need_fill;
int pos;
} rand_state_;
} rand_state_;
/* pan_state_ struct -- corresponds to respective fortran module in panphasia_routines.f
* data structure that contains all panphasia state variables
* it needs to get passed between the fortran routines to enable
* thread-safe execution.
*/
typedef struct {
typedef struct
{
int base_state[5], base_lev_start[5][maxdim + 1];
rand_offset_ poweroffset[maxpow + 1], superjump;
rand_state_ current_state[maxpow + 2];
@ -45,29 +47,31 @@ typedef struct {
} pan_state_;
extern "C" {
void start_panphasia_(pan_state_ *lstate, const char *descriptor, int *ngrid, int *bverbose);
extern "C"
{
void start_panphasia_(pan_state_ *lstate, const char *descriptor, int *ngrid, int *bverbose);
void parse_descriptor_(const char *descriptor, int16_t *l, int32_t *ix, int32_t *iy, int32_t *iz, int16_t *side1,
int16_t *side2, int16_t *side3, int32_t *check_int, char *name);
void parse_descriptor_(const char *descriptor, int16_t *l, int32_t *ix, int32_t *iy, int32_t *iz, int16_t *side1,
int16_t *side2, int16_t *side3, int32_t *check_int, char *name);
void panphasia_cell_properties_(pan_state_ *lstate, int *ixcell, int *iycell, int *izcell, double *cell_prop);
void panphasia_cell_properties_(pan_state_ *lstate, int *ixcell, int *iycell, int *izcell, double *cell_prop);
void adv_panphasia_cell_properties_(pan_state_ *lstate, int *ixcell, int *iycell, int *izcell, int *layer_min,
int *layer_max, int *indep_field, double *cell_prop);
void adv_panphasia_cell_properties_(pan_state_ *lstate, int *ixcell, int *iycell, int *izcell, int *layer_min,
int *layer_max, int *indep_field, double *cell_prop);
void set_phases_and_rel_origin_(pan_state_ *lstate, const char *descriptor, int *lev, long long *ix_rel,
long long *iy_rel, long long *iz_rel, int *VERBOSE);
/*void set_local_box_( pan_state_ *lstate, int lev, int8_t ix_abs, int8_t iy_abs, int8_t iz_abs,
int8_t ix_per, int8_t iy_per, int8_t iz_per, int8_t ix_rel, int8_t iy_rel,
int8_t iz_rel, int wn_level_base, int8_t check_rand, char *phase_name, int MYID);*/
/*extern struct {
int layer_min, layer_max, hoswitch;
}oct_range_;
*/
void set_phases_and_rel_origin_(pan_state_ *lstate, const char *descriptor, int *lev, long long *ix_rel,
long long *iy_rel, long long *iz_rel, int *VERBOSE);
/*void set_local_box_( pan_state_ *lstate, int lev, int8_t ix_abs, int8_t iy_abs, int8_t iz_abs,
int8_t ix_per, int8_t iy_per, int8_t iz_per, int8_t ix_rel, int8_t iy_rel,
int8_t iz_rel, int wn_level_base, int8_t check_rand, char *phase_name, int MYID);*/
/*extern struct {
int layer_min, layer_max, hoswitch;
}oct_range_;
*/
}
class RNG_panphasia : public RNG_plugin {
class RNG_panphasia : public RNG_plugin
{
private:
void forward_transform_field(real_t *field, int n0, int n1, int n2);
void forward_transform_field(real_t *field, int n) { forward_transform_field(field, n, n, n); }
@ -83,26 +87,29 @@ protected:
double inter_grid_phase_adjustment_;
// double translation_phase_;
pan_state_ *lstate;
int grid_p_,grid_m_;
int grid_p_, grid_m_;
double grid_rescale_fac_;
int coordinate_system_shift_[3];
int ix_abs_[3], ix_per_[3], ix_rel_[3], level_p_, lextra_;
const refinement_hierarchy *prefh_;
struct panphasia_descriptor {
struct panphasia_descriptor
{
int16_t wn_level_base;
int32_t i_xorigin_base, i_yorigin_base, i_zorigin_base;
int16_t i_base, i_base_y, i_base_z;
int32_t check_rand;
std::string name;
explicit panphasia_descriptor(std::string dstring) {
explicit panphasia_descriptor(std::string dstring)
{
char tmp[100];
memset(tmp, ' ', 100);
parse_descriptor_(dstring.c_str(), &wn_level_base, &i_xorigin_base, &i_yorigin_base, &i_zorigin_base, &i_base,
&i_base_y, &i_base_z, &check_rand, tmp);
for (int i = 0; i < 100; i++)
if (tmp[i] == ' ') {
if (tmp[i] == ' ')
{
tmp[i] = '\0';
break;
}
@ -111,8 +118,10 @@ protected:
}
};
void clear_panphasia_thread_states(void) {
for (int i = 0; i < num_threads_; ++i) {
void clear_panphasia_thread_states(void)
{
for (int i = 0; i < num_threads_; ++i)
{
lstate[i].init = 0;
lstate[i].init_cell_props = 0;
lstate[i].init_lecuyer_state = 0;
@ -120,7 +129,8 @@ protected:
}
// greatest common divisor
int gcd(int a, int b) {
int gcd(int a, int b)
{
if (b == 0)
return a;
return gcd(b, a % b);
@ -129,19 +139,22 @@ protected:
// least common multiple
int lcm(int a, int b) { return abs(a * b) / gcd(a, b); }
// Two or largest power of 2 less than the argument
int largest_power_two_lte(int b) {
// Two or largest power of 2 less than the argument
int largest_power_two_lte(int b)
{
int a = 1;
if (b<=a) return a;
while (2*a < b) a = 2*a;
if (b <= a)
return a;
while (2 * a < b)
a = 2 * a;
return a;
}
}
panphasia_descriptor *pdescriptor_;
public:
explicit RNG_panphasia(config_file &cf) : RNG_plugin(cf) {
explicit RNG_panphasia(config_file &cf) : RNG_plugin(cf)
{
descriptor_string_ = pcf_->getValue<std::string>("random", "descriptor");
#ifdef _OPENMP
@ -160,16 +173,17 @@ public:
// write panphasia base size into config file for the grid construction
// as the gridding unit we use the least common multiple of 2 and i_base
std::stringstream ss;
//ARJ ss << lcm(2, pdescriptor_->i_base);
//ss << two_or_largest_power_two_less_than(pdescriptor_->i_base);//ARJ
ss << 2; //ARJ - set gridding unit to two
// ARJ ss << lcm(2, pdescriptor_->i_base);
// ss << two_or_largest_power_two_less_than(pdescriptor_->i_base);//ARJ
ss << 2; // ARJ - set gridding unit to two
pcf_->insertValue("setup", "gridding_unit", ss.str());
ss.str(std::string());
ss << pdescriptor_->i_base ;
pcf_->insertValue("random","base_unit", ss.str());
ss << pdescriptor_->i_base;
pcf_->insertValue("random", "base_unit", ss.str());
}
void initialize_for_grid_structure(const refinement_hierarchy &refh) {
void initialize_for_grid_structure(const refinement_hierarchy &refh)
{
prefh_ = &refh;
levelmin_ = prefh_->levelmin();
levelmin_final_ = pcf_->getValue<unsigned>("setup", "levelmin");
@ -180,20 +194,21 @@ public:
// if ngrid is not a multiple of i_base, then we need to enlarge and then sample down
ngrid_ = 1 << levelmin_;
grid_p_ = pdescriptor_->i_base;
grid_m_ = largest_power_two_lte(grid_p_);
lextra_ = (log10((double)ngrid_ / (double)pdescriptor_->i_base) + 0.001) / log10(2.0);
int ratio = 1 << lextra_;
int ratio = 1 << lextra_;
grid_rescale_fac_ = 1.0;
coordinate_system_shift_[0] = -pcf_->getValue<int>("setup", "shift_x");
coordinate_system_shift_[1] = -pcf_->getValue<int>("setup", "shift_y");
coordinate_system_shift_[2] = -pcf_->getValue<int>("setup", "shift_z");
incongruent_fields_ = false;
if (ngrid_ != ratio * pdescriptor_->i_base) {
if (ngrid_ != ratio * pdescriptor_->i_base)
{
incongruent_fields_ = true;
ngrid_ = 2 * ratio * pdescriptor_->i_base;
grid_rescale_fac_ = (double)ngrid_ / (1 << levelmin_);
@ -201,7 +216,7 @@ public:
" (%d -> %d) * 2**ref compatible with PANPHASIA\n"
" will Fourier interpolate after.",
grid_m_, grid_p_);
}
}
}
~RNG_panphasia() { delete[] lstate; }
@ -211,7 +226,8 @@ public:
bool is_multiscale() const { return true; }
};
void RNG_panphasia::forward_transform_field(real_t *field, int nx, int ny, int nz) {
void RNG_panphasia::forward_transform_field(real_t *field, int nx, int ny, int nz)
{
fftw_real *rfield = reinterpret_cast<fftw_real *>(field);
fftw_complex *cfield = reinterpret_cast<fftw_complex *>(field);
@ -241,7 +257,8 @@ void RNG_panphasia::forward_transform_field(real_t *field, int nx, int ny, int n
#endif
}
void RNG_panphasia::backward_transform_field(real_t *field, int nx, int ny, int nz) {
void RNG_panphasia::backward_transform_field(real_t *field, int nx, int ny, int nz)
{
fftw_real *rfield = reinterpret_cast<fftw_real *>(field);
fftw_complex *cfield = reinterpret_cast<fftw_complex *>(field);
@ -272,7 +289,8 @@ void RNG_panphasia::backward_transform_field(real_t *field, int nx, int ny, int
}
#include <sys/time.h>
inline double get_wtime(void) {
inline double get_wtime(void)
{
#ifdef _OPENMP
return omp_get_wtime();
#else
@ -280,63 +298,69 @@ inline double get_wtime(void) {
#endif
}
void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R) {
void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
{
fftw_real *pr0, *pr1, *pr2, *pr3, *pr4;
fftw_complex *pc0, *pc1, *pc2, *pc3, *pc4;
// determine resolution and offset so that we can do proper resampling
int ileft[3], ileft_corner[3], nx[3], nxremap[3];
int iexpand_left[3];
for (int k = 0; k < 3; ++k) {
ileft[k] = prefh_->offset_abs(level, k);
for (int k = 0; k < 3; ++k)
{
ileft[k] = prefh_->offset_abs(level, k);
nx[k] = R.size(k);
assert(nx[k] % 4 == 0);
if (level == levelmin_) {
ileft_corner[k] = ileft[k]; // Top level - periodic
}else{
ileft_corner[k] = (ileft[k] - nx[k]/4 + (1 << level))%(1 << level); // Isolated
assert(nx[k] % 4 == 0);
if (level == levelmin_)
{
ileft_corner[k] = ileft[k]; // Top level - periodic
}
iexpand_left[k] = (ileft_corner[k]%grid_m_ ==0) ? 0 : ileft_corner[k]%grid_m_;
fprintf(stderr, "dim=%c : ileft = %d, ileft_corner %d, nx = %d\n", 'x' + k, ileft[k],ileft_corner[k],nx[k]);
else
{
ileft_corner[k] = (ileft[k] - nx[k] / 4 + (1 << level)) % (1 << level); // Isolated
}
iexpand_left[k] = (ileft_corner[k] % grid_m_ == 0) ? 0 : ileft_corner[k] % grid_m_;
// fprintf(stderr, "dim=%c : ileft = %d, ileft_corner %d, nx = %d\n", 'x' + k, ileft[k],ileft_corner[k],nx[k]);
};
int ileft_corner_m[3], ileft_corner_p[3],nx_m[3];
int ileft_max_expand = std::max(iexpand_left[0],std::max(iexpand_left[1],iexpand_left[2]));
int ileft_corner_m[3], ileft_corner_p[3], nx_m[3];
int ileft_max_expand = std::max(iexpand_left[0], std::max(iexpand_left[1], iexpand_left[2]));
for (int k = 0; k < 3; ++k) {
ileft_corner_m[k] = ((ileft_corner[k] - iexpand_left[k]) +
coordinate_system_shift_[k] * (1 << (level - levelmin_final_)) + (1 << level)) % (1 << level);
for (int k = 0; k < 3; ++k)
{
ileft_corner_m[k] = ((ileft_corner[k] - iexpand_left[k]) +
coordinate_system_shift_[k] * (1 << (level - levelmin_final_)) + (1 << level)) %
(1 << level);
ileft_corner_p[k] = grid_p_ * ileft_corner_m[k]/grid_m_;
nx_m[k] = (ileft_max_expand!=0)? nx[k] + ileft_max_expand: nx[k];
if (nx_m[k]%grid_m_ !=0) nx_m[k] = nx_m[k] + grid_m_ - nx_m[k]%grid_m_;
nxremap[k] = grid_p_ * nx_m[k]/grid_m_;
if (nxremap[k]%2==1){
ileft_corner_p[k] = grid_p_ * ileft_corner_m[k] / grid_m_;
nx_m[k] = (ileft_max_expand != 0) ? nx[k] + ileft_max_expand : nx[k];
if (nx_m[k] % grid_m_ != 0)
nx_m[k] = nx_m[k] + grid_m_ - nx_m[k] % grid_m_;
nxremap[k] = grid_p_ * nx_m[k] / grid_m_;
if (nxremap[k] % 2 == 1)
{
nx_m[k] = nx_m[k] + grid_m_;
nxremap[k] = grid_p_ * nx_m[k]/grid_m_;
}
nxremap[k] = grid_p_ * nx_m[k] / grid_m_;
}
}
if ( (nx_m[0]!=nx_m[1]) || (nx_m[0]!=nx_m[2])) LOGERR("Fatal error: non-cubic refinement being requested");
if ((nx_m[0] != nx_m[1]) || (nx_m[0] != nx_m[2]))
LOGERR("Fatal error: non-cubic refinement being requested");
inter_grid_phase_adjustment_ = M_PI * (1.0 / (double)nx_m[0] - 1.0 / (double)nxremap[0]);
LOGINFO("The value of the phase adjustement is %f\n", inter_grid_phase_adjustment_);
LOGUSER("The value of the phase adjustement is %f\n", inter_grid_phase_adjustment_);
// LOGINFO("ileft[0],ileft[1],ileft[2] %d %d %d", ileft[0], ileft[1], ileft[2]);
// LOGINFO("ileft_corner[0,1,2] %d %d %d", ileft_corner[0], ileft_corner[1], ileft_corner[2]);
LOGINFO("ileft[0],ileft[1],ileft[2] %d %d %d", ileft[0], ileft[1], ileft[2]);
LOGINFO("ileft_corner[0,1,2] %d %d %d", ileft_corner[0], ileft_corner[1], ileft_corner[2]);
// LOGINFO("iexpand_left[1,2,3] = (%d, %d, %d) Max %d ",iexpand_left[0],iexpand_left[1],iexpand_left[2], ileft_max_expand);
LOGINFO("iexpand_left[1,2,3] = (%d, %d, %d) Max %d ",iexpand_left[0],iexpand_left[1],iexpand_left[2],
ileft_max_expand);
LOGINFO("ileft_corner_m[0,1,2] = (%d,%d,%d)",ileft_corner_m[0],ileft_corner_m[1],ileft_corner_m[2]);
LOGINFO("grid_m_ %d grid_p_ %d",grid_m_,grid_p_);
LOGINFO("nx_m[0,1,2] = (%d,%d,%d)",nx_m[0],nx_m[1],nx_m[2]);
LOGINFO("ileft_corner_p[0,1,2] = (%d,%d,%d)",ileft_corner_p[0],ileft_corner_p[1],ileft_corner_p[2]);
LOGINFO("nxremap[0,1,2] = (%d,%d,%d)",nxremap[0],nxremap[1],nxremap[2]);
// LOGINFO("ileft_corner_m[0,1,2] = (%d,%d,%d)",ileft_corner_m[0],ileft_corner_m[1],ileft_corner_m[2]);
// LOGINFO("grid_m_ %d grid_p_ %d",grid_m_,grid_p_);
// LOGINFO("nx_m[0,1,2] = (%d,%d,%d)",nx_m[0],nx_m[1],nx_m[2]);
// LOGINFO("ileft_corner_p[0,1,2] = (%d,%d,%d)",ileft_corner_p[0],ileft_corner_p[1],ileft_corner_p[2]);
// LOGINFO("nxremap[0,1,2] = (%d,%d,%d)",nxremap[0],nxremap[1],nxremap[2]);
size_t ngp = nxremap[0] * nxremap[1] * (nxremap[2] + 2);
@ -358,7 +382,6 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R) {
double t1 = get_wtime();
double tp = t1;
#pragma omp parallel
{
#ifdef _OPENMP
@ -368,13 +391,14 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R) {
#endif
int odd_x, odd_y, odd_z;
int ng_level = ngrid_ * (1 << (level - levelmin_)); // full resolution of current level
int verbosity = (mythread == 0);
char descriptor[100];
memset(descriptor, 0, 100);
memcpy(descriptor, descriptor_string_.c_str(), descriptor_string_.size());
if (level == levelmin_) {
if (level == levelmin_)
{
start_panphasia_(&lstate[mythread], descriptor, &ng_level, &verbosity);
}
@ -388,16 +412,11 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R) {
int ratio = 1 << lextra;
assert(ng_level == ratio * d.i_base);
ix_rel[0] = ileft_corner_p[0];
ix_rel[1] = ileft_corner_p[1];
ix_rel[2] = ileft_corner_p[2];
ix_rel[0] = ileft_corner_p[0];
ix_rel[1] = ileft_corner_p[1];
ix_rel[2] = ileft_corner_p[2];
// Code above ignores the coordinate_system_shift_ - but currently this is set to zero //
// Code above ignores the coordinate_system_shift_ - but currently this is set to zero //
lstate[mythread].layer_min = 0;
lstate[mythread].layer_max = level_p;
@ -417,30 +436,34 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R) {
if (verbosity)
t1 = get_wtime();
//***************************************************************
// Process Panphasia values: p000, p001, p010, p100 and indep field
//****************************************************************
// START //
//***************************************************************
// Process Panphasia values: p000, p001, p010, p100 and indep field
//****************************************************************
// START //
#pragma omp for //nowait
for (int i = 0; i < nxremap[0] / 2 + odd_x; ++i) {
#pragma omp for // nowait
for (int i = 0; i < nxremap[0] / 2 + odd_x; ++i)
{
double cell_prop[9];
pan_state_ *ps = &lstate[mythread];
for (int j = 0; j < nxremap[1] / 2 + odd_y; ++j)
for (int k = 0; k < nxremap[2] / 2 + odd_z; ++k) {
for (int k = 0; k < nxremap[2] / 2 + odd_z; ++k)
{
// ARJ - added inner set of loops to speed up evaluation of Panphasia
for (int ix = 0; ix < 2; ++ix)
for (int iy = 0; iy < 2; ++iy)
for (int iz = 0; iz < 2; ++iz) {
for (int iz = 0; iz < 2; ++iz)
{
int ii = 2 * i + ix - odd_x;
int jj = 2 * j + iy - odd_y;
int kk = 2 * k + iz - odd_z;
if (((ii >= 0) && (ii < nxremap[0])) && ((jj >= 0) && (jj < nxremap[1])) &&
((kk >= 0) && (kk < nxremap[2]))) {
((kk >= 0) && (kk < nxremap[2])))
{
size_t idx = ((size_t)ii * nxremap[1] + (size_t)jj) * (nxremap[2] + 2) + (size_t)kk;
adv_panphasia_cell_properties_(ps, &ii, &jj, &kk, &ps->layer_min, &ps->layer_max, &ps->indep_field,
@ -456,14 +479,14 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R) {
}
}
}
LOGINFO("time for calculating PANPHASIA for level %d : %f s, %f µs/cell", level, get_wtime() - t1,
LOGUSER("time for calculating PANPHASIA for level %d : %f s, %f µs/cell", level, get_wtime() - t1,
1e6 * (get_wtime() - t1) / ((double)nxremap[2] * (double)nxremap[1] * (double)nxremap[0]));
LOGINFO("time for calculating PANPHASIA for level %d : %f s, %f µs/cell", level, get_wtime() - t1,
LOGUSER("time for calculating PANPHASIA for level %d : %f s, %f µs/cell", level, get_wtime() - t1,
1e6 * (get_wtime() - t1) / ((double)nxremap[2] * (double)nxremap[1] * (double)nxremap[0]));
//////////////////////////////////////////////////////////////////////////////////////////////
LOGINFO("\033[31mtiming level %d [adv_panphasia_cell_properties]: %f s\033[0m", level, get_wtime() - tp);
LOGUSER("\033[31mtiming level %d [adv_panphasia_cell_properties]: %f s\033[0m", level, get_wtime() - tp);
tp = get_wtime();
/////////////////////////////////////////////////////////////////////////
@ -478,7 +501,8 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R) {
#pragma omp parallel for
for (int i = 0; i < nxremap[0]; i++)
for (int j = 0; j < nxremap[1]; j++)
for (int k = 0; k < nxremap[2] / 2 + 1; k++) {
for (int k = 0; k < nxremap[2] / 2 + 1; k++)
{
size_t idx = ((size_t)i * nxremap[1] + (size_t)j) * (nxremap[2] / 2 + 1) + (size_t)k;
double fx(1.0), fy(1.0), fz(1.0), arg = 0.;
@ -492,30 +516,38 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R) {
// int kkmax = std::max(abs(ii),std::max(abs(jj),abs(kk)));
if (ii != 0) {
if (ii != 0)
{
arg = M_PI * (double)ii / (double)nxremap[0];
fx = sin(arg) / arg;
gx = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
} else {
}
else
{
fx = 1.0;
gx = 0.0;
}
if (jj != 0) {
if (jj != 0)
{
arg = M_PI * (double)jj / (double)nxremap[1];
fy = sin(arg) / arg;
gy = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
} else {
}
else
{
fy = 1.0;
gy = 0.0;
}
if (kk != 0) {
if (kk != 0)
{
arg = M_PI * (double)kk / (double)nxremap[2];
fz = sin(arg) / arg;
gz = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
} else {
}
else
{
fz = 1.0;
gz = 0.0;
}
@ -524,7 +556,8 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R) {
double magnitude = sqrt(1.0 - std::abs(temp_comp * temp_comp));
if (abs(ii) != nxremap[0] / 2 && abs(jj) != nxremap[1] / 2 &&
abs(kk) != nxremap[2] / 2) { // kkmax != nxremap[2]/2 ){
abs(kk) != nxremap[2] / 2)
{ // kkmax != nxremap[2]/2 ){
complex x, y0(RE(pc0[idx]), IM(pc0[idx])), y1(RE(pc1[idx]), IM(pc1[idx])), y2(RE(pc2[idx]), IM(pc2[idx])),
y3(RE(pc3[idx]), IM(pc3[idx])), y4(RE(pc4[idx]), IM(pc4[idx]));
@ -538,12 +571,12 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R) {
// END
LOGINFO("\033[31mtiming level %d [build panphasia field]: %f s\033[0m", level, get_wtime() - tp);
LOGUSER("\033[31mtiming level %d [build panphasia field]: %f s\033[0m", level, get_wtime() - tp);
tp = get_wtime();
//***************************************************************
// Process Panphasia values: p000, p001, p010, p100 and indep field
//****************************************************************
//***************************************************************
// Process Panphasia values: p000, p001, p010, p100 and indep field
//****************************************************************
#pragma omp parallel
{
@ -559,7 +592,8 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R) {
memset(descriptor, 0, 100);
memcpy(descriptor, descriptor_string_.c_str(), descriptor_string_.size());
if (level == levelmin_) {
if (level == levelmin_)
{
start_panphasia_(&lstate[mythread], descriptor, &ng_level, &verbosity);
}
@ -577,7 +611,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R) {
ix_rel[1] = ileft_corner_p[1];
ix_rel[2] = ileft_corner_p[2];
// Code above ignores the coordinate_system_shift_ - but currently this is set to zero //
// Code above ignores the coordinate_system_shift_ - but currently this is set to zero //
lstate[mythread].layer_min = 0;
lstate[mythread].layer_max = level_p;
@ -601,25 +635,29 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R) {
//***************************************************************
// Process Panphasia values: p110, p011, p101, p111
//****************************************************************
#pragma omp for //nowait
for (int i = 0; i < nxremap[0] / 2 + odd_x; ++i) {
#pragma omp for // nowait
for (int i = 0; i < nxremap[0] / 2 + odd_x; ++i)
{
double cell_prop[9];
pan_state_ *ps = &lstate[mythread];
for (int j = 0; j < nxremap[1] / 2 + odd_y; ++j)
for (int k = 0; k < nxremap[2] / 2 + odd_z; ++k) {
for (int k = 0; k < nxremap[2] / 2 + odd_z; ++k)
{
// ARJ - added inner set of loops to speed up evaluation of Panphasia
for (int ix = 0; ix < 2; ++ix)
for (int iy = 0; iy < 2; ++iy)
for (int iz = 0; iz < 2; ++iz) {
for (int iz = 0; iz < 2; ++iz)
{
int ii = 2 * i + ix - odd_x;
int jj = 2 * j + iy - odd_y;
int kk = 2 * k + iz - odd_z;
if (((ii >= 0) && (ii < nxremap[0])) && ((jj >= 0) && (jj < nxremap[1])) &&
((kk >= 0) && (kk < nxremap[2]))) {
((kk >= 0) && (kk < nxremap[2])))
{
size_t idx = ((size_t)ii * nxremap[1] + (size_t)jj) * (nxremap[2] + 2) + (size_t)kk;
adv_panphasia_cell_properties_(ps, &ii, &jj, &kk, &ps->layer_min, &ps->layer_max, &ps->indep_field,
@ -637,7 +675,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R) {
LOGINFO("time for calculating PANPHASIA for level %d : %f s, %f µs/cell", level, get_wtime() - t1,
1e6 * (get_wtime() - t1) / ((double)nxremap[2] * (double)nxremap[1] * (double)nxremap[0]));
LOGINFO("\033[31mtiming level %d [adv_panphasia_cell_properties2]: %f s \033[0m", level, get_wtime() - tp);
LOGUSER("\033[31mtiming level %d [adv_panphasia_cell_properties2]: %f s \033[0m", level, get_wtime() - tp);
tp = get_wtime();
/////////////////////////////////////////////////////////////////////////
@ -651,7 +689,8 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R) {
#pragma omp parallel for
for (int i = 0; i < nxremap[0]; i++)
for (int j = 0; j < nxremap[1]; j++)
for (int k = 0; k < nxremap[2] / 2 + 1; k++) {
for (int k = 0; k < nxremap[2] / 2 + 1; k++)
{
size_t idx = ((size_t)i * nxremap[1] + (size_t)j) * (nxremap[2] / 2 + 1) + (size_t)k;
double fx(1.0), fy(1.0), fz(1.0), arg = 0.;
@ -665,35 +704,45 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R) {
// int kkmax = std::max(abs(ii),std::max(abs(jj),abs(kk)));
if (ii != 0) {
if (ii != 0)
{
arg = M_PI * (double)ii / (double)nxremap[0];
fx = sin(arg) / arg;
gx = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
} else {
}
else
{
fx = 1.0;
gx = 0.0;
}
if (jj != 0) {
if (jj != 0)
{
arg = M_PI * (double)jj / (double)nxremap[1];
fy = sin(arg) / arg;
gy = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
} else {
}
else
{
fy = 1.0;
gy = 0.0;
}
if (kk != 0) {
if (kk != 0)
{
arg = M_PI * (double)kk / (double)nxremap[2];
fz = sin(arg) / arg;
gz = complex(0.0, (arg * cos(arg) - sin(arg)) / (arg * arg));
} else {
}
else
{
fz = 1.0;
gz = 0.0;
}
if (abs(ii) != nxremap[0] / 2 && abs(jj) != nxremap[1] / 2 &&
abs(kk) != nxremap[2] / 2) { // kkmax != nxremap[2]/2 ){
abs(kk) != nxremap[2] / 2)
{ // kkmax != nxremap[2]/2 ){
complex x, y1(RE(pc1[idx]), IM(pc1[idx])), y2(RE(pc2[idx]), IM(pc2[idx])), y3(RE(pc3[idx]), IM(pc3[idx])),
y4(RE(pc4[idx]), IM(pc4[idx]));
@ -704,7 +753,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R) {
}
}
LOGINFO("\033[31mtiming level %d [build panphasia field2]: %f s\033[0m", level, get_wtime() - tp);
LOGUSER("\033[31mtiming level %d [build panphasia field2]: %f s\033[0m", level, get_wtime() - tp);
tp = get_wtime();
// END
@ -717,21 +766,24 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R) {
// do we need to cut off the small scales?
// int nn = 1<<level;
if (incongruent_fields_) {
if (incongruent_fields_)
{
LOGINFO("Remapping fields from dimension %d -> %d",nxremap[0],nx_m[0]);
LOGINFO("Remapping fields from dimension %d -> %d", nxremap[0], nx_m[0]);
memset(pr1, 0, ngp * sizeof(fftw_real));
#pragma omp parallel for
#pragma omp parallel for
for (int i = 0; i < nxremap[0]; i++)
for (int j = 0; j < nxremap[1]; j++)
for (int k = 0; k < nxremap[2] / 2 + 1; k++) {
for (int k = 0; k < nxremap[2] / 2 + 1; k++)
{
int ii = (i > nxremap[0] / 2) ? i - nxremap[0] : i, jj = (j > nxremap[1] / 2) ? j - nxremap[1] : j, kk = k;
int ia(abs(ii)), ja(abs(jj)), ka(abs(kk));
if (ia < nx_m[0] / 2 && ja < nx_m[1] / 2 && ka < nx_m[2] / 2) {
if (ia < nx_m[0] / 2 && ja < nx_m[1] / 2 && ka < nx_m[2] / 2)
{
size_t idx = ((size_t)(i)*nxremap[1] + (size_t)(j)) * (nxremap[2] / 2 + 1) + (size_t)(k);
@ -739,27 +791,26 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R) {
size_t idx2 = ((size_t)ir * nx_m[1] + (size_t)jr) * ((size_t)nx_m[2] / 2 + 1) + (size_t)kr;
complex x(RE(pc0[idx]),IM(pc0[idx]));
double total_phase_shift;
total_phase_shift = inter_grid_phase_adjustment_ * (double)(ii+jj+kk);
x = x * exp(complex(0.0, total_phase_shift));
RE(pc1[idx2]) = x.real();
IM(pc1[idx2]) = x.imag();
}
complex x(RE(pc0[idx]), IM(pc0[idx]));
double total_phase_shift;
total_phase_shift = inter_grid_phase_adjustment_ * (double)(ii + jj + kk);
x = x * exp(complex(0.0, total_phase_shift));
RE(pc1[idx2]) = x.real();
IM(pc1[idx2]) = x.imag();
}
}
memcpy(pr0, pr1, ngp * sizeof(fftw_real));
}
// if (level == 9)
// {
// LOGUSER("DC mode of level is %g", RE(pc0[0]));
// // RE(pc0[0]) = 1e8;
// // IM(pc0[0]) = 0.0;
// }
if( level == 9 ){
LOGINFO("DC mode of level is %g",RE(pc0[0]));
//RE(pc0[0]) = 1e8;
//IM(pc0[0]) = 0.0;
}
LOGINFO("\033[31mtiming level %d [remap noncongruent]: %f s\033[0m", level, get_wtime() - tp);
LOGUSER("\033[31mtiming level %d [remap noncongruent]: %f s\033[0m", level, get_wtime() - tp);
tp = get_wtime();
/////////////////////////////////////////////////////////////////////////
// transform back
@ -774,7 +825,7 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R) {
delete[] pr4;
LOGINFO("Copying random field data %d,%d,%d -> %d,%d,%d", nxremap[0], nxremap[1], nxremap[2], nx[0], nx[1], nx[2]);
// n = 1<<level;
// ng = n;
// ngp = ng*ng*2*(ng/2+1);
@ -785,15 +836,16 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R) {
/*double norm = 1.0 / sqrt((double)nxremap[0] * (double)nxremap[1] * (double)nxremap[2] * (double)nx[0] *
(double)nx[1] * (double)nx[2]);*/
double norm = 1.0 / sqrt((double)nxremap[0] * (double)nxremap[1] * (double)nxremap[2] * (double)nx_m[0] *
(double)nx_m[1] * (double)nx_m[2]);
double norm = 1.0 / sqrt((double)nxremap[0] * (double)nxremap[1] * (double)nxremap[2] * (double)nx_m[0] *
(double)nx_m[1] * (double)nx_m[2]);
#pragma omp parallel for reduction(+ : sum, sum2, count)
#pragma omp parallel for reduction(+ \
: sum, sum2, count)
for (int k = 0; k < nx[2]; ++k) // ARJ - swapped roles of i,k, and reverse ordered loops
for (int j = 0; j < nx[1]; ++j)
for (int i = 0; i < nx[0]; ++i) {
size_t idx = ((size_t)(i+iexpand_left[0])*nx_m[1] + (size_t)(j+iexpand_left[1])) * (nx_m[2] + 2)
+ (size_t)(k+iexpand_left[2]);
for (int i = 0; i < nx[0]; ++i)
{
size_t idx = ((size_t)(i + iexpand_left[0]) * nx_m[1] + (size_t)(j + iexpand_left[1])) * (nx_m[2] + 2) + (size_t)(k + iexpand_left[2]);
R(i, j, k) = pr0[idx] * norm;
sum += R(i, j, k);
@ -814,8 +866,9 @@ void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R) {
LOGINFO("PANPHASIA level %d mean and variance are\n <p> = %g | var(p) = %g", level, sum, sum2);
}
namespace {
RNG_plugin_creator_concrete<RNG_panphasia> creator("PANPHASIA");
namespace
{
RNG_plugin_creator_concrete<RNG_panphasia> creator("PANPHASIA");
}
#endif

1837
src/random.cc
View file

File diff suppressed because it is too large Load diff

371
src/random.hh
View file

@ -4,7 +4,7 @@
a code to generate multi-scale initial conditions
for cosmological simulations
Copyright (C) 2010 Oliver Hahn
Copyright (C) 2010-23 by Oliver Hahn
*/
@ -13,15 +13,17 @@
// #define DEGRADE_RAND2
//.....................................
#ifndef __RANDOM_HH
#define __RANDOM_HH
#pragma once
#define DEF_RAN_CUBE_SIZE 32
#include <fstream>
#include <algorithm>
#include <map>
#ifdef _OPENMP
#include <omp.h>
#endif
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
@ -30,18 +32,23 @@
#include "mesh.hh"
#include "mg_operators.hh"
#include "constraints.hh"
// #include "convolution_kernel.hh"
#include "density_grid.hh"
class RNG_plugin
{
protected:
config_file *pcf_; //!< pointer to config_file from which to read parameters
config_file *pcf_; //!< pointer to config_file from which to read parameters
const refinement_hierarchy *prefh_; //!< pointer to refinement hierarchy structure containing the grid sizes
public:
explicit RNG_plugin(config_file &cf)
: pcf_(&cf)
explicit RNG_plugin(config_file &cf) //, const refinement_hierarchy& refh )
: pcf_(&cf) //, prefh_( & refh )
{
}
virtual ~RNG_plugin() {}
virtual bool is_multiscale() const = 0;
virtual void fill_grid(int level, DensityGrid<real_t> &R) = 0;
virtual void initialize_for_grid_structure(const refinement_hierarchy &refh) = 0;
};
struct RNG_plugin_creator
@ -65,356 +72,54 @@ struct RNG_plugin_creator_concrete : public RNG_plugin_creator
}
//! create an instance of the plugin
RNG_plugin *create(config_file &cf) const
RNG_plugin *create(config_file &cf) const //, const refinement_hierarchy& refh ) const
{
return new Derived(cf);
return new Derived(cf); //, refh );
}
};
typedef RNG_plugin RNG_instance;
RNG_plugin *select_RNG_plugin(config_file &cf);
/*!
* @brief encapsulates all things random number generator related
*/
template <typename T>
class random_numbers
{
public:
unsigned
res_, //!< resolution of the full mesh
cubesize_, //!< size of one independent random number cube
ncubes_; //!< number of random number cubes to cover the full mesh
long baseseed_; //!< base seed from which cube seeds are computed
protected:
//! vector of 3D meshes (the random number cubes) with random numbers
std::vector<Meshvar<T> *> rnums_;
//! map of 3D indices to cube index
std::map<size_t, size_t> cubemap_;
typedef std::map<size_t, size_t>::iterator cubemap_iterator;
protected:
//! register a cube with the hash map
void register_cube(int i, int j, int k);
//! fills a subcube with random numbers
double fill_cube(int i, int j, int k);
//! subtract a constant from an entire cube
void subtract_from_cube(int i, int j, int k, double val);
//! copy random numbers from a cube to a full grid array
template <class C>
void copy_cube(int i, int j, int k, C &dat)
{
int offi, offj, offk;
offi = i * cubesize_;
offj = j * cubesize_;
offk = k * cubesize_;
i = (i + ncubes_) % ncubes_;
j = (j + ncubes_) % ncubes_;
k = (k + ncubes_) % ncubes_;
size_t icube = (i * ncubes_ + j) * ncubes_ + k;
cubemap_iterator it = cubemap_.find(icube);
if (it == cubemap_.end())
{
LOGERR("attempting to copy data from non-existing RND cube %d,%d,%d", i, j, k);
throw std::runtime_error("attempting to copy data from non-existing RND cube");
}
size_t cubeidx = it->second;
for (int ii = 0; ii < (int)cubesize_; ++ii)
for (int jj = 0; jj < (int)cubesize_; ++jj)
for (int kk = 0; kk < (int)cubesize_; ++kk)
dat(offi + ii, offj + jj, offk + kk) = (*rnums_[cubeidx])(ii, jj, kk);
}
//! free the memory associated with a subcube
void free_cube(int i, int j, int k);
//! initialize member variables and allocate memory
void initialize(void);
//! fill a cubic subvolume of the full grid with random numbers
double fill_subvolume(int *i0, int *n);
//! fill an entire grid with random numbers
double fill_all(void);
//! fill an external array instead of the internal field
template <class C>
double fill_all(C &dat)
{
double sum = 0.0;
for (int i = 0; i < (int)ncubes_; ++i)
for (int j = 0; j < (int)ncubes_; ++j)
for (int k = 0; k < (int)ncubes_; ++k)
{
int ii(i), jj(j), kk(k);
register_cube(ii, jj, kk);
}
#pragma omp parallel for reduction(+ \
: sum)
for (int i = 0; i < (int)ncubes_; ++i)
for (int j = 0; j < (int)ncubes_; ++j)
for (int k = 0; k < (int)ncubes_; ++k)
{
int ii(i), jj(j), kk(k);
ii = (ii + ncubes_) % ncubes_;
jj = (jj + ncubes_) % ncubes_;
kk = (kk + ncubes_) % ncubes_;
sum += fill_cube(ii, jj, kk);
copy_cube(ii, jj, kk, dat);
free_cube(ii, jj, kk);
}
return sum / (ncubes_ * ncubes_ * ncubes_);
}
//! write the number of allocated random number cubes to stdout
void print_allocated(void);
public:
//! constructor
random_numbers(unsigned res, unsigned cubesize, long baseseed, int *x0, int *lx);
//! constructor for constrained fine field
random_numbers(random_numbers<T> &rc, unsigned cubesize, long baseseed,
bool kspace = false, bool isolated = false, int *x0_ = NULL, int *lx_ = NULL, bool zeromean = true);
//! constructor
random_numbers(unsigned res, unsigned cubesize, long baseseed, bool zeromean = true);
//! constructor to read white noise from file
random_numbers(unsigned res, std::string randfname, bool rndsign);
//! copy constructor for averaged field (not copying) hence explicit!
explicit random_numbers(/*const*/ random_numbers<T> &rc, bool kdegrade = true);
//! destructor
~random_numbers()
{
for (unsigned i = 0; i < rnums_.size(); ++i)
if (rnums_[i] != NULL)
delete rnums_[i];
rnums_.clear();
}
//! access a random number, this allocates a cube and fills it with consistent random numbers
inline T &operator()(int i, int j, int k, bool fillrand = true)
{
int ic, jc, kc, is, js, ks;
if (ncubes_ == 0)
throw std::runtime_error("random_numbers: internal error, not properly initialized");
//... determine cube
ic = (int)((double)i / cubesize_ + ncubes_) % ncubes_;
jc = (int)((double)j / cubesize_ + ncubes_) % ncubes_;
kc = (int)((double)k / cubesize_ + ncubes_) % ncubes_;
size_t icube = ((size_t)ic * ncubes_ + (size_t)jc) * ncubes_ + (size_t)kc;
cubemap_iterator it = cubemap_.find(icube);
if (it == cubemap_.end())
{
LOGERR("Attempting to copy data from non-existing RND cube %d,%d,%d @ %d,%d,%d", ic, jc, kc, i, j, k);
throw std::runtime_error("attempting to copy data from non-existing RND cube");
}
size_t cubeidx = it->second;
if (rnums_[cubeidx] == NULL)
{
LOGERR("Attempting to access data from non-allocated RND cube %d,%d,%d", ic, jc, kc);
throw std::runtime_error("attempting to access data from non-allocated RND cube");
}
//... determine cell in cube
is = (i - ic * cubesize_ + cubesize_) % cubesize_;
js = (j - jc * cubesize_ + cubesize_) % cubesize_;
ks = (k - kc * cubesize_ + cubesize_) % cubesize_;
return (*rnums_[cubeidx])(is, js, ks);
}
//! free all cubes
void free_all_mem(void)
{
for (unsigned i = 0; i < rnums_.size(); ++i)
if (rnums_[i] != NULL)
{
delete rnums_[i];
rnums_[i] = NULL;
}
}
};
RNG_plugin *select_RNG_plugin(config_file &cf); //, const refinement_hierarchy& refh );
/*!
* @brief encapsulates all things for multi-scale white noise generation
*/
template <typename rng, typename T>
template <typename T>
class random_number_generator
{
protected:
config_file *pcf_;
refinement_hierarchy *prefh_;
constraint_set constraints;
int levelmin_,
levelmax_,
levelmin_seed_;
std::vector<long> rngseeds_;
std::vector<std::string> rngfnames_;
bool disk_cached_;
bool restart_;
std::vector<std::vector<T> *> mem_cache_;
unsigned ran_cube_size_;
protected:
//! checks if the specified string is numeric
bool is_number(const std::string &s);
//! parses the random number parameters in the conf file
void parse_rand_parameters(void);
//! correct coarse grid averages for the change in small scale when using Fourier interpolation
void correct_avg(int icoarse, int ifine);
//! the main driver routine for multi-scale white noise generation
void compute_random_numbers(void);
//! store the white noise fields in memory or on disk
void store_rnd(int ilevel, rng *prng);
// const refinement_hierarchy * prefh_;
RNG_plugin *generator_;
int levelmin_, levelmax_;
public:
//! constructor
random_number_generator(config_file &cf, refinement_hierarchy &refh, transfer_function *ptf = NULL);
random_number_generator(config_file &cf, transfer_function *ptf = NULL)
: pcf_(&cf) //, prefh_( &refh )
{
levelmin_ = pcf_->getValue<int>("setup", "levelmin");
levelmax_ = pcf_->getValue<int>("setup", "levelmax");
generator_ = select_RNG_plugin(cf);
}
//! destructor
~random_number_generator();
~random_number_generator()
{
}
//! initialize_for_grid_structure
void initialize_for_grid_structure(const refinement_hierarchy &refh)
{
generator_->initialize_for_grid_structure(refh);
}
//! load random numbers to a new array
template <typename array>
void load(array &A, int ilevel)
{
if (restart_)
LOGINFO("Attempting to restart using random numbers for level %d\n from file \'wnoise_%04d.bin\'.", ilevel, ilevel);
if (disk_cached_)
{
char fname[128];
sprintf(fname, "wnoise_%04d.bin", ilevel);
LOGUSER("Loading white noise from file \'%s\'...", fname);
std::ifstream ifs(fname, std::ios::binary);
if (!ifs.good())
{
LOGERR("White noise file \'%s\'was not found.", fname);
throw std::runtime_error("A white noise file was not found. This is an internal inconsistency and bad.");
}
int nx, ny, nz;
ifs.read(reinterpret_cast<char *>(&nx), sizeof(int));
ifs.read(reinterpret_cast<char *>(&ny), sizeof(int));
ifs.read(reinterpret_cast<char *>(&nz), sizeof(int));
if (nx != (int)A.size(0) || ny != (int)A.size(1) || nz != (int)A.size(2))
{
if (nx == (int)A.size(0) * 2 && ny == (int)A.size(1) * 2 && nz == (int)A.size(2) * 2)
{
std::cerr << "CHECKPOINT" << std::endl;
int ox = nx / 4, oy = ny / 4, oz = nz / 4;
std::vector<T> slice(ny * nz, 0.0);
for (int i = 0; i < nx; ++i)
{
ifs.read(reinterpret_cast<char *>(&slice[0]), ny * nz * sizeof(T));
if (i < ox)
continue;
if (i >= 3 * ox)
break;
#pragma omp parallel for
for (int j = oy; j < 3 * oy; ++j)
for (int k = oz; k < 3 * oz; ++k)
A(i - ox, j - oy, k - oz) = slice[j * nz + k];
}
ifs.close();
}
else
{
LOGERR("White noise file is not aligned with array. File: [%d,%d,%d]. Mem: [%d,%d,%d].",
nx, ny, nz, A.size(0), A.size(1), A.size(2));
throw std::runtime_error("White noise file is not aligned with array. This is an internal inconsistency and bad.");
}
}
else
{
for (int i = 0; i < nx; ++i)
{
std::vector<T> slice(ny * nz, 0.0);
ifs.read(reinterpret_cast<char *>(&slice[0]), ny * nz * sizeof(T));
#pragma omp parallel for
for (int j = 0; j < ny; ++j)
for (int k = 0; k < nz; ++k)
A(i, j, k) = slice[j * nz + k];
}
ifs.close();
}
}
else
{
LOGUSER("Copying white noise from memory cache...");
if (mem_cache_[ilevel - levelmin_] == NULL)
LOGERR("Tried to access mem-cached random numbers for level %d. But these are not available!\n", ilevel);
int nx(A.size(0)), ny(A.size(1)), nz(A.size(2));
if ((size_t)nx * (size_t)ny * (size_t)nz != mem_cache_[ilevel - levelmin_]->size())
{
LOGERR("White noise file is not aligned with array. File: [%d,%d,%d]. Mem: [%d,%d,%d].", nx, ny, nz, A.size(0), A.size(1), A.size(2));
throw std::runtime_error("White noise file is not aligned with array. This is an internal inconsistency and bad");
}
#pragma omp parallel for
for (int i = 0; i < nx; ++i)
for (int j = 0; j < ny; ++j)
for (int k = 0; k < nz; ++k)
A(i, j, k) = (*mem_cache_[ilevel - levelmin_])[((size_t)i * ny + (size_t)j) * nz + (size_t)k];
std::vector<T>().swap(*mem_cache_[ilevel - levelmin_]);
delete mem_cache_[ilevel - levelmin_];
mem_cache_[ilevel - levelmin_] = NULL;
}
generator_->fill_grid(ilevel, A);
}
};
typedef random_numbers<real_t> rand_nums;
typedef random_number_generator<rand_nums, real_t> rand_gen;
#endif //__RANDOM_HH
using noise_generator = random_number_generator<real_t>;