cmbfast help

 Posts: 13
 Joined: March 30 2005
 Affiliation: Center For Theoretical Physics, Jamia Millia Islamia, new delhi
cmbfast help
Hi,
I am trying to use cmbfast for a matter+dark energy model taking into account the dark energy perturbation. The dark energy I am studying is a perfect fluid one given by [tex]p = f(\rho)[/tex].
Now in the cmbfast version 4.5.1 one can give the dark energy equation of state input as a file with scale factor and [tex]w[/tex] as the two colomns.
My question is that when I give this file (which I generate from my dark energy equation of state) as input and choose the dark energy perturbation option in cmbfast, will it take the speed of sound = 1 as in standrad scalar field, or will it calculate the speed of sound accordingly?
If it always takes [tex]c_s^2= 1[/tex] then is there any way one can change this and incorporate the actual [tex]c_s^2[/tex] for a pefect fluid type dark energy?
Thanks
anjan sen
I am trying to use cmbfast for a matter+dark energy model taking into account the dark energy perturbation. The dark energy I am studying is a perfect fluid one given by [tex]p = f(\rho)[/tex].
Now in the cmbfast version 4.5.1 one can give the dark energy equation of state input as a file with scale factor and [tex]w[/tex] as the two colomns.
My question is that when I give this file (which I generate from my dark energy equation of state) as input and choose the dark energy perturbation option in cmbfast, will it take the speed of sound = 1 as in standrad scalar field, or will it calculate the speed of sound accordingly?
If it always takes [tex]c_s^2= 1[/tex] then is there any way one can change this and incorporate the actual [tex]c_s^2[/tex] for a pefect fluid type dark energy?
Thanks
anjan sen

 Posts: 43
 Joined: November 11 2004
 Affiliation: LUTH, Observatoire de ParisMeudon
 Contact:
cmbfast help
Hi,
the original version of CMBFAST 4.5.1 uses the input tabulated w(a) to compute an effective scalar field potential and solve the scalar field perturbed equations (namely the second order perturbed KleinGordon equation). This force the dark energy perturbations to behave as a generic fluid with an effective equation of state (I stress effective) =1. If you are interested in generic fluids where c_s^2 is a free parameter you need to change the dark energy perturbation equations.
I have implemented in CROSS_CMBFAST the equations for a generic dark energy fluid specified by w, or w(a)=w_0+w_1 a and c_s^2, which maybe is what you are looking for.
If you are interested in models for which p= w(a) \rho you just need to modify the file subcross.F, more specifically the function wdyn_func(a) (if I remember the name correctly). In that sub you can just add a spline function which interpolate your tabulated w(a), which you must have stored in a matrix before the CMBFAST computation starts. The latter can be done by adding a bunch of lines in drivercross.F before the sub cmbflat(..) is called.
Then you are free of varying c_s^2 as you like since is an input parameter of the dark energy perturbation equations which are in cmbcrossflat.F (sub fderivs).
By the way if your models are not of the form p=w(a)\rho but like p=f(\rho) remeber to modify the dark energy perturbation equations accordingly in fderivs, you don't certainly want to compute the perturbations of another model.
All the best,
Pier Stefano
the original version of CMBFAST 4.5.1 uses the input tabulated w(a) to compute an effective scalar field potential and solve the scalar field perturbed equations (namely the second order perturbed KleinGordon equation). This force the dark energy perturbations to behave as a generic fluid with an effective equation of state (I stress effective) =1. If you are interested in generic fluids where c_s^2 is a free parameter you need to change the dark energy perturbation equations.
I have implemented in CROSS_CMBFAST the equations for a generic dark energy fluid specified by w, or w(a)=w_0+w_1 a and c_s^2, which maybe is what you are looking for.
If you are interested in models for which p= w(a) \rho you just need to modify the file subcross.F, more specifically the function wdyn_func(a) (if I remember the name correctly). In that sub you can just add a spline function which interpolate your tabulated w(a), which you must have stored in a matrix before the CMBFAST computation starts. The latter can be done by adding a bunch of lines in drivercross.F before the sub cmbflat(..) is called.
Then you are free of varying c_s^2 as you like since is an input parameter of the dark energy perturbation equations which are in cmbcrossflat.F (sub fderivs).
By the way if your models are not of the form p=w(a)\rho but like p=f(\rho) remeber to modify the dark energy perturbation equations accordingly in fderivs, you don't certainly want to compute the perturbations of another model.
All the best,
Pier Stefano

 Posts: 43
 Joined: November 11 2004
 Affiliation: LUTH, Observatoire de ParisMeudon
 Contact:
