Mercurial > repos > public > sbplib_julia

comparison test/testSbpOperators.jl @ 628:316dbfd31d35 feature/volume_and_boundary_operators

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add tests for VolumeOperator
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 31 Dec 2020 08:22:56 +0100
parents f799678357df
children 22a0971f7f84
comparison
equal deleted inserted replaced
627:9f27f451d0a0 628:316dbfd31d35
6 using LinearAlgebra 6 using LinearAlgebra
7 using TOML 7 using TOML
8 8
9 import Sbplib.SbpOperators.Stencil 9 import Sbplib.SbpOperators.Stencil
10 import Sbplib.SbpOperators.VolumeOperator 10 import Sbplib.SbpOperators.VolumeOperator
11 import Sbplib.SbpOperators.volume_operator
11 import Sbplib.SbpOperators.BoundaryOperator 12 import Sbplib.SbpOperators.BoundaryOperator
12 import Sbplib.SbpOperators.boundary_operator 13 import Sbplib.SbpOperators.boundary_operator
14 import Sbplib.SbpOperators.Parity
15
13 16
14 @testset "SbpOperators" begin 17 @testset "SbpOperators" begin
15 18
16 @testset "Stencil" begin 19 @testset "Stencil" begin
17 s = Stencil((-2,2), (1.,2.,2.,3.,4.)) 20 s = Stencil((-2,2), (1.,2.,2.,3.,4.))
129 # v = [2.,3.,2.,4.,5.,4.,3.,4.,5.,4.5] 132 # v = [2.,3.,2.,4.,5.,4.,3.,4.,5.,4.5]
130 # for i ∈ 1:N 133 # for i ∈ 1:N
131 # @test apply_quadrature(op, h, v[i], i, N) == q[i]*v[i] 134 # @test apply_quadrature(op, h, v[i], i, N) == q[i]*v[i]
132 # end 135 # end
133 # end 136 # end
137
138 @testset "VolumeOperator" begin
139 inner_stencil = Stencil(1/4 .* (1.,2.,1.),center=2)
140 closure_stencils = (Stencil(1/2 .* (1.,1.),center=1),Stencil((0.,1.),center=2))
141 g_1D = EquidistantGrid(11,0.,1.)
142 g_2D = EquidistantGrid((11,12),(0.,0.),(1.,1.))
143 g_3D = EquidistantGrid((11,12,10),(0.,0.,0.),(1.,1.,1.))
144 @testset "Constructors" begin
145 #TODO: How are even and odd in SbpOperators.Parity exposed? Currently constructing even as Parity(1)
146 @testset "1D" begin
147 op = VolumeOperator(inner_stencil,closure_stencils,(11,),Parity(1))
148 @test op == VolumeOperator(g_1D,inner_stencil,closure_stencils,Parity(1))
149 @test op == volume_operator(g_1D,inner_stencil,closure_stencils,Parity(1),1)
150 @test op isa TensorMapping{T,1,1} where T
151 end
152 @testset "2D" begin
153 op_x = volume_operator(g_2D,inner_stencil,closure_stencils,Parity(1),1)
154 op_y = volume_operator(g_2D,inner_stencil,closure_stencils,Parity(1),2)
155 Ix = IdentityMapping{Float64}((11,))
156 Iy = IdentityMapping{Float64}((12,))
157 @test op_x == VolumeOperator(inner_stencil,closure_stencils,(11,),Parity(1))⊗Iy
158 @test op_y == Ix⊗VolumeOperator(inner_stencil,closure_stencils,(12,),Parity(1))
159 @test op_x isa TensorMapping{T,2,2} where T
160 @test op_y isa TensorMapping{T,2,2} where T
161 end
162 @testset "3D" begin
163 op_x = volume_operator(g_3D,inner_stencil,closure_stencils,Parity(1),1)
164 op_y = volume_operator(g_3D,inner_stencil,closure_stencils,Parity(1),2)
165 op_z = volume_operator(g_3D,inner_stencil,closure_stencils,Parity(1),3)
166 Ix = IdentityMapping{Float64}((11,))
167 Iy = IdentityMapping{Float64}((12,))
168 Iz = IdentityMapping{Float64}((10,))
169 @test op_x == VolumeOperator(inner_stencil,closure_stencils,(11,),Parity(1))⊗Iy⊗Iz
170 @test op_y == Ix⊗VolumeOperator(inner_stencil,closure_stencils,(12,),Parity(1))⊗Iz
171 @test op_z == Ix⊗Iy⊗VolumeOperator(inner_stencil,closure_stencils,(10,),Parity(1))
172 @test op_x isa TensorMapping{T,3,3} where T
173 @test op_y isa TensorMapping{T,3,3} where T
174 @test op_z isa TensorMapping{T,3,3} where T
175 end
176 end
177
178 @testset "Sizes" begin
179 @testset "1D" begin
180 op = volume_operator(g_1D,inner_stencil,closure_stencils,Parity(1),1)
181 @test range_size(op) == domain_size(op) == size(g_1D)
182 end
183
184 @testset "2D" begin
185 op_x = volume_operator(g_2D,inner_stencil,closure_stencils,Parity(1),1)
186 op_y = volume_operator(g_2D,inner_stencil,closure_stencils,Parity(1),2)
187 @test range_size(op_y) == domain_size(op_y) ==
188 range_size(op_x) == domain_size(op_x) == size(g_2D)
189 end
190 @testset "3D" begin
191 op_x = volume_operator(g_3D,inner_stencil,closure_stencils,Parity(1),1)
192 op_y = volume_operator(g_3D,inner_stencil,closure_stencils,Parity(1),2)
193 op_z = volume_operator(g_3D,inner_stencil,closure_stencils,Parity(1),3)
194 @test range_size(op_z) == domain_size(op_z) ==
195 range_size(op_y) == domain_size(op_y) ==
196 range_size(op_x) == domain_size(op_x) == size(g_3D)
197 end
198 end
199
200 # TODO: Test for other dimensions?
201 op_x = volume_operator(g_2D,inner_stencil,closure_stencils,Parity(1),1)
202 op_y = volume_operator(g_2D,inner_stencil,closure_stencils,Parity(1),2)
203 v = zeros(size(g_2D))
204 Nx = size(g_2D)[1]
205 for i = 1:Nx
206 v[i,:] .= i
207 end
208 rx = copy(v)
209 rx[1,:] .= 1.5
210 rx[end,:] .= (2*Nx-1)/2
211 ry = copy(v)
212
213 @testset "Application" begin
214 @test op_x*v ≈ rx rtol = 1e-14
215 @test op_y*v ≈ ry rtol = 1e-14
216 end
217
218 # TODO: Test for other dimensions?
219 @testset "Regions" begin
220 @test (op_x*v)[Index(1,Lower),Index(3,Interior)] ≈ rx[1,3] rtol = 1e-14
221 @test (op_x*v)[Index(2,Lower),Index(3,Interior)] ≈ rx[2,3] rtol = 1e-14
222 @test (op_x*v)[Index(6,Interior),Index(3,Interior)] ≈ rx[6,3] rtol = 1e-14
223 @test (op_x*v)[Index(10,Upper),Index(3,Interior)] ≈ rx[10,3] rtol = 1e-14
224 @test (op_x*v)[Index(11,Upper),Index(3,Interior)] ≈ rx[11,3] rtol = 1e-14
225
226 @test_throws BoundsError (op_x*v)[Index(3,Lower),Index(3,Interior)]
227 @test_throws BoundsError (op_x*v)[Index(9,Upper),Index(3,Interior)]
228
229 @test (op_y*v)[Index(3,Interior),Index(1,Lower)] ≈ ry[3,1] rtol = 1e-14
230 @test (op_y*v)[Index(3,Interior),Index(2,Lower)] ≈ ry[3,2] rtol = 1e-14
231 @test (op_y*v)[Index(3,Interior),Index(6,Interior)] ≈ ry[3,6] rtol = 1e-14
232 @test (op_y*v)[Index(3,Interior),Index(11,Upper)] ≈ ry[3,11] rtol = 1e-14
233 @test (op_y*v)[Index(3,Interior),Index(12,Upper)] ≈ ry[3,12] rtol = 1e-14
234
235 @test_throws BoundsError (op_y*v)[Index(3,Interior),Index(10,Upper)]
236 @test_throws BoundsError (op_y*v)[Index(3,Interior),Index(3,Lower)]
237 end
238
239 # TODO: Test for other dimensions?
240 @testset "Inferred" begin
241 @inferred apply(op_x, v,1,1)
242 @inferred apply(op_x, v, Index(1,Lower),Index(1,Lower))
243 @inferred apply(op_x, v, Index(6,Interior),Index(1,Lower))
244 @inferred apply(op_x, v, Index(11,Upper),Index(1,Lower))
245
246 @inferred apply(op_y, v,1,1)
247 @inferred apply(op_y, v, Index(1,Lower),Index(1,Lower))
248 @inferred apply(op_y, v, Index(1,Lower),Index(6,Interior))
249 @inferred apply(op_y, v, Index(1,Lower),Index(11,Upper))
250 end
251
252 end
134 253
135 @testset "SecondDerivative" begin 254 @testset "SecondDerivative" begin
136 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 255 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
137 Lx = 3.5 256 Lx = 3.5
138 Ly = 3. 257 Ly = 3.