Mercurial > repos > public > sbplib_julia

comparison test/testSbpOperators.jl @ 621:60d7c1752cfa feature/volume_and_boundary_operators

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Restructure tests for SecondDerivative
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 16 Dec 2020 16:42:29 +0100
parents bfc893d03cbf
children f799678357df
comparison
equal deleted inserted replaced
620:bfc893d03cbf 621:60d7c1752cfa
5 using Sbplib.LazyTensors 5 using Sbplib.LazyTensors
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.BoundaryOperator 11 import Sbplib.SbpOperators.BoundaryOperator
11 import Sbplib.SbpOperators.boundary_operator 12 import Sbplib.SbpOperators.boundary_operator
12 13
13 @testset "SbpOperators" begin 14 @testset "SbpOperators" begin
14 15
131 # end 132 # end
132 # end 133 # end
133 134
134 @testset "SecondDerivative" begin 135 @testset "SecondDerivative" begin
135 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 136 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
136 L = 3.5 137 Lx = 3.5
137 g = EquidistantGrid(101, 0.0, L) 138 Ly = 3.
138 Dₓₓ = SecondDerivative(g,op.innerStencil,op.closureStencils) 139 g_1D = EquidistantGrid(121, 0.0, Lx)
139 140 g_2D = EquidistantGrid((121,123), (0.0, 0.0), (Lx, Ly))
140 f0(x) = 1. 141
141 f1(x) = x 142 @testset "Constructors" begin
142 f2(x) = 1/2*x^2 143 @testset "1D" begin
143 f3(x) = 1/6*x^3 144 Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils)
144 f4(x) = 1/24*x^4 145 Dₓₓ == SecondDerivative(g_1D,op.innerStencil,op.closureStencils,1)
145 f5(x) = sin(x) 146 Dₓₓ isa VolumeOperator
146 f5ₓₓ(x) = -f5(x) 147 end
147 148 @testset "2D" begin
148 v0 = evalOn(g,f0) 149 Dₓₓ = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,1)
149 v1 = evalOn(g,f1) 150 Dₓₓ isa TensorMappingComposition
150 v2 = evalOn(g,f2) 151 Dₓₓ isa TensorMapping{2,2}
151 v3 = evalOn(g,f3) 152 end
152 v4 = evalOn(g,f4) 153 end
153 v5 = evalOn(g,f5) 154
154 155 @testset "Accuracy" begin
155 @test Dₓₓ isa TensorMapping{T,1,1} where T 156
156 @test Dₓₓ' isa TensorMapping{T,1,1} where T 157 @testset "1D" begin
157 158 v0 = evalOn(g_1D,x->0.)
158 # 4th order interior stencil, 2nd order boundary stencil, 159 monomials = ()
159 # implies that L*v should be exact for v - monomial up to order 3. 160 maxOrder = 4;
160 # Exact differentiation is measured point-wise. For other grid functions 161 for i = 0:maxOrder-1
161 # the error is measured in the l2-norm. 162 f_i(x) = 1/factorial(i)*x^i
162 @test norm(Dₓₓ*v0) ≈ 0.0 atol=5e-10 163 monomials = (monomials...,evalOn(g_1D,f_i))
163 @test norm(Dₓₓ*v1) ≈ 0.0 atol=5e-10 164 end
164 @test Dₓₓ*v2 ≈ v0 atol=5e-11 165 l2(v) = sqrt(spacing(g_1D)[1]*sum(v.^2));
165 @test Dₓₓ*v3 ≈ v1 atol=5e-11 166
166 167 #TODO: Error when reading second order stencil!
167 h = spacing(g)[1]; 168 # @testset "2nd order" begin
168 l2(v) = sqrt(h*sum(v.^2)) 169 # op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
169 @test Dₓₓ*v4 ≈ v2 atol=5e-4 norm=l2 170 # Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils)
170 @test Dₓₓ*v5 ≈ -v5 atol=5e-4 norm=l2 171 # @test Dₓₓ*monomials[1] ≈ v0 atol = 5e-13
172 # @test Dₓₓ*monomials[2] ≈ v0 atol = 5e-13
173 # end
174
175 # 4th order interior stencil, 2nd order boundary stencil,
176 # implies that L*v should be exact for monomials up to order 3.
177 # Exact differentiation is measured point-wise. For other grid functions
178 # the error is measured in the l2-norm.
179 @testset "4th order" begin
180 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
181 Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils)
182 # TODO: high tolerances for checking the "exact" differentiation
183 # due to accumulation of round-off errors/cancellation errors?
184 @test Dₓₓ*monomials[1] ≈ v0 atol = 5e-10
185 @test Dₓₓ*monomials[2] ≈ v0 atol = 5e-10
186 @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10
187 @test Dₓₓ*monomials[4] ≈ monomials[2] atol = 5e-10
188 @test Dₓₓ*evalOn(g_1D,x -> sin(x)) ≈ evalOn(g_1D,x -> -sin(x)) rtol = 5e-4 norm = l2
189 end
190 end
191
192 #TODO: Error when reading second order stencil!
193 @testset "2D" begin
194 binomials = ()
195 maxOrder = 4;
196 for i = 0:maxOrder-1
197 f_i(x,y) = 1/factorial(i)*y^i + x^i
198 binomials = (binomials...,evalOn(g_2D,f_i))
199 end
200 l2(v) = sqrt(prod(spacing(g_2D))*sum(v.^2));
201 # @testset "2nd order" begin
202 # op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
203 # Dyy = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,2)
204 # @test Dyy*binomials[1] ≈ evalOn(g_2D,(x,y)->0.) atol = 5e-12
205 # @test Dyy*binomials[2] ≈ evalOn(g_2D,(x,y)->0.) atol = 5e-12
206 # end
207
208 # 4th order interior stencil, 2nd order boundary stencil,
209 # implies that L*v should be exact for binomials up to order 3.
210 # Exact differentiation is measured point-wise. For other grid functions
211 # the error is measured in the l2-norm.
212 @testset "4th order" begin
213 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
214 Dyy = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,2)
215 # TODO: high tolerances for checking the "exact" differentiation
216 # due to accumulation of round-off errors/cancellation errors?
217 @test Dyy*binomials[1] ≈ evalOn(g_2D,(x,y)->0.) atol = 5e-9
218 @test Dyy*binomials[2] ≈ evalOn(g_2D,(x,y)->0.) atol = 5e-9
219 @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9
220 @test Dyy*binomials[4] ≈ evalOn(g_2D,(x,y)->y) atol = 5e-9
221 @test Dyy*evalOn(g_2D, (x,y) -> sin(x)+cos(y)) ≈ evalOn(g_2D,(x,y) -> -cos(y)) rtol = 5e-4 norm = l2
222 end
223 end
224 end
171 end 225 end
172 226
173 227
174 @testset "Laplace2D" begin 228 @testset "Laplace2D" begin
175 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 229 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)