Mercurial > repos > public > sbplib_julia

comparison test/testSbpOperators.jl @ 634:fb5ac62563aa feature/volume_and_boundary_operators

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Integrate feature/quadrature_as_outer_product into branch, before closing feature/quadrature_as_outer_product. (It is now obsolete apart from tests)
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 01 Jan 2021 16:39:57 +0100
parents bf8b66c596f7 a78bda7084f6
children 08b2c7a2d063
comparison
equal deleted inserted replaced
632:bf8b66c596f7 634:fb5ac62563aa
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 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
404 end 404 end
405 end 405 end
406 end 406 end
407 407
408 @testset "DiagonalInnerProduct" begin 408 @testset "DiagonalQuadrature" begin
409 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 409 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
410 L = 2.3 410 Lx = π/2.
411 g = EquidistantGrid(77, 0.0, L) 411 Ly = Float64(π)
412 H = DiagonalInnerProduct(g,op.quadratureClosure) 412 g_1D = EquidistantGrid(77, 0.0, Lx)
413 v = ones(Float64, size(g)) 413 g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
414 414 integral(H,v) = sum(H*v)
415 @test H isa TensorMapping{T,1,1} where T 415 @testset "Constructors" begin
416 @test H' isa TensorMapping{T,1,1} where T 416 # 1D
417 @test sum(H*v) ≈ L 417 H_x = DiagonalQuadrature(spacing(g_1D)[1],op.quadratureClosure,size(g_1D));
418 @test H*v == H'*v 418 @test H_x == DiagonalQuadrature(g_1D,op.quadratureClosure)
419 end 419 @test H_x == diagonal_quadrature(g_1D,op.quadratureClosure)
420 420 @test H_x isa TensorMapping{T,1,1} where T
421 @testset "Quadrature" begin 421 @test H_x' isa TensorMapping{T,1,1} where T
422 # 2D
423 H_xy = diagonal_quadrature(g_2D,op.quadratureClosure)
424 @test H_xy isa TensorMappingComposition
425 @test H_xy isa TensorMapping{T,2,2} where T
426 @test H_xy' isa TensorMapping{T,2,2} where T
427 end
428
429 @testset "Sizes" begin
430 # 1D
431 H_x = diagonal_quadrature(g_1D,op.quadratureClosure)
432 @test domain_size(H_x) == size(g_1D)
433 @test range_size(H_x) == size(g_1D)
434 # 2D
435 H_xy = diagonal_quadrature(g_2D,op.quadratureClosure)
436 @test domain_size(H_xy) == size(g_2D)
437 @test range_size(H_xy) == size(g_2D)
438 end
439
440 @testset "Application" begin
441 # 1D
442 H_x = diagonal_quadrature(g_1D,op.quadratureClosure)
443 a = 3.2
444 v_1D = a*ones(Float64, size(g_1D))
445 u_1D = evalOn(g_1D,x->sin(x))
446 @test integral(H_x,v_1D) ≈ a*Lx rtol = 1e-13
447 @test integral(H_x,u_1D) ≈ 1. rtol = 1e-8
448 @test H_x*v_1D == H_x'*v_1D
449 # 2D
450 H_xy = diagonal_quadrature(g_2D,op.quadratureClosure)
451 b = 2.1
452 v_2D = b*ones(Float64, size(g_2D))
453 u_2D = evalOn(g_2D,(x,y)->sin(x)+cos(y))
454 @test integral(H_xy,v_2D) ≈ b*Lx*Ly rtol = 1e-13
455 @test integral(H_xy,u_2D) ≈ π rtol = 1e-8
456 @test H_xy*v_2D ≈ H_xy'*v_2D rtol = 1e-16 #Failed for exact equality. Must differ in operation order for some reason?
457 end
458
459 @testset "Accuracy" begin
460 v = ()
461 for i = 0:4
462 f_i(x) = 1/factorial(i)*x^i
463 v = (v...,evalOn(g_1D,f_i))
464 end
465 # TODO: Bug in readOperator for 2nd order
466 # # 2nd order
467 # op2 = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
468 # H2 = diagonal_quadrature(g_1D,op2.quadratureClosure)
469 # for i = 1:3
470 # @test integral(H2,v[i]) ≈ v[i+1] rtol = 1e-14
471 # end
472
473 # 4th order
474 op4 = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
475 H4 = diagonal_quadrature(g_1D,op4.quadratureClosure)
476 for i = 1:4
477 @test integral(H4,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14
478 end
479 end
480
481 @testset "Inferred" begin
482 H_x = diagonal_quadrature(g_1D,op.quadratureClosure)
483 H_xy = diagonal_quadrature(g_2D,op.quadratureClosure)
484 v_1D = ones(Float64, size(g_1D))
485 v_2D = ones(Float64, size(g_2D))
486 @inferred H_x*v_1D
487 @inferred H_x'*v_1D
488 @inferred H_xy*v_2D
489 @inferred H_xy'*v_2D
490 end
491 end
492
493 @testset "InverseDiagonalQuadrature" begin
422 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 494 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
423 Lx = 2.3 495 Lx = π/2.
424 Ly = 5.2 496 Ly = Float64(π)
425 g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) 497 g_1D = EquidistantGrid(77, 0.0, Lx)
426 498 g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
427 Q = Quadrature(g, op.quadratureClosure) 499 @testset "Constructors" begin
428 500 # 1D
429 @test Q isa TensorMapping{T,2,2} where T 501 Hi_x = InverseDiagonalQuadrature(inverse_spacing(g_1D)[1], 1. ./ op.quadratureClosure, size(g_1D));
430 @test Q' isa TensorMapping{T,2,2} where T 502 @test Hi_x == InverseDiagonalQuadrature(g_1D,op.quadratureClosure)
431 503 @test Hi_x == inverse_diagonal_quadrature(g_1D,op.quadratureClosure)
432 v = ones(Float64, size(g)) 504 @test Hi_x isa TensorMapping{T,1,1} where T
433 @test sum(Q*v) ≈ Lx*Ly 505 @test Hi_x' isa TensorMapping{T,1,1} where T
434 506
435 v = 2*ones(Float64, size(g)) 507 # 2D
436 @test_broken sum(Q*v) ≈ 2*Lx*Ly 508 Hi_xy = inverse_diagonal_quadrature(g_2D,op.quadratureClosure)
437 509 @test Hi_xy isa TensorMappingComposition
438 @test Q*v == Q'*v 510 @test Hi_xy isa TensorMapping{T,2,2} where T
439 end 511 @test Hi_xy' isa TensorMapping{T,2,2} where T
440 512 end
441 @testset "InverseDiagonalInnerProduct" begin 513
442 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 514 @testset "Sizes" begin
443 L = 2.3 515 # 1D
444 g = EquidistantGrid(77, 0.0, L) 516 Hi_x = inverse_diagonal_quadrature(g_1D,op.quadratureClosure)
445 H = DiagonalInnerProduct(g, op.quadratureClosure) 517 @test domain_size(Hi_x) == size(g_1D)
446 Hi = InverseDiagonalInnerProduct(g,op.quadratureClosure) 518 @test range_size(Hi_x) == size(g_1D)
447 v = evalOn(g, x->sin(x)) 519 # 2D
448 520 Hi_xy = inverse_diagonal_quadrature(g_2D,op.quadratureClosure)
449 @test Hi isa TensorMapping{T,1,1} where T 521 @test domain_size(Hi_xy) == size(g_2D)
450 @test Hi' isa TensorMapping{T,1,1} where T 522 @test range_size(Hi_xy) == size(g_2D)
451 @test Hi*H*v ≈ v 523 end
452 @test Hi*v == Hi'*v 524
453 end 525 @testset "Application" begin
454 526 # 1D
455 @testset "InverseQuadrature" begin 527 H_x = diagonal_quadrature(g_1D,op.quadratureClosure)
456 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 528 Hi_x = inverse_diagonal_quadrature(g_1D,op.quadratureClosure)
457 Lx = 7.3 529 v_1D = evalOn(g_1D,x->sin(x))
458 Ly = 8.2 530 u_1D = evalOn(g_1D,x->x^3-x^2+1)
459 g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) 531 @test Hi_x*H_x*v_1D ≈ v_1D rtol = 1e-15
460 532 @test Hi_x*H_x*u_1D ≈ u_1D rtol = 1e-15
461 Q = Quadrature(g, op.quadratureClosure) 533 @test Hi_x*v_1D == Hi_x'*v_1D
462 Qinv = InverseQuadrature(g, op.quadratureClosure) 534 # 2D
463 v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y) 535 H_xy = diagonal_quadrature(g_2D,op.quadratureClosure)
464 536 Hi_xy = inverse_diagonal_quadrature(g_2D,op.quadratureClosure)
465 @test Qinv isa TensorMapping{T,2,2} where T 537 v_2D = evalOn(g_2D,(x,y)->sin(x)+cos(y))
466 @test Qinv' isa TensorMapping{T,2,2} where T 538 u_2D = evalOn(g_2D,(x,y)->x*y + x^5 - sqrt(y))
467 @test_broken Qinv*(Q*v) ≈ v 539 @test Hi_xy*H_xy*v_2D ≈ v_2D rtol = 1e-15
468 @test Qinv*v == Qinv'*v 540 @test Hi_xy*H_xy*u_2D ≈ u_2D rtol = 1e-15
541 @test Hi_xy*v_2D ≈ Hi_xy'*v_2D rtol = 1e-16 #Failed for exact equality. Must differ in operation order for some reason?
542 end
543
544 @testset "Inferred" begin
545 Hi_x = inverse_diagonal_quadrature(g_1D,op.quadratureClosure)
546 Hi_xy = inverse_diagonal_quadrature(g_2D,op.quadratureClosure)
547 v_1D = ones(Float64, size(g_1D))
548 v_2D = ones(Float64, size(g_2D))
549 @inferred Hi_x*v_1D
550 @inferred Hi_x'*v_1D
551 @inferred Hi_xy*v_2D
552 @inferred Hi_xy'*v_2D
553 end
469 end 554 end
470 555
471 @testset "BoundaryOperator" begin 556 @testset "BoundaryOperator" begin
472 closure_stencil = Stencil((0,2), (2.,1.,3.)) 557 closure_stencil = Stencil((0,2), (2.,1.,3.))
473 g_1D = EquidistantGrid(11, 0.0, 1.0) 558 g_1D = EquidistantGrid(11, 0.0, 1.0)