Guide to EFTCAMB_Oct14


There might exist some typos in this guide, please check it carefully with our original source codes. If you encounter any problems during your modifications, please contact us.



----------------------------------------------------------------------------------------------------------------------------
Add: EFT_main.f90
----------------------------------------------------------------------------------------------------------------------------
Add: EFTstabilitySpace.f90
----------------------------------------------------------------------------------------------------------------------------
Modify: Makefile_main

Replace:
EQUATIONS ?= equations_EFT
#… Original:
#EQUATIONS ?= equations

Add: After #DRIVER ?= tester.f90
# EFTCAMB: optional program to test stability.
#DRIVER ?= EFTstabilitySpace_V3.f90


Replace:
# EFTCAMB MOD START
CAMBOBJ = constants.o utils.o $(EXTCAMBFILES) subroutines.o inifile.o $(POWERSPECTRUM).o $(RECOMBINATION).o \
$(REIONIZATION).o modules.o bessels.o EFT_main.o $(EQUATIONS).o $(NONLINEAR).o lensing.o $(BISPECTRUM).o cmbmain.o camb.o
# EFTCAMB MOD END

#… Original:
#CAMBOBJ = constants.o utils.o $(EXTCAMBFILES) subroutines.o inifile.o $(POWERSPECTRUM).o $(RECOMBINATION).o \
# $(REIONIZATION).o modules.o bessels.o $(EQUATIONS).o $(NONLINEAR).o lensing.o $(BISPECTRUM).o cmbmain.o camb.o


----------------------------------------------------------------------------------------------------------------------------
Modify: params.ini

Add: At the beginning
####### Model selection flags for EFTCAMB #######

# EFT flags: set up in which mode EFTCAMB is used.
# We refer to the documentation (EFTCAMB:numerical notes) for a thorough
# explanation of the effect of these flags.
#
# 1) Main EFT flag:
# EFTflag = 0 : GR code. Every EFT modification is ignored.
# EFTflag = 1 : Pure EFT code.
# EFTflag = 2 : Designer matching EFT.
# EFTflag = 3 : Matching EFT.

EFTflag = 0

# 2) Background Dark Energy equation of state flag:
# EFTwDE = 0 : Cosmological constant
# EFTwDE = 1 : DE with constant Eos determined by EFTw0
# EFTwDE = 2 : CPL parametrization
# EFTwDE = 3 : JBP parametrization
# EFTwDE = 4 : turning point parametrization
# EFTwDE = 5 : Taylor expansion
# EFTwDE = 6 : User defined

EFTwDE = 0

# 3) Pure EFT model selection flag:
# The following structure applies to every operator and can be specified for every operator
# separately.
#
# PureEFTmodel___ = 0 : Zero (operator ignored)
# PureEFTmodel___ = 1 : Constant model
# PureEFTmodel___ = 2 : Linear model
# PureEFTmodel___ = 3 : Power law model
# PureEFTmodel___ = 4 : Exponential model
# PureEFTmodel___ = 5 : User defined

PureEFTmodelOmega = 1
PureEFTmodelAlpha1 = 0
PureEFTmodelAlpha2 = 0
PureEFTmodelAlpha3 = 0
PureEFTmodelAlpha4 = 0
PureEFTmodelAlpha5 = 0
PureEFTmodelAlpha6 = 0

# 4) Designer matching EFT model selection flag:
# DesignerEFTmodel = 1 : designer f(R)
# DesignerEFTmodel = 2 : designer minimally coupled quintessence

DesignerEFTmodel = 1

# 5) Matching EFT model selection flag:
# No model implemented so far.

MatchingEFTmodel = 1

####### Model parameters for EFTCAMB #######

# Notice that if the model is not selected via the model selection flags then
# the values of the parameters are automatically ignored.

# 1) Background Dark Energy equation of state parameters:

EFTw0 = -0.7
EFTwa = -0.3
EFTwn = 2
EFTwat = 0.8
EFtw2 = 0.1
EFTw3 = 0.1

