libMesh
eigenproblems_ex2.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // <h1>Eigenproblems Example 2 - Solving a generalized Eigen Problem</h1>
21 // \author Steffen Petersen
22 // \date 2006
23 //
24 // This example shows how the previous EigenSolver example
25 // can be adapted to solve generalized eigenvalue problems.
26 //
27 // For solving eigen problems, libMesh interfaces
28 // SLEPc (www.grycap.upv.es/slepc/) which again is based on PETSc.
29 // Hence, this example will only work if the library is compiled
30 // with SLEPc support enabled.
31 //
32 // In this example some eigenvalues for a generalized symmetric
33 // eigenvalue problem A*x=lambda*B*x are computed, where the
34 // matrices A and B are assembled according to stiffness and
35 // mass matrix, respectively.
36 
37 // libMesh include files.
38 #include "libmesh/libmesh.h"
39 #include "libmesh/mesh.h"
40 #include "libmesh/mesh_generation.h"
41 #include "libmesh/exodusII_io.h"
42 #include "libmesh/eigen_system.h"
43 #include "libmesh/equation_systems.h"
44 #include "libmesh/slepc_eigen_solver.h"
45 #include "libmesh/fe.h"
46 #include "libmesh/quadrature_gauss.h"
47 #include "libmesh/dense_matrix.h"
48 #include "libmesh/sparse_matrix.h"
49 #include "libmesh/numeric_vector.h"
50 #include "libmesh/dof_map.h"
51 
52 // Bring in everything from the libMesh namespace
53 using namespace libMesh;
54 
55 
56 // Function prototype. This is the function that will assemble
57 // the eigen system. Here, we will simply assemble a mass matrix.
59  const std::string & system_name);
60 
61 
62 
63 int main (int argc, char ** argv)
64 {
65  // Initialize libMesh and the dependent libraries.
66  LibMeshInit init (argc, argv);
67 
68  // This example is designed for the SLEPc eigen solver interface.
69 #ifndef LIBMESH_HAVE_SLEPC
70  if (init.comm().rank() == 0)
71  libMesh::err << "ERROR: This example requires libMesh to be\n"
72  << "compiled with SLEPc eigen solvers support!"
73  << std::endl;
74 
75  return 0;
76 #else
77 
78 #ifdef LIBMESH_DEFAULT_SINGLE_PRECISION
79  // SLEPc currently gives us a nasty crash with Real==float
80  libmesh_example_requires(false, "--disable-singleprecision");
81 #endif
82 
83  // Check for proper usage.
84  if (argc < 3)
85  libmesh_error_msg("\nUsage: " << argv[0] << " -n <number of eigen values>");
86 
87  // Tell the user what we are doing.
88  else
89  {
90  libMesh::out << "Running " << argv[0];
91 
92  for (int i=1; i<argc; i++)
93  libMesh::out << " " << argv[i];
94 
95  libMesh::out << std::endl << std::endl;
96  }
97 
98  // Get the number of eigen values to be computed from argv[2]
99  const unsigned int nev = std::atoi(argv[2]);
100 
101  // Skip this 2D example if libMesh was compiled as 1D-only.
102  libmesh_example_requires(2 <= LIBMESH_DIM, "2D support");
103 
104  // Create a mesh, with dimension to be overridden later, on the
105  // default MPI communicator.
106  Mesh mesh(init.comm());
107 
108  // Use the internal mesh generator to create a uniform
109  // 2D grid on a square.
111  20, 20,
112  -1., 1.,
113  -1., 1.,
114  QUAD4);
115 
116  // Print information about the mesh to the screen.
117  mesh.print_info();
118 
119  // Create an equation systems object.
120  EquationSystems equation_systems (mesh);
121 
122  // Create a EigenSystem named "Eigensystem" and (for convenience)
123  // use a reference to the system we create.
124  EigenSystem & eigen_system =
125  equation_systems.add_system<EigenSystem> ("Eigensystem");
126 
127  // Declare the system variables.
128  // Adds the variable "p" to "Eigensystem". "p"
129  // will be approximated using second-order approximation.
130  eigen_system.add_variable("p", FIRST);
131 
132  // Give the system a pointer to the matrix assembly
133  // function defined below.
135 
136  // Set necessary parameters used in EigenSystem::solve(),
137  // i.e. the number of requested eigenpairs nev and the number
138  // of basis vectors ncv used in the solution algorithm. Note that
139  // ncv >= nev must hold and ncv >= 2*nev is recommended.
140  equation_systems.parameters.set<unsigned int>("eigenpairs") = nev;
141  equation_systems.parameters.set<unsigned int>("basis vectors") = nev*3;
142 
143  // You may optionally change the default eigensolver used by SLEPc.
144  // The Krylov-Schur method is mathematically equivalent to implicitly
145  // restarted Arnoldi, the method of Arpack, so there is currently no
146  // point in using SLEPc with Arpack.
147  // ARNOLDI = default in SLEPc 2.3.1 and earlier
148  // KRYLOVSCHUR default in SLEPc 2.3.2 and later
149  // eigen_system.eigen_solver->set_eigensolver_type(KRYLOVSCHUR);
150 
151  // Set the solver tolerance and the maximum number of iterations.
152  equation_systems.parameters.set<Real>("linear solver tolerance") = pow(TOLERANCE, 5./3.);
153  equation_systems.parameters.set<unsigned int>("linear solver maximum iterations") = 1000;
154 
155  // Set the type of the problem, here we deal with
156  // a generalized Hermitian problem.
157  eigen_system.set_eigenproblem_type(GHEP);
158 
159  // Set the eigenvalues to be computed. Note that not
160  // all solvers in SLEPc support this capability.
161  eigen_system.eigen_solver->set_position_of_spectrum(2.3);
162 
163  // Initialize the data structures for the equation system.
164  equation_systems.init();
165 
166  // Prints information about the system to the screen.
167  equation_systems.print_info();
168 
169 #if SLEPC_VERSION_LESS_THAN(3,1,0)
170  libmesh_error_msg("SLEPc 3.1 is required to call EigenSolver::set_initial_space()");
171 #else
172  // Get the SLEPc solver object and set initial guess for one basis vector
173  // this has to be done _after_ the EquationSystems object is initialized
174  UniquePtr<EigenSolver<Number>> & slepc_eps = eigen_system.eigen_solver;
175  NumericVector<Number> & initial_space = eigen_system.add_vector("initial_space");
176  initial_space.add(1.0);
177  slepc_eps->set_initial_space(initial_space);
178 #endif
179 
180  // Solve the system "Eigensystem".
181  eigen_system.solve();
182 
183  // Get the number of converged eigen pairs.
184  unsigned int nconv = eigen_system.get_n_converged();
185 
186  libMesh::out << "Number of converged eigenpairs: "
187  << nconv
188  << "\n"
189  << std::endl;
190 
191  // Get the last converged eigenpair
192  if (nconv != 0)
193  {
194  eigen_system.get_eigenpair(nconv-1);
195 
196 #ifdef LIBMESH_HAVE_EXODUS_API
197  // Write the eigen vector to file.
198  ExodusII_IO (mesh).write_equation_systems ("out.e", equation_systems);
199 #endif // #ifdef LIBMESH_HAVE_EXODUS_API
200  }
201  else
202  {
203  libMesh::out << "WARNING: Solver did not converge!\n" << nconv << std::endl;
204  }
205 
206 #endif // LIBMESH_HAVE_SLEPC
207 
208  // All done.
209  return 0;
210 }
211 
212 
213 
215  const std::string & libmesh_dbg_var(system_name))
216 {
217 
218  // It is a good idea to make sure we are assembling
219  // the proper system.
220  libmesh_assert_equal_to (system_name, "Eigensystem");
221 
222 #ifdef LIBMESH_HAVE_SLEPC
223 
224  // Get a constant reference to the mesh object.
225  const MeshBase & mesh = es.get_mesh();
226 
227  // The dimension that we are running.
228  const unsigned int dim = mesh.mesh_dimension();
229 
230  // Get a reference to our system.
231  EigenSystem & eigen_system = es.get_system<EigenSystem> ("Eigensystem");
232 
233  // Get a constant reference to the Finite Element type
234  // for the first (and only) variable in the system.
235  FEType fe_type = eigen_system.get_dof_map().variable_type(0);
236 
237  // A reference to the two system matrices
238  SparseMatrix<Number> & matrix_A = *eigen_system.matrix_A;
239  SparseMatrix<Number> & matrix_B = *eigen_system.matrix_B;
240 
241  // Build a Finite Element object of the specified type. Since the
242  // FEBase::build() member dynamically creates memory we will
243  // store the object as a UniquePtr<FEBase>. This can be thought
244  // of as a pointer that will clean up after itself.
245  UniquePtr<FEBase> fe (FEBase::build(dim, fe_type));
246 
247  // A Gauss quadrature rule for numerical integration.
248  // Use the default quadrature order.
249  QGauss qrule (dim, fe_type.default_quadrature_order());
250 
251  // Tell the finite element object to use our quadrature rule.
252  fe->attach_quadrature_rule (&qrule);
253 
254  // The element Jacobian * quadrature weight at each integration point.
255  const std::vector<Real> & JxW = fe->get_JxW();
256 
257  // The element shape functions evaluated at the quadrature points.
258  const std::vector<std::vector<Real>> & phi = fe->get_phi();
259 
260  // The element shape function gradients evaluated at the quadrature
261  // points.
262  const std::vector<std::vector<RealGradient>> & dphi = fe->get_dphi();
263 
264  // A reference to the DofMap object for this system. The DofMap
265  // object handles the index translation from node and element numbers
266  // to degree of freedom numbers.
267  const DofMap & dof_map = eigen_system.get_dof_map();
268 
269  // The element mass and stiffness matrices.
272 
273  // This vector will hold the degree of freedom indices for
274  // the element. These define where in the global system
275  // the element degrees of freedom get mapped.
276  std::vector<dof_id_type> dof_indices;
277 
278  // Now we will loop over all the elements in the mesh that
279  // live on the local processor. We will compute the element
280  // matrix and right-hand-side contribution. In case users
281  // later modify this program to include refinement, we will
282  // be safe and will only consider the active elements;
283  // hence we use a variant of the active_elem_iterator.
284  for (const auto & elem : mesh.active_local_element_ptr_range())
285  {
286  // Get the degree of freedom indices for the
287  // current element. These define where in the global
288  // matrix and right-hand-side this element will
289  // contribute to.
290  dof_map.dof_indices (elem, dof_indices);
291 
292  // Compute the element-specific data for the current
293  // element. This involves computing the location of the
294  // quadrature points (q_point) and the shape functions
295  // (phi, dphi) for the current element.
296  fe->reinit (elem);
297 
298  // Zero the element matrices before
299  // summing them. We use the resize member here because
300  // the number of degrees of freedom might have changed from
301  // the last element. Note that this will be the case if the
302  // element type is different (i.e. the last element was a
303  // triangle, now we are on a quadrilateral).
304  Ke.resize (dof_indices.size(), dof_indices.size());
305  Me.resize (dof_indices.size(), dof_indices.size());
306 
307  // Now loop over the quadrature points. This handles
308  // the numeric integration.
309  //
310  // We will build the element matrix. This involves
311  // a double loop to integrate the test functions (i) against
312  // the trial functions (j).
313  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
314  for (std::size_t i=0; i<phi.size(); i++)
315  for (std::size_t j=0; j<phi.size(); j++)
316  {
317  Me(i,j) += JxW[qp]*phi[i][qp]*phi[j][qp];
318  Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
319  }
320 
321  // On an unrefined mesh, constrain_element_matrix does
322  // nothing. If this assembly function is ever repurposed to
323  // run on a refined mesh, getting the hanging node constraints
324  // right will be important. Note that, even with
325  // asymmetric_constraint_rows = false, the constrained dof
326  // diagonals still exist in the matrix, with diagonal entries
327  // that are there to ensure non-singular matrices for linear
328  // solves but which would generate positive non-physical
329  // eigenvalues for eigensolves.
330  // dof_map.constrain_element_matrix(Ke, dof_indices, false);
331  // dof_map.constrain_element_matrix(Me, dof_indices, false);
332 
333  // Finally, simply add the element contribution to the
334  // overall matrices A and B.
335  matrix_A.add_matrix (Ke, dof_indices);
336  matrix_B.add_matrix (Me, dof_indices);
337  } // end of element loop
338 
339 
340 #else
341  // Avoid compiler warnings
342  libmesh_ignore(es);
343 #endif // LIBMESH_HAVE_SLEPC
344 }
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
OStreamProxy err
virtual std::pair< Real, Real > get_eigenpair(dof_id_type i)
Definition: eigen_system.C:257
This is the EquationSystems class.
const FEType & variable_type(const unsigned int c) const
Definition: dof_map.h:1697
unsigned int dim
The ExodusII_IO class implements reading meshes in the ExodusII file format from Sandia National Labs...
Definition: exodusII_io.h:52
void assemble_mass(EquationSystems &es, const std::string &system_name)
unsigned int add_variable(const std::string &var, const FEType &type, const std::set< subdomain_id_type > *const active_subdomains=libmesh_nullptr)
Adds the variable var to the list of variables for this system.
Definition: system.C:1101
UniquePtr< SparseMatrix< Number > > matrix_B
A second system matrix for generalized eigenvalue problems.
Definition: eigen_system.h:155
MeshBase & mesh
unsigned int get_n_converged() const
Definition: eigen_system.h:124
void build_square(UnstructuredMesh &mesh, const unsigned int nx, const unsigned int ny, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 2D meshes.
static const Real TOLERANCE
The LibMeshInit class, when constructed, initializes the dependent libraries (e.g.
Definition: libmesh.h:62
The libMesh namespace provides an interface to certain functionality in the library.
UniquePtr< SparseMatrix< Number > > matrix_A
The system matrix for standard eigenvalue problems.
Definition: eigen_system.h:150
This is the MeshBase class.
Definition: mesh_base.h:68
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual SimpleRange< element_iterator > active_local_element_ptr_range()=0
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols)=0
Add the full matrix dm to the SparseMatrix.
NumericVector< Number > & add_vector(const std::string &vec_name, const bool projections=true, const ParallelType type=PARALLEL)
Adds the additional vector vec_name to this system.
Definition: system.C:681
const DofMap & get_dof_map() const
Definition: system.h:2030
virtual void solve() libmesh_override
Assembles & solves the eigen system.
Definition: eigen_system.C:199
void print_info(std::ostream &os=libMesh::out) const
Prints information about the equation systems, by default to libMesh::out.
void attach_assemble_function(void fptr(EquationSystems &es, const std::string &name))
Register a user function to use in assembling the system matrix and RHS.
Definition: system.C:1795
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
virtual System & add_system(const std::string &system_type, const std::string &name)
Add the system of type system_type named name to the systems array.
const Parallel::Communicator & comm() const
Definition: libmesh.h:81
PetscErrorCode Vec Mat libmesh_dbg_var(j)
double pow(double a, int b)
virtual void write_equation_systems(const std::string &, const EquationSystems &, const std::set< std::string > *system_names=libmesh_nullptr)
This method implements writing a mesh with data to a specified file where the data is taken from the ...
Definition: mesh_output.C:31
void libmesh_ignore(const T &)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
int main(int argc, char **argv)
UniquePtr< EigenSolver< Number > > eigen_solver
The EigenSolver, defining which interface, i.e solver package to use.
Definition: eigen_system.h:161
T & set(const std::string &)
Definition: parameters.h:469
OStreamProxy out
unsigned int mesh_dimension() const
Definition: mesh_base.C:148
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
This class implements specific orders of Gauss quadrature.
Parameters parameters
Data structure holding arbitrary parameters.
const MeshBase & get_mesh() const
void set_eigenproblem_type(EigenProblemType ept)
Sets the type of the current eigen problem.
Definition: eigen_system.C:78
unsigned int rank() const
Definition: parallel.h:724
const T_sys & get_system(const std::string &name) const
virtual void init()
Initialize all the systems.
virtual void add(const numeric_index_type i, const T value)=0
Adds value to each entry of the vector.
This class provides a specific system class.
Definition: eigen_system.h:50
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
void print_info(std::ostream &os=libMesh::out) const
Prints relevant information about the mesh.
Definition: mesh_base.C:448
void dof_indices(const Elem *const elem, std::vector< dof_id_type > &di) const
Fills the vector di with the global degree of freedom indices for the element.
Definition: dof_map.C:1917
static UniquePtr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.