Error with WENO Scheme

Specific Domains Modelling & Simulations
1

I’m having trouble using the WENO Scheme for my system of PDEs. My system is given as

Screen Shot 2023-05-17 at 12.08.30 PM
Screen Shot 2023-05-17 at 12.08.30 PM884×158 12.9 KB

and I have the following code to solve the numerical solution for each equation:

using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets, DifferentialEquations, Plots

# Method of Manufactured Solutions: exact solution

# model parameters

Q = 0.5
p = 0.5
Dc = 0.05
mu = 6.0
Yu = 0.63
ku = 4.0
X_inf = 10000.0
d = 0.25
a = 0.4
Y_lambda = 0.63
eta = Y_lambda/X_inf
kd = 0.4
flux(S,L) = mu*((S/(p-2*L))/(ku+(S/(p-2*L))))*L 

# Parameters, variables, and derivatives
@parameters t x
@variables S(..), W(..), L(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2

# 1D PDE and boundary conditions
eq  = [Dt(S(t, x)) ~ (-Q*(p-2*L(t,x))*Dx(S(t,x)) + Q*S(t,x)*(-2)*Dx((L(t,x))))/(p-2*L(t,x))^2 + (Dc*(p-2*L(t,x))*Dxx(S(t,x)) - Dc*(-2)*Dxx(L(t,x))*(S(t,x)))/(p-2*L(t,x))^2 - 2*flux(S(t,x),L(t,x)) - (mu*S(t,x)*W(t,x))/(Yu*(ku*(p-2*L(t,x))+S(t,x))),
     Dt(W(t, x)) ~ (-Q*(p-2*L(t,x))*Dx(W(t,x)) + Q*W(t,x)*(-2)*Dx((L(t,x))))/(p-2*L(t,x))^2 + (Dc*(p-2*L(t,x))*Dxx(W(t,x)) - Dc*(-2)*Dxx(L(t,x))*(W(t,x)))/(p-2*L(t,x))^2 + 2*X_inf*d*L(t,x) - 2*a*W(t,x) + (mu*S(t,x)*W(t,x))/(ku*(p-2*L(t,x))+S(t,x)),
     Dt(L(t, x)) ~ eta*flux(S(t,x),L(t,x)) - (kd+d)L(t,x) + (a*W(t,x))/(X_inf)]

bcs = [S(0, x) ~ 25,
       S(t, 0) ~ 25,
       Dx(S(t,1)) ~ 0.0,
        W(0, x) ~ 0,
       W(t, 0) ~ 0.0,
       Dx(W(t,1)) ~ 0.0,
        L(0, x) ~ 0.0025]

# Space and time domains
domains = [t ∈ Interval(0.0, 4.0),
           x ∈ Interval(0.0, 1.0)]

# PDE system
@named pdesys = PDESystem(eq, bcs, domains, [t, x], [S(t, x), W(t,x), L(t,x)])

# Method of lines discretization
dx = 0.01
order = 2
discretization = MOLFiniteDifference([x => dx], t, advection_scheme = WENOScheme()

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

# Solve ODE problem
sol = solve(prob,Rosenbrock23(), saveat=0.4)

# Plot results and compare with exact solution
discrete_x = sol[x]
discrete_t = sol[t]

solS = sol[S(t,x)]
solW = sol[W(t,x)]
solL = sol[L(t,x)]

p1 = plot()
for i in 1:length(discrete_t)
    plot!(discrete_x, solS[i, :], title="Numerical S", label="t=$(discrete_t[i])")
end
p2 = plot()
for i in 1:length(discrete_t)
	plot!(discrete_x, solW[i,:], title="Numerical W", label="t=$(discrete_t[i])")
end
p3 = plot()
for i in 1:length(discrete_t)
    plot!(discrete_x, solL[i, :], title="Numerical L", label="t=$(discrete_t[i])")
end


plot(p1, p2, p3, legend=:bottomleft)

However, when I try this I get the error:

LoadError: Scheme WENO applied to S(t, x) in direction of x at point CartesianIndex(1,) is not defined.
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] function_scheme(F::FunctionalScheme{typeof(MethodOfLines.weno_f), Vector{Nothing}, Vector{Nothing}, Vector{Float64}}, II::CartesianIndex{1}, s::MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}, bs::Vector{Any}, jx::Tuple{Int64, SymbolicUtils.BasicSymbolic{Real}}, u::SymbolicUtils.BasicSymbolic{Real}, ufunc::MethodOfLines.var"#central_ufunc#213"{MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}})
    @ MethodOfLines ~/.julia/packages/MethodOfLines/Egrbc/src/discretization/schemes/function_scheme/function_scheme.jl:33
  [3] (::MethodOfLines.var"#212#215"{SymbolicUtils.BasicSymbolic{Real}, FunctionalScheme{typeof(MethodOfLines.weno_f), Vector{Nothing}, Vector{Nothing}, Vector{Float64}}, CartesianIndex{1}, MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}, Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}, MethodOfLines.var"#central_ufunc#213"{MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}}})(x::SymbolicUtils.BasicSymbolic{Real})
    @ MethodOfLines ./none:0
  [4] iterate
    @ ./generator.jl:47 [inlined]
  [5] collect(itr::Base.Generator{Vector{Any}, MethodOfLines.var"#212#215"{SymbolicUtils.BasicSymbolic{Real}, FunctionalScheme{typeof(MethodOfLines.weno_f), Vector{Nothing}, Vector{Nothing}, Vector{Float64}}, CartesianIndex{1}, MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}, Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}, MethodOfLines.var"#central_ufunc#213"{MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}}}})
    @ Base ./array.jl:782
  [6] (::MethodOfLines.var"#211#214"{FunctionalScheme{typeof(MethodOfLines.weno_f), Vector{Nothing}, Vector{Nothing}, Vector{Float64}}, CartesianIndex{1}, MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}, Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}, MethodOfLines.var"#central_ufunc#213"{MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}}})(u::SymbolicUtils.BasicSymbolic{Real})
    @ MethodOfLines ./none:0
  [7] iterate
    @ ./generator.jl:47 [inlined]
  [8] collect(itr::Base.Generator{Vector{Any}, MethodOfLines.var"#211#214"{FunctionalScheme{typeof(MethodOfLines.weno_f), Vector{Nothing}, Vector{Nothing}, Vector{Float64}}, CartesianIndex{1}, MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}, Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}, MethodOfLines.var"#central_ufunc#213"{MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}}}})
    @ Base ./array.jl:782
  [9] generate_advection_rules
    @ ~/.julia/packages/MethodOfLines/Egrbc/src/discretization/schemes/function_scheme/function_scheme.jl:55 [inlined]
 [10] generate_finite_difference_rules(II::CartesianIndex{1}, s::MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}, depvars::Vector{Any}, pde::Equation, derivweights::MethodOfLines.DifferentialDiscretizer{Pair{SymbolicUtils.BasicSymbolic{Real}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Float64, StaticArraysCore.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}, FunctionalScheme{typeof(MethodOfLines.weno_f), Vector{Nothing}, Vector{Nothing}, Vector{Float64}}}, bmap::Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, indexmap::Dict{SymbolicUtils.BasicSymbolic{Real}, Int64})
    @ MethodOfLines ~/.julia/packages/MethodOfLines/Egrbc/src/discretization/generate_finite_difference_rules.jl:58
 [11] discretize_equation_at_point(II::CartesianIndex{1}, s::MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}, depvars::Vector{Any}, pde::Equation, derivweights::MethodOfLines.DifferentialDiscretizer{Pair{SymbolicUtils.BasicSymbolic{Real}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Float64, StaticArraysCore.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}, FunctionalScheme{typeof(MethodOfLines.weno_f), Vector{Nothing}, Vector{Nothing}, Vector{Float64}}}, bcmap::Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, eqvar::SymbolicUtils.BasicSymbolic{Real}, indexmap::Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}, boundaryvalfuncs::Vector{Any})
    @ MethodOfLines ~/.julia/packages/MethodOfLines/Egrbc/src/scalar_discretization.jl:34
 [12] #439
    @ ~/.julia/packages/MethodOfLines/Egrbc/src/scalar_discretization.jl:24 [inlined]
 [13] iterate
    @ ./generator.jl:47 [inlined]
 [14] _collect(c::CartesianIndices{1, Tuple{UnitRange{Int64}}}, itr::Base.Generator{CartesianIndices{1, Tuple{UnitRange{Int64}}}, MethodOfLines.var"#439#441"{Equation, SymbolicUtils.BasicSymbolic{Real}, Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, Vector{Any}, MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}, MethodOfLines.DifferentialDiscretizer{Pair{SymbolicUtils.BasicSymbolic{Real}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Float64, StaticArraysCore.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}, FunctionalScheme{typeof(MethodOfLines.weno_f), Vector{Nothing}, Vector{Nothing}, Vector{Float64}}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}, Vector{Any}}}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})
    @ Base ./array.jl:802
 [15] collect_similar(cont::CartesianIndices{1, Tuple{UnitRange{Int64}}}, itr::Base.Generator{CartesianIndices{1, Tuple{UnitRange{Int64}}}, MethodOfLines.var"#439#441"{Equation, SymbolicUtils.BasicSymbolic{Real}, Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, Vector{Any}, MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}, MethodOfLines.DifferentialDiscretizer{Pair{SymbolicUtils.BasicSymbolic{Real}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Float64, StaticArraysCore.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}, FunctionalScheme{typeof(MethodOfLines.weno_f), Vector{Nothing}, Vector{Nothing}, Vector{Float64}}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}, Vector{Any}}})
    @ Base ./array.jl:711
 [16] map(f::Function, A::CartesianIndices{1, Tuple{UnitRange{Int64}}})
    @ Base ./abstractarray.jl:3261
 [17] discretize_equation!(disc_state::PDEBase.EquationState, pde::Equation, interiormap::MethodOfLines.InteriorMap, eqvar::SymbolicUtils.BasicSymbolic{Real}, bcmap::Dict{SymbolicUtils.BasicSymbolic{SymbolicUtils.FnType{Tuple, Real}}, Dict{SymbolicUtils.BasicSymbolic{Real}, Vector{Any}}}, depvars::Vector{Any}, s::MethodOfLines.DiscreteSpace{1, 3, MethodOfLines.CenterAlignedGrid}, derivweights::MethodOfLines.DifferentialDiscretizer{Pair{SymbolicUtils.BasicSymbolic{Real}, Vector{Int64}}, Dict{Function, MethodOfLines.DerivativeOperator{Float64, Nothing, false, Float64, StaticArraysCore.SVector{3, Float64}, S2, S3, Nothing, Nothing} where {S2, S3}}, FunctionalScheme{typeof(MethodOfLines.weno_f), Vector{Nothing}, Vector{Nothing}, Vector{Float64}}}, indexmap::Dict{SymbolicUtils.BasicSymbolic{Real}, Int64}, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
    @ MethodOfLines ~/.julia/packages/MethodOfLines/Egrbc/src/scalar_discretization.jl:23
 [18] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
    @ PDEBase ~/.julia/packages/PDEBase/Dg0WI/src/symbolic_discretize.jl:79
 [19] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}; analytic::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ PDEBase ~/.julia/packages/PDEBase/Dg0WI/src/discretization_state.jl:58
 [20] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization})
    @ PDEBase ~/.julia/packages/PDEBase/Dg0WI/src/discretization_state.jl:55
 [21] top-level scope
    @ ~/Desktop/Emma_Apr5.jl:54
 [22] include(fname::String)
    @ Base.MainInclude ./client.jl:478
 [23] top-level scope
    @ REPL[5]:1
in expression starting at /Users/emmabottomley/Desktop/Emma_Apr5.jl:54

Wondering if anyone knows how to solve this problem?
Thank you!

2

Hey, I’m the maintainer, sorry about this, could you post the full stacktrace with the latest version? Also, we’d love to have your system in PDESystemLibrary.jl, what’s it for? It would be amazing if you’d submit a PR with a brief description.

3

Hey thanks! I’m doing my master’s right now and I wrote this code for my research in biofilms. I’m hoping to publish it so I don’t really want it to be part of a library until then!