DeaglanBartlett
commented
2 months ago The first thing to note here is that this does not appear to be an issue with BORG's implementation/wrapper of CLASS, but with the CLASS code itself. Using this cosmology using CLASS's python wrapper (v3.2.0):
from classy import Class
import numpy as np
import camb
redshift = 0.0
tau = 0.0561
Om, Ob, h, ns, sigma8 = 0.4791, 0.06759, 0.8131, 0.9573, 0.7391
kmin = 9e-3
kmax = 4.9
nk = 200
k = np.logspace(np.log10(kmin), np.log10(kmax), nk)
# Convert As to sigma8 (at z=0)
As = 2.1e-09 # initial guess
pars = camb.CAMBparams()
pars.set_cosmology(H0=h*100,
ombh2=Ob * h ** 2,
omch2=(Om - Ob) * h ** 2,
mnu=0.0,
omk=0,
tau=tau,)
pars.InitPower.set_params(As=As, ns=ns, r=0)
pars.set_matter_power(redshifts=[0.], kmax=k[-1])
pars.NonLinear = camb.model.NonLinear_none
results = camb.get_results(pars)
sigma8_new = results.get_sigma8()[0]
As *= (sigma8 / sigma8_new) ** 2
print('Using As:', As)
# Run camb
num_massive_neutrinos = 0
Oc = Om - Ob
pars = camb.CAMBparams()
pars.set_cosmology(H0=h*100,
ombh2=Ob * h ** 2,
omch2=(Om - Ob) * h ** 2,
mnu=0.0,
omk=0,
tau=tau,)
pars.set_dark_energy(w=-1.0, wa=0.0)
pars.InitPower.set_params(ns=ns, As=As, r=0)
pars.set_matter_power(redshifts=[redshift], kmax=k[-1])
pars.NonLinear = camb.model.NonLinear_none
results = camb.get_results(pars)
index = 6
_, _, plin_camb = results.get_matter_power_spectrum(
var1=(1 + index),
var2=(1 + index),
minkh=k.min(),
maxkh=k.max(),
npoints=len(k)
)
plin_camb = plin_camb[0]
print('CAMB sigma8', results.get_sigma8()[0])
# Run class
Oc = Om - Ob
class_params = {
'h': h,
'omega_b': Ob * h**2,
'omega_cdm': Oc * h**2,
'A_s': As,
'n_s': ns,
'output': 'mPk',
'P_k_max_1/Mpc': k.max() * h,
'w0_fld': -1.0,
'wa_fld': 0.0,
'Omega_Lambda': 0, # Set to 0 because we're using w0_fld and wa_fld instead
'tau_reio':tau,
'z_max_pk': 3.0, # Max redshift for P(k) output
}
cosmo = Class()
cosmo.set(class_params)
cosmo.compute()
sigma8_class = cosmo.sigma8()
plin_class = np.array([cosmo.pk_lin(kk*h, redshift) for kk in k]) * h ** 3
# Get sigma8 from class
sigma8_class = cosmo.sigma8()
print('CLASS sigma8:', sigma8_class)
# Memory cleanup for class
cosmo.struct_cleanup()
cosmo.empty()also gives an error at the line cosmo.compute():
File classy.pyx:389, in classy.Class.compute()
CosmoComputationError:
Error in Class: thermodynamics_init(L:344) :error in thermodynamics_helium_from_bbn(ppr,pba,pth);
=>thermodynamics_helium_from_bbn(L:615) :condition (omega_b > omegab[num_omegab-1]) is true; You have asked for an unrealistic high value omega_b = 4.468589e-02. The corresponding value of the primordial helium fraction cannot be found in the interpolation table. If you really want this value, you should fix YHe to a given value rather than to BBNI will investigate whether it is possible to follow the advice If you really want this value, you should fix YHe to a given value rather than to BBN.
ludvigdoeser
maho3
The CLASS implementation in BORG seems to break for high values of Omega_b.
For example, running
lhid=15at the cosmology:[Omega_m, Omega_b, h, n_s, sigma8] = [0.4791 0.06759 0.8131 0.9573 0.7391]results in...
Note, it complained about
omega_b = 4.468589e-02which is different than the config cosmology. I assume this is the omega_b at z=127, or whenever the transfer function is applied.This seems to break for ~40% of Quijote cosmologies.
Can we fix this, or will we need to rely on the CAMB implementation in FastPM to run reliable PM sims?
EDIT: I am also getting a different error with BORG+CLASS, which I think may be related to this^ problem. For example, for lhid=408, i.e.
[Omega_m, Omega_b, h, n_s, sigma8] = [0.1219, 0.05495, 0.8649, 1.0065, 0.6183]I get...