Unitful tolerances

[devmotion]

The following code fails for the in-place method but works for the out-of-place version:

using OrdinaryDiffEq, Unitful

f_oop(y,p,t) = 0.5*y / 3.0u"s"
f_ip(dy,y,p,t) = (dy .= 0.5.*y ./ 3.0u"s")

# works
function test_oop()
    u0 = 1.0u"N"
    prob = ODEProblem(f_oop, u0, (0.0u"s", 1.0u"s"))

    solve(prob, Tsit5(); abstol = 1e-6, reltol = 1e-3)
end

# does not work
function test_ip()
    u0 = [1.0u"N"]
    prob = ODEProblem(f_ip, u0, (0.0u"s", 1.0u"s"))

    solve(prob, Tsit5(); abstol = 1e-6, reltol = 1e-3)
end

The inplace method fails in https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/initdt.jl#L11.

Actually I tend to think that the out-of-place method should fail as well since abstol is specified without units (alternatively, one could try to convert them to the correct units during setup of the integrator?). The out-of-place method does not fail in https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/initdt.jl#L130 at the moment since internalnorm does not preserve units:

julia> using DiffEqBase, Unitful
julia> DiffEqBase.ODE_DEFAULT_NORM(1u"N", 2.0u"s")
1

However, I think it should respect the units (as least in this case) since in my opinion the relative tolerance should be unitless but the units of the product of the norm and the relative tolerance has to have the same units as abstol (and additionally abs(1u"N") == 1u"N").

Since in my opinion the relative tolerance should be unitless, additionally I think that the logic in https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/solve.jl#L124-L134 is not correct since it automatically creates a unitful relative tolerance if none is specified (however, before changing that one would have to change internalnorm first).

[ChrisRackauckas]

Yes, abstol should respect units, but I didn't handle this for now because it ended up being really difficult so I got something mostly working where no units were specified in the tolerances, and that seemed good enough. It would be a good idea to spend some time getting this all correct. The real difficulty comes with heterogeneous units. I am just not sure how to handle this correctly.

[wnoise]

The current code for initializing abstol seems to respect units, even for heterogeneous arrays. The problem as mentioned above is is that reltol does not:

https://github.com/SciML/OrdinaryDiffEq.jl/blob/5fef4209b30b1461feac162b29a2e08eeb6920b5/src/solve.jl#L169-L177

Changing the norm for Unitful quantities to preserve the units (just dropping the call to value in the method), and reltol to something like reltol_internal = real.(DiffEqBase.value.(oneunit.(u)).*1//10^3) now has a snag on calculate_residuals, which returns a Vector{Any}[] all of whose elements are Unitful quantities with the no units. I don't understand why the type is this instead of Vector{Quantity{Float64, D, U} where {D, U}}, though presumably it has something to do with what FastBroadcast.@.. is doing to the returns from the per-element calculate_residuals. Replacing @.. with @. gets the even more specific type Vector{Quantity{Float64, NoDims, Unitful.FreeUnits{(), NoDims, nothing}}}, which proceeds to actually work, and finish solving.

But of course, there is no overload that handles interpolation for these arrays...