search

SciML / ModelingToolkit.jl

An acausal modeling framework for automatically parallelized scientific machine learning (SciML) in Julia. A computer algebra system for integrated symbolics for physics-informed machine learning and automated transformations of differential equations
https://mtk.sciml.ai/dev/
Other
1.41k stars 204 forks source link

Solution of ODEProblem errors for u0 given as dictionary, but not for u0 given as vector #1502

Open rschiemann1 opened 2 years ago

rschiemann1 commented 2 years ago

I am constructing an ODEProblem with u0 defined as dictionary, see function build_problem in below code. When doing so, solving the constructed problem errors with no method matching oneunit(::Type{Any}). If instead I build u0 as a plain, boring vector of Floats, solution of system works.

using OrdinaryDiffEq, ModelingToolkit
using DifferentialEquations

function build_system(n)

    @parameters t
    vars = @variables (Tgas[1:n](t) = fill(100,n)), 
                      (Tsol[1:n](t) = fill(20, n)), 
                      (cR[1:n](t) = fill(3000, n)), 
                      (mflux(t) = 0.5)

    @parameters params[1:3]
    D = Differential(t)

    dp, Tin, h = params 

    dx = h / n

    ############################################################
    # Build equations ##########################################
    ############################################################

    eqns = Array{Equation}(undef, 3 * n + 1)
    count = 1

    # specification as function or intermediate variabel does not matter. I could also make this an 
    # actual @variable, then I could retrieve it later.
    rR = 5.5e6 .* exp.(.-5680 ./(273 .+ Tsol)) .* (1 .- exp.(.-cR./50))

    q = 1e4 * (Tgas - Tsol)

    # gas energy balance
    eqns[count] = 0 ~  mflux * 1000 * (Tin - Tgas[1]) - q[1] * dx
    count += 1
    for i in 2:n
        #eqns[count] = 0 ~ mflux * 1000 * (Tgas[i-1] - Tgas[i]) - q(Tsol[i], Tgas[i]) * dx
        eqns[count] = 0 ~ mflux * 1000 * (Tgas[i-1] - Tgas[i]) - q[i] * dx
        count += 1
    end

    # solids energy balance
    for i in 1:n
        #eqns[count] = D(Tsol[i]) ~ q(Tsol[i], Tgas[i]) / (1000 * 2000)
        eqns[count] = D(Tsol[i]) ~ q[i] / (1000 * 2000)
        count += 1
    end

    # gas momentum balance
    eqns[count] = 0 ~ - mflux[1] + 5e-3 * sqrt(dp)
    count += 1

    # reactant concentration equations
    for i in 1:n
        #eqns[count] = D(cR[i]) ~ -rR(Tsol[i], cR[i])
        eqns[count] = D(cR[i]) ~ -rR[i]

        count += 1
    end

    @named sys = ODESystem(eqns, t, vcat(vars...), params)

    return sys

end

function build_problem(sys; tspan = (0.0, 1.5), param_values = [100e2, 200, 0.5])

    u0map =   [ 
        (sys.Tgas .=> fill(100.0, N))...
        (sys.Tsol .=> fill(20.0, N))...
        (sys.cR .=> fill(3000.0, N))...
        (sys.mflux => 0.5)
    ]

    prob = ODEProblem(sys, u0map, tspan, param_values, jac=true, sparse=true)
    return prob

end

dp = 100e2
Tin = 200
h = 0.5
parameters = [dp, Tin, h]

N = 100
tspan = (0.0, 1500.0)

sys = build_system(N)
prob = build_problem(sys, tspan=tspan, param_values = parameters)

# error: no method matching oneunit(::Type{Any})
# error does not happen if a plain vector is used as u0 in ODEProblem construction.
sol = solve(prob, QBDF(), reltol = 1e-3)
ChrisRackauckas commented 2 years ago

It's probably https://github.com/SciML/ModelingToolkit.jl/pull/1490/files#diff-47c27891e951c8cd946b850dc2df31082624afdf57446c21cb6992f5f4b74aa2R475. I don't think that method was made safe for symbolic arrays @ValentinKaisermayer

rschiemann1 commented 2 years ago

Then how does everyone else specify their initial values for array variables? I don't think I am the first one to try this.

ChrisRackauckas commented 2 years ago

You change the variables to arrays of symbolics via collect(Tgas) etc.

ValentinKaisermayer commented 2 years ago 
function build_problem(sys; tspan = (0.0, 1.5), param_values = [100e2, 200, 0.5])
    u0map = [ 
        (collect(sys.Tgas) .=> fill(100.0, N))...
        (collect(sys.Tsol) .=> fill(20.0, N))...
        (collect(sys.cR) .=> fill(3000.0, N))...
        (sys.mflux => 0.5)
    ] 
    return ODEProblem(sys, u0map, tspan, param_values)
end

N = 3
tspan = (0.0, 1500.0)
sys = build_system(N)
prob = build_problem(sys, tspan=tspan, param_values = pars)

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 1500.0)
u0: 10-element Vector{Float64}:
  100.0
  100.0
  100.0
   20.0
   20.0
   20.0
 3000.0
 3000.0
 3000.0
    0.5

However, I get an error if I have ODEProblem(sys, u0map, tspan, param_values, jac=true, sparse=true)

