Problem with sigma_8

[miguelzuma]

I am experiencing a problem when the code computes sigma_8 in certain models:

Error in spectra_init =>spectra_init(L:1268) :error in spectra_pk(pba,ppt,ppm,pnl,psp); =>spectra_pk(L:2711) :error in spectra_sigma(pba,ppm,psp,8./pba->h,0.,&(psp->sigma8)); =>spectra_sigma(L:2801) :error in spectra_pk_at_k_and_z(pba,ppm,psp,k,z,&pk,pk_ic); =>spectra_pk_at_k_and_z(L:568) :condition ((k < 0.) || (k > exp(psp->ln_k[psp->ln_k_size-1]))) is true; k=5.505146e+00 out of bounds [0.000000e+00:5.505146e+00]

(the parameters.ini file is copied at the end) There is no problem for other models, e.g. when fixing omega_cdm = 0.110616 there is no problem.

I've tried to substitute (i=0;iln_k_size;i++) -> (i=0;iln_k_size -1 ;i++) in the loop for the integration of sigma (in spectra.c). That fixes the error, but gives a very bad value of sigma_8:

-> sigma8=1.48732 (computed till k = 7.86449 h/Mpc)

Here is the input file:

List of input/precision parameters actually read (all other parameters set to default values) Obtained with CLASS v2.1.2 (for developpers: svn version 6142M) This file can be used as the input file of another run h = 0.7 N_ur = 3.046 N_ncdm = 0 Omega_k = 0. Omega_fld = 0. YHe = BBN recombination = RECFAST reio_parametrization = reio_camb z_reio = 10. output = mPk,mTk non linear = halofit modes = s ic = ad gauge = synchronous P_k_ini type = analytic_Pk k_pivot = 0.002 A_s = 2.3e-9 n_s = 1. alpha_s = 0. P_k_max_h/Mpc = 1. zpk = 0., 1., 3.5 root = output/planck headers = yes format = class write background = yup write thermodynamics = no write parameters = yeap background_verbose = 1 thermodynamics_verbose = 1 perturbations_verbose = 1 transfer_verbose = 1 primordial_verbose = 1 spectra_verbose = 1 nonlinear_verbose = 1 lensing_verbose = 1 output_verbose = 1

[lesgourg]

Dear Miguel, as you could figure out this comes from a rounding error done by the computer (the operation k=exp(psp->ln_k[i]) is done in two different place and returns a different rounded number). It is computer-dependent. What you did is a possible way out, but after changing psp->ln_k_size into psp->ln_k_size-1 in the loop ( for (i=0;iln_k_size;i++) {...}, line 2795 of spectra.c), you must do the same for constancy in the list of arguments of the functions array_spline() and array_integrate_all_spline() which are just below. Finally, in line 2716, you should do it again:

        exp(psp->ln_k[psp->ln_k_size-2])/pba->h);

Then, the result will be correct, although sigma8 will be computed up to a slightly smaller k.

Another fix (that would preserve the calculation till the same k) would be to insert only one line in 2797, just after

for (i=0;iln_k_size;i++) { k=exp(psp->ln_k[i]);

The new line would be:

if (i == psp->ln_k_size-1) k*= 0.9999999;

It should work well but it is not terribly elegant. I think that I will use the first fix in the next release. I'll think about it a few more days.

[lesgourg]

Fixed in new release v2.3.0 using last method described above