petsc_precondition.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 2004 - 2017 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#include <deal.II/lac/petsc_precondition.h>

#ifdef DEAL_II_WITH_PETSC

#  include <deal.II/base/utilities.h>
#  include <deal.II/lac/exceptions.h>
#  include <deal.II/lac/petsc_compatibility.h>
#  include <deal.II/lac/petsc_matrix_base.h>
#  include <deal.II/lac/petsc_vector_base.h>
#  include <deal.II/lac/petsc_solver.h>
#  include <petscconf.h>
#  include <cmath>

DEAL_II_NAMESPACE_OPEN

namespace PETScWrappers
{
  PreconditionerBase::PreconditionerBase ()
    :
    pc(nullptr), matrix(nullptr)
  {}

  PreconditionerBase::~PreconditionerBase ()
  {
    try
      {
        clear();
      }
    catch (...)
      {}
  }

  void
  PreconditionerBase::clear ()
  {
    matrix = nullptr;

    if (pc!=nullptr)
      {
        PetscErrorCode ierr = PCDestroy(&pc);
        pc = nullptr;
        AssertThrow (ierr == 0, ExcPETScError(ierr));
      }
  }


  void
  PreconditionerBase::vmult (VectorBase       &dst,
                             const VectorBase &src) const
  {
    AssertThrow (pc != nullptr, StandardExceptions::ExcInvalidState ());

    const PetscErrorCode ierr = PCApply(pc, src, dst);
    AssertThrow (ierr == 0, ExcPETScError(ierr));
  }


  void
  PreconditionerBase::create_pc ()
  {
    // only allow the creation of the
    // preconditioner once
    AssertThrow (pc == nullptr, StandardExceptions::ExcInvalidState ());

    MPI_Comm comm;
    // this ugly cast is necessary because the
    // type Mat and PETScObject are
    // unrelated.
    PetscErrorCode ierr = PetscObjectGetComm(reinterpret_cast<PetscObject>(matrix),
                                             &comm);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = PCCreate(comm, &pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

#if DEAL_II_PETSC_VERSION_LT(3, 5, 0)
    ierr = PCSetOperators(pc, matrix, matrix, SAME_PRECONDITIONER);
#else
    ierr = PCSetOperators(pc, matrix, matrix);
#endif
    AssertThrow (ierr == 0, ExcPETScError(ierr));
  }


  const PC &
  PreconditionerBase::get_pc () const
  {
    return pc;
  }


  PreconditionerBase::operator Mat () const
  {
    return matrix;
  }


  /* ----------------- PreconditionJacobi -------------------- */
  PreconditionJacobi::PreconditionJacobi (const MPI_Comm comm,
                                          const AdditionalData &additional_data_)
  {
    additional_data = additional_data_;

    PetscErrorCode ierr = PCCreate(comm, &pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    initialize();
  }



  PreconditionJacobi::PreconditionJacobi (const MatrixBase     &matrix,
                                          const AdditionalData &additional_data)
  {
    initialize(matrix, additional_data);
  }

  void
  PreconditionJacobi::initialize()
  {
    AssertThrow (pc != nullptr, StandardExceptions::ExcInvalidState ());

    PetscErrorCode ierr = PCSetType (pc, const_cast<char *>(PCJACOBI));
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = PCSetFromOptions (pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));
  }

  void
  PreconditionJacobi::initialize (const MatrixBase     &matrix_,
                                  const AdditionalData &additional_data_)
  {
    clear ();

    matrix = static_cast<Mat>(matrix_);
    additional_data = additional_data_;

    create_pc();
    initialize();

    PetscErrorCode ierr = PCSetUp (pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));
  }


  /* ----------------- PreconditionBlockJacobi -------------------- */
  PreconditionBlockJacobi::PreconditionBlockJacobi (const MPI_Comm comm,
                                                    const AdditionalData &additional_data_)
  {
    additional_data = additional_data_;

    PetscErrorCode ierr = PCCreate(comm, &pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    initialize();
  }



  PreconditionBlockJacobi::
  PreconditionBlockJacobi (const MatrixBase     &matrix,
                           const AdditionalData &additional_data)
  {
    initialize(matrix, additional_data);
  }

  void
  PreconditionBlockJacobi::initialize()
  {
    PetscErrorCode ierr = PCSetType (pc, const_cast<char *>(PCBJACOBI));
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = PCSetFromOptions (pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));
  }


