Mercurial > repos > public > sbplib_julia

comparison test/SbpOperators/volumeops/laplace/laplace_test.jl @ 866:1784b1c0af3e feature/laplace_opset

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