libMesh::PetscLinearSolver< T > Class Template Reference

#include <petsc_linear_solver.h>

Inheritance diagram for libMesh::PetscLinearSolver< T >:

Public Member Functions

 PetscLinearSolver (const libMesh::Parallel::Communicator &comm_in LIBMESH_CAN_DEFAULT_TO_COMMWORLD)
 
 ~PetscLinearSolver ()
 
virtual void clear () libmesh_override
 
virtual void init (const char *name=libmesh_nullptr) libmesh_override
 
void init (PetscMatrix< T > *matrix, const char *name=libmesh_nullptr)
 
virtual void init_names (const System &) libmesh_override
 
virtual void restrict_solve_to (const std::vector< unsigned int > *const dofs, const SubsetSolveMode subset_solve_mode=SUBSET_ZERO) libmesh_override
 
virtual std::pair< unsigned int, Realsolve (SparseMatrix< T > &matrix_in, NumericVector< T > &solution_in, NumericVector< T > &rhs_in, const double tol, const unsigned int m_its) libmesh_override
 
virtual std::pair< unsigned int, Realadjoint_solve (SparseMatrix< T > &matrix_in, NumericVector< T > &solution_in, NumericVector< T > &rhs_in, const double tol, const unsigned int m_its) libmesh_override
 
virtual std::pair< unsigned int, Realsolve (SparseMatrix< T > &matrix, SparseMatrix< T > &preconditioner, NumericVector< T > &solution, NumericVector< T > &rhs, const double tol, const unsigned int m_its) libmesh_override
 
virtual std::pair< unsigned int, Realsolve (const ShellMatrix< T > &shell_matrix, NumericVector< T > &solution_in, NumericVector< T > &rhs_in, const double tol, const unsigned int m_its) libmesh_override
 
virtual std::pair< unsigned int, Realsolve (const ShellMatrix< T > &shell_matrix, const SparseMatrix< T > &precond_matrix, NumericVector< T > &solution_in, NumericVector< T > &rhs_in, const double tol, const unsigned int m_its) libmesh_override
 
PC pc ()
 
KSP ksp ()
 
void get_residual_history (std::vector< double > &hist)
 
Real get_initial_residual ()
 
virtual LinearConvergenceReason get_converged_reason () const libmesh_override
 
bool initialized () const
 
SolverType solver_type () const
 
void set_solver_type (const SolverType st)
 
PreconditionerType preconditioner_type () const
 
void set_preconditioner_type (const PreconditionerType pct)
 
void attach_preconditioner (Preconditioner< T > *preconditioner)
 
virtual void reuse_preconditioner (bool)
 
bool get_same_preconditioner ()
 
std::pair< unsigned int, Realsolve (SparseMatrix< T > &matrix, SparseMatrix< T > *precond_matrix, NumericVector< T > &, NumericVector< T > &, const double, const unsigned int)
 
std::pair< unsigned int, Realsolve (const ShellMatrix< T > &matrix, const SparseMatrix< T > *precond_matrix, NumericVector< T > &, NumericVector< T > &, const double, const unsigned int)
 
virtual void print_converged_reason () const
 
void set_solver_configuration (SolverConfiguration &solver_configuration)
 
const Parallel::Communicatorcomm () const
 
processor_id_type n_processors () const
 
processor_id_type processor_id () const
 

Static Public Member Functions

static std::unique_ptr< LinearSolver< T > > build (const libMesh::Parallel::Communicator &comm_in, const SolverPackage solver_package=libMesh::default_solver_package())
 
static std::string get_info ()
 
static void print_info (std::ostream &out=libMesh::out)
 
static unsigned int n_objects ()
 
static void enable_print_counter_info ()
 
static void disable_print_counter_info ()
 

Protected Types

typedef std::map< std::string, std::pair< unsigned int, unsigned int > > Counts
 

Protected Member Functions

void increment_constructor_count (const std::string &name)
 
void increment_destructor_count (const std::string &name)
 

Protected Attributes

SolverType _solver_type
 
PreconditionerType _preconditioner_type
 
bool _is_initialized
 
Preconditioner< T > * _preconditioner
 
bool same_preconditioner
 
SolverConfiguration_solver_configuration
 
const Parallel::Communicator_communicator
 

Static Protected Attributes

static Counts _counts
 
static Threads::atomic< unsigned int > _n_objects
 
static Threads::spin_mutex _mutex
 
static bool _enable_print_counter = true
 

Private Member Functions

void set_petsc_solver_type ()
 
PetscInt _restrict_solve_to_is_local_size () const
 
void _create_complement_is (const NumericVector< T > &vec_in)
 

Static Private Member Functions

static PetscErrorCode _petsc_shell_matrix_mult (Mat mat, Vec arg, Vec dest)
 
static PetscErrorCode _petsc_shell_matrix_mult_add (Mat mat, Vec arg, Vec add, Vec dest)
 
static PetscErrorCode _petsc_shell_matrix_get_diagonal (Mat mat, Vec dest)
 

Private Attributes

PC _pc
 
KSP _ksp
 
IS _restrict_solve_to_is
 
IS _restrict_solve_to_is_complement
 
SubsetSolveMode _subset_solve_mode
 

Detailed Description

template<typename T>
class libMesh::PetscLinearSolver< T >

This class provides an interface to PETSc iterative solvers that is compatible with the libMesh LinearSolver<>

Author
Benjamin Kirk
Date
2002-2007

Definition at line 82 of file petsc_linear_solver.h.

Member Typedef Documentation

typedef std::map<std::string, std::pair<unsigned int, unsigned int> > libMesh::ReferenceCounter::Counts
protectedinherited

Data structure to log the information. The log is identified by the class name.

Definition at line 119 of file reference_counter.h.

Constructor & Destructor Documentation

template<typename T >
libMesh::PetscLinearSolver< T >::PetscLinearSolver ( const libMesh::Parallel::Communicator &comm_in  LIBMESH_CAN_DEFAULT_TO_COMMWORLD)
inline

Constructor. Initializes Petsc data structures

Definition at line 313 of file petsc_linear_solver.h.

References libMesh::LinearSolver< T >::_preconditioner_type, libMesh::BLOCK_JACOBI_PRECOND, libMesh::ILU_PRECOND, and libMesh::ParallelObject::n_processors().

313  :
314  LinearSolver<T>(comm_in),
318 {
319  if (this->n_processors() == 1)
321  else
323 }
Set dofs outside the subset to zero.
processor_id_type n_processors() const
const class libmesh_nullptr_t libmesh_nullptr
PreconditionerType _preconditioner_type
template<typename T >
libMesh::PetscLinearSolver< T >::~PetscLinearSolver ( )
inline

Destructor.

Definition at line 329 of file petsc_linear_solver.h.

References libMesh::PetscLinearSolver< T >::clear().

330 {
331  this->clear ();
332 }
virtual void clear() libmesh_override

Member Function Documentation

template<typename T >
void libMesh::PetscLinearSolver< T >::_create_complement_is ( const NumericVector< T > &  vec_in)
private

Creates _restrict_solve_to_is_complement to contain all indices that are local in vec_in, except those that are contained in _restrict_solve_to_is.

Definition at line 353 of file petsc_linear_solver.h.

References libMesh::PetscLinearSolver< T >::_restrict_solve_to_is, libMesh::PetscLinearSolver< T >::_restrict_solve_to_is_complement, ierr, and libmesh_nullptr.

360 {
361  libmesh_assert(_restrict_solve_to_is);
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 }
const class libmesh_nullptr_t libmesh_nullptr
PetscErrorCode ierr
template<typename T >
PetscErrorCode libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_get_diagonal ( Mat  mat,
Vec  dest 
)
staticprivate

Internal function if shell matrix mode is used.

Definition at line 1731 of file petsc_linear_solver.C.

References libMesh::ParallelObject::comm(), libMesh::Parallel::Communicator::get(), libMesh::ShellMatrix< T >::get_diagonal(), and ierr.

Referenced by libMesh::PetscLinearSolver< T >::ksp().

1732 {
1733  /* Get the matrix context. */
1734  PetscErrorCode ierr=0;
1735  void * ctx;
1736  ierr = MatShellGetContext(mat,&ctx);
1737 
1738  /* Get user shell matrix object. */
1739  const ShellMatrix<T> & shell_matrix = *static_cast<const ShellMatrix<T> *>(ctx);
1740  CHKERRABORT(shell_matrix.comm().get(), ierr);
1741 
1742  /* Make \p NumericVector instances around the vector. */
1743  PetscVector<T> dest_global(dest, shell_matrix.comm());
1744 
1745  /* Call the user function. */
1746  shell_matrix.get_diagonal(dest_global);
1747 
1748  return ierr;
1749 }
PetscErrorCode ierr
template<typename T >
PetscErrorCode libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_mult ( Mat  mat,
Vec  arg,
Vec  dest 
)
staticprivate

Internal function if shell matrix mode is used.

Definition at line 1677 of file petsc_linear_solver.C.

References libMesh::ParallelObject::comm(), libMesh::Parallel::Communicator::get(), ierr, and libMesh::ShellMatrix< T >::vector_mult().

Referenced by libMesh::PetscLinearSolver< T >::ksp().

1678 {
1679  /* Get the matrix context. */
1680  PetscErrorCode ierr=0;
1681  void * ctx;
1682  ierr = MatShellGetContext(mat,&ctx);
1683 
1684  /* Get user shell matrix object. */
1685  const ShellMatrix<T> & shell_matrix = *static_cast<const ShellMatrix<T> *>(ctx);
1686  CHKERRABORT(shell_matrix.comm().get(), ierr);
1687 
1688  /* Make \p NumericVector instances around the vectors. */
1689  PetscVector<T> arg_global(arg, shell_matrix.comm());
1690  PetscVector<T> dest_global(dest, shell_matrix.comm());
1691 
1692  /* Call the user function. */
1693  shell_matrix.vector_mult(dest_global,arg_global);
1694 
1695  return ierr;
1696 }
PetscErrorCode ierr
template<typename T >
PetscErrorCode libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_mult_add ( Mat  mat,
Vec  arg,
Vec  add,
Vec  dest 
)
staticprivate

Internal function if shell matrix mode is used.

Definition at line 1701 of file petsc_linear_solver.C.

References libMesh::ParallelObject::comm(), libMesh::Parallel::Communicator::get(), ierr, and libMesh::ShellMatrix< T >::vector_mult_add().

Referenced by libMesh::PetscLinearSolver< T >::ksp().

1702 {
1703  /* Get the matrix context. */
1704  PetscErrorCode ierr=0;
1705  void * ctx;
1706  ierr = MatShellGetContext(mat,&ctx);
1707 
1708  /* Get user shell matrix object. */
1709  const ShellMatrix<T> & shell_matrix = *static_cast<const ShellMatrix<T> *>(ctx);
1710  CHKERRABORT(shell_matrix.comm().get(), ierr);
1711 
1712  /* Make \p NumericVector instances around the vectors. */
1713  PetscVector<T> arg_global(arg, shell_matrix.comm());
1714  PetscVector<T> dest_global(dest, shell_matrix.comm());
1715  PetscVector<T> add_global(add, shell_matrix.comm());
1716 
1717  if (add!=arg)
1718  {
1719  arg_global = add_global;
1720  }
1721 
1722  /* Call the user function. */
1723  shell_matrix.vector_mult_add(dest_global,arg_global);
1724 
1725  return ierr;
1726 }
PetscErrorCode ierr
template<typename T >
PetscInt libMesh::PetscLinearSolver< T >::_restrict_solve_to_is_local_size ( ) const
inlineprivate
Returns
The local size of _restrict_solve_to_is.

Definition at line 338 of file petsc_linear_solver.h.

References libMesh::PetscLinearSolver< T >::_restrict_solve_to_is, and ierr.

339 {
340  libmesh_assert(_restrict_solve_to_is);
341 
342  PetscInt s;
343  int ierr = ISGetLocalSize(_restrict_solve_to_is, &s);
344  LIBMESH_CHKERR(ierr);
345 
346  return s;
347 }
PetscErrorCode ierr
template<typename T >
std::pair< unsigned int, Real > libMesh::PetscLinearSolver< T >::adjoint_solve ( SparseMatrix< T > &  matrix_in,
NumericVector< T > &  solution_in,
NumericVector< T > &  rhs_in,
const double  tol,
const unsigned int  m_its 
)
virtual

Call the Petsc solver. It calls the method below, using the same matrix for the system and preconditioner matrices.

Reimplemented from libMesh::LinearSolver< T >.

Definition at line 658 of file petsc_linear_solver.C.

References libMesh::PetscMatrix< T >::close(), ierr, libMesh::TriangleWrapper::init(), libmesh_nullptr, libMesh::NumericVector< T >::local_size(), libMesh::PetscMatrix< T >::mat(), libMesh::SUBSET_COPY_RHS, libMesh::SUBSET_DONT_TOUCH, and libMesh::SUBSET_ZERO.

Referenced by libMesh::PetscLinearSolver< T >::solve().

