petsc_linear_solver.h
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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"
29 
30 // Petsc include files.
31 #include <petscksp.h>
32 
33 // Local includes
34 #include "libmesh/linear_solver.h"
35 
36 // C++ includes
37 #include <cstddef>
38 #include <vector>
39 
40 //--------------------------------------------------------------------
41 // Functions with C linkage to pass to PETSc. PETSc will call these
42 // methods as needed for preconditioning
43 //
44 // Since they must have C linkage they have no knowledge of a namespace.
45 // Give them an obscure name to avoid namespace pollution.
46 extern "C"
47 {
48 #if PETSC_RELEASE_LESS_THAN(3,0,1)
49 
53  PetscErrorCode libmesh_petsc_preconditioner_setup (void * ctx);
54 
59  PetscErrorCode libmesh_petsc_preconditioner_apply(void * ctx, Vec x, Vec y);
60 #else
61  PetscErrorCode libmesh_petsc_preconditioner_setup (PC);
62  PetscErrorCode libmesh_petsc_preconditioner_apply(PC, Vec x, Vec y);
63 #endif
64 
65 #if LIBMESH_ENABLE_DEPRECATED
66 #if PETSC_RELEASE_LESS_THAN(3,0,1)
67 
71  PetscErrorCode __libmesh_petsc_preconditioner_setup (void * ctx);
72 
77  PetscErrorCode __libmesh_petsc_preconditioner_apply(void * ctx, Vec x, Vec y);
78 #else
79  PetscErrorCode __libmesh_petsc_preconditioner_setup (PC);
80  PetscErrorCode __libmesh_petsc_preconditioner_apply(PC, Vec x, Vec y);
81 #endif
82 #endif
83 } // end extern "C"
84 
85 
86 namespace libMesh
87 {
88 
89 // forward declarations
90 template <typename T> class PetscMatrix;
91 
100 template <typename T>
102 {
103 public:
108 
113 
117  virtual void clear () override;
118 
124  virtual void init (const char * name = nullptr) override;
125 
129  void init (PetscMatrix<T> * matrix,
130  const char * name = nullptr);
131 
140  virtual void init_names (const System &) override;
141 
149  virtual void restrict_solve_to (const std::vector<unsigned int> * const dofs,
150  const SubsetSolveMode subset_solve_mode=SUBSET_ZERO) override;
151 
156  virtual std::pair<unsigned int, Real>
157  solve (SparseMatrix<T> & matrix_in,
158  NumericVector<T> & solution_in,
159  NumericVector<T> & rhs_in,
160  const double tol,
161  const unsigned int m_its) override
162  {
163  return this->solve(matrix_in, matrix_in, solution_in, rhs_in, tol, m_its);
164  }
165 
166 
171  virtual std::pair<unsigned int, Real>
172  adjoint_solve (SparseMatrix<T> & matrix_in,
173  NumericVector<T> & solution_in,
174  NumericVector<T> & rhs_in,
175  const double tol,
176  const unsigned int m_its) override;
177 
194  virtual std::pair<unsigned int, Real>
195  solve (SparseMatrix<T> & matrix,
196  SparseMatrix<T> & preconditioner,
197  NumericVector<T> & solution,
198  NumericVector<T> & rhs,
199  const double tol,
200  const unsigned int m_its) override;
201 
205  virtual std::pair<unsigned int, Real>
206  solve (const ShellMatrix<T> & shell_matrix,
207  NumericVector<T> & solution_in,
208  NumericVector<T> & rhs_in,
209  const double tol,
210  const unsigned int m_its) override;
211 
217  virtual std::pair<unsigned int, Real>
218  solve (const ShellMatrix<T> & shell_matrix,
219  const SparseMatrix<T> & precond_matrix,
220  NumericVector<T> & solution_in,
221  NumericVector<T> & rhs_in,
222  const double tol,
223  const unsigned int m_its) override;
224 
232  PC pc() { this->init(); return _pc; }
233 
240  KSP ksp() { this->init(); return _ksp; }
241 
246  void get_residual_history(std::vector<double> & hist);
247 
256 
260  virtual LinearConvergenceReason get_converged_reason() const override;
261 
262 private:
263 
268  void set_petsc_solver_type ();
269 
273  static PetscErrorCode _petsc_shell_matrix_mult(Mat mat, Vec arg, Vec dest);
274 
278  static PetscErrorCode _petsc_shell_matrix_mult_add(Mat mat, Vec arg, Vec add, Vec dest);
279 
283  static PetscErrorCode _petsc_shell_matrix_get_diagonal(Mat mat, Vec dest);
284 
288  PC _pc;
289 
293  KSP _ksp;
294 
300 
307 
311  PetscInt _restrict_solve_to_is_local_size() const;
312 
318  void _create_complement_is (const NumericVector<T> & vec_in);
319 
325 };
326 
327 
328 /*----------------------- functions ----------------------------------*/
329 template <typename T>
330 inline
332 {
333  this->clear ();
334 }
335 
336 
337 
338 template <typename T>
339 inline PetscInt
341 {
342  libmesh_assert(_restrict_solve_to_is);
343 
344  PetscInt s;
345  int ierr = ISGetLocalSize(_restrict_solve_to_is, &s);
346  LIBMESH_CHKERR(ierr);
347 
348  return s;
349 }
350 
351 
352 
353 template <typename T>
354 void
356 #if PETSC_VERSION_LESS_THAN(3,0,0)
357  // unnamed to avoid compiler "unused parameter" warning
358 #else
359  vec_in
360 #endif
361  )
362 {
363  libmesh_assert(_restrict_solve_to_is);
364 #if PETSC_VERSION_LESS_THAN(3,0,0)
365  // No ISComplement in PETSc 2.3.3
366  libmesh_not_implemented();
367 #else
368  if (_restrict_solve_to_is_complement==nullptr)
369  {
370  int ierr = ISComplement(_restrict_solve_to_is,
371  vec_in.first_local_index(),
372  vec_in.last_local_index(),
373  &_restrict_solve_to_is_complement);
374  LIBMESH_CHKERR(ierr);
375  }
376 #endif
377 }
378 
379 
380 
381 } // namespace libMesh
382 
383 
384 #endif // #ifdef LIBMESH_HAVE_PETSC
385 #endif // LIBMESH_PETSC_LINEAR_SOLVER_H
std::string name(const ElemQuality q)
Definition: elem_quality.C:42
void get_residual_history(std::vector< double > &hist)
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) override
virtual void init(const char *name=nullptr) override
void _create_complement_is(const NumericVector< T > &vec_in)
static PetscErrorCode _petsc_shell_matrix_mult_add(Mat mat, Vec arg, Vec add, Vec dest)
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
Definition: diff_context.h:40
PetscErrorCode __libmesh_petsc_preconditioner_apply(void *ctx, Vec x, Vec y)
PetscLinearSolver(const libMesh::Parallel::Communicator &comm_in)
Manages consistently variables, degrees of freedom, and coefficient vectors.
Definition: system.h:92
PetscErrorCode __libmesh_petsc_preconditioner_setup(void *ctx)
virtual void restrict_solve_to(const std::vector< unsigned int > *const dofs, const SubsetSolveMode subset_solve_mode=SUBSET_ZERO) override
virtual std::pair< unsigned int, Real > adjoint_solve(SparseMatrix< T > &matrix_in, NumericVector< T > &solution_in, NumericVector< T > &rhs_in, const double tol, const unsigned int m_its) override
virtual void init_names(const System &) override
PetscErrorCode libmesh_petsc_preconditioner_apply(void *ctx, Vec x, Vec y)
static PetscErrorCode _petsc_shell_matrix_mult(Mat mat, Vec arg, Vec dest)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
PetscErrorCode ierr
PetscErrorCode libmesh_petsc_preconditioner_setup(void *ctx)
SparseMatrix interface to PETSc Mat.
PetscInt _restrict_solve_to_is_local_size() const
virtual LinearConvergenceReason get_converged_reason() const override
static PetscErrorCode _petsc_shell_matrix_get_diagonal(Mat mat, Vec dest)
virtual void clear() override