ISTL_Preconditioner.hpp

#pragma once

#include <dune/istl/preconditioners.hh>

#include <amdis/CreatorInterface.hpp>
#include <amdis/linearalgebra/istl/AMGPrecon.hpp>
#include <amdis/linearalgebra/istl/Fwd.hpp>
#include <amdis/linearalgebra/istl/ISTLPreconCreatorInterface.hpp>

namespace AMDiS
{
  template <class> struct Type {};

  // creators for preconditioners, depending on matrix-type
  // ---------------------------------------------------------------------------

  template <class Precon, class Matrix>
  struct ISTLPreconCreator
      : public ISTLPreconCreatorInterface<Matrix, typename Precon::domain_type, typename Precon::range_type>
  {
    using Sol = typename Precon::domain_type;
    using Rhs = typename Precon::range_type;
    using Super = ISTLPreconCreatorInterface<Matrix, Sol, Rhs>;
    using Self = ISTLPreconCreator;

    struct Creator : CreatorInterfaceName<Super>
    {
      std::unique_ptr<Super> createWithString(std::string prefix) override
      {
        return std::make_unique<Self>(prefix);
      }
    };

    ISTLPreconCreator(std::string const& prefix)
    {
      Parameters::get(prefix + "->relaxation", w_);
      Parameters::get(prefix + "->iterations", iter_);
    }

    using PreconBase = Dune::Preconditioner<Sol, Rhs>;
    std::unique_ptr<PreconBase> create(Matrix const& A) const override
    {
      return createImpl(A, Type<Precon>{});
    }

  private:
    template <class P>
    std::unique_ptr<P> createImpl(Matrix const& A, Type<P>) const
    {
      return std::make_unique<P>(A, iter_, w_);
    }

    std::unique_ptr<Dune::SeqILU<Matrix, Sol, Rhs>>
    createImpl(Matrix const& A, Type<Dune::SeqILU<Matrix, Sol, Rhs>>) const
    {
      return std::make_unique<Dune::SeqILU<Matrix, Sol, Rhs>>(A, iter_, w_);
    }

    std::unique_ptr<Dune::Richardson<Sol, Rhs>>
    createImpl(Matrix const& /*A*/, Type<Dune::Richardson<Sol, Rhs>>) const
    {
      return std::make_unique<Dune::Richardson<Sol, Rhs>>(w_);
    }

    double w_ = 1.0;
    int iter_ = 0;
  };


  /// Adds default creators for linear solvers based on `Dune::BCRSMatrix`.
  /**
   * Adds creators for full-matrix aware solvers.
   * - *cg*: conjugate gradient method, \see Dune::CGSolver
   * - *bcgs*: stabilized bi-conjugate gradient method, \see Dune::BiCGSTABSolver
   * - *minres*: Minimal residul method, \see Dune::MINRESSolver
   * - *gmres*: Generalized minimal residula method, \see Dune::RestartedGMResSolver
   * - *umfpack*: external UMFPACK solver, \see Dune::UMFPack
   **/
  template <class Mat, class Sol, class Rhs>
  class DefaultCreators<ISTLPreconCreatorInterface<Mat, Sol, Rhs>>
  {
    template <class Precon>
    using PreconCreator
      = typename ISTLPreconCreator<Precon, Mat>::Creator;

    using Map = CreatorMap<ISTLPreconCreatorInterface<Mat, Sol, Rhs>>;

  public:
    static void init()
    {
      auto jacobi = new PreconCreator<Dune::SeqJac<Mat, Sol, Rhs>>;
      Map::addCreator("diag", jacobi);
      Map::addCreator("jacobi", jacobi);

      auto gs = new PreconCreator<Dune::SeqGS<Mat, Sol, Rhs>>;
      Map::addCreator("gs", gs);
      Map::addCreator("gauss_seidel", gs);

      auto sor = new PreconCreator<Dune::SeqSOR<Mat, Sol, Rhs>>;
      Map::addCreator("sor", sor);

      auto ssor = new PreconCreator<Dune::SeqSSOR<Mat, Sol, Rhs>>;
      Map::addCreator("ssor", ssor);

      auto ilu = new PreconCreator<Dune::SeqILU<Mat, Sol, Rhs>>;
      Map::addCreator("ilu", ilu);
      Map::addCreator("ilu0", ilu);

      auto amg = new AMGPreconCreator<Dune::Amg::AMG, Mat, Sol, Rhs>;
      Map::addCreator("amg", amg);
      auto fastamg = new AMGPreconCreator<Dune::Amg::FastAMG, Mat, Sol, Rhs>;
      Map::addCreator("fastamg", fastamg);

      auto richardson = new PreconCreator<Dune::Richardson<Sol, Rhs>>;
      Map::addCreator("richardson", richardson);
      Map::addCreator("default", richardson);
    }
  };

} // end namespace AMDiS