Solving models

First create a symbolic cosmological model:

using SymBoltz
M = SymBoltz.ΛCDM()

Creating the problem

Once the symbolic cosmological model M has been constructed, it can be turned into a numerical problem: For example:

pars = SymBoltz.parameters_Planck18(M)
prob = CosmologyProblem(M, pars)

Cosmology problem for model ΛCDM
Stages:
  Background: 5 unknowns, 0 splines
  Perturbations: 61 unknowns, 5 splines
Parameters & initial conditions:
  c₊Ω₀ = 0.26447041034523616
  b₊Ω₀ = 0.04936780993111075
  γ₊T₀ = 2.7255
  b₊YHe = 0.2454
  h₊m_eV = 0.06
  h = 0.6736
  I₊ns = 0.965
  ν₊Neff = 2.99
  I₊ln_As1e10 = 3.0440461338325417

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...
)

Updating the parameters

Constructing a CosmologyProblem is an expensive operation that compiles all the symbolics down to numerics. It is not necessary to repeat this just to update parameter values. To do so, use the function parameter_updater that returns a function that quickly creates new problems with updated parameter values:

probmaker = parameter_updater(prob, [M.g.h, M.c.Ω₀]) # fast factory function
prob = probmaker([0.70, 0.27]) # create updated problem

Cosmology problem for model ΛCDM
Stages:
  Background: 5 unknowns, 0 splines
  Perturbations: 61 unknowns, 5 splines
Parameters & initial conditions:
  c₊Ω₀ = 0.27
  b₊Ω₀ = 0.04936780993111075
  γ₊T₀ = 2.7255
  b₊YHe = 0.2454
  h₊m_eV = 0.06
  h = 0.7
  I₊ns = 0.965
  ν₊Neff = 2.99
  I₊ln_As1e10 = 3.0440461338325417

parameter_updater(prob::CosmologyProblem, idxs; kwargs...)

Solving the problem

The (updated) problem can now be solved for some wavenumbers:

using Unitful, UnitfulAstro
ks = 10 .^ range(-5, +1, length=100) / u"Mpc"
sol = solve(prob, ks)

Cosmology solution for model ΛCDM
Stages:
  Background: return code Terminated; solved with Rodas4P; 1311 points
  Perturbations: return codes Success; solved with KenCarp4; 81-826 points; x100 k ∈ [0.042827494000000015, 42827.49400000001] H₀/c (interpolation in-between)

CommonSolve.solve(args...; kwargs...)

CommonSolve.solve(prob::ProblemType,alg::SolverType; kwargs...)::SolutionType

solve(args...; kwargs...) = solve!(init(args...; kwargs...))

function solve(
    prob::CosmologyProblem, ks::Union{Nothing, AbstractArray} = nothing;
    bgopts = (alg = Rodas4P(), reltol = 1e-9, abstol = 1e-9),
    ptopts = (alg = KenCarp4(), reltol = 1e-8, abstol = 1e-8),
    shootopts = (alg = NewtonRaphson(), abstol = 1e-5),
    thread = true, verbose = false, kwargs...
)

Accessing the solution

The returned solution sol can be conveniently accessed to obtain any variable y of the model M:

• sol(y, τ) returns the background variable(s) $y(τ)$ as a function of conformal time(s) $τ$. It interpolates between time points using the ODE solver's custom-tailored interpolator.
• sol(y, τ, k) returns the perturbation variable(s) $y(τ,k)$ as a function of the wavenumber(s) $k$ and conformal time(s) $τ$. It also interpolates linearly between the logarithms of the wavenumbers passed to solve.

Note that y can be any symbolic variables in the model M, and even expressions thereof. Unknown variables are part of the state vector integrated by the ODE solver, and are returned directly from its solution. Observed variables or expressions are functions of the unknowns, and are automatically calculated from the equations that define them in the symbolic model. For example:

τs = sol[M.τ] # get time points used in the background solution
ks = [1e-3, 1e-2, 1e-1, 1e0] / u"Mpc" # wavenumbers
as = sol(M.g.a, τs) # scale factors
Ωms = sol((M.b.ρ + M.c.ρ) / M.G.ρ, τs) # matter-to-total density ratios
κs = sol(M.b.κ, τs) # optical depths
Φs = sol(M.g.Φ, τs, ks) # metric potentials
Φs_over_Ψs = sol(M.g.Φ / M.g.Ψ, τs, ks) # ratio between metric potentials

Plotting the solution

SymBoltz.jl includes plot recipes for easily visualizing the solution. It works similarly to the solution accessing: call plot(sol, [wavenumber(s),] x_expr, y_expr) to plot y_expr as a function of x_expr. For example, to plot some of the same quantities that we obtained above:

using Plots
p1 = plot(sol, log10(M.g.a), (M.b.ρ + M.c.ρ) / M.G.ρ)
p2 = plot(sol, log10(M.g.a), log10(abs(M.b.κ)))
p3 = plot(sol, log10(M.g.a), M.g.Φ / M.g.Ψ, ks[1:3]) # exclude last k, where Φ and Ψ cross 0
plot(p1, p2, p3, layout=(3, 1), size=(600, 800))

More examples are shown on the models page.

Solve background and perturbations directly

For lower-level control, you can solve the background and perturbations separately:

solvebg(bgprob::ODEProblem; alg = Rodas4P(), reltol = 1e-9, abstol = 1e-9, verbose = false, kwargs...)

solvept(ptprob::ODEProblem, bgsol::ODESolution, ks::AbstractArray, var2spl::Dict; alg = KenCarp4(), reltol = 1e-8, abstol = 1e-8, output_func = (sol, i) -> (sol, false), thread = true, verbose = false, kwargs...)

Choice of ODE solver

In principle, models can be solved with any DifferentialEquations.jl ODE solver. But most cosmological models have very stiff Einstein-Boltzmann equations that can only be solved by implicit solvers, while explicit solvers usually fail. For the standard ΛCDM model, some good solvers are:

• Rodas4P: Slow for large systems. Very accurate. Handles extreme stiffness. Default background solver.
• Rodas5P: Slow for large systems. Very accurate. Handles severe stiffness.
• KenCarp4 (and KenCarp47): Fast. Handles medium stiffness. Default perturbation solver.
• Kvaerno5: Behaves similar to KenCarp4. Slightly more accurate. Slightly slower.
• TRBDF2: Very fast. Decent accuracy. Handles severe stiffness.

See also the solver benchmarks.