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_PETSC_LINEAR_SOLVER_H
21 #define LIBMESH_PETSC_LINEAR_SOLVER_H
22
23 #include "libmesh/libmesh_config.h"
24
25 #ifdef LIBMESH_HAVE_PETSC
26
27 #include "libmesh/petsc_macro.h"
28 #include "libmesh/petsc_solver_exception.h"
29
30 // Petsc include files.
31 #ifdef I
32 # define LIBMESH_SAW_I
33 #endif
34 #include <petscksp.h>
35 #ifndef LIBMESH_SAW_I
36 # undef I // Avoid complex.h contamination
37 #endif
38
39 // Local includes
40 #include "libmesh/linear_solver.h"
41
42 // C++ includes
43 #include <cstddef>
44 #include <vector>
45
46 //--------------------------------------------------------------------
47 // Functions with C linkage to pass to PETSc. PETSc will call these
48 // methods as needed for preconditioning
49 //
50 // Since they must have C linkage they have no knowledge of a namespace.
51 // Give them an obscure name to avoid namespace pollution.
52 extern "C"
53 {
54 /**
55 * This function is called by PETSc to initialize the preconditioner.
56 * ctx will hold the Preconditioner.
57 */
58 PetscErrorCode libmesh_petsc_preconditioner_setup (PC);
59
60 /**
61 * This function is called by PETSc to actually apply the preconditioner.
62 * ctx will hold the Preconditioner.
63 */
64 PetscErrorCode libmesh_petsc_preconditioner_apply(PC, Vec x, Vec y);
65
66 #if LIBMESH_ENABLE_DEPRECATED
67 /**
68 * This function is called by PETSc to initialize the preconditioner.
69 * ctx will hold the Preconditioner.
70 */
71 PetscErrorCode __libmesh_petsc_preconditioner_setup (PC);
72
73 /**
74 * This function is called by PETSc to actually apply the preconditioner.
75 * ctx will hold the Preconditioner.
76 */
77 PetscErrorCode __libmesh_petsc_preconditioner_apply(PC, Vec x, Vec y);
78 #endif
79 } // end extern "C"
80
81
82 namespace libMesh
83 {
84
85 // forward declarations
86 template <typename T> class PetscMatrix;
87
88 /**
89 * This class provides an interface to PETSc
90 * iterative solvers that is compatible with the \p libMesh
91 * \p LinearSolver<>
92 *
93 * \author Benjamin Kirk
94 * \date 2002-2007
95 */
96 template <typename T>
97 class PetscLinearSolver : public LinearSolver<T>
98 {
99 public:
100 /**
101 * Constructor. Initializes Petsc data structures
102 */
103 PetscLinearSolver (const libMesh::Parallel::Communicator & comm_in);
104
105 /**
106 * Destructor.
107 */
108 ~PetscLinearSolver ();
109
110 /**
111 * Release all memory and clear data structures.
112 */
113 virtual void clear () override;
114
115 /**
116 * Initialize data structures if not done so already.
117 * Assigns a name, which is turned into an underscore-separated
118 * prefix for the underlying KSP object.
119 */
120 virtual void init (const char * name = nullptr) override;
121
122 /**
123 * Initialize data structures if not done so already plus much more
124 */
125 void init (PetscMatrix<T> * matrix,
126 const char * name = nullptr);
127
128 /**
129 * Apply names to the system to be solved. This sets an option
130 * prefix from the system name and sets field names from the
131 * system's variable names.
132 *
133 * Since field names are applied to DoF numberings, this method must
134 * be called again after any System reinit.
135 */
136 virtual void init_names (const System &) override;
137
138 /**
139 * After calling this method, all successive solves will be
140 * restricted to the given set of dofs, which must contain local
141 * dofs on each processor only and not contain any duplicates. This
142 * mode can be disabled by calling this method with \p dofs being a
143 * \p nullptr.
144 */
145 virtual void restrict_solve_to (const std::vector<unsigned int> * const dofs,
146 const SubsetSolveMode subset_solve_mode=SUBSET_ZERO) override;
147
148 /**
149 * Call the Petsc solver. It calls the method below, using the
150 * same matrix for the system and preconditioner matrices.
151 */
152 virtual std::pair<unsigned int, Real>
solve(SparseMatrix<T> & matrix_in,NumericVector<T> & solution_in,NumericVector<T> & rhs_in,const double tol,const unsigned int m_its)153 solve (SparseMatrix<T> & matrix_in,
154 NumericVector<T> & solution_in,
155 NumericVector<T> & rhs_in,
156 const double tol,
157 const unsigned int m_its) override
158 {
159 return this->solve(matrix_in, matrix_in, solution_in, rhs_in, tol, m_its);
160 }
161
162
163 /**
164 * Call the Petsc solver. It calls the method below, using the
165 * same matrix for the system and preconditioner matrices.
166 */
167 virtual std::pair<unsigned int, Real>
168 adjoint_solve (SparseMatrix<T> & matrix_in,
169 NumericVector<T> & solution_in,
170 NumericVector<T> & rhs_in,
171 const double tol,
172 const unsigned int m_its) override;
173
174 /**
175 * This method allows you to call a linear solver while specifying
176 * the matrix to use as the (left) preconditioning matrix.
177 *
178 * \note The linear solver will not compute a preconditioner in this
179 * case, and will instead premultiply by the matrix you provide. In
180 * PETSc, this is accomplished by calling
181 * \code
182 * PCSetType(_pc, PCMAT);
183 * \endcode
184 * before invoking KSPSolve().
185 *
186 * \note This functionality is not implemented in the LinearSolver
187 * class since there is not a built-in analog to this method for
188 * LASPACK. You could probably implement it by hand if you wanted.
189 */
190 virtual std::pair<unsigned int, Real>
191 solve (SparseMatrix<T> & matrix,
192 SparseMatrix<T> & preconditioner,
193 NumericVector<T> & solution,
194 NumericVector<T> & rhs,
195 const double tol,
196 const unsigned int m_its) override;
197
198 /**
199 * This function solves a system whose matrix is a shell matrix.
200 */
201 virtual std::pair<unsigned int, Real>
202 solve (const ShellMatrix<T> & shell_matrix,
203 NumericVector<T> & solution_in,
204 NumericVector<T> & rhs_in,
205 const double tol,
206 const unsigned int m_its) override;
207
208 /**
209 * This function solves a system whose matrix is a shell matrix, but
210 * a sparse matrix is used as preconditioning matrix, this allowing
211 * other preconditioners than JACOBI.
212 */
213 virtual std::pair<unsigned int, Real>
214 solve (const ShellMatrix<T> & shell_matrix,
215 const SparseMatrix<T> & precond_matrix,
216 NumericVector<T> & solution_in,
217 NumericVector<T> & rhs_in,
218 const double tol,
219 const unsigned int m_its) override;
220
221 /**
222 * \returns The raw PETSc preconditioner context pointer.
223 *
224 * This allows you to specify the PCShellSetApply() and
225 * PCShellSetSetUp() functions if you desire. Just don't do
226 * anything crazy like calling libMeshPCDestroy() on the pointer!
227 */
pc()228 PC pc() { this->init(); return _pc; }
229
230 /**
231 * \returns The raw PETSc ksp context pointer.
232 *
233 * This is useful if you are for example setting a custom
234 * convergence test with KSPSetConvergenceTest().
235 */
ksp()236 KSP ksp() { this->init(); return _ksp; }
237
238 /**
239 * Fills the input vector with the sequence of residual norms
240 * from the latest iterative solve.
241 */
242 void get_residual_history(std::vector<double> & hist);
243
244 /**
245 * \returns Just the initial residual for the solve just
246 * completed with this interface.
247 *
248 * Use this method instead of the one above if you just want the
249 * starting residual and not the entire history.
250 */
251 Real get_initial_residual();
252
253 /**
254 * \returns The solver's convergence flag
255 */
256 virtual LinearConvergenceReason get_converged_reason() const override;
257
258 private:
259
260 /**
261 * Tells PETSC to use the user-specified solver stored in
262 * \p _solver_type
263 */
264 void set_petsc_solver_type ();
265
266 /**
267 * Internal function if shell matrix mode is used.
268 */
269 static PetscErrorCode _petsc_shell_matrix_mult(Mat mat, Vec arg, Vec dest);
270
271 /**
272 * Internal function if shell matrix mode is used.
273 */
274 static PetscErrorCode _petsc_shell_matrix_mult_add(Mat mat, Vec arg, Vec add, Vec dest);
275
276 /**
277 * Internal function if shell matrix mode is used.
278 */
279 static PetscErrorCode _petsc_shell_matrix_get_diagonal(Mat mat, Vec dest);
280
281 /**
282 * Preconditioner context
283 */
284 PC _pc;
285
286 /**
287 * Krylov subspace context
288 */
289 KSP _ksp;
290
291 /**
292 * PETSc index set containing the dofs on which to solve (\p nullptr
293 * means solve on all dofs).
294 */
295 IS _restrict_solve_to_is;
296
297 /**
298 * PETSc index set, complement to \p _restrict_solve_to_is. This
299 * will be created on demand by the method \p
300 * _create_complement_is().
301 */
302 IS _restrict_solve_to_is_complement;
303
304 /**
305 * \returns The local size of \p _restrict_solve_to_is.
306 */
307 PetscInt _restrict_solve_to_is_local_size() const;
308
309 /**
310 * Creates \p _restrict_solve_to_is_complement to contain all
311 * indices that are local in \p vec_in, except those that are
312 * contained in \p _restrict_solve_to_is.
313 */
314 void _create_complement_is (const NumericVector<T> & vec_in);
315
316 /**
317 * If restrict-solve-to-subset mode is active, this member decides
318 * what happens with the dofs outside the subset.
319 */
320 SubsetSolveMode _subset_solve_mode;
321 };
322
323
324 /*----------------------- functions ----------------------------------*/
325 template <typename T>
326 inline
~PetscLinearSolver()327 PetscLinearSolver<T>::~PetscLinearSolver ()
328 {
329 this->clear ();
330 }
331
332
333
334 template <typename T>
335 inline PetscInt
_restrict_solve_to_is_local_size()336 PetscLinearSolver<T>::_restrict_solve_to_is_local_size() const
337 {
338 libmesh_assert(_restrict_solve_to_is);
339
340 PetscInt s;
341 int ierr = ISGetLocalSize(_restrict_solve_to_is, &s);
342 LIBMESH_CHKERR(ierr);
343
344 return s;
345 }
346
347
348
349 template <typename T>
350 void
_create_complement_is(const NumericVector<T> & vec_in)351 PetscLinearSolver<T>::_create_complement_is (const NumericVector<T> & vec_in)
352 {
353 libmesh_assert(_restrict_solve_to_is);
354 if (_restrict_solve_to_is_complement==nullptr)
355 {
356 int ierr = ISComplement(_restrict_solve_to_is,
357 vec_in.first_local_index(),
358 vec_in.last_local_index(),
359 &_restrict_solve_to_is_complement);
360 LIBMESH_CHKERR(ierr);
361 }
362 }
363
364
365
366 } // namespace libMesh
367
368
369 #endif // #ifdef LIBMESH_HAVE_PETSC
370 #endif // LIBMESH_PETSC_LINEAR_SOLVER_H
371