New full_limber scheme. Also, cubic interpolation of lnPk, and minor …

…bug fixes in fourier_pk_at_z and in get_pk_all()
miguelzuma · Feb 3, 2024 · a7304a5 · a7304a5
2 parents 997d1ac + 8df566c
commit a7304a5
Show file tree

Hide file tree

Showing 12 changed files with 738 additions and 308 deletions.
diff --git a/explanatory.ini b/explanatory.ini
@@ -1151,7 +1151,17 @@ A_L =
 #lcmb_tilt = 0
 #lcmb_pivot = 0.1
 
-
+# In general, do we want to use the full Limber scheme introduced in
+# v3.2.2? With this full Limber scheme, the calculation of the CMB
+# lensing potential spectrum C_l^phiphi for l > ppr->l_switch_limber
+# is based on a new integration scheme. Compared to the previous
+# scheme, which can be recovered by setting this parameter to 'no',
+# the new scheme uses a larger k_max and a coarser k-grid (or q-grid)
+# than the CMB transfer function. The new scheme is used by default,
+# because the old one is inaccurate at large l due to the too small
+# k_max. Set to anything starting with the letter 'y' or
+# 'n'. (default: 'yes')
+want_lcmb_full_limber = yes
 
 # -------------------------------------
 # ----> Distortions parameters:

diff --git a/include/common.h b/include/common.h
@@ -15,7 +15,7 @@
 #ifndef __COMMON__
 #define __COMMON__
 
-#define _VERSION_ "v3.2.1"
+#define _VERSION_ "v3.2.2"
 
 /* @cond INCLUDE_WITH_DOXYGEN */
 

