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