diff --git a/cs_util/cosmo.py b/cs_util/cosmo.py index f777541..62673db 100644 --- a/cs_util/cosmo.py +++ b/cs_util/cosmo.py @@ -371,13 +371,37 @@ def _ccl_to_camb(cosmo): """ h = cosmo["h"] + masses = np.atleast_1d(cosmo["m_nu"]) + masses = masses[masses > 0] + mnu = float(np.sum(masses)) + n_massive = len(masses) camb_params = { "H0": h * 100, "ombh2": cosmo["Omega_b"] * h**2, "omch2": cosmo["Omega_c"] * h**2, "ns": cosmo["n_s"], + "mnu": mnu, + "num_massive_neutrinos": n_massive, + "nnu": cosmo["Neff"], + "TCMB": cosmo["T_CMB"], + "omk": cosmo["Omega_k"], + "w": cosmo["w0"], + "wa": cosmo["wa"], } + # Match CCL's resolved species (including explicit mass arrays), rather + # than CAMB's default single 0.06-eV species. As in pyccl.boltzmann, + # degeneracies encode T_ncdm relative to the standard neutrino temperature. + degeneracy = (cosmo["T_ncdm"] / (4 / 11) ** (1 / 3)) ** 4 + camb_params.update( + omnuh2=cosmo["Omega_nu_mass"] * h**2, + num_nu_massless=cosmo["N_nu_rel"], + nu_mass_eigenstates=n_massive, + nu_mass_numbers=np.ones(n_massive, dtype=int), + nu_mass_fractions=masses / mnu if n_massive else [], + nu_mass_degeneracies=np.full(n_massive, degeneracy), + ) + # Handle normalization: prefer As, but convert sigma8 to As if needed As_val = cosmo.__getitem__("A_s") sigma8_val = cosmo.__getitem__("sigma8") @@ -422,11 +446,6 @@ def _ccl_to_camb(cosmo): # No normalization specified, use CAMB default pass - # Add dark energy parameters if they exist - for camb_key, cosmo_key in [("w", "w0"), ("wa", "wa")]: - if hasattr(cosmo._params, cosmo_key): - camb_params[camb_key] = cosmo[cosmo_key] - return camb_params @@ -711,10 +730,7 @@ def get_theo_c_ell( NonLinear=camb.model.NonLinear_both, ) - # Adjust for neutrino contribution - if "mnu" in camb_kwargs and camb_kwargs["mnu"] > 0: - omch2_adj = camb_kwargs["omch2"] - pars.omeganu * (pars.H0 / 100) ** 2 - pars.set_cosmology(omch2=omch2_adj) + # CCL's Omega_c and CAMB's omch2 both exclude massive neutrinos. # Set up lensing source window. CAMB's min_l is the scalar-C_ell floor # (1 or 2 only), so leave it at its default; ell below ell.min() are diff --git a/cs_util/tests/test_camb_neutrinos.py b/cs_util/tests/test_camb_neutrinos.py new file mode 100644 index 0000000..99a7211 --- /dev/null +++ b/cs_util/tests/test_camb_neutrinos.py @@ -0,0 +1,121 @@ +"""Regression tests for the CCL-to-CAMB cosmology bridge.""" + +import camb +import numpy as np +import numpy.testing as npt +import pyccl as ccl +import pytest + +from cs_util import cosmo + + +@pytest.mark.parametrize( + "masses, split", + [ + (0.0, "normal"), + (0.06, "normal"), + (0.3, "inverted"), + (0.3, "equal"), + ([0.01, 0.03, 0.08], "list"), + ], +) +def test_ccl_to_camb_species(masses, split): + """Preserve resolved masses, temperatures and CDM, not CAMB defaults.""" + c = ccl.Cosmology( + Omega_c=0.25, + Omega_b=0.05, + h=0.7, + n_s=0.96, + A_s=2.1e-9, + m_nu=masses, + mass_split=split, + Neff=3.2, + T_CMB=2.73, + Omega_k=0.01, + w0=-0.9, + wa=0.1, + ) + kwargs = cosmo._ccl_to_camb(c) + p = camb.set_params(**kwargs) + assert kwargs["mnu"] == pytest.approx(np.sum(c["m_nu"])) + assert p.num_nu_massive == c["N_nu_mass"] + assert p.omch2 == pytest.approx(c["Omega_c"] * c["h"] ** 2) + assert p.ombh2 == pytest.approx(c["Omega_b"] * c["h"] ** 2) + assert p.omnuh2 == pytest.approx(c["Omega_nu_mass"] * c["h"] ** 2) + assert p.num_nu_massless == pytest.approx(c["N_nu_rel"]) + assert p.N_eff == pytest.approx(c["Neff"]) + assert p.TCMB == c["T_CMB"] + assert p.omk == c["Omega_k"] + assert p.DarkEnergy.w == c["w0"] + assert p.DarkEnergy.wa == c["wa"] + if kwargs["mnu"]: + npt.assert_allclose(p.nu_mass_fractions, np.array(c["m_nu"]) / kwargs["mnu"]) + + +def test_ccl_to_camb_scalar_mass(): + """Also accept a scalar mass representation.""" + c = ccl.Cosmology( + Omega_c=0.25, + Omega_b=0.05, + h=0.7, + n_s=0.96, + A_s=2.1e-9, + m_nu=[0.3], + mass_split="list", + ) + keys = ( + "h", + "Omega_b", + "Omega_c", + "n_s", + "Neff", + "T_CMB", + "Omega_k", + "w0", + "wa", + "T_ncdm", + "Omega_nu_mass", + "N_nu_rel", + "A_s", + "sigma8", + ) + params = {key: c[key] for key in keys} + params["m_nu"] = 0.3 + kwargs = cosmo._ccl_to_camb(params) + assert kwargs["mnu"] == pytest.approx(0.3) + assert kwargs["num_massive_neutrinos"] == 1 + + +def test_ccl_to_camb_sigma8_background(): + """Normalize on the requested neutrino and dark-energy background.""" + c = ccl.Cosmology( + Omega_c=0.25, + Omega_b=0.05, + h=0.7, + n_s=0.96, + sigma8=0.8, + m_nu=0.3, + w0=-0.9, + wa=0.1, + ) + p = camb.set_params(**cosmo._ccl_to_camb(c)) + p.set_matter_power(redshifts=[0], kmax=2) + assert camb.get_results(p).get_sigma8_0() == pytest.approx(0.8, rel=1e-3) + + +def test_camb_neutrino_response(): + """The mass response agrees even though the nonlinear models differ.""" + ell = np.array([100, 500, 1000]) + z = np.linspace(0.01, 3, 300) + nz = np.exp(-0.5 * ((z - 0.7) / 0.2) ** 2) + spectra = {} + for mass in (0.06, 0.3): + c = cosmo.get_cosmo(mnu=mass) + spectra[mass] = { + backend: cosmo.get_theo_c_ell(ell, z, nz, backend=backend, cosmo=c)["W1xW1"] + for backend in ("camb", "ccl") + } + response_camb = spectra[0.3]["camb"] / spectra[0.06]["camb"] + response_ccl = spectra[0.3]["ccl"] / spectra[0.06]["ccl"] + assert abs(response_camb[-1] - 1) > 0.005 + npt.assert_allclose(response_camb, response_ccl, rtol=0.02)