1
0
Fork 0
mirror of https://github.com/cosmo-sims/monofonIC.git synced 2024-09-19 17:03:45 +02:00

added pre-defined sets of cosmological parameters

This commit is contained in:
Oliver Hahn 2020-09-10 23:08:00 +02:00
parent 4ecc025ec7
commit db5ce865a4
7 changed files with 229 additions and 59 deletions

View file

@ -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
#########################################################################################

View file

@ -21,15 +21,21 @@
#include <physical_constants.hh>
#include <config_file.hh>
#include <general.hh>
namespace cosmology
{
//! singleton structure for cosmological parameters
//! structure for cosmological parameters
class parameters
{
public:
using defaultmmap_t = std::map<std::string,std::map<std::string,real_t>>;
private:
std::map<std::string, double> 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<double>("cosmology", "Tcmb", 2.7255);
pmap_["YHe"] = cf.get_value_safe<double>("cosmology", "YHe", 0.2454006);
// H0
pmap_["H0"] = cf.get_value<double>("cosmology", "H0");
pmap_["h"] = cf.get_value<double>("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<double>("cosmology", "nspec");
else
pmap_["n_s"] = cf.get_value<double>("cosmology", "n_s");
pmap_["A_s"] = cf.get_value_safe<double>("cosmology", "A_s", -1.0);
pmap_["k_p"] = cf.get_value_safe<double>("cosmology", "k_p", 0.05);
pmap_["sigma_8"] = cf.get_value_safe<double>("cosmology", "sigma_8", -1.0);
music::ilog << "-------------------------------------------------------------------------------" << std::endl;
// baryon and non-relativistic matter content
pmap_["Omega_b"] = cf.get_value<double>("cosmology", "Omega_b");
pmap_["Omega_m"] = cf.get_value<double>("cosmology", "Omega_m");
if( cf.get_value_safe<std::string>("cosmology","ParameterSet","none") == std::string("none"))
{
//-------------------------------------------------------------------------------
// read parameters from config file
//-------------------------------------------------------------------------------
// massive neutrino species
pmap_["m_nu1"] = cf.get_value_safe<double>("cosmology", "m_nu1", 0.06);
pmap_["m_nu2"] = cf.get_value_safe<double>("cosmology", "m_nu2", 0.0);
pmap_["m_nu3"] = cf.get_value_safe<double>("cosmology", "m_nu3", 0.0);
auto defaultp = default_pmaps_.find("Planck2018EE+BAO+SN")->second;
// CMB
pmap_["Tcmb"] = cf.get_value_safe<double>("cosmology", "Tcmb", defaultp["Tcmb"]);
pmap_["YHe"] = cf.get_value_safe<double>("cosmology", "YHe", defaultp["YHe"]);
// H0
pmap_["h"] = cf.get_value_safe<double>("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<double>("cosmology", "nspec", defaultp["n_s"]);
else
pmap_["n_s"] = cf.get_value_safe<double>("cosmology", "n_s", defaultp["n_s"]);
pmap_["A_s"] = cf.get_value_safe<double>("cosmology", "A_s", cf.contains_key("cosmology/sigma_8")? -1.0 : defaultp["A_s"]);
pmap_["k_p"] = cf.get_value_safe<double>("cosmology", "k_p", defaultp["k_p"]);
pmap_["sigma_8"] = cf.get_value_safe<double>("cosmology", "sigma_8", -1.0);
// baryon and non-relativistic matter content
pmap_["Omega_b"] = cf.get_value_safe<double>("cosmology", "Omega_b", defaultp["Omega_b"]);
pmap_["Omega_m"] = cf.get_value_safe<double>("cosmology", "Omega_m", defaultp["Omega_m"]);
// massive neutrino species
pmap_["m_nu1"] = cf.get_value_safe<double>("cosmology", "m_nu1", defaultp["m_nu1"]);
pmap_["m_nu2"] = cf.get_value_safe<double>("cosmology", "m_nu2", defaultp["m_nu2"]);
pmap_["m_nu3"] = cf.get_value_safe<double>("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<double>("cosmology", "N_ur", 3.046 - N_nu_massive);
// dark energy
pmap_["Omega_DE"] = cf.get_value_safe<double>("cosmology", "Omega_L", defaultp["Omega_DE"]);
pmap_["w_0"] = cf.get_value_safe<double>("cosmology", "w_0", defaultp["w_0"]);
pmap_["w_a"] = cf.get_value_safe<double>("cosmology", "w_a", defaultp["w_a"]);
}else{
//-------------------------------------------------------------------------------
// load predefined parameter set
//-------------------------------------------------------------------------------
auto pset_name = cf.get_value<std::string>("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<double>("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<double>("cosmology", "Omega_L");
pmap_["w_0"] = cf.get_value_safe<double>("cosmology", "w_0", -1.0);
pmap_["w_a"] = cf.get_value_safe<double>("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<bool>("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

View file

@ -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<real_t,false> D_xx_, D_xy_, D_xz_, D_yy_, D_yz_, D_zz_;
Grid_FFT<real_t,false> grad_x_, grad_y_, grad_z_;
std::vector<vec3_t<real_t>> vectk_;
@ -522,7 +522,7 @@ public:
aini_ ( 1.0/(1.0+the_config.get_value<double>("setup", "zstart")) ),
ngmapto_( the_config.get_value<size_t>("setup", "GridRes") ),
ngrid_( ngridself ), ngrid32_( std::pow(ngrid_, 1.5) ), mapratio_(real_t(ngrid_)/real_t(ngmapto_)),
XmL_ ( the_config.get_value<double>("cosmology", "Omega_L") / the_config.get_value<double>("cosmology", "Omega_m") ),
//XmL_ ( the_config.get_value<double>("cosmology", "Omega_L") / the_config.get_value<double>("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}),

111
src/cosmology_parameters.cc Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
#include <cosmology_parameters.hh>
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

View file

@ -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<double>("cosmology", "Omega_m");
double Ob = the_config.get_value<double>("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<double>("cosmology", "Omega_m");
double Om = the_cosmo_calc->cosmo_param_["Omega_m"];
Omega[cosmo_species::dm] = Om;
Omega[cosmo_species::baryon] = 0.0;
}

View file

@ -28,6 +28,7 @@
#include <general.hh>
#include <ic_generator.hh>
#include <cosmology_parameters.hh>
#include <particle_plt.hh>
@ -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();

View file

@ -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<double>("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;