libMesh
Public Member Functions | Static Public Member Functions | Protected Types | Protected Member Functions | Protected Attributes | Static Protected Attributes | Private Member Functions | Static Private Member Functions | Private Attributes | List of all members
libMesh::PetscLinearSolver< T > Class Template Reference

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

#include <petsc_linear_solver.h>

Inheritance diagram for libMesh::PetscLinearSolver< T >:
[legend]

Public Member Functions

 PetscLinearSolver (const libMesh::Parallel::Communicator &comm_in)
 Constructor. More...
 
 ~PetscLinearSolver ()
 Destructor. More...
 
virtual void clear () override
 Release all memory and clear data structures. More...
 
virtual void init (const char *name=nullptr) override
 Initialize data structures if not done so already. More...
 
void init (PetscMatrix< T > *matrix, const char *name=nullptr)
 Initialize data structures if not done so already plus much more. More...
 
virtual void init_names (const System &) override
 Apply names to the system to be solved. More...
 
virtual void restrict_solve_to (const std::vector< unsigned int > *const dofs, const SubsetSolveMode subset_solve_mode=SUBSET_ZERO) override
 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. More...
 
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
 Call the Petsc solver. More...
 
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
 Call the Petsc solver. More...
 
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
 This method allows you to call a linear solver while specifying the matrix to use as the (left) preconditioning matrix. More...
 
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
 This function solves a system whose matrix is a shell matrix. More...
 
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
 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. More...
 
PC pc ()
 
KSP ksp ()
 
void get_residual_history (std::vector< double > &hist)
 Fills the input vector with the sequence of residual norms from the latest iterative solve. More...
 
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)
 Sets the type of solver to use. More...
 
PreconditionerType preconditioner_type () const
 
void set_preconditioner_type (const PreconditionerType pct)
 Sets the type of preconditioner to use. More...
 
void attach_preconditioner (Preconditioner< T > *preconditioner)
 Attaches a Preconditioner object to be used. More...
 
virtual void reuse_preconditioner (bool)
 Set the same_preconditioner flag, which indicates if we reuse the same preconditioner for subsequent solves. More...
 
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)
 This function calls the solver "_solver_type" preconditioned with the "_preconditioner_type" preconditioner. More...
 
std::pair< unsigned int, Realsolve (const ShellMatrix< T > &matrix, const SparseMatrix< T > *precond_matrix, NumericVector< T > &, NumericVector< T > &, const double, const unsigned int)
 This function solves a system whose matrix is a shell matrix, but an optional sparse matrix may be used as preconditioning matrix. More...
 
virtual void print_converged_reason () const
 Prints a useful message about why the latest linear solve con(di)verged. More...
 
void set_solver_configuration (SolverConfiguration &solver_configuration)
 Set the solver configuration object. More...
 
const Parallel::Communicator & comm () 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())
 Builds a LinearSolver using the linear solver package specified by solver_package. More...
 
static std::string get_info ()
 Gets a string containing the reference information. More...
 
static void print_info (std::ostream &out=libMesh::out)
 Prints the reference information, by default to libMesh::out. More...
 
static unsigned int n_objects ()
 Prints the number of outstanding (created, but not yet destroyed) objects. More...
 
static void enable_print_counter_info ()
 Methods to enable/disable the reference counter output from print_info() More...
 
static void disable_print_counter_info ()
 

Protected Types

typedef std::map< std::string, std::pair< unsigned int, unsigned int > > Counts
 Data structure to log the information. More...
 

Protected Member Functions

void increment_constructor_count (const std::string &name)
 Increments the construction counter. More...
 
void increment_destructor_count (const std::string &name)
 Increments the destruction counter. More...
 

Protected Attributes

SolverType _solver_type
 Enum stating which type of iterative solver to use. More...
 
PreconditionerType _preconditioner_type
 Enum stating with type of preconditioner to use. More...
 
bool _is_initialized
 Flag indicating if the data structures have been initialized. More...
 
Preconditioner< T > * _preconditioner
 Holds the Preconditioner object to be used for the linear solves. More...
 
bool same_preconditioner
 Boolean flag to indicate whether we want to use an identical preconditioner to the previous solve. More...
 
SolverConfiguration_solver_configuration
 Optionally store a SolverOptions object that can be used to set parameters like solver type, tolerances and iteration limits. More...
 
const Parallel::Communicator & _communicator
 

