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
|
|
|
|
|
|
|
|
|
|
|
|
#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
|
|
|
|
|
|
|
|
|
|
|
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);
|
|
|
|
|
|
|
|
|
|
|
|
#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();
|
|
|
|
}
|
|
|
|
|
|
|
|
#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 *********************************************
|
|
|
|
\*****************************************************************************************/
|
|
|
|
|
|
|
|
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-08-22 04:43:21 +02:00
|
|
|
kernel_real( const parameters& cp, tf_type type )
|
|
|
|
: kernel( cp, type )
|
2010-07-02 20:49:30 +02:00
|
|
|
{
|
|
|
|
kdata_.assign( cparam_.nx*cparam_.ny*2*(cparam_.nz/2+1), 0.0 );
|
2010-08-22 04:43:21 +02:00
|
|
|
compute_kernel( type );
|
2010-07-02 20:49:30 +02:00
|
|
|
}
|
|
|
|
|
|
|
|
void *get_ptr()
|
|
|
|
{ return reinterpret_cast<void*> (&kdata_[0]); }
|
|
|
|
|
|
|
|
~kernel_real()
|
|
|
|
{ std::vector<real_t>().swap( kdata_ ); }
|
|
|
|
|
|
|
|
};
|
|
|
|
|
|
|
|
|
|
|
|
template< typename real_t >
|
2010-08-22 04:43:21 +02:00
|
|
|
void kernel_real<real_t>::compute_kernel( tf_type type )
|
2010-07-02 20:49:30 +02:00
|
|
|
{
|
|
|
|
double
|
|
|
|
kny = cparam_.nx*M_PI/cparam_.lx,
|
|
|
|
fac = cparam_.lx*cparam_.ly*cparam_.lz/pow(2.0*M_PI,3)/(cparam_.nx*cparam_.ny*cparam_.nz),
|
|
|
|
dx = cparam_.lx/cparam_.nx,
|
|
|
|
dy = cparam_.ly/cparam_.ny,
|
|
|
|
dz = cparam_.lz/cparam_.nz,
|
|
|
|
boxlength = cparam_.pcf->getValue<double>("setup","boxlength"),
|
2010-09-01 06:59:31 +02:00
|
|
|
boxlength2 = 0.5*boxlength,
|
2010-07-21 10:10:29 +02:00
|
|
|
nspec = cparam_.pcf->getValue<double>("cosmology","nspec"),
|
|
|
|
pnorm = cparam_.pcf->getValue<double>("cosmology","pnorm"),
|
|
|
|
dplus = cparam_.pcf->getValue<double>("cosmology","dplus");
|
2010-07-02 20:49:30 +02:00
|
|
|
|
|
|
|
bool
|
2010-09-08 10:06:18 +02:00
|
|
|
bavgfine = cparam_.pcf->getValueSafe<bool>("setup","avg_fine",true),
|
2010-09-01 06:59:31 +02:00
|
|
|
bperiodic = cparam_.pcf->getValueSafe<bool>("setup","periodic_TF",true),
|
|
|
|
kspacepoisson = (cparam_.pcf->getValueSafe<bool>("poisson","fft_fine",true)|
|
|
|
|
cparam_.pcf->getValueSafe<bool>("poisson","kspace",false));
|
2010-08-22 04:43:21 +02:00
|
|
|
unsigned
|
|
|
|
levelmax = cparam_.pcf->getValue<unsigned>("setup","levelmax");
|
|
|
|
|
|
|
|
|
|
|
|
if( bavgfine )
|
|
|
|
cparam_.coarse_fact++;
|
2010-09-01 06:59:31 +02:00
|
|
|
|
2010-08-15 01:52:25 +02:00
|
|
|
const int ref_fac = (int)(pow(2,cparam_.coarse_fact)+0.5), 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-08-22 04:43:21 +02:00
|
|
|
new TransferFunction_real(type, cparam_.ptf,nspec,pnorm,dplus,
|
2010-07-21 10:10:29 +02:00
|
|
|
0.25*cparam_.lx/(double)cparam_.nx,2.0*boxlength,kny, (int)pow(2,levelmax+2));
|
|
|
|
|
|
|
|
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-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)
|
|
|
|
for( int i=0; i<=cparam_.nx/2; ++i )
|
|
|
|
for( int j=0; j<=cparam_.ny/2; ++j )
|
|
|
|
for( int k=0; k<=cparam_.nz/2; ++k )
|
2010-07-02 20:49:30 +02:00
|
|
|
{
|
|
|
|
int iix(i), iiy(j), iiz(k);
|
|
|
|
double rr[3], rr2;
|
|
|
|
|
|
|
|
if( iix > (int)cparam_.nx/2 ) iix -= cparam_.nx;
|
|
|
|
if( iiy > (int)cparam_.ny/2 ) iiy -= cparam_.ny;
|
|
|
|
if( iiz > (int)cparam_.nz/2 ) iiz -= cparam_.nz;
|
|
|
|
|
2010-08-15 01:52:25 +02:00
|
|
|
|
|
|
|
//... speed up 8x by copying data to other octants
|
|
|
|
int idx[8];
|
|
|
|
int nx=cparam_.nx, ny = cparam_.ny, nz = cparam_.nz;
|
|
|
|
|
|
|
|
idx[0] = ((i)*cparam_.ny + (j)) * 2*(cparam_.nz/2+1) + (k);
|
|
|
|
idx[1] = ((nx-i)*cparam_.ny + (j)) * 2*(cparam_.nz/2+1) + (k);
|
|
|
|
idx[2] = ((i)*cparam_.ny + (ny-j)) * 2*(cparam_.nz/2+1) + (k);
|
|
|
|
idx[3] = ((nx-i)*cparam_.ny + (ny-j)) * 2*(cparam_.nz/2+1) + (k);
|
|
|
|
idx[4] = ((i)*cparam_.ny + (j)) * 2*(cparam_.nz/2+1) + (nz-k);
|
|
|
|
idx[5] = ((nx-i)*cparam_.ny + (j)) * 2*(cparam_.nz/2+1) + (nz-k);
|
|
|
|
idx[6] = ((i)*cparam_.ny + (ny-j)) * 2*(cparam_.nz/2+1) + (nz-k);
|
|
|
|
idx[7] = ((nx-i)*cparam_.ny + (ny-j)) * 2*(cparam_.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;
|
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;
|
|
|
|
rr[1] = ((double)iiy ) * dy + jj*boxlength;
|
|
|
|
rr[2] = ((double)iiz ) * dz + kk*boxlength;
|
|
|
|
|
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
|
|
|
|
for( int i=0; i<cparam_.nx; ++i )
|
|
|
|
for( int j=0; j<cparam_.ny; ++j )
|
|
|
|
for( int k=0; k<cparam_.nz; ++k )
|
|
|
|
{
|
|
|
|
int iix(i), iiy(j), iiz(k);
|
|
|
|
double rr[3], rr2;
|
|
|
|
|
|
|
|
if( iix > (int)cparam_.nx/2 ) iix -= cparam_.nx;
|
|
|
|
if( iiy > (int)cparam_.ny/2 ) iiy -= cparam_.ny;
|
|
|
|
if( iiz > (int)cparam_.nz/2 ) iiz -= cparam_.nz;
|
|
|
|
|
|
|
|
unsigned idx = (i*cparam_.ny + j) * 2*(cparam_.nz/2+1) + k;
|
|
|
|
|
|
|
|
rr[0] = ((double)iix ) * dx;
|
|
|
|
rr[1] = ((double)iiy ) * dy;
|
|
|
|
rr[2] = ((double)iiz ) * dz;
|
|
|
|
|
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
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
}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-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
|
|
|
|
|
|
|
if( cparam_.is_finest )
|
2010-09-01 06:59:31 +02:00
|
|
|
{
|
2010-08-22 04:43:21 +02:00
|
|
|
kdata_[0] = tfr->compute_real(0.0)*fac;
|
2010-09-01 06:59:31 +02:00
|
|
|
//std::cerr << "T(0) = " << kdata_[0] << std::endl;
|
|
|
|
//kdata_[0] = sqrt(8.0)*tfr->compute_real(0.5*sqrt(3.0)*dx/4)*fac;
|
|
|
|
//std::cerr << "T(0) = " << kdata_[0] << std::endl;
|
|
|
|
}
|
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-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-07-08 07:15:33 +02:00
|
|
|
for( int i=0; i<cparam_.nx; ++i )
|
2010-07-02 20:49:30 +02:00
|
|
|
for( int j=0; j<cparam_.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;
|
|
|
|
|
|
|
|
if( kx > cparam_.nx/2 ) kx -= cparam_.nx;
|
|
|
|
if( ky > cparam_.ny/2 ) ky -= cparam_.ny;
|
|
|
|
|
|
|
|
|
|
|
|
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
|
|
|
|
|
|
|
unsigned q = (i*cparam_.ny+j)*nzp+k;
|
|
|
|
|
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-07-21 10:10:29 +02:00
|
|
|
if( k==0 || k==cparam_.nz/2 )
|
|
|
|
{
|
|
|
|
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-07-28 00:21:50 +02:00
|
|
|
if( cparam_.is_finest )
|
|
|
|
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-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)
|
|
|
|
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;
|
|
|
|
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-01 06:59:31 +02:00
|
|
|
|
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-08-22 04:43:21 +02:00
|
|
|
kernel_k( const parameters& cp, tf_type type )
|
|
|
|
: kernel( cp, type )
|
2010-07-02 20:49:30 +02:00
|
|
|
{
|
|
|
|
kdata_.assign( cparam_.nx*cparam_.ny*2*(cparam_.nz/2+1), 0.0 );
|
2010-08-22 04:43:21 +02:00
|
|
|
compute_kernel(type);
|
2010-07-02 20:49:30 +02:00
|
|
|
}
|
|
|
|
|
|
|
|
void *get_ptr()
|
|
|
|
{ return reinterpret_cast<void*> (&kdata_[0]); }
|
|
|
|
|
|
|
|
~kernel_k()
|
|
|
|
{ std::vector<real_t>().swap( kdata_ ); }
|
|
|
|
|
|
|
|
};
|
|
|
|
|
|
|
|
template< typename real_t >
|
2010-08-22 04:43:21 +02:00
|
|
|
void kernel_k<real_t>::compute_kernel( tf_type type )
|
2010-07-02 20:49:30 +02:00
|
|
|
{
|
|
|
|
double
|
2010-07-21 10:10:29 +02:00
|
|
|
fac = cparam_.lx*cparam_.ly*cparam_.lz/pow(2.0*M_PI,3),
|
2010-07-02 20:49:30 +02:00
|
|
|
boxlength = cparam_.pcf->getValue<double>("setup","boxlength"),
|
|
|
|
nspec = cparam_.pcf->getValue<double>("cosmology","nspec"),
|
|
|
|
pnorm = cparam_.pcf->getValue<double>("cosmology","pnorm"),
|
|
|
|
dplus = cparam_.pcf->getValue<double>("cosmology","dplus");
|
|
|
|
|
2010-08-22 04:43:21 +02:00
|
|
|
TransferFunction_k *tfk = new TransferFunction_k(type,cparam_.ptf,nspec,pnorm,dplus);
|
2010-07-02 20:49:30 +02:00
|
|
|
|
|
|
|
fftw_complex *kdata = reinterpret_cast<fftw_complex*> ( this->get_ptr() );
|
|
|
|
|
|
|
|
unsigned nx = cparam_.nx, ny = cparam_.ny, nz = cparam_.nz, 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-07-02 20:49:30 +02:00
|
|
|
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 )
|
|
|
|
{
|
|
|
|
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-07-28 00:21:50 +02:00
|
|
|
if( k==0 || k==cparam_.nz/2 )
|
|
|
|
{
|
|
|
|
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-08-11 05:35:29 +02:00
|
|
|
}
|
|
|
|
|
2010-08-22 04:43:21 +02:00
|
|
|
|
|
|
|
/**************************************************************************************/
|
|
|
|
/**************************************************************************************/
|
|
|
|
|
|
|
|
|
2010-07-02 20:49:30 +02:00
|
|
|
namespace{
|
|
|
|
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_k<double> > creator_kd("tf_kernel_k_double");
|
|
|
|
convolution::kernel_creator_concrete< convolution::kernel_k<float> > creator_kf("tf_kernel_k_float");
|
|
|
|
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
|
|
|