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/include/cosmology_calculator.hh b/include/cosmology_calculator.hh index 9565feb..7ac9a5c 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, 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(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; diff --git a/src/ic_generator.cc b/src/ic_generator.cc index 3064355..f1f9ba9 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) { @@ -412,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 @@ -428,13 +440,14 @@ 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(); 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 } diff --git a/src/main.cc b/src/main.cc index db319dc..9fe9de3 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;