From 31290837183f457567084ee70122e857346025e0 Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Tue, 29 Oct 2013 09:45:47 +0100 Subject: [PATCH] working backup commit to repository. maybe this one should not be used. --- convolution_kernel.cc | 131 ++++++++++++++++++++++++---- densities.cc | 7 +- random.cc | 194 +++++++++++++++++++++++++++++++++++++++++- transfer_function.hh | 14 ++- 4 files changed, 324 insertions(+), 22 deletions(-) diff --git a/convolution_kernel.cc b/convolution_kernel.cc index 84cd7ff..e9fbe19 100644 --- a/convolution_kernel.cc +++ b/convolution_kernel.cc @@ -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::precompute_kernel( transfer_function* ptf, tf_type type, const refinement_hierarchy& refh ) @@ -604,6 +672,8 @@ namespace convolution{ 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; + + 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 (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,19 +818,26 @@ 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; val *= fac; @@ -771,9 +848,16 @@ namespace convolution{ } } - - rkernel[0] = tfr->compute_real(0.0)*fac; - + { +#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; ixsubtract_oct_mean(); convolution::perform( the_tf_kernel, reinterpret_cast (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; } @@ -738,6 +738,9 @@ void normalize_density( grid_hierarchy& delta ) void coarsen_density( const refinement_hierarchy& rh, GridHierarchy& 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 ) diff --git a/random.cc b/random.cc index e530d13..14d834d 100644 --- a/random.cc +++ b/random.cc @@ -1375,6 +1375,191 @@ void random_number_generator::parse_rand_parameters( void ) } +template< typename rng, typename T > +void random_number_generator::correct_avg( int icoarse, int ifine ) +{ + int shift[3], levelmin_poisson; + shift[0] = pcf_->getValue("setup","shift_x"); + shift[1] = pcf_->getValue("setup","shift_y"); + shift[2] = pcf_->getValue("setup","shift_z"); + + levelmin_poisson = pcf_->getValue("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 (&nxf), sizeof(unsigned) ); + iffine.read( reinterpret_cast (&nyf), sizeof(unsigned) ); + iffine.read( reinterpret_cast (&nzf), sizeof(unsigned) ); + + ifcoarse.read( reinterpret_cast (&nxc), sizeof(unsigned) ); + ifcoarse.read( reinterpret_cast (&nyc), sizeof(unsigned) ); + ifcoarse.read( reinterpret_cast (&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 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 fine_rand( 2*nyf*nzf, 0.0 ); + iffine.read( reinterpret_cast (&fine_rand[0]), 2*nyf*nzf*sizeof(T) ); + +#pragma omp parallel for + for( int j=0; j coarse_rand(nxc*nyc*nzc,0.0); + ifcoarse.read( reinterpret_cast (&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= 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 (&nxc), sizeof(unsigned) ); + ofcoarse.write( reinterpret_cast (&nyc), sizeof(unsigned) ); + ofcoarse.write( reinterpret_cast (&nzc), sizeof(unsigned) ); + ofcoarse.write( reinterpret_cast (&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= 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::compute_random_numbers( void ) { @@ -1389,7 +1574,7 @@ void random_number_generator::compute_random_numbers( void ) if( levelmin_seed_ < levelmin_ ) { if( rngfnames_[levelmin_seed_].size() > 0 ) - randc[levelmin_seed_] + randc[levelmin_seed_] = new rng( 1<::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] ); @@ -1497,6 +1682,11 @@ void random_number_generator::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_ ) diff --git a/transfer_function.hh b/transfer_function.hh index 6ada43c..ff053b8 100644 --- a/transfer_function.hh +++ b/transfer_function.hh @@ -532,7 +532,7 @@ public: for( unsigned i=0; 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; } };