Merge pull request #1110 from LSSTDESC/baryons_vd

Baryons as in van Daalen+ 2019

CHANGELOG.md

@@ -3,6 +3,7 @@
 # v2.8.0 Changes
 
 ## Python library
+- van Daalen+ baryons model (#1110)
 - Refactor pyccl under an abstract CCLObject (#934).
 - Bump Python version to 3.8 (#1031).
 - Optimize HM autocorrelation power spectra (#891).

doc/0000-ccl_note/main.bib

@@ -1,3 +1,23 @@
+@ARTICLE{vanDaalen19,
+       author = {{van Daalen}, Marcel P. and {McCarthy}, Ian G. and {Schaye}, Joop},
+        title = "{Exploring the effects of galaxy formation on matter clustering through a library of simulation power spectra}",
+      journal = {\mnras},
+     keywords = {gravitational lensing: weak, surveys, galaxies: formation, large-scale structure of Universe, cosmology: theory, Astrophysics - Cosmology and Nongalactic Astrophysics, Astrophysics - Astrophysics of Galaxies},
+         year = 2020,
+        month = jan,
+       volume = {491},
+       number = {2},
+        pages = {2424-2446},
+          doi = {10.1093/mnras/stz3199},
+archivePrefix = {arXiv},
+       eprint = {1906.00968},
+ primaryClass = {astro-ph.CO},
+       adsurl = {https://ui.adsabs.harvard.edu/abs/2020MNRAS.491.2424V},
+      adsnote = {Provided by the SAO/NASA Astrophysics Data System}
+}
+
+
+
 @ARTICLE{isitgr1,
        author = {{Dossett}, Jason N. and {Ishak}, Mustapha and {Moldenhauer}, Jacob},
         title = "{Testing general relativity at cosmological scales: Implementation and parameter correlations}",

doc/0000-ccl_note/main.tex

@@ -518,8 +518,9 @@ \subsubsection{Cosmic emulator}
 \subsubsection{Impact of baryons}
 \label{ss:baryons}
 
-\ccl incorporates the impact of baryons on the total matter power spectrum via
-the ``baryonic correction model'' (BCM) of \citet{Schneider15}. When this
+\ccl incorporates several models of the impact of baryons on the total matter power spectrum.
+
+One of the modelling options is the ``baryonic correction model'' (BCM) of \citet{Schneider15}. When this
 correction is in use, the nonlinear matter power spectrum (whichever the
 method that provides it) is multiplied by a correction factor $F(k,z)$ which
 models the impact of baryons.
@@ -541,6 +542,14 @@ \subsubsection{Impact of baryons}
 spectrum with baryonic feedback included, \ccl will assume the default parameters
 of \citet{Schneider15}.
 
+A second option is the fitting function presented in \citet{vanDaalen19}. This model
+uses only one free parameter, which relates to the fraction of baryons in halos of mass
+$10^{14}\,h^{-1}\rm M_\odot$, normalised to $\Omega_{\rm b}/\Omega_{\rm m}$, the universal
+baryon fraction. There are two available parametrisations depending on whether the
+baryon fraction is measured within 500 or 200 times the critical density.
+Notice the model has only been validated against cosmological simulations at $z=0$
+and can be applied up to $k\leq 1\,h/\rm Mpc$ for the time being. 
+
 \subsubsection{Modified gravity ($\mu, \Sigma$)}
 
 \ccl supports the quasistatic parameterization of modified gravity with

pyccl/baryons/__init__.py

@@ -1,3 +1,4 @@
 from .baryons_base import *
 from .schneider15 import *
 from .baccoemu_baryons import *
+from .vandaalen19 import *

pyccl/baryons/vandaalen19.py

@@ -0,0 +1,120 @@
+__all__ = ("BaryonsvanDaalen19",)
+
+import numpy as np
+
+from .. import Pk2D
+from . import Baryons
+
+
+class BaryonsvanDaalen19(Baryons):
+    """The "van Daalen+ 2019" model boost factor for baryons.
+
+    .. note:: First presented in
+              `van Daalen et al., <https://arxiv.org/abs/1906.00968>`_. See the
+              `DESC Note <https://github.com/LSSTDESC/CCL/blob/master/doc\
+/0000-ccl_note/main.pdf>`_
+              for details.
+
+              The boost factor is applied multiplicatively so that
+              :math:`P_{\\rm bar.}(k, a) = P_{\\rm DMO}(k, a)\\,
+              f_{\\rm BCM}(k, a)`.
+
+              Notice the model has only been tested at z=0 and is valid for
+              :math:`k\\leq 1 \\,h/{\\rm Mpc}`.
+
+    Args:
+        fbar (:obj:`float`): the fraction of baryons in a halo within
+            an overdensity of X times the critical density, given in units
+            of the ratio of :math:`\\Omega_b` to :math:`\\Omega_m`.
+            Default to 0.7 which is approximately compatible with observations.
+            See Figure 16 of the paper.
+        mass_def (:obj:`string`): whether the mass definition corresponds to
+            X=500 or 200 critical. Options are "500c" or "200c".
+
+    """
+    name = 'vanDaalen19'
+    __repr_attrs__ = __eq_attrs__ = ("fbar", "mass_def",)
+
+    def __init__(self, fbar=0.7, mass_def='500c'):
+        self.fbar = fbar
+        self.mass_def = mass_def
+        if mass_def not in ["500c", "200c"]:
+            raise ValueError(f"Mass definition {mass_def} not supported "
+                             "for van Daalen 2019 model.")
+
+    def boost_factor(self, cosmo, k, a):
+        """The vd19 model boost factor for baryons.
+
+        Args:
+            cosmo (:class:`~pyccl.cosmology.Cosmology`): Cosmological parameters.
+            k (:obj:`float` or `array`): Wavenumber (in :math:`{\\rm Mpc}^{-1}`).
+            a (:obj:`float` or `array`): Scale factor.
+
+        Returns:
+            :obj:`float` or `array`: Correction factor to apply to
+                the power spectrum.
+
+        """ # noqa
+        a_use, k_use = map(np.atleast_1d, [a, k])
+        # [k]=1/Mpc [k_use]=h/Mpc
+        a_use, k_use = a_use[:, None], k_use[None, :]/cosmo['h']
+
+        if self.mass_def == '500c':
+            avD19 = 2.215
+            bvD19 = 0.1276
+            cvD19 = 1.309
+            dvD19 = -5.99
+            evD19 = -0.5107
+        else:
+            avD19 = 2.111
+            bvD19 = 0.0038
+            cvD19 = 1.371
+            dvD19 = -5.816
+            evD19 = -0.4005
+
+        numf = (2**avD19 +
+                2**bvD19*(cvD19*self.fbar)**(bvD19-avD19))
+        expf = np.exp(dvD19*self.fbar+evD19)
+        denf1 = (cvD19*self.fbar)**(bvD19-avD19)
+        denf = k_use**(-avD19)+k_use**(-bvD19)*denf1
+        fka = 1-numf/denf*expf
+
+        if np.ndim(k) == 0:
+            fka = np.squeeze(fka, axis=-1)
+        if np.ndim(a) == 0:
+            fka = np.squeeze(fka, axis=0)
+        return fka
+
+    def update_parameters(self, fbar=None, mass_def=None):
+        """Update van Daalen 2019 parameters. All parameters set to
+        ``None`` will be left untouched.
+
+        Args:
+            fbar (:obj:`float`): baryonic fraction in halos
+                within X times the critical density.
+
+            mass_def (:obj:`string`): mass definition: whether 500 ("500c")
+                or 200 critical ("200c").
+        """
+        if fbar is not None:
+            self.fbar = fbar
+        if mass_def is not None:
+            self.mass_def = mass_def
+            if mass_def not in ["500c", "200c"]:
+                raise ValueError(f"Mass definition {mass_def} not supported "
+                                 "for van Daalen 2019 model.")
+
+    def _include_baryonic_effects(self, cosmo, pk):
+        # Applies boost factor
+        a_arr, lk_arr, pk_arr = pk.get_spline_arrays()
+        k_arr = np.exp(lk_arr)
+        fka = self.boost_factor(cosmo, k_arr, a_arr)
+        pk_arr *= fka
+
+        if pk.psp.is_log:
+            np.log(pk_arr, out=pk_arr)  # in-place log
+
+        return Pk2D(a_arr=a_arr, lk_arr=lk_arr, pk_arr=pk_arr,
+                    is_logp=pk.psp.is_log,
+                    extrap_order_lok=pk.extrap_order_lok,
+                    extrap_order_hik=pk.extrap_order_hik)

pyccl/tests/test_baryons_vd19.py

@@ -0,0 +1,66 @@
+import numpy as np
+import pyccl as ccl
+import pytest
+
+# Set tolerances
+BOOST_TOLERANCE = 1e-5
+
+# Set up the cosmological parameters to be used
+cosmo = ccl.Cosmology(Omega_c=0.27, Omega_b=0.045, h=0.67,
+                      A_s=2.1e-9, n_s=0.96,
+                      transfer_function='boltzmann_class')
+k = 0.5*cosmo['h']
+a = 1.
+
+# Set up power without baryons
+pk2D_no_baryons = cosmo.get_nonlin_power()
+
+fbarcvec = np.linspace(0.25, 1, 20)
+mdef_vec = ['500c', '200c']
+
+
+def compare_boost():
+    vdboost = []
+    vdboost_expect = []
+    pk_nl = pk2D_no_baryons(k, a)
+    for mdef in mdef_vec:
+        for f in fbarcvec:
+            # Takes ftilde as argument
+            vD19 = ccl.BaryonsvanDaalen19(fbar=f, mass_def=mdef)
+            pk2D_with_baryons = vD19.include_baryonic_effects(
+                cosmo, pk2D_no_baryons)
+            # Takes k in units of 1/Mpc as argument
+            pk_nl_bar = pk2D_with_baryons(k, a)
+            vdboost.append(pk_nl_bar/pk_nl-1)
+            if mdef == '500c':
+                vdboost_expect.append(-np.exp(-5.99*f-0.5107))
+            else:
+                vdboost_expect.append(-np.exp(-5.816*f-0.4005))
+    assert np.allclose(
+        vdboost, vdboost_expect, atol=1e-5, rtol=BOOST_TOLERANCE)
+
+
+def test_boost_model():
+    compare_boost()
+
+
+def test_baryons_from_name():
+    baryons = ccl.BaryonsvanDaalen19(fbar=0.7, mass_def='500c')
+    bar2 = ccl.Baryons.from_name('vanDaalen19')
+    assert baryons.name == bar2.name
+    assert baryons.name == 'vanDaalen19'
+
+
+def test_baryons_vd19_raises():
+    with pytest.raises(ValueError):
+        ccl.BaryonsvanDaalen19(fbar=0.7, mass_def='blah')
+    b = ccl.BaryonsvanDaalen19(fbar=0.7, mass_def='500c')
+    with pytest.raises(ValueError):
+        b.update_parameters(mass_def='blah')
+
+
+def test_update_params():
+    b = ccl.BaryonsvanDaalen19(fbar=0.7, mass_def='500c')
+    b.update_parameters(fbar=0.6, mass_def='200c')
+    assert b.mass_def == '200c'
+    assert b.fbar == 0.6

readthedocs/api/pyccl.baryons.rst

@@ -10,6 +10,7 @@ Submodules
    pyccl.baryons.baccoemu_baryons
    pyccl.baryons.baryons_base
    pyccl.baryons.schneider15
+   pyccl.baryons.vandaalen19
 
 Module contents
 ---------------