Mercurial > repos > public > sbplib_julia

changeset 559:08e27dee76c3 feature/quadrature_as_outer_product

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Restructure and extend tests for InverseDiagonalQuadrature. Also minor cleanup in tests for DiagonalQuadrature.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Sun, 29 Nov 2020 22:42:23 +0100
parents 9b5710ae6587
children 04d7b4eb63ef
files test/testSbpOperators.jl
diffstat 1 files changed, 64 insertions(+), 31 deletions(-) [+]
line wrap: on
line diff
--- a/test/testSbpOperators.jl	Sun Nov 29 22:06:53 2020 +0100
+++ b/test/testSbpOperators.jl	Sun Nov 29 22:42:23 2020 +0100
@@ -124,7 +124,6 @@
         @test H_x == diagonal_quadrature(g_1D,op.quadratureClosure)
         @test H_x isa TensorMapping{T,1,1} where T
         @test H_x' isa TensorMapping{T,1,1} where T
-
         # 2D
         H_xy = diagonal_quadrature(g_2D,op.quadratureClosure)
         @test H_xy isa TensorMappingComposition
@@ -133,30 +132,30 @@
     end
 
     @testset "Sizes" begin
+        # 1D
         H_x = diagonal_quadrature(g_1D,op.quadratureClosure)
-        H_xy = diagonal_quadrature(g_2D,op.quadratureClosure)
-        # 1D
         @test domain_size(H_x) == size(g_1D)
         @test range_size(H_x) == size(g_1D)
         # 2D
+        H_xy = diagonal_quadrature(g_2D,op.quadratureClosure)
         @test domain_size(H_xy) == size(g_2D)
         @test range_size(H_xy) == size(g_2D)
     end
 
     @testset "Application" begin
+        # 1D
         H_x = diagonal_quadrature(g_1D,op.quadratureClosure)
-        H_xy = diagonal_quadrature(g_2D,op.quadratureClosure)
         a = 3.2
-        b = 2.1
         v_1D = a*ones(Float64, size(g_1D))
-        v_2D = b*ones(Float64, size(g_2D))
         u_1D = evalOn(g_1D,x->sin(x))
-        u_2D = evalOn(g_2D,(x,y)->sin(x)+cos(y))
-        # 1D
         @test integral(H_x,v_1D) ≈ a*Lx rtol = 1e-13
         @test integral(H_x,u_1D) ≈ 1. rtol = 1e-8
         @test H_x*v_1D == H_x'*v_1D
         # 2D
+        H_xy = diagonal_quadrature(g_2D,op.quadratureClosure)
+        b = 2.1
+        v_2D = b*ones(Float64, size(g_2D))
+        u_2D = evalOn(g_2D,(x,y)->sin(x)+cos(y))
         @test integral(H_xy,v_2D) ≈ b*Lx*Ly rtol = 1e-13
         @test integral(H_xy,u_2D) ≈ π rtol = 1e-8
         @test H_xy*v_2D ≈ H_xy'*v_2D rtol = 1e-16 #Failed for exact equality. Must differ in operation order for some reason?
@@ -173,7 +172,7 @@
         # op2 = readOperator(sbp_operators_path()*"d2_2nd.txt",sbp_operators_path()*"h_2nd.txt")
         # H2 = diagonal_quadrature(g_1D,op2.quadratureClosure)
         # for i = 1:3
-        #     @test integral(H4,v[i]) ≈ v[i+1] rtol = 1e-14
+        #     @test integral(H2,v[i]) ≈ v[i+1] rtol = 1e-14
         # end
 
         # 4th order
@@ -198,31 +197,65 @@
 
 @testset "InverseDiagonalQuadrature" begin
     op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
