10 #include "libmesh/libmesh_config.h" 24 #include "libmesh/system.h" 25 #include "libmesh/eigen_solver.h" 26 #include "libmesh/enum_eigen_solver_type.h" 35 params.
addParam<
bool>(
"negative_sign_eigen_kernel",
37 "Whether or not to use a negative sign for eigenvalue kernels. " 38 "Using a negative sign makes eigenvalue kernels consistent with " 39 "a nonlinear solver");
44 "Which eigenvector is used to compute residual and also associated to nonlinear variable");
45 params.
addParam<PostprocessorName>(
"bx_norm",
"A postprocessor describing the norm of Bx");
52 #ifdef LIBMESH_HAVE_SLEPC
55 _n_eigen_pairs_required(1),
56 _generalized_eigenvalue_problem(false),
57 _negative_sign_eigen_kernel(getParam<bool>(
"negative_sign_eigen_kernel")),
58 _active_eigen_index(getParam<unsigned
int>(
"active_eigen_index")),
59 _do_free_power_iteration(false),
60 _output_inverse_eigenvalue(false),
61 _on_linear_solver(false),
62 _matrices_formed(false),
63 _constant_matrices(false),
64 _has_normalization(false),
66 _first_solve(declareRestartableData<bool>(
"first_solve", true)),
67 _bx_norm_name(isParamValid(
"bx_norm")
68 ?
std::make_optional(getParam<PostprocessorName>(
"bx_norm"))
72 #ifdef LIBMESH_HAVE_SLEPC 75 "eigen problems do not currently support multiple nonlinear eigen systems");
81 nl = std::make_shared<NonlinearEigenSystem>(*
this, sys_name);
86 _aux = std::make_shared<AuxiliarySystem>(*
this,
"aux0");
94 mooseError(
"Need to install SLEPc to solve eigenvalue problems, please reconfigure\n");
98 #if PETSC_RELEASE_LESS_THAN(3, 13, 0) 100 "Please use SLEPc-3.13.0 or higher. Old versions of SLEPc likely produce bad convergence");
109 #ifdef LIBMESH_HAVE_SLEPC 113 switch (eigen_problem_type)
141 mooseError(
"libMesh does not support EPT_POS_GEN_NON_HERMITIAN currently \n");
169 TIME_SECTION(
"computeJacobianTag", 3);
173 _nl_eigen->disassociateDefaultMatrixTags();
181 for (
const auto & matrix_tag : matrix_tags)
182 if (
_nl_eigen->hasMatrix(matrix_tag.second))
187 _nl_eigen->associateMatrixToTag(jacobian, tag);
191 _nl_eigen->disassociateMatrixFromTag(jacobian, tag);
198 const std::set<TagID> & tags)
200 TIME_SECTION(
"computeMatricesTags", 3);
202 if (jacobians.size() != tags.size())
205 " does not equal the number of tags ",
210 _nl_eigen->disassociateDefaultMatrixTags();
217 for (
auto tag : tags)
218 _nl_eigen->associateMatrixToTag(*(jacobians[i++]), tag);
223 for (
auto tag : tags)
224 _nl_eigen->disassociateMatrixFromTag(*(jacobians[i++]), tag);
229 const unsigned int nl_sys_num)
231 TIME_SECTION(
"computeJacobianBlocks", 3);
253 TIME_SECTION(
"computeJacobianAB", 3);
257 _nl_eigen->disassociateDefaultMatrixTags();
266 for (
const auto & matrix_tag : matrix_tags)
267 if (
_nl_eigen->hasMatrix(matrix_tag.second))
272 _nl_eigen->associateMatrixToTag(jacobianA, tagA);
273 _nl_eigen->associateMatrixToTag(jacobianB, tagB);
277 _nl_eigen->disassociateMatrixFromTag(jacobianA, tagA);
278 _nl_eigen->disassociateMatrixFromTag(jacobianB, tagB);
286 TIME_SECTION(
"computeResidualTag", 3);
290 _nl_eigen->disassociateDefaultVectorTags();
298 for (
const auto & vector_tag : residual_vector_tags)
299 if (
_nl_eigen->hasVector(vector_tag._id))
302 _nl_eigen->associateVectorToTag(residual, tag);
308 _nl_eigen->disassociateVectorFromTag(residual, tag);
318 TIME_SECTION(
"computeResidualAB", 3);
322 _nl_eigen->disassociateDefaultVectorTags();
331 for (
const auto & vector_tag : residual_vector_tags)
332 if (
_nl_eigen->hasVector(vector_tag._id))
335 _nl_eigen->associateVectorToTag(residualA, tagA);
336 _nl_eigen->associateVectorToTag(residualB, tagB);
342 _nl_eigen->disassociateVectorFromTag(residualA, tagA);
343 _nl_eigen->disassociateVectorFromTag(residualB, tagB);
355 Real eigenvalue = 1.0;
357 if (_active_eigen_index < _nl_eigen->getNumConvergedEigenvalues())
361 _nl_eigen->residualVectorBX() *= eigenvalue;
369 return _nl_eigen->residualVectorAX().l2_norm();
376 for (
auto & vn : var_names)
385 for (
unsigned int vc = 0; vc < var->
count(); ++vc)
387 std::set<dof_id_type> var_indices;
388 _nl_eigen->system().local_dof_indices(var->
number() + vc, var_indices);
389 for (
const auto & dof : var_indices)
424 Real v =
std::sqrt(eig.first * eig.first + eig.second * eig.second);
441 mooseError(
"Number of converged eigenvalues ",
443 " but you required eigenvalue ",
449 v = 1 /
std::sqrt(eig.first * eig.first + eig.second * eig.second);
458 mooseAssert(v != 0.,
"normal factor can not be zero");
460 unsigned int itr = 0;
466 mooseError(
"Can not scale eigenvector to the required factor ",
468 " please check if postprocessor is defined on only eigen variables");
470 mooseAssert(c != 0.,
"postprocessor value used for scaling can not be zero");
509 mooseError(
"There is no executioner for this moose app");
517 #if !PETSC_RELEASE_LESS_THAN(3, 12, 0) 527 TIME_SECTION(
"solve", 1);
535 if (_active_eigen_index < _nl_eigen->getNumConvergedEigenvalues())
549 _console << std::endl <<
" -------------------------------" << std::endl;
550 _console <<
" Free power iteration starts ..." << std::endl;
551 _console <<
" -------------------------------" << std::endl << std::endl;
559 _console << std::endl <<
" --------------------------------------" << std::endl;
560 _console <<
" Extra Free power iteration starts ..." << std::endl;
561 _console <<
" --------------------------------------" << std::endl << std::endl;
569 _console << std::endl <<
" -------------------------------------" << std::endl;
572 _console <<
" Nonlinear Newton iteration starts ..." << std::endl;
574 _console <<
" Nonlinear power iteration starts ..." << std::endl;
576 _console <<
" -------------------------------------" << std::endl << std::endl;
590 #if !PETSC_RELEASE_LESS_THAN(3, 12, 0) 614 #if !PETSC_RELEASE_LESS_THAN(3, 13, 0) 654 "We should not get here unless a bx_norm postprocessor has been provided");
Nonlinear eigenvalue system to be solved.
virtual void update(bool update_libmesh_system=true)
Update the system (doing libMesh magic)
Generalized Non-Hermitian.
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.
static InputParameters validParams()
void computeJacobianBlocks(std::vector< JacobianBlock *> &blocks)
Computes several Jacobian blocks simultaneously, summing their contributions into smaller preconditio...
bool _matrices_formed
Whether or not matrices had been formed.
Generalized Hermitian indefinite.
Newton-based eigensolver with an assembled Jacobian matrix (fully coupled by default) ...
Moose::EigenSolveType _eigen_solve_type
bool isUltimateMaster() const
Whether or not this app is the ultimate master app.
virtual void initPetscOutputAndSomeSolverSettings() override
Hook up monitors for SNES and KSP.
SolverParams & solverParams()
Get the solver parameters.
PostprocessorName _normalization
Postprocessor used to compute a factor from eigenvector.
bool absoluteFuzzyEqual(const T &var1, const T2 &var2, const T3 &tol=libMesh::TOLERANCE *libMesh::TOLERANCE)
Function to check whether two variables are equal within an absolute tolerance.
void mooseDeprecated(Args &&... args) const
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...
unsigned int number() const
Get variable number coming from libMesh.
bool _currently_computing_jacobian
Flag to determine whether the problem is currently computing Jacobian.
virtual void init() override
The same as PJFNK except that matrix-vector multiplication is employed to replace residual evaluation...
void initEigenvector(const Real initial_value)
For nonlinear eigen solver, a good initial value can help convergence.
void clearFreeNonlinearPowerIterations(const InputParameters ¶ms)
bool _has_normalization
Whether or not we normalize eigenvector.
EigenProblem(const InputParameters ¶meters)
unsigned int count() const
Get the number of components Note: For standard and vector variables, the number is one...
virtual void execute(const ExecFlagType &exec_type) override
Convenience function for performing execution of MOOSE systems.
PetscOptions _petsc_option_data_base
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 bool hasScalarVariable(const std::string &var_name) const override
Returns a Boolean indicating whether any system contains a variable with the name provided...
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.
use whatever SLPEC has by default
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
bool _negative_sign_eigen_kernel
Whether or not use negative sign for Bx.
Generalized Non-Hermitian with positive (semi-)definite B.
Specialization of SubProblem for solving nonlinear equations plus auxiliary equations.
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...
MooseApp & getMooseApp() const
Get the MooseApp this class is associated with.
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 computeResidualTags(const std::set< TagID > &tags)
Form multiple residual vectors and each is associated with one tag.
auto max(const L &left, const R &right)
virtual void newAssemblyArray(std::vector< std::shared_ptr< NonlinearSystemBase >> &nl)
std::shared_ptr< NonlinearEigenSystem > _nl_eigen
bool eigen() const
Whether or not this variable operates on an eigen kernel.
bool _generalized_eigenvalue_problem
virtual EquationSystems & es() override
bool isNonlinearEigenvalueSolver() const
Real formNorm()
Form the Bx norm.
bool & _first_solve
A flag to indicate if it is the first time calling the solve.
std::vector< std::shared_ptr< NonlinearSystemBase > > _nl
The nonlinear systems.
void createTagSolutions()
Create extra tagged solution vectors.
virtual void execute(const ExecFlagType &exec_type)
Convenience function for performing execution of MOOSE systems.
virtual std::vector< VariableName > getVariableNames()
Returns a list of all the variables in the problem (both from the NL and Aux systems.
virtual Real computeResidualL2Norm() override
Compute the residual of Ax - Bx.
Preconditioned Jacobian-free Newton Krylov.
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...
const bool & _solve
Whether or not to actually solve the nonlinear system.
NonlinearSystemBase * _current_nl_sys
The current nonlinear system that we are solving.
void setCurrentNonlinearSystem(const unsigned int nl_sys_num)
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
void createTagVectors()
Create extra tagged vectors and matrices.
virtual void init() override
std::shared_ptr< AuxiliarySystem > _aux
The auxiliary system.
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 MooseVariableScalar & getScalarVariable(const THREAD_ID tid, const std::string &var_name) override
Returns the scalar variable reference from whichever system contains it.
void computeMatricesTags(const NumericVector< Number > &soln, const std::vector< std::unique_ptr< SparseMatrix< Number >>> &jacobians, const std::set< TagID > &tags)
Form several matrices simultaneously.
std::vector< VectorTag > getVectorTags(const std::set< TagID > &tag_ids) const
registerMooseObject("MooseApp", EigenProblem)
MooseApp & _app
The MOOSE application this is associated with.
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 ...
const ExecFlagType EXEC_LINEAR
virtual void checkProblemIntegrity()
Method called to perform a series of sanity checks before a simulation is run.
Real _initial_eigenvalue
A value used for initial normalization.
unsigned int _active_eigen_index
Which eigenvalue is used to compute residual.
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.
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...
virtual void initNullSpaceVectors(const InputParameters ¶meters, std::vector< std::shared_ptr< NonlinearSystemBase >> &nl)
Executioner * getExecutioner() const
Retrieve the Executioner for this App.
virtual void solve() override=0
Solve the system (using libMesh magic)
void solveSetup()
Calls the timestepSetup function for each of the output objects.
const ExecFlagType EXEC_PRE_DISPLACE
const std::vector< NonlinearSystemName > _nl_sys_names
The nonlinear system names.
const ExecFlagType EXEC_NONLINEAR
std::set< TagID > _fe_matrix_tags
void scaleEigenvector(const Real scaling_factor)
Scale eigenvector.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Class for containing MooseEnum item information.
void setFreeNonlinearPowerIterations(unsigned int free_power_iterations)
Set SLEPc/PETSc options to trigger free power iteration.
void postScaleEigenVector()
Normalize eigen vector.
bool relativeFuzzyEqual(const T &var1, const T2 &var2, const T3 &tol=libMesh::TOLERANCE *libMesh::TOLERANCE)
Function to check whether two variables are equal within a relative tolerance.
Jacobian-free Newton Krylov.
virtual bool nlConverged(const unsigned int nl_sys_num) override
std::set< TagID > _fe_vector_tags
void mooseError(Args &&... args) const
Emits an error prefixed with object name and type.
std::shared_ptr< DisplacedProblem > _displaced_problem
const InputParameters & parameters() const
Get the parameters of the object.
Real _normal_factor
Postprocessor target value.
static InputParameters validParams()
const ConsoleStream _console
An instance of helper class to write streams to the Console objects.
std::optional< PostprocessorName > _bx_norm_name
The name of the Postprocessor providing the Bx norm.
Problem for solving eigenvalue problems.
bool bxNormProvided() const
Whether a Bx norm postprocessor has been provided.
virtual void checkProblemIntegrity() override
Method called to perform a series of sanity checks before a simulation is run.
EigenProblemType
Type of the eigen problem.
bool doFreePowerIteration() const
Whether or not we are doing free power iteration.
void setEigenproblemType(Moose::EigenProblemType eigen_problem_type)
Set eigen problem type.
Number initial_value(const Point &, const Parameters &, const std::string &, const std::string &)
void ErrorVector unsigned int
virtual void solve(const unsigned int nl_sys_num) override
auto index_range(const T &sizable)
OutputWarehouse & getOutputWarehouse()
Get the OutputWarehouse objects.
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 computeJacobianTags(const std::set< TagID > &tags)
Form multiple matrices, and each is associated with a tag.
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.
const ExecFlagType EXEC_INITIAL