diff --git a/include/harmonic.h b/include/harmonic.h
@@ -181,6 +181,7 @@ extern "C" {
                        );
 
   int harmonic_cls(
+                   struct precision * ppr,
                    struct background * pba,
                    struct perturbations * ppt,
                    struct transfer * ptr,
@@ -189,6 +190,7 @@ extern "C" {
                    );
 
   int harmonic_compute_cl(
+                          struct precision * ppr,
                           struct background * pba,
                           struct perturbations * ppt,
                           struct transfer * ptr,
@@ -200,6 +202,7 @@ extern "C" {
                           int index_l,
                           int cl_integrand_num_columns,
                           double * cl_integrand,
+                          double * cl_integrand_limber,
                           double * primordial_pk,
                           double * transfer_ic1,
                           double * transfer_ic2

diff --git a/include/perturbations.h b/include/perturbations.h
@@ -150,6 +150,8 @@ struct perturbations
   int l_lss_max; /**< maximum l value for LSS \f$ C_l \f$'s (density and lensing potential in  bins) */
   double k_max_for_pk; /**< maximum value of k in 1/Mpc required for the output of P(k,z) and T(k,z) */
 
+  short want_lcmb_full_limber; /**< In general, do we want to use the full Limber scheme introduced in v3.2.2? With this full Limber scheme, the calculation of the CMB lensing potential spectrum C_l^phiphi for l > ppr->l_switch_limber is based on a new integration scheme. Compared to the previous scheme, which can be recovered by switching this parameter to _FALSE_, the new scheme uses a larger k_max and a coarser k-grid (or q-grid) than the CMB transfer function. The new scheme is used by default, because the old one is inaccurate at large l due to the too small k_max. */
+
   int selection_num;                            /**< number of selection functions
                                                    (i.e. bins) for matter density \f$ C_l \f$'s */
   enum selection_type selection;                /**< type of selection functions */

diff --git a/include/precisions.h b/include/precisions.h
@@ -247,7 +247,7 @@ class_precision_parameter(z_start_chi_approx,double,2.0e3) /**< Switching redshi
 
 class_precision_parameter(k_min_tau0,double,0.1) /**< number defining k_min for the computation of Cl's and P(k)'s (dimensionless): (k_min tau_0), usually chosen much smaller than one */
 
-class_precision_parameter(k_max_tau0_over_l_max,double,2.4) /**< number defining k_max for the computation of Cl's (dimensionless): (k_max tau_0)/l_max, usually chosen around two */
+class_precision_parameter(k_max_tau0_over_l_max,double,1.8) /**< number defining k_max for the computation of Cl's (dimensionless): (k_max tau_0)/l_max, usually chosen around two. In v3.2.2: lowered from 2.4 to 1.8, because the high value 2.4 was needed to keep CMB lensing accurate enough. With the new Limber scheme, this will be the case anyway, and k_max can be lowered in other observables in order to speed up the code. */
 class_precision_parameter(k_step_sub,double,0.05) /**< step in k space, in units of one period of acoustic oscillation at decoupling, for scales inside sound horizon at decoupling */
 class_precision_parameter(k_step_super,double,0.002) /**< step in k space, in units of one period of acoustic oscillation at decoupling, for scales above sound horizon at decoupling */
 class_precision_parameter(k_step_transition,double,0.2) /**< dimensionless number regulating the transition from 'sub' steps to 'super' steps. Decrease for more precision. */
@@ -325,6 +325,15 @@ class_precision_parameter(perturbations_integration_stepsize,double,0.5)
  * default step \f$ d \tau \f$ for sampling the source function, in units of the timescale involved in the sources: \f$ (\dot{\kappa}- \ddot{\kappa}/\dot{\kappa})^{-1} \f$
  */
 class_precision_parameter(perturbations_sampling_stepsize,double,0.1)
+/**
+ * added in v 3.2.2: age fraction (between 0 and 1 ) such that, when
+ * tau > conformal_age * age_fraction, the time sampling of sources is
+ * twice finer, in order to boost the accuracy of the lensing
+ * line-of-sight integrals (for l < l_switch_limber) without changing
+ * that of unlensed CMB observables. Setting to 1.0 disables this
+ * functionality.
+*/
+class_precision_parameter(perturbations_sampling_boost_above_age_fraction, double, 0.9)
 /**
  * control parameter for the precision of the perturbation integration,
  * IMPORTANT FOR SETTING THE STEPSIZE OF NDF15
@@ -461,6 +470,9 @@ class_precision_parameter(q_numstep_transition,double,250.0) /**< number of step
                                  q_logstep_spline steps (transition
                                  must be smooth for spline) */
 
+class_precision_parameter(q_logstep_limber,double,1.025) /**< new in v3.2.2: in the new 'full limber' scheme, logarithmic step for the k-grid (and q-grid) */
+class_precision_parameter(k_max_limber_over_l_max_scalars,double,0.001) /**< new in v3.2.2: in the new 'full limber' scheme, the integral runs up to k_max = l_max_scalars times this parameter (units of 1/Mpc) */
+
 class_precision_parameter(transfer_neglect_delta_k_S_t0,double,0.15) /**< for temperature source function T0 of scalar mode, range of k values (in 1/Mpc) taken into account in transfer function: for l < (k-delta_k)*tau0, ie for k > (l/tau0 + delta_k), the transfer function is set to zero */
 class_precision_parameter(transfer_neglect_delta_k_S_t1,double,0.04) /**< same for temperature source function T1 of scalar mode */
 class_precision_parameter(transfer_neglect_delta_k_S_t2,double,0.15) /**< same for temperature source function T2 of scalar mode */

diff --git a/include/transfer.h b/include/transfer.h
@@ -171,14 +171,22 @@ struct transfer {
 
   //@{
 
-  size_t q_size; /**< number of wavenumber values */
+  size_t q_size; /**< number of wavenumber values corresponding to k up to k_max_cl */
 
   double * q;  /**< list of wavenumber values, q[index_q] */
 
   double ** k; /**< list of wavenumber values for each requested mode, k[index_md][index_q]. In flat universes k=q. In non-flat universes q and k differ through q2 = k2 + K(1+m), where m=0,1,2 for scalar, vector, tensor. q should be used throughout the transfer module, excepted when interpolating or manipulating the source functions S(k,tau): for a given value of q this should be done in k(q). */
 
   int index_q_flat_approximation; /**< index of the first q value using the flat rescaling approximation */
 
+  short do_lcmb_full_limber; /**< in this particular run, will we use the full Limber scheme? */
+
+  size_t q_size_limber; /**< number of wavenumber values corresponding to k up to k_max */
+
+  double * q_limber;  /**< list of wavenumber values used in full limber scheme, q_limber[index_q] */
+
+  double ** k_limber; /**< list of wavenumber values used in full limber scheme */
+
   //@}
 
   /** @name - transfer functions */
@@ -187,6 +195,8 @@ struct transfer {
 
   double ** transfer; /**< table of transfer functions for each mode, initial condition, type, multipole and wavenumber, with argument transfer[index_md][((index_ic * ptr->tt_size[index_md] + index_tt) * ptr->l_size[index_md] + index_l) * ptr->q_size + index_q] */
 
+  double ** transfer_limber; /**< table of transfer functions used in full limber scheme */
+
   //@}
 
   /** @name - technical parameters */
@@ -367,14 +377,13 @@ extern "C" {
                           int sgnK
                           );
 
-  int transfer_get_q_list_v1(
-                             struct precision * ppr,
-                             struct perturbations * ppt,
-                             struct transfer * ptr,
-                             double q_period,
-                             double K,
-                             int sgnK
-                             );
+  int transfer_get_q_limber_list(
+                                 struct precision * ppr,
+                                 struct perturbations * ppt,
+                                 struct transfer * ptr,
+                                 double K,
+                                 int sgnK
+                                 );
 
   int transfer_get_k_list(
                           struct perturbations * ppt,
@@ -427,7 +436,8 @@ extern "C" {
                                   double *** sources,
                                   double *** sources_spline,
                                   double * window,
-                                  struct transfer_workspace * ptw
+                                  struct transfer_workspace * ptw,
+                                  short use_full_limber
                                   );
 
   int transfer_radial_coordinates(
@@ -440,7 +450,7 @@ extern "C" {
   int transfer_interpolate_sources(
                                    struct perturbations * ppt,
                                    struct transfer * ptr,
-                                   int index_q,
+                                   double k,
                                    int index_md,
                                    int index_ic,
                                    int index_type,
@@ -456,7 +466,7 @@ extern "C" {
                        struct transfer * ptr,
                        double * interpolated_sources,
                        double tau_rec,
-                       int index_q,
+                       double k,
                        int index_md,
                        int index_tt,
                        double * sources,
@@ -548,7 +558,8 @@ extern "C" {
                                   int index_l,
                                   double l,
                                   double q_max_bessel,
-                                  radial_function_type radial_type
+                                  radial_function_type radial_type,
+                                  short use_full_limber
                                   );
 
   int transfer_use_limber(

diff --git a/python/classy.pyx b/python/classy.pyx
@@ -894,7 +894,7 @@ cdef class Class:
                     pk_cb[index_k,index_z,index_mu] = self.pk_cb_lin(k[index_k,index_z,index_mu],z[index_z])
         return pk_cb
 
-    def get_pk_all(self, k, z, nonlinear = True, cdmbar = False, z_axis_in_k_arr = 0):
+    def get_pk_all(self, k, z, nonlinear = True, cdmbar = False, z_axis_in_k_arr = 0, interpolation_kind='cubic'):
         """ General function to get the P(k,z) for ARBITRARY shapes of k,z
             Additionally, it includes the functionality of selecting wether to use the non-linear parts or not,
             and wether to use the cdm baryon power spectrum only
@@ -929,16 +929,16 @@ cdef class Class:
 
         # Only get the nonlinear function where the nonlinear treatment is possible
         def _islinear(z):
-          if z > z_max_nonlinear:
+          if z > z_max_nonlinear or (self.fo.method == nl_none):
             return True
           else:
             return False
 
         # A simple wrapper for writing the P(k) in the given location and interpolating it
         def _interpolate_pk_at_z(karr,z):
           _write_pk(z,_islinear(z),ispkcb)
-          interp_func = interp1d(k_out,pk_out,kind='linear',copy=True)
-          return interp_func(karr)
+          interp_func = interp1d(k_out,np.log(pk_out),kind=interpolation_kind,copy=True)
+          return np.exp(interp_func(karr))
 
         # 2) Check if z array, or z value
         if not isinstance(z,(list,np.ndarray)):

diff --git a/source/fourier.c b/source/fourier.c
@@ -86,6 +86,8 @@ int fourier_pk_at_z(
   if ((pk_output == pk_linear) && (pfo->ic_size > 1) && (out_pk_ic != NULL))
     do_ic = _TRUE_;
 
+  class_test(pk_output == pk_nonlinear && pfo->method == nl_none, pfo->error_message, "Cannot get nonlinear power spectrum when no nonlinear method is employed");
+
   /** - case z=0 requiring no interpolation in z */
   if (z == 0) {