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