Static Protected Attributes

static Counts _counts
 Actually holds the data. More...
 
static Threads::atomic< unsigned int_n_objects
 The number of objects. More...
 
static Threads::spin_mutex _mutex
 Mutual exclusion object to enable thread-safe reference counting. More...
 
static bool _enable_print_counter = true
 Flag to control whether reference count information is printed when print_info is called. More...
 

Private Member Functions

void set_petsc_solver_type ()
 Tells PETSC to use the user-specified solver stored in _solver_type. More...
 
PetscInt _restrict_solve_to_is_local_size () const
 
void _create_complement_is (const NumericVector< T > &vec_in)
 Creates _restrict_solve_to_is_complement to contain all indices that are local in vec_in, except those that are contained in _restrict_solve_to_is. More...
 

Static Private Member Functions

static PetscErrorCode _petsc_shell_matrix_mult (Mat mat, Vec arg, Vec dest)
 Internal function if shell matrix mode is used. More...
 
static PetscErrorCode _petsc_shell_matrix_mult_add (Mat mat, Vec arg, Vec add, Vec dest)
 Internal function if shell matrix mode is used. More...
 
static PetscErrorCode _petsc_shell_matrix_get_diagonal (Mat mat, Vec dest)
 Internal function if shell matrix mode is used. More...
 

Private Attributes

PC _pc
 Preconditioner context. More...
 
KSP _ksp
 Krylov subspace context. More...
 
IS _restrict_solve_to_is
 PETSc index set containing the dofs on which to solve (nullptr means solve on all dofs). More...
 
IS _restrict_solve_to_is_complement
 PETSc index set, complement to _restrict_solve_to_is. More...
 
SubsetSolveMode _subset_solve_mode
 If restrict-solve-to-subset mode is active, this member decides what happens with the dofs outside the subset. More...
 

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 107 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()

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

Constructor.

Initializes Petsc data structures

Definition at line 135 of file petsc_linear_solver.C.

135  :
136  LinearSolver<T>(comm_in),
137  _restrict_solve_to_is(nullptr),
140 {
141  if (this->n_processors() == 1)
143  else
145 }

◆ ~PetscLinearSolver()

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

Destructor.

Definition at line 337 of file petsc_linear_solver.h.

338 {
339  this->clear ();
340 }

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

368 {
370 #if PETSC_VERSION_LESS_THAN(3,0,0)
371  // No ISComplement in PETSc 2.3.3
372  libmesh_not_implemented();
373 #else
375  {
376  int ierr = ISComplement(_restrict_solve_to_is,
377  vec_in.first_local_index(),
378  vec_in.last_local_index(),
380  LIBMESH_CHKERR(ierr);
381  }
382 #endif
383 }

◆ _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 1788 of file petsc_linear_solver.C.

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

◆ _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 1734 of file petsc_linear_solver.C.

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

◆ _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 1758 of file petsc_linear_solver.C.

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

◆ _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 346 of file petsc_linear_solver.h.

347 {
349 
350  PetscInt s;
351  int ierr = ISGetLocalSize(_restrict_solve_to_is, &s);
352  LIBMESH_CHKERR(ierr);
353 
354  return s;
355 }

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

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

◆ 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 }

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

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 }

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

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

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

◆ comm()

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

Definition at line 94 of file parallel_object.h.

