Dimensionality-related error

[schmidtjonathan]

Goal

Benchmark the MOL-discretized Laplace operator on a large-ish grid against my own implementation.
Therefore, simply set up a PDESystem for the Laplace operator and use MethodOfLines.jl to discretize it.

Code

using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets
@parameters t x y
@variables u(..)
Dt = Differential(t)
Dxx = Differential(x)^2
Dyy = Differential(y)^2

eq  = Dt(u(t, x, y)) ~ 0.01 * (Dxx(u(t, x, y)) + Dyy(u(t, x, y)))

domains = [x ∈ Interval(0.0, 1.0),
        y ∈ Interval(0.0, 1.0),
        t ∈ Interval(0.0, 1.0)]

initf(t, x, y) = x   # not relevant

bcs = [u(0.0, x, y) ~ initf(t, x, y),
        u(t, 0.0, y) ~ 0.0,
        u(t, x, 0.0) ~ 0.0,
        u(t, 1.0, y) ~ 0.0,
        u(t, x, 1.0) ~ 0.0,
]

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

N = 100
dx = 1.0 / N
discretization = MOLFiniteDifference([x => dx, y => dx], t, approx_order=2)

prob = discretize(pdesys, discretization)

Error

Then, calling

out = similar(prob.u0)
prob.f(out, prob.u0, nothing, 0.0)

yields

