SolverCreator.hpp

#pragma once

#include <memory>

#include <dune/common/classname.hh>
#include <dune/common/version.hh>

#include <dune/common/ftraits.hh>

#include <amdis/CreatorMap.hpp>
#include <amdis/Environment.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/Output.hpp>

#include <amdis/linearalgebra/LinearSolver.hpp>
#include <amdis/linearalgebra/istl/CreatorInterfaces.hpp>
#include <amdis/linearalgebra/istl/ISTLSolver.hpp>
#include <amdis/linearalgebra/istl/SolverWrapper.hpp>
#include <amdis/linearalgebra/istl/precompiled/Solvers.hpp>

namespace AMDiS
{
  namespace tag
  {
    template <class Solver> struct pcg {};
    template <class Solver> struct gmres {};
  }


  template <class Model, class Traits>
  struct ISTLSolverCreator
      : public ISTLSolverCreatorInterface<Traits>
  {
    using Interface = ISTLSolverCreatorInterface<Traits>;
    struct Creator : CreatorInterfaceName<Interface>
    {
      std::unique_ptr<Interface> createWithString(std::string prefix) override
      {
        return std::make_unique<Model>(prefix);
      }
    };

    explicit ISTLSolverCreator(std::string const& prefix)
    {
      if (Environment::mpiRank() == 0)
        Parameters::get(prefix + "->info", info_);
    }

  protected:
    int info_ = 0;
  };



  template <class Creator, class Traits>
  class ISTLIterativeSolverCreator
      : public ISTLSolverCreator<Creator, Traits>
  {
    using Super = ISTLSolverCreator<Creator, Traits>;

    using LinOp = typename Traits::LinOp;
    using ScalProd = typename Traits::ScalProd;
    using Prec = typename Traits::Prec;

    using real_type = typename Dune::FieldTraits<typename Traits::M::field_type>::real_type;

  public:
    explicit ISTLIterativeSolverCreator(std::string const& prefix)
      : Super(prefix)
    {
      Parameters::get(prefix + "->max iteration", maxIter_);
      Parameters::get(prefix + "->relative tolerance", rTol_);

      std::string precon = "default";
      Parameters::get(prefix + "->precon", precon);

      using CreatorMap = CreatorMap<ISTLPreconCreatorInterface<Traits>>;
      auto* creator = named(CreatorMap::getCreator(precon, prefix + "->precon"));
      preconCreator_ = creator->createWithString(prefix + "->precon");
      assert(preconCreator_);
    }

  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);

      auto precon = preconCreator_->create(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;
    std::shared_ptr<ISTLPreconCreatorInterface<Traits>> preconCreator_;
  };