663 {
664  LOG_SCOPE("solve()", "PetscLinearSolver");
665 
666  // Make sure the data passed in are really of Petsc types
667  PetscMatrix<T> * matrix = cast_ptr<PetscMatrix<T> *>(&matrix_in);
668  // Note that the matrix and precond matrix are the same
669  PetscMatrix<T> * precond = cast_ptr<PetscMatrix<T> *>(&matrix_in);
670  PetscVector<T> * solution = cast_ptr<PetscVector<T> *>(&solution_in);
671  PetscVector<T> * rhs = cast_ptr<PetscVector<T> *>(&rhs_in);
672 
673  this->init (matrix);
674 
675  PetscErrorCode ierr=0;
676  PetscInt its=0, max_its = static_cast<PetscInt>(m_its);
677  PetscReal final_resid=0.;
678 
679  // Close the matrices and vectors in case this wasn't already done.
680  matrix->close ();
681  precond->close ();
682  solution->close ();
683  rhs->close ();
684 
685  Mat submat = libmesh_nullptr;
686  Mat subprecond = libmesh_nullptr;
687  Vec subrhs = libmesh_nullptr;
688  Vec subsolution = libmesh_nullptr;
689  VecScatter scatter = libmesh_nullptr;
690  std::unique_ptr<PetscMatrix<Number>> subprecond_matrix;
691 
692  // Set operators. Also restrict rhs and solution vector to
693  // subdomain if necessary.
695  {
696  PetscInt is_local_size = this->_restrict_solve_to_is_local_size();
697 
698  ierr = VecCreate(this->comm().get(),&subrhs);
699  LIBMESH_CHKERR(ierr);
700  ierr = VecSetSizes(subrhs,is_local_size,PETSC_DECIDE);
701  LIBMESH_CHKERR(ierr);
702  ierr = VecSetFromOptions(subrhs);
703  LIBMESH_CHKERR(ierr);
704 
705  ierr = VecCreate(this->comm().get(),&subsolution);
706  LIBMESH_CHKERR(ierr);
707  ierr = VecSetSizes(subsolution,is_local_size,PETSC_DECIDE);
708  LIBMESH_CHKERR(ierr);
709  ierr = VecSetFromOptions(subsolution);
710  LIBMESH_CHKERR(ierr);
711 
712  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is, subrhs, libmesh_nullptr, &scatter);
713  LIBMESH_CHKERR(ierr);
714 
715  ierr = VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
716  LIBMESH_CHKERR(ierr);
717  ierr = VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
718  LIBMESH_CHKERR(ierr);
719 
720  ierr = VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
721  LIBMESH_CHKERR(ierr);
722  ierr = VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
723  LIBMESH_CHKERR(ierr);
724 
725  ierr = LibMeshCreateSubMatrix(matrix->mat(),
728 #if PETSC_VERSION_LESS_THAN(3,1,0)
729  PETSC_DECIDE,
730 #endif
731  MAT_INITIAL_MATRIX,
732  &submat);
733  LIBMESH_CHKERR(ierr);
734 
735  ierr = LibMeshCreateSubMatrix(precond->mat(),
738 #if PETSC_VERSION_LESS_THAN(3,1,0)
739  PETSC_DECIDE,
740 #endif
741  MAT_INITIAL_MATRIX,
742  &subprecond);
743  LIBMESH_CHKERR(ierr);
744 
745  /* Since removing columns of the matrix changes the equation
746  system, we will now change the right hand side to compensate
747  for this. Note that this is not necessary if \p SUBSET_ZERO
748  has been selected. */
750  {
751  _create_complement_is(rhs_in);
752  PetscInt is_complement_local_size =
753  cast_int<PetscInt>(rhs_in.local_size()-is_local_size);
754 
755  Vec subvec1 = libmesh_nullptr;
756  Mat submat1 = libmesh_nullptr;
757  VecScatter scatter1 = libmesh_nullptr;
758 
759  ierr = VecCreate(this->comm().get(),&subvec1);
760  LIBMESH_CHKERR(ierr);
761  ierr = VecSetSizes(subvec1,is_complement_local_size,PETSC_DECIDE);
762  LIBMESH_CHKERR(ierr);
763  ierr = VecSetFromOptions(subvec1);
764  LIBMESH_CHKERR(ierr);
765 
766  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is_complement, subvec1, libmesh_nullptr, &scatter1);
767  LIBMESH_CHKERR(ierr);
768 
769  ierr = VecScatterBegin(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
770  LIBMESH_CHKERR(ierr);
771  ierr = VecScatterEnd(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
772  LIBMESH_CHKERR(ierr);
773 
774  ierr = VecScale(subvec1,-1.0);
775  LIBMESH_CHKERR(ierr);
776 
777  ierr = LibMeshCreateSubMatrix(matrix->mat(),
780 #if PETSC_VERSION_LESS_THAN(3,1,0)
781  PETSC_DECIDE,
782 #endif
783  MAT_INITIAL_MATRIX,
784  &submat1);
785  LIBMESH_CHKERR(ierr);
786 
787  ierr = MatMultAdd(submat1,subvec1,subrhs,subrhs);
788  LIBMESH_CHKERR(ierr);
789 
790  ierr = LibMeshVecScatterDestroy(&scatter1);
791  LIBMESH_CHKERR(ierr);
792  ierr = LibMeshVecDestroy(&subvec1);
793  LIBMESH_CHKERR(ierr);
794  ierr = LibMeshMatDestroy(&submat1);
795  LIBMESH_CHKERR(ierr);
796  }
797 #if PETSC_RELEASE_LESS_THAN(3,5,0)
798  ierr = KSPSetOperators(_ksp, submat, subprecond,
799  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
800 #else
801  ierr = KSPSetOperators(_ksp, submat, subprecond);
802 
803  PetscBool ksp_reuse_preconditioner = this->same_preconditioner ? PETSC_TRUE : PETSC_FALSE;
804  ierr = KSPSetReusePreconditioner(_ksp, ksp_reuse_preconditioner);
805 #endif
806  LIBMESH_CHKERR(ierr);
807 
808  if (this->_preconditioner)
809  {
810  subprecond_matrix.reset(new PetscMatrix<Number>(subprecond, this->comm()));
811  this->_preconditioner->set_matrix(*subprecond_matrix);
812  this->_preconditioner->init();
813  }
814  }
815  else
816  {
817 #if PETSC_RELEASE_LESS_THAN(3,5,0)
818  ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
819  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
820 #else
821  ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat());
822 
823  PetscBool ksp_reuse_preconditioner = this->same_preconditioner ? PETSC_TRUE : PETSC_FALSE;
824  ierr = KSPSetReusePreconditioner(_ksp, ksp_reuse_preconditioner);
825 #endif
826  LIBMESH_CHKERR(ierr);
827 
828  if (this->_preconditioner)
829  {
830  this->_preconditioner->set_matrix(matrix_in);
831  this->_preconditioner->init();
832  }
833  }
834 
835  // Set the tolerances for the iterative solver. Use the user-supplied
836  // tolerance for the relative residual & leave the others at default values.
837  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
838  PETSC_DEFAULT, max_its);
839  LIBMESH_CHKERR(ierr);
840 
841  // Allow command line options to override anything set programmatically.
842  ierr = KSPSetFromOptions(_ksp);
843  LIBMESH_CHKERR(ierr);
844 
845  // Solve the linear system
846  if (_restrict_solve_to_is != libmesh_nullptr)
847  {
848  ierr = KSPSolveTranspose (_ksp, subrhs, subsolution);
849  LIBMESH_CHKERR(ierr);
850  }
851  else
852  {
853  ierr = KSPSolveTranspose (_ksp, rhs->vec(), solution->vec());
854  LIBMESH_CHKERR(ierr);
855  }
856 
857  // Get the number of iterations required for convergence
858  ierr = KSPGetIterationNumber (_ksp, &its);
859  LIBMESH_CHKERR(ierr);
860 
861  // Get the norm of the final residual to return to the user.
862  ierr = KSPGetResidualNorm (_ksp, &final_resid);
863  LIBMESH_CHKERR(ierr);
864 
865  if (_restrict_solve_to_is != libmesh_nullptr)
866  {
867  switch(_subset_solve_mode)
868  {
869  case SUBSET_ZERO:
870  ierr = VecZeroEntries(solution->vec());
871  LIBMESH_CHKERR(ierr);
872  break;
873 
874  case SUBSET_COPY_RHS:
875  ierr = VecCopy(rhs->vec(),solution->vec());
876  LIBMESH_CHKERR(ierr);
877  break;
878 
879  case SUBSET_DONT_TOUCH:
880  /* Nothing to do here. */
881  break;
882 
883  default:
884  libmesh_error_msg("Invalid subset solve mode = " << _subset_solve_mode);
885  }
886  ierr = VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
887  LIBMESH_CHKERR(ierr);
888  ierr = VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
889  LIBMESH_CHKERR(ierr);
890 
891  ierr = LibMeshVecScatterDestroy(&scatter);
892  LIBMESH_CHKERR(ierr);
893 
894  if (this->_preconditioner)
895  {
896  // Before subprecond_matrix gets cleaned up, we should give
897  // the _preconditioner a different matrix.
898  this->_preconditioner->set_matrix(matrix_in);
899  this->_preconditioner->init();
900  }
901 
902  ierr = LibMeshVecDestroy(&subsolution);
903  LIBMESH_CHKERR(ierr);
904  ierr = LibMeshVecDestroy(&subrhs);
905  LIBMESH_CHKERR(ierr);
906  ierr = LibMeshMatDestroy(&submat);
907  LIBMESH_CHKERR(ierr);
908  ierr = LibMeshMatDestroy(&subprecond);
909  LIBMESH_CHKERR(ierr);
910  }
911 
912  // return the # of its. and the final residual norm.
913  return std::make_pair(its, final_resid);
914 }
Leaves dofs outside the subset unchanged. This is fastest, but also most confusing because it abandon...
Preconditioner< T > * _preconditioner
Set dofs outside the subset to zero.
void _create_complement_is(const NumericVector< T > &vec_in)
PetscInt _restrict_solve_to_is_local_size() const
const class libmesh_nullptr_t libmesh_nullptr
virtual void init(const char *name=libmesh_nullptr) libmesh_override
Set dofs outside the subset to the value of the corresponding dofs of the right hand side...
PetscTruth PetscBool
Definition: petsc_macro.h:64
PetscErrorCode ierr
const Parallel::Communicator & comm() const
template<typename T>
void libMesh::LinearSolver< T >::attach_preconditioner ( Preconditioner< T > *  preconditioner)
inherited

Attaches a Preconditioner object to be used

Definition at line 101 of file linear_solver.C.

References libMesh::libMeshPrivateData::_is_initialized, and libMesh::SHELL_PRECOND.

Referenced by libMesh::LinearSolver< Number >::set_solver_type().

102 {
103  if (this->_is_initialized)
104  libmesh_error_msg("Preconditioner must be attached before the solver is initialized!");
105 
107  _preconditioner = preconditioner;
108 }
Preconditioner< T > * _preconditioner
PreconditionerType _preconditioner_type
template<typename T >
std::unique_ptr< LinearSolver< T > > libMesh::LinearSolver< T >::build ( const libMesh::Parallel::Communicator comm_in,
const SolverPackage  solver_package = libMesh::default_solver_package() 
)
staticinherited

Builds a LinearSolver using the linear solver package specified by solver_package

Definition at line 40 of file linear_solver.C.

References libMesh::EIGEN_SOLVERS, libMesh::LASPACK_SOLVERS, libMesh::libmesh_ignore(), libMesh::PETSC_SOLVERS, and libMesh::TRILINOS_SOLVERS.

Referenced by libMesh::ImplicitSystem::get_linear_solver(), and libMesh::TimeSolver::init().

42 {
43  // Avoid unused parameter warnings when no solver packages are enabled.
45 
46  // Build the appropriate solver
47  switch (solver_package)
48  {
49 #ifdef LIBMESH_HAVE_LASPACK
50  case LASPACK_SOLVERS:
51  return libmesh_make_unique<LaspackLinearSolver<T>>(comm);
52 #endif
53 
54 
55 #ifdef LIBMESH_HAVE_PETSC
56  case PETSC_SOLVERS:
57  return libmesh_make_unique<PetscLinearSolver<T>>(comm);
58 #endif
59 
60 
61 #ifdef LIBMESH_TRILINOS_HAVE_AZTECOO
62  case TRILINOS_SOLVERS:
63  return libmesh_make_unique<AztecLinearSolver<T>>(comm);
64 #endif
65 
66 
67 #ifdef LIBMESH_HAVE_EIGEN
68  case EIGEN_SOLVERS:
69  return libmesh_make_unique<EigenSparseLinearSolver<T>>(comm);
70 #endif
71 
72  default:
73  libmesh_error_msg("ERROR: Unrecognized solver package: " << solver_package);
74  }
75 
76  return std::unique_ptr<LinearSolver<T>>();
77 }
EIGEN_SOLVERS
Definition: libmesh.C:261
TRILINOS_SOLVERS
Definition: libmesh.C:259
LASPACK_SOLVERS
Definition: libmesh.C:263
void libmesh_ignore(const T &)
const Parallel::Communicator & comm() const
template<typename T >
void libMesh::PetscLinearSolver< T >::clear ( )
virtual

Release all memory and clear data structures.

Reimplemented from libMesh::LinearSolver< T >.

Definition at line 102 of file petsc_linear_solver.C.

References libMesh::libMeshPrivateData::_is_initialized, libMesh::BLOCK_JACOBI_PRECOND, libMesh::GMRES, ierr, libMesh::ILU_PRECOND, libMesh::initialized(), libmesh_nullptr, and libMesh::n_processors().

Referenced by libMesh::PetscLinearSolver< T >::~PetscLinearSolver().

103 {
104  if (this->initialized())
105  {
106  /* If we were restricted to some subset, this restriction must
107  be removed and the subset index set destroyed. */
109  {
110  PetscErrorCode ierr = LibMeshISDestroy(&_restrict_solve_to_is);
111  LIBMESH_CHKERR(ierr);
113  }
114 
116  {
117  PetscErrorCode ierr = LibMeshISDestroy(&_restrict_solve_to_is_complement);
118  LIBMESH_CHKERR(ierr);
120  }
121 
122  this->_is_initialized = false;
123 
124  PetscErrorCode ierr=0;
125 
126  ierr = LibMeshKSPDestroy(&_ksp);
127  LIBMESH_CHKERR(ierr);
128 
129  // Mimic PETSc default solver and preconditioner
130  this->_solver_type = GMRES;
131 
132  if (!this->_preconditioner)
133  {
134  if (this->n_processors() == 1)
136  else
138  }
139  }
140 }
bool initialized() const
Definition: linear_solver.h:86
Preconditioner< T > * _preconditioner
processor_id_type n_processors() const
const class libmesh_nullptr_t libmesh_nullptr
PetscErrorCode ierr
PreconditionerType _preconditioner_type
const Parallel::Communicator& libMesh::ParallelObject::comm ( ) const
inlineinherited
Returns
A reference to the Parallel::Communicator object used by this mesh.

Definition at line 87 of file parallel_object.h.

References libMesh::ParallelObject::_communicator.

