MethodOfLines.jl icon indicating copy to clipboard operation
MethodOfLines.jl copied to clipboard

Crashes when trying to solve the BBM equation

Open ranocha opened this issue 2 years ago • 6 comments

Trying to solve the BBM equation yields

julia> using Pkg; Pkg.activate(temp=true); Pkg.add(["MethodOfLines", "ModelingToolkit", "DomainSets"]); using MethodOfLines, ModelingToolkit, DomainSets
  Activating new project at `/tmp/jl_4Y3LmJ`
    Updating registry at `~/.julia/registries/General.toml`
   Resolving package versions...
    Updating `/tmp/jl_4Y3LmJ/Project.toml`
  [5b8099bc] + DomainSets v0.5.11
  [94925ecb] + MethodOfLines v0.3.2
  [961ee093] + ModelingToolkit v8.18.1
...

julia> using ModelingToolkit, MethodOfLines, DomainSets

julia> @parameters t x
2-element Vector{Num}:
 t
 x

julia> @variables u(..)
1-element Vector{Symbolics.CallWithMetadata{SymbolicUtils.FnType{Tuple, Real}, Base.ImmutableDict{DataType, Any}}}:
 u⋆

julia> Dt = Differential(t)
(::Differential) (generic function with 2 methods)

julia> Dx = Differential(x)
(::Differential) (generic function with 2 methods)

julia> Dxx = Differential(x)^2
Differential(x) ∘ Differential(x)

julia> Dtxx = Dt * Dxx
Differential(t) ∘ Differential(x) ∘ Differential(x)

julia> eq  = Dt(u(t, x)) + Dx(u(t, x)) - (Dt*Dxx)(u(t, x)) ~ 0
Differential(t)(u(t, x)) + Differential(x)(u(t, x)) - Differential(t)(Differential(x)(Differential(x)(u(t, x)))) ~ 0

julia> ic_bc = [u(0, x) ~ cospi(x), u(t, 0) ~ u(t, 1)]
2-element Vector{Equation}:
 u(0, x) ~ cospi(x)
 u(t, 0) ~ u(t, 1)

julia> domains = [t ∈ Interval(0.0, 1.0),
                  x ∈ Interval(0.0, 2.0)]
2-element Vector{Symbolics.VarDomainPairing}:
 Symbolics.VarDomainPairing(t, 0.0..1.0)
 Symbolics.VarDomainPairing(x, 0.0..2.0)

julia> @named pdesys = PDESystem(eq, ic_bc, domains, [t, x], [u(t, x)])
PDESystem
Equations: Equation[Differential(t)(u(t, x)) + Differential(x)(u(t, x)) - Differential(t)(Differential(x)(Differential(x)(u(t, x)))) ~ 0]
Boundary Conditions: Equation[u(0, x) ~ cospi(x), u(t, 0) ~ u(t, 1)]
Domain: Symbolics.VarDomainPairing[Symbolics.VarDomainPairing(t, 0.0..1.0), Symbolics.VarDomainPairing(x, 0.0..2.0)]
Dependent Variables: Num[u(t, x)]
Independent Variables: Num[t, x]
Parameters: SciMLBase.NullParameters()
Default Parameter ValuesDict{Any, Any}()

julia> discretization = MOLFiniteDifference([x => 0.1], t)
MOLFiniteDifference{MethodOfLines.CenterAlignedGrid}(Pair{Num, Float64}[x => 0.1], t, 2, MethodOfLines.UpwindScheme(), MethodOfLines.CenterAlignedGrid(), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}())

julia> prob = discretize(pdesys,discretization)
ERROR: MethodError: no method matching deleteat!(::OffsetArrays.OffsetVector{Num, Vector{Num}}, ::UnitRange{Int64})
Closest candidates are:
  deleteat!(::Vector, ::UnitRange{<:Integer}) at ~/Software/julia-1.7.3/share/julia/base/array.jl:1427
  deleteat!(::Vector, ::AbstractVector) at ~/Software/julia-1.7.3/share/julia/base/array.jl:1469
  deleteat!(::Vector, ::Any) at ~/Software/julia-1.7.3/share/julia/base/array.jl:1468
  ...
