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

Static Public Member Functions

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

Protected Types

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

Protected Member Functions

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

Protected Attributes

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

Static Protected Attributes

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

Private Member Functions

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

Static Private Member Functions

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

Private Attributes

PC _pc
 
KSP _ksp
 
IS _restrict_solve_to_is
 
IS _restrict_solve_to_is_complement
 
SubsetSolveMode _subset_solve_mode
 

Detailed Description

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

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

Author
Benjamin Kirk
Date
2002-2007

Definition at line 101 of file petsc_linear_solver.h.

Member Typedef Documentation

◆ Counts

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 117 of file reference_counter.h.

Constructor & Destructor Documentation

◆ PetscLinearSolver()

◆ ~PetscLinearSolver()

template<typename T >
libMesh::PetscLinearSolver< T >::~PetscLinearSolver ( )
inline

Destructor.

Definition at line 331 of file petsc_linear_solver.h.

332 {
333  this->clear ();
334 }
virtual void clear() override

Member Function Documentation

◆ _create_complement_is()

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 355 of file petsc_linear_solver.h.

References ierr.

362 {
363  libmesh_assert(_restrict_solve_to_is);
364 #if PETSC_VERSION_LESS_THAN(3,0,0)
365  // No ISComplement in PETSc 2.3.3
366  libmesh_not_implemented();
367 #else
369  {
370  int ierr = ISComplement(_restrict_solve_to_is,
371  vec_in.first_local_index(),
372  vec_in.last_local_index(),
374  LIBMESH_CHKERR(ierr);
375  }
376 #endif
377 }
PetscErrorCode ierr

◆ _petsc_shell_matrix_get_diagonal()

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 1781 of file petsc_linear_solver.C.

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

1782 {
1783  /* Get the matrix context. */
1784  PetscErrorCode ierr=0;
1785  void * ctx;
1786  ierr = MatShellGetContext(mat,&ctx);
1787 
1788  /* Get user shell matrix object. */
1789  const ShellMatrix<T> & shell_matrix = *static_cast<const ShellMatrix<T> *>(ctx);
1790  CHKERRABORT(shell_matrix.comm().get(), ierr);
1791 
1792  /* Make \p NumericVector instances around the vector. */
1793  PetscVector<T> dest_global(dest, shell_matrix.comm());
1794 
1795  /* Call the user function. */
1796  shell_matrix.get_diagonal(dest_global);
1797 
1798  return ierr;
1799 }
PetscErrorCode ierr

◆ _petsc_shell_matrix_mult()

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 1727 of file petsc_linear_solver.C.

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

1728 {
1729  /* Get the matrix context. */
1730  PetscErrorCode ierr=0;
1731  void * ctx;
1732  ierr = MatShellGetContext(mat,&ctx);
1733 
1734  /* Get user shell matrix object. */
1735  const ShellMatrix<T> & shell_matrix = *static_cast<const ShellMatrix<T> *>(ctx);
1736  CHKERRABORT(shell_matrix.comm().get(), ierr);
1737 
1738  /* Make \p NumericVector instances around the vectors. */
1739  PetscVector<T> arg_global(arg, shell_matrix.comm());
1740  PetscVector<T> dest_global(dest, shell_matrix.comm());
1741 
1742  /* Call the user function. */
1743  shell_matrix.vector_mult(dest_global,arg_global);
1744 
1745  return ierr;
1746 }
PetscErrorCode ierr

◆ _petsc_shell_matrix_mult_add()

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 1751 of file petsc_linear_solver.C.

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

1752 {
1753  /* Get the matrix context. */
1754  PetscErrorCode ierr=0;
1755  void * ctx;
1756  ierr = MatShellGetContext(mat,&ctx);
1757 
1758  /* Get user shell matrix object. */
1759  const ShellMatrix<T> & shell_matrix = *static_cast<const ShellMatrix<T> *>(ctx);
1760  CHKERRABORT(shell_matrix.comm().get(), ierr);
1761 
1762  /* Make \p NumericVector instances around the vectors. */
1763  PetscVector<T> arg_global(arg, shell_matrix.comm());
1764  PetscVector<T> dest_global(dest, shell_matrix.comm());
1765  PetscVector<T> add_global(add, shell_matrix.comm());
1766 
1767  if (add!=arg)
1768  {
1769  arg_global = add_global;
1770  }
1771 
1772  /* Call the user function. */
1773  shell_matrix.vector_mult_add(dest_global,arg_global);
1774 
1775  return ierr;
1776 }
PetscErrorCode ierr

◆ _restrict_solve_to_is_local_size()

template<typename T >
PetscInt libMesh::PetscLinearSolver< T >::_restrict_solve_to_is_local_size ( ) const
inlineprivate
Returns
The local size of _restrict_solve_to_is.

Definition at line 340 of file petsc_linear_solver.h.

References ierr.

341 {
342  libmesh_assert(_restrict_solve_to_is);
343 
344  PetscInt s;
345  int ierr = ISGetLocalSize(_restrict_solve_to_is, &s);
346  LIBMESH_CHKERR(ierr);
347 
348  return s;
349 }
PetscErrorCode ierr

◆ adjoint_solve()

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 
)
overridevirtual

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 708 of file petsc_linear_solver.C.

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

