Mercurial > repos > public > sbplib_julia

comparison test/testSbpOperators.jl @ 633:a78bda7084f6 feature/quadrature_as_outer_product

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge w. default
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 01 Jan 2021 16:34:55 +0100
parents 08e27dee76c3 eaa8c852ddf2
children fb5ac62563aa
comparison
equal deleted inserted replaced
561:04d7b4eb63ef 633:a78bda7084f6
2 using Sbplib.SbpOperators 2 using Sbplib.SbpOperators
3 using Sbplib.Grids 3 using Sbplib.Grids
4 using Sbplib.RegionIndices 4 using Sbplib.RegionIndices
5 using Sbplib.LazyTensors 5 using Sbplib.LazyTensors
6 using LinearAlgebra 6 using LinearAlgebra
7 using TOML
8
9 import Sbplib.SbpOperators.Stencil
7 10
8 @testset "SbpOperators" begin 11 @testset "SbpOperators" begin
9 12
13 @testset "Stencil" begin
14 s = Stencil((-2,2), (1.,2.,2.,3.,4.))
15 @test s isa Stencil{Float64, 5}
16
17 @test eltype(s) == Float64
18 @test SbpOperators.scale(s, 2) == Stencil((-2,2), (2.,4.,4.,6.,8.))
19
20 @test Stencil((1,2,3,4), center=1) == Stencil((0, 3),(1,2,3,4))
21 @test Stencil((1,2,3,4), center=2) == Stencil((-1, 2),(1,2,3,4))
22 @test Stencil((1,2,3,4), center=4) == Stencil((-3, 0),(1,2,3,4))
23 end
24
25 @testset "parse_rational" begin
26 @test SbpOperators.parse_rational("1") isa Rational
27 @test SbpOperators.parse_rational("1") == 1//1
28 @test SbpOperators.parse_rational("1/2") isa Rational
29 @test SbpOperators.parse_rational("1/2") == 1//2
30 @test SbpOperators.parse_rational("37/13") isa Rational
31 @test SbpOperators.parse_rational("37/13") == 37//13
32 end
33
34 @testset "readoperator" begin
35 toml_str = """
36 [meta]
37 type = "equidistant"
38
39 [order2]
40 H.inner = ["1"]
41
42 D1.inner_stencil = ["-1/2", "0", "1/2"]
43 D1.closure_stencils = [
44 ["-1", "1"],
45 ]
46
47 d1.closure = ["-3/2", "2", "-1/2"]
48
49 [order4]
50 H.closure = ["17/48", "59/48", "43/48", "49/48"]
51
52 D2.inner_stencil = ["-1/12","4/3","-5/2","4/3","-1/12"]
53 D2.closure_stencils = [
54 [ "2", "-5", "4", "-1", "0", "0"],
55 [ "1", "-2", "1", "0", "0", "0"],
56 [ "-4/43", "59/43", "-110/43", "59/43", "-4/43", "0"],
57 [ "-1/49", "0", "59/49", "-118/49", "64/49", "-4/49"],
58 ]
59 """
60
61 parsed_toml = TOML.parse(toml_str)
62 @testset "get_stencil" begin
63 @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil") == Stencil((-1/2, 0., 1/2), center=2)
64 @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil", center=1) == Stencil((-1/2, 0., 1/2); center=1)
65 @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil", center=3) == Stencil((-1/2, 0., 1/2); center=3)
66
67 @test get_stencil(parsed_toml, "order2", "H", "inner") == Stencil((1.,), center=1)
68
69 @test_throws AssertionError get_stencil(parsed_toml, "meta", "type")
70 @test_throws AssertionError get_stencil(parsed_toml, "order2", "D1", "closure_stencils")
71 end
72
73 @testset "get_stencils" begin
74 @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=(1,)) == (Stencil((-1., 1.), center=1),)
75 @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=(2,)) == (Stencil((-1., 1.), center=2),)
76 @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=[2]) == (Stencil((-1., 1.), center=2),)
77
78 @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=[1,1,1,1]) == (
79 Stencil(( 2., -5., 4., -1., 0., 0.), center=1),
80 Stencil(( 1., -2., 1., 0., 0., 0.), center=1),
81 Stencil(( -4/43, 59/43, -110/43, 59/43, -4/43, 0.), center=1),
82 Stencil(( -1/49, 0., 59/49, -118/49, 64/49, -4/49), center=1),
83 )
84
85 @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(4,2,3,1)) == (
86 Stencil(( 2., -5., 4., -1., 0., 0.), center=4),
87 Stencil(( 1., -2., 1., 0., 0., 0.), center=2),
88 Stencil(( -4/43, 59/43, -110/43, 59/43, -4/43, 0.), center=3),
89 Stencil(( -1/49, 0., 59/49, -118/49, 64/49, -4/49), center=1),
90 )
91
92 @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=1:4) == (
93 Stencil(( 2., -5., 4., -1., 0., 0.), center=1),
94 Stencil(( 1., -2., 1., 0., 0., 0.), center=2),
95 Stencil(( -4/43, 59/43, -110/43, 59/43, -4/43, 0.), center=3),
96 Stencil(( -1/49, 0., 59/49, -118/49, 64/49, -4/49), center=4),
97 )
98
99 @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(1,2,3))
100 @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(1,2,3,5,4))
101 @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "inner_stencil",centers=(1,2))
102 end
103
104 @testset "get_tuple" begin
105 @test get_tuple(parsed_toml, "order2", "d1", "closure") == (-3/2, 2, -1/2)
106
107 @test_throws AssertionError get_tuple(parsed_toml, "meta", "type")
108 end
109 end
110
10 # @testset "apply_quadrature" begin 111 # @testset "apply_quadrature" begin
11 # op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 112 # op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
12 # h = 0.5 113 # h = 0.5
13 # 114 #
14 # @test apply_quadrature(op, h, 1.0, 10, 100) == h 115 # @test apply_quadrature(op, h, 1.0, 10, 100) == h
15 # 116 #
16 # N = 10 117 # N = 10
27 # @test apply_quadrature(op, h, v[i], i, N) == q[i]*v[i] 128 # @test apply_quadrature(op, h, v[i], i, N) == q[i]*v[i]
28 # end 129 # end
29 # end 130 # end
30 131
31 @testset "SecondDerivative" begin 132 @testset "SecondDerivative" begin
32 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 133 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
33 L = 3.5 134 L = 3.5
34 g = EquidistantGrid(101, 0.0, L) 135 g = EquidistantGrid(101, 0.0, L)
35 Dₓₓ = SecondDerivative(g,op.innerStencil,op.closureStencils) 136 Dₓₓ = SecondDerivative(g,op.innerStencil,op.closureStencils)
36 137
37 f0(x) = 1. 138 f0(x) = 1.
67 @test Dₓₓ*v5 ≈ -v5 atol=5e-4 norm=l2 168 @test Dₓₓ*v5 ≈ -v5 atol=5e-4 norm=l2
68 end 169 end
69 170
70 171
71 @testset "Laplace2D" begin 172 @testset "Laplace2D" begin
72 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 173 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
73 Lx = 1.5 174 Lx = 1.5
74 Ly = 3.2 175 Ly = 3.2
75 g = EquidistantGrid((102,131), (0.0, 0.0), (Lx,Ly)) 176 g = EquidistantGrid((102,131), (0.0, 0.0), (Lx,Ly))
76 L = Laplace(g, op.innerStencil, op.closureStencils) 177 L = Laplace(g, op.innerStencil, op.closureStencils)
77 178
109 @test L*v4 ≈ v2 atol=5e-4 norm=l2 210 @test L*v4 ≈ v2 atol=5e-4 norm=l2
110 @test L*v5 ≈ v5ₓₓ atol=5e-4 norm=l2 211 @test L*v5 ≈ v5ₓₓ atol=5e-4 norm=l2
111 end 212 end
112 213
113 @testset "DiagonalQuadrature" begin 214 @testset "DiagonalQuadrature" begin
114 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 215 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
115 Lx = π/2. 216 Lx = π/2.
116 Ly = Float64(π) 217 Ly = Float64(π)
117 g_1D = EquidistantGrid(77, 0.0, Lx) 218 g_1D = EquidistantGrid(77, 0.0, Lx)
118 g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) 219 g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
119 integral(H,v) = sum(H*v) 220 integral(H,v) = sum(H*v)
167 f_i(x) = 1/factorial(i)*x^i 268 f_i(x) = 1/factorial(i)*x^i
168 v = (v...,evalOn(g_1D,f_i)) 269 v = (v...,evalOn(g_1D,f_i))
169 end 270 end
170 # TODO: Bug in readOperator for 2nd order 271 # TODO: Bug in readOperator for 2nd order
171 # # 2nd order 272 # # 2nd order
172 # op2 = readOperator(sbp_operators_path()*"d2_2nd.txt",sbp_operators_path()*"h_2nd.txt") 273 # op2 = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
173 # H2 = diagonal_quadrature(g_1D,op2.quadratureClosure) 274 # H2 = diagonal_quadrature(g_1D,op2.quadratureClosure)
174 # for i = 1:3 275 # for i = 1:3
175 # @test integral(H2,v[i]) ≈ v[i+1] rtol = 1e-14 276 # @test integral(H2,v[i]) ≈ v[i+1] rtol = 1e-14
176 # end 277 # end
177 278
178 # 4th order 279 # 4th order
179 op4 = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 280 op4 = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
180 H4 = diagonal_quadrature(g_1D,op4.quadratureClosure) 281 H4 = diagonal_quadrature(g_1D,op4.quadratureClosure)
181 for i = 1:4 282 for i = 1:4
182 @test integral(H4,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14 283 @test integral(H4,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14
183 end 284 end
184 end 285 end
194 @inferred H_xy'*v_2D 295 @inferred H_xy'*v_2D
195 end 296 end
196 end 297 end
197 298
198 @testset "InverseDiagonalQuadrature" begin 299 @testset "InverseDiagonalQuadrature" begin
199 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 300 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
200 Lx = π/2. 301 Lx = π/2.
201 Ly = Float64(π) 302 Ly = Float64(π)
202 g_1D = EquidistantGrid(77, 0.0, Lx) 303 g_1D = EquidistantGrid(77, 0.0, Lx)
203 g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) 304 g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
204 @testset "Constructors" begin 305 @testset "Constructors" begin
255 @inferred Hi_x'*v_1D 356 @inferred Hi_x'*v_1D
256 @inferred Hi_xy*v_2D 357 @inferred Hi_xy*v_2D
257 @inferred Hi_xy'*v_2D 358 @inferred Hi_xy'*v_2D
258 end 359 end
259 end 360 end
260 # 361
261 # @testset "BoundaryValue" begin 362 @testset "BoundaryRestrictrion" begin
262 # op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 363 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
263 # g = EquidistantGrid((4,5), (0.0, 0.0), (1.0,1.0)) 364 g_1D = EquidistantGrid(11, 0.0, 1.0)
264 # 365 g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))
265 # e_w = BoundaryValue(op, g, CartesianBoundary{1,Lower}()) 366
266 # e_e = BoundaryValue(op, g, CartesianBoundary{1,Upper}()) 367 @testset "Constructors" begin
267 # e_s = BoundaryValue(op, g, CartesianBoundary{2,Lower}()) 368 @testset "1D" begin
268 # e_n = BoundaryValue(op, g, CartesianBoundary{2,Upper}()) 369 e_l = BoundaryRestriction{Lower}(op.eClosure,size(g_1D)[1])
269 # 370 @test e_l == BoundaryRestriction(g_1D,op.eClosure,Lower())
270 # v = zeros(Float64, 4, 5) 371 @test e_l == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Lower}())
271 # v[:,5] = [1, 2, 3,4] 372 @test e_l isa TensorMapping{T,0,1} where T
272 # v[:,4] = [1, 2, 3,4] 373
273 # v[:,3] = [4, 5, 6, 7] 374 e_r = BoundaryRestriction{Upper}(op.eClosure,size(g_1D)[1])
274 # v[:,2] = [7, 8, 9, 10] 375 @test e_r == BoundaryRestriction(g_1D,op.eClosure,Upper())
275 # v[:,1] = [10, 11, 12, 13] 376 @test e_r == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Upper}())
276 # 377 @test e_r isa TensorMapping{T,0,1} where T
277 # @test e_w isa TensorMapping{T,2,1} where T 378 end
278 # @test e_w' isa TensorMapping{T,1,2} where T 379
279 # 380 @testset "2D" begin
280 # @test domain_size(e_w, (3,2)) == (2,) 381 e_w = boundary_restriction(g_2D,op.eClosure,CartesianBoundary{1,Upper}())
281 # @test domain_size(e_e, (3,2)) == (2,) 382 @test e_w isa InflatedTensorMapping
282 # @test domain_size(e_s, (3,2)) == (3,) 383 @test e_w isa TensorMapping{T,1,2} where T
283 # @test domain_size(e_n, (3,2)) == (3,) 384 end
284 # 385 end
285 # @test size(e_w'*v) == (5,) 386
286 # @test size(e_e'*v) == (5,) 387 e_l = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Lower}())
287 # @test size(e_s'*v) == (4,) 388 e_r = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Upper}())
288 # @test size(e_n'*v) == (4,) 389
289 # 390 e_w = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{1,Lower}())
290 # @test e_w'*v == [10,7,4,1.0,1] 391 e_e = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{1,Upper}())
291 # @test e_e'*v == [13,10,7,4,4.0] 392 e_s = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{2,Lower}())
292 # @test e_s'*v == [10,11,12,13.0] 393 e_n = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{2,Upper}())
293 # @test e_n'*v == [1,2,3,4.0] 394
294 # 395 @testset "Sizes" begin
295 # g_x = [1,2,3,4.0] 396 @testset "1D" begin
296 # g_y = [5,4,3,2,1.0] 397 @test domain_size(e_l) == (11,)
297 # 398 @test domain_size(e_r) == (11,)
298 # G_w = zeros(Float64, (4,5)) 399
299 # G_w[1,:] = g_y 400 @test range_size(e_l) == ()
300 # 401 @test range_size(e_r) == ()
301 # G_e = zeros(Float64, (4,5)) 402 end
302 # G_e[4,:] = g_y 403
303 # 404 @testset "2D" begin
304 # G_s = zeros(Float64, (4,5)) 405 @test domain_size(e_w) == (11,15)
305 # G_s[:,1] = g_x 406 @test domain_size(e_e) == (11,15)
306 # 407 @test domain_size(e_s) == (11,15)
307 # G_n = zeros(Float64, (4,5)) 408 @test domain_size(e_n) == (11,15)
308 # G_n[:,5] = g_x 409
309 # 410 @test range_size(e_w) == (15,)
310 # @test size(e_w*g_y) == (UnknownDim,5) 411 @test range_size(e_e) == (15,)
311 # @test size(e_e*g_y) == (UnknownDim,5) 412 @test range_size(e_s) == (11,)
312 # @test size(e_s*g_x) == (4,UnknownDim) 413 @test range_size(e_n) == (11,)
313 # @test size(e_n*g_x) == (4,UnknownDim) 414 end
314 # 415 end
315 # # These tests should be moved to where they are possible (i.e we know what the grid should be) 416
316 # @test_broken e_w*g_y == G_w 417
317 # @test_broken e_e*g_y == G_e 418 @testset "Application" begin
318 # @test_broken e_s*g_x == G_s 419 @testset "1D" begin
319 # @test_broken e_n*g_x == G_n 420 v = evalOn(g_1D,x->1+x^2)
320 # end 421 u = fill(3.124)
422 @test (e_l*v)[] == v[1]
423 @test (e_r*v)[] == v[end]
424 @test (e_r*v)[1] == v[end]
425 @test e_l'*u == [u[]; zeros(10)]
426 @test e_r'*u == [zeros(10); u[]]
427 end
428
429 @testset "2D" begin
430 v = rand(11, 15)
431 u = fill(3.124)
432
433 @test e_w*v == v[1,:]
434 @test e_e*v == v[end,:]
435 @test e_s*v == v[:,1]
436 @test e_n*v == v[:,end]
437
438
439 g_x = rand(11)
440 g_y = rand(15)
441
442 G_w = zeros(Float64, (11,15))
443 G_w[1,:] = g_y
444
445 G_e = zeros(Float64, (11,15))
446 G_e[end,:] = g_y
447
448 G_s = zeros(Float64, (11,15))
449 G_s[:,1] = g_x
450
451 G_n = zeros(Float64, (11,15))
452 G_n[:,end] = g_x
453
454 @test e_w'*g_y == G_w
455 @test e_e'*g_y == G_e
456 @test e_s'*g_x == G_s
457 @test e_n'*g_x == G_n
458 end
459
460 @testset "Regions" begin
461 u = fill(3.124)
462 @test (e_l'*u)[Index(1,Lower)] == 3.124
463 @test (e_l'*u)[Index(2,Lower)] == 0
464 @test (e_l'*u)[Index(6,Interior)] == 0
465 @test (e_l'*u)[Index(10,Upper)] == 0
466 @test (e_l'*u)[Index(11,Upper)] == 0
467
468 @test (e_r'*u)[Index(1,Lower)] == 0
469 @test (e_r'*u)[Index(2,Lower)] == 0
470 @test (e_r'*u)[Index(6,Interior)] == 0
471 @test (e_r'*u)[Index(10,Upper)] == 0
472 @test (e_r'*u)[Index(11,Upper)] == 3.124
473 end
474 end
475
476 @testset "Inferred" begin
477 v = ones(Float64, 11)
478 u = fill(1.)
479
480 @inferred apply(e_l, v)
481 @inferred apply(e_r, v)
482
483 @inferred apply_transpose(e_l, u, 4)
484 @inferred apply_transpose(e_l, u, Index(1,Lower))
485 @inferred apply_transpose(e_l, u, Index(2,Lower))
486 @inferred apply_transpose(e_l, u, Index(6,Interior))
487 @inferred apply_transpose(e_l, u, Index(10,Upper))
488 @inferred apply_transpose(e_l, u, Index(11,Upper))
489
490 @inferred apply_transpose(e_r, u, 4)
491 @inferred apply_transpose(e_r, u, Index(1,Lower))
492 @inferred apply_transpose(e_r, u, Index(2,Lower))
493 @inferred apply_transpose(e_r, u, Index(6,Interior))
494 @inferred apply_transpose(e_r, u, Index(10,Upper))
495 @inferred apply_transpose(e_r, u, Index(11,Upper))
496 end
497
498 end
321 # 499 #
322 # @testset "NormalDerivative" begin 500 # @testset "NormalDerivative" begin
323 # op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 501 # op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
324 # g = EquidistantGrid((5,6), (0.0, 0.0), (4.0,5.0)) 502 # g = EquidistantGrid((5,6), (0.0, 0.0), (4.0,5.0))
325 # 503 #
326 # d_w = NormalDerivative(op, g, CartesianBoundary{1,Lower}()) 504 # d_w = NormalDerivative(op, g, CartesianBoundary{1,Lower}())
327 # d_e = NormalDerivative(op, g, CartesianBoundary{1,Upper}()) 505 # d_e = NormalDerivative(op, g, CartesianBoundary{1,Upper}())
328 # d_s = NormalDerivative(op, g, CartesianBoundary{2,Lower}()) 506 # d_s = NormalDerivative(op, g, CartesianBoundary{2,Lower}())
395 # @test_broken d_s*g_x .≈ G_s 573 # @test_broken d_s*g_x .≈ G_s
396 # @test_broken d_n*g_x .≈ G_n 574 # @test_broken d_n*g_x .≈ G_n
397 # end 575 # end
398 # 576 #
399 # @testset "BoundaryQuadrature" begin 577 # @testset "BoundaryQuadrature" begin
400 # op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 578 # op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
401 # g = EquidistantGrid((10,11), (0.0, 0.0), (1.0,1.0)) 579 # g = EquidistantGrid((10,11), (0.0, 0.0), (1.0,1.0))
402 # 580 #
403 # H_w = BoundaryQuadrature(op, g, CartesianBoundary{1,Lower}()) 581 # H_w = BoundaryQuadrature(op, g, CartesianBoundary{1,Lower}())
404 # H_e = BoundaryQuadrature(op, g, CartesianBoundary{1,Upper}()) 582 # H_e = BoundaryQuadrature(op, g, CartesianBoundary{1,Upper}())
405 # H_s = BoundaryQuadrature(op, g, CartesianBoundary{2,Lower}()) 583 # H_s = BoundaryQuadrature(op, g, CartesianBoundary{2,Lower}())