petsc_precondition.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2004 - 2026 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#ifndef dealii_petsc_precondition_h
#define dealii_petsc_precondition_h
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/enable_observer_pointer.h>
#include <deal.II/base/mpi_stub.h>
#include <deal.II/base/point.h>
 
#ifdef DEAL_II_WITH_PETSC
 
#  include <deal.II/lac/exceptions.h>
 
#  include <petscpc.h>
 
#  include <functional>
 
#endif // DEAL_II_WITH_PETSC
 
DEAL_II_NAMESPACE_OPEN
 
#ifdef DEAL_II_WITH_PETSC
namespace PETScWrappers
{
  // forward declarations
#  ifndef DOXYGEN
  class MatrixBase;
  class VectorBase;
#  endif
 
  class PreconditionBase : public EnableObserverPointer
  {
  public:
    explicit PreconditionBase(const MPI_Comm mpi_communicator);
 
    PreconditionBase();
 
    virtual ~PreconditionBase();
 
    void
    clear();
 
    void
    vmult(VectorBase &dst, const VectorBase &src) const;
 
    void
    Tvmult(VectorBase &dst, const VectorBase &src) const;
 
    void
    setup();
 
    const PC &
    get_pc() const;
 
    MPI_Comm
    get_mpi_communicator() const;
 
  protected:
    PC pc;
 
    void
    create_pc_with_mat(const MatrixBase &);
 
    void
    create_pc_with_comm(const MPI_Comm);
  };
 
 
 
  class PreconditionJacobi : public PreconditionBase
  {
  public:
    struct AdditionalData
    {};
 
    PreconditionJacobi();
 
    PreconditionJacobi(
      const MatrixBase     &matrix,
      const AdditionalData &additional_data = AdditionalData());
 
    PreconditionJacobi(
      const MPI_Comm        communicator,
      const AdditionalData &additional_data = AdditionalData());
 
    void
    initialize(const MatrixBase     &matrix,
               const AdditionalData &additional_data = AdditionalData());
 
  protected:
    AdditionalData additional_data;
 
    void
    initialize();
  };
 
 
 
  class PreconditionBlockJacobi : public PreconditionBase
  {
  public:
    struct AdditionalData
    {};
 
    PreconditionBlockJacobi();
 
    PreconditionBlockJacobi(
      const MatrixBase     &matrix,
      const AdditionalData &additional_data = AdditionalData());
 
    PreconditionBlockJacobi(
      const MPI_Comm        communicator,
      const AdditionalData &additional_data = AdditionalData());
 
 
    void
    initialize(const MatrixBase     &matrix,
               const AdditionalData &additional_data = AdditionalData());
 
  protected:
    AdditionalData additional_data;
 
    void
    initialize();
  };
 
 
 
  class PreconditionSOR : public PreconditionBase
  {
  public:
    struct AdditionalData
    {
      AdditionalData(const double omega = 1);
 
      double omega;
    };
 
    PreconditionSOR();
 
    PreconditionSOR(const MatrixBase     &matrix,
                    const AdditionalData &additional_data = AdditionalData());
 
    void
    initialize(const MatrixBase     &matrix,
               const AdditionalData &additional_data = AdditionalData());
 
  protected:
    AdditionalData additional_data;
  };
 
 
 
  class PreconditionSSOR : public PreconditionBase
  {
  public:
    struct AdditionalData
    {
      explicit AdditionalData(const double omega = 1);
 
      double omega;
    };
 
    PreconditionSSOR();
 
    PreconditionSSOR(const MatrixBase     &matrix,
                     const AdditionalData &additional_data = AdditionalData());
 
    void
    initialize(const MatrixBase     &matrix,
               const AdditionalData &additional_data = AdditionalData());
 
  protected:
    AdditionalData additional_data;
  };
 
  class PreconditionICC : public PreconditionBase
  {
  public:
    struct AdditionalData
    {
      AdditionalData(const unsigned int levels = 0);
 
      unsigned int levels;
    };
 
    PreconditionICC();
 
    PreconditionICC(const MatrixBase     &matrix,
                    const AdditionalData &additional_data = AdditionalData());
 