ERROR: UndefVarError: (u(t))[1, 2] not defined
Stacktrace:
  [1] macro expansion
    @ ~/.julia/packages/SymbolicUtils/qulQp/src/code.jl:394 [inlined]
  [2] macro expansion
    @ ~/.julia/packages/Symbolics/J8IHJ/src/build_function.jl:504 [inlined]
  [3] macro expansion
    @ ~/.julia/packages/SymbolicUtils/qulQp/src/code.jl:351 [inlined]
  [4] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/KrkGo/src/RuntimeGeneratedFunctions.jl:129 [inlined]
  [5] macro expansion
    @ ./none:0 [inlined]
  [6] generated_callfunc
    @ ./none:0 [inlined]
  [7] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), Symbolics.var"#_RGF_ModTag", Symbolics.var"#_RGF_ModTag", (0x31798052, 0xa8339a20, 0xe7ba35d5, 0xa26036e4, 0x09c47050)})(::Vector{Float64}, ::Vector{Float64}, ::Vector{Any}, ::Float64)
    @ RuntimeGeneratedFunctions ~/.julia/packages/RuntimeGeneratedFunctions/KrkGo/src/RuntimeGeneratedFunctions.jl:117
  [8] macro expansion
    @ ~/.julia/packages/Symbolics/J8IHJ/src/build_function.jl:367 [inlined]
  [9] macro expansion
    @ ~/.julia/packages/SymbolicUtils/qulQp/src/code.jl:351 [inlined]
 [10] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/KrkGo/src/RuntimeGeneratedFunctions.jl:129 [inlined]
 [11] macro expansion
    @ ./none:0 [inlined]
 [12] generated_callfunc(::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), Symbolics.var"#_RGF_ModTag", Symbolics.var"#_RGF_ModTag", (0x65002b91, 0x50030b6e, 0x0e77a634, 0x5a060fbd, 0xeec46872)}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Any}, ::Float64)
    @ Symbolics ./none:0
 [13] RuntimeGeneratedFunction
    @ ~/.julia/packages/RuntimeGeneratedFunctions/KrkGo/src/RuntimeGeneratedFunctions.jl:117 [inlined]
 [14] macro expansion
    @ ~/.julia/packages/Symbolics/J8IHJ/src/build_function.jl:367 [inlined]
 [15] macro expansion
    @ ~/.julia/packages/SymbolicUtils/qulQp/src/code.jl:351 [inlined]
 [16] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/KrkGo/src/RuntimeGeneratedFunctions.jl:129 [inlined]
 [17] macro expansion
    @ ./none:0 [inlined]
 [18] generated_callfunc(::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), Symbolics.var"#_RGF_ModTag", Symbolics.var"#_RGF_ModTag", (0x0d72d0ee, 0x26305dfb, 0x3508a7a1, 0xc775c143, 0x2deaad10)}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Any}, ::Float64)
    @ Symbolics ./none:0
 [19] RuntimeGeneratedFunction
    @ ~/.julia/packages/RuntimeGeneratedFunctions/KrkGo/src/RuntimeGeneratedFunctions.jl:117 [inlined]
 [20] macro expansion
    @ ~/.julia/packages/Symbolics/J8IHJ/src/build_function.jl:367 [inlined]
 [21] macro expansion
    @ ~/.julia/packages/SymbolicUtils/qulQp/src/code.jl:351 [inlined]
 [22] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/KrkGo/src/RuntimeGeneratedFunctions.jl:129 [inlined]
 [23] macro expansion
    @ ./none:0 [inlined]
 [24] generated_callfunc(::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), Symbolics.var"#_RGF_ModTag", Symbolics.var"#_RGF_ModTag", (0x7c5f3805, 0xc53675ae, 0x37f82c0c, 0xcee36cba, 0xb23316b7)}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{Any}, ::Float64)
    @ Symbolics ./none:0
 [25] RuntimeGeneratedFunction
    @ ~/.julia/packages/RuntimeGeneratedFunctions/KrkGo/src/RuntimeGeneratedFunctions.jl:117 [inlined]
 [26] macro expansion
    @ ~/.julia/packages/Symbolics/J8IHJ/src/build_function.jl:367 [inlined]
 [27] macro expansion
    @ ~/.julia/packages/SymbolicUtils/qulQp/src/code.jl:351 [inlined]
 [28] macro expansion
    @ ~/.julia/packages/RuntimeGeneratedFunctions/KrkGo/src/RuntimeGeneratedFunctions.jl:129 [inlined]
 [29] macro expansion
    @ ./none:0 [inlined]
 [30] generated_callfunc
    @ ./none:0 [inlined]
 [31] RuntimeGeneratedFunction
    @ ~/.julia/packages/RuntimeGeneratedFunctions/KrkGo/src/RuntimeGeneratedFunctions.jl:117 [inlined]
 [32] (::ModelingToolkit.var"#f#464"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x4fbc4e5b, 0xd6540a24, 0x3047d7cc, 0x8bbe05e8, 0x19bb90d6)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf0c39447, 0x60f33fcc, 0xa84bcb8a, 0x3e197216, 0x49215fb8)}})(du::Vector{Float64}, u::Vector{Float64}, p::Vector{Any}, t::Float64)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/Dukjb/src/systems/diffeqs/abstractodesystem.jl:279
 [33] (::ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#464"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x4fbc4e5b, 0xd6540a24, 0x3047d7cc, 0x8bbe05e8, 0x19bb90d6)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf0c39447, 0x60f33fcc, 0xa84bcb8a, 0x3e197216, 0x49215fb8)}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, ModelingToolkit.var"#472#generated_observed#471"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem})(::Vector{Float64}, ::Vararg{Any})
    @ SciMLBase ~/.julia/packages/SciMLBase/bsPLM/src/scimlfunctions.jl:1704
 [34] top-level scope
    @ REPL[12]:1

What is really interesting:

For N = 30, i.e. a smaller grid, it works.

System

Julia Version 1.8.1
Commit afb6c60d69a (2022-09-06 15:09 UTC)
Platform Info:
  OS: macOS (arm64-apple-darwin21.5.0)
  CPU: 10 × Apple M1 Pro
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, apple-m1)
  Threads: 1 on 8 virtual cores

Please let me know, if you need more information.
Thanks a lot in advance for your help!

[xtalax]

Ok, since the defaults are an array of pairs of symbols to values, similar won't work this way, please extract the shape of u0 and use it with zeros for example.

Extremely weird that it works with small systems though.