Commit a3907e03 authored by lesgourg's avatar lesgourg
Browse files

final polishing of pk_eq method

parent 61bfe207
......@@ -780,7 +780,10 @@ struct precision
whcih defines the wavenumber of
non-linearity, k_nl=1./R_nl */
double pk_eq_z_max; /**< Maximum z until the pk_eq method tries to calculate */
double pk_eq_z_max; /**< Maximum z until which the Pk_equal method of 0810.0190 and 1601.07230 is used */
double pk_eq_tol; /**< tolerance for finding the equivalent models of the pk_equal method */
//@}
/** @name - parameters related to lensing */
......
......@@ -54,15 +54,15 @@ struct nonlinear {
short has_pk_eq; /**< flag: will we use the pk_eq method? */
int index_eq_w; /**< index of w in table eq_w_and_Omega */
int index_eq_Omega_m; /**< index of Omega_m in table eq_w_and_Omega */
int eq_size; /**< number of indices in table eq_w_and_Omega */
int index_pk_eq_w; /**< index of w in table pk_eq_w_and_Omega */
int index_pk_eq_Omega_m; /**< index of Omega_m in table pk_eq_w_and_Omega */
int pk_eq_size; /**< number of indices in table pk_eq_w_and_Omega */
int eq_tau_size; /**< number of times (and raws in table eq_w_and_Omega) */
int pk_eq_tau_size; /**< number of times (and raws in table pk_eq_w_and_Omega) */
double * eq_tau; /**< table of time values */
double * eq_w_and_Omega; /**< table of background quantites */
double * eq_ddw_and_ddOmega; /**< table of second derivatives */
double * pk_eq_tau; /**< table of time values */
double * pk_eq_w_and_Omega; /**< table of background quantites */
double * pk_eq_ddw_and_ddOmega; /**< table of second derivatives */
//@{
......
......@@ -2783,6 +2783,7 @@ int input_read_parameters(
class_read_double("halofit_sigma_precision",ppr->halofit_sigma_precision);
class_read_double("halofit_tol_sigma",ppr->halofit_tol_sigma);
class_read_double("pk_eq_z_max",ppr->pk_eq_z_max);
class_read_double("pk_eq_tol",ppr->pk_eq_tol);
/** - (h.7.) parameter related to lensing */
......@@ -2862,8 +2863,8 @@ int input_read_parameters(
}
/** - (i.5) special buisness if we want Halofit with wa_fld non-zero:
so-called "Pk-equal method" of 0810.0190 and 1601.07230 */
/** - (i.5) special steps if we want Halofit with wa_fld non-zero:
so-called "Pk_equal method" of 0810.0190 and 1601.07230 */
if ((pnl->method == nl_halofit) && (pba->Omega0_fld != 0.) && (pba->wa_fld != 0.))
pnl->has_pk_eq = _TRUE_;
......@@ -2872,7 +2873,7 @@ int input_read_parameters(
if (input_verbose > 0) {
printf(" -> since you want to use Halofit with a non-zero wa_fld, calling background module to\n");
printf(" extract the effective w(tau), Omega_m(tau) parameters required by the Pk-equal method\n");
printf(" extract the effective w(tau), Omega_m(tau) parameters required by the Pk_equal method\n");
}
class_call(input_prepare_pk_eq(ppr,pba,pth,pnl,input_verbose,errmsg),
errmsg,
......@@ -3462,7 +3463,8 @@ int input_default_precision ( struct precision * ppr ) {
ppr->halofit_sigma_precision = 0.05;
ppr->halofit_tol_sigma = 1.e-6;
ppr->pk_eq_z_max = 5.;
ppr->pk_eq_tol = 1.e-7;
/**
* - parameter related to lensing
*/
......@@ -4111,6 +4113,26 @@ int compare_doubles(const void *a,const void *b) {
return 0;
}
/**
* Perform preliminary steps fur using the method called Pk_equal,
* described in 0810.0190 and 1601.07230, extending the range of
* validity of HALOFIT from constant w to (w0,wa) models. In that
* case, one must compute here some effective values of w0_eff(z_i)
* and Omega_m_eff(z_i), that will be interpolated later at arbitrary
* redshift in the non-linear module.
*
* Returns table of values [z_i, tau_i, w0_eff_i, Omega_m_eff_i]
* stored in nonlinear structure.
*
* @param ppr Input: pointer to precision structure
* @param pba Input: pointer to background structure
* @param pth Input: pointer to thermodynamics structure
* @param pnl Input/Output: pointer to nonlinear structure
* @param input_verbose Input: verbosity of this input module
* @param errmsg Input/Ouput: error message
*/
int input_prepare_pk_eq(
struct precision * ppr,
struct background *pba,
......@@ -4120,13 +4142,17 @@ int input_prepare_pk_eq(
ErrorMsg errmsg
) {
/** Summary: */
/** - define local variables */
double tau_of_z;
double delta_tau;
double error;
double delta_tau_eq;
double * pvecback;
int last_index=0;
int index_eq_z;
int index_pk_eq_z;
int index_eq;
int true_background_verbose;
int true_thermodynamics_verbose;
......@@ -4134,39 +4160,59 @@ int input_prepare_pk_eq(
double true_wa_fld;
double * z;
/** - store the true cosmological parameters (w0, wa) somwhere before using temporarily some fake ones in this function */
true_background_verbose = pba->background_verbose;
true_thermodynamics_verbose = pth->thermodynamics_verbose;
true_w0_fld = pba->w0_fld;
true_wa_fld = pba->wa_fld;
////
/** - the fake calls of the background and thermodynamics module will be done in non-verbose mode */
pba->background_verbose = 0;
pth->thermodynamics_verbose = 0;
pnl->eq_tau_size = 10;
class_alloc(pnl->eq_tau,pnl->eq_tau_size*sizeof(double),errmsg);
class_alloc(z,pnl->eq_tau_size*sizeof(double),errmsg);
/** - allocate indices and arrays for storing the results */
pnl->pk_eq_tau_size = 10;
class_alloc(pnl->pk_eq_tau,pnl->pk_eq_tau_size*sizeof(double),errmsg);
class_alloc(z,pnl->pk_eq_tau_size*sizeof(double),errmsg);
index_eq = 0;
class_define_index(pnl->index_eq_w,_TRUE_,index_eq,1);
class_define_index(pnl->index_eq_Omega_m,_TRUE_,index_eq,1);
pnl->eq_size = index_eq;
class_alloc(pnl->eq_w_and_Omega,pnl->eq_tau_size*pnl->eq_size*sizeof(double),errmsg);
class_alloc(pnl->eq_ddw_and_ddOmega,pnl->eq_tau_size*pnl->eq_size*sizeof(double),errmsg);
class_define_index(pnl->index_pk_eq_w,_TRUE_,index_eq,1);
class_define_index(pnl->index_pk_eq_Omega_m,_TRUE_,index_eq,1);
pnl->pk_eq_size = index_eq;
class_alloc(pnl->pk_eq_w_and_Omega,pnl->pk_eq_tau_size*pnl->pk_eq_size*sizeof(double),errmsg);
class_alloc(pnl->pk_eq_ddw_and_ddOmega,pnl->pk_eq_tau_size*pnl->pk_eq_size*sizeof(double),errmsg);
/** - call the background module in order to fill a table of tau_i[z_i] */
class_call(background_init(ppr,pba), pba->error_message, errmsg);
for (index_eq_z=0; index_eq_z<pnl->eq_tau_size; index_eq_z++) {
z[index_eq_z] = exp(log(1.+ppr->pk_eq_z_max)/(pnl->eq_tau_size-1)*index_eq_z)-1.;
class_call(background_tau_of_z(pba,z[index_eq_z],&tau_of_z),
for (index_pk_eq_z=0; index_pk_eq_z<pnl->pk_eq_tau_size; index_pk_eq_z++) {
z[index_pk_eq_z] = exp(log(1.+ppr->pk_eq_z_max)/(pnl->pk_eq_tau_size-1)*index_pk_eq_z)-1.;
class_call(background_tau_of_z(pba,z[index_pk_eq_z],&tau_of_z),
pba->error_message, errmsg);
pnl->eq_tau[index_eq_z] = tau_of_z;
pnl->pk_eq_tau[index_pk_eq_z] = tau_of_z;
}
class_call(background_free_noinput(pba), pba->error_message, errmsg);
for (index_eq_z=0; index_eq_z<pnl->eq_tau_size; index_eq_z++) {
/** - loop over z_i values. For each of them, we will call the
background and thermodynamics module for fake models. The goal is
to find, for each z_i, and effective w0_eff[z_i] and
Omega_m_eff{z_i], such that: the true model with (w0,wa) and the
equivalent model with (w0_eff[z_i],0) have the same conformal
distance between z_i and z_recombination, namely chi = tau[z_i] -
tau_rec. It is thus necessary to call both the background and
thermodynamics module for each fake model and to re-compute
tau_rec for each of them. Once the eqauivalent model is found we
compute and store Omega_m_effa(z_i) of the equivalent model */
for (index_pk_eq_z=0; index_pk_eq_z<pnl->pk_eq_tau_size; index_pk_eq_z++) {
if (input_verbose > 2)
printf(" * computing PK-equal parameters at z=%e\n",z[index_eq_z]);
printf(" * computing Pk_equal parameters at z=%e\n",z[index_pk_eq_z]);
/* get chi = (tau[z_i] - tau_rec) in true model */
pba->w0_fld = true_w0_fld;
pba->wa_fld = true_wa_fld;
......@@ -4174,9 +4220,9 @@ int input_prepare_pk_eq(
class_call(background_init(ppr,pba), pba->error_message, errmsg);
class_call(thermodynamics_init(ppr,pba,pth), pth->error_message, errmsg);
delta_tau = pnl->eq_tau[index_eq_z] - pth->tau_rec;
delta_tau = pnl->pk_eq_tau[index_pk_eq_z] - pth->tau_rec;
///////
/* launch iterations in order to coverge to effective model with wa=0 but the same chi = (tau[z_i] - tau_rec) */
pba->wa_fld=0.;
......@@ -4185,7 +4231,7 @@ int input_prepare_pk_eq(
class_call(thermodynamics_free(pth), pth->error_message, errmsg);
class_call(background_init(ppr,pba), pba->error_message, errmsg);
class_call(background_tau_of_z(pba,z[index_eq_z],&tau_of_z), pba->error_message, errmsg);
class_call(background_tau_of_z(pba,z[index_pk_eq_z],&tau_of_z), pba->error_message, errmsg);
class_call(thermodynamics_init(ppr,pba,pth), pth->error_message, errmsg);
delta_tau_eq = tau_of_z - pth->tau_rec;
......@@ -4194,9 +4240,11 @@ int input_prepare_pk_eq(
pba->w0_fld = pba->w0_fld*pow(1.+error,10.);
}
while(fabs(error) > 1.e-7);
while(fabs(error) > ppr->pk_eq_tol);
/* Equivalent model found. Store w0(z) in that model. Find Omega_m(z) in that model and store it. */
pnl->eq_w_and_Omega[pnl->eq_size*index_eq_z+pnl->index_eq_w] = pba->w0_fld;
pnl->pk_eq_w_and_Omega[pnl->pk_eq_size*index_pk_eq_z+pnl->index_pk_eq_w] = pba->w0_fld;
class_alloc(pvecback,pba->bg_size*sizeof(double),pba->error_message);
class_call(background_at_tau(pba,
......@@ -4206,7 +4254,7 @@ int input_prepare_pk_eq(
&last_index,
pvecback),
pba->error_message, errmsg);
pnl->eq_w_and_Omega[pnl->eq_size*index_eq_z+pnl->index_eq_Omega_m] = pvecback[pba->index_bg_Omega_m];
pnl->pk_eq_w_and_Omega[pnl->pk_eq_size*index_pk_eq_z+pnl->index_pk_eq_Omega_m] = pvecback[pba->index_bg_Omega_m];
free(pvecback);
class_call(background_free_noinput(pba), pba->error_message, errmsg);
......@@ -4214,34 +4262,40 @@ int input_prepare_pk_eq(
}
/** - restore cosmological parameters (w0, wa) to their true values before main call to CLASS modules */
pba->background_verbose = true_background_verbose;
pth->thermodynamics_verbose = true_thermodynamics_verbose;
pba->w0_fld = true_w0_fld;
pba->wa_fld = true_wa_fld;
/* in verbose mode, report the results */
if (input_verbose > 1) {
fprintf(stdout," Effective parameters for Pk-equal:\n");
fprintf(stdout," Effective parameters for Pk_equal:\n");
for (index_eq_z=0; index_eq_z<pnl->eq_tau_size; index_eq_z++) {
for (index_pk_eq_z=0; index_pk_eq_z<pnl->pk_eq_tau_size; index_pk_eq_z++) {
fprintf(stdout," * at z=%e, tau=%e w=%e Omega_m=%e\n",
z[index_eq_z],
pnl->eq_tau[index_eq_z],
pnl->eq_w_and_Omega[pnl->eq_size*index_eq_z+pnl->index_eq_w],
pnl->eq_w_and_Omega[pnl->eq_size*index_eq_z+pnl->index_eq_Omega_m]
z[index_pk_eq_z],
pnl->pk_eq_tau[index_pk_eq_z],
pnl->pk_eq_w_and_Omega[pnl->pk_eq_size*index_pk_eq_z+pnl->index_pk_eq_w],
pnl->pk_eq_w_and_Omega[pnl->pk_eq_size*index_pk_eq_z+pnl->index_pk_eq_Omega_m]
);
}
}
free(z);
/** - spline the table for later interpolation */
class_call(array_spline_table_lines(
pnl->eq_tau,
pnl->eq_tau_size,
pnl->eq_w_and_Omega,
pnl->eq_size,
pnl->eq_ddw_and_ddOmega,
pnl->pk_eq_tau,
pnl->pk_eq_tau_size,
pnl->pk_eq_w_and_Omega,
pnl->pk_eq_size,
pnl->pk_eq_ddw_and_ddOmega,
_SPLINE_NATURAL_,
errmsg),
errmsg,errmsg);
......
......@@ -180,7 +180,7 @@ int nonlinear_init(
class_alloc(lnk_l[index_pk],pnl->k_size*sizeof(double),pnl->error_message);//this is not really necessary
class_alloc(lnpk_l[index_pk],pnl->k_size*sizeof(double),pnl->error_message);
class_alloc(ddlnpk_l[index_pk],pnl->k_size*sizeof(double),pnl->error_message);
}
print_warning=_FALSE_;
......@@ -220,7 +220,7 @@ int nonlinear_init(
/* get P_NL(k) at this time */
if (print_warning == _FALSE_) {
class_call(nonlinear_halofit(ppr,
pba,
ppt,
......@@ -241,7 +241,7 @@ int nonlinear_init(
if (halofit_found_k_max == _TRUE_) {
// for debugging:
/*if ((index_tau == pnl->tau_size-1)){
/*if ((index_tau == pnl->tau_size-1)){
for (index_k=0; index_k<pnl->k_size; index_k++) {
fprintf(stdout,"%d %e %e %e\n",index_pk,pnl->k[index_k],pk_l[index_pk][index_k],pk_nl[index_pk][index_k]);
}
......@@ -285,7 +285,7 @@ int nonlinear_init(
}
}//end loop over pk_type
}//end loop over tau
for (index_pk=0; index_pk<pnl->pk_size; index_pk++){
free(pk_l[index_pk]);
......@@ -313,7 +313,7 @@ int nonlinear_free(
struct nonlinear *pnl
) {
int index_pk;
if (pnl->method > nl_none) {
if (pnl->method == nl_halofit) {
......@@ -329,9 +329,9 @@ int nonlinear_free(
}
if (pnl->has_pk_eq == _TRUE_) {
free(pnl->eq_tau);
free(pnl->eq_w_and_Omega);
free(pnl->eq_ddw_and_ddOmega);
free(pnl->pk_eq_tau);
free(pnl->pk_eq_w_and_Omega);
free(pnl->pk_eq_ddw_and_ddOmega);
}
return _SUCCESS_;
......@@ -358,7 +358,7 @@ int nonlinear_pk_l(
double source_ic1,source_ic2;
index_md = ppt->index_md_scalars;
// Initialize first, then assign correct value
index_delta = ppt->index_tp_delta_m;
if(index_pk == pnl->index_pk_m){
......@@ -506,7 +506,7 @@ int nonlinear_halofit(
double * w_and_Omega;
class_alloc(pvecback,pba->bg_size*sizeof(double),pnl->error_message);
Omega0_m = (pba->Omega0_cdm + pba->Omega0_b + pba->Omega0_ncdm_tot + pba->Omega0_dcdm);
//Initialize first, then assign correct value
......@@ -519,7 +519,7 @@ int nonlinear_halofit(
}
else {
class_stop(pnl->error_message,"P(k) is set neither to total matter nor to cold dark matter + baryons, index_pk=%d \n",index_pk);
}
}
if (pnl->has_pk_eq == _FALSE_) {
......@@ -538,33 +538,33 @@ int nonlinear_halofit(
}
else {
/* alternative method called PK-equal, described in 0810.0190 and
1601.0723, extending the range of validity of
HALOFIT from constant w to w0-wa models. In that
/* alternative method called Pk_equal, described in 0810.0190 and
1601.07230, extending the range of validity of
HALOFIT from constant w to (w0,wa) models. In that
case, some effective values of w0(tau_i) and
Omega_m(tau_i) have been pre-computed in the input
module, and we just ned to interpolate within
tabulated arrays, to get them at the current tau
value. */
Omega_m(tau_i) have been pre-computed in the
input module, and we just ned to interpolate
within tabulated arrays, to get them at the
current tau value. */
class_alloc(w_and_Omega,pnl->eq_size*sizeof(double),pnl->error_message);
class_alloc(w_and_Omega,pnl->pk_eq_size*sizeof(double),pnl->error_message);
class_call(array_interpolate_spline(
pnl->eq_tau,
pnl->eq_tau_size,
pnl->eq_w_and_Omega,
pnl->eq_ddw_and_ddOmega,
pnl->eq_size,
pnl->pk_eq_tau,
pnl->pk_eq_tau_size,
pnl->pk_eq_w_and_Omega,
pnl->pk_eq_ddw_and_ddOmega,
pnl->pk_eq_size,
tau,
&last_index,
w_and_Omega,
pnl->eq_size,
pnl->pk_eq_size,
pnl->error_message),
pnl->error_message,
pnl->error_message);
w0 = w_and_Omega[pnl->index_eq_w];
Omega_m = w_and_Omega[pnl->index_eq_Omega_m];
w0 = w_and_Omega[pnl->index_pk_eq_w];
Omega_m = w_and_Omega[pnl->index_pk_eq_Omega_m];
Omega_v = 1.-Omega_m;
free(w_and_Omega);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment