Problems to check a new parameterisation for the Equation of State

[MarisolJC]

Hi!

I am using CLASS code to check a new parameterisation for the Equation of State (EoS) and I'm trying to link it with Montepython. I wrote the new equation of state in background.c together with the derivative and the corresponding integral, I also defined the new parameters of this new EoS in the background.h, in the file input.c and cclassy.pxd. When I go to the CLASS directory and compile the code using make clean, make and make-j seems that the changes works, after I type ./class explanatory.ini and again looks that everything works well, the problem is when I try to make a run with MontePython, I get the next error message:

Error in Class: background_init(L:593) :condition (pba->shooting_failed == _TRUE_) is true; Shooting failed, try optimising input_get_guess(). Error message:

input_init(L:330) :error in input_find_root(&xzero, &fevals, &fzw, errmsg); =>input_find_root(L:3937) :error in input_fzerofun_1d(x1, pfzw, &f1, errmsg); =>input_fzerofun_1d(L:3474) :error in input_try_unknown_parameters(&input, 1, pfzw, output, error_message); =>input_try_unknown_parameters(L:3651) :error in background_init(&pr,&ba); =>background_init(L:632) :condition (w_fld >= 1./3.) is true; Your choice for w(a--->0)=1.17471 is suspicious, since it is bigger than -1/3 there cannot be radiation domination at early times

# 5000 steps done, acceptance rate: 0.0 /!\ The acceptance rate is below 0.05. You might want to set the jumping factor to a lower value than the default (2.4), with the option-f 1.5for instance.

I don't get anything good, regardless of the number of steps I indicate in Montepyhton, at the end of the run, only one of them is accepted.

¿What does that errors mean? ¿How could I solve it? I'm new in this area, I'm trying to understand how the code works.

Marisol.

[Stefan-Heimersheim]

Hi Marisol,

The error message says your $w(a=0)$ is 1.17 -- does this make sense in your model?

When you run ./class explanatory.ini, did you modify explanatory.ini to include your new parameters? Otherwise this just runs the normal standard model. Or did you change the code such that it has a changed EoS without you giving extra parameters?

Edit: If ./class works but the Montepython doesn't, use this python script to test if the wrapper works for a simple case (if it doesn't output anything it worked, details here):

from classy import Class
cosmo = Class()
params = {
    'output': 'tCl lCl',
    'l_max_scalars': 2000,
    'lensing': 'yes',
    'A_s': 2.3e-9,
    'n_s': 0.9624, 
    'h': 0.6711,
    'omega_b': 0.022068,
    'omega_cdm': 0.12029}

cosmo.set(params)
cosmo.compute()

Best, Stefan

[MarisolJC]

Hi again!

Thaks for your answer. I noted that I was not modifying the explanatory.ini, in order to not to complicate things so much, now I working with a model that only use one parameter (I'm using w0_fld to doesn't add an extraparameter in class_public/Python), I set Omega_Lambda=0.0, wa_fld=0, commented Omega_fld=0.0 and modified the value of w0_fld, then I type make clean, make and make -j without problem, but when I run ./class explanatory.ini I get the following:

Reading input parameters
 -> matched budget equations by adjusting Omega_fld = 6.878622e-01
Running CLASS version v2.6.3
Computing background

Error running background_init 
=>background_init(L:682) :error in background_solve(ppr,pba);
=>background_solve(L:1630) :error in generic_integrator(background_derivs, tau_start, tau_end, pvecback_integration, &bpaw, ppr->tol_background_integration, ppr->smallest_allowed_variation, &gi);
=>generic_integrator(L:112) :error in rkqs(&x, h, eps, &hdid, &hnext, derivs, parameters_and_workspace_for_derivs, pgi);
=>rkqs(L:156) :condition (xnew == *x) is true; stepsize underflow at x=1.561831e+04

I don't understand the problem, what is that xnew== x*? Seems that there is a problem with the generic_integrator. For my model w_fld doesn't lead to a simple analytic solution of the integral necessary, so I implemented a numerical calculation of this integral using mathematica and added this result in background_initial_conditions() instead of calling background_w_fld and after setting integral_fld =0. What else could I do?

Marisol.

[Stefan-Heimersheim]

