# PDE discretization error in MethodOfLines.jl

**URL:** <https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566>\
**Category:** Modelling & Simulations\
**Created:** [March 27, 2022, 6:53pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566 "2022-03-27T18:53:50Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![RickDW](https://avatars.discourse-cdn.com/v4/letter/r/97f17d/32.png) [@RickDW](https://discourse.julialang.org/u/RickDW)\
**Post date:** [March 27, 2022, 6:53pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/1 "2022-03-27T18:53:50Z")

</div>

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:

```julia
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**
>
> ```julia
> ## 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)
> 
> ```

---

<div class="post-metadata">

**Author:** ![RickDW](https://avatars.discourse-cdn.com/v4/letter/r/97f17d/32.png) [@RickDW](https://discourse.julialang.org/u/RickDW)\
**Post date:** [March 27, 2022, 8:02pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/2 "2022-03-27T20:02:54Z")

</div>

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)

```julia
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**
>
> ```julia
> 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)
> 
> ```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 27, 2022, 8:34pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/3 "2022-03-27T20:34:14Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![RickDW](https://avatars.discourse-cdn.com/v4/letter/r/97f17d/32.png) [@RickDW](https://discourse.julialang.org/u/RickDW)\
**Post date:** [March 28, 2022, 1:25pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/4 "2022-03-28T13:25:08Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 28, 2022, 1:31pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/5 "2022-03-28T13:31:46Z")

</div>

