1
0
Fork 0
mirror of https://github.com/cosmo-sims/monofonIC.git synced 2024-09-19 17:03:45 +02:00

added CLASS transfer function plugin

This commit is contained in:
Oliver Hahn 2019-05-15 21:03:10 +02:00
parent 80867c6fd6
commit e989fbc1fb
7 changed files with 237 additions and 29 deletions

55
.gitignore vendored
View file

@ -1,2 +1,55 @@
build
.vscode
.vscode
src/CMakeFiles/3.12.2/CompilerIdC/CMakeCCompilerId.c
src/CMakeFiles/feature_tests.c
src/CMakeFiles/feature_tests.cxx
src/CMakeFiles/progress.marks
src/CMakeFiles/3.12.2/CMakeCCompiler.cmake
src/CMakeFiles/3.12.2/CMakeCXXCompiler.cmake
src/CMakeFiles/3.12.2/CMakeDetermineCompilerABI_C.bin
src/CMakeFiles/3.12.2/CMakeDetermineCompilerABI_CXX.bin
src/CMakeFiles/3.12.2/CMakeSystem.cmake
src/CMakeFiles/fastLPT.dir/build.make
src/CMakeFiles/FindMPI/test_mpi.cpp
src/CMakeFiles/FindMPI/test_mpi_C.bin
src/CMakeFiles/FindMPI/test_mpi_CXX.bin
src/CMakeFiles/FindOpenMP/OpenMPCheckVersion.c
src/CMakeFiles/FindOpenMP/OpenMPCheckVersion.cpp
src/CMakeFiles/FindOpenMP/OpenMPTryFlag.c
src/CMakeFiles/FindOpenMP/OpenMPTryFlag.cpp
src/CMakeFiles/FindOpenMP/ompver_C.bin
src/CMakeFiles/FindOpenMP/ompver_CXX.bin
src/CMakeFiles/fastLPT.dir/CXX.includecache
src/CMakeFiles/fastLPT.dir/DependInfo.cmake
src/CMakeFiles/fastLPT.dir/plugins/transfer_eisenstein.cc.o
src/CMakeFiles/3.12.2/CompilerIdCXX/a.out
src/CMakeFiles/fastLPT.dir/cmake_clean.cmake
src/CMakeFiles/fastLPT.dir/depend.internal
src/CMakeFiles/fastLPT.dir/depend.make
src/CMakeFiles/fastLPT.dir/flags.make
src/CMakeFiles/fastLPT.dir/grid_fft.cc.o
src/CMakeFiles/fastLPT.dir/link.txt
src/CMakeFiles/fastLPT.dir/logger.cc.o
src/CMakeFiles/fastLPT.dir/main.cc.o
src/CMakeFiles/fastLPT.dir/progress.make
src/CMakeFiles/fastLPT.dir/random_plugin.cc.o
src/CMakeFiles/fastLPT.dir/transfer_function_plugin.cc.o
src/CMakeFiles/fastLPT.dir/plugins/random_music.cc.o
src/CMakeFiles/fastLPT.dir/plugins/random_music_wnoise_generator.cc.o
src/CMakeFiles/feature_tests.bin
src/CMakeFiles/CMakeDirectoryInformation.cmake
src/CMakeFiles/CMakeOutput.log
src/CMakeFiles/Makefile.cmake
src/CMakeFiles/Makefile2
src/CMakeFiles/TargetDirectories.txt
src/CMakeFiles/cmake.check_cache
src/CMakeFiles/3.12.2/CompilerIdC/a.out
src/CMakeFiles/3.12.2/CompilerIdCXX/CMakeCXXCompilerId.cpp
src/CMakeFiles/hdf5/cmake_hdf5_test.c
src/fastLPT.dSYM/Contents/Info.plist
src/fastLPT.dSYM/Contents/Resources/DWARF/fastLPT
src/cmake_install.cmake
src/CMakeCache.txt
src/fastLPT
src/input_powerspec.txt
src/Makefile