-    Lx = 2.3
-    g = EquidistantGrid(77, 0.0, Lx)
-    H = DiagonalQuadrature(g, op.quadratureClosure)
-    Hi = InverseDiagonalQuadrature(g,op.quadratureClosure)
-    v = evalOn(g, x->sin(x))
+    Lx = π/2.
+    Ly = Float64(π)
+    g_1D = EquidistantGrid(77, 0.0, Lx)
+    g_2D = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
+    @testset "Constructors" begin
+        # 1D
+        Hi_x = InverseDiagonalQuadrature(inverse_spacing(g_1D)[1], 1. ./ op.quadratureClosure, size(g_1D));
+        @test Hi_x == InverseDiagonalQuadrature(g_1D,op.quadratureClosure)
+        @test Hi_x == inverse_diagonal_quadrature(g_1D,op.quadratureClosure)
+        @test Hi_x isa TensorMapping{T,1,1} where T
+        @test Hi_x' isa TensorMapping{T,1,1} where T
 
-    @test Hi isa TensorMapping{T,1,1} where T
-    @test Hi' isa TensorMapping{T,1,1} where T
-    @test Hi == inverse_diagonal_quadrature(g,op.quadratureClosure)
-    @test Hi*H*v ≈ v
-    @test Hi*v == Hi'*v
-    @inferred Hi*v
+        # 2D
+        Hi_xy = inverse_diagonal_quadrature(g_2D,op.quadratureClosure)
+        @test Hi_xy isa TensorMappingComposition
+        @test Hi_xy isa TensorMapping{T,2,2} where T
+        @test Hi_xy' isa TensorMapping{T,2,2} where T
+    end
+
+    @testset "Sizes" begin
+        # 1D
+        Hi_x = inverse_diagonal_quadrature(g_1D,op.quadratureClosure)
+        @test domain_size(Hi_x) == size(g_1D)
+        @test range_size(Hi_x) == size(g_1D)
+        # 2D
+        Hi_xy = inverse_diagonal_quadrature(g_2D,op.quadratureClosure)
+        @test domain_size(Hi_xy) == size(g_2D)
+        @test range_size(Hi_xy) == size(g_2D)
+    end
 
-    # Test multiple dimension
-    Ly = 5.2
-    g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
+    @testset "Application" begin
+        # 1D
+        H_x = diagonal_quadrature(g_1D,op.quadratureClosure)
+        Hi_x = inverse_diagonal_quadrature(g_1D,op.quadratureClosure)
+        v_1D = evalOn(g_1D,x->sin(x))
+        u_1D = evalOn(g_1D,x->x^3-x^2+1)
+        @test Hi_x*H_x*v_1D ≈ v_1D rtol = 1e-15
+        @test Hi_x*H_x*u_1D ≈ u_1D rtol = 1e-15
+        @test Hi_x*v_1D == Hi_x'*v_1D
+        # 2D
+        H_xy = diagonal_quadrature(g_2D,op.quadratureClosure)
+        Hi_xy = inverse_diagonal_quadrature(g_2D,op.quadratureClosure)
+        v_2D = evalOn(g_2D,(x,y)->sin(x)+cos(y))
+        u_2D = evalOn(g_2D,(x,y)->x*y + x^5 - sqrt(y))
+        @test Hi_xy*H_xy*v_2D ≈ v_2D rtol = 1e-15
+        @test Hi_xy*H_xy*u_2D ≈ u_2D rtol = 1e-15
+        @test Hi_xy*v_2D ≈ Hi_xy'*v_2D rtol = 1e-16 #Failed for exact equality. Must differ in operation order for some reason?
+    end
 
-    H = diagonal_quadrature(g, op.quadratureClosure)
-    Hi = inverse_diagonal_quadrature(g, op.quadratureClosure)
-    v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
-
-    @test Hi isa TensorMapping{T,2,2} where T
-    @test Hi' isa TensorMapping{T,2,2} where T
-    @test Hi*H*v ≈ v
-    @test_broken Hi*v == Hi'*v # apply_transpose not implemented for InflatedTensorMapping!
+    @testset "Inferred" begin
+        Hi_x = inverse_diagonal_quadrature(g_1D,op.quadratureClosure)
+        Hi_xy = inverse_diagonal_quadrature(g_2D,op.quadratureClosure)
+        v_1D = ones(Float64, size(g_1D))
+        v_2D = ones(Float64, size(g_2D))
+        @inferred Hi_x*v_1D
+        @inferred Hi_x'*v_1D
+        @inferred Hi_xy*v_2D
+        @inferred Hi_xy'*v_2D
+    end
 end
 #
 # @testset "BoundaryValue" begin