> [@RickDW](#):
>
> Good call, thanks. I’m currently using v0.0.1. Do you think v0.2 has a fix for this?

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

---

<div class="post-metadata">

**Author:** ![RickDW](https://avatars.discourse-cdn.com/v4/letter/r/97f17d/32.png) [@RickDW](https://discourse.julialang.org/u/RickDW)\
**Post date:** [March 28, 2022, 1:33pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/6 "2022-03-28T13:33:49Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![RickDW](https://avatars.discourse-cdn.com/v4/letter/r/97f17d/32.png) [@RickDW](https://discourse.julialang.org/u/RickDW)\
**Post date:** [March 28, 2022, 1:53pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/7 "2022-03-28T13:53:39Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 28, 2022, 3:03pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/8 "2022-03-28T15:03:20Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![RickDW](https://avatars.discourse-cdn.com/v4/letter/r/97f17d/32.png) [@RickDW](https://discourse.julialang.org/u/RickDW)\
**Post date:** [March 28, 2022, 3:24pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/9 "2022-03-28T15:24:03Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 28, 2022, 3:32pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/10 "2022-03-28T15:32:44Z")

</div>

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 😓

---

<div class="post-metadata">

**Author:** ![RickDW](https://avatars.discourse-cdn.com/v4/letter/r/97f17d/32.png) [@RickDW](https://discourse.julialang.org/u/RickDW)\
**Post date:** [March 28, 2022, 3:37pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/11 "2022-03-28T15:37:38Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [March 29, 2022, 5:20pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/12 "2022-03-29T17:20:56Z")

</div>

And it’s fixed, next please!

---

<div class="post-metadata">

**Author:** ![mkretzschmar](https://avatars.discourse-cdn.com/v4/letter/m/ba8739/32.png) [@mkretzschmar](https://discourse.julialang.org/u/mkretzschmar)\
**Post date:** [August 1, 2023, 9:51am UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/13 "2023-08-01T09:51:00Z")

</div>

I am using the MethodOfLines to solve a system of 3 first order pde’s. It worked fine with a two variable function. Now I extended to 3 variables and got into problems with the discretization. The error message was “Differential w.r.t. multiple variables are not allowed”.  
This is how my equations are defined:

# Parameters, variables, and derivatives

```
@parameters t a b
@variables S(..) I(..) R(..) 
Dt = Differential(t) 
Da = Differential(a) 
    Db = Differential(b)
Ia = Integral(a in DomainSets.ClosedInterval(amin,amax))

# PDE and boundary conditions
eqs = [ 
Dt(S(t,a,b)) + Da(S(t,a,b)) + Db(frailty(a)*S(t,a,b)) ~ - β * S(t,a,b) * Ia(I(t,a,b))/N + alpha * R(t,a,b) - mu(a)*S(t,a,b) - vacc(a,5.0, 10.0)*S(t,a,b), 
Dt(I(t,a,b)) + Da(I(t,a,b)) + Db(frailty(a)*I(t,a,b)) ~ β * S(t,a,b) * Ia(I(t,a,b))/N - γ*I(t,a,b) - mu(a)*I(t,a,b), 
Dt(R(t,a,b)) + Da(R(t,a,b)) + Db(frailty(a)*R(t,a,b)) ~ γ * I(t,a,b) - alpha * R(t,a,b) - mu(a)*R(t,a,b)+ vacc(a,5.0, 10.0)*S(t,a,b)] 

S₀ = N * 0.99 #initial number susceptible
I₀ = N * 0.01 #initial number infected
R₀ = 0.0 #initial number recovered
min_age = 0.0 #minimum age
max_age = 100.0 #maximum age 
min_b = 0.0 #minimum frailty
max_b = 1.0 #maximum frailty 
bcs = [ 
	S(0.0,a,b) ~ S₀/(max_age - min_age) * theta(max_age - a) * theta(a - min_age) * theta(max_b -b) * theta(b - min_b), 
	I(0.0,a,b) ~ I₀/(max_age - min_age) * theta(max_age - a) * theta(a - min_age) * theta(max_b -b) * theta(b - min_b), 
	R(0.0,a,b) ~ 0.0, 
	S(t,0.0,b) ~ N/(max_age - min_age),  
	I(t,0.0,b) ~ 0.0, 
	R(t,0.0,b) ~ 0.0,
    S(t,a,0.0) ~ N/(max_age - min_age),  
	I(t,a,0.0) ~ 0.0, 
	R(t,a,0.0) ~ 0.0 
]

# Space, age, and frailty domains
domains = [t ∈ Interval(tmin,tmax), a ∈ Interval(amin,amax), b ∈ Interval(bmin,bmax)]	

# PDE system
@named SIRmodel = PDESystem(eqs, bcs, domains, [t, a, b], [S(t, a, b), I(t,a,b), R(t,a,b)]) 

# Method of lines discretization
da = 1.0 
db = 0.01
discretization = MOLFiniteDifference([a=>da, b=>db],t)

```

I got the error message when applying prob = discretize(SIRmodel, discretization). I don’t know how to interpret the error message. I would be grateful for any hints!

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [August 8, 2023, 5:59pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/14 "2023-08-08T17:59:20Z")

</div>

> [@mkretzschmar](#):
>
> `Db(frailty(a)*S(t,a,b))`

Can you rewrite these terms as `frailty(a)*Db(S(t,a,b))`?

---

<div class="post-metadata">

**Author:** ![mkretzschmar](https://avatars.discourse-cdn.com/v4/letter/m/ba8739/32.png) [@mkretzschmar](https://discourse.julialang.org/u/mkretzschmar)\
**Post date:** [August 16, 2023, 2:26pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/15 "2023-08-16T14:26:04Z")

</div>

Thanks for the suggestion. I rewrote the terms as you suggested. Now when I try to discretize it takes so long that I have not been able to finish, and I am not sure whether it gets stuck somewhere or is still working on it. I tried to reduce the grid resolution to make it faster, but it did not help. Eventually I would like to add a diffusion term to the equation, so it would get even more complex and I am wondering whether it can be done at all with the method of lines. If I make any progress, I will post it.

---

<div class="post-metadata">

**Author:** ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)\
**Post date:** [August 17, 2023, 12:11pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/16 "2023-08-17T12:11:06Z")

</div>

I expect it’s getting stuck in compilation. We are working on methods to make this much faster, but this may not be possible at present

---

<div class="post-metadata">

**Author:** ![mkretzschmar](https://avatars.discourse-cdn.com/v4/letter/m/ba8739/32.png) [@mkretzschmar](https://discourse.julialang.org/u/mkretzschmar)\
**Post date:** [August 17, 2023, 4:17pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/17 "2023-08-17T16:17:10Z")

</div>

After 2 days of calculating a discretization, I got the following message, and I have no idea what it means:

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true

Is there a place where error messages are explained? I could not find anything in the manuals.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [August 17, 2023, 4:23pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/18 "2023-08-17T16:23:20Z")

</div>

That’s not an error message. That’s just the print of a problem type.

---

<div class="post-metadata">

**Author:** ![mkretzschmar](https://avatars.discourse-cdn.com/v4/letter/m/ba8739/32.png) [@mkretzschmar](https://discourse.julialang.org/u/mkretzschmar)\
**Post date:** [August 17, 2023, 4:49pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/19 "2023-08-17T16:49:36Z")

</div>

I see, thanks!

---

<div class="post-metadata">

**Author:** ![mkretzschmar](https://avatars.discourse-cdn.com/v4/letter/m/ba8739/32.png) [@mkretzschmar](https://discourse.julialang.org/u/mkretzschmar)\
**Post date:** [September 1, 2023, 12:43pm UTC](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566/20 "2023-09-01T12:43:34Z")

</div>

To make it faster I have now drastically reduced the resolution of the discretization. Discretization seems to work.  
When solving, I now get the error message

Cannot `convert` an object of type SymbolicUtils.BasicSymbolic{Float64} to an object of type Float64

I don’t understand where the problem is.

This is the code:  
using DifferentialEquations  
using OrdinaryDiffEq  
using ModelingToolkit  
using MethodOfLines  
using DomainSets  
using DiffEqCallbacks  
using Interpolations  
using Plots  
using Distributions  
#using SpecialFunctions  
#to find the exponential integral of x  
#expinti(x)

```
tmin = 0.0
tmax = 50.0
amin = 0.0 
amax = 100.0
bmin = 0.0
bmax = 1.0

#Model parameters
β = 3.0 #transmission rate
γ = 1.0 # recovery rate
alpha = 0.05 #waning rate
N = 100.0 #population size
k = 10.0
lambda = 90.0
 

function mu(age) #mortality function
#	0.01
   pdf(Weibull(k,lambda),age)/(1-cdf(Weibull(k, lambda),age))
    end
 @register_symbolic mu(age)

#Heaviside function 
theta(a) = (a > 0) ? 1 : 0 
@register_symbolic theta(a)

#Vaccination function
function vacc(age, alower, aupper)
    coverage = 0.9 #vaccination coverage
    parambeta = 2.0 #parameter betadistribution
    coverage * pdf(Beta(parambeta),(age - alower)/(aupper - alower))*theta(aupper - age) * theta(age - alower)
end
@register_symbolic vacc(age, alower, aupper)

# Parameters, variables, and derivatives
@parameters t a b
@variables S(..) I(..) R(..) 
Dt = Differential(t)
Da = Differential(a) 
Db = Differential(b)
Ia = Integral(a in DomainSets.ClosedInterval(amin,amax))
Ib = Integral(b in DomainSets.ClosedInterval(bmin,bmax))

# PDE and boundary conditions
eqs = [ 
Dt(S(t,a,b)) + Da(S(t,a,b)) + Db(S(t,a,b)) ~ - β * S(t,a,b) * Ia(Ib(I(t,a,b)))/N + alpha * R(t,a,b) - mu(a)*S(t,a,b), 
Dt(I(t,a,b)) + Da(I(t,a,b)) + Db(I(t,a,b)) ~ β * S(t,a,b) * Ia(Ib(I(t,a,b)))/N - γ*I(t,a,b) - mu(a)*I(t,a,b), 
Dt(R(t,a,b)) + Da(R(t,a,b)) + Db(R(t,a,b)) ~ γ * I(t,a,b) - alpha * R(t,a,b) - mu(a)*R(t,a,b)] 

S₀ = N * 0.99 #initial number susceptible
I₀ = N * 0.01 #initial number infected
R₀ = 0.0 #initial number recovered
min_age = 0.0 #minimum age
max_age = 100.0 #maximum age 
min_b = 0.0 #minimum frailty
max_b = 1.0 #maximum frailty 
bcs = [ 
	S(0.0,a,b) ~ S₀/(max_age - min_age) * theta(max_age - a) * theta(a - min_age) * theta(max_b - b) * theta(b - min_b), 
	I(0.0,a,b) ~ I₀/(max_age - min_age) * theta(max_age - a) * theta(a - min_age) * theta(max_b - b) * theta(b - min_b), 
	R(0.0,a,b) ~ 0.0, 
	S(t,0.0,b) ~ S₀/(max_age - min_age),  
	I(t,0.0,b) ~ I₀/(max_age - min_age), 
	R(t,0.0,b) ~ 0.0,
    S(t,a,0.0) ~ N/(max_age - min_age),  
	I(t,a,0.0) ~ 0.0, 
	R(t,a,0.0) ~ 0.0 
]

# Space, age, and frailty domains
domains = [t ∈ Interval(tmin,tmax),a ∈ Interval(amin,amax), b ∈ Interval(bmin,bmax)]	

# PDE system
@named SIRmodel = PDESystem(eqs, bcs, domains, [t,a,b], [S(t,a,b), I(t,a,b), R(t,a,b)]) 

# Method of lines discretization 
da = 10.0
db = 0.1
discretization = MOLFiniteDifference([a=>da,b=>db],t)

# Convert the PDE problem into an ODE problem
prob = discretize(SIRmodel, discretization) 

# Solve ODE problem
dt = tmax/100.0
sol_pde = solve(prob, Tsit5(), saveat = dt)

```

[Next page](https://discourse.julialang.org/t/pde-discretization-error-in-methodoflines-jl/78566.md?page=2)