petsc_linear_solver.c

一个用来实现偏微分方程中网格的计算库

// $Id: petsc_linear_solver.C 2789 2008-04-13 02:24:40Z roystgnr $

// The libMesh Finite Element Library.
// Copyright (C) 2002-2007  Benjamin S. Kirk, John W. Peterson
  
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
  
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



#include "libmesh_common.h"

#ifdef HAVE_PETSC


// C++ includes

// Local Includes
#include "libmesh_logging.h"
#include "petsc_linear_solver.h"



/*----------------------- functions ----------------------------------*/
template <typename T>
void PetscLinearSolver<T>::clear ()
{
  if (this->initialized())
    {
      this->_is_initialized = false;

      int ierr=0;

#if PETSC_VERSION_LESS_THAN(2,2,0)

  // 2.1.x & earlier style      
  ierr = SLESDestroy(_sles);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

#else

  // 2.2.0 & newer style
  ierr = KSPDestroy(_ksp);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

#endif

	     
      // Mimic PETSc default solver and preconditioner
      this->_solver_type           = GMRES;

      if (libMesh::n_processors() == 1)
	this->_preconditioner_type = ILU_PRECOND;
      else
	this->_preconditioner_type = BLOCK_JACOBI_PRECOND;
    }
}



template <typename T>
void PetscLinearSolver<T>::init ()
{
  // Initialize the data structures if not done so already.
  if (!this->initialized())
    {
      this->_is_initialized = true;
      
      int ierr=0;
  
// 2.1.x & earlier style      
#if PETSC_VERSION_LESS_THAN(2,2,0)
      
      // Create the linear solver context
      ierr = SLESCreate (libMesh::COMM_WORLD, &_sles);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
      
      // Create the Krylov subspace & preconditioner contexts
      ierr = SLESGetKSP       (_sles, &_ksp);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
      ierr = SLESGetPC        (_sles, &_pc);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
      
      // Have the Krylov subspace method use our good initial guess rather than 0
      ierr = KSPSetInitialGuessNonzero (_ksp, PETSC_TRUE);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
      
      // Set user-specified  solver and preconditioner types
      this->set_petsc_solver_type();
      this->set_petsc_preconditioner_type();
      
      // Set the options from user-input
      // Set runtime options, e.g.,
      //      -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
      //  These options will override those specified above as long as
      //  SLESSetFromOptions() is called _after_ any other customization
      //  routines.
      
      ierr = SLESSetFromOptions (_sles);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);

// 2.2.0 & newer style
#else
      
      // Create the linear solver context
      ierr = KSPCreate (libMesh::COMM_WORLD, &_ksp);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
      
      // Create the preconditioner context
      ierr = KSPGetPC        (_ksp, &_pc);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);

      // Have the Krylov subspace method use our good initial guess rather than 0
      ierr = KSPSetInitialGuessNonzero (_ksp, PETSC_TRUE);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
      
      // Set user-specified  solver and preconditioner types
      this->set_petsc_solver_type();
      this->set_petsc_preconditioner_type();
      
      // Set the options from user-input
      // Set runtime options, e.g.,
      //      -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
      //  These options will override those specified above as long as
      //  KSPSetFromOptions() is called _after_ any other customization
      //  routines.
      
      ierr = KSPSetFromOptions (_ksp);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);

      // Not sure if this is necessary, or if it is already handled by KSPSetFromOptions?
      // NOT NECESSARY!!!!
      //ierr = PCSetFromOptions (_pc);
      //CHKERRABORT(libMesh::COMM_WORLD,ierr);

	       
#endif
	     
      // Notify PETSc of location to store residual history.
      // This needs to be called before any solves, since
      // it sets the residual history length to zero.  The default
      // behavior is for PETSc to allocate (internally) an array
      // of size 1000 to hold the residual norm history.
      ierr = KSPSetResidualHistory(_ksp,
				   PETSC_NULL,   // pointer to the array which holds the history
				   PETSC_DECIDE, // size of the array holding the history
				   PETSC_TRUE);  // Whether or not to reset the history for each solve. 
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
    }
}


