2ndOrderODEProblem + Unitful, Type issues

[louisponet]

I ran into issues while trying to solve 2nd order diffeqs in combination with Unitful. MWE:

using DifferentialEquations
using Unitful
import Unitful: m, s

const ω = 0.5/s

ddy(y) = -ω^2/2 * y

prob = SecondOrderODEProblem((dy, y, p, t) -> ddy(y), 0.0m/s,0.0m, (0.0s, 10.0s))
sol = solve(prob, Tsit5()) #error
sol = solve(prob, Tsit5(), abstol=1.0e-8m/s^2, reltol = 1.0e-8m/s^2)
sol = solve(prob, abstol=1.0e-8m/s^2, reltol = 1.0e-8m/s^2) #error

So it's not only dependent on parameters, but also on which solver is being used. What this feels like to me is that there are no promotion mechanisms from standard Julia types which would be used in the various parameters, and the Number types used in Unitful, giving some issues during various operations performed by different solvers.

[devmotion]

Unfortunately, not all solvers support general number types (such as units) yet, so probably at the moment it's best to specify a suitable solver explicitly (such as in solve(prob, Tsit5())). Second order problems are even more complicated to handle since there's a mix of different units in the state vector.

I made the following observations when I tried to reproduce your example on on my computer (Julia 1.0.3 with up-to-date packages):

The lines

sol = solve(prob, Tsit5())
sol = solve(prob, Tsit5(), abstol=1.0e-8m/s^2, reltol = 1.0e-8m/s^2)

do not error, in contrast to what you experienced, it seems. There was a fix regarding units introduced to OrdinaryDiffEq (the package in which all the solvers are defined) in version 4.13.0 (current version is 4.18.3) - so maybe you used an old version? By the way, I'm actually wondering why (and whether) the second line should work at all - the units are incorrect for both the derivative dy (has units m/s) and y (has units m); IMO something like

sol = solve(prob, Tsit5(), abstol=ArrayPartition([1.0e-8m/s], [1.0e-8m]), reltol = ArrayPartition([1.0e-8m/s], [1.0e-8m]))

would be more appropriate.

Next I investigated the line

sol = solve(prob, abstol=1.0e-8m/s^2, reltol = 1.0e-8m/s^2) #error

(BTW, also here I think the units are not correct):

• The first problem is that the comparison in https://github.com/JuliaDiffEq/DifferentialEquations.jl/blob/master/src/default_arg_parsing.jl#L9 (and similar lines in that function) is not defined for units (and not for element-wise tolerances either). So probably reltol should be defined rather as minimum(x -> x / oneunit(x), o[:reltol]) (or something similar).
• The next problem are comparisons as in https://github.com/JuliaDiffEq/DifferentialEquations.jl/blob/master/src/ode_default_alg.jl#L46 (which probably are also suboptimal for other use cases such as nested arrays). Due to this check, in the example the AutoVern9 algorithm is chosen even for high tolerances.

[louisponet]

Thanks for investigating this! In my mwe, I meant that only the lines with #error behind them actually give an error. So the second line actually works as you also indicate. Sorry for not being clear about that.

[devmotion]

I meant that only the lines with #error behind them actually give an error. So the second line actually works as you also indicate. Sorry for not being clear about that.

That's what I assumed and hence I was surprised that the line

sol = solve(prob, Tsit5()) #error

does not give an error although you marked it with #error. My best guess was that the problem you experienced is already fixed. Are your packages up-to-date?

[ChrisRackauckas]

I ran it and got an error. I am not sure it's fixed?

[ChrisRackauckas]

It looks like there's a lot of other issues in this example though. (dy, y, p, t) -> ddy(y) doesn't make sense because then dy=0 for the whole evolution so the solution is constant. I assume it meant to update dy. The units look correct though. But yes, the original post definitely is using an older version since unitless tolerances are supported in this case as of that commit.

[devmotion]

I ran it and got an error. I am not sure it's fixed?

Sorry for the confusion, what I meant was: I get an error but ONLY for the last line - the first two calls of solve do not produce any errors although the first call is marked with #error.

[jonniedie]

ComponentArrays can have a heterogeneous mix of types now, so mixed units can work. Here is a similar example:

using ComponentArrays
using OrdinaryDiffEq
using LinearAlgebra
using Unitful

