MFEM  v4.6.0
Finite element discretization library
ex13p.cpp
Go to the documentation of this file.
1 // MFEM Example 13 - Parallel Version
2 //
3 // Compile with: make ex13p
4 //
5 // Sample runs: mpirun -np 4 ex13p -m ../data/star.mesh
6 // mpirun -np 4 ex13p -m ../data/square-disc.mesh -o 2 -n 4
7 // mpirun -np 4 ex13p -m ../data/beam-tet.mesh
8 // mpirun -np 4 ex13p -m ../data/beam-hex.mesh
9 // mpirun -np 4 ex13p -m ../data/escher.mesh
10 // mpirun -np 4 ex13p -m ../data/fichera.mesh
11 // mpirun -np 4 ex13p -m ../data/fichera-q2.vtk
12 // mpirun -np 4 ex13p -m ../data/fichera-q3.mesh
13 // mpirun -np 4 ex13p -m ../data/square-disc-nurbs.mesh
14 // mpirun -np 4 ex13p -m ../data/beam-hex-nurbs.mesh
15 // mpirun -np 4 ex13p -m ../data/amr-quad.mesh -o 2
16 // mpirun -np 4 ex13p -m ../data/amr-hex.mesh
17 // mpirun -np 4 ex13p -m ../data/mobius-strip.mesh -n 8 -o 2
18 // mpirun -np 4 ex13p -m ../data/klein-bottle.mesh -n 10 -o 2
19 //
20 // Description: This example code solves the Maxwell (electromagnetic)
21 // eigenvalue problem curl curl E = lambda E with homogeneous
22 // Dirichlet boundary conditions E x n = 0.
23 //
24 // We compute a number of the lowest nonzero eigenmodes by
25 // discretizing the curl curl operator using a Nedelec FE space of
26 // the specified order in 2D or 3D.
27 //
28 // The example highlights the use of the AME subspace eigenvalue
29 // solver from HYPRE, which uses LOBPCG and AMS internally.
30 // Reusing a single GLVis visualization window for multiple
31 // eigenfunctions is also illustrated.
32 //
33 // We recommend viewing examples 3 and 11 before viewing this
34 // example.
35 
36 #include "mfem.hpp"
37 #include <fstream>
38 #include <iostream>
39 
40 using namespace std;
41 using namespace mfem;
42 
43 int main(int argc, char *argv[])
44 {
45  // 1. Initialize MPI and HYPRE.
46  Mpi::Init(argc, argv);
47  int num_procs = Mpi::WorldSize();
48  int myid = Mpi::WorldRank();
49  Hypre::Init();
50 
51  // 2. Parse command-line options.
52  const char *mesh_file = "../data/beam-tet.mesh";
53  int ser_ref_levels = 2;
54  int par_ref_levels = 1;
55  int order = 1;
56  int nev = 5;
57  bool visualization = 1;
58  const char *device_config = "cpu";
59 
60  OptionsParser args(argc, argv);
61  args.AddOption(&mesh_file, "-m", "--mesh",
62  "Mesh file to use.");
63  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
64  "Number of times to refine the mesh uniformly in serial.");
65  args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
66  "Number of times to refine the mesh uniformly in parallel.");
67  args.AddOption(&order, "-o", "--order",
68  "Finite element order (polynomial degree) or -1 for"
69  " isoparametric space.");
70  args.AddOption(&nev, "-n", "--num-eigs",
71  "Number of desired eigenmodes.");
72  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
73  "--no-visualization",
74  "Enable or disable GLVis visualization.");
75  args.AddOption(&device_config, "-d", "--device",
76  "Device configuration string, see Device::Configure().");
77  args.Parse();
78  if (!args.Good())
79  {
80  if (myid == 0)
81  {
82  args.PrintUsage(cout);
83  }
84  return 1;
85  }
86  if (myid == 0)
87  {
88  args.PrintOptions(cout);
89  }
90 
91  // 3. Enable hardware devices such as GPUs, and programming models such as
92  // CUDA, OCCA, RAJA and OpenMP based on command line options.
93  Device device(device_config);
94  if (myid == 0) { device.Print(); }
95 
96  // 4. Read the (serial) mesh from the given mesh file on all processors. We
97  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
98  // and volume meshes with the same code.
99  Mesh *mesh = new Mesh(mesh_file, 1, 1);
100  int dim = mesh->Dimension();
101 
102  // 5. Refine the serial mesh on all processors to increase the resolution. In
103  // this example we do 'ref_levels' of uniform refinement (2 by default, or
104  // specified on the command line with -rs).
105  for (int lev = 0; lev < ser_ref_levels; lev++)
106  {
107  mesh->UniformRefinement();
108  }
109 
110  // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
111  // this mesh further in parallel to increase the resolution (1 time by
112  // default, or specified on the command line with -rp). Once the parallel
113  // mesh is defined, the serial mesh can be deleted.
114  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
115  delete mesh;
116  for (int lev = 0; lev < par_ref_levels; lev++)
117  {
118  pmesh->UniformRefinement();
119  }
120 
121  // 7. Define a parallel finite element space on the parallel mesh. Here we
122  // use the Nedelec finite elements of the specified order.
123  FiniteElementCollection *fec = new ND_FECollection(order, dim);
124  ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
125  HYPRE_BigInt size = fespace->GlobalTrueVSize();
126  if (myid == 0)
127  {
128  cout << "Number of unknowns: " << size << endl;
129  }
130 
131  // 8. Set up the parallel bilinear forms a(.,.) and m(.,.) on the finite
132  // element space. The first corresponds to the curl curl, while the second
133  // is a simple mass matrix needed on the right hand side of the
134  // generalized eigenvalue problem below. The boundary conditions are
135  // implemented by marking all the boundary attributes from the mesh as
136  // essential. The corresponding degrees of freedom are eliminated with
137  // special values on the diagonal to shift the Dirichlet eigenvalues out
138  // of the computational range. After serial and parallel assembly we
139  // extract the corresponding parallel matrices A and M.
140  ConstantCoefficient one(1.0);
141  Array<int> ess_bdr;
142  if (pmesh->bdr_attributes.Size())
143  {
144  ess_bdr.SetSize(pmesh->bdr_attributes.Max());
145  ess_bdr = 1;
146  }
147 
148  ParBilinearForm *a = new ParBilinearForm(fespace);
149  a->AddDomainIntegrator(new CurlCurlIntegrator(one));
150  if (pmesh->bdr_attributes.Size() == 0)
151  {
152  // Add a mass term if the mesh has no boundary, e.g. periodic mesh or
153  // closed surface.
154  a->AddDomainIntegrator(new VectorFEMassIntegrator(one));
155  }
156  a->Assemble();
157  a->EliminateEssentialBCDiag(ess_bdr, 1.0);
158  a->Finalize();
159 
160  ParBilinearForm *m = new ParBilinearForm(fespace);
162  m->Assemble();
163  // shift the eigenvalue corresponding to eliminated dofs to a large value
164  m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<double>::min());
165  m->Finalize();
166 
167  HypreParMatrix *A = a->ParallelAssemble();
169 
170  delete a;
171  delete m;
172 
173  // 9. Define and configure the AME eigensolver and the AMS preconditioner for
174  // A to be used within the solver. Set the matrices which define the
175  // generalized eigenproblem A x = lambda M x.
176  HypreAMS *ams = new HypreAMS(*A,fespace);
177  ams->SetPrintLevel(0);
178  ams->SetSingularProblem();
179 
180  HypreAME *ame = new HypreAME(MPI_COMM_WORLD);
181  ame->SetNumModes(nev);
182  ame->SetPreconditioner(*ams);
183  ame->SetMaxIter(100);
184  ame->SetTol(1e-8);
185  ame->SetPrintLevel(1);
186  ame->SetMassMatrix(*M);
187  ame->SetOperator(*A);
188 
189  // 10. Compute the eigenmodes and extract the array of eigenvalues. Define a
190  // parallel grid function to represent each of the eigenmodes returned by
191  // the solver.
192  Array<double> eigenvalues;
193  ame->Solve();
194  ame->GetEigenvalues(eigenvalues);
195  ParGridFunction x(fespace);
196 
197  // 11. Save the refined mesh and the modes in parallel. This output can be
198  // viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
199  {
200  ostringstream mesh_name, mode_name;
201  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
202 
203  ofstream mesh_ofs(mesh_name.str().c_str());
204  mesh_ofs.precision(8);
205  pmesh->Print(mesh_ofs);
206 
207  for (int i=0; i<nev; i++)
208  {
209  // convert eigenvector from HypreParVector to ParGridFunction
210  x = ame->GetEigenvector(i);
211 
212  mode_name << "mode_" << setfill('0') << setw(2) << i << "."
213  << setfill('0') << setw(6) << myid;
214 
215  ofstream mode_ofs(mode_name.str().c_str());
216  mode_ofs.precision(8);
217  x.Save(mode_ofs);
218  mode_name.str("");
219  }
220  }
221 
222  // 12. Send the solution by socket to a GLVis server.
223  if (visualization)
224  {
225  char vishost[] = "localhost";
226  int visport = 19916;
227  socketstream mode_sock(vishost, visport);
228  mode_sock.precision(8);
229 
230  for (int i=0; i<nev; i++)
231  {
232  if ( myid == 0 )
233  {
234  cout << "Eigenmode " << i+1 << '/' << nev
235  << ", Lambda = " << eigenvalues[i] << endl;
236  }
237 
238  // convert eigenvector from HypreParVector to ParGridFunction
239  x = ame->GetEigenvector(i);
240 
241  mode_sock << "parallel " << num_procs << " " << myid << "\n"
242  << "solution\n" << *pmesh << x << flush
243  << "window_title 'Eigenmode " << i+1 << '/' << nev
244  << ", Lambda = " << eigenvalues[i] << "'" << endl;
245 
246  char c;
247  if (myid == 0)
248  {
249  cout << "press (q)uit or (c)ontinue --> " << flush;
250  cin >> c;
251  }
252  MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
253 
254  if (c != 'c')
255  {
256  break;
257  }
258  }
259  mode_sock.close();
260  }
261 
262  // 13. Free the used memory.
263  delete ame;
264  delete ams;
265  delete M;
266  delete A;
267 
268  delete fespace;
269  delete fec;
270  delete pmesh;
271 
272  return 0;
273 }
int visport
The Auxiliary-space Maxwell Solver in hypre.
Definition: hypre.hpp:1743
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
void SetMassMatrix(const HypreParMatrix &M)
Definition: hypre.cpp:6471
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
Integrator for (curl u, curl v) for Nedelec elements.
void GetEigenvalues(Array< double > &eigenvalues) const
Collect the converged eigenvalues.
Definition: hypre.cpp:6490
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
void SetPreconditioner(HypreSolver &precond)
Definition: hypre.cpp:6450
void SetTol(double tol)
Definition: hypre.cpp:6419
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:279
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
Abstract parallel finite element space.
Definition: pfespace.hpp:28
STL namespace.
int main(int argc, char *argv[])
Definition: ex13p.cpp:43
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
int close()
Close the socketstream.
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
char vishost[]
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:10232
void Assemble(int skip_zeros=1)
Assemble the local matrix.
void EliminateEssentialBCDiag(const Array< int > &bdr_attr_is_ess, double value)
Perform elimination and set the diagonal entry to the given value.
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:281
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:5645
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
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
void Solve()
Solve the eigenproblem.
Definition: hypre.cpp:6478
HYPRE_Int HYPRE_BigInt
void SetNumModes(int num_eigs)
Definition: hypre.cpp:6411
void SetMaxIter(int max_iter)
Definition: hypre.cpp:6435
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:909
double a
Definition: lissajous.cpp:41
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
int dim
Definition: ex24.cpp:53
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
void SetPrintLevel(int logging)
Definition: hypre.cpp:6441
Class for parallel bilinear form.
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
Arbitrary order H(curl)-conforming Nedelec finite elements.
Definition: fe_coll.hpp:454
void SetOperator(const HypreParMatrix &A)
Definition: hypre.cpp:6456
void Print(std::ostream &out=mfem::out) const override
Definition: pmesh.cpp:4825
Class for parallel grid function.
Definition: pgridfunc.hpp:32
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:121
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
Class for parallel meshes.
Definition: pmesh.hpp:32
void SetSingularProblem()
Set this option when solving a curl-curl problem with zero mass term.
Definition: hypre.hpp:1802
const HypreParVector & GetEigenvector(unsigned int i) const
Extract a single eigenvector.
Definition: hypre.cpp:6514