From a1464cc90497fe1ca474f874cb3977d1c48f4904 Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Mon, 17 Jan 2011 18:07:20 +0000 Subject: [PATCH] Restructured the internals of PETScWrappers::Precondition* to allow a PETSc PC object to exist without a solver. New: Precondition*::vmult(), ::initialize(), and default constructors. git-svn-id: https://svn.dealii.org/trunk@23205 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 11 +- .../include/deal.II/lac/petsc_precondition.h | 349 ++++++++++++------ deal.II/source/lac/petsc_precondition.cc | 273 ++++++++++---- deal.II/source/lac/petsc_solver.cc | 7 +- 4 files changed, 449 insertions(+), 191 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 1f0fd34a7d..982a507044 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -40,13 +40,22 @@ inconvenience this causes.

Specific improvements

    +
  1. Restructured the internals of PETScWrappers::Precondition* to allow a +PETSc PC object to exist without a solver. New: use Precondition*::vmult() to +apply the preconditioner once. Preconditioners now have a default constructor +and an initialize() function and are no longer initialized in the solver call, +but in the constructor or initialize(). +
    +(Timo Heister, 2011/01/17) +
  2. +
  3. Fixed: Boundary conditions in the step-23 tutorial program are now applied correctly. Matrix columns get eliminated with the used method and introduce some contribution to the right hand side coming from inhomogeneous boundary values. The old implementation did not reset the matrix columns before applying new boundary values.
    (Martin Stoll, Martin Kronbichler, 2011/01/14) -