Stacktrace:
 [1] filter!(f::Base.var"#321#322"{Set{Num}}, a::OffsetArrays.OffsetVector{Num, Vector{Num}})
   @ Base ./array.jl:2536
 [2] _shrink!(shrinker!::Function, v::OffsetArrays.OffsetVector{Num, Vector{Num}}, itrs::Tuple{Vector{Num}})
   @ Base ./array.jl:2622
 [3] setdiff!(v::OffsetArrays.OffsetVector{Num, Vector{Num}}, itrs::Vector{Num})
   @ Base ./array.jl:2626
 [4] error_analysis(sys::ODESystem, e::ModelingToolkit.ExtraVariablesSystemException)
   @ MethodOfLines ~/.julia/packages/MethodOfLines/5wITU/src/error_analysis.jl:13
 [5] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/5wITU/src/MOL_discretization.jl:157
 [6] top-level scope
   @ REPL[16]:1

caused by: ExtraVariablesSystemException: The system is unbalanced. There are 41 highest order derivative variables and 20 equations.
More variables than equations, here are the potential extra variable(s):
 Differential(t)(99.99999999999999u[1](t) + 99.99999999999999u[3](t) - 199.99999999999997u[2](t))
 Differential(t)(99.99999999999999u[2](t) + 99.99999999999999u[4](t) - 199.99999999999997u[3](t))
 Differential(t)(99.99999999999999u[3](t) + 99.99999999999999u[5](t) - 199.99999999999997u[4](t))
 Differential(t)(99.99999999999999u[4](t) + 99.99999999999999u[6](t) - 199.99999999999997u[5](t))
 Differential(t)(99.99999999999999u[5](t) + 99.99999999999999u[7](t) - 199.99999999999997u[6](t))
 Differential(t)(99.99999999999999u[6](t) + 99.99999999999999u[8](t) - 199.99999999999997u[7](t))
 Differential(t)(99.99999999999999u[7](t) + 99.99999999999999u[9](t) - 199.99999999999997u[8](t))
 Differential(t)(99.99999999999999u[8](t) + 99.99999999999999u[10](t) - 199.99999999999997u[9](t))
 Differential(t)(99.99999999999999u[9](t) + 99.99999999999999u[11](t) - 199.99999999999997u[10](t))
 Differential(t)(99.99999999999999u[10](t) + 99.99999999999999u[12](t) - 199.99999999999997u[11](t))
 Differential(t)(99.99999999999999u[11](t) + 99.99999999999999u[13](t) - 199.99999999999997u[12](t))
 Differential(t)(99.99999999999999u[12](t) + 99.99999999999999u[14](t) - 199.99999999999997u[13](t))
 Differential(t)(99.99999999999999u[13](t) + 99.99999999999999u[15](t) - 199.99999999999997u[14](t))
 Differential(t)(99.99999999999999u[14](t) + 99.99999999999999u[16](t) - 199.99999999999997u[15](t))
 Differential(t)(99.99999999999999u[15](t) + 99.99999999999999u[17](t) - 199.99999999999997u[16](t))
 Differential(t)(99.99999999999999u[16](t) + 99.99999999999999u[18](t) - 199.99999999999997u[17](t))
 Differential(t)(99.99999999999999u[17](t) + 99.99999999999999u[19](t) - 199.99999999999997u[18](t))
 Differential(t)(99.99999999999999u[18](t) + 99.99999999999999u[20](t) - 199.99999999999997u[19](t))
 Differential(t)(99.99999999999999u[19](t) + 99.99999999999999u[21](t) - 199.99999999999997u[20](t))
 Differential(t)(399.99999999999994u[19](t) + 199.99999999999997u[21](t) - 99.99999999999999u[18](t) - 499.99999999999994u[20](t))
 u[2](t)
 99.99999999999999u[1](t) + 99.99999999999999u[3](t) - 199.99999999999997u[2](t)
 u[3](t)
 99.99999999999999u[2](t) + 99.99999999999999u[4](t) - 199.99999999999997u[3](t)
 u[4](t)
 99.99999999999999u[3](t) + 99.99999999999999u[5](t) - 199.99999999999997u[4](t)
 u[5](t)
 99.99999999999999u[4](t) + 99.99999999999999u[6](t) - 199.99999999999997u[5](t)
 u[6](t)
 99.99999999999999u[5](t) + 99.99999999999999u[7](t) - 199.99999999999997u[6](t)
 u[7](t)
 99.99999999999999u[6](t) + 99.99999999999999u[8](t) - 199.99999999999997u[7](t)
 u[8](t)
 99.99999999999999u[7](t) + 99.99999999999999u[9](t) - 199.99999999999997u[8](t)
 u[9](t)
 99.99999999999999u[8](t) + 99.99999999999999u[10](t) - 199.99999999999997u[9](t)
 u[10](t)
 99.99999999999999u[9](t) + 99.99999999999999u[11](t) - 199.99999999999997u[10](t)
 u[11](t)
 99.99999999999999u[10](t) + 99.99999999999999u[12](t) - 199.99999999999997u[11](t)
 u[12](t)
 99.99999999999999u[11](t) + 99.99999999999999u[13](t) - 199.99999999999997u[12](t)
 u[13](t)
 99.99999999999999u[12](t) + 99.99999999999999u[14](t) - 199.99999999999997u[13](t)
 u[14](t)
 99.99999999999999u[13](t) + 99.99999999999999u[15](t) - 199.99999999999997u[14](t)
 u[15](t)
 99.99999999999999u[14](t) + 99.99999999999999u[16](t) - 199.99999999999997u[15](t)
 u[16](t)
 99.99999999999999u[15](t) + 99.99999999999999u[17](t) - 199.99999999999997u[16](t)
 u[17](t)
 99.99999999999999u[16](t) + 99.99999999999999u[18](t) - 199.99999999999997u[17](t)
 u[18](t)
 99.99999999999999u[17](t) + 99.99999999999999u[19](t) - 199.99999999999997u[18](t)
 u[19](t)
 99.99999999999999u[18](t) + 99.99999999999999u[20](t) - 199.99999999999997u[19](t)
 u[20](t)
 99.99999999999999u[19](t) + 99.99999999999999u[21](t) - 199.99999999999997u[20](t)
 u[21](t)
 399.99999999999994u[19](t) + 199.99999999999997u[21](t) - 99.99999999999999u[18](t) - 499.99999999999994u[20](t)
 u[1](t)
