Mercurial > repos > public > sbplib_julia

comparison test/testSbpOperators.jl @ 697:1b3b8f82349e refactor/operator_naming

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Update tests for inner_product and inverse_inner_product
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Sun, 14 Feb 2021 13:49:44 +0100
parents 6ab473e0ea80
children 5ddf28ddee18
comparison
equal deleted inserted replaced
696:0bec3c4e78c0 697:1b3b8f82349e
391 @test L*v ≈ Δv rtol = 5e-4 norm = l2 391 @test L*v ≈ Δv rtol = 5e-4 norm = l2
392 end 392 end
393 end 393 end
394 end 394 end
395 395
396 @testset "Quadrature diagonal" begin 396 @testset "Diagonal-stencil inner_product" begin
397 Lx = π/2. 397 Lx = π/2.
398 Ly = Float64(π) 398 Ly = Float64(π)
399 Lz = 1. 399 Lz = 1.
400 g_1D = EquidistantGrid(77, 0.0, Lx) 400 g_1D = EquidistantGrid(77, 0.0, Lx)
401 g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) 401 g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
402 g_3D = EquidistantGrid((10,10, 10), (0.0, 0.0, 0.0), (Lx,Ly,Lz)) 402 g_3D = EquidistantGrid((10,10, 10), (0.0, 0.0, 0.0), (Lx,Ly,Lz))
403 integral(H,v) = sum(H*v) 403 integral(H,v) = sum(H*v)
404 @testset "quadrature" begin 404 @testset "inner_product" begin
405 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 405 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
406 @testset "0D" begin 406 @testset "0D" begin
407 H = quadrature(EquidistantGrid{Float64}(),op.quadratureClosure) 407 H = inner_product(EquidistantGrid{Float64}(),op.quadratureClosure)
408 @test H == IdentityMapping{Float64}() 408 @test H == IdentityMapping{Float64}()
409 @test H isa TensorMapping{T,0,0} where T 409 @test H isa TensorMapping{T,0,0} where T
410 end 410 end
411 @testset "1D" begin 411 @testset "1D" begin
412 H = quadrature(g_1D,op.quadratureClosure) 412 H = inner_product(g_1D,op.quadratureClosure)
413 inner_stencil = CenteredStencil(1.) 413 inner_stencil = CenteredStencil(1.)
414 @test H == quadrature(g_1D,op.quadratureClosure,inner_stencil) 414 @test H == inner_product(g_1D,op.quadratureClosure,inner_stencil)
415 @test H isa TensorMapping{T,1,1} where T 415 @test H isa TensorMapping{T,1,1} where T
416 end 416 end
417 @testset "2D" begin 417 @testset "2D" begin
418 H = quadrature(g_2D,op.quadratureClosure) 418 H = inner_product(g_2D,op.quadratureClosure)
419 H_x = quadrature(restrict(g_2D,1),op.quadratureClosure) 419 H_x = inner_product(restrict(g_2D,1),op.quadratureClosure)
420 H_y = quadrature(restrict(g_2D,2),op.quadratureClosure) 420 H_y = inner_product(restrict(g_2D,2),op.quadratureClosure)
421 @test H == H_x⊗H_y 421 @test H == H_x⊗H_y
422 @test H isa TensorMapping{T,2,2} where T 422 @test H isa TensorMapping{T,2,2} where T
423 end 423 end
424 end 424 end
425 425
426 @testset "Sizes" begin 426 @testset "Sizes" begin
427 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 427 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
428 @testset "1D" begin 428 @testset "1D" begin
429 H = quadrature(g_1D,op.quadratureClosure) 429 H = inner_product(g_1D,op.quadratureClosure)
430 @test domain_size(H) == size(g_1D) 430 @test domain_size(H) == size(g_1D)
431 @test range_size(H) == size(g_1D) 431 @test range_size(H) == size(g_1D)
432 end 432 end
433 @testset "2D" begin 433 @testset "2D" begin
434 H = quadrature(g_2D,op.quadratureClosure) 434 H = inner_product(g_2D,op.quadratureClosure)
435 @test domain_size(H) == size(g_2D) 435 @test domain_size(H) == size(g_2D)
436 @test range_size(H) == size(g_2D) 436 @test range_size(H) == size(g_2D)
437 end 437 end
438 end 438 end
439 439
446 end 446 end
447 u = evalOn(g_1D,x->sin(x)) 447 u = evalOn(g_1D,x->sin(x))
448 448
449 @testset "2nd order" begin 449 @testset "2nd order" begin
450 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) 450 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
451 H = quadrature(g_1D,op.quadratureClosure) 451 H = inner_product(g_1D,op.quadratureClosure)
452 for i = 1:2 452 for i = 1:2
453 @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14 453 @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14
454 end 454 end
455 @test integral(H,u) ≈ 1. rtol = 1e-4 455 @test integral(H,u) ≈ 1. rtol = 1e-4
456 end 456 end
457 457
458 @testset "4th order" begin 458 @testset "4th order" begin
459 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 459 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
460 H = quadrature(g_1D,op.quadratureClosure) 460 H = inner_product(g_1D,op.quadratureClosure)
461 for i = 1:4 461 for i = 1:4
462 @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14 462 @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14
463 end 463 end
464 @test integral(H,u) ≈ 1. rtol = 1e-8 464 @test integral(H,u) ≈ 1. rtol = 1e-8
465 end 465 end
469 b = 2.1 469 b = 2.1
470 v = b*ones(Float64, size(g_2D)) 470 v = b*ones(Float64, size(g_2D))
471 u = evalOn(g_2D,(x,y)->sin(x)+cos(y)) 471 u = evalOn(g_2D,(x,y)->sin(x)+cos(y))
472 @testset "2nd order" begin 472 @testset "2nd order" begin
473 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) 473 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
474 H = quadrature(g_2D,op.quadratureClosure) 474 H = inner_product(g_2D,op.quadratureClosure)
475 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13 475 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13
476 @test integral(H,u) ≈ π rtol = 1e-4 476 @test integral(H,u) ≈ π rtol = 1e-4
477 end 477 end
478 @testset "4th order" begin 478 @testset "4th order" begin
479 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 479 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
480 H = quadrature(g_2D,op.quadratureClosure) 480 H = inner_product(g_2D,op.quadratureClosure)
481 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13 481 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13
482 @test integral(H,u) ≈ π rtol = 1e-8 482 @test integral(H,u) ≈ π rtol = 1e-8
483 end 483 end
484 end 484 end
485 end 485 end
486 end 486 end
487 487
488 @testset "InverseDiagonalQuadrature" begin 488 @testset "Diagonal-stencil inverse_inner_product" begin
489 Lx = π/2. 489 Lx = π/2.
490 Ly = Float64(π) 490 Ly = Float64(π)
491 g_1D = EquidistantGrid(77, 0.0, Lx) 491 g_1D = EquidistantGrid(77, 0.0, Lx)
492 g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) 492 g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
493 @testset "Constructors" begin 493 @testset "inverse_inner_product" begin
494 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)
495 @testset "1D" begin 495 @testset "1D" begin
496 Hi = InverseDiagonalQuadrature(g_1D, op.quadratureClosure); 496 Hi = inverse_inner_product(g_1D, op.quadratureClosure);
497 inner_stencil = CenteredStencil(1.) 497 inner_stencil = CenteredStencil(1.)
498 closures = () 498 closures = ()
499 for i = 1:length(op.quadratureClosure) 499 for i = 1:length(op.quadratureClosure)
500 closures = (closures...,Stencil(op.quadratureClosure[i].range,1.0./op.quadratureClosure[i].weights)) 500 closures = (closures...,Stencil(op.quadratureClosure[i].range,1.0./op.quadratureClosure[i].weights))
501 end 501 end
502 @test Hi == InverseQuadrature(g_1D,inner_stencil,closures) 502 @test Hi == inverse_inner_product(g_1D,closures,inner_stencil)
503 @test Hi isa TensorMapping{T,1,1} where T 503 @test Hi isa TensorMapping{T,1,1} where T
504 end 504 end
505 @testset "2D" begin 505 @testset "2D" begin
506 Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure) 506 Hi = inverse_inner_product(g_2D,op.quadratureClosure)
507 Hi_x = InverseDiagonalQuadrature(restrict(g_2D,1),op.quadratureClosure) 507 Hi_x = inverse_inner_product(restrict(g_2D,1),op.quadratureClosure)
508 Hi_y = InverseDiagonalQuadrature(restrict(g_2D,2),op.quadratureClosure) 508 Hi_y = inverse_inner_product(restrict(g_2D,2),op.quadratureClosure)
509 @test Hi == Hi_x⊗Hi_y 509 @test Hi == Hi_x⊗Hi_y
510 @test Hi isa TensorMapping{T,2,2} where T 510 @test Hi isa TensorMapping{T,2,2} where T
511 end 511 end
512 end 512 end
513 513
514 @testset "Sizes" begin 514 @testset "Sizes" begin
515 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 515 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
516 @testset "1D" begin 516 @testset "1D" begin
517 Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure) 517 Hi = inverse_inner_product(g_1D,op.quadratureClosure)
518 @test domain_size(Hi) == size(g_1D) 518 @test domain_size(Hi) == size(g_1D)
519 @test range_size(Hi) == size(g_1D) 519 @test range_size(Hi) == size(g_1D)
520 end 520 end
521 @testset "2D" begin 521 @testset "2D" begin
522 Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure) 522 Hi = inverse_inner_product(g_2D,op.quadratureClosure)
523 @test domain_size(Hi) == size(g_2D) 523 @test domain_size(Hi) == size(g_2D)
524 @test range_size(Hi) == size(g_2D) 524 @test range_size(Hi) == size(g_2D)
525 end 525 end
526 end 526 end
527 527
529 @testset "1D" begin 529 @testset "1D" begin
530 v = evalOn(g_1D,x->sin(x)) 530 v = evalOn(g_1D,x->sin(x))
531 u = evalOn(g_1D,x->x^3-x^2+1) 531 u = evalOn(g_1D,x->x^3-x^2+1)
532 @testset "2nd order" begin 532 @testset "2nd order" begin
533 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) 533 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
534 H = quadrature(g_1D,op.quadratureClosure) 534 H = inner_product(g_1D,op.quadratureClosure)
535 Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure) 535 Hi = inverse_inner_product(g_1D,op.quadratureClosure)
536 @test Hi*H*v ≈ v rtol = 1e-15 536 @test Hi*H*v ≈ v rtol = 1e-15
537 @test Hi*H*u ≈ u rtol = 1e-15 537 @test Hi*H*u ≈ u rtol = 1e-15
538 end 538 end
539 @testset "4th order" begin 539 @testset "4th order" begin
540 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 540 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
541 H = quadrature(g_1D,op.quadratureClosure) 541 H = inner_product(g_1D,op.quadratureClosure)
542 Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure) 542 Hi = inverse_inner_product(g_1D,op.quadratureClosure)
543 @test Hi*H*v ≈ v rtol = 1e-15 543 @test Hi*H*v ≈ v rtol = 1e-15
544 @test Hi*H*u ≈ u rtol = 1e-15 544 @test Hi*H*u ≈ u rtol = 1e-15
545 end 545 end
546 end 546 end
547 @testset "2D" begin 547 @testset "2D" begin
548 v = evalOn(g_2D,(x,y)->sin(x)+cos(y)) 548 v = evalOn(g_2D,(x,y)->sin(x)+cos(y))
549 u = evalOn(g_2D,(x,y)->x*y + x^5 - sqrt(y)) 549 u = evalOn(g_2D,(x,y)->x*y + x^5 - sqrt(y))
550 @testset "2nd order" begin 550 @testset "2nd order" begin
551 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) 551 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
552 H = quadrature(g_2D,op.quadratureClosure) 552 H = inner_product(g_2D,op.quadratureClosure)
553 Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure) 553 Hi = inverse_inner_product(g_2D,op.quadratureClosure)
554 @test Hi*H*v ≈ v rtol = 1e-15 554 @test Hi*H*v ≈ v rtol = 1e-15
555 @test Hi*H*u ≈ u rtol = 1e-15 555 @test Hi*H*u ≈ u rtol = 1e-15
556 end 556 end
557 @testset "4th order" begin 557 @testset "4th order" begin
558 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 558 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
559 H = quadrature(g_2D,op.quadratureClosure) 559 H = inner_product(g_2D,op.quadratureClosure)
560 Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure) 560 Hi = inverse_inner_product(g_2D,op.quadratureClosure)
561 @test Hi*H*v ≈ v rtol = 1e-15 561 @test Hi*H*v ≈ v rtol = 1e-15
562 @test Hi*H*u ≈ u rtol = 1e-15 562 @test Hi*H*u ≈ u rtol = 1e-15
563 end 563 end
564 end 564 end
565 end 565 end