Changed ISTL solver interface to use parallel solver; Added empty communicator...

Changed ISTL solver interface to use parallel solver; Added empty communicator for sequential use of parallel solvers

src/amdis/LinearAlgebra.hpp

@@ -10,6 +10,7 @@
 #include <amdis/linearalgebra/mtl/DOFVector.hpp>
 #include <amdis/linearalgebra/mtl/ITL_Solver.hpp>
 #include <amdis/linearalgebra/mtl/ITL_Preconditioner.hpp>
+#include <amdis/linearalgebra/mtl/Traits.hpp>
 
 #elif HAVE_EIGEN
 
@@ -17,6 +18,7 @@
 #include <amdis/linearalgebra/eigen/DOFMatrix.hpp>
 #include <amdis/linearalgebra/eigen/DOFVector.hpp>
 #include <amdis/linearalgebra/eigen/SolverCreator.hpp>
+#include <amdis/linearalgebra/eigen/Traits.hpp>
 
 #else // ISTL
 
@@ -25,5 +27,6 @@
 #include <amdis/linearalgebra/istl/DOFVector.hpp>
 #include <amdis/linearalgebra/istl/ISTL_Solver.hpp>
 #include <amdis/linearalgebra/istl/ISTL_Preconditioner.hpp>
+#include <amdis/linearalgebra/istl/Traits.hpp>
 
 #endif
\ No newline at end of file

src/amdis/ProblemStat.hpp

@@ -11,7 +11,6 @@
 #include <dune/common/fvector.hh>
 #include <dune/common/fmatrix.hh>
 #include <dune/grid/common/grid.hh>
-
 #include <amdis/AdaptInfo.hpp>
 #include <amdis/CreatorInterface.hpp>
 #include <amdis/CreatorMap.hpp>
@@ -72,8 +71,13 @@ namespace AMDiS
 
     using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, typename Traits::CoefficientType>;
     using SystemVector = DOFVector<GlobalBasis, typename Traits::CoefficientType>;
+    using LinearSolverTraits = SolverTraits<typename SystemMatrix::BaseMatrix,
+                                            typename SystemVector::BaseVector,
+                                            typename SystemVector::BaseVector,
+                                            GlobalBasis>;
 
-    using LinearSolverType = LinearSolverInterface<typename SystemMatrix::BaseMatrix, typename SystemVector::BaseVector>;
+    using LinearSolverType = LinearSolverInterface<LinearSolverTraits>;
+    using CommInfo = typename LinearSolverTraits::Comm;
 
   public:
     /**

src/amdis/ProblemStat.inc.hpp

@@ -309,10 +309,10 @@ solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
   solverInfo.setCreateMatrixData(createMatrixData);
   solverInfo.setStoreMatrixData(storeMatrixData);
 
-  solution_->compress();
+  auto commInfo = CommInfo::create(*globalBasis_, name_ + "->solver");
 
   linearSolver_->solve(systemMatrix_->matrix(), solution_->vector(), rhs_->vector(),
-                      solverInfo);
+                      *commInfo, solverInfo);
 
   if (solverInfo.info() > 0) {
     msg("solution of discrete system needed {} seconds", t.elapsed());

src/amdis/functions/GlobalIdSet.hpp

@@ -81,10 +81,23 @@ namespace AMDiS
     using GridView = typename GlobalBasis::GridView;
     using Grid = typename GridView::Grid;
     using Element = typename GridView::template Codim<0>::Entity;
-    using EntityIdType = typename Grid::GlobalIdSet::IdType; 
-    using IdType = std::pair<EntityIdType, std::size_t>;
+    using EntityIdType = typename Grid::GlobalIdSet::IdType;
     using PartitionType = Dune::PartitionType;
 
+    struct IdType : std::pair<EntityIdType, std::size_t>
+    {
+      using Super = std::pair<EntityIdType, std::size_t>;
+
+      IdType(int i = 0) : Super() {};
+      using Super::Super;
+
+      friend std::ostream& operator<<(std::ostream& os, IdType const& id)
+      {
+        os << "(" << id.first << "," << id.second << ")";
+        return os;
+      }
+    };
+
     using PreBasis = typename GlobalBasis::PreBasis;
     using TreePath = typename GlobalBasis::PrefixPath;
     using NodeIdSet = AMDiS::NodeIdSet<PreBasis, TreePath>;

src/amdis/linearalgebra/Common.hpp

 #pragma once
 
 #include <cstddef>
+#include <memory>
+#include <string>
 
 namespace AMDiS
 {
@@ -11,4 +13,30 @@ namespace AMDiS
     T value;
   };
 
+  template <class Basis>
+  class DefaultCommunication
+  {
+  public:
+    static std::unique_ptr<DefaultCommunication> create(Basis const& basis, std::string const& prefix)
+    {
+      DUNE_UNUSED_PARAMETER(basis);
+      DUNE_UNUSED_PARAMETER(prefix);
+      return std::make_unique<DefaultCommunication>();
+    }
+  };
+
+  /** Base traits class for a linear solver for the system AX=B using an FE space described by a
+   *  dune-functions Basis. This defines the general interface typedefs, all implementations are
+   *  required to provide the typedefs listed here, by e.g. inheriting from this.
+   */
+  template <class A, class X, class B, class Basis>
+  class SolverTraitsBase
+  {
+  public:
+    using Mat = A;
+    using Sol = X;
+    using Rhs = B;
+    using Comm = DefaultCommunication<Basis>;
+  };
+
 } // end namespace AMDiS

