1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "catch.hpp"
13 #include "mfem.hpp"
14 #include <fstream>
15 #include <iostream>
16 
17 using namespace mfem;
18 
19 namespace ceed_test
20 {
21 
22 #ifdef MFEM_USE_CEED
23 
24 enum class CeedCoeffType { Const, Grid, Quad, VecConst, VecGrid, VecQuad };
25 
coeff_function(const Vector & x)26 double coeff_function(const Vector &x)
27 {
28    return 1.0 + x[0]*x[0];
29 }
30 
31 // Velocity coefficient
velocity_function(const Vector & x,Vector & v)32 void velocity_function(const Vector &x, Vector &v)
33 {
34    int dim = x.Size();
35    const double w = 1.0 + x[0]*x[0];
36    switch (dim)
37    {
38       case 1: v(0) = w; break;
39       case 2: v(0) = w*sqrt(2./3.); v(1) = w*sqrt(1./3.); break;
40       case 3: v(0) = w*sqrt(3./6.); v(1) = w*sqrt(2./6.); v(2) = w*sqrt(1./6.); break;
41    }
42 }
43 
getString(AssemblyLevel assembly)44 std::string getString(AssemblyLevel assembly)
45 {
46    switch (assembly)
47    {
48       case AssemblyLevel::NONE:
49          return "NONE";
50          break;
51       case AssemblyLevel::PARTIAL:
52          return "PARTIAL";
53          break;
54       case AssemblyLevel::ELEMENT:
55          return "ELEMENT";
56          break;
57       case AssemblyLevel::FULL:
58          return "FULL";
59          break;
60       case AssemblyLevel::LEGACY:
61          return "LEGACY";
62          break;
63    }
64    MFEM_ABORT("Unknown AssemblyLevel.");
65    return "";
66 }
67 
getString(CeedCoeffType coeff_type)68 std::string getString(CeedCoeffType coeff_type)
69 {
70    switch (coeff_type)
71    {
72       case CeedCoeffType::Const:
73          return "Const";
74          break;
75       case CeedCoeffType::Grid:
76          return "Grid";
77          break;
78       case CeedCoeffType::Quad:
79          return "Quad";
80          break;
81       case CeedCoeffType::VecConst:
82          return "VecConst";
83          break;
84       case CeedCoeffType::VecGrid:
85          return "VecGrid";
86          break;
87       case CeedCoeffType::VecQuad:
88          return "VecQuad";
89          break;
90    }
91    MFEM_ABORT("Unknown CeedCoeffType.");
92    return "";
93 }
94 
95 enum class Problem { Mass,
96                      Convection,
97                      Diffusion,
98                      VectorMass,
99                      VectorDiffusion,
100                      MassDiffusion
101                    };
102 
getString(Problem pb)103 std::string getString(Problem pb)
104 {
105    switch (pb)
106    {
107       case Problem::Mass:
108          return "Mass";
109          break;
110       case Problem::Convection:
111          return "Convection";
112          break;
113       case Problem::Diffusion:
114          return "Diffusion";
115          break;
116       case Problem::VectorMass:
117          return "VectorMass";
118          break;
119       case Problem::VectorDiffusion:
120          return "VectorDiffusion";
121          break;
122       case Problem::MassDiffusion:
123          return "MassDiffusion";
124          break;
125    }
126    MFEM_ABORT("Unknown Problem.");
127    return "";
128 }
129 
130 enum class NLProblem {Convection};
131 
getString(NLProblem pb)132 std::string getString(NLProblem pb)
133 {
134    switch (pb)
135    {
136       case NLProblem::Convection:
137          return "Convection";
138          break;
139    }
140    MFEM_ABORT("Unknown Problem.");
141    return "";
142 }
143 
InitCoeff(Mesh & mesh,FiniteElementCollection & fec,const int dim,const CeedCoeffType coeff_type,GridFunction * & gf,FiniteElementSpace * & coeff_fes,Coefficient * & coeff,VectorCoefficient * & vcoeff)144 void InitCoeff(Mesh &mesh, FiniteElementCollection &fec, const int dim,
145                const CeedCoeffType coeff_type, GridFunction *&gf,
146                FiniteElementSpace *& coeff_fes,
147                Coefficient *&coeff, VectorCoefficient *&vcoeff)
148 {
149    switch (coeff_type)
150    {
151       case CeedCoeffType::Const:
152          coeff = new ConstantCoefficient(1.0);
153          break;
154       case CeedCoeffType::Grid:
155       {
156          FunctionCoefficient f_coeff(coeff_function);
157          coeff_fes = new FiniteElementSpace(&mesh, &fec);
158          gf = new GridFunction(coeff_fes);
159          gf->ProjectCoefficient(f_coeff);
160          coeff = new GridFunctionCoefficient(gf);
161          break;
162       }
163       case CeedCoeffType::Quad:
164          coeff = new FunctionCoefficient(coeff_function);
165          break;
166       case CeedCoeffType::VecConst:
167       {
168          Vector val(dim);
169          for (int i = 0; i < dim; i++)
170          {
171             val(i) = 1.0;
172          }
173          vcoeff = new VectorConstantCoefficient(val);
174          break;
175       }
176       case CeedCoeffType::VecGrid:
177       {
178          VectorFunctionCoefficient f_vcoeff(dim, velocity_function);
179          coeff_fes = new FiniteElementSpace(&mesh, &fec, dim);
180          gf = new GridFunction(coeff_fes);
181          gf->ProjectCoefficient(f_vcoeff);
182          vcoeff = new VectorGridFunctionCoefficient(gf);
183          break;
184       }
185       case CeedCoeffType::VecQuad:
186          vcoeff = new VectorFunctionCoefficient(dim, velocity_function);
187          break;
188    }
189 }
190 
test_ceed_operator(const char * input,int order,const CeedCoeffType coeff_type,const Problem pb,const AssemblyLevel assembly)191 void test_ceed_operator(const char* input, int order,
192                         const CeedCoeffType coeff_type, const Problem pb,
193                         const AssemblyLevel assembly)
194 {
195    std::string section = "assembly: " + getString(assembly) + "\n" +
196                          "coeff_type: " + getString(coeff_type) + "\n" +
197                          "pb: " + getString(pb) + "\n" +
198                          "order: " + std::to_string(order) + "\n" +
199                          "mesh: " + input;
200    INFO(section);
201    Mesh mesh(input, 1, 1);
202    mesh.EnsureNodes();
203    int dim = mesh.Dimension();
204    H1_FECollection fec(order, dim);
205 
206    // Coefficient Initialization
207    GridFunction *gf = nullptr;
208    FiniteElementSpace *coeff_fes = nullptr;
209    Coefficient *coeff = nullptr;
210    VectorCoefficient *vcoeff = nullptr;
211    InitCoeff(mesh, fec, dim, coeff_type, gf, coeff_fes, coeff, vcoeff);
212 
213    // Build the BilinearForm
214    bool vecOp = pb == Problem::VectorMass || pb == Problem::VectorDiffusion;
215    const int vdim = vecOp ? dim : 1;
216    FiniteElementSpace fes(&mesh, &fec, vdim);
217 
218    BilinearForm k_test(&fes);
219    BilinearForm k_ref(&fes);
220    switch (pb)
221    {
222       case Problem::Mass:
223          k_ref.AddDomainIntegrator(new MassIntegrator(*coeff));
224          k_test.AddDomainIntegrator(new MassIntegrator(*coeff));
225          break;
226       case Problem::Convection:
227          k_ref.AddDomainIntegrator(new ConvectionIntegrator(*vcoeff,-1));
228          k_test.AddDomainIntegrator(new ConvectionIntegrator(*vcoeff,-1));
229          break;
230       case Problem::Diffusion:
231          k_ref.AddDomainIntegrator(new DiffusionIntegrator(*coeff));
232          k_test.AddDomainIntegrator(new DiffusionIntegrator(*coeff));
233          break;
234       case Problem::VectorMass:
235          k_ref.AddDomainIntegrator(new VectorMassIntegrator(*coeff));
236          k_test.AddDomainIntegrator(new VectorMassIntegrator(*coeff));
237          break;
238       case Problem::VectorDiffusion:
239          k_ref.AddDomainIntegrator(new VectorDiffusionIntegrator(*coeff));
240          k_test.AddDomainIntegrator(new VectorDiffusionIntegrator(*coeff));
241          break;
242       case Problem::MassDiffusion:
243          k_ref.AddDomainIntegrator(new MassIntegrator(*coeff));
244          k_test.AddDomainIntegrator(new MassIntegrator(*coeff));
245          k_ref.AddDomainIntegrator(new DiffusionIntegrator(*coeff));
246          k_test.AddDomainIntegrator(new DiffusionIntegrator(*coeff));
247          break;
248    }
249 
250    k_ref.Assemble();
251    k_ref.Finalize();
252 
253    k_test.SetAssemblyLevel(assembly);
254    k_test.Assemble();
255 
256    // Compare ceed with mfem.
257    GridFunction x(&fes), y_ref(&fes), y_test(&fes);
258 
259    x.Randomize(1);
260 
261    k_ref.Mult(x,y_ref);
262    k_test.Mult(x,y_test);
263 
264    y_test -= y_ref;
265 
266    REQUIRE(y_test.Norml2() < 1.e-12);
267    delete gf;
268    delete coeff_fes;
269    delete coeff;
270    delete vcoeff;
271 }
272 
test_ceed_nloperator(const char * input,int order,const CeedCoeffType coeff_type,const NLProblem pb,const AssemblyLevel assembly)273 void test_ceed_nloperator(const char* input, int order,
274                           const CeedCoeffType coeff_type,
275                           const NLProblem pb, const AssemblyLevel assembly)
276 {
277    std::string section = "assembly: " + getString(assembly) + "\n" +
278                          "coeff_type: " + getString(coeff_type) + "\n" +
279                          "pb: " + getString(pb) + "\n" +
280                          "order: " + std::to_string(order) + "\n" +
281                          "mesh: " + input;
282    INFO(section);
283    Mesh mesh(input, 1, 1);
284    mesh.EnsureNodes();
285    int dim = mesh.Dimension();
286    H1_FECollection fec(order, dim);
287 
288    // Coefficient Initialization
289    GridFunction *gf = nullptr;
290    FiniteElementSpace *coeff_fes = nullptr;
291    Coefficient *coeff = nullptr;
292    VectorCoefficient *vcoeff = nullptr;
293    InitCoeff(mesh, fec, dim, coeff_type, gf, coeff_fes, coeff, vcoeff);
294 
295    // Build the NonlinearForm
296    bool vecOp = pb == NLProblem::Convection;
297    const int vdim = vecOp ? dim : 1;
298    FiniteElementSpace fes(&mesh, &fec, vdim);
299 
300    NonlinearForm k_test(&fes);
301    NonlinearForm k_ref(&fes);
302    switch (pb)
303    {
304       case NLProblem::Convection:
305          k_ref.AddDomainIntegrator(new VectorConvectionNLFIntegrator(*coeff));
306          k_test.AddDomainIntegrator(new VectorConvectionNLFIntegrator(*coeff));
307          break;
308    }
309 
310    k_test.SetAssemblyLevel(assembly);
311    k_test.Setup();
312    k_ref.Setup();
313 
314    // Compare ceed with mfem.
315    GridFunction x(&fes), y_ref(&fes), y_test(&fes);
316 
317    x.Randomize(1);
318 
319    k_ref.Mult(x,y_ref);
320    k_test.Mult(x,y_test);
321 
322    y_test -= y_ref;
323 
324    REQUIRE(y_test.Norml2() < 1.e-12);
325    delete gf;
326    delete coeff_fes;
327    delete coeff;
328    delete vcoeff;
329 }
330 
331 TEST_CASE("CEED mass & diffusion", "[CEED]")
332 {
333    auto assembly = GENERATE(AssemblyLevel::PARTIAL,AssemblyLevel::NONE);
334    auto coeff_type = GENERATE(CeedCoeffType::Const,CeedCoeffType::Grid,
335                               CeedCoeffType::Quad);
336    auto pb = GENERATE(Problem::Mass,Problem::Diffusion,Problem::MassDiffusion,
337                       Problem::VectorMass,Problem::VectorDiffusion);
338    auto order = GENERATE(1);
339    auto mesh = GENERATE("../../data/inline-quad.mesh","../../data/inline-hex.mesh",
340                         "../../data/periodic-square.mesh",
341                         "../../data/star-q2.mesh","../../data/fichera-q2.mesh",
342                         "../../data/amr-quad.mesh","../../data/fichera-amr.mesh");
343    test_ceed_operator(mesh, order, coeff_type, pb, assembly);
344 } // test case
345 
346 TEST_CASE("CEED convection", "[CEED],[Convection]")
347 {
348    auto assembly = GENERATE(AssemblyLevel::PARTIAL,AssemblyLevel::NONE);
349    auto coeff_type = GENERATE(CeedCoeffType::VecConst,CeedCoeffType::VecGrid,
350                               CeedCoeffType::VecQuad);
351    auto pb = GENERATE(Problem::Convection);
352    auto order = GENERATE(1);
353    auto mesh = GENERATE("../../data/inline-quad.mesh",
354                         "../../data/inline-hex.mesh",
355                         "../../data/star-q2.mesh",
356                         "../../data/fichera-q2.mesh",
357                         "../../data/amr-quad.mesh",
358                         "../../data/fichera-amr.mesh");
359    test_ceed_operator(mesh, order, coeff_type, pb, assembly);
360 } // test case
361 
362 TEST_CASE("CEED non-linear convection", "[CEED],[NLConvection]")
363 {
364    auto assembly = GENERATE(AssemblyLevel::PARTIAL,AssemblyLevel::NONE);
365    auto coeff_type = GENERATE(CeedCoeffType::Const,CeedCoeffType::Grid,
366                               CeedCoeffType::Quad);
367    auto pb = GENERATE(NLProblem::Convection);
368    auto order = GENERATE(1);
369    auto mesh = GENERATE("../../data/inline-quad.mesh",
370                         "../../data/inline-hex.mesh",
371                         "../../data/periodic-square.mesh",
372                         "../../data/star-q2.mesh",
373                         "../../data/fichera.mesh");
374    test_ceed_nloperator(mesh, order, coeff_type, pb, assembly);
375 } // test case
376 
377 #endif
378 
379 } // namespace ceed_test
380