[feature][solver] This adds a FemSolverBackend and the

AMGXSolverBackend and various smaller fixes to accommodate this.

Author: Robert Kloefkorn

ewoms-prereqs.cmake

@@ -13,8 +13,12 @@ set (ewoms_CONFIG_VAR
   HAVE_DUNE_ISTL
   HAVE_DUNE_ALUGRID
   HAVE_DUNE_FEM
+  HAVE_PETSC
+  HAVE_AMGXSOLVER
   HAVE_ECL_INPUT
   HAVE_ECL_OUTPUT
+  HAVE_VIENNACL
+  HAVE_OPENCL
   DUNE_AVOID_CAPABILITIES_IS_PARALLEL_DEPRECATION_WARNING
   )
 
@@ -36,6 +40,12 @@ set (ewoms_DEPS
   "dune-alugrid"
   "dune-fem"
   "opm-grid"
+  # PETSc numerical backend
+  "PETSc"
+  # AMGX wrapper using PETSc
+  "AmgXSolver"
+  # ViennaCL
+  "ViennaCL"
   # valgrind client requests
   "Valgrind"
   # quadruple precision floating point calculations

ewoms/common/basicproperties.hh

@@ -78,6 +78,10 @@ NEW_PROP_TAG(GridView);
 
 #if HAVE_DUNE_FEM
 NEW_PROP_TAG(GridPart);
+//! The class describing the discrete function space when dune-fem is used, otherwise it points to the stencil class
+NEW_PROP_TAG(DiscreteFunctionSpace);
+NEW_PROP_TAG(DiscreteFunction);
+NEW_PROP_TAG(LinearOperator);
 #endif
 
 //! Property which tells the Vanguard how often the grid should be refined

ewoms/disc/common/fvbasediscretization.hh

@@ -76,7 +76,13 @@
 #include <dune/fem/space/common/restrictprolongtuple.hh>
 #include <dune/fem/function/blockvectorfunction.hh>
 #include <dune/fem/misc/capabilities.hh>
+
+#if HAVE_PETSC
+#include <dune/fem/operator/linear/petscoperator.hh>
 #endif
+#include <dune/fem/operator/linear/istloperator.hh>
+#include <dune/fem/operator/linear/spoperator.hh>
+#endif // endif HAVE_DUNE_FEM
 
 #include <limits>
 #include <list>
@@ -130,7 +136,66 @@ SET_TYPE_PROP(FvBaseDiscretization, DiscExtensiveQuantities, Ewoms::FvBaseExtens
 //! Calculates the gradient of any quantity given the index of a flux approximation point
 SET_TYPE_PROP(FvBaseDiscretization, GradientCalculator, Ewoms::FvBaseGradientCalculator<TypeTag>);
 
-//! Set the type of a global jacobian matrix from the solution types
+//! Set the type of a global Jacobian matrix from the solution types
+#if HAVE_DUNE_FEM
+SET_PROP(FvBaseDiscretization, DiscreteFunction)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunctionSpace) DiscreteFunctionSpace;
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables)      PrimaryVariables;
+public:
+    // discrete function storing solution data
+    typedef Dune::Fem::ISTLBlockVectorDiscreteFunction<DiscreteFunctionSpace, PrimaryVariables> type;
+};
+#endif
+
+#if USE_DUNE_FEM_SOLVERS
+SET_PROP(FvBaseDiscretization, SparseMatrixAdapter)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunctionSpace) DiscreteFunctionSpace;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar)                Scalar;
+    // discrete function storing solution data
+    typedef Dune::Fem::ISTLBlockVectorDiscreteFunction<DiscreteFunctionSpace> DiscreteFunction;
+
+#if USE_DUNE_FEM_PETSC_SOLVERS
+#warning "Using Dune-Fem PETSc solvers"
+    typedef Dune::Fem::PetscLinearOperator< DiscreteFunction, DiscreteFunction > LinearOperator;
+#elif USE_DUNE_FEM_VIENNACL_SOLVERS
+#warning "Using Dune-Fem ViennaCL solvers"
+    typedef Dune::Fem::SparseRowLinearOperator < DiscreteFunction, DiscreteFunction > LinearOperator;
+#else
+#warning "Using Dune-Fem ISTL solvers"
+    typedef Dune::Fem::ISTLLinearOperator < DiscreteFunction, DiscreteFunction > LinearOperator;
+#endif
+
+    struct FemMatrixBackend : public LinearOperator
+    {
+        typedef LinearOperator  ParentType;
+        typedef typename LinearOperator :: MatrixType    Matrix;
+        typedef typename ParentType :: MatrixBlockType   MatrixBlock;
+        template <class Simulator>
+        FemMatrixBackend( const Simulator& simulator )
+            : LinearOperator("eWoms::Jacobian", simulator.model().space(), simulator.model().space() )
+        {}
+
+        void commit()
+        {
+          this->flushAssembly();
+        }
+
+        template< class LocalBlock >
+        void addToBlock ( const size_t row, const size_t col, const LocalBlock& block )
+        {
+          this->addBlock( row, col, block );
+        }
+
+        void clearRow( const size_t row, const Scalar diag = 1.0 ) { this->unitRow( row ); }
+    };
+public:
+    typedef FemMatrixBackend type;
+};
+#else
 SET_PROP(FvBaseDiscretization, SparseMatrixAdapter)
 {
 private:
@@ -141,6 +206,7 @@ private:
 public:
     typedef typename Ewoms::Linear::IstlSparseMatrixAdapter<Block> type;
 };