Referenced by libMesh::__libmesh_petsc_diff_solver_jacobian(), libMesh::__libmesh_petsc_diff_solver_monitor(), libMesh::__libmesh_petsc_diff_solver_residual(), libMesh::__libmesh_petsc_snes_fd_residual(), libMesh::__libmesh_petsc_snes_jacobian(), libMesh::__libmesh_petsc_snes_mffd_residual(), libMesh::__libmesh_petsc_snes_postcheck(), libMesh::__libmesh_petsc_snes_residual(), libMesh::__libmesh_tao_equality_constraints(), libMesh::__libmesh_tao_equality_constraints_jacobian(), libMesh::__libmesh_tao_gradient(), libMesh::__libmesh_tao_hessian(), libMesh::__libmesh_tao_inequality_constraints(), libMesh::__libmesh_tao_inequality_constraints_jacobian(), libMesh::__libmesh_tao_objective(), libMesh::MeshRefinement::_coarsen_elements(), libMesh::ExactSolution::_compute_error(), libMesh::ParmetisPartitioner::_do_repartition(), libMesh::UniformRefinementEstimator::_estimate_error(), libMesh::BoundaryInfo::_find_id_maps(), libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_get_diagonal(), libMesh::SlepcEigenSolver< T >::_petsc_shell_matrix_get_diagonal(), libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_mult(), libMesh::SlepcEigenSolver< T >::_petsc_shell_matrix_mult(), libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_mult_add(), libMesh::EquationSystems::_read_impl(), libMesh::MeshRefinement::_refine_elements(), libMesh::MeshRefinement::_smooth_flags(), libMesh::ImplicitSystem::add_matrix(), libMesh::System::add_vector(), libMesh::EigenSparseLinearSolver< T >::adjoint_solve(), libMesh::UnstructuredMesh::all_second_order(), libMesh::MeshTools::Modification::all_tri(), libMesh::LaplaceMeshSmoother::allgather_graph(), libMesh::FEMSystem::assemble_qoi(), libMesh::MeshCommunication::assign_global_indices(), libMesh::ParmetisPartitioner::assign_partitioning(), libMesh::DofMap::attach_matrix(), libMesh::Parallel::BinSorter< KeyType, IdxType >::binsort(), libMesh::Parallel::Sort< KeyType, IdxType >::binsort(), libMesh::MeshCommunication::broadcast(), libMesh::SparseMatrix< T >::build(), libMesh::MeshTools::Generation::build_extrusion(), libMesh::Parallel::Histogram< KeyType, IdxType >::build_histogram(), libMesh::PetscNonlinearSolver< T >::build_mat_null_space(), libMesh::BoundaryInfo::build_node_list_from_side_list(), libMesh::EquationSystems::build_parallel_solution_vector(), libMesh::MeshBase::cache_elem_dims(), libMesh::System::calculate_norm(), libMesh::DofMap::check_dirichlet_bcid_consistency(), libMesh::DistributedVector< T >::clone(), libMesh::EigenSparseVector< T >::clone(), libMesh::LaspackVector< T >::clone(), libMesh::EpetraVector< T >::clone(), libMesh::PetscVector< T >::clone(), libMesh::EpetraVector< T >::close(), libMesh::Parallel::Sort< KeyType, IdxType >::communicate_bins(), libMesh::Nemesis_IO_Helper::compute_num_global_elem_blocks(), libMesh::Nemesis_IO_Helper::compute_num_global_nodesets(), libMesh::Nemesis_IO_Helper::compute_num_global_sidesets(), libMesh::Problem_Interface::computeF(), libMesh::Problem_Interface::computeJacobian(), libMesh::Problem_Interface::computePreconditioner(), libMesh::ExodusII_IO::copy_elemental_solution(), libMesh::MeshTools::correct_node_proc_ids(), libMesh::MeshTools::create_bounding_box(), libMesh::MeshTools::create_nodal_bounding_box(), libMesh::MeshRefinement::create_parent_error_vector(), libMesh::MeshTools::create_processor_bounding_box(), libMesh::MeshTools::create_subdomain_bounding_box(), libMesh::MeshCommunication::delete_remote_elements(), libMesh::DofMap::distribute_dofs(), DMlibMeshFunction(), DMlibMeshJacobian(), DMlibMeshSetSystem_libMesh(), DMVariableBounds_libMesh(), libMesh::MeshRefinement::eliminate_unrefined_patches(), libMesh::EpetraVector< T >::EpetraVector(), libMesh::WeightedPatchRecoveryErrorEstimator::estimate_error(), libMesh::PatchRecoveryErrorEstimator::estimate_error(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::ExactErrorEstimator::estimate_error(), libMesh::MeshRefinement::flag_elements_by_elem_fraction(), libMesh::MeshRefinement::flag_elements_by_error_fraction(), libMesh::MeshRefinement::flag_elements_by_nelem_target(), libMesh::MeshCommunication::gather(), libMesh::MeshCommunication::gather_neighboring_elements(), libMesh::CondensedEigenSystem::get_eigenpair(), libMesh::DofMap::get_info(), libMesh::ImplicitSystem::get_linear_solver(), libMesh::EquationSystems::get_solution(), libMesh::LocationMap< T >::init(), libMesh::TopologyMap::init(), libMesh::TimeSolver::init(), libMesh::TaoOptimizationSolver< T >::init(), libMesh::PetscNonlinearSolver< T >::init(), libMesh::DistributedVector< T >::init(), libMesh::EpetraVector< T >::init(), libMesh::PetscVector< T >::init(), libMesh::SystemSubsetBySubdomain::init(), libMesh::EigenSystem::init_data(), libMesh::EigenSystem::init_matrices(), libMesh::ParmetisPartitioner::initialize(), libMesh::OptimizationSystem::initialize_equality_constraints_storage(), libMesh::OptimizationSystem::initialize_inequality_constraints_storage(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Elem >(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Node >(), libMesh::MeshTools::libmesh_assert_topology_consistent_procids< Node >(), libMesh::MeshTools::libmesh_assert_valid_boundary_ids(), libMesh::MeshTools::libmesh_assert_valid_dof_ids(), libMesh::MeshTools::libmesh_assert_valid_neighbors(), libMesh::DistributedMesh::libmesh_assert_valid_parallel_flags(), libMesh::DistributedMesh::libmesh_assert_valid_parallel_object_ids(), libMesh::DistributedMesh::libmesh_assert_valid_parallel_p_levels(), libMesh::MeshTools::libmesh_assert_valid_refinement_flags(), libMesh::MeshTools::libmesh_assert_valid_unique_ids(), libMesh::libmesh_petsc_snes_residual_helper(), libMesh::MeshRefinement::limit_level_mismatch_at_edge(), libMesh::MeshRefinement::limit_level_mismatch_at_node(), libMesh::MeshRefinement::limit_overrefined_boundary(), libMesh::MeshRefinement::limit_underrefined_boundary(), libMesh::MeshRefinement::make_coarsening_compatible(), libMesh::MeshCommunication::make_elems_parallel_consistent(), libMesh::MeshRefinement::make_flags_parallel_consistent(), libMesh::MeshCommunication::make_new_node_proc_ids_parallel_consistent(), libMesh::MeshCommunication::make_new_nodes_parallel_consistent(), libMesh::MeshCommunication::make_node_ids_parallel_consistent(), libMesh::MeshCommunication::make_node_proc_ids_parallel_consistent(), libMesh::MeshCommunication::make_node_unique_ids_parallel_consistent(), libMesh::MeshCommunication::make_nodes_parallel_consistent(), libMesh::MeshCommunication::make_p_levels_parallel_consistent(), libMesh::MeshRefinement::make_refinement_compatible(), libMesh::DistributedVector< T >::max(), libMesh::FEMSystem::mesh_position_set(), libMesh::MeshSerializer::MeshSerializer(), libMesh::DistributedVector< T >::min(), libMesh::DistributedMesh::n_active_elem(), libMesh::MeshTools::n_active_levels(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::BoundaryInfo::n_edge_conds(), libMesh::CondensedEigenSystem::n_global_non_condensed_dofs(), libMesh::MeshTools::n_levels(), libMesh::BoundaryInfo::n_nodeset_conds(), libMesh::MeshTools::n_p_levels(), libMesh::BoundaryInfo::n_shellface_conds(), libMesh::DistributedMesh::parallel_max_elem_id(), libMesh::DistributedMesh::parallel_max_node_id(), libMesh::ReplicatedMesh::parallel_max_unique_id(), libMesh::DistributedMesh::parallel_max_unique_id(), libMesh::DistributedMesh::parallel_n_elem(), libMesh::DistributedMesh::parallel_n_nodes(), libMesh::SparsityPattern::Build::parallel_sync(), libMesh::MeshTools::paranoid_n_levels(), libMesh::Partitioner::partition(), libMesh::LinearPartitioner::partition_range(), libMesh::MetisPartitioner::partition_range(), libMesh::Partitioner::partition_unpartitioned_elements(), libMesh::petsc_auto_fieldsplit(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::MeshBase::prepare_for_use(), libMesh::SparseMatrix< T >::print(), libMesh::Nemesis_IO::read(), libMesh::XdrIO::read(), libMesh::CheckpointIO::read_header(), libMesh::XdrIO::read_header(), libMesh::XdrIO::read_serialized_bc_names(), libMesh::XdrIO::read_serialized_bcs_helper(), libMesh::XdrIO::read_serialized_connectivity(), libMesh::XdrIO::read_serialized_nodes(), libMesh::XdrIO::read_serialized_nodesets(), libMesh::XdrIO::read_serialized_subdomain_names(), libMesh::MeshBase::recalculate_n_partitions(), libMesh::MeshCommunication::redistribute(), libMesh::MeshRefinement::refine_and_coarsen_elements(), libMesh::DistributedMesh::renumber_dof_objects(), libMesh::CheckpointIO::select_split_config(), libMesh::MeshCommunication::send_coarse_ghosts(), libMesh::Partitioner::set_node_processor_ids(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::Partitioner::set_parent_processor_ids(), libMesh::PetscDiffSolver::setup_petsc_data(), libMesh::LaplaceMeshSmoother::smooth(), libMesh::Parallel::Sort< KeyType, IdxType >::sort(), libMesh::split_mesh(), libMesh::MeshBase::subdomain_ids(), libMesh::BoundaryInfo::sync(), libMesh::Parallel::sync_element_data_by_parent_id(), libMesh::Parallel::sync_node_data_by_element_id(), libMesh::MeshRefinement::test_level_one(), libMesh::MeshRefinement::test_unflagged(), libMesh::MeshTools::total_weight(), libMesh::MeshRefinement::uniformly_coarsen(), libMesh::NameBasedIO::write(), libMesh::XdrIO::write(), libMesh::XdrIO::write_serialized_bcs_helper(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::XdrIO::write_serialized_nodes(), libMesh::XdrIO::write_serialized_nodesets(), libMesh::DistributedVector< T >::zero_clone(), libMesh::EigenSparseVector< T >::zero_clone(), libMesh::LaspackVector< T >::zero_clone(), libMesh::EpetraVector< T >::zero_clone(), and libMesh::PetscVector< T >::zero_clone().

88  { return _communicator; }
const Parallel::Communicator & _communicator
void libMesh::ReferenceCounter::disable_print_counter_info ( )
staticinherited
void libMesh::ReferenceCounter::enable_print_counter_info ( )
staticinherited

Methods to enable/disable the reference counter output from print_info()

Definition at line 101 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

Referenced by libMesh::ReferenceCounter::n_objects().

102 {
103  _enable_print_counter = true;
104  return;
105 }
template<typename T >
LinearConvergenceReason libMesh::PetscLinearSolver< T >::get_converged_reason ( ) const
virtual
Returns
The solver's convergence flag

Implements libMesh::LinearSolver< T >.

Definition at line 1634 of file petsc_linear_solver.C.

References libMesh::CONVERGED_ATOL, libMesh::CONVERGED_ATOL_NORMAL, libMesh::CONVERGED_CG_CONSTRAINED, libMesh::CONVERGED_CG_NEG_CURVE, libMesh::CONVERGED_HAPPY_BREAKDOWN, libMesh::CONVERGED_ITERATING, libMesh::CONVERGED_ITS, libMesh::CONVERGED_RTOL, libMesh::CONVERGED_RTOL_NORMAL, libMesh::CONVERGED_STEP_LENGTH, libMesh::DIVERGED_BREAKDOWN, libMesh::DIVERGED_BREAKDOWN_BICG, libMesh::DIVERGED_DTOL, libMesh::DIVERGED_INDEFINITE_MAT, libMesh::DIVERGED_INDEFINITE_PC, libMesh::DIVERGED_ITS, libMesh::DIVERGED_NAN, libMesh::DIVERGED_NONSYMMETRIC, libMesh::DIVERGED_NULL, libMesh::DIVERGED_PCSETUP_FAILED, libMesh::err, and libMesh::UNKNOWN_FLAG.

Referenced by libMesh::PetscLinearSolver< T >::ksp().

1635 {
1636  KSPConvergedReason reason;
1637  KSPGetConvergedReason(_ksp, &reason);
1638 
1639  switch(reason)
1640  {
1641 #if !PETSC_VERSION_LESS_THAN(3,2,0)
1642  case KSP_CONVERGED_RTOL_NORMAL : return CONVERGED_RTOL_NORMAL;
1643  case KSP_CONVERGED_ATOL_NORMAL : return CONVERGED_ATOL_NORMAL;
1644 #endif
1645  case KSP_CONVERGED_RTOL : return CONVERGED_RTOL;
1646  case KSP_CONVERGED_ATOL : return CONVERGED_ATOL;
1647  case KSP_CONVERGED_ITS : return CONVERGED_ITS;
1648  case KSP_CONVERGED_CG_NEG_CURVE : return CONVERGED_CG_NEG_CURVE;
1649  case KSP_CONVERGED_CG_CONSTRAINED : return CONVERGED_CG_CONSTRAINED;
1650  case KSP_CONVERGED_STEP_LENGTH : return CONVERGED_STEP_LENGTH;
1651  case KSP_CONVERGED_HAPPY_BREAKDOWN : return CONVERGED_HAPPY_BREAKDOWN;
1652  case KSP_DIVERGED_NULL : return DIVERGED_NULL;
1653  case KSP_DIVERGED_ITS : return DIVERGED_ITS;
1654  case KSP_DIVERGED_DTOL : return DIVERGED_DTOL;
1655  case KSP_DIVERGED_BREAKDOWN : return DIVERGED_BREAKDOWN;
1656  case KSP_DIVERGED_BREAKDOWN_BICG : return DIVERGED_BREAKDOWN_BICG;
1657  case KSP_DIVERGED_NONSYMMETRIC : return DIVERGED_NONSYMMETRIC;
1658  case KSP_DIVERGED_INDEFINITE_PC : return DIVERGED_INDEFINITE_PC;
1659 #if PETSC_VERSION_LESS_THAN(3,4,0)
1660  case KSP_DIVERGED_NAN : return DIVERGED_NAN;
1661 #else
1662  case KSP_DIVERGED_NANORINF : return DIVERGED_NAN;
1663 #endif
1664  case KSP_DIVERGED_INDEFINITE_MAT : return DIVERGED_INDEFINITE_MAT;
1665  case KSP_CONVERGED_ITERATING : return CONVERGED_ITERATING;
1666 #if !PETSC_VERSION_LESS_THAN(3,7,0)
1667  case KSP_DIVERGED_PCSETUP_FAILED : return DIVERGED_PCSETUP_FAILED;
1668 #endif
1669  default :
1670  libMesh::err << "Unknown convergence flag!" << std::endl;
1671  return UNKNOWN_FLAG;
1672  }
1673 }
OStreamProxy err(std::cerr)
std::string libMesh::ReferenceCounter::get_info ( )
staticinherited

Gets a string containing the reference information.

Definition at line 47 of file reference_counter.C.

References libMesh::ReferenceCounter::_counts, and libMesh::Quality::name().

Referenced by libMesh::ReferenceCounter::print_info().

48 {
49 #if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
50 
51  std::ostringstream oss;
52 
53  oss << '\n'
54  << " ---------------------------------------------------------------------------- \n"
55  << "| Reference count information |\n"
56  << " ---------------------------------------------------------------------------- \n";
57 
58  for (Counts::iterator it = _counts.begin();
59  it != _counts.end(); ++it)
60  {
61  const std::string name(it->first);
62  const unsigned int creations = it->second.first;
63  const unsigned int destructions = it->second.second;
64 
65  oss << "| " << name << " reference count information:\n"
66  << "| Creations: " << creations << '\n'
67  << "| Destructions: " << destructions << '\n';
68  }
69 
70  oss << " ---------------------------------------------------------------------------- \n";
71 
72  return oss.str();
73 
74 #else
75 
76  return "";
77 
78 #endif
79 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:39
template<typename T >
Real libMesh::PetscLinearSolver< T >::get_initial_residual ( )
Returns
Just the initial residual for the solve just completed with this interface.

Use this method instead of the one above if you just want the starting residual and not the entire history.

Definition at line 1520 of file petsc_linear_solver.C.

References libMesh::err, and ierr.

Referenced by libMesh::PetscLinearSolver< T >::ksp().

1521 {
1522  PetscErrorCode ierr = 0;
1523  PetscInt its = 0;
1524 
1525  // Fill the residual history vector with the residual norms
1526  // Note that GetResidualHistory() does not copy any values, it
1527  // simply sets the pointer p. Note that for some Krylov subspace
1528  // methods, the number of residuals returned in the history
1529  // vector may be different from what you are expecting. For
1530  // example, TFQMR returns two residual values per iteration step.
1531  PetscReal * p;
1532  ierr = KSPGetResidualHistory(_ksp, &p, &its);
1533  LIBMESH_CHKERR(ierr);
1534 
1535  // Check no residual history
1536  if (its == 0)
1537  {
1538  libMesh::err << "No iterations have been performed, returning 0." << std::endl;
1539  return 0.;
1540  }
1541 
1542  // Otherwise, return the value pointed to by p.
1543  return *p;
1544 }
OStreamProxy err(std::cerr)
PetscErrorCode ierr
template<typename T >
void libMesh::PetscLinearSolver< T >::get_residual_history ( std::vector< double > &  hist)

Fills the input vector with the sequence of residual norms from the latest iterative solve.

Definition at line 1487 of file petsc_linear_solver.C.

References ierr.

Referenced by libMesh::PetscLinearSolver< T >::ksp().

1488 {
1489  PetscErrorCode ierr = 0;
1490  PetscInt its = 0;
1491 
1492  // Fill the residual history vector with the residual norms
1493  // Note that GetResidualHistory() does not copy any values, it
1494  // simply sets the pointer p. Note that for some Krylov subspace
1495  // methods, the number of residuals returned in the history
1496  // vector may be different from what you are expecting. For
1497  // example, TFQMR returns two residual values per iteration step.
1498  PetscReal * p;
1499  ierr = KSPGetResidualHistory(_ksp, &p, &its);
1500  LIBMESH_CHKERR(ierr);
1501 
1502  // Check for early return
1503  if (its == 0) return;
1504 
1505  // Create space to store the result
1506  hist.resize(its);
1507 
1508  // Copy history into the vector provided by the user.
1509  for (PetscInt i=0; i<its; ++i)
1510  {
1511  hist[i] = *p;
1512  p++;
1513  }
1514 }
PetscErrorCode ierr
template<typename T >
bool libMesh::LinearSolver< T >::get_same_preconditioner ( )
inlineinherited
Returns
same_preconditioner, which indicates if we reuse the same preconditioner for subsequent solves.

Definition at line 327 of file linear_solver.h.

References libMesh::LinearSolver< T >::same_preconditioner.

Referenced by libMesh::LinearSolver< Number >::set_solver_type().

328 {
329  return same_preconditioner;
330 }
void libMesh::ReferenceCounter::increment_constructor_count ( const std::string &  name)
inlineprotectedinherited

Increments the construction counter. Should be called in the constructor of any derived class that will be reference counted.

Definition at line 185 of file reference_counter.h.

References libMesh::ReferenceCounter::_counts, libMesh::Quality::name(), and libMesh::Threads::spin_mtx.

Referenced by libMesh::ReferenceCounter::n_objects(), and libMesh::ReferenceCountedObject< RBParametrized >::ReferenceCountedObject().

186 {
187  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
188  std::pair<unsigned int, unsigned int> & p = _counts[name];
189 
190  p.first++;
191 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:39
spin_mutex spin_mtx
Definition: threads.C:29
void libMesh::ReferenceCounter::increment_destructor_count ( const std::string &  name)
inlineprotectedinherited

Increments the destruction counter. Should be called in the destructor of any derived class that will be reference counted.

Definition at line 198 of file reference_counter.h.

References libMesh::ReferenceCounter::_counts, libMesh::Quality::name(), and libMesh::Threads::spin_mtx.

Referenced by libMesh::ReferenceCounter::n_objects(), and libMesh::ReferenceCountedObject< RBParametrized >::~ReferenceCountedObject().

199 {
200  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
201  std::pair<unsigned int, unsigned int> & p = _counts[name];
202 
203  p.second++;
204 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:39
spin_mutex spin_mtx
Definition: threads.C:29
template<typename T >
void libMesh::PetscLinearSolver< T >::init ( const char *  name = libmesh_nullptr)
virtual

Initialize data structures if not done so already. Assigns a name, which is turned into an underscore-separated prefix for the underlying KSP object.

Implements libMesh::LinearSolver< T >.

Definition at line 145 of file petsc_linear_solver.C.

References libMesh::__libmesh_petsc_preconditioner_apply(), libMesh::__libmesh_petsc_preconditioner_setup(), libMesh::libMeshPrivateData::_is_initialized, ierr, libMesh::initialized(), and libMesh::PetscPreconditioner< T >::set_petsc_preconditioner_type().

Referenced by libMesh::PetscLinearSolver< T >::ksp(), and libMesh::PetscLinearSolver< T >::pc().

146 {
147  // Initialize the data structures if not done so already.
148  if (!this->initialized())
149  {
150  this->_is_initialized = true;
151 
152  PetscErrorCode ierr=0;
153 
154  // Create the linear solver context
155  ierr = KSPCreate (this->comm().get(), &_ksp);
156  LIBMESH_CHKERR(ierr);
157 
158  if (name)
159  {
160  ierr = KSPSetOptionsPrefix(_ksp, name);
161  LIBMESH_CHKERR(ierr);
162  }
163 
164  // Create the preconditioner context
165  ierr = KSPGetPC (_ksp, &_pc);
166  LIBMESH_CHKERR(ierr);
167 
168  // Set user-specified solver and preconditioner types
169  this->set_petsc_solver_type();
170 
171  // If the SolverConfiguration object is provided, use it to set
172  // options during solver initialization.
173  if (this->_solver_configuration)
174  {
176  }
177 
178  // Set the options from user-input
179  // Set runtime options, e.g.,
180  // -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
181  // These options will override those specified above as long as
182  // KSPSetFromOptions() is called _after_ any other customization
183  // routines.
184 
185  ierr = KSPSetFromOptions (_ksp);
186  LIBMESH_CHKERR(ierr);
187 
188  // Not sure if this is necessary, or if it is already handled by KSPSetFromOptions?
189  // NOT NECESSARY!!!!
190  //ierr = PCSetFromOptions (_pc);
191  //LIBMESH_CHKERR(ierr);
192 
193  // Have the Krylov subspace method use our good initial guess
194  // rather than 0, unless the user requested a KSPType of
195  // preonly, which complains if asked to use initial guesses.
196 #if PETSC_VERSION_LESS_THAN(3,0,0) || !PETSC_RELEASE_LESS_THAN(3,4,0)
197  // Pre-3.0 and petsc-dev (as of October 2012) use non-const versions
198  KSPType ksp_type;
199 #else
200  const KSPType ksp_type;
201 #endif
202 
203  ierr = KSPGetType (_ksp, &ksp_type);
204  LIBMESH_CHKERR(ierr);
205 
206  if (strcmp(ksp_type, "preonly"))
207  {
208  ierr = KSPSetInitialGuessNonzero (_ksp, PETSC_TRUE);
209  LIBMESH_CHKERR(ierr);
210  }
211 
212  // Notify PETSc of location to store residual history.
213  // This needs to be called before any solves, since
214  // it sets the residual history length to zero. The default
215  // behavior is for PETSc to allocate (internally) an array
216  // of size 1000 to hold the residual norm history.
217  ierr = KSPSetResidualHistory(_ksp,
218  PETSC_NULL, // pointer to the array which holds the history
219  PETSC_DECIDE, // size of the array holding the history
220  PETSC_TRUE); // Whether or not to reset the history for each solve.
221  LIBMESH_CHKERR(ierr);
222 
224 
225  //If there is a preconditioner object we need to set the internal setup and apply routines
226  if (this->_preconditioner)
227  {
228  this->_preconditioner->init();
229  PCShellSetContext(_pc,(void *)this->_preconditioner);
230  PCShellSetSetUp(_pc,__libmesh_petsc_preconditioner_setup);
231  PCShellSetApply(_pc,__libmesh_petsc_preconditioner_apply);
232  }
233  }
234 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:39
bool initialized() const
Definition: linear_solver.h:86
Preconditioner< T > * _preconditioner
SolverConfiguration * _solver_configuration
PetscErrorCode __libmesh_petsc_preconditioner_setup(void *ctx)
static void set_petsc_preconditioner_type(const PreconditionerType &preconditioner_type, PC &pc)
PetscErrorCode ierr
const Parallel::Communicator & comm() const
PreconditionerType _preconditioner_type
PetscErrorCode __libmesh_petsc_preconditioner_apply(void *ctx, Vec x, Vec y)
template<typename T >
void libMesh::PetscLinearSolver< T >::init ( PetscMatrix< T > *  matrix,
const char *  name = libmesh_nullptr 
)

Initialize data structures if not done so already plus much more

Definition at line 238 of file petsc_linear_solver.C.

References libMesh::__libmesh_petsc_preconditioner_apply(), libMesh::__libmesh_petsc_preconditioner_setup(), libMesh::libMeshPrivateData::_is_initialized, ierr, libMesh::initialized(), libMesh::PetscMatrix< T >::mat(), and libMesh::PetscPreconditioner< T >::set_petsc_preconditioner_type().

240 {
241  // Initialize the data structures if not done so already.
242  if (!this->initialized())
243  {
244  this->_is_initialized = true;
245 
246  PetscErrorCode ierr=0;
247 
248  // Create the linear solver context
249  ierr = KSPCreate (this->comm().get(), &_ksp);
250  LIBMESH_CHKERR(ierr);
251 
252  if (name)
253  {
254  ierr = KSPSetOptionsPrefix(_ksp, name);
255  LIBMESH_CHKERR(ierr);
256  }
257 
258  //ierr = PCCreate (this->comm().get(), &_pc);
259  // LIBMESH_CHKERR(ierr);
260 
261  // Create the preconditioner context
262  ierr = KSPGetPC (_ksp, &_pc);
263  LIBMESH_CHKERR(ierr);
264 
265  // Set operators. The input matrix works as the preconditioning matrix
266 #if PETSC_RELEASE_LESS_THAN(3,5,0)
267  ierr = KSPSetOperators(_ksp, matrix->mat(), matrix->mat(),DIFFERENT_NONZERO_PATTERN);
268 #else
269  ierr = KSPSetOperators(_ksp, matrix->mat(), matrix->mat());
270 #endif
271  LIBMESH_CHKERR(ierr);
272 
273  // Set user-specified solver and preconditioner types
274  this->set_petsc_solver_type();
275 
276  // If the SolverConfiguration object is provided, use it to set
277  // options during solver initialization.
278  if (this->_solver_configuration)
279  {
281  }
282 
283  // Set the options from user-input
284  // Set runtime options, e.g.,
285  // -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
286  // These options will override those specified above as long as
287  // KSPSetFromOptions() is called _after_ any other customization
288  // routines.
289 
290  ierr = KSPSetFromOptions (_ksp);
291  LIBMESH_CHKERR(ierr);
292 
293  // Not sure if this is necessary, or if it is already handled by KSPSetFromOptions?
294  // NOT NECESSARY!!!!
295  //ierr = PCSetFromOptions (_pc);
296  //LIBMESH_CHKERR(ierr);
297 
298  // Have the Krylov subspace method use our good initial guess
299  // rather than 0, unless the user requested a KSPType of
300  // preonly, which complains if asked to use initial guesses.
301 #if PETSC_VERSION_LESS_THAN(3,0,0) || !PETSC_RELEASE_LESS_THAN(3,4,0)
302  KSPType ksp_type;
303 #else
304  const KSPType ksp_type;
305 #endif
306 
307  ierr = KSPGetType (_ksp, &ksp_type);
308  LIBMESH_CHKERR(ierr);
309 
310  if (strcmp(ksp_type, "preonly"))
311  {
312  ierr = KSPSetInitialGuessNonzero (_ksp, PETSC_TRUE);
313  LIBMESH_CHKERR(ierr);
314  }
315 
316  // Notify PETSc of location to store residual history.
317  // This needs to be called before any solves, since
318  // it sets the residual history length to zero. The default
319  // behavior is for PETSc to allocate (internally) an array
320  // of size 1000 to hold the residual norm history.
321  ierr = KSPSetResidualHistory(_ksp,
322  PETSC_NULL, // pointer to the array which holds the history
323  PETSC_DECIDE, // size of the array holding the history
324  PETSC_TRUE); // Whether or not to reset the history for each solve.
325  LIBMESH_CHKERR(ierr);
326 
328  if (this->_preconditioner)
329  {
330  this->_preconditioner->set_matrix(*matrix);
331  this->_preconditioner->init();
332  PCShellSetContext(_pc,(void *)this->_preconditioner);
333  PCShellSetSetUp(_pc,__libmesh_petsc_preconditioner_setup);
334  PCShellSetApply(_pc,__libmesh_petsc_preconditioner_apply);
335  }
336  }
337 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:39
bool initialized() const
Definition: linear_solver.h:86
Preconditioner< T > * _preconditioner
SolverConfiguration * _solver_configuration
PetscErrorCode __libmesh_petsc_preconditioner_setup(void *ctx)
static void set_petsc_preconditioner_type(const PreconditionerType &preconditioner_type, PC &pc)
PetscErrorCode ierr
const Parallel::Communicator & comm() const
PreconditionerType _preconditioner_type
PetscErrorCode __libmesh_petsc_preconditioner_apply(void *ctx, Vec x, Vec y)
template<typename T >
void libMesh::PetscLinearSolver< T >::init_names ( const System sys)
virtual

Apply names to the system to be solved. This sets an option prefix from the system name and sets field names from the system's variable names.

Since field names are applied to DoF numberings, this method must be called again after any System reinit.

Reimplemented from libMesh::LinearSolver< T >.

Definition at line 343 of file petsc_linear_solver.C.

References libMesh::petsc_auto_fieldsplit().

344 {
345  petsc_auto_fieldsplit(this->pc(), sys);
346 }
void petsc_auto_fieldsplit(PC my_pc, const System &sys)
template<typename T>
bool libMesh::LinearSolver< T >::initialized ( ) const
inlineinherited
Returns
true if the data structures are initialized, false otherwise.

Definition at line 86 of file linear_solver.h.

Referenced by libMesh::EigenSparseLinearSolver< T >::clear(), and libMesh::EigenSparseLinearSolver< T >::init().

86 { return _is_initialized; }
static unsigned int libMesh::ReferenceCounter::n_objects ( )
inlinestaticinherited
processor_id_type libMesh::ParallelObject::n_processors ( ) const
inlineinherited
Returns
The number of processors in the group.

Definition at line 93 of file parallel_object.h.

References libMesh::ParallelObject::_communicator, and libMesh::Parallel::Communicator::size().

Referenced by libMesh::ParmetisPartitioner::_do_repartition(), libMesh::BoundaryInfo::_find_id_maps(), libMesh::DistributedMesh::add_elem(), libMesh::DistributedMesh::add_node(), libMesh::LaplaceMeshSmoother::allgather_graph(), libMesh::FEMSystem::assembly(), libMesh::ParmetisPartitioner::assign_partitioning(), libMesh::AztecLinearSolver< T >::AztecLinearSolver(), libMesh::MeshCommunication::broadcast(), libMesh::BoundaryInfo::build_node_list_from_side_list(), libMesh::DistributedMesh::clear(), libMesh::Nemesis_IO_Helper::compute_border_node_ids(), libMesh::Nemesis_IO_Helper::construct_nemesis_filename(), libMesh::MeshTools::correct_node_proc_ids(), libMesh::UnstructuredMesh::create_pid_mesh(), libMesh::MeshTools::create_processor_bounding_box(), libMesh::DofMap::distribute_dofs(), libMesh::DofMap::distribute_local_dofs_node_major(), libMesh::DofMap::distribute_local_dofs_var_major(), libMesh::DistributedMesh::DistributedMesh(), libMesh::EnsightIO::EnsightIO(), libMesh::MeshCommunication::gather(), libMesh::MeshCommunication::gather_neighboring_elements(), libMesh::MeshBase::get_info(), libMesh::EquationSystems::get_solution(), libMesh::DistributedVector< T >::init(), libMesh::SystemSubsetBySubdomain::init(), libMesh::ParmetisPartitioner::initialize(), libMesh::Nemesis_IO_Helper::initialize(), libMesh::DistributedMesh::insert_elem(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Elem >(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Node >(), libMesh::MeshTools::libmesh_assert_topology_consistent_procids< Node >(), libMesh::MeshTools::libmesh_assert_valid_boundary_ids(), libMesh::MeshTools::libmesh_assert_valid_dof_ids(), libMesh::MeshTools::libmesh_assert_valid_neighbors(), libMesh::MeshTools::libmesh_assert_valid_refinement_flags(), libMesh::DofMap::local_variable_indices(), libMesh::MeshRefinement::make_coarsening_compatible(), libMesh::MeshBase::n_active_elem_on_proc(), libMesh::MeshBase::n_elem_on_proc(), libMesh::MeshBase::n_nodes_on_proc(), libMesh::SparsityPattern::Build::parallel_sync(), libMesh::Partitioner::partition(), libMesh::MeshBase::partition(), libMesh::Partitioner::partition_unpartitioned_elements(), libMesh::PetscLinearSolver< T >::PetscLinearSolver(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::SparseMatrix< T >::print(), libMesh::NameBasedIO::read(), libMesh::Nemesis_IO::read(), libMesh::CheckpointIO::read(), libMesh::CheckpointIO::read_connectivity(), libMesh::XdrIO::read_header(), libMesh::CheckpointIO::read_nodes(), libMesh::MeshCommunication::redistribute(), libMesh::DistributedMesh::renumber_dof_objects(), libMesh::Partitioner::repartition(), libMesh::MeshCommunication::send_coarse_ghosts(), libMesh::Partitioner::set_node_processor_ids(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::Parallel::Sort< KeyType, IdxType >::sort(), libMesh::MeshRefinement::uniformly_coarsen(), libMesh::DistributedMesh::update_parallel_id_counts(), libMesh::GMVIO::write_binary(), libMesh::GMVIO::write_discontinuous_gmv(), libMesh::XdrIO::write_serialized_bcs_helper(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::XdrIO::write_serialized_nodes(), and libMesh::XdrIO::write_serialized_nodesets().

94  { return cast_int<processor_id_type>(_communicator.size()); }
unsigned int size() const
Definition: parallel.h:725
const Parallel::Communicator & _communicator
template<typename T >
PC libMesh::PetscLinearSolver< T >::pc ( )
inline
Returns
The raw PETSc preconditioner context pointer.

This allows you to specify the PCShellSetApply() and PCShellSetSetUp() functions if you desire. Just don't do anything crazy like calling libMeshPCDestroy() on the pointer!

Definition at line 214 of file petsc_linear_solver.h.

References libMesh::PetscLinearSolver< T >::_pc, and libMesh::PetscLinearSolver< T >::init().

214 { this->init(); return _pc; }
virtual void init(const char *name=libmesh_nullptr) libmesh_override
template<typename T >
PreconditionerType libMesh::LinearSolver< T >::preconditioner_type ( ) const
inherited
Returns
The type of preconditioner to use.

Definition at line 81 of file linear_solver.C.

Referenced by libMesh::LinearSolver< Number >::set_solver_type().

82 {
83  if (_preconditioner)
84  return _preconditioner->type();
85 
86  return _preconditioner_type;
87 }
Preconditioner< T > * _preconditioner
PreconditionerType _preconditioner_type
template<typename T >
void libMesh::LinearSolver< T >::print_converged_reason ( ) const
virtualinherited

Prints a useful message about why the latest linear solve con(di)verged.

Reimplemented in libMesh::AztecLinearSolver< T >, and libMesh::LaspackLinearSolver< T >.

Definition at line 153 of file linear_solver.C.

References libMesh::Utility::enum_to_string(), and libMesh::out.

Referenced by libMesh::LinearSolver< Number >::set_solver_type().

154 {
156  libMesh::out << "Linear solver convergence/divergence reason: " << Utility::enum_to_string(reason) << std::endl;
157 }
virtual LinearConvergenceReason get_converged_reason() const =0
std::string enum_to_string(const T e)
OStreamProxy out(std::cout)
void libMesh::ReferenceCounter::print_info ( std::ostream &  out = libMesh::out)
staticinherited

Prints the reference information, by default to libMesh::out.

Definition at line 88 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter, and libMesh::ReferenceCounter::get_info().

Referenced by libMesh::LibMeshInit::LibMeshInit().

89 {
91  out_stream << ReferenceCounter::get_info();
92 }
static std::string get_info()
processor_id_type libMesh::ParallelObject::processor_id ( ) const
inlineinherited
Returns
The rank of this processor in the group.

Definition at line 99 of file parallel_object.h.

References libMesh::ParallelObject::_communicator, and libMesh::Parallel::Communicator::rank().

Referenced by libMesh::BoundaryInfo::_find_id_maps(), libMesh::EquationSystems::_read_impl(), libMesh::DistributedMesh::add_elem(), libMesh::BoundaryInfo::add_elements(), libMesh::DofMap::add_neighbors_to_send_list(), libMesh::DistributedMesh::add_node(), libMesh::MeshRefinement::add_node(), libMesh::MeshTools::Modification::all_tri(), libMesh::FEMSystem::assembly(), libMesh::ParmetisPartitioner::assign_partitioning(), libMesh::MeshCommunication::broadcast(), libMesh::EquationSystems::build_discontinuous_solution_vector(), libMesh::Nemesis_IO_Helper::build_element_and_node_maps(), libMesh::ParmetisPartitioner::build_graph(), libMesh::InfElemBuilder::build_inf_elem(), libMesh::BoundaryInfo::build_node_list_from_side_list(), libMesh::DofMap::build_sparsity(), libMesh::DistributedMesh::clear(), libMesh::ExodusII_IO_Helper::close(), libMesh::Nemesis_IO_Helper::compute_border_node_ids(), libMesh::Nemesis_IO_Helper::compute_communication_map_parameters(), libMesh::Nemesis_IO_Helper::compute_internal_and_border_elems_and_internal_nodes(), libMesh::Nemesis_IO_Helper::compute_node_communication_maps(), libMesh::Nemesis_IO_Helper::compute_num_global_elem_blocks(), libMesh::Nemesis_IO_Helper::compute_num_global_nodesets(), libMesh::Nemesis_IO_Helper::compute_num_global_sidesets(), libMesh::Nemesis_IO_Helper::construct_nemesis_filename(), libMesh::MeshTools::correct_node_proc_ids(), libMesh::ExodusII_IO_Helper::create(), libMesh::DistributedMesh::delete_elem(), libMesh::DistributedMesh::delete_node(), libMesh::MeshCommunication::delete_remote_elements(), libMesh::DofMap::distribute_dofs(), libMesh::DofMap::distribute_local_dofs_node_major(), libMesh::DofMap::distribute_local_dofs_var_major(), libMesh::DistributedMesh::DistributedMesh(), libMesh::EnsightIO::EnsightIO(), libMesh::MeshFunction::find_element(), libMesh::MeshFunction::find_elements(), libMesh::UnstructuredMesh::find_neighbors(), libMesh::MeshCommunication::gather(), libMesh::MeshCommunication::gather_neighboring_elements(), libMesh::Nemesis_IO_Helper::get_cmap_params(), libMesh::Nemesis_IO_Helper::get_eb_info_global(), libMesh::Nemesis_IO_Helper::get_elem_cmap(), libMesh::Nemesis_IO_Helper::get_elem_map(), libMesh::MeshBase::get_info(), libMesh::DofMap::get_info(), libMesh::Nemesis_IO_Helper::get_init_global(), libMesh::Nemesis_IO_Helper::get_init_info(), libMesh::Nemesis_IO_Helper::get_loadbal_param(), libMesh::Nemesis_IO_Helper::get_node_cmap(), libMesh::Nemesis_IO_Helper::get_node_map(), libMesh::Nemesis_IO_Helper::get_ns_param_global(), libMesh::EquationSystems::get_solution(), libMesh::Nemesis_IO_Helper::get_ss_param_global(), libMesh::SparsityPattern::Build::handle_vi_vj(), libMesh::DistributedVector< T >::init(), libMesh::SystemSubsetBySubdomain::init(), libMesh::ParmetisPartitioner::initialize(), libMesh::ExodusII_IO_Helper::initialize(), libMesh::ExodusII_IO_Helper::initialize_element_variables(), libMesh::ExodusII_IO_Helper::initialize_global_variables(), libMesh::ExodusII_IO_Helper::initialize_nodal_variables(), libMesh::DistributedMesh::insert_elem(), libMesh::DofMap::is_evaluable(), libMesh::SparsityPattern::Build::join(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Elem >(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Node >(), libMesh::MeshTools::libmesh_assert_valid_neighbors(), libMesh::DistributedMesh::libmesh_assert_valid_parallel_object_ids(), libMesh::DofMap::local_variable_indices(), libMesh::MeshRefinement::make_coarsening_compatible(), libMesh::MeshBase::n_active_local_elem(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::BoundaryInfo::n_edge_conds(), libMesh::System::n_local_dofs(), libMesh::MeshBase::n_local_elem(), libMesh::MeshBase::n_local_nodes(), libMesh::BoundaryInfo::n_nodeset_conds(), libMesh::BoundaryInfo::n_shellface_conds(), libMesh::WeightedPatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::SparsityPattern::Build::operator()(), libMesh::PatchRecoveryErrorEstimator::EstimateError::operator()(), libMesh::SparsityPattern::Build::parallel_sync(), libMesh::MetisPartitioner::partition_range(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::SparseMatrix< T >::print(), libMesh::NumericVector< T >::print_global(), libMesh::Nemesis_IO_Helper::put_cmap_params(), libMesh::Nemesis_IO_Helper::put_elem_cmap(), libMesh::Nemesis_IO_Helper::put_elem_map(), libMesh::Nemesis_IO_Helper::put_loadbal_param(), libMesh::Nemesis_IO_Helper::put_node_cmap(), libMesh::Nemesis_IO_Helper::put_node_map(), libMesh::NameBasedIO::read(), libMesh::Nemesis_IO::read(), libMesh::XdrIO::read(), libMesh::CheckpointIO::read(), libMesh::ExodusII_IO_Helper::read_elem_num_map(), libMesh::ExodusII_IO_Helper::read_global_values(), libMesh::CheckpointIO::read_header(), libMesh::XdrIO::read_header(), libMesh::ExodusII_IO_Helper::read_node_num_map(), libMesh::XdrIO::read_serialized_bc_names(), libMesh::XdrIO::read_serialized_bcs_helper(), libMesh::XdrIO::read_serialized_connectivity(), libMesh::XdrIO::read_serialized_nodes(), libMesh::XdrIO::read_serialized_nodesets(), libMesh::XdrIO::read_serialized_subdomain_names(), libMesh::MeshCommunication::redistribute(), libMesh::DistributedMesh::renumber_dof_objects(), libMesh::CheckpointIO::select_split_config(), libMesh::MeshCommunication::send_coarse_ghosts(), libMesh::Partitioner::set_node_processor_ids(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::LaplaceMeshSmoother::smooth(), libMesh::MeshTools::total_weight(), libMesh::MeshRefinement::uniformly_coarsen(), libMesh::Parallel::Packing< Node * >::unpack(), libMesh::Parallel::Packing< Elem * >::unpack(), libMesh::DistributedMesh::update_parallel_id_counts(), libMesh::MeshTools::weight(), libMesh::NameBasedIO::write(), libMesh::XdrIO::write(), libMesh::CheckpointIO::write(), libMesh::EquationSystems::write(), libMesh::GMVIO::write_discontinuous_gmv(), libMesh::ExodusII_IO::write_element_data(), libMesh::ExodusII_IO_Helper::write_element_values(), libMesh::ExodusII_IO_Helper::write_elements(), libMesh::ExodusII_IO::write_global_data(), libMesh::ExodusII_IO_Helper::write_global_values(), libMesh::ExodusII_IO::write_information_records(), libMesh::ExodusII_IO_Helper::write_information_records(), libMesh::ExodusII_IO_Helper::write_nodal_coordinates(), libMesh::UCDIO::write_nodal_data(), libMesh::ExodusII_IO::write_nodal_data(), libMesh::ExodusII_IO::write_nodal_data_discontinuous(), libMesh::ExodusII_IO_Helper::write_nodal_values(), libMesh::Nemesis_IO_Helper::write_nodesets(), libMesh::ExodusII_IO_Helper::write_nodesets(), libMesh::XdrIO::write_serialized_bc_names(), libMesh::XdrIO::write_serialized_bcs_helper(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::XdrIO::write_serialized_nodes(), libMesh::XdrIO::write_serialized_nodesets(), libMesh::XdrIO::write_serialized_subdomain_names(), libMesh::Nemesis_IO_Helper::write_sidesets(), libMesh::ExodusII_IO_Helper::write_sidesets(), libMesh::ExodusII_IO::write_timestep(), libMesh::ExodusII_IO_Helper::write_timestep(), and libMesh::ExodusII_IO::write_timestep_discontinuous().

100  { return cast_int<processor_id_type>(_communicator.rank()); }
const Parallel::Communicator & _communicator
unsigned int rank() const
Definition: parallel.h:723
template<typename T >
void libMesh::PetscLinearSolver< T >::restrict_solve_to ( const std::vector< unsigned int > *const  dofs,
const SubsetSolveMode  subset_solve_mode = SUBSET_ZERO 
)
virtual

After calling this method, all successive solves will be restricted to the given set of dofs, which must contain local dofs on each processor only and not contain any duplicates. This mode can be disabled by calling this method with dofs being a NULL pointer.

Reimplemented from libMesh::LinearSolver< T >.

Definition at line 352 of file petsc_linear_solver.C.

References ierr, libmesh_nullptr, and PETSC_OWN_POINTER.

354 {
355  PetscErrorCode ierr=0;
356 
357  /* The preconditioner (in particular if a default preconditioner)
358  will have to be reset. We call this->clear() to do that. This
359  call will also remove and free any previous subset that may have
360  been set before. */
361  this->clear();
362 
363  _subset_solve_mode = subset_solve_mode;
364 
365  if (dofs != libmesh_nullptr)
366  {
367  PetscInt * petsc_dofs = libmesh_nullptr;
368  ierr = PetscMalloc(dofs->size()*sizeof(PetscInt), &petsc_dofs);
369  LIBMESH_CHKERR(ierr);
370 
371  for (std::size_t i=0; i<dofs->size(); i++)
372  petsc_dofs[i] = (*dofs)[i];
373 
374  ierr = ISCreateLibMesh(this->comm().get(),dofs->size(),petsc_dofs,PETSC_OWN_POINTER,&_restrict_solve_to_is);
375  LIBMESH_CHKERR(ierr);
376  }
377 }
virtual void clear() libmesh_override
unsigned int size() const
Definition: parallel.h:725
const class libmesh_nullptr_t libmesh_nullptr
PetscErrorCode ierr
const Parallel::Communicator & comm() const
template<typename T >
void libMesh::LinearSolver< T >::reuse_preconditioner ( bool  reuse_flag)
virtualinherited

Set the same_preconditioner flag, which indicates if we reuse the same preconditioner for subsequent solves.

Definition at line 112 of file linear_solver.C.

Referenced by libMesh::ImplicitSystem::disable_cache(), and libMesh::LinearSolver< Number >::set_solver_type().

113 {
114  same_preconditioner = reuse_flag;
115 }
template<typename T >
void libMesh::PetscLinearSolver< T >::set_petsc_solver_type ( )
private

Tells PETSC to use the user-specified solver stored in _solver_type

Definition at line 1550 of file petsc_linear_solver.C.

References libMesh::BICG, libMesh::BICGSTAB, libMesh::CG, libMesh::CGS, libMesh::CHEBYSHEV, libMesh::CR, libMesh::Utility::enum_to_string(), libMesh::err, libMesh::GMRES, ierr, libMesh::LSQR, libMesh::MINRES, libMesh::RICHARDSON, libMesh::TCQMR, and libMesh::TFQMR.

Referenced by libMesh::PetscLinearSolver< T >::ksp().

1551 {
1552  PetscErrorCode ierr = 0;
1553 
1554  switch (this->_solver_type)
1555  {
1556 
1557  case CG:
1558  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPCG));
1559  LIBMESH_CHKERR(ierr);
1560  return;
1561 
1562  case CR:
1563  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPCR));
1564  LIBMESH_CHKERR(ierr);
1565  return;
1566 
1567  case CGS:
1568  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPCGS));
1569  LIBMESH_CHKERR(ierr);
1570  return;
1571 
1572  case BICG:
1573  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPBICG));
1574  LIBMESH_CHKERR(ierr);
1575  return;
1576 
1577  case TCQMR:
1578  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPTCQMR));
1579  LIBMESH_CHKERR(ierr);
1580  return;
1581 
1582  case TFQMR:
1583  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPTFQMR));
1584  LIBMESH_CHKERR(ierr);
1585  return;
1586 
1587  case LSQR:
1588  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPLSQR));
1589  LIBMESH_CHKERR(ierr);
1590  return;
1591 
1592  case BICGSTAB:
1593  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPBCGS));
1594  LIBMESH_CHKERR(ierr);
1595  return;
1596 
1597  case MINRES:
1598  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPMINRES));
1599  LIBMESH_CHKERR(ierr);
1600  return;
1601 
1602  case GMRES:
1603  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPGMRES));
1604  LIBMESH_CHKERR(ierr);
1605  return;
1606 
1607  case RICHARDSON:
1608  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPRICHARDSON));
1609  LIBMESH_CHKERR(ierr);
1610  return;
1611 
1612  case CHEBYSHEV:
1613 #if defined(LIBMESH_HAVE_PETSC) && PETSC_VERSION_LESS_THAN(3,3,0)
1614  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPCHEBYCHEV));
1615  LIBMESH_CHKERR(ierr);
1616  return;
1617 #else
1618  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPCHEBYSHEV));
1619  LIBMESH_CHKERR(ierr);
1620  return;
1621 #endif
1622 
1623 
1624  default:
1625  libMesh::err << "ERROR: Unsupported PETSC Solver: "
1626  << Utility::enum_to_string(this->_solver_type) << std::endl
1627  << "Continuing with PETSC defaults" << std::endl;
1628  }
1629 }
OStreamProxy err(std::cerr)
std::string enum_to_string(const T e)
PetscErrorCode ierr
template<typename T >
void libMesh::LinearSolver< T >::set_preconditioner_type ( const PreconditionerType  pct)
inherited

