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[1]=x[2]
    dx[3]=x[4]
    dx[2]=(-x[2]^2/2*(3-x[2]/c) +(1+(1-3*κ)*x[2]/c)*((P_stat-P_v)/ρ+2*σ/(ρ*p[1]))*abs(complex(p[1]/x[1]))^(3.0*κ)
        -2*σ/(ρ*x[1]) -4*µ*x[2]/(ρ*x[1]) -(1+x[2]/c)*(P_stat-P_v+p[3]*sin(2*pi*p[4]*(t+x[1]/c)))/ρ -
        1/p[5]*(2*x[3]*x[4]^2+x[3]^2*dx[4]) )/((1-x[2]/c)*x[1]+4*µ/(ρ*c))

    dx[4]=(-x[4]^2/2*(3-x[4]/c) +(1+(1-3*κ)*x[4]/c)*((P_stat-P_v)/ρ+2*σ/(ρ*p[2]))*abs(complex(p[2]/x[3]))^(3.0*κ)
        -2*σ/(ρ*x[3]) -4*µ*x[4]/(ρ*x[3]) -(1+x[4]/c)*(P_stat-P_v+p[3]*sin(2*pi*p[4]*(t+x[3]/c)))/ρ -
        1/p[5]*(2*x[1]*x[2]^2+x[1]^2*dx[2]) )/((1-x[4]/c)*x[3]+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[1], 0.0, stats[2], 0.0]
    tspan = (0.0,stats[6])
    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=>RadauIIA5()),
          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[10],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:
  JULIA_NUM_THREADS = 2

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

Thank you for your effort.
I see that you are fixing the caching problem.
But in terms of stiffness detection/estimation are there different methods?
I would like to show in which regions it behaves stiff.

I worked around it a little bit with trying to solve with a nonstiff alg and when integration wasn’t a success trying again with a stiff alg. But I don’t know how reliable this is to conclude stiffness.

Not yet. Only one detection method for now, but open an issue. If you open an issue with some good examples to work on, then I’ll make this one of the student projects over the next year or summer. We know of some alternative methods we could use but it’s really hard to get a good set of test problems for this, so we’d be extremely grateful to have that it would really help us push forwards.