@@ -8,138 +8,163 @@ author: Chris Rackauckas
88using Distributed
99addprocs(2)
1010
11- p1 = Vector{Any}(undef,3)
12- p2 = Vector{Any}(undef,3)
13- p3 = Vector{Any}(undef,3)
11+ p1 = Vector{Any}(undef, 3)
12+ p2 = Vector{Any}(undef, 3)
13+ p3 = Vector{Any}(undef, 3)
1414
1515@everywhere begin
16- using StochasticDiffEq, SDEProblemLibrary, DiffEqNoiseProcess, Plots, ParallelDataTransfer
17- import SDEProblemLibrary: prob_sde_additive,
18- prob_sde_linear, prob_sde_wave
16+ using StochasticDiffEq, SDEProblemLibrary, DiffEqNoiseProcess, Plots,
17+ ParallelDataTransfer
18+ import SDEProblemLibrary: prob_sde_additive,
19+ prob_sde_linear, prob_sde_wave
1920end
2021
2122using StochasticDiffEq, SDEProblemLibrary, DiffEqNoiseProcess, Plots, ParallelDataTransfer
2223import SDEProblemLibrary: prob_sde_additive,
23- prob_sde_linear, prob_sde_wave
24+ prob_sde_linear, prob_sde_wave
2425
25- probs = Matrix{SDEProblem}(undef,3, 3)
26+ probs = Matrix{SDEProblem}(undef, 3, 3)
2627## Problem 1
2728prob = prob_sde_linear
28- probs[1,1] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM1)))
29- probs[1,2] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM2)))
30- probs[1,3] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM3)))
29+ probs[1, 1] = SDEProblem(prob.f, prob.g, prob.u0, prob.tspan, prob.p,
30+ noise = WienerProcess(0.0, 0.0, 0.0, rswm = RSWM(adaptivealg = :RSwM1)))
31+ probs[1, 2] = SDEProblem(prob.f, prob.g, prob.u0, prob.tspan, prob.p,
32+ noise = WienerProcess(0.0, 0.0, 0.0, rswm = RSWM(adaptivealg = :RSwM2)))
33+ probs[1, 3] = SDEProblem(prob.f, prob.g, prob.u0, prob.tspan, prob.p,
34+ noise = WienerProcess(0.0, 0.0, 0.0, rswm = RSWM(adaptivealg = :RSwM3)))
3135## Problem 2
3236prob = prob_sde_wave
33- probs[2,1] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM1)))
34- probs[2,2] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM2)))
35- probs[2,3] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM3)))
37+ probs[2, 1] = SDEProblem(prob.f, prob.g, prob.u0, prob.tspan, prob.p,
38+ noise = WienerProcess(0.0, 0.0, 0.0, rswm = RSWM(adaptivealg = :RSwM1)))
39+ probs[2, 2] = SDEProblem(prob.f, prob.g, prob.u0, prob.tspan, prob.p,
40+ noise = WienerProcess(0.0, 0.0, 0.0, rswm = RSWM(adaptivealg = :RSwM2)))
41+ probs[2, 3] = SDEProblem(prob.f, prob.g, prob.u0, prob.tspan, prob.p,
42+ noise = WienerProcess(0.0, 0.0, 0.0, rswm = RSWM(adaptivealg = :RSwM3)))
3643## Problem 3
3744prob = prob_sde_additive
38- probs[3,1] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM1)))
39- probs[3,2] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM2)))
40- probs[3,3] = SDEProblem(prob.f,prob.g,prob.u0,prob.tspan,prob.p,noise=WienerProcess(0.0,0.0,0.0,rswm=RSWM(adaptivealg=:RSwM3)))
41-
42- fullMeans = Vector{Array}(undef,3)
43- fullMedians = Vector{Array}(undef,3)
44- fullElapsed = Vector{Array}(undef,3)
45- fullTols = Vector{Array}(undef,3)
45+ probs[3, 1] = SDEProblem(prob.f, prob.g, prob.u0, prob.tspan, prob.p,
46+ noise = WienerProcess(0.0, 0.0, 0.0, rswm = RSWM(adaptivealg = :RSwM1)))
47+ probs[3, 2] = SDEProblem(prob.f, prob.g, prob.u0, prob.tspan, prob.p,
48+ noise = WienerProcess(0.0, 0.0, 0.0, rswm = RSWM(adaptivealg = :RSwM2)))
49+ probs[3, 3] = SDEProblem(prob.f, prob.g, prob.u0, prob.tspan, prob.p,
50+ noise = WienerProcess(0.0, 0.0, 0.0, rswm = RSWM(adaptivealg = :RSwM3)))
51+
52+ fullMeans = Vector{Array}(undef, 3)
53+ fullMedians = Vector{Array}(undef, 3)
54+ fullElapsed = Vector{Array}(undef, 3)
55+ fullTols = Vector{Array}(undef, 3)
4656offset = 0
4757
48- Ns = [17,23,
49- 17]
58+ Ns = [17, 23,
59+ 17]
5060```
5161
5262Timings are only valid if no workers die. Workers die if you run out of memory.
5363
5464```julia
55- for k in 1:size(probs,1)
56- global probs, Ns, fullMeans, fullMedians, fullElapsed, fullTols
57- println("Problem $k")
58- ## Setup
59- N = Ns[k]
60-
61- msims = Vector{Any}(undef,N)
62- elapsed = Array{Float64}(undef,N,3)
63- medians = Array{Float64}(undef,N,3)
64- means = Array{Float64}(undef,N,3)
65- tols = Array{Float64}(undef,N,3)
66-
67- #Compile
68- prob = probs[k,1]
69- ParallelDataTransfer.sendto(workers(), prob=prob)
70- monte_prob = EnsembleProblem(prob)
71- solve(monte_prob,SRIW1(),dt=1/2^(4),adaptive=true,trajectories=1000,abstol=2.0^(-1),reltol=0)
72-
73- println("RSwM1")
74- for i=1+offset:N+offset
75- tols[i-offset,1] = 2.0^(-i-1)
76- msims[i-offset] = DiffEqBase.calculate_monte_errors(solve(monte_prob,SRIW1(),
77- trajectories=1000,abstol=2.0^(-i-1),
78- reltol=0,force_dtmin=true))
79- elapsed[i-offset,1] = msims[i-offset].elapsedTime
80- medians[i-offset,1] = msims[i-offset].error_medians[:final]
81- means[i-offset,1] = msims[i-offset].error_means[:final]
82- end
83-
84- println("RSwM2")
85- prob = probs[k,2]
86-
87- ParallelDataTransfer.sendto(workers(), prob=prob)
88- monte_prob = EnsembleProblem(prob)
89- solve(monte_prob,SRIW1(),dt=1/2^(4),adaptive=true,trajectories=1000,abstol=2.0^(-1),reltol=0)
90-
91- for i=1+offset:N+offset
92- tols[i-offset,2] = 2.0^(-i-1)
93- msims[i-offset] = DiffEqBase.calculate_monte_errors(solve(monte_prob,SRIW1(),
94- trajectories=1000,abstol=2.0^(-i-1),
95- reltol=0,force_dtmin=true))
96- elapsed[i-offset,2] = msims[i-offset].elapsedTime
97- medians[i-offset,2] = msims[i-offset].error_medians[:final]
98- means[i-offset,2] = msims[i-offset].error_means[:final]
99- end
100-
101- println("RSwM3")
102- prob = probs[k,3]
103- ParallelDataTransfer.sendto(workers(), prob=prob)
104- monte_prob = EnsembleProblem(prob)
105- solve(monte_prob,SRIW1(),dt=1/2^(4),adaptive=true,trajectories=1000,abstol=2.0^(-1),reltol=0)
106-
107- for i=1+offset:N+offset
108- tols[i-offset,3] = 2.0^(-i-1)
109- msims[i-offset] = DiffEqBase.calculate_monte_errors(solve(monte_prob,SRIW1(),
110- adaptive=true,trajectories=1000,abstol=2.0^(-i-1),
111- reltol=0,force_dtmin=true))
112- elapsed[i-offset,3] = msims[i-offset].elapsedTime
113- medians[i-offset,3] = msims[i-offset].error_medians[:final]
114- means[i-offset,3] = msims[i-offset].error_means[:final]
115- end
116-
117- fullMeans[k] = means
118- fullMedians[k] =medians
119- fullElapsed[k] = elapsed
120- fullTols[k] = tols
65+ for k in 1:size(probs, 1)
66+ global probs, Ns, fullMeans, fullMedians, fullElapsed, fullTols
67+ println("Problem $k")
68+ ## Setup
69+ N = Ns[k]
70+
71+ msims = Vector{Any}(undef, N)
72+ elapsed = Array{Float64}(undef, N, 3)
73+ medians = Array{Float64}(undef, N, 3)
74+ means = Array{Float64}(undef, N, 3)
75+ tols = Array{Float64}(undef, N, 3)
76+
77+ #Compile
78+ prob = probs[k, 1]
79+ ParallelDataTransfer.sendto(workers(), prob = prob)
80+ monte_prob = EnsembleProblem(prob)
81+ solve(monte_prob, SRIW1(), dt = 1/2^(4), adaptive = true,
82+ trajectories = 1000, abstol = 2.0^(-1), reltol = 0)
83+
84+ println("RSwM1")
85+ for i in (1 + offset):(N + offset)
86+ tols[i - offset, 1] = 2.0^(-i-1)
87+ msims[i - offset] = DiffEqBase.calculate_monte_errors(solve(monte_prob, SRIW1(),
88+ trajectories = 1000, abstol = 2.0^(-i-1),
89+ reltol = 0, force_dtmin = true))
90+ elapsed[i - offset, 1] = msims[i - offset].elapsedTime
91+ medians[i - offset, 1] = msims[i - offset].error_medians[:final]
92+ means[i - offset, 1] = msims[i - offset].error_means[:final]
93+ end
94+
95+ println("RSwM2")
96+ prob = probs[k, 2]
97+
98+ ParallelDataTransfer.sendto(workers(), prob = prob)
99+ monte_prob = EnsembleProblem(prob)
100+ solve(monte_prob, SRIW1(), dt = 1/2^(4), adaptive = true,
101+ trajectories = 1000, abstol = 2.0^(-1), reltol = 0)
102+
103+ for i in (1 + offset):(N + offset)
104+ tols[i - offset, 2] = 2.0^(-i-1)
105+ msims[i - offset] = DiffEqBase.calculate_monte_errors(solve(monte_prob, SRIW1(),
106+ trajectories = 1000, abstol = 2.0^(-i-1),
107+ reltol = 0, force_dtmin = true))
108+ elapsed[i - offset, 2] = msims[i - offset].elapsedTime
109+ medians[i - offset, 2] = msims[i - offset].error_medians[:final]
110+ means[i - offset, 2] = msims[i - offset].error_means[:final]
111+ end
112+
113+ println("RSwM3")
114+ prob = probs[k, 3]
115+ ParallelDataTransfer.sendto(workers(), prob = prob)
116+ monte_prob = EnsembleProblem(prob)
117+ solve(monte_prob, SRIW1(), dt = 1/2^(4), adaptive = true,
118+ trajectories = 1000, abstol = 2.0^(-1), reltol = 0)
119+
120+ for i in (1 + offset):(N + offset)
121+ tols[i - offset, 3] = 2.0^(-i-1)
122+ msims[i - offset] = DiffEqBase.calculate_monte_errors(solve(monte_prob, SRIW1(),
123+ adaptive = true, trajectories = 1000, abstol = 2.0^(-i-1),
124+ reltol = 0, force_dtmin = true))
125+ elapsed[i - offset, 3] = msims[i - offset].elapsedTime
126+ medians[i - offset, 3] = msims[i - offset].error_medians[:final]
127+ means[i - offset, 3] = msims[i - offset].error_means[:final]
128+ end
129+
130+ fullMeans[k] = means
131+ fullMedians[k] = medians
132+ fullElapsed[k] = elapsed
133+ fullTols[k] = tols
121134end
122135```
123136
124137```julia
125- gr(fmt= :svg)
138+ gr(fmt = :svg)
126139lw=3
127- leg=String["RSwM1","RSwM2","RSwM3"]
140+ leg=String["RSwM1", "RSwM2", "RSwM3"]
128141
129142titleFontSize = 16
130143guideFontSize = 14
131- legendFontSize= 14
132- tickFontSize = 12
133-
134- for k in 1:size(probs,1)
135- global probs, Ns, fullMeans, fullMedians, fullElapsed, fullTols
136- p1[k] = Plots.plot(fullTols[k],fullMeans[k],xscale=:log10,yscale=:log10, xguide="Absolute Tolerance",yguide="Mean Final Error",title="Example $k" ,linewidth=lw,grid=false,lab=leg,titlefont=font(titleFontSize),legendfont=font(legendFontSize),tickfont=font(tickFontSize),guidefont=font(guideFontSize))
137- p2[k] = Plots.plot(fullTols[k],fullMedians[k],xscale=:log10,yscale=:log10,xguide="Absolute Tolerance",yguide="Median Final Error",title="Example $k",linewidth=lw,grid=false,lab=leg,titlefont=font(titleFontSize),legendfont=font(legendFontSize),tickfont=font(tickFontSize),guidefont=font(guideFontSize))
138- p3[k] = Plots.plot(fullTols[k],fullElapsed[k],xscale=:log10,yscale=:log10,xguide="Absolute Tolerance",yguide="Elapsed Time",title="Example $k" ,linewidth=lw,grid=false,lab=leg,titlefont=font(titleFontSize),legendfont=font(legendFontSize),tickfont=font(tickFontSize),guidefont=font(guideFontSize))
144+ legendFontSize = 14
145+ tickFontSize = 12
146+
147+ for k in 1:size(probs, 1)
148+ global probs, Ns, fullMeans, fullMedians, fullElapsed, fullTols
149+ p1[k] = Plots.plot(fullTols[k], fullMeans[k], xscale = :log10, yscale = :log10,
150+ xguide = "Absolute Tolerance", yguide = "Mean Final Error",
151+ title = "Example $k", linewidth = lw, grid = false, lab = leg,
152+ titlefont = font(titleFontSize), legendfont = font(legendFontSize),
153+ tickfont = font(tickFontSize), guidefont = font(guideFontSize))
154+ p2[k] = Plots.plot(fullTols[k], fullMedians[k], xscale = :log10, yscale = :log10,
155+ xguide = "Absolute Tolerance", yguide = "Median Final Error",
156+ title = "Example $k", linewidth = lw, grid = false, lab = leg,
157+ titlefont = font(titleFontSize), legendfont = font(legendFontSize),
158+ tickfont = font(tickFontSize), guidefont = font(guideFontSize))
159+ p3[k] = Plots.plot(fullTols[k], fullElapsed[k], xscale = :log10, yscale = :log10,
160+ xguide = "Absolute Tolerance", yguide = "Elapsed Time",
161+ title = "Example $k", linewidth = lw, grid = false, lab = leg,
162+ titlefont = font(titleFontSize), legendfont = font(legendFontSize),
163+ tickfont = font(tickFontSize), guidefont = font(guideFontSize))
139164end
140165
141166Plots.plot!(p1[1])
142- Plots.plot(p1[1],p1[2],p1[3],layout= (3,1),size= (1000,800))
167+ Plots.plot(p1[1], p1[2], p1[3], layout = (3, 1), size = (1000, 800))
143168```
144169
145170```julia
@@ -148,17 +173,17 @@ Plots.plot(p1[1],p1[2],p1[3],layout=(3,1),size=(1000,800))
148173```
149174
150175```julia
151- plot(p3[1],p3[2],p3[3],layout= (3,1),size= (1000,800))
176+ plot(p3[1], p3[2], p3[3], layout = (3, 1), size = (1000, 800))
152177#savefig("timevstol.png")
153178#savefig("timevstol.pdf")
154179```
155180
156181```julia
157- plot(p1[1],p3[1],p1[2],p3[2],p1[3],p3[3],layout= (3,2),size= (1000,800))
182+ plot(p1[1], p3[1], p1[2], p3[2], p1[3], p3[3], layout = (3, 2), size = (1000, 800))
158183```
159184
160185```julia
161186
162187using SciMLBenchmarks
163- SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])
188+ SciMLBenchmarks.bench_footer(WEAVE_ARGS[:folder], WEAVE_ARGS[:file])
164189```
0 commit comments