search

jonniedie / ComponentArrays.jl

Arrays with arbitrarily nested named components.
MIT License
287 stars 34 forks source link

Component Arrays with DAE's #128

Open aabills opened 2 years ago

aabills commented 2 years ago

Hi! My team and I stumbled upon this package and have found it really useful in our current project. Unfortunately though, we're having trouble getting it to play nicely with much of the DAE ecosystem in DifferentialEquations.jl, and in particular with linear solvers and nonlinear solvers. Here're a few MW examples:

Top of file:

using OrdinaryDiffEq
using ComponentArrays
using Parameters
using SciMLNLSolve
using LinearSolve

p = ComponentArray(a=0.04,b=1e4,c=3e7,d=1.0)
u0 = ComponentArray(x1=1.0,x2=0.0,x3=0.0)

function f(resid,du,u,p,t)
    @unpack a,b,c,d=p
    resid[1] = - a*u[1] + b*u[2]*u[3] - du[1]
    resid[2] = + a*u[1] - c*u[2]^2 - b*u[2]*u[3] - du[2]
    resid[3] = u[1] + u[2] + u[3] - d
end

du0 = similar(u0)
fill!(du0,0.0)

tspan = (0.0,1e5)

prob = DAEProblem(f,du0,u0,tspan,p,differential_vars=[true,true,false])

This all works fine. Unfortunately, when passed into a solver, sol = solve(prob,DFBDF()) The default linear solver in this case will be KrylovJL_GMRES() and it will error with this stack trace:

ERROR: MethodError: Cannot `convert` an object of type UndefInitializer to an object of type Vector{Float64}
Closest candidates are:
  convert(::Type{Array{T, N}}, ::FillArrays.Zeros{V, N}) where {T, V, N} at ~/.julia/packages/FillArrays/5Arin/src/FillArrays.jl:440
  convert(::Type{<:Array}, ::LabelledArrays.LArray) at ~/.julia/packages/LabelledArrays/qbHni/src/larray.jl:104
  convert(::Type{Array{T, N}}, ::StaticArrays.SizedArray{S, T, N, N, Array{T, N}}) where {S, T, N} at ~/.julia/packages/StaticArrays/58yy1/src/SizedArray.jl:118
  ...
