Multi-unit handling

[ChrisRackauckas]

We want codes like this to work:

using DifferentialEquations, Unitful

function g(dydt, y, p, t)
dydt[1] = -p[1] / p[2] * sin(y[2])
dydt[2] = y[1]
end

tspan = (0.0u"s",4.0u"s")
Y0=[0. * u"rad/s", .5u"rad"]
p = [9.8u"m/s^2",1u"m"]
prob = ODEProblem(g,Y0, tspan,p)

#Pass to solvers
sol = solve(prob,Tsit5())

I'm not sure this is very easy since it requires the units in the RK coefficients to be multi-united dependent on the part of the linear combination they are in. This may need a new array type and a multi-unit element type scalar. xref: https://github.com/JuliaDiffEq/DiffEqTutorials.jl/issues/27

Right now we can actually handle this case via an ArrayPartition and using a SecondOrderODEProblem solver. So if the ODE partitions into state units and state/time units (for the velocities), it's handled fine. But in general, allowing arbitrary units for the terms of the "vector" would be ideal. A standard Julia vector cannot do this well though, so we need new types to handle the arbitrary (more than 2 units) case.

[ChrisRackauckas]

@ajkeller34 an MWE for what I am talking about goes like this. We have united arrays A and B which exactly match types. We want to define scalars a and b where typeof(a) == typeof(b) == eltype(A), so that

@. dt*(a*A + b*B)

works. If A all has the same units, this already works. If A is multi-united, we need for a to somehow be a scalar / abstract array mix s.t. a[i] has the right units (but always the same value, it's just a scalar with more type information), and so that broadcast works. So it's a scalar with type information and some overloads to make array-like behavior work, but array-like behavior only effects the units.

As for the data structure of A, I said array but we should also consider a tuple. A tuple would be strictly typed but immutable and not scale well to a larger size. Ideally we could also have this as a heterogeneously typed array, backed just by a single contiugous array of non-united T (since it's just type information tagged on). If it is just a Julia array it could still work, relying heavily on Julia's small union optimizations.

[ChrisRackauckas]

In the example above a and b should actually be unitless, so it's easier in this case. But I'm still not sure how to represent the Array properly. Use union of units?

[LebedevRI]

Apologies for necroposting.. A number of issues were closed as being duplicate of this bug, and then this bug has been closed with no comments, no linked changes, no linked pull requests. I'm now also hitting this issue, and i'm having trouble deciphering the outcome of this issue? Should this be now working and i'm doing it wrong, or the issue was closed as wontfix? Perhaps the docs (http://tutorials.juliadiffeq.org/html/type_handling/03-unitful.html) could use some more words.

[ChrisRackauckas]

That tutorial was written in 2016 (the time stamps are just newer because of the auto-updating infrastructure). In 202(1), I think the way that Unitful.jl expresses units is a really bad idea and can never actually be used. A completely different unitful array implementation is required. If it wants to do it all in the type domain, it needs to commit to having its own array type that specializes on the element types of each individual portion, using similar machinery to https://github.com/SciML/LabelledArrays.jl . The current implementation of unitful arrays has some fundamental issues there, like the fact that its eltype is not well-defined (because each element can have different units and those units are types), so it has no zero element, it has no linear algebra defined, etc. That's not the kind of library that can robustly work with a differential equation solver beyond what is shown in that tutorial.

I do plan to, as somewhat of a sideshow, make my own implementation of a correct unitful array so I will reopen this, but I don't have a timeline on it since it doesn't intersect with any of my real work.

[LebedevRI]

Thank you for a speedy response! It is good to know that this simply still doesn't really work, and the problem is mostly not on my side!

[RomeoV]

Is there any progress on this by any chance? It has come up again in a recent project of mine.

[ChrisRackauckas]

Use ModelingToolkit https://mtk.sciml.ai/dev/basics/Validation/