Merge pull request #116 from ErikQQY/qqy/bvp_working_precision

Add BVP precision diagram

Author: ChrisRackauckas

Project.toml

@@ -1,7 +1,7 @@
 name = "DiffEqDevTools"
 uuid = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "2.36.0"
+version = "2.37.0"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
@@ -28,6 +28,8 @@ SciMLBase = "1.74"
 julia = "1.6"
 
 [extras]
+BoundaryValueDiffEq = "764a87c0-6b3e-53db-9096-fe964310641d"
+BVProblemLibrary = "ded0fc24-dfea-4565-b1d9-79c027d14d84"
 DelayDiffEq = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb"
 DDEProblemLibrary = "f42792ee-6ffc-4e2a-ae83-8ee2f22de800"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
@@ -42,4 +44,4 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["DelayDiffEq", "DDEProblemLibrary", "NonlinearSolve", "ODEProblemLibrary", "SDEProblemLibrary", "OrdinaryDiffEq", "ParameterizedFunctions", "Random", "StochasticDiffEq", "StochasticDelayDiffEq", "Plots", "Test"]
+test = ["BoundaryValueDiffEq", "BVProblemLibrary", "DelayDiffEq", "DDEProblemLibrary", "NonlinearSolve", "ODEProblemLibrary", "SDEProblemLibrary", "OrdinaryDiffEq", "ParameterizedFunctions", "Random", "StochasticDiffEq", "StochasticDelayDiffEq", "Plots", "Test"]

src/DiffEqDevTools.jl

@@ -11,7 +11,7 @@ import Base: length
 import DiffEqBase: AbstractODEProblem, AbstractDDEProblem, AbstractDDEAlgorithm,
                    AbstractODESolution, AbstractRODEProblem, AbstractSDEProblem,
                    AbstractSDDEProblem, AbstractEnsembleProblem,
-                   AbstractDAEProblem, @def, ConvergenceSetup, DEAlgorithm,
+                   AbstractDAEProblem, AbstractBVProblem, @def, ConvergenceSetup, DEAlgorithm,
                    ODERKTableau, AbstractTimeseriesSolution, ExplicitRKTableau,
                    ImplicitRKTableau
 

src/benchmark.jl

