2010-07-02 20:49:30 +02:00
|
|
|
/*
|
|
|
|
|
|
|
|
convolution_kernel.cc - This file is part of MUSIC -
|
|
|
|
a code to generate multi-scale initial conditions
|
|
|
|
for cosmological simulations
|
|
|
|
|
|
|
|
Copyright (C) 2010 Oliver Hahn
|
|
|
|
|
|
|
|
This program is free software: you can redistribute it and/or modify
|
|
|
|
it under the terms of the GNU General Public License as published by
|
|
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
|
|
(at your option) any later version.
|
|
|
|
|
|
|
|
This program is distributed in the hope that it will be useful,
|
|
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
|
|
GNU General Public License for more details.
|
|
|
|
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
|
|
along with this program. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
|
|
|
|
*/
|
|
|
|
|
|
|
|
#ifdef SINGLE_PRECISION
|
|
|
|
#ifdef SINGLETHREAD_FFTW
|
|
|
|
#include <srfftw.h>
|
|
|
|
#else
|
|
|
|
#include <srfftw_threads.h>
|
|
|
|
#endif
|
|
|
|
#else
|
|
|
|
#ifdef SINGLETHREAD_FFTW
|
|
|
|
#include <drfftw.h>
|
|
|
|
#else
|
|
|
|
#include <drfftw_threads.h>
|
|
|
|
#endif
|
|
|
|
#endif
|
|
|
|
|
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
|
2010-07-02 20:49:30 +02:00
|
|
|
#include "densities.hh"
|
|
|
|
#include "convolution_kernel.hh"
|
|
|
|
|
2010-07-28 00:21:50 +02:00
|
|
|
double T0 = 1.0;
|
2010-07-02 20:49:30 +02:00
|
|
|
|
|
|
|
namespace convolution{
|
|
|
|
|
|
|
|
std::map< std::string, kernel_creator *>&
|
|
|
|
get_kernel_map()
|
|
|
|
{
|
|
|
|
static std::map< std::string, kernel_creator* > kernel_map;
|
|
|
|
return kernel_map;
|
|
|
|
}
|
|
|
|
|
|
|
|
template< typename real_t >
|
|
|
|
void perform( kernel * pk, void *pd )
|
|
|
|
{
|
2010-09-08 10:06:18 +02:00
|
|
|
//return;
|
2010-07-02 20:49:30 +02:00
|
|
|
parameters cparam_ = pk->cparam_;
|
|
|
|
double fftnorm = pow(2.0*M_PI,1.5)/sqrt(cparam_.lx*cparam_.ly*cparam_.lz)/sqrt((double)(cparam_.nx*cparam_.ny*cparam_.nz));
|
|
|
|
|
|
|
|
fftw_complex *cdata,*ckernel;
|
|
|
|
fftw_real *data;
|
|
|
|
|
|
|
|
data = reinterpret_cast<fftw_real*>(pd);
|
|
|
|
cdata = reinterpret_cast<fftw_complex*>(data);
|
|
|
|
ckernel = reinterpret_cast<fftw_complex*>( pk->get_ptr() );
|
|
|
|
|
|
|
|
rfftwnd_plan iplan, plan;
|
|
|
|
|
2010-07-31 14:38:55 +02:00
|
|
|
std::cout << " - Performing density convolution... ("
|
|
|
|
<< cparam_.nx << ", " << cparam_.ny << ", " << cparam_.nz << ")\n";
|
2010-07-02 20:49:30 +02:00
|
|
|
|
|
|
|
|
2010-09-30 00:42:07 +02:00
|
|
|
LOGUSER("Performing kernel convolution on (%5d,%5d,%5d) grid",cparam_.nx ,cparam_.ny ,cparam_.nz );
|
|
|
|
|
2010-07-02 20:49:30 +02:00
|
|
|
plan = rfftw3d_create_plan( cparam_.nx, cparam_.ny, cparam_.nz,
|
|
|
|
FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE|FFTW_IN_PLACE);
|
|
|
|
|
|
|
|
iplan = rfftw3d_create_plan( cparam_.nx, cparam_.ny, cparam_.nz,
|
|
|
|
FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE|FFTW_IN_PLACE);
|
|
|
|
|
|
|
|
|
2010-09-30 00:42:07 +02:00
|
|
|
LOGUSER("Performing forward FFT...");
|
2010-07-02 20:49:30 +02:00
|
|
|
#ifndef SINGLETHREAD_FFTW
|
|
|
|
rfftwnd_threads_one_real_to_complex( omp_get_max_threads(), plan, data, NULL );
|
|
|
|
#else
|
|
|
|
rfftwnd_one_real_to_complex( plan, data, NULL );
|
|
|
|
#endif
|
|
|
|
|
|
|
|
#pragma omp parallel for
|
|
|
|
for( int i=0; i<cparam_.nx; ++i )
|
|
|
|
for( int j=0; j<cparam_.ny; ++j )
|
|
|
|
for( int k=0; k<cparam_.nz/2+1; ++k )
|
|
|
|
{
|
|
|
|
unsigned ii = (i*cparam_.ny + j) * (cparam_.nz/2+1) + k;
|
2010-09-01 06:59:31 +02:00
|
|
|
|
2010-08-22 04:43:21 +02:00
|
|
|
std::complex<double> ccdata(cdata[ii].re,cdata[ii].im), cckernel(ckernel[ii].re,ckernel[ii].im);
|
2010-07-21 10:10:29 +02:00
|
|
|
ccdata = ccdata * cckernel *fftnorm;
|
2010-09-01 06:59:31 +02:00
|
|
|
|
2010-07-02 20:49:30 +02:00
|
|
|
cdata[ii].re = ccdata.real();
|
|
|
|
cdata[ii].im = ccdata.imag();
|
|
|
|
}
|
2010-09-30 00:42:07 +02:00
|
|
|
|
|
|
|
LOGUSER("Performing backward FFT...");
|
2010-07-02 20:49:30 +02:00
|
|
|
#ifndef SINGLETHREAD_FFTW
|
|
|
|
rfftwnd_threads_one_complex_to_real( omp_get_max_threads(), iplan, cdata, NULL);
|
|
|
|
#else
|
|
|
|
rfftwnd_one_complex_to_real(iplan, cdata, NULL);
|
|
|
|
#endif
|
|
|
|
|
|
|
|
rfftwnd_destroy_plan(plan);
|
|
|
|
rfftwnd_destroy_plan(iplan);
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
template void perform<double>( kernel* pk, void *pd );
|
|
|
|
template void perform<float>( kernel* pk, void *pd );
|
|
|
|
|
|
|
|
/*****************************************************************************************\
|
|
|
|
*** SPECIFIC KERNEL IMPLEMENTATIONS *********************************************
|
|
|
|
\*****************************************************************************************/
|
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
2010-07-02 20:49:30 +02:00
|
|
|
template< typename real_t >
|
|
|
|
class kernel_real : public kernel
|
|
|
|
{
|
|
|
|
protected:
|
|
|
|
std::vector<real_t> kdata_;
|
|
|
|
|
2010-08-22 04:43:21 +02:00
|
|
|
void compute_kernel( tf_type type );
|
2010-07-02 20:49:30 +02:00
|
|
|
|
|
|
|
public:
|
2010-09-25 01:14:09 +02:00
|
|
|
kernel_real( config_file& cf, transfer_function* ptf, refinement_hierarchy& refh, tf_type type )
|
|
|
|
: kernel( cf, ptf, refh, type )
|
|
|
|
{ }
|
2010-07-02 20:49:30 +02:00
|
|
|
|
|
|
|
void *get_ptr()
|
|
|
|
{ return reinterpret_cast<void*> (&kdata_[0]); }
|
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
kernel* fetch_kernel( int ilevel, bool isolated=false );
|
|
|
|
|
2010-07-02 20:49:30 +02:00
|
|
|
~kernel_real()
|
2010-09-25 01:14:09 +02:00
|
|
|
{ deallocate(); }
|
|
|
|
|
|
|
|
void deallocate()
|
2010-07-02 20:49:30 +02:00
|
|
|
{ std::vector<real_t>().swap( kdata_ ); }
|
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
|
|
|
|
|
2010-07-02 20:49:30 +02:00
|
|
|
};
|
|
|
|
|
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
|
|
|
|
|
2010-07-02 20:49:30 +02:00
|
|
|
template< typename real_t >
|
2010-09-25 01:14:09 +02:00
|
|
|
kernel* kernel_real<real_t>::fetch_kernel( int ilevel, bool isolated )
|
2010-07-02 20:49:30 +02:00
|
|
|
{
|
2010-09-25 01:14:09 +02:00
|
|
|
double
|
|
|
|
boxlength = pcf_->getValue<double>("setup","boxlength"),
|
|
|
|
boxlength2 = 0.5*boxlength;
|
|
|
|
|
|
|
|
int
|
|
|
|
levelmax = prefh_->levelmax(),
|
|
|
|
levelmin = prefh_->levelmin(),
|
2010-09-29 00:58:41 +02:00
|
|
|
nx = prefh_->size(ilevel,0),
|
|
|
|
ny = prefh_->size(ilevel,1),
|
|
|
|
nz = prefh_->size(ilevel,2);
|
2010-09-25 01:14:09 +02:00
|
|
|
|
2010-09-29 00:58:41 +02:00
|
|
|
/*if( isolated )
|
|
|
|
{ nx *= 2; ny *= 2; nz *= 2; }*/
|
2010-09-25 01:14:09 +02:00
|
|
|
if( ilevel != levelmin )
|
2010-09-29 00:58:41 +02:00
|
|
|
{ nx*=2; ny*=2; nz*=2; }
|
2010-09-25 01:14:09 +02:00
|
|
|
|
|
|
|
kdata_.assign( nx*ny*2*(nz/2+1), 0.0 );
|
|
|
|
|
|
|
|
double
|
2010-09-29 00:58:41 +02:00
|
|
|
dx = boxlength / (1<<ilevel),
|
2010-09-25 01:14:09 +02:00
|
|
|
lx = dx * nx,
|
|
|
|
ly = dx * ny,
|
|
|
|
lz = dx * nz;
|
|
|
|
|
2010-07-02 20:49:30 +02:00
|
|
|
double
|
2010-09-25 01:14:09 +02:00
|
|
|
kny = nx*M_PI/lx,
|
|
|
|
fac = lx*ly*lz/pow(2.0*M_PI,3)/(nx*ny*nz),
|
|
|
|
nspec = pcf_->getValue<double>("cosmology","nspec"),
|
|
|
|
pnorm = pcf_->getValue<double>("cosmology","pnorm"),
|
|
|
|
dplus = pcf_->getValue<double>("cosmology","dplus");
|
2010-07-02 20:49:30 +02:00
|
|
|
|
|
|
|
bool
|
2010-09-25 01:14:09 +02:00
|
|
|
bavgfine = pcf_->getValueSafe<bool>("setup","avg_fine",true),
|
|
|
|
bperiodic = pcf_->getValueSafe<bool>("setup","periodic_TF",true),
|
2010-09-29 00:58:41 +02:00
|
|
|
kspacepoisson = (pcf_->getValueSafe<bool>("poisson","fft_fine",true)
|
2010-10-05 01:34:01 +02:00
|
|
|
|pcf_->getValueSafe<bool>("poisson","kspace",false));
|
2010-09-25 01:14:09 +02:00
|
|
|
|
|
|
|
cparam_.nx = nx;
|
|
|
|
cparam_.ny = ny;
|
|
|
|
cparam_.nz = nz;
|
|
|
|
cparam_.lx = dx * cparam_.nx;
|
|
|
|
cparam_.ly = dx * cparam_.ny;
|
|
|
|
cparam_.lz = dx * cparam_.nz;
|
|
|
|
cparam_.pcf = pcf_;
|
|
|
|
|
2010-09-29 00:58:41 +02:00
|
|
|
if( pcf_->getValueSafe<bool>("setup","deconvolve",true) )
|
|
|
|
cparam_.deconvolve = true;
|
|
|
|
else
|
|
|
|
cparam_.deconvolve = false;
|
2010-08-22 04:43:21 +02:00
|
|
|
|
2010-09-29 00:58:41 +02:00
|
|
|
cparam_.coarse_fact = (prefh_->levelmax()-ilevel);
|
|
|
|
|
2010-08-22 04:43:21 +02:00
|
|
|
if( bavgfine )
|
|
|
|
cparam_.coarse_fact++;
|
2010-09-01 06:59:31 +02:00
|
|
|
|
2010-09-29 00:58:41 +02:00
|
|
|
const int ref_fac = 1<<cparam_.coarse_fact, ql = -ref_fac/2+1, qr=ql+ref_fac, rf8=(int)pow(ref_fac,3);
|
2010-08-05 21:14:41 +02:00
|
|
|
|
2010-07-21 10:10:29 +02:00
|
|
|
std::cout << " - Computing transfer function kernel...\n";
|
2010-07-02 20:49:30 +02:00
|
|
|
|
2010-07-21 10:10:29 +02:00
|
|
|
if(! cparam_.is_finest )
|
2010-09-01 06:59:31 +02:00
|
|
|
kny *= pow(2,cparam_.coarse_fact);
|
|
|
|
|
2010-07-21 10:10:29 +02:00
|
|
|
|
|
|
|
TransferFunction_real *tfr =
|
2010-10-05 01:34:01 +02:00
|
|
|
new TransferFunction_real( boxlength, 1<<levelmax, type_, ptf_,nspec,pnorm,dplus,
|
2010-09-25 01:14:09 +02:00
|
|
|
0.25*lx/(double)nx,2.0*boxlength,kny, (int)pow(2,levelmax+2));
|
2010-07-21 10:10:29 +02:00
|
|
|
|
|
|
|
fftw_real *rkernel = reinterpret_cast<fftw_real*>( &kdata_[0] );
|
|
|
|
fftw_complex *kkernel = reinterpret_cast<fftw_complex*>( &kdata_[0] );
|
|
|
|
|
|
|
|
rfftwnd_plan plan = rfftw3d_create_plan( cparam_.nx, cparam_.ny, cparam_.nz,
|
|
|
|
FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE|FFTW_IN_PLACE);
|
|
|
|
|
|
|
|
int nzp = cparam_.nz/2+1;
|
|
|
|
double kmax = 0.5*M_PI/std::max(cparam_.nx,std::max(cparam_.ny,cparam_.nz));
|
|
|
|
|
2010-08-15 01:52:25 +02:00
|
|
|
T0 = tfr->compute_real(0.0);
|
|
|
|
|
2010-09-01 06:59:31 +02:00
|
|
|
|
|
|
|
|
|
|
|
if(bavgfine)
|
|
|
|
cparam_.coarse_fact--;
|
2010-09-25 01:14:09 +02:00
|
|
|
|
2010-07-31 14:38:55 +02:00
|
|
|
//... obtain a grid-version of the real space transfer function
|
2010-07-21 10:10:29 +02:00
|
|
|
if( bperiodic )
|
2010-08-15 01:52:25 +02:00
|
|
|
{
|
|
|
|
#pragma omp parallel for //schedule(static,4)
|
2010-09-25 01:14:09 +02:00
|
|
|
for( int i=0; i<=nx/2; ++i )
|
|
|
|
for( int j=0; j<=ny/2; ++j )
|
|
|
|
for( int k=0; k<=nz/2; ++k )
|
2010-07-02 20:49:30 +02:00
|
|
|
{
|
|
|
|
int iix(i), iiy(j), iiz(k);
|
|
|
|
double rr[3], rr2;
|
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
if( iix > (int)nx/2 ) iix -= nx;
|
|
|
|
if( iiy > (int)ny/2 ) iiy -= ny;
|
|
|
|
if( iiz > (int)nz/2 ) iiz -= nz;
|
2010-07-02 20:49:30 +02:00
|
|
|
|
2010-08-15 01:52:25 +02:00
|
|
|
|
|
|
|
//... speed up 8x by copying data to other octants
|
|
|
|
int idx[8];
|
2010-09-25 01:14:09 +02:00
|
|
|
//int nx=cparam_.nx, ny = cparam_.ny, nz = cparam_.nz;
|
2010-08-15 01:52:25 +02:00
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
idx[0] = ((i)*ny + (j)) * 2*(nz/2+1) + (k);
|
|
|
|
idx[1] = ((nx-i)*ny + (j)) * 2*(nz/2+1) + (k);
|
|
|
|
idx[2] = ((i)*ny + (ny-j)) * 2*(nz/2+1) + (k);
|
|
|
|
idx[3] = ((nx-i)*ny + (ny-j)) * 2*(nz/2+1) + (k);
|
|
|
|
idx[4] = ((i)*ny + (j)) * 2*(nz/2+1) + (nz-k);
|
|
|
|
idx[5] = ((nx-i)*ny + (j)) * 2*(nz/2+1) + (nz-k);
|
|
|
|
idx[6] = ((i)*ny + (ny-j)) * 2*(nz/2+1) + (nz-k);
|
|
|
|
idx[7] = ((nx-i)*ny + (ny-j)) * 2*(nz/2+1) + (nz-k);
|
2010-08-15 01:52:25 +02:00
|
|
|
|
|
|
|
if(i==0||i==nx/2){ idx[1]=idx[3]=idx[5]=idx[7]=-1;}
|
|
|
|
if(j==0||j==ny/2){ idx[2]=idx[3]=idx[6]=idx[7]=-1;}
|
|
|
|
if(k==0||k==nz/2){ idx[4]=idx[5]=idx[6]=idx[7]=-1;}
|
|
|
|
|
|
|
|
double val = 0.0;
|
2010-07-02 20:49:30 +02:00
|
|
|
|
|
|
|
for( int ii=-1; ii<=1; ++ii )
|
|
|
|
for( int jj=-1; jj<=1; ++jj )
|
|
|
|
for( int kk=-1; kk<=1; ++kk )
|
|
|
|
{
|
|
|
|
rr[0] = ((double)iix ) * dx + ii*boxlength;
|
2010-09-25 01:14:09 +02:00
|
|
|
rr[1] = ((double)iiy ) * dx + jj*boxlength;
|
|
|
|
rr[2] = ((double)iiz ) * dx + kk*boxlength;
|
2010-07-02 20:49:30 +02:00
|
|
|
|
2010-08-15 01:52:25 +02:00
|
|
|
if( rr[0] > -boxlength && rr[0] < boxlength
|
|
|
|
&& rr[1] > -boxlength && rr[1] < boxlength
|
|
|
|
&& rr[2] > -boxlength && rr[2] < boxlength )
|
2010-07-02 20:49:30 +02:00
|
|
|
{
|
2010-08-15 01:52:25 +02:00
|
|
|
if( ref_fac > 1 )//&& fabs(tfr->get_grad(rr2)*dx/T0) > 1e-4 )
|
2010-08-05 21:14:41 +02:00
|
|
|
{
|
|
|
|
double rrr[3];
|
|
|
|
register double rrr2[3];
|
|
|
|
for( int iii=ql; iii<qr; ++iii )
|
|
|
|
{
|
|
|
|
rrr[0] = rr[0]+(double)iii*0.5*dx - 0.25*dx;
|
|
|
|
rrr2[0]= rrr[0]*rrr[0];
|
|
|
|
for( int jjj=ql; jjj<qr; ++jjj )
|
|
|
|
{
|
|
|
|
rrr[1] = rr[1]+(double)jjj*0.5*dx - 0.25*dx;
|
|
|
|
rrr2[1]= rrr[1]*rrr[1];
|
|
|
|
for( int kkk=ql; kkk<qr; ++kkk )
|
|
|
|
{
|
|
|
|
rrr[2] = rr[2]+(double)kkk*0.5*dx - 0.25*dx;
|
|
|
|
rrr2[2]= rrr[2]*rrr[2];
|
|
|
|
rr2 = rrr2[0]+rrr2[1]+rrr2[2];
|
2010-08-15 01:52:25 +02:00
|
|
|
val += tfr->compute_real(rr2)/rf8;
|
2010-08-05 21:14:41 +02:00
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
}else{
|
|
|
|
rr2 = rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2];
|
2010-08-15 01:52:25 +02:00
|
|
|
val += tfr->compute_real(rr2);
|
2010-08-05 21:14:41 +02:00
|
|
|
}
|
2010-07-02 20:49:30 +02:00
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2010-08-15 01:52:25 +02:00
|
|
|
val *= fac;
|
|
|
|
|
|
|
|
for(int q=0;q<8;++q)
|
|
|
|
if(idx[q]!=-1)
|
2010-08-22 04:43:21 +02:00
|
|
|
kdata_[idx[q]] = val;
|
2010-07-02 20:49:30 +02:00
|
|
|
}
|
2010-08-22 04:43:21 +02:00
|
|
|
|
2010-07-02 20:49:30 +02:00
|
|
|
}else{
|
|
|
|
#pragma omp parallel for
|
2010-09-25 01:14:09 +02:00
|
|
|
for( int i=0; i<nx; ++i )
|
|
|
|
for( int j=0; j<ny; ++j )
|
|
|
|
for( int k=0; k<nz; ++k )
|
2010-07-02 20:49:30 +02:00
|
|
|
{
|
|
|
|
int iix(i), iiy(j), iiz(k);
|
|
|
|
double rr[3], rr2;
|
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
if( iix > (int)nx/2 ) iix -= nx;
|
|
|
|
if( iiy > (int)ny/2 ) iiy -= ny;
|
|
|
|
if( iiz > (int)nz/2 ) iiz -= nz;
|
2010-07-02 20:49:30 +02:00
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
unsigned idx = (i*ny + j) * 2*(nz/2+1) + k;
|
2010-07-02 20:49:30 +02:00
|
|
|
|
|
|
|
rr[0] = ((double)iix ) * dx;
|
2010-09-25 01:14:09 +02:00
|
|
|
rr[1] = ((double)iiy ) * dx;
|
|
|
|
rr[2] = ((double)iiz ) * dx;
|
|
|
|
|
|
|
|
kdata_[idx] = 0.0;
|
2010-07-02 20:49:30 +02:00
|
|
|
|
2010-08-05 21:14:41 +02:00
|
|
|
if( ref_fac > 1 )
|
|
|
|
{
|
|
|
|
double rrr[3];
|
|
|
|
register double rrr2[3];
|
|
|
|
for( int iii=ql; iii<qr; ++iii )
|
|
|
|
{
|
|
|
|
rrr[0] = rr[0]+(double)iii*0.5*dx - 0.25*dx;
|
|
|
|
rrr2[0]= rrr[0]*rrr[0];
|
|
|
|
for( int jjj=ql; jjj<qr; ++jjj )
|
|
|
|
{
|
|
|
|
rrr[1] = rr[1]+(double)jjj*0.5*dx - 0.25*dx;
|
|
|
|
rrr2[1]= rrr[1]*rrr[1];
|
|
|
|
for( int kkk=ql; kkk<qr; ++kkk )
|
|
|
|
{
|
|
|
|
rrr[2] = rr[2]+(double)kkk*0.5*dx - 0.25*dx;
|
|
|
|
rrr2[2]= rrr[2]*rrr[2];
|
|
|
|
rr2 = rrr2[0]+rrr2[1]+rrr2[2];
|
2010-09-01 06:59:31 +02:00
|
|
|
if( fabs(rr[0])<=boxlength2||fabs(rr[1])<=boxlength2||fabs(rr[2])<=boxlength2 )
|
|
|
|
kdata_[idx] += (fftw_real)(tfr->compute_real(rr2)/rf8)*fac;
|
2010-08-05 21:14:41 +02:00
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
2010-09-25 01:14:09 +02:00
|
|
|
if( false )//cparam_.deconvolve )//&& cparam_.is_finest )
|
|
|
|
{
|
|
|
|
double rx((double)iix*dx/boxlength),ry((double)iiy*dx/boxlength),rz((double)iiz*dx/boxlength), ico;
|
|
|
|
ico = 1.0;
|
|
|
|
if( i>0 && fabs(dx*(double)iix) < 0.5*boxlength)//&& i<cparam_.nx/2)
|
|
|
|
ico *= cos(M_PI*rx);//sin(M_PI*rx)/(M_PI*rx);
|
|
|
|
if( j>0 && fabs(dx*(double)iiy) < 0.5*boxlength)//&& j<cparam_.ny/2 )
|
|
|
|
ico *= cos(M_PI*ry);//sin(M_PI*ry)/(M_PI*ry);
|
|
|
|
if( k>0 && fabs(dx*(double)iiz) < 0.5*boxlength)//&& k<cparam_.nz/2)
|
|
|
|
ico *= cos(M_PI*rz);//sin(M_PI*rz)/(M_PI*rz);
|
|
|
|
//kdata_[idx] /= ico;
|
|
|
|
//if( fabs(ico)>1e-8 );
|
|
|
|
kdata_[idx] /= ico;
|
|
|
|
|
|
|
|
}
|
2010-08-05 21:14:41 +02:00
|
|
|
}else{
|
|
|
|
rr2 = rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2];
|
2010-09-01 06:59:31 +02:00
|
|
|
if( fabs(rr[0])<=boxlength2||fabs(rr[1])<=boxlength2||fabs(rr[2])<=boxlength2 )
|
|
|
|
kdata_[idx] += (fftw_real)tfr->compute_real(rr2)*fac;
|
2010-09-25 01:14:09 +02:00
|
|
|
|
2010-08-05 21:14:41 +02:00
|
|
|
}
|
2010-07-02 20:49:30 +02:00
|
|
|
}
|
2010-07-21 10:10:29 +02:00
|
|
|
}
|
2010-08-22 04:43:21 +02:00
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
|
|
|
|
if( ilevel == levelmax )
|
2010-08-22 04:43:21 +02:00
|
|
|
kdata_[0] = tfr->compute_real(0.0)*fac;
|
2010-09-29 00:58:41 +02:00
|
|
|
|
2010-07-28 00:21:50 +02:00
|
|
|
T0 = kdata_[0];
|
2010-07-21 10:10:29 +02:00
|
|
|
delete tfr;
|
2010-07-02 20:49:30 +02:00
|
|
|
|
2010-07-21 10:10:29 +02:00
|
|
|
double k0 = kdata_[0];
|
2010-09-01 06:59:31 +02:00
|
|
|
|
2010-08-05 21:14:41 +02:00
|
|
|
//... subtract white noise component before deconvolution
|
|
|
|
if( cparam_.deconvolve )//&& cparam_.is_finest)
|
2010-07-21 10:10:29 +02:00
|
|
|
kdata_[0] = 0.0;
|
2010-07-02 20:49:30 +02:00
|
|
|
|
2010-07-21 10:10:29 +02:00
|
|
|
#ifndef SINGLETHREAD_FFTW
|
|
|
|
rfftwnd_threads_one_real_to_complex( omp_get_max_threads(), plan, rkernel, NULL );
|
|
|
|
#else
|
|
|
|
rfftwnd_one_real_to_complex( plan, rkernel, NULL );
|
|
|
|
#endif
|
2010-09-25 01:14:09 +02:00
|
|
|
|
|
|
|
/*#define XI_SAMPLE
|
|
|
|
#ifdef XI_SAMPLE
|
|
|
|
#pragma omp parallel for
|
|
|
|
for( int i=0; i<cparam_.nx; ++i )
|
|
|
|
for( int j=0; j<cparam_.ny; ++j )
|
|
|
|
for( int k=0; k<nzp; ++k )
|
|
|
|
{
|
|
|
|
unsigned q = (i*cparam_.ny+j)*nzp+k;
|
|
|
|
std::complex<double> z(kkernel[q].re,kkernel[q].im);
|
|
|
|
z=sqrt(z);
|
|
|
|
kkernel[q].re = z.real();
|
|
|
|
kkernel[q].im = z.imag();
|
|
|
|
}
|
|
|
|
|
|
|
|
#endif*/
|
|
|
|
|
|
|
|
|
2010-07-02 20:49:30 +02:00
|
|
|
|
2010-09-01 06:59:31 +02:00
|
|
|
//... do some k-space filter correction stuff
|
2010-07-31 14:38:55 +02:00
|
|
|
if( cparam_.deconvolve || cparam_.smooth )
|
2010-07-02 20:49:30 +02:00
|
|
|
{
|
|
|
|
|
2010-07-21 10:10:29 +02:00
|
|
|
double ksum = 0.0;
|
|
|
|
unsigned kcount = 0;
|
|
|
|
|
|
|
|
#pragma omp parallel for reduction(+:ksum,kcount)
|
2010-09-25 01:14:09 +02:00
|
|
|
for( int i=0; i<nx; ++i )
|
|
|
|
for( int j=0; j<ny; ++j )
|
2010-07-21 10:10:29 +02:00
|
|
|
for( int k=0; k<nzp; ++k )
|
2010-07-02 20:49:30 +02:00
|
|
|
{
|
2010-07-21 10:10:29 +02:00
|
|
|
double kx,ky,kz;
|
|
|
|
|
|
|
|
kx = (double)i;
|
|
|
|
ky = (double)j;
|
|
|
|
kz = (double)k;
|
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
if( kx > nx/2 ) kx -= nx;
|
|
|
|
if( ky > ny/2 ) ky -= ny;
|
2010-07-21 10:10:29 +02:00
|
|
|
|
|
|
|
|
|
|
|
double ipix = 1.0;
|
|
|
|
|
2010-07-31 14:38:55 +02:00
|
|
|
|
2010-08-11 00:51:49 +02:00
|
|
|
double kkmax = kmax;// / pow(2,cparam_.coarse_fact);
|
2010-07-31 14:38:55 +02:00
|
|
|
|
2010-08-11 00:51:49 +02:00
|
|
|
if( true )//cparam_.is_finest||cparam_.smooth )
|
2010-08-05 21:14:41 +02:00
|
|
|
{
|
|
|
|
|
|
|
|
//... deconvolve with grid-cell (NGP) kernel
|
|
|
|
if( i > 0 )
|
|
|
|
ipix /= sin(kx*2.0*kkmax)/(kx*2.0*kkmax);
|
|
|
|
if( j > 0 )
|
|
|
|
ipix /= sin(ky*2.0*kkmax)/(ky*2.0*kkmax);
|
|
|
|
if( k > 0 )
|
|
|
|
ipix /= sin(kz*2.0*kkmax)/(kz*2.0*kkmax);
|
|
|
|
}
|
2010-07-21 10:10:29 +02:00
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
unsigned q = (i*ny+j)*nzp+k;
|
2010-07-21 10:10:29 +02:00
|
|
|
|
2010-07-31 14:38:55 +02:00
|
|
|
if( !cparam_.smooth )
|
|
|
|
{
|
2010-09-01 06:59:31 +02:00
|
|
|
if( cparam_.is_finest )
|
2010-08-22 04:43:21 +02:00
|
|
|
{
|
2010-09-01 06:59:31 +02:00
|
|
|
if( kspacepoisson )
|
|
|
|
{
|
|
|
|
kkernel[q].re /= cos(kx*kkmax)*cos(ky*kkmax)*cos(kz*kkmax);
|
|
|
|
kkernel[q].im /= cos(kx*kkmax)*cos(ky*kkmax)*cos(kz*kkmax);
|
|
|
|
///////kkernel[q].re *= ipix;
|
|
|
|
///////kkernel[q].im *= ipix;
|
2010-09-08 10:06:18 +02:00
|
|
|
//kkernel[q].re *= ipix;
|
|
|
|
//kkernel[q].im *= ipix;
|
2010-09-01 06:59:31 +02:00
|
|
|
}else{
|
|
|
|
|
|
|
|
////////kkernel[q].re /= cos(kx*kkmax)*cos(ky*kkmax)*cos(kz*kkmax);
|
|
|
|
////////kkernel[q].im /= cos(kx*kkmax)*cos(ky*kkmax)*cos(kz*kkmax);
|
|
|
|
|
|
|
|
kkernel[q].re *= ipix;
|
|
|
|
kkernel[q].im *= ipix;
|
|
|
|
}
|
|
|
|
|
2010-08-22 04:43:21 +02:00
|
|
|
}else{
|
|
|
|
kkernel[q].re /= cos(kx*kkmax)*cos(ky*kkmax)*cos(kz*kkmax);
|
|
|
|
kkernel[q].im /= cos(kx*kkmax)*cos(ky*kkmax)*cos(kz*kkmax);
|
2010-09-01 06:59:31 +02:00
|
|
|
////////kkernel[q].re *= ipix;
|
|
|
|
////////kkernel[q].im *= ipix;
|
2010-08-22 04:43:21 +02:00
|
|
|
}
|
2010-09-01 06:59:31 +02:00
|
|
|
|
2010-07-31 14:38:55 +02:00
|
|
|
}else{
|
|
|
|
//... if smooth==true, convolve with
|
|
|
|
//... NGP kernel to get CIC smoothness
|
|
|
|
kkernel[q].re /= ipix;
|
|
|
|
kkernel[q].im /= ipix;
|
|
|
|
}
|
|
|
|
|
|
|
|
//... store k-space average
|
2010-09-25 01:14:09 +02:00
|
|
|
if( k==0 || k==nz/2 )
|
2010-07-21 10:10:29 +02:00
|
|
|
{
|
|
|
|
ksum += kkernel[q].re;
|
|
|
|
kcount++;
|
|
|
|
}else{
|
|
|
|
ksum += 2.0*(kkernel[q].re);
|
|
|
|
kcount+=2;
|
|
|
|
}
|
2010-07-08 07:15:33 +02:00
|
|
|
}
|
2010-07-21 10:10:29 +02:00
|
|
|
|
2010-07-28 00:21:50 +02:00
|
|
|
double dk;
|
|
|
|
|
2010-07-31 14:38:55 +02:00
|
|
|
//... re-add white noise component for finest grid
|
2010-09-25 01:14:09 +02:00
|
|
|
if( ilevel == levelmax )
|
2010-07-28 00:21:50 +02:00
|
|
|
dk = k0-ksum/kcount;
|
|
|
|
else
|
2010-08-05 21:14:41 +02:00
|
|
|
dk = k0-ksum/kcount;
|
2010-07-21 10:10:29 +02:00
|
|
|
|
2010-07-31 14:38:55 +02:00
|
|
|
//... set white noise component to zero if smoothing is enabled
|
|
|
|
if( cparam_.smooth )
|
|
|
|
dk = 0.0;
|
2010-09-25 01:14:09 +02:00
|
|
|
|
|
|
|
//dk = 0.0;
|
2010-09-01 06:59:31 +02:00
|
|
|
|
2010-07-21 10:10:29 +02:00
|
|
|
//... enforce the r=0 component by adjusting the k-space mean
|
|
|
|
#pragma omp parallel for reduction(+:ksum,kcount)
|
2010-09-25 01:14:09 +02:00
|
|
|
for( int i=0; i<nx; ++i )
|
|
|
|
for( int j=0; j<ny; ++j )
|
2010-07-21 10:10:29 +02:00
|
|
|
for( int k=0; k<nzp; ++k )
|
|
|
|
{
|
2010-09-25 01:14:09 +02:00
|
|
|
unsigned q = (i*ny+j)*nzp+k;
|
2010-07-21 10:10:29 +02:00
|
|
|
kkernel[q].re += dk;
|
|
|
|
}
|
|
|
|
}
|
2010-07-02 20:49:30 +02:00
|
|
|
|
|
|
|
rfftwnd_destroy_plan(plan);
|
2010-08-22 04:43:21 +02:00
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
return this;
|
2010-07-02 20:49:30 +02:00
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
template< typename real_t >
|
|
|
|
class kernel_k : public kernel
|
|
|
|
{
|
|
|
|
protected:
|
|
|
|
std::vector<real_t> kdata_;
|
|
|
|
|
2010-08-22 04:43:21 +02:00
|
|
|
void compute_kernel( tf_type type );
|
2010-07-02 20:49:30 +02:00
|
|
|
|
|
|
|
public:
|
2010-09-25 01:14:09 +02:00
|
|
|
kernel_k( config_file& cf, transfer_function* ptf, refinement_hierarchy& refh, tf_type type )
|
|
|
|
: kernel( cf, ptf, refh, type )
|
|
|
|
{ }
|
|
|
|
|
|
|
|
kernel* fetch_kernel( int ilevel, bool isolated=false );
|
2010-07-02 20:49:30 +02:00
|
|
|
|
|
|
|
void *get_ptr()
|
|
|
|
{ return reinterpret_cast<void*> (&kdata_[0]); }
|
|
|
|
|
|
|
|
~kernel_k()
|
2010-09-25 01:14:09 +02:00
|
|
|
{ deallocate(); }
|
|
|
|
|
|
|
|
void deallocate()
|
2010-07-02 20:49:30 +02:00
|
|
|
{ std::vector<real_t>().swap( kdata_ ); }
|
|
|
|
|
|
|
|
};
|
|
|
|
|
|
|
|
template< typename real_t >
|
2010-09-25 01:14:09 +02:00
|
|
|
kernel* kernel_k<real_t>::fetch_kernel( int ilevel, bool isolated )
|
2010-07-02 20:49:30 +02:00
|
|
|
{
|
|
|
|
double
|
2010-09-25 01:14:09 +02:00
|
|
|
boxlength = pcf_->getValue<double>("setup","boxlength"),
|
|
|
|
nspec = pcf_->getValue<double>("cosmology","nspec"),
|
|
|
|
pnorm = pcf_->getValue<double>("cosmology","pnorm"),
|
|
|
|
dplus = pcf_->getValue<double>("cosmology","dplus"),
|
|
|
|
fac = pow(boxlength,3)/pow(2.0*M_PI,3);
|
|
|
|
|
|
|
|
TransferFunction_k *tfk = new TransferFunction_k(type_,ptf_,nspec,pnorm,dplus);
|
2010-07-02 20:49:30 +02:00
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
int
|
|
|
|
nx = prefh_->size(prefh_->levelmax(),0),
|
|
|
|
ny = prefh_->size(prefh_->levelmax(),1),
|
|
|
|
nz = prefh_->size(prefh_->levelmax(),2);
|
|
|
|
|
|
|
|
cparam_.nx = nx;
|
|
|
|
cparam_.ny = ny;
|
|
|
|
cparam_.nz = nz;
|
|
|
|
cparam_.lx = boxlength;
|
|
|
|
cparam_.ly = boxlength;
|
|
|
|
cparam_.lz = boxlength;
|
|
|
|
cparam_.pcf = pcf_;
|
|
|
|
|
|
|
|
kdata_.assign( nx*ny*2*(nz/2+1), 0.0 );
|
2010-07-02 20:49:30 +02:00
|
|
|
|
|
|
|
fftw_complex *kdata = reinterpret_cast<fftw_complex*> ( this->get_ptr() );
|
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
unsigned nzp = (nz/2+1);
|
2010-07-31 14:38:55 +02:00
|
|
|
fac =1.0;
|
2010-07-02 20:49:30 +02:00
|
|
|
|
2010-07-28 00:21:50 +02:00
|
|
|
double kfac = 2.0*M_PI/boxlength, ksum = 0.0;
|
|
|
|
unsigned q=0, kcount = 0;
|
|
|
|
|
|
|
|
#pragma omp parallel for reduction(+:ksum,kcount)
|
2010-09-25 01:14:09 +02:00
|
|
|
for( int i=0; i<nx; ++i )
|
|
|
|
for( int j=0; j<ny; ++j )
|
|
|
|
for( int k=0; k<nz/2+1; ++k )
|
2010-07-02 20:49:30 +02:00
|
|
|
{
|
|
|
|
double kx,ky,kz;
|
|
|
|
|
|
|
|
kx = (double)i;
|
|
|
|
ky = (double)j;
|
|
|
|
kz = (double)k;
|
|
|
|
|
|
|
|
if( kx > nx/2 ) kx -= nx;
|
|
|
|
if( ky > ny/2 ) ky -= ny;
|
|
|
|
|
|
|
|
q = (i*ny+j)*nzp+k;
|
|
|
|
kdata[q].re = fac*tfk->compute(kfac*sqrt(kx*kx+ky*ky+kz*kz));
|
|
|
|
kdata[q].im = 0.0;
|
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
if( k==0 || k==nz/2 )
|
2010-07-28 00:21:50 +02:00
|
|
|
{
|
|
|
|
ksum += kdata[q].re;
|
|
|
|
kcount++;
|
|
|
|
}else{
|
|
|
|
ksum += 2.0*(kdata[q].re);
|
|
|
|
kcount+=2;
|
|
|
|
}
|
|
|
|
|
2010-07-02 20:49:30 +02:00
|
|
|
}
|
|
|
|
|
|
|
|
delete tfk;
|
2010-09-25 01:14:09 +02:00
|
|
|
|
|
|
|
return this;
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
////////////////////////////////////////////////////////////////////////////
|
|
|
|
|
|
|
|
template< typename real_t >
|
|
|
|
class kernel_real_cached : public kernel
|
|
|
|
{
|
|
|
|
protected:
|
|
|
|
std::vector<real_t> kdata_;
|
|
|
|
void precompute_kernel( transfer_function* ptf, tf_type type, const refinement_hierarchy& refh );
|
|
|
|
|
|
|
|
public:
|
|
|
|
kernel_real_cached( config_file& cf, transfer_function* ptf, refinement_hierarchy& refh, tf_type type )
|
|
|
|
: kernel( cf, ptf, refh, type )
|
|
|
|
{
|
|
|
|
precompute_kernel(ptf, type, refh);
|
|
|
|
}
|
|
|
|
|
|
|
|
kernel* fetch_kernel( int ilevel, bool isolated=false );
|
|
|
|
|
|
|
|
void *get_ptr()
|
|
|
|
{ return reinterpret_cast<void*> (&kdata_[0]); }
|
|
|
|
|
|
|
|
~kernel_real_cached()
|
|
|
|
{
|
|
|
|
deallocate();
|
|
|
|
}
|
|
|
|
|
|
|
|
void deallocate()
|
|
|
|
{
|
|
|
|
std::vector<real_t>().swap( kdata_ );
|
|
|
|
}
|
|
|
|
|
|
|
|
};
|
|
|
|
|
|
|
|
template< typename real_t >
|
|
|
|
kernel* kernel_real_cached<real_t>::fetch_kernel( int ilevel, bool isolated )
|
|
|
|
{
|
|
|
|
char cachefname[128];
|
|
|
|
sprintf(cachefname,"temp_kernel_level%03d.tmp",ilevel);
|
|
|
|
FILE *fp = fopen(cachefname,"r");
|
|
|
|
|
|
|
|
|
2010-09-29 00:58:41 +02:00
|
|
|
std::cout << " - Fetching kernel for level " << ilevel << std::endl;
|
2010-09-25 01:14:09 +02:00
|
|
|
|
2010-09-30 00:42:07 +02:00
|
|
|
LOGUSER("Loading kernel for level %3d from file \'%s\'...",ilevel,cachefname);
|
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
if( fp == NULL )
|
2010-09-30 00:42:07 +02:00
|
|
|
{
|
|
|
|
LOGERR("Could not open kernel file \'%s\'.",cachefname);
|
2010-09-25 01:14:09 +02:00
|
|
|
throw std::runtime_error("Internal error: cached convolution kernel does not exist on disk!");
|
2010-09-30 00:42:07 +02:00
|
|
|
}
|
2010-09-25 01:14:09 +02:00
|
|
|
|
|
|
|
unsigned nx,ny,nz;
|
|
|
|
|
|
|
|
fread( reinterpret_cast<void*> (&nx), sizeof(unsigned), 1, fp );
|
|
|
|
fread( reinterpret_cast<void*> (&ny), sizeof(unsigned), 1, fp );
|
|
|
|
fread( reinterpret_cast<void*> (&nz), sizeof(unsigned), 1, fp );
|
|
|
|
|
|
|
|
kdata_.assign(nx*ny*nz,0.0);
|
|
|
|
|
|
|
|
fread( reinterpret_cast<void*>(&kdata_[0]), sizeof(fftw_real), nx*ny*nz, fp );
|
|
|
|
fclose(fp);
|
|
|
|
|
|
|
|
//... set parameters
|
|
|
|
|
|
|
|
double boxlength = pcf_->getValue<double>("setup","boxlength");
|
|
|
|
double dx = boxlength / (1<<ilevel);
|
|
|
|
|
|
|
|
cparam_.nx = nx;
|
|
|
|
cparam_.ny = ny;
|
|
|
|
cparam_.nz = nz-2;
|
|
|
|
cparam_.lx = dx * cparam_.nx;
|
|
|
|
cparam_.ly = dx * cparam_.ny;
|
|
|
|
cparam_.lz = dx * cparam_.nz;
|
|
|
|
cparam_.pcf = pcf_;
|
|
|
|
|
|
|
|
fftw_real *rkernel = reinterpret_cast<fftw_real*>( &kdata_[0] );
|
|
|
|
|
|
|
|
rfftwnd_plan plan = rfftw3d_create_plan( cparam_.nx, cparam_.ny, cparam_.nz,
|
|
|
|
FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE|FFTW_IN_PLACE);
|
|
|
|
|
|
|
|
|
|
|
|
#ifndef SINGLETHREAD_FFTW
|
|
|
|
rfftwnd_threads_one_real_to_complex( omp_get_max_threads(), plan, rkernel, NULL );
|
|
|
|
#else
|
|
|
|
rfftwnd_one_real_to_complex( plan, rkernel, NULL );
|
|
|
|
#endif
|
|
|
|
|
|
|
|
rfftwnd_destroy_plan( plan );
|
|
|
|
return this;
|
2010-07-02 20:49:30 +02:00
|
|
|
}
|
2010-09-25 01:14:09 +02:00
|
|
|
|
|
|
|
|
|
|
|
template< typename real_t >
|
|
|
|
void kernel_real_cached<real_t>::precompute_kernel( transfer_function* ptf, tf_type type, const refinement_hierarchy& refh )
|
|
|
|
{
|
|
|
|
//... compute kernel for finest level
|
|
|
|
int nx,ny,nz;
|
|
|
|
double dx,lx,ly,lz;
|
|
|
|
|
|
|
|
double
|
2010-09-29 00:58:41 +02:00
|
|
|
boxlength = pcf_->getValue<double>("setup","boxlength"),
|
|
|
|
boxlength2 = 0.5*boxlength;
|
2010-09-25 01:14:09 +02:00
|
|
|
|
2010-09-29 00:58:41 +02:00
|
|
|
int
|
|
|
|
levelmax = refh.levelmax(),
|
|
|
|
levelmin = refh.levelmin();
|
2010-09-30 00:42:07 +02:00
|
|
|
|
|
|
|
LOGUSER("Precomputing transfer function kernels...");
|
2010-09-29 00:58:41 +02:00
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
nx = refh.size(refh.levelmax(),0);
|
|
|
|
ny = refh.size(refh.levelmax(),1);
|
|
|
|
nz = refh.size(refh.levelmax(),2);
|
|
|
|
|
|
|
|
if( levelmax != levelmin )
|
|
|
|
{
|
|
|
|
nx *= 2;
|
|
|
|
ny *= 2;
|
|
|
|
nz *= 2;
|
|
|
|
}
|
|
|
|
|
|
|
|
dx = boxlength / (1<<refh.levelmax());
|
|
|
|
lx = dx * nx;
|
|
|
|
ly = dx * ny;
|
|
|
|
lz = dx * nz;
|
|
|
|
|
|
|
|
double
|
|
|
|
kny = M_PI/dx,
|
|
|
|
fac = lx*ly*lz/pow(2.0*M_PI,3)/(nx*ny*nz),
|
|
|
|
|
|
|
|
nspec = pcf_->getValue<double>("cosmology","nspec"),
|
|
|
|
pnorm = pcf_->getValue<double>("cosmology","pnorm"),
|
|
|
|
dplus = pcf_->getValue<double>("cosmology","dplus");
|
|
|
|
|
|
|
|
bool
|
|
|
|
bperiodic = pcf_->getValueSafe<bool>("setup","periodic_TF",true);
|
|
|
|
|
|
|
|
std::cout << " - Computing transfer function kernel...\n";
|
|
|
|
|
|
|
|
TransferFunction_real *tfr =
|
2010-10-05 01:34:01 +02:00
|
|
|
new TransferFunction_real( boxlength, 1<<levelmax, type, ptf,nspec,pnorm,dplus,
|
2010-09-29 00:58:41 +02:00
|
|
|
0.25*dx,2.0*boxlength,kny, (int)pow(2,levelmax+2));
|
2010-09-25 01:14:09 +02:00
|
|
|
|
|
|
|
|
|
|
|
fftw_real *rkernel = new fftw_real[nx*ny*2*(nz/2+1)], *rkernel_coarse;
|
2010-09-29 00:58:41 +02:00
|
|
|
|
|
|
|
for( int i=0; i<nx; ++i )
|
|
|
|
for( int j=0; j<ny; ++j )
|
|
|
|
for( int k=0; k<nz; ++k )
|
|
|
|
{
|
|
|
|
int q=((i)*ny + (j)) * 2*(nz/2+1) + (k);
|
|
|
|
rkernel[q] = 0.0;
|
|
|
|
|
|
|
|
}
|
|
|
|
|
2010-09-30 00:42:07 +02:00
|
|
|
LOGUSER("Computing fine kernel (level %d)...", levelmax);
|
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
if( bperiodic )
|
|
|
|
{
|
2010-09-30 00:42:07 +02:00
|
|
|
|
|
|
|
|
2010-09-29 00:58:41 +02:00
|
|
|
#pragma omp parallel for
|
2010-09-25 01:14:09 +02:00
|
|
|
for( int i=0; i<=nx/2; ++i )
|
|
|
|
for( int j=0; j<=ny/2; ++j )
|
|
|
|
for( int k=0; k<=nz/2; ++k )
|
|
|
|
{
|
|
|
|
int iix(i), iiy(j), iiz(k);
|
|
|
|
double rr[3], rr2;
|
|
|
|
|
|
|
|
if( iix > (int)nx/2 ) iix -= nx;
|
|
|
|
if( iiy > (int)ny/2 ) iiy -= ny;
|
|
|
|
if( iiz > (int)nz/2 ) iiz -= nz;
|
|
|
|
|
|
|
|
|
|
|
|
//... speed up 8x by copying data to other octants
|
|
|
|
int idx[8];
|
|
|
|
|
|
|
|
idx[0] = ((i)*ny + (j)) * 2*(nz/2+1) + (k);
|
|
|
|
idx[1] = ((nx-i)*ny + (j)) * 2*(nz/2+1) + (k);
|
|
|
|
idx[2] = ((i)*ny + (ny-j)) * 2*(nz/2+1) + (k);
|
|
|
|
idx[3] = ((nx-i)*ny + (ny-j)) * 2*(nz/2+1) + (k);
|
|
|
|
idx[4] = ((i)*ny + (j)) * 2*(nz/2+1) + (nz-k);
|
|
|
|
idx[5] = ((nx-i)*ny + (j)) * 2*(nz/2+1) + (nz-k);
|
|
|
|
idx[6] = ((i)*ny + (ny-j)) * 2*(nz/2+1) + (nz-k);
|
|
|
|
idx[7] = ((nx-i)*ny + (ny-j)) * 2*(nz/2+1) + (nz-k);
|
|
|
|
|
|
|
|
if(i==0||i==nx/2){ idx[1]=idx[3]=idx[5]=idx[7]=-1;}
|
|
|
|
if(j==0||j==ny/2){ idx[2]=idx[3]=idx[6]=idx[7]=-1;}
|
|
|
|
if(k==0||k==nz/2){ idx[4]=idx[5]=idx[6]=idx[7]=-1;}
|
|
|
|
|
|
|
|
double val = 0.0;
|
|
|
|
|
|
|
|
for( int ii=-1; ii<=1; ++ii )
|
|
|
|
for( int jj=-1; jj<=1; ++jj )
|
|
|
|
for( int kk=-1; kk<=1; ++kk )
|
|
|
|
{
|
|
|
|
rr[0] = ((double)iix ) * dx + ii*boxlength;
|
|
|
|
rr[1] = ((double)iiy ) * dx + jj*boxlength;
|
|
|
|
rr[2] = ((double)iiz ) * dx + kk*boxlength;
|
|
|
|
|
|
|
|
if( rr[0] > -boxlength && rr[0] < boxlength
|
|
|
|
&& rr[1] > -boxlength && rr[1] < boxlength
|
|
|
|
&& rr[2] > -boxlength && rr[2] < boxlength )
|
|
|
|
{
|
|
|
|
rr2 = rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2];
|
|
|
|
val += tfr->compute_real(rr2);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
val *= fac;
|
|
|
|
|
|
|
|
for(int q=0;q<8;++q)
|
|
|
|
if(idx[q]!=-1)
|
|
|
|
rkernel[idx[q]] = val;
|
|
|
|
}
|
|
|
|
|
|
|
|
}else{
|
2010-09-29 00:58:41 +02:00
|
|
|
#pragma omp parallel for
|
2010-09-25 01:14:09 +02:00
|
|
|
for( int i=0; i<nx; ++i )
|
|
|
|
for( int j=0; j<ny; ++j )
|
|
|
|
for( int k=0; k<nz; ++k )
|
|
|
|
{
|
|
|
|
int iix(i), iiy(j), iiz(k);
|
|
|
|
double rr[3], rr2;
|
|
|
|
|
|
|
|
if( iix > (int)nx/2 ) iix -= nx;
|
|
|
|
if( iiy > (int)ny/2 ) iiy -= ny;
|
|
|
|
if( iiz > (int)nz/2 ) iiz -= nz;
|
|
|
|
|
|
|
|
unsigned idx = (i*ny + j) * 2*(nz/2+1) + k;
|
|
|
|
|
|
|
|
rr[0] = ((double)iix ) * dx;
|
|
|
|
rr[1] = ((double)iiy ) * dx;
|
|
|
|
rr[2] = ((double)iiz ) * dx;
|
|
|
|
|
|
|
|
rkernel[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 )
|
|
|
|
//if( rr2 <= boxlength2*boxlength2 )
|
|
|
|
rkernel[idx] += (fftw_real)tfr->compute_real(rr2)*fac;
|
|
|
|
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
2010-09-29 00:58:41 +02:00
|
|
|
rkernel[0] = tfr->compute_real(0.0)*fac;
|
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
/*************************************************************************************/
|
|
|
|
/*************************************************************************************/
|
|
|
|
/******* perform deconvolution *******************************************************/
|
|
|
|
|
2010-09-25 01:16:31 +02:00
|
|
|
if( pcf_->getValueSafe<bool>("setup","deconvolve",true) )
|
2010-09-25 01:14:09 +02:00
|
|
|
{
|
2010-09-25 01:16:31 +02:00
|
|
|
|
2010-09-30 00:42:07 +02:00
|
|
|
LOGUSER("Deconvolving fine kernel...");
|
2010-09-29 00:58:41 +02:00
|
|
|
std::cout << " - Deconvoling density kernel...\n";
|
|
|
|
|
2010-09-25 01:16:31 +02:00
|
|
|
bool kspacepoisson = (pcf_->getValueSafe<bool>("poisson","fft_fine",true)|
|
|
|
|
pcf_->getValueSafe<bool>("poisson","kspace",false));
|
|
|
|
fftw_complex *kkernel = reinterpret_cast<fftw_complex*>( &rkernel[0] );
|
|
|
|
rfftwnd_plan plan = rfftw3d_create_plan( nx,ny,nz, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE|FFTW_IN_PLACE),
|
|
|
|
iplan = rfftw3d_create_plan( nx,ny,nz, FFTW_COMPLEX_TO_REAL, FFTW_ESTIMATE|FFTW_IN_PLACE);
|
|
|
|
|
|
|
|
double fftnorm = 1.0/(nx*ny*nz);
|
|
|
|
double k0 = rkernel[0];
|
|
|
|
|
|
|
|
//... subtract white noise component before deconvolution
|
2010-09-29 00:58:41 +02:00
|
|
|
rkernel[0] = 0.0;
|
2010-09-25 01:16:31 +02:00
|
|
|
|
|
|
|
#ifndef SINGLETHREAD_FFTW
|
|
|
|
rfftwnd_threads_one_real_to_complex( omp_get_max_threads(), plan, rkernel, NULL );
|
|
|
|
#else
|
|
|
|
rfftwnd_one_real_to_complex( plan, rkernel, NULL );
|
|
|
|
#endif
|
|
|
|
|
|
|
|
{
|
|
|
|
|
|
|
|
double ksum = 0.0;
|
|
|
|
unsigned kcount = 0;
|
|
|
|
double kmax = 0.5*M_PI/std::max(nx,std::max(ny,nz));
|
|
|
|
|
2010-09-29 00:58:41 +02:00
|
|
|
#pragma omp parallel for reduction(+:ksum,kcount)
|
2010-09-25 01:16:31 +02:00
|
|
|
for( int i=0; i<nx; ++i )
|
|
|
|
for( int j=0; j<ny; ++j )
|
|
|
|
for( int k=0; k<nz/2+1; ++k )
|
2010-09-25 01:14:09 +02:00
|
|
|
{
|
2010-09-25 01:16:31 +02:00
|
|
|
double kx,ky,kz;
|
|
|
|
|
|
|
|
kx = (double)i;
|
|
|
|
ky = (double)j;
|
|
|
|
kz = (double)k;
|
|
|
|
|
|
|
|
if( kx > nx/2 ) kx -= nx;
|
|
|
|
if( ky > ny/2 ) ky -= ny;
|
|
|
|
|
|
|
|
|
2010-09-29 00:58:41 +02:00
|
|
|
double kkmax = kmax;
|
2010-09-25 01:16:31 +02:00
|
|
|
unsigned q = (i*ny+j)*(nz/2+1)+k;
|
|
|
|
|
|
|
|
if( true )//!cparam_.smooth )
|
2010-09-25 01:14:09 +02:00
|
|
|
{
|
2010-09-25 01:16:31 +02:00
|
|
|
if( kspacepoisson )
|
|
|
|
{
|
2010-09-29 00:58:41 +02:00
|
|
|
//... Use child average response function to emulate sub-sampling
|
|
|
|
double ipix = cos(kx*kkmax)*cos(ky*kkmax)*cos(kz*kkmax);
|
|
|
|
kkernel[q].re /= ipix;
|
|
|
|
kkernel[q].im /= ipix;
|
2010-09-25 01:16:31 +02:00
|
|
|
}else{
|
|
|
|
|
2010-09-29 00:58:41 +02:00
|
|
|
//... Use piecewise constant response function (NGP-kernel)
|
2010-09-25 01:16:31 +02:00
|
|
|
//... for finite difference methods
|
|
|
|
kkmax = kmax;//*2.0;
|
|
|
|
double ipix = 1.0;
|
2010-09-29 00:58:41 +02:00
|
|
|
if( i > 0 )
|
2010-09-25 01:16:31 +02:00
|
|
|
ipix /= sin(kx*2.0*kkmax)/(kx*2.0*kkmax);
|
2010-09-29 00:58:41 +02:00
|
|
|
if( j > 0 )
|
2010-09-25 01:16:31 +02:00
|
|
|
ipix /= sin(ky*2.0*kkmax)/(ky*2.0*kkmax);
|
2010-09-29 00:58:41 +02:00
|
|
|
if( k > 0 )
|
2010-09-25 01:16:31 +02:00
|
|
|
ipix /= sin(kz*2.0*kkmax)/(kz*2.0*kkmax);
|
|
|
|
|
|
|
|
kkernel[q].re *= ipix;
|
|
|
|
kkernel[q].im *= ipix;
|
|
|
|
}
|
2010-09-25 01:14:09 +02:00
|
|
|
}else{
|
2010-09-25 01:16:31 +02:00
|
|
|
//... if smooth==true, convolve with
|
|
|
|
//... NGP kernel to get CIC smoothness
|
2010-09-25 01:14:09 +02:00
|
|
|
double ipix = 1.0;
|
2010-09-25 01:16:31 +02:00
|
|
|
if( i > 0 )
|
2010-09-25 01:14:09 +02:00
|
|
|
ipix /= sin(kx*2.0*kkmax)/(kx*2.0*kkmax);
|
2010-09-25 01:16:31 +02:00
|
|
|
if( j > 0 )
|
2010-09-25 01:14:09 +02:00
|
|
|
ipix /= sin(ky*2.0*kkmax)/(ky*2.0*kkmax);
|
2010-09-25 01:16:31 +02:00
|
|
|
if( k > 0 )
|
2010-09-25 01:14:09 +02:00
|
|
|
ipix /= sin(kz*2.0*kkmax)/(kz*2.0*kkmax);
|
|
|
|
|
2010-09-25 01:16:31 +02:00
|
|
|
kkernel[q].re /= ipix;
|
|
|
|
kkernel[q].im /= ipix;
|
2010-09-25 01:14:09 +02:00
|
|
|
}
|
|
|
|
|
2010-09-25 01:16:31 +02:00
|
|
|
//... store k-space average
|
|
|
|
if( k==0 || k==nz/2 )
|
|
|
|
{
|
|
|
|
ksum += kkernel[q].re;
|
|
|
|
kcount++;
|
|
|
|
}else{
|
|
|
|
ksum += 2.0*(kkernel[q].re);
|
|
|
|
kcount+=2;
|
|
|
|
}
|
2010-09-25 01:14:09 +02:00
|
|
|
}
|
2010-09-25 01:16:31 +02:00
|
|
|
|
|
|
|
double dk;
|
|
|
|
|
|
|
|
//... re-add white noise component for finest grid
|
|
|
|
dk = k0-ksum/kcount;
|
|
|
|
|
|
|
|
//... set white noise component to zero if smoothing is enabled
|
|
|
|
if( false )//cparam_.smooth )
|
|
|
|
dk = 0.0;
|
|
|
|
|
|
|
|
//... enforce the r=0 component by adjusting the k-space mean
|
|
|
|
#pragma omp parallel for reduction(+:ksum,kcount)
|
|
|
|
for( int i=0; i<nx; ++i )
|
|
|
|
for( int j=0; j<ny; ++j )
|
|
|
|
for( int k=0; k<(nz/2+1); ++k )
|
2010-09-25 01:14:09 +02:00
|
|
|
{
|
2010-09-25 01:16:31 +02:00
|
|
|
unsigned q = (i*ny+j)*(nz/2+1)+k;
|
|
|
|
kkernel[q].re += dk;
|
|
|
|
|
|
|
|
kkernel[q].re *= fftnorm;
|
|
|
|
kkernel[q].im *= fftnorm;
|
2010-09-25 01:14:09 +02:00
|
|
|
}
|
2010-09-25 01:16:31 +02:00
|
|
|
}
|
2010-09-25 01:14:09 +02:00
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
2010-09-25 01:16:31 +02:00
|
|
|
#ifndef SINGLETHREAD_FFTW
|
|
|
|
rfftwnd_threads_one_complex_to_real( omp_get_max_threads(), iplan, kkernel, NULL );
|
|
|
|
#else
|
|
|
|
rfftwnd_one_complex_to_real( iplan, kkernel, NULL );
|
|
|
|
#endif
|
|
|
|
rfftwnd_destroy_plan(plan);
|
|
|
|
rfftwnd_destroy_plan(iplan);
|
2010-09-29 00:58:41 +02:00
|
|
|
}
|
2010-09-25 01:14:09 +02:00
|
|
|
|
|
|
|
/*************************************************************************************/
|
|
|
|
/*************************************************************************************/
|
|
|
|
/*************************************************************************************/
|
|
|
|
|
|
|
|
|
|
|
|
char cachefname[128];
|
|
|
|
sprintf(cachefname,"temp_kernel_level%03d.tmp",levelmax);
|
2010-09-30 00:42:07 +02:00
|
|
|
LOGUSER("Storing kernel in temp file \'%s\'.",cachefname);
|
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
FILE *fp = fopen(cachefname,"w+");
|
|
|
|
unsigned q = nx;
|
|
|
|
fwrite( reinterpret_cast<void*> (&q), sizeof(unsigned), 1, fp );
|
|
|
|
q = ny;
|
|
|
|
fwrite( reinterpret_cast<void*> (&q), sizeof(unsigned), 1, fp );
|
|
|
|
q = 2*(nz/2+1);
|
|
|
|
fwrite( reinterpret_cast<void*> (&q), sizeof(unsigned), 1, fp );
|
|
|
|
fwrite( reinterpret_cast<void*>(&rkernel[0]), sizeof(fftw_real), nx*ny*2*(nz/2+1), fp );
|
|
|
|
fclose(fp);
|
|
|
|
|
|
|
|
//... average and fill for other levels
|
|
|
|
for( int ilevel=levelmax-1; ilevel>=levelmin; ilevel-- )
|
|
|
|
{
|
2010-09-30 00:42:07 +02:00
|
|
|
LOGUSER("Computing coarse kernel (level %d)...",ilevel);
|
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
int nxc, nyc, nzc;
|
|
|
|
double dxc,lxc,lyc,lzc;
|
|
|
|
|
|
|
|
nxc = refh.size(ilevel,0);
|
|
|
|
nyc = refh.size(ilevel,1);
|
|
|
|
nzc = refh.size(ilevel,2);
|
|
|
|
|
|
|
|
if( ilevel!=levelmin )
|
|
|
|
{
|
|
|
|
nxc *= 2;
|
|
|
|
nyc *= 2;
|
|
|
|
nzc *= 2;
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
dxc = boxlength / (1<<ilevel);
|
|
|
|
lxc = dxc * nxc;
|
|
|
|
lyc = dxc * nyc;
|
|
|
|
lzc = dxc * nzc;
|
|
|
|
|
|
|
|
rkernel_coarse = new fftw_real[nxc*nyc*2*(nzc/2+1)];
|
|
|
|
fac = lxc*lyc*lzc/pow(2.0*M_PI,3)/(nxc*nyc*nzc);
|
|
|
|
|
|
|
|
if( bperiodic )
|
|
|
|
{
|
2010-09-29 00:58:41 +02:00
|
|
|
#pragma omp parallel for
|
2010-09-25 01:14:09 +02:00
|
|
|
for( int i=0; i<=nxc/2; ++i )
|
|
|
|
for( int j=0; j<=nyc/2; ++j )
|
|
|
|
for( int k=0; k<=nzc/2; ++k )
|
|
|
|
{
|
|
|
|
int iix(i), iiy(j), iiz(k);
|
|
|
|
double rr[3], rr2;
|
|
|
|
|
|
|
|
if( iix > (int)nxc/2 ) iix -= nxc;
|
|
|
|
if( iiy > (int)nyc/2 ) iiy -= nyc;
|
|
|
|
if( iiz > (int)nzc/2 ) iiz -= nzc;
|
|
|
|
|
|
|
|
//... speed up 8x by copying data to other octants
|
|
|
|
int idx[8];
|
|
|
|
|
|
|
|
idx[0] = ((i)*nyc + (j)) * 2*(nzc/2+1) + (k);
|
|
|
|
idx[1] = ((nxc-i)*nyc + (j)) * 2*(nzc/2+1) + (k);
|
|
|
|
idx[2] = ((i)*nyc + (nyc-j)) * 2*(nzc/2+1) + (k);
|
|
|
|
idx[3] = ((nxc-i)*nyc + (nyc-j)) * 2*(nzc/2+1) + (k);
|
|
|
|
idx[4] = ((i)*nyc + (j)) * 2*(nzc/2+1) + (nzc-k);
|
|
|
|
idx[5] = ((nxc-i)*nyc + (j)) * 2*(nzc/2+1) + (nzc-k);
|
|
|
|
idx[6] = ((i)*nyc + (nyc-j)) * 2*(nzc/2+1) + (nzc-k);
|
|
|
|
idx[7] = ((nxc-i)*nyc + (nyc-j)) * 2*(nzc/2+1) + (nzc-k);
|
|
|
|
|
|
|
|
if(i==0||i==nxc/2){ idx[1]=idx[3]=idx[5]=idx[7]=-1;}
|
|
|
|
if(j==0||j==nyc/2){ idx[2]=idx[3]=idx[6]=idx[7]=-1;}
|
|
|
|
if(k==0||k==nzc/2){ idx[4]=idx[5]=idx[6]=idx[7]=-1;}
|
|
|
|
|
|
|
|
double val = 0.0;
|
|
|
|
|
|
|
|
for( int ii=-1; ii<=1; ++ii )
|
|
|
|
for( int jj=-1; jj<=1; ++jj )
|
|
|
|
for( int kk=-1; kk<=1; ++kk )
|
|
|
|
{
|
|
|
|
rr[0] = ((double)iix ) * dxc + ii*boxlength;
|
|
|
|
rr[1] = ((double)iiy ) * dxc + jj*boxlength;
|
|
|
|
rr[2] = ((double)iiz ) * dxc + kk*boxlength;
|
|
|
|
|
|
|
|
if( rr[0] > -boxlength && rr[0] < boxlength
|
|
|
|
&& rr[1] > -boxlength && rr[1] < boxlength
|
|
|
|
&& rr[2] > -boxlength && rr[2] < boxlength )
|
|
|
|
{
|
|
|
|
rr2 = rr[0]*rr[0]+rr[1]*rr[1]+rr[2]*rr[2];
|
|
|
|
val += tfr->compute_real(rr2);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
val *= fac;
|
|
|
|
|
|
|
|
for(int qq=0;qq<8;++qq)
|
|
|
|
if(idx[qq]!=-1)
|
|
|
|
rkernel_coarse[idx[qq]] = val;
|
|
|
|
}
|
|
|
|
|
|
|
|
}else{
|
2010-09-29 00:58:41 +02:00
|
|
|
#pragma omp parallel for
|
2010-09-25 01:14:09 +02:00
|
|
|
for( int i=0; i<nxc; ++i )
|
|
|
|
for( int j=0; j<nyc; ++j )
|
|
|
|
for( int k=0; k<nzc; ++k )
|
|
|
|
{
|
|
|
|
int iix(i), iiy(j), iiz(k);
|
|
|
|
double rr[3], rr2;
|
|
|
|
|
|
|
|
if( iix > (int)nxc/2 ) iix -= nxc;
|
|
|
|
if( iiy > (int)nyc/2 ) iiy -= nyc;
|
|
|
|
if( iiz > (int)nzc/2 ) iiz -= nzc;
|
|
|
|
|
|
|
|
unsigned idx = (i*nyc + j) * 2*(nzc/2+1) + k;
|
|
|
|
|
|
|
|
rr[0] = ((double)iix ) * dxc;
|
|
|
|
rr[1] = ((double)iiy ) * dxc;
|
|
|
|
rr[2] = ((double)iiz ) * dxc;
|
|
|
|
|
|
|
|
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;
|
|
|
|
|
|
|
|
}
|
|
|
|
}
|
2010-09-29 00:58:41 +02:00
|
|
|
|
2010-09-30 00:42:07 +02:00
|
|
|
LOGUSER("Averaging fine kernel to coarse kernel...");
|
2010-09-29 00:58:41 +02:00
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
//... copy averaged and convolved fine kernel to coarse kernel
|
2010-09-29 00:58:41 +02:00
|
|
|
for( int ix=0; ix<nx; ix+=2 )
|
|
|
|
for( int iy=0; iy<ny; iy+=2 )
|
|
|
|
for( int iz=0; iz<nz; iz+=2 )
|
|
|
|
{
|
|
|
|
int iix(ix/2),iiy(iy/2),iiz(iz/2);
|
|
|
|
if( ix>nx/2 ) iix+=nxc-nx/2;
|
|
|
|
if( iy>ny/2 ) iiy+=nyc-ny/2;
|
|
|
|
if( iz>nz/2 ) iiz+=nzc-nz/2;
|
|
|
|
|
|
|
|
if( ix==nx/2||iy==ny/2||iz==nz/2 ) continue;
|
|
|
|
|
|
|
|
for( int i=0; i<=1; ++i )
|
|
|
|
for( int j=0; j<=1; ++j )
|
|
|
|
for( int k=0; k<=1; ++k )
|
|
|
|
if(i==0&&k==0&&j==0 )
|
|
|
|
rkernel_coarse[ ACC_RC(iix,iiy,iiz) ] =
|
|
|
|
0.125 * ( rkernel[ACC_RF(ix-i,iy-j,iz-k)] + rkernel[ACC_RF(ix-i+1,iy-j,iz-k)]
|
|
|
|
+rkernel[ACC_RF(ix-i,iy-j+1,iz-k)] + rkernel[ACC_RF(ix-i,iy-j,iz-k+1)]
|
|
|
|
+rkernel[ACC_RF(ix-i+1,iy-j+1,iz-k)] + rkernel[ACC_RF(ix-i+1,iy-j,iz-k+1)]
|
|
|
|
+rkernel[ACC_RF(ix-i,iy-j+1,iz-k+1)] + rkernel[ACC_RF(ix-i+1,iy-j+1,iz-k+1)]);
|
|
|
|
|
|
|
|
else
|
|
|
|
{
|
|
|
|
|
|
|
|
rkernel_coarse[ ACC_RC(iix,iiy,iiz) ] +=
|
|
|
|
0.125 * ( rkernel[ACC_RF(ix-i,iy-j,iz-k)] + rkernel[ACC_RF(ix-i+1,iy-j,iz-k)]
|
|
|
|
+rkernel[ACC_RF(ix-i,iy-j+1,iz-k)] + rkernel[ACC_RF(ix-i,iy-j,iz-k+1)]
|
|
|
|
+rkernel[ACC_RF(ix-i+1,iy-j+1,iz-k)] + rkernel[ACC_RF(ix-i+1,iy-j,iz-k+1)]
|
|
|
|
+rkernel[ACC_RF(ix-i,iy-j+1,iz-k+1)] + rkernel[ACC_RF(ix-i+1,iy-j+1,iz-k+1)]);
|
|
|
|
}
|
|
|
|
|
|
|
|
}
|
|
|
|
|
2010-09-25 01:14:09 +02:00
|
|
|
|
|
|
|
sprintf(cachefname,"temp_kernel_level%03d.tmp",ilevel);
|
2010-09-30 00:42:07 +02:00
|
|
|
LOGUSER("Storing kernel in temp file \'%s\'.",cachefname);
|
2010-09-25 01:14:09 +02:00
|
|
|
fp = fopen(cachefname,"w+");
|
|
|
|
q = nxc;
|
|
|
|
fwrite( reinterpret_cast<void*> (&q), sizeof(unsigned), 1, fp );
|
|
|
|
q = nyc;
|
|
|
|
fwrite( reinterpret_cast<void*> (&q), sizeof(unsigned), 1, fp );
|
|
|
|
q = 2*(nzc/2+1);
|
|
|
|
fwrite( reinterpret_cast<void*> (&q), sizeof(unsigned), 1, fp );
|
|
|
|
fwrite( reinterpret_cast<void*>(&rkernel_coarse[0]), sizeof(fftw_real), nxc*nyc*2*(nzc/2+1), fp );
|
|
|
|
fclose(fp);
|
|
|
|
|
|
|
|
delete[] rkernel;
|
|
|
|
|
|
|
|
//... prepare for next iteration
|
|
|
|
nx = nxc;
|
|
|
|
ny = nyc;
|
|
|
|
nz = nzc;
|
|
|
|
lx = lxc;
|
|
|
|
ly = lyc;
|
|
|
|
lz = lzc;
|
|
|
|
dx = dxc;
|
|
|
|
rkernel = rkernel_coarse;
|
|
|
|
}
|
|
|
|
|
|
|
|
//... clean up
|
|
|
|
delete[] rkernel;
|
|
|
|
}
|
|
|
|
|
2010-08-11 05:35:29 +02:00
|
|
|
}
|
|
|
|
|
2010-08-22 04:43:21 +02:00
|
|
|
|
|
|
|
/**************************************************************************************/
|
|
|
|
/**************************************************************************************/
|
|
|
|
|
|
|
|
|
2010-07-02 20:49:30 +02:00
|
|
|
namespace{
|
2010-09-25 01:14:09 +02:00
|
|
|
//convolution::kernel_creator_concrete< convolution::kernel_real<double> > creator_d("tf_kernel_real_double");
|
|
|
|
//convolution::kernel_creator_concrete< convolution::kernel_real<float> > creator_f("tf_kernel_real_float");
|
|
|
|
convolution::kernel_creator_concrete< convolution::kernel_real_cached<double> > creator_d("tf_kernel_real_double");
|
|
|
|
convolution::kernel_creator_concrete< convolution::kernel_real_cached<float> > creator_f("tf_kernel_real_float");
|
2010-07-02 20:49:30 +02:00
|
|
|
convolution::kernel_creator_concrete< convolution::kernel_k<double> > creator_kd("tf_kernel_k_double");
|
|
|
|
convolution::kernel_creator_concrete< convolution::kernel_k<float> > creator_kf("tf_kernel_k_float");
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|