Hello all! I’m new to Julia and am trying to use the Method of Lines and MTK packages to solve a system of PDEs. I am running into repeated method errors, and the stacktrace keeps pointing to a line from Macro Expanded code, and I’m not sure where to start with troubleshooting. Here’s the code and stack trace.
using ModelingToolkit
using DifferentialEquations
using MethodOfLines
using DomainSets
using Symbolics
#Parameters and Variables:
JJ = 3
N = 45
P = 20
@parameters r z t
@parameters D ε εp K kI u_avg qm Cf rp[1:JJ] psd[1:JJ] kf[1:JJ] zL
@variables c(..)[1:JJ] q(..)[1:JJ] Csl(..) Csh(..) #Pore Fluid Concentration, Sorption Concentration, Bulk Concentration of Slice, Bulk Concentration of Column
@variables c_surface(..)[1:JJ] #Helper Variable for cj(r=rpj, t)
@variables Q(..) #Total Pore Concentration
@variables cstar(..)[1:JJ] #Equilibrium Concentration
#Differential Operators:
Dt = Differential(t)
Dr = Differential(r)
Drr = Differential(r)^2
Dz = Differential(z)
#Equations
eqs = []
bcs = []
rp_vals = [50,62,75]
kf_vals = 1.54e-5./rp_vals
frac = [0.22, 0.54, 0.24]
psd_vals = kf_vals./rp_vals .*frac
push!(eqs, Dt(εp *c(r,z,t)[j]) + Dt(q(r,z,t)[j]) ~ D*(2/r*Dr(c(r,z,t)[j]) + Drr(c(r,z,t)[j])) for j in 1:JJ ) #Diffusion into a pore
push!(eqs, c_surface(z, t)[j] ~ c(rp[j],z, t)[j] for j in 1:JJ) #Surface Flux
push!(eqs, cstar(r,z,t)[j] ~ q(r,z,t)[j] / (K * (qm - q(r,z,t)[j])) for j in 1:JJ) #Langmuir Isotherm Equilibrium Relationship
push!(eqs, Dt(q(r,z,t)[j] ~ kI * (c_surface(z,t)[j] - cstar(r,z,t)[j])) for j in 1:JJ) #LDF Model for Equilibrium
push!(eqs, Q(z,t) ~ (εp * c(r,z,t)[j] + q(r,z,t)[j]) for j in 1:JJ) #Total Concentration entering Pore
push!(eqs, Dt(Q(z,t) * -(1/3) * (1-ε / ε)) ~ sum(psd[j] * (Csl(z,t) - c_surface(z,t)[j]) for j in 1:JJ) ) # Column Material Balance
push!(eqs, Dt(Csh(z,t)) ~ (-(1-ε)/ε) * Dt(Q(z,t)) - (u_avg / ε) * Dz(Csl(z,t))) #Column Balance
#Boundary Conditions
push!(bcs, Dr(c(0, z, t)[j]) ~ 0 for j in JJ) #Inside of Particle Neumann
push!(bcs, Dr(c(rp[j],z, t)[j]) ~ kf[j] * (Csl(z, t) - c_surface(z,t)[j]) for j in JJ ) #Surface of Particle Robin
push!(bcs, Dz(Csl(zL, t)) ~ 0) #Bottom of Column Neumann
push!(bcs, Csl(0, 0) ~ Cf) #Feed Concentration at top of column Dirichlet at both z = 0, t = 0
#Space and Time Domains
dmns = []
push!(dmns, r ∈ Interval(0, maximum(rp_vals))) #Interval for particle radii
push!(dmns, z ∈ Interval(0, zL)) #Interval for axial dimension from 0 to column length zL
push!(dmns, t ∈ Interval(0, Inf)) #Time Interval, based on spatial bounds
#Defining the PDE Systems
@named pdesys = PDESystem(eqs, bcs, dmns, [r,z,t], [c(r,z,t)[1:JJ], q(r,z,t)[1:JJ], Csl(z,t), Csh(z,t), Q(z,t), c_surface(z,t)[1:JJ], cstar(r,z,t)[1:JJ]],
[D => 5e-6, ε => 0.25, εp => 0.1, K => 367, kI => 4000, u_avg => 2.5, qm => 242, Cf => 3.0, rp[1:JJ] .=> rp_vals, psd[1:JJ] .=> psd_vals, kf[1:JJ] .=> kf_vals, zL => 2 ])
#Method of Lines Discretizations
order = 2
discretization = MOLFiniteDifference([r => N, z => P], t, approx_order = order)
#Convert the Problem from a PDE to an ODE
@time problem = discretize(pdesys, discretization)
#Solve the ODE
using DifferentialEquations
@time solution = solve(problem, FBDF(), saveat = 0.1)
#Plot Final Concentration vs Time
using Plots
plot(solution[Csh(z,t)], linewidth = 2, xaxis = "Time(t)", yaxis = "Concentration")
and here’s the stack trace
ERROR: MethodError: no method matching real(::Base.Generator{UnitRange{Int64}, var"#11#12"})
The function `real` exists, but no method is defined for this combination of argument types.
Closest candidates are:
real(::Type{Union{}}, Any...)
@ Base complex.jl:123
real(::Symbolics.SemiMonomial)
@ Symbolics C:\Users\U1087908\.julia\packages\Symbolics\lJiB2\src\semipoly.jl:123
real(::Missing)
@ Base missing.jl:101
...
Stacktrace:
[1] hascomplex(term::Base.Generator{UnitRange{Int64}, var"#11#12"})
@ PDEBase C:\Users\U1087908\.julia\packages\PDEBase\CpWUp\src\symbolic_utils.jl:292
[2] (::PDEBase.var"#125#133")(eq::Base.Generator{UnitRange{Int64}, var"#11#12"})
@ PDEBase C:\Users\U1087908\.julia\packages\PDEBase\CpWUp\src\make_pdesys_compatible.jl:100
[3] _any(f::PDEBase.var"#125#133", itr::Vector{Any}, ::Colon)
@ Base .\reduce.jl:1237
[4] any(f::Function, a::Vector{Any}; dims::Function)
@ Base .\reducedim.jl:992
[5] handle_complex(pdesys::PDESystem)
@ PDEBase C:\Users\U1087908\.julia\packages\PDEBase\CpWUp\src\make_pdesys_compatible.jl:100
[6] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
@ PDEBase C:\Users\U1087908\.julia\packages\PDEBase\CpWUp\src\symbolic_discretize.jl:10
[7] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…}; analytic::Nothing, kwargs::@Kwargs{})
@ PDEBase C:\Users\U1087908\.julia\packages\PDEBase\CpWUp\src\discretization_state.jl:57
[8] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
@ PDEBase C:\Users\U1087908\.julia\packages\PDEBase\CpWUp\src\discretization_state.jl:54
[9] macro expansion
@ .\timing.jl:581 [inlined]
[10] top-level scope
@ c:\Users\U1087908\OneDrive - CompanyName\Vscode\ChromatographyModeling\BreakthroughPredict.jl:315
@ Base .\reduce.jl:1237
my guess is it’s an issue with the syntax for the indexing, but I’m not sure how to troubleshoot that. Thank you for any help.