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

working backup commit to repository. maybe this one should not be used.

Browse source
This commit is contained in:
Oliver Hahn 2013-10-29 09:45:47 +01:00
parent 5f467e81d5
commit 3129083718
4 changed files with 324 additions and 22 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

127
convolution_kernel.cc
View file

@ -530,6 +530,74 @@ namespace convolution{
return this;
}
template< typename real_t >
inline real_t sqr( real_t x )
{ return x*x; }
template< typename real_t >
inline real_t eval_split_recurse( const TransferFunction_real* tfr, real_t *xmid, real_t dx, real_t prevval, int nsplit )
{
const real_t abs_err = 1e-10, rel_err = 1e-6;
const int nmaxsplits = 10;
real_t dxnew = dx/2, dxnew2 = dx/4;
real_t dV = dxnew*dxnew*dxnew;
real_t xl[3] = {xmid[0]-dxnew2,xmid[1]-dxnew2,xmid[2]-dxnew2};
real_t xr[3] = {xmid[0]+dxnew2,xmid[1]+dxnew2,xmid[2]+dxnew2};
real_t xc[8][3] =
{
{xl[0],xl[1],xl[2]},
{xl[0],xl[1],xr[2]},
{xl[0],xr[1],xl[2]},
{xl[0],xr[1],xr[2]},
{xr[0],xl[1],xl[2]},
{xr[0],xl[1],xr[2]},
{xr[0],xr[1],xl[2]},
{xr[0],xr[1],xr[2]},
};
real_t rr2, res[8], ressum = 0.;
for( int i=0; i<8; ++i )
{
rr2 = sqr(xc[i][0])+sqr(xc[i][1])+sqr(xc[i][2]);
res[i] = tfr->compute_real(rr2)*dV;
if( res[i] != res[i] ){ LOGERR("NaN encountered at r=%f, dx=%f, dV=%f : TF = %f",sqrt(rr2),dx,dV,res[i]); abort(); }
ressum += res[i];
}
real_t ae = fabs((prevval-ressum));
real_t re = fabs(ae/ressum);
if( ae < abs_err || re < rel_err ) return ressum;
if( nsplit > nmaxsplits )
{
LOGWARN("reached maximum number of supdivisions in eval_split_recurse. Ending recursion... : abs. err.=%f, rel. err.=%f",ae, re);
return ressum;
}
//otherwise keep splitting
ressum = 0;
for( int i=0; i<8; ++i )
ressum += eval_split_recurse( tfr, xc[i], dxnew, res[i], nsplit+1 );
return ressum;
}
template< typename real_t >
inline real_t eval_split_recurse( const TransferFunction_real *tfr, real_t *xmid, real_t dx, int nsplit = 0 )
{
//sqr(xmid[0])+sqr(xmid[1])+sqr(xmid[2])
real_t rr2 = sqr(xmid[0]) + sqr(xmid[1]) + sqr(xmid[2]);
real_t prevval = tfr->compute_real(rr2)*dx*dx*dx;
return eval_split_recurse( tfr, xmid, dx, prevval, nsplit );
}
//#define OLD_KERNEL_SAMPLING
template< typename real_t >
void kernel_real_cached<real_t>::precompute_kernel( transfer_function* ptf, tf_type type, const refinement_hierarchy& refh )
@ -605,6 +673,8 @@ namespace convolution{
const double rf8 = pow(ref_fac,3);
const double dx05 = 0.5*dx, dx025 = 0.25*dx;
std::cerr << ">>>>>>>>>>>> " << ref_fac << " <<<<<<<<<<<<<<<<" << std::endl;
if( bperiodic )
{
@ -615,7 +685,7 @@ namespace convolution{
for( int k=0; k<=nz/2; ++k )
{
int iix(i), iiy(j), iiz(k);
double rr[3], rr2;
double rr[3];
if( iix > (int)nx/2 ) iix -= nx;
if( iiy > (int)ny/2 ) iiy -= ny;
@ -668,16 +738,14 @@ namespace convolution{
{
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;
val += tfr->compute_real(rrr2[0]+rrr2[1]+rrr2[2])/rf8;
}
}
}
}
else
{
rr2 = rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2];
val += tfr->compute_real(rr2);
val += tfr->compute_real(rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2]);
}
@ -698,7 +766,7 @@ namespace convolution{
for( int k=0; k<nz; ++k )
{
int iix(i), iiy(j), iiz(k);
double rr[3], rr2;
double rr[3];
if( iix > (int)nx/2 ) iix -= nx;
if( iiy > (int)ny/2 ) iiy -= ny;
@ -732,6 +800,8 @@ namespace convolution{
double val = 0.0;//(fftw_real)tfr->compute_real(rr2)*fac;
#ifdef OLD_KERNEL_SAMPLING
if( ref_fac > 0 )
{
double rrr[3];
@ -748,18 +818,25 @@ namespace convolution{
{
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;
val += tfr->compute_real(rrr2[0]+rrr2[1]+rrr2[2])/rf8;
}
}
}
}
else
{
rr2 = rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2];
val = tfr->compute_real(rr2);
val = tfr->compute_real(rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2]);
}
#else
if( i == 0 && j == 0 && k == 0 ) continue;
// use new exact volume integration scheme
real_t xmid[3] = { rr[0], rr[1], rr[2] };
real_t ddx = dx;
val = eval_split_recurse( tfr, xmid, ddx ) / (ddx*ddx*ddx);
#endif
//if( rr2 <= boxlength2*boxlength2 )
//rkernel[idx] += (fftw_real)tfr->compute_real(rr2)*fac;
@ -771,9 +848,16 @@ namespace convolution{
}
}
{
#ifdef OLD_KERNEL_SAMPLING
rkernel[0] = tfr->compute_real(0.0)*fac;
#else
real_t xmid[3] = {0.0,0.0,0.0};
real_t ddx = dx;
rkernel[0] = fac*eval_split_recurse( tfr, xmid, ddx ) / (ddx*ddx*ddx);
deconv = false;
#endif
}
/*************************************************************************************/
/*************************************************************************************/
/******* perform deconvolution *******************************************************/
@ -1096,11 +1180,26 @@ namespace convolution{
rr[1] = ((double)iiy ) * dxc;
rr[2] = ((double)iiz ) * dxc;
#ifdef OLD_KERNEL_SAMPLING
rkernel_coarse[idx] = 0.0;
rr2 = rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2];
if( fabs(rr[0])<=boxlength2||fabs(rr[1])<=boxlength2||fabs(rr[2])<=boxlength2 )
rkernel_coarse[idx] += (fftw_real)tfr->compute_real(rr2)*fac;
#else
rkernel_coarse[idx] = 0.0;
//if( i==0 && j==0 && k==0 ) continue;
real_t xmid[3] = { rr[0], rr[1], rr[2] };
real_t ddx = dxc;
real_t val = eval_split_recurse( tfr, xmid, ddx ) / (ddx*ddx*ddx);
if( fabs(rr[0])<=boxlength2||fabs(rr[1])<=boxlength2||fabs(rr[2])<=boxlength2 )
rkernel_coarse[idx] += val * fac;
#endif
}
}
@ -1108,7 +1207,7 @@ namespace convolution{
LOGUSER("Averaging fine kernel to coarse kernel...");
//... copy averaged and convolved fine kernel to coarse kernel
#pragma omp parallel for
/*#pragma omp parallel for
for( int ix=0; ix<nx; ix+=2 )
for( int iy=0; iy<ny; iy+=2 )
for( int iz=0; iz<nz; iz+=2 )
@ -1140,7 +1239,7 @@ namespace convolution{
+rkernel[ACC_RF(ix-i,iy-j+1,iz-k+1)] + rkernel[ACC_RF(ix-i+1,iy-j+1,iz-k+1)]);
}
}
}*/
sprintf(cachefname,"temp_kernel_level%03d.tmp",ilevel);

7
densities.cc
View file

@ -615,7 +615,7 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type
coarse->subtract_oct_mean();
convolution::perform<real_t>( the_tf_kernel, reinterpret_cast<void*> (coarse->get_data_ptr()), shift );
coarse->subtract_mean();
coarse->upload_bnd_add( *delta.get_grid(levelmin+i-1) );
//coarse->upload_bnd_add( *delta.get_grid(levelmin+i-1) );
//... clean up
the_tf_kernel->deallocate();
@ -667,7 +667,7 @@ void GenerateDensityHierarchy( config_file& cf, transfer_function *ptf, tf_type
coarse->subtract_mean();
//... upload data to coarser grid
coarse->upload_bnd_add( *delta.get_grid(levelmax-1) );
//coarse->upload_bnd_add( *delta.get_grid(levelmax-1) );
delete coarse;
}
@ -739,6 +739,9 @@ void coarsen_density( const refinement_hierarchy& rh, GridHierarchy<real_t>& u )
{
unsigned levelmin_TF = rh.levelmin();
/* for( int i=rh.levelmax(); i>0; --i )
mg_straight().restrict( *(u.get_grid(i)), *(u.get_grid(i-1)) );*/
//for( unsigned i=levelmin_TF+1; i<=rh.levelmax(); ++i )
for( unsigned i=1; i<=rh.levelmax(); ++i )
{

192
random.cc
View file

@ -1375,6 +1375,191 @@ 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 )
{
@ -1489,7 +1674,7 @@ void random_number_generator<rng,T>::compute_random_numbers( void )
//... apply constraints to this level, if any
//if( ilevel == levelmax_ )
constraints.apply( ilevel, x0, lx, randc[ilevel] );
//constraints.apply( ilevel, x0, lx, randc[ilevel] );
//... store numbers
store_rnd( ilevel, randc[ilevel] );
@ -1498,6 +1683,11 @@ 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_ )
{

14
transfer_function.hh
View file

@ -532,7 +532,7 @@ public:
for( unsigned i=0; i<r.size(); ++i )
{
if( r[i] > rmin && r[i] < rmax )
if( r[i] > rmin/512. && r[i] < rmax )
{
xsp.push_back( log10(r[i]) );
ysp.push_back( T[i]*r[i]*r[i] );
@ -676,7 +676,17 @@ public:
y2 = m_ytable[i+1];
//divide by r**2 because r^2 T is tabulated
return (real_t)((y1 + (y2-y1)*(ii-(double)i))/r2);
//return (real_t)((y1 + (y2-y1)*(ii-(double)i))/r2);
real_t retval = (real_t)((y1 + (y2-y1)*(ii-(double)i))/r2);
if( retval != retval ){
std::cerr << "FAILURE FAILURE FAILURE" << std::endl;
fprintf(stderr,"r2 = %f, r = %f, i = %d, y1 = %f, y2 = %f, xtable[i]=%f",r2,r,i,y1,y2, m_xtable[i]);
abort();
}
return retval;
}
};