From 9913c1db850304f544ae053a936f1886e431e7c1 Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Mon, 17 Jun 2019 15:07:04 +0100 Subject: [PATCH] working commit for baryon back-scaling --- example.conf | 9 +- external/class | 2 +- include/cosmology_calculator.hh | 1 + src/main.cc | 14 +- src/plugins/transfer_CLASS.cc | 248 ++++++++++++++++++++++++++++++-- 5 files changed, 251 insertions(+), 23 deletions(-) diff --git a/example.conf b/example.conf index c9c5de7..76e4df6 100644 --- a/example.conf +++ b/example.conf @@ -1,8 +1,9 @@ [setup] -GridRes = 64 +GridRes = 128 BoxLength = 300 -zstart = 49.0 -LPTorder = 3 +zstart = 1.0 +ztarget = 1.0 +LPTorder = 2 SymplecticPT = no DoFixing = no @@ -20,7 +21,7 @@ generator = NGENIC seed = 9001 [cosmology] -transfer = CLASS #eisenstein +transfer = eisenstein # CLASS #eisenstein Omega_m = 0.302 Omega_b = 0.045 Omega_L = 0.698 diff --git a/external/class b/external/class index b34d7f6..0b9f62d 160000 --- a/external/class +++ b/external/class @@ -1 +1 @@ -Subproject commit b34d7f6c2b72eab3a347c28e62298d62ca9dd69b +Subproject commit 0b9f62dabdbbc5815652bc172966510512cb082a diff --git a/include/cosmology_calculator.hh b/include/cosmology_calculator.hh index ad5a38b..4b5a7c6 100644 --- a/include/cosmology_calculator.hh +++ b/include/cosmology_calculator.hh @@ -106,6 +106,7 @@ public: csoca::ilog << "Wrote power spectrum at a=" << a << " to file \'" << fname << "\'" << std::endl; } + const CosmologyParameters &GetParams(void) const { return cosmo_param_; diff --git a/src/main.cc b/src/main.cc index 5f892c7..3041bbc 100644 --- a/src/main.cc +++ b/src/main.cc @@ -14,13 +14,13 @@ // initialise with "default" values namespace CONFIG{ -int MPI_thread_support = -1; -int MPI_task_rank = 0; -int MPI_task_size = 1; -bool MPI_ok = false; -bool MPI_threads_ok = false; -bool FFTW_threads_ok = false; -int num_threads = 1; + int MPI_thread_support = -1; + int MPI_task_rank = 0; + int MPI_task_size = 1; + bool MPI_ok = false; + bool MPI_threads_ok = false; + bool FFTW_threads_ok = false; + int num_threads = 1; } diff --git a/src/plugins/transfer_CLASS.cc b/src/plugins/transfer_CLASS.cc index 85b65b8..d23dbdb 100644 --- a/src/plugins/transfer_CLASS.cc +++ b/src/plugins/transfer_CLASS.cc @@ -9,7 +9,7 @@ #include #include #include - +#include #include #include @@ -19,13 +19,28 @@ #include #include +struct fluid_params{ + double H0, Omega_m, fb, fc; +}; + +// constexpr int nvar{4}; + +// define some operators on arrays +using vec4 = std::array; +vec4 operator*( double x, vec4 y ){ return {x*y[0],x*y[1],x*y[2],x*y[3]}; } +vec4 operator/( vec4 y, double x ){ return {y[0]/x,y[1]/x,y[2]/x,y[3]/x}; } +vec4 operator+( vec4 x, vec4 y ){ return {x[0]+y[0],x[1]+y[1],x[2]+y[2],x[3]+y[3]}; } +vec4 operator-( vec4 x, vec4 y ){ return {x[0]-y[0],x[1]-y[1],x[2]-y[2],x[3]-y[3]}; } + +void backscale_twofluid( double astart, double aend, vec4& instate, fluid_params& p, vec4& outstate ); + class transfer_CLASS_plugin : public TransferFunction_plugin { 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_; + double Omega_m_, Omega_b_, Omega_r_, N_ur_, zstart_, ztarget_, kmax_, kmin_, H0_, h_, fb_, fc_, gfac_; void ClassEngine_get_data( void ){ std::vector d_ncdm, t_ncdm, phi, psi; @@ -39,32 +54,159 @@ private: 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_Lambda",1.0-Omega_m_); + pars.add("Omega_k",0.0); + 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); // tnis 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"); +// lensing = no +// ic = ad +// gauge = synchronous +// P_k_ini type = analytic_Pk +// k_pivot = 0.05 +// A_s = 2.215e-9 +// n_s = 0.96 +// alpha_s = 0. pars.add("k_per_decade_for_pk",50); pars.add("k_per_decade_for_bao",50); + pars.add("tol_perturb_integration",1e-8); + pars.add("tol_background_integration",1e-9); pars.add("compute damping scale","yes"); - pars.add("z_reio",-1.0); // make sure reionisation is not included + // pars.add("evolver",1); + + pars.add("z_reio",10.0); // make sure reionisation is not included std::unique_ptr CE = std::make_unique(pars, false); CE->getTk(ztarget_, tab_lnk_, tab_dc_, tab_db_, d_ncdm, tab_dtot_, tab_tc_, tab_tb_, t_ncdm, tab_ttot_, phi, psi ); + double astart = 1.0/(1.0+ztarget_); + // assume dtot = fc * dc + fb * db -> dc = (dtot-fB * db)/fc + for( size_t i=0; i ztarget_ ){ + while( a > afinal ){ + state = RK4_step(state, a, -da*a); + a-=da*a; + da = std::min(da*a,a-afinal)/a; + ++istep; + } + }else{ + while( a < afinal ){ + state = RK4_step(state, a, da*a); + a+=da*a; + da = std::min(da*a,afinal-a)/a; + ++istep; + } + } + // std::cerr << "istep=" << istep << std::endl; + ofs << std::exp(tab_lnk_[i]) << " "; + for( int j=0; j<4; ++j ) + ofs << state0[j] << " " << state[j] << " "; + ofs << std::endl; + } + } + + ///////////////////////////////////////////////////////////////////////////////////////////// + ///////////////////////////////////////////////////////////////////////////////////////////// + + + // void BackScaleTF( void ){ + // vec4 state, endstate; + // fluid_params p; + // p.H0 = 100.0*h_; + // p.Omega_m = Omega_m_; + // p.fb = Omega_b_/Omega_m_; + // p.fc = 1.0-p.fb; + + // csoca::ilog << "Back-Scaling two-fluid TF from z=" << ztarget_ << " to z=" << zstart_ << "..." << std::endl; + // std::ofstream ofs("backsc.txt"); + // for( size_t i=0; iGetValue("cosmology","Omega_m"); Omega_b_ = pcf_->GetValue("cosmology","Omega_b"); N_ur_ = pcf_->GetValueSafe("cosmology","N_ur", 3.046); - ztarget_ = pcf_->GetValueSafe("cosmology","ztarget",0.0); + ztarget_ = pcf_->GetValueSafe("setup","ztarget",0.0); zstart_ = pcf_->GetValue("setup","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 + double Tcmb = 2.726; + double Omega_g = 4.48146636e-7 * std::pow(Tcmb,4) /h_/h_; + double Omega_nu = N_ur_ * (7./8) * std::pow(4./11,4./3) * Omega_g; + + Omega_r_ = 0.0; //Omega_g+Omega_nu; + std::cerr << "Omega_r_ = " << Omega_r_ << std::endl; + kmax_ = 2.0*M_PI/lbox * nres/2 * std::sqrt(3.0) * 2.0; // 200% of spatial diagonal + H0_ = 100.0*h_; + fb_ = Omega_b_/Omega_m_; + fc_ = 1.0-Omega_b_/Omega_m_; + gfac_ = 1.5*H0_*H0_*Omega_m_; + this->ClassEngine_get_data(); + + this->BackScaleTF(); gsl_ia_dtot_ = gsl_interp_accel_alloc(); gsl_ia_dc_ = gsl_interp_accel_alloc(); @@ -151,4 +306,75 @@ namespace { TransferFunction_plugin_creator_concrete creator("CLASS"); } + +// #include +// #include +// #include + + + + +// int fluid_eq( double a, const double y[], double f[], void* params ) +// { +// fluid_params *p = reinterpret_cast(params); +// double delta_c = y[0]; +// double theta_c = y[1]; +// double delta_b = y[2]; +// double theta_b = y[3]; +// // a=-a; + +// double adot = p->H0 * std::sqrt( p->Omega_m / a + (1.0-p->Omega_m) * a * a ); +// double Gfac = 1.5*p->H0*p->H0*p->Omega_m; + +// f[0] = (-theta_c / (a*adot)); +// f[1] = (-theta_c / a - Gfac * delta_c / (a*a*adot) ); +// // f[1] = (-theta_c / a - Gfac * delta_c / (a*a*adot) ); +// //f[1] = -(theta_c / a + Gfac * (p->fc * delta_c + p->fb * delta_b) / (a*a*adot)); +// //f[2] = -(theta_b / (a*adot)); +// //f[3] = -(theta_b / a + Gfac * (p->fc * delta_c + p->fb * delta_b) / (a*a*adot)); + +// return GSL_SUCCESS; +// } + +// void backscale_twofluid( double astart, double aend, vec4& instate, fluid_params& p, vec4& outstate ) +// { +// static const gsl_odeiv2_step_type * T = gsl_odeiv2_step_rk4; +// const size_t nvar = instate.size(); + +// gsl_odeiv2_step *s = gsl_odeiv2_step_alloc (T, 4); +// gsl_odeiv2_control *c = gsl_odeiv2_control_yp_new (1e-10,1e-10); +// gsl_odeiv2_evolve *e = gsl_odeiv2_evolve_alloc (4); + +// gsl_odeiv2_system sys = {fluid_eq, nullptr, nvar, reinterpret_cast (&p)}; + +// double yk[nvar]; +// for( size_t ivar=0;ivar