MFEM  v4.5.2
Finite element discretization library
ex11p.cpp
Go to the documentation of this file.
1 // MFEM Example 11 - Parallel Version
2 // PETSc Modification
3 //
4 // Compile with: make ex11p
5 //
6 // Sample runs: mpirun -np 4 ex11p -m ../../data/star.mesh
7 // mpirun -np 4 ex11p -m ../../data/star.mesh --slepcopts rc_ex11p_lobpcg
8 // mpirun -np 4 ex11p -m ../../data/star.mesh --slepcopts rc_ex11p_gd
9 //
10 // Description: This example code demonstrates the use of MFEM to solve the
11 // eigenvalue problem -Delta u = lambda u with homogeneous
12 // Dirichlet boundary conditions.
13 //
14 // We compute a number of the lowest eigenmodes by discretizing
15 // the Laplacian and Mass operators using a FE space of the
16 // specified order, or an isoparametric/isogeometric space if
17 // order < 1 (quadratic for quadratic curvilinear mesh, NURBS for
18 // NURBS mesh, etc.)
19 //
20 // The example demonstrates the use of the SLEPc eigensolver as an
21 // alternative to the LOBPCG eigenvalue solver. The shift and
22 // invert spectral transformation is used to help the convergence
23 // to the smaller eigenvalues. Alternative solver parameters can
24 // be passed in a file with "-slepcopts".
25 //
26 // Reusing a single GLVis visualization window for multiple
27 // eigenfunctions is also illustrated.
28 //
29 // We recommend viewing Example 1 before viewing this example.
30 
31 #include "mfem.hpp"
32 #include <fstream>
33 #include <iostream>
34 
35 #ifndef MFEM_USE_SLEPC
36 #error This examples requires that MFEM is build with MFEM_USE_SLEPC=YES
37 #endif
38 
39 using namespace std;
40 using namespace mfem;
41 
42 int main(int argc, char *argv[])
43 {
44  // 1. Initialize MPI and HYPRE.
45  Mpi::Init(argc, argv);
46  int num_procs = Mpi::WorldSize();
47  int myid = Mpi::WorldRank();
48  Hypre::Init();
49 
50  // 2. Parse command-line options.
51  const char *mesh_file = "../../data/star.mesh";
52  int ser_ref_levels = 2;
53  int par_ref_levels = 1;
54  int order = 1;
55  int nev = 5;
56  int seed = 75;
57  bool slu_solver = false;
58  bool sp_solver = false;
59  bool visualization = 1;
60  bool use_slepc = true;
61  const char *slepcrc_file = "";
62  const char *device_config = "cpu";
63 
64  OptionsParser args(argc, argv);
65  args.AddOption(&mesh_file, "-m", "--mesh",
66  "Mesh file to use.");
67  args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
68  "Number of times to refine the mesh uniformly in serial.");
69  args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
70  "Number of times to refine the mesh uniformly in parallel.");
71  args.AddOption(&order, "-o", "--order",
72  "Finite element order (polynomial degree) or -1 for"
73  " isoparametric space.");
74  args.AddOption(&nev, "-n", "--num-eigs",
75  "Number of desired eigenmodes.");
76  args.AddOption(&seed, "-s", "--seed",
77  "Random seed used to initialize LOBPCG.");
78 #ifdef MFEM_USE_SUPERLU
79  args.AddOption(&slu_solver, "-slu", "--superlu", "-no-slu",
80  "--no-superlu", "Use the SuperLU Solver.");
81 #endif
82 #ifdef MFEM_USE_STRUMPACK
83  args.AddOption(&sp_solver, "-sp", "--strumpack", "-no-sp",
84  "--no-strumpack", "Use the STRUMPACK Solver.");
85 #endif
86  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
87  "--no-visualization",
88  "Enable or disable GLVis visualization.");
89  args.AddOption(&use_slepc, "-useslepc","--useslepc","-no-slepc",
90  "--no-slepc","Use or not SLEPc to solve the eigenvalue problem");
91  args.AddOption(&slepcrc_file, "-slepcopts", "--slepcopts",
92  "SlepcOptions file to use.");
93  args.AddOption(&device_config, "-d", "--device",
94  "Device configuration string, see Device::Configure().");
95  args.Parse();
96  if (slu_solver && sp_solver)
97  {
98  if (myid == 0)
99  cout << "WARNING: Both SuperLU and STRUMPACK have been selected,"
100  << " please choose either one." << endl
101  << " Defaulting to SuperLU." << endl;
102  sp_solver = false;
103  }
104  // The command line options are also passed to the STRUMPACK
105  // solver. So do not exit if some options are not recognized.
106  if (!sp_solver)
107  {
108  if (!args.Good())
109  {
110  if (myid == 0)
111  {
112  args.PrintUsage(cout);
113  }
114  return 1;
115  }
116  }
117  if (myid == 0)
118  {
119  args.PrintOptions(cout);
120  }
121 
122  // 2b. Enable hardware devices such as GPUs, and programming models such as
123  // CUDA, OCCA, RAJA and OpenMP based on command line options.
124  Device device(device_config);
125  if (myid == 0) { device.Print(); }
126 
127  // 2c. We initialize SLEPc. This internally initializes PETSc as well.
128  MFEMInitializeSlepc(NULL,NULL,slepcrc_file,NULL);
129 
130  // 3. Read the (serial) mesh from the given mesh file on all processors. We
131  // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
132  // and volume meshes with the same code.
133  Mesh *mesh = new Mesh(mesh_file, 1, 1);
134  int dim = mesh->Dimension();
135 
136  // 4. Refine the serial mesh on all processors to increase the resolution. In
137  // this example we do 'ref_levels' of uniform refinement (2 by default, or
138  // specified on the command line with -rs).
139  for (int lev = 0; lev < ser_ref_levels; lev++)
140  {
141  mesh->UniformRefinement();
142  }
143 
144  // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
145  // this mesh further in parallel to increase the resolution (1 time by
146  // default, or specified on the command line with -rp). Once the parallel
147  // mesh is defined, the serial mesh can be deleted.
148  ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
149  delete mesh;
150  for (int lev = 0; lev < par_ref_levels; lev++)
151  {
152  pmesh->UniformRefinement();
153  }
154 
155  // 6. Define a parallel finite element space on the parallel mesh. Here we
156  // use continuous Lagrange finite elements of the specified order. If
157  // order < 1, we instead use an isoparametric/isogeometric space.
159  if (order > 0)
160  {
161  fec = new H1_FECollection(order, dim);
162  }
163  else if (pmesh->GetNodes())
164  {
165  fec = pmesh->GetNodes()->OwnFEC();
166  }
167  else
168  {
169  fec = new H1_FECollection(order = 1, dim);
170  }
171  ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
172  HYPRE_BigInt size = fespace->GlobalTrueVSize();
173  if (myid == 0)
174  {
175  cout << "Number of unknowns: " << size << endl;
176  }
177 
178  // 7. Set up the parallel bilinear forms a(.,.) and m(.,.) on the finite
179  // element space. The first corresponds to the Laplacian operator -Delta,
180  // while the second is a simple mass matrix needed on the right hand side
181  // of the generalized eigenvalue problem below. The boundary conditions
182  // are implemented by elimination with special values on the diagonal to
183  // shift the Dirichlet eigenvalues out of the computational range. After
184  // serial and parallel assembly we extract the corresponding parallel
185  // matrices A and M.
186  ConstantCoefficient one(1.0);
187  Array<int> ess_bdr;
188  if (pmesh->bdr_attributes.Size())
189  {
190  ess_bdr.SetSize(pmesh->bdr_attributes.Max());
191  ess_bdr = 1;
192  }
193 
194  ParBilinearForm *a = new ParBilinearForm(fespace);
195  a->AddDomainIntegrator(new DiffusionIntegrator(one));
196  if (pmesh->bdr_attributes.Size() == 0)
197  {
198  // Add a mass term if the mesh has no boundary, e.g. periodic mesh or
199  // closed surface.
200  a->AddDomainIntegrator(new MassIntegrator(one));
201  }
202  a->Assemble();
203  a->EliminateEssentialBCDiag(ess_bdr, 1.0);
204  a->Finalize();
205 
206  ParBilinearForm *m = new ParBilinearForm(fespace);
207  m->AddDomainIntegrator(new MassIntegrator(one));
208  m->Assemble();
209  // shift the eigenvalue corresponding to eliminated dofs to a large value
210  m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<double>::min());
211  m->Finalize();
212 
213  PetscParMatrix *pA = NULL, *pM = NULL;
214  HypreParMatrix *A = NULL, *M = NULL;
215  Operator::Type tid =
216  !use_slepc ? Operator::Hypre_ParCSR : Operator::PETSC_MATAIJ;
217  OperatorHandle Ah(tid), Mh(tid);
218 
219  a->ParallelAssemble(Ah);
220  if (!use_slepc) { Ah.Get(A); }
221  else { Ah.Get(pA); }
222  Ah.SetOperatorOwner(false);
223 
224  m->ParallelAssemble(Mh);
225  if (!use_slepc) {Mh.Get(M); }
226  else {Mh.Get(pM); }
227  Mh.SetOperatorOwner(false);
228 
229 #if defined(MFEM_USE_SUPERLU) || defined(MFEM_USE_STRUMPACK)
230  Operator * Arow = NULL;
231 #ifdef MFEM_USE_SUPERLU
232  if (slu_solver)
233  {
234  Arow = new SuperLURowLocMatrix(*A);
235  }
236 #endif
237 #ifdef MFEM_USE_STRUMPACK
238  if (sp_solver)
239  {
240  Arow = new STRUMPACKRowLocMatrix(*A);
241  }
242 #endif
243 #endif
244 
245  delete a;
246  delete m;
247 
248  // 8. Define and configure the LOBPCG eigensolver and the BoomerAMG
249  // preconditioner for A to be used within the solver. Set the matrices
250  // which define the generalized eigenproblem A x = lambda M x.
251  Solver * precond = NULL;
252  if (!use_slepc)
253  {
254  if (!slu_solver && !sp_solver)
255  {
256  HypreBoomerAMG * amg = new HypreBoomerAMG(*A);
257  amg->SetPrintLevel(0);
258  precond = amg;
259  }
260  else
261  {
262 #ifdef MFEM_USE_SUPERLU
263  if (slu_solver)
264  {
265  SuperLUSolver * superlu = new SuperLUSolver(MPI_COMM_WORLD);
266  superlu->SetPrintStatistics(false);
267  superlu->SetSymmetricPattern(true);
269  superlu->SetOperator(*Arow);
270  precond = superlu;
271  }
272 #endif
273 #ifdef MFEM_USE_STRUMPACK
274  if (sp_solver)
275  {
276  STRUMPACKSolver * strumpack = new STRUMPACKSolver(argc, argv, MPI_COMM_WORLD);
277  strumpack->SetPrintFactorStatistics(true);
278  strumpack->SetPrintSolveStatistics(false);
279  strumpack->SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
280  strumpack->SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
281  strumpack->DisableMatching();
282  strumpack->SetOperator(*Arow);
283  strumpack->SetFromCommandLine();
284  precond = strumpack;
285  }
286 #endif
287  }
288  }
289 
290  HypreLOBPCG * lobpcg = NULL;
291  SlepcEigenSolver * slepc = NULL;
292  if (!use_slepc)
293  {
294 
295  lobpcg = new HypreLOBPCG(MPI_COMM_WORLD);
296  lobpcg->SetNumModes(nev);
297  lobpcg->SetRandomSeed(seed);
298  lobpcg->SetPreconditioner(*precond);
299  lobpcg->SetMaxIter(200);
300  lobpcg->SetTol(1e-8);
301  lobpcg->SetPrecondUsageMode(1);
302  lobpcg->SetPrintLevel(1);
303  lobpcg->SetMassMatrix(*M);
304  lobpcg->SetOperator(*A);
305  }
306  else
307  {
308  slepc = new SlepcEigenSolver(MPI_COMM_WORLD);
309  slepc->SetNumModes(nev);
310  slepc->SetWhichEigenpairs(SlepcEigenSolver::TARGET_REAL);
311  slepc->SetTarget(0.0);
312  slepc->SetSpectralTransformation(SlepcEigenSolver::SHIFT_INVERT);
313  slepc->SetOperators(*pA,*pM);
314  }
315 
316  // 9. Compute the eigenmodes and extract the array of eigenvalues. Define a
317  // parallel grid function to represent each of the eigenmodes returned by
318  // the solver.
319  Array<double> eigenvalues;
320  if (!use_slepc)
321  {
322  lobpcg->Solve();
323  lobpcg->GetEigenvalues(eigenvalues);
324  }
325  else
326  {
327  slepc->Solve();
328  eigenvalues.SetSize(nev);
329  for (int i=0; i<nev; i++)
330  {
331  slepc->GetEigenvalue(i,eigenvalues[i]);
332  }
333  }
334  Vector temp(fespace->GetTrueVSize());
335  ParGridFunction x(fespace);
336 
337  // 10. Save the refined mesh and the modes in parallel. This output can be
338  // viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
339  {
340  ostringstream mesh_name, mode_name;
341  mesh_name << "mesh." << setfill('0') << setw(6) << myid;
342 
343  ofstream mesh_ofs(mesh_name.str().c_str());
344  mesh_ofs.precision(8);
345  pmesh->Print(mesh_ofs);
346 
347  for (int i=0; i<nev; i++)
348  {
349  // convert eigenvector from HypreParVector to ParGridFunction
350  if (!use_slepc)
351  {
352  x = lobpcg->GetEigenvector(i);
353  }
354  else
355  {
356  slepc->GetEigenvector(i,temp);
357  x.Distribute(temp);
358  }
359 
360  mode_name << "mode_" << setfill('0') << setw(2) << i << "."
361  << setfill('0') << setw(6) << myid;
362 
363  ofstream mode_ofs(mode_name.str().c_str());
364  mode_ofs.precision(8);
365  x.Save(mode_ofs);
366  mode_name.str("");
367  }
368  }
369 
370  // 11. Send the solution by socket to a GLVis server.
371  if (visualization)
372  {
373  char vishost[] = "localhost";
374  int visport = 19916;
375  socketstream mode_sock(vishost, visport);
376  mode_sock.precision(8);
377 
378  for (int i=0; i<nev; i++)
379  {
380  if ( myid == 0 )
381  {
382  cout << "Eigenmode " << i+1 << '/' << nev
383  << ", Lambda = " << eigenvalues[i] << endl;
384  }
385 
386  // convert eigenvector from HypreParVector to ParGridFunction
387  if (!use_slepc)
388  {
389  x = lobpcg->GetEigenvector(i);
390  }
391  else
392  {
393  slepc->GetEigenvector(i,temp);
394  x.Distribute(temp);
395  }
396 
397  mode_sock << "parallel " << num_procs << " " << myid << "\n"
398  << "solution\n" << *pmesh << x << flush
399  << "window_title 'Eigenmode " << i+1 << '/' << nev
400  << ", Lambda = " << eigenvalues[i] << "'" << endl;
401 
402  char c;
403  if (myid == 0)
404  {
405  cout << "press (q)uit or (c)ontinue --> " << flush;
406  cin >> c;
407  }
408  MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
409 
410  if (c != 'c')
411  {
412  break;
413  }
414  }
415  mode_sock.close();
416  }
417 
418  // 12. Free the used memory.
419  delete lobpcg;
420  delete slepc;
421  delete precond;
422  delete M;
423  delete A;
424  delete pA;
425  delete pM;
426 #if defined(MFEM_USE_SUPERLU) || defined(MFEM_USE_STRUMPACK)
427  delete Arow;
428 #endif
429  delete fespace;
430  if (order > 0)
431  {
432  delete fec;
433  }
434  delete pmesh;
435 
436  // We finalize SLEPc
438 
439  return 0;
440 }
void GetEigenvalues(Array< double > &eigenvalues) const
Collect the converged eigenvalues.
Definition: hypre.cpp:6230
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
void MFEMInitializeSlepc()
Definition: slepc.cpp:29
int Dimension() const
Definition: mesh.hpp:1047
Wrapper for PETSc&#39;s matrix class.
Definition: petsc.hpp:315
const HypreParVector & GetEigenvector(unsigned int i) const
Extract a single eigenvector.
Definition: hypre.cpp:6242
void GetEigenvalue(unsigned int i, double &lr) const
Get the corresponding eigenvalue.
Definition: slepc.cpp:147
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
void SetPreconditioner(Solver &precond)
Definition: hypre.cpp:6164
void GetEigenvector(unsigned int i, Vector &vr) const
Get the corresponding eigenvector.
Definition: slepc.cpp:158
void Print(std::ostream &out=mfem::out)
Print the configuration of the MFEM virtual device object.
Definition: device.cpp:279
void SetMassMatrix(Operator &M)
Definition: hypre.cpp:6220
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
void Get(OpType *&A) const
Return the Operator pointer statically cast to a given OpType.
Definition: handle.hpp:114
void SetPrintLevel(int logging)
Definition: hypre.cpp:6149
STL namespace.
void SetTarget(double target)
Definition: slepc.cpp:229
void SetSymmetricPattern(bool sym)
Definition: superlu.cpp:425
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
void SetTol(double tol)
Definition: hypre.cpp:6127
int close()
Close the socketstream.
void SetOperator(const Operator &op)
Set/update the solver for the given operator.
Definition: superlu.cpp:591
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1590
void SetSpectralTransformation(SpectralTransformation transformation)
Definition: slepc.cpp:234
void SetMaxIter(int max_iter)
Definition: hypre.cpp:6143
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[]
void SetColumnPermutation(superlu::ColPerm col_perm)
Definition: superlu.cpp:345
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition: mesh.cpp:9878
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1670
constexpr int visport
void Assemble(int skip_zeros=1)
Assemble the local matrix.
void SetNumModes(int num_eigs)
Set the number of required eigenmodes.
Definition: slepc.cpp:124
void EliminateEssentialBCDiag(const Array< int > &bdr_attr_is_ess, double value)
Perform elimination and set the diagonal entry to the given value.
void SetPrintStatistics(bool print_stat)
Definition: superlu.cpp:327
void Solve()
Solve the eigenvalue problem for the specified number of eigenvalues.
Definition: slepc.cpp:130
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:285
void MFEMFinalizeSlepc()
Definition: slepc.cpp:53
Type
Enumeration defining IDs for some classes derived from Operator.
Definition: operator.hpp:283
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 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:684
HYPRE_Int HYPRE_BigInt
void SetRandomSeed(int s)
Definition: hypre.hpp:2004
virtual int GetTrueVSize() const
Return the number of local vector true dofs.
Definition: pfespace.hpp:289
double a
Definition: lissajous.cpp:41
void SetOperatorOwner(bool own=true)
Set the ownership flag for the held Operator.
Definition: handle.hpp:120
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void SetWhichEigenpairs(Which which)
Definition: slepc.cpp:195
int dim
Definition: ex24.cpp:53
void SetOperator(Operator &A)
Definition: hypre.cpp:6173
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
Class for parallel bilinear form.
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
Vector data type.
Definition: vector.hpp:60
void SetPrecondUsageMode(int pcg_mode)
Definition: hypre.cpp:6158
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:252
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7949
void Print(std::ostream &out=mfem::out) const override
Definition: pmesh.cpp:4839
Base class for solvers.
Definition: operator.hpp:682
Class for parallel grid function.
Definition: pgridfunc.hpp:32
int main(int argc, char *argv[])
Definition: ex11p.cpp:58
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:121
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
void SetNumModes(int num_eigs)
Definition: hypre.hpp:2002
Class for parallel meshes.
Definition: pmesh.hpp:32
void Solve()
Solve the eigenproblem.
Definition: hypre.cpp:6287
void SetOperators(const PetscParMatrix &op, const PetscParMatrix &opB)
Set operator for generalized eigenvalue problem.
Definition: slepc.cpp:93