Add back the AMGPrecon

amdis/linearalgebra/istl/AMGPrecon.hpp

@@ -14,7 +14,7 @@
 #include <amdis/Initfile.hpp>
 #include <amdis/common/ConceptsBase.hpp>
 #include <amdis/linearalgebra/istl/Communication.hpp>
-#include <amdis/linearalgebra/istl/CreatorInterfaces.hpp>
+#include <amdis/linearalgebra/istl/ISTLPreconCreator.hpp>
 #include <amdis/linearalgebra/istl/precompiled/Preconditioners.hpp>
 #include <amdis/linearalgebra/istl/precompiled/Solvers.hpp>
 
@@ -56,10 +56,21 @@ namespace AMDiS
   {
     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>
-    create(std::string prefix, LinOp const& linOp, Criterion const& criterion, [[maybe_unused]] SmootherArgs<Smoother> const& smootherArgs,
-           Dune::Amg::SequentialInformation const& comm)
+    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);
 
@@ -67,10 +78,10 @@ namespace AMDiS
       return std::make_unique<Solver>(linOp, criterion, criterion, symmetric, comm);
     }
 
-    template <class Smoother, class LinOp, class Criterion, class SmootherArgs, class Comm,
-      REQUIRES(not std::is_same_v<Comm, Dune::Amg::SequentialInformation>)>
+    template <class Smoother, class LinOp, class Criterion, class SmootherArgs, class Comm>
     static std::unique_ptr<PrecBase>
-    create(std::string, LinOp const&, Criterion const&, SmootherArgs const&, Comm const&)
+    createImpl1(std::string, LinOp const&, Criterion const&, SmootherArgs const&, Comm const&,
+                Dune::PriorityTag<1>)
     {
       error_exit("Currently only sequential FastAMG supported.");
       return nullptr;
@@ -96,7 +107,8 @@ namespace AMDiS
 
     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)
+    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");
 
@@ -105,8 +117,7 @@ namespace AMDiS
       double minDefectReduction = 1.e-1;
       Parameters::get(prefix + "->krylov solver->minDefectReduction", minDefectReduction);
 
-      if (solver == "pcg" ||
-          solver == "default")
+      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);
@@ -123,8 +134,7 @@ namespace AMDiS
         return std::make_unique<Solver>(linOp, criterion, smootherArgs, maxLevelKrylovSteps, minDefectReduction, comm);
       }
 #endif
-      else if (solver == "bicgstab" ||
-               solver == "bcgs")
+      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);
@@ -152,82 +162,7 @@ namespace AMDiS
   };
 
 
