Mismatch in Symbolic names, when using symbolic arrays as state variables

[bgctw]

Describe the bug 🐞

When using symbolic arrays as state variables, the name of the corresponding discretized state currently use several stacked var"". When trying to relate the discretized variables back to the original pde variables, this name is hard to reconstruct.
This leads to special cases and bugs, e.g. when trying to extract grid locations for a pde-variable as in the example below.

Expected behavior

Either avoid generating the names with several "var" or modify function _get_gridloc to deal with this case.
SciMLBase.ODEFunction(pdesys) should not throw an error.

Minimal Reproducible Example 👇

using OrdinaryDiffEq, ModelingToolkit, DomainSets
using MethodOfLines

@parameters t z 
z_m = +0.3 # maximum depth in m, change to positive depth
n_z = 16
z_grid = collect(range(0,z_m, length=n_z))
dz = z_m/(n_z-1) # control of same grid represented by a single number
discretization = MOLFiniteDifference([z => z_grid], t;
    advection_scheme = UpwindScheme(), approx_order = 2)

n_comp = 2
@parameters k_Y Y0 ω 
@variables Y(..)[1:n_comp] 
∂_t = Differential(t)
∂_z = Differential(z)
params = [
        k_Y => 2.0,
        Y0 => 200.0,
        ω => 0.01,
]

eqs_i = [[
    ∂_t(Y(t, z)[i]) ~ - k_Y * Y(t, z)[i] + -ω * ∂_z(Y(t, z)[i]) + 20*i,
] for i in 1:n_comp]
eqs = vcat(eqs_i...)

bcs_i = [[
    #Y(0, z) ~ Dz_exp(z, z_m, 2.0) * Y0, # initial exponential distribution with depth
    Y(0, z)[i] ~ Y0/z_m, #Dz_lin(z, z_m) * Y0, # initial constant with depth, modified by set problem.u0
    ω * Y(t, 0)[i] ~ 10.0*i, # specified flux at upper boundary
    ∂_z(Y(t, z_m)[i]) ~ 0, # negligible change in concentration at lower boundary
] for i in 1:n_comp]
bcs = vcat(bcs_i...)

# Space and time domains
domains = [t ∈ Interval(0.0, 500.0), z ∈ Interval(0.0, z_m), ]

# PDE system
state_vars = [Y(t, z)]
@named pdesys = PDESystem(eqs, bcs, domains, [t, z], state_vars, params)
# Convert the PDE problem into an ODE problem
prob = discretize(pdesys, discretization) 

analytic_u(p, t, x) = x / (t + p[1]) # dummy: not correct 
ODEFunction(pdesys, discretization; 
    analytic = vcat([Symbolics.scalarize(Y(t, z)) .=>  analytic_u]...))

# from _get_gridloc
# (ut, u, operation(s.ū[1])) = (var"var\"Y_Any[1]\""(t), var"var\"Y_Any[1]\"", var"Y_Any[1]")

Error & Stacktrace ⚠️

Discretization failed, please post an issue on https://github.com/SciML/MethodOfLines.jl with the failing code and system at low point count.

ERROR: ArgumentError: invalid index: nothing of type Nothing
Stacktrace:
  [1] to_index(i::Nothing)
    @ Base ./indices.jl:300
  [2] to_index(A::Vector{Any}, i::Nothing)
    @ Base ./indices.jl:277
  [3] _to_indices1(A::Vector{Any}, inds::Tuple{Base.OneTo{Int64}}, I1::Nothing)
    @ Base ./indices.jl:359
  [4] to_indices
    @ ./indices.jl:354 [inlined]
  [5] to_indices
    @ ./indices.jl:345 [inlined]
  [6] getindex(A::Vector{Any}, I::Nothing)
    @ Base ./abstractarray.jl:1294
  [7] _get_gridloc(s::MethodOfLines.DiscreteSpace{1, 2, MethodOfLines.CenterAlignedGrid}, ut::SymbolicUtils.BasicSymbolic{Vector{Real}}, is::Int64)
    @ MethodOfLines /Net/Groups/BGI/people/twutz/julia/dev/MethodOfLines/src/MOL_utils.jl:24
  [8] get_gridloc(u::SymbolicUtils.BasicSymbolic{Real}, s::MethodOfLines.DiscreteSpace{1, 2, MethodOfLines.CenterAlignedGrid})
    @ MethodOfLines /Net/Groups/BGI/people/twutz/julia/dev/MethodOfLines/src/MOL_utils.jl:35
  [9] _broadcast_getindex_evalf
    @ ./broadcast.jl:683 [inlined]
 [10] _broadcast_getindex
    @ ./broadcast.jl:656 [inlined]
 [11] getindex
    @ ./broadcast.jl:610 [inlined]
 [12] copy(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Tuple{Base.OneTo{Int64}}, typeof(MethodOfLines.get_gridloc), Tuple{Vector{Any}, Tuple{MethodOfLines.DiscreteSpace{1, 2, MethodOfLines.CenterAlignedGrid}}}})
    @ Base.Broadcast ./broadcast.jl:912
 [13] materialize
    @ ./broadcast.jl:873 [inlined]
 [14] ODEFunction(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}; analytic::Vector{Pair{Num, typeof(analytic_u)}}, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ MethodOfLines /Net/Groups/BGI/people/twutz/julia/dev/MethodOfLines/src/MOL_discretization.jl:97
 [15] top-level scope
    @ REPL[81]:1

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

