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