Merge remote-tracking branch 'class_public/master' into gravity_dream

cl_ref.pre

@@ -27,8 +27,14 @@ l_max_pol_g=25
 l_max_ur=50
 l_max_ncdm=50
 
-tol_perturb_integration=1.e-6
-perturb_sampling_stepsize=0.01
+#Maximum multipoles in the tensor hierarchy
+#This requires hierarchy = tam in the .ini file, otherwise the precision is too stiff for the optimal hierarchy.
+l_max_g_ten=50
+l_max_pol_g_ten=25
+l_max_ur_ten=50
+
+tol_perturbations_integration=1.e-6
+perturbations_sampling_stepsize=0.01
 
 radiation_streaming_approximation = 2
 radiation_streaming_trigger_tau_over_tau_k = 240.

explanatory.ini

@@ -151,14 +151,14 @@ gauge = synchronous
 #      (default: set to 'no')
 #nbody_gauge_transfer_functions = yes
 
-# 5) Hubble parameter : either 'H0' in km/s/Mpc or 'h' or '100*theta_s', where
+# 5) Hubble parameter : either 'H0' in km/s/Mpc or 'h' (or 'theta_s_100'), where
 #    the latter is the peak scale parameter defined exactly as 100(ds_dec/da_dec)
 #    with a decoupling time given by maximum of visibility function (quite different
 #    from theta_MC of CosmoMC and slightly different from theta_* of CAMB)
 #    (default: 'h' set to 0.67810 such that 100*theta_s = 1.041783 like in Planck 2018)
 h = 0.67810
 #H0 = 67.810
-#100*theta_s = 1.041783
+#theta_s_100 = 1.041783
 
 
 # 6) Primordial Helium fraction 'YHe', e.g. 0.25; if set to 'BBN' or 'bbn',
@@ -725,7 +725,7 @@ PBH_accretion_eigenvalue = 0.1
 #      user can pass explicitely the value of f_eff. (default: f_eff=1)
 #
 #    - for 'from_file' the code reads a precomputed function in an external file
-#      with a path set by the user (default set to "/external/heating/example_f_eff_file.dat")
+#      with a path set by the user (default set to "external/heating/example_f_eff_file.dat")
 f_eff_type = on_the_spot
 #f_eff =
 #f_eff_file = external/heating/example_f_eff_file.dat
@@ -751,9 +751,9 @@ chi_type = CK_2004
 #      Two files 'example_chix_file.dat' and 'example_chiz_file.dat' are given
 #      as example in external/heating. Note that 'example_chix_file.dat' has
 #      been computed according to the approximations of Galli et al. 2013.
-#      (default set to "/external/heating/example_f_eff_file.dat")
-#chi_z_file = external/heating/example_chiz_file.dat;
-#chi_x_file = external/heating/example_chix_file.dat;
+#      (default set to "external/heating/example_f_eff_file.dat")
+#chi_file = external/heating/example_chiz_file.dat
+#chi_file = external/heating/example_chix_file.dat
 
 
 
@@ -831,11 +831,12 @@ k_pivot = 0.05
 # 1.b) For type 'analytic_Pk':
 # 1.b.1) For scalar perturbations
 #        curvature power spectrum value at pivot scale ('A_s' or
-#        'ln10^{10}A_s') OR 'sigma8' (found by iterations using a shooting
+#        'ln_A_s_1e10') OR one of 'sigma8' or 'S8' (found by iterations using a shooting
 #        method). (default: set 'A_s' to 2.100549e-09)
 A_s = 2.100549e-09
-#ln10^{10}A_s = 3.04478383
-# sigma8 = 0.824398
+#ln_A_s_1e10 = 3.04478383
+#sigma8 = 0.824398
+#S8 = 0.837868
 
 # 1.b.1.1) Adiabatic perturbations:
 #          tilt at the same scale 'n_s', and tilt running 'alpha_s'
@@ -1150,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:

include/common.h

@@ -15,7 +15,7 @@
 #ifndef __COMMON__
 #define __COMMON__
 
-#define _VERSION_ "v3.2.0"
+#define _VERSION_ "v3.2.2"
 
 /* @cond INCLUDE_WITH_DOXYGEN */
 

include/fourier.h

@@ -105,7 +105,6 @@ struct fourier {
                           condition z(tau) < z_max_pk */
 
   int ln_tau_size;     /**< total number of values in this array */
-  int index_ln_tau_pk; /**< first index relevant for output of P(k,z) and T(k,z) */
 
   double ** ln_pk_ic_l;   /**< Matter power spectrum (linear).
                              Depends on indices index_pk, index_ic1_ic2, index_k, index_tau as:

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

include/input.h

@@ -203,9 +203,9 @@
  * For shooting method: definition of the possible targets
  */
 
- enum target_names {theta_s, Omega_dcdmdr, omega_dcdmdr, Omega_scf, Omega_ini_dcdm, omega_ini_dcdm, Omega_smg, M2_today_smg, sigma8};
+ enum target_names {theta_s, theta_s_100, Omega_dcdmdr, omega_dcdmdr, Omega_scf, Omega_ini_dcdm, omega_ini_dcdm, Omega_smg, M2_today_smg, sigma8, S8};
  /* Important: Keep this number equal to the number of target_names (except sigma8) */
- #define _NUM_TARGETS_ 8
+ #define _NUM_TARGETS_ 9
 /* Important: add one for each new target_names */
 enum computation_stage {cs_background, cs_thermodynamics, cs_perturbations, cs_primordial, cs_nonlinear, cs_transfer, cs_spectra};
 

include/perturbations.h

@@ -170,6 +170,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 */
@@ -414,7 +416,6 @@ struct perturbations
                                final time range required for the output of
                                Fourier transfer functions (used for interpolations) */
   int ln_tau_size;         /**< total number of values in this array */
-  int index_ln_tau_pk;     /**< first index relevant for output of P(k,z) and T(k,z) */
 
   double *** late_sources; /**< Pointer towards the source interpolation table
                               late_sources[index_md]

include/precisions.h

@@ -12,7 +12,7 @@ class_precision_parameter(a_ini_over_a_today_default,double,1.e-14)
 /**
  * Number of background integration steps that are stored in the output vector
  */
-class_precision_parameter(background_Nloga,int,3000)
+class_precision_parameter(background_Nloga,int,40000)
 /**
  * Evolver to be used for thermodynamics (rk, ndf15)
  */
@@ -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 */

include/thermodynamics.h

@@ -519,6 +519,8 @@ extern "C" {
 
   int thermodynamics_free(struct thermodynamics * pth);
 
+  int thermodynamics_free_input(struct thermodynamics * pth);
+
   /* internal functions of the module */
 
   int thermodynamics_helium_from_bbn(struct precision * ppr,

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(