Mercurial > repos > public > sbplib_julia

comparison test/testSbpOperators.jl @ 638:08b2c7a2d063 feature/volume_and_boundary_operators

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