libMesh
condensed_eigen_system.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2025 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 
26 #include "libmesh/dof_map.h"
27 #include "libmesh/equation_systems.h"
28 #include "libmesh/int_range.h"
29 #include "libmesh/libmesh_logging.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_dofs_initialized(false),
41  _have_condensed_dofs(false),
42  _create_submatrices_in_solve(true)
43 {
44 }
45 
47 
48 void
49 CondensedEigenSystem::initialize_condensed_dofs(const std::set<dof_id_type> & global_dirichlet_dofs_set)
50 {
51  const DofMap & dof_map = this->get_dof_map();
52  if (global_dirichlet_dofs_set.empty() && !dof_map.n_constrained_dofs())
53  {
54  _have_condensed_dofs = false;
56  return;
57  }
58 
59  // First, put all unconstrained local dofs into non_dirichlet_dofs_set
60  std::set<dof_id_type> local_non_condensed_dofs_set;
61  for (auto i : make_range(dof_map.first_dof(), dof_map.end_dof()))
62 #if LIBMESH_ENABLE_CONSTRAINTS
63  if (!dof_map.is_constrained_dof(i))
64 #endif
65  local_non_condensed_dofs_set.insert(i);
66 
67  // Now erase the condensed dofs
68  for (const auto dof : global_dirichlet_dofs_set)
69  if ((dof_map.first_dof() <= dof) && (dof < dof_map.end_dof()))
70  local_non_condensed_dofs_set.erase(dof);
71 
72  // Finally, move local_non_condensed_dofs_set over to a vector for convenience in solve()
73  this->local_non_condensed_dofs_vector.clear();
74 
75  for (const auto dof : local_non_condensed_dofs_set)
76  this->local_non_condensed_dofs_vector.push_back(dof);
77 
78  _have_condensed_dofs = true;
80 }
81 
82 void
84 {
87 }
88 
89 void
91 {
94 
96  return;
97 
98  const auto m = cast_int<numeric_index_type>(this->local_non_condensed_dofs_vector.size());
99  auto M = m;
100  _communicator.sum(M);
101  if (this->has_condensed_matrix_A())
102  {
103  this->get_condensed_matrix_A().init(M, M, m, m);
104  // A bit ludicrously MatCopy requires the matrix being copied to to be assembled
105  this->get_condensed_matrix_A().close();
106  }
107  if (this->has_condensed_matrix_B())
108  {
109  this->get_condensed_matrix_B().init(M, M, m, m);
110  this->get_condensed_matrix_B().close();
111  }
112  if (this->has_condensed_precond_matrix())
113  {
114  this->get_condensed_precond_matrix().init(M, M, m, m);
116  }
117 }
118 
120 {
122  return this->n_dofs();
123  else
124  {
126  cast_int<dof_id_type>(local_non_condensed_dofs_vector.size());
127  this->comm().sum(n_global_non_condensed_dofs);
128 
130  }
131 }
132 
133 #ifdef LIBMESH_ENABLE_DEPRECATED
134 void
136 {
139 }
140 #endif
141 
142 void
144 {
145  Parent::clear();
146  _condensed_matrix_A = nullptr;
147  _condensed_matrix_B = nullptr;
148  _condensed_precond_matrix = nullptr;
149 #ifdef LIBMESH_ENABLE_DEPRECATED
151 #endif
152 }
153 
154 void
156 {
158 
159  if (!this->use_shell_matrices())
160  {
161  if (!_condensed_matrix_A)
163  if (!_condensed_matrix_B)
165 
166  // When not using shell matrices we use A for P as well so we don't need to add P
167  }
168  // we *are* using shell matrices but we might not be using a shell matrix for P
169  else if (!this->use_shell_precond_matrix())
170  {
173  }
174 
175 #ifdef LIBMESH_ENABLE_DEPRECATED
177 #endif
178 }
179 
181 {
182  LOG_SCOPE("solve()", "CondensedEigenSystem");
183 
184  // If we haven't initialized any condensed dofs,
185  // just use the default eigen_system
187  {
188  Parent::solve();
189  return;
190  }
191 
192  // check that necessary parameters have been set
194  this->get_equation_systems().parameters.have_parameter<unsigned int>("eigenpairs"));
196  this->get_equation_systems().parameters.have_parameter<unsigned int>("basis vectors"));
197 
198  if (this->assemble_before_solve)
199  {
200  // Assemble the linear system
201  this->assemble();
202 
203  // And close the assembled matrices; using a non-closed matrix
204  // with create_submatrix() is deprecated.
205  if (matrix_A)
206  matrix_A->close();
207  if (generalized() && matrix_B)
208  matrix_B->close();
209  if (precond_matrix)
211  }
212 
213  // If we reach here, then there should be some non-condensed dofs
215 
217  {
218  if (matrix_A)
221  if (generalized() && matrix_B)
224  if (precond_matrix)
228  }
229 
230  solve_helper(
232 }
233 
235 {
236  LOG_SCOPE("get_eigenpair()", "CondensedEigenSystem");
237 
238  // If we haven't initialized any condensed dofs,
239  // just use the default eigen_system
241  return Parent::get_eigenpair(i);
242 
243  // If we reach here, then there should be some non-condensed dofs
245 
246  // This function assumes that condensed_solve has just been called.
247  // If this is not the case, then we will trip an asset in get_eigenpair
248  std::unique_ptr<NumericVector<Number>> temp = NumericVector<Number>::build(this->comm());
249  const dof_id_type n_local =
250  cast_int<dof_id_type>(local_non_condensed_dofs_vector.size());
251  dof_id_type n = n_local;
252  this->comm().sum(n);
253 
254  temp->init (n, n_local, false, PARALLEL);
255 
256  std::pair<Real, Real> eval = eigen_solver->get_eigenpair (i, *temp);
257 
258  // Now map temp to solution. Loop over local entries of local_non_condensed_dofs_vector
259  libmesh_assert(this->comm().verify(this->solution->closed()));
260  if (!this->solution->closed())
261  this->solution->close();
262  this->solution->zero();
263  for (auto j : make_range(n_local))
264  {
266  solution->set(index,(*temp)(temp->first_local_index()+j));
267  }
268 
269  // Enforcing constraints requires creating a ghosted version of the solution, which requires the
270  // solution be assembled
271  solution->close();
273  this->update();
274 
275  return eval;
276 }
277 
278 
279 
281 {
283  return *_condensed_matrix_A;
284 }
285 
287 {
289  return *_condensed_matrix_B;
290 }
291 
293 {
296 }
297 
298 void
300  NumericVector<Number> & super)
301 {
302  libmesh_assert_equal_to(sub.local_size(), local_non_condensed_dofs_vector.size());
303  libmesh_assert_equal_to(sub.local_size() + this->get_dof_map().n_local_constrained_dofs(),
304  super.local_size());
305  auto super_sub_view = super.get_subvector(local_non_condensed_dofs_vector);
306  *super_sub_view = sub;
307  super.restore_subvector(std::move(super_sub_view), local_non_condensed_dofs_vector);
308 }
309 
310 void
312  NumericVector<Number> & sub)
313 {
314  libmesh_assert_equal_to(sub.local_size(), local_non_condensed_dofs_vector.size());
315  libmesh_assert_equal_to(sub.local_size() + this->get_dof_map().n_local_constrained_dofs(),
316  super.local_size());
317  auto super_sub_view = super.get_subvector(local_non_condensed_dofs_vector);
318  sub = *super_sub_view;
319  super.restore_subvector(std::move(super_sub_view), local_non_condensed_dofs_vector);
320 }
321 
322 void
324  SparseMatrix<Number> & sub)
325 {
326  parallel_object_only();
327 
328  libmesh_assert_equal_to(sub.local_m(), local_non_condensed_dofs_vector.size());
329  libmesh_assert_equal_to(sub.local_m() + this->get_dof_map().n_local_constrained_dofs(),
330  super.local_m());
331  auto super_sub_view = SparseMatrix<Number>::build(super.comm());
332  super.create_submatrix(
334  sub = *super_sub_view;
335 }
336 
337 } // namespace libMesh
338 
339 
340 #endif // LIBMESH_HAVE_SLEPC
virtual void reinit() override
Reinitializes the member data fields associated with the system, so that, e.g., assemble() may be use...
Definition: eigen_system.C:189
virtual std::pair< Real, Real > get_eigenpair(dof_id_type i)
Definition: eigen_system.C:308
void solve_helper(SparseMatrix< Number > *const A, SparseMatrix< Number > *const B, SparseMatrix< Number > *const P)
Definition: eigen_system.C:215
bool generalized() const
Definition: eigen_system.C:325
virtual void create_submatrix(SparseMatrix< T > &submatrix, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols) const
This function creates a matrix called "submatrix" which is defined by the row and column indices give...
void initialize_condensed_matrices()
Initializes the condensed matrices.
dof_id_type end_dof(const processor_id_type proc) const
Definition: dof_map_base.h:191
This is the EquationSystems class.
bool _condensed_dofs_initialized
A private flag to indicate whether the condensed dofs have been initialized.
bool have_parameter(std::string_view) const
Definition: parameters.h:395
SparseMatrix< Number > * matrix_B
A second system matrix for generalized eigenvalue problems.
Definition: eigen_system.h:321
SparseMatrix< Number > & get_condensed_precond_matrix()
virtual numeric_index_type local_m() const
Get the number of rows owned by this process.
dof_id_type n_local_constrained_dofs() const
const EquationSystems & get_equation_systems() const
Definition: system.h:721
void sum(T &r) const
virtual void assemble()
Prepares matrix and _dof_map for matrix assembly.
Definition: system.C:566
CondensedEigenSystem(EquationSystems &es, const std::string &name_in, const unsigned int number_in)
Constructor.
SparseMatrix< Number > * condensed_matrix_B
A second (condensed) system matrix for generalized eigenvalue problems.
const Parallel::Communicator & comm() const
bool use_shell_precond_matrix() const
Definition: eigen_system.h:171
const Parallel::Communicator & _communicator
virtual void clear() override
Clear all the data structures associated with the system.
The libMesh namespace provides an interface to certain functionality in the library.
void copy_super_to_sub(NumericVector< Number > &super, NumericVector< Number > &sub)
Copy a logically super-vector into a sub-vector.
SparseMatrix< Number > * condensed_matrix_A
The (condensed) system matrix for standard eigenvalue problems.
Parameters parameters
Parameters for the system. If a parameter is not provided, it should be retrieved from the EquationSy...
Definition: system.h:1526
static std::unique_ptr< SparseMatrix< T > > build(const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package(), const MatrixBuildType matrix_build_type=MatrixBuildType::AUTOMATIC)
Builds a SparseMatrix<T> using the linear solver package specified by solver_package.
virtual void solve() override
Override to solve the condensed eigenproblem with the dofs in local_non_condensed_dofs_vector strippe...
dof_id_type n_dofs() const
Definition: system.C:121
This class handles the numbering of degrees of freedom on a mesh.
Definition: dof_map.h:179
SparseMatrix< Number > * precond_matrix
A preconditioning matrix.
Definition: eigen_system.h:345
SparseMatrix< Number > & get_condensed_matrix_A()
bool use_shell_matrices() const
Definition: eigen_system.h:161
dof_id_type n_global_non_condensed_dofs() const
void copy_sub_to_super(const NumericVector< Number > &sub, NumericVector< Number > &super)
Copy a logically sub-vector into a super-vector.
bool _have_condensed_dofs
Whether there are any condensed degrees of freedom.
std::unique_ptr< NumericVector< Number > > solution
Data structure to hold solution values.
Definition: system.h:1593
libmesh_assert(ctx)
bool is_constrained_dof(const dof_id_type dof) const
Definition: dof_map.h:2258
SparseMatrix< Number > * matrix_A
The system matrix for standard eigenvalue problems.
Definition: eigen_system.h:313
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 std::unique_ptr< NumericVector< T > > get_subvector(const std::vector< numeric_index_type > &)
Creates a view into this vector using the indices in rows.
virtual void solve() override
Assembles & solves the eigen system.
Definition: eigen_system.C:287
virtual void add_matrices() override
Adds the necessary matrices and shell matrices.
virtual void update()
Update the local values to reflect the solution on neighboring processors.
Definition: system.C:510
virtual void close()=0
Calls the SparseMatrix&#39;s internal assembly routines, ensuring that the values are consistent across p...
dof_id_type n_constrained_dofs() const
std::unique_ptr< SparseMatrix< Number > > _condensed_matrix_B
A second (condensed) system matrix for generalized eigenvalue problems.
std::unique_ptr< EigenSolver< Number > > eigen_solver
The EigenSolver, defining which interface, i.e solver package to use.
Definition: eigen_system.h:362
virtual numeric_index_type local_size() const =0
virtual void add_matrices() override
Adds the necessary matrices and shell matrices.
Definition: eigen_system.C:126
virtual std::pair< Real, Real > get_eigenpair(dof_id_type i) override
Override get_eigenpair() to retrieve the eigenpair for the condensed eigensolve.
std::vector< dof_id_type > local_non_condensed_dofs_vector
Vector storing the local dof indices that will not be condensed.
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:140
static std::unique_ptr< NumericVector< T > > build(const Parallel::Communicator &comm, SolverPackage solver_package=libMesh::default_solver_package(), ParallelType parallel_type=AUTOMATIC)
Builds a NumericVector on the processors in communicator comm using the linear solver package specifi...
bool _create_submatrices_in_solve
Denotes whether to create the condensed submatrices from the global matrices in the solve...
virtual void reinit() override
Reinitializes the member data fields associated with the system, so that, e.g., assemble() may be use...
SparseMatrix< Number > & get_condensed_matrix_B()
dof_id_type first_dof(const processor_id_type proc) const
Definition: dof_map_base.h:185
virtual void clear() override
Clear all the data structures associated with the system.
Definition: eigen_system.C:62
Manages consistently variables, degrees of freedom, and coefficient vectors for eigenvalue problems...
Definition: eigen_system.h:55
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:1547
const DofMap & get_dof_map() const
Definition: system.h:2374
std::unique_ptr< SparseMatrix< Number > > _condensed_precond_matrix
The condensed preconditioning matrix.
std::unique_ptr< SparseMatrix< Number > > _condensed_matrix_A
The (condensed) system matrix for standard eigenvalue problems.
uint8_t dof_id_type
Definition: id_types.h:67
virtual void init(const numeric_index_type m, const numeric_index_type n, const numeric_index_type m_l, const numeric_index_type n_l, const numeric_index_type nnz=30, const numeric_index_type noz=10, const numeric_index_type blocksize=1)=0
Initialize SparseMatrix with the specified sizes.
void enforce_constraints_exactly(const System &system, NumericVector< Number > *v=nullptr, bool homogeneous=false) const
Constrains the numeric vector v, which represents a solution defined on the mesh. ...
Definition: dof_map.h:2350
virtual void restore_subvector(std::unique_ptr< NumericVector< T >>, const std::vector< numeric_index_type > &)
Restores a view into this vector using the indices in rows.