Mercurial > repos > public > sbplib_julia

comparison test/testSbpOperators.jl @ 610:e40e7439d1b4 feature/volume_and_boundary_operators

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add a general boundary operator and make BoundaryRestriction a specialization of it.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Sat, 05 Dec 2020 18:12:31 +0100
parents 8e4f86c4bf75
children 1db945cba3a2
comparison
equal deleted inserted replaced
606:8f9b3eac128a 610:e40e7439d1b4
179 @test Qinv' isa TensorMapping{T,2,2} where T 179 @test Qinv' isa TensorMapping{T,2,2} where T
180 @test_broken Qinv*(Q*v) ≈ v 180 @test_broken Qinv*(Q*v) ≈ v
181 @test Qinv*v == Qinv'*v 181 @test Qinv*v == Qinv'*v
182 end 182 end
183 183
184 @testset "BoundaryRestrictrion" begin 184 @testset "BoundaryOperator" begin
185 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 185 closure_stencil = SbpOperators.Stencil((0,2), (2.,1.,3.))
186 g_1D = EquidistantGrid(11, 0.0, 1.0) 186 g_1D = EquidistantGrid(11, 0.0, 1.0)
187 g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0)) 187 g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))
188 188
189 @testset "Constructors" begin 189 @testset "Constructors" begin
190 @testset "1D" begin 190 @testset "1D" begin
191 e_l = BoundaryRestriction{Lower}(op.eClosure,size(g_1D)[1]) 191 op_l = BoundaryOperator{Lower}(closure_stencil,size(g_1D)[1])
192 @test e_l == BoundaryRestriction(g_1D,op.eClosure,Lower()) 192 @test op_l == BoundaryOperator(g_1D,closure_stencil,Lower())
193 @test e_l == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Lower}()) 193 @test op_l == boundary_operator(g_1D,closure_stencil,CartesianBoundary{1,Lower}())
194 @test e_l isa TensorMapping{T,0,1} where T 194 @test op_l isa TensorMapping{T,0,1} where T
195 195
196 e_r = BoundaryRestriction{Upper}(op.eClosure,size(g_1D)[1]) 196 op_r = BoundaryOperator{Upper}(closure_stencil,size(g_1D)[1])
197 @test e_r == BoundaryRestriction(g_1D,op.eClosure,Upper()) 197 @test op_r == BoundaryRestriction(g_1D,closure_stencil,Upper())
198 @test e_r == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Upper}()) 198 @test op_r == boundary_operator(g_1D,closure_stencil,CartesianBoundary{1,Upper}())
199 @test e_r isa TensorMapping{T,0,1} where T 199 @test op_r isa TensorMapping{T,0,1} where T
200 end 200 end
201 201
202 @testset "2D" begin 202 @testset "2D" begin
203 e_w = boundary_restriction(g_2D,op.eClosure,CartesianBoundary{1,Upper}()) 203 e_w = boundary_operator(g_2D,closure_stencil,CartesianBoundary{1,Upper}())
204 @test e_w isa InflatedTensorMapping 204 @test e_w isa InflatedTensorMapping
205 @test e_w isa TensorMapping{T,1,2} where T 205 @test e_w isa TensorMapping{T,1,2} where T
206 end 206 end
207 end 207 end
208 208
209 e_l = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Lower}()) 209 op_l = boundary_operator(g_1D, closure_stencil, CartesianBoundary{1,Lower}())
210 e_r = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Upper}()) 210 op_r = boundary_operator(g_1D, closure_stencil, CartesianBoundary{1,Upper}())
211 211
212 e_w = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{1,Lower}()) 212 op_w = boundary_operator(g_2D, closure_stencil, CartesianBoundary{1,Lower}())
213 e_e = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{1,Upper}()) 213 op_e = boundary_operator(g_2D, closure_stencil, CartesianBoundary{1,Upper}())
214 e_s = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{2,Lower}()) 214 op_s = boundary_operator(g_2D, closure_stencil, CartesianBoundary{2,Lower}())
215 e_n = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{2,Upper}()) 215 op_n = boundary_operator(g_2D, closure_stencil, CartesianBoundary{2,Upper}())
216 216
217 @testset "Sizes" begin 217 @testset "Sizes" begin
218 @testset "1D" begin 218 @testset "1D" begin
219 @test domain_size(e_l) == (11,) 219 @test domain_size(op_l) == (11,)
220 @test domain_size(e_r) == (11,) 220 @test domain_size(op_r) == (11,)
221 221
222 @test range_size(e_l) == () 222 @test range_size(op_l) == ()
223 @test range_size(e_r) == () 223 @test range_size(op_r) == ()
224 end 224 end
225 225
226 @testset "2D" begin 226 @testset "2D" begin
227 @test domain_size(e_w) == (11,15) 227 @test domain_size(op_w) == (11,15)
228 @test domain_size(e_e) == (11,15) 228 @test domain_size(op_e) == (11,15)
229 @test domain_size(e_s) == (11,15) 229 @test domain_size(op_s) == (11,15)
230 @test domain_size(e_n) == (11,15) 230 @test domain_size(op_n) == (11,15)
231 231
232 @test range_size(e_w) == (15,) 232 @test range_size(op_w) == (15,)
233 @test range_size(e_e) == (15,) 233 @test range_size(op_e) == (15,)
234 @test range_size(e_s) == (11,) 234 @test range_size(op_s) == (11,)
235 @test range_size(e_n) == (11,) 235 @test range_size(op_n) == (11,)
236 end 236 end
237 end 237 end
238 238
239 239
240 @testset "Application" begin 240 @testset "Application" begin
241 @testset "1D" begin 241 @testset "1D" begin
242 v = evalOn(g_1D,x->1+x^2) 242 v = evalOn(g_1D,x->1+x^2)
243 u = fill(3.124) 243 u = fill(3.124)
244 @test (op_l*v)[] == 2*v[1] + v[2] + 3*v[3]
245 @test (op_r*v)[] == 2*v[end] + v[end-1] + 3*v[end-2]
246 @test (op_r*v)[1] == 2*v[end] + v[end-1] + 3*v[end-2]
247 @test op_l'*u == [2*u[]; u[]; 3*u[]; zeros(8)]
248 @test op_r'*u == [zeros(8); 3*u[]; u[]; 2*u[]]
249 end
250
251 @testset "2D" begin
252 v = rand(size(g_2D)...)
253 u = fill(3.124)
254 @test op_w*v ≈ 2*v[1,:] + v[2,:] + 3*v[3,:] rtol = 1e-14
255 @test op_e*v ≈ 2*v[end,:] + v[end-1,:] + 3*v[end-2,:] rtol = 1e-14
256 @test op_s*v ≈ 2*v[:,1] + v[:,2] + 3*v[:,3] rtol = 1e-14
257 @test op_n*v ≈ 2*v[:,end] + v[:,end-1] + 3*v[:,end-2] rtol = 1e-14
258
259
260 g_x = rand(size(g_2D)[1])
261 g_y = rand(size(g_2D)[2])
262
263 G_w = zeros(Float64, size(g_2D)...)
264 G_w[1,:] = 2*g_y
265 G_w[2,:] = g_y
266 G_w[3,:] = 3*g_y
267
268 G_e = zeros(Float64, size(g_2D)...)
269 G_e[end,:] = 2*g_y
270 G_e[end-1,:] = g_y
271 G_e[end-2,:] = 3*g_y
272
273 G_s = zeros(Float64, size(g_2D)...)
274 G_s[:,1] = 2*g_x
275 G_s[:,2] = g_x
276 G_s[:,3] = 3*g_x
277
278 G_n = zeros(Float64, size(g_2D)...)
279 G_n[:,end] = 2*g_x
280 G_n[:,end-1] = g_x
281 G_n[:,end-2] = 3*g_x
282
283 @test op_w'*g_y == G_w
284 @test op_e'*g_y == G_e
285 @test op_s'*g_x == G_s
286 @test op_n'*g_x == G_n
287 end
288
289 @testset "Regions" begin
290 u = fill(3.124)
291 @test (op_l'*u)[Index(1,Lower)] == 2*u[]
292 @test (op_l'*u)[Index(2,Lower)] == u[]
293 @test (op_l'*u)[Index(6,Interior)] == 0
294 @test (op_l'*u)[Index(10,Upper)] == 0
295 @test (op_l'*u)[Index(11,Upper)] == 0
296
297 @test (op_r'*u)[Index(1,Lower)] == 0
298 @test (op_r'*u)[Index(2,Lower)] == 0
299 @test (op_r'*u)[Index(6,Interior)] == 0
300 @test (op_r'*u)[Index(10,Upper)] == u[]
301 @test (op_r'*u)[Index(11,Upper)] == 2*u[]
302 end
303 end
304
305 @testset "Inferred" begin
306 v = ones(Float64, 11)
307 u = fill(1.)
308
309 @inferred apply(op_l, v)
310 @inferred apply(op_r, v)
311
312 @inferred apply_transpose(op_l, u, 4)
313 @inferred apply_transpose(op_l, u, Index(1,Lower))
314 @inferred apply_transpose(op_l, u, Index(2,Lower))
315 @inferred apply_transpose(op_l, u, Index(6,Interior))
316 @inferred apply_transpose(op_l, u, Index(10,Upper))
317 @inferred apply_transpose(op_l, u, Index(11,Upper))
318
319 @inferred apply_transpose(op_r, u, 4)
320 @inferred apply_transpose(op_r, u, Index(1,Lower))
321 @inferred apply_transpose(op_r, u, Index(2,Lower))
322 @inferred apply_transpose(op_r, u, Index(6,Interior))
323 @inferred apply_transpose(op_r, u, Index(10,Upper))
324 @inferred apply_transpose(op_r, u, Index(11,Upper))
325 end
326
327 end
328
329 @testset "BoundaryRestriction" begin
330 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
331 g_1D = EquidistantGrid(11, 0.0, 1.0)
332 g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))
333
334 @testset "Constructors" begin
335 @testset "1D" begin
336 e_l = BoundaryRestriction(g_1D,op.eClosure,Lower())
337 @test e_l == BoundaryRestriction(g_1D,op.eClosure,CartesianBoundary{1,Lower}())
338 @test e_l == BoundaryOperator(g_1D,op.eClosure,Lower())
339 @test e_l isa BoundaryOperator{T,Lower} where T
340 @test e_l isa TensorMapping{T,0,1} where T
341
342 e_r = BoundaryRestriction(g_1D,op.eClosure,Upper())
343 @test e_r == BoundaryRestriction(g_1D,op.eClosure,CartesianBoundary{1,Upper}())
344 @test e_r == BoundaryOperator(g_1D,op.eClosure,Upper())
345 @test e_r isa BoundaryOperator{T,Upper} where T
346 @test e_r isa TensorMapping{T,0,1} where T
347 end
348
349 @testset "2D" begin
350 e_w = BoundaryRestriction(g_2D,op.eClosure,CartesianBoundary{1,Upper}())
351 @test e_w isa InflatedTensorMapping
352 @test e_w isa TensorMapping{T,1,2} where T
353 end
354 end
355
356 @testset "Application" begin
357 @testset "1D" begin
358 e_l = BoundaryRestriction(g_1D, op.eClosure, CartesianBoundary{1,Lower}())
359 e_r = BoundaryRestriction(g_1D, op.eClosure, CartesianBoundary{1,Upper}())
360
361 v = evalOn(g_1D,x->1+x^2)
362 u = fill(3.124)
363
244 @test (e_l*v)[] == v[1] 364 @test (e_l*v)[] == v[1]
245 @test (e_r*v)[] == v[end] 365 @test (e_r*v)[] == v[end]
246 @test (e_r*v)[1] == v[end] 366 @test (e_r*v)[1] == v[end]
247 @test e_l'*u == [u[]; zeros(10)] 367 @test e_l'*u == [u[]; zeros(10)]
248 @test e_r'*u == [zeros(10); u[]] 368 @test e_r'*u == [zeros(10); u[]]
249 end 369 end
250 370
251 @testset "2D" begin 371 @testset "2D" begin
372 e_w = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{1,Lower}())
373 e_e = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{1,Upper}())
374 e_s = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{2,Lower}())
375 e_n = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{2,Upper}())
376
252 v = rand(11, 15) 377 v = rand(11, 15)
253 u = fill(3.124) 378 u = fill(3.124)
254 379
255 @test e_w*v == v[1,:] 380 @test e_w*v == v[1,:]
256 @test e_e*v == v[end,:] 381 @test e_e*v == v[end,:]
276 @test e_w'*g_y == G_w 401 @test e_w'*g_y == G_w
277 @test e_e'*g_y == G_e 402 @test e_e'*g_y == G_e
278 @test e_s'*g_x == G_s 403 @test e_s'*g_x == G_s
279 @test e_n'*g_x == G_n 404 @test e_n'*g_x == G_n
280 end 405 end
281
282 @testset "Regions" begin
283 u = fill(3.124)
284 @test (e_l'*u)[Index(1,Lower)] == 3.124
285 @test (e_l'*u)[Index(2,Lower)] == 0
286 @test (e_l'*u)[Index(6,Interior)] == 0
287 @test (e_l'*u)[Index(10,Upper)] == 0
288 @test (e_l'*u)[Index(11,Upper)] == 0
289
290 @test (e_r'*u)[Index(1,Lower)] == 0
291 @test (e_r'*u)[Index(2,Lower)] == 0
292 @test (e_r'*u)[Index(6,Interior)] == 0
293 @test (e_r'*u)[Index(10,Upper)] == 0
294 @test (e_r'*u)[Index(11,Upper)] == 3.124
295 end
296 end 406 end
297
298 @testset "Inferred" begin
299 v = ones(Float64, 11)
300 u = fill(1.)
301
302 @inferred apply(e_l, v)
303 @inferred apply(e_r, v)
304
305 @inferred apply_transpose(e_l, u, 4)
306 @inferred apply_transpose(e_l, u, Index(1,Lower))
307 @inferred apply_transpose(e_l, u, Index(2,Lower))
308 @inferred apply_transpose(e_l, u, Index(6,Interior))
309 @inferred apply_transpose(e_l, u, Index(10,Upper))
310 @inferred apply_transpose(e_l, u, Index(11,Upper))
311
312 @inferred apply_transpose(e_r, u, 4)
313 @inferred apply_transpose(e_r, u, Index(1,Lower))
314 @inferred apply_transpose(e_r, u, Index(2,Lower))
315 @inferred apply_transpose(e_r, u, Index(6,Interior))
316 @inferred apply_transpose(e_r, u, Index(10,Upper))
317 @inferred apply_transpose(e_r, u, Index(11,Upper))
318 end
319
320 end 407 end
321 # 408 #
322 # @testset "NormalDerivative" begin 409 # @testset "NormalDerivative" begin
323 # op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 410 # op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
324 # g = EquidistantGrid((5,6), (0.0, 0.0), (4.0,5.0)) 411 # g = EquidistantGrid((5,6), (0.0, 0.0), (4.0,5.0))