+#endif
 
 //! The maximum allowed number of timestep divisions for the
 //! Newton solver
@@ -342,13 +408,15 @@ class FvBaseDiscretization
 
     typedef typename LocalResidual::LocalEvalBlockVector LocalEvalBlockVector;
 
+    typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunctionSpace)    DiscreteFunctionSpace;
+
     class BlockVectorWrapper
     {
     protected:
         SolutionVector blockVector_;
     public:
-        BlockVectorWrapper(const std::string& name OPM_UNUSED, const size_t size)
-            : blockVector_(size)
+        BlockVectorWrapper(const std::string& name OPM_UNUSED, const DiscreteFunctionSpace& space)
+            : blockVector_(space.size())
         {}
 
         SolutionVector& blockVector()
@@ -358,14 +426,12 @@ class FvBaseDiscretization
     };
 
 #if HAVE_DUNE_FEM
-    typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunctionSpace)    DiscreteFunctionSpace;
-
     // discrete function storing solution data
-    typedef Dune::Fem::ISTLBlockVectorDiscreteFunction<DiscreteFunctionSpace, PrimaryVariables> DiscreteFunction;
+    typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunction)         DiscreteFunction;
 
     // problem restriction and prolongation operator for adaptation
-    typedef typename GET_PROP_TYPE(TypeTag, Problem)   Problem;
-    typedef typename Problem :: RestrictProlongOperator  ProblemRestrictProlongOperator;
+    typedef typename GET_PROP_TYPE(TypeTag, Problem)                  Problem;
+    typedef typename Problem :: RestrictProlongOperator               ProblemRestrictProlongOperator;
 
     // discrete function restriction and prolongation operator for adaptation
     typedef Dune::Fem::RestrictProlongDefault< DiscreteFunction > DiscreteFunctionRestrictProlong;
@@ -374,7 +440,6 @@ class FvBaseDiscretization
     typedef Dune::Fem::AdaptationManager<Grid, RestrictProlong  > AdaptationManager;
 #else
     typedef BlockVectorWrapper  DiscreteFunction;
-    typedef size_t              DiscreteFunctionSpace;
 #endif
 
     // copying a discretization object is not a good idea
@@ -394,14 +459,14 @@ public:
         , elementMapper_(gridView_)
         , vertexMapper_(gridView_)
 #endif
-        , newtonMethod_(simulator)
-        , localLinearizer_(ThreadManager::maxThreads())
-        , linearizer_(new Linearizer())
 #if HAVE_DUNE_FEM
-        , space_( simulator.vanguard().gridPart() )
+        , discreteFunctionSpace_( simulator.vanguard().gridPart() )
 #else
-        , space_( asImp_().numGridDof() )
+        , discreteFunctionSpace_( asImp_().numGridDof() )
 #endif
