MFEM  v4.3.0 Finite element discretization library
ex1p.cpp
Go to the documentation of this file.
1 // MFEM Example 1 - Parallel Version
2 // PETSc Modification
3 //
4 // Compile with: make ex1p
5 //
6 // Sample runs: mpirun -np 4 ex1p -m ../../data/amr-quad.mesh
7 // mpirun -np 4 ex1p -m ../../data/amr-quad.mesh --petscopts rc_ex1p
8 //
9 // Device sample runs:
10 // mpirun -np 4 ex1p -pa -d cuda --petscopts rc_ex1p_cuda
11 //
12 // Description: This example code demonstrates the use of MFEM to define a
13 // simple finite element discretization of the Laplace problem
14 // -Delta u = 1 with homogeneous Dirichlet boundary conditions.
15 // Specifically, we discretize using a FE space of the specified
16 // order, or if order < 1 using an isoparametric/isogeometric
18 // NURBS mesh, etc.)
19 //
20 // The example highlights the use of mesh refinement, finite
21 // element grid functions, as well as linear and bilinear forms
22 // corresponding to the left-hand side and right-hand side of the
23 // discrete linear system. We also cover the explicit elimination
24 // of essential boundary conditions, static condensation, and the
25 // optional connection to the GLVis tool for visualization.
26 // The example also shows how PETSc Krylov solvers can be used by
27 // wrapping a HypreParMatrix (or not) and a Solver, together with
28 // customization using an options file (see rc_ex1p) We also
29 // provide an example on how to visualize the iterative solution
30 // inside a PETSc solver.
31
32 #include "mfem.hpp"
33 #include <fstream>
34 #include <iostream>
35
36 #ifndef MFEM_USE_PETSC
37 #error This example requires that MFEM is built with MFEM_USE_PETSC=YES
38 #endif
39
40 using namespace std;
41 using namespace mfem;
42
43 class UserMonitor : public PetscSolverMonitor
44 {
45 private:
48
49 public:
50  UserMonitor(ParBilinearForm *a_, ParLinearForm *b_)
51  : PetscSolverMonitor(true,false), a(a_), b(b_) {}
52
53  void MonitorSolution(PetscInt it, PetscReal norm, const Vector &X)
54  {
55  // we plot the first 5 iterates
56  if (!it || it > 5) { return; }
57  ParFiniteElementSpace *fespace = a->ParFESpace();
58  ParMesh *mesh = fespace->GetParMesh();
59  ParGridFunction x(fespace);
60  a->RecoverFEMSolution(X, *b, x);
61
62  char vishost[] = "localhost";
63  int visport = 19916;
64  int num_procs, myid;
65
66  MPI_Comm_size(mesh->GetComm(),&num_procs);
67  MPI_Comm_rank(mesh->GetComm(),&myid);
68  socketstream sol_sock(vishost, visport);
69  sol_sock << "parallel " << num_procs << " " << myid << "\n";
70  sol_sock.precision(8);
71  sol_sock << "solution\n" << *mesh << x
72  << "window_title 'Iteration no " << it << "'" << flush;
73  }
74 };
75
76 int main(int argc, char *argv[])
77 {
78  // 1. Initialize MPI.
79  int num_procs, myid;
80  MPI_Init(&argc, &argv);
81  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
82  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
83
84  // 2. Parse command-line options.
85  const char *mesh_file = "../../data/star.mesh";
86  int order = 1;
87  bool static_cond = false;
88  bool pa = false;
89  bool visualization = false;
90  const char *device_config = "cpu";
91  bool use_petsc = true;
92  const char *petscrc_file = "";
93  bool petscmonitor = false;
94  bool forcewrap = false;
95  bool useh2 = false;
96
97  OptionsParser args(argc, argv);
99  "Mesh file to use.");
101  "Finite element order (polynomial degree) or -1 for"
102  " isoparametric space.");
104  "--no-static-condensation", "Enable static condensation.");
106  "--no-partial-assembly", "Enable Partial Assembly.");
108  "Device configuration string, see Device::Configure().");
110  "--no-visualization",
111  "Enable or disable GLVis visualization.");
113  "--no-petsc",
114  "Use or not PETSc to solve the linear system.");
116  "PetscOptions file to use.");
118  "-no-petscmonitor", "--no-petscmonitor",
119  "Enable or disable GLVis visualization of residual.");
121  "-noforce-wrap", "--noforce-wrap",
122  "Force matrix-free.");
124  "--no-h2",
125  "Use or not the H2 matrix solver.");
126  args.Parse();
127  if (!args.Good())
128  {
129  if (myid == 0)
130  {
131  args.PrintUsage(cout);
132  }
133  MPI_Finalize();
134  return 1;
135  }
136  if (myid == 0)
137  {
138  args.PrintOptions(cout);
139  }
140
141  // 3. Enable hardware devices such as GPUs, and programming models such as
142  // CUDA, OCCA, RAJA and OpenMP based on command line options.
143  Device device(device_config);
144  if (myid == 0) { device.Print(); }
145
146  // 3b. We initialize PETSc
147  MFEMInitializePetsc(NULL,NULL,petscrc_file,NULL);
148
149  // 4. Read the (serial) mesh from the given mesh file on all processors. We
150  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
151  // and volume meshes with the same code.
152  Mesh *mesh = new Mesh(mesh_file, 1, 1);
153  int dim = mesh->Dimension();
154
155  // 5. Refine the serial mesh on all processors to increase the resolution. In
156  // this example we do 'ref_levels' of uniform refinement. We choose
157  // 'ref_levels' to be the largest number that gives a final mesh with no
158  // more than 10,000 elements.
159  {
160  int ref_levels =
161  (int)floor(log(10000./mesh->GetNE())/log(2.)/dim);
162  for (int l = 0; l < ref_levels; l++)
163  {
164  mesh->UniformRefinement();
165  }
166  }
167
168  // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
169  // this mesh further in parallel to increase the resolution. Once the
170  // parallel mesh is defined, the serial mesh can be deleted.
171  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
172  delete mesh;
173  {
174  int par_ref_levels = 2;
175  for (int l = 0; l < par_ref_levels; l++)
176  {
177  pmesh->UniformRefinement();
178  }
179  }
180
181  // 7. Define a parallel finite element space on the parallel mesh. Here we
182  // use continuous Lagrange finite elements of the specified order. If
183  // order < 1, we instead use an isoparametric/isogeometric space.
185  bool delete_fec;
186  if (order > 0)
187  {
188  fec = new H1_FECollection(order, dim);
189  delete_fec = true;
190  }
191  else if (pmesh->GetNodes())
192  {
193  fec = pmesh->GetNodes()->OwnFEC();
194  delete_fec = false;
195  if (myid == 0)
196  {
197  cout << "Using isoparametric FEs: " << fec->Name() << endl;
198  }
199  }
200  else
201  {
202  fec = new H1_FECollection(order = 1, dim);
203  delete_fec = true;
204  }
205  ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
206  HYPRE_BigInt size = fespace->GlobalTrueVSize();
207  if (myid == 0)
208  {
209  cout << "Number of finite element unknowns: " << size << endl;
210  }
211
212  // 8. Determine the list of true (i.e. parallel conforming) essential
213  // boundary dofs. In this example, the boundary conditions are defined
214  // by marking all the boundary attributes from the mesh as essential
215  // (Dirichlet) and converting them to a list of true dofs.
216  Array<int> ess_tdof_list;
217  if (pmesh->bdr_attributes.Size())
218  {
219  Array<int> ess_bdr(pmesh->bdr_attributes.Max());
220  ess_bdr = 1;
221  fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
222  }
223
224  // 9. Set up the parallel linear form b(.) which corresponds to the
225  // right-hand side of the FEM linear system, which in this case is
226  // (1,phi_i) where phi_i are the basis functions in fespace.
227  ParLinearForm *b = new ParLinearForm(fespace);
228  ConstantCoefficient one(1.0);
230  b->Assemble();
231
232  // 10. Define the solution vector x as a parallel finite element grid function
233  // corresponding to fespace. Initialize x with initial guess of zero,
234  // which satisfies the boundary conditions.
235  ParGridFunction x(fespace);
236  x = 0.0;
237
238  // 11. Set up the parallel bilinear form a(.,.) on the finite element space
239  // corresponding to the Laplacian operator -Delta, by adding the Diffusion
240  // domain integrator.
241  ParBilinearForm *a = new ParBilinearForm(fespace);
242  if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
244
245  // 12. Assemble the parallel bilinear form and the corresponding linear
246  // system, applying any necessary transformations such as: parallel
247  // assembly, eliminating boundary conditions, applying conforming
248  // constraints for non-conforming AMR, static condensation, etc.
249  if (static_cond) { a->EnableStaticCondensation(); }
250  a->Assemble();
251
252  OperatorPtr A;
253  Vector B, X;
254  a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
255
256  // 13. Solve the linear system A X = B.
257  // If using MFEM with HYPRE
258  // * With full assembly, use the BoomerAMG preconditioner from hypre.
259  // * With partial assembly, use Jacobi smoothing, for now.
260  // If using MFEM with PETSc
261  // * With full assembly, use command line options or H2 matrix solver
262  // * With partial assembly, wrap Jacobi smoothing, for now.
263  Solver *prec = NULL;
264  if (pa)
265  {
266  if (UsesTensorBasis(*fespace))
267  {
268  prec = new OperatorJacobiSmoother(*a, ess_tdof_list);
269  }
270  }
271  else
272  {
273  prec = new HypreBoomerAMG;
274  }
275
276  if (!use_petsc)
277  {
278  CGSolver *pcg = new CGSolver(MPI_COMM_WORLD);
279  if (prec) { pcg->SetPreconditioner(*prec); }
280  pcg->SetOperator(*A);
281  pcg->SetRelTol(1e-12);
282  pcg->SetMaxIter(200);
283  pcg->SetPrintLevel(1);
284  pcg->Mult(B, X);
285  delete pcg;
286  }
287  else
288  {
289  PetscPCGSolver *pcg;
290  // If petscrc_file has been given, we convert the HypreParMatrix to a
291  // PetscParMatrix; the user can then experiment with PETSc command line
292  // options unless forcewrap is true.
293  bool wrap = forcewrap ? true : (pa ? true : !strlen(petscrc_file));
294  if (wrap)
295  {
296  pcg = new PetscPCGSolver(MPI_COMM_WORLD);
297  pcg->SetOperator(*A);
298  if (useh2)
299  {
300  delete prec;
301  prec = new PetscH2Solver(*A.Ptr(),fespace);
302  }
303  else if (!pa) // We need to pass the preconditioner constructed from the HypreParMatrix
304  {
305  delete prec;
306  HypreParMatrix *hA = A.As<HypreParMatrix>();
307  prec = new HypreBoomerAMG(*hA);
308  }
309  if (prec) { pcg->SetPreconditioner(*prec); }
310  }
311  else // Not wrapping, pass the HypreParMatrix so that users can experiment with command line
312  {
313  HypreParMatrix *hA = A.As<HypreParMatrix>();
314  pcg = new PetscPCGSolver(*hA, false);
315  if (useh2)
316  {
317  delete prec;
318  prec = new PetscH2Solver(*hA,fespace);
319  }
320  }
321  pcg->iterative_mode = true; // iterative_mode is true by default with CGSolver
322  pcg->SetRelTol(1e-12);
323  pcg->SetAbsTol(1e-12);
324  pcg->SetMaxIter(200);
325  pcg->SetPrintLevel(1);
326
327  UserMonitor mymon(a,b);
328  if (visualization && petscmonitor)
329  {
330  pcg->SetMonitor(&mymon);
331  pcg->iterative_mode = true;
332  X.Randomize();
333  }
334  pcg->Mult(B, X);
335  delete pcg;
336  }
337
338  // 14. Recover the parallel grid function corresponding to X. This is the
339  // local finite element solution on each processor.
340  a->RecoverFEMSolution(X, *b, x);
341
342  // 15. Save the refined mesh and the solution in parallel. This output can
343  // be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
344  {
345  ostringstream mesh_name, sol_name;
346  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
347  sol_name << "sol." << setfill('0') << setw(6) << myid;
348
349  ofstream mesh_ofs(mesh_name.str().c_str());
350  mesh_ofs.precision(8);
351  pmesh->Print(mesh_ofs);
352
353  ofstream sol_ofs(sol_name.str().c_str());
354  sol_ofs.precision(8);
355  x.Save(sol_ofs);
356  }
357
358  // 16. Send the solution by socket to a GLVis server.
359  if (visualization)
360  {
361  char vishost[] = "localhost";
362  int visport = 19916;
363  socketstream sol_sock(vishost, visport);
364  sol_sock << "parallel " << num_procs << " " << myid << "\n";
365  sol_sock.precision(8);
366  sol_sock << "solution\n" << *pmesh << x << flush;
367  }
368
369  // 17. Free the used memory.
370  if (delete_fec)
371  {
372  delete fec;
373  }
374  delete a;
375  delete b;
376  delete fespace;
377  delete pmesh;
378  delete prec;
379
380  // We finalize PETSc
382
383  MPI_Finalize();
384
385  return 0;
386 }
int Size() const
Return the logical size of the array.
Definition: array.hpp:134
void SetPreconditioner(Solver &precond)
Definition: petsc.cpp:3025
void SetPrintLevel(int plev)
Definition: petsc.cpp:2374
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:97
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:790
Definition: solvers.hpp:316
OpType * As() const
Return the Operator pointer statically cast to a specified OpType. Similar to the method Get()...
Definition: handle.hpp:99
ParMesh * GetParMesh() const
Definition: pfespace.hpp:267
virtual void SetOperator(const Operator &op)
Definition: petsc.cpp:2880
A coefficient that is constant across space and time.
Definition: coefficient.hpp:78
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(...
HYPRE_Int PetscInt
Definition: petsc.hpp:47
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:616
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:275
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:846
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:841
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.
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void SetRelTol(double tol)
Definition: petsc.cpp:2297
bool iterative_mode
If true, use the second argument of Mult() as an initial guess.
Definition: operator.hpp:652
void Randomize(int seed=0)
Set random values in the vector.
Definition: vector.cpp:763
bool UsesTensorBasis(const FiniteElementSpace &fes)
Definition: fespace.hpp:957
void SetMonitor(PetscSolverMonitor *ctx)
Sets user-defined monitoring routine.
Definition: petsc.cpp:2459
void SetMaxIter(int max_iter)
Definition: petsc.cpp:2347
double PetscReal
Definition: petsc.hpp:49
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1387
Class for parallel linear form.
Definition: plinearform.hpp:26
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:71
Abstract class for monitoring PETSc&#39;s solvers.
Definition: petsc.hpp:955
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:150
constexpr char vishost[]
Jacobi smoothing for a given bilinear form (no matrix necessary).
Definition: solvers.hpp:130
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9143
constexpr int visport
void SetMaxIter(int max_it)
Definition: solvers.hpp:100
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
void Assemble(int skip_zeros=1)
Assemble the local matrix.
virtual void Print(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:4382
int Dimension() const
Definition: mesh.hpp:911
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:457
Operator * Ptr() const
Access the underlying Operator pointer.
Definition: handle.hpp:82
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:204
void SetRelTol(double rtol)
Definition: solvers.hpp:98
MPI_Comm GetComm() const
Definition: pmesh.hpp:276
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
HYPRE_Int HYPRE_BigInt
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
double a
Definition: lissajous.cpp:41
void MFEMFinalizePetsc()
Definition: petsc.cpp:192
virtual const char * Name() const
Definition: fe_coll.hpp:61
void MFEMInitializePetsc()
Convenience functions to initialize/finalize PETSc.
Definition: petsc.cpp:168
int dim
Definition: ex24.cpp:53
Adds new Domain Integrator. Assumes ownership of bfi.
void SetAbsTol(double tol)
Definition: petsc.cpp:2322
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:327
Class for parallel bilinear form.
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.hpp:330
Vector data type.
Definition: vector.hpp:60
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7343
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:92
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
Base class for solvers.
Definition: operator.hpp:648
Class for parallel grid function.
Definition: pgridfunc.hpp:32
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:121
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:277
virtual void Mult(const Vector &b, Vector &x) const
Application of the solver.
Definition: petsc.cpp:3105
Class for parallel meshes.
Definition: pmesh.hpp:32
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