Sets the type of preconditioner to use.

Definition at line 91 of file linear_solver.C.

Referenced by libMesh::LinearSolver< Number >::set_solver_type().

92 {
93  if (_preconditioner)
94  _preconditioner->set_type(pct);
95  else
97 }
Preconditioner< T > * _preconditioner
PreconditionerType _preconditioner_type
template<typename T >
void libMesh::LinearSolver< T >::set_solver_configuration ( SolverConfiguration solver_configuration)
inherited

Set the solver configuration object.

Definition at line 160 of file linear_solver.C.

Referenced by libMesh::LinearSolver< Number >::set_solver_type().

161 {
162  _solver_configuration = &solver_configuration;
163 }
SolverConfiguration * _solver_configuration
template<typename T>
void libMesh::LinearSolver< T >::set_solver_type ( const SolverType  st)
inlineinherited

Sets the type of solver to use.

Definition at line 117 of file linear_solver.h.

118  { _solver_type = st; }
template<typename T >
virtual std::pair<unsigned int, Real> libMesh::PetscLinearSolver< T >::solve ( SparseMatrix< T > &  matrix_in,
NumericVector< T > &  solution_in,
NumericVector< T > &  rhs_in,
const double  tol,
const unsigned int  m_its 
)
inlinevirtual

Call the Petsc solver. It calls the method below, using the same matrix for the system and preconditioner matrices.

