Mercurial > repos > public > sbplib_julia

comparison benchmark/benchmark_laplace.jl @ 1376:2ad2de55061a feature/variable_derivatives

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add temp. benchmark script for constant and variable coeff Laplace
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 30 May 2023 17:27:45 +0200
parents
children 43aaf710463e
comparison
equal deleted inserted replaced
1375:127d2558926e 1376:2ad2de55061a
1 using Sbplib
2 using Sbplib.SbpOperators
3 using Sbplib.Grids
4 using Sbplib.RegionIndices
5 using BenchmarkTools
6
7 # TODO: Move the below benchmarks into the benchmark suite
8
9 const operator_path = sbp_operators_path()*"standard_diagonal.toml"
10
11 function benchmark_const_coeff_1d(;N = 100, order = 4)
12 stencil_set = read_stencil_set(operator_path; order=order)
13 g = equidistant_grid(N, 0., 1.)
14 D = second_derivative(g, stencil_set)
15 u = rand(size(g)...)
16 u_xx = rand(size(g)...)
17
18 b_naive = @benchmark $u_xx .= $D*$u
19 b_reg = @benchmark $apply_region_1d!($u_xx,$u,$D)
20 b_thrd = @benchmark $apply_region_threaded_1d!($u_xx,$u,$D)
21 print_benchmark_result("Naive apply",b_naive)
22 print_benchmark_result("Region apply",b_reg)
23 print_benchmark_result("Threaded region apply",b_thrd)
24 end
25
26 function benchmark_var_coeff_1d(;N = 100, order = 4)
27 stencil_set = read_stencil_set(operator_path; order=order)
28 g = equidistant_grid(N, 0., 1.)
29 c = rand(size(g)...)
30 c_lz = eval_on(g, x -> 0.5)
31 D = second_derivative_variable(g, c, stencil_set)
32 D_lz = second_derivative_variable(g, c_lz, stencil_set)
33 u = rand(size(g)...)
34 u_xx = rand(size(g)...)
35
36 b_naive = @benchmark $u_xx .= $D*$u
37 b_naive_lz = @benchmark $u_xx .= $D_lz*$u
38 b_reg = @benchmark $apply_region_1d!($u_xx,$u,$D)
39 b_reg_lz = @benchmark $apply_region_1d!($u_xx,$u,$D_lz)
40 b_thrd = @benchmark $apply_region_threaded_1d!($u_xx,$u,$D)
41 b_thrd_lz = @benchmark $apply_region_threaded_1d!($u_xx,$u,$D_lz)
42 print_benchmark_result("Naive apply",b_naive)
43 print_benchmark_result("Naive apply lazy coeff",b_naive_lz)
44 print_benchmark_result("Region apply",b_reg)
45 print_benchmark_result("Region apply lazy coeff",b_reg_lz)
46 print_benchmark_result("Threaded region apply",b_thrd)
47 print_benchmark_result("Threaded region apply lazy coeff",b_thrd_lz)
48 end
49
50 function benchmark_const_coeff_2d(;N = 100, order = 4)
51 stencil_set = read_stencil_set(operator_path; order=order)
52 g = equidistant_grid((N,N), (0.,0.,),(1.,1.))
53 D = Laplace(g, stencil_set)
54 u = rand(size(g)...)
55 u_xx = rand(size(g)...)
56 if order == 2
57 clz_sz = 1
58 elseif order == 4
59 clz_sz = 4
60 else
61 error()
62 end
63
64 b_naive = @benchmark $u_xx .= $D*$u
65 b_reg = @benchmark $apply_region_2d!($u_xx,$u,$D,$clz_sz)
66 b_thrd = @benchmark $apply_region_threaded_2d!($u_xx,$u,$D,$clz_sz)
67 print_benchmark_result("Naive apply",b_naive)
68 print_benchmark_result("Region apply",b_reg)
69 print_benchmark_result("Threaded region apply",b_thrd)
70 end
71
72 function benchmark_var_coeff_2d(;N = 100, order = 4)
73 stencil_set = read_stencil_set(operator_path; order=order)
74 g = equidistant_grid((N,N), (0.,0.,),(1.,1.))
75 c = rand(size(g)...)
76 c_lz = eval_on(g, x-> 0.5)
77 D = second_derivative_variable(g, c, stencil_set, 1) + second_derivative_variable(g, c, stencil_set, 2)
78 D_lz = second_derivative_variable(g, c_lz, stencil_set, 1) + second_derivative_variable(g, c_lz, stencil_set, 2)
79 u = rand(size(g)...)
80 u_xx = rand(size(g)...)
81
82 if order == 2
83 clz_sz = 1
84 elseif order == 4
85 clz_sz = 6
86 else
87 error()
88 end
89
90 # Check correctnesss
91 # u_xx .= D*u
92 # u_xx_tmp = zeros(size(u_xx)...)
93 # u_xx_tmp .= u_xx
94 # apply_region_threaded_2d!(u_xx, u, D, clz_sz)
95
96 # @show sum(abs.(u_xx_tmp .- u_xx))
97 # @show pointer(u_xx_tmp) == pointer(u_xxs
98
99
100 b_naive = @benchmark $u_xx .= $D*$u
101 b_naive_lz = @benchmark $u_xx .= $D_lz*$u
102 b_reg = @benchmark $apply_region_2d!($u_xx,$u,$D, $clz_sz)
103 b_reg_lz = @benchmark $apply_region_2d!($u_xx,$u,$D_lz, $clz_sz)
104 b_thrd = @benchmark $apply_region_threaded_2d!($u_xx,$u,$D, $clz_sz)
105 b_thrd_lz = @benchmark $apply_region_threaded_2d!($u_xx,$u,$D_lz, $clz_sz)
106 print_benchmark_result("Naive apply",b_naive)
107 print_benchmark_result("Naive apply lazy coeff",b_naive_lz)
108 print_benchmark_result("Region apply",b_reg)
109 print_benchmark_result("Region apply lazy coeff",b_reg_lz)
110 print_benchmark_result("Threaded region apply",b_thrd)
111 print_benchmark_result("Threaded region apply lazy coeff",b_thrd_lz)
112 end
113
114 function print_benchmark_result(title_str,res)
115 if title_str[1] != ' '
116 title_str = lpad(title_str,length(title_str)+1, " ")
117 end
118 if title_str[end] != ' '
119 title_str = rpad(title_str,length(title_str)+1, " ")
120 end
121 tot_len = 76
122 pad_len = Int(tot_len/2)
123 header = lpad(title_str,pad_len,"*")
124 header = rpad(header,tot_len,"*")
125 bottom = repeat("*",tot_len)
126 println(header)
127 display(res)
128 println(bottom)
129 return
130 end
131
132 function apply_region_1d!(u_xx, u, D)
133 clz_sz = SbpOperators.closure_size(D)
134 tm = D*u
135 for i ∈ @view eachindex(u)[1:clz_sz]
136 u_xx[i] = tm[Index{Lower}(i)]
137 end
138 for i ∈ @view eachindex(u)[clz_sz+1:end-clz_sz]
139 u_xx[i] = tm[Index{Interior}(i)]
140 end
141 for i ∈ @view eachindex(u)[end-clz_sz+1:end]
142 u_xx[i] = tm[Index{Upper}(i)]
143 end
144 end
145
146 function apply_region_threaded_1d!(u_xx, u, D)
147 clz_sz = SbpOperators.closure_size(D)
148 tm = D*u
149 for i ∈ @view eachindex(u)[1:clz_sz]
150 u_xx[i] = tm[Index{Lower}(i)]
151 end
152 Threads.@threads for i ∈ @view eachindex(u)[clz_sz+1:end-clz_sz]
153 u_xx[i] = tm[Index{Interior}(i)]
154 end
155 for i ∈ @view eachindex(u)[end-clz_sz+1:end]
156 u_xx[i] = tm[Index{Upper}(i)]
157 end
158 end
159
160 function apply_region_2d!(u_xx, u, D, clz_sz)
161 tm = D*u
162 for I ∈ @view CartesianIndices(u)[1:clz_sz,1:clz_sz]
163 u_xx[I] = tm[Index{Lower}(I[1]),Index{Lower}(I[2])]
164 end
165 for I ∈ @view CartesianIndices(u)[1:clz_sz,clz_sz+1:end-clz_sz]
166 u_xx[I] = tm[Index{Lower}(I[1]),Index{Interior}(I[2])]
167 end
168 for I ∈ @view CartesianIndices(u)[1:clz_sz,end-clz_sz+1:end]
169 u_xx[I] = tm[Index{Lower}(I[1]),Index{Upper}(I[2])]
170 end
171 for I ∈ @view CartesianIndices(u)[clz_sz+1:end-clz_sz,1:clz_sz]
172 u_xx[I] = tm[Index{Interior}(I[1]),Index{Lower}(I[2])]
173 end
174 for I ∈ @view CartesianIndices(u)[clz_sz+1:end-clz_sz,clz_sz+1:end-clz_sz]
175 u_xx[I] = tm[Index{Interior}(I[1]),Index{Interior}(I[2])]
176 end
177 for I ∈ @view CartesianIndices(u)[clz_sz+1:end-clz_sz,end-clz_sz+1:end]
178 u_xx[I] = tm[Index{Interior}(I[1]),Index{Upper}(I[2])]
179 end
180 for I ∈ @view CartesianIndices(u)[end-clz_sz+1:end,1:clz_sz]
181 u_xx[I] = tm[Index{Upper}(I[1]),Index{Lower}(I[2])]
182 end
183 for I ∈ @view CartesianIndices(u)[end-clz_sz+1:end,clz_sz+1:end-clz_sz]
184 u_xx[I] = tm[Index{Upper}(I[1]),Index{Interior}(I[2])]
185 end
186 for I ∈ @view CartesianIndices(u)[end-clz_sz+1:end,end-clz_sz+1:end]
187 u_xx[I] = tm[Index{Upper}(I[1]),Index{Upper}(I[2])]
188 end
189 end
190
191 function apply_region_threaded_2d!(u_xx, u, D, clz_sz)
192 tm = D*u
193 for I ∈ @view CartesianIndices(u)[1:clz_sz,1:clz_sz]
194 u_xx[I] = tm[Index{Lower}(I[1]),Index{Lower}(I[2])]
195 end
196 for I ∈ @view CartesianIndices(u)[1:clz_sz,clz_sz+1:end-clz_sz]
197 u_xx[I] = tm[Index{Lower}(I[1]),Index{Interior}(I[2])]
198 end
199 for I ∈ @view CartesianIndices(u)[1:clz_sz,end-clz_sz+1:end]
200 u_xx[I] = tm[Index{Lower}(I[1]),Index{Upper}(I[2])]
201 end
202 for I ∈ @view CartesianIndices(u)[clz_sz+1:end-clz_sz,1:clz_sz]
203 u_xx[I] = tm[Index{Interior}(I[1]),Index{Lower}(I[2])]
204 end
205 Threads.@threads for I ∈ @view CartesianIndices(u)[clz_sz+1:end-clz_sz,clz_sz+1:end-clz_sz]
206 u_xx[I] = tm[Index{Interior}(I[1]),Index{Interior}(I[2])]
207 end
208 for I ∈ @view CartesianIndices(u)[clz_sz+1:end-clz_sz,end-clz_sz+1:end]
209 u_xx[I] = tm[Index{Interior}(I[1]),Index{Upper}(I[2])]
210 end
211 for I ∈ @view CartesianIndices(u)[end-clz_sz+1:end,1:clz_sz]
212 u_xx[I] = tm[Index{Upper}(I[1]),Index{Lower}(I[2])]
213 end
214 for I ∈ @view CartesianIndices(u)[end-clz_sz+1:end,clz_sz+1:end-clz_sz]
215 u_xx[I] = tm[Index{Upper}(I[1]),Index{Interior}(I[2])]
216 end
217 for I ∈ @view CartesianIndices(u)[end-clz_sz+1:end,end-clz_sz+1:end]
218 u_xx[I] = tm[Index{Upper}(I[1]),Index{Upper}(I[2])]
219 end
220 end