From fa55c07f7465b7260886a16ef2d5fb6a99f7b375 Mon Sep 17 00:00:00 2001 From: Emilio Bellini Date: Tue, 23 Apr 2024 15:15:03 +0200 Subject: [PATCH] fixed todos --- gravity_smg/gravity_models_smg.c | 18 ++++++------------ gravity_smg/input_smg.c | 24 ++++++++++++------------ hi_class.ini | 18 ++++++++++++------ include/background.h | 1 + source/background.c | 9 ++++----- 5 files changed, 35 insertions(+), 35 deletions(-) diff --git a/gravity_smg/gravity_models_smg.c b/gravity_smg/gravity_models_smg.c index 690867f7..aaabb1a5 100644 --- a/gravity_smg/gravity_models_smg.c +++ b/gravity_smg/gravity_models_smg.c @@ -43,13 +43,8 @@ int gravity_models_gravity_properties_smg( * -> real models: "parameters_smg" to pba->parameters_smg * -> parameterizations: "parameters_smg" to pba->parameters_2_smg * "expansion_smg" to pba->parameters_smg - * NOTE: can change class_read_list_of_doubles_or_default <-> class_read_list_of_doubles - * to make it mandatory or allow for default values */ - // TODO_EB: if parameters_smg is not specified it gives segfault. There are two cases: - // - gravity_model is not specified either. In this case we redirect to propto_omega. We should have default values also for the parameters - // - gravity_model is specified. In this case it should return a documented error rather than segfault. if (strcmp(string1,"propto_omega") == 0) { pba->gravity_model_smg = propto_omega; pba->field_evolution_smg = _FALSE_; @@ -525,7 +520,7 @@ int gravity_models_expansion_properties_smg( flag2=_TRUE_; pba->parameters_size_smg = 1; pba->rho_evolution_smg=_FALSE_; - class_read_list_of_doubles_or_default("expansion_smg",pba->parameters_smg,0.0,pba->parameters_size_smg); + class_read_list_of_doubles("expansion_smg",pba->parameters_smg,pba->parameters_size_smg); } //accept different names @@ -534,7 +529,7 @@ int gravity_models_expansion_properties_smg( flag2=_TRUE_; pba->parameters_size_smg = 3; pba->rho_evolution_smg=_FALSE_; - class_read_list_of_doubles_or_default("expansion_smg",pba->parameters_smg,0.0,pba->parameters_size_smg); + class_read_list_of_doubles("expansion_smg",pba->parameters_smg,pba->parameters_size_smg); } if (strcmp(string1,"wowa_w") == 0 || strcmp(string1,"w0wa_w") == 0 || strcmp(string1,"cpl_w") == 0 ) { @@ -542,7 +537,7 @@ int gravity_models_expansion_properties_smg( flag2=_TRUE_; pba->parameters_size_smg = 3; pba->rho_evolution_smg=_TRUE_; - class_read_list_of_doubles_or_default("expansion_smg",pba->parameters_smg,0.0,pba->parameters_size_smg); + class_read_list_of_doubles("expansion_smg",pba->parameters_smg,pba->parameters_size_smg); } if (strcmp(string1,"wede") == 0) { @@ -550,7 +545,8 @@ int gravity_models_expansion_properties_smg( pba->expansion_model_smg = wede; flag2=_TRUE_; pba->parameters_size_smg = 3; - class_read_list_of_doubles_or_default("expansion_smg",pba->parameters_smg,0.0,pba->parameters_size_smg); + class_read_list_of_doubles("expansion_smg",pba->parameters_smg,pba->parameters_size_smg); + class_read_double("wede_Omega_e_regularizer_smg",pba->wede_Omega_e_regularizer_smg); // //optimize the guessing BUG: eventually leads to problem in the MCMC, perhaps the guess is too good? // if(pba->tuning_index_smg == 0){ // pba->parameters_smg[0] = pba->Omega0_smg; @@ -787,9 +783,7 @@ int gravity_models_get_back_par_smg( double Om0 = pba->parameters_smg[0]; double w0 = pba->parameters_smg[1]; - // TODO_EB: add a precision parameter for 1e-10 - //NOTE: I've regularized the expression adding a tiny Omega_e - double Ome = pba->parameters_smg[2] + 1e-10; + double Ome = pba->parameters_smg[2] + pba->wede_Omega_e_regularizer_smg; double Om = ((Om0 - Ome*(1.-pow(a,-3.*w0)))/(Om0 + (1.-Om0)*pow(a,3*w0)) + Ome*(1.-pow(a,-3*w0))); double dOm_da = (3*pow(a,-1 - 3*w0)*(-1 + Om0)*(-2*pow(a,3*w0)*(-1 + Om0)*Ome + Om0*Ome + pow(a,6*w0)*(Om0 - 2*Ome + Om0*Ome))*w0)/pow(-(pow(a,3*w0)*(-1 + Om0)) + Om0,2); //from Mathematica diff --git a/gravity_smg/input_smg.c b/gravity_smg/input_smg.c index 0ce7d194..439811f9 100644 --- a/gravity_smg/input_smg.c +++ b/gravity_smg/input_smg.c @@ -165,16 +165,14 @@ int input_read_parameters_smg( class_call(parser_read_string(pfc,"gravity_model",&string1,&flag1,errmsg), errmsg, errmsg); - - if (flag1 == _FALSE_) { - printf(" gravity_model not read, default will be used \n"); - } - else { - /** Fix flags for each gravity model */ - class_call(gravity_models_gravity_properties_smg(pfc, pba, string1, has_tuning_index_smg, has_dxdy_guess_smg, errmsg), - errmsg, - errmsg); - } + class_test(flag1 == _FALSE_, + errmsg, + "gravity_model not read, you should specify one!\n"); + + /** Fix flags for each gravity model */ + class_call(gravity_models_gravity_properties_smg(pfc, pba, string1, has_tuning_index_smg, has_dxdy_guess_smg, errmsg), + errmsg, + errmsg); // end of loop over models if(pba->field_evolution_smg == _TRUE_){ @@ -192,8 +190,9 @@ int input_read_parameters_smg( class_call(parser_read_string(pfc,"expansion_model",&string1,&flag1,errmsg), errmsg, errmsg); - if (flag1 == _FALSE_) - printf("No expansion model specified, will take default one \n"); + class_test(flag1 == _FALSE_, + errmsg, + "expansion_model not read, you should specify one!\n"); /** Fix flags for each expansion model */ class_call(gravity_models_expansion_properties_smg(pfc, pba, string1, errmsg), @@ -337,6 +336,7 @@ int input_default_params_smg( pba->hubble_evolution = _TRUE_; /** dynamical evolution of Friedmann eq. */ pba->hubble_friction = 3.; /** friction coefficient in H' equation: H' = ... + H_friction*(H^2 - rho_crit) [NOT ONLY IN SMG!] */ pba->rho_evolution_smg = _FALSE_; /*does the model require to evolve the background energy density? (only for parameterizations)*/ + pba->wede_Omega_e_regularizer_smg = 1.e-10; /** regularize adding a tiny Omega_e */ pba->kineticity_safe_smg = 0; /* value added to the kineticity, useful to cure perturbations at early time in some models */ pba->cs2_safe_smg = 0; /* threshold to consider the sound speed of scalars negative in the stability check */ diff --git a/hi_class.ini b/hi_class.ini index 007e63e1..3cb3e149 100644 --- a/hi_class.ini +++ b/hi_class.ini @@ -292,7 +292,7 @@ n_max_qs_smg = 1e4 start_small_k_at_tau_c_over_tau_h = 1e-4 start_large_k_at_tau_h_over_tau_k = 1e-4 -perturb_sampling_stepsize = 0.05 +perturbations_sampling_stepsize = 0.05 l_logstep = 1.045 l_linstep = 50 @@ -310,20 +310,26 @@ output = tCl,pCl,lCl,mPk 8b) file name root 'root' for all output files -root = output/test_ +root = output/test +overwrite_root = yes 8c) Do you want to write a file with the input parameters you used? Do you want to write a table of background quantitites in a file? -write parameters = yeap -write background = yeah +write_parameters = yeap +write_background = yeah k_output_values = 0.1, 0.01, 0.0001 8d) Verbose parameters input_verbose = 1 background_verbose = 1 -output_verbose = 1 thermodynamics_verbose = 1 perturbations_verbose = 1 -spectra_verbose = 1 +transfer_verbose = 1 +primordial_verbose = 1 +harmonic_verbose = 1 +fourier_verbose = 1 +lensing_verbose = 1 +distortions_verbose = 1 +output_verbose = 1 diff --git a/include/background.h b/include/background.h index 4ac78156..1280d4d3 100644 --- a/include/background.h +++ b/include/background.h @@ -149,6 +149,7 @@ struct background double hubble_friction; /** friction coefficient in H' equation: H' = ... + H_friction*(H^2 - rho_crit) [NOT ONLY IN SMG!] */ int hubble_evolution; /** whether to evolve H' from the equation */ + double wede_Omega_e_regularizer_smg; /** regularize adding a tiny Omega_e */ enum gravity_model gravity_model_smg; /** Horndeski model */ // enum gravity_model_subclass gravity_submodel_smg; /** Horndeski model */ diff --git a/source/background.c b/source/background.c index 9dcb0c6a..b759ffa6 100644 --- a/source/background.c +++ b/source/background.c @@ -624,8 +624,7 @@ int background_functions( /** - compute critical density */ rho_crit = rho_tot-pba->K/a/a; - // TODO_EB: here the condition in hi_class was (rho_crit <= 0.) && (pba->has_smg == _FALSE_). Why? - class_test(rho_crit <= 0., + class_test((rho_crit <= 0.) && (pba->initial_conditions_set_smg == _TRUE_), pba->error_message, "rho_crit = %e instead of strictly positive",rho_crit); @@ -1975,9 +1974,10 @@ int background_solve( class_alloc(pvecback_integration,pba->bi_size*sizeof(double),pba->error_message); /** - impose initial conditions with background_initial_conditions() */ - class_call(background_initial_conditions(ppr,pba,pvecback,pvecback_integration,&(loga_ini)), + class_call_except(background_initial_conditions(ppr,pba,pvecback,pvecback_integration,&(loga_ini)), pba->error_message, - pba->error_message); + pba->error_message, + free(pvecback);free(pvecback_integration);); /** - Determine output vector */ loga_final = 0.; // with our conventions, loga is in fact log(a/a_0); we integrate until today, when log(a/a_0) = 0 @@ -2403,7 +2403,6 @@ int background_initial_conditions( /* Just checking that our initial time indeed is deep enough in the radiation dominated regime (_smg) */ - // TODO_EB: There was a class_test_except with free(pvecback);free(pvecback_integration);background_free(pba) class_test(fabs(pvecback[pba->index_bg_Omega_r]+pvecback[pba->index_bg_Omega_de]-1.) > ppr->tol_initial_Omega_r, pba->error_message, "Omega_r = %e, Omega_de = %e, not close enough to 1. Decrease a_ini_over_a_today_default in order to start from radiation domination.",