Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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"
27 : #include "libmesh/equation_systems.h"
28 : #include "libmesh/sparse_matrix.h"
29 : #include "libmesh/shell_matrix.h"
30 : #include "libmesh/eigen_solver.h"
31 : #include "libmesh/dof_map.h"
32 : #include "libmesh/mesh_base.h"
33 : #include "libmesh/enum_eigen_solver_type.h"
34 :
35 : namespace libMesh
36 : {
37 :
38 :
39 26 : EigenSystem::EigenSystem (EquationSystems & es,
40 : const std::string & name_in,
41 : const unsigned int number_in
42 26 : ) :
43 : Parent (es, name_in, number_in),
44 10 : matrix_A (nullptr),
45 10 : matrix_B (nullptr),
46 10 : precond_matrix (nullptr),
47 10 : eigen_solver (EigenSolver<Number>::build(es.comm())),
48 10 : _n_converged_eigenpairs (0),
49 10 : _n_iterations (0),
50 10 : _eigen_problem_type (NHEP),
51 10 : _use_shell_matrices (false),
52 26 : _use_shell_precond_matrix(false)
53 : {
54 26 : }
55 :
56 :
57 :
58 14 : EigenSystem::~EigenSystem () = default;
59 :
60 :
61 :
62 0 : void EigenSystem::clear ()
63 : {
64 : // Clear the parent data
65 0 : Parent::clear();
66 :
67 : // The SparseMatrices are contained in _matrices
68 : // and are cleared with the above call
69 0 : matrix_A = nullptr;
70 0 : matrix_B = nullptr;
71 0 : precond_matrix = nullptr;
72 :
73 : // Operators
74 0 : shell_matrix_A.reset();
75 0 : shell_matrix_B.reset();
76 :
77 : // Preconditioning matrices
78 0 : shell_precond_matrix.reset();
79 :
80 : // clear the solver
81 0 : eigen_solver->clear();
82 0 : }
83 :
84 :
85 20 : void EigenSystem::set_eigenproblem_type (EigenProblemType ept)
86 : {
87 20 : if (!can_add_matrices())
88 0 : libmesh_error_msg("ERROR: Cannot change eigen problem type after system initialization");
89 :
90 20 : _eigen_problem_type = ept;
91 :
92 6 : eigen_solver->set_eigenproblem_type(ept);
93 :
94 : // libMesh::out << "The Problem type is set to be: " << std::endl;
95 :
96 20 : switch (_eigen_problem_type)
97 : {
98 0 : case HEP:
99 : // libMesh::out << "Hermitian" << std::endl;
100 0 : break;
101 :
102 0 : case NHEP:
103 : // libMesh::out << "Non-Hermitian" << std::endl;
104 0 : break;
105 :
106 6 : case GHEP:
107 : // libMesh::out << "Generalized Hermitian" << std::endl;
108 6 : break;
109 :
110 0 : case GNHEP:
111 : // libMesh::out << "Generalized Non-Hermitian" << std::endl;
112 0 : break;
113 :
114 0 : case GHIEP:
115 : // libMesh::out << "Generalized indefinite Hermitian" << std::endl;
116 0 : break;
117 :
118 0 : default:
119 : // libMesh::out << "not properly specified" << std::endl;
120 0 : libmesh_error_msg("Unrecognized _eigen_problem_type = " << _eigen_problem_type);
121 : }
122 20 : }
123 :
124 :
125 :
126 26 : void EigenSystem::add_matrices ()
127 : {
128 26 : Parent::add_matrices();
129 :
130 26 : if (_use_shell_matrices)
131 : {
132 2 : if (!shell_matrix_A)
133 4 : shell_matrix_A = ShellMatrix<Number>::build(this->comm());
134 :
135 2 : if (generalized() && !shell_matrix_B)
136 4 : shell_matrix_B = ShellMatrix<Number>::build(this->comm());
137 :
138 2 : if (_use_shell_precond_matrix)
139 : {
140 0 : if (!shell_precond_matrix)
141 0 : shell_precond_matrix = ShellMatrix<Number>::build(this->comm());
142 : }
143 2 : else if (!precond_matrix)
144 2 : precond_matrix = &(this->add_matrix("Eigen Preconditioner"));
145 : }
146 : else
147 : {
148 24 : if (!matrix_A)
149 24 : matrix_A = &(this->add_matrix("Eigen Matrix A"));
150 :
151 24 : if (generalized() && !matrix_B)
152 18 : matrix_B = &(this->add_matrix("Eigen Matrix B"));
153 : }
154 26 : }
155 :
156 :
157 :
158 26 : void EigenSystem::init_matrices ()
159 : {
160 26 : Parent::init_matrices();
161 :
162 26 : const bool condense_constraints = condense_constrained_dofs();
163 26 : if (shell_matrix_A)
164 : {
165 0 : shell_matrix_A->attach_dof_map(this->get_dof_map());
166 2 : if (condense_constraints)
167 0 : shell_matrix_A->omit_constrained_dofs();
168 2 : shell_matrix_A->init();
169 : }
170 :
171 26 : if (shell_matrix_B)
172 : {
173 0 : shell_matrix_B->attach_dof_map(this->get_dof_map());
174 2 : if (condense_constraints)
175 0 : shell_matrix_B->omit_constrained_dofs();
176 2 : shell_matrix_B->init();
177 : }
178 :
179 26 : if (shell_precond_matrix)
180 : {
181 0 : shell_precond_matrix->attach_dof_map(this->get_dof_map());
182 0 : if (condense_constraints)
183 0 : shell_precond_matrix->omit_constrained_dofs();
184 0 : shell_precond_matrix->init();
185 : }
186 26 : }
187 :
188 :
189 0 : void EigenSystem::reinit ()
190 : {
191 : // initialize parent data
192 : // this calls reinit on matrix_A, matrix_B, and precond_matrix (if any)
193 0 : Parent::reinit();
194 :
195 0 : if (shell_matrix_A)
196 : {
197 0 : shell_matrix_A->clear();
198 0 : shell_matrix_A->init();
199 : }
200 :
201 0 : if (shell_matrix_B)
202 : {
203 0 : shell_matrix_B->clear();
204 0 : shell_matrix_B->init();
205 : }
206 :
207 0 : if (shell_precond_matrix)
208 : {
209 0 : shell_precond_matrix->clear();
210 0 : shell_precond_matrix->init();
211 : }
212 0 : }
213 :
214 : void
215 26 : EigenSystem::solve_helper(SparseMatrix<Number> * const A,
216 : SparseMatrix<Number> * const B,
217 : SparseMatrix<Number> * const P)
218 : {
219 : // A reference to the EquationSystems
220 16 : EquationSystems & es = this->get_equation_systems();
221 :
222 : // Get the tolerance for the solver and the maximum
223 : // number of iterations. Here, we simply adopt the linear solver
224 : // specific parameters.
225 44 : const double tol = parameters.have_parameter<Real>("linear solver tolerance") ?
226 0 : double(parameters.get<Real>("linear solver tolerance")) :
227 44 : double(es.parameters.get<Real>("linear solver tolerance"));
228 :
229 44 : const unsigned int maxits = parameters.have_parameter<unsigned int>("linear solver maximum iterations") ?
230 0 : parameters.get<unsigned int>("linear solver maximum iterations") :
231 44 : es.parameters.get<unsigned int>("linear solver maximum iterations");
232 :
233 44 : const unsigned int nev = parameters.have_parameter<unsigned int>("eigenpairs") ?
234 0 : parameters.get<unsigned int>("eigenpairs") :
235 44 : es.parameters.get<unsigned int>("eigenpairs");
236 :
237 44 : const unsigned int ncv = parameters.have_parameter<unsigned int>("basis vectors") ?
238 0 : parameters.get<unsigned int>("basis vectors") :
239 26 : es.parameters.get<unsigned int>("basis vectors");
240 :
241 8 : std::pair<unsigned int, unsigned int> solve_data;
242 :
243 : // call the solver depending on the type of eigenproblem
244 :
245 26 : if (_use_shell_matrices)
246 : {
247 : // Generalized eigenproblem
248 2 : if (generalized())
249 : // Shell preconditioning matrix
250 2 : if (_use_shell_precond_matrix)
251 0 : solve_data = eigen_solver->solve_generalized(
252 0 : *shell_matrix_A, *shell_matrix_B, *shell_precond_matrix, nev, ncv, tol, maxits);
253 : else
254 2 : solve_data = eigen_solver->solve_generalized(
255 0 : *shell_matrix_A, *shell_matrix_B, *P, nev, ncv, tol, maxits);
256 :
257 : // Standard eigenproblem
258 : else
259 : {
260 0 : libmesh_assert(!shell_matrix_B);
261 : // Shell preconditioning matrix
262 0 : if (_use_shell_precond_matrix)
263 0 : solve_data = eigen_solver->solve_standard(
264 0 : *shell_matrix_A, *shell_precond_matrix, nev, ncv, tol, maxits);
265 : else
266 0 : solve_data = eigen_solver->solve_standard(*shell_matrix_A, *P, nev, ncv, tol, maxits);
267 : }
268 : }
269 : else
270 : {
271 : // Generalized eigenproblem
272 24 : if (generalized())
273 18 : solve_data = eigen_solver->solve_generalized(*A, *B, nev, ncv, tol, maxits);
274 :
275 : // Standard eigenproblem
276 : else
277 : {
278 2 : libmesh_assert(!matrix_B);
279 6 : solve_data = eigen_solver->solve_standard(*A, nev, ncv, tol, maxits);
280 : }
281 : }
282 :
283 26 : this->_n_converged_eigenpairs = solve_data.first;
284 26 : this->_n_iterations = solve_data.second;
285 26 : }
286 :
287 12 : void EigenSystem::solve ()
288 : {
289 12 : if (get_dof_map().n_constrained_dofs())
290 : libmesh_warning("EigenSystem does not have first-class support for constrained degrees of "
291 : "freedom. You may see spurious effects of the constrained degrees of freedom "
292 : "in a given eigenvector. If you wish to perform a reliable solve on a system "
293 : "with constraints, please use the CondensedEigenSystem class instead");
294 :
295 : // check that necessary parameters have been set
296 4 : libmesh_assert(
297 : this->get_equation_systems().parameters.have_parameter<unsigned int>("eigenpairs"));
298 4 : libmesh_assert(
299 : this->get_equation_systems().parameters.have_parameter<unsigned int>("basis vectors"));
300 :
301 12 : if (this->assemble_before_solve)
302 : // Assemble the linear system
303 12 : this->assemble ();
304 :
305 12 : solve_helper(matrix_A, matrix_B, precond_matrix);
306 12 : }
307 :
308 12 : std::pair<Real, Real> EigenSystem::get_eigenpair (dof_id_type i)
309 : {
310 : // call the eigen_solver get_eigenpair method
311 16 : return eigen_solver->get_eigenpair (i, *solution);
312 : }
313 :
314 0 : std::pair<Real, Real> EigenSystem::get_eigenvalue (dof_id_type i)
315 : {
316 0 : return eigen_solver->get_eigenvalue (i);
317 : }
318 :
319 0 : void EigenSystem::set_initial_space (NumericVector<Number> & initial_space_in)
320 : {
321 0 : eigen_solver->set_initial_space (initial_space_in);
322 0 : }
323 :
324 :
325 82 : bool EigenSystem::generalized () const
326 : {
327 112 : return _eigen_problem_type == GNHEP ||
328 86 : _eigen_problem_type == GHEP ||
329 86 : _eigen_problem_type == GHIEP;
330 : }
331 :
332 :
333 0 : const SparseMatrix<Number> & EigenSystem::get_matrix_A() const
334 : {
335 0 : libmesh_assert(matrix_A);
336 0 : libmesh_assert(!_use_shell_matrices);
337 0 : libmesh_assert(has_matrix_A());
338 0 : return *matrix_A;
339 : }
340 :
341 :
342 24 : SparseMatrix<Number> & EigenSystem::get_matrix_A()
343 : {
344 8 : libmesh_assert(matrix_A);
345 8 : libmesh_assert(!_use_shell_matrices);
346 8 : libmesh_assert(has_matrix_A());
347 24 : return *matrix_A;
348 : }
349 :
350 :
351 0 : const SparseMatrix<Number> & EigenSystem::get_matrix_B() const
352 : {
353 0 : libmesh_assert(matrix_B);
354 0 : libmesh_assert(!_use_shell_matrices);
355 0 : libmesh_assert(generalized());
356 0 : libmesh_assert(has_matrix_B());
357 0 : return *matrix_B;
358 : }
359 :
360 :
361 18 : SparseMatrix<Number> & EigenSystem::get_matrix_B()
362 : {
363 6 : libmesh_assert(matrix_B);
364 6 : libmesh_assert(!_use_shell_matrices);
365 6 : libmesh_assert(generalized());
366 6 : libmesh_assert(has_matrix_B());
367 18 : return *matrix_B;
368 : }
369 :
370 :
371 0 : const SparseMatrix<Number> & EigenSystem::get_precond_matrix() const
372 : {
373 0 : libmesh_assert(precond_matrix);
374 0 : libmesh_assert(_use_shell_matrices);
375 0 : libmesh_assert(!_use_shell_precond_matrix);
376 0 : libmesh_assert(has_precond_matrix());
377 0 : return *precond_matrix;
378 : }
379 :
380 :
381 8 : SparseMatrix<Number> & EigenSystem::get_precond_matrix()
382 : {
383 0 : libmesh_assert(precond_matrix);
384 0 : libmesh_assert(_use_shell_matrices);
385 0 : libmesh_assert(!_use_shell_precond_matrix);
386 0 : libmesh_assert(has_precond_matrix());
387 8 : return *precond_matrix;
388 : }
389 :
390 :
391 0 : const ShellMatrix<Number> & EigenSystem::get_shell_matrix_A() const
392 : {
393 0 : libmesh_assert(_use_shell_matrices);
394 0 : libmesh_assert(has_shell_matrix_A());
395 0 : return *shell_matrix_A;
396 : }
397 :
398 :
399 2 : ShellMatrix<Number> & EigenSystem::get_shell_matrix_A()
400 : {
401 0 : libmesh_assert(_use_shell_matrices);
402 0 : libmesh_assert(has_shell_matrix_A());
403 2 : return *shell_matrix_A;
404 : }
405 :
406 :
407 0 : const ShellMatrix<Number> & EigenSystem::get_shell_matrix_B() const
408 : {
409 0 : libmesh_assert(_use_shell_matrices);
410 0 : libmesh_assert(generalized());
411 0 : libmesh_assert(has_shell_matrix_B());
412 0 : return *shell_matrix_B;
413 : }
414 :
415 :
416 2 : ShellMatrix<Number> & EigenSystem::get_shell_matrix_B()
417 : {
418 0 : libmesh_assert(_use_shell_matrices);
419 0 : libmesh_assert(generalized());
420 0 : libmesh_assert(has_shell_matrix_B());
421 2 : return *shell_matrix_B;
422 : }
423 :
424 :
425 0 : const ShellMatrix<Number> & EigenSystem::get_shell_precond_matrix() const
426 : {
427 0 : libmesh_assert(_use_shell_matrices);
428 0 : libmesh_assert(_use_shell_precond_matrix);
429 0 : libmesh_assert(has_shell_precond_matrix());
430 0 : return *shell_precond_matrix;
431 : }
432 :
433 :
434 0 : ShellMatrix<Number> & EigenSystem::get_shell_precond_matrix()
435 : {
436 0 : libmesh_assert(_use_shell_matrices);
437 0 : libmesh_assert(_use_shell_precond_matrix);
438 0 : libmesh_assert(has_shell_precond_matrix());
439 0 : return *shell_precond_matrix;
440 : }
441 :
442 :
443 8 : bool EigenSystem::has_matrix_A() const
444 : {
445 8 : libmesh_assert_equal_to(request_matrix("Eigen Matrix A"), matrix_A);
446 8 : return matrix_A;
447 : }
448 :
449 :
450 6 : bool EigenSystem::has_matrix_B() const
451 : {
452 6 : libmesh_assert_equal_to(request_matrix("Eigen Matrix B"), matrix_B);
453 6 : return matrix_B;
454 : }
455 :
456 :
457 0 : bool EigenSystem::has_precond_matrix() const
458 : {
459 0 : libmesh_assert_equal_to(request_matrix("Eigen Preconditioner"), precond_matrix);
460 0 : return precond_matrix;
461 : }
462 :
463 :
464 0 : bool EigenSystem::has_shell_matrix_A() const
465 : {
466 0 : return shell_matrix_A.get();
467 : }
468 :
469 :
470 0 : bool EigenSystem::has_shell_matrix_B() const
471 : {
472 0 : return shell_matrix_B.get();
473 : }
474 :
475 :
476 0 : bool EigenSystem::has_shell_precond_matrix() const
477 : {
478 0 : return shell_precond_matrix.get();
479 : }
480 :
481 :
482 0 : const EigenSolver<Number> & EigenSystem::get_eigen_solver() const
483 : {
484 0 : libmesh_assert(eigen_solver);
485 0 : return *eigen_solver;
486 : }
487 :
488 :
489 28 : EigenSolver<Number> & EigenSystem::get_eigen_solver()
490 : {
491 8 : libmesh_assert(eigen_solver);
492 28 : return *eigen_solver;
493 : }
494 :
495 :
496 : } // namespace libMesh
497 :
498 : #endif // LIBMESH_HAVE_SLEPC
|
