Bound error on Discretization when specifying a uniform grid as a vector

[bgctw]

Describe the bug 🐞

When I provide a uniform grid as vector to MOLFiniteDifference, the discretization
of a pde results in bounds errors.
Contrary, the discretization works, if the same grid is specified by a single delta-number and the domain.

Expected behavior

I expect providing the grid as a vector would yield the same results as providing the grid by a delta-number. In the future, I want to specify a non-uniform grid with closer spacing near the top of the depth domain.

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 = grid_exp(n_z, z_m, 3)
z_grid = collect(range(0,z_m, length=n_z))
# z_grid = 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)
    #advection_scheme=WENOScheme(), approx_order = 2)
dzs = diff(z_grid)
dzsl = vcat(dzs[1]/2, dzs[2:end], dzs[end]/2) # assume first and last layer only half 

@parameters k_Y Y0 i_Y ω i_Y_agr i_Y_agr_pulse
@variables Y(..) i_Yo(..) adv_Yo(..) dec_Y(..) Y_t(..) Y_zz(..) Yi(..) i_Yi(..) Y_z(..) adv_Yi(..) 
∂_t = Differential(t)
∂_z = Differential(z)
params = [
        k_Y => 2.0,
        Y0 => 200.0,
        i_Y => 50.0,
        ω => 0.01,
        i_Y_agr => 10.0,
        i_Y_agr_pulse => 0.0,        
]


i_Yz(t,z,i_Y) = 20.0 #Dz_exp(z, z_m,4.0) * i_Y  # inputs that integrate to i_Y independent of time
@register_symbolic i_Yz(t, z, i_Y)

z_min = 0.0  # directly using 0.0 in Integral causes error in solution wrapping
Iz = Integral(z in DomainSets.ClosedInterval(z_min, z))


eqs = [
    ∂_t(Y(t, z)) ~ Y_t(t, z),
    Y_t(t, z) ~ i_Yo(t, z) - dec_Y(t, z) + adv_Yo(t, z),
    i_Yo(t, z) ~ i_Yz(t, z, i_Y),         # observable input Y
    dec_Y(t, z) ~ k_Y * Y(t, z),          # observable of decomposition 
    adv_Yo(t, z) ~ -ω * ∂_z(Y(t, z)),  # observable advective flux of Y
    # further observables that are not used in eqs
    Y_z(t, z) ~ ∂_z(Y(t, z)),
    #Y_zz(t,z) ~ ∂_zz(Y(t, z)),
    Yi(t, z) ~ Iz(Y(t, z)),
    i_Yi(t, z) ~ Iz(i_Yo(t, z)), # integral over litter inputs to check
    adv_Yi(t,z) ~ Iz(adv_Yo(t, z)), # integral over advection inputs
]

fgauss(x, μ, σ2) = 1/sqrt(2*pi*σ2) * exp(-(x-μ)^2/(2*σ2))
fagr(t, i_Y_agr, i_Y_agr_pulse) = i_Y_agr + i_Y_agr_pulse*fgauss(t, 80, 2^2)
@register_symbolic fagr(t, i_Y_agr, i_Y_agr_pulse)

bcs = [
    Y(0, z) ~ 1/z_m * Y0, 
    ∂_z(Y(t, 0)) ~ (Y(t, 0) - fagr(t,i_Y_agr, i_Y_agr_pulse)/ω) / dzs[1], 
    ∂_z(Y(t, z_m)) ~ 0.0, 
    ]

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

