Index: include/numerics/petsc_preconditioner.h =================================================================== --- include/numerics/petsc_preconditioner.h (revision 0) +++ include/numerics/petsc_preconditioner.h (revision 0) @@ -0,0 +1,133 @@ +// $Id: petsc_preconditioner.h 3158 2009-02-02 17:44:57Z drgasto $ + +// The libMesh Finite Element Library. +// Copyright (C) 2002-2008 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner + +// 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 + + + +#ifndef __petsc_preconditioner_h__ +#define __petsc_preconditioner_h__ + +#include "libmesh_config.h" + +#ifdef LIBMESH_HAVE_PETSC + +// C++ includes + +// Local includes +#include "preconditioner.h" +#include "libmesh_common.h" +#include "enum_solver_package.h" +#include "enum_preconditioner_type.h" +#include "reference_counted_object.h" +#include "libmesh.h" + +// Petsc includes +#include "petscpc.h" + +// forward declarations +template class AutoPtr; +template class SparseMatrix; +template class NumericVector; +template class ShellMatrix; + + +/** + * This class provides an interface to the suite of preconditioners available + * from Petsc. + * + * @author Derek Gaston, 2009 + */ + +template +class PetscPreconditioner : public Preconditioner +{ +public: + + /** + * Constructor. Initializes PetscPreconditioner data structures + */ + PetscPreconditioner (); + + /** + * Destructor. + */ + virtual ~PetscPreconditioner (); + + /** + * Computes the preconditioned vector "y" based on input "x". + * Usually by solving Py=x to get the action of P^-1 x. + */ + virtual void apply(const NumericVector & x, NumericVector & y); + + /** + * Release all memory and clear data structures. + */ + virtual void clear () {} + + /** + * Initialize data structures if not done so already. + */ + virtual void init (); + + /** + * Returns the actual Petsc PC struct. Useful for more advanced + * purposes + */ + PC pc() { return _pc; } + + /** + * Tells PETSC to use the user-specified preconditioner + */ + static void set_petsc_preconditioner_type (const PreconditionerType & preconditioner_type, PC & pc); + +protected: + + /** + * Preconditioner context + */ + PC _pc; + + /** + * Petsc Matrix that's been pulled out of the _matrix object. + * This happens during init... + */ + Mat _mat; +}; + + + + +/*----------------------- inline functions ----------------------------------*/ +template +inline +PetscPreconditioner::PetscPreconditioner () : + Preconditioner() +{ +} + + + +template +inline +PetscPreconditioner::~PetscPreconditioner () +{ + this->clear (); +} + +#endif // #ifdef LIBMESH_HAVE_PETSC +#endif // #ifdef __petsc_linear_solver_h__ Index: include/numerics/petsc_linear_solver.h =================================================================== --- include/numerics/petsc_linear_solver.h (revision 3235) +++ include/numerics/petsc_linear_solver.h (working copy) @@ -187,12 +187,6 @@ void set_petsc_solver_type (); /** - * Tells PETSC to use the user-specified preconditioner stored in - * \p _preconditioner_type - */ - void set_petsc_preconditioner_type (); - - /** * Internal function if shell matrix mode is used. */ static PetscErrorCode _petsc_shell_matrix_mult(Mat mat, Vec arg, Vec dest); Index: include/numerics/preconditioner.h =================================================================== --- include/numerics/preconditioner.h (revision 0) +++ include/numerics/preconditioner.h (revision 0) @@ -0,0 +1,178 @@ +// $Id: preconditioner.h 3158 2009-02-02 17:44:57Z drgasto $ + +// The libMesh Finite Element Library. +// Copyright (C) 2002-2008 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner + +// 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 + + + +#ifndef __preconditioner_h__ +#define __preconditioner_h__ + + +// C++ includes + +// Local includes +#include "libmesh_common.h" +#include "enum_solver_package.h" +#include "enum_solver_type.h" +#include "enum_preconditioner_type.h" +#include "reference_counted_object.h" +#include "libmesh.h" + +// forward declarations +template class AutoPtr; +template class SparseMatrix; +template class NumericVector; +template class ShellMatrix; + + + + + +/** + * This class provides a uniform interface for preconditioners. This base + * class is overloaded to provide linear solvers from different packages + * like PETSC or Trilinos. + * + * In the below comments P is the matrix to be preconditioned with Apply() + * performing the equivalent of the matrix vector product P^-1 x. This + * can also be thought of as (usually approximately) solving for Py=x. + * + * @author Derek Gaston, 2009 + */ + +template +class Preconditioner : public ReferenceCountedObject > +{ +public: + + /** + * Constructor. Initializes Preconditioner data structures + */ + Preconditioner (); + + /** + * Destructor. + */ + virtual ~Preconditioner (); + + /** + * Builds a \p Preconditioner using the linear solver package specified by + * \p solver_package + */ + static AutoPtr > build(const SolverPackage solver_package = + libMesh::default_solver_package()); + + /** + * @returns true if the data structures are + * initialized, false otherwise. + */ + bool initialized () const { return _is_initialized; } + + /** + * Computes the preconditioned vector "y" based on input "x". + * Usually by solving Py=x to get the action of P^-1 x. + */ + virtual void apply(const NumericVector & x, NumericVector & y) = 0; + + /** + * Release all memory and clear data structures. + */ + virtual void clear () {} + + /** + * Initialize data structures if not done so already. + */ + virtual void init () = 0; + + /** + * Sets the matrix P to be preconditioned. + */ + void set_matrix(SparseMatrix & mat); + + /** + * Returns the type of preconditioner to use. + */ + PreconditionerType preconditioner_type () const + { return _preconditioner_type; } + + /** + * Sets the type of preconditioner to use. + */ + void set_preconditioner_type (const PreconditionerType pct); + +protected: + + /** + * The matrix P... ie the matrix to be preconditioned. + * This is often the actual system matrix of a linear sytem. + */ + SparseMatrix * _matrix; + + /** + * Enum statitng with type of preconditioner to use. + */ + PreconditionerType _preconditioner_type; + + /** + * Flag indicating if the data structures have been initialized. + */ + bool _is_initialized; +}; + + + + +/*----------------------- inline functions ----------------------------------*/ +template +inline +Preconditioner::Preconditioner () : + _matrix(NULL), + _preconditioner_type (ILU_PRECOND), + _is_initialized (false) +{ +} + + + +template +inline +Preconditioner::~Preconditioner () +{ + this->clear (); +} + +template +void +Preconditioner::set_matrix(SparseMatrix & mat) +{ + //If the matrix is changing then we (probably) need to reinitialize. + _is_initialized = false; + _matrix = &mat; +} + +template +void +Preconditioner::set_preconditioner_type (const PreconditionerType pct) +{ + //If the preconditioner type changes we (probably) need to reinitialize. + _is_initialized = false; + _preconditioner_type = pct; +} + + +#endif // #ifdef __preconditioner_h__ Index: include/enums/enum_preconditioner_type.h =================================================================== --- include/enums/enum_preconditioner_type.h (revision 3235) +++ include/enums/enum_preconditioner_type.h (working copy) @@ -47,6 +47,7 @@ LU_PRECOND, USER_PRECOND, SHELL_PRECOND, + AMG_PRECOND, INVALID_PRECONDITIONER}; } Index: include/solvers/system.h =================================================================== --- include/solvers/system.h (revision 3235) +++ include/solvers/system.h (working copy) @@ -653,6 +653,12 @@ */ AutoPtr > current_local_solution; + /** + * Fills the std::set with the degrees of freedom on the local + * processor corresponding the the variable number passed in. + */ + void local_dof_indices (const unsigned int var, + std::set & var_indices) const; protected: @@ -676,19 +682,10 @@ */ void project_vector (const NumericVector&, NumericVector&) const; - private: /** - * Fills the std::set with the degrees of freedom on the local - * processor corresponding the the variable number passed in. - */ - void local_dof_indices (const unsigned int var, - std::set & var_indices) const; - - - /** * Finds the discrete norm for the entries in the vector * corresponding to Dofs associated with var. */ Index: src/numerics/petsc_linear_solver.C =================================================================== --- src/numerics/petsc_linear_solver.C (revision 3235) +++ src/numerics/petsc_linear_solver.C (working copy) @@ -30,9 +30,9 @@ #include "libmesh_logging.h" #include "petsc_linear_solver.h" #include "shell_matrix.h" +#include "petsc_preconditioner.h" - /*----------------------- functions ----------------------------------*/ template void PetscLinearSolver::clear () @@ -99,7 +99,7 @@ // Set user-specified solver and preconditioner types this->set_petsc_solver_type(); - this->set_petsc_preconditioner_type(); + PetscPreconditioner::set_petsc_preconditioner_type(this->_preconditioner_type,_pc); // Set the options from user-input // Set runtime options, e.g., @@ -128,7 +128,7 @@ // Set user-specified solver and preconditioner types this->set_petsc_solver_type(); - this->set_petsc_preconditioner_type(); + PetscPreconditioner::set_petsc_preconditioner_type(this->_preconditioner_type,_pc); // Set the options from user-input // Set runtime options, e.g., @@ -191,7 +191,7 @@ // Set user-specified solver and preconditioner types this->set_petsc_solver_type(); - this->set_petsc_preconditioner_type(); + PetscPreconditioner::set_petsc_preconditioner_type(this->_preconditioner_type,_pc); // Set the options from user-input // Set runtime options, e.g., @@ -228,7 +228,7 @@ // Set user-specified solver and preconditioner types this->set_petsc_solver_type(); - this->set_petsc_preconditioner_type(); + PetscPreconditioner::set_petsc_preconditioner_type(this->_preconditioner_type,_pc); // Set the options from user-input // Set runtime options, e.g., @@ -709,67 +709,7 @@ } } - - - - - - - template -void PetscLinearSolver::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 void PetscLinearSolver::print_converged_reason() { #if PETSC_VERSION_LESS_THAN(2,3,1) Index: src/numerics/preconditioner.C =================================================================== --- src/numerics/preconditioner.C (revision 0) +++ src/numerics/preconditioner.C (revision 0) @@ -0,0 +1,84 @@ +// $Id: preconditioner.C 3158 2009-02-02 17:44:57Z drgasto $ + +// The libMesh Finite Element Library. +// Copyright (C) 2002-2008 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner + +// 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 + + + +// C++ includes + +// Local Includes +#include "auto_ptr.h" +#include "preconditioner.h" +#include "petsc_preconditioner.h" + +//------------------------------------------------------------------ +// Preconditioner members +template +AutoPtr > +Preconditioner::build(const SolverPackage solver_package) +{ + // Build the appropriate solver + switch (solver_package) + { + +/* +#ifdef LIBMESH_HAVE_LASPACK + case LASPACK_SOLVERS: + { + AutoPtr > ap(new LaspackPreconditioner); + return ap; + } +#endif +*/ + +#ifdef LIBMESH_HAVE_PETSC + case PETSC_SOLVERS: + { + AutoPtr > ap(new PetscPreconditioner); + return ap; + } +#endif + +/* +#ifdef LIBMESH_HAVE_TRILINOS + case TRILINOS_SOLVERS: + { + AutoPtr > ap(new AztecPreconditioner); + return ap; + } +#endif +*/ + default: + std::cerr << "ERROR: Unrecognized solver package: " + << solver_package + << std::endl; + libmesh_error(); + } + + AutoPtr > ap(NULL); + return ap; +} + + + +//------------------------------------------------------------------ +// Explicit instantiations +template class Preconditioner; + + + Index: src/numerics/petsc_preconditioner.C =================================================================== --- src/numerics/petsc_preconditioner.C (revision 0) +++ src/numerics/petsc_preconditioner.C (revision 0) @@ -0,0 +1,145 @@ +// $Id: linear_solver.C 3158 2009-02-02 17:44:57Z drgasto $ + +// The libMesh Finite Element Library. +// Copyright (C) 2002-2008 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner + +// 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 LIBMESH_HAVE_PETSC + +// C++ includes + +// Local Includes +#include "auto_ptr.h" +#include "petsc_preconditioner.h" +#include "petsc_macro.h" +#include "petsc_matrix.h" +#include "petsc_vector.h" +#include "petsc_macro.h" + +#include "libmesh_common.h" + +template +void +PetscPreconditioner::apply(const NumericVector & x, NumericVector & y) +{ + PetscVector & x_pvec = libmesh_cast_ref&, NumericVector >(const_cast&>(x)); + PetscVector & y_pvec = libmesh_cast_ref&, NumericVector >(const_cast&>(y)); + + Vec x_vec = x_pvec.vec(); + Vec y_vec = y_pvec.vec(); + + PCApply(_pc,x_vec,y_vec); +} + +template +void +PetscPreconditioner::init () +{ + if(this->_preconditioner_type == IDENTITY_PRECOND) + return; + + if(!this->_matrix) + { + std::cerr << "ERROR: No matrix set for PetscPreconditioner, but init() called" << std::endl; + libmesh_error(); + } + + //Clear the preconditioner in case it has been created in the past + if(this->_is_initialized) + PCDestroy(_pc); + + //Create the preconditioning object + PCCreate(libMesh::COMM_WORLD,&_pc); + + //Set the PCType + set_petsc_preconditioner_type(this->_preconditioner_type, _pc); + + if(this->_preconditioner_type == AMG_PRECOND) + PCHYPRESetType(this->_pc, "boomeramg"); + + PetscMatrix * pmatrix = libmesh_cast_ptr*, SparseMatrix >(this->_matrix); + + _mat = pmatrix->mat(); + + PCSetOperators(_pc,_mat,_mat,SAME_NONZERO_PATTERN); + + this->_is_initialized = true; +} + +template +void +PetscPreconditioner::set_petsc_preconditioner_type (const PreconditionerType & preconditioner_type, PC & pc) +{ + int ierr = 0; + + switch (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; + + case AMG_PRECOND: + ierr = PCSetType (pc, (char*) PCHYPRE); 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: " + << preconditioner_type << std::endl + << "Continuing with PETSC defaults" << std::endl; + } +} + +//------------------------------------------------------------------ +// Explicit instantiations +template class PetscPreconditioner; + +#endif // #ifdef LIBMESH_HAVE_PETSC