Mercurial > repos > public > sbplib_julia

comparison test/SbpOperators/volumeops/laplace/laplace_test.jl @ 924:12e8e431b43c feature/laplace_opset

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Start restructuring Laplace making it more minimal.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 21 Feb 2022 13:12:47 +0100
parents 6a4d36eccf39
children 47425442bbc5
comparison
equal deleted inserted replaced
923:88bf50821cf5 924:12e8e431b43c
2 2
3 using Sbplib.SbpOperators 3 using Sbplib.SbpOperators
4 using Sbplib.Grids 4 using Sbplib.Grids
5 using Sbplib.LazyTensors 5 using Sbplib.LazyTensors
6 using Sbplib.RegionIndices 6 using Sbplib.RegionIndices
7 using Sbplib.StaticDicts
8 7
8 # Default stencils (4th order)
9 operator_path = sbp_operators_path()*"standard_diagonal.toml" 9 operator_path = sbp_operators_path()*"standard_diagonal.toml"
10 # Default stencils (4th order)
11 stencil_set = read_stencil_set(operator_path; order=4) 10 stencil_set = read_stencil_set(operator_path; order=4)
12 inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"]) 11 inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"])
13 closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"]) 12 closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
14 e_closure = parse_stencil(stencil_set["e"]["closure"])
15 d_closure = parse_stencil(stencil_set["d1"]["closure"])
16 quadrature_interior = parse_scalar(stencil_set["H"]["inner"])
17 quadrature_closure = parse_tuple(stencil_set["H"]["closure"])
18 13
19 @testset "Laplace" begin 14 @testset "Laplace" begin
20 g_1D = EquidistantGrid(101, 0.0, 1.) 15 g_1D = EquidistantGrid(101, 0.0, 1.)
21 g_3D = EquidistantGrid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.)) 16 g_3D = EquidistantGrid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.))
22 @testset "Constructors" begin 17 @testset "Constructors" begin
23 18
24 @testset "1D" begin 19 @testset "1D" begin
25 Δ = laplace(g_1D, inner_stencil, closure_stencils) 20 Δ = laplace(g_1D, inner_stencil, closure_stencils)
26 H = inner_product(g_1D, quadrature_interior, quadrature_closure) 21 @test Laplace(g_1D, stencil_set) == Laplace(Δ, stencil_set)
27 Hi = inverse_inner_product(g_1D, quadrature_interior, quadrature_closure) 22 @test Laplace(g_1D, stencil_set) isa TensorMapping{T,1,1} where T
28
29 (id_l, id_r) = boundary_identifiers(g_1D)
30
31 e_l = boundary_restriction(g_1D, e_closure,id_l)
32 e_r = boundary_restriction(g_1D, e_closure,id_r)
33 e_dict = StaticDict(id_l => e_l, id_r => e_r)
34
35 d_l = normal_derivative(g_1D, d_closure,id_l)
36 d_r = normal_derivative(g_1D, d_closure,id_r)
37 d_dict = StaticDict(id_l => d_l, id_r => d_r)
38
39 H_l = inner_product(boundary_grid(g_1D,id_l), quadrature_interior, quadrature_closure)
40 H_r = inner_product(boundary_grid(g_1D,id_r), quadrature_interior, quadrature_closure)
41 Hb_dict = StaticDict(id_l => H_l, id_r => H_r)
42
43 L = Laplace(g_1D, operator_path; order=4)
44 @test L == Laplace(Δ, H, Hi, e_dict, d_dict, Hb_dict)
45 @test L isa TensorMapping{T,1,1} where T
46 @inferred Laplace(Δ, H, Hi, e_dict, d_dict, Hb_dict)
47 # REVIEW: The tests above seem very tied to the implementation. Is
48 # it important that the components of the operator set are stored
49 # in static dicts? Is something like below better?
50 #
51 # ```
52 # L = Laplace(g_1D, operator_path; order=4)
53 # @test L isa TensorMapping{T,1,1} where T
54 # @test boundary_restriction(L,id_l) == boundary_restriction(g_1D, e_closure,id_l)
55 # ...
56 # ```
57 # I guess this is more or less simply a reorganization of the test and skipping testing for the struct layout
58 end 23 end
59 @testset "3D" begin 24 @testset "3D" begin
60 Δ = laplace(g_3D, inner_stencil, closure_stencils) 25 Δ = laplace(g_3D, inner_stencil, closure_stencils)
61 H = inner_product(g_3D, quadrature_interior, quadrature_closure) 26 @test Laplace(g_3D, stencil_set) == Laplace(Δ,stencil_set)
62 Hi = inverse_inner_product(g_3D, quadrature_interior, quadrature_closure) 27 @test Laplace(g_3D, stencil_set) isa TensorMapping{T,3,3} where T
63
64 (id_l, id_r, id_s, id_n, id_b, id_t) = boundary_identifiers(g_3D)
65 e_l = boundary_restriction(g_3D, e_closure,id_l)
66 e_r = boundary_restriction(g_3D, e_closure,id_r)
67 e_s = boundary_restriction(g_3D, e_closure,id_s)
68 e_n = boundary_restriction(g_3D, e_closure,id_n)
69 e_b = boundary_restriction(g_3D, e_closure,id_b)
70 e_t = boundary_restriction(g_3D, e_closure,id_t)
71 e_dict = StaticDict(id_l => e_l, id_r => e_r,
72 id_s => e_s, id_n => e_n,
73 id_b => e_b, id_t => e_t)
74
75 d_l = normal_derivative(g_3D, d_closure,id_l)
76 d_r = normal_derivative(g_3D, d_closure,id_r)
77 d_s = normal_derivative(g_3D, d_closure,id_s)
78 d_n = normal_derivative(g_3D, d_closure,id_n)
79 d_b = normal_derivative(g_3D, d_closure,id_b)
80 d_t = normal_derivative(g_3D, d_closure,id_t)
81 d_dict = StaticDict(id_l => d_l, id_r => d_r,
82 id_s => d_s, id_n => d_n,
83 id_b => d_b, id_t => d_t)
84
85 H_l = inner_product(boundary_grid(g_3D,id_l), quadrature_interior, quadrature_closure)
86 H_r = inner_product(boundary_grid(g_3D,id_r), quadrature_interior, quadrature_closure)
87 H_s = inner_product(boundary_grid(g_3D,id_s), quadrature_interior, quadrature_closure)
88 H_n = inner_product(boundary_grid(g_3D,id_n), quadrature_interior, quadrature_closure)
89 H_b = inner_product(boundary_grid(g_3D,id_b), quadrature_interior, quadrature_closure)
90 H_t = inner_product(boundary_grid(g_3D,id_t), quadrature_interior, quadrature_closure)
91 Hb_dict = StaticDict(id_l => H_l, id_r => H_r,
92 id_s => H_s, id_n => H_n,
93 id_b => H_b, id_t => H_t)
94
95 L = Laplace(g_3D, operator_path; order=4)
96 @test L == Laplace(Δ,H,Hi,e_dict,d_dict,Hb_dict)
97 @test L isa TensorMapping{T,3,3} where T
98 @inferred Laplace(Δ,H,Hi,e_dict,d_dict,Hb_dict)
99 end 28 end
100 end
101
102 # REVIEW: Is this testset misplaced? Should it really be inside the "Laplace" testset?
103 @testset "laplace" begin
104 @testset "1D" begin
105 L = laplace(g_1D, inner_stencil, closure_stencils)
106 @test L == second_derivative(g_1D, inner_stencil, closure_stencils)
107 @test L isa TensorMapping{T,1,1} where T
108 end
109 @testset "3D" begin
110 L = laplace(g_3D, inner_stencil, closure_stencils)
111 @test L isa TensorMapping{T,3,3} where T
112 Dxx = second_derivative(g_3D, inner_stencil, closure_stencils, 1)
113 Dyy = second_derivative(g_3D, inner_stencil, closure_stencils, 2)
114 Dzz = second_derivative(g_3D, inner_stencil, closure_stencils, 3)
115 @test L == Dxx + Dyy + Dzz
116 @test L isa TensorMapping{T,3,3} where T
117 end
118 end
119
120 @testset "inner_product" begin
121 L = Laplace(g_3D, operator_path; order=4)
122 @test inner_product(L) == inner_product(g_3D, quadrature_interior, quadrature_closure)
123 end
124
125 @testset "inverse_inner_product" begin
126 L = Laplace(g_3D, operator_path; order=4)
127 @test inverse_inner_product(L) == inverse_inner_product(g_3D, quadrature_interior, quadrature_closure)
128 end
129
130 @testset "boundary_restriction" begin
131 L = Laplace(g_3D, operator_path; order=4)
132 (id_l, id_r, id_s, id_n, id_b, id_t) = boundary_identifiers(g_3D)
133 @test boundary_restriction(L, id_l) == boundary_restriction(g_3D, e_closure,id_l)
134 @test boundary_restriction(L, id_r) == boundary_restriction(g_3D, e_closure,id_r)
135 @test boundary_restriction(L, id_s) == boundary_restriction(g_3D, e_closure,id_s)
136 @test boundary_restriction(L, id_n) == boundary_restriction(g_3D, e_closure,id_n)
137 @test boundary_restriction(L, id_b) == boundary_restriction(g_3D, e_closure,id_b)
138 @test boundary_restriction(L, id_t) == boundary_restriction(g_3D, e_closure,id_t)
139
140 ids = boundary_identifiers(g_3D)
141 es = boundary_restriction(L, ids)
142 @test es == (boundary_restriction(L, id_l),
143 boundary_restriction(L, id_r),
144 boundary_restriction(L, id_s),
145 boundary_restriction(L, id_n),
146 boundary_restriction(L, id_b),
147 boundary_restriction(L, id_t));
148 @test es == boundary_restriction(L, ids...)
149 end
150
151 @testset "normal_derivative" begin
152 L = Laplace(g_3D, operator_path; order=4)
153 (id_l, id_r, id_s, id_n, id_b, id_t) = boundary_identifiers(g_3D)
154 @test normal_derivative(L, id_l) == normal_derivative(g_3D, d_closure,id_l)
155 @test normal_derivative(L, id_r) == normal_derivative(g_3D, d_closure,id_r)
156 @test normal_derivative(L, id_s) == normal_derivative(g_3D, d_closure,id_s)
157 @test normal_derivative(L, id_n) == normal_derivative(g_3D, d_closure,id_n)
158 @test normal_derivative(L, id_b) == normal_derivative(g_3D, d_closure,id_b)
159 @test normal_derivative(L, id_t) == normal_derivative(g_3D, d_closure,id_t)
160
161 ids = boundary_identifiers(g_3D)
162 ds = normal_derivative(L, ids)
163 @test ds == (normal_derivative(L, id_l),
164 normal_derivative(L, id_r),
165 normal_derivative(L, id_s),
166 normal_derivative(L, id_n),
167 normal_derivative(L, id_b),
168 normal_derivative(L, id_t));
169 @test ds == normal_derivative(L, ids...)
170 end
171
172 @testset "boundary_quadrature" begin
173 L = Laplace(g_3D, operator_path; order=4)
174 (id_l, id_r, id_s, id_n, id_b, id_t) = boundary_identifiers(g_3D)
175 @test boundary_quadrature(L, id_l) == inner_product(boundary_grid(g_3D, id_l), quadrature_interior, quadrature_closure)
176 @test boundary_quadrature(L, id_r) == inner_product(boundary_grid(g_3D, id_r), quadrature_interior, quadrature_closure)
177 @test boundary_quadrature(L, id_s) == inner_product(boundary_grid(g_3D, id_s), quadrature_interior, quadrature_closure)
178 @test boundary_quadrature(L, id_n) == inner_product(boundary_grid(g_3D, id_n), quadrature_interior, quadrature_closure)
179 @test boundary_quadrature(L, id_b) == inner_product(boundary_grid(g_3D, id_b), quadrature_interior, quadrature_closure)
180 @test boundary_quadrature(L, id_t) == inner_product(boundary_grid(g_3D, id_t), quadrature_interior, quadrature_closure)
181
182 ids = boundary_identifiers(g_3D)
183 H_gammas = boundary_quadrature(L, ids)
184 @test H_gammas == (boundary_quadrature(L, id_l),
185 boundary_quadrature(L, id_r),
186 boundary_quadrature(L, id_s),
187 boundary_quadrature(L, id_n),
188 boundary_quadrature(L, id_b),
189 boundary_quadrature(L, id_t));
190 @test H_gammas == boundary_quadrature(L, ids...)
191 end 29 end
192 30
193 # Exact differentiation is measured point-wise. In other cases 31 # Exact differentiation is measured point-wise. In other cases
194 # the error is measured in the l2-norm. 32 # the error is measured in the l2-norm.
195 @testset "Accuracy" begin 33 @testset "Accuracy" begin
205 43
206 # 2nd order interior stencil, 1st order boundary stencil, 44 # 2nd order interior stencil, 1st order boundary stencil,
207 # implies that L*v should be exact for binomials up to order 2. 45 # implies that L*v should be exact for binomials up to order 2.
208 @testset "2nd order" begin 46 @testset "2nd order" begin
209 stencil_set = read_stencil_set(operator_path; order=2) 47 stencil_set = read_stencil_set(operator_path; order=2)
210 inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"]) 48 Δ = Laplace(g_3D, stencil_set)
211 closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"]) 49 @test Δ*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
212 L = laplace(g_3D, inner_stencil, closure_stencils) 50 @test Δ*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
213 @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9 51 @test Δ*polynomials[3] ≈ polynomials[1] atol = 5e-9
214 @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9 52 @test Δ*v ≈ Δv rtol = 5e-2 norm = l2
215 @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9
216 @test L*v ≈ Δv rtol = 5e-2 norm = l2
217 end 53 end
218 54
219 # 4th order interior stencil, 2nd order boundary stencil, 55 # 4th order interior stencil, 2nd order boundary stencil,
220 # implies that L*v should be exact for binomials up to order 3. 56 # implies that L*v should be exact for binomials up to order 3.
221 @testset "4th order" begin 57 @testset "4th order" begin
222 stencil_set = read_stencil_set(operator_path; order=4) 58 stencil_set = read_stencil_set(operator_path; order=4)
223 inner_stencil = parse_stencil(stencil_set["D2"]["inner_stencil"]) 59 Δ = Laplace(g_3D, stencil_set)
224 closure_stencils = parse_stencil.(stencil_set["D2"]["closure_stencils"])
225 L = laplace(g_3D, inner_stencil, closure_stencils)
226 # NOTE: high tolerances for checking the "exact" differentiation 60 # NOTE: high tolerances for checking the "exact" differentiation
227 # due to accumulation of round-off errors/cancellation errors? 61 # due to accumulation of round-off errors/cancellation errors?
228 @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9 62 @test Δ*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
229 @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9 63 @test Δ*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
230 @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9 64 @test Δ*polynomials[3] ≈ polynomials[1] atol = 5e-9
231 @test L*polynomials[4] ≈ polynomials[2] atol = 5e-9 65 @test Δ*polynomials[4] ≈ polynomials[2] atol = 5e-9
232 @test L*v ≈ Δv rtol = 5e-4 norm = l2 66 @test Δ*v ≈ Δv rtol = 5e-4 norm = l2
233 end 67 end
234 end 68 end
235 end 69 end
70
71 @testset "laplace" begin
72 @testset "1D" begin
73 Δ = laplace(g_1D, inner_stencil, closure_stencils)
74 @test Δ == second_derivative(g_1D, inner_stencil, closure_stencils)
75 @test Δ isa TensorMapping{T,1,1} where T
76 end
77 @testset "3D" begin
78 Δ = laplace(g_3D, inner_stencil, closure_stencils)
79 @test Δ isa TensorMapping{T,3,3} where T
80 Dxx = second_derivative(g_3D, inner_stencil, closure_stencils, 1)
81 Dyy = second_derivative(g_3D, inner_stencil, closure_stencils, 2)
82 Dzz = second_derivative(g_3D, inner_stencil, closure_stencils, 3)
83 @test Δ == Dxx + Dyy + Dzz
84 @test Δ isa TensorMapping{T,3,3} where T
85 end
86 end
87