Few bug fixes and mnior improvements, especially on Halofit precision

* 2.6: changed version number to 2.6.2 add sigma8 as an input parameter using shooting method with target A_s Added the new quadrature functionality Added check for mvec library in glibc to avoid issue lesgourg#99 changed default value of k_max_for_pk to 1/Mpc for precise sigma8 even with default settings fixed new precision settings for halofit just debugged previous commit changes in spectra module: store pk_nl in table with size independent from pk_l, for more correct interpolation improved coding style inside nonlinear_halofit added tau_reio() in classy now possible to ask for tensors only improved estimate of diffusion damping scale using R corrections from W. Hu added Omega_g() to classy (to get Omega_gamma today) corrected new BBN file: uses Yp, not Yp^BBN; added marcucci et al. version as an alternative corrected new BBN file: uses Yp, not Yp^BBN; added marcucci et al. version as an alternative updated BBN file
abacusorg · Oct 23, 2017 · db1589b · db1589b
2 parents 5908e1c + 8e8b35d
commit db1589b
Show file tree

Hide file tree

Showing 22 changed files with 1,949 additions and 472 deletions.
diff --git a/bbn/sBBN_2017.dat b/bbn/sBBN_2017.dat
diff --git a/bbn/sBBN_2017_marcucci.dat b/bbn/sBBN_2017_marcucci.dat
diff --git a/cl_ref.pre b/cl_ref.pre
@@ -40,7 +40,8 @@ ur_fluid_trigger_tau_over_tau_k = 50.
 ncdm_fluid_approximation = 3
 ncdm_fluid_trigger_tau_over_tau_k = 51.
 
-tol_ncdm = 1.e-10
+tol_ncdm_synchronous = 1.e-10
+tol_ncdm_newtonian = 1.e-10
 
 l_logstep=1.026
 l_linstep=25
@@ -54,7 +55,7 @@ hyper_x_tol = 1.e-4
 hyper_flat_approximation_nu = 1.e6
 
 q_linstep=0.20
-q_logstep_spline= 20. 
+q_logstep_spline= 20.
 q_logstep_trapzd = 0.5
 q_numstep_transition = 250
 
@@ -73,6 +74,8 @@ transfer_neglect_delta_k_T_b = 100.
 neglect_CMB_sources_below_visibility = 1.e-30
 transfer_neglect_late_source = 3000.
 
+halofit_k_per_decade = 3000.
+
 l_switch_limber = 40.
 accurate_lensing=1
 num_mu_minus_lmax = 1000.

diff --git a/explanatory.ini b/explanatory.ini
@@ -138,6 +138,19 @@ ksi_ncdm =
 
 deg_ncdm =
 
+  -> 'Quadrature strategy' is the method used for the momentum sampling of the
+  ncdm distribution function. 0 is the automatic method, 1 is Gauss-Laguerre
+  quadrature, 2 is the trapezoidal rule on [0,Infinity] using the transformation
+  q->1/t-1. 3 is the trapezoidal rule on [0,q_max] where q_max is the next input.
+  (default: set to 0)
+
+Quadrature strategy =
+
+  -> 'Maximum q' is the maximum q relevant only for Quadrature strategy 3.
+  (default: set to 15)
+
+Maximum q =
+
 7) curvature: 'Omega_k' (default: 'Omega_k' set to 0)
 
 Omega_k = 0.
@@ -675,7 +688,7 @@ l_max_tensors = 500
 2) maximum k in P(k), 'P_k_max_h/Mpc' in units of h/Mpc or 'P_k_max_1/Mpc' in
    units of 1/Mpc. If scalar Cls are also requested, a minimum value is
    automatically imposed (the same as in scalar Cls computation)
-   (default: set to 0.1h/Mpc)
+   (default: set to 1 1/Mpc)
 
 P_k_max_h/Mpc = 1.
 #P_k_max_1/Mpc = 0.7

diff --git a/include/background.h b/include/background.h
@@ -280,7 +280,9 @@ struct background
 
 
   //@{
