Hi all, I’m currently working on implementing a relatively simple reaction-diffusion-convection model of a 2D hydrogen flame. The PDESystem is created successfully, but the discretization step using methodoflines.jl fails with an error that I do not understand:
ERROR: ArgumentError: collection must be non-empty
Stacktrace:
[1] first
@ ./abstractarray.jl:387 [inlined]
[2] MethodOfLines.BoundaryHandler(bcs::Vector{Equation}, s::MethodOfLines.DiscreteSpace{2, 1, MethodOfLines.CenterAlignedGrid}, depvar_ops::Vector{Sym{SymbolicUtils.FnType{Tuple, Real}, Nothing}}, tspan::Tuple{Float64, Float64}, derivweights::MethodOfLines.DifferentialDiscretizer{Any, Dict{Function, DiffEqOperators.DerivativeOperator{Float64, Nothing, false, Float64, StaticArrays.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}})
@ MethodOfLines ~/.julia/packages/MethodOfLines/2RBAY/src/bcs/boundary_types.jl:142
[3] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
@ MethodOfLines ~/.julia/packages/MethodOfLines/2RBAY/src/discretization/MOL_discretization.jl:41
[4] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
@ MethodOfLines ~/.julia/packages/MethodOfLines/2RBAY/src/discretization/MOL_discretization.jl:120
[5] top-level scope
@ ./timing.jl:210 [inlined]
[6] top-level scope
@ ~/Museum/uni/msc5/UQDDM/project/scripts/hydrogenFlame.jl:0
If you know what’s going on here I’d love to know. In case it’s useful, I’ve added the model definition below (the error occurs at the very last line)
Summary
## Define variables and derivative operators
# the two spatial dimensions and time
@parameters x y t
# the mass fractions MF (respectively H₂, O₂, and H₂O), temperature,
# and the source terms of the mass fractions and temperature
@variables MF[1:3](..) T(..) Smf[1:3](..) ST(..)
Dt = Differential(t)
Dx = Differential(x)
Dy = Differential(y)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
# laplace operator
∇²(u) = Dxx(u) + Dyy(u)
# gradient operator
grad(u) = [Dx(u), Dy(u)]
## Define domains
xmin = ymin = 0.0 # cm
xmax = 1.8 # cm
ymax = 0.9 # cm
tmin = 0.0
tmax = 0.06 # s
domains = [
x ∈ Interval(xmin, xmax),
y ∈ Interval(ymin, ymax),
t ∈ Interval(tmin, tmax)
]
## Define parameters
# diffusivity coefficient (for temperature and mass fractions)
κ = 2.0 #* cm^2/s
# constant (divergence-free) velocity field
U = [50, 0] #.* cm/s
# density of the mixture
ρ = 1.39e-3 #* g/cm^3
# molecular weights (respectively H₂, O₂, and H₂O)
W = [2.016, 31.9, 18] #.* g/mol
# stoichiometric coefficients
ν = [2, 1, 2]
# heat of the reaction
Q = 9800 #* K
# universal gas constant
R = 8.314472 * 100 #* 1e-2 J/mol/K -> because J = Nm = 100 N cm
# y coordinates of the inlet area on the left side of the domain
inlety = (0.3, 0.6) # mm
# TODO try out different values (systematically)
# pre-exponential factor in source term
A = 5.5e11 # dimensionless
# activation energy in source term
E = 5.5e13 * 100 # 1e-2 J/mol -> J = Nm = 100 N cm
## Define the model
# diffusion and convection terms for the mass fractions
# TODO check if you can vectorize these equations by joining them
eqs = [
Dt( MF[i](x,y,t) ) ~ κ * ∇²( MF[i](x,y,t) ) - U' * grad(MF[i](x,y,t))
for i in 1:3
]
# diffusion and convection terms for the temperature
push!(eqs, Dt(T(x,y,t)) ~ κ * ∇²( T(x,y,t) ) - U' * grad(T(x,y,t)))
# TODO continue here please
function atInlet(x,y,t)
return (inlety[1] < y) * (y < inlety[2])
end
bcs = [
#
# initial conditions
#
T(x, y, 0) ~ 300.0, # K; around 26 C
# domain is empty at the start
MF[1](x,y,0) ~ 0.0,
MF[2](x,y,0) ~ 0.0,
MF[3](x,y,0) ~ 0.0,
#
# boundary conditions
#
# left side
T(xmin,y,t) ~ atInlet(xmin,y,t) * 950 + (1-atInlet(xmin,y,t)) * 300, # K
MF[1](xmin,y,t) ~ atInlet(xmin,y,t) * 0.0282, # mass fraction
MF[2](xmin,y,t) ~ atInlet(xmin,y,t) * 0.2259,
MF[3](xmin,y,t) ~ 0.0,
# right side
Dt( T(xmax,y,t) ) ~ 0.0, # K/s
Dt( MF[1](xmax,y,t) ) ~ 0.0, # mass fraction/s
Dt( MF[2](xmax,y,t) ) ~ 0.0,
Dt( MF[3](xmax,y,t) ) ~ 0.0,
# top
Dt( T(x,ymax,t) ) ~ 0.0, # K/s
Dt( MF[1](x,ymax,t) ) ~ 0.0, # mass fraction/s
Dt( MF[2](x,ymax,t) ) ~ 0.0,
Dt( MF[3](x,ymax,t) ) ~ 0.0,
# bottom
Dt( T(x,ymin,t) ) ~ 0.0, # K/s
Dt( MF[1](x,ymin,t) ) ~ 0.0, # mass fraction/s
Dt( MF[2](x,ymin,t) ) ~ 0.0,
Dt( MF[3](x,ymin,t) ) ~ 0.0,
]
@named pdesys = PDESystem(eqs, bcs, domains, [x,y,t],[T(x, y, t)])
## Discretize the system
N = 32
dx = 1/N
dy = 1/N
order = 2
discretization = MOLFiniteDifference(
[x=>dx, y=>dy], t, approx_order=order
)
# this creates an ODEProblem or a NonlinearProblem, depending on the system
problem = discretize(pdesys, discretization)
. I would just try v0.2.
If I could make a recommendation, I think it would be good to make a note of the still-early-in-development stage in the documentation in that case. At the moment it seems like the MOL package is stable and that can make for an interesting surprise.
I understand that you’d have confidence in the package though, the Brusselator model seems quite complex.