search

SciML / MethodOfLines.jl

Automatic Finite Difference PDE solving with Julia SciML
https://docs.sciml.ai/MethodOfLines/stable/
MIT License
160 stars 28 forks source link

Advection Equation with time varying velocity: UndefVarError #194

Closed verflieger closed 1 year ago

verflieger commented 1 year ago

If i model a simple advection equation with time varying velocity, I run into an UndefVarError C not defined when simulating the discretized model. Maybe this is related to #189, but I am not sure. The problem does not occur if I use a constant velocity term or if I replace the zero boundary condition at the inlet with a nonzero bc.

The varying velocity is introduced via DataInterpolations.jl (the problem seems also to occur with Interpolations.jl). A MWE is

using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets, DataInterpolations

@parameters t z
@variables C(..)

Dt = Differential(t)
Dz = Differential(z)

v = 1:11;
Time = 0:0.1:1

tvip =  DataInterpolations.LinearInterpolation(vcat(v...),vcat(Time...))
vip1(t) = tvip(t);

@register vip1(t)

# define PDEs
# eq = [Dt(C(z,t)) ~ -1*Dz(C(z,t)), ]
eq = [Dt(C(z,t)) ~ -vip1(t)*Dz(C(z,t)), ]

bcs = [ C(0,t) ~ 0,         # boundary condition
        C(z,0) ~ 1,         # initial condition   
    ]

domains = [t ∈ Interval(0.0,1),
           z ∈ Interval(0.0,1)]

@named pdesys = PDESystem(eq,bcs,domains,[z,t],[
    C(z,t), 
    ], 
);

N = 5;
dz = 1/N;
order = 2;
discretization = MOLFiniteDifference([z => dz],t, approx_order=order, grid_align=center_align)

@time prob = discretize(pdesys,discretization);
solve(prob,Tsit5());

Here is the complete

