From 75ff08eac426685793c843e98dd760c0610c9be6 Mon Sep 17 00:00:00 2001 From: Stefano Zampini Date: Tue, 18 Apr 2023 16:48:50 +0300 Subject: [PATCH] PETScWrappers::Solver Bring the class into the new century --- include/deal.II/lac/petsc_compatibility.h | 19 - include/deal.II/lac/petsc_solver.h | 127 ++--- source/lac/petsc_solver.cc | 491 ++++++++------------ source/lac/slepc_spectral_transformation.cc | 2 +- 4 files changed, 244 insertions(+), 395 deletions(-) diff --git a/include/deal.II/lac/petsc_compatibility.h b/include/deal.II/lac/petsc_compatibility.h index 96bd230629..46d61f7822 100644 --- a/include/deal.II/lac/petsc_compatibility.h +++ b/include/deal.II/lac/petsc_compatibility.h @@ -53,25 +53,6 @@ namespace PETScWrappers - /** - * Destroy a Krylov Subspace (KSP) PETSc solver. This function wraps - * KSPDestroy with a version check (the signature of this function changed - * in PETSc 3.2.0). - * - * @warning Since the primary intent of this function is to enable RAII - * semantics in the PETSc wrappers, this function will not throw an - * exception if an error occurs, but instead just returns the error code - * given by MatDestroy. - */ - inline PetscErrorCode - destroy_krylov_solver(KSP &krylov_solver) - { - // PETSc will check whether or not matrix is nullptr. - return KSPDestroy(&krylov_solver); - } - - - /** * Set a PETSc matrix option. This function wraps MatSetOption with a * version check. diff --git a/include/deal.II/lac/petsc_solver.h b/include/deal.II/lac/petsc_solver.h index a687966709..eaee6fca71 100644 --- a/include/deal.II/lac/petsc_solver.h +++ b/include/deal.II/lac/petsc_solver.h @@ -21,13 +21,13 @@ #ifdef DEAL_II_WITH_PETSC +# include + # include # include # include -# include - # ifdef DEAL_II_WITH_SLEPC # include # endif @@ -77,11 +77,15 @@ namespace PETScWrappers * linear solver and can be obtained from the documentation and manual pages. * + * @note Command line options relative to convergence tolerances are ignored + * when the solver is instantiated with a SolverControl. + * * @note Repeated calls to solve() on a solver object with a Preconditioner * must be used with care. The preconditioner is initialized in the first * call to solve() and subsequent calls reuse the solver and preconditioner * object. This is done for performance reasons. The solver and - * preconditioner can be reset by calling reset(). + * preconditioner can be reset by calling reset() or by using the command + * line option "-ksp_reuse_preconditioner false". * * @ingroup PETScWrappers */ @@ -89,14 +93,20 @@ namespace PETScWrappers { public: /** - * Constructor. + * Default constructor without SolverControl. The resulting solver will + * use PETSc's default convergence testing routines. + */ + SolverBase(); + + /** + * Constructor with a SolverControl instance. */ SolverBase(SolverControl &cn); /** * Destructor. */ - virtual ~SolverBase() = default; + virtual ~SolverBase(); /** * Solve the linear system Ax=b. Depending on the information @@ -111,7 +121,6 @@ namespace PETScWrappers const VectorBase & b, const PreconditionBase &preconditioner); - /** * Resets the contained preconditioner and solver object. See class * description for more details. @@ -119,7 +128,6 @@ namespace PETScWrappers virtual void reset(); - /** * Sets a prefix name for the solver object. Useful when customizing the * PETSc KSP object with command-line options. @@ -127,7 +135,6 @@ namespace PETScWrappers void set_prefix(const std::string &prefix); - /** * Access to object that controls convergence. */ @@ -141,21 +148,54 @@ namespace PETScWrappers void initialize(const PreconditionBase &preconditioner); + /** + * Return the PETSc KSP object. + */ + KSP + petsc_ksp(); + + /** + * Conversion operator to gain access to the underlying PETSc type. If you + * do this, you cut this class off some information it may need, so this + * conversion operator should only be used if you know what you do. + */ + operator KSP() const; + protected: + /** + * The PETSc KSP object. + */ + KSP ksp; + /** * Reference to the object that controls convergence of the iterative * solver. In fact, for these PETSc wrappers, PETSc does so itself, but we * copy the data from this object before starting the solution process, * and copy the data back into it afterwards. */ - SolverControl &solver_control; + SmartPointer solver_control; + + /** + * Utility to create the KSP object and attach convergence test. + */ + void + initialize_ksp_with_comm(const MPI_Comm &comm); /** * %Function that takes a Krylov Subspace Solver context object, and sets * the type of solver that is requested by the derived class. */ virtual void - set_solver_type(KSP &ksp) const = 0; + set_solver_type(KSP &ksp) const; + + /** + * Utility to use deal.II convergence testing. + * + * This call changes the convergence criterion when the instance of the + * class has a SolverControl object associated. + */ + void + perhaps_set_convergence_test() const; /** * Solver prefix name to qualify options specific to the PETSc KSP object @@ -179,41 +219,6 @@ namespace PETScWrappers KSPConvergedReason *reason, void * solver_control); - /** - * A structure that contains the PETSc solver and preconditioner objects. - * This object is preserved between subsequent calls to the solver if the - * same preconditioner is used as in the previous solver step. This may - * save some computation time, if setting up a preconditioner is - * expensive, such as in the case of an ILU for example. - * - * The actual declaration of this class is complicated by the fact that - * PETSc changed its solver interface completely and incompatibly between - * versions 2.1.6 and 2.2.0 :-( - * - * Objects of this type are explicitly created, but are destroyed when the - * surrounding solver object goes out of scope, or when we assign a new - * value to the pointer to this object. The respective *Destroy functions - * are therefore written into the destructor of this object, even though - * the object does not have a constructor. - */ - struct SolverData - { - /** - * Destructor - */ - ~SolverData(); - - /** - * Object for Krylov subspace solvers. - */ - KSP ksp; - }; - - /** - * Pointer to an object that stores the solver context. This is recreated - * in the main solver routine if necessary. - */ - std::unique_ptr solver_data; # ifdef DEAL_II_WITH_SLEPC // Make the transformation class a friend, since it needs to set the KSP @@ -945,38 +950,6 @@ namespace PETScWrappers virtual void set_solver_type(KSP &ksp) const override; - private: - /** - * A function that is used in PETSc as a callback to check convergence. It - * takes the information provided from PETSc and checks it against - * deal.II's own SolverControl objects to see if convergence has been - * reached. - */ - static PetscErrorCode - convergence_test(KSP ksp, - const PetscInt iteration, - const PetscReal residual_norm, - KSPConvergedReason *reason, - void * solver_control); - - /** - * A structure that contains the PETSc solver and preconditioner objects. - * Since the solve member function in the base is not used here, the - * private SolverData struct located in the base could not be used either. - */ - struct SolverDataMUMPS - { - /** - * Destructor - */ - ~SolverDataMUMPS(); - - KSP ksp; - PC pc; - }; - - std::unique_ptr solver_data; - /** * Flag specifies whether matrix being factorized is symmetric or not. It * influences the type of the used preconditioner (PCLU or PCCHOLESKY) diff --git a/source/lac/petsc_solver.cc b/source/lac/petsc_solver.cc index 834a3b0c9f..060525f3e0 100644 --- a/source/lac/petsc_solver.cc +++ b/source/lac/petsc_solver.cc @@ -26,25 +26,57 @@ # include # include -# include - -# include -# include +// Shorthand notation for PETSc error codes. +# define AssertPETSc(code) \ + do \ + { \ + PetscErrorCode ierr = (code); \ + AssertThrow(ierr == 0, ExcPETScError(ierr)); \ + } \ + while (0) DEAL_II_NAMESPACE_OPEN namespace PETScWrappers { - SolverBase::SolverData::~SolverData() + SolverBase::SolverBase() + : ksp(nullptr) + , solver_control(nullptr) + {} + + + SolverBase::SolverBase(SolverControl &cn) + : ksp(nullptr) + , solver_control(&cn) + {} + + + + void + SolverBase::set_solver_type(KSP &) const + {} + + + + SolverBase::~SolverBase() { - destroy_krylov_solver(ksp); + AssertPETSc(KSPDestroy(&ksp)); } - SolverBase::SolverBase(SolverControl &cn) - : solver_control(cn) - {} + KSP + SolverBase::petsc_ksp() + { + return ksp; + } + + + + SolverBase::operator KSP() const + { + return ksp; + } @@ -54,86 +86,54 @@ namespace PETScWrappers const VectorBase & b, const PreconditionBase &preconditioner) { - /* - TODO: PETSc duplicates communicators, so this does not work (you put - MPI_COMM_SELF in, but get something other out when you ask PETSc for the - communicator. This mainly fails due to the MatrixFree classes, that can not - ask PETSc for a communicator. //Timo Heister - Assert(A.get_mpi_communicator()==mpi_communicator, ExcMessage("PETSc Solver - and Matrix need to use the same MPI_Comm.")); - Assert(x.get_mpi_communicator()==mpi_communicator, ExcMessage("PETSc Solver - and Vector need to use the same MPI_Comm.")); - Assert(b.get_mpi_communicator()==mpi_communicator, ExcMessage("PETSc Solver - and Vector need to use the same MPI_Comm.")); - */ - // first create a solver object if this // is necessary - if (solver_data.get() == nullptr) + if (ksp == nullptr) { - solver_data = std::make_unique(); - - PetscErrorCode ierr = - KSPCreate(A.get_mpi_communicator(), &solver_data->ksp); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + initialize_ksp_with_comm(A.get_mpi_communicator()); // let derived classes set the solver // type, and the preconditioning // object set the type of // preconditioner - set_solver_type(solver_data->ksp); - - ierr = KSPSetPC(solver_data->ksp, preconditioner.get_pc()); - AssertThrow(ierr == 0, ExcPETScError(ierr)); - - // setting the preconditioner overwrites the used matrices. - // hence, we need to set the matrices after the 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 - // function. that function simply - // checks with the solver_control - // object we have in this object for - // convergence - ierr = KSPSetConvergenceTest(solver_data->ksp, - &convergence_test, - reinterpret_cast(&solver_control), - nullptr); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + set_solver_type(ksp); + + AssertPETSc(KSPSetPC(ksp, preconditioner.get_pc())); + + /* + * by default we set up the preconditioner only once. + * this can be overriden by command line. + */ + AssertPETSc(KSPSetReusePreconditioner(ksp, PETSC_TRUE)); } + // setting the preconditioner overwrites the used matrices. + // hence, we need to set the matrices after the preconditioner. + Mat B; + AssertPETSc(KSPGetOperators(ksp, nullptr, &B)); + AssertPETSc(KSPSetOperators(ksp, A, B)); + // set the command line option prefix name - PetscErrorCode ierr = - KSPSetOptionsPrefix(solver_data->ksp, prefix_name.c_str()); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetOptionsPrefix(ksp, prefix_name.c_str())); // set the command line options provided // by the user to override the defaults - ierr = KSPSetFromOptions(solver_data->ksp); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetFromOptions(ksp)); // then do the real work: set up solver // internal data and solve the // system. - ierr = KSPSetUp(solver_data->ksp); - AssertThrow(ierr == 0, ExcPETScError(ierr)); - - ierr = KSPSolve(solver_data->ksp, b, x); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetUp(ksp)); - // do not destroy solver object - // solver_data.reset (); + AssertPETSc(KSPSolve(ksp, b, x)); // in case of failure: throw // exception - if (solver_control.last_check() != SolverControl::success) + if (solver_control && + solver_control->last_check() != SolverControl::success) AssertThrow(false, - SolverControl::NoConvergence(solver_control.last_step(), - solver_control.last_value())); + SolverControl::NoConvergence(solver_control->last_step(), + solver_control->last_value())); // otherwise exit as normal } @@ -148,14 +148,18 @@ namespace PETScWrappers void SolverBase::reset() { - solver_data.reset(); + AssertPETSc(KSPDestroy(&ksp)); } SolverControl & SolverBase::control() const { - return solver_control; + AssertThrow( + solver_control, + ExcMessage( + "You need to create the solver with a SolverControl object if you want to call the function that returns it.")); + return *solver_control; } @@ -179,7 +183,7 @@ namespace PETScWrappers break; case ::dealii::SolverControl::success: - *reason = static_cast(1); + *reason = KSP_CONVERGED_RTOL; break; case ::dealii::SolverControl::failure: @@ -200,39 +204,44 @@ namespace PETScWrappers void - SolverBase::initialize(const PreconditionBase &preconditioner) + SolverBase::initialize_ksp_with_comm(const MPI_Comm &comm) + { + // Create the PETSc KSP object + AssertPETSc(KSPCreate(comm, &ksp)); + + // then a convergence monitor + // function that simply + // checks with the solver_control + // object we have in this object for + // convergence + perhaps_set_convergence_test(); + } + + + + void + SolverBase::perhaps_set_convergence_test() const { - PetscErrorCode ierr; + if (ksp && solver_control) + AssertPETSc( + KSPSetConvergenceTest(ksp, &convergence_test, solver_control, nullptr)); + } - solver_data = std::make_unique(); - ierr = KSPCreate(preconditioner.get_mpi_communicator(), &solver_data->ksp); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + void + SolverBase::initialize(const PreconditionBase &preconditioner) + { + initialize_ksp_with_comm(preconditioner.get_mpi_communicator()); // let derived classes set the solver // type, and the preconditioning // object set the type of // preconditioner - set_solver_type(solver_data->ksp); - - ierr = KSPSetPC(solver_data->ksp, preconditioner.get_pc()); - AssertThrow(ierr == 0, ExcPETScError(ierr)); - - // then a convergence monitor - // function. that function simply - // checks with the solver_control - // object we have in this object for - // convergence - ierr = KSPSetConvergenceTest(solver_data->ksp, - &convergence_test, - reinterpret_cast(&solver_control), - nullptr); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + set_solver_type(ksp); // set the command line options provided // by the user to override the defaults - ierr = KSPSetFromOptions(solver_data->ksp); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetFromOptions(ksp)); } @@ -263,18 +272,15 @@ namespace PETScWrappers void SolverRichardson::set_solver_type(KSP &ksp) const { - PetscErrorCode ierr = KSPSetType(ksp, KSPRICHARDSON); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetType(ksp, KSPRICHARDSON)); // set the damping factor from the data - ierr = KSPRichardsonSetScale(ksp, additional_data.omega); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPRichardsonSetScale(ksp, additional_data.omega)); // in the deal.II solvers, we always // honor the initial guess in the // solution vector. do so here as well: - ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE)); // Hand over the absolute // tolerance and the maximum @@ -289,12 +295,11 @@ namespace PETScWrappers // the Richardson iteration, // where no residual is // available. - ierr = KSPSetTolerances(ksp, - PETSC_DEFAULT, - this->solver_control.tolerance(), - PETSC_DEFAULT, - this->solver_control.max_steps() + 1); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetTolerances(ksp, + PETSC_DEFAULT, + this->solver_control->tolerance(), + PETSC_DEFAULT, + this->solver_control->max_steps() + 1)); } @@ -317,14 +322,12 @@ namespace PETScWrappers void SolverChebychev::set_solver_type(KSP &ksp) const { - PetscErrorCode ierr = KSPSetType(ksp, KSPCHEBYSHEV); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetType(ksp, KSPCHEBYSHEV)); // in the deal.II solvers, we always // honor the initial guess in the // solution vector. do so here as well: - ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE)); } @@ -346,14 +349,12 @@ namespace PETScWrappers void SolverCG::set_solver_type(KSP &ksp) const { - PetscErrorCode ierr = KSPSetType(ksp, KSPCG); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetType(ksp, KSPCG)); // in the deal.II solvers, we always // honor the initial guess in the // solution vector. do so here as well: - ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE)); } @@ -375,14 +376,12 @@ namespace PETScWrappers void SolverBiCG::set_solver_type(KSP &ksp) const { - PetscErrorCode ierr = KSPSetType(ksp, KSPBICG); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetType(ksp, KSPBICG)); // in the deal.II solvers, we always // honor the initial guess in the // solution vector. do so here as well: - ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE)); } @@ -413,24 +412,20 @@ namespace PETScWrappers void SolverGMRES::set_solver_type(KSP &ksp) const { - PetscErrorCode ierr = KSPSetType(ksp, KSPGMRES); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetType(ksp, KSPGMRES)); - ierr = KSPGMRESSetRestart(ksp, additional_data.restart_parameter); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPGMRESSetRestart(ksp, additional_data.restart_parameter)); // Set preconditioning side to right if (additional_data.right_preconditioning) { - ierr = KSPSetPCSide(ksp, PC_RIGHT); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetPCSide(ksp, PC_RIGHT)); } // in the deal.II solvers, we always // honor the initial guess in the // solution vector. do so here as well: - ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE)); } @@ -452,14 +447,12 @@ namespace PETScWrappers void SolverBicgstab::set_solver_type(KSP &ksp) const { - PetscErrorCode ierr = KSPSetType(ksp, KSPBCGS); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetType(ksp, KSPBCGS)); // in the deal.II solvers, we always // honor the initial guess in the // solution vector. do so here as well: - ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE)); } @@ -481,14 +474,12 @@ namespace PETScWrappers void SolverCGS::set_solver_type(KSP &ksp) const { - PetscErrorCode ierr = KSPSetType(ksp, KSPCGS); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetType(ksp, KSPCGS)); // in the deal.II solvers, we always // honor the initial guess in the // solution vector. do so here as well: - ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE)); } @@ -510,14 +501,12 @@ namespace PETScWrappers void SolverTFQMR::set_solver_type(KSP &ksp) const { - PetscErrorCode ierr = KSPSetType(ksp, KSPTFQMR); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetType(ksp, KSPTFQMR)); // in the deal.II solvers, we always // honor the initial guess in the // solution vector. do so here as well: - ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE)); } @@ -539,14 +528,12 @@ namespace PETScWrappers void SolverTCQMR::set_solver_type(KSP &ksp) const { - PetscErrorCode ierr = KSPSetType(ksp, KSPTCQMR); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetType(ksp, KSPTCQMR)); // in the deal.II solvers, we always // honor the initial guess in the // solution vector. do so here as well: - ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE)); } @@ -568,14 +555,12 @@ namespace PETScWrappers void SolverCR::set_solver_type(KSP &ksp) const { - PetscErrorCode ierr = KSPSetType(ksp, KSPCR); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetType(ksp, KSPCR)); // in the deal.II solvers, we always // honor the initial guess in the // solution vector. do so here as well: - ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE)); } @@ -599,14 +584,19 @@ namespace PETScWrappers void SolverLSQR::set_solver_type(KSP &ksp) const { - PetscErrorCode ierr = KSPSetType(ksp, KSPLSQR); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetType(ksp, KSPLSQR)); // in the deal.II solvers, we always // honor the initial guess in the // solution vector. do so here as well: - ierr = KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE)); + + // The KSPLSQR implementation overwrites the user-defined + // convergence test at creation (i.e. KSPSetType) time. + // This is probably a bad design decision in PETSc. + // Anyway, here we make sure we use our own convergence + // test. + perhaps_set_convergence_test(); } @@ -628,8 +618,7 @@ namespace PETScWrappers void SolverPreOnly::set_solver_type(KSP &ksp) const { - PetscErrorCode ierr = KSPSetType(ksp, KSPPREONLY); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetType(ksp, KSPPREONLY)); // The KSPPREONLY solver of // PETSc never calls the convergence @@ -639,13 +628,12 @@ namespace PETScWrappers // is set to some nice values, which // guarantee a nice result at the end // of the solution process. - solver_control.check(1, 0.0); + solver_control->check(1, 0.0); // Using the PREONLY solver with // a nonzero initial guess leads // PETSc to produce some error messages. - ierr = KSPSetInitialGuessNonzero(ksp, PETSC_FALSE); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_FALSE)); } @@ -668,14 +656,6 @@ namespace PETScWrappers - SparseDirectMUMPS::SolverDataMUMPS::~SolverDataMUMPS() - { - destroy_krylov_solver(ksp); - // the 'pc' object is owned by the 'ksp' object, and consequently - // does not have to be destroyed explicitly here - } - - void SparseDirectMUMPS::set_solver_type(KSP &ksp) const { @@ -684,8 +664,7 @@ namespace PETScWrappers * preconditioner. Its use is due to SparseDirectMUMPS being a direct * (rather than iterative) solver */ - PetscErrorCode ierr = KSPSetType(ksp, KSPPREONLY); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetType(ksp, KSPPREONLY)); /* * The KSPPREONLY solver of PETSc never calls the convergence monitor, @@ -693,14 +672,13 @@ namespace PETScWrappers * SolverControl status is set to some nice values, which guarantee a * nice result at the end of the solution process. */ - solver_control.check(1, 0.0); + solver_control->check(1, 0.0); /* * Using a PREONLY solver with a nonzero initial guess leads PETSc to * produce some error messages. */ - ierr = KSPSetInitialGuessNonzero(ksp, PETSC_FALSE); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetInitialGuessNonzero(ksp, PETSC_FALSE)); } void @@ -709,154 +687,110 @@ namespace PETScWrappers const VectorBase &b) { # ifdef DEAL_II_PETSC_WITH_MUMPS - /* - * factorization matrix to be obtained from MUMPS - */ - Mat F; - - /* - * setting MUMPS integer control parameters ICNTL to be passed to - * MUMPS. Setting entry 7 of MUMPS ICNTL array (of size 40) to a value - * of 2. This sets use of Approximate Minimum Fill (AMF) - */ - PetscInt ival = 2, icntl = 7; - /* - * number of iterations to solution (should be 1) for a direct solver - */ - PetscInt its; - /* - * norm of residual - */ - PetscReal rnorm; - /* * creating a solver object if this is necessary */ - if (solver_data == nullptr) + if (ksp == nullptr) { - solver_data = std::make_unique(); + initialize_ksp_with_comm(A.get_mpi_communicator()); /* - * creates the default KSP context and puts it in the location - * solver_data->ksp + * setting the solver type */ - PetscErrorCode ierr = - KSPCreate(A.get_mpi_communicator(), &solver_data->ksp); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + set_solver_type(ksp); /* * set the matrices involved. the last argument is irrelevant here, * since we use the solver only once anyway */ - ierr = KSPSetOperators(solver_data->ksp, A, A); - AssertThrow(ierr == 0, ExcPETScError(ierr)); - - /* - * setting the solver type - */ - set_solver_type(solver_data->ksp); + AssertPETSc(KSPSetOperators(ksp, A, A)); /* * getting the associated preconditioner context */ - ierr = KSPGetPC(solver_data->ksp, &solver_data->pc); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + PC pc; + AssertPETSc(KSPGetPC(ksp, &pc)); /* * build PETSc PC for particular PCLU or PCCHOLESKY preconditioner * depending on whether the symmetric mode has been set */ if (symmetric_mode) - ierr = PCSetType(solver_data->pc, PCCHOLESKY); + AssertPETSc(PCSetType(pc, PCCHOLESKY)); else - ierr = PCSetType(solver_data->pc, PCLU); - AssertThrow(ierr == 0, ExcPETScError(ierr)); - - /* - * convergence monitor function that checks with the solver_control - * object for convergence - */ - ierr = KSPSetConvergenceTest(solver_data->ksp, - &convergence_test, - reinterpret_cast(&solver_control), - nullptr); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(PCSetType(pc, PCLU)); - /* - * set the software that is to be used to perform the lu - * factorization here we start to see differences with the base - * class solve function - */ + /* + * set the software that is to be used to perform the lu + * factorization here we start to see differences with the base + * class solve function + */ # if DEAL_II_PETSC_VERSION_LT(3, 9, 0) - ierr = PCFactorSetMatSolverPackage(solver_data->pc, MATSOLVERMUMPS); + AssertPETSc(PCFactorSetMatSolverPackage(pc, MATSOLVERMUMPS)); # else - ierr = PCFactorSetMatSolverType(solver_data->pc, MATSOLVERMUMPS); + AssertPETSc(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS)); # endif - AssertThrow(ierr == 0, ExcPETScError(ierr)); /* * set up the package to call for the factorization */ # if DEAL_II_PETSC_VERSION_LT(3, 9, 0) - ierr = PCFactorSetUpMatSolverPackage(solver_data->pc); + AssertPETSc(PCFactorSetUpMatSolverPackage(pc)); # else - ierr = PCFactorSetUpMatSolverType(solver_data->pc); + AssertPETSc(PCFactorSetUpMatSolverType(pc)); # endif - AssertThrow(ierr == 0, ExcPETScError(ierr)); - - /* - * get the factored matrix F from the preconditioner context. This - * routine is valid only for LU, ILU, Cholesky, and incomplete - * Cholesky - */ - ierr = PCFactorGetMatrix(solver_data->pc, &F); - AssertThrow(ierr == 0, ExcPETScError(ierr)); /* - * Passing the control parameters to MUMPS + * get the factored matrix F from the preconditioner context. */ - ierr = MatMumpsSetIcntl(F, icntl, ival); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + Mat F; + AssertPETSc(PCFactorGetMatrix(pc, &F)); /* - * set the command line option prefix name + * pass control parameters to MUMPS. + * Setting entry 7 of MUMPS ICNTL array to a value + * of 2. This sets use of Approximate Minimum Fill (AMF) */ - ierr = KSPSetOptionsPrefix(solver_data->ksp, prefix_name.c_str()); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(MatMumpsSetIcntl(F, 7, 2)); /* - * set the command line options provided by the user to override - * the defaults + * by default we set up the preconditioner only once. + * this can be overriden by command line. */ - ierr = KSPSetFromOptions(solver_data->ksp); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSetReusePreconditioner(ksp, PETSC_TRUE)); } + /* + * set the matrices involved. the last argument is irrelevant here, + * since we use the solver only once anyway + */ + AssertPETSc(KSPSetOperators(ksp, A, A)); + + /* + * set the command line option prefix name + */ + AssertPETSc(KSPSetOptionsPrefix(ksp, prefix_name.c_str())); + + /* + * set the command line options provided by the user to override + * the defaults + */ + AssertPETSc(KSPSetFromOptions(ksp)); + /* * solve the linear system */ - PetscErrorCode ierr = KSPSolve(solver_data->ksp, b, x); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + AssertPETSc(KSPSolve(ksp, b, x)); /* * in case of failure throw exception */ - if (solver_control.last_check() != SolverControl::success) + if (solver_control && + solver_control->last_check() != SolverControl::success) { AssertThrow(false, - SolverControl::NoConvergence(solver_control.last_step(), - solver_control.last_value())); - } - else - { - /* - * obtain convergence information. obtain the number of iterations - * and residual norm - */ - ierr = KSPGetIterationNumber(solver_data->ksp, &its); - AssertThrow(ierr == 0, ExcPETScError(ierr)); - ierr = KSPGetResidualNorm(solver_data->ksp, &rnorm); - AssertThrow(ierr == 0, ExcPETScError(ierr)); + SolverControl::NoConvergence(solver_control->last_step(), + solver_control->last_value())); } # else // DEAL_II_PETSC_WITH_MUMPS @@ -877,45 +811,6 @@ namespace PETScWrappers - PetscErrorCode - SparseDirectMUMPS::convergence_test(KSP /*ksp*/, - const PetscInt iteration, - const PetscReal residual_norm, - KSPConvergedReason *reason, - void * solver_control_x) - { - SolverControl &solver_control = - *reinterpret_cast(solver_control_x); - - const SolverControl::State state = - solver_control.check(iteration, residual_norm); - - switch (state) - { - case ::dealii::SolverControl::iterate: - *reason = KSP_CONVERGED_ITERATING; - break; - - case ::dealii::SolverControl::success: - *reason = static_cast(1); - break; - - case ::dealii::SolverControl::failure: - if (solver_control.last_step() > solver_control.max_steps()) - *reason = KSP_DIVERGED_ITS; - else - *reason = KSP_DIVERGED_DTOL; - break; - - default: - Assert(false, ExcNotImplemented()); - } - - return 0; - } - - - void SparseDirectMUMPS::set_symmetric_mode(const bool flag) { diff --git a/source/lac/slepc_spectral_transformation.cc b/source/lac/slepc_spectral_transformation.cc index dc68c3c8c7..3a8cf0321b 100644 --- a/source/lac/slepc_spectral_transformation.cc +++ b/source/lac/slepc_spectral_transformation.cc @@ -55,7 +55,7 @@ namespace SLEPcWrappers void TransformationBase::set_solver(const PETScWrappers::SolverBase &solver) { - PetscErrorCode ierr = STSetKSP(st, solver.solver_data->ksp); + PetscErrorCode ierr = STSetKSP(st, solver); AssertThrow(ierr == 0, SolverBase::ExcSLEPcError(ierr)); } -- 2.39.5