ERROR: Differentiation of expressions involving arrays and array variables is not yet supported.
Stacktrace:
  [1] error(s::String)
    @ Base .\error.jl:33
  [2] occursin_info(x::Term{Real, Base.ImmutableDict{DataType, Any}}, expr::Sym{Vector{Real}, Base.ImmutableDict{DataType, Any}})
    @ Symbolics .julia\packages\Symbolics\1OrKJ\src\diff.jl:85
  [3] (::Symbolics.var"#133#135"{Term{Real, Base.ImmutableDict{DataType, Any}}})(a::Sym{Vector{Real}, Base.ImmutableDict{DataType, Any}})
    @ Symbolics .julia\packages\Symbolics\1OrKJ\src\diff.jl:75
  [4] iterate
    @ .\generator.jl:47 [inlined]
  [5] _collect(c::Vector{Any}, itr::Base.Generator{Vector{Any}, Symbolics.var"#133#135"{Term{Real, Base.ImmutableDict{DataType, Any}}}}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})
    @ Base .\array.jl:744
  [6] collect_similar(cont::Vector{Any}, itr::Base.Generator{Vector{Any}, Symbolics.var"#133#135"{Term{Real, Base.ImmutableDict{DataType, Any}}}})
    @ Base .\array.jl:653
  [7] map(f::Function, A::Vector{Any})
    @ Base .\abstractarray.jl:2849
  [8] occursin_info(x::Term{Real, Base.ImmutableDict{DataType, Any}}, expr::Term{Real, Base.ImmutableDict{DataType, Any}})
    @ Symbolics .julia\packages\Symbolics\1OrKJ\src\diff.jl:75
  [9] (::Symbolics.var"#133#135"{Term{Real, Base.ImmutableDict{DataType, Any}}})(a::Term{Real, Base.ImmutableDict{DataType, Any}})
    @ Symbolics .julia\packages\Symbolics\1OrKJ\src\diff.jl:75
 [10] iterate
    @ .\generator.jl:47 [inlined]
 [11] collect_to!(dest::Vector{Any}, itr::Base.Generator{Vector{Any}, Symbolics.var"#133#135"{Term{Real, Base.ImmutableDict{DataType, Any}}}}, offs::Int64, st::Int64)
    @ Base .\array.jl:782
 [12] collect_to!(dest::Vector{Bool}, itr::Base.Generator{Vector{Any}, Symbolics.var"#133#135"{Term{Real, Base.ImmutableDict{DataType, Any}}}}, offs::Int64, st::Int64)
    @ Base .\array.jl:790
 [13] collect_to_with_first!(dest::Vector{Bool}, v1::Bool, itr::Base.Generator{Vector{Any}, Symbolics.var"#133#135"{Term{Real, Base.ImmutableDict{DataType, Any}}}}, st::Int64)
    @ Base .\array.jl:760
 [14] _collect(c::Vector{Any}, itr::Base.Generator{Vector{Any}, Symbolics.var"#133#135"{Term{Real, Base.ImmutableDict{DataType, Any}}}}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})
    @ Base .\array.jl:754
 [15] collect_similar(cont::Vector{Any}, itr::Base.Generator{Vector{Any}, Symbolics.var"#133#135"{Term{Real, Base.ImmutableDict{DataType, Any}}}})
    @ Base .\array.jl:653
 [16] map(f::Function, A::Vector{Any})
    @ Base .\abstractarray.jl:2849
 [17] occursin_info(x::Term{Real, Base.ImmutableDict{DataType, Any}}, expr::SymbolicUtils.Mul{Real, Rational{Int64}, Dict{Any, Number}, Nothing})
    @ Symbolics .julia\packages\Symbolics\1OrKJ\src\diff.jl:75
 [18] (::Symbolics.var"#133#135"{Term{Real, Base.ImmutableDict{DataType, Any}}})(a::SymbolicUtils.Mul{Real, Rational{Int64}, Dict{Any, Number}, Nothing})
    @ Symbolics .julia\packages\Symbolics\1OrKJ\src\diff.jl:75
 [19] iterate
    @ .\generator.jl:47 [inlined]
 [20] _collect(c::Vector{SymbolicUtils.Mul{Real, T, Dict{Any, Number}, Nothing} where T<:Number}, itr::Base.Generator{Vector{SymbolicUtils.Mul{Real, T, Dict{Any, Number}, Nothing} where T<:Number}, Symbolics.var"#133#135"{Term{Real, Base.ImmutableDict{DataType, Any}}}}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})
    @ Base .\array.jl:744
 [21] collect_similar(cont::Vector{SymbolicUtils.Mul{Real, T, Dict{Any, Number}, Nothing} where T<:Number}, itr::Base.Generator{Vector{SymbolicUtils.Mul{Real, T, Dict{Any, Number}, Nothing} where T<:Number}, Symbolics.var"#133#135"{Term{Real, Base.ImmutableDict{DataType, Any}}}})
    @ Base .\array.jl:653
 [22] map(f::Function, A::Vector{SymbolicUtils.Mul{Real, T, Dict{Any, Number}, Nothing} where T<:Number})
    @ Base .\abstractarray.jl:2849
 [23] occursin_info(x::Term{Real, Base.ImmutableDict{DataType, Any}}, expr::SymbolicUtils.Add{Real, Int64, Dict{Any, Number}, Nothing})
    @ Symbolics .julia\packages\Symbolics\1OrKJ\src\diff.jl:75
 [24] expand_derivatives(O::Term{Real, Nothing}, simplify::Bool; occurances::Nothing)
    @ Symbolics .julia\packages\Symbolics\1OrKJ\src\diff.jl:137
ChrisRackauckas commented 2 years ago

This isolates it:

using OrdinaryDiffEq, ModelingToolkit
using DifferentialEquations

function build_system(n)

    @parameters t
    vars = collect.(@variables (Tgas[1:n](t) = fill(100,n)), 
                      (Tsol[1:n](t) = fill(20, n)), 
                      (cR[1:n](t) = fill(3000, n)))
    push!(vars, (@variables mflux(t) = 0.5))
    params = collect(@parameters(params[1:3])[1])
    D = Differential(t)

    dp, Tin, h = params 

    dx = h / n

    ############################################################
    # Build equations ##########################################
    ############################################################

    eqns = Array{Equation}(undef, 3 * n + 1)
    count = 1

    # specification as function or intermediate variabel does not matter. I could also make this an 
    # actual @variable, then I could retrieve it later.
    rR = 5.5e6 .* exp.(.-5680 ./(273 .+ Tsol)) .* (1 .- exp.(.-cR./50))

    q = 1e4 * (Tgas - Tsol)

    # gas energy balance
    eqns[count] = 0 ~  mflux * 1000 * (Tin - Tgas[1]) - q[1] * dx
    count += 1
    for i in 2:n
        #eqns[count] = 0 ~ mflux * 1000 * (Tgas[i-1] - Tgas[i]) - q(Tsol[i], Tgas[i]) * dx
        eqns[count] = 0 ~ mflux * 1000 * (Tgas[i-1] - Tgas[i]) - q[i] * dx
        count += 1
    end

    # solids energy balance
    for i in 1:n
        #eqns[count] = D(Tsol[i]) ~ q(Tsol[i], Tgas[i]) / (1000 * 2000)
        eqns[count] = D(Tsol[i]) ~ q[i] / (1000 * 2000)
        count += 1
    end

    # gas momentum balance
    eqns[count] = 0 ~ - mflux[1] + 5e-3 * sqrt(dp)
    count += 1

    # reactant concentration equations
    for i in 1:n
        #eqns[count] = D(cR[i]) ~ -rR(Tsol[i], cR[i])
        eqns[count] = D(cR[i]) ~ -rR[i]

        count += 1
    end

    @named sys = ODESystem(eqns, t, vcat(vars...), params)

    return sys

end

function build_problem(sys; tspan = (0.0, 1.5), param_values = [100e2, 200, 0.5])

    u0map =   [ 
        (sys.Tgas .=> fill(100.0, N))...
        (sys.Tsol .=> fill(20.0, N))...
        (sys.cR .=> fill(3000.0, N))...
        (sys.mflux => 0.5)
    ]

    prob = ODEProblem(sys, u0map, tspan, param_values, jac=true, sparse=true)
    return prob

end

dp = 100e2
Tin = 200
h = 0.5
parameters = [dp, Tin, h]

N = 100
tspan = (0.0, 1500.0)

sys = build_system(N)
prob = build_problem(sys, tspan=tspan, param_values=parameters)

Which gave the error:

ERROR: Differentiation of expressions involving arrays and array variables is not yet supported.
Stacktrace:
  [1] error(s::String)
    @ Base .\error.jl:33
  [2] occursin_info(x::Term{Real, Base.ImmutableDict{DataType, Any}}, expr::Sym{Vector{Real}, Base.ImmutableDict{DataType, Any}})
    @ Symbolics C:\Users\accou\.julia\packages\Symbolics\1OrKJ\src\diff.jl:85
  [3] (::Symbolics.var"#884#886"{Term{Real, Base.ImmutableDict{DataType, Any}}})(a::Sym{Vector{Real}, Base.ImmutableDict{DataType, Any}})
    @ Symbolics c:\Users\accou\.julia\packages\Symbolics\1OrKJ\src\diff.jl:76
  [4] iterate
    @ .\generator.jl:47 [inlined]
  [5] _collect(c::Vector{Any}, itr::Base.Generator{Vector{Any}, Symbolics.var"#884#886"{Term{Real, Base.ImmutableDict{DataType, Any}}}}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})
    @ Base .\array.jl:744
  [6] collect_similar(cont::Vector{Any}, itr::Base.Generator{Vector{Any}, Symbolics.var"#884#886"{Term{Real, Base.ImmutableDict{DataType, Any}}}})
    @ Base .\array.jl:653
  [7] map(f::Function, A::Vector{Any})
    @ Base .\abstractarray.jl:2849
  [8] occursin_info(x::Term{Real, Base.ImmutableDict{DataType, Any}}, expr::Term{Real, Base.ImmutableDict{DataType, Any}})
    @ Symbolics c:\Users\accou\.julia\packages\Symbolics\1OrKJ\src\diff.jl:76
  [9] (::Symbolics.var"#884#886"{Term{Real, Base.ImmutableDict{DataType, Any}}})(a::Term{Real, Base.ImmutableDict{DataType, Any}})
    @ Symbolics c:\Users\accou\.julia\packages\Symbolics\1OrKJ\src\diff.jl:76
 [10] iterate
    @ .\generator.jl:47 [inlined]
 [11] collect_to!(dest::Vector{Any}, itr::Base.Generator{Vector{Any}, Symbolics.var"#884#886"{Term{Real, Base.ImmutableDict{DataType, Any}}}}, offs::Int64, st::Int64)     
    @ Base .\array.jl:782
 [12] collect_to!(dest::Vector{Bool}, itr::Base.Generator{Vector{Any}, Symbolics.var"#884#886"{Term{Real, Base.ImmutableDict{DataType, Any}}}}, offs::Int64, st::Int64)    
    @ Base .\array.jl:790
 [13] collect_to_with_first!(dest::Vector{Bool}, v1::Bool, itr::Base.Generator{Vector{Any}, Symbolics.var"#884#886"{Term{Real, Base.ImmutableDict{DataType, Any}}}}, st::Int64)
    @ Base .\array.jl:760
 [14] _collect(c::Vector{Any}, itr::Base.Generator{Vector{Any}, Symbolics.var"#884#886"{Term{Real, Base.ImmutableDict{DataType, Any}}}}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})
    @ Base .\array.jl:754
 [15] collect_similar(cont::Vector{Any}, itr::Base.Generator{Vector{Any}, Symbolics.var"#884#886"{Term{Real, Base.ImmutableDict{DataType, Any}}}})
    @ Base .\array.jl:653
 [16] map(f::Function, A::Vector{Any})
    @ Base .\abstractarray.jl:2849
 [17] occursin_info(x::Term{Real, Base.ImmutableDict{DataType, Any}}, expr::SymbolicUtils.Mul{Real, Rational{Int64}, Dict{Any, Number}, Nothing})
    @ Symbolics c:\Users\accou\.julia\packages\Symbolics\1OrKJ\src\diff.jl:76
 [18] (::Symbolics.var"#884#886"{Term{Real, Base.ImmutableDict{DataType, Any}}})(a::SymbolicUtils.Mul{Real, Rational{Int64}, Dict{Any, Number}, Nothing})
    @ Symbolics c:\Users\accou\.julia\packages\Symbolics\1OrKJ\src\diff.jl:76
 [19] iterate
    @ .\generator.jl:47 [inlined]
 [20] _collect(c::Vector{SymbolicUtils.Mul{Real, T, Dict{Any, Number}, Nothing} where T<:Number}, itr::Base.Generator{Vector{SymbolicUtils.Mul{Real, T, Dict{Any, Number}, Nothing} where T<:Number}, Symbolics.var"#884#886"{Term{Real, Base.ImmutableDict{DataType, Any}}}}, #unused#::Base.EltypeUnknown, isz::Base.HasShape{1})
    @ Base .\array.jl:744
 [21] collect_similar(cont::Vector{SymbolicUtils.Mul{Real, T, Dict{Any, Number}, Nothing} where T<:Number}, itr::Base.Generator{Vector{SymbolicUtils.Mul{Real, T, Dict{Any, Number}, Nothing} where T<:Number}, Symbolics.var"#884#886"{Term{Real, Base.ImmutableDict{DataType, Any}}}})
    @ Base .\array.jl:653
 [22] map(f::Function, A::Vector{SymbolicUtils.Mul{Real, T, Dict{Any, Number}, Nothing} where T<:Number})
    @ Base .\abstractarray.jl:2849
 [23] occursin_info(x::Term{Real, Base.ImmutableDict{DataType, Any}}, expr::SymbolicUtils.Add{Real, Int64, Dict{Any, Number}, Nothing})
    @ Symbolics c:\Users\accou\.julia\packages\Symbolics\1OrKJ\src\diff.jl:76
 [24] expand_derivatives(O::Term{Real, Nothing}, simplify::Bool; occurances::Nothing)
    @ Symbolics C:\Users\accou\.julia\packages\Symbolics\1OrKJ\src\diff.jl:137
 [25] expand_derivatives(O::Term{Real, Nothing}, simplify::Bool)
    @ Symbolics C:\Users\accou\.julia\packages\Symbolics\1OrKJ\src\diff.jl:132
 [26] sparsejacobian(ops::Vector{SymbolicUtils.Symbolic{Real}}, vars::Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}; simplify::Bool)
    @ Symbolics C:\Users\accou\.julia\packages\Symbolics\1OrKJ\src\diff.jl:435
 [27] calculate_jacobian(sys::ODESystem; sparse::Bool, simplify::Bool, dvs::Vector{Term{Real, Base.ImmutableDict{DataType, Any}}})
    @ ModelingToolkit c:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:35
 [28] generate_jacobian(sys::ODESystem, dvs::Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}, ps::Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}; simplify::Bool, sparse::Bool, kwargs::Base.Pairs{Symbol, Any, NTuple{7, Symbol}, NamedTuple{(:expression, :expression_module, :checkbounds, :ddvs, :linenumbers, :parallel, :has_difference), Tuple{DataType, Module, Bool, Nothing, Bool, Symbolics.SerialForm, Bool}}})
    @ ModelingToolkit C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:70
 [29] (ODEFunction{true})(sys::ODESystem, dvs::Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}, ps::Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}, u0::Vector{Any}; version::Nothing, tgrad::Bool, jac::Bool, eval_expression::Bool, sparse::Bool, simplify::Bool, eval_module::Module, steady_state::Bool, checkbounds::Bool, sparsity::Bool, kwargs::Base.Pairs{Symbol, Any, NTuple{4, Symbol}, NamedTuple{(:ddvs, :linenumbers, :parallel, :has_difference), Tuple{Nothing, Bool, Symbolics.SerialForm, Bool}}})
    @ ModelingToolkit C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:357
 [30] process_DEProblem(constructor::Type, sys::ODESystem, u0map::Vector{Any}, parammap::Vector{Float64}; implicit_dae::Bool, du0map::Nothing, version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, kwargs::Base.Pairs{Symbol, Bool, Tuple{Symbol}, NamedTuple{(:has_difference,), Tuple{Bool}}})
    @ ModelingToolkit C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:593
 [31] (ODEProblem{true})(sys::ODESystem, u0map::Vector{Any}, tspan::Tuple{Float64, Float64}, parammap::Vector{Float64}; callback::Nothing, kwargs::Base.Pairs{Symbol, Bool, Tuple{Symbol, Symbol}, NamedTuple{(:jac, :sparse), Tuple{Bool, Bool}}})
    @ ModelingToolkit C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:671
 [32] #ODEProblem#334
    @ C:\Users\accou\.julia\dev\ModelingToolkit\src\systems\diffeqs\abstractodesystem.jl:649 [inlined]
 [33] build_problem(sys::ODESystem; tspan::Tuple{Float64, Float64}, param_values::Vector{Float64})
    @ Main c:\Users\accou\OneDrive\Computer\Desktop\test.jl:75
 [34] top-level scope
    @ c:\Users\accou\OneDrive\Computer\Desktop\test.jl:91

