diff --git a/example.conf b/example.conf index 6b57a7c..db8fa73 100644 --- a/example.conf +++ b/example.conf @@ -30,18 +30,18 @@ ParticleLoad = sc # particle load, can be 'sc' (1x), 'bcc' (2x) or 'fcc [cosmology] ## transfer = ... specifies the Einstein-Boltzmann plugin module -# main cosmological parameters -Omega_m = 0.3158 -Omega_b = 0.0494 -Omega_L = 0.6842 -H0 = 67.321 -n_s = 0.9661 -sigma_8 = 0.8102 +ParameterSet = Planck2018EE+BAO+SN # specify a pre-defined parameter set, or set to 'none' and set manually below -ZeroRadiation = false # For Back-scaling only: set to true if your simulation code - # cannot deal with Omega_r!=0 in its background FLRW model - -## Additional cosmological parameters (set by default to the given values) +## cosmological parameters, to set, choose ParameterSet = none, +## default values (those not specified) are set to the values +## from 'Planck2018EE+BAO+SN', we currently assume flatness +# Omega_m = 0.3158 +# Omega_b = 0.0494 +# Omega_L = 0.6842 +# H0 = 67.321 +# n_s = 0.9661 +# sigma_8 = 0.8102 +# A_s = 2.148752e-09 # can use A_s instead of sigma_8 when using CLASS # Tcmb = 2.7255 # k_p = 0.05 # N_ur = 2.046 @@ -51,6 +51,8 @@ ZeroRadiation = false # For Back-scaling only: set to true if your simulation # w_0 = -1.0 # not supported yet! # w_a = 0.0 # not supported yet! +ZeroRadiation = false # For Back-scaling only: set to true if your simulation code + # cannot deal with Omega_r!=0 in its background FLRW model ## Use below for anisotropic large scale tidal field ICs up to 2LPT ## see Stuecker+2020 (https://arxiv.org/abs/2003.06427) @@ -78,7 +80,6 @@ ZeroRadiation = false # For Back-scaling only: set to true if your simulation transfer = CLASS ztarget = 2.5 # target redshift for CLASS module, output at ztarget will be back-scaled to zstart -# A_s = 2.148752e-09 # can use A_s instead of sigma_8 when using CLASS ######################################################################################### diff --git a/include/cosmology_parameters.hh b/include/cosmology_parameters.hh index 0505893..5a46b95 100644 --- a/include/cosmology_parameters.hh +++ b/include/cosmology_parameters.hh @@ -21,15 +21,21 @@ #include #include +#include namespace cosmology { - //! singleton structure for cosmological parameters + //! structure for cosmological parameters class parameters { + public: + using defaultmmap_t = std::map>; + private: std::map pmap_; //!< All parameters are stored here as key-value pairs + static defaultmmap_t default_pmaps_; //!< Holds pre-defined parameter sets, see src/cosmology_parameters.cc + public: //! get routine for cosmological parameter key-value pairs double get(const std::string &key) const @@ -63,8 +69,8 @@ namespace cosmology //! shortcut get routine for cosmological parameter key-value pairs through bracket operator inline double operator[](const std::string &key) const { return this->get(key); } - //! no default constructor - parameters() = delete; + //! default constructor does nothing + parameters() {} //! default copy constructor parameters(const parameters &) = default; @@ -72,44 +78,87 @@ namespace cosmology //! main constructor for explicit construction from input config file explicit parameters( config_file &cf ) { - // CMB - pmap_["Tcmb"] = cf.get_value_safe("cosmology", "Tcmb", 2.7255); - pmap_["YHe"] = cf.get_value_safe("cosmology", "YHe", 0.2454006); - - // H0 - pmap_["H0"] = cf.get_value("cosmology", "H0"); - pmap_["h"] = cf.get_value("cosmology", "H0") / 100.0; - const double h = pmap_["h"]; - - // primordial and normalisation - if(!cf.contains_key("cosmology/n_s")) - pmap_["n_s"] = cf.get_value("cosmology", "nspec"); - else - pmap_["n_s"] = cf.get_value("cosmology", "n_s"); - - pmap_["A_s"] = cf.get_value_safe("cosmology", "A_s", -1.0); - pmap_["k_p"] = cf.get_value_safe("cosmology", "k_p", 0.05); - - pmap_["sigma_8"] = cf.get_value_safe("cosmology", "sigma_8", -1.0); + music::ilog << "-------------------------------------------------------------------------------" << std::endl; - // baryon and non-relativistic matter content - pmap_["Omega_b"] = cf.get_value("cosmology", "Omega_b"); - pmap_["Omega_m"] = cf.get_value("cosmology", "Omega_m"); + if( cf.get_value_safe("cosmology","ParameterSet","none") == std::string("none")) + { + //------------------------------------------------------------------------------- + // read parameters from config file + //------------------------------------------------------------------------------- - // massive neutrino species - pmap_["m_nu1"] = cf.get_value_safe("cosmology", "m_nu1", 0.06); - pmap_["m_nu2"] = cf.get_value_safe("cosmology", "m_nu2", 0.0); - pmap_["m_nu3"] = cf.get_value_safe("cosmology", "m_nu3", 0.0); + auto defaultp = default_pmaps_.find("Planck2018EE+BAO+SN")->second; + + // CMB + pmap_["Tcmb"] = cf.get_value_safe("cosmology", "Tcmb", defaultp["Tcmb"]); + pmap_["YHe"] = cf.get_value_safe("cosmology", "YHe", defaultp["YHe"]); + + // H0 + pmap_["h"] = cf.get_value_safe("cosmology", "H0", defaultp["h"]*100.0) / 100.0; + + // primordial and normalisation + if(!cf.contains_key("cosmology/n_s")) + pmap_["n_s"] = cf.get_value_safe("cosmology", "nspec", defaultp["n_s"]); + else + pmap_["n_s"] = cf.get_value_safe("cosmology", "n_s", defaultp["n_s"]); + + pmap_["A_s"] = cf.get_value_safe("cosmology", "A_s", cf.contains_key("cosmology/sigma_8")? -1.0 : defaultp["A_s"]); + pmap_["k_p"] = cf.get_value_safe("cosmology", "k_p", defaultp["k_p"]); + + pmap_["sigma_8"] = cf.get_value_safe("cosmology", "sigma_8", -1.0); + + // baryon and non-relativistic matter content + pmap_["Omega_b"] = cf.get_value_safe("cosmology", "Omega_b", defaultp["Omega_b"]); + pmap_["Omega_m"] = cf.get_value_safe("cosmology", "Omega_m", defaultp["Omega_m"]); + + // massive neutrino species + pmap_["m_nu1"] = cf.get_value_safe("cosmology", "m_nu1", defaultp["m_nu1"]); + pmap_["m_nu2"] = cf.get_value_safe("cosmology", "m_nu2", defaultp["m_nu2"]); + pmap_["m_nu3"] = cf.get_value_safe("cosmology", "m_nu3", defaultp["m_nu3"]); + int N_nu_massive = int(this->get("m_nu1") > 1e-9) + int(this->get("m_nu2") > 1e-9) + int(this->get("m_nu3") > 1e-9);; + + // number ultrarelativistic neutrinos + pmap_["N_ur"] = cf.get_value_safe("cosmology", "N_ur", 3.046 - N_nu_massive); + + // dark energy + pmap_["Omega_DE"] = cf.get_value_safe("cosmology", "Omega_L", defaultp["Omega_DE"]); + pmap_["w_0"] = cf.get_value_safe("cosmology", "w_0", defaultp["w_0"]); + pmap_["w_a"] = cf.get_value_safe("cosmology", "w_a", defaultp["w_a"]); + + }else{ + + //------------------------------------------------------------------------------- + // load predefined parameter set + //------------------------------------------------------------------------------- + auto pset_name = cf.get_value("cosmology","ParameterSet"); + auto it = default_pmaps_.find( pset_name ); + if( it == default_pmaps_.end() ){ + music::elog << "Unknown cosmology parameter set \'" << pset_name << "\'!" << std::endl; + music::ilog << "Valid pre-defined sets are: " << std::endl; + for( auto ii : default_pmaps_ ){ + music::ilog << " " << ii.first << std::endl; + } + throw std::runtime_error("Invalid value for cosmology/ParameterSet"); + }else{ + music::ilog << "Loading cosmological parameter set \'" << it->first << "\'..." << std::endl; + } + + // copy pre-defined set in + for( auto entry : it->second ){ + pmap_[entry.first] = entry.second; + } + } + + //------------------------------------------------------------------------------- + // derived parameters + //------------------------------------------------------------------------------- + double h = this->get("h"); + pmap_["H0"] = 100.0 * h; + + // massive neutrinos pmap_["N_nu_massive"] = int(this->get("m_nu1") > 1e-9) + int(this->get("m_nu2") > 1e-9) + int(this->get("m_nu3") > 1e-9); const double sum_m_nu = this->get("m_nu1") + this->get("m_nu2") + this->get("m_nu3"); - - // number ultrarelativistic neutrinos - pmap_["N_ur"] = cf.get_value_safe("cosmology", "N_ur", 3.046 - this->get("N_nu_massive")); pmap_["Omega_nu_massive"] = sum_m_nu / (93.14 * h * h); // Omega_nu_m = \sum_i m_i / (93.14 eV h^2) - // compute amount of cold dark matter as the rest - pmap_["Omega_c"] = this->get("Omega_m") - this->get("Omega_b") - this->get("Omega_nu_massive"); - // calculate energy density in ultrarelativistic species from Tcmb and Neff // photons pmap_["Omega_gamma"] = 4 * phys_const::sigma_SI / std::pow(phys_const::c_SI, 3) * std::pow(this->get("Tcmb"), 4.0) @@ -119,10 +168,8 @@ namespace cosmology // total relativistic pmap_["Omega_r"] = this->get("Omega_gamma") + this->get("Omega_nu_massless"); - // dark energy - pmap_["Omega_DE"] = cf.get_value("cosmology", "Omega_L"); - pmap_["w_0"] = cf.get_value_safe("cosmology", "w_0", -1.0); - pmap_["w_a"] = cf.get_value_safe("cosmology", "w_a", 0.0); + // compute amount of cold dark matter as the rest + pmap_["Omega_c"] = this->get("Omega_m") - this->get("Omega_b") - this->get("Omega_nu_massive"); if (cf.get_value_safe("cosmology", "ZeroRadiation", false)) { @@ -147,7 +194,6 @@ namespace cosmology pmap_["sqrtpnorm"] = 0.0; pmap_["vfact"] = 0.0; - music::ilog << "-------------------------------------------------------------------------------" << std::endl; music::ilog << "Cosmological parameters are: " << std::endl; music::ilog << " h = " << std::setw(16) << this->get("h"); if( this->get("A_s") > 0.0 ) @@ -156,7 +202,7 @@ namespace cosmology music::ilog << "sigma_8 = " << std::setw(16) << this->get("sigma_8"); music::ilog << "n_s = " << std::setw(16) << this->get("n_s") << std::endl; music::ilog << " Omega_c = " << std::setw(16) << this->get("Omega_c") << "Omega_b = " << std::setw(16) << this->get("Omega_b") << "Omega_m = " << std::setw(16) << this->get("Omega_m") << std::endl; - music::ilog << " Omega_r = " << std::setw(16) << this->get("Omega_r") << "Omega_nu = " << std::setw(16) << this->get("Omega_nu_massive") << "∑m_nu = " << std::setw(11) << sum_m_nu << "eV" << std::endl; + music::ilog << " Omega_r = " << std::setw(16) << this->get("Omega_r") << "Omega_nu = " << std::setw(16) << this->get("Omega_nu_massive") << "∑m_nu = " << sum_m_nu << "eV" << std::endl; music::ilog << " Omega_DE = " << std::setw(16) << this->get("Omega_DE") << "w_0 = " << std::setw(16) << this->get("w_0") << "w_a = " << std::setw(16) << this->get("w_a") << std::endl; //music::ilog << " Omega_k = " << 1.0 - this->get("Omega_m") - this->get("Omega_r") - this->get("Omega_DE") << std::endl; if (this->get("Omega_r") > 0.0) @@ -165,5 +211,15 @@ namespace cosmology music::wlog << " Make sure your sim code supports this, otherwise set [cosmology] / ZeroRadiation=true." << std::endl; } } + + //! print the names of all pre-defined parameter sets + void print_available_sets( void ) const + { + for( auto ii : default_pmaps_ ){ + music::ilog << "\t\'" << ii.first << "\'" << std::endl; + } + } }; + + void print_ParameterSets( void ); } // namespace cosmology \ No newline at end of file diff --git a/include/particle_plt.hh b/include/particle_plt.hh index cb4946f..2216340 100644 --- a/include/particle_plt.hh +++ b/include/particle_plt.hh @@ -51,7 +51,7 @@ class lattice_gradient{ private: const real_t boxlen_, aini_; const size_t ngmapto_, ngrid_, ngrid32_; - const real_t mapratio_, XmL_; + const real_t mapratio_;//, XmL_; Grid_FFT D_xx_, D_xy_, D_xz_, D_yy_, D_yz_, D_zz_; Grid_FFT grad_x_, grad_y_, grad_z_; std::vector> vectk_; @@ -522,7 +522,7 @@ public: aini_ ( 1.0/(1.0+the_config.get_value("setup", "zstart")) ), ngmapto_( the_config.get_value("setup", "GridRes") ), ngrid_( ngridself ), ngrid32_( std::pow(ngrid_, 1.5) ), mapratio_(real_t(ngrid_)/real_t(ngmapto_)), - XmL_ ( the_config.get_value("cosmology", "Omega_L") / the_config.get_value("cosmology", "Omega_m") ), + //XmL_ ( the_config.get_value("cosmology", "Omega_L") / the_config.get_value("cosmology", "Omega_m") ), D_xx_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}), D_xy_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}), D_xz_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}), D_yy_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}), D_yz_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}), D_zz_({ngrid_, ngrid_, ngrid_}, {1.0,1.0,1.0}), diff --git a/src/cosmology_parameters.cc b/src/cosmology_parameters.cc new file mode 100644 index 0000000..088c967 --- /dev/null +++ b/src/cosmology_parameters.cc @@ -0,0 +1,111 @@ +// This file is part of monofonIC (MUSIC2) +// A software package to generate ICs for cosmological simulations +// Copyright (C) 2020 by Oliver Hahn +// +// monofonIC is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// monofonIC is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see . + +#include + +namespace cosmology{ + +//! we store here the preset cosmological paramters +parameters::defaultmmap_t parameters::default_pmaps_ +{ + //============================================================================= + // Planck 2018 baseline cosmologies + // cf. https://wiki.cosmos.esa.int/planck-legacy-archive/images/b/be/Baseline_params_table_2018_68pc.pdf + //============================================================================= + + // baseline 2.17 base_plikHM_TTTEEE_lowl_lowE_lensing + {"Planck2018EE", { + {"h", 0.67321}, + {"Omega_m", 0.3158}, + {"Omega_b", 0.04938898}, + {"Omega_DE", 0.6842}, + {"w_0", -1.0}, + {"w_a", 0.0}, + {"n_s", 0.96605}, + {"A_s", 2.1005e-9}, + {"k_p", 0.05}, + {"YHe", 0.245401}, + {"N_ur", 2.046}, + {"m_nu1", 0.06}, + {"m_nu2", 0.0}, + {"m_nu3", 0.0}, + {"Tcmb", 2.7255}}}, + + // baseline 2.18 base_plikHM_TTTEEE_lowl_lowE_lensing_post_BAO + {"Planck2018EE+BAO", { + {"h", 0.67702}, + {"Omega_m", 0.3106}, + {"Omega_b", 0.04897284}, + {"Omega_DE", 0.6894}, + {"w_0", -1.0}, + {"w_a", 0.0}, + {"n_s", 0.96824}, + {"A_s", 2.1073e-9}, + {"k_p", 0.05}, + {"YHe", 0.245425}, + {"N_ur", 2.046}, + {"m_nu1", 0.06}, + {"m_nu2", 0.0}, + {"m_nu3", 0.0}, + {"Tcmb", 2.7255}}}, + + // baseline 2.19 base_plikHM_TTTEEE_lowl_lowE_lensing_post_Pantheon + {"Planck2018EE+SN", { + {"h", 0.6749}, + {"Omega_m", 0.3134}, + {"Omega_b", 0.04919537}, + {"Omega_DE", 0.6866}, + {"w_0", -1.0}, + {"w_a", 0.0}, + {"n_s", 0.96654}, + {"A_s", 2.1020e-9}, + {"k_p", 0.05}, + {"YHe", 0.245411}, + {"N_ur", 2.046}, + {"m_nu1", 0.06}, + {"m_nu2", 0.0}, + {"m_nu3", 0.0}, + {"Tcmb", 2.7255}}}, + + // baseline 2.20 base_plikHM_TTTEEE_lowl_lowE_lensing_post_BAO_Pantheon + {"Planck2018EE+BAO+SN", { + {"h", 0.67742}, + {"Omega_m", 0.3099}, + {"Omega_b", 0.048891054}, + {"Omega_DE", 0.6901}, + {"w_0", -1.0}, + {"w_a", 0.0}, + {"n_s", 0.96822}, + {"A_s", 2.1064e-9}, + {"k_p", 0.05}, + {"YHe", 0.245421}, + {"N_ur", 2.046}, + {"m_nu1", 0.06}, + {"m_nu2", 0.0}, + {"m_nu3", 0.0}, + {"Tcmb", 2.7255}}} +}; + + +void print_ParameterSets( void ){ + music::ilog << "Available cosmology parameter sets:" << std::endl; + cosmology::parameters p; + p.print_available_sets(); + music::ilog << std::endl; +} + +}// end namespace cosmology \ No newline at end of file diff --git a/src/ic_generator.cc b/src/ic_generator.cc index 4a63648..1088956 100644 --- a/src/ic_generator.cc +++ b/src/ic_generator.cc @@ -101,12 +101,12 @@ int run( config_file& the_config ) // compute mass fractions std::map< cosmo_species, double > Omega; if( bDoBaryons ){ - double Om = the_config.get_value("cosmology", "Omega_m"); - double Ob = the_config.get_value("cosmology", "Omega_b"); + double Om = the_cosmo_calc->cosmo_param_["Omega_m"]; + double Ob = the_cosmo_calc->cosmo_param_["Omega_b"]; Omega[cosmo_species::dm] = Om-Ob; Omega[cosmo_species::baryon] = Ob; }else{ - double Om = the_config.get_value("cosmology", "Omega_m"); + double Om = the_cosmo_calc->cosmo_param_["Omega_m"]; Omega[cosmo_species::dm] = Om; Omega[cosmo_species::baryon] = 0.0; } diff --git a/src/main.cc b/src/main.cc index fb22277..1b10ca7 100644 --- a/src/main.cc +++ b/src/main.cc @@ -28,6 +28,7 @@ #include #include +#include #include @@ -123,6 +124,7 @@ int main( int argc, char** argv ) if (argc != 2) { // print_region_generator_plugins(); + cosmology::print_ParameterSets(); print_TransferFunction_plugins(); print_RNG_plugins(); print_output_plugins(); diff --git a/src/plugins/output_arepo.cc b/src/plugins/output_arepo.cc index b8bd6c8..b847608 100644 --- a/src/plugins/output_arepo.cc +++ b/src/plugins/output_arepo.cc @@ -118,7 +118,7 @@ public: // initial gas temperature double Tcmb0 = pcc->cosmo_param_["Tcmb"]; double Omegab = pcc->cosmo_param_["Omega_b"]; - double h = cf_.get_value("cosmology", "H0") / 100.0, h2 = h*h; + double h = pcc->cosmo_param_["h"], h2 = h*h; double adec = 1.0 / (160.0 * pow(Omegab * h2 / 0.022, 2.0 / 5.0)); Tini_ = astart < adec ? Tcmb0 / astart : Tcmb0 / astart / astart * adec;