# 2) Pure EFT parameters:

EFTOmega0 = 0.1
EFTOmegaExp = 2.8

EFTAlpha10 = 0
EFTAlpha1Exp = 0

EFTAlpha20 = 0.1
EFTAlpha2Exp = 0

EFTAlpha30 = 0.001
EFTAlpha3Exp = 0

EFTAlpha40 = 0.001
EFTAlpha4Exp = 0

EFTAlpha50 = 0.1
EFTAlpha5Exp = 0

EFTAlpha60 = 0.1
EFTAlpha6Exp = 0

# 3) Designer matching parameters:
# Model 1: designer f(R) theories

EFTB0 = 0.1

####### Parameters for CAMB #######

----------------------------------------------------------------------------------------------------------------------------
Modify: cmbmain.f90

::::: Module CAMBmain :::::

Add: After use Errors
! EFTCAMB MOD START
use EFTinitialization
! EFTCAMB MOD END

Replace:
! EFTCAMB MOD START:
! To debug the code it should be launched on just one thread by
! changing the appropriate entry in the parameter file.
real(dl) :: fixq = 0._dl !Debug output of one q
! EFTCAMB MOD END.
#… Original:
real(dl) :: fixq = 0._dl !Debug output of one q

::::: subroutine cmbmain :::::

Add: After integer q_ix
! EFTCAMB MOD START
logical EFTsuccess
! EFTCAMB MOD END

Add: After WantLateTime = CP%DoLensing .or. num_redshiftwindows > 0
! EFTCAMB MOD START
if (EFTCAMBuseinCOSMOMC==1.and.CP%EFTflag/=0) then
call EFTCAMB_initialization(EFTsuccess)
if (.not.EFTsuccess) stop
end if
! EFTCAMB MOD END

::::: subroutine CalcScalarSources(EV,taustart) :::::

Replace:
! EFTCAMB MOD START: It's always nice to plot the evolution of perturbations in DE/MG models!
! In order to print the user needs to choose the scale by means of fixq and then turno on the
! flag that prints quantities from the output routine.
! It is suggested to run the code on just one core.

if (fixq/=0._dl) then
tol1=tol/exp(AccuracyBoost-1)

write(*,*) 'EFTCAMB: start printing.'
call CreateTxtFile('Results/Debug_Evolution/Files/1_FRW.dat',1)
call CreateTxtFile('Results/Debug_Evolution/Files/2_EFTfunctions.dat',2)
call CreateTxtFile('Results/Debug_Evolution/Files/3_EFTBackground.dat',3)
call CreateTxtFile('Results/Debug_Evolution/Files/4_PiFieldSolution.dat',4)
call CreateTxtFile('Results/Debug_Evolution/Files/5_MetricSolution.dat',5)
call CreateTxtFile('Results/Debug_Evolution/Files/6_DensitySolution.dat',6)
call CreateTxtFile('Results/Debug_Evolution/Files/7_EinsteinEqFactors.dat',7)
call CreateTxtFile('Results/Debug_Evolution/Files/8_PiEqFactors.dat',8)
call CreateTxtFile('Results/Debug_Evolution/Files/9_Sources.dat',9)
call CreateTxtFile('Results/Debug_Evolution/Files/10_MuGamma.dat',10)

call CreateTxtFile('Results/Debug_File/Debug.dat',11)

do j=1,10000
tauend = taustart +REAL(j-1)*(CP%tau0-taustart)/REAL(10000-1)
call GaugeInterface_EvolveScal(EV,tau,y,tauend,tol1,ind,c,w)
yprime = 0
call derivs(EV,EV%ScalEqsToPropagate,tau,y,yprime)
call output(EV,y,j,tau,sources)
end do

close(1)
close(2)
close(3)
close(4)
close(5)
close(6)
close(7)
close(8)
close(9)
close(10)
close(11)
write(*,*) 'Stop printing'
stop
end if

