solvers are reworked and a lot of thingt corrected in the parallel petsc backend

AMDiS/src/AMDiS.cc

@@ -29,6 +29,10 @@
 #endif
 #include "boost/program_options.hpp"
 
+#if defined HAVE_PETSC || defined HAVE_SEQ_PETSC || defined HAVE_PARALLEL_PETSC
+#include <petsc.h>
+#endif
+
 namespace AMDiS {
 
   using namespace std;

AMDiS/src/AMDiS.h

@@ -154,9 +154,7 @@
 #include "parallel/StdMpi.h"
 #include "parallel/ParallelProblemStat.h"
 
-#if HAVE_PARALLEL_MTL4
-// #include "parallel/PMTL4Solver.h"
-#else
+#if !HAVE_PARALLEL_MTL4
 #include "parallel/PetscSolver.h"
 #include "parallel/PetscSolverNavierStokes.h"
 #endif

AMDiS/src/AMDiS_fwd.h

@@ -35,6 +35,7 @@ namespace AMDiS {
   class AdaptStationary;
   class Assembler;
   class BasisFunction; 
+  struct BlockMTLMatrix;
   class BoundaryManager;
 //   class CGSolver;
   class CoarseningManager;

AMDiS/src/CreatorMap.cc

@@ -25,6 +25,7 @@
 #include "MTL4Types.h"
 #include "solver/LinearSolver.h"
 #include "solver/ITL_Solver.h"
+#include "solver/BITL_Solver.h"
 #include "solver/ITL_Preconditioner.h"
 #include "solver/UmfPackSolver.h"
 #include "solver/KrylovPreconditioner.h"
@@ -61,6 +62,9 @@ namespace AMDiS {
   {
     LinearSolverCreator *creator;
 
+    // creators for MTL4 krylov solvers
+    // _________________________________________________________________________
+    
     creator = new CGSolver::Creator;
     addCreator("mtl_cg", creator);
 
@@ -115,11 +119,64 @@ namespace AMDiS {
     addCreator("mtl_hypre", creator);
 #endif 
     
+    // creators for block krylov solvers
+    // _________________________________________________________________________
+    
+    creator = new B_CGSolver::Creator;
+    addCreator("bmtl_cg", creator);
+
+    creator = new B_CGSSolver::Creator;
+    addCreator("bmtl_cgs", creator);
+
+//     creator = new B_BiCGSolver::Creator;
+//     addCreator("mtl_bicg", creator);
+
+    creator = new B_BiCGStabSolver::Creator;
+    addCreator("bmtl_bicgstab", creator);
+
+//     creator = new B_BiCGStab2Solver::Creator;
+//     addCreator("bmtl_bicgstab2", creator);      // MTL Error
+
+//     creator = new B_BiCGStabEllSolver::Creator;
+//     addCreator("bmtl_bicgstab_ell", creator);   // MTL Error
+
+//     creator = new QMRSolver::Creator;
+//     addCreator("mtl_qmr", creator);
+
+    creator = new B_TFQMRSolver::Creator;
+    addCreator("bmtl_tfqmr", creator);
+
+    creator = new B_GMResSolver::Creator;
+    addCreator("bmtl_gmres", creator);
+
+    creator = new B_GcrSolver::Creator;
+    addCreator("bmtl_gcr", creator);    
+
+    creator = new B_FGMResSolver::Creator;
+    addCreator("bmtl_fgmres", creator);   
+    
+//     creator = new B_IDRsSolver::Creator;
+//     addCreator("bmtl_idr_s", creator);     // MTL Error
+
+    creator = new B_MinResSolver::Creator;
+    addCreator("bmtl_minres", creator);
+
+    creator = new B_PreOnly::Creator;
+    addCreator("bmtl_preonly", creator);
+    addCreator("bmtl_richardson", creator);
+    
+    // creators for PETSc solvers
+    // _________________________________________________________________________
+    
 #if defined HAVE_SEQ_PETSC
     // sequential PETSc-Solver
-    creator = new PetscSolver<>::Creator;
+    creator = new PetscSolver::Creator;
     addCreator("petsc_petsc", creator); // standard creator for petsc solver
     addCreator("petsc", creator);
+    
+    LinearSolverCreator* creator2 = new PetscSolverNested::Creator;
+    addCreator("bpetsc_petsc", creator2); // standard creator for petsc solver
+    addCreator("bpetsc", creator2);
         
     std::map<std::string,std::string>::iterator it;
     PetscParameters params;
@@ -127,7 +184,8 @@ namespace AMDiS {
 	 it!= params.solverMap.end();
 	 it++) {
       CreatorMap< LinearSolver >::addCreator("petsc_" + it->first, creator);
-    }  
+      CreatorMap< LinearSolver >::addCreator("bpetsc_" + it->first, creator2);
+    }
 #endif
   }
 
@@ -159,6 +217,22 @@ namespace AMDiS {
   }
 
 
+  template<>
+  void CreatorMap<ITL_BasePreconditioner<BlockMTLMatrix, MTLTypes::MTLVector> >::addDefaultCreators()
+  {
+    addCreator("no", new BlockIdentityPreconditioner::Creator);
+    addCreator("diag", new BlockDiagonalPreconditioner::Creator);
+  }
+
+
+#if defined HAVE_SEQ_PETSC
+  template<>
+  void CreatorMap<PetscPreconditioner>::addDefaultCreators() { }
+  
+  template<>
+  void CreatorMap<PetscPreconditionerNested>::addDefaultCreators() { }
+#endif
+
   template<>
   void CreatorMap<Estimator>::addDefaultCreators()
   {

AMDiS/src/CreatorMap.h

@@ -43,13 +43,16 @@ namespace AMDiS {
   template<typename BaseClass>
   class CreatorMap
   {
+  public:
+    typedef std::map< std::string, CreatorInterface<BaseClass>* > CreatorMapType;
+    
   public:
     /// Adds a new creator together with the given key to the map.
     static void addCreator(std::string key, CreatorInterface<BaseClass>* creator) 
     {
       FUNCNAME("CreatorMap::addCreator()");
       init();
-      TEST_EXIT(creatorMap[key] == nullptr)
+      TEST_EXIT(creatorMap[key] == NULL)
 	("there is already a creator for key %s\n",key.c_str());
       creatorMap[key] = creator;
     }
@@ -66,8 +69,10 @@ namespace AMDiS {
       FUNCNAME("CreatorMap::getCreator()");
       init();
       CreatorInterface<BaseClass> *creator = creatorMap[key];
-      TEST_EXIT(creator)("No creator for key \"%s\" defined in init file for parameter \"%s\"\n", 
-			 key.c_str(), initFileStr.c_str());
+      TEST_EXIT(creator)
+	("No creator for key \"%s\" defined in init file for parameter \"%s\"\n", 
+	  key.c_str(), initFileStr.c_str());
+	
       return creator;
     }
 
@@ -89,13 +94,13 @@ namespace AMDiS {
 
   protected:
     /// STL map containing the pairs of keys and creators.
-    static std::map< std::string, CreatorInterface<BaseClass>* > creatorMap;
+    static CreatorMapType creatorMap;
 
     static bool initialized;
   };
 
   template<typename BaseClass>
-  std::map< std::string, CreatorInterface<BaseClass>* > CreatorMap<BaseClass>::creatorMap;
+  typename CreatorMap<BaseClass>::CreatorMapType CreatorMap<BaseClass>::creatorMap;
 
   template<typename BaseClass>
   bool CreatorMap<BaseClass>::initialized = false;

AMDiS/src/DOFMatrix.cc

@@ -67,10 +67,6 @@ namespace AMDiS {
 
     if (!colFeSpace)
       colFeSpace = rowFeSpace;
-#if 0
-    if (rowFeSpace && rowFeSpace->getAdmin())
-      (const_cast<DOFAdmin*>(rowFeSpace->getAdmin()))->addDOFIndexed(this);
-#endif
     boundaryManager = new BoundaryManager(rowFeSpace);
 
     nRow = rowFeSpace->getBasisFcts()->getNumber();
@@ -87,10 +83,6 @@ namespace AMDiS {
     FUNCNAME("DOFMatrix::DOFMatrix()");
 
     *this = rhs;
-#if 0
-    if (rowFeSpace && rowFeSpace->getAdmin())
-      (const_cast<DOFAdmin*>( rowFeSpace->getAdmin()))->addDOFIndexed(this);
-#endif
     TEST_EXIT(rhs.inserter == 0)("Cannot copy during insertion!\n");
     inserter = 0;
   }
@@ -98,10 +90,6 @@ namespace AMDiS {
 
   DOFMatrix::~DOFMatrix()
   {
-#if 0
-    if (rowFeSpace && rowFeSpace->getAdmin())
-      (const_cast<DOFAdmin*>(rowFeSpace->getAdmin()))->removeDOFIndexed(this);
-#endif
     if (boundaryManager) 
       delete boundaryManager;
     if (inserter) 
@@ -116,43 +104,6 @@ namespace AMDiS {
   }
 
 
-  bool DOFMatrix::symmetric()
-  {
-    double tol = 1e-5;
-
-    using mtl::tag::major; using mtl::tag::nz; using mtl::begin; using mtl::end;
-    namespace traits= mtl::traits;
-    typedef base_matrix_type   Matrix;
-
-    traits::row<Matrix>::type                                 row(matrix);
-    traits::col<Matrix>::type                                 col(matrix);
-    traits::const_value<Matrix>::type                         value(matrix);
-
-    typedef traits::range_generator<major, Matrix>::type      cursor_type;
-    typedef traits::range_generator<nz, cursor_type>::type    icursor_type;
-    
-    for (cursor_type cursor = begin<major>(matrix), cend = end<major>(matrix); cursor != cend; ++cursor)
-      for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); icursor != icend; ++icursor)
-	// Compare each non-zero entry with its transposed
-	if (abs(value(*icursor) - matrix[col(*icursor)][row(*icursor)]) > tol)
-	  return false;
-    return true;
-  }
-
-
-  void DOFMatrix::test()
-  {
-    FUNCNAME("DOFMatrix::test()");
-
-    int non_symmetric = !symmetric();
-
-    if (non_symmetric)
-      MSG("Matrix `%s' not symmetric.\n", name.data());
-    else
-      MSG("Matrix `%s' is symmetric.\n", name.data());
-  }
-
-
   DOFMatrix& DOFMatrix::operator=(const DOFMatrix& rhs)
   {
     rowFeSpace = rhs.rowFeSpace;
@@ -290,14 +241,11 @@ namespace AMDiS {
   }
 
 
+#if 0
   double DOFMatrix::logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const
   {
     return matrix[a][b];
   }
-
-#if 0
-  void DOFMatrix::freeDOFContent(int index)
-  {}
 #endif
 
   void DOFMatrix::assemble(double factor, 
@@ -418,7 +366,7 @@ namespace AMDiS {
 
   void DOFMatrix::finishAssembling()
   {
-    // call the operatos cleanup procedures
+    // call the operators cleanup procedures
     for (std::vector<Operator*>::iterator it = operators.begin();
 	 it != operators.end(); ++it)
       (*it)->finishAssembling();
@@ -436,7 +384,7 @@ namespace AMDiS {
     return fillFlag;
   }
 
-
+#if 0
   void DOFMatrix::axpy(double a, const DOFMatrix& x, const DOFMatrix& y)
   {
     matrix += a * x.matrix + y.matrix;
@@ -447,7 +395,7 @@ namespace AMDiS {
   {
     matrix *= b;
   }
-
+#endif
 
   void DOFMatrix::addOperator(Operator *op, double* factor, double* estFactor) 
   { 

AMDiS/src/DOFMatrix.h

@@ -69,6 +69,9 @@ namespace AMDiS {
     /// Symbolic constant for an unused entry which is not followed by any
     /// used entry in this row
     static const int NO_MORE_ENTRIES = -2;
+    
+    /// Minimal slot-size used for the inserter. See \ref startInsertion for details.
+    static const int MIN_NNZ_PER_ROW = 20;
 
   public:    
     DOFMatrix();
@@ -100,56 +103,7 @@ namespace AMDiS {
     {
 	return matrix;
     }
-#if 0
-    // Only to get rid of the abstract functions, I hope they are not used
-    std::vector<bool>::iterator begin() 
-    {
-      ERROR_EXIT("Shouldn't be used, only fake."); std::vector<bool> v; 
-      return v.begin();
-    }
-    
-    std::vector<bool>::iterator end() 
-    {
-      ERROR_EXIT("Shouldn't be used, only fake."); std::vector<bool> v; 
-      return v.end();
-    }
-
-    // Only to get rid of the abstract functions, I hope they are not used
-    std::vector<bool>::const_iterator begin() const
-    {
-      ERROR_EXIT("Shouldn't be used, only fake."); std::vector<bool> v;
-      return v.begin();
-    }
-
-    std::vector<bool>::const_iterator end() const
-    {
-      ERROR_EXIT("Shouldn't be used, only fake."); std::vector<bool> v;
-      return v.end();
-    }
-    
-    bool dummy; // Must be deleted later
-    reference operator[](DegreeOfFreedom i) 
-    {
-      ERROR_EXIT("Shouldn't be used, only fake."); 
-      std::vector<bool> d; d.push_back(dummy);
-      return d[0];
-    }
-
-    const_reference operator[](DegreeOfFreedom i) const 
-    {
-      ERROR_EXIT("Shouldn't be used, only fake."); 
-      std::vector<bool> d; d.push_back(dummy);
-      return d[0];
-    }
- 
-    /// DOFMatrix does not need to be compressed before assembling, when using MTL4.
-    void compressDOFIndexed(int first, int last, std::vector<DegreeOfFreedom> &newDOF) 
-    {}
-    
-    /// Implements DOFIndexedBase::freeDOFContent()
-    virtual void freeDOFContent(int index);
 
-#endif
     /// Returns \ref coupleMatrix.
     inline bool isCoupleMatrix() 
     { 
@@ -162,11 +116,13 @@ namespace AMDiS {
       coupleMatrix = c; 
     }
 
+#if 0
     /// a * x + y
     void axpy(double a, const DOFMatrix& x, const DOFMatrix& y);
 
     /// Multiplication with a scalar.
     void scal(double s);
+#endif
 
     /// Adds an operator to the DOFMatrix. A factor, that is multipled to the 
     /// operator, and a multilier factor for the estimator may be also given.
@@ -267,7 +223,7 @@ namespace AMDiS {
     /// entries per row (per column for CCS matrices).  Choosing this parameter
     /// too small can induce perceivable overhead for compressed matrices. Thus,
     /// it's better to choose a bit too large than too small.
-    void startInsertion(int nnz_per_row = 10);
+    void startInsertion(int nnz_per_row = MIN_NNZ_PER_ROW);
 
     /// Finishes insertion. For compressed matrix types, this is where the
     /// compression happens.
@@ -281,15 +237,6 @@ namespace AMDiS {
       inserter= 0;
     }
 
-#if 0
-    /// Returns whether restriction should be performed after coarsening
-    /// (false by default)
-    virtual bool coarseRestrict() 
-    {
-      return false;
-    }
-#endif
-
     /// Returns const \ref rowFeSpace
     const FiniteElemSpace* getRowFeSpace() const 
     { 
@@ -332,12 +279,6 @@ namespace AMDiS {
       return name; 
     }
 #if 0
-    /// Resizes \ref matrix to n rows
-    inline void resize(int n) 
-    { 
-      TEST_EXIT_DBG(n >= 0)("Can't resize DOFMatrix to negative size\n");
-    }
-#endif
     /// Returns value at logical indices a,b
     double logAcc(DegreeOfFreedom a, DegreeOfFreedom b) const;
 
@@ -349,6 +290,7 @@ namespace AMDiS {
     void addSparseDOFEntry(double sign,
 			   int irow, int jcol, double entry,
 			   bool add = true);
+#endif
 
     void clearDirichletRows();
 
@@ -361,11 +303,6 @@ namespace AMDiS {
 	set_to_zero(matrix);
     }
 
-    /// Test whether \ref matrix is symmetric. Exits if not.
-    void test();
-
-    bool symmetric();
-
     inline std::vector<Operator*>& getOperators() 
     { 
       return operators; 
@@ -398,8 +335,8 @@ namespace AMDiS {
 
       if (num_rows(matrix) != 0)
 	nnzPerRow = int(double(matrix.nnz()) / num_rows(matrix) * 1.2); 
-      if (nnzPerRow < 20) 
-	nnzPerRow = 20;
+      if (nnzPerRow < MIN_NNZ_PER_ROW) 
+	nnzPerRow = MIN_NNZ_PER_ROW;
     }
 
     /// Returns \ref nnzPerRow.

AMDiS/src/DOFVector.cc

@@ -110,7 +110,7 @@ namespace AMDiS {
       delete oldElInfo;
     } else
       inside = mesh->findElInfoAtPoint(p, elInfo, lambda, nullptr, nullptr, nullptr);
-
+    
     if (oldElInfo)
       oldElInfo = elInfo;
 
@@ -122,7 +122,7 @@ namespace AMDiS {
       value = basFcts->evalUh(lambda, uh);
     } else {
 #ifdef HAVE_PARALLEL_DOMAIN_AMDIS
-      return 0.0;
+      value = 0.0;
 #else
       ERROR_EXIT("Can not eval DOFVector at point p, because point is outside geometry.");
 #endif
@@ -132,6 +132,9 @@ namespace AMDiS {
     if (oldElInfo == nullptr)
       delete elInfo;
 
+#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
+    Parallel::mpi::globalAdd(value);
+#endif
     return value;
   }
 

AMDiS/src/GenericOperatorTerm.h

@@ -315,7 +315,7 @@ struct GenericSecondOrderTerm_1 : public SecondOrderTerm
       l1lt(grdLambda, LALt[iq], term(iq));
   }
 
-  /// Implenetation of SecondOrderTerm::eval().
+  /// Implemetation of SecondOrderTerm::eval().
   void eval(int nPoints,
 	    const mtl::dense_vector<double>& uhAtQP,
 	    const mtl::dense_vector<WorldVector<double> >& grdUhAtQP,
@@ -336,7 +336,7 @@ struct GenericSecondOrderTerm_1 : public SecondOrderTerm
     }
   }
 
-  /// Implenetation of SecondOrderTerm::weakEval().
+  /// Implemetation of SecondOrderTerm::weakEval().
   void weakEval(const std::vector<WorldVector<double> > &grdUhAtQP,
 		std::vector<WorldVector<double> > &result) 
   {
@@ -521,25 +521,6 @@ inline void addZOT(Operator& op, double term)
   addZOT(op, constant(term));
 }
 
-
-// .............................................................................
-// Rosenbrock method
-
-// template <typename Term>
-// void addZOT_(Operator* op, const Term& term)
-// {
-//   if (op->getUhOld()) {
-//     addZOT(op, jacobian(term,0) * ...
-//     
-// }
-
-// template <typename Term>
-// void addZOT(Operator& op, const Term& term)
-// {
-//   op.addZeroOrderTerm(new GenericZeroOrderTerm<Term>(term));
-// }
-
-
 // _____________________________________________________________________________
 
 // first order term using FirstOrderTerm::l1
@@ -775,7 +756,7 @@ inline expressions::Function1<expressions::Wrapper<TOut,WorldVector<double> >, e
 eval(AbstractFunction<TOut, WorldVector<double> >* fct) { return function_(wrap(fct), X()); }
 
 /// Integrate an expression over a domain.
-/** IF the expression does not contain any DOFVector, a mesh must be given as second argument */
+/** If the expression does not contain any DOFVector, a mesh must be given as second argument */
 template<typename Term>
 inline typename boost::enable_if<typename traits::is_expr<Term>::type, typename Term::value_type>::type
 integrate(Term term, Mesh* mesh_ = NULL);

AMDiS/src/Global.cc

@@ -188,7 +188,11 @@ namespace AMDiS {
     PRINT_LINE((*error), buff);
     va_end(arg);
 #if defined HAVE_PARALLEL_DOMAIN_AMDIS && !defined HAVE_PARALLEL_MTL4 && (DEBUG == 0 || defined NDEBUG)
+  #if (PETSC_VERSION_MINOR >= 5)
+    PetscError(MPI_COMM_WORLD, __LINE__, "Msg::print_error_exit", "Global.cc", 1, PETSC_ERROR_INITIAL, buff);
+  #else
     PetscError(MPI_COMM_WORLD, __LINE__, "Msg::print_error_exit", "Global.cc", "src/", 1, PETSC_ERROR_INITIAL, buff);
+  #endif
 #else
     throw std::runtime_error(buff);
 #endif

AMDiS/src/LinearAlgebra.h

+/******************************************************************************
+ *
+ * AMDiS - Adaptive multidimensional simulations
+ *
+ * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
+ * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
+ *
+ * Authors: 
+ * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
+ *
+ * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
+ * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
+ *
+ *
+ * This file is part of AMDiS
+ *
+ * See also license.opensource.txt in the distribution.
+ * 
+ ******************************************************************************/
+
+ /// This file is work in progress and may replace the assembling now part of
+ /// DOFMatrix directly. The idea is to assemble matrices to a data-structure
+ /// provided by a linear algebra backend. You have to specify how an element-
+ /// matrix can be added to the global matrix, how a Dirichlet-row is inserted
+ /// and how the matrix is initialized/finalized.
+
+
+/** \file LinearAlgebra.h */
+
+#ifndef AMDIS_LINEAR_ALGEBRA_H
+#define AMDIS_LINEAR_ALGEBRA_H
+
+#include "utility/calculate_nnz.hpp"
+
+namespace AMDiS {
+    
+  namespace tag 
+  {
+    struct distributed {};
+    struct non_distributed {};
+    
+    // backends
+    struct mtl {};
+    struct pmtl {};
+    struct petsc {};
+    struct p_petsc {};
+  }
+  
+  namespace traits
+  {
+    template< class T >
+    struct distributed_tag : boost::mpl::if_< 
+      mtl::traits::is_distributed< T >,
+      tag::distributed,
+      tag::non_distributed
+    > {};
+  }
+  
+  namespace detail
+  {
+    template< class Dummy = void >
+    struct DefaultParameters {};
+  }
+  
+  /** fill for each backend a linear algebra class that
+   *  contains the following typedefs and functions
+   **/
+  template< class Backend = tag::mtl, class Value = double, class Parameters = detail::DefaultParameters<> >
+  struct LinearAlgebra 
+  {
+  public: // typedefs
+    
+    typedef Value         ValueType;
+    typedef unsigned int  SizeType;
+    
+    typedef mtl::matrix::parameters<mtl::row_major, mtl::index::c_index, mtl::non_fixed::dimensions, false, SizeType> MatrixPara;
+    typedef mtl::matrix::compressed2D<ValueType, MatrixPara>   MatrixBase;
+    typedef mtl::vector::dense_vector<ValueType>               VectorBase;
+    
+    typedef mtl::matrix::inserter<MatrixBase, mtl::operations::update_plus<ValueType> >  Inserter;
+    typedef MPI::Intracomm   Communicator;
+    
+    static const int minSlotSize = 20;
+    
+  public: // methods
+    
+    /// constructor
+    LinearAlgebra(Communicator& comm)
+      : comm(comm), ins(NULL)
+    { }
+