Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Crashes when trying to solve the BBM equation #145

Closed
ranocha opened this issue Jul 28, 2022 · 6 comments · Fixed by #168
Closed

Crashes when trying to solve the BBM equation #145

ranocha opened this issue Jul 28, 2022 · 6 comments · Fixed by #168
Labels
bug Something isn't working

Comments

@ranocha
Copy link
Member

ranocha commented Jul 28, 2022

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
@xtalax
Copy link
Member

xtalax commented Jul 28, 2022

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 xtalax mentioned this issue Jul 28, 2022
Closed
@ranocha
Copy link
Member Author

ranocha commented Jul 28, 2022

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

@xtalax
Copy link
Member

xtalax commented Jul 29, 2022

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

@ranocha
Copy link
Member Author

ranocha commented Jul 29, 2022

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

@xtalax
Copy link
Member

xtalax commented Jul 30, 2022

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

@xtalax xtalax added the bug Something isn't working label Jul 30, 2022
@ranocha
Copy link
Member Author

ranocha commented Jul 30, 2022 •
edited
Loading

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

@xtalax xtalax linked a pull request Oct 7, 2022 that will close this issue
System Transformation - Automatically make PDESystems compatible #168
Merged
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
bug Something isn't working
Projects
None yet
Milestone
No milestone
Development

Successfully merging a pull request may close this issue.

System Transformation - Automatically make PDESystems compatible xtalax/MethodOfLines.jl
2 participants
@xtalax @ranocha