Mercurial > repos > public > sbplib_julia

comparison test/testSbpOperators.jl @ 553:d5c45785d318 feature/quadrature_as_outer_product

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add test of accuracy conditions for DiagonalQuadrature
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Sun, 29 Nov 2020 18:57:29 +0100
parents 9bfecd6b6aac
children dab9df9c4d66
comparison
equal deleted inserted replaced
552:9bfecd6b6aac 553:d5c45785d318
114 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 114 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
115 Lx = π/2. 115 Lx = π/2.
116 Ly = Float64(π) 116 Ly = Float64(π)
117 g_1D = EquidistantGrid(77, 0.0, Lx) 117 g_1D = EquidistantGrid(77, 0.0, Lx)
118 g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) 118 g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
119 integral(H,v) = sum(H*v)
119 @testset "Constructors" begin 120 @testset "Constructors" begin
120 # 1D 121 # 1D
121 H_x = DiagonalQuadrature(spacing(g_1D)[1],op.quadratureClosure,size(g_1D)); 122 H_x = DiagonalQuadrature(spacing(g_1D)[1],op.quadratureClosure,size(g_1D));
122 @test H_x == DiagonalQuadrature(g_1D,op.quadratureClosure) 123 @test H_x == DiagonalQuadrature(g_1D,op.quadratureClosure)
123 @test H_x == diagonal_quadrature(g_1D,op.quadratureClosure) 124 @test H_x == diagonal_quadrature(g_1D,op.quadratureClosure)
129 @test H_xy isa TensorMappingComposition 130 @test H_xy isa TensorMappingComposition
130 @test H_xy isa TensorMapping{T,2,2} where T 131 @test H_xy isa TensorMapping{T,2,2} where T
131 @test H_xy' isa TensorMapping{T,2,2} where T 132 @test H_xy' isa TensorMapping{T,2,2} where T
132 end 133 end
133 134
134 H_x = diagonal_quadrature(g_1D,op.quadratureClosure)
135 H_xy = diagonal_quadrature(g_2D,op.quadratureClosure)
136 @testset "Sizes" begin 135 @testset "Sizes" begin
136 H_x = diagonal_quadrature(g_1D,op.quadratureClosure)
137 H_xy = diagonal_quadrature(g_2D,op.quadratureClosure)
137 # 1D 138 # 1D
138 @test domain_size(H_x) == size(g_1D) 139 @test domain_size(H_x) == size(g_1D)
139 @test range_size(H_x) == size(g_1D) 140 @test range_size(H_x) == size(g_1D)
140 # 2D 141 # 2D
141 @test domain_size(H_xy) == size(g_2D) 142 @test domain_size(H_xy) == size(g_2D)
142 @test range_size(H_xy) == size(g_2D) 143 @test range_size(H_xy) == size(g_2D)
143 end 144 end
144 a = 3.2 145
145 b = 2.1
146 v_1D = a*ones(Float64, size(g_1D))
147 v_2D = b*ones(Float64, size(g_2D))
148 u_1D = evalOn(g_1D,x->sin(x))
149 u_2D = evalOn(g_2D,(x,y)->sin(x)+cos(y))
150 integral(H,v) = sum(H*v)
151 @testset "Application" begin 146 @testset "Application" begin
147 H_x = diagonal_quadrature(g_1D,op.quadratureClosure)
148 H_xy = diagonal_quadrature(g_2D,op.quadratureClosure)
149 a = 3.2
150 b = 2.1
151 v_1D = a*ones(Float64, size(g_1D))
152 v_2D = b*ones(Float64, size(g_2D))
153 u_1D = evalOn(g_1D,x->sin(x))
154 u_2D = evalOn(g_2D,(x,y)->sin(x)+cos(y))
152 # 1D 155 # 1D
153 @test integral(H_x,v_1D) ≈ a*Lx rtol = 1e-13 156 @test integral(H_x,v_1D) ≈ a*Lx rtol = 1e-13
154 @test integral(H_x,u_1D) ≈ 1. rtol = 1e-8 157 @test integral(H_x,u_1D) ≈ 1. rtol = 1e-8
155 @test H_x*v_1D == H_x'*v_1D 158 @test H_x*v_1D == H_x'*v_1D
156 # 2D 159 # 2D
157 @test integral(H_xy,v_2D) ≈ b*Lx*Ly rtol = 1e-13 160 @test integral(H_xy,v_2D) ≈ b*Lx*Ly rtol = 1e-13
158 @test integral(H_xy,u_2D) ≈ π rtol = 1e-8 161 @test integral(H_xy,u_2D) ≈ π rtol = 1e-8
159 @test H_xy*v_2D ≈ H_xy'*v_2D rtol = 1e-16 #Failed for exact equality. Must differ in operation order for some reason? 162 @test H_xy*v_2D ≈ H_xy'*v_2D rtol = 1e-16 #Failed for exact equality. Must differ in operation order for some reason?
160 end 163 end
161 164
165 @testset "Accuracy" begin
166 v = ()
167 for i = 0:4
168 f_i(x) = 1/factorial(i)*x^i
169 v = (v...,evalOn(g_1D,f_i))
170 end
171 # # 2nd order
172 # op2 = readOperator(sbp_operators_path()*"d2_2nd.txt",sbp_operators_path()*"h_2nd.txt")
173 # H2 = diagonal_quadrature(g_1D,op2.quadratureClosure)
174 # for i = 1:3
175 # @test integral(H4,v[i]) ≈ v[i+1] rtol = 1e-14
176 # end
177
178 # 4th order
179 op4 = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
180 H4 = diagonal_quadrature(g_1D,op4.quadratureClosure)
181 for i = 1:4
182 @test integral(H4,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14
183 end
184 end
185
162 @testset "Inferred" begin 186 @testset "Inferred" begin
187 H_x = diagonal_quadrature(g_1D,op.quadratureClosure)
188 H_xy = diagonal_quadrature(g_2D,op.quadratureClosure)
189 v_1D = ones(Float64, size(g_1D))
190 v_2D = ones(Float64, size(g_2D))
163 @inferred H_x*v_1D 191 @inferred H_x*v_1D
164 @inferred H_x'*v_1D 192 @inferred H_x'*v_1D
165 @inferred H_xy*v_2D 193 @inferred H_xy*v_2D
166 @inferred H_xy'*v_2D 194 @inferred H_xy'*v_2D
167 end 195 end