# PDE system
state_vars = [Y(t, z), i_Yo(t,z), adv_Yo(t,z), dec_Y(t,z), Y_t(t,z), Yi(t,z), i_Yi(t,z), Y_z(t,z), adv_Yi(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) 

Error & Stacktrace ⚠️

ERROR: BoundsError: attempt to access 16-element Vector{SymbolicUtils.BasicSymbolic{Real}} at index [CartesianIndex{1}[CartesianIndex(0,), CartesianIndex(1,)]]
Stacktrace:
  [1] throw_boundserror(A::Vector{SymbolicUtils.BasicSymbolic{Real}}, I::Tuple{Vector{CartesianIndex{1}}})
    @ Base ./abstractarray.jl:744
  [2] checkbounds
    @ ./abstractarray.jl:709 [inlined]
  [3] _getindex
    @ ./multidimensional.jl:860 [inlined]
  [4] getindex(A::Vector{SymbolicUtils.BasicSymbolic{Real}}, I::Vector{CartesianIndex{1}})
    @ Base ./abstractarray.jl:1294
  [5] (::MethodOfLines.var"#wind_ufunc#279"{MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}})(v::SymbolicUtils.BasicSymbolic{Real}, I::Vector{CartesianIndex{1}}, x::SymbolicUtils.BasicSymbolic{Real})
    @ MethodOfLines ~/julia/dev/MethodOfLines/src/discretization/schemes/upwind_difference/upwind_difference.jl:83
  [6] upwind_difference(d::Int64, II::CartesianIndex{1}, s::MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}, bs::Vector{Any}, derivweights::MethodOfLines.DifferentialDiscretizer{Pair{SymbolicUtils.BasicSymbolic{Real}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Vector{Float64}, Vector{StaticArraysCore.SVector{3, Float64}}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}, jx::Tuple{Int64, SymbolicUtils.BasicSymbolic{Real}}, u::SymbolicUtils.BasicSymbolic{Real}, ufunc::MethodOfLines.var"#wind_ufunc#279"{MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}}, ispositive::Bool)
    @ MethodOfLines ~/julia/dev/MethodOfLines/src/discretization/schemes/upwind_difference/upwind_difference.jl:71
  [7] (::MethodOfLines.var"#278#290"{SymbolicUtils.BasicSymbolic{Real}, Int64, SymbolicUtils.BasicSymbolic{Real}, CartesianIndex{1}, MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}, MethodOfLines.DifferentialDiscretizer{Pair{SymbolicUtils.BasicSymbolic{Real}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Vector{Float64}, Vector{StaticArraysCore.SVector{3, Float64}}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}, Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}, MethodOfLines.var"#wind_ufunc#279"{MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}}})(d::Int64)
    @ MethodOfLines ~/julia/dev/MethodOfLines/src/discretization/schemes/upwind_difference/upwind_difference.jl:119
  [8] iterate
    @ ./generator.jl:47 [inlined]
  [9] _collect(c::Vector{Any}, itr::Base.Generator{Vector{Any}, MethodOfLines.var"#278#290"{SymbolicUtils.BasicSymbolic{Real}, Int64, SymbolicUtils.BasicSymbolic{Real}, CartesianIndex{1}, MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}, MethodOfLines.DifferentialDiscretizer{Pair{SymbolicUtils.BasicSymbolic{Real}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Vector{Float64}, Vector{StaticArraysCore.SVector{3, Float64}}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}, Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}, MethodOfLines.var"#wind_ufunc#279"{MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}}}}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})
    @ Base ./array.jl:802
 [10] collect_similar(cont::Vector{Any}, itr::Base.Generator{Vector{Any}, MethodOfLines.var"#278#290"{SymbolicUtils.BasicSymbolic{Real}, Int64, SymbolicUtils.BasicSymbolic{Real}, CartesianIndex{1}, MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}, MethodOfLines.DifferentialDiscretizer{Pair{SymbolicUtils.BasicSymbolic{Real}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Vector{Float64}, Vector{StaticArraysCore.SVector{3, Float64}}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}, Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}, MethodOfLines.var"#wind_ufunc#279"{MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}}}})
    @ Base ./array.jl:711
 [11] map(f::Function, A::Vector{Any})
    @ Base ./abstractarray.jl:3261
 [12] (::MethodOfLines.var"#277#289"{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}, CartesianIndex{1}, MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}, MethodOfLines.DifferentialDiscretizer{Pair{SymbolicUtils.BasicSymbolic{Real}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Vector{Float64}, Vector{StaticArraysCore.SVector{3, Float64}}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}, Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}, MethodOfLines.var"#wind_ufunc#279"{MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}}})(x::SymbolicUtils.BasicSymbolic{Real})
    @ MethodOfLines ~/julia/dev/MethodOfLines/src/discretization/schemes/upwind_difference/upwind_difference.jl:118
 [13] MappingRF
    @ ./reduce.jl:95 [inlined]
 [14] _foldl_impl(op::Base.MappingRF{MethodOfLines.var"#277#289"{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}, CartesianIndex{1}, MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}, MethodOfLines.DifferentialDiscretizer{Pair{SymbolicUtils.BasicSymbolic{Real}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Vector{Float64}, Vector{StaticArraysCore.SVector{3, Float64}}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}, Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}, MethodOfLines.var"#wind_ufunc#279"{MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}}}, Base.BottomRF{typeof(MethodOfLines.safe_vcat)}}, init::Vector{Any}, itr::Vector{Any})
    @ Base ./reduce.jl:58
 [15] foldl_impl
    @ ./reduce.jl:48 [inlined]
 [16] mapfoldl_impl(f::MethodOfLines.var"#277#289"{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}, CartesianIndex{1}, MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}, MethodOfLines.DifferentialDiscretizer{Pair{SymbolicUtils.BasicSymbolic{Real}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Vector{Float64}, Vector{StaticArraysCore.SVector{3, Float64}}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}, Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}, MethodOfLines.var"#wind_ufunc#279"{MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}}}, op::typeof(MethodOfLines.safe_vcat), nt::Vector{Any}, itr::Vector{Any})
    @ Base ./reduce.jl:44
 [17] _mapreduce_dim(f::Function, op::Function, nt::Vector{Any}, A::Vector{Any}, #unused#::Colon)
    @ Base ./reducedim.jl:362
 [18] #mapreduce#800
    @ ./reducedim.jl:357 [inlined]
 [19] (::MethodOfLines.var"#276#288"{Vector{Any}, CartesianIndex{1}, MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}, MethodOfLines.DifferentialDiscretizer{Pair{SymbolicUtils.BasicSymbolic{Real}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Vector{Float64}, Vector{StaticArraysCore.SVector{3, Float64}}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}, Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}, MethodOfLines.var"#wind_ufunc#279"{MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}}})(u::SymbolicUtils.BasicSymbolic{Real})
    @ MethodOfLines ~/julia/dev/MethodOfLines/src/discretization/schemes/upwind_difference/upwind_difference.jl:112
 [20] MappingRF
    @ ./reduce.jl:95 [inlined]
 [21] _foldl_impl(op::Base.MappingRF{MethodOfLines.var"#276#288"{Vector{Any}, CartesianIndex{1}, MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}, MethodOfLines.DifferentialDiscretizer{Pair{SymbolicUtils.BasicSymbolic{Real}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Vector{Float64}, Vector{StaticArraysCore.SVector{3, Float64}}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}, Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}, MethodOfLines.var"#wind_ufunc#279"{MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}}}, Base.BottomRF{typeof(MethodOfLines.safe_vcat)}}, init::Vector{Any}, itr::Vector{Any})
    @ Base ./reduce.jl:58
 [22] foldl_impl
    @ ./reduce.jl:48 [inlined]
 [23] mapfoldl_impl
    @ ./reduce.jl:44 [inlined]
 [24] _mapreduce_dim
    @ ./reducedim.jl:362 [inlined]
 [25] #mapreduce#800
    @ ./reducedim.jl:357 [inlined]
 [26] mapreduce
    @ ./reducedim.jl:357 [inlined]
 [27] generate_winding_rules(II::CartesianIndex{1}, s::MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}, depvars::Vector{Any}, derivweights::MethodOfLines.DifferentialDiscretizer{Pair{SymbolicUtils.BasicSymbolic{Real}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Vector{Float64}, Vector{StaticArraysCore.SVector{3, Float64}}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}, bcmap::Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, indexmap::Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}, terms::Vector{Any}; skip::Vector{Any})
    @ MethodOfLines ~/julia/dev/MethodOfLines/src/discretization/schemes/upwind_difference/upwind_difference.jl:111
 [28] generate_winding_rules
    @ ~/julia/dev/MethodOfLines/src/discretization/schemes/upwind_difference/upwind_difference.jl:82 [inlined]
 [29] generate_finite_difference_rules(II::CartesianIndex{1}, s::MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}, depvars::Vector{Any}, pde::Equation, derivweights::MethodOfLines.DifferentialDiscretizer{Pair{SymbolicUtils.BasicSymbolic{Real}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Vector{Float64}, Vector{StaticArraysCore.SVector{3, Float64}}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}, bmap::Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, indexmap::Dict{SymbolicUtils.BasicSymbolic{Real}, Int64})
    @ MethodOfLines ~/julia/dev/MethodOfLines/src/discretization/generate_finite_difference_rules.jl:57
 [30] discretize_equation_at_point(II::CartesianIndex{1}, s::MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}, depvars::Vector{Any}, pde::Equation, derivweights::MethodOfLines.DifferentialDiscretizer{Pair{SymbolicUtils.BasicSymbolic{Real}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Vector{Float64}, Vector{StaticArraysCore.SVector{3, Float64}}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}, bcmap::Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, eqvar::SymbolicUtils.BasicSymbolic{Real}, indexmap::Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}, boundaryvalfuncs::Vector{Any})
    @ MethodOfLines ~/julia/dev/MethodOfLines/src/scalar_discretization.jl:34
 [31] #470
    @ ~/julia/dev/MethodOfLines/src/scalar_discretization.jl:24 [inlined]
 [32] iterate
    @ ./generator.jl:47 [inlined]
 [33] _collect(c::CartesianIndices{1, Tuple{UnitRange{Int64}}}, itr::Base.Generator{CartesianIndices{1, Tuple{UnitRange{Int64}}}, MethodOfLines.var"#470#472"{Equation, SymbolicUtils.BasicSymbolic{Real}, Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, Vector{Any}, MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}, MethodOfLines.DifferentialDiscretizer{Pair{SymbolicUtils.BasicSymbolic{Real}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Vector{Float64}, Vector{StaticArraysCore.SVector{3, Float64}}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}, Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}, Vector{Any}}}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})
    @ Base ./array.jl:802
 [34] collect_similar(cont::CartesianIndices{1, Tuple{UnitRange{Int64}}}, itr::Base.Generator{CartesianIndices{1, Tuple{UnitRange{Int64}}}, MethodOfLines.var"#470#472"{Equation, SymbolicUtils.BasicSymbolic{Real}, Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, Vector{Any}, MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}, MethodOfLines.DifferentialDiscretizer{Pair{SymbolicUtils.BasicSymbolic{Real}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Vector{Float64}, Vector{StaticArraysCore.SVector{3, Float64}}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}, Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}, Vector{Any}}})
    @ Base ./array.jl:711
 [35] map(f::Function, A::CartesianIndices{1, Tuple{UnitRange{Int64}}})
    @ Base ./abstractarray.jl:3261
 [36] discretize_equation!(disc_state::PDEBase.EquationState, pde::Equation, interiormap::MethodOfLines.InteriorMap, eqvar::SymbolicUtils.BasicSymbolic{Real}, bcmap::Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, depvars::Vector{Any}, s::MethodOfLines.DiscreteSpace{1, 9, MethodOfLines.CenterAlignedGrid}, derivweights::MethodOfLines.DifferentialDiscretizer{Pair{SymbolicUtils.BasicSymbolic{Real}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Vector{Float64}, Vector{StaticArraysCore.SVector{3, Float64}}, S2, S3, Nothing, Nothing} where {S2, S3}}, UpwindScheme}, indexmap::Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
    @ MethodOfLines ~/julia/dev/MethodOfLines/src/scalar_discretization.jl:23
 [37] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
    @ PDEBase ~/scratch/twutz/julia_cluster_depots/packages/PDEBase/jX5Yp/src/symbolic_discretize.jl:79
 [38] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}; analytic::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ PDEBase ~/scratch/twutz/julia_cluster_depots/packages/PDEBase/jX5Yp/src/discretization_state.jl:58
 [39] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
    @ PDEBase ~/scratch/twutz/julia_cluster_depots/packages/PDEBase/jX5Yp/src/discretization_state.jl:55
 [40] top-level scope
    @ REPL[46]:1

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