This prompted the question, what expressions were causing the error? So I went to the spot in the stacktrace and just made it print:

function occursin_info(x, expr)
    @show x, expr

and then reran:

#=
(x, expr) = (Tgas[1](t), 1000(params[2] - Tgas[1](t))*mflux(t) - (1//100)*(10000.0Tgas[1](t) - 10000.0Tsol[1](t))*params[3])
(x, expr) = (Tgas[1](t), (-1//100)*(10000.0Tgas[1](t) - 10000.0Tsol[1](t))*params[3])
(x, expr) = (Tgas[1](t), -1//100)
(x, expr) = (Tgas[1](t), 10000.0Tgas[1](t) - 10000.0Tsol[1](t))
(x, expr) = (Tgas[1](t), 10000.0Tgas[1](t))
(x, expr) = (Tgas[1](t), 10000.0)
(x, expr) = (Tgas[1](t), Tgas[1](t))
(x, expr) = (Tgas[1](t), -10000.0Tsol[1](t))
(x, expr) = (Tgas[1](t), -10000.0)
(x, expr) = (Tgas[1](t), Tsol[1](t))
(x, expr) = (Tgas[1](t), params[3])
=#

This isolates the error to:

Symbolics.derivative(ModelingToolkit.parameters(sys)[3], sys.Tgas[1])

which shows there's very specific cases which trigger it (not all array scalar differentiation fails, many work before it hits the error!). This then constructs the MWE:

using Symbolics

n = 10
vars = collect(@variables(Tgas[1:n](t))[1])
ps = collect(@parameters(ps[1:3])[1])
Symbolics.derivative(ps[3], vars[1])

@shashi do you know why this would be a bad case? It's just the derivative of a constant by a variable (so... zero).

rschiemann1 commented 2 years ago

Thank you two for looking into this. Just for the record, this is (in my view) the same issue as in #1501.

ChrisRackauckas commented 2 years ago

No that's unrelated.

rschiemann1 commented 2 years ago

Oh. Of course I have almost no insight into the inner workings here, but @ValentinKaisermayer reported a running code with ODEProblem(sys, u0map, tspan, param_values) and reported the same error as in #1501 with ODEProblem(sys, u0map, tspan, param_values, jac=true, sparse=true). So from an outside point of view, the exact same error messages surface under the same conditions. A possible difference between the two cases may be that in #1501, the system was structurally simplified.

rschiemann1 commented 2 years ago 
function build_problem(sys; tspan = (0.0, 1.5), param_values = [100e2, 200, 0.5])
    u0map = [ 
        (collect(sys.Tgas) .=> fill(100.0, N))...
        (collect(sys.Tsol) .=> fill(20.0, N))...
        (collect(sys.cR) .=> fill(3000.0, N))...
        (sys.mflux => 0.5)
    ] 
    return ODEProblem(sys, u0map, tspan, param_values)
end

N = 3
tspan = (0.0, 1500.0)
sys = build_system(N)
prob = build_problem(sys, tspan=tspan, param_values = pars)

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 1500.0)
u0: 10-element Vector{Float64}:
  100.0
  100.0
  100.0
   20.0
   20.0
   20.0
 3000.0
 3000.0
 3000.0
    0.5

Is that a solution to my original issue? The code in my original post errors on the solve, not on the ODEProblem construction. It also seems to me like all following disussions about ODEProblem building. I am getting a bit lost, are we still discussing my original issue?

ChrisRackauckas commented 2 years ago

Yes, the workaround here should be collect, but it's not working because of the recent derivative check. So fixing that at least makes the workaround work, and then the full symbolic array stuff is step two.

rschiemann1 commented 2 years ago

Okay, thank you for the clarification. The workaround with collect errors on the solve on my side as well - both with and without jac and sparse in ODEProblem construction. Not sure if that is already on the table.

edit: No, I was wrong. I see the same behavior as @ValentinKaisermayer: With the collect workaround, everything is working without jac and sparse. With jac and sparse, I see the same issue as he does.

ValentinKaisermayer commented 2 years ago 
using ModelingToolkit, OrdinaryDiffEq

function build_system(n)

    @parameters t
    vars = @variables (Tgas[1:n](t) = fill(100,n)), 
                      (Tsol[1:n](t) = fill(20, n)), 
                      (cR[1:n](t) = fill(3000, n)), 
                      (mflux(t) = 0.5)

    @parameters params[1:3]
    D = Differential(t)

    dp, Tin, h = params 

    dx = h / n

    ############################################################
    # Build equations ##########################################
    ############################################################

    eqns = Array{Equation}(undef, 3 * n + 1)
    count = 1

    # specification as function or intermediate variabel does not matter. I could also make this an 
    # actual @variable, then I could retrieve it later.
    rR = 5.5e6 .* exp.(.-5680 ./(273 .+ Tsol)) .* (1 .- exp.(.-cR./50))

    q = 1e4 * (Tgas - Tsol)

    # gas energy balance
    eqns[count] = 0 ~  mflux * 1000 * (Tin - Tgas[1]) - q[1] * dx
    count += 1
    for i in 2:n
        #eqns[count] = 0 ~ mflux * 1000 * (Tgas[i-1] - Tgas[i]) - q(Tsol[i], Tgas[i]) * dx
        eqns[count] = 0 ~ mflux * 1000 * (Tgas[i-1] - Tgas[i]) - q[i] * dx
        count += 1
    end

    # solids energy balance
    for i in 1:n
        #eqns[count] = D(Tsol[i]) ~ q(Tsol[i], Tgas[i]) / (1000 * 2000)
        eqns[count] = D(Tsol[i]) ~ q[i] / (1000 * 2000)
        count += 1
    end

    # gas momentum balance
    eqns[count] = 0 ~ - mflux[1] + 5e-3 * sqrt(dp)
    count += 1

    # reactant concentration equations
    for i in 1:n
        #eqns[count] = D(cR[i]) ~ -rR(Tsol[i], cR[i])
        eqns[count] = D(cR[i]) ~ -rR[i]

        count += 1
    end

    @named sys = ODESystem(eqns, t, vcat(vars...), params)

    return sys

end

function build_problem(sys; tspan = (0.0, 1.5), param_values = [100e2, 200, 0.5])

    u0map = [ 
        (collect(sys.Tgas) .=> fill(100.0, N))...
        (collect(sys.Tsol) .=> fill(20.0, N))...
        (collect(sys.cR) .=> fill(3000.0, N))...
        (sys.mflux => 0.5)
    ]

    prob = ODEProblem(sys, u0map, tspan, param_values, 
        # jac=true, sparse=true,
    )
    return prob

end

dp = 100e2
Tin = 200
h = 0.5
pars = [dp, Tin, h]

N = 100
tspan = (0.0, 1500.0)
sys = build_system(N)
prob = build_problem(sys, tspan=tspan, param_values=pars)

sol = solve(prob, QBDF(), reltol = 1e-3)
retcode: Success
Interpolation: 3rd order Hermite
t: 644-element Vector{Float64}:
    0.0
    1.0e-6
    2.0e-6
    3.0e-6
    1.3e-5
    2.2999999999999997e-5
    ⋮
 1492.9298837024442
 1495.1034815980252
 1497.2770794936062
 1499.4506773891871
 1500.0
u: 644-element Vector{Vector{Float64}}:
 [183.63636363636363, 168.7603305785124, 155.236664162284, 142.94242196571273, 131.76583815064794, 121.60530740967994, 112.3684612815272, 103.971328437752, 96.33757130704728, 89.39779209731572  …  3000.0, 3000.0, 3000.0, 3000.0, 3000.0, 3000.0, 3000.0, 3000.0, 3000.0, 0.5]
 [183.6363637107438, 168.76033071374906, 155.23666434669764, 142.9424221892444, 131.7658384046612, 121.60530768678532, 112.36846157542686, 103.97132874310229, 96.33757161933735, 89.39779241276023  …  2999.9999999790452, 2999.9999999790452, 2999.9999999790452, 2999.9999999790452, 2999.9999999790452, 2999.9999999790452, 2999.9999999790452, 2999.9999999790452, 2999.9999999790452, 0.5]
 [183.63636378512396, 168.76033084898572, 155.23666453111127, 142.94242241277607, 131.76583865867445, 121.60530796389071, 112.36846186932651, 103.97132904845259, 96.33757193162741, 89.39779272820473  …  2999.999999958091, 2999.999999958091, 2999.999999958091, 2999.999999958091, 2999.999999958091, 2999.999999958091, 2999.999999958091, 2999.999999958091, 2999.999999958091, 0.5]
 [183.63636385950414, 168.76033098422238, 155.23666471552488, 142.94242263630775, 131.76583891268774, 121.60530824099608, 112.36846216322616, 103.97132935380286, 96.33757224391748, 89.39779304364924  …  2999.9999999371366, 2999.9999999371366, 2999.9999999371366, 2999.9999999371366, 2999.9999999371366, 2999.9999999371366, 2999.9999999371366, 2999.9999999371366, 2999.9999999371366, 0.5]
 [183.63636460330576, 168.760332336589, 155.23666655966116, 142.94242487162444, 131.76584145282035, 121.60531101204985, 112.36846510222259, 103.97133240730565, 96.33757536681806, 89.39779619809428  …  2999.9999997275927, 2999.9999997275927, 2999.9999997275927, 2999.9999997275927, 2999.9999997275927, 2999.9999997275927, 2999.9999997275927, 2999.9999997275927, 2999.9999997275927, 0.5]
 [183.63636534710736, 168.7603336889555, 155.23666840379735, 142.94242710694104, 131.76584399295288, 121.60531378310353, 112.36846804121892, 103.97133546080835, 96.33757848971857, 89.39779935253924  …  2999.9999995180488, 2999.9999995180488, 2999.9999995180488, 2999.9999995180488, 2999.9999995180488, 2999.9999995180488, 2999.9999995180488, 2999.9999995180488, 2999.9999995180488, 0.5]
 ⋮
 [199.98134292690412, 199.95294196513478, 199.91320479204953, 199.8605062751289, 199.7931984613385, 199.7096203863243, 199.60810764602462, 199.48700167916658, 199.3446587148666, 199.17945834514524  …  2874.3654884312796, 2880.46887651642, 2886.1323601254185, 2891.390084042567, 2896.2734568590317, 2900.8113674401307, 2905.030385455698, 2908.954946939496, 2912.6075258156047, 0.5]
 [199.98152635100192, 199.9533880161691, 199.91400116277046, 199.8617487027554, 199.79499022420163, 199.7120717496069, 199.6113352635379, 199.491127957983, 199.3498111494614, 199.18576882749298  …  2873.450285045084, 2879.6120671316075, 2885.3297111103516, 2890.637676790367, 2895.567665104997, 2900.148836019693, 2904.4080104419854, 2908.369857097457, 2912.0570653084487, 0.5]
 [199.9817079717944, 199.95382984505076, 199.91479024916475, 199.8629801196625, 199.79676657168056, 199.71450260723552, 199.6145365987342, 199.49522149695574, 199.35492371739795, 199.1920316643912  …  2872.527516747898, 2878.7482046985174, 2884.520485690924, 2889.8791368358347, 2894.8561535673743, 2899.4809689932576, 2903.780657186879, 2907.7801213829134, 2911.502268011835, 0.5]
 [199.98188780701048, 199.95426749168544, 199.91557211763336, 199.86420062287993, 199.7985276353314, 199.7169131289097, 199.61771186275416, 199.49928255161163, 199.35999672115398, 199.19824720743117  …  2871.5971231446715, 2877.877232644986, 2883.704630900183, 2889.114414607855, 2894.1388758700928, 2898.8077229877517, 2903.148285137084, 2907.1857018879437, 2910.9430984975597, 0.5]
 [199.98193297531046, 199.95437743861544, 199.91576858062655, 199.86450735844332, 199.79897029634742, 199.71751912899782, 199.61851022940806, 199.50030376980197, 199.36127256904774, 199.1998105858311  …  2871.3607751012896, 2877.655984003059, 2883.497388223221, 2888.920165589697, 2893.956683524072, 2898.6367202010847, 2902.987669575532, 2907.034731235199, 2910.8010860134027, 0.5]

(test) pkg> st
  [961ee093] ModelingToolkit v8.6.0
  [1dea7af3] OrdinaryDiffEq v6.8.0
rschiemann1 commented 2 years ago

With yesterday's release of Symbolics (closing this issue), constructing the ODEProblem with jac=true, sparse = true works, which is one step closer.

Unfortunately, solving the constructed ODEProblem fails with a LinearAlgebra.SingularException error. But I will create a separate issue (edit: see #1517) for this.

Anyway, I think this issue should remain open because hunting down the actual bug of this issue (as isolated in Chris' longer post above) is still ongoing.

ChrisRackauckas commented 2 years ago

@YingboMa could you check would could be missing here?

YingboMa commented 2 years ago

Symbolics seems to give the wrong Jacobian. With https://github.com/SciML/ModelingToolkit.jl/pull/1521 I get

julia> prob_simplified_jacSparse = ODEProblem(sys_simplified, Pair[], tspan, params, jac=true, sparse=true);

julia> sol_simplified_jacSparse = solve(prob_simplified_jacSparse, QBDF());
ERROR: LinearAlgebra.SingularException(0)
Stacktrace:
...

julia> prob_simplified_jac = ODEProblem(sys_simplified, Pair[], tspan, params, jac=true, sparse=false);

julia> sol_simplified_jac = solve(prob_simplified_jac, QBDF());
ERROR: LinearAlgebra.SingularException(12)
Stacktrace:
  [1] checknonsingular
...

@shashi

AayushSabharwal commented 6 months ago

This works with #2605 

julia> function build_system(n)
           vars = @variables (Tgas(t)[1:n] = fill(100,n)),
                             (Tsol(t)[1:n] = fill(20, n)),
                             (cR(t)[1:n] = fill(3000, n)),
                             (mflux(t) = 0.5)

           @parameters params[1:3]
           dp, Tin, h = params
           dx = h / n
           eqns = Equation[]
           count = 1
           # specification as function or intermediate variabel does not matter. I could also make this an
           # actual @variable, then I could retrieve it later.
           rR = (5.5e6 .* exp.(.-5680 ./(273 .+ Tsol)) .* (1 .- exp.(.-cR./50)))
           q = 1e4 * (Tgas - Tsol)
           # gas energy balance
           push!(eqns, 0 ~  mflux * 1000 * (Tin - Tgas[1]) - q[1] * dx)
           count += 1
           for i in 2:n
               #eqns[count] = 0 ~ mflux * 1000 * (Tgas[i-1] - Tgas[i]) - q(Tsol[i], Tgas[i]) * dx
               push!(eqns, 0 ~ mflux * 1000 * (Tgas[i-1] - Tgas[i]) - q[i] * dx)
               count += 1
           end
           # solids energy balance
           for i in 1:n
               #eqns[count] = D(Tsol[i]) ~ q(Tsol[i], Tgas[i]) / (1000 * 2000)
               push!(eqns, D(Tsol[i]) ~ q[i] / (1000 * 2000))
               count += 1
           end
           # gas momentum balance
           push!(eqns, 0 ~ - mflux[1] + 5e-3 * sqrt(dp))
           count += 1
           # reactant concentration equations
           for i in 1:n
               #eqns[count] = D(cR[i]) ~ -rR(Tsol[i], cR[i])
               push!(eqns, D(cR[i]) ~ -rR[i])
               count += 1
           end
           @mtkbuild sys = ODESystem(eqns, t)
       end
julia> sys = build_system(10)
Model sys with 30 equations
Unknowns (30):
  (Tsol(t))[1] [defaults to 20]
  (Tsol(t))[2] [defaults to 20]
  (Tsol(t))[3] [defaults to 20]
  (Tsol(t))[4] [defaults to 20]
  (Tsol(t))[5] [defaults to 20]
  (Tsol(t))[6] [defaults to 20]
  (Tsol(t))[7] [defaults to 20]
  (Tsol(t))[8] [defaults to 20]
  (Tsol(t))[9] [defaults to 20]
  (Tsol(t))[10] [defaults to 20]
⋮
Parameters (1):
  params
Incidence matrix:30×50 SparseArrays.SparseMatrixCSC{Num, Int64} with 89 stored entries:
⎡⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⎤
⎢⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⎥
⎢⢄⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠑⎥
⎢⠀⠑⢄⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠑⢄⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⎥
⎢⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠳⣄⠀⠀⠀⎥
⎢⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⎥
⎣⠀⠀⠀⠀⠑⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠓⎦
julia> prob = ODEProblem(sys, [sys.Tgas => 100ones(10), sys.Tsol => 20ones(10), sys.cR => 3000ones(10), sys.mflux => 0.5], (0.0, 1.5), [sys.params => [100e2, 200, 0.5]])
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 1.5)
u0: 30-element Vector{Float64}:
   20.0
   20.0
   20.0
   20.0
   20.0
   20.0
   20.0
   20.0
   20.0
   20.0
 3000.0
 3000.0
 3000.0
 3000.0
 3000.0
 3000.0
 3000.0
 3000.0
 3000.0
 3000.0
  100.0
  100.0
  100.0
  100.0
  100.0
  100.0
  100.0
  100.0
  100.0
  100.0
AayushSabharwal commented 6 months ago

The fix in #2605 is actually incorrect and breaks other things. Upon further investigation, this looks like a Symbolics codegen bug. I have not been successful in reducing to an MWE, but the above example leads to the following error when the ODEProblem is constructed:

ERROR: type NameState has no field symbolify
Stacktrace:
  [1] getproperty
    @ ./Base.jl:37 [inlined]
  [2] toexpr(x::Symbolics.ArrayOp{AbstractVector{Real}}, st::SymbolicUtils.Code.NameState)
    @ Symbolics ~/Julia/SciML/Symbolics.jl/src/arrays.jl:993
  [3] (::SymbolicUtils.Code.var"#4#5"{SymbolicUtils.Code.NameState})(x::Symbolics.ArrayOp{AbstractVector{Real}})
    @ SymbolicUtils.Code ~/.julia/packages/SymbolicUtils/c0xQb/src/code.jl:187
  [4] iterate
    @ ./generator.jl:47 [inlined]
  [5] _collect(c::Vector{…}, itr::Base.Generator{…}, ::Base.EltypeUnknown, isz::Base.HasShape{…})
    @ Base ./array.jl:854
  [6] collect_similar(cont::Vector{Any}, itr::Base.Generator{Vector{…}, SymbolicUtils.Code.var"#4#5"{…}})
    @ Base ./array.jl:763
  [7] map(f::Function, A::Vector{Any})
    @ Base ./abstractarray.jl:3285
  [8] toexpr(O::SymbolicUtils.BasicSymbolic{Real}, st::SymbolicUtils.Code.NameState)
    @ SymbolicUtils.Code ~/.julia/packages/SymbolicUtils/c0xQb/src/code.jl:187
  [9] (::Base.Fix2{typeof(toexpr), SymbolicUtils.Code.NameState})(y::SymbolicUtils.BasicSymbolic{Real})

This can be "fixed" by Symbolics.scalarize-ing rR and q in the system constructor function.

AayushSabharwal commented 6 months ago

I also noticed that if we don't scalarize the broadcasts, the incidence matrix is different.

With scalarization of rR and q:

Incidence matrix:30×50 SparseArrays.SparseMatrixCSC{Num, Int64} with 89 stored entries:
⎡⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⎤
⎢⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⎥
⎢⢄⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠑⎥
⎢⠀⠑⢄⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠑⢄⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⎥
⎢⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠳⣄⠀⠀⠀⎥
⎢⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⣄⠀⎥
⎣⠀⠀⠀⠀⠑⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠓⎦

Without scalarization of rR and q:

Incidence matrix:30×50 SparseArrays.SparseMatrixCSC{Num, Int64} with 620 stored entries:
⎡⣿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⠀⣿⣿⣿⣿⣿⎤
⎢⣿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⣿⣿⣿⣿⣿⎥
⎢⣿⣿⣿⣿⣿⣤⣤⣤⣤⣤⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠛⠛⠛⠛⠛⎥
⎢⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⎥
⎢⣿⣿⣿⣿⣿⣿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⎥
⎢⣿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣿⣿⣿⣿⣿⎥
⎢⣿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣿⣿⣿⣿⣿⎥
⎣⠛⠛⠛⠛⠛⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠛⠛⠛⠛⠛⎦

We need a method in symbolics that does this scalarization without scalarizing symbolic arrays (i.e. turn broadcast(+, x, 1)[1] into x[1] + 1 but don't turn my_registered_function(x) into my_registered_function([x[1], x[2], x[3]])). Alternatively, linear_expansion needs to work for such broadcasted expressions.