    void
    initialize(const MatrixBase     &matrix,
               const AdditionalData &additional_data = AdditionalData());
 
  protected:
    AdditionalData additional_data;
  };
 
 
 
  class PreconditionILU : public PreconditionBase
  {
  public:
    struct AdditionalData
    {
      AdditionalData(const unsigned int levels = 0);
 
      unsigned int levels;
    };
 
    PreconditionILU();
 
    PreconditionILU(const MatrixBase     &matrix,
                    const AdditionalData &additional_data = AdditionalData());
 
    void
    initialize(const MatrixBase     &matrix,
               const AdditionalData &additional_data = AdditionalData());
 
  protected:
    AdditionalData additional_data;
  };
 
 
 
  class PreconditionLU : public PreconditionBase
  {
  public:
    struct AdditionalData
    {
      AdditionalData(const double pivoting   = 1.e-6,
                     const double zero_pivot = 1.e-12,
                     const double damping    = 0.0);
 
      double pivoting;
 
      double zero_pivot;
 
      double damping;
    };
 
    PreconditionLU();
 
    PreconditionLU(const MatrixBase     &matrix,
                   const AdditionalData &additional_data = AdditionalData());
 
    void
    initialize(const MatrixBase     &matrix,
               const AdditionalData &additional_data = AdditionalData());
 
  protected:
    AdditionalData additional_data;
  };
 
 
 
  class PreconditionBoomerAMG : public PreconditionBase
  {
  public:
    struct AdditionalData
    {
      enum class RelaxationType
      {
        Jacobi,
        sequentialGaussSeidel,
        seqboundaryGaussSeidel,
        SORJacobi,
        backwardSORJacobi,
        symmetricSORJacobi,
        l1scaledSORJacobi,
        GaussianElimination,
        l1GaussSeidel,
        backwardl1GaussSeidel,
        CG,
        Chebyshev,
        FCFJacobi,
        l1scaledJacobi,
        None
      };
 
      AdditionalData(
        const bool           symmetric_operator               = false,
        const double         strong_threshold                 = 0.25,
        const double         max_row_sum                      = 0.9,
        const unsigned int   aggressive_coarsening_num_levels = 0,
        const bool           output_details                   = false,
        const RelaxationType relaxation_type_up   = RelaxationType::SORJacobi,
        const RelaxationType relaxation_type_down = RelaxationType::SORJacobi,
        const RelaxationType relaxation_type_coarse =
          RelaxationType::GaussianElimination,
        const unsigned int n_sweeps_coarse = 1,
        const double       tol             = 0.0,
        const unsigned int max_iter        = 1,
        const bool         w_cycle         = false);
 
      bool symmetric_operator;
 
      double strong_threshold;
 
      double max_row_sum;
 
      unsigned int aggressive_coarsening_num_levels;
 
      bool output_details;
 
      RelaxationType relaxation_type_up;
 
      RelaxationType relaxation_type_down;
 
      RelaxationType relaxation_type_coarse;
 
      unsigned int n_sweeps_coarse;
 
      double tol;
 
      unsigned int max_iter;
 
      bool w_cycle;
    };
 
    PreconditionBoomerAMG();
 
    PreconditionBoomerAMG(
      const MatrixBase     &matrix,
      const AdditionalData &additional_data = AdditionalData());
 
    PreconditionBoomerAMG(
      const MPI_Comm        communicator,
      const AdditionalData &additional_data = AdditionalData());
 
 
    void
    initialize(const MatrixBase     &matrix,
               const AdditionalData &additional_data = AdditionalData());
 
  protected:
    AdditionalData additional_data;
 
    void
    initialize();
  };
 
 
 
  class PreconditionParaSails : public PreconditionBase
  {
  public:
    struct AdditionalData
    {
      AdditionalData(const unsigned int symmetric      = 1,
                     const unsigned int n_levels       = 1,
                     const double       threshold      = 0.1,
                     const double       filter         = 0.05,
                     const bool         output_details = false);
 
      unsigned int symmetric;
 
      unsigned int n_levels;
 
      double threshold;
 
      double filter;
 