! Original code:
!!Example code for plotting out variable evolution
!!if (fixq/=0._dl) then
!! tol1=tol/exp(AccuracyBoost-1)
!! call CreateTxtFile('evolve_q005.txt',1)
!! do j=1,1000
!! tauend = taustart+(j-1)*(CP%tau0-taustart)/1000
!! call GaugeInterface_EvolveScal(EV,tau,y,tauend,tol1,ind,c,w)
!! yprime = 0
!! call derivs(EV,EV%ScalEqsToPropagate,tau,y,yprime)
!! adotoa = 1/(y(1)*dtauda(y(1)))
!! ddelta= (yprime(3)*grhoc+yprime(4)*grhob)/(grhob+grhoc)
!! delta=(grhoc*y(3)+grhob*y(4))/(grhob+grhoc)
!! growth= ddelta/delta/adotoa
!! write (1,'(7E15.5)') tau, delta, growth, y(3), y(4), y(EV%g_ix), y(1)
!! end do
!! close(1)
!! stop
!!end if
! EFTCAMB MOD END

::::: subroutine GetTransfer(EV,tau) :::::

Replace:
! EFTCAMB MOD START: this seems to create some problems. It should be safe to turn it off.
if (CP%EFTflag/=0) then
if (CP%Transfer%high_precision) atol=atol
else
if (CP%Transfer%high_precision) atol=atol/10000
end if

! Original code:
! if (CP%Transfer%high_precision) atol=atol/10000
! EFTCAMB MOD END.

----------------------------------------------------------------------------------------------------------------------------
Modify: modules.f90

::::: module ModelParams :::::

Add: After integer :: Nu_mass_numbers(max_nu) !physical number per eigenstate
! EFTCAMB MOD START
! 1) Definition of flags:
integer :: EFTflag
integer :: EFTwDE
integer :: PureEFTmodelOmega
integer :: PureEFTmodelAlpha1, PureEFTmodelAlpha2, PureEFTmodelAlpha3
integer :: PureEFTmodelAlpha4, PureEFTmodelAlpha5, PureEFTmodelAlpha6
integer :: DesignerEFTmodel
integer :: MatchingEFTmodel
! 2) Definition of model parameters:
real(dl) :: EFTw0, EFTwa, EFTwn, EFTwat, EFtw2, EFTw3
real(dl) :: EFTOmega0, EFTOmegaExp
real(dl) :: EFTAlpha10, EFTAlpha1Exp, EFTAlpha20, EFTAlpha2Exp
real(dl) :: EFTAlpha30, EFTAlpha3Exp, EFTAlpha40, EFTAlpha4Exp
real(dl) :: EFTAlpha50, EFTAlpha5Exp, EFTAlpha60, EFTAlpha6Exp
real(dl) :: EFTB0
! EFTCAMB MOD END

::::: module MassiveNu :::::

Replace:
! EFTCAMB MOD START: compatibility with massive neutrinos
real(dl), dimension(:), allocatable :: r1,p1,dr1,dp1,ddr1, ddp1, dddp1
! Original code:
! real(dl), dimension(:), allocatable :: r1,p1,dr1,dp1,ddr1
! EFTCAMB MOD END.

Replace:
! EFTCAMB MOD START: compatibility with massive neutrinos
public const,Nu_Init,Nu_background, Nu_rho, Nu_drho, nqmax0, nqmax, &
nu_int_kernel, nu_q, Nu_pidot, Nu_pidotdot
! Original code:
! public const,Nu_Init,Nu_background, Nu_rho, Nu_drho, nqmax0, nqmax, &
! nu_int_kernel, nu_q
! EFTCAMB MOD END.

::::: subroutine Nu_init :::::