cmbfast help
forgot to add that you also need to modify the function dynrho(a) in subcross.F that is the \rho_DE(a) of your model.
PS
PS

 Posts: 13
 Joined: March 30 2005
 Affiliation: Center For Theoretical Physics, Jamia Millia Islamia, new delhi
cmbfast help
Thanks a lot Pier.
I shall try this.
anjan
I shall try this.
anjan

 Posts: 13
 Joined: March 30 2005
 Affiliation: Center For Theoretical Physics, Jamia Millia Islamia, new delhi
cmbfast help
Hi Pier,
while running make for the program it is generating some errors which I listed below:
g77 ffastmath c o jlgen.o jlgen.F
jlgen.F: In program `jlgen':
jlgen.F:14:
include 'cmbfast.inc'
^
Unable to open INCLUDE file `cmbfast.inc' at (^)
jlgen.F:18:
integer l(lmax),i,j,lmo
^
Invalid declaration of or reference to symbol `lmax' at (^) [initially seen at (^)]
jlgen.F:18:
integer l(lmax),i,j,lmo
1
jlgen.F:21: (continued):
common /lvalues1/ l,l0,lmo
2
Invalid declaration of or reference to symbol `l' at (2) [initially seen at (1)]jlgen.F: In subroutine `initlval':
jlgen.F:121:
include 'cmbfast.inc'
^
Unable to open INCLUDE file `cmbfast.inc' at (^)
jlgen.F:123:
integer l(lmax),l0
1
jlgen.F:125: (continued):
common /lvalues1/ l,l0,lmo
2
Invalid declaration of or reference to symbol `l' at (2) [initially seen at (1)]jlgen.F:133:
include 'lvalues.inc'
^
Unable to open INCLUDE file `lvalues.inc' at (^)
jlgen.F:123:
integer l(lmax),l0
1
jlgen.F:136: (continued):
end
2
Invalid declaration of or reference to symbol `lmax' at (2) [initially seen at (1)]
jlgen.F: In subroutine `bjl':
jlgen.F:96: warning:
call bjl(l(j),x,ajl)
1
jlgen.F:140: (continued):
subroutine bjl(L,x,jl)
2
Argument #1 (named `l') of `bjl' is one type at (2) but is some other type at (1) [info f g77 M GLOBALS]
make: *** [jlgen.o] Error 1
can you help me in this regard?
anjan
while running make for the program it is generating some errors which I listed below:
g77 ffastmath c o jlgen.o jlgen.F
jlgen.F: In program `jlgen':
jlgen.F:14:
include 'cmbfast.inc'
^
Unable to open INCLUDE file `cmbfast.inc' at (^)
jlgen.F:18:
integer l(lmax),i,j,lmo
^
Invalid declaration of or reference to symbol `lmax' at (^) [initially seen at (^)]
jlgen.F:18:
integer l(lmax),i,j,lmo
1
jlgen.F:21: (continued):
common /lvalues1/ l,l0,lmo
2
Invalid declaration of or reference to symbol `l' at (2) [initially seen at (1)]jlgen.F: In subroutine `initlval':
jlgen.F:121:
include 'cmbfast.inc'
^
Unable to open INCLUDE file `cmbfast.inc' at (^)
jlgen.F:123:
integer l(lmax),l0
1
jlgen.F:125: (continued):
common /lvalues1/ l,l0,lmo
2
Invalid declaration of or reference to symbol `l' at (2) [initially seen at (1)]jlgen.F:133:
include 'lvalues.inc'
^
Unable to open INCLUDE file `lvalues.inc' at (^)
jlgen.F:123:
integer l(lmax),l0
1
jlgen.F:136: (continued):
end
2
Invalid declaration of or reference to symbol `lmax' at (2) [initially seen at (1)]
jlgen.F: In subroutine `bjl':
jlgen.F:96: warning:
call bjl(l(j),x,ajl)
1
jlgen.F:140: (continued):
subroutine bjl(L,x,jl)
2
Argument #1 (named `l') of `bjl' is one type at (2) but is some other type at (1) [info f g77 M GLOBALS]
make: *** [jlgen.o] Error 1
can you help me in this regard?
anjan

 Posts: 13
 Joined: March 30 2005
 Affiliation: Center For Theoretical Physics, Jamia Millia Islamia, new delhi
