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!

4

Hey, I’m wondering if you can help me with another issue I’m having with MOL. I posted this issue yesterday here: Comparing Plot Outputs for Different Versions of Julia

I’ve come to realize that it’s an issue with one of the packages I use and not the version specifically as I’ve downloaded v-1.8.1 on a different workstation and using the current packages I get the same issue. So the code works properly with packages from before April but creates faults with the current packages.
Given that MOL is the newest and likely less used package in my code, I wonder if you could look into this and see if there is something I should be doing differently / If it’s an issue with a different package.

Thank you!