Creators.hpp

#pragma once
 
#include <memory>
#include <utility>
 
#include <dune/common/typeutilities.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/scalarproducts.hh>
#include <dune/istl/solvercategory.hh>
 
#if HAVE_MPI
#include <dune/istl/novlpschwarz.hh>
#include <dune/istl/schwarz.hh>
#endif
 
#include <amdis/Output.hpp>
#include <amdis/linearalgebra/istl/PreconWrapper.hpp>
 
namespace AMDiS
{
  template <class X>
  class ISTLScalarProductCreator
  {
  public:
    using Interface = Dune::ScalarProduct<X>;
 
    template <class Comm>
    static std::unique_ptr<Interface> create(Dune::SolverCategory::Category cat, Comm const& comm)
    {
      return create_impl(cat, comm.impl(), Dune::PriorityTag<10>{});
    }
 
  private:
    static std::unique_ptr<Interface>
    create_impl(Dune::SolverCategory::Category cat, Dune::Amg::SequentialInformation, Dune::PriorityTag<2>)
    {
      assert(cat == Dune::SolverCategory::sequential);
      return std::make_unique< Dune::SeqScalarProduct<X> >();
    }
 
    template <class Comm>
    static std::unique_ptr<Interface>
    create_impl(Dune::SolverCategory::Category cat, Comm const& comm, Dune::PriorityTag<1>)
    {
      switch (cat)
      {
        case Dune::SolverCategory::sequential:
          return std::make_unique< Dune::SeqScalarProduct<X> >();
#if HAVE_MPI
        case Dune::SolverCategory::overlapping:
          return std::make_unique< Dune::OverlappingSchwarzScalarProduct<X,Comm> >(comm);
        case Dune::SolverCategory::nonoverlapping:
          return std::make_unique< Dune::NonoverlappingSchwarzScalarProduct<X,Comm> >(comm);
#endif
        default:
          error_exit("Invalid solver category for scalar product\n");
          return nullptr; // avoid warnings
      }
    }
  };
 
  template <class M, class X, class Y>
  class ISTLLinearOperatorCreator
  {
  public:
    using Interface = Dune::AssembledLinearOperator<M,X,Y>;
 
    template <class Comm>
    static std::unique_ptr<Interface> create(Dune::SolverCategory::Category cat, M const& mat, Comm const& comm)
    {
      return create_impl(cat, mat, comm.impl(), Dune::PriorityTag<10>{});
    }
 
  private:
    static std::unique_ptr<Interface>
    create_impl(Dune::SolverCategory::Category cat, M const& mat, Dune::Amg::SequentialInformation, Dune::PriorityTag<2>)
    {
      assert(cat == Dune::SolverCategory::sequential);
      return std::make_unique< Dune::MatrixAdapter<M,X,Y> >(mat);
    }
 
    template <class Comm>
    static std::unique_ptr<Interface>
    create_impl(Dune::SolverCategory::Category cat, M const& mat, Comm const& comm, Dune::PriorityTag<1>)
    {
      switch (cat)
      {
        case Dune::SolverCategory::sequential:
          return std::make_unique< Dune::MatrixAdapter<M,X,Y> >(mat);
#if HAVE_MPI
        case Dune::SolverCategory::overlapping:
          return std::make_unique< Dune::OverlappingSchwarzOperator<M,X,Y,Comm> >(mat, comm);
        case Dune::SolverCategory::nonoverlapping:
          return std::make_unique< Dune::NonoverlappingSchwarzOperator<M,X,Y,Comm> >(mat, comm);
#endif
        default:
          error_exit("Invalid solver category for linear operator\n");
          return nullptr; // avoid warnings
      }
    }
  };
 
 
  template <class X, class Y>
  class ISTLParallelPreconditionerCreator
  {
  public:
    using Interface = Dune::Preconditioner<X,Y>;
 
    template <class Comm>
    static std::unique_ptr<Interface> create(Dune::SolverCategory::Category cat, std::unique_ptr<Interface> prec, Comm const& comm)
    {
      return create_impl(cat, std::move(prec), comm.impl(), Dune::PriorityTag<10>{});
    }
 
  private:
    static std::unique_ptr<Interface>
    create_impl(Dune::SolverCategory::Category cat, std::unique_ptr<Interface> prec, Dune::Amg::SequentialInformation, Dune::PriorityTag<2>)
    {
      assert(cat == Dune::SolverCategory::sequential);
      return prec;
    }
 
    template <class Comm>
    static std::unique_ptr<Interface>
    create_impl(Dune::SolverCategory::Category cat, std::unique_ptr<Interface> prec, Comm const& comm, Dune::PriorityTag<1>)
    {
      switch (cat)
      {
        case Dune::SolverCategory::sequential:
          return prec;
#if HAVE_MPI
        case Dune::SolverCategory::overlapping:
        {
          using BP = Dune::BlockPreconditioner<X,Y,Comm>;
          return std::make_unique< PreconWrapper<BP,Interface> >(std::move(prec), comm);
        }
        case Dune::SolverCategory::nonoverlapping:
        {
          using BP = Dune::NonoverlappingBlockPreconditioner<Comm,Interface>;
          return std::make_unique< PreconWrapper<BP,Interface> >(std::move(prec), comm);
        }
#endif
        default:
          error_exit("Invalid solver category for parallel preconditioner\n");
          return nullptr; // avoid warnings
      }
    }
  };
 
} // end namespace AMDiS