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

added missing random number plugin

Browse source
This commit is contained in:
Oliver Hahn 2013-10-24 11:27:09 +02:00
parent ccaf6762bd
commit c74e6dbd55
2 changed files with 19 additions and 196 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

19
plugins/random_music.cc Normal file
View file

@ -0,0 +1,19 @@
#include "random.hh"
class RNG_music : public RNG_plugin{
public:
explicit RNG_music( config_file& cf )
: RNG_plugin( cf )
{ }
~RNG_music() { }
bool is_multiscale() const
{
}
};
namespace{
RNG_plugin_creator_concrete< RNG_music > creator("MUSIC");
}

196
random.cc
View file

@ -1375,193 +1375,6 @@ void random_number_generator<rng,T>::parse_rand_parameters( void )
}
template< typename rng, typename T >
void random_number_generator<rng,T>::correct_avg( int icoarse, int ifine )
{
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 lfacc = 1<<(icoarse-levelmin_poisson);
int nc[3], i0c[3], nf[3], i0f[3];
if( icoarse != levelmin_ )
{
nc[0] = 2*prefh_->size(icoarse, 0);
nc[1] = 2*prefh_->size(icoarse, 1);
nc[2] = 2*prefh_->size(icoarse, 2);
i0c[0] = prefh_->offset_abs(icoarse, 0) - lfacc*shift[0] - nc[0]/4;
i0c[1] = prefh_->offset_abs(icoarse, 1) - lfacc*shift[1] - nc[1]/4;
i0c[2] = prefh_->offset_abs(icoarse, 2) - lfacc*shift[2] - nc[2]/4;
}
else
{
nc[0] = prefh_->size(icoarse, 0);
nc[1] = prefh_->size(icoarse, 1);
nc[2] = prefh_->size(icoarse, 2);
i0c[0] = - lfacc*shift[0];
i0c[1] = - lfacc*shift[1];
i0c[2] = - lfacc*shift[2];
}
nf[0] = 2*prefh_->size(ifine, 0);
nf[1] = 2*prefh_->size(ifine, 1);
nf[2] = 2*prefh_->size(ifine, 2);
i0f[0] = prefh_->offset_abs(ifine, 0) - 2*lfacc*shift[0] - nf[0]/4;
i0f[1] = prefh_->offset_abs(ifine, 1) - 2*lfacc*shift[1] - nf[1]/4;
i0f[2] = prefh_->offset_abs(ifine, 2) - 2*lfacc*shift[2] - nf[2]/4;
//.................................
if( disk_cached_ )
{
char fncoarse[128], fnfine[128];
sprintf(fncoarse,"wnoise_%04d.bin",icoarse);
sprintf(fnfine,"wnoise_%04d.bin",ifine);
std::ifstream
iffine( fnfine, std::ios::binary ),
ifcoarse( fncoarse, std::ios::binary );
int nxc,nyc,nzc,nxf,nyf,nzf;
iffine.read( reinterpret_cast<char*> (&nxf), sizeof(unsigned) );
iffine.read( reinterpret_cast<char*> (&nyf), sizeof(unsigned) );
iffine.read( reinterpret_cast<char*> (&nzf), sizeof(unsigned) );
ifcoarse.read( reinterpret_cast<char*> (&nxc), sizeof(unsigned) );
ifcoarse.read( reinterpret_cast<char*> (&nyc), sizeof(unsigned) );
ifcoarse.read( reinterpret_cast<char*> (&nzc), sizeof(unsigned) );
if( nxf!=nf[0] || nyf!=nf[1] || nzf!=nf[2] || nxc!=nc[0] || nyc!=nc[1] || nzc!=nc[2] )
{
LOGERR("White noise file mismatch. This should not happen. Notify a developer!");
throw std::runtime_error("White noise file mismatch. This should not happen. Notify a developer!");
}
int nxd(nxf/2),nyd(nyf/2),nzd(nzf/2);
std::vector<T> deg_rand( (size_t)nxd*(size_t)nyd*(size_t)nzd, 0.0 );
double fac = 1.0/sqrt(8.0);
for( int i=0, ic=0; i<nxf; i+=2, ic++ )
{
std::vector<T> fine_rand( 2*nyf*nzf, 0.0 );
iffine.read( reinterpret_cast<char*> (&fine_rand[0]), 2*nyf*nzf*sizeof(T) );
#pragma omp parallel for
for( int j=0; j<nyf; j+=2 )
for( int k=0; k<nzf; k+=2 )
{
int jc = j/2, kc = k/2;
//size_t qc = (((size_t)i/2)*(size_t)nyd+((size_t)j/2))*(size_t)nzd+((size_t)k/2);
size_t qc = ((size_t)(ic*nyd+jc))*(size_t)nzd+(size_t)kc;
size_t qf[8];
qf[0] = (0*(size_t)nyf+(size_t)j+0)*(size_t)nzf+(size_t)k+0;
qf[1] = (0*(size_t)nyf+(size_t)j+0)*(size_t)nzf+(size_t)k+1;
qf[2] = (0*(size_t)nyf+(size_t)j+1)*(size_t)nzf+(size_t)k+0;
qf[3] = (0*(size_t)nyf+(size_t)j+1)*(size_t)nzf+(size_t)k+1;
qf[4] = (1*(size_t)nyf+(size_t)j+0)*(size_t)nzf+(size_t)k+0;
qf[5] = (1*(size_t)nyf+(size_t)j+0)*(size_t)nzf+(size_t)k+1;
qf[6] = (1*(size_t)nyf+(size_t)j+1)*(size_t)nzf+(size_t)k+0;
qf[7] = (1*(size_t)nyf+(size_t)j+1)*(size_t)nzf+(size_t)k+1;
double d = 0.0;
for( int q=0; q<8; ++q )
d += fac*fine_rand[qf[q]];
//deg_rand[qc] += d;
deg_rand[qc] = d;
}
}
//... now deg_rand holds the oct-averaged fine field, store this in the coarse field
std::vector<T> coarse_rand(nxc*nyc*nzc,0.0);
ifcoarse.read( reinterpret_cast<char*> (&coarse_rand[0]), nxc*nyc*nzc*sizeof(T) );
int di,dj,dk;
di = i0f[0]/2-i0c[0];
dj = i0f[1]/2-i0c[1];
dk = i0f[2]/2-i0c[2];
#pragma omp parallel for
for( int i=0; i<nxd; i++ )
for( int j=0; j<nyd; j++ )
for( int k=0; k<nzd; k++ )
{
//unsigned qc = (((i+di+nxc)%nxc)*nyc+(((j+dj+nyc)%nyc)))*nzc+((k+dk+nzc)%nzc);
if( i+di < 0 || i+di >= nxc || j+dj < 0 || j+dj >= nyc || k+dk < 0 || k+dk >= nzc )
continue;
size_t qc = (((size_t)i+(size_t)di)*(size_t)nyc+((size_t)j+(size_t)dj))*(size_t)nzc+(size_t)(k+dk);
size_t qcd = (size_t)(i*nyd+j)*(size_t)nzd+(size_t)k;
coarse_rand[qc] = deg_rand[qcd];
}
deg_rand.clear();
ifcoarse.close();
std::ofstream ofcoarse( fncoarse, std::ios::binary|std::ios::trunc );
ofcoarse.write( reinterpret_cast<char*> (&nxc), sizeof(unsigned) );
ofcoarse.write( reinterpret_cast<char*> (&nyc), sizeof(unsigned) );
ofcoarse.write( reinterpret_cast<char*> (&nzc), sizeof(unsigned) );
ofcoarse.write( reinterpret_cast<char*> (&coarse_rand[0]), nxc*nyc*nzc*sizeof(T) );
ofcoarse.close();
}
else
{
int nxc,nyc,nzc,nxf,nyf,nzf;
nxc = nc[0]; nyc = nc[1]; nzc = nc[2];
nxf = nf[0]; nyf = nf[1]; nzf = nf[2];
int nxd(nxf/2),nyd(nyf/2),nzd(nzf/2);
int di,dj,dk;
di = i0f[0]/2-i0c[0];
dj = i0f[1]/2-i0c[1];
dk = i0f[2]/2-i0c[2];
double fac = 1.0/sqrt(8.0);
#pragma omp parallel for
for( int i=0; i<nxd; i++ )
for( int j=0; j<nyd; j++ )
for( int k=0; k<nzd; k++ )
{
if( i+di < 0 || i+di >= nxc || j+dj < 0 || j+dj >= nyc || k+dk < 0 || k+dk >= nzc )
continue;
size_t qf[8];
qf[0] = (size_t)((2*i+0)*nyf+2*j+0)*(size_t)nzf+(size_t)(2*k+0);
qf[1] = (size_t)((2*i+0)*nyf+2*j+0)*(size_t)nzf+(size_t)(2*k+1);
qf[2] = (size_t)((2*i+0)*nyf+2*j+1)*(size_t)nzf+(size_t)(2*k+0);
qf[3] = (size_t)((2*i+0)*nyf+2*j+1)*(size_t)nzf+(size_t)(2*k+1);
qf[4] = (size_t)((2*i+1)*nyf+2*j+0)*(size_t)nzf+(size_t)(2*k+0);
qf[5] = (size_t)((2*i+1)*nyf+2*j+0)*(size_t)nzf+(size_t)(2*k+1);
qf[6] = (size_t)((2*i+1)*nyf+2*j+1)*(size_t)nzf+(size_t)(2*k+0);
qf[7] = (size_t)((2*i+1)*nyf+2*j+1)*(size_t)nzf+(size_t)(2*k+1);
double finesum = 0.0;
for( int q=0; q<8; ++q )
finesum += fac*(*mem_cache_[ifine-levelmin_])[qf[q]];
size_t qc = ((size_t)(i+di)*nyc+(size_t)(j+dj))*(size_t)nzc+(size_t)(k+dk);
(*mem_cache_[icoarse-levelmin_])[qc] = finesum;
}
}
}
template< typename rng, typename T >
void random_number_generator<rng,T>::compute_random_numbers( void )
{
@ -1685,15 +1498,6 @@ void random_number_generator<rng,T>::compute_random_numbers( void )
delete randc[levelmax_];
randc[levelmax_] = NULL;
//... make sure that the coarse grid contains oct averages where it overlaps with a fine grid
//... this also ensures that constraints enforced on fine grids are carried to the coarser grids
for( int ilevel=levelmax_; ilevel>levelmin_; --ilevel )
correct_avg( ilevel-1, ilevel );
//.. we do not have random numbers for a coarse level, generate them
/*if( levelmax_rand_ >= (int)levelmin_ )
{