Remove h from CosmologySolution (get it as any other parameter)

Author: hersle

src/constants.jl

@@ -44,12 +44,7 @@ const EHet2p2s = EHet2p1s - EHet2s1s
 
 δkron(i, j) = (i == j ? 1 : 0) # Kronecker delta
 
-function k_dimensionless(k::Number, h)
-    if unit(k) == NoUnits
-        return k
-    else
-        H₀ = h * H100 # s⁻¹
-        k0 = (H₀ / c) / u"m"
-        return NoUnits(k / k0)
-    end
-end
+k_dimensionless(k::Number, h) = k
+k_dimensionless(k::Quantity, h) = NoUnits(k / (h*H100 / c / u"m"))
+k_dimensionless(k::Number, bgsol::ODESolution) = k
+k_dimensionless(k::Quantity, bgsol::ODESolution) = k_dimensionless(k, getsym(bgsol, :h)(bgsol))

src/models/inflation.jl

@@ -11,7 +11,7 @@ function harrison_zeldovich(g; kwargs...)
         P ~ 2*Num(π)^2 / k^3 * As * (k/kpivot)^(ns-1)
     ]
     defaults = [
-        kpivot => 0.05 / (H100*g.h/c) / Mpc # k = 0.05/Mpc ≠ 0.05/(Mpc/h)
+        kpivot => 0.05 / Mpc / (H100/c) / g.h # k = 0.05/Mpc ≠ 0.05/(Mpc/h)
         ln_As1e10 => NaN # keep uninitialized if not needed
         ns => NaN
     ]

src/observables/fourier.jl

@@ -3,20 +3,21 @@ using MatterPower
 using TwoFAST
 
 """
-    spectrum_primordial(k, h, As, ns=1.0; kp = k_dimensionless(0.05 / u"Mpc", h))
+    spectrum_primordial(k, h, As, ns=1.0; kp = 0.05/u"Mpc")
 
 Compute the primordial power spectrum
 ```math
 P₀(k) = 2π² Aₛ (k/kₚ)^{nₛ-1} / k³
 ```
 with spectral amplitude `As`, spectral index `ns` and pivot scale wavenumber `kp` at the wavenumber(s) `k`.
 """
-function spectrum_primordial(k, h, As, ns=1.0; kp = k_dimensionless(0.05 / u"Mpc", h))
-    P = @. 2*π^2 / k^3 * As
+function spectrum_primordial(k, h, As, ns=1.0; kp = 0.05/u"Mpc")
+    P = 2*π^2 * As ./ k.^3
 
-    # compute k / kpivot with equal wavenumber units
+    # ensure both k and kp to dimensionless wavenumbers k/(H₀/c) before taking ratio
     k = k_dimensionless.(k, h)
-    @. P *= (k/kp)^(ns-1)
+    kp = k_dimensionless(kp, h)
+    P .*= (k./kp).^(ns-1)
 
     return P
 end
@@ -121,7 +122,7 @@ function variance_matter(sol::CosmologySolution, R)
     P = spectrum_matter(sol, k)
     lgPspl = spline(log.(P), log.(k))
     Pf(k) = exp(lgPspl(log(k)))
-    R = 1 / k_dimensionless(1 / R, sol[M.g.h]) # make dimensionless
+    R = 1 / k_dimensionless(1 / R, sol.bg) # make dimensionless
     return MatterPower.sigma2(Pf, R)
 end
 """

src/solve.jl

@@ -31,12 +31,11 @@ struct CosmologyProblem{Tbg <: ODEProblem, Tpt <: Union{ODEProblem, Nothing}, Tb
     bgspline::Tbgspline
 end
 
-struct CosmologySolution{Tbg <: ODESolution, Tpts <: Union{Nothing, EnsembleSolution}, Tks <: Union{Nothing, AbstractVector}, Th <: Number}
+struct CosmologySolution{Tbg <: ODESolution, Tpts <: Union{Nothing, EnsembleSolution}, Tks <: Union{Nothing, AbstractVector}}
     prob::CosmologyProblem # problem which is solved
     bg::Tbg # background solution
     ks::Tks # perturbation wavenumbers
     pts::Tpts # perturbation solutions
-    h::Th # reduced Hubble parameter h = H/(100 km/s/Mpc)
 end
 
 algname(alg) = string(nameof(typeof(alg)))
@@ -293,16 +292,15 @@ function solve(
         bgsol = solvebg(prob.bg; verbose, bgopts..., kwargs...)
     end
 
-    h = getsym(bgsol, :h)(bgsol) # TODO: remove symbolic getter
     if isnothing(ks) || isempty(ks)
         ks = nothing
         ptsol = nothing
     else
-        ks = k_dimensionless.(ks, h)
+        ks = k_dimensionless.(ks, Ref(bgsol))
         ptsol = solvept(prob.pt, bgsol, ks, prob.bgspline; thread, verbose, ptopts..., kwargs...)
     end
 
-    return CosmologySolution(prob, bgsol, ks, ptsol, h)
+    return CosmologySolution(prob, bgsol, ks, ptsol)
 end
 function solve(prob::CosmologyProblem, k::Number; kwargs...)
     return solve(prob, [k]; kwargs...)
@@ -570,7 +568,7 @@ function getsym(provider::Union{CosmologyProblem, CosmologySolution}, p)
 end
 
 function neighboring_modes_indices(sol::CosmologySolution, k)
-    k = k_dimensionless.(k, sol.h)
+    k = k_dimensionless.(k, Ref(sol.bg))
     if k == sol.ks[begin] # k == kmin
         i1 = i2 = 1
     elseif k == sol.ks[end] # k == kmax
@@ -588,7 +586,7 @@ function (sol::CosmologySolution)(out::AbstractArray, is::AbstractArray, ts::Abs
     if isnothing(sol.ks) || isempty(sol.ks)
         throw(error("No perturbations solved for. Pass ks to solve()."))
     end
-    ks = k_dimensionless.(ks, sol.h)
+    ks = k_dimensionless.(ks, Ref(sol.bg))
     kmin, kmax = extrema(sol.ks)
     minimum(ks) >= kmin || throw("Requested wavenumber k = $(minimum(ks)) is below the minimum solved wavenumber $kmin")
     maximum(ks) <= kmax || throw("Requested wavenumber k = $(maximum(ks)) is above the maximum solved wavenumber $kmax")

test/runtests.jl

@@ -199,11 +199,23 @@ end
     @test_throws "not an unknown" SymBoltz.mtkcompile_spline(M, [M.g.H])
 end
 
-@testset "Primordial power spectrum pivot scale" begin
+@testset "Wavenumber units and primordial power spectrum pivot scale" begin
     h = pars[M.g.h]
     k = 0.05 / u"Mpc" # ≠ 0.05/(Mpc/h)
     sol = solve(prob, k)
-    @test sol[M.I.kpivot] ≈ SymBoltz.k_dimensionless(k, h)
+    @test sol[M.I.kpivot] ≈ SymBoltz.k_dimensionless(k, sol.bg)
+
+    ks = 1.0:100.0
+    sol = solve(prob, ks)
+    P1 = sol(M.I.P, sol.bg.t[begin], ks)
+    P2 = spectrum_primordial(ks, sol)
+    @test all(isapprox.(P1, P2))
+
+    ks = collect(1.0:100.0) / u"Mpc"
+    sol = solve(prob, ks)
+    P1 = sol(M.I.P, sol.bg.t[begin], ks)
+    P2 = spectrum_primordial(SymBoltz.k_dimensionless.(ks, h), sol)
+    @test all(isapprox.(P1, P2))
 end
 
 @testset "Time and optical depth today" begin