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