@@ -2,130 +2,130 @@ using GPUArrays
22using ReadVTK, WriteVTK
33
44@info " Test backends: $(join (arrays," , " )) "
5- # @testset "util.jl" begin
6- # I = CartesianIndex(1,2,3,4)
7- # @test I+δ(3,I) == CartesianIndex(1,2,4,4)
8- # @test WaterLily.CI(I,5)==CartesianIndex(1,2,3,4,5)
9- # @test WaterLily.CIj(3,I,5)==CartesianIndex(1,2,5,4)
10- # @test WaterLily.CIj(2,CartesianIndex(16,16,16,3),14)==CartesianIndex(16,14,16,3)
11-
12- # @test loc(3,CartesianIndex(3,4,5)) == SVector(3,4,4.5) .- 1.5
13- # I = CartesianIndex(rand(2:10,3)...)
14- # @test loc(0,I) == SVector(I.I...) .- 1.5
15-
16- # ex,sym = :(a[I,i] = Math.add(p.b[I],func(I,q))),[]
17- # WaterLily.grab!(sym,ex)
18- # @test ex == :(a[I, i] = Math.add(b[I], func(I, q)))
19- # @test sym == [:a, :I, :i, :(p.b), :q]
20-
21- # for f ∈ arrays
22- # p = zeros(4,5) |> f
23- # apply!(x->x[1]+x[2]+3,p) # add 2×1.5 to move edge to origin
24- # @test inside(p) == CartesianIndices((2:3,2:4))
25- # @test inside(p,buff=0) == CartesianIndices(p)
26- # @test L₂(p) == 187
27-
28- # u = zeros(5,5,2) |> f
29- # apply!((i,x)->x[i],u)
30- # @test GPUArrays.@allowscalar [u[i,j,1].-(i-2) for i in 1:3, j in 1:3]==zeros(3,3)
31-
32- # Ng, D, U = (6, 6), 2, (1.0, 0.5)
33- # u = rand(Ng..., D) |> f # vector
34- # σ = rand(Ng...) |> f # scalar
35- # BC!(u, U)
36- # @test GPUArrays.@allowscalar all(u[1, :, 1] .== U[1]) && all(u[2, :, 1] .== U[1]) && all(u[end, :, 1] .== U[1]) &&
37- # all(u[3:end-1, 1, 1] .== u[3:end-1, 2, 1]) && all(u[3:end-1, end, 1] .== u[3:end-1, end-1, 1])
38- # @test GPUArrays.@allowscalar all(u[:, 1, 2] .== U[2]) && all(u[:, 2, 2] .== U[2]) && all(u[:, end, 2] .== U[2]) &&
39- # all(u[1, 3:end-1, 2] .== u[2, 3:end-1, 2]) && all(u[end, 3:end-1, 2] .== u[end-1, 3:end-1, 2])
40-
41- # GPUArrays.@allowscalar u[end,:,1] .= 3
42- # BC!(u,U,true) # save exit values
43- # @test GPUArrays.@allowscalar all(u[end, :, 1] .== 3)
44-
45- # WaterLily.exitBC!(u,u,0) # conservative exit check
46- # @test GPUArrays.@allowscalar all(u[end,2:end-1, 1] .== U[1])
47-
48- # # test BC with function
49- # Ubc(i,x,t) = i==1 ? 1.0 : 0.5
50- # v = rand(Ng..., D) |> f # vector
51- # BC!(v,Ubc,false); BC!(u,U,false) # make sure we apply the same
52- # @test all(v[1, :, 1] .≈ u[1, :, 1]) && all(v[2, :, 1] .≈ u[2, :, 1]) && all(v[end, :, 1] .≈ u[end, :, 1])
53- # @test all(v[:, 1, 2] .≈ u[:, 1, 2]) && all(v[:, 2, 2] .≈ u[:, 2, 2]) && all(v[:, end, 2] .≈ u[:, end, 2])
54- # # test exit bc
55- # GPUArrays.@allowscalar v[end,:,1] .= 3
56- # BC!(v,Ubc,true) # save exit values
57- # @test GPUArrays.@allowscalar all(v[end, :, 1] .== 3)
58-
59- # BC!(u,U,true,(2,)) # periodic in y and save exit values
60- # @test GPUArrays.@allowscalar all(u[:, 1:2, 1] .== u[:, end-1:end, 1]) && all(u[:, 1:2, 1] .== u[:,end-1:end,1])
61- # WaterLily.perBC!(σ,(1,2)) # periodic in two directions
62- # @test GPUArrays.@allowscalar all(σ[1, 2:end-1] .== σ[end-1, 2:end-1]) && all(σ[2:end-1, 1] .== σ[2:end-1, end-1])
63-
64- # u = rand(Ng..., D) |> f # vector
65- # BC!(u,U,true,(1,)) #saveexit has no effect here as x-periodic
66- # @test GPUArrays.@allowscalar all(u[1:2, :, 1] .== u[end-1:end, :, 1]) && all(u[1:2, :, 2] .== u[end-1:end, :, 2]) &&
67- # all(u[:, 1, 2] .== U[2]) && all(u[:, 2, 2] .== U[2]) && all(u[:, end, 2] .== U[2])
68-
69- # # test interpolation
70- # a = zeros(5,5,2) |> f; b = zeros(5,5) |> f
71- # apply!((i,x)->x[i]+1.5,a); apply!(x->x[1]+1.5,b) # offset for start of grid
72- # @test GPUArrays.@allowscalar all(WaterLily.interp(SVector(2.5,1),a) .≈ [2.5,1.])
73- # @test GPUArrays.@allowscalar all(WaterLily.interp(SVector(3.5,3),a) .≈ [3.5,3.])
74- # @test GPUArrays.@allowscalar WaterLily.interp(SVector(2.5,1),b) ≈ 2.5
75- # @test GPUArrays.@allowscalar WaterLily.interp(SVector(3.5,3),b) ≈ 3.5
76- # end
77- # end
78-
79- # function Poisson_setup(poisson,N::NTuple{D};f=Array,T=Float32) where D
80- # c = ones(T,N...,D) |> f; BC!(c, ntuple(zero,D))
81- # x = zeros(T,N) |> f; z = copy(x)
82- # pois = poisson(x,c,z)
83- # soln = map(I->T(I.I[1]),CartesianIndices(N)) |> f
84- # I = first(inside(x))
85- # GPUArrays.@allowscalar @. soln -= soln[I]
86- # z = mult!(pois,soln)
87- # solver!(pois)
88- # GPUArrays.@allowscalar @. x -= x[I]
89- # return L₂(x-soln)/L₂(soln),pois
90- # end
91-
92- # @testset "Poisson.jl" begin
93- # for f ∈ arrays
94- # err,pois = Poisson_setup(Poisson,(5,5);f)
95- # @test GPUArrays.@allowscalar parent(pois.D)==f(Float32[0 0 0 0 0; 0 -2 -3 -2 0; 0 -3 -4 -3 0; 0 -2 -3 -2 0; 0 0 0 0 0])
96- # @test GPUArrays.@allowscalar parent(pois.iD)≈f(Float32[0 0 0 0 0; 0 -1/2 -1/3 -1/2 0; 0 -1/3 -1/4 -1/3 0; 0 -1/2 -1/3 -1/2 0; 0 0 0 0 0])
97- # @test err < 1e-5
98- # err,pois = Poisson_setup(Poisson,(2^6+2,2^6+2);f)
99- # @test err < 1e-6
100- # @test pois.n[] < 310
101- # err,pois = Poisson_setup(Poisson,(2^4+2,2^4+2,2^4+2);f)
102- # @test err < 1e-6
103- # @test pois.n[] < 35
104- # end
105- # end
106-
107- # @testset "MultiLevelPoisson.jl" begin
108- # I = CartesianIndex(4,3,2)
109- # @test all(WaterLily.down(J)==I for J ∈ WaterLily.up(I))
110- # @test_throws AssertionError("MultiLevelPoisson requires size=a2ⁿ, where n>2") Poisson_setup(MultiLevelPoisson,(15+2,3^4+2))
111-
112- # err,pois = Poisson_setup(MultiLevelPoisson,(10,10))
113- # @test pois.levels[3].D == Float32[0 0 0 0; 0 -2 -2 0; 0 -2 -2 0; 0 0 0 0]
114- # @test err < 1e-5
115-
116- # pois.levels[1].L[5:6,:,1].=0
117- # WaterLily.update!(pois)
118- # @test pois.levels[3].D == Float32[0 0 0 0; 0 -1 -1 0; 0 -1 -1 0; 0 0 0 0]
119-
120- # for f ∈ arrays
121- # err,pois = Poisson_setup(MultiLevelPoisson,(2^6+2,2^6+2);f)
122- # @test err < 1e-6
123- # @test pois.n[] ≤ 3
124- # err,pois = Poisson_setup(MultiLevelPoisson,(2^4+2,2^4+2,2^4+2);f)
125- # @test err < 1e-6
126- # @test pois.n[] ≤ 3
127- # end
128- # end
5+ @testset " util.jl" begin
6+ I = CartesianIndex (1 ,2 ,3 ,4 )
7+ @test I+ δ (3 ,I) == CartesianIndex (1 ,2 ,4 ,4 )
8+ @test WaterLily. CI (I,5 )== CartesianIndex (1 ,2 ,3 ,4 ,5 )
9+ @test WaterLily. CIj (3 ,I,5 )== CartesianIndex (1 ,2 ,5 ,4 )
10+ @test WaterLily. CIj (2 ,CartesianIndex (16 ,16 ,16 ,3 ),14 )== CartesianIndex (16 ,14 ,16 ,3 )
11+
12+ @test loc (3 ,CartesianIndex (3 ,4 ,5 )) == SVector (3 ,4 ,4.5 ) .- 1.5
13+ I = CartesianIndex (rand (2 : 10 ,3 )... )
14+ @test loc (0 ,I) == SVector (I. I... ) .- 1.5
15+
16+ ex,sym = :(a[I,i] = Math. add (p. b[I],func (I,q))),[]
17+ WaterLily. grab! (sym,ex)
18+ @test ex == :(a[I, i] = Math. add (b[I], func (I, q)))
19+ @test sym == [:a , :I , :i , :(p. b), :q ]
20+
21+ for f ∈ arrays
22+ p = zeros (4 ,5 ) |> f
23+ apply! (x-> x[1 ]+ x[2 ]+ 3 ,p) # add 2×1.5 to move edge to origin
24+ @test inside (p) == CartesianIndices ((2 : 3 ,2 : 4 ))
25+ @test inside (p,buff= 0 ) == CartesianIndices (p)
26+ @test L₂ (p) == 187
27+
28+ u = zeros (5 ,5 ,2 ) |> f
29+ apply! ((i,x)-> x[i],u)
30+ @test GPUArrays. @allowscalar [u[i,j,1 ]. - (i- 2 ) for i in 1 : 3 , j in 1 : 3 ]== zeros (3 ,3 )
31+
32+ Ng, D, U = (6 , 6 ), 2 , (1.0 , 0.5 )
33+ u = rand (Ng... , D) |> f # vector
34+ σ = rand (Ng... ) |> f # scalar
35+ BC! (u, U)
36+ @test GPUArrays. @allowscalar all (u[1 , :, 1 ] .== U[1 ]) && all (u[2 , :, 1 ] .== U[1 ]) && all (u[end , :, 1 ] .== U[1 ]) &&
37+ all (u[3 : end - 1 , 1 , 1 ] .== u[3 : end - 1 , 2 , 1 ]) && all (u[3 : end - 1 , end , 1 ] .== u[3 : end - 1 , end - 1 , 1 ])
38+ @test GPUArrays. @allowscalar all (u[:, 1 , 2 ] .== U[2 ]) && all (u[:, 2 , 2 ] .== U[2 ]) && all (u[:, end , 2 ] .== U[2 ]) &&
39+ all (u[1 , 3 : end - 1 , 2 ] .== u[2 , 3 : end - 1 , 2 ]) && all (u[end , 3 : end - 1 , 2 ] .== u[end - 1 , 3 : end - 1 , 2 ])
40+
41+ GPUArrays. @allowscalar u[end ,:,1 ] .= 3
42+ BC! (u,U,true ) # save exit values
43+ @test GPUArrays. @allowscalar all (u[end , :, 1 ] .== 3 )
44+
45+ WaterLily. exitBC! (u,u,0 ) # conservative exit check
46+ @test GPUArrays. @allowscalar all (u[end ,2 : end - 1 , 1 ] .== U[1 ])
47+
48+ # test BC with function
49+ Ubc (i,x,t) = i== 1 ? 1.0 : 0.5
50+ v = rand (Ng... , D) |> f # vector
51+ BC! (v,Ubc,false ); BC! (u,U,false ) # make sure we apply the same
52+ @test all (v[1 , :, 1 ] .≈ u[1 , :, 1 ]) && all (v[2 , :, 1 ] .≈ u[2 , :, 1 ]) && all (v[end , :, 1 ] .≈ u[end , :, 1 ])
53+ @test all (v[:, 1 , 2 ] .≈ u[:, 1 , 2 ]) && all (v[:, 2 , 2 ] .≈ u[:, 2 , 2 ]) && all (v[:, end , 2 ] .≈ u[:, end , 2 ])
54+ # test exit bc
55+ GPUArrays. @allowscalar v[end ,:,1 ] .= 3
56+ BC! (v,Ubc,true ) # save exit values
57+ @test GPUArrays. @allowscalar all (v[end , :, 1 ] .== 3 )
58+
59+ BC! (u,U,true ,(2 ,)) # periodic in y and save exit values
60+ @test GPUArrays. @allowscalar all (u[:, 1 : 2 , 1 ] .== u[:, end - 1 : end , 1 ]) && all (u[:, 1 : 2 , 1 ] .== u[:,end - 1 : end ,1 ])
61+ WaterLily. perBC! (σ,(1 ,2 )) # periodic in two directions
62+ @test GPUArrays. @allowscalar all (σ[1 , 2 : end - 1 ] .== σ[end - 1 , 2 : end - 1 ]) && all (σ[2 : end - 1 , 1 ] .== σ[2 : end - 1 , end - 1 ])
63+
64+ u = rand (Ng... , D) |> f # vector
65+ BC! (u,U,true ,(1 ,)) # saveexit has no effect here as x-periodic
66+ @test GPUArrays. @allowscalar all (u[1 : 2 , :, 1 ] .== u[end - 1 : end , :, 1 ]) && all (u[1 : 2 , :, 2 ] .== u[end - 1 : end , :, 2 ]) &&
67+ all (u[:, 1 , 2 ] .== U[2 ]) && all (u[:, 2 , 2 ] .== U[2 ]) && all (u[:, end , 2 ] .== U[2 ])
68+
69+ # test interpolation
70+ a = zeros (5 ,5 ,2 ) |> f; b = zeros (5 ,5 ) |> f
71+ apply! ((i,x)-> x[i]+ 1.5 ,a); apply! (x-> x[1 ]+ 1.5 ,b) # offset for start of grid
72+ @test GPUArrays. @allowscalar all (WaterLily. interp (SVector (2.5 ,1 ),a) .≈ [2.5 ,1. ])
73+ @test GPUArrays. @allowscalar all (WaterLily. interp (SVector (3.5 ,3 ),a) .≈ [3.5 ,3. ])
74+ @test GPUArrays. @allowscalar WaterLily. interp (SVector (2.5 ,1 ),b) ≈ 2.5
75+ @test GPUArrays. @allowscalar WaterLily. interp (SVector (3.5 ,3 ),b) ≈ 3.5
76+ end
77+ end
78+
79+ function Poisson_setup (poisson,N:: NTuple{D} ;f= Array,T= Float32) where D
80+ c = ones (T,N... ,D) |> f; BC! (c, ntuple (zero,D))
81+ x = zeros (T,N) |> f; z = copy (x)
82+ pois = poisson (x,c,z)
83+ soln = map (I-> T (I. I[1 ]),CartesianIndices (N)) |> f
84+ I = first (inside (x))
85+ GPUArrays. @allowscalar @. soln -= soln[I]
86+ z = mult! (pois,soln)
87+ solver! (pois)
88+ GPUArrays. @allowscalar @. x -= x[I]
89+ return L₂ (x- soln)/ L₂ (soln),pois
90+ end
91+
92+ @testset " Poisson.jl" begin
93+ for f ∈ arrays
94+ err,pois = Poisson_setup (Poisson,(5 ,5 );f)
95+ @test GPUArrays. @allowscalar parent (pois. D)== f (Float32[0 0 0 0 0 ; 0 - 2 - 3 - 2 0 ; 0 - 3 - 4 - 3 0 ; 0 - 2 - 3 - 2 0 ; 0 0 0 0 0 ])
96+ @test GPUArrays. @allowscalar parent (pois. iD)≈ f (Float32[0 0 0 0 0 ; 0 - 1 / 2 - 1 / 3 - 1 / 2 0 ; 0 - 1 / 3 - 1 / 4 - 1 / 3 0 ; 0 - 1 / 2 - 1 / 3 - 1 / 2 0 ; 0 0 0 0 0 ])
97+ @test err < 1e-5
98+ err,pois = Poisson_setup (Poisson,(2 ^ 6 + 2 ,2 ^ 6 + 2 );f)
99+ @test err < 1e-6
100+ @test pois. n[] < 310
101+ err,pois = Poisson_setup (Poisson,(2 ^ 4 + 2 ,2 ^ 4 + 2 ,2 ^ 4 + 2 );f)
102+ @test err < 1e-6
103+ @test pois. n[] < 35
104+ end
105+ end
106+
107+ @testset " MultiLevelPoisson.jl" begin
108+ I = CartesianIndex (4 ,3 ,2 )
109+ @test all (WaterLily. down (J)== I for J ∈ WaterLily. up (I))
110+ @test_throws AssertionError (" MultiLevelPoisson requires size=a2ⁿ, where n>2" ) Poisson_setup (MultiLevelPoisson,(15 + 2 ,3 ^ 4 + 2 ))
111+
112+ err,pois = Poisson_setup (MultiLevelPoisson,(10 ,10 ))
113+ @test pois. levels[3 ]. D == Float32[0 0 0 0 ; 0 - 2 - 2 0 ; 0 - 2 - 2 0 ; 0 0 0 0 ]
114+ @test err < 1e-5
115+
116+ pois. levels[1 ]. L[5 : 6 ,:,1 ]. = 0
117+ WaterLily. update! (pois)
118+ @test pois. levels[3 ]. D == Float32[0 0 0 0 ; 0 - 1 - 1 0 ; 0 - 1 - 1 0 ; 0 0 0 0 ]
119+
120+ for f ∈ arrays
121+ err,pois = Poisson_setup (MultiLevelPoisson,(2 ^ 6 + 2 ,2 ^ 6 + 2 );f)
122+ @test err < 1e-6
123+ @test pois. n[] ≤ 3
124+ err,pois = Poisson_setup (MultiLevelPoisson,(2 ^ 4 + 2 ,2 ^ 4 + 2 ,2 ^ 4 + 2 );f)
125+ @test err < 1e-6
126+ @test pois. n[] ≤ 3
127+ end
128+ end
129129
130130@testset " Flow.jl" begin
131131 # test than vanLeer behaves correctly
0 commit comments