mirror of
https://github.com/cosmo-sims/MUSIC.git
synced 2024-09-19 17:03:46 +02:00
WIP transition to monofonic background cosmology module
This commit is contained in:
parent
0fa33b9ea9
commit
0877b8d1f1
4 changed files with 870 additions and 0 deletions
452
src/cosmology_calculator.hh
Normal file
452
src/cosmology_calculator.hh
Normal file
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
#pragma once
|
||||
|
||||
#include <array>
|
||||
#include <vec.hh>
|
||||
|
||||
#include <cosmology_parameters.hh>
|
||||
#include <physical_constants.hh>
|
||||
#include <transfer_function_plugin.hh>
|
||||
#include <math/ode_integrate.hh>
|
||||
#include <logger.hh>
|
||||
|
||||
#include <math/interpolate.hh>
|
||||
|
||||
#include <gsl/gsl_integration.h>
|
||||
#include <gsl/gsl_errno.h>
|
||||
|
||||
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<TransferFunction_plugin> transfer_function_;
|
||||
|
||||
private:
|
||||
static constexpr double REL_PRECISION = 1e-10;
|
||||
interpolated_function_1d<true,true,false> 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<real_t>(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<double>& tab_a, std::vector<double>& tab_D, std::vector<double>& 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<double>("setup","zstart")) ),
|
||||
atarget_( 1.0/(1.0+cf.get_value_safe<double>("cosmology","ztarget",0.0)) )
|
||||
{
|
||||
// pre-compute growth factors and store for interpolation
|
||||
std::vector<double> 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<cosmology::calculator *>(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<cosmology::calculator *>(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<double>(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<double>(kmin), static_cast<double>(kmax), this);
|
||||
else{
|
||||
sigma0 = 4.0 * M_PI * integrate(&dSigma8_0, static_cast<double>(kmin), static_cast<double>(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
|
118
src/cosmology_parameters.cc
Normal file
118
src/cosmology_parameters.cc
Normal file
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
#include <cosmology_parameters.hh>
|
||||
|
||||
/**
|
||||
* @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
|
225
src/cosmology_parameters.hh
Normal file
225
src/cosmology_parameters.hh
Normal file
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
#pragma once
|
||||
|
||||
#include <map>
|
||||
#include <string>
|
||||
#include <cmath>
|
||||
#include <physical_constants.hh>
|
||||
#include <config_file.hh>
|
||||
#include <general.hh>
|
||||
|
||||
namespace cosmology_new
|
||||
{
|
||||
//! 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
|
||||
{
|
||||
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<std::string>("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<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");
|
||||
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<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
|
||||
|
||||
#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
|
75
src/physical_constants.hh
Normal file
75
src/physical_constants.hh
Normal file
|
@ -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 <http://www.gnu.org/licenses/>.
|
||||
#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
|
Loading…
Reference in a new issue