Add bench script

Author: tansongchen

bench/Project.toml

@@ -0,0 +1,8 @@
+[deps]
+DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqTaylorSeries = "9c7f1690-dd92-42a3-8318-297ee24d8d39"
+ParameterizedFunctions = "65888b18-ceab-5e60-b2b9-181511a3b968"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+TaylorIntegration = "92b13dbe-c966-51a2-8445-caca9f8a7d42"

bench/adaptive_order.jl

@@ -0,0 +1,40 @@
+using OrdinaryDiffEq, ParameterizedFunctions, DiffEqDevTools, StaticArrays,
+      OrdinaryDiffEqTaylorSeries
+using Plots, BenchmarkTools
+gr();
+
+abstols = 1.0 ./ 10.0 .^ (6:15)
+reltols = 1.0 ./ 10.0 .^ (3:12);
+
+function lotka(du, u, p, t)
+    x = u[1]
+    y = u[2]
+    du[1] = p[1] * x - p[2] * x * y
+    du[2] = -p[3] * y + p[4] * x * y
+end
+prob = ODEProblem{true, SciMLBase.FullSpecialize}(
+    lotka, [1.0, 1.0], (0.0, 10.0), [1.5, 1.0, 3.0, 1.0])
+sol = solve(prob, Vern7(), abstol = 1 / 10^14, reltol = 1 / 10^14)
+
+# Profile
+cache1 = init(prob, ExplicitTaylorAdaptiveOrder(min_order = Val(6), max_order = Val(12)))
+OrdinaryDiffEqTaylorSeries.perform_step!(cache1, cache1.cache, false)
+@profview for _ in 1:10000000
+    OrdinaryDiffEqTaylorSeries.perform_step!(cache1, cache1.cache, false)
+end
+@btime OrdinaryDiffEqTaylorSeries.perform_step!(cache1, cache1.cache, false)
+
+# Make work-precision plot
+setups = [Dict(:alg => DP5())
+          Dict(:alg => Tsit5())
+          Dict(:alg => Vern6())
+          Dict(:alg => Vern8())
+          Dict(:alg => ExplicitTaylorAdaptiveOrder(min_order = Val(8), max_order = Val(12)))];
+names = ["DP5", "Tsit5", "Vern6", "Vern8", "Taylor"]
+for order in 6:2:12
+    push!(names, "Taylor $(order)")
+    push!(setups, Dict(:alg => ExplicitTaylor(order = Val(order + 1))))
+end
+wp = WorkPrecisionSet([prob], abstols, reltols, setups; names = names, appxsol = [sol],
+    save_everystep = false, numruns = 100, maxiters = 10000)
+plot(wp)

lib/OrdinaryDiffEqTaylorSeries/src/TaylorSeries_perform_step.jl

@@ -110,6 +110,7 @@ end
 
     min_order_value = get_value(min_order)
     max_order_value = get_value(max_order)
+    # println("current order: ", current_order[])
     jet_index = current_order[] - min_order_value + 1
     # compute one additional order for adaptive order
     jet = jets[jet_index + 1]
@@ -122,10 +123,12 @@ end
         min_work = Inf
         start_order = max(min_order_value, current_order[] - 1)
         end_order = min(max_order_value - 1, current_order[] + 1)
+        dtn = dt^start_order
         for i in start_order:end_order
             A = i * i
-            @.. broadcast=false thread=thread utilde=TaylorDiff.get_coefficient(
-                utaylor, i) * dt^i
+            for ix in eachindex(utaylor)
+                utilde[ix] = TaylorDiff.partials(utaylor[ix])[i] * dtn
+            end
             calculate_residuals!(atmp, utilde, uprev, u, integrator.opts.abstol,
                 integrator.opts.reltol, integrator.opts.internalnorm, t)
             EEst = integrator.opts.internalnorm(atmp, t)
@@ -147,6 +150,7 @@ end
                 min_work = work
                 integrator.EEst = EEst
             end
+            dtn *= dt
         end
     end
     return nothing