Discretize MOLFiniteDifference with more than one dependent variable

[michaelsachs]

I am trying to use ModelingToolkit for a 1D diffusion simulation with more than one dependent variable, but am getting a BoundsError when discretizing using MOLFiniteDifference.

The following code with one dependent variable works fine:

using OrdinaryDiffEq, ModelingToolkit, DiffEqOperators

# Parameters, variables, and derivatives
@parameters t x
@variables n(..) 
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

# 1D PDE and boundary conditions
eqs  = Dt(n(t,x)) ~ Dxx(n(t,x))

bcs = [n(0,x) ~ 2x,
        Dx(n(t,0)) ~ 0,
        Dx(n(t,10)) ~ 0]

# Space and time domains
domains = [t ∈ IntervalDomain(0.0,100.0),
           x ∈ IntervalDomain(0.0,10.0)]

# PDE system
pdesys = PDESystem(eqs,bcs,domains,[t,x],[n(t,x)])

# Method of lines discretization
Δx = 0.1
discretization = MOLFiniteDifference([x=>Δx],t)

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

# Solve ODE problem
@time sol = solve(prob,Tsit5(),saveat=1.0)

However, when adding a second dependent variable (while keeping the system 1D) I am getting a BoundsError in the discretize(pdesys,discretization) step:

using OrdinaryDiffEq, ModelingToolkit, DiffEqOperators

# Parameters, variables, and derivatives
@parameters t x
@variables n(..) p(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

# 1D PDE and boundary conditions
eqs  = [Dt(n(t,x)) ~ Dxx(n(t,x)),
        Dt(p(t,x)) ~ Dxx(p(t,x))]

bcs = [n(0,x) ~ 2x,
        Dx(n(t,0)) ~ 0,
        Dx(n(t,10)) ~ 0,
        p(0,x) ~ 2x,
        Dx(p(t,0)) ~ 0,
        Dx(p(t,10)) ~ 0]

# Space and time domains
domains = [t ∈ IntervalDomain(0.0,100.0),
           x ∈ IntervalDomain(0.0,10.0)]

# PDE system
pdesys = PDESystem(eqs,bcs,domains,[t,x],[n(t,x),p(t,x)])

# Method of lines discretization
Δx = 0.1
discretization = MOLFiniteDifference([x=>Δx],t)

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

# Solve ODE problem
@time sol = solve(prob,Tsit5(),saveat=1.0)

With the stacktrace being as follows:

ERROR: LoadError: BoundsError: attempt to access 1-element Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}} at index [2]
Stacktrace:
  [1] getindex
    @ .\array.jl:801 [inlined]
  [2] left_weights
    @ ~\.julia\packages\DiffEqOperators\Fmo6r\src\MOLFiniteDifference\MOL_discretization.jl:57 [inlined]
  [3] (::DiffEqOperators.var"#204#243"{DiffEqOperators.var"#right_idxs#242"{Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}}, Vector{CartesianIndex{1}}, DiffEqOperators.var"#right_weights#238"{MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num}, Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}}, DiffEqOperators.var"#left_weights#237"{MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num}, Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}}})(::Tuple{Int64, Vector{Num}})
    @ DiffEqOperators .\none:0
  [4] iterate
    @ .\generator.jl:47 [inlined]
  [5] collect_to!(dest::Vector{Vector{Num}}, itr::Base.Generator{Base.Iterators.Enumerate{Vector{Vector{Num}}}, DiffEqOperators.var"#204#243"{DiffEqOperators.var"#right_idxs#242"{Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}}, Vector{CartesianIndex{1}}, DiffEqOperators.var"#right_weights#238"{MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num}, Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}}, DiffEqOperators.var"#left_weights#237"{MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num}, Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}}}}, offs::Int64, st::Tuple{Int64, Int64})
    @ Base .\array.jl:724
  [6] collect_to_with_first!(dest::Vector{Vector{Num}}, v1::Vector{Num}, itr::Base.Generator{Base.Iterators.Enumerate{Vector{Vector{Num}}}, DiffEqOperators.var"#204#243"{DiffEqOperators.var"#right_idxs#242"{Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}}, Vector{CartesianIndex{1}}, DiffEqOperators.var"#right_weights#238"{MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num}, Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}}, DiffEqOperators.var"#left_weights#237"{MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num}, Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}}}}, st::Tuple{Int64, Int64})
    @ Base .\array.jl:702
  [7] collect(itr::Base.Generator{Base.Iterators.Enumerate{Vector{Vector{Num}}}, DiffEqOperators.var"#204#243"{DiffEqOperators.var"#right_idxs#242"{Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}}, Vector{CartesianIndex{1}}, DiffEqOperators.var"#right_weights#238"{MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num}, Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}}, DiffEqOperators.var"#left_weights#237"{MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num}, Vector{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}}}}})
    @ Base .\array.jl:683
  [8] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num})
    @ DiffEqOperators ~\.julia\packages\DiffEqOperators\Fmo6r\src\MOLFiniteDifference\MOL_discretization.jl:63
  [9] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num})
    @ DiffEqOperators ~\.julia\packages\DiffEqOperators\Fmo6r\src\MOLFiniteDifference\MOL_discretization.jl:152

Looking into MOL_discretization.jl, the problem seems to be that space is a one element vector when it should contain two elements for the following piece of code to run successfully (as depvars also contains two elements). Not sure if this is a bug or just me getting the syntax wrong?

# 1D system
left_weights(j) = DiffEqOperators.calculate_weights(discretization.upwind_order, space[j][1], space[j][1:2])
right_weights(j) = DiffEqOperators.calculate_weights(discretization.upwind_order, space[j][end], space[j][end-1:end])
central_neighbor_idxs(i,j) = [i+CartesianIndex([ifelse(l==j,-1,0) for l in 1:nspace]...),i,i+CartesianIndex([ifelse(l==j,1,0) for l in 1:nspace]...)]
left_idxs = central_neighbor_idxs(CartesianIndex(2),1)[1:2]
right_idxs(j) = central_neighbor_idxs(CartesianIndex(length(space[j])-1),1)[end-1:end]
# Constructs symbolic spatially discretized terms of the form e.g. au₂ - bu₁
derivars = [[dot(left_weights(j),depvar[left_idxs]), dot(right_weights(j),depvar[right_idxs(j)])]
            for (j, depvar) in enumerate(depvars)]

[valentinsulzer]

What's missing is #350 to handle the higher-dimensional Neumann BCs. It should currently work with Dirichlet BCs, see https://github.com/SciML/DiffEqOperators.jl/blob/master/test/MOL/MOL_2D_Diffusion.jl

[michaelsachs]

Thanks for the reply! I tried Dirichlet BCs and actually encountered the same error.

Is the dimensionality not related to the numer of spatial dimensions? Both of my examples should be 1-dimensional, as they only have one x dimension. In the example that throws the BoundsError I am just trying to add another variable (with the aim of eventually making n(t,x) and p(t,x) interact) but no additional spatial dimension.

[valentinsulzer]

Oh my bad, I misread the question, yeah that is a genuine bug.

Since both p and n are defined on the same dependent variable x, space should only have 1 element, but I'll try to see what else is going wrong.