Stacktrace:
 [1] error_reporting(state::TearingState{ODESystem}, bad_idxs::Vector{Int64}, n_highest_vars::Int64, iseqs::Bool)
   @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/mQkI6/src/structural_transformation/utils.jl:35
 [2] check_consistency(state::TearingState{ODESystem})
   @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/mQkI6/src/structural_transformation/utils.jl:65
 [3] structural_simplify(sys::ODESystem, io::Nothing; simplify::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/mQkI6/src/systems/abstractsystem.jl:969
 [4] structural_simplify (repeats 2 times)
   @ ~/.julia/packages/ModelingToolkit/mQkI6/src/systems/abstractsystem.jl:960 [inlined]
 [5] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/5wITU/src/MOL_discretization.jl:150
 [6] top-level scope
   @ REPL[16]:1

ranocha avatar Jul 28 '22 09:07 ranocha

Explicit mixed derivatives are unfortunately not yet supported, though may be solved by defining an auxiliary variable like Dxx_u(t, x) ~ Dxx(u(t,x)) and doing Dt(Dxx_u(t, x)) instead.

xtalax avatar Jul 28 '22 16:07 xtalax

Like this?

julia> using Pkg; Pkg.activate(temp=true); Pkg.add(["MethodOfLines", "ModelingToolkit", "DomainSets"]); using MethodOfLines, ModelingToolkit, DomainSets
  Activating new project at `/tmp/jl_1KdU3o`
    Updating registry at `~/.julia/registries/General.toml`
   Resolving package versions...
   Installed AbstractAlgebra ─ v0.27.1
    Updating `/tmp/jl_1KdU3o/Project.toml`
  [5b8099bc] + DomainSets v0.5.11
  [94925ecb] + MethodOfLines v0.3.2
  [961ee093] + ModelingToolkit v8.18.1
...

julia> using ModelingToolkit, MethodOfLines, DomainSets

julia> @parameters t x
2-element Vector{Num}:
 t
 x

julia> @variables u(..)
1-element Vector{Symbolics.CallWithMetadata{SymbolicUtils.FnType{Tuple, Real}, Base.ImmutableDict{DataType, Any}}}:
 u⋆

julia> Dt = Differential(t)
(::Differential) (generic function with 2 methods)

julia> Dx = Differential(x)
(::Differential) (generic function with 2 methods)

julia> Dxx = Differential(x)^2
Differential(x) ∘ Differential(x)

julia> @variables Dxx_u(..)
1-element Vector{Symbolics.CallWithMetadata{SymbolicUtils.FnType{Tuple, Real}, Base.ImmutableDict{DataType, Any}}}:
 Dxx_u⋆

julia> eq  = [Dt(u(t, x)) + Dx(u(t, x)) - Dt(Dxx_u(t, x)) ~ 0, Dxx_u(t, x) ~ Dxx(u(t, x))]
2-element Vector{Equation}:
 Differential(t)(u(t, x)) + Differential(x)(u(t, x)) - Differential(t)(Dxx_u(t, x)) ~ 0
 Dxx_u(t, x) ~ Differential(x)(Differential(x)(u(t, x)))

julia> ic_bc = [u(0, x) ~ cospi(x), u(t, 0) ~ u(t, 1)]
2-element Vector{Equation}:
 u(0, x) ~ cospi(x)
 u(t, 0) ~ u(t, 1)

julia> domains = [t ∈ Interval(0.0, 1.0),
                  x ∈ Interval(0.0, 2.0)]
2-element Vector{Symbolics.VarDomainPairing}:
 Symbolics.VarDomainPairing(t, 0.0..1.0)
 Symbolics.VarDomainPairing(x, 0.0..2.0)

julia> @named pdesys = PDESystem(eq, ic_bc, domains, [t, x], [u(t, x), Dxx_u(t, x)])
PDESystem
Equations: Equation[Differential(t)(u(t, x)) + Differential(x)(u(t, x)) - Differential(t)(Dxx_u(t, x)) ~ 0, Dxx_u(t, x) ~ Differential(x)(Differential(x)(u(t, x)))]
Boundary Conditions: Equation[u(0, x) ~ cospi(x), u(t, 0) ~ u(t, 1)]
Domain: Symbolics.VarDomainPairing[Symbolics.VarDomainPairing(t, 0.0..1.0), Symbolics.VarDomainPairing(x, 0.0..2.0)]
Dependent Variables: Num[u(t, x), Dxx_u(t, x)]
Independent Variables: Num[t, x]
Parameters: SciMLBase.NullParameters()
Default Parameter ValuesDict{Any, Any}()

julia> discretization = MOLFiniteDifference([x => 0.1], t)
MOLFiniteDifference{MethodOfLines.CenterAlignedGrid}(Pair{Num, Float64}[x => 0.1], t, 2, MethodOfLines.UpwindScheme(), MethodOfLines.CenterAlignedGrid(), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}())

