@@ -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
@@ -170,17 +170,19 @@ end
170170 N = 4 ; a = zeros (N,N,2 ) |> f
171171 WaterLily. accelerate! (a,1 ,nothing ,())
172172 @test all (a .== 0 )
173- WaterLily. accelerate! (a,1 ,(i,x,t)-> i== 1 ? t : 2 * t,())
173+ WaterLily. accelerate! (a,1. ,(i,x,t)-> i== 1 ? t : 2 * t,())
174174 @test all (a[:,:,1 ] .== 1 ) && all (a[:,:,2 ] .== 2 )
175- WaterLily. accelerate! (a,1 ,nothing ,(i,x,t) -> i== 1 ? - t : - 2 * t)
175+ WaterLily. accelerate! (a,1. ,nothing ,(i,x,t) -> i== 1 ? - t : - 2 * t)
176176 @test all (a[:,:,1 ] .== 0 ) && all (a[:,:,2 ] .== 0 )
177- WaterLily. accelerate! (a,1 ,(i,x,t) -> i== 1 ? t : 2 * t,(i,x,t) -> i== 1 ? - t : - 2 * t)
177+ WaterLily. accelerate! (a,1. ,(i,x,t) -> i== 1 ? t : 2 * t,(i,x,t) -> i== 1 ? - t : - 2 * t)
178178 @test all (a[:,:,1 ] .== 0 ) && all (a[:,:,2 ] .== 0 )
179179 # check applying body force (changes in x but not t)
180180 b = zeros (N,N,2 ) |> f
181- WaterLily. accelerate! (b,0 ,(i,x,t)-> 1 ,nothing )
181+ WaterLily. accelerate! (b,0. ,(i,x,t)-> 1 ,nothing )
182182 @test all (b .== 1 )
183- a .= 0 # reset and accelerate using a non-uniform velocity field
183+ WaterLily. accelerate! (b,1. ,(i,x,t)-> 0 ,(i,x,t)-> t)
184+ @test all (b .== 2 )
185+ a .= 0 ; b .= 1 # reset and accelerate using a non-uniform velocity field
184186 WaterLily. accelerate! (a,0. ,nothing ,(i,x,t)-> t* (x[i]+ 1.0 ))
185187 WaterLily. accelerate! (b,0 ,(i,x,t)-> x[i],nothing )
186188 @test all (b .== a)
0 commit comments