sbplib_julia: test/testSbpOperators.jl comparison

comparison test/testSbpOperators.jl @ 649:351937390162 feature/volume_and_boundary_operators

Clean up testSbpOperators (remove some redundat tests, remove todos and fix use of Parity)

author	Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date	Fri, 15 Jan 2021 14:09:53 +0100
parents	27b4bcb22154
children	538ccbaeb1f8 e14627e79a54 5ff162f3ed72

comparison

equal deleted inserted replaced

-:d6edde60909b
+:351937390162
 import Sbplib.SbpOperators.Stencil
 import Sbplib.SbpOperators.VolumeOperator
 import Sbplib.SbpOperators.volume_operator
 import Sbplib.SbpOperators.BoundaryOperator
 import Sbplib.SbpOperators.boundary_operator
-import Sbplib.SbpOperators.Parity
+import Sbplib.SbpOperators.even
+import Sbplib.SbpOperators.odd
 @testset "SbpOperators" begin
 @testset "Stencil" begin
 closure_stencils = (Stencil(1/2 .* (1.,1.),center=1),Stencil((0.,1.),center=2))
 g_1D = EquidistantGrid(11,0.,1.)
 g_2D = EquidistantGrid((11,12),(0.,0.),(1.,1.))
 g_3D = EquidistantGrid((11,12,10),(0.,0.,0.),(1.,1.,1.))
 @testset "Constructors" begin
-#TODO: How are even and odd in SbpOperators.Parity exposed? Currently constructing even as Parity(1)
+@testset "1D" begin
-@testset "1D" begin
+op = VolumeOperator(inner_stencil,closure_stencils,(11,),even)
-op = VolumeOperator(inner_stencil,closure_stencils,(11,),Parity(1))
+@test op == VolumeOperator(g_1D,inner_stencil,closure_stencils,even)
-@test op == VolumeOperator(g_1D,inner_stencil,closure_stencils,Parity(1))
+@test op == volume_operator(g_1D,inner_stencil,closure_stencils,even,1)
-@test op == volume_operator(g_1D,inner_stencil,closure_stencils,Parity(1),1)
 @test op isa TensorMapping{T,1,1} where T
 end
 @testset "2D" begin
-op_x = volume_operator(g_2D,inner_stencil,closure_stencils,Parity(1),1)
+op_x = volume_operator(g_2D,inner_stencil,closure_stencils,even,1)
-op_y = volume_operator(g_2D,inner_stencil,closure_stencils,Parity(1),2)
+op_y = volume_operator(g_2D,inner_stencil,closure_stencils,even,2)
 Ix = IdentityMapping{Float64}((11,))
 Iy = IdentityMapping{Float64}((12,))
-@test op_x == VolumeOperator(inner_stencil,closure_stencils,(11,),Parity(1))⊗Iy
+@test op_x == VolumeOperator(inner_stencil,closure_stencils,(11,),even)⊗Iy
-@test op_y == Ix⊗VolumeOperator(inner_stencil,closure_stencils,(12,),Parity(1))
+@test op_y == Ix⊗VolumeOperator(inner_stencil,closure_stencils,(12,),even)
 @test op_x isa TensorMapping{T,2,2} where T
 @test op_y isa TensorMapping{T,2,2} where T
 end
 @testset "3D" begin
-op_x = volume_operator(g_3D,inner_stencil,closure_stencils,Parity(1),1)
+op_x = volume_operator(g_3D,inner_stencil,closure_stencils,even,1)
-op_y = volume_operator(g_3D,inner_stencil,closure_stencils,Parity(1),2)
+op_y = volume_operator(g_3D,inner_stencil,closure_stencils,even,2)
-op_z = volume_operator(g_3D,inner_stencil,closure_stencils,Parity(1),3)
+op_z = volume_operator(g_3D,inner_stencil,closure_stencils,even,3)
 Ix = IdentityMapping{Float64}((11,))
 Iy = IdentityMapping{Float64}((12,))
 Iz = IdentityMapping{Float64}((10,))
