Preserve CCL neutrino masses and species in CAMB predictions

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

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)