ISTLSolverCreator.hpp

#pragma once
 
#include <memory>
 
#include <dune/common/version.hh>
 
#include <dune/common/ftraits.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/solvertype.hh>
 
#if HAVE_SUITESPARSE_UMFPACK
#include <dune/istl/umfpack.hh>
#endif
#if HAVE_SUITESPARSE_LDL
#include <dune/istl/ldl.hh>
#endif
#if HAVE_SUITESPARSE_SPQR
#include <dune/istl/spqr.hh>
#endif
#if HAVE_SUPERLU
#include <dune/istl/superlu.hh>
#endif
 
#include <amdis/CreatorMap.hpp>
#include <amdis/Environment.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/Output.hpp>
 
#include <amdis/linearalgebra/istl/ISTLPreconCreator.hpp>
#include <amdis/linearalgebra/istl/SolverWrapper.hpp>
#include <amdis/linearalgebra/istl/Traits.hpp>
 
namespace AMDiS
{
  namespace tag
  {
    template <class Traits>
    struct solver {};
 
    // forward declaration
    template <class Traits>
    struct preconditioner;
  }
 
 
  template <class Traits>
  class ISTLSolverCreatorBase
      : public CreatorInterface<tag::solver<Traits>>
  {
    using X = typename Traits::X;
    using Y = typename Traits::Y;
 
    std::unique_ptr<tag::solver<Traits>> create() final { return {}; };
 
  public:
    virtual void init(std::string const& prefix)
    {
      if (Environment::mpiRank() == 0)
        Parameters::get(prefix + "->info", info_);
    }
 
    virtual std::unique_ptr<Dune::InverseOperator<X,Y>>
    createSolver(typename Traits::M const& A, typename Traits::Comm const& comm) const = 0;
 
  protected:
    int info_ = 0;
  };
 
 
 
  template <class Traits>
  class ISTLIterativeSolverCreatorBase
      : public ISTLSolverCreatorBase<Traits>
  {
    using Super = ISTLSolverCreatorBase<Traits>;
 
    using real_type = typename Dune::FieldTraits<typename Traits::M::field_type>::real_type;
 
  public:
    void init(std::string const& prefix) override
    {
      Super::init(prefix);
 
      maxIter_ = Parameters::get<int>(prefix + "->maxit").value_or(
                 Parameters::get<int>(prefix + "->max iteration").value_or(maxIter_));
      rTol_ = Parameters::get<real_type>(prefix + "->reduction").value_or(
              Parameters::get<real_type>(prefix + "->relative tolerance").value_or(rTol_));
 
      std::string precon = "default";
      Parameters::get(prefix + "->precon", precon);
 
      using CreatorMap = CreatorMap<tag::preconditioner<Traits>>;
      auto* creator = CreatorMap::get(precon, prefix + "->precon");
      preconCreator_ = dynamic_cast<ISTLPreconCreatorBase<Traits>*>(creator);
      assert(preconCreator_ != nullptr);
      preconCreator_->init(prefix + "->precon");
    }
 
  protected:
    template <class Solver, class... Args>
    auto create_impl(typename Traits::M const& mat, typename Traits::Comm const& comm, Args&&... args) const
    {
      auto cat = Dune::SolverCategory::category(comm);
      auto sp = Traits::ScalProdCreator::create(cat, comm);
      auto linOp = Traits::LinOpCreator::create(cat, mat, comm);
 
      assert(preconCreator_ != nullptr);
      auto precon = preconCreator_->createPrecon(mat, comm);
      return std::make_unique<IterativeSolverWrapper<Solver>>(
        std::move(linOp), std::move(sp), std::move(precon), FWD(args)...);
    }
 
  protected:
    int maxIter_ = 500;
    real_type rTol_ = 1.e-6;
    ISTLPreconCreatorBase<Traits>* preconCreator_ = nullptr;
  };
 
 
 
