From f872e5dc1993703446f1fcc38f7ee714e9e400a0 Mon Sep 17 00:00:00 2001 From: Nathanael Bosch Date: Thu, 12 Oct 2023 13:19:00 +0200 Subject: [PATCH] Add Pleiades as a second-order ODE benchmark --- benchmarks/pleiades.jmd | 131 ++++++++ benchmarks/runall.jl | 1 + docs/make.jl | 1 + docs/src/benchmarks/figures/pleiades_2_1.svg | 156 +++++++++ docs/src/benchmarks/figures/pleiades_3_1.svg | 314 +++++++++++++++++++ docs/src/benchmarks/figures/pleiades_4_1.svg | 156 +++++++++ docs/src/benchmarks/figures/pleiades_5_1.svg | 231 ++++++++++++++ docs/src/benchmarks/pleiades.md | 161 ++++++++++ 8 files changed, 1151 insertions(+) create mode 100644 benchmarks/pleiades.jmd create mode 100644 docs/src/benchmarks/figures/pleiades_2_1.svg create mode 100644 docs/src/benchmarks/figures/pleiades_3_1.svg create mode 100644 docs/src/benchmarks/figures/pleiades_4_1.svg create mode 100644 docs/src/benchmarks/figures/pleiades_5_1.svg create mode 100644 docs/src/benchmarks/pleiades.md diff --git a/benchmarks/pleiades.jmd b/benchmarks/pleiades.jmd new file mode 100644 index 000000000..004999a68 --- /dev/null +++ b/benchmarks/pleiades.jmd @@ -0,0 +1,131 @@ +# Pleiades benchmark + +```julia +using LinearAlgebra, Statistics, InteractiveUtils +using DiffEqDevTools, ParameterizedFunctions, SciMLBase, OrdinaryDiffEq, Sundials, Plots +using ModelingToolkit +using ProbNumDiffEq + +# Plotting theme +theme(:dao; + markerstrokewidth=0.5, + legend=:outertopright, + bottom_margin=5Plots.mm, + size = (1000, 400), +) +``` + +### Pleiades problem definition + +```julia +# first-order ODE +@fastmath function pleiades(du, u, p, t) + v = view(u, 1:7) # x + w = view(u, 8:14) # y + x = view(u, 15:21) # x′ + y = view(u, 22:28) # y′ + du[15:21] .= v + du[22:28] .= w + @inbounds @simd ivdep for i in 1:14 + du[i] = zero(eltype(u)) + end + @inbounds @simd ivdep for i in 1:7 + @inbounds @simd ivdep for j in 1:7 + if i != j + r = ((x[i] - x[j])^2 + (y[i] - y[j])^2)^(3 / 2) + du[i] += j * (x[j] - x[i]) / r + du[7+i] += j * (y[j] - y[i]) / r + end + end + end +end +x0 = [3.0, 3.0, -1.0, -3.0, 2.0, -2.0, 2.0] +y0 = [3.0, -3.0, 2.0, 0, 0, -4.0, 4.0] +dx0 = [0, 0, 0, 0, 0, 1.75, -1.5] +dy0 = [0, 0, 0, -1.25, 1, 0, 0] +u0 = [dx0; dy0; x0; y0] +tspan = (0.0, 3.0) +prob1 = ODEProblem(pleiades, u0, tspan) + +# second-order ODE +function pleiades2(ddu, du, u, p, t) + x = view(u, 1:7) + y = view(u, 8:14) + for i in 1:14 + ddu[i] = zero(eltype(u)) + end + for i in 1:7, j in 1:7 + if i != j + r = ((x[i] - x[j])^2 + (y[i] - y[j])^2)^(3 / 2) + ddu[i] += j * (x[j] - x[i]) / r + ddu[7+i] += j * (y[j] - y[i]) / r + end + end +end +u0 = [x0; y0] +du0 = [dx0; dy0] +prob2 = SecondOrderODEProblem(pleiades2, du0, u0, tspan) +probs = [prob1, prob2] + +ref_sol1 = solve(prob1, Vern9(), abstol=1/10^14, reltol=1/10^14, dense=false) +ref_sol2 = solve(prob2, Vern9(), abstol=1/10^14, reltol=1/10^14, dense=false) +ref_sols = [ref_sol1, ref_sol2] + +plot(ref_sol1, idxs=[(14+i,21+i) for i in 1:7], title="Pleiades Solution", legend=false) +scatter!(ref_sol1.u[end][15:21], ref_sol1.u[end][22:end], color=1:7) +``` + +## First-order ODE vs. second-order ODE +```julia +DENSE = false; +SAVE_EVERYSTEP = false; + +_setups = [ + "EK0(3) (1st order ODE)" => Dict(:alg => EK0(order=3, smooth=DENSE), :prob_choice => 1) + "EK0(5) (1st order ODE)" => Dict(:alg => EK0(order=5, smooth=DENSE), :prob_choice => 1) + "EK1(3) (1st order ODE)" => Dict(:alg => EK1(order=3, smooth=DENSE), :prob_choice => 1) + "EK1(5) (1st order ODE)" => Dict(:alg => EK1(order=5, smooth=DENSE), :prob_choice => 1) + "EK0(4) (2nd order ODE)" => Dict(:alg => EK0(order=4, smooth=DENSE), :prob_choice => 2) + "EK0(6) (2nd order ODE)" => Dict(:alg => EK0(order=6, smooth=DENSE), :prob_choice => 2) + "EK1(4) (2nd order ODE)" => Dict(:alg => EK1(order=4, smooth=DENSE), :prob_choice => 2) + "EK1(6) (2nd order ODE)" => Dict(:alg => EK1(order=6, smooth=DENSE), :prob_choice => 2) +] + +labels = first.(_setups) +setups = last.(_setups) + +abstols = 1.0 ./ 10.0 .^ (6:11) +reltols = 1.0 ./ 10.0 .^ (3:8) + +wp = WorkPrecisionSet( + probs, abstols, reltols, setups; + names = labels, + #print_names = true, + appxsol = ref_sols, + dense = DENSE, + save_everystep = SAVE_EVERYSTEP, + numruns = 10, + maxiters = Int(1e7), + timeseries_errors = false, + verbose = false, +) + +plot(wp, color=[1 1 2 2 3 3 4 4], + # xticks = 10.0 .^ (-16:1:5) +) +``` + + +## Conclusion + +- If the problem is a second-order ODE, _implement it as a second-order ODE_! + Just use `SecondOrderODEProblem`. + + +## Appendix + +Computer information: + +```julia +InteractiveUtils.versioninfo() +``` diff --git a/benchmarks/runall.jl b/benchmarks/runall.jl index 5790315bb..d086e17d0 100644 --- a/benchmarks/runall.jl +++ b/benchmarks/runall.jl @@ -4,6 +4,7 @@ FILES = [ "lotkavolterra.jmd", "vanderpol.jmd", "rober.jmd", + "pleiades.jmd", "multi-language-wrappers.jmd", ] diff --git a/docs/make.jl b/docs/make.jl index 41131646f..4c3c48fee 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -50,6 +50,7 @@ makedocs( "Multi-Language Wrapper Benchmark" => "benchmarks/multi-language-wrappers.md", "Non-stiff ODE: Lotka-Volterra" => "benchmarks/lotkavolterra.md", "Stiff ODE: Van der Pol" => "benchmarks/vanderpol.md", + "Second order ODE: Pleiades" => "benchmarks/pleiades.md", "DAE: ROBER" => "benchmarks/rober.md", ], "Internals" => [ diff --git a/docs/src/benchmarks/figures/pleiades_2_1.svg b/docs/src/benchmarks/figures/pleiades_2_1.svg new file mode 100644 index 000000000..81f8b19e6 --- /dev/null +++ b/docs/src/benchmarks/figures/pleiades_2_1.svg @@ -0,0 +1,156 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/benchmarks/figures/pleiades_3_1.svg b/docs/src/benchmarks/figures/pleiades_3_1.svg new file mode 100644 index 000000000..cdbddeec1 --- /dev/null +++ b/docs/src/benchmarks/figures/pleiades_3_1.svg @@ -0,0 +1,314 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/benchmarks/figures/pleiades_4_1.svg b/docs/src/benchmarks/figures/pleiades_4_1.svg new file mode 100644 index 000000000..0779f282a --- /dev/null +++ b/docs/src/benchmarks/figures/pleiades_4_1.svg @@ -0,0 +1,156 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/benchmarks/figures/pleiades_5_1.svg b/docs/src/benchmarks/figures/pleiades_5_1.svg new file mode 100644 index 000000000..e26844737 --- /dev/null +++ b/docs/src/benchmarks/figures/pleiades_5_1.svg @@ -0,0 +1,231 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/benchmarks/pleiades.md b/docs/src/benchmarks/pleiades.md new file mode 100644 index 000000000..2c8cae29c --- /dev/null +++ b/docs/src/benchmarks/pleiades.md @@ -0,0 +1,161 @@ +# Pleiades benchmark + +```julia +using LinearAlgebra, Statistics, InteractiveUtils +using DiffEqDevTools, ParameterizedFunctions, SciMLBase, OrdinaryDiffEq, Sundials, Plots +using ModelingToolkit +using ProbNumDiffEq + +# Plotting theme +theme(:dao; + markerstrokewidth=0.5, + legend=:outertopright, + bottom_margin=5Plots.mm, + size = (1000, 400), +) +``` + + + + +### Pleiades problem definition + +```julia +# first-order ODE +@fastmath function pleiades(du, u, p, t) + v = view(u, 1:7) # x + w = view(u, 8:14) # y + x = view(u, 15:21) # x′ + y = view(u, 22:28) # y′ + du[15:21] .= v + du[22:28] .= w + @inbounds @simd ivdep for i in 1:14 + du[i] = zero(eltype(u)) + end + @inbounds @simd ivdep for i in 1:7 + @inbounds @simd ivdep for j in 1:7 + if i != j + r = ((x[i] - x[j])^2 + (y[i] - y[j])^2)^(3 / 2) + du[i] += j * (x[j] - x[i]) / r + du[7+i] += j * (y[j] - y[i]) / r + end + end + end +end +x0 = [3.0, 3.0, -1.0, -3.0, 2.0, -2.0, 2.0] +y0 = [3.0, -3.0, 2.0, 0, 0, -4.0, 4.0] +dx0 = [0, 0, 0, 0, 0, 1.75, -1.5] +dy0 = [0, 0, 0, -1.25, 1, 0, 0] +u0 = [dx0; dy0; x0; y0] +tspan = (0.0, 3.0) +prob1 = ODEProblem(pleiades, u0, tspan) + +# second-order ODE +function pleiades2(ddu, du, u, p, t) + x = view(u, 1:7) + y = view(u, 8:14) + for i in 1:14 + ddu[i] = zero(eltype(u)) + end + for i in 1:7, j in 1:7 + if i != j + r = ((x[i] - x[j])^2 + (y[i] - y[j])^2)^(3 / 2) + ddu[i] += j * (x[j] - x[i]) / r + ddu[7+i] += j * (y[j] - y[i]) / r + end + end +end +u0 = [x0; y0] +du0 = [dx0; dy0] +prob2 = SecondOrderODEProblem(pleiades2, du0, u0, tspan) +probs = [prob1, prob2] + +ref_sol1 = solve(prob1, Vern9(), abstol=1/10^14, reltol=1/10^14, dense=false) +ref_sol2 = solve(prob2, Vern9(), abstol=1/10^14, reltol=1/10^14, dense=false) +ref_sols = [ref_sol1, ref_sol2] + +plot(ref_sol1, idxs=[(14+i,21+i) for i in 1:7], title="Pleiades Solution", legend=false) +scatter!(ref_sol1.u[end][15:21], ref_sol1.u[end][22:end], color=1:7) +``` + +![](figures/pleiades_2_1.svg) + + + +## First-order ODE vs. second-order ODE +```julia +DENSE = false; +SAVE_EVERYSTEP = false; + +_setups = [ + "EK0(3) (1st order ODE)" => Dict(:alg => EK0(order=3, smooth=DENSE), :prob_choice => 1) + "EK0(5) (1st order ODE)" => Dict(:alg => EK0(order=5, smooth=DENSE), :prob_choice => 1) + "EK1(3) (1st order ODE)" => Dict(:alg => EK1(order=3, smooth=DENSE), :prob_choice => 1) + "EK1(5) (1st order ODE)" => Dict(:alg => EK1(order=5, smooth=DENSE), :prob_choice => 1) + "EK0(4) (2nd order ODE)" => Dict(:alg => EK0(order=4, smooth=DENSE), :prob_choice => 2) + "EK0(6) (2nd order ODE)" => Dict(:alg => EK0(order=6, smooth=DENSE), :prob_choice => 2) + "EK1(4) (2nd order ODE)" => Dict(:alg => EK1(order=4, smooth=DENSE), :prob_choice => 2) + "EK1(6) (2nd order ODE)" => Dict(:alg => EK1(order=6, smooth=DENSE), :prob_choice => 2) +] + +labels = first.(_setups) +setups = last.(_setups) + +abstols = 1.0 ./ 10.0 .^ (6:11) +reltols = 1.0 ./ 10.0 .^ (3:8) + +wp = WorkPrecisionSet( + probs, abstols, reltols, setups; + names = labels, + #print_names = true, + appxsol = ref_sols, + dense = DENSE, + save_everystep = SAVE_EVERYSTEP, + numruns = 10, + maxiters = Int(1e7), + timeseries_errors = false, + verbose = false, +) + +plot(wp, color=[1 1 2 2 3 3 4 4], + # xticks = 10.0 .^ (-16:1:5) +) +``` + +![](figures/pleiades_3_1.svg) + + + + +## Conclusion + +- If the problem is a second-order ODE, _implement it as a second-order ODE_! + Just use `SecondOrderODEProblem`. + + +## Appendix + +Computer information: + +```julia +InteractiveUtils.versioninfo() +``` + +``` +Julia Version 1.9.3 +Commit bed2cd540a1 (2023-08-24 14:43 UTC) +Build Info: + Official https://julialang.org/ release +Platform Info: + OS: Linux (x86_64-linux-gnu) + CPU: 12 × Intel(R) Core(TM) i7-6800K CPU @ 3.40GHz + WORD_SIZE: 64 + LIBM: libopenlibm + LLVM: libLLVM-14.0.6 (ORCJIT, broadwell) + Threads: 12 on 12 virtual cores +Environment: + JULIA_NUM_THREADS = auto + JULIA_STACKTRACE_MINIMAL = true +``` + +