From fb0c546672a9f1bc1f7d793841b79e3a2ff2793c Mon Sep 17 00:00:00 2001 From: Stefano Zampini Date: Wed, 23 Nov 2022 20:14:55 +0300 Subject: [PATCH] PETScWrappers: improve preconditioner class Added support for SHELL preconditioning --- include/deal.II/lac/petsc_precondition.h | 101 ++++++--- source/lac/petsc_precondition.cc | 271 +++++++++++++---------- source/lac/petsc_solver.cc | 11 +- tests/petsc/solver_03_mf.cc | 7 + tests/petsc/solver_03_mf.output | 1 + 5 files changed, 238 insertions(+), 153 deletions(-) diff --git a/include/deal.II/lac/petsc_precondition.h b/include/deal.II/lac/petsc_precondition.h index 53b708ba39..321c250d48 100644 --- a/include/deal.II/lac/petsc_precondition.h +++ b/include/deal.II/lac/petsc_precondition.h @@ -28,6 +28,8 @@ # include +# include + DEAL_II_NAMESPACE_OPEN @@ -38,7 +40,6 @@ namespace PETScWrappers # ifndef DOXYGEN class MatrixBase; class VectorBase; - class SolverBase; # endif /** @@ -65,11 +66,9 @@ namespace PETScWrappers explicit PreconditionBase(const MPI_Comm &mpi_communicator); /** - * Constructor. This constructor is deprecated. + * Constructor. * - * @deprecated */ - DEAL_II_DEPRECATED PreconditionBase(); /** @@ -79,7 +78,7 @@ namespace PETScWrappers /** * Destroys the preconditioner, leaving an object like just after having - * called the constructor. + * called the default constructor. */ void clear(); @@ -96,6 +95,12 @@ namespace PETScWrappers void Tvmult(VectorBase &dst, const VectorBase &src) const; + /** + * Explictly call setup. This is usually not needed since PETSc will + * automatically call the setup function when needed. + */ + void + setup(); /** * Give access to the underlying PETSc object. @@ -106,42 +111,27 @@ namespace PETScWrappers /** * Return the MPI communicator object used by this preconditioner. */ - MPI_Comm + const MPI_Comm & get_mpi_communicator() const; protected: - /** - * The communicator to be used for this preconditioner. - */ - MPI_Comm mpi_communicator; - /** * The PETSc preconditioner object */ PC pc; - /** - * A pointer to the matrix that acts as a preconditioner. - */ - Mat matrix; - /** * Internal function to create the PETSc preconditioner object. Fails if * called twice. */ void - create_pc(); + create_pc_with_mat(const MatrixBase &); /** - * Conversion operator to get a representation of the matrix that - * represents this preconditioner. We use this inside the actual solver, - * where we need to pass this matrix to the PETSc solvers. + * Internal function to create the PETSc preconditioner object. */ - operator Mat() const; - - // Make the solver class a friend, since it needs to call the conversion - // operator. - friend class SolverBase; + void + create_pc_with_comm(const MPI_Comm &); }; @@ -1090,6 +1080,67 @@ namespace PETScWrappers initialize(); }; + class PreconditionShell : public PreconditionBase + { + public: + /** + * Empty Constructor. You need to call initialize() before using this + * object. + */ + PreconditionShell() = default; + + /** + * Constructor. Take the matrix which is used to form the preconditioner. + */ + PreconditionShell(const MatrixBase &matrix); + + /** + * Same as above but without setting a matrix to form the preconditioner. + */ + PreconditionShell(const MPI_Comm &communicator); + + + /** + * The callback for the application of the preconditioner. Defaults to + * a copy operation if not provided. + */ + std::function apply; + + /** + * The callback for the application of the transposed preconditioner. + * Defaults to the non-transpose operation if not provided. + */ + std::function applyT; + + protected: + /** + * Initialize the preconditioner object without knowing a particular + * matrix. This function sets up the PCSHELL preconditioner + */ + void + initialize(const MPI_Comm &comm); + + /** + * Initialize the preconditioner object with a particular + * matrix. This function sets up the PCSHELL preconditioner + */ + void + initialize(const MatrixBase &matrix); + + private: + /** + * Callback-function invoked by PCApply + */ + static int + pcapply(PC pc, Vec src, Vec dst); + + /** + * Callback-function invoked by PCApplyTranspose + */ + static int + pcapply_transpose(PC pc, Vec src, Vec dst); + }; + /** * Alias for backwards-compatibility. * @deprecated Use PETScWrappers::PreconditionBase instead. diff --git a/source/lac/petsc_precondition.cc b/source/lac/petsc_precondition.cc index a39d482ca4..67fe466393 100644 --- a/source/lac/petsc_precondition.cc +++ b/source/lac/petsc_precondition.cc @@ -34,19 +34,15 @@ DEAL_II_NAMESPACE_OPEN namespace PETScWrappers { PreconditionBase::PreconditionBase(const MPI_Comm &comm) - : mpi_communicator(comm) - , pc(nullptr) - , matrix(nullptr) - {} - - + : pc(nullptr) + { + create_pc_with_comm(comm); + } PreconditionBase::PreconditionBase() - : PreconditionBase(MPI_COMM_NULL) + : pc(nullptr) {} - - PreconditionBase::~PreconditionBase() { try @@ -57,23 +53,16 @@ namespace PETScWrappers {} } - void PreconditionBase::clear() { - mpi_communicator = MPI_COMM_NULL; - - matrix = nullptr; - - if (pc != nullptr) + if (pc) { PetscErrorCode ierr = PCDestroy(&pc); - pc = nullptr; AssertThrow(ierr == 0, ExcPETScError(ierr)); } } - void PreconditionBase::vmult(VectorBase &dst, const VectorBase &src) const { @@ -83,8 +72,6 @@ namespace PETScWrappers AssertThrow(ierr == 0, ExcPETScError(ierr)); } - - void PreconditionBase::Tvmult(VectorBase &dst, const VectorBase &src) const { @@ -94,28 +81,39 @@ namespace PETScWrappers AssertThrow(ierr == 0, ExcPETScError(ierr)); } + void + PreconditionBase::setup() + { + AssertThrow(pc != nullptr, StandardExceptions::ExcInvalidState()); + PetscErrorCode ierr = PCSetUp(pc); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + } - MPI_Comm + const MPI_Comm & PreconditionBase::get_mpi_communicator() const { - return mpi_communicator; + static MPI_Comm comm = PETSC_COMM_SELF; + MPI_Comm pcomm = PetscObjectComm(reinterpret_cast(pc)); + if (pcomm != MPI_COMM_NULL) + comm = pcomm; + return comm; } - - void - PreconditionBase::create_pc() + PreconditionBase::create_pc_with_mat(const MatrixBase &matrix) { // only allow the creation of the // preconditioner once AssertThrow(pc == nullptr, StandardExceptions::ExcInvalidState()); - MPI_Comm comm = get_mpi_communicator(); - - PetscErrorCode ierr = PCCreate(comm, &pc); + MPI_Comm comm; + PetscErrorCode ierr = PetscObjectGetComm( + reinterpret_cast(static_cast(matrix)), &comm); AssertThrow(ierr == 0, ExcPETScError(ierr)); + create_pc_with_comm(comm); + # if DEAL_II_PETSC_VERSION_LT(3, 5, 0) ierr = PCSetOperators(pc, matrix, matrix, SAME_PRECONDITIONER); # else @@ -124,6 +122,13 @@ namespace PETScWrappers AssertThrow(ierr == 0, ExcPETScError(ierr)); } + void + PreconditionBase::create_pc_with_comm(const MPI_Comm &comm) + { + clear(); + PetscErrorCode ierr = PCCreate(comm, &pc); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + } const PC & PreconditionBase::get_pc() const @@ -132,16 +137,10 @@ namespace PETScWrappers } - PreconditionBase::operator Mat() const - { - return matrix; - } - - /* ----------------- PreconditionJacobi -------------------- */ PreconditionJacobi::PreconditionJacobi() - : PreconditionBase(MPI_COMM_NULL) + : PreconditionBase() {} @@ -189,23 +188,17 @@ namespace PETScWrappers { clear(); - mpi_communicator = matrix_.get_mpi_communicator(); - - matrix = static_cast(matrix_); additional_data = additional_data_; - create_pc(); + create_pc_with_mat(matrix_); initialize(); - - PetscErrorCode ierr = PCSetUp(pc); - AssertThrow(ierr == 0, ExcPETScError(ierr)); } /* ----------------- PreconditionBlockJacobi -------------------- */ PreconditionBlockJacobi::PreconditionBlockJacobi() - : PreconditionBase(MPI_COMM_NULL) + : PreconditionBase() {} PreconditionBlockJacobi::PreconditionBlockJacobi( @@ -251,23 +244,17 @@ namespace PETScWrappers { clear(); - mpi_communicator = matrix_.get_mpi_communicator(); - - matrix = static_cast(matrix_); additional_data = additional_data_; - create_pc(); + create_pc_with_mat(matrix_); initialize(); - - PetscErrorCode ierr = PCSetUp(pc); - AssertThrow(ierr == 0, ExcPETScError(ierr)); } /* ----------------- PreconditionSOR -------------------- */ PreconditionSOR::PreconditionSOR() - : PreconditionBase(MPI_COMM_NULL) + : PreconditionBase() {} @@ -292,12 +279,9 @@ namespace PETScWrappers { clear(); - mpi_communicator = matrix_.get_mpi_communicator(); - - matrix = static_cast(matrix_); additional_data = additional_data_; - create_pc(); + create_pc_with_mat(matrix_); PetscErrorCode ierr = PCSetType(pc, const_cast(PCSOR)); AssertThrow(ierr == 0, ExcPETScError(ierr)); @@ -308,16 +292,13 @@ namespace PETScWrappers ierr = PCSetFromOptions(pc); AssertThrow(ierr == 0, ExcPETScError(ierr)); - - ierr = PCSetUp(pc); - AssertThrow(ierr == 0, ExcPETScError(ierr)); } /* ----------------- PreconditionSSOR -------------------- */ PreconditionSSOR::PreconditionSSOR() - : PreconditionBase(MPI_COMM_NULL) + : PreconditionBase() {} @@ -342,12 +323,9 @@ namespace PETScWrappers { clear(); - mpi_communicator = matrix_.get_mpi_communicator(); - - matrix = static_cast(matrix_); additional_data = additional_data_; - create_pc(); + create_pc_with_mat(matrix_); PetscErrorCode ierr = PCSetType(pc, const_cast(PCSOR)); AssertThrow(ierr == 0, ExcPETScError(ierr)); @@ -362,16 +340,13 @@ namespace PETScWrappers ierr = PCSetFromOptions(pc); AssertThrow(ierr == 0, ExcPETScError(ierr)); - - ierr = PCSetUp(pc); - AssertThrow(ierr == 0, ExcPETScError(ierr)); } /* ----------------- PreconditionICC -------------------- */ PreconditionICC::PreconditionICC() - : PreconditionBase(MPI_COMM_NULL) + : PreconditionBase() {} @@ -396,12 +371,9 @@ namespace PETScWrappers { clear(); - mpi_communicator = matrix_.get_mpi_communicator(); - - matrix = static_cast(matrix_); additional_data = additional_data_; - create_pc(); + create_pc_with_mat(matrix_); PetscErrorCode ierr = PCSetType(pc, const_cast(PCICC)); AssertThrow(ierr == 0, ExcPETScError(ierr)); @@ -412,16 +384,13 @@ namespace PETScWrappers ierr = PCSetFromOptions(pc); AssertThrow(ierr == 0, ExcPETScError(ierr)); - - ierr = PCSetUp(pc); - AssertThrow(ierr == 0, ExcPETScError(ierr)); } /* ----------------- PreconditionILU -------------------- */ PreconditionILU::PreconditionILU() - : PreconditionBase(MPI_COMM_NULL) + : PreconditionBase() {} @@ -446,12 +415,9 @@ namespace PETScWrappers { clear(); - mpi_communicator = matrix_.get_mpi_communicator(); - - matrix = static_cast(matrix_); additional_data = additional_data_; - create_pc(); + create_pc_with_mat(matrix_); PetscErrorCode ierr = PCSetType(pc, const_cast(PCILU)); AssertThrow(ierr == 0, ExcPETScError(ierr)); @@ -462,9 +428,6 @@ namespace PETScWrappers ierr = PCSetFromOptions(pc); AssertThrow(ierr == 0, ExcPETScError(ierr)); - - ierr = PCSetUp(pc); - AssertThrow(ierr == 0, ExcPETScError(ierr)); } @@ -579,7 +542,7 @@ namespace PETScWrappers PreconditionBoomerAMG::PreconditionBoomerAMG() - : PreconditionBase(MPI_COMM_NULL) + : PreconditionBase() {} @@ -733,17 +696,11 @@ namespace PETScWrappers # ifdef DEAL_II_PETSC_WITH_HYPRE clear(); - mpi_communicator = matrix_.get_mpi_communicator(); - - matrix = static_cast(matrix_); additional_data = additional_data_; - create_pc(); + create_pc_with_mat(matrix_); initialize(); - PetscErrorCode ierr = PCSetUp(pc); - AssertThrow(ierr == 0, ExcPETScError(ierr)); - # else // DEAL_II_PETSC_WITH_HYPRE (void)matrix_; (void)additional_data_; @@ -772,7 +729,7 @@ namespace PETScWrappers PreconditionParaSails::PreconditionParaSails() - : PreconditionBase(MPI_COMM_NULL) + : PreconditionBase() {} @@ -792,13 +749,10 @@ namespace PETScWrappers { clear(); - mpi_communicator = matrix_.get_mpi_communicator(); - - matrix = static_cast(matrix_); additional_data = additional_data_; # ifdef DEAL_II_PETSC_WITH_HYPRE - create_pc(); + create_pc_with_mat(matrix_); PetscErrorCode ierr = PCSetType(pc, const_cast(PCHYPRE)); AssertThrow(ierr == 0, ExcPETScError(ierr)); @@ -861,9 +815,6 @@ namespace PETScWrappers ierr = PCSetFromOptions(pc); AssertThrow(ierr == 0, ExcPETScError(ierr)); - ierr = PCSetUp(pc); - AssertThrow(ierr == 0, ExcPETScError(ierr)); - # else // DEAL_II_PETSC_WITH_HYPRE (void)pc; Assert(false, @@ -876,7 +827,7 @@ namespace PETScWrappers /* ----------------- PreconditionNone ------------------------- */ PreconditionNone::PreconditionNone() - : PreconditionBase(MPI_COMM_NULL) + : PreconditionBase() {} @@ -895,21 +846,15 @@ namespace PETScWrappers { clear(); - mpi_communicator = matrix_.get_mpi_communicator(); - - matrix = static_cast(matrix_); additional_data = additional_data_; - create_pc(); + create_pc_with_mat(matrix_); PetscErrorCode ierr = PCSetType(pc, const_cast(PCNONE)); AssertThrow(ierr == 0, ExcPETScError(ierr)); ierr = PCSetFromOptions(pc); AssertThrow(ierr == 0, ExcPETScError(ierr)); - - ierr = PCSetUp(pc); - AssertThrow(ierr == 0, ExcPETScError(ierr)); } @@ -926,7 +871,7 @@ namespace PETScWrappers PreconditionLU::PreconditionLU() - : PreconditionBase(MPI_COMM_NULL) + : PreconditionBase() {} @@ -945,12 +890,9 @@ namespace PETScWrappers { clear(); - mpi_communicator = matrix_.get_mpi_communicator(); - - matrix = static_cast(matrix_); additional_data = additional_data_; - create_pc(); + create_pc_with_mat(matrix_); PetscErrorCode ierr = PCSetType(pc, const_cast(PCLU)); AssertThrow(ierr == 0, ExcPETScError(ierr)); @@ -967,9 +909,6 @@ namespace PETScWrappers ierr = PCSetFromOptions(pc); AssertThrow(ierr == 0, ExcPETScError(ierr)); - - ierr = PCSetUp(pc); - AssertThrow(ierr == 0, ExcPETScError(ierr)); } /* ----------------- PreconditionBDDC -------------------- */ @@ -992,7 +931,7 @@ namespace PETScWrappers template PreconditionBDDC::PreconditionBDDC() - : PreconditionBase(MPI_COMM_NULL) + : PreconditionBase() {} @@ -1034,8 +973,16 @@ namespace PETScWrappers // The matrix must be of IS type. We check for this to avoid the PETSc error // in order to suggest the correct matrix reinit method. { - MatType current_type; - ierr = MatGetType(matrix, ¤t_type); + MatType current_type; + Mat A, P; + PetscBool flg; + + ierr = PCGetOperators(pc, &A, &P); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + ierr = PCGetUseAmat(pc, &flg); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + + ierr = MatGetType(flg ? A : P, ¤t_type); AssertThrow(ierr == 0, ExcPETScError(ierr)); AssertThrow( strcmp(current_type, MATIS) == 0, @@ -1105,18 +1052,100 @@ namespace PETScWrappers { clear(); - mpi_communicator = matrix_.get_mpi_communicator(); - - matrix = static_cast(matrix_); additional_data = additional_data_; - create_pc(); + create_pc_with_mat(matrix_); initialize(); + } - PetscErrorCode ierr = PCSetUp(pc); + /* ----------------- PreconditionShell -------------------- */ + + PreconditionShell::PreconditionShell(const MatrixBase &matrix) + { + initialize(matrix); + } + + PreconditionShell::PreconditionShell(const MPI_Comm &comm) + { + initialize(comm); + } + + void + PreconditionShell::initialize(const MPI_Comm &comm) + { + PetscErrorCode ierr; + if (pc) + { + ierr = PCDestroy(&pc); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + } + create_pc_with_comm(comm); + + ierr = PCSetType(pc, PCSHELL); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + ierr = PCShellSetContext(pc, static_cast(this)); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + ierr = PCShellSetApply(pc, PreconditionShell::pcapply); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + ierr = PCShellSetApplyTranspose(pc, PreconditionShell::pcapply_transpose); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + ierr = PCShellSetName(pc, "deal.II user solve"); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + } + + void + PreconditionShell::initialize(const MatrixBase &matrix) + { + initialize(matrix.get_mpi_communicator()); + PetscErrorCode ierr; +# if DEAL_II_PETSC_VERSION_LT(3, 5, 0) + ierr = PCSetOperators(pc, matrix, matrix, DIFFERENT_NONZERO_PATTERN); +# else + ierr = PCSetOperators(pc, matrix, matrix); +# endif AssertThrow(ierr == 0, ExcPETScError(ierr)); } + int + PreconditionShell::pcapply(PC ppc, Vec x, Vec y) + { + void *ctx; + + PetscFunctionBeginUser; + PetscErrorCode ierr = PCShellGetContext(ppc, &ctx); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + + VectorBase src(x); + VectorBase dst(y); + auto user = static_cast(ctx); + if (user->apply) + user->apply(dst, src); + else + dst.sadd(0, src); + PetscFunctionReturn(0); + } + + int + PreconditionShell::pcapply_transpose(PC ppc, Vec x, Vec y) + { + void *ctx; + + PetscFunctionBeginUser; + PetscErrorCode ierr = PCShellGetContext(ppc, &ctx); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + + auto user = static_cast(ctx); + VectorBase src(x); + VectorBase dst(y); + if (user->applyT) + user->applyT(dst, src); + else if (user->apply) + user->apply(dst, src); + else + dst.sadd(0, src); + PetscFunctionReturn(0); + } + } // namespace PETScWrappers diff --git a/source/lac/petsc_solver.cc b/source/lac/petsc_solver.cc index 33b9e23292..8a4ec5c471 100644 --- a/source/lac/petsc_solver.cc +++ b/source/lac/petsc_solver.cc @@ -86,15 +86,12 @@ namespace PETScWrappers ierr = KSPSetPC(solver_data->ksp, preconditioner.get_pc()); AssertThrow(ierr == 0, ExcPETScError(ierr)); - // make sure the preconditioner has an associated matrix set - const Mat B = preconditioner; - AssertThrow(B != nullptr, - ExcMessage("PETSc preconditioner should have an " - "associated matrix set to be used in solver.")); - // setting the preconditioner overwrites the used matrices. // hence, we need to set the matrices after the preconditioner. - ierr = KSPSetOperators(solver_data->ksp, A, preconditioner); + Mat B; + ierr = PCGetOperators(preconditioner.get_pc(), nullptr, &B); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + ierr = KSPSetOperators(solver_data->ksp, A, B); AssertThrow(ierr == 0, ExcPETScError(ierr)); // then a convergence monitor diff --git a/tests/petsc/solver_03_mf.cc b/tests/petsc/solver_03_mf.cc index b162e33740..31a65e9ed1 100644 --- a/tests/petsc/solver_03_mf.cc +++ b/tests/petsc/solver_03_mf.cc @@ -62,5 +62,12 @@ main(int argc, char **argv) control.last_step(), 42, 44); + + u = 0.; + PETScWrappers::PreconditionShell preconditioner_user(A); + check_solver_within_range(solver.solve(A, u, f, preconditioner_user), + control.last_step(), + 42, + 44); } } diff --git a/tests/petsc/solver_03_mf.output b/tests/petsc/solver_03_mf.output index a6685358ae..636735f86a 100644 --- a/tests/petsc/solver_03_mf.output +++ b/tests/petsc/solver_03_mf.output @@ -2,3 +2,4 @@ DEAL::Size 32 Unknowns 961 DEAL::Solver type: N6dealii13PETScWrappers8SolverCGE DEAL::Solver stopped within 42 - 44 iterations +DEAL::Solver stopped within 42 - 44 iterations -- 2.39.5