julia> prob = discretize(pdesys,discretization)
ERROR: MethodError: no method matching deleteat!(::OffsetArrays.OffsetVector{Num, Vector{Num}}, ::UnitRange{Int64})
Closest candidates are:
  deleteat!(::Vector, ::UnitRange{<:Integer}) at ~/Software/julia-1.7.3/share/julia/base/array.jl:1427
  deleteat!(::Vector, ::AbstractVector) at ~/Software/julia-1.7.3/share/julia/base/array.jl:1469
  deleteat!(::Vector, ::Any) at ~/Software/julia-1.7.3/share/julia/base/array.jl:1468
  ...
Stacktrace:
 [1] filter!(f::Base.var"#321#322"{Set{Num}}, a::OffsetArrays.OffsetVector{Num, Vector{Num}})
   @ Base ./array.jl:2536
 [2] _shrink!(shrinker!::Function, v::OffsetArrays.OffsetVector{Num, Vector{Num}}, itrs::Tuple{Vector{Num}})
   @ Base ./array.jl:2622
 [3] setdiff!(v::OffsetArrays.OffsetVector{Num, Vector{Num}}, itrs::Vector{Num})
   @ Base ./array.jl:2626
 [4] error_analysis(sys::ODESystem, e::ModelingToolkit.ExtraVariablesSystemException)
   @ MethodOfLines ~/.julia/packages/MethodOfLines/5wITU/src/error_analysis.jl:13
 [5] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/5wITU/src/MOL_discretization.jl:157
 [6] top-level scope
   @ REPL[23]:1

