Mercurial > repos > public > sbplib_julia

comparison SbpOperators/test/runtests.jl @ 314:accb0876da12

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Move tests for Laplace from DiffOps/test to SbpOperators/test and add test for Secondderivative
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 09 Sep 2020 21:21:35 +0200
parents 4ca3794fffef
children 9cc5d1498b2d
comparison
equal deleted inserted replaced
313:d1004b881da1 314:accb0876da12
1 using Test
1 using SbpOperators 2 using SbpOperators
2 using Test 3 using Grids
3 4 using RegionIndices
4 @testset "apply_quadrature" begin 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
5 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 29 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
6 h = 0.5 30 L = 3.5
7 31 g = EquidistantGrid((101,), (0.0,), (L,))
8 @test apply_quadrature(op, h, 1.0, 10, 100) == h 32 h_inv = inverse_spacing(g)
9 33 h = 1/h_inv[1];
10 N = 10 34 Dₓₓ = SecondDerivative(h_inv[1],op.innerStencil,op.closureStencils,op.parity)
11 qc = op.quadratureClosure 35
12 q = h.*(qc..., ones(N-2*closuresize(op))..., reverse(qc)...) 36 f0(x::Float64) = 1.
13 @assert length(q) == N 37 f1(x::Float64) = x
14 38 f2(x::Float64) = 1/2*x^2
15 for i ∈ 1:N 39 f3(x::Float64) = 1/6*x^3
16 @test apply_quadrature(op, h, 1.0, i, N) == q[i] 40 f4(x::Float64) = 1/24*x^4
17 end 41 f5(x::Float64) = sin(x)
18 42 f5ₓₓ(x::Float64) = -f5(x)
19 v = [2.,3.,2.,4.,5.,4.,3.,4.,5.,4.5] 43
20 for i ∈ 1:N 44 v0 = evalOn(g,f0)
21 @test apply_quadrature(op, h, v[i], i, N) == q[i]*v[i] 45 v1 = evalOn(g,f1)
22 end 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
23 end 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 "Quadrature" begin
124 # op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
125 # Lx = 2.3
126 # Ly = 5.2
127 # g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
128 # H = Quadrature(op,g)
129 # v = ones(Float64, size(g))
130 #
131 # @test H isa TensorOperator{T,2} where T
132 # @test H' isa TensorMapping{T,2,2} where T
133 # @test sum(collect(H*v)) ≈ (Lx*Ly)
134 # @test collect(H*v) == collect(H'*v)
135 # end
136 #
137 # @testset "InverseQuadrature" begin
138 # op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
139 # Lx = 7.3
140 # Ly = 8.2
141 # g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
142 # H = Quadrature(op,g)
143 # Hinv = InverseQuadrature(op,g)
144 # v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
145 #
146 # @test Hinv isa TensorOperator{T,2} where T
147 # @test Hinv' isa TensorMapping{T,2,2} where T
148 # @test collect(Hinv*H*v) ≈ v
149 # @test collect(Hinv*v) == collect(Hinv'*v)
150 # end
151 #
152 # @testset "BoundaryValue" begin
153 # op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
154 # g = EquidistantGrid((4,5), (0.0, 0.0), (1.0,1.0))
155 #
156 # e_w = BoundaryValue(op, g, CartesianBoundary{1,Lower}())
157 # e_e = BoundaryValue(op, g, CartesianBoundary{1,Upper}())
158 # e_s = BoundaryValue(op, g, CartesianBoundary{2,Lower}())
159 # e_n = BoundaryValue(op, g, CartesianBoundary{2,Upper}())
160 #
161 # v = zeros(Float64, 4, 5)
162 # v[:,5] = [1, 2, 3,4]
163 # v[:,4] = [1, 2, 3,4]
164 # v[:,3] = [4, 5, 6, 7]
165 # v[:,2] = [7, 8, 9, 10]
166 # v[:,1] = [10, 11, 12, 13]
167 #
168 # @test e_w isa TensorMapping{T,2,1} where T
169 # @test e_w' isa TensorMapping{T,1,2} where T
170 #
171 # @test domain_size(e_w, (3,2)) == (2,)
172 # @test domain_size(e_e, (3,2)) == (2,)
173 # @test domain_size(e_s, (3,2)) == (3,)
174 # @test domain_size(e_n, (3,2)) == (3,)
175 #
176 # @test size(e_w'*v) == (5,)
177 # @test size(e_e'*v) == (5,)
178 # @test size(e_s'*v) == (4,)
179 # @test size(e_n'*v) == (4,)
180 #
181 # @test collect(e_w'*v) == [10,7,4,1.0,1]
182 # @test collect(e_e'*v) == [13,10,7,4,4.0]
183 # @test collect(e_s'*v) == [10,11,12,13.0]
184 # @test collect(e_n'*v) == [1,2,3,4.0]
185 #
186 # g_x = [1,2,3,4.0]
187 # g_y = [5,4,3,2,1.0]
188 #
189 # G_w = zeros(Float64, (4,5))
190 # G_w[1,:] = g_y
191 #
192 # G_e = zeros(Float64, (4,5))
193 # G_e[4,:] = g_y
194 #
195 # G_s = zeros(Float64, (4,5))
196 # G_s[:,1] = g_x
197 #
198 # G_n = zeros(Float64, (4,5))
199 # G_n[:,5] = g_x
200 #
201 # @test size(e_w*g_y) == (UnknownDim,5)
202 # @test size(e_e*g_y) == (UnknownDim,5)
203 # @test size(e_s*g_x) == (4,UnknownDim)
204 # @test size(e_n*g_x) == (4,UnknownDim)
205 #
206 # # These tests should be moved to where they are possible (i.e we know what the grid should be)
207 # @test_broken collect(e_w*g_y) == G_w
208 # @test_broken collect(e_e*g_y) == G_e
209 # @test_broken collect(e_s*g_x) == G_s
210 # @test_broken collect(e_n*g_x) == G_n
211 # end
212 #
213 # @testset "NormalDerivative" begin
214 # op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
215 # g = EquidistantGrid((5,6), (0.0, 0.0), (4.0,5.0))
216 #
217 # d_w = NormalDerivative(op, g, CartesianBoundary{1,Lower}())
218 # d_e = NormalDerivative(op, g, CartesianBoundary{1,Upper}())
219 # d_s = NormalDerivative(op, g, CartesianBoundary{2,Lower}())
220 # d_n = NormalDerivative(op, g, CartesianBoundary{2,Upper}())
221 #
222 #
223 # v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
224 # v∂x = evalOn(g, (x,y)-> 2*x + y)
225 # v∂y = evalOn(g, (x,y)-> 2*(y-1) + x)
226 #
227 # @test d_w isa TensorMapping{T,2,1} where T
228 # @test d_w' isa TensorMapping{T,1,2} where T
229 #
230 # @test domain_size(d_w, (3,2)) == (2,)
231 # @test domain_size(d_e, (3,2)) == (2,)
232 # @test domain_size(d_s, (3,2)) == (3,)
233 # @test domain_size(d_n, (3,2)) == (3,)
234 #
235 # @test size(d_w'*v) == (6,)
236 # @test size(d_e'*v) == (6,)
237 # @test size(d_s'*v) == (5,)
238 # @test size(d_n'*v) == (5,)
239 #
240 # @test collect(d_w'*v) ≈ v∂x[1,:]
241 # @test collect(d_e'*v) ≈ v∂x[5,:]
242 # @test collect(d_s'*v) ≈ v∂y[:,1]
243 # @test collect(d_n'*v) ≈ v∂y[:,6]
244 #
245 #
246 # d_x_l = zeros(Float64, 5)
247 # d_x_u = zeros(Float64, 5)
248 # for i ∈ eachindex(d_x_l)
249 # d_x_l[i] = op.dClosure[i-1]
250 # d_x_u[i] = -op.dClosure[length(d_x_u)-i]
251 # end
252 #
253 # d_y_l = zeros(Float64, 6)
254 # d_y_u = zeros(Float64, 6)
255 # for i ∈ eachindex(d_y_l)
256 # d_y_l[i] = op.dClosure[i-1]
257 # d_y_u[i] = -op.dClosure[length(d_y_u)-i]
258 # end
259 #
260 # function prod_matrix(x,y)
261 # G = zeros(Float64, length(x), length(y))
262 # for I ∈ CartesianIndices(G)
263 # G[I] = x[I[1]]*y[I[2]]
264 # end
265 #
266 # return G
267 # end
268 #
269 # g_x = [1,2,3,4.0,5]
270 # g_y = [5,4,3,2,1.0,11]
271 #
272 # G_w = prod_matrix(d_x_l, g_y)
273 # G_e = prod_matrix(d_x_u, g_y)
274 # G_s = prod_matrix(g_x, d_y_l)
275 # G_n = prod_matrix(g_x, d_y_u)
276 #
277 #
278 # @test size(d_w*g_y) == (UnknownDim,6)
279 # @test size(d_e*g_y) == (UnknownDim,6)
280 # @test size(d_s*g_x) == (5,UnknownDim)
281 # @test size(d_n*g_x) == (5,UnknownDim)
282 #
283 # # These tests should be moved to where they are possible (i.e we know what the grid should be)
284 # @test_broken collect(d_w*g_y) ≈ G_w
285 # @test_broken collect(d_e*g_y) ≈ G_e
286 # @test_broken collect(d_s*g_x) ≈ G_s
287 # @test_broken collect(d_n*g_x) ≈ G_n
288 # end
289 #
290 # @testset "BoundaryQuadrature" begin
291 # op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
292 # g = EquidistantGrid((10,11), (0.0, 0.0), (1.0,1.0))
293 #
294 # H_w = BoundaryQuadrature(op, g, CartesianBoundary{1,Lower}())
295 # H_e = BoundaryQuadrature(op, g, CartesianBoundary{1,Upper}())
296 # H_s = BoundaryQuadrature(op, g, CartesianBoundary{2,Lower}())
297 # H_n = BoundaryQuadrature(op, g, CartesianBoundary{2,Upper}())
298 #
299 # v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
300 #
301 # function get_quadrature(N)
302 # qc = op.quadratureClosure
303 # q = (qc..., ones(N-2*closuresize(op))..., reverse(qc)...)
304 # @assert length(q) == N
305 # return q
306 # end
307 #
308 # v_w = v[1,:]
309 # v_e = v[10,:]
310 # v_s = v[:,1]
311 # v_n = v[:,11]
312 #
313 # q_x = spacing(g)[1].*get_quadrature(10)
314 # q_y = spacing(g)[2].*get_quadrature(11)
315 #
316 # @test H_w isa TensorOperator{T,1} where T
317 #
318 # @test domain_size(H_w, (3,)) == (3,)
319 # @test domain_size(H_n, (3,)) == (3,)
320 #
321 # @test range_size(H_w, (3,)) == (3,)
322 # @test range_size(H_n, (3,)) == (3,)
323 #
324 # @test size(H_w*v_w) == (11,)
325 # @test size(H_e*v_e) == (11,)
326 # @test size(H_s*v_s) == (10,)
327 # @test size(H_n*v_n) == (10,)
328 #
329 # @test collect(H_w*v_w) ≈ q_y.*v_w
330 # @test collect(H_e*v_e) ≈ q_y.*v_e
331 # @test collect(H_s*v_s) ≈ q_x.*v_s
332 # @test collect(H_n*v_n) ≈ q_x.*v_n
333 #
334 # @test collect(H_w'*v_w) == collect(H_w'*v_w)
335 # @test collect(H_e'*v_e) == collect(H_e'*v_e)
336 # @test collect(H_s'*v_s) == collect(H_s'*v_s)
337 # @test collect(H_n'*v_n) == collect(H_n'*v_n)
338 # end