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 203 forks source link

`promote_rule_constant` error on `ODEFunction(system; sparse = true)` #1085

Open cadojo opened 3 years ago

cadojo commented 3 years ago

I'm modeling the Circular Restricted Three-body Problem within AstrodynamicalModels.jl. I can create an ODEFunction from the ODESystem for a dense output, but only 14 elements are non-zero in the Jacobian, so I'd like to sparsify.

I think this is a bug? Sorry, I'm not super familiar with sparse math and could be missing something in the implementation / modeling syntax. Thanks for any guidance anyone can provide. 


"""
A `ModelingToolkit.ODESystem` for the Circular Restricted Three-body Problem,
with the local linearization included in the state vector and dynamics.
"""
CR3BPWithSTM = let 

    @parameters t μ 
    @variables x(t) y(t) z(t) ẋ(t) ẏ(t) ż(t) Φ[1:6,1:6](t)
    δ = Differential(t)
    r = @SVector [x,y,z]
    v = @SVector [ẋ,ẏ,ż]

    # Potential energy of spacecraft (this code was symbolically generated with Symbolics.jl, then copied and pasted)
    U = μ*(sqrt((μ + x - 1)^2 + y^2 + z^2)^-1) + (1//2) * (x^2) + (1//2) * (y^2) + (sqrt(y^2 + z^2 + (μ + x)^2)^-1)*(1 - μ)

    # Hessian of potential energy
    H = Symbolics.hessian(U, r)

    eqs = let
        LHS = [δ(Φ[i,j]) for j in 1:size(Φ,2) for i in 1:size(Φ,1)]
        RHS = vec(vcat(
            Matrix{Float64}(hcat(zeros(3,3), I(3))),
            hcat(H, Float64[0 2 0; -2 0 0; 0 0 0])
        ) * Φ)

        @assert length(LHS) == length(RHS) == 36 "If this assertion fails, please file an issue at https://github.com/cadojo/AstrodynamicalModels.jl!"
        [LHS[i] ~ RHS[i] for i in 1:length(LHS)]
    end

    @named CR3BPWithSTM = ODESystem(vcat(equations(CR3BP)..., eqs...), t, vcat(r,v,Φ...), [μ])
end

"""
A `DifferentialEquations`-compatible `ODEFunction` for R2BP dynamics.
Note that this function has several methods, including an in-place 
method! Function signatures follow `ModelingToolkit` and `DifferentialEquations`
conventions.
"""
CR3BPWithSTMVectorField = ODEFunction(CR3BPWithSTM; jac = true, tgrad = false, sparse = true)
ERROR: MethodError: no method matching promote_rule_constant(::Type{Int64}, ::Type{MultivariatePolynomials.AbstractPolynomialLike})
Closest candidates are:
  promote_rule_constant(::Type{T}, ::Type{MultivariatePolynomials.RationalPoly{NT, DT}}) where {T, NT, DT} at /Users/joey/.julia/packages/MultivariatePolynomials/XkHHT/src/promote.jl:33
  promote_rule_constant(::Type{S}, ::Type{DynamicPolynomials.Term{C, T}}) where {S, C, T} at /Users/joey/.julia/packages/DynamicPolynomials/5p8Ga/src/promote.jl:2
  promote_rule_constant(::Type{S}, ::Type{var"#s20"} where var"#s20"<:Union{DynamicPolynomials.Polynomial{C, T}, DynamicPolynomials.Term{C, T}}) where {S, C, T} at /Users/joey/.julia/packages/DynamicPolynomials/5p8Ga/src/promote.jl:3
  ...
Stacktrace:
  [1] promote_rule(#unused#::Type{Int64}, #unused#::Type{MultivariatePolynomials.AbstractPolynomialLike})
    @ MultivariatePolynomials ~/.julia/packages/MultivariatePolynomials/XkHHT/src/promote.jl:30
  [2] promote_type
    @ ./promotion.jl:233 [inlined]
  [3] promote_eltypeof(v1::Int64, vs::Vector{MultivariatePolynomials.AbstractPolynomialLike})
    @ Base ./abstractarray.jl:1460
  [4] _cat(::Val{1}, ::Int64, ::Vararg{Any, N} where N)
    @ Base ./abstractarray.jl:1641
  [5] #cat#129
    @ ./abstractarray.jl:1781 [inlined]
  [6] vcat
    @ ./abstractarray.jl:1710 [inlined]
  [7] arguments(a::SymbolicUtils.Add{Real, Int64, Dict{Any, Number}, Nothing})
    @ SymbolicUtils ~/.julia/packages/SymbolicUtils/BKm0J/src/types.jl:630
  [8] (::SymbolicUtils.Rewriters.Walk{:post, SymbolicUtils.Rewriters.PassThrough{SymbolicUtils.Rewriters.RestartedChain{Vector{SymbolicUtils.AbstractRule}}}, SymbolicUtils.var"#simterm#125"{SymbolicUtils.var"#simterm#114#126"}, false})(x::SymbolicUtils.Add{Real, Int64, Dict{Any, Number}, Nothing})
    @ SymbolicUtils.Rewriters ~/.julia/packages/SymbolicUtils/BKm0J/src/rewriters.jl:142
  [9] (::SymbolicUtils.Rewriters.Fixpoint{SymbolicUtils.Rewriters.Walk{:post, SymbolicUtils.Rewriters.PassThrough{SymbolicUtils.Rewriters.RestartedChain{Vector{SymbolicUtils.AbstractRule}}}, SymbolicUtils.var"#simterm#125"{SymbolicUtils.var"#simterm#114#126"}, false}})(x::SymbolicUtils.Add{Real, Int64, Dict{Any, Number}, Nothing})
    @ SymbolicUtils.Rewriters ~/.julia/packages/SymbolicUtils/BKm0J/src/rewriters.jl:104
 [10] to_mpoly(t::Term{Real, Nothing}, variable_type::Type, dicts::Tuple{OrderedCollections.OrderedDict{Sym, Any}, OrderedCollections.OrderedDict{Any, Sym}})
    @ SymbolicUtils ~/.julia/packages/SymbolicUtils/BKm0J/src/abstractalgebra.jl:91
 [11] to_mpoly(t::Term{Real, Nothing})
    @ SymbolicUtils ~/.julia/packages/SymbolicUtils/BKm0J/src/abstractalgebra.jl:78
 [12] iszero(x::Num)
    @ Symbolics ~/.julia/packages/Symbolics/AuRtL/src/num.jl:114
 [13] _iszero(x::Num)
    @ SparseArrays.HigherOrderFns /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/higherorderfns.jl:202
 [14] _map_zeropres!(f::Symbolics.var"#256#257"{Symbol}, C::SparseArrays.SparseMatrixCSC{Any, Int64}, A::SparseArrays.SparseMatrixCSC{Num, Int64})
    @ SparseArrays.HigherOrderFns /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/higherorderfns.jl:248
 [15] _noshapecheck_map(::Symbolics.var"#256#257"{Symbol}, ::SparseArrays.SparseMatrixCSC{Num, Int64})
    @ SparseArrays.HigherOrderFns /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/higherorderfns.jl:166
 [16] map
    @ /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.6/SparseArrays/src/higherorderfns.jl:1154 [inlined]
 [17] _make_array(rhss::SparseArrays.SparseMatrixCSC{Num, Int64}, similarto::Symbol)
    @ Symbolics ~/.julia/packages/Symbolics/AuRtL/src/build_function.jl:285
 [18] make_array(s::Symbolics.SerialForm, dargs::Vector{Any}, arr::SparseArrays.SparseMatrixCSC{Num, Int64}, similarto::Symbol)
    @ Symbolics ~/.julia/packages/Symbolics/AuRtL/src/build_function.jl:246
 [19] _build_function(::Symbolics.JuliaTarget, ::SparseArrays.SparseMatrixCSC{Num, Int64}, ::Vector{Term{Real, Nothing}}, ::Vararg{Any, N} where N; expression::Type, expression_module::Module, checkbounds::Bool, linenumbers::Bool, outputidxs::Nothing, skipzeros::Bool, wrap_code::Tuple{Nothing, Nothing}, fillzeros::Bool, parallel::Symbolics.SerialForm, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ Symbolics ~/.julia/packages/Symbolics/AuRtL/src/build_function.jl:219
 [20] build_function(::SparseArrays.SparseMatrixCSC{Num, Int64}, ::Vararg{Any, N} where N; target::Symbolics.JuliaTarget, kwargs::Base.Iterators.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:expression, :expression_module, :checkbounds), Tuple{DataType, Module, Bool}}})
    @ Symbolics ~/.julia/packages/Symbolics/AuRtL/src/build_function.jl:52
 [21] generate_jacobian(sys::ODESystem, dvs::Vector{Term{Real, Nothing}}, ps::Vector{Sym{Real, Base.ImmutableDict{DataType, Any}}}; simplify::Bool, sparse::Bool, kwargs::Base.Iterators.Pairs{Symbol, Any, Tuple{Symbol, Symbol, Symbol}, NamedTuple{(:expression, :expression_module, :checkbounds), Tuple{DataType, Module, Bool}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/MYBNE/src/systems/diffeqs/abstractodesystem.jl:50
 [22] (ODEFunction{true, F, TMM, Ta, Tt, TJ, JVP, VJP, JP, SP, TW, TWt, TPJ, S, S2, O, TCV} where {F, TMM, Ta, Tt, TJ, JVP, VJP, JP, SP, TW, TWt, TPJ, S, S2, O, TCV})(sys::ODESystem, dvs::Vector{Term{Real, Nothing}}, ps::Vector{Sym{Real, Base.ImmutableDict{DataType, Any}}}, u0::Nothing; version::Nothing, tgrad::Bool, jac::Bool, eval_expression::Bool, sparse::Bool, simplify::Bool, eval_module::Module, steady_state::Bool, checkbounds::Bool, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/MYBNE/src/systems/diffeqs/abstractodesystem.jl:187
 [23] #ODEFunction#143
    @ ~/.julia/packages/ModelingToolkit/MYBNE/src/systems/diffeqs/abstractodesystem.jl:140 [inlined]
 [24] top-level scope
    @ REPL[48]:1
ChrisRackauckas commented 3 years ago

@shashi @YingboMa where are these polynomials coming in?