condensed_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 #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 
25 
26 #include "libmesh/dof_map.h"
28 #include "libmesh/int_range.h"
30 #include "libmesh/numeric_vector.h"
31 #include "libmesh/parallel.h"
32 
33 namespace libMesh
34 {
35 
37  const std::string & name_in,
38  const unsigned int number_in)
39  : Parent(es, name_in, number_in),
40  condensed_matrix_A(SparseMatrix<Number>::build(es.comm())),
41  condensed_matrix_B(SparseMatrix<Number>::build(es.comm())),
42  condensed_dofs_initialized(false)
43 {
44 }
45 
46 void
47 CondensedEigenSystem::initialize_condensed_dofs(const std::set<dof_id_type> & global_dirichlet_dofs_set)
48 {
49  const DofMap & dof_map = this->get_dof_map();
50 
51  // First, put all unconstrained local dofs into non_dirichlet_dofs_set
52  std::set<dof_id_type> local_non_condensed_dofs_set;
53  for (dof_id_type i=this->get_dof_map().first_dof(); i<this->get_dof_map().end_dof(); i++)
54  if (!dof_map.is_constrained_dof(i))
55  local_non_condensed_dofs_set.insert(i);
56 
57  // Now erase the condensed dofs
58  for (const auto & dof : global_dirichlet_dofs_set)
59  if ((this->get_dof_map().first_dof() <= dof) && (dof < this->get_dof_map().end_dof()))
60  local_non_condensed_dofs_set.erase(dof);
61 
62  // Finally, move local_non_condensed_dofs_set over to a vector for convenience in solve()
63  this->local_non_condensed_dofs_vector.clear();
64 
65  for (const auto & dof : local_non_condensed_dofs_set)
66  this->local_non_condensed_dofs_vector.push_back(dof);
67 
69 }
70 
72 {
74  {
75  return this->n_dofs();
76  }
77  else
78  {
80  cast_int<dof_id_type>(local_non_condensed_dofs_vector.size());
81  this->comm().sum(n_global_non_condensed_dofs);
82 
84  }
85 }
86 
87 
89 {
90  LOG_SCOPE("solve()", "CondensedEigenSystem");
91 
92  // If we haven't initialized any condensed dofs,
93  // just use the default eigen_system
95  {
96  Parent::solve();
97  return;
98  }
99 
100  // A reference to the EquationSystems
101  EquationSystems & es = this->get_equation_systems();
102 
103  // check that necessary parameters have been set
104  libmesh_assert (es.parameters.have_parameter<unsigned int>("eigenpairs"));
105  libmesh_assert (es.parameters.have_parameter<unsigned int>("basis vectors"));
106 
107  if (this->assemble_before_solve)
108  {
109  // Assemble the linear system
110  this->assemble ();
111 
112  // And close the assembled matrices; using a non-closed matrix
113  // with create_submatrix() is deprecated.
114  matrix_A->close();
115  if (generalized())
116  matrix_B->close();
117  }
118 
119  // If we reach here, then there should be some non-condensed dofs
120  libmesh_assert(!local_non_condensed_dofs_vector.empty());
121 
122  // Now condense the matrices
123  matrix_A->create_submatrix(*condensed_matrix_A,
126 
127  if (generalized())
128  {
129  matrix_B->create_submatrix(*condensed_matrix_B,
132  }
133 
134 
135  // Get the tolerance for the solver and the maximum
136  // number of iterations. Here, we simply adopt the linear solver
137  // specific parameters.
138  const Real tol =
139  es.parameters.get<Real>("linear solver tolerance");
140 
141  const unsigned int maxits =
142  es.parameters.get<unsigned int>("linear solver maximum iterations");
143 
144  const unsigned int nev =
145  es.parameters.get<unsigned int>("eigenpairs");
146 
147  const unsigned int ncv =
148  es.parameters.get<unsigned int>("basis vectors");
149 
150  std::pair<unsigned int, unsigned int> solve_data;
151 
152  // call the solver depending on the type of eigenproblem
153  if (generalized())
154  {
155  //in case of a generalized eigenproblem
156  solve_data = eigen_solver->solve_generalized
157  (*condensed_matrix_A,*condensed_matrix_B, nev, ncv, tol, maxits);
158  }
159 
160  else
161  {
162  libmesh_assert (!matrix_B);
163 
164  //in case of a standard eigenproblem
165  solve_data = eigen_solver->solve_standard (*condensed_matrix_A, nev, ncv, tol, maxits);
166  }
167 
168  set_n_converged(solve_data.first);
169  set_n_iterations(solve_data.second);
170 }
171 
172 
173 
175 {
176  LOG_SCOPE("get_eigenpair()", "CondensedEigenSystem");
177 
178  // If we haven't initialized any condensed dofs,
179  // just use the default eigen_system
181  return Parent::get_eigenpair(i);
182 
183  // If we reach here, then there should be some non-condensed dofs
184  libmesh_assert(!local_non_condensed_dofs_vector.empty());
185 
186  // This function assumes that condensed_solve has just been called.
187  // If this is not the case, then we will trip an asset in get_eigenpair
188  std::unique_ptr<NumericVector<Number>> temp = NumericVector<Number>::build(this->comm());
189  const dof_id_type n_local =
190  cast_int<dof_id_type>(local_non_condensed_dofs_vector.size());
191  dof_id_type n = n_local;
192  this->comm().sum(n);
193 
194  temp->init (n, n_local, false, PARALLEL);
195 
196  std::pair<Real, Real> eval = eigen_solver->get_eigenpair (i, *temp);
197 
198  // Now map temp to solution. Loop over local entries of local_non_condensed_dofs_vector
199  this->solution->zero();
200  for (auto j : IntRange<dof_id_type>(0, n_local))
201  {
203  solution->set(index,(*temp)(temp->first_local_index()+j));
204  }
205 
206  solution->close();
207  this->update();
208 
209  return eval;
210 }
211 
212 
213 
214 
215 } // namespace libMesh
216 
217 
218 #endif // LIBMESH_HAVE_SLEPC
virtual std::pair< Real, Real > get_eigenpair(dof_id_type i)
Definition: eigen_system.C:256
bool generalized() const
Definition: eigen_system.h:149
Manages multiples systems of equations.
std::unique_ptr< SparseMatrix< Number > > condensed_matrix_A
const EquationSystems & get_equation_systems() const
Definition: system.h:712
CondensedEigenSystem(EquationSystems &es, const std::string &name_in, const unsigned int number_in)
const Parallel::Communicator & comm() const
dof_id_type n_dofs() const
Definition: system.C:150
virtual void assemble() override
Definition: eigen_system.C:247
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
dof_id_type n_global_non_condensed_dofs() const
void set_n_iterations(unsigned int its)
Definition: eigen_system.h:193
std::unique_ptr< NumericVector< Number > > solution
Definition: system.h:1523
bool is_constrained_dof(const dof_id_type dof) const
Definition: dof_map.h:1829
void set_n_converged(unsigned int nconv)
Definition: eigen_system.h:186
void initialize_condensed_dofs(const std::set< dof_id_type > &global_condensed_dofs_set=std::set< dof_id_type >())
static std::unique_ptr< NumericVector< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package())
virtual void solve() override
Definition: eigen_system.C:198
virtual void update()
Definition: system.C:408
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::unique_ptr< EigenSolver< Number > > eigen_solver
Definition: eigen_system.h:165
std::unique_ptr< SparseMatrix< Number > > condensed_matrix_B
virtual std::pair< Real, Real > get_eigenpair(dof_id_type i) override
std::vector< dof_id_type > local_non_condensed_dofs_vector
const T & get(const std::string &) const
Definition: parameters.h:425
std::unique_ptr< SparseMatrix< Number > > matrix_B
Definition: eigen_system.h:159
bool have_parameter(const std::string &) const
Definition: parameters.h:406
dof_id_type end_dof(const processor_id_type proc) const
Definition: dof_map.h:641
Manages consistently variables, degrees of freedom, and coefficient vectors for eigenvalue problems...
Definition: eigen_system.h:54
bool assemble_before_solve
Definition: system.h:1477
const DofMap & get_dof_map() const
Definition: system.h:2049
uint8_t dof_id_type
Definition: id_types.h:64