Solution interface does not work when using init() and step!()

[aplund]

It seems that one cannot use the solution using init() and step!() from the integrator interface.

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets

@variables x y t
@variables u(..) v(..)
Dt = Differential(t)
Dx = Differential(x)
Dy = Differential(y)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

∇²(u) = Dxx(u) + Dyy(u)

brusselator_f(x, y, t) = (((x-0.3)^2 + (y-0.6)^2) <= 0.1^2) * (t >= 1.1) * 5.

x_min = y_min = t_min = 0.0
x_max = y_max = 1.0
t_max = 11.5

α = 10.

u0(x,y,t) = 22(y*(1-y))^(3/2)
v0(x,y,t) = 27(x*(1-x))^(3/2)

eq = [Dt(u(x,y,t)) ~ 1. + v(x,y,t)*u(x,y,t)^2 - 4.4*u(x,y,t) + α*∇²(u(x,y,t)) + brusselator_f(x, y, t),
       Dt(v(x,y,t)) ~ 3.4*u(x,y,t) - v(x,y,t)*u(x,y,t)^2 + α*∇²(v(x,y,t))]

domains = [x ∈ Interval(x_min, x_max),
              y ∈ Interval(y_min, y_max),
              t ∈ Interval(t_min, t_max)]

# Periodic BCs
bcs = [u(x,y,0) ~ u0(x,y,0),
       u(0,y,t) ~ u(1,y,t),
       u(x,0,t) ~ u(x,1,t),

       v(x,y,0) ~ v0(x,y,0),
       v(0,y,t) ~ v(1,y,t),
       v(x,0,t) ~ v(x,1,t)]

@named pdesys = PDESystem(eq,bcs,domains,[x,y,t],[u(x,y,t),v(x,y,t)])

N = 32

order = 2 # This may be increased to improve accuracy of some schemes

# Integers for x and y are interpreted as number of points. Use a Float to directtly specify stepsizes dx and dy.
discretization = MOLFiniteDifference([x=>N, y=>N], t, approx_order=order)

# Convert the PDE problem into an ODE problem
println("Discretization:")
@time prob = discretize(pdesys,discretization)

state = init(prob, TRBDF2())
step!(state, 0.1)

discrete_x = state.sol[x]
discrete_y = state.sol[y]
discrete_t = state.sol[t]

solu = state.sol[u(x, y, t)]
solv = state.sol[v(x, y, t)]

The error given is

ERROR: LoadError: ArgumentError: x is neither an observed nor a state variable.

with an extraordinarily long stack trace starting at build_explicit_observed_function.

[xtalax]

There is no hook to wrap the solution object in step!, this is by design as it isn't an optimal store for ODE solving. If you can construct an ODESolution from your state and call SciMLBase.wrap_sol on it with the system metadata you can get this functionality though

[ChrisRackauckas]

This is something we should address at some point. In theory the integrator indexing should be handled too.