Mercurial > repos > public > sbplib_julia

comparison test/SbpOperators/volumeops/volume_operator_test.jl @ 1858:4a9be96f2569 feature/documenter_logo

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge default
author Jonatan Werpers <jonatan@werpers.com>
date Sun, 12 Jan 2025 21:18:44 +0100
parents 0656b46a1a74
children
comparison
equal deleted inserted replaced
1857:ffde7dad9da5 1858:4a9be96f2569
1 using Test 1 using Test
2 2
3 using Sbplib.SbpOperators 3 using Diffinitive.SbpOperators
4 using Sbplib.Grids 4 using Diffinitive.Grids
5 using Sbplib.RegionIndices 5 using Diffinitive.RegionIndices
6 using Sbplib.LazyTensors 6 using Diffinitive.LazyTensors
7 7
8 import Sbplib.SbpOperators.Stencil 8 import Diffinitive.SbpOperators.Stencil
9 import Sbplib.SbpOperators.VolumeOperator 9 import Diffinitive.SbpOperators.VolumeOperator
10 import Sbplib.SbpOperators.volume_operator 10 import Diffinitive.SbpOperators.odd
11 import Sbplib.SbpOperators.odd 11 import Diffinitive.SbpOperators.even
12 import Sbplib.SbpOperators.even 12
13 13
14 @testset "VolumeOperator" begin 14 @testset "VolumeOperator" begin
15 inner_stencil = CenteredStencil(1/4, 2/4, 1/4) 15 inner_stencil = CenteredStencil(1/4, 2/4, 1/4)
16 closure_stencils = (Stencil(1/2, 1/2; center=1), Stencil(0.,1.; center=2)) 16 closure_stencils = (Stencil(1/2, 1/2; center=1), Stencil(2.,1.; center=2))
17 g_1D = EquidistantGrid(11,0.,1.) 17 g = equidistant_grid(0.,1., 11)
18 g_2D = EquidistantGrid((11,12),(0.,0.),(1.,1.)) 18
19 g_3D = EquidistantGrid((11,12,10),(0.,0.,0.),(1.,1.,1.))
20 @testset "Constructors" begin 19 @testset "Constructors" begin
21 @testset "1D" begin 20 op = VolumeOperator(inner_stencil, closure_stencils, 11, even)
22 op = VolumeOperator(inner_stencil,closure_stencils,(11,),even) 21 @test op == VolumeOperator(g,inner_stencil,closure_stencils,even)
23 @test op == VolumeOperator(g_1D,inner_stencil,closure_stencils,even) 22 @test op isa LazyTensor{T,1,1} where T
24 @test op == volume_operator(g_1D,inner_stencil,closure_stencils,even,1)
25 @test op isa TensorMapping{T,1,1} where T
26 end
27 @testset "2D" begin
28 op_x = volume_operator(g_2D,inner_stencil,closure_stencils,even,1)
29 op_y = volume_operator(g_2D,inner_stencil,closure_stencils,even,2)
30 Ix = IdentityMapping{Float64}((11,))
31 Iy = IdentityMapping{Float64}((12,))
32 @test op_x == VolumeOperator(inner_stencil,closure_stencils,(11,),even)⊗Iy
33 @test op_y == Ix⊗VolumeOperator(inner_stencil,closure_stencils,(12,),even)
34 @test op_x isa TensorMapping{T,2,2} where T
35 @test op_y isa TensorMapping{T,2,2} where T
36 end
37 @testset "3D" begin
38 op_x = volume_operator(g_3D,inner_stencil,closure_stencils,even,1)
39 op_y = volume_operator(g_3D,inner_stencil,closure_stencils,even,2)
40 op_z = volume_operator(g_3D,inner_stencil,closure_stencils,even,3)
41 Ix = IdentityMapping{Float64}((11,))
42 Iy = IdentityMapping{Float64}((12,))
43 Iz = IdentityMapping{Float64}((10,))
44 @test op_x == VolumeOperator(inner_stencil,closure_stencils,(11,),even)⊗Iy⊗Iz
45 @test op_y == Ix⊗VolumeOperator(inner_stencil,closure_stencils,(12,),even)⊗Iz
46 @test op_z == Ix⊗Iy⊗VolumeOperator(inner_stencil,closure_stencils,(10,),even)
47 @test op_x isa TensorMapping{T,3,3} where T
48 @test op_y isa TensorMapping{T,3,3} where T
49 @test op_z isa TensorMapping{T,3,3} where T
50 end
51 end 23 end
52 24
53 @testset "Sizes" begin 25 @testset "Sizes" begin
54 @testset "1D" begin 26 op = VolumeOperator(g,inner_stencil,closure_stencils,even)
55 op = volume_operator(g_1D,inner_stencil,closure_stencils,even,1) 27 @test range_size(op) == domain_size(op) == size(g)
56 @test range_size(op) == domain_size(op) == size(g_1D)
57 end
58
59 @testset "2D" begin
60 op_x = volume_operator(g_2D,inner_stencil,closure_stencils,even,1)
61 op_y = volume_operator(g_2D,inner_stencil,closure_stencils,even,2)
62 @test range_size(op_y) == domain_size(op_y) ==
63 range_size(op_x) == domain_size(op_x) == size(g_2D)
64 end
65 @testset "3D" begin
66 op_x = volume_operator(g_3D,inner_stencil,closure_stencils,even,1)
67 op_y = volume_operator(g_3D,inner_stencil,closure_stencils,even,2)
68 op_z = volume_operator(g_3D,inner_stencil,closure_stencils,even,3)
69 @test range_size(op_z) == domain_size(op_z) ==
70 range_size(op_y) == domain_size(op_y) ==
71 range_size(op_x) == domain_size(op_x) == size(g_3D)
72 end
73 end 28 end
74 29
75 op_x = volume_operator(g_2D,inner_stencil,closure_stencils,even,1) 30
76 op_y = volume_operator(g_2D,inner_stencil,closure_stencils,odd,2) 31 op_even = VolumeOperator(g, inner_stencil, closure_stencils, even)
77 v = zeros(size(g_2D)) 32 op_odd = VolumeOperator(g, inner_stencil, closure_stencils, odd)
78 Nx = size(g_2D)[1] 33
79 Ny = size(g_2D)[2] 34 N = size(g)[1]
80 for i = 1:Nx 35 v = rand(N)
81 v[i,:] .= i 36
37 r_even = copy(v)
38 r_odd = copy(v)
39
40 r_even[1] = (v[1] + v[2])/2
41 r_odd[1] = (v[1] + v[2])/2
42
43 r_even[2] = 2v[1] + v[2]
44 r_odd[2] = 2v[1] + v[2]
45
46 for i ∈ 3:N-2
47 r_even[i] = (v[i-1] + 2v[i] + v[i+1])/4
48 r_odd[i] = (v[i-1] + 2v[i] + v[i+1])/4
82 end 49 end
83 rx = copy(v) 50
84 rx[1,:] .= 1.5 51 r_even[N-1] = v[N-1] + 2v[N]
85 rx[Nx,:] .= (2*Nx-1)/2 52 r_odd[N-1] = -v[N-1] - 2v[N]
86 ry = copy(v) 53
87 ry[:,Ny-1:Ny] = -v[:,Ny-1:Ny] 54 r_even[N] = (v[N-1] + v[N])/2
55 r_odd[N] = -(v[N-1] + v[N])/2
56
88 57
89 @testset "Application" begin 58 @testset "Application" begin
90 @test op_x*v ≈ rx rtol = 1e-14 59 @test op_even*v ≈ r_even
91 @test op_y*v ≈ ry rtol = 1e-14 60 @test op_odd*v ≈ r_odd
61
62 @test (op_even*rand(ComplexF64,size(g)))[2] isa ComplexF64
92 end 63 end
93 64
94 @testset "Regions" begin 65 @testset "Regions" begin
95 @test (op_x*v)[Index(1,Lower),Index(3,Interior)] ≈ rx[1,3] rtol = 1e-14 66 @test (op_even*v)[Index(1,Lower)] ≈ r_even[1]
96 @test (op_x*v)[Index(2,Lower),Index(3,Interior)] ≈ rx[2,3] rtol = 1e-14 67 @test (op_even*v)[Index(2,Lower)] ≈ r_even[2]
97 @test (op_x*v)[Index(6,Interior),Index(3,Interior)] ≈ rx[6,3] rtol = 1e-14 68 @test (op_even*v)[Index(6,Interior)] ≈ r_even[6]
98 @test (op_x*v)[Index(10,Upper),Index(3,Interior)] ≈ rx[10,3] rtol = 1e-14 69 @test (op_even*v)[Index(10,Upper)] ≈ r_even[10]
99 @test (op_x*v)[Index(11,Upper),Index(3,Interior)] ≈ rx[11,3] rtol = 1e-14 70 @test (op_even*v)[Index(11,Upper)] ≈ r_even[11]
100 71
101 @test_throws BoundsError (op_x*v)[Index(3,Lower),Index(3,Interior)] 72 @test_throws BoundsError (op_even*v)[Index(3,Lower)]
102 @test_throws BoundsError (op_x*v)[Index(9,Upper),Index(3,Interior)] 73 @test_throws BoundsError (op_even*v)[Index(9,Upper)]
103
104 @test (op_y*v)[Index(3,Interior),Index(1,Lower)] ≈ ry[3,1] rtol = 1e-14
105 @test (op_y*v)[Index(3,Interior),Index(2,Lower)] ≈ ry[3,2] rtol = 1e-14
106 @test (op_y*v)[Index(3,Interior),Index(6,Interior)] ≈ ry[3,6] rtol = 1e-14
107 @test (op_y*v)[Index(3,Interior),Index(11,Upper)] ≈ ry[3,11] rtol = 1e-14
108 @test (op_y*v)[Index(3,Interior),Index(12,Upper)] ≈ ry[3,12] rtol = 1e-14
109
110 @test_throws BoundsError (op_y*v)[Index(3,Interior),Index(10,Upper)]
111 @test_throws BoundsError (op_y*v)[Index(3,Interior),Index(3,Lower)]
112 end 74 end
113 75
114 @testset "Inferred" begin 76 @testset "Inferred" begin
115 @test_skip @inferred apply(op_x, v,1,1) 77 @inferred apply(op_even, v, 1)
116 @inferred apply(op_x, v, Index(1,Lower),Index(1,Lower)) 78 @inferred apply(op_even, v, Index(1,Lower))
117 @inferred apply(op_x, v, Index(6,Interior),Index(1,Lower)) 79 @inferred apply(op_even, v, Index(6,Interior))
118 @inferred apply(op_x, v, Index(11,Upper),Index(1,Lower)) 80 @inferred apply(op_even, v, Index(11,Upper))
119 @test_skip @inferred apply(op_y, v,1,1)
120 @inferred apply(op_y, v, Index(1,Lower),Index(1,Lower))
121 @inferred apply(op_y, v, Index(1,Lower),Index(6,Interior))
122 @inferred apply(op_y, v, Index(1,Lower),Index(11,Upper))
123 end 81 end
124 end 82 end