eigen_system.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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 // Local includes
26 #include "libmesh/eigen_system.h"
28 #include "libmesh/sparse_matrix.h"
29 #include "libmesh/eigen_solver.h"
30 #include "libmesh/dof_map.h"
31 #include "libmesh/mesh_base.h"
33 
34 namespace libMesh
35 {
36 
37 
38 // ------------------------------------------------------------
39 // EigenSystem implementation
41  const std::string & name_in,
42  const unsigned int number_in
43  ) :
44  Parent (es, name_in, number_in),
45  eigen_solver (EigenSolver<Number>::build(es.comm())),
46  _n_converged_eigenpairs (0),
47  _n_iterations (0),
48  _is_generalized_eigenproblem (false),
49  _eigen_problem_type (NHEP)
50 {
51 }
52 
53 
54 
56 {
57  // clear data
58  this->clear();
59 }
60 
61 
62 
64 {
65  // Clear the parent data
66  Parent::clear();
67 
68  // Clean up the matrices
69  matrix_A.reset();
70  matrix_B.reset();
71 
72  // clear the solver
73  eigen_solver->clear();
74 }
75 
76 
78 {
79  _eigen_problem_type = ept;
80 
81  eigen_solver->set_eigenproblem_type(ept);
82 
83  // libMesh::out << "The Problem type is set to be: " << std::endl;
84 
85  switch (_eigen_problem_type)
86  {
87  case HEP:
88  // libMesh::out << "Hermitian" << std::endl;
89  break;
90 
91  case NHEP:
92  // libMesh::out << "Non-Hermitian" << std::endl;
93  break;
94 
95  case GHEP:
96  // libMesh::out << "Generalized Hermitian" << std::endl;
97  break;
98 
99  case GNHEP:
100  // libMesh::out << "Generalized Non-Hermitian" << std::endl;
101  break;
102 
103  case GHIEP:
104  // libMesh::out << "Generalized indefinite Hermitian" << std::endl;
105  break;
106 
107  default:
108  // libMesh::out << "not properly specified" << std::endl;
109  libmesh_error_msg("Unrecognized _eigen_problem_type = " << _eigen_problem_type);
110  }
111 }
112 
113 
114 
115 
117 {
118  // initialize parent data
120 
121  // define the type of eigenproblem
122  if (_eigen_problem_type == GNHEP ||
126 
127  // build the system matrix
129 
130  this->init_matrices();
131 }
132 
133 
134 
136 {
137  DofMap & dof_map = this->get_dof_map();
138 
139  dof_map.attach_matrix(*matrix_A);
140 
141  // build matrix_B only in case of a
142  // generalized problem
144  {
146  dof_map.attach_matrix(*matrix_B);
147  }
148 
149  dof_map.compute_sparsity(this->get_mesh());
150 
151  // initialize and zero system matrix
152  matrix_A->init();
153  matrix_A->zero();
154 
155  // eventually initialize and zero system matrix_B
157  {
158  matrix_B->init();
159  matrix_B->zero();
160  }
161 }
162 
163 
164 
166 {
167  // initialize parent data
168  Parent::reinit();
169 
170  // Clear the matrices
171  matrix_A->clear();
172 
174  matrix_B->clear();
175 
176  DofMap & dof_map = this->get_dof_map();
177 
178  // Clear the sparsity pattern
179  dof_map.clear_sparsity();
180 
181  // Compute the sparsity pattern for the current
182  // mesh and DOF distribution. This also updates
183  // both matrices, \p DofMap now knows them
184  dof_map.compute_sparsity(this->get_mesh());
185 
186  matrix_A->init();
187  matrix_A->zero();
188 
190  {
191  matrix_B->init();
192  matrix_B->zero();
193  }
194 }
195 
196 
197 
199 {
200 
201  // A reference to the EquationSystems
202  EquationSystems & es = this->get_equation_systems();
203 
204  // check that necessary parameters have been set
205  libmesh_assert (es.parameters.have_parameter<unsigned int>("eigenpairs"));
206  libmesh_assert (es.parameters.have_parameter<unsigned int>("basis vectors"));
207 
208  if (this->assemble_before_solve)
209  // Assemble the linear system
210  this->assemble ();
211 
212  // Get the tolerance for the solver and the maximum
213  // number of iterations. Here, we simply adopt the linear solver
214  // specific parameters.
215  const Real tol =
216  es.parameters.get<Real>("linear solver tolerance");
217 
218  const unsigned int maxits =
219  es.parameters.get<unsigned int>("linear solver maximum iterations");
220 
221  const unsigned int nev =
222  es.parameters.get<unsigned int>("eigenpairs");
223 
224  const unsigned int ncv =
225  es.parameters.get<unsigned int>("basis vectors");
226 
227  std::pair<unsigned int, unsigned int> solve_data;
228 
229  // call the solver depending on the type of eigenproblem
230 
231  // Generalized eigenproblem
233  solve_data = eigen_solver->solve_generalized (*matrix_A, *matrix_B, nev, ncv, tol, maxits);
234 
235  // Standard eigenproblem
236  else
237  {
238  libmesh_assert (!matrix_B);
239  solve_data = eigen_solver->solve_standard (*matrix_A, nev, ncv, tol, maxits);
240  }
241 
242  this->_n_converged_eigenpairs = solve_data.first;
243  this->_n_iterations = solve_data.second;
244 }
245 
246 
248 {
249 
250  // Assemble the linear system
251  Parent::assemble ();
252 
253 }
254 
255 
256 std::pair<Real, Real> EigenSystem::get_eigenpair (dof_id_type i)
257 {
258  // call the eigen_solver get_eigenpair method
259  return eigen_solver->get_eigenpair (i, *solution);
260 }
261 
262 } // namespace libMesh
263 
264 #endif // LIBMESH_HAVE_SLEPC
static std::unique_ptr< SparseMatrix< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
virtual void reinit() override
Definition: eigen_system.C:165
virtual std::pair< Real, Real > get_eigenpair(dof_id_type i)
Definition: eigen_system.C:256
Manages multiples systems of equations.
virtual void clear()
Definition: system.C:205
unsigned int _n_converged_eigenpairs
Definition: eigen_system.h:202
virtual void reinit()
Definition: system.C:390
const EquationSystems & get_equation_systems() const
Definition: system.h:712
virtual void assemble()
Definition: system.C:462
void attach_matrix(SparseMatrix< Number > &matrix)
Definition: dof_map.C:262
virtual void init_data()
Definition: system.C:262
const Parallel::Communicator & comm() const
const MeshBase & get_mesh() const
Definition: system.h:2033
virtual void assemble() override
Definition: eigen_system.C:247
unsigned int _n_iterations
Definition: eigen_system.h:207
Manages the degrees of freedom (DOFs) in a simulation.
Definition: dof_map.h:176
std::unique_ptr< SparseMatrix< Number > > matrix_A
Definition: eigen_system.h:154
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:92
std::unique_ptr< NumericVector< Number > > solution
Definition: system.h:1523
virtual ~EigenSystem()
Definition: eigen_system.C:55
virtual void solve() override
Definition: eigen_system.C:198
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
EigenSystem(EquationSystems &es, const std::string &name_in, const unsigned int number_in)
Definition: eigen_system.C:40
std::unique_ptr< EigenSolver< Number > > eigen_solver
Definition: eigen_system.h:165
virtual void init_matrices()
Definition: eigen_system.C:135
EigenProblemType _eigen_problem_type
Definition: eigen_system.h:218
void clear_sparsity()
Definition: dof_map.C:1771
const T & get(const std::string &) const
Definition: parameters.h:425
void set_eigenproblem_type(EigenProblemType ept)
Definition: eigen_system.C:77
std::unique_ptr< SparseMatrix< Number > > matrix_B
Definition: eigen_system.h:159
bool have_parameter(const std::string &) const
Definition: parameters.h:406
virtual void clear() override
Definition: eigen_system.C:63
void compute_sparsity(const MeshBase &)
Definition: dof_map.C:1739
bool assemble_before_solve
Definition: system.h:1477
const DofMap & get_dof_map() const
Definition: system.h:2049
virtual void init_data() override
Definition: eigen_system.C:116
uint8_t dof_id_type
Definition: id_types.h:64
Base class which defines the interface for solving eigenproblems.
Definition: eigen_solver.h:68