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"
36 params.
addParam<
bool>(
"negative_sign_eigen_kernel",
38 "Whether or not to use a negative sign for eigenvalue kernels. "
39 "Using a negative sign makes eigenvalue kernels consistent with "
40 "a nonlinear solver");
45 "Which eigenvector is used to compute residual and also associated to nonlinear variable");
46 params.
addParam<PostprocessorName>(
"bx_norm",
"A postprocessor describing the norm of Bx");
56#ifdef LIBMESH_HAVE_SLEPC
59 _n_eigen_pairs_required(1),
60 _generalized_eigenvalue_problem(false),
61 _negative_sign_eigen_kernel(getParam<bool>(
"negative_sign_eigen_kernel")),
62 _active_eigen_index(getParam<unsigned
int>(
"active_eigen_index")),
63 _do_free_power_iteration(false),
64 _output_inverse_eigenvalue(false),
65 _on_linear_solver(false),
66 _matrices_formed(false),
67 _constant_matrices(false),
68 _has_normalization(false),
70 _first_solve(declareRestartableData<bool>(
"first_solve", true)),
71 _bx_norm_name(isParamValid(
"bx_norm")
72 ?
std::make_optional(getParam<PostprocessorName>(
"bx_norm"))
76#ifdef LIBMESH_HAVE_SLEPC
79 "eigen problems do not currently support multiple nonlinear eigen systems");
81 paramError(
"linear_sys_names",
"EigenProblem only works with a single nonlinear eigen system");
87 nl = std::make_shared<NonlinearEigenSystem>(*
this, sys_name);
88 _nl_eigen = std::dynamic_pointer_cast<NonlinearEigenSystem>(nl);
94 _aux = std::make_shared<AuxiliarySystem>(*
this,
"aux0");
102 mooseError(
"Need to install SLEPc to solve eigenvalue problems, please reconfigure libMesh\n");
106#if PETSC_RELEASE_LESS_THAN(3, 13, 0)
108 "Please use SLEPc-3.13.0 or higher. Old versions of SLEPc likely produce bad convergence");
117#ifdef LIBMESH_HAVE_SLEPC
121 switch (eigen_problem_type)
149 mooseError(
"libMesh does not support EPT_POS_GEN_NON_HERMITIAN currently \n");
174 SparseMatrix<Number> & jacobian,
177 TIME_SECTION(
"computeJacobianTag", 3);
181 _nl_eigen->disassociateDefaultMatrixTags();
189 for (
const auto & matrix_tag : matrix_tags)
190 if (
_nl_eigen->hasMatrix(matrix_tag.second))
195 _nl_eigen->associateMatrixToTag(jacobian, tag);
200 _nl_eigen->disassociateMatrixFromTag(jacobian, tag);
205 const std::vector<SparseMatrix<Number> *> & jacobians,
206 const std::set<TagID> & tags)
208 TIME_SECTION(
"computeMatricesTags", 3);
210 if (jacobians.size() != tags.size())
213 " does not equal the number of tags ",
218 _nl_eigen->disassociateDefaultMatrixTags();
225 for (
auto tag : tags)
226 _nl_eigen->associateMatrixToTag(*(jacobians[i++]), tag);
232 for (
auto tag : tags)
233 _nl_eigen->disassociateMatrixFromTag(*(jacobians[i++]), tag);
238 const unsigned int nl_sys_num)
240 TIME_SECTION(
"computeJacobianBlocks", 3);
257 SparseMatrix<Number> & jacobianA,
258 SparseMatrix<Number> & jacobianB,
262 TIME_SECTION(
"computeJacobianAB", 3);
266 _nl_eigen->disassociateDefaultMatrixTags();
275 for (
const auto & matrix_tag : matrix_tags)
276 if (
_nl_eigen->hasMatrix(matrix_tag.second))
281 _nl_eigen->associateMatrixToTag(jacobianA, tagA);
282 _nl_eigen->associateMatrixToTag(jacobianB, tagB);
287 _nl_eigen->disassociateMatrixFromTag(jacobianA, tagA);
288 _nl_eigen->disassociateMatrixFromTag(jacobianB, tagB);
293 NumericVector<Number> & residual,
296 TIME_SECTION(
"computeResidualTag", 3);
300 _nl_eigen->disassociateDefaultVectorTags();
303 mooseAssert(
_fe_vector_tags.empty(),
"This should be empty indicating a clean starting state");
308 for (
const auto & vector_tag : residual_vector_tags)
309 if (
_nl_eigen->hasVector(vector_tag._id))
312 _nl_eigen->associateVectorToTag(residual, tag);
320 _nl_eigen->disassociateVectorFromTag(residual, tag);
325 NumericVector<Number> & residualA,
326 NumericVector<Number> & residualB,
330 TIME_SECTION(
"computeResidualAB", 3);
334 _nl_eigen->disassociateDefaultVectorTags();
337 mooseAssert(
_fe_vector_tags.empty(),
"This should be empty indicating a clean starting state");
343 for (
const auto & vector_tag : residual_vector_tags)
344 if (
_nl_eigen->hasVector(vector_tag._id))
347 _nl_eigen->associateVectorToTag(residualA, tagA);
348 _nl_eigen->associateVectorToTag(residualB, tagB);
355 _nl_eigen->disassociateVectorFromTag(residualA, tagA);
356 _nl_eigen->disassociateVectorFromTag(residualB, tagB);
368 Real eigenvalue = 1.0;
370 if (_active_eigen_index < _nl_eigen->getNumConvergedEigenvalues())
374 _nl_eigen->residualVectorBX() *= eigenvalue;
382 return _nl_eigen->residualVectorAX().l2_norm();
389 for (
auto & vn : var_names)
398 for (
unsigned int vc = 0; vc < var->
count(); ++vc)
400 std::set<dof_id_type> var_indices;
401 _nl_eigen->system().local_dof_indices(var->
number() + vc, var_indices);
402 for (
const auto & dof : var_indices)
403 _nl_eigen->solution().set(dof, scaling ? (
_nl_eigen->solution()(dof) * value) : value);
437 Real v = std::sqrt(eig.first * eig.first + eig.second * eig.second);
441 if (!MooseUtils::absoluteFuzzyEqual(factor, 1))
454 mooseError(
"Number of converged eigenvalues ",
456 " but you required eigenvalue ",
462 v = 1 / std::sqrt(eig.first * eig.first + eig.second * eig.second);
471 mooseAssert(v != 0.,
"normal factor can not be zero");
473 unsigned int itr = 0;
475 while (!MooseUtils::relativeFuzzyEqual(v, c))
479 mooseError(
"Can not scale eigenvector to the required factor ",
481 " please check if postprocessor is defined on only eigen variables");
483 mooseAssert(c != 0.,
"postprocessor value used for scaling can not be zero");
505 paramWarning(
"bx_norm",
"This parameter is only used for nonlinear solve types");
508 pp.paramError(
"execute_on",
509 "If providing the Bx norm, this postprocessor must execute on linear e.g. "
510 "during residual evaluations");
532 mooseError(
"There is no executioner for this moose app");
540#if !PETSC_RELEASE_LESS_THAN(3, 12, 0)
550 TIME_SECTION(
"solve", 1);
558 if (_active_eigen_index < _nl_eigen->getNumConvergedEigenvalues())
572 _console << std::endl <<
" -------------------------------" << std::endl;
573 _console <<
" Free power iteration starts ..." << std::endl;
574 _console <<
" -------------------------------" << std::endl << std::endl;
582 _console << std::endl <<
" --------------------------------------" << std::endl;
583 _console <<
" Extra Free power iteration starts ..." << std::endl;
584 _console <<
" --------------------------------------" << std::endl << std::endl;
592 _console << std::endl <<
" -------------------------------------" << std::endl;
595 _console <<
" Nonlinear Newton iteration starts ..." << std::endl;
597 _console <<
" Nonlinear power iteration starts ..." << std::endl;
599 _console <<
" -------------------------------------" << std::endl << std::endl;
613#if !PETSC_RELEASE_LESS_THAN(3, 12, 0)
615 LibmeshPetscCall(PetscOptionsPop());
637#if PETSC_RELEASE_LESS_THAN(3, 13, 0)
644 "We should have errored during construction if we had more than one nonlinear system");
646 "We should have errored during construction if we had any linear systems");
649 _nl_eigen->sys().assemble_before_solve =
false;
675 const auto & solver_params =
solverParams(eigen_sys_num);
693 "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)