SciML / SciMLBase.jl

The Base interface of the SciML ecosystem
https://docs.sciml.ai/SciMLBase/stable
MIT License
118 stars 91 forks source link

Symbolic indexing fails for certain `SteadyStateProblem`s #660

Closed TorkelE closed 2 months ago

TorkelE commented 2 months ago

Certain types of indexing fails for steady-state problems. It only seems to be a case when some variables are eliminated as observables, and then one tries to access these:

using OrdinaryDiffEq, SteadyStateDiffEq, ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using ModelingToolkit: getu, getp, setu, setp
@parameters p d
@variables X(t) X2(t)
eqs = [
    D(X) ~ p - d*X,
    X2 ~ 2*X
]
@mtkbuild osys = ODESystem(eqs, t)

u0 = [X => 0.1]
p = [p => 1.0, d => 0.2]
prob = SteadyStateProblem(osys, u0, p)

# Get u values (including observables).
prob[X] == prob[osys.X] == prob[:X] == 0.1
prob[X2] == prob[osys.X2] == prob[:X2] == 0.2
prob[[X,X2]] == prob[[osys.X,osys.X2]] == prob[[:X,:X2]] == [0.1, 0.2]
# prob[(X,X2)] == prob[(osys.X,osys.X2)] == prob[(:X,:X2)] == (0.1, 0.2) # Does not work for any problem types.
getu(prob, X)(prob) == getu(prob, osys.X)(prob) == getu(prob, :X)(prob) == 0.1
getu(prob, X2)(prob) == getu(prob, osys.X2)(prob) == getu(prob, :X2)(prob) == 0.2 # ERROR: type SteadyStateProblem has no field tspan
getu(prob, [X,X2])(prob) == getu(prob, [osys.X,osys.X2])(prob) == getu(prob, [:X,:X2])(prob) == [0.1, 0.2] # ERROR: type SteadyStateProblem has no field tspan
getu(prob, (X,X2))(prob) == getu(prob, (osys.X,osys.X2))(prob) == getu(prob, (:X,:X2))(prob) == (0.1, 0.2) # ERROR: type SteadyStateProblem has no field tspan

For some cases it works fine though.