From 0877b8d1f12804b7c070653f36610b79bf6e3bad Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Tue, 14 Feb 2023 10:37:29 -0800 Subject: [PATCH] WIP transition to monofonic background cosmology module --- src/cosmology_calculator.hh | 452 ++++++++++++++++++++++++++++++++++++ src/cosmology_parameters.cc | 118 ++++++++++ src/cosmology_parameters.hh | 225 ++++++++++++++++++ src/physical_constants.hh | 75 ++++++ 4 files changed, 870 insertions(+) create mode 100644 src/cosmology_calculator.hh create mode 100644 src/cosmology_parameters.cc create mode 100644 src/cosmology_parameters.hh create mode 100644 src/physical_constants.hh diff --git a/src/cosmology_calculator.hh b/src/cosmology_calculator.hh new file mode 100644 index 0000000..fb48eab --- /dev/null +++ b/src/cosmology_calculator.hh @@ -0,0 +1,452 @@ +// 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 . +#pragma once + +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +namespace cosmology_new +{ + +/*! + * @class cosmology::calculator + * @brief provides functions to compute cosmological quantities + * + * This class provides member functions to compute cosmological quantities + * related to the Friedmann equations and linear perturbation theory, it also + * provides the functionality to work with back-scaled cosmological fields + */ +class calculator +{ +public: + //! data structure to store cosmological parameters + cosmology::parameters cosmo_param_; + + //! pointer to an instance of a transfer function plugin + std::unique_ptr transfer_function_; + +private: + static constexpr double REL_PRECISION = 1e-10; + interpolated_function_1d D_of_a_, f_of_a_, a_of_D_; + double Dnow_, Dplus_start_, Dplus_target_, astart_, atarget_; + + double m_n_s_, m_sqrtpnorm_; + + //! wrapper for GSL adaptive integration routine, do not use if many integrations need to be done as it allocates and deallocates memory + //! set to 61-point Gauss-Kronrod and large workspace, used for sigma_8 normalisation + real_t integrate(double (*func)(double x, void *params), double a, double b, void *params) const + { + constexpr size_t wspace_size{100000}; + + double result{0.0}; + double error{0.0}; + + gsl_function F; + F.function = func; + F.params = params; + + auto errh = gsl_set_error_handler_off(); + gsl_integration_workspace *wspace = gsl_integration_workspace_alloc(wspace_size); + gsl_integration_qag(&F, a, b, 0, REL_PRECISION, wspace_size, GSL_INTEG_GAUSS61, wspace, &result, &error); + gsl_integration_workspace_free(wspace); + gsl_set_error_handler(errh); + + if (error / result > REL_PRECISION) + music::wlog << "no convergence in function 'integrate', rel. error=" << error / result << std::endl; + + return static_cast(result); + } + + //! compute the linear theory growth factor D+ by solving the single fluid ODE, returns tables D(a), f(a) + /*! + * @param tab_a reference to STL vector for values of a at which table entries exist + * @param tab_D reference to STL vector for values D(a) with a from tab_a + * @param tab_f reference to STL vector for values f(a)=dlog D / dlog a with a from tab_a + */ + void compute_growth( std::vector& tab_a, std::vector& tab_D, std::vector& tab_f ) + { + using v_t = vec_t<3, double>; + + // set ICs, very deep in radiation domination + const double a0 = 1e-10; + const double D0 = a0; + const double Dprime0 = 2.0 * D0 * H_of_a(a0) / std::pow(phys_const::c_SI, 2); + const double t0 = 1.0 / (a0 * H_of_a(a0)); + + v_t y0({a0, D0, Dprime0}); + + // set up integration + double dt = 1e-9; + double dtdid, dtnext; + const double amax = 2.0; + + v_t yy(y0); + double t = t0; + const double eps = 1e-10; + + const double Omega_m( cosmo_param_["Omega_m"] ), H0( cosmo_param_["H0"] ); + + while (yy[0] < amax) + { + // RHS of ODEs + auto rhs = [&](double t, v_t y) -> v_t { + auto a = y[0]; + auto D = y[1]; + auto Dprime = y[2]; + v_t dy; + // da/dtau = a^2 H(a) + dy[0] = a * a * H_of_a(a); + // d D/dtau + dy[1] = Dprime; + // d^2 D / dtau^2 + dy[2] = -a * H_of_a(a) * Dprime + 3.0 / 2.0 * Omega_m * std::pow(H0, 2) * D / a; + return dy; + }; + + // scale by predicted value to get approx. constant fractional errors + v_t yyscale = yy.abs() + dt * rhs(t, yy).abs(); + + // call integrator + ode_integrate::rk_step_qs(dt, t, yy, yyscale, rhs, eps, dtdid, dtnext); + + tab_a.push_back(yy[0]); + tab_D.push_back(yy[1]); + tab_f.push_back(yy[2]); // temporarily store D' in table + + dt = dtnext; + } + + // compute f, before we stored here D' + for (size_t i = 0; i < tab_a.size(); ++i) + { + tab_f[i] = tab_f[i] / (tab_a[i] * H_of_a(tab_a[i]) * tab_D[i]); + tab_D[i] = tab_D[i]; + tab_a[i] = tab_a[i]; + } + } + +public: + + calculator() = delete; + + calculator(const calculator& c) = delete; + + //! constructor for a cosmology calculator object + /*! + * @param acosmo a cosmological parameters structure + * @param pTransferFunction pointer to an instance of a transfer function object + */ + explicit calculator(config_file &cf) + : cosmo_param_(cf), astart_( 1.0/(1.0+cf.get_value("setup","zstart")) ), + atarget_( 1.0/(1.0+cf.get_value_safe("cosmology","ztarget",0.0)) ) + { + // pre-compute growth factors and store for interpolation + std::vector tab_a, tab_D, tab_f; + this->compute_growth(tab_a, tab_D, tab_f); + D_of_a_.set_data(tab_a,tab_D); + f_of_a_.set_data(tab_a,tab_f); + a_of_D_.set_data(tab_D,tab_a); + Dnow_ = D_of_a_(1.0); + + Dplus_start_ = D_of_a_( astart_ ) / Dnow_; + Dplus_target_ = D_of_a_( atarget_ ) / Dnow_; + + music::ilog << "Linear growth factors: D+_target = " << Dplus_target_ << ", D+_start = " << Dplus_start_ << std::endl; + + // set up transfer functions and compute normalisation + transfer_function_ = select_TransferFunction_plugin(cf, cosmo_param_); + transfer_function_->intialise(); + if( !transfer_function_->tf_isnormalised_ ){ + cosmo_param_.set("pnorm", this->compute_pnorm_from_sigma8() ); + }else{ + cosmo_param_.set("pnorm", 1.0/Dplus_target_/Dplus_target_); + auto sigma8 = this->compute_sigma8(); + music::ilog << "Measured sigma_8 for given PS normalisation is " << sigma8 << std::endl; + } + cosmo_param_.set("sqrtpnorm", std::sqrt(cosmo_param_["pnorm"])); + + music::ilog << std::setw(32) << std::left << "TF supports distinct CDM+baryons" + << " : " << (transfer_function_->tf_is_distinct() ? "yes" : "no") << std::endl; + music::ilog << std::setw(32) << std::left << "TF maximum wave number" + << " : " << transfer_function_->get_kmax() << " h/Mpc" << std::endl; + + m_n_s_ = cosmo_param_["n_s"]; + m_sqrtpnorm_ = cosmo_param_["sqrtpnorm"]; + } + + ~calculator() { } + + //! Write out a correctly scaled power spectrum at time a + void write_powerspectrum(real_t a, std::string fname) const + { + // const real_t Dplus0 = this->get_growth_factor(a); + + if (CONFIG::MPI_task_rank == 0) + { + double kmin = std::max(1e-4, transfer_function_->get_kmin()); + + // write power spectrum to a file + std::ofstream ofs(fname.c_str()); + std::stringstream ss; + ss << " ,ap=" << a << ""; + 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::endl; + for (double k = kmin; k < transfer_function_->get_kmax(); k *= 1.01) + { + 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; + } + + //! Write out a correctly scaled power spectrum at starting time + void write_transfer( std::string fname ) const + { + // const real_t Dplus0 = this->get_growth_factor(a); + + if (CONFIG::MPI_task_rank == 0) + { + double kmin = std::max(1e-4, transfer_function_->get_kmin()); + + // write power spectrum to a file + std::ofstream ofs(fname.c_str()); + std::stringstream ss; + ss << " ,ap=" << astart_ << ""; + ofs << "# " << std::setw(18) << "k [h/Mpc]" + << std::setw(20) << ("delta_c(k,a=ap)") + << std::setw(20) << ("delta_b(k,a=ap)") + << std::setw(20) << ("delta_m(k,a=ap)") + << std::setw(20) << ("delta_bc(k,a=ap)") + << std::endl; + double fb = cosmo_param_["f_b"], fc = cosmo_param_["f_c"]; + 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) << dc + << std::setw(20) << std::setprecision(10) << db + << std::setw(20) << std::setprecision(10) << dm + << std::setw(20) << std::setprecision(10) << dbc + 2 * tbc * (std::sqrt( Dplus_target_ / Dplus_start_ ) - 1.0) + << std::setw(20) << std::setprecision(10) << tc / std::pow( Dplus_start_ / Dplus_target_, 0.5 ) + << std::setw(20) << std::setprecision(10) << tb / std::pow( Dplus_start_ / Dplus_target_, 0.5 ) + << std::setw(20) << std::setprecision(10) << tm / std::pow( Dplus_start_ / Dplus_target_, 0.5 ) + << std::setw(20) << std::setprecision(10) << tbc / std::pow( Dplus_start_ / Dplus_target_, 0.5 ) + << std::endl; + } + } + music::ilog << "Wrote input transfer functions at a=" << astart_ << " to file \'" << fname << "\'" << std::endl; + } + + const cosmology::parameters &get_parameters(void) const noexcept + { + return cosmo_param_; + } + + //! return the value of the Hubble function H(a) = dloga/dt + inline double H_of_a(double a) const noexcept + { + double HH2 = 0.0; + HH2 += cosmo_param_["Omega_r"] / (a * a * a * a); + HH2 += cosmo_param_["Omega_m"] / (a * a * a); + HH2 += cosmo_param_["Omega_k"] / (a * a); + HH2 += cosmo_param_["Omega_DE"] * std::pow(a, -3. * (1. + cosmo_param_["w_0"] + cosmo_param_["w_a"])) * exp(-3. * (1.0 - a) * cosmo_param_["w_a"]); + return cosmo_param_["H0"] * std::sqrt(HH2); + } + + //! Computes the linear theory growth factor D+, normalised to D+(a=1)=1 + real_t get_growth_factor(real_t a) const noexcept + { + return D_of_a_(a) / Dnow_; + } + + //! Computes the inverse of get_growth_factor + real_t get_a( real_t Dplus ) const noexcept + { + return a_of_D_( Dplus * Dnow_ ); + } + + //! Computes the linear theory growth rate f + /*! Function computes (by interpolating on precalculated table) + * f = dlog D+ / dlog a + */ + real_t get_f(real_t a) const noexcept + { + return f_of_a_(a); + } + + //! Compute the factor relating particle displacement and velocity + /*! Function computes + * vfac = a * (H(a)/h) * dlogD+ / dlog a + */ + real_t get_vfact(real_t a) const noexcept + { + return f_of_a_(a) * a * H_of_a(a) / cosmo_param_["h"]; + } + + //! Integrand for the sigma_8 normalization of the power spectrum + /*! Returns the value of the primordial power spectrum multiplied with + the transfer function and the window function of 8 Mpc/h at wave number k */ + static double dSigma8(double k, void *pParams) + { + cosmology::calculator *pcc = reinterpret_cast(pParams); + + const double x = k * 8.0; + const double w = (x < 0.001)? 1.0-0.1*x*x : 3.0 * (std::sin(x) - x * std::cos(x)) / (x * x * x); + + static double nspect = (double)pcc->cosmo_param_["n_s"]; + double tf = pcc->transfer_function_->compute(k, delta_matter); + + //... no growth factor since we compute at z=0 and normalize so that D+(z=0)=1 + return k * k * w * w * pow((double)k, (double)nspect) * tf * tf; + } + + //! Integrand for the sigma_8 normalization of the power spectrum + /*! Returns the value of the primordial power spectrum multiplied with + the transfer function and the window function of 8 Mpc/h at wave number k */ + static double dSigma8_0(double k, void *pParams) + { + cosmology::calculator *pcc = reinterpret_cast(pParams); + + const double x = k * 8.0; + const double w = (x < 0.001)? 1.0-0.1*x*x : 3.0 * (std::sin(x) - x * std::cos(x)) / (x * x * x); + + static double nspect = static_cast(pcc->cosmo_param_["n_s"]); + double tf = pcc->transfer_function_->compute(k, delta_matter0); + + //... no growth factor since we compute at z=0 and normalize so that D+(z=0)=1 + return k * k * w * w * std::pow(k, nspect) * tf * tf; + } + + //! Computes the amplitude of a mode from the power spectrum + /*! Function evaluates the supplied transfer function transfer_function_ + * and returns the amplitude of fluctuations at wave number k (in h/Mpc) back-scaled to z=z_start + * @param k wave number at which to evaluate + * @param type one of the species: {delta,theta}_{matter,cdm,baryon,neutrino} + */ + inline real_t get_amplitude( const real_t k, const tf_type type) const + { + return std::pow(k, 0.5 * m_n_s_) * transfer_function_->compute(k, type) * m_sqrtpnorm_; + } + + //! Compute amplitude of the back-scaled delta_bc mode, with decaying velocity v_bc included or not (in which case delta_bc=const) + inline real_t get_amplitude_delta_bc( const real_t k, bool withvbc ) const + { + const real_t Dratio = Dplus_target_ / Dplus_start_; + const real_t dbc = transfer_function_->compute(k, delta_bc) + (withvbc? 2 * transfer_function_->compute(k, theta_bc) * (std::sqrt(Dratio) - 1.0) : 0.0); + // need to multiply with Dplus_target since sqrtpnorm rescales like that + return std::pow(k, 0.5 * m_n_s_) * dbc * (m_sqrtpnorm_ * Dplus_target_); + } + + //! Compute amplitude of the back-scaled relative velocity theta_bc mode if withvbc==true, otherwise return zero + inline real_t get_amplitude_theta_bc( const real_t k, bool withvbc ) const + { + const real_t Dratio = Dplus_target_ / Dplus_start_; + const real_t tbc = transfer_function_->compute(k, theta_bc) * std::sqrt(Dratio); + // need to multiply with Dplus_target since sqrtpnorm rescales like that + return withvbc ? std::pow(k, 0.5 * m_n_s_) * tbc * (m_sqrtpnorm_ * Dplus_target_) : 0.0; + } + + + //! Computes the normalization for the power spectrum + /*! + * integrates the power spectrum to fix the normalization to that given + * by the sigma_8 parameter + */ + real_t compute_sigma8(void) + { + real_t sigma0, kmin, kmax; + kmax = transfer_function_->get_kmax(); + kmin = transfer_function_->get_kmin(); + + if (!transfer_function_->tf_has_total0()) + sigma0 = 4.0 * M_PI * integrate(&dSigma8, static_cast(kmin), static_cast(kmax), this); + else{ + sigma0 = 4.0 * M_PI * integrate(&dSigma8_0, static_cast(kmin), static_cast(kmax), this); + } + + return std::sqrt(sigma0); + } + + //! Computes the normalization for the power spectrum + /*! + * integrates the power spectrum to fix the normalization to that given + * by the sigma_8 parameter + */ + real_t compute_pnorm_from_sigma8(void) + { + auto measured_sigma8 = this->compute_sigma8(); + return cosmo_param_["sigma_8"] * cosmo_param_["sigma_8"] / (measured_sigma8 * measured_sigma8); + } +}; + +//! compute the jeans sound speed +/*! given a density in g/cm^-3 and a mass in g it gives back the sound + * speed in cm/s for which the input mass is equal to the jeans mass + * @param rho density + * @param mass mass scale + * @returns jeans sound speed + */ +// inline double jeans_sound_speed(double rho, double mass) +// { +// const double G = 6.67e-8; +// return pow(6.0 * mass / M_PI * std::sqrt(rho) * std::pow(G, 1.5), 1.0 / 3.0); +// } + +} // namespace cosmology diff --git a/src/cosmology_parameters.cc b/src/cosmology_parameters.cc new file mode 100644 index 0000000..c3a9ba0 --- /dev/null +++ b/src/cosmology_parameters.cc @@ -0,0 +1,118 @@ +// 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 + +/** + * @brief namespace encapsulating all things cosmology + * + */ +namespace cosmology_new{ + +//! 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}}} +}; + +/** + * @brief Output all sets of cosmological parameters that we store internally + * + */ +void print_ParameterSets( void ){ + music::ilog << "Available cosmology parameter sets:" << std::endl; + parameters p; + p.print_available_sets(); + music::ilog << std::endl; +} + +}// end namespace cosmology \ No newline at end of file diff --git a/src/cosmology_parameters.hh b/src/cosmology_parameters.hh new file mode 100644 index 0000000..848b714 --- /dev/null +++ b/src/cosmology_parameters.hh @@ -0,0 +1,225 @@ +// 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 . +#pragma once + +#include +#include +#include +#include +#include +#include + +namespace cosmology_new +{ + //! 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 + { + 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; + } + + //! set routine for cosmological parameter key-value pairs + 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()); + } + } + + //! shortcut get routine for cosmological parameter key-value pairs through bracket operator + inline double operator[](const std::string &key) const { return this->get(key); } + + //! default constructor does nothing + parameters() {} + + //! default copy constructor + parameters(const parameters &) = default; + + //! main constructor for explicit construction from input config file + explicit parameters( config_file &cf ) + { + music::ilog << "-------------------------------------------------------------------------------" << std::endl; + + if( cf.get_value_safe("cosmology","ParameterSet","none") == std::string("none")) + { + //------------------------------------------------------------------------------- + // read parameters from config file + //------------------------------------------------------------------------------- + + 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"); + 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"); + + // 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)) + { + 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 + +#if 1 + // 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; +#else + // allow for curvature + Omega_k = 1.0 - Omega_m - Omega_DE - Omega_r; +#endif + + 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"); + 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 = " << 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) + { + 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; + } + } + + //! 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/src/physical_constants.hh b/src/physical_constants.hh new file mode 100644 index 0000000..ceec634 --- /dev/null +++ b/src/physical_constants.hh @@ -0,0 +1,75 @@ +// 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 . +#pragma once + +//------------------------------------------------------------------------------------ +// physical constants for convenience, all values have been taken from +// the 2018 edition of the Particle Data Group Booklet, +// http://pdg.lbl.gov/2019/mobile/reviews/pdf/rpp2018-rev-phys-constants-m.pdf + +namespace phys_const +{ +// helper value of pi so that we don't need to include any other header just for this +static constexpr double pi_ = 3.141592653589793115997963468544185161590576171875; + +//--- unit conversions --------------------------------------------------- + +// 1 Mpc in m +static constexpr double Mpc_SI = 3.0857e22; + +// 1 Gyr in s +static constexpr double Gyr_SI = 3.1536e16; + +// 1 eV in J +static constexpr double eV_SI = 1.602176487e-19; + +// 1 erg in J +static constexpr double erg_SI = 1e-7; + +//--- physical constants ------------------------------------------------ + +// speed of light c in m/s +static constexpr double c_SI = 2.99792458e8; + +// gravitational constant G in m^3/s^2/kg +static constexpr double G_SI = 6.6740800e-11; + +// Boltzmann constant k_B in kg m^2/s^2/K +static constexpr double kB_SI = 1.38064852e-23; + +// reduced Planck's quantum \hbar in kg m^2/s +static constexpr double hbar_SI = 1.054571800e-34; + +// Stefan-Boltzmann constant sigma in J/m^2/s/K^-4 +static constexpr double sigma_SI = (pi_ * pi_) * (kB_SI * kB_SI * kB_SI * kB_SI) / 60. / (hbar_SI * hbar_SI * hbar_SI) / (c_SI * c_SI); + +// electron mass in kg +static constexpr double me_SI = 9.10938356e-31; + +// proton mass in kg +static constexpr double mp_SI = 1.672621898e-27; + +// unified atomic mass unit (u) in kg +static constexpr double u_SI = 1.660539040e-27; + +// critical density of the Universe in h^2 kg/m^3 +static constexpr double rhocrit_h2_SI = 3 * 1e10 / (8 * pi_ * G_SI) / Mpc_SI / Mpc_SI; + +// solar mass in kg +static constexpr double Msun_SI = 1.98847e30; + +} // namespace phys_const \ No newline at end of file