+        , newtonMethod_(simulator)
+        , localLinearizer_(ThreadManager::maxThreads())
+        , linearizer_(new Linearizer( ))
         , enableGridAdaptation_( EWOMS_GET_PARAM(TypeTag, bool, EnableGridAdaptation) )
         , enableIntensiveQuantityCache_(EWOMS_GET_PARAM(TypeTag, bool, EnableIntensiveQuantityCache))
         , enableStorageCache_(EWOMS_GET_PARAM(TypeTag, bool, EnableStorageCache))
@@ -426,7 +491,7 @@ public:
 
         size_t numDof = asImp_().numGridDof();
         for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
-            solution_[timeIdx].reset(new DiscreteFunction("solution", space_));
+            solution_[timeIdx].reset(new DiscreteFunction("solution", discreteFunctionSpace_));
 
             if (storeIntensiveQuantities()) {
                 intensiveQuantityCache_[timeIdx].resize(numDof);
@@ -1091,6 +1156,16 @@ public:
     SolutionVector& solution(unsigned timeIdx)
     { return solution_[timeIdx]->blockVector(); }
 
+    template <class BVector>
+    void communicate( BVector& x ) const
+    {
+#if HAVE_DUNE_FEM
+        typedef Dune::Fem::ISTLBlockVectorDiscreteFunction<DiscreteFunctionSpace, typename BVector::block_type> DF;
+        DF tmpX("temp-x", discreteFunctionSpace_, x);
+        tmpX.communicate();
+#endif
+    }
+
   protected:
     /*!
      * \copydoc solution(int) const
@@ -1508,7 +1583,7 @@ public:
     void resetLinearizer ()
     {
         delete linearizer_;
-        linearizer_ = new Linearizer;
+        linearizer_ = new Linearizer();
         linearizer_->init(simulator_);
     }
 
@@ -1742,6 +1817,11 @@ public:
             solution(timeIdx).resize(numDof);
 
         auxMod->applyInitial();
+
+#if DUNE_VERSION_NEWER( DUNE_FEM, 2, 7 )
+        discreteFunctionSpace_.extendSize( asImp_().numAuxiliaryDof() );
+#endif
+
     }
 
     /*!
@@ -1808,6 +1888,11 @@ public:
     const Ewoms::Timer& updateTimer() const
     { return updateTimer_; }
 
+#if HAVE_DUNE_FEM
+    const DiscreteFunctionSpace& space() const { return discreteFunctionSpace_; }
+#endif
+
+
 protected:
     void resizeAndResetIntensiveQuantitiesCache_()
     {
@@ -1878,9 +1963,18 @@ protected:
     ElementMapper elementMapper_;
     VertexMapper vertexMapper_;
 
+    DiscreteFunctionSpace discreteFunctionSpace_;
+
     // a vector with all auxiliary equations to be considered
     std::vector<BaseAuxiliaryModule<TypeTag>*> auxEqModules_;
 
+    mutable std::array< std::unique_ptr< DiscreteFunction >, historySize > solution_;
+
+#if HAVE_DUNE_FEM
+    std::unique_ptr< RestrictProlong  > restrictProlong_;
+    std::unique_ptr< AdaptationManager> adaptationManager_;
+#endif
+
     NewtonMethod newtonMethod_;
 
     Ewoms::Timer prePostProcessTimer_;
@@ -1899,15 +1993,6 @@ protected:
     mutable IntensiveQuantitiesVector intensiveQuantityCache_[historySize];
     mutable std::vector<bool> intensiveQuantityCacheUpToDate_[historySize];
 
-    DiscreteFunctionSpace space_;
-    mutable std::array< std::unique_ptr< DiscreteFunction >, historySize > solution_;
-
-#if HAVE_DUNE_FEM
-    std::unique_ptr<RestrictProlong> restrictProlong_;
-    std::unique_ptr<AdaptationManager> adaptationManager_;
-#endif
-
-
     std::list<BaseOutputModule<TypeTag>*> outputModules_;
 
     Scalar gridTotalVolume_;

ewoms/disc/common/fvbaselinearizer.hh

@@ -54,6 +54,7 @@ namespace Ewoms {
 template<class TypeTag>
 class EcfvDiscretization;
 
+
 /*!
  * \ingroup FiniteVolumeDiscretizations
  *
@@ -95,8 +96,6 @@ class FvBaseLinearizer
 
     typedef GlobalEqVector Vector;
 
-    typedef typename SparseMatrixAdapter::IstlMatrix IstlMatrix;
-
     enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
     enum { historySize = GET_PROP_VALUE(TypeTag, TimeDiscHistorySize) };
 
@@ -229,9 +228,6 @@ public:
      */
     void linearizeAuxiliaryEquations()
     {
-        // flush possible local caches into matrix structure
-        jacobian_->commit();
-
         auto& model = model_();
         const auto& comm = simulator_().gridView().comm();
         for (unsigned auxModIdx = 0; auxModIdx < model.numAuxiliaryModules(); ++auxModIdx) {
@@ -262,6 +258,9 @@ public:
             if (!succeeded)
                 throw Opm::NumericalIssue("linearization of an auxilary equation failed");
         }
+
+        // flush possible local caches into matrix structure
+        jacobian_->commit();
     }
 
     /*!
@@ -499,6 +498,9 @@ private:
             std::rethrow_exception(exceptionPtr);
         }
 
+        // flush possible local caches into matrix structure
+        jacobian_->commit();
+
         applyConstraintsToLinearization_();
     }
 

ewoms/disc/common/fvbaseproperties.hh

@@ -37,6 +37,9 @@
 #include <ewoms/common/basicproperties.hh>
 #include <ewoms/io/vtkprimaryvarsmodule.hh>
 #include <ewoms/linear/parallelbicgstabbackend.hh>
+#include <ewoms/linear/parallelistlbackend.hh>
+#include <ewoms/linear/femsolverbackend.hh>
+#include <ewoms/linear/amgxsolverbackend.hh>
 
 BEGIN_PROPERTIES
 
@@ -54,10 +57,19 @@ NEW_PROP_TAG(ParallelBiCGStabLinearSolver);
 NEW_PROP_TAG(LocalLinearizerSplice);
 NEW_PROP_TAG(FiniteDifferenceLocalLinearizer);
 
+NEW_PROP_TAG(DiscreteFunctionSpace);
+
 SET_SPLICES(FvBaseDiscretization, LinearSolverSplice, LocalLinearizerSplice);
 
 //! use a parallel BiCGStab linear solver by default
+#if USE_AMGX_SOLVERS
+SET_TAG_PROP(FvBaseDiscretization, LinearSolverSplice, AmgXSolverBackend);
+#elif USE_DUNE_FEM_SOLVERS
+SET_TAG_PROP(FvBaseDiscretization, LinearSolverSplice, FemSolverBackend);
+#else
 SET_TAG_PROP(FvBaseDiscretization, LinearSolverSplice, ParallelBiCGStabLinearSolver);
+//SET_TAG_PROP(FvBaseDiscretization, LinearSolverSplice, ParallelIstlLinearSolver);
+#endif
 
 //! by default, use finite differences to linearize the system of PDEs
 SET_TAG_PROP(FvBaseDiscretization, LocalLinearizerSplice, FiniteDifferenceLocalLinearizer);
@@ -80,9 +92,6 @@ NEW_PROP_TAG(GridView);
 //! The class describing the stencil of the spatial discretization
 NEW_PROP_TAG(Stencil);
 
-//! The class describing the discrete function space when dune-fem is used, otherwise it points to the stencil class
-NEW_PROP_TAG(DiscreteFunctionSpace);
-
 //! The type of the problem
 NEW_PROP_TAG(Problem);
 //! The type of the base class for all problems which use this model

ewoms/disc/ecfv/ecfvdiscretization.hh

@@ -91,6 +91,24 @@ private:
 public:
     typedef Dune::Fem::FiniteVolumeSpace< FunctionSpace, GridPart, 0 > type;
 };
+#else
+SET_PROP(EcfvDiscretization, DiscreteFunctionSpace)
+{
+private:
+    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+    struct DiscreteFunctionSpace
+    {
+        static const int dimRange = numEq ;
+        size_t numGridDofs_;
+        size_t extension_;
+        DiscreteFunctionSpace( const size_t numGridDofs )
+            : numGridDofs_( numGridDofs ), extension_(0) {}
+        void extendSize( const size_t extension ) { extension_ = extension; }
+        size_t size() const { return dimRange * (numGridDofs_ + extension_); }
+    };
+public:
+    typedef DiscreteFunctionSpace type;
+};
 #endif
 
 //! Set the border list creator for to the one of an element based