libMesh
condensed_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 #include "libmesh/libmesh_config.h"
19 
20 // Currently, the EigenSystem should only be available
21 // if SLEPc support is enabled.
22 #if defined(LIBMESH_HAVE_SLEPC)
23 
24 #include "libmesh/condensed_eigen_system.h"
25 #include "libmesh/libmesh_logging.h"
26 #include "libmesh/numeric_vector.h"
27 #include "libmesh/equation_systems.h"
28 #include "libmesh/dof_map.h"
29 #include "libmesh/parallel.h"
30 
31 namespace libMesh
32 {
33 
35  const std::string & name_in,
36  const unsigned int number_in)
37  : Parent(es, name_in, number_in),
38  condensed_matrix_A(SparseMatrix<Number>::build(es.comm())),
39  condensed_matrix_B(SparseMatrix<Number>::build(es.comm())),
40  condensed_dofs_initialized(false)
41 {
42 }
43 
44 void
45 CondensedEigenSystem::initialize_condensed_dofs(const std::set<dof_id_type> & global_dirichlet_dofs_set)
46 {
47  const DofMap & dof_map = this->get_dof_map();
48 
49  // First, put all unconstrained local dofs into non_dirichlet_dofs_set
50  std::set<dof_id_type> local_non_condensed_dofs_set;
51  for (dof_id_type i=this->get_dof_map().first_dof(); i<this->get_dof_map().end_dof(); i++)
52  if (!dof_map.is_constrained_dof(i))
53  local_non_condensed_dofs_set.insert(i);
54 
55  // Now erase the condensed dofs
56  std::set<dof_id_type>::iterator iter = global_dirichlet_dofs_set.begin();
57  std::set<dof_id_type>::iterator iter_end = global_dirichlet_dofs_set.end();
58 
59  for ( ; iter != iter_end ; ++iter)
60  {
61  dof_id_type condensed_dof_index = *iter;
62  if ( (this->get_dof_map().first_dof() <= condensed_dof_index) &&
63  (condensed_dof_index < this->get_dof_map().end_dof()) )
64  {
65  local_non_condensed_dofs_set.erase(condensed_dof_index);
66  }
67  }
68 
69  // Finally, move local_non_condensed_dofs_set over to a vector for convenience in solve()
70  iter = local_non_condensed_dofs_set.begin();
71  iter_end = local_non_condensed_dofs_set.end();
72 
73  this->local_non_condensed_dofs_vector.clear();
74 
75  for ( ; iter != iter_end; ++iter)
76  {
77  dof_id_type non_condensed_dof_index = *iter;
78 
79  this->local_non_condensed_dofs_vector.push_back(non_condensed_dof_index);
80  }
81 
83 }
84 
86 {
88  {
89  return this->n_dofs();
90  }
91  else
92  {
94  this->comm().sum(n_global_non_condensed_dofs);
95 
97  }
98 }
99 
100 
102 {
103  LOG_SCOPE("solve()", "CondensedEigenSystem");
104 
105  // If we haven't initialized any condensed dofs,
106  // just use the default eigen_system
108  {
109  Parent::solve();
110  return;
111  }
112 
113  // A reference to the EquationSystems
114  EquationSystems & es = this->get_equation_systems();
115 
116  // check that necessary parameters have been set
117  libmesh_assert (es.parameters.have_parameter<unsigned int>("eigenpairs"));
118  libmesh_assert (es.parameters.have_parameter<unsigned int>("basis vectors"));
119 
120  if (this->assemble_before_solve)
121  {
122  // Assemble the linear system
123  this->assemble ();
124 
125  // And close the assembled matrices; using a non-closed matrix
126  // with create_submatrix() is deprecated.
127  matrix_A->close();
128  if (generalized())
129  matrix_B->close();
130  }
131 
132  // If we reach here, then there should be some non-condensed dofs
134 
135  // Now condense the matrices
136  matrix_A->create_submatrix(*condensed_matrix_A,
139 
140  if (generalized())
141  {
142  matrix_B->create_submatrix(*condensed_matrix_B,
145  }
146 
147 
148  // Get the tolerance for the solver and the maximum
149  // number of iterations. Here, we simply adopt the linear solver
150  // specific parameters.
151  const Real tol =
152  es.parameters.get<Real>("linear solver tolerance");
153 
154  const unsigned int maxits =
155  es.parameters.get<unsigned int>("linear solver maximum iterations");
156 
157  const unsigned int nev =
158  es.parameters.get<unsigned int>("eigenpairs");
159 
160  const unsigned int ncv =
161  es.parameters.get<unsigned int>("basis vectors");
162 
163  std::pair<unsigned int, unsigned int> solve_data;
164 
165  // call the solver depending on the type of eigenproblem
166  if (generalized())
167  {
168  //in case of a generalized eigenproblem
169  solve_data = eigen_solver->solve_generalized
170  (*condensed_matrix_A,*condensed_matrix_B, nev, ncv, tol, maxits);
171  }
172 
173  else
174  {
176 
177  //in case of a standard eigenproblem
178  solve_data = eigen_solver->solve_standard (*condensed_matrix_A, nev, ncv, tol, maxits);
179  }
180 
181  set_n_converged(solve_data.first);
182  set_n_iterations(solve_data.second);
183 }
184 
185 
186 
188 {
189  LOG_SCOPE("get_eigenpair()", "CondensedEigenSystem");
190 
191  // If we haven't initialized any condensed dofs,
192  // just use the default eigen_system
194  return Parent::get_eigenpair(i);
195 
196  // If we reach here, then there should be some non-condensed dofs
198 
199  // This function assumes that condensed_solve has just been called.
200  // If this is not the case, then we will trip an asset in get_eigenpair
203  dof_id_type n = n_local;
204  this->comm().sum(n);
205 
206  temp->init (n, n_local, false, PARALLEL);
207 
208  std::pair<Real, Real> eval = eigen_solver->get_eigenpair (i, *temp);
209 
210  // Now map temp to solution. Loop over local entries of local_non_condensed_dofs_vector
211  this->solution->zero();
212  for (std::size_t j=0; j<local_non_condensed_dofs_vector.size(); j++)
213  {
215  solution->set(index,(*temp)(temp->first_local_index()+j));
216  }
217 
218  solution->close();
219  this->update();
220 
221  return eval;
222 }
223 
224 
225 
226 
227 } // namespace libMesh
228 
229 
230 #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.
UniquePtr< SparseMatrix< Number > > condensed_matrix_B
A second (condensed) system matrix for generalized eigenvalue problems.
UniquePtr< SparseMatrix< Number > > matrix_B
A second system matrix for generalized eigenvalue problems.
Definition: eigen_system.h:155
virtual void solve() libmesh_override
Override to solve the condensed eigenproblem with the dofs in local_non_condensed_dofs_vector strippe...
CondensedEigenSystem(EquationSystems &es, const std::string &name_in, const unsigned int number_in)
Constructor.
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
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
static UniquePtr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
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
dof_id_type end_dof(const processor_id_type proc) const
Definition: dof_map.h:577
bool is_constrained_dof(const dof_id_type dof) const
Definition: dof_map.h:1734
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 set_n_iterations(unsigned int its)
Set the _n_iterations member, useful for subclasses of EigenSystem.
Definition: eigen_system.h:189
bool generalized() const
Definition: eigen_system.h:145
void set_n_converged(unsigned int nconv)
Set the _n_converged_eigenpairs member, useful for subclasses of EigenSystem.
Definition: eigen_system.h:182
UniquePtr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1523
void initialize_condensed_dofs(const std::set< dof_id_type > &global_condensed_dofs_set=std::set< dof_id_type >())
Loop over the dofs on each processor to initialize the list of non-condensed dofs.
virtual void assemble() libmesh_override
Assembles the system matrix.
Definition: eigen_system.C:248
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:425
UniquePtr< SparseMatrix< Number > > condensed_matrix_A
The (condensed) system matrix for standard eigenvalue problems.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
UniquePtr< EigenSolver< Number > > eigen_solver
The EigenSolver, defining which interface, i.e solver package to use.
Definition: eigen_system.h:161
const Parallel::Communicator & comm() const
std::vector< dof_id_type > local_non_condensed_dofs_vector
Vector storing the local dof indices that will not be condensed.
const EquationSystems & get_equation_systems() const
Definition: system.h:712
Parameters parameters
Data structure holding arbitrary parameters.
bool condensed_dofs_initialized
A private flag to indicate whether the condensed dofs have been initialized.
bool have_parameter(const std::string &) const
Definition: parameters.h:411
dof_id_type n_dofs() const
Definition: system.C:148
This class provides a specific system class.
Definition: eigen_system.h:50
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
virtual std::pair< Real, Real > get_eigenpair(dof_id_type i) libmesh_override
Override get_eigenpair() to retrieve the eigenpair for the condensed eigensolve.
void sum(T &r) const
Take a local variable and replace it with the sum of it&#39;s values on all processors.
dof_id_type n_global_non_condensed_dofs() const
uint8_t dof_id_type
Definition: id_types.h:64