713 {
714  LOG_SCOPE("solve()", "PetscLinearSolver");
715 
716  // Make sure the data passed in are really of Petsc types
717  PetscMatrix<T> * matrix = cast_ptr<PetscMatrix<T> *>(&matrix_in);
718  // Note that the matrix and precond matrix are the same
719  PetscMatrix<T> * precond = cast_ptr<PetscMatrix<T> *>(&matrix_in);
720  PetscVector<T> * solution = cast_ptr<PetscVector<T> *>(&solution_in);
721  PetscVector<T> * rhs = cast_ptr<PetscVector<T> *>(&rhs_in);
722 
723  this->init (matrix);
724 
725  PetscErrorCode ierr=0;
726  PetscInt its=0, max_its = static_cast<PetscInt>(m_its);
727  PetscReal final_resid=0.;
728 
729  // Close the matrices and vectors in case this wasn't already done.
730  matrix->close ();
731  precond->close ();
732  solution->close ();
733  rhs->close ();
734 
735  Mat submat = nullptr;
736  Mat subprecond = nullptr;
737  Vec subrhs = nullptr;
738  Vec subsolution = nullptr;
739  VecScatter scatter = nullptr;
740  std::unique_ptr<PetscMatrix<Number>> subprecond_matrix;
741 
742  // Set operators. Also restrict rhs and solution vector to
743  // subdomain if necessary.
744  if (_restrict_solve_to_is != nullptr)
745  {
746  PetscInt is_local_size = this->_restrict_solve_to_is_local_size();
747 
748  ierr = VecCreate(this->comm().get(),&subrhs);
749  LIBMESH_CHKERR(ierr);
750  ierr = VecSetSizes(subrhs,is_local_size,PETSC_DECIDE);
751  LIBMESH_CHKERR(ierr);
752  ierr = VecSetFromOptions(subrhs);
753  LIBMESH_CHKERR(ierr);
754 
755  ierr = VecCreate(this->comm().get(),&subsolution);
756  LIBMESH_CHKERR(ierr);
757  ierr = VecSetSizes(subsolution,is_local_size,PETSC_DECIDE);
758  LIBMESH_CHKERR(ierr);
759  ierr = VecSetFromOptions(subsolution);
760  LIBMESH_CHKERR(ierr);
761 
762  ierr = LibMeshVecScatterCreate(rhs->vec(),_restrict_solve_to_is, subrhs, nullptr, &scatter);
763  LIBMESH_CHKERR(ierr);
764 
765  ierr = VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
766  LIBMESH_CHKERR(ierr);
767  ierr = VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
768  LIBMESH_CHKERR(ierr);
769 
770  ierr = VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
771  LIBMESH_CHKERR(ierr);
772  ierr = VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
773  LIBMESH_CHKERR(ierr);
774 
775  ierr = LibMeshCreateSubMatrix(matrix->mat(),
778 #if PETSC_VERSION_LESS_THAN(3,1,0)
779  PETSC_DECIDE,
780 #endif
781  MAT_INITIAL_MATRIX,
782  &submat);
783  LIBMESH_CHKERR(ierr);
784 
785  ierr = LibMeshCreateSubMatrix(precond->mat(),
788 #if PETSC_VERSION_LESS_THAN(3,1,0)
789  PETSC_DECIDE,
790 #endif
791  MAT_INITIAL_MATRIX,
792  &subprecond);
793  LIBMESH_CHKERR(ierr);
794 
795  /* Since removing columns of the matrix changes the equation
796  system, we will now change the right hand side to compensate
797  for this. Note that this is not necessary if \p SUBSET_ZERO
798  has been selected. */
800  {
801  _create_complement_is(rhs_in);
802  PetscInt is_complement_local_size =
803  cast_int<PetscInt>(rhs_in.local_size()-is_local_size);
804 
805  Vec subvec1 = nullptr;
806  Mat submat1 = nullptr;
807  VecScatter scatter1 = nullptr;
808 
809  ierr = VecCreate(this->comm().get(),&subvec1);
810  LIBMESH_CHKERR(ierr);
811  ierr = VecSetSizes(subvec1,is_complement_local_size,PETSC_DECIDE);
812  LIBMESH_CHKERR(ierr);
813  ierr = VecSetFromOptions(subvec1);
814  LIBMESH_CHKERR(ierr);
815 
816  ierr = LibMeshVecScatterCreate(rhs->vec(),_restrict_solve_to_is_complement, subvec1, nullptr, &scatter1);
817  LIBMESH_CHKERR(ierr);
818 
819  ierr = VecScatterBegin(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
820  LIBMESH_CHKERR(ierr);
821  ierr = VecScatterEnd(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
822  LIBMESH_CHKERR(ierr);
823 
824  ierr = VecScale(subvec1,-1.0);
825  LIBMESH_CHKERR(ierr);
826 
827  ierr = LibMeshCreateSubMatrix(matrix->mat(),
830 #if PETSC_VERSION_LESS_THAN(3,1,0)
831  PETSC_DECIDE,
832 #endif
833  MAT_INITIAL_MATRIX,
834  &submat1);
835  LIBMESH_CHKERR(ierr);
836 
837  ierr = MatMultAdd(submat1,subvec1,subrhs,subrhs);
838  LIBMESH_CHKERR(ierr);
839 
840  ierr = LibMeshVecScatterDestroy(&scatter1);
841  LIBMESH_CHKERR(ierr);
842  ierr = LibMeshVecDestroy(&subvec1);
843  LIBMESH_CHKERR(ierr);
844  ierr = LibMeshMatDestroy(&submat1);
845  LIBMESH_CHKERR(ierr);
846  }
847 #if PETSC_RELEASE_LESS_THAN(3,5,0)
848  ierr = KSPSetOperators(_ksp, submat, subprecond,
849  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
850 #else
851  ierr = KSPSetOperators(_ksp, submat, subprecond);
852 
853  PetscBool ksp_reuse_preconditioner = this->same_preconditioner ? PETSC_TRUE : PETSC_FALSE;
854  ierr = KSPSetReusePreconditioner(_ksp, ksp_reuse_preconditioner);
855 #endif
856  LIBMESH_CHKERR(ierr);
857 
858  if (this->_preconditioner)
859  {
860  subprecond_matrix.reset(new PetscMatrix<Number>(subprecond, this->comm()));
861  this->_preconditioner->set_matrix(*subprecond_matrix);
862  this->_preconditioner->init();
863  }
864  }
865  else
866  {
867 #if PETSC_RELEASE_LESS_THAN(3,5,0)
868  ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
869  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
870 #else
871  ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat());
872 
873  PetscBool ksp_reuse_preconditioner = this->same_preconditioner ? PETSC_TRUE : PETSC_FALSE;
874  ierr = KSPSetReusePreconditioner(_ksp, ksp_reuse_preconditioner);
875 #endif
876  LIBMESH_CHKERR(ierr);
877 
878  if (this->_preconditioner)
879  {
880  this->_preconditioner->set_matrix(matrix_in);
881  this->_preconditioner->init();
882  }
883  }
884 
885  // Set the tolerances for the iterative solver. Use the user-supplied
886  // tolerance for the relative residual & leave the others at default values.
887  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
888  PETSC_DEFAULT, max_its);
889  LIBMESH_CHKERR(ierr);
890 
891  // Allow command line options to override anything set programmatically.
892  ierr = KSPSetFromOptions(_ksp);
893  LIBMESH_CHKERR(ierr);
894 
895  // Solve the linear system
896  if (_restrict_solve_to_is != nullptr)
897  {
898  ierr = KSPSolveTranspose (_ksp, subrhs, subsolution);
899  LIBMESH_CHKERR(ierr);
900  }
901  else
902  {
903  ierr = KSPSolveTranspose (_ksp, rhs->vec(), solution->vec());
904  LIBMESH_CHKERR(ierr);
905  }
906 
907  // Get the number of iterations required for convergence
908  ierr = KSPGetIterationNumber (_ksp, &its);
909  LIBMESH_CHKERR(ierr);
910 
911  // Get the norm of the final residual to return to the user.
912  ierr = KSPGetResidualNorm (_ksp, &final_resid);
913  LIBMESH_CHKERR(ierr);
914 
915  if (_restrict_solve_to_is != nullptr)
916  {
917  switch(_subset_solve_mode)
918  {
919  case SUBSET_ZERO:
920  ierr = VecZeroEntries(solution->vec());
921  LIBMESH_CHKERR(ierr);
922  break;
923 
924  case SUBSET_COPY_RHS:
925  ierr = VecCopy(rhs->vec(),solution->vec());
926  LIBMESH_CHKERR(ierr);
927  break;
928 
929  case SUBSET_DONT_TOUCH:
930  /* Nothing to do here. */
931  break;
932 
933  default:
934  libmesh_error_msg("Invalid subset solve mode = " << _subset_solve_mode);
935  }
936  ierr = VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
937  LIBMESH_CHKERR(ierr);
938  ierr = VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
939  LIBMESH_CHKERR(ierr);
940 
941  ierr = LibMeshVecScatterDestroy(&scatter);
942  LIBMESH_CHKERR(ierr);
943 
944  if (this->_preconditioner)
945  {
946  // Before subprecond_matrix gets cleaned up, we should give
947  // the _preconditioner a different matrix.
948  this->_preconditioner->set_matrix(matrix_in);
949  this->_preconditioner->init();
950  }
951 
952  ierr = LibMeshVecDestroy(&subsolution);
953  LIBMESH_CHKERR(ierr);
954  ierr = LibMeshVecDestroy(&subrhs);
955  LIBMESH_CHKERR(ierr);
956  ierr = LibMeshMatDestroy(&submat);
957  LIBMESH_CHKERR(ierr);
958  ierr = LibMeshMatDestroy(&subprecond);
959  LIBMESH_CHKERR(ierr);
960  }
961 
962  // return the # of its. and the final residual norm.
963  return std::make_pair(its, final_resid);
964 }
Preconditioner< T > * _preconditioner
virtual void init(const char *name=nullptr) override
void _create_complement_is(const NumericVector< T > &vec_in)
const Parallel::Communicator & comm() const
PetscTruth PetscBool
Definition: petsc_macro.h:67
PetscErrorCode ierr
PetscInt _restrict_solve_to_is_local_size() const

◆ attach_preconditioner()

template<typename T>
void libMesh::LinearSolver< T >::attach_preconditioner ( Preconditioner< T > *  preconditioner)
inherited

Attaches a Preconditioner object to be used

Definition at line 118 of file linear_solver.C.

119 {
120  if (this->_is_initialized)
121  libmesh_error_msg("Preconditioner must be attached before the solver is initialized!");
122 
124  _preconditioner = preconditioner;
125 }
Preconditioner< T > * _preconditioner
PreconditionerType _preconditioner_type

◆ build()

template<typename T >
std::unique_ptr< LinearSolver< T > > libMesh::LinearSolver< T >::build ( const libMesh::Parallel::Communicator comm_in,
const SolverPackage  solver_package = libMesh::default_solver_package() 
)
staticinherited

Builds a LinearSolver using the linear solver package specified by solver_package

Definition at line 57 of file linear_solver.C.

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

59 {
60  // Avoid unused parameter warnings when no solver packages are enabled.
62 
63  // Build the appropriate solver
64  switch (solver_package)
65  {
66 #ifdef LIBMESH_HAVE_LASPACK
67  case LASPACK_SOLVERS:
68  return libmesh_make_unique<LaspackLinearSolver<T>>(comm);
69 #endif
70 
71 
72 #ifdef LIBMESH_HAVE_PETSC
73  case PETSC_SOLVERS:
74  return libmesh_make_unique<PetscLinearSolver<T>>(comm);
75 #endif
76 
77 
78 #ifdef LIBMESH_TRILINOS_HAVE_AZTECOO
79  case TRILINOS_SOLVERS:
80  return libmesh_make_unique<AztecLinearSolver<T>>(comm);
81 #endif
82 
83 
84 #ifdef LIBMESH_HAVE_EIGEN
85  case EIGEN_SOLVERS:
86  return libmesh_make_unique<EigenSparseLinearSolver<T>>(comm);
87 #endif
88 
89  default:
90  libmesh_error_msg("ERROR: Unrecognized solver package: " << solver_package);
91  }
92 
93  return std::unique_ptr<LinearSolver<T>>();
94 }
EIGEN_SOLVERS
Definition: libmesh.C:246
TRILINOS_SOLVERS
Definition: libmesh.C:244
const Parallel::Communicator & comm() const
void libmesh_ignore(const Args &...)
LASPACK_SOLVERS
Definition: libmesh.C:248

◆ clear()

template<typename T >
void libMesh::PetscLinearSolver< T >::clear ( )
overridevirtual

Release all memory and clear data structures.

Reimplemented from libMesh::LinearSolver< T >.

Definition at line 149 of file petsc_linear_solver.C.

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

150 {
151  if (this->initialized())
152  {
153  /* If we were restricted to some subset, this restriction must
154  be removed and the subset index set destroyed. */
155  if (_restrict_solve_to_is != nullptr)
156  {
157  PetscErrorCode ierr = LibMeshISDestroy(&_restrict_solve_to_is);
158  LIBMESH_CHKERR(ierr);
159  _restrict_solve_to_is = nullptr;
160  }
161 
162  if (_restrict_solve_to_is_complement != nullptr)
163  {
164  PetscErrorCode ierr = LibMeshISDestroy(&_restrict_solve_to_is_complement);
165  LIBMESH_CHKERR(ierr);
167  }
168 
169  this->_is_initialized = false;
170 
171  PetscErrorCode ierr=0;
172 
173  ierr = LibMeshKSPDestroy(&_ksp);
174  LIBMESH_CHKERR(ierr);
175 
176  // Mimic PETSc default solver and preconditioner
177  this->_solver_type = GMRES;
178 
179  if (!this->_preconditioner)
180  {
181  if (this->n_processors() == 1)
183  else
185  }
186  }
187 }
Preconditioner< T > * _preconditioner
processor_id_type n_processors() const
bool initialized() const
Definition: linear_solver.h:96
PetscErrorCode ierr
PreconditionerType _preconditioner_type

◆ comm()

const Parallel::Communicator& libMesh::ParallelObject::comm ( ) const
inlineinherited
Returns
A reference to the Parallel::Communicator object used by this mesh.

Definition at line 89 of file parallel_object.h.

References libMesh::ParallelObject::_communicator.

Referenced by libMesh::__libmesh_petsc_diff_solver_jacobian(), libMesh::__libmesh_petsc_diff_solver_monitor(), libMesh::__libmesh_petsc_diff_solver_residual(), libMesh::__libmesh_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::UniformRefinementEstimator::_estimate_error(), libMesh::BoundaryInfo::_find_id_maps(), libMesh::SlepcEigenSolver< T >::_petsc_shell_matrix_get_diagonal(), libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_get_diagonal(), libMesh::SlepcEigenSolver< T >::_petsc_shell_matrix_mult(), libMesh::PetscLinearSolver< 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::PetscDMWrapper::add_dofs_helper(), libMesh::PetscDMWrapper::add_dofs_to_section(), libMesh::ImplicitSystem::add_matrix(), libMesh::System::add_vector(), libMesh::UnstructuredMesh::all_second_order(), libMesh::MeshTools::Modification::all_tri(), libMesh::LaplaceMeshSmoother::allgather_graph(), libMesh::FEMSystem::assemble_qoi(), libMesh::MeshCommunication::assign_global_indices(), libMesh::DofMap::attach_matrix(), libMesh::MeshTools::Generation::build_extrusion(), libMesh::BoundaryInfo::build_node_list_from_side_list(), libMesh::EquationSystems::build_parallel_elemental_solution_vector(), libMesh::EquationSystems::build_parallel_solution_vector(), libMesh::PetscDMWrapper::build_section(), libMesh::PetscDMWrapper::build_sf(), libMesh::MeshBase::cache_elem_dims(), libMesh::System::calculate_norm(), libMesh::DofMap::check_dirichlet_bcid_consistency(), libMesh::PetscDMWrapper::check_section_n_dofs(), libMesh::Nemesis_IO_Helper::compute_num_global_elem_blocks(), libMesh::Nemesis_IO_Helper::compute_num_global_nodesets(), libMesh::Nemesis_IO_Helper::compute_num_global_sidesets(), libMesh::Problem_Interface::computeF(), libMesh::Problem_Interface::computeJacobian(), libMesh::Problem_Interface::computePreconditioner(), libMesh::ExodusII_IO::copy_elemental_solution(), libMesh::MeshTools::correct_node_proc_ids(), libMesh::MeshTools::create_bounding_box(), libMesh::MeshTools::create_nodal_bounding_box(), libMesh::MeshRefinement::create_parent_error_vector(), libMesh::MeshTools::create_processor_bounding_box(), libMesh::MeshTools::create_subdomain_bounding_box(), libMesh::MeshCommunication::delete_remote_elements(), libMesh::DofMap::distribute_dofs(), DMlibMeshFunction(), DMlibMeshJacobian(), DMlibMeshSetSystem_libMesh(), DMVariableBounds_libMesh(), libMesh::MeshRefinement::eliminate_unrefined_patches(), libMesh::EpetraVector< T >::EpetraVector(), libMesh::WeightedPatchRecoveryErrorEstimator::estimate_error(), libMesh::PatchRecoveryErrorEstimator::estimate_error(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::ExactErrorEstimator::estimate_error(), libMesh::MeshRefinement::flag_elements_by_elem_fraction(), libMesh::MeshRefinement::flag_elements_by_error_fraction(), libMesh::MeshRefinement::flag_elements_by_nelem_target(), libMesh::CondensedEigenSystem::get_eigenpair(), libMesh::DofMap::get_info(), libMesh::ImplicitSystem::get_linear_solver(), libMesh::LocationMap< T >::init(), libMesh::TimeSolver::init(), libMesh::SystemSubsetBySubdomain::init(), libMesh::PetscDMWrapper::init_and_attach_petscdm(), libMesh::EigenSystem::init_data(), libMesh::EigenSystem::init_matrices(), libMesh::OptimizationSystem::initialize_equality_constraints_storage(), libMesh::OptimizationSystem::initialize_inequality_constraints_storage(), libMesh::MeshTools::libmesh_assert_consistent_distributed(), libMesh::MeshTools::libmesh_assert_consistent_distributed_nodes(), libMesh::MeshTools::libmesh_assert_contiguous_dof_ids(), libMesh::MeshTools::libmesh_assert_parallel_consistent_new_node_procids(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Elem >(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Node >(), libMesh::MeshTools::libmesh_assert_topology_consistent_procids< Node >(), libMesh::MeshTools::libmesh_assert_valid_boundary_ids(), libMesh::MeshTools::libmesh_assert_valid_dof_ids(), libMesh::MeshTools::libmesh_assert_valid_neighbors(), libMesh::DistributedMesh::libmesh_assert_valid_parallel_flags(), libMesh::DistributedMesh::libmesh_assert_valid_parallel_object_ids(), libMesh::DistributedMesh::libmesh_assert_valid_parallel_p_levels(), libMesh::MeshTools::libmesh_assert_valid_refinement_flags(), libMesh::MeshTools::libmesh_assert_valid_unique_ids(), libMesh::libmesh_petsc_snes_fd_residual(), libMesh::libmesh_petsc_snes_jacobian(), libMesh::libmesh_petsc_snes_mffd_residual(), libMesh::libmesh_petsc_snes_postcheck(), libMesh::libmesh_petsc_snes_residual(), libMesh::libmesh_petsc_snes_residual_helper(), libMesh::MeshRefinement::limit_level_mismatch_at_edge(), libMesh::MeshRefinement::limit_level_mismatch_at_node(), libMesh::MeshRefinement::limit_overrefined_boundary(), libMesh::MeshRefinement::limit_underrefined_boundary(), libMesh::MeshRefinement::make_coarsening_compatible(), libMesh::MeshCommunication::make_elems_parallel_consistent(), libMesh::MeshRefinement::make_flags_parallel_consistent(), libMesh::MeshCommunication::make_new_node_proc_ids_parallel_consistent(), libMesh::MeshCommunication::make_new_nodes_parallel_consistent(), libMesh::MeshCommunication::make_node_ids_parallel_consistent(), libMesh::MeshCommunication::make_node_proc_ids_parallel_consistent(), libMesh::MeshCommunication::make_node_unique_ids_parallel_consistent(), libMesh::MeshCommunication::make_nodes_parallel_consistent(), libMesh::MeshCommunication::make_p_levels_parallel_consistent(), libMesh::MeshRefinement::make_refinement_compatible(), libMesh::FEMSystem::mesh_position_set(), 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::petsc_auto_fieldsplit(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::MeshBase::prepare_for_use(), libMesh::Nemesis_IO::read(), libMesh::XdrIO::read(), libMesh::CheckpointIO::read_header(), libMesh::XdrIO::read_header(), libMesh::System::read_header(), libMesh::System::read_legacy_data(), libMesh::System::read_SCALAR_dofs(), libMesh::XdrIO::read_serialized_bc_names(), libMesh::XdrIO::read_serialized_bcs_helper(), libMesh::System::read_serialized_blocked_dof_objects(), libMesh::XdrIO::read_serialized_connectivity(), libMesh::XdrIO::read_serialized_nodes(), libMesh::XdrIO::read_serialized_nodesets(), libMesh::XdrIO::read_serialized_subdomain_names(), libMesh::System::read_serialized_vector(), libMesh::MeshBase::recalculate_n_partitions(), libMesh::MeshRefinement::refine_and_coarsen_elements(), libMesh::DistributedMesh::renumber_dof_objects(), libMesh::CheckpointIO::select_split_config(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::PetscDMWrapper::set_point_range_in_section(), libMesh::PetscDiffSolver::setup_petsc_data(), libMesh::LaplaceMeshSmoother::smooth(), libMesh::split_mesh(), libMesh::MeshBase::subdomain_ids(), libMesh::BoundaryInfo::sync(), libMesh::MeshRefinement::test_level_one(), libMesh::MeshRefinement::test_unflagged(), libMesh::MeshTools::total_weight(), libMesh::MeshRefinement::uniformly_coarsen(), libMesh::NameBasedIO::write(), libMesh::XdrIO::write(), libMesh::System::write_SCALAR_dofs(), libMesh::XdrIO::write_serialized_bcs_helper(), libMesh::System::write_serialized_blocked_dof_objects(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::XdrIO::write_serialized_nodes(), and libMesh::XdrIO::write_serialized_nodesets().

90  { return _communicator; }
const Parallel::Communicator & _communicator

◆ disable_print_counter_info()

void libMesh::ReferenceCounter::disable_print_counter_info ( )
staticinherited

Definition at line 106 of file reference_counter.C.

References libMesh::ReferenceCounter::_enable_print_counter.

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

107 {
108  _enable_print_counter = false;
109  return;
110 }

◆ enable_print_counter_info()

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.

101 {
102  _enable_print_counter = true;
103  return;
104 }

◆ get_converged_reason()

template<typename T >
LinearConvergenceReason libMesh::PetscLinearSolver< T >::get_converged_reason ( ) const
overridevirtual
Returns
The solver's convergence flag

Implements libMesh::LinearSolver< T >.

Definition at line 1684 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.

1685 {
1686  KSPConvergedReason reason;
1687  KSPGetConvergedReason(_ksp, &reason);
1688 
1689  switch(reason)
1690  {
1691 #if !PETSC_VERSION_LESS_THAN(3,2,0)
1692  case KSP_CONVERGED_RTOL_NORMAL : return CONVERGED_RTOL_NORMAL;
1693  case KSP_CONVERGED_ATOL_NORMAL : return CONVERGED_ATOL_NORMAL;
1694 #endif
1695  case KSP_CONVERGED_RTOL : return CONVERGED_RTOL;
1696  case KSP_CONVERGED_ATOL : return CONVERGED_ATOL;
1697  case KSP_CONVERGED_ITS : return CONVERGED_ITS;
1698  case KSP_CONVERGED_CG_NEG_CURVE : return CONVERGED_CG_NEG_CURVE;
1699  case KSP_CONVERGED_CG_CONSTRAINED : return CONVERGED_CG_CONSTRAINED;
1700  case KSP_CONVERGED_STEP_LENGTH : return CONVERGED_STEP_LENGTH;
1701  case KSP_CONVERGED_HAPPY_BREAKDOWN : return CONVERGED_HAPPY_BREAKDOWN;
1702  case KSP_DIVERGED_NULL : return DIVERGED_NULL;
1703  case KSP_DIVERGED_ITS : return DIVERGED_ITS;
1704  case KSP_DIVERGED_DTOL : return DIVERGED_DTOL;
1705  case KSP_DIVERGED_BREAKDOWN : return DIVERGED_BREAKDOWN;
1706  case KSP_DIVERGED_BREAKDOWN_BICG : return DIVERGED_BREAKDOWN_BICG;
1707  case KSP_DIVERGED_NONSYMMETRIC : return DIVERGED_NONSYMMETRIC;
1708  case KSP_DIVERGED_INDEFINITE_PC : return DIVERGED_INDEFINITE_PC;
1709 #if PETSC_VERSION_LESS_THAN(3,4,0)
1710  case KSP_DIVERGED_NAN : return DIVERGED_NAN;
1711 #else
1712  case KSP_DIVERGED_NANORINF : return DIVERGED_NAN;
1713 #endif
1714  case KSP_DIVERGED_INDEFINITE_MAT : return DIVERGED_INDEFINITE_MAT;
1715  case KSP_CONVERGED_ITERATING : return CONVERGED_ITERATING;
1716 #if !PETSC_VERSION_LESS_THAN(3,7,0)
1717  case KSP_DIVERGED_PCSETUP_FAILED : return DIVERGED_PCSETUP_FAILED;
1718 #endif
1719  default :
1720  libMesh::err << "Unknown convergence flag!" << std::endl;
1721  return UNKNOWN_FLAG;
1722  }
1723 }
OStreamProxy err(std::cerr)

◆ get_info()

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 (const auto & pr : _counts)
59  {
60  const std::string name(pr.first);
61  const unsigned int creations = pr.second.first;
62  const unsigned int destructions = pr.second.second;
63 
64  oss << "| " << name << " reference count information:\n"
65  << "| Creations: " << creations << '\n'
66  << "| Destructions: " << destructions << '\n';
67  }
68 
69  oss << " ---------------------------------------------------------------------------- \n";
70 
71  return oss.str();
72 
73 #else
74 
75  return "";
76 
77 #endif
78 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:42

◆ get_initial_residual()

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 1570 of file petsc_linear_solver.C.

References libMesh::err, and ierr.

1571 {
1572  PetscErrorCode ierr = 0;
1573  PetscInt its = 0;
1574 
1575  // Fill the residual history vector with the residual norms
1576  // Note that GetResidualHistory() does not copy any values, it
1577  // simply sets the pointer p. Note that for some Krylov subspace
1578  // methods, the number of residuals returned in the history
1579  // vector may be different from what you are expecting. For
1580  // example, TFQMR returns two residual values per iteration step.
1581  PetscReal * p;
1582  ierr = KSPGetResidualHistory(_ksp, &p, &its);
1583  LIBMESH_CHKERR(ierr);
1584 
1585  // Check no residual history
1586  if (its == 0)
1587  {
1588  libMesh::err << "No iterations have been performed, returning 0." << std::endl;
1589  return 0.;
1590  }
1591 
1592  // Otherwise, return the value pointed to by p.
1593  return *p;
1594 }
OStreamProxy err(std::cerr)
PetscErrorCode ierr

◆ get_residual_history()

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 1537 of file petsc_linear_solver.C.

References ierr.

1538 {
1539  PetscErrorCode ierr = 0;
1540  PetscInt its = 0;
1541 
1542  // Fill the residual history vector with the residual norms
1543  // Note that GetResidualHistory() does not copy any values, it
1544  // simply sets the pointer p. Note that for some Krylov subspace
1545  // methods, the number of residuals returned in the history
1546  // vector may be different from what you are expecting. For
1547  // example, TFQMR returns two residual values per iteration step.
1548  PetscReal * p;
1549  ierr = KSPGetResidualHistory(_ksp, &p, &its);
1550  LIBMESH_CHKERR(ierr);
1551 
1552  // Check for early return
1553  if (its == 0) return;
1554 
1555  // Create space to store the result
1556  hist.resize(its);
1557 
1558  // Copy history into the vector provided by the user.
1559  for (PetscInt i=0; i<its; ++i)
1560  {
1561  hist[i] = *p;
1562  p++;
1563  }
1564 }
PetscErrorCode ierr

◆ get_same_preconditioner()

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 322 of file linear_solver.h.

323 {
324  return same_preconditioner;
325 }

◆ increment_constructor_count()

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 181 of file reference_counter.h.

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

Referenced by libMesh::ReferenceCountedObject< RBParametrized >::ReferenceCountedObject().

182 {
183  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
184  std::pair<unsigned int, unsigned int> & p = _counts[name];
185 
186  p.first++;
187 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:42
spin_mutex spin_mtx
Definition: threads.C:29

◆ increment_destructor_count()

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 194 of file reference_counter.h.

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

Referenced by libMesh::ReferenceCountedObject< RBParametrized >::~ReferenceCountedObject().

195 {
196  Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
197  std::pair<unsigned int, unsigned int> & p = _counts[name];
198 
199  p.second++;
200 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:42
spin_mutex spin_mtx
Definition: threads.C:29

◆ init() [1/2]

template<typename T >
void libMesh::PetscLinearSolver< T >::init ( const char *  name = nullptr)
overridevirtual

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 192 of file petsc_linear_solver.C.

References libMesh::libMeshPrivateData::_is_initialized, ierr, libMesh::initialized(), libMesh::libmesh_petsc_preconditioner_apply(), libMesh::libmesh_petsc_preconditioner_setup(), libMesh::Quality::name(), and libMesh::PetscPreconditioner< T >::set_petsc_preconditioner_type().

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

193 {
194  // Initialize the data structures if not done so already.
195  if (!this->initialized())
196  {
197  this->_is_initialized = true;
198 
199  PetscErrorCode ierr=0;
200 
201  // Create the linear solver context
202  ierr = KSPCreate (this->comm().get(), &_ksp);
203  LIBMESH_CHKERR(ierr);
204 
205  if (name)
206  {
207  ierr = KSPSetOptionsPrefix(_ksp, name);
208  LIBMESH_CHKERR(ierr);
209  }
210 
211  // Create the preconditioner context
212  ierr = KSPGetPC (_ksp, &_pc);
213  LIBMESH_CHKERR(ierr);
214 
215  // Set user-specified solver and preconditioner types
216  this->set_petsc_solver_type();
217 
218  // If the SolverConfiguration object is provided, use it to set
219  // options during solver initialization.
220  if (this->_solver_configuration)
221  {
223  }
224 
225  // Set the options from user-input
226  // Set runtime options, e.g.,
227  // -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
228  // These options will override those specified above as long as
229  // KSPSetFromOptions() is called _after_ any other customization
230  // routines.
231 
232  ierr = KSPSetFromOptions (_ksp);
233  LIBMESH_CHKERR(ierr);
234 
235  // Not sure if this is necessary, or if it is already handled by KSPSetFromOptions?
236  // NOT NECESSARY!!!!
237  //ierr = PCSetFromOptions (_pc);
238  //LIBMESH_CHKERR(ierr);
239 
240  // Have the Krylov subspace method use our good initial guess
241  // rather than 0, unless the user requested a KSPType of
242  // preonly, which complains if asked to use initial guesses.
243 #if PETSC_VERSION_LESS_THAN(3,0,0) || !PETSC_RELEASE_LESS_THAN(3,4,0)
244  // Pre-3.0 and petsc-dev (as of October 2012) use non-const versions
245  KSPType ksp_type;
246 #else
247  const KSPType ksp_type;
248 #endif
249 
250  ierr = KSPGetType (_ksp, &ksp_type);
251  LIBMESH_CHKERR(ierr);
252 
253  if (strcmp(ksp_type, "preonly"))
254  {
255  ierr = KSPSetInitialGuessNonzero (_ksp, PETSC_TRUE);
256  LIBMESH_CHKERR(ierr);
257  }
258 
259  // Notify PETSc of location to store residual history.
260  // This needs to be called before any solves, since
261  // it sets the residual history length to zero. The default
262  // behavior is for PETSc to allocate (internally) an array
263  // of size 1000 to hold the residual norm history.
264  ierr = KSPSetResidualHistory(_ksp,
265  PETSC_NULL, // pointer to the array which holds the history
266  PETSC_DECIDE, // size of the array holding the history
267  PETSC_TRUE); // Whether or not to reset the history for each solve.
268  LIBMESH_CHKERR(ierr);
269 
271 
272  //If there is a preconditioner object we need to set the internal setup and apply routines
273  if (this->_preconditioner)
274  {
275  this->_preconditioner->init();
276  PCShellSetContext(_pc,(void *)this->_preconditioner);
277  PCShellSetSetUp(_pc,libmesh_petsc_preconditioner_setup);
278  PCShellSetApply(_pc,libmesh_petsc_preconditioner_apply);
279  }
280  }
281 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:42
Preconditioner< T > * _preconditioner
const Parallel::Communicator & comm() const
bool initialized() const
Definition: linear_solver.h:96
PetscErrorCode libmesh_petsc_preconditioner_setup(void *ctx)
SolverConfiguration * _solver_configuration
static void set_petsc_preconditioner_type(const PreconditionerType &preconditioner_type, PC &pc)
PetscErrorCode ierr
PetscErrorCode libmesh_petsc_preconditioner_apply(void *ctx, Vec x, Vec y)
PreconditionerType _preconditioner_type

◆ init() [2/2]

template<typename T >
void libMesh::PetscLinearSolver< T >::init ( PetscMatrix< T > *  matrix,
const char *  name = nullptr 
)

Initialize data structures if not done so already plus much more

Definition at line 285 of file petsc_linear_solver.C.

References libMesh::libMeshPrivateData::_is_initialized, ierr, libMesh::initialized(), libMesh::libmesh_petsc_preconditioner_apply(), libMesh::libmesh_petsc_preconditioner_setup(), libMesh::PetscMatrix< T >::mat(), libMesh::Quality::name(), and libMesh::PetscPreconditioner< T >::set_petsc_preconditioner_type().

287 {
288  // Initialize the data structures if not done so already.
289  if (!this->initialized())
290  {
291  this->_is_initialized = true;
292 
293  PetscErrorCode ierr=0;
294 
295  // Create the linear solver context
296  ierr = KSPCreate (this->comm().get(), &_ksp);
297  LIBMESH_CHKERR(ierr);
298 
299  if (name)
300  {
301  ierr = KSPSetOptionsPrefix(_ksp, name);
302  LIBMESH_CHKERR(ierr);
303  }
304 
305  //ierr = PCCreate (this->comm().get(), &_pc);
306  // LIBMESH_CHKERR(ierr);
307 
308  // Create the preconditioner context
309  ierr = KSPGetPC (_ksp, &_pc);
310  LIBMESH_CHKERR(ierr);
311 
312  // Set operators. The input matrix works as the preconditioning matrix
313 #if PETSC_RELEASE_LESS_THAN(3,5,0)
314  ierr = KSPSetOperators(_ksp, matrix->mat(), matrix->mat(),DIFFERENT_NONZERO_PATTERN);
315 #else
316  ierr = KSPSetOperators(_ksp, matrix->mat(), matrix->mat());
317 #endif
318  LIBMESH_CHKERR(ierr);
319 
320  // Set user-specified solver and preconditioner types
321  this->set_petsc_solver_type();
322 
323  // If the SolverConfiguration object is provided, use it to set
324  // options during solver initialization.
325  if (this->_solver_configuration)
326  {
328  }
329 
330  // Set the options from user-input
331  // Set runtime options, e.g.,
332  // -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
333  // These options will override those specified above as long as
334  // KSPSetFromOptions() is called _after_ any other customization
335  // routines.
336 
337  ierr = KSPSetFromOptions (_ksp);
338  LIBMESH_CHKERR(ierr);
339 
340  // Not sure if this is necessary, or if it is already handled by KSPSetFromOptions?
341  // NOT NECESSARY!!!!
342  //ierr = PCSetFromOptions (_pc);
343  //LIBMESH_CHKERR(ierr);
344 
345  // Have the Krylov subspace method use our good initial guess
346  // rather than 0, unless the user requested a KSPType of
347  // preonly, which complains if asked to use initial guesses.
348 #if PETSC_VERSION_LESS_THAN(3,0,0) || !PETSC_RELEASE_LESS_THAN(3,4,0)
349  KSPType ksp_type;
350 #else
351  const KSPType ksp_type;
352 #endif
353 
354  ierr = KSPGetType (_ksp, &ksp_type);
355  LIBMESH_CHKERR(ierr);
356 
357  if (strcmp(ksp_type, "preonly"))
358  {
359  ierr = KSPSetInitialGuessNonzero (_ksp, PETSC_TRUE);
360  LIBMESH_CHKERR(ierr);
361  }
362 
363  // Notify PETSc of location to store residual history.
364  // This needs to be called before any solves, since
365  // it sets the residual history length to zero. The default
366  // behavior is for PETSc to allocate (internally) an array
367  // of size 1000 to hold the residual norm history.
368  ierr = KSPSetResidualHistory(_ksp,
369  PETSC_NULL, // pointer to the array which holds the history
370  PETSC_DECIDE, // size of the array holding the history
371  PETSC_TRUE); // Whether or not to reset the history for each solve.
372  LIBMESH_CHKERR(ierr);
373 
375  if (this->_preconditioner)
376  {
377  this->_preconditioner->set_matrix(*matrix);
378  this->_preconditioner->init();
379  PCShellSetContext(_pc,(void *)this->_preconditioner);
380  PCShellSetSetUp(_pc,libmesh_petsc_preconditioner_setup);
381  PCShellSetApply(_pc,libmesh_petsc_preconditioner_apply);
382  }
383  }
384 }
std::string name(const ElemQuality q)
Definition: elem_quality.C:42
Preconditioner< T > * _preconditioner
const Parallel::Communicator & comm() const
bool initialized() const
Definition: linear_solver.h:96
PetscErrorCode libmesh_petsc_preconditioner_setup(void *ctx)
SolverConfiguration * _solver_configuration
static void set_petsc_preconditioner_type(const PreconditionerType &preconditioner_type, PC &pc)
PetscErrorCode ierr
PetscErrorCode libmesh_petsc_preconditioner_apply(void *ctx, Vec x, Vec y)
PreconditionerType _preconditioner_type

◆ init_names()

template<typename T >
void libMesh::PetscLinearSolver< T >::init_names ( const System sys)
overridevirtual

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 390 of file petsc_linear_solver.C.

References libMesh::petsc_auto_fieldsplit().

391 {
392  petsc_auto_fieldsplit(this->pc(), sys);
393 }
void petsc_auto_fieldsplit(PC my_pc, const System &sys)

◆ initialized()

template<typename T>
bool libMesh::LinearSolver< T >::initialized ( ) const
inlineinherited
Returns
true if the data structures are initialized, false otherwise.

Definition at line 96 of file linear_solver.h.

96 { return _is_initialized; }

◆ ksp()

template<typename T >
KSP libMesh::PetscLinearSolver< T >::ksp ( )
inline
Returns
The raw PETSc ksp context pointer.

This is useful if you are for example setting a custom convergence test with KSPSetConvergenceTest().

Definition at line 240 of file petsc_linear_solver.h.

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

240 { this->init(); return _ksp; }
virtual void init(const char *name=nullptr) override

◆ n_objects()

static unsigned int libMesh::ReferenceCounter::n_objects ( )
inlinestaticinherited

Prints the number of outstanding (created, but not yet destroyed) objects.

Definition at line 83 of file reference_counter.h.

References libMesh::ReferenceCounter::_n_objects.

84  { return _n_objects; }
static Threads::atomic< unsigned int > _n_objects

◆ n_processors()

processor_id_type libMesh::ParallelObject::n_processors ( ) const
inlineinherited
Returns
The number of processors in the group.

Definition at line 95 of file parallel_object.h.

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

Referenced by libMesh::BoundaryInfo::_find_id_maps(), libMesh::PetscDMWrapper::add_dofs_to_section(), libMesh::DistributedMesh::add_elem(), libMesh::DistributedMesh::add_node(), libMesh::LaplaceMeshSmoother::allgather_graph(), libMesh::FEMSystem::assembly(), libMesh::AztecLinearSolver< T >::AztecLinearSolver(), libMesh::BoundaryInfo::build_node_list_from_side_list(), libMesh::EquationSystems::build_parallel_elemental_solution_vector(), libMesh::DistributedMesh::clear(), libMesh::Nemesis_IO_Helper::compute_border_node_ids(), libMesh::Nemesis_IO_Helper::construct_nemesis_filename(), libMesh::UnstructuredMesh::create_pid_mesh(), libMesh::MeshTools::create_processor_bounding_box(), libMesh::DofMap::distribute_dofs(), libMesh::DofMap::distribute_local_dofs_node_major(), libMesh::DofMap::distribute_local_dofs_var_major(), libMesh::EnsightIO::EnsightIO(), libMesh::MeshBase::get_info(), libMesh::SystemSubsetBySubdomain::init(), libMesh::PetscDMWrapper::init_and_attach_petscdm(), libMesh::Nemesis_IO_Helper::initialize(), libMesh::DistributedMesh::insert_elem(), libMesh::MeshTools::libmesh_assert_contiguous_dof_ids(), libMesh::MeshTools::libmesh_assert_parallel_consistent_new_node_procids(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Elem >(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Node >(), libMesh::MeshTools::libmesh_assert_topology_consistent_procids< Node >(), libMesh::MeshTools::libmesh_assert_valid_boundary_ids(), libMesh::MeshTools::libmesh_assert_valid_dof_ids(), libMesh::MeshTools::libmesh_assert_valid_neighbors(), libMesh::MeshTools::libmesh_assert_valid_refinement_flags(), libMesh::DofMap::local_variable_indices(), libMesh::MeshRefinement::make_coarsening_compatible(), libMesh::MeshBase::n_active_elem_on_proc(), libMesh::MeshBase::n_elem_on_proc(), libMesh::MeshBase::n_nodes_on_proc(), libMesh::MeshBase::partition(), libMesh::PetscLinearSolver< T >::PetscLinearSolver(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::NameBasedIO::read(), libMesh::Nemesis_IO::read(), libMesh::CheckpointIO::read(), libMesh::CheckpointIO::read_connectivity(), libMesh::XdrIO::read_header(), libMesh::CheckpointIO::read_nodes(), libMesh::System::read_parallel_data(), libMesh::System::read_SCALAR_dofs(), libMesh::System::read_serialized_blocked_dof_objects(), libMesh::System::read_serialized_vector(), libMesh::DistributedMesh::renumber_dof_objects(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::PetscDMWrapper::set_point_range_in_section(), libMesh::MeshRefinement::uniformly_coarsen(), libMesh::DistributedMesh::update_parallel_id_counts(), libMesh::GMVIO::write_binary(), libMesh::GMVIO::write_discontinuous_gmv(), libMesh::System::write_parallel_data(), libMesh::System::write_SCALAR_dofs(), libMesh::XdrIO::write_serialized_bcs_helper(), libMesh::System::write_serialized_blocked_dof_objects(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::XdrIO::write_serialized_nodes(), and libMesh::XdrIO::write_serialized_nodesets().

96  { return cast_int<processor_id_type>(_communicator.size()); }
processor_id_type size() const
Definition: communicator.h:175
const Parallel::Communicator & _communicator

◆ pc()

template<typename T >
PC libMesh::PetscLinearSolver< T >::pc ( )
inline
Returns
The raw PETSc preconditioner context pointer.

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

Definition at line 232 of file petsc_linear_solver.h.

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

232 { this->init(); return _pc; }
virtual void init(const char *name=nullptr) override

◆ preconditioner_type()

template<typename T >
PreconditionerType libMesh::LinearSolver< T >::preconditioner_type ( ) const
inherited
Returns
The type of preconditioner to use.

Definition at line 98 of file linear_solver.C.

99 {
100  if (_preconditioner)
101  return _preconditioner->type();
102 
103  return _preconditioner_type;
104 }
Preconditioner< T > * _preconditioner
PreconditionerType _preconditioner_type

◆ print_converged_reason()

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 170 of file linear_solver.C.

171 {
173  libMesh::out << "Linear solver convergence/divergence reason: " << Utility::enum_to_string(reason) << std::endl;
174 }
std::string enum_to_string(const T e)
virtual LinearConvergenceReason get_converged_reason() const =0
OStreamProxy out(std::cout)

◆ print_info()

void libMesh::ReferenceCounter::print_info ( std::ostream &  out = libMesh::out)
staticinherited

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

Definition at line 87 of file reference_counter.C.

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

88 {
90  out_stream << ReferenceCounter::get_info();
91 }
static std::string get_info()

◆ processor_id()

processor_id_type libMesh::ParallelObject::processor_id ( ) const
inlineinherited
Returns
The rank of this processor in the group.

Definition at line 101 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::PetscDMWrapper::add_dofs_to_section(), libMesh::DistributedMesh::add_elem(), libMesh::BoundaryInfo::add_elements(), libMesh::DofMap::add_neighbors_to_send_list(), libMesh::DistributedMesh::add_node(), libMesh::UnstructuredMesh::all_second_order(), libMesh::MeshTools::Modification::all_tri(), libMesh::FEMSystem::assembly(), libMesh::EquationSystems::build_discontinuous_solution_vector(), libMesh::Nemesis_IO_Helper::build_element_and_node_maps(), libMesh::InfElemBuilder::build_inf_elem(), libMesh::BoundaryInfo::build_node_list_from_side_list(), libMesh::EquationSystems::build_parallel_elemental_solution_vector(), 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::DofMap::end_dof(), libMesh::DofMap::end_old_dof(), libMesh::EnsightIO::EnsightIO(), libMesh::MeshFunction::find_element(), libMesh::MeshFunction::find_elements(), libMesh::UnstructuredMesh::find_neighbors(), libMesh::DofMap::first_dof(), libMesh::DofMap::first_old_dof(), 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::Nemesis_IO_Helper::get_ss_param_global(), libMesh::SparsityPattern::Build::handle_vi_vj(), libMesh::SystemSubsetBySubdomain::init(), 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::DofMap::last_dof(), libMesh::MeshTools::libmesh_assert_consistent_distributed(), libMesh::MeshTools::libmesh_assert_consistent_distributed_nodes(), libMesh::MeshTools::libmesh_assert_contiguous_dof_ids(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Elem >(), 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::DofMap::n_local_dofs(), 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::SparsityPattern::Build::operator()(), libMesh::DistributedMesh::own_node(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::Nemesis_IO_Helper::put_cmap_params(), libMesh::Nemesis_IO_Helper::put_elem_cmap(), libMesh::Nemesis_IO_Helper::put_elem_map(), libMesh::Nemesis_IO_Helper::put_loadbal_param(), libMesh::Nemesis_IO_Helper::put_node_cmap(), libMesh::Nemesis_IO_Helper::put_node_map(), libMesh::NameBasedIO::read(), libMesh::Nemesis_IO::read(), libMesh::XdrIO::read(), libMesh::CheckpointIO::read(), libMesh::ExodusII_IO_Helper::read_elem_num_map(), libMesh::ExodusII_IO_Helper::read_global_values(), libMesh::CheckpointIO::read_header(), libMesh::XdrIO::read_header(), libMesh::System::read_header(), libMesh::System::read_legacy_data(), libMesh::ExodusII_IO_Helper::read_node_num_map(), libMesh::System::read_parallel_data(), libMesh::System::read_SCALAR_dofs(), libMesh::XdrIO::read_serialized_bc_names(), libMesh::XdrIO::read_serialized_bcs_helper(), libMesh::System::read_serialized_blocked_dof_objects(), libMesh::XdrIO::read_serialized_connectivity(), libMesh::System::read_serialized_data(), libMesh::XdrIO::read_serialized_nodes(), libMesh::XdrIO::read_serialized_nodesets(), libMesh::XdrIO::read_serialized_subdomain_names(), libMesh::System::read_serialized_vector(), libMesh::System::read_serialized_vectors(), libMesh::DistributedMesh::renumber_dof_objects(), libMesh::CheckpointIO::select_split_config(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::PetscDMWrapper::set_point_range_in_section(), libMesh::LaplaceMeshSmoother::smooth(), libMesh::MeshTools::total_weight(), libMesh::MeshRefinement::uniformly_coarsen(), libMesh::Parallel::Packing< T >::unpack(), libMesh::DistributedMesh::update_parallel_id_counts(), libMesh::NameBasedIO::write(), libMesh::XdrIO::write(), libMesh::CheckpointIO::write(), libMesh::EquationSystems::write(), libMesh::GMVIO::write_discontinuous_gmv(), libMesh::ExodusII_IO::write_element_data(), libMesh::ExodusII_IO_Helper::write_element_values(), libMesh::ExodusII_IO_Helper::write_elements(), libMesh::ExodusII_IO::write_global_data(), libMesh::ExodusII_IO_Helper::write_global_values(), libMesh::System::write_header(), 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::System::write_parallel_data(), libMesh::System::write_SCALAR_dofs(), libMesh::XdrIO::write_serialized_bc_names(), libMesh::XdrIO::write_serialized_bcs_helper(), libMesh::System::write_serialized_blocked_dof_objects(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::System::write_serialized_data(), libMesh::XdrIO::write_serialized_nodes(), libMesh::XdrIO::write_serialized_nodesets(), libMesh::XdrIO::write_serialized_subdomain_names(), libMesh::System::write_serialized_vector(), libMesh::System::write_serialized_vectors(), libMesh::Nemesis_IO_Helper::write_sidesets(), libMesh::ExodusII_IO_Helper::write_sidesets(), libMesh::ExodusII_IO::write_timestep(), libMesh::ExodusII_IO_Helper::write_timestep(), and libMesh::ExodusII_IO::write_timestep_discontinuous().

102  { return cast_int<processor_id_type>(_communicator.rank()); }
const Parallel::Communicator & _communicator
processor_id_type rank() const
Definition: communicator.h:173

◆ restrict_solve_to()

template<typename T >
void libMesh::PetscLinearSolver< T >::restrict_solve_to ( const std::vector< unsigned int > *const  dofs,
const SubsetSolveMode  subset_solve_mode = SUBSET_ZERO 
)
overridevirtual

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 nullptr.

Reimplemented from libMesh::LinearSolver< T >.

Definition at line 399 of file petsc_linear_solver.C.

References ierr, and PETSC_OWN_POINTER.

401 {
402  PetscErrorCode ierr=0;
403 
404  /* The preconditioner (in particular if a default preconditioner)
405  will have to be reset. We call this->clear() to do that. This
406  call will also remove and free any previous subset that may have
407  been set before. */
408  this->clear();
409 
410  _subset_solve_mode = subset_solve_mode;
411 
412  if (dofs != nullptr)
413  {
414  PetscInt * petsc_dofs = nullptr;
415  ierr = PetscMalloc(dofs->size()*sizeof(PetscInt), &petsc_dofs);
416  LIBMESH_CHKERR(ierr);
417 
418  for (std::size_t i=0; i<dofs->size(); i++)
419  petsc_dofs[i] = (*dofs)[i];
420 
421  ierr = ISCreateLibMesh(this->comm().get(),
422  cast_int<PetscInt>(dofs->size()),
423  petsc_dofs, PETSC_OWN_POINTER,
425  LIBMESH_CHKERR(ierr);
426  }
427 }
const Parallel::Communicator & comm() const
PetscErrorCode ierr
virtual void clear() override

◆ reuse_preconditioner()

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 129 of file linear_solver.C.

Referenced by libMesh::ImplicitSystem::disable_cache().

130 {
131  same_preconditioner = reuse_flag;
132 }

◆ set_petsc_solver_type()

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 1600 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.

1601 {
1602  PetscErrorCode ierr = 0;
1603 
1604  switch (this->_solver_type)
1605  {
1606 
1607  case CG:
1608  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPCG));
1609  LIBMESH_CHKERR(ierr);
1610  return;
1611 
1612  case CR:
1613  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPCR));
1614  LIBMESH_CHKERR(ierr);
1615  return;
1616 
1617  case CGS:
1618  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPCGS));
1619  LIBMESH_CHKERR(ierr);
1620  return;
1621 
1622  case BICG:
1623  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPBICG));
1624  LIBMESH_CHKERR(ierr);
1625  return;
1626 
1627  case TCQMR:
1628  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPTCQMR));
1629  LIBMESH_CHKERR(ierr);
1630  return;
1631 
1632  case TFQMR:
1633  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPTFQMR));
1634  LIBMESH_CHKERR(ierr);
1635  return;
1636 
1637  case LSQR:
1638  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPLSQR));
1639  LIBMESH_CHKERR(ierr);
1640  return;
1641 
1642  case BICGSTAB:
1643  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPBCGS));
1644  LIBMESH_CHKERR(ierr);
1645  return;
1646 
1647  case MINRES:
1648  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPMINRES));
1649  LIBMESH_CHKERR(ierr);
1650  return;
1651 
1652  case GMRES:
1653  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPGMRES));
1654  LIBMESH_CHKERR(ierr);
1655  return;
1656 
1657  case RICHARDSON:
1658  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPRICHARDSON));
1659  LIBMESH_CHKERR(ierr);
1660  return;
1661 
1662  case CHEBYSHEV:
1663 #if defined(LIBMESH_HAVE_PETSC) && PETSC_VERSION_LESS_THAN(3,3,0)
1664  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPCHEBYCHEV));
1665  LIBMESH_CHKERR(ierr);
1666  return;
1667 #else
1668  ierr = KSPSetType (_ksp, const_cast<KSPType>(KSPCHEBYSHEV));
1669  LIBMESH_CHKERR(ierr);
1670  return;
1671 #endif
1672 
1673 
1674  default:
1675  libMesh::err << "ERROR: Unsupported PETSC Solver: "
1676  << Utility::enum_to_string(this->_solver_type) << std::endl
1677  << "Continuing with PETSC defaults" << std::endl;
1678  }
1679 }
OStreamProxy err(std::cerr)
std::string enum_to_string(const T e)
PetscErrorCode ierr