-  template <class Creator, class Smoother, class Traits>
-  class AMGPrecon;
-
-  /// A creator for \ref AMGPrecon, reads the smoother type from initfile:
-  /**
-   * Initfile parameters:
-   * - `[PRECON]->smoother:  The Smoother type to use. [jacobi]
-   *
-   * Note: possible values for the smoother types: [diag|jacobi], ssor.
-   */
-  template <template <class...> class AMGSolver, class Traits>
-  class AMGPreconCreator
-      : public CreatorInterfaceName<ISTLPreconCreatorInterface<Traits>>
-  {
-    using M = typename Traits::M;
-    using X = typename Traits::X;
-    using Y = typename Traits::Y;
-
-    using PreconCreatorBase = ISTLPreconCreatorInterface<Traits>;
-    using Creator = AMGCreator<AMGSolver,Traits>;
-
-    template <class Smoother>
-    using Precon = AMGPrecon<Creator, Smoother, Traits>;
-
-  public:
-    std::unique_ptr<PreconCreatorBase> createWithString(std::string prefix) override
-    {
-      // get creator string for preconditioner
-      std::string smoother = "default";
-      Parameters::get(prefix + "->smoother", smoother);
-
-      if (smoother == "diag" ||
-          smoother == "jacobi" ||
-          smoother == "default")
-      {
-        auto creator = typename Precon<Dune::SeqJac<M,X,Y>>::Creator{};
-        return creator.createWithString(prefix);
-      }
-#if 0
-      // NOTE: apply<forward> not implemented correctly in BlockPreconditioner and
-      // NonoverlappingBlockPreconditioner. See !302 in dune-istl
-      else if (smoother == "sor")
-      {
-        auto creator = typename Precon<Dune::SeqSOR<M,X,Y>>::Creator{};
-        return creator.createWithString(prefix);
-      }
-#endif
-#if 0
-      // NOTE: ConstructionTraits not implemented for SeqGS. See !303 in dune-istl
-      else if (smoother == "gs" ||
-               smoother == "gauss_seidel")
-      {
-        auto creator = typename Precon<Dune::SeqGS<M,X,Y>>::Creator{};
-        return creator.createWithString(prefix);
-      }
-#endif
-#if 0
-      // NOTE: ConstructionTraits not implemented for Richardson. See !304 in dune-istl
-      else if (smoother == "richardson") {
-        auto creator = typename Precon<Dune::Richardson<X,Y>>::Creator{};
-        return creator.createWithString(prefix);
-      }
-#endif
-      else if (smoother == "ssor") {
-        auto creator = typename Precon<Dune::SeqSSOR<M,X,Y>>::Creator{};
-        return creator.createWithString(prefix);
-      }
-      else {
-        error_exit("Unknown smoother type '{}' given for parameter '{}'", smoother, prefix + "->smoother");
-        return nullptr;
-      }
-    }
-  };
-
-
-  /// Implementation of \ref ISTLPreconCreatorInterface to be used in the
+  /// Implementation of \ref ISTLPreconCreatorBase to be used in the
   /// \ref DefaultCreators.
   /**
    * Read several parameters of the AMG preconditioner and calls the corresponding
@@ -249,17 +184,13 @@ namespace AMDiS
    * - `[PRECON]->aggregation->...`: Parameters for the coarsening by aggregation
    * - `[PRECON]->dependency->...`: Parameters for the characterization for variable dependencies.
    */
-  template <class AMGCreator, class Smoother, class Traits>
-  class AMGPrecon
-      : public ISTLPreconCreatorInterface<Traits>
+  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 Comm = typename Traits::Comm;
-
-    using Super = ISTLPreconCreatorInterface<Traits>;
-    using Self = AMGPrecon;
 
     using Interface = Dune::Preconditioner<X,Y>;
     using LinOperator = typename Traits::LinOpCreator::Interface;
@@ -273,36 +204,79 @@ namespace AMDiS
     using UnSymCriterion = Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<M,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)
-      : prefix_(prefix)
+    void init(std::string const& prefix) override
     {
+      prefix_ = prefix;
       initParameters(prefix);
     }
 
-    /// Implements \ref ISTLPreconCreatorInterface::create
-    std::unique_ptr<Interface> create(M const& mat, typename Traits::Comm const& comm) const override
+    /// Implements \ref ISTLPreconCreatorBase::create
+    std::unique_ptr<Interface>
+    create(M const& mat, typename Traits::Comm const& comm) const override
     {
-      return createImpl1(Dune::SolverCategory::category(comm),mat,comm.get());
+      return createImpl0(Dune::SolverCategory::category(comm), mat, comm.get());
     }
 
-  private:
-    template <class Comm,
-      REQUIRES(not std::is_same_v<Comm, Dune::Amg::SequentialInformation>)>
+  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(SolverCategory cat, M const& mat, Comm const& comm) const
+    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(cat, mat, Dune::Amg::SequentialInformation{});
+          return createImpl1<Smoother>(cat, mat, Dune::Amg::SequentialInformation{}, Dune::PriorityTag<5>{});
         }
 #if HAVE_MPI
         case SolverCategory::overlapping:
@@ -328,31 +302,37 @@ namespace AMDiS
       }
     }
 
-    std::unique_ptr<Interface>
-    createImpl1([[maybe_unused]] SolverCategory cat, M const& mat, Dune::Amg::SequentialInformation const& comm) 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);
-    }
+  private: // step 2:
 