Stacktrace:
  [1] ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}(data::UndefInitializer, axes::Int64)
    @ ComponentArrays ~/.julia/packages/ComponentArrays/4wSaG/src/componentarray.jl:35
  [2] Krylov.GmresSolver(n::Int64, m::Int64, memory::Int64, S::Type{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}})
    @ Krylov ~/.julia/packages/Krylov/cjUQT/src/krylov_solvers.jl:1476
  [3] Krylov.GmresSolver(A::ComponentMatrix{Float64}, b::ComponentVector{Float64}, memory::Int64)
    @ Krylov ~/.julia/packages/Krylov/cjUQT/src/krylov_solvers.jl:1494
  [4] init_cacheval(alg::LinearSolve.KrylovJL{typeof(Krylov.gmres!), Int64, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, A::ComponentMatrix{Float64}, b::ComponentVector{Float64}, u::ComponentVector{Float64}, Pl::LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}}}, Pr::LinearAlgebra.Diagonal{Float64, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}}, maxiters::Int64, abstol::Float64, reltol::Float64, verbose::Bool)
    @ LinearSolve ~/.julia/packages/LinearSolve/aLfeI/src/iterative_wrappers.jl:78
  [5] init_cacheval(alg::Nothing, A::ComponentMatrix{Float64}, b::ComponentVector{Float64}, u::ComponentVector{Float64}, Pl::LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}}}, Pr::LinearAlgebra.Diagonal{Float64, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}}, maxiters::Int64, abstol::Float64, reltol::Float64, verbose::Bool)
    @ LinearSolve ~/.julia/packages/LinearSolve/aLfeI/src/default.jl:169
  [6] init(::LinearProblem{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, true, ComponentMatrix{Float64, Matrix{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}, Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, SciMLBase.NullParameters, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::Nothing; alias_A::Bool, alias_b::Bool, abstol::Float64, reltol::Float64, maxiters::Int64, verbose::Bool, Pl::LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}}}, Pr::LinearAlgebra.Diagonal{Float64, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}}, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ LinearSolve ~/.julia/packages/LinearSolve/aLfeI/src/common.jl:97
  [7] build_nlsolver(alg::DFBDF{5, 3, true, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}, nlalg::NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, u::ComponentVector{Float64}, uprev::ComponentVector{Float64}, p::ComponentVector{Float64}, t::Float64, dt::Float64, f::Function, rate_prototype::ComponentVector{Float64}, #unused#::Type{Float64}, #unused#::Type{Float64}, #unused#::Type{Float64}, γ::Float64, c::Float64, α::Int64, #unused#::Val{true})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vGDGq/src/nlsolve/utils.jl:169
  [8] build_nlsolver
    @ ~/.julia/packages/OrdinaryDiffEq/vGDGq/src/nlsolve/utils.jl:124 [inlined]
  [9] build_nlsolver(alg::DFBDF{5, 3, true, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}, u::ComponentVector{Float64}, uprev::ComponentVector{Float64}, p::ComponentVector{Float64}, t::Float64, dt::Float64, f::Function, rate_prototype::ComponentVector{Float64}, #unused#::Type{Float64}, #unused#::Type{Float64}, #unused#::Type{Float64}, γ::Float64, c::Float64, iip::Val{true})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vGDGq/src/nlsolve/utils.jl:118
 [10] alg_cache(alg::DFBDF{5, 3, true, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}, du::ComponentVector{Float64}, u::ComponentVector{Float64}, res_prototype::ComponentVector{Float64}, rate_prototype::ComponentVector{Float64}, uEltypeNoUnits::Type, uBottomEltypeNoUnits::Type, tTypeNoUnits::Type{Float64}, uprev::ComponentVector{Float64}, uprev2::ComponentVector{Float64}, f::DAEFunction{true, typeof(f), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, t::Float64, dt::Float64, reltol::Float64, p::ComponentVector{Float64}, calck::Bool, #unused#::Val{true})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vGDGq/src/caches/dae_caches.jl:183
 [11] __init(prob::DAEProblem{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(a = 1, b = 2, c = 3, d = 4)}}}, DAEFunction{true, typeof(f), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, alg::DFBDF{5, 3, true, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}, 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 ~/.julia/packages/OrdinaryDiffEq/vGDGq/src/solve.jl:293
 [12] __init(prob::DAEProblem{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(a = 1, b = 2, c = 3, d = 4)}}}, DAEFunction{true, typeof(f), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, alg::DFBDF{5, 3, true, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}) (repeats 5 times)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vGDGq/src/solve.jl:67
 [13] __solve(::DAEProblem{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(a = 1, b = 2, c = 3, d = 4)}}}, DAEFunction{true, typeof(f), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, ::DFBDF{5, 3, true, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vGDGq/src/solve.jl:4
 [14] __solve(::DAEProblem{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(a = 1, b = 2, c = 3, d = 4)}}}, DAEFunction{true, typeof(f), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, ::DFBDF{5, 3, true, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vGDGq/src/solve.jl:4
 [15] solve_call(_prob::DAEProblem{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(a = 1, b = 2, c = 3, d = 4)}}}, DAEFunction{true, typeof(f), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, args::DFBDF{5, 3, true, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}; merge_callbacks::Bool, kwargshandle::KeywordArgError, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/202zR/src/solve.jl:384
 [16] solve_call(_prob::DAEProblem{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(a = 1, b = 2, c = 3, d = 4)}}}, DAEFunction{true, typeof(f), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, args::DFBDF{5, 3, true, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/202zR/src/solve.jl:366
 [17] solve_up(prob::DAEProblem{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(a = 1, b = 2, c = 3, d = 4)}}}, DAEFunction{true, typeof(f), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, sensealg::Nothing, u0::ComponentVector{Float64}, p::ComponentVector{Float64}, args::DFBDF{5, 0, true, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/202zR/src/solve.jl:412
 [18] solve_up
    @ ~/.julia/packages/DiffEqBase/202zR/src/solve.jl:401 [inlined]
 [19] #solve#29
    @ ~/.julia/packages/DiffEqBase/202zR/src/solve.jl:396 [inlined]
 [20] solve(prob::DAEProblem{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(a = 1, b = 2, c = 3, d = 4)}}}, DAEFunction{true, typeof(f), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, args::DFBDF{5, 0, true, Nothing, NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/202zR/src/solve.jl:391
 [21] top-level scope
    @ ~/Sandbox/solving_probs/nlsolve_CA_mwe.jl:25

This error occurs when using any solver from Krylov.jl.

Taking a look at the first line of the stack trace, it looks like Krylov is attempting to create an empty array, by a syntax which is not compatible with ComponentArrays.

I also wanted to try IterativeSolvers.jl, but unfortunately it doesn't seem to work at all with DAEProblem (see https://github.com/SciML/OrdinaryDiffEq.jl/issues/1673

Moving to different choices in nonlinear solvers, we have a different issue: solve(prob,DFBDF(NLSolveJL()))

Results in:

ERROR: MethodError: no method matching build_nlsolver(::DFBDF{5, 3, true, Nothing, NLSolveJL{Static, SciMLNLSolve.var"#2#4"}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}, ::NLSolveJL{Static, SciMLNLSolve.var"#2#4"}, ::ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, ::ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, ::ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(a = 1, b = 2, c = 3, d = 4)}}}, ::Float64, ::Float64, ::DAEFunction{true, typeof(f), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, ::ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, ::Type{Float64}, ::Type{Float64}, ::Type{Float64}, ::Float64, ::Float64, ::Int64, ::Val{true})
Closest candidates are:
  build_nlsolver(::Any, ::Union{NLAnderson, NLFunctional, NLNewton}, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, ::Any, ::Any, ::Any, ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} at ~/.julia/packages/OrdinaryDiffEq/vGDGq/src/nlsolve/utils.jl:128
  build_nlsolver(::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, ::Any, ::Any, ::Any, ::Any) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} at ~/.julia/packages/OrdinaryDiffEq/vGDGq/src/nlsolve/utils.jl:122
  build_nlsolver(::Any, ::Union{NLAnderson, NLFunctional, NLNewton}, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, ::Type{uEltypeNoUnits}, ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, ::Any, ::Any, ::Any, ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits} at ~/.julia/packages/OrdinaryDiffEq/vGDGq/src/nlsolve/utils.jl:202
  ...
Stacktrace:
  [1] build_nlsolver(alg::DFBDF{5, 3, true, Nothing, NLSolveJL{Static, SciMLNLSolve.var"#2#4"}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}, u::ComponentVector{Float64}, uprev::ComponentVector{Float64}, p::ComponentVector{Float64}, t::Float64, dt::Float64, f::Function, rate_prototype::ComponentVector{Float64}, #unused#::Type{Float64}, #unused#::Type{Float64}, #unused#::Type{Float64}, γ::Float64, c::Float64, α::Int64, iip::Val{true})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vGDGq/src/nlsolve/utils.jl:124
  [2] build_nlsolver(alg::DFBDF{5, 3, true, Nothing, NLSolveJL{Static, SciMLNLSolve.var"#2#4"}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}, u::ComponentVector{Float64}, uprev::ComponentVector{Float64}, p::ComponentVector{Float64}, t::Float64, dt::Float64, f::Function, rate_prototype::ComponentVector{Float64}, #unused#::Type{Float64}, #unused#::Type{Float64}, #unused#::Type{Float64}, γ::Float64, c::Float64, iip::Val{true})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vGDGq/src/nlsolve/utils.jl:118
  [3] alg_cache(alg::DFBDF{5, 3, true, Nothing, NLSolveJL{Static, SciMLNLSolve.var"#2#4"}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}, du::ComponentVector{Float64}, u::ComponentVector{Float64}, res_prototype::ComponentVector{Float64}, rate_prototype::ComponentVector{Float64}, uEltypeNoUnits::Type, uBottomEltypeNoUnits::Type, tTypeNoUnits::Type{Float64}, uprev::ComponentVector{Float64}, uprev2::ComponentVector{Float64}, f::DAEFunction{true, typeof(f), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, t::Float64, dt::Float64, reltol::Float64, p::ComponentVector{Float64}, calck::Bool, #unused#::Val{true})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vGDGq/src/caches/dae_caches.jl:183
  [4] __init(prob::DAEProblem{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(a = 1, b = 2, c = 3, d = 4)}}}, DAEFunction{true, typeof(f), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, alg::DFBDF{5, 3, true, Nothing, NLSolveJL{Static, SciMLNLSolve.var"#2#4"}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}, 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 ~/.julia/packages/OrdinaryDiffEq/vGDGq/src/solve.jl:293
  [5] __init(prob::DAEProblem{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(a = 1, b = 2, c = 3, d = 4)}}}, DAEFunction{true, typeof(f), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, alg::DFBDF{5, 3, true, Nothing, NLSolveJL{Static, SciMLNLSolve.var"#2#4"}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}, timeseries_init::Tuple{}, ts_init::Tuple{}, ks_init::Tuple{}, recompile::Type{Val{true}}) (repeats 5 times)
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vGDGq/src/solve.jl:67
  [6] __solve(::DAEProblem{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(a = 1, b = 2, c = 3, d = 4)}}}, DAEFunction{true, typeof(f), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, ::DFBDF{5, 3, true, Nothing, NLSolveJL{Static, SciMLNLSolve.var"#2#4"}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vGDGq/src/solve.jl:4
  [7] __solve(::DAEProblem{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(a = 1, b = 2, c = 3, d = 4)}}}, DAEFunction{true, typeof(f), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, ::DFBDF{5, 3, true, Nothing, NLSolveJL{Static, SciMLNLSolve.var"#2#4"}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing})
    @ OrdinaryDiffEq ~/.julia/packages/OrdinaryDiffEq/vGDGq/src/solve.jl:4
  [8] solve_call(_prob::DAEProblem{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(a = 1, b = 2, c = 3, d = 4)}}}, DAEFunction{true, typeof(f), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, args::DFBDF{5, 3, true, Nothing, NLSolveJL{Static, SciMLNLSolve.var"#2#4"}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}; merge_callbacks::Bool, kwargshandle::KeywordArgError, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/202zR/src/solve.jl:384
  [9] solve_call(_prob::DAEProblem{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(a = 1, b = 2, c = 3, d = 4)}}}, DAEFunction{true, typeof(f), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, args::DFBDF{5, 3, true, Nothing, NLSolveJL{Static, SciMLNLSolve.var"#2#4"}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/202zR/src/solve.jl:366
 [10] solve_up(prob::DAEProblem{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(a = 1, b = 2, c = 3, d = 4)}}}, DAEFunction{true, typeof(f), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, sensealg::Nothing, u0::ComponentVector{Float64}, p::ComponentVector{Float64}, args::DFBDF{5, 0, true, Nothing, NLSolveJL{Static, SciMLNLSolve.var"#2#4"}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/202zR/src/solve.jl:412
 [11] solve_up
    @ ~/.julia/packages/DiffEqBase/202zR/src/solve.jl:401 [inlined]
 [12] #solve#29
    @ ~/.julia/packages/DiffEqBase/202zR/src/solve.jl:396 [inlined]
 [13] solve(prob::DAEProblem{ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(x1 = 1, x2 = 2, x3 = 3)}}}, Tuple{Float64, Float64}, true, ComponentVector{Float64, Vector{Float64}, Tuple{Axis{(a = 1, b = 2, c = 3, d = 4)}}}, DAEFunction{true, typeof(f), Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, Vector{Bool}}, args::DFBDF{5, 0, true, Nothing, NLSolveJL{Static, SciMLNLSolve.var"#2#4"}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/202zR/src/solve.jl:391
 [14] top-level scope
    @ ~/Sandbox/solving_probs/nlsolve_CA_mwe.jl:25

This is true for any NLSolve solver. Some things do work! For example, solve(prob,DFBDF(linsolve=LUFactorization()) works fine. However, it'd be amazing to have more compatibility of DAE's and ComponentArrays.

jonniedie commented 2 years ago

Hey Alec. Yeah, I think I need to explicitly add some of the factorizations. A lot of them don’t work because they call a similar method that is fundamentally type-unstable with ComponentArrays, StaticArrays, and similar types of arrays that encode their size (in one way or another) in the type domain. The approach for ComponentArrays and StaticArrays is just to fall back to a plain Matrix in these cases, which tends to break things down the line that expect them to be the same type as what went into the factorization. So I think that’s what’s happening here.

I’ll take a look at it when I get the chance. Unfortunately I’ve kinda got a backlog of issues to fix this week and haven’t had time to commit to them.