      bool output_details;
    };
 
 
 
    PreconditionParaSails();
 
    PreconditionParaSails(
      const MatrixBase     &matrix,
      const AdditionalData &additional_data = AdditionalData());
 
    void
    initialize(const MatrixBase     &matrix,
               const AdditionalData &additional_data = AdditionalData());
 
  private:
    AdditionalData additional_data;
  };
 
 
 
  class PreconditionNone : public PreconditionBase
  {
  public:
    struct AdditionalData
    {};
 
    PreconditionNone();
 
    PreconditionNone(const MatrixBase     &matrix,
                     const AdditionalData &additional_data = AdditionalData());
 
    void
    initialize(const MatrixBase     &matrix,
               const AdditionalData &additional_data = AdditionalData());
 
  private:
    AdditionalData additional_data;
  };
 
  template <int dim>
  class PreconditionBDDC : public PreconditionBase
  {
  public:
    struct AdditionalData
    {
      AdditionalData(const bool                    use_vertices = true,
                     const bool                    use_edges    = false,
                     const bool                    use_faces    = false,
                     const bool                    symmetric    = false,
                     const std::vector<Point<dim>> coords       = {});
 
      bool use_vertices;
 
      bool use_edges;
 
      bool use_faces;
 
      bool symmetric;
 
      std::vector<Point<dim>> coords;
    };
 
    PreconditionBDDC();
 
    PreconditionBDDC(const MatrixBase     &matrix,
                     const AdditionalData &additional_data = AdditionalData());
 
    PreconditionBDDC(const MPI_Comm        communicator,
                     const AdditionalData &additional_data = AdditionalData());
 
    void
    initialize(const MatrixBase     &matrix,
               const AdditionalData &additional_data = AdditionalData());
 
  protected:
    AdditionalData additional_data;
 
    void
    initialize();
  };
 
 
 
  class PreconditionASM : public PreconditionBase
  {
  public:
    enum class ASMType
    {
      basic,
      restrict,
      interpolate,
      none
    };
 
    struct AdditionalData
    {
      struct LocalSubdomains
      {
        LocalSubdomains(
          const std::vector<std::vector<PetscInt>> &subdomain_indices = {},
          const std::vector<std::vector<PetscInt>> &local_indices     = {})
        {
          this->subdomain_indices = subdomain_indices;
          this->local_indices     = local_indices;
        }
 
        std::vector<std::vector<PetscInt>> subdomain_indices;
 
        std::vector<std::vector<PetscInt>> local_indices;
      };
 
      AdditionalData(
        const unsigned int     overlap          = 1,
        const ASMType          type             = ASMType::restrict,
        const LocalSubdomains &local_subdomains = LocalSubdomains({}, {}));
 
      unsigned int overlap;
 
      ASMType type;
 
      LocalSubdomains local_subdomains;
    };
 
    PreconditionASM();
 
    PreconditionASM(const MatrixBase     &matrix,
                    const AdditionalData &additional_data = AdditionalData());
 
    void
    initialize(const MatrixBase     &matrix,
               const AdditionalData &additional_data = AdditionalData());
 
 
  protected:
    AdditionalData additional_data;
  };
 
 
 
  class PreconditionShell : public PreconditionBase
  {
  public:
    PreconditionShell() = default;
 
    PreconditionShell(const MatrixBase &matrix);
 
    PreconditionShell(const MPI_Comm communicator);
 
    std::function<void(VectorBase &dst, const VectorBase &src)> vmult;
 
    std::function<void(VectorBase &dst, const VectorBase &src)> vmultT;
 
  protected:
    void
    initialize(const MPI_Comm comm);
 
    void
    initialize(const MatrixBase &matrix);
 
  private:
    static PetscErrorCode
    pcapply(PC pc, Vec src, Vec dst);
 
    static PetscErrorCode
    pcapply_transpose(PC pc, Vec src, Vec dst);
 
    static PetscErrorCode
    pcsetup(PC pc);
  };
} // namespace PETScWrappers
 
#endif // DEAL_II_WITH_PETSC
 
DEAL_II_NAMESPACE_CLOSE
 
#endif