Error trying to discretize a PDAE

[valentinsulzer]

When one of the equations is of the form

0 ~ Dxx(u(t,x))

MOLFiniteDifference does not work because it tries to do SymbolicUtils.operation(0).

Then, ODESystem is hardcoded - how should the resulting DAE system be created? Edit: maybe it's fine and the problem is structural_simplify is being too smart for me

[valentinsulzer]

Complete example:

using ModelingToolkit
using IfElse

# MTK model
# ('negative electrode', 'separator', 'positive electrode') -> x
@parameters t x
# 'Electrolyte concentration ' -> c_e
# 'Electrolyte potential' -> phi_e
@variables c_e(..) phi_e(..)
Dt = Differential(t)
Dx = Differential(x)

# 'Electrolyte concentration ' equation

function concatenation(n, s, p)
   # A concatenation in the electrolyte domain
  f= (x) -> IfElse.ifelse(
      x < 0.4444444444444445, n, IfElse.ifelse(
         x < 0.5555555555555556, s, p
      )
   )
end

cache_5101060308695467050(x) = concatenation(2.25, 0.0, -2.25)(x)
@register cache_5101060308695467050(x)
cache_8771569224475106856 = (Dx(Dx(c_e(t, x)))) + cache_5101060308695467050(x)

# 'Electrolyte potential' equation
cache_5101060308695467050(x) = concatenation(2.25, 0.0, -2.25)(x)
@register cache_5101060308695467050(x)
# cache_4483180157687090897 = (Dx(1 / c_e(t, x) * Dx(c_e(t, x))) - Dx(Dx(phi_e(t, x)))) - cache_5101060308695467050(x)
cache_4483180157687090897 = - Dx(Dx(phi_e(t, x))) - cache_5101060308695467050(x)


eqs = [
   Dt(c_e(t, x)) ~ cache_8771569224475106856,
   0 ~ cache_4483180157687090897,
]

ics_bcs = [
   # initial conditions
   c_e(0, x) ~ 1.0,
   phi_e(0, x) ~ 0.0,
   # boundary conditions
   Dx(c_e(t, 0.0)) ~ 0.0,
   Dx(c_e(t, 1.0)) ~ 0.0,
   phi_e(t, 0.0) ~ 0.0,
   Dx(phi_e(t, 1.0)) ~ 0.0,
]

t_domain = IntervalDomain(0.0, 3600.0)
x_domain = IntervalDomain(0.0, 1.0)

domains = [
   t in t_domain,
   x in x_domain,
]
ind_vars = [t, x]
dep_vars = [c_e(t,x), phi_e(t,x)]

reduced_c_phi_pde_system = PDESystem(eqs, ics_bcs, domains, ind_vars, dep_vars)


# Finite difference solution
using DiffEqOperators
using OrdinaryDiffEq

ln = 4/9
ls = 1/9
dx = vcat(range(0,ln,length=11)[1:end-1], range(ln,ln+ls,length=11)[1:end-1], range(ln+ls,1,length=11))
discretization = MOLFiniteDifference([x=>dx],t; grid_align=edge_align)
prob = discretize(reduced_c_phi_pde_system,discretization)

[valentinsulzer]

Looks like everything is working fine and it's just structural_simplify magic:
after discretization

