Mercurial > repos > public > sbplib_julia

comparison test/SbpOperators/volumeops/laplace/laplace_test.jl @ 750:f88b2117dc69 feature/laplace_opset

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge in default
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Fri, 19 Mar 2021 16:52:53 +0100
parents 6114274447f5
children f94feb005e7d
comparison
equal deleted inserted replaced
723:c16abc564b82 750:f88b2117dc69
1 using Test
2
3 using Sbplib.SbpOperators
4 using Sbplib.Grids
5 using Sbplib.LazyTensors
6 using Sbplib.RegionIndices
7
8 """
9 cmp_fields(s1,s2)
10
11 Compares the fields of two structs s1, s2, using the == operator.
12 """
13 function cmp_fields(s1::T,s2::T) where T
14 f = fieldnames(T)
15 return getfield.(Ref(s1),f) == getfield.(Ref(s2),f)
16 end
17
18 @testset "Laplace" begin
19 g_1D = EquidistantGrid(101, 0.0, 1.)
20 g_3D = EquidistantGrid((51,101,52), (0.0, -1.0, 0.0), (1., 1., 1.))
21 op = read_D2_operator(sbp_operators_path()*"standard_diagonal.toml"; order=4)
22 @testset "Constructors" begin
23
24 @testset "1D" begin
25 # Create all tensor mappings included in Laplace
26 Δ = laplace(g_1D, op.innerStencil, op.closureStencils)
27 H = inner_product(g_1D, op.quadratureClosure)
28 Hi = inverse_inner_product(g_1D, op.quadratureClosure)
29
30 (id_l, id_r) = boundary_identifiers(g_1D)
31
32 e_l = boundary_restriction(g_1D,op.eClosure,id_l)
33 e_r = boundary_restriction(g_1D,op.eClosure,id_r)
34 e_dict = Dict(Pair(id_l,e_l),Pair(id_r,e_r))
35
36 d_l = normal_derivative(g_1D,op.dClosure,id_l)
37 d_r = normal_derivative(g_1D,op.dClosure,id_r)
38 d_dict = Dict(Pair(id_l,d_l),Pair(id_r,d_r))
39
40 H_l = inner_product(boundary_grid(g_1D,id_l),op.quadratureClosure)
41 H_r = inner_product(boundary_grid(g_1D,id_r),op.quadratureClosure)
42 Hb_dict = Dict(Pair(id_l,H_l),Pair(id_r,H_r))
43
44 L = Laplace(g_1D, sbp_operators_path()*"standard_diagonal.toml"; order=4)
45 @test cmp_fields(L,Laplace(Δ,H,Hi,e_dict,d_dict,Hb_dict))
46 @test L isa TensorMapping{T,1,1} where T
47 @inferred Laplace(Δ,H,Hi,e_dict,d_dict,Hb_dict)
48 end
49 @testset "3D" begin
50 # Create all tensor mappings included in Laplace
51 Δ = laplace(g_3D, op.innerStencil, op.closureStencils)
52 H = inner_product(g_3D, op.quadratureClosure)
53 Hi = inverse_inner_product(g_3D, op.quadratureClosure)
54
55 (id_l, id_r, id_s, id_n, id_b, id_t) = boundary_identifiers(g_3D)
56
57 e_l = boundary_restriction(g_3D,op.eClosure,id_l)
58 e_r = boundary_restriction(g_3D,op.eClosure,id_r)
59 e_s = boundary_restriction(g_3D,op.eClosure,id_s)
60 e_n = boundary_restriction(g_3D,op.eClosure,id_n)
61 e_b = boundary_restriction(g_3D,op.eClosure,id_b)
62 e_t = boundary_restriction(g_3D,op.eClosure,id_t)
63 e_dict = Dict(Pair(id_l,e_l),Pair(id_r,e_r),
64 Pair(id_s,e_s),Pair(id_n,e_n),
65 Pair(id_b,e_b),Pair(id_t,e_t))
66
67 d_l = normal_derivative(g_3D,op.dClosure,id_l)
68 d_r = normal_derivative(g_3D,op.dClosure,id_r)
69 d_s = normal_derivative(g_3D,op.dClosure,id_s)
70 d_n = normal_derivative(g_3D,op.dClosure,id_n)
71 d_b = normal_derivative(g_3D,op.dClosure,id_b)
72 d_t = normal_derivative(g_3D,op.dClosure,id_t)
73 d_dict = Dict(Pair(id_l,d_l),Pair(id_r,d_r),
74 Pair(id_s,d_s),Pair(id_n,d_n),
75 Pair(id_b,d_b),Pair(id_t,d_t))
76
77 H_l = inner_product(boundary_grid(g_3D,id_l),op.quadratureClosure)
78 H_r = inner_product(boundary_grid(g_3D,id_r),op.quadratureClosure)
79 H_s = inner_product(boundary_grid(g_3D,id_s),op.quadratureClosure)
80 H_n = inner_product(boundary_grid(g_3D,id_n),op.quadratureClosure)
81 H_b = inner_product(boundary_grid(g_3D,id_b),op.quadratureClosure)
82 H_t = inner_product(boundary_grid(g_3D,id_t),op.quadratureClosure)
83 Hb_dict = Dict(Pair(id_l,H_l),Pair(id_r,H_r),
84 Pair(id_s,H_s),Pair(id_n,H_n),
85 Pair(id_b,H_b),Pair(id_t,H_t))
86
87 L = Laplace(g_3D, sbp_operators_path()*"standard_diagonal.toml"; order=4)
88 @test cmp_fields(L,Laplace(Δ,H,Hi,e_dict,d_dict,Hb_dict))
89 @test L isa TensorMapping{T,3,3} where T
90 @inferred Laplace(Δ,H,Hi,e_dict,d_dict,Hb_dict)
91 end
92 end
93
94 @testset "laplace" begin
95 @testset "1D" begin
96 L = laplace(g_1D, op.innerStencil, op.closureStencils)
97 @test L == second_derivative(g_1D, op.innerStencil, op.closureStencils)
98 @test L isa TensorMapping{T,1,1} where T
99 end
100 @testset "3D" begin
101 L = laplace(g_3D, op.innerStencil, op.closureStencils)
102 @test L isa TensorMapping{T,3,3} where T
103 Dxx = second_derivative(g_3D, op.innerStencil, op.closureStencils,1)
104 Dyy = second_derivative(g_3D, op.innerStencil, op.closureStencils,2)
105 Dzz = second_derivative(g_3D, op.innerStencil, op.closureStencils,3)
106 @test L == Dxx + Dyy + Dzz
107 @test L isa TensorMapping{T,3,3} where T
108 end
109 end
110
111 @testset "inner_product" begin
112 L = Laplace(g_3D, sbp_operators_path()*"standard_diagonal.toml"; order=4)
113 @test inner_product(L) == inner_product(g_3D,op.quadratureClosure)
114 end
115
116 @testset "inverse_inner_product" begin
117 L = Laplace(g_3D, sbp_operators_path()*"standard_diagonal.toml"; order=4)
118 @test inverse_inner_product(L) == inverse_inner_product(g_3D,op.quadratureClosure)
119 end
120
121 @testset "boundary_restriction" begin
122 L = Laplace(g_3D, sbp_operators_path()*"standard_diagonal.toml"; order=4)
123 id_l = CartesianBoundary{1,Lower}()
124 id_r = CartesianBoundary{1,Upper}()
125 id_s = CartesianBoundary{2,Lower}()
126 id_n = CartesianBoundary{2,Upper}()
127 id_b = CartesianBoundary{3,Lower}()
128 id_t = CartesianBoundary{3,Upper}()
129 @test boundary_restriction(L,id_l) == boundary_restriction(g_3D,op.eClosure,id_l)
130 @test boundary_restriction(L,id_r) == boundary_restriction(g_3D,op.eClosure,id_r)
131 @test boundary_restriction(L,id_s) == boundary_restriction(g_3D,op.eClosure,id_s)
132 @test boundary_restriction(L,id_n) == boundary_restriction(g_3D,op.eClosure,id_n)
133 @test boundary_restriction(L,id_b) == boundary_restriction(g_3D,op.eClosure,id_b)
134 @test boundary_restriction(L,id_t) == boundary_restriction(g_3D,op.eClosure,id_t)
135 end
136
137 @testset "normal_derivative" begin
138 L = Laplace(g_3D, sbp_operators_path()*"standard_diagonal.toml"; order=4)
139 id_l = CartesianBoundary{1,Lower}()
140 id_r = CartesianBoundary{1,Upper}()
141 id_s = CartesianBoundary{2,Lower}()
142 id_n = CartesianBoundary{2,Upper}()
143 id_b = CartesianBoundary{3,Lower}()
144 id_t = CartesianBoundary{3,Upper}()
145 @test normal_derivative(L,id_l) == normal_derivative(g_3D,op.dClosure,id_l)
146 @test normal_derivative(L,id_r) == normal_derivative(g_3D,op.dClosure,id_r)
147 @test normal_derivative(L,id_s) == normal_derivative(g_3D,op.dClosure,id_s)
148 @test normal_derivative(L,id_n) == normal_derivative(g_3D,op.dClosure,id_n)
149 @test normal_derivative(L,id_b) == normal_derivative(g_3D,op.dClosure,id_b)
150 @test normal_derivative(L,id_t) == normal_derivative(g_3D,op.dClosure,id_t)
151 end
152
153 @testset "boundary_quadrature" begin
154 L = Laplace(g_3D, sbp_operators_path()*"standard_diagonal.toml"; order=4)
155 id_l = CartesianBoundary{1,Lower}()
156 id_r = CartesianBoundary{1,Upper}()
157 id_s = CartesianBoundary{2,Lower}()
158 id_n = CartesianBoundary{2,Upper}()
159 id_b = CartesianBoundary{3,Lower}()
160 id_t = CartesianBoundary{3,Upper}()
161 @test boundary_quadrature(L,id_l) == inner_product(boundary_grid(g_3D,id_l),op.quadratureClosure)
162 @test boundary_quadrature(L,id_r) == inner_product(boundary_grid(g_3D,id_r),op.quadratureClosure)
163 @test boundary_quadrature(L,id_s) == inner_product(boundary_grid(g_3D,id_s),op.quadratureClosure)
164 @test boundary_quadrature(L,id_n) == inner_product(boundary_grid(g_3D,id_n),op.quadratureClosure)
165 @test boundary_quadrature(L,id_b) == inner_product(boundary_grid(g_3D,id_b),op.quadratureClosure)
166 @test boundary_quadrature(L,id_t) == inner_product(boundary_grid(g_3D,id_t),op.quadratureClosure)
167 end
168
169 # Exact differentiation is measured point-wise. In other cases
170 # the error is measured in the l2-norm.
171 @testset "Accuracy" begin
172 l2(v) = sqrt(prod(spacing(g_3D))*sum(v.^2));
173 polynomials = ()
174 maxOrder = 4;
175 for i = 0:maxOrder-1
176 f_i(x,y,z) = 1/factorial(i)*(y^i + x^i + z^i)
177 polynomials = (polynomials...,evalOn(g_3D,f_i))
178 end
179 v = evalOn(g_3D, (x,y,z) -> sin(x) + cos(y) + exp(z))
180 Δv = evalOn(g_3D,(x,y,z) -> -sin(x) - cos(y) + exp(z))
181
182 # 2nd order interior stencil, 1st order boundary stencil,
183 # implies that L*v should be exact for binomials up to order 2.
184 @testset "2nd order" begin
185 L = Laplace(g_3D, sbp_operators_path()*"standard_diagonal.toml"; order=2)
186 @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
187 @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
188 @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9
189 @test L*v ≈ Δv rtol = 5e-2 norm = l2
190 end
191
192 # 4th order interior stencil, 2nd order boundary stencil,
193 # implies that L*v should be exact for binomials up to order 3.
194 @testset "4th order" begin
195 L = Laplace(g_3D, sbp_operators_path()*"standard_diagonal.toml"; order=4)
196 # NOTE: high tolerances for checking the "exact" differentiation
197 # due to accumulation of round-off errors/cancellation errors?
198 @test L*polynomials[1] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
199 @test L*polynomials[2] ≈ zeros(Float64, size(g_3D)...) atol = 5e-9
200 @test L*polynomials[3] ≈ polynomials[1] atol = 5e-9
201 @test L*polynomials[4] ≈ polynomials[2] atol = 5e-9
202 @test L*v ≈ Δv rtol = 5e-4 norm = l2
203 end
204 end
205 end