Method Error when using Method of Lines to solve PDEs

Specific Domains Modelling & Simulations
, , , ,
1

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.

2

The issue stems from improper use of generator expressions when building the equation array, i.e.:

  push!(eqs, Dt(εp * c(r,z,t)[j]) + ... for j in 1:JJ)

This creates a generator object rather than expanding into individual equations. When PDESystem constructor tries to process the equations (specifically checking for complex numbers in PDEBase.jl), it receives a generator object instead of actual symbolic equations, causing the method error. Instead, just expand it yourself:

  for j in 1:JJ
      push!(eqs, equation)  # Correct - adds individual equations
  end

Or collect the generator first:

  append!(eqs, [equation for j in 1:JJ])  # Also correct

so you just get a normal array.

3

Thank you so much for the response! I’ve managed to get it working, but I did have a couple questions regarding MOL/MTK specifically.

I moved away from the use of generator expressions, since MOL doesn’t currently support calling Boundary Conditions not at the edges of the domain.

#This got me the equations I needed
rp_vals = [50, 60, 70]
equationvec = Equation[]
for j in 1:J
    eqn1 = ep*Dt(c(r,t)[j]) ~ D*(Drr(c(r,t)[j]) + (2/r * Dr(c(r,t)[j])))
    global equationvec = push!(equationvec, eqn1)
end

#This caused an error because rp_vals[1] isn't on the boundary of the domain for r
boundaryvec = Equation[]
for j in 1:J
    bc1 = Dr(c(0,t)[j]) ~ 0
    global boundaryvec = push!(boundaryvec, bc1)
    bc2 = c(rp_vals[j],t)[j] ~ Cb
    global boundaryvec = push!(boundaryvec, bc2)
    bc3 = c(r, 0)[j] ~ 0
    global boundaryvec = push!(boundaryvec, bc3)
end

dom1 = r ∈ Interval(0, maximum(rp_vals))
dom2 = t ∈ Interval(0, 50)
domainvec = [dom1, dom2]

I ended up switching to defining each r as a separate variable like this:

PoreDiff1 = ep * Dt(c1(r1,z,t)) ~ De * ((2/r1 *Dr1(c1(r1,z,t))) + Drr1(c1(r1,z,t))) - Dt(q1(r1,z,t))

But that is a bit tedious to manually copy paste each equation for each r. On the docs it says that support for calling BCs in the middle of the domain is a feature in progress, is that something that will be added soon?

My other question is about compile times for MOL. I really like the straightforward syntax and readability. But right now, it’s slow to compile a system with 360 ODES (approx 7 mins). The docs mention JuliaSimCompiler, but the link there leads to nothing, and I can’t use the commercial software.

I also found this tutorial, which has a section on improving the process, but there is a warning and adding the following drastically increases compile time.

FullModel = discretize(ColModel, Discrete_ColModel; jac = true, sparse = true)

I’m guessing it’s still deprecated for now, but what’s the status on decreasing compile times and increasing scalability in general? Other than manually discretizing the system, are there any options that increase speed at high ODE count that currently work with MOL?

Thanks for your time!

4

We’re doing a bit of a rewrite to make this better, and it doesn’t require JuliaSimCompiler. More coming out on this soon.

1 Like