Rodas5 instability with current massive neutrinos

[jmsull]

Following up from mass_nu PR: """ Regarding the instability you encountered -- let's turn that into an issue, I imagine it's easy to make an minimal demonstration. My current set of questions is

1. Does this happen if you remove the effect of massive neutrinos on the metric? (i.e. they're just free streaming the entire time)
2. Does this happen if you truncate the massive neutrinos to just the monopole / dipole?
3. Does this happen when you change the q grid?

However, we can deal with it later.

Originally posted by @xzackli in https://github.com/xzackli/Bolt.jl/issues/39#issuecomment-795041642 """

[jmsull]

For 3., I did try adding 7x as many q-points and it didn't help, but this is when I was using a regular grid instead of quadrature points. Changing from the regular grid to quadrature points did not seem to help either.

I have not tried 1, or 2 but will check on those.

[jmsull]

For 1., removing the massive neutrino contributions to the metric does not help. I still get the dt<=dtmin error.

For 2. there is still an instability, but the error is slightly different with no dtmin warning, just: "Warning: Instability detected. Aborting". This is when I truncate at the quadrupole (\ellmax=2). Instead truncating at the dipole (\ellmax=1), I get the old "dt<=dtmin" error.

For 3. Changing the number of quadrature q points (down to 5 or up to 150) does not seem to help either (with or without the metric fluctuations).

[xzackli]

Just briefly poking this problem, haven't looked at equations yet. Actually, I can make it successfully solve if I comment out the impact of relativistic neutrinos on the metric.

# metric perturbations (00 and ij FRW Einstein eqns)
    Ψ = -Φ - 12H₀² / k^2 / a^2 * (Ω_r * Θ[2]
                                #   + Ω_ν * 𝒩[2] #add rel quadrupole
                                  + σℳ / bg.ρ_crit/ norm𝒩′) #add mnu integrated quadrupole

    # println("x= ",x, " so a = ", exp(x))
    # println("Size of terms in i neq j eqn. Ω_ν: ", Ω_ν * 𝒩[2], " and σℳ ", σℳ / bg.ρ_crit / norm𝒩′)

    Φ′ = Ψ - k^2 / (3ℋₓ^2) * Φ + H₀² / (2ℋₓ^2) * (
        Ω_m * a^(-1) * δ + Ω_b * a^(-1) * δ_b + 4Ω_r * a^(-2) * Θ[0]
        # + 4Ω_ν * a^(-2) * 𝒩[0] #add rel monopole on this line
        + a^(-2) * ρℳ  / bg.ρ_crit / norm𝒩′) #add mnu integrated monopole

Here's my little test script

using Bolt
using BenchmarkTools
using OrdinaryDiffEq

par = CosmoParams()
bg = Background(par)
ih = IonizationHistory(Peebles(), par, bg)

k_grid = quadratic_k(0.1bg.H₀, 1000bg.H₀, 100)
k_i = 10
hierarchy = Hierarchy(BasicNewtonian(), par, bg, ih, k_grid[k_i])

# xᵢ = log(1e-7)
xᵢ = (hierarchy.bg.x_grid)[10]
u₀ = Bolt.initial_conditions(xᵢ, hierarchy)
prob = ODEProblem{true}(Bolt.hierarchy!, u₀, (xᵢ, 0.0), hierarchy)
@time sol = solve(prob, Rodas5(), reltol=1e-4)

[jmsull]

This is surprising - I reproduce as you say when I try this example. Similarly, in the Einstein equations when I leave both species of neutrinos in but comment out the photons it also successfully solves...

[xzackli]

The minimal change that makes it stable is removing the relativistic contribution to Ψ,

Ψ = -Φ - 12H₀² / k^2 / a^2 * (Ω_r * Θ[2]
                                #   + Ω_ν * 𝒩[2] #add rel quadrupole
                                  + σℳ / bg.ρ_crit/ norm𝒩′) #add mnu integrated quadrupole

[xzackli]

By the way, I might be misunderstanding, but is there a difference in f_nu convention?

https://github.com/xzackli/Bolt.jl/blob/80255b4b3df2b5338f95e60452757d3f7f27081c/src/perturbations.jl#L212

[jmsull]

No sorry that's a mistake on my part - it should be 5/(2f\nu) not 5f\nu/2. Thanks for catching that.

[xzackli]

All good, I didn't catch it when I merged! I think the time has come for me to write some unit tests, in which I make up some nice varied numbers for inputs, and then work out the expected output using Mathematica or something. That'll help with general debugging and development efforts, especially if we do stuff like change notation.

Here's my list of suspects:

1. We could have made some user errors and are solving an incorrect ODE or initial condition. Also maybe units?
2. The solver is encountering instability at extremely early times when the fluid is tightly coupled. Maybe we really do need TCA, and a modern stiff solver isn't enough in the presence of neutrinos?
3. The current hierarchy is based on Callin06, which honestly hasn't been extensively tested. It's very easy to make typos (I feel like I encounter errors in papers for every project I work on).

The plans of attack are therefore to (1) write detailed tests, (2) try the TCA in Callin, and (3) implement the Newtonian gauge equations from MB instead. These are all things I need to do at some point anyway, in pursuit of correctness. I won't be able to start on this until later this week though.

[jmsull]

The instability error seems to have gone away in #44, likely the solver was struggling with the incorrect implementation of the massive neutrino hierarchy