condensed_eigen_system.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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 
26 #include "libmesh/numeric_vector.h"
28 #include "libmesh/dof_map.h"
29 #include "libmesh/parallel.h"
30 
31 namespace libMesh
32 {
33 
35  const std::string& name,
36  const unsigned int number)
37  : Parent(es, name, number),
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 CondensedEigenSystem::initialize_condensed_dofs(std::set<unsigned int>& global_dirichlet_dofs_set)
45 {
46  // First, put all local dofs into non_dirichlet_dofs_set and
47  std::set<unsigned int> local_non_condensed_dofs_set;
48  for(unsigned int i=this->get_dof_map().first_dof(); i<this->get_dof_map().end_dof(); i++)
49  local_non_condensed_dofs_set.insert(i);
50 
51  // Now erase the condensed dofs
52  std::set<unsigned int>::iterator iter = global_dirichlet_dofs_set.begin();
53  std::set<unsigned int>::iterator iter_end = global_dirichlet_dofs_set.end();
54 
55  for ( ; iter != iter_end ; iter++)
56  {
57  unsigned int condensed_dof_index = *iter;
58  if ( (this->get_dof_map().first_dof() <= condensed_dof_index) &&
59  (condensed_dof_index < this->get_dof_map().end_dof()) )
60  {
61  local_non_condensed_dofs_set.erase(condensed_dof_index);
62  }
63  }
64 
65  // Finally, move local_non_condensed_dofs_set over to a vector for convenience in solve()
66  iter = local_non_condensed_dofs_set.begin();
67  iter_end = local_non_condensed_dofs_set.end();
68 
69  this->local_non_condensed_dofs_vector.clear();
70 
71  for ( ; iter != iter_end; ++iter)
72  {
73  unsigned int non_condensed_dof_index = *iter;
74 
75  this->local_non_condensed_dofs_vector.push_back(non_condensed_dof_index);
76  }
77 
79 }
80 
82 {
84  {
85  return this->n_dofs();
86  }
87  else
88  {
90  this->comm().sum(n_global_non_condensed_dofs);
91 
93  }
94 }
95 
96 
98 {
99  START_LOG("solve()", "CondensedEigenSystem");
100 
101  // If we haven't initialized any condensed dofs,
102  // just use the default eigen_system
104  {
105  STOP_LOG("solve()", "CondensedEigenSystem");
106  Parent::solve();
107  return;
108  }
109 
110  // A reference to the EquationSystems
112 
113  // check that necessary parameters have been set
114  libmesh_assert (es.parameters.have_parameter<unsigned int>("eigenpairs"));
115  libmesh_assert (es.parameters.have_parameter<unsigned int>("basis vectors"));
116 
117  if (this->assemble_before_solve)
118  // Assemble the linear system
119  this->assemble ();
120 
121  // If we reach here, then there should be some non-condensed dofs
123 
124  // Now condense the matrices
128 
129  if(generalized())
130  {
134  }
135 
136 
137  // Get the tolerance for the solver and the maximum
138  // number of iterations. Here, we simply adopt the linear solver
139  // specific parameters.
140  const Real tol =
141  es.parameters.get<Real>("linear solver tolerance");
142 
143  const unsigned int maxits =
144  es.parameters.get<unsigned int>("linear solver maximum iterations");
145 
146  const unsigned int nev =
147  es.parameters.get<unsigned int>("eigenpairs");
148 
149  const unsigned int ncv =
150  es.parameters.get<unsigned int>("basis vectors");
151 
152  std::pair<unsigned int, unsigned int> solve_data;
153 
154  // call the solver depending on the type of eigenproblem
155  if ( generalized() )
156  {
157  //in case of a generalized eigenproblem
158  solve_data = eigen_solver->solve_generalized
159  (*condensed_matrix_A,*condensed_matrix_B, nev, ncv, tol, maxits);
160  }
161 
162  else
163  {
165 
166  //in case of a standard eigenproblem
167  solve_data = eigen_solver->solve_standard (*condensed_matrix_A, nev, ncv, tol, maxits);
168  }
169 
170  set_n_converged(solve_data.first);
171  set_n_iterations(solve_data.second);
172 
173  STOP_LOG("solve()", "CondensedEigenSystem");
174 }
175 
176 
177 
178 std::pair<Real, Real> CondensedEigenSystem::get_eigenpair(unsigned int i)
179 {
180  START_LOG("get_eigenpair()", "CondensedEigenSystem");
181 
182  // If we haven't initialized any condensed dofs,
183  // just use the default eigen_system
185  {
186  STOP_LOG("get_eigenpair()", "CondensedEigenSystem");
187  return Parent::get_eigenpair(i);
188  }
189 
190  // If we reach here, then there should be some non-condensed dofs
192 
193  // This function assumes that condensed_solve has just been called.
194  // If this is not the case, then we will trip an asset in get_eigenpair
196  unsigned int n_local = local_non_condensed_dofs_vector.size();
197  unsigned int n = n_local;
198  this->comm().sum(n);
199 
200  temp->init (n, n_local, false, libMeshEnums::PARALLEL);
201 
202  std::pair<Real, Real> eval = eigen_solver->get_eigenpair (i, *temp);
203 
204  // Now map temp to solution. Loop over local entries of local_non_condensed_dofs_vector
205  this->solution->zero();
206  for (unsigned int j=0; j<local_non_condensed_dofs_vector.size(); j++)
207  {
208  unsigned int index = local_non_condensed_dofs_vector[j];
209  solution->set(index,(*temp)(temp->first_local_index()+j));
210  }
211 
212  solution->close();
213  this->update();
214 
215  STOP_LOG("get_eigenpair()", "CondensedEigenSystem");
216 
217  return eval;
218 }
219 
220 
221 
222 
223 } // namespace libMesh
224 
225 
226 #endif // LIBMESH_HAVE_SLEPC

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:04 UTC

Hosted By:
SourceForge.net Logo