2020-08-21 17:03:33 +02:00
// This file is part of monofonIC (MUSIC2)
// A software package to generate ICs for cosmological simulations
// Copyright (C) 2020 by Oliver Hahn
2020-09-07 22:55:39 +02:00
//
2020-08-21 17:03:33 +02:00
// 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.
2020-09-07 22:55:39 +02:00
//
2020-08-21 17:03:33 +02:00
// 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.
2020-09-07 22:55:39 +02:00
//
2020-08-21 17:03:33 +02:00
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
2019-05-07 01:05:16 +02:00
# pragma once
2020-09-01 15:29:40 +02:00
# include <map>
# include <string>
2020-03-29 14:49:17 +02:00
# include <physical_constants.hh>
2019-05-07 01:05:16 +02:00
# include <config_file.hh>
2020-09-10 23:08:00 +02:00
# include <general.hh>
2019-05-07 01:05:16 +02:00
2020-03-29 14:49:17 +02:00
namespace cosmology
{
2020-09-10 23:08:00 +02:00
//! structure for cosmological parameters
2020-09-07 22:55:39 +02:00
class parameters
{
2020-09-10 23:08:00 +02:00
public :
using defaultmmap_t = std : : map < std : : string , std : : map < std : : string , real_t > > ;
2020-09-07 22:55:39 +02:00
private :
2020-09-07 23:20:28 +02:00
std : : map < std : : string , double > pmap_ ; //!< All parameters are stored here as key-value pairs
2020-09-01 15:29:40 +02:00
2020-09-10 23:08:00 +02:00
static defaultmmap_t default_pmaps_ ; //!< Holds pre-defined parameter sets, see src/cosmology_parameters.cc
2020-09-07 22:55:39 +02:00
public :
2020-09-10 05:02:43 +02:00
//! get routine for cosmological parameter key-value pairs
2020-09-07 22:55:39 +02:00
double get ( const std : : string & key ) const
{
auto it = pmap_ . find ( key ) ;
if ( it = = pmap_ . end ( ) )
{
auto errmsg = std : : string ( " Cosmological parameter \' " ) + key + std : : string ( " \' does not exist in internal list. " ) ;
music : : elog < < errmsg < < std : : endl ;
throw std : : runtime_error ( errmsg . c_str ( ) ) ;
}
return it - > second ;
2020-09-01 15:29:40 +02:00
}
2020-03-29 14:49:17 +02:00
2020-09-10 05:02:43 +02:00
//! set routine for cosmological parameter key-value pairs
2020-09-07 22:55:39 +02:00
void set ( const std : : string & key , const double value )
{
auto it = pmap_ . find ( key ) ;
if ( it ! = pmap_ . end ( ) )
{
pmap_ [ key ] = value ;
}
else
{
auto errmsg = std : : string ( " Cosmological parameter \' " ) + key + std : : string ( " \' does not exist in internal list. Needs to be defaulted before it can be set! " ) ;
music : : elog < < errmsg < < std : : endl ;
throw std : : runtime_error ( errmsg . c_str ( ) ) ;
}
}
2020-03-29 14:49:17 +02:00
2020-09-10 05:02:43 +02:00
//! shortcut get routine for cosmological parameter key-value pairs through bracket operator
2020-09-07 22:55:39 +02:00
inline double operator [ ] ( const std : : string & key ) const { return this - > get ( key ) ; }
2019-05-07 01:05:16 +02:00
2020-09-10 23:08:00 +02:00
//! default constructor does nothing
parameters ( ) { }
2020-03-29 14:49:17 +02:00
2020-09-10 05:02:43 +02:00
//! default copy constructor
2020-09-07 22:55:39 +02:00
parameters ( const parameters & ) = default ;
2020-03-29 14:49:17 +02:00
2020-09-10 05:02:43 +02:00
//! main constructor for explicit construction from input config file
2020-09-07 23:20:28 +02:00
explicit parameters ( config_file & cf )
2020-03-29 14:49:17 +02:00
{
2020-09-10 23:08:00 +02:00
music : : ilog < < " ------------------------------------------------------------------------------- " < < std : : endl ;
if ( cf . get_value_safe < std : : string > ( " cosmology " , " ParameterSet " , " none " ) = = std : : string ( " none " ) )
{
//-------------------------------------------------------------------------------
// read parameters from config file
//-------------------------------------------------------------------------------
2020-09-07 23:20:28 +02:00
2020-09-10 23:08:00 +02:00
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 ) ;
2020-09-07 22:55:39 +02:00
2020-09-10 23:08:00 +02:00
// 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
2020-09-07 22:55:39 +02:00
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 " ) ;
pmap_ [ " Omega_nu_massive " ] = sum_m_nu / ( 93.14 * h * h ) ; // Omega_nu_m = \sum_i m_i / (93.14 eV h^2)
// 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 )
/ phys_const : : rhocrit_h2_SI / ( this - > get ( " h " ) * this - > get ( " h " ) ) ;
// massless neutrinos
pmap_ [ " Omega_nu_massless " ] = this - > get ( " N_ur " ) * this - > get ( " Omega_gamma " ) * 7. / 8. * std : : pow ( 4. / 11. , 4. / 3. ) ;
// total relativistic
pmap_ [ " Omega_r " ] = this - > get ( " Omega_gamma " ) + this - > get ( " Omega_nu_massless " ) ;
2020-09-10 23:08:00 +02:00
// 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 " ) ;
2020-09-07 22:55:39 +02:00
if ( cf . get_value_safe < bool > ( " cosmology " , " ZeroRadiation " , false ) )
{
pmap_ [ " Omega_r " ] = 0.0 ;
}
pmap_ [ " f_b " ] = this - > get ( " Omega_b " ) / this - > get ( " Omega_m " ) ;
pmap_ [ " f_c " ] = 1.0 - this - > get ( " f_b " ) ; // this means we add massive neutrinos to CDM here
2020-09-01 15:29:40 +02:00
2020-04-02 12:45:24 +02:00
# if 1
2020-09-07 22:55:39 +02:00
// assume zero curvature, take difference from dark energy
pmap_ [ " Omega_DE " ] + = 1.0 - this - > get ( " Omega_m " ) - this - > get ( " Omega_DE " ) - this - > get ( " Omega_r " ) ;
// Omega_DE += 1.0 - Omega_m - Omega_DE - Omega_r;
pmap_ [ " Omega_k " ] = 0.0 ;
2020-04-02 12:45:24 +02:00
# else
2020-09-07 22:55:39 +02:00
// allow for curvature
Omega_k = 1.0 - Omega_m - Omega_DE - Omega_r ;
2020-04-02 12:45:24 +02:00
# endif
2020-09-07 22:55:39 +02:00
pmap_ [ " dplus " ] = 0.0 ;
pmap_ [ " pnorm " ] = 0.0 ;
pmap_ [ " sqrtpnorm " ] = 0.0 ;
pmap_ [ " vfact " ] = 0.0 ;
music : : ilog < < " Cosmological parameters are: " < < std : : endl ;
music : : ilog < < " h = " < < std : : setw ( 16 ) < < this - > get ( " h " ) ;
if ( this - > get ( " A_s " ) > 0.0 )
music : : ilog < < " A_s = " < < std : : setw ( 16 ) < < this - > get ( " A_s " ) ;
else
music : : ilog < < " sigma_8 = " < < std : : setw ( 16 ) < < this - > get ( " sigma_8 " ) ;
2020-09-10 05:02:43 +02:00
music : : ilog < < " n_s = " < < std : : setw ( 16 ) < < this - > get ( " n_s " ) < < std : : endl ;
2020-09-07 22:55:39 +02:00
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 ;
2020-09-10 23:08:00 +02:00
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 ;
2020-09-07 22:55:39 +02:00
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 ;
2020-09-07 23:20:28 +02:00
//music::ilog << " Omega_k = " << 1.0 - this->get("Omega_m") - this->get("Omega_r") - this->get("Omega_DE") << std::endl;
2020-09-07 22:55:39 +02:00
if ( this - > get ( " Omega_r " ) > 0.0 )
{
music : : wlog < < " Radiation enabled, using Omega_r= " < < this - > get ( " Omega_r " ) < < " internally for backscaling. " < < std : : endl ;
music : : wlog < < " Make sure your sim code supports this, otherwise set [cosmology] / ZeroRadiation=true. " < < std : : endl ;
}
2020-04-04 01:24:05 +02:00
}
2020-09-10 23:08:00 +02:00
//! 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 ;
}
}
2020-09-07 22:55:39 +02:00
} ;
2020-09-10 23:08:00 +02:00
void print_ParameterSets ( void ) ;
2020-03-29 14:49:17 +02:00
} // namespace cosmology