From 03f9de975dcd62c8bffecf2095fec088c478f8fa Mon Sep 17 00:00:00 2001 From: Emilio Bellini Date: Fri, 19 Apr 2024 14:28:38 +0200 Subject: [PATCH] solved some todo --- gravity_smg/background_smg.c | 11 ++---- gravity_smg/fourier_smg.c | 2 +- gravity_smg/gravity_functions_smg.c | 12 +++--- gravity_smg/gravity_models_smg.c | 12 +++--- gravity_smg/include/hi_class.h | 4 +- gravity_smg/include/perturbations_smg.h | 2 - gravity_smg/input_smg.c | 5 +-- gravity_smg/perturbations_smg.c | 52 +++++++------------------ source/background.c | 2 +- source/input.c | 1 - source/perturbations.c | 8 ---- 11 files changed, 34 insertions(+), 77 deletions(-) diff --git a/gravity_smg/background_smg.c b/gravity_smg/background_smg.c index f360cbb0..38a830d2 100644 --- a/gravity_smg/background_smg.c +++ b/gravity_smg/background_smg.c @@ -1,6 +1,6 @@ /** @file background_smg.c Documented background_smg module * - * Emilio Bellini, Ignacy Sawicki, Miguel Zumalacarregui, TODO_EB: date here xx.xx.xxxx + * Emilio Bellini, Ignacy Sawicki, Miguel Zumalacarregui, 2024 * * Additional functions for the background module. * It contains all the hi_class related functions (_smg) @@ -30,7 +30,6 @@ * and phi is dimensionless (we can think of phi as given in units of the Planck mass * - A default module to numerically compute the derivatives when no analytic functions are given should be added. * Numerical derivatives may further serve as a consistency check. -* TODO: Add a background_write_alpha_primes * * @param pba Input: pointer to background structure * @param a Input: scale factor @@ -57,7 +56,6 @@ int background_gravity_functions_smg( pvecback[pba->index_bg_p_tot_wo_smg] = *ptr_p_tot; // scalar field + curvature not yet implemented - // TODO_EB: look for a better place to put this test class_test(pba->K !=0 , pba->error_message, "has_smg with curvature K = %e not yet implemented",pba->K); @@ -333,7 +331,7 @@ int background_gravity_functions_smg( *ptr_p_tot += pvecback[pba->index_bg_p_smg]; //divide relativistic & nonrelativistic (not very meaningful for oscillatory models) - //TODO: need to define menaingfully -> separate early universe (IC, BBN...) from late (Halofit...) + //TODO: need to define meaningfully -> separate early universe (IC, BBN...) from late (Halofit...) //BUG: causes problem with halofit!, if not, causes bug with Brans-Dicke *ptr_rho_de += pvecback[pba->index_bg_rho_smg]; @@ -731,7 +729,7 @@ int background_solve_smg( pba->error_message, pba->error_message); - // TODO_EB: note that the derivative is now calculated w.r.t. loga, while our _prime are w.r.t. tau + // NOTE: note that the derivative is now calculated w.r.t. loga, while our _prime are w.r.t. tau double factor = pvecback[pba->index_bg_a]*pvecback[pba->index_bg_H]; double d_over_dtau; @@ -967,7 +965,6 @@ int background_store_columntitles_smg( class_store_columntitle(titles,"slip_eff_smg",_TRUE_); } - //TODO: add in output background trigger if (pba->output_background_smg >= 3){ class_store_columntitle(titles,"kineticity_prime_smg",_TRUE_); class_store_columntitle(titles,"braiding_prime_smg",_TRUE_); @@ -1286,7 +1283,7 @@ int derivatives_alphas_smg( /* necessary for calling array_interpolate(), but never used */ int last_index; - // TODO_EB: note that the derivative is now calculated w.r.t. loga, while our _prime are w.r.t. tau + // NOTE: note that the derivative is now calculated w.r.t. log(a), while our _prime are w.r.t. tau double factor = pvecback[pba->index_bg_a]*pvecback[pba->index_bg_H]; double d_over_dtau; diff --git a/gravity_smg/fourier_smg.c b/gravity_smg/fourier_smg.c index 7bc34512..2e5f54e7 100644 --- a/gravity_smg/fourier_smg.c +++ b/gravity_smg/fourier_smg.c @@ -1,6 +1,6 @@ /** @file fourier_smg.c Documented fourier_smg module * - * Emilio Bellini, Ignacy Sawicki, Miguel Zumalacarregui, TODO_EB: date here xx.xx.xxxx + * Emilio Bellini, Ignacy Sawicki, Miguel Zumalacarregui, 2024 * * Additional functions for the fourier module. * It contains all the hi_class related functions (_smg) diff --git a/gravity_smg/gravity_functions_smg.c b/gravity_smg/gravity_functions_smg.c index 2ade65a1..fdf394da 100644 --- a/gravity_smg/gravity_functions_smg.c +++ b/gravity_smg/gravity_functions_smg.c @@ -1,6 +1,6 @@ /** @file gravity_functions_smg.c Documented gravity_functions_smg module * - * Emilio Bellini, Ignacy Sawicki, Miguel Zumalacarregui, TODO_EB: date here xx.xx.xxxx + * Emilio Bellini, Ignacy Sawicki, Miguel Zumalacarregui, 2024 * * This module contains all the complicated expressions * used in hi_class. In particular, they are casted in @@ -21,7 +21,6 @@ #include "gravity_functions_smg.h" -// TODO_EB: reread all the documentation and check in particular Input/Output /** * Get gravity functions Es from Gs. These are the functions * entering the Friedmann time-time constraint, which is @@ -56,7 +55,6 @@ int gravity_functions_Es_from_Gs_smg( double rho_tot = pvecback[pba->index_bg_rho_tot_wo_smg]; double p_tot = pvecback[pba->index_bg_p_tot_wo_smg]; - // TODO_EB: decide if it is better to keep it like this or add the pointers to the equations double G2=pgf->G2; double G2_X=pgf->G2_X, G2_XX=pgf->G2_XX, G2_XXX=pgf->G2_XXX; double G2_phi=pgf->G2_phi, G2_Xphi=pgf->G2_Xphi, G2_XXphi=pgf->G2_XXphi; @@ -658,7 +656,7 @@ int gravity_functions_As_from_alphas_smg( )*pow(H,-2); - // TODO_EB: remove below this when IC are recalculated + // NOTE: this could be written in terms of the A's and B's variables. pvecback[pba->index_bg_lambda_1_smg] = (run + (-1.)*ten)*(-3.)*bra + (1. + ten)*kin; @@ -690,7 +688,8 @@ int gravity_functions_As_from_alphas_smg( + (1. + beh)*bra_p/a/H + (2. - bra)*beh_p/a/H; - // TODO_EB: check if there is a better alternative to regularizing these quantities + // NOTE: this is to regularize cs2 when both the numerator and denominator are + // below numerical precision if (pvecback[pba->index_bg_cs2num_smg] == pvecback[pba->index_bg_kinetic_D_smg]) { pvecback[pba->index_bg_cs2_smg] = 1.; } @@ -698,7 +697,8 @@ int gravity_functions_As_from_alphas_smg( pvecback[pba->index_bg_cs2_smg] = pvecback[pba->index_bg_cs2num_smg]/pvecback[pba->index_bg_kinetic_D_smg]; } - // TODO_EB: rewrite Geff and slip for beyond Horndeski (calculate them in the hi_class.nb Mathematica notebook) + // NOTE: Geff and slip are calculated for Horndeski. Consider extending them to + // beyond Horndeski double beta_1 = (run + (-1.)*ten)*2. + (1. + ten)*bra; double beta_2 = 2.*beta_1 + (2. + (-2.)*M2 + bra*M2)*(rho_tot + p_tot)*(-3.)*pow(H,-2)*pow(M2,-1) + ((-2.) + bra)*(rho_smg + p_smg)*(-3.)*pow(H,-2) + 2.*pow(H,-1)*bra_p*pow(a,-1); diff --git a/gravity_smg/gravity_models_smg.c b/gravity_smg/gravity_models_smg.c index 6973346b..690867f7 100644 --- a/gravity_smg/gravity_models_smg.c +++ b/gravity_smg/gravity_models_smg.c @@ -641,7 +641,6 @@ int gravity_models_get_Gs_smg( } else if(pba->gravity_model_smg == galileon){ - //TODO: change name with ccc_galileon double M3 = pow(pba->H0,2); @@ -784,14 +783,14 @@ int gravity_models_get_back_par_smg( //Doran-Robbers model astro-ph/0601544 //as implemented in Pettorino et al. 1301.5279 - //TODO: check these expressions, they probably assume the standard evolution/friedmann eqs, etc... - //TODO: rewrite the expressions integrating the equation of state + //NOTE: rewrite the expressions integrating the equation of state 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; - //NOTE: I've regularized the expression adding a tiny Omega_e 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 //I took a_eq = a*rho_r/rho_m, with rho_r = 3*p_tot_wo_smg @@ -1128,7 +1127,7 @@ int gravity_models_initial_conditions_smg( double n = pba->parameters_smg[1]; double Rshift0 = pba->parameters_smg[2]; - double H = sqrt(rho_rad); //TODO: for low n -> need to solve tracker + Constraint simultaneously. Increasing H (e.g. double H = 10*sqrt(rho_rad);) works for n~0.65. + double H = sqrt(rho_rad); //NOTE: for low n -> need to solve tracker + Constraint simultaneously. Increasing H (e.g. double H = 10*sqrt(rho_rad);) works for n~0.65. double H0 = pba->H0; double signg = copysign(1.,g); @@ -1187,8 +1186,7 @@ int gravity_models_initial_conditions_smg( /* expansion_smg when w is parametrized and rho obtained with integration */ if (pba->rho_evolution_smg == _TRUE_){ - //DT: moved here from the definition of wowa_w model, to avoid pvecback[pba->index_bg_rho_smg] mix up - // TODO_EB: rethink this + //NOTE: moved here from the definition of wowa_w model, to avoid pvecback[pba->index_bg_rho_smg] mix up if (pba->expansion_model_smg == wowa_w) { double Omega_const_smg = pba->parameters_smg[0]; double w0 = pba->parameters_smg[1]; diff --git a/gravity_smg/include/hi_class.h b/gravity_smg/include/hi_class.h index e9d568ce..04b34383 100644 --- a/gravity_smg/include/hi_class.h +++ b/gravity_smg/include/hi_class.h @@ -1,9 +1,7 @@ -// TODO_EB: write some overall description - #ifndef __HI_CLASS__ #define __HI_CLASS__ -#define _HI_CLASS_VERSION_ "TODO_EB:write version" +#define _HI_CLASS_VERSION_ "3.0" #include "gravity_functions_smg.h" #include "gravity_models_smg.h" diff --git a/gravity_smg/include/perturbations_smg.h b/gravity_smg/include/perturbations_smg.h index b1609494..cfe569a6 100644 --- a/gravity_smg/include/perturbations_smg.h +++ b/gravity_smg/include/perturbations_smg.h @@ -59,8 +59,6 @@ int perturbations_get_h_prime_ic_from_00_smg(struct background * pba, double eta, double delta_rho_tot); -int perturbations_get_x_x_prime_newtonian_smg(struct perturbations_workspace * ppw); - int perturbations_einstein_scalar_smg(struct precision * ppr, struct background * pba, struct thermodynamics * pth, diff --git a/gravity_smg/input_smg.c b/gravity_smg/input_smg.c index ae0af746..0ce7d194 100644 --- a/gravity_smg/input_smg.c +++ b/gravity_smg/input_smg.c @@ -1,6 +1,6 @@ /** @file input_smg.c Documented input_smg module * - * Emilio Bellini, Ignacy Sawicki, Miguel Zumalacarregui, TODO_EB: date here xx.xx.xxxx + * Emilio Bellini, Ignacy Sawicki, Miguel Zumalacarregui, 2024 * * Additional functions for the input module. * It contains all the hi_class related functions (_smg) @@ -35,9 +35,6 @@ int input_warnings_smg( int input_verbose ) { - // TODO_EB: rethink this test. The default should have no warnings, and - // we have to be sure that get_h_from_trace == _TRUE_ is working properly - /* Here we put a warning as we want to encourage hi_class users to get h_prime from the trace of the Einstein ij equation than from the Einstein 00 equation. This is because the Einstein 00 equation has a gauge dependent diff --git a/gravity_smg/perturbations_smg.c b/gravity_smg/perturbations_smg.c index 8136d991..ace4030a 100644 --- a/gravity_smg/perturbations_smg.c +++ b/gravity_smg/perturbations_smg.c @@ -1,6 +1,6 @@ /** @file perturbations_smg.c Documented perturbations_smg module * - * Emilio Bellini, Ignacy Sawicki, Miguel Zumalacarregui, TODO_EB: date here xx.xx.xxxx + * Emilio Bellini, Ignacy Sawicki, Miguel Zumalacarregui, 2024 * * Additional functions for the perturbation module. * It contains all the hi_class related functions (_smg) @@ -364,14 +364,14 @@ int perturbations_prepare_k_output_smg( ) { if (ppt->use_pert_var_deltaphi_smg==_TRUE_) { - class_store_columntitle(ppt->scalar_titles, "delta_phi_smg", _TRUE_); - class_store_columntitle(ppt->scalar_titles, "delta_phi_prime_smg", _TRUE_); - class_store_columntitle(ppt->scalar_titles, "delta_phi_prime_prime_smg", _TRUE_); + class_store_columntitle(ppt->scalar_titles, "delta_phi_smg (sync)", _TRUE_); + class_store_columntitle(ppt->scalar_titles, "delta_phi_prime_smg (sync)", _TRUE_); + class_store_columntitle(ppt->scalar_titles, "delta_phi_prime_prime_smg (sync)", _TRUE_); } else { - class_store_columntitle(ppt->scalar_titles, "V_x_smg", _TRUE_); - class_store_columntitle(ppt->scalar_titles, "V_x_prime_smg", _TRUE_); - class_store_columntitle(ppt->scalar_titles, "V_x_prime_prime_smg", _TRUE_); + class_store_columntitle(ppt->scalar_titles, "V_x_smg (sync)", _TRUE_); + class_store_columntitle(ppt->scalar_titles, "V_x_prime_smg (sync)", _TRUE_); + class_store_columntitle(ppt->scalar_titles, "V_x_prime_prime_smg (sync)", _TRUE_); } /* Quasi-static functions smg*/ @@ -486,35 +486,13 @@ int perturbations_get_h_prime_ic_from_00_smg( double p_smg = ppw->pvecback[pba->index_bg_p_smg]; double kin = ppw->pvecback[pba->index_bg_kineticity_smg]; double bra = ppw->pvecback[pba->index_bg_braiding_smg]; - /* TODO_EB: rewrite this equation with new variables */ + /* NOTE: this could be rewritten this equation with A and B variables */ ppw->pv->y[ppw->pv->index_pt_h_prime_from_trace] = (-4. * pow(H, -1) * pow(k, 2) * eta / a - 6. * pow(H, -1) * pow(M2, -1) * delta_rho_tot * a + 2. * H * (3. * bra + kin) * ppw->pv->y[ppw->pv->index_pt_x_prime_smg] * a + (2. * bra * pow(k, 2) + (-18. + 15. * bra + 2. * kin) * rho_smg * pow(a, 2) + (-18. * DelM2 + 15. * bra * M2 + 2. * kin * M2) * rho_tot * pow(M2, -1) * pow(a, 2) + (-2. * DelM2 + bra * M2) * 9. * pow(M2, -1) * p_tot * pow(a, 2) + 9. * (-2. + bra) * p_smg * pow(a, 2)) * ppw->pv->y[ppw->pv->index_pt_x_smg]) * pow(-2. + bra, -1); return _SUCCESS_; } -/** - * Transform synchronous gauge scalar field to newtonian. - * - * @param ppw Input/Output: pointer to perturbation workspace structure - * @return the error status - */ -int perturbations_get_x_x_prime_newtonian_smg( - struct perturbations_workspace * ppw - ) { - - int qs_array_smg[] = _VALUES_QS_SMG_FLAGS_; - - /* scalar field: TODO_EB: add gauge transformations (when we add Newtonian gauge) */ - if (qs_array_smg[ppw->approx[ppw->index_ap_qs_smg]] == 0) { - ppw->pv->y[ppw->pv->index_pt_x_smg] += 0.; - ppw->pv->y[ppw->pv->index_pt_x_prime_smg] += 0.; - } - - return _SUCCESS_; -} - - /** * Scalar Einstein equations. * @@ -2027,7 +2005,7 @@ int perturbations_adiabatic_ic_smg( a_prime_over_a = ppw->pvecback[pba->index_bg_H]*a; - H = ppw->pvecback[pba->index_bg_H];//TODO_EB + H = ppw->pvecback[pba->index_bg_H]; Hprime = ppw->pvecback[pba->index_bg_H_prime]; a = ppw->pvecback[pba->index_bg_a]; rho_tot = ppw->pvecback[pba->index_bg_rho_tot_wo_smg]; @@ -2068,7 +2046,7 @@ int perturbations_adiabatic_ic_smg( DelM2 = ppw->pvecback[pba->index_bg_delta_M2_smg];//M2-1 - /* TODO_EB: revisit initial conditions for beyond horndeski and oscillations */ + /* NOTE: initial conditions are computed for EFT Horndeski. Beyond horndeski and background scalar field oscillations could be added */ if (qs_array_smg[ppw->approx[ppw->index_ap_qs_smg]] == 0) { @@ -2412,7 +2390,7 @@ int perturbations_adiabatic_ic_smg( shear_ur = A_sigma_nu_smg/A1_eta_smg* ktau_two * ppr->curvature_ini; // /(45.+12.*fracnu) * (3.*s2_squared-1.) * (1.+(4.*fracnu-5.)/4./(2.*fracnu+15.)*tau*om) * ppr->curvature_ini;//TBC /s2_squared; /* shear of ultra-relativistic neutrinos/relics */ //TBC:0 - //TODO: needs to be modified? + //TODO: needs to be modified? (ask ILS) l3_ur = ktau_three/A1_eta_smg*2./7./(12.*fracnu+45.)* ppr->curvature_ini;//ILS if(ppt->perturbations_verbose > 8) @@ -2633,7 +2611,7 @@ int perturbations_isocurvature_cdm_ic_smg( //only set ICs for smg if have smg, we are in exernal field attractor and we are *not* quasi-static if((ppt->pert_initial_conditions_smg==ext_field_attr) && (qs_array_smg[ppw->approx[ppw->index_ap_qs_smg]] == 0)) { - /* TODO_EB: revisit isocurvature initial conditions for beyond horndeski and oscillations */ + /* NOTE: initial conditions are computed for EFT Horndeski. Beyond horndeski and background scalar field oscillations could be added */ double coeff_isocurv_smg; @@ -2718,7 +2696,7 @@ int perturbations_isocurvature_b_ic_smg( //only set ICs for smg if have smg, we are in exernal field attractor and we are *not* quasi-static if((ppt->pert_initial_conditions_smg==ext_field_attr)&&(qs_array_smg[ppw->approx[ppw->index_ap_qs_smg]] == 0)) { - /* TODO_EB: revisit isocurvature initial conditions for beyond horndeski and oscillations */ + /* NOTE: initial conditions are computed for EFT Horndeski. Beyond horndeski and background scalar field oscillations could be added */ double coeff_isocurv_smg; int nexpo= 1; @@ -2808,7 +2786,7 @@ int perturbations_isocurvature_urd_ic_smg( //only set ICs for smg if have smg, we are in exernal field attractor and we are *not* quasi-static if((ppt->pert_initial_conditions_smg==ext_field_attr)&&(qs_array_smg[ppw->approx[ppw->index_ap_qs_smg]] == 0)) { - /* TODO_EB: revisit isocurvature initial conditions for beyond horndeski and oscillations */ + /* NOTE: initial conditions are computed for EFT Horndeski. Beyond horndeski and background scalar field oscillations could be added */ double coeff_isocurv_smg; double a = ppw->pvecback[pba->index_bg_a]; @@ -2913,7 +2891,7 @@ int perturbations_isocurvature_urv_ic_smg( //only set ICs for smg if have smg, we are in exernal field attractor and we are *not* quasi-static if((pba->has_smg == _TRUE_)&&(ppt->pert_initial_conditions_smg==ext_field_attr)&&(qs_array_smg[ppw->approx[ppw->index_ap_qs_smg]] == 0)) { - /* TODO_EB: revisit isocurvature initial conditions for beyond horndeski and oscillations */ + /* NOTE: initial conditions are computed for EFT Horndeski. Beyond horndeski and background scalar field oscillations could be added */ double coeff_isocurv_smg; diff --git a/source/background.c b/source/background.c index 11a47ed7..9dcb0c6a 100644 --- a/source/background.c +++ b/source/background.c @@ -2403,7 +2403,7 @@ int background_initial_conditions( /* Just checking that our initial time indeed is deep enough in the radiation dominated regime (_smg) */ - // TODO_EB: rethink this test (r+(de?)). There was a class_test_except with free(pvecback);free(pvecback_integration);background_free(pba) + // 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.", diff --git a/source/input.c b/source/input.c index 4f532d05..11babe9d 100644 --- a/source/input.c +++ b/source/input.c @@ -555,7 +555,6 @@ int input_shooting(struct file_content * pfc, /* for each target, module up to which we need to run CLASS in order to compute the targetted quantities (not running the whole code each time to saves a lot of time) */ - // TODO_EB: here we didn't add any computation stage in the previous version of hi_class, I think we should. Check it enum computation_stage target_cs[] = {cs_thermodynamics, /* computation stage for target '100*theta_s' */ cs_thermodynamics, /* computation stage for target 'theta_s_100' */ cs_background, /* computation stage for target 'Omega_dcdmdr' */ diff --git a/source/perturbations.c b/source/perturbations.c index 5de8ff36..5436d6fa 100644 --- a/source/perturbations.c +++ b/source/perturbations.c @@ -6060,14 +6060,6 @@ int perturbations_initial_conditions(struct precision * ppr, +ppw->pvecback[pba->index_bg_phi_prime_scf]*alpha_prime); } - if (pba->has_smg == _TRUE_) { - class_call( - perturbations_get_x_x_prime_newtonian_smg(ppw), - ppt->error_message, - ppt->error_message - ); - } - if ((pba->has_ur == _TRUE_) || (pba->has_ncdm == _TRUE_) || (pba->has_dr == _TRUE_) || (pba->has_idr == _TRUE_)) { delta_ur -= 4.*a_prime_over_a*alpha;