AMGPrecon.hpp

#pragma once

#include <memory>

#include <dune/common/unused.hh>
#include <dune/common/version.hh>

#include <dune/istl/novlpschwarz.hh>
#include <dune/istl/schwarz.hh>
#include <dune/istl/paamg/amg.hh>
#include <dune/istl/paamg/fastamg.hh>
#include <dune/istl/paamg/kamg.hh>

#include <amdis/Initfile.hpp>
#include <amdis/common/ConceptsBase.hpp>
#include <amdis/linearalgebra/istl/IndexDistribution.hpp>
#include <amdis/linearalgebra/istl/ISTLPreconCreator.hpp>
#include <amdis/linearalgebra/istl/precompiled/Preconditioners.hpp>
#include <amdis/linearalgebra/istl/precompiled/Solvers.hpp>

namespace AMDiS
{
  template <template <class...> class AMGSolver, class Traits>
  struct AMGCreator;

  template <class Smoother>
  using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;


  // Specialization for \ref Dune::Amg::AMG preconditioner
  template <class Traits>
  struct AMGCreator<Dune::Amg::AMG, Traits>
  {
    using PrecBase = typename Traits::Prec;

    template <class Smoother, class LinOp, class Criterion, class Comm>
    static std::unique_ptr<PrecBase>
    create([[maybe_unused]] std::string prefix, LinOp const& linOp, Criterion const& criterion, SmootherArgs<Smoother> const& smootherArgs, Comm const& comm)
    {
      using Solver = Dune::Amg::AMG<LinOp, typename Traits::X, Smoother, Comm>;
      return std::make_unique<Solver>(linOp, criterion, smootherArgs, comm);
    }
  };


  // Specialization for \ref Dune::Amg::FastAMG preconditioner
  template <class Traits>
  struct AMGCreator<Dune::Amg::FastAMG, Traits>
  {
    using PrecBase = typename Traits::Prec;

    template <class Smoother, class LinOp, class Criterion, class SmootherArgs, class Comm>
    static std::unique_ptr<PrecBase>
    create(std::string const& prefix, LinOp const& linOp, Criterion const& criterion,
           SmootherArgs const& smootherArgs, Comm const& comm)
    {
      return createImpl1<Smoother>(prefix, linOp, criterion, smootherArgs, comm, Dune::PriorityTag<5>{});
    }

  private: // step 1:

    template <class Smoother, class LinOp, class Criterion>
    static std::unique_ptr<PrecBase>
    createImpl1(std::string prefix, LinOp const& linOp, Criterion const& criterion,
                [[maybe_unused]] SmootherArgs<Smoother> const& smootherArgs,
                Dune::Amg::SequentialInformation const& comm, Dune::PriorityTag<2>)
    {
      bool symmetric = Parameters::get<bool>(prefix + "->symmetric").value_or(true);

      using Solver = Dune::Amg::FastAMG<LinOp, typename Traits::X, Dune::Amg::SequentialInformation>;
      return std::make_unique<Solver>(linOp, criterion, criterion, symmetric, comm);
    }

    template <class Smoother, class LinOp, class Criterion, class SmootherArgs, class Comm>
    static std::unique_ptr<PrecBase>
    createImpl1(std::string, LinOp const&, Criterion const&, SmootherArgs const&, Comm const&,
                Dune::PriorityTag<1>)
    {
      error_exit("Currently only sequential FastAMG supported.");
      return nullptr;
    }
  };


  // Specialization for \ref Dune::Amg::KAMG preconditioner
  template <class Traits>
  struct AMGCreator<Dune::Amg::KAMG, Traits>
  {
    using X = typename Traits::X;
    using PrecBase = typename Traits::Prec;