95  { return _communicator; }

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::PetscLinearSolver< Number >::_petsc_shell_matrix_get_diagonal(), libMesh::SlepcEigenSolver< libMesh::Number >::_petsc_shell_matrix_get_diagonal(), libMesh::PetscLinearSolver< Number >::_petsc_shell_matrix_mult(), libMesh::SlepcEigenSolver< libMesh::Number >::_petsc_shell_matrix_mult(), libMesh::PetscLinearSolver< Number >::_petsc_shell_matrix_mult_add(), libMesh::EquationSystems::_read_impl(), libMesh::MeshRefinement::_refine_elements(), libMesh::MeshRefinement::_smooth_flags(), libMesh::DofMap::add_constraints_to_send_list(), add_cube_convex_hull_to_mesh(), libMesh::PetscDMWrapper::add_dofs_helper(), libMesh::PetscDMWrapper::add_dofs_to_section(), libMesh::TransientRBConstruction::add_IC_to_RB_space(), libMesh::ImplicitSystem::add_matrix(), libMesh::RBConstruction::add_scaled_matrix_and_vector(), libMesh::DynaIO::add_spline_constraints(), libMesh::System::add_vector(), libMesh::UnstructuredMesh::all_second_order(), libMesh::MeshTools::Modification::all_tri(), libMesh::LaplaceMeshSmoother::allgather_graph(), libMesh::DofMap::allgather_recursive_constraints(), libMesh::TransientRBConstruction::allocate_data_structures(), libMesh::RBConstruction::allocate_data_structures(), libMesh::TransientRBConstruction::assemble_affine_expansion(), 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::System::calculate_norm(), libMesh::DofMap::check_dirichlet_bcid_consistency(), libMesh::RBConstruction::compute_Fq_representor_innerprods(), libMesh::RBConstruction::compute_max_error_bound(), 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::RBConstruction::compute_output_dual_innerprods(), libMesh::RBConstruction::compute_residual_dual_norm_slow(), libMesh::RBSCMConstruction::compute_SCM_bounds_on_training_set(), libMesh::Problem_Interface::computeF(), libMesh::Problem_Interface::computeJacobian(), libMesh::Problem_Interface::computePreconditioner(), libMesh::ExodusII_IO::copy_elemental_solution(), libMesh::ExodusII_IO::copy_scalar_solution(), libMesh::MeshTools::correct_node_proc_ids(), libMesh::MeshTools::create_bounding_box(), libMesh::DofMap::create_dof_constraints(), 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::DTKSolutionTransfer::DTKSolutionTransfer(), libMesh::MeshRefinement::eliminate_unrefined_patches(), libMesh::RBEIMConstruction::enrich_RB_space(), libMesh::TransientRBConstruction::enrich_RB_space(), libMesh::RBConstruction::enrich_RB_space(), libMesh::EpetraVector< T >::EpetraVector(), AssembleOptimization::equality_constraints(), libMesh::WeightedPatchRecoveryErrorEstimator::estimate_error(), libMesh::PatchRecoveryErrorEstimator::estimate_error(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::ExactErrorEstimator::estimate_error(), libMesh::RBEIMConstruction::evaluate_mesh_function(), libMesh::MeshRefinement::flag_elements_by_elem_fraction(), libMesh::MeshRefinement::flag_elements_by_error_fraction(), libMesh::MeshRefinement::flag_elements_by_error_tolerance(), libMesh::MeshRefinement::flag_elements_by_mean_stddev(), libMesh::MeshRefinement::flag_elements_by_nelem_target(), libMesh::DofMap::gather_constraints(), libMesh::MeshfreeInterpolation::gather_remote_data(), libMesh::CondensedEigenSystem::get_eigenpair(), libMesh::DofMap::get_info(), libMesh::ImplicitSystem::get_linear_solver(), AssembleOptimization::inequality_constraints(), AssembleOptimization::inequality_constraints_jacobian(), libMesh::LocationMap< T >::init(), libMesh::TimeSolver::init(), libMesh::SystemSubsetBySubdomain::init(), libMesh::PetscDMWrapper::init_and_attach_petscdm(), libMesh::EigenSystem::init_matrices(), libMesh::OptimizationSystem::initialize_equality_constraints_storage(), libMesh::OptimizationSystem::initialize_inequality_constraints_storage(), libMesh::RBEIMConstruction::initialize_rb_construction(), integrate_function(), 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_preconditioner_apply(), 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(), main(), 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::TransientRBConstruction::mass_matrix_scaled_matvec(), libMesh::FEMSystem::mesh_position_set(), LinearElasticityWithContact::move_mesh(), libMesh::DistributedMesh::n_active_elem(), libMesh::MeshTools::n_active_levels(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::DofMap::n_constrained_dofs(), 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::DofMap::print_dof_constraints(), FEMParameters::read(), libMesh::Nemesis_IO::read(), libMesh::XdrIO::read(), libMesh::CheckpointIO::read_header(), libMesh::XdrIO::read_header(), libMesh::System::read_header(), libMesh::RBEvaluation::read_in_vectors_from_multiple_files(), libMesh::System::read_legacy_data(), libMesh::TransientRBConstruction::read_riesz_representors_from_files(), libMesh::RBConstruction::read_riesz_representors_from_files(), 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::MeshRefinement::refine_and_coarsen_elements(), libMesh::DistributedMesh::renumber_dof_objects(), LinearElasticityWithContact::residual_and_jacobian(), OverlappingAlgebraicGhostingTest::run_ghosting_test(), OverlappingCouplingGhostingTest::run_sparsity_pattern_test(), libMesh::DofMap::scatter_constraints(), libMesh::CheckpointIO::select_split_config(), libMesh::TransientRBConstruction::set_error_temporal_data(), libMesh::RBEIMConstruction::set_explicit_sys_subvector(), 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::BoundaryInfo::sync(), libMesh::MeshRefinement::test_level_one(), MeshfunctionDFEM::test_mesh_function_dfem(), MeshfunctionDFEM::test_mesh_function_dfem_grad(), MeshFunctionTest::test_p_level(), libMesh::MeshRefinement::test_unflagged(), SystemsTest::testBlockRestrictedVarNDofs(), PointLocatorTest::testLocator(), BoundaryInfoTest::testMesh(), SystemsTest::testProjectCubeWithMeshFunction(), CheckpointIOTest::testSplitter(), libMesh::MeshTools::total_weight(), libMesh::MeshFunctionSolutionTransfer::transfer(), libMesh::MeshfreeSolutionTransfer::transfer(), libMesh::TransientRBConstruction::truth_assembly(), libMesh::RBConstruction::truth_assembly(), libMesh::MeshRefinement::uniformly_coarsen(), libMesh::TransientRBConstruction::update_RB_initial_condition_all_N(), libMesh::RBEIMConstruction::update_RB_system_matrices(), libMesh::TransientRBConstruction::update_RB_system_matrices(), libMesh::RBConstruction::update_RB_system_matrices(), libMesh::TransientRBConstruction::update_residual_terms(), libMesh::RBConstruction::update_residual_terms(), libMesh::NameBasedIO::write(), libMesh::XdrIO::write(), libMesh::VTKIO::write_nodal_data(), libMesh::RBEvaluation::write_out_vectors(), libMesh::TransientRBConstruction::write_riesz_representors_to_files(), libMesh::RBConstruction::write_riesz_representors_to_files(), 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(), libMesh::XdrIO::write_serialized_nodesets(), libMesh::RBDataSerialization::RBEvaluationSerialization::write_to_file(), libMesh::RBDataSerialization::TransientRBEvaluationSerialization::write_to_file(), libMesh::RBDataSerialization::RBEIMEvaluationSerialization::write_to_file(), and libMesh::RBDataSerialization::RBSCMEvaluationSerialization::write_to_file().

◆ disable_print_counter_info()

void libMesh::ReferenceCounter::disable_print_counter_info ( )
staticinherited

Definition at line 106 of file reference_counter.C.

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

References libMesh::ReferenceCounter::_enable_print_counter.

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

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

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

References libMesh::ReferenceCounter::_enable_print_counter.

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

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

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

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 }

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

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

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

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

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

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

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

322 {
323  return same_preconditioner;
324 }

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

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 }

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

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

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

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 }

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

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

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

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

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

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

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

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

