mirror of
https://github.com/cosmo-sims/monofonIC.git
synced 2024-09-19 17:03:45 +02:00
add back decaying relative velocity mode at linear order only, cleanup of cosmology calculator module
This commit is contained in:
parent
f7a59c6634
commit
da7dc2b210
8 changed files with 145 additions and 129 deletions
|
@ -14,6 +14,7 @@ zstart = 24.0 # starting redshift
|
||||||
LPTorder = 3 # order of the LPT to be used (1,2 or 3)
|
LPTorder = 3 # order of the LPT to be used (1,2 or 3)
|
||||||
|
|
||||||
DoBaryons = no # also do baryon ICs?
|
DoBaryons = no # also do baryon ICs?
|
||||||
|
DoBaryonVrel = no # if doing baryons, incl. also relative velocity to linear order?
|
||||||
|
|
||||||
DoFixing = yes # do mode fixing à la Angulo&Pontzen (https://arxiv.org/abs/1603.05253)
|
DoFixing = yes # do mode fixing à la Angulo&Pontzen (https://arxiv.org/abs/1603.05253)
|
||||||
DoInversion = no # invert phases (for paired simulations)
|
DoInversion = no # invert phases (for paired simulations)
|
||||||
|
|
|
@ -38,7 +38,8 @@ namespace cosmology
|
||||||
* @brief provides functions to compute cosmological quantities
|
* @brief provides functions to compute cosmological quantities
|
||||||
*
|
*
|
||||||
* This class provides member functions to compute cosmological quantities
|
* This class provides member functions to compute cosmological quantities
|
||||||
* related to the Friedmann equations and linear perturbation theory
|
* related to the Friedmann equations and linear perturbation theory, it also
|
||||||
|
* provides the functionality to work with back-scaled cosmological fields
|
||||||
*/
|
*/
|
||||||
class calculator
|
class calculator
|
||||||
{
|
{
|
||||||
|
@ -54,34 +55,42 @@ private:
|
||||||
interpolated_function_1d<true,true,false> D_of_a_, f_of_a_, a_of_D_;
|
interpolated_function_1d<true,true,false> D_of_a_, f_of_a_, a_of_D_;
|
||||||
double Dnow_, Dplus_start_, Dplus_target_, astart_, atarget_;
|
double Dnow_, Dplus_start_, Dplus_target_, astart_, atarget_;
|
||||||
|
|
||||||
|
//! 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
|
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;
|
gsl_function F;
|
||||||
F.function = func;
|
F.function = func;
|
||||||
F.params = params;
|
F.params = params;
|
||||||
|
|
||||||
double result;
|
auto errh = gsl_set_error_handler_off();
|
||||||
double error;
|
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_set_error_handler_off();
|
gsl_integration_workspace_free(wspace);
|
||||||
gsl_integration_workspace *w = gsl_integration_workspace_alloc(100000);
|
gsl_set_error_handler(errh);
|
||||||
gsl_integration_qag(&F, a, b, 0, REL_PRECISION, 100000, 6, w, &result, &error);
|
|
||||||
|
|
||||||
gsl_integration_workspace_free(w);
|
|
||||||
|
|
||||||
gsl_set_error_handler(NULL);
|
|
||||||
|
|
||||||
if (error / result > REL_PRECISION)
|
if (error / result > REL_PRECISION)
|
||||||
music::wlog << "no convergence in function 'integrate', rel. error=" << error / result << std::endl;
|
music::wlog << "no convergence in function 'integrate', rel. error=" << error / result << std::endl;
|
||||||
|
|
||||||
return (real_t)result;
|
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 )
|
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>;
|
using v_t = vec_t<3, double>;
|
||||||
|
|
||||||
// set ICs
|
// set ICs, very deep in radiation domination
|
||||||
const double a0 = 1e-10;
|
const double a0 = 1e-10;
|
||||||
const double D0 = a0;
|
const double D0 = a0;
|
||||||
const double Dprime0 = 2.0 * D0 * H_of_a(a0) / std::pow(phys_const::c_SI, 2);
|
const double Dprime0 = 2.0 * D0 * H_of_a(a0) / std::pow(phys_const::c_SI, 2);
|
||||||
|
@ -125,7 +134,7 @@ private:
|
||||||
|
|
||||||
tab_a.push_back(yy[0]);
|
tab_a.push_back(yy[0]);
|
||||||
tab_D.push_back(yy[1]);
|
tab_D.push_back(yy[1]);
|
||||||
tab_f.push_back(yy[2]);
|
tab_f.push_back(yy[2]); // temporarily store D' in table
|
||||||
|
|
||||||
dt = dtnext;
|
dt = dtnext;
|
||||||
}
|
}
|
||||||
|
@ -162,7 +171,7 @@ public:
|
||||||
a_of_D_.set_data(tab_D,tab_a);
|
a_of_D_.set_data(tab_D,tab_a);
|
||||||
Dnow_ = D_of_a_(1.0);
|
Dnow_ = D_of_a_(1.0);
|
||||||
|
|
||||||
Dplus_start_ = D_of_a_( astart_ ) / Dnow_;
|
Dplus_start_ = D_of_a_( astart_ ) / Dnow_;
|
||||||
Dplus_target_ = D_of_a_( atarget_ ) / Dnow_;
|
Dplus_target_ = D_of_a_( atarget_ ) / Dnow_;
|
||||||
|
|
||||||
music::ilog << "Linear growth factors: D+_target = " << Dplus_target_ << ", D+_start = " << Dplus_start_ << std::endl;
|
music::ilog << "Linear growth factors: D+_target = " << Dplus_target_ << ", D+_start = " << Dplus_start_ << std::endl;
|
||||||
|
@ -185,9 +194,7 @@ public:
|
||||||
<< " : " << transfer_function_->get_kmax() << " h/Mpc" << std::endl;
|
<< " : " << transfer_function_->get_kmax() << " h/Mpc" << std::endl;
|
||||||
}
|
}
|
||||||
|
|
||||||
~calculator()
|
~calculator() { }
|
||||||
{
|
|
||||||
}
|
|
||||||
|
|
||||||
//! Write out a correctly scaled power spectrum at time a
|
//! Write out a correctly scaled power spectrum at time a
|
||||||
void write_powerspectrum(real_t a, std::string fname) const
|
void write_powerspectrum(real_t a, std::string fname) const
|
||||||
|
@ -213,22 +220,23 @@ public:
|
||||||
<< std::setw(20) << ("P_dbar(k,a=1)")
|
<< std::setw(20) << ("P_dbar(k,a=1)")
|
||||||
<< std::setw(20) << ("P_tcdm(k,a=1)")
|
<< std::setw(20) << ("P_tcdm(k,a=1)")
|
||||||
<< std::setw(20) << ("P_tbar(k,a=1)")
|
<< std::setw(20) << ("P_tbar(k,a=1)")
|
||||||
<< std::setw(20) << ("P_dtot(K,a=1)")
|
|
||||||
<< std::endl;
|
<< std::endl;
|
||||||
for (double k = kmin; k < transfer_function_->get_kmax(); k *= 1.01)
|
for (double k = kmin; k < transfer_function_->get_kmax(); k *= 1.01)
|
||||||
{
|
{
|
||||||
ofs << std::setw(20) << std::setprecision(10) << k
|
ofs << std::setw(20) << std::setprecision(10) << k
|
||||||
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, total)*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, cdm)*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, baryon)*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, vcdm)*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, vbaryon)*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, total0), 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, cdm0), 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, baryon0), 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, vcdm0), 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, vbaryon0), 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, vtotal), 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;
|
<< std::endl;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
@ -257,15 +265,24 @@ public:
|
||||||
double fb = cosmo_param_["f_b"], fc = cosmo_param_["f_c"];
|
double fb = cosmo_param_["f_b"], fc = cosmo_param_["f_c"];
|
||||||
for (double k = kmin; k < transfer_function_->get_kmax(); k *= 1.01)
|
for (double k = kmin; k < transfer_function_->get_kmax(); k *= 1.01)
|
||||||
{
|
{
|
||||||
double dm = this->get_amplitude(k, total) * Dplus_start_ / Dplus_target_;
|
const double dm = this->get_amplitude(k, delta_matter) * Dplus_start_ / Dplus_target_;
|
||||||
double dbc = this->get_amplitude(k, baryon) - this->get_amplitude(k, cdm);
|
const double dbc = this->get_amplitude(k, delta_bc);
|
||||||
double db = dm + fc * dbc;
|
const double db = dm + fc * dbc;
|
||||||
double dc = dm - fb * 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
|
ofs << std::setw(20) << std::setprecision(10) << k
|
||||||
<< std::setw(20) << std::setprecision(10) << dc
|
<< std::setw(20) << std::setprecision(10) << dc
|
||||||
<< std::setw(20) << std::setprecision(10) << db
|
<< std::setw(20) << std::setprecision(10) << db
|
||||||
<< std::setw(20) << std::setprecision(10) << dm
|
<< std::setw(20) << std::setprecision(10) << dm
|
||||||
<< std::setw(20) << std::setprecision(10) << dbc
|
<< 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;
|
<< std::endl;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
@ -327,9 +344,9 @@ public:
|
||||||
|
|
||||||
const double x = k * 8.0;
|
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);
|
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"];
|
static double nspect = (double)pcc->cosmo_param_["n_s"];
|
||||||
double tf = pcc->transfer_function_->compute(k, total);
|
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
|
//... 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;
|
return k * k * w * w * pow((double)k, (double)nspect) * tf * tf;
|
||||||
|
@ -346,47 +363,42 @@ public:
|
||||||
const double w = (x < 0.001)? 1.0-0.1*x*x : 3.0 * (std::sin(x) - x * std::cos(x)) / (x * x * x);
|
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"]);
|
static double nspect = static_cast<double>(pcc->cosmo_param_["n_s"]);
|
||||||
double tf = pcc->transfer_function_->compute(k, total0);
|
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
|
//... 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;
|
return k * k * w * w * std::pow(k, nspect) * tf * tf;
|
||||||
}
|
}
|
||||||
|
|
||||||
// static double dSigma_bc(double k, void *pParams)
|
|
||||||
// {
|
|
||||||
// if (k <= 0.0)
|
|
||||||
// return 0.0f;
|
|
||||||
|
|
||||||
// cosmology::calculator *pcc = reinterpret_cast<cosmology::calculator *>(pParams);
|
|
||||||
|
|
||||||
// static double nspect = static_cast<double>(pcc->cosmo_param_.nspect);
|
|
||||||
// double tf = pcc->transfer_function_->compute(k, deltabc);
|
|
||||||
|
|
||||||
// //... no growth factor since we compute at z=0 and normalize so that D+(z=0)=1
|
|
||||||
// return k * k * std::pow(k, nspect) * tf * tf; // * cosmo_param_.sqrtpnorm * Dplus_target_;
|
|
||||||
// }
|
|
||||||
|
|
||||||
//! Computes the amplitude of a mode from the power spectrum
|
//! Computes the amplitude of a mode from the power spectrum
|
||||||
/*! Function evaluates the supplied transfer function ptransfer_fun_
|
/*! Function evaluates the supplied transfer function transfer_function_
|
||||||
* and returns the amplitude of fluctuations at wave number k at z=0
|
* 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 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
|
inline real_t get_amplitude( const real_t k, const tf_type type) const
|
||||||
{
|
{
|
||||||
return std::pow(k, 0.5 * cosmo_param_["n_s"]) * transfer_function_->compute(k, type) * cosmo_param_["sqrtpnorm"];// * ((type!=deltabc)? 1.0 : 1.0/Dplus_target_);
|
return std::pow(k, 0.5 * cosmo_param_["n_s"]) * transfer_function_->compute(k, type) * cosmo_param_["sqrtpnorm"];
|
||||||
}
|
}
|
||||||
|
|
||||||
inline real_t get_amplitude_phibc( const real_t k ) const
|
//! 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
|
// need to multiply with Dplus_target since sqrtpnorm rescales like that
|
||||||
return -std::pow(k, 0.5 * cosmo_param_["n_s"]-2.0) * transfer_function_->compute(k, deltabc) * cosmo_param_["sqrtpnorm"] * Dplus_target_;
|
return std::pow(k, 0.5 * cosmo_param_["n_s"]) * dbc * (cosmo_param_["sqrtpnorm"] * Dplus_target_);
|
||||||
}
|
}
|
||||||
inline real_t get_amplitude_rhobc( const real_t k ) const
|
|
||||||
|
//! 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
|
// need to multiply with Dplus_target since sqrtpnorm rescales like that
|
||||||
return std::pow(k, 0.5 * cosmo_param_["n_s"]) * transfer_function_->compute(k, deltabc) * cosmo_param_["sqrtpnorm"] * Dplus_target_;
|
return withvbc ? std::pow(k, 0.5 * cosmo_param_["n_s"]) * tbc * (cosmo_param_["sqrtpnorm"] * Dplus_target_) : 0.0;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
//! Computes the normalization for the power spectrum
|
//! Computes the normalization for the power spectrum
|
||||||
/*!
|
/*!
|
||||||
* integrates the power spectrum to fix the normalization to that given
|
* integrates the power spectrum to fix the normalization to that given
|
||||||
|
@ -407,19 +419,6 @@ public:
|
||||||
return std::sqrt(sigma0);
|
return std::sqrt(sigma0);
|
||||||
}
|
}
|
||||||
|
|
||||||
// real_t compute_sigma_bc( void )
|
|
||||||
// {
|
|
||||||
// real_t sigma0, kmin, kmax;
|
|
||||||
// kmax = 100.0; //transfer_function_->get_kmax();
|
|
||||||
// kmin = transfer_function_->get_kmin();
|
|
||||||
|
|
||||||
// sigma0 = 4.0 * M_PI * integrate(&dSigma_bc, static_cast<double>(kmin), static_cast<double>(kmax), this);
|
|
||||||
// sigma0 = std::sqrt(sigma0);
|
|
||||||
// sigma0 *= cosmo_param_.sqrtpnorm * Dplus_target_;
|
|
||||||
// std::cout << "kmin = " << kmin << ", kmax = " << kmax << ", sigma_bc = " << sigma0 << std::endl;
|
|
||||||
// return sigma0;
|
|
||||||
// }
|
|
||||||
|
|
||||||
//! Computes the normalization for the power spectrum
|
//! Computes the normalization for the power spectrum
|
||||||
/*!
|
/*!
|
||||||
* integrates the power spectrum to fix the normalization to that given
|
* integrates the power spectrum to fix the normalization to that given
|
||||||
|
|
|
@ -31,7 +31,7 @@ namespace cosmology
|
||||||
std::map<std::string, double> pmap_; //!< All parameters are stored here as key-value pairs
|
std::map<std::string, double> pmap_; //!< All parameters are stored here as key-value pairs
|
||||||
|
|
||||||
public:
|
public:
|
||||||
//!< get routine for cosmological parameter key-value pairs
|
//! get routine for cosmological parameter key-value pairs
|
||||||
double get(const std::string &key) const
|
double get(const std::string &key) const
|
||||||
{
|
{
|
||||||
auto it = pmap_.find(key);
|
auto it = pmap_.find(key);
|
||||||
|
@ -44,7 +44,7 @@ namespace cosmology
|
||||||
return it->second;
|
return it->second;
|
||||||
}
|
}
|
||||||
|
|
||||||
//!< set routine for cosmological parameter key-value pairs
|
//! set routine for cosmological parameter key-value pairs
|
||||||
void set(const std::string &key, const double value)
|
void set(const std::string &key, const double value)
|
||||||
{
|
{
|
||||||
auto it = pmap_.find(key);
|
auto it = pmap_.find(key);
|
||||||
|
@ -60,16 +60,16 @@ namespace cosmology
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
//!< shortcut get routine for cosmological parameter key-value pairs through bracket operator
|
//! shortcut get routine for cosmological parameter key-value pairs through bracket operator
|
||||||
inline double operator[](const std::string &key) const { return this->get(key); }
|
inline double operator[](const std::string &key) const { return this->get(key); }
|
||||||
|
|
||||||
//!< no default constructor
|
//! no default constructor
|
||||||
parameters() = delete;
|
parameters() = delete;
|
||||||
|
|
||||||
//!< default copy constructor
|
//! default copy constructor
|
||||||
parameters(const parameters &) = default;
|
parameters(const parameters &) = default;
|
||||||
|
|
||||||
//!< main constructor for explicit construction from input config file
|
//! main constructor for explicit construction from input config file
|
||||||
explicit parameters( config_file &cf )
|
explicit parameters( config_file &cf )
|
||||||
{
|
{
|
||||||
// CMB
|
// CMB
|
||||||
|
@ -154,7 +154,7 @@ namespace cosmology
|
||||||
music::ilog << "A_s = " << std::setw(16) << this->get("A_s");
|
music::ilog << "A_s = " << std::setw(16) << this->get("A_s");
|
||||||
else
|
else
|
||||||
music::ilog << "sigma_8 = " << std::setw(16) << this->get("sigma_8");
|
music::ilog << "sigma_8 = " << std::setw(16) << this->get("sigma_8");
|
||||||
music::ilog << "nspec = " << std::setw(16) << this->get("n_s") << std::endl;
|
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_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 = " << std::setw(11) << 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_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;
|
||||||
|
|
|
@ -25,19 +25,20 @@
|
||||||
|
|
||||||
enum tf_type
|
enum tf_type
|
||||||
{
|
{
|
||||||
total,
|
delta_matter,
|
||||||
cdm,
|
delta_cdm,
|
||||||
baryon,
|
delta_baryon,
|
||||||
vtotal,
|
theta_matter,
|
||||||
vcdm,
|
theta_cdm,
|
||||||
vbaryon,
|
theta_baryon,
|
||||||
deltabc,
|
delta_bc,
|
||||||
total0,
|
theta_bc,
|
||||||
cdm0,
|
delta_matter0,
|
||||||
baryon0,
|
delta_cdm0,
|
||||||
vtotal0,
|
delta_baryon0,
|
||||||
vcdm0,
|
theta_matter0,
|
||||||
vbaryon0,
|
theta_cdm0,
|
||||||
|
theta_baryon0,
|
||||||
};
|
};
|
||||||
|
|
||||||
class TransferFunction_plugin
|
class TransferFunction_plugin
|
||||||
|
|
|
@ -45,7 +45,7 @@ int initialise( config_file& the_config )
|
||||||
the_random_number_generator = std::move(select_RNG_plugin(the_config));
|
the_random_number_generator = std::move(select_RNG_plugin(the_config));
|
||||||
the_cosmo_calc = std::make_unique<cosmology::calculator>(the_config);
|
the_cosmo_calc = std::make_unique<cosmology::calculator>(the_config);
|
||||||
the_output_plugin = std::move(select_output_plugin(the_config, the_cosmo_calc));
|
the_output_plugin = std::move(select_output_plugin(the_config, the_cosmo_calc));
|
||||||
|
|
||||||
return 0;
|
return 0;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -96,6 +96,9 @@ int run( config_file& the_config )
|
||||||
//--------------------------------------------------------------------------------------------------------
|
//--------------------------------------------------------------------------------------------------------
|
||||||
//! do baryon ICs?
|
//! do baryon ICs?
|
||||||
const bool bDoBaryons = the_config.get_value_safe<bool>("setup", "DoBaryons", false );
|
const bool bDoBaryons = the_config.get_value_safe<bool>("setup", "DoBaryons", false );
|
||||||
|
//! enable also back-scaled decaying relative velocity mode? only first order!
|
||||||
|
const bool bDoLinearBCcorr = the_config.get_value_safe<bool>("cosmology", "DoBaryonVrel", false);
|
||||||
|
// compute mass fractions
|
||||||
std::map< cosmo_species, double > Omega;
|
std::map< cosmo_species, double > Omega;
|
||||||
if( bDoBaryons ){
|
if( bDoBaryons ){
|
||||||
double Om = the_config.get_value<double>("cosmology", "Omega_m");
|
double Om = the_config.get_value<double>("cosmology", "Omega_m");
|
||||||
|
@ -164,8 +167,6 @@ int run( config_file& the_config )
|
||||||
const double g1 = -Dplus0;
|
const double g1 = -Dplus0;
|
||||||
const double g2 = ((LPTorder>1)? -3.0/7.0*Dplus0*Dplus0 : 0.0);
|
const double g2 = ((LPTorder>1)? -3.0/7.0*Dplus0*Dplus0 : 0.0);
|
||||||
const double g3 = ((LPTorder>2)? 1.0/3.0*Dplus0*Dplus0*Dplus0 : 0.0);
|
const double g3 = ((LPTorder>2)? 1.0/3.0*Dplus0*Dplus0*Dplus0 : 0.0);
|
||||||
// const double g3a = ((LPTorder>2)? 1.0/3.0*Dplus0*Dplus0*Dplus0 : 0.0);
|
|
||||||
// const double g3b = ((LPTorder>2)? -10.0/21.*Dplus0*Dplus0*Dplus0 : 0.0);
|
|
||||||
const double g3c = ((LPTorder>2)? 1.0/7.0*Dplus0*Dplus0*Dplus0 : 0.0);
|
const double g3c = ((LPTorder>2)? 1.0/7.0*Dplus0*Dplus0*Dplus0 : 0.0);
|
||||||
|
|
||||||
// vfac = d log D+ / dt
|
// vfac = d log D+ / dt
|
||||||
|
@ -328,7 +329,7 @@ int run( config_file& the_config )
|
||||||
phi.FourierTransformForward(false);
|
phi.FourierTransformForward(false);
|
||||||
phi.assign_function_of_grids_kdep([&](auto k, auto wn) {
|
phi.assign_function_of_grids_kdep([&](auto k, auto wn) {
|
||||||
real_t kmod = k.norm();
|
real_t kmod = k.norm();
|
||||||
ccomplex_t delta = wn * the_cosmo_calc->get_amplitude(kmod, total);
|
ccomplex_t delta = wn * the_cosmo_calc->get_amplitude(kmod, delta_matter);
|
||||||
return -delta / (kmod * kmod);
|
return -delta / (kmod * kmod);
|
||||||
}, wnoise);
|
}, wnoise);
|
||||||
|
|
||||||
|
@ -525,7 +526,7 @@ int run( config_file& the_config )
|
||||||
wnoise.FourierTransformForward();
|
wnoise.FourierTransformForward();
|
||||||
rho.FourierTransformForward(false);
|
rho.FourierTransformForward(false);
|
||||||
rho.assign_function_of_grids_kdep( [&]( auto k, auto wn ){
|
rho.assign_function_of_grids_kdep( [&]( auto k, auto wn ){
|
||||||
return wn * the_cosmo_calc->get_amplitude_rhobc(k.norm());;
|
return wn * the_cosmo_calc->get_amplitude_delta_bc(k.norm(),bDoLinearBCcorr);
|
||||||
}, wnoise );
|
}, wnoise );
|
||||||
rho.zero_DC_mode();
|
rho.zero_DC_mode();
|
||||||
rho.FourierTransformBackward();
|
rho.FourierTransformBackward();
|
||||||
|
@ -556,7 +557,7 @@ int run( config_file& the_config )
|
||||||
wnoise.FourierTransformForward();
|
wnoise.FourierTransformForward();
|
||||||
rho.FourierTransformForward(false);
|
rho.FourierTransformForward(false);
|
||||||
rho.assign_function_of_grids_kdep( [&]( auto k, auto wn ){
|
rho.assign_function_of_grids_kdep( [&]( auto k, auto wn ){
|
||||||
return wn * the_cosmo_calc->get_amplitude_rhobc(k.norm());;
|
return wn * the_cosmo_calc->get_amplitude_delta_bc(k.norm(), false);
|
||||||
}, wnoise );
|
}, wnoise );
|
||||||
rho.zero_DC_mode();
|
rho.zero_DC_mode();
|
||||||
rho.FourierTransformBackward();
|
rho.FourierTransformBackward();
|
||||||
|
@ -682,6 +683,7 @@ int run( config_file& the_config )
|
||||||
A3[1]->FourierTransformForward();
|
A3[1]->FourierTransformForward();
|
||||||
A3[2]->FourierTransformForward();
|
A3[2]->FourierTransformForward();
|
||||||
}
|
}
|
||||||
|
wnoise.FourierTransformForward();
|
||||||
|
|
||||||
// write out positions
|
// write out positions
|
||||||
for( int idim=0; idim<3; ++idim ){
|
for( int idim=0; idim<3; ++idim ){
|
||||||
|
@ -769,6 +771,12 @@ int run( config_file& the_config )
|
||||||
tmp.kelem(idx) += vfac3 * (lg.gradient(idimp,tmp.get_k3(i,j,k)) * A3[idimpp]->kelem(idx) - lg.gradient(idimpp,tmp.get_k3(i,j,k)) * A3[idimp]->kelem(idx));
|
tmp.kelem(idx) += vfac3 * (lg.gradient(idimp,tmp.get_k3(i,j,k)) * A3[idimpp]->kelem(idx) - lg.gradient(idimpp,tmp.get_k3(i,j,k)) * A3[idimp]->kelem(idx));
|
||||||
}
|
}
|
||||||
|
|
||||||
|
// if multi-species, then add vbc component backwards
|
||||||
|
if( bDoBaryons & bDoLinearBCcorr ){
|
||||||
|
real_t knorm = wnoise.get_k<real_t>(i,j,k).norm();
|
||||||
|
tmp.kelem(idx) -= vfac1 * C_species * the_cosmo_calc->get_amplitude_theta_bc(knorm, bDoLinearBCcorr) * wnoise.kelem(i,j,k) * lg.gradient(idim,tmp.get_k3(i,j,k)) / (knorm*knorm);
|
||||||
|
}
|
||||||
|
|
||||||
// correct with interpolation kernel if we used interpolation to read out the positions (for glasses)
|
// correct with interpolation kernel if we used interpolation to read out the positions (for glasses)
|
||||||
if( the_output_plugin->write_species_as( this_species ) == output_type::particles && lattice_type == particle::lattice_glass){
|
if( the_output_plugin->write_species_as( this_species ) == output_type::particles && lattice_type == particle::lattice_glass){
|
||||||
tmp.kelem(idx) *= interp.compensation_kernel( tmp.get_k<real_t>(i,j,k) );
|
tmp.kelem(idx) *= interp.compensation_kernel( tmp.get_k<real_t>(i,j,k) );
|
||||||
|
|
|
@ -44,7 +44,7 @@ private:
|
||||||
{
|
{
|
||||||
music::ilog << "Reading tabulated transfer function data from file:" << std::endl
|
music::ilog << "Reading tabulated transfer function data from file:" << std::endl
|
||||||
<< " \'" << filename << "\'" << std::endl;
|
<< " \'" << filename << "\'" << std::endl;
|
||||||
|
|
||||||
std::string line;
|
std::string line;
|
||||||
std::ifstream ifs(filename.c_str());
|
std::ifstream ifs(filename.c_str());
|
||||||
|
|
||||||
|
@ -183,32 +183,35 @@ public:
|
||||||
|
|
||||||
switch (type)
|
switch (type)
|
||||||
{
|
{
|
||||||
case total0:
|
case delta_matter0:
|
||||||
case total:
|
case delta_matter:
|
||||||
return delta_m_(k);
|
return delta_m_(k);
|
||||||
|
|
||||||
case cdm0:
|
case delta_cdm0:
|
||||||
case cdm:
|
case delta_cdm:
|
||||||
return delta_c_(k);
|
return delta_c_(k);
|
||||||
|
|
||||||
case baryon0:
|
case delta_baryon0:
|
||||||
case baryon:
|
case delta_baryon:
|
||||||
return delta_b_(k);
|
return delta_b_(k);
|
||||||
|
|
||||||
case vtotal0:
|
case theta_matter0:
|
||||||
case vtotal:
|
case theta_matter:
|
||||||
return theta_m_(k);
|
return theta_m_(k);
|
||||||
|
|
||||||
case vcdm0:
|
case theta_cdm0:
|
||||||
case vcdm:
|
case theta_cdm:
|
||||||
return theta_c_(k);
|
return theta_c_(k);
|
||||||
|
|
||||||
case vbaryon0:
|
case theta_baryon0:
|
||||||
case vbaryon:
|
case theta_baryon:
|
||||||
return theta_b_(k);
|
return theta_b_(k);
|
||||||
|
|
||||||
case deltabc:
|
case delta_bc:
|
||||||
return delta_b_(k)-delta_c_(k);
|
return delta_b_(k)-delta_c_(k);
|
||||||
|
|
||||||
|
case theta_bc:
|
||||||
|
return theta_b_(k)-theta_c_(k);
|
||||||
|
|
||||||
default:
|
default:
|
||||||
throw std::runtime_error("Invalid type requested in transfer function evaluation");
|
throw std::runtime_error("Invalid type requested in transfer function evaluation");
|
||||||
|
|
|
@ -178,8 +178,6 @@ private:
|
||||||
the_ClassEngine_->getTk(z, k, dc, db, dn, dm, tc, tb, tn, tm);
|
the_ClassEngine_->getTk(z, k, dc, db, dn, dm, tc, tb, tn, tm);
|
||||||
|
|
||||||
const double h = cosmo_params_.get("h");
|
const double h = cosmo_params_.get("h");
|
||||||
const double fc = cosmo_params_.get("f_c");
|
|
||||||
const double fb = cosmo_params_.get("f_b");
|
|
||||||
|
|
||||||
for (size_t i = 0; i < k.size(); ++i)
|
for (size_t i = 0; i < k.size(); ++i)
|
||||||
{
|
{
|
||||||
|
@ -288,33 +286,35 @@ public:
|
||||||
switch (type)
|
switch (type)
|
||||||
{
|
{
|
||||||
// values at ztarget:
|
// values at ztarget:
|
||||||
case total:
|
case delta_matter:
|
||||||
val = delta_m_(k); break;
|
val = delta_m_(k); break;
|
||||||
case cdm:
|
case delta_cdm:
|
||||||
val = delta_c_(k); break;
|
val = delta_c_(k); break;
|
||||||
case baryon:
|
case delta_baryon:
|
||||||
val = delta_b_(k); break;
|
val = delta_b_(k); break;
|
||||||
case vtotal:
|
case theta_matter:
|
||||||
val = theta_m_(k); break;
|
val = theta_m_(k); break;
|
||||||
case vcdm:
|
case theta_cdm:
|
||||||
val = theta_c_(k); break;
|
val = theta_c_(k); break;
|
||||||
case vbaryon:
|
case theta_baryon:
|
||||||
val = theta_b_(k); break;
|
val = theta_b_(k); break;
|
||||||
case deltabc:
|
case delta_bc:
|
||||||
val = delta_b_(k)-delta_c_(k); break;
|
val = delta_b_(k)-delta_c_(k); break;
|
||||||
|
case theta_bc:
|
||||||
|
val = theta_b_(k)-theta_c_(k); break;
|
||||||
|
|
||||||
// values at zstart:
|
// values at zstart:
|
||||||
case total0:
|
case delta_matter0:
|
||||||
val = delta_m0_(k); break;
|
val = delta_m0_(k); break;
|
||||||
case cdm0:
|
case delta_cdm0:
|
||||||
val = delta_c0_(k); break;
|
val = delta_c0_(k); break;
|
||||||
case baryon0:
|
case delta_baryon0:
|
||||||
val = delta_b0_(k); break;
|
val = delta_b0_(k); break;
|
||||||
case vtotal0:
|
case theta_matter0:
|
||||||
val = theta_m0_(k); break;
|
val = theta_m0_(k); break;
|
||||||
case vcdm0:
|
case theta_cdm0:
|
||||||
val = theta_c0_(k); break;
|
val = theta_c0_(k); break;
|
||||||
case vbaryon0:
|
case theta_baryon0:
|
||||||
val = theta_b0_(k); break;
|
val = theta_b0_(k); break;
|
||||||
default:
|
default:
|
||||||
throw std::runtime_error("Invalid type requested in transfer function evaluation");
|
throw std::runtime_error("Invalid type requested in transfer function evaluation");
|
||||||
|
|
|
@ -237,7 +237,8 @@ public:
|
||||||
//! Computes the transfer function for k in Mpc/h by calling TFfit_onek
|
//! Computes the transfer function for k in Mpc/h by calling TFfit_onek
|
||||||
inline double compute(double k, tf_type type) const
|
inline double compute(double k, tf_type type) const
|
||||||
{
|
{
|
||||||
return (type!=deltabc)? etf_.at_k(k) : 0.0;
|
if( type == theta_bc || type == delta_bc ) return 0.0;
|
||||||
|
return etf_.at_k(k);
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double get_kmin(void) const
|
inline double get_kmin(void) const
|
||||||
|
@ -332,6 +333,7 @@ public:
|
||||||
|
|
||||||
inline double compute(double k, tf_type type) const
|
inline double compute(double k, tf_type type) const
|
||||||
{
|
{
|
||||||
|
if( type == theta_bc || type == delta_bc ) return 0.0;
|
||||||
return etf_.at_k(k) * pow(1.0 + pow(m_WDMalpha * k, 2.0 * wdmnu_), -5.0 / wdmnu_);
|
return etf_.at_k(k) * pow(1.0 + pow(m_WDMalpha * k, 2.0 * wdmnu_), -5.0 / wdmnu_);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -384,6 +386,8 @@ public:
|
||||||
double kkfs2 = kkfs * kkfs;
|
double kkfs2 = kkfs * kkfs;
|
||||||
double kkd2 = (k / kd_) * (k / kd_);
|
double kkd2 = (k / kd_) * (k / kd_);
|
||||||
|
|
||||||
|
if( type == theta_bc || type == delta_bc ) return 0.0;
|
||||||
|
|
||||||
// in principle the Green et al. (2004) works only up to k/k_fs < 1
|
// in principle the Green et al. (2004) works only up to k/k_fs < 1
|
||||||
// the fit crosses zero at (k/k_fs)**2 = 3/2, we just zero it there...
|
// the fit crosses zero at (k/k_fs)**2 = 3/2, we just zero it there...
|
||||||
if (kkfs2 < 1.5)
|
if (kkfs2 < 1.5)
|
||||||
|
|
Loading…
Reference in a new issue