CCL Halofit power spectrum significantly different from CLASS result

[philbull]

In certain circumstances, the CCL Halofit power spectrum is significantly different from the CLASS Halofit power spectrum, even though CCL is calling CLASS to calculate it in the first place. This seems to be most significant at higher k and higher z when Omega_cdm is relatively low, although this is not necessarily the only time that the issue occurs.

The plot below demonstrates the problem. The solid red/blue lines are the CLASS Halofit power spectrum, obtained using the classy Python wrapper (I also confirmed that the CLASS executable produces the same results). The dashed yellow/cyan lines are the CCL Halofit power spectrum, obtained using ccl_nonlin_matter_power(). The agreement seems to be pretty good at all k for z <= 1, but degrades at higher redshifts at higher k.

This particular example was produced with the following CLASS parameter values:

    'h':                5.537500e-01,
    'n_s':              9.514000e-01,
    'A_s':              1.955000e-09,
    'Omega_cdm':        1.770000e-01,
    'Omega_b':          3.347000e-02,
    'z_pk':             '0.000,0.500,1.000,1.500,2.000,2.500',
    'non linear':       'halofit',
    'Omega_k':          0.0,
    'N_ur':             3.046000,
    'N_ncdm':           0.000000,
    'T_cmb':            2.725,

[philbull]

Here's the code I used to produce the plot (you need classy installed; it's distributed with CLASS):

import numpy as np
import pylab as P
from classy import Class
import pyccl as ccl

# Define CLASS cosmology
pclass = {
    'output':           'mPk',
    'lensing':          'no',

    'h':                5.537500e-01,
    'n_s':              9.514000e-01,
    'A_s':              1.955000e-09,
    'Omega_cdm':        1.770000e-01,
    'Omega_b':          3.347000e-02,

    'z_pk':             '0.000,0.500,1.000,1.500,2.000,2.500',
    'non linear':       'halofit',
    'Omega_k':          0.0,
    'N_ur':             3.046000,
    'N_ncdm':           0.000000,
    'T_cmb':            2.725,
}

# Define CCL cosmology
pccl = {key: pclass[key] for key in ['A_s', 'n_s', 'h', 'Omega_b']}
pccl['Omega_c'] = pclass['Omega_cdm']
pccl['transfer_function'] = 'boltzmann'
pccl['matter_power_spectrum'] = 'halofit'

# Define k array and redshift
k = np.logspace(-4., np.log10(2.), 300)

# Calculate CLASS power spectrum
cclass = Class()
cclass.set(pclass)
cclass.compute()

pk_class_nl0 = np.array([cclass.pk(_k, 0.) for _k in k])
pk_class_nl1 = np.array([cclass.pk(_k, 0.5) for _k in k])
pk_class_nl2 = np.array([cclass.pk(_k, 1.) for _k in k])
#pk_class_nl2a = np.array([cclass.pk(_k, 1.2) for _k in k])
pk_class_nl3 = np.array([cclass.pk(_k, 1.5) for _k in k])
pk_class_nl4 = np.array([cclass.pk(_k, 2.) for _k in k])

cclass.struct_cleanup()

# Calculate CCL power spectrum
cosmo = ccl.Cosmology(**pccl)
pk_ccl_nl0 = ccl.nonlin_matter_power(cosmo, k, 1./(1.+0.))
pk_ccl_nl1 = ccl.nonlin_matter_power(cosmo, k, 1./(1.+0.5))
pk_ccl_nl2 = ccl.nonlin_matter_power(cosmo, k, 1./(1.+1.))
#pk_ccl_nl2a = ccl.nonlin_matter_power(cosmo, k, 1./(1.+1.2))
pk_ccl_nl3 = ccl.nonlin_matter_power(cosmo, k, 1./(1.+1.5))
pk_ccl_nl4 = ccl.nonlin_matter_power(cosmo, k, 1./(1.+2.))

# Plot power spectrum
P.subplot(111)
P.plot(k, pk_class_nl0, 'r-', lw=1.8, label='CLASS Halofit z=0.0')
P.plot(k, pk_class_nl1, 'b-', lw=1.8, label='CLASS Halofit z=0.5')
P.plot(k, pk_class_nl2, 'r-', lw=1.8, label='CLASS Halofit z=1.0')
#P.plot(k, pk_class_nl2a, 'g-', lw=1.8, label='CLASS Halofit z=1.2')
P.plot(k, pk_class_nl3, 'b-', lw=1.8, label='CLASS Halofit z=1.5')
P.plot(k, pk_class_nl4, 'r-', lw=1.8, label='CLASS Halofit z=2.0')

