AutoSwitch algorithm become a bit too sensitive

[LEquinoxy]

When solving the same ODE (a 100-dimensional strongly slow-fast system, where the non-smooth component is approximated using hyperbolic tangent function ) problem using DifferentialEquations.jl , I noticed that ODE solving is fast on [email protected] but significantly slower on later versions

Here are the test results:
[email protected]:

AutoVern6(Rodas5P(),nonstifftol = 15/ 10):
BenchmarkTools.Trial: 10 samples with 1 evaluation.
Range (min … max):  1.418 s …    2.346 s  ┊ GC (min … max):  0.00% … 41.36%
Memory estimate: 661.43 MiB, allocs estimate: 2418966.

Vern6():
BenchmarkTools.Trial: 10 samples with 1 evaluation.
Range (min … max):  1.451 s …    1.982 s  ┊ GC (min … max):  0.00% … 26.63%
Memory estimate: 659.45 MiB, allocs estimate: 2418817.

Rodas5P():
BenchmarkTools.Trial: 4 samples with 1 evaluation.
Range (min … max):  22.905 s …   24.300 s  ┊ GC (min … max): 2.90% … 2.82%
Memory estimate: 1.88 GiB, allocs estimate: 21146606.

[email protected]:

AutoVern6(Rodas5P(),nonstifftol = 15/ 10):
BenchmarkTools.Trial: 2 samples with 1 evaluation.
Range (min … max):  20.355 s …  20.452 s  ┊ GC (min … max): 7.00% … 7.22%
 Memory estimate: 1.52 GiB, allocs estimate: 15735599.

Vern6():
BenchmarkTools.Trial: 10 samples with 1 evaluation.
Range (min … max):  1.651 s …    2.611 s  ┊ GC (min … max):  9.16% … 44.33%
Memory estimate: 660.53 MiB, allocs estimate: 3311381.

Rodas5P():
BenchmarkTools.Trial: 1 sample with 1 evaluation.
Single result which took 26.859 s (5.53% GC) to evaluate,
with a memory estimate of 1.81 GiB, over 19996952 allocations.

Here is the code, the ODE ODE_SOC_PO and initial state u0 can be found in the attachment.

using Revise, Parameters, Plots, BifurcationKit, JLD2, DifferentialEquations,ParameterizedFunctions,LinearAlgebra,MATLAB
const BK = BifurcationKit
using Arpack,BenchmarkTools
#ODE_Setting
RT=1e-8;AT=1e-10
include("ODE_SOC_PO.jl")
file = jldopen("u0_new2.jld2", "r")
u0 = file["u0_new"]
close(file)
# u0=[zeros(25,1);1;0;zeros(8,1);zeros(38,1)];
# tsidas_alg = AutoVern6(Rodas5P(),nonstifftol = 15/ 10)# Decreasing nonstifftol makes switching to the stiff algorithm more likely
tsidas_alg = Vern6()
# tsidas_alg = Rodas5P()
par_ODE = (AA_i=5e4, AA_firing=5e5, AA_ctrl=5e2, Ihold_pu=2.5e-4, 
           Us_pu_rec=1.0, Kp_dc=1.0989, Ki_dc=1/0.01092, Kp_PLL_rec=10.0, Ki_PLL_rec=1/0.02,
           Us_pu_inv=1.0, Kp_CC=0.5,    Ki_CC=10.0,      Kp_PLL_inv=10.0, Ki_PLL_inv=1/0.02,
           Kp_udc=2.0,Ki_udc=1/0.01,Kp_cQ=0.0,Ki_cQ=1/0.05,
           Lfault_real=0.1,Uref=1.0)
prob_ODE = ODEProblem(ODE_SOC_PO, u0, (3.0,4), par_ODE, reltol=RT, abstol=AT)
# sol = DifferentialEquations.solve(prob_ODE,progress = true,alg=tsidas_alg)
# sol = DifferentialEquations.solve(prob_ODE,progress = true,alg=tsidas_alg)
sol =@time DifferentialEquations.solve(prob_ODE,progress = true,alg=tsidas_alg)

