ModelingToolkit.jl Solution of ODEProblem errors for u0 given as dictionary, but not for u0 given as vector

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)

Mar 25 '22 14:03 rschiemann1

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

Mar 25 '22 18:03 ChrisRackauckas

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

Apr 05 '22 15:04 rschiemann1

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

Apr 06 '22 00:04 ChrisRackauckas

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

Apr 06 '22 13:04 ValentinKaisermayer

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).

Apr 06 '22 16:04 ChrisRackauckas

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

Apr 06 '22 17:04 rschiemann1

No that's unrelated.

Apr 06 '22 17:04 ChrisRackauckas

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.

Apr 06 '22 17:04 rschiemann1

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?

Apr 06 '22 18:04 rschiemann1

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.

Apr 06 '22 18:04 ChrisRackauckas

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.

Apr 06 '22 19:04 rschiemann1

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

Apr 07 '22 14:04 ValentinKaisermayer

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.

Apr 08 '22 07:04 rschiemann1

@YingboMa could you check would could be missing here?

Apr 08 '22 14:04 ChrisRackauckas

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

Apr 09 '22 00:04 YingboMa

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

Apr 04 '24 13:04 AayushSabharwal

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.

Apr 05 '24 05:04 AayushSabharwal

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.

Apr 05 '24 06:04 AayushSabharwal

ModelingToolkit.jl ModelingToolkit.jl copied to clipboard

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

ModelingToolkit.jl
ModelingToolkit.jl copied to clipboard