Implements libMesh::LinearSolver< T >.

Definition at line 139 of file petsc_linear_solver.h.

References libMesh::PetscLinearSolver< T >::adjoint_solve().

144  {
145  return this->solve(matrix_in, matrix_in, solution_in, rhs_in, tol, m_its);
146  }
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
template<typename T >
std::pair< unsigned int, Real > libMesh::PetscLinearSolver< T >::solve ( SparseMatrix< T > &  matrix,
SparseMatrix< T > &  preconditioner,
NumericVector< T > &  solution,
NumericVector< T > &  rhs,
const double  tol,
const unsigned int  m_its 
)
virtual

This method allows you to call a linear solver while specifying the matrix to use as the (left) preconditioning matrix.

Note
The linear solver will not compute a preconditioner in this case, and will instead premultiply by the matrix you provide. In PETSc, this is accomplished by calling
PCSetType(_pc, PCMAT);
before invoking KSPSolve().
This functionality is not implemented in the LinearSolver class since there is not a built-in analog to this method for LASPACK. You could probably implement it by hand if you wanted.

Implements libMesh::LinearSolver< T >.

Definition at line 383 of file petsc_linear_solver.C.

References libMesh::PetscMatrix< T >::close(), ierr, libMesh::TriangleWrapper::init(), libmesh_nullptr, libMesh::NumericVector< T >::local_size(), libMesh::PetscMatrix< T >::mat(), libMesh::SUBSET_COPY_RHS, libMesh::SUBSET_DONT_TOUCH, and libMesh::SUBSET_ZERO.