Tracelog ``` ERROR: UndefVarError: C not defined Stacktrace: [1] macro expansion @ C:\Users\user\.julia\packages\SymbolicUtils\qulQp\src\code.jl:394 [inlined] [2] macro expansion @ C:\Users\user\.julia\packages\Symbolics\FGTCH\src\build_function.jl:519 [inlined] [3] macro expansion @ C:\Users\user\.julia\packages\SymbolicUtils\qulQp\src\code.jl:351 [inlined] [4] macro expansion @ C:\Users\user\.julia\packages\RuntimeGeneratedFunctions\6v5Gn\src\RuntimeGeneratedFunctions.jl:137 [inlined] [5] macro expansion @ .\none:0 [inlined] [6] generated_callfunc @ .\none:0 [inlined] [7] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xfd86900f, 0x3d09c415, 0x3a70905f, 0xc3f8b7e6, 0x2731f34e)})(::Vector{Float64}, ::Vector{Float64}, ::SciMLBase.NullParameters, ::Float64) @ RuntimeGeneratedFunctions C:\Users\user\.julia\packages\RuntimeGeneratedFunctions\6v5Gn\src\RuntimeGeneratedFunctions.jl:124 [8] f @ C:\Users\user\.julia\packages\ModelingToolkit\e7Li7\src\systems\diffeqs\abstractodesystem.jl:279 [inlined] [9] Void @ C:\Users\user\.julia\packages\SciMLBase\AwWUI\src\utils.jl:462 [inlined] [10] (::FunctionWrappers.CallWrapper{Nothing})(f::SciMLBase.Void{ModelingToolkit.var"#f#465"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf6985fb9, 0xb0545c9a, 0xb6be35b7, 0x1ca248be, 0x19479b71)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xfd86900f, 0x3d09c415, 0x3a70905f, 0xc3f8b7e6, 0x2731f34e)}}}, arg1::Vector{Float64}, arg2::Vector{Float64}, arg3::SciMLBase.NullParameters, arg4::Float64) @ FunctionWrappers C:\Users\user\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:65 [11] macro expansion @ C:\Users\user\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:137 [inlined] [12] do_ccall @ C:\Users\user\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:125 [inlined] [13] FunctionWrapper @ C:\Users\user\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:144 [inlined] [14] _call @ C:\Users\user\.julia\packages\FunctionWrappersWrappers\YyXoO\src\FunctionWrappersWrappers.jl:12 [inlined] [15] FunctionWrappersWrapper @ C:\Users\user\.julia\packages\FunctionWrappersWrappers\YyXoO\src\FunctionWrappersWrappers.jl:10 [inlined] [16] ODEFunction @ C:\Users\user\.julia\packages\SciMLBase\AwWUI\src\scimlfunctions.jl:1964 [inlined] [17] initialize!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Nothing, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, SciMLBase.NullParameters, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, SciMLBase.NullParameters, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, SciMLBase.NullParameters, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Any}, ModelingToolkit.var"#488#generated_observed#472"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, MethodOfLines.ScalarizedDiscretization}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, SciMLBase.NullParameters, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, SciMLBase.NullParameters, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, SciMLBase.NullParameters, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Any}, ModelingToolkit.var"#488#generated_observed#472"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats, Nothing}, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, SciMLBase.NullParameters, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, SciMLBase.NullParameters, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, SciMLBase.NullParameters, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Any}, ModelingToolkit.var"#488#generated_observed#472"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}) @ OrdinaryDiffEq C:\Users\user\.julia\packages\OrdinaryDiffEq\4qO6L\src\perform_step\low_order_rk_perform_step.jl:736 [18] __init(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, SciMLBase.NullParameters, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, SciMLBase.NullParameters, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, SciMLBase.NullParameters, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Any}, ModelingToolkit.var"#488#generated_observed#472"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, MethodOfLines.ScalarizedDiscretization}}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; 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::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, 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), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) @ OrdinaryDiffEq C:\Users\user\.julia\packages\OrdinaryDiffEq\4qO6L\src\solve.jl:493 [19] __init (repeats 5 times) @ C:\Users\user\.julia\packages\OrdinaryDiffEq\4qO6L\src\solve.jl:10 [inlined] [20] #__solve#562 @ C:\Users\user\.julia\packages\OrdinaryDiffEq\4qO6L\src\solve.jl:5 [inlined] [21] __solve @ C:\Users\user\.julia\packages\OrdinaryDiffEq\4qO6L\src\solve.jl:1 [inlined] [22] #solve_call#26 @ C:\Users\user\.julia\packages\DiffEqBase\g4OeQ\src\solve.jl:473 [inlined] [23] solve_call @ C:\Users\user\.julia\packages\DiffEqBase\g4OeQ\src\solve.jl:443 [inlined] [24] #solve_up#32 @ C:\Users\user\.julia\packages\DiffEqBase\g4OeQ\src\solve.jl:835 [inlined] [25] solve_up @ C:\Users\user\.julia\packages\DiffEqBase\g4OeQ\src\solve.jl:808 [inlined] [26] #solve#31 @ C:\Users\user\.julia\packages\DiffEqBase\g4OeQ\src\solve.jl:802 [inlined] [27] solve(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#465"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf6985fb9, 0xb0545c9a, 0xb6be35b7, 0x1ca248be, 0x19479b71)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xfd86900f, 0x3d09c415, 0x3a70905f, 0xc3f8b7e6, 0x2731f34e)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Any}, ModelingToolkit.var"#488#generated_observed#472"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, MethodOfLines.ScalarizedDiscretization}}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}) @ DiffEqBase C:\Users\user\.julia\packages\DiffEqBase\g4OeQ\src\solve.jl:792 [28] top-level scope @ c:\Users\user\KineticIso\simpleadvection.jl:42 ```

Package versions: Julia v1.8.1 [961ee093] ModelingToolkit v8.33.1 [94925ecb] MethodOfLines v0.6.1 [1dea7af3] OrdinaryDiffEq v6.31.2 [82cc6244] DataInterpolations v3.10.1

xtalax commented 1 year ago

Is this giving something like undef var error - C not defined?

verflieger commented 1 year ago

Yes exactly, sorry for not being exact there. Here is the complete