template <typename T>
void PetscLinearSolver<T>::init ( PetscMatrix<T>* matrix )
{
  // Initialize the data structures if not done so already.
  if (!this->initialized())
    {
      this->_is_initialized = true;
      
      int ierr=0;

// 2.1.x & earlier style      
#if PETSC_VERSION_LESS_THAN(2,2,0)
      
      // Create the linear solver context
      ierr = SLESCreate (libMesh::COMM_WORLD, &_sles);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
      
      // Create the Krylov subspace & preconditioner contexts
      ierr = SLESGetKSP       (_sles, &_ksp);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
      ierr = SLESGetPC        (_sles, &_pc);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
      
      // Have the Krylov subspace method use our good initial guess rather than 0
      ierr = KSPSetInitialGuessNonzero (_ksp, PETSC_TRUE);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
      
      // Set user-specified  solver and preconditioner types
      this->set_petsc_solver_type();
      this->set_petsc_preconditioner_type();
      
      // Set the options from user-input
      // Set runtime options, e.g.,
      //      -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
      //  These options will override those specified above as long as
      //  SLESSetFromOptions() is called _after_ any other customization
      //  routines.
      
      ierr = SLESSetFromOptions (_sles);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);

// 2.2.0 & newer style
#else
      
      // Create the linear solver context
      ierr = KSPCreate (libMesh::COMM_WORLD, &_ksp);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
      
    
      //ierr = PCCreate (libMesh::COMM_WORLD, &_pc);
      //     CHKERRABORT(libMesh::COMM_WORLD,ierr);
      
      // Create the preconditioner context
      ierr = KSPGetPC        (_ksp, &_pc);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
      
      // Set operators. The input matrix works as the preconditioning matrix
      ierr = KSPSetOperators(_ksp, matrix->mat(), matrix->mat(),SAME_NONZERO_PATTERN);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);

      // Have the Krylov subspace method use our good initial guess rather than 0
      ierr = KSPSetInitialGuessNonzero (_ksp, PETSC_TRUE);
             CHKERRABORT(libMesh::COMM_WORLD,ierr);
      
      // Set user-specified  solver and preconditioner types
      this->set_petsc_solver_type();
      this->set_petsc_preconditioner_type();
      
      // Set the options from user-input
      // Set runtime options, e.g.,
      //      -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
      //  These options will override those specified above as long as
      //  KSPSetFromOptions() is called _after_ any other customization
      //  routines.
      
      ierr = KSPSetFromOptions (_ksp);
      CHKERRABORT(libMesh::COMM_WORLD,ierr);

      // Not sure if this is necessary, or if it is already handled by KSPSetFromOptions?
      // NOT NECESSARY!!!!
      //ierr = PCSetFromOptions (_pc);
      //CHKERRABORT(libMesh::COMM_WORLD,ierr);

	       
#endif
	     
      // Notify PETSc of location to store residual history.
      // This needs to be called before any solves, since
      // it sets the residual history length to zero.  The default
      // behavior is for PETSc to allocate (internally) an array
      // of size 1000 to hold the residual norm history.
      ierr = KSPSetResidualHistory(_ksp,
				   PETSC_NULL,   // pointer to the array which holds the history
				   PETSC_DECIDE, // size of the array holding the history
				   PETSC_TRUE);  // Whether or not to reset the history for each solve. 
      CHKERRABORT(libMesh::COMM_WORLD,ierr);
    }
}






template <typename T>
std::pair<unsigned int, Real> 
PetscLinearSolver<T>::solve (SparseMatrix<T>&  matrix_in,
			     SparseMatrix<T>&  precond_in,
			     NumericVector<T>& solution_in,
			     NumericVector<T>& rhs_in,
			     const double tol,
			     const unsigned int m_its)
{
  START_LOG("solve()", "PetscLinearSolver");
  
  PetscMatrix<T>* matrix   = dynamic_cast<PetscMatrix<T>*>(&matrix_in);
  PetscMatrix<T>* precond  = dynamic_cast<PetscMatrix<T>*>(&precond_in);
  PetscVector<T>* solution = dynamic_cast<PetscVector<T>*>(&solution_in);
  PetscVector<T>* rhs      = dynamic_cast<PetscVector<T>*>(&rhs_in);

  // We cast to pointers so we can be sure that they succeeded
  // by comparing the result against NULL.
  libmesh_assert(matrix   != NULL);
  libmesh_assert(precond  != NULL);
  libmesh_assert(solution != NULL);
  libmesh_assert(rhs      != NULL);
  
  this->init (matrix);

  int ierr=0;
  int its=0, max_its = static_cast<int>(m_its);
  PetscReal final_resid=0.;

  // Close the matrices and vectors in case this wasn't already done.
  matrix->close ();
  precond->close ();
  solution->close ();
  rhs->close ();

//   // If matrix != precond, then this means we have specified a
//   // special preconditioner, so reset preconditioner type to PCMAT.
//   if (matrix != precond)
//     {
//       this->_preconditioner_type = USER_PRECOND;
//       this->set_petsc_preconditioner_type ();
//     }
  
// 2.1.x & earlier style      
#if PETSC_VERSION_LESS_THAN(2,2,0)
      
  // Set operators. The input matrix works as the preconditioning matrix
  ierr = SLESSetOperators(_sles, matrix->mat(), precond->mat(),
			  SAME_NONZERO_PATTERN);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);


  // Set the tolerances for the iterative solver.  Use the user-supplied
  // tolerance for the relative residual & leave the others at default values.
  ierr = KSPSetTolerances (_ksp, tol, PETSC_DEFAULT,
 			   PETSC_DEFAULT, max_its);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);


  // Solve the linear system
  ierr = SLESSolve (_sles, rhs->vec(), solution->vec(), &its);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);