@@ -273,6 +273,103 @@ function WorkPrecision(prob, alg, abstols, reltols, dts = nothing;
     return WorkPrecision(prob, abstols, reltols, errors, times, name, N)
 end
 
+# Work precision information for a BVP
+function WorkPrecision(prob::AbstractBVProblem, alg, abstols, reltols, dts = nothing;
+    name = nothing, appxsol = nothing, error_estimate = :final,
+    numruns = 20, seconds = 2, kwargs...)
+    N = length(abstols)
+    errors = Vector{Float64}(undef, N)
+    times = Vector{Float64}(undef, N)
+    if name === nothing
+        name = "WP-Alg"
+    end
+
+    if haskey(kwargs, :prob_choice)
+        _prob = prob[kwargs[:prob_choice]]
+    elseif prob isa AbstractArray
+        _prob = prob[1]
+    else
+        _prob = prob
+    end
+
+    let _prob = _prob
+        timeseries_errors = error_estimate ∈ TIMESERIES_ERRORS
+        dense_errors = error_estimate ∈ DENSE_ERRORS
+        for i in 1:N
+            if dts === nothing
+                sol = solve(_prob, alg; kwargs..., abstol = abstols[i],
+                    reltol = reltols[i], timeseries_errors = timeseries_errors,
+                    dense_errors = dense_errors)
+            else
+                sol = DiffEqBase.solve(_prob, alg; kwargs..., abstol = abstols[i],
+                    reltol = reltols[i], dt = dts[i],
+                    timeseries_errors = timeseries_errors,
+                    dense_errors = dense_errors)
+            end
+
+            if haskey(kwargs, :prob_choice)
+                cur_appxsol = appxsol[kwargs[:prob_choice]]
+            elseif prob isa AbstractArray
+                cur_appxsol = appxsol[1]
+            else
+                cur_appxsol = appxsol
+            end
+
+            if cur_appxsol !== nothing
+                errsol = appxtrue(sol, cur_appxsol)
+                errors[i] = mean(errsol.errors[error_estimate])
+            else
+                errors[i] = mean(sol.errors[error_estimate])
+            end
+
+            benchmark_f = let dts = dts, _prob = _prob, alg = alg, sol = sol,
+                abstols = abstols, reltols = reltols, kwargs = kwargs
+
+                if dts === nothing
+                    if typeof(_prob) <: DAEProblem
+                        () -> @elapsed solve(_prob, alg;
+                            abstol = abstols[i],
+                            reltol = reltols[i],
+                            timeseries_errors = false,
+                            dense_errors = false, kwargs...)
+                    else
+                        () -> @elapsed solve(_prob, alg;
+                            abstol = abstols[i],
+                            reltol = reltols[i],
+                            timeseries_errors = false,
+                            dense_errors = false, kwargs...)
+                    end
+                else
+                    if typeof(_prob) <: DAEProblem
+                        () -> @elapsed solve(_prob, alg;
+                            abstol = abstols[i],
+                            reltol = reltols[i],
+                            dt = dts[i],
+                            timeseries_errors = false,
+                            dense_errors = false, kwargs...)
+                    else
+                        () -> @elapsed solve(_prob, alg;
+                            abstol = abstols[i],
+                            reltol = reltols[i],
+                            dt = dts[i],
+                            timeseries_errors = false,
+                            dense_errors = false, kwargs...)
+                    end
+                end
+            end
+            benchmark_f() # pre-compile
+
+            b_t = benchmark_f()
+            if b_t > seconds
+                times[i] = b_t
+            else
+                times[i] = mapreduce(i -> benchmark_f(), min, 2:numruns; init = b_t)
+            end
+        end
+    end
+    return WorkPrecision(prob, abstols, reltols, errors, times, name, N)
+end
+
 # Work precision information for a nonlinear problem.
 function WorkPrecision(prob::NonlinearProblem, alg, abstols, reltols, dts = nothing; name = nothing, appxsol = nothing, error_estimate = :l2, numruns = 20, seconds = 2, kwargs...)
     isnothing(appxsol) && error("Must provide the real value as the \"appxsol\" kwarg.")
@@ -587,6 +684,33 @@ function WorkPrecisionSet(prob::AbstractEnsembleProblem, abstols, reltols, setup
                      Int(trajectories))
 end
 
+function WorkPrecisionSet(prob::AbstractBVProblem,
+    abstols, reltols, setups;
+    print_names = false, names = nothing, appxsol = nothing,
+    error_estimate = :final,
+    test_dt = nothing, kwargs...)
+    N = length(setups)
+    @assert names === nothing || length(setups) == length(names)
+    wps = Vector{WorkPrecision}(undef, N)
+    if names === nothing
+        names = [_default_name(setup[:alg]) for setup in setups]
+    end
+    for i in 1:N
+        print_names && println(names[i])
+        _abstols = get(setups[i], :abstols, abstols)
+        _reltols = get(setups[i], :reltols, reltols)
+        _dts = get(setups[i], :dts, nothing)
+        filtered_setup = filter(p -> p.first in DiffEqBase.allowedkeywords, setups[i])
+
+        wps[i] = WorkPrecision(prob, setups[i][:alg], _abstols, _reltols, _dts;
+            appxsol = appxsol,
+            error_estimate = error_estimate,
+            name = names[i], kwargs..., filtered_setup...)
+    end
+    return WorkPrecisionSet(wps, N, abstols, reltols, prob, setups, names, error_estimate,
+        nothing)
+end
+
 function get_sample_errors(prob::AbstractRODEProblem, setup, test_dt = nothing;
                            appxsol_setup = nothing,
                            numruns, error_estimate = :final,

test/benchmark_tests.jl

@@ -1,6 +1,7 @@
-using OrdinaryDiffEq, DelayDiffEq, DiffEqDevTools, DiffEqBase, Test
+using OrdinaryDiffEq, DelayDiffEq, BoundaryValueDiffEq, DiffEqDevTools, DiffEqBase, Test
 using ODEProblemLibrary: prob_ode_2Dlinear, prob_ode_linear
 using DDEProblemLibrary: prob_dde_constant_1delay_ip
+using BVProblemLibrary: prob_bvp_linear_1
 
 ## Setup Tests
 
@@ -176,3 +177,25 @@ setups = [Dict(:alg => RadauIIA5()),
     Dict(:alg => RosShamp4())]
 shoot = Shootout(prob, setups; appxsol = TestSolution(sol))
 @test shoot.names == ["RadauIIA5", "RosShamp4"]
+
+# BVP problem
+prob = prob_bvp_linear_1
+abstols = 1.0 ./ 10.0 .^ (2:3)
+reltols = 1.0 ./ 10.0 .^ (2:3)
+sol = solve(prob, Shooting(Tsit5()), abstol = 1e-14, reltol = 1e-14)
+test_sol = TestSolution(sol.t, sol.u)
+
+setups = [Dict(:alg => MIRK4(), :dts => 1.0 ./ 5.0 .^ ((1:length(reltols)) .+ 1))
+    Dict(:alg => MIRK5(), :dts => 1.0 ./ 5.0 .^ ((1:length(reltols)) .+ 1))]
+labels = ["MIRK4", "MIRK5"]
+
+println("Test MIRK4 and MIRK5")
+wp = WorkPrecisionSet(prob,
+    abstols,
+    reltols,
+    setups;
+    approxsol = sol,
+    names = labels,
+    print_names = true)
+@test wp.names == ["MIRK4", "MIRK5"]
+println("BVP Done")