389 {
390  LOG_SCOPE("solve()", "PetscLinearSolver");
391 
392  // Make sure the data passed in are really of Petsc types
393  PetscMatrix<T> * matrix = cast_ptr<PetscMatrix<T> *>(&matrix_in);
394  PetscMatrix<T> * precond = cast_ptr<PetscMatrix<T> *>(&precond_in);
395  PetscVector<T> * solution = cast_ptr<PetscVector<T> *>(&solution_in);
396  PetscVector<T> * rhs = cast_ptr<PetscVector<T> *>(&rhs_in);
397 
398  this->init (matrix);
399 
400  PetscErrorCode ierr=0;
401  PetscInt its=0, max_its = static_cast<PetscInt>(m_its);
402  PetscReal final_resid=0.;
403 
404  // Close the matrices and vectors in case this wasn't already done.
405  matrix->close ();
406  precond->close ();
407  solution->close ();
408  rhs->close ();
409 
410  // // If matrix != precond, then this means we have specified a
411  // // special preconditioner, so reset preconditioner type to PCMAT.
412  // if (matrix != precond)
413  // {
414  // this->_preconditioner_type = USER_PRECOND;
415  // this->set_petsc_preconditioner_type ();
416  // }
417 
418  Mat submat = libmesh_nullptr;
419  Mat subprecond = libmesh_nullptr;
420  Vec subrhs = libmesh_nullptr;
421  Vec subsolution = libmesh_nullptr;
422  VecScatter scatter = libmesh_nullptr;
423  std::unique_ptr<PetscMatrix<Number>> subprecond_matrix;
424 
425  // Set operators. Also restrict rhs and solution vector to
426  // subdomain if necessary.
428  {
429  PetscInt is_local_size = this->_restrict_solve_to_is_local_size();
430 
431  ierr = VecCreate(this->comm().get(),&subrhs);
432  LIBMESH_CHKERR(ierr);
433  ierr = VecSetSizes(subrhs,is_local_size,PETSC_DECIDE);
434  LIBMESH_CHKERR(ierr);
435  ierr = VecSetFromOptions(subrhs);
436  LIBMESH_CHKERR(ierr);
437 
438  ierr = VecCreate(this->comm().get(),&subsolution);
439  LIBMESH_CHKERR(ierr);
440  ierr = VecSetSizes(subsolution,is_local_size,PETSC_DECIDE);
441  LIBMESH_CHKERR(ierr);
442  ierr = VecSetFromOptions(subsolution);
443  LIBMESH_CHKERR(ierr);
444 
445  ierr = VecScatterCreate(rhs->vec(), _restrict_solve_to_is, subrhs, libmesh_nullptr, &scatter);
446  LIBMESH_CHKERR(ierr);
447 
448  ierr = VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
449  LIBMESH_CHKERR(ierr);
450  ierr = VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
451  LIBMESH_CHKERR(ierr);
452 
453  ierr = VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
454  LIBMESH_CHKERR(ierr);
455  ierr = VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
456  LIBMESH_CHKERR(ierr);
457 
458  ierr = LibMeshCreateSubMatrix(matrix->mat(),
461 #if PETSC_VERSION_LESS_THAN(3,1,0)
462  PETSC_DECIDE,
463 #endif
464  MAT_INITIAL_MATRIX,
465  &submat);
466  LIBMESH_CHKERR(ierr);
467 
468  ierr = LibMeshCreateSubMatrix(precond->mat(),
471 #if PETSC_VERSION_LESS_THAN(3,1,0)
472  PETSC_DECIDE,
473 #endif
474  MAT_INITIAL_MATRIX,
475  &subprecond);
476  LIBMESH_CHKERR(ierr);
477 
478  /* Since removing columns of the matrix changes the equation
479  system, we will now change the right hand side to compensate
480  for this. Note that this is not necessary if \p SUBSET_ZERO
481  has been selected. */
483  {
484  _create_complement_is(rhs_in);
485  PetscInt is_complement_local_size =
486  cast_int<PetscInt>(rhs_in.local_size()-is_local_size);
487 
488  Vec subvec1 = libmesh_nullptr;
489  Mat submat1 = libmesh_nullptr;
490  VecScatter scatter1 = libmesh_nullptr;
491 
492  ierr = VecCreate(this->comm().get(),&subvec1);
493  LIBMESH_CHKERR(ierr);
494  ierr = VecSetSizes(subvec1,is_complement_local_size,PETSC_DECIDE);
495  LIBMESH_CHKERR(ierr);
496  ierr = VecSetFromOptions(subvec1);
497  LIBMESH_CHKERR(ierr);
498 
499  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is_complement, subvec1, libmesh_nullptr, &scatter1);
500  LIBMESH_CHKERR(ierr);
501 
502  ierr = VecScatterBegin(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
503  LIBMESH_CHKERR(ierr);
504  ierr = VecScatterEnd(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
505  LIBMESH_CHKERR(ierr);
506 
507  ierr = VecScale(subvec1,-1.0);
508  LIBMESH_CHKERR(ierr);
509 
510  ierr = LibMeshCreateSubMatrix(matrix->mat(),
513 #if PETSC_VERSION_LESS_THAN(3,1,0)
514  PETSC_DECIDE,
515 #endif
516  MAT_INITIAL_MATRIX,
517  &submat1);
518  LIBMESH_CHKERR(ierr);
519 
520  ierr = MatMultAdd(submat1,subvec1,subrhs,subrhs);
521  LIBMESH_CHKERR(ierr);
522 
523  ierr = LibMeshVecScatterDestroy(&scatter1);
524  LIBMESH_CHKERR(ierr);
525  ierr = LibMeshVecDestroy(&subvec1);
526  LIBMESH_CHKERR(ierr);
527  ierr = LibMeshMatDestroy(&submat1);
528  LIBMESH_CHKERR(ierr);
529  }
530 #if PETSC_RELEASE_LESS_THAN(3,5,0)
531  ierr = KSPSetOperators(_ksp, submat, subprecond,
532  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
533 #else
534  ierr = KSPSetOperators(_ksp, submat, subprecond);
535 
536  PetscBool ksp_reuse_preconditioner = this->same_preconditioner ? PETSC_TRUE : PETSC_FALSE;
537  ierr = KSPSetReusePreconditioner(_ksp, ksp_reuse_preconditioner);
538 #endif
539  LIBMESH_CHKERR(ierr);
540 
541  if (this->_preconditioner)
542  {
543  subprecond_matrix.reset(new PetscMatrix<Number>(subprecond, this->comm()));
544  this->_preconditioner->set_matrix(*subprecond_matrix);
545  this->_preconditioner->init();
546  }
547  }
548  else
549  {
550 #if PETSC_RELEASE_LESS_THAN(3,5,0)
551  ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
552  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
553 #else
554  ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat());
555 
556  PetscBool ksp_reuse_preconditioner = this->same_preconditioner ? PETSC_TRUE : PETSC_FALSE;
557  ierr = KSPSetReusePreconditioner(_ksp, ksp_reuse_preconditioner);
558 #endif
559  LIBMESH_CHKERR(ierr);
560 
561  if (this->_preconditioner)
562  {
563  this->_preconditioner->set_matrix(matrix_in);
564  this->_preconditioner->init();
565  }
566  }
567 
568  // Set the tolerances for the iterative solver. Use the user-supplied
569  // tolerance for the relative residual & leave the others at default values.
570  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
571  PETSC_DEFAULT, max_its);
572  LIBMESH_CHKERR(ierr);
573 
574  // Allow command line options to override anything set programmatically.
575  ierr = KSPSetFromOptions(_ksp);
576  LIBMESH_CHKERR(ierr);
577 
578  // If the SolverConfiguration object is provided, use it to override
579  // solver options.
580  if (this->_solver_configuration)
581  {
583  }
584 
585  // Solve the linear system
586  if (_restrict_solve_to_is != libmesh_nullptr)
587  {
588  ierr = KSPSolve (_ksp, subrhs, subsolution);
589  LIBMESH_CHKERR(ierr);
590  }
591  else
592  {
593  ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
594  LIBMESH_CHKERR(ierr);
595  }
596 
597  // Get the number of iterations required for convergence
598  ierr = KSPGetIterationNumber (_ksp, &its);
599  LIBMESH_CHKERR(ierr);
600 
601  // Get the norm of the final residual to return to the user.
602  ierr = KSPGetResidualNorm (_ksp, &final_resid);
603  LIBMESH_CHKERR(ierr);
604 
605  if (_restrict_solve_to_is != libmesh_nullptr)
606  {
607  switch(_subset_solve_mode)
608  {
609  case SUBSET_ZERO:
610  ierr = VecZeroEntries(solution->vec());
611  LIBMESH_CHKERR(ierr);
612  break;
613 
614  case SUBSET_COPY_RHS:
615  ierr = VecCopy(rhs->vec(),solution->vec());
616  LIBMESH_CHKERR(ierr);
617  break;
618 
619  case SUBSET_DONT_TOUCH:
620  /* Nothing to do here. */
621  break;
622 
623  default:
624  libmesh_error_msg("Invalid subset solve mode = " << _subset_solve_mode);
625  }
626  ierr = VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
627  LIBMESH_CHKERR(ierr);
628  ierr = VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
629  LIBMESH_CHKERR(ierr);
630 
631  ierr = LibMeshVecScatterDestroy(&scatter);
632  LIBMESH_CHKERR(ierr);
633 
634  if (this->_preconditioner)
635  {
636  // Before subprecond_matrix gets cleaned up, we should give
637  // the _preconditioner a different matrix.
638  this->_preconditioner->set_matrix(matrix_in);
639  this->_preconditioner->init();
640  }
641 
642  ierr = LibMeshVecDestroy(&subsolution);
643  LIBMESH_CHKERR(ierr);
644  ierr = LibMeshVecDestroy(&subrhs);
645  LIBMESH_CHKERR(ierr);
646  ierr = LibMeshMatDestroy(&submat);
647  LIBMESH_CHKERR(ierr);
648  ierr = LibMeshMatDestroy(&subprecond);
649  LIBMESH_CHKERR(ierr);
650  }
651 
652  // return the # of its. and the final residual norm.
653  return std::make_pair(its, final_resid);
654 }
Leaves dofs outside the subset unchanged. This is fastest, but also most confusing because it abandon...
Preconditioner< T > * _preconditioner
Set dofs outside the subset to zero.
void _create_complement_is(const NumericVector< T > &vec_in)
PetscInt _restrict_solve_to_is_local_size() const
const class libmesh_nullptr_t libmesh_nullptr
virtual void configure_solver()=0
SolverConfiguration * _solver_configuration
virtual void init(const char *name=libmesh_nullptr) libmesh_override
Set dofs outside the subset to the value of the corresponding dofs of the right hand side...
PetscTruth PetscBool
Definition: petsc_macro.h:64
PetscErrorCode ierr
const Parallel::Communicator & comm() const
template<typename T >
std::pair< unsigned int, Real > libMesh::PetscLinearSolver< T >::solve ( const ShellMatrix< T > &  shell_matrix,
NumericVector< T > &  solution_in,
NumericVector< T > &  rhs_in,
const double  tol,
const unsigned int  m_its 
)
virtual

This function solves a system whose matrix is a shell matrix.

Implements libMesh::LinearSolver< T >.

Definition at line 919 of file petsc_linear_solver.C.