-
+  int * ncdm_quadrature_strategy; /**< Vector of integers according to quadrature strategy. */
+  int * ncdm_input_q_size; /**< Vector of numbers of q bins */
+  double * ncdm_qmax;   /**< Vector of maximum value of q */
   double ** q_ncdm_bg;  /**< Pointers to vectors of background sampling in q */
   double ** w_ncdm_bg;  /**< Pointers to vectors of corresponding quadrature weights w */
   double ** q_ncdm;     /**< Pointers to vectors of perturbation sampling in q */

diff --git a/include/common.h b/include/common.h
@@ -15,7 +15,7 @@
 #ifndef __COMMON__
 #define __COMMON__
 
-#define _VERSION_ "v2.6.1"
+#define _VERSION_ "v2.6.2"
 /* @cond INCLUDE_WITH_DOXYGEN */
 
 #define _TRUE_ 1 /**< integer associated to true statement */
@@ -744,33 +744,36 @@ struct precision
 
   /** parameters relevant for HALOFIT computation */
 
-  double halofit_dz; /**< spacing in redshift space defining values of z
-			at which HALOFIT will be used. Intermediate
-			values will be obtained by
-			interpolation. Decrease for more precise
-			interpolations, at the expense of increasing
-			time spent in nonlinear_init() */
+  double halofit_min_k_nonlinear; /**< value of k in 1/Mpc below which
+				     non-linear corrections will be neglected */
 
-  double halofit_min_k_nonlinear; /**< value of k in 1/Mpc above
-				     which non-linear corrections will
-				     be computed */
-
-  double halofit_sigma_precision; /**< a smaller value will lead to a
-				      more precise halofit result at
-				      the highest requested redshift,
-				      at the expense of requiring a
-				      larger k_max */
-
-  double halofit_min_k_max; /**< when halofit is used, k_max must be at
-                               least equal to this value (otherwise
+  double halofit_min_k_max; /**< when halofit is used, k_max must be
+                               at least equal to this value (otherwise
                                halofit could not find the scale of
-                               non-linearity) */
+                               non-linearity). Calculations are done
+                               internally until this k_max, but the
+                               output is still controlled by
+                               P_k_max_1/Mpc or P_k_max_h/Mpc even if
+                               they are smaller */
 
   double halofit_k_per_decade; /**< halofit needs to evalute integrals
                                   (linear power spectrum times some
                                   kernels). They are sampled using
                                   this logarithmic step size. */
 
+  double halofit_sigma_precision; /**< a smaller value will lead to a
+				      more precise halofit result at the *highest*
+				      redshift at which halofit can make computations,
+				      at the expense of requiring a larger k_max; but
+				      this parameter is not relevant for the
+				      precision on P_nl(k,z) at other redshifts, so
+				      there is normally no need to change it */
+
+  double halofit_tol_sigma; /**< tolerance required on sigma(R) when
+                               matching the condition sigma(R_nl)=1,
+                               whcih defines the wavenumber of
+                               non-linearity, k_nl=1./R_nl */
+
   //@}
 
   /** @name - parameters related to lensing */

diff --git a/include/input.h b/include/input.h
@@ -134,10 +134,10 @@
  * temporary parameters for background fzero function
  */
 
-enum target_names {theta_s, Omega_dcdmdr, omega_dcdmdr, Omega_scf, Omega_ini_dcdm, omega_ini_dcdm};
+enum target_names {theta_s, Omega_dcdmdr, omega_dcdmdr, Omega_scf, Omega_ini_dcdm, omega_ini_dcdm, sigma8};
 enum computation_stage {cs_background, cs_thermodynamics, cs_perturbations,
                         cs_primordial, cs_nonlinear, cs_transfer, cs_spectra};
-#define _NUM_TARGETS_ 6 //Keep this number as number of target_names
+#define _NUM_TARGETS_ 7 //Keep this number as number of target_names
 
 struct input_pprpba {
   struct precision * ppr;
@@ -155,7 +155,7 @@ struct fzerofun_workspace {
 
 
 /**************************************************************/
-/* @cond INCLUDE_WITH_DOXYGEN */ 
+/* @cond INCLUDE_WITH_DOXYGEN */
 /*
  * Boilerplate for C++
  */
@@ -284,4 +284,4 @@ extern "C" {
 /**************************************************************/
 
 #endif
-/* @endcond */
+/* @endcond */