◆ set_preconditioner_type()

template<typename T >
void libMesh::LinearSolver< T >::set_preconditioner_type ( const PreconditionerType  pct)
inherited

Sets the type of preconditioner to use.

Definition at line 108 of file linear_solver.C.

109 {
110  if (_preconditioner)
111  _preconditioner->set_type(pct);
112  else
113  _preconditioner_type = pct;
114 }
Preconditioner< T > * _preconditioner
PreconditionerType _preconditioner_type

◆ set_solver_configuration()

template<typename T >
void libMesh::LinearSolver< T >::set_solver_configuration ( SolverConfiguration solver_configuration)
inherited

Set the solver configuration object.

Definition at line 177 of file linear_solver.C.

178 {
179  _solver_configuration = &solver_configuration;
180 }
SolverConfiguration * _solver_configuration

◆ set_solver_type()

template<typename T>
void libMesh::LinearSolver< T >::set_solver_type ( const SolverType  st)
inlineinherited

Sets the type of solver to use.

Definition at line 127 of file linear_solver.h.

128  { _solver_type = st; }

◆ solve() [1/6]

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 
)
inlineoverridevirtual

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 157 of file petsc_linear_solver.h.

162  {
163  return this->solve(matrix_in, matrix_in, solution_in, rhs_in, tol, m_its);
164  }
virtual std::pair< unsigned int, Real > solve(SparseMatrix< T > &matrix_in, NumericVector< T > &solution_in, NumericVector< T > &rhs_in, const double tol, const unsigned int m_its) override

