mirror of
https://github.com/cosmo-sims/monofonIC.git
synced 2024-09-19 17:03:45 +02:00
Merge branch 'random' into 'master'
Merge in new RNG structures and other changes See merge request ohahn/fastlpt!1
This commit is contained in:
Oliver Hahn
commit
9e04da5c92
9 changed files with 433 additions and 131 deletions
|
|
@ -10,6 +10,10 @@ fbase_analysis = output
|
|||
format = gadget2
|
||||
filename = ics_gadget.dat
|
||||
|
||||
[random]
|
||||
generator = NGENIC
|
||||
seed = 9001
|
||||
|
||||
[cosmology]
|
||||
transfer = CLASS #eisenstein
|
||||
Omega_m = 0.302
|
||||
|
|
|
|||
|
|
@ -5,16 +5,215 @@
|
|||
#include <general.hh>
|
||||
#include <grid_fft.hh>
|
||||
|
||||
//! convolution class, respecting Orszag's 3/2 rule
|
||||
template< typename data_t >
|
||||
class OrszagConvolver
|
||||
template< typename data_t, typename derived_t >
|
||||
class BaseConvolver
|
||||
{
|
||||
protected:
|
||||
std::array<size_t,3> np_;
|
||||
std::array<real_t,3> length_;
|
||||
public:
|
||||
BaseConvolver( const std::array<size_t, 3> &N, const std::array<real_t, 3> &L )
|
||||
: np_( N ), length_( L ) {}
|
||||
|
||||
virtual ~BaseConvolver() {}
|
||||
|
||||
// implements convolution of two Fourier-space fields
|
||||
template< typename kfunc1, typename kfunc2, typename opp >
|
||||
void convolve2( kfunc1 kf1, kfunc2 kf2, Grid_FFT<data_t> & res, opp op ) {}
|
||||
|
||||
// implements convolution of three Fourier-space fields
|
||||
template< typename kfunc1, typename kfunc2, typename kfunc3, typename opp >
|
||||
void convolve3( kfunc1 kf1, kfunc2 kf2, kfunc3 kf3, Grid_FFT<data_t> & res, opp op ) {}
|
||||
|
||||
public:
|
||||
|
||||
template< typename opp >
|
||||
void convolve_Hessians( Grid_FFT<data_t> & inl, const std::array<int,2>& d2l, Grid_FFT<data_t> & inr, const std::array<int,2>& d2r, Grid_FFT<data_t> & res, opp op ){
|
||||
// transform to FS in case fields are not
|
||||
inl.FourierTransformForward();
|
||||
inr.FourierTransformForward();
|
||||
// perform convolution of Hessians
|
||||
static_cast<derived_t&>(*this).convolve2(
|
||||
[&]( size_t i, size_t j, size_t k ) -> ccomplex_t{
|
||||
auto kk = inl.template get_k<real_t>(i,j,k);
|
||||
return -kk[d2l[0]] * kk[d2l[1]] * inl.kelem(i,j,k);
|
||||
},
|
||||
[&]( size_t i, size_t j, size_t k ){
|
||||
auto kk = inr.template get_k<real_t>(i,j,k);
|
||||
return -kk[d2r[0]] * kk[d2r[1]] * inr.kelem(i,j,k);
|
||||
}, res, op );
|
||||
}
|
||||
|
||||
template< typename opp >
|
||||
void convolve_Hessians( Grid_FFT<data_t> & inl, const std::array<int,2>& d2l,
|
||||
Grid_FFT<data_t> & inm, const std::array<int,2>& d2m,
|
||||
Grid_FFT<data_t> & inr, const std::array<int,2>& d2r,
|
||||
Grid_FFT<data_t> & res, opp op )
|
||||
{
|
||||
// transform to FS in case fields are not
|
||||
inl.FourierTransformForward();
|
||||
inm.FourierTransformForward();
|
||||
inr.FourierTransformForward();
|
||||
// perform convolution of Hessians
|
||||
static_cast<derived_t&>(*this).convolve3(
|
||||
[&inl,&d2l]( size_t i, size_t j, size_t k ) -> ccomplex_t{
|
||||
auto kk = inl.template get_k<real_t>(i,j,k);
|
||||
return -kk[d2l[0]] * kk[d2l[1]] * inl.kelem(i,j,k);
|
||||
},
|
||||
[&inm,&d2m]( size_t i, size_t j, size_t k ) -> ccomplex_t{
|
||||
auto kk = inm.template get_k<real_t>(i,j,k);
|
||||
return -kk[d2m[0]] * kk[d2m[1]] * inm.kelem(i,j,k);
|
||||
},
|
||||
[&inr,&d2r]( size_t i, size_t j, size_t k ) -> ccomplex_t{
|
||||
auto kk = inr.template get_k<real_t>(i,j,k);
|
||||
return -kk[d2r[0]] * kk[d2r[1]] * inr.kelem(i,j,k);
|
||||
}, res, op );
|
||||
}
|
||||
|
||||
template< typename opp >
|
||||
void convolve_SumHessians( Grid_FFT<data_t> & inl, const std::array<int,2>& d2l, Grid_FFT<data_t> & inr, const std::array<int,2>& d2r1,
|
||||
const std::array<int,2>& d2r2, Grid_FFT<data_t> & res, opp op ){
|
||||
// transform to FS in case fields are not
|
||||
inl.FourierTransformForward();
|
||||
inr.FourierTransformForward();
|
||||
// perform convolution of Hessians
|
||||
static_cast<derived_t&>(*this).convolve2(
|
||||
[&inl,&d2l]( size_t i, size_t j, size_t k ) -> ccomplex_t{
|
||||
auto kk = inl.template get_k<real_t>(i,j,k);
|
||||
return -kk[d2l[0]] * kk[d2l[1]] * inl.kelem(i,j,k);
|
||||
},
|
||||
[&inr,&d2r1,&d2r2]( size_t i, size_t j, size_t k ) -> ccomplex_t{
|
||||
auto kk = inr.template get_k<real_t>(i,j,k);
|
||||
return (-kk[d2r1[0]] * kk[d2r1[1]] -kk[d2r2[0]] * kk[d2r2[1]]) * inr.kelem(i,j,k);
|
||||
}, res, op );
|
||||
}
|
||||
};
|
||||
|
||||
//! naive convolution class, disrespecting aliasing
|
||||
template< typename data_t >
|
||||
class NaiveConvolver : public BaseConvolver<data_t, NaiveConvolver<data_t>>
|
||||
{
|
||||
protected:
|
||||
Grid_FFT<data_t> *fbuf1_, *fbuf2_;
|
||||
|
||||
using BaseConvolver<data_t, NaiveConvolver<data_t>>::np_;
|
||||
using BaseConvolver<data_t, NaiveConvolver<data_t>>::length_;
|
||||
|
||||
public:
|
||||
NaiveConvolver( const std::array<size_t, 3> &N, const std::array<real_t, 3> &L )
|
||||
: BaseConvolver<data_t, NaiveConvolver<data_t>>( N, L )
|
||||
{
|
||||
fbuf1_ = new Grid_FFT<data_t>(N, length_, kspace_id);
|
||||
fbuf2_ = new Grid_FFT<data_t>(N, length_, kspace_id);
|
||||
}
|
||||
|
||||
~NaiveConvolver()
|
||||
{
|
||||
delete fbuf1_;
|
||||
delete fbuf2_;
|
||||
}
|
||||
|
||||
template< typename kfunc1, typename kfunc2, typename opp >
|
||||
void convolve2( kfunc1 kf1, kfunc2 kf2, Grid_FFT<data_t> & res, opp op )
|
||||
{
|
||||
//... prepare data 1
|
||||
fbuf1_->FourierTransformForward(false);
|
||||
this->copy_in( kf1, *fbuf1_ );
|
||||
|
||||
//... prepare data 2
|
||||
fbuf2_->FourierTransformForward(false);
|
||||
this->copy_in( kf2, *fbuf2_ );
|
||||
|
||||
//... convolve
|
||||
fbuf1_->FourierTransformBackward();
|
||||
fbuf2_->FourierTransformBackward();
|
||||
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < fbuf1_->ntot_; ++i){
|
||||
(*fbuf2_).relem(i) *= (*fbuf1_).relem(i);
|
||||
}
|
||||
fbuf2_->FourierTransformForward();
|
||||
//... copy data back
|
||||
res.FourierTransformForward();
|
||||
copy_out( *fbuf2_, op, res );
|
||||
}
|
||||
|
||||
template< typename kfunc1, typename kfunc2, typename kfunc3, typename opp >
|
||||
void convolve3( kfunc1 kf1, kfunc2 kf2, kfunc3 kf3, Grid_FFT<data_t> & res, opp op )
|
||||
{
|
||||
//... prepare data 1
|
||||
fbuf1_->FourierTransformForward(false);
|
||||
this->copy_in( kf1, *fbuf1_ );
|
||||
|
||||
//... prepare data 2
|
||||
fbuf2_->FourierTransformForward(false);
|
||||
this->copy_in( kf2, *fbuf2_ );
|
||||
|
||||
//... convolve
|
||||
fbuf1_->FourierTransformBackward();
|
||||
fbuf2_->FourierTransformBackward();
|
||||
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < fbuf1_->ntot_; ++i){
|
||||
(*fbuf2_).relem(i) *= (*fbuf1_).relem(i);
|
||||
}
|
||||
|
||||
//... prepare data 2
|
||||
fbuf1_->FourierTransformForward(false);
|
||||
this->copy_in( kf3, *fbuf1_ );
|
||||
|
||||
//... convolve
|
||||
fbuf1_->FourierTransformBackward();
|
||||
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < fbuf1_->ntot_; ++i){
|
||||
(*fbuf2_).relem(i) *= (*fbuf1_).relem(i);
|
||||
}
|
||||
|
||||
fbuf2_->FourierTransformForward();
|
||||
//... copy data back
|
||||
res.FourierTransformForward();
|
||||
copy_out( *fbuf2_, op, res );
|
||||
}
|
||||
|
||||
private:
|
||||
template< typename kfunc >
|
||||
void copy_in( kfunc kf, Grid_FFT<data_t>& g )
|
||||
{
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < g.size(0); ++i){
|
||||
for (size_t j = 0; j < g.size(1); ++j){
|
||||
for (size_t k = 0; k < g.size(2); ++k){
|
||||
g.kelem(i, j, k) = kf(i, j, k);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template< typename opp >
|
||||
void copy_out( Grid_FFT<data_t>& f, opp op, Grid_FFT<data_t>& res )
|
||||
{
|
||||
#pragma omp parallel for
|
||||
for (size_t i = 0; i < f.size(0); ++i){
|
||||
for (size_t j = 0; j < f.size(1); ++j){
|
||||
for (size_t k = 0; k < f.size(2); ++k){
|
||||
res.kelem(i, j, k) = op(f.kelem(i, j, k), res.kelem(i, j, k));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
//! convolution class, respecting Orszag's 3/2 rule
|
||||
template< typename data_t >
|
||||
class OrszagConvolver : public BaseConvolver<data_t, OrszagConvolver<data_t>>
|
||||
{
|
||||
private:
|
||||
Grid_FFT<data_t> *f1p_, *f2p_;
|
||||
Grid_FFT<data_t> *fbuf_, *fbuf2_;
|
||||
|
||||
std::array<size_t,3> np_;
|
||||
std::array<real_t,3> length_;
|
||||
using BaseConvolver<data_t, OrszagConvolver<data_t>>::np_;
|
||||
using BaseConvolver<data_t, OrszagConvolver<data_t>>::length_;
|
||||
|
||||
ccomplex_t *crecvbuf_;
|
||||
real_t *recvbuf_;
|
||||
|
|
@ -22,8 +221,6 @@ protected:
|
|||
std::vector<ptrdiff_t> offsets_, offsetsp_;
|
||||
std::vector<size_t> sizes_, sizesp_;
|
||||
|
||||
private:
|
||||
|
||||
int get_task(ptrdiff_t index, const std::vector<ptrdiff_t>& offsets, const std::vector<size_t>& sizes, const int ntasks )
|
||||
{
|
||||
int itask = 0;
|
||||
|
|
@ -35,7 +232,8 @@ public:
|
|||
|
||||
|
||||
OrszagConvolver( const std::array<size_t, 3> &N, const std::array<real_t, 3> &L )
|
||||
: np_({3*N[0]/2,3*N[1]/2,3*N[2]/2}), length_(L)
|
||||
: BaseConvolver<data_t,OrszagConvolver<data_t>>( {3*N[0]/2,3*N[1]/2,3*N[2]/2}, L )
|
||||
//: np_({3*N[0]/2,3*N[1]/2,3*N[2]/2}), length_(L)
|
||||
{
|
||||
//... create temporaries
|
||||
f1p_ = new Grid_FFT<data_t>(np_, length_, kspace_id);
|
||||
|
|
@ -80,67 +278,6 @@ public:
|
|||
#endif
|
||||
}
|
||||
|
||||
template< typename opp >
|
||||
void convolve_Hessians( Grid_FFT<data_t> & inl, const std::array<int,2>& d2l, Grid_FFT<data_t> & inr, const std::array<int,2>& d2r, Grid_FFT<data_t> & res, opp op ){
|
||||
// transform to FS in case fields are not
|
||||
inl.FourierTransformForward();
|
||||
inr.FourierTransformForward();
|
||||
// perform convolution of Hessians
|
||||
this->convolve2(
|
||||
[&]( size_t i, size_t j, size_t k ) -> ccomplex_t{
|
||||
auto kk = inl.template get_k<real_t>(i,j,k);
|
||||
return -kk[d2l[0]] * kk[d2l[1]] * inl.kelem(i,j,k);
|
||||
},
|
||||
[&]( size_t i, size_t j, size_t k ){
|
||||
auto kk = inr.template get_k<real_t>(i,j,k);
|
||||
return -kk[d2r[0]] * kk[d2r[1]] * inr.kelem(i,j,k);
|
||||
}, res, op );
|
||||
}
|
||||
|
||||
template< typename opp >
|
||||
void convolve_Hessians( Grid_FFT<data_t> & inl, const std::array<int,2>& d2l,
|
||||
Grid_FFT<data_t> & inm, const std::array<int,2>& d2m,
|
||||
Grid_FFT<data_t> & inr, const std::array<int,2>& d2r,
|
||||
Grid_FFT<data_t> & res, opp op )
|
||||
{
|
||||
// transform to FS in case fields are not
|
||||
inl.FourierTransformForward();
|
||||
inm.FourierTransformForward();
|
||||
inr.FourierTransformForward();
|
||||
// perform convolution of Hessians
|
||||
this->convolve3(
|
||||
[&inl,&d2l]( size_t i, size_t j, size_t k ) -> ccomplex_t{
|
||||
auto kk = inl.template get_k<real_t>(i,j,k);
|
||||
return -kk[d2l[0]] * kk[d2l[1]] * inl.kelem(i,j,k);
|
||||
},
|
||||
[&inm,&d2m]( size_t i, size_t j, size_t k ) -> ccomplex_t{
|
||||
auto kk = inm.template get_k<real_t>(i,j,k);
|
||||
return -kk[d2m[0]] * kk[d2m[1]] * inm.kelem(i,j,k);
|
||||
},
|
||||
[&inr,&d2r]( size_t i, size_t j, size_t k ) -> ccomplex_t{
|
||||
auto kk = inr.template get_k<real_t>(i,j,k);
|
||||
return -kk[d2r[0]] * kk[d2r[1]] * inr.kelem(i,j,k);
|
||||
}, res, op );
|
||||
}
|
||||
|
||||
template< typename opp >
|
||||
void convolve_SumHessians( Grid_FFT<data_t> & inl, const std::array<int,2>& d2l, Grid_FFT<data_t> & inr, const std::array<int,2>& d2r1,
|
||||
const std::array<int,2>& d2r2, Grid_FFT<data_t> & res, opp op ){
|
||||
// transform to FS in case fields are not
|
||||
inl.FourierTransformForward();
|
||||
inr.FourierTransformForward();
|
||||
// perform convolution of Hessians
|
||||
this->convolve2(
|
||||
[&inl,&d2l]( size_t i, size_t j, size_t k ) -> ccomplex_t{
|
||||
auto kk = inl.template get_k<real_t>(i,j,k);
|
||||
return -kk[d2l[0]] * kk[d2l[1]] * inl.kelem(i,j,k);
|
||||
},
|
||||
[&inr,&d2r1,&d2r2]( size_t i, size_t j, size_t k ) -> ccomplex_t{
|
||||
auto kk = inr.template get_k<real_t>(i,j,k);
|
||||
return (-kk[d2r1[0]] * kk[d2r1[1]] -kk[d2r2[0]] * kk[d2r2[1]]) * inr.kelem(i,j,k);
|
||||
}, res, op );
|
||||
}
|
||||
|
||||
template< typename kfunc1, typename kfunc2, typename opp >
|
||||
void convolve2( kfunc1 kf1, kfunc2 kf2, Grid_FFT<data_t> & res, opp op )
|
||||
{
|
||||
|
|
|
|||
|
|
@ -90,6 +90,11 @@ public:
|
|||
|
||||
size_t size(size_t i) const { return sizes_[i]; }
|
||||
|
||||
const bounding_box<size_t>& get_global_range( void ) const
|
||||
{
|
||||
return global_range_;
|
||||
}
|
||||
|
||||
void zero()
|
||||
{
|
||||
#pragma omp parallel for
|
||||
|
|
@ -148,8 +153,7 @@ public:
|
|||
}
|
||||
|
||||
void cell_pos( int ilevel, size_t i, size_t j, size_t k, double* x ) const {
|
||||
#warning needs to be fixed for MPI
|
||||
x[0] = double(i)/size(0);
|
||||
x[0] = double(i+local_0_start_)/size(0);
|
||||
x[1] = double(j)/size(1);
|
||||
x[2] = double(k)/size(2);
|
||||
}
|
||||
|
|
|
|||
|
|
@ -2,6 +2,8 @@
|
|||
|
||||
#include <map>
|
||||
#include <config_file.hh>
|
||||
#include <general.hh>
|
||||
#include <grid_fft.hh>
|
||||
|
||||
#define DEF_RAN_CUBE_SIZE 32
|
||||
|
||||
|
|
@ -16,6 +18,7 @@ class RNG_plugin
|
|||
}
|
||||
virtual ~RNG_plugin() {}
|
||||
virtual bool isMultiscale() const = 0;
|
||||
virtual void Fill_Grid( Grid_FFT<real_t>& g ) const = 0;
|
||||
//virtual void FillGrid(int level, DensityGrid<real_t> &R) = 0;
|
||||
};
|
||||
|
||||
|
|
|
|||
|
|
@ -87,6 +87,8 @@ void Grid_FFT<data_t>::Setup(void)
|
|||
|
||||
local_0_size_ = n_[0];
|
||||
local_1_size_ = n_[1];
|
||||
local_0_start_= 0;
|
||||
local_1_start_= 0;
|
||||
|
||||
if (space_ == rspace_id)
|
||||
{
|
||||
|
|
@ -562,9 +564,7 @@ void Grid_FFT<data_t>::Write_PowerSpectrum( std::string ofname )
|
|||
{
|
||||
#endif
|
||||
std::ofstream ofs( ofname.c_str());
|
||||
std::size_t numcells = size(0) * size(1) * size(2);
|
||||
|
||||
//ofs << "# a = " << aexp << std::endl;
|
||||
|
||||
ofs << "# " << std::setw(14) << "k" << std::setw(16) << "P(k)" << std::setw(16) << "err. P(k)"
|
||||
<< std::setw(16) <<"#modes" << "\n";
|
||||
|
||||
|
|
|
|||
81
src/main.cc
81
src/main.cc
|
|
@ -71,7 +71,7 @@ int main( int argc, char** argv )
|
|||
{
|
||||
// print_region_generator_plugins();
|
||||
print_TransferFunction_plugins();
|
||||
// print_RNG_plugins();
|
||||
print_RNG_plugins();
|
||||
print_output_plugins();
|
||||
|
||||
csoca::elog << "In order to run, you need to specify a parameter file!" << std::endl;
|
||||
|
|
@ -146,8 +146,19 @@ int main( int argc, char** argv )
|
|||
Grid_FFT<real_t> phi3b({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
|
||||
OrszagConvolver<real_t> Conv({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
// NaiveConvolver<real_t> Conv({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
|
||||
//--------------------------------------------------------------------
|
||||
// Some operators to add or subtract terms
|
||||
auto assign_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return res; };
|
||||
auto add_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return val+res; };
|
||||
auto add2_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return val+2.0*res; };
|
||||
auto sub_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return val-res; };
|
||||
auto sub2_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return val-2.0*res; };
|
||||
//--------------------------------------------------------------------
|
||||
|
||||
phi.FillRandomReal(6519);
|
||||
//phi.FillRandomReal(6519);
|
||||
the_random_number_generator->Fill_Grid( phi );
|
||||
|
||||
//======================================================================
|
||||
//... compute 1LPT displacement potential ....
|
||||
|
|
@ -163,19 +174,15 @@ int main( int argc, char** argv )
|
|||
ccomplex_t delta = x * the_cosmo_calc->GetAmplitude(kmod, total);
|
||||
return -delta / (kmod * kmod) * phifac / volfac;
|
||||
});
|
||||
|
||||
phi.zero_DC_mode();
|
||||
csoca::ilog << " took " << get_wtime()-wtime << "s" << std::endl;
|
||||
|
||||
//////////////////////////////////////////////////////////////////////////////////////////////////
|
||||
auto assign_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return res; };
|
||||
auto add_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return val+res; };
|
||||
auto add2_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return val+2.0*res; };
|
||||
auto sub_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return val-res; };
|
||||
auto sub2_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return val-2.0*res; };
|
||||
|
||||
//======================================================================
|
||||
//... compute 2LPT displacement potential ....
|
||||
|
||||
wtime = get_wtime();
|
||||
csoca::ilog << "Computing phi(2) term..." << std::flush;
|
||||
// Compute the source term for phi(2)
|
||||
Conv.convolve_SumHessians( phi, {0,0}, phi, {1,1}, {2,2}, phi2, assign_op );
|
||||
Conv.convolve_Hessians( phi, {1,1}, phi, {2,2}, phi2, add_op );
|
||||
Conv.convolve_Hessians( phi, {0,1}, phi, {0,1}, phi2, sub_op );
|
||||
|
|
@ -184,58 +191,36 @@ int main( int argc, char** argv )
|
|||
phi2.apply_InverseLaplacian();
|
||||
phi2 /= phifac;
|
||||
csoca::ilog << " took " << get_wtime()-wtime << "s" << std::endl;
|
||||
|
||||
// phi2.FourierTransformForward();
|
||||
// phi2.apply_function_k_dep([&](auto x, auto k) {
|
||||
// real_t kmod2 = k.norm_squared();
|
||||
// return x * (-1.0 / kmod2) * phifac / phifac / phifac;
|
||||
// });
|
||||
// phi2.zero_DC_mode();
|
||||
|
||||
|
||||
//======================================================================
|
||||
//... compute 3LPT displacement potential
|
||||
|
||||
wtime = get_wtime();
|
||||
csoca::ilog << "Computing phi(3a) term..." << std::flush;
|
||||
Conv.convolve_SumHessians( phi, {0,0}, phi2, {1,1}, {2,2}, phi3a, assign_op );
|
||||
Conv.convolve_SumHessians( phi, {1,1}, phi2, {2,2}, {0,0}, phi3a, add_op );
|
||||
Conv.convolve_SumHessians( phi, {2,2}, phi2, {0,0}, {1,1}, phi3a, add_op );
|
||||
Conv.convolve_Hessians( phi, {0,1}, phi2, {0,1}, phi3a, sub2_op );
|
||||
Conv.convolve_Hessians( phi, {0,2}, phi2, {0,2}, phi3a, sub2_op );
|
||||
Conv.convolve_Hessians( phi, {1,2}, phi2, {1,2}, phi3a, sub2_op );
|
||||
Conv.convolve_Hessians( phi, {0,0}, phi, {1,1}, phi, {2,2}, phi3a, assign_op );
|
||||
Conv.convolve_Hessians( phi, {0,1}, phi, {0,2}, phi, {1,2}, phi3a, add2_op );
|
||||
Conv.convolve_Hessians( phi, {1,2}, phi, {1,2}, phi, {0,0}, phi3a, sub_op );
|
||||
Conv.convolve_Hessians( phi, {0,2}, phi, {0,2}, phi, {1,1}, phi3a, sub_op );
|
||||
Conv.convolve_Hessians( phi, {0,1}, phi, {0,1}, phi, {2,2}, phi3a, sub_op );
|
||||
phi3a.apply_InverseLaplacian();
|
||||
phi3a *= 0.5/phifac; // factor 1/2 from definition of phi(3a)!
|
||||
// phi3a.apply_function_k_dep([&](auto x, auto k) {
|
||||
// return 0.5 * x;
|
||||
// });
|
||||
phi3a /= phifac*phifac;
|
||||
csoca::ilog << " took " << get_wtime()-wtime << "s" << std::endl;
|
||||
|
||||
|
||||
// phi3a.FourierTransformForward();
|
||||
// phi3a.apply_function_k_dep([&](auto x, auto k) {
|
||||
// real_t kmod2 = k.norm_squared();
|
||||
// return x * (-1.0 / kmod2) * phifac / phifac / phifac;
|
||||
// });
|
||||
// phi3a.zero_DC_mode();
|
||||
//...
|
||||
|
||||
wtime = get_wtime();
|
||||
csoca::ilog << "Computing phi(3b) term..." << std::flush;
|
||||
Conv.convolve_Hessians( phi, {0,0}, phi, {1,1}, phi, {2,2}, phi3b, assign_op );
|
||||
Conv.convolve_Hessians( phi, {0,1}, phi, {0,2}, phi, {1,2}, phi3b, add2_op );
|
||||
Conv.convolve_Hessians( phi, {1,2}, phi, {1,2}, phi, {0,0}, phi3b, sub_op );
|
||||
Conv.convolve_Hessians( phi, {0,2}, phi, {0,2}, phi, {1,1}, phi3b, sub_op );
|
||||
Conv.convolve_Hessians( phi, {0,1}, phi, {0,1}, phi, {2,2}, phi3b, sub_op );
|
||||
|
||||
// phi3b.FourierTransformForward();
|
||||
// phi3b.apply_function_k_dep([&](auto x, auto k) {
|
||||
// real_t kmod2 = k.norm_squared();
|
||||
// return x * (-1.0 / kmod2) * phifac / phifac / phifac/phifac;
|
||||
// });
|
||||
// phi3b.zero_DC_mode();
|
||||
Conv.convolve_SumHessians( phi, {0,0}, phi2, {1,1}, {2,2}, phi3b, assign_op );
|
||||
Conv.convolve_SumHessians( phi, {1,1}, phi2, {2,2}, {0,0}, phi3b, add_op );
|
||||
Conv.convolve_SumHessians( phi, {2,2}, phi2, {0,0}, {1,1}, phi3b, add_op );
|
||||
Conv.convolve_Hessians( phi, {0,1}, phi2, {0,1}, phi3b, sub2_op );
|
||||
Conv.convolve_Hessians( phi, {0,2}, phi2, {0,2}, phi3b, sub2_op );
|
||||
Conv.convolve_Hessians( phi, {1,2}, phi2, {1,2}, phi3b, sub2_op );
|
||||
phi3b.apply_InverseLaplacian();
|
||||
phi3b /= phifac*phifac;
|
||||
phi3b *= 0.5/phifac; // factor 1/2 from definition of phi(3a)!
|
||||
csoca::ilog << " took " << get_wtime()-wtime << "s" << std::endl;
|
||||
|
||||
|
||||
///////////////////////////////////////////////////////////////////////
|
||||
// we store the densities here if we compute them
|
||||
const bool compute_densities = true;
|
||||
|
|
|
|||
|
|
@ -40,6 +40,8 @@ public:
|
|||
|
||||
bool isMultiscale() const { return true; }
|
||||
|
||||
void Fill_Grid( Grid_FFT<real_t>& g ) const { }
|
||||
|
||||
void initialize_for_grid_structure()//const refinement_hierarchy &refh)
|
||||
{
|
||||
//prefh_ = &refh;
|
||||
|
|
|
|||
162
src/plugins/random_ngenic.cc
Normal file
162
src/plugins/random_ngenic.cc
Normal file
|
|
@ -0,0 +1,162 @@
|
|||
|
||||
#include <general.hh>
|
||||
#include <random_plugin.hh>
|
||||
#include <config_file.hh>
|
||||
|
||||
#include <vector>
|
||||
#include <cmath>
|
||||
|
||||
#include <grid_fft.hh>
|
||||
#include <gsl/gsl_rng.h>
|
||||
|
||||
class RNG_ngenic : public RNG_plugin
|
||||
{
|
||||
private:
|
||||
gsl_rng *pRandomGenerator_;
|
||||
long RandomSeed_;
|
||||
size_t nres_;
|
||||
std::vector<unsigned int> SeedTable_;
|
||||
|
||||
public:
|
||||
explicit RNG_ngenic(ConfigFile &cf) : RNG_plugin(cf)
|
||||
{
|
||||
|
||||
RandomSeed_ = cf.GetValue<long>("random", "seed");
|
||||
nres_ = cf.GetValue<size_t>("setup", "GridRes");
|
||||
pRandomGenerator_ = gsl_rng_alloc(gsl_rng_ranlxd1);
|
||||
|
||||
SeedTable_.assign(nres_ * nres_, 0u);
|
||||
|
||||
for (size_t i = 0; i < nres_ / 2; ++i)
|
||||
{
|
||||
for (size_t j = 0; j < i; j++)
|
||||
SeedTable_[i * nres_ + j] = 0x7fffffff * gsl_rng_uniform(pRandomGenerator_);
|
||||
|
||||
for (size_t j = 0; j < i + 1; j++)
|
||||
SeedTable_[j * nres_ + i] = 0x7fffffff * gsl_rng_uniform(pRandomGenerator_);
|
||||
|
||||
for (size_t j = 0; j < i; j++)
|
||||
SeedTable_[(nres_ - 1 - i) * nres_ + j] = 0x7fffffff * gsl_rng_uniform(pRandomGenerator_);
|
||||
|
||||
for (size_t j = 0; j < i + 1; j++)
|
||||
SeedTable_[(nres_ - 1 - j) * nres_ + i] = 0x7fffffff * gsl_rng_uniform(pRandomGenerator_);
|
||||
|
||||
for (size_t j = 0; j < i; j++)
|
||||
SeedTable_[i * nres_ + (nres_ - 1 - j)] = 0x7fffffff * gsl_rng_uniform(pRandomGenerator_);
|
||||
|
||||
for (size_t j = 0; j < i + 1; j++)
|
||||
SeedTable_[j * nres_ + (nres_ - 1 - i)] = 0x7fffffff * gsl_rng_uniform(pRandomGenerator_);
|
||||
|
||||
for (size_t j = 0; j < i; j++)
|
||||
SeedTable_[(nres_ - 1 - i) * nres_ + (nres_ - 1 - j)] = 0x7fffffff * gsl_rng_uniform(pRandomGenerator_);
|
||||
|
||||
for (size_t j = 0; j < i + 1; j++)
|
||||
SeedTable_[(nres_ - 1 - j) * nres_ + (nres_ - 1 - i)] = 0x7fffffff * gsl_rng_uniform(pRandomGenerator_);
|
||||
}
|
||||
}
|
||||
|
||||
virtual ~RNG_ngenic()
|
||||
{
|
||||
gsl_rng_free(pRandomGenerator_);
|
||||
}
|
||||
|
||||
bool isMultiscale() const { return false; }
|
||||
|
||||
void Fill_Grid(Grid_FFT<real_t> &g) const
|
||||
{
|
||||
double fnorm = std::pow((double)nres_, -1.5);
|
||||
g.FourierTransformForward(false);
|
||||
|
||||
#ifdef USE_MPI
|
||||
// transform is transposed!
|
||||
for (size_t j = 0; j < g.size(0); ++j)
|
||||
{
|
||||
for (size_t i = 0; i < g.size(1); ++i)
|
||||
{
|
||||
#else
|
||||
for (size_t i = 0; i < g.size(0); ++i)
|
||||
{
|
||||
for (size_t j = 0; j < g.size(1); ++j)
|
||||
{
|
||||
#endif
|
||||
ptrdiff_t ii = (i>0)? g.size(1) - i : 0;
|
||||
gsl_rng_set( pRandomGenerator_, SeedTable_[i * nres_ + j]);
|
||||
for (size_t k = 0; k < g.size(2); ++k)
|
||||
{
|
||||
double phase = gsl_rng_uniform(pRandomGenerator_) * 2 * M_PI;
|
||||
double ampl;
|
||||
do
|
||||
{
|
||||
ampl = gsl_rng_uniform(pRandomGenerator_);
|
||||
} while (ampl == 0);
|
||||
if (i == nres_ / 2 || j == nres_ / 2 || k == nres_ / 2)
|
||||
continue;
|
||||
if (i == 0 && j == 0 && k == 0)
|
||||
continue;
|
||||
|
||||
real_t rp = -std::sqrt(-std::log(ampl)) * std::cos(phase);// * fnorm;
|
||||
real_t ip = -std::sqrt(-std::log(ampl)) * std::sin(phase);// * fnorm;
|
||||
ccomplex_t zrand(rp,ip);
|
||||
|
||||
if (k > 0)
|
||||
{
|
||||
g.kelem(i,j,k) = zrand;
|
||||
// RE(knoise[(i * res + j) * (res / 2 + 1) + k]) = rp;
|
||||
// IM(knoise[(i * res + j) * (res / 2 + 1) + k]) = ip;
|
||||
}
|
||||
else /* k=0 plane needs special treatment */
|
||||
{
|
||||
if (i == 0)
|
||||
{
|
||||
if (j >= nres_ / 2)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
else
|
||||
{
|
||||
int jj = (int)nres_ - (int)j; /* note: j!=0 surely holds at this point */
|
||||
g.kelem(i,j,k) = zrand;
|
||||
g.kelem(i,jj,k) = std::conj(zrand);
|
||||
|
||||
// RE(knoise[(i * res + j) * (res / 2 + 1) + k]) = rp;
|
||||
// IM(knoise[(i * res + j) * (res / 2 + 1) + k]) = ip;
|
||||
|
||||
// RE(knoise[(i * res + jj) * (res / 2 + 1) + k]) = rp;
|
||||
// IM(knoise[(i * res + jj) * (res / 2 + 1) + k]) = -ip;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if (i >= nres_ / 2)
|
||||
{
|
||||
continue;
|
||||
}
|
||||
else
|
||||
{
|
||||
ptrdiff_t ii = (i>0)? nres_ - i : 0;
|
||||
ptrdiff_t jj = (j>0)? nres_ - j : 0;
|
||||
|
||||
g.kelem(i,j,k) = zrand;
|
||||
|
||||
// RE(knoise[(i * res + j) * (res / 2 + 1) + k]) = rp;
|
||||
// IM(knoise[(i * res + j) * (res / 2 + 1) + k]) = ip;
|
||||
|
||||
if (ii >= 0 && ii < (int)nres_)
|
||||
{
|
||||
// RE(knoise[(ii * res + jj) * (res / 2 + 1) + k]) = rp;
|
||||
// IM(knoise[(ii * res + jj) * (res / 2 + 1) + k]) = -ip;
|
||||
g.kelem(ii,jj,k) = std::conj(zrand);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
namespace
|
||||
{
|
||||
RNG_plugin_creator_concrete<RNG_ngenic> creator("NGENIC");
|
||||
}
|
||||
|
|
@ -25,10 +25,10 @@ private:
|
|||
std::vector<double> tab_lnk_, tab_dtot_, tab_dc_, tab_db_, tab_ttot_, tab_tc_, tab_tb_;
|
||||
gsl_interp_accel *gsl_ia_dtot_, *gsl_ia_dc_, *gsl_ia_db_, *gsl_ia_ttot_, *gsl_ia_tc_, *gsl_ia_tb_;
|
||||
gsl_spline *gsl_sp_dtot_, *gsl_sp_dc_, *gsl_sp_db_, *gsl_sp_ttot_, *gsl_sp_tc_, *gsl_sp_tb_;
|
||||
double Omega_m_, Omega_b_, zstart_, ztarget_, kmax_, kmin_, h_;
|
||||
double Omega_m_, Omega_b_, N_ur_, zstart_, ztarget_, kmax_, kmin_, h_;
|
||||
|
||||
void ClassEngine_get_data( void ){
|
||||
std::vector<double> d_g, d_ur, t_g, t_ur, phi, psi;
|
||||
std::vector<double> d_ncdm, t_ncdm, phi, psi;
|
||||
|
||||
csoca::ilog << "Computing transfer function via ClassEngine..." << std::flush;
|
||||
double wtime = get_wtime();
|
||||
|
|
@ -41,9 +41,13 @@ private:
|
|||
pars.add("P_k_max_h/Mpc", kmax_);
|
||||
pars.add("h",h_);
|
||||
pars.add("Omega_b",Omega_b_);
|
||||
pars.add("Omega_ur",0.0);
|
||||
// pars.add("Omega_k",0.0);
|
||||
// pars.add("Omega_ur",0.0);
|
||||
pars.add("N_ur",N_ur_);
|
||||
pars.add("Omega_cdm",Omega_m_-Omega_b_);
|
||||
pars.add("Omega_Lambda",1.0-Omega_m_);
|
||||
// pars.add("Omega_fld",0.0);
|
||||
// pars.add("Omega_scf",0.0);
|
||||
pars.add("A_s",2.42e-9);
|
||||
pars.add("n_s",.96);
|
||||
pars.add("output","dTk,vTk");
|
||||
|
|
@ -51,8 +55,8 @@ private:
|
|||
|
||||
std::unique_ptr<ClassEngine> CE = std::make_unique<ClassEngine>(pars, false);
|
||||
|
||||
CE->getTk(target_redshift, tab_lnk_, d_g, tab_db_, tab_dc_, d_ur, tab_dtot_,
|
||||
phi, psi, t_g, tab_tb_, t_ur, tab_ttot_);
|
||||
CE->getTk(target_redshift, tab_lnk_, tab_dc_, tab_db_, d_ncdm, tab_dtot_,
|
||||
tab_tc_, tab_tb_, t_ncdm, tab_ttot_, phi, psi );
|
||||
tab_tc_ = tab_ttot_;
|
||||
#warning need to fix CDM velocities
|
||||
|
||||
|
|
@ -67,8 +71,9 @@ public:
|
|||
h_ = cf.GetValue<double>("cosmology","H0") / 100.0;
|
||||
Omega_m_ = cf.GetValue<double>("cosmology","Omega_m");
|
||||
Omega_b_ = cf.GetValue<double>("cosmology","Omega_b");
|
||||
zstart_ = cf.GetValue<double>("setup","zstart");
|
||||
N_ur_ = cf.GetValueSafe<double>("cosmology","N_ur", 3.046);
|
||||
ztarget_ = cf.GetValueSafe<double>("cosmology","ztarget",0.0);
|
||||
zstart_ = cf.GetValue<double>("setup","zstart");
|
||||
kmax_ = 1000.0;
|
||||
|
||||
this->ClassEngine_get_data();
|
||||
|
|
@ -130,12 +135,12 @@ public:
|
|||
}
|
||||
|
||||
double d = (k<=kmin_)? gsl_spline_eval(splineT, std::log(kmin_), accT)
|
||||
: gsl_spline_eval(splineT, std::log(k), accT);
|
||||
: gsl_spline_eval(splineT, std::log(k*h_), accT);
|
||||
return -d/(k*k);
|
||||
}
|
||||
|
||||
inline double get_kmin(void) const { return std::exp(tab_lnk_[0]); }
|
||||
inline double get_kmax(void) const { return std::exp(tab_lnk_[tab_lnk_.size()-1]); }
|
||||
inline double get_kmin(void) const { return std::exp(tab_lnk_[0])/h_; }
|
||||
inline double get_kmax(void) const { return std::exp(tab_lnk_[tab_lnk_.size()-1])/h_; }
|
||||
};
|
||||
|
||||
namespace {
|
||||
|
|
|
|||
Loading…
Reference in a new issue