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_LINEAR_SOLVER_H
21 #define LIBMESH_LINEAR_SOLVER_H
22
23 // Local includes
24 #include "libmesh/libmesh_common.h"
25 #include "libmesh/enum_subset_solve_mode.h" // SUBSET_ZERO
26 #include "libmesh/reference_counted_object.h"
27 #include "libmesh/libmesh.h"
28 #include "libmesh/parallel_object.h"
29
30 #ifdef LIBMESH_FORWARD_DECLARE_ENUMS
31 namespace libMesh
32 {
33 enum SolverPackage : int;
34 enum PreconditionerType : int;
35 enum SolverType : int;
36 enum LinearConvergenceReason : int;
37 }
38 #else
39 #include "libmesh/enum_solver_package.h"
40 #include "libmesh/enum_preconditioner_type.h"
41 #include "libmesh/enum_solver_type.h"
42 #include "libmesh/enum_convergence_flags.h"
43 #endif
44
45 // C++ includes
46 #include <cstddef>
47 #include <vector>
48 #include <memory>
49
50 namespace libMesh
51 {
52
53 // forward declarations
54 template <typename T> class SparseMatrix;
55 template <typename T> class NumericVector;
56 template <typename T> class ShellMatrix;
57 template <typename T> class Preconditioner;
58 class System;
59 class SolverConfiguration;
60
61 /**
62 * This base class can be inherited from to provide interfaces to
63 * linear solvers from different packages like PETSc and LASPACK.
64 *
65 * \author Benjamin Kirk
66 * \date 2003
67 */
68 template <typename T>
69 class LinearSolver : public ReferenceCountedObject<LinearSolver<T>>,
70 public ParallelObject
71 {
72 public:
73
74 /**
75 * Constructor. Initializes Solver data structures
76 */
77 LinearSolver (const libMesh::Parallel::Communicator & comm_in);
78
79 /**
80 * Destructor.
81 */
82 virtual ~LinearSolver ();
83
84 /**
85 * Builds a \p LinearSolver using the linear solver package specified by
86 * \p solver_package
87 */
88 static std::unique_ptr<LinearSolver<T>> build(const libMesh::Parallel::Communicator & comm_in,
89 const SolverPackage solver_package = libMesh::default_solver_package());
90
91 /**
92 * \returns \p true if the data structures are
93 * initialized, false otherwise.
94 */
initialized()95 bool initialized () const { return _is_initialized; }
96
97 /**
98 * Release all memory and clear data structures.
99 */
clear()100 virtual void clear () {}
101
102 /**
103 * Initialize data structures if not done so already.
104 * May assign a name to the solver in some implementations
105 */
106 virtual void init (const char * name = nullptr) = 0;
107
108 /**
109 * Apply names to the system to be solved. For most packages this
110 * is a no-op; for PETSc this sets an option prefix from the system
111 * name and sets field names from the system's variable names.
112 *
113 * Since field names are applied to DoF numberings, this method must
114 * be called again after any System reinit.
115 */
init_names(const System &)116 virtual void init_names (const System &) {}
117
118 /**
119 * \returns The type of solver to use.
120 */
solver_type()121 SolverType solver_type () const { return _solver_type; }
122
123 /**
124 * Sets the type of solver to use.
125 */
set_solver_type(const SolverType st)126 void set_solver_type (const SolverType st)
127 { _solver_type = st; }
128
129 /**
130 * \returns The type of preconditioner to use.
131 */
132 PreconditionerType preconditioner_type () const;
133
134 /**
135 * Sets the type of preconditioner to use.
136 */
137 void set_preconditioner_type (const PreconditionerType pct);
138
139 /**
140 * Attaches a Preconditioner object to be used
141 */
142 void attach_preconditioner(Preconditioner<T> * preconditioner);
143
144 /**
145 * Set the same_preconditioner flag, which indicates if we reuse the
146 * same preconditioner for subsequent solves.
147 */
148 virtual void reuse_preconditioner(bool );
149
150 /**
151 * \returns \p same_preconditioner, which indicates if we reuse the
152 * same preconditioner for subsequent solves.
153 */
154 bool get_same_preconditioner();
155
156 /**
157 * After calling this method, all successive solves will be
158 * restricted to the given set of dofs, which must contain local
159 * dofs on each processor only and not contain any duplicates. This
160 * mode can be disabled by calling this method with \p dofs being a
161 * \p nullptr.
162 */
163 virtual void restrict_solve_to (const std::vector<unsigned int> * const dofs,
164 const SubsetSolveMode subset_solve_mode=SUBSET_ZERO);
165
166 /**
167 * This function calls the solver \p _solver_type preconditioned
168 * with the \p _preconditioner_type preconditioner.
169 *
170 * \note This method will compute the preconditioner from the system
171 * matrix.
172 */
173 virtual std::pair<unsigned int, Real> solve (SparseMatrix<T> &, // System Matrix
174 NumericVector<T> &, // Solution vector
175 NumericVector<T> &, // RHS vector
176 const double, // Stopping tolerance
177 const unsigned int) = 0; // N. Iterations
178
179 /**
180 * Function to solve the adjoint system.
181 *
182 * \note This method will compute the preconditioner from the system
183 * matrix. This is not a pure virtual function and is defined
184 * linear_solver.C
185 */
186 virtual std::pair<unsigned int, Real> adjoint_solve (SparseMatrix<T> &, // System Matrix
187 NumericVector<T> &, // Solution vector
188 NumericVector<T> &, // RHS vector
189 const double, // Stopping tolerance
190 const unsigned int); // N. Iterations
191
192 /**
193 * This function calls the solver
194 * "_solver_type" preconditioned with the
195 * "_preconditioner_type" preconditioner.
196 */
197 virtual std::pair<unsigned int, Real> solve (SparseMatrix<T> &, // System Matrix
198 SparseMatrix<T> &, // Preconditioning Matrix
199 NumericVector<T> &, // Solution vector
200 NumericVector<T> &, // RHS vector
201 const double, // Stopping tolerance
202 const unsigned int) = 0; // N. Iterations
203
204 /**
205 * This function calls the solver "_solver_type" preconditioned with
206 * the "_preconditioner_type" preconditioner. The preconditioning
207 * matrix is used if it is provided, or the system matrix is used if
208 * \p precond_matrix is null
209 */
210 std::pair<unsigned int, Real> solve (SparseMatrix<T> & matrix,
211 SparseMatrix<T> * precond_matrix,
212 NumericVector<T> &, // Solution vector
213 NumericVector<T> &, // RHS vector
214 const double, // Stopping tolerance
215 const unsigned int); // N. Iterations
216
217
218
219 /**
220 * This function solves a system whose matrix is a shell matrix.
221 */
222 virtual std::pair<unsigned int, Real> solve (const ShellMatrix<T> & shell_matrix,
223 NumericVector<T> &, // Solution vector
224 NumericVector<T> &, // RHS vector
225 const double, // Stopping tolerance
226 const unsigned int) = 0; // N. Iterations
227
228
229 /**
230 * This function solves a system whose matrix is a shell matrix, but
231 * a sparse matrix is used as preconditioning matrix, this allowing
232 * other preconditioners than JACOBI.
233 */
234 virtual std::pair<unsigned int, Real> solve (const ShellMatrix<T> & shell_matrix,
235 const SparseMatrix<T> & precond_matrix,
236 NumericVector<T> &, // Solution vector
237 NumericVector<T> &, // RHS vector
238 const double, // Stopping tolerance
239 const unsigned int) = 0; // N. Iterations
240
241
242 /**
243 * This function solves a system whose matrix is a shell matrix, but
244 * an optional sparse matrix may be used as preconditioning matrix.
245 */
246 std::pair<unsigned int, Real> solve (const ShellMatrix<T> & matrix,
247 const SparseMatrix<T> * precond_matrix,
248 NumericVector<T> &, // Solution vector
249 NumericVector<T> &, // RHS vector
250 const double, // Stopping tolerance
251 const unsigned int); // N. Iterations
252
253
254 /**
255 * Prints a useful message about why the latest linear solve
256 * con(di)verged.
257 */
258 virtual void print_converged_reason() const;
259
260 /**
261 * \returns The solver's convergence flag
262 */
263 virtual LinearConvergenceReason get_converged_reason() const = 0;
264
265 /**
266 * Set the solver configuration object.
267 */
268 void set_solver_configuration(SolverConfiguration & solver_configuration);
269
270 protected:
271
272
273 /**
274 * Enum stating which type of iterative solver to use.
275 */
276 SolverType _solver_type;
277
278 /**
279 * Enum stating with type of preconditioner to use.
280 */
281 PreconditionerType _preconditioner_type;
282
283 /**
284 * Flag indicating if the data structures have been initialized.
285 */
286 bool _is_initialized;
287
288 /**
289 * Holds the Preconditioner object to be used for the linear solves.
290 */
291 Preconditioner<T> * _preconditioner;
292
293 /**
294 * Boolean flag to indicate whether we want to use an identical
295 * preconditioner to the previous solve. This can save
296 * substantial work in the cases where the system matrix is
297 * the same for successive solves.
298 */
299 bool same_preconditioner;
300
301 /**
302 * Optionally store a SolverOptions object that can be used
303 * to set parameters like solver type, tolerances and iteration limits.
304 */
305 SolverConfiguration * _solver_configuration;
306 };
307
308
309
310
311 /*----------------------- inline functions ----------------------------------*/
312 template <typename T>
313 inline
~LinearSolver()314 LinearSolver<T>::~LinearSolver ()
315 {
316 this->clear ();
317 }
318
319 template <typename T>
320 inline
get_same_preconditioner()321 bool LinearSolver<T>::get_same_preconditioner()
322 {
323 return same_preconditioner;
324 }
325
326 template <typename T>
327 inline
328 std::pair<unsigned int, Real>
solve(SparseMatrix<T> & mat,SparseMatrix<T> * pc_mat,NumericVector<T> & sol,NumericVector<T> & rhs,const double tol,const unsigned int n_iter)329 LinearSolver<T>::solve (SparseMatrix<T> & mat,
330 SparseMatrix<T> * pc_mat,
331 NumericVector<T> & sol,
332 NumericVector<T> & rhs,
333 const double tol,
334 const unsigned int n_iter)
335 {
336 if (pc_mat)
337 return this->solve(mat, *pc_mat, sol, rhs, tol, n_iter);
338 else
339 return this->solve(mat, sol, rhs, tol, n_iter);
340 }
341
342
343 template <typename T>
344 inline
345 std::pair<unsigned int, Real>
solve(const ShellMatrix<T> & mat,const SparseMatrix<T> * pc_mat,NumericVector<T> & sol,NumericVector<T> & rhs,const double tol,const unsigned int n_iter)346 LinearSolver<T>::solve (const ShellMatrix<T> & mat,
347 const SparseMatrix<T> * pc_mat,
348 NumericVector<T> & sol,
349 NumericVector<T> & rhs,
350 const double tol,
351 const unsigned int n_iter)
352 {
353 if (pc_mat)
354 return this->solve(mat, *pc_mat, sol, rhs, tol, n_iter);
355 else
356 return this->solve(mat, sol, rhs, tol, n_iter);
357 }
358
359 } // namespace libMesh
360
361
362 #endif // LIBMESH_LINEAR_SOLVER_H
363