petsc_linear_solver.c

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

  // Get the norm of the final residual to return to the user.
  ierr = KSPGetResidualNorm (_ksp, &final_resid);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);

// 2.2.0
#elif PETSC_VERSION_LESS_THAN(2,2,1)
      
  // Set operators. The input matrix works as the preconditioning matrix
  //ierr = KSPSetOperators(_ksp, 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.
  // Convergence is detected at iteration k if
  // ||r_k||_2 < max(rtol*||b||_2 , abstol)
  // where r_k is the residual vector and b is the right-hand side.  Note that
  // it is the *maximum* of the two values, the larger of which will almost
  // always be rtol*||b||_2. 
  ierr = KSPSetTolerances (_ksp,
			   tol,           // rtol   = relative decrease in residual  (1.e-5)
			   PETSC_DEFAULT, // abstol = absolute convergence tolerance (1.e-50)
 			   PETSC_DEFAULT, // dtol   = divergence tolerance           (1.e+5)
			   max_its);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);


  // Set the solution vector to use
  ierr = KSPSetSolution (_ksp, solution->vec());
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
	 
  // Set the RHS vector to use
  ierr = KSPSetRhs (_ksp, rhs->vec());
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
   
  // Solve the linear system
  ierr = KSPSolve (_ksp);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
	 
  // Get the number of iterations required for convergence
  ierr = KSPGetIterationNumber (_ksp, &its);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
	 
  // Get the norm of the final residual to return to the user.
  ierr = KSPGetResidualNorm (_ksp, &final_resid);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
	 
// 2.2.1 & newer style
#else
      
  // Set operators. The input matrix works as the preconditioning matrix
  ierr = KSPSetOperators(_ksp, 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 = KSPSolve (_ksp, rhs->vec(), solution->vec());
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
	 
  // Get the number of iterations required for convergence
  ierr = KSPGetIterationNumber (_ksp, &its);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
	 
  // Get the norm of the final residual to return to the user.
  ierr = KSPGetResidualNorm (_ksp, &final_resid);
         CHKERRABORT(libMesh::COMM_WORLD,ierr);
	 
#endif

  STOP_LOG("solve()", "PetscLinearSolver");
  // return the # of its. and the final residual norm.
  return std::make_pair(its, final_resid);
}



template <typename T>
void PetscLinearSolver<T>::get_residual_history(std::vector<double>& hist)
{
  int ierr = 0;
  int its  = 0;

  // Fill the residual history vector with the residual norms
  // Note that GetResidualHistory() does not copy any values, it
  // simply sets the pointer p.  Note that for some Krylov subspace
  // methods, the number of residuals returned in the history
  // vector may be different from what you are expecting.  For
  // example, TFQMR returns two residual values per iteration step.
  PetscReal* p;
  ierr = KSPGetResidualHistory(_ksp, &p, &its);
  CHKERRABORT(libMesh::COMM_WORLD,ierr);

  // Check for early return
  if (its == 0) return;
  
  // Create space to store the result
  hist.resize(its);

  // Copy history into the vector provided by the user.
  for (int i=0; i<its; ++i)
    {
      hist[i] = *p;
      p++;
    }
}




template <typename T>
Real PetscLinearSolver<T>::get_initial_residual()
{
  int ierr = 0;
  int its  = 0;

  // Fill the residual history vector with the residual norms
  // Note that GetResidualHistory() does not copy any values, it
  // simply sets the pointer p.  Note that for some Krylov subspace
  // methods, the number of residuals returned in the history
  // vector may be different from what you are expecting.  For
  // example, TFQMR returns two residual values per iteration step.
  PetscReal* p;
  ierr = KSPGetResidualHistory(_ksp, &p, &its);
  CHKERRABORT(libMesh::COMM_WORLD,ierr);

  // Check no residual history
  if (its == 0)
    {
      std::cerr << "No iterations have been performed, returning 0." << std::endl;
      return 0.;
    }

  // Otherwise, return the value pointed to by p.
  return *p;
}




template <typename T>
void PetscLinearSolver<T>::set_petsc_solver_type()
{
  int ierr = 0;
  
  switch (this->_solver_type)
    {

    case CG:
      ierr = KSPSetType (_ksp, (char*) KSPCG);         CHKERRABORT(libMesh::COMM_WORLD,ierr); return;

    case CR:
      ierr = KSPSetType (_ksp, (char*) KSPCR);         CHKERRABORT(libMesh::COMM_WORLD,ierr); return;

    case CGS:
      ierr = KSPSetType (_ksp, (char*) KSPCGS);        CHKERRABORT(libMesh::COMM_WORLD,ierr); return;

    case BICG:
      ierr = KSPSetType (_ksp, (char*) KSPBICG);       CHKERRABORT(libMesh::COMM_WORLD,ierr); return;

    case TCQMR:
      ierr = KSPSetType (_ksp, (char*) KSPTCQMR);      CHKERRABORT(libMesh::COMM_WORLD,ierr); return;
 
    case TFQMR:
      ierr = KSPSetType (_ksp, (char*) KSPTFQMR);      CHKERRABORT(libMesh::COMM_WORLD,ierr); return;

    case LSQR:
      ierr = KSPSetType (_ksp, (char*) KSPLSQR);       CHKERRABORT(libMesh::COMM_WORLD,ierr); return;

    case BICGSTAB:
      ierr = KSPSetType (_ksp, (char*) KSPBCGS);       CHKERRABORT(libMesh::COMM_WORLD,ierr); return;

    case MINRES:
      ierr = KSPSetType (_ksp, (char*) KSPMINRES);     CHKERRABORT(libMesh::COMM_WORLD,ierr); return;

    case GMRES:
      ierr = KSPSetType (_ksp, (char*) KSPGMRES);      CHKERRABORT(libMesh::COMM_WORLD,ierr); return;

    case RICHARDSON:
      ierr = KSPSetType (_ksp, (char*) KSPRICHARDSON); CHKERRABORT(libMesh::COMM_WORLD,ierr); return;

    case CHEBYSHEV: 
      ierr = KSPSetType (_ksp, (char*) KSPCHEBYCHEV);  CHKERRABORT(libMesh::COMM_WORLD,ierr); return;

    default:
      std::cerr << "ERROR:  Unsupported PETSC Solver: "
		<< this->_solver_type               << std::endl
		<< "Continuing with PETSC defaults" << std::endl;
    }
}








