Mercurial > repos > public > sbplib_julia

comparison test/testSbpOperators.jl @ 333:01b851161018 refactor/combine_to_one_package

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Start converting to one package by moving all the files to their correct location
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 25 Sep 2020 13:06:02 +0200
parents SbpOperators/test/runtests.jl@408c37b295c2
children f4e3e71a4ff4
comparison
equal deleted inserted replaced
332:535f1bff4bcc 333:01b851161018
1 using Test
2 using SbpOperators
3 using Grids
4 using RegionIndices
5 using LazyTensors
6
7 # @testset "apply_quadrature" begin
8 # op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
9 # h = 0.5
10 #
11 # @test apply_quadrature(op, h, 1.0, 10, 100) == h
12 #
13 # N = 10
14 # qc = op.quadratureClosure
15 # q = h.*(qc..., ones(N-2*closuresize(op))..., reverse(qc)...)
16 # @assert length(q) == N
17 #
18 # for i ∈ 1:N
19 # @test apply_quadrature(op, h, 1.0, i, N) == q[i]
20 # end
21 #
22 # v = [2.,3.,2.,4.,5.,4.,3.,4.,5.,4.5]
23 # for i ∈ 1:N
24 # @test apply_quadrature(op, h, v[i], i, N) == q[i]*v[i]
25 # end
26 # end
27
28 @testset "SecondDerivative" begin
29 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
30 L = 3.5
31 g = EquidistantGrid((101,), (0.0,), (L,))
32 h_inv = inverse_spacing(g)
33 h = 1/h_inv[1];
34 Dₓₓ = SecondDerivative(h_inv[1],op.innerStencil,op.closureStencils,op.parity)
35
36 f0(x::Float64) = 1.
37 f1(x::Float64) = x
38 f2(x::Float64) = 1/2*x^2
39 f3(x::Float64) = 1/6*x^3
40 f4(x::Float64) = 1/24*x^4
41 f5(x::Float64) = sin(x)
42 f5ₓₓ(x::Float64) = -f5(x)
43
44 v0 = evalOn(g,f0)
45 v1 = evalOn(g,f1)
46 v2 = evalOn(g,f2)
47 v3 = evalOn(g,f3)
48 v4 = evalOn(g,f4)
49 v5 = evalOn(g,f5)
50
51 @test Dₓₓ isa TensorOperator{T,1} where T
52 @test Dₓₓ' isa TensorMapping{T,1,1} where T
53
54 # TODO: Should perhaps set tolerance level for isapporx instead?
55 # Are these tolerance levels resonable or should tests be constructed
56 # differently?
57 equalitytol = 0.5*1e-10
58 accuracytol = 0.5*1e-3
59 # 4th order interior stencil, 2nd order boundary stencil,
60 # implies that L*v should be exact for v - monomial up to order 3.
61 # Exact differentiation is measured point-wise. For other grid functions
62 # the error is measured in the l2-norm.
63 @test all(abs.(collect(Dₓₓ*v0)) .<= equalitytol)
64 @test all(abs.(collect(Dₓₓ*v1)) .<= equalitytol)
65 @test all(abs.((collect(Dₓₓ*v2) - v0)) .<= equalitytol)
66 @test all(abs.((collect(Dₓₓ*v3) - v1)) .<= equalitytol)
67 e4 = collect(Dₓₓ*v4) - v2
68 e5 = collect(Dₓₓ*v5) + v5
69 @test sqrt(h*sum(collect(e4.^2))) <= accuracytol
70 @test sqrt(h*sum(collect(e5.^2))) <= accuracytol
71 end
72
73 @testset "Laplace2D" begin
74 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
75 Lx = 1.5
76 Ly = 3.2
77 g = EquidistantGrid((102,131), (0.0, 0.0), (Lx,Ly))
78
79 h_inv = inverse_spacing(g)
80 h = spacing(g)
81 D_xx = SecondDerivative(h_inv[1],op.innerStencil,op.closureStencils,op.parity)
82 D_yy = SecondDerivative(h_inv[2],op.innerStencil,op.closureStencils,op.parity)
83 L = Laplace((D_xx,D_yy))
84
85 f0(x::Float64,y::Float64) = 2.
86 f1(x::Float64,y::Float64) = x+y
87 f2(x::Float64,y::Float64) = 1/2*x^2 + 1/2*y^2
88 f3(x::Float64,y::Float64) = 1/6*x^3 + 1/6*y^3
89 f4(x::Float64,y::Float64) = 1/24*x^4 + 1/24*y^4
90 f5(x::Float64,y::Float64) = sin(x) + cos(y)
91 f5ₓₓ(x::Float64,y::Float64) = -f5(x,y)
92
93 v0 = evalOn(g,f0)
94 v1 = evalOn(g,f1)
95 v2 = evalOn(g,f2)
96 v3 = evalOn(g,f3)
97 v4 = evalOn(g,f4)
98 v5 = evalOn(g,f5)
99 v5ₓₓ = evalOn(g,f5ₓₓ)
100
101 @test L isa TensorOperator{T,2} where T
102 @test L' isa TensorMapping{T,2,2} where T
103
104 # TODO: Should perhaps set tolerance level for isapporx instead?
105 # Are these tolerance levels resonable or should tests be constructed
106 # differently?
107 equalitytol = 0.5*1e-10
108 accuracytol = 0.5*1e-3
109 # 4th order interior stencil, 2nd order boundary stencil,
110 # implies that L*v should be exact for v - monomial up to order 3.
111 # Exact differentiation is measured point-wise. For other grid functions
112 # the error is measured in the H-norm.
113 @test all(abs.(collect(L*v0)) .<= equalitytol)
114 @test all(abs.(collect(L*v1)) .<= equalitytol)
115 @test all(collect(L*v2) .≈ v0) # Seems to be more accurate
116 @test all(abs.((collect(L*v3) - v1)) .<= equalitytol)
117 e4 = collect(L*v4) - v2
118 e5 = collect(L*v5) - v5ₓₓ
119 @test sqrt(prod(h)*sum(collect(e4.^2))) <= accuracytol
120 @test sqrt(prod(h)*sum(collect(e5.^2))) <= accuracytol
121 end
122
123 @testset "DiagonalInnerProduct" begin
124 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
125 L = 2.3
126 g = EquidistantGrid((77,), (0.0,), (L,))
127 h = spacing(g)
128 H = DiagonalInnerProduct(h[1],op.quadratureClosure)
129 v = ones(Float64, size(g))
130
131 @test H isa TensorOperator{T,1} where T
132 @test H' isa TensorMapping{T,1,1} where T
133 @test sum(collect(H*v)) ≈ L
134 @test collect(H*v) == collect(H'*v)
135 end
136
137 @testset "Quadrature" begin
138 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
139 Lx = 2.3
140 Ly = 5.2
141 g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
142
143 h = spacing(g)
144 Hx = DiagonalInnerProduct(h[1],op.quadratureClosure);
145 Hy = DiagonalInnerProduct(h[2],op.quadratureClosure);
146 Q = Quadrature((Hx,Hy))
147
148 v = ones(Float64, size(g))
149
150 @test Q isa TensorOperator{T,2} where T
151 @test Q' isa TensorMapping{T,2,2} where T
152 @test sum(collect(Q*v)) ≈ (Lx*Ly)
153 @test collect(Q*v) == collect(Q'*v)
154 end
155
156 @testset "InverseDiagonalInnerProduct" begin
157 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
158 L = 2.3
159 g = EquidistantGrid((77,), (0.0,), (L,))
160 h = spacing(g)
161 H = DiagonalInnerProduct(h[1],op.quadratureClosure)
162
163 h_i = inverse_spacing(g)
164 Hi = InverseDiagonalInnerProduct(h_i[1],1 ./ op.quadratureClosure)
165 v = evalOn(g, x->sin(x))
166
167 @test Hi isa TensorOperator{T,1} where T
168 @test Hi' isa TensorMapping{T,1,1} where T
169 @test collect(Hi*H*v) ≈ v
170 @test collect(Hi*v) == collect(Hi'*v)
171 end
172
173 @testset "InverseQuadrature" begin
174 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
175 Lx = 7.3
176 Ly = 8.2
177 g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
178
179 h = spacing(g)
180 Hx = DiagonalInnerProduct(h[1], op.quadratureClosure);
181 Hy = DiagonalInnerProduct(h[2], op.quadratureClosure);
182 Q = Quadrature((Hx,Hy))
183
184 hi = inverse_spacing(g)
185 Hix = InverseDiagonalInnerProduct(hi[1], 1 ./ op.quadratureClosure);
186 Hiy = InverseDiagonalInnerProduct(hi[2], 1 ./ op.quadratureClosure);
187 Qinv = InverseQuadrature((Hix,Hiy))
188 v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
189
190 @test Qinv isa TensorOperator{T,2} where T
191 @test Qinv' isa TensorMapping{T,2,2} where T
192 @test collect(Qinv*Q*v) ≈ v
193 @test collect(Qinv*v) == collect(Qinv'*v)
194 end
195 #
196 # @testset "BoundaryValue" begin
197 # op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
198 # g = EquidistantGrid((4,5), (0.0, 0.0), (1.0,1.0))
199 #
200 # e_w = BoundaryValue(op, g, CartesianBoundary{1,Lower}())
201 # e_e = BoundaryValue(op, g, CartesianBoundary{1,Upper}())
202 # e_s = BoundaryValue(op, g, CartesianBoundary{2,Lower}())
203 # e_n = BoundaryValue(op, g, CartesianBoundary{2,Upper}())
204 #
205 # v = zeros(Float64, 4, 5)
206 # v[:,5] = [1, 2, 3,4]
207 # v[:,4] = [1, 2, 3,4]
208 # v[:,3] = [4, 5, 6, 7]
209 # v[:,2] = [7, 8, 9, 10]
210 # v[:,1] = [10, 11, 12, 13]
211 #
212 # @test e_w isa TensorMapping{T,2,1} where T
213 # @test e_w' isa TensorMapping{T,1,2} where T
214 #
215 # @test domain_size(e_w, (3,2)) == (2,)
216 # @test domain_size(e_e, (3,2)) == (2,)
217 # @test domain_size(e_s, (3,2)) == (3,)
218 # @test domain_size(e_n, (3,2)) == (3,)
219 #
220 # @test size(e_w'*v) == (5,)
221 # @test size(e_e'*v) == (5,)
222 # @test size(e_s'*v) == (4,)
223 # @test size(e_n'*v) == (4,)
224 #
225 # @test collect(e_w'*v) == [10,7,4,1.0,1]
226 # @test collect(e_e'*v) == [13,10,7,4,4.0]
227 # @test collect(e_s'*v) == [10,11,12,13.0]
228 # @test collect(e_n'*v) == [1,2,3,4.0]
229 #
230 # g_x = [1,2,3,4.0]
231 # g_y = [5,4,3,2,1.0]
232 #
233 # G_w = zeros(Float64, (4,5))
234 # G_w[1,:] = g_y
235 #
236 # G_e = zeros(Float64, (4,5))
237 # G_e[4,:] = g_y
238 #
239 # G_s = zeros(Float64, (4,5))
240 # G_s[:,1] = g_x
241 #
242 # G_n = zeros(Float64, (4,5))
243 # G_n[:,5] = g_x
244 #
245 # @test size(e_w*g_y) == (UnknownDim,5)
246 # @test size(e_e*g_y) == (UnknownDim,5)
247 # @test size(e_s*g_x) == (4,UnknownDim)
248 # @test size(e_n*g_x) == (4,UnknownDim)
249 #
250 # # These tests should be moved to where they are possible (i.e we know what the grid should be)
251 # @test_broken collect(e_w*g_y) == G_w
252 # @test_broken collect(e_e*g_y) == G_e
253 # @test_broken collect(e_s*g_x) == G_s
254 # @test_broken collect(e_n*g_x) == G_n
255 # end
256 #
257 # @testset "NormalDerivative" begin
258 # op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
259 # g = EquidistantGrid((5,6), (0.0, 0.0), (4.0,5.0))
260 #
261 # d_w = NormalDerivative(op, g, CartesianBoundary{1,Lower}())
262 # d_e = NormalDerivative(op, g, CartesianBoundary{1,Upper}())
263 # d_s = NormalDerivative(op, g, CartesianBoundary{2,Lower}())
264 # d_n = NormalDerivative(op, g, CartesianBoundary{2,Upper}())
265 #
266 #
267 # v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
268 # v∂x = evalOn(g, (x,y)-> 2*x + y)
269 # v∂y = evalOn(g, (x,y)-> 2*(y-1) + x)
270 #
271 # @test d_w isa TensorMapping{T,2,1} where T
272 # @test d_w' isa TensorMapping{T,1,2} where T
273 #
274 # @test domain_size(d_w, (3,2)) == (2,)
275 # @test domain_size(d_e, (3,2)) == (2,)
276 # @test domain_size(d_s, (3,2)) == (3,)
277 # @test domain_size(d_n, (3,2)) == (3,)
278 #
279 # @test size(d_w'*v) == (6,)
280 # @test size(d_e'*v) == (6,)
281 # @test size(d_s'*v) == (5,)
282 # @test size(d_n'*v) == (5,)
283 #
284 # @test collect(d_w'*v) ≈ v∂x[1,:]
285 # @test collect(d_e'*v) ≈ v∂x[5,:]
286 # @test collect(d_s'*v) ≈ v∂y[:,1]
287 # @test collect(d_n'*v) ≈ v∂y[:,6]
288 #
289 #
290 # d_x_l = zeros(Float64, 5)
291 # d_x_u = zeros(Float64, 5)
292 # for i ∈ eachindex(d_x_l)
293 # d_x_l[i] = op.dClosure[i-1]
294 # d_x_u[i] = -op.dClosure[length(d_x_u)-i]
295 # end
296 #
297 # d_y_l = zeros(Float64, 6)
298 # d_y_u = zeros(Float64, 6)
299 # for i ∈ eachindex(d_y_l)
300 # d_y_l[i] = op.dClosure[i-1]
301 # d_y_u[i] = -op.dClosure[length(d_y_u)-i]
302 # end
303 #
304 # function prod_matrix(x,y)
305 # G = zeros(Float64, length(x), length(y))
306 # for I ∈ CartesianIndices(G)
307 # G[I] = x[I[1]]*y[I[2]]
308 # end
309 #
310 # return G
311 # end
312 #
313 # g_x = [1,2,3,4.0,5]
314 # g_y = [5,4,3,2,1.0,11]
315 #
316 # G_w = prod_matrix(d_x_l, g_y)
317 # G_e = prod_matrix(d_x_u, g_y)
318 # G_s = prod_matrix(g_x, d_y_l)
319 # G_n = prod_matrix(g_x, d_y_u)
320 #
321 #
322 # @test size(d_w*g_y) == (UnknownDim,6)
323 # @test size(d_e*g_y) == (UnknownDim,6)
324 # @test size(d_s*g_x) == (5,UnknownDim)
325 # @test size(d_n*g_x) == (5,UnknownDim)
326 #
327 # # These tests should be moved to where they are possible (i.e we know what the grid should be)
328 # @test_broken collect(d_w*g_y) ≈ G_w
329 # @test_broken collect(d_e*g_y) ≈ G_e
330 # @test_broken collect(d_s*g_x) ≈ G_s
331 # @test_broken collect(d_n*g_x) ≈ G_n
332 # end
333 #
334 # @testset "BoundaryQuadrature" begin
335 # op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
336 # g = EquidistantGrid((10,11), (0.0, 0.0), (1.0,1.0))
337 #
338 # H_w = BoundaryQuadrature(op, g, CartesianBoundary{1,Lower}())
339 # H_e = BoundaryQuadrature(op, g, CartesianBoundary{1,Upper}())
340 # H_s = BoundaryQuadrature(op, g, CartesianBoundary{2,Lower}())
341 # H_n = BoundaryQuadrature(op, g, CartesianBoundary{2,Upper}())
342 #
343 # v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
344 #
345 # function get_quadrature(N)
346 # qc = op.quadratureClosure
347 # q = (qc..., ones(N-2*closuresize(op))..., reverse(qc)...)
348 # @assert length(q) == N
349 # return q
350 # end
351 #
352 # v_w = v[1,:]
353 # v_e = v[10,:]
354 # v_s = v[:,1]
355 # v_n = v[:,11]
356 #
357 # q_x = spacing(g)[1].*get_quadrature(10)
358 # q_y = spacing(g)[2].*get_quadrature(11)
359 #
360 # @test H_w isa TensorOperator{T,1} where T
361 #
362 # @test domain_size(H_w, (3,)) == (3,)
363 # @test domain_size(H_n, (3,)) == (3,)
364 #
365 # @test range_size(H_w, (3,)) == (3,)
366 # @test range_size(H_n, (3,)) == (3,)
367 #
368 # @test size(H_w*v_w) == (11,)
369 # @test size(H_e*v_e) == (11,)
370 # @test size(H_s*v_s) == (10,)
371 # @test size(H_n*v_n) == (10,)
372 #
373 # @test collect(H_w*v_w) ≈ q_y.*v_w
374 # @test collect(H_e*v_e) ≈ q_y.*v_e
375 # @test collect(H_s*v_s) ≈ q_x.*v_s
376 # @test collect(H_n*v_n) ≈ q_x.*v_n
377 #
378 # @test collect(H_w'*v_w) == collect(H_w'*v_w)
379 # @test collect(H_e'*v_e) == collect(H_e'*v_e)
380 # @test collect(H_s'*v_s) == collect(H_s'*v_s)
381 # @test collect(H_n'*v_n) == collect(H_n'*v_n)
382 # end