From 7e8f921be6019d8416e24df83ecc3ae47d2072cc Mon Sep 17 00:00:00 2001 From: adrigut10 Date: Tue, 9 Apr 2024 13:42:07 +0200 Subject: [PATCH 1/3] Test on pushing stuff to the online repo --- test | 1 + 1 file changed, 1 insertion(+) create mode 100644 test diff --git a/test b/test new file mode 100644 index 0000000..980a0d5 --- /dev/null +++ b/test @@ -0,0 +1 @@ +Hello World! From 6a677a08f92399b26e7b6e2b745f4b5f596024d4 Mon Sep 17 00:00:00 2001 From: adrigut10 Date: Tue, 9 Apr 2024 13:45:13 +0200 Subject: [PATCH 2/3] test removing a file --- test | 1 - 1 file changed, 1 deletion(-) delete mode 100644 test diff --git a/test b/test deleted file mode 100644 index 980a0d5..0000000 --- a/test +++ /dev/null @@ -1 +0,0 @@ -Hello World! From edb7c08b8611bc5ddb7bcd5ba1072018c03860bb Mon Sep 17 00:00:00 2001 From: adrigut10 Date: Mon, 15 Apr 2024 11:34:09 +0200 Subject: [PATCH 3/3] Added fNL and gNL part --- example.conf | 13 ++--- include/convolution.hh | 19 +++++++- include/cosmology_calculator.hh | 12 ++++- src/ic_generator.cc | 84 +++++++++++++++++++++++++++++---- 4 files changed, 111 insertions(+), 17 deletions(-) diff --git a/example.conf b/example.conf index 76f3a1d..e5a7db4 100644 --- a/example.conf +++ b/example.conf @@ -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 diff --git a/include/convolution.hh b/include/convolution.hh index 56675b4..1eff731 100644 --- a/include/convolution.hh +++ b/include/convolution.hh @@ -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 // TOMA + void multiply_field(Grid_FFT &inl, Grid_FFT &inr, opp output_op) + { + // transform to FS in case fields are not + inl.FourierTransformForward(); + inr.FourierTransformForward(); + // perform convolution of Hessians + static_cast(*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); } }; diff --git a/include/cosmology_calculator.hh b/include/cosmology_calculator.hh index b11b1d0..9565feb 100644 --- a/include/cosmology_calculator.hh +++ b/include/cosmology_calculator.hh @@ -55,7 +55,7 @@ private: interpolated_function_1d 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 { diff --git a/src/ic_generator.cc b/src/ic_generator.cc index cfd9003..c0c4f2a 100644 --- a/src/ic_generator.cc +++ b/src/ic_generator.cc @@ -28,6 +28,7 @@ #include // for unlink + /** * @brief the possible species of fluids * @@ -62,9 +63,9 @@ std::unique_ptr 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(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("setup","LPTorder",100); - + const double fnl = the_config.get_value_safe("cosmology","fnl",0); + const double nf = the_config.get_value_safe("cosmology","nf",0); + const double k0 = the_config.get_value_safe("cosmology","k0",0); + const double gnl = the_config.get_value_safe("cosmology","gnl",0); + const double norm = the_config.get_value_safe("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 phi({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); Grid_FFT phi2({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}, false); // do not allocate these unless needed + Grid_FFT delta_power({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}, false); // TOMA + Grid_FFT phi3({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}, false); // .. Grid_FFT A3x({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}, false); // .. Grid_FFT 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;