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