1 // The libMesh Finite Element Library. 2 // Copyright (C) 2002-2020 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner 3 4 // This library is free software; you can redistribute it and/or 5 // modify it under the terms of the GNU Lesser General Public 6 // License as published by the Free Software Foundation; either 7 // version 2.1 of the License, or (at your option) any later version. 8 9 // This library is distributed in the hope that it will be useful, 10 // but WITHOUT ANY WARRANTY; without even the implied warranty of 11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12 // Lesser General Public License for more details. 13 14 // You should have received a copy of the GNU Lesser General Public 15 // License along with this library; if not, write to the Free Software 16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 17 18 19 20 #ifndef LIBMESH_EIGEN_TIME_SOLVER_H 21 #define LIBMESH_EIGEN_TIME_SOLVER_H 22 23 #include "libmesh/libmesh_config.h" 24 #ifdef LIBMESH_HAVE_SLEPC 25 26 // Local includes 27 #include "libmesh/time_solver.h" 28 29 // C++ includes 30 31 namespace libMesh 32 { 33 34 // Forward declarations 35 template <typename T> class EigenSolver; 36 37 38 /** 39 * The name of this class is confusing...it's meant to refer to the 40 * base class (TimeSolver) while still telling one that it's for solving 41 * (generalized) EigenValue problems that arise from finite element 42 * discretizations. For a time-dependent problem du/dt=F(u), with a 43 * steady solution 0=F(u_0), we look at the time evolution of a small 44 * perturbation, p=u-u_0, for which the (linearized) governing equation is 45 * 46 * dp/dt = F'(u_0)p 47 * 48 * where F'(u_0) is the Jacobian. The generalized eigenvalue problem arises 49 * by considering perturbations of the general form p = exp(lambda*t)x, which 50 * leads to 51 * 52 * Ax = lambda*Bx 53 * 54 * where A is the (discretized by FEM) Jacobian matrix and B is the 55 * (discretized by FEM) mass matrix. 56 * 57 * The EigenSystem class (by Steffen Petersen) is related but does not 58 * fall under the FEMSystem paradigm invented by Roy Stogner. The EigenSolver 59 * class (also by Steffen) is meant to provide a generic "linear solver" 60 * interface for EigenValue software. The only current concrete implementation 61 * is a SLEPc-based eigensolver class, which we make use of here as well. 62 * 63 * \author John W. Peterson 64 * \date 2007 65 */ 66 class EigenTimeSolver : public TimeSolver 67 { 68 public: 69 /** 70 * The type of system 71 */ 72 typedef DifferentiableSystem sys_type; 73 74 /** 75 * The parent class 76 */ 77 typedef TimeSolver Parent; 78 79 /** 80 * Constructor. Requires a reference to the system 81 * to be solved. 82 */ 83 explicit 84 EigenTimeSolver (sys_type & s); 85 86 /** 87 * Destructor. 88 */ 89 virtual ~EigenTimeSolver (); 90 91 /** 92 * The initialization function. This method is used to 93 * initialize internal data structures before a simulation begins. 94 */ 95 virtual void init () override; 96 97 /** 98 * The reinitialization function. This method is used after 99 * changes in the mesh 100 */ 101 virtual void reinit () override; 102 103 /** 104 * Implements the assembly of both matrices A and B, and calls 105 * the EigenSolver to compute the eigenvalues. 106 */ 107 virtual void solve () override; 108 109 /** 110 * It doesn't make sense to advance the timestep, so we shouldn't call this. 111 */ advance_timestep()112 virtual void advance_timestep () override {} 113 114 /** 115 * error convergence order against deltat is 116 * not applicable to an eigenvalue problem. 117 */ error_order()118 Real error_order() const { return 0.; } 119 120 /** 121 * Forms either the spatial (Jacobian) or mass matrix part of the 122 * operator, depending on which is requested. 123 */ 124 virtual bool element_residual (bool get_jacobian, 125 DiffContext &) override; 126 127 /** 128 * Forms the jacobian of the boundary terms. 129 */ 130 virtual bool side_residual (bool get_jacobian, 131 DiffContext &) override; 132 133 /** 134 * Forms the jacobian of the nonlocal terms. 135 */ 136 virtual bool nonlocal_residual (bool get_jacobian, 137 DiffContext &) override; 138 139 /** 140 * \returns 0, but derived classes should override this function to 141 * compute the size of the difference between successive solution 142 * iterates ||u^{n+1} - u^{n}|| in some norm. 143 */ du(const SystemNorm &)144 virtual Real du (const SystemNorm &) const override { return 0.; } 145 146 /** 147 * This is effectively a steady-state solver. 148 */ is_steady()149 virtual bool is_steady() const override { return true; } 150 151 /** 152 * The EigenSolver object. This is what actually 153 * makes the calls to SLEPc. 154 */ 155 std::unique_ptr<EigenSolver<Number>> eigen_solver; 156 157 /** 158 * The linear solver tolerance to be used when solving the 159 * eigenvalue problem. FIXME: need more info... 160 */ 161 Real tol; 162 163 /** 164 * The maximum number of iterations allowed to solve the problem. 165 */ 166 unsigned int maxits; 167 168 /** 169 * The number of eigenvectors/values to be computed. 170 */ 171 unsigned int n_eigenpairs_to_compute; 172 173 /** 174 * The number of basis vectors to use in the computation. According 175 * to ex16, the number of basis vectors must be >= the number of 176 * eigenpairs requested, and ncv >= 2*nev is recommended. 177 * Increasing this number, even by a little bit, can _greatly_ 178 * reduce the number of (EigenSolver) iterations required to compute 179 * the desired number of eigenpairs, but the _cost per iteration_ 180 * goes up drastically as well. 181 */ 182 unsigned int n_basis_vectors_to_use; 183 184 /** 185 * After a solve, holds the number of eigenpairs successfully 186 * converged. 187 */ 188 unsigned int n_converged_eigenpairs; 189 190 /** 191 * After a solve, holds the number of iterations required to converge 192 * the requested number of eigenpairs. 193 */ 194 unsigned int n_iterations_reqd; 195 196 private: 197 198 enum NowAssembling { 199 /** 200 * The matrix associated with the spatial part of the operator. 201 */ 202 Matrix_A, 203 204 /** 205 * The matrix associated with the time derivative (mass matrix). 206 */ 207 Matrix_B, 208 209 /** 210 * The enum is in an invalid state. 211 */ 212 Invalid_Matrix 213 }; 214 215 /** 216 * Flag which controls the internals of element_residual() and side_residual(). 217 */ 218 NowAssembling now_assembling; 219 }; 220 221 } // namespace libMesh 222 223 224 #endif // LIBMESH_HAVE_SLEPC 225 #endif // LIBMESH_EIGEN_TIME_SOLVER_H 226