Desktop.zip
https://discourse.julialang.org/t/ode-solving-is-fast-on-differentialequation-7-9-0-but-significantly-slower-on-later-versions/124377/11

[oscardssmith]

One possibility of what's going on here is that the stiffness calculation takes into account the dt which at the start of the solve may not be settled in to the correct value.

[ChrisRackauckas]

Maybe forcing the solver to take a few non-stiff steps before starting the check could make sense, so it "gets into the flow" a bit and then asks "okay, is it better to make the switch given what I have seen in a few steps?" Since indeed with the initial dt at the initial point it might not be complete information to make the switch, and switching even just once is a big effect on performance so you really want to make sure you actually need to.

[oscardssmith]

What computer are you running this on? I'm seeing fairly different timings from you. Specifically, I'm seeing similar times for Rodas5P, but for me, Vern6 is 10 seconds rather than 2 seconds. Also, FBDF seems to be doing quite well here:

julia> @time sol =solve(prob, Vern6());
 11.217420 seconds (24.43 M allocations: 4.860 GiB, 17.26% gc time, 0.41% compilation time)

julia> @time sol =solve(prob, Rodas5P());
 28.194117 seconds (21.08 M allocations: 1.973 GiB, 3.08% gc time)

julia> @time sol =solve(prob, FBDF());
  6.642380 seconds (6.23 M allocations: 982.093 MiB, 9.08% gc time, 0.68% compilation time)

[oscardssmith]

Also using sparsity to make the jacobian generation faster as discussed in https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/#Choosing-Jacobian-Types seems to help a decent amount.

julia> using SparseConnectivityTracer, ADTypes, LinearSolve

julia> jac_sparsity = ADTypes.jacobian_sparsity(
           (du, u) -> ODE_SOC_PO(du, u, par_ODE, 0.0), du0, u0, detector)
73×73 SparseArrays.SparseMatrixCSC{Bool, Int64} with 638 stored entries:
⎡⠑⠀⢄⠀⠀⠀⢕⠀⠀⠀⠀⠀⠘⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠑⠑⢄⠀⢄⢕⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⢄⢄⢕⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠑⠑⠑⠑⠀⠑⠑⣄⠀⠑⠘⠃⠜⠂⠀⠀⠀⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⣿⣤⢺⣤⡄⠀⣼⡇⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡄⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠶⣿⢸⠉⢰⡶⣿⡇⠀⠀⠰⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠉⠈⠉⠉⠉⢹⣇⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠁⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠚⠄⡼⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠒⠐⠂⢐⡒⠒⠋⠛⢄⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⢀⣀⣀⣀⣀⡀⠀⢀⡀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⡇⠀⠀⠀⠈⢹⣿⠛⠭⣉⠉⣿⣿⣿⣿⢸⣿⣿⣿⣿⣷⠀⠈⠁⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⡇⠀⠀⠀⠘⣼⣿⠀⠀⠀⠀⣿⣿⣿⣿⢸⣿⣿⣿⣿⡇⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⠃⠀⠀⠀⠸⡛⠛⠀⠀⠀⢄⠛⠛⠛⠛⠘⠛⠛⠛⠛⠃⠀⠀⠄⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⡇⠀⠀⠀⢸⡇⠀⠀⠀⠀⣿⣵⢟⠀⠀⢸⠀⠀⠀⠀⠀⠀⢸⡇⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⡇⠀⠀⠀⢸⡇⠀⠀⠀⠀⠀⠀⠀⣵⢟⢸⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⠃⠀⠀⠀⠘⠣⠤⠀⠀⠀⠛⠀⠛⠀⠀⠚⣄⠀⠀⠀⠠⠀⠘⠃⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠅⠀⠪⡀⠀⠒⠅⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠰⠆⠀⠀⠀⠀⢰⣶⠀⠀⠀⠀⠀⠀⠀⠈⠰⠀⠀⠀⠀⠀⡀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⠀⠀⠀⠀⢀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠫⢄⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⎦