Tracelog ``` ERROR: UndefVarError: C not defined Stacktrace: [1] macro expansion @ C:\Users\user\.julia\packages\SymbolicUtils\qulQp\src\code.jl:394 [inlined] [2] macro expansion @ C:\Users\user\.julia\packages\Symbolics\FGTCH\src\build_function.jl:519 [inlined] [3] macro expansion @ C:\Users\user\.julia\packages\SymbolicUtils\qulQp\src\code.jl:351 [inlined] [4] macro expansion @ C:\Users\user\.julia\packages\RuntimeGeneratedFunctions\6v5Gn\src\RuntimeGeneratedFunctions.jl:137 [inlined] [5] macro expansion @ .\none:0 [inlined] [6] generated_callfunc @ .\none:0 [inlined] [7] (::RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xfd86900f, 0x3d09c415, 0x3a70905f, 0xc3f8b7e6, 0x2731f34e)})(::Vector{Float64}, ::Vector{Float64}, ::SciMLBase.NullParameters, ::Float64) @ RuntimeGeneratedFunctions C:\Users\user\.julia\packages\RuntimeGeneratedFunctions\6v5Gn\src\RuntimeGeneratedFunctions.jl:124 [8] f @ C:\Users\user\.julia\packages\ModelingToolkit\e7Li7\src\systems\diffeqs\abstractodesystem.jl:279 [inlined] [9] Void @ C:\Users\user\.julia\packages\SciMLBase\AwWUI\src\utils.jl:462 [inlined] [10] (::FunctionWrappers.CallWrapper{Nothing})(f::SciMLBase.Void{ModelingToolkit.var"#f#465"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf6985fb9, 0xb0545c9a, 0xb6be35b7, 0x1ca248be, 0x19479b71)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xfd86900f, 0x3d09c415, 0x3a70905f, 0xc3f8b7e6, 0x2731f34e)}}}, arg1::Vector{Float64}, arg2::Vector{Float64}, arg3::SciMLBase.NullParameters, arg4::Float64) @ FunctionWrappers C:\Users\user\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:65 [11] macro expansion @ C:\Users\user\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:137 [inlined] [12] do_ccall @ C:\Users\user\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:125 [inlined] [13] FunctionWrapper @ C:\Users\user\.julia\packages\FunctionWrappers\Q5cBx\src\FunctionWrappers.jl:144 [inlined] [14] _call @ C:\Users\user\.julia\packages\FunctionWrappersWrappers\YyXoO\src\FunctionWrappersWrappers.jl:12 [inlined] [15] FunctionWrappersWrapper @ C:\Users\user\.julia\packages\FunctionWrappersWrappers\YyXoO\src\FunctionWrappersWrappers.jl:10 [inlined] [16] ODEFunction @ C:\Users\user\.julia\packages\SciMLBase\AwWUI\src\scimlfunctions.jl:1964 [inlined] [17] initialize!(integrator::OrdinaryDiffEq.ODEIntegrator{Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, true, Vector{Float64}, Nothing, Float64, SciMLBase.NullParameters, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, SciMLBase.NullParameters, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, SciMLBase.NullParameters, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, SciMLBase.NullParameters, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Any}, ModelingToolkit.var"#488#generated_observed#472"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, MethodOfLines.ScalarizedDiscretization}}, Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.InterpolationData{ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, SciMLBase.NullParameters, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, SciMLBase.NullParameters, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, SciMLBase.NullParameters, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Any}, ModelingToolkit.var"#488#generated_observed#472"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}}, DiffEqBase.DEStats, Nothing}, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, SciMLBase.NullParameters, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, SciMLBase.NullParameters, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, SciMLBase.NullParameters, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Any}, ModelingToolkit.var"#488#generated_observed#472"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, OrdinaryDiffEq.DEOptions{Float64, Float64, Float64, Float64, PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Nothing, CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Int64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEq.DefaultInit}, cache::OrdinaryDiffEq.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, OrdinaryDiffEq.Tsit5ConstantCache{Float64, Float64}, typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}) @ OrdinaryDiffEq C:\Users\user\.julia\packages\OrdinaryDiffEq\4qO6L\src\perform_step\low_order_rk_perform_step.jl:736 [18] __init(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.AutoSpecialize, FunctionWrappersWrappers.FunctionWrappersWrapper{Tuple{FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, SciMLBase.NullParameters, Float64}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{Float64}, SciMLBase.NullParameters, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}, FunctionWrappers.FunctionWrapper{Nothing, Tuple{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}, SciMLBase.NullParameters, ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 1}}}}, false}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Any}, ModelingToolkit.var"#488#generated_observed#472"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, MethodOfLines.ScalarizedDiscretization}}, alg::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}; 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::Nothing, dtmax::Float64, force_dtmin::Bool, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmin::Rational{Int64}, qmax::Int64, qsteady_min::Int64, qsteady_max::Int64, beta1::Nothing, beta2::Nothing, qoldinit::Rational{Int64}, 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), userdata::Nothing, allow_extrapolation::Bool, initialize_integrator::Bool, alias_u0::Bool, alias_du0::Bool, initializealg::OrdinaryDiffEq.DefaultInit, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) @ OrdinaryDiffEq C:\Users\user\.julia\packages\OrdinaryDiffEq\4qO6L\src\solve.jl:493 [19] __init (repeats 5 times) @ C:\Users\user\.julia\packages\OrdinaryDiffEq\4qO6L\src\solve.jl:10 [inlined] [20] #__solve#562 @ C:\Users\user\.julia\packages\OrdinaryDiffEq\4qO6L\src\solve.jl:5 [inlined] [21] __solve @ C:\Users\user\.julia\packages\OrdinaryDiffEq\4qO6L\src\solve.jl:1 [inlined] [22] #solve_call#26 @ C:\Users\user\.julia\packages\DiffEqBase\g4OeQ\src\solve.jl:473 [inlined] [23] solve_call @ C:\Users\user\.julia\packages\DiffEqBase\g4OeQ\src\solve.jl:443 [inlined] [24] #solve_up#32 @ C:\Users\user\.julia\packages\DiffEqBase\g4OeQ\src\solve.jl:835 [inlined] [25] solve_up @ C:\Users\user\.julia\packages\DiffEqBase\g4OeQ\src\solve.jl:808 [inlined] [26] #solve#31 @ C:\Users\user\.julia\packages\DiffEqBase\g4OeQ\src\solve.jl:802 [inlined] [27] solve(prob::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, ODEFunction{true, SciMLBase.AutoSpecialize, ModelingToolkit.var"#f#465"{RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xf6985fb9, 0xb0545c9a, 0xb6be35b7, 0x1ca248be, 0x19479b71)}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :ˍ₋arg1, :ˍ₋arg2, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xfd86900f, 0x3d09c415, 0x3a70905f, 0xc3f8b7e6, 0x2731f34e)}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Vector{Symbol}, Symbol, Vector{Any}, ModelingToolkit.var"#488#generated_observed#472"{Bool, ODESystem, Dict{Any, Any}}, Nothing, ODESystem}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, MethodOfLines.MOLMetadata{Val{true}(), MethodOfLines.DiscreteSpace{1, 1, MethodOfLines.CenterAlignedGrid}, MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}, PDESystem, Base.RefValue{Any}, MethodOfLines.ScalarizedDiscretization}}, args::Tsit5{typeof(OrdinaryDiffEq.trivial_limiter!), typeof(OrdinaryDiffEq.trivial_limiter!), Static.False}) @ DiffEqBase C:\Users\user\.julia\packages\DiffEqBase\g4OeQ\src\solve.jl:792 [28] top-level scope @ c:\Users\user\KineticIso\simpleadvection.jl:42 ```