template <typename T>
void PetscLinearSolver<T>::set_petsc_preconditioner_type()
{
  int ierr = 0;
 
  switch (this->_preconditioner_type)
    {
    case IDENTITY_PRECOND:
      ierr = PCSetType (_pc, (char*) PCNONE);      CHKERRABORT(libMesh::COMM_WORLD,ierr); return;
	
    case CHOLESKY_PRECOND:
      ierr = PCSetType (_pc, (char*) PCCHOLESKY);  CHKERRABORT(libMesh::COMM_WORLD,ierr); return;

    case ICC_PRECOND:
      ierr = PCSetType (_pc, (char*) PCICC);       CHKERRABORT(libMesh::COMM_WORLD,ierr); return;

    case ILU_PRECOND:
      ierr = PCSetType (_pc, (char*) PCILU);       CHKERRABORT(libMesh::COMM_WORLD,ierr); return;

    case LU_PRECOND:
      ierr = PCSetType (_pc, (char*) PCLU);        CHKERRABORT(libMesh::COMM_WORLD,ierr); return;
      
    case ASM_PRECOND:
      ierr = PCSetType (_pc, (char*) PCASM);       CHKERRABORT(libMesh::COMM_WORLD,ierr); return;

    case JACOBI_PRECOND:
      ierr = PCSetType (_pc, (char*) PCJACOBI);    CHKERRABORT(libMesh::COMM_WORLD,ierr); return;

    case BLOCK_JACOBI_PRECOND:
      ierr = PCSetType (_pc, (char*) PCBJACOBI);   CHKERRABORT(libMesh::COMM_WORLD,ierr); return;

    case SOR_PRECOND:
      ierr = PCSetType (_pc, (char*) PCSOR);       CHKERRABORT(libMesh::COMM_WORLD,ierr); return;

    case EISENSTAT_PRECOND:
      ierr = PCSetType (_pc, (char*) PCEISENSTAT); CHKERRABORT(libMesh::COMM_WORLD,ierr); return;

#if !(PETSC_VERSION_LESS_THAN(2,1,2))
    // Only available for PETSC >= 2.1.2      
    case USER_PRECOND:
      ierr = PCSetType (_pc, (char*) PCMAT);       CHKERRABORT(libMesh::COMM_WORLD,ierr); return;
#endif

    case SHELL_PRECOND:
      ierr = PCSetType (_pc, (char*) PCSHELL);     CHKERRABORT(libMesh::COMM_WORLD,ierr); return;

    default:
      std::cerr << "ERROR:  Unsupported PETSC Preconditioner: "
		<< this->_preconditioner_type       << std::endl
		<< "Continuing with PETSC defaults" << std::endl;
    }
}

template <typename T>
void PetscLinearSolver<T>::print_converged_reason()
{
#if PETSC_VERSION_LESS_THAN(2,3,1)
  std::cout << "This method is currently not supported "
	    << "(but may work!) for Petsc 2.3.0 and earlier." << std::endl;
#else
  KSPConvergedReason reason;
  KSPGetConvergedReason(_ksp, &reason);
  
  //  KSP_CONVERGED_RTOL (residual 2-norm decreased by a factor of rtol, from 2-norm of right hand side)
  //  KSP_CONVERGED_ATOL (residual 2-norm less than abstol)
  //  KSP_CONVERGED_ITS (used by the preonly preconditioner that always uses ONE iteration) 
  //  KSP_CONVERGED_STEP_LENGTH
  //  KSP_DIVERGED_ITS  (required more than its to reach convergence)
  //  KSP_DIVERGED_DTOL (residual norm increased by a factor of divtol)
  //  KSP_DIVERGED_NAN (residual norm became Not-a-number likely do to 0/0)
  //  KSP_DIVERGED_BREAKDOWN (generic breakdown in method)

  switch (reason)
    {
    case KSP_CONVERGED_RTOL:
       {
	std::cout << "Linear solver converged, relative tolerance reached." << std::endl;
	break;
       }
    case KSP_CONVERGED_ATOL:
       {
	 std::cout << "Linear solver converged, absolute tolerance reached." << std::endl;
	 break;
       }

      // Divergence
    case KSP_DIVERGED_ITS:
       {
	 std::cout << "Linear solver diverged, max no. of iterations reached." << std::endl;
	 break;
       }
    case KSP_DIVERGED_DTOL:
       {
	 std::cout << "Linear solver diverged, residual norm increase by dtol (default 1.e5)." << std::endl;
	 break;
       }
    case KSP_DIVERGED_NAN:
       {
	 std::cout << "Linear solver diverged, residual norm is NaN." << std::endl;
	 break;
       }
    case KSP_DIVERGED_BREAKDOWN:
       {
	 std::cout << "Linear solver diverged, generic breakdown in the method." << std::endl;
	 break;
       }
    default:
      {
	std::cout << "Unknown/unsupported con(di)vergence reason: " << reason << std::endl;
      }
    }
#endif
}


//------------------------------------------------------------------
// Explicit instantiations
template class PetscLinearSolver<Number>;
 


#endif // #ifdef HAVE_PETSC