◆ solve() [2/6]

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 
)
overridevirtual

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

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

Implements libMesh::LinearSolver< T >.

Definition at line 433 of file petsc_linear_solver.C.

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

439 {
440  LOG_SCOPE("solve()", "PetscLinearSolver");
441 
442  // Make sure the data passed in are really of Petsc types
443  PetscMatrix<T> * matrix = cast_ptr<PetscMatrix<T> *>(&matrix_in);
444  PetscMatrix<T> * precond = cast_ptr<PetscMatrix<T> *>(&precond_in);
445  PetscVector<T> * solution = cast_ptr<PetscVector<T> *>(&solution_in);
446  PetscVector<T> * rhs = cast_ptr<PetscVector<T> *>(&rhs_in);
447 
448  this->init (matrix);
449 
450  PetscErrorCode ierr=0;
451  PetscInt its=0, max_its = static_cast<PetscInt>(m_its);
452  PetscReal final_resid=0.;
453 
454  // Close the matrices and vectors in case this wasn't already done.
455  matrix->close ();
456  precond->close ();
457  solution->close ();
458  rhs->close ();
459 
460  // // If matrix != precond, then this means we have specified a
461  // // special preconditioner, so reset preconditioner type to PCMAT.
462  // if (matrix != precond)
463  // {
464  // this->_preconditioner_type = USER_PRECOND;
465  // this->set_petsc_preconditioner_type ();
466  // }
467 
468  Mat submat = nullptr;
469  Mat subprecond = nullptr;
470  Vec subrhs = nullptr;
471  Vec subsolution = nullptr;
472  VecScatter scatter = nullptr;
473  std::unique_ptr<PetscMatrix<Number>> subprecond_matrix;
474 
475  // Set operators. Also restrict rhs and solution vector to
476  // subdomain if necessary.
477  if (_restrict_solve_to_is != nullptr)
478  {
479  PetscInt is_local_size = this->_restrict_solve_to_is_local_size();
480 
481  ierr = VecCreate(this->comm().get(),&subrhs);
482  LIBMESH_CHKERR(ierr);
483  ierr = VecSetSizes(subrhs,is_local_size,PETSC_DECIDE);
484  LIBMESH_CHKERR(ierr);
485  ierr = VecSetFromOptions(subrhs);
486  LIBMESH_CHKERR(ierr);
487 
488  ierr = VecCreate(this->comm().get(),&subsolution);
489  LIBMESH_CHKERR(ierr);
490  ierr = VecSetSizes(subsolution,is_local_size,PETSC_DECIDE);
491  LIBMESH_CHKERR(ierr);
492  ierr = VecSetFromOptions(subsolution);
493  LIBMESH_CHKERR(ierr);
494 
495  ierr = LibMeshVecScatterCreate(rhs->vec(), _restrict_solve_to_is, subrhs, nullptr, &scatter);
496  LIBMESH_CHKERR(ierr);
497 
498  ierr = VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
499  LIBMESH_CHKERR(ierr);
500  ierr = VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
501  LIBMESH_CHKERR(ierr);
502 
503  ierr = VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
504  LIBMESH_CHKERR(ierr);
505  ierr = VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
506  LIBMESH_CHKERR(ierr);
507 
508  ierr = LibMeshCreateSubMatrix(matrix->mat(),
511 #if PETSC_VERSION_LESS_THAN(3,1,0)
512  PETSC_DECIDE,
513 #endif
514  MAT_INITIAL_MATRIX,
515  &submat);
516  LIBMESH_CHKERR(ierr);
517 
518  ierr = LibMeshCreateSubMatrix(precond->mat(),
521 #if PETSC_VERSION_LESS_THAN(3,1,0)
522  PETSC_DECIDE,
523 #endif
524  MAT_INITIAL_MATRIX,
525  &subprecond);
526  LIBMESH_CHKERR(ierr);
527 
528  /* Since removing columns of the matrix changes the equation
529  system, we will now change the right hand side to compensate
530  for this. Note that this is not necessary if \p SUBSET_ZERO
531  has been selected. */
533  {
534  _create_complement_is(rhs_in);
535  PetscInt is_complement_local_size =
536  cast_int<PetscInt>(rhs_in.local_size()-is_local_size);
537 
538  Vec subvec1 = nullptr;
539  Mat submat1 = nullptr;
540  VecScatter scatter1 = nullptr;
541 
542  ierr = VecCreate(this->comm().get(),&subvec1);
543  LIBMESH_CHKERR(ierr);
544  ierr = VecSetSizes(subvec1,is_complement_local_size,PETSC_DECIDE);
545  LIBMESH_CHKERR(ierr);
546  ierr = VecSetFromOptions(subvec1);
547  LIBMESH_CHKERR(ierr);
548 
549  ierr = LibMeshVecScatterCreate(rhs->vec(),_restrict_solve_to_is_complement, subvec1, nullptr, &scatter1);
550  LIBMESH_CHKERR(ierr);
551 
552  ierr = VecScatterBegin(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
553  LIBMESH_CHKERR(ierr);
554  ierr = VecScatterEnd(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
555  LIBMESH_CHKERR(ierr);
556 
557  ierr = VecScale(subvec1,-1.0);
558  LIBMESH_CHKERR(ierr);
559 
560  ierr = LibMeshCreateSubMatrix(matrix->mat(),
563 #if PETSC_VERSION_LESS_THAN(3,1,0)
564  PETSC_DECIDE,
565 #endif
566  MAT_INITIAL_MATRIX,
567  &submat1);
568  LIBMESH_CHKERR(ierr);
569 
570  ierr = MatMultAdd(submat1,subvec1,subrhs,subrhs);
571  LIBMESH_CHKERR(ierr);
572 
573  ierr = LibMeshVecScatterDestroy(&scatter1);
574  LIBMESH_CHKERR(ierr);
575  ierr = LibMeshVecDestroy(&subvec1);
576  LIBMESH_CHKERR(ierr);
577  ierr = LibMeshMatDestroy(&submat1);
578  LIBMESH_CHKERR(ierr);
579  }
580 #if PETSC_RELEASE_LESS_THAN(3,5,0)
581  ierr = KSPSetOperators(_ksp, submat, subprecond,
582  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
583 #else
584  ierr = KSPSetOperators(_ksp, submat, subprecond);
585 
586  PetscBool ksp_reuse_preconditioner = this->same_preconditioner ? PETSC_TRUE : PETSC_FALSE;
587  ierr = KSPSetReusePreconditioner(_ksp, ksp_reuse_preconditioner);
588 #endif
589  LIBMESH_CHKERR(ierr);
590 
591  if (this->_preconditioner)
592  {
593  subprecond_matrix.reset(new PetscMatrix<Number>(subprecond, this->comm()));
594  this->_preconditioner->set_matrix(*subprecond_matrix);
595  this->_preconditioner->init();
596  }
597  }
598  else
599  {
600 #if PETSC_RELEASE_LESS_THAN(3,5,0)
601  ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat(),
602  this->same_preconditioner ? SAME_PRECONDITIONER : DIFFERENT_NONZERO_PATTERN);
603 #else
604  ierr = KSPSetOperators(_ksp, matrix->mat(), precond->mat());
605 
606  PetscBool ksp_reuse_preconditioner = this->same_preconditioner ? PETSC_TRUE : PETSC_FALSE;
607  ierr = KSPSetReusePreconditioner(_ksp, ksp_reuse_preconditioner);
608 #endif
609  LIBMESH_CHKERR(ierr);
610 
611  if (this->_preconditioner)
612  {
613  this->_preconditioner->set_matrix(matrix_in);
614  this->_preconditioner->init();
615  }
616  }
617 
618  // Set the tolerances for the iterative solver. Use the user-supplied
619  // tolerance for the relative residual & leave the others at default values.
620  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
621  PETSC_DEFAULT, max_its);
622  LIBMESH_CHKERR(ierr);
623 
624  // Allow command line options to override anything set programmatically.
625  ierr = KSPSetFromOptions(_ksp);
626  LIBMESH_CHKERR(ierr);
627 
628  // If the SolverConfiguration object is provided, use it to override
629  // solver options.
630  if (this->_solver_configuration)
631  {
633  }
634 
635  // Solve the linear system
636  if (_restrict_solve_to_is != nullptr)
637  {
638  ierr = KSPSolve (_ksp, subrhs, subsolution);
639  LIBMESH_CHKERR(ierr);
640  }
641  else
642  {
643  ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
644  LIBMESH_CHKERR(ierr);
645  }
646 
647  // Get the number of iterations required for convergence
648  ierr = KSPGetIterationNumber (_ksp, &its);
649  LIBMESH_CHKERR(ierr);
650 
651  // Get the norm of the final residual to return to the user.
652  ierr = KSPGetResidualNorm (_ksp, &final_resid);
653  LIBMESH_CHKERR(ierr);
654 
655  if (_restrict_solve_to_is != nullptr)
656  {
657  switch(_subset_solve_mode)
658  {
659  case SUBSET_ZERO:
660  ierr = VecZeroEntries(solution->vec());
661  LIBMESH_CHKERR(ierr);
662  break;
663 
664  case SUBSET_COPY_RHS:
665  ierr = VecCopy(rhs->vec(),solution->vec());
666  LIBMESH_CHKERR(ierr);
667  break;
668 
669  case SUBSET_DONT_TOUCH:
670  /* Nothing to do here. */
671  break;
672 
673  default:
674  libmesh_error_msg("Invalid subset solve mode = " << _subset_solve_mode);
675  }
676  ierr = VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
677  LIBMESH_CHKERR(ierr);
678  ierr = VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
679  LIBMESH_CHKERR(ierr);
680 
681  ierr = LibMeshVecScatterDestroy(&scatter);
682  LIBMESH_CHKERR(ierr);
683 
684  if (this->_preconditioner)
685  {
686  // Before subprecond_matrix gets cleaned up, we should give
687  // the _preconditioner a different matrix.
688  this->_preconditioner->set_matrix(matrix_in);
689  this->_preconditioner->init();
690  }
691 
692  ierr = LibMeshVecDestroy(&subsolution);
693  LIBMESH_CHKERR(ierr);
694  ierr = LibMeshVecDestroy(&subrhs);
695  LIBMESH_CHKERR(ierr);
696  ierr = LibMeshMatDestroy(&submat);
697  LIBMESH_CHKERR(ierr);
698  ierr = LibMeshMatDestroy(&subprecond);
699  LIBMESH_CHKERR(ierr);
700  }
701 
702  // return the # of its. and the final residual norm.
703  return std::make_pair(its, final_resid);
704 }
Preconditioner< T > * _preconditioner
virtual void init(const char *name=nullptr) override
void _create_complement_is(const NumericVector< T > &vec_in)
const Parallel::Communicator & comm() const
virtual void configure_solver()=0
SolverConfiguration * _solver_configuration
PetscTruth PetscBool
Definition: petsc_macro.h:67
PetscErrorCode ierr
PetscInt _restrict_solve_to_is_local_size() const

◆ solve() [3/6]

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 
)
overridevirtual

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