There is no error if I either

xtalax commented 1 year ago

Useful debugging fodder, thanks - this is a really weird one

verflieger commented 1 year ago

The problem is also not there if I add a diffusion term and a Neumann boundary condition

eq = [Dt(C(z,t)) ~ Dzz(C(z,t))-vip1(t)*Dz(C(z,t)), ]

bcs = [ C(0,t) ~ 0,         # boundary condition
        Dz(C(1,t)) ~ 0,
        C(z,0) ~ 1,         # initial condition   
    ]

I have some issues with debugging this error, so I can't really narrow it down further

xtalax commented 1 year ago

Smells like a bug in structural_simplify, can you try solving with the ODESysttem you get from symbolic_discrretize?

verflieger commented 1 year ago

I think that you are right. With symbolic_discretize I get a DAE problem with 6 states (before applying structural_simplify), see here and I can solve this with 

prob1 = DAEProblem(prob_symb,repeat([0.0],6),repeat([1.0],6), tspan; )
sol = solve(prob1, DFBDF() );

If I apply structural_simplify to this problem I get this output.

xtalax commented 1 year ago

Can you create an MWE generating the first system using only MTK, then we can see if this bug still shows in this case - and have something to give to @YingboMa. The Metadata shouldn't be relevant so don't worry about that

verflieger commented 1 year ago

At the moment, I cannot reproduce the problem with MTK alone Here is my MWE:

using ModelingToolkit, DifferentialEquations, DataInterpolations, OrdinaryDiffEq, DomainSets, MethodOfLines

@parameters t z
@variables C(..)

