diff --git a/camb/cambdll.dll b/camb/cambdll.dll index 10c163d0..7d8659be 100644 Binary files a/camb/cambdll.dll and b/camb/cambdll.dll differ diff --git a/camb/dark_energy.py b/camb/dark_energy.py index e2e0c257..04c1c174 100644 --- a/camb/dark_energy.py +++ b/camb/dark_energy.py @@ -105,17 +105,18 @@ class AxionEffectiveFluid(DarkEnergyModel): Not well tested, but should serve to demonstrate how to make your own custom classes. """ _fields_ = [ - ("w_n", c_double), - ("om", c_double), - ("a_c", c_double), - ("theta_i", c_double)] + ("w_n", c_double, "effective equation of state parameter"), + ("fde_zc", c_double, "energy density fraction at z=zc"), + ("zc", c_double, "decay transition redshift (not same as peak of energy density fraction)"), + ("theta_i", c_double, "initial condition field value")] + _fortran_class_name_ = 'TAxionEffectiveFluid' _fortran_class_module_ = 'DarkEnergyFluid' - def set_params(self, w_n, om, a_c, theta_i=None): + def set_params(self, w_n, fde_zc, zc, theta_i=None): self.w_n = w_n - self.om = om - self.a_c = a_c + self.fde_zc = fde_zc + self.zc = zc if theta_i is not None: self.theta_i = theta_i diff --git a/camb/tests/camb_test.py b/camb/tests/camb_test.py index f168365e..1dbd02a2 100644 --- a/camb/tests/camb_test.py +++ b/camb/tests/camb_test.py @@ -92,7 +92,7 @@ def testAssigments(self): 'cosmomc_theta', 'YHe', 'wa', 'cs2', 'H0', 'mnu', 'Alens', 'TCMB', 'ns', 'nrun', 'As', 'nt', 'r', 'w', 'omch2'}) params2 = camb.get_valid_numerical_params(dark_energy_model='AxionEffectiveFluid') - self.assertEqual(params2.difference(params), {'om', 'w_n', 'a_c', 'theta_i'}) + self.assertEqual(params2.difference(params), {'fde_zc', 'w_n', 'zc', 'theta_i'}) def testBackground(self): pars = camb.CAMBparams() @@ -704,5 +704,5 @@ def test_memory(self): def test_quintessence(self): pars = camb.set_params( ombh2=0.022, omch2=0.122, - dark_energy_model='Quintessence', H0=67) + dark_energy_model='EarlyQuintessence', H0=67) camb.get_background(pars) diff --git a/fortran/DarkEnergyFluid.f90 b/fortran/DarkEnergyFluid.f90 index 89e8fef3..0849e790 100644 --- a/fortran/DarkEnergyFluid.f90 +++ b/fortran/DarkEnergyFluid.f90 @@ -23,10 +23,12 @@ module DarkEnergyFluid !This is an example, it's not supposed to be a rigorous model! (not very well tested) type, extends(TDarkEnergyModel) :: TAxionEffectiveFluid real(dl) :: w_n = 1._dl !Effective equation of state when oscillating - real(dl) :: Om = 0._dl !Omega of the early DE component today (assumed to be negligible compared to omega_lambda) - real(dl) :: a_c !transition scale factor + real(dl) :: fde_zc = 0._dl ! energy density fraction at a_c (not the same as peak dark energy fraction) + real(dl) :: zc !transition redshift (scale factor a_c) real(dl) :: theta_i = const_pi/2 !Initial value - real(dl), private :: pow, omL, acpow, freq, n !cached internally + !om is Omega of the early DE component today (assumed to be negligible compared to omega_lambda) + !omL is the lambda component of the total dark energy omega + real(dl), private :: a_c, pow, om, omL, acpow, freq, n !cached internally contains procedure :: ReadParams => TAxionEffectiveFluid_ReadParams procedure, nopass :: PythonClass => TAxionEffectiveFluid_PythonClass @@ -151,9 +153,12 @@ subroutine TAxionEffectiveFluid_ReadParams(this, Ini) class(TIniFile), intent(in) :: Ini call this%TDarkEnergyModel%ReadParams(Ini) + if (Ini%HasKey('AxionEffectiveFluid_a_c')) then + error stop 'AxionEffectiveFluid inputs changed to AxionEffectiveFluid_fde_zc and AxionEffectiveFluid_zc' + end if this%w_n = Ini%Read_Double('AxionEffectiveFluid_w_n') - this%om = Ini%Read_Double('AxionEffectiveFluid_om') - this%a_c = Ini%Read_Double('AxionEffectiveFluid_a_c') + this%fde_zc = Ini%Read_Double('AxionEffectiveFluid_fde_zc') + this%zc = Ini%Read_Double('AxionEffectiveFluid_zc') call Ini%Read('AxionEffectiveFluid_theta_i', this%theta_i) end subroutine TAxionEffectiveFluid_ReadParams @@ -184,10 +189,14 @@ subroutine TAxionEffectiveFluid_Init(this, State) select type(State) class is (CAMBdata) - this%is_cosmological_constant = this%om==0 + this%is_cosmological_constant = this%fde_zc==0 this%pow = 3*(1+this%w_n) - this%omL = State%Omega_de - this%om !Omega_de is total dark energy density today + this%a_c = 1/(1+this%zc) this%acpow = this%a_c**this%pow + !Omega in early de at z=0 + this%om = 2*this%fde_zc/(1-this%fde_zc)*& + (State%grho_no_de(this%a_c)/this%a_c**4/State%grhocrit + State%Omega_de)/(1 + 1/this%acpow) + this%omL = State%Omega_de - this%om !Omega_de is total dark energy density today this%num_perturb_equations = 2 if (this%w_n < 0.9999) then ! n <> infinity diff --git a/fortran/DarkEnergyQuintessence.f90 b/fortran/DarkEnergyQuintessence.f90 index 4c76b5eb..49b62f5d 100644 --- a/fortran/DarkEnergyQuintessence.f90 +++ b/fortran/DarkEnergyQuintessence.f90 @@ -267,7 +267,7 @@ subroutine TQuintessence_PerturbationEvolve(this, ayprime, w, w_ix, & ayprime(w_ix)= vq ayprime(w_ix+1) = - 2*adotoa*vq - k*z*phidot - k**2*clxq - a**2*clxq*this%Vofphi(phi,2) - end subroutine TQuintessence_PerturbationEvolve + end subroutine TQuintessence_PerturbationEvolve ! Early Quintessence example, axion potential from e.g. arXiv: 1908.06995 @@ -297,7 +297,6 @@ end function TEarlyQuintessence_VofPhi subroutine TEarlyQuintessence_Init(this, State) - use FileUtils use Powell class(TEarlyQuintessence), intent(inout) :: this class(TCAMBdata), intent(in), target :: State @@ -308,7 +307,6 @@ subroutine TEarlyQuintessence_Init(this, State) real(dl), parameter :: splZero = 0._dl real(dl) lastsign, da_osc, last_a, a_c real(dl) initial_phi, initial_phidot, a2 - Type(TTextFile) Fout real(dl), dimension(:), allocatable :: sampled_a, phi_a, phidot_a, fde integer npoints, tot_points, max_ix logical has_peak @@ -348,30 +346,32 @@ subroutine TEarlyQuintessence_Init(this, State) end do call Timer%WriteTime('Timing for fitting') end if - if (.True.) then - if (this%DebugLevel>0) call Timer%Start() - !Minimize in log f, log m - param_min(1) = log(0.001_dl) - param_min(2) = log(1d-57) - param_max(1) = log(1e5_dl) - param_max(2) = log(1d-50) - if (Minimize%BOBYQA(this, match_fde_zc, 2, 5, log_params,param_min, & - param_max, 0.8_dl,1e-4_dl,this%DebugLevel,2000)) then - if (Minimize%Last_bestfit > 1e-3) then - print *, 'ERROR finding good solution for fde, zc - best fit = ', Minimize%Last_bestfit - end if - this%f = exp(log_params(1)) - this%m = exp(log_params(2)) - if (this%DebugLevel>0) then - call this%calc_zc_fde(fzero, xzero) - print *, 'matched outputs Bobyqa zc, fde = ', fzero, xzero - end if - else - print *, 'ERROR finding solution for fde, zc' + if (this%DebugLevel>0) call Timer%Start() + !Minimize in log f, log m + param_min(1) = log(0.001_dl) + param_min(2) = log(1d-57) + param_max(1) = log(1e5_dl) + param_max(2) = log(1d-50) + if (Minimize%BOBYQA(this, match_fde_zc, 2, 5, log_params,param_min, & + param_max, 0.8_dl,1e-4_dl,this%DebugLevel,2000)) then + if (Minimize%Last_bestfit > 1e-3) then + global_error_flag = error_darkenergy + global_error_message= 'TEarlyQuintessence ERROR converging solution for fde, zc' + print *, 'last-bestfit= ', Minimize%Last_bestfit + return + end if + this%f = exp(log_params(1)) + this%m = exp(log_params(2)) + if (this%DebugLevel>0) then + call this%calc_zc_fde(fzero, xzero) + print *, 'matched outputs Bobyqa zc, fde = ', fzero, xzero end if - if (this%DebugLevel>0) call Timer%WriteTime('Timing for BOBYQA') + else + global_error_flag = error_darkenergy + global_error_message= 'TEarlyQuintessence ERROR finding solution for fde, zc' + return end if - + if (this%DebugLevel>0) call Timer%WriteTime('Timing for BOBYQA') end if this%dloga = (-this%log_astart)/(this%npoints-1) @@ -449,7 +449,6 @@ subroutine TEarlyQuintessence_Init(this, State) ind=1 afrom=this%log_astart - !call Fout%CreateOpenFile('z:\test.txt') do i=1, npoints-1 aend = this%log_astart + this%dloga*i ix = i+1 @@ -468,15 +467,12 @@ subroutine TEarlyQuintessence_Init(this, State) lastsign= y(2) end if - !Define fde is ratio of energy density to total assuming the neutrinos fully relativistic - !adotrad = sqrt((this%grhog+this%grhornomass+sum(this%grhormass(1:this%CP%Nu_mass_eigenstates)))/3) - + !Define fde as ratio of early dark energy density to total fde(ix) = 1/((this%state%grho_no_de(sampled_a(ix)) + this%frac_lambda0*this%State%grhov*a2**2) & /(a2*(0.5d0* phidot_a(ix)**2 + a2*this%Vofphi(y(1),0))) + 1) if (max_ix==0 .and. ix > 2 .and. fde(ix)< fde(ix-1)) then max_ix = ix-1 end if - !call Fout%Write(exp(aend), phi_a(ix), phidot_a(ix), fde(ix)) if (sampled_a(ix)*(exp(this%dloga)-1)*this%min_steps_per_osc > da_osc) then !Step size getting too big to sample oscillations well exit @@ -518,11 +514,8 @@ subroutine TEarlyQuintessence_Init(this, State) if (max_ix==0 .and. this%fde(ix)< this%fde(ix-1)) then max_ix = ix-1 end if - !call Fout%Write(aend, this%phi_a(ix), this%phidot_a(ix), this%fde(ix)) end do - !call Fout%Close() - call spline(this%sampled_a,this%phi_a,tot_points,splZero,splZero,this%ddphi_a) call spline(this%sampled_a,this%phidot_a,tot_points,splZero,splZero,this%ddphidot_a) call spline(this%sampled_a,this%fde,tot_points,splZero,splZero,this%ddfde) @@ -540,10 +533,12 @@ subroutine TEarlyQuintessence_Init(this, State) this%zc = 1/a_c-1 this%fde_zc = this%fdeAta(a_c) else - print *, 'NO PEAK FINAL' + if (this%DebugLevel>0) write(*,*) 'TEarlyQuintessence: NO PEAK ' this%zc = -1 end if - print *, 'zc, fde used', this%zc, this%fde_zc + if (this%DebugLevel>0) then + write(*,*) 'TEarlyQuintessence zc, fde used', this%zc, this%fde_zc + end if end subroutine TEarlyQuintessence_Init @@ -610,8 +605,10 @@ function match_fde_zc(this, x) this%f = exp(x(1)) this%m = exp(x(2)) call this%calc_zc_fde(zc, fde_zc) - match_fde_zc = (log(this%fde_zc)-log(fde_zc))**2 + (zc/this%zc-1)**2 - + + match_fde_zc = (log(this%fde_zc)-log(fde_zc))**2 + (log(zc)-log(this%zc))**2 + !print *, 'f, m, zc, fde_zc', this%f, this%m, zc, fde_zc, match_fde_zc + end function match_fde_zc subroutine calc_zc_fde(this, z_c, fde_zc) diff --git a/fortran/PowellConstrainedMinimize.f90 b/fortran/PowellConstrainedMinimize.f90 index e378f137..684acf99 100644 --- a/fortran/PowellConstrainedMinimize.f90 +++ b/fortran/PowellConstrainedMinimize.f90 @@ -2,7 +2,7 @@ !F90 translated version of plato.asu.edu/ftp/other_software/bobyqa.zip !Powell 2009 Method http://www.damtp.cam.ac.uk/user/na/NA_papers/NA2009_06.pdf - !AL Sept 2012: translated to F90, changed to input "funkk" class function argument + !AL Sept 2012/2020: translated to F90, changed to input "funkk" class function argument !to minimize rather than calfun subroutine as original MODULE Powell @@ -41,7 +41,7 @@ function BOBYQA (this, obj, funk, N,NPT,X,XL,XU,RHOBEG,RHOEND,IPRINT,MAXFUN) use iso_c_binding class(TBOBYQA) this class(*), target :: obj - real(dp), external :: funk !a class function f(obj,x) + real(dp), external :: funk !a class function f(obj,x) where obj is any class instance integer, intent(in)::n, npt, maxfun, iPrint real(dp), intent(in) :: RHOBEG, RHOEND logical BOBYQA diff --git a/fortran/inidriver.f90 b/fortran/inidriver.f90 index 51225cd7..bdb2c5de 100644 --- a/fortran/inidriver.f90 +++ b/fortran/inidriver.f90 @@ -13,6 +13,7 @@ program driver if (InputFile == '') error stop 'No parameter input file' call CAMB_CommandLineRun(InputFile) + deallocate(InputFile) ! Just so no memory leaks in valgrind end program driver