libMesh
eigen_time_solver.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 #include "libmesh/libmesh_config.h"
21 #ifdef LIBMESH_HAVE_SLEPC
22 
23 #include "libmesh/diff_system.h"
24 #include "libmesh/eigen_time_solver.h"
25 #include "libmesh/eigen_solver.h"
26 #include "libmesh/sparse_matrix.h"
27 
28 namespace libMesh
29 {
30 
31 // Constructor
33  : Parent(s),
34  eigen_solver (EigenSolver<Number>::build(s.comm())),
35  tol(pow(TOLERANCE, 5./3.)),
36  maxits(1000),
37  n_eigenpairs_to_compute(5),
38  n_basis_vectors_to_use(3*n_eigenpairs_to_compute),
39  n_converged_eigenpairs(0),
40  n_iterations_reqd(0)
41 {
42  libmesh_experimental();
43  eigen_solver->set_eigenproblem_type(GHEP);//or GNHEP
44  eigen_solver->set_position_of_spectrum(LARGEST_MAGNITUDE);
45 }
46 
48 {
49 }
50 
52 {
53  // empty...
54 }
55 
57 {
58  // Add matrix "B" to _system if not already there.
59  // The user may have already added a matrix "B" before
60  // calling the System initialization. This would be
61  // necessary if e.g. the System originally started life
62  // with a different type of TimeSolver and only later
63  // had its TimeSolver changed to an EigenTimeSolver.
64  if (!_system.have_matrix("B"))
65  _system.add_matrix("B");
66 }
67 
69 {
70  // The standard implementation is basically to call:
71  // _diff_solver->solve();
72  // which internally assembles (when necessary) matrices and vectors
73  // and calls linear solver software while also doing Newton steps (see newton_solver.C)
74  //
75  // The element_residual and side_residual functions below control
76  // what happens in the interior of the element assembly loops.
77  // We have a system reference, so it's possible to call _system.assembly()
78  // ourselves if we want to...
79  //
80  // Interestingly, for the EigenSolver we don't need residuals...just Jacobians.
81  // The Jacobian should therefore always be requested, and always return
82  // jacobian_computed as being true.
83 
84  // The basic plan of attack is:
85  // .) Construct the Jacobian using _system.assembly(true,true) as we
86  // would for a steady system. Use a flag in this class to
87  // control behavior in element_residual and side_residual
88  // .) Swap _system.matrix to matrix "B" (be sure to add this extra matrix during init)
89  // .) Call _system.assembly(true,true) again, using the flag in element_residual
90  // and side_residual to only get the mass matrix terms.
91  // .) Send A and B to Steffen's EigenSolver interface.
92 
93  // Assemble the spatial part (matrix A) of the operator
94  if (!this->quiet)
95  libMesh::out << "Assembling matrix A." << std::endl;
96  _system.matrix = &( _system.get_matrix ("System Matrix") );
97  this->now_assembling = Matrix_A;
98  _system.assembly(true, true);
99  //_system.matrix->print_matlab("matrix_A.m");
100 
101  // Point the system's matrix at B, call assembly again.
102  if (!this->quiet)
103  libMesh::out << "Assembling matrix B." << std::endl;
104  _system.matrix = &( _system.get_matrix ("B") );
105  this->now_assembling = Matrix_B;
106  _system.assembly(true, true);
107  //_system.matrix->print_matlab("matrix_B.m");
108 
109  // Send matrices A, B to Steffen's SlepcEigenSolver interface
110  //libmesh_here();
111  if (!this->quiet)
112  libMesh::out << "Calling the EigenSolver." << std::endl;
113  std::pair<unsigned int, unsigned int> solve_data =
114  eigen_solver->solve_generalized (_system.get_matrix ("System Matrix"),
115  _system.get_matrix ("B"),
118  tol,
119  maxits);
120 
121  this->n_converged_eigenpairs = solve_data.first;
122  this->n_iterations_reqd = solve_data.second;
123 }
124 
125 
126 
127 bool EigenTimeSolver::element_residual(bool request_jacobian,
128  DiffContext & context)
129 {
130  // The EigenTimeSolver always computes jacobians!
131  libmesh_assert (request_jacobian);
132 
133  context.elem_solution_rate_derivative = 1;
134  context.elem_solution_derivative = 1;
135 
136  // Assemble the operator for the spatial part.
137  if (now_assembling == Matrix_A)
138  {
139  bool jacobian_computed =
140  _system.get_physics()->element_time_derivative(request_jacobian, context);
141 
142  // The user shouldn't compute a jacobian unless requested
143  libmesh_assert(request_jacobian || !jacobian_computed);
144 
145  bool jacobian_computed2 =
146  _system.get_physics()->element_constraint(jacobian_computed, context);
147 
148  // The user shouldn't compute a jacobian unless requested
149  libmesh_assert (jacobian_computed || !jacobian_computed2);
150 
151  return jacobian_computed && jacobian_computed2;
152 
153  }
154 
155  // Assemble the mass matrix operator
156  else if (now_assembling == Matrix_B)
157  {
158  bool mass_jacobian_computed =
159  _system.get_physics()->mass_residual(request_jacobian, context);
160 
161  // Scale Jacobian by -1 to get positive matrix from new negative
162  // mass_residual convention
163  context.get_elem_jacobian() *= -1.0;
164 
165  return mass_jacobian_computed;
166  }
167 
168  else
169  libmesh_error_msg("Unrecognized value now_assembling = " << now_assembling);
170 }
171 
172 
173 
174 bool EigenTimeSolver::side_residual(bool request_jacobian,
175  DiffContext & context)
176 {
177  // The EigenTimeSolver always requests jacobians?
178  //libmesh_assert (request_jacobian);
179 
180  context.elem_solution_rate_derivative = 1;
181  context.elem_solution_derivative = 1;
182 
183  // Assemble the operator for the spatial part.
184  if (now_assembling == Matrix_A)
185  {
186  bool jacobian_computed =
187  _system.get_physics()->side_time_derivative(request_jacobian, context);
188 
189  // The user shouldn't compute a jacobian unless requested
190  libmesh_assert (request_jacobian || !jacobian_computed);
191 
192  bool jacobian_computed2 =
193  _system.get_physics()->side_constraint(jacobian_computed, context);
194 
195  // The user shouldn't compute a jacobian unless requested
196  libmesh_assert (jacobian_computed || !jacobian_computed2);
197 
198  return jacobian_computed && jacobian_computed2;
199 
200  }
201 
202  // There is now a "side" equivalent for the mass matrix
203  else if (now_assembling == Matrix_B)
204  {
205  bool mass_jacobian_computed =
206  _system.get_physics()->side_mass_residual(request_jacobian, context);
207 
208  // Scale Jacobian by -1 to get positive matrix from new negative
209  // mass_residual convention
210  context.get_elem_jacobian() *= -1.0;
211 
212  return mass_jacobian_computed;
213  }
214 
215  else
216  libmesh_error_msg("Unrecognized value now_assembling = " << now_assembling);
217 }
218 
219 
220 
221 bool EigenTimeSolver::nonlocal_residual(bool request_jacobian,
222  DiffContext & context)
223 {
224  // The EigenTimeSolver always requests jacobians?
225  //libmesh_assert (request_jacobian);
226 
227  // Assemble the operator for the spatial part.
228  if (now_assembling == Matrix_A)
229  {
230  bool jacobian_computed =
231  _system.get_physics()->nonlocal_time_derivative(request_jacobian, context);
232 
233  // The user shouldn't compute a jacobian unless requested
234  libmesh_assert (request_jacobian || !jacobian_computed);
235 
236  bool jacobian_computed2 =
237  _system.get_physics()->nonlocal_constraint(jacobian_computed, context);
238 
239  // The user shouldn't compute a jacobian unless requested
240  libmesh_assert (jacobian_computed || !jacobian_computed2);
241 
242  return jacobian_computed && jacobian_computed2;
243 
244  }
245 
246  // There is now a "side" equivalent for the mass matrix
247  else if (now_assembling == Matrix_B)
248  {
249  bool mass_jacobian_computed =
250  _system.get_physics()->nonlocal_mass_residual(request_jacobian, context);
251 
252  // Scale Jacobian by -1 to get positive matrix from new negative
253  // mass_residual convention
254  context.get_elem_jacobian() *= -1.0;
255 
256  return mass_jacobian_computed;
257  }
258 
259  else
260  libmesh_error_msg("Unrecognized value now_assembling = " << now_assembling);
261 }
262 
263 } // namespace libMesh
264 
265 #endif // LIBMESH_HAVE_SLEPC
bool quiet
Print extra debugging information if quiet == false.
Definition: time_solver.h:191
virtual bool side_constraint(bool request_jacobian, DiffContext &)
Adds the constraint contribution on side of elem to elem_residual.
Definition: diff_physics.h:190
This class provides all data required for a physics package (e.g.
Definition: diff_context.h:54
virtual bool element_time_derivative(bool request_jacobian, DiffContext &)
Adds the time derivative contribution on elem to elem_residual.
Definition: diff_physics.h:123
This is a generic class that defines a solver to handle time integration of DifferentiableSystems.
Definition: time_solver.h:58
virtual void solve() libmesh_override
Implements the assembly of both matrices A and B, and calls the EigenSolver to compute the eigenvalue...
virtual void reinit() libmesh_override
The reinitialization function.
SparseMatrix< Number > & add_matrix(const std::string &mat_name)
Adds the additional matrix mat_name to this system.
virtual bool nonlocal_time_derivative(bool request_jacobian, DiffContext &)
Adds any nonlocal time derivative contributions (e.g.
Definition: diff_physics.h:208
static const Real TOLERANCE
The libMesh namespace provides an interface to certain functionality in the library.
virtual bool nonlocal_constraint(bool request_jacobian, DiffContext &)
Adds any nonlocal constraint contributions (e.g.
Definition: diff_physics.h:226
virtual bool mass_residual(bool request_jacobian, DiffContext &)
Subtracts a mass vector contribution on elem from elem_residual.
Definition: diff_physics.h:317
libmesh_assert(j)
virtual void assembly(bool get_residual, bool get_jacobian, bool apply_heterogeneous_constraints=false, bool apply_no_constraints=false) libmesh_override=0
Assembles a residual in rhs and/or a jacobian in matrix, as requested.
This class provides a specific system class.
Definition: diff_system.h:53
const DenseMatrix< Number > & get_elem_jacobian() const
Const accessor for element Jacobian.
Definition: diff_context.h:282
sys_type & _system
A reference to the system we are solving.
Definition: time_solver.h:257
virtual bool nonlocal_residual(bool get_jacobian, DiffContext &) libmesh_override
Forms the jacobian of the nonlocal terms.
Real elem_solution_rate_derivative
The derivative of elem_solution_rate with respect to the current nonlinear solution, for use by systems with non default mass_residual terms.
Definition: diff_context.h:490
bool have_matrix(const std::string &mat_name) const
Real elem_solution_derivative
The derivative of elem_solution with respect to the current nonlinear solution.
Definition: diff_context.h:483
UniquePtr< EigenSolver< Number > > eigen_solver
The EigenSolver object.
const DifferentiablePhysics * get_physics() const
Definition: diff_system.h:169
virtual void init() libmesh_override
The initialization function.
virtual bool element_residual(bool get_jacobian, DiffContext &) libmesh_override
Forms either the spatial (Jacobian) or mass matrix part of the operator, depending on which is reques...
double pow(double a, int b)
unsigned int maxits
The maximum number of iterations allowed to solve the problem.
unsigned int n_iterations_reqd
After a solve, holds the number of iterations required to converge the requested number of eigenpairs...
const SparseMatrix< Number > & get_matrix(const std::string &mat_name) const
virtual bool side_residual(bool get_jacobian, DiffContext &) libmesh_override
Forms the jacobian of the boundary terms.
virtual ~EigenTimeSolver()
Destructor.
NowAssembling now_assembling
Flag which controls the internals of element_residual() and side_residual().
OStreamProxy out
SparseMatrix< Number > * matrix
The system matrix.
unsigned int n_basis_vectors_to_use
The number of basis vectors to use in the computation.
Real tol
The linear solver tolerance to be used when solving the eigenvalue problem.
EigenTimeSolver(sys_type &s)
Constructor.
virtual bool side_mass_residual(bool request_jacobian, DiffContext &)
Subtracts a mass vector contribution on side of elem from elem_residual.
Definition: diff_physics.h:334
virtual bool side_time_derivative(bool request_jacobian, DiffContext &)
Adds the time derivative contribution on side of elem to elem_residual.
Definition: diff_physics.h:170
unsigned int n_converged_eigenpairs
After a solve, holds the number of eigenpairs successfully converged.
virtual bool element_constraint(bool request_jacobian, DiffContext &)
Adds the constraint contribution on elem to elem_residual.
Definition: diff_physics.h:141
unsigned int n_eigenpairs_to_compute
The number of eigenvectors/values to be computed.
The matrix associated with the spatial part of the operator.
The matrix associated with the time derivative (mass matrix).
virtual bool nonlocal_mass_residual(bool request_jacobian, DiffContext &c)
Subtracts any nonlocal mass vector contributions (e.g.
Definition: diff_physics.C:63
This class provides an interface to solvers for eigenvalue problems.
Definition: eigen_solver.h:54