ode_order_lowering equation ordering

[agerlach]

I used ode_order_lowering for the first time and was caught a little off guard by the resulting equation order, i.e. descending order instead of ascending. I was essentially doing the following. Which works just fine.

using ModelingToolkit

@variables t x(t)
D = Differential(t)
eqs = [(D∘D)(x) ~ -1]
sys = ode_order_lowering(ODESystem(eqs))

tspan = (0.0, 1.0)
u0 = [x => 100.0, D(x) =>0.0]
prob = ODEProblem(sys,u0,tspan,jac=true)

However, I then went to update prob via remake(prob, u0=[x_value, dx_value]) which is incorrect. It needs to be u0=[dx_value, x_value].

A couple thoughts:

1. The order seems unnatural. I don't know if I've ever seen an ODE lowered in that order in a textbook. Obviously, the order doesn't matter in principle, but it does impact the interface for updating the resulting numerical problem. At a minimum this should be documented.
2. Possibly add a "sort order" kwarg to ode_order_lowering. However, I think below is a better option
3. Treat u0 = [x => 100.0, D(x) =>0.0] as a template for the interface, e.g. reorder the equations based on the order in the array. However, I can see situations where the order of the equations impacts performance and using a symbolic transformation to put it in a favorable form would conflict w/ this, e.g. grouping uncertain states for expectation() to prevent cache miss.

[ChrisRackauckas]

Note that in general users shouldn't rely on ordering. You can even now do sol[D(x)], so in the end I don't think users should be expecting an order and instead we should give interfaces that are all order-independent. In fact, with many passes variables can even be deleted, in which case we give back values which were never computed in the solver itself!

[agerlach]

Agreed completely. I was of the misconception that all symbolic was discarded after ODEProblem. What do you envision for a remake interface then? I don't think a user would want to form arrays of pairs to implement a prob_func for ensembles. Or does something already exist?

However, with that said, if a user ever builds the functions to use directly, then order will matter to the user.

[ChrisRackauckas]

We do need a nice way to handle remake. That's one of the issues coming up this week.

We are doing something of a defaults system, i.e.

remake(prob,p = [delta => 2.0])

which would then change only delta, and reuse the others. We were just discussing this the other day actually. There's a way to get to this but it will require #777

[agerlach]

Cool! I can see the verbosity of using pairs getting annoying though for large systems. Does it make any sense have the option to pass the variables and values separately, e.g.

u0s =  # Array of ic
u0vars = [x, D(x), D(D(x))]
function prob_func(prob,i,repeat) 
    remake(prob,u0=u0s[i], u0vars = u0vars )
end

or would you just propose

remake(prob,u0=Pair.(u0vars, u0s[i] ))

Seems like needing to copy to a Pair could lead to a big performance hit.

[ChrisRackauckas]

Pairs are stack allocated so that wouldn't be a huge issue. And in this form you'd only have to provide the u0's that you're changing, so the list shouldn't be huge? Though I do want to expose an index map somehow to allow the user to get the permutation from the original variables.