Mercurial > repos > public > sbplib_julia

comparison test/testSbpOperators.jl @ 643:0928bbc3ee8b feature/volume_and_boundary_operators

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add tests for SecondDerivative and Laplace for 2nd order accurate case
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 04 Jan 2021 17:42:42 +0100
parents f4a16b403487
children e3fd4b7aa789
comparison
equal deleted inserted replaced
642:f4a16b403487 643:0928bbc3ee8b
274 @test Dₓₓ == D2⊗I 274 @test Dₓₓ == D2⊗I
275 @test Dₓₓ isa TensorMapping{T,2,2} where T 275 @test Dₓₓ isa TensorMapping{T,2,2} where T
276 end 276 end
277 end 277 end
278 278
279 # Exact differentiation is measured point-wise. In other cases
280 # the error is measured in the l2-norm.
279 @testset "Accuracy" begin 281 @testset "Accuracy" begin
280 @testset "1D" begin 282 @testset "1D" begin
283 l2(v) = sqrt(spacing(g_1D)[1]*sum(v.^2));
281 monomials = () 284 monomials = ()
282 maxOrder = 4; 285 maxOrder = 4;
283 for i = 0:maxOrder-1 286 for i = 0:maxOrder-1
284 f_i(x) = 1/factorial(i)*x^i 287 f_i(x) = 1/factorial(i)*x^i
285 monomials = (monomials...,evalOn(g_1D,f_i)) 288 monomials = (monomials...,evalOn(g_1D,f_i))
286 end 289 end
287 l2(v) = sqrt(spacing(g_1D)[1]*sum(v.^2)); 290 v = evalOn(g_1D,x -> sin(x))
288 291 vₓₓ = evalOn(g_1D,x -> -sin(x))
289 #TODO: Error when reading second order stencil! 292
290 # @testset "2nd order" begin 293 # 2nd order interior stencil, 1nd order boundary stencil,
291 # op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) 294 # implies that L*v should be exact for monomials up to order 2.
292 # Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils) 295 @testset "2nd order" begin
293 # @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-13 296 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
294 # @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-13 297 Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils)
295 # end 298 @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
299 @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
300 @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10
301 @test Dₓₓ*v ≈ vₓₓ rtol = 5e-2 norm = l2
302 end
296 303
297 # 4th order interior stencil, 2nd order boundary stencil, 304 # 4th order interior stencil, 2nd order boundary stencil,
298 # implies that L*v should be exact for monomials up to order 3. 305 # implies that L*v should be exact for monomials up to order 3.
299 # Exact differentiation is measured point-wise. For other grid functions
300 # the error is measured in the l2-norm.
301 @testset "4th order" begin 306 @testset "4th order" begin
302 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 307 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
303 Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils) 308 Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils)
304 # TODO: high tolerances for checking the "exact" differentiation 309 # TODO: high tolerances for checking the "exact" differentiation
305 # due to accumulation of round-off errors/cancellation errors? 310 # due to accumulation of round-off errors/cancellation errors?
306 @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10 311 @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
307 @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10 312 @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
308 @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10 313 @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10
309 @test Dₓₓ*monomials[4] ≈ monomials[2] atol = 5e-10 314 @test Dₓₓ*monomials[4] ≈ monomials[2] atol = 5e-10
310 @test Dₓₓ*evalOn(g_1D,x -> sin(x)) ≈ evalOn(g_1D,x -> -sin(x)) rtol = 5e-4 norm = l2 315 @test Dₓₓ*v ≈ vₓₓ rtol = 5e-4 norm = l2
311 end 316 end
312 end 317 end
313 318
314 @testset "2D" begin 319 @testset "2D" begin
320 l2(v) = sqrt(prod(spacing(g_2D))*sum(v.^2));
315 binomials = () 321 binomials = ()
316 maxOrder = 4; 322 maxOrder = 4;
317 for i = 0:maxOrder-1 323 for i = 0:maxOrder-1
318 f_i(x,y) = 1/factorial(i)*y^i + x^i 324 f_i(x,y) = 1/factorial(i)*y^i + x^i
319 binomials = (binomials...,evalOn(g_2D,f_i)) 325 binomials = (binomials...,evalOn(g_2D,f_i))
320 end 326 end
321 l2(v) = sqrt(prod(spacing(g_2D))*sum(v.^2)); 327 v = evalOn(g_2D, (x,y) -> sin(x)+cos(y))
322 328 v_yy = evalOn(g_2D,(x,y) -> -cos(y))
323 #TODO: Error when reading second order stencil! 329
324 # @testset "2nd order" begin 330 # 2nd order interior stencil, 1st order boundary stencil,
325 # op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) 331 # implies that L*v should be exact for binomials up to order 2.
326 # Dyy = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,2) 332 @testset "2nd order" begin
327 # @test Dyy*binomials[1] ≈ evalOn(g_2D,(x,y)->0.) atol = 5e-12 333 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
328 # @test Dyy*binomials[2] ≈ evalOn(g_2D,(x,y)->0.) atol = 5e-12 334 Dyy = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,2)
329 # end 335 @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
336 @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
337 @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9
338 @test Dyy*v ≈ v_yy rtol = 5e-2 norm = l2
339 end
330 340
331 # 4th order interior stencil, 2nd order boundary stencil, 341 # 4th order interior stencil, 2nd order boundary stencil,
332 # implies that L*v should be exact for binomials up to order 3. 342 # implies that L*v should be exact for binomials up to order 3.
333 # Exact differentiation is measured point-wise. For other grid functions
334 # the error is measured in the l2-norm.
335 @testset "4th order" begin 343 @testset "4th order" begin
336 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 344 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
337 Dyy = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,2) 345 Dyy = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,2)
338 # TODO: high tolerances for checking the "exact" differentiation 346 # TODO: high tolerances for checking the "exact" differentiation
339 # due to accumulation of round-off errors/cancellation errors? 347 # due to accumulation of round-off errors/cancellation errors?
340 @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9 348 @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
341 @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9 349 @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
342 @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9 350 @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9
343 @test Dyy*binomials[4] ≈ evalOn(g_2D,(x,y)->y) atol = 5e-9 351 @test Dyy*binomials[4] ≈ evalOn(g_2D,(x,y)->y) atol = 5e-9
344 @test Dyy*evalOn(g_2D, (x,y) -> sin(x)+cos(y)) ≈ evalOn(g_2D,(x,y) -> -cos(y)) rtol = 5e-4 norm = l2 352 @test Dyy*v ≈ v_yy rtol = 5e-4 norm = l2
345 end 353 end
346 end 354 end
347 end 355 end
348 end 356 end
349 357
350 @testset "Laplace" begin 358 @testset "Laplace" begin
351 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
352 g_1D = EquidistantGrid(101, 0.0, 1.) 359 g_1D = EquidistantGrid(101, 0.0, 1.)
353 #TODO: It's nice to verify that 3D works somewhere at least, but perhaps should keep 3D tests to a minimum to avoid 360 #TODO: It's nice to verify that 3D works somewhere at least, but perhaps should keep 3D tests to a minimum to avoid
354 # long run time for test? 361 # long run time for test?
355 g_3D = EquidistantGrid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.)) 362 g_3D = EquidistantGrid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.))
356 # TODO: These areant really constructors. Better name? 363 # TODO: These areant really constructors. Better name?
357 @testset "Constructors" begin 364 @testset "Constructors" begin
365 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
358 @testset "1D" begin 366 @testset "1D" begin
359 L = Laplace(g_1D, op.innerStencil, op.closureStencils) 367 L = Laplace(g_1D, op.innerStencil, op.closureStencils)
360 @test L == SecondDerivative(g_1D, op.innerStencil, op.closureStencils) 368 @test L == SecondDerivative(g_1D, op.innerStencil, op.closureStencils)
361 @test L isa TensorMapping{T,1,1} where T 369 @test L isa TensorMapping{T,1,1} where T
362 end 370 end
368 Dzz = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,3) 376 Dzz = SecondDerivative(g_3D, op.innerStencil, op.closureStencils,3)
369 @test L == Dxx + Dyy + Dzz 377 @test L == Dxx + Dyy + Dzz
370 end 378 end
371 end 379 end
372 380
381 # Exact differentiation is measured point-wise. In other cases
382 # the error is measured in the l2-norm.
373 @testset "Accuracy" begin 383 @testset "Accuracy" begin
384 l2(v) = sqrt(prod(spacing(g_3D))*sum(v.^2));
374 polynomials = () 385 polynomials = ()
375 maxOrder = 4; 386 maxOrder = 4;
376 for i = 0:maxOrder-1 387 for i = 0:maxOrder-1
377 f_i(x,y,z) = 1/factorial(i)*(y^i + x^i + z^i) 388 f_i(x,y,z) = 1/factorial(i)*(y^i + x^i + z^i)
378 polynomials = (polynomials...,evalOn(g_3D,f_i)) 389 polynomials = (polynomials...,evalOn(g_3D,f_i))
379 end 390 end
380 l2(v) = sqrt(prod(spacing(g_3D))*sum(v.^2)); 391 v = evalOn(g_3D, (x,y,z) -> sin(x) + cos(y) + exp(z))
381 392 Δv = evalOn(g_3D,(x,y,z) -> -sin(x) - cos(y) + exp(z))
382 #TODO: Error when reading second order stencil! 393
383 # @testset "2nd order" begin 394 # 2nd order interior stencil, 1st order boundary stencil,
384 # op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) 395 # implies that L*v should be exact for binomials up to order 2.
385 # Dyy = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,2) 396 @testset "2nd order" begin
386 # @test Dyy*binomials[1] ≈ evalOn(g_2D,(x,y)->0.) atol = 5e-12 397 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
387 # @test Dyy*binomials[2] ≈ evalOn(g_2D,(x,y)->0.) atol = 5e-12 398 L = Laplace(g_3D,op.innerStencil,op.closureStencils)
388 # end 399 @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
400 @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
401 @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9
402 @test L*v ≈ Δv rtol = 5e-2 norm = l2
403 end
389 404
390 # 4th order interior stencil, 2nd order boundary stencil, 405 # 4th order interior stencil, 2nd order boundary stencil,
391 # implies that L*v should be exact for binomials up to order 3. 406 # implies that L*v should be exact for binomials up to order 3.
392 # Exact differentiation is measured point-wise. For other grid functions
393 # the error is measured in the l2-norm.
394 @testset "4th order" begin 407 @testset "4th order" begin
395 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 408 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
396 L = Laplace(g_3D,op.innerStencil,op.closureStencils) 409 L = Laplace(g_3D,op.innerStencil,op.closureStencils)
397 # TODO: high tolerances for checking the "exact" differentiation 410 # TODO: high tolerances for checking the "exact" differentiation
398 # due to accumulation of round-off errors/cancellation errors? 411 # due to accumulation of round-off errors/cancellation errors?
399 @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9 412 @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
400 @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9 413 @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
401 @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9 414 @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9
402 @test L*polynomials[4] ≈ polynomials[2] atol = 5e-9 415 @test L*polynomials[4] ≈ polynomials[2] atol = 5e-9
403 @test L*evalOn(g_3D, (x,y,z) -> sin(x) + cos(y) + exp(z)) ≈ evalOn(g_3D,(x,y,z) -> -sin(x)-cos(y) + exp(z)) rtol = 5e-4 norm = l2 416 @test L*v ≈ Δv rtol = 5e-4 norm = l2
404 end 417 end
405 end 418 end
406 end 419 end
407 420
408 @testset "DiagonalQuadrature" begin 421 @testset "DiagonalQuadrature" begin
409 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
410 Lx = π/2. 422 Lx = π/2.
411 Ly = Float64(π) 423 Ly = Float64(π)
412 g_1D = EquidistantGrid(77, 0.0, Lx) 424 g_1D = EquidistantGrid(77, 0.0, Lx)
413 g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) 425 g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
414 integral(H,v) = sum(H*v) 426 integral(H,v) = sum(H*v)
415 @testset "Constructors" begin 427 @testset "Constructors" begin
428 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
416 @testset "1D" begin 429 @testset "1D" begin
417 H = DiagonalQuadrature(g_1D,op.quadratureClosure) 430 H = DiagonalQuadrature(g_1D,op.quadratureClosure)
418 inner_stencil = Stencil((1.,),center=1) 431 inner_stencil = Stencil((1.,),center=1)
419 @test H == Quadrature(g_1D,inner_stencil,op.quadratureClosure) 432 @test H == Quadrature(g_1D,inner_stencil,op.quadratureClosure)
420 @test H isa TensorMapping{T,1,1} where T 433 @test H isa TensorMapping{T,1,1} where T
427 @test H isa TensorMapping{T,2,2} where T 440 @test H isa TensorMapping{T,2,2} where T
428 end 441 end
429 end 442 end
430 443
431 @testset "Sizes" begin 444 @testset "Sizes" begin
445 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
432 @testset "1D" begin 446 @testset "1D" begin
433 H = DiagonalQuadrature(g_1D,op.quadratureClosure) 447 H = DiagonalQuadrature(g_1D,op.quadratureClosure)
434 @test domain_size(H) == size(g_1D) 448 @test domain_size(H) == size(g_1D)
435 @test range_size(H) == size(g_1D) 449 @test range_size(H) == size(g_1D)
436 end 450 end