Mercurial > repos > public > sbplib_julia

comparison test/testSbpOperators.jl @ 617:f59e1732eacc feature/volume_and_boundary_operators

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 07 Dec 2020 12:07:29 +0100
parents 1db945cba3a2 eaa8c852ddf2
children 332f65c1abf3
comparison
equal deleted inserted replaced
616:d9324671b412 617:f59e1732eacc
2 using Sbplib.SbpOperators 2 using Sbplib.SbpOperators
3 using Sbplib.Grids 3 using Sbplib.Grids
4 using Sbplib.RegionIndices 4 using Sbplib.RegionIndices
5 using Sbplib.LazyTensors 5 using Sbplib.LazyTensors
6 using LinearAlgebra 6 using LinearAlgebra
7 using TOML
8
9 import Sbplib.SbpOperators.Stencil
7 10
8 @testset "SbpOperators" begin 11 @testset "SbpOperators" begin
9 12
10 @testset "Stencil" begin 13 @testset "Stencil" begin
11 s = SbpOperators.Stencil((-2,2), (1.,2.,2.,3.,4.)) 14 s = Stencil((-2,2), (1.,2.,2.,3.,4.))
12 @test s isa SbpOperators.Stencil{Float64, 5} 15 @test s isa Stencil{Float64, 5}
13 16
14 @test eltype(s) == Float64 17 @test eltype(s) == Float64
15 @test SbpOperators.scale(s, 2) == SbpOperators.Stencil((-2,2), (2.,4.,4.,6.,8.)) 18 @test SbpOperators.scale(s, 2) == Stencil((-2,2), (2.,4.,4.,6.,8.))
19
20 @test Stencil((1,2,3,4), center=1) == Stencil((0, 3),(1,2,3,4))
21 @test Stencil((1,2,3,4), center=2) == Stencil((-1, 2),(1,2,3,4))
22 @test Stencil((1,2,3,4), center=4) == Stencil((-3, 0),(1,2,3,4))
23 end
24
25 @testset "parse_rational" begin
26 @test SbpOperators.parse_rational("1") isa Rational
27 @test SbpOperators.parse_rational("1") == 1//1
28 @test SbpOperators.parse_rational("1/2") isa Rational
29 @test SbpOperators.parse_rational("1/2") == 1//2
30 @test SbpOperators.parse_rational("37/13") isa Rational
31 @test SbpOperators.parse_rational("37/13") == 37//13
32 end
33
34 @testset "readoperator" begin
35 toml_str = """
36 [meta]
37 type = "equidistant"
38
39 [order2]
40 H.inner = ["1"]
41
42 D1.inner_stencil = ["-1/2", "0", "1/2"]
43 D1.closure_stencils = [
44 ["-1", "1"],
45 ]
46
47 d1.closure = ["-3/2", "2", "-1/2"]
48
49 [order4]
50 H.closure = ["17/48", "59/48", "43/48", "49/48"]
51
52 D2.inner_stencil = ["-1/12","4/3","-5/2","4/3","-1/12"]
53 D2.closure_stencils = [
54 [ "2", "-5", "4", "-1", "0", "0"],
55 [ "1", "-2", "1", "0", "0", "0"],
56 [ "-4/43", "59/43", "-110/43", "59/43", "-4/43", "0"],
57 [ "-1/49", "0", "59/49", "-118/49", "64/49", "-4/49"],
58 ]
59 """
60
61 parsed_toml = TOML.parse(toml_str)
62 @testset "get_stencil" begin
63 @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil") == Stencil((-1/2, 0., 1/2), center=2)
64 @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil", center=1) == Stencil((-1/2, 0., 1/2); center=1)
65 @test get_stencil(parsed_toml, "order2", "D1", "inner_stencil", center=3) == Stencil((-1/2, 0., 1/2); center=3)
66
67 @test get_stencil(parsed_toml, "order2", "H", "inner") == Stencil((1.,), center=1)
68
69 @test_throws AssertionError get_stencil(parsed_toml, "meta", "type")
70 @test_throws AssertionError get_stencil(parsed_toml, "order2", "D1", "closure_stencils")
71 end
72
73 @testset "get_stencils" begin
74 @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=(1,)) == (Stencil((-1., 1.), center=1),)
75 @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=(2,)) == (Stencil((-1., 1.), center=2),)
76 @test get_stencils(parsed_toml, "order2", "D1", "closure_stencils", centers=[2]) == (Stencil((-1., 1.), center=2),)
77
78 @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=[1,1,1,1]) == (
79 Stencil(( 2., -5., 4., -1., 0., 0.), center=1),
80 Stencil(( 1., -2., 1., 0., 0., 0.), center=1),
81 Stencil(( -4/43, 59/43, -110/43, 59/43, -4/43, 0.), center=1),
82 Stencil(( -1/49, 0., 59/49, -118/49, 64/49, -4/49), center=1),
83 )
84
85 @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(4,2,3,1)) == (
86 Stencil(( 2., -5., 4., -1., 0., 0.), center=4),
87 Stencil(( 1., -2., 1., 0., 0., 0.), center=2),
88 Stencil(( -4/43, 59/43, -110/43, 59/43, -4/43, 0.), center=3),
89 Stencil(( -1/49, 0., 59/49, -118/49, 64/49, -4/49), center=1),
90 )
91
92 @test get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=1:4) == (
93 Stencil(( 2., -5., 4., -1., 0., 0.), center=1),
94 Stencil(( 1., -2., 1., 0., 0., 0.), center=2),
95 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=4),
97 )
98
99 @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(1,2,3))
100 @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "closure_stencils",centers=(1,2,3,5,4))
101 @test_throws AssertionError get_stencils(parsed_toml, "order4", "D2", "inner_stencil",centers=(1,2))
102 end
103
104 @testset "get_tuple" begin
105 @test get_tuple(parsed_toml, "order2", "d1", "closure") == (-3/2, 2, -1/2)
106
107 @test_throws AssertionError get_tuple(parsed_toml, "meta", "type")
108 end
16 end 109 end
17 110
18 # @testset "apply_quadrature" begin 111 # @testset "apply_quadrature" begin
19 # op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 112 # op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
20 # h = 0.5 113 # h = 0.5
21 # 114 #
22 # @test apply_quadrature(op, h, 1.0, 10, 100) == h 115 # @test apply_quadrature(op, h, 1.0, 10, 100) == h
23 # 116 #
24 # N = 10 117 # N = 10
35 # @test apply_quadrature(op, h, v[i], i, N) == q[i]*v[i] 128 # @test apply_quadrature(op, h, v[i], i, N) == q[i]*v[i]
36 # end 129 # end
37 # end 130 # end
38 131
39 @testset "SecondDerivative" begin 132 @testset "SecondDerivative" begin
40 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 133 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
41 L = 3.5 134 L = 3.5
42 g = EquidistantGrid(101, 0.0, L) 135 g = EquidistantGrid(101, 0.0, L)
43 Dₓₓ = SecondDerivative(g,op.innerStencil,op.closureStencils) 136 Dₓₓ = SecondDerivative(g,op.innerStencil,op.closureStencils)
44 137
45 f0(x) = 1. 138 f0(x) = 1.
75 @test Dₓₓ*v5 ≈ -v5 atol=5e-4 norm=l2 168 @test Dₓₓ*v5 ≈ -v5 atol=5e-4 norm=l2
76 end 169 end
77 170
78 171
79 @testset "Laplace2D" begin 172 @testset "Laplace2D" begin
80 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 173 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
81 Lx = 1.5 174 Lx = 1.5
82 Ly = 3.2 175 Ly = 3.2
83 g = EquidistantGrid((102,131), (0.0, 0.0), (Lx,Ly)) 176 g = EquidistantGrid((102,131), (0.0, 0.0), (Lx,Ly))
84 L = Laplace(g, op.innerStencil, op.closureStencils) 177 L = Laplace(g, op.innerStencil, op.closureStencils)
85 178
117 @test L*v4 ≈ v2 atol=5e-4 norm=l2 210 @test L*v4 ≈ v2 atol=5e-4 norm=l2
118 @test L*v5 ≈ v5ₓₓ atol=5e-4 norm=l2 211 @test L*v5 ≈ v5ₓₓ atol=5e-4 norm=l2
119 end 212 end
120 213
121 @testset "DiagonalInnerProduct" begin 214 @testset "DiagonalInnerProduct" begin
122 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 215 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
123 L = 2.3 216 L = 2.3
124 g = EquidistantGrid(77, 0.0, L) 217 g = EquidistantGrid(77, 0.0, L)
125 H = DiagonalInnerProduct(g,op.quadratureClosure) 218 H = DiagonalInnerProduct(g,op.quadratureClosure)
126 v = ones(Float64, size(g)) 219 v = ones(Float64, size(g))
127 220
130 @test sum(H*v) ≈ L 223 @test sum(H*v) ≈ L
131 @test H*v == H'*v 224 @test H*v == H'*v
132 end 225 end
133 226
134 @testset "Quadrature" begin 227 @testset "Quadrature" begin
135 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 228 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
136 Lx = 2.3 229 Lx = 2.3
137 Ly = 5.2 230 Ly = 5.2
138 g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) 231 g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
139 232
140 Q = Quadrature(g, op.quadratureClosure) 233 Q = Quadrature(g, op.quadratureClosure)
150 243
151 @test Q*v == Q'*v 244 @test Q*v == Q'*v
152 end 245 end
153 246
154 @testset "InverseDiagonalInnerProduct" begin 247 @testset "InverseDiagonalInnerProduct" begin
155 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 248 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
156 L = 2.3 249 L = 2.3
157 g = EquidistantGrid(77, 0.0, L) 250 g = EquidistantGrid(77, 0.0, L)
158 H = DiagonalInnerProduct(g, op.quadratureClosure) 251 H = DiagonalInnerProduct(g, op.quadratureClosure)
159 Hi = InverseDiagonalInnerProduct(g,op.quadratureClosure) 252 Hi = InverseDiagonalInnerProduct(g,op.quadratureClosure)
160 v = evalOn(g, x->sin(x)) 253 v = evalOn(g, x->sin(x))
164 @test Hi*H*v ≈ v 257 @test Hi*H*v ≈ v
165 @test Hi*v == Hi'*v 258 @test Hi*v == Hi'*v
166 end 259 end
167 260
168 @testset "InverseQuadrature" begin 261 @testset "InverseQuadrature" begin
169 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 262 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
170 Lx = 7.3 263 Lx = 7.3
171 Ly = 8.2 264 Ly = 8.2
172 g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) 265 g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
173 266
174 Q = Quadrature(g, op.quadratureClosure) 267 Q = Quadrature(g, op.quadratureClosure)
180 @test_broken Qinv*(Q*v) ≈ v 273 @test_broken Qinv*(Q*v) ≈ v
181 @test Qinv*v == Qinv'*v 274 @test Qinv*v == Qinv'*v
182 end 275 end
183 276
184 @testset "BoundaryOperator" begin 277 @testset "BoundaryOperator" begin
185 closure_stencil = SbpOperators.Stencil((0,2), (2.,1.,3.)) 278 closure_stencil = Stencil((0,2), (2.,1.,3.))
186 g_1D = EquidistantGrid(11, 0.0, 1.0) 279 g_1D = EquidistantGrid(11, 0.0, 1.0)
187 g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0)) 280 g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))
188 281
189 @testset "Constructors" begin 282 @testset "Constructors" begin
190 @testset "1D" begin 283 @testset "1D" begin
325 end 418 end
326 419
327 end 420 end
328 421
329 @testset "BoundaryRestriction" begin 422 @testset "BoundaryRestriction" begin
330 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 423 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
331 g_1D = EquidistantGrid(11, 0.0, 1.0) 424 g_1D = EquidistantGrid(11, 0.0, 1.0)
332 g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0)) 425 g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))
333 426
334 @testset "Constructors" begin 427 @testset "Constructors" begin
335 @testset "1D" begin 428 @testset "1D" begin
405 end 498 end
406 end 499 end
407 end 500 end
408 # 501 #
409 # @testset "NormalDerivative" begin 502 # @testset "NormalDerivative" begin
410 # op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 503 # op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
411 # g = EquidistantGrid((5,6), (0.0, 0.0), (4.0,5.0)) 504 # g = EquidistantGrid((5,6), (0.0, 0.0), (4.0,5.0))
412 # 505 #
413 # d_w = NormalDerivative(op, g, CartesianBoundary{1,Lower}()) 506 # d_w = NormalDerivative(op, g, CartesianBoundary{1,Lower}())
414 # d_e = NormalDerivative(op, g, CartesianBoundary{1,Upper}()) 507 # d_e = NormalDerivative(op, g, CartesianBoundary{1,Upper}())
415 # d_s = NormalDerivative(op, g, CartesianBoundary{2,Lower}()) 508 # d_s = NormalDerivative(op, g, CartesianBoundary{2,Lower}())
482 # @test_broken d_s*g_x .≈ G_s 575 # @test_broken d_s*g_x .≈ G_s
483 # @test_broken d_n*g_x .≈ G_n 576 # @test_broken d_n*g_x .≈ G_n
484 # end 577 # end
485 # 578 #
486 # @testset "BoundaryQuadrature" begin 579 # @testset "BoundaryQuadrature" begin
487 # op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 580 # op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
488 # g = EquidistantGrid((10,11), (0.0, 0.0), (1.0,1.0)) 581 # g = EquidistantGrid((10,11), (0.0, 0.0), (1.0,1.0))
489 # 582 #
490 # H_w = BoundaryQuadrature(op, g, CartesianBoundary{1,Lower}()) 583 # H_w = BoundaryQuadrature(op, g, CartesianBoundary{1,Lower}())
491 # H_e = BoundaryQuadrature(op, g, CartesianBoundary{1,Upper}()) 584 # H_e = BoundaryQuadrature(op, g, CartesianBoundary{1,Upper}())
492 # H_s = BoundaryQuadrature(op, g, CartesianBoundary{2,Lower}()) 585 # H_s = BoundaryQuadrature(op, g, CartesianBoundary{2,Lower}())