Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       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             : #pragma once
      11             : 
      12             : // MOOSE Includes
      13             : #include "FEProblemBase.h"
      14             : #include "Eigenvalue.h"
      15             : 
      16             : class NonlinearEigenSystem;
      17             : 
      18             : /**
      19             :  * Problem for solving eigenvalue problems
      20             :  */
      21             : class EigenProblem : public FEProblemBase
      22             : {
      23             : public:
      24             :   static InputParameters validParams();
      25             : 
      26             :   EigenProblem(const InputParameters & parameters);
      27             : 
      28             :   virtual std::string solverTypeString(unsigned int solver_sys_num = 0) override;
      29             : 
      30             : #ifdef LIBMESH_HAVE_SLEPC
      31             :   virtual void solve(const unsigned int nl_sys_num) override;
      32             : 
      33             :   virtual void init() override;
      34             : 
      35             :   virtual bool solverSystemConverged(const unsigned int solver_sys_num) override;
      36             : 
      37        1894 :   unsigned int getNEigenPairsRequired() const { return _n_eigen_pairs_required; }
      38         588 :   void setNEigenPairsRequired(unsigned int n_eigen_pairs)
      39             :   {
      40         588 :     _n_eigen_pairs_required = n_eigen_pairs;
      41         588 :   }
      42             :   bool isGeneralizedEigenvalueProblem() const { return _generalized_eigenvalue_problem; }
      43             :   bool isNonlinearEigenvalueSolver(unsigned int eigen_sys_num) const;
      44             :   // silences warning in debug mode about the other computeJacobian signature being hidden
      45             :   using FEProblemBase::computeJacobian;
      46             : 
      47             :   NonlinearEigenSystem & getNonlinearEigenSystem(const unsigned int nl_sys_num);
      48             :   NonlinearEigenSystem & getCurrentNonlinearEigenSystem();
      49             : 
      50             :   virtual void checkProblemIntegrity() override;
      51             : 
      52             :   /**
      53             :    * A flag indicates if a negative sign is used in eigen kernels.
      54             :    * If the negative sign is used, eigen kernels are consistent in nonlinear solver.
      55             :    * In nonlinear solver, RHS kernels always have a negative sign.
      56             :    */
      57       19230 :   bool negativeSignEigenKernel() const { return _negative_sign_eigen_kernel; }
      58             : 
      59             :   /**
      60             :    * Set postprocessor and normalization factor
      61             :    * 'Postprocessor' is often used to compute an integral of physics variables
      62             :    */
      63             :   void setNormalization(const PostprocessorName & pp,
      64             :                         const Real value = std::numeric_limits<Real>::max());
      65             : 
      66             :   /**
      67             :    * Set an initial eigenvalue for initial normalization
      68             :    */
      69         584 :   void setInitialEigenvalue(const Real initial_eigenvalue)
      70             :   {
      71         584 :     _initial_eigenvalue = initial_eigenvalue;
      72         584 :   }
      73             : 
      74             :   /**
      75             :    * Whether or not we are doing free power iteration. It is used in convergence check.
      76             :    * We need to mark the solver as "converged" when doing free power to retrieve
      77             :    * the final solution from SLPEc
      78             :    */
      79        2988 :   bool doFreePowerIteration() const { return _do_free_power_iteration; }
      80             : 
      81             :   /**
      82             :    * Set a flag to indicate whether or not we are doing free power iterations.
      83             :    */
      84        1116 :   void doFreePowerIteration(bool do_power) { _do_free_power_iteration = do_power; }
      85             : 
      86             :   /**
      87             :    * Eigenvector need to be scaled back if it was scaled in an earlier stage
      88             :    * Scaling eigen vector does not affect solution (eigenvalue, eigenvector),
      89             :    * but it does affect the convergence rate. To have an optimal convergence rate,
      90             :    * We pre-scale eigen vector using the same factor as the one computed in
      91             :    * "postScaleEigenVector"
      92             :    */
      93             :   void preScaleEigenVector(const std::pair<Real, Real> & eig);
      94             : 
      95             :   /**
      96             :    * Normalize eigen vector. Scale eigen vector such as ||x|| = _normal_factor
      97             :    * This might be useful when coupling to other physics
      98             :    */
      99             :   void postScaleEigenVector();
     100             : 
     101             :   /**
     102             :    * Scale eigenvector. Scaling_factor is often computed based on physics.
     103             :    */
     104             :   void scaleEigenvector(const Real scaling_factor);
     105             : 
     106             :   /**
     107             :    * Set eigen problem type. Don't need to use this if we use Newton eigenvalue solver.
     108             :    */
     109             :   void setEigenproblemType(Moose::EigenProblemType eigen_problem_type);
     110             : 
     111             :   /**
     112             :    * Compute the residual of Ax - \lambda Bx. If there is no \lambda yet, "1" will
     113             :    * be used.
     114             :    */
     115             :   virtual Real computeResidualL2Norm() override;
     116             : 
     117             :   /**
     118             :    * Form a Jacobian matrix for all kernels and BCs with a given tag
     119             :    */
     120             :   virtual void computeJacobianTag(const NumericVector<Number> & soln,
     121             :                                   SparseMatrix<Number> & jacobian,
     122             :                                   TagID tag) override;
     123             : 
     124             :   /**
     125             :    * Form several matrices simultaneously
     126             :    */
     127             :   void computeMatricesTags(const NumericVector<Number> & soln,
     128             :                            const std::vector<SparseMatrix<Number> *> & jacobians,
     129             :                            const std::set<TagID> & tags);
     130             : 
     131             :   /**
     132             :    * Form two Jacobian matrices, where each is associated with one tag, through one
     133             :    * element-loop.
     134             :    */
     135             :   void computeJacobianAB(const NumericVector<Number> & soln,
     136             :                          SparseMatrix<Number> & jacobianA,
     137             :                          SparseMatrix<Number> & jacobianB,
     138             :                          TagID tagA,
     139             :                          TagID tagB);
     140             : 
     141             :   virtual void computeJacobianBlocks(std::vector<JacobianBlock *> & blocks,
     142             :                                      const unsigned int nl_sys_num) override;
     143             : 
     144             :   /**
     145             :    * Form a vector for all kernels and BCs with a given tag
     146             :    */
     147             :   virtual void computeResidualTag(const NumericVector<Number> & soln,
     148             :                                   NumericVector<Number> & residual,
     149             :                                   TagID tag) override;
     150             : 
     151             :   /**
     152             :    * Form two vetors, where each is associated with one tag, through one
     153             :    * element-loop.
     154             :    */
     155             :   void computeResidualAB(const NumericVector<Number> & soln,
     156             :                          NumericVector<Number> & residualA,
     157             :                          NumericVector<Number> & residualB,
     158             :                          TagID tagA,
     159             :                          TagID tagB);
     160             : 
     161             :   /**
     162             :    * Convenience function for performing execution of MOOSE systems.
     163             :    * We override this function to perform an initial scaling.
     164             :    */
     165             :   virtual void execute(const ExecFlagType & exec_type) override;
     166             : 
     167             :   /**
     168             :    * For nonlinear eigen solver, a good initial value can help convergence.
     169             :    * Should set initial values for only eigen variables.
     170             :    */
     171             :   void initEigenvector(const Real initial_value);
     172             : 
     173             :   /**
     174             :    * Which eigenvalue is active
     175             :    */
     176        1296 :   unsigned int activeEigenvalueIndex() const { return _active_eigen_index; }
     177             : 
     178             :   /**
     179             :    * Hook up monitors for SNES and KSP
     180             :    */
     181             :   virtual void initPetscOutputAndSomeSolverSettings() override;
     182             : 
     183             :   /**
     184             :    * Whether or not to output eigenvalue inverse. The inverse is useful for
     185             :    * neutronics community
     186             :    */
     187        5112 :   bool outputInverseEigenvalue() const { return _output_inverse_eigenvalue; }
     188             : 
     189             :   /**
     190             :    * Set a flag to indicate whether or not to output eigenvalue inverse.
     191             :    */
     192          10 :   void outputInverseEigenvalue(bool inverse) { _output_inverse_eigenvalue = inverse; }
     193             : 
     194             :   /**
     195             :    * Whether or not we are in a linear solver iteration
     196             :    */
     197        9128 :   bool onLinearSolver() const { return _on_linear_solver; }
     198             : 
     199             :   /**
     200             :    * Set a flag to indicate whether or not we are in a linear solver iteration
     201             :    */
     202        9272 :   void onLinearSolver(bool ols) { _on_linear_solver = ols; }
     203             : 
     204             :   /**
     205             :    * Whether or not matrices are constant
     206             :    */
     207        8448 :   bool constantMatrices() const { return _constant_matrices; }
     208             : 
     209             :   /**
     210             :    * Set a flag to indicate whether or not we use constant matrices
     211             :    */
     212         588 :   void constantMatrices(bool cm) { _constant_matrices = cm; }
     213             : 
     214             :   /**
     215             :    * Whether or not constant matrices were already formed
     216             :    */
     217        3412 :   bool wereMatricesFormed() const { return _matrices_formed; }
     218             : 
     219             :   /**
     220             :    * Set a flag to indicate whether or not constant matrices were already formed
     221             :    */
     222         380 :   void wereMatricesFormed(bool mf) { _matrices_formed = mf; }
     223             : 
     224             :   /**
     225             :    * Form the Bx norm
     226             :    */
     227             :   Real formNorm();
     228             : 
     229             :   /**
     230             :    * Whether a Bx norm postprocessor has been provided
     231             :    */
     232        8056 :   bool bxNormProvided() const { return _bx_norm_name.has_value(); }
     233             : 
     234             :   /**
     235             :    * Set the Bx norm postprocessor programatically
     236             :    */
     237             :   void setBxNorm(const PostprocessorName & bx_norm) { _bx_norm_name = bx_norm; }
     238             : 
     239             : protected:
     240             :   unsigned int _n_eigen_pairs_required;
     241             :   bool _generalized_eigenvalue_problem;
     242             :   std::shared_ptr<NonlinearEigenSystem> _nl_eigen;
     243             : 
     244             :   /// Whether or not use negative sign for Bx. Use negative sign by default to
     245             :   /// make the eigen system consistent with nonlinear system
     246             :   bool _negative_sign_eigen_kernel;
     247             :   /// Which eigenvalue is used to compute residual. By default the zeroth eigenvalue
     248             :   /// is used.
     249             :   unsigned int _active_eigen_index;
     250             : 
     251             :   /// Whether or not we are doing free power iteration. Free power iteration is
     252             :   /// often used to compute initial guess for Newton eigen solver. It is automatically
     253             :   /// triggered by Eigenvalue Executioner
     254             :   bool _do_free_power_iteration;
     255             :   /// Whether or not output eigenvalue as its inverse. By default, we output regular eigenvalue.
     256             :   bool _output_inverse_eigenvalue;
     257             :   /// Whether or not we are in linear solver
     258             :   bool _on_linear_solver;
     259             :   /// Whether or not matrices had been formed
     260             :   bool _matrices_formed;
     261             :   /// Whether or not require constant matrices
     262             :   bool _constant_matrices;
     263             :   /// Whether or not we normalize eigenvector
     264             :   bool _has_normalization;
     265             :   /// Postprocessor used to compute a factor from eigenvector
     266             :   PostprocessorName _normalization;
     267             :   /// Postprocessor target value. The value of postprocessor should equal to
     268             :   /// '_normal_factor' by adjusting eigenvector
     269             :   Real _normal_factor;
     270             :   /// A flag to indicate if it is the first time calling the solve
     271             :   bool & _first_solve;
     272             :   /// A value used for initial normalization
     273             :   Real _initial_eigenvalue;
     274             : 
     275             : private:
     276             :   /**
     277             :    * Do some free/extra power iterations
     278             :    */
     279             :   void doFreeNonlinearPowerIterations(unsigned int free_power_iterations);
     280             : 
     281             :   /**
     282             :    * Adjust eigen vector by either scaling the existing values or setting new values
     283             :    * The operations are applied for only eigen variables
     284             :    */
     285             :   void adjustEigenVector(const Real value, bool scaling);
     286             : 
     287             :   /// The name of the Postprocessor providing the Bx norm. This may be empty in which case the
     288             :   /// default L2 norm of Bx will be used as the Bx norm
     289             :   std::optional<PostprocessorName> _bx_norm_name;
     290             : 
     291             : #endif
     292             : 
     293             :   using FEProblemBase::_nl;
     294             : };
     295             : 
     296             : #ifdef LIBMESH_HAVE_SLEPC
     297             : 
     298             : inline NonlinearEigenSystem &
     299         644 : EigenProblem::getNonlinearEigenSystem(const unsigned int nl_sys_num)
     300             : {
     301         644 :   if (nl_sys_num > 0)
     302           0 :     mooseError("eigen problems do not currently support multiple nonlinear eigen systems");
     303         644 :   return *_nl_eigen;
     304             : }
     305             : 
     306             : inline NonlinearEigenSystem &
     307       69374 : EigenProblem::getCurrentNonlinearEigenSystem()
     308             : {
     309       69374 :   return *_nl_eigen;
     310             : }
     311             : 
     312             : #endif