istl AMG preconditioner added

The dune-istl AMG method is available as precondition. The corresponding initfile parameters are

[...]->solver->precon: [amg,fastamg]
[...]->solver->precon->smoother: [jacobi,sor,ssor] % later more smoothers can be added

[...]->solver->precon->[...] % general amg parameters, see Dune::Amg::Parameters
[...]->solver->precon->coarsening->[...] % Parameters of Dune::Amg::CoarseningParameters
[...]->solver->precon->aggregation->[...] % Parameters of Dune::Amg::AggregationParameters
[...]->solver->precon->dependency->[...] % Parameters of Dune::Amg::DependencyParameters
[...]->solver->precon->smoother->[...] % Parameters of Dune::Amg::DefaultSmootherArgs

See amdis/linearalgebra/istl/AMGPrecon.hpp for a detailed list of parameters. The parameter names follow the variable names in the Dune Parameters class.

src/amdis/linearalgebra/istl/AMGPrecon.hpp

+#pragma once
+
+#include <dune/istl/paamg/amg.hh>
+#include <dune/istl/paamg/fastamg.hh>
+
+#include <amdis/linearalgebra/istl/ISTLPreconCreatorInterface.hpp>
+
+namespace AMDiS
+{
+  template <template <class...> class AMGSolver, class Mat, class Sol, class Rhs>
+  struct AMGTraits;
+
+  template <class Mat, class Sol, class Rhs>
+  struct AMGTraits<Dune::Amg::AMG, Mat, Sol, Rhs>
+  {
+    using Operator = Dune::MatrixAdapter<Mat, Sol, Rhs>;
+
+    template <class Smoother>
+    using Solver = Dune::Amg::AMG<Operator, Sol, Smoother>;
+
+    template <class Smoother, class Criterion, class SmootherArgs>
+    static auto create(Operator const& fineOperator, Criterion const& criterion, SmootherArgs const& smootherArgs)
+    {
+      return std::make_unique<Solver<Smoother>>(fineOperator, criterion, smootherArgs);
+    }
+  };
+
+  template <class Mat, class Sol, class Rhs>
+  struct AMGTraits<Dune::Amg::FastAMG, Mat, Sol, Rhs>
+  {
+    using Operator = Dune::MatrixAdapter<Mat, Sol, Rhs>;
+
+    template <class>
+    using Solver = Dune::Amg::FastAMG<Operator, Sol>;
+
+    template <class Smoother, class Criterion, class SmootherArgs>
+    static auto create(Operator const& fineOperator, Criterion const& criterion, SmootherArgs const& /*smootherArgs*/)
+    {
+      return std::make_unique<Solver<Smoother>>(fineOperator, criterion, criterion);
+    }
+  };
+
+
+  template <class Traits, class Smoother, class Mat, class Sol, class Rhs>
+  class AMGPrecon;
+
+  template <template <class...> class AMGSolver, class Mat, class Sol, class Rhs>
+  class AMGPreconCreator
+      : public CreatorInterfaceName<ISTLPreconCreatorInterface<Mat, Sol, Rhs>>
+  {
+    using PreconCreatorBase = ISTLPreconCreatorInterface<Mat, Sol, Rhs>;
+    using Traits = AMGTraits<AMGSolver,Mat,Sol,Rhs>;
+
+    template <class Smoother>
+    using Precon = AMGPrecon<Traits, Smoother, Mat, Sol, Rhs>;
+
+  public:
+    std::unique_ptr<PreconCreatorBase> createWithString(std::string prefix) override
+    {
+      // get creator string for preconditioner
+      std::string smoother = "no";
+      Parameters::get(prefix + "->smoother", smoother);
+
+      if (smoother == "diag" ||
+          smoother == "jacobi")
+      {
+        auto creator = typename Precon<Dune::SeqJac<Mat, Sol, Rhs>>::Creator{};
+        return creator.createWithString(prefix);
+      }
+      // else if (smoother == "gs" ||
+      //          smoother == "hauss_seidel")
+      // {
+      //   auto creator = typename Precon<Dune::SeqGS<Mat, Sol, Rhs>>::Creator{};
+      //   return creator.createWithString(prefix);
+      // }
+      else if (smoother == "sor")
+      {
+        auto creator = typename Precon<Dune::SeqSOR<Mat, Sol, Rhs>>::Creator{};
+        return creator.createWithString(prefix);
+      }
+      else if (smoother == "ssor") {
+        auto creator = typename Precon<Dune::SeqSSOR<Mat, Sol, Rhs>>::Creator{};
+        return creator.createWithString(prefix);
+      }
+      // else if (smoother == "richardson" ||
+      //          smoother == "default") {
+      //   auto creator = typename Precon<Dune::Richardson<Sol, Rhs>>::Creator{};
+      //   return creator.createWithString(prefix);
+      // }
+      else {
+        error_exit("Unknown smoother type '{}' given for parameter '{}'", smoother, prefix + "->smoother");
+        return nullptr;
+      }
+    }
+  };
+
+
+  template <class Traits, class Smoother, class Mat, class Sol, class Rhs>
+  class AMGPrecon
+      : public ISTLPreconCreatorInterface<Mat, Sol, Rhs>
+  {
+    using Super = ISTLPreconCreatorInterface<Mat, Sol, Rhs>;
+    using Self = AMGPrecon;
+
+    using PreconBase = Dune::Preconditioner<Sol, Rhs>;
+
+    using LinOperator = typename Traits::Operator;
+    using SmootherArgs = typename Dune::Amg::SmootherTraits<Smoother>::Arguments;
+
+    using Norm = std::conditional_t<std::is_arithmetic<typename Dune::FieldTraits<Mat>::field_type>::value,
+      Dune::Amg::FirstDiagonal, Dune::Amg::RowSum>;
+    // TODO: make criterion flexible
+    using Criterion = Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<Mat,Norm>>;
+
+  public:
+    struct Creator : CreatorInterfaceName<Super>
+    {
+      std::unique_ptr<Super> createWithString(std::string prefix) override
+      {
+        return std::make_unique<Self>(prefix);
+      }
+    };
+
+    AMGPrecon(std::string const& prefix)
+    {
+      initParameters(prefix);
+    }
+
+    /// Implementes \ref ISTLPreconCreatorInterface::create
+    std::unique_ptr<PreconBase> create(Mat const& A) const override
+    {
+      linOperator_ = std::make_shared<LinOperator>(A);
+      criterion_ = std::make_shared<Criterion>(parameters_);
+
+      return Traits::template create<Smoother>(*linOperator_, *criterion_, smootherArgs_);
+    }
+
+  protected:
+    void initParameters(std::string const& prefix)
+    {
+      // 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));
+
+      initCoarseningParameters(prefix + "->coarsening");
+      initAggregationParameters(prefix + "->aggregation");
+      initDependencyParameters(prefix + "->dependency");
+      initSmootherParameters(prefix + "->smoother");
+    }
+
+    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 prologated 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));
+    }
+
+    void initSmootherParameters(std::string const& prefix)
+    {
+      Parameters::get(prefix + "->iterations", smootherArgs_.iterations);
+      Parameters::get(prefix + "->relaxationFactor", smootherArgs_.relaxationFactor);
+    }
+
+  private:
+    SmootherArgs smootherArgs_;
+    Dune::Amg::Parameters parameters_;
+
+    mutable std::shared_ptr<Criterion>    criterion_;
+    mutable std::shared_ptr<LinOperator>  linOperator_;
+  };
+
+} // end namespace AMDiS