mirror of
https://github.com/cosmo-sims/MUSIC.git
synced 2024-09-19 17:03:46 +02:00
Re-enabled explicit averaging when using hybrid+deconvolution
This commit is contained in:
parent
b730cd5b5f
commit
39d679685c
3 changed files with 175 additions and 13 deletions
|
@ -403,7 +403,13 @@ namespace convolution{
|
||||||
dplus = pcf_->getValue<double>("cosmology","dplus");
|
dplus = pcf_->getValue<double>("cosmology","dplus");
|
||||||
|
|
||||||
bool
|
bool
|
||||||
bperiodic = pcf_->getValueSafe<bool>("setup","periodic_TF",true);
|
bperiodic = pcf_->getValueSafe<bool>("setup","periodic_TF",true),
|
||||||
|
deconv = pcf_->getValueSafe<bool>("setup","deconvolve",true);
|
||||||
|
// bool deconv_baryons = true;//pcf_->getValueSafe<bool>("setup","deconvolve_baryons",false) || do_SPH;
|
||||||
|
bool bsmooth_baryons = false;//type==baryon && !deconv_baryons;
|
||||||
|
//bool bbaryons = type==baryon | type==vbaryon;
|
||||||
|
bool kspacepoisson = ((pcf_->getValueSafe<bool>("poisson","fft_fine",true)|
|
||||||
|
pcf_->getValueSafe<bool>("poisson","kspace",false)));// & !(type==baryon&!do_SPH);//&!baryons ;
|
||||||
|
|
||||||
std::cout << " - Computing transfer function kernel...\n";
|
std::cout << " - Computing transfer function kernel...\n";
|
||||||
|
|
||||||
|
@ -425,6 +431,11 @@ namespace convolution{
|
||||||
|
|
||||||
LOGUSER("Computing fine kernel (level %d)...", levelmax);
|
LOGUSER("Computing fine kernel (level %d)...", levelmax);
|
||||||
|
|
||||||
|
int ref_fac = (deconv&&kspacepoisson)? 2 : 0;
|
||||||
|
const int ql = -ref_fac/2+1, qr = ql+ref_fac;
|
||||||
|
const double rf8 = pow(ref_fac,3);
|
||||||
|
const double dx05 = 0.5*dx, dx025 = 0.25*dx;
|
||||||
|
|
||||||
if( bperiodic )
|
if( bperiodic )
|
||||||
{
|
{
|
||||||
|
|
||||||
|
@ -472,8 +483,35 @@ namespace convolution{
|
||||||
&& rr[1] > -boxlength && rr[1] <= boxlength
|
&& rr[1] > -boxlength && rr[1] <= boxlength
|
||||||
&& rr[2] > -boxlength && rr[2] <= boxlength )
|
&& rr[2] > -boxlength && rr[2] <= boxlength )
|
||||||
{
|
{
|
||||||
rr2 = rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2];
|
if( ref_fac > 0 )
|
||||||
val += tfr->compute_real(rr2);
|
{
|
||||||
|
double rrr[3];
|
||||||
|
register double rrr2[3];
|
||||||
|
for( int iii=ql; iii<qr; ++iii )
|
||||||
|
{
|
||||||
|
rrr[0] = rr[0]+(double)iii*dx05 - dx025;
|
||||||
|
rrr2[0]= rrr[0]*rrr[0];
|
||||||
|
for( int jjj=ql; jjj<qr; ++jjj )
|
||||||
|
{
|
||||||
|
rrr[1] = rr[1]+(double)jjj*dx05 - dx025;
|
||||||
|
rrr2[1]= rrr[1]*rrr[1];
|
||||||
|
for( int kkk=ql; kkk<qr; ++kkk )
|
||||||
|
{
|
||||||
|
rrr[2] = rr[2]+(double)kkk*dx05 - dx025;
|
||||||
|
rrr2[2]= rrr[2]*rrr[2];
|
||||||
|
rr2 = rrr2[0]+rrr2[1]+rrr2[2];
|
||||||
|
val += tfr->compute_real(rr2)/rf8;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
rr2 = rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2];
|
||||||
|
val += tfr->compute_real(rr2);
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -505,7 +543,7 @@ namespace convolution{
|
||||||
|
|
||||||
//rkernel[idx] = 0.0;
|
//rkernel[idx] = 0.0;
|
||||||
|
|
||||||
rr2 = rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2];
|
//rr2 = rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2];
|
||||||
|
|
||||||
//... speed up 8x by copying data to other octants
|
//... speed up 8x by copying data to other octants
|
||||||
size_t idx[8];
|
size_t idx[8];
|
||||||
|
@ -523,8 +561,37 @@ namespace convolution{
|
||||||
if(j==0||j==ny/2){ idx[2]=idx[3]=idx[6]=idx[7]=(size_t)-1;}
|
if(j==0||j==ny/2){ idx[2]=idx[3]=idx[6]=idx[7]=(size_t)-1;}
|
||||||
if(k==0||k==nz/2){ idx[4]=idx[5]=idx[6]=idx[7]=(size_t)-1;}
|
if(k==0||k==nz/2){ idx[4]=idx[5]=idx[6]=idx[7]=(size_t)-1;}
|
||||||
|
|
||||||
|
double val = 0.0;//(fftw_real)tfr->compute_real(rr2)*fac;
|
||||||
|
|
||||||
|
if( ref_fac > 0 )
|
||||||
|
{
|
||||||
|
double rrr[3];
|
||||||
|
register double rrr2[3];
|
||||||
|
for( int iii=ql; iii<qr; ++iii )
|
||||||
|
{
|
||||||
|
rrr[0] = rr[0]+(double)iii*dx05 - dx025;
|
||||||
|
rrr2[0]= rrr[0]*rrr[0];
|
||||||
|
for( int jjj=ql; jjj<qr; ++jjj )
|
||||||
|
{
|
||||||
|
rrr[1] = rr[1]+(double)jjj*dx05 - dx025;
|
||||||
|
rrr2[1]= rrr[1]*rrr[1];
|
||||||
|
for( int kkk=ql; kkk<qr; ++kkk )
|
||||||
|
{
|
||||||
|
rrr[2] = rr[2]+(double)kkk*dx05 - dx025;
|
||||||
|
rrr2[2]= rrr[2]*rrr[2];
|
||||||
|
rr2 = rrr2[0]+rrr2[1]+rrr2[2];
|
||||||
|
val += tfr->compute_real(rr2)/rf8;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
rr2 = rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2];
|
||||||
|
val = tfr->compute_real(rr2);
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
double val = (fftw_real)tfr->compute_real(rr2)*fac;
|
|
||||||
//if( rr2 <= boxlength2*boxlength2 )
|
//if( rr2 <= boxlength2*boxlength2 )
|
||||||
//rkernel[idx] += (fftw_real)tfr->compute_real(rr2)*fac;
|
//rkernel[idx] += (fftw_real)tfr->compute_real(rr2)*fac;
|
||||||
|
|
||||||
|
@ -541,10 +608,7 @@ namespace convolution{
|
||||||
/*************************************************************************************/
|
/*************************************************************************************/
|
||||||
/*************************************************************************************/
|
/*************************************************************************************/
|
||||||
/******* perform deconvolution *******************************************************/
|
/******* perform deconvolution *******************************************************/
|
||||||
bool deconv = pcf_->getValueSafe<bool>("setup","deconvolve",true);
|
|
||||||
// bool deconv_baryons = true;//pcf_->getValueSafe<bool>("setup","deconvolve_baryons",false) || do_SPH;
|
|
||||||
bool bsmooth_baryons = false;//type==baryon && !deconv_baryons;
|
|
||||||
//bool bbaryons = type==baryon | type==vbaryon;
|
|
||||||
|
|
||||||
|
|
||||||
//bool baryons = type==baryon||type==vbaryon;
|
//bool baryons = type==baryon||type==vbaryon;
|
||||||
|
@ -554,8 +618,7 @@ namespace convolution{
|
||||||
LOGUSER("Deconvolving fine kernel...");
|
LOGUSER("Deconvolving fine kernel...");
|
||||||
std::cout << " - Deconvoling density kernel...\n";
|
std::cout << " - Deconvoling density kernel...\n";
|
||||||
|
|
||||||
bool kspacepoisson = ((pcf_->getValueSafe<bool>("poisson","fft_fine",true)|
|
|
||||||
pcf_->getValueSafe<bool>("poisson","kspace",false)));// & !(type==baryon&!do_SPH);//&!baryons ;
|
|
||||||
|
|
||||||
double fftnorm = 1.0/((size_t)nx*(size_t)ny*(size_t)nz);
|
double fftnorm = 1.0/((size_t)nx*(size_t)ny*(size_t)nz);
|
||||||
double k0 = rkernel[0];
|
double k0 = rkernel[0];
|
||||||
|
|
102
random.cc
102
random.cc
|
@ -11,6 +11,8 @@
|
||||||
#include "random.hh"
|
#include "random.hh"
|
||||||
|
|
||||||
|
|
||||||
|
//#define DEGRADE_RAND1
|
||||||
|
|
||||||
template< typename T >
|
template< typename T >
|
||||||
random_numbers<T>::random_numbers( unsigned res, unsigned cubesize, long baseseed, int *x0, int *lx )
|
random_numbers<T>::random_numbers( unsigned res, unsigned cubesize, long baseseed, int *x0, int *lx )
|
||||||
: res_( res ), cubesize_( cubesize ), ncubes_( 1 ), baseseed_( baseseed )
|
: res_( res ), cubesize_( cubesize ), ncubes_( 1 ), baseseed_( baseseed )
|
||||||
|
@ -1288,6 +1290,104 @@ void random_number_generator<rng,T>:: store_rnd( int ilevel, rng* prng )
|
||||||
|
|
||||||
int lfac = 1<<(ilevel-levelmin_poisson);
|
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 );
|
||||||
|
|
||||||
|
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_ )
|
if( disk_cached_ )
|
||||||
{
|
{
|
||||||
std::vector<T> data;
|
std::vector<T> data;
|
||||||
|
@ -1314,7 +1414,7 @@ void random_number_generator<rng,T>:: store_rnd( int ilevel, rng* prng )
|
||||||
data.assign( N*N, 0.0 );
|
data.assign( N*N, 0.0 );
|
||||||
for( int i=0; i<N; ++i )
|
for( int i=0; i<N; ++i )
|
||||||
{
|
{
|
||||||
#pragma omp parallel for
|
#pragma omp parallel for
|
||||||
for( int j=0; j<N; ++j )
|
for( int j=0; j<N; ++j )
|
||||||
for( int k=0; k<N; ++k )
|
for( int k=0; k<N; ++k )
|
||||||
data[j*N+k] = (*prng)(i+i0,j+j0,k+k0);
|
data[j*N+k] = (*prng)(i+i0,j+j0,k+k0);
|
||||||
|
|
|
@ -146,7 +146,6 @@ public:
|
||||||
if( type == total )
|
if( type == total )
|
||||||
fname = "input_powerspec_total.txt";
|
fname = "input_powerspec_total.txt";
|
||||||
|
|
||||||
|
|
||||||
if( type == cdm || type == baryon || type == total )
|
if( type == cdm || type == baryon || type == total )
|
||||||
{
|
{
|
||||||
std::ofstream ofs(fname.c_str());
|
std::ofstream ofs(fname.c_str());
|
||||||
|
|
Loading…
Reference in a new issue