ggebbie / TMItransient.jl

Transient simulations with the TMI circulation matrix
MIT License
0 stars 1 forks source link

`vintagedistribution` encounters an instability for the 2deg TMI #21

Closed b-r-hamilton closed 1 year ago

b-r-hamilton commented 1 year ago

The following line runs forever, and then encounters an instability

julia> @time g = vintagedistribution(TMIversion, γ, L, B, 2000, 2010, tmodern = 2015) 
┌ Warning: Instability detected. Aborting
└ @ SciMLBase ~/.julia/packages/SciMLBase/VdcHg/src/integrator_interface.jl:596
ERROR: BoundsError: attempt to access 1-element Vector{Vector{Float64}} at index [2]
Stacktrace:
 [1] getindex
   @ ./essentials.jl:13 [inlined]
 [2] getindex(A::SciMLBase.ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, SciMLBase.ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, SciMLBase.NullParameters, SciMLBase.ODEFunction{true, SciMLBase.FullSpecialize, TMItransient.var"#f#18"{SparseMatrixCSC{Float64, Int64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, OrdinaryDiffEq.QNDF{5, 12, true, LinearSolve.UMFPACKFactorization, OrdinaryDiffEq.NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing, Nothing, Nothing, NTuple{5, Float64}}, OrdinaryDiffEq.InterpolationData{SciMLBase.ODEFunction{true, SciMLBase.FullSpecialize, TMItransient.var"#f#18"{SparseMatrixCSC{Float64, Int64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, OrdinaryDiffEq.QNDFCache{5, StaticArraysCore.SMatrix{5, 5, Float64, 25}, Matrix{Float64}, Vector{Float64}, OrdinaryDiffEq.NLSolver{OrdinaryDiffEq.NLNewton{Rational{Int64}, Rational{Int64}, Rational{Int64}, Rational{Int64}}, true, Vector{Float64}, Float64, Nothing, Float64, OrdinaryDiffEq.NLNewtonCache{Vector{Float64}, Float64, Float64, Vector{Float64}, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, SciMLBase.UJacobianWrapper{SciMLBase.ODEFunction{true, SciMLBase.FullSpecialize, TMItransient.var"#f#18"{SparseMatrixCSC{Float64, Int64}}, UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, SparseMatrixCSC{Float64, Int64}, SparseMatrixCSC{Float64, Int64}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing, Nothing}, Float64, SciMLBase.NullParameters}, SparseDiffTools.ForwardColorJacCache{Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqBase.OrdinaryDiffEqTag, Float64}, Float64, 12}}, Vector{Float64}, Vector{Vector{NTuple{12, Float64}}}, Vector{Int64}, SparseMatrixCSC{Float64, Int64}}, LinearSolve.LinearCache{SparseMatrixCSC{Float64, Int64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, LinearSolve.UMFPACKFactorization, SparseArrays.UMFPACK.UmfpackLU{Float64, Int64}, LinearSolve.InvPreconditioner{Diagonal{Float64, Vector{Float64}}}, Diagonal{Float64, Vector{Float64}}, Float64, true}}}, Matrix{Float64}, Float64, Float64, StaticArraysCore.SVector{5, Float64}, Vector{Float64}, Vector{Float64}}}, DiffEqBase.DEStats, Nothing}, I::Int64)
   @ SciMLBase ~/.julia/packages/SciMLBase/VdcHg/src/solutions/solution_interface.jl:42
 [3] vintagedistribution(TMIversion::String, γ::Grid, L::SparseMatrixCSC{Float64, Int64}, B::SparseMatrixCSC{Float64, Int64}, t₀::Int64, tf::Int64; tmodern::Int64)
   @ TMItransient ~/Code/TMItransient.jl/src/TMItransient.jl:374
 [4] top-level scope
   @ ./timing.jl:273 [inlined]
 [5] top-level scope
   @ ./REPL[31]:0
ggebbie commented 1 year ago

is it related to Issue #17 ?

The solution was to pin OrdinaryDiffEq to 6.33.3.

I am trying to automate the solution in the Manifest now.

b-r-hamilton commented 1 year ago

Hm this might make sense. I can run this line with no issues in the TMItransient.jl project, which has ODE 6.33.3, but I tried to add TMItransient to a different project and it won't work there.

ggebbie commented 1 year ago

Try updating to latest main branch

b-r-hamilton commented 1 year ago

Revisiting this - I've tried updating to main branch, still throws the error. My current soluiton is to install OrdinaryDiffEq into my package (even though it is not actually needed) and pin it to 6.33.3

ggebbie commented 1 year ago

My attempted solution was to add a compat entry like this: https://github.com/ggebbie/TMItransient.jl/blob/62675cb3ff95a4f165a081ff4f3c2a0354076a50/Project.toml#L22

If that doesn't solve it, here's another idea. Is your environment using DifferentialEquations.jl? I think it may conflict with OrdinaryDiffEq.jl and that only the latter is probably necessary.

b-r-hamilton commented 1 year ago

Hm, I didn't have DiffEq installed, but I did have Sundials installed - removed that and it seems to work fine now. Not 100% that was really the issue, but it's working fine now.