    template <class Smoother, class LinOp, class Criterion, class Comm>
    static std::unique_ptr<PrecBase>
    create(std::string prefix, LinOp const& linOp, Criterion const& criterion,
           [[maybe_unused]] SmootherArgs<Smoother> const& smootherArgs, Comm const& comm)
    {
      std::string solver = Parameters::get<std::string>(prefix + "->krylov solver").value_or("default");

      std::size_t maxLevelKrylovSteps = 3;
      Parameters::get(prefix + "->krylov solver->maxLevelKrylovSteps", maxLevelKrylovSteps);
      double minDefectReduction = 1.e-1;
      Parameters::get(prefix + "->krylov solver->minDefectReduction", minDefectReduction);

      if (solver == "pcg" || solver == "default")
      {
        using Solver = Dune::Amg::KAMG<LinOp, X, Smoother, Comm, Dune::GeneralizedPCGSolver<X>>;
        return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
      }
      else if (solver == "fcg")
      {
        using Solver = Dune::Amg::KAMG<LinOp, X, Smoother, Comm, Dune::RestartedFCGSolver<X>>;
        return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
      }
      else if (solver == "cfcg")
      {
        using Solver = Dune::Amg::KAMG<LinOp, X, Smoother, Comm, Dune::CompleteFCGSolver<X>>;
        return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
      }
      else if (solver == "bicgstab" || solver == "bcgs")
      {
        using Solver = Dune::Amg::KAMG<LinOp, X, Smoother, Comm, Dune::BiCGSTABSolver<X>>;
        return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
      }
      else if (solver == "minres")
      {
        using Solver = Dune::Amg::KAMG<LinOp, X, Smoother, Comm, Dune::MINRESSolver<X>>;
        return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
      }
#if 0
      // NOTE: can not be constructed inside the KAMG precon, since additional constructor argument.
      // Needs something like ConstructionTraits for solvers.
      else if (solver == "gmres")
      {
        using Solver = Dune::Amg::KAMG<LinOp, X, Smoother, Comm, Dune::RestartedGMResSolver<X>>;
        return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
      }
#endif
      else
      {
        error_exit("Unknown coarse space solver {} given for parameter `{}`.", solver, prefix + "->coarse space solver");
        return nullptr;
      }
    }
  };



