Commit 1b0882f6 authored by Julien Lesgourgues's avatar Julien Lesgourgues
Browse files

Corrected a few bugs affecting Halofit precision with massive nu or decaying DM

* 2.4:
  version number changed to v2.4.5
  included dcdm in definition of fnu in halofit
  series of fixes to remove numerical errors from halofit implementation
  relevant only for quintessence: corrected an abs() into fabs()
  renoved selection bias from shear spectrum, only relevant for number count (thanks Aurelien Benoit-Levy)
  corrected bug in interface to halofit, problematic only with decaying CDM (Omega_drdcdm>0)
  improved reco-reio matching, helps to get more precise z_reio(tau_reio) search
  Improved error handling for shooting method not to generate input errors for MontePython
parents 54f22882 049b2f85
......@@ -262,7 +262,18 @@ int array_integrate_all_trapzd_or_spline(
int * last_index,
double * result,
int result_size,
ErrorMsg errmsg); /** from 1 to n_columns */
ErrorMsg errmsg);
int array_interpolate_one_growing_closeby(
double * array,
int n_columns,
int n_lines,
int index_x, /** from 0 to (n_columns-1) */
double x,
int * last_index,
int index_y,
double * result,
ErrorMsg errmsg);
int array_interpolate_spline_growing_closeby(
double * x_array,
......
......@@ -303,6 +303,10 @@ struct background
//@{
short shooting_failed; /**< flag is set to true if shooting failed. */
ErrorMsg shooting_error; /**< Error message from shooting failed. */
short background_verbose; /**< flag regulating the amount of information sent to standard output (none if set to zero) */
ErrorMsg error_message; /**< zone for writing error messages */
......
......@@ -15,7 +15,7 @@
#ifndef __COMMON__
#define __COMMON__
#define _VERSION_ "v2.4.4"
#define _VERSION_ "v2.4.5"
#define _TRUE_ 1 /**< integer associated to true statement */
#define _FALSE_ 0 /**< integer associated to false statement */
......@@ -100,6 +100,14 @@ int get_number_of_titles(char * titlestring);
} \
}
/* macro for trying to call function */
#define class_call_try(function, error_message_from_function, error_message_output,list_of_commands) { \
if (function == _FAILURE_) { \
class_call_message(error_message_output,#function,error_message_from_function); \
list_of_commands; \
} \
}
/* macro for calling function and returning error if it failed */
#define class_call(function, error_message_from_function, error_message_output) \
class_call_except(function, error_message_from_function,error_message_output,)
......@@ -737,6 +745,11 @@ struct precision
halofit could not find the scale of
non-linearity) */
double halofit_k_per_decade; /* halofit needs to evalute integrals
(linear power spectrum times some
kernels). They are sampled using
this logarithmic step size. */
//@}
/** @name - parameters related to lensing */
......
......@@ -259,6 +259,11 @@ extern "C" {
struct fzerofun_workspace * pfzw,
ErrorMsg errmsg);
int input_find_root(double *xzero,
int *fevals,
struct fzerofun_workspace *pfzw,
ErrorMsg errmsg);
int file_exists(const char *fname);
int input_auxillary_target_conditions(struct file_content * pfc,
......
......@@ -88,7 +88,10 @@ extern "C" {
struct primordial *ppm,
struct nonlinear *pnl,
int index_tau,
double *pk_l);
double *pk_l,
double *lnk,
double *lnpk,
double *ddlnpk);
int nonlinear_halofit(
struct precision *ppr,
......@@ -97,6 +100,9 @@ extern "C" {
struct nonlinear *pnl,
double tau,
double *pk_l,
double *lnk,
double *lnpk,
double *ddlnpk,
double *pk_nl,
double *k_nl
);
......
......@@ -520,6 +520,12 @@ int background_init(
}
}
/** - if shooting failed during input, catch the error here: */
class_test(pba->shooting_failed == _TRUE_,
pba->error_message,
"Shooting failed, try optimising input_get_guess(). Error message:\n\n%s",
pba->shooting_error);
/** - assign values to all indices in vectors of background quantities with background_indices()*/
class_call(background_indices(pba),
pba->error_message,
......
......@@ -233,7 +233,7 @@ int input_init(
enum computation_stage target_cs[] = {cs_thermodynamics, cs_background, cs_background,
cs_background, cs_background, cs_background};
int input_verbose = 0, int1, aux_flag;
int input_verbose = 0, int1, aux_flag, shooting_failed=_FALSE_;
class_read_int("input_verbose",input_verbose);
......@@ -318,68 +318,14 @@ int input_init(
if (unknown_parameters_size == 1){
/* We can do 1 dimensional root finding */
/** Here is our guess: */
class_call(input_get_guess(&x1, &dxdy, &fzw, errmsg),
errmsg, errmsg);
// printf("x1= %g\n",x1);
class_call(input_fzerofun_1d(x1,
&fzw,
&f1,
errmsg),
errmsg, errmsg);
fevals++;
//printf("x1= %g, f1= %g\n",x1,f1);
dx = 1.5*f1*dxdy;
/** Do linear hunt for boundaries: */
for (iter=1; iter<=15; iter++){
//x2 = x1 + search_dir*dx;
x2 = x1 - dx;
for (iter2=1; iter2 <= 3; iter2++) {
return_function = input_fzerofun_1d(x2,&fzw,&f2,errmsg);
fevals++;
//printf("x2= %g, f2= %g\n",x2,f2);
//fprintf(stderr,"iter2=%d\n",iter2);
if (return_function ==_SUCCESS_) {
break;
}
else if (iter2 < 3) {
dx*=0.5;
x2 = x1-dx;
}
else {
//fprintf(stderr,"get here\n");
class_stop(errmsg,errmsg);
}
}
if (f1*f2<0.0){
/** root has been bracketed */
if (0==1){//(pba->background_verbose > 4){
printf("Root has been bracketed after %d iterations: [%g, %g].\n",iter,x1,x2);
}
break;
}
x1 = x2;
f1 = f2;
}
/** Find root using Ridders method. (Exchange for bisection if you are old-school.) */
class_call(class_fzero_ridder(input_fzerofun_1d,
x1,
x2,
1e-5*MAX(fabs(x1),fabs(x2)),
&fzw,
&f1,
&f2,
&xzero,
&fevals,
errmsg),
errmsg, errmsg);
/* If shooting fails, postpone error to background module to play nice with MontePython. */
class_call_try(input_find_root(&xzero,
&fevals,
&fzw,
errmsg),
errmsg,
pba->shooting_error,
shooting_failed=_TRUE_);
/* Store xzero */
sprintf(fzw.fc.value[fzw.unknown_parameters_index[0]],"%e",xzero);
......@@ -403,16 +349,16 @@ int input_init(
errmsg),
errmsg, errmsg);
class_call(fzero_Newton(input_try_unknown_parameters,
x_inout,
dxdF,
unknown_parameters_size,
1e-4,
1e-6,
&fzw,
&fevals,
errmsg),
errmsg,errmsg);
class_call_try(fzero_Newton(input_try_unknown_parameters,
x_inout,
dxdF,
unknown_parameters_size,
1e-4,
1e-6,
&fzw,
&fevals,
errmsg),
errmsg, pba->shooting_error,shooting_failed=_TRUE_);
if (input_verbose > 0) {
fprintf(stdout,"Computing unknown input parameters\n");
......@@ -454,6 +400,9 @@ int input_init(
errmsg,
errmsg);
/** Set status of shooting: */
pba->shooting_failed = shooting_failed;
/* all parameters read in fzw must be considered as read in
pfc. At the same time the parameters read before in pfc (like
theta_s,...) must still be considered as read (hence we could
......@@ -1046,7 +995,7 @@ int input_read_parameters(
class_read_double("scf_shooting_parameter",pba->scf_parameters[pba->scf_tuning_index]);
scf_lambda = pba->scf_parameters[0];
if ((abs(scf_lambda) <3.)&&(pba->background_verbose>1))
if ((fabs(scf_lambda) <3.)&&(pba->background_verbose>1))
printf("lambda = %e <3 won't be tracking (for exp quint) unless overwritten by tuning function\n",scf_lambda);
class_call(parser_read_string(pfc,
......@@ -2488,6 +2437,7 @@ int input_read_parameters(
class_read_double("halofit_dz",ppr->halofit_dz);
class_read_double("halofit_min_k_nonlinear",ppr->halofit_min_k_nonlinear);
class_read_double("halofit_k_per_decade",ppr->halofit_k_per_decade);
class_read_double("halofit_sigma_precision",ppr->halofit_sigma_precision);
/** h.8. parameter related to lensing */
......@@ -2658,6 +2608,8 @@ int input_default_params(
pba->wa_fld=0.;
pba->cs2_fld=1.;
pba->shooting_failed = _FALSE_;
/** - thermodynamics structure */
pth->YHe=_BBN_;
......@@ -3117,6 +3069,7 @@ int input_default_precision ( struct precision * ppr ) {
ppr->halofit_dz=0.1;
ppr->halofit_min_k_nonlinear=0.0035;
ppr->halofit_k_per_decade = 80.;
ppr->halofit_sigma_precision=0.05;
ppr->halofit_min_k_max=5.;
......@@ -3593,6 +3546,79 @@ int input_get_guess(double *xguess,
return _SUCCESS_;
}
int input_find_root(double *xzero,
int *fevals,
struct fzerofun_workspace *pfzw,
ErrorMsg errmsg){
double x1, x2, f1, f2, dxdy, dx;
int iter, iter2;
int return_function;
/** Here is our guess: */
class_call(input_get_guess(&x1, &dxdy, pfzw, errmsg),
errmsg, errmsg);
// printf("x1= %g\n",x1);
class_call(input_fzerofun_1d(x1,
pfzw,
&f1,
errmsg),
errmsg, errmsg);
(*fevals)++;
//printf("x1= %g, f1= %g\n",x1,f1);
dx = 1.5*f1*dxdy;
/** Do linear hunt for boundaries: */
for (iter=1; iter<=15; iter++){
//x2 = x1 + search_dir*dx;
x2 = x1 - dx;
for (iter2=1; iter2 <= 3; iter2++) {
return_function = input_fzerofun_1d(x2,pfzw,&f2,errmsg);
(*fevals)++;
//printf("x2= %g, f2= %g\n",x2,f2);
//fprintf(stderr,"iter2=%d\n",iter2);
if (return_function ==_SUCCESS_) {
break;
}
else if (iter2 < 3) {
dx*=0.5;
x2 = x1-dx;
}
else {
//fprintf(stderr,"get here\n");
class_stop(errmsg,errmsg);
}
}
if (f1*f2<0.0){
/** root has been bracketed */
if (0==1){
printf("Root has been bracketed after %d iterations: [%g, %g].\n",iter,x1,x2);
}
break;
}
x1 = x2;
f1 = f2;
}
/** Find root using Ridders method. (Exchange for bisection if you are old-school.)*/
class_call(class_fzero_ridder(input_fzerofun_1d,
x1,
x2,
1e-5*MAX(fabs(x1),fabs(x2)),
pfzw,
&f1,
&f2,
xzero,
fevals,
errmsg),
errmsg,errmsg);
return _SUCCESS_;
}
int file_exists(const char *fname){
FILE *file = fopen(fname, "r");
if (file != NULL){
......
......@@ -63,6 +63,9 @@ int nonlinear_init(
int index_tau;
double *pk_l;
double *pk_nl;
double *lnk_l;
double *lnpk_l;
double *ddlnpk_l;
/** (a) if non non-linear corrections requested */
......@@ -102,12 +105,16 @@ int nonlinear_init(
class_alloc(pk_l,pnl->k_size*sizeof(double),pnl->error_message);
class_alloc(pk_nl,pnl->k_size*sizeof(double),pnl->error_message);
class_alloc(lnk_l,pnl->k_size*sizeof(double),pnl->error_message);
class_alloc(lnpk_l,pnl->k_size*sizeof(double),pnl->error_message);
class_alloc(ddlnpk_l,pnl->k_size*sizeof(double),pnl->error_message);
/** - loop over time */
for (index_tau = pnl->tau_size-1; index_tau>=0; index_tau--) {
/* get P_L(k) at this time */
class_call(nonlinear_pk_l(ppt,ppm,pnl,index_tau,pk_l),
class_call(nonlinear_pk_l(ppt,ppm,pnl,index_tau,pk_l,lnk_l,lnpk_l,ddlnpk_l),
pnl->error_message,
pnl->error_message);
......@@ -127,6 +134,9 @@ int nonlinear_init(
pnl->tau[index_tau],
pk_l,
pk_nl,
lnk_l,
lnpk_l,
ddlnpk_l,
&(pnl->k_nl[index_tau])) == _SUCCESS_) {
/*
......@@ -158,6 +168,10 @@ int nonlinear_init(
free(pk_l);
free(pk_nl);
free(lnk_l);
free(lnpk_l);
free(ddlnpk_l);
}
else {
......@@ -192,8 +206,10 @@ int nonlinear_pk_l(
struct primordial *ppm,
struct nonlinear *pnl,
int index_tau,
double *pk_l
) {
double *pk_l,
double *lnk,
double *lnpk,
double *ddlnpk) {
int index_md;
int index_k;
......@@ -255,8 +271,20 @@ int nonlinear_pk_l(
}
}
lnk[index_k] = log(pnl->k[index_k]);
lnpk[index_k] = log(pk_l[index_k]);
}
class_call(array_spline_table_columns(lnk,
pnl->k_size,
lnpk,
1,
ddlnpk,
_SPLINE_NATURAL_,
pnl->error_message),
pnl->error_message,
pnl->error_message);
free(primordial_pk);
return _SUCCESS_;
......@@ -271,6 +299,9 @@ int nonlinear_halofit(
double tau,
double *pk_l,
double *pk_nl,
double *lnk_l,
double *lnpk_l,
double *ddlnpk_l,
double *k_nl
) {
......@@ -296,20 +327,83 @@ int nonlinear_halofit(
double anorm;
double *integrand_array;
int integrand_size;
int index_ia_k;
int index_ia_pk;
int index_ia_sum1;
int index_ia_ddsum1;
int index_ia_sum2;
int index_ia_ddsum2;
int index_ia_sum3;
int index_ia_ddsum3;
int ia_size;
int index_ia;
double k_integrand;
double lnpk_integrand;
double x2,R;
class_alloc(pvecback,pba->bg_size*sizeof(double),pnl->error_message);
Omega0_m = (pba->Omega0_cdm + pba->Omega0_b + pba->Omega0_ncdm_tot);
Omega0_m = (pba->Omega0_cdm + pba->Omega0_b + pba->Omega0_ncdm_tot + pba->Omega0_dcdm);
w0 = pba->w0_fld;
fnu = pba->Omega0_ncdm_tot/(pba->Omega0_b + pba->Omega0_cdm+pba->Omega0_ncdm_tot);
fnu = pba->Omega0_ncdm_tot/Omega0_m;
anorm = 1./(2*pow(_PI_,2));
class_alloc(integrand_array,pnl->k_size*7*sizeof(double),pnl->error_message);
/* Until the 17.02.2015 the values of k used for integrating sigma(R) quantities needed by Halofit where the same as in the perturbation module.
Since then, we sample these integrals on more values, in order to get more precise integrals (thanks Matteo Zennaro for noticing the need for this).
We create a temporary integrand_array which columns will be:
- k in 1/Mpc
- just linear P(k) in Mpc**3
- 1/(2(pi**2)) P(k) k**2 exp(-(kR)**2)
- second derivative of previous line with spline
- 1/(2(pi**2)) P(k) k**2 2 (kR) exp(-(kR)**2)
- second derivative of previous line with spline
- 1/(2(pi**2)) P(k) k**2 4 (kR)(1-kR) exp(-(kR)**2)
- second derivative of previous line with spline
*/
index_ia=0;
class_define_index(index_ia_k, _TRUE_,index_ia,1);
class_define_index(index_ia_pk, _TRUE_,index_ia,1);
class_define_index(index_ia_sum1, _TRUE_,index_ia,1);
class_define_index(index_ia_ddsum1,_TRUE_,index_ia,1);
class_define_index(index_ia_sum2, _TRUE_,index_ia,1);
class_define_index(index_ia_ddsum2,_TRUE_,index_ia,1);
class_define_index(index_ia_sum3, _TRUE_,index_ia,1);
class_define_index(index_ia_ddsum3,_TRUE_,index_ia,1);
ia_size = index_ia;
integrand_size=(int)(log(pnl->k[pnl->k_size-1]/pnl->k[0])/log(10.)*ppr->halofit_k_per_decade)+1;
class_alloc(integrand_array,integrand_size*ia_size*sizeof(double),pnl->error_message);
/* we fill integrand_array with values of k and P(k) using interpolation */
last_index=0;
for (index_k=0; index_k < integrand_size; index_k++) {
k_integrand=pnl->k[0]*pow(10.,index_k/ppr->halofit_k_per_decade);
class_call(array_interpolate_spline(lnk_l,
pnl->k_size,
lnpk_l,
ddlnpk_l,
1,
log(k_integrand),
&last_index,
&lnpk_integrand,
1,
pnl->error_message),
pnl->error_message,
pnl->error_message);
integrand_array[index_k*ia_size + index_ia_k] = k_integrand;
integrand_array[index_k*ia_size + index_ia_pk] = exp(lnpk_integrand);
for (index_k=0; index_k < pnl->k_size; index_k++) {
integrand_array[index_k*7 + 0] = pnl->k[index_k];
}
class_call(background_at_tau(pba,tau,pba->long_info,pba->inter_normal,&last_index,pvecback),
......@@ -320,20 +414,35 @@ int nonlinear_halofit(
Omega_v = 1.-pvecback[pba->index_bg_Omega_m]-pvecback[pba->index_bg_Omega_r];
/* minimum value of R such that the integral giving sigma_R is converged */
R=sqrt(-log(ppr->halofit_sigma_precision))/pnl->k[pnl->k_size-1];
R=sqrt(-log(ppr->halofit_sigma_precision))/integrand_array[(integrand_size-1)*ia_size + index_ia_k];
/* corresponding value of sigma_R */
sum1=0.;
for (index_k=0; index_k < pnl->k_size; index_k++) {
x2 = pow(pnl->k[index_k]*R,2);
integrand_array[index_k*7 + 1] = pk_l[index_k]*pow(pnl->k[index_k],2)*anorm*exp(-x2);
for (index_k=0; index_k < integrand_size; index_k++) {
x2 = pow(integrand_array[index_k*ia_size + index_ia_k]*R,2);
integrand_array[index_k*ia_size + index_ia_sum1] = integrand_array[index_k*ia_size + index_ia_pk]
*pow(integrand_array[index_k*ia_size + index_ia_k],2)*anorm*exp(-x2);
}
/* fill in second derivatives */
class_call(array_spline(integrand_array,7,pnl->k_size,0,1,2,_SPLINE_EST_DERIV_,pnl->error_message),
class_call(array_spline(integrand_array,
ia_size,
integrand_size,
index_ia_k,
index_ia_sum1,
index_ia_ddsum1,
_SPLINE_EST_DERIV_,
pnl->error_message),
pnl->error_message,
pnl->error_message);
/* integrate */
class_call(array_integrate_all_spline(integrand_array,7,pnl->k_size,0,1,2,&sum1,pnl->error_message),
class_call(array_integrate_all_spline(integrand_array,
ia_size,
integrand_size,
index_ia_k,
index_ia_sum1,
index_ia_ddsum1,
&sum1,
pnl->error_message),
pnl->error_message,
pnl->error_message);
sigma = sqrt(sum1);
......@@ -352,16 +461,31 @@ int nonlinear_halofit(
/* corresponding value of sigma_R */
sum1=0.;
for (index_k=0; index_k < pnl->k_size; index_k++) {
x2 = pow(pnl->k[index_k]*R,2);
integrand_array[index_k*7 + 1] = pk_l[index_k]*pow(pnl->k[index_k],2)*anorm*exp(-x2);
for (index_k=0; index_k < integrand_size; index_k++) {
x2 = pow(integrand_array[index_k*ia_size + index_ia_k]*R,2);
integrand_array[index_k*ia_size + index_ia_sum1] = integrand_array[index_k*ia_size + index_ia_pk]
*pow(integrand_array[index_k*ia_size + index_ia_k],2)*anorm*exp(-x2);
}
/* fill in second derivatives */
class_call(array_spline(integrand_array,7,pnl->k_size,0,1,2,_SPLINE_EST_DERIV_,pnl->error_message),
class_call(array_spline(integrand_array,
ia_size,
integrand_size,
index_ia_k,
index_ia_sum1,
index_ia_ddsum1,
_SPLINE_EST_DERIV_,
pnl->error_message),
pnl->error_message,
pnl->error_message);
/* integrate */
class_call(array_integrate_all_spline(integrand_array,7,pnl->k_size,0,1,2,&sum1,pnl->error_message),
class_call(array_integrate_all_spline(integrand_array,
ia_size,
integrand_size,
index_ia_k,
index_ia_sum1,
index_ia_ddsum1,
&sum1,
pnl->error_message),
pnl->error_message,
pnl->error_message);
sigma = sqrt(sum1);
......@@ -384,31 +508,81 @@ int nonlinear_halofit(
sum2=0.;
sum3=0.;
for (index_k=0; index_k < pnl->k_size; index_k++) {
x2 = pnl->k[index_k]*pnl->k[index_k]*rmid*rmid;
integrand_array[index_k*7 + 1] = pk_l[index_k]*pow(pnl->k[index_k],2)*anorm*exp(-x2);
integrand_array[index_k*7 + 3] = pk_l[index_k]*pow(pnl->k[index_k],2)*anorm*2.*x2*exp(-x2);
integrand_array[index_k*7 + 5] = pk_l[index_k]*pow(pnl->k[index_k],2)*anorm*4.*x2*(1.-x2)*exp(-x2);
for (index_k=0; index_k < integrand_size; index_k++) {
x2 = pow(integrand_array[index_k*ia_size + index_ia_k],2)*rmid*rmid;
integrand_array[index_k*ia_size + index_ia_sum1] = integrand_array[index_k*ia_size + index_ia_pk]
*pow(integrand_array[index_k*ia_size + index_ia_k],2)*anorm*exp(-x2);
integrand_array[index_k*ia_size + index_ia_sum2] = integrand_array[index_k*