+
  • Extended: base/tensor.h has an additional collection of contractions between three tensors (ie. contract3). diff --git a/deal.II/include/deal.II/lac/petsc_precondition.h b/deal.II/include/deal.II/lac/petsc_precondition.h index 6b0b26c2e2..47903a6b8e 100644 --- a/deal.II/include/deal.II/lac/petsc_precondition.h +++ b/deal.II/include/deal.II/lac/petsc_precondition.h @@ -47,31 +47,54 @@ namespace PETScWrappers * for parallel jobs, such as for example the ILU preconditioner. * * @ingroup PETScWrappers - * @author Wolfgang Bangerth, 2004 + * @author Wolfgang Bangerth, Timo Heister, 2004, 2011 */ class PreconditionerBase { public: /** - * Constructor. Take a pointer to the - * matrix from which the preconditioner - * shall be constructed. + * Constructor. */ - PreconditionerBase (const MatrixBase &matrix); + PreconditionerBase (); /** * Destructor. */ virtual ~PreconditionerBase (); + /** + * Apply the preconditioner once to the + * given src vector. + */ + void vmult (VectorBase &dst, + const VectorBase &src) const; + + /** + * Gives access to the underlying PETSc + * object. + */ + const PC & get_pc () const; + protected: + /** + * the PETSc preconditioner object + */ + PC pc; + /** * A pointer to the matrix that acts as * a preconditioner. */ - const Mat matrix; + Mat matrix; + /** + * Internal function to create the + * PETSc preconditioner object. Fails + * if called twice. + */ + void create_pc (); + /** * Conversion operator to get a * representation of the matrix that @@ -81,16 +104,7 @@ namespace PETScWrappers * the PETSc solvers. */ operator Mat () const; - - /** - * Function that takes a Krylov - * Subspace Preconditioner context - * object, and sets the type of - * preconditioner that is requested by - * the derived class. - */ - virtual void set_preconditioner_type (PC &pc) const = 0; - + /** * Make the solver class a friend, * since it needs to call the @@ -109,7 +123,7 @@ namespace PETScWrappers * preconditioner may or may not work. * * @ingroup PETScWrappers - * @author Wolfgang Bangerth, 2004 + * @author Wolfgang Bangerth, Timo Heister, 2004, 2011 */ class PreconditionJacobi : public PreconditionerBase { @@ -122,6 +136,14 @@ namespace PETScWrappers struct AdditionalData {}; + /** + * Empty Constructor. You need to call + * initialize() before using this + * object. + */ + PreconditionJacobi (); + + /** * Constructor. Take the matrix which * is used to form the preconditioner, @@ -131,21 +153,26 @@ namespace PETScWrappers PreconditionJacobi (const MatrixBase &matrix, const AdditionalData &additional_data = AdditionalData()); + /** + * Initializes the preconditioner + * object and calculate all data that + * is necessary for applying it in a + * solver. This function is + * automatically called when calling + * the constructor with the same + * arguments and is only used if you + * create the preconditioner without + * arguments. + */ + void initialize (const MatrixBase &matrix, + const AdditionalData &additional_data = AdditionalData()); + protected: /** * Store a copy of the flags for this * particular preconditioner. */ - const AdditionalData additional_data; - - /** - * Function that takes a Krylov - * Subspace Preconditioner context - * object, and sets the type of - * preconditioner that is appropriate - * for the present class. - */ - virtual void set_preconditioner_type (PC &pc) const; + AdditionalData additional_data; }; @@ -165,7 +192,7 @@ namespace PETScWrappers * preconditioner may or may not work. * * @ingroup PETScWrappers - * @author Wolfgang Bangerth, 2004 + * @author Wolfgang Bangerth, Timo Heister, 2004, 2011 */ class PreconditionBlockJacobi : public PreconditionerBase { @@ -177,6 +204,13 @@ namespace PETScWrappers */ struct AdditionalData {}; + + /** + * Empty Constructor. You need to call + * initialize() before using this + * object. + */ + PreconditionBlockJacobi (); /** * Constructor. Take the matrix which @@ -187,21 +221,26 @@ namespace PETScWrappers PreconditionBlockJacobi (const MatrixBase &matrix, const AdditionalData &additional_data = AdditionalData()); + /** + * Initializes the preconditioner + * object and calculate all data that + * is necessary for applying it in a + * solver. This function is + * automatically called when calling + * the constructor with the same + * arguments and is only used if you + * create the preconditioner without + * arguments. + */ + void initialize (const MatrixBase &matrix, + const AdditionalData &additional_data = AdditionalData()); + protected: /** * Store a copy of the flags for this * particular preconditioner. */ - const AdditionalData additional_data; - - /** - * Function that takes a Krylov - * Subspace Preconditioner context - * object, and sets the type of - * preconditioner that is appropriate - * for the present class. - */ - virtual void set_preconditioner_type (PC &pc) const; + AdditionalData additional_data; }; @@ -214,7 +253,7 @@ namespace PETScWrappers * preconditioner may or may not work. * * @ingroup PETScWrappers - * @author Wolfgang Bangerth, 2004 + * @author Wolfgang Bangerth, Timo Heister, 2004, 2011 */ class PreconditionSOR : public PreconditionerBase { @@ -239,6 +278,13 @@ namespace PETScWrappers double omega; }; + /** + * Empty Constructor. You need to call + * initialize() before using this + * object. + */ + PreconditionSOR (); + /** * Constructor. Take the matrix which * is used to form the preconditioner, @@ -248,21 +294,26 @@ namespace PETScWrappers PreconditionSOR (const MatrixBase &matrix, const AdditionalData &additional_data = AdditionalData()); + /** + * Initializes the preconditioner + * object and calculate all data that + * is necessary for applying it in a + * solver. This function is + * automatically called when calling + * the constructor with the same + * arguments and is only used if you + * create the preconditioner without + * arguments. + */ + void initialize (const MatrixBase &matrix, + const AdditionalData &additional_data = AdditionalData()); + protected: /** * Store a copy of the flags for this * particular preconditioner. */ - const AdditionalData additional_data; - - /** - * Function that takes a Krylov - * Subspace Preconditioner context - * object, and sets the type of - * preconditioner that is appropriate - * for the present class. - */ - virtual void set_preconditioner_type (PC &pc) const; + AdditionalData additional_data; }; @@ -275,7 +326,7 @@ namespace PETScWrappers * preconditioner may or may not work. * * @ingroup PETScWrappers - * @author Wolfgang Bangerth, 2004 + * @author Wolfgang Bangerth, Timo Heister, 2004, 2011 */ class PreconditionSSOR : public PreconditionerBase { @@ -300,6 +351,13 @@ namespace PETScWrappers double omega; }; + /** + * Empty Constructor. You need to call + * initialize() before using this + * object. + */ + PreconditionSSOR (); + /** * Constructor. Take the matrix which * is used to form the preconditioner, @@ -309,21 +367,26 @@ namespace PETScWrappers PreconditionSSOR (const MatrixBase &matrix, const AdditionalData &additional_data = AdditionalData()); + /** + * Initializes the preconditioner + * object and calculate all data that + * is necessary for applying it in a + * solver. This function is + * automatically called when calling + * the constructor with the same + * arguments and is only used if you + * create the preconditioner without + * arguments. + */ + void initialize (const MatrixBase &matrix, + const AdditionalData &additional_data = AdditionalData()); + protected: /** * Store a copy of the flags for this * particular preconditioner. */ - const AdditionalData additional_data; - - /** - * Function that takes a Krylov - * Subspace Preconditioner context - * object, and sets the type of - * preconditioner that is appropriate - * for the present class. - */ - virtual void set_preconditioner_type (PC &pc) const; + AdditionalData additional_data; }; @@ -336,7 +399,7 @@ namespace PETScWrappers * preconditioner may or may not work. * * @ingroup PETScWrappers - * @author Wolfgang Bangerth, 2004 + * @author Wolfgang Bangerth, Timo Heister, 2004, 2011 */ class PreconditionEisenstat : public PreconditionerBase { @@ -361,6 +424,13 @@ namespace PETScWrappers double omega; }; + /** + * Empty Constructor. You need to call + * initialize() before using this + * object. + */ + PreconditionEisenstat (); + /** * Constructor. Take the matrix which * is used to form the preconditioner, @@ -370,21 +440,26 @@ namespace PETScWrappers PreconditionEisenstat (const MatrixBase &matrix, const AdditionalData &additional_data = AdditionalData()); + /** + * Initializes the preconditioner + * object and calculate all data that + * is necessary for applying it in a + * solver. This function is + * automatically called when calling + * the constructor with the same + * arguments and is only used if you + * create the preconditioner without + * arguments. + */ + void initialize (const MatrixBase &matrix, + const AdditionalData &additional_data = AdditionalData()); + protected: /** * Store a copy of the flags for this * particular preconditioner. */ - const AdditionalData additional_data; - - /** - * Function that takes a Krylov - * Subspace Preconditioner context - * object, and sets the type of - * preconditioner that is appropriate - * for the present class. - */ - virtual void set_preconditioner_type (PC &pc) const; + AdditionalData additional_data; }; @@ -397,7 +472,7 @@ namespace PETScWrappers * preconditioner may or may not work. * * @ingroup PETScWrappers - * @author Wolfgang Bangerth, 2004 + * @author Wolfgang Bangerth, Timo Heister, 2004, 2011 */ class PreconditionICC : public PreconditionerBase { @@ -422,6 +497,13 @@ namespace PETScWrappers unsigned int levels; }; + /** + * Empty Constructor. You need to call + * initialize() before using this + * object. + */ + PreconditionICC (); + /** * Constructor. Take the matrix which * is used to form the preconditioner, @@ -431,21 +513,26 @@ namespace PETScWrappers PreconditionICC (const MatrixBase &matrix, const AdditionalData &additional_data = AdditionalData()); + /** + * Initializes the preconditioner + * object and calculate all data that + * is necessary for applying it in a + * solver. This function is + * automatically called when calling + * the constructor with the same + * arguments and is only used if you + * create the preconditioner without + * arguments. + */ + void initialize (const MatrixBase &matrix, + const AdditionalData &additional_data = AdditionalData()); + protected: /** * Store a copy of the flags for this * particular preconditioner. */ - const AdditionalData additional_data; - - /** - * Function that takes a Krylov - * Subspace Preconditioner context - * object, and sets the type of - * preconditioner that is appropriate - * for the present class. - */ - virtual void set_preconditioner_type (PC &pc) const; + AdditionalData additional_data; }; @@ -458,7 +545,7 @@ namespace PETScWrappers * preconditioner may or may not work. * * @ingroup PETScWrappers - * @author Wolfgang Bangerth, 2004 + * @author Wolfgang Bangerth, Timo Heister, 2004, 2011 */ class PreconditionILU : public PreconditionerBase { @@ -483,6 +570,13 @@ namespace PETScWrappers unsigned int levels; }; + /** + * Empty Constructor. You need to call + * initialize() before using this + * object. + */ + PreconditionILU (); + /** * Constructor. Take the matrix which * is used to form the preconditioner, @@ -492,21 +586,26 @@ namespace PETScWrappers PreconditionILU (const MatrixBase &matrix, const AdditionalData &additional_data = AdditionalData()); + /** + * Initializes the preconditioner + * object and calculate all data that + * is necessary for applying it in a + * solver. This function is + * automatically called when calling + * the constructor with the same + * arguments and is only used if you + * create the preconditioner without + * arguments. + */ + void initialize (const MatrixBase &matrix, + const AdditionalData &additional_data = AdditionalData()); + protected: /** * Store a copy of the flags for this * particular preconditioner. */ - const AdditionalData additional_data; - - /** - * Function that takes a Krylov - * Subspace Preconditioner context - * object, and sets the type of - * preconditioner that is appropriate - * for the present class. - */ - virtual void set_preconditioner_type (PC &pc) const; + AdditionalData additional_data; }; @@ -568,6 +667,13 @@ namespace PETScWrappers double damping; }; + /** + * Empty Constructor. You need to call + * initialize() before using this + * object. + */ + PreconditionLU (); + /** * Constructor. Take the matrix which * is used to form the preconditioner, @@ -577,21 +683,26 @@ namespace PETScWrappers PreconditionLU (const MatrixBase &matrix, const AdditionalData &additional_data = AdditionalData()); + /** + * Initializes the preconditioner + * object and calculate all data that + * is necessary for applying it in a + * solver. This function is + * automatically called when calling + * the constructor with the same + * arguments and is only used if you + * create the preconditioner without + * arguments. + */ + void initialize (const MatrixBase &matrix, + const AdditionalData &additional_data = AdditionalData()); + protected: /** * Store a copy of the flags for this * particular preconditioner. */ - const AdditionalData additional_data; - - /** - * Function that takes a Krylov - * Subspace Preconditioner context - * object, and sets the type of - * preconditioner that is appropriate - * for the present class. - */ - virtual void set_preconditioner_type (PC &pc) const; + AdditionalData additional_data; }; @@ -696,6 +807,13 @@ namespace PETScWrappers + /** + * Empty Constructor. You need to call + * initialize() before using this + * object. + */ + PreconditionBoomerAMG (); + /** * Constructor. Take the matrix which * is used to form the preconditioner, @@ -705,21 +823,26 @@ namespace PETScWrappers PreconditionBoomerAMG (const MatrixBase &matrix, const AdditionalData &additional_data = AdditionalData()); + /** + * Initializes the preconditioner + * object and calculate all data that + * is necessary for applying it in a + * solver. This function is + * automatically called when calling + * the constructor with the same + * arguments and is only used if you + * create the preconditioner without + * arguments. + */ + void initialize (const MatrixBase &matrix, + const AdditionalData &additional_data = AdditionalData()); + protected: /** * Store a copy of the flags for this * particular preconditioner. */ - const AdditionalData additional_data; - - /** - * Function that takes a Krylov - * Subspace Preconditioner context - * object, and sets the type of - * preconditioner that is appropriate - * for the present class. - */ - virtual void set_preconditioner_type (PC &pc) const; + AdditionalData additional_data; }; } diff --git a/deal.II/source/lac/petsc_precondition.cc b/deal.II/source/lac/petsc_precondition.cc index d1b0e496ea..33d54111e1 100644 --- a/deal.II/source/lac/petsc_precondition.cc +++ b/deal.II/source/lac/petsc_precondition.cc @@ -19,6 +19,7 @@ # include # include # include +# include # include # include @@ -26,18 +27,64 @@ DEAL_II_NAMESPACE_OPEN namespace PETScWrappers { - PreconditionerBase::PreconditionerBase (const MatrixBase &matrix) - : - matrix (matrix) + PreconditionerBase::PreconditionerBase () + : + pc(NULL), matrix(NULL) {} - PreconditionerBase::~PreconditionerBase () - {} + { + if (pc!=NULL) + { + int ierr = PCDestroy(pc); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + } + } + + + void + PreconditionerBase::vmult (VectorBase &dst, + const VectorBase &src) const + { + AssertThrow (pc != NULL, StandardExceptions::ExcInvalidState ()); + + int ierr; + ierr = PCApply(pc, src, dst); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + } + + void + PreconditionerBase::create_pc () + { + // only allow the creation of the + // preconditioner once + AssertThrow (pc == NULL, StandardExceptions::ExcInvalidState ()); + + MPI_Comm comm; + int ierr; + // this ugly cast is necessay because the + // type Mat and PETScObject are + // unrelated. + ierr = PetscObjectGetComm(reinterpret_cast(matrix), &comm); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + + ierr = PCCreate(comm, &pc); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + ierr = PCSetOperators(pc , matrix, matrix, SAME_PRECONDITIONER); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + } + + const PC & + PreconditionerBase::get_pc () const + { + return pc; + } + + PreconditionerBase::operator Mat () const { return matrix; @@ -46,52 +93,69 @@ namespace PETScWrappers /* ----------------- PreconditionJacobi -------------------- */ - + PreconditionJacobi::PreconditionJacobi () + {} + + PreconditionJacobi::PreconditionJacobi (const MatrixBase &matrix, const AdditionalData &additional_data) - : - PreconditionerBase (matrix), - additional_data (additional_data) - {} + { + initialize(matrix, additional_data); + } void - PreconditionJacobi::set_preconditioner_type (PC &pc) const + PreconditionJacobi::initialize (const MatrixBase &matrix_, + const AdditionalData &additional_data_) { - // set the right type for the - // preconditioner + matrix = static_cast(matrix_); + additional_data = additional_data_; + + create_pc(); + int ierr; ierr = PCSetType (pc, const_cast(PCJACOBI)); AssertThrow (ierr == 0, ExcPETScError(ierr)); ierr = PCSetFromOptions (pc); AssertThrow (ierr == 0, ExcPETScError(ierr)); + + ierr = PCSetUp (pc); + AssertThrow (ierr == 0, ExcPETScError(ierr)); } -/* ----------------- PreconditionJacobi -------------------- */ - +/* ----------------- PreconditionBlockJacobi -------------------- */ + PreconditionBlockJacobi::PreconditionBlockJacobi () + {} + + PreconditionBlockJacobi:: PreconditionBlockJacobi (const MatrixBase &matrix, const AdditionalData &additional_data) - : - PreconditionerBase (matrix), - additional_data (additional_data) - {} - + { + initialize(matrix, additional_data); + } void - PreconditionBlockJacobi::set_preconditioner_type (PC &pc) const + PreconditionBlockJacobi::initialize (const MatrixBase &matrix_, + const AdditionalData &additional_data_) { - // set the right type for the - // preconditioner + matrix = static_cast(matrix_); + additional_data = additional_data_; + + create_pc(); + int ierr; ierr = PCSetType (pc, const_cast(PCBJACOBI)); AssertThrow (ierr == 0, ExcPETScError(ierr)); ierr = PCSetFromOptions (pc); AssertThrow (ierr == 0, ExcPETScError(ierr)); + + ierr = PCSetUp (pc); + AssertThrow (ierr == 0, ExcPETScError(ierr)); } @@ -104,20 +168,26 @@ namespace PETScWrappers {} + PreconditionSOR::PreconditionSOR () + {} + PreconditionSOR::PreconditionSOR (const MatrixBase &matrix, const AdditionalData &additional_data) - : - PreconditionerBase (matrix), - additional_data (additional_data) - {} + { + initialize(matrix, additional_data); + } void - PreconditionSOR::set_preconditioner_type (PC &pc) const + PreconditionSOR::initialize (const MatrixBase &matrix_, + const AdditionalData &additional_data_) { - // set the right type for the - // preconditioner + matrix = static_cast(matrix_); + additional_data = additional_data_; + + create_pc(); + int ierr; ierr = PCSetType (pc, const_cast(PCSOR)); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -128,6 +198,9 @@ namespace PETScWrappers ierr = PCSetFromOptions (pc); AssertThrow (ierr == 0, ExcPETScError(ierr)); + + ierr = PCSetUp (pc); + AssertThrow (ierr == 0, ExcPETScError(ierr)); } @@ -139,21 +212,27 @@ namespace PETScWrappers omega (omega) {} + + PreconditionSSOR::PreconditionSSOR () + {} PreconditionSSOR::PreconditionSSOR (const MatrixBase &matrix, const AdditionalData &additional_data) - : - PreconditionerBase (matrix), - additional_data (additional_data) - {} + { + initialize(matrix, additional_data); + } void - PreconditionSSOR::set_preconditioner_type (PC &pc) const + PreconditionSSOR::initialize (const MatrixBase &matrix_, + const AdditionalData &additional_data_) { - // set the right type for the - // preconditioner + matrix = static_cast(matrix_); + additional_data = additional_data_; + + create_pc(); + int ierr; ierr = PCSetType (pc, const_cast(PCSOR)); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -168,6 +247,9 @@ namespace PETScWrappers ierr = PCSetFromOptions (pc); AssertThrow (ierr == 0, ExcPETScError(ierr)); + + ierr = PCSetUp (pc); + AssertThrow (ierr == 0, ExcPETScError(ierr)); } @@ -180,20 +262,26 @@ namespace PETScWrappers {} + PreconditionEisenstat::PreconditionEisenstat () + {} + PreconditionEisenstat::PreconditionEisenstat (const MatrixBase &matrix, const AdditionalData &additional_data) - : - PreconditionerBase (matrix), - additional_data (additional_data) - {} + { + initialize(matrix, additional_data); + } void - PreconditionEisenstat::set_preconditioner_type (PC &pc) const + PreconditionEisenstat::initialize (const MatrixBase &matrix_, + const AdditionalData &additional_data_) { - // set the right type for the - // preconditioner + matrix = static_cast(matrix_); + additional_data = additional_data_; + + create_pc(); + int ierr; ierr = PCSetType (pc, const_cast(PCEISENSTAT)); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -204,6 +292,9 @@ namespace PETScWrappers ierr = PCSetFromOptions (pc); AssertThrow (ierr == 0, ExcPETScError(ierr)); + + ierr = PCSetUp (pc); + AssertThrow (ierr == 0, ExcPETScError(ierr)); } @@ -217,20 +308,26 @@ namespace PETScWrappers {} + PreconditionICC::PreconditionICC () + {} + PreconditionICC::PreconditionICC (const MatrixBase &matrix, const AdditionalData &additional_data) - : - PreconditionerBase (matrix), - additional_data (additional_data) - {} + { + initialize(matrix, additional_data); + } void - PreconditionICC::set_preconditioner_type (PC &pc) const + PreconditionICC::initialize (const MatrixBase &matrix_, + const AdditionalData &additional_data_) { - // set the right type for the - // preconditioner + matrix = static_cast(matrix_); + additional_data = additional_data_; + + create_pc(); + int ierr; ierr = PCSetType (pc, const_cast(PCICC)); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -245,6 +342,9 @@ namespace PETScWrappers ierr = PCSetFromOptions (pc); AssertThrow (ierr == 0, ExcPETScError(ierr)); + + ierr = PCSetUp (pc); + AssertThrow (ierr == 0, ExcPETScError(ierr)); } @@ -257,20 +357,26 @@ namespace PETScWrappers {} + PreconditionILU::PreconditionILU () + {} + PreconditionILU::PreconditionILU (const MatrixBase &matrix, const AdditionalData &additional_data) - : - PreconditionerBase (matrix), - additional_data (additional_data) - {} + { + initialize(matrix, additional_data); + } void - PreconditionILU::set_preconditioner_type (PC &pc) const + PreconditionILU::initialize (const MatrixBase &matrix_, + const AdditionalData &additional_data_) { - // set the right type for the - // preconditioner + matrix = static_cast(matrix_); + additional_data = additional_data_; + + create_pc(); + int ierr; ierr = PCSetType (pc, const_cast(PCILU)); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -285,6 +391,9 @@ namespace PETScWrappers ierr = PCSetFromOptions (pc); AssertThrow (ierr == 0, ExcPETScError(ierr)); + + ierr = PCSetUp (pc); + AssertThrow (ierr == 0, ExcPETScError(ierr)); } @@ -306,20 +415,27 @@ namespace PETScWrappers {} + PreconditionBoomerAMG::PreconditionBoomerAMG () + {} + + PreconditionBoomerAMG::PreconditionBoomerAMG (const MatrixBase &matrix, const AdditionalData &additional_data) - : - PreconditionerBase (matrix), - additional_data (additional_data) - {} + { + initialize(matrix, additional_data); + } void - PreconditionBoomerAMG::set_preconditioner_type (PC &pc) const + PreconditionBoomerAMG::initialize (const MatrixBase &matrix_, + const AdditionalData &additional_data_) { + matrix = static_cast(matrix_); + additional_data = additional_data_; + #ifdef PETSC_HAVE_HYPRE - // set the right type for the - // preconditioner + create_pc(); + int ierr; ierr = PCSetType (pc, const_cast(PCHYPRE)); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -358,6 +474,10 @@ namespace PETScWrappers ierr = PCSetFromOptions (pc); AssertThrow (ierr == 0, ExcPETScError(ierr)); + + ierr = PCSetUp (pc); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + #else // PETSC_HAVE_HYPRE (void)pc; Assert (false, @@ -380,20 +500,26 @@ namespace PETScWrappers {} + PreconditionLU::PreconditionLU () + {} + PreconditionLU::PreconditionLU (const MatrixBase &matrix, const AdditionalData &additional_data) - : - PreconditionerBase (matrix), - additional_data (additional_data) - {} + { + initialize(matrix, additional_data); + } void - PreconditionLU::set_preconditioner_type (PC &pc) const + PreconditionLU::initialize (const MatrixBase &matrix_, + const AdditionalData &additional_data_) { - // set the right type for the - // preconditioner + matrix = static_cast(matrix_); + additional_data = additional_data_; + + create_pc(); + int ierr; ierr = PCSetType (pc, const_cast(PCLU)); AssertThrow (ierr == 0, ExcPETScError(ierr)); @@ -427,6 +553,9 @@ namespace PETScWrappers ierr = PCSetFromOptions (pc); AssertThrow (ierr == 0, ExcPETScError(ierr)); + + ierr = PCSetUp (pc); + AssertThrow (ierr == 0, ExcPETScError(ierr)); } diff --git a/deal.II/source/lac/petsc_solver.cc b/deal.II/source/lac/petsc_solver.cc index 8f5af44aa6..2c8fbdaef6 100644 --- a/deal.II/source/lac/petsc_solver.cc +++ b/deal.II/source/lac/petsc_solver.cc @@ -68,7 +68,6 @@ namespace PETScWrappers const PreconditionerBase &preconditioner) { int ierr; - // first create a solver object if this // is necessary if (solver_data.get() == 0) @@ -92,10 +91,8 @@ namespace PETScWrappers // preconditioner set_solver_type (solver_data->ksp); - ierr = KSPGetPC (solver_data->ksp, &solver_data->pc); - AssertThrow (ierr == 0, ExcPETScError(ierr)); - - preconditioner.set_preconditioner_type (solver_data->pc); + ierr = KSPSetPC (solver_data->ksp, preconditioner.get_pc()); + AssertThrow (ierr == 0, ExcPETScError(ierr)); // then a convergence monitor // function. that function simply -- 2.39.5