Preconditioners.hpp

#pragma once
 
// MTL4 headers
#include <boost/numeric/itl/itl.hpp>
#include <boost/numeric/itl/pc/ilu_0.hpp>
#include <boost/numeric/itl/pc/ic_0.hpp>
#include <boost/numeric/mtl/vector/assigner.hpp>
 
#include <amdis/CreatorMap.hpp>
#include <amdis/linearalgebra/mtl/HyprePrecon.hpp>
#include <amdis/linearalgebra/mtl/Preconditioner.hpp>
#include <amdis/linearalgebra/mtl/SolverPrecon.hpp>
#include <amdis/linearalgebra/mtl/Traits.hpp>
#include <amdis/linearalgebra/mtl/itl/masslumping.hpp>
 
namespace AMDiS
{
  template <class Matrix>
  using DiagonalPreconditioner = itl::pc::diagonal<Matrix>;
 
  template <class Matrix>
  using MassLumpingPreconditioner = itl::pc::masslumping<Matrix>;
 
 
  template <class Matrix>
  using IdentityPreconditioner = itl::pc::identity<Matrix>;
 
 
  template <class Matrix>
  using ILUPreconditioner = itl::pc::ilu_0<Matrix>;
 
 
  template <class Matrix>
  using ICPreconditioner = itl::pc::ic_0<Matrix>;
 
 
  template <class M, class X, class Y>
  class DefaultCreators< PreconditionerInterface<M,X,Y> >
  {
    using PreconBase = PreconditionerInterface<M,X,Y>;
 
    template <template <class> class ITLPrecon>
    using PreconCreator = typename Preconditioner<M,X,Y,ITLPrecon>::Creator;
 
    using Map = CreatorMap<PreconBase>;
 
  public:
    static void init()
    {
      auto pc_diag = new PreconCreator<DiagonalPreconditioner>;
      Map::addCreator("diag", pc_diag);
      Map::addCreator("jacobi", pc_diag);
 
      auto pc_mass = new PreconCreator<MassLumpingPreconditioner>;
      Map::addCreator("masslumping", pc_mass);
 
      auto pc_ilu = new PreconCreator<ILUPreconditioner>;
      Map::addCreator("ilu", pc_ilu);
      Map::addCreator("ilu0", pc_ilu);
 
      auto pc_ic = new PreconCreator<ICPreconditioner>;
      Map::addCreator("ic", pc_ic);
      Map::addCreator("ic0", pc_ic);
 
      auto pc_id = new PreconCreator<IdentityPreconditioner>;
      Map::addCreator("identity", pc_id);
      Map::addCreator("no", pc_id);
 
      auto pc_solver = new typename SolverPrecon<M,X,Y>::Creator;
      Map::addCreator("solver", pc_solver);
 
      Map::addCreator("default", pc_id);
 
      init_hypre(std::is_same<typename Dune::FieldTraits<typename M::value_type>::real_type, double>{});
    }
 
    static void init_hypre(std::false_type) {}
    static void init_hypre(std::true_type)
    {
#if AMDIS_HAS_HYPRE && HAVE_MPI
      auto pc_hypre = new typename HyprePrecon<M,X,Y>::Creator;
      Map::addCreator("hypre", pc_hypre);
#endif
    }
  };
 
  template <class M, class X, class Y>
  itl::pc::solver<PreconditionerInterface<M,X,Y>, X, false>
  solve(PreconditionerInterface<M,X,Y> const& P, X const& vin)
  {
    return {P, vin};
  }
 
  template <class M, class X, class Y>
  itl::pc::solver<PreconditionerInterface<M,X,Y>, X, true>
  adjoint_solve(PreconditionerInterface<M,X,Y> const& P, X const& vin)
  {
    return {P, vin};
  }
 
} // namespace AMDiS