xref: /libCEED/julia/LibCEED.jl/test/runtests.jl (revision 04085da31fa6846cbcc143b2bb9b87a1be852ac9)
1using Test, LibCEED, LinearAlgebra, StaticArrays
2
3showstr(x) = sprint(show, MIME("text/plain"), x)
4summarystr(x) = sprint(summary, x)
5getoutput(fname) = chomp(read(joinpath(@__DIR__, "output", fname), String))
6
7mutable struct CtxData
8    io::IOBuffer
9    x::Vector{Float64}
10end
11
12const run_dev_tests = !isrelease() || ("--run-dev-tests" in ARGS)
13
14if run_dev_tests
15    include("rundevtests.jl")
16end
17
18if LibCEED.minimum_libceed_version > ceedversion() && !run_dev_tests
19    @warn "Skipping tests because of incompatible libCEED versions."
20else
21    @testset "LibCEED Release Tests" begin
22        @testset "LibCEED" begin
23            @test ceedversion() isa VersionNumber
24            @test isrelease() isa Bool
25            @test isfile(get_libceed_path())
26        end
27
28        @testset "Ceed" begin
29            res = "/cpu/self/ref/serial"
30            c = Ceed(res)
31            @test isdeterministic(c)
32            @test getresource(c) == res
33            @test !iscuda(c)
34            @test get_preferred_memtype(c) == MEM_HOST
35            @test_throws LibCEED.CeedError create_interior_qfunction(c, "")
36            @test showstr(c) == """
37                Ceed
38                  Ceed Resource: $res
39                  Preferred MemType: host"""
40        end
41
42        @testset "Context" begin
43            c = Ceed()
44            data = zeros(3)
45            ctx = Context(c, data)
46            @test showstr(ctx) == """
47                CeedQFunctionContext
48                  Context Data Size: $(sizeof(data))"""
49            @test_throws Exception set_data!(ctx, MEM_HOST, OWN_POINTER, data)
50        end
51
52        @testset "CeedVector" begin
53            n = 10
54            c = Ceed()
55            v = CeedVector(c, n)
56            @test size(v) == (n,)
57            @test length(v) == n
58            @test axes(v) == (1:n,)
59            @test ndims(v) == 1
60            @test ndims(CeedVector) == 1
61
62            v[] = 0.0
63            @test @witharray(a = v, all(a .== 0.0))
64
65            v1 = rand(n)
66            v2 = CeedVector(c, v1)
67            @test @witharray_read(a = v2, mtype = MEM_HOST, a == v1)
68            @test Vector(v2) == v1
69            v[] = v1
70            for p ∈ [1, 2, Inf]
71                @test norm(v, p) ≈ norm(v1, p)
72            end
73            @test_throws Exception norm(v, 3)
74            @test witharray_read(sum, v) == sum(v1)
75            reciprocal!(v)
76            @test @witharray(a = v, mtype = MEM_HOST, all(a .== 1.0 ./ v1))
77
78            witharray(x -> x .= 1.0, v)
79            @test @witharray(a = v, all(a .== 1.0))
80
81            @test summarystr(v) == "$n-element CeedVector"
82            @test sprint(show, v) == @witharray_read(a = v, sprint(show, a))
83            io = IOBuffer()
84            summary(io, v)
85            println(io, ":")
86            @witharray_read(a = v, Base.print_array(io, a))
87            s1 = String(take!(io))
88            @test showstr(v) == s1
89
90            setarray!(v, MEM_HOST, USE_POINTER, v1)
91            syncarray!(v, MEM_HOST)
92            @test @witharray_read(a = v, a == v1)
93            p = takearray!(v, MEM_HOST)
94            @test p == pointer(v1)
95
96            m = rand(10, 10)
97            vm = CeedVector(c, vec(m))
98            @test @witharray_read(a = vm, size = size(m), a == m)
99
100            @test CeedVectorActive()[] == LibCEED.C.CEED_VECTOR_ACTIVE[]
101            @test CeedVectorNone()[] == LibCEED.C.CEED_VECTOR_NONE[]
102
103            w1 = rand(n)
104            w2 = rand(n)
105            w3 = rand(n)
106
107            cv1 = CeedVector(c, w1)
108            cv2 = CeedVector(c, w2)
109            cv3 = CeedVector(c, w3)
110
111            alpha = rand()
112
113            scale!(cv1, alpha)
114            w1 .*= alpha
115            @test @witharray_read(a = cv1, a == w1)
116
117            pointwisemult!(cv1, cv2, cv3)
118            w1 .= w2.*w3
119            @test @witharray_read(a = cv1, a == w1)
120
121            axpy!(alpha, cv2, cv1)
122            axpy!(alpha, w2, w1)
123            @test @witharray_read(a = cv1, a ≈ w1)
124        end
125
126        @testset "Basis" begin
127            c = Ceed()
128            dim = 3
129            ncomp = 1
130            p = 4
131            q = 6
132            b1 = create_tensor_h1_lagrange_basis(c, dim, ncomp, p, q, GAUSS_LOBATTO)
133
134            @test showstr(b1) == getoutput("b1.out")
135            @test getdimension(b1) == 3
136            @test gettopology(b1) == HEX
137            @test getnumcomponents(b1) == ncomp
138            @test getnumnodes(b1) == p^dim
139            @test getnumnodes1d(b1) == p
140            @test getnumqpts(b1) == q^dim
141            @test getnumqpts1d(b1) == q
142
143            q1d, w1d = lobatto_quadrature(3, AbscissaAndWeights)
144            @test q1d ≈ [-1.0, 0.0, 1.0]
145            @test w1d ≈ [1/3, 4/3, 1/3]
146
147            q1d, w1d = gauss_quadrature(3)
148            @test q1d ≈ [-sqrt(3/5), 0.0, sqrt(3/5)]
149            @test w1d ≈ [5/9, 8/9, 5/9]
150
151            b1d = [1.0 0.0; 0.5 0.5; 0.0 1.0]
152            d1d = [-0.5 0.5; -0.5 0.5; -0.5 0.5]
153            q1d = [-1.0, 0.0, 1.0]
154            w1d = [1/3, 4/3, 1/3]
155            q, p = size(b1d)
156            d2d = zeros(2, q*q, p*p)
157            d2d[1, :, :] = kron(b1d, d1d)
158            d2d[2, :, :] = kron(d1d, b1d)
159
160            dim2 = 2
161            b2 = create_tensor_h1_basis(c, dim2, 1, p, q, b1d, d1d, q1d, w1d)
162            @test getinterp(b2) == kron(b1d, b1d)
163            @test getinterp1d(b2) == b1d
164            @test getgrad(b2) == d2d
165            @test getgrad1d(b2) == d1d
166            @test showstr(b2) == getoutput("b2.out")
167
168            b3 = create_h1_basis(c, LINE, 1, p, q, b1d, reshape(d1d, 1, q, p), q1d, w1d)
169            @test getqref(b3) == q1d
170            @test getqweights(b3) == w1d
171            @test showstr(b3) == getoutput("b3.out")
172
173            v = rand(2)
174            vq = apply(b3, v)
175            vd = apply(b3, v; emode=EVAL_GRAD)
176            @test vq ≈ b1d*v
177            @test vd ≈ d1d*v
178
179            @test BasisCollocated()[] == LibCEED.C.CEED_BASIS_COLLOCATED[]
180        end
181
182        @testset "Request" begin
183            @test RequestImmediate()[] == LibCEED.C.CEED_REQUEST_IMMEDIATE[]
184            @test RequestOrdered()[] == LibCEED.C.CEED_REQUEST_ORDERED[]
185        end
186
187        @testset "Misc" begin
188            for dim = 1:3
189                D = CeedDim(dim)
190                J = rand(dim, dim)
191                @test det(J, D) ≈ det(J)
192                J = J + J' # make symmetric
193                @test setvoigt(SMatrix{dim,dim}(J)) == setvoigt(J, D)
194                @test getvoigt(setvoigt(J, D)) == J
195                V = zeros(dim*(dim + 1)÷2)
196                setvoigt!(V, J, D)
197                @test V == setvoigt(J, D)
198                J2 = zeros(dim, dim)
199                getvoigt!(J2, V, D)
200                @test J2 == J
201            end
202        end
203
204        @testset "QFunction" begin
205            c = Ceed()
206
207            id = create_identity_qfunction(c, 1, EVAL_INTERP, EVAL_INTERP)
208            Q = 10
209            v = rand(Q)
210            v1 = CeedVector(c, v)
211            v2 = CeedVector(c, Q)
212            apply!(id, Q, [v1], [v2])
213            @test @witharray(a = v2, a == v)
214
215            @test showstr(create_interior_qfunction(c, "Poisson3DApply")) == """
216                Gallery CeedQFunction Poisson3DApply
217                  2 Input Fields:
218                    Input Field [0]:
219                      Name: "du"
220                      Size: 3
221                      EvalMode: "gradient"
222                    Input Field [1]:
223                      Name: "qdata"
224                      Size: 6
225                      EvalMode: "none"
226                  1 Output Field:
227                    Output Field [0]:
228                      Name: "dv"
229                      Size: 3
230                      EvalMode: "gradient\""""
231
232            @interior_qf id2 = (c, (a, :in, EVAL_INTERP), (b, :out, EVAL_INTERP), b.=a)
233            v2[] = 0.0
234            apply!(id2, Q, [v1], [v2])
235            @test @witharray(a = v2, a == v)
236
237            ctxdata = CtxData(IOBuffer(), rand(3))
238            ctx = Context(c, ctxdata)
239            dim = 3
240            @interior_qf qf = (
241                c,
242                dim=dim,
243                ctxdata::CtxData,
244                (a, :in, EVAL_GRAD, dim),
245                (b, :in, EVAL_NONE),
246                (c, :out, EVAL_INTERP),
247                begin
248                    c[] = b*sum(a)
249                    show(ctxdata.io, MIME("text/plain"), ctxdata.x)
250                end,
251            )
252            set_context!(qf, ctx)
253            in_sz, out_sz = LibCEED.get_field_sizes(qf)
254            @test in_sz == [dim, 1]
255            @test out_sz == [1]
256            v1 = rand(dim)
257            v2 = rand(1)
258            cv1 = CeedVector(c, v1)
259            cv2 = CeedVector(c, v2)
260            cv3 = CeedVector(c, 1)
261            apply!(qf, 1, [cv1, cv2], [cv3])
262            @test String(take!(ctxdata.io)) == showstr(ctxdata.x)
263            @test @witharray_read(v3 = cv3, v3[1] == v2[1]*sum(v1))
264
265            @test QFunctionNone()[] == LibCEED.C.CEED_QFUNCTION_NONE[]
266        end
267
268        @testset "Operator" begin
269            c = Ceed()
270            @interior_qf id = (
271                c,
272                (input, :in, EVAL_INTERP),
273                (output, :out, EVAL_INTERP),
274                begin
275                    output[] = input
276                end,
277            )
278            b = create_tensor_h1_lagrange_basis(c, 3, 1, 3, 3, GAUSS_LOBATTO)
279            n = getnumnodes(b)
280            offsets = Vector{CeedInt}(0:n-1)
281            r = create_elem_restriction(c, 1, n, 1, 1, n, offsets)
282            op = Operator(
283                c;
284                qf=id,
285                fields=[
286                    (:input, r, b, CeedVectorActive()),
287                    (:output, r, b, CeedVectorActive()),
288                ],
289            )
290            @test showstr(op) == """
291                CeedOperator
292                  2 Fields
293                  1 Input Field:
294                    Input Field [0]:
295                      Name: "input"
296                      Active vector
297                  1 Output Field:
298                    Output Field [0]:
299                      Name: "output"
300                      Active vector"""
301
302            v = rand(n)
303            v1 = CeedVector(c, v)
304            v2 = CeedVector(c, n)
305            apply!(op, v1, v2)
306            @test @witharray_read(a1 = v1, @witharray_read(a2 = v2, a1 == a2))
307            apply_add!(op, v1, v2)
308            @test @witharray_read(a1 = v1, @witharray_read(a2 = v2, a1 + a1 == a2))
309
310            diag_vector = create_lvector(r)
311            LibCEED.assemble_diagonal!(op, diag_vector)
312            @test @witharray_read(a = diag_vector, a == ones(n))
313            # TODO: change this test after bug-fix in libCEED
314            diag_vector[] = 0.0
315            LibCEED.assemble_add_diagonal!(op, diag_vector)
316            @test @witharray(a = diag_vector, a == fill(1.0, n))
317
318            comp_op = create_composite_operator(c, [op])
319            apply!(comp_op, v1, v2)
320            @test @witharray_read(a1 = v1, @witharray_read(a2 = v2, a1 == a2))
321        end
322
323        @testset "ElemRestriction" begin
324            c = Ceed()
325            n = 10
326            offsets = Vector{CeedInt}([0:n-1; n-1:2*n-2])
327            lsize = 2*n - 1
328            r = create_elem_restriction(c, 2, n, 1, lsize, lsize, offsets)
329            @test getcompstride(r) == lsize
330            @test getnumelements(r) == 2
331            @test getelementsize(r) == n
332            @test getlvectorsize(r) == lsize
333            @test getnumcomponents(r) == 1
334            @test length(create_lvector(r)) == lsize
335            @test length(create_evector(r)) == 2*n
336            lv, ev = create_vectors(r)
337            @test length(lv) == lsize
338            @test length(ev) == 2*n
339            mult = getmultiplicity(r)
340            mult2 = ones(lsize)
341            mult2[n] = 2
342            @test mult == mult2
343            rand_lv = rand(lsize)
344            rand_ev = [rand_lv[1:n]; rand_lv[n:end]]
345            @test apply(r, rand_lv) == rand_ev
346            @test apply(r, rand_ev; tmode=TRANSPOSE) == rand_lv.*mult
347            @test showstr(r) == string(
348                "CeedElemRestriction from (19, 1) to 2 elements ",
349                "with 10 nodes each and component stride 19",
350            )
351
352            strides = CeedInt[1, n, n]
353            rs = create_elem_restriction_strided(c, 1, n, 1, n, strides)
354            @test showstr(rs) == string(
355                "CeedElemRestriction from (10, 1) to 1 elements ",
356                "with 10 nodes each and strides [1, $n, $n]",
357            )
358
359            @test ElemRestrictionNone()[] == LibCEED.C.CEED_ELEMRESTRICTION_NONE[]
360        end
361    end
362end
363