-  private:
-    template <class S, class LinOp, class Comm>
+    template <class Smoother, class LinOp, class Comm>
     std::unique_ptr<Interface> createImpl2(LinOp const& linOp, Comm const& comm) const
     {
-      typename Dune::Amg::SmootherTraits<S>::Arguments smootherArgs;
+      SmootherArgs<Smoother> smootherArgs;
       initSmootherParameters(prefix_ + "->smoother", smootherArgs);
 
       return symmetricAggregation_
-        ? AMGCreator::template create<S>(prefix_, linOp, SymCriterion(parameters_), smootherArgs, comm)
-        : AMGCreator::template create<S>(prefix_, linOp, UnSymCriterion(parameters_), smootherArgs, comm);
+        ? 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
@@ -414,6 +394,7 @@ namespace AMDiS
 
   private:
     std::string prefix_;
+    std::string smoother_;
     Dune::Amg::Parameters parameters_;
     bool symmetricAggregation_ = true;
 

amdis/linearalgebra/istl/CMakeLists.txt

@@ -3,7 +3,7 @@ target_sources(amdis PRIVATE
     Solvers.cpp)
 
 install(FILES
-    #AMGPrecon.hpp
+    AMGPrecon.hpp
     Communication.hpp
     Communication.inc.hpp
     Constraints.hpp

amdis/linearalgebra/istl/Preconditioners.hpp

@@ -5,7 +5,7 @@
 #include <amdis/CreatorInterface.hpp>
 #include <amdis/CreatorMap.hpp>
 #include <amdis/common/TypeTraits.hpp>
-// #include <amdis/linearalgebra/istl/AMGPrecon.hpp>
+#include <amdis/linearalgebra/istl/AMGPrecon.hpp>
 #include <amdis/linearalgebra/istl/ISTLPreconCreator.hpp>
 #include <amdis/linearalgebra/istl/precompiled/Preconditioners.hpp>
 
@@ -57,7 +57,7 @@ namespace AMDiS
       Map::addCreator("ssor", ssor);
 
       init_ilu(std::is_arithmetic<typename FTraits::field_type>{});
-      // init_amg(std::is_same<typename FTraits::real_type, double>{});
+      init_amg(std::is_same<typename FTraits::real_type, double>{});
 
       auto richardson = new Creator<Dune::Richardson<X,Y>>;
       Map::addCreator("richardson", richardson);
@@ -85,21 +85,21 @@ namespace AMDiS
       Map::addCreator("ildl", ildl);
     }
 
-    // static void init_amg(std::false_type)
-    // {
-    //   warning("AMG preconditioners not created for the matrix with real_type = {}.",
-    //     Dune::className<typename FTraits::real_type>());
-    // }
-
-    // static void init_amg(std::true_type)
-    // {
-    //   auto amg = new AMGPreconCreator<Dune::Amg::AMG, Traits>;
-    //   Map::addCreator("amg", amg);
-    //   auto fastamg = new AMGPreconCreator<Dune::Amg::FastAMG, Traits>;
-    //   Map::addCreator("fastamg", fastamg);
-    //   auto kamg = new AMGPreconCreator<Dune::Amg::KAMG, Traits>;
-    //   Map::addCreator("kamg", kamg);
-    // }
+    static void init_amg(std::false_type)
+    {
+      warning("AMG preconditioners not created for the matrix with real_type = {}.",
+        Dune::className<typename FTraits::real_type>());
+    }
+
+    static void init_amg(std::true_type)
+    {
+      auto amg = new AMGPreconCreator<Dune::Amg::AMG, Traits>;
+      Map::addCreator("amg", amg);
+      auto fastamg = new AMGPreconCreator<Dune::Amg::FastAMG, Traits>;
+      Map::addCreator("fastamg", fastamg);
+      auto kamg = new AMGPreconCreator<Dune::Amg::KAMG, Traits>;
+      Map::addCreator("kamg", kamg);
+    }
 
     static void init_bjacobi(Types<Dune::Amg::SequentialInformation>, Dune::PriorityTag<2>) {}
 

amdis/linearalgebra/petsc/LinearSolver.hpp

@@ -90,7 +90,7 @@ namespace AMDiS
     }
 
     /// Implements \ref Dune::InverOperator::apply()
-    void apply(Vector& x, Vector& b, Dune::InverseOperatorResult& stat) override
+    void apply(Vector& x, Vector const& b, Dune::InverseOperatorResult& stat) override
     {
       KSPSolve(ksp_, b.vector(), x.vector());