Project MethodOfLines v0.10.2
Status `/Net/Groups/BGI/people/twutz/julia/dev/MethodOfLines/Project.toml`
  [861a8166] Combinatorics v1.0.2
  [2b5f629d] DiffEqBase v6.141.0
⌅ [5b8099bc] DomainSets v0.6.7
  [615f187c] IfElse v0.1.1
⌅ [a98d9a8b] Interpolations v0.14.7
  [23fbe1c1] Latexify v0.16.1
  [961ee093] ModelingToolkit v8.73.1
  [1dea7af3] OrdinaryDiffEq v6.59.3
  [a7812802] PDEBase v0.1.6
  [aea7be01] PrecompileTools v1.2.0
  [7e49a35a] RuntimeGeneratedFunctions v0.5.12
  [0bca4576] SciMLBase v2.9.1
  [90137ffa] StaticArrays v1.7.0
  [d1185830] SymbolicUtils v1.4.0
  [0c5d862f] Symbolics v5.10.0
  [8ea1fca8] TermInterface v0.3.3
  [37e2e46d] LinearAlgebra

• Output of using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)

Project MethodOfLines v0.10.2
Status `/Net/Groups/BGI/people/twutz/julia/dev/MethodOfLines/Project.toml`
  [861a8166] Combinatorics v1.0.2
  [2b5f629d] DiffEqBase v6.141.0
⌅ [5b8099bc] DomainSets v0.6.7
  [615f187c] IfElse v0.1.1