Replace:
! EFTCAMB MOD START: compatibility with massive neutrinos
allocate(r1(nrhopn),p1(nrhopn),dr1(nrhopn),dp1(nrhopn),ddr1(nrhopn),ddp1(nrhopn),dddp1(nrhopn))
! Original code:
! allocate(r1(nrhopn),p1(nrhopn),dr1(nrhopn),dp1(nrhopn),ddr1(nrhopn))
! EFTCAMB MOD END.

Add: After call splder(dr1,ddr1,nrhopn,spline_data)
! EFTCAMB MOD START: compatibility with massive neutrinos
call splder(dp1,ddp1,nrhopn,spline_data)
call splder(ddp1,dddp1,nrhopn,spline_data)
! EFTCAMB MOD END.

Add: After end function Nu_drho
! EFTCAMB MOD START: compatibility with massive neutrinos
function Nu_pidot(am,adotoa,presnu) result (presnudot)
use precision
use ModelParams

real(dl) adotoa,presnu,presnudot
real(dl) d
real(dl), intent(IN) :: am
integer i

if (am< am_minp) then
presnudot = -2*const2*am**2*adotoa/3._dl
else if (am>am_maxp) then
presnudot = -((15._dl*(4._dl*am**2*zeta5 -189._dl*Zeta7))/(8._dl*am**3*const))*adotoa
else
d=log(am/am_min)/dlnam+1._dl
i=int(d)
d=d-i

presnudot = dp1(i)+d*(ddp1(i)+d*(3._dl*(dp1(i+1)-dp1(i))-2._dl*ddp1(i) &
-ddp1(i+1)+d*(ddp1(i)+ddp1(i+1)+2._dl*(dp1(i)-dp1(i+1)))))

presnudot=presnu*adotoa*presnudot/dlnam
end if

end function Nu_pidot

function Nu_pidotdot(am,adotoa,Hdot,presnu,presnudot) result (presnudotdot)
use precision
use ModelParams

real(dl) adotoa,Hdot,presnu,presnudot,presnudotdot
real(dl) d
real(dl), intent(in) :: am
integer i

if (am< am_minp) then
presnudotdot = presnudot*(adotoa +Hdot/adotoa) +am**2*adotoa**2*(-2._dl*const2/3._dl)
else if (am>am_maxp) then
presnudotdot = presnudot*(adotoa +Hdot/adotoa) +am**2*adotoa**2*(&
&-((15._dl*zeta5)/(am**3*const)) + (15._dl*(4._dl*am**2*zeta5 -189._dl*Zeta7))/(2._dl*am**5*const))
else

d=log(am/am_min)/dlnam+1._dl
i=int(d)
d=d-i

presnudotdot = ddp1(i)+d*(dddp1(i)+d*(3._dl*(ddp1(i+1)-ddp1(i))-2._dl*dddp1(i) &
-dddp1(i+1)+d*(dddp1(i)+dddp1(i+1)+2._dl*(ddp1(i)-ddp1(i+1)))))

presnudotdot = +adotoa**2*presnu*presnudotdot/dlnam +Hdot/adotoa*presnudot +presnudot**2/presnu
end if

end function Nu_pidotdot
! EFTCAMB MOD END.

----------------------------------------------------------------------------------------------------------------------------
Modify: inidriver.F90

::::: program driver :::::

Add: After
#ifdef NAGF95
use F90_UNIX
#endif
! EFTCAMB MOD START
use EFTdef
! EFTCAMB MOD END

Add: After P%tcmb = Ini_Read_Double('temp_cmb',COBE_CMBTemp)
! EFTCAMB MOD START

! 1) Initialization of EFTCAMB flags.

P%EFTflag = Ini_Read_Int('EFTflag',0)

P%EFTwDE = Ini_Read_Int('EFTwDE',0)

