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 "unit_tests.hpp"
13 #include "mfem.hpp"
14 
15 using namespace mfem;
16 
17 #ifdef MFEM_USE_SUITESPARSE
18 #define DIRECT_SOLVE_SERIAL
19 #endif
20 #ifdef MFEM_USE_MUMPS
21 #define DIRECT_SOLVE_PARALLEL
22 #endif
23 #ifdef MFEM_USE_SUPERLU
24 #define DIRECT_SOLVE_PARALLEL
25 #endif
26 
27 #if defined(DIRECT_SOLVE_SERIAL) || defined(DIRECT_SOLVE_PARALLEL)
28 
29 int dim;
uexact(const Vector & x)30 double uexact(const Vector& x)
31 {
32    double u;
33    switch (dim)
34    {
35       case 1:
36          u  = 3.0 + 2.0 * x(0) - 0.5 * x(0) * x(0);
37          break;
38       case 2:
39          u = 1.0 + 0.2 * x(0) - 0.9 * x(0) * x(1) + x(1) * x(1) * x(0);
40          break;
41       default:
42          u = x(2) * x(2) * x(2) - 5.0 * x(0) * x(0) * x(1) * x(2);
43          break;
44    }
45    return u;
46 }
47 
gradexact(const Vector & x,Vector & grad)48 void gradexact(const Vector& x, Vector & grad)
49 {
50    grad.SetSize(dim);
51    switch (dim)
52    {
53       case 1:
54          grad[0] = 2.0 - x(0);
55          break;
56       case 2:
57          grad[0] = 0.2 - 0.9 * x(1) + x(1) * x (1);
58          grad[1] = - 0.9 * x(0) + 2.0 * x(0) * x(1);
59          break;
60       default:
61          grad[0] = -10.0 * x(0) * x(1) * x(2);
62          grad[1] = - 5.0 * x(0) * x(0) * x(2);
63          grad[2] = 3.0 * x(2) * x(2) - 5.0 * x(0) * x(0) * x(1);
64          break;
65    }
66 }
67 
d2uexact(const Vector & x)68 double d2uexact(const Vector& x) // returns \Delta u
69 {
70    double d2u;
71    switch (dim)
72    {
73       case 1:
74          d2u  = -1.0;
75          break;
76       case 2:
77          d2u = 2.0 * x(0);
78          break;
79       default:
80          d2u = -10.0 * x(1) * x(2) + 6.0 * x(2);
81          break;
82    }
83    return d2u;
84 }
85 
fexact(const Vector & x)86 double fexact(const Vector& x) // returns -\Delta u
87 {
88    double d2u = d2uexact(x);
89    return -d2u;
90 }
91 
92 #endif
93 
94 #ifdef DIRECT_SOLVE_SERIAL
95 
96 TEST_CASE("direct-serial","[CUDA]")
97 {
98    const int ne = 2;
99    for (dim = 1; dim < 4; ++dim)
100    {
101       Mesh mesh;
102       if (dim == 1)
103       {
104          mesh = Mesh::MakeCartesian1D(ne,  1.0);
105       }
106       else if (dim == 2)
107       {
108          mesh = Mesh::MakeCartesian2D(
109                    ne, ne, Element::QUADRILATERAL, 1, 1.0, 1.0);
110       }
111       else
112       {
113          mesh = Mesh::MakeCartesian3D(
114                    ne, ne, ne, Element::HEXAHEDRON, 1.0, 1.0, 1.0);
115       }
116       int order = 3;
117       FiniteElementCollection* fec = new H1_FECollection(order, dim);
118       FiniteElementSpace fespace(&mesh, fec);
119       Array<int> ess_tdof_list;
120       Array<int> ess_bdr(mesh.bdr_attributes.Max());
121       ess_bdr = 1;
122       fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
123 
124       FunctionCoefficient f(fexact);
125       LinearForm b(&fespace);
126       b.AddDomainIntegrator(new DomainLFIntegrator(f));
127       b.Assemble();
128 
129       BilinearForm a(&fespace);
130       ConstantCoefficient one(1.0);
131       a.AddDomainIntegrator(new DiffusionIntegrator(one));
132       a.Assemble();
133 
134       GridFunction x(&fespace);
135       FunctionCoefficient uex(uexact);
136       x = 0.0;
137       x.ProjectBdrCoefficient(uex,ess_bdr);
138 
139       OperatorPtr A;
140       Vector B, X;
141       a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
142 
143       UMFPackSolver umf_solver;
144       umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
145       umf_solver.SetOperator(*A);
146       umf_solver.Mult(B, X);
147 
148       Vector Y(X.Size());
149       A->Mult(X,Y);
150       Y-=B;
151       REQUIRE(Y.Norml2() < 1.e-12);
152 
153       a.RecoverFEMSolution(X, b, x);
154       VectorFunctionCoefficient grad(dim,gradexact);
155       double err = x.ComputeH1Error(&uex,&grad);
156       REQUIRE(err < 1.e-12);
157       delete fec;
158    }
159 }
160 
161 #endif
162 
163 #ifdef DIRECT_SOLVE_PARALLEL
164 
165 TEST_CASE("direct-parallel", "[Parallel], [CUDA]")
166 {
167    int rank;
168    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
169    const int ne = 2;
170    for (dim = 1; dim < 4; ++dim)
171    {
172       Mesh mesh;
173       if (dim == 1)
174       {
175          mesh = Mesh::MakeCartesian1D(ne,  1.0);
176       }
177       else if (dim == 2)
178       {
179          mesh = Mesh::MakeCartesian2D(
180                    ne, ne, Element::QUADRILATERAL, 1, 1.0, 1.0);
181       }
182       else
183       {
184          mesh = Mesh::MakeCartesian3D(
185                    ne, ne, ne, Element::HEXAHEDRON, 1.0, 1.0, 1.0);
186       }
187 
188       ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, mesh);
189       mesh.Clear();
190       int order = 3;
191       FiniteElementCollection* fec = new H1_FECollection(order, dim);
192       ParFiniteElementSpace fespace(pmesh, fec);
193       Array<int> ess_tdof_list;
194       Array<int> ess_bdr;
195       if (pmesh->bdr_attributes.Size())
196       {
197          ess_bdr.SetSize(pmesh->bdr_attributes.Max());
198          ess_bdr = 1;
199          fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
200       }
201       FunctionCoefficient f(fexact);
202       ParLinearForm b(&fespace);
203       b.AddDomainIntegrator(new DomainLFIntegrator(f));
204       b.Assemble();
205 
206       ParBilinearForm a(&fespace);
207       ConstantCoefficient one(1.0);
208       a.AddDomainIntegrator(new DiffusionIntegrator(one));
209       a.Assemble();
210 
211       ParGridFunction x(&fespace);
212       FunctionCoefficient uex(uexact);
213       x = 0.0;
214       x.ProjectBdrCoefficient(uex,ess_bdr);
215 
216       OperatorPtr A;
217       Vector B, X;
218       a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
219 
220 #ifdef MFEM_USE_MUMPS
221       {
222          MUMPSSolver mumps;
223          mumps.SetPrintLevel(0);
224          mumps.SetOperator(*A.As<HypreParMatrix>());
225          mumps.Mult(B,X);
226          Vector Y(X.Size());
227          A->Mult(X,Y);
228          Y-=B;
229          REQUIRE(Y.Norml2() < 1.e-12);
230 
231          a.RecoverFEMSolution(X, b, x);
232          VectorFunctionCoefficient grad(dim,gradexact);
233          double err = x.ComputeH1Error(&uex,&grad);
234          REQUIRE(err < 1.e-12);
235       }
236 #endif
237 #ifdef MFEM_USE_SUPERLU
238       // Transform to monolithic HypreParMatrix
239       {
240          SuperLURowLocMatrix SA(*A.As<HypreParMatrix>());
241          SuperLUSolver superlu(MPI_COMM_WORLD);
242          superlu.SetPrintStatistics(false);
243          superlu.SetSymmetricPattern(false);
244          superlu.SetColumnPermutation(superlu::METIS_AT_PLUS_A);
245          superlu.SetOperator(SA);
246          superlu.Mult(B, X);
247          Vector Y(X.Size());
248          A->Mult(X,Y);
249          Y-=B;
250          REQUIRE(Y.Norml2() < 1.e-12);
251          a.RecoverFEMSolution(X, b, x);
252          VectorFunctionCoefficient grad(dim,gradexact);
253          double err = x.ComputeH1Error(&uex,&grad);
254          REQUIRE(err < 1.e-12);
255       }
256 #endif
257       delete fec;
258       delete pmesh;
259    }
260 }
261 
262 #endif
263