PDE discretization error

[RickDW]

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:
(I'm using method of lines v0.2.0 and modeling toolkit v8.8.5)

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)

## 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)

[xtalax]

Please retry with MethodOfLines.jl v0.2.0, there was a bug with array element variables that should be fixed. You also need to pass in MF[i](x, y, t) as dependent variables. I have not tried vector valued equation components like the output of grad(u), could you retry with a simple system using this if the error persists after the above changes?

[RickDW]

I've tried to run this in v0.2.0 and I unfortunately received the same error. I fixed the mistake of not passing dependent variables into the system definition and scalarized the equation component, and I now get a different error, so there's progress. I'll provide the script below that produces this error:

ERROR: AssertionError: Boundary condition Differential(t)(T(1.8, y, t)) ~ 0.0 is not on a boundary of the domain, or is not a valid boundary condition
Stacktrace:
 [1] BoundaryHandler(bcs::Vector{Equation}, s::MethodOfLines.DiscreteSpace{2, 4, 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/9MLPx/src/bcs/boundary_types.jl:204
 [2] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:48
 [3] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:146
 [4] top-level scope
   @ ~/Museum/uni/msc5/UQDDM/project/scripts/hydrogenFlame2.jl:135

I also ran into this error when I tried to implement a simple 2D heat diffusion model with non-homogeneous Dirichlet boundary conditions. If you think it could be helpful, I posted that error and the related script on discourse: https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/2?u=rickdw

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 (respectively H₂, O₂, and H₂O), temperature,
# and the source terms of the mass fractions and temperature
@variables YF(..) YO(..) YP(..) 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)


## 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, 0.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.0] #.* 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 and temperature
eqs = [
    Dt( YF(x,y,t) ) ~ κ * ∇²( YF(x,y,t) ) - U[1] * Dx( YF(x,y,t) ) - U[2] * Dy( YF(x,y,t) )
    Dt( YO(x,y,t) ) ~ κ * ∇²( YO(x,y,t) ) - U[1] * Dx( YO(x,y,t) ) - U[2] * Dy( YO(x,y,t) )
    Dt( YP(x,y,t) ) ~ κ * ∇²( YP(x,y,t) ) - U[1] * Dx( YP(x,y,t) ) - U[2] * Dy( YP(x,y,t) )
    Dt( T(x,y,t) )  ~ κ * ∇²( T(x,y,t) )  - U[1] * Dx( T(x,y,t) )  - U[2] * Dy( T(x,y,t) )
]

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
    YF(x,y,0) ~ 0.0,
    YO(x,y,0) ~ 0.0, 
    YP(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
    YF(xmin,y,t) ~ atInlet(xmin,y,t) * 0.0282, # mass fraction
    YO(xmin,y,t) ~ atInlet(xmin,y,t) * 0.2259,
    YP(xmin,y,t) ~ 0.0,

    # right side
    Dt( T(xmax,y,t) ) ~ 0.0, # K/s
    Dt( YF(xmax,y,t) ) ~ 0.0, # mass fraction/s
    Dt( YO(xmax,y,t) ) ~ 0.0,
    Dt( YP(xmax,y,t) ) ~ 0.0,

    # top
    Dt( T(x,ymax,t) ) ~ 0.0, # K/s
    Dt( YF(x,ymax,t) ) ~ 0.0, # mass fraction/s
    Dt( YO(x,ymax,t) ) ~ 0.0,
    Dt( YP(x,ymax,t) ) ~ 0.0,

    # bottom
    Dt( T(x,ymin,t) ) ~ 0.0, # K/s
    Dt( YF(x,ymin,t) ) ~ 0.0, # mass fraction/s
    Dt( YO(x,ymin,t) ) ~ 0.0,
    Dt( YP(x,ymin,t) ) ~ 0.0,
]

@named pdesys = PDESystem(eqs, bcs, domains, [x,y,t],[YF(x,y,t), YO(x,y,t), YP(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)

[xtalax]

Ok, that is a bug - thank you for finding it. I will take a look and try to get this working. We don't yet have interface tests for time derivatives in boundary conditions, I suspect the handling for that is the culprit.

[RickDW]

Does that mean Neumann boundary conditions are not officially supported yet, or is this just a regression?

[xtalax]

Neumann boundary conditions are supported, just time derivative neumann is untested and looks like it needs fixing

EDIT: the problem turned out to be something else, Dt bcs works

[RickDW]

Because of your comment I just realized that I've misinterpreted the boundary conditions of my model, and that I should've been using spatial derivatives at the boundary conditions. However, I get the same kind of error when I try that, so it's probably related to this bug.
To be specific, I replaced the Dt's with Dx's in the 'right side' boundary conditions, and I replaced the Dt's in the 'top' and 'bottom' BCs by Dy's. Here's the error which is pretty much identical to the other one:

ERROR: AssertionError: Boundary condition Differential(x)(T(1.8, y, t)) ~ 0.0 is not on a boundary of the domain, or is not a valid boundary condition
Stacktrace:
 [1] BoundaryHandler(bcs::Vector{Equation}, s::MethodOfLines.DiscreteSpace{2, 4, 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/9MLPx/src/bcs/boundary_types.jl:204
 [2] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:48
 [3] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:146
 [4] top-level scope
   @ ~/Museum/uni/msc5/UQDDM/project/scripts/hydrogenFlame2.jl:135

[xtalax]

Hmm, this is a serious bug, this has my full attention.

[xtalax]

@RickDW
This has revealed something I'm amazed that we've missed, a workaround is to specify a step size that goes in to the interval for its variable exactly,
e.g.

dx = (xmax-xmin)/N
dy = (ymax-ymin)/N

[RickDW]

It must have taken a while to figure out that's the issue, thank you so much for giving me a workaround! The discretization runs fine on my end as well, I'm just waiting to see what the actual simulation results look like but I doubt there'll be any more surprise there. Thanks again!

[RickDW]

If you don't mind me asking, is there any reason to think it's normal that this simulation with N=32 should run for over 20 minutes? I'm using this command to solve the ODE:

sol = solve(problem, TRBDF2(), saveat=0.01)

[xtalax]

Compile times are quite slow for the generated functions at the moment, and do not scale well with dimension and point count, when we have this they will be significantly (orders of magnitude) reduced and should not scale too much with these parameters

[RickDW]

Okay thanks for the heads up, wouldn't have expected compilation to be the issue here.

[xtalax]

Over 99% of the total time for solve is compile time right now, a manual looping semidiscretization which will be what is generated in the future compiles in thousandths of that total