From d06e18fba8a314aa6fbd7bc405c7e23168f09096 Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Fri, 7 Jun 2024 13:23:33 +0200 Subject: [PATCH 1/4] allowed separate selection of convolver for PNG terms --- CMakeLists.txt | 11 ++++++++++- include/cmake_config.hh.in | 7 +++++++ src/ic_generator.cc | 19 +++++++++++++++---- src/main.cc | 3 ++- 4 files changed, 34 insertions(+), 6 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 1e091ee..a867d8a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -99,13 +99,22 @@ set_property ( # convolver type, right now only orszag or naive set ( CONVOLVER_TYPE "ORSZAG" - CACHE STRING "Convolution algorithm to be used (Naive=no dealiasing, Orszag=dealiased)" + CACHE STRING "Convolution algorithm to be used for LPT terms (Naive=no dealiasing, Orszag=dealiased)" ) set_property ( CACHE CONVOLVER_TYPE PROPERTY STRINGS ORSZAG NAIVE ) +set ( + CONVOLVER_TYPE_PNG "ORSZAG" + CACHE STRING "Convolution algorithm to be used for PNG terms (Naive=no dealiasing, Orszag=dealiased)" +) +set_property ( + CACHE CONVOLVER_TYPE_PNG + PROPERTY STRINGS ORSZAG NAIVE +) + ######################################################################################################################## # PLT options, right now only on/off option(ENABLE_PLT "Enable PLT (particle linear theory) corrections" OFF) diff --git a/include/cmake_config.hh.in b/include/cmake_config.hh.in index cc0da5e..9af13e8 100644 --- a/include/cmake_config.hh.in +++ b/include/cmake_config.hh.in @@ -20,6 +20,13 @@ constexpr char CMAKE_BUILDTYPE_STR[] = "${CMAKE_BUILD_TYPE}"; constexpr char CMAKE_CONVOLVER_STR[] = "Aliased"; #endif +#define USE_CONVOLVER_${CONVOLVER_TYPE_PNG}_PNG +#if defined(USE_CONVOLVER_ORSZAG_PNG) + constexpr char CMAKE_CONVOLVER_PNG_STR[] = "Orszag3/2"; +#elif defined(USE_CONVOLVER_NAIVE_PNG) + constexpr char CMAKE_CONVOLVER_PNG_STR[] = "Aliased"; +#endif + #if defined(ENABLE_PLT) constexpr char CMAKE_PLT_STR[] = "PLT corr. on"; #else diff --git a/src/ic_generator.cc b/src/ic_generator.cc index c0c4f2a..fffd33a 100644 --- a/src/ic_generator.cc +++ b/src/ic_generator.cc @@ -113,18 +113,22 @@ int run( config_file& the_config ) 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) + #if defined(USE_CONVOLVER_ORSZAG) | defined(USE_CONVOLVER_ORSZAG_PNG) //! check if grid size is even for Orszag convolver if( (ngrid%2 != 0) && (LPTorder>1) ){ music::elog << "ERROR: Orszag convolver for LPTorder>1 requires even grid size!" << std::endl; throw std::runtime_error("Orszag convolver for LPTorder>1 requires even grid size!"); return 0; } - #else + #elif !defined(USE_CONVOLVER_ORSZAG) //! warn if Orszag convolver is not used if( LPTorder>1 ){ music::wlog << "WARNING: LPTorder>1 requires USE_CONVOLVER_ORSZAG to be enabled to avoid aliased results!" << std::endl; } + #elif !defined(USE_CONVOLVER_ORSZAG_PNG) + if( fnl != 0 || gnl != 0 ){ + music::wlog << "WARNING: fNL!=0 or gNL!=0 requires USE_CONVOLVER_ORSZAG_PNG to be enabled to avoid aliased results!" << std::endl; + } #endif //-------------------------------------------------------------------------------------------------------- @@ -358,6 +362,13 @@ int run( config_file& the_config ) #elif defined(USE_CONVOLVER_NAIVE) NaiveConvolver Conv({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); #endif +#if defined(USE_CONVOLVER_ORSZAG_PNG) + OrszagConvolver Conv_png({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); +#elif defined(USE_CONVOLVER_NAIVE_PNG) + NaiveConvolver Conv_png({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); +#endif + + //-------------------------------------------------------------------- //-------------------------------------------------------------------- @@ -400,7 +411,7 @@ int run( config_file& the_config ) delta_power.allocate(); delta_power.FourierTransformForward(false); - Conv.multiply_field(phi, phi , op::assign_to(delta_power)); // phi2 = zeta^2 + Conv_png.multiply_field(phi, phi , op::assign_to(delta_power)); // phi2 = zeta^2 if (nf != 0) { @@ -428,7 +439,7 @@ int run( config_file& the_config ) music::ilog << "\n>>> Computing gnl term.... <<<\n" << std::endl; - Conv.multiply_field(delta_power, phi , op::assign_to(delta_power)); // delta3 = delta^3 + Conv_png.multiply_field(delta_power, phi , op::assign_to(delta_power)); // delta3 = delta^3 delta_power.FourierTransformBackward(); phi.FourierTransformBackward(); diff --git a/src/main.cc b/src/main.cc index 4668964..d26ba8e 100644 --- a/src/main.cc +++ b/src/main.cc @@ -125,7 +125,8 @@ int main( int argc, char** argv ) music::ilog << "-------------------------------------------------------------------------------" << std::endl; music::ilog << "Compile time options : " << std::endl; music::ilog << " Precision : " << CMAKE_PRECISION_STR << std::endl; - music::ilog << " Convolutions : " << CMAKE_CONVOLVER_STR << std::endl; + music::ilog << " Convolutions (LPT): " << CMAKE_CONVOLVER_STR << std::endl; + music::ilog << " Convolutions (PNG): " << CMAKE_CONVOLVER_PNG_STR << std::endl; music::ilog << " PLT : " << CMAKE_PLT_STR << std::endl; music::ilog << "-------------------------------------------------------------------------------" << std::endl; From 03e5fb5e8435b624fb9ff7a9f2ea0cd6e334cd9f Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Tue, 11 Jun 2024 16:34:29 +0200 Subject: [PATCH 2/4] fixed subtraction of phi in gNL term --- src/ic_generator.cc | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/ic_generator.cc b/src/ic_generator.cc index fffd33a..562acbf 100644 --- a/src/ic_generator.cc +++ b/src/ic_generator.cc @@ -423,12 +423,13 @@ int run( config_file& the_config ) delta_power.FourierTransformBackward(); phi.FourierTransformBackward(); + real_t var_phi = delta_power.mean(); // = 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 - var_phi)*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 @@ -442,10 +443,11 @@ int run( config_file& the_config ) Conv_png.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); + return norm*(delta1 - gnl*(delta_power - 3*var_phi*delta1)*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 } From 7598bcc07f5c77405b18377398ef934b859608e0 Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Mon, 8 Jul 2024 23:26:06 +0200 Subject: [PATCH 3/4] more output columns for input_powerspec --- include/cosmology_calculator.hh | 72 +++++++++++++++++++++------------ 1 file changed, 47 insertions(+), 25 deletions(-) diff --git a/include/cosmology_calculator.hh b/include/cosmology_calculator.hh index 9565feb..dd37cc0 100644 --- a/include/cosmology_calculator.hh +++ b/include/cosmology_calculator.hh @@ -222,40 +222,62 @@ public: if (CONFIG::MPI_task_rank == 0) { double kmin = std::max(1e-4, transfer_function_->get_kmin()); + double fb = cosmo_param_["f_b"], fc = cosmo_param_["f_c"]; // write power spectrum to a file std::ofstream ofs(fname.c_str()); - std::stringstream ss; - ss << " ,ap=" << a << ""; + // std::stringstream ss; + // ss << " ,ap=" << a << ""; + ofs << "# astart = " << a << std::endl; + ofs << "# atarget = " << atarget_ << std::endl; + ofs << "# D+start = " << this->get_growth_factor(a) << std::endl; ofs << "# " << std::setw(18) << "k [h/Mpc]" - << std::setw(20) << ("P_dtot(k,a=ap)") - << std::setw(20) << ("P_dcdm(k,a=ap)") - << std::setw(20) << ("P_dbar(k,a=ap)") - << std::setw(20) << ("P_tcdm(k,a=ap)") - << std::setw(20) << ("P_tbar(k,a=ap)") - << std::setw(20) << ("P_dtot(k,a=1)") - << std::setw(20) << ("P_dcdm(k,a=1)") - << std::setw(20) << ("P_dbar(k,a=1)") - << std::setw(20) << ("P_tcdm(k,a=1)") - << std::setw(20) << ("P_tbar(k,a=1)") + << std::setw(20) << ("P_deltam") + << std::setw(20) << ("P_deltac") + << std::setw(20) << ("P_deltab") + << std::setw(20) << ("P_deltabc") + << std::setw(20) << ("P_thetam") + << std::setw(20) << ("P_thetac") + << std::setw(20) << ("P_thetab") + << std::setw(20) << ("P_thetabc") << std::endl; + for (double k = kmin; k < transfer_function_->get_kmax(); k *= 1.01) { + const double dm = this->get_amplitude(k, delta_matter) * Dplus_start_ / Dplus_target_; + const double dbc = this->get_amplitude(k, delta_bc); + const double db = dm + fc * dbc; + const double dc = dm - fb * dbc; + const double tm = this->get_amplitude(k, delta_matter) * Dplus_start_ / Dplus_target_; + const double tbc = this->get_amplitude(k, theta_bc); + const double tb = dm + fc * dbc; + const double tc = dm - fb * dbc; + ofs << std::setw(20) << std::setprecision(10) << k - << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_matter)*Dplus_start_, 2.0) - << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_cdm)*Dplus_start_, 2.0) - << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_baryon)*Dplus_start_, 2.0) - // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_matter)*Dplus_start_, 2.0) - // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_cdm)*Dplus_start_, 2.0) - // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_baryon)*Dplus_start_, 2.0) - // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_cdm)*Dplus_start_, 2.0) - // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_baryon)*Dplus_start_, 2.0) - // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_matter0)* Dplus_start_ / Dplus_target_, 2.0) - // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_cdm0)* Dplus_start_ / Dplus_target_, 2.0) - // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_baryon0)* Dplus_start_ / Dplus_target_, 2.0) - // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_cdm0)* Dplus_start_ / Dplus_target_, 2.0) - // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_baryon0)* Dplus_start_ / Dplus_target_, 2.0) + << std::setw(20) << std::setprecision(10) << std::pow(dm,2) + << std::setw(20) << std::setprecision(10) << std::pow(dc,2) + << std::setw(20) << std::setprecision(10) << std::pow(db,2) + << std::setw(20) << std::setprecision(10) << std::pow(dbc + 2 * tbc * (std::sqrt( Dplus_target_ / Dplus_start_ ) - 1.0),2) + << std::setw(20) << std::setprecision(10) << std::pow(tm / std::pow( Dplus_start_ / Dplus_target_, 0.5 ),2) + << std::setw(20) << std::setprecision(10) << std::pow(tc / std::pow( Dplus_start_ / Dplus_target_, 0.5 ),2) + << std::setw(20) << std::setprecision(10) << std::pow(tb / std::pow( Dplus_start_ / Dplus_target_, 0.5 ),2) + << std::setw(20) << std::setprecision(10) << std::pow(tbc / std::pow( Dplus_start_ / Dplus_target_, 0.5 ),2) << std::endl; + // ofs << std::setw(20) << std::setprecision(10) << k + // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_matter)*Dplus_start_, 2.0) + // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_cdm)*Dplus_start_, 2.0) + // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_baryon)*Dplus_start_, 2.0) + // // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_matter)*Dplus_start_, 2.0) + // // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_cdm)*Dplus_start_, 2.0) + // // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_baryon)*Dplus_start_, 2.0) + // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_cdm)*Dplus_start_, 2.0) + // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_baryon)*Dplus_start_, 2.0) + // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_matter0)* Dplus_start_ / Dplus_target_, 2.0) + // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_cdm0)* Dplus_start_ / Dplus_target_, 2.0) + // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_baryon0)* Dplus_start_ / Dplus_target_, 2.0) + // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_cdm0)* Dplus_start_ / Dplus_target_, 2.0) + // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_baryon0)* Dplus_start_ / Dplus_target_, 2.0) + // << std::endl; } } music::ilog << "Wrote power spectrum at a=" << a << " to file \'" << fname << "\'" << std::endl; From f64cc6d2e142bb14b3d60dbdc1ac790e9a8e9bfc Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Thu, 22 Aug 2024 11:06:10 +0200 Subject: [PATCH 4/4] fixed wrong normalisation of input_powerspec --- include/cosmology_calculator.hh | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/include/cosmology_calculator.hh b/include/cosmology_calculator.hh index dd37cc0..7ac9a5c 100644 --- a/include/cosmology_calculator.hh +++ b/include/cosmology_calculator.hh @@ -235,33 +235,33 @@ public: << std::setw(20) << ("P_deltam") << std::setw(20) << ("P_deltac") << std::setw(20) << ("P_deltab") - << std::setw(20) << ("P_deltabc") - << std::setw(20) << ("P_thetam") - << std::setw(20) << ("P_thetac") - << std::setw(20) << ("P_thetab") - << std::setw(20) << ("P_thetabc") + // << std::setw(20) << ("P_deltabc") + // << std::setw(20) << ("P_thetam") + // << std::setw(20) << ("P_thetac") + // << std::setw(20) << ("P_thetab") + // << std::setw(20) << ("P_thetabc") << std::endl; for (double k = kmin; k < transfer_function_->get_kmax(); k *= 1.01) { - const double dm = this->get_amplitude(k, delta_matter) * Dplus_start_ / Dplus_target_; + const double dm = this->get_amplitude(k, delta_matter) * Dplus_start_; // / Dplus_target_; const double dbc = this->get_amplitude(k, delta_bc); const double db = dm + fc * dbc; const double dc = dm - fb * dbc; - const double tm = this->get_amplitude(k, delta_matter) * Dplus_start_ / Dplus_target_; - const double tbc = this->get_amplitude(k, theta_bc); - const double tb = dm + fc * dbc; - const double tc = dm - fb * dbc; + // const double tm = this->get_amplitude(k, theta_matter) * Dplus_start_; // / Dplus_target_; + // const double tbc = this->get_amplitude(k, theta_bc); + // const double tb = dm + fc * dbc; + // const double tc = dm - fb * dbc; ofs << std::setw(20) << std::setprecision(10) << k << std::setw(20) << std::setprecision(10) << std::pow(dm,2) << std::setw(20) << std::setprecision(10) << std::pow(dc,2) << std::setw(20) << std::setprecision(10) << std::pow(db,2) - << std::setw(20) << std::setprecision(10) << std::pow(dbc + 2 * tbc * (std::sqrt( Dplus_target_ / Dplus_start_ ) - 1.0),2) - << std::setw(20) << std::setprecision(10) << std::pow(tm / std::pow( Dplus_start_ / Dplus_target_, 0.5 ),2) - << std::setw(20) << std::setprecision(10) << std::pow(tc / std::pow( Dplus_start_ / Dplus_target_, 0.5 ),2) - << std::setw(20) << std::setprecision(10) << std::pow(tb / std::pow( Dplus_start_ / Dplus_target_, 0.5 ),2) - << std::setw(20) << std::setprecision(10) << std::pow(tbc / std::pow( Dplus_start_ / Dplus_target_, 0.5 ),2) + // << std::setw(20) << std::setprecision(10) << std::pow(dbc + 2 * tbc * (std::sqrt( Dplus_target_ / Dplus_start_ ) - 1.0),2) + // << std::setw(20) << std::setprecision(10) << std::pow(tm / std::pow( Dplus_start_ / Dplus_target_, 0.5 ),2) + // << std::setw(20) << std::setprecision(10) << std::pow(tc / std::pow( Dplus_start_ / Dplus_target_, 0.5 ),2) + // << std::setw(20) << std::setprecision(10) << std::pow(tb / std::pow( Dplus_start_ / Dplus_target_, 0.5 ),2) + // << std::setw(20) << std::setprecision(10) << std::pow(tbc / std::pow( Dplus_start_ / Dplus_target_, 0.5 ),2) << std::endl; // ofs << std::setw(20) << std::setprecision(10) << k // << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_matter)*Dplus_start_, 2.0)