comparison test/SbpOperators/boundaryops/boundary_operator_test.jl @ 1207:f1c2a4fa0ee1 performance/get_region_type_inference

Merge default
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 03 Feb 2023 22:14:47 +0100
parents ae006e844870
children 1cc45207817e 102ebdaf7c11
comparison
equal deleted inserted replaced
919:b41180efb6c2 1207:f1c2a4fa0ee1
4 using Sbplib.SbpOperators 4 using Sbplib.SbpOperators
5 using Sbplib.Grids 5 using Sbplib.Grids
6 using Sbplib.RegionIndices 6 using Sbplib.RegionIndices
7 import Sbplib.SbpOperators.Stencil 7 import Sbplib.SbpOperators.Stencil
8 import Sbplib.SbpOperators.BoundaryOperator 8 import Sbplib.SbpOperators.BoundaryOperator
9 import Sbplib.SbpOperators.boundary_operator 9
10 10
11 @testset "BoundaryOperator" begin 11 @testset "BoundaryOperator" begin
12 closure_stencil = Stencil((0,2), (2.,1.,3.)) 12 closure_stencil = Stencil(2.,1.,3.; center = 1)
13 g_1D = EquidistantGrid(11, 0.0, 1.0) 13 g_1D = EquidistantGrid(11, 0.0, 1.0)
14 g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0)) 14 g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))
15 15
16 @testset "Constructors" begin 16 @testset "Constructors" begin
17 @testset "1D" begin 17 @test BoundaryOperator(g_1D, closure_stencil, Lower()) isa LazyTensor{T,0,1} where T
18 op_l = BoundaryOperator{Lower}(closure_stencil,size(g_1D)[1]) 18 @test BoundaryOperator(g_1D, closure_stencil, Upper()) isa LazyTensor{T,0,1} where T
19 @test op_l == BoundaryOperator(g_1D,closure_stencil,Lower())
20 @test op_l == boundary_operator(g_1D,closure_stencil,CartesianBoundary{1,Lower}())
21 @test op_l isa TensorMapping{T,0,1} where T
22
23 op_r = BoundaryOperator{Upper}(closure_stencil,size(g_1D)[1])
24 @test op_r == BoundaryOperator(g_1D,closure_stencil,Upper())
25 @test op_r == boundary_operator(g_1D,closure_stencil,CartesianBoundary{1,Upper}())
26 @test op_r isa TensorMapping{T,0,1} where T
27 end
28
29 @testset "2D" begin
30 e_w = boundary_operator(g_2D,closure_stencil,CartesianBoundary{1,Upper}())
31 @test e_w isa InflatedTensorMapping
32 @test e_w isa TensorMapping{T,1,2} where T
33 end
34 end 19 end
35 20
36 op_l = boundary_operator(g_1D, closure_stencil, CartesianBoundary{1,Lower}()) 21 op_l = BoundaryOperator(g_1D, closure_stencil, Lower())
37 op_r = boundary_operator(g_1D, closure_stencil, CartesianBoundary{1,Upper}()) 22 op_r = BoundaryOperator(g_1D, closure_stencil, Upper())
38
39 op_w = boundary_operator(g_2D, closure_stencil, CartesianBoundary{1,Lower}())
40 op_e = boundary_operator(g_2D, closure_stencil, CartesianBoundary{1,Upper}())
41 op_s = boundary_operator(g_2D, closure_stencil, CartesianBoundary{2,Lower}())
42 op_n = boundary_operator(g_2D, closure_stencil, CartesianBoundary{2,Upper}())
43 23
44 @testset "Sizes" begin 24 @testset "Sizes" begin
45 @testset "1D" begin 25 @test domain_size(op_l) == (11,)
46 @test domain_size(op_l) == (11,) 26 @test domain_size(op_r) == (11,)
47 @test domain_size(op_r) == (11,)
48 27
49 @test range_size(op_l) == () 28 @test range_size(op_l) == ()
50 @test range_size(op_r) == () 29 @test range_size(op_r) == ()
51 end
52
53 @testset "2D" begin
54 @test domain_size(op_w) == (11,15)
55 @test domain_size(op_e) == (11,15)
56 @test domain_size(op_s) == (11,15)
57 @test domain_size(op_n) == (11,15)
58
59 @test range_size(op_w) == (15,)
60 @test range_size(op_e) == (15,)
61 @test range_size(op_s) == (11,)
62 @test range_size(op_n) == (11,)
63 end
64 end 30 end
65 31
66 @testset "Application" begin 32 @testset "Application" begin
67 @testset "1D" begin 33 v = evalOn(g_1D,x->1+x^2)
68 v = evalOn(g_1D,x->1+x^2) 34 u = fill(3.124)
69 u = fill(3.124) 35 @test (op_l*v)[] == 2*v[1] + v[2] + 3*v[3]
70 @test (op_l*v)[] == 2*v[1] + v[2] + 3*v[3] 36 @test (op_r*v)[] == 2*v[end] + v[end-1] + 3*v[end-2]
71 @test (op_r*v)[] == 2*v[end] + v[end-1] + 3*v[end-2] 37 @test (op_r*v)[1] == 2*v[end] + v[end-1] + 3*v[end-2]
72 @test (op_r*v)[1] == 2*v[end] + v[end-1] + 3*v[end-2] 38 @test op_l'*u == [2*u[]; u[]; 3*u[]; zeros(8)]
73 @test op_l'*u == [2*u[]; u[]; 3*u[]; zeros(8)] 39 @test op_r'*u == [zeros(8); 3*u[]; u[]; 2*u[]]
74 @test op_r'*u == [zeros(8); 3*u[]; u[]; 2*u[]]
75 end
76 40
77 @testset "2D" begin 41 v = evalOn(g_1D, x->1. +x*im)
78 v = rand(size(g_2D)...) 42 @test (op_l*v)[] isa ComplexF64
79 u = fill(3.124)
80 @test op_w*v ≈ 2*v[1,:] + v[2,:] + 3*v[3,:] rtol = 1e-14
81 @test op_e*v ≈ 2*v[end,:] + v[end-1,:] + 3*v[end-2,:] rtol = 1e-14
82 @test op_s*v ≈ 2*v[:,1] + v[:,2] + 3*v[:,3] rtol = 1e-14
83 @test op_n*v ≈ 2*v[:,end] + v[:,end-1] + 3*v[:,end-2] rtol = 1e-14
84 43
44 u = fill(1. +im)
45 @test (op_l'*u)[1] isa ComplexF64
46 @test (op_l'*u)[5] isa ComplexF64
47 @test (op_l'*u)[11] isa ComplexF64
85 48
86 g_x = rand(size(g_2D)[1]) 49 u = fill(3.124)
87 g_y = rand(size(g_2D)[2]) 50 @test (op_l'*u)[Index(1,Lower)] == 2*u[]
51 @test (op_l'*u)[Index(2,Lower)] == u[]
52 @test (op_l'*u)[Index(6,Interior)] == 0
53 @test (op_l'*u)[Index(10,Upper)] == 0
54 @test (op_l'*u)[Index(11,Upper)] == 0
88 55
89 G_w = zeros(Float64, size(g_2D)...) 56 @test (op_r'*u)[Index(1,Lower)] == 0
90 G_w[1,:] = 2*g_y 57 @test (op_r'*u)[Index(2,Lower)] == 0
91 G_w[2,:] = g_y 58 @test (op_r'*u)[Index(6,Interior)] == 0
92 G_w[3,:] = 3*g_y 59 @test (op_r'*u)[Index(10,Upper)] == u[]
93 60 @test (op_r'*u)[Index(11,Upper)] == 2*u[]
94 G_e = zeros(Float64, size(g_2D)...)
95 G_e[end,:] = 2*g_y
96 G_e[end-1,:] = g_y
97 G_e[end-2,:] = 3*g_y
98
99 G_s = zeros(Float64, size(g_2D)...)
100 G_s[:,1] = 2*g_x
101 G_s[:,2] = g_x
102 G_s[:,3] = 3*g_x
103
104 G_n = zeros(Float64, size(g_2D)...)
105 G_n[:,end] = 2*g_x
106 G_n[:,end-1] = g_x
107 G_n[:,end-2] = 3*g_x
108
109 @test op_w'*g_y == G_w
110 @test op_e'*g_y == G_e
111 @test op_s'*g_x == G_s
112 @test op_n'*g_x == G_n
113 end
114
115 @testset "Regions" begin
116 u = fill(3.124)
117 @test (op_l'*u)[Index(1,Lower)] == 2*u[]
118 @test (op_l'*u)[Index(2,Lower)] == u[]
119 @test (op_l'*u)[Index(6,Interior)] == 0
120 @test (op_l'*u)[Index(10,Upper)] == 0
121 @test (op_l'*u)[Index(11,Upper)] == 0
122
123 @test (op_r'*u)[Index(1,Lower)] == 0
124 @test (op_r'*u)[Index(2,Lower)] == 0
125 @test (op_r'*u)[Index(6,Interior)] == 0
126 @test (op_r'*u)[Index(10,Upper)] == u[]
127 @test (op_r'*u)[Index(11,Upper)] == 2*u[]
128 end
129 end 61 end
130 62
131 @testset "Inferred" begin 63 @testset "Inferred" begin
132 v = ones(Float64, 11) 64 v = ones(Float64, 11)
133 u = fill(1.) 65 u = fill(1.)