P%PureEFTmodelOmega = Ini_Read_Int('PureEFTmodelOmega',0)
P%PureEFTmodelAlpha1 = Ini_Read_Int('PureEFTmodelAlpha1',0)
P%PureEFTmodelAlpha2 = Ini_Read_Int('PureEFTmodelAlpha2',0)
P%PureEFTmodelAlpha3 = Ini_Read_Int('PureEFTmodelAlpha3',0)
P%PureEFTmodelAlpha4 = Ini_Read_Int('PureEFTmodelAlpha4',0)
P%PureEFTmodelAlpha5 = Ini_Read_Int('PureEFTmodelAlpha5',0)
P%PureEFTmodelAlpha6 = Ini_Read_Int('PureEFTmodelAlpha6',0)

P%DesignerEFTmodel = Ini_Read_Int('DesignerEFTmodel',1)
P%MatchingEFTmodel = Ini_Read_Int('MatchingEFTmodel',1)

! 2) Initialization of EFTCAMB model parameters.

P%EFTw0 = Ini_Read_Double('EFTw0',-1._dl)
P%EFTwa = Ini_Read_Double('EFTwa',0._dl)
P%EFTwn = Ini_Read_Double('EFTwn',2._dl)
P%EFTwat = Ini_Read_Double('EFTwat',1._dl)
P%EFtw2 = Ini_Read_Double('EFtw2',0._dl)
P%EFTw3 = Ini_Read_Double('EFTw3',0._dl)

P%EFTOmega0 = Ini_Read_Double('EFTOmega0')
P%EFTOmegaExp = Ini_Read_Double('EFTOmegaExp')
P%EFTAlpha10 = Ini_Read_Double('EFTAlpha10')
P%EFTAlpha1Exp = Ini_Read_Double('EFTAlpha1Exp')
P%EFTAlpha20 = Ini_Read_Double('EFTAlpha20')
P%EFTAlpha2Exp = Ini_Read_Double('EFTAlpha2Exp')
P%EFTAlpha30 = Ini_Read_Double('EFTAlpha30')
P%EFTAlpha3Exp = Ini_Read_Double('EFTAlpha3Exp')
P%EFTAlpha40 = Ini_Read_Double('EFTAlpha40')
P%EFTAlpha4Exp = Ini_Read_Double('EFTAlpha4Exp')
P%EFTAlpha50 = Ini_Read_Double('EFTAlpha50')
P%EFTAlpha5Exp = Ini_Read_Double('EFTAlpha5Exp')
P%EFTAlpha60 = Ini_Read_Double('EFTAlpha60')
P%EFTAlpha6Exp = Ini_Read_Double('EFTAlpha6Exp')

P%EFTB0 = Ini_Read_Double('EFTB0')

! EFTCAMB MOD END

----------------------------------------------------------------------------------------------------------------------------
Modify: equations.f90 => equations_EFT.f90

::::: function dtauda(a) :::::

Add: After use LambdaGeneral
! EFTCAMB MOD START
use EFTfunctions
! EFTCAMB MOD END

Replace:
! EFTCAMB MOD START
if (CP%EFTflag==0) then ! Original CAMB code.
if (w_lam == -1._dl) then
grhoa2=grhoa2+grhov*a2**2
else
grhoa2=grhoa2+grhov*a**(1-3*w_lam)
end if
else if (CP%EFTflag/=0) then
grhoa2=grhoa2+grhov*EFTw(a,3)*a2
end if
! EFTCAMB MOD END

::::: module GaugeInterface :::::

Add: After use Errors
! EFTCAMB MOD START
use EFTfunctions
! EFTCAMB MOD END

Add: After real(dl) denlkt(4,max_l_evolve),Kft(max_l_evolve)
! EFTCAMB MOD START
logical EFTCAMBactive
real(dl) EFTturnOnTime
! EFTCAMB MOD END

::::: recursive subroutine GaugeInterface_EvolveScal(EV,tau,y,tauend,tol1,ind,c,w) :::::

Add: After tau_switch_nu_massive= noSwitch

! EFTCAMB MOD START: add the dark energy switch.
if (CP%EFTflag==0) EV%EFTturnOnTime = noSwitch
! EFTCAMB MOD END

