Unable to Index into a steadystate solution using symbols

[yonatanwesen]

I am unable to index into the steadystate solution using symbols.

Here is MWE:

using ModelingToolkit,MethodOfLines, DomainSets
using OrdinaryDiffEq, SteadyStateDiffEq

@parameters t x y
@parameters D b
@variables A(..) S(..)
Dx = Differential(x)
Dy = Differential(y)
Dt = Differential(t)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
#2d PDE and bcs
eqs = [Dt(A(t,x,y)) ~ Dxx(A(t,x,y)) + Dyy(A(t,x,y))  + A(t,x,y)^ 2 /S(t,x,y) - b*A(t,x,y), #+ Dyy(A(t,x))
       Dt(S(t,x,y)) ~ D*Dxx(S(t,x,y)) + D*Dyy(S(t,x,y)) + A(t,x,y)^ 2 - S(t,x,y)] #+ D*Dyy(S(t,x))
A0, S0 = homo_turing(0.35)
d = Normal(0.0,1.0)
bcs = [A(0, x,y) ~ A0 * (1 + 0.1*rand(d)) + 0.1*exp(-1.0*x^2*y^2),
S(0, x,y) ~ S0*(1 +  0.1*rand(d))  + 0.1* exp(-1.0*x^2*y^2) , #
Dx(A(t, 0,y)) ~ 0.0,
Dx(A(t, 10,y)) ~ 0.0,
Dx(S(t, 0,y)) ~ 0.0,
Dx(S(t, 10,y)) ~ 0.0,
Dy(A(t, x,0)) ~ 0.0,
Dy(A(t, x,10)) ~ 0.0,
Dy(S(t, x,0)) ~ 0.0,
Dy(S(t, x,10)) ~ 0.0
]
#domain
domains = [t ∈ Interval(0, 1.0), x ∈ Interval(0.0, 10.0),y ∈ Interval(0.0, 10.0)]
dx = 1.0
dy = 1.0
order = 2
#construct the pde sys
@named pdesys = PDESystem(eqs,
    bcs,
    domains,
    [t, x,y],
    [A(t, x,y), S(t, x,y)],
    [D => 30.0, b=> 0.35])
discretization = MOLFiniteDifference([x => dx,y =>dy], t, approx_order = order)
prob = discretize(pdesys, discretization)
steadystateprob = SteadyStateProblem(prob)
@time steadystate = solve(steadystateprob, DynamicSS(AutoTsit5(Rosenbrock23())))

steadystate[A(t,x,y)] # broken
steadystate[A(x,y)]  # broken

Doing steadystate[A(t,x,y)] or steadystate[A(x,y)] results in the following error message

Error & Stacktrace ⚠️

ERROR: ArgumentError: A(t, x, y) is neither an observed nor a state variable.
Stacktrace:
 [1] build_explicit_observed_function(sys::ODESystem, ts::Num; inputs::Nothing, expression::Bool, output_type::Type, checkbounds::Bool, drop_expr::typeof(RuntimeGeneratedFunctions.drop_expr), ps::Vector{…}, throw::Bool)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/N1HdY/src/systems/diffeqs/odesystem.jl:375
 [2] build_explicit_observed_function
   @ ~/.julia/packages/ModelingToolkit/N1HdY/src/systems/diffeqs/odesystem.jl:315 [inlined]
 [3] #541
   @ ~/.julia/packages/ModelingToolkit/N1HdY/src/systems/diffeqs/abstractodesystem.jl:476 [inlined]
 [4] get!(default::ModelingToolkit.var"#541#556"{…}, h::Dict{…}, key::SymbolicUtils.BasicSymbolic{…})
   @ Base ./dict.jl:479
 [5] (::ModelingToolkit.var"#649#generated_observed#555"{…})(::Num, ::Vector{…}, ::Vararg{…})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/N1HdY/src/systems/diffeqs/abstractodesystem.jl:475
 [6] observed
   @ SciMLBase ~/.julia/packages/SciMLBase/XEPyX/src/solutions/solution_interface.jl:184 [inlined]
 [7] getindex(A::SciMLBase.NonlinearSolution{…}, sym::Num)
   @ SciMLBase ~/.julia/packages/SciMLBase/XEPyX/src/solutions/solution_interface.jl:174
 [8] top-level scope
   @ REPL[41]:1
Some type information was truncated. Use `show(err)` to see complete types.

Environment (please complete the following information):

• The version of the packages I am using is below. I am using older version of ModelingToolkit due to latest version being broken as show in https://github.com/SciML/PDEBase.jl/issues/19

[5b8099bc] DomainSets v0.6.7
⌃ [94925ecb] MethodOfLines v0.10.1
⌃ [961ee093] ModelingToolkit v8.71.2
    [9672c7b4] SteadyStateDiffEq v1.16.1
⌃ [1dea7af3] OrdinaryDiffEq v6.59.0

• Output of versioninfo()

Julia Version 1.10.0-rc1
Commit 5aaa9485436 (2023-11-03 07:44 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 12 × Intel(R) Xeon(R) CPU E5-2603 v4 @ 1.70GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, broadwell)
  Threads: 1 on 12 virtual cores
Environment:
  JULIA_EDITOR = code
  JULIA_NUM_THREADS = 

cc @ChrisRackauckas

[xtalax]

If you wanted to solve steady state, you need to write your problem independent of time and solve with a NonlinearProblem, see the docs: https://docs.sciml.ai/MethodOfLines/stable/tutorials/heatss/

[yonatanwesen]

The reason I formulated the problem as steadystate problem rather than NonlinearProblem is because PDEbase doesn't allow setting initial conditions for u0. Even if I reformulate steadystate as NonlinearProblem, I run into the same issue

[xtalax]

I'll add that feature then, thanks for the heads up!