Commit 0bf59c07 authored by Julien Lesgourgues's avatar Julien Lesgourgues
Browse files

Removed few bugs pointed out in Munich workshop, most important was error in...

Removed few bugs pointed out in Munich workshop, most important was error in baryon drag redshift introduced in 2.1.0

* devel:
  updated version number to 2.1.2
  can ask for lensed polarisation without temperature
  removed obsolete warning for tensor polarisation
  fixed case r=0, and removed useless lines in background
  changed comment for halofit
  Update the model of Halofit used to match that of arXiv:1208.2701, which works better on small scales and incorporates a correction for dark energy.
  always write a # in headers
  fixed bug in clauclation of z_d
  when printing perturbation evolution, print also a
  corrected bug in print_variables: was printed only in last approx interval
  fixed bug appearing when tau_d is very small and z_d is above 10000
  get warning if pass file which is neither .ini nor .pre
  q smapling robust in range -0.5+0.5, and omegak now limited to this range
  z values in recfast derived in more robust way, avoiding compiler-dependent rounding issues
parents 2a867cda 5cb277ea
......@@ -445,6 +445,8 @@ extern "C" {
#define _TCMB_BIG_ 2.8 /**< maximal \f$ T_{cmb} \f$ in K */
#define _TCMB_SMALL_ 2.7 /**< minimal \f$ T_{cmb} \f$ in K */
#define _TOLERANCE_ON_CURVATURE_ 1.e-5 /**< if \f$ | \Omega_k | \f$ smaller than this, considered as flat */
#define _OMEGAK_BIG_ 0.5 /**< maximal \f$ Omega_k \f$ */
#define _OMEGAK_SMALL_ -0.5 /**< minimal \f$ Omega_k \f$ */
//@}
......
......@@ -15,7 +15,7 @@
#ifndef __COMMON__
#define __COMMON__
#define _VERSION_ "v2.1.1"
#define _VERSION_ "v2.1.2"
#define _TRUE_ 1 /**< integer associated to true statement */
#define _FALSE_ 0 /**< integer associated to false statement */
......
......@@ -242,7 +242,7 @@ int background_init(
/** - local variables : */
int n_ncdm;
double Omega0_tot,rho_ncdm_rel,rho_nu_rel;
double rho_ncdm_rel,rho_nu_rel;
int filenum=0;
/** - in verbose mode, provide some information */
......@@ -315,6 +315,12 @@ int background_init(
pba->error_message,
"T_cmb=%g out of bounds (%g<T_cmb<%g)",pba->T_cmb,_TCMB_SMALL_,_TCMB_BIG_);
/* H0 in Mpc^{-1} */
class_test((pba->Omega0_k < _OMEGAK_SMALL_)||(pba->Omega0_k > _OMEGAK_BIG_),
pba->error_message,
"Omegak = %g out of bounds (%g<Omegak<%g) \n",pba->Omega0_k,_OMEGAK_SMALL_,_OMEGAK_BIG_);
/* fluid equation of state */
if (pba->has_fld == _TRUE_) {
class_test(pba->w0_fld+pba->wa_fld>=1./3.,
pba->error_message,
......@@ -334,24 +340,6 @@ int background_init(
}
}
/* compute Omega0_tot and check curvature */
Omega0_tot = pba->Omega0_g + pba->Omega0_b;
if (pba->has_cdm == _TRUE_) {
Omega0_tot += pba->Omega0_cdm;
}
if (pba->has_ncdm == _TRUE_) {
Omega0_tot += pba->Omega0_ncdm_tot;
}
if (pba->has_lambda == _TRUE_) {
Omega0_tot += pba->Omega0_lambda;
}
if (pba->has_fld == _TRUE_) {
Omega0_tot += pba->Omega0_fld;
}
if (pba->has_ur == _TRUE_) {
Omega0_tot += pba->Omega0_ur;
}
/* check other quantities which would lead to segmentation fault if zero */
class_test(pba->a_today <= 0,
pba->error_message,
......
......@@ -68,12 +68,15 @@ int input_init_from_arguments(
"You have passed more than one input file with extension '.ini', choose one.");
strcpy(input_file,argv[i]);
}
if (strcmp(extension,".pre") == 0) {
else if (strcmp(extension,".pre") == 0) {
class_test(precision_file[0] != '\0',
errmsg,
"You have passed more than one precision with extension '.pre', choose one.");
strcpy(precision_file,argv[i]);
}
else {
fprintf(stdout,"Warning: the file %s has an extension different from .ini and .pre, so it has been ignored\n",argv[i]);
}
}
}
......@@ -1134,34 +1137,39 @@ int input_init(
class_read_double("r",ppm->r);
class_call(parser_read_string(pfc,"n_t",&string1,&flag1,errmsg),
errmsg,
errmsg);
if (ppm->r <= 0) {
ppt->has_tensors = _FALSE_;
}
else {
if (flag1 == _TRUE_) {
class_call(parser_read_string(pfc,"n_t",&string1,&flag1,errmsg),
errmsg,
errmsg);
if ((strstr(string1,"SCC") != NULL) || (strstr(string1,"scc") != NULL)) {
ppm->n_t = -ppm->r/8.*(2.-ppm->r/8.-ppm->n_s);
}
else {
class_read_double("n_t",ppm->n_t);
}
if (flag1 == _TRUE_) {
}
if ((strstr(string1,"SCC") != NULL) || (strstr(string1,"scc") != NULL)) {
ppm->n_t = -ppm->r/8.*(2.-ppm->r/8.-ppm->n_s);
}
else {
class_read_double("n_t",ppm->n_t);
}
class_call(parser_read_string(pfc,"alpha_t",&string1,&flag1,errmsg),
errmsg,
errmsg);
}
if (flag1 == _TRUE_) {
class_call(parser_read_string(pfc,"alpha_t",&string1,&flag1,errmsg),
errmsg,
errmsg);
if ((strstr(string1,"SCC") != NULL) || (strstr(string1,"scc") != NULL)) {
ppm->alpha_t = ppm->r/8.*(ppm->r/8.+ppm->n_s-1.);
}
else {
class_read_double("alpha_t",ppm->alpha_t);
}
if (flag1 == _TRUE_) {
if ((strstr(string1,"SCC") != NULL) || (strstr(string1,"scc") != NULL)) {
ppm->alpha_t = ppm->r/8.*(ppm->r/8.+ppm->n_s-1.);
}
else {
class_read_double("alpha_t",ppm->alpha_t);
}
}
}
}
}
......
......@@ -691,13 +691,15 @@ int lensing_init(
/** - compute lensed Cls by integration */
//debut = omp_get_wtime();
class_call(lensing_lensed_cl_tt(ksi,d00,w8,num_mu-1,ple),
ple->error_message,
ple->error_message);
if (ppr->accurate_lensing == _FALSE_) {
class_call(lensing_addback_cl_tt(ple,cl_tt),
if (ple->has_tt==_TRUE_) {
class_call(lensing_lensed_cl_tt(ksi,d00,w8,num_mu-1,ple),
ple->error_message,
ple->error_message);
if (ppr->accurate_lensing == _FALSE_) {
class_call(lensing_addback_cl_tt(ple,cl_tt),
ple->error_message,
ple->error_message);
}
}
if (ple->has_te==_TRUE_) {
......
......@@ -75,7 +75,7 @@ int nonlinear_init(
else if (pnl->method == nl_halofit) {
if (pnl->nonlinear_verbose > 0)
printf("Computing non-linear matter power spectrum with Halofit (including update by Bird et al 2011)\n");
printf("Computing non-linear matter power spectrum with Halofit (including update Takahashi et al. 2012 and Bird 2014)\n");
if (pba->has_ncdm) {
for (index_ncdm=0;index_ncdm < pba->N_ncdm; index_ncdm++){
......@@ -274,7 +274,7 @@ int nonlinear_halofit(
double *k_nl
) {
double Omega_m,Omega_v,fnu,Omega0_m;
double Omega_m,Omega_v,fnu,Omega0_m, w0;
/** determine non linear ratios (from pk) **/
......@@ -283,7 +283,7 @@ int nonlinear_halofit(
double sigma,rknl,rneff,rncur,d1,d2;
double diff,xlogr1,xlogr2,rmid;
double extragam,gam,a,b,c,xmu,xnu,alpha,beta,f1,f2,f3;
double gam,a,b,c,xmu,xnu,alpha,beta,f1,f2,f3;
double pk_linaa;
double y;
double f1a,f2a,f3a,f1b,f2b,f3b,frac;
......@@ -302,6 +302,7 @@ int nonlinear_halofit(
class_alloc(pvecback,pba->bg_size*sizeof(double),pnl->error_message);
Omega0_m = (pba->Omega0_cdm + pba->Omega0_b + pba->Omega0_ncdm_tot);
w0 = pba->w0_fld;
fnu = pba->Omega0_ncdm_tot/(pba->Omega0_b + pba->Omega0_cdm+pba->Omega0_ncdm_tot);
anorm = 1./(2*pow(_PI_,2));
......@@ -454,18 +455,19 @@ int nonlinear_halofit(
* scales by a factor of two. Add an extra correction from the
* simulations in Bird, Viel,Haehnelt 2011 which partially accounts for
* this.*/
extragam = 0.3159 -0.0765*rneff -0.8350*rncur;
gam=extragam+0.86485+0.2989*rneff+0.1631*rncur;
a=1.4861+1.83693*rneff+1.67618*rneff*rneff+0.7940*rneff*rneff*rneff+0.1670756*rneff*rneff*rneff*rneff-0.620695*rncur;
/*SPB14: This version of halofit is an updated version of the fit to the massive neutrinos
* based on the results of Takahashi 2012, (arXiv:1208.2701).
*/
gam=0.1971-0.0843*rneff+0.8460*rncur;
a=1.5222+2.8553*rneff+2.3706*rneff*rneff+0.9903*rneff*rneff*rneff+ 0.2250*rneff*rneff*rneff*rneff-0.6038*rncur+0.1749*Omega_v*(1.+w0);
a=pow(10,a);
b=pow(10,(0.9463+0.9466*rneff+0.3084*rneff*rneff-0.940*rncur));
c=pow(10,(-0.2807+0.6669*rneff+0.3214*rneff*rneff-0.0793*rncur));
xmu = pow(10,(-3.54419+0.19086*rneff));
xnu = pow(10,(0.95897+1.2857*rneff));
alpha = 1.38848+0.3701*rneff-0.1452*rneff*rneff;
beta = 0.8291+0.9854*rneff+0.3400*pow(rneff,2)+fnu*(-6.4868+1.4373*pow(rneff,2));
b=pow(10, (-0.5642+0.5864*rneff+0.5716*rneff*rneff-1.5474*rncur+0.2279*Omega_v*(1.+w0)));
c=pow(10, 0.3698+2.0404*rneff+0.8161*rneff*rneff+0.5869*rncur);
xmu=0.;
xnu=pow(10,5.2105+3.6902*rneff);
alpha=abs(6.0835+1.3373*rneff-0.1959*rneff*rneff-5.5274*rncur);
beta=2.0379-0.7354*rneff+0.3157*pow(rneff,2)+1.2490*pow(rneff,3)+0.3980*pow(rneff,4)-0.1682*rncur + fnu*(1.081 + 0.395*pow(rneff,2));
if(fabs(1-Omega_m)>0.01) { /*then omega evolution */
f1a=pow(Omega_m,(-0.0732));
f2a=pow(Omega_m,(-0.1423));
......@@ -486,8 +488,8 @@ int nonlinear_halofit(
y=(rk/rknl);
pk_halo = a*pow(y,f1*3.)/(1.+b*pow(y,f2)+pow(f3*c*y,3.-gam));
pk_halo=pk_halo/(1+xmu*pow(y,-1)+xnu*pow(y,-2))*(1+fnu*(2.080-12.39*(Omega0_m-0.3))/(1+1.201e-03*pow(y,3)));
pk_linaa=pk_lin*(1+fnu*26.29*pow(rk,2)/(1+1.5*pow(rk,2)));
pk_halo=pk_halo/(1+xmu*pow(y,-1)+xnu*pow(y,-2))*(1+fnu*(0.977-18.015*(Omega0_m-0.3)));
pk_linaa=pk_lin*(1+fnu*47.48*pow(rk,2)/(1+1.5*pow(rk,2)));
pk_quasi=pk_lin*pow((1+pk_linaa),beta)/(1+pk_linaa*alpha)*exp(-y/4.0-pow(y,2)/8.0);
pk_nl[index_k] = (pk_halo+pk_quasi)/pow(pnl->k[index_k],3)/anorm;
......
......@@ -1764,6 +1764,7 @@ int output_open_background_file(
fprintf(*backfile,"# All densities are mutiplied by (8piG/3) (below, shortcut notation (.) for this factor) \n");
/** Length of the columntitle should be less than _OUTPUTPRECISION_+6 to be indented correctly,
but it can be as long as . */
fprintf(*backfile,"#");
class_fprintf_columntitle(*backfile,"z",_TRUE_);
class_fprintf_columntitle(*backfile,"proper time [Gyr]",_TRUE_);
class_fprintf_columntitle(*backfile,"conf. time * c [Mpc]",_TRUE_);
......@@ -1808,7 +1809,7 @@ int output_one_line_of_background(
int n;
fprintf(backfile," ");
class_fprintf_double(backfile,pba->a_today/pvecback[pba->index_bg_a]-1.,_TRUE_);
class_fprintf_double(backfile,pvecback[pba->index_bg_time]/_Gyr_over_Mpc_,_TRUE_);
class_fprintf_double(backfile,pba->conformal_age-pvecback[pba->index_bg_conf_distance],_TRUE_);
......@@ -1869,6 +1870,7 @@ int output_open_thermodynamics_file(
/** Length of the columntitle should be less than _OUTPUTPRECISION_+6 to be indented correctly,
but it can be as long as _COLUMNTITLE_. */
fprintf(*thermofile,"#");
class_fprintf_columntitle(*thermofile,"z",_TRUE_);
class_fprintf_columntitle(*thermofile,"conf. time [Mpc]",_TRUE_);
class_fprintf_columntitle(*thermofile,"x_e",_TRUE_);
......@@ -1908,6 +1910,7 @@ int output_one_line_of_thermodynamics(
double * pvecthermo
) {
fprintf(thermofile," ");
class_fprintf_double(thermofile,z,_TRUE_);
class_fprintf_double(thermofile,tau,_TRUE_);
class_fprintf_double(thermofile,pvecthermo[pth->index_th_xe],_TRUE_);
......
......@@ -199,11 +199,6 @@ int perturb_init(
if ((ppt->has_tensors == _TRUE_) && (ppt->perturbations_verbose > 0))
printf("Warning: neglect coupling between gravity waves and neutrino shear:\n leads to typically 40 per cent error, but only for l > 150 (so irrelevant for current experiments)\n");
if ((ppt->has_cl_cmb_polarization == _TRUE_) &&
(ppt->has_tensors == _TRUE_) && (ppt->perturbations_verbose > 0)) {
printf("Warning: yet unexplained difference with CAMB for impact of reionization on polarized tensors,\n affecting ClEE, BB, TE for typically l<10\n");
}
/** - initialize all indices and lists in perturbs structure using perturb_indices_of_perturbs() */
class_call(perturb_indices_of_perturbs(ppr,
......@@ -1825,6 +1820,18 @@ int perturb_solve(
ppaw.ppw = ppw;
ppaw.ppw->inter_mode = pba->inter_closeby;
/** - check whether we need to print perturbations to a file for this wavenumber */
perhaps_print_variables = NULL;
for (index_ikout=0; index_ikout<ppt->k_output_values_num; index_ikout++){
if (ppt->index_k_output_values[index_ikout] == index_k){
perhaps_print_variables = perturb_print_variables;
class_call(perturb_prepare_output_file(pba,ppt,ppw,index_ikout,index_md),
ppt->error_message,
ppt->error_message);
}
}
/** - loop over intervals over which approximatiomn scheme is uniform. For each interval: */
for (index_interval=0; index_interval<interval_number; index_interval++) {
......@@ -1876,16 +1883,6 @@ int perturb_solve(
generic_evolver = evolver_ndf15;
}
perhaps_print_variables = NULL;
for (index_ikout=0; index_ikout<ppt->k_output_values_num; index_ikout++){
if (ppt->index_k_output_values[index_ikout] == index_k){
perhaps_print_variables = perturb_print_variables;
class_call(perturb_prepare_output_file(pba,ppt,ppw,index_ikout,index_md),
ppt->error_message,
ppt->error_message);
}
}
class_call(generic_evolver(perturb_derivs,
interval_limit[index_interval],
interval_limit[index_interval+1],
......@@ -1905,11 +1902,13 @@ int perturb_solve(
ppt->error_message,
ppt->error_message);
if (perhaps_print_variables != NULL)
fclose(ppw->perturb_output_file);
}
/** - if perturbations were printed in a file, close the file */
if (perhaps_print_variables != NULL)
fclose(ppw->perturb_output_file);
/** fill the source terms array with zeros for all times between
then last integrated time tau_max and tau_today. */
......@@ -1962,7 +1961,9 @@ int perturb_prepare_output_file(struct background * pba,
fprintf(ppw->perturb_output_file,
"#scalar perturbations for mode k = %.*e Mpc^(-1)\n",
_OUTPUTPRECISION_,k);
class_fprintf_columntitle(ppw->perturb_output_file,"tau [Mpc^-1]",_TRUE_);
fprintf(ppw->perturb_output_file,"#");
class_fprintf_columntitle(ppw->perturb_output_file,"tau [Mpc]",_TRUE_);
class_fprintf_columntitle(ppw->perturb_output_file,"a",_TRUE_);
class_fprintf_columntitle(ppw->perturb_output_file,"delta_g",_TRUE_);
class_fprintf_columntitle(ppw->perturb_output_file,"theta_g",_TRUE_);
class_fprintf_columntitle(ppw->perturb_output_file,"shear_g",_TRUE_);
......@@ -2001,7 +2002,9 @@ int perturb_prepare_output_file(struct background * pba,
fprintf(ppw->perturb_output_file,
"#tensor perturbations for mode k = %.*e Mpc^(-1)\n",
_OUTPUTPRECISION_,k);
class_fprintf_columntitle(ppw->perturb_output_file,"tau [Mpc^-1]",_TRUE_);
fprintf(ppw->perturb_output_file,"#");
class_fprintf_columntitle(ppw->perturb_output_file,"tau [Mpc]",_TRUE_);
class_fprintf_columntitle(ppw->perturb_output_file,"a",_TRUE_);
class_fprintf_columntitle(ppw->perturb_output_file,"delta_g",_TRUE_);
class_fprintf_columntitle(ppw->perturb_output_file,"shear_g",_TRUE_);
class_fprintf_columntitle(ppw->perturb_output_file,"l4_g",_TRUE_);
......@@ -5194,7 +5197,9 @@ int perturb_print_variables(double tau,
}
fprintf(ppw->perturb_output_file," ");
class_fprintf_double(ppw->perturb_output_file, tau, _TRUE_);
class_fprintf_double(ppw->perturb_output_file, pvecback[pba->index_bg_a], _TRUE_);
class_fprintf_double(ppw->perturb_output_file, delta_g, _TRUE_);
class_fprintf_double(ppw->perturb_output_file, theta_g, _TRUE_);
class_fprintf_double(ppw->perturb_output_file, shear_g, _TRUE_);
......@@ -5257,7 +5262,9 @@ int perturb_print_variables(double tau,
pol4_g = 0.;
}
fprintf(ppw->perturb_output_file," ");
class_fprintf_double(ppw->perturb_output_file, tau, _TRUE_);
class_fprintf_double(ppw->perturb_output_file, pvecback[pba->index_bg_a], _TRUE_);
class_fprintf_double(ppw->perturb_output_file, delta_g, _TRUE_);
class_fprintf_double(ppw->perturb_output_file, shear_g, _TRUE_);
class_fprintf_double(ppw->perturb_output_file, l4_g, _TRUE_);
......
......@@ -413,7 +413,7 @@ int thermodynamics_init(
class_call(background_at_tau(pba,
tau_table[index_tau],
pba->short_info,
pba->normal_info,
pba->inter_closeby,
&last_index_back,
pvecback),
......@@ -663,7 +663,7 @@ int thermodynamics_init(
interpolation) and sound horizon at that time */
index_tau=0;
while (pth->thermodynamics_table[(index_tau)*pth->th_size+pth->index_th_tau_d] < 1.)
while ((pth->thermodynamics_table[(index_tau)*pth->th_size+pth->index_th_tau_d] < 1.) && (index_tau < pth->tt_size))
index_tau++;
pth->z_d = pth->z_table[index_tau-1]+
......@@ -2562,8 +2562,9 @@ int thermodynamics_recombination_with_recfast(
for(i=0; i <Nz; i++) {
zstart = zinitial * (1. - (double)i / (double)Nz);
zend = zinitial * (1. - (double)(i+1) / (double)Nz);
zstart = zinitial * (double)(Nz-i) / (double)Nz;
zend = zinitial * (double)(Nz-i-1) / (double)Nz;
z = zend;
/** -> first approximation: H and Helium fully ionized */
......
......@@ -1153,10 +1153,10 @@ int transfer_get_q_list(
/* max contribution from non-integer nu values */
q_step = 1.+q_period*ppr->q_logstep_spline;
q_size_max = 2*(int)(log(q_max/q_min)/log(q_step));
q_size_max = 5*(int)(log(q_max/q_min)/log(q_step));
q_step = q_period*ppr->q_linstep;
q_size_max += 2*(int)((q_max-q_min)/q_step);
q_size_max += 5*(int)((q_max-q_min)/q_step);
}
......
/** @file class.c
* Julien Lesgourgues, 17.04.2011
/** @file class.c
* Julien Lesgourgues, 17.04.2011
*/
/* this main runs only the background, thermodynamics and perturbation part */
......@@ -21,7 +21,7 @@ int main(int argc, char **argv) {
ErrorMsg errmsg; /* for error messages */
if (input_init_from_arguments(argc, argv,&pr,&ba,&th,&pt,&tr,&pm,&sp,&nl,&le,&op,errmsg) == _FAILURE_) {
printf("\n\nError running input_init_from_arguments \n=>%s\n",errmsg);
printf("\n\nError running input_init_from_arguments \n=>%s\n",errmsg);
return _FAILURE_;
}
......@@ -40,40 +40,44 @@ int main(int argc, char **argv) {
return _FAILURE_;
}
/*********************************************************************/
/* here you can output the source function S(k,tau) of your choice */
/*********************************************************************/
if (pt.has_perturbations == _TRUE_) {
FILE * output;
int index_k,index_tau;
/*********************************************************************/
/* here you can output the source function S(k,tau) of your choice */
/*********************************************************************/
/* choose a mode (scalar, tensor, ...) */
int index_mode=pt.index_md_scalars;
FILE * output;
int index_k,index_tau;
/* choose a type (temperature, polarization, grav. pot., ...) */
int index_type=pt.index_tp_t0;
/* choose a mode (scalar, tensor, ...) */
int index_mode=pt.index_md_scalars;
/* choose an initial condition (ad, bi, cdi, nid, niv, ...) */
int index_ic=pt.index_ic_ad;
/* choose a type (temperature, polarization, grav. pot., ...) */
int index_type=pt.index_tp_t0;
output=fopen("output/source.dat","w");
fprintf(output,"# k tau S\n");
/* choose an initial condition (ad, bi, cdi, nid, niv, ...) */
int index_ic=pt.index_ic_ad;
for (index_k=0; index_k < pt.k_size; index_k++) {
for (index_tau=0; index_tau < pt.tau_size; index_tau++) {
fprintf(output,"%e %e %e\n",
pt.k[index_k],
pt.tau_sampling[index_tau],
pt.sources[index_mode]
[index_ic * pt.tp_size[index_mode] + index_type]
[index_tau * pt.k_size + index_k]
);
output=fopen("output/source.dat","w");
fprintf(output,"# k tau S\n");
for (index_k=0; index_k < pt.k_size; index_k++) {
for (index_tau=0; index_tau < pt.tau_size; index_tau++) {
fprintf(output,"%e %e %e\n",
pt.k[index_k],
pt.tau_sampling[index_tau],
pt.sources[index_mode]
[index_ic * pt.tp_size[index_mode] + index_type]
[index_tau * pt.k_size + index_k]
);
}
fprintf(output,"\n");
}
fprintf(output,"\n");
fclose(output);
}
fclose(output);
/****** all calculations done, now free the structures ******/
......
Markdown is supported
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