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 UniquePtr< 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 110 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 309 of file petsc_linear_solver.h.

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

309  :
310  LinearSolver<T>(comm_in),
314 {
315  if (this->n_processors() == 1)
317  else
319 }
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 325 of file petsc_linear_solver.h.

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

326 {
327  this->clear ();
328 }
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 349 of file petsc_linear_solver.h.

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

356 {
358 #if PETSC_VERSION_LESS_THAN(3,0,0)
359  // No ISComplement in PETSc 2.3.3
360  libmesh_not_implemented();
361 #else
363  {
364  int ierr = ISComplement(_restrict_solve_to_is,
365  vec_in.first_local_index(),
366  vec_in.last_local_index(),
368  LIBMESH_CHKERR(ierr);
369  }
370 #endif
371 }
const class libmesh_nullptr_t libmesh_nullptr
libmesh_assert(j)
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 1736 of file petsc_linear_solver.C.

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

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

1737 {
1738  /* Get the matrix context. */
1739  PetscErrorCode ierr=0;
1740  void * ctx;
1741  ierr = MatShellGetContext(mat,&ctx);
1742 
1743  /* Get user shell matrix object. */
1744  const ShellMatrix<T> & shell_matrix = *static_cast<const ShellMatrix<T> *>(ctx);
1745  CHKERRABORT(shell_matrix.comm().get(), ierr);
1746 
1747  /* Make \p NumericVector instances around the vector. */
1748  PetscVector<T> dest_global(dest, shell_matrix.comm());
1749 
1750  /* Call the user function. */
1751  shell_matrix.get_diagonal(dest_global);
1752 
1753  return ierr;
1754 }
PetscErrorCode Vec Mat Mat void * ctx
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 1682 of file petsc_linear_solver.C.

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

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

1683 {
1684  /* Get the matrix context. */
1685  PetscErrorCode ierr=0;
1686  void * ctx;
1687  ierr = MatShellGetContext(mat,&ctx);
1688 
1689  /* Get user shell matrix object. */
1690  const ShellMatrix<T> & shell_matrix = *static_cast<const ShellMatrix<T> *>(ctx);
1691  CHKERRABORT(shell_matrix.comm().get(), ierr);
1692 
1693  /* Make \p NumericVector instances around the vectors. */
1694  PetscVector<T> arg_global(arg, shell_matrix.comm());
1695  PetscVector<T> dest_global(dest, shell_matrix.comm());
1696 
1697  /* Call the user function. */
1698  shell_matrix.vector_mult(dest_global,arg_global);
1699 
1700  return ierr;
1701 }
PetscErrorCode Vec Mat Mat void * ctx
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 1706 of file petsc_linear_solver.C.

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

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

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

Internal method that returns the local size of _restrict_solve_to_is.

Definition at line 334 of file petsc_linear_solver.h.

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

335 {
337 
338  PetscInt s;
339  int ierr = ISGetLocalSize(_restrict_solve_to_is, &s);
340  LIBMESH_CHKERR(ierr);
341 
342  return s;
343 }
libmesh_assert(j)
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 660 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().

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

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

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

103 {
104  if (this->_is_initialized)
105  libmesh_error_msg("Preconditioner must be attached before the solver is initialized!");
106 
108  _preconditioner = preconditioner;
109 }
Preconditioner< T > * _preconditioner
PreconditionerType _preconditioner_type
template<typename T >
UniquePtr< 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 42 of file linear_solver.C.

References libMesh::EIGEN_SOLVERS, libMesh::LASPACK_SOLVERS, libMesh::PETSC_SOLVERS, and libMesh::TRILINOS_SOLVERS.

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

44 {
45  // Build the appropriate solver
46  switch (solver_package)
47  {
48 
49 
50 #ifdef LIBMESH_HAVE_LASPACK
51  case LASPACK_SOLVERS:
52  return UniquePtr<LinearSolver<T> >(new LaspackLinearSolver<T>(comm));
53 #endif
54 
55 
56 #ifdef LIBMESH_HAVE_PETSC
57  case PETSC_SOLVERS:
58  return UniquePtr<LinearSolver<T> >(new PetscLinearSolver<T>(comm));
59 #endif
60 
61 
62 #ifdef LIBMESH_TRILINOS_HAVE_AZTECOO
63  case TRILINOS_SOLVERS:
64  return UniquePtr<LinearSolver<T> >(new AztecLinearSolver<T>(comm));
65 #endif
66 
67 
68 #ifdef LIBMESH_HAVE_EIGEN
69  case EIGEN_SOLVERS:
70  return UniquePtr<LinearSolver<T> >(new EigenSparseLinearSolver<T>(comm));
71 #endif
72 
73  default:
74  libmesh_error_msg("ERROR: Unrecognized solver package: " << solver_package);
75  }
76 
77  return UniquePtr<LinearSolver<T> >();
78 }
EIGEN_SOLVERS
Definition: libmesh.C:260
TRILINOS_SOLVERS
Definition: libmesh.C:258
LASPACK_SOLVERS
Definition: libmesh.C:262
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_monitor(), libMesh::__libmesh_petsc_diff_solver_residual(), libMesh::__libmesh_petsc_snes_jacobian(), 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::MeshTools::bounding_box(), 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::MeshTools::correct_node_proc_ids(), libMesh::MeshRefinement::create_parent_error_vector(), 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::PetscDiffSolver::init(), libMesh::TimeSolver::init(), libMesh::TopologyMap::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::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::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::MeshTools::processor_bounding_box(), libMesh::Nemesis_IO::read(), libMesh::XdrIO::read(), 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::DistributedMesh::renumber_dof_objects(), libMesh::MeshCommunication::send_coarse_ghosts(), libMesh::Partitioner::set_node_processor_ids(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::Partitioner::set_parent_processor_ids(), libMesh::LaplaceMeshSmoother::smooth(), libMesh::Parallel::Sort< KeyType, IdxType >::sort(), libMesh::MeshTools::subdomain_bounding_box(), 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::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::LaspackVector< T >::zero_clone(), libMesh::EigenSparseVector< 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 100 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

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

101 {
102  _enable_print_counter = true;
103  return;
104 }
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 1639 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().

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

References libMesh::err, and ierr.

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

1526 {
1527  PetscErrorCode ierr = 0;
1528  PetscInt its = 0;
1529 
1530  // Fill the residual history vector with the residual norms
1531  // Note that GetResidualHistory() does not copy any values, it
1532  // simply sets the pointer p. Note that for some Krylov subspace
1533  // methods, the number of residuals returned in the history
1534  // vector may be different from what you are expecting. For
1535  // example, TFQMR returns two residual values per iteration step.
1536  PetscReal * p;
1537  ierr = KSPGetResidualHistory(_ksp, &p, &its);
1538  LIBMESH_CHKERR(ierr);
1539 
1540  // Check no residual history
1541  if (its == 0)
1542  {
1543  libMesh::err << "No iterations have been performed, returning 0." << std::endl;
1544  return 0.;
1545  }
1546 
1547  // Otherwise, return the value pointed to by p.
1548  return *p;
1549 }
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 1492 of file petsc_linear_solver.C.

References ierr.

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

1493 {
1494  PetscErrorCode ierr = 0;
1495  PetscInt its = 0;
1496 
1497  // Fill the residual history vector with the residual norms
1498  // Note that GetResidualHistory() does not copy any values, it
1499  // simply sets the pointer p. Note that for some Krylov subspace
1500  // methods, the number of residuals returned in the history
1501  // vector may be different from what you are expecting. For
1502  // example, TFQMR returns two residual values per iteration step.
1503  PetscReal * p;
1504  ierr = KSPGetResidualHistory(_ksp, &p, &its);
1505  LIBMESH_CHKERR(ierr);
1506 
1507  // Check for early return
1508  if (its == 0) return;
1509 
1510  // Create space to store the result
1511  hist.resize(its);
1512 
1513  // Copy history into the vector provided by the user.
1514  for (PetscInt i=0; i<its; ++i)
1515  {
1516  hist[i] = *p;
1517  p++;
1518  }
1519 }
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 324 of file linear_solver.h.

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

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

325 {
326  return same_preconditioner;
327 }
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 160 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().

161 {
162  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
163  std::pair<unsigned int, unsigned int> & p = _counts[name];
164 
165  p.first++;
166 }
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 173 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().

174 {
175  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
176  std::pair<unsigned int, unsigned int> & p = _counts[name];
177 
178  p.second++;
179 }
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::pc, and libMesh::petsc_auto_fieldsplit().

344 {
345  petsc_auto_fieldsplit(this->pc(), sys);
346 }
void petsc_auto_fieldsplit(PC my_pc, const System &sys)
ImplicitSystem & 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::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::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::MeshTools::processor_bounding_box(), libMesh::NameBasedIO::read(), libMesh::Nemesis_IO::read(), libMesh::XdrIO::read(), 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::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:679
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 PCDestroy()!

Definition at line 211 of file petsc_linear_solver.h.

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

211 { 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 82 of file linear_solver.C.

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

83 {
84  if(_preconditioner)
85  return _preconditioner->type();
86 
87  return _preconditioner_type;
88 }
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 154 of file linear_solver.C.

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

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

155 {
157  libMesh::out << "Linear solver convergence/divergence reason: " << Utility::enum_to_string(reason) << std::endl;
158 }
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 }
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::ExodusII_IO_Helper::read_elem_num_map(), 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::MeshCommunication::send_coarse_ghosts(), libMesh::Partitioner::set_node_processor_ids(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::LaplaceMeshSmoother::smooth(), libMesh::MeshTools::total_weight(), 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::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(), and libMesh::ExodusII_IO_Helper::write_timestep().

100  { return cast_int<processor_id_type>(_communicator.rank()); }
const Parallel::Communicator & _communicator
unsigned int rank() const
Definition: parallel.h:677
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:679
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 113 of file linear_solver.C.

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

114 {
115  same_preconditioner = reuse_flag;
116 }
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 1555 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().

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

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

93 {
94  if(_preconditioner)
95  _preconditioner->set_type(pct);
96  else
98 }
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 161 of file linear_solver.C.

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

162 {
163  _solver_configuration = &solver_configuration;
164 }
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 that 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(). Note: 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  UniquePtr<PetscMatrix<Number> > subprecond_matrix;
424 
425  // Set operators. Also restrict rhs and solution vector to
426  // subdomain if neccessary.
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 #if PETSC_VERSION_LESS_THAN(3,1,0)
459  ierr = MatGetSubMatrix(matrix->mat(),
461  PETSC_DECIDE,MAT_INITIAL_MATRIX,&submat);
462  LIBMESH_CHKERR(ierr);
463  ierr = MatGetSubMatrix(precond->mat(),
465  PETSC_DECIDE,MAT_INITIAL_MATRIX,&subprecond);
466  LIBMESH_CHKERR(ierr);
467 #else
468  ierr = MatGetSubMatrix(matrix->mat(),
470  MAT_INITIAL_MATRIX,&submat);
471  LIBMESH_CHKERR(ierr);
472  ierr = MatGetSubMatrix(precond->mat(),
474  MAT_INITIAL_MATRIX,&subprecond);
475  LIBMESH_CHKERR(ierr);
476 #endif
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 #if PETSC_VERSION_LESS_THAN(3,1,0)
511  ierr = MatGetSubMatrix(matrix->mat(),
513  PETSC_DECIDE,MAT_INITIAL_MATRIX,&submat1);
514  LIBMESH_CHKERR(ierr);
515 #else
516  ierr = MatGetSubMatrix(matrix->mat(),
518  MAT_INITIAL_MATRIX,&submat1);
519  LIBMESH_CHKERR(ierr);
520 #endif
521 
522  ierr = MatMultAdd(submat1,subvec1,subrhs,subrhs);
523  LIBMESH_CHKERR(ierr);
524 
525  ierr = LibMeshVecScatterDestroy(&scatter1);
526  LIBMESH_CHKERR(ierr);
527  ierr = LibMeshVecDestroy(&subvec1);
528  LIBMESH_CHKERR(ierr);
529  ierr = LibMeshMatDestroy(&submat1);
530  LIBMESH_CHKERR(ierr);
531  }
532 #if PETSC_RELEASE_LESS_THAN(3,5,0)
533  ierr = KSPSetOperators(_ksp, submat, subprecond,
534  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
535 #else
536  ierr = KSPSetOperators(_ksp, submat, subprecond);
537 
538  PetscBool ksp_reuse_preconditioner = this->same_preconditioner ? PETSC_TRUE : PETSC_FALSE;
539  ierr = KSPSetReusePreconditioner(_ksp, ksp_reuse_preconditioner);
540 #endif
541  LIBMESH_CHKERR(ierr);
542 
543  if (this->_preconditioner)
544  {
545  subprecond_matrix.reset(new PetscMatrix<Number>(subprecond, this->comm()));
546  this->_preconditioner->set_matrix(*subprecond_matrix);
547  this->_preconditioner->init();
548  }
549  }
550  else
551  {
552 #if PETSC_RELEASE_LESS_THAN(3,5,0)
553  ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
554  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
555 #else
556  ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat());
557 
558  PetscBool ksp_reuse_preconditioner = this->same_preconditioner ? PETSC_TRUE : PETSC_FALSE;
559  ierr = KSPSetReusePreconditioner(_ksp, ksp_reuse_preconditioner);
560 #endif
561  LIBMESH_CHKERR(ierr);
562 
563  if(this->_preconditioner)
564  {
565  this->_preconditioner->set_matrix(matrix_in);
566  this->_preconditioner->init();
567  }
568  }
569 
570  // Set the tolerances for the iterative solver. Use the user-supplied
571  // tolerance for the relative residual & leave the others at default values.
572  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
573  PETSC_DEFAULT, max_its);
574  LIBMESH_CHKERR(ierr);
575 
576  // Allow command line options to override anything set programmatically.
577  ierr = KSPSetFromOptions(_ksp);
578  LIBMESH_CHKERR(ierr);
579 
580  // If the SolverConfiguration object is provided, use it to override
581  // solver options.
582  if(this->_solver_configuration)
583  {
585  }
586 
587  // Solve the linear system
588  if (_restrict_solve_to_is != libmesh_nullptr)
589  {
590  ierr = KSPSolve (_ksp, subrhs, subsolution);
591  LIBMESH_CHKERR(ierr);
592  }
593  else
594  {
595  ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
596  LIBMESH_CHKERR(ierr);
597  }
598 
599  // Get the number of iterations required for convergence
600  ierr = KSPGetIterationNumber (_ksp, &its);
601  LIBMESH_CHKERR(ierr);
602 
603  // Get the norm of the final residual to return to the user.
604  ierr = KSPGetResidualNorm (_ksp, &final_resid);
605  LIBMESH_CHKERR(ierr);
606 
607  if (_restrict_solve_to_is != libmesh_nullptr)
608  {
609  switch(_subset_solve_mode)
610  {
611  case SUBSET_ZERO:
612  ierr = VecZeroEntries(solution->vec());
613  LIBMESH_CHKERR(ierr);
614  break;
615 
616  case SUBSET_COPY_RHS:
617  ierr = VecCopy(rhs->vec(),solution->vec());
618  LIBMESH_CHKERR(ierr);
619  break;
620 
621  case SUBSET_DONT_TOUCH:
622  /* Nothing to do here. */
623  break;
624 
625  default:
626  libmesh_error_msg("Invalid subset solve mode = " << _subset_solve_mode);
627  }
628  ierr = VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
629  LIBMESH_CHKERR(ierr);
630  ierr = VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
631  LIBMESH_CHKERR(ierr);
632 
633  ierr = LibMeshVecScatterDestroy(&scatter);
634  LIBMESH_CHKERR(ierr);
635 
636  if (this->_preconditioner)
637  {
638  // Before subprecond_matrix gets cleaned up, we should give
639  // the _preconditioner a different matrix.
640  this->_preconditioner->set_matrix(matrix_in);
641  this->_preconditioner->init();
642  }
643 
644  ierr = LibMeshVecDestroy(&subsolution);
645  LIBMESH_CHKERR(ierr);
646  ierr = LibMeshVecDestroy(&subrhs);
647  LIBMESH_CHKERR(ierr);
648  ierr = LibMeshMatDestroy(&submat);
649  LIBMESH_CHKERR(ierr);
650  ierr = LibMeshMatDestroy(&subprecond);
651  LIBMESH_CHKERR(ierr);
652  }
653 
654  // return the # of its. and the final residual norm.
655  return std::make_pair(its, final_resid);
656 }
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 923 of file petsc_linear_solver.C.

References libMesh::PetscVector< T >::close(), ierr, libMesh::TriangleWrapper::init(), libMesh::libmesh_assert(), 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().

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

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

338 {
339  if (pc_mat)
340  return this->solve(mat, *pc_mat, sol, rhs, tol, n_iter);
341  else
342  return this->solve(mat, sol, rhs, tol, n_iter);
343 }
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::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 1192 of file petsc_linear_solver.C.

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

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

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

355 {
356  if (pc_mat)
357  return this->solve(mat, *pc_mat, sol, rhs, tol, n_iter);
358  else
359  return this->solve(mat, sol, rhs, tol, n_iter);
360 }
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 134 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 274 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 270 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 128 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 123 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 265 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 279 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 276 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 283 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 293 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 302 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 287 of file linear_solver.h.

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


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