-@test op_x == VolumeOperator(inner_stencil,closure_stencils,(11,),Parity(1))⊗Iy⊗Iz
+@test op_x == VolumeOperator(inner_stencil,closure_stencils,(11,),even)⊗Iy⊗Iz
-@test op_y == Ix⊗VolumeOperator(inner_stencil,closure_stencils,(12,),Parity(1))⊗Iz
+@test op_y == Ix⊗VolumeOperator(inner_stencil,closure_stencils,(12,),even)⊗Iz
-@test op_z == Ix⊗Iy⊗VolumeOperator(inner_stencil,closure_stencils,(10,),Parity(1))
+@test op_z == Ix⊗Iy⊗VolumeOperator(inner_stencil,closure_stencils,(10,),even)
 @test op_x isa TensorMapping{T,3,3} where T
 @test op_y isa TensorMapping{T,3,3} where T
 @test op_z isa TensorMapping{T,3,3} where T
 end
 end
 @testset "Sizes" begin
 @testset "1D" begin
-op = volume_operator(g_1D,inner_stencil,closure_stencils,Parity(1),1)
+op = volume_operator(g_1D,inner_stencil,closure_stencils,even,1)
 @test range_size(op) == domain_size(op) == size(g_1D)
 end
 @testset "2D" begin
-op_x = volume_operator(g_2D,inner_stencil,closure_stencils,Parity(1),1)
+op_x = volume_operator(g_2D,inner_stencil,closure_stencils,even,1)
-op_y = volume_operator(g_2D,inner_stencil,closure_stencils,Parity(1),2)
+op_y = volume_operator(g_2D,inner_stencil,closure_stencils,even,2)
 @test range_size(op_y) == domain_size(op_y) ==
 range_size(op_x) == domain_size(op_x) == size(g_2D)
 end
 @testset "3D" begin
-op_x = volume_operator(g_3D,inner_stencil,closure_stencils,Parity(1),1)
+op_x = volume_operator(g_3D,inner_stencil,closure_stencils,even,1)
-op_y = volume_operator(g_3D,inner_stencil,closure_stencils,Parity(1),2)
+op_y = volume_operator(g_3D,inner_stencil,closure_stencils,even,2)
-op_z = volume_operator(g_3D,inner_stencil,closure_stencils,Parity(1),3)
+op_z = volume_operator(g_3D,inner_stencil,closure_stencils,even,3)
 @test range_size(op_z) == domain_size(op_z) ==
 range_size(op_y) == domain_size(op_y) ==
 range_size(op_x) == domain_size(op_x) == size(g_3D)
 end
 end
-# TODO: Test for other dimensions?
+op_x = volume_operator(g_2D,inner_stencil,closure_stencils,even,1)
-op_x = volume_operator(g_2D,inner_stencil,closure_stencils,Parity(1),1)
+op_y = volume_operator(g_2D,inner_stencil,closure_stencils,odd,2)
-op_y = volume_operator(g_2D,inner_stencil,closure_stencils,Parity(-1),2)
 v = zeros(size(g_2D))
 Nx = size(g_2D)[1]
 Ny = size(g_2D)[2]
 for i = 1:Nx
 v[i,:] .= i
 @testset "Application" begin
 @test op_x*v ≈ rx rtol = 1e-14
 @test op_y*v ≈ ry rtol = 1e-14
 end
-# TODO: Test for other dimensions?
 @testset "Regions" begin
 @test (op_x*v)[Index(1,Lower),Index(3,Interior)] ≈ rx[1,3] rtol = 1e-14
 @test (op_x*v)[Index(2,Lower),Index(3,Interior)] ≈ rx[2,3] rtol = 1e-14
 @test (op_x*v)[Index(6,Interior),Index(3,Interior)] ≈ rx[6,3] rtol = 1e-14
 @test (op_x*v)[Index(10,Upper),Index(3,Interior)] ≈ rx[10,3] rtol = 1e-14
 @test_throws BoundsError (op_y*v)[Index(3,Interior),Index(10,Upper)]
 @test_throws BoundsError (op_y*v)[Index(3,Interior),Index(3,Lower)]
 end
