mirror/monofonIC
1
0
Fork 0
mirror of https://github.com/cosmo-sims/monofonIC.git synced 2024-09-16 13:33:45 +02:00
Code Releases Activity

Added fNL and gNL part

Browse source
This commit is contained in:
adrigut10 2024-04-15 11:34:09 +02:00
parent 6a677a08f9
commit edb7c08b86
4 changed files with 111 additions and 17 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

13
example.conf
View file

@ -54,7 +54,8 @@ ParameterSet = Planck2018EE+BAO+SN # specify a pre-defined parameter set, or
# m_nu3 = 0.0
# w_0 = -1.0 # not supported yet!
# w_a = 0.0 # not supported yet!
# fnl = 100.0
# gnl = 0.0
ZeroRadiation = false # For Back-scaling only: set to true if your simulation code
# cannot deal with Omega_r!=0 in its background FLRW model
@ -132,8 +133,8 @@ NumThreads = 8
# grafic_use_SPT = no # if no then uses PPT, otherwise linear SPT
##> Gadget-2/3 'fortran unformatted binary'-style format
# format = gadget2
# filename = ics_gadget.dat
#format = gadget2
#filename = ics_gadget.dat
# UseLongids = false
##> Gadget-2/3 HDF5 format
@ -156,6 +157,6 @@ NumThreads = 8
# UseLongids = true
##> Generic HDF5 output format for testing or PT-based calculations
# format = generic
# filename = debug.hdf5
# generic_out_eulerian = yes # if yes then uses PPT for output
format = generic
filename = debug.hdf5
generic_out_eulerian = yes # if yes then uses PPT for output

19
include/convolution.hh
View file

@ -270,7 +270,24 @@ public:
auto grad22 = inr.gradient(d2r2[1],{i,j,k});
return (grad11*grad12-grad21*grad22)*inr.kelem(i, j, k);
},
// -> output operator
output_op);
}
/// @brief performs the multiplication of two fields by convolving them in Fourier space
/// @tparam opp output operator type
/// @param inl left input field a)
/// @param inr right input field b
/// @param output_op output operator
template <typename opp> // TOMA
void multiply_field(Grid_FFT<data_t> &inl, Grid_FFT<data_t> &inr, opp output_op)
{
// transform to FS in case fields are not
inl.FourierTransformForward();
inr.FourierTransformForward();
// perform convolution of Hessians
static_cast<derived_t &>(*this).convolve2(
[&inl](size_t i, size_t j, size_t k) -> ccomplex_t {return inl.kelem(i, j, k); },
[&inr](size_t i, size_t j, size_t k) -> ccomplex_t {return inr.kelem(i, j, k); },
output_op);
}
};

12
include/cosmology_calculator.hh
View file

