Mercurial > repos > public > sbplib_julia

comparison test/testSbpOperators.jl @ 642:f4a16b403487 feature/volume_and_boundary_operators

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Implement the inverse quadrature operator as a volume operator and update tests.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 04 Jan 2021 17:17:40 +0100
parents 5e50e9815732
children 0928bbc3ee8b
comparison
equal deleted inserted replaced
641:5e50e9815732 642:f4a16b403487
487 end 487 end
488 end 488 end
489 end 489 end
490 end 490 end
491 491
492 # @testset "InverseDiagonalQuadrature" begin 492 @testset "InverseDiagonalQuadrature" begin
493 # op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 493 Lx = π/2.
494 # Lx = π/2. 494 Ly = Float64(π)
495 # Ly = Float64(π) 495 g_1D = EquidistantGrid(77, 0.0, Lx)
496 # g_1D = EquidistantGrid(77, 0.0, Lx) 496 g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
497 # g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) 497 @testset "Constructors" begin
498 # @testset "Constructors" begin 498 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
499 # # 1D 499 @testset "1D" begin
500 # Hi_x = InverseDiagonalQuadrature(inverse_spacing(g_1D)[1], 1. ./ op.quadratureClosure, size(g_1D)); 500 Hi = InverseDiagonalQuadrature(g_1D, op.quadratureClosure);
501 # @test Hi_x == InverseDiagonalQuadrature(g_1D,op.quadratureClosure) 501 inner_stencil = Stencil((1.,),center=1)
502 # @test Hi_x == inverse_diagonal_quadrature(g_1D,op.quadratureClosure) 502 closures = ()
503 # @test Hi_x isa TensorMapping{T,1,1} where T 503 for i = 1:length(op.quadratureClosure)
504 # @test Hi_x' isa TensorMapping{T,1,1} where T 504 closures = (closures...,Stencil(op.quadratureClosure[i].range,1.0./op.quadratureClosure[i].weights))
505 # 505 end
506 # # 2D 506 @test Hi == InverseQuadrature(g_1D,inner_stencil,closures)
507 # Hi_xy = inverse_diagonal_quadrature(g_2D,op.quadratureClosure) 507 @test Hi isa TensorMapping{T,1,1} where T
508 # @test Hi_xy isa TensorMappingComposition 508 end
509 # @test Hi_xy isa TensorMapping{T,2,2} where T 509 @testset "2D" begin
510 # @test Hi_xy' isa TensorMapping{T,2,2} where T 510 Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure)
511 # end 511 Hi_x = InverseDiagonalQuadrature(restrict(g_2D,1),op.quadratureClosure)
512 # 512 Hi_y = InverseDiagonalQuadrature(restrict(g_2D,2),op.quadratureClosure)
513 # @testset "Sizes" begin 513 @test Hi == Hi_x⊗Hi_y
514 # # 1D 514 @test Hi isa TensorMapping{T,2,2} where T
515 # Hi_x = inverse_diagonal_quadrature(g_1D,op.quadratureClosure) 515 end
516 # @test domain_size(Hi_x) == size(g_1D) 516 end
517 # @test range_size(Hi_x) == size(g_1D) 517
518 # # 2D 518 @testset "Sizes" begin
519 # Hi_xy = inverse_diagonal_quadrature(g_2D,op.quadratureClosure) 519 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
520 # @test domain_size(Hi_xy) == size(g_2D) 520 @testset "1D" begin
521 # @test range_size(Hi_xy) == size(g_2D) 521 Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure)
522 # end 522 @test domain_size(Hi) == size(g_1D)
523 # 523 @test range_size(Hi) == size(g_1D)
524 # @testset "Application" begin 524 end
525 # # 1D 525 @testset "2D" begin
526 # H_x = diagonal_quadrature(g_1D,op.quadratureClosure) 526 Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure)
527 # Hi_x = inverse_diagonal_quadrature(g_1D,op.quadratureClosure) 527 @test domain_size(Hi) == size(g_2D)
528 # v_1D = evalOn(g_1D,x->sin(x)) 528 @test range_size(Hi) == size(g_2D)
529 # u_1D = evalOn(g_1D,x->x^3-x^2+1) 529 end
530 # @test Hi_x*H_x*v_1D ≈ v_1D rtol = 1e-15 530 end
531 # @test Hi_x*H_x*u_1D ≈ u_1D rtol = 1e-15 531
532 # @test Hi_x*v_1D == Hi_x'*v_1D 532 @testset "Accuracy" begin
533 # # 2D 533 @testset "1D" begin
534 # H_xy = diagonal_quadrature(g_2D,op.quadratureClosure) 534 v = evalOn(g_1D,x->sin(x))
535 # Hi_xy = inverse_diagonal_quadrature(g_2D,op.quadratureClosure) 535 u = evalOn(g_1D,x->x^3-x^2+1)
536 # v_2D = evalOn(g_2D,(x,y)->sin(x)+cos(y)) 536 @testset "2nd order" begin
537 # u_2D = evalOn(g_2D,(x,y)->x*y + x^5 - sqrt(y)) 537 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
538 # @test Hi_xy*H_xy*v_2D ≈ v_2D rtol = 1e-15 538 H = DiagonalQuadrature(g_1D,op.quadratureClosure)
539 # @test Hi_xy*H_xy*u_2D ≈ u_2D rtol = 1e-15 539 Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure)
540 # @test Hi_xy*v_2D ≈ Hi_xy'*v_2D rtol = 1e-16 #Failed for exact equality. Must differ in operation order for some reason? 540 @test Hi*H*v ≈ v rtol = 1e-15
541 # end 541 @test Hi*H*u ≈ u rtol = 1e-15
542 # 542 end
543 # @testset "Inferred" begin 543 @testset "4th order" begin
544 # Hi_x = inverse_diagonal_quadrature(g_1D,op.quadratureClosure) 544 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
545 # Hi_xy = inverse_diagonal_quadrature(g_2D,op.quadratureClosure) 545 H = DiagonalQuadrature(g_1D,op.quadratureClosure)
546 # v_1D = ones(Float64, size(g_1D)) 546 Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure)
547 # v_2D = ones(Float64, size(g_2D)) 547 @test Hi*H*v ≈ v rtol = 1e-15
548 # @inferred Hi_x*v_1D 548 @test Hi*H*u ≈ u rtol = 1e-15
549 # @inferred Hi_x'*v_1D 549 end
550 # @inferred Hi_xy*v_2D 550 end
551 # @inferred Hi_xy'*v_2D 551 @testset "2D" begin
552 # end 552 v = evalOn(g_2D,(x,y)->sin(x)+cos(y))
553 # end 553 u = evalOn(g_2D,(x,y)->x*y + x^5 - sqrt(y))
554 @testset "2nd order" begin
555 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
556 H = DiagonalQuadrature(g_2D,op.quadratureClosure)
557 Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure)
558 @test Hi*H*v ≈ v rtol = 1e-15
559 @test Hi*H*u ≈ u rtol = 1e-15
560 end
561 @testset "4th order" begin
562 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
563 H = DiagonalQuadrature(g_2D,op.quadratureClosure)
564 Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure)
565 @test Hi*H*v ≈ v rtol = 1e-15
566 @test Hi*H*u ≈ u rtol = 1e-15
567 end
568 end
569 end
570 end
554 571
555 @testset "BoundaryOperator" begin 572 @testset "BoundaryOperator" begin
556 closure_stencil = Stencil((0,2), (2.,1.,3.)) 573 closure_stencil = Stencil((0,2), (2.,1.,3.))
557 g_1D = EquidistantGrid(11, 0.0, 1.0) 574 g_1D = EquidistantGrid(11, 0.0, 1.0)
558 g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0)) 575 g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))