View file

@ -1,6 +1,6 @@
[setup]
LPTorder = 2
GridRes = 128
LPTorder = 3
GridRes = 256
BoxLength = 300
zstart = 1.0
H0 = 70.0
@ -15,7 +15,7 @@ format = gadget2
filename = ics_gadget.dat
[cosmology]
transfer = eisenstein
transfer = CLASS #eisenstein
Omega_m = 0.302
Omega_b = 0.045
Omega_L = 0.698

View file

@ -8,6 +8,7 @@ else()
endif()
set(CLASS_INCLUDE_DIR ${CMAKE_CURRENT_LIST_DIR}/class/include)
set(CLASS_INCLUDE_CPP_DIR ${CMAKE_CURRENT_LIST_DIR}/class/cpp)
# list of object files generated by class
set(CLASS_OBJECT_FILES
@ -60,6 +61,8 @@ target_include_directories(class_cpp
macro(target_setup_class target_name)
target_include_directories(${target_name}
PRIVATE ${CLASS_INCLUDE_DIR})
target_include_directories(${target_name}
PRIVATE ${CLASS_INCLUDE_CPP_DIR})
target_link_libraries(${target_name} ${CLASS_OBJECT_FILES})
target_link_libraries(${target_name} class_cpp)
add_dependencies(${target_name} class_objects)

View file

@ -45,13 +45,13 @@ class TransferFunction_plugin
virtual ~TransferFunction_plugin(){};
//! 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
virtual double get_kmax(void) = 0;
virtual double get_kmax(void) const = 0;
//! 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
bool tf_is_distinct(void)

View file

@ -132,7 +132,7 @@ int main( int argc, char** argv )
{
// write power spectrum to a file
std::ofstream ofs("input_powerspec.txt");
for( double k=1e-4; k<1e4; k*=1.1 ){
for( double k=the_transfer_function->get_kmin(); k<the_transfer_function->get_kmax(); k*=1.1 ){
ofs << std::setw(16) << k
<< std::setw(16) << std::pow(the_cosmo_calc->GetAmplitude(k, total) * Dplus0, 2.0)
<< std::setw(16) << std::pow(the_cosmo_calc->GetAmplitude(k, total), 2.0)
@ -164,7 +164,8 @@ int main( int argc, char** argv )
//======================================================================
//... compute 1LPT displacement potential ....
// phi = - delta / k^2
csoca::ilog << "Computing phi(1) term..." << std::endl;
double wtime = get_wtime();
csoca::ilog << "Computing phi(1) term..." << std::flush;
phi.FourierTransformForward();
phi.apply_function_k_dep([&](auto x, auto k) -> ccomplex_t {
@ -174,23 +175,25 @@ int main( int argc, char** argv )
ccomplex_t delta = x * the_cosmo_calc->GetAmplitude(kmod, total);
return -delta / (kmod * kmod) * phifac / volfac;
});
phi.zero_DC_mode();
csoca::ilog << " took " << get_wtime()-wtime << "s" << std::endl;
//////////////////////////////////////////////////////////////////////////////////////////////////
auto assign_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return res; };
auto add_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return val+res; };
auto add2_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return val+2.0*res; };
auto sub_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return val-res; };
auto sub2_op = []( ccomplex_t res, ccomplex_t val ) -> ccomplex_t{ return val-2.0*res; };
csoca::ilog << "Computing phi(2) term..." << std::endl;
wtime = get_wtime();
csoca::ilog << "Computing phi(2) term..." << std::flush;
// Compute the source term for phi(2)
Conv.convolve_SumHessians( phi, {0,0}, phi, {1,1}, {2,2}, phi2, assign_op );
Conv.convolve_Hessians( phi, {1,1}, phi, {2,2}, phi2, add_op );
Conv.convolve_Hessians( phi, {0,1}, phi, {0,1}, phi2, sub_op );
Conv.convolve_Hessians( phi, {0,2}, phi, {0,2}, phi2, sub_op );
Conv.convolve_Hessians( phi, {1,2}, phi, {1,2}, phi2, sub_op );
csoca::ilog << " took " << get_wtime()-wtime << "s" << std::endl;
phi2.FourierTransformForward();
phi2.apply_function_k_dep([&](auto x, auto k) {
@ -202,7 +205,8 @@ int main( int argc, char** argv )
//======================================================================
//... compute 3LPT displacement potential
csoca::ilog << "Computing phi(3a) term..." << std::endl;
wtime = get_wtime();
csoca::ilog << "Computing phi(3a) term..." << std::flush;
Conv.convolve_SumHessians( phi, {0,0}, phi2, {1,1}, {2,2}, phi3a, assign_op );
Conv.convolve_SumHessians( phi, {1,1}, phi2, {2,2}, {0,0}, phi3a, add_op );
Conv.convolve_SumHessians( phi, {2,2}, phi2, {0,0}, {1,1}, phi3a, add_op );
@ -213,6 +217,8 @@ int main( int argc, char** argv )
phi3a.apply_function_k_dep([&](auto x, auto k) {
return 0.5 * x;
});
csoca::ilog << " took " << get_wtime()-wtime << "s" << std::endl;
phi3a.FourierTransformForward();
phi3a.apply_function_k_dep([&](auto x, auto k) {
@ -221,7 +227,8 @@ int main( int argc, char** argv )
});
phi3a.zero_DC_mode();
csoca::ilog << "Computing phi(3b) term..." << std::endl;
wtime = get_wtime();
csoca::ilog << "Computing phi(3b) term..." << std::flush;
Conv.convolve_Hessians( phi, {0,0}, phi, {1,1}, phi, {2,2}, phi3b, assign_op );
Conv.convolve_Hessians( phi, {0,1}, phi, {0,2}, phi, {1,2}, phi3b, add2_op );
Conv.convolve_Hessians( phi, {1,2}, phi, {1,2}, phi, {0,0}, phi3b, sub_op );
@ -234,6 +241,7 @@ int main( int argc, char** argv )
return x * (-1.0 / kmod2) * phifac / phifac / phifac/phifac;
});
phi3b.zero_DC_mode();
csoca::ilog << " took " << get_wtime()-wtime << "s" << std::endl;
///////////////////////////////////////////////////////////////////////
// we store the densities here if we compute them

View file

@ -0,0 +1,144 @@
// transfer_CLASS.cc - This file is part of MUSIC -
// a code to generate multi-scale initial conditions for cosmological simulations
// Copyright (C) 2019 Oliver Hahn
#include <cmath>
#include <string>
#include <vector>
#include <memory>
#include <ClassEngine.hh>
#include <general.hh>
#include <config_file.hh>
#include <transfer_function_plugin.hh>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
const double tiny = 1e-30;
class transfer_CLASS_plugin : public TransferFunction_plugin {
private:
std::vector<double> tab_lnk_, tab_dtot_, tab_dc_, tab_db_, tab_ttot_, tab_tc_, tab_tb_;
gsl_interp_accel *gsl_ia_dtot_, *gsl_ia_dc_, *gsl_ia_db_, *gsl_ia_ttot_, *gsl_ia_tc_, *gsl_ia_tb_;
gsl_spline *gsl_sp_dtot_, *gsl_sp_dc_, *gsl_sp_db_, *gsl_sp_ttot_, *gsl_sp_tc_, *gsl_sp_tb_;
double Omega_m_, Omega_b_, zstart_, ztarget_, kmax_, kmin_, h_;
void ClassEngine_get_data( void ){
std::vector<double> d_g, d_ur, t_g, t_ur, phi, psi;
csoca::ilog << "Computing transfer function via ClassEngine..." << std::flush;
double wtime = get_wtime();
ClassParams pars;
double target_redshift = 0.0;
pars.add("z_pk",target_redshift);
pars.add("P_k_max_h/Mpc", kmax_);
pars.add("h",h_);
pars.add("Omega_b",Omega_b_);
pars.add("Omega_ur",0.0);
pars.add("Omega_cdm",Omega_m_-Omega_b_);
pars.add("Omega_Lambda",1.0-Omega_m_);
pars.add("A_s",2.42e-9);
pars.add("n_s",.96);
pars.add("output","dTk,vTk");
pars.add("YHe",0.25);
std::unique_ptr<ClassEngine> CE = std::make_unique<ClassEngine>(pars, false);
CE->getTk(target_redshift, tab_lnk_, d_g, tab_db_, tab_dc_, d_ur, tab_dtot_,
phi, psi, t_g, tab_tb_, t_ur, tab_ttot_);
tab_tc_ = tab_ttot_;
#warning need to fix CDM velocities
wtime = get_wtime() - wtime;
csoca::ilog << " took " << wtime << " s / " << tab_lnk_.size() << " modes." << std::endl;
}
public:
explicit transfer_CLASS_plugin( ConfigFile &cf)
: TransferFunction_plugin(cf)
{
h_ = cf.GetValue<double>("cosmology","H0") / 100.0;
Omega_m_ = cf.GetValue<double>("cosmology","Omega_m");
Omega_b_ = cf.GetValue<double>("cosmology","Omega_b");
zstart_ = cf.GetValue<double>("setup","zstart");
ztarget_ = cf.GetValueSafe<double>("cosmology","ztarget",0.0);
kmax_ = 1000.0;
this->ClassEngine_get_data();
gsl_ia_dtot_ = gsl_interp_accel_alloc();
gsl_ia_dc_ = gsl_interp_accel_alloc();
gsl_ia_db_ = gsl_interp_accel_alloc();
gsl_ia_ttot_ = gsl_interp_accel_alloc();
gsl_ia_tc_ = gsl_interp_accel_alloc();
gsl_ia_tb_ = gsl_interp_accel_alloc();
gsl_sp_dtot_ = gsl_spline_alloc(gsl_interp_cspline, tab_lnk_.size());
gsl_sp_dc_ = gsl_spline_alloc(gsl_interp_cspline, tab_lnk_.size());
gsl_sp_db_ = gsl_spline_alloc(gsl_interp_cspline, tab_lnk_.size());
gsl_sp_ttot_ = gsl_spline_alloc(gsl_interp_cspline, tab_lnk_.size());
gsl_sp_tc_ = gsl_spline_alloc(gsl_interp_cspline, tab_lnk_.size());
gsl_sp_tb_ = gsl_spline_alloc(gsl_interp_cspline, tab_lnk_.size());
gsl_spline_init(gsl_sp_dtot_, &tab_lnk_[0], &tab_dtot_[0], tab_lnk_.size());
gsl_spline_init(gsl_sp_dc_, &tab_lnk_[0], &tab_dc_[0], tab_lnk_.size());
gsl_spline_init(gsl_sp_db_, &tab_lnk_[0], &tab_db_[0], tab_lnk_.size());
gsl_spline_init(gsl_sp_ttot_, &tab_lnk_[0], &tab_ttot_[0], tab_lnk_.size());
gsl_spline_init(gsl_sp_tc_, &tab_lnk_[0], &tab_tc_[0], tab_lnk_.size());
gsl_spline_init(gsl_sp_tb_, &tab_lnk_[0], &tab_tb_[0], tab_lnk_.size());
kmin_ = std::exp(tab_lnk_[0]);
tf_distinct_ = true; // [150612SH: different density between CDM v.s. Baryon]
tf_withvel_ = true; // [150612SH: using velocity transfer function]
}
~transfer_CLASS_plugin(){
gsl_spline_free(gsl_sp_dtot_);
gsl_spline_free(gsl_sp_dc_);
gsl_spline_free(gsl_sp_db_);
gsl_spline_free(gsl_sp_ttot_);
gsl_spline_free(gsl_sp_tc_);
gsl_spline_free(gsl_sp_tb_);
gsl_interp_accel_free(gsl_ia_dtot_);
gsl_interp_accel_free(gsl_ia_dc_);
gsl_interp_accel_free(gsl_ia_db_);
gsl_interp_accel_free(gsl_ia_ttot_);
gsl_interp_accel_free(gsl_ia_tc_);
gsl_interp_accel_free(gsl_ia_tb_);
}
inline double compute(double k, tf_type type) const {
gsl_spline *splineT = nullptr;
gsl_interp_accel *accT = nullptr;
switch(type){
case total: splineT = gsl_sp_dtot_; accT = gsl_ia_dtot_; break;
case cdm: splineT = gsl_sp_dc_; accT = gsl_ia_dc_; break;
case baryon: splineT = gsl_sp_db_; accT = gsl_ia_db_; break;
case vtotal: splineT = gsl_sp_ttot_; accT = gsl_ia_ttot_; break;
case vcdm: splineT = gsl_sp_tc_; accT = gsl_ia_tc_; break;
case vbaryon: splineT = gsl_sp_tb_; accT = gsl_ia_tb_; break;
default:
throw std::runtime_error("Invalid type requested in transfer function evaluation");
}
double d = (k<=kmin_)? gsl_spline_eval(splineT, std::log(kmin_), accT)
: gsl_spline_eval(splineT, std::log(k), accT);
return -d/(k*k);
}
inline double get_kmin(void) const { return std::exp(tab_lnk_[0]); }
inline double get_kmax(void) const { return std::exp(tab_lnk_[tab_lnk_.size()-1]); }
};
namespace {
TransferFunction_plugin_creator_concrete<transfer_CLASS_plugin> creator("CLASS");
}

View file

@ -36,13 +36,13 @@ struct eisenstein_transfer
alpha_gamma; /* Gamma suppression in approximate TF */
template <typename T>
inline T SQR(T x) { return x * x; }
inline T SQR(T x) const { return x * x; }
template <typename T>
inline T CUBE(T x) { return x * SQR(x); }
inline T CUBE(T x) const { return x * SQR(x); }
template<typename T>
inline T POW4( T x ){ return SQR(SQR(x)); }
inline T POW4( T x ) const { return SQR(SQR(x)); }
//! private member function: sets internal quantities for Eisenstein & Hu fitting
void TFset_parameters(double omega0hh, double f_baryon, double Tcmb)
@ -111,7 +111,7 @@ struct eisenstein_transfer
}
//! 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.
*tf_baryon, *tf_cdm -- Input value not used; replaced on output if
the input was not NULL. */
@ -177,7 +177,7 @@ struct eisenstein_transfer
fc_ = (Omega_m-Omega_b)/(Omega_m) ;
}
inline double at_k( double k )
inline double at_k( double k ) const
{
double tfb, tfcdm;
TFfit_onek( k*m_h0, &tfb, &tfcdm );
@ -218,15 +218,15 @@ public:
}
//! 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 );
}
inline double get_kmin( void ){
inline double get_kmin( void ) const{
return 1e-4;
}
inline double get_kmax( void ){
inline double get_kmax( void ) const{
return 1.e4;
}
@ -313,16 +313,16 @@ public:
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_);
}
inline double get_kmin( void ){
inline double get_kmin( void ) const{
return 1e-4;
}
inline double get_kmax( void ){
inline double get_kmax( void ) const{
return 1.e4;
}
@ -358,7 +358,7 @@ public:
//LOGINFO(" 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 kkfs2 = kkfs*kkfs;
@ -372,11 +372,11 @@ public:
return 0.0;
}
inline double get_kmin( void ){
inline double get_kmin( void ) const{
return 1e-4;
}
inline double get_kmax( void ){
inline double get_kmax( void ) const{
return 1.e8;
}