Commit 54f22882 authored by Julien Lesgourgues's avatar Julien Lesgourgues
Browse files

Few bug fixes in thermodynamics and transfer

* 2.4:
  change version number to 2.4.4
  rewritten small part related to annihilation to avoid rare bug
  slightly revised the expression and sampling of the jumps in the reio_bins_tanh scheme, to make it more robust and avoid rare interpolation errors
  fixed two bugs concerning reio_bins_tanh
  fixed bug in total velocity transfer function (rho_cdm contribution in sync was missing)
  fixed bug when using hyrec + annihilation
parents 129d024b c4525f19
......@@ -103,8 +103,14 @@ double rec_Tmss(double xe, double Tr, double H, double fHe, double nH, double en
double chi_heat;
//chi_heat = (1.+2.*xe)/3.; // old approximation from Chen and Kamionkowski
chi_heat = 0.996857*(1.-pow(1.-pow(xe,0.300134),1.51035)); // coefficient as revised by Galli et al. 2013 (in fact it is a fit by Vivian Poulin of columns 1 and 2 in Table V of Galli et al. 2013)
if (chi_heat > 1.) chi_heat = 1.;
// coefficient as revised by Galli et al. 2013 (in fact it is a fit by Vivian Poulin of columns 1 and 2 in Table V of Galli et al. 2013)
if (xe < 1.) {
chi_heat = 0.996857*(1.-pow(1.-pow(xe,0.300134),1.51035));
if (chi_heat > 1.) chi_heat = 1.;
}
else
chi_heat = 1.;
return Tr/(1.+H/4.91466895548409e-22/Tr/Tr/Tr/Tr*(1.+xe+fHe)/xe)
+2./3./kBoltz*chi_heat/nH*energy_rate/(4.91466895548409e-22*pow(Tr,4)*xe);
......@@ -121,8 +127,14 @@ double rec_dTmdlna(double xe, double Tm, double Tr, double H, double fHe, double
double chi_heat;
//chi_heat = (1.+2.*xe)/3.; // old approximation from Chen and Kamionkowski
chi_heat = 0.996857*(1.-pow(1.-pow(xe,0.300134),1.51035)); // coefficient as revised by Galli et al. 2013 (in fact it is a fit by Vivian Poulin of columns 1 and 2 in Table V of Galli et al. 2013)
if (chi_heat > 1.) chi_heat = 1.;
// coefficient as revised by Galli et al. 2013 (in fact it is a fit by Vivian Poulin of columns 1 and 2 in Table V of Galli et al. 2013)
if (xe < 1.) {
chi_heat = 0.996857*(1.-pow(1.-pow(xe,0.300134),1.51035));
if (chi_heat > 1.) chi_heat = 1.;
}
else
chi_heat = 1.;
return -2.*Tm + 4.91466895548409e-22*Tr*Tr*Tr*Tr*xe/(1.+xe+fHe)*(Tr-Tm)/H
+2./3./kBoltz*chi_heat/nH*energy_rate/(1.+xe+fHe)/H;
......
......@@ -15,7 +15,7 @@
#ifndef __COMMON__
#define __COMMON__
#define _VERSION_ "v2.4.3"
#define _VERSION_ "v2.4.4"
#define _TRUE_ 1 /**< integer associated to true statement */
#define _FALSE_ 0 /**< integer associated to false statement */
......
......@@ -4699,6 +4699,14 @@ int perturb_approximations(
/** (b.2.a) compute recombination time scale for photons, \f$ \tau_{\gamma} = 1/ \kappa' \f$ */
tau_c = 1./ppw->pvecthermo[pth->index_th_dkappa];
class_test(tau_c < 0.,
ppt->error_message,
"tau_c = 1/kappa' should always be positive unless there is something wrong in the thermodynamics module. However you have here tau_c=%e at z=%e, conformal time=%e x_e=%e. (This could come from the interpolation of a too poorly sampled reionisation history?).\n",
tau_c,
1./ppw->pvecback[pba->index_bg_a]-1.,
tau,
ppw->pvecthermo[pth->index_th_xe]);
/** (b.2.b) check whether tight-coupling approximation should be on */
if ((tau_c/tau_h < ppr->tight_coupling_trigger_tau_c_over_tau_h) &&
(tau_c/tau_k < ppr->tight_coupling_trigger_tau_c_over_tau_k)) {
......
......@@ -2915,10 +2915,10 @@ int spectra_matter_transfers(
delta_rho_tot += rho_i * delta_i;
rho_tot += rho_i;
}
rho_tot += rho_i;
if (ppt->has_source_theta_b == _TRUE_) {
theta_i = ppt->sources[index_md]
......@@ -2929,10 +2929,10 @@ int spectra_matter_transfers(
rho_plus_p_theta_tot += rho_i * theta_i;
rho_plus_p_tot += rho_i;
}
rho_plus_p_tot += rho_i;
/* T_cdm(k,tau) */
if (pba->has_cdm == _TRUE_) {
......@@ -2949,10 +2949,10 @@ int spectra_matter_transfers(
delta_rho_tot += rho_i * delta_i;
rho_tot += rho_i;
}
rho_tot += rho_i;
if (ppt->has_source_theta_cdm == _TRUE_) {
theta_i = ppt->sources[index_md]
......@@ -2963,10 +2963,10 @@ int spectra_matter_transfers(
rho_plus_p_theta_tot += rho_i * theta_i;
rho_plus_p_tot += rho_i;
}
rho_plus_p_tot += rho_i;
}
/* T_dcdm(k,tau) */
......@@ -2985,10 +2985,10 @@ int spectra_matter_transfers(
delta_rho_tot += rho_i * delta_i;
rho_tot += rho_i;
}
rho_tot += rho_i;
if (ppt->has_source_theta_dcdm == _TRUE_) {
theta_i = ppt->sources[index_md]
......@@ -2999,10 +2999,10 @@ int spectra_matter_transfers(
rho_plus_p_theta_tot += rho_i * theta_i;
rho_plus_p_tot += rho_i;
}
rho_plus_p_tot += rho_i;
}
/* T_scf(k,tau) */
......@@ -3021,9 +3021,10 @@ int spectra_matter_transfers(
delta_rho_tot += rho_i * delta_i;
rho_tot += rho_i;
}
rho_tot += rho_i;
if (ppt->has_source_theta_scf == _TRUE_) {
theta_i = ppt->sources[index_md]
......@@ -3034,9 +3035,10 @@ int spectra_matter_transfers(
rho_plus_p_theta_tot += (rho_i + pvecback_sp_long[pba->index_bg_p_scf]) * theta_i;
rho_plus_p_tot += (rho_i + pvecback_sp_long[pba->index_bg_p_scf]);
}
rho_plus_p_tot += (rho_i + pvecback_sp_long[pba->index_bg_p_scf]);
}
......@@ -3056,10 +3058,10 @@ int spectra_matter_transfers(
delta_rho_tot += rho_i * delta_i;
rho_tot += rho_i;
}
rho_tot += rho_i;
if (ppt->has_source_theta_fld == _TRUE_) {
theta_i = ppt->sources[index_md]
......@@ -3070,10 +3072,10 @@ int spectra_matter_transfers(
rho_plus_p_theta_tot += (1. + pba->w0_fld + pba->wa_fld * (1. - pvecback_sp_long[pba->index_bg_a] / pba->a_today)) * rho_i * theta_i;
rho_plus_p_tot += (1. + pba->w0_fld + pba->wa_fld * (1. - pvecback_sp_long[pba->index_bg_a] / pba->a_today)) * rho_i;
}
rho_plus_p_tot += (1. + pba->w0_fld + pba->wa_fld * (1. - pvecback_sp_long[pba->index_bg_a] / pba->a_today)) * rho_i;
}
/* T_ur(k,tau) */
......@@ -3092,10 +3094,10 @@ int spectra_matter_transfers(
delta_rho_tot += rho_i * delta_i;
rho_tot += rho_i;
}
rho_tot += rho_i;
if (ppt->has_source_theta_ur == _TRUE_) {
theta_i = ppt->sources[index_md]
......@@ -3106,10 +3108,10 @@ int spectra_matter_transfers(
rho_plus_p_theta_tot += 4./3. * rho_i * theta_i;
rho_plus_p_tot += 4./3. * rho_i;
}
rho_plus_p_tot += 4./3. * rho_i;
}
/* T_dr(k,tau) */
......@@ -3128,10 +3130,10 @@ int spectra_matter_transfers(
delta_rho_tot += rho_i * delta_i;
rho_tot += rho_i;
}
rho_tot += rho_i;
if (ppt->has_source_theta_dr == _TRUE_) {
theta_i = ppt->sources[index_md]
......@@ -3142,10 +3144,10 @@ int spectra_matter_transfers(
rho_plus_p_theta_tot += 4./3. * rho_i * theta_i;
rho_plus_p_tot += 4./3. * rho_i;
}
rho_plus_p_tot += 4./3. * rho_i;
}
/* T_ncdm_i(k,tau) */
......@@ -3166,9 +3168,10 @@ int spectra_matter_transfers(
delta_rho_tot += rho_i * delta_i;
rho_tot += rho_i;
}
rho_tot += rho_i;
if (ppt->has_source_theta_ncdm == _TRUE_) {
theta_i = ppt->sources[index_md]
......@@ -3179,9 +3182,10 @@ int spectra_matter_transfers(
rho_plus_p_theta_tot += (rho_i + pvecback_sp_long[pba->index_bg_p_ncdm1+n_ncdm]) * theta_i;
rho_plus_p_tot += (rho_i + pvecback_sp_long[pba->index_bg_p_ncdm1+n_ncdm]);
}
rho_plus_p_tot += (rho_i + pvecback_sp_long[pba->index_bg_p_ncdm1+n_ncdm]);
}
}
......
......@@ -1322,6 +1322,7 @@ int thermodynamics_reionization_function(
/** - define local variables */
double argument;
int i;
double z_jump;
/** - implementation of ionization function similar to the one in CAMB */
......@@ -1396,6 +1397,13 @@ int thermodynamics_reionization_function(
i = 0;
while (preio->reionization_parameters[preio->index_reio_first_z+i+1]<z) i++;
/* This is the expression of the tanh-like jumps of the
reio_bins_tanh scheme until the 10.06.2015. It appeared to be
not robust enough. It could lead to a kink in xe(z) near the
maximum value of z at which reionisation is sampled. It has
been replaced by the simpler and more robust expression
below.
*xe = preio->reionization_parameters[preio->index_reio_first_xe+i]
+0.5*(tanh((2.*(z-preio->reionization_parameters[preio->index_reio_first_z+i])
/(preio->reionization_parameters[preio->index_reio_first_z+i+1]
......@@ -1404,7 +1412,27 @@ int thermodynamics_reionization_function(
/tanh(1./preio->reionization_parameters[preio->index_reio_step_sharpness])+1.)
*(preio->reionization_parameters[preio->index_reio_first_xe+i+1]
-preio->reionization_parameters[preio->index_reio_first_xe+i]);
*/
/* compute the central redshift value of the tanh jump */
if (i == preio->reio_num_z-2) {
z_jump = preio->reionization_parameters[preio->index_reio_first_z+i]
+ 0.5*(preio->reionization_parameters[preio->index_reio_first_z+i]
-preio->reionization_parameters[preio->index_reio_first_z+i-1]);
}
else {
z_jump = 0.5*(preio->reionization_parameters[preio->index_reio_first_z+i+1]
+ preio->reionization_parameters[preio->index_reio_first_z+i]);
}
/* implementation of the tanh jump */
*xe = preio->reionization_parameters[preio->index_reio_first_xe+i]
+0.5*(tanh((z-z_jump)
/preio->reionization_parameters[preio->index_reio_step_sharpness])+1.)
*(preio->reionization_parameters[preio->index_reio_first_xe+i+1]
-preio->reionization_parameters[preio->index_reio_first_xe+i]);
}
......@@ -1689,9 +1717,11 @@ int thermodynamics_reionization(
/* check that this input can be interpreted by the code */
for (bin=1; bin<pth->binned_reio_num; bin++) {
class_test(pth->binned_reio_z[bin]<pth->binned_reio_z[bin],
class_test(pth->binned_reio_z[bin-1]>=pth->binned_reio_z[bin],
pth->error_message,
"value of reionization bin centers z_i expected to be passed in growing order");
"value of reionization bin centers z_i expected to be passed in growing order: %e, %e",
pth->binned_reio_z[bin-1],
pth->binned_reio_z[bin]);
}
/* the code will not only copy here the "bin centers" passed in
......@@ -1705,11 +1735,12 @@ int thermodynamics_reionization(
/* find largest value of z in the array. We choose to define it as
z_(i_max) + (the distance between z_(i_max) and z_(i_max-1)). E.g. if
the bins are in 10,12,14, the largest z will be 16. */
z_(i_max) + 2*(the distance between z_(i_max) and z_(i_max-1)). E.g. if
the bins are in 10,12,14, the largest z will be 18. */
preio->reionization_parameters[preio->index_reio_first_z+preio->reio_num_z-1] =
2.*preio->reionization_parameters[preio->index_reio_first_z+preio->reio_num_z-2]
-preio->reionization_parameters[preio->index_reio_first_z+preio->reio_num_z-3];
preio->reionization_parameters[preio->index_reio_first_z+preio->reio_num_z-2]
+2.*(preio->reionization_parameters[preio->index_reio_first_z+preio->reio_num_z-2]
-preio->reionization_parameters[preio->index_reio_first_z+preio->reio_num_z-3]);
/* copy this value in reio_start */
preio->reionization_parameters[preio->index_reio_start] = preio->reionization_parameters[preio->index_reio_first_z+preio->reio_num_z-1];
......@@ -1730,9 +1761,15 @@ int thermodynamics_reionization(
-preio->reionization_parameters[preio->index_reio_first_z+2];
/* check it's not too small */
/* 6.06.2015: changed this test to simply imposing that the first z is at least zero */
/*
class_test(preio->reionization_parameters[preio->index_reio_first_z] < 0,
pth->error_message,
"final redshift for reionization = %e, you must change the binning or redefine the way in which the code extrapolates below the first value of z_i",preio->reionization_parameters[preio->index_reio_first_z]);
*/
if (preio->reionization_parameters[preio->index_reio_first_z] < 0) {
preio->reionization_parameters[preio->index_reio_first_z] = 0.;
}
/* infer xe before reio */
class_call(thermodynamics_get_xe_before_reionization(ppr,
......@@ -2034,19 +2071,33 @@ int thermodynamics_reionization_sample(
/** - derivative of baryon temperature */
class_call(thermodynamics_energy_injection(ppr,pba,preco,z,&energy_rate,pth->error_message),
pth->error_message,
pth->error_message);
//chi_heat = (1.+2.*preio->reionization_table[i*preio->re_size+preio->index_re_xe])/3.; // old approximation from Chen and Kamionkowski
chi_heat = MIN(0.996857*(1.-pow(1.-pow(preio->reionization_table[i*preio->re_size+preio->index_re_xe],0.300134),1.51035)),1); // coefficient as revised by Slatyer et al. 2013 (in fact it is a fit by Vivian Poulin of columns 1 and 2 in Table V of Slatyer et al. 2013)
dTdz=2./(1+z)*preio->reionization_table[i*preio->re_size+preio->index_re_Tb]
-2.*mu/_m_e_*4.*pvecback[pba->index_bg_rho_g]/3./pvecback[pba->index_bg_rho_b]*opacity*
(pba->T_cmb * (1.+z)-preio->reionization_table[i*preio->re_size+preio->index_re_Tb])/pvecback[pba->index_bg_H]
-2./(3.*_k_B_)*energy_rate*chi_heat
/(preco->Nnow*pow(1.+z,3))/(1.+preco->fHe+preio->reionization_table[i*preio->re_size+preio->index_re_xe])
/(pvecback[pba->index_bg_H]*_c_/_Mpc_over_m_*(1.+z)); /* energy injection */
(pba->T_cmb * (1.+z)-preio->reionization_table[i*preio->re_size+preio->index_re_Tb])/pvecback[pba->index_bg_H];
if (preco->annihilation > 0) {
class_call(thermodynamics_energy_injection(ppr,pba,preco,z,&energy_rate,pth->error_message),
pth->error_message,
pth->error_message);
// old approximation from Chen and Kamionkowski:
// chi_heat = (1.+2.*preio->reionization_table[i*preio->re_size+preio->index_re_xe])/3.;
// coefficient as revised by Slatyer et al. 2013
// (in fact it is a fit by Vivian Poulin of columns 1 and 2 in Table V
// of Slatyer et al. 2013):
xe = preio->reionization_table[i*preio->re_size+preio->index_re_xe];
if (xe < 1.)
chi_heat = MIN(0.996857*(1.-pow(1.-pow(xe,0.300134),1.51035)),1);
else
chi_heat = 1.;
dTdz+= -2./(3.*_k_B_)*energy_rate*chi_heat
/(preco->Nnow*pow(1.+z,3))/(1.+preco->fHe+preio->reionization_table[i*preio->re_size+preio->index_re_xe])
/(pvecback[pba->index_bg_H]*_c_/_Mpc_over_m_*(1.+z)); /* energy injection */
}
/** - increment baryon temperature */
......@@ -3079,7 +3130,12 @@ int thermodynamics_derivs_with_recfast(
/* equations modified to take into account energy injection from dark matter */
//chi_heat = (1.+2.*preio->reionization_table[i*preio->re_size+preio->index_re_xe])/3.; // old approximation from Chen and Kamionkowski
chi_heat = MIN(0.996857*(1.-pow(1.-pow(x,0.300134),1.51035)),1.); // coefficient as revised by Galli et al. 2013 (in fact it is a fit by Vivian Poulin of columns 1 and 2 in Table V of Galli et al. 2013)
// coefficient as revised by Slatyer et al. 2013 (in fact it is a fit by Vivian Poulin of columns 1 and 2 in Table V of Slatyer et al. 2013)
if (x < 1.)
chi_heat = MIN(0.996857*(1.-pow(1.-pow(x,0.300134),1.51035)),1);
else
chi_heat = 1.;
dy[2]= preco->CT * pow(Trad,4) * x / (1.+x+preco->fHe) * (Tmat-Trad) / (Hz*(1.+z)) + 2.*Tmat/(1.+z)
-2./(3.*_k_B_)*energy_rate*chi_heat/n/(1.+preco->fHe+x)/(Hz*(1.+z)); /* energy injection */
......
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