Remove allocations from `reinit!` with `reset_dt = true`

[jClugstor]

Is your feature request related to a problem? Please describe.
If you reinit! an integrator with reset_dt = true it always allocates. I would like to try to see if there's a way to calculate the initial dt value without allocating.

using OrdinaryDiffEq

# Simple exponential decay: du/dt = -u (in-place)
f!(du, u, _, _) = (du[1] = -u[1]; nothing)

u0 = [1.0]
tspan = (0.0, 1.0)
prob = ODEProblem(f!, u0, tspan)
integrator = init(prob, Tsit5())

# Warm up
reinit!(integrator)
solve!(integrator)

# reinit! allocates even on a plain ODE with no DAE initialization
println("Allocations from reinit!:")
allocs = @allocated reinit!(integrator)
@profview_allocs reinit!(integrator) sample_rate = 1.0

[jClugstor]

@oscardssmith any ideas on the best way to go about adding a better caching interface, specifically for the pieces needed for initdt?

It seems like every algorithm needs the same things for the dt initialization? In that case the cache for initdt should go in the integrators, and not the algorithm caches?

Claude suggested creating an InitDtCache which would hold all of the space needed for initdt, and putting that as a field to the integrator, but I'm not sure how we would feel about adding another field just for this.

[oscardssmith]

The integrator already has a cache of u values. You should just be able to use that...

[jClugstor]

Are you talking about u, uprev, uprev2?

[jClugstor]

In _ode_initdt_iip there are three uType arrays, one of which is from get_tmp_cache, the other two are newly allocated in _ode_initdt_iip. Then there are two rateType arrays, one sometimes uses integrator.fsallast, and one is always allocated with zero.

I could use u, uprev, uprev2, du, duprev but that would change those values, which I don't think is a good idea?

[oscardssmith]

IIUC, get_tmp_cache can get mutliple cache vectors.

[jClugstor]

Yes, but what you get depends on the type of the cache, from OrdinaryDiffEqCore, integrator_interface.jl:

#TODO: Bigger caches for most algorithms
@inline function SciMLBase.get_tmp_cache(integrator::ODEIntegrator)
    return get_tmp_cache(integrator::ODEIntegrator, integrator.alg, integrator.cache)
end

# the ordering of the cache arrays is important!!!
@inline function SciMLBase.get_tmp_cache(
        integrator, alg::OrdinaryDiffEqAlgorithm,
        cache::OrdinaryDiffEqConstantCache
    )
    return nothing
end
@inline function SciMLBase.get_tmp_cache(
        integrator, alg::OrdinaryDiffEqAlgorithm,
        cache::OrdinaryDiffEqMutableCache
    )
    return (cache.tmp,)
end
@inline function SciMLBase.get_tmp_cache(
        integrator,
        alg::OrdinaryDiffEqNewtonAdaptiveAlgorithm,
        cache::OrdinaryDiffEqMutableCache
    )
    return (cache.nlsolver.tmp, cache.atmp)
end

and none of the get_tmp_cache currently returns enough, in fact for constant caches it returns nothing.

[oscardssmith]

for constant caches you will just do the allocations (constant cache is for out of place where we assume that you can't mutate the us.

[jClugstor]

Right, my main point was that get_tmp_cache doesn't have enough space for everything that _ode_initdt_iip needs. At worst it returns cache.tmp which is just a single uType. At best it returns two uTypes.

I guess ideally we would want to be able to call get_tmp_cache and get all of the space needed to do initdt, but in order to do that I would need to add more fields to every single cache type. Or another option would be to make it so that every cache.tmp is a type that holds a few copies of uType, uNoUnitsType (forgot to mention that we need a few of these), and rateType.

[oscardssmith]

yeah this is fairly annoying. IMO it would be really nice if all solvers stored their vectors in compatible ways to enable easier reuse, but we aren't in that situation currently.