inverse_operator.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2022 by Luca Heltai
//
// This file is part of the FSI-suite platform, based on the deal.II library.
//
// The FSI-suite platform 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 3.0 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 FSI-suite platform
// distribution.
//
// ---------------------------------------------------------------------
 
#ifndef inverse_opeator_h
#define inverse_opeator_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/algorithms/general_data_storage.h>
 
#include <deal.II/base/parameter_acceptor.h>
 
#include <deal.II/lac/linear_operator_tools.h>
#include <deal.II/lac/solver.h>
#include <deal.II/lac/solver_bicgstab.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_control.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/solver_minres.h>
#include <deal.II/lac/solver_qmrs.h>
#include <deal.II/lac/solver_richardson.h>
 
namespace ParsedLAC
{
  enum class SolverControlType
  {
    tolerance              = 1 << 0, 
    consecutive_iterations = 1 << 1, 
    iteration_number       = 1 << 2, 
    reduction              = 1 << 3, 
  };
 
  class InverseOperator : public dealii::ParameterAcceptor
  {
  public:
    InverseOperator(
      const std::string       &section_name   = "",
      const std::string       &default_solver = "cg",
      const SolverControlType &control_type   = SolverControlType::tolerance,
      const unsigned int       max_iterations = 1000,
      const double             tolerance      = 1e-12,
      const double             reduction      = 1e-6,
      const unsigned int       consecutive_iterations = 2,
      const bool              &log_history            = false,
      const bool              &log_result             = false);
 
    template <typename Domain,
              typename Payload,
              typename PreconditionerType,
              typename Range = Domain>
    dealii::LinearOperator<Domain, Range, Payload>
    operator()(const dealii::LinearOperator<Range, Domain, Payload> &op,
               const PreconditionerType                             &prec,
               const double       abs_tol = 0.0,
               const std::string &prefix  = "") const;
 
    template <typename Range,
              typename Payload =
                dealii::internal::LinearOperatorImplementation::EmptyPayload,
              typename Domain = Range,
              typename MatrixType,
              typename PreconditionerType>
    dealii::LinearOperator<Range, Domain, Payload>
    solver(const MatrixType         &op,
           const PreconditionerType &prec,
           const double              abs_tol = 0.0,
           const std::string        &prefix  = "") const;
 
 
    template <typename MatrixType,
              typename PreconditionerType,
              typename VectorType>
    void
    solve(const MatrixType         &matrix,
          const PreconditionerType &preconditioner,
          const VectorType         &src,
          VectorType               &dst,
          const double              abs_tol = 0.0) const;
 
 
    std::string
    get_solver_name() const;
 
    std::unique_ptr<dealii::SolverControl>
    setup_new_solver_control(const double abs_tol = 0.0) const;
 
    template <typename Range, typename Payload>
    std::shared_ptr<dealii::SolverBase<Range>>
    setup_new_solver(const double abs_tol = 0.0) const;
 
  private:
    SolverControlType control_type;
 
    mutable std::unique_ptr<dealii::SolverControl> control;
 
    std::string solver_name;
 
    unsigned int max_iterations;
 
    unsigned int consecutive_iterations;
 
    double tolerance;
 
    double reduction;
 
    bool log_history;
 
    bool log_result;
 
    mutable dealii::GeneralDataStorage storage;
  };
 
 
#ifndef DOXYGEN
  // ============================================================
  // Explicit template instantiation
  // ============================================================
  template <typename MatrixType,
            typename PreconditionerType,
            typename VectorType>
  void
  InverseOperator::solve(const MatrixType         &matrix,
                         const PreconditionerType &preconditioner,
                         const VectorType         &src,
                         VectorType               &dst,
                         const double              abs_tol) const
  {
    control = setup_new_solver_control(abs_tol);
    if (solver_name == "cg")
      {
        dealii::SolverCG<VectorType> solver(*control);
        solver.solve(matrix, dst, src, preconditioner);
      }
    else if (solver_name == "bicgstab")
      {
        dealii::SolverBicgstab<VectorType> solver(*control);
        solver.solve(matrix, dst, src, preconditioner);
      }
    else if (solver_name == "gmres")
      {
        dealii::SolverGMRES<VectorType> solver(*control);
        solver.solve(matrix, dst, src, preconditioner);
      }
    else if (solver_name == "fgmres")
      {
        dealii::SolverFGMRES<VectorType> solver(*control);
        solver.solve(matrix, dst, src, preconditioner);
      }
    else if (solver_name == "minres")
      {
        dealii::SolverMinRes<VectorType> solver(*control);
        solver.solve(matrix, dst, src, preconditioner);
      }
    else if (solver_name == "qmrs")
      {
        dealii::SolverQMRS<VectorType> solver(*control);
        solver.solve(matrix, dst, src, preconditioner);
      }
    else if (solver_name == "richardson")
      {
        dealii::SolverRichardson<VectorType> solver(*control);
        solver.solve(matrix, dst, src, preconditioner);
      }
    else
      {
        Assert(false,
               dealii::ExcInternalError("Solver should not be unknonw."));
      }
  }
 
 
  template <typename Range,
            typename Payload,
            typename Domain,
            typename MatrixType,
            typename PreconditionerType>
  dealii::LinearOperator<Range, Domain, Payload>
  InverseOperator::solver(const MatrixType         &op,
                          const PreconditionerType &prec,
                          const double              abs_tol,
                          const std::string        &prefix) const
  {
    control = setup_new_solver_control(abs_tol);
 
    // Make sure the solver itself is left around until the object is
    // destroyed we need a general storage class, since we have no
    // idea what types are used when calling this function.
    using SolverType = std::shared_ptr<dealii::SolverBase<Range>>;
    auto &solver     = storage.template get_or_add_object_with_name<SolverType>(
      prefix + "solver");
 
    dealii::LinearOperator<Range, Domain, Payload> inverse;
 
    auto initialize_solver = [&](auto *s) {
      solver.reset(s);
      inverse = dealii::inverse_operator(op, *s, prec);
    };
 
    if (solver_name == "cg")
      {
        initialize_solver(new dealii::SolverCG<Range>(*control));
      }
    else if (solver_name == "bicgstab")
      {
        initialize_solver(new dealii::SolverBicgstab<Range>(*control));
      }
    else if (solver_name == "gmres")
      {
        initialize_solver(new dealii::SolverGMRES<Range>(*control));
      }
    else if (solver_name == "fgmres")
      {
        initialize_solver(new dealii::SolverFGMRES<Range>(*control));
      }
    else if (solver_name == "minres")
      {
        initialize_solver(new dealii::SolverMinRes<Range>(*control));
      }
    else if (solver_name == "qmrs")
      {
        initialize_solver(new dealii::SolverQMRS<Range>(*control));
      }
    else if (solver_name == "richardson")
      {
        initialize_solver(new dealii::SolverRichardson<Range>(*control));
      }
    else
      {
        Assert(false,
               dealii::ExcInternalError("Solver should not be unknonw."));
      }
    return inverse;
  }
 
 
 
  template <typename Domain,
            typename Payload,
            typename PreconditionerType,
            typename Range>
  dealii::LinearOperator<Domain, Range, Payload>
  InverseOperator::operator()(
    const dealii::LinearOperator<Range, Domain, Payload> &op,
    const PreconditionerType                             &prec,
    const double                                          abs_tol,
    const std::string                                    &prefix) const
  {
    return solver<Range, Payload>(op, prec, abs_tol, prefix);
  }
#endif
} // namespace ParsedLAC
 
 
#endif