Use nonzero initial conformal time and radiation-dominated initial scale factor

Update cosmologies.jl

Author: hersle

docs/src/comparison.md

@@ -199,21 +199,21 @@ nothing # hide
 ```@example class
 a1 = (1 ./ (sol1["background"][:,"z"] .+ 1))
 a2 = sol2[M.g.a]
-χ1 = sol1["background"][:,"conf. time [Mpc]"][end].-sol1["background"][:,"conf. time [Mpc]"]
-χ2 = sol2[M.χ] / (h*SymBoltz.k0)
-plot_compare(a1, a2, χ1, χ2, "a", "χ"; tol = 2e-3)
+τ1 = sol1["background"][:,"conf. time [Mpc]"]
+τ2 = sol2[M.τ] / (h*SymBoltz.k0)
+plot_compare(a1, a2, τ1, τ2, "a", "χ"; tol = 2e-2)
 ```
 ### Hubble function
 ```@example class
 E1 = sol1["background"][:,"H [1/Mpc]"]./sol1["background"][end,"H [1/Mpc]"]
 E2 = sol2[M.g.E]
-plot_compare(a1, a2, E1, E2, "a", "E"; lgx=true, lgy=true, tol = 2e8)
+plot_compare(a1, a2, E1, E2, "a", "E"; lgx=true, lgy=true, tol = 8e6)
 ```
 ### Energy densities
 ```@example class
 ρ1 = map(s -> sol1["background"][:,"(.)rho_$s"], ["g", "ur", "cdm", "b", "fld", "ncdm[0]"])
 ρ2 = map(s -> sol2[s.ρ] * 8π/3*(h*SymBoltz.k0)^2, [M.γ, M.ν, M.c, M.b, M.X, M.h])
-plot_compare(a1, a2, ρ1, ρ2, "a", ["ργ", "ρb", "ρc", "ρX", "ρν", "ρh"]; lgx=true, lgy=true, tol = 2e15)
+plot_compare(a1, a2, ρ1, ρ2, "a", ["ργ", "ρb", "ρc", "ρX", "ρν", "ρh"]; lgx=true, lgy=true, tol = 6e12)
 ```
 ### Equations of state
 ```@example class
@@ -227,13 +227,13 @@ plot_compare(a1, a2, [wh1, wX1], [wh2, wX2], "a", ["wh", "wX"]; lgx=true, tol =
 ```@example class
 rs1 = sol1["background"][:,"comov.snd.hrz."]
 rs2 = sound_horizon(sol2) ./ (h*SymBoltz.k0)
-plot_compare(a1, a2, rs1, rs2, "a", "rₛ"; lgx = true, tol = 2e-2)
+plot_compare(a1, a2, rs1, rs2, "a", "rₛ"; lgx = true, tol = 1e-2)
 ```
 ### Luminosity distance
 ```@example class
 dL1 = sol1["background"][:,"lum. dist."]
 dL2 = SymBoltz.distance_luminosity(sol2) / SymBoltz.Mpc
-plot_compare(a1, a2, dL1, dL2, "a", "dL"; lgx=true, lgy=true, tol = 1e6)
+plot_compare(a1, a2, dL1, dL2, "a", "dL"; lgx=true, lgy=true, tol = 3e5)
 ```
 
 ## Thermodynamics
@@ -262,15 +262,15 @@ plot_compare(a1, a2, v1, v2, "a", "v"; lgx=true, lgy=false, tol = 2e-6)
 ```@example class
 Xe1 = reverse(sol1["thermodynamics"][:,"x_e"])
 Xe2 = sol2[M.b.Xe]
-plot_compare(a1, a2, Xe1, Xe2, "a", "Xe"; lgx=true, lgy=false, tol = 5e-4)
+plot_compare(a1, a2, Xe1, Xe2, "a", "Xe"; lgx=true, lgy=false, tol = 4e-4)
 ```
 ### Baryon temperature
 ```@example class
 Tb1 = reverse(sol1["thermodynamics"][:,"Tb [K]"])
 Tb2 = sol2[M.b.T]
 dTb1 = reverse(sol1["thermodynamics"][:,"dTb [K]"])
 dTb2 = sol2[M.b.DT] ./ -sol2[M.g.E] # convert my dT/dt̂ to CLASS' dT/dz = -1/H * dT/dt
