PETScKrylov.hh

Go to the documentation of this file.

#pragma once
 
// external libraries
#include <deal.II/base/config.h>
 
#ifdef DEAL_II_WITH_PETSC
 
#include <deal.II/lac/petsc_precondition.h>
#include <deal.II/lac/petsc_solver.h>
#include <deal.II/lac/solver_control.h>
 
// standard library
#include <stdexcept>
 
// DiFfRG
#include <DiFfRG/timestepping/linear_solver/abstract_linear_solver.hh>
 
namespace DiFfRG
{
  namespace internal
  {
#ifdef DEAL_II_PETSC_WITH_HYPRE
    using PETScDefaultPreconditioner = dealii::PETScWrappers::PreconditionBoomerAMG;
#else
    using PETScDefaultPreconditioner = dealii::PETScWrappers::PreconditionBlockJacobi;
#endif
 
    template <typename PreconditionerType> typename PreconditionerType::AdditionalData precondition_settings()
    {
      return typename PreconditionerType::AdditionalData();
    }
 
#ifdef DEAL_II_PETSC_WITH_HYPRE
    template <>
    inline dealii::PETScWrappers::PreconditionBoomerAMG::AdditionalData
    precondition_settings<dealii::PETScWrappers::PreconditionBoomerAMG>()
    {
      using AMG = dealii::PETScWrappers::PreconditionBoomerAMG;
      AMG::AdditionalData data;
      data.relaxation_type_coarse = AMG::AdditionalData::RelaxationType::SORJacobi;
      return data;
    }
#endif
  } // namespace internal
 
  template <typename SparseMatrixType, typename VectorType,
            typename PreconditionerType = internal::PETScDefaultPreconditioner>
  class PETScKrylov : public AbstractLinearSolver<SparseMatrixType, VectorType>
  {
  public:
    static constexpr bool performs_factorization = false;
 
    PETScKrylov() : matrix(nullptr) {}
 
    void init(const SparseMatrixType &matrix)
    {
      this->matrix = &matrix;
      // Rebuilding the preconditioner on every Jacobian is deliberate: IDA hands us a new
      // Jacobian only when it has decided the old one is stale, so reusing a preconditioner
      // across init() calls would be reusing it across exactly the changes it exists to track.
      preconditioner.initialize(matrix, internal::precondition_settings<PreconditionerType>());
    }
 
    bool invert() { return false; }
 
    int solve(const VectorType &src, VectorType &dst, const double tol)
    {
      if (!matrix) throw std::runtime_error("PETScKrylov::solve: matrix not initialized");
 
      // Same iteration cap as the serial GMRES wrapper. src.size() is the *global* size for a
      // PETSc vector, which is what we want here -- the cap should not depend on rank count.
      dealii::SolverControl solver_control(std::max<std::size_t>(1000, src.size() / 10), tol);
      dealii::PETScWrappers::SolverGMRES solver(solver_control);
 
      try {
        solver.solve(*matrix, dst, src, preconditioner);
      } catch (std::exception &e) {
        std::cerr << "PETSc GMRES linear solver failed: " << e.what() << std::endl;
        throw;
      }
 
      return solver_control.last_step();
    }
 
  private:
    const SparseMatrixType *matrix;
    PreconditionerType preconditioner;
  };
} // namespace DiFfRG
 
#endif // DEAL_II_WITH_PETSC