WENOScheme seems to require support

[hurak]

I am trying to solve a simple linear advection equation. I am trying to explore the performance of the WENOScheme as described at https://docs.sciml.ai/MethodOfLines/stable/advection_schemes/. Since the proposed way of calling the function with the keyword argument epsilon does not work, as I report here at GitHub in another issue, I tried to omit the epsilon and call it just without any arguments: discretization = MOLFiniteDifference([x => Δx], t, advection_scheme = WENOScheme()). This seems to pass but then if I use it to discrete the problem and then solve it using some ODE solver, I get error message that "This solver is not able to use mass matrices" (see below). And besides the default Tsit5 I also tried those referenced on the above webpage as "Strong-Stability-Preserving (SSP) solver" such as SSPRK22 etc to no avail.

julia> sol = solve(prob, Tsit5())
ERROR: This solver is not able to use mass matrices. For compatible solvers see https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] __init(prob::ODEProblem{…}, alg::Tsit5{…}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}; saveat::Tuple{}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_on::Bool, save_start::Bool, save_end::Nothing, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, dtmin::Float64, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{…}, abstol::Nothing, reltol::Nothing, qmin::Rational{…}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{…}, controller::Nothing, fullnormalize::Bool, failfactor::Int64, maxiters::Int64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), internalopnorm::typeof(LinearAlgebra.opnorm), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), progress_id::Symbol, userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias::ODEAliasSpecifier, initializealg::OrdinaryDiffEqCore.DefaultInit, kwargs::@Kwargs{})
    @ OrdinaryDiffEqCore ~/.julia/packages/OrdinaryDiffEqCore/nLyB1/src/solve.jl:95
  [3] __init (repeats 2 times)
    @ ~/.julia/packages/OrdinaryDiffEqCore/nLyB1/src/solve.jl:11 [inlined]
  [4] #__solve#62
    @ ~/.julia/packages/OrdinaryDiffEqCore/nLyB1/src/solve.jl:6 [inlined]
  [5] __solve
    @ ~/.julia/packages/OrdinaryDiffEqCore/nLyB1/src/solve.jl:1 [inlined]
  [6] #solve_call#31
    @ ~/.julia/packages/DiffEqBase/yEe9v/src/solve.jl:127 [inlined]
  [7] solve_call
    @ ~/.julia/packages/DiffEqBase/yEe9v/src/solve.jl:84 [inlined]
  [8] #solve_up#40
    @ ~/.julia/packages/DiffEqBase/yEe9v/src/solve.jl:683 [inlined]
  [9] solve_up
    @ ~/.julia/packages/DiffEqBase/yEe9v/src/solve.jl:660 [inlined]
 [10] #solve#38
    @ ~/.julia/packages/DiffEqBase/yEe9v/src/solve.jl:555 [inlined]
 [11] solve(prob::ODEProblem{…}, args::Tsit5{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/yEe9v/src/solve.jl:545
 [12] top-level scope
    @ REPL[3]:1
Some type information was truncated. Use `show(err)` to see complete types.

For reproducibility I include the full code

using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets

@independent_variables t x
@variables T(..)
Dₜ = Differential(t)
Dₓ = Differential(x)

eq  = Dₜ(T(t, x)) ~ -Dₓ(T(t, x))

T₀ = 0.0

bcs = [T(0, x) ~ T₀+exp(-((x-2.5)^2)/0.5),
       T(t, 0) ~ T₀+exp(-((-t-2.5)^2)/0.5)]

domains = [t ∈ Interval(0.0, 7.0),
           x ∈ Interval(0.0, 5.0)]

@named pdesys = PDESystem(eq, bcs, domains, [t, x], [T(t, x)])

Δx = 0.2

discretization = MOLFiniteDifference([x => Δx], t, advection_scheme = WENOScheme())

prob = discretize(pdesys,discretization)

sol = solve(prob, Tsit5())
#sol = solve(prob, SSPRK22())