Matter power spectrum for small scales

[aslanyan]

I’m having a problem running class to get the matter power spectrum for very small scales. I’m running the latest version 2.4.2. I just changed the explanatory.ini file to include the matter power spectrum in the output (line 292) and the k_max value for matter power spectrum (line 598). Currently using P_k_max_h/Mpc = 1e3 but I’ll need to go to even smaller scales. Running class I get the following error:

Running CLASS version v2.4.2 Computing background -> age = 13.795359 Gyr -> conformal age = 14165.045412 Mpc Computing thermodynamics with Y_He=0.2477 -> recombination at z = 1089.267451 corresponding to conformal time = 280.576042 Mpc with comoving sound horizon = 144.695940 Mpc angular diameter distance = 12.734921 Mpc and sound horizon angle 100*theta_s = 1.042142 -> baryon drag stops at z = 1059.171266 corresponding to conformal time = 286.488461 Mpc with comoving sound horizon rs = 147.376600 Mpc -> reionization with optical depth = 0.092473 corresponding to conformal time = 4255.316282 Mpc Computing sources

Error in perturb_init =>perturb_init(L:363) :error in perturb_solve(ppr, pba, pth, ppt, index_md, index_ic, index_k, pppw[thread]); =>perturb_solve(L:2372) :error in generic_evolver(perturb_derivs, interval_limit[index_interval], interval_limit[index_interval+1], ppw->pv->y, ppw->pv->used_in_sources, ppw->pv->pt_size, &ppaw, ppr->tol_perturb_integration, ppr->smallest_allowed_variation, perturb_timescale, ppr->perturb_integration_stepsize, ppt->tau_sampling, tau_actual_size, perturb_sources, perhaps_print_variables, ppt->error_message); =>evolver_ndf15(L:461) :condition (absh <= hmin) is true; Step size too small: step:2.22045e-16, minimum:2.22045e-16, in interval: [1.14856:364.763]

[mcataneo]

Hi,

I ran precisely in the same issue lately and if speed is not one of your concerns you may set the

'evolver' = 0

in either the ini file or the pre file. I just started using CLASS so I'm not an expert, maybe there's a more efficient solution I'm not aware of.

Cheers, Matteo

[lesgourg]

Hi guys, thanks for pointing at the problem, and at this cheap way to avoid it. However the fact that there can be an error when using the default 'evolver=1' clearly suggests that there is some kind of bug, which shows up only when computing the P(k) at very large k (>1000h/Mpc) like you are doing.

I just found one such bug, and updated the branch 2.4 on github accordingly. This fix will go to a tagged release of the master branch very soon, but not immediately. However you could immediately try the new version of the branch 2.4 - I guess that you know how to do that, but let me recall it just in case:

If you are working with the master branch (or your own customised branch), you need to commit all your changes in order to be allowed to leave that branch. Then you can type:

git checkout 2.4

In principle git should tell you that it will now track the existing branch 'origin/2.4' Then you can type

git pull origin 2.4

You now have the latest 2.4 branch (last change this afternoon). After a 'make clean', 'make -j', you have the new code ready. You can check whether you still have the problem when you run with your own input file computing the P(k) up to 1000, with the default 'evolver=1'. I hope that now, things will work, but I am not completely sure, because we have just found a second bug (causing an error on some machines but not all) that we will try to solve it in the next hours/days. I don't know if your previous error came from the bug that has already been solved, or from the other one. In any case we will tell you when the second bug is fixed.

Fortunately in the meantime you can work with 'evolver=0' as suggested by Matteo. Sorry for these issues (they were difficult to find because we rarely use such high k values when we do debugging tests - but we should from now on).

[aslanyan]

Thanks for the comments! I tried branch 2.4 but it didn't seem to help - I still get the same error. Setting "evolver=0" seems to work up to certain values of k but at some point (around k=10^5) the resulting power spectrum looks wrong. After the peak the power spectrum is decreasing as a function of k but for large enough values it starts sharply increasing again. This must be a numerical issue. Also, "evolver=0" makes it extremely slow, especially with very large values of k.

I basically just need the matter power spectrum up to k~10^7 or possibly beyond, and so far I haven't been able to do so with CLASS. I understand that this value is very large and CLASS isn't necessarily designed to handle this case so not sure if I'd call this a bug, but it would definitely be nice to have that feature if possible.

[ThomasTram]

When I use branch 2.4, I do not get any errors using the standard evolver, and I can compute P(k) up to k=10^4 h/Mpc in a few minutes. I tried running for 10 hours on 16 cores using P_k_max_h/Mpc = 1e7, which did not finish but did not throw any small step size error. Such a large value of k poses several fundamental problems: 1) The frequency of oscillations in the perturbations scale with k, so increasing k by a factor 10 also increases the number of time steps by 10. 2) Initial conditions for the perturbations has to be imposed when the baryons and photons are tightly coupled. For large k, this is extremely early. 3) The memory requirements of the code also scales with k, since it needs to store the source functions as a function of tau which sufficient resolution. This is in principle not necessary for the P(k), so that could be fixed.

In the end, I think you are better off extending the matter power spectrum to small values by analytic means, for instance by using Eqs. (2.30, B.3, C.2-C.5) in http://arxiv.org/pdf/1312.5301.pdf.

[JinsuKim]

Hi,

I ran the code with the same setup. Have you tried "primordial=yes" option (L708)? It seems that this option allows me to avoid the bug.

Cheers, Jinsu