lesgourg / class_public

Public repository of the Cosmic Linear Anisotropy Solving System (master for the most recent version of the standard code; GW_CLASS to include Cosmic Gravitational Wave Background anisotropies; classnet branch for acceleration with neutral networks; ExoCLASS branch for exotic energy injection; class_matter branch for FFTlog)
220 stars 291 forks source link

HMCode for Decaying dark matter #492

Open ShazAlvi opened 2 years ago

ShazAlvi commented 2 years ago

Hi,

I tried to use the HMCode non-linear model to produce matter power spectra in the decaying dark matter model. I used the following code to produce the spectra.

from classy import Class
import matplotlib.pyplot as plt
import numpy as np

k = np.logspace(-3,1,1000)
Pk_LCDM = np.zeros((len(k)))
Pk_LDDM = np.zeros((len(k)))
params = {'output': 'mPk','omega_b': 0.02237,'omega_cdm':0.1201,'ln10^{10}A_s': 3.054, 'n_s': 0.9651,
                '100*theta_s': 1.04090, 'tau_reio': 0.0543, 'non linear': 'HMCODE','P_k_max_1/Mpc': 15.47,
                'z_max_pk':2.5}

InsCLASS = Class()
InsCLASS.set(params)
InsCLASS.compute()
for i, k_i in enumerate(k):
    Pk_LCDM[i] = InsCLASS.pk(k_i,z=0)
InsCLASS.empty()
InsCLASS.struct_cleanup()

params = {'output': 'mPk','omega_b': 0.02237,'omega_ini_dcdm':0.1201,'omega_cdm':0.0, 'Gamma_dcdm': 294.12,'ln10^{10}A_s': 3.054, 'n_s': 0.9651,
                '100*theta_s': 1.04090, 'tau_reio': 0.0543, 'non linear': 'HMCODE','P_k_max_1/Mpc': 15.47,'z_max_pk':2.5}

for i, k_i in enumerate(k):
    Pk_LDDM[i] = InsCLASS.pk(k_i,z=0)
InsCLASS = Class()
InsCLASS.set(params)
InsCLASS.compute()

InsCLASS.empty()
InsCLASS.struct_cleanup()

plt.semilogx(k, (1-(Pk_LDDM/Pk_LCDM))*100)
plt.show()

While the code works for the LCDM model, in the decaying dark matter model CLASS gives different errors (even if the same code is run at two different times without any changes). Sometimes the error is,

Traceback (most recent call last):
  File "DCDM_pk_HMCode.py", line 24, in <module>
    Pk_LDDM[i] = InsCLASS.pk(k_i,z=0)
  File "classy.pyx", line 775, in classy.Class.pk
classy.CosmoSevereError: 

Error in Class: nonlinear_pk_at_k_and_z(L:413) :condition ((k < 0.) || (k > exp(pnl->ln_k[pnl->k_size-1]))) is true; k=1.006939e+00 out of bounds [0.000000e+00:1.000000e+00]

Other times there is a "segmentation fault (core dumped)" error or "free(): invalid pointer Aborted (core dumped)".

Other times there is no error and the plot is produced, however, the percentage difference between the two power spectra is close to zero. I suspect that the HMCode is not written for a decaying dark matter model (?). It is mentioned in the code that for warm dark matter or for interacting dark matter the HMCode is not reliable. Perhaps due to a similar reason, it has not been programmed for a decaying dark matter model (?).

Any feedback would helpful! Thanks a bunch!