Use setuppt/ptprobgen in adaptive source grid function

Author: hersle

src/observables/fourier.jl

@@ -262,11 +262,12 @@ The options `bgopts` and `ptopts` are passed to the background and perturbation
 """
 function source_grid_adaptive(prob::CosmologyProblem, Ss::AbstractVector, τs, ks; bgopts = (), ptopts = (), kwargs...)
     bgsol = solvebg(prob.bg; bgopts...)
+    ptprob0, ptprobgen = setuppt(prob.pt, bgsol, prob.bgspline)
 
     getSs = map(S -> getsym(prob.pt, S), Ss)
     function Sk(k)
-        ptsols = solvept(prob.pt, bgsol, [k], prob.bgspline; saveat = τs, thread = false, ptopts...)
-        ptsol = only(ptsols)
+        ptprob = ptprobgen(ptprob0, k)
+        ptsol = solvept(ptprob; saveat = τs, ptopts...)
         S = transpose(stack(map(getS -> getS(ptsol), getSs)))
         S[:, end] .= 0.0 # avoid NaN/Infs from 1/χ today; fine in LOS integration, since always weighted by 0-valued Bessel function # TODO: avoid
         return S

src/solve.jl

@@ -427,6 +427,27 @@ function solvept(ptprob::ODEProblem, bgsol::ODESolution, ks::AbstractArray, bgsp
     verbose && println()
     return ptsols
 end
+"""
+    solvept(ptprob::ODEProblem; alg = DEFAULT_PTALG, reltol = 1e-8, abstol = 1e-8, kwargs...)
+
+Solve the perturbation problem `ptprob` and return the solution.
+Its wavenumber and background spline must already be initialized, for example with `setuppt`.
+
+# Examples
+
+```julia
+# ...
+prob = CosmologyProblem(M, pars)
+bgsol = solvebg(prob.bg)
+ptprob0, ptprobgen = SymBoltz.setuppt(prob.pt, bgsol, prob.bgspline)
+k = 1.0
+ptprob = ptprobgen(ptprob0, k)
+ptsol = solvept(ptprob)
+```
+"""
+function solvept(ptprob::ODEProblem; alg = DEFAULT_PTALG, reltol = 1e-8, abstol = 1e-8, kwargs...)
+    return solve(ptprob, alg; reltol, abstol, kwargs...)
+end
 
 function time_today(prob::CosmologyProblem)
     getτ0 = SymBoltz.getsym(prob.bg, :τ0)