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