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>
#include <amdis/linearalgebra/istl/precompiled/Solvers.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