Commit 80c0d527 authored by Julia Stadler's avatar Julia Stadler
Browse files

first implementation, gives bugs

parent e65e4565
......@@ -525,6 +525,11 @@ struct precision
double start_large_k_at_tau_h_over_tau_k; /**< largest wavelengths start being sampled when mode is sufficiently outside Hubble scale. This is quantified in terms of the ratio of hubble time scale to wavenumber time scale, \f$ \tau_h/\tau_k \f$ which is roughly equal to (k*tau). Start when this ratio equals start_large_k_at_tau_k_over_tau_h. Decrease this value to start integrating the wavenumbers earlier in time. */
/** if we have neutrino interactions we need three new precision parameters to controll how we set initial conditions **/
int has_nuDM_initially;
double start_small_k_at_dmu_nuDM_over_aH;
double start_large_k_at_aH_over_dmu_nuDM;
/**
* when to switch off tight-coupling approximation: first condition:
* \f$ \tau_c/\tau_H \f$ > tight_coupling_trigger_tau_c_over_tau_h.
......
......@@ -712,6 +712,7 @@ extern "C" {
int perturb_initial_conditions(
struct precision * ppr,
struct background * pba,
struct thermo * pth,
struct perturbs * ppt,
int index_md,
int index_ic,
......
......@@ -20,7 +20,7 @@ N_ur = 3.046
---> g-cdm coupling
-------------------
--> u_nuDM(T) is given by u_nuDM_0 * a^n_nuDM
u_nuDM_0 = 0.
u_nuDM_0 = 1.e-4
n_nuDM = 0.
--> for alpha_nuDM give a list of values which apply to increasing moltipoles
......@@ -34,16 +34,23 @@ gauge = newtonian
k_pivot = 0.05
sBBN file = ./bbn/sBBN.dat
--------------------------------------------
---> interacting neutrino initial conditions
--------------------------------------------
has_nuDM_initially = 0
start_small_k_at_dmu_nuDM_over_aH = 0.001
start_large_k_at_aH_over_dmu_nuDM = 0.001
-----------
---> OUTPUT
-----------
output = tCl, pCl, mPk, dTk
extra metric transfer functions = yes
k_output_values = 0.3366, 2.0196, 3.366, 13.464, 20.196, 26.926, 53.806
k_output_values = 0.3366, 2.0196, 3.366, 13.464, 20.196, 26.926
k_per_decade_for_pk = 100
P_k_max_h/Mpc = 100
P_k_max_h/Mpc = 50
write parameters = yes
root = output/2019-05-21_jcap-review/lcdm_
root = output/2019-05-24_initial-conditions/u1em4_old_
---------------------------
---> some precision settings
......@@ -51,12 +58,12 @@ root = output/2019-05-21_jcap-review/lcdm_
ur_fluid_trigger_tau_over_tau_k = 200
l_max_ur = 200
tol_perturb_integration = 1.e-9 # default: 1.e-5
tol_background_integration = 1.e-4 # default: 1.e-2
tol_tau_approx=1.e-11 # default: 1.e-10
perturb_sampling_stepsize=0.01 # default: 0.1
tol_thermo_integration=1.e-5 # default: 1.e-2
recfast_Nz0=100000 # default: 20,000
#tol_perturb_integration = 1.e-9 # default: 1.e-5
#tol_background_integration = 1.e-4 # default: 1.e-2
#tol_tau_approx=1.e-11 # default: 1.e-10
#perturb_sampling_stepsize=0.01 # default: 0.1
#tol_thermo_integration=1.e-5 # default: 1.e-2
#recfast_Nz0=100000 # default: 20,000
-----------------------
---> breake the silence
......
......@@ -500,7 +500,12 @@ int input_init(
printf(" alpha_nuDM = ");
for (i=0; i<ppr->l_max_ur; i++)
printf("%f, ", ppt->alpha_nuDM[i]);
printf("\n");
printf("\n");
printf("has_nuDM_initially = %d\n", ppr->has_nuDM_initially);
if (ppr->has_nuDM_initially == _TRUE_){
printf("small k threshold = %f\n", ppr->start_small_k_at_dmu_nuDM_over_aH);
printf("large k threshold = %f\n", ppr->start_large_k_at_aH_over_dmu_nuDM);
}
}
return _SUCCESS_;
......@@ -1308,7 +1313,12 @@ int input_read_parameters(
class_read_double("u_nuDM_0", pth->u_nuDM_0);
if (pth->u_nuDM_0>0.){
pth->has_coupling_nuDM = _TRUE_;
class_read_double("n_nuDM", pth->n_nuDM)
class_read_double("n_nuDM", pth->n_nuDM);
class_read_int("has_nuDM_initially", ppr->has_nuDM_initially);
if(ppr->has_nuDM_initially){
class_read_double("start_small_k_at_dmu_nuDM_over_aH", ppr->start_small_k_at_dmu_nuDM_over_aH);
class_read_double("start_large_k_at_aH_over_dmu_nuDM", ppr->start_large_k_at_aH_over_dmu_nuDM);
}
}
......@@ -3394,6 +3404,11 @@ int input_default_precision ( struct precision * ppr ) {
ppr->start_sources_at_tau_c_over_tau_h = 0.008; /* decrease to start earlier in time */
ppr->tight_coupling_approximation=(int)compromise_CLASS;
/** parameters for nuDM initial conditions **/
ppr->has_nuDM_initially=_TRUE_;
ppr->start_small_k_at_dmu_nuDM_over_aH=1.;
ppr->start_large_k_at_aH_over_dmu_nuDM=1.;
ppr->l_max_g=12;
ppr->l_max_pol_g=10;
ppr->l_max_dr=17;
......
......@@ -2158,7 +2158,7 @@ int perturb_solve(
/* approximation scheme within previous interval: previous_approx[index_ap] */
int * previous_approx;
int n_ncdm,is_early_enough;
int n_ncdm,is_early_enough,account_for_nuDM_interactions;
/* function pointer to ODE evolver and names of possible evolvers */
......@@ -2249,6 +2249,25 @@ int perturb_solve(
}
}
/** check if we have to impose a condition on (aH/kappa_nuDM) **/
account_for_nuDM_interactions == _FALSE_;
if(pth->has_coupling_nuDM==_TRUE_ && ppr->has_nuDM_initially==_TRUE_){
if(ppw->pvecthermo[pth->index_th_dmu_nuDM]/
ppw->pvecback[pba->index_bg_a]/
ppw->pvecback[pba->index_bg_H]
> ppr->start_small_k_at_dmu_nuDM_over_aH)
account_for_nuDM_interactions = _TRUE_;
}
if(account_for_nuDM_interactions = _TRUE_){
class_test(ppw->pvecback[pba->index_bg_a]*
ppw->pvecback[pba->index_bg_H]/
ppw->pvecthermo[pth->index_th_dmu_nuDM] > ppr->start_large_k_at_aH_over_dmu_nuDM,
ppt->error_message, "your choice of initial time for integrating wavenumbers is inappropriate: it corresponds to a time before that at which the background has been integrated. You should increase 'start_large_k_at_aH_over_dmu_nuDM' up to at least %g, or decrease 'a_ini_over_a_today_default'\n",
ppw->pvecback[pba->index_bg_a]*
ppw->pvecback[pba->index_bg_H]/
ppw->pvecthermo[pth->index_th_dmu_nuDM]);
}
/* is at most the time at which sources must be sampled */
tau_upper = ppt->tau_sampling[0];
......@@ -2276,6 +2295,25 @@ int perturb_solve(
}
}
/* check if the neutrino interactions need to be taken into account */
if(pth->has_coupling_nuDM==_TRUE_ && ppr->has_nuDM_initially==_TRUE_){
if(ppw->pvecthermo[pth->index_th_dmu_nuDM]/
ppw->pvecback[pba->index_bg_H]/
ppw->pvecback[pba->index_bg_a]
< ppr->start_small_k_at_dmu_nuDM_over_aH)
account_for_nuDM_interactions = _FALSE_;
}
/* check if integration starts early enough for nuDM initial conditions */
if(account_for_nuDM_interactions == _TRUE_){
if(ppw->pvecback[pba->index_bg_H]*
ppw->pvecback[pba->index_bg_a]/
ppw->pvecthermo[pth->index_th_dmu_nuDM]
> ppr->start_large_k_at_aH_over_dmu_nuDM)
is_early_enough = _FALSE_;
}
/* also check that the two conditions on (aH/kappa') and (aH/k) are fulfilled */
if (is_early_enough == _TRUE_) {
......@@ -3450,6 +3488,7 @@ int perturb_vector_init(
class_call(perturb_initial_conditions(ppr,
pba,
pth,
ppt,
index_md,
index_ic,
......@@ -4066,6 +4105,7 @@ int perturb_vector_free(
int perturb_initial_conditions(struct precision * ppr,
struct background * pba,
struct thermo * pth,
struct perturbs * ppt,
int index_md,
int index_ic,
......@@ -4092,6 +4132,8 @@ int perturb_initial_conditions(struct precision * ppr,
double velocity_tot;
double s2_squared;
int account_for_nuDM_interactions;
/** --> For scalars */
if (_scalars_) {
......@@ -4186,6 +4228,27 @@ int perturb_initial_conditions(struct precision * ppr,
s2_squared = 1.-3.*pba->K/k/k;
//////////////////////////////////////////////////
/* If neutrinos interact with dark matter this supresses their shear and changes the initial
conditions.Here we check if we have to take this into account */
account_for_nuDM_interactions = _FALSE_;
if (pth->has_coupling_nuDM && ppr->has_nuDM_initially){
class_call(thermodynamics_at_z(pba,
pth,
1./ppw->pvecback[pba->index_bg_a]-1., /* redshift z=1/a-1 */
pth->inter_normal,
&(ppw->last_index_thermo),
ppw->pvecback,
ppw->pvecthermo),
pth->error_message,
ppt->error_message);
if(ppw->pvecthermo[pth->index_th_dmu_nuDM]/a_prime_over_a < ppr->start_small_k_at_dmu_nuDM_over_aH)
account_for_nuDM_interactions = _TRUE_;
}
//////////////////////////////////////////////////
/** - (b) starts by setting everything in synchronous gauge. If
another gauge is needed, we will perform a gauge
transformation below. */
......@@ -4278,13 +4341,30 @@ int perturb_initial_conditions(struct precision * ppr,
if (pba->has_dr == _TRUE_) delta_dr = delta_ur;
//////////////////////////////////////////////////
if (account_for_nuDM_interactions == _TRUE_)
{
/* same as photon velocity */
theta_ur = ppw->pv->y[ppw->pv->index_pt_theta_g];
shear_ur = 0.;
/* tighly-coupled dark matter */
ppw->pv->y[ppw->pv->index_pt_theta_cdm] = theta_ur;
}
//////////////////////////////////////////////////
}
/* synchronous metric perturbation eta */
//eta = ppr->curvature_ini * (1.-ktau_two/12./(15.+4.*fracnu)*(5.+4.*fracnu - (16.*fracnu*fracnu+280.*fracnu+325)/10./(2.*fracnu+15.)*tau*om)) / s2_squared;
//eta = ppr->curvature_ini * s2_squared * (1.-ktau_two/12./(15.+4.*fracnu)*(15.*s2_squared-10.+4.*s2_squared*fracnu - (16.*fracnu*fracnu+280.*fracnu+325)/10./(2.*fracnu+15.)*tau*om));
eta = ppr->curvature_ini * (1.-ktau_two/12./(15.+4.*fracnu)*(5.+4.*s2_squared*fracnu - (16.*fracnu*fracnu+280.*fracnu+325)/10./(2.*fracnu+15.)*tau*om));
//////////////////////////////////////////////////
if (account_for_nuDM_interactions == _TRUE_){
/* ToDo: this is only the lowest order expression as in Ma/Bertschinger
Obtain higher order corrections. */
eta = 1 - 0.5/18.*ktau_two;
}
//////////////////////////////////////////////////
}
/* isocurvature initial conditions taken from Bucher, Moodely,
......@@ -4424,7 +4504,7 @@ int perturb_initial_conditions(struct precision * ppr,
if (ppt->gauge == newtonian) {
/* alpha is like in Ma & Bertschinger: (h'+6 eta')/(2k^2). We obtain it from the first two Einstein equations:
alpha = [eta + 3/2 (a'/a)^2 (delta_rho/rho_c) / k^2 /s_2^2 + 3/2 (a'/a)^3 3 ((rho+p)theta/rho_c) / k^4 / s_2^2] / (a'/a)
= [eta + 3/2 (a'/a)^2 / k^2 /s_2^2 {delta_tot + 3 (a'/a) /k^2 velocity_tot}] / (a'/a)
......
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