MFEM  v4.2.0 Finite element discretization library
ex14p.cpp
Go to the documentation of this file.
1 // MFEM Example 14 - Parallel Version
2 //
3 // Compile with: make ex14p
4 //
5 // Sample runs: mpirun -np 4 ex14p -m ../data/inline-quad.mesh -o 0
6 // mpirun -np 4 ex14p -m ../data/star.mesh -o 2
7 // mpirun -np 4 ex14p -m ../data/star-mixed.mesh -o 2
8 // mpirun -np 4 ex14p -m ../data/escher.mesh -s 1
9 // mpirun -np 4 ex14p -m ../data/fichera.mesh -s 1 -k 1
10 // mpirun -np 4 ex14p -m ../data/fichera-mixed.mesh -s 1 -k 1
11 // mpirun -np 4 ex14p -m ../data/square-disc-p2.vtk -o 2
12 // mpirun -np 4 ex14p -m ../data/square-disc-p3.mesh -o 3
13 // mpirun -np 4 ex14p -m ../data/square-disc-nurbs.mesh -o 1
14 // mpirun -np 4 ex14p -m ../data/disc-nurbs.mesh -rs 4 -o 2 -s 1 -k 0
15 // mpirun -np 4 ex14p -m ../data/pipe-nurbs.mesh -o 1
16 // mpirun -np 4 ex14p -m ../data/inline-segment.mesh -rs 5
17 // mpirun -np 4 ex14p -m ../data/amr-quad.mesh -rs 3
18 // mpirun -np 4 ex14p -m ../data/amr-hex.mesh
19 //
20 // Description: This example code demonstrates the use of MFEM to define a
21 // discontinuous Galerkin (DG) finite element discretization of
22 // the Laplace problem -Delta u = 1 with homogeneous Dirichlet
23 // boundary conditions. Finite element spaces of any order,
24 // including zero on regular grids, are supported. The example
25 // highlights the use of discontinuous spaces and DG-specific face
26 // integrators.
27 //
28 // We recommend viewing examples 1 and 9 before viewing this
29 // example.
30
31 #include "mfem.hpp"
32 #include <fstream>
33 #include <iostream>
34
35 using namespace std;
36 using namespace mfem;
37
38 class CustomSolverMonitor : public IterativeSolverMonitor
39 {
40 public:
41  CustomSolverMonitor(const ParMesh *m,
42  ParGridFunction *f) :
43  pmesh(m),
44  pgf(f) {}
45
46  void MonitorSolution(int i, double norm, const Vector &x, bool final)
47  {
48  char vishost[] = "localhost";
49  int visport = 19916;
50  int num_procs, myid;
51
52  MPI_Comm_size(pmesh->GetComm(),&num_procs);
53  MPI_Comm_rank(pmesh->GetComm(),&myid);
54
55  pgf->SetFromTrueDofs(x);
56
57  socketstream sol_sock(vishost, visport);
58  sol_sock << "parallel " << num_procs << " " << myid << "\n";
59  sol_sock.precision(8);
60  sol_sock << "solution\n" << *pmesh << *pgf
61  << "window_title 'Iteration no " << i << "'"
62  << "keys rRjlc\n" << flush;
63  }
64
65 private:
66  const ParMesh *pmesh;
67  ParGridFunction *pgf;
68 };
69
70 int main(int argc, char *argv[])
71 {
72  // 1. Initialize MPI.
73  int num_procs, myid;
74  MPI_Init(&argc, &argv);
75  MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
76  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
77
78  // 2. Parse command-line options.
79  const char *mesh_file = "../data/star.mesh";
80  int ser_ref_levels = -1;
81  int par_ref_levels = 2;
82  int order = 1;
83  double sigma = -1.0;
84  double kappa = -1.0;
85  bool visualization = 1;
86
87  OptionsParser args(argc, argv);
89  "Mesh file to use.");
91  "Number of times to refine the mesh uniformly in serial,"
92  " -1 for auto.");
94  "Number of times to refine the mesh uniformly in parallel.");
96  "Finite element order (polynomial degree) >= 0.");
98  "One of the two DG penalty parameters, typically +1/-1."
99  " See the documentation of class DGDiffusionIntegrator.");
101  "One of the two DG penalty parameters, should be positive."
102  " Negative values are replaced with (order+1)^2.");
104  "--no-visualization",
105  "Enable or disable GLVis visualization.");
106  args.Parse();
107  if (!args.Good())
108  {
109  if (myid == 0)
110  {
111  args.PrintUsage(cout);
112  }
113  MPI_Finalize();
114  return 1;
115  }
116  if (kappa < 0)
117  {
118  kappa = (order+1)*(order+1);
119  }
120  if (myid == 0)
121  {
122  args.PrintOptions(cout);
123  }
124
125  // 3. Read the (serial) mesh from the given mesh file on all processors. We
126  // can handle triangular, quadrilateral, tetrahedral and hexahedral meshes
127  // with the same code. NURBS meshes are projected to second order meshes.
128  Mesh *mesh = new Mesh(mesh_file, 1, 1);
129  int dim = mesh->Dimension();
130
131  // 4. Refine the serial mesh on all processors to increase the resolution. In
132  // this example we do 'ser_ref_levels' of uniform refinement. By default,
133  // or if ser_ref_levels < 0, we choose it to be the largest number that
134  // gives a final mesh with no more than 50,000 elements.
135  {
136  if (ser_ref_levels < 0)
137  {
138  ser_ref_levels = (int)floor(log(10000./mesh->GetNE())/log(2.)/dim);
139  }
140  for (int l = 0; l < ser_ref_levels; l++)
141  {
142  mesh->UniformRefinement();
143  }
144  }
145  if (mesh->NURBSext)
146  {
147  mesh->SetCurvature(max(order, 1));
148  }
149
150  // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
151  // this mesh further in parallel to increase the resolution. Once the
152  // parallel mesh is defined, the serial mesh can be deleted.
153  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
154  delete mesh;
155  {
156  for (int l = 0; l < par_ref_levels; l++)
157  {
158  pmesh->UniformRefinement();
159  }
160  }
161
162  // 6. Define a parallel finite element space on the parallel mesh. Here we
163  // use discontinuous finite elements of the specified order >= 0.
164  FiniteElementCollection *fec = new DG_FECollection(order, dim);
165  ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
166  HYPRE_Int size = fespace->GlobalTrueVSize();
167  if (myid == 0)
168  {
169  cout << "Number of unknowns: " << size << endl;
170  }
171
172  // 7. Set up the parallel linear form b(.) which corresponds to the
173  // right-hand side of the FEM linear system.
174  ParLinearForm *b = new ParLinearForm(fespace);
175  ConstantCoefficient one(1.0);
176  ConstantCoefficient zero(0.0);
179  new DGDirichletLFIntegrator(zero, one, sigma, kappa));
180  b->Assemble();
181
182  // 8. Define the solution vector x as a parallel finite element grid function
183  // corresponding to fespace. Initialize x with initial guess of zero.
184  ParGridFunction x(fespace);
185  x = 0.0;
186
187  // 9. Set up the bilinear form a(.,.) on the finite element space
188  // corresponding to the Laplacian operator -Delta, by adding the Diffusion
189  // domain integrator and the interior and boundary DG face integrators.
190  // Note that boundary conditions are imposed weakly in the form, so there
191  // is no need for dof elimination. After serial and parallel assembly we
192  // extract the corresponding parallel matrix A.
193  ParBilinearForm *a = new ParBilinearForm(fespace);
197  a->Assemble();
198  a->Finalize();
199
200  // 10. Define the parallel (hypre) matrix and vectors representing a(.,.),
201  // b(.) and the finite element approximation.
205
206  delete a;
207  delete b;
208
209  // 11. Depending on the symmetry of A, define and apply a parallel PCG or
210  // GMRES solver for AX=B using the BoomerAMG preconditioner from hypre.
211  HypreSolver *amg = new HypreBoomerAMG(*A);
212  if (sigma == -1.0)
213  {
214  HyprePCG pcg(*A);
215  pcg.SetTol(1e-12);
216  pcg.SetMaxIter(200);
217  pcg.SetPrintLevel(2);
218  pcg.SetPreconditioner(*amg);
219  pcg.Mult(*B, *X);
220  }
221  else
222  {
223  CustomSolverMonitor monitor(pmesh, &x);
224  GMRESSolver gmres(MPI_COMM_WORLD);
225  gmres.SetAbsTol(0.0);
226  gmres.SetRelTol(1e-12);
227  gmres.SetMaxIter(200);
228  gmres.SetKDim(10);
229  gmres.SetPrintLevel(1);
230  gmres.SetOperator(*A);
231  gmres.SetPreconditioner(*amg);
232  gmres.SetMonitor(monitor);
233  gmres.Mult(*B, *X);
234  }
235  delete amg;
236
237  // 12. Extract the parallel grid function corresponding to the finite element
238  // approximation X. This is the local solution on each processor.
239  x = *X;
240
241  // 13. Save the refined mesh and the solution in parallel. This output can
242  // be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
243  {
244  ostringstream mesh_name, sol_name;
245  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
246  sol_name << "sol." << setfill('0') << setw(6) << myid;
247
248  ofstream mesh_ofs(mesh_name.str().c_str());
249  mesh_ofs.precision(8);
250  pmesh->Print(mesh_ofs);
251
252  ofstream sol_ofs(sol_name.str().c_str());
253  sol_ofs.precision(8);
254  x.Save(sol_ofs);
255  }
256
257  // 14. Send the solution by socket to a GLVis server.
258  if (visualization)
259  {
260  char vishost[] = "localhost";
261  int visport = 19916;
262  socketstream sol_sock(vishost, visport);
263  sol_sock << "parallel " << num_procs << " " << myid << "\n";
264  sol_sock.precision(8);
265  sol_sock << "solution\n" << *pmesh << x << flush;
266  }
267
268  // 15. Free the used memory.
269  delete X;
270  delete B;
271  delete A;
272
273  delete fespace;
274  delete fec;
275  delete pmesh;
276
277  MPI_Finalize();
278
279  return 0;
280 }
void SetTol(double tol)
Definition: hypre.cpp:2619
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:93
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:79
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:737
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:819
double kappa
Definition: ex24.cpp:54
Abstract parallel finite element space.
Definition: pfespace.hpp:28
int main(int argc, char *argv[])
Definition: ex1.cpp:66
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:2634
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
virtual void Mult(const Vector &b, Vector &x) const
Operator application: y=A(x).
Definition: solvers.cpp:798
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1079
Class for parallel linear form.
Definition: plinearform.hpp:26
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:70
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition: solvers.cpp:97
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[]
double b
Definition: lissajous.cpp:42
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:8382
constexpr int visport
void SetMaxIter(int max_it)
Definition: solvers.hpp:98
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Definition: mesh.cpp:4165
HYPRE_Int GlobalTrueVSize() const
Definition: pfespace.hpp:251
void Assemble(int skip_zeros=1)
Assemble the local matrix.
void SetKDim(int dim)
Set the number of iteration to perform between restarts, default is 50.
Definition: solvers.hpp:303
virtual void Print(std::ostream &out=mfem::out) const
Definition: pmesh.cpp:4169
Adds new Boundary Face Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:66
int Dimension() const
Definition: mesh.hpp:788
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:434
void SetMaxIter(int max_iter)
Definition: hypre.cpp:2624
L2_FECollection DG_FECollection
Declare an alternative name for L2_FECollection = DG_FECollection.
Definition: fe_coll.hpp:269
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:70
Abstract base class for an iterative solver monitor.
Definition: solvers.hpp:34
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:39
void SetAbsTol(double atol)
Definition: solvers.hpp:97
void SetRelTol(double rtol)
Definition: solvers.hpp:96
PCG solver in hypre.
Definition: hypre.hpp:759
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
GMRES method.
Definition: solvers.hpp:290
void SetMonitor(IterativeSolverMonitor &m)
Set the iterative solver monitor.
Definition: solvers.hpp:112
double a
Definition: lissajous.cpp:41
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:203
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
int dim
Definition: ex24.cpp:53
Adds new interior Face Integrator. Assumes ownership of bfi.
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:2639
Adds new Domain Integrator. Assumes ownership of bfi.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:304
Class for parallel bilinear form.
Abstract class for hypre&#39;s solvers and preconditioners.
Definition: hypre.hpp:700
Vector data type.
Definition: vector.hpp:51
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition: solvers.cpp:91
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:181
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre&#39;s PCG.
Definition: hypre.cpp:2662
Class for parallel meshes.
Definition: pmesh.hpp:32
void ParallelAssemble(Vector &tv)
Assemble the vector on the true dofs, i.e. P^t v.
Definition: plinearform.cpp:46