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