My apologies if this is something really simple. I have looked at the docs and whatever examples I could find, but I am beginning to confuse the trees for the forest. Any directions would be greatly appreciated.
I am trying to learn ModelingToolkit by way of building a simple, isothermal axial dispersion reactor model with one reaction and three species. I started off by trying a dynamic model, but it looks like the PDE discretization tools are still WIP, so I changed this to a steady-state model.
using ModelingToolkit
#=
Differential molar balance for i
∂Cᵢ/∂t = 0 = Dₐ*∂²Cᵢ/∂z² - ∂(uCᵢ)/∂z + ∑(reaction j)νᵢⱼrᵢ
Peₐ = ul/Dₐ (l = characteristic length)
Dankwerts boundary conditions
z = 0 uₛCᵢ(z = 0) - Dₐ*dCₜₒₜ/dz(z=0) = uₛCᵢ(feed)
z = L dCᵢ/dz(z=L) = 0
Cₜₒₜ = P/RT
Cᵢ = nᵢ/nₜₒₜ*Cₜₒₜ
Reaction:
A + 2B -> C
r = k*pA*pC / (1 + Ka*pA + Kb*pB)
rA = -r
rB = -2r
rC = +r
=#
const N = 3 # Number of species
const ν = [-1.0, -2.0, 1.0]
const R₀ = 8.314
const Cfeed = [1.0, 1.0, 0.0]
const L = 2.0
@parameters T Dₐ k Ka Kb
@variables z u(z) C[1:N](z) p[1:N](z) r[1:N](z)
Dz = Differential(z)
Dzz = Differential(z)^2
# 1D ADPFR with Dankwerts boundary conditions
# Cᵢ = Pᵢ/R₀T
presconc = [C[i] ~ p[i]/R₀/T for i in 1:N]
# A + 2B -> C
# r = k*pA*pC / (1 + Ka*pA + Kb*pB)
reactions = [r[i] ~ ν[i] * (k * p[1] * p[2] / (1 + Ka*p[1] + Kb*p[2])) for i in 1:N]
# ∂Cᵢ/∂t = 0 = Dₐ*∂²Cᵢ/∂z² - ∂(uCᵢ)/∂z + ∑(reaction j)νᵢⱼrᵢ
axdisp = [0 ~ Dₐ*Dzz(C[i]) - Dz(u*C[i]) + r[i] for i in 1:N]
eqs = vcat(presconc, reactions, axdisp)
sys = ODESystem(eqs)
Each line executes as expected, up to the last, which returns an error:
ERROR: LoadError: MethodError: isequal(::Missing, ::Sym{Real, Base.ImmutableDict{DataType, Any}}) is ambiguous. Candidates:
isequal(x, ::SymbolicUtils.Symbolic) in SymbolicUtils at C:\Users\Braam\.julia\packages\SymbolicUtils\4n4Gz\src\types.jl:123
isequal(::Missing, ::Any) in Base at missing.jl:81
Possible fix, define
isequal(::Missing, ::SymbolicUtils.Symbolic)
Stacktrace:
[1] collect_var!(states::OrderedCollections.OrderedSet{Any}, parameters::OrderedCollections.OrderedSet{Any}, var::Missing, iv::Sym{Real, Base.ImmutableDict{DataType, Any}})
@ ModelingToolkit C:\Users\Braam\.julia\packages\ModelingToolkit\3q7jk\src\utils.jl:303
[2] collect_vars!(states::OrderedCollections.OrderedSet{Any}, parameters::OrderedCollections.OrderedSet{Any}, expr::SymbolicUtils.Add{Real, Int64, Dict{Any, Number}, Nothing}, iv::Sym{Real, Base.ImmutableDict{DataType, Any}})
@ ModelingToolkit C:\Users\Braam\.julia\packages\ModelingToolkit\3q7jk\src\utils.jl:296
[3] ODESystem(eqs::Vector{Equation}, iv::Nothing; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ ModelingToolkit C:\Users\Braam\.julia\packages\ModelingToolkit\3q7jk\src\systems\diffeqs\odesystem.jl:172
[4] ODESystem(eqs::Vector{Equation}, iv::Nothing) (repeats 2 times)
@ ModelingToolkit C:\Users\Braam\.julia\packages\ModelingToolkit\3q7jk\src\systems\diffeqs\odesystem.jl:151
[5] top-level scope
@ d:\JuliaCode\ADPFR\ADPFD.jl:49
in expression starting at d:\JuliaCode\ADPFR\ADPFD.jl:49
I am sure that I am doing something wrong, but after staring at this for too long, I can no longer see even obvious problems. If anyone could point me in a direction to resolve this, I would be extremely appreciative.