From a33a5870cbdc599b71ce3eac617bc948ef3de1f8 Mon Sep 17 00:00:00 2001 From: Emilio Bellini Date: Tue, 3 May 2016 12:58:30 +0200 Subject: [PATCH 1/2] Modified a bit stability test for D. Moved stability tests before the isnan test. Modified hi_class.ini --- hi_class.ini | 79 +++++++++++++++++++++++---------------------- source/background.c | 51 +++++++++++++++-------------- 2 files changed, 67 insertions(+), 63 deletions(-) diff --git a/hi_class.ini b/hi_class.ini index 57232541..d74892e4 100644 --- a/hi_class.ini +++ b/hi_class.ini @@ -6,6 +6,8 @@ > with detailed comments. The idea here is that the parameters used also in CLASS > are usually not reported, i.e. they assume their standard value. For a complete > list of these parameters look at the "explanatory.ini" file. +> All the quantities related to the scalar field introduced in hi_class have been +> marked with the "_smg" label (scalar modified gravity). > Only lines containing an equal sign not preceded by a sharp sign "#" are > considered by the code. Hence, do not write an equal sign within a comment, > the whole line would be interpreted as relevant input. Input files must have @@ -16,25 +18,23 @@ ----> background parameters: ---------------------------- -1a) Omega_smg is the fractional density of the scalar field today. It can be specified or - found automatically if Omega_Lambda and Omega_fld are given and Omega_smg set to -1. - This is to avoid someone modifying gravity by mistake. The code will then use Omega_smg - to satisfy the closure equation (sum_i Omega_i) equals (1 + Omega_k). If Omega_smg has - a positive (and smaller than 1) value and either Omega_Lambda or Omega_fld are left - unspecified, the code will use the closure equation to determine the unspecified - fractional density. +1a) Omega_smg is the fractional density of the scalar field today. There are three possibilities + i) Omega_smg unspecified or equivalently set to 0. In this case the code will + ignore the scalar field equations and will use the standard CLASS equations + ii) Omega_smg has a value larger than 0 but smaller than 1. In this case you should + leave either Omega_Lambda or Omega_fld unspecified. Then, hi_class will run + with the scalar field equations, and Omega_Lambda or Omega_fld will be inferred + using the closure equation (sum_i Omega_i) equals (1 + Omega_k) + iii) Omega_smg has a negqative value. In this case the equations for the scalar field + will be used, you have to specify both Omega_Lambda and Omega_fld, and Omega_smg will + be inferred by the code using the closure equation Omega_Lambda = 0 # fractional density of the cosmological constant today Omega_fld = 0 # fractional density of a perfect fluid today, see the "explanatory.ini" file for details on this fluid Omega_smg = -1 # fractional density of the DE field today. If it is specified to a given value, then need to comment Omega_Lambda or Omega_fld -1b) For debugging purposes one can specify Omega_smg_debug AND unspecify Omega_smg,fld,Lambda. - In this case the closure equation is not tested - -# Omega_smg_debug = 0. - -1c) For now, there are two possible expansion histories of the universe +1b) For now, there are two possible expansion histories of the universe i) "lcdm": fixes the expansion history to be the one predicted by a cosmological constant ii) "wowa": introduces a fluid with equation of state p/rho equal to w0+wa(1-a/a0) This information is stored in the "expansion_model" variable. After specifying it, you have @@ -45,8 +45,10 @@ Omega_smg = -1 # fractional density of the DE field today. If it is specified to expansion_model = lcdm expansion_smg = 0.5 #this value will be overwritten using the closure equation. See below +# expansion_model = wowa +# expansion_smg = 0.5, -1., 0. -1d) Given a vector containing the values of the parameters needed to fully fix the background +1c) Given a vector containing the values of the parameters needed to fully fix the background evolution you chose, you should specify which parameter you want to vary to satisfy the closure equation (tuning_index_smg, default: 0) and how much you want to vary it (tuning_dxdy_guess_smg, default: 1) @@ -55,39 +57,17 @@ expansion_smg = 0.5 #this value will be overwritten using the closure equation. # tuning_dxdy_guess_smg = 1. ---------------------------- -----> stability parameters: ---------------------------- - -2a) The absence of ghost and gradient instabilities is tested in the code for both scalars - and tensors. It is possible to avoid these tests using - -skip_stability_tests_smg = no # default: no - -2b) If skip_stability_tests_smg is set to no, it is possible to relax the stability conditions - (moslty to avoid numerical noise around 0). Indeed, the code internally verify that - i) cs2_smg > - abs(cs2_safe_smg). Sound speed of the scalar sector - ii) ct2_smg > - abs(ct2_safe_smg). Sound speed of the tensor sector - iii) D_smg > - abs(D_safe_smg). Kinetic term of the scalar sector - iv) M2_smg > - abs(M2_safe_smg). Kinetic term of the tensor sector - -cs2_safe_smg = 0. # threshold to consider the sound speed of scalars negative in the stability check -D_safe_smg = 0. # threshold to consider the kinetic term of scalars negative in the stability check -ct2_safe_smg = 0. # threshold to consider the sound speed of tensors negative in the stability check -M2_safe_smg = 0. # threshold to consider the kinetic term of tensors (M2) negative in the stability check - - ------------------------------ ----> perturbation parameters: ------------------------------ -3a) There are two possible initial conditions for the perturbations of the scalar field +2a) There are two possible initial conditions for the perturbations of the scalar field i) "single_clock": single_clock IC given with respect to photons (default) ii) "zero": set the initial value of the scalar field and its first derivative to 0 pert_initial_conditions_smg = single_clock -3b) At the perturbation level, there are four possible parametrizations of the alpha functions +2b) At the perturbation level, there are four possible parametrizations of the alpha functions i) "propto_omega": all the alphas are proportional to the fractional density of the Dark Energy fluid at the background. Then, you have to provide a vector containg the factor of proportionality of each alpha and the initial value of the Planck Mass @@ -110,6 +90,29 @@ pert_initial_conditions_smg = single_clock gravity_model = propto_omega parameters_smg = 1., 0., 0., 0., 1. + +--------------------------- +----> stability parameters: +--------------------------- + +3a) The absence of ghost and gradient instabilities is tested in the code for both scalars + and tensors. It is possible to avoid these tests using + +skip_stability_tests_smg = no # default: no + +3b) If skip_stability_tests_smg is set to no, it is possible to relax the stability conditions + (moslty to avoid numerical noise around 0). Indeed, the code internally verify that + i) cs2_smg > - abs(cs2_safe_smg). Sound speed of the scalar sector + ii) ct2_smg > - abs(ct2_safe_smg). Sound speed of the tensor sector + iii) D_smg > - abs(D_safe_smg). Kinetic term of the scalar sector + iv) M2_smg > - abs(M2_safe_smg). Kinetic term of the tensor sector + +cs2_safe_smg = 0. # threshold to consider the sound speed of scalars negative in the stability check +D_safe_smg = 0. # threshold to consider the kinetic term of scalars negative in the stability check +ct2_safe_smg = 0. # threshold to consider the sound speed of tensors negative in the stability check +M2_safe_smg = 0. # threshold to consider the kinetic term of tensors (M2) negative in the stability check + + --------------------------- ----> precision parameters: --------------------------- diff --git a/source/background.c b/source/background.c index d7ca7f6a..ff2e93b7 100755 --- a/source/background.c +++ b/source/background.c @@ -1964,7 +1964,32 @@ int background_solve( }//end of has_smg } - /* Yet another (third!) loop to make sure the background table makes sense + + /* Horndeski stability tests + * only if not overriden + * and model is tuned! + * TODO: commented out till properly checked! + */ + if ((pba->has_smg == _TRUE_) && + (pba->parameters_tuned_smg == _TRUE_) && + (pba->skip_stability_tests_smg == _FALSE_)){ + + class_test(pba->min_D_smg <= -fabs(pba->D_safe_smg), + pba->error_message, + "Ghost instability for scalar field perturbations with minimum D=%g \n",pba->min_D_smg); + class_test(pba->min_cs2_smg < -fabs(pba->cs2_safe_smg), + pba->error_message, + "Gradient instability for scalar field perturbations with minimum c_s^2=%g \n",pba->min_cs2_smg); + class_test(pba->min_M2_smg < -fabs(pba->M2_safe_smg), + pba->error_message, + "Ghost instability for metric tensor perturbations with minimum M*^2=%g \n",pba->min_M2_smg); + class_test(pba->min_ct2_smg < -fabs(pba->ct2_safe_smg), + pba->error_message, + "Gradient instability for metric tensor perturbations with minimum c_t^2=%g \n",pba->min_ct2_smg); + + } + + /* Yet another (third!) loop to make sure the background table makes sense */ for (i=0; i < pba->bt_size; i++) { @@ -2034,30 +2059,6 @@ int background_solve( free(pvecback); free(pvecback_integration); - /* Horndeski stability tests - * only if not overriden - * and model is tuned! - * TODO: commented out till properly checked! - */ - if ((pba->has_smg == _TRUE_) && - (pba->parameters_tuned_smg == _TRUE_) && - (pba->skip_stability_tests_smg == _FALSE_)){ - - class_test(pba->min_D_smg < -fabs(pba->D_safe_smg), - pba->error_message, - "Ghost instability for scalar field perturbations with minimum D=%g \n",pba->min_D_smg); - class_test(pba->min_cs2_smg < -fabs(pba->cs2_safe_smg), - pba->error_message, - "Gradient instability for scalar field perturbations with minimum c_s^2=%g \n",pba->min_cs2_smg); - class_test(pba->min_M2_smg < -fabs(pba->M2_safe_smg), - pba->error_message, - "Ghost instability for metric tensor perturbations with minimum M*^2=%g \n",pba->min_M2_smg); - class_test(pba->min_ct2_smg < -fabs(pba->ct2_safe_smg), - pba->error_message, - "Gradient instability for metric tensor perturbations with minimum c_t^2=%g \n",pba->min_ct2_smg); - - } - return _SUCCESS_; } From 84d5f74db013aa6fa62f5aaf4ad5bacc270408e6 Mon Sep 17 00:00:00 2001 From: Emilio Bellini Date: Wed, 4 May 2016 13:10:24 +0200 Subject: [PATCH 2/2] Corrected bug in the tensor equation --- source/perturbations.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/source/perturbations.c b/source/perturbations.c index 56a42aff..11f67dc7 100755 --- a/source/perturbations.c +++ b/source/perturbations.c @@ -5347,10 +5347,11 @@ int perturb_einstein( /* modified version if gravity is non-standard. Note that no curvature is allowed in this case */ else{ + double M2 = ppw->pvecback[pba->index_bg_M2_smg]; double run = ppw->pvecback[pba->index_bg_mpl_running_smg]; double c_t2 = (1. + ppw->pvecback[pba->index_bg_tensor_excess_smg]); - ppw->pvecmetric[ppw->index_mt_gw_prime_prime] = -(2. + run)*a_prime_over_a*y[ppw->pv->index_pt_gwdot]-k2*c_t2*y[ppw->pv->index_pt_gw]+ppw->gw_source; + ppw->pvecmetric[ppw->index_mt_gw_prime_prime] = -(2. + run)*a_prime_over_a*y[ppw->pv->index_pt_gwdot]-k2*c_t2*y[ppw->pv->index_pt_gw]+ppw->gw_source/M2; } }