References libMesh::PetscVector< T >::close(), ierr, libMesh::TriangleWrapper::init(), libmesh_nullptr, libMesh::NumericVector< T >::local_size(), libMesh::NumericVector< T >::size(), libMesh::SUBSET_COPY_RHS, libMesh::SUBSET_DONT_TOUCH, libMesh::SUBSET_ZERO, and libMesh::PetscVector< T >::vec().

924 {
925 
926 #if PETSC_VERSION_LESS_THAN(3,1,0)
928  libmesh_error_msg("The current implementation of subset solves with " \
929  << "shell matrices requires PETSc version 3.1 or above. " \
930  << "Older PETSc version do not support automatic " \
931  << "submatrix generation of shell matrices.");
932 #endif
933 
934  LOG_SCOPE("solve()", "PetscLinearSolver");
935 
936  // Make sure the data passed in are really of Petsc types
937  PetscVector<T> * solution = cast_ptr<PetscVector<T> *>(&solution_in);
938  PetscVector<T> * rhs = cast_ptr<PetscVector<T> *>(&rhs_in);
939 
940  this->init ();
941 
942  PetscErrorCode ierr=0;
943  PetscInt its=0, max_its = static_cast<PetscInt>(m_its);
944  PetscReal final_resid=0.;
945 
946  Mat submat = libmesh_nullptr;
947  Vec subrhs = libmesh_nullptr;
948  Vec subsolution = libmesh_nullptr;
949  VecScatter scatter = libmesh_nullptr;
950 
951  // Close the matrices and vectors in case this wasn't already done.
952  solution->close ();
953  rhs->close ();
954 
955  // Prepare the matrix.
956  Mat mat;
957  ierr = MatCreateShell(this->comm().get(),
958  rhs_in.local_size(),
959  solution_in.local_size(),
960  rhs_in.size(),
961  solution_in.size(),
962  const_cast<void *>(static_cast<const void *>(&shell_matrix)),
963  &mat);
964  /* Note that the const_cast above is only necessary because PETSc
965  does not accept a const void *. Inside the member function
966  _petsc_shell_matrix() below, the pointer is casted back to a
967  const ShellMatrix<T> *. */
968 
969  LIBMESH_CHKERR(ierr);
970  ierr = MatShellSetOperation(mat,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
971  LIBMESH_CHKERR(ierr);
972  ierr = MatShellSetOperation(mat,MATOP_MULT_ADD,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult_add));
973  LIBMESH_CHKERR(ierr);
974  ierr = MatShellSetOperation(mat,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
975  LIBMESH_CHKERR(ierr);
976 
977  // Restrict rhs and solution vectors and set operators. The input
978  // matrix works as the preconditioning matrix.
980  {
981  PetscInt is_local_size = this->_restrict_solve_to_is_local_size();
982 
983  ierr = VecCreate(this->comm().get(),&subrhs);
984  LIBMESH_CHKERR(ierr);
985  ierr = VecSetSizes(subrhs,is_local_size,PETSC_DECIDE);
986  LIBMESH_CHKERR(ierr);
987  ierr = VecSetFromOptions(subrhs);
988  LIBMESH_CHKERR(ierr);
989 
990  ierr = VecCreate(this->comm().get(),&subsolution);
991  LIBMESH_CHKERR(ierr);
992  ierr = VecSetSizes(subsolution,is_local_size,PETSC_DECIDE);
993  LIBMESH_CHKERR(ierr);
994  ierr = VecSetFromOptions(subsolution);
995  LIBMESH_CHKERR(ierr);
996 
997  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is, subrhs, libmesh_nullptr, &scatter);
998  LIBMESH_CHKERR(ierr);
999 
1000  ierr = VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
1001  LIBMESH_CHKERR(ierr);
1002  ierr = VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
1003  LIBMESH_CHKERR(ierr);
1004 
1005  ierr = VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
1006  LIBMESH_CHKERR(ierr);
1007  ierr = VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
1008  LIBMESH_CHKERR(ierr);
1009 
1010 #if !PETSC_VERSION_LESS_THAN(3,1,0)
1011  ierr = LibMeshCreateSubMatrix(mat,
1014  MAT_INITIAL_MATRIX,
1015  &submat);
1016  LIBMESH_CHKERR(ierr);
1017 #endif
1018 
1019  /* Since removing columns of the matrix changes the equation
1020  system, we will now change the right hand side to compensate
1021  for this. Note that this is not necessary if \p SUBSET_ZERO
1022  has been selected. */
1024  {
1025  _create_complement_is(rhs_in);
1026  PetscInt is_complement_local_size =
1027  cast_int<PetscInt>(rhs_in.local_size()-is_local_size);
1028 
1029  Vec subvec1 = libmesh_nullptr;
1030  Mat submat1 = libmesh_nullptr;
1031  VecScatter scatter1 = libmesh_nullptr;
1032 
1033  ierr = VecCreate(this->comm().get(),&subvec1);
1034  LIBMESH_CHKERR(ierr);
1035  ierr = VecSetSizes(subvec1,is_complement_local_size,PETSC_DECIDE);
1036  LIBMESH_CHKERR(ierr);
1037  ierr = VecSetFromOptions(subvec1);
1038  LIBMESH_CHKERR(ierr);
1039 
1040  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is_complement, subvec1, libmesh_nullptr, &scatter1);
1041  LIBMESH_CHKERR(ierr);
1042 
1043  ierr = VecScatterBegin(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
1044  LIBMESH_CHKERR(ierr);
1045  ierr = VecScatterEnd(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
1046  LIBMESH_CHKERR(ierr);
1047 
1048  ierr = VecScale(subvec1,-1.0);
1049  LIBMESH_CHKERR(ierr);
1050 
1051 #if !PETSC_VERSION_LESS_THAN(3,1,0)
1052  ierr = LibMeshCreateSubMatrix(mat,
1055  MAT_INITIAL_MATRIX,
1056  &submat1);
1057  LIBMESH_CHKERR(ierr);
1058 #endif
1059 
1060  // The following lines would be correct, but don't work
1061  // correctly in PETSc up to 3.1.0-p5. See discussion in
1062  // petsc-users of Nov 9, 2010.
1063  //
1064  // ierr = MatMultAdd(submat1,subvec1,subrhs,subrhs);
1065  // LIBMESH_CHKERR(ierr);
1066  //
1067  // We workaround by using a temporary vector. Note that the
1068  // fix in PETsc 3.1.0-p6 uses a temporary vector internally,
1069  // so this is no effective performance loss.
1070  Vec subvec2 = libmesh_nullptr;
1071  ierr = VecCreate(this->comm().get(),&subvec2);
1072  LIBMESH_CHKERR(ierr);
1073  ierr = VecSetSizes(subvec2,is_local_size,PETSC_DECIDE);
1074  LIBMESH_CHKERR(ierr);
1075  ierr = VecSetFromOptions(subvec2);
1076  LIBMESH_CHKERR(ierr);
1077  ierr = MatMult(submat1,subvec1,subvec2);
1078  LIBMESH_CHKERR(ierr);
1079  ierr = VecAXPY(subrhs,1.0,subvec2);
1080 
1081  ierr = LibMeshVecScatterDestroy(&scatter1);
1082  LIBMESH_CHKERR(ierr);
1083  ierr = LibMeshVecDestroy(&subvec1);
1084  LIBMESH_CHKERR(ierr);
1085  ierr = LibMeshMatDestroy(&submat1);
1086  LIBMESH_CHKERR(ierr);
1087  }
1088 #if PETSC_RELEASE_LESS_THAN(3,5,0)
1089  ierr = KSPSetOperators(_ksp, submat, submat,
1090  DIFFERENT_NONZERO_PATTERN);
1091 #else
1092  ierr = KSPSetOperators(_ksp, submat, submat);
1093 #endif
1094  LIBMESH_CHKERR(ierr);
1095  }
1096  else
1097  {
1098 #if PETSC_RELEASE_LESS_THAN(3,5,0)
1099  ierr = KSPSetOperators(_ksp, mat, mat,
1100  DIFFERENT_NONZERO_PATTERN);
1101 #else
1102  ierr = KSPSetOperators(_ksp, mat, mat);
1103 #endif
1104  LIBMESH_CHKERR(ierr);
1105  }
1106 
1107  // Set the tolerances for the iterative solver. Use the user-supplied
1108  // tolerance for the relative residual & leave the others at default values.
1109  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
1110  PETSC_DEFAULT, max_its);
1111  LIBMESH_CHKERR(ierr);
1112 
1113  // Allow command line options to override anything set programmatically.
1114  ierr = KSPSetFromOptions(_ksp);
1115  LIBMESH_CHKERR(ierr);
1116 
1117  // Solve the linear system
1118  if (_restrict_solve_to_is != libmesh_nullptr)
1119  {
1120  ierr = KSPSolve (_ksp, subrhs, subsolution);
1121  LIBMESH_CHKERR(ierr);
1122  }
1123  else
1124  {
1125  ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
1126  LIBMESH_CHKERR(ierr);
1127  }
1128 
1129  // Get the number of iterations required for convergence
1130  ierr = KSPGetIterationNumber (_ksp, &its);
1131  LIBMESH_CHKERR(ierr);
1132 
1133  // Get the norm of the final residual to return to the user.
1134  ierr = KSPGetResidualNorm (_ksp, &final_resid);
1135  LIBMESH_CHKERR(ierr);
1136 
1137  if (_restrict_solve_to_is != libmesh_nullptr)
1138  {
1139  switch(_subset_solve_mode)
1140  {
1141  case SUBSET_ZERO:
1142  ierr = VecZeroEntries(solution->vec());
1143  LIBMESH_CHKERR(ierr);
1144  break;
1145 
1146  case SUBSET_COPY_RHS:
1147  ierr = VecCopy(rhs->vec(),solution->vec());
1148  LIBMESH_CHKERR(ierr);
1149  break;
1150 
1151  case SUBSET_DONT_TOUCH:
1152  /* Nothing to do here. */
1153  break;
1154 
1155  default:
1156  libmesh_error_msg("Invalid subset solve mode = " << _subset_solve_mode);
1157  }
1158  ierr = VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1159  LIBMESH_CHKERR(ierr);
1160  ierr = VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1161  LIBMESH_CHKERR(ierr);
1162 
1163  ierr = LibMeshVecScatterDestroy(&scatter);
1164  LIBMESH_CHKERR(ierr);
1165 
1166  ierr = LibMeshVecDestroy(&subsolution);
1167  LIBMESH_CHKERR(ierr);
1168  ierr = LibMeshVecDestroy(&subrhs);
1169  LIBMESH_CHKERR(ierr);
1170  ierr = LibMeshMatDestroy(&submat);
1171  LIBMESH_CHKERR(ierr);
1172  }
1173 
1174  // Destroy the matrix.
1175  ierr = LibMeshMatDestroy(&mat);
1176  LIBMESH_CHKERR(ierr);
1177 
1178  // return the # of its. and the final residual norm.
1179  return std::make_pair(its, final_resid);
1180 }
Leaves dofs outside the subset unchanged. This is fastest, but also most confusing because it abandon...
unsigned int size() const
Definition: parallel.h:725
Set dofs outside the subset to zero.
void _create_complement_is(const NumericVector< T > &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)
virtual void init(const char *name=libmesh_nullptr) libmesh_override
Set dofs outside the subset to the value of the corresponding dofs of the right hand side...
static PetscErrorCode _petsc_shell_matrix_mult(Mat mat, Vec arg, Vec dest)
PetscErrorCode ierr
const Parallel::Communicator & comm() const
static PetscErrorCode _petsc_shell_matrix_get_diagonal(Mat mat, Vec dest)
template<typename T >
std::pair< unsigned int, Real > libMesh::PetscLinearSolver< T >::solve ( const ShellMatrix< T > &  shell_matrix,
const SparseMatrix< T > &  precond_matrix,
NumericVector< T > &  solution_in,
NumericVector< T > &  rhs_in,
const double  tol,
const unsigned int  m_its 
)
virtual

This function solves a system whose matrix is a shell matrix, but a sparse matrix is used as preconditioning matrix, this allowing other preconditioners than JACOBI.

Implements libMesh::LinearSolver< T >.

Definition at line 1186 of file petsc_linear_solver.C.

References ierr, libMesh::TriangleWrapper::init(), libmesh_nullptr, libMesh::NumericVector< T >::local_size(), libMesh::NumericVector< T >::size(), libMesh::SUBSET_COPY_RHS, libMesh::SUBSET_DONT_TOUCH, and libMesh::SUBSET_ZERO.

