GSL interpolation error when using a WDM transfer function

[robmost]

Hi,

I'm getting the following error when I run genetIC with a WDM transfer function:

gsl: interp.c:37: ERROR: insufficient number of points for interpolation type
Default GSL error handler invoked.
[1]    55914 IOT instruction (core dumped)  ./genetIC 

The genetIC parameter file I used is the following:

# Output parameters
outdir   .
outformat gadget3
outname B75_N128_WDM_genetIC_ICS-parent

# Cosmology, based on Planck 2018 (CMB + lensing + BAO)
Om     0.3111
TCMB   2.7255
Ol     0.6889
#Ob     0.0489
hubble 0.6766
s8     0.8102
ns     0.9665 
z      127

# Import transfer function data.
camb    ./LWDM_planck18+_genetIC_tk.dat

# Seed the field using the default algorithm (parallel Fourier space):
random_seed 192837465

# Specify the base-level grid. Basegrid 75 Mpc/h, 128^3
base_grid 75.0 128
gadget_particle_type 1

centre_output

done

The WDM transfer function I'm using was generated by taking a CDM transfer function from CLASS (using the CAMB output setting, so it should have the correct normalization) and then multiplying it by the Viel et al. 05 (Eq.6 here) WDM cutoff parametrization for a certain WDM particle mass. It has 7 columns, where only the first two, corresponding to the wavenumber and the DM transfer function, are non-zero.

I initially thought that my WDM transfer function has too few points for the interpolation (140 wavenumbers) but I tried with a CDM one that has the same number of values and for that setup it does work. Here are some values from the CDM transfer function:

#    1:k (h/Mpc)        2:-T_cdm/k2        3:-T_b/k2           4:-T_g/k2          5:-T_ur/k2         6:-T_ncdm/k2       7:-T_tot/k2
1.043719779506e-05  1.983212541770e+07  1.983212496349e+07  2.642898281863e+07  2.642868383264e+07  0.000000000000e+00  1.982339280717e+07
1.055805484953e-05  1.983212369680e+07  1.983212323202e+07  2.642865796929e+07  2.642835202515e+07  0.000000000000e+00  1.982339108082e+07
...
5.053402767380e+02  1.105144958511e+00  5.747035482245e-01 -6.041156819334e-14 -6.143209445162e-14  0.000000000000e+00  1.021467887905e+00

And from the WDM one:

#    1:k (h/Mpc)        2:-T_cdm/k2        3:-T_b/k2           4:-T_g/k2          5:-T_ur/k2         6:-T_ncdm/k2       7:-T_tot/k2
1.043719779506e-05 1.983212541749e+07 0.000000000000e+00 0.000000000000e+00 0.000000000000e+00 0.000000000000e+00 0.000000000000e+00
1.055805484953e-05 1.983212369658e+07 0.000000000000e+00 0.000000000000e+00 0.000000000000e+00 0.000000000000e+00 0.000000000000e+00
...
5.053402767380e+02 1.205880041635e-25 0.000000000000e+00 0.000000000000e+00 0.000000000000e+00 0.000000000000e+00 0.000000000000e+00

[apontzen]

Thanks - could you provide the full files needed to reproduce this? Otherwise it's hard to guess at what might be happening.

EDIT: at a guess, this may be because genetIC interpolates in log-log space, so is it possible some of your transfer function cut-off actually sets the power to zero? If so, it'd be -inf in the log-log space and I can see why this might cause weird errors.

[robmost]

Thanks - could you provide the full files needed to reproduce this? Otherwise it's hard to guess at what might be happening.

Sure, I created this link with the CDM and WDM transfer function and the genetIC param file.

EDIT: at a guess, this may be because genetIC interpolates in log-log space, so is it possible some of your transfer function cut-off actually sets the power to zero? If so, it'd be -inf in the log-log space and I can see why this might cause weird errors.

Makes sense... Is there a way of accounting for such transfer functions in the interpolation step? I guess I could use a transfer function up to smaller k so that it doesn't reach such low powers. This particular model uses WDM mass cutoff of 0.1 keV just so I can test it on my laptop but for a more realistic WDM mass, e.g. 6 keV, I would definitely need higher ks to probe such scales, which would probably put me back in this same situation.

[apontzen]

Ahah, so actually it turns out I already thought of this issue and there is no intrinsic problem dealing with zeros (or close to zeros) in the transfer function.

Instead what seems what is going on is the error is coming when genetIC tries to read the 'baryon' and 'all' transfer functions. These are both zero in your file, at every k. So there is no valid k at which they can be evaluated, and the GSL error results. At the moment genetIC uses the 'all' transfer function but reads all three.

I think you can just apply the same Viel suppression to the 'baryon' and 'all' transfer functions rather than only the DM column, and you should be fine?

[robmost]

That is in fact what was happening. Thank you! Now I was able to generate the ICs.

Can you clarify if genetIC uses the 'all' transfer function (i.e. 7th column) or the 'dark matter' transfer function (i.e. 2nd column) instead to generate the ICs? I thought it used only the first 3 columns (i.e. k, dark matter, baryons) but it needed all of the 7 CAMB columns (hence the zero padding), and at the moment it only supported the 'dark matter' transfer function for both dm and gas components so I didn't bother setting T_baryons.

From the genetIC manual in Sec. 4.2:

If transfer functions are available from another source, they should be prepared in the same format as a CAMB file. Specifically: space delimited columns and new-line delimited rows, with column 1 containing k values, column 2 T_DarkMatter , and column 3 containing T_Baryons . 7 or 13 columns must be present in total, but can be padded with zeros as only the first three columns are actually used by the code.

[apontzen]

Thanks for highlighting -- this is an error in the manual. If T_all is available, it will default to using that. (Of course it is extremely similar to T_DM)

[robmost]

Got it. Thank you!