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" 34 #ifdef LIBMESH_HAVE_SLEPC 43 EigenSystem & eigen_system = es.get_system<EigenSystem>(system_name);
54 #if PETSC_RELEASE_LESS_THAN(3, 13, 0) 55 if (eigen_system.has_matrix_B())
57 eigen_system.get_matrix_B(),
58 eigen_nl.eigenMatrixTag());
63 #if !PETSC_RELEASE_LESS_THAN(3, 13, 0) 66 if (eigen_system.use_shell_matrices() && !eigen_system.use_shell_precond_matrix())
69 eigen_system.get_precond_matrix(),
70 eigen_nl.precondMatrixTag());
76 if (eigen_system.generalized())
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 85 MatScale(
static_cast<PetscMatrix<Number> &
>(eigen_system.get_matrix_B()).mat(), -1.0);
93 eigen_system.get_matrix_A(),
94 eigen_nl.nonEigenMatrixTag());
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)
114 SlepcEigenSolver<Number> * solver =
115 cast_ptr<SlepcEigenSolver<Number> *>(
_eigen_sys.eigen_solver.get());
118 mooseError(
"A slepc eigen solver is required");
125 _Ax_tag = eigen_problem.addVectorTag(
"Ax_tag");
127 _Bx_tag = eigen_problem.addVectorTag(
"Eigen");
129 _A_tag = eigen_problem.addMatrixTag(
"A_tag");
131 _B_tag = eigen_problem.addMatrixTag(
"Eigen");
137 _precond_tag = eigen_problem.addMatrixTag(
"Eigen_precond");
147 !dynamic_cast<EigenArrayDirichletBC *>(&
object))
150 auto & vtags =
object.getVectorTags({});
151 auto & mtags =
object.getMatrixTags({});
153 if (vtags.find(
_Bx_tag) != vtags.end() || mtags.find(
_B_tag) != mtags.end())
157 auto vname =
object.variable().
name();
159 sys->getScalarVariable(0, vname).eigen(
true);
161 sys->getVariable(0, vname).eigen(
true);
167 if (vtags.find(
_Bx_tag) == vtags.end() && mtags.find(
_B_tag) == mtags.end())
170 object.useVectorTag(
_Ax_tag, {});
172 object.useMatrixTag(
_A_tag, {});
180 object.useMatrixTag(
_B_tag, {});
181 object.useVectorTag(
_Bx_tag, {});
188 const bool presolve_succeeded =
preSolve();
189 if (!presolve_succeeded)
196 #if DEBUG && !PETSC_RELEASE_LESS_THAN(3, 13, 0) 197 if (
sys().has_matrix_B())
198 sys().get_matrix_B().close();
224 for (
unsigned int n = 0; n < n_converged_eigenvalues; n++)
228 if (n_converged_eigenvalues)
236 _eigen_sys.get_eigen_solver().set_close_matrix_before_solve(
false);
241 Mat mat =
static_cast<PetscMatrix<Number> &
>(
_eigen_sys.get_matrix_A()).mat();
249 Mat mat =
static_cast<PetscMatrix<Number> &
>(
_eigen_sys.get_matrix_B()).mat();
257 Mat mat =
static_cast<PetscShellMatrix<Number> &
>(
_eigen_sys.get_shell_matrix_A()).mat();
269 Mat mat =
static_cast<PetscShellMatrix<Number> &
>(
_eigen_sys.get_shell_matrix_B()).mat();
280 Mat mat =
static_cast<PetscMatrix<Number> &
>(
_eigen_sys.get_precond_matrix()).mat();
288 Mat mat =
static_cast<PetscShellMatrix<Number> &
>(
_eigen_sys.get_shell_precond_matrix()).mat();
356 mooseError(
"There is no SNES in linear eigen solver");
362 SlepcEigenSolver<Number> * solver =
363 cast_ptr<SlepcEigenSolver<Number> *>(&(*
_eigen_sys.eigen_solver));
366 mooseError(
"Unable to retrieve eigen solver");
368 return solver->eps();
377 for (
const auto & nodal_bc : nodal_bcs)
388 if (nbc && nbc->variable().eigen() && nbc->getParam<
Real>(
"value"))
390 "Can't set an inhomogeneous Dirichlet boundary condition for eigenvalue problems.");
395 for (MooseIndex(values) i = 0; i < values.size(); i++)
398 mooseError(
"Can't set an inhomogeneous array Dirichlet boundary condition for " 399 "eigenvalue problems.");
402 else if (!nbc && !eigen_nbc && !anbc && !aeigen_nbc)
404 "Invalid NodalBC for eigenvalue problems, please use homogeneous (array) Dirichlet.");
409 std::pair<Real, Real>
413 if (n >= n_converged_eigenvalues)
414 mooseError(n,
" not in [0, ", n_converged_eigenvalues,
")");
418 std::pair<Real, Real>
423 if (n >= n_converged_eigenvalues)
424 mooseError(n,
" not in [0, ", n_converged_eigenvalues,
")");
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 " 493 const std::string & )
494 :
libMesh::ParallelObject(eigen_problem)
496 mooseError(
"Need to install SLEPc to solve eigenvalue problems, please reconfigure libMesh\n");
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.
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.
NumericVector< Number > & solution()
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application...
virtual void postAddResidualObject(ResidualObject &object) override
Called after any ResidualObject-derived objects are added to the system.
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.
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
The following methods are specializations for using the libMesh::Parallel::packed_range_* routines fo...
Base class for a system (of equations)
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
EigenProblem & _eigen_problem
unsigned int getNumConvergedEigenvalues() const
Get the number of converged eigenvalues.
Nonlinear system to be solved.
virtual const std::string & name() const
unsigned int getNEigenPairsRequired() const
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.
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.
virtual bool converged() override
Returns the convergence state.
Boundary condition of a Dirichlet type.
virtual void attachPreconditioner(Preconditioner< Number > *preconditioner) override
Attach a customized preconditioner that requires physics knowledge.
bool _precond_matrix_includes_eigen
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.
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)
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.
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.
virtual SNES getSNES() override
Retrieve snes from slepc eigen solver.
std::shared_ptr< TimeIntegrator > _time_integrator
Time integrator.
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
void attachSLEPcCallbacks()
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.
std::vector< std::pair< Real, Real > > _eigen_values
TagID precondMatrixTag() const
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
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.
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.