The code tries to numerically integrate (generic_integrator) the evolution of the everything in background_solve(). Lines 410 to 426 deal with dark energy (fld). They call background_w_fld() (line 519) which contains the EoS etc. and quite a few comments that might be helpful [e.g. "Obsolete: at the beginning, we had here the analytic integral solution ...", "so far, HyRec explicitely assumes that w(a)= w0 + wa (1-a/a0)"].

I'm not sure how to implement the analytic solution again (but maybe see the comment how it was done before). You could try implementing your EoS directly in background_w_fld(), use Recfast (if HyRec still assumes the standard EoS) and see if that works.

I should add that I never played around with changing the EoS so there might be things I'm missing; although you can search for "fld" in the code, most things look fairly general and probably don't need adapting.

Best, Stefan

[Uendert]

Hello @MarisolJC and @Stefan-Heimersheim

I've been experience a related issue. I had made a couple of changes in the code and it works fine when I am using a MAC computer. However, when I use the same code in Linux ubuntu I got the same error as Marisol, namely,

Error running background_init 
=>background_init(L:682) :error in background_solve(ppr,pba);
=>background_solve(L:1630) :error in generic_integrator(background_derivs, tau_start, tau_end, pvecback_integration, &bpaw, ppr->tol_background_integration, ppr->smallest_allowed_variation, &gi);
=>generic_integrator(L:112) :error in rkqs(&x, h, eps, &hdid, &hnext, derivs, parameters_and_workspace_for_derivs, pgi);
=>rkqs(L:156) :condition (xnew == *x) is true; stepsize underflow at x=1.561831e+04

Just for reference, I also included analytic integral solution to the fluid equation for the CLP parameterisation, but with wa_fld=0.

The weirdest thing is that it works fine on MAC, but fails when I use Linux Ubuntu. I've tested this in four different Mac and Linux and the same trends hold. Any thoughts?

Best, Uendert

[Uendert]

Hello @MarisolJC and @Stefan-Heimersheim

I've been experience a related issue. I had made a couple of changes in the code and it works fine when I am using a MAC computer. However, when I use the same code in Linux ubuntu I got the same error as Marisol, namely,

Error running background_init 
=>background_init(L:682) :error in background_solve(ppr,pba);
=>background_solve(L:1630) :error in generic_integrator(background_derivs, tau_start, tau_end, pvecback_integration, &bpaw, ppr->tol_background_integration, ppr->smallest_allowed_variation, &gi);
=>generic_integrator(L:112) :error in rkqs(&x, h, eps, &hdid, &hnext, derivs, parameters_and_workspace_for_derivs, pgi);
=>rkqs(L:156) :condition (xnew == *x) is true; stepsize underflow at x=1.561831e+04

Just for reference, I also included analytic integral solution to the fluid equation for the CLP parameterisation, but with wa_fld=0.

The weirdest thing is that it works fine on MAC, but fails when I use Linux Ubuntu. I've tested this in four different Mac and Linux and the same trends hold. Any thoughts?

Best, Uendert

As I was conjecturing, the main difference between Mac and Linux would be the compilers. The former uses clang while the latter uses gcc. So I installed clang on Linux and then the error just disappeared.

[MarisolJC]

Hi @Uendert!

I have installed clang on Linux, but in my case, the error message just have changed:

Reading input parameters
 -> matched budget equations by adjusting Omega_fld = 6.878622e-01
Running CLASS version v2.6.3
Computing background

Error running background_init 
=>background_init(L:691) :error in background_solve(ppr,pba);
=>background_solve(L:1608) :error in background_functions(pba,pvecback_integration, pba->short_info, pvecback);
=>background_functions(L:286) :condition (a_rel <= 0.) is true; a = -nan instead of strictly positive

Now, I don't know what is that a_rel ... :(

[Uendert]

Hi @MarisolJC,

I think it is promising now. a_rel stands for the normalized scale factor, namely, a_rel = a / a_0. Your integration is leading to negative values of the scale factor when integrating the Friedmann equation.

If I followed your messages correctly, you said you added an analytic solution into background_initial_conditions(), but you didn't mention that you had added it into background_functions(), where it defined the Friedmann equation.

So a fix could be adding the analytic solution into the background_functions() when it is taken into account the fluid contribution into the energy density budget. Just for reference, this is around

 409   /* fluid with w(a) and constant cs2 */
 410   if (pba->has_fld == _TRUE_) {
 411 
 412     /* get rho_fld from vector of integrated variables */
 413     pvecback[pba->index_bg_rho_fld] = pvecback_B[pba->index_bi_rho_fld];
 414