caused by: ExtraVariablesSystemException: The system is unbalanced. There are 42 highest order derivative variables and 41 equations.
More variables than equations, here are the potential extra variable(s):
 Differential(t)(Dxx_u[1](t))
 Differential(t)(Dxx_u[2](t))
 Differential(t)(Dxx_u[3](t))
 Differential(t)(Dxx_u[4](t))
 Differential(t)(Dxx_u[5](t))
 Differential(t)(Dxx_u[6](t))
 Differential(t)(Dxx_u[7](t))
 Differential(t)(Dxx_u[8](t))
 Differential(t)(Dxx_u[9](t))
 Differential(t)(Dxx_u[10](t))
 Differential(t)(Dxx_u[11](t))
 Differential(t)(Dxx_u[12](t))
 Differential(t)(Dxx_u[13](t))
 Differential(t)(Dxx_u[14](t))
 Differential(t)(Dxx_u[15](t))
 Differential(t)(Dxx_u[16](t))
 Differential(t)(Dxx_u[17](t))
 Differential(t)(Dxx_u[18](t))
 Differential(t)(Dxx_u[19](t))
 Differential(t)(Dxx_u[20](t))
 Differential(t)(Dxx_u[21](t))
 u[1](t)
 Dxx_u[1](t)
 u[2](t)
 Dxx_u[2](t)
 u[3](t)
 Dxx_u[3](t)
 u[4](t)
 Dxx_u[4](t)
 u[5](t)
 Dxx_u[5](t)
 u[6](t)
 Dxx_u[6](t)
 u[7](t)
 Dxx_u[7](t)
 u[8](t)
 Dxx_u[8](t)
 u[9](t)
 Dxx_u[9](t)
 u[10](t)
 Dxx_u[10](t)
 u[11](t)
 Dxx_u[11](t)
 u[12](t)
 Dxx_u[12](t)
 u[13](t)
 Dxx_u[13](t)
 u[14](t)
 Dxx_u[14](t)
 u[15](t)
 Dxx_u[15](t)
 u[16](t)
 Dxx_u[16](t)
 u[17](t)
 Dxx_u[17](t)
 u[18](t)
 Dxx_u[18](t)
 u[19](t)
 Dxx_u[19](t)
 u[20](t)
 Dxx_u[20](t)
 u[21](t)
 Dxx_u[21](t)
Stacktrace:
 [1] error_reporting(state::TearingState{ODESystem}, bad_idxs::Vector{Int64}, n_highest_vars::Int64, iseqs::Bool)
   @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/mQkI6/src/structural_transformation/utils.jl:35
 [2] check_consistency(state::TearingState{ODESystem})
   @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/mQkI6/src/structural_transformation/utils.jl:65
 [3] structural_simplify(sys::ODESystem, io::Nothing; simplify::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/mQkI6/src/systems/abstractsystem.jl:969
 [4] structural_simplify (repeats 2 times)
   @ ~/.julia/packages/ModelingToolkit/mQkI6/src/systems/abstractsystem.jl:960 [inlined]
 [5] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/5wITU/src/MOL_discretization.jl:150
 [6] top-level scope
   @ REPL[23]:1

ranocha avatar Jul 28 '22 17:07 ranocha

This variable also needs an initial condition, but then I would expect this to work.

xtalax avatar Jul 29 '22 13:07 xtalax

Like this?

julia> using Pkg; Pkg.activate(temp=true); Pkg.add(["MethodOfLines", "ModelingToolkit", "DomainSets"]); using MethodOfLines, ModelingToolkit, DomainSets
  Activating new project at `/tmp/jl_PDoHJg`
    Updating registry at `~/.julia/registries/General.toml`
   Resolving package versions...
   Installed Symbolics ─ v4.10.3
    Updating `/tmp/jl_PDoHJg/Project.toml`
  [5b8099bc] + DomainSets v0.5.11
  [94925ecb] + MethodOfLines v0.3.2
  [961ee093] + ModelingToolkit v8.18.1
...

julia> using ModelingToolkit, MethodOfLines, DomainSets

julia> @parameters t x
2-element Vector{Num}:
 t
 x

julia> @variables u(..)
1-element Vector{Symbolics.CallWithMetadata{SymbolicUtils.FnType{Tuple, Real}, Base.ImmutableDict{DataType, Any}}}:
 u⋆

julia> Dt = Differential(t)
(::Differential) (generic function with 2 methods)

julia> Dx = Differential(x)
(::Differential) (generic function with 2 methods)

julia> Dxx = Differential(x)^2
Differential(x) ∘ Differential(x)

julia> @variables Dxx_u(..)
1-element Vector{Symbolics.CallWithMetadata{SymbolicUtils.FnType{Tuple, Real}, Base.ImmutableDict{DataType, Any}}}:
 Dxx_u⋆

julia> eq  = [Dt(u(t, x)) + Dx(u(t, x)) - Dt(Dxx_u(t, x)) ~ 0, Dxx_u(t, x) ~ Dxx(u(t, x))]
2-element Vector{Equation}:
 Differential(t)(u(t, x)) + Differential(x)(u(t, x)) - Differential(t)(Dxx_u(t, x)) ~ 0
 Dxx_u(t, x) ~ Differential(x)(Differential(x)(u(t, x)))

julia> ic_bc = [u(0, x) ~ cospi(x), u(t, 0) ~ u(t, 1), Dxx_u(0, x) ~ -π^2 * cospi(x)]
3-element Vector{Equation}:
 u(0, x) ~ cospi(x)
 u(t, 0) ~ u(t, 1)
 Dxx_u(0, x) ~ -9.869604401089358cospi(x)

julia> domains = [t ∈ Interval(0.0, 1.0),
                  x ∈ Interval(0.0, 2.0)]
2-element Vector{Symbolics.VarDomainPairing}:
 Symbolics.VarDomainPairing(t, 0.0..1.0)
 Symbolics.VarDomainPairing(x, 0.0..2.0)

julia> @named pdesys = PDESystem(eq, ic_bc, domains, [t, x], [u(t, x), Dxx_u(t, x)])
PDESystem
Equations: Equation[Differential(t)(u(t, x)) + Differential(x)(u(t, x)) - Differential(t)(Dxx_u(t, x)) ~ 0, Dxx_u(t, x) ~ Differential(x)(Differential(x)(u(t, x)))]
Boundary Conditions: Equation[u(0, x) ~ cospi(x), u(t, 0) ~ u(t, 1), Dxx_u(0, x) ~ -9.869604401089358cospi(x)]
Domain: Symbolics.VarDomainPairing[Symbolics.VarDomainPairing(t, 0.0..1.0), Symbolics.VarDomainPairing(x, 0.0..2.0)]
Dependent Variables: Num[u(t, x), Dxx_u(t, x)]
Independent Variables: Num[t, x]
Parameters: SciMLBase.NullParameters()
Default Parameter ValuesDict{Any, Any}()

julia> discretization = MOLFiniteDifference([x => 0.1], t)
MOLFiniteDifference{MethodOfLines.CenterAlignedGrid}(Pair{Num, Float64}[x => 0.1], t, 2, MethodOfLines.UpwindScheme(), MethodOfLines.CenterAlignedGrid(), Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}())

julia> prob = discretize(pdesys,discretization)
ERROR: MethodError: no method matching deleteat!(::OffsetArrays.OffsetVector{Num, Vector{Num}}, ::UnitRange{Int64})
Closest candidates are:
  deleteat!(::Vector, ::UnitRange{<:Integer}) at ~/Software/julia-1.7.3/share/julia/base/array.jl:1427
  deleteat!(::Vector, ::AbstractVector) at ~/Software/julia-1.7.3/share/julia/base/array.jl:1469
  deleteat!(::Vector, ::Any) at ~/Software/julia-1.7.3/share/julia/base/array.jl:1468
  ...
Stacktrace:
 [1] filter!(f::Base.var"#321#322"{Set{Num}}, a::OffsetArrays.OffsetVector{Num, Vector{Num}})
   @ Base ./array.jl:2536
 [2] _shrink!(shrinker!::Function, v::OffsetArrays.OffsetVector{Num, Vector{Num}}, itrs::Tuple{Vector{Num}})
   @ Base ./array.jl:2622
 [3] setdiff!(v::OffsetArrays.OffsetVector{Num, Vector{Num}}, itrs::Vector{Num})
   @ Base ./array.jl:2626
 [4] error_analysis(sys::ODESystem, e::ModelingToolkit.ExtraVariablesSystemException)
   @ MethodOfLines ~/.julia/packages/MethodOfLines/5wITU/src/error_analysis.jl:13
 [5] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/5wITU/src/MOL_discretization.jl:157
 [6] top-level scope
   @ REPL[14]:1

caused by: ExtraVariablesSystemException: The system is unbalanced. There are 42 highest order derivative variables and 41 equations.
More variables than equations, here are the potential extra variable(s):
 Differential(t)(Dxx_u[1](t))
 Differential(t)(Dxx_u[2](t))
 Differential(t)(Dxx_u[3](t))
 Differential(t)(Dxx_u[4](t))
 Differential(t)(Dxx_u[5](t))
 Differential(t)(Dxx_u[6](t))
 Differential(t)(Dxx_u[7](t))
 Differential(t)(Dxx_u[8](t))
 Differential(t)(Dxx_u[9](t))
 Differential(t)(Dxx_u[10](t))
 Differential(t)(Dxx_u[11](t))
 Differential(t)(Dxx_u[12](t))
 Differential(t)(Dxx_u[13](t))
 Differential(t)(Dxx_u[14](t))
 Differential(t)(Dxx_u[15](t))
 Differential(t)(Dxx_u[16](t))
 Differential(t)(Dxx_u[17](t))
 Differential(t)(Dxx_u[18](t))
 Differential(t)(Dxx_u[19](t))
 Differential(t)(Dxx_u[20](t))
 Differential(t)(Dxx_u[21](t))
 u[1](t)
 Dxx_u[1](t)
 u[2](t)
 Dxx_u[2](t)
 u[3](t)
 Dxx_u[3](t)
 u[4](t)
 Dxx_u[4](t)
 u[5](t)
 Dxx_u[5](t)
 u[6](t)
 Dxx_u[6](t)
 u[7](t)
 Dxx_u[7](t)
 u[8](t)
 Dxx_u[8](t)
 u[9](t)
 Dxx_u[9](t)
 u[10](t)
 Dxx_u[10](t)
 u[11](t)
 Dxx_u[11](t)
 u[12](t)
 Dxx_u[12](t)
 u[13](t)
 Dxx_u[13](t)
 u[14](t)
 Dxx_u[14](t)
 u[15](t)
 Dxx_u[15](t)
 u[16](t)
 Dxx_u[16](t)
 u[17](t)
 Dxx_u[17](t)
 u[18](t)
 Dxx_u[18](t)
 u[19](t)
 Dxx_u[19](t)
 u[20](t)
 Dxx_u[20](t)
 u[21](t)
 Dxx_u[21](t)
Stacktrace:
 [1] error_reporting(state::TearingState{ODESystem}, bad_idxs::Vector{Int64}, n_highest_vars::Int64, iseqs::Bool)
   @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/mQkI6/src/structural_transformation/utils.jl:35
 [2] check_consistency(state::TearingState{ODESystem})
   @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/mQkI6/src/structural_transformation/utils.jl:65
 [3] structural_simplify(sys::ODESystem, io::Nothing; simplify::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/mQkI6/src/systems/abstractsystem.jl:969
 [4] structural_simplify (repeats 2 times)
   @ ~/.julia/packages/ModelingToolkit/mQkI6/src/systems/abstractsystem.jl:960 [inlined]
 [5] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/5wITU/src/MOL_discretization.jl:150
 [6] top-level scope
   @ REPL[14]:1

With another BC for Dxx_u(t, x) (continuing the REPL session above):

julia> ic_bc = [u(0, x) ~ cospi(x), u(t, 0) ~ u(t, 1), Dxx_u(0, x) ~ -π^2 * cospi(x), Dxx_u(t, 0) ~ Dxx_u(t, 1)]
4-element Vector{Equation}:
 u(0, x) ~ cospi(x)
 u(t, 0) ~ u(t, 1)
 Dxx_u(0, x) ~ -9.869604401089358cospi(x)
 Dxx_u(t, 0) ~ Dxx_u(t, 1)

julia> @named pdesys = PDESystem(eq, ic_bc, domains, [t, x], [u(t, x), Dxx_u(t, x)])
PDESystem
Equations: Equation[Differential(t)(u(t, x)) + Differential(x)(u(t, x)) - Differential(t)(Dxx_u(t, x)) ~ 0, Dxx_u(t, x) ~ Differential(x)(Differential(x)(u(t, x)))]
Boundary Conditions: Equation[u(0, x) ~ cospi(x), u(t, 0) ~ u(t, 1), Dxx_u(0, x) ~ -9.869604401089358cospi(x), Dxx_u(t, 0) ~ Dxx_u(t, 1)]
Domain: Symbolics.VarDomainPairing[Symbolics.VarDomainPairing(t, 0.0..1.0), Symbolics.VarDomainPairing(x, 0.0..2.0)]
Dependent Variables: Num[u(t, x), Dxx_u(t, x)]
Independent Variables: Num[t, x]
Parameters: SciMLBase.NullParameters()
Default Parameter ValuesDict{Any, Any}()

julia> prob = discretize(pdesys,discretization)
ERROR: MethodError: no method matching deleteat!(::OffsetArrays.OffsetVector{Num, Vector{Num}}, ::UnitRange{Int64})
Closest candidates are:
  deleteat!(::Vector, ::UnitRange{<:Integer}) at ~/Software/julia-1.7.3/share/julia/base/array.jl:1427
  deleteat!(::Vector, ::AbstractVector) at ~/Software/julia-1.7.3/share/julia/base/array.jl:1469
  deleteat!(::Vector, ::Any) at ~/Software/julia-1.7.3/share/julia/base/array.jl:1468
  ...
Stacktrace:
 [1] filter!(f::Base.var"#321#322"{Set{Num}}, a::OffsetArrays.OffsetVector{Num, Vector{Num}})
   @ Base ./array.jl:2536
 [2] _shrink!(shrinker!::Function, v::OffsetArrays.OffsetVector{Num, Vector{Num}}, itrs::Tuple{Vector{Num}})
   @ Base ./array.jl:2622
 [3] setdiff!(v::OffsetArrays.OffsetVector{Num, Vector{Num}}, itrs::Vector{Num})
   @ Base ./array.jl:2626
 [4] error_analysis(sys::ODESystem, e::ModelingToolkit.ExtraVariablesSystemException)
   @ MethodOfLines ~/.julia/packages/MethodOfLines/5wITU/src/error_analysis.jl:13
 [5] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/5wITU/src/MOL_discretization.jl:157
 [6] top-level scope
   @ REPL[17]:1

caused by: ExtraVariablesSystemException: The system is unbalanced. There are 42 highest order derivative variables and 41 equations.
More variables than equations, here are the potential extra variable(s):
 Differential(t)(Dxx_u[2](t))
 Differential(t)(Dxx_u[3](t))
 Differential(t)(Dxx_u[4](t))
 Differential(t)(Dxx_u[5](t))
 Differential(t)(Dxx_u[6](t))
 Differential(t)(Dxx_u[7](t))
 Differential(t)(Dxx_u[8](t))
 Differential(t)(Dxx_u[9](t))
 Differential(t)(Dxx_u[10](t))
 Differential(t)(Dxx_u[11](t))
 Differential(t)(Dxx_u[12](t))
 Differential(t)(Dxx_u[13](t))
 Differential(t)(Dxx_u[14](t))
 Differential(t)(Dxx_u[15](t))
 Differential(t)(Dxx_u[16](t))
 Differential(t)(Dxx_u[17](t))
 Differential(t)(Dxx_u[18](t))
 Differential(t)(Dxx_u[19](t))
 Differential(t)(Dxx_u[20](t))
 Differential(t)(Dxx_u[21](t))
 u[2](t)
 Dxx_u[2](t)
 u[3](t)
 Dxx_u[3](t)
 u[4](t)
 Dxx_u[4](t)
 u[5](t)
 Dxx_u[5](t)
 u[6](t)
 Dxx_u[6](t)
 u[7](t)
 Dxx_u[7](t)
 u[8](t)
 Dxx_u[8](t)
 u[9](t)
 Dxx_u[9](t)
 u[10](t)
 Dxx_u[10](t)
 u[11](t)
 Dxx_u[11](t)
 u[12](t)
 Dxx_u[12](t)
 u[13](t)
 Dxx_u[13](t)
 u[14](t)
 Dxx_u[14](t)
 u[15](t)
 Dxx_u[15](t)
 u[16](t)
 Dxx_u[16](t)
 u[17](t)
 Dxx_u[17](t)
 u[18](t)
 Dxx_u[18](t)
 u[19](t)
 Dxx_u[19](t)
 u[20](t)
 Dxx_u[20](t)
 u[21](t)
 Dxx_u[21](t)
Stacktrace:
 [1] error_reporting(state::TearingState{ODESystem}, bad_idxs::Vector{Int64}, n_highest_vars::Int64, iseqs::Bool)
   @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/mQkI6/src/structural_transformation/utils.jl:35
 [2] check_consistency(state::TearingState{ODESystem})
   @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/mQkI6/src/structural_transformation/utils.jl:65
 [3] structural_simplify(sys::ODESystem, io::Nothing; simplify::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ ModelingToolkit ~/.julia/packages/ModelingToolkit/mQkI6/src/systems/abstractsystem.jl:969
 [4] structural_simplify (repeats 2 times)
   @ ~/.julia/packages/ModelingToolkit/mQkI6/src/systems/abstractsystem.jl:960 [inlined]
 [5] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
   @ MethodOfLines ~/.julia/packages/MethodOfLines/5wITU/src/MOL_discretization.jl:150
 [6] top-level scope
   @ REPL[17]:1

ranocha avatar Jul 29 '22 13:07 ranocha

Ok I'm sorry about this, I don't know where this is coming from. I will take a closer look soon.

xtalax avatar Jul 30 '22 12:07 xtalax

No worries - you asked to crash MethodOfLines.jl some time ago so here I am 😉

ranocha avatar Jul 30 '22 12:07 ranocha