Commit 6e3a3f2a authored by Julien Lesgourgues's avatar Julien Lesgourgues
Browse files

Added massless and massive neutrinos to tensor calculation

* devel:
  Test class restored
  Update the test to test tensor modes
  replaced input parameter N_eff by N_ur (N_eff deprecated but still works)
  three options for tensor method now working
  almost done with neutrino tensors
  still in this process, tm_photons_only and tm_exact work, not tm_massless_approximation
  in precess of implementing massless neutrino approx in tensors, still buggy
  removed warning for neutrinos in tensors
  ncdm in tensors works
  ndm still wrong
  in course of impelmenting massive neutrinos in tensors; compiles, but not tested and not finsihed
  Added massless neutrino contribution to tensor modes
parents 0bf59c07 2d0336ec
......@@ -27,9 +27,9 @@ T_cmb = 2.726
Omega_b = 0.05
#omega_b = 0.0266691
4) ultra-relativistic species / massless neutrino density: either 'N_eff' or 'Omega_ur' or 'omega_ur' (default: 'N_eff' set to 3.046)
4) ultra-relativistic species / massless neutrino density: either 'N_ur' or 'Omega_ur' or 'omega_ur' (default: 'N_ur' set to 3.046) (note: instead of 'N_ur' you can pass equivalently 'N_eff', although this syntax is deprecated)
N_eff = 3.046
N_ur = 3.046
#Omega_ur = 3.4869665e-05
#omega_ur = 1.70861334e-5
......@@ -209,12 +209,16 @@ modes = s
lensing = yes
6) list of initial conditions for scalars ('ad' for adiabatic, 'bi' for baryon isocurvature, 'cdi' for CDM isocurvature, 'nid' for neutrino density isocurvature, 'niv' for neutrino velocity isocurvature). More than one of these allowed, can be attached or separated by arbitrary characters; letters can be small or capital. (default: set to 'ad')
6) which perturbations should be included in tensor calculations? write 'exact' to include photons, ultra-relativistic species 'ur' and all non-cold dark matter species 'ncdm'; write 'massless' to appriximate 'ncdm' as extra relativistic species (good approximation if ncdm is still relativistic at the time of recombination); write 'photons' to include only photons (default: 'massless')
tensor method =
7) list of initial conditions for scalars ('ad' for adiabatic, 'bi' for baryon isocurvature, 'cdi' for CDM isocurvature, 'nid' for neutrino density isocurvature, 'niv' for neutrino velocity isocurvature). More than one of these allowed, can be attached or separated by arbitrary characters; letters can be small or capital. (default: set to 'ad')
ic = ad
#ic = ad&bi&nid
7) gauge in which calculations are performed: 'sync' or 'synchronous' or 'Synchronous' for synchronous, 'new' or 'newtonian' or 'Newtonian' for Newtonian/longitudinal gauge (default: set to synchronous)
8) gauge in which calculations are performed: 'sync' or 'synchronous' or 'Synchronous' for synchronous, 'new' or 'newtonian' or 'Newtonian' for Newtonian/longitudinal gauge (default: set to synchronous)
gauge = synchronous
......
......@@ -43,6 +43,7 @@ enum tca_method {first_order_MB,first_order_CAMB,first_order_CLASS,second_order_
enum rsa_method {rsa_null,rsa_MD,rsa_MD_with_reio,rsa_none};
enum ufa_method {ufa_mb,ufa_hu,ufa_CLASS,ufa_none};
enum ncdmfa_method {ncdmfa_mb,ncdmfa_hu,ncdmfa_CLASS,ncdmfa_none};
enum tensor_methods {tm_photons_only,tm_massless_approximation,tm_exact};
//@}
......@@ -114,9 +115,14 @@ struct perturbs
short has_nid; /**< do we need isocurvature nid mode? */
short has_niv; /**< do we need isocurvature niv mode? */
/* perturbed recombination */
/* Do we want to consider perturbed temperature and ionization fraction? */
short has_perturbed_recombination;
/* perturbed recombination */
/* Do we want to consider perturbed temperature and ionization fraction? */
short has_perturbed_recombination;
/** Neutrino contribution to tensors */
enum tensor_methods tensor_method; /**< way to treat neutrinos in tensor perturbations(neglect, approximate as massless, take exact equations) */
short evolve_tensor_ur; /**< will we evolve ur tensor perturbations (either becasue we have ur species, or we have ncdm species with massless approximation) ? */
short evolve_tensor_ncdm; /**< will we evolve ncdm tensor perturbations (if we have ncdm species and we use the exact method) ? */
short has_cl_cmb_temperature; /**< do we need Cl's for CMB temperature? */
short has_cl_cmb_polarization; /**< do we need Cl's for CMB polarization? */
......@@ -411,6 +417,7 @@ struct perturb_workspace
int index_mt_eta_prime; /**< eta' (wrt conf. time) in synchronous gauge */
int index_mt_alpha; /**< \alpha = (h' + 6 \eta') / (2 k^2) \f$ in synchronous gauge */
int index_mt_alpha_prime; /**< alpha' wrt conf. time) in synchronous gauge */
int index_mt_gw_prime_prime;/**< second derivative wrt confromal time of gravitational wave field, often called h */
int mt_size; /**< size of metric perturbation vector */
//@}
......@@ -432,6 +439,7 @@ struct perturb_workspace
double rho_plus_p_theta;
double rho_plus_p_shear;
double delta_p;
double gw_source;
double tca_shear_g; /**< photon shear in tight-coupling approximation */
double tca_slip; /**< photon-baryon slip in tight-coupling approximation */
......@@ -676,6 +684,7 @@ extern "C" {
struct background * pba,
struct thermo * pth,
struct perturbs * ppt,
int index_md,
double k,
double * y,
struct perturb_workspace * ppw
......
......@@ -40,7 +40,8 @@ class TestClass(unittest.TestCase):
self.scenario = {'lensing':'yes'}
def tearDown(self):
self.cosmo.cleanup()
self.cosmo.struct_cleanup()
self.cosmo.empty()
del self.scenario
@parameterized.expand(
......@@ -132,5 +133,34 @@ class TestClass(unittest.TestCase):
print '~~~~~~~~ passed ? '
@parameterized.expand(
itertools.product(
('massless', 'massive', 'both'),
('photons', 'massless', 'exact'),
('t', 's, t')))
def test_tensors(self, scenario, method, modes):
"""Test the new tensor mode implementation"""
self.scenario = {}
if scenario == 'massless':
self.scenario.update({'N_eff': 3.046, 'N_ncdm':0})
elif scenario == 'massiv':
self.scenario.update(
{'N_eff': 0, 'N_ncdm': 2, 'm_ncdm': '0.03, 0.04',
'deg_ncdm': '2, 1'})
elif scenario == 'both':
self.scenario.update(
{'N_eff': 1.5, 'N_ncdm': 2, 'm_ncdm': '0.03, 0.04',
'deg_ncdm': '1, 0.5'})
self.scenario.update({
'tensor method': method, 'modes': modes,
'output': 'tCl, pCl'})
for key, value in self.scenario.iteritems():
sys.stderr.write("%s = %s\n" % (key, value))
sys.stderr.write("\n")
self.cosmo.set(
dict(self.verbose.items()+self.scenario.items()))
self.cosmo.compute()
if __name__ == '__main__':
unittest.main()
......@@ -295,15 +295,37 @@ int input_init(
Omega_tot += pba->Omega0_b;
/* Omega_0_ur (ultra-relativistic species / massless neutrino) */
class_call(parser_read_double(pfc,"N_eff",&param1,&flag1,errmsg),
/* (a) try to read N_ur */
class_call(parser_read_double(pfc,"N_ur",&param1,&flag1,errmsg),
errmsg,
errmsg);
/* these lines have been added for coimpatibility with deprecated syntax 'N_eff' instead of 'N_ur', in the future they could be supressed */
class_call(parser_read_double(pfc,"N_eff",&param2,&flag2,errmsg),
errmsg,
errmsg);
class_test((flag1 == _TRUE_) && (flag2 == _TRUE_),
errmsg,
"In input file, you can only enter one of N_eff (deprecated syntax) or N_ur (up-to-date syntax), since they botgh describe the same, i.e. the contribution ukltra-relativistic species to the effective neutrino number");
if (flag2 == _TRUE_) {
param1 = param2;
flag1 = _TRUE_;
flag2 = _FALSE_;
}
/* end of lines for deprecated syntax */
/* (b) try to read Omega_ur */
class_call(parser_read_double(pfc,"Omega_ur",&param2,&flag2,errmsg),
errmsg,
errmsg);
/* (c) try to read omega_ur */
class_call(parser_read_double(pfc,"omega_ur",&param3,&flag3,errmsg),
errmsg,
errmsg);
/* (d) infer the unpassed ones from the passed one */
class_test(class_at_least_two_of_three(flag1,flag2,flag3),
errmsg,
"In input file, you can only enter one of N_eff, Omega_ur or omega_ur, choose one");
......@@ -1649,6 +1671,21 @@ int input_init(
class_read_int("tight_coupling_approximation",ppr->tight_coupling_approximation);
if (ppt->has_tensors == _TRUE_) {
/** Include ur and ncdm shear in tensor computation? */
class_call(parser_read_string(pfc,"tensor method",&string1,&flag1,errmsg),
errmsg,
errmsg);
if (flag1 == _TRUE_) {
if (strstr(string1,"photons") != NULL)
ppt->tensor_method = tm_photons_only;
if (strstr(string1,"massless") != NULL)
ppt->tensor_method = tm_massless_approximation;
if (strstr(string1,"exact") != NULL)
ppt->tensor_method = tm_exact;
}
}
/** derivatives of baryon sound speed only computed if some non-minimal tight-coupling schemes is requested */
if ((ppr->tight_coupling_approximation == (int)first_order_CLASS) || (ppr->tight_coupling_approximation == (int)second_order_CLASS)) {
pth->compute_cb2_derivatives = _TRUE_;
......@@ -2008,6 +2045,9 @@ int input_default_params(
ppt->has_niv=_FALSE_;
ppt->has_perturbed_recombination=_FALSE_;
ppt->tensor_method = tm_massless_approximation;
ppt->evolve_tensor_ur = _FALSE_;
ppt->evolve_tensor_ncdm = _FALSE_;
ppt->has_scalars=_TRUE_;
ppt->has_vectors=_FALSE_;
......
This diff is collapsed.
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