eigen_time_solver.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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"
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  // Assemble the operator for the spatial part.
134  if (now_assembling == Matrix_A)
135  {
136  bool jacobian_computed =
137  _system.element_time_derivative(request_jacobian, context);
138 
139  // The user shouldn't compute a jacobian unless requested
140  libmesh_assert(request_jacobian || !jacobian_computed);
141 
142  bool jacobian_computed2 =
143  _system.element_constraint(jacobian_computed, context);
144 
145  // The user shouldn't compute a jacobian unless requested
146  libmesh_assert (jacobian_computed || !jacobian_computed2);
147 
148  return jacobian_computed && jacobian_computed2;
149 
150  }
151 
152  // Assemble the mass matrix operator
153  else if (now_assembling == Matrix_B)
154  {
155  bool mass_jacobian_computed =
156  _system.mass_residual(request_jacobian, context);
157 
158  // Scale Jacobian by -1?
159  //context.elem_jacobian *= -1.0;
160 
161  return mass_jacobian_computed;
162  }
163 
164  else
165  {
166  libmesh_error();
167  return false;
168  }
169 }
170 
171 
172 
173 bool EigenTimeSolver::side_residual(bool request_jacobian,
174  DiffContext &context)
175 {
176  // The EigenTimeSolver always requests jacobians?
177  //libmesh_assert (request_jacobian);
178 
179  // Assemble the operator for the spatial part.
180  if (now_assembling == Matrix_A)
181  {
182  bool jacobian_computed =
183  _system.side_time_derivative(request_jacobian, context);
184 
185  // The user shouldn't compute a jacobian unless requested
186  libmesh_assert (request_jacobian || !jacobian_computed);
187 
188  bool jacobian_computed2 =
189  _system.side_constraint(jacobian_computed, context);
190 
191  // The user shouldn't compute a jacobian unless requested
192  libmesh_assert (jacobian_computed || !jacobian_computed2);
193 
194  return jacobian_computed && jacobian_computed2;
195 
196  }
197 
198  // There is now a "side" equivalent for the mass matrix
199  else if (now_assembling == Matrix_B)
200  {
201  bool mass_jacobian_computed =
202  _system.side_mass_residual(request_jacobian, context);
203 
204  return mass_jacobian_computed;
205  }
206 
207  else
208  {
209  libmesh_error();
210  return false;
211  }
212 }
213 
214 } // namespace libMesh
215 
216 #endif // LIBMESH_HAVE_SLEPC

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:04 UTC

Hosted By:
SourceForge.net Logo