Replace:
! EFTCAMB MOD START: add the dark energy switch.
next_switch = min(tau_switch_ktau, tau_switch_nu_massless,EV%TightSwitchoffTime, tau_switch_nu_massive, &
tau_switch_no_nu_multpoles, tau_switch_no_phot_multpoles, tau_switch_nu_nonrel,noSwitch, EV%EFTturnOnTime)
! Original code:
! next_switch = min(tau_switch_ktau, tau_switch_nu_massless,EV%TightSwitchoffTime, tau_switch_nu_massive, &
! tau_switch_no_nu_multpoles, tau_switch_no_phot_multpoles, tau_switch_nu_nonrel,noSwitch)
! EFTCAMB MOD END

Add: Inside else if (next_switch==tau_switch_no_phot_multpoles) then After EV=EVout
! EFTCAMB MOD START: effect of the DE switch.
else if (next_switch==EV%EFTturnOnTime) then
ind=1
EVout%EFTCAMBactive = .true.
EVout%EFTturnOnTime = noSwitch
call SetupScalarArrayIndices(EVout)
call CopyScalarVariableArray(y,yout, EV, EVout)
y=yout
EV=EVout
call EFTpiInitialConditions(y,EV,tau)
! EFTCAMB MOD STOP

::::: subroutine SetupScalarArrayIndices(EV, max_num_eqns) :::::

Replace:
! EFTCAMB MOD START
if ((w_lam /= -1 .and. w_Perturb).or.CP%EFTflag/=0) then
EV%w_ix = neq+1
neq=neq+2
maxeq=maxeq+2
else
EV%w_ix=0
end if
! Original code:
! !Dark energy
! if (w_lam /= -1 .and. w_Perturb) then
! EV%w_ix = neq+1
! neq=neq+2
! maxeq=maxeq+2
! else
! EV%w_ix=0
! end if
! EFTCAMB MOD END

::::: subroutine CopyScalarVariableArray(y,yout, EV, EVout) :::::

Replace:
! EFTCAMB MOD START
if ((w_lam /= -1.and.w_Perturb).or.&
&(CP%EFTflag/=0.and.EVout%EFTCAMBactive).or.&
&(CP%EFTflag/=0.and.EV%EFTCAMBactive)) then

yout(EVout%w_ix)=y(EV%w_ix)
yout(EVout%w_ix+1)=y(EV%w_ix+1)
end if
! Original code:
! if (w_lam /= -1 .and. w_Perturb) then
! yout(EVout%w_ix)=y(EV%w_ix)
! yout(EVout%w_ix+1)=y(EV%w_ix+1)
! end if
! EFTCAMB MOD END

::::: subroutine GetNumEqns(EV) :::::

Add: After EV%MassiveNuApprox=.false.
! EFTCAMB MOD START:
EV%EFTCAMBactive = .false.
! EFTCAMB MOD END

::::: subroutine MassiveNuVars(EV,y,a,grho,gpres,dgrho,dgq, wnu_arr) :::::

Replace:
! EFTCAMB MOD START: compatibility with massive neutrinos
subroutine MassiveNuVars(EV,y,a,grho,gpres,dgrho,dgq, wnu_arr, dgp)
! Original code:
! subroutine MassiveNuVars(EV,y,a,grho,gpres,dgrho,dgq, wnu_arr)
! EFTCAMB MOD END.
implicit none
type(EvolutionVars) EV
real(dl) :: y(EV%nvar), a, grho,gpres,dgrho,dgq
real(dl), intent(out), optional :: wnu_arr(max_nu)
! EFTCAMB MOD START: compatibility with massive neutrinos
real(dl), intent(out), optional :: dgp
! EFTCAMB MOD END.

Add: After dgq = dgq + grhonu_t*qnu
! EFTCAMB MOD START: compatibility with massive neutrinos
if (present(dgp)) then
dgp = dgp + gpnu_t*clxnu
end if
! EFTCAMB MOD END.

-1-