cmbfast help
Hi pier,
I have solved the previous problem.
Now I have checked in your subcross.F file in the cs2_ad(a) subroutine that it is calculating the adiabatic sound speed for the dark energy perturbation given the w(a).
But while running, what is the cs^2 it is asking for?
regards
anjan
I have solved the previous problem.
Now I have checked in your subcross.F file in the cs2_ad(a) subroutine that it is calculating the adiabatic sound speed for the dark energy perturbation given the w(a).
But while running, what is the cs^2 it is asking for?
regards
anjan

 Posts: 43
 Joined: November 11 2004
 Affiliation: LUTH, Observatoire de ParisMeudon
 Contact:
cmbfast help
One is the adiabatic sound speed c_a^2 which entirely determined by w(a) and the other is the effective sound speed, the two are differently defined. If you like the latter parameterize how close the perturbations are to a pure canonical scalar field case (for which c_eff^2=1). You can have a better insight looking at the literature, see for instance Bean and Dore' astroph/0307100.
PS
PS

 Posts: 13
 Joined: March 30 2005
 Affiliation: Center For Theoretical Physics, Jamia Millia Islamia, new delhi
cmbfast help
Yes I understand.
But if one takes a perfect fluid for dark energy then there is no dissipative effects and in principle these two should be equal. Then in this case how can one modify the code?
anjan
But if one takes a perfect fluid for dark energy then there is no dissipative effects and in principle these two should be equal. Then in this case how can one modify the code?
anjan

 Posts: 43
 Joined: November 11 2004
 Affiliation: LUTH, Observatoire de ParisMeudon
 Contact:
cmbfast help
In this case you can impose c_e^2=c_a^2 by changing accordingly the dark energy perturbation equations in the subroutine fderivs in cmbcrossflat.F. Simply substitute cs2l with c2adia at line 2275 in dgpres_phi, at line
2453 in yprime(n1) and 2457 in yprime(n).
Let me add a comment: in this case the dark energy sound speed is only the adiabatic one, hence w(a) specifies everything and the evolution of the perturbations is equivalent to solving the perturbed KG equation for a scalar field (as in in the original version of CMBFAST) with adiabatic initial conditions.
For this you can just use the original version of CMBFAST4.5.1, but be sure to properly modify the scalar field initial conditions in the subroutine finitial(..). If I remember [tex]\delta{\phi}[/tex] and [tex]\delta\dot{\varphi}[/tex] are set =0 by default, but these are not adiabatic initial conditions for the scalar field perturbations as you may realize by writing the relative entropy perturbation S_{\phi i} (with i being whatever other cosmic fluid) and the intrinsic scalar field entropy perturbation \Gamma_{\phi} in terms of \delta\phi and [tex]\delta\dot\varphi[/tex].
Hope this clarify the point.
PS
2453 in yprime(n1) and 2457 in yprime(n).
Let me add a comment: in this case the dark energy sound speed is only the adiabatic one, hence w(a) specifies everything and the evolution of the perturbations is equivalent to solving the perturbed KG equation for a scalar field (as in in the original version of CMBFAST) with adiabatic initial conditions.
For this you can just use the original version of CMBFAST4.5.1, but be sure to properly modify the scalar field initial conditions in the subroutine finitial(..). If I remember [tex]\delta{\phi}[/tex] and [tex]\delta\dot{\varphi}[/tex] are set =0 by default, but these are not adiabatic initial conditions for the scalar field perturbations as you may realize by writing the relative entropy perturbation S_{\phi i} (with i being whatever other cosmic fluid) and the intrinsic scalar field entropy perturbation \Gamma_{\phi} in terms of \delta\phi and [tex]\delta\dot\varphi[/tex].
Hope this clarify the point.
PS

 Posts: 13
 Joined: March 30 2005
 Affiliation: Center For Theoretical Physics, Jamia Millia Islamia, new delhi
