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