ISTLPreconCreator.hpp

#pragma once

#include <memory>
#include <string>

#include <dune/istl/preconditioners.hh>
#include <dune/istl/novlpschwarz.hh>
#include <dune/istl/schwarz.hh>

#include <amdis/CreatorInterface.hpp>
#include <amdis/CreatorMap.hpp>
#include <amdis/linearalgebra/istl/PreconWrapper.hpp>
#include <amdis/linearalgebra/istl/Traits.hpp>

namespace AMDiS
{
  namespace tag
  {
    struct bjacobi {};

    // forward declaration
    template <class Traits>
    struct solver;

    template <class Traits>
    struct preconditioner {};
  }



  template <class Traits>
  class ISTLPreconCreatorBase
      : public CreatorInterface<tag::preconditioner<Traits>>
  {
    using X = typename Traits::X;
    using Y = typename Traits::Y;

    std::unique_ptr<tag::preconditioner<Traits>> create() final { return {}; };

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

    virtual std::unique_ptr<Dune::Preconditioner<X,Y>>
    createPrecon(typename Traits::M const& A, typename Traits::Comm const& comm) const = 0;

  protected:
    double w_ = 1.0;
    int iter_ = 1;
  };



  template <class Precon, class Traits>
  struct ISTLPreconCreator
      : public ISTLPreconCreatorBase<Traits>
  {
    std::unique_ptr<typename Traits::Prec>
    createPrecon(typename Traits::M const& mat, typename Traits::Comm const& /*comm*/) const override
    {
      return std::make_unique<Precon>(mat, this->iter_, this->w_);
    }
  };


  template <class X, class Y, class Traits>
  struct ISTLPreconCreator<Dune::Richardson<X, Y>, Traits>
      : public ISTLPreconCreatorBase<Traits>
  {
    std::unique_ptr<typename Traits::Prec>
    createPrecon(typename Traits::M const& /*mat*/, typename Traits::Comm const& /*comm*/) const override
    {
      using Precon = Dune::Richardson<X, Y>;
      return std::make_unique<Precon>(this->w_);
    }
  };


  template <class M, class X, class Y, class Traits>
  struct ISTLPreconCreator<Dune::SeqILDL<M, X, Y>, Traits>
      : public ISTLPreconCreatorBase<Traits>
  {
    std::unique_ptr<typename Traits::Prec>
    createPrecon(typename Traits::M const& mat, typename Traits::Comm const& /*comm*/) const override
    {
      using Precon = Dune::SeqILDL<M, X, Y>;
      return std::make_unique<Precon>(mat, this->w_);
    }
  };



  template <class M, class X, class Y, class Comm, class Traits>
  struct ISTLPreconCreator<Dune::ParSSOR<M, X, Y, Comm>, Traits>
      : public ISTLPreconCreatorBase<Traits>
  {
    std::unique_ptr<typename Traits::Prec>
    createPrecon(typename Traits::M const& mat, typename Traits::Comm const& comm) const override
    {
      test_exit(Dune::SolverCategory::category(comm) == Dune::SolverCategory::overlapping,
        "Dune::ParSSOR preconditioner can be used with overlapping domain decomposition.");

      using Precon = Dune::ParSSOR<M,X,Y,Comm>;
      return std::make_unique<Precon>(mat, this->iter_, this->w_, comm.impl());
    }
  };


  // forward declaration
  template <class Traits>
  class ISTLSolverCreatorBase;



  template <class Traits>
  struct ISTLPreconCreator<tag::solver<Traits>, Traits>
      : public ISTLPreconCreatorBase<Traits>
  {
    using Super = ISTLPreconCreatorBase<Traits>;

    void init(std::string const& prefix) override
    {
      Super::init(prefix);

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

      using CreatorMap = CreatorMap<tag::solver<Traits>>;
      auto* creator = CreatorMap::get(solver, prefix + "->solver");
      solverCreator_ = dynamic_cast<ISTLSolverCreatorBase<Traits>*>(creator);
      assert(solverCreator_ != nullptr);

      solverCreator_->init(prefix + "->solver");
    }

    std::unique_ptr<typename Traits::Prec>
    createPrecon(typename Traits::M const& mat, typename Traits::Comm const& comm) const override
    {
      using InverseOp = Dune::InverseOperator<typename Traits::X, typename Traits::Y>;
      using Precon = Dune::InverseOperator2Preconditioner<InverseOp>;
      using Wrapper = PreconWrapper<Precon, InverseOp>;

      assert(solverCreator_ != nullptr);
      return std::make_unique<Wrapper>(solverCreator_->createSolver(mat, comm));
    }

  private:
    ISTLSolverCreatorBase<Traits>* solverCreator_ = nullptr;
  };



  template <class Traits>
  struct ISTLPreconCreator<tag::bjacobi, Traits>
      : public ISTLPreconCreatorBase<Traits>
  {
    using Super = ISTLPreconCreatorBase<Traits>;
    using SeqTraits = SeqSolverTraits<Traits>;

    void init(std::string const& prefix) override
    {
      Super::init(prefix);

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

      using CreatorMap = CreatorMap<tag::preconditioner<SeqTraits>>;
      auto* creator = CreatorMap::get(subPrecon, prefix + "->sub precon");
      subPreconCreator_ = dynamic_cast<ISTLPreconCreatorBase<SeqTraits>*>(creator);
      assert(subPreconCreator_ != nullptr);
      subPreconCreator_->init(prefix + "->sub precon");
    }

    std::unique_ptr<typename Traits::Prec>
    createPrecon(typename Traits::M const& mat, typename Traits::Comm const& comm) const override
    {
      assert(subPreconCreator_ != nullptr);
      return Traits::ParPrecCreator::create(Dune::SolverCategory::category(comm),
        subPreconCreator_->createPrecon(mat, comm.sequential()),
        comm);
    }

  private:
    ISTLPreconCreatorBase<SeqTraits>* subPreconCreator_ = nullptr;
  };

} // end namespace AMDiS