⌅ [a98d9a8b] Interpolations v0.14.7
  [23fbe1c1] Latexify v0.16.1
  [961ee093] ModelingToolkit v8.73.1
  [1dea7af3] OrdinaryDiffEq v6.59.3
  [a7812802] PDEBase v0.1.6
  [aea7be01] PrecompileTools v1.2.0
  [7e49a35a] RuntimeGeneratedFunctions v0.5.12
  [0bca4576] SciMLBase v2.9.1
  [90137ffa] StaticArrays v1.7.0
  [d1185830] SymbolicUtils v1.4.0
  [0c5d862f] Symbolics v5.10.0
  [8ea1fca8] TermInterface v0.3.3
  [37e2e46d] LinearAlgebra
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated`

julia> Pkg.status(; mode = PKGMODE_MANIFEST)
Project MethodOfLines v0.10.2
Status `/Net/Groups/BGI/people/twutz/julia/dev/MethodOfLines/Manifest.toml`
  [47edcb42] ADTypes v0.2.5
⌅ [c3fe647b] AbstractAlgebra v0.32.5
  [1520ce14] AbstractTrees v0.4.4
  [7d9f7c33] Accessors v0.1.33
  [79e6a3ab] Adapt v3.7.1
  [ec485272] ArnoldiMethod v0.2.0
  [4fba245c] ArrayInterface v7.6.1
  [13072b0f] AxisAlgorithms v1.0.1
  [e2ed5e7c] Bijections v0.1.6
  [62783981] BitTwiddlingConvenienceFunctions v0.1.5
  [2a0fbf3d] CPUSummary v0.2.4
  [00ebfdb7] CSTParser v3.3.6
  [49dc2e85] Calculus v0.5.1
  [d360d2e6] ChainRulesCore v1.18.0
  [fb6a15b2] CloseOpenIntervals v0.1.12
  [861a8166] Combinatorics v1.0.2
  [a80b9123] CommonMark v0.8.12
  [38540f10] CommonSolve v0.2.4
  [bbf7d656] CommonSubexpressions v0.3.0
  [34da2185] Compat v4.10.0
  [b152e2b5] CompositeTypes v0.1.3
  [a33af91c] CompositionsBase v0.1.2
  [2569d6c7] ConcreteStructs v0.2.3
  [187b0558] ConstructionBase v1.5.4
  [adafc99b] CpuId v0.3.1
  [a8cc5b0e] Crayons v4.1.1
  [9a962f9c] DataAPI v1.15.0
  [864edb3b] DataStructures v0.18.15
  [e2d170a0] DataValueInterfaces v1.0.0
  [8bb1440f] DelimitedFiles v1.9.1
  [2b5f629d] DiffEqBase v6.141.0
  [459566f4] DiffEqCallbacks v2.34.0
  [163ba53b] DiffResults v1.1.0
  [b552c78f] DiffRules v1.15.1
  [b4f34e82] Distances v0.10.10
  [31c24e10] Distributions v0.25.103
  [ffbed154] DocStringExtensions v0.9.3
⌅ [5b8099bc] DomainSets v0.6.7
  [fa6b7ba4] DualNumbers v0.6.8
  [7c1d4256] DynamicPolynomials v0.5.3
  [4e289a0a] EnumX v1.0.4
  [f151be2c] EnzymeCore v0.6.4
  [d4d017d3] ExponentialUtilities v1.25.0
  [e2ba6199] ExprTools v0.1.10
  [7034ab61] FastBroadcast v0.2.8
  [9aa1b823] FastClosures v0.3.2
  [29a986be] FastLapackInterface v2.0.0
  [1a297f60] FillArrays v1.8.0
  [6a86dc24] FiniteDiff v2.21.1
  [59287772] Formatting v0.4.2
  [f6369f11] ForwardDiff v0.10.36
  [069b7b12] FunctionWrappers v1.1.3
  [77dc65aa] FunctionWrappersWrappers v0.1.3
  [d9f16b24] Functors v0.4.5
  [46192b85] GPUArraysCore v0.1.5
  [c145ed77] GenericSchur v0.5.3
  [c27321d9] Glob v1.3.1
  [86223c79] Graphs v1.9.0
⌅ [0b43b601] Groebner v0.4.4
  [d5909c97] GroupsCore v0.4.0
  [3e5b6fbb] HostCPUFeatures v0.1.16
  [34004b35] HypergeometricFunctions v0.3.23
  [615f187c] IfElse v0.1.1
  [d25df0c9] Inflate v0.1.4
  [18e54dd8] IntegerMathUtils v0.1.2
⌅ [a98d9a8b] Interpolations v0.14.7
  [8197267c] IntervalSets v0.7.8
  [3587e190] InverseFunctions v0.1.12
  [92d709cd] IrrationalConstants v0.2.2
  [82899510] IteratorInterfaceExtensions v1.0.0
  [692b3bcd] JLLWrappers v1.5.0
  [682c06a0] JSON v0.21.4
  [98e50ef6] JuliaFormatter v1.0.43
  [ccbc3e58] JumpProcesses v9.8.0
  [ef3ab10e] KLU v0.4.1
  [ba0b0d4f] Krylov v0.9.4
  [b964fa9f] LaTeXStrings v1.3.1
  [2ee39098] LabelledArrays v1.14.0
  [984bce1d] LambertW v0.4.6
  [23fbe1c1] Latexify v0.16.1
  [73f95e8e] LatticeRules v0.0.1
  [10f19ff3] LayoutPointers v0.1.15
  [50d2b5c4] Lazy v0.15.1
  [d3d80556] LineSearches v7.2.0
  [7ed4a6bd] LinearSolve v2.20.0
  [2ab3a3ac] LogExpFunctions v0.3.26
  [bdcacae8] LoopVectorization v0.12.166
  [d8e11817] MLStyle v0.4.17
  [1914dd2f] MacroTools v0.5.11
  [d125e4d3] ManualMemory v0.1.8
  [e1d29d7a] Missings v1.1.0
  [961ee093] ModelingToolkit v8.73.1
  [46d2c3a1] MuladdMacro v0.2.4
  [102ac46a] MultivariatePolynomials v0.5.3
  [d8a4904e] MutableArithmetics v1.3.3
  [d41bc354] NLSolversBase v7.8.3
  [2774e3e8] NLsolve v4.5.1
  [77ba4419] NaNMath v1.0.2
  [8913a72c] NonlinearSolve v2.8.2
  [6fe1bfb0] OffsetArrays v1.12.10
  [bac558e1] OrderedCollections v1.6.3
  [1dea7af3] OrdinaryDiffEq v6.59.3
  [a7812802] PDEBase v0.1.6
  [90014a1f] PDMats v0.11.30
  [65ce6f38] PackageExtensionCompat v1.0.2
  [d96e819e] Parameters v0.12.3
  [69de0a69] Parsers v2.8.0
  [e409e4f3] PoissonRandom v0.4.4
  [f517fe37] Polyester v0.7.9
  [1d0040c9] PolyesterWeave v0.2.1
  [d236fae5] PreallocationTools v0.4.12
  [aea7be01] PrecompileTools v1.2.0
  [21216c6a] Preferences v1.4.1
  [27ebfcd6] Primes v0.5.5
  [1fd47b50] QuadGK v2.9.1
  [8a4e6c94] QuasiMonteCarlo v0.3.3
  [fb686558] RandomExtensions v0.4.4
  [e6cf234a] RandomNumbers v1.5.3
  [c84ed2f1] Ratios v0.4.5
  [3cdcf5f2] RecipesBase v1.3.4
  [731186ca] RecursiveArrayTools v2.38.10
  [f2c3362d] RecursiveFactorization v0.2.21
  [189a3867] Reexport v1.2.2
  [ae029012] Requires v1.3.0
  [79098fc4] Rmath v0.7.1
  [7e49a35a] RuntimeGeneratedFunctions v0.5.12
  [fdea26ae] SIMD v3.4.6
  [94e857df] SIMDTypes v0.1.0
  [476501e8] SLEEFPirates v0.6.42
  [0bca4576] SciMLBase v2.9.1
  [e9a6253c] SciMLNLSolve v0.1.9
  [c0aeaf25] SciMLOperators v0.3.7
  [efcf1570] Setfield v1.1.1
  [727e6d20] SimpleNonlinearSolve v0.1.25
  [699a6c99] SimpleTraits v0.9.4
  [ce78b400] SimpleUnPack v1.1.0
  [66db9d55] SnoopPrecompile v1.0.3
  [ed01d8cd] Sobol v1.5.0
  [a2af1166] SortingAlgorithms v1.2.0
  [47a9eef4] SparseDiffTools v2.13.0
  [e56a9233] Sparspak v0.3.9
  [276daf66] SpecialFunctions v2.3.1
  [aedffcd0] Static v0.8.8
  [0d7ed370] StaticArrayInterface v1.4.1
  [90137ffa] StaticArrays v1.7.0
  [1e83bf80] StaticArraysCore v1.4.2
  [82ae8749] StatsAPI v1.7.0
  [2913bbd2] StatsBase v0.34.2
  [4c63d2b9] StatsFuns v1.3.0
  [7792a7ef] StrideArraysCore v0.5.2
  [2efcf032] SymbolicIndexingInterface v0.2.2
  [d1185830] SymbolicUtils v1.4.0
  [0c5d862f] Symbolics v5.10.0
  [3783bdb8] TableTraits v1.0.1
  [bd369af6] Tables v1.11.1
  [8ea1fca8] TermInterface v0.3.3
  [8290d209] ThreadingUtilities v0.5.2
  [a759f4b9] TimerOutputs v0.5.23
  [0796e94c] Tokenize v0.5.26
  [a2a6695c] TreeViews v0.3.0
  [d5829a12] TriangularSolve v0.1.20
  [410a4b4d] Tricks v0.1.8
  [781d530d] TruncatedStacktraces v1.4.0
  [5c2747f8] URIs v1.5.1
  [3a884ed6] UnPack v1.0.2
  [1986cc42] Unitful v1.19.0
  [a7c27f48] Unityper v0.1.5
  [3d5dd08c] VectorizationBase v0.21.65
  [19fa3120] VertexSafeGraphs v0.2.0
  [efce3f68] WoodburyMatrices v0.5.6
  [1d5cc7b8] IntelOpenMP_jll v2023.2.0+0
  [856f044c] MKL_jll v2023.2.0+0
  [efe28fd5] OpenSpecFun_jll v0.5.5+0
  [f50d1b31] Rmath_jll v0.4.0+0
  [0dad84c5] ArgTools v1.1.1
  [56f22d72] Artifacts
  [2a0f44e3] Base64
  [ade2ca70] Dates
  [8ba89e20] Distributed
  [f43a241f] Downloads v1.6.0
  [7b1f6079] FileWatching
  [9fa8497b] Future
  [b77e0a4c] InteractiveUtils
  [4af54fe1] LazyArtifacts
  [b27032c2] LibCURL v0.6.3
  [76f85450] LibGit2
  [8f399da3] Libdl
  [37e2e46d] LinearAlgebra
  [56ddb016] Logging
  [d6f4376e] Markdown
  [a63ad114] Mmap
  [ca575930] NetworkOptions v1.2.0
  [44cfe95a] Pkg v1.9.2
  [de0858da] Printf
  [3fa0cd96] REPL
  [9a3f8284] Random
  [ea8e919c] SHA v0.7.0
  [9e88b42a] Serialization
  [1a1011a3] SharedArrays
  [6462fe0b] Sockets
  [2f01184e] SparseArrays
  [10745b16] Statistics v1.9.0
  [4607b0f0] SuiteSparse
  [fa267f1f] TOML v1.0.3
  [a4e569a6] Tar v1.10.0
  [8dfed614] Test
  [cf7118a7] UUIDs
  [4ec0a83e] Unicode
  [e66e0078] CompilerSupportLibraries_jll v1.0.5+0
  [deac9b47] LibCURL_jll v7.84.0+0
  [29816b5a] LibSSH2_jll v1.10.2+0
  [c8ffd9c3] MbedTLS_jll v2.28.2+0
  [14a3606d] MozillaCACerts_jll v2022.10.11
  [4536629a] OpenBLAS_jll v0.3.21+4
  [05823500] OpenLibm_jll v0.8.1+0
  [bea87d4a] SuiteSparse_jll v5.10.1+6
  [83775a58] Zlib_jll v1.2.13+0
  [8e850b90] libblastrampoline_jll v5.8.0+0
  [8e850ede] nghttp2_jll v1.48.0+0
  [3f19e933] p7zip_jll v17.4.0+0
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated -m`

