libMesh
eigen_system.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 #include "libmesh/libmesh_config.h"
20 
21 // Currently, the EigenSystem should only be available
22 // if SLEPc support is enabled.
23 #if defined(LIBMESH_HAVE_SLEPC)
24 
25 // C++ includes
26 
27 // Local includes
28 #include "libmesh/eigen_system.h"
29 #include "libmesh/equation_systems.h"
30 #include "libmesh/sparse_matrix.h"
31 #include "libmesh/eigen_solver.h"
32 #include "libmesh/dof_map.h"
33 #include "libmesh/mesh_base.h"
34 
35 namespace libMesh
36 {
37 
38 
39 // ------------------------------------------------------------
40 // EigenSystem implementation
42  const std::string & name_in,
43  const unsigned int number_in
44  ) :
45  Parent (es, name_in, number_in),
46  eigen_solver (EigenSolver<Number>::build(es.comm())),
47  _n_converged_eigenpairs (0),
48  _n_iterations (0),
49  _is_generalized_eigenproblem (false),
50  _eigen_problem_type (NHEP)
51 {
52 }
53 
54 
55 
57 {
58  // clear data
59  this->clear();
60 }
61 
62 
63 
65 {
66  // Clear the parent data
67  Parent::clear();
68 
69  // Clean up the matrices
70  matrix_A.reset();
71  matrix_B.reset();
72 
73  // clear the solver
74  eigen_solver->clear();
75 }
76 
77 
79 {
80  _eigen_problem_type = ept;
81 
82  eigen_solver->set_eigenproblem_type(ept);
83 
84  // libMesh::out << "The Problem type is set to be: " << std::endl;
85 
86  switch (_eigen_problem_type)
87  {
88  case HEP:
89  // libMesh::out << "Hermitian" << std::endl;
90  break;
91 
92  case NHEP:
93  // libMesh::out << "Non-Hermitian" << std::endl;
94  break;
95 
96  case GHEP:
97  // libMesh::out << "Generalized Hermitian" << std::endl;
98  break;
99 
100  case GNHEP:
101  // libMesh::out << "Generalized Non-Hermitian" << std::endl;
102  break;
103 
104  case GHIEP:
105  // libMesh::out << "Generalized indefinite Hermitian" << std::endl;
106  break;
107 
108  default:
109  // libMesh::out << "not properly specified" << std::endl;
110  libmesh_error_msg("Unrecognized _eigen_problem_type = " << _eigen_problem_type);
111  }
112 }
113 
114 
115 
116 
118 {
119  // initialize parent data
121 
122  // define the type of eigenproblem
123  if (_eigen_problem_type == GNHEP ||
127 
128  // build the system matrix
129  matrix_A.reset(SparseMatrix<Number>::build(this->comm()).release());
130 
131  this->init_matrices();
132 }
133 
134 
135 
137 {
138  DofMap & dof_map = this->get_dof_map();
139 
140  dof_map.attach_matrix(*matrix_A);
141 
142  // build matrix_B only in case of a
143  // generalized problem
145  {
146  matrix_B.reset(SparseMatrix<Number>::build(this->comm()).release());
147  dof_map.attach_matrix(*matrix_B);
148  }
149 
150  dof_map.compute_sparsity(this->get_mesh());
151 
152  // initialize and zero system matrix
153  matrix_A->init();
154  matrix_A->zero();
155 
156  // eventually initialize and zero system matrix_B
158  {
159  matrix_B->init();
160  matrix_B->zero();
161  }
162 }
163 
164 
165 
167 {
168  // initialize parent data
169  Parent::reinit();
170 
171  // Clear the matrices
172  matrix_A->clear();
173 
175  matrix_B->clear();
176 
177  DofMap & dof_map = this->get_dof_map();
178 
179  // Clear the sparsity pattern
180  dof_map.clear_sparsity();
181 
182  // Compute the sparsity pattern for the current
183  // mesh and DOF distribution. This also updates
184  // both matrices, \p DofMap now knows them
185  dof_map.compute_sparsity(this->get_mesh());
186 
187  matrix_A->init();
188  matrix_A->zero();
189 
191  {
192  matrix_B->init();
193  matrix_B->zero();
194  }
195 }
196 
197 
198 
200 {
201 
202  // A reference to the EquationSystems
203  EquationSystems & es = this->get_equation_systems();
204 
205  // check that necessary parameters have been set
206  libmesh_assert (es.parameters.have_parameter<unsigned int>("eigenpairs"));
207  libmesh_assert (es.parameters.have_parameter<unsigned int>("basis vectors"));
208 
209  if (this->assemble_before_solve)
210  // Assemble the linear system
211  this->assemble ();
212 
213  // Get the tolerance for the solver and the maximum
214  // number of iterations. Here, we simply adopt the linear solver
215  // specific parameters.
216  const Real tol =
217  es.parameters.get<Real>("linear solver tolerance");
218 
219  const unsigned int maxits =
220  es.parameters.get<unsigned int>("linear solver maximum iterations");
221 
222  const unsigned int nev =
223  es.parameters.get<unsigned int>("eigenpairs");
224 
225  const unsigned int ncv =
226  es.parameters.get<unsigned int>("basis vectors");
227 
228  std::pair<unsigned int, unsigned int> solve_data;
229 
230  // call the solver depending on the type of eigenproblem
231 
232  // Generalized eigenproblem
234  solve_data = eigen_solver->solve_generalized (*matrix_A, *matrix_B, nev, ncv, tol, maxits);
235 
236  // Standard eigenproblem
237  else
238  {
240  solve_data = eigen_solver->solve_standard (*matrix_A, nev, ncv, tol, maxits);
241  }
242 
243  this->_n_converged_eigenpairs = solve_data.first;
244  this->_n_iterations = solve_data.second;
245 }
246 
247 
249 {
250 
251  // Assemble the linear system
252  Parent::assemble ();
253 
254 }
255 
256 
257 std::pair<Real, Real> EigenSystem::get_eigenpair (dof_id_type i)
258 {
259  // call the eigen_solver get_eigenpair method
260  return eigen_solver->get_eigenpair (i, *solution);
261 }
262 
263 } // namespace libMesh
264 
265 #endif // LIBMESH_HAVE_SLEPC
virtual std::pair< Real, Real > get_eigenpair(dof_id_type i)
Definition: eigen_system.C:257
This is the EquationSystems class.
virtual void clear()
Clear all the data structures associated with the system.
Definition: system.C:203
unsigned int _n_converged_eigenpairs
The number of converged eigenpairs.
Definition: eigen_system.h:198
virtual void reinit()
Reinitializes degrees of freedom and other required data on the current mesh.
Definition: system.C:385
EigenProblemType
Defines an enum for eigenproblem types.
virtual void assemble()
Prepares matrix and _dof_map for matrix assembly.
Definition: system.C:479
void attach_matrix(SparseMatrix< Number > &matrix)
Additional matrices may be handled with this DofMap.
Definition: dof_map.C:246
UniquePtr< SparseMatrix< Number > > matrix_B
A second system matrix for generalized eigenvalue problems.
Definition: eigen_system.h:155
virtual void init_data()
Initializes the data for the system.
Definition: system.C:260
const T & get(const std::string &) const
Definition: parameters.h:430
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
virtual void init_data() libmesh_override
Initializes the member data fields associated with the system, so that, e.g., assemble() may be used...
Definition: eigen_system.C:117
libmesh_assert(j)
unsigned int _n_iterations
The number of iterations of the eigen solver algorithm.
Definition: eigen_system.h:203
Generic sparse matrix.
Definition: dof_map.h:65
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:167
const MeshBase & get_mesh() const
Definition: system.h:2014
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
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
bool _is_generalized_eigenproblem
A boolean flag to indicate whether we are dealing with a generalized eigenvalue problem.
Definition: eigen_system.h:209
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
virtual ~EigenSystem()
Destructor.
Definition: eigen_system.C:56
virtual void clear() libmesh_override
Clear all the data structures associated with the system.
Definition: eigen_system.C:64
virtual void assemble() libmesh_override
Assembles the system matrix.
Definition: eigen_system.C:248
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
EigenSystem(EquationSystems &es, const std::string &name_in, const unsigned int number_in)
Constructor.
Definition: eigen_system.C:41
UniquePtr< EigenSolver< Number > > eigen_solver
The EigenSolver, defining which interface, i.e solver package to use.
Definition: eigen_system.h:161
virtual void init_matrices()
Initializes the matrices associated with the system.
Definition: eigen_system.C:136
const Parallel::Communicator & comm() const
EigenProblemType _eigen_problem_type
The type of the eigenvalue problem.
Definition: eigen_system.h:214
void clear_sparsity()
Clears the sparsity pattern.
Definition: dof_map.C:1772
const EquationSystems & get_equation_systems() const
Definition: system.h:712
Parameters parameters
Data structure holding arbitrary parameters.
void set_eigenproblem_type(EigenProblemType ept)
Sets the type of the current eigen problem.
Definition: eigen_system.C:78
bool have_parameter(const std::string &) const
Definition: parameters.h:411
virtual void reinit() libmesh_override
Reinitializes the member data fields associated with the system, so that, e.g., assemble() may be use...
Definition: eigen_system.C:166
void compute_sparsity(const MeshBase &)
Computes the sparsity pattern for the matrices corresponding to proc_id and sends that data to Linear...
Definition: dof_map.C:1735
bool assemble_before_solve
Flag which tells the system to whether or not to call the user assembly function during each call to ...
Definition: system.h:1477
uint8_t dof_id_type
Definition: id_types.h:64
This class provides an interface to solvers for eigenvalue problems.
Definition: eigen_solver.h:54