Mercurial > repos > public > sbplib_julia

comparison test/testSbpOperators.jl @ 370:8e55dee6a1a1

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge branch refactor/remove_dynamic_size_tensormapping
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 28 Sep 2020 22:56:54 +0200
parents 0844069ab5ff
children de4746d6d126
comparison
equal deleted inserted replaced
366:0af6da238d95 370:8e55dee6a1a1
31 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 31 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
32 L = 3.5 32 L = 3.5
33 g = EquidistantGrid((101,), (0.0,), (L,)) 33 g = EquidistantGrid((101,), (0.0,), (L,))
34 h_inv = inverse_spacing(g) 34 h_inv = inverse_spacing(g)
35 h = 1/h_inv[1]; 35 h = 1/h_inv[1];
36 Dₓₓ = SecondDerivative(h_inv[1],op.innerStencil,op.closureStencils,op.parity) 36 Dₓₓ = SecondDerivative(g,op.innerStencil,op.closureStencils)
37 37
38 f0(x::Float64) = 1. 38 f0(x::Float64) = 1.
39 f1(x::Float64) = x 39 f1(x::Float64) = x
40 f2(x::Float64) = 1/2*x^2 40 f2(x::Float64) = 1/2*x^2
41 f3(x::Float64) = 1/6*x^3 41 f3(x::Float64) = 1/6*x^3
48 v2 = evalOn(g,f2) 48 v2 = evalOn(g,f2)
49 v3 = evalOn(g,f3) 49 v3 = evalOn(g,f3)
50 v4 = evalOn(g,f4) 50 v4 = evalOn(g,f4)
51 v5 = evalOn(g,f5) 51 v5 = evalOn(g,f5)
52 52
53 @test Dₓₓ isa TensorOperator{T,1} where T 53 @test Dₓₓ isa TensorMapping{T,1,1} where T
54 @test Dₓₓ' isa TensorMapping{T,1,1} where T 54 @test Dₓₓ' isa TensorMapping{T,1,1} where T
55 55
56 # TODO: Should perhaps set tolerance level for isapporx instead? 56 # TODO: Should perhaps set tolerance level for isapporx instead?
57 # Are these tolerance levels resonable or should tests be constructed 57 # Are these tolerance levels resonable or should tests be constructed
58 # differently? 58 # differently?
75 @testset "Laplace2D" begin 75 @testset "Laplace2D" begin
76 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 76 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
77 Lx = 1.5 77 Lx = 1.5
78 Ly = 3.2 78 Ly = 3.2
79 g = EquidistantGrid((102,131), (0.0, 0.0), (Lx,Ly)) 79 g = EquidistantGrid((102,131), (0.0, 0.0), (Lx,Ly))
80 80 L = Laplace(g, op.innerStencil, op.closureStencils)
81 h_inv = inverse_spacing(g) 81
82 h = spacing(g)
83 D_xx = SecondDerivative(h_inv[1],op.innerStencil,op.closureStencils,op.parity)
84 D_yy = SecondDerivative(h_inv[2],op.innerStencil,op.closureStencils,op.parity)
85 L = Laplace((D_xx,D_yy))
86 82
87 f0(x::Float64,y::Float64) = 2. 83 f0(x::Float64,y::Float64) = 2.
88 f1(x::Float64,y::Float64) = x+y 84 f1(x::Float64,y::Float64) = x+y
89 f2(x::Float64,y::Float64) = 1/2*x^2 + 1/2*y^2 85 f2(x::Float64,y::Float64) = 1/2*x^2 + 1/2*y^2
90 f3(x::Float64,y::Float64) = 1/6*x^3 + 1/6*y^3 86 f3(x::Float64,y::Float64) = 1/6*x^3 + 1/6*y^3
98 v3 = evalOn(g,f3) 94 v3 = evalOn(g,f3)
99 v4 = evalOn(g,f4) 95 v4 = evalOn(g,f4)
100 v5 = evalOn(g,f5) 96 v5 = evalOn(g,f5)
101 v5ₓₓ = evalOn(g,f5ₓₓ) 97 v5ₓₓ = evalOn(g,f5ₓₓ)
102 98
103 @test L isa TensorOperator{T,2} where T 99 @test L isa TensorMapping{T,2,2} where T
104 @test L' isa TensorMapping{T,2,2} where T 100 @test L' isa TensorMapping{T,2,2} where T
105 101
106 # TODO: Should perhaps set tolerance level for isapporx instead? 102 # TODO: Should perhaps set tolerance level for isapporx instead?
107 # Are these tolerance levels resonable or should tests be constructed 103 # Are these tolerance levels resonable or should tests be constructed
108 # differently? 104 # differently?
116 @test all(abs.(collect(L*v1)) .<= equalitytol) 112 @test all(abs.(collect(L*v1)) .<= equalitytol)
117 @test all(collect(L*v2) .≈ v0) # Seems to be more accurate 113 @test all(collect(L*v2) .≈ v0) # Seems to be more accurate
118 @test all(abs.((collect(L*v3) - v1)) .<= equalitytol) 114 @test all(abs.((collect(L*v3) - v1)) .<= equalitytol)
119 e4 = collect(L*v4) - v2 115 e4 = collect(L*v4) - v2
120 e5 = collect(L*v5) - v5ₓₓ 116 e5 = collect(L*v5) - v5ₓₓ
117
118 h = spacing(g)
121 @test sqrt(prod(h)*sum(collect(e4.^2))) <= accuracytol 119 @test sqrt(prod(h)*sum(collect(e4.^2))) <= accuracytol
122 @test sqrt(prod(h)*sum(collect(e5.^2))) <= accuracytol 120 @test sqrt(prod(h)*sum(collect(e5.^2))) <= accuracytol
123 end 121 end
124 122
125 @testset "DiagonalInnerProduct" begin 123 @testset "DiagonalInnerProduct" begin
126 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 124 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
127 L = 2.3 125 L = 2.3
128 g = EquidistantGrid((77,), (0.0,), (L,)) 126 g = EquidistantGrid((77,), (0.0,), (L,))
129 h = spacing(g) 127 H = DiagonalInnerProduct(g,op.quadratureClosure)
130 H = DiagonalInnerProduct(h[1],op.quadratureClosure)
131 v = ones(Float64, size(g)) 128 v = ones(Float64, size(g))
132 129
133 @test H isa TensorOperator{T,1} where T 130 @test H isa TensorMapping{T,1,1} where T
134 @test H' isa TensorMapping{T,1,1} where T 131 @test H' isa TensorMapping{T,1,1} where T
135 @test sum(collect(H*v)) ≈ L 132 @test sum(collect(H*v)) ≈ L
136 @test collect(H*v) == collect(H'*v) 133 @test collect(H*v) == collect(H'*v)
137 end 134 end
138 135
140 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 137 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
141 Lx = 2.3 138 Lx = 2.3
142 Ly = 5.2 139 Ly = 5.2
143 g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) 140 g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
144 141
145 h = spacing(g) 142 Q = Quadrature(g, op.quadratureClosure)
146 Hx = DiagonalInnerProduct(h[1],op.quadratureClosure);
147 Hy = DiagonalInnerProduct(h[2],op.quadratureClosure);
148 Q = Quadrature((Hx,Hy))
149 143
150 v = ones(Float64, size(g)) 144 v = ones(Float64, size(g))
151 145
152 @test Q isa TensorOperator{T,2} where T 146 @test Q isa TensorMapping{T,2,2} where T
153 @test Q' isa TensorMapping{T,2,2} where T 147 @test Q' isa TensorMapping{T,2,2} where T
154 @test sum(collect(Q*v)) ≈ (Lx*Ly) 148 @test sum(collect(Q*v)) ≈ (Lx*Ly)
155 @test collect(Q*v) == collect(Q'*v) 149 @test collect(Q*v) == collect(Q'*v)
156 end 150 end
157 151
158 @testset "InverseDiagonalInnerProduct" begin 152 @testset "InverseDiagonalInnerProduct" begin
159 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 153 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
160 L = 2.3 154 L = 2.3
161 g = EquidistantGrid((77,), (0.0,), (L,)) 155 g = EquidistantGrid((77,), (0.0,), (L,))
162 h = spacing(g) 156 H = DiagonalInnerProduct(g, op.quadratureClosure)
163 H = DiagonalInnerProduct(h[1],op.quadratureClosure) 157 Hi = InverseDiagonalInnerProduct(g,op.quadratureClosure)
164
165 h_i = inverse_spacing(g)
166 Hi = InverseDiagonalInnerProduct(h_i[1],1 ./ op.quadratureClosure)
167 v = evalOn(g, x->sin(x)) 158 v = evalOn(g, x->sin(x))
168 159
169 @test Hi isa TensorOperator{T,1} where T 160 @test Hi isa TensorMapping{T,1,1} where T
170 @test Hi' isa TensorMapping{T,1,1} where T 161 @test Hi' isa TensorMapping{T,1,1} where T
171 @test collect(Hi*H*v) ≈ v 162 @test collect(Hi*H*v) ≈ v
172 @test collect(Hi*v) == collect(Hi'*v) 163 @test collect(Hi*v) == collect(Hi'*v)
173 end 164 end
174 165
176 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 167 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
177 Lx = 7.3 168 Lx = 7.3
178 Ly = 8.2 169 Ly = 8.2
179 g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) 170 g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
180 171
181 h = spacing(g) 172 Q = Quadrature(g, op.quadratureClosure)
182 Hx = DiagonalInnerProduct(h[1], op.quadratureClosure); 173 Qinv = InverseQuadrature(g, op.quadratureClosure)
183 Hy = DiagonalInnerProduct(h[2], op.quadratureClosure);
184 Q = Quadrature((Hx,Hy))
185
186 hi = inverse_spacing(g)
187 Hix = InverseDiagonalInnerProduct(hi[1], 1 ./ op.quadratureClosure);
188 Hiy = InverseDiagonalInnerProduct(hi[2], 1 ./ op.quadratureClosure);
189 Qinv = InverseQuadrature((Hix,Hiy))
190 v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y) 174 v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
191 175
192 @test Qinv isa TensorOperator{T,2} where T 176 @test Qinv isa TensorMapping{T,2,2} where T
193 @test Qinv' isa TensorMapping{T,2,2} where T 177 @test Qinv' isa TensorMapping{T,2,2} where T
194 @test collect(Qinv*Q*v) ≈ v 178 @test collect(Qinv*(Q*v)) ≈ v
195 @test collect(Qinv*v) == collect(Qinv'*v) 179 @test collect(Qinv*v) == collect(Qinv'*v)
196 end 180 end
197 # 181 #
198 # @testset "BoundaryValue" begin 182 # @testset "BoundaryValue" begin
199 # op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 183 # op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")