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