Status `/Net/Groups/BGI/people/twutz/julia/dev/MTKHelpers/Project.toml`
⌃ [b0b7db55] ComponentArrays v0.15.0
⌃ [2b5f629d] DiffEqBase v6.128.4
⌃ [31c24e10] Distributions v0.25.100
⌅ [5b8099bc] DomainSets v0.6.7
  [842dd82b] InlineStrings v1.4.0
  [2ee39098] LabelledArrays v1.14.0
  [94925ecb] MethodOfLines v0.10.2 `~/julia/dev/MethodOfLines`
⌃ [961ee093] ModelingToolkit v8.69.1
  [86f7a689] NamedArrays v0.10.0
⌃ [1dea7af3] OrdinaryDiffEq v6.58.2
  [731186ca] RecursiveArrayTools v2.38.10
  [ae029012] Requires v1.3.0
⌅ [0bca4576] SciMLBase v1.98.1
⌃ [90137ffa] StaticArrays v1.6.3
  [d1185830] SymbolicUtils v1.4.0

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

Status `/Net/Groups/BGI/people/twutz/julia/dev/MTKHelpers/Manifest.toml`
  [47edcb42] ADTypes v0.2.5
⌅ [c3fe647b] AbstractAlgebra v0.32.5
  [1520ce14] AbstractTrees v0.4.4
  [79e6a3ab] Adapt v3.7.1
  [ec485272] ArnoldiMethod v0.2.0
  [4fba245c] ArrayInterface v7.5.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