@ -55,7 +55,7 @@ private:
interpolated_function_1d<true,true,false> D_of_a_, f_of_a_, a_of_D_;
double Dnow_, Dplus_start_, Dplus_target_, astart_, atarget_;
double m_n_s_, m_sqrtpnorm_;
double m_n_s_, m_sqrtpnorm_, tnorm_;
//! wrapper for GSL adaptive integration routine, do not use if many integrations need to be done as it allocates and deallocates memory
//! set to 61-point Gauss-Kronrod and large workspace, used for sigma_8 normalisation
@ -181,7 +181,7 @@ public:
music::ilog << "Linear growth factors: D+_target = " << Dplus_target_ << ", D+_start = " << Dplus_start_ << std::endl;
// set up transfer functions and compute normalisation
transfer_function_ = select_TransferFunction_plugin(cf, cosmo_param_);
transfer_function_ = std::move(select_TransferFunction_plugin(cf, cosmo_param_));
transfer_function_->intialise();
if( !transfer_function_->tf_isnormalised_ ){
cosmo_param_.set("pnorm", this->compute_pnorm_from_sigma8() );
@ -204,6 +204,9 @@ public:
m_n_s_ = cosmo_param_["n_s"];
m_sqrtpnorm_ = cosmo_param_["sqrtpnorm"];
double k_p = cosmo_param_["k_p"] / cosmo_param_["h"];
tnorm_ = std::sqrt(2.0 * M_PI * M_PI * cosmo_param_["A_s"] * std::pow(1.0 / k_p, cosmo_param_["n_s"] - 1) / std::pow(2.0 * M_PI, 3.0));
}
//! destructor
@ -404,6 +407,11 @@ public:
return std::pow(k, 0.5 * m_n_s_) * transfer_function_->compute(k, type) * m_sqrtpnorm_;
}
inline real_t get_transfer( const real_t k, const tf_type type) const
{
return -transfer_function_->compute(k, type)*k*k / tnorm_ * m_sqrtpnorm_;
}
//! Compute amplitude of the back-scaled delta_bc mode, with decaying velocity v_bc included or not (in which case delta_bc=const)
inline real_t get_amplitude_delta_bc( const real_t k, bool withvbc ) const
{

84
src/ic_generator.cc
View file

@ -28,6 +28,7 @@
#include <unistd.h> // for unlink
/**
* @brief the possible species of fluids
*
@ -62,9 +63,9 @@ std::unique_ptr<cosmology::calculator> the_cosmo_calc;
*/
int initialise( config_file& the_config )
{
the_random_number_generator = select_RNG_plugin(the_config);
the_random_number_generator = std::move(select_RNG_plugin(the_config));
the_cosmo_calc = std::make_unique<cosmology::calculator>(the_config);
the_output_plugin = select_output_plugin(the_config, the_cosmo_calc);
the_output_plugin = std::move(select_output_plugin(the_config, the_cosmo_calc));
return 0;
}
@ -107,7 +108,11 @@ int run( config_file& the_config )
//--------------------------------------------------------------------------------------------------------
//! order of the LPT approximation
const int LPTorder = the_config.get_value_safe<double>("setup","LPTorder",100);
const double fnl = the_config.get_value_safe<double>("cosmology","fnl",0);
const double nf = the_config.get_value_safe<double>("cosmology","nf",0);
const double k0 = the_config.get_value_safe<double>("cosmology","k0",0);
const double gnl = the_config.get_value_safe<double>("cosmology","gnl",0);
const double norm = the_config.get_value_safe<double>("cosmology","norm",1);
#if defined(USE_CONVOLVER_ORSZAG)
//! check if grid size is even for Orszag convolver
if( (ngrid%2 != 0) && (LPTorder>1) ){
@ -254,6 +259,8 @@ int run( config_file& the_config )
//... Next, declare LPT related arrays, allocated only as needed by order
Grid_FFT<real_t> phi({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
Grid_FFT<real_t> phi2({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}, false); // do not allocate these unless needed
Grid_FFT<real_t> delta_power({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}, false); // TOMA
Grid_FFT<real_t> phi3({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}, false); // ..
Grid_FFT<real_t> A3x({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}, false); // ..
Grid_FFT<real_t> A3y({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}, false); // ..
@ -380,12 +387,73 @@ int run( config_file& the_config )
music::ilog << std::setw(79) << std::setfill('.') << std::left << ">> Computing phi(1) term" << std::endl;
phi.FourierTransformForward(false);
phi.assign_function_of_grids_kdep([&](auto k, auto wn) {
real_t kmod = k.norm();
ccomplex_t delta = wn * the_cosmo_calc->get_amplitude(kmod, delta_matter);
return -delta / (kmod * kmod);
}, wnoise);
if (fnl != 0 || gnl != 0) {
phi.assign_function_of_grids_kdep([&](auto k, auto wn) {
real_t kmod = k.norm();
ccomplex_t zeta = wn * the_cosmo_calc->get_amplitude(kmod, delta_matter) / the_cosmo_calc->get_transfer(kmod, delta_matter);
return zeta; // zeta is temporarely stored in phi
}, wnoise);
phi.zero_DC_mode();
delta_power.allocate();
delta_power.FourierTransformForward(false);
Conv.multiply_field(phi, phi , op::assign_to(delta_power)); // phi2 = zeta^2
if (nf != 0)
{
delta_power.assign_function_of_grids_kdep([&](auto k, auto delta_power) {
real_t kmod = k.norm();
return std::pow(kmod/k0, nf)*delta_power;
}, delta_power);
}
delta_power.FourierTransformBackward();
phi.FourierTransformBackward();
if (fnl != 0)
{
music::ilog << "\n>>> Computing fnl term.... <<<\n" << std::endl;
phi.assign_function_of_grids_r([&](auto delta1, auto delta_power ){
return norm*(delta1 -fnl*delta_power*3.0/5.0) ;}, phi, delta_power);
// return norm*(delta1 - fnl*delta_power) ;}, phi, delta_power);
// the -3/5 factor is to match the usual fnl in terms of phi
// 3/5 fnl_zeta = fnl_phi
// no need to sustract the mean since latter
//is assured that the mean of delta is zero
}
else{
music::ilog << "\n>>> Computing gnl term.... <<<\n" << std::endl;
Conv.multiply_field(delta_power, phi , op::assign_to(delta_power)); // delta3 = delta^3
delta_power.FourierTransformBackward();
phi.FourierTransformBackward();
phi.assign_function_of_grids_r([&](auto delta1, auto delta_power ){
return norm*(delta1 - gnl*delta_power*9.0/25.0) ;}, phi, delta_power);
// the -9/25 factor is to match the usual gnl in terms of phi
// 9/25 gnl_zeta = gnl_phi
}
phi.FourierTransformForward();
phi.assign_function_of_grids_kdep([&](auto k, auto delta) {
real_t kmod = k.norm();
return - delta * the_cosmo_calc->get_transfer(kmod, delta_matter) / kmod /kmod ;
}, phi);
} else {
phi.assign_function_of_grids_kdep([&](auto k, auto wn) {
real_t kmod = k.norm();
ccomplex_t delta = wn * the_cosmo_calc->get_amplitude(kmod, delta_matter);
return -delta / (kmod * kmod);
}, wnoise);
}
phi.zero_DC_mode();
music::ilog << std::setw(70) << std::setfill(' ') << std::right << "took : " << std::setw(8) << get_wtime() - wtime << "s" << std::endl;