  template <class Solver, class Traits>
  class IterativeSolverCreator
      : public ISTLIterativeSolverCreator<IterativeSolverCreator<Solver,Traits>, Traits>
  {
    using Super = ISTLIterativeSolverCreator<IterativeSolverCreator,Traits>;
    using Interface = typename Traits::Solver;

  public:
    using Super::Super;

    std::unique_ptr<Interface> create(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 IterativeSolverCreator<tag::gmres<Solver>, Traits>
      : public ISTLIterativeSolverCreator<IterativeSolverCreator<tag::gmres<Solver>,Traits>, Traits>
  {
    using Super = ISTLIterativeSolverCreator<IterativeSolverCreator, Traits>;
    using Interface = typename Traits::Solver;

  public:
    explicit IterativeSolverCreator(std::string const& prefix)
      : Super(prefix)
    {
      Parameters::get(prefix + "->restart", restart_);
    }

    std::unique_ptr<Interface> create(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 Solver, class Traits>
  struct IterativeSolverCreator<tag::pcg<Solver>, Traits>
      : public ISTLIterativeSolverCreator<IterativeSolverCreator<tag::pcg<Solver>,Traits>, Traits>
  {
    using Super = ISTLIterativeSolverCreator<IterativeSolverCreator, Traits>;
    using Interface = typename Traits::Solver;

  public:
    explicit IterativeSolverCreator(std::string const& prefix)
      : Super(prefix)
    {
      Parameters::get(prefix + "->restart", restart_);
    }

    std::unique_ptr<Interface> create(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 Solver, class Traits>
  struct DirectSolverCreator
      : public ISTLSolverCreator<DirectSolverCreator<Solver,Traits>, Traits>
  {
    using Super = ISTLSolverCreator<DirectSolverCreator,Traits>;

    explicit DirectSolverCreator(std::string const& prefix)
      : Super(prefix)
    {
      Parameters::get(prefix + "->reuse vector", reuseVector_);
    }

    std::unique_ptr<typename Traits::Solver>
    create(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;
  };

#if HAVE_SUITESPARSE_CHOLMOD && DUNE_VERSION_GTE(DUNE_ISTL,2,7)
  template <class X, class Traits>
  struct DirectSolverCreator<Dune::Cholmod<X>,Traits>
      : public ISTLSolverCreator<DirectSolverCreator<Dune::Cholmod<X>,Traits>, Traits>
  {
    using Super = ISTLSolverCreator<DirectSolverCreator,Traits>;
    using Super::Super;

    std::unique_ptr<typename Traits::Solver>
    create(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.");
      auto solver = std::make_unique<Dune::Cholmod<X>>();
      solver->setMatrix(mat);
      return std::move(solver);
    }
  };
#endif



  template <class Traits, class Interface>
  class DefaultSolverCreators
  {
    using M = typename Traits::M;
    using X = typename Traits::X;
    using Y = typename Traits::Y;

    using FTraits = Dune::FieldTraits<typename M::field_type>;

    template <class Solver>
    using IterativeSolver = typename IterativeSolverCreator<Solver,Traits>::Creator;

    template <class Solver>
    using DirectSolver = typename DirectSolverCreator<Solver,Traits>::Creator;

    using Map = CreatorMap<Interface>;

  public:
    static void init()
    {
      auto cg = new IterativeSolver<Dune::CGSolver<X>>;
      addCreator("cg", cg);

      // Generalized preconditioned conjugate gradient solver.
      auto pcg = new IterativeSolver<tag::pcg<Dune::GeneralizedPCGSolver<X>>>;
      addCreator("pcg", pcg);

#if DUNE_VERSION_GTE(DUNE_ISTL,2,7)
      auto fcg = new IterativeSolver<tag::pcg<Dune::RestartedFCGSolver<X>>>;
      addCreator("fcg", fcg);

      auto cfcg = new IterativeSolver<tag::pcg<Dune::CompleteFCGSolver<X>>>;
      addCreator("cfcg", cfcg);
#endif

      auto bicgstab = new IterativeSolver<Dune::BiCGSTABSolver<X>>;
      addCreator("bicgstab", bicgstab);
      addCreator("bcgs", bicgstab);
      addCreator("default", bicgstab);

      auto minres = new IterativeSolver<Dune::MINRESSolver<X>>;
      addCreator("minres", minres);

      auto gmres = new IterativeSolver<tag::gmres<Dune::RestartedGMResSolver<X,Y>>>;
      addCreator("gmres", gmres);

#if DUNE_VERSION_GTE(DUNE_ISTL,2,7)
      auto fgmres = new IterativeSolver<tag::gmres<Dune::RestartedFlexibleGMResSolver<X,Y>>>;
      addCreator("fgmres", fgmres);
#endif

      init_direct(std::is_same<typename FTraits::real_type, double>{});
    }

    static void init_direct(std::false_type)
    {
      warning("Direct solvers not created for the matrix with real_type = {}.",
        Dune::className<typename FTraits::real_type>());
    }

    static void init_direct(std::true_type)
    {
#if HAVE_SUITESPARSE_UMFPACK
      auto umfpack = new DirectSolver<Dune::UMFPack<M>>;
      addCreator("umfpack", umfpack);
#endif

#if HAVE_SUITESPARSE_LDL
      auto ldl = new DirectSolver<Dune::LDL<M>>;
      addCreator("ldl", ldl);
#endif

#if HAVE_SUITESPARSE_SPQR
      auto spqr = new DirectSolver<Dune::SPQR<M>>;
      addCreator("spqr", spqr);
#endif

#if HAVE_SUITESPARSE_CHOLMOD && DUNE_VERSION_GTE(DUNE_ISTL,2,7)
      auto cholmod = new DirectSolver<Dune::Cholmod<X>>;
      addCreator("cholmod", cholmod);
#endif

#if HAVE_SUPERLU
      auto superlu = new DirectSolver<Dune::SuperLU<M>>;
      addCreator("superlu", superlu);
#endif

      // default direct solvers
#if HAVE_SUITESPARSE_UMFPACK
      addCreator("direct", umfpack);
#elif HAVE_SUPERLU
      addCreator("direct", superlu);
#endif
    }

  private:
    template <class T> struct Type {};

    template <class Creator>
    static void addCreator(std::string name, Creator* creator)
    {
      addCreatorImpl(name, creator, Type<Interface>{});
    }

    template <class Creator>
    static void addCreatorImpl(std::string name, Creator* creator, Type<ISTLSolverCreatorInterface<Traits>>)
    {
      Map::addCreator(name, creator);
    }

    template <class Creator, class Mat, class Vec>
    static void addCreatorImpl(std::string name, Creator* creator, Type<LinearSolverInterface<Mat,Vec>>)
    {
      using LinearSolverCreator = typename ISTLSolver<Mat, Vec>::Creator;
      Map::addCreator(name, new LinearSolverCreator);
    }
  };


  template <class Traits>
  class DefaultCreators< ISTLSolverCreatorInterface<Traits> >
      : public DefaultSolverCreators<Traits, ISTLSolverCreatorInterface<Traits>>
  {};

  template <class Mat, class Vec>
  class DefaultCreators< LinearSolverInterface<Mat,Vec> >
      : public DefaultSolverCreators<SolverTraits<Mat, Vec>, LinearSolverInterface<Mat,Vec>>
  {};


  // extern template declarations:
  extern template class DefaultCreators< ISTLSolverCreatorInterface<Precompiled::SolverTraits> >;

} // end namespace AMDiS