  template <class Solver, class Traits>
  class IterativeSolverCreator
      : public ISTLIterativeSolverCreatorBase<Traits>
  {
    using Super = ISTLIterativeSolverCreatorBase<Traits>;
    using Interface = typename Traits::Solver;
 
  public:
    std::unique_ptr<Interface>
    createSolver(typename Traits::M const& mat, typename Traits::Comm const& comm) const override
    {
      return this->template create_impl<Solver>(mat, comm, this->rTol_, this->maxIter_, this->info_);
    }
  };
 
 
 
  template <class Solver, class Traits>
  struct GMResSolverCreator
      : public ISTLIterativeSolverCreatorBase<Traits>
  {
    using Super = ISTLIterativeSolverCreatorBase<Traits>;
    using Interface = typename Traits::Solver;
 
  public:
    void init(std::string const& prefix) override
    {
      Super::init(prefix);
      Parameters::get(prefix + "->restart", restart_);
    }
 
    std::unique_ptr<Interface>
    createSolver(typename Traits::M const& mat, typename Traits::Comm const& comm) const override
    {
      return this->template create_impl<Solver>(mat, comm, this->rTol_, restart_, this->maxIter_, this->info_);
    }
 
  private:
    int restart_ = 30;
  };
 
  template <class X, class Y, class Traits>
  struct IterativeSolverCreator<Dune::RestartedGMResSolver<X,Y>, Traits>
      : public GMResSolverCreator<Dune::RestartedGMResSolver<X,Y>, Traits>
  {};
 
  template <class X, class Y, class Traits>
  struct IterativeSolverCreator<Dune::RestartedFlexibleGMResSolver<X,Y>, Traits>
      : public GMResSolverCreator<Dune::RestartedFlexibleGMResSolver<X,Y>, Traits>
  {};
 
 
 
  template <class Solver, class Traits>
  struct PCGSolverCreator
      : public ISTLIterativeSolverCreatorBase<Traits>
  {
    using Super = ISTLIterativeSolverCreatorBase<Traits>;
    using Interface = typename Traits::Solver;
 
  public:
    void init(std::string const& prefix) override
    {
      Super::init(prefix);
      Parameters::get(prefix + "->restart", restart_);
    }
 
    std::unique_ptr<Interface>
    createSolver(typename Traits::M const& mat, typename Traits::Comm const& comm) const override
    {
      return this->template create_impl<Solver>(mat, comm, this->rTol_, this->maxIter_, this->info_, restart_);
    }
 
  private:
    int restart_ = 30;
  };
 
  template <class X, class Traits>
  struct IterativeSolverCreator<Dune::GeneralizedPCGSolver<X>, Traits>
      : public PCGSolverCreator<Dune::GeneralizedPCGSolver<X>, Traits>
  {};
 
  template <class X, class Traits>
  struct IterativeSolverCreator<Dune::RestartedFCGSolver<X>, Traits>
      : public PCGSolverCreator<Dune::RestartedFCGSolver<X>, Traits>
  {};
 
  template <class X, class Traits>
  struct IterativeSolverCreator<Dune::CompleteFCGSolver<X>, Traits>
      : public PCGSolverCreator<Dune::CompleteFCGSolver<X>, Traits>
  {};
 
 
 
  template <class Solver, class Traits>
  struct DirectSolverCreator
      : public ISTLSolverCreatorBase<Traits>
  {
    using Super = ISTLSolverCreatorBase<Traits>;
 
    void init(std::string const& prefix) override
    {
      Super::init(prefix);
      Parameters::get(prefix + "->reuse vector", reuseVector_);
    }
 
    std::unique_ptr<typename Traits::Solver>
    createSolver(typename Traits::M const& mat, typename Traits::Comm const& comm) const override
    {
      test_exit(Dune::SolverCategory::category(comm) == Dune::SolverCategory::sequential,
        "Direct solver can be used as sequential solver only.");
      return std::make_unique<Solver>(mat, this->info_, reuseVector_);
    }
 
  protected:
    bool reuseVector_ = true;
  };
 
 
  template <class Solver, class Traits>
  using ISTLSolverCreator = std::conditional_t<Dune::IsDirectSolver<Solver>::value,
    DirectSolverCreator<Solver,Traits>,
    IterativeSolverCreator<Solver,Traits>>;
 
} // end namespace AMDiS