Implements libMesh::LinearSolver< T >.

Definition at line 969 of file petsc_linear_solver.C.

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

974 {
975 
976 #if PETSC_VERSION_LESS_THAN(3,1,0)
977  if (_restrict_solve_to_is != nullptr)
978  libmesh_error_msg("The current implementation of subset solves with " \
979  << "shell matrices requires PETSc version 3.1 or above. " \
980  << "Older PETSc version do not support automatic " \
981  << "submatrix generation of shell matrices.");
982 #endif
983 
984  LOG_SCOPE("solve()", "PetscLinearSolver");
985 
986  // Make sure the data passed in are really of Petsc types
987  PetscVector<T> * solution = cast_ptr<PetscVector<T> *>(&solution_in);
988  PetscVector<T> * rhs = cast_ptr<PetscVector<T> *>(&rhs_in);
989 
990  this->init ();
991 
992  PetscErrorCode ierr=0;
993  PetscInt its=0, max_its = static_cast<PetscInt>(m_its);
994  PetscReal final_resid=0.;
995 
996  Mat submat = nullptr;
997  Vec subrhs = nullptr;
998  Vec subsolution = nullptr;
999  VecScatter scatter = nullptr;
1000 
1001  // Close the matrices and vectors in case this wasn't already done.
1002  solution->close ();
1003  rhs->close ();
1004 
1005  // Prepare the matrix.
1006  Mat mat;
1007  ierr = MatCreateShell(this->comm().get(),
1008  rhs_in.local_size(),
1009  solution_in.local_size(),
1010  rhs_in.size(),
1011  solution_in.size(),
1012  const_cast<void *>(static_cast<const void *>(&shell_matrix)),
1013  &mat);
1014  /* Note that the const_cast above is only necessary because PETSc
1015  does not accept a const void *. Inside the member function
1016  _petsc_shell_matrix() below, the pointer is casted back to a
1017  const ShellMatrix<T> *. */
1018 
1019  LIBMESH_CHKERR(ierr);
1020  ierr = MatShellSetOperation(mat,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
1021  LIBMESH_CHKERR(ierr);
1022  ierr = MatShellSetOperation(mat,MATOP_MULT_ADD,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult_add));
1023  LIBMESH_CHKERR(ierr);
1024  ierr = MatShellSetOperation(mat,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
1025  LIBMESH_CHKERR(ierr);
1026 
1027  // Restrict rhs and solution vectors and set operators. The input
1028  // matrix works as the preconditioning matrix.
1029  if (_restrict_solve_to_is != nullptr)
1030  {
1031  PetscInt is_local_size = this->_restrict_solve_to_is_local_size();
1032 
1033  ierr = VecCreate(this->comm().get(),&subrhs);
1034  LIBMESH_CHKERR(ierr);
1035  ierr = VecSetSizes(subrhs,is_local_size,PETSC_DECIDE);
1036  LIBMESH_CHKERR(ierr);
1037  ierr = VecSetFromOptions(subrhs);
1038  LIBMESH_CHKERR(ierr);
1039 
1040  ierr = VecCreate(this->comm().get(),&subsolution);
1041  LIBMESH_CHKERR(ierr);
1042  ierr = VecSetSizes(subsolution,is_local_size,PETSC_DECIDE);
1043  LIBMESH_CHKERR(ierr);
1044  ierr = VecSetFromOptions(subsolution);
1045  LIBMESH_CHKERR(ierr);
1046 
1047  ierr = LibMeshVecScatterCreate(rhs->vec(),_restrict_solve_to_is, subrhs, nullptr, &scatter);
1048  LIBMESH_CHKERR(ierr);
1049 
1050  ierr = VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
1051  LIBMESH_CHKERR(ierr);
1052  ierr = VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
1053  LIBMESH_CHKERR(ierr);
1054 
1055  ierr = VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
1056  LIBMESH_CHKERR(ierr);
1057  ierr = VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
1058  LIBMESH_CHKERR(ierr);
1059 
1060 #if !PETSC_VERSION_LESS_THAN(3,1,0)
1061  ierr = LibMeshCreateSubMatrix(mat,
1064  MAT_INITIAL_MATRIX,
1065  &submat);
1066  LIBMESH_CHKERR(ierr);
1067 #endif
1068 
1069  /* Since removing columns of the matrix changes the equation
1070  system, we will now change the right hand side to compensate
1071  for this. Note that this is not necessary if \p SUBSET_ZERO
1072  has been selected. */
1074  {
1075  _create_complement_is(rhs_in);
1076  PetscInt is_complement_local_size =
1077  cast_int<PetscInt>(rhs_in.local_size()-is_local_size);
1078 
1079  Vec subvec1 = nullptr;
1080  Mat submat1 = nullptr;
1081  VecScatter scatter1 = nullptr;
1082 
1083  ierr = VecCreate(this->comm().get(),&subvec1);
1084  LIBMESH_CHKERR(ierr);
1085  ierr = VecSetSizes(subvec1,is_complement_local_size,PETSC_DECIDE);
1086  LIBMESH_CHKERR(ierr);
1087  ierr = VecSetFromOptions(subvec1);
1088  LIBMESH_CHKERR(ierr);
1089 
1090  ierr = LibMeshVecScatterCreate(rhs->vec(),_restrict_solve_to_is_complement, subvec1, nullptr, &scatter1);
1091  LIBMESH_CHKERR(ierr);
1092 
1093  ierr = VecScatterBegin(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
1094  LIBMESH_CHKERR(ierr);
1095  ierr = VecScatterEnd(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
1096  LIBMESH_CHKERR(ierr);
1097 
1098  ierr = VecScale(subvec1,-1.0);
1099  LIBMESH_CHKERR(ierr);
1100 
1101 #if !PETSC_VERSION_LESS_THAN(3,1,0)
1102  ierr = LibMeshCreateSubMatrix(mat,
1105  MAT_INITIAL_MATRIX,
1106  &submat1);
1107  LIBMESH_CHKERR(ierr);
1108 #endif
1109 
1110  // The following lines would be correct, but don't work
1111  // correctly in PETSc up to 3.1.0-p5. See discussion in
1112  // petsc-users of Nov 9, 2010.
1113  //
1114  // ierr = MatMultAdd(submat1,subvec1,subrhs,subrhs);
1115  // LIBMESH_CHKERR(ierr);
1116  //
1117  // We workaround by using a temporary vector. Note that the
1118  // fix in PETsc 3.1.0-p6 uses a temporary vector internally,
1119  // so this is no effective performance loss.
1120  Vec subvec2 = nullptr;
1121  ierr = VecCreate(this->comm().get(),&subvec2);
1122  LIBMESH_CHKERR(ierr);
1123  ierr = VecSetSizes(subvec2,is_local_size,PETSC_DECIDE);
1124  LIBMESH_CHKERR(ierr);
1125  ierr = VecSetFromOptions(subvec2);
1126  LIBMESH_CHKERR(ierr);
1127  ierr = MatMult(submat1,subvec1,subvec2);
1128  LIBMESH_CHKERR(ierr);
1129  ierr = VecAXPY(subrhs,1.0,subvec2);
1130 
1131  ierr = LibMeshVecScatterDestroy(&scatter1);
1132  LIBMESH_CHKERR(ierr);
1133  ierr = LibMeshVecDestroy(&subvec1);
1134  LIBMESH_CHKERR(ierr);
1135  ierr = LibMeshMatDestroy(&submat1);
1136  LIBMESH_CHKERR(ierr);
1137  }
1138 #if PETSC_RELEASE_LESS_THAN(3,5,0)
1139  ierr = KSPSetOperators(_ksp, submat, submat,
1140  DIFFERENT_NONZERO_PATTERN);
1141 #else
1142  ierr = KSPSetOperators(_ksp, submat, submat);
1143 #endif
1144  LIBMESH_CHKERR(ierr);
1145  }
1146  else
1147  {
1148 #if PETSC_RELEASE_LESS_THAN(3,5,0)
1149  ierr = KSPSetOperators(_ksp, mat, mat,
1150  DIFFERENT_NONZERO_PATTERN);
1151 #else
1152  ierr = KSPSetOperators(_ksp, mat, mat);
1153 #endif
1154  LIBMESH_CHKERR(ierr);
1155  }
1156 
1157  // Set the tolerances for the iterative solver. Use the user-supplied
1158  // tolerance for the relative residual & leave the others at default values.
1159  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
1160  PETSC_DEFAULT, max_its);
1161  LIBMESH_CHKERR(ierr);
1162 
1163  // Allow command line options to override anything set programmatically.
1164  ierr = KSPSetFromOptions(_ksp);
1165  LIBMESH_CHKERR(ierr);
1166 
1167  // Solve the linear system
1168  if (_restrict_solve_to_is != nullptr)
1169  {
1170  ierr = KSPSolve (_ksp, subrhs, subsolution);
1171  LIBMESH_CHKERR(ierr);
1172  }
1173  else
1174  {
1175  ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
1176  LIBMESH_CHKERR(ierr);
1177  }
1178 
1179  // Get the number of iterations required for convergence
1180  ierr = KSPGetIterationNumber (_ksp, &its);
1181  LIBMESH_CHKERR(ierr);
1182 
1183  // Get the norm of the final residual to return to the user.
1184  ierr = KSPGetResidualNorm (_ksp, &final_resid);
1185  LIBMESH_CHKERR(ierr);
1186 
1187  if (_restrict_solve_to_is != nullptr)
1188  {
1189  switch(_subset_solve_mode)
1190  {
1191  case SUBSET_ZERO:
1192  ierr = VecZeroEntries(solution->vec());
1193  LIBMESH_CHKERR(ierr);
1194  break;
1195 
1196  case SUBSET_COPY_RHS:
1197  ierr = VecCopy(rhs->vec(),solution->vec());
1198  LIBMESH_CHKERR(ierr);
1199  break;
1200 
1201  case SUBSET_DONT_TOUCH:
1202  /* Nothing to do here. */
1203  break;
1204 
1205  default:
1206  libmesh_error_msg("Invalid subset solve mode = " << _subset_solve_mode);
1207  }
1208  ierr = VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1209  LIBMESH_CHKERR(ierr);
1210  ierr = VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1211  LIBMESH_CHKERR(ierr);
1212 
1213  ierr = LibMeshVecScatterDestroy(&scatter);
1214  LIBMESH_CHKERR(ierr);
1215 
1216  ierr = LibMeshVecDestroy(&subsolution);
1217  LIBMESH_CHKERR(ierr);
1218  ierr = LibMeshVecDestroy(&subrhs);
1219  LIBMESH_CHKERR(ierr);
1220  ierr = LibMeshMatDestroy(&submat);
1221  LIBMESH_CHKERR(ierr);
1222  }
1223 
1224  // Destroy the matrix.
1225  ierr = LibMeshMatDestroy(&mat);
1226  LIBMESH_CHKERR(ierr);
1227 
1228  // return the # of its. and the final residual norm.
1229  return std::make_pair(its, final_resid);
1230 }
processor_id_type size() const
Definition: communicator.h:175
virtual void init(const char *name=nullptr) override
void _create_complement_is(const NumericVector< T > &vec_in)
static PetscErrorCode _petsc_shell_matrix_mult_add(Mat mat, Vec arg, Vec add, Vec dest)
const Parallel::Communicator & comm() const
static PetscErrorCode _petsc_shell_matrix_mult(Mat mat, Vec arg, Vec dest)
PetscErrorCode ierr
PetscInt _restrict_solve_to_is_local_size() const
static PetscErrorCode _petsc_shell_matrix_get_diagonal(Mat mat, Vec dest)

◆ solve() [4/6]

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 330 of file linear_solver.h.

336 {
337  if (pc_mat)
338  return this->solve(mat, *pc_mat, sol, rhs, tol, n_iter);
339  else
340  return this->solve(mat, sol, rhs, tol, n_iter);
341 }
virtual std::pair< unsigned int, Real > solve(SparseMatrix< T > &, NumericVector< T > &, NumericVector< T > &, const double, const unsigned int)=0

◆ solve() [5/6]

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 
)
overridevirtual

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 1236 of file petsc_linear_solver.C.

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

1242 {
1243 
1244 #if PETSC_VERSION_LESS_THAN(3,1,0)
1245  if (_restrict_solve_to_is != nullptr)
1246  libmesh_error_msg("The current implementation of subset solves with " \
1247  << "shell matrices requires PETSc version 3.1 or above. " \
1248  << "Older PETSc version do not support automatic " \
1249  << "submatrix generation of shell matrices.");
1250 #endif
1251 
1252  LOG_SCOPE("solve()", "PetscLinearSolver");
1253 
1254  // Make sure the data passed in are really of Petsc types
1255  const PetscMatrix<T> * precond = cast_ptr<const PetscMatrix<T> *>(&precond_matrix);
1256  PetscVector<T> * solution = cast_ptr<PetscVector<T> *>(&solution_in);
1257  PetscVector<T> * rhs = cast_ptr<PetscVector<T> *>(&rhs_in);
1258 
1259  this->init ();
1260 
1261  PetscErrorCode ierr=0;
1262  PetscInt its=0, max_its = static_cast<PetscInt>(m_its);
1263  PetscReal final_resid=0.;
1264 
1265  Mat submat = nullptr;
1266  Mat subprecond = nullptr;
1267  Vec subrhs = nullptr;
1268  Vec subsolution = nullptr;
1269  VecScatter scatter = nullptr;
1270  std::unique_ptr<PetscMatrix<Number>> subprecond_matrix;
1271 
1272  // Close the matrices and vectors in case this wasn't already done.
1273  solution->close ();
1274  rhs->close ();
1275 
1276  // Prepare the matrix.
1277  Mat mat;
1278  ierr = MatCreateShell(this->comm().get(),
1279  rhs_in.local_size(),
1280  solution_in.local_size(),
1281  rhs_in.size(),
1282  solution_in.size(),
1283  const_cast<void *>(static_cast<const void *>(&shell_matrix)),
1284  &mat);
1285  /* Note that the const_cast above is only necessary because PETSc
1286  does not accept a const void *. Inside the member function
1287  _petsc_shell_matrix() below, the pointer is casted back to a
1288  const ShellMatrix<T> *. */
1289 
1290  LIBMESH_CHKERR(ierr);
1291  ierr = MatShellSetOperation(mat,MATOP_MULT,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult));
1292  LIBMESH_CHKERR(ierr);
1293  ierr = MatShellSetOperation(mat,MATOP_MULT_ADD,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_mult_add));
1294  LIBMESH_CHKERR(ierr);
1295  ierr = MatShellSetOperation(mat,MATOP_GET_DIAGONAL,reinterpret_cast<void(*)(void)>(_petsc_shell_matrix_get_diagonal));
1296  LIBMESH_CHKERR(ierr);
1297 
1298  // Restrict rhs and solution vectors and set operators. The input
1299  // matrix works as the preconditioning matrix.
1300  if (_restrict_solve_to_is != nullptr)
1301  {
1302  PetscInt is_local_size = this->_restrict_solve_to_is_local_size();
1303 
1304  ierr = VecCreate(this->comm().get(),&subrhs);
1305  LIBMESH_CHKERR(ierr);
1306  ierr = VecSetSizes(subrhs,is_local_size,PETSC_DECIDE);
1307  LIBMESH_CHKERR(ierr);
1308  ierr = VecSetFromOptions(subrhs);
1309  LIBMESH_CHKERR(ierr);
1310 
1311  ierr = VecCreate(this->comm().get(),&subsolution);
1312  LIBMESH_CHKERR(ierr);
1313  ierr = VecSetSizes(subsolution,is_local_size,PETSC_DECIDE);
1314  LIBMESH_CHKERR(ierr);
1315  ierr = VecSetFromOptions(subsolution);
1316  LIBMESH_CHKERR(ierr);
1317 
1318  ierr = LibMeshVecScatterCreate(rhs->vec(),_restrict_solve_to_is, subrhs, nullptr, &scatter);
1319  LIBMESH_CHKERR(ierr);
1320 
1321  ierr = VecScatterBegin(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
1322  LIBMESH_CHKERR(ierr);
1323  ierr = VecScatterEnd(scatter,rhs->vec(),subrhs,INSERT_VALUES,SCATTER_FORWARD);
1324  LIBMESH_CHKERR(ierr);
1325 
1326  ierr = VecScatterBegin(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
1327  LIBMESH_CHKERR(ierr);
1328  ierr = VecScatterEnd(scatter,solution->vec(),subsolution,INSERT_VALUES,SCATTER_FORWARD);
1329  LIBMESH_CHKERR(ierr);
1330 
1331 #if !PETSC_VERSION_LESS_THAN(3,1,0)
1332  ierr = LibMeshCreateSubMatrix(mat,
1335  MAT_INITIAL_MATRIX,
1336  &submat);
1337  LIBMESH_CHKERR(ierr);
1338 
1339  ierr = LibMeshCreateSubMatrix(const_cast<PetscMatrix<T> *>(precond)->mat(),
1342  MAT_INITIAL_MATRIX,
1343  &subprecond);
1344  LIBMESH_CHKERR(ierr);
1345 #endif
1346 
1347  /* Since removing columns of the matrix changes the equation
1348  system, we will now change the right hand side to compensate
1349  for this. Note that this is not necessary if \p SUBSET_ZERO
1350  has been selected. */
1352  {
1353  _create_complement_is(rhs_in);
1354  PetscInt is_complement_local_size = rhs_in.local_size()-is_local_size;
1355 
1356  Vec subvec1 = nullptr;
1357  Mat submat1 = nullptr;
1358  VecScatter scatter1 = nullptr;
1359 
1360  ierr = VecCreate(this->comm().get(),&subvec1);
1361  LIBMESH_CHKERR(ierr);
1362  ierr = VecSetSizes(subvec1,is_complement_local_size,PETSC_DECIDE);
1363  LIBMESH_CHKERR(ierr);
1364  ierr = VecSetFromOptions(subvec1);
1365  LIBMESH_CHKERR(ierr);
1366 
1367  ierr = LibMeshVecScatterCreate(rhs->vec(),_restrict_solve_to_is_complement, subvec1, nullptr, &scatter1);
1368  LIBMESH_CHKERR(ierr);
1369 
1370  ierr = VecScatterBegin(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
1371  LIBMESH_CHKERR(ierr);
1372  ierr = VecScatterEnd(scatter1,_subset_solve_mode==SUBSET_COPY_RHS ? rhs->vec() : solution->vec(),subvec1,INSERT_VALUES,SCATTER_FORWARD);
1373  LIBMESH_CHKERR(ierr);
1374 
1375  ierr = VecScale(subvec1,-1.0);
1376  LIBMESH_CHKERR(ierr);
1377 
1378 #if !PETSC_VERSION_LESS_THAN(3,1,0)
1379  ierr = LibMeshCreateSubMatrix(mat,
1382  MAT_INITIAL_MATRIX,
1383  &submat1);
1384  LIBMESH_CHKERR(ierr);
1385 #endif
1386 
1387  // The following lines would be correct, but don't work
1388  // correctly in PETSc up to 3.1.0-p5. See discussion in
1389  // petsc-users of Nov 9, 2010.
1390  //
1391  // ierr = MatMultAdd(submat1,subvec1,subrhs,subrhs);
1392  // LIBMESH_CHKERR(ierr);
1393  //
1394  // We workaround by using a temporary vector. Note that the
1395  // fix in PETsc 3.1.0-p6 uses a temporary vector internally,
1396  // so this is no effective performance loss.
1397  Vec subvec2 = nullptr;
1398  ierr = VecCreate(this->comm().get(),&subvec2);
1399  LIBMESH_CHKERR(ierr);
1400  ierr = VecSetSizes(subvec2,is_local_size,PETSC_DECIDE);
1401  LIBMESH_CHKERR(ierr);
1402  ierr = VecSetFromOptions(subvec2);
1403  LIBMESH_CHKERR(ierr);
1404  ierr = MatMult(submat1,subvec1,subvec2);
1405  LIBMESH_CHKERR(ierr);
1406  ierr = VecAXPY(subrhs,1.0,subvec2);
1407  LIBMESH_CHKERR(ierr);
1408 
1409  ierr = LibMeshVecScatterDestroy(&scatter1);
1410  LIBMESH_CHKERR(ierr);
1411  ierr = LibMeshVecDestroy(&subvec1);
1412  LIBMESH_CHKERR(ierr);
1413  ierr = LibMeshMatDestroy(&submat1);
1414  LIBMESH_CHKERR(ierr);
1415  }
1416 
1417 #if PETSC_RELEASE_LESS_THAN(3,5,0)
1418  ierr = KSPSetOperators(_ksp, submat, subprecond,
1419  DIFFERENT_NONZERO_PATTERN);
1420 #else
1421  ierr = KSPSetOperators(_ksp, submat, subprecond);
1422 #endif
1423  LIBMESH_CHKERR(ierr);
1424 
1425  if (this->_preconditioner)
1426  {
1427  subprecond_matrix.reset(new PetscMatrix<Number>(subprecond, this->comm()));
1428  this->_preconditioner->set_matrix(*subprecond_matrix);
1429  this->_preconditioner->init();
1430  }
1431  }
1432  else
1433  {
1434 #if PETSC_RELEASE_LESS_THAN(3,5,0)
1435  ierr = KSPSetOperators(_ksp, mat, const_cast<PetscMatrix<T> *>(precond)->mat(),
1436  DIFFERENT_NONZERO_PATTERN);
1437 #else
1438  ierr = KSPSetOperators(_ksp, mat, const_cast<PetscMatrix<T> *>(precond)->mat());
1439 #endif
1440  LIBMESH_CHKERR(ierr);
1441 
1442  if (this->_preconditioner)
1443  {
1444  this->_preconditioner->set_matrix(const_cast<SparseMatrix<Number> &>(precond_matrix));
1445  this->_preconditioner->init();
1446  }
1447  }
1448 
1449  // Set the tolerances for the iterative solver. Use the user-supplied
1450  // tolerance for the relative residual & leave the others at default values.
1451  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
1452  PETSC_DEFAULT, max_its);
1453  LIBMESH_CHKERR(ierr);
1454 
1455  // Allow command line options to override anything set programmatically.
1456  ierr = KSPSetFromOptions(_ksp);
1457  LIBMESH_CHKERR(ierr);
1458 
1459  // Solve the linear system
1460  if (_restrict_solve_to_is != nullptr)
1461  {
1462  ierr = KSPSolve (_ksp, subrhs, subsolution);
1463  LIBMESH_CHKERR(ierr);
1464  }
1465  else
1466  {
1467  ierr = KSPSolve (_ksp, rhs->vec(), solution->vec());
1468  LIBMESH_CHKERR(ierr);
1469  }
1470 
1471  // Get the number of iterations required for convergence
1472  ierr = KSPGetIterationNumber (_ksp, &its);
1473  LIBMESH_CHKERR(ierr);
1474 
1475  // Get the norm of the final residual to return to the user.
1476  ierr = KSPGetResidualNorm (_ksp, &final_resid);
1477  LIBMESH_CHKERR(ierr);
1478 
1479  if (_restrict_solve_to_is != nullptr)
1480  {
1481  switch(_subset_solve_mode)
1482  {
1483  case SUBSET_ZERO:
1484  ierr = VecZeroEntries(solution->vec());
1485  LIBMESH_CHKERR(ierr);
1486  break;
1487 
1488  case SUBSET_COPY_RHS:
1489  ierr = VecCopy(rhs->vec(),solution->vec());
1490  LIBMESH_CHKERR(ierr);
1491  break;
1492 
1493  case SUBSET_DONT_TOUCH:
1494  /* Nothing to do here. */
1495  break;
1496 
1497  default:
1498  libmesh_error_msg("Invalid subset solve mode = " << _subset_solve_mode);
1499  }
1500  ierr = VecScatterBegin(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1501  LIBMESH_CHKERR(ierr);
1502  ierr = VecScatterEnd(scatter,subsolution,solution->vec(),INSERT_VALUES,SCATTER_REVERSE);
1503  LIBMESH_CHKERR(ierr);
1504 
1505  ierr = LibMeshVecScatterDestroy(&scatter);
1506  LIBMESH_CHKERR(ierr);
1507 
1508  if (this->_preconditioner)
1509  {
1510  // Before subprecond_matrix gets cleaned up, we should give
1511  // the _preconditioner a different matrix.
1512  this->_preconditioner->set_matrix(const_cast<SparseMatrix<Number> &>(precond_matrix));
1513  this->_preconditioner->init();
1514  }
1515 
1516  ierr = LibMeshVecDestroy(&subsolution);
1517  LIBMESH_CHKERR(ierr);
1518  ierr = LibMeshVecDestroy(&subrhs);
1519  LIBMESH_CHKERR(ierr);
1520  ierr = LibMeshMatDestroy(&submat);
1521  LIBMESH_CHKERR(ierr);
1522  ierr = LibMeshMatDestroy(&subprecond);
1523  LIBMESH_CHKERR(ierr);
1524  }
1525 
1526  // Destroy the matrix.
1527  ierr = LibMeshMatDestroy(&mat);
1528  LIBMESH_CHKERR(ierr);
1529 
1530  // return the # of its. and the final residual norm.
1531  return std::make_pair(its, final_resid);
1532 }
Preconditioner< T > * _preconditioner
processor_id_type size() const
Definition: communicator.h:175
virtual void init(const char *name=nullptr) override
void _create_complement_is(const NumericVector< T > &vec_in)
static PetscErrorCode _petsc_shell_matrix_mult_add(Mat mat, Vec arg, Vec add, Vec dest)
const Parallel::Communicator & comm() const
static PetscErrorCode _petsc_shell_matrix_mult(Mat mat, Vec arg, Vec dest)
PetscErrorCode ierr
PetscInt _restrict_solve_to_is_local_size() const
static PetscErrorCode _petsc_shell_matrix_get_diagonal(Mat mat, Vec dest)

◆ solve() [6/6]

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 347 of file linear_solver.h.

353 {
354  if (pc_mat)
355  return this->solve(mat, *pc_mat, sol, rhs, tol, n_iter);
356  else
357  return this->solve(mat, sol, rhs, tol, n_iter);
358 }
virtual std::pair< unsigned int, Real > solve(SparseMatrix< T > &, NumericVector< T > &, NumericVector< T > &, const double, const unsigned int)=0

◆ solver_type()

template<typename T>
SolverType libMesh::LinearSolver< T >::solver_type ( ) const
inlineinherited
Returns
The type of solver to use.

Definition at line 122 of file linear_solver.h.

122 { return _solver_type; }

Member Data Documentation

◆ _communicator

◆ _counts

ReferenceCounter::Counts libMesh::ReferenceCounter::_counts
staticprotectedinherited

◆ _enable_print_counter

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 141 of file reference_counter.h.

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

◆ _is_initialized

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

Flag indicating if the data structures have been initialized.

Definition at line 287 of file linear_solver.h.

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

◆ _ksp

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

Krylov subspace context

Definition at line 293 of file petsc_linear_solver.h.

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

◆ _mutex

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

Mutual exclusion object to enable thread-safe reference counting.

Definition at line 135 of file reference_counter.h.

◆ _n_objects

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 130 of file reference_counter.h.

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

◆ _pc

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

Preconditioner context

Definition at line 288 of file petsc_linear_solver.h.

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

◆ _preconditioner

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

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

Definition at line 292 of file linear_solver.h.

◆ _preconditioner_type

template<typename T>
PreconditionerType libMesh::LinearSolver< T >::_preconditioner_type
protectedinherited

Enum stating with type of preconditioner to use.

Definition at line 282 of file linear_solver.h.

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

◆ _restrict_solve_to_is

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

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

Definition at line 299 of file petsc_linear_solver.h.

◆ _restrict_solve_to_is_complement

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 306 of file petsc_linear_solver.h.

◆ _solver_configuration

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 306 of file linear_solver.h.

◆ _solver_type

template<typename T>
SolverType libMesh::LinearSolver< T >::_solver_type
protectedinherited

◆ _subset_solve_mode

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 324 of file petsc_linear_solver.h.

◆ same_preconditioner

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 300 of file linear_solver.h.


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