Mercurial > repos > public > sbplib_julia

comparison test/testDiffOps.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 DiffOps/test/runtests.jl@e21dcda55163
children f4e3e71a4ff4
comparison
equal deleted inserted replaced
332:535f1bff4bcc 333:01b851161018
1 using Test
2 using Sbplib
3 using DiffOps
4 using Grids
5 using SbpOperators
6 using RegionIndices
7 using LazyTensors
8
9 @testset "Laplace2D" begin
10 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
11 Lx = 3.5
12 Ly = 7.2
13 g = EquidistantGrid((42,41), (0.0, 0.0), (Lx,Ly))
14 L = Laplace(g, 1., op)
15 H = quadrature(L)
16
17 f0(x::Float64,y::Float64) = 2.
18 f1(x::Float64,y::Float64) = x+y
19 f2(x::Float64,y::Float64) = 1/2*x^2 + 1/2*y^2
20 f3(x::Float64,y::Float64) = 1/6*x^3 + 1/6*y^3
21 f4(x::Float64,y::Float64) = 1/24*x^4 + 1/24*y^4
22 f5(x::Float64,y::Float64) = sin(x) + cos(y)
23 f5ₓₓ(x::Float64,y::Float64) = -f5(x,y)
24
25 v0 = evalOn(g,f0)
26 v1 = evalOn(g,f1)
27 v2 = evalOn(g,f2)
28 v3 = evalOn(g,f3)
29 v4 = evalOn(g,f4)
30 v5 = evalOn(g,f5)
31 v5ₓₓ = evalOn(g,f5ₓₓ)
32
33 @test L isa TensorOperator{T,2} where T
34 @test L' isa TensorMapping{T,2,2} where T
35
36 # TODO: Should perhaps set tolerance level for isapporx instead?
37 # Are these tolerance levels resonable or should tests be constructed
38 # differently?
39 equalitytol = 0.5*1e-10
40 accuracytol = 0.5*1e-3
41 # 4th order interior stencil, 2nd order boundary stencil,
42 # implies that L*v should be exact for v - monomial up to order 3.
43 # Exact differentiation is measured point-wise. For other grid functions
44 # the error is measured in the H-norm.
45 @test all(abs.(collect(L*v0)) .<= equalitytol)
46 @test all(abs.(collect(L*v1)) .<= equalitytol)
47 @test all(collect(L*v2) .≈ v0) # Seems to be more accurate
48 @test all(abs.((collect(L*v3) - v1)) .<= equalitytol)
49 e4 = collect(L*v4) - v2
50 e5 = collect(L*v5) - v5ₓₓ
51 @test sum(collect(H*e4.^2)) <= accuracytol
52 @test sum(collect(H*e5.^2)) <= accuracytol
53 end
54
55 @testset "Quadrature" begin
56 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
57 Lx = 2.3
58 Ly = 5.2
59 g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
60 H = Quadrature(op,g)
61 v = ones(Float64, size(g))
62
63 @test H isa TensorOperator{T,2} where T
64 @test H' isa TensorMapping{T,2,2} where T
65 @test sum(collect(H*v)) ≈ (Lx*Ly)
66 @test collect(H*v) == collect(H'*v)
67 end
68
69 @testset "InverseQuadrature" begin
70 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
71 Lx = 7.3
72 Ly = 8.2
73 g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
74 H = Quadrature(op,g)
75 Hinv = InverseQuadrature(op,g)
76 v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
77
78 @test Hinv isa TensorOperator{T,2} where T
79 @test Hinv' isa TensorMapping{T,2,2} where T
80 @test collect(Hinv*H*v) ≈ v
81 @test collect(Hinv*v) == collect(Hinv'*v)
82 end
83
84 @testset "BoundaryValue" begin
85 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
86 g = EquidistantGrid((4,5), (0.0, 0.0), (1.0,1.0))
87
88 e_w = BoundaryValue(op, g, CartesianBoundary{1,Lower}())
89 e_e = BoundaryValue(op, g, CartesianBoundary{1,Upper}())
90 e_s = BoundaryValue(op, g, CartesianBoundary{2,Lower}())
91 e_n = BoundaryValue(op, g, CartesianBoundary{2,Upper}())
92
93 v = zeros(Float64, 4, 5)
94 v[:,5] = [1, 2, 3,4]
95 v[:,4] = [1, 2, 3,4]
96 v[:,3] = [4, 5, 6, 7]
97 v[:,2] = [7, 8, 9, 10]
98 v[:,1] = [10, 11, 12, 13]
99
100 @test e_w isa TensorMapping{T,2,1} where T
101 @test e_w' isa TensorMapping{T,1,2} where T
102
103 @test domain_size(e_w, (3,2)) == (2,)
104 @test domain_size(e_e, (3,2)) == (2,)
105 @test domain_size(e_s, (3,2)) == (3,)
106 @test domain_size(e_n, (3,2)) == (3,)
107
108 @test size(e_w'*v) == (5,)
109 @test size(e_e'*v) == (5,)
110 @test size(e_s'*v) == (4,)
111 @test size(e_n'*v) == (4,)
112
113 @test collect(e_w'*v) == [10,7,4,1.0,1]
114 @test collect(e_e'*v) == [13,10,7,4,4.0]
115 @test collect(e_s'*v) == [10,11,12,13.0]
116 @test collect(e_n'*v) == [1,2,3,4.0]
117
118 g_x = [1,2,3,4.0]
119 g_y = [5,4,3,2,1.0]
120
121 G_w = zeros(Float64, (4,5))
122 G_w[1,:] = g_y
123
124 G_e = zeros(Float64, (4,5))
125 G_e[4,:] = g_y
126
127 G_s = zeros(Float64, (4,5))
128 G_s[:,1] = g_x
129
130 G_n = zeros(Float64, (4,5))
131 G_n[:,5] = g_x
132
133 @test size(e_w*g_y) == (UnknownDim,5)
134 @test size(e_e*g_y) == (UnknownDim,5)
135 @test size(e_s*g_x) == (4,UnknownDim)
136 @test size(e_n*g_x) == (4,UnknownDim)
137
138 # These tests should be moved to where they are possible (i.e we know what the grid should be)
139 @test_broken collect(e_w*g_y) == G_w
140 @test_broken collect(e_e*g_y) == G_e
141 @test_broken collect(e_s*g_x) == G_s
142 @test_broken collect(e_n*g_x) == G_n
143 end
144
145 @testset "NormalDerivative" begin
146 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
147 g = EquidistantGrid((5,6), (0.0, 0.0), (4.0,5.0))
148
149 d_w = NormalDerivative(op, g, CartesianBoundary{1,Lower}())
150 d_e = NormalDerivative(op, g, CartesianBoundary{1,Upper}())
151 d_s = NormalDerivative(op, g, CartesianBoundary{2,Lower}())
152 d_n = NormalDerivative(op, g, CartesianBoundary{2,Upper}())
153
154
155 v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
156 v∂x = evalOn(g, (x,y)-> 2*x + y)
157 v∂y = evalOn(g, (x,y)-> 2*(y-1) + x)
158
159 @test d_w isa TensorMapping{T,2,1} where T
160 @test d_w' isa TensorMapping{T,1,2} where T
161
162 @test domain_size(d_w, (3,2)) == (2,)
163 @test domain_size(d_e, (3,2)) == (2,)
164 @test domain_size(d_s, (3,2)) == (3,)
165 @test domain_size(d_n, (3,2)) == (3,)
166
167 @test size(d_w'*v) == (6,)
168 @test size(d_e'*v) == (6,)
169 @test size(d_s'*v) == (5,)
170 @test size(d_n'*v) == (5,)
171
172 @test collect(d_w'*v) ≈ v∂x[1,:]
173 @test collect(d_e'*v) ≈ v∂x[5,:]
174 @test collect(d_s'*v) ≈ v∂y[:,1]
175 @test collect(d_n'*v) ≈ v∂y[:,6]
176
177
178 d_x_l = zeros(Float64, 5)
179 d_x_u = zeros(Float64, 5)
180 for i ∈ eachindex(d_x_l)
181 d_x_l[i] = op.dClosure[i-1]
182 d_x_u[i] = -op.dClosure[length(d_x_u)-i]
183 end
184
185 d_y_l = zeros(Float64, 6)
186 d_y_u = zeros(Float64, 6)
187 for i ∈ eachindex(d_y_l)
188 d_y_l[i] = op.dClosure[i-1]
189 d_y_u[i] = -op.dClosure[length(d_y_u)-i]
190 end
191
192 function prod_matrix(x,y)
193 G = zeros(Float64, length(x), length(y))
194 for I ∈ CartesianIndices(G)
195 G[I] = x[I[1]]*y[I[2]]
196 end
197
198 return G
199 end
200
201 g_x = [1,2,3,4.0,5]
202 g_y = [5,4,3,2,1.0,11]
203
204 G_w = prod_matrix(d_x_l, g_y)
205 G_e = prod_matrix(d_x_u, g_y)
206 G_s = prod_matrix(g_x, d_y_l)
207 G_n = prod_matrix(g_x, d_y_u)
208
209
210 @test size(d_w*g_y) == (UnknownDim,6)
211 @test size(d_e*g_y) == (UnknownDim,6)
212 @test size(d_s*g_x) == (5,UnknownDim)
213 @test size(d_n*g_x) == (5,UnknownDim)
214
215 # These tests should be moved to where they are possible (i.e we know what the grid should be)
216 @test_broken collect(d_w*g_y) ≈ G_w
217 @test_broken collect(d_e*g_y) ≈ G_e
218 @test_broken collect(d_s*g_x) ≈ G_s
219 @test_broken collect(d_n*g_x) ≈ G_n
220 end
221
222 @testset "BoundaryQuadrature" begin
223 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
224 g = EquidistantGrid((10,11), (0.0, 0.0), (1.0,1.0))
225
226 H_w = BoundaryQuadrature(op, g, CartesianBoundary{1,Lower}())
227 H_e = BoundaryQuadrature(op, g, CartesianBoundary{1,Upper}())
228 H_s = BoundaryQuadrature(op, g, CartesianBoundary{2,Lower}())
229 H_n = BoundaryQuadrature(op, g, CartesianBoundary{2,Upper}())
230
231 v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
232
233 function get_quadrature(N)
234 qc = op.quadratureClosure
235 q = (qc..., ones(N-2*closuresize(op))..., reverse(qc)...)
236 @assert length(q) == N
237 return q
238 end
239
240 v_w = v[1,:]
241 v_e = v[10,:]
242 v_s = v[:,1]
243 v_n = v[:,11]
244
245 q_x = spacing(g)[1].*get_quadrature(10)
246 q_y = spacing(g)[2].*get_quadrature(11)
247
248 @test H_w isa TensorOperator{T,1} where T
249
250 @test domain_size(H_w, (3,)) == (3,)
251 @test domain_size(H_n, (3,)) == (3,)
252
253 @test range_size(H_w, (3,)) == (3,)
254 @test range_size(H_n, (3,)) == (3,)
255
256 @test size(H_w*v_w) == (11,)
257 @test size(H_e*v_e) == (11,)
258 @test size(H_s*v_s) == (10,)
259 @test size(H_n*v_n) == (10,)
260
261 @test collect(H_w*v_w) ≈ q_y.*v_w
262 @test collect(H_e*v_e) ≈ q_y.*v_e
263 @test collect(H_s*v_s) ≈ q_x.*v_s
264 @test collect(H_n*v_n) ≈ q_x.*v_n
265
266 @test collect(H_w'*v_w) == collect(H_w'*v_w)
267 @test collect(H_e'*v_e) == collect(H_e'*v_e)
268 @test collect(H_s'*v_s) == collect(H_s'*v_s)
269 @test collect(H_n'*v_n) == collect(H_n'*v_n)
270 end