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