Mercurial > repos > public > sbplib_julia

comparison test/testSbpOperators.jl @ 702:3cd582257072 feature/laplace_opset

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge in default
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 15 Feb 2021 11:30:34 +0100
parents 54ce3f9771e5 5ddf28ddee18
children 988e9cfcd58d
comparison
equal deleted inserted replaced
679:54ce3f9771e5 702:3cd582257072
22 @test s isa Stencil{Float64, 5} 22 @test s isa Stencil{Float64, 5}
23 23
24 @test eltype(s) == Float64 24 @test eltype(s) == Float64
25 @test SbpOperators.scale(s, 2) == Stencil((-2,2), (2.,4.,4.,6.,8.)) 25 @test SbpOperators.scale(s, 2) == Stencil((-2,2), (2.,4.,4.,6.,8.))
26 26
27 @test Stencil((1,2,3,4), center=1) == Stencil((0, 3),(1,2,3,4)) 27 @test Stencil(1,2,3,4; center=1) == Stencil((0, 3),(1,2,3,4))
28 @test Stencil((1,2,3,4), center=2) == Stencil((-1, 2),(1,2,3,4)) 28 @test Stencil(1,2,3,4; center=2) == Stencil((-1, 2),(1,2,3,4))
29 @test Stencil((1,2,3,4), center=4) == Stencil((-3, 0),(1,2,3,4)) 29 @test Stencil(1,2,3,4; center=4) == Stencil((-3, 0),(1,2,3,4))
30
31 @test CenteredStencil(1,2,3,4,5) == Stencil((-2, 2), (1,2,3,4,5))
32 @test_throws ArgumentError CenteredStencil(1,2,3,4)
30 end 33 end
31 34
32 @testset "parse_rational" begin 35 @testset "parse_rational" begin
33 @test SbpOperators.parse_rational("1") isa Rational 36 @test SbpOperators.parse_rational("1") isa Rational
34 @test SbpOperators.parse_rational("1") == 1//1 37 @test SbpOperators.parse_rational("1") == 1//1
65 ] 68 ]
66 """ 69 """
67 70
68 parsed_toml = TOML.parse(toml_str) 71 parsed_toml = TOML.parse(toml_str)
69 @testset "get_stencil" begin 72 @testset "get_stencil" begin
70 @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil") == Stencil((-1/2, 0., 1/2), center=2) 73 @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil") == Stencil(-1/2, 0., 1/2, center=2)
71 @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil", center=1) == Stencil((-1/2, 0., 1/2); center=1) 74 @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil", center=1) == Stencil(-1/2, 0., 1/2; center=1)
72 @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil", center=3) == Stencil((-1/2, 0., 1/2); center=3) 75 @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil", center=3) == Stencil(-1/2, 0., 1/2; center=3)
73 76
74 @test get_stencil(parsed_toml, "order2", "H", "inner") == Stencil((1.,), center=1) 77 @test get_stencil(parsed_toml, "order2", "H", "inner") == Stencil(1.; center=1)
75 78
76 @test_throws AssertionError get_stencil(parsed_toml, "meta", "type") 79 @test_throws AssertionError get_stencil(parsed_toml, "meta", "type")
77 @test_throws AssertionError get_stencil(parsed_toml, "order2", "D1", "closure_stencils") 80 @test_throws AssertionError get_stencil(parsed_toml, "order2", "D1", "closure_stencils")
78 end 81 end
79 82
80 @testset "get_stencils" begin 83 @testset "get_stencils" begin
81 @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=(1,)) == (Stencil((-1., 1.), center=1),) 84 @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=(1,)) == (Stencil(-1., 1., center=1),)
82 @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=(2,)) == (Stencil((-1., 1.), center=2),) 85 @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=(2,)) == (Stencil(-1., 1., center=2),)
83 @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=[2]) == (Stencil((-1., 1.), center=2),) 86 @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=[2]) == (Stencil(-1., 1., center=2),)
84 87
85 @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=[1,1,1,1]) == ( 88 @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=[1,1,1,1]) == (
86 Stencil(( 2., -5., 4., -1., 0., 0.), center=1), 89 Stencil( 2., -5., 4., -1., 0., 0., center=1),
87 Stencil(( 1., -2., 1., 0., 0., 0.), center=1), 90 Stencil( 1., -2., 1., 0., 0., 0., center=1),
88 Stencil(( -4/43, 59/43, -110/43, 59/43, -4/43, 0.), center=1), 91 Stencil( -4/43, 59/43, -110/43, 59/43, -4/43, 0., center=1),
89 Stencil(( -1/49, 0., 59/49, -118/49, 64/49, -4/49), center=1), 92 Stencil( -1/49, 0., 59/49, -118/49, 64/49, -4/49, center=1),
90 ) 93 )
91 94
92 @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(4,2,3,1)) == ( 95 @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(4,2,3,1)) == (
93 Stencil(( 2., -5., 4., -1., 0., 0.), center=4), 96 Stencil( 2., -5., 4., -1., 0., 0., center=4),
94 Stencil(( 1., -2., 1., 0., 0., 0.), center=2), 97 Stencil( 1., -2., 1., 0., 0., 0., center=2),
95 Stencil(( -4/43, 59/43, -110/43, 59/43, -4/43, 0.), center=3), 98 Stencil( -4/43, 59/43, -110/43, 59/43, -4/43, 0., center=3),
96 Stencil(( -1/49, 0., 59/49, -118/49, 64/49, -4/49), center=1), 99 Stencil( -1/49, 0., 59/49, -118/49, 64/49, -4/49, center=1),
97 ) 100 )
98 101
99 @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=1:4) == ( 102 @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=1:4) == (
100 Stencil(( 2., -5., 4., -1., 0., 0.), center=1), 103 Stencil( 2., -5., 4., -1., 0., 0., center=1),
101 Stencil(( 1., -2., 1., 0., 0., 0.), center=2), 104 Stencil( 1., -2., 1., 0., 0., 0., center=2),
102 Stencil(( -4/43, 59/43, -110/43, 59/43, -4/43, 0.), center=3), 105 Stencil( -4/43, 59/43, -110/43, 59/43, -4/43, 0., center=3),
103 Stencil(( -1/49, 0., 59/49, -118/49, 64/49, -4/49), center=4), 106 Stencil( -1/49, 0., 59/49, -118/49, 64/49, -4/49, center=4),
104 ) 107 )
105 108
106 @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(1,2,3)) 109 @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(1,2,3))
107 @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(1,2,3,5,4)) 110 @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(1,2,3,5,4))
108 @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "inner_stencil",centers=(1,2)) 111 @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "inner_stencil",centers=(1,2))
114 @test_throws AssertionError get_tuple(parsed_toml, "meta", "type") 117 @test_throws AssertionError get_tuple(parsed_toml, "meta", "type")
115 end 118 end
116 end 119 end
117 120
118 @testset "VolumeOperator" begin 121 @testset "VolumeOperator" begin
119 inner_stencil = Stencil(1/4 .* (1.,2.,1.),center=2) 122 inner_stencil = CenteredStencil(1/4, 2/4, 1/4)
120 closure_stencils = (Stencil(1/2 .* (1.,1.),center=1),Stencil((0.,1.),center=2)) 123 closure_stencils = (Stencil(1/2, 1/2; center=1), Stencil(0.,1.; center=2))
121 g_1D = EquidistantGrid(11,0.,1.) 124 g_1D = EquidistantGrid(11,0.,1.)
122 g_2D = EquidistantGrid((11,12),(0.,0.),(1.,1.)) 125 g_2D = EquidistantGrid((11,12),(0.,0.),(1.,1.))
123 g_3D = EquidistantGrid((11,12,10),(0.,0.,0.),(1.,1.,1.)) 126 g_3D = EquidistantGrid((11,12,10),(0.,0.,0.),(1.,1.,1.))
124 @testset "Constructors" begin 127 @testset "Constructors" begin
125 @testset "1D" begin 128 @testset "1D" begin
236 g_1D = EquidistantGrid(121, 0.0, Lx) 239 g_1D = EquidistantGrid(121, 0.0, Lx)
237 g_2D = EquidistantGrid((121,123), (0.0, 0.0), (Lx, Ly)) 240 g_2D = EquidistantGrid((121,123), (0.0, 0.0), (Lx, Ly))
238 241
239 @testset "Constructors" begin 242 @testset "Constructors" begin
240 @testset "1D" begin 243 @testset "1D" begin
241 Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils) 244 Dₓₓ = second_derivative(g_1D,op.innerStencil,op.closureStencils)
242 @test Dₓₓ == SecondDerivative(g_1D,op.innerStencil,op.closureStencils,1) 245 @test Dₓₓ == second_derivative(g_1D,op.innerStencil,op.closureStencils,1)
243 @test Dₓₓ isa VolumeOperator 246 @test Dₓₓ isa VolumeOperator
244 end 247 end
245 @testset "2D" begin 248 @testset "2D" begin
246 Dₓₓ = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,1) 249 Dₓₓ = second_derivative(g_2D,op.innerStencil,op.closureStencils,1)
247 D2 = SecondDerivative(g_1D,op.innerStencil,op.closureStencils) 250 D2 = second_derivative(g_1D,op.innerStencil,op.closureStencils)
248 I = IdentityMapping{Float64}(size(g_2D)[2]) 251 I = IdentityMapping{Float64}(size(g_2D)[2])
249 @test Dₓₓ == D2⊗I 252 @test Dₓₓ == D2⊗I
250 @test Dₓₓ isa TensorMapping{T,2,2} where T 253 @test Dₓₓ isa TensorMapping{T,2,2} where T
251 end 254 end
252 end 255 end
267 270
268 # 2nd order interior stencil, 1nd order boundary stencil, 271 # 2nd order interior stencil, 1nd order boundary stencil,
269 # implies that L*v should be exact for monomials up to order 2. 272 # implies that L*v should be exact for monomials up to order 2.
270 @testset "2nd order" begin 273 @testset "2nd order" begin
271 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) 274 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
272 Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils) 275 Dₓₓ = second_derivative(g_1D,op.innerStencil,op.closureStencils)
273 @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10 276 @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
274 @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10 277 @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
275 @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10 278 @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10
276 @test Dₓₓ*v ≈ vₓₓ rtol = 5e-2 norm = l2 279 @test Dₓₓ*v ≈ vₓₓ rtol = 5e-2 norm = l2
277 end 280 end
278 281
279 # 4th order interior stencil, 2nd order boundary stencil, 282 # 4th order interior stencil, 2nd order boundary stencil,
280 # implies that L*v should be exact for monomials up to order 3. 283 # implies that L*v should be exact for monomials up to order 3.
281 @testset "4th order" begin 284 @testset "4th order" begin
282 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 285 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
283 Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils) 286 Dₓₓ = second_derivative(g_1D,op.innerStencil,op.closureStencils)
284 # NOTE: high tolerances for checking the "exact" differentiation 287 # NOTE: high tolerances for checking the "exact" differentiation
285 # due to accumulation of round-off errors/cancellation errors? 288 # due to accumulation of round-off errors/cancellation errors?
286 @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10 289 @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
287 @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10 290 @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
288 @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10 291 @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10
304 307
305 # 2nd order interior stencil, 1st order boundary stencil, 308 # 2nd order interior stencil, 1st order boundary stencil,
306 # implies that L*v should be exact for binomials up to order 2. 309 # implies that L*v should be exact for binomials up to order 2.
307 @testset "2nd order" begin 310 @testset "2nd order" begin
308 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) 311 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
309 Dyy = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,2) 312 Dyy = second_derivative(g_2D,op.innerStencil,op.closureStencils,2)
310 @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9 313 @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
311 @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9 314 @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
312 @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9 315 @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9
313 @test Dyy*v ≈ v_yy rtol = 5e-2 norm = l2 316 @test Dyy*v ≈ v_yy rtol = 5e-2 norm = l2
314 end 317 end
315 318
316 # 4th order interior stencil, 2nd order boundary stencil, 319 # 4th order interior stencil, 2nd order boundary stencil,
317 # implies that L*v should be exact for binomials up to order 3. 320 # implies that L*v should be exact for binomials up to order 3.
318 @testset "4th order" begin 321 @testset "4th order" begin
319 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 322 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
320 Dyy = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,2) 323 Dyy = second_derivative(g_2D,op.innerStencil,op.closureStencils,2)
321 # NOTE: high tolerances for checking the "exact" differentiation 324 # NOTE: high tolerances for checking the "exact" differentiation
322 # due to accumulation of round-off errors/cancellation errors? 325 # due to accumulation of round-off errors/cancellation errors?
323 @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9 326 @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
324 @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9 327 @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
325 @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9 328 @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9
336 @testset "Constructors" begin 339 @testset "Constructors" begin
337 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 340 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
338 @testset "1D" begin 341 @testset "1D" begin
339 # Create all tensor mappings included in Laplace 342 # Create all tensor mappings included in Laplace
340 Δ = laplace(g_1D, op.innerStencil, op.closureStencils) 343 Δ = laplace(g_1D, op.innerStencil, op.closureStencils)
341 H = quadrature(g_1D, op.quadratureClosure) 344 H = inner_product(g_1D, op.quadratureClosure)
342 Hi = InverseDiagonalQuadrature(g_1D, op.quadratureClosure) 345 Hi = inverse_inner_product(g_1D, op.quadratureClosure)
343 346
344 (id_l, id_r) = boundary_identifiers(g_1D) 347 (id_l, id_r) = boundary_identifiers(g_1D)
345 348
346 e_l = BoundaryRestriction(g_1D,op.eClosure,id_l) 349 e_l = boundary_restriction(g_1D,op.eClosure,id_l)
347 e_r = BoundaryRestriction(g_1D,op.eClosure,id_r) 350 e_r = boundary_restriction(g_1D,op.eClosure,id_r)
348 e_dict = Dict(Pair(id_l,e_l),Pair(id_r,e_r)) 351 e_dict = Dict(Pair(id_l,e_l),Pair(id_r,e_r))
349 352
350 d_l = NormalDerivative(g_1D,op.dClosure,id_l) 353 d_l = normal_derivative(g_1D,op.dClosure,id_l)
351 d_r = NormalDerivative(g_1D,op.dClosure,id_r) 354 d_r = normal_derivative(g_1D,op.dClosure,id_r)
352 d_dict = Dict(Pair(id_l,d_l),Pair(id_r,d_r)) 355 d_dict = Dict(Pair(id_l,d_l),Pair(id_r,d_r))
353 356
354 H_l = boundary_quadrature(g_1D,op.quadratureClosure,id_r) 357 H_l = inner_product(boundary_grid(g_1D,id_l),op.quadratureClosure)
355 H_r = boundary_quadrature(g_1D,op.quadratureClosure,id_r) 358 H_r = inner_product(boundary_grid(g_1D,id_r),op.quadratureClosure)
356 Hb_dict = Dict(Pair(id_l,H_l),Pair(id_r,H_r)) 359 Hb_dict = Dict(Pair(id_l,H_l),Pair(id_r,H_r))
357 360
358 # TODO: Not sure why this doesnt work? Comparing the fields of 361 # TODO: Not sure why this doesnt work? Comparing the fields of
359 # Laplace seems to work. Reformulate below once solved. 362 # Laplace seems to work. Reformulate below once solved.
360 @test_broken Laplace(Δ,H,Hi,e_dict,d_dict,Hb_dict) == Laplace(g_1D, sbp_operators_path()*"standard_diagonal.toml"; order=4) 363 @test_broken Laplace(Δ,H,Hi,e_dict,d_dict,Hb_dict) == Laplace(g_1D, sbp_operators_path()*"standard_diagonal.toml"; order=4)
370 @inferred Laplace(Δ,H,Hi,e_dict,d_dict,Hb_dict) 373 @inferred Laplace(Δ,H,Hi,e_dict,d_dict,Hb_dict)
371 end 374 end
372 @testset "3D" begin 375 @testset "3D" begin
373 # Create all tensor mappings included in Laplace 376 # Create all tensor mappings included in Laplace
374 Δ = laplace(g_3D, op.innerStencil, op.closureStencils) 377 Δ = laplace(g_3D, op.innerStencil, op.closureStencils)
375 H = quadrature(g_3D, op.quadratureClosure) 378 H = inner_product(g_3D, op.quadratureClosure)
376 Hi = InverseDiagonalQuadrature(g_3D, op.quadratureClosure) 379 Hi = inverse_inner_product(g_3D, op.quadratureClosure)
377 380
378 (id_l, id_r, id_s, id_n, id_b, id_t) = boundary_identifiers(g_3D) 381 (id_l, id_r, id_s, id_n, id_b, id_t) = boundary_identifiers(g_3D)
379 382
380 e_l = BoundaryRestriction(g_3D,op.eClosure,id_l) 383 e_l = boundary_restriction(g_3D,op.eClosure,id_l)
381 e_r = BoundaryRestriction(g_3D,op.eClosure,id_r) 384 e_r = boundary_restriction(g_3D,op.eClosure,id_r)
382 e_s = BoundaryRestriction(g_3D,op.eClosure,id_s) 385 e_s = boundary_restriction(g_3D,op.eClosure,id_s)
383 e_n = BoundaryRestriction(g_3D,op.eClosure,id_n) 386 e_n = boundary_restriction(g_3D,op.eClosure,id_n)
384 e_b = BoundaryRestriction(g_3D,op.eClosure,id_b) 387 e_b = boundary_restriction(g_3D,op.eClosure,id_b)
385 e_t = BoundaryRestriction(g_3D,op.eClosure,id_t) 388 e_t = boundary_restriction(g_3D,op.eClosure,id_t)
386 e_dict = Dict(Pair(id_l,e_l),Pair(id_r,e_r), 389 e_dict = Dict(Pair(id_l,e_l),Pair(id_r,e_r),
387 Pair(id_s,e_s),Pair(id_n,e_n), 390 Pair(id_s,e_s),Pair(id_n,e_n),
388 Pair(id_b,e_b),Pair(id_t,e_t)) 391 Pair(id_b,e_b),Pair(id_t,e_t))
389 392
390 d_l = NormalDerivative(g_3D,op.dClosure,id_l) 393 d_l = normal_derivative(g_3D,op.dClosure,id_l)
391 d_r = NormalDerivative(g_3D,op.dClosure,id_r) 394 d_r = normal_derivative(g_3D,op.dClosure,id_r)
392 d_s = NormalDerivative(g_3D,op.dClosure,id_s) 395 d_s = normal_derivative(g_3D,op.dClosure,id_s)
393 d_n = NormalDerivative(g_3D,op.dClosure,id_n) 396 d_n = normal_derivative(g_3D,op.dClosure,id_n)
394 d_b = NormalDerivative(g_3D,op.dClosure,id_b) 397 d_b = normal_derivative(g_3D,op.dClosure,id_b)
395 d_t = NormalDerivative(g_3D,op.dClosure,id_t) 398 d_t = normal_derivative(g_3D,op.dClosure,id_t)
396 d_dict = Dict(Pair(id_l,d_l),Pair(id_r,d_r), 399 d_dict = Dict(Pair(id_l,d_l),Pair(id_r,d_r),
397 Pair(id_s,d_s),Pair(id_n,d_n), 400 Pair(id_s,d_s),Pair(id_n,d_n),
398 Pair(id_b,d_b),Pair(id_t,d_t)) 401 Pair(id_b,d_b),Pair(id_t,d_t))
399 402
400 H_l = boundary_quadrature(g_3D,op.quadratureClosure,id_r) 403 H_l = inner_product(boundary_grid(g_3D,id_l),op.quadratureClosure)
401 H_r = boundary_quadrature(g_3D,op.quadratureClosure,id_r) 404 H_r = inner_product(boundary_grid(g_3D,id_r),op.quadratureClosure)
402 H_s = boundary_quadrature(g_3D,op.quadratureClosure,id_s) 405 H_s = inner_product(boundary_grid(g_3D,id_s),op.quadratureClosure)
403 H_n = boundary_quadrature(g_3D,op.quadratureClosure,id_n) 406 H_n = inner_product(boundary_grid(g_3D,id_n),op.quadratureClosure)
404 H_b = boundary_quadrature(g_3D,op.quadratureClosure,id_b) 407 H_b = inner_product(boundary_grid(g_3D,id_b),op.quadratureClosure)
405 H_t = boundary_quadrature(g_3D,op.quadratureClosure,id_t) 408 H_t = inner_product(boundary_grid(g_3D,id_t),op.quadratureClosure)
406 Hb_dict = Dict(Pair(id_l,H_l),Pair(id_r,H_r), 409 Hb_dict = Dict(Pair(id_l,H_l),Pair(id_r,H_r),
407 Pair(id_s,H_s),Pair(id_n,H_n), 410 Pair(id_s,H_s),Pair(id_n,H_n),
408 Pair(id_b,H_b),Pair(id_t,H_t)) 411 Pair(id_b,H_b),Pair(id_t,H_t))
409 412
410 # TODO: Not sure why this doesnt work? Comparing the fields of 413 # TODO: Not sure why this doesnt work? Comparing the fields of
424 427
425 @testset "laplace" begin 428 @testset "laplace" begin
426 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 429 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
427 @testset "1D" begin 430 @testset "1D" begin
428 L = laplace(g_1D, op.innerStencil, op.closureStencils) 431 L = laplace(g_1D, op.innerStencil, op.closureStencils)
429 @test L == SecondDerivative(g_1D, op.innerStencil, op.closureStencils) 432 @test L == second_derivative(g_1D, op.innerStencil, op.closureStencils)
430 @test L isa TensorMapping{T,1,1} where T 433 @test L isa TensorMapping{T,1,1} where T
431 end 434 end
432 @testset "3D" begin 435 @testset "3D" begin
433 L = laplace(g_3D, op.innerStencil, op.closureStencils) 436 L = laplace(g_3D, op.innerStencil, op.closureStencils)
434 Dxx = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,1) 437 @test L isa TensorMapping{T,3,3} where T
435 Dyy = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,2) 438 Dxx = second_derivative(g_3D, op.innerStencil, op.closureStencils,1)
436 Dzz = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,3) 439 Dyy = second_derivative(g_3D, op.innerStencil, op.closureStencils,2)
440 Dzz = second_derivative(g_3D, op.innerStencil, op.closureStencils,3)
437 @test L == Dxx + Dyy + Dzz 441 @test L == Dxx + Dyy + Dzz
438 @test L isa TensorMapping{T,3,3} where T 442 @test L isa TensorMapping{T,3,3} where T
439 end 443 end
440 end 444 end
441 445
468 Δv = evalOn(g_3D,(x,y,z) -> -sin(x) - cos(y) + exp(z)) 472 Δv = evalOn(g_3D,(x,y,z) -> -sin(x) - cos(y) + exp(z))
469 473
470 # 2nd order interior stencil, 1st order boundary stencil, 474 # 2nd order interior stencil, 1st order boundary stencil,
471 # implies that L*v should be exact for binomials up to order 2. 475 # implies that L*v should be exact for binomials up to order 2.
472 @testset "2nd order" begin 476 @testset "2nd order" begin
473 L = Laplace(g_3D,sbp_operators_path()*"standard_diagonal.toml"; order=2) 477 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
478 L = laplace(g_3D,op.innerStencil,op.closureStencils)
474 @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9 479 @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
475 @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9 480 @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
476 @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9 481 @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9
477 @test L*v ≈ Δv rtol = 5e-2 norm = l2 482 @test L*v ≈ Δv rtol = 5e-2 norm = l2
478 end 483 end
479 484
480 # 4th order interior stencil, 2nd order boundary stencil, 485 # 4th order interior stencil, 2nd order boundary stencil,
481 # implies that L*v should be exact for binomials up to order 3. 486 # implies that L*v should be exact for binomials up to order 3.
482 @testset "4th order" begin 487 @testset "4th order" begin
483 L = Laplace(g_3D,sbp_operators_path()*"standard_diagonal.toml"; order=4) 488 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
489 L = laplace(g_3D,op.innerStencil,op.closureStencils)
484 # NOTE: high tolerances for checking the "exact" differentiation 490 # NOTE: high tolerances for checking the "exact" differentiation
485 # due to accumulation of round-off errors/cancellation errors? 491 # due to accumulation of round-off errors/cancellation errors?
486 @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9 492 @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
487 @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9 493 @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
488 @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9 494 @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9
490 @test L*v ≈ Δv rtol = 5e-4 norm = l2 496 @test L*v ≈ Δv rtol = 5e-4 norm = l2
491 end 497 end
492 end 498 end
493 end 499 end
494 500
495 @testset "Quadrature diagonal" begin 501 @testset "Diagonal-stencil inner_product" begin
496 Lx = π/2. 502 Lx = π/2.
497 Ly = Float64(π) 503 Ly = Float64(π)
498 Lz = 1. 504 Lz = 1.
499 g_1D = EquidistantGrid(77, 0.0, Lx) 505 g_1D = EquidistantGrid(77, 0.0, Lx)
500 g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) 506 g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
501 g_3D = EquidistantGrid((10,10, 10), (0.0, 0.0, 0.0), (Lx,Ly,Lz)) 507 g_3D = EquidistantGrid((10,10, 10), (0.0, 0.0, 0.0), (Lx,Ly,Lz))
502 integral(H,v) = sum(H*v) 508 integral(H,v) = sum(H*v)
503 @testset "quadrature" begin 509 @testset "inner_product" begin
504 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 510 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
505 @testset "1D" begin 511 @testset "0D" begin
506 H = quadrature(g_1D,op.quadratureClosure) 512 H = inner_product(EquidistantGrid{Float64}(),op.quadratureClosure)
507 inner_stencil = Stencil((1.,),center=1) 513 @test H == IdentityMapping{Float64}()
508 @test H == quadrature(g_1D,inner_stencil,op.quadratureClosure) 514 @test H isa TensorMapping{T,0,0} where T
515 end
516 @testset "1D" begin
517 H = inner_product(g_1D,op.quadratureClosure)
518 inner_stencil = CenteredStencil(1.)
519 @test H == inner_product(g_1D,op.quadratureClosure,inner_stencil)
509 @test H isa TensorMapping{T,1,1} where T 520 @test H isa TensorMapping{T,1,1} where T
510 end 521 end
511 @testset "2D" begin 522 @testset "2D" begin
512 H = quadrature(g_2D,op.quadratureClosure) 523 H = inner_product(g_2D,op.quadratureClosure)
513 H_x = quadrature(restrict(g_2D,1),op.quadratureClosure) 524 H_x = inner_product(restrict(g_2D,1),op.quadratureClosure)
514 H_y = quadrature(restrict(g_2D,2),op.quadratureClosure) 525 H_y = inner_product(restrict(g_2D,2),op.quadratureClosure)
515 @test H == H_x⊗H_y 526 @test H == H_x⊗H_y
516 @test H isa TensorMapping{T,2,2} where T 527 @test H isa TensorMapping{T,2,2} where T
517 end 528 end
518 end 529 end
519 530
520 @testset "boundary_quadrature" begin
521 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
522 @testset "1D" begin
523 (id_l, id_r) = boundary_identifiers(g_1D)
524 @test boundary_quadrature(g_1D,op.quadratureClosure,id_l) == IdentityMapping{Float64}()
525 @test boundary_quadrature(g_1D,op.quadratureClosure,id_r) == IdentityMapping{Float64}()
526
527 end
528 @testset "2D" begin
529 (id_w, id_e, id_s, id_n) = boundary_identifiers(g_2D)
530 H_x = quadrature(restrict(g_2D,1),op.quadratureClosure)
531 H_y = quadrature(restrict(g_2D,2),op.quadratureClosure)
532 @test boundary_quadrature(g_2D,op.quadratureClosure,id_w) == H_y
533 @test boundary_quadrature(g_2D,op.quadratureClosure,id_e) == H_y
534 @test boundary_quadrature(g_2D,op.quadratureClosure,id_s) == H_x
535 @test boundary_quadrature(g_2D,op.quadratureClosure,id_n) == H_x
536 end
537 @testset "3D" begin
538 (id_w, id_e,
539 id_s, id_n,
540 id_t, id_b) = boundary_identifiers(g_3D)
541 H_xy = quadrature(restrict(g_3D,[1,2]),op.quadratureClosure)
542 H_xz = quadrature(restrict(g_3D,[1,3]),op.quadratureClosure)
543 H_yz = quadrature(restrict(g_3D,[2,3]),op.quadratureClosure)
544 @test boundary_quadrature(g_3D,op.quadratureClosure,id_w) == H_yz
545 @test boundary_quadrature(g_3D,op.quadratureClosure,id_e) == H_yz
546 @test boundary_quadrature(g_3D,op.quadratureClosure,id_s) == H_xz
547 @test boundary_quadrature(g_3D,op.quadratureClosure,id_n) == H_xz
548 @test boundary_quadrature(g_3D,op.quadratureClosure,id_t) == H_xy
549 @test boundary_quadrature(g_3D,op.quadratureClosure,id_b) == H_xy
550 end
551 end
552
553 @testset "Sizes" begin 531 @testset "Sizes" begin
554 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 532 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
555 @testset "1D" begin 533 @testset "1D" begin
556 H = quadrature(g_1D,op.quadratureClosure) 534 H = inner_product(g_1D,op.quadratureClosure)
557 @test domain_size(H) == size(g_1D) 535 @test domain_size(H) == size(g_1D)
558 @test range_size(H) == size(g_1D) 536 @test range_size(H) == size(g_1D)
559 end 537 end
560 @testset "2D" begin 538 @testset "2D" begin
561 H = quadrature(g_2D,op.quadratureClosure) 539 H = inner_product(g_2D,op.quadratureClosure)
562 @test domain_size(H) == size(g_2D) 540 @test domain_size(H) == size(g_2D)
563 @test range_size(H) == size(g_2D) 541 @test range_size(H) == size(g_2D)
564 end 542 end
565 end 543 end
566 544
573 end 551 end
574 u = evalOn(g_1D,x->sin(x)) 552 u = evalOn(g_1D,x->sin(x))
575 553
576 @testset "2nd order" begin 554 @testset "2nd order" begin
577 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) 555 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
578 H = quadrature(g_1D,op.quadratureClosure) 556 H = inner_product(g_1D,op.quadratureClosure)
579 for i = 1:2 557 for i = 1:2
580 @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14 558 @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14
581 end 559 end
582 @test integral(H,u) ≈ 1. rtol = 1e-4 560 @test integral(H,u) ≈ 1. rtol = 1e-4
583 end 561 end
584 562
585 @testset "4th order" begin 563 @testset "4th order" begin
586 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 564 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
587 H = quadrature(g_1D,op.quadratureClosure) 565 H = inner_product(g_1D,op.quadratureClosure)
588 for i = 1:4 566 for i = 1:4
589 @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14 567 @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14
590 end 568 end
591 @test integral(H,u) ≈ 1. rtol = 1e-8 569 @test integral(H,u) ≈ 1. rtol = 1e-8
592 end 570 end
596 b = 2.1 574 b = 2.1
597 v = b*ones(Float64, size(g_2D)) 575 v = b*ones(Float64, size(g_2D))
598 u = evalOn(g_2D,(x,y)->sin(x)+cos(y)) 576 u = evalOn(g_2D,(x,y)->sin(x)+cos(y))
599 @testset "2nd order" begin 577 @testset "2nd order" begin
600 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) 578 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
601 H = quadrature(g_2D,op.quadratureClosure) 579 H = inner_product(g_2D,op.quadratureClosure)
602 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13 580 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13
603 @test integral(H,u) ≈ π rtol = 1e-4 581 @test integral(H,u) ≈ π rtol = 1e-4
604 end 582 end
605 @testset "4th order" begin 583 @testset "4th order" begin
606 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 584 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
607 H = quadrature(g_2D,op.quadratureClosure) 585 H = inner_product(g_2D,op.quadratureClosure)
608 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13 586 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13
609 @test integral(H,u) ≈ π rtol = 1e-8 587 @test integral(H,u) ≈ π rtol = 1e-8
610 end 588 end
611 end 589 end
612 end 590 end
613 end 591 end
614 592
615 @testset "InverseDiagonalQuadrature" begin 593 @testset "Diagonal-stencil inverse_inner_product" begin
616 Lx = π/2. 594 Lx = π/2.
617 Ly = Float64(π) 595 Ly = Float64(π)
618 g_1D = EquidistantGrid(77, 0.0, Lx) 596 g_1D = EquidistantGrid(77, 0.0, Lx)
619 g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) 597 g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
620 @testset "Constructors" begin 598 @testset "inverse_inner_product" begin
621 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 599 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
622 @testset "1D" begin 600 @testset "0D" begin
623 Hi = InverseDiagonalQuadrature(g_1D, op.quadratureClosure); 601 Hi = inverse_inner_product(EquidistantGrid{Float64}(),op.quadratureClosure)
624 inner_stencil = Stencil((1.,),center=1) 602 @test Hi == IdentityMapping{Float64}()
603 @test Hi isa TensorMapping{T,0,0} where T
604 end
605 @testset "1D" begin
606 Hi = inverse_inner_product(g_1D, op.quadratureClosure);
607 inner_stencil = CenteredStencil(1.)
625 closures = () 608 closures = ()
626 for i = 1:length(op.quadratureClosure) 609 for i = 1:length(op.quadratureClosure)
627 closures = (closures...,Stencil(op.quadratureClosure[i].range,1.0./op.quadratureClosure[i].weights)) 610 closures = (closures...,Stencil(op.quadratureClosure[i].range,1.0./op.quadratureClosure[i].weights))
628 end 611 end
629 @test Hi == InverseQuadrature(g_1D,inner_stencil,closures) 612 @test Hi == inverse_inner_product(g_1D,closures,inner_stencil)
630 @test Hi isa TensorMapping{T,1,1} where T 613 @test Hi isa TensorMapping{T,1,1} where T
631 end 614 end
632 @testset "2D" begin 615 @testset "2D" begin
633 Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure) 616 Hi = inverse_inner_product(g_2D,op.quadratureClosure)
634 Hi_x = InverseDiagonalQuadrature(restrict(g_2D,1),op.quadratureClosure) 617 Hi_x = inverse_inner_product(restrict(g_2D,1),op.quadratureClosure)
635 Hi_y = InverseDiagonalQuadrature(restrict(g_2D,2),op.quadratureClosure) 618 Hi_y = inverse_inner_product(restrict(g_2D,2),op.quadratureClosure)
636 @test Hi == Hi_x⊗Hi_y 619 @test Hi == Hi_x⊗Hi_y
637 @test Hi isa TensorMapping{T,2,2} where T 620 @test Hi isa TensorMapping{T,2,2} where T
638 end 621 end
639 end 622 end
640 623
641 @testset "Sizes" begin 624 @testset "Sizes" begin
642 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 625 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
643 @testset "1D" begin 626 @testset "1D" begin
644 Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure) 627 Hi = inverse_inner_product(g_1D,op.quadratureClosure)
645 @test domain_size(Hi) == size(g_1D) 628 @test domain_size(Hi) == size(g_1D)
646 @test range_size(Hi) == size(g_1D) 629 @test range_size(Hi) == size(g_1D)
647 end 630 end
648 @testset "2D" begin 631 @testset "2D" begin
649 Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure) 632 Hi = inverse_inner_product(g_2D,op.quadratureClosure)
650 @test domain_size(Hi) == size(g_2D) 633 @test domain_size(Hi) == size(g_2D)
651 @test range_size(Hi) == size(g_2D) 634 @test range_size(Hi) == size(g_2D)
652 end 635 end
653 end 636 end
654 637
656 @testset "1D" begin 639 @testset "1D" begin
657 v = evalOn(g_1D,x->sin(x)) 640 v = evalOn(g_1D,x->sin(x))
658 u = evalOn(g_1D,x->x^3-x^2+1) 641 u = evalOn(g_1D,x->x^3-x^2+1)
659 @testset "2nd order" begin 642 @testset "2nd order" begin
660 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) 643 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
661 H = quadrature(g_1D,op.quadratureClosure) 644 H = inner_product(g_1D,op.quadratureClosure)
662 Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure) 645 Hi = inverse_inner_product(g_1D,op.quadratureClosure)
663 @test Hi*H*v ≈ v rtol = 1e-15 646 @test Hi*H*v ≈ v rtol = 1e-15
664 @test Hi*H*u ≈ u rtol = 1e-15 647 @test Hi*H*u ≈ u rtol = 1e-15
665 end 648 end
666 @testset "4th order" begin 649 @testset "4th order" begin
667 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 650 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
668 H = quadrature(g_1D,op.quadratureClosure) 651 H = inner_product(g_1D,op.quadratureClosure)
669 Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure) 652 Hi = inverse_inner_product(g_1D,op.quadratureClosure)
670 @test Hi*H*v ≈ v rtol = 1e-15 653 @test Hi*H*v ≈ v rtol = 1e-15
671 @test Hi*H*u ≈ u rtol = 1e-15 654 @test Hi*H*u ≈ u rtol = 1e-15
672 end 655 end
673 end 656 end
674 @testset "2D" begin 657 @testset "2D" begin
675 v = evalOn(g_2D,(x,y)->sin(x)+cos(y)) 658 v = evalOn(g_2D,(x,y)->sin(x)+cos(y))
676 u = evalOn(g_2D,(x,y)->x*y + x^5 - sqrt(y)) 659 u = evalOn(g_2D,(x,y)->x*y + x^5 - sqrt(y))
677 @testset "2nd order" begin 660 @testset "2nd order" begin
678 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) 661 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
679 H = quadrature(g_2D,op.quadratureClosure) 662 H = inner_product(g_2D,op.quadratureClosure)
680 Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure) 663 Hi = inverse_inner_product(g_2D,op.quadratureClosure)
681 @test Hi*H*v ≈ v rtol = 1e-15 664 @test Hi*H*v ≈ v rtol = 1e-15
682 @test Hi*H*u ≈ u rtol = 1e-15 665 @test Hi*H*u ≈ u rtol = 1e-15
683 end 666 end
684 @testset "4th order" begin 667 @testset "4th order" begin
685 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 668 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
686 H = quadrature(g_2D,op.quadratureClosure) 669 H = inner_product(g_2D,op.quadratureClosure)
687 Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure) 670 Hi = inverse_inner_product(g_2D,op.quadratureClosure)
688 @test Hi*H*v ≈ v rtol = 1e-15 671 @test Hi*H*v ≈ v rtol = 1e-15
689 @test Hi*H*u ≈ u rtol = 1e-15 672 @test Hi*H*u ≈ u rtol = 1e-15
690 end 673 end
691 end 674 end
692 end 675 end
703 @test op_l == BoundaryOperator(g_1D,closure_stencil,Lower()) 686 @test op_l == BoundaryOperator(g_1D,closure_stencil,Lower())
704 @test op_l == boundary_operator(g_1D,closure_stencil,CartesianBoundary{1,Lower}()) 687 @test op_l == boundary_operator(g_1D,closure_stencil,CartesianBoundary{1,Lower}())
705 @test op_l isa TensorMapping{T,0,1} where T 688 @test op_l isa TensorMapping{T,0,1} where T
706 689
707 op_r = BoundaryOperator{Upper}(closure_stencil,size(g_1D)[1]) 690 op_r = BoundaryOperator{Upper}(closure_stencil,size(g_1D)[1])
708 @test op_r == BoundaryRestriction(g_1D,closure_stencil,Upper()) 691 @test op_r == BoundaryOperator(g_1D,closure_stencil,Upper())
709 @test op_r == boundary_operator(g_1D,closure_stencil,CartesianBoundary{1,Upper}()) 692 @test op_r == boundary_operator(g_1D,closure_stencil,CartesianBoundary{1,Upper}())
710 @test op_r isa TensorMapping{T,0,1} where T 693 @test op_r isa TensorMapping{T,0,1} where T
711 end 694 end
712 695
713 @testset "2D" begin 696 @testset "2D" begin
834 @inferred apply_transpose(op_r, u, Index(11,Upper)) 817 @inferred apply_transpose(op_r, u, Index(11,Upper))
835 end 818 end
836 819
837 end 820 end
838 821
839 @testset "BoundaryRestriction" begin 822 @testset "boundary_restriction" begin
840 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 823 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
841 g_1D = EquidistantGrid(11, 0.0, 1.0) 824 g_1D = EquidistantGrid(11, 0.0, 1.0)
842 g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0)) 825 g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))
843 826
844 @testset "Constructors" begin 827 @testset "boundary_restriction" begin
845 @testset "1D" begin 828 @testset "1D" begin
846 e_l = BoundaryRestriction(g_1D,op.eClosure,Lower()) 829 e_l = boundary_restriction(g_1D,op.eClosure,Lower())
847 @test e_l == BoundaryRestriction(g_1D,op.eClosure,CartesianBoundary{1,Lower}()) 830 @test e_l == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Lower}())
848 @test e_l == BoundaryOperator(g_1D,op.eClosure,Lower()) 831 @test e_l == BoundaryOperator(g_1D,op.eClosure,Lower())
849 @test e_l isa BoundaryOperator{T,Lower} where T 832 @test e_l isa BoundaryOperator{T,Lower} where T
850 @test e_l isa TensorMapping{T,0,1} where T 833 @test e_l isa TensorMapping{T,0,1} where T
851 834
852 e_r = BoundaryRestriction(g_1D,op.eClosure,Upper()) 835 e_r = boundary_restriction(g_1D,op.eClosure,Upper())
853 @test e_r == BoundaryRestriction(g_1D,op.eClosure,CartesianBoundary{1,Upper}()) 836 @test e_r == boundary_restriction(g_1D,op.eClosure,CartesianBoundary{1,Upper}())
854 @test e_r == BoundaryOperator(g_1D,op.eClosure,Upper()) 837 @test e_r == BoundaryOperator(g_1D,op.eClosure,Upper())
855 @test e_r isa BoundaryOperator{T,Upper} where T 838 @test e_r isa BoundaryOperator{T,Upper} where T
856 @test e_r isa TensorMapping{T,0,1} where T 839 @test e_r isa TensorMapping{T,0,1} where T
857 end 840 end
858 841
859 @testset "2D" begin 842 @testset "2D" begin
860 e_w = BoundaryRestriction(g_2D,op.eClosure,CartesianBoundary{1,Upper}()) 843 e_w = boundary_restriction(g_2D,op.eClosure,CartesianBoundary{1,Upper}())
861 @test e_w isa InflatedTensorMapping 844 @test e_w isa InflatedTensorMapping
862 @test e_w isa TensorMapping{T,1,2} where T 845 @test e_w isa TensorMapping{T,1,2} where T
863 end 846 end
864 end 847 end
865 848
866 @testset "Application" begin 849 @testset "Application" begin
867 @testset "1D" begin 850 @testset "1D" begin
868 e_l = BoundaryRestriction(g_1D, op.eClosure, CartesianBoundary{1,Lower}()) 851 e_l = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Lower}())
869 e_r = BoundaryRestriction(g_1D, op.eClosure, CartesianBoundary{1,Upper}()) 852 e_r = boundary_restriction(g_1D, op.eClosure, CartesianBoundary{1,Upper}())
870 853
871 v = evalOn(g_1D,x->1+x^2) 854 v = evalOn(g_1D,x->1+x^2)
872 u = fill(3.124) 855 u = fill(3.124)
873 856
874 @test (e_l*v)[] == v[1] 857 @test (e_l*v)[] == v[1]
875 @test (e_r*v)[] == v[end] 858 @test (e_r*v)[] == v[end]
876 @test (e_r*v)[1] == v[end] 859 @test (e_r*v)[1] == v[end]
877 end 860 end
878 861
879 @testset "2D" begin 862 @testset "2D" begin
880 e_w = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{1,Lower}()) 863 e_w = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{1,Lower}())
881 e_e = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{1,Upper}()) 864 e_e = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{1,Upper}())
882 e_s = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{2,Lower}()) 865 e_s = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{2,Lower}())
883 e_n = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{2,Upper}()) 866 e_n = boundary_restriction(g_2D, op.eClosure, CartesianBoundary{2,Upper}())
884 867
885 v = rand(11, 15) 868 v = rand(11, 15)
886 u = fill(3.124) 869 u = fill(3.124)
887 870
888 @test e_w*v == v[1,:] 871 @test e_w*v == v[1,:]
891 @test e_n*v == v[:,end] 874 @test e_n*v == v[:,end]
892 end 875 end
893 end 876 end
894 end 877 end
895 878
896 @testset "NormalDerivative" begin 879 @testset "normal_derivative" begin
897 g_1D = EquidistantGrid(11, 0.0, 1.0) 880 g_1D = EquidistantGrid(11, 0.0, 1.0)
898 g_2D = EquidistantGrid((11,12), (0.0, 0.0), (1.0,1.0)) 881 g_2D = EquidistantGrid((11,12), (0.0, 0.0), (1.0,1.0))
899 @testset "Constructors" begin 882 @testset "normal_derivative" begin
900 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 883 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
901 @testset "1D" begin 884 @testset "1D" begin
902 d_l = NormalDerivative(g_1D, op.dClosure, Lower()) 885 d_l = normal_derivative(g_1D, op.dClosure, Lower())
903 @test d_l == NormalDerivative(g_1D, op.dClosure, CartesianBoundary{1,Lower}()) 886 @test d_l == normal_derivative(g_1D, op.dClosure, CartesianBoundary{1,Lower}())
904 @test d_l isa BoundaryOperator{T,Lower} where T 887 @test d_l isa BoundaryOperator{T,Lower} where T
905 @test d_l isa TensorMapping{T,0,1} where T 888 @test d_l isa TensorMapping{T,0,1} where T
906 end 889 end
907 @testset "2D" begin 890 @testset "2D" begin
908 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 891 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
909 d_w = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}()) 892 d_w = normal_derivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}())
910 d_n = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}()) 893 d_n = normal_derivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}())
911 Ix = IdentityMapping{Float64}((size(g_2D)[1],)) 894 Ix = IdentityMapping{Float64}((size(g_2D)[1],))
912 Iy = IdentityMapping{Float64}((size(g_2D)[2],)) 895 Iy = IdentityMapping{Float64}((size(g_2D)[2],))
913 d_l = NormalDerivative(restrict(g_2D,1),op.dClosure,Lower()) 896 d_l = normal_derivative(restrict(g_2D,1),op.dClosure,Lower())
914 d_r = NormalDerivative(restrict(g_2D,2),op.dClosure,Upper()) 897 d_r = normal_derivative(restrict(g_2D,2),op.dClosure,Upper())
915 @test d_w == d_l⊗Iy 898 @test d_w == d_l⊗Iy
916 @test d_n == Ix⊗d_r 899 @test d_n == Ix⊗d_r
917 @test d_w isa TensorMapping{T,1,2} where T 900 @test d_w isa TensorMapping{T,1,2} where T
918 @test d_n isa TensorMapping{T,1,2} where T 901 @test d_n isa TensorMapping{T,1,2} where T
919 end 902 end
923 v∂x = evalOn(g_2D, (x,y)-> 2*x + y) 906 v∂x = evalOn(g_2D, (x,y)-> 2*x + y)
924 v∂y = evalOn(g_2D, (x,y)-> 2*(y-1) + x) 907 v∂y = evalOn(g_2D, (x,y)-> 2*(y-1) + x)
925 # TODO: Test for higher order polynomials? 908 # TODO: Test for higher order polynomials?
926 @testset "2nd order" begin 909 @testset "2nd order" begin
927 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) 910 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
928 d_w = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}()) 911 d_w = normal_derivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}())
929 d_e = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{1,Upper}()) 912 d_e = normal_derivative(g_2D, op.dClosure, CartesianBoundary{1,Upper}())
930 d_s = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{2,Lower}()) 913 d_s = normal_derivative(g_2D, op.dClosure, CartesianBoundary{2,Lower}())
931 d_n = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}()) 914 d_n = normal_derivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}())
932 915
933 @test d_w*v ≈ v∂x[1,:] atol = 1e-13 916 @test d_w*v ≈ v∂x[1,:] atol = 1e-13
934 @test d_e*v ≈ -v∂x[end,:] atol = 1e-13 917 @test d_e*v ≈ -v∂x[end,:] atol = 1e-13
935 @test d_s*v ≈ v∂y[:,1] atol = 1e-13 918 @test d_s*v ≈ v∂y[:,1] atol = 1e-13
936 @test d_n*v ≈ -v∂y[:,end] atol = 1e-13 919 @test d_n*v ≈ -v∂y[:,end] atol = 1e-13
937 end 920 end
938 921
939 @testset "4th order" begin 922 @testset "4th order" begin
940 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 923 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
941 d_w = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}()) 924 d_w = normal_derivative(g_2D, op.dClosure, CartesianBoundary{1,Lower}())
942 d_e = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{1,Upper}()) 925 d_e = normal_derivative(g_2D, op.dClosure, CartesianBoundary{1,Upper}())
943 d_s = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{2,Lower}()) 926 d_s = normal_derivative(g_2D, op.dClosure, CartesianBoundary{2,Lower}())
944 d_n = NormalDerivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}()) 927 d_n = normal_derivative(g_2D, op.dClosure, CartesianBoundary{2,Upper}())
945 928
946 @test d_w*v ≈ v∂x[1,:] atol = 1e-13 929 @test d_w*v ≈ v∂x[1,:] atol = 1e-13
947 @test d_e*v ≈ -v∂x[end,:] atol = 1e-13 930 @test d_e*v ≈ -v∂x[end,:] atol = 1e-13
948 @test d_s*v ≈ v∂y[:,1] atol = 1e-13 931 @test d_s*v ≈ v∂y[:,1] atol = 1e-13
949 @test d_n*v ≈ -v∂y[:,end] atol = 1e-13 932 @test d_n*v ≈ -v∂y[:,end] atol = 1e-13