Merge branch 'feature/istl_amg' into 'master'

istl AMG preconditioner added

See merge request !32

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

src/amdis/linearalgebra/istl/ISTLPreconCreatorInterface.hpp

+#pragma once
+
+#include <memory>
+
+#include <dune/istl/preconditioners.hh>
+
+namespace AMDiS
+{
+  template <class Mat, class Sol, class Rhs>
+  struct ISTLPreconCreatorInterface
+  {
+    virtual ~ISTLPreconCreatorInterface() = default;
+
+    using PreconBase = Dune::Preconditioner<Sol, Rhs>;
+    virtual std::unique_ptr<PreconBase> create(Mat const& matrix) const = 0;
+  };
+
+} // end namespace AMDiS
\ No newline at end of file

src/amdis/linearalgebra/istl/ISTLRunner.hpp

@@ -5,7 +5,7 @@
 #include <amdis/linearalgebra/RunnerInterface.hpp>
 #include <amdis/linearalgebra/SolverInfo.hpp>
 #include <amdis/linearalgebra/istl/Fwd.hpp>
-#include <amdis/linearalgebra/istl/ISTL_Preconditioner.hpp>
+#include <amdis/linearalgebra/istl/ISTLPreconCreatorInterface.hpp>
 
 namespace AMDiS
 {
@@ -16,7 +16,7 @@ namespace AMDiS
     using LinOperator = Dune::MatrixAdapter<Matrix, VectorX, VectorB>;
     using BaseSolver  = Dune::InverseOperator<VectorX, VectorB>;
     using BasePrecon  = Dune::Preconditioner<VectorX, VectorB>;
-    using ISTLPrecon  = ISTLPreconInterface<Matrix, VectorX, VectorB>;
+    using ISTLPreconCreator  = ISTLPreconCreatorInterface<Matrix, VectorX, VectorB>;
 
   public:
     ISTLRunner(std::string const& prefix)
@@ -64,10 +64,10 @@ namespace AMDiS
       // get creator string for preconditioner
       std::string preconName = "no";
       std::string initFileStr = "";
-      for (std::string postfix : {"left precon", "right precon", "precon->name"}) {
+      for (std::string postfix : {"left precon", "right precon", "precon", "precon->name"}) {
         Parameters::get(prefix + "->" + postfix, preconName);
         if (!preconName.empty() && preconName != "no") {
-          initFileStr = prefix + "->" + postfix;
+          initFileStr = postfix == "precon->name" ? prefix + "->precon" : prefix + "->" + postfix;
           break;
         }
       }
@@ -75,14 +75,14 @@ namespace AMDiS
       if (preconName.empty() || preconName == "no")
         preconName = "default";
 
-      auto* creator = named(CreatorMap<ISTLPrecon>::getCreator(preconName, initFileStr));
+      auto* creator = named(CreatorMap<ISTLPreconCreator>::getCreator(preconName, initFileStr));
       preconCreator_ = creator->createWithString(initFileStr);
       assert(preconCreator_);
     }
 
   private:
     ISTLSolverCreator<ISTLSolver> solverCreator_;
-    std::shared_ptr<ISTLPrecon> preconCreator_;
+    std::shared_ptr<ISTLPreconCreator> preconCreator_;
 
     std::shared_ptr<BasePrecon>   precon_;
     std::shared_ptr<LinOperator>  linOperator_;

src/amdis/linearalgebra/istl/ISTL_Preconditioner.hpp

@@ -3,30 +3,25 @@
 #include <dune/istl/preconditioners.hh>
 
 #include <amdis/CreatorInterface.hpp>
+#include <amdis/linearalgebra/istl/AMGPrecon.hpp>
+#include <amdis/linearalgebra/istl/Fwd.hpp>
+#include <amdis/linearalgebra/istl/ISTLPreconCreatorInterface.hpp>
 
 namespace AMDiS
 {
-  template <class Matrix, class VectorX, class VectorB>
-  struct ISTLPreconInterface
-  {
-    virtual ~ISTLPreconInterface() = default;
-    using PreconBase = Dune::Preconditioner<VectorX, VectorB>;
-    virtual std::unique_ptr<PreconBase> create(Matrix const& matrix) const = 0;
-  };
-
   template <class> struct Type {};
 
   // creators for preconditioners, depending on matrix-type
   // ---------------------------------------------------------------------------
 
   template <class Precon, class Matrix>
-  struct ISTLPrecon
-      : public ISTLPreconInterface<Matrix, typename Precon::domain_type, typename Precon::range_type>
+  struct ISTLPreconCreator
+      : public ISTLPreconCreatorInterface<Matrix, typename Precon::domain_type, typename Precon::range_type>
   {
-    using VectorX = typename Precon::domain_type;
-    using VectorB = typename Precon::range_type;
-    using Super = ISTLPreconInterface<Matrix, VectorX, VectorB>;
-    using Self = ISTLPrecon;
+    using Sol = typename Precon::domain_type;
+    using Rhs = typename Precon::range_type;
+    using Super = ISTLPreconCreatorInterface<Matrix, Sol, Rhs>;
+    using Self = ISTLPreconCreator;
 
     struct Creator : CreatorInterfaceName<Super>
     {
@@ -36,13 +31,13 @@ namespace AMDiS
       }
     };
 
-    ISTLPrecon(std::string const& prefix)
+    ISTLPreconCreator(std::string const& prefix)
     {
       Parameters::get(prefix + "->relaxation", w_);
       Parameters::get(prefix + "->iterations", iter_);
     }
 
-    using PreconBase = Dune::Preconditioner<VectorX, VectorB>;
+    using PreconBase = Dune::Preconditioner<Sol, Rhs>;
     std::unique_ptr<PreconBase> create(Matrix const& A) const override
     {
       return createImpl(A, Type<Precon>{});
@@ -55,16 +50,16 @@ namespace AMDiS
       return std::make_unique<P>(A, iter_, w_);
     }
 
-    std::unique_ptr<Dune::SeqILU<Matrix, VectorX, VectorB>>
-    createImpl(Matrix const& A, Type<Dune::SeqILU<Matrix, VectorX, VectorB>>) const
+    std::unique_ptr<Dune::SeqILU<Matrix, Sol, Rhs>>
+    createImpl(Matrix const& A, Type<Dune::SeqILU<Matrix, Sol, Rhs>>) const
     {
-      return std::make_unique<Dune::SeqILU<Matrix, VectorX, VectorB>>(A, iter_, w_);
+      return std::make_unique<Dune::SeqILU<Matrix, Sol, Rhs>>(A, iter_, w_);
     }
 
-    std::unique_ptr<Dune::Richardson<VectorX, VectorB>>
-    createImpl(Matrix const& /*A*/, Type<Dune::Richardson<VectorX, VectorB>>) const
+    std::unique_ptr<Dune::Richardson<Sol, Rhs>>
+    createImpl(Matrix const& /*A*/, Type<Dune::Richardson<Sol, Rhs>>) const
     {
-      return std::make_unique<Dune::Richardson<VectorX, VectorB>>(w_);
+      return std::make_unique<Dune::Richardson<Sol, Rhs>>(w_);
     }
 
     double w_ = 1.0;
@@ -81,37 +76,42 @@ namespace AMDiS
    * - *gmres*: Generalized minimal residula method, \see Dune::RestartedGMResSolver
    * - *umfpack*: external UMFPACK solver, \see Dune::UMFPack
    **/
-  template <class Matrix, class VectorX, class VectorB>
-  class DefaultCreators<ISTLPreconInterface<Matrix, VectorX, VectorB>>
+  template <class Mat, class Sol, class Rhs>
+  class DefaultCreators<ISTLPreconCreatorInterface<Mat, Sol, Rhs>>
   {
     template <class Precon>
     using PreconCreator
-      = typename ISTLPrecon<Precon, Matrix>::Creator;
+      = typename ISTLPreconCreator<Precon, Mat>::Creator;
 
-    using Map = CreatorMap<ISTLPreconInterface<Matrix, VectorX, VectorB>>;
+    using Map = CreatorMap<ISTLPreconCreatorInterface<Mat, Sol, Rhs>>;
 
   public:
     static void init()
     {
-      auto jacobi = new PreconCreator<Dune::SeqJac<Matrix, VectorX, VectorB>>;
+      auto jacobi = new PreconCreator<Dune::SeqJac<Mat, Sol, Rhs>>;
       Map::addCreator("diag", jacobi);
       Map::addCreator("jacobi", jacobi);
 
-      auto gs = new PreconCreator<Dune::SeqGS<Matrix, VectorX, VectorB>>;
+      auto gs = new PreconCreator<Dune::SeqGS<Mat, Sol, Rhs>>;
       Map::addCreator("gs", gs);
       Map::addCreator("gauss_seidel", gs);
 
-      auto sor = new PreconCreator<Dune::SeqSOR<Matrix, VectorX, VectorB>>;
+      auto sor = new PreconCreator<Dune::SeqSOR<Mat, Sol, Rhs>>;
       Map::addCreator("sor", sor);
 
-      auto ssor = new PreconCreator<Dune::SeqSSOR<Matrix, VectorX, VectorB>>;
+      auto ssor = new PreconCreator<Dune::SeqSSOR<Mat, Sol, Rhs>>;
       Map::addCreator("ssor", ssor);
 
-      auto ilu = new PreconCreator<Dune::SeqILU<Matrix, VectorX, VectorB>>;
+      auto ilu = new PreconCreator<Dune::SeqILU<Mat, Sol, Rhs>>;
       Map::addCreator("ilu", ilu);
       Map::addCreator("ilu0", ilu);
 
-      auto richardson = new PreconCreator<Dune::Richardson<VectorX, VectorB>>;
+      auto amg = new AMGPreconCreator<Dune::Amg::AMG, Mat, Sol, Rhs>;
+      Map::addCreator("amg", amg);
+      auto fastamg = new AMGPreconCreator<Dune::Amg::FastAMG, Mat, Sol, Rhs>;
+      Map::addCreator("fastamg", fastamg);
+
+      auto richardson = new PreconCreator<Dune::Richardson<Sol, Rhs>>;
       Map::addCreator("richardson", richardson);
       Map::addCreator("default", richardson);
     }