Preconditioners.hpp

#pragma once

#include <type_traits>

#include <amdis/CreatorInterface.hpp>
#include <amdis/CreatorMap.hpp>
#include <amdis/common/TypeTraits.hpp>
#include <amdis/linearalgebra/istl/AMGPrecon.hpp>
#include <amdis/linearalgebra/istl/ISTLPreconCreator.hpp>
#include <amdis/linearalgebra/istl/precompiled/Preconditioners.hpp>

namespace AMDiS
{

  template <class Traits>
  class DefaultCreators<tag::preconditioner<Traits>>
  {
    using M = typename Traits::M;
    using X = typename Traits::X;
    using Y = typename Traits::Y;

    template <class Precon>
    using Creator = ISTLPreconCreator<Precon, Traits>;

    using Map = CreatorMap<tag::preconditioner<Traits>>;
    using FTraits = Dune::FieldTraits<typename M::field_type>;

  public:
    static void init()
    {
      auto jacobi = new Creator<Dune::SeqJac<M,X,Y>>;
      Map::addCreator("diag", jacobi);
      Map::addCreator("jacobi", jacobi);

      auto gs = new Creator<Dune::SeqGS<M,X,Y>>;
      Map::addCreator("gs", gs);
      Map::addCreator("gauss_seidel", gs);

      auto sor = new Creator<Dune::SeqSOR<M,X,Y>>;
      Map::addCreator("sor", sor);

      auto ssor = new Creator<Dune::SeqSSOR<M,X,Y>>;
      Map::addCreator("ssor", ssor);

      init_ilu(std::is_arithmetic<typename FTraits::field_type>{});
      init_amg(std::is_same<typename FTraits::real_type, double>{});

      auto richardson = new Creator<Dune::Richardson<X,Y>>;
      Map::addCreator("richardson", richardson);
      Map::addCreator("default", richardson);

      auto solver = new Creator<tag::solver<Traits>>;
      Map::addCreator("solver", solver);

      init_bjacobi(Types<TYPEOF(std::declval<typename Traits::Comm>().impl())>{}, Dune::PriorityTag<10>{});
    }

    static void init_ilu(std::false_type)
    {
      warning("ILU preconditioners not created for the matrix with field_type = {}.",
        Dune::className<typename FTraits::field_type>());
    }

    static void init_ilu(std::true_type)
    {
      auto ilu = new Creator<Dune::SeqILU<M,X,Y>>;
      Map::addCreator("ilu", ilu);
      Map::addCreator("ilu0", ilu);

      auto ildl = new Creator<Dune::SeqILDL<M,X,Y>>;
      Map::addCreator("ildl", ildl);
    }

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

    static void init_amg(std::true_type)
    {
      auto amg = new AMGPreconCreator<Dune::Amg::AMG, Traits>;
      Map::addCreator("amg", amg);
      auto fastamg = new AMGPreconCreator<Dune::Amg::FastAMG, Traits>;
      Map::addCreator("fastamg", fastamg);
      auto kamg = new AMGPreconCreator<Dune::Amg::KAMG, Traits>;
      Map::addCreator("kamg", kamg);
    }

    static void init_bjacobi(Types<Dune::Amg::SequentialInformation>, Dune::PriorityTag<2>) {}

    template <class Comm>
    static void init_bjacobi(Types<Comm>, Dune::PriorityTag<1>)
    {
      auto pssor = new Creator<Dune::ParSSOR<M,X,Y,Comm>>;
      Map::addCreator("pssor", pssor);

      auto bjacobi = new Creator<tag::bjacobi>;
      Map::addCreator("bjacobi", bjacobi);
    }
  };

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

} // end namespace AMDiS