P.plot(k, pk_ccl_nl0, 'y--', lw=1.8, label='CCL Halofit z=0.0')
P.plot(k, pk_ccl_nl1, 'c--', lw=1.8, label='CCL Halofit z=0.5')
P.plot(k, pk_ccl_nl2, 'y--', lw=1.8, label='CCL Halofit z=1.0')
#P.plot(k, pk_ccl_nl2a, 'm--', lw=1.8, label='CCL Halofit z=1.2')
P.plot(k, pk_ccl_nl3, 'c--', lw=1.8, label='CCL Halofit z=1.5')
P.plot(k, pk_ccl_nl4, 'y--', lw=1.8, label='CCL Halofit z=2.0')

P.ylabel("$P(k)$", fontsize=16)
P.xlabel("k", fontsize=16)
P.xscale('log')
P.yscale('log')
P.legend(loc='lower left', ncol=2)

P.tight_layout()
P.show()

[alexander-mead]

In your plot it looks like the CLASS P(k) curves have no Halofit correction for z=1.5 and z=2.0. They look like they are just linear theory.

[alexander-mead]

When I try to run your code I get the following error:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-4-8832e85db1d8> in <module>()
     48 # Calculate CCL power spectrum
     49 cosmo = ccl.Cosmology(**pccl)
---> 50 pk_ccl_nl0 = ccl.nonlin_matter_power(cosmo, k, 1./(1.+0.))
     51 pk_ccl_nl1 = ccl.nonlin_matter_power(cosmo, k, 1./(1.+0.5))
     52 pk_ccl_nl2 = ccl.nonlin_matter_power(cosmo, k, 1./(1.+1.))

/usr/local/lib/python2.7/site-packages/pyccl/power.pyc in nonlin_matter_power(cosmo, k, a)
     31     """
     32     return _vectorize_fn2(lib.nonlin_matter_power, 
---> 33                           lib.nonlin_matter_power_vec, cosmo, k, a)
     34 
     35 def sigmaR(cosmo, R):

/usr/local/lib/python2.7/site-packages/pyccl/pyutils.pyc in _vectorize_fn2(fn, fn_vec, cosmo, x, z, returns_status)
    185 
    186     # Check result and return
--> 187     check(status, cosmo_in)
    188     if scalar:
    189         return f[0]

/usr/local/lib/python2.7/site-packages/pyccl/pyutils.pyc in check(status, cosmo)
     22     # Check for known error status
     23     if status in pyccl.core.error_types.keys():
---> 24         raise RuntimeError("Error %s: %s" % (pyccl.core.error_types[status], msg))
     25 
     26     # Check for unknown error

RuntimeError: Error CCL_ERROR_CLASS: ccl_power.c: ccl_cosmology_compute_power_class(): Error running CLASS thermodynamics:thermodynamics_init(L:304) :error in thermodynamics_helium_from_bbn(ppr,pba,pth);
=>thermodynamics_helium_from_bbn(L:1143) :could not open fA with name /private/var/folders/vz/3cw9xc9940d1f75qgkmvd65c0000gn/T/pip-install-LFr7Cs/pyccl/build/CLASS/src/CLASS/bbn/sBBN_2017.dat and mode "r"

[philbull]

@alexander-mead Yes, that's the weird thing. It's as if halofit is just being switched off by CLASS at z >~ 1.5 or so. This happens if you use the CLASS binary or classy. But, weirdly, it doesn't happen if you call CLASS through CCL...

The error you're getting suggests that CLASS can't load a BBN-related data file that it depends on. I'm not sure why this would be happening -- are you using a version of CCL with Francois' new build system? Perhaps that might have missed these files or something.

[alexander-mead]

I did use the new build system.

[philbull]

The plot thickens. There is an automated switch inside CLASS that switches off halofit corrections if an internal parameter called kmax is too low. This is normally a silent failure, but you can see it if you set nonlinear_verbose = 1 in the CLASS .ini file. When you do that, you get the following error:

-> [WARNING:] Halofit non-linear corrections could not be computed at redshift z= 1.50 and higher.
    This is because k_max is too small for Halofit to be able to compute the scale k_NL at this redshift.
    If non-linear corrections at such high redshift really matter for you,
    just try to increase one of the parameters P_k_max_h/Mpc or P_k_max_1/Mpc or halofit_min_k_max (the code will take the max of these parameters) until reaching desired z.```

[philbull]

OK, it turns out that the Halofit corrections can be prevented from switching off by setting P_k_max_h/Mpc to a large-ish value, as recommended by the error message.

[philbull]

I think I can mark this as closed now -- the behaviour we get is expected, although confusing. CCL seems to be doing the right thing though.