• Output of versioninfo()

Julia Version 1.9.2
Commit e4ee485e909 (2023-07-05 09:39 UTC)
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 64 × Intel(R) Xeon(R) Silver 4216 CPU @ 2.10GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-14.0.6 (ORCJIT, cascadelake)
  Threads: 8 on 64 virtual cores
Environment:
  LD_LIBRARY_PATH = /Net/Groups/Services/HPC_22/apps/julia/julia-1.9.2/lib
  JULIA_PKG_DEVDIR = /User/homes/twutz/julia/dev
  JULIA_DEPOT_PATH = /User/homes/twutz/scratch/twutz/julia_cluster_depots
  JULIA_PKG_USE_CLI_GIT = true
  JULIA_NUM_THREADS = 8

Additional context

We use symbolic arrays as state variables, so that we can concisely generate similar equations for several states.
We also want to use non-equidistant spacial 1-d grid.
For modifying the problem state (reapplying the boundary condition function with different parameters that are optimized) I need to get the grid locations from the system or problem.
This works for scalar symbolic states, but fails for the symbolic arrays.

This is not urgent. I just want to support development and testing with this MWE.

[xtalax]

You can try scalarizing yourself before solving, we try to do this internally to make systems compatible. please be aware that ODEFunction is intended to be internal for MOL, you should try getting this with prob.f instead.

Non equidistant grid is supported, see: https://docs.sciml.ai/MethodOfLines/stable/nonuniform/. you can also write[x => xpoints] where xpoints is a vector of grid points.

you can use

grid = get_discrete(pdesys, discretization)
discrete_x = grid[x]

to get the grid points, nb it is recommended to use remake to update system parameters to avoid re-discretizing each time

I will look at this bug in more detail soon