-plot_compare(a1, a2, [Tb1, dTb1], [Tb2, dTb2], "a", ["Tb", "dTb"]; lgx=true, lgy=true, tol = 6e0)
+plot_compare(a1, a2, [Tb1, dTb1], [Tb2, dTb2], "a", ["Tb", "dTb"]; lgx=true, lgy=true, tol = 5e0)
 ```
 ### Baryon equation of state
 ```@example class
@@ -294,7 +294,7 @@ a1 = sol1["perturbations_k0_s"][:,"a"]
 a2 = sol2[1, M.g.a]
 Φ1, Ψ1 = sol1["perturbations_k0_s"][:,"phi"], sol1["perturbations_k0_s"][:,"psi"]
 Φ2, Ψ2 = sol2[1, M.g.Φ], sol2[1, M.g.Ψ]
-plot_compare(a1, a2, [Φ1, Ψ1], [Φ2, Ψ2], "a", ["Ψ", "Φ"]; lgx=true, tol = 2e-4)
+plot_compare(a1, a2, [Φ1, Ψ1], [Φ2, Ψ2], "a", ["Ψ", "Φ"]; lgx=true, tol = 2e-3)
 ```
 ### Energy overdensities
 ```@example class
@@ -306,25 +306,25 @@ plot_compare(a1, a2, δ1, δ2, "a", ["δb", "δc", "δγ", "δν", "δh"]; lgx=t
 ```@example class
 θ1 = map(s -> sol1["perturbations_k0_s"][:,"theta_$s"], ["b", "cdm", "g", "ur", "ncdm[0]"])
 θ2 = map(s -> sol2[1, s.θ] * (h*SymBoltz.k0), [M.b, M.c, M.γ, M.ν, M.h])
-plot_compare(a1, a2, θ1, θ2, "a", ["θb", "θc", "θγ", "θν", "θh"]; lgx=true, lgy=true, tol = 2e-2)
+plot_compare(a1, a2, θ1, θ2, "a", ["θb", "θc", "θγ", "θν", "θh"]; lgx=true, lgy=true, tol = 3e-2)
 ```
 ### Dark energy overdensity
 ```@example class
 δρX1 = sol1["perturbations_k0_s"][:,"delta_rho_fld"]
 δρX2 = sol2[1, M.X.δ*M.X.ρ] * 8π/3*(h*SymBoltz.k0)^2
-plot_compare(a1, a2, δρX1, δρX2, "a", "δρX"; lgx=true, lgy=true, tol = 2e-7)
+plot_compare(a1, a2, δρX1, δρX2, "a", "δρX"; lgx=true, lgy=true, tol = 4e-7)
 ```
 ### Dark energy momentum
 ```@example class
 pX1 = sol1["perturbations_k0_s"][:,"rho_plus_p_theta_fld"]
 pX2 = sol2[1, (M.X.ρ+M.X.P)*M.X.θ * 8π/3*(h*SymBoltz.k0)^3]
-plot_compare(a1, a2, pX1, pX2, "a", "pX"; lgx=true, lgy=true, tol = 7e-8)
+plot_compare(a1, a2, pX1, pX2, "a", "pX"; lgx=true, lgy=true, tol = 3e-7)
 ```
 ### Shear stresses
 ```@example class
 σ1 = [sol1["perturbations_k0_s"][:,"shear_g"], sol1["perturbations_k0_s"][:,"shear_ur"]]
 σ2 = [sol2[1, M.γ.σ], sol2[1, M.ν.F[2]/2]]
-plot_compare(a1, a2, σ1, σ2, "a", ["σγ", "σν"]; lgx=true, tol = 6e-4)
+plot_compare(a1, a2, σ1, σ2, "a", ["σγ", "σν"]; lgx=true, tol = 7e-4)
 ```
 ### Polarization
 ```@example class
@@ -354,7 +354,7 @@ function P_symboltz(k, pars)
 end
 k, P1 = P_class(pars)
 P2 = P_symboltz(k, pars)
-plot_compare(k, k, P1, P2, "k/Mpc⁻¹", "P/Mpc³"; lgx = true, lgy = true, tol = 4e1)
+plot_compare(k, k, P1, P2, "k/Mpc⁻¹", "P/Mpc³"; lgx = true, lgy = true, tol = 5e1)
 ```
 ```@example class
 using ForwardDiff, FiniteDiff

docs/src/observables.md

@@ -135,7 +135,7 @@ sol = solve(prob)
 zs = 0.0:1.0:10.0
 τs = SymBoltz.timeseries(sol, M.g.z, zs) # times at given redshifts
 dLs = SymBoltz.distance_luminosity(sol, τs) / SymBoltz.Gpc
-@assert dLs[begin] == 0.0 || zs[begin] != 0.0 # ensure bug does not reappear # hide
+@assert isapprox(dLs[begin], 0.0; atol = 1e-14) || zs[begin] != 0.0 # ensure bug does not reappear # hide
 plot(zs, dLs; marker=:dot, xlabel="z", ylabel="dL / Gpc", label=nothing)
 ```
 

src/models/cosmologies.jl

@@ -59,7 +59,8 @@ function ΛCDM(;
     end
     defs = Dict(
         C => 1//2,
-        τ0 => NaN
+        τ0 => NaN,
+        g.a => √(sum(s.Ω₀ for s in [ν, h, γ] if have(s))) * τ # default initial scale factor
     )
     ics = [
         g.Ψ ~ 20C / (15 + 4fν) # Φ found from solving initialization system
@@ -152,7 +153,8 @@ function RMΛ(;
     ]
     defs = Dict(
         g.Ψ => 20 // 15,
-        τ0 => NaN
+        τ0 => NaN,
+        g.a => √(r.Ω₀) * τ # default initial scale factor
     )
     append!(eqs, Ω₀_eqs(G, species)) # parameter equations
     connections = System(eqs, τ, vars, [pars; k]; defaults = defs, name)

src/models/metric.jl

@@ -20,9 +20,6 @@ function metric(; name = :g, kwargs...)
     vars = a, ℰ, E, ℋ, H, Ψ, Φ, Ψ̇, Φ̇, z, ż = GlobalScope.(vars)
 
     pars = h, = GlobalScope.(@parameters h [description = "Dimensionless Hubble parameter today (H₀/(100km/s/Mpc))"])
-    defaults = Dict(
-        a => 1e-8 # default initial scale factor # TODO: move elsewhere?
-    )
     description = "Spacetime FLRW metric in Newtonian gauge"
     return System([
         z ~ 1/a - 1
@@ -31,5 +28,5 @@ function metric(; name = :g, kwargs...)
         E ~ ℰ / a # E = H/H₀
         ℋ ~ ℰ * H100 * h
         H ~ E * H100 * h
-    ], τ, vars, pars; defaults, name, description, kwargs...)
+    ], τ, vars, pars; name, description, kwargs...)
 end

src/observables/angular.jl

@@ -166,12 +166,13 @@ function spectrum_cmb(modes::AbstractVector, prob::CosmologyProblem, ls::Abstrac
     τ0 = getsym(sol, prob.M.τ0)(sol)
     ks_coarse = kτ0s_coarse ./ τ0
     τs = sol.bg.t # by default, use background (thermodynamics) time points for line of sight integration
+    τi = τs[begin]
     if Nlos != 0 # instead choose Nlos time points τ = τ(u) corresponding to uniformly spaced u
         τmin, τmax = extrema(τs)
         umin, umax = u(τmin), u(τmax)
         us = range(umin, umax, length = Nlos)
         τs = u⁻¹.(us)
-        τs[begin] = 0.0
+        τs[begin] = τi
         τs[end] = τ0
     end
 

src/solve.jl

@@ -104,7 +104,7 @@ end
 """
     CosmologyProblem(
         M::System, pars::Dict, shoot_pars = Dict(), shoot_conditions = [];
-        ivspan = (0.0, 100.0), bg = true, pt = true, spline = true, debug = false, fully_determined = true, jac = false, kwargs...
+        ivspan = (1e-6, 100.0), bg = true, pt = true, spline = true, debug = false, fully_determined = true, jac = false, kwargs...
     )
 
 Create a numerical cosmological problem from the model `M` with parameters `pars`.
@@ -116,7 +116,7 @@ If `spline` is a `Vector`, it rather decides which (unknown and observed) variab
 """
 function CosmologyProblem(
     M::System, pars::Dict, shoot_pars = Dict(), shoot_conditions = [];
-    ivspan = (0.0, 100.0), bg = true, pt = true, spline = true, debug = false, fully_determined = true, jac = false, kwargs... # TODO: restore jac = true, not working with massive neutrinos?
+    ivspan = (1e-6, 100.0), bg = true, pt = true, spline = true, debug = false, fully_determined = true, jac = false, kwargs... # TODO: restore jac = true, not working with massive neutrinos?
 )
     pars = merge(pars, shoot_pars) # save full dictionary for constructor
     parsk = merge(pars, Dict(M.k => NaN)) # k is unused, but must be set

test/runtests.jl

@@ -141,15 +141,16 @@ end
 end
 
 @testset "Initial conditions" begin
+    τini = prob.bg.tspan[1]
     ks = [1e-1, 1e0] / u"Mpc"
     sol = solve(prob, ks)
 
     # Check that Fₗ(0) ∝ kˡ
-    Fls = sol([M.γ.F0; collect(M.γ.F)], 0.0, ks)
+    Fls = sol([M.γ.F0; collect(M.γ.F)], τini, ks)
     @test all(Fls[:,1] ./ Fls[:,2] .≈ map(l -> (ks[1]/ks[2])^l, 0:size(Fls)[1]-1))
 
     # Check initial ratio of metric potentials
-    @test all(isapprox.(sol(M.g.Φ / M.g.Ψ, 0.0, ks), sol((1+2/5*M.fν), 0.0); atol = 1e-5))
+    @test all(isapprox.(sol(M.g.Φ / M.g.Ψ, τini, ks), sol((1+2/5*M.fν), τini); atol = 1e-4))
 end
 
 @testset "Automatic background/thermodynamics splining" begin
@@ -336,8 +337,9 @@ end
 
     # custom output_func for e.g. source function
     getS = SymBoltz.getsym(prob.pt, M.ST0)
+    τi = bgsol.t[begin]
     τ0 = bgsol.t[end]
-    τs = range(0.0, τ0, length = 768)
+    τs = range(τi, τ0, length = 768)
     ks = range(1.0, 1000.0, length = 1000)
     ptsol = solvept(prob.pt, bgsol, ks, prob.var2spl; saveat = τs, output_func = (ptsol, _) -> (getS(ptsol), false))
     Ss = stack(ptsol.u)
@@ -348,15 +350,16 @@ end
 # TODO: define closure based on this?
 @testset "Type stability" begin
     getτ0 = SymBoltz.getsym(prob.bg, M.τ0)
-    sol = solve(prob; save_everystep = false, save_start = false, save_end = true)
+    sol = solve(prob; save_everystep = false, save_start = true, save_end = true)
+    τi = sol.bg.t[begin]
     τ0 = @inferred getτ0(sol.bg)
-    @test sol.bg.t[begin] == sol.bg.t[end] == τ0
+    @test sol.bg.t[end] == τ0
 
     ls = 20:20:500
     kτ0s_coarse, kτ0s_fine = SymBoltz.cmb_kτ0s(ls[begin], ls[end])
     ks_coarse, ks_fine = kτ0s_coarse / τ0, kτ0s_fine / τ0
     ks_fine = clamp.(ks_fine, ks_coarse[begin], ks_coarse[end]) # numerics can change endpoint slightly
-    τs = range(0.0, τ0, length = 768)
+    τs = range(τi, τ0, length = 768)
     #sol = solve(prob, ks_coarse)
     #Ss = sol(ks_fine, τs, M.ST0)
     sol = solve(prob, ks_coarse; ptopts = (alg = SymBoltz.KenCarp4(), reltol = 1e-8, abstol = 1e-8, saveat = τs,))