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 #ifdef MFEM_USE_MPI
20
21 namespace sparse_matrix_test
22 {
23
coeff_function(const Vector & x)24 double coeff_function(const Vector &x)
25 {
26 return 1.0 + x[0]*x[0];
27 }
28
velocity_function(const Vector & x,Vector & v)29 void velocity_function(const Vector &x, Vector &v)
30 {
31 int dim = x.Size();
32 switch (dim)
33 {
34 case 1: v(0) = 1.0; break;
35 case 2: v(0) = x(1); v(1) = -x(0); break;
36 case 3: v(0) = x(1); v(1) = -x(0); v(2) = x(0); break;
37 }
38 }
39
40 enum class Coeff {Const, Grid, Quad};
41
getString(Coeff coeff_type)42 static std::string getString(Coeff coeff_type)
43 {
44 switch (coeff_type)
45 {
46 case Coeff::Const:
47 return "Const";
48 break;
49 case Coeff::Grid:
50 return "Grid";
51 break;
52 case Coeff::Quad:
53 return "Quad";
54 break;
55 }
56 mfem_error("Unknown CeedCoeff.");
57 return "";
58 }
59
60 enum class Problem {Mass, Convection, Diffusion};
61
getString(Problem pb)62 static std::string getString(Problem pb)
63 {
64 switch (pb)
65 {
66 case Problem::Mass:
67 return "Mass";
68 break;
69 case Problem::Convection:
70 return "Convection";
71 break;
72 case Problem::Diffusion:
73 return "Diffusion";
74 break;
75 }
76 mfem_error("Unknown Problem.");
77 return "";
78 }
79
test_sparse_matrix(const char * input,int order,const Coeff coeff_type,const Problem pb,const bool keep_nbr_block,int basis)80 void test_sparse_matrix(const char* input, int order, const Coeff coeff_type,
81 const Problem pb, const bool keep_nbr_block,
82 int basis)
83 {
84 std::string knb = keep_nbr_block ? "ON" : "OFF";
85 std::string section = "keep_nbr_block: " + knb + "\n" +
86 "coeff_type: " + getString(coeff_type) + "\n" +
87 "pb: " + getString(pb) + "\n" +
88 "order: " + std::to_string(order) + "\n" +
89 "mesh: " + input;
90 INFO(section);
91 Mesh mesh(input, 1, 1);
92 mesh.EnsureNodes();
93 if (mesh.GetNE() < 16) { mesh.UniformRefinement(); }
94 ParMesh pmesh(MPI_COMM_WORLD, mesh);
95 int dim = mesh.Dimension();
96
97 FiniteElementCollection *fec;
98 if (pb == Problem::Convection)
99 {
100 fec = new L2_FECollection(order, dim, basis);
101 }
102 else
103 {
104 fec = new H1_FECollection(order, dim, basis);
105 }
106 ParFiniteElementSpace fes(&pmesh, fec);
107
108 ParBilinearForm k_test(&fes);
109 ParBilinearForm k_ref(&fes);
110
111 ParFiniteElementSpace coeff_fes(&pmesh, fec);
112 ParGridFunction gf(&coeff_fes);
113
114 Coefficient *coeff = nullptr;
115 ConstantCoefficient rho(1.0);
116 VectorFunctionCoefficient velocity(dim, velocity_function);
117 switch (coeff_type)
118 {
119 case Coeff::Const:
120 coeff = new ConstantCoefficient(1.0);
121 break;
122 case Coeff::Grid:
123 {
124 FunctionCoefficient f_coeff(coeff_function);
125 gf.ProjectCoefficient(f_coeff);
126 coeff = new GridFunctionCoefficient(&gf);
127 break;
128 }
129 case Coeff::Quad:
130 coeff = new FunctionCoefficient(coeff_function);
131 break;
132 }
133
134 switch (pb)
135 {
136 case Problem::Mass:
137 k_ref.AddDomainIntegrator(new MassIntegrator(*coeff));
138 k_test.AddDomainIntegrator(new MassIntegrator(*coeff));
139 break;
140 case Problem::Convection:
141 k_ref.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
142 k_ref.AddInteriorFaceIntegrator(
143 new TransposeIntegrator(new DGTraceIntegrator(rho, velocity, 1.0, -0.5)));
144 k_ref.AddBdrFaceIntegrator(
145 new TransposeIntegrator(new DGTraceIntegrator(rho, velocity, 1.0, -0.5)));
146 k_test.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
147 k_test.AddInteriorFaceIntegrator(
148 new TransposeIntegrator(new DGTraceIntegrator(rho, velocity, 1.0, -0.5)));
149 k_test.AddBdrFaceIntegrator(
150 new TransposeIntegrator(new DGTraceIntegrator(rho, velocity, 1.0, -0.5)));
151 break;
152 case Problem::Diffusion:
153 k_ref.AddDomainIntegrator(new DiffusionIntegrator(*coeff));
154 k_test.AddDomainIntegrator(new DiffusionIntegrator(*coeff));
155 break;
156 }
157
158 if (keep_nbr_block) { k_ref.KeepNbrBlock(); }
159 k_ref.Assemble();
160 k_ref.Finalize();
161
162 k_test.SetAssemblyLevel(AssemblyLevel::FULL);
163 if (keep_nbr_block) { k_test.KeepNbrBlock(); }
164 k_test.Assemble();
165
166 const int sizeIn = pb == Problem::Convection ?
167 fes.GetVSize() + fes.GetFaceNbrVSize() :
168 fes.GetVSize();
169 const int sizeOut = (pb == Problem::Convection && keep_nbr_block)?
170 fes.GetVSize() + fes.GetFaceNbrVSize() :
171 fes.GetVSize();
172 const int sizeEnd = fes.GetVSize();
173 Vector x(sizeIn), y_test(sizeOut), y_ref(sizeOut);
174 x.Randomize(1);
175
176 k_test.SpMat().Mult(x,y_test);
177 k_ref.SpMat().Mult(x,y_ref);
178
179 y_test -= y_ref;
180
181 Vector result(y_test.HostReadWrite(), sizeEnd);
182
183 REQUIRE(result.Norml2() < 1.e-12);
184 delete coeff;
185 delete fec;
186 }
187
188 TEST_CASE("Sparse Matrix", "[Parallel]")
189 {
190 auto basis = GENERATE(BasisType::GaussLobatto,BasisType::Positive);
191 auto keep_nbr_block = GENERATE(false);
192 auto coeff_type = GENERATE(Coeff::Const,Coeff::Grid,Coeff::Quad);
193 auto pb = GENERATE(Problem::Mass,Problem::Convection,Problem::Diffusion);
194 auto order = GENERATE(1,2,3);
195 auto mesh = GENERATE("../../data/inline-quad.mesh",
196 "../../data/inline-hex.mesh",
197 "../../data/star-q2.mesh",
198 "../../data/fichera-q2.mesh");
199 test_sparse_matrix(mesh, order, coeff_type, pb, keep_nbr_block, basis);
200 } // test case
201
202 } // namespace sparse_matrix_test
203
204 #endif // MFEM_USE_MPI
205