mirror of
https://github.com/cosmo-sims/MUSIC.git
synced 2024-09-19 17:03:46 +02:00
wip commit, memory optimization
This commit is contained in:
parent
33d784a137
commit
437e77f74f
6 changed files with 364 additions and 341 deletions
|
@ -298,9 +298,15 @@ public:
|
||||||
}
|
}
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
|
if( prefh_->get_margin() < 0 ){
|
||||||
cparam_.nx = 2 * prefh_->size(ilevel, 0);
|
cparam_.nx = 2 * prefh_->size(ilevel, 0);
|
||||||
cparam_.ny = 2 * prefh_->size(ilevel, 1);
|
cparam_.ny = 2 * prefh_->size(ilevel, 1);
|
||||||
cparam_.nz = 2 * prefh_->size(ilevel, 2);
|
cparam_.nz = 2 * prefh_->size(ilevel, 2);
|
||||||
|
}else{
|
||||||
|
cparam_.nx = prefh_->size(ilevel, 0) + 2*prefh_->get_margin();
|
||||||
|
cparam_.ny = prefh_->size(ilevel, 1) + 2*prefh_->get_margin();
|
||||||
|
cparam_.nz = prefh_->size(ilevel, 2) + 2*prefh_->get_margin();
|
||||||
|
}
|
||||||
|
|
||||||
cparam_.lx = (double)cparam_.nx / (double)(1 << ilevel) * boxlength_;
|
cparam_.lx = (double)cparam_.nx / (double)(1 << ilevel) * boxlength_;
|
||||||
cparam_.ly = (double)cparam_.ny / (double)(1 << ilevel) * boxlength_;
|
cparam_.ly = (double)cparam_.ny / (double)(1 << ilevel) * boxlength_;
|
||||||
|
|
|
@ -632,12 +632,22 @@ void GenerateDensityHierarchy(config_file &cf, transfer_function *ptf, tf_type t
|
||||||
LOGUSER(" size =(%5d,%5d,%5d)", refh.size(levelmin + i, 0),
|
LOGUSER(" size =(%5d,%5d,%5d)", refh.size(levelmin + i, 0),
|
||||||
refh.size(levelmin + i, 1), refh.size(levelmin + i, 2));
|
refh.size(levelmin + i, 1), refh.size(levelmin + i, 2));
|
||||||
|
|
||||||
|
if( refh.get_margin() > 0 )
|
||||||
|
fine = new PaddedDensitySubGrid<real_t>(refh.offset(levelmin + i, 0),
|
||||||
|
refh.offset(levelmin + i, 1),
|
||||||
|
refh.offset(levelmin + i, 2),
|
||||||
|
refh.size(levelmin + i, 0),
|
||||||
|
refh.size(levelmin + i, 1),
|
||||||
|
refh.size(levelmin + i, 2),
|
||||||
|
refh.get_margin(), refh.get_margin(), refh.get_margin() );
|
||||||
|
else
|
||||||
fine = new PaddedDensitySubGrid<real_t>(refh.offset(levelmin + i, 0),
|
fine = new PaddedDensitySubGrid<real_t>(refh.offset(levelmin + i, 0),
|
||||||
refh.offset(levelmin + i, 1),
|
refh.offset(levelmin + i, 1),
|
||||||
refh.offset(levelmin + i, 2),
|
refh.offset(levelmin + i, 2),
|
||||||
refh.size(levelmin + i, 0),
|
refh.size(levelmin + i, 0),
|
||||||
refh.size(levelmin + i, 1),
|
refh.size(levelmin + i, 1),
|
||||||
refh.size(levelmin + i, 2));
|
refh.size(levelmin + i, 2));
|
||||||
|
|
||||||
/////////////////////////////////////////////////////////////////////////
|
/////////////////////////////////////////////////////////////////////////
|
||||||
|
|
||||||
// load white noise for patch
|
// load white noise for patch
|
||||||
|
|
|
@ -230,7 +230,7 @@ void store_grid_structure( config_file& cf, const refinement_hierarchy& rh )
|
||||||
for( int j=0; j<3; ++j )
|
for( int j=0; j<3; ++j )
|
||||||
{
|
{
|
||||||
sprintf(str1,"offset(%d,%d)",i,j);
|
sprintf(str1,"offset(%d,%d)",i,j);
|
||||||
sprintf(str2,"%d",rh.offset(i,j));
|
sprintf(str2,"%ld",rh.offset(i,j));
|
||||||
cf.insertValue("setup",str1,str2);
|
cf.insertValue("setup",str1,str2);
|
||||||
|
|
||||||
sprintf(str1,"size(%d,%d)",i,j);
|
sprintf(str1,"size(%d,%d)",i,j);
|
||||||
|
|
35
src/mesh.hh
35
src/mesh.hh
|
@ -1207,17 +1207,6 @@ class refinement_hierarchy
|
||||||
std::vector<index3_t> absoffsets_;
|
std::vector<index3_t> absoffsets_;
|
||||||
std::vector<index3_t> len_;
|
std::vector<index3_t> len_;
|
||||||
|
|
||||||
// std::vector<unsigned>
|
|
||||||
// ox_, //!< relative x-coordinates of grid origins (in coarser grid cells)
|
|
||||||
// oy_, //!< relative y-coordinates of grid origins (in coarser grid cells)
|
|
||||||
// oz_, //!< relative z-coordinates of grid origins (in coarser grid cells)
|
|
||||||
// oax_, //!< absolute x-coordinates of grid origins (in fine grid cells)
|
|
||||||
// oay_, //!< absolute y-coordinates of grid origins (in fine grid cells)
|
|
||||||
// oaz_, //!< absolute z-coordinates of grid origins (in fine grid cells)
|
|
||||||
// nx_, //!< x-extent of grids (in fine grid cells)
|
|
||||||
// ny_, //!< y-extent of grids (in fine grid cells)
|
|
||||||
// nz_; //!< z-extent of grids (in fine grid cells)
|
|
||||||
|
|
||||||
unsigned
|
unsigned
|
||||||
levelmin_, //!< minimum grid level for Poisson solver
|
levelmin_, //!< minimum grid level for Poisson solver
|
||||||
levelmax_, //!< maximum grid level for all operations
|
levelmax_, //!< maximum grid level for all operations
|
||||||
|
@ -1225,16 +1214,14 @@ class refinement_hierarchy
|
||||||
padding_, //!< padding in number of coarse cells between refinement levels
|
padding_, //!< padding in number of coarse cells between refinement levels
|
||||||
blocking_factor_;
|
blocking_factor_;
|
||||||
|
|
||||||
|
int margin_; //!< number of cells used for additional padding for convolutions with isolated boundaries (-1 = double padding)
|
||||||
|
|
||||||
config_file &cf_; //!< reference to config_file
|
config_file &cf_; //!< reference to config_file
|
||||||
|
|
||||||
bool align_top_, //!< bool whether to align all grids with coarsest grid cells
|
bool align_top_, //!< bool whether to align all grids with coarsest grid cells
|
||||||
preserve_dims_, //!< bool whether to preserve user-specified grid dimensions
|
preserve_dims_, //!< bool whether to preserve user-specified grid dimensions
|
||||||
equal_extent_; //!< bool whether the simulation code requires squared refinement regions (e.g. RAMSES)
|
equal_extent_; //!< bool whether the simulation code requires squared refinement regions (e.g. RAMSES)
|
||||||
|
|
||||||
// double
|
|
||||||
// x0ref_[3], //!< coordinates of refinement region origin (in [0..1[)
|
|
||||||
// lxref_[3]; //!< extent of refinement region (int [0..1[)
|
|
||||||
|
|
||||||
vec3_t x0ref_; //!< coordinates of refinement region origin (in [0..1[)
|
vec3_t x0ref_; //!< coordinates of refinement region origin (in [0..1[)
|
||||||
vec3_t lxref_; //!< extent of refinement region (int [0..1[)
|
vec3_t lxref_; //!< extent of refinement region (int [0..1[)
|
||||||
|
|
||||||
|
@ -1264,6 +1251,7 @@ public:
|
||||||
preserve_dims_ = cf_.getValueSafe<bool>("setup", "preserve_dims", false);
|
preserve_dims_ = cf_.getValueSafe<bool>("setup", "preserve_dims", false);
|
||||||
equal_extent_ = cf_.getValueSafe<bool>("setup", "force_equal_extent", false);
|
equal_extent_ = cf_.getValueSafe<bool>("setup", "force_equal_extent", false);
|
||||||
blocking_factor_ = cf.getValueSafe<unsigned>("setup", "blocking_factor", 0);
|
blocking_factor_ = cf.getValueSafe<unsigned>("setup", "blocking_factor", 0);
|
||||||
|
margin_ = cf.getValueSafe<int>("setup","convolution_margin",32);
|
||||||
|
|
||||||
bool bnoshift = cf_.getValueSafe<bool>("setup", "no_shift", false);
|
bool bnoshift = cf_.getValueSafe<bool>("setup", "no_shift", false);
|
||||||
bool force_shift = cf_.getValueSafe<bool>("setup", "force_shift", false);
|
bool force_shift = cf_.getValueSafe<bool>("setup", "force_shift", false);
|
||||||
|
@ -1289,14 +1277,8 @@ public:
|
||||||
// if not doing any refinement levels, set extent to full box
|
// if not doing any refinement levels, set extent to full box
|
||||||
if (levelmin_ == levelmax_)
|
if (levelmin_ == levelmax_)
|
||||||
{
|
{
|
||||||
x0ref_ = {0.0, 0.0, 0.0};
|
x0ref_ = { 0.0, 0.0, 0.0 };
|
||||||
// x0ref_[0] = 0.0;
|
lxref_ = { 1.0, 1.0, 1.0 };
|
||||||
// x0ref_[1] = 0.0;
|
|
||||||
// x0ref_[2] = 0.0;
|
|
||||||
|
|
||||||
lxref_[0] = 1.0;
|
|
||||||
lxref_[1] = 1.0;
|
|
||||||
lxref_[2] = 1.0;
|
|
||||||
}
|
}
|
||||||
|
|
||||||
unsigned ncoarse = 1 << levelmin_;
|
unsigned ncoarse = 1 << levelmin_;
|
||||||
|
@ -1658,6 +1640,7 @@ public:
|
||||||
offsets_ = o.offsets_;
|
offsets_ = o.offsets_;
|
||||||
absoffsets_ = o.absoffsets_;
|
absoffsets_ = o.absoffsets_;
|
||||||
len_ = o.len_;
|
len_ = o.len_;
|
||||||
|
margin_ = o.margin_;
|
||||||
// ox_ = o.ox_;
|
// ox_ = o.ox_;
|
||||||
// oy_ = o.oy_;
|
// oy_ = o.oy_;
|
||||||
// oz_ = o.oz_;
|
// oz_ = o.oz_;
|
||||||
|
@ -1776,6 +1759,12 @@ public:
|
||||||
return xshift_[idim];
|
return xshift_[idim];
|
||||||
}
|
}
|
||||||
|
|
||||||
|
//! get the margin reserved for isolated convolutions (-1=double pad)
|
||||||
|
int get_margin() const
|
||||||
|
{
|
||||||
|
return margin_;
|
||||||
|
}
|
||||||
|
|
||||||
//! get the total shift of the coordinate system in box coordinates
|
//! get the total shift of the coordinate system in box coordinates
|
||||||
const double *get_coord_shift(void) const
|
const double *get_coord_shift(void) const
|
||||||
{
|
{
|
||||||
|
|
|
@ -1697,12 +1697,26 @@ void random_number_generator<rng, T>::compute_random_numbers(void)
|
||||||
|
|
||||||
int lfac = 1 << (ilevel - levelmin_poisson);
|
int lfac = 1 << (ilevel - levelmin_poisson);
|
||||||
|
|
||||||
lx[0] = 2 * prefh_->size(ilevel, 0);
|
int margin[3];
|
||||||
lx[1] = 2 * prefh_->size(ilevel, 1);
|
if (prefh_->get_margin()>0){
|
||||||
lx[2] = 2 * prefh_->size(ilevel, 2);
|
margin[0] = prefh_->get_margin();
|
||||||
x0[0] = prefh_->offset_abs(ilevel, 0) - lfac * shift[0] - lx[0] / 4;
|
margin[1] = prefh_->get_margin();
|
||||||
x0[1] = prefh_->offset_abs(ilevel, 1) - lfac * shift[1] - lx[1] / 4;
|
margin[2] = prefh_->get_margin();
|
||||||
x0[2] = prefh_->offset_abs(ilevel, 2) - lfac * shift[2] - lx[2] / 4;
|
}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];//lx[0] / 4;
|
||||||
|
x0[1] = prefh_->offset_abs(ilevel, 1) - lfac * shift[1] - margin[1];//lx[1] / 4;
|
||||||
|
x0[2] = prefh_->offset_abs(ilevel, 2) - lfac * shift[2] - margin[2];//lx[2] / 4;
|
||||||
|
|
||||||
|
LOGINFO("margin=%ld",prefh_->get_margin());
|
||||||
|
LOGINFO("x0=(%ld,%ld,%ld) lx=(%ld,%ld,%ld)",x0[0],x0[1],x0[2],lx[0],lx[1],lx[2]);
|
||||||
|
|
||||||
if (randc[ilevel] == NULL)
|
if (randc[ilevel] == NULL)
|
||||||
randc[ilevel] = new rng(*randc[ilevel - 1], ran_cube_size_, rngseeds_[ilevel], kavg, ilevel == levelmin_ + 1, x0, lx);
|
randc[ilevel] = new rng(*randc[ilevel - 1], ran_cube_size_, rngseeds_[ilevel], kavg, ilevel == levelmin_ + 1, x0, lx);
|
||||||
|
@ -1812,7 +1826,7 @@ void random_number_generator<rng, T>::store_rnd(int ilevel, rng *prng)
|
||||||
char fname[128];
|
char fname[128];
|
||||||
sprintf(fname, "grafic_wnoise_%04d.bin", ilevel);
|
sprintf(fname, "grafic_wnoise_%04d.bin", ilevel);
|
||||||
|
|
||||||
LOGUSER("Storing white noise field for grafic in file \'%s\'...", fname);
|
LOGINFO("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);
|
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);
|
std::ofstream ofs(fname, std::ios::binary | std::ios::trunc);
|
||||||
|
@ -1884,12 +1898,23 @@ void random_number_generator<rng, T>::store_rnd(int ilevel, rng *prng)
|
||||||
int nx, ny, nz;
|
int nx, ny, nz;
|
||||||
int i0, j0, k0;
|
int i0, j0, k0;
|
||||||
|
|
||||||
nx = 2 * prefh_->size(ilevel, 0);
|
int margin[3];
|
||||||
ny = 2 * prefh_->size(ilevel, 1);
|
if (prefh_->get_margin()>0){
|
||||||
nz = 2 * prefh_->size(ilevel, 2);
|
margin[0] = prefh_->get_margin();
|
||||||
i0 = prefh_->offset_abs(ilevel, 0) - lfac * shift[0] - nx / 4;
|
margin[1] = prefh_->get_margin();
|
||||||
j0 = prefh_->offset_abs(ilevel, 1) - lfac * shift[1] - ny / 4; // was nx/4
|
margin[2] = prefh_->get_margin();
|
||||||
k0 = prefh_->offset_abs(ilevel, 2) - lfac * shift[2] - nz / 4; // was nx/4
|
}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]; //nx / 4;
|
||||||
|
j0 = prefh_->offset_abs(ilevel, 1) - lfac * shift[1] - margin[1]; //ny / 4; // was nx/4
|
||||||
|
k0 = prefh_->offset_abs(ilevel, 2) - lfac * shift[2] - margin[2]; //nz / 4; // was nx/4
|
||||||
|
|
||||||
char fname[128];
|
char fname[128];
|
||||||
sprintf(fname, "wnoise_%04d.bin", ilevel);
|
sprintf(fname, "wnoise_%04d.bin", ilevel);
|
||||||
|
@ -1930,12 +1955,22 @@ void random_number_generator<rng, T>::store_rnd(int ilevel, rng *prng)
|
||||||
}
|
}
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
nx = 2 * prefh_->size(ilevel, 0);
|
int margin[3];
|
||||||
ny = 2 * prefh_->size(ilevel, 1);
|
if (prefh_->get_margin()>0){
|
||||||
nz = 2 * prefh_->size(ilevel, 2);
|
margin[0] = prefh_->get_margin();
|
||||||
i0 = prefh_->offset_abs(ilevel, 0) - lfac * shift[0] - nx / 4;
|
margin[1] = prefh_->get_margin();
|
||||||
j0 = prefh_->offset_abs(ilevel, 1) - lfac * shift[1] - ny / 4; // was nx/4
|
margin[2] = prefh_->get_margin();
|
||||||
k0 = prefh_->offset_abs(ilevel, 2) - lfac * shift[2] - nz / 4; // was nx/4
|
}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];//nx / 4;
|
||||||
|
j0 = prefh_->offset_abs(ilevel, 1) - lfac * shift[1] - margin[1];//ny / 4; // was nx/4
|
||||||
|
k0 = prefh_->offset_abs(ilevel, 2) - lfac * shift[2] - margin[2];//nz / 4; // was nx/4
|
||||||
}
|
}
|
||||||
|
|
||||||
mem_cache_[ilevel - levelmin_] = new std::vector<T>(nx * ny * nz, 0.0);
|
mem_cache_[ilevel - levelmin_] = new std::vector<T>(nx * ny * nz, 0.0);
|
||||||
|
|
321
src/random.hh
321
src/random.hh
|
@ -9,8 +9,8 @@
|
||||||
*/
|
*/
|
||||||
|
|
||||||
//... for testing purposes.............
|
//... for testing purposes.............
|
||||||
//#define DEGRADE_RAND1
|
// #define DEGRADE_RAND1
|
||||||
//#define DEGRADE_RAND2
|
// #define DEGRADE_RAND2
|
||||||
//.....................................
|
//.....................................
|
||||||
|
|
||||||
#ifndef __RANDOM_HH
|
#ifndef __RANDOM_HH
|
||||||
|
@ -31,54 +31,53 @@
|
||||||
#include "mg_operators.hh"
|
#include "mg_operators.hh"
|
||||||
#include "constraints.hh"
|
#include "constraints.hh"
|
||||||
|
|
||||||
|
class RNG_plugin
|
||||||
class RNG_plugin{
|
{
|
||||||
protected:
|
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
|
||||||
public:
|
public:
|
||||||
explicit RNG_plugin( config_file& cf )
|
explicit RNG_plugin(config_file &cf)
|
||||||
: pcf_( &cf )
|
: pcf_(&cf)
|
||||||
{ }
|
{
|
||||||
virtual ~RNG_plugin() { }
|
}
|
||||||
|
virtual ~RNG_plugin() {}
|
||||||
virtual bool is_multiscale() const = 0;
|
virtual bool is_multiscale() const = 0;
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
||||||
struct RNG_plugin_creator
|
struct RNG_plugin_creator
|
||||||
{
|
{
|
||||||
virtual RNG_plugin * create( config_file& cf ) const = 0;
|
virtual RNG_plugin *create(config_file &cf) const = 0;
|
||||||
virtual ~RNG_plugin_creator() { }
|
virtual ~RNG_plugin_creator() {}
|
||||||
};
|
};
|
||||||
|
|
||||||
std::map< std::string, RNG_plugin_creator *>&
|
std::map<std::string, RNG_plugin_creator *> &
|
||||||
get_RNG_plugin_map();
|
get_RNG_plugin_map();
|
||||||
|
|
||||||
void print_RNG_plugins( void );
|
void print_RNG_plugins(void);
|
||||||
|
|
||||||
template< class Derived >
|
template <class Derived>
|
||||||
struct RNG_plugin_creator_concrete : public RNG_plugin_creator
|
struct RNG_plugin_creator_concrete : public RNG_plugin_creator
|
||||||
{
|
{
|
||||||
//! register the plugin by its name
|
//! register the plugin by its name
|
||||||
RNG_plugin_creator_concrete( const std::string& plugin_name )
|
RNG_plugin_creator_concrete(const std::string &plugin_name)
|
||||||
{
|
{
|
||||||
get_RNG_plugin_map()[ plugin_name ] = this;
|
get_RNG_plugin_map()[plugin_name] = this;
|
||||||
}
|
}
|
||||||
|
|
||||||
//! create an instance of the plugin
|
//! create an instance of the plugin
|
||||||
RNG_plugin* create( config_file& cf ) const
|
RNG_plugin *create(config_file &cf) const
|
||||||
{
|
{
|
||||||
return new Derived( cf );
|
return new Derived(cf);
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
|
|
||||||
typedef RNG_plugin RNG_instance;
|
typedef RNG_plugin RNG_instance;
|
||||||
RNG_plugin *select_RNG_plugin( config_file& cf );
|
RNG_plugin *select_RNG_plugin(config_file &cf);
|
||||||
|
|
||||||
|
|
||||||
/*!
|
/*!
|
||||||
* @brief encapsulates all things random number generator related
|
* @brief encapsulates all things random number generator related
|
||||||
*/
|
*/
|
||||||
template< typename T >
|
template <typename T>
|
||||||
class random_numbers
|
class random_numbers
|
||||||
{
|
{
|
||||||
public:
|
public:
|
||||||
|
@ -88,164 +87,159 @@ public:
|
||||||
ncubes_; //!< number of random number cubes to cover the full mesh
|
ncubes_; //!< number of random number cubes to cover the full mesh
|
||||||
long baseseed_; //!< base seed from which cube seeds are computed
|
long baseseed_; //!< base seed from which cube seeds are computed
|
||||||
|
|
||||||
|
|
||||||
protected:
|
protected:
|
||||||
//! vector of 3D meshes (the random number cubes) with random numbers
|
//! vector of 3D meshes (the random number cubes) with random numbers
|
||||||
std::vector< Meshvar<T>* > rnums_;
|
std::vector<Meshvar<T> *> rnums_;
|
||||||
|
|
||||||
//! map of 3D indices to cube index
|
//! map of 3D indices to cube index
|
||||||
std::map<size_t,size_t> cubemap_;
|
std::map<size_t, size_t> cubemap_;
|
||||||
|
|
||||||
typedef std::map<size_t,size_t>::iterator cubemap_iterator;
|
typedef std::map<size_t, size_t>::iterator cubemap_iterator;
|
||||||
|
|
||||||
protected:
|
protected:
|
||||||
|
|
||||||
//! register a cube with the hash map
|
//! register a cube with the hash map
|
||||||
void register_cube( int i, int j, int k);
|
void register_cube(int i, int j, int k);
|
||||||
|
|
||||||
//! fills a subcube with random numbers
|
//! fills a subcube with random numbers
|
||||||
double fill_cube( int i, int j, int k);
|
double fill_cube(int i, int j, int k);
|
||||||
|
|
||||||
//! subtract a constant from an entire cube
|
//! subtract a constant from an entire cube
|
||||||
void subtract_from_cube( int i, int j, int k, double val );
|
void subtract_from_cube(int i, int j, int k, double val);
|
||||||
|
|
||||||
//! copy random numbers from a cube to a full grid array
|
//! copy random numbers from a cube to a full grid array
|
||||||
template< class C >
|
template <class C>
|
||||||
void copy_cube( int i, int j, int k, C& dat )
|
void copy_cube(int i, int j, int k, C &dat)
|
||||||
{
|
{
|
||||||
int offi, offj, offk;
|
int offi, offj, offk;
|
||||||
|
|
||||||
offi = i*cubesize_;
|
offi = i * cubesize_;
|
||||||
offj = j*cubesize_;
|
offj = j * cubesize_;
|
||||||
offk = k*cubesize_;
|
offk = k * cubesize_;
|
||||||
|
|
||||||
i = (i+ncubes_)%ncubes_;
|
i = (i + ncubes_) % ncubes_;
|
||||||
j = (j+ncubes_)%ncubes_;
|
j = (j + ncubes_) % ncubes_;
|
||||||
k = (k+ncubes_)%ncubes_;
|
k = (k + ncubes_) % ncubes_;
|
||||||
|
|
||||||
size_t icube = (i*ncubes_+j)*ncubes_+k;
|
size_t icube = (i * ncubes_ + j) * ncubes_ + k;
|
||||||
cubemap_iterator it = cubemap_.find( icube );
|
cubemap_iterator it = cubemap_.find(icube);
|
||||||
|
|
||||||
if( it == cubemap_.end() )
|
if (it == cubemap_.end())
|
||||||
{
|
{
|
||||||
LOGERR("attempting to copy data from non-existing RND cube %d,%d,%d",i,j,k);
|
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");
|
throw std::runtime_error("attempting to copy data from non-existing RND cube");
|
||||||
}
|
}
|
||||||
|
|
||||||
size_t cubeidx = it->second;
|
size_t cubeidx = it->second;
|
||||||
|
|
||||||
for( int ii=0; ii<(int)cubesize_; ++ii )
|
for (int ii = 0; ii < (int)cubesize_; ++ii)
|
||||||
for( int jj=0; jj<(int)cubesize_; ++jj )
|
for (int jj = 0; jj < (int)cubesize_; ++jj)
|
||||||
for( int kk=0; kk<(int)cubesize_; ++kk )
|
for (int kk = 0; kk < (int)cubesize_; ++kk)
|
||||||
dat(offi+ii,offj+jj,offk+kk) = (*rnums_[cubeidx])(ii,jj,kk);
|
dat(offi + ii, offj + jj, offk + kk) = (*rnums_[cubeidx])(ii, jj, kk);
|
||||||
}
|
}
|
||||||
|
|
||||||
//! free the memory associated with a subcube
|
//! free the memory associated with a subcube
|
||||||
void free_cube( int i, int j, int k );
|
void free_cube(int i, int j, int k);
|
||||||
|
|
||||||
//! initialize member variables and allocate memory
|
//! initialize member variables and allocate memory
|
||||||
void initialize( void );
|
void initialize(void);
|
||||||
|
|
||||||
//! fill a cubic subvolume of the full grid with random numbers
|
//! fill a cubic subvolume of the full grid with random numbers
|
||||||
double fill_subvolume( int *i0, int *n );
|
double fill_subvolume(int *i0, int *n);
|
||||||
|
|
||||||
//! fill an entire grid with random numbers
|
//! fill an entire grid with random numbers
|
||||||
double fill_all( void );
|
double fill_all(void);
|
||||||
|
|
||||||
//! fill an external array instead of the internal field
|
//! fill an external array instead of the internal field
|
||||||
template< class C >
|
template <class C>
|
||||||
double fill_all( C& dat )
|
double fill_all(C &dat)
|
||||||
{
|
{
|
||||||
double sum = 0.0;
|
double sum = 0.0;
|
||||||
|
|
||||||
for( int i=0; i<(int)ncubes_; ++i )
|
for (int i = 0; i < (int)ncubes_; ++i)
|
||||||
for( int j=0; j<(int)ncubes_; ++j )
|
for (int j = 0; j < (int)ncubes_; ++j)
|
||||||
for( int k=0; k<(int)ncubes_; ++k )
|
for (int k = 0; k < (int)ncubes_; ++k)
|
||||||
{
|
{
|
||||||
int ii(i),jj(j),kk(k);
|
int ii(i), jj(j), kk(k);
|
||||||
register_cube(ii,jj,kk);
|
register_cube(ii, jj, kk);
|
||||||
}
|
}
|
||||||
|
|
||||||
#pragma omp parallel for reduction(+:sum)
|
#pragma omp parallel for reduction(+ \
|
||||||
for( int i=0; i<(int)ncubes_; ++i )
|
: sum)
|
||||||
for( int j=0; j<(int)ncubes_; ++j )
|
for (int i = 0; i < (int)ncubes_; ++i)
|
||||||
for( int k=0; k<(int)ncubes_; ++k )
|
for (int j = 0; j < (int)ncubes_; ++j)
|
||||||
|
for (int k = 0; k < (int)ncubes_; ++k)
|
||||||
{
|
{
|
||||||
int ii(i),jj(j),kk(k);
|
int ii(i), jj(j), kk(k);
|
||||||
|
|
||||||
ii = (ii+ncubes_)%ncubes_;
|
ii = (ii + ncubes_) % ncubes_;
|
||||||
jj = (jj+ncubes_)%ncubes_;
|
jj = (jj + ncubes_) % ncubes_;
|
||||||
kk = (kk+ncubes_)%ncubes_;
|
kk = (kk + ncubes_) % ncubes_;
|
||||||
|
|
||||||
sum+=fill_cube(ii, jj, kk);
|
sum += fill_cube(ii, jj, kk);
|
||||||
copy_cube(ii,jj,kk,dat);
|
copy_cube(ii, jj, kk, dat);
|
||||||
free_cube(ii, jj, kk);
|
free_cube(ii, jj, kk);
|
||||||
}
|
}
|
||||||
|
|
||||||
return sum/(ncubes_*ncubes_*ncubes_);
|
return sum / (ncubes_ * ncubes_ * ncubes_);
|
||||||
}
|
}
|
||||||
|
|
||||||
//! write the number of allocated random number cubes to stdout
|
//! write the number of allocated random number cubes to stdout
|
||||||
void print_allocated( void );
|
void print_allocated(void);
|
||||||
|
|
||||||
public:
|
public:
|
||||||
|
|
||||||
//! constructor
|
//! constructor
|
||||||
random_numbers( unsigned res, unsigned cubesize, long baseseed, int *x0, int *lx );
|
random_numbers(unsigned res, unsigned cubesize, long baseseed, int *x0, int *lx);
|
||||||
|
|
||||||
//! constructor for constrained fine field
|
//! constructor for constrained fine field
|
||||||
random_numbers( random_numbers<T>& rc, unsigned cubesize, long baseseed,
|
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 );
|
bool kspace = false, bool isolated = false, int *x0_ = NULL, int *lx_ = NULL, bool zeromean = true);
|
||||||
|
|
||||||
//! constructor
|
//! constructor
|
||||||
random_numbers( unsigned res, unsigned cubesize, long baseseed, bool zeromean=true );
|
random_numbers(unsigned res, unsigned cubesize, long baseseed, bool zeromean = true);
|
||||||
|
|
||||||
|
|
||||||
//! constructor to read white noise from file
|
//! constructor to read white noise from file
|
||||||
random_numbers( unsigned res, std::string randfname, bool rndsign );
|
random_numbers(unsigned res, std::string randfname, bool rndsign);
|
||||||
|
|
||||||
|
|
||||||
//! copy constructor for averaged field (not copying) hence explicit!
|
//! copy constructor for averaged field (not copying) hence explicit!
|
||||||
explicit random_numbers( /*const*/ random_numbers <T>& rc, bool kdegrade = true );
|
explicit random_numbers(/*const*/ random_numbers<T> &rc, bool kdegrade = true);
|
||||||
|
|
||||||
//! destructor
|
//! destructor
|
||||||
~random_numbers()
|
~random_numbers()
|
||||||
{
|
{
|
||||||
for( unsigned i=0; i<rnums_.size(); ++i )
|
for (unsigned i = 0; i < rnums_.size(); ++i)
|
||||||
if( rnums_[i] != NULL )
|
if (rnums_[i] != NULL)
|
||||||
delete rnums_[i];
|
delete rnums_[i];
|
||||||
rnums_.clear();
|
rnums_.clear();
|
||||||
}
|
}
|
||||||
|
|
||||||
//! access a random number, this allocates a cube and fills it with consistent random numbers
|
//! 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 )
|
inline T &operator()(int i, int j, int k, bool fillrand = true)
|
||||||
{
|
{
|
||||||
int ic, jc, kc, is, js, ks;
|
int ic, jc, kc, is, js, ks;
|
||||||
|
|
||||||
if( ncubes_ == 0 )
|
if (ncubes_ == 0)
|
||||||
throw std::runtime_error("random_numbers: internal error, not properly initialized");
|
throw std::runtime_error("random_numbers: internal error, not properly initialized");
|
||||||
|
|
||||||
//... determine cube
|
//... determine cube
|
||||||
ic = (int)((double)i/cubesize_ + ncubes_) % ncubes_;
|
ic = (int)((double)i / cubesize_ + ncubes_) % ncubes_;
|
||||||
jc = (int)((double)j/cubesize_ + ncubes_) % ncubes_;
|
jc = (int)((double)j / cubesize_ + ncubes_) % ncubes_;
|
||||||
kc = (int)((double)k/cubesize_ + ncubes_) % ncubes_;
|
kc = (int)((double)k / cubesize_ + ncubes_) % ncubes_;
|
||||||
|
|
||||||
size_t icube = ((size_t)ic*ncubes_+(size_t)jc)*ncubes_+(size_t)kc;
|
size_t icube = ((size_t)ic * ncubes_ + (size_t)jc) * ncubes_ + (size_t)kc;
|
||||||
|
|
||||||
cubemap_iterator it = cubemap_.find( icube );
|
cubemap_iterator it = cubemap_.find(icube);
|
||||||
|
|
||||||
if( it == cubemap_.end() )
|
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);
|
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");
|
throw std::runtime_error("attempting to copy data from non-existing RND cube");
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
size_t cubeidx = it->second;
|
size_t cubeidx = it->second;
|
||||||
|
|
||||||
if( rnums_[ cubeidx ] == NULL )
|
if (rnums_[cubeidx] == NULL)
|
||||||
{
|
{
|
||||||
LOGERR("Attempting to access data from non-allocated RND cube %d,%d,%d",ic,jc,kc);
|
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");
|
throw std::runtime_error("attempting to access data from non-allocated RND cube");
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -254,33 +248,30 @@ public:
|
||||||
js = (j - jc * cubesize_ + cubesize_) % cubesize_;
|
js = (j - jc * cubesize_ + cubesize_) % cubesize_;
|
||||||
ks = (k - kc * cubesize_ + cubesize_) % cubesize_;
|
ks = (k - kc * cubesize_ + cubesize_) % cubesize_;
|
||||||
|
|
||||||
return (*rnums_[ cubeidx ])(is,js,ks);
|
return (*rnums_[cubeidx])(is, js, ks);
|
||||||
}
|
}
|
||||||
|
|
||||||
//! free all cubes
|
//! free all cubes
|
||||||
void free_all_mem( void )
|
void free_all_mem(void)
|
||||||
{
|
{
|
||||||
for( unsigned i=0; i<rnums_.size(); ++i )
|
for (unsigned i = 0; i < rnums_.size(); ++i)
|
||||||
if( rnums_[i] != NULL )
|
if (rnums_[i] != NULL)
|
||||||
{
|
{
|
||||||
delete rnums_[i];
|
delete rnums_[i];
|
||||||
rnums_[i] = NULL;
|
rnums_[i] = NULL;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
||||||
/*!
|
/*!
|
||||||
* @brief encapsulates all things for multi-scale white noise generation
|
* @brief encapsulates all things for multi-scale white noise generation
|
||||||
*/
|
*/
|
||||||
template< typename rng, typename T >
|
template <typename rng, typename T>
|
||||||
class random_number_generator
|
class random_number_generator
|
||||||
{
|
{
|
||||||
protected:
|
protected:
|
||||||
config_file * pcf_;
|
config_file *pcf_;
|
||||||
refinement_hierarchy * prefh_;
|
refinement_hierarchy *prefh_;
|
||||||
constraint_set constraints;
|
constraint_set constraints;
|
||||||
|
|
||||||
int levelmin_,
|
int levelmin_,
|
||||||
|
@ -291,86 +282,82 @@ protected:
|
||||||
|
|
||||||
bool disk_cached_;
|
bool disk_cached_;
|
||||||
bool restart_;
|
bool restart_;
|
||||||
std::vector< std::vector<T>* > mem_cache_;
|
std::vector<std::vector<T> *> mem_cache_;
|
||||||
|
|
||||||
unsigned ran_cube_size_;
|
unsigned ran_cube_size_;
|
||||||
|
|
||||||
|
|
||||||
protected:
|
protected:
|
||||||
|
|
||||||
//! checks if the specified string is numeric
|
//! checks if the specified string is numeric
|
||||||
bool is_number(const std::string& s);
|
bool is_number(const std::string &s);
|
||||||
|
|
||||||
//! parses the random number parameters in the conf file
|
//! parses the random number parameters in the conf file
|
||||||
void parse_rand_parameters( void );
|
void parse_rand_parameters(void);
|
||||||
|
|
||||||
//! correct coarse grid averages for the change in small scale when using Fourier interpolation
|
//! correct coarse grid averages for the change in small scale when using Fourier interpolation
|
||||||
void correct_avg( int icoarse, int ifine );
|
void correct_avg(int icoarse, int ifine);
|
||||||
|
|
||||||
//! the main driver routine for multi-scale white noise generation
|
//! the main driver routine for multi-scale white noise generation
|
||||||
void compute_random_numbers( void );
|
void compute_random_numbers(void);
|
||||||
|
|
||||||
//! store the white noise fields in memory or on disk
|
//! store the white noise fields in memory or on disk
|
||||||
void store_rnd( int ilevel, rng* prng );
|
void store_rnd(int ilevel, rng *prng);
|
||||||
|
|
||||||
|
|
||||||
public:
|
public:
|
||||||
|
|
||||||
//! constructor
|
//! constructor
|
||||||
random_number_generator( config_file& cf, refinement_hierarchy& refh, transfer_function *ptf = NULL );
|
random_number_generator(config_file &cf, refinement_hierarchy &refh, transfer_function *ptf = NULL);
|
||||||
|
|
||||||
//! destructor
|
//! destructor
|
||||||
~random_number_generator();
|
~random_number_generator();
|
||||||
|
|
||||||
//! load random numbers to a new array
|
//! load random numbers to a new array
|
||||||
template< typename array >
|
template <typename array>
|
||||||
void load( array& A, int ilevel )
|
void load(array &A, int ilevel)
|
||||||
{
|
{
|
||||||
if( restart_ )
|
if (restart_)
|
||||||
LOGINFO("Attempting to restart using random numbers for level %d\n from file \'wnoise_%04d.bin\'.",ilevel,ilevel);
|
LOGINFO("Attempting to restart using random numbers for level %d\n from file \'wnoise_%04d.bin\'.", ilevel, ilevel);
|
||||||
|
|
||||||
if( disk_cached_ )
|
if (disk_cached_)
|
||||||
{
|
{
|
||||||
char fname[128];
|
char fname[128];
|
||||||
sprintf(fname,"wnoise_%04d.bin",ilevel);
|
sprintf(fname, "wnoise_%04d.bin", ilevel);
|
||||||
|
|
||||||
LOGUSER("Loading white noise from file \'%s\'...",fname);
|
LOGUSER("Loading white noise from file \'%s\'...", fname);
|
||||||
|
|
||||||
std::ifstream ifs( fname, std::ios::binary );
|
std::ifstream ifs(fname, std::ios::binary);
|
||||||
if( !ifs.good() )
|
if (!ifs.good())
|
||||||
{
|
{
|
||||||
LOGERR("White noise file \'%s\'was not found.",fname);
|
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.");
|
throw std::runtime_error("A white noise file was not found. This is an internal inconsistency and bad.");
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
int nx,ny,nz;
|
int nx, ny, nz;
|
||||||
ifs.read( reinterpret_cast<char*> (&nx), sizeof(int) );
|
ifs.read(reinterpret_cast<char *>(&nx), sizeof(int));
|
||||||
ifs.read( reinterpret_cast<char*> (&ny), sizeof(int) );
|
ifs.read(reinterpret_cast<char *>(&ny), sizeof(int));
|
||||||
ifs.read( reinterpret_cast<char*> (&nz), 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) || 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 )
|
if (nx == (int)A.size(0) * 2 && ny == (int)A.size(1) * 2 && nz == (int)A.size(2) * 2)
|
||||||
{
|
{
|
||||||
std::cerr << "CHECKPOINT" << std::endl;
|
std::cerr << "CHECKPOINT" << std::endl;
|
||||||
|
|
||||||
|
int ox = nx / 4, oy = ny / 4, oz = nz / 4;
|
||||||
|
std::vector<T> slice(ny * nz, 0.0);
|
||||||
|
|
||||||
int ox = nx/4, oy = ny/4, oz = nz/4;
|
for (int i = 0; i < nx; ++i)
|
||||||
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) );
|
ifs.read(reinterpret_cast<char *>(&slice[0]), ny * nz * sizeof(T));
|
||||||
|
|
||||||
if( i<ox ) continue;
|
if (i < ox)
|
||||||
if( i>=3*ox ) break;
|
continue;
|
||||||
|
if (i >= 3 * ox)
|
||||||
|
break;
|
||||||
|
|
||||||
#pragma omp parallel for
|
#pragma omp parallel for
|
||||||
for( int j=oy; j<3*oy; ++j )
|
for (int j = oy; j < 3 * oy; ++j)
|
||||||
for( int k=oz; k<3*oz; ++k )
|
for (int k = oz; k < 3 * oz; ++k)
|
||||||
A(i-ox,j-oy,k-oz) = slice[j*nz+k];
|
A(i - ox, j - oy, k - oz) = slice[j * nz + k];
|
||||||
}
|
}
|
||||||
|
|
||||||
ifs.close();
|
ifs.close();
|
||||||
|
@ -378,21 +365,22 @@ public:
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
LOGERR("White noise file is not aligned with array. File: [%d,%d,%d]. Mem: [%d,%d,%d].",
|
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));
|
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.");
|
throw std::runtime_error("White noise file is not aligned with array. This is an internal inconsistency and bad.");
|
||||||
}
|
}
|
||||||
}else{
|
}
|
||||||
|
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 i = 0; i < nx; ++i)
|
||||||
for( int j=0; j<ny; ++j )
|
{
|
||||||
for( int k=0; k<nz; ++k )
|
std::vector<T> slice(ny * nz, 0.0);
|
||||||
A(i,j,k) = slice[j*nz+k];
|
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();
|
ifs.close();
|
||||||
|
@ -402,36 +390,31 @@ public:
|
||||||
{
|
{
|
||||||
LOGUSER("Copying white noise from memory cache...");
|
LOGUSER("Copying white noise from memory cache...");
|
||||||
|
|
||||||
if( mem_cache_[ilevel-levelmin_] == NULL )
|
if (mem_cache_[ilevel - levelmin_] == NULL)
|
||||||
LOGERR("Tried to access mem-cached random numbers for level %d. But these are not available!\n",ilevel);
|
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) );
|
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() )
|
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));
|
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");
|
throw std::runtime_error("White noise file is not aligned with array. This is an internal inconsistency and bad");
|
||||||
}
|
}
|
||||||
|
|
||||||
#pragma omp parallel for
|
#pragma omp parallel for
|
||||||
for( int i=0; i<nx; ++i )
|
for (int i = 0; i < nx; ++i)
|
||||||
for( int j=0; j<ny; ++j )
|
for (int j = 0; j < ny; ++j)
|
||||||
for( int k=0; k<nz; ++k )
|
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];
|
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;
|
|
||||||
|
|
||||||
|
std::vector<T>().swap(*mem_cache_[ilevel - levelmin_]);
|
||||||
|
delete mem_cache_[ilevel - levelmin_];
|
||||||
|
mem_cache_[ilevel - levelmin_] = NULL;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
|
|
||||||
typedef random_numbers<real_t> rand_nums;
|
typedef random_numbers<real_t> rand_nums;
|
||||||
typedef random_number_generator< rand_nums,real_t> rand_gen;
|
typedef random_number_generator<rand_nums, real_t> rand_gen;
|
||||||
|
|
||||||
|
|
||||||
#endif //__RANDOM_HH
|
#endif //__RANDOM_HH
|
||||||
|
|
||||||
|
|
Loading…
Reference in a new issue