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