392 {
393  petsc_auto_fieldsplit(this->pc(), sys);
394 }

◆ initialized()

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

Definition at line 95 of file linear_solver.h.

95 { 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 246 of file petsc_linear_solver.h.

246 { this->init(); return _ksp; }

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

84  { return _n_objects; }

References libMesh::ReferenceCounter::_n_objects.

◆ n_processors()

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

Definition at line 100 of file parallel_object.h.

101  { return cast_int<processor_id_type>(_communicator.size()); }

References libMesh::ParallelObject::_communicator.

Referenced by libMesh::BoundaryInfo::_find_id_maps(), libMesh::DofMap::add_constraints_to_send_list(), libMesh::PetscDMWrapper::add_dofs_to_section(), libMesh::DistributedMesh::add_elem(), libMesh::DofMap::add_neighbors_to_send_list(), libMesh::DistributedMesh::add_node(), libMesh::LaplaceMeshSmoother::allgather_graph(), libMesh::DofMap::allgather_recursive_constraints(), 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::ExodusII_IO::copy_scalar_solution(), 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::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::partition(), libMesh::System::point_gradient(), libMesh::System::point_hessian(), libMesh::System::point_value(), libMesh::DofMap::prepare_send_list(), libMesh::DofMap::print_dof_constraints(), 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(), OverlappingFunctorTest::run_partitioner_test(), libMesh::DofMap::scatter_constraints(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::PetscDMWrapper::set_point_range_in_section(), CheckpointIOTest::testSplitter(), WriteVecAndScalar::testWrite(), libMesh::MeshRefinement::uniformly_coarsen(), libMesh::DistributedMesh::update_parallel_id_counts(), libMesh::GMVIO::write_binary(), libMesh::GMVIO::write_discontinuous_gmv(), libMesh::VTKIO::write_nodal_data(), 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().

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

238 { this->init(); return _pc; }

◆ 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 }

◆ 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 }

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

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

