Getting an error when running discretize on a PDESystem with MOLFiniteDifference

[killah-t-cell]

Hi everyone,

I am getting an ArgumentError: reducing over an empty collection is not allowed when running discretize on a PDESystem with MOLFiniteDifference. Though I spent some time combing through the code before asking, I am not sure what I am doing wrong at this point.

Here is the code for reference. I am a beginner, so there might be some silly mistake. But the model works when I discretize it with PhysicsInformedNN

@parameters x, y, z, t
@variables A1(..), A2(..), A3(..), V(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
Dzz = Differential(z)^2
Dtt = Differential(t)^2
Dx = Differential(x)
Dy = Differential(y)
Dz = Differential(z)
Dt = Differential(t)

# Physics
μ_0 = 1.25663706212e-6
ε_0 = 8.8541878128e-12

# Helpers
divergenceA = Dx(A1(x, y, z, t)) + Dy(A2(x, y, z, t)) + Dz(A3(x, y, z, t))
laplacianV  = Dxx(V(x, y, z, t)) + Dyy(V(x, y, z, t)) + Dzz(V(x, y, z, t))
laplacianA1 = Dxx(A1(x, y, z, t)) + Dyy(A1(x, y, z, t)) + Dzz(A1(x, y, z, t))
laplacianA2 = Dxx(A2(x, y, z, t)) + Dyy(A2(x, y, z, t)) + Dzz(A2(x, y, z, t))
laplacianA3 = Dxx(A3(x, y, z, t)) + Dyy(A3(x, y, z, t)) + Dzz(A3(x, y, z, t))

S = [ρ/ε_0, μ_0*jx, μ_0*jy, μ_0*jz]
f = [V(x, y, z, t), A1(x, y, z, t), A2(x, y, z, t), A3(x, y, z, t)]
laplacians = [laplacianV, laplacianA1, laplacianA2, laplacianA3]

# Problem setup
eqs = Dtt.(f) .~ (S + laplacians)/(μ_0*ε_0)
# Boundary Conditions
bcs = [Dt(V(x, y, z, t)) ~ divergenceA / (μ_0 * ε_0), # set this as a condition
       A1(0,x, y, z) ~ sin(4*pi*z),
       A2(0,x, y, z) ~ cos(pi*x),
       A3(0,x, y, z) ~ sin(pi*x),
       V(0, x, y, z) ~ 1.0,
       A1(t,0, y, z) ~ exp(-t),
       A2(t,0, y, z) ~ exp(-t),
       A3(t,0, y, z) ~ exp(-t),
       V(t, 0, y, z) ~ exp(-t),
       A1(t,1.0, y, z) ~ exp(-t) * cos(1),
       A2(t,1.0, y, z) ~ exp(-t) * cos(1),
       A3(t,1.0, y, z) ~ exp(-t) * cos(1),
       V(t, 1.0, y, z) ~ exp(-t) * cos(1),
       A1(t,x, 0, z) ~ exp(-t),
       A2(t,x, 0, z) ~ exp(-t),
       A3(t,x, 0, z) ~ exp(-t),
       V(t, x, 0, z) ~ exp(-t),
       A1(t,x, 1.0, z) ~ exp(-t) * cos(1),
       A2(t,x, 1.0, z) ~ exp(-t) * cos(1),
       A3(t,x, 1.0, z) ~ exp(-t) * cos(1),
       V(t, x, 1.0, z) ~ exp(-t) * cos(1),
       A1(t,x, y, 0) ~ exp(-t),
       A2(t,x, y, 0) ~ exp(-t),
       A3(t,x, y, 0) ~ exp(-t),
       V(t, x, y, 0) ~ exp(-t),
       A1(t,x, y, 1.0) ~ exp(-t) * cos(1),
       A2(t,x, y, 1.0) ~ exp(-t) * cos(1),
       A3(t,x, y, 1.0) ~ exp(-t) * cos(1),
       V(t, x, y, 1.0) ~ exp(-t) * cos(1)]
# Space and time domains
domains = [t ∈ Interval(0.0,1.0),
           x ∈ Interval(0.0,1.0),
           y ∈ Interval(0.0,1.0),
           z ∈ Interval(0.0,1.0)]
pde_system = PDESystem(eqs,bcs,domains,[t, x, y, z],[V(x, y, z, t), A1(x, y, z, t), A2(x, y, z, t), A3(x, y, z, t)],[ρ=>0.5, jx=>2, jy=>1.2, jz=>4.8])
# Method of lines discretization
dx, dy, dz = 0.1, 0.1, 0.1
order = 2
discretization = MOLFiniteDifference([x=>dx, y=>dy, z=>dz],t)
prob = discretize(pde_system,discretization)
# Solve ODE problem
sol = solve(prob,Tsit5(),saveat=0.2)

[bgroenks96]

Possibly related to #415.

But you should make sure your operation Dtt.(f) .~ (S + laplacians)/(μ_0*ε_0) is properly broadcasted; i.e. it should be (I think) Dtt.(f) .~ (S .+ laplacians)./(μ_0*ε_0)

[killah-t-cell]

@bgroenks96 thanks! Unfortunately, I ran it again with the proper broadcasting and still ran into the same issue

[bgroenks96]

I'm not sure what the issues is... but I know I have seen the exact same error before when trying to set up a PDE problem.

[valentinsulzer]

Haven't tried this locally yet but the order of your independent variables (currently A(x,y,z,t)) might need to be the same as the order in the BCs (currently A(t,x,y,z)). We need way better error messages but better to wait for the function to be more stable

[killah-t-cell]

@tinosulzer nice catch! That in fact changed the error message.

Now I get

ERROR: LoadError: MethodError: no method matching ~(::Matrix{Num}, ::Float64)
Closest candidates are:
 ~(::Complex, ::Number) at /Users/gabrielbirnbaum/.julia/packages/Symbolics/sITWZ/src/equations.jl:71
 ~(::SymbolicUtils.Symbolic, ::Any) at /Users/gabrielbirnbaum/.julia/packages/Symbolics/sITWZ/src/equations.jl:67
 ~(::Num, ::Number) at /Users/gabrielbirnbaum/.julia/packages/Symbolics/sITWZ/src/equations.jl:64
 ...

Which I guess is promising because I can try to debug again. Will update here if I manage to fix it.

Re: better error messages – is there anything I can do to help? I'm still a bit of a Julia noob but can gladly try to participate. An error/warning like "you inverted your independent variables, are you sure about that?" would be super helpful.

[valentinsulzer]

Another thing is I don't think you can specify 2nd order (in time) systems yet, so you might need to explicitly rewrite it as a first-order system for now, but this is definitely functionality we should add

[killah-t-cell]

@tinosulzer Oh! Thanks for heads up. Beginner question: what's the best way to rewrite the system as a first order system?

[valentinsulzer]

Dtt(u) ~ f(u) --> Dt(u) ~ v, Dt(v) ~ f(u) (and make v a new variable with same dependent variables as u)

[killah-t-cell]

Quick update: I just rewrote it into a first-order system and tried to debug for a bit. Sadly, no change. It seems like the compiler throws an error before it reads eqs

For convenience here is an update version of the code.

using OrdinaryDiffEq, ModelingToolkit, DiffEqOperators, DomainSets

@parameters t, x, y, z 
@variables A1(..), A2(..), A3(..), V(..), w1(..), w2(..), w3(..), w4(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
Dzz = Differential(z)^2
Dtt = Differential(t)^2
Dx = Differential(x)
Dy = Differential(y)
Dz = Differential(z)
Dt = Differential(t)

# Physics
μ_0 = 1.25663706212e-6 # N A⁻²
ε_0 = 8.8541878128e-12 # F ms⁻¹

# Helpers
divergenceA = Dx(A1(t, x, y, z)) + Dy(A2(t, x, y, z)) + Dz(A3(t, x, y, z))
laplacianV  = Dxx(V(t, x, y, z)) + Dyy(V(t, x, y, z)) + Dzz(V(t, x, y, z))
laplacianA1 = Dxx(A1(t, x, y, z)) + Dyy(A1(t, x, y, z)) + Dzz(A1(t, x, y, z))
laplacianA2 = Dxx(A2(t, x, y, z)) + Dyy(A2(t, x, y, z)) + Dzz(A2(t, x, y, z))
laplacianA3 = Dxx(A3(t, x, y, z)) + Dyy(A3(t, x, y, z)) + Dzz(A3(t, x, y, z))

jx, jy, jz = 8.32, 3.1, 4.0 
ρ = 3.21

S = [ρ/ε_0, μ_0*jx, μ_0*jy, μ_0*jz]
f = [V(t, x, y, z), A1(t, x, y, z), A2(t, x, y, z), A3(t, x, y, z)]
w = [w1(t, x, y, z), w2(t, x, y, z), w3(t, x, y, z), w4(t, x, y, z)]
laplacians = [laplacianV, laplacianA1, laplacianA2, laplacianA3]

# Problem setup – Lorenz gauge electrodynamic potential
# eqs = Dtt.(f) .~ (S .+ laplacians)./(μ_0*ε_0)
first_order = Dt.(f) .~ w
eq = Dt.(w) .~ (S .+ laplacians)./(μ_0*ε_0)
eqs = [first_order; eq]

# Boundary Conditions
bcs = [Dt(V(t, x, y, z)) ~ divergenceA / (μ_0 * ε_0), # set this as a condition
       A1(0,x, y, z) ~ 0,
       A2(0,x, y, z) ~ 0,
       A3(0,x, y, z) ~ 0,
       V(0, x, y, z) ~ 1,
       A1(t,0, y, z) ~ 0,
       A2(t,0, y, z) ~ 0,
       A3(t,0, y, z) ~ 0,
       V(t, 0, y, z) ~ 0,
       A1(t,1.0, y, z) ~ 1,
       A2(t,1.0, y, z) ~ 1,
       A3(t,1.0, y, z) ~ 1,
       V(t, 1.0, y, z) ~ 1,
       A1(t,x, 0, z) ~ 0,
       A2(t,x, 0, z) ~ 0,
       A3(t,x, 0, z) ~ 0,
       V(t, x, 0, z) ~ 0,
       A1(t,x, 1.0, z) ~ 1,
       A2(t,x, 1.0, z) ~ 1,
       A3(t,x, 1.0, z) ~ 1,
       V(t, x, 1.0, z) ~ 1,
       A1(t,x, y, 0) ~ 0,
       A2(t,x, y, 0) ~ 0,
       A3(t,x, y, 0) ~ 0,
       V(t, x, y, 0) ~ 0,
       A1(t,x, y, 1.0) ~ 0,
       A2(t,x, y, 1.0) ~ 0,
       A3(t,x, y, 1.0) ~ 0,
       V(t, x, y, 1.0) ~ 0]

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

pde_system = PDESystem(eqs,bcs,domains,[t, x, y, z],[V(t, x, y, z), A1(t, x, y, z), A2(t, x, y, z), A3(t, x, y, z), w1(t, x, y, z), w2(t, x, y, z), w3(t, x, y, z), w4(t, x, y, z)])

# Method of lines discretization
dx, dy, dz = 0.1, 0.1, 0.1
order = 2
discretization = MOLFiniteDifference([x=>dx, y=>dy, z=>dz],t)

prob = discretize(pde_system,discretization) 

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

The error is the same ERROR: LoadError: MethodError: no method matching ~(::Matrix{Num}, ::Float64) as before. I will try to define this numerically with DiffEqOperators for now.

[valentinsulzer]

Looks like the problem could be broadcasting the ~ when defining your eqns. I don't know whether this is supported. Could you define your equations one by one in a for loop for now? This won't come at any performance cost since discretize rewrites the equations optimally anyway

[killah-t-cell]

That's a good guess but unfortunately it didn't do anything. the compiler threw an error instantly (which I took it as a cache-hint that the problem is not with eqs). The problem seems to be in the boundary conditions.

I rewrote the boundary into Dirichlet boundary conditions (just like in this example https://diffeqoperators.sciml.ai/stable/symbolic_tutorials/mol_heat/#Dirichlet-boundary-conditions). The compiler had to evaluate for a bit, and then I got the most interesting of errors:

LoadError: MethodError: no method matching -(::Term{Real, Nothing}, ::Matrix{Num})
Closest candidates are:
  -(::VectorizationBase.CartesianVIndex, ::Any) at /Users/gabrielbirnbaum/.julia/packages/VectorizationBase/EjKep/src/cartesianvindex.jl:58
  -(::StaticArrays.StaticArray, ::AbstractArray) at /Users/gabrielbirnbaum/.julia/packages/StaticArrays/aMe0r/src/linalg.jl:18
  -(::ChainRulesCore.AbstractThunk, ::Any) at /Users/gabrielbirnbaum/.julia/packages/ChainRulesCore/e5hAX/src/differentials/thunks.jl:30
  ...
Stacktrace:
 [1] initialize_system_structure(sys::ODESystem)
   @ ModelingToolkit.SystemStructures ~/.julia/packages/ModelingToolkit/MuQfD/src/systems/systemstructure.jl:112
 [2] alias_elimination(sys::ODESystem)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/MuQfD/src/systems/alias_elimination.jl:6
 [3] structural_simplify(sys::ODESystem)
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/MuQfD/src/systems/abstractsystem.jl:642
 [4] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{Vector{Pair{Num, Float64}}, Num})
   @ DiffEqOperators ~/.julia/packages/DiffEqOperators/FnmOt/src/MOLFiniteDifference/MOL_discretization.jl:373
 [5] top-level scope
   @ ~/.julia/dev/MaxwellEqSolver/src/numerical.jl:96

See the boundary array below for reference. I think discretize is somehow converting the @parameter t or each equation in the array into a matrix during evaluation? May be something that it needs to do for a system of equations or higher dimensions?

bcs = [A1(0,x, y, z) ~ cos(x),
       A2(0,x, y, z) ~ cos(x),
       A3(0,x, y, z) ~ sin(z),
       V(0, x, y, z) ~ 2,
       A1(t,0, y, z) ~ exp(-t),
       A2(t,0, y, z) ~ exp(-t),
       A3(t,0, y, z) ~ exp(-t),
       V(t, 0, y, z) ~ exp(-t),
       A1(t,1.0, y, z) ~ exp(-t) * cos(1),
       A2(t,1.0, y, z) ~ exp(-t) * cos(1),
       A3(t,1.0, y, z) ~ exp(-t) * cos(1),
       V(t, 1.0, y, z) ~ exp(-t) * cos(1),
       A1(t,x, 0, z) ~ exp(-t),
       A2(t,x, 0, z) ~ exp(-t),
       A3(t,x, 0, z) ~ exp(-t),
       V(t, x, 0, z) ~ exp(-t),
       A1(t,x, 1.0, z) ~ exp(-t) * cos(1),
       A2(t,x, 1.0, z) ~ exp(-t) * cos(1),
       A3(t,x, 1.0, z) ~ exp(-t) * cos(1),
       V(t, x, 1.0, z) ~ exp(-t) * cos(1),
       A1(t,x, y, 0) ~ exp(-t),
       A2(t,x, y, 0) ~ exp(-t),
       A3(t,x, y, 0) ~ exp(-t),
       V(t, x, y, 0) ~ exp(-t),
       A1(t,x, y, 1.0) ~ exp(-t) * cos(1),
       A2(t,x, y, 1.0) ~ exp(-t) * cos(1),
       A3(t,x, y, 1.0) ~ exp(-t) * cos(1),
       V(t, x, y, 1.0) ~ exp(-t) * cos(1)]

[valentinsulzer]

Have identified the problem but don't fully understand it or know how to fix it. It's a bug in

DiffEqOperators.jl/src/MOLFiniteDifference/MOL_discretization.jl

Lines 208 to 211 in 40aa4b9

	lhs = substitute(lhs,depvarbcmaps[i])
	rhs = substitute.((bc.rhs,),edgemaps[bcargs])
	lhs = lhs isa Vector ? lhs : [lhs] # handle 1D
	push!(bceqs,lhs .~ rhs)

with the boundary conditions in 3D.

2D example

> depvarbcmaps[i]
u(t, 0.0, y) => Num[u₁ˏ₁(t), u₁ˏ₂(t), u₁ˏ₃(t), u₁ˏ₄(t), u₁ˏ₅(t), u₁ˏ₆(t), u₁ˏ₇(t), u₁ˏ₈(t), u₁ˏ₉(t), u₁ˏ₁₀(t), u₁ˏ₁₁(t)]
> lhs
Num[u₁ˏ₁(t), u₁ˏ₂(t), u₁ˏ₃(t), u₁ˏ₄(t), u₁ˏ₅(t), u₁ˏ₆(t), u₁ˏ₇(t), u₁ˏ₈(t), u₁ˏ₉(t), u₁ˏ₁₀(t), u₁ˏ₁₁(t)]
> rhs
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

3D example

> depvarbcmaps[i]
V(t, 0.0, y, z) => Num[V₁ˏ₁ˏ₁(t) V₁ˏ₁ˏ₂(t) V₁ˏ₁ˏ₃(t) V₁ˏ₁ˏ₄(t) V₁ˏ₁ˏ₅(t) V₁ˏ₁ˏ₆(t) V₁ˏ₁ˏ₇(t) V₁ˏ₁ˏ₈(t) V₁ˏ₁ˏ₉(t) V₁ˏ₁ˏ₁₀(t) V₁ˏ₁ˏ₁₁(t); V₁ˏ₂ˏ₁(t) V₁ˏ₂ˏ₂(t) V₁ˏ₂ˏ₃(t) V₁ˏ₂ˏ₄(t) V₁ˏ₂ˏ₅(t) V₁ˏ₂ˏ₆(t) V₁ˏ₂ˏ₇(t) V₁ˏ₂ˏ₈(t) V₁ˏ₂ˏ₉(t) V₁ˏ₂ˏ₁₀(t) V₁ˏ₂ˏ₁₁(t); V₁ˏ₃ˏ₁(t) V₁ˏ₃ˏ₂(t) V₁ˏ₃ˏ₃(t) V₁ˏ₃ˏ₄(t) V₁ˏ₃ˏ₅(t) V₁ˏ₃ˏ₆(t) V₁ˏ₃ˏ₇(t) V₁ˏ₃ˏ₈(t) V₁ˏ₃ˏ₉(t) V₁ˏ₃ˏ₁₀(t) V₁ˏ₃ˏ₁₁(t); V₁ˏ₄ˏ₁(t) V₁ˏ₄ˏ₂(t) V₁ˏ₄ˏ₃(t) V₁ˏ₄ˏ₄(t) V₁ˏ₄ˏ₅(t) V₁ˏ₄ˏ₆(t) V₁ˏ₄ˏ₇(t) V₁ˏ₄ˏ₈(t) V₁ˏ₄ˏ₉(t) V₁ˏ₄ˏ₁₀(t) V₁ˏ₄ˏ₁₁(t); V₁ˏ₅ˏ₁(t) V₁ˏ₅ˏ₂(t) V₁ˏ₅ˏ₃(t) V₁ˏ₅ˏ₄(t) V₁ˏ₅ˏ₅(t) V₁ˏ₅ˏ₆(t) V₁ˏ₅ˏ₇(t) V₁ˏ₅ˏ₈(t) V₁ˏ₅ˏ₉(t) V₁ˏ₅ˏ₁₀(t) V₁ˏ₅ˏ₁₁(t); V₁ˏ₆ˏ₁(t) V₁ˏ₆ˏ₂(t) V₁ˏ₆ˏ₃(t) V₁ˏ₆ˏ₄(t) V₁ˏ₆ˏ₅(t) V₁ˏ₆ˏ₆(t) V₁ˏ₆ˏ₇(t) V₁ˏ₆ˏ₈(t) V₁ˏ₆ˏ₉(t) V₁ˏ₆ˏ₁₀(t) V₁ˏ₆ˏ₁₁(t); V₁ˏ₇ˏ₁(t) V₁ˏ₇ˏ₂(t) V₁ˏ₇ˏ₃(t) V₁ˏ₇ˏ₄(t) V₁ˏ₇ˏ₅(t) V₁ˏ₇ˏ₆(t) V₁ˏ₇ˏ₇(t) V₁ˏ₇ˏ₈(t) V₁ˏ₇ˏ₉(t) V₁ˏ₇ˏ₁₀(t) V₁ˏ₇ˏ₁₁(t); V₁ˏ₈ˏ₁(t) V₁ˏ₈ˏ₂(t) V₁ˏ₈ˏ₃(t) V₁ˏ₈ˏ₄(t) V₁ˏ₈ˏ₅(t) V₁ˏ₈ˏ₆(t) V₁ˏ₈ˏ₇(t) V₁ˏ₈ˏ₈(t) V₁ˏ₈ˏ₉(t) V₁ˏ₈ˏ₁₀(t) V₁ˏ₈ˏ₁₁(t); V₁ˏ₉ˏ₁(t) V₁ˏ₉ˏ₂(t) V₁ˏ₉ˏ₃(t) V₁ˏ₉ˏ₄(t) V₁ˏ₉ˏ₅(t) V₁ˏ₉ˏ₆(t) V₁ˏ₉ˏ₇(t) V₁ˏ₉ˏ₈(t) V₁ˏ₉ˏ₉(t) V₁ˏ₉ˏ₁₀(t) V₁ˏ₉ˏ₁₁(t); V₁ˏ₁₀ˏ₁(t) V₁ˏ₁₀ˏ₂(t) V₁ˏ₁₀ˏ₃(t) V₁ˏ₁₀ˏ₄(t) V₁ˏ₁₀ˏ₅(t) V₁ˏ₁₀ˏ₆(t) V₁ˏ₁₀ˏ₇(t) V₁ˏ₁₀ˏ₈(t) V₁ˏ₁₀ˏ₉(t) V₁ˏ₁₀ˏ₁₀(t) V₁ˏ₁₀ˏ₁₁(t); V₁ˏ₁₁ˏ₁(t) V₁ˏ₁₁ˏ₂(t) V₁ˏ₁₁ˏ₃(t) V₁ˏ₁₁ˏ₄(t) V₁ˏ₁₁ˏ₅(t) V₁ˏ₁₁ˏ₆(t) V₁ˏ₁₁ˏ₇(t) V₁ˏ₁₁ˏ₈(t) V₁ˏ₁₁ˏ₉(t) V₁ˏ₁₁ˏ₁₀(t) V₁ˏ₁₁ˏ₁₁(t)]
> lhs
Matrix{Num}[[V₁ˏ₁ˏ₁(t) V₁ˏ₁ˏ₂(t) V₁ˏ₁ˏ₃(t) V₁ˏ₁ˏ₄(t) V₁ˏ₁ˏ₅(t) V₁ˏ₁ˏ₆(t) V₁ˏ₁ˏ₇(t) V₁ˏ₁ˏ₈(t) V₁ˏ₁ˏ₉(t) V₁ˏ₁ˏ₁₀(t) V₁ˏ₁ˏ₁₁(t); V₁ˏ₂ˏ₁(t) V₁ˏ₂ˏ₂(t) V₁ˏ₂ˏ₃(t) V₁ˏ₂ˏ₄(t) V₁ˏ₂ˏ₅(t) V₁ˏ₂ˏ₆(t) V₁ˏ₂ˏ₇(t) V₁ˏ₂ˏ₈(t) V₁ˏ₂ˏ₉(t) V₁ˏ₂ˏ₁₀(t) V₁ˏ₂ˏ₁₁(t); V₁ˏ₃ˏ₁(t) V₁ˏ₃ˏ₂(t) V₁ˏ₃ˏ₃(t) V₁ˏ₃ˏ₄(t) V₁ˏ₃ˏ₅(t) V₁ˏ₃ˏ₆(t) V₁ˏ₃ˏ₇(t) V₁ˏ₃ˏ₈(t) V₁ˏ₃ˏ₉(t) V₁ˏ₃ˏ₁₀(t) V₁ˏ₃ˏ₁₁(t); V₁ˏ₄ˏ₁(t) V₁ˏ₄ˏ₂(t) V₁ˏ₄ˏ₃(t) V₁ˏ₄ˏ₄(t) V₁ˏ₄ˏ₅(t) V₁ˏ₄ˏ₆(t) V₁ˏ₄ˏ₇(t) V₁ˏ₄ˏ₈(t) V₁ˏ₄ˏ₉(t) V₁ˏ₄ˏ₁₀(t) V₁ˏ₄ˏ₁₁(t); V₁ˏ₅ˏ₁(t) V₁ˏ₅ˏ₂(t) V₁ˏ₅ˏ₃(t) V₁ˏ₅ˏ₄(t) V₁ˏ₅ˏ₅(t) V₁ˏ₅ˏ₆(t) V₁ˏ₅ˏ₇(t) V₁ˏ₅ˏ₈(t) V₁ˏ₅ˏ₉(t) V₁ˏ₅ˏ₁₀(t) V₁ˏ₅ˏ₁₁(t); V₁ˏ₆ˏ₁(t) V₁ˏ₆ˏ₂(t) V₁ˏ₆ˏ₃(t) V₁ˏ₆ˏ₄(t) V₁ˏ₆ˏ₅(t) V₁ˏ₆ˏ₆(t) V₁ˏ₆ˏ₇(t) V₁ˏ₆ˏ₈(t) V₁ˏ₆ˏ₉(t) V₁ˏ₆ˏ₁₀(t) V₁ˏ₆ˏ₁₁(t); V₁ˏ₇ˏ₁(t) V₁ˏ₇ˏ₂(t) V₁ˏ₇ˏ₃(t) V₁ˏ₇ˏ₄(t) V₁ˏ₇ˏ₅(t) V₁ˏ₇ˏ₆(t) V₁ˏ₇ˏ₇(t) V₁ˏ₇ˏ₈(t) V₁ˏ₇ˏ₉(t) V₁ˏ₇ˏ₁₀(t) V₁ˏ₇ˏ₁₁(t); V₁ˏ₈ˏ₁(t) V₁ˏ₈ˏ₂(t) V₁ˏ₈ˏ₃(t) V₁ˏ₈ˏ₄(t) V₁ˏ₈ˏ₅(t) V₁ˏ₈ˏ₆(t) V₁ˏ₈ˏ₇(t) V₁ˏ₈ˏ₈(t) V₁ˏ₈ˏ₉(t) V₁ˏ₈ˏ₁₀(t) V₁ˏ₈ˏ₁₁(t); V₁ˏ₉ˏ₁(t) V₁ˏ₉ˏ₂(t) V₁ˏ₉ˏ₃(t) V₁ˏ₉ˏ₄(t) V₁ˏ₉ˏ₅(t) V₁ˏ₉ˏ₆(t) V₁ˏ₉ˏ₇(t) V₁ˏ₉ˏ₈(t) V₁ˏ₉ˏ₉(t) V₁ˏ₉ˏ₁₀(t) V₁ˏ₉ˏ₁₁(t); V₁ˏ₁₀ˏ₁(t) V₁ˏ₁₀ˏ₂(t) V₁ˏ₁₀ˏ₃(t) V₁ˏ₁₀ˏ₄(t) V₁ˏ₁₀ˏ₅(t) V₁ˏ₁₀ˏ₆(t) V₁ˏ₁₀ˏ₇(t) V₁ˏ₁₀ˏ₈(t) V₁ˏ₁₀ˏ₉(t) V₁ˏ₁₀ˏ₁₀(t) V₁ˏ₁₀ˏ₁₁(t); V₁ˏ₁₁ˏ₁(t) V₁ˏ₁₁ˏ₂(t) V₁ˏ₁₁ˏ₃(t) V₁ˏ₁₁ˏ₄(t) V₁ˏ₁₁ˏ₅(t) V₁ˏ₁₁ˏ₆(t) V₁ˏ₁₁ˏ₇(t) V₁ˏ₁₁ˏ₈(t) V₁ˏ₁₁ˏ₉(t) V₁ˏ₁₁ˏ₁₀(t) V₁ˏ₁₁ˏ₁₁(t)]]
> rhs
[0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0]

So the problem is that in 3D lhs is a Matrix{Num} instead of a Num, and then substitute doesn't work. This is despite depvarbcmaps[i] still being a Num.

[ChrisRackauckas]

@shashi

[shashi]

Actually I'm not sure if I'm going to be more helpful than pointing out what you already know haha:

 lhs = substitute(lhs,depvarbcmaps[i]) 

Will only work if lhs is a single thing. Otherwise it treats whatever you give it as a literal and returns it without modification. This is by design.

You may want to check to see if lhs isa AbstractArray and broadcast.

Also:

 rhs = substitute.((bc.rhs,),edgemaps[bcargs]) 

You may get a better performance if you did: substitute(bc.rhs, merge(edgemaps[bcargs]...)).