Merge pull request #977 from SciML/defaults

A new round of ODE defaults

Author: ChrisRackauckas

src/ode_default_alg.jl

@@ -23,67 +23,46 @@ function default_algorithm(prob::DiffEqBase.AbstractODEProblem{uType, tType, inp
     elseif typeof(prob.f) <: DynamicalODEFunction
         if tol_level == :low_tol || tol_level == :med_tol
             alg = Tsit5()
-        elseif tol_level == :high_tol
-            alg = Vern7(lazy = !callbacks)
         else
-            alg = Vern9(lazy = !callbacks)
+            alg = Vern7(lazy = !callbacks)
         end
     else # Standard ODE
-        if :nonstiff ∈ alg_hints || length(prob.u0) > 10000
+        if :nonstiff ∈ alg_hints
             # Don't default to implicit here because of memory requirements
             # And because the linear system gets unruly
             if (!(uEltype <: Float64) && !(uEltype <: Float32) && !(uEltype <: Complex)) ||
-               tol_level == :extreme_tol
+               tol_level == :extreme_tol || tol_level == :low_tol
                 # Most likely higher precision, so use a higher order method
-                alg = Vern9(lazy = !callbacks)
-            elseif tol_level == :low_tol
                 alg = Vern7(lazy = !callbacks)
-            elseif tol_level == :med_tol
+            else
                 alg = Tsit5()
-            else # tol_level == :high_tol
-                alg = BS3()
             end
         elseif :stiff ∈ alg_hints || mm # The problem is stiff
             if length(prob.u0) > 500
                 # Use Krylov method when huge!
-                if !mm
-                    alg = CVODE_BDF(linear_solver = :GMRES)
-                elseif mm
-                    alg = Rodas4(autodiff = false, linsolve = LinearSolve.KrylovJL_GMRES())
-                end
-            elseif length(prob.u0) > 50
-                if callbacks && !mm
-                    alg = CVODE_BDF()
-                elseif !callbacks
-                    alg = QNDF(autodiff = false)
-                else
-                    alg = Rodas4(autodiff = false)
-                end
+                alg = FBDF(autodiff = false, linsolve = LinearSolve.KrylovJL_GMRES())
+            elseif length(prob.u0) > 500
+                alg = FBDF(autodiff = false)
             elseif tol_level == :high_tol
                 alg = Rosenbrock23(autodiff = false)
             else
-                alg = Rodas4(autodiff = false)
+                alg = Rodas5P(autodiff = false)
             end
         else # :auto ∈ alg_hints
-            if (!(uEltype <: Float64) && !(uEltype <: Float32)) || tol_level == :extreme_tol
+            if (!(uEltype <: Float64) && !(uEltype <: Float32)) || tol_level == :extreme_tol  ||
+                tol_level == :low_tol
                 # Most likely higher precision, so use a higher order method
-                if length(prob.u0) > 100
-                    alg = AutoVern9(KenCarp47(autodiff = false), lazy = !callbacks)
-                else
-                    alg = AutoVern9(Rodas5P(autodiff = false), lazy = !callbacks)
-                end
-            elseif tol_level == :low_tol
                 if length(prob.u0) > 500
-                    alg = AutoVern7(Rodas4(autodiff = false,
-                            linsolve = LinearSolve.KrylovJL_GMRES()),
-                        lazy = !callbacks)
+                    alg = AutoVern7(KenCarp47(autodiff = false, linsolve = LinearSolve.KrylovJL_GMRES()), lazy = !callbacks)
                 elseif length(prob.u0) > 50
-                    alg = AutoVern7(TRBDF2(autodiff = false), lazy = !callbacks)
+                    alg = AutoVern7(KenCarp47(autodiff = false), lazy = !callbacks)
                 else
-                    alg = AutoVern7(Rodas4(autodiff = false), lazy = !callbacks)
+                    alg = AutoVern7(Rodas5P(autodiff = false), lazy = !callbacks)
                 end
             else # :med or low
-                if length(prob.u0) > 100
+                if length(prob.u0) > 500
+                    alg = AutoTsit5(TRBDF2(autodiff = false, linsolve = LinearSolve.KrylovJL_GMRES()))
+                elseif length(prob.u0) > 50
                     alg = AutoTsit5(TRBDF2(autodiff = false))
                 else
                     alg = AutoTsit5(Rosenbrock23(autodiff = false))

test/default_ode_alg_test.jl

@@ -20,16 +20,16 @@ sol = solve(prob_ode_2Dlinear; reltol = 1e-1)
 sol = solve(prob_ode_2Dlinear; reltol = 1e-7)
 
 @test typeof(sol.alg.algs[1]) <: Vern7
-@test typeof(sol.alg.algs[2]) <: Rodas4
+@test typeof(sol.alg.algs[2]) <: Rodas5P
 
 sol = solve(prob_ode_2Dlinear; reltol = 1e-10)
 
-@test typeof(sol.alg.algs[1]) <: Vern9
+@test typeof(sol.alg.algs[1]) <: Vern7
 @test typeof(sol.alg.algs[2]) <: Rodas5P
 
 sol = solve(prob_ode_2Dlinear; alg_hints = [:stiff])
 
-@test typeof(sol.alg) <: Rodas4
+@test typeof(sol.alg) <: Rodas5P
 
 sol = solve(prob_ode_2Dlinear; alg_hints = [:stiff], reltol = 1e-1)
 
@@ -46,18 +46,18 @@ prob_ode_bigfloat2Dlinear = ODEProblem(f, map(BigFloat, rand(4, 2)) .* ones(4, 2
     (0.0, 1.0))
 
 sol = solve(prob_ode_bigfloat2Dlinear; dt = 1 // 2^(4))
-@test typeof(sol.alg.algs[1]) <: Vern9
+@test typeof(sol.alg.algs[1]) <: Vern7
 @test typeof(sol.alg.algs[2]) <: Rodas5P
 
 default_algorithm(prob_ode_bigfloat2Dlinear; alg_hints = [:stiff])
 
 sol = solve(prob_ode_bigfloat2Dlinear; alg_hints = [:stiff])
 
-@test typeof(sol.alg) <: Rodas4
+@test typeof(sol.alg) <: Rodas5P
 
 sol = solve(prob_ode_bigfloat2Dlinear, nothing; alg_hints = [:stiff])
 
-@test typeof(sol.alg) <: Rodas4
+@test typeof(sol.alg) <: Rodas5P
 
 struct FooAlg end
 
@@ -73,21 +73,21 @@ prob = ODEProblem(f, rand(4, 2) .* ones(4, 2) / 2, (0.0, 1.0))
 
 sol = solve(prob; alg_hints = [:stiff])
 
-@test typeof(sol.alg) <: Rodas4
+@test typeof(sol.alg) <: Rodas5P
 
 sol = solve(prob; alg_hints = [:stiff], reltol = 1e-1)
 
 @test typeof(sol.alg) <: Rosenbrock23
 
 sol = solve(prob; alg_hints = [:stiff], callback = CallbackSet())
 
-@test typeof(sol.alg) <: Rodas4
+@test typeof(sol.alg) <: Rodas5P
 
 prob = ODEProblem(f, rand(4, 2) .* ones(4, 2) / 2, (0.0, 1.0))
 
 alg, kwargs = default_algorithm(prob; alg_hints = [:stiff])
 
-@test typeof(alg) <: Rodas4
+@test typeof(alg) <: Rodas5P
 
 m = 1.0
 ω = 1.0