Mercurial > repos > public > sbplib_julia

comparison test/Grids/mapped_grid_test.jl @ 1751:f3d7e2d7a43f feature/sbp_operators/laplace_curvilinear

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge feature/grids/manifolds
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 11 Sep 2024 16:26:19 +0200
parents 03894fd7b132
children 819ab806960f
comparison
equal deleted inserted replaced
1731:3684db043add 1751:f3d7e2d7a43f
1 using Sbplib.Grids 1 using Diffinitive.Grids
2 using Sbplib.RegionIndices 2 using Diffinitive.RegionIndices
3 using Test 3 using Test
4 using StaticArrays 4 using StaticArrays
5 using LinearAlgebra 5 using LinearAlgebra
6 6
7 7
26 26
27 return x̄, J 27 return x̄, J
28 end 28 end
29 29
30 @testset "MappedGrid" begin 30 @testset "MappedGrid" begin
31 lg = equidistant_grid((0,0), (1,1), 11, 11) # TODO: Change dims of the grid to be different 31 @testset "Constructor" begin
32 x̄ = map(ξ̄ -> 2ξ̄, lg) 32 lg = equidistant_grid((0,0), (1,1), 11, 21)
33 J = map(ξ̄ -> @SArray(fill(2., 2, 2)), lg) 33
34 mg = MappedGrid(lg, x̄, J) 34 x̄ = map(ξ̄ -> 2ξ̄, lg)
35 35 J = map(ξ̄ -> @SArray(fill(2., 2, 2)), lg)
36 # TODO: Test constructor for different dims of range and domain for the coordinates 36 mg = MappedGrid(lg, x̄, J)
37 # TODO: Test constructor with different type than TensorGrid. a dummy type? 37
38 38 @test mg isa Grid{SVector{2, Float64},2}
39 @test_broken false # @test_throws ArgumentError("Sizes must match") MappedGrid(lg, map(ξ̄ -> @SArray[ξ̄[1], ξ̄[2], -ξ̄[1]], lg), rand(SMatrix{2,3,Float64},15,11)) 39 @test jacobian(mg) isa Array{<:AbstractMatrix}
40 40 @test logicalgrid(mg) isa Grid
41 41
42 @test mg isa Grid{SVector{2, Float64},2} 42 @test collect(mg) == x̄
43 43 @test jacobian(mg) == J
44 @test jacobian(mg) isa Array{<:AbstractMatrix} 44 @test logicalgrid(mg) == lg
45 @test logicalgrid(mg) isa Grid 45
46
47 x̄ = map(ξ̄ -> @SVector[ξ̄[1],ξ̄[2], ξ̄[1] + ξ̄[2]], lg)
48 J = map(ξ̄ -> @SMatrix[1 0; 0 1; 1 1], lg)
49 mg = MappedGrid(lg, x̄, J)
50
51 @test mg isa Grid{SVector{3, Float64},2}
52 @test jacobian(mg) isa Array{<:AbstractMatrix}
53 @test logicalgrid(mg) isa Grid
54
55 @test collect(mg) == x̄
56 @test jacobian(mg) == J
57 @test logicalgrid(mg) == lg
58
59 sz1 = (10,11)
60 sz2 = (10,12)
61 @test_throws ArgumentError("Sizes must match") MappedGrid(
62 equidistant_grid((0,0), (1,1), sz2...),
63 rand(SVector{2},sz1...),
64 rand(SMatrix{2,2},sz1...),
65 )
66
67 @test_throws ArgumentError("Sizes must match") MappedGrid(
68 equidistant_grid((0,0), (1,1), sz1...),
69 rand(SVector{2},sz2...),
70 rand(SMatrix{2,2},sz1...),
71 )
72
73 @test_throws ArgumentError("Sizes must match") MappedGrid(
74 equidistant_grid((0,0), (1,1), sz1...),
75 rand(SVector{2},sz1...),
76 rand(SMatrix{2,2},sz2...),
77 )
78
79 err_str = "The size of the jacobian must match the dimensions of the grid and coordinates"
80 @test_throws ArgumentError(err_str) MappedGrid(
81 equidistant_grid((0,0), (1,1), 10, 11),
82 rand(SVector{3}, 10, 11),
83 rand(SMatrix{3,4}, 10, 11),
84 )
85
86 @test_throws ArgumentError(err_str) MappedGrid(
87 equidistant_grid((0,0), (1,1), 10, 11),
88 rand(SVector{3}, 10, 11),
89 rand(SMatrix{4,2}, 10, 11),
90 )
91 end
46 92
47 @testset "Indexing Interface" begin 93 @testset "Indexing Interface" begin
94 lg = equidistant_grid((0,0), (1,1), 11, 21)
95 x̄ = map(ξ̄ -> 2ξ̄, lg)
96 J = map(ξ̄ -> @SArray(fill(2., 2, 2)), lg)
48 mg = MappedGrid(lg, x̄, J) 97 mg = MappedGrid(lg, x̄, J)
49 @test mg[1,1] == [0.0, 0.0] 98 @test mg[1,1] == [0.0, 0.0]
50 @test mg[4,2] == [0.6, 0.2] 99 @test mg[4,2] == [0.6, 0.1]
51 @test mg[6,10] == [1., 1.8] 100 @test mg[6,10] == [1., 0.9]
52 101
53 @test mg[begin, begin] == [0.0, 0.0] 102 @test mg[begin, begin] == [0.0, 0.0]
54 @test mg[end,end] == [2.0, 2.0] 103 @test mg[end,end] == [2.0, 2.0]
55 @test mg[begin,end] == [0., 2.] 104 @test mg[begin,end] == [0., 2.]
56 105
57 @test eachindex(mg) == CartesianIndices((11,11)) 106 @test axes(mg) == (1:11, 1:21)
58 107
59 @testset "cartesian indexing" begin 108 @testset "cartesian indexing" begin
60 cases = [ 109 cases = [
61 (1,1) , 110 (1,1) ,
62 (3,5) , 111 (3,5) ,
69 @test mg[CartesianIndex(is...)] == mg[is...] 118 @test mg[CartesianIndex(is...)] == mg[is...]
70 end 119 end
71 end 120 end
72 121
73 @testset "eachindex" begin 122 @testset "eachindex" begin
74 @test eachindex(mg) == CartesianIndices((11,11)) 123 @test eachindex(mg) == CartesianIndices((11,21))
75 end 124 end
76 125
77 @testset "firstindex" begin 126 @testset "firstindex" begin
78 @test firstindex(mg, 1) == 1 127 @test firstindex(mg, 1) == 1
79 @test firstindex(mg, 2) == 1 128 @test firstindex(mg, 2) == 1
80 end 129 end
81 130
82 @testset "lastindex" begin 131 @testset "lastindex" begin
83 @test lastindex(mg, 1) == 11 132 @test lastindex(mg, 1) == 11
84 @test lastindex(mg, 2) == 11 133 @test lastindex(mg, 2) == 21
85 end 134 end
86 end 135 end
87 # TODO: Test with different types of logical grids
88 136
89 @testset "Iterator interface" begin 137 @testset "Iterator interface" begin
138 lg = equidistant_grid((0,0), (1,1), 11, 21)
139 x̄ = map(ξ̄ -> 2ξ̄, lg)
140 J = map(ξ̄ -> @SArray(fill(2., 2, 2)), lg)
141
142 mg = MappedGrid(lg, x̄, J)
143
144 lg2 = equidistant_grid((0,0), (1,1), 15, 11)
90 sg = MappedGrid( 145 sg = MappedGrid(
91 equidistant_grid((0,0), (1,1), 15, 11), 146 equidistant_grid((0,0), (1,1), 15, 11),
92 map(ξ̄ -> @SArray[ξ̄[1], ξ̄[2], -ξ̄[1]], lg), rand(SMatrix{2,3,Float64},15,11) 147 map(ξ̄ -> @SArray[ξ̄[1], ξ̄[2], -ξ̄[1]], lg2), rand(SMatrix{3,2,Float64},15,11)
93 ) 148 )
94 149
95 @test eltype(mg) == SVector{2,Float64} 150 @test eltype(mg) == SVector{2,Float64}
96 @test eltype(sg) == SVector{3,Float64} 151 @test eltype(sg) == SVector{3,Float64}
97 152
98 @test eltype(typeof(mg)) == SVector{2,Float64} 153 @test eltype(typeof(mg)) == SVector{2,Float64}
99 @test eltype(typeof(sg)) == SVector{3,Float64} 154 @test eltype(typeof(sg)) == SVector{3,Float64}
100 155
101 @test size(mg) == (11,11) 156 @test size(mg) == (11,21)
102 @test size(sg) == (15,11) 157 @test size(sg) == (15,11)
103 158
104 @test size(mg,2) == 11 159 @test size(mg,2) == 21
105 @test size(sg,2) == 11 160 @test size(sg,2) == 11
106 161
107 @test length(mg) == 121 162 @test length(mg) == 231
108 @test length(sg) == 165 163 @test length(sg) == 165
109 164
110 @test Base.IteratorSize(mg) == Base.HasShape{2}() 165 @test Base.IteratorSize(mg) == Base.HasShape{2}()
111 @test Base.IteratorSize(typeof(mg)) == Base.HasShape{2}() 166 @test Base.IteratorSize(typeof(mg)) == Base.HasShape{2}()
112 167
125 180
126 @test collect(mg) == 2 .* lg 181 @test collect(mg) == 2 .* lg
127 end 182 end
128 183
129 @testset "Base" begin 184 @testset "Base" begin
185 lg = equidistant_grid((0,0), (1,1), 11, 21)
186 x̄ = map(ξ̄ -> 2ξ̄, lg)
187 J = map(ξ̄ -> @SArray(fill(2., 2, 2)), lg)
188 mg = MappedGrid(lg, x̄, J)
189
130 @test ndims(mg) == 2 190 @test ndims(mg) == 2
131 end 191 end
132 192
193 @testset "==" begin
194 sz = (15,11)
195 lg = equidistant_grid((0,0), (1,1), sz...)
196 x = rand(SVector{3,Float64}, sz...)
197 J = rand(SMatrix{3,2,Float64}, sz...)
198
199 sg = MappedGrid(lg, x, J)
200
201 sg1 = MappedGrid(equidistant_grid((0,0), (1,1), sz...), copy(x), copy(J))
202
203 sz2 = (15,12)
204 lg2 = equidistant_grid((0,0), (1,1), sz2...)
205 x2 = rand(SVector{3,Float64}, sz2...)
206 J2 = rand(SMatrix{3,2,Float64}, sz2...)
207 sg2 = MappedGrid(lg2, x2, J2)
208
209 sg3 = MappedGrid(lg, rand(SVector{3,Float64}, sz...), J)
210 sg4 = MappedGrid(lg, x, rand(SMatrix{3,2,Float64}, sz...))
211
212 @test sg == sg1
213 @test sg != sg2 # Different size
214 @test sg != sg3 # Different coordinates
215 @test sg != sg4 # Different jacobian
216 end
217
133 @testset "boundary_identifiers" begin 218 @testset "boundary_identifiers" begin
219 lg = equidistant_grid((0,0), (1,1), 11, 15)
220 x̄ = map(ξ̄ -> 2ξ̄, lg)
221 J = map(ξ̄ -> @SArray(fill(2., 2, 2)), lg)
222 mg = MappedGrid(lg, x̄, J)
134 @test boundary_identifiers(mg) == boundary_identifiers(lg) 223 @test boundary_identifiers(mg) == boundary_identifiers(lg)
135 end 224 end
136 225
137 @testset "boundary_indices" begin 226 @testset "boundary_indices" begin
138 @test boundary_indices(mg, CartesianBoundary{1,Lower}()) == boundary_indices(lg,CartesianBoundary{1,Lower}()) 227 lg = equidistant_grid((0,0), (1,1), 11, 15)
139 @test boundary_indices(mg, CartesianBoundary{2,Lower}()) == boundary_indices(lg,CartesianBoundary{2,Lower}()) 228 x̄ = map(ξ̄ -> 2ξ̄, lg)
140 @test boundary_indices(mg, CartesianBoundary{1,Upper}()) == boundary_indices(lg,CartesianBoundary{1,Upper}()) 229 J = map(ξ̄ -> @SArray(fill(2., 2, 2)), lg)
230 mg = MappedGrid(lg, x̄, J)
231
232 @test boundary_indices(mg, CartesianBoundary{1,LowerBoundary}()) == boundary_indices(lg,CartesianBoundary{1,LowerBoundary}())
233 @test boundary_indices(mg, CartesianBoundary{2,LowerBoundary}()) == boundary_indices(lg,CartesianBoundary{2,LowerBoundary}())
234 @test boundary_indices(mg, CartesianBoundary{1,UpperBoundary}()) == boundary_indices(lg,CartesianBoundary{1,UpperBoundary}())
141 end 235 end
142 236
143 @testset "boundary_grid" begin 237 @testset "boundary_grid" begin
144 x̄, J = _partially_curved_mapping() 238 x̄, J = _partially_curved_mapping()
145 mg = mapped_grid(x̄, J, 10, 11) 239 mg = mapped_grid(x̄, J, 10, 11)
150 J2((ξ, η)) = @SMatrix[ 244 J2((ξ, η)) = @SMatrix[
151 0; 245 0;
152 1+ξ*(ξ-1); 246 1+ξ*(ξ-1);
153 ] 247 ]
154 248
155 function test_boundary_grid(mg, bId, Jb) 249 function expected_bg(mg, bId, Jb)
156 bg = boundary_grid(mg, bId)
157
158 lg = logicalgrid(mg) 250 lg = logicalgrid(mg)
159 expected_bg = MappedGrid( 251 return MappedGrid(
160 boundary_grid(lg, bId), 252 boundary_grid(lg, bId),
161 map(x̄, boundary_grid(lg, bId)), 253 map(x̄, boundary_grid(lg, bId)),
162 map(Jb, boundary_grid(lg, bId)), 254 map(Jb, boundary_grid(lg, bId)),
163 ) 255 )
164 256 end
165 @testset let bId=bId, bg=bg, expected_bg=expected_bg 257
166 @test collect(bg) == collect(expected_bg) 258 let bid = TensorGridBoundary{1, LowerBoundary}()
167 @test logicalgrid(bg) == logicalgrid(expected_bg) 259 @test boundary_grid(mg, bid) == expected_bg(mg, bid, J2)
168 @test jacobian(bg) == jacobian(expected_bg) 260 end
169 # TODO: Implement equality of a curvilinear grid and simlify the above 261
170 end 262 let bid = TensorGridBoundary{1, UpperBoundary}()
171 end 263 @test boundary_grid(mg, bid) == expected_bg(mg, bid, J2)
172 264 end
173 @testset test_boundary_grid(mg, TensorGridBoundary{1, Lower}(), J2) 265
174 @testset test_boundary_grid(mg, TensorGridBoundary{1, Upper}(), J2) 266 let bid = TensorGridBoundary{2, LowerBoundary}()
175 @testset test_boundary_grid(mg, TensorGridBoundary{2, Lower}(), J1) 267 @test boundary_grid(mg, bid) == expected_bg(mg, bid, J1)
176 @testset test_boundary_grid(mg, TensorGridBoundary{2, Upper}(), J1) 268 end
269
270 let bid = TensorGridBoundary{2, UpperBoundary}()
271 @test boundary_grid(mg, bid) == expected_bg(mg, bid, J1)
272 end
177 end 273 end
178 end 274 end
179 275
180 @testset "mapped_grid" begin 276 @testset "mapped_grid" begin
181 x̄, J = _partially_curved_mapping() 277 x̄, J = _partially_curved_mapping()
183 @test mg isa MappedGrid{SVector{2,Float64}, 2} 279 @test mg isa MappedGrid{SVector{2,Float64}, 2}
184 280
185 lg = equidistant_grid((0,0), (1,1), 10, 11) 281 lg = equidistant_grid((0,0), (1,1), 10, 11)
186 @test logicalgrid(mg) == lg 282 @test logicalgrid(mg) == lg
187 @test collect(mg) == map(x̄, lg) 283 @test collect(mg) == map(x̄, lg)
284
285 @test mapped_grid(lg, x̄, J) == mg
188 end 286 end
189 287
190 @testset "jacobian_determinant" begin 288 @testset "jacobian_determinant" begin
191 @test_broken false 289 x̄((ξ, η)) = @SVector[ξ*η, ξ + η^2]
290 J((ξ, η)) = @SMatrix[
291 η ξ;
292 1 2η;
293 ]
294
295 g = mapped_grid(x̄, J, 10, 11)
296 J = map(logicalgrid(g)) do (ξ,η)
297 2η^2 - ξ
298 end
299 @test jacobian_determinant(g) ≈ J
300
301
302 lg = equidistant_grid((0,0), (1,1), 11, 21)
303 x̄ = map(ξ̄ -> @SVector[ξ̄[1],ξ̄[2], ξ̄[1] + ξ̄[2]], lg)
304 J = map(ξ̄ -> @SMatrix[1 0; 0 1; 1 1], lg)
305 mg = MappedGrid(lg, x̄, J)
306
307 @test_broken jacobian(mg) isa AbstractArray{2,Float64}
192 end 308 end
193 309
194 @testset "metric_tensor" begin 310 @testset "metric_tensor" begin
195 @test_broken false 311 x̄((ξ, η)) = @SVector[ξ*η, ξ + η^2]
312 J((ξ, η)) = @SMatrix[
313 η ξ;
314 1 2η;
315 ]
316
317 g = mapped_grid(x̄, J, 10, 11)
318 G = map(logicalgrid(g)) do (ξ,η)
319 @SMatrix[
320 1+η^2 ξ*η+2η;
321 ξ*η+2η ξ^2 + 4η^2;
322 ]
323 end
324 @test metric_tensor(g) ≈ G
196 end 325 end
197 326
198 @testset "metric_tensor_inverse" begin 327 @testset "metric_tensor_inverse" begin
199 @test_broken false 328 x̄((ξ, η)) = @SVector[ξ + ξ^2/2, η + η^2 + ξ^2/2]
329 J((ξ, η)) = @SMatrix[
330 1+ξ 0;
331 ξ 1+η;
332 ]
333
334 g = mapped_grid(x̄, J, 10, 11)
335 G⁻¹ = map(logicalgrid(g)) do (ξ,η)
336 @SMatrix[
337 (1+η)^2 -ξ*(1+η);
338 -ξ*(1+η) (1+ξ)^2+ξ^2;
339 ]/(((1+ξ)^2+ξ^2)*(1+η)^2 - ξ^2*(1+η)^2)
340
341 end
342
343 @test metric_tensor_inverse(g) ≈ G⁻¹
200 end 344 end
201 345
202 @testset "min_spacing" begin 346 @testset "min_spacing" begin
203 let g = mapped_grid(identity, x->@SMatrix[1], 11) 347 let g = mapped_grid(identity, x->@SMatrix[1], 11)
204 @test min_spacing(g) ≈ 0.1 348 @test min_spacing(g) ≈ 0.1
235 end 379 end
236 380
237 @testset "normal" begin 381 @testset "normal" begin
238 g = mapped_grid(_partially_curved_mapping()...,10, 11) 382 g = mapped_grid(_partially_curved_mapping()...,10, 11)
239 383
240 @test normal(g, CartesianBoundary{1,Lower}()) == fill(@SVector[-1,0], 11) 384 @test normal(g, CartesianBoundary{1,LowerBoundary}()) == fill(@SVector[-1,0], 11)
241 @test normal(g, CartesianBoundary{1,Upper}()) == fill(@SVector[1,0], 11) 385 @test normal(g, CartesianBoundary{1,UpperBoundary}()) == fill(@SVector[1,0], 11)
242 @test normal(g, CartesianBoundary{2,Lower}()) == fill(@SVector[0,-1], 10) 386 @test normal(g, CartesianBoundary{2,LowerBoundary}()) == fill(@SVector[0,-1], 10)
243 @test normal(g, CartesianBoundary{2,Upper}()) ≈ map(boundary_grid(g,CartesianBoundary{2,Upper}())|>logicalgrid) do ξ̄ 387 @test normal(g, CartesianBoundary{2,UpperBoundary}()) ≈ map(boundary_grid(g,CartesianBoundary{2,UpperBoundary}())|>logicalgrid) do ξ̄
244 α = 1-2ξ̄[1] 388 α = 1-2ξ̄[1]
245 @SVector[α,1]/√(α^2 + 1) 389 @SVector[α,1]/√(α^2 + 1)
246 end 390 end
247 391
248 g = mapped_grid(_fully_curved_mapping()...,5,4) 392 g = mapped_grid(_fully_curved_mapping()...,5,4)
249 393
250 unit(v) = v/norm(v) 394 unit(v) = v/norm(v)
251 @testset let bId = CartesianBoundary{1,Lower}() 395 @testset let bId = CartesianBoundary{1,LowerBoundary}()
252 lbg = boundary_grid(logicalgrid(g), bId) 396 lbg = boundary_grid(logicalgrid(g), bId)
253 @test normal(g, bId) ≈ map(lbg) do (ξ, η) 397 @test normal(g, bId) ≈ map(lbg) do (ξ, η)
254 -unit(@SVector[1/2, η/3-1/6]) 398 -unit(@SVector[1/2, η/3-1/6])
255 end 399 end
256 end 400 end
257 401
258 @testset let bId = CartesianBoundary{1,Upper}() 402 @testset let bId = CartesianBoundary{1,UpperBoundary}()
259 lbg = boundary_grid(logicalgrid(g), bId) 403 lbg = boundary_grid(logicalgrid(g), bId)
260 @test normal(g, bId) ≈ map(lbg) do (ξ, η) 404 @test normal(g, bId) ≈ map(lbg) do (ξ, η)
261 unit(@SVector[7/2, 2η-1]/(5 + 3η + 2η^2)) 405 unit(@SVector[7/2, 2η-1]/(5 + 3η + 2η^2))
262 end 406 end
263 end 407 end
264 408
265 @testset let bId = CartesianBoundary{2,Lower}() 409 @testset let bId = CartesianBoundary{2,LowerBoundary}()
266 lbg = boundary_grid(logicalgrid(g), bId) 410 lbg = boundary_grid(logicalgrid(g), bId)
267 @test normal(g, bId) ≈ map(lbg) do (ξ, η) 411 @test normal(g, bId) ≈ map(lbg) do (ξ, η)
268 -unit(@SVector[-2ξ, 2]/(6 + ξ^2 - 2ξ)) 412 -unit(@SVector[-2ξ, 2]/(6 + ξ^2 - 2ξ))
269 end 413 end
270 end 414 end
271 415
272 @testset let bId = CartesianBoundary{2,Upper}() 416 @testset let bId = CartesianBoundary{2,UpperBoundary}()
273 lbg = boundary_grid(logicalgrid(g), bId) 417 lbg = boundary_grid(logicalgrid(g), bId)
274 @test normal(g, bId) ≈ map(lbg) do (ξ, η) 418 @test normal(g, bId) ≈ map(lbg) do (ξ, η)
275 unit(@SVector[-3ξ, 2]/(6 + ξ^2 + 3ξ)) 419 unit(@SVector[-3ξ, 2]/(6 + ξ^2 + 3ξ))
276 end 420 end
277 end 421 end