# ODE AutoTsit5(Rodas5()) switching bug (depends from tolerance history)

I want to solve my coupled DiffEq with different algorithms, because i don’t know the specific characteristics of it. But i know that it can exhibit stiff behaviour in some cases. I need to solve it many times in a parameter space, so efficiency is key and it would be nice to have a good stiffness detection and automated switching. This is some of the code:

``````using Plots, DifferentialEquations;
##  Liquid Properties   ##
const c=1500.0        #[m/s] sound velocity in the liquid
const κ=1.4           #polytropic exponent of the gas in the bubble
const P_stat=100e+3   #[Pa] static pressure
const P_v=0.          #[Pa] vapor pressure
const ρ=998.0         #[kg/m^3] density of liquid
const σ=0.0725        #[N/m] surface tension
const µ= 0.001        #[N s/m^2] viscosity

function bubbleKME2w(dx,x,p,t) #mymodel
dx=x
dx=x
dx=(-x^2/2*(3-x/c) +(1+(1-3*κ)*x/c)*((P_stat-P_v)/ρ+2*σ/(ρ*p))*abs(complex(p/x))^(3.0*κ)
-2*σ/(ρ*x) -4*µ*x/(ρ*x) -(1+x/c)*(P_stat-P_v+p*sin(2*pi*p*(t+x/c)))/ρ -
1/p*(2*x*x^2+x^2*dx) )/((1-x/c)*x+4*µ/(ρ*c))

dx=(-x^2/2*(3-x/c) +(1+(1-3*κ)*x/c)*((P_stat-P_v)/ρ+2*σ/(ρ*p))*abs(complex(p/x))^(3.0*κ)
-2*σ/(ρ*x) -4*µ*x/(ρ*x) -(1+x/c)*(P_stat-P_v+p*sin(2*pi*p*(t+x/c)))/ρ -
1/p*(2*x*x^2+x^2*dx) )/((1-x/c)*x+4*µ/(ρ*c))
end
function plot_probw(stats, alg, tol)
#=  R1_n, R2_n, P_s, v, d , tend= stats
1     2      3   4  5 ,  6  =#
x0 = [stats, 0.0, stats, 0.0]
tspan = (0.0,stats)
prob = ODEProblem(bubbleKME2w, x0, tspan, stats)
sol = solve(prob, alg[:alg], reltol=tol, abstol=tol, maxiters=1e7)
return sol
end
setups = [
Dict(:alg=>Rosenbrock23()),#1
Dict(:alg=>Tsit5()),#3
Dict(:alg=>AutoTsit5(Rodas4(),nonstifftol=7/10,stifftol=9/10)),
Dict(:alg=>AutoVern7(KenCarp4(),nonstifftol=7/10,stifftol=9/10)),#5
Dict(:alg=>VCAB4()),
Dict(:alg=>ABDF2()),#7
Dict(:alg=>Feagin12()),
Dict(:alg=>AutoTsit5(ABDF2(),nonstifftol=7/10,stifftol=9/10)),#9
Dict(:alg=>AutoTsit5(Rodas5(),nonstifftol=7/10,stifftol=9/10))
]
sol=plot_probw([0.5e-6, 1.4e-6, 132e+3,20e+3, 200e-6,300e-6],setups,1e-8)
stifftimes=length(findall(x-> x >1, sol.alg_choice)) #gets the number of steps with the stiff method
``````

In this case it runs with the `AutoTsit5(Rodas5(),nonstifftol=7/10,stifftol=9/10))` algorithm with tolerance `1e-8` and it is not switching to the stiff method. I get same behaviour for tolerances `<=1e-6`. For a tolerance of `1e-5` it detects stiffness and switches to `Rodas5`.
After that I get stiffness detects and switches for the lower tolerances and even other starting values, which were previously classified as non-stiff, too!
After reassigning my dictionary of algorithms I get again no stiff detects for tolerances `<=1e-6` until I run with tolerance `>=1e-5` once.

So it happens that the algorithm is switching to the stiff method even cases where it’s not necessary. This slows down the computation a lot.

I am working with Atom and:

``````Julia Version 1.5.0
Commit 96786e22cc (2020-08-01 23:44 UTC)
Platform Info:
OS: Windows (x86_64-w64-mingw32)
CPU: Intel(R) Core(TM) i5-7440HQ CPU @ 2.80GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-9.0.1 (ORCJIT, skylake)
Environment:
``````

See https://github.com/SciML/OrdinaryDiffEq.jl/blob/a9162c40904ebdde1d5a435137211a0495500fbe/src/composite_algs.jl#L31-L40. It looks like the stiffness detection has state.

Could you open an issue? This seems like it’s “mostly” tuning. The current stiffness detection method that we’re using is an extremely cheap heuristic, so it works best as a default method switch once or twice in a whole integration and get close to the efficiency of having directly used a method without detection. There are more intrusive methods that can more accurately calculate the cutoff to make switching back and forth easier, but of course they have a higher cost in the eigenvalue estimates. Still, this seems like an interesting problem to test it on so it would be good to get an issue in OrdinaryDiffEq.jl.

@YingboMa we should move that caching to the integrator.

Yes, that should be an easy fix.

Does this mean that the cache will be reset for each solve? Because if I’ve understood it correctly, the problem is that all those different solve calls that uses AutoTsit5 ends up sharing counters etc. Very dangerous if you put thread a number of solve calls!

Yes, I don’t know how this got closed: https://github.com/SciML/OrdinaryDiffEq.jl/issues/325 that was an accident. But @YingboMa is finally getting around to it: https://github.com/SciML/OrdinaryDiffEq.jl/pull/1258

1 Like