  template <template <class...> class AMGSolver, class Traits>
  class AMGPreconCreator
      : public ISTLPreconCreatorBase<Traits>
  {
    using M = typename Traits::M;
    using X = typename Traits::X;
    using Y = typename Traits::Y;

    using Interface = Dune::Preconditioner<X,Y>;
    using LinOperator = typename Traits::LinOpCreator::Interface;

    using SolverCategory = typename Dune::SolverCategory::Category;

    static const bool is_arithmetic = std::is_arithmetic_v<typename Dune::FieldTraits<M>::field_type>;
    using Norm = std::conditional_t<is_arithmetic, Dune::Amg::FirstDiagonal, Dune::Amg::RowSum>;

    using SymCriterion = Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<M,Norm>>;
    using UnSymCriterion = Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<M,Norm>>;

  public:

    void init(std::string const& prefix) override
    {
      prefix_ = prefix;
      initParameters(prefix);
    }

    std::unique_ptr<Interface>
    createPrecon(M const& mat, typename Traits::Comm const& comm) const override
    {
      return createImpl0(Dune::SolverCategory::category(comm), mat, comm.impl());
    }

  private: // step 0:

    template <class Comm>
    std::unique_ptr<Interface>
    createImpl0(SolverCategory cat, M const& mat, Comm const& comm) const
    {
      if (smoother_ == "diag" || smoother_ == "jacobi" || smoother_ == "default") {
        return createImpl1<Dune::SeqJac<M,X,Y>>(cat,mat,comm, Dune::PriorityTag<5>{});
      }
#if 0
      // NOTE: apply<forward> not implemented correctly in BlockPreconditioner and
      // NonoverlappingBlockPreconditioner. See !302 in dune-istl
      else if (smoother_ == "sor") {
        return createImpl1<Dune::SeqSOR<M,X,Y>>(cat,mat,comm, Dune::PriorityTag<5>{});
      }
#endif
#if 0
      // NOTE: ConstructionTraits not implemented for SeqGS. See !303 in dune-istl
      else if (smoother_ == "gs" || smoother_ == "gauss_seidel") {
        return createImpl1<Dune::SeqGS<M,X,Y>>(cat,mat,comm, Dune::PriorityTag<5>{});
      }
#endif
#if 0
      // NOTE: ConstructionTraits not implemented for Richardson. See !304 in dune-istl
      else if (smoother_ == "richardson") {
        return createImpl1<Dune::Richardson<X,Y>>(cat,mat,comm, Dune::PriorityTag<5>{});
      }
#endif
      else if (smoother_ == "ssor") {
        return createImpl1<Dune::SeqSSOR<M,X,Y>>(cat,mat,comm, Dune::PriorityTag<5>{});
      }
      else {
        error_exit("Unknown smoother type '{}' given for parameter '{}'", smoother_, prefix_ + "->smoother");
        return nullptr;
      }
    }

  private: // step 1:

    template <class Smoother>
    std::unique_ptr<Interface>
    createImpl1([[maybe_unused]] SolverCategory cat, M const& mat,
                Dune::Amg::SequentialInformation const& comm, Dune::PriorityTag<2>) const
    {
      assert(cat == SolverCategory::sequential);
      using LinOp = Dune::MatrixAdapter<M,X,Y>;
      LinOp* linOpPtr = new LinOp(mat);
      linOperator_.reset(linOpPtr);
      return createImpl2<Smoother>(*linOpPtr, comm);
    }

    template <class Smoother, class Comm>
    std::unique_ptr<Interface>
    createImpl1(SolverCategory cat, M const& mat, Comm const& comm, Dune::PriorityTag<1>) const
    {
      switch (cat) {
        case SolverCategory::sequential:
        {
          return createImpl1<Smoother>(cat, mat, Dune::Amg::SequentialInformation{}, Dune::PriorityTag<5>{});
        }
#if HAVE_MPI
        case SolverCategory::overlapping:
        {
          using LinOp = Dune::OverlappingSchwarzOperator<M,X,Y,Comm>;
          using ParSmoother = Dune::BlockPreconditioner<X,Y,Comm,Smoother>;
          LinOp* linOpPtr = new LinOp(mat, comm);
          linOperator_.reset(linOpPtr);
          return createImpl2<ParSmoother>(*linOpPtr, comm);
        }
        case SolverCategory::nonoverlapping:
        {
          using LinOp = Dune::NonoverlappingSchwarzOperator<M,X,Y,Comm>;
          using ParSmoother = Dune::NonoverlappingBlockPreconditioner<Comm,Smoother>;
          LinOp* linOpPtr = new LinOp(mat, comm);
          linOperator_.reset(linOpPtr);
          return createImpl2<ParSmoother>(*linOpPtr, comm);
        }
#endif
        default:
          error_exit("Invalid solver category for AMGSolver\n");
          return nullptr; // avoid warnings
      }
    }

  private: // step 2:

    template <class Smoother, class LinOp, class Comm>
    std::unique_ptr<Interface> createImpl2(LinOp const& linOp, Comm const& comm) const
    {
      SmootherArgs<Smoother> smootherArgs;
      initSmootherParameters(prefix_ + "->smoother", smootherArgs);

      return symmetricAggregation_
        ? createImpl3<Smoother>(linOp, SymCriterion(parameters_), smootherArgs, comm)
        : createImpl3<Smoother>(linOp, UnSymCriterion(parameters_), smootherArgs, comm);
    }

  private: // step 3:

    template <class Smoother, class LinOp, class Criterion, class Comm>
    std::unique_ptr<Interface>
    createImpl3(LinOp const& linOp, Criterion const& criterion,
                SmootherArgs<Smoother> const& smootherArgs, Comm const& comm) const
    {
      using Creator = AMGCreator<AMGSolver,Traits>;
      return Creator::template create<Smoother>(prefix_,linOp,criterion,smootherArgs,comm);
    }

  protected:
    void initParameters(std::string const& prefix)
    {
      // get creator string for preconditioner
      smoother_ = "default";
      Parameters::get(prefix + "->smoother", smoother_);

      // The debugging level. If 0 no debugging output will be generated.
      parameters_.setDebugLevel(Parameters::get<int>(prefix + "->debugLevel").value_or(2));
      // The number of presmoothing steps to apply
      parameters_.setNoPreSmoothSteps(Parameters::get<std::size_t>(prefix + "->preSmoothSteps").value_or(2));
      // The number of postsmoothing steps to apply
      parameters_.setNoPostSmoothSteps(Parameters::get<std::size_t>(prefix + "->postSmoothSteps").value_or(2));
      // Value of gamma; 1 for V-cycle, 2 for W-cycle
      parameters_.setGamma(Parameters::get<std::size_t>(prefix + "->gamma").value_or(1));
      // Whether to use additive multigrid.
      parameters_.setAdditive(Parameters::get<bool>(prefix + "->additive").value_or(false));
      // Whether to use symmetric or unsymmetric aggregation criterion
      symmetricAggregation_ = Parameters::get<bool>(prefix + "->symmetric aggregation").value_or(true);

      initCoarseningParameters(prefix + "->coarsening");
      initAggregationParameters(prefix + "->aggregation");
      initDependencyParameters(prefix + "->dependency");
    }

    void initCoarseningParameters(std::string const& prefix)
    {
      // The maximum number of levels allowed in the hierarchy.
      parameters_.setMaxLevel(Parameters::get<int>(prefix + "->maxLevel").value_or(100));
      // The maximum number of unknowns allowed on the coarsest level.
      parameters_.setCoarsenTarget(Parameters::get<int>(prefix + "->coarsenTarget").value_or(1000));
      // The minimum coarsening rate to be achieved.
      parameters_.setMinCoarsenRate(Parameters::get<double>(prefix + "->minCoarsenRate").value_or(1.2));
      // The damping factor to apply to the prolongated correction.
      parameters_.setProlongationDampingFactor(Parameters::get<double>(prefix + "->dampingFactor").value_or(1.6));
    }

    void initAggregationParameters(std::string const& prefix)
    {
      // Tthe maximal distance allowed between to nodes in a aggregate.
      parameters_.setMaxDistance(Parameters::get<std::size_t>(prefix + "->maxDistance").value_or(2));
      // The minimum number of nodes a aggregate has to consist of.
      parameters_.setMinAggregateSize(Parameters::get<std::size_t>(prefix + "->minAggregateSize").value_or(4));
      // The maximum number of nodes a aggregate is allowed to have.
      parameters_.setMaxAggregateSize(Parameters::get<std::size_t>(prefix + "->maxAggregateSize").value_or(6));
      // The maximum number of connections a aggregate is allowed to have.
      parameters_.setMaxConnectivity(Parameters::get<std::size_t>(prefix + "->maxConnectivity").value_or(15));
      // The maximum number of connections a aggregate is allowed to have.
      parameters_.setSkipIsolated(Parameters::get<bool>(prefix + "->skipIsolated").value_or(false));
    }

    void initDependencyParameters(std::string const& prefix)
    {
      // The scaling value for marking connections as strong.
      parameters_.setAlpha(Parameters::get<double>(prefix + "->alpha").value_or(1.0/3.0));
      // The threshold for marking nodes as isolated.
      parameters_.setBeta(Parameters::get<double>(prefix + "->beta").value_or(1.e-5));
    }

    template <class SA>
    void initSmootherParameters(std::string const& prefix, SA& smootherArgs) const
    {
      Parameters::get(prefix + "->iterations", smootherArgs.iterations);
      Parameters::get(prefix + "->relaxationFactor", smootherArgs.relaxationFactor);
    }

  private:
    std::string prefix_;
    std::string smoother_;
    Dune::Amg::Parameters parameters_;
    bool symmetricAggregation_ = true;

    mutable std::shared_ptr<LinOperator> linOperator_;
  };

} // end namespace AMDiS