⌃ [b0b7db55] ComponentArrays v0.15.0
  [b152e2b5] CompositeTypes v0.1.3
  [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.128.4
  [459566f4] DiffEqCallbacks v2.34.0
  [163ba53b] DiffResults v1.1.0
  [b552c78f] DiffRules v1.15.1
  [b4f34e82] Distances v0.10.10
⌃ [31c24e10] Distributions v0.25.100
  [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.3
  [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.7.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
  [842dd82b] InlineStrings v1.4.0
  [18e54dd8] IntegerMathUtils v0.1.2
  [a98d9a8b] Interpolations v0.14.7
  [8197267c] IntervalSets v0.7.8
  [41ab1584] InvertedIndices v1.3.0
  [92d709cd] IrrationalConstants v0.2.2
  [82899510] IteratorInterfaceExtensions v1.0.0
  [692b3bcd] JLLWrappers v1.5.0
  [682c06a0] JSON v0.21.4
  [98e50ef6] JuliaFormatter v1.0.42
  [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
  [10f19ff3] LayoutPointers v0.1.15
  [50d2b5c4] Lazy v0.15.1
  [d3d80556] LineSearches v7.2.0
⌃ [7ed4a6bd] LinearSolve v2.16.2
  [2ab3a3ac] LogExpFunctions v0.3.26
  [bdcacae8] LoopVectorization v0.12.166
  [33e6dc65] MKL v0.6.1
  [d8e11817] MLStyle v0.4.17
  [1914dd2f] MacroTools v0.5.11
  [d125e4d3] ManualMemory v0.1.8
  [94925ecb] MethodOfLines v0.10.2 `~/julia/dev/MethodOfLines`
  [e1d29d7a] Missings v1.1.0
⌃ [961ee093] ModelingToolkit v8.69.1
  [46d2c3a1] MuladdMacro v0.2.4
  [102ac46a] MultivariatePolynomials v0.5.2
  [d8a4904e] MutableArithmetics v1.3.3
  [d41bc354] NLSolversBase v7.8.3
  [2774e3e8] NLsolve v4.5.1
  [77ba4419] NaNMath v1.0.2
  [86f7a689] NamedArrays v0.10.0
⌃ [8913a72c] NonlinearSolve v1.10.1
  [6fe1bfb0] OffsetArrays v1.12.10
  [bac558e1] OrderedCollections v1.6.2
⌃ [1dea7af3] OrdinaryDiffEq v6.58.2
  [a7812802] PDEBase v0.1.6
  [90014a1f] PDMats v0.11.29
  [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
  [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 v1.98.1
  [e9a6253c] SciMLNLSolve v0.1.9
  [c0aeaf25] SciMLOperators v0.3.7
  [efcf1570] Setfield v1.1.1
⌃ [727e6d20] SimpleNonlinearSolve v0.1.23
  [699a6c99] SimpleTraits v0.9.4
  [ce78b400] SimpleUnPack v1.1.0
  [66db9d55] SnoopPrecompile v1.0.3
  [a2af1166] SortingAlgorithms v1.2.0
  [47a9eef4] SparseDiffTools v2.12.0
  [e56a9233] Sparspak v0.3.9
  [276daf66] SpecialFunctions v2.3.1
  [aedffcd0] Static v0.8.8
  [0d7ed370] StaticArrayInterface v1.4.1
⌃ [90137ffa] StaticArrays v1.6.3
  [1e83bf80] StaticArraysCore v1.4.2
  [82ae8749] StatsAPI v1.7.0
  [2913bbd2] StatsBase v0.34.2
  [4c63d2b9] StatsFuns v1.3.0
  [7792a7ef] StrideArraysCore v0.5.1
  [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.18.0
  [a7c27f48] Unityper v0.1.5
  [3d5dd08c] VectorizationBase v0.21.64
  [19fa3120] VertexSafeGraphs v0.2.0
  [efce3f68] WoodburyMatrices v0.5.5
  [700de1a5] ZygoteRules v0.2.4
  [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

• 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

When trying to debug the issue I observed method _upwind_difference(...){Number} is going to
a different branch than _upwind_difference(...){AbstractVector}.
This was because D.offside was 1 in the Number-method, but 0 in the Vector method.

[xtalax]

I have seen this one recently, I will look closer soon

[xtalax]

I see the problem, you need to supply bcs for any of your observables that contain derivatives else they are not well specified. This needs a better error though.

[bgctw]

Thanks for looking into it so quickly.

Why do I need to specify additional boundary conditions (bcs) when specifying the same grid by a vector compared to specifying the same grid by a grid spacing? I thought the problem description did not change semantically.
For me its not quite intuitive to specify bcs for observables, such as a local sink decomposition flux, that are not related to transport.

Would it be sensible to assume zero-fluxes by default? Could this be implemented as a standard in the package that is only changed if bcs are explicitly specified for some observables?

For the delta-grid control case, this seems to be already the case. Then, is the advection-observable in the example for the first layer actually taking into account the specified bc for Y(t,0) or does it assume a default no-flux bc? Although the overall solution is correct, the user should be made aware that observables related to transport are not fully corresponding to the solution at the boundary.

[bgctw]

Thinking a bit more about it: its not a zero-flux condition that is applied by default with the delta-grid with the upwind scheme, but the flux at the border is the same as between the adjacent two cells/layers. Hence, the default condition corresponds rather to a zero-gradient in the flux, i.e. a zero second derivative of the amount.

[bgctw]

Based on @xtalax advice, I constructed a modified working vector-grid example with specifying addtional bcs. It does not throw the Boundary error. It skips most of the observed variables.

The MWE works only, if I specify bcs for the non-transport related dec_Y. The actual formulation of the boundary condition should not matter for the solution of the problem, but only for inspecting the observed variables at the borders.

When inspecting the simplified corresponding ODESystem (using symbolic_discretize):

(dec_Y(t))[15] ~ (Y(t))[15]*k_Y   # expected
 (dec_Y(t))[1] ~ 0.0002(12500.000000000002(dec_Y(t))[2] - 10000.000000000004(dec_Y(t))[3] + 2500.0000000000014(dec_Y(t))[4])
 (dec_Y(t))[16] ~ 0.00019999999999999963(2499.9999999999936(dec_Y(t))[13] - 9999.999999999982(dec_Y(t))[14] + 12499.999999999996(dec_Y(t))[15])

I would expect an observed variable: (dec_Y(t))[16] ~ (Y(t))[16]*k_Y

The example demonstrates that equations of observed variables unrelated to transport are currently strange at the borders - should be treated differently, and the user should not be required to specify boundary conditions for them.
Further, a default bc corresponding to the dz-number based discretization for the transport-related observed variables would be really helpful.

discretization_grid = MOLFiniteDifference([z => z_grid], t;
    advection_scheme = UpwindScheme(), approx_order = 2)
∂_zz = Differential(z)^2
eqs2 = [
    ∂_t(Y(t, z)) ~ i_Yz(t, z, i_Y) - dec_Y(t, z) + adv_Yo(t, z), 
    dec_Y(t, z) ~ k_Y * Y(t, z),          # observable of decomposition 
    adv_Yo(t, z) ~ -ω * ∂_z(Y(t, z)),  # observable advective flux of Y
]
bcs2 = [
    Y(0, z) ~ 1/z_m * Y0, 
    ∂_z(Y(t, 0)) ~ (Y(t, 0) - fagr(t,i_Y_agr, i_Y_agr_pulse)/ω) / dzs[1], 
    ∂_zz(Y(t, z_m)) ~ 0.0, 
    ∂_zz(dec_Y(t, 0)) ~ 0.0, 
    ∂_zz(dec_Y(t, z_m)) ~ 0.0, 
    #∂_zz(adv_Yo(t, 0)) ~ 0.0, 
    ∂_z(adv_Yo(t, 0)) ~ (fagr(t,i_Y_agr, i_Y_agr_pulse) - ω * Y(t, 0))/dzs[1], 
    ∂_zz(adv_Yo(t, z_m)) ~ 0.0, 
    ]
state_vars2 = [Y(t, z), dec_Y(t, z), adv_Yo(t,z)]
@named pdesys2 = PDESystem(eqs2, bcs2, domains, [t, z], state_vars2, params)
# Convert the PDE problem into an ODE problem
prob2 = discretize(pdesys2, discretization_grid)