julia> sys.eqs
64-element Vector{Equation}:
 Differential(t)(c_e₂(t)) ~ 2.25 + 506.25c_e₁(t) + 506.25c_e₃(t) - (1012.5c_e₂(t))
 Differential(t)(c_e₃(t)) ~ 2.25 + 506.2500000000001c_e₂(t) + 506.2500000000001c_e₄(t) - (1012.5000000000002c_e₃(t))
 Differential(t)(c_e₄(t)) ~ 2.25 + 506.25000000000006c_e₃(t) + 506.2499999999999c_e₅(t) - (1012.5c_e₄(t))
 Differential(t)(c_e₅(t)) ~ 2.25 + 506.24999999999983c_e₄(t) + 506.2499999999998c_e₆(t) - (1012.4999999999995c_e₅(t))
 Differential(t)(c_e₆(t)) ~ 2.25 + 506.25c_e₅(t) + 506.25000000000034c_e₇(t) - (1012.5000000000003c_e₆(t))
 Differential(t)(c_e₇(t)) ~ 2.25 + 506.25000000000045c_e₆(t) + 506.25000000000045c_e₈(t) - (1012.5000000000009c_e₇(t))
 ⋮
 0 ~ 2.25 + 1012.5000000000022phi_e₃₁(t) - (506.24999999999983phi_e₃₀(t)) - (506.25000000000233phi_e₃₂(t))
 22.5c_e₂(t) - (22.5c_e₁(t)) ~ 0.0
 22.50000000000008c_e₃₂(t) - (22.50000000000008c_e₃₁(t)) ~ 0.0
 0.5phi_e₁(t) + 0.5phi_e₂(t) ~ 0.0
 22.50000000000008phi_e₃₂(t) - (22.50000000000008phi_e₃₁(t)) ~ 0.0

after structural_simplify

julia> simpsys.eqs
32-element Vector{Equation}:
 Differential(t)(c_e₂(t)) ~ 2.25 + 506.25c_e₃(t) - (506.25c_e₂(t))
 Differential(t)(c_e₃(t)) ~ 2.25 + 506.2500000000001c_e₂(t) + 506.2500000000001c_e₄(t) - (1012.5000000000002c_e₃(t))
 Differential(t)(c_e₄(t)) ~ 2.25 + 506.25000000000006c_e₃(t) + 506.2499999999999c_e₅(t) - (1012.5c_e₄(t))
 Differential(t)(c_e₅(t)) ~ 2.25 + 506.24999999999983c_e₄(t) + 506.2499999999998c_e₆(t) - (1012.4999999999995c_e₅(t))
 Differential(t)(c_e₆(t)) ~ 2.25 + 506.25c_e₅(t) + 506.25000000000034c_e₇(t) - (1012.5000000000003c_e₆(t))
 Differential(t)(c_e₇(t)) ~ 2.25 + 506.25000000000045c_e₆(t) + 506.25000000000045c_e₈(t) - (1012.5000000000009c_e₇(t))
 ⋮
 Differential(t)(c_e₂₉(t)) ~ 506.25000000000045c_e₂₈(t) + 506.2499999999992c_e₃₀(t) - 2.25 - (1012.4999999999997c_e₂₉(t))
 Differential(t)(c_e₃₀(t)) ~ 506.2499999999985c_e₂₉(t) + 506.2499999999985c_e₃₁(t) - 2.25 - (1012.499999999997c_e₃₀(t))
 Differential(t)(c_e₃₁(t)) ~ 506.24999999999983c_e₃₀(t) - 2.25 - (506.2499999999999c_e₃₁(t))
 0 ~ 0.0022222222222222222 + 0.5phi_e₃(t) - (1.5phi_e₂(t))
 0 ~ 22.499881744384766phi_e₂(t) - 0.10000050067901611 - (22.499717712402344phi_e₃(t))

How would I recover sol[phi_e] in this case?

[ChrisRackauckas]

sol[phi_e[i]]

[valentinsulzer]

That gives me

