10#include "libmesh/libmesh_config.h"
25#include "libmesh/system.h"
26#include "libmesh/eigen_solver.h"
27#include "libmesh/enum_eigen_solver_type.h"
39 params.
addParam<
bool>(
"negative_sign_eigen_kernel",
41 "Whether or not to use a negative sign for eigenvalue kernels. "
42 "Using a negative sign makes eigenvalue kernels consistent with "
43 "a nonlinear solver");
48 "Which eigenvector is used to compute residual and also associated to nonlinear variable");
49 params.
addParam<PostprocessorName>(
"bx_norm",
"A postprocessor describing the norm of Bx");
59#ifdef LIBMESH_HAVE_SLEPC
62 _n_eigen_pairs_required(1),
63 _generalized_eigenvalue_problem(false),
64 _negative_sign_eigen_kernel(getParam<bool>(
"negative_sign_eigen_kernel")),
65 _active_eigen_index(getParam<unsigned
int>(
"active_eigen_index")),
66 _do_free_power_iteration(false),
67 _output_inverse_eigenvalue(false),
68 _on_linear_solver(false),
69 _matrices_formed(false),
70 _constant_matrices(false),
71 _has_normalization(false),
73 _first_solve(declareRestartableData<bool>(
"first_solve", true)),
74 _bx_norm_name(isParamValid(
"bx_norm")
75 ?
std::make_optional(getParam<PostprocessorName>(
"bx_norm"))
79#ifdef LIBMESH_HAVE_SLEPC
82 "eigen problems do not currently support multiple nonlinear eigen systems");
84 paramError(
"linear_sys_names",
"EigenProblem only works with a single nonlinear eigen system");
90 nl = std::make_shared<NonlinearEigenSystem>(*
this, sys_name);
91 _nl_eigen = std::dynamic_pointer_cast<NonlinearEigenSystem>(nl);
97 _aux = std::make_shared<AuxiliarySystem>(*
this,
"aux0");
105 mooseError(
"Need to install SLEPc to solve eigenvalue problems, please reconfigure libMesh\n");
109#if PETSC_RELEASE_LESS_THAN(3, 13, 0)
111 "Please use SLEPc-3.13.0 or higher. Old versions of SLEPc likely produce bad convergence");
120#ifdef LIBMESH_HAVE_SLEPC
124 switch (eigen_problem_type)
152 mooseError(
"libMesh does not support EPT_POS_GEN_NON_HERMITIAN currently \n");
177 SparseMatrix<Number> & jacobian,
180 TIME_SECTION(
"computeJacobianTag", 3);
184 _nl_eigen->disassociateDefaultMatrixTags();
192 for (
const auto & matrix_tag : matrix_tags)
193 if (
_nl_eigen->hasMatrix(matrix_tag.second))
198 _nl_eigen->associateMatrixToTag(jacobian, tag);
203 _nl_eigen->disassociateMatrixFromTag(jacobian, tag);
208 const std::vector<SparseMatrix<Number> *> & jacobians,
209 const std::set<TagID> & tags)
211 TIME_SECTION(
"computeMatricesTags", 3);
213 if (jacobians.size() != tags.size())
216 " does not equal the number of tags ",
221 _nl_eigen->disassociateDefaultMatrixTags();
228 for (
auto tag : tags)
229 _nl_eigen->associateMatrixToTag(*(jacobians[i++]), tag);
235 for (
auto tag : tags)
236 _nl_eigen->disassociateMatrixFromTag(*(jacobians[i++]), tag);
241 const unsigned int nl_sys_num)
243 TIME_SECTION(
"computeJacobianBlocks", 3);
260 SparseMatrix<Number> & jacobianA,
261 SparseMatrix<Number> & jacobianB,
265 TIME_SECTION(
"computeJacobianAB", 3);
269 _nl_eigen->disassociateDefaultMatrixTags();
278 for (
const auto & matrix_tag : matrix_tags)
279 if (
_nl_eigen->hasMatrix(matrix_tag.second))
284 _nl_eigen->associateMatrixToTag(jacobianA, tagA);
285 _nl_eigen->associateMatrixToTag(jacobianB, tagB);
290 _nl_eigen->disassociateMatrixFromTag(jacobianA, tagA);
291 _nl_eigen->disassociateMatrixFromTag(jacobianB, tagB);
296 NumericVector<Number> & residual,
299 TIME_SECTION(
"computeResidualTag", 3);
303 _nl_eigen->disassociateDefaultVectorTags();
306 mooseAssert(
_fe_vector_tags.empty(),
"This should be empty indicating a clean starting state");
311 for (
const auto & vector_tag : residual_vector_tags)
312 if (
_nl_eigen->hasVector(vector_tag._id))
315 _nl_eigen->associateVectorToTag(residual, tag);
323 _nl_eigen->disassociateVectorFromTag(residual, tag);
328 NumericVector<Number> & residualA,
329 NumericVector<Number> & residualB,
333 TIME_SECTION(
"computeResidualAB", 3);
337 _nl_eigen->disassociateDefaultVectorTags();
340 mooseAssert(
_fe_vector_tags.empty(),
"This should be empty indicating a clean starting state");
346 for (
const auto & vector_tag : residual_vector_tags)
347 if (
_nl_eigen->hasVector(vector_tag._id))
350 _nl_eigen->associateVectorToTag(residualA, tagA);
351 _nl_eigen->associateVectorToTag(residualB, tagB);
358 _nl_eigen->disassociateVectorFromTag(residualA, tagA);
359 _nl_eigen->disassociateVectorFromTag(residualB, tagB);
371 Real eigenvalue = 1.0;
373 if (_active_eigen_index < _nl_eigen->getNumConvergedEigenvalues())
377 _nl_eigen->residualVectorBX() *= eigenvalue;
385 return _nl_eigen->residualVectorAX().l2_norm();
392 for (
auto & vn : var_names)
401 for (
unsigned int vc = 0; vc < var->
count(); ++vc)
403 std::set<dof_id_type> var_indices;
404 _nl_eigen->system().local_dof_indices(var->
number() + vc, var_indices);
405 for (
const auto & dof : var_indices)
406 _nl_eigen->solution().set(dof, scaling ? (
_nl_eigen->solution()(dof) * value) : value);
440 Real v = std::sqrt(eig.first * eig.first + eig.second * eig.second);
444 if (!MooseUtils::absoluteFuzzyEqual(factor, 1))
457 mooseError(
"Number of converged eigenvalues ",
459 " but you required eigenvalue ",
465 v = 1 / std::sqrt(eig.first * eig.first + eig.second * eig.second);
474 mooseAssert(v != 0.,
"normal factor can not be zero");
476 unsigned int itr = 0;
478 while (!MooseUtils::relativeFuzzyEqual(v, c))
482 mooseError(
"Can not scale eigenvector to the required factor ",
484 " please check if postprocessor is defined on only eigen variables");
486 mooseAssert(c != 0.,
"postprocessor value used for scaling can not be zero");
508 paramWarning(
"bx_norm",
"This parameter is only used for nonlinear solve types");
511 pp.paramError(
"execute_on",
512 "If providing the Bx norm, this postprocessor must execute on linear e.g. "
513 "during residual evaluations");
535 mooseError(
"There is no executioner for this moose app");
543#if !PETSC_RELEASE_LESS_THAN(3, 12, 0)
553 TIME_SECTION(
"solve", 1);
561 if (_active_eigen_index < _nl_eigen->getNumConvergedEigenvalues())
575 _console << std::endl <<
" -------------------------------" << std::endl;
576 _console <<
" Free power iteration starts ..." << std::endl;
577 _console <<
" -------------------------------" << std::endl << std::endl;
585 _console << std::endl <<
" --------------------------------------" << std::endl;
586 _console <<
" Extra Free power iteration starts ..." << std::endl;
587 _console <<
" --------------------------------------" << std::endl << std::endl;
595 _console << std::endl <<
" -------------------------------------" << std::endl;
598 _console <<
" Nonlinear Newton iteration starts ..." << std::endl;
600 _console <<
" Nonlinear power iteration starts ..." << std::endl;
602 _console <<
" -------------------------------------" << std::endl << std::endl;
616#if !PETSC_RELEASE_LESS_THAN(3, 12, 0)
618 LibmeshPetscCall(PetscOptionsPop());
640#if PETSC_RELEASE_LESS_THAN(3, 13, 0)
647 "We should have errored during construction if we had more than one nonlinear system");
649 "We should have errored during construction if we had any linear systems");
652 _nl_eigen->sys().assemble_before_solve =
false;
678 const auto & solver_params =
solverParams(eigen_sys_num);
696 "We should not get here unless a bx_norm postprocessor has been provided");
registerMooseObject("MooseApp", EigenProblem)
void mooseError(Args &&... args)
Emit an error message with the given stringified, concatenated args and terminate the application.
void mooseDeprecated(Args &&... args)
Emit a deprecated code/feature message with the given stringified, concatenated args.
const ExecFlagType EXEC_INITIAL
const ExecFlagType EXEC_LINEAR
const ExecFlagType EXEC_NONLINEAR
const ExecFlagType EXEC_PRE_DISPLACE
Number initial_value(const Point &, const Parameters &, const std::string &, const std::string &)
void ErrorVector unsigned int
const ConsoleStream _console
An instance of helper class to write streams to the Console objects.
Problem for solving eigenvalue problems.
virtual std::string solverTypeString(unsigned int solver_sys_num=0) override
Return solver type as a human readable string.
void computeMatricesTags(const NumericVector< Number > &soln, const std::vector< SparseMatrix< Number > * > &jacobians, const std::set< TagID > &tags)
Form several matrices simultaneously.
static InputParameters validParams()
bool & _first_solve
A flag to indicate if it is the first time calling the solve.
virtual void computeJacobianBlocks(std::vector< JacobianBlock * > &blocks, const unsigned int nl_sys_num) override
Computes several Jacobian blocks simultaneously, summing their contributions into smaller preconditio...
virtual bool solverSystemConverged(const unsigned int solver_sys_num) override
void setNormalization(const PostprocessorName &pp, const Real value=std::numeric_limits< Real >::max())
Set postprocessor and normalization factor 'Postprocessor' is often used to compute an integral of ph...
EigenProblem(const InputParameters ¶meters)
unsigned int _active_eigen_index
Which eigenvalue is used to compute residual.
void computeResidualAB(const NumericVector< Number > &soln, NumericVector< Number > &residualA, NumericVector< Number > &residualB, TagID tagA, TagID tagB)
Form two vetors, where each is associated with one tag, through one element-loop.
virtual void checkProblemIntegrity() override
Method called to perform a series of sanity checks before a simulation is run.
virtual void execute(const ExecFlagType &exec_type) override
Convenience function for performing execution of MOOSE systems.
std::shared_ptr< NonlinearEigenSystem > _nl_eigen
std::optional< PostprocessorName > _bx_norm_name
The name of the Postprocessor providing the Bx norm.
bool _negative_sign_eigen_kernel
Whether or not use negative sign for Bx.
bool bxNormProvided() const
Whether a Bx norm postprocessor has been provided.
bool doFreePowerIteration() const
Whether or not we are doing free power iteration.
Real _normal_factor
Postprocessor target value.
void initEigenvector(const Real initial_value)
For nonlinear eigen solver, a good initial value can help convergence.
virtual void solve(const unsigned int nl_sys_num) override
Real _initial_eigenvalue
A value used for initial normalization.
void scaleEigenvector(const Real scaling_factor)
Scale eigenvector.
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.
bool _has_normalization
Whether or not we normalize eigenvector.
virtual Real computeResidualL2Norm() override
Compute the residual of Ax - \lambda Bx.
bool isNonlinearEigenvalueSolver(unsigned int eigen_sys_num) const
virtual void init() override
bool _generalized_eigenvalue_problem
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.
std::vector< std::shared_ptr< NonlinearSystemBase > > _nl
The nonlinear systems.
Real formNorm()
Form the Bx norm.
PostprocessorName _normalization
Postprocessor used to compute a factor from eigenvector.
bool _matrices_formed
Whether or not matrices had been formed.
void postScaleEigenVector()
Normalize eigen vector.
void preScaleEigenVector(const std::pair< Real, Real > &eig)
Eigenvector need to be scaled back if it was scaled in an earlier stage Scaling eigen vector does not...
void adjustEigenVector(const Real value, bool scaling)
Adjust eigen vector by either scaling the existing values or setting new values The operations are ap...
void doFreeNonlinearPowerIterations(unsigned int free_power_iterations)
Do some free/extra power iterations.
virtual void initPetscOutputAndSomeSolverSettings() override
Hook up monitors for SNES and KSP.
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.
void setEigenproblemType(Moose::EigenProblemType eigen_problem_type)
Set eigen problem type.
Specialization of SubProblem for solving nonlinear equations plus auxiliary equations.
std::shared_ptr< AuxiliarySystem > _aux
The auxiliary system.
std::set< TagID > _fe_vector_tags
virtual void init() override
virtual void checkProblemIntegrity()
Method called to perform a series of sanity checks before a simulation is run.
void createTagVectors()
Create extra tagged vectors and matrices.
const bool & _solve
Whether or not to actually solve the nonlinear system.
virtual bool hasScalarVariable(const std::string &var_name) const override
Returns a Boolean indicating whether any system contains a variable with the name provided.
virtual void computeResidualTags(const std::set< TagID > &tags)
Form multiple residual vectors and each is associated with one tag.
virtual std::size_t numLinearSystems() const override
virtual MooseVariableScalar & getScalarVariable(const THREAD_ID tid, const std::string &var_name) override
Returns the scalar variable reference from whichever system contains it.
virtual libMesh::EquationSystems & es() override
virtual std::size_t numNonlinearSystems() const override
void computeSystems(const ExecFlagType &type)
Do generic system computations.
const PostprocessorValue & getPostprocessorValueByName(const PostprocessorName &name, std::size_t t_index=0) const
Get a read-only reference to the value associated with a Postprocessor that exists.
const bool _use_hash_table_matrix_assembly
Whether to assemble matrices using hash tables instead of preallocating matrix memory.
virtual void newAssemblyArray(std::vector< std::shared_ptr< SolverSystem > > &solver_systems)
std::set< TagID > _fe_matrix_tags
void setCurrentNonlinearSystem(const unsigned int nl_sys_num)
const std::vector< LinearSystemName > _linear_sys_names
The linear system names.
SolverParams & solverParams(unsigned int solver_sys_num=0)
Get the solver parameters.
virtual const MooseVariableFieldBase & getVariable(const THREAD_ID tid, const std::string &var_name, Moose::VarKindType expected_var_type=Moose::VarKindType::VAR_ANY, Moose::VarFieldType expected_var_field_type=Moose::VarFieldType::VAR_FIELD_ANY) const override
Returns the variable reference for requested variable which must be of the expected_var_type (Nonline...
virtual void initNullSpaceVectors(const InputParameters ¶meters, std::vector< std::shared_ptr< NonlinearSystemBase > > &nl)
PetscOptions _petsc_option_data_base
const UserObject & getUserObjectBase(const std::string &name, const THREAD_ID tid=0) const
Get the user object by its name.
NonlinearSystemBase * _current_nl_sys
The current nonlinear system that we are solving.
std::shared_ptr< DisplacedProblem > _displaced_problem
virtual std::vector< VariableName > getVariableNames()
Returns a list of all the variables in the problem (both from the NL and Aux systems.
virtual void execute(const ExecFlagType &exec_type)
Convenience function for performing execution of MOOSE systems.
virtual void computeJacobianTags(const std::set< TagID > &tags)
Form multiple matrices, and each is associated with a tag.
void createTagSolutions()
Create extra tagged solution vectors.
std::vector< std::shared_ptr< SolverSystem > > _solver_systems
Combined container to base pointer of every solver system.
const std::vector< NonlinearSystemName > _nl_sys_names
The nonlinear system names.
static InputParameters validParams()
OutputWarehouse & getOutputWarehouse()
Get the OutputWarehouse objects.
bool isRestarting() const
Whether or not this is a "restart" calculation.
Executioner * getExecutioner() const
Retrieve the Executioner for this App.
bool isUltimateMaster() const
Whether or not this app is the ultimate master app.
const InputParameters & parameters() const
Get the parameters of the object.
void paramError(const std::string ¶m, Args... args) const
Emits an error prefixed with the file and line number of the given param (from the input file) along ...
MooseApp & getMooseApp() const
Get the MooseApp this class is associated with.
Class for containing MooseEnum item information.
MooseApp & _app
The MOOSE application this is associated with.
bool eigen() const
Whether or not this variable operates on an eigen kernel.
unsigned int number() const
Get variable number coming from libMesh.
unsigned int count() const
Get the number of components Note: For standard and vector variables, the number is one.
void computeJacobianBlocks(std::vector< JacobianBlock * > &blocks)
Computes several Jacobian blocks simultaneously, summing their contributions into smaller preconditio...
virtual void solve() override=0
Solve the system (using libMesh magic)
void solveSetup()
Calls the timestepSetup function for each of the output objects.
void paramWarning(const std::string ¶m, Args... args) const
std::vector< VectorTag > getVectorTags(const std::set< TagID > &tag_ids) const
bool _currently_computing_jacobian
Flag to determine whether the problem is currently computing Jacobian.
virtual std::map< TagName, TagID > & getMatrixTags()
Return all matrix tags in the system, where a tag is represented by a map from name to ID.
void update()
Update the system (doing libMesh magic)
T & set(const std::string &)
void setFreeNonlinearPowerIterations(unsigned int free_power_iterations)
Set SLEPc/PETSc options to trigger free power iteration.
void clearFreeNonlinearPowerIterations(const InputParameters ¶ms)
std::string stringify(const T &t)
conversion to string
@ EST_PJFNKMO
The same as PJFNK except that matrix-vector multiplication is employed to replace residual evaluation...
@ EST_JFNK
Jacobian-free Newton Krylov.
@ EST_NEWTON
Newton-based eigensolver with an assembled Jacobian matrix (fully coupled by default)
@ EST_NONLINEAR_POWER
Nonlinear inverse power.
@ EST_PJFNK
Preconditioned Jacobian-free Newton Krylov.
EigenProblemType
Type of the eigen problem.
@ EPT_GEN_INDEFINITE
Generalized Hermitian indefinite.
@ EPT_NON_HERMITIAN
Non-Hermitian.
@ EPT_GEN_HERMITIAN
Generalized Hermitian.
@ EPT_HERMITIAN
Hermitian.
@ EPT_GEN_NON_HERMITIAN
Generalized Non-Hermitian.
@ EPT_POS_GEN_NON_HERMITIAN
Generalized Non-Hermitian with positive (semi-)definite B.
@ EPT_SLEPC_DEFAULT
use whatever SLPEC has by default
void assemble_matrix(EquationSystems &es, const std::string &system_name)