julia> f = ODEFunction(ODE_SOC_PO, jac_prototype = float.(jac_sparsity));

julia> prob_sparse = ODEProblem(f, u0, (3,4), par_ODE, reltol=1e-8, abstol=1e-10);

julia> @time sol2 =solve(prob_sparse, FBDF(linsolve=KLUFactorization()));
  5.204844 seconds (4.74 M allocations: 1.409 GiB, 14.34% gc time)

[LEquinoxy]

a HP omen16 laptop with i7-14650HX, RAM=16GB. I also run the code in my desktops([email protected] RAM=16GB), the results are quite the same.

https://www.amazon.com/HP-Display-i7-14650HX-Dedicated-16-ae0100nr/dp/B0D1575NN4

[LEquinoxy]

Also using sparsity to make the jacobian generation faster as discussed in https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/#Choosing-Jacobian-Types seems to help a decent amount.

julia> using SparseConnectivityTracer, ADTypes

julia> jac_sparsity = ADTypes.jacobian_sparsity(
           (du, u) -> ODE_SOC_PO(du, u, par_ODE, 0.0), du0, u0, detector)
73×73 SparseArrays.SparseMatrixCSC{Bool, Int64} with 638 stored entries:
⎡⠑⠀⢄⠀⠀⠀⢕⠀⠀⠀⠀⠀⠘⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠑⠑⢄⠀⢄⢕⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⢄⢄⢕⢄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠑⠑⠑⠑⠀⠑⠑⣄⠀⠑⠘⠃⠜⠂⠀⠀⠀⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⣿⣤⢺⣤⡄⠀⣼⡇⠀⠀⢸⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡄⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠶⣿⢸⠉⢰⡶⣿⡇⠀⠀⠰⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠉⠈⠉⠉⠉⢹⣇⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠁⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠚⠄⡼⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠒⠐⠂⢐⡒⠒⠋⠛⢄⣀⣀⣀⣀⣀⣀⣀⣀⣀⣀⢀⣀⣀⣀⣀⡀⠀⢀⡀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⡇⠀⠀⠀⠈⢹⣿⠛⠭⣉⠉⣿⣿⣿⣿⢸⣿⣿⣿⣿⣷⠀⠈⠁⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⡇⠀⠀⠀⠘⣼⣿⠀⠀⠀⠀⣿⣿⣿⣿⢸⣿⣿⣿⣿⡇⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⠃⠀⠀⠀⠸⡛⠛⠀⠀⠀⢄⠛⠛⠛⠛⠘⠛⠛⠛⠛⠃⠀⠀⠄⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⡇⠀⠀⠀⢸⡇⠀⠀⠀⠀⣿⣵⢟⠀⠀⢸⠀⠀⠀⠀⠀⠀⢸⡇⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢸⡇⠀⠀⠀⢸⡇⠀⠀⠀⠀⠀⠀⠀⣵⢟⢸⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⠃⠀⠀⠀⠘⠣⠤⠀⠀⠀⠛⠀⠛⠀⠀⠚⣄⠀⠀⠀⠠⠀⠘⠃⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠅⠀⠪⡀⠀⠒⠅⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠰⠆⠀⠀⠀⠀⢰⣶⠀⠀⠀⠀⠀⠀⠀⠈⠰⠀⠀⠀⠀⠀⡀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⠀⠀⠀⠀⢀⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠫⢄⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠁⎦

julia> f = ODEFunction(ODE_SOC_PO, jac_prototype = float.(jac_sparsity));

julia> prob_sparse = ODEProblem(f, u0, (3,4), par_ODE, reltol=1e-8, abstol=1e-10);

julia> @time sol2 =solve(prob_sparse, FBDF(linsolve=KLUFactorization()));
  5.204844 seconds (4.74 M allocations: 1.409 GiB, 14.34% gc time)

Wow! I will try this later, 😀