src/amdis/linearalgebra/LinearSolver.hpp

@@ -20,13 +20,18 @@ namespace AMDiS
    * solvers where the backend provides an interface, can be assigned
    * by different Runner objects.
    **/
-  template <class Mat, class Sol, class Rhs, class Runner>
+  template <class Traits, class Runner>
   class LinearSolver
-      : public LinearSolverInterface<Mat, Sol, Rhs>
+      : public LinearSolverInterface<Traits>
   {
     using Self = LinearSolver;
-    using Super = LinearSolverInterface<Mat, Sol, Rhs>;
+    using Super = LinearSolverInterface<Traits>;
+
     using RunnerBase = typename Super::RunnerBase;
+    using Mat = typename Traits::Mat;
+    using Sol = typename Traits::Sol;
+    using Rhs = typename Traits::Rhs;
+    using Comm = typename Traits::Comm;
 
   public:
     /// A creator to be used instead of the constructor.
@@ -52,12 +57,12 @@ namespace AMDiS
 
   private:
     /// Implements \ref LinearSolverInterface::solveSystemImpl()
-    void solveImpl(Mat const& A, Sol& x, Rhs const& b, SolverInfo& solverInfo) override
+    void solveImpl(Mat const& A, Sol& x, Rhs const& b, Comm& comm, SolverInfo& solverInfo) override
     {
       Dune::Timer t;
       if (solverInfo.doCreateMatrixData()) {
         // init matrix or wrap block-matrix or ...
-        runner_->init(A);
+        runner_->init(A, comm);
       }
 
       if (solverInfo.info() > 0)

src/amdis/linearalgebra/LinearSolverInterface.hpp

@@ -20,11 +20,15 @@
 namespace AMDiS
 {
   /// Abstract base class for linear solvers
-  template <class Mat, class Sol, class Rhs = Sol>
+  template <class Traits>
   class LinearSolverInterface
   {
   protected:
-    using RunnerBase = RunnerInterface<Mat, Sol, Rhs>;
+    using RunnerBase = RunnerInterface<Traits>;
+    using Mat = typename Traits::Mat;
+    using Sol = typename Traits::Sol;
+    using Rhs = typename Traits::Rhs;
+    using Comm = typename Traits::Comm;
 
   public:
     /// Destructor.
@@ -40,9 +44,9 @@ namespace AMDiS
      *  \p x     A [block-]vector for the unknown components.
      *  \p b     A [block-]vector for the right-hand side of the linear system.
      **/
-    void solve(Mat const& A, Sol& x, Rhs const& b, SolverInfo& solverInfo)
+    void solve(Mat const& A, Sol& x, Rhs const& b, Comm& comm, SolverInfo& solverInfo)
     {
-      solveImpl(A, x, b, solverInfo);
+      solveImpl(A, x, b, comm, solverInfo);
     }
 
     // return the runner/worker corresponding to this solver (optional)
@@ -50,7 +54,7 @@ namespace AMDiS
 
   private:
     /// main methods that all solvers must implement
-    virtual void solveImpl(Mat const& A, Sol& x, Rhs const& b, SolverInfo& solverInfo) = 0;
+    virtual void solveImpl(Mat const& A, Sol& x, Rhs const& b, Comm& comm, SolverInfo& solverInfo) = 0;
   };
 
 

src/amdis/linearalgebra/PreconditionerInterface.hpp

@@ -5,9 +5,13 @@
 namespace AMDiS
 {
   /// Interface for Preconditioner types
-  template <class Mat, class Sol, class Rhs = Sol>
+  template <class Traits>
   struct PreconditionerInterface
   {
+    using Mat = typename Traits::Mat;
+    using Sol = typename Traits::Sol;
+    using Rhs = typename Traits::Rhs;
+
     /// Virtual destructor.
     virtual ~PreconditionerInterface() = default;
 

src/amdis/linearalgebra/RunnerInterface.hpp

@@ -8,18 +8,22 @@ namespace AMDiS
   class SolverInfo;
 
   /// Interface for Runner / Worker types used in solver classes
-  template <class Mat, class Sol, class Rhs = Sol>
+  template <class Traits>
   class RunnerInterface
   {
   protected:
-    using PreconBase = PreconditionerInterface<Mat, Sol, Rhs>;
+    using PreconBase = PreconditionerInterface<Traits>;
+    using Mat = typename Traits::Mat;
+    using Sol = typename Traits::Sol;
+    using Rhs = typename Traits::Rhs;
+    using Comm = typename Traits::Comm;
 
   public:
     /// virtual destructor
     virtual ~RunnerInterface() = default;
 
     /// Is called at the beginning of a solution procedure
-    virtual void init(Mat const& A) = 0;
+    virtual void init(Mat const& A, Comm& comm) = 0;
 
     /// Is called at the end of a solution procedure
     virtual void exit() = 0;

src/amdis/linearalgebra/eigen/CMakeLists.txt

@@ -7,4 +7,5 @@ install(FILES
     PreconConfig.hpp
     SolverConfig.hpp
     SolverCreator.hpp
+    Traits.hpp
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/linearalgebra/eigen)

src/amdis/linearalgebra/eigen/DirectRunner.hpp

@@ -14,13 +14,17 @@ namespace AMDiS
   * \class AMDiS::DirectRunner
   * \brief \implements RunnerInterface for the (external) direct solvers
   */
-  template <template <class> class Solver, class Matrix, class VectorX, class VectorB>
+  template <template <class> class Solver, class Traits>
   class DirectRunner
-      : public RunnerInterface<Matrix, VectorX, VectorB>
+      : public RunnerInterface<Traits>
   {
   protected:
-    using Super = RunnerInterface<Matrix, VectorX, VectorB>;
-    using EigenSolver = Solver<Matrix>;
+    using Super = RunnerInterface<Traits>;
+    using Mat = typename Traits::Mat;
+    using Sol = typename Traits::Sol;
+    using Rhs = typename Traits::Rhs;
+    using Comm = typename Traits::Comm;
+    using EigenSolver = Solver<Mat>;
 
   public:
     /// Constructor.
@@ -31,9 +35,11 @@ namespace AMDiS
       Parameters::get(prefix + "->reuse pattern", reusePattern_);
     }
 
-    /// Implementes \ref RunnerInterface::init()
-    void init(Matrix const& A) override
+    /// Implements \ref RunnerInterface::init()
+    void init(Mat const& A, Comm& comm) override
     {
+      DUNE_UNUSED_PARAMETER(comm);
+
       if (!reusePattern_ || !initialized_) {
         solver_.analyzePattern(A);
         initialized_ = true;
@@ -45,16 +51,17 @@ namespace AMDiS
     }
 
 
-    /// Implementes \ref RunnerInterface::exit()
+    /// Implements \ref RunnerInterface::exit()
     void exit() override {}
 
-    /// Implementes \ref RunnerInterface::solve()
-    int solve(Matrix const& A, VectorX& x, VectorB const& b,
-              SolverInfo& solverInfo) override
+    /// Implements \ref RunnerInterface::solve()
+    int solve(Mat const& A, Sol& x, Rhs const& b, SolverInfo& solverInfo) override
     {
+      DUNE_UNUSED_PARAMETER(A);
+
       x = solver_.solve(b);
 
-      auto r = VectorB(b);
+      auto r = Rhs(b);
       if (x.norm() != 0)
         r -= A * x;
 

src/amdis/linearalgebra/eigen/IterativeRunner.hpp

@@ -9,12 +9,16 @@
 
 namespace AMDiS
 {
-  template <class IterativeSolver, class Matrix, class VectorX, class VectorB>
+  template <class IterativeSolver, class Traits>
   class IterativeRunner
-      : public RunnerInterface<Matrix, VectorX, VectorB>
+      : public RunnerInterface<Traits>
   {
     using SolverCfg = SolverConfig<IterativeSolver>;
     using PreconCfg = PreconConfig<typename IterativeSolver::Preconditioner>;
+    using Mat = typename Traits::Mat;
+    using Sol = typename Traits::Sol;
+    using Rhs = typename Traits::Rhs;
+    using Comm = typename Traits::Comm;
 
   public:
     IterativeRunner(std::string const& prefix)
@@ -26,9 +30,11 @@ namespace AMDiS
     }
 
 
-    /// Implementes \ref RunnerInterface::init()
-    void init(Matrix const& A) override
+    /// Implements \ref RunnerInterface::init()
+    void init(Mat const& A, Comm& comm) override
     {
+      DUNE_UNUSED_PARAMETER(comm);
+
       if (!reusePattern_ || !initialized_) {
         solver_.analyzePattern(A);
         initialized_ = true;
@@ -39,17 +45,18 @@ namespace AMDiS
         "Error in solver.compute(matrix)");
     }
 
-    /// Implementes \ref RunnerInterface::exit()
+    /// Implements \ref RunnerInterface::exit()
     void exit() override {}
 
-    /// Implementes \ref RunnerInterface::solve()
-    int solve(Matrix const& A, VectorX& x, VectorB const& b,
-              SolverInfo& solverInfo) override
+    /// Implements \ref RunnerInterface::solve()
+    int solve(Mat const& A, Sol& x, Rhs const& b, SolverInfo& solverInfo) override
     {
+      DUNE_UNUSED_PARAMETER(A);
+
       solver_.setTolerance(solverInfo.relTolerance());
       x = solver_.solveWithGuess(b, x);
 
-      auto r = VectorB(b);
+      auto r = Rhs(b);
       if (x.norm() != 0)
         r -= A * x;
 

src/amdis/linearalgebra/eigen/SolverCreator.hpp

@@ -11,24 +11,27 @@
 
 #include <amdis/CreatorMap.hpp>
 #include <amdis/Initfile.hpp>
+#include <amdis/Output.hpp>
 
 #include <amdis/linearalgebra/LinearSolver.hpp>
 #include <amdis/linearalgebra/eigen/DirectRunner.hpp>
 #include <amdis/linearalgebra/eigen/IterativeRunner.hpp>
+#include <amdis/linearalgebra/eigen/Traits.hpp>
 
 namespace AMDiS
 {
-  template <class Matrix, class VectorX, class VectorB, template <class> class IterativeSolver>
+  template <class Traits, template <class> class IterativeSolver>
   struct IterativeSolverCreator
-      : public CreatorInterfaceName<LinearSolverInterface<Matrix, VectorX, VectorB>>
+      : public CreatorInterfaceName<LinearSolverInterface<Traits>>
   {
     template <class Precon>
     using SolverCreator
-      = typename LinearSolver<Matrix, VectorX, VectorB,
-          IterativeRunner<IterativeSolver<Precon>, Matrix, VectorX, VectorB>>::Creator;
+      = typename LinearSolver<Traits,
+          IterativeRunner<IterativeSolver<Precon>, Traits>>::Creator;
 
-    using SolverBase = LinearSolverInterface<Matrix, VectorX, VectorB>;
-    using Scalar = typename Matrix::Scalar;
+    using SolverBase = LinearSolverInterface<Traits>;
+    using Mat = typename Traits::Mat;
+    using Scalar = typename Mat::Scalar;
 
     std::unique_ptr<SolverBase> createWithString(std::string prefix) override
     {
@@ -68,17 +71,16 @@ namespace AMDiS
       Parameters::get(prefix + "->precon->ordering", ordering);
 
       if (ordering == "amd") {
-        using AMD = Eigen::AMDOrdering<typename Matrix::StorageIndex>;
+        using AMD = Eigen::AMDOrdering<typename Mat::StorageIndex>;
         return IncompleteCholesky<AMD>{}.createWithString(prefix);
       }
       else {
-        using NATURAL = Eigen::NaturalOrdering<typename Matrix::StorageIndex>;
+        using NATURAL = Eigen::NaturalOrdering<typename Mat::StorageIndex>;
         return IncompleteCholesky<NATURAL>{}.createWithString(prefix);
       }
     }
   };
 
-
   /// Adds default creators for linear solvers based on `Eigen::SparseMatrix`.
   /**
    * Adds creators for full-matrix aware solvers.
@@ -90,36 +92,37 @@ namespace AMDiS
    * - *umfpack*: external UMFPACK solver, \see Eigen::UmfPackLU
    * - *superlu*: external SparseLU solver, \see Eigen::SparseLU
    **/
-  template <class T, int O, class VectorX, class VectorB>
-  class DefaultCreators<LinearSolverInterface<Eigen::SparseMatrix<T,O>, VectorX, VectorB>>
+  template <class Traits>
+  class DefaultCreators< LinearSolverInterface<Traits> >
   {
-    using Matrix = Eigen::SparseMatrix<T,O>;
-    using SolverBase = LinearSolverInterface<Matrix, VectorX, VectorB>;
+    using Mat = typename Traits::Mat;
+
+    using SolverBase = LinearSolverInterface<Traits>;
 
     template <template <class> class IterativeSolver>
     using SolverCreator
-      = IterativeSolverCreator<Matrix, VectorX, VectorB, IterativeSolver>;
+      = IterativeSolverCreator<Traits, IterativeSolver>;
 
     template <template <class> class DirectSolver>
     using DirectSolverCreator
-      = typename LinearSolver<Matrix, VectorX, VectorB,
-          DirectRunner<DirectSolver, Matrix, VectorX, VectorB>>::Creator;
+      = typename LinearSolver<Traits,
+          DirectRunner<DirectSolver, Traits>>::Creator;
 
     //@{
     template <class Precon>
-    using CG = Eigen::ConjugateGradient<Matrix, Eigen::Lower|Eigen::Upper, Precon>;
+    using CG = Eigen::ConjugateGradient<Mat, Eigen::Lower|Eigen::Upper, Precon>;
 
     template <class Precon>
-    using BCGS = Eigen::BiCGSTAB<Matrix, Precon>;
+    using BCGS = Eigen::BiCGSTAB<Mat, Precon>;
 
     template <class Precon>
-    using MINRES = Eigen::MINRES<Matrix, Eigen::Lower|Eigen::Upper, Precon>;
+    using MINRES = Eigen::MINRES<Mat, Eigen::Lower|Eigen::Upper, Precon>;
 
     template <class Precon>
-    using GMRES = Eigen::GMRES<Matrix, Precon>;
+    using GMRES = Eigen::GMRES<Mat, Precon>;
 
     template <class Precon>
-    using DGMRES = Eigen::DGMRES<Matrix, Precon>;
+    using DGMRES = Eigen::DGMRES<Mat, Precon>;
     // @}
 
     using Map = CreatorMap<SolverBase>;

src/amdis/linearalgebra/eigen/Traits.hpp

+#pragma once
+
+#include <amdis/linearalgebra/Common.hpp>
+
+namespace AMDiS
+{
+  /** Traits class for a linear solver for the system AX=B using an FE space described by a dune-functions Basis
+   *  Contains typedefs specific to the EIGEN backend.
+   */
+  template <class A, class X, class B, class Basis>
+  class SolverTraits
+      : public SolverTraitsBase<A,X,B,Basis>
+  {};
+
+} // end namespace AMDiS

src/amdis/linearalgebra/istl/CMakeLists.txt

 install(FILES
+    Communication.hpp
+    Communication.inc.hpp
     Constraints.hpp
+    Creators.hpp
     DirectRunner.hpp
     DOFMatrix.hpp
     DOFVector.hpp
@@ -7,4 +10,5 @@ install(FILES
     ISTL_Preconditioner.hpp
     ISTL_Solver.hpp
     ISTLRunner.hpp
+    Traits.hpp
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/linearalgebra/istl)

src/amdis/linearalgebra/istl/Communication.hpp

+#pragma once
+
+#include <memory>
+#include <ostream>
+#include <set>
+#include <string>
+#include <utility>
+
+#include <dune/common/version.hh>
+#include <dune/istl/owneroverlapcopy.hh>
+#include <dune/istl/solvercategory.hh>
+
+#include <amdis/functions/GlobalIdSet.hpp>
+
+namespace AMDiS
+{
+
+#if HAVE_MPI
+  /// Implementation class for ISTL-specific communication to be used in parallel solvers
+  template <class Basis>
+  class ISTLCommunication
+      : public Dune::OwnerOverlapCopyCommunication<typename GlobalBasisIdSet<Basis>::IdType, typename Basis::size_type>
+  {
+    using GlobalCommIdType = typename GlobalBasisIdSet<Basis>::IdType;
+    using LocalCommIdType = typename Basis::size_type;
+
+    using Self = ISTLCommunication;
+    using Super = Dune::OwnerOverlapCopyCommunication<GlobalCommIdType, LocalCommIdType>;
+
+    using GridView = typename Basis::GridView;
+    using Grid = typename GridView::Grid;
+    using Element = typename GridView::template Codim<0>::Entity;
+    using IdSet = typename Grid::GlobalIdSet;
+    using IdType = typename IdSet::IdType;
+    using EntitySet = std::set<IdType>;
+    using Comm = typename Dune::MPIHelper::MPICommunicator;
+    using SolverType = typename Dune::SolverCategory::Category;
+
+  public:
+    /**
+     * Constructor. Takes the global basis, an MPI communicator and the ISTL solver type
+     *
+     * \param basis   Global basis associated to the solvers' matrices and vectors
+     * \param comm    MPI communicator of the group containing the distributed matrices and vectors
+     * \param cat     Category of the solver used with this communicator, one of
+     *                Dune::SolverCategory::Category::sequential, overlapping, nonoverlapping
+     **/
+    ISTLCommunication(Basis const& basis, Comm const& comm, SolverType cat)
+      : Super(comm, cat)
+    {
+      if (cat != SolverType::sequential)
+      {
+        auto& pis = this->indexSet();
+        buildParallelIndexSet(basis, pis);
+        auto& ris = this->remoteIndices();
+        ris.setIndexSets(pis,pis,comm);
+        ris.template rebuild<true>();
+      }
+    }
+
+  private:
+    /// Builds the parallel index set containing all DoFs of the basis
+    void buildParallelIndexSet(Basis const& basis, typename Super::PIS& pis);
+
+  public:
+    /**
+     * ISTLCommunication creator function. Takes the global basis and the solver's initfile prefix
+     * and returns a unique pointer to the ISTLCommunication object.
+     * Deduces the category of the solver used from the initfile and basis
+     *
+     * \param basis   Global basis associated to the solvers' matrices and vectors
+     * \param prefix  Solver prefix used in the initfile, usually "<problem name>->solver"
+     **/
+    static std::unique_ptr<Self> create(Basis const& basis, std::string const& prefix);
+  };
+
+#else // !HAVE_MPI
+
+  /// Dummy implementation for ISTL-specific communication when no MPI is found
+  template <class Basis>
+  class ISTLCommunication
+      : public DefaultCommunication<Basis>
+  {
+    using Self = ISTLCommunication;
+
+  public:
+    ISTLCommunication() = default;
+
+    Dune::SolverCategory::Category category() const