diff --git a/doc/0000-ccl_note/0000-ccl_note.pdf b/doc/0000-ccl_note/0000-ccl_note.pdf index b25cd133c..542328247 100644 Binary files a/doc/0000-ccl_note/0000-ccl_note.pdf and b/doc/0000-ccl_note/0000-ccl_note.pdf differ diff --git a/doc/0000-ccl_note/main.pdf b/doc/0000-ccl_note/main.pdf index b25cd133c..542328247 100644 Binary files a/doc/0000-ccl_note/main.pdf and b/doc/0000-ccl_note/main.pdf differ diff --git a/doc/0000-ccl_note/main.tex b/doc/0000-ccl_note/main.tex index f2d139c3f..29184603e 100644 --- a/doc/0000-ccl_note/main.tex +++ b/doc/0000-ccl_note/main.tex @@ -476,7 +476,7 @@ \subsubsection{Spline parameters \& the INI file} \item {\tt A$\_$SPLINE$\_$MIN$\_$PK}: maximum value for the logarithmic part and minimum value for the linear part of the scale factor spline. The default is {\tt A$\_$SPLINE$\_$MIN$\_$PK}$=0.1$. This only applies for the 2D power spectrum splines. \item {\tt K$\_$MAX$\_$SPLINE}: The maximum value of wavenumber considered in the spline. This is explained in more detail in the coming subsections. The default is {\tt K$\_$MAX$\_$SPLINE}$=50/$Mpc. In the case of the emulator's nonlinear power spectrum, the maximum wavenumber is {\tt K$\_$MAX$\_$SPLINE}$=5/$Mpc. \item {\tt K$\_$MAX}: The maximum value of wavenumber when integrating over $k$. The default is {\tt K$\_$MAX}$=1000$/Mpc. -\item {\tt K$\_$MIN$\_$DEFAULT}: The minimum value of wavenumber when integrating over $k$. The default is {\tt K$\_$MIN$\_$DEFAULT}$=5 \times 10^{-5}$/Mpc. +\item {\tt K$\_$MIN}: The minimum value of wavenumber when integrating over $k$. The default is {\tt K$\_$MIN}$=5 \times 10^{-5}$/Mpc. \item {\tt N$\_$K}: Number of bins per decade for the wavenumber. The default is {\tt N$\_$K}$=167$. \item {\tt N$\_$K$\_$3DCOR}: Number of bins per decade for the wavenumber in the 3d correlation function calculation. The default is {\tt N$\_$K$\_$3DCOR}~=~100,000. @@ -511,11 +511,11 @@ \subsubsection{Spline parameters \& the INI file} \subsubsection{Extrapolation for the nonlinear power spectrum} \label{sec:NLextrapol} -The computation of the power spectrum from {\tt CLASS} can be significantly sped up by extrapolating in the range $k>$~{\tt K$\_$MAX$\_$SPLINE} and $k<$~{\tt K$\_$MIN}. In this section, we describe the implementation of the extrapolation and the accuracy attained. These tests are performed in a flat $\Lambda$CDM cosmology with $\Omega_c=0.25$, $\Omega_b=0.05$, $A_s=2.1\times10^{-9}$, $h=0.7$ and $n_s=0.96$. +The computation of the power spectrum from {\tt CLASS} can be significantly sped up by extrapolating in the range $k>$~{\tt K$\_$MAX$\_$SPLINE} and $k<$~{\tt K$\_$MIN$\_$SPLINE}. In this section, we describe the implementation of the extrapolation and the accuracy attained. These tests are performed in a flat $\Lambda$CDM cosmology with $\Omega_c=0.25$, $\Omega_b=0.05$, $A_s=2.1\times10^{-9}$, $h=0.7$ and $n_s=0.96$. We first describe the extrapolation at high wavenumbers. The introduction of the parameter {\tt K$\_$MAX$\_$SPLINE} allows us to spline the matter power spectrum within the {\tt cosmo} structure up to that value of $k$ (in units of $1/$Mpc). A separate {\tt K$\_$MAX} parameter sets the limit for evaluation of the matter power spectrum. The range between {\tt K$\_$MAX$\_$SPLINE}$${\tt K$\_$MAX$\_$SPLINE}. The Taylor expansion gives +First, we compute the first and second derivative of the $\ln P(k,z)$ at $k_0={\rm \tt K\_MAX\_SPLINE}-2\Delta\ln k$ by computing the numerical derivatives by finite differences using GSL. The fiducial choice for $\Delta\ln k$ is $10^{-2}$. We then apply a second order Taylor expansion to extrapolate the matter power spectrum to $k>${\tt K$\_$MAX$\_$SPLINE}. The Taylor expansion gives % \begin{equation} \ln P(k,z) \simeq \ln P(k_0,z) + \frac{d\ln P}{d\ln k}(\ln k_0,z) (\ln k-\ln k_0) + \frac{1}{2} \frac{d^2\ln P}{d\ln k^2}(\ln k_0,z) (\ln k-\ln k_0)^2. @@ -524,11 +524,13 @@ \subsubsection{Extrapolation for the nonlinear power spectrum} The accuracy of this approximation is shown in Figure \ref{fig:NLextrapol} for redshifts $z=0$, $z=3$ and $z=20$. We compare the nonlinear matter power spectrum at these redshifts computed with the previously described approximation, to the matter power spectrum obtained by setting the power spectrum splines to high values. We find that for typical values of $\Delta \ln k=10^{-2}$ and {\tt K$\_$MAX$\_$SPLINE}$=50$/Mpc, $\ln P$ has converged to an accuracy that surpasses the expected impact of baryonic effects on the matter power spectrum at $k>10/$Mpc. (For an estimate of the impact of baryons on the total matter power spectrum, see \citealt{Schneider15}.) The lower {\tt K$\_$MAX$\_$SPLINE} is, the faster \ccl will run. The optimum choice of {\tt K$\_$MAX$\_$SPLINE} is left to the user for their particular application. -We also extrapolate the power spectrum at small wavenumbers within the static routine {\tt ccl\_power\_extrapol\_lowk}. In this case, the power spectrum below {\tt K$\_$MIN} is obtained by a power-law extrapolation with index $n_s$: +We also extrapolate the power spectrum at small wavenumbers within the static routine {\tt ccl\_power\_extrapol\_lowk}. In this case, the power spectrum below {\tt K$\_$MIN$\_$SPLINE} is obtained by a power-law extrapolation with index $n_s$: \begin{equation} - \log P(k<{\tt K\_MIN},z) = \log P({\tt K\_MIN},z) + n_s (\log k-\log{\tt K\_MIN}) + \log P(k<{\tt K\_MIN\_SPLINE},z) = \log P({\tt K\_MIN\_SPLINE},z) + n_s (\log k-\log{\tt K\_MIN\_SPLINE}) \end{equation} -The value adopted for {\tt K\_MIN} depends on the choice of power spectrum method. For CLASS and the nonlinear power spectrum, we adopt {\tt K\_MIN} that coincides with the smallest wavenumber output by CLASS, {\tt K\_MIN}$=7\times 10^{-6}$/Mpc.\footnote{For BBKS, the power spectrum is computed analytically at all $k$, there is no extrapolation. For the Eisenstein \& Hu implementation, the splines of the power spectrum span {\tt K\_MIN\_DEFAULT}$ 10^{-5}$, and to $A_{\mathrm s}$ normalization otherwise. -In the {\tt python} implementation, {\tt CCL} allows for either $\sigma 8$ or {\tt A\_s} to be passed as parameters. +In the {\tt python} implementation, {\tt CCL} allows for either $\sigma_8$ or {\tt A\_s} to be passed as parameters. \subsection{Angular power spectra} \label{sec:cl} diff --git a/examples/ccl_sample_baryons.c b/examples/ccl_sample_baryons.c index de9bc96f8..4330e7930 100644 --- a/examples/ccl_sample_baryons.c +++ b/examples/ccl_sample_baryons.c @@ -32,7 +32,7 @@ int main(int argc, char * argv[]) //The linear power spectrum is not changed when baryons are passed /*printf("Linear matter PS\n"); printf("# k [1/Mpc],P(k,z=0),P(k,z=1),P(k,z=2),P(k,z=3)\n"); - for (k = ccl_splines->K_MIN_DEFAULT; kK_MAX; k*=1.05) { + for (k = ccl_splines->K_MIN; kK_MAX; k*=1.05) { p = ccl_linear_matter_power(cosmo, k,1.0, &status); p1 = ccl_linear_matter_power(cosmo,k, a_at_z1,&status); p2 = ccl_linear_matter_power(cosmo,k, a_at_z2,&status); @@ -41,7 +41,7 @@ int main(int argc, char * argv[]) }*/ printf("# Total matter power spectrum\n"); printf("# k [1/Mpc],P(k,z=0),P(k,z=1),P(k,z=2),P(k,z=3)\n"); - for (k = ccl_splines->K_MIN_DEFAULT; kK_MAX; k*=1.05) { + for (k = ccl_splines->K_MIN; kK_MAX; k*=1.05) { p = ccl_nonlin_matter_power(cosmo, k,1.0,&status); p1 = ccl_nonlin_matter_power(cosmo,k, a_at_z1,&status); p2 = ccl_nonlin_matter_power(cosmo,k, a_at_z2,&status); diff --git a/examples/ccl_sample_power.c b/examples/ccl_sample_power.c index 05e90bfc8..612213cca 100644 --- a/examples/ccl_sample_power.c +++ b/examples/ccl_sample_power.c @@ -25,7 +25,7 @@ int main(int argc, char * argv[]) double a_at_z2=1./3.; double a_at_z3=0.25; if(cosmo->config.matter_power_spectrum_method==ccl_linear) { - for (k = ccl_splines->K_MIN_DEFAULT; kK_MAX; k*=1.05) { + for (k = ccl_splines->K_MIN; kK_MAX; k*=1.05) { p = ccl_linear_matter_power(cosmo, k,1.0, &status); p1 = ccl_linear_matter_power(cosmo,k, a_at_z1,&status); p2 = ccl_linear_matter_power(cosmo,k, a_at_z2,&status); @@ -35,7 +35,7 @@ int main(int argc, char * argv[]) } else { if(cosmo->config.matter_power_spectrum_method==ccl_halofit) { - for (k = ccl_splines->K_MIN_DEFAULT; kK_MAX; k*=1.05) { + for (k = ccl_splines->K_MIN; kK_MAX; k*=1.05) { p = ccl_nonlin_matter_power(cosmo, k,1.0,&status); p1 = ccl_nonlin_matter_power(cosmo,k, a_at_z1,&status); p2 = ccl_nonlin_matter_power(cosmo,k, a_at_z2,&status); diff --git a/include/ccl_params.h b/include/ccl_params.h index 40e97e9d8..08044c4dc 100644 --- a/include/ccl_params.h +++ b/include/ccl_params.h @@ -32,7 +32,7 @@ typedef struct ccl_spline_params { //k-splines and integrals double K_MAX_SPLINE; double K_MAX; - double K_MIN_DEFAULT; + double K_MIN; int N_K; int N_K_3DCOR; diff --git a/include/ccl_params.ini b/include/ccl_params.ini index 21551c177..a87a0acd8 100644 --- a/include/ccl_params.ini +++ b/include/ccl_params.ini @@ -35,7 +35,7 @@ A_SPLINE_NLOG_PK=11 ;These are in units of Mpc (no factor of h) K_MAX_SPLINE=50 K_MAX=1e3 -K_MIN_DEFAULT=5e-5 +K_MIN=5e-5 ; This is now the number of k per decade N_K=167 ; Number of k per decade for 3D-correlation diff --git a/src/ccl_cls.c b/src/ccl_cls.c index 896ec8822..2fb09fd1d 100644 --- a/src/ccl_cls.c +++ b/src/ccl_cls.c @@ -896,13 +896,13 @@ static double *get_lkarr(ccl_cosmology *cosmo,CCL_ClWorkspace *w, //First compute relevant k-range for this ell double kmin,kmax,lkmin,lkmax; if(l>w->l_limber) { - kmin=CCL_MAX(ccl_splines->K_MIN_DEFAULT,0.8*(l+0.5)/chimax); + kmin=CCL_MAX(ccl_splines->K_MIN,0.8*(l+0.5)/chimax); kmax=CCL_MIN(ccl_splines->K_MAX,1.2*(l+0.5)/chimin); } else { double xmin,xmax; limits_bessel(l,CCL_FRAC_RELEVANT,&xmin,&xmax); - kmin=CCL_MAX(ccl_splines->K_MIN_DEFAULT,xmin/chimax); + kmin=CCL_MAX(ccl_splines->K_MIN,xmin/chimax); kmax=CCL_MIN(ccl_splines->K_MAX,xmax/chimin); //Cap by maximum meaningful argument of the Bessel function kmax=CCL_MIN(kmax,2*(w->l_arr[w->n_ls-1]+0.5)/chimin); //Cap by 2 x inverse scale corresponding to l_max @@ -1079,7 +1079,7 @@ static void get_k_interval(ccl_cosmology *cosmo,CCL_ClWorkspace *w, if(ll_limber) { chimin=2*(l+0.5)/ccl_splines->K_MAX; - chimax=0.5*(l+0.5)/ccl_splines->K_MIN_DEFAULT; + chimax=0.5*(l+0.5)/ccl_splines->K_MIN; } else { if(cut_low_1) { @@ -1098,7 +1098,7 @@ static void get_k_interval(ccl_cosmology *cosmo,CCL_ClWorkspace *w, } else { chimin=0.5*(l+0.5)/ccl_splines->K_MAX; - chimax=2*(l+0.5)/ccl_splines->K_MIN_DEFAULT; + chimax=2*(l+0.5)/ccl_splines->K_MIN; } } diff --git a/src/ccl_core.c b/src/ccl_core.c index 01836ce23..5f997af33 100644 --- a/src/ccl_core.c +++ b/src/ccl_core.c @@ -64,7 +64,6 @@ void ccl_cosmology_read_config(void) // Use default ini file param_file = EXPAND_STR(__CCL_DATA_DIR__) "/ccl_params.ini"; } - if ((fconfig=fopen(param_file, "r")) == NULL) { char msg[256]; snprintf(msg, 256, "ccl_core.c: Failed to open config file: %s", param_file); @@ -103,7 +102,7 @@ void ccl_cosmology_read_config(void) if(strcmp(var_name,"A_SPLINE_NLOG_PK")==0) ccl_splines->A_SPLINE_NLOG_PK=(int) var_dbl; if(strcmp(var_name,"K_MAX_SPLINE")==0) ccl_splines->K_MAX_SPLINE=var_dbl; if(strcmp(var_name,"K_MAX")==0) ccl_splines->K_MAX=var_dbl; - if(strcmp(var_name,"K_MIN_DEFAULT")==0) ccl_splines->K_MIN_DEFAULT=var_dbl; + if(strcmp(var_name,"K_MIN")==0) ccl_splines->K_MIN=var_dbl; if(strcmp(var_name,"N_K")==0) ccl_splines->N_K=(int) var_dbl; // 3dcorr parameters if(strcmp(var_name,"N_K_3DCOR")==0) ccl_splines->N_K_3DCOR=(int) var_dbl; diff --git a/src/ccl_correlation.c b/src/ccl_correlation.c index 47cca7525..2533772a9 100644 --- a/src/ccl_correlation.c +++ b/src/ccl_correlation.c @@ -452,9 +452,9 @@ void ccl_correlation_3d(ccl_cosmology *cosmo, double a, double *k_arr,*pk_arr,*r_arr,*xi_arr; //number of data points for k and pk array - N_ARR=(int)(ccl_splines->N_K_3DCOR*log10(ccl_splines->K_MAX/ccl_splines->K_MIN_DEFAULT)); + N_ARR=(int)(ccl_splines->N_K_3DCOR*log10(ccl_splines->K_MAX/ccl_splines->K_MIN)); - k_arr=ccl_log_spacing(ccl_splines->K_MIN_DEFAULT,ccl_splines->K_MAX,N_ARR); + k_arr=ccl_log_spacing(ccl_splines->K_MIN,ccl_splines->K_MAX,N_ARR); if(k_arr==NULL) { *status=CCL_ERROR_MEMORY; strcpy(cosmo->status_message,"ccl_correlation.c: ccl_correlation_3d ran out of memory\n"); diff --git a/src/ccl_power.c b/src/ccl_power.c index d3629745d..156bfc11e 100644 --- a/src/ccl_power.c +++ b/src/ccl_power.c @@ -421,6 +421,7 @@ static void ccl_cosmology_compute_power_class(ccl_cosmology * cosmo, int * statu return; } + //These are the limits of the splining range cosmo->data.k_min_lin=2*exp(sp.ln_k[0]); cosmo->data.k_max_lin=ccl_splines->K_MAX_SPLINE; @@ -493,6 +494,7 @@ static void ccl_cosmology_compute_power_class(ccl_cosmology * cosmo, int * statu strcpy(cosmo->status_message ,"ccl_power.c: ccl_cosmology_compute_power_class(): K_MIN is less than CLASS's kmin. Not yet supported for nonlinear P(k).\n"); } + //These are the limits of the splining range cosmo->data.k_min_nl=2*exp(sp.ln_k[0]); cosmo->data.k_max_nl=ccl_splines->K_MAX_SPLINE; @@ -802,8 +804,9 @@ static double eh_power(ccl_parameters *params,eh_struct *eh,double k,int wiggled static void ccl_cosmology_compute_power_eh(ccl_cosmology * cosmo, int * status) { - cosmo->data.k_min_lin = ccl_splines->K_MIN_DEFAULT; - cosmo->data.k_min_nl = ccl_splines->K_MIN_DEFAULT; + //These are the limits of the splining range + cosmo->data.k_min_lin = ccl_splines->K_MIN; + cosmo->data.k_min_nl = ccl_splines->K_MIN; cosmo->data.k_max_lin = ccl_splines->K_MAX; cosmo->data.k_max_nl = ccl_splines->K_MAX; double kmin = cosmo->data.k_min_lin; @@ -979,8 +982,9 @@ TASK: provide spline for the BBKS power spectrum with baryonic correction static void ccl_cosmology_compute_power_bbks(ccl_cosmology * cosmo, int * status) { - cosmo->data.k_min_lin=ccl_splines->K_MIN_DEFAULT; - cosmo->data.k_min_nl=ccl_splines->K_MIN_DEFAULT; + //These are the limits of the splining range + cosmo->data.k_min_lin=ccl_splines->K_MIN; + cosmo->data.k_min_nl=ccl_splines->K_MIN; cosmo->data.k_max_lin=ccl_splines->K_MAX; cosmo->data.k_max_nl=ccl_splines->K_MAX; double kmin = cosmo->data.k_min_lin; @@ -1222,9 +1226,10 @@ static void ccl_cosmology_compute_power_emu(ccl_cosmology * cosmo, int * status) return; } + //These are the limits of the splining range cosmo->data.k_min_lin=2*exp(sp.ln_k[0]); cosmo->data.k_max_lin=ccl_splines->K_MAX_SPLINE; -//CLASS calculations done - now allocate CCL splines + //CLASS calculations done - now allocate CCL splines double kmin = cosmo->data.k_min_lin; double kmax = ccl_splines->K_MAX_SPLINE; //Compute nk from number of decades and N_K = # k per decade @@ -1286,6 +1291,7 @@ static void ccl_cosmology_compute_power_emu(ccl_cosmology * cosmo, int * status) } //Now start the NL computation with the emulator + //These are the limits of the splining range cosmo->data.k_min_nl=K_MIN_EMU; cosmo->data.k_max_nl=K_MAX_EMU; amin = A_MIN_EMU; //limit of the emulator @@ -1394,7 +1400,7 @@ INPUT: ccl_cosmology * cosmo, a, k [1/Mpc] TASK: extrapolate power spectrum at high k */ static double ccl_power_extrapol_highk(ccl_cosmology * cosmo, double k, double a, - gsl_spline2d * powerspl, double kmax, int * status) + gsl_spline2d * powerspl, double kmax_spline, int * status) { double log_p_1; double deltak=1e-2; //step for numerical derivative; @@ -1402,7 +1408,7 @@ static double ccl_power_extrapol_highk(ccl_cosmology * cosmo, double k, double a double lkmid; double lpk_kmid; - lkmid = log(kmax)-2*deltak; + lkmid = log(kmax_spline)-2*deltak; int gslstatus = gsl_spline2d_eval_e(powerspl, lkmid,a,NULL ,NULL ,&lpk_kmid); if(gslstatus != GSL_SUCCESS) { @@ -1432,16 +1438,16 @@ static double ccl_power_extrapol_highk(ccl_cosmology * cosmo, double k, double a } -/*------ ROUTINE: ccl_power_extrapol_hxighk ----- +/*------ ROUTINE: ccl_power_extrapol_lowk ----- INPUT: ccl_cosmology * cosmo, a, k [1/Mpc] TASK: extrapolate power spectrum at low k */ static double ccl_power_extrapol_lowk(ccl_cosmology * cosmo, double k, double a, - gsl_spline2d * powerspl, double kmin, int * status) + gsl_spline2d * powerspl, double kmin_spline, int * status) { double log_p_1; double deltak=1e-2; //safety step - double lkmin=log(kmin)+deltak; + double lkmin=log(kmin_spline)+deltak; double lpk_kmin; int gslstatus = gsl_spline2d_eval_e(powerspl,lkmin,a,NULL,NULL,&lpk_kmin); @@ -1669,7 +1675,7 @@ double ccl_sigmaR(ccl_cosmology *cosmo,double R, int *status) F.function=&sigmaR_integrand; F.params=∥ double sigma_R; - int gslstatus = gsl_integration_cquad(&F, log10(ccl_splines->K_MIN_DEFAULT), log10(ccl_splines->K_MAX), + int gslstatus = gsl_integration_cquad(&F, log10(ccl_splines->K_MIN), log10(ccl_splines->K_MAX), 0.0, ccl_gsl->INTEGRATION_SIGMAR_EPSREL, workspace,&sigma_R,NULL,NULL); if(gslstatus != GSL_SUCCESS) {