mirror of
https://github.com/cosmo-sims/MUSIC.git
synced 2024-09-19 17:03:46 +02:00
WIP added CLASS transfer functions from monofonic, and CLASS as a project module
This commit is contained in:
parent
ef9afa5f7c
commit
526af6affa
11 changed files with 420 additions and 106 deletions
|
@ -179,6 +179,9 @@ elseif(CODE_PRECISION STREQUAL "LONGDOUBLE")
|
||||||
target_link_libraries(${PRGNAME} PRIVATE FFTW3::FFTW3_LONGDOUBLE_SERIAL)
|
target_link_libraries(${PRGNAME} PRIVATE FFTW3::FFTW3_LONGDOUBLE_SERIAL)
|
||||||
endif()
|
endif()
|
||||||
|
|
||||||
|
target_include_directories(${PRGNAME} PRIVATE ${GSL_INCLUDE_DIR})
|
||||||
|
target_link_libraries(${PRGNAME} PRIVATE GSL::gsl)
|
||||||
|
|
||||||
if(HDF5_FOUND)
|
if(HDF5_FOUND)
|
||||||
target_link_libraries(${PRGNAME} PRIVATE ${HDF5_LIBRARIES})
|
target_link_libraries(${PRGNAME} PRIVATE ${HDF5_LIBRARIES})
|
||||||
target_include_directories(${PRGNAME} PRIVATE ${HDF5_INCLUDE_DIRS})
|
target_include_directories(${PRGNAME} PRIVATE ${HDF5_INCLUDE_DIRS})
|
||||||
|
@ -192,10 +195,22 @@ if(TIRPC_FOUND)
|
||||||
target_compile_options(${PRGNAME} PRIVATE "-DHAVE_TIRPC")
|
target_compile_options(${PRGNAME} PRIVATE "-DHAVE_TIRPC")
|
||||||
endif(TIRPC_FOUND)
|
endif(TIRPC_FOUND)
|
||||||
|
|
||||||
|
########################################################################################################################
|
||||||
|
# include PANPHASIA
|
||||||
if(ENABLE_PANPHASIA)
|
if(ENABLE_PANPHASIA)
|
||||||
target_compile_options(${PRGNAME} PRIVATE "-DHAVE_PANPHASIA")
|
target_compile_options(${PRGNAME} PRIVATE "-DHAVE_PANPHASIA")
|
||||||
endif(ENABLE_PANPHASIA)
|
endif(ENABLE_PANPHASIA)
|
||||||
|
########################################################################################################################
|
||||||
|
|
||||||
target_include_directories(${PRGNAME} PRIVATE ${GSL_INCLUDE_DIR})
|
|
||||||
target_link_libraries(${PRGNAME} PRIVATE GSL::gsl)
|
########################################################################################################################
|
||||||
|
# include CLASS
|
||||||
|
option(ENABLE_CLASS "Enable CLASS support (as a submodule)." ON)
|
||||||
|
|
||||||
|
if(ENABLE_CLASS)
|
||||||
|
include(${CMAKE_CURRENT_SOURCE_DIR}/ext/class.cmake)
|
||||||
|
target_link_libraries(${PRGNAME} PRIVATE class::libclass_cpp)
|
||||||
|
target_compile_definitions(${PRGNAME} PRIVATE "USE_CLASS")
|
||||||
|
endif()
|
||||||
|
########################################################################################################################
|
||||||
|
|
||||||
|
|
17
ext/class.cmake
Normal file
17
ext/class.cmake
Normal file
|
@ -0,0 +1,17 @@
|
||||||
|
cmake_minimum_required(VERSION 3.11)
|
||||||
|
include(FetchContent)
|
||||||
|
FetchContent_Declare(
|
||||||
|
class
|
||||||
|
GIT_REPOSITORY https://github.com/ohahn/class_public.git
|
||||||
|
GIT_TAG master
|
||||||
|
GIT_SHALLOW YES
|
||||||
|
GIT_PROGRESS TRUE
|
||||||
|
USES_TERMINAL_DOWNLOAD TRUE # <---- this is needed only for Ninja
|
||||||
|
)
|
||||||
|
|
||||||
|
FetchContent_GetProperties(class)
|
||||||
|
if(NOT class_POPULATED)
|
||||||
|
set(FETCHCONTENT_QUIET OFF)
|
||||||
|
FetchContent_Populate(class)
|
||||||
|
add_subdirectory(${class_SOURCE_DIR} ${class_BINARY_DIR})
|
||||||
|
endif()
|
|
@ -80,70 +80,26 @@ inline T POW4( T a ){
|
||||||
//return a*a*a*a;
|
//return a*a*a*a;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
#if defined(_OPENMP)
|
||||||
|
#include <omp.h>
|
||||||
|
inline double get_wtime()
|
||||||
|
{
|
||||||
|
return omp_get_wtime();
|
||||||
|
}
|
||||||
|
#else
|
||||||
|
#include <ctime>
|
||||||
|
inline double get_wtime()
|
||||||
|
{
|
||||||
|
return std::clock() / double(CLOCKS_PER_SEC);
|
||||||
|
}
|
||||||
|
#endif
|
||||||
|
|
||||||
//! structure for cosmological parameters
|
// //! basic box/grid/refinement structure parameters
|
||||||
/*
|
// typedef struct {
|
||||||
typedef struct cosmology_old{
|
// unsigned levelmin, levelmax;
|
||||||
double
|
// double boxlength;
|
||||||
Omega_m, //!< baryon+dark matter density
|
// std::vector<unsigned> offx,offy,offz,llx,lly,llz;
|
||||||
Omega_b, //!< baryon matter density
|
// }Parameters;
|
||||||
Omega_DE, //!< dark energy density (cosmological constant or parameterised)
|
|
||||||
Omega_r, //!< photon + relativistic particle density
|
|
||||||
Omega_k, //!< curvature density
|
|
||||||
H0, //!< Hubble constant in km/s/Mpc
|
|
||||||
nspect, //!< long-wave spectral index (scale free is nspect=1)
|
|
||||||
sigma8, //!< power spectrum normalization
|
|
||||||
w_0, //!< dark energy equation of state parameter 1: w = w0 + a * wa
|
|
||||||
w_a, //!< dark energy equation of state parameter 2: w = w0 + a * wa
|
|
||||||
|
|
||||||
//Gamma, //!< shape parameter (of historical interest, as a free parameter)
|
|
||||||
//fnl, //!< non-gaussian contribution parameter
|
|
||||||
//w0, //!< dark energy equation of state parameter (not implemented, i.e. =1 at the moment)
|
|
||||||
//wa, //!< dark energy equation of state parameter (not implemented, i.e. =1 at the moment)
|
|
||||||
dplus, //!< linear perturbation growth factor
|
|
||||||
pnorm, //!< actual power spectrum normalisation factor
|
|
||||||
vfact, //!< velocity<->displacement conversion factor in Zel'dovich approx.
|
|
||||||
WDMmass, //!< Warm DM particle mass
|
|
||||||
WDMg_x, //!< Warm DM particle degrees of freedom
|
|
||||||
astart; //!< expansion factor a for which to generate initial conditions
|
|
||||||
|
|
||||||
cosmology( config_file cf )
|
|
||||||
{
|
|
||||||
double zstart = cf.get_value<double>( "setup", "zstart" );
|
|
||||||
|
|
||||||
astart = 1.0/(1.0+zstart);
|
|
||||||
Omega_b = cf.get_value<double>( "cosmology", "Omega_b" );
|
|
||||||
Omega_m = cf.get_value<double>( "cosmology", "Omega_m" );
|
|
||||||
Omega_DE = cf.get_value<double>( "cosmology", "Omega_L" );
|
|
||||||
w_0 = cf.get_value_safe<double>( "cosmology", "w0", -1.0 );
|
|
||||||
w_a = cf.get_value_safe<double>( "cosmology", "wa", 0.0 );
|
|
||||||
|
|
||||||
Omega_r = cf.get_value_safe<double>( "cosmology", "Omega_r", 0.0 ); // no longer default to nonzero (8.3e-5)
|
|
||||||
Omega_k = 1.0 - Omega_m - Omega_DE - Omega_r;
|
|
||||||
|
|
||||||
H0 = cf.get_value<double>( "cosmology", "H0" );
|
|
||||||
sigma8 = cf.get_value<double>( "cosmology", "sigma_8" );
|
|
||||||
nspect = cf.get_value<double>( "cosmology", "nspec" );
|
|
||||||
WDMg_x = cf.get_value_safe<double>( "cosmology", "WDMg_x", 1.5 );
|
|
||||||
WDMmass = cf.get_value_safe<double>( "cosmology", "WDMmass", 0.0 );
|
|
||||||
|
|
||||||
dplus = 0.0;
|
|
||||||
pnorm = 0.0;
|
|
||||||
vfact = 0.0;
|
|
||||||
}
|
|
||||||
|
|
||||||
cosmology( void )
|
|
||||||
{
|
|
||||||
|
|
||||||
}
|
|
||||||
}Cosmology;*/
|
|
||||||
|
|
||||||
//! basic box/grid/refinement structure parameters
|
|
||||||
typedef struct {
|
|
||||||
unsigned levelmin, levelmax;
|
|
||||||
double boxlength;
|
|
||||||
std::vector<unsigned> offx,offy,offz,llx,lly,llz;
|
|
||||||
}Parameters;
|
|
||||||
|
|
||||||
//! measure elapsed wallclock time
|
//! measure elapsed wallclock time
|
||||||
inline double time_seconds( void )
|
inline double time_seconds( void )
|
||||||
|
|
|
@ -261,16 +261,6 @@ void RNG_panphasia::backward_transform_field(real_t *field, int nx, int ny, int
|
||||||
FFTW_API(destroy_plan(ipf));
|
FFTW_API(destroy_plan(ipf));
|
||||||
}
|
}
|
||||||
|
|
||||||
#include <sys/time.h>
|
|
||||||
inline double get_wtime(void)
|
|
||||||
{
|
|
||||||
#ifdef _OPENMP
|
|
||||||
return omp_get_wtime();
|
|
||||||
#else
|
|
||||||
return (double)clock() / CLOCKS_PER_SEC;
|
|
||||||
#endif
|
|
||||||
}
|
|
||||||
|
|
||||||
void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
|
void RNG_panphasia::fill_grid(int level, DensityGrid<real_t> &R)
|
||||||
{
|
{
|
||||||
real_t *pr0, *pr1, *pr2, *pr3, *pr4;
|
real_t *pr0, *pr1, *pr2, *pr3, *pr4;
|
||||||
|
|
336
src/plugins/transfer_CLASS.cc
Normal file
336
src/plugins/transfer_CLASS.cc
Normal file
|
@ -0,0 +1,336 @@
|
||||||
|
// 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/>.
|
||||||
|
|
||||||
|
#ifdef USE_CLASS
|
||||||
|
|
||||||
|
#include <cmath>
|
||||||
|
#include <string>
|
||||||
|
#include <vector>
|
||||||
|
#include <memory>
|
||||||
|
#include <sstream>
|
||||||
|
|
||||||
|
#include <ClassEngine.hh>
|
||||||
|
|
||||||
|
#include <general.hh>
|
||||||
|
#include <config_file.hh>
|
||||||
|
#include <transfer_function.hh>
|
||||||
|
// #include <ic_generator.hh>
|
||||||
|
|
||||||
|
#include <math/interpolate.hh>
|
||||||
|
|
||||||
|
class transfer_CLASS_plugin : public transfer_function_plugin
|
||||||
|
{
|
||||||
|
private:
|
||||||
|
|
||||||
|
using transfer_function_plugin::cosmo_params_;
|
||||||
|
|
||||||
|
interpolated_function_1d<true, true, false> delta_c_, delta_b_, delta_n_, delta_m_, theta_c_, theta_b_, theta_n_, theta_m_;
|
||||||
|
interpolated_function_1d<true, true, false> delta_c0_, delta_b0_, delta_n0_, delta_m0_, theta_c0_, theta_b0_, theta_n0_, theta_m0_;
|
||||||
|
|
||||||
|
double zstart_, ztarget_, astart_, atarget_, kmax_, kmin_, h_, tnorm_;
|
||||||
|
|
||||||
|
ClassParams pars_;
|
||||||
|
std::unique_ptr<ClassEngine> the_ClassEngine_;
|
||||||
|
std::ofstream ofs_class_input_;
|
||||||
|
|
||||||
|
template <typename T>
|
||||||
|
void add_class_parameter(std::string parameter_name, const T parameter_value)
|
||||||
|
{
|
||||||
|
pars_.add(parameter_name, parameter_value);
|
||||||
|
ofs_class_input_ << parameter_name << " = " << parameter_value << std::endl;
|
||||||
|
}
|
||||||
|
|
||||||
|
//! Set up class parameters from MUSIC cosmological parameters
|
||||||
|
void init_ClassEngine(void)
|
||||||
|
{
|
||||||
|
//--- general parameters ------------------------------------------
|
||||||
|
add_class_parameter("z_max_pk", std::max(std::max(zstart_, ztarget_),199.0)); // use 1.2 as safety
|
||||||
|
add_class_parameter("P_k_max_h/Mpc", std::max(2.0,kmax_));
|
||||||
|
add_class_parameter("output", "dTk,vTk");
|
||||||
|
add_class_parameter("extra metric transfer functions","yes");
|
||||||
|
// add_class_parameter("lensing", "no");
|
||||||
|
|
||||||
|
//--- choose gauge ------------------------------------------------
|
||||||
|
// add_class_parameter("extra metric transfer functions", "yes");
|
||||||
|
add_class_parameter("gauge", "synchronous");
|
||||||
|
|
||||||
|
//--- cosmological parameters, densities --------------------------
|
||||||
|
add_class_parameter("h", cosmo_params_.get("h"));
|
||||||
|
|
||||||
|
add_class_parameter("Omega_b", cosmo_params_.get("Omega_b"));
|
||||||
|
add_class_parameter("Omega_cdm", cosmo_params_.get("Omega_c"));
|
||||||
|
add_class_parameter("Omega_k", cosmo_params_.get("Omega_k"));
|
||||||
|
add_class_parameter("Omega_fld", 0.0);
|
||||||
|
add_class_parameter("Omega_scf", 0.0);
|
||||||
|
|
||||||
|
|
||||||
|
// add_class_parameter("fluid_equation_of_state","CLP");
|
||||||
|
// add_class_parameter("w0_fld", -1 );
|
||||||
|
// add_class_parameter("wa_fld", 0. );
|
||||||
|
// add_class_parameter("cs2_fld", 1);
|
||||||
|
|
||||||
|
//--- massive neutrinos -------------------------------------------
|
||||||
|
#if 0
|
||||||
|
//default off
|
||||||
|
// add_class_parameter("Omega_ur",0.0);
|
||||||
|
add_class_parameter("N_ur", cosmo_params_.get("N_ur"));
|
||||||
|
add_class_parameter("N_ncdm", 0);
|
||||||
|
|
||||||
|
#else
|
||||||
|
|
||||||
|
add_class_parameter("N_ur", cosmo_params_.get("N_ur"));
|
||||||
|
add_class_parameter("N_ncdm", cosmo_params_.get("N_nu_massive"));
|
||||||
|
if( cosmo_params_.get("N_nu_massive") > 0 ){
|
||||||
|
std::stringstream sstr;
|
||||||
|
if( cosmo_params_.get("m_nu1") > 1e-9 ) sstr << cosmo_params_.get("m_nu1");
|
||||||
|
if( cosmo_params_.get("m_nu2") > 1e-9 ) sstr << ", " << cosmo_params_.get("m_nu2");
|
||||||
|
if( cosmo_params_.get("m_nu3") > 1e-9 ) sstr << ", " << cosmo_params_.get("m_nu3");
|
||||||
|
add_class_parameter("m_ncdm", sstr.str().c_str());
|
||||||
|
}
|
||||||
|
|
||||||
|
// change above to enable
|
||||||
|
//add_class_parameter("omega_ncdm", 0.0006451439);
|
||||||
|
//add_class_parameter("m_ncdm", "0.4");
|
||||||
|
//add_class_parameter("T_ncdm", 0.71611);
|
||||||
|
#endif
|
||||||
|
|
||||||
|
//--- cosmological parameters, primordial -------------------------
|
||||||
|
add_class_parameter("P_k_ini type", "analytic_Pk");
|
||||||
|
|
||||||
|
if( cosmo_params_.get("A_s") > 0.0 ){
|
||||||
|
add_class_parameter("A_s", cosmo_params_.get("A_s"));
|
||||||
|
}else{
|
||||||
|
add_class_parameter("sigma8", cosmo_params_.get("sigma_8"));
|
||||||
|
}
|
||||||
|
add_class_parameter("n_s", cosmo_params_.get("n_s"));
|
||||||
|
add_class_parameter("alpha_s", 0.0);
|
||||||
|
add_class_parameter("T_cmb", cosmo_params_.get("Tcmb"));
|
||||||
|
add_class_parameter("YHe", cosmo_params_.get("YHe"));
|
||||||
|
|
||||||
|
// additional parameters
|
||||||
|
add_class_parameter("reio_parametrization", "reio_none");
|
||||||
|
|
||||||
|
// precision parameters
|
||||||
|
add_class_parameter("k_per_decade_for_pk", 100);
|
||||||
|
add_class_parameter("k_per_decade_for_bao", 100);
|
||||||
|
add_class_parameter("compute damping scale", "yes");
|
||||||
|
add_class_parameter("tol_perturb_integration", 1.e-8);
|
||||||
|
add_class_parameter("tol_background_integration", 1e-9);
|
||||||
|
|
||||||
|
// high precision options from cl_permille.pre:
|
||||||
|
// precision file to be passed as input in order to achieve at least percent precision on scalar Cls
|
||||||
|
add_class_parameter("hyper_flat_approximation_nu", 7000.);
|
||||||
|
add_class_parameter("transfer_neglect_delta_k_S_t0", 0.17);
|
||||||
|
add_class_parameter("transfer_neglect_delta_k_S_t1", 0.05);
|
||||||
|
add_class_parameter("transfer_neglect_delta_k_S_t2", 0.17);
|
||||||
|
add_class_parameter("transfer_neglect_delta_k_S_e", 0.13);
|
||||||
|
add_class_parameter("delta_l_max", 1000);
|
||||||
|
|
||||||
|
int class_verbosity = 0;
|
||||||
|
|
||||||
|
add_class_parameter("background_verbose", class_verbosity);
|
||||||
|
add_class_parameter("thermodynamics_verbose", class_verbosity);
|
||||||
|
add_class_parameter("perturbations_verbose", class_verbosity);
|
||||||
|
add_class_parameter("transfer_verbose", class_verbosity);
|
||||||
|
add_class_parameter("primordial_verbose", class_verbosity);
|
||||||
|
add_class_parameter("spectra_verbose", class_verbosity);
|
||||||
|
add_class_parameter("nonlinear_verbose", class_verbosity);
|
||||||
|
add_class_parameter("lensing_verbose", class_verbosity);
|
||||||
|
add_class_parameter("output_verbose", class_verbosity);
|
||||||
|
|
||||||
|
// output parameters, only needed for the control CLASS .ini file that we output
|
||||||
|
std::stringstream zlist;
|
||||||
|
if (ztarget_ == zstart_)
|
||||||
|
zlist << ztarget_ << ((ztarget_!=0.0)? ", 0.0" : "");
|
||||||
|
else
|
||||||
|
zlist << std::max(ztarget_, zstart_) << ", " << std::min(ztarget_, zstart_) << ", 0.0";
|
||||||
|
add_class_parameter("z_pk", zlist.str());
|
||||||
|
|
||||||
|
music::ilog << "Computing transfer function via ClassEngine..." << std::endl;
|
||||||
|
double wtime = get_wtime();
|
||||||
|
|
||||||
|
the_ClassEngine_ = std::make_unique<ClassEngine>(pars_, false);
|
||||||
|
|
||||||
|
wtime = get_wtime() - wtime;
|
||||||
|
music::ilog << "CLASS took " << wtime << " s." << std::endl;
|
||||||
|
}
|
||||||
|
|
||||||
|
//! run ClassEngine with parameters set up
|
||||||
|
void run_ClassEngine(double z, std::vector<double> &k, std::vector<double> &dc, std::vector<double> &tc, std::vector<double> &db, std::vector<double> &tb,
|
||||||
|
std::vector<double> &dn, std::vector<double> &tn, std::vector<double> &dm, std::vector<double> &tm)
|
||||||
|
{
|
||||||
|
k.clear();
|
||||||
|
dc.clear(); db.clear(); dn.clear(); dm.clear();
|
||||||
|
tc.clear(); tb.clear(); tn.clear(); tm.clear();
|
||||||
|
|
||||||
|
the_ClassEngine_->getTk(z, k, dc, db, dn, dm, tc, tb, tn, tm);
|
||||||
|
|
||||||
|
const double h = cosmo_params_.get("h");
|
||||||
|
|
||||||
|
for (size_t i = 0; i < k.size(); ++i)
|
||||||
|
{
|
||||||
|
// convert to 'CAMB' format, since we interpolate loglog and
|
||||||
|
// don't want negative numbers...
|
||||||
|
auto ik2 = 1.0 / (k[i] * k[i]) * h * h;
|
||||||
|
dc[i] = -dc[i] * ik2;
|
||||||
|
db[i] = -db[i] * ik2;
|
||||||
|
dn[i] = -dn[i] * ik2;
|
||||||
|
dm[i] = -dm[i] * ik2;
|
||||||
|
tc[i] = -tc[i] * ik2;
|
||||||
|
tb[i] = -tb[i] * ik2;
|
||||||
|
tn[i] = -tn[i] * ik2;
|
||||||
|
tm[i] = -tm[i] * ik2;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
public:
|
||||||
|
explicit transfer_CLASS_plugin(config_file &cf, const cosmology::parameters& cosmo_params)
|
||||||
|
: transfer_function_plugin(cf,cosmo_params)
|
||||||
|
{
|
||||||
|
this->tf_isnormalised_ = true;
|
||||||
|
|
||||||
|
ofs_class_input_.open("input_class_parameters.ini", std::ios::trunc);
|
||||||
|
|
||||||
|
// all cosmological parameters need to be passed through the_cosmo_calc
|
||||||
|
|
||||||
|
ztarget_ = pcf_->get_value_safe<double>("cosmology", "ztarget", 0.0);
|
||||||
|
atarget_ = 1.0 / (1.0 + ztarget_);
|
||||||
|
zstart_ = pcf_->get_value<double>("setup", "zstart");
|
||||||
|
astart_ = 1.0 / (1.0 + zstart_);
|
||||||
|
|
||||||
|
h_ = cosmo_params_["h"];
|
||||||
|
|
||||||
|
if (cosmo_params_["A_s"] > 0.0) {
|
||||||
|
music::ilog << "CLASS: Using A_s=" << cosmo_params_["A_s"] << " to normalise the transfer function." << std::endl;
|
||||||
|
}else{
|
||||||
|
double sigma8 = cosmo_params_["sigma_8"];
|
||||||
|
if( sigma8 < 0 ){
|
||||||
|
throw std::runtime_error("Need to specify either A_s or sigma_8 for CLASS plugin...");
|
||||||
|
}
|
||||||
|
music::ilog << "CLASS: Using sigma8_ =" << sigma8<< " to normalise the transfer function." << std::endl;
|
||||||
|
}
|
||||||
|
|
||||||
|
// determine highest k we will need for the resolution selected
|
||||||
|
double lbox = pcf_->get_value<double>("setup", "BoxLength");
|
||||||
|
int nres = pcf_->get_value<double>("setup", "GridRes");
|
||||||
|
kmax_ = std::max(20.0, 2.0 * M_PI / lbox * nres / 2 * sqrt(3) * 2.0); // 120% of spatial diagonal, or k=10h Mpc-1
|
||||||
|
|
||||||
|
// initialise CLASS and get the normalisation
|
||||||
|
this->init_ClassEngine();
|
||||||
|
double A_s_ = the_ClassEngine_->get_A_s(); // this either the input one, or the one computed from sigma8
|
||||||
|
|
||||||
|
// compute the normalisation to interface with MUSIC
|
||||||
|
double k_p = cosmo_params["k_p"] / cosmo_params["h"];
|
||||||
|
tnorm_ = std::sqrt(2.0 * M_PI * M_PI * A_s_ * std::pow(1.0 / k_p, cosmo_params["n_s"] - 1) / std::pow(2.0 * M_PI, 3.0));
|
||||||
|
|
||||||
|
// compute the transfer function at z=0 using CLASS engine
|
||||||
|
std::vector<double> k, dc, tc, db, tb, dn, tn, dm, tm;
|
||||||
|
this->run_ClassEngine(0.0, k, dc, tc, db, tb, dn, tn, dm, tm);
|
||||||
|
|
||||||
|
delta_c0_.set_data(k, dc);
|
||||||
|
theta_c0_.set_data(k, tc);
|
||||||
|
delta_b0_.set_data(k, db);
|
||||||
|
theta_b0_.set_data(k, tb);
|
||||||
|
delta_n0_.set_data(k, dn);
|
||||||
|
theta_n0_.set_data(k, tn);
|
||||||
|
delta_m0_.set_data(k, dm);
|
||||||
|
theta_m0_.set_data(k, tm);
|
||||||
|
|
||||||
|
// compute the transfer function at z=z_target using CLASS engine
|
||||||
|
this->run_ClassEngine(ztarget_, k, dc, tc, db, tb, dn, tn, dm, tm);
|
||||||
|
delta_c_.set_data(k, dc);
|
||||||
|
theta_c_.set_data(k, tc);
|
||||||
|
delta_b_.set_data(k, db);
|
||||||
|
theta_b_.set_data(k, tb);
|
||||||
|
delta_n_.set_data(k, dn);
|
||||||
|
theta_n_.set_data(k, tn);
|
||||||
|
delta_m_.set_data(k, dm);
|
||||||
|
theta_m_.set_data(k, tm);
|
||||||
|
|
||||||
|
kmin_ = k[0];
|
||||||
|
kmax_ = k.back();
|
||||||
|
|
||||||
|
music::ilog << "CLASS table contains k = " << this->get_kmin() << " to " << this->get_kmax() << " h Mpc-1." << std::endl;
|
||||||
|
|
||||||
|
tf_distinct_ = true;
|
||||||
|
tf_withvel_ = true;
|
||||||
|
tf_withtotal0_ = true;
|
||||||
|
}
|
||||||
|
|
||||||
|
~transfer_CLASS_plugin()
|
||||||
|
{
|
||||||
|
}
|
||||||
|
|
||||||
|
inline double compute(double k, tf_type type) const
|
||||||
|
{
|
||||||
|
k *= h_;
|
||||||
|
|
||||||
|
if (k < kmin_ || k > kmax_)
|
||||||
|
{
|
||||||
|
return 0.0;
|
||||||
|
}
|
||||||
|
|
||||||
|
real_t val(0.0);
|
||||||
|
switch (type)
|
||||||
|
{
|
||||||
|
// values at ztarget:
|
||||||
|
case delta_matter:
|
||||||
|
val = delta_m_(k); break;
|
||||||
|
case delta_cdm:
|
||||||
|
val = delta_c_(k); break;
|
||||||
|
case delta_baryon:
|
||||||
|
val = delta_b_(k); break;
|
||||||
|
case theta_matter:
|
||||||
|
val = theta_m_(k); break;
|
||||||
|
case theta_cdm:
|
||||||
|
val = theta_c_(k); break;
|
||||||
|
case theta_baryon:
|
||||||
|
val = theta_b_(k); break;
|
||||||
|
case delta_bc:
|
||||||
|
val = delta_b_(k)-delta_c_(k); break;
|
||||||
|
case theta_bc:
|
||||||
|
val = theta_b_(k)-theta_c_(k); break;
|
||||||
|
|
||||||
|
// values at zstart:
|
||||||
|
case delta_matter0:
|
||||||
|
val = delta_m0_(k); break;
|
||||||
|
case delta_cdm0:
|
||||||
|
val = delta_c0_(k); break;
|
||||||
|
case delta_baryon0:
|
||||||
|
val = delta_b0_(k); break;
|
||||||
|
case theta_matter0:
|
||||||
|
val = theta_m0_(k); break;
|
||||||
|
case theta_cdm0:
|
||||||
|
val = theta_c0_(k); break;
|
||||||
|
case theta_baryon0:
|
||||||
|
val = theta_b0_(k); break;
|
||||||
|
default:
|
||||||
|
throw std::runtime_error("Invalid type requested in transfer function evaluation");
|
||||||
|
}
|
||||||
|
return val * tnorm_;
|
||||||
|
}
|
||||||
|
|
||||||
|
inline double get_kmin(void) const { return kmin_ / h_; }
|
||||||
|
inline double get_kmax(void) const { return kmax_ / h_; }
|
||||||
|
};
|
||||||
|
|
||||||
|
namespace
|
||||||
|
{
|
||||||
|
transfer_function_plugin_creator_concrete<transfer_CLASS_plugin> creator("CLASS");
|
||||||
|
}
|
||||||
|
|
||||||
|
#endif // USE_CLASS
|
|
@ -57,7 +57,7 @@ public:
|
||||||
}
|
}
|
||||||
|
|
||||||
//! computes the value of the BBKS transfer function for mode k (in h/Mpc)
|
//! computes the value of the BBKS transfer function for mode k (in h/Mpc)
|
||||||
inline double compute( double k, tf_type type ){
|
inline double compute( double k, tf_type type ) const{
|
||||||
double q, f1, f2;
|
double q, f1, f2;
|
||||||
|
|
||||||
if(k < 1e-7 )
|
if(k < 1e-7 )
|
||||||
|
@ -71,11 +71,11 @@ public:
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double get_kmin( void ){
|
inline double get_kmin( void ) const{
|
||||||
return 1e-4;
|
return 1e-4;
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double get_kmax( void ){
|
inline double get_kmax( void ) const{
|
||||||
return 1.e4;
|
return 1.e4;
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
|
|
|
@ -221,7 +221,7 @@ public:
|
||||||
}
|
}
|
||||||
|
|
||||||
// linear interpolation in log-log
|
// linear interpolation in log-log
|
||||||
inline double extrap_right(double k, const tf_type &type) {
|
inline double extrap_right(double k, const tf_type &type) const {
|
||||||
int n = m_tab_k.size() - 1, n1 = n - 1;
|
int n = m_tab_k.size() - 1, n1 = n - 1;
|
||||||
|
|
||||||
double v1(1.0), v2(1.0);
|
double v1(1.0), v2(1.0);
|
||||||
|
@ -267,7 +267,7 @@ public:
|
||||||
return 0.0;
|
return 0.0;
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double compute(double k, tf_type type) {
|
inline double compute(double k, tf_type type) const{
|
||||||
// use constant interpolation on the left side of the tabulated values
|
// use constant interpolation on the left side of the tabulated values
|
||||||
if (k < m_kmin) {
|
if (k < m_kmin) {
|
||||||
switch (type) {
|
switch (type) {
|
||||||
|
@ -320,9 +320,9 @@ public:
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double get_kmin(void) { return pow(10.0, m_tab_k[1]); }
|
inline double get_kmin(void) const { return pow(10.0, m_tab_k[1]); }
|
||||||
|
|
||||||
inline double get_kmax(void) {
|
inline double get_kmax(void) const {
|
||||||
return pow(10.0, m_tab_k[m_tab_k.size() - 2]);
|
return pow(10.0, m_tab_k[m_tab_k.size() - 2]);
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
|
|
|
@ -102,7 +102,7 @@ struct eisenstein_transfer
|
||||||
}
|
}
|
||||||
|
|
||||||
//! private member function: computes transfer function for mode k (k in Mpc)
|
//! private member function: computes transfer function for mode k (k in Mpc)
|
||||||
inline double TFfit_onek(double k, double *tf_baryon, double *tf_cdm)
|
inline double TFfit_onek(double k, double *tf_baryon, double *tf_cdm) const
|
||||||
/* Input: k -- Wavenumber at which to calculate transfer function, in Mpc^-1.
|
/* Input: k -- Wavenumber at which to calculate transfer function, in Mpc^-1.
|
||||||
*tf_baryon, *tf_cdm -- Input value not used; replaced on output if
|
*tf_baryon, *tf_cdm -- Input value not used; replaced on output if
|
||||||
the input was not NULL. */
|
the input was not NULL. */
|
||||||
|
@ -173,7 +173,7 @@ struct eisenstein_transfer
|
||||||
fc_ = (cp["Omega_m"] - cp["Omega_b"]) / cp["Omega_m"] ;
|
fc_ = (cp["Omega_m"] - cp["Omega_b"]) / cp["Omega_m"] ;
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double at_k(double k)
|
inline double at_k(double k) const
|
||||||
{
|
{
|
||||||
double tfb, tfcdm;
|
double tfb, tfcdm;
|
||||||
TFfit_onek(k * m_h0, &tfb, &tfcdm);
|
TFfit_onek(k * m_h0, &tfb, &tfcdm);
|
||||||
|
@ -211,17 +211,17 @@ 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)
|
inline double compute(double k, tf_type type) const
|
||||||
{
|
{
|
||||||
return etf_.at_k(k);
|
return etf_.at_k(k);
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double get_kmin(void)
|
inline double get_kmin(void) const
|
||||||
{
|
{
|
||||||
return 1e-5;
|
return 1e-5;
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double get_kmax(void)
|
inline double get_kmax(void) const
|
||||||
{
|
{
|
||||||
return 1.e5;
|
return 1.e5;
|
||||||
}
|
}
|
||||||
|
@ -303,17 +303,17 @@ public:
|
||||||
std::cerr << "WDM alpha = " << m_WDMalpha << std::endl;
|
std::cerr << "WDM alpha = " << m_WDMalpha << std::endl;
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double compute(double k, tf_type type)
|
inline double compute(double k, tf_type type) const
|
||||||
{
|
{
|
||||||
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_);
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double get_kmin(void)
|
inline double get_kmin(void) const
|
||||||
{
|
{
|
||||||
return 1e-4;
|
return 1e-4;
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double get_kmax(void)
|
inline double get_kmax(void) const
|
||||||
{
|
{
|
||||||
return 1.e4;
|
return 1.e4;
|
||||||
}
|
}
|
||||||
|
@ -348,7 +348,7 @@ public:
|
||||||
music::ilog.Print(" bino CDM: k_fs = %g, k_d = %g", kfs_, kd_);
|
music::ilog.Print(" bino CDM: k_fs = %g, k_d = %g", kfs_, kd_);
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double compute(double k, tf_type type)
|
inline double compute(double k, tf_type type) const
|
||||||
{
|
{
|
||||||
double kkfs = k / kfs_;
|
double kkfs = k / kfs_;
|
||||||
double kkfs2 = kkfs * kkfs;
|
double kkfs2 = kkfs * kkfs;
|
||||||
|
@ -362,12 +362,12 @@ public:
|
||||||
return 0.0;
|
return 0.0;
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double get_kmin(void)
|
inline double get_kmin(void) const
|
||||||
{
|
{
|
||||||
return 1e-4;
|
return 1e-4;
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double get_kmax(void)
|
inline double get_kmax(void) const
|
||||||
{
|
{
|
||||||
return 1.e8;
|
return 1.e8;
|
||||||
}
|
}
|
||||||
|
|
|
@ -197,7 +197,7 @@ public:
|
||||||
gsl_interp_accel_free(acc_dtot0);
|
gsl_interp_accel_free(acc_dtot0);
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double extrap_left(double k, const tf_type &type)
|
inline double extrap_left(double k, const tf_type &type) const
|
||||||
{
|
{
|
||||||
if (k < 1e-8)
|
if (k < 1e-8)
|
||||||
return 1.0;
|
return 1.0;
|
||||||
|
@ -246,7 +246,7 @@ public:
|
||||||
return pow(10.0, (v2 - v1) / dk * (delk) + v1);
|
return pow(10.0, (v2 - v1) / dk * (delk) + v1);
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double extrap_right(double k, const tf_type &type)
|
inline double extrap_right(double k, const tf_type &type) const
|
||||||
{
|
{
|
||||||
double v1(1.0), v2(1.0);
|
double v1(1.0), v2(1.0);
|
||||||
|
|
||||||
|
@ -294,7 +294,7 @@ public:
|
||||||
return pow(10.0, (v2 - v1) / dk * (delk) + v2);
|
return pow(10.0, (v2 - v1) / dk * (delk) + v2);
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double compute(double k, tf_type type)
|
inline double compute(double k, tf_type type) const
|
||||||
{
|
{
|
||||||
|
|
||||||
double lk = log10(k);
|
double lk = log10(k);
|
||||||
|
@ -335,12 +335,12 @@ public:
|
||||||
return 1.0;
|
return 1.0;
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double get_kmin(void)
|
inline double get_kmin(void) const
|
||||||
{
|
{
|
||||||
return pow(10.0, m_tab_k[0]);
|
return pow(10.0, m_tab_k[0]);
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double get_kmax(void)
|
inline double get_kmax(void) const
|
||||||
{
|
{
|
||||||
return pow(10.0, m_tab_k[m_tab_k.size() - 1]);
|
return pow(10.0, m_tab_k[m_tab_k.size() - 1]);
|
||||||
}
|
}
|
||||||
|
|
|
@ -161,7 +161,7 @@ public:
|
||||||
gsl_interp_accel_free (acc_theta_baryon);
|
gsl_interp_accel_free (acc_theta_baryon);
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double extrap_left( double k, const tf_type& type )
|
inline double extrap_left( double k, const tf_type& type ) const
|
||||||
{
|
{
|
||||||
if( k<1e-8 )
|
if( k<1e-8 )
|
||||||
return 1.0;
|
return 1.0;
|
||||||
|
@ -202,7 +202,7 @@ public:
|
||||||
return pow(10.0,(v2-v1)/dk*(delk)+v1);
|
return pow(10.0,(v2-v1)/dk*(delk)+v1);
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double extrap_right( double k, const tf_type& type )
|
inline double extrap_right( double k, const tf_type& type ) const
|
||||||
{
|
{
|
||||||
double v1(1.0), v2(1.0);
|
double v1(1.0), v2(1.0);
|
||||||
|
|
||||||
|
@ -242,7 +242,7 @@ public:
|
||||||
return pow(10.0,(v2-v1)/dk*(delk)+v2);
|
return pow(10.0,(v2-v1)/dk*(delk)+v2);
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double compute( double k, tf_type type ){
|
inline double compute( double k, tf_type type ) const{
|
||||||
|
|
||||||
double lk = log10(k);
|
double lk = log10(k);
|
||||||
|
|
||||||
|
@ -279,11 +279,11 @@ public:
|
||||||
return 1.0;
|
return 1.0;
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double get_kmin( void ){
|
inline double get_kmin( void ) const{
|
||||||
return pow(10.0,m_tab_k[0]);
|
return pow(10.0,m_tab_k[0]);
|
||||||
}
|
}
|
||||||
|
|
||||||
inline double get_kmax( void ){
|
inline double get_kmax( void ) const{
|
||||||
return pow(10.0,m_tab_k[m_tab_k.size()-1]);
|
return pow(10.0,m_tab_k[m_tab_k.size()-1]);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
|
@ -74,13 +74,13 @@ public:
|
||||||
virtual void intialise(void) {}
|
virtual void intialise(void) {}
|
||||||
|
|
||||||
//! compute value of transfer function at waven umber
|
//! compute value of transfer function at waven umber
|
||||||
virtual double compute(double k, tf_type type) = 0;
|
virtual double compute(double k, tf_type type) const = 0;
|
||||||
|
|
||||||
//! return maximum wave number allowed
|
//! return maximum wave number allowed
|
||||||
virtual double get_kmax(void) = 0;
|
virtual double get_kmax(void) const = 0;
|
||||||
|
|
||||||
//! return minimum wave number allowed
|
//! return minimum wave number allowed
|
||||||
virtual double get_kmin(void) = 0;
|
virtual double get_kmin(void) const = 0;
|
||||||
|
|
||||||
//! return if density transfer function is distinct for baryons and DM
|
//! return if density transfer function is distinct for baryons and DM
|
||||||
bool tf_is_distinct(void) const
|
bool tf_is_distinct(void) const
|
||||||
|
|
Loading…
Reference in a new issue