MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
ex1.cpp
Go to the documentation of this file.
1 // MFEM Example 1
2 // GINKGO Modification
3 //
4 // Compile with: make ex1
5 //
6 // Sample runs: ex1 -m ../data/square-disc.mesh
7 // ex1 -m ../data/star.mesh
8 // ex1 -m ../data/star-mixed.mesh
9 // ex1 -m ../data/escher.mesh
10 // ex1 -m ../data/fichera.mesh
11 // ex1 -m ../data/fichera-mixed.mesh
12 // ex1 -m ../data/toroid-wedge.mesh
13 // ex1 -m ../data/square-disc-p2.vtk -o 2
14 // ex1 -m ../data/square-disc-p3.mesh -o 3
15 // ex1 -m ../data/square-disc-nurbs.mesh -o -1
16 // ex1 -m ../data/star-mixed-p2.mesh -o 2
17 // ex1 -m ../data/disc-nurbs.mesh -o -1
18 // ex1 -m ../data/pipe-nurbs.mesh -o -1
19 // ex1 -m ../data/fichera-mixed-p2.mesh -o 2
20 // ex1 -m ../data/star-surf.mesh
21 // ex1 -m ../data/square-disc-surf.mesh
22 // ex1 -m ../data/inline-segment.mesh
23 // ex1 -m ../data/amr-quad.mesh
24 // ex1 -m ../data/amr-hex.mesh
25 // ex1 -m ../data/fichera-amr.mesh
26 // ex1 -m ../data/mobius-strip.mesh
27 // ex1 -m ../data/mobius-strip.mesh -o -1 -sc
28 //
29 // Device sample runs:
30 // ex1 -pa -d cuda
31 // ex1 -pa -d raja-cuda
32 // ex1 -pa -d occa-cuda
33 // ex1 -pa -d raja-omp
34 // ex1 -pa -d occa-omp
35 // ex1 -m ../data/beam-hex.mesh -pa -d cuda
36 //
37 // Description: This example code demonstrates the use of MFEM to define a
38 // simple finite element discretization of the Laplace problem
39 // -Delta u = 1 with homogeneous Dirichlet boundary conditions.
40 // Specifically, we discretize using a FE space of the specified
41 // order, or if order < 1 using an isoparametric/isogeometric
42 // space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
43 // NURBS mesh, etc.)
44 //
45 // The example highlights the use of mesh refinement, finite
46 // element grid functions, as well as linear and bilinear forms
47 // corresponding to the left-hand side and right-hand side of the
48 // discrete linear system. We also cover the explicit elimination
49 // of essential boundary conditions, static condensation, and the
50 // optional connection to the GLVis tool for visualization.
51 
52 #include "mfem.hpp"
53 #include <fstream>
54 #include <iostream>
55 
56 #ifndef MFEM_USE_GINKGO
57 #error This example requires that MFEM is built with MFEM_USE_GINKGO=YES
58 #endif
59 
60 using namespace std;
61 using namespace mfem;
62 
63 int main(int argc, char *argv[])
64 {
65  // 1. Parse command-line options.
66  const char *mesh_file = "../../data/star.mesh";
67  int order = 1;
68  bool static_cond = false;
69  bool pa = false;
70  const char *device_config = "cpu";
71  bool visualization = true;
72  int solver_config = 0;
73  int print_lvl = 1;
74 
75  OptionsParser args(argc, argv);
76  args.AddOption(&mesh_file, "-m", "--mesh",
77  "Mesh file to use.");
78  args.AddOption(&order, "-o", "--order",
79  "Finite element order (polynomial degree) or -1 for"
80  " isoparametric space.");
81  args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
82  "--no-static-condensation", "Enable static condensation.");
83  args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
84  "--no-partial-assembly", "Enable Partial Assembly.");
85  args.AddOption(&device_config, "-d", "--device",
86  "Device configuration string, see Device::Configure().");
87  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
88  "--no-visualization",
89  "Enable or disable GLVis visualization.");
90  args.AddOption(&solver_config, "-s", "--solver-config",
91  "Solver and preconditioner combination: \n\t"
92  " 0 - Ginkgo solver and Ginkgo preconditioner, \n\t"
93  " 1 - Ginkgo solver and MFEM preconditioner, \n\t"
94  " 2 - MFEM solver and Ginkgo preconditioner, \n\t"
95  " 3 - MFEM solver and MFEM preconditioner.");
96  args.AddOption(&print_lvl, "-pl", "--print-level",
97  "Print level for iterative solver (1 prints every iteration).");
98  args.Parse();
99  if (!args.Good())
100  {
101  args.PrintUsage(cout);
102  return 1;
103  }
104  args.PrintOptions(cout);
105 
106  // 2. Enable hardware devices such as GPUs, and programming models such as
107  // CUDA, OCCA, RAJA and OpenMP based on command line options.
108  Device device(device_config);
109  device.Print();
110 
111  // 3. Read the mesh from the given mesh file. We can handle triangular,
112  // quadrilateral, tetrahedral, hexahedral, surface and volume meshes with
113  // the same code.
114  Mesh *mesh = new Mesh(mesh_file, 1, 1);
115  int dim = mesh->Dimension();
116 
117  // 4. Refine the mesh to increase the resolution. In this example we do
118  // 'ref_levels' of uniform refinement. We choose 'ref_levels' to be the
119  // largest number that gives a final mesh with no more than 50,000
120  // elements.
121  {
122  int ref_levels =
123  (int)floor(log(50000./mesh->GetNE())/log(2.)/dim);
124  for (int l = 0; l < ref_levels; l++)
125  {
126  mesh->UniformRefinement();
127  }
128  }
129 
130  // 5. Define a finite element space on the mesh. Here we use continuous
131  // Lagrange finite elements of the specified order. If order < 1, we
132  // instead use an isoparametric/isogeometric space.
134  if (order > 0)
135  {
136  fec = new H1_FECollection(order, dim);
137  }
138  else if (mesh->GetNodes())
139  {
140  fec = mesh->GetNodes()->OwnFEC();
141  cout << "Using isoparametric FEs: " << fec->Name() << endl;
142  }
143  else
144  {
145  fec = new H1_FECollection(order = 1, dim);
146  }
147  FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
148  cout << "Number of finite element unknowns: "
149  << fespace->GetTrueVSize() << endl;
150 
151  // 6. Determine the list of true (i.e. conforming) essential boundary dofs.
152  // In this example, the boundary conditions are defined by marking all
153  // the boundary attributes from the mesh as essential (Dirichlet) and
154  // converting them to a list of true dofs.
155  Array<int> ess_tdof_list;
156  if (mesh->bdr_attributes.Size())
157  {
158  Array<int> ess_bdr(mesh->bdr_attributes.Max());
159  ess_bdr = 1;
160  fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
161  }
162 
163  // 7. Set up the linear form b(.) which corresponds to the right-hand side of
164  // the FEM linear system, which in this case is (1,phi_i) where phi_i are
165  // the basis functions in the finite element fespace.
166  LinearForm *b = new LinearForm(fespace);
167  ConstantCoefficient one(1.0);
169  b->Assemble();
170 
171  // 8. Define the solution vector x as a finite element grid function
172  // corresponding to fespace. Initialize x with initial guess of zero,
173  // which satisfies the boundary conditions.
174  GridFunction x(fespace);
175  x = 0.0;
176 
177  // 9. Set up the bilinear form a(.,.) on the finite element space
178  // corresponding to the Laplacian operator -Delta, by adding the Diffusion
179  // domain integrator.
180  BilinearForm *a = new BilinearForm(fespace);
181  if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
183 
184  // 10. Assemble the bilinear form and the corresponding linear system,
185  // applying any necessary transformations such as: eliminating boundary
186  // conditions, applying conforming constraints for non-conforming AMR,
187  // static condensation, etc.
188  if (static_cond) { a->EnableStaticCondensation(); }
189  a->Assemble();
190 
191  OperatorPtr A;
192  Vector B, X;
193  a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
194 
195  cout << "Size of linear system: " << A->Height() << endl;
196 
197  // 11. Solve the linear system A X = B.
198  if (!pa)
199  {
200  switch (solver_config)
201  {
202  // Solve the linear system with CG + IC from Ginkgo
203  case 0:
204  {
205  cout << "Using Ginkgo solver + preconditioner...\n";
206  Ginkgo::GinkgoExecutor exec(device);
207  Ginkgo::IcPreconditioner ginkgo_precond(exec, "paric", 30);
208  Ginkgo::CGSolver ginkgo_solver(exec, ginkgo_precond);
209  ginkgo_solver.SetPrintLevel(print_lvl);
210  ginkgo_solver.SetRelTol(1e-12);
211  ginkgo_solver.SetAbsTol(0.0);
212  ginkgo_solver.SetMaxIter(400);
213  ginkgo_solver.SetOperator(*(A.Ptr()));
214  ginkgo_solver.Mult(B, X);
215  break;
216  }
217 
218  // Solve the linear system with CG from Ginkgo + MFEM preconditioner
219  case 1:
220  {
221  cout << "Using Ginkgo solver + MFEM preconditioner...\n";
222  Ginkgo::GinkgoExecutor exec(device);
223  //Create MFEM preconditioner and wrap it for Ginkgo's use.
224  DSmoother M((SparseMatrix&)(*A));
225  Ginkgo::MFEMPreconditioner gko_M(exec, M);
226  Ginkgo::CGSolver ginkgo_solver(exec, gko_M);
227  ginkgo_solver.SetPrintLevel(print_lvl);
228  ginkgo_solver.SetRelTol(1e-12);
229  ginkgo_solver.SetAbsTol(0.0);
230  ginkgo_solver.SetMaxIter(400);
231  ginkgo_solver.SetOperator(*(A.Ptr()));
232  ginkgo_solver.Mult(B, X);
233  break;
234  }
235 
236  // Ginkgo IC preconditioner + MFEM CG solver
237  case 2:
238  {
239  cout << "Using MFEM solver + Ginkgo preconditioner...\n";
240  Ginkgo::GinkgoExecutor exec(device);
241  Ginkgo::IcPreconditioner M(exec, "paric", 30);
242  M.SetOperator(*(A.Ptr())); // Generate the preconditioner for the matrix A.
243  PCG(*A, M, B, X, print_lvl, 400, 1e-12, 0.0);
244  break;
245  }
246 
247  // MFEM solver + MFEM preconditioner
248  case 3:
249  {
250  cout << "Using MFEM solver + MFEM preconditioner...\n";
251  // Use a simple Jacobi preconditioner with PCG.
252  DSmoother M((SparseMatrix&)(*A));
253  PCG(*A, M, B, X, print_lvl, 400, 1e-12, 0.0);
254  break;
255  }
256  } // End switch on solver_config
257  }
258  // Partial assembly mode. Cannot use Ginkgo preconditioners, but can use Ginkgo
259  // solvers.
260  else
261  {
262  if (UsesTensorBasis(*fespace))
263  {
264  // Use Jacobi preconditioning in partial assembly mode.
265  OperatorJacobiSmoother M(*a, ess_tdof_list);
266  switch (solver_config)
267  {
268  // No Ginkgo preconditioners work with matrix-free; error
269  case 0:
270  {
271  cout << "Using Ginkgo solver + preconditioner...\n";
272  MFEM_ABORT("Cannot use Ginkgo preconditioner in partial assembly mode.\n"
273  " Try -s 1 to test Ginkgo solver with an MFEM preconditioner.");
274  break;
275  }
276 
277  // Use Ginkgo solver with MFEM preconditioner
278  case 1:
279  {
280  cout << "Using Ginkgo solver + MFEM preconditioner...\n";
281  Ginkgo::GinkgoExecutor exec(device);
282  // Wrap MFEM preconditioner for Ginkgo's use.
283  Ginkgo::MFEMPreconditioner gko_M(exec, M);
284  Ginkgo::CGSolver ginkgo_solver(exec, gko_M);
285  ginkgo_solver.SetPrintLevel(print_lvl);
286  ginkgo_solver.SetRelTol(1e-12);
287  ginkgo_solver.SetAbsTol(0.0);
288  ginkgo_solver.SetMaxIter(400);
289  ginkgo_solver.SetOperator(*(A.Ptr()));
290  ginkgo_solver.Mult(B, X);
291  break;
292  }
293 
294  // No Ginkgo preconditioners work with matrix-free; error
295  case 2:
296  {
297  cout << "Using MFEM solver + Ginkgo preconditioner...\n";
298  MFEM_ABORT("Cannot use Ginkgo preconditioner in partial assembly mode.\n"
299  " Try -s 1 to test Ginkgo solver with an MFEM preconditioner.");
300  break;
301  }
302 
303  // Use MFEM solver and preconditioner
304  case 3:
305  {
306  cout << "Using MFEM solver + MFEM preconditioner...\n";
307  PCG(*A, M, B, X, print_lvl, 400, 1e-12, 0.0);
308  break;
309  }
310  } // End switch on solver_config
311  }
312  else // CG with no preconditioning
313  {
314  cout << "Using MFEM solver + no preconditioner...\n";
315  CG(*A, B, X, print_lvl, 400, 1e-12, 0.0);
316  }
317  }
318 
319  // 12. Recover the solution as a finite element grid function.
320  a->RecoverFEMSolution(X, *b, x);
321 
322  // 13. Save the refined mesh and the solution. This output can be viewed later
323  // using GLVis: "glvis -m refined.mesh -g sol.gf".
324  ofstream mesh_ofs("refined.mesh");
325  mesh_ofs.precision(8);
326  mesh->Print(mesh_ofs);
327  ofstream sol_ofs("sol.gf");
328  sol_ofs.precision(8);
329  x.Save(sol_ofs);
330 
331  // 14. Send the solution by socket to a GLVis server.
332  if (visualization)
333  {
334  char vishost[] = "localhost";
335  int visport = 19916;
336  socketstream sol_sock(vishost, visport);
337  sol_sock.precision(8);
338  sol_sock << "solution\n" << *mesh << x << flush;
339  }
340 
341  // 15. Free the used memory.
342  delete a;
343  delete b;
344  delete fespace;
345  if (order > 0) { delete fec; }
346  delete mesh;
347 
348  return 0;
349 }
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:97
void SetRelTol(double rtol)
Definition: ginkgo.hpp:776
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
Data type for scaled Jacobi-type smoother of sparse matrix.
virtual void SetOperator(const Operator &op)
Definition: ginkgo.cpp:410
void Assemble(int skip_zeros=1)
Assembles the form i.e. sums over all domain/bdr integrators.
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:102
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this bilinear form and the linear form b(...
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
Recover the solution of a linear system formed with FormLinearSystem().
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Get a list of essential true dofs, ess_tdof_list, corresponding to the boundary attributes marked in ...
Definition: fespace.cpp:564
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:928
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:279
void SetAssemblyLevel(AssemblyLevel assembly_level)
Set the desired assembly level.
void SetMaxIter(int max_it)
Definition: ginkgo.hpp:800
bool UsesTensorBasis(const FiniteElementSpace &fes)
Definition: fespace.hpp:983
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3619
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
Data type sparse matrix.
Definition: sparsemat.hpp:46
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
constexpr char vishost[]
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:274
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9737
constexpr int visport
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
void CG(const Operator &A, const Vector &b, Vector &x, int print_iter, int max_num_iter, double RTOLERANCE, double ATOLERANCE)
Conjugate gradient method. (tolerances are squared)
Definition: solvers.cpp:889
void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x, int print_iter, int max_num_iter, double RTOLERANCE, double ATOLERANCE)
Preconditioned conjugate gradient method. (tolerances are squared)
Definition: solvers.cpp:904
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:566
int Dimension() const
Definition: mesh.hpp:999
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
Operator * Ptr() const
Access the underlying Operator pointer.
Definition: handle.hpp:87
void SetPrintLevel(int print_lvl)
Definition: ginkgo.hpp:609
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:39
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:270
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:88
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
FDualNumber< tbase > log(const FDualNumber< tbase > &f)
log([dual number])
Definition: fdual.hpp:515
double a
Definition: lissajous.cpp:41
A &quot;square matrix&quot; operator for the associated FE space and BLFIntegrators The sum of all the BLFInteg...
virtual const char * Name() const
Definition: fe_coll.hpp:61
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:1646
int dim
Definition: ex24.cpp:53
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
virtual void Mult(const Vector &x, Vector &y) const
Definition: ginkgo.cpp:282
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
void SetAbsTol(double atol)
Definition: ginkgo.hpp:788
Vector data type.
Definition: vector.hpp:60
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7841
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
Vector with associated FE space and LinearFormIntegrators.
Definition: linearform.hpp:23
virtual void SetOperator(const Operator &op)
Definition: ginkgo.cpp:879
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:121
int main()
void EnableStaticCondensation()
Enable the use of static condensation. For details see the description for class StaticCondensation i...
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150