ODE systems with heterogenous units

[helgee]

As discussed here: https://discourse.julialang.org/t/differentialequations-jl-and-systems-with-heterogenous-units/2952/4

The following example does not work:

function newton!(t, y, dy, μ)
    r = norm(y[1:3])
    dy[1:3] = y[4:6]
    dy[4:6] = -μ * y[1:3] / r^3
end

function test2()
    r0 = [1131.340, -2282.343, 6672.423]u"km"
    v0 = [-5.64305, 4.30333, 2.42879]u"km/s"
    Δt = 86400.0*365u"s"
    mu = 398600.4418u"km^3/s^2"
    rv0 = [r0; v0]

    prob = ODEProblem((t, y, dy) -> newton!(t, y, dy, mu), rv0, (0.0u"s", Δt))
    solve(prob, Vern8())[end]
end

[ChrisRackauckas]

Few things to work out here. First of all, there is a problem due to heterogeneous typing. A simpler version of this problem is illustrated by an ODE where

u0 = Vector{Real}(2)
u0[1] = 1; u0[2] = 2.0

Put that into any ODEProblem and it will fail. The reason is because many constants try to generate themselves as eltype(u0). But that will fail if eltype(u0) is not a concrete type. So first things first: it needs to find out how to handle non-concrete eltypes.

But then secondly, it would have to handle units. Some things which would error here is, what is the norm of a vector which has different units? That will cause adaptive timestepping to fail unless the user gives a dispatch. I feel like a new library HeterogeneousArrays.jl could be built to handle this. It would also have to handle the fact that the array is heterogeneously typed, and thus is very slow to access. It could be made faster if an iterator for broadcast is made using a generated function and it treats the separately typed parts as different fields.

But this is a ways off. These are some pretty fundamental problems with multiply-typed initial conditions, and it will take some sweat to make this compatible. I don't plan on having this for awhile. Sorry if I got your hopes up earlier. But, I can forsee a future where this could work.

[helgee]

A few more things to consider and notes to myself:

Not only has y heterogenous units but dy also has different heterogenous units, e.g.

unit(r) == u"km"
unit(v) == u"km/s"

unit(dr) == u"km/s"
unit(dv) == u"km/s^2"

Another issue is that the HeterogenousArray type as a subtype of AbstractArray approach does not work because AbstractArray expects a consistent eltype.

In my specific case I would be possible to partition the problem somehow because the velocity equations dy[1:3] = y[4:6] are always the same and only the acceleration equation dy[4:6] = ... need to be defined in the RHS.

[ChrisRackauckas]

Two things are folded together here. First of all, this would be fixed by using a method for 2nd order ODEs. In that case, the velocity variables could get their units as the position units divided by the time units. The integrator would be better off keeping these separate, and so this problem would be avoided. That's in the works:

SciML/DiffEqBase.jl#14

but isn't a priority (for me. I'm sorry, I just don't have any 2nd order ODEs in my research to justify it!). I could help you navigate the system to make it a reality though.

But secondly, this is starting to make me question the "units as type" approach, at least for this application. The whole point of the approach is that, since it's all type information, unit-checking is done by dispatch and so it has zero runtime cost. However, putting different units on different variables seems like a very sensible thing here. But if the moment you do that, you lose concreteness of the eltype and thus get bad performance, that makes it not a good idea. I think it would be nice to have a

immutable NumberWithUnits{T} <: Number
  x::T
  units::Symbol
end

and so you'd have Vector{NumberWithUnits} with the units stored as a field. Then the arithmetic x+y would check x.units and y.units for compatibility and then do the arithmetic (it can even do conversions automatically!). This does have a runtime cost, but it is more flexible. For this application it would work better.

[ChrisRackauckas]

I have a few different solutions here.

First, with the next release of RecursiveArrayTools.jl, you'll be able to create an ArrayPartition. You can use this to create an array where different parts have different types:

u0 = 3ones(4)
v0 = ones(4)
u = ArrayPartition((u0,v0))
f = function (t,u,du)
  du.x[1] .= u.x[2]
  du.x[2] .= -2u.x[1]
end
prob = ODEProblem(f,u,(0.0,5.0))
sol = solve(prob,Euler(),dt=1/100)

It's just an array formed by tuples of arrays.

But secondly

u0 = 3ones(4)
v0 = ones(4)
f2 = function (t,u,v,dv)
  dv .= -2u
end
prob = SecondOrderODEProblem(f2,u0,v0,(0.0,5.0))
sol2 = solve(prob,SymplecticEuler(),dt=1/100)

with the next release you'll be able to directly solve 2nd order problems (underneath, it is using ArrayPartition). More generally, there is PartitionedODEProblem. See

http://docs.juliadiffeq.org/latest/types/refined_ode_types.html#Mathematical-Specification-of-a-Partitioned-ODE-Problem-1
http://docs.juliadiffeq.org/latest/solvers/refined_ode_solve.html#Functional-Forms-1

Now the PartitionedODEProblem and SecondOrderODEProblem solvers "know" about this form, and so they use it with full efficiency. But there aren't many solvers there yet. The other ODE solvers will work, but not with the same efficiency. However, they will be fully efficient if I figure out one last problem:

https://discourse.julialang.org/t/efficient-broadcast-of-heterogeneous-abstractarrays/3244

So this is completed, but I'll leave it open until this is all released.

[ChrisRackauckas]

https://discourse.julialang.org/t/efficient-broadcast-of-heterogeneous-abstractarrays/3244/9?u=chrisrackauckas

Using ArrayPartitions is now just as efficient with broadcasting as an array, even with heterogeneous typing. Thus this is now a really good solution. RecursiveArrayTools.jl needs some docs on it.

[ChrisRackauckas]

function newton!(t, y, dy, μ)
    r = norm(y.x[1])
    dy.x[1] .= y.x[2]
    dy.x[2] .= -μ .* y.x[1] / r^3
end

using Unitful, RecursiveArrayTools, DiffEqBase, OrdinaryDiffEq

r0 = [1131.340, -2282.343, 6672.423]u"km"
v0 = [-5.64305, 4.30333, 2.42879]u"km/s"
Δt = 86400.0*365u"s"
mu = 398600.4418u"km^3/s^2"
rv0 = ArrayPartition(r0,v0)

prob = ODEProblem((t, y, dy) -> newton!(t, y, dy, mu), rv0, (0.0u"s", Δt))
sol = solve(prob, Vern8(),dt=100u"s",adaptive=false)

[ChrisRackauckas]

If you use an ArrayPartition you'll be good. Two things still need to be done:

1. Initial dt calculation needs to be compatible with heterogeneous units
2. Adaptivity needs to be compatible with heterogeneous units

(2) will come for free with v0.6 when everything changes to broadcast. Can't do that right now though because it would massively slow down v0.5. (1) I'll finish fixing up soon. This example will be added to the docs in a second on AbstractArray usage.

Note this is slow on v0.5 because of the lack of broadcast (the indexing is really bad) but will end up fine on v0.6 after the broadcast change.

[ChrisRackauckas]

Initial dt is fixed. Adaptivity will be free with v0.6 changes. For now, set adaptive=false and give a dt.