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             : #include "libmesh/libmesh_config.h"
      13             : 
      14             : #include "NonlinearSystemBase.h"
      15             : #include "SlepcEigenSolverConfiguration.h"
      16             : 
      17             : #include "libmesh/condensed_eigen_system.h"
      18             : 
      19             : // forward declarations
      20             : class EigenProblem;
      21             : class KernelBase;
      22             : class ResidualObject;
      23             : 
      24             : #ifdef LIBMESH_HAVE_SLEPC
      25             : 
      26             : #include <slepceps.h>
      27             : 
      28             : namespace Moose
      29             : {
      30             : void assemble_matrix(EquationSystems & es, const std::string & system_name);
      31             : }
      32             : 
      33             : /**
      34             :  * Nonlinear eigenvalue system to be solved
      35             :  */
      36             : class NonlinearEigenSystem : public NonlinearSystemBase
      37             : {
      38             : public:
      39             :   NonlinearEigenSystem(EigenProblem & problem, const std::string & name);
      40             : 
      41             :   virtual void solve() override;
      42             : 
      43             :   /**
      44             :    * Quit the current solve as soon as possible.
      45             :    */
      46             :   virtual void stopSolve(const ExecFlagType & exec_flag,
      47             :                          const std::set<TagID> & vector_tags_to_close) override;
      48             : 
      49             :   /**
      50             :    * Returns the current nonlinear iteration number.  In libmesh, this is
      51             :    * updated during the nonlinear solve, so it should be up-to-date.
      52             :    */
      53             :   virtual unsigned int getCurrentNonlinearIterationNumber() override;
      54             : 
      55             :   virtual void setupFiniteDifferencedPreconditioner() override;
      56             : 
      57             :   /**
      58             :    * Returns the convergence state
      59             :    *
      60             :    * @return true if converged, otherwise false
      61             :    */
      62             :   virtual bool converged() override;
      63             : 
      64             :   virtual NumericVector<Number> & RHS() override;
      65             : 
      66             :   /*
      67             :    *  A residual vector at MOOSE side for AX
      68             :    */
      69             :   NumericVector<Number> & residualVectorAX();
      70             : 
      71             :   /*
      72             :    * A residual vector at MOOSE side for BX
      73             :    */
      74             :   NumericVector<Number> & residualVectorBX();
      75             : 
      76             :   void attachSLEPcCallbacks();
      77             : 
      78             :   /**
      79             :    * Get the number of converged eigenvalues
      80             :    *
      81             :    * @return The number of converged eigenvalues
      82             :    */
      83        4918 :   unsigned int getNumConvergedEigenvalues() const { return _eigen_sys.get_n_converged(); };
      84             : 
      85             :   virtual unsigned int nNonlinearIterations() const override;
      86             :   virtual unsigned int nLinearIterations() const override;
      87             :   virtual Real finalNonlinearResidual() const override;
      88             : 
      89             :   virtual libMesh::NonlinearSolver<Number> * nonlinearSolver() override;
      90             : 
      91             :   /**
      92             :    * Retrieve snes from slepc eigen solver. It is valid for only nonlinear eigen solver.
      93             :    * You should see a big error if you do this for linear solver.
      94             :    */
      95             :   virtual SNES getSNES() override;
      96             : 
      97             :   /**
      98             :    * Retrieve EPS (SLEPc eigen solver)
      99             :    */
     100             :   virtual EPS getEPS();
     101             : 
     102       42861 :   libMesh::CondensedEigenSystem & sys() { return _eigen_sys; }
     103             : 
     104             :   /**
     105             :    * For eigenvalue problems (including standard and generalized), inhomogeneous (Dirichlet or
     106             :    * Neumann)
     107             :    * boundary conditions are not allowed.
     108             :    */
     109             :   void checkIntegrity();
     110             : 
     111             :   /**
     112             :    * Return the Nth converged eigenvalue.
     113             :    *
     114             :    * @return The Nth converged eigenvalue as a complex number, i.e. the first and the second number
     115             :    * is the real and the imaginary part of
     116             :    * the eigenvalue, respectively.
     117             :    */
     118             :   std::pair<Real, Real> getConvergedEigenvalue(dof_id_type n) const;
     119             : 
     120             :   /**
     121             :    * Return the Nth converged eigenvalue and copies the respective eigen vector to the solution
     122             :    * vector.
     123             :    *
     124             :    * @return The Nth converged eigenvalue as a complex number, i.e. the first and the second number
     125             :    * is the real and the imaginary part of
     126             :    * the eigenvalue, respectively.
     127             :    */
     128             :   std::pair<Real, Real> getConvergedEigenpair(dof_id_type n) const;
     129             : 
     130             :   /**
     131             :    * Get the number of converged eigenvalues
     132             :    *
     133             :    * @return all converged eigenvalues as complex numbers
     134             :    */
     135         668 :   const std::vector<std::pair<Real, Real>> & getAllConvergedEigenvalues() const
     136             :   {
     137         668 :     return _eigen_values;
     138             :   }
     139             : 
     140             :   /**
     141             :    * Vector tag ID of right hand side
     142             :    */
     143       47302 :   TagID eigenVectorTag() const { return _Bx_tag; }
     144             : 
     145             :   /**
     146             :    * Vector tag ID of left hand side
     147             :    */
     148       57358 :   TagID nonEigenVectorTag() const { return _Ax_tag; }
     149             : 
     150             :   /**
     151             :    * Matrix tag ID of right hand side
     152             :    */
     153        4307 :   TagID eigenMatrixTag() const { return _B_tag; }
     154             : 
     155             :   /**
     156             :    * Matrix tag ID of left hand side
     157             :    */
     158        4340 :   TagID nonEigenMatrixTag() const { return _A_tag; }
     159             : 
     160             :   std::set<TagID> defaultVectorTags() const override;
     161             :   std::set<TagID> defaultMatrixTags() const override;
     162             : 
     163             :   /**
     164             :    * If the preconditioning matrix includes eigen kernels
     165             :    */
     166         562 :   void precondMatrixIncludesEigenKernels(bool precond_matrix_includes_eigen)
     167             :   {
     168         562 :     _precond_matrix_includes_eigen = precond_matrix_includes_eigen;
     169         562 :   }
     170             : 
     171             :   bool precondMatrixIncludesEigenKernels() const { return _precond_matrix_includes_eigen; }
     172             : 
     173        3858 :   TagID precondMatrixTag() const { return _precond_tag; }
     174             : 
     175             :   virtual void attachPreconditioner(libMesh::Preconditioner<Number> * preconditioner) override;
     176             : 
     177        1431 :   libMesh::Preconditioner<Number> * preconditioner() const { return _preconditioner; }
     178             : 
     179             :   virtual void turnOffJacobian() override;
     180             : 
     181             :   void residualAndJacobianTogether() override;
     182             : 
     183             :   /**
     184             :    * Initialize the condensed matrices. This is a no-op if there are no
     185             :    * constraints in the DofMap
     186             :    */
     187             :   void initializeCondensedMatrices();
     188             : 
     189             :   virtual void postInit() override;
     190             :   virtual void reinit() override;
     191             : 
     192             : protected:
     193             :   virtual void postAddResidualObject(ResidualObject & object) override;
     194             : 
     195             :   void computeScalingJacobian() override;
     196             :   void computeScalingResidual() override;
     197             : 
     198             :   libMesh::CondensedEigenSystem & _eigen_sys;
     199             :   EigenProblem & _eigen_problem;
     200             :   std::unique_ptr<SlepcEigenSolverConfiguration> _solver_configuration;
     201             :   std::vector<std::pair<Real, Real>> _eigen_values;
     202             :   unsigned int _n_eigen_pairs_required;
     203             :   NumericVector<Number> & _work_rhs_vector_AX;
     204             :   NumericVector<Number> & _work_rhs_vector_BX;
     205             :   TagID _Ax_tag;
     206             :   TagID _Bx_tag;
     207             :   TagID _A_tag;
     208             :   TagID _B_tag;
     209             :   TagID _precond_tag;
     210             :   bool _precond_matrix_includes_eigen;
     211             :   // Libmesh preconditioner
     212             :   libMesh::Preconditioner<Number> * _preconditioner;
     213             : 
     214             :   /// The number of degrees of freedom constrained at the libMesh level, e.g. via hanging node or
     215             :   /// periodic boundary constraints
     216             :   dof_id_type _num_constrained_dofs;
     217             : };
     218             : 
     219             : #else
     220             : 
     221             : class NonlinearEigenSystem : public libMesh::ParallelObject
     222             : {
     223             : public:
     224             :   NonlinearEigenSystem(EigenProblem & problem, const std::string & name);
     225             : 
     226             :   /**
     227             :    * Returns the convergence state
     228             :    * @return true if converged, otherwise false
     229             :    */
     230             :   bool converged() { return false; }
     231             : 
     232             :   void checkIntegrity() {}
     233             : };
     234             : 
     235             : #endif