function newton(du, u, p, t)
    mu = 398600.4418u"km^3/s^2"
    r = norm(u.r)
    du.r = u.v
    du.v = -mu .* u.r / r^3
end

r0 = [1131.340, -2282.343, 6672.423]u"km"
v0 = [-5.64305, 4.30333, 2.42879]u"km/s"
Δt = 86400.0*365u"s"
rv0 = ComponentArray(r=r0, v=v0)

prob = ODEProblem(newton, rv0, (0.0u"s", Δt))
sol = solve(prob, Vern8(), dt=100u"s", adaptive=false)

[MatFi]

ComponentArrays can have a heterogeneous mix of types now, so mixed units can work.

Sadly one(::ComponentVector) is not implemented so it fails with implicit solvers

[jonniedie]

Sadly one(::ComponentVector) is not implemented so it fails with implicit solvers

I don't think one is supposed to work with arrays:

julia> a = [1,2,3,4]
4-element Array{Int64,1}:
 1
 2
 3
 4

julia> one(a)
ERROR: MethodError: no method matching one(::Array{Int64,1})
Closest candidates are:
  one(::Type{Missing}) at missing.jl:103
  one(::BitArray{2}) at bitarray.jl:426
  one(::Missing) at missing.jl:100
  ...
Stacktrace:
 [1] top-level scope at REPL[17]:1
 [2] include_string(::Function, ::Module, ::String, ::String) at .\loading.jl:1088

There is no reason ComponentArrays should be failing with implicit solvers. If you've got something that is failing with a ComponentVector that works with a regular Vector, go ahead and open up and issue with a MWE over here and I'll try to track down what's going on.

[cadojo]

ComponentArrays can have a heterogeneous mix of types now, so mixed units can work. Here is a similar example:

using ComponentArrays
using OrdinaryDiffEq
using LinearAlgebra
using Unitful

function newton(du, u, p, t)
    mu = 398600.4418u"km^3/s^2"
    r = norm(u.r)
    du.r = u.v
    du.v = -mu .* u.r / r^3
end

r0 = [1131.340, -2282.343, 6672.423]u"km"
v0 = [-5.64305, 4.30333, 2.42879]u"km/s"
Δt = 86400.0*365u"s"
rv0 = ComponentArray(r=r0, v=v0)

prob = ODEProblem(newton, rv0, (0.0u"s", Δt))
sol = solve(prob, Vern8(), adaptive=true, abstol=1e-12, reltol=1e-12)> 

@jonniedie Does this ComponentArray solution work with adaptive=true? I modified the last line of the block quoted code to show what I'm talking about.

If the ComponentArray fields have different units, such as in the example you gave, setting abstol=1e-12 produces DimensionError: and km are not dimensionally compatible. Is there a way around this? I saw in the documentation that you can set abstol and reltol as vectors with the same size as u0, but this does not seem to work either.

[jonniedie]

@cadojo It's trying to compare the unitless value 1e-12 with the values in the state vector and Unitful (rightly) says they're dimensionally incompatible. You can get it to work by swapping out that last line with

tol = ComponentArray(r=fill(1e-12u"km", 3), v=fill(1e-12u"km/s", 3))
sol = solve(prob, Vern8(), adaptive=true, abstol=tol, reltol=tol)

[jonniedie]

I should mention, though, that ComponentArrays with mixed types like this will be slow because arrays in general with non-concrete eltypes are slow. For this particular problem, since there are only two eltypes in the state vector, you'll get some speed gains by explicitly annotating a Union eltype. So this would look like:

r0 = [1131.340, -2282.343, 6672.423]u"km"
v0 = [-5.64305, 4.30333, 2.42879]u"km/s"
Δt = 86400.0*365u"s"

RV = Union{eltype(r0), eltype(v0)}

rv0 = ComponentArray{RV}(r=r0, v=v0)
tol = ComponentArray{RV}(r=fill(1e-12u"km", 3), v=fill(1e-12u"km/s", 3))

That makes things better, but it's still not going to be as good as calling ustrip on everything before the solve. If it's just unit checking that's important, you could run a dummy call on your derivative function before the solve call to make sure the units work. So like:

drv0 = copy(rv0)/Δt
newton(drv0, rv0, nothing, Δt)

[jonniedie]

On my machine, the times for this problem look like:

Baseline:        41.2 seconds
Union eltype:    16.0 seconds
No units:        11.3 seconds