libMesh
petsc_linear_solver.h
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 #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 #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 } // end extern "C"
65 
66 
67 namespace libMesh
68 {
69 
70 // forward declarations
71 template <typename T> class PetscMatrix;
72 
81 template <typename T>
83 {
84 public:
89  LIBMESH_CAN_DEFAULT_TO_COMMWORLD);
90 
95 
99  virtual void clear () libmesh_override;
100 
106  virtual void init (const char * name = libmesh_nullptr) libmesh_override;
107 
111  void init (PetscMatrix<T> * matrix,
112  const char * name = libmesh_nullptr);
113 
122  virtual void init_names (const System &) libmesh_override;
123 
131  virtual void restrict_solve_to (const std::vector<unsigned int> * const dofs,
132  const SubsetSolveMode subset_solve_mode=SUBSET_ZERO) libmesh_override;
133 
138  virtual std::pair<unsigned int, Real>
139  solve (SparseMatrix<T> & matrix_in,
140  NumericVector<T> & solution_in,
141  NumericVector<T> & rhs_in,
142  const double tol,
143  const unsigned int m_its) libmesh_override
144  {
145  return this->solve(matrix_in, matrix_in, solution_in, rhs_in, tol, m_its);
146  }
147 
148 
153  virtual std::pair<unsigned int, Real>
154  adjoint_solve (SparseMatrix<T> & matrix_in,
155  NumericVector<T> & solution_in,
156  NumericVector<T> & rhs_in,
157  const double tol,
158  const unsigned int m_its) libmesh_override;
159 
176  virtual std::pair<unsigned int, Real>
177  solve (SparseMatrix<T> & matrix,
178  SparseMatrix<T> & preconditioner,
179  NumericVector<T> & solution,
180  NumericVector<T> & rhs,
181  const double tol,
182  const unsigned int m_its) libmesh_override;
183 
187  virtual std::pair<unsigned int, Real>
188  solve (const ShellMatrix<T> & shell_matrix,
189  NumericVector<T> & solution_in,
190  NumericVector<T> & rhs_in,
191  const double tol,
192  const unsigned int m_its) libmesh_override;
193 
199  virtual std::pair<unsigned int, Real>
200  solve (const ShellMatrix<T> & shell_matrix,
201  const SparseMatrix<T> & precond_matrix,
202  NumericVector<T> & solution_in,
203  NumericVector<T> & rhs_in,
204  const double tol,
205  const unsigned int m_its) libmesh_override;
206 
214  PC pc() { this->init(); return _pc; }
215 
222  KSP ksp() { this->init(); return _ksp; }
223 
228  void get_residual_history(std::vector<double> & hist);
229 
238 
242  virtual LinearConvergenceReason get_converged_reason() const libmesh_override;
243 
244 private:
245 
250  void set_petsc_solver_type ();
251 
255  static PetscErrorCode _petsc_shell_matrix_mult(Mat mat, Vec arg, Vec dest);
256 
260  static PetscErrorCode _petsc_shell_matrix_mult_add(Mat mat, Vec arg, Vec add, Vec dest);
261 
265  static PetscErrorCode _petsc_shell_matrix_get_diagonal(Mat mat, Vec dest);
266 
270  PC _pc;
271 
275  KSP _ksp;
276 
282 
289 
293  PetscInt _restrict_solve_to_is_local_size() const;
294 
300  void _create_complement_is (const NumericVector<T> & vec_in);
301 
307 };
308 
309 
310 /*----------------------- functions ----------------------------------*/
311 template <typename T>
312 inline
314  LinearSolver<T>(comm_in),
318 {
319  if (this->n_processors() == 1)
321  else
323 }
324 
325 
326 
327 template <typename T>
328 inline
330 {
331  this->clear ();
332 }
333 
334 
335 
336 template <typename T>
337 inline PetscInt
339 {
341 
342  PetscInt s;
343  int ierr = ISGetLocalSize(_restrict_solve_to_is, &s);
344  LIBMESH_CHKERR(ierr);
345 
346  return s;
347 }
348 
349 
350 
351 template <typename T>
352 void
354 #if PETSC_VERSION_LESS_THAN(3,0,0)
355  // unnamed to avoid compiler "unused parameter" warning
356 #else
357  vec_in
358 #endif
359  )
360 {
362 #if PETSC_VERSION_LESS_THAN(3,0,0)
363  // No ISComplement in PETSc 2.3.3
364  libmesh_not_implemented();
365 #else
367  {
368  int ierr = ISComplement(_restrict_solve_to_is,
369  vec_in.first_local_index(),
370  vec_in.last_local_index(),
372  LIBMESH_CHKERR(ierr);
373  }
374 #endif
375 }
376 
377 
378 
379 } // namespace libMesh
380 
381 
382 #endif // #ifdef LIBMESH_HAVE_PETSC
383 #endif // LIBMESH_PETSC_LINEAR_SOLVER_H
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:39
void get_residual_history(std::vector< double > &hist)
Fills the input vector with the sequence of residual norms from the latest iterative solve...
Encapsulates the MPI_Comm object.
Definition: parallel.h:657
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) libmesh_override
Call the Petsc solver.
virtual void clear() libmesh_override
Release all memory and clear data structures.
virtual LinearConvergenceReason get_converged_reason() const libmesh_override
Set dofs outside the subset to zero.
processor_id_type n_processors() const
void _create_complement_is(const NumericVector< T > &vec_in)
Creates _restrict_solve_to_is_complement to contain all indices that are local in vec_in...
PetscInt _restrict_solve_to_is_local_size() const
const class libmesh_nullptr_t libmesh_nullptr
static PetscErrorCode _petsc_shell_matrix_mult_add(Mat mat, Vec arg, Vec add, Vec dest)
Internal function if shell matrix mode is used.
Numeric vector.
Definition: dof_map.h:66
The libMesh namespace provides an interface to certain functionality in the library.
libmesh_assert(j)
Generic sparse matrix.
Definition: dof_map.h:65
This base class can be inherited from to provide interfaces to linear solvers from different packages...
Definition: linear_solver.h:58
virtual void init_names(const System &) libmesh_override
Apply names to the system to be solved.
PetscErrorCode __libmesh_petsc_preconditioner_apply(void *ctx, Vec x, Vec y)
This function is called by PETSc to actually apply the preconditioner.
PetscErrorCode Vec x
PetscLinearSolver(const libMesh::Parallel::Communicator &comm_in LIBMESH_CAN_DEFAULT_TO_COMMWORLD)
Constructor.
IS _restrict_solve_to_is
PETSc index set containing the dofs on which to solve (NULL means solve on all dofs).
This is the base class for classes which contain information related to any physical process that mig...
Definition: system.h:76
void set_petsc_solver_type()
Tells PETSC to use the user-specified solver stored in _solver_type.
PetscErrorCode __libmesh_petsc_preconditioner_setup(void *ctx)
This function is called by PETSc to initialize the preconditioner.
PetscErrorCode Vec Mat Mat void * ctx
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) libmesh_override
Call the Petsc solver.
virtual void init(const char *name=libmesh_nullptr) libmesh_override
Initialize data structures if not done so already.
PC _pc
Preconditioner context.
static PetscErrorCode _petsc_shell_matrix_mult(Mat mat, Vec arg, Vec dest)
Internal function if shell matrix mode is used.
IS _restrict_solve_to_is_complement
PETSc index set, complement to _restrict_solve_to_is.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
PetscErrorCode ierr
This class provides an interface to PETSc iterative solvers that is compatible with the libMesh Linea...
PreconditionerType _preconditioner_type
Enum stating with type of preconditioner to use.
LinearConvergenceReason
Linear solver convergence flags (taken from the PETSc flags)
SubsetSolveMode _subset_solve_mode
If restrict-solve-to-subset mode is active, this member decides what happens with the dofs outside th...
SubsetSolveMode
defines an enum for the question what happens to the dofs outside the given subset when a system is s...
KSP _ksp
Krylov subspace context.
Generic shell matrix, i.e.
static PetscErrorCode _petsc_shell_matrix_get_diagonal(Mat mat, Vec dest)
Internal function if shell matrix mode is used.
virtual void restrict_solve_to(const std::vector< unsigned int > *const dofs, const SubsetSolveMode subset_solve_mode=SUBSET_ZERO) libmesh_override
After calling this method, all successive solves will be restricted to the given set of dofs...