julia> sol[phi_e[1]]
ERROR: Cannot convert getindex(phi_e(::Vararg{Any, N} where N)::Real, 1) to a symbol
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:33
  [2] tosymbol(t::Term{Real, Nothing}; states::Nothing, escape::Bool)
    @ Symbolics ~/.julia/packages/Symbolics/h8kPL/src/utils.jl:152
  [3] tosymbol(t::Term{Real, Nothing})
    @ Symbolics ~/.julia/packages/Symbolics/h8kPL/src/utils.jl:136
  [4] tosymbol(t::Num; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Symbolics ~/.julia/packages/Symbolics/h8kPL/src/utils.jl:76
  [5] tosymbol
    @ ~/.julia/packages/Symbolics/h8kPL/src/utils.jl:76 [inlined]
  [6] Symbol
    @ ~/.julia/packages/Symbolics/h8kPL/src/utils.jl:73 [inlined]
  [7] sym_to_index(sym::Num, syms::Vector{Symbol})
    @ SciMLBase ~/.julia/packages/SciMLBase/EFFG1/src/solutions/solution_interface.jl:293
  [8] sym_to_index
    @ ~/.julia/packages/SciMLBase/EFFG1/src/solutions/solution_interface.jl:292 [inlined]
  [9] getindex(A::ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Nothing, ODEFunction{true, ModelingToolkit.var"#f#151"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#540"), Symbol("##arg#541"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x11dc4677, 0xcfefcb0b, 0xf8d8e25e, 0xbcc9106f, 0xb11a8eaa)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#542"), Symbol("##arg#540"), Symbol("##arg#541"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6951d4c3, 0x17288a49, 0x7682651f, 0x663e1546, 0x0aabaf53)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, ModelingToolkit.var"#124#generated_observed#158"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, KenCarp47{0, true, DefaultLinSolve, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, DataType}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, ModelingToolkit.var"#f#151"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#540"), Symbol("##arg#541"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x11dc4677, 0xcfefcb0b, 0xf8d8e25e, 0xbcc9106f, 0xb11a8eaa)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#542"), Symbol("##arg#540"), Symbol("##arg#541"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6951d4c3, 0x17288a49, 0x7682651f, 0x663e1546, 0x0aabaf53)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, ModelingToolkit.var"#124#generated_observed#158"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.KenCarp47Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, Matrix{Float64}, Matrix{Float64}, SciMLBase.UJacobianWrapper{ODEFunction{true, ModelingToolkit.var"#f#151"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#540"), Symbol("##arg#541"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x11dc4677, 0xcfefcb0b, 0xf8d8e25e, 0xbcc9106f, 0xb11a8eaa)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#542"), Symbol("##arg#540"), Symbol("##arg#541"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6951d4c3, 0x17288a49, 0x7682651f, 0x663e1546, 0x0aabaf53)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, ModelingToolkit.var"#124#generated_observed#158"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Float64, Nothing}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, ModelingToolkit.var"#f#151"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#540"), Symbol("##arg#541"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x11dc4677, 0xcfefcb0b, 0xf8d8e25e, 0xbcc9106f, 0xb11a8eaa)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#542"), Symbol("##arg#540"), Symbol("##arg#541"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6951d4c3, 0x17288a49, 0x7682651f, 0x663e1546, 0x0aabaf53)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, ModelingToolkit.var"#124#generated_observed#158"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Float64, Nothing}, Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{SciMLBase.UJacobianWrapper{ODEFunction{true, ModelingToolkit.var"#f#151"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##arg#540"), Symbol("##arg#541"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x11dc4677, 0xcfefcb0b, 0xf8d8e25e, 0xbcc9106f, 0xb11a8eaa)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(Symbol("##out#542"), Symbol("##arg#540"), Symbol("##arg#541"), :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6951d4c3, 0x17288a49, 0x7682651f, 0x663e1546, 0x0aabaf53)}}, Matrix{Float64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, ModelingToolkit.var"#124#generated_observed#158"{Bool, ODESystem, Dict{Any, Any}}, Nothing}, Float64, Nothing}, Float64}, Float64, 12}}, Vector{Float64}, Vector{Vector{NTuple{12, Bool}}}, UnitRange{Int64}, Nothing}, DefaultLinSolve}}, OrdinaryDiffEq.KenCarp47Tableau{Float64, Float64}, Nothing}}, DiffEqBase.DEStats}, sym::Num)
    @ SciMLBase ~/.julia/packages/SciMLBase/EFFG1/src/solutions/solution_interface.jl:39
 [10] top-level scope

[ChrisRackauckas]

@variables phi_e[1:n]

[valentinsulzer]

@variables phi_e[1:n](..)
sol[phi_e[1](t)]

seems to be working

[ChrisRackauckas]

Yeah, we can setup sol[phi_e], but that's one of the symbolic array things.