Merge pull request #1 from andlaus/fix_fem_solver_backend

Fixes for the dune-fem linear solver backend

Author: dr-robertk

ewoms/common/basicproperties.hh

@@ -78,10 +78,6 @@ 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

@@ -50,7 +50,6 @@
 #include <ewoms/parallel/gridcommhandles.hh>
 #include <ewoms/parallel/threadmanager.hh>
 #include <ewoms/linear/nullborderlistmanager.hh>
-#include <ewoms/linear/istlsparsematrixadapter.hh>
 #include <ewoms/common/simulator.hh>
 #include <ewoms/common/alignedallocator.hh>
 #include <ewoms/common/timer.hh>
@@ -76,13 +75,7 @@
 #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>
@@ -136,78 +129,6 @@ 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
-#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:
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
-    typedef Ewoms::MatrixBlock<Scalar, numEq, numEq> Block;
-
-public:
-    typedef typename Ewoms::Linear::IstlSparseMatrixAdapter<Block> type;
-};
-#endif
-
 //! The maximum allowed number of timestep divisions for the
 //! Newton solver
 SET_INT_PROP(FvBaseDiscretization, MaxTimeStepDivisions, 10);
@@ -408,15 +329,13 @@ 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 DiscreteFunctionSpace& space)
-            : blockVector_(space.size())
+        BlockVectorWrapper(const std::string& name OPM_UNUSED, const size_t size)
+            : blockVector_(size)
         {}
 
         SolutionVector& blockVector()
@@ -426,12 +345,14 @@ class FvBaseDiscretization
     };
 
 #if HAVE_DUNE_FEM
+    typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunctionSpace)    DiscreteFunctionSpace;
+
     // discrete function storing solution data
-    typedef typename GET_PROP_TYPE(TypeTag, DiscreteFunction)         DiscreteFunction;
+    typedef Dune::Fem::ISTLBlockVectorDiscreteFunction<DiscreteFunctionSpace, PrimaryVariables> 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;
@@ -440,6 +361,7 @@ 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
@@ -459,14 +381,14 @@ public:
         , elementMapper_(gridView_)
         , vertexMapper_(gridView_)
 #endif
+        , newtonMethod_(simulator)
+        , localLinearizer_(ThreadManager::maxThreads())
+        , linearizer_(new Linearizer())
 #if HAVE_DUNE_FEM
-        , discreteFunctionSpace_( simulator.vanguard().gridPart() )
+        , space_( simulator.vanguard().gridPart() )
 #else
-        , discreteFunctionSpace_( asImp_().numGridDof() )
+        , space_( 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))
@@ -491,7 +413,7 @@ public:
 
         size_t numDof = asImp_().numGridDof();
         for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
-            solution_[timeIdx].reset(new DiscreteFunction("solution", discreteFunctionSpace_));
+            solution_[timeIdx].reset(new DiscreteFunction("solution", space_));
 
             if (storeIntensiveQuantities()) {
                 intensiveQuantityCache_[timeIdx].resize(numDof);
@@ -1156,16 +1078,6 @@ 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
@@ -1583,7 +1495,7 @@ public:
     void resetLinearizer ()
     {
         delete linearizer_;
-        linearizer_ = new Linearizer();
+        linearizer_ = new Linearizer;
         linearizer_->init(simulator_);
     }
 
@@ -1817,11 +1729,6 @@ public:
             solution(timeIdx).resize(numDof);
 
         auxMod->applyInitial();
-
-#if DUNE_VERSION_NEWER( DUNE_FEM, 2, 7 )
-        discreteFunctionSpace_.extendSize( asImp_().numAuxiliaryDof() );
-#endif
-
     }
 
     /*!
@@ -1888,11 +1795,6 @@ public:
     const Ewoms::Timer& updateTimer() const
     { return updateTimer_; }
 
-#if HAVE_DUNE_FEM
-    const DiscreteFunctionSpace& space() const { return discreteFunctionSpace_; }
-#endif
-
-
 protected:
     void resizeAndResetIntensiveQuantitiesCache_()
     {
@@ -1963,18 +1865,9 @@ 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_;
@@ -1993,6 +1886,15 @@ 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,7 +54,6 @@ namespace Ewoms {
 template<class TypeTag>
 class EcfvDiscretization;
 
-
 /*!
  * \ingroup FiniteVolumeDiscretizations
  *

ewoms/disc/common/fvbaseproperties.hh

@@ -37,9 +37,6 @@
 #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
 
@@ -57,19 +54,10 @@ 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);
@@ -92,6 +80,9 @@ 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,24 +91,6 @@ 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

ewoms/disc/vcfv/vcfvdiscretization.hh

@@ -101,24 +101,6 @@ public:
     // Lagrange discrete function space with unknowns at the cell vertices
     typedef Dune::Fem::LagrangeDiscreteFunctionSpace< FunctionSpace, GridPart, 1 > type;
 };
-#else
-SET_PROP(VcfvDiscretization, 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 vertices