  void
  PreconditionBlockJacobi::initialize (const MatrixBase     &matrix_,
                                       const AdditionalData &additional_data_)
  {
    clear ();

    matrix = static_cast<Mat>(matrix_);
    additional_data = additional_data_;

    create_pc();
    initialize();

    PetscErrorCode ierr = PCSetUp (pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));
  }


  /* ----------------- PreconditionSOR -------------------- */

  PreconditionSOR::AdditionalData::
  AdditionalData (const double omega)
    :
    omega (omega)
  {}



  PreconditionSOR::PreconditionSOR (const MatrixBase     &matrix,
                                    const AdditionalData &additional_data)
  {
    initialize(matrix, additional_data);
  }


  void
  PreconditionSOR::initialize (const MatrixBase     &matrix_,
                               const AdditionalData &additional_data_)
  {
    clear ();

    matrix = static_cast<Mat>(matrix_);
    additional_data = additional_data_;

    create_pc();

    PetscErrorCode ierr = PCSetType (pc, const_cast<char *>(PCSOR));
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    // then set flags as given
    ierr = PCSORSetOmega (pc, additional_data.omega);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = PCSetFromOptions (pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = PCSetUp (pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));
  }


  /* ----------------- PreconditionSSOR -------------------- */

  PreconditionSSOR::AdditionalData::
  AdditionalData (const double omega)
    :
    omega (omega)
  {}



  PreconditionSSOR::PreconditionSSOR (const MatrixBase     &matrix,
                                      const AdditionalData &additional_data)
  {
    initialize(matrix, additional_data);
  }


  void
  PreconditionSSOR::initialize (const MatrixBase     &matrix_,
                                const AdditionalData &additional_data_)
  {
    clear ();

    matrix = static_cast<Mat>(matrix_);
    additional_data = additional_data_;

    create_pc();

    PetscErrorCode ierr = PCSetType (pc, const_cast<char *>(PCSOR));
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    // then set flags as given
    ierr = PCSORSetOmega (pc, additional_data.omega);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    // convert SOR to SSOR
    ierr = PCSORSetSymmetric (pc, SOR_SYMMETRIC_SWEEP);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = PCSetFromOptions (pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = PCSetUp (pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));
  }


  /* ----------------- PreconditionEisenstat -------------------- */

  PreconditionEisenstat::AdditionalData::
  AdditionalData (const double omega)
    :
    omega (omega)
  {}



  PreconditionEisenstat::PreconditionEisenstat (const MatrixBase     &matrix,
                                                const AdditionalData &additional_data)
  {
    initialize(matrix, additional_data);
  }


  void
  PreconditionEisenstat::initialize (const MatrixBase     &matrix_,
                                     const AdditionalData &additional_data_)
  {
    clear ();

    matrix = static_cast<Mat>(matrix_);
    additional_data = additional_data_;

    create_pc();

    PetscErrorCode ierr = PCSetType (pc, const_cast<char *>(PCEISENSTAT));
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    // then set flags as given
    ierr = PCEisenstatSetOmega (pc, additional_data.omega);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = PCSetFromOptions (pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = PCSetUp (pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));
  }


  /* ----------------- PreconditionICC -------------------- */


  PreconditionICC::AdditionalData::
  AdditionalData (const unsigned int levels)
    :
    levels (levels)
  {}



  PreconditionICC::PreconditionICC (const MatrixBase     &matrix,
                                    const AdditionalData &additional_data)
  {
    initialize(matrix, additional_data);
  }


  void
  PreconditionICC::initialize (const MatrixBase     &matrix_,
                               const AdditionalData &additional_data_)
  {
    clear ();

    matrix = static_cast<Mat>(matrix_);
    additional_data = additional_data_;

    create_pc();

    PetscErrorCode ierr = PCSetType (pc, const_cast<char *>(PCICC));
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    // then set flags
    ierr = PCFactorSetLevels (pc, additional_data.levels);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = PCSetFromOptions (pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = PCSetUp (pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));
  }


  /* ----------------- PreconditionILU -------------------- */

  PreconditionILU::AdditionalData::
  AdditionalData (const unsigned int levels)
    :
    levels (levels)
  {}



  PreconditionILU::PreconditionILU (const MatrixBase     &matrix,
                                    const AdditionalData &additional_data)
  {
    initialize(matrix, additional_data);
  }


  void
  PreconditionILU::initialize (const MatrixBase     &matrix_,
                               const AdditionalData &additional_data_)
  {
    clear ();

    matrix = static_cast<Mat>(matrix_);
    additional_data = additional_data_;

    create_pc();

    PetscErrorCode ierr = PCSetType (pc, const_cast<char *>(PCILU));
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    // then set flags
    ierr = PCFactorSetLevels (pc, additional_data.levels);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = PCSetFromOptions (pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = PCSetUp (pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));
  }


  /* ----------------- PreconditionBoomerAMG -------------------- */

  PreconditionBoomerAMG::AdditionalData::
  AdditionalData(const bool symmetric_operator,
                 const double strong_threshold,
                 const double max_row_sum,
                 const unsigned int aggressive_coarsening_num_levels,
                 const bool output_details
                )
    :
    symmetric_operator(symmetric_operator),
    strong_threshold(strong_threshold),
    max_row_sum(max_row_sum),
    aggressive_coarsening_num_levels(aggressive_coarsening_num_levels),
    output_details(output_details)
  {}



  PreconditionBoomerAMG::PreconditionBoomerAMG (const MPI_Comm comm,
                                                const AdditionalData &additional_data_)
  {
    additional_data = additional_data_;

    PetscErrorCode ierr = PCCreate(comm, &pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

#ifdef DEAL_II_PETSC_WITH_HYPRE
    initialize();
#else // DEAL_II_PETSC_WITH_HYPRE
    (void)pc;
    Assert (false,
            ExcMessage ("Your PETSc installation does not include a copy of "
                        "the hypre package necessary for this preconditioner."));
#endif
  }


  PreconditionBoomerAMG::PreconditionBoomerAMG (const MatrixBase     &matrix,
                                                const AdditionalData &additional_data)
  {
    initialize(matrix, additional_data);
  }

  void
  PreconditionBoomerAMG::initialize ()
  {
#ifdef DEAL_II_PETSC_WITH_HYPRE
    PetscErrorCode ierr = PCSetType (pc, const_cast<char *>(PCHYPRE));
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = PCHYPRESetType(pc, "boomeramg");
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    if (additional_data.output_details)
      {
        set_option_value("-pc_hypre_boomeramg_print_statistics", "1");
      }

    set_option_value("-pc_hypre_boomeramg_agg_nl",
                     Utilities::to_string(additional_data.aggressive_coarsening_num_levels));

    std::stringstream ssStream;
    ssStream << additional_data.max_row_sum;
    set_option_value("-pc_hypre_boomeramg_max_row_sum", ssStream.str());

    ssStream.str(""); // empty the stringstream
    ssStream << additional_data.strong_threshold;
    set_option_value("-pc_hypre_boomeramg_strong_threshold", ssStream.str());

    if (additional_data.symmetric_operator)
      {
        set_option_value("-pc_hypre_boomeramg_relax_type_up", "symmetric-SOR/Jacobi");
        set_option_value("-pc_hypre_boomeramg_relax_type_down", "symmetric-SOR/Jacobi");
        set_option_value("-pc_hypre_boomeramg_relax_type_coarse", "Gaussian-elimination");
      }
    else
      {
        set_option_value("-pc_hypre_boomeramg_relax_type_up", "SOR/Jacobi");
        set_option_value("-pc_hypre_boomeramg_relax_type_down", "SOR/Jacobi");
        set_option_value("-pc_hypre_boomeramg_relax_type_coarse", "Gaussian-elimination");
      }

    ierr = PCSetFromOptions (pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));
#else
    Assert (false,
            ExcMessage ("Your PETSc installation does not include a copy of "
                        "the hypre package necessary for this preconditioner."));
#endif
  }

  void
  PreconditionBoomerAMG::initialize (const MatrixBase     &matrix_,
                                     const AdditionalData &additional_data_)
  {
#ifdef DEAL_II_PETSC_WITH_HYPRE
    clear ();

    matrix = static_cast<Mat>(matrix_);
    additional_data = additional_data_;

    create_pc();
    initialize ();

    PetscErrorCode ierr = PCSetUp (pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

#else // DEAL_II_PETSC_WITH_HYPRE
    (void)matrix_;
    (void)additional_data_;
    Assert (false,
            ExcMessage ("Your PETSc installation does not include a copy of "
                        "the hypre package necessary for this preconditioner."));
#endif
  }


  /* ----------------- PreconditionParaSails -------------------- */

  PreconditionParaSails::AdditionalData::
  AdditionalData(const unsigned int symmetric,
                 const unsigned int n_levels,
                 const double threshold,
                 const double filter,
                 const bool output_details)
    :
    symmetric(symmetric),
    n_levels(n_levels),
    threshold(threshold),
    filter(filter),
    output_details(output_details)
  {}



  PreconditionParaSails::PreconditionParaSails (const MatrixBase     &matrix,
                                                const AdditionalData &additional_data)
  {
    initialize(matrix, additional_data);
  }


  void
  PreconditionParaSails::initialize (const MatrixBase     &matrix_,
                                     const AdditionalData &additional_data_)
  {
    clear ();

    matrix = static_cast<Mat>(matrix_);
    additional_data = additional_data_;

#ifdef DEAL_II_PETSC_WITH_HYPRE
    create_pc();

    PetscErrorCode ierr = PCSetType (pc, const_cast<char *>(PCHYPRE));
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = PCHYPRESetType(pc, "parasails");
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    if (additional_data.output_details)
      {
        set_option_value("-pc_hypre_parasails_logging","1");
      }

    Assert ((additional_data.symmetric == 0 ||
             additional_data.symmetric == 1 ||
             additional_data.symmetric == 2),
            ExcMessage("ParaSails parameter symmetric can only be equal to 0, 1, 2!"));

    std::stringstream ssStream;

    switch (additional_data.symmetric)
      {
      case 0:
      {
        ssStream << "nonsymmetric";
        break;
      }

      case 1:
      {
        ssStream << "SPD";
        break;
      }

      case 2:
      {
        ssStream << "nonsymmetric,SPD";
        break;
      }

      default:
        Assert (false,
                ExcMessage("ParaSails parameter symmetric can only be equal to 0, 1, 2!"));
      };

    set_option_value("-pc_hypre_parasails_sym",ssStream.str());

    set_option_value ("-pc_hypre_parasails_nlevels",
                      Utilities::to_string(additional_data.n_levels));

    ssStream.str(""); // empty the stringstream
    ssStream << additional_data.threshold;
    set_option_value("-pc_hypre_parasails_thresh", ssStream.str());

    ssStream.str(""); // empty the stringstream
    ssStream << additional_data.filter;
    set_option_value("-pc_hypre_parasails_filter", ssStream.str());

    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,
            ExcMessage ("Your PETSc installation does not include a copy of "
                        "the hypre package necessary for this preconditioner."));
#endif
  }


  /* ----------------- PreconditionNone ------------------------- */

  PreconditionNone::PreconditionNone (const MatrixBase     &matrix,
                                      const AdditionalData &additional_data)
  {
    initialize (matrix, additional_data);
  }


  void
  PreconditionNone::initialize (const MatrixBase     &matrix_,
                                const AdditionalData &additional_data_)
  {
    clear ();

    matrix = static_cast<Mat>(matrix_);
    additional_data = additional_data_;

    create_pc();

    PetscErrorCode ierr = PCSetType (pc, const_cast<char *>(PCNONE));
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = PCSetFromOptions (pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = PCSetUp (pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));
  }


  /* ----------------- PreconditionLU -------------------- */

  PreconditionLU::AdditionalData::
  AdditionalData (const double pivoting,
                  const double zero_pivot,
                  const double damping)
    :
    pivoting (pivoting),
    zero_pivot (zero_pivot),
    damping (damping)
  {}



  PreconditionLU::PreconditionLU (const MatrixBase     &matrix,
                                  const AdditionalData &additional_data)
  {
    initialize(matrix, additional_data);
  }


  void
  PreconditionLU::initialize (const MatrixBase     &matrix_,
                              const AdditionalData &additional_data_)
  {
    clear ();

    matrix = static_cast<Mat>(matrix_);
    additional_data = additional_data_;

    create_pc();

    PetscErrorCode ierr = PCSetType (pc, const_cast<char *>(PCLU));
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    // set flags as given
    ierr = PCFactorSetColumnPivot (pc, additional_data.pivoting);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = PCFactorSetZeroPivot (pc, additional_data.zero_pivot);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = PCFactorSetShiftAmount (pc, additional_data.damping);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = PCSetFromOptions (pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));

    ierr = PCSetUp (pc);
    AssertThrow (ierr == 0, ExcPETScError(ierr));
  }

}

DEAL_II_NAMESPACE_CLOSE

#endif // DEAL_II_WITH_PETSC