References libMesh::ReferenceCounter::_enable_print_counter, and libMesh::ReferenceCounter::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 106 of file parallel_object.h.

107  { return cast_int<processor_id_type>(_communicator.rank()); }

References libMesh::ParallelObject::_communicator.

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::DofMap::allgather_recursive_constraints(), libMesh::FEMSystem::assembly(), 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::RBConstruction::compute_max_error_bound(), 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::ExodusII_IO::copy_scalar_solution(), 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::RBEIMConstruction::evaluate_mesh_function(), 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::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::DofMap::get_local_constraints(), 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::PetscDMWrapper::init_and_attach_petscdm(), HeatSystem::init_data(), 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::TransientRBEvaluation::legacy_write_offline_data_to_files(), libMesh::RBEIMEvaluation::legacy_write_offline_data_to_files(), libMesh::RBEvaluation::legacy_write_offline_data_to_files(), libMesh::RBSCMEvaluation::legacy_write_offline_data_to_files(), libMesh::RBEIMEvaluation::legacy_write_out_interpolation_points_elem(), 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(), main(), libMesh::MeshRefinement::make_coarsening_compatible(), AugmentSparsityOnInterface::mesh_reinit(), 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::DofMap::print_dof_constraints(), 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::RBEvaluation::read_in_vectors_from_multiple_files(), libMesh::System::read_legacy_data(), libMesh::ExodusII_IO_Helper::read_node_num_map(), libMesh::System::read_parallel_data(), libMesh::TransientRBConstruction::read_riesz_representors_from_files(), libMesh::RBConstruction::read_riesz_representors_from_files(), 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::DofMap::scatter_constraints(), libMesh::CheckpointIO::select_split_config(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::PetscDMWrapper::set_point_range_in_section(), libMesh::LaplaceMeshSmoother::smooth(), DefaultCouplingTest::testCoupling(), PointNeighborCouplingTest::testCoupling(), MeshInputTest::testDynaReadElem(), MeshInputTest::testDynaReadPatch(), MeshInputTest::testExodusCopyElementSolution(), MeshInputTest::testExodusWriteElementDataFromDiscontinuousNodalData(), SystemsTest::testProjectMatrix1D(), SystemsTest::testProjectMatrix2D(), SystemsTest::testProjectMatrix3D(), BoundaryInfoTest::testShellFaceConstraints(), CheckpointIOTest::testSplitter(), WriteVecAndScalar::testWrite(), libMesh::MeshTools::total_weight(), libMesh::MeshRefinement::uniformly_coarsen(), libMesh::Parallel::Packing< Node * >::unpack(), libMesh::Parallel::Packing< Elem * >::unpack(), libMesh::DistributedMesh::update_parallel_id_counts(), libMesh::DTKAdapter::update_variable_values(), 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_element_values_element_major(), 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::VTKIO::write_nodal_data(), 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::RBEvaluation::write_out_vectors(), write_output_solvedata(), libMesh::System::write_parallel_data(), libMesh::RBConstruction::write_riesz_representors_to_files(), 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::ExodusII_IO_Helper::write_sideset_data(), 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().

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

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

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

130 {
131  same_preconditioner = reuse_flag;
132 }

Referenced by libMesh::ImplicitSystem::disable_cache(), libMesh::RBConstruction::initialize_rb_construction(), and main().

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

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

◆ 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 }

Referenced by TimeSolverTestImplementation< NewmarkSolver >::run_test_with_exact_soln(), and SystemsTest::testDofCouplingWithVarGroups().

◆ 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 }

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

127  { _solver_type = st; }

Referenced by main(), TimeSolverTestImplementation< NewmarkSolver >::run_test_with_exact_soln(), and SystemsTest::testDofCouplingWithVarGroups().

