PDE discretization error in MethodOfLines.jl

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’m not entirely sure whether this is related, but I ran into a different error while trying to fix this. The same code runs perfectly fine if I use xmax = ymax = 1.0 instead for the domain. Admittedly I’m not a PDE expert, but this does not seem like an invalid domain/boundary condition pairing to me. Does someone understand what is going on here?

(the script that produces this error can be found below)

ERROR: AssertionError: Boundary condition T(1.8, y, t) ~ -0.9510565162951538sin(6.283185307179586y) is not on a boundary of the domain, or is not a valid boundary condition
Stacktrace:
 [1] 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:189
 [2] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/2RBAY/src/discretization/MOL_discretization.jl:41
 [3] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/2RBAY/src/discretization/MOL_discretization.jl:120
 [4] top-level scope
   @ ~/Museum/uni/msc5/UQDDM/project/scripts/hydrogenFlame copy.jl:140
Summary
using ModelingToolkit
using MethodOfLines
using DomainSets
using DifferentialEquations
using GLMakie


## Define variables and derivative operators

# the two spatial dimensions and time
@parameters x y t
# the mass fractions MF, MO, and MP (respectively H₂, O₂, and H₂O), temperature
# @variables YF(..) YO(..) YP(..) T(..)
@variables T(..)
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( YF(x,y,t) ) ~ κ * ∇²( YF(x,y,t) ) #- U' * grad(MF(x,y,t))
    # Dt( YO(x,y,t) ) ~ κ * ∇²( YO(x,y,t) ) #- U' * grad(MO(x,y,t))
    # Dt( YP(x,y,t) ) ~ κ * ∇²( YP(x,y,t) ) #- U' * grad(MP(x,y,t))
    Dt( T(x,y,t) )  ~ κ * ∇²( T(x,y,t) )  #- U' * grad(T(x,y,t)))
]

# T0(x,y,t) = 300.0
T0(x, y, t) = sin(2π * x) * sin(2π * y)
# YF0(x,y,t) = 0.0
# YO0(x,y,t) = 0.0
# YP0(x,y,t) = 0.0

bcs = [
    #
    # initial conditions
    #

    T(x,y,0) ~ T0(x,y,0), # K; around 26 C
    # domain is empty at the start
    # YF(x,y,0) ~ YF0(x,y,0),
    # YO(x,y,0) ~ YO0(x,y,0), 
    # YP(x,y,0) ~ YP0(x,y,0),

    #
    # boundary conditions
    #

    # left side
    T(xmin,y,t) ~ T0(xmin,y,t), #, #atInlet(xmin,y,t) * 950 + (1-atInlet(xmin,y,t)) *  300, # K
    # YF(xmin,y,t) ~ YF0(xmin,y,t), #atInlet(xmin,y,t) * 0.0282, # mass fraction
    # YO(xmin,y,t) ~ YO0(xmin,y,t), #atInlet(xmin,y,t) * 0.2259,
    # YP(xmin,y,t) ~ YP0(xmin,y,t),

    # right side
    T(xmax,y,t) ~ T0(xmax,y,t), 
    # YF(xmax,y,t) ~ YF0(xmax,y,t),
    # YO(xmax,y,t) ~ YO0(xmax,y,t),
    # YP(xmax,y,t) ~ YP0(xmax,y,t),

    # top
    T(x,ymax,t) ~ T0(x,ymax,t), 
    # YF(x,ymax,t) ~ YF0(x,ymax,t),
    # YO(x,ymax,t) ~ YO0(x,ymax,t),
    # YP(x,ymax,t) ~ YP0(x,ymax,t),

    # bottom
    T(x,ymin,t) ~ T0(x,ymin,t), 
    # YF(x,ymin,t) ~ YF0(x,ymin,t),
    # YO(x,ymin,t) ~ YO0(x,ymin,t),
    # YP(x,ymin,t) ~ YP0(x,ymin,t),
]

@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 changed the title since this isn’t about ModelingToolkit.jl it seems. Is this with MethodOfLines v0.2? Can you open an issue there?

Good call, thanks. I’m currently using v0.0.1 v0.1.0. Do you think v0.2 has a fix for this?
I’ll create a github issue, maybe someone happens to know what’s going on.

0.0.1 effectively does nothing. It was registered early just to start the new package timer :sweat_smile:. I would just try v0.2.

Oops, sorry I meant v0.1.0 :sweat_smile:. I’ll give v0.2.0 a go in an hour and post an update afterwards. Thanks for the pointer :smiley:

The newer version raises the same error, that’s a shame. In the meantime I’ll go back to trying to implement this using DiffEqOperators. I tried that before but also had trouble implementing the boundary conditions, hopefully that goes more smoothly this time around.

Alex should have time to look at it today. MethodOfLines is still very early but it’s progressing quickly. DEO is more stalled.

Sounds good, thanks for all the quick replies :slight_smile: 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.

We had hubris because we got Brusselator and such working, so we thought it might be good enough now. Alex just asked on Slack yesterday to try and find ways to break it. It sounds like you won the prize :sweat:

Oof, I’m not sure whether I should be happy about this then :laughing: I understand that you’d have confidence in the package though, the Brusselator model seems quite complex.

And it’s fixed, next please!

2 Likes