Vector Equations with ModelingToolkit and MethodOfLines (Updated Syntax)

I have attempted to update some code to the preferred updated syntax u(..)[1:n] since u[1:n](..) is slated for deprecation. When I modify this in the code, I get a new error.

I have attached an example problem from adapted from the tutorial with the fixes from a previous discussion (Vector equations LHS/RHS missing with MTK and MethodOfLines)

The error seems to have something to do with the variable map step, and specifically (I think) to something going wrong with getting the properties of the independent variable t in this case. Is there any workaround for this issue?

In the meantime, is the syntax deprecation of u(..)[1:n] likely to happen soon, and are there any resources I can use to learn how to figure out what is going wrong deep in packages like this?

using DifferentialEquations, ModelingToolkit, MethodOfLines, DomainSets

# Parameters, variables, and derivatives
n_comp = 2
@parameters t, x, p[1:n_comp], q[1:n_comp]
@variables u(..)[1:n_comp]
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2
params = Symbolics.scalarize(reduce(vcat,[p .=> [1.5, 2.0], q .=> [1.2, 1.8]]))
# 1D PDE and boundary conditions

eqs  = [Dt(u(t, x)[i]) ~ p[i] * Dxx(u(t, x)[i]) for i in 1:n_comp]

bcs = [[u(0, x)[i] ~ q[i] * cos(x),
        u(t, 0)[i] ~ sin(t),
        u(t, 1)[i] ~ exp(-t) * cos(1),
        Dx(u(t,0)[i]) ~ 0.0] for i in 1:n_comp]
bcs_collected = reduce(vcat, bcs)

# Space and time domains
domains = [t ∈ Interval(0.0, 1.0),
           x ∈ Interval(0.0, 1.0)]

# PDE system

@named pdesys = PDESystem(eqs, bcs_collected, domains, [t, x], [u(t, x)[i] for i in 1:n_comp], Symbolics.scalarize(params))


# Method of lines discretization
dx = 0.1
order = 2
discretization = MOLFiniteDifference([x => dx], t; approx_order = order)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization) #error occurs here

# Solve ODE problem
    sol = solve(prob, Tsit5(), saveat=0.2)

The error I get from attempting discretization is:

ERROR: ArgumentError: invalid index: nothing of type Nothing
Stacktrace:
  [1] to_index(i::Nothing)
    @ Base .\indices.jl:300
  [2] to_index(A::Vector{Symbolics.VarDomainPairing}, i::Nothing)
    @ Base .\indices.jl:277
  [3] to_indices
    @ .\indices.jl:333 [inlined]
  [4] to_indices
    @ .\indices.jl:325 [inlined]
  [5] getindex
    @ .\abstractarray.jl:1241 [inlined]
  [6] (::PDEBase.var"#48#57"{Vector{Symbolics.VarDomainPairing}})(x::SymbolicUtils.BasicSymbolic{Vector{Real}})
    @ PDEBase \.julia\packages\PDEBase\rQehP\src\variable_map.jl:35
  [7] iterate
    @ .\generator.jl:47 [inlined]
  [8] _collect(c::Vector{Any}, itr::Base.Generator{Vector{Any}, PDEBase.var"#48#57"{Vector{Symbolics.VarDomainPairing}}}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})
    @ Base .\array.jl:807
  [9] collect_similar(cont::Vector{Any}, itr::Base.Generator{Vector{Any}, PDEBase.var"#48#57"{Vector{Symbolics.VarDomainPairing}}})
    @ Base .\array.jl:716
 [10] map(f::Function, A::Vector{Any})
    @ Base .\abstractarray.jl:2933
 [11] VariableMap(pdesys::PDESystem, disc::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
    @ PDEBase \.julia\packages\PDEBase\rQehP\src\variable_map.jl:34
 [12] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
    @ PDEBase \.julia\packages\PDEBase\rQehP\src\symbolic_discretize.jl:15
 [13] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}; analytic::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ PDEBase \.julia\packages\PDEBase\rQehP\src\discretization_state.jl:58
 [14] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
    @ PDEBase \.julia\packages\PDEBase\rQehP\src\discretization_state.jl:55
 [15] top-level scope
    @ \test_for_debug.jl:38
1 Like

I don’t believe that’s supported yet. @xtalax ?

Ah, okay. Thanks! Is the legacy syntax from Symbolics expected to be fully deprecated soon, or should I be safe for a while still?

I thought this was supported, it’s definitely possible to get working. Let me take a closer look.

EDIT: Alright I’m going to push a fix so this works, in the meantime I think the old syntax will be fine for a while.