cmbfast help
Hi pier,
I just tried to change the cs2l by c2adia in the lines you suggested. It makes the dark energy being a perfect fluid completely determined by its eos. I then used rest of your conditions that is the eos parametrization, input variables etc etc exactly the same. But now running the code is not successful as it the dverk function generates errors aparently due to some blowup. I am copying the message:
computation stopped in dverk with the followingvalues 
ind = 3 tol = 4.599721E10 x = 1.192250269439929E+04
n = 43 hmin = 1.000000E08 xend = 1.193284983313287E+04
nw = 451 hmax = 2.000000E+00 prev xend = 1.191905854892556E+04
no of successful steps = 13608.
no of successive failures = 25.
no of function evals = 109411.
the components of y are
0.875591193635206E+00 NAN INF 0.486776136956397+302 0.000000000000000E+00
0.486776135995405+302 0.575591201110352+292 0.560556019509850+302 0.402089594204575+300 INF
0.666915840569582+301 0.264843740374802+301 0.107621821176573+301 0.444249500465392+300 0.185367882478533+300
0.781958390964206+299 0.328098351305341+299 0.146813690766940+299 0.478611438405599+298 0.526225602102048+298
0.711832717182278+296 0.130135132888999+296 0.182856633063571+296 0.703064335912279+295 0.280627382892827+295
0.114619230203352+295 0.475558556171797+294 0.199443366060634+294 0.845713413770606+293 0.356548169422175+293
0.160531002932601+293 0.524220704522618+292 0.580731822887188+292 0.560518164022928+302 0.402261629273723+300
INF 0.668105614780505+301 0.262519675558030+301 0.112315377983086+301 0.348021316543802+300
0.385270502704955+300 INF INF
do you have any clue?
sorry for bothering you much but I am just new to these things and need help.
Thanks
anjan sen
I just tried to change the cs2l by c2adia in the lines you suggested. It makes the dark energy being a perfect fluid completely determined by its eos. I then used rest of your conditions that is the eos parametrization, input variables etc etc exactly the same. But now running the code is not successful as it the dverk function generates errors aparently due to some blowup. I am copying the message:
computation stopped in dverk with the followingvalues 
ind = 3 tol = 4.599721E10 x = 1.192250269439929E+04
n = 43 hmin = 1.000000E08 xend = 1.193284983313287E+04
nw = 451 hmax = 2.000000E+00 prev xend = 1.191905854892556E+04
no of successful steps = 13608.
no of successive failures = 25.
no of function evals = 109411.
the components of y are
0.875591193635206E+00 NAN INF 0.486776136956397+302 0.000000000000000E+00
0.486776135995405+302 0.575591201110352+292 0.560556019509850+302 0.402089594204575+300 INF
0.666915840569582+301 0.264843740374802+301 0.107621821176573+301 0.444249500465392+300 0.185367882478533+300
0.781958390964206+299 0.328098351305341+299 0.146813690766940+299 0.478611438405599+298 0.526225602102048+298
0.711832717182278+296 0.130135132888999+296 0.182856633063571+296 0.703064335912279+295 0.280627382892827+295
0.114619230203352+295 0.475558556171797+294 0.199443366060634+294 0.845713413770606+293 0.356548169422175+293
0.160531002932601+293 0.524220704522618+292 0.580731822887188+292 0.560518164022928+302 0.402261629273723+300
INF 0.668105614780505+301 0.262519675558030+301 0.112315377983086+301 0.348021316543802+300
0.385270502704955+300 INF INF
do you have any clue?
sorry for bothering you much but I am just new to these things and need help.
Thanks
anjan sen

 Posts: 43
 Joined: November 11 2004
 Affiliation: LUTH, Observatoire de ParisMeudon
 Contact:
cmbfast help
It looks like a problem of error convergence, you can try increasing the accuracy of the integrator. If you have used a w(a) which crosses 1, wierd things can happen as the code is not optimized for it, well the equations are simply not valid for such cases. I ll just try to check what happen when setting cs2l=c2adia. Thanks for pointing out this bug.
Regards,
PS
Regards,
PS

 Posts: 13
 Joined: March 30 2005
 Affiliation: Center For Theoretical Physics, Jamia Millia Islamia, new delhi
cmbfast help
Dear Pier,
I used the same wvar.ini and wcons.ini that you provided. The only thing I did is to replace cs2l by c2adia.
Please let me know if you find anything.
regards
anjan
I used the same wvar.ini and wcons.ini that you provided. The only thing I did is to replace cs2l by c2adia.
Please let me know if you find anything.
regards
anjan