-# TODO: Test for other dimensions?
 @testset "Inferred" begin
 @inferred apply(op_x, v,1,1)
 @inferred apply(op_x, v, Index(1,Lower),Index(1,Lower))
 @inferred apply(op_x, v, Index(6,Interior),Index(1,Lower))
 @inferred apply(op_x, v, Index(11,Upper),Index(1,Lower))
 Lx = 3.5
 Ly = 3.
 g_1D = EquidistantGrid(121, 0.0, Lx)
 g_2D = EquidistantGrid((121,123), (0.0, 0.0), (Lx, Ly))
-# TODO: These areant really constructors. Better name?
 @testset "Constructors" begin
 @testset "1D" begin
 Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils)
 @test Dₓₓ == SecondDerivative(g_1D,op.innerStencil,op.closureStencils,1)
 @test Dₓₓ isa VolumeOperator
 # 4th order interior stencil, 2nd order boundary stencil,
 # implies that L*v should be exact for monomials up to order 3.
 @testset "4th order" begin
 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
 Dₓₓ = SecondDerivative(g_1D,op.innerStencil,op.closureStencils)
-# TODO: high tolerances for checking the "exact" differentiation
+# NOTE: high tolerances for checking the "exact" differentiation
 # due to accumulation of round-off errors/cancellation errors?
 @test Dₓₓ*monomials[1] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
 @test Dₓₓ*monomials[2] ≈ zeros(Float64,size(g_1D)...) atol = 5e-10
 @test Dₓₓ*monomials[3] ≈ monomials[1] atol = 5e-10
 @test Dₓₓ*monomials[4] ≈ monomials[2] atol = 5e-10
 @test Dₓₓ*v ≈ vₓₓ rtol = 5e-4 norm = l2
 end
 end
 @testset "2D" begin
 l2(v) = sqrt(prod(spacing(g_2D))*sum(v.^2));
 binomials = ()
 maxOrder = 4;
 for i = 0:maxOrder-1
 # 4th order interior stencil, 2nd order boundary stencil,
 # implies that L*v should be exact for binomials up to order 3.
 @testset "4th order" begin
 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
 Dyy = SecondDerivative(g_2D,op.innerStencil,op.closureStencils,2)
-# TODO: high tolerances for checking the "exact" differentiation
+# NOTE: high tolerances for checking the "exact" differentiation
 # due to accumulation of round-off errors/cancellation errors?
 @test Dyy*binomials[1] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
 @test Dyy*binomials[2] ≈ zeros(Float64,size(g_2D)...) atol = 5e-9
 @test Dyy*binomials[3] ≈ evalOn(g_2D,(x,y)->1.) atol = 5e-9
 @test Dyy*binomials[4] ≈ evalOn(g_2D,(x,y)->y) atol = 5e-9
 end
 end
 @testset "Laplace" begin
 g_1D = EquidistantGrid(101, 0.0, 1.)
-#TODO: It's nice to verify that 3D works somewhere at least, but perhaps should keep 3D tests to a minimum to avoid
-# long run time for test?
 g_3D = EquidistantGrid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.))
-# TODO: These areant really constructors. Better name?
 @testset "Constructors" begin
 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
 @testset "1D" begin
 L = Laplace(g_1D, op.innerStencil, op.closureStencils)
 @test L == SecondDerivative(g_1D, op.innerStencil, op.closureStencils)
 # 4th order interior stencil, 2nd order boundary stencil,
 # implies that L*v should be exact for binomials up to order 3.
 @testset "4th order" begin
 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
 L = Laplace(g_3D,op.innerStencil,op.closureStencils)
