diff --git a/example.conf b/example.conf index 58fc969..2ebaeac 100644 --- a/example.conf +++ b/example.conf @@ -4,7 +4,7 @@ GridRes = 128 # length of the box in Mpc/h BoxLength = 250 # starting redshift -zstart = 49.0 +zstart = 100.0 # order of the LPT to be used (1,2 or 3) LPTorder = 3 # also do baryon ICs? @@ -15,8 +15,9 @@ DoFixing = no ParticleLoad = sc [cosmology] -#transfer = CLASS -transfer = eisenstein +transfer = CLASS +ztarget = 100.0 +#transfer = eisenstein Omega_m = 0.302 Omega_b = 0.045 Omega_L = 0.698 diff --git a/src/plugins/transfer_CLASS.cc b/src/plugins/transfer_CLASS.cc index 85b65b8..3469b7e 100644 --- a/src/plugins/transfer_CLASS.cc +++ b/src/plugins/transfer_CLASS.cc @@ -25,7 +25,13 @@ private: std::vector 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_, N_ur_, zstart_, ztarget_, kmax_, kmin_, h_; + + // single fluid growing/decaying mode decomposition + gsl_interp_accel *gsl_ia_Cplus_, *gsl_ia_Cminus_; + gsl_spline *gsl_sp_Cplus_, *gsl_sp_Cminus_; + std::vector tab_Cplus_, tab_Cminus_; + + double Omega_m_, Omega_b_, N_ur_, zstart_, ztarget_, kmax_, kmin_, h_, astart_, atarget_; void ClassEngine_get_data( void ){ std::vector d_ncdm, t_ncdm, phi, psi; @@ -37,25 +43,46 @@ private: pars.add("extra metric transfer functions", "yes"); pars.add("z_pk",ztarget_); pars.add("P_k_max_h/Mpc", kmax_); + pars.add("h",h_); pars.add("Omega_b",Omega_b_); // pars.add("Omega_k",0.0); // pars.add("Omega_ur",0.0); pars.add("N_ur",N_ur_); pars.add("Omega_cdm",Omega_m_-Omega_b_); - pars.add("Omega_Lambda",1.0-Omega_m_); - // pars.add("Omega_fld",0.0); - // pars.add("Omega_scf",0.0); + pars.add("Omega_k",0.0); + // pars.add("Omega_Lambda",1.0-Omega_m_); + pars.add("Omega_fld",0.0); + pars.add("Omega_scf",0.0); + pars.add("A_s",2.42e-9); - pars.add("n_s",.96); // tnis doesn't matter for TF + pars.add("n_s",.961); // this doesn't matter for TF pars.add("output","dTk,vTk"); pars.add("YHe",0.248); + pars.add("lensing","no"); + pars.add("alpha_s",0.0); + pars.add("P_k_ini type","analytic_Pk"); + pars.add("gauge","synchronous"); + + pars.add("k_per_decade_for_pk",100); + pars.add("k_per_decade_for_bao",100); - pars.add("k_per_decade_for_pk",50); - pars.add("k_per_decade_for_bao",50); pars.add("compute damping scale","yes"); pars.add("z_reio",-1.0); // make sure reionisation is not included + pars.add("tol_perturb_integration",1.e-8); + pars.add("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 + pars.add("hyper_flat_approximation_nu", 7000. ); + pars.add("transfer_neglect_delta_k_S_t0", 0.17 ); + pars.add("transfer_neglect_delta_k_S_t1", 0.05 ); + pars.add("transfer_neglect_delta_k_S_t2", 0.17 ); + pars.add("transfer_neglect_delta_k_S_e", 0.13 ); + pars.add("delta_l_max", 1000 ); + + std::unique_ptr CE = std::make_unique(pars, false); CE->getTk(ztarget_, tab_lnk_, tab_dc_, tab_db_, d_ncdm, tab_dtot_, @@ -74,7 +101,9 @@ public: Omega_b_ = pcf_->GetValue("cosmology","Omega_b"); N_ur_ = pcf_->GetValueSafe("cosmology","N_ur", 3.046); ztarget_ = pcf_->GetValueSafe("cosmology","ztarget",0.0); + atarget_ = 1.0/(1.0+ztarget_); zstart_ = pcf_->GetValue("setup","zstart"); + astart_ = 1.0/(1.0+zstart_); double lbox = pcf_->GetValue("setup","BoxLength"); int nres = pcf_->GetValue("setup","GridRes"); kmax_ = 2.0*M_PI/lbox * nres/2 * sqrt(3) * 2.0; // 120% of spatial diagonal @@ -102,6 +131,33 @@ public: 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()); + //-------------------------------------------------------------------------- + // single fluid growing/decaying mode decomposition + //-------------------------------------------------------------------------- + gsl_ia_Cplus_ = gsl_interp_accel_alloc(); + gsl_ia_Cminus_ = gsl_interp_accel_alloc(); + + gsl_sp_Cplus_ = gsl_spline_alloc(gsl_interp_cspline, tab_lnk_.size()); + gsl_sp_Cminus_ = gsl_spline_alloc(gsl_interp_cspline, tab_lnk_.size()); + + tab_Cplus_.assign(tab_lnk_.size(),0); + tab_Cminus_.assign(tab_lnk_.size(),0); + + std::ofstream ofs("grow_decay.txt"); + + for( size_t i=0; i