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