1192 {
1193 
1194 #if PETSC_VERSION_LESS_THAN(3,1,0)
1196  libmesh_error_msg("The current implementation of subset solves with " \
1197  << "shell matrices requires PETSc version 3.1 or above. " \
1198  << "Older PETSc version do not support automatic " \
1199  << "submatrix generation of shell matrices.");
1200 #endif
1201 
1202  LOG_SCOPE("solve()", "PetscLinearSolver");
1203 
1204  // Make sure the data passed in are really of Petsc types
1205  const PetscMatrix<T> * precond = cast_ptr<const PetscMatrix<T> *>(&precond_matrix);
1206  PetscVector<T> * solution = cast_ptr<PetscVector<T> *>(&solution_in);
1207  PetscVector<T> * rhs = cast_ptr<PetscVector<T> *>(&rhs_in);
1208 
1209  this->init ();
1210 
1211  PetscErrorCode ierr=0;
1212  PetscInt its=0, max_its = static_cast<PetscInt>(m_its);
1213  PetscReal final_resid=0.;
1214 
1215  Mat submat = libmesh_nullptr;
1216  Mat subprecond = libmesh_nullptr;
1217  Vec subrhs = libmesh_nullptr;
1218  Vec subsolution = libmesh_nullptr;
1219  VecScatter scatter = libmesh_nullptr;
1220  std::unique_ptr<PetscMatrix<Number>> subprecond_matrix;
1221 
1222  // Close the matrices and vectors in case this wasn't already done.
1223  solution->close ();
1224  rhs->close ();
1225 
1226  // Prepare the matrix.
1227  Mat mat;
1228  ierr = MatCreateShell(this->comm().get(),
1229  rhs_in.local_size(),
1230  solution_in.local_size(),
1231  rhs_in.size(),
1232  solution_in.size(),
1233  const_cast<void *>(static_cast<const void *>(&shell_matrix)),
1234  &mat);
1235  /* Note that the const_cast above is only necessary because PETSc
1236  does not accept a const void *. Inside the member function
1237  _petsc_shell_matrix() below, the pointer is casted back to a
1238  const ShellMatrix<T> *. */
1239 
1240  LIBMESH_CHKERR(ierr);
1241  ierr = MatShellSetOperation(mat,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
1242  LIBMESH_CHKERR(ierr);
1243  ierr = MatShellSetOperation(mat,MATOP_MULT_ADD,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult_add));
1244  LIBMESH_CHKERR(ierr);
1245  ierr = MatShellSetOperation(mat,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
1246  LIBMESH_CHKERR(ierr);
1247 
1248  // Restrict rhs and solution vectors and set operators. The input
1249  // matrix works as the preconditioning matrix.
1251  {
1252  PetscInt is_local_size = this->_restrict_solve_to_is_local_size();
1253 
1254  ierr = VecCreate(this->comm().get(),&subrhs);
1255  LIBMESH_CHKERR(ierr);
1256  ierr = VecSetSizes(subrhs,is_local_size,PETSC_DECIDE);
1257  LIBMESH_CHKERR(ierr);
1258  ierr = VecSetFromOptions(subrhs);
1259  LIBMESH_CHKERR(ierr);
1260 
1261  ierr = VecCreate(this->comm().get(),&subsolution);
1262  LIBMESH_CHKERR(ierr);
1263  ierr = VecSetSizes(subsolution,is_local_size,PETSC_DECIDE);
1264  LIBMESH_CHKERR(ierr);
1265  ierr = VecSetFromOptions(subsolution);
1266  LIBMESH_CHKERR(ierr);
1267 
1268  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is, subrhs, libmesh_nullptr, &scatter);
1269  LIBMESH_CHKERR(ierr);
1270 
1271  ierr = VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
1272  LIBMESH_CHKERR(ierr);
1273  ierr = VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
1274  LIBMESH_CHKERR(ierr);
1275 
1276  ierr = VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
1277  LIBMESH_CHKERR(ierr);
1278  ierr = VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
1279  LIBMESH_CHKERR(ierr);
1280 
1281 #if !PETSC_VERSION_LESS_THAN(3,1,0)
1282  ierr = LibMeshCreateSubMatrix(mat,
1285  MAT_INITIAL_MATRIX,
1286  &submat);
1287  LIBMESH_CHKERR(ierr);
1288 
1289  ierr = LibMeshCreateSubMatrix(const_cast<PetscMatrix<T> *>(precond)->mat(),
1292  MAT_INITIAL_MATRIX,
1293  &subprecond);
1294  LIBMESH_CHKERR(ierr);
1295 #endif
1296 
1297  /* Since removing columns of the matrix changes the equation
1298  system, we will now change the right hand side to compensate
1299  for this. Note that this is not necessary if \p SUBSET_ZERO
1300  has been selected. */
1302  {
1303  _create_complement_is(rhs_in);
1304  PetscInt is_complement_local_size = rhs_in.local_size()-is_local_size;
1305 
1306  Vec subvec1 = libmesh_nullptr;
1307  Mat submat1 = libmesh_nullptr;
1308  VecScatter scatter1 = libmesh_nullptr;
1309 
1310  ierr = VecCreate(this->comm().get(),&subvec1);
1311  LIBMESH_CHKERR(ierr);
1312  ierr = VecSetSizes(subvec1,is_complement_local_size,PETSC_DECIDE);
1313  LIBMESH_CHKERR(ierr);
1314  ierr = VecSetFromOptions(subvec1);
1315  LIBMESH_CHKERR(ierr);
1316 
1317  ierr = VecScatterCreate(rhs->vec(),_restrict_solve_to_is_complement, subvec1, libmesh_nullptr, &scatter1);
1318  LIBMESH_CHKERR(ierr);
1319 
1320  ierr = VecScatterBegin(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
1321  LIBMESH_CHKERR(ierr);
1322  ierr = VecScatterEnd(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
1323  LIBMESH_CHKERR(ierr);
1324 
1325  ierr = VecScale(subvec1,-1.0);
1326  LIBMESH_CHKERR(ierr);
1327 
1328 #if !PETSC_VERSION_LESS_THAN(3,1,0)
1329  ierr = LibMeshCreateSubMatrix(mat,
1332  MAT_INITIAL_MATRIX,
1333  &submat1);
1334  LIBMESH_CHKERR(ierr);
1335 #endif
1336 
1337  // The following lines would be correct, but don't work
1338  // correctly in PETSc up to 3.1.0-p5. See discussion in
1339  // petsc-users of Nov 9, 2010.
1340  //
1341  // ierr = MatMultAdd(submat1,subvec1,subrhs,subrhs);
1342  // LIBMESH_CHKERR(ierr);
1343  //
1344  // We workaround by using a temporary vector. Note that the
1345  // fix in PETsc 3.1.0-p6 uses a temporary vector internally,
1346  // so this is no effective performance loss.
1347  Vec subvec2 = libmesh_nullptr;
1348  ierr = VecCreate(this->comm().get(),&subvec2);
1349  LIBMESH_CHKERR(ierr);
1350  ierr = VecSetSizes(subvec2,is_local_size,PETSC_DECIDE);
1351  LIBMESH_CHKERR(ierr);
1352  ierr = VecSetFromOptions(subvec2);
1353  LIBMESH_CHKERR(ierr);
1354  ierr = MatMult(submat1,subvec1,subvec2);
1355  LIBMESH_CHKERR(ierr);
1356  ierr = VecAXPY(subrhs,1.0,subvec2);
1357  LIBMESH_CHKERR(ierr);
1358 
1359  ierr = LibMeshVecScatterDestroy(&scatter1);
1360  LIBMESH_CHKERR(ierr);
1361  ierr = LibMeshVecDestroy(&subvec1);
1362  LIBMESH_CHKERR(ierr);
1363  ierr = LibMeshMatDestroy(&submat1);
1364  LIBMESH_CHKERR(ierr);
1365  }
1366 
1367 #if PETSC_RELEASE_LESS_THAN(3,5,0)
1368  ierr = KSPSetOperators(_ksp, submat, subprecond,
1369  DIFFERENT_NONZERO_PATTERN);
1370 #else
1371  ierr = KSPSetOperators(_ksp, submat, subprecond);
1372 #endif
1373  LIBMESH_CHKERR(ierr);
1374 
1375  if (this->_preconditioner)
1376  {
1377  subprecond_matrix.reset(new PetscMatrix<Number>(subprecond, this->comm()));
1378  this->_preconditioner->set_matrix(*subprecond_matrix);
1379  this->_preconditioner->init();
1380  }
1381  }
1382  else
1383  {
1384 #if PETSC_RELEASE_LESS_THAN(3,5,0)
1385  ierr = KSPSetOperators(_ksp, mat, const_cast<PetscMatrix<T> *>(precond)->mat(),
1386  DIFFERENT_NONZERO_PATTERN);
1387 #else
1388  ierr = KSPSetOperators(_ksp, mat, const_cast<PetscMatrix<T> *>(precond)->mat());
1389 #endif
1390  LIBMESH_CHKERR(ierr);
1391 
1392  if (this->_preconditioner)
1393  {
1394  this->_preconditioner->set_matrix(const_cast<SparseMatrix<Number> &>(precond_matrix));
1395  this->_preconditioner->init();
1396  }
1397  }
1398 
1399  // Set the tolerances for the iterative solver. Use the user-supplied
1400  // tolerance for the relative residual & leave the others at default values.
1401  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
1402  PETSC_DEFAULT, max_its);
1403  LIBMESH_CHKERR(ierr);
1404 
1405  // Allow command line options to override anything set programmatically.
1406  ierr = KSPSetFromOptions(_ksp);
1407  LIBMESH_CHKERR(ierr);
1408 
1409  // Solve the linear system
1410  if (_restrict_solve_to_is != libmesh_nullptr)
1411  {
1412  ierr = KSPSolve (_ksp, subrhs, subsolution);
1413  LIBMESH_CHKERR(ierr);
1414  }
1415  else
1416  {
1417  ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
1418  LIBMESH_CHKERR(ierr);
1419  }
1420 
1421  // Get the number of iterations required for convergence
1422  ierr = KSPGetIterationNumber (_ksp, &its);
1423  LIBMESH_CHKERR(ierr);
1424 
1425  // Get the norm of the final residual to return to the user.
1426  ierr = KSPGetResidualNorm (_ksp, &final_resid);
1427  LIBMESH_CHKERR(ierr);
1428 
1429  if (_restrict_solve_to_is != libmesh_nullptr)
1430  {
1431  switch(_subset_solve_mode)
1432  {
1433  case SUBSET_ZERO:
1434  ierr = VecZeroEntries(solution->vec());
1435  LIBMESH_CHKERR(ierr);
1436  break;
1437 
1438  case SUBSET_COPY_RHS:
1439  ierr = VecCopy(rhs->vec(),solution->vec());
1440  LIBMESH_CHKERR(ierr);
1441  break;
1442 
1443  case SUBSET_DONT_TOUCH:
1444  /* Nothing to do here. */
1445  break;
1446 
1447  default:
1448  libmesh_error_msg("Invalid subset solve mode = " << _subset_solve_mode);
1449  }
1450  ierr = VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1451  LIBMESH_CHKERR(ierr);
1452  ierr = VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1453  LIBMESH_CHKERR(ierr);
1454 
1455  ierr = LibMeshVecScatterDestroy(&scatter);
1456  LIBMESH_CHKERR(ierr);
1457 
1458  if (this->_preconditioner)
1459  {
1460  // Before subprecond_matrix gets cleaned up, we should give
1461  // the _preconditioner a different matrix.
1462  this->_preconditioner->set_matrix(const_cast<SparseMatrix<Number> &>(precond_matrix));
1463  this->_preconditioner->init();
1464  }
1465 
1466  ierr = LibMeshVecDestroy(&subsolution);
1467  LIBMESH_CHKERR(ierr);
1468  ierr = LibMeshVecDestroy(&subrhs);
1469  LIBMESH_CHKERR(ierr);
1470  ierr = LibMeshMatDestroy(&submat);
1471  LIBMESH_CHKERR(ierr);
1472  ierr = LibMeshMatDestroy(&subprecond);
1473  LIBMESH_CHKERR(ierr);
1474  }
1475 
1476  // Destroy the matrix.
1477  ierr = LibMeshMatDestroy(&mat);
1478  LIBMESH_CHKERR(ierr);
1479 
1480  // return the # of its. and the final residual norm.
1481  return std::make_pair(its, final_resid);
1482 }
Leaves dofs outside the subset unchanged. This is fastest, but also most confusing because it abandon...
unsigned int size() const
Definition: parallel.h:725
Preconditioner< T > * _preconditioner
Set dofs outside the subset to zero.
void _create_complement_is(const NumericVector< T > &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)
virtual void init(const char *name=libmesh_nullptr) libmesh_override
Set dofs outside the subset to the value of the corresponding dofs of the right hand side...
static PetscErrorCode _petsc_shell_matrix_mult(Mat mat, Vec arg, Vec dest)
PetscErrorCode ierr
const Parallel::Communicator & comm() const
static PetscErrorCode _petsc_shell_matrix_get_diagonal(Mat mat, Vec dest)
template<typename T>
std::pair< unsigned int, Real > libMesh::LinearSolver< T >::solve ( SparseMatrix< T > &  matrix,
SparseMatrix< T > *  precond_matrix,
NumericVector< T > &  sol,
NumericVector< T > &  rhs,
const double  tol,
const unsigned int  n_iter 
)
inlineinherited

This function calls the solver "_solver_type" preconditioned with the "_preconditioner_type" preconditioner. The preconditioning matrix is used if it is provided, or the system matrix is used if precond_matrix is null

Definition at line 335 of file linear_solver.h.

References libMesh::LinearSolver< T >::solve().

341 {
342  if (pc_mat)
343  return this->solve(mat, *pc_mat, sol, rhs, tol, n_iter);
344  else
345  return this->solve(mat, sol, rhs, tol, n_iter);
346 }
virtual std::pair< unsigned int, Real > solve(SparseMatrix< T > &, NumericVector< T > &, NumericVector< T > &, const double, const unsigned int)=0
template<typename T>
std::pair< unsigned int, Real > libMesh::LinearSolver< T >::solve ( const ShellMatrix< T > &  matrix,
const SparseMatrix< T > *  precond_matrix,
NumericVector< T > &  sol,
NumericVector< T > &  rhs,
const double  tol,
const unsigned int  n_iter 
)
inlineinherited

This function solves a system whose matrix is a shell matrix, but an optional sparse matrix may be used as preconditioning matrix.

Definition at line 352 of file linear_solver.h.

References libMesh::LinearSolver< T >::solve().

358 {
359  if (pc_mat)
360  return this->solve(mat, *pc_mat, sol, rhs, tol, n_iter);
361  else
362  return this->solve(mat, sol, rhs, tol, n_iter);
363 }
virtual std::pair< unsigned int, Real > solve(SparseMatrix< T > &, NumericVector< T > &, NumericVector< T > &, const double, const unsigned int)=0
template<typename T>
SolverType libMesh::LinearSolver< T >::solver_type ( ) const
inlineinherited
Returns
The type of solver to use.

Definition at line 112 of file linear_solver.h.

112 { return _solver_type; }

Member Data Documentation

ReferenceCounter::Counts libMesh::ReferenceCounter::_counts
staticprotectedinherited
bool libMesh::ReferenceCounter::_enable_print_counter = true
staticprotectedinherited

Flag to control whether reference count information is printed when print_info is called.

Definition at line 143 of file reference_counter.h.

Referenced by libMesh::ReferenceCounter::disable_print_counter_info(), libMesh::ReferenceCounter::enable_print_counter_info(), and libMesh::ReferenceCounter::print_info().

template<typename T>
bool libMesh::LinearSolver< T >::_is_initialized
protectedinherited

Flag indicating if the data structures have been initialized.

Definition at line 277 of file linear_solver.h.

Referenced by libMesh::EigenSparseLinearSolver< T >::clear(), libMesh::EigenSparseLinearSolver< T >::init(), and libMesh::LinearSolver< Number >::initialized().

template<typename T >
KSP libMesh::PetscLinearSolver< T >::_ksp
private

Krylov subspace context

Definition at line 275 of file petsc_linear_solver.h.

Referenced by libMesh::PetscLinearSolver< T >::ksp().

Threads::spin_mutex libMesh::ReferenceCounter::_mutex
staticprotectedinherited

Mutual exclusion object to enable thread-safe reference counting.

Definition at line 137 of file reference_counter.h.

Threads::atomic< unsigned int > libMesh::ReferenceCounter::_n_objects
staticprotectedinherited

The number of objects. Print the reference count information when the number returns to 0.

Definition at line 132 of file reference_counter.h.

Referenced by libMesh::ReferenceCounter::n_objects(), libMesh::ReferenceCounter::ReferenceCounter(), and libMesh::ReferenceCounter::~ReferenceCounter().

template<typename T >
PC libMesh::PetscLinearSolver< T >::_pc
private

Preconditioner context

Definition at line 270 of file petsc_linear_solver.h.

Referenced by libMesh::PetscLinearSolver< T >::pc().

template<typename T>
Preconditioner<T>* libMesh::LinearSolver< T >::_preconditioner
protectedinherited

Holds the Preconditioner object to be used for the linear solves.

Definition at line 282 of file linear_solver.h.

template<typename T>
PreconditionerType libMesh::LinearSolver< T >::_preconditioner_type
protectedinherited
template<typename T >
IS libMesh::PetscLinearSolver< T >::_restrict_solve_to_is
private

PETSc index set containing the dofs on which to solve (NULL means solve on all dofs).

Definition at line 281 of file petsc_linear_solver.h.

Referenced by libMesh::PetscLinearSolver< T >::_create_complement_is(), and libMesh::PetscLinearSolver< T >::_restrict_solve_to_is_local_size().

template<typename T >
IS libMesh::PetscLinearSolver< T >::_restrict_solve_to_is_complement
private

PETSc index set, complement to _restrict_solve_to_is. This will be created on demand by the method _create_complement_is().

Definition at line 288 of file petsc_linear_solver.h.

Referenced by libMesh::PetscLinearSolver< T >::_create_complement_is().

template<typename T>
SolverConfiguration* libMesh::LinearSolver< T >::_solver_configuration
protectedinherited

Optionally store a SolverOptions object that can be used to set parameters like solver type, tolerances and iteration limits.

Definition at line 296 of file linear_solver.h.

Referenced by libMesh::EigenSparseLinearSolver< T >::solve().

template<typename T >
SubsetSolveMode libMesh::PetscLinearSolver< T >::_subset_solve_mode
private

If restrict-solve-to-subset mode is active, this member decides what happens with the dofs outside the subset.

Definition at line 306 of file petsc_linear_solver.h.

template<typename T>
bool libMesh::LinearSolver< T >::same_preconditioner
protectedinherited

Boolean flag to indicate whether we want to use an identical preconditioner to the previous solve. This can save substantial work in the cases where the system matrix is the same for successive solves.

Definition at line 290 of file linear_solver.h.

Referenced by libMesh::LinearSolver< T >::get_same_preconditioner().


The documentation for this class was generated from the following files: