Mercurial > repos > public > sbplib_julia

comparison test/testSbpOperators.jl @ 670:538ccbaeb1f8 feature/boundary_quads

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Update quadrature tests
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Sat, 06 Feb 2021 12:04:06 +0100
parents 351937390162
children 1ce3a104afc8 730565f7cc2e
comparison
equal deleted inserted replaced
669:2a95beb9ef1d 670:538ccbaeb1f8
288 @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10 288 @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10
289 @test Dₓₓ*monomials[4] ≈ monomials[2] atol = 5e-10 289 @test Dₓₓ*monomials[4] ≈ monomials[2] atol = 5e-10
290 @test Dₓₓ*v ≈ vₓₓ rtol = 5e-4 norm = l2 290 @test Dₓₓ*v ≈ vₓₓ rtol = 5e-4 norm = l2
291 end 291 end
292 end 292 end
293 293
294 @testset "2D" begin 294 @testset "2D" begin
295 l2(v) = sqrt(prod(spacing(g_2D))*sum(v.^2)); 295 l2(v) = sqrt(prod(spacing(g_2D))*sum(v.^2));
296 binomials = () 296 binomials = ()
297 maxOrder = 4; 297 maxOrder = 4;
298 for i = 0:maxOrder-1 298 for i = 0:maxOrder-1
388 @test L*v ≈ Δv rtol = 5e-4 norm = l2 388 @test L*v ≈ Δv rtol = 5e-4 norm = l2
389 end 389 end
390 end 390 end
391 end 391 end
392 392
393 @testset "DiagonalQuadrature" begin 393 @testset "Quadrature diagonal" begin
394 Lx = π/2. 394 Lx = π/2.
395 Ly = Float64(π) 395 Ly = Float64(π)
396 Lz = 1.
396 g_1D = EquidistantGrid(77, 0.0, Lx) 397 g_1D = EquidistantGrid(77, 0.0, Lx)
397 g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) 398 g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
399 g_3D = EquidistantGrid((10,10, 10), (0.0, 0.0, 0.0), (Lx,Ly,Lz))
398 integral(H,v) = sum(H*v) 400 integral(H,v) = sum(H*v)
399 @testset "Constructors" begin 401 @testset "quadrature" begin
400 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 402 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
401 @testset "1D" begin 403 @testset "1D" begin
402 H = DiagonalQuadrature(g_1D,op.quadratureClosure) 404 H = quadrature(g_1D,op.quadratureClosure)
403 inner_stencil = Stencil((1.,),center=1) 405 inner_stencil = Stencil((1.,),center=1)
404 @test H == Quadrature(g_1D,inner_stencil,op.quadratureClosure) 406 @test H == quadrature(g_1D,inner_stencil,op.quadratureClosure)
405 @test H isa TensorMapping{T,1,1} where T 407 @test H isa TensorMapping{T,1,1} where T
406 end 408 end
407 @testset "1D" begin 409 @testset "2D" begin
408 H = DiagonalQuadrature(g_2D,op.quadratureClosure) 410 H = quadrature(g_2D,op.quadratureClosure)
409 H_x = DiagonalQuadrature(restrict(g_2D,1),op.quadratureClosure) 411 H_x = quadrature(restrict(g_2D,1),op.quadratureClosure)
410 H_y = DiagonalQuadrature(restrict(g_2D,2),op.quadratureClosure) 412 H_y = quadrature(restrict(g_2D,2),op.quadratureClosure)
411 @test H == H_x⊗H_y 413 @test H == H_x⊗H_y
412 @test H isa TensorMapping{T,2,2} where T 414 @test H isa TensorMapping{T,2,2} where T
413 end 415 end
414 end 416 end
415 417
418 @testset "boundary_quadrature" begin
419 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
420 @testset "1D" begin
421 (id_l, id_r) = boundary_identifiers(g_1D)
422 @test boundary_quadrature(g_1D,op.quadratureClosure,id_l) == IdentityMapping{Float64}()
423 @test boundary_quadrature(g_1D,op.quadratureClosure,id_r) == IdentityMapping{Float64}()
424
425 end
426 @testset "2D" begin
427 (id_w, id_e, id_s, id_n) = boundary_identifiers(g_2D)
428 H_x = quadrature(restrict(g_2D,1),op.quadratureClosure)
429 H_y = quadrature(restrict(g_2D,2),op.quadratureClosure)
430 @test boundary_quadrature(g_2D,op.quadratureClosure,id_w) == H_y
431 @test boundary_quadrature(g_2D,op.quadratureClosure,id_e) == H_y
432 @test boundary_quadrature(g_2D,op.quadratureClosure,id_s) == H_x
433 @test boundary_quadrature(g_2D,op.quadratureClosure,id_n) == H_x
434 end
435 @testset "3D" begin
436 (id_w, id_e,
437 id_s, id_n,
438 id_t, id_b) = boundary_identifiers(g_3D)
439 H_xy = quadrature(restrict(g_3D,[1,2]),op.quadratureClosure)
440 H_xz = quadrature(restrict(g_3D,[1,3]),op.quadratureClosure)
441 H_yz = quadrature(restrict(g_3D,[2,3]),op.quadratureClosure)
442 @test boundary_quadrature(g_3D,op.quadratureClosure,id_w) == H_yz
443 @test boundary_quadrature(g_3D,op.quadratureClosure,id_e) == H_yz
444 @test boundary_quadrature(g_3D,op.quadratureClosure,id_s) == H_xz
445 @test boundary_quadrature(g_3D,op.quadratureClosure,id_n) == H_xz
446 @test boundary_quadrature(g_3D,op.quadratureClosure,id_t) == H_xy
447 @test boundary_quadrature(g_3D,op.quadratureClosure,id_b) == H_xy
448 end
449 end
450
416 @testset "Sizes" begin 451 @testset "Sizes" begin
417 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 452 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
418 @testset "1D" begin 453 @testset "1D" begin
419 H = DiagonalQuadrature(g_1D,op.quadratureClosure) 454 H = quadrature(g_1D,op.quadratureClosure)
420 @test domain_size(H) == size(g_1D) 455 @test domain_size(H) == size(g_1D)
421 @test range_size(H) == size(g_1D) 456 @test range_size(H) == size(g_1D)
422 end 457 end
423 @testset "2D" begin 458 @testset "2D" begin
424 H = DiagonalQuadrature(g_2D,op.quadratureClosure) 459 H = quadrature(g_2D,op.quadratureClosure)
425 @test domain_size(H) == size(g_2D) 460 @test domain_size(H) == size(g_2D)
426 @test range_size(H) == size(g_2D) 461 @test range_size(H) == size(g_2D)
427 end 462 end
428 end 463 end
429 464
436 end 471 end
437 u = evalOn(g_1D,x->sin(x)) 472 u = evalOn(g_1D,x->sin(x))
438 473
439 @testset "2nd order" begin 474 @testset "2nd order" begin
440 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) 475 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
441 H = DiagonalQuadrature(g_1D,op.quadratureClosure) 476 H = quadrature(g_1D,op.quadratureClosure)
442 for i = 1:2 477 for i = 1:2
443 @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14 478 @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14
444 end 479 end
445 @test integral(H,u) ≈ 1. rtol = 1e-4 480 @test integral(H,u) ≈ 1. rtol = 1e-4
446 end 481 end
447 482
448 @testset "4th order" begin 483 @testset "4th order" begin
449 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 484 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
450 H = DiagonalQuadrature(g_1D,op.quadratureClosure) 485 H = quadrature(g_1D,op.quadratureClosure)
451 for i = 1:4 486 for i = 1:4
452 @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14 487 @test integral(H,v[i]) ≈ v[i+1][end] - v[i+1][1] rtol = 1e-14
453 end 488 end
454 @test integral(H,u) ≈ 1. rtol = 1e-8 489 @test integral(H,u) ≈ 1. rtol = 1e-8
455 end 490 end
459 b = 2.1 494 b = 2.1
460 v = b*ones(Float64, size(g_2D)) 495 v = b*ones(Float64, size(g_2D))
461 u = evalOn(g_2D,(x,y)->sin(x)+cos(y)) 496 u = evalOn(g_2D,(x,y)->sin(x)+cos(y))
462 @testset "2nd order" begin 497 @testset "2nd order" begin
463 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) 498 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
464 H = DiagonalQuadrature(g_2D,op.quadratureClosure) 499 H = quadrature(g_2D,op.quadratureClosure)
465 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13 500 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13
466 @test integral(H,u) ≈ π rtol = 1e-4 501 @test integral(H,u) ≈ π rtol = 1e-4
467 end 502 end
468 @testset "4th order" begin 503 @testset "4th order" begin
469 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 504 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
470 H = DiagonalQuadrature(g_2D,op.quadratureClosure) 505 H = quadrature(g_2D,op.quadratureClosure)
471 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13 506 @test integral(H,v) ≈ b*Lx*Ly rtol = 1e-13
472 @test integral(H,u) ≈ π rtol = 1e-8 507 @test integral(H,u) ≈ π rtol = 1e-8
473 end 508 end
474 end 509 end
475 end 510 end
519 @testset "1D" begin 554 @testset "1D" begin
520 v = evalOn(g_1D,x->sin(x)) 555 v = evalOn(g_1D,x->sin(x))
521 u = evalOn(g_1D,x->x^3-x^2+1) 556 u = evalOn(g_1D,x->x^3-x^2+1)
522 @testset "2nd order" begin 557 @testset "2nd order" begin
523 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) 558 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
524 H = DiagonalQuadrature(g_1D,op.quadratureClosure) 559 H = quadrature(g_1D,op.quadratureClosure)
525 Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure) 560 Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure)
526 @test Hi*H*v ≈ v rtol = 1e-15 561 @test Hi*H*v ≈ v rtol = 1e-15
527 @test Hi*H*u ≈ u rtol = 1e-15 562 @test Hi*H*u ≈ u rtol = 1e-15
528 end 563 end
529 @testset "4th order" begin 564 @testset "4th order" begin
530 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 565 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
531 H = DiagonalQuadrature(g_1D,op.quadratureClosure) 566 H = quadrature(g_1D,op.quadratureClosure)
532 Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure) 567 Hi = InverseDiagonalQuadrature(g_1D,op.quadratureClosure)
533 @test Hi*H*v ≈ v rtol = 1e-15 568 @test Hi*H*v ≈ v rtol = 1e-15
534 @test Hi*H*u ≈ u rtol = 1e-15 569 @test Hi*H*u ≈ u rtol = 1e-15
535 end 570 end
536 end 571 end
537 @testset "2D" begin 572 @testset "2D" begin
538 v = evalOn(g_2D,(x,y)->sin(x)+cos(y)) 573 v = evalOn(g_2D,(x,y)->sin(x)+cos(y))
539 u = evalOn(g_2D,(x,y)->x*y + x^5 - sqrt(y)) 574 u = evalOn(g_2D,(x,y)->x*y + x^5 - sqrt(y))
540 @testset "2nd order" begin 575 @testset "2nd order" begin
541 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2) 576 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=2)
542 H = DiagonalQuadrature(g_2D,op.quadratureClosure) 577 H = quadrature(g_2D,op.quadratureClosure)
543 Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure) 578 Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure)
544 @test Hi*H*v ≈ v rtol = 1e-15 579 @test Hi*H*v ≈ v rtol = 1e-15
545 @test Hi*H*u ≈ u rtol = 1e-15 580 @test Hi*H*u ≈ u rtol = 1e-15
546 end 581 end
547 @testset "4th order" begin 582 @testset "4th order" begin
548 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4) 583 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
549 H = DiagonalQuadrature(g_2D,op.quadratureClosure) 584 H = quadrature(g_2D,op.quadratureClosure)
550 Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure) 585 Hi = InverseDiagonalQuadrature(g_2D,op.quadratureClosure)
551 @test Hi*H*v ≈ v rtol = 1e-15 586 @test Hi*H*v ≈ v rtol = 1e-15
552 @test Hi*H*u ≈ u rtol = 1e-15 587 @test Hi*H*u ≈ u rtol = 1e-15
553 end 588 end
554 end 589 end