◆ solve() [1/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 346 of file linear_solver.h.

352 {
353  if (pc_mat)
354  return this->solve(mat, *pc_mat, sol, rhs, tol, n_iter);
355  else
356  return this->solve(mat, sol, rhs, tol, n_iter);
357 }

◆ solve() [2/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 1237 of file petsc_linear_solver.C.

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

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

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

◆ solve() [4/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 434 of file petsc_linear_solver.C.

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

◆ solve() [5/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 329 of file linear_solver.h.

335 {
336  if (pc_mat)
337  return this->solve(mat, *pc_mat, sol, rhs, tol, n_iter);
338  else
339  return this->solve(mat, sol, rhs, tol, n_iter);
340 }

◆ solve() [6/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 163 of file petsc_linear_solver.h.

168  {
169  return this->solve(matrix_in, matrix_in, solution_in, rhs_in, tol, m_its);
170  }

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

◆ solver_type()

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

Definition at line 121 of file linear_solver.h.

121 { return _solver_type; }

Member Data Documentation

◆ _communicator

const Parallel::Communicator& libMesh::ParallelObject::_communicator
protectedinherited

◆ _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 286 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 299 of file petsc_linear_solver.h.

Referenced by libMesh::PetscLinearSolver< Number >::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 294 of file petsc_linear_solver.h.

Referenced by libMesh::PetscLinearSolver< Number >::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 291 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 281 of file linear_solver.h.

Referenced by libMesh::AztecLinearSolver< T >::AztecLinearSolver().

◆ _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 305 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 312 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 305 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 330 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 299 of file linear_solver.h.


The documentation for this class was generated from the following files:
libMesh::DIVERGED_BREAKDOWN_BICG
Definition: enum_convergence_flags.h:49
libMesh::PetscLinearSolver::_restrict_solve_to_is_local_size
PetscInt _restrict_solve_to_is_local_size() const
Definition: petsc_linear_solver.h:346
libMesh::CONVERGED_STEP_LENGTH
Definition: enum_convergence_flags.h:42
libMesh::LinearSolver::get_converged_reason
virtual LinearConvergenceReason get_converged_reason() const =0
libMesh::UNKNOWN_FLAG
Definition: enum_convergence_flags.h:58
libMesh::PetscLinearSolver::_petsc_shell_matrix_get_diagonal
static PetscErrorCode _petsc_shell_matrix_get_diagonal(Mat mat, Vec dest)
Internal function if shell matrix mode is used.
Definition: petsc_linear_solver.C:1788
libMesh::DIVERGED_INDEFINITE_MAT
Definition: enum_convergence_flags.h:53
libMesh::PETSC_SOLVERS
Definition: enum_solver_package.h:36
libMesh::PetscLinearSolver::_petsc_shell_matrix_mult_add
static PetscErrorCode _petsc_shell_matrix_mult_add(Mat mat, Vec arg, Vec add, Vec dest)
Internal function if shell matrix mode is used.
Definition: petsc_linear_solver.C:1758
libMesh::BICGSTAB
Definition: enum_solver_type.h:42
libMesh::DIVERGED_BREAKDOWN
Definition: enum_convergence_flags.h:48
libMesh::DIVERGED_DTOL
Definition: enum_convergence_flags.h:47
libMesh::LinearSolver::solve
virtual std::pair< unsigned int, Real > solve(SparseMatrix< T > &, NumericVector< T > &, NumericVector< T > &, const double, const unsigned int)=0
This function calls the solver _solver_type preconditioned with the _preconditioner_type precondition...
libMesh::DIVERGED_ITS
Definition: enum_convergence_flags.h:46
libMesh::BICG
Definition: enum_solver_type.h:41
libMesh::CR
Definition: enum_solver_type.h:37
libMesh::index_range
IntRange< std::size_t > index_range(const std::vector< T > &vec)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:106
libMesh::LinearSolver::_solver_type
SolverType _solver_type
Enum stating which type of iterative solver to use.
Definition: linear_solver.h:276
libMesh::RICHARDSON
Definition: enum_solver_type.h:50
libMesh::ReferenceCounter::_counts
static Counts _counts
Actually holds the data.
Definition: reference_counter.h:122
libMesh::MINRES
Definition: enum_solver_type.h:43
libMesh::DIVERGED_NULL
Definition: enum_convergence_flags.h:45
libMesh::ParallelObject::comm
const Parallel::Communicator & comm() const
Definition: parallel_object.h:94
libMesh::ReferenceCounter::_n_objects
static Threads::atomic< unsigned int > _n_objects
The number of objects.
Definition: reference_counter.h:130
libMesh::CONVERGED_CG_CONSTRAINED
Definition: enum_convergence_flags.h:41
libMesh::SUBSET_COPY_RHS
Definition: enum_subset_solve_mode.h:40
libMesh::PetscLinearSolver::pc
PC pc()
Definition: petsc_linear_solver.h:238
libMesh::ierr
ierr
Definition: petsc_dm_wrapper.C:72
libMesh::CONVERGED_ITERATING
Definition: enum_convergence_flags.h:56
libMesh::ReferenceCounter::get_info
static std::string get_info()
Gets a string containing the reference information.
Definition: reference_counter.C:47
libMesh::DIVERGED_NONSYMMETRIC
Definition: enum_convergence_flags.h:50
libMesh::LinearSolver::_preconditioner_type
PreconditionerType _preconditioner_type
Enum stating with type of preconditioner to use.
Definition: linear_solver.h:281
libMesh::SUBSET_ZERO
Definition: enum_subset_solve_mode.h:37
PetscBool
PetscTruth PetscBool
Definition: petsc_macro.h:73
libMesh::PetscLinearSolver::_ksp
KSP _ksp
Krylov subspace context.
Definition: petsc_linear_solver.h:299
libMesh::CONVERGED_RTOL
Definition: enum_convergence_flags.h:37
libMesh::PetscLinearSolver::_restrict_solve_to_is
IS _restrict_solve_to_is
PETSc index set containing the dofs on which to solve (nullptr means solve on all dofs).
Definition: petsc_linear_solver.h:305
libMesh::PetscLinearSolver::_petsc_shell_matrix_mult
static PetscErrorCode _petsc_shell_matrix_mult(Mat mat, Vec arg, Vec dest)
Internal function if shell matrix mode is used.
Definition: petsc_linear_solver.C:1734
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::CG
Definition: enum_solver_type.h:34
libMesh::LinearConvergenceReason
LinearConvergenceReason
Linear solver convergence flags (taken from the PETSc flags).
Definition: enum_convergence_flags.h:33
libMesh::TFQMR
Definition: enum_solver_type.h:40
libMesh::libmesh_petsc_preconditioner_setup
PetscErrorCode libmesh_petsc_preconditioner_setup(void *ctx)
This function is called by PETSc to initialize the preconditioner.
Definition: petsc_linear_solver.C:49
libMesh::CGS
Definition: enum_solver_type.h:36
libMesh::ParallelObject::n_processors
processor_id_type n_processors() const
Definition: parallel_object.h:100
libMesh::PetscLinearSolver::set_petsc_solver_type
void set_petsc_solver_type()
Tells PETSC to use the user-specified solver stored in _solver_type.
Definition: petsc_linear_solver.C:1601
libMesh::PetscLinearSolver::clear
virtual void clear() override
Release all memory and clear data structures.
Definition: petsc_linear_solver.C:150
libMesh::Threads::spin_mtx
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
Definition: threads.C:29
libMesh::libmesh_ignore
void libmesh_ignore(const Args &...)
Definition: libmesh_common.h:526
libMesh::LinearSolver::_preconditioner
Preconditioner< T > * _preconditioner
Holds the Preconditioner object to be used for the linear solves.
Definition: linear_solver.h:291
libMesh::ParallelObject::_communicator
const Parallel::Communicator & _communicator
Definition: parallel_object.h:112
libMesh::TCQMR
Definition: enum_solver_type.h:39
libMesh::BLOCK_JACOBI_PRECOND
Definition: enum_preconditioner_type.h:36
libMesh::PetscLinearSolver::_subset_solve_mode
SubsetSolveMode _subset_solve_mode
If restrict-solve-to-subset mode is active, this member decides what happens with the dofs outside th...
Definition: petsc_linear_solver.h:330
libMesh::PetscLinearSolver::solve
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
Call the Petsc solver.
Definition: petsc_linear_solver.h:163
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh::LASPACK_SOLVERS
Definition: enum_solver_package.h:38
libMesh::CONVERGED_HAPPY_BREAKDOWN
Definition: enum_convergence_flags.h:43
libMesh::GMRES
Definition: enum_solver_type.h:44
libMesh::LinearSolver::_is_initialized
bool _is_initialized
Flag indicating if the data structures have been initialized.
Definition: linear_solver.h:286
libMesh::CONVERGED_RTOL_NORMAL
Definition: enum_convergence_flags.h:35
libMesh::DIVERGED_INDEFINITE_PC
Definition: enum_convergence_flags.h:51
PETSC_OWN_POINTER
Definition: petsc_macro.h:103
libMesh::TRILINOS_SOLVERS
Definition: enum_solver_package.h:37
libMesh::CONVERGED_CG_NEG_CURVE
Definition: enum_convergence_flags.h:40
libMesh::SolverConfiguration::set_options_during_init
virtual void set_options_during_init()
Apply options during initialization of a solver.
Definition: solver_configuration.h:58
libMesh::LinearSolver::same_preconditioner
bool same_preconditioner
Boolean flag to indicate whether we want to use an identical preconditioner to the previous solve.
Definition: linear_solver.h:299
libMesh::CONVERGED_ATOL_NORMAL
Definition: enum_convergence_flags.h:36
libMesh::ReferenceElem::get
const Elem & get(const ElemType type_in)
Definition: reference_elem.C:237
libMesh::PetscLinearSolver::init
virtual void init(const char *name=nullptr) override
Initialize data structures if not done so already.
Definition: petsc_linear_solver.C:193
libMesh::SUBSET_DONT_TOUCH
Definition: enum_subset_solve_mode.h:45
libMesh::CONVERGED_ITS
Definition: enum_convergence_flags.h:39
libMesh::SHELL_PRECOND
Definition: enum_preconditioner_type.h:46
libMesh::LinearSolver::initialized
bool initialized() const
Definition: linear_solver.h:95
libMesh::CONVERGED_ATOL
Definition: enum_convergence_flags.h:38
libMesh::LSQR
Definition: enum_solver_type.h:45
libMesh::PetscLinearSolver::_pc
PC _pc
Preconditioner context.
Definition: petsc_linear_solver.h:294
libMesh::DIVERGED_PCSETUP_FAILED
Definition: enum_convergence_flags.h:54
libMesh::LinearSolver::_solver_configuration
SolverConfiguration * _solver_configuration
Optionally store a SolverOptions object that can be used to set parameters like solver type,...
Definition: linear_solver.h:305
libMesh::DIVERGED_NAN
Definition: enum_convergence_flags.h:52
libMesh::ctx
void * ctx
Definition: petsc_dm_wrapper.C:71
libMesh::libmesh_petsc_preconditioner_apply
PetscErrorCode libmesh_petsc_preconditioner_apply(void *ctx, Vec x, Vec y)
This function is called by PETSc to actually apply the preconditioner.
Definition: petsc_linear_solver.C:62
libMesh::PetscLinearSolver::_create_complement_is
void _create_complement_is(const NumericVector< T > &vec_in)
Creates _restrict_solve_to_is_complement to contain all indices that are local in vec_in,...
Definition: petsc_linear_solver.h:361
libMesh::PetscLinearSolver::_restrict_solve_to_is_complement
IS _restrict_solve_to_is_complement
PETSc index set, complement to _restrict_solve_to_is.
Definition: petsc_linear_solver.h:312
libMesh::err
OStreamProxy err
libMesh::CHEBYSHEV
Definition: enum_solver_type.h:51
libMesh::ReferenceCounter::_enable_print_counter
static bool _enable_print_counter
Flag to control whether reference count information is printed when print_info is called.
Definition: reference_counter.h:141
libMesh::petsc_auto_fieldsplit
void petsc_auto_fieldsplit(PC my_pc, const System &sys)
Definition: petsc_auto_fieldsplit.C:58
libMesh::PetscPreconditioner::set_petsc_preconditioner_type
static void set_petsc_preconditioner_type(const PreconditionerType &preconditioner_type, PC &pc)
Tells PETSc to use the user-specified preconditioner.
Definition: petsc_preconditioner.C:115
libMesh::out
OStreamProxy out
libMesh::EIGEN_SOLVERS
Definition: enum_solver_package.h:40
libMesh::Quality::name
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition: elem_quality.C:42
libMesh::SolverConfiguration::configure_solver
virtual void configure_solver()=0
Apply solver options to a particular solver.
libMesh::ILU_PRECOND
Definition: enum_preconditioner_type.h:43