Mercurial > repos > public > sbplib_julia

comparison DiffOps/test/runtests.jl @ 282:ce6a2f3f732a boundary_conditions

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make Laplace a TensorOperator and add tests. NOTE: Two of the tests for Laplace2D are currently failing.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 09 Jan 2020 10:54:24 +0100
parents f1e90a92ad74
children 12a12a5cd973
comparison
equal deleted inserted replaced
281:1eefaefdd0c7 282:ce6a2f3f732a
3 using Grids 3 using Grids
4 using SbpOperators 4 using SbpOperators
5 using RegionIndices 5 using RegionIndices
6 using LazyTensors 6 using LazyTensors
7 7
8 @testset "Laplace2D" begin
9 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
10 Lx = 3.5
11 Ly = 7.2
12 g = EquidistantGrid((21,42), (0.0, 0.0), (Lx,Ly))
13 L = Laplace(g, 1., op)
14
15 f0(x::Float64,y::Float64) = 2.
16 f1(x::Float64,y::Float64) = x+y
17 f2(x::Float64,y::Float64) = 1/2*x^2 + 1/2*y^2
18 f3(x::Float64,y::Float64) = 1/6*x^3 + 1/6*y^3
19 f4(x::Float64,y::Float64) = 1/24*x^4 + 1/24*y^4
20 f5(x::Float64,y::Float64) = x^5 + 2*y^3 + 3/2*x^2 + y + 7
21 f5ₓₓ(x::Float64,y::Float64) = 20*x^3 + 12*y + 3
22
23 v0 = evalOn(g,f0)
24 v1 = evalOn(g,f1)
25 v2 = evalOn(g,f2)
26 v3 = evalOn(g,f3)
27 v4 = evalOn(g,f4)
28 v5 = evalOn(g,f5)
29 v5ₓₓ = evalOn(g,f5ₓₓ)
30
31 @test L isa TensorOperator{T,2} where T
32 @test L' isa TensorMapping{T,2,2} where T
33 # TODO: Should perhaps set tolerance level for isapporx instead?
34 equalitytol = 0.5*1e-11
35 accuracytol = 1e-4
36 @test all(collect(L*v0) .<= equalitytol)
37 @test all(collect(L*v1) .<= equalitytol)
38 @test all(collect(L*v2) .≈ v0) # Seems to be more accurate
39 @test all((collect(L*v3) - v1) .<= equalitytol)
40 @show maximum(abs.(collect(L*v4) - v2))
41 @show maximum(abs.(collect(L*v5) - v5ₓₓ))
42 @test_broken all((collect(L*v4) - v2) .<= accuracytol) #TODO: Should be equality?
43 @test_broken all((collect(L*v5) - v5ₓₓ) .<= accuracytol) #TODO: This is just wrong
44 end
45
8 @testset "Quadrature" begin 46 @testset "Quadrature" begin
9 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt") 47 op = readOperator(sbp_operators_path()*"d2_4th.txt",sbp_operators_path()*"h_4th.txt")
10 Lx = 2.3 48 Lx = 2.3
11 Ly = 5.2 49 Ly = 5.2
12 g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) 50 g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
13 H = Quadrature(op,g) 51 H = Quadrature(op,g)
14 v = ones(Float64, size(g)) 52 v = ones(Float64, size(g))
15 53
16 @test H isa TensorMapping{T,2,2} where T 54 @test H isa TensorOperator{T,2} where T
17 @test H' isa TensorMapping{T,2,2} where T 55 @test H' isa TensorMapping{T,2,2} where T
18 @test sum(collect(H*v)) ≈ (Lx*Ly) 56 @test sum(collect(H*v)) ≈ (Lx*Ly)
19 @test collect(H*v) == collect(H'*v) 57 @test collect(H*v) == collect(H'*v)
20 end 58 end
21 59
26 g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly)) 64 g = EquidistantGrid((77,66), (0.0, 0.0), (Lx,Ly))
27 H = Quadrature(op,g) 65 H = Quadrature(op,g)
28 Hinv = InverseQuadrature(op,g) 66 Hinv = InverseQuadrature(op,g)
29 v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y) 67 v = evalOn(g, (x,y)-> x^2 + (y-1)^2 + x*y)
30 68
31 @test Hinv isa TensorMapping{T,2,2} where T 69 @test Hinv isa TensorOperator{T,2} where T
32 @test Hinv' isa TensorMapping{T,2,2} where T 70 @test Hinv' isa TensorMapping{T,2,2} where T
33 @test collect(Hinv*H*v) ≈ v 71 @test collect(Hinv*H*v) ≈ v
34 @test collect(Hinv*v) == collect(Hinv'*v) 72 @test collect(Hinv*v) == collect(Hinv'*v)
35 end 73 end
36 74
139 for i ∈ eachindex(d_y_l) 177 for i ∈ eachindex(d_y_l)
140 d_y_l[i] = op.dClosure[i-1] 178 d_y_l[i] = op.dClosure[i-1]
141 d_y_u[i] = -op.dClosure[length(d_y_u)-i] 179 d_y_u[i] = -op.dClosure[length(d_y_u)-i]
142 end 180 end
143 181
144 function ❓(x,y) 182 function prod_matrix(x,y)
145 G = zeros(Float64, length(x), length(y)) 183 G = zeros(Float64, length(x), length(y))
146 for I ∈ CartesianIndices(G) 184 for I ∈ CartesianIndices(G)
147 G[I] = x[I[1]]*y[I[2]] 185 G[I] = x[I[1]]*y[I[2]]
148 end 186 end
149 187
151 end 189 end
152 190
153 g_x = [1,2,3,4.0,5] 191 g_x = [1,2,3,4.0,5]
154 g_y = [5,4,3,2,1.0,11] 192 g_y = [5,4,3,2,1.0,11]
155 193
156 G_w = ❓(d_x_l, g_y) 194 G_w = prod_matrix(d_x_l, g_y)
157 G_e = ❓(d_x_u, g_y) 195 G_e = prod_matrix(d_x_u, g_y)
158 G_s = ❓(g_x, d_y_l) 196 G_s = prod_matrix(g_x, d_y_l)
159 G_n = ❓(g_x, d_y_u) 197 G_n = prod_matrix(g_x, d_y_u)
160 198
161 199
162 @test size(d_w*g_y) == (5,6) 200 @test size(d_w*g_y) == (5,6)
163 @test size(d_e*g_y) == (5,6) 201 @test size(d_e*g_y) == (5,6)
164 @test size(d_s*g_x) == (5,6) 202 @test size(d_s*g_x) == (5,6)