Dt = Differential(t)
Dz = Differential(z)

v = 1:11;
Time = 0:0.1:1

tvip =  DataInterpolations.LinearInterpolation(vcat(v...),vcat(Time...))
vip1(t) = tvip(t);
@register vip1(t)

# define PDEs
eq = [Dt(C(z,t)) ~ -vip1(t)*Dz(C(z,t)), ]

bcs = [ C(0,t) ~ 0,         # boundary condition
        C(z,0) ~ 1,         # initial condition   
    ]

domains = [t ∈ Interval(0.0,1),
           z ∈ Interval(0.0,1)]

@named pdesys = PDESystem(eq,bcs,domains,[z,t],[C(z,t), ], );

# Method of lines discretization
N = 5;
dz = 1/N;
order = 2;
discretization = MOLFiniteDifference([z => dz],t, approx_order=order, grid_align=center_align)
@time prob = discretize(pdesys,discretization);

symbmol, tspan = symbolic_discretize(pdesys,discretization);
prob1 = DAEProblem(symbmol,repeat([0.0],6),repeat([1.0],6), tspan; )
sol = solve(prob1, DFBDF() );

# define DAE using MTK
@variables (C(t))[1:6]

eqsmtk =[ 
    Dt(C[2]) + ifelse(vip1(t) > 0, (5.0(C)[2]-5.0(C)[1])*vip1(t), (5.0(C)[3] - 5.0(C)[2])*vip1(t)) ~ 0,
    Dt((C)[3]) + ifelse(vip1(t) > 0, (5.0(C)[3] - 5.0(C)[2])*vip1(t), (5.0(C)[4] - 5.0(C)[3])*vip1(t)) ~ 0,
    Dt((C)[4]) + ifelse(vip1(t) > 0, (5.0(C)[4] - 5.0(C)[3])*vip1(t), (5.0(C)[5] - 5.0(C)[4])*vip1(t)) ~ 0,
    Dt((C)[5]) + ifelse(vip1(t) > 0, (5.0(C)[5] - 5.0(C)[4])*vip1(t), (5.0(C)[6] - 5.0(C)[5])*vip1(t)) ~ 0,
    (C)[1] ~ 0, 
    (C)[6] ~ 10.0(C)[3] + 5.0(C)[5] + (C)[1] - 5.0(C)[2] - 10.0(C)[4],
    ]

@named symbmtk = ODESystem(eqsmtk,t,C,[],defaults=Dict(C[1]=>1.0, C[2]=>1.0,C[3]=>1.0,C[4]=>1.0,C[5]=>1.0,C[6]=>1.0,))

daeprob = DAEProblem(symbmtk,repeat([0.0],6),repeat([1.0],6), tspan; )
soldae = solve(daeprob, DFBDF() );

err = maximum(abs.(vcat((soldae.u - sol.u)...)))
equations(symbmol) == equations(symbmtk)
equations(structural_simplify(symbmol)) == equations(structural_simplify(symbmtk))

simpsysmol = structural_simplify(symbmol)
simpsysmtk = structural_simplify(symbmtk)

prob2 = ODEProblem(simpsysmtk,Pair[], tspan; )
solve(prob2,Rodas3())

Here I have both versions, the MOL and the MTK implementation. The solutions of both resulting DAE systems are the same. The symbolic equations before simplification are the same ( equations(structural_simplify(symb_mol)) == equations(structural_simplify(symbmtk)) returns true), but equations(structural_simplify(symb_mol)) == equations(structural_simplify(symbmtk)) is false. I cannot spot the difference between the ODESystem returned by MethodofLines.jl and the ODESystem which I defined via MTK.

xtalax commented 1 year ago

Oh great, it's one of those...

I will write a script to find the problem equations and determine the difference

xtalax commented 1 year ago

Hey, I realized that your problem can be represented with the new InterfaceBCs

https://docs.sciml.ai/MethodOfLines/stable/boundary_conditions/#Interfaces

Is this helpful? Of course i will try to fix the original bug too.

verflieger commented 1 year ago

At the moment I cannot see how this helps... For now, my dirty workaround is to use a very small but nonzero boundary condition or to add a small diffusion term. I stick to the first approach to keep the model complexity low for my real problem (parameter identification for a plug-flow reactor with several reactions and species)...

YingboMa commented 1 year ago

See https://github.com/SciML/ModelingToolkit.jl/issues/1948