-# TODO: high tolerances for checking the "exact" differentiation
+# NOTE: high tolerances for checking the "exact" differentiation
 # due to accumulation of round-off errors/cancellation errors?
 @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
 @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
 @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9
 @test L*polynomials[4] ≈ polynomials[2] atol = 5e-9
 @testset "BoundaryRestriction" begin
 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
 g_1D = EquidistantGrid(11, 0.0, 1.0)
 g_2D = EquidistantGrid((11,15), (0.0, 0.0), (1.0,1.0))
-# TODO: These areant really constructors. Better name?
 @testset "Constructors" begin
 @testset "1D" begin
 e_l = BoundaryRestriction(g_1D,op.eClosure,Lower())
 @test e_l == BoundaryRestriction(g_1D,op.eClosure,CartesianBoundary{1,Lower}())
 @test e_l == BoundaryOperator(g_1D,op.eClosure,Lower())
 u = fill(3.124)
 @test (e_l*v)[] == v[1]
 @test (e_r*v)[] == v[end]
 @test (e_r*v)[1] == v[end]
-@test e_l'*u == [u[]; zeros(10)]
-@test e_r'*u == [zeros(10); u[]]
 end
 @testset "2D" begin
 e_w = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{1,Lower}())
 e_e = BoundaryRestriction(g_2D, op.eClosure, CartesianBoundary{1,Upper}())
 @test e_w*v == v[1,:]
 @test e_e*v == v[end,:]
 @test e_s*v == v[:,1]
 @test e_n*v == v[:,end]
-g_x = rand(11)
-g_y = rand(15)
-G_w = zeros(Float64, (11,15))
-G_w[1,:] = g_y
-G_e = zeros(Float64, (11,15))
-G_e[end,:] = g_y
-G_s = zeros(Float64, (11,15))
-G_s[:,1] = g_x
-G_n = zeros(Float64, (11,15))
-G_n[:,end] = g_x
-@test e_w'*g_y == G_w
-@test e_e'*g_y == G_e
-@test e_s'*g_x == G_s
-@test e_n'*g_x == G_n
 end
 end
 end
 @testset "NormalDerivative" begin
 @test d_e*v ≈ -v∂x[end,:] atol = 1e-13
 @test d_s*v ≈ v∂y[:,1] atol = 1e-13
 @test d_n*v ≈ -v∂y[:,end] atol = 1e-13
 end
 end
+end
-#TODO: Need to check this? Tested in BoundaryOperator. If so, the old test
-# below should be fixed.
+end
-@testset "Transpose" begin
-# op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
-# d_x_l = zeros(Float64, size(g_2D)[1])
-# d_x_u = zeros(Float64, size(g_2D)[1])
-# for i ∈ eachindex(d_x_l)
-#     d_x_l[i] = op.dClosure[i-1]
-#     d_x_u[i] = op.dClosure[length(d_x_u)-i]
-# end
-# d_y_l = zeros(Float64, size(g_2D)[2])
-# d_y_u = zeros(Float64, size(g_2D)[2])
-# for i ∈ eachindex(d_y_l)
-#     d_y_l[i] = op.dClosure[i-1]
-#     d_y_u[i] = op.dClosure[length(d_y_u)-i]
-# end
-# function prod_matrix(x,y)
-#     G = zeros(Float64, length(x), length(y))
-#     for I ∈ CartesianIndices(G)
-#         G[I] = x[I[1]]*y[I[2]]
-#     end
-#
-#     return G
-# end
-# g_x = [1,2,3,4.0,5]
-# g_y = [5,4,3,2,1.0,11]
-#
-# G_w = prod_matrix(d_x_l, g_y)
-# G_e = prod_matrix(d_x_u, g_y)
-# G_s = prod_matrix(g_x, d_y_l)
-# G_n = prod_matrix(g_x, d_y_u)
-#
-# @test d_w'*g_y ≈ G_w
-# @test_broken d_e'*g_y ≈ G_e
-# @test d_s'*g_x ≈ G_s
-# @test_broken d_n'*g_x ≈ G_n
-end
-end
-end

Mercurial > repos > public > sbplib_julia

comparison test/testSbpOperators.jl @ 649:351937390162 feature/volume_and_boundary_operators