www.mooseframework.org
NonlinearEigenSystem.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
3 //*
4 //* All rights reserved, see COPYRIGHT for full restrictions
5 //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 //*
7 //* Licensed under LGPL 2.1, please see LICENSE for details
8 //* https://www.gnu.org/licenses/lgpl-2.1.html
9 
10 #include "NonlinearEigenSystem.h"
11 
12 // MOOSE includes
13 #include "DirichletBC.h"
14 #include "EigenDirichletBC.h"
15 #include "ArrayDirichletBC.h"
16 #include "EigenArrayDirichletBC.h"
17 #include "EigenProblem.h"
18 #include "IntegratedBC.h"
19 #include "KernelBase.h"
20 #include "NodalBC.h"
21 #include "TimeIntegrator.h"
22 #include "SlepcSupport.h"
23 #include "DGKernelBase.h"
24 #include "ScalarKernelBase.h"
25 #include "MooseVariableScalar.h"
26 #include "ResidualObject.h"
27 
28 #include "libmesh/libmesh_config.h"
29 #include "libmesh/petsc_matrix.h"
30 #include "libmesh/sparse_matrix.h"
31 #include "libmesh/diagonal_matrix.h"
32 #include "libmesh/petsc_shell_matrix.h"
33 
34 #ifdef LIBMESH_HAVE_SLEPC
35 
36 namespace Moose
37 {
38 
39 void
40 assemble_matrix(EquationSystems & es, const std::string & system_name)
41 {
42  EigenProblem * p = es.parameters.get<EigenProblem *>("_eigen_problem");
43  EigenSystem & eigen_system = es.get_system<EigenSystem>(system_name);
44  NonlinearEigenSystem & eigen_nl = p->getNonlinearEigenSystem(/*nl_sys_num=*/0);
45 
46  // If this is a nonlinear eigenvalue problem,
47  // we do not need to assemble anything
49  {
50  // If you want an efficient eigensolver,
51  // please use PETSc 3.13 or newer.
52  // We need to do an unnecessary assembly,
53  // if you use PETSc that is older than 3.13.
54 #if PETSC_RELEASE_LESS_THAN(3, 13, 0)
55  if (eigen_system.has_matrix_B())
56  p->computeJacobianTag(*eigen_system.current_local_solution,
57  eigen_system.get_matrix_B(),
58  eigen_nl.eigenMatrixTag());
59 #endif
60  return;
61  }
62 
63 #if !PETSC_RELEASE_LESS_THAN(3, 13, 0)
64  // If we use shell matrices and do not use a shell preconditioning matrix,
65  // we only need to form a preconditioning matrix
66  if (eigen_system.use_shell_matrices() && !eigen_system.use_shell_precond_matrix())
67  {
68  p->computeJacobianTag(*eigen_system.current_local_solution,
69  eigen_system.get_precond_matrix(),
70  eigen_nl.precondMatrixTag());
71  return;
72  }
73 #endif
74  // If it is a linear generalized eigenvalue problem,
75  // we assemble A and B together
76  if (eigen_system.generalized())
77  {
78  p->computeJacobianAB(*eigen_system.current_local_solution,
79  eigen_system.get_matrix_A(),
80  eigen_system.get_matrix_B(),
81  eigen_nl.nonEigenMatrixTag(),
82  eigen_nl.eigenMatrixTag());
83 #if LIBMESH_HAVE_SLEPC
84  if (p->negativeSignEigenKernel())
85  MatScale(static_cast<PetscMatrix<Number> &>(eigen_system.get_matrix_B()).mat(), -1.0);
86 #endif
87  return;
88  }
89 
90  // If it is a linear eigenvalue problem, we assemble matrix A
91  {
92  p->computeJacobianTag(*eigen_system.current_local_solution,
93  eigen_system.get_matrix_A(),
94  eigen_nl.nonEigenMatrixTag());
95 
96  return;
97  }
98 }
99 }
100 
101 NonlinearEigenSystem::NonlinearEigenSystem(EigenProblem & eigen_problem, const std::string & name)
102  : NonlinearSystemBase(eigen_problem, eigen_problem.es().add_system<EigenSystem>(name), name),
103  _eigen_sys(eigen_problem.es().get_system<EigenSystem>(name)),
104  _eigen_problem(eigen_problem),
105  _solver_configuration(nullptr),
106  _n_eigen_pairs_required(eigen_problem.getNEigenPairsRequired()),
107  _work_rhs_vector_AX(addVector("work_rhs_vector_Ax", false, PARALLEL)),
108  _work_rhs_vector_BX(addVector("work_rhs_vector_Bx", false, PARALLEL)),
109  _precond_matrix_includes_eigen(false),
110  _preconditioner(nullptr)
111 {
112  sys().attach_assemble_function(Moose::assemble_matrix);
113 
114  SlepcEigenSolver<Number> * solver =
115  cast_ptr<SlepcEigenSolver<Number> *>(_eigen_sys.eigen_solver.get());
116 
117  if (!solver)
118  mooseError("A slepc eigen solver is required");
119 
120  // setup of our class @SlepcSolverConfiguration
121  _solver_configuration = std::make_unique<SlepcEigenSolverConfiguration>(eigen_problem, *solver);
122 
123  solver->set_solver_configuration(*_solver_configuration);
124 
125  _Ax_tag = eigen_problem.addVectorTag("Ax_tag");
126 
127  _Bx_tag = eigen_problem.addVectorTag("Eigen");
128 
129  _A_tag = eigen_problem.addMatrixTag("A_tag");
130 
131  _B_tag = eigen_problem.addMatrixTag("Eigen");
132 
133  // By default, _precond_tag and _A_tag will share the same
134  // objects. If we want to include eigen contributions to
135  // the preconditioning matrix, and then _precond_tag will
136  // point to part of "B" objects
137  _precond_tag = eigen_problem.addMatrixTag("Eigen_precond");
138 }
139 
140 void
142 {
143  // If it is an eigen dirichlet boundary condition, we should skip it because their
144  // contributions should be zero. If we do not skip it, preconditioning matrix will
145  // be singular because boundary elements are zero.
146  if (_precond_matrix_includes_eigen && !dynamic_cast<EigenDirichletBC *>(&object) &&
147  !dynamic_cast<EigenArrayDirichletBC *>(&object))
148  object.useMatrixTag(_precond_tag, {});
149 
150  auto & vtags = object.getVectorTags({});
151  auto & mtags = object.getMatrixTags({});
152  // If it is an eigen kernel, mark its variable as eigen
153  if (vtags.find(_Bx_tag) != vtags.end() || mtags.find(_B_tag) != mtags.end())
154  {
155  // Note: the object may be on the displaced system
156  auto sys = object.parameters().get<SystemBase *>("_sys");
157  auto vname = object.variable().name();
158  if (hasScalarVariable(vname))
159  sys->getScalarVariable(0, vname).eigen(true);
160  else
161  sys->getVariable(0, vname).eigen(true);
162  }
163 
164  // If this is not an eigen kernel
165  // If there is no vector eigen tag and no matrix eigen tag,
166  // then we consider this as noneigen kernel
167  if (vtags.find(_Bx_tag) == vtags.end() && mtags.find(_B_tag) == mtags.end())
168  {
169  // Noneigen Vector tag
170  object.useVectorTag(_Ax_tag, {});
171  // Noneigen Matrix tag
172  object.useMatrixTag(_A_tag, {});
173  // Noneigen Kernels
174  object.useMatrixTag(_precond_tag, {});
175  }
176  else
177  {
178  // Associate the eigen matrix tag and the vector tag
179  // if this is a eigen kernel
180  object.useMatrixTag(_B_tag, {});
181  object.useVectorTag(_Bx_tag, {});
182  }
183 }
184 
185 void
187 {
188  const bool presolve_succeeded = preSolve();
189  if (!presolve_succeeded)
190  return;
191 
192 // In DEBUG mode, Libmesh will check the residual automatically. This may cause
193 // an error because B does not need to assembly by default.
194 // When PETSc is older than 3.13, we always need to do an extra assembly,
195 // so we do not do "close" here
196 #if DEBUG && !PETSC_RELEASE_LESS_THAN(3, 13, 0)
197  if (sys().has_matrix_B())
198  sys().get_matrix_B().close();
199 #endif
200 
201  // We apply initial guess for only nonlinear solver
203  _eigen_sys.set_initial_space(solution());
204 
205  // Solve the transient problem if we have a time integrator; the
206  // steady problem if not.
207  if (_time_integrator)
208  {
209  _time_integrator->solve();
210  _time_integrator->postSolve();
211  }
212  else
213  system().solve();
214 
215  // store eigenvalues
216  unsigned int n_converged_eigenvalues = getNumConvergedEigenvalues();
217 
219 
220  if (_n_eigen_pairs_required < n_converged_eigenvalues)
221  n_converged_eigenvalues = _n_eigen_pairs_required;
222 
223  _eigen_values.resize(n_converged_eigenvalues);
224  for (unsigned int n = 0; n < n_converged_eigenvalues; n++)
226 
227  // Update the solution vector to the active eigenvector
228  if (n_converged_eigenvalues)
230 }
231 
232 void
234 {
235  // Tell libmesh not close matrices before solve
236  _eigen_sys.get_eigen_solver().set_close_matrix_before_solve(false);
237 
238  // Matrix A
239  if (_eigen_sys.has_matrix_A())
240  {
241  Mat mat = static_cast<PetscMatrix<Number> &>(_eigen_sys.get_matrix_A()).mat();
242 
244  }
245 
246  // Matrix B
247  if (_eigen_sys.has_matrix_B())
248  {
249  Mat mat = static_cast<PetscMatrix<Number> &>(_eigen_sys.get_matrix_B()).mat();
250 
252  }
253 
254  // Shell matrix A
255  if (_eigen_sys.has_shell_matrix_A())
256  {
257  Mat mat = static_cast<PetscShellMatrix<Number> &>(_eigen_sys.get_shell_matrix_A()).mat();
258 
259  // Attach callbacks for nonlinear eigenvalue solver
261 
262  // Set MatMult operations for shell
264  }
265 
266  // Shell matrix B
267  if (_eigen_sys.has_shell_matrix_B())
268  {
269  Mat mat = static_cast<PetscShellMatrix<Number> &>(_eigen_sys.get_shell_matrix_B()).mat();
270 
272 
273  // Set MatMult operations for shell
275  }
276 
277  // Preconditioning matrix
278  if (_eigen_sys.has_precond_matrix())
279  {
280  Mat mat = static_cast<PetscMatrix<Number> &>(_eigen_sys.get_precond_matrix()).mat();
281 
283  }
284 
285  // Shell preconditioning matrix
286  if (_eigen_sys.has_shell_precond_matrix())
287  {
288  Mat mat = static_cast<PetscShellMatrix<Number> &>(_eigen_sys.get_shell_precond_matrix()).mat();
289 
291  }
292 }
293 
294 void
296 {
297  mooseError("did not implement yet \n");
298 }
299 
300 void
302 {
303  mooseError("did not implement yet \n");
304 }
305 
306 bool
308 {
309  return _eigen_sys.get_n_converged();
310 }
311 
312 unsigned int
314 {
315  mooseError("did not implement yet \n");
316  return 0;
317 }
318 
321 {
322  return _work_rhs_vector_BX;
323 }
324 
327 {
328  return _work_rhs_vector_AX;
329 }
330 
333 {
334  return _work_rhs_vector_BX;
335 }
336 
339 {
340  mooseError("did not implement yet \n");
341  return NULL;
342 }
343 
344 SNES
346 {
347  EPS eps = getEPS();
348 
350  {
351  SNES snes = nullptr;
353  return snes;
354  }
355  else
356  mooseError("There is no SNES in linear eigen solver");
357 }
358 
359 EPS
361 {
362  SlepcEigenSolver<Number> * solver =
363  cast_ptr<SlepcEigenSolver<Number> *>(&(*_eigen_sys.eigen_solver));
364 
365  if (!solver)
366  mooseError("Unable to retrieve eigen solver");
367 
368  return solver->eps();
369 }
370 
371 void
373 {
375  {
376  const auto & nodal_bcs = _nodal_bcs.getActiveObjects();
377  for (const auto & nodal_bc : nodal_bcs)
378  {
379  // If this is a dirichlet boundary condition
380  auto nbc = std::dynamic_pointer_cast<DirichletBC>(nodal_bc);
381  // If this is a eigen Dirichlet boundary condition
382  auto eigen_nbc = std::dynamic_pointer_cast<EigenDirichletBC>(nodal_bc);
383  // ArrayDirichletBC
384  auto anbc = std::dynamic_pointer_cast<ArrayDirichletBC>(nodal_bc);
385  // EigenArrayDirichletBC
386  auto aeigen_nbc = std::dynamic_pointer_cast<EigenArrayDirichletBC>(nodal_bc);
387  // If it is a Dirichlet boundary condition, then value has to be zero
388  if (nbc && nbc->variable().eigen() && nbc->getParam<Real>("value"))
389  mooseError(
390  "Can't set an inhomogeneous Dirichlet boundary condition for eigenvalue problems.");
391  // If it is an array Dirichlet boundary condition, all values should be zero
392  else if (anbc)
393  {
394  auto & values = anbc->getParam<RealEigenVector>("values");
395  for (MooseIndex(values) i = 0; i < values.size(); i++)
396  {
397  if (values(i))
398  mooseError("Can't set an inhomogeneous array Dirichlet boundary condition for "
399  "eigenvalue problems.");
400  }
401  }
402  else if (!nbc && !eigen_nbc && !anbc && !aeigen_nbc)
403  mooseError(
404  "Invalid NodalBC for eigenvalue problems, please use homogeneous (array) Dirichlet.");
405  }
406  }
407 }
408 
409 std::pair<Real, Real>
411 {
412  unsigned int n_converged_eigenvalues = getNumConvergedEigenvalues();
413  if (n >= n_converged_eigenvalues)
414  mooseError(n, " not in [0, ", n_converged_eigenvalues, ")");
415  return _eigen_sys.get_eigenvalue(n);
416 }
417 
418 std::pair<Real, Real>
420 {
421  unsigned int n_converged_eigenvalues = getNumConvergedEigenvalues();
422 
423  if (n >= n_converged_eigenvalues)
424  mooseError(n, " not in [0, ", n_converged_eigenvalues, ")");
425 
426  return _eigen_sys.get_eigenpair(n);
427 }
428 
429 void
431 {
433 
434  // If we have a customized preconditioner,
435  // We need to let PETSc know that
436  if (_preconditioner)
437  {
439  // Mark this, and then we can setup correct petsc options
442  }
443 }
444 
445 void
447 {
448  // Let us do nothing at the current moment
449 }
450 
451 void
453 {
454  mooseError(
455  "NonlinearEigenSystem::residualAndJacobianTogether is not implemented. It might even be "
456  "nonsensical. If it is sensical and you want this capability, please contact a MOOSE "
457  "developer.");
458 }
459 
460 void
462 {
464 }
465 
466 void
468 {
470 }
471 
472 std::set<TagID>
474 {
476  tags.insert(eigenVectorTag());
477  tags.insert(nonEigenVectorTag());
478  return tags;
479 }
480 
481 std::set<TagID>
483 {
485  tags.insert(eigenMatrixTag());
486  tags.insert(nonEigenMatrixTag());
487  return tags;
488 }
489 
490 #else
491 
493  const std::string & /*name*/)
494  : libMesh::ParallelObject(eigen_problem)
495 {
496  mooseError("Need to install SLEPc to solve eigenvalue problems, please reconfigure libMesh\n");
497 }
498 
499 #endif /* LIBMESH_HAVE_SLEPC */
std::string name(const ElemQuality q)
Nonlinear eigenvalue system to be solved.
virtual void computeResidualTag(const NumericVector< Number > &soln, NumericVector< Number > &residual, TagID tag) override
Form a vector for all kernels and BCs with a given tag.
Definition: EigenProblem.C:282
std::unique_ptr< DiagonalMatrix< Number > > _scaling_matrix
A diagonal matrix used for computing scaling.
NumericVector< Number > & residualVectorAX()
int eps(unsigned int i, unsigned int j)
2D version
std::unique_ptr< SlepcEigenSolverConfiguration > _solver_configuration
SolverParams & solverParams()
Get the solver parameters.
unsigned int activeEigenvalueIndex() const
Which eigenvalue is active.
Definition: EigenProblem.h:179
NumericVector< Number > & solution()
Definition: SystemBase.h:176
PARALLEL
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
Definition: MooseError.h:284
virtual void postAddResidualObject(ResidualObject &object) override
Called after any ResidualObject-derived objects are added to the system.
std::vector< std::pair< R1, R2 > > get(const std::string &param1, const std::string &param2) const
Combine two vector parameters into a single vector of pairs.
TagID nonEigenVectorTag() const
Vector tag ID of left hand side.
virtual NonlinearSolver< Number > * nonlinearSolver() override
TagID eigenVectorTag() const
Vector tag ID of right hand side.
MooseObjectTagWarehouse< NodalBCBase > _nodal_bcs
std::set< TagID > defaultVectorTags() const override
Get the default vector tags associated with this system.
bool negativeSignEigenKernel() const
A flag indicates if a negative sign is used in eigen kernels.
Definition: EigenProblem.h:60
NumericVector< Number > & _work_rhs_vector_AX
std::unique_ptr< T_DEST, T_DELETER > dynamic_pointer_cast(std::unique_ptr< T_SRC, T_DELETER > &src)
These are reworked from https://stackoverflow.com/a/11003103.
bool _customized_pc_for_eigen
Definition: SolverParams.h:29
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
Base class for a system (of equations)
Definition: SystemBase.h:84
NonlinearEigenSystem(EigenProblem &problem, const std::string &name)
Preconditioner< Number > * _preconditioner
NumericVector< Number > & _work_rhs_vector_BX
virtual void stopSolve(const ExecFlagType &exec_flag) override
Quit the current solve as soon as possible.
bool isNonlinearEigenvalueSolver() const
Definition: EigenProblem.C:635
unsigned int getNumConvergedEigenvalues() const
Get the number of converged eigenvalues.
Nonlinear system to be solved.
virtual const std::string & name() const
Definition: SystemBase.C:1297
unsigned int getNEigenPairsRequired() const
Definition: EigenProblem.h:40
virtual EPS getEPS()
Retrieve EPS (SLEPc eigen solver)
const std::vector< std::shared_ptr< T > > & getActiveObjects(THREAD_ID tid=0) const
Retrieve complete vector to the active all/block/boundary restricted objects for a given thread...
void attachCallbacksToMat(EigenProblem &eigen_problem, Mat mat, bool eigen)
Attach call backs to mat.
Definition: SlepcSupport.C:982
NumericVector< Number > & residualVectorBX()
virtual void turnOffJacobian() override
Turn off the Jacobian (must be called before equation system initialization)
std::pair< Real, Real > getConvergedEigenpair(dof_id_type n) const
Return the Nth converged eigenvalue and copies the respective eigen vector to the solution vector...
Jacobian-Free Newton Krylov.
Definition: MooseTypes.h:758
virtual bool converged() override
Returns the convergence state.
Boundary condition of a Dirichlet type.
Definition: DirichletBC.h:19
virtual void attachPreconditioner(Preconditioner< Number > *preconditioner) override
Attach a customized preconditioner that requires physics knowledge.
Moose::SolveType _type
Definition: SolverParams.h:19
void assemble_matrix(EquationSystems &es, const std::string &system_name)
Set Dirichlet boundary condition for eigenvalue problems.
virtual void setupFiniteDifferencedPreconditioner() override
virtual std::set< TagID > defaultMatrixTags() const
Get the default matrix tags associted with this system.
Definition: SystemBase.h:287
void computeScalingResidual() override
Compute a "residual" for automatic scaling purposes.
virtual unsigned int getCurrentNonlinearIterationNumber() override
Returns the current nonlinear iteration number.
NonlinearEigenSystem & getNonlinearEigenSystem(const unsigned int nl_sys_num)
Definition: EigenProblem.h:302
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
This is the common base class for objects that give residual contributions.
PetscErrorCode mooseSlepcEPSGetSNES(EPS eps, SNES *snes)
Retrieve SNES from EPS.
Class for containing MooseEnum item information.
Definition: MooseEnumItem.h:18
bool hasActiveObjects(THREAD_ID tid=0) const
virtual System & system() override
Get the reference to the libMesh system.
Preconditioner< Number > * preconditioner() const
void residualAndJacobianTogether() override
Call this method if you want the residual and Jacobian to be computed simultaneously.
virtual std::set< TagID > defaultVectorTags() const
Get the default vector tags associated with this system.
Definition: SystemBase.h:280
virtual SNES getSNES() override
Retrieve snes from slepc eigen solver.
std::shared_ptr< TimeIntegrator > _time_integrator
Time integrator.
Definition: SystemBase.h:973
TagID eigenMatrixTag() const
Matrix tag ID of right hand side.
const NumericVector< Number > * _current_solution
solution vector from nonlinear solver
void computeScalingJacobian() override
Compute a "Jacobian" for automatic scaling purposes.
Boundary condition of a Dirichlet type for the eigen side.
const InputParameters & parameters() const
Get the parameters of the object.
void setOperationsForShellMat(EigenProblem &eigen_problem, Mat mat, bool eigen)
Set operations to shell mat.
Eigen::Matrix< Real, Eigen::Dynamic, 1 > RealEigenVector
Definition: MooseTypes.h:138
std::pair< Real, Real > getConvergedEigenvalue(dof_id_type n) const
Return the Nth converged eigenvalue.
TagID nonEigenMatrixTag() const
Matrix tag ID of left hand side.
MOOSE now contains C++17 code, so give a reasonable error message stating what the user can do to add...
bool preSolve()
Perform some steps to get ready for the solver.
Problem for solving eigenvalue problems.
Definition: EigenProblem.h:21
std::vector< std::pair< Real, Real > > _eigen_values
PETSC_EXTERN PetscErrorCode registerPCToPETSc()
Let PETSc know there is a preconditioner.
virtual NumericVector< Number > & RHS() override
unsigned int _n_eigen_pairs_required
virtual void solve() override
Solve the system (using libMesh magic)
virtual bool hasScalarVariable(const std::string &var_name) const
Definition: SystemBase.C:825
Boundary condition of a Dirichlet type.
std::set< TagID > defaultMatrixTags() const override
Get the default matrix tags associted with this system.
void checkIntegrity()
For eigenvalue problems (including standard and generalized), inhomogeneous (Dirichlet or Neumann) bo...
void computeJacobianAB(const NumericVector< Number > &soln, SparseMatrix< Number > &jacobianA, SparseMatrix< Number > &jacobianB, TagID tagA, TagID tagB)
Form two Jacobian matrices, where each is associated with one tag, through one element-loop.
Definition: EigenProblem.C:247
virtual void computeJacobianTag(const NumericVector< Number > &soln, SparseMatrix< Number > &jacobian, TagID tag) override
Form a Jacobian matrix for all kernels and BCs with a given tag.
Definition: EigenProblem.C:165
uint8_t dof_id_type