MethodOfLines.jl
MethodOfLines.jl copied to clipboard
Crashes when trying to solve the BBM equation
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
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.
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
This variable also needs an initial condition, but then I would expect this to work.
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
Ok I'm sorry